/[ascend]/trunk/models/johnpye/fprops/helmholtz.c
ViewVC logotype

Contents of /trunk/models/johnpye/fprops/helmholtz.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1909 - (show annotations) (download) (as text)
Sun Sep 28 07:41:24 2008 UTC (13 years, 9 months ago) by jpye
File MIME type: text/x-csrc
File size: 15728 byte(s)
The helmholtz_p function now provides its own derivatives to ASCEND.
Restored full test suite for Hydrogen.
1 /* ASCEND modelling environment
2 Copyright (C) 2008 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 *//** @file
19 Implementation of the reduced molar Helmholtz free energy equation of state.
20
21 For nomenclature see Tillner-Roth, Harms-Watzenberg and Baehr, Eine neue
22 Fundamentalgleichung f端r Ammoniak.
23
24 John Pye, 29 Jul 2008.
25 */
26
27 #include <math.h>
28
29 #include "helmholtz.h"
30 #include "ideal_impl.h"
31
32 #ifdef TEST
33 #include <assert.h>
34 #include <stdlib.h>
35 #include <stdio.h>
36 #endif
37
38 #define SQ(X) ((X)*(X))
39
40 /* forward decls */
41
42 static double helm_resid(double tau, double delta, const HelmholtzData *data);
43 static double helm_resid_del(double tau, double delta, const HelmholtzData *data);
44 static double helm_resid_tau(double tau, double delta, const HelmholtzData *data);
45 static double helm_resid_deltau(double tau, double delta, const HelmholtzData *data);
46 static double helm_resid_deldel(double tau, double delta, const HelmholtzData *data);
47
48 /**
49 Function to calculate pressure from Helmholtz free energy EOS, given temperature
50 and mass density.
51
52 @param T temperature in K
53 @param rho mass density in kg/m続
54 @return pressure in Pa???
55 */
56 double helmholtz_p(double T, double rho, const HelmholtzData *data){
57
58 double tau = data->T_star / T;
59 double delta = rho / data->rho_star;
60
61 #ifdef TEST
62 assert(data->rho_star!=0);
63 assert(T!=0);
64 assert(!isnan(tau));
65 assert(!isnan(delta));
66 assert(!isnan(data->R));
67
68 //fprintf(stderr,"p calc: T = %f\n",T);
69 //fprintf(stderr,"p calc: tau = %f\n",tau);
70 //fprintf(stderr,"p calc: rho = %f\n",rho);
71 //fprintf(stderr,"p calc: delta = %f\n",delta);
72 //fprintf(stderr,"p calc: R*T*rho = %f\n",data->R * T * rho);
73
74 //fprintf(stderr,"T = %f\n", T);
75 //fprintf(stderr,"rhob = %f, rhob* = %f, delta = %f\n", rho/data->M, data->rho_star/data->M, delta);
76 #endif
77
78 return data->R * T * rho * (1 + delta * helm_resid_del(tau,delta,data));
79 }
80
81 /**
82 Function to calculate internal energy from Helmholtz free energy EOS, given
83 temperature and mass density.
84
85 @param T temperature in K
86 @param rho mass density in kg/m続
87 @return internal energy in ???
88 */
89 double helmholtz_u(double T, double rho, const HelmholtzData *data){
90
91 double tau = data->T_star / T;
92 double delta = rho / data->rho_star;
93
94 #ifdef TEST
95 assert(data->rho_star!=0);
96 assert(T!=0);
97 assert(!isnan(tau));
98 assert(!isnan(delta));
99 assert(!isnan(data->R));
100 #endif
101
102 #if 0
103 fprintf(stderr,"ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
104 fprintf(stderr,"resid_tau = %f\n",helm_resid_tau(tau,delta,data));
105 fprintf(stderr,"R T = %f\n",data->R * data->T_star);
106 #endif
107
108 return data->R * data->T_star * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data));
109 }
110
111 /**
112 Function to calculate enthalpy from Helmholtz free energy EOS, given
113 temperature and mass density.
114
115 @param T temperature in K
116 @param rho mass density in kg/m続
117 @return enthalpy in J/kg
118 */
119 double helmholtz_h(double T, double rho, const HelmholtzData *data){
120
121 double tau = data->T_star / T;
122 double delta = rho / data->rho_star;
123
124 #ifdef TEST
125 assert(data->rho_star!=0);
126 assert(T!=0);
127 assert(!isnan(tau));
128 assert(!isnan(delta));
129 assert(!isnan(data->R));
130 #endif
131
132 return data->R * T * (1 + tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data)) + delta*helm_resid_del(tau,delta,data));
133 }
134
135 /**
136 Function to calculate entropy from Helmholtz free energy EOS, given
137 temperature and mass density.
138
139 @param T temperature in K
140 @param rho mass density in kg/m続
141 @return entropy in J/kgK
142 */
143 double helmholtz_s(double T, double rho, const HelmholtzData *data){
144
145 double tau = data->T_star / T;
146 double delta = rho / data->rho_star;
147
148 #ifdef ENTROPY_DEBUG
149 assert(data->rho_star!=0);
150 assert(T!=0);
151 assert(!isnan(tau));
152 assert(!isnan(delta));
153 assert(!isnan(data->R));
154
155 fprintf(stderr,"helm_ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
156 fprintf(stderr,"helm_resid_tau = %f\n",helm_resid_tau(tau,delta,data));
157 fprintf(stderr,"helm_ideal = %f\n",helm_ideal(tau,delta,data->ideal));
158 fprintf(stderr,"helm_resid = %f\n",helm_resid(tau,delta,data));
159 #endif
160 return data->R * (
161 tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data))
162 - (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data))
163 );
164 }
165
166 /**
167 Function to calculate Helmholtz energy from the Helmholtz free energy EOS,
168 given temperature and mass density.
169
170 @param T temperature in K
171 @param rho mass density in kg/m続
172 @return Helmholtz energy 'a', in J/kg
173 */
174 double helmholtz_a(double T, double rho, const HelmholtzData *data){
175
176 double tau = data->T_star / T;
177 double delta = rho / data->rho_star;
178
179 #ifdef TEST
180 assert(data->rho_star!=0);
181 assert(T!=0);
182 assert(!isnan(tau));
183 assert(!isnan(delta));
184 assert(!isnan(data->R));
185 #endif
186
187 #ifdef HELMHOLTZ_DEBUG
188 fprintf(stderr,"helmholtz_a: T = %f, rho = %f\n",T,rho);
189 fprintf(stderr,"multiplying by RT = %f\n",data->R*T);
190 #endif
191
192 return data->R * T * (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data));
193 }
194
195
196 /**
197 Calculation zero-pressure specific heat capacity
198 */
199 double helmholtz_cp0(double T, const HelmholtzData *data){
200 double val = helm_cp0(T,data->ideal);
201 #if 0
202 fprintf(stderr,"val = %f\n",val);
203 #endif
204 return val;
205 }
206
207 /**
208 Calculate partial derivative of p with respect to T, with rho constant
209 */
210 double helmholtz_dpdT_rho(double T, double rho, const HelmholtzData *data){
211 double tau = data->T_star / T;
212 double delta = rho / data->rho_star;
213
214 double phir_del = helm_resid_del(tau,delta,data);
215 double phir_deltau = helm_resid_deltau(tau,delta,data);
216 #ifdef TEST
217 assert(!isinf(phir_del));
218 assert(!isinf(phir_deltau));
219 assert(!isnan(phir_del));
220 assert(!isnan(phir_deltau));
221 assert(!isnan(data->R));
222 assert(!isnan(rho));
223 assert(!isnan(tau));
224 #endif
225
226 double res = data->R * rho * (1 + delta*phir_del - delta*tau*phir_deltau);
227
228 #ifdef TEST
229 assert(!isnan(res));
230 assert(!isinf(res));
231 #endif
232 return res;
233 }
234
235 /**
236 Calculate partial derivative of p with respect to rho, with T constant
237 */
238 double helmholtz_dpdrho_T(double T, double rho, const HelmholtzData *data){
239 double tau = data->T_star / T;
240 double delta = rho / data->rho_star;
241
242 double phir_del = helm_resid_del(tau,delta,data);
243 double phir_deldel = helm_resid_deldel(tau,delta,data);
244 #ifdef TEST
245 assert(!isinf(phir_del));
246 assert(!isinf(phir_deldel));
247 #endif
248 return data->R * T * (1 + 2*delta*phir_del + delta*delta* phir_deldel);
249 }
250
251 /*---------------------------------------------
252 UTILITY FUNCTION(S)
253 */
254
255 /* ipow: public domain by Mark Stephen with suggestions by Keiichi Nakasato */
256 static double ipow(double x, int n){
257 double t = 1.0;
258
259 if(!n)return 1.0; /* At the top. x^0 = 1 */
260
261 if (n < 0){
262 n = -n;
263 x = 1.0/x; /* error if x == 0. Good */
264 } /* ZTC/SC returns inf, which is even better */
265
266 if (x == 0.0)return 0.0;
267
268 do{
269 if(n & 1)t *= x;
270 n /= 2; /* KN prefers if (n/=2) x*=x; This avoids an */
271 x *= x; /* unnecessary but benign multiplication on */
272 }while(n); /* the last pass, but the comparison is always
273 true _except_ on the last pass. */
274
275 return t;
276 }
277
278 //#define RESID_DEBUG
279
280 /**
281 Residual part of helmholtz function.
282 */
283 double helm_resid(double tau, double delta, const HelmholtzData *data){
284 double dell,ldell, term, sum, res = 0;
285 unsigned n, i;
286 const HelmholtzPowTerm *pt;
287 const HelmholtzGausTerm *gt;
288
289 n = data->np;
290 pt = &(data->pt[0]);
291
292 #ifdef RESID_DEBUG
293 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
294 #endif
295
296 /* power terms */
297 sum = 0;
298 dell = ipow(delta,pt->l);
299 ldell = pt->l * dell;
300 unsigned oldl;
301 for(i=0; i<n; ++i){
302 term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
303 sum += term;
304 #ifdef RESID_DEBUG
305 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
306 if(pt->l==0){
307 fprintf(stderr,",row=%e\n",term);
308 }else{
309 fprintf(stderr,",row=%e\n,",term*exp(-dell));
310 }
311 #endif
312 oldl = pt->l;
313 ++pt;
314 if(i+1==n || oldl != pt->l){
315 if(oldl == 0){
316 #ifdef RESID_DEBUG
317 fprintf(stderr,"linear ");
318 #endif
319 res += sum;
320 }else{
321 #ifdef RESID_DEBUG
322 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
323 #endif
324 res += sum * exp(-dell);
325 }
326 #ifdef RESID_DEBUG
327 fprintf(stderr,"i = %d, res = %f\n",i,res);
328 #endif
329 sum = 0;
330 dell = ipow(delta,pt->l);
331 ldell = pt->l*dell;
332 }
333 }
334
335 /* gaussian terms */
336 n = data->ng;
337 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
338 gt = &(data->gt[0]);
339 for(i=0; i<n; ++i){
340 #ifdef RESID_DEBUG
341 fprintf(stderr,"i = %d, GAUSSIAN, n = %e, t = %f, d = %f, alpha = %f, beta = %f, gamma = %f, epsilon = %f\n",i+1, gt->n, gt->t, gt->d, gt->alpha, gt->beta, gt->gamma, gt->epsilon);
342 #endif
343 double d1 = delta - gt->epsilon;
344 double t1 = tau - gt->gamma;
345 double e1 = -gt->alpha*d1*d1 - gt->beta*t1*t1;
346 sum = gt->n * pow(tau,gt->t) * pow(delta,gt->d) * exp(e1);
347 //fprintf(stderr,"sum = %f\n",sum);
348 res += sum;
349 ++gt;
350 }
351
352 #ifdef RESID_DEBUG
353 fprintf(stderr,"phir = %f\n",res);
354 #endif
355 return res;
356 }
357
358 /**
359 Derivative of the helmholtz residual function with respect to
360 delta.
361 */
362 double helm_resid_del(double tau,double delta, const HelmholtzData *data){
363 double sum = 0, res = 0;
364 double dell, ldell;
365 unsigned n, i;
366 const HelmholtzPowTerm *pt;
367 const HelmholtzGausTerm *gt;
368
369
370 #ifdef RESID_DEBUG
371 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
372 #endif
373
374 /* power terms */
375 n = data->np;
376 pt = &(data->pt[0]);
377 dell = ipow(delta,pt->l);
378 ldell = pt->l * dell;
379 unsigned oldl;
380 for(i=0; i<n; ++i){
381 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
382 oldl = pt->l;
383 ++pt;
384 if(i+1==n || oldl != pt->l){
385 if(oldl == 0){
386 res += sum;
387 }else{
388 res += sum * exp(-dell);
389 }
390 sum = 0;
391 dell = ipow(delta,pt->l);
392 ldell = pt->l*dell;
393 }
394 }
395
396 /* gaussian terms */
397 n = data->ng;
398 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
399 gt = &(data->gt[0]);
400 for(i=0; i<n; ++i){
401 #ifdef RESID_DEBUG
402 fprintf(stderr,"i = %d, GAUSSIAN, n = %e, t = %f, d = %f, alpha = %f, beta = %f, gamma = %f, epsilon = %f\n",i+1, gt->n, gt->t, gt->d, gt->alpha, gt->beta, gt->gamma, gt->epsilon);
403 #endif
404 double val2;
405 val2 = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
406 * (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
407 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
408 res += val2;
409 #ifdef RESID_DEBUG
410 fprintf(stderr,"val2 = %f --> res = %f\n",val2,res);
411 #endif
412 ++gt;
413 }
414
415 return res;
416 }
417
418 /**
419 Derivative of the helmholtz residual function with respect to
420 tau.
421 */
422 double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
423
424 double sum;
425 double res = 0;
426 double delX;
427 unsigned l;
428 unsigned n, i;
429 const HelmholtzPowTerm *pt;
430 const HelmholtzGausTerm *gt;
431
432 n = data->np;
433 pt = &(data->pt[0]);
434
435 delX = 1;
436
437 l = 0;
438 sum = 0;
439 for(i=0; i<n; ++i){
440 if(pt->t){
441 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
442 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
443 }
444 ++pt;
445 //fprintf(stderr,"l = %d\n",l);
446 if(i+1==n || l != pt->l){
447 if(l==0){
448 //fprintf(stderr,"Adding non-exp term\n");
449 res += sum;
450 }else{
451 //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
452 res += sum * exp(-delX);
453 }
454 /* set l to new value */
455 if(i+1!=n){
456 l = pt->l;
457 //fprintf(stderr,"New l = %d\n",l);
458 delX = ipow(delta,l);
459 sum = 0;
460 }
461 }
462 }
463
464 //#define RESID_DEBUG
465 /* gaussian terms */
466 n = data->ng;
467 gt = &(data->gt[0]);
468 for(i=0; i<n; ++i){
469 #ifdef RESID_DEBUG
470 fprintf(stderr,"i = %d, GAUSSIAN, n = %e, t = %f, d = %f, alpha = %f, beta = %f, gamma = %f, epsilon = %f\n",i+1, gt->n, gt->t, gt->d, gt->alpha, gt->beta, gt->gamma, gt->epsilon);
471 #endif
472
473 double val2;
474 val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
475 * (2. * gt->beta * tau * (tau - gt->gamma) - gt->t)
476 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
477 res += val2;
478 #ifdef RESID_DEBUG
479 fprintf(stderr,"res = %f\n",res);
480 #endif
481
482 ++gt;
483 }
484
485 return res;
486 }
487
488
489
490 /**
491 Mixed derivative of the helmholtz residual function with respect to
492 delta and tau.
493 */
494 double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
495 double dell,ldell, term, sum = 0, res = 0;
496 unsigned n, i;
497 const HelmholtzPowTerm *pt;
498 const HelmholtzGausTerm *gt;
499
500 /* power terms */
501 n = data->np;
502 pt = &(data->pt[0]);
503 dell = ipow(delta,pt->l);
504 ldell = pt->l * dell;
505 unsigned oldl;
506 for(i=0; i<n; ++i){
507 sum += pt->a * pt->t * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * (pt->d - ldell);
508 oldl = pt->l;
509 ++pt;
510 if(i+1==n || oldl != pt->l){
511 if(oldl == 0){
512 res += sum;
513 }else{
514 res += sum * exp(-dell);
515 }
516 sum = 0;
517 dell = ipow(delta,pt->l);
518 ldell = pt->l*dell;
519 }
520 }
521
522 #ifdef TEST
523 assert(!isinf(res));
524 #endif
525
526 /* gaussian terms */
527 n = data->ng;
528 gt = &(data->gt[0]);
529 for(i=0; i<n; ++i){
530 #ifdef RESID_DEBUG
531 fprintf(stderr,"i = %d, GAUSSIAN, n = %e, t = %f, d = %f, alpha = %f, beta = %f, gamma = %f, epsilon = %f\n",i+1, gt->n, gt->t, gt->d, gt->alpha, gt->beta, gt->gamma, gt->epsilon);
532 #endif
533 double d1 = delta - gt->epsilon;
534 double t1 = tau - gt->gamma;
535 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
536
537 double f1 = gt->t - 2*gt->beta*tau*(tau - gt->gamma);
538 double g1 = gt->d - 2*gt->alpha*delta*(delta - gt->epsilon);
539
540 sum = gt->n * f1 * pow(tau,gt->t-1) * g1 * pow(delta,gt->d-1) * exp(e1);
541
542 //fprintf(stderr,"sum = %f\n",sum);
543 res += sum;
544 #ifdef TEST
545 assert(!isinf(res));
546 #endif
547 ++gt;
548 }
549
550 #ifdef RESID_DEBUG
551 fprintf(stderr,"phir = %f\n",res);
552 #endif
553
554 #ifdef TEST
555 assert(!isnan(res));
556 assert(!isinf(res));
557 #endif
558 return res;
559 }
560
561 /**
562 Second derivative of helmholtz residual function with respect to
563 delta (twice).
564
565 FIXME this function is WRONG.
566 */
567 double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
568 double sum = 0, res = 0;
569 double dell, ldell;
570 unsigned n, i;
571 const HelmholtzPowTerm *pt;
572 const HelmholtzGausTerm *gt;
573
574
575 #ifdef RESID_DEBUG
576 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
577 #endif
578
579 /* power terms */
580 n = data->np;
581 pt = &(data->pt[0]);
582 dell = ipow(delta,pt->l);
583 ldell = pt->l * dell;
584 unsigned oldl;
585 for(i=0; i<n; ++i){
586 double lpart = pt->l ? SQ(ldell) + ldell*(1. - 2*pt->d - pt->l) : 0;
587 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (pt->d*(pt->d - 1) + lpart);
588 oldl = pt->l;
589 ++pt;
590 if(i+1==n || oldl != pt->l){
591 if(oldl == 0){
592 res += sum;
593 }else{
594 res += sum * exp(-dell);
595 }
596 sum = 0;
597 dell = ipow(delta,pt->l);
598 ldell = pt->l*dell;
599 }
600 }
601
602 /* gaussian terms */
603 n = data->ng;
604 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
605 gt = &(data->gt[0]);
606 for(i=0; i<n; ++i){
607 double s1 = SQ(delta - gt->epsilon);
608 double f1 = gt->d*(gt->d - 1)
609 + 2.*gt->alpha*delta * (
610 delta * (2. * gt->alpha * s1 - 1)
611 - 2. * gt->d * (delta - gt->epsilon)
612 );
613 res += gt->n * pow(tau,gt->t) * pow(delta, gt->d - 2.)
614 * f1
615 * exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
616 ++gt;
617 }
618
619 return res;
620 }
621

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22