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

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