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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1836 - (hide annotations) (download) (as text)
Mon Aug 25 08:32:38 2008 UTC (12 years, 1 month ago) by jpye
File MIME type: text/x-csrc
File size: 11207 byte(s)
Fixed helm_resid and helm_resid_tau.
1 jpye 1810 /* ASCEND modelling environment
2 jpye 1831 Copyright (C) 2008 Carnegie Mellon University
3 jpye 1810
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 jpye 1822 For nomenclature see Tillner-Roth, Harms-Watzenberg and Baehr, Eine neue
22     Fundamentalgleichung f端r Ammoniak.
23    
24 jpye 1810 John Pye, 29 Jul 2008.
25     */
26    
27 jpye 1822 #include <math.h>
28 jpye 1810
29 jpye 1822 #include "helmholtz.h"
30    
31 jpye 1826 #ifdef TEST
32     #include <assert.h>
33     #include <stdlib.h>
34     #include <stdio.h>
35     #endif
36    
37 jpye 1822 /* forward decls */
38    
39 jpye 1835 static double helm_ideal(double tau, double delta, const IdealData *data);
40     static double helm_ideal_tau(double tau, double delta, const IdealData *data);
41 jpye 1826 static double helm_resid(double tau, double delta, const HelmholtzData *data);
42     static double helm_resid_del(double tau, double delta, const HelmholtzData *data);
43     static double helm_resid_tau(double tau, double delta, const HelmholtzData *data);
44     static double helm_resid_deltau(double tau, double delta, const HelmholtzData *data);
45     static double helm_resid_deldel(double tau, double delta, const HelmholtzData *data);
46 jpye 1822
47 jpye 1810 /**
48     Function to calculate pressure from Helmholtz free energy EOS, given temperature
49 jpye 1824 and mass density.
50 jpye 1810
51     @param T temperature in K
52 jpye 1824 @param rho mass density in kg/m続
53     @return pressure in Pa???
54 jpye 1810 */
55 jpye 1826 double helmholtz_p(double T, double rho, const HelmholtzData *data){
56 jpye 1822
57     double tau = data->T_star / T;
58     double delta = rho / data->rho_star;
59 jpye 1810
60 jpye 1822 return data->R * T * rho * (1. + delta * helm_resid_del(tau,delta,data));
61     }
62 jpye 1810
63 jpye 1824 /**
64     Function to calculate internal energy from Helmholtz free energy EOS, given
65     temperature and mass density.
66    
67     @param T temperature in K
68     @param rho mass density in kg/m続
69     @return internal energy in ???
70     */
71 jpye 1826 double helmholtz_u(double T, double rho, const HelmholtzData *data){
72 jpye 1824
73     double tau = data->T_star / T;
74     double delta = rho / data->rho_star;
75    
76 jpye 1832 #ifdef TEST
77 jpye 1835 fprintf(stderr,"ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
78 jpye 1827 fprintf(stderr,"resid_tau = %f\n",helm_resid_tau(tau,delta,data));
79     fprintf(stderr,"R T = %f\n",data->R * data->T_star);
80 jpye 1832 #endif
81 jpye 1827
82 jpye 1835 return data->R * data->T_star * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data));
83 jpye 1824 }
84    
85     /**
86     Function to calculate enthalpy from Helmholtz free energy EOS, given
87     temperature and mass density.
88    
89     @param T temperature in K
90     @param rho mass density in kg/m続
91 jpye 1829 @return enthalpy in J/kg
92 jpye 1824 */
93 jpye 1826 double helmholtz_h(double T, double rho, const HelmholtzData *data){
94 jpye 1824
95     double tau = data->T_star / T;
96     double delta = rho / data->rho_star;
97    
98 jpye 1835 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));
99 jpye 1824 }
100    
101 jpye 1829 /**
102     Function to calculate entropy from Helmholtz free energy EOS, given
103     temperature and mass density.
104    
105     @param T temperature in K
106     @param rho mass density in kg/m続
107     @return entropy in J/kgK
108     */
109     double helmholtz_s(double T, double rho, const HelmholtzData *data){
110    
111     double tau = data->T_star / T;
112     double delta = rho / data->rho_star;
113    
114     return data->R * (
115 jpye 1835 tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data))
116     - helm_ideal(tau,delta,data->ideal) - helm_resid(tau,delta,data)
117 jpye 1829 );
118     }
119    
120 jpye 1824 /*---------------------------------------------
121 jpye 1826 UTILITY FUNCTION(S)
122     */
123    
124     /* ipow: public domain by Mark Stephen with suggestions by Keiichi Nakasato */
125     static double ipow(double x, int n){
126     double t = 1.0;
127    
128 jpye 1829 if(!n)return 1.0; /* At the top. x^0 = 1 */
129 jpye 1826
130     if (n < 0){
131     n = -n;
132     x = 1.0/x; /* error if x == 0. Good */
133     } /* ZTC/SC returns inf, which is even better */
134    
135     if (x == 0.0)return 0.0;
136    
137     do{
138     if(n & 1)t *= x;
139     n /= 2; /* KN prefers if (n/=2) x*=x; This avoids an */
140     x *= x; /* unnecessary but benign multiplication on */
141     }while(n); /* the last pass, but the comparison is always
142     true _except_ on the last pass. */
143    
144     return t;
145     }
146    
147     /*---------------------------------------------
148 jpye 1824 IDEAL COMPONENT RELATIONS
149     */
150    
151     /**
152     Ideal component of helmholtz function
153     */
154 jpye 1835 double helm_ideal(double tau, double delta, const IdealData *data){
155 jpye 1810
156 jpye 1835 const IdealPowTerm *pt;
157     const IdealExpTerm *et;
158 jpye 1822
159 jpye 1835 unsigned i;
160     double sum = log(delta) + data->c + data->m * tau - log(tau);
161     double term;
162 jpye 1827
163 jpye 1835 //fprintf(stderr,"constant = %f, linear = %f", data->c, data->m);
164     //fprintf(stderr,"initial terms = %f\n",sum);
165     pt = &(data->pt[0]);
166     for(i = 0; i<data->np; ++i, ++pt){
167     term = pt->a0 * pow(tau, pt->t0);
168     //fprintf(stderr,"i = %d: a0 = %f, t0 = %f, term = %f\n",i,pt->a0, pt->t0, term);
169     sum += pt->a0 * pow(tau, pt->t0);
170     }
171    
172     #if 0
173     et = &(data->et[0]);
174     for(i=0; i<data->ne; ++i, ++et){
175     sum += et->b * log( 1 - exp(- et->B * tau));
176     }
177     #endif
178    
179     return sum;
180 jpye 1822 }
181    
182     /**
183 jpye 1824 Partial dervivative of ideal component of helmholtz residual function with
184     respect to tau.
185     */
186 jpye 1835 double helm_ideal_tau(double tau, double delta, const IdealData *data){
187     const IdealPowTerm *pt;
188     const IdealExpTerm *et;
189 jpye 1824
190 jpye 1835 unsigned i;
191     double sum = -1./tau + data->m;
192 jpye 1824
193 jpye 1835 pt = &(data->pt[0]);
194     for(i = 0; i<data->np; ++i, ++pt){
195     sum += pt->a0 * pt->t0 * pow(tau, pt->t0 - 1);
196     }
197 jpye 1827
198 jpye 1835 #if 0
199     et = &(data->et[0]);
200     for(i=0; i<data->ne; ++i, ++et){
201     sum += et->b * log( 1 - exp(- et->B * tau));
202     }
203     #endif
204     return sum;
205 jpye 1824 }
206    
207     /**
208 jpye 1822 Residual part of helmholtz function. Note: we have NOT prematurely
209     optimised here ;-)
210     */
211 jpye 1826 double helm_resid(double tau, double delta, const HelmholtzData *data){
212 jpye 1822 double sum;
213 jpye 1836 double res = 0;
214     double delX;
215     unsigned l;
216     unsigned nr, i;
217     const HelmholtzATDL *atdl;
218 jpye 1822
219 jpye 1836 nr = data->nr;
220     atdl = &(data->atdl[0]);
221 jpye 1810
222 jpye 1836 delX = 1;
223 jpye 1810
224 jpye 1836 l = 0;
225 jpye 1822 sum = 0;
226 jpye 1836 for(i=0; i<nr; ++i){
227     //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, atdl->a, atdl->t, atdl->d, atdl->l);
228 jpye 1832 sum += atdl->a * pow(tau, atdl->t) * ipow(delta, atdl->d);
229     ++atdl;
230 jpye 1836 //fprintf(stderr,"l = %d\n",l);
231     if(i+1==nr || l != atdl->l){
232     if(l==0){
233     //fprintf(stderr,"Adding non-exp term\n");
234     res += sum;
235     }else{
236     //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
237     res += sum * exp(-delX);
238     }
239     /* set l to new value */
240     if(i+1!=nr){
241     l = atdl->l;
242     //fprintf(stderr,"New l = %d\n",l);
243     delX = ipow(delta,l);
244     sum = 0;
245     }
246     }
247 jpye 1822 }
248    
249 jpye 1836 return res;
250 jpye 1810 }
251    
252 jpye 1822 /**
253     Derivative of the helmholtz residual function with respect to
254     delta.
255 jpye 1830
256 jpye 1836 THERE APPEARS TO BE AN ERROR IN THIS FUNCTION.
257 jpye 1822 */
258 jpye 1826 double helm_resid_del(double tau,double delta, const HelmholtzData *data){
259 jpye 1822 double sum;
260 jpye 1836 double res = 0;
261     double delX, XdelX;
262     unsigned l;
263     unsigned nr, i;
264     const HelmholtzATDL *atdl;
265 jpye 1810
266 jpye 1836 nr = data->nr;
267     atdl = &(data->atdl[0]);
268 jpye 1822
269 jpye 1836 delX = 1;
270 jpye 1822
271 jpye 1836 l = 0;
272 jpye 1822 sum = 0;
273 jpye 1836 XdelX = 0;
274     for(i=0; i<nr; ++i){
275     //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, atdl->a, atdl->t, atdl->d, atdl->l);
276     sum += atdl->a * pow(tau, atdl->t) * ipow(delta, atdl->d - 1) * (atdl->d - XdelX);
277 jpye 1832 ++atdl;
278 jpye 1836 //fprintf(stderr,"l = %d\n",l);
279     if(i+1==nr || l != atdl->l){
280     if(l==0){
281     //fprintf(stderr,"Adding non-exp term\n");
282     //fprintf(stderr,"sum = %f\n",sum);
283     res += sum;
284     }else{
285     //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
286     //fprintf(stderr,"sum = %f\n",sum);
287     res += sum * exp(-delX);
288     }
289     /* set l to new value */
290     if(i+1!=nr){
291     l = atdl->l;
292     delX = ipow(delta,l);
293     XdelX = l * delX;
294     //fprintf(stderr,"New l = %d, XdelX = %f\n",l,XdelX);
295     sum = 0;
296     }
297     }
298 jpye 1822 }
299    
300 jpye 1836 return res;
301 jpye 1822 }
302    
303 jpye 1824 /**
304     Derivative of the helmholtz residual function with respect to
305     tau.
306     */
307 jpye 1826 double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
308 jpye 1824
309     double sum;
310 jpye 1827 double res = 0;
311 jpye 1832 double delX;
312     unsigned l;
313     unsigned nr, i;
314     const HelmholtzATDL *atdl;
315 jpye 1822
316 jpye 1832 nr = data->nr;
317     atdl = &(data->atdl[0]);
318    
319     delX = 1;
320    
321     l = 0;
322     sum = 0;
323     for(i=0; i<nr; ++i){
324     if(atdl->t){
325     //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, atdl->a, atdl->t, atdl->d, atdl->l);
326     sum += atdl->a * pow(tau, atdl->t - 1) * ipow(delta, atdl->d) * atdl->t;
327     }
328     ++atdl;
329     //fprintf(stderr,"l = %d\n",l);
330     if(i+1==nr || l != atdl->l){
331     if(l==0){
332     //fprintf(stderr,"Adding non-exp term\n");
333     res += sum;
334     }else{
335     //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
336     res += sum * exp(-delX);
337     }
338     /* set l to new value */
339     if(i+1!=nr){
340     l = atdl->l;
341     //fprintf(stderr,"New l = %d\n",l);
342     delX = ipow(delta,l);
343     sum = 0;
344     }
345     }
346     }
347    
348 jpye 1827 return res;
349 jpye 1824 }
350    
351    
352    
353 jpye 1826 /**
354     Mixed derivative of the helmholtz residual function with respect to
355     delta and tau
356     */
357     double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
358    
359     double sum;
360     double phir = 0;
361     unsigned i;
362     double XdelX;
363 jpye 1824
364 jpye 1832 const HelmholtzATDL *atdl = &(data->atdl[0]);
365 jpye 1826
366     for(i=0; i<5; ++i){
367 jpye 1832 phir += atdl->a * pow(tau, atdl->t - 1) * ipow(delta, atdl->d - 1) * atdl->d * atdl->t;
368     ++atdl;
369 jpye 1826 }
370    
371     sum = 0;
372     XdelX = delta;
373     for(i=5; i<10; ++i){
374 jpye 1832 sum += atdl->a * pow(tau, atdl->t - 1) * ipow(delta, atdl->d - 1) * atdl->t *(atdl->d - XdelX);
375     ++atdl;
376 jpye 1826 }
377     phir += exp(-delta) * sum;
378    
379     sum = 0;
380     XdelX = 2*delta*delta;
381     for(i=10; i<17; ++i){
382 jpye 1832 sum += atdl->a * pow(tau, atdl->t - 1) * ipow(delta, atdl->d - 1) * atdl->t *(atdl->d - XdelX);
383     ++atdl;
384 jpye 1826 }
385     phir += exp(-delta*delta) * sum;
386    
387     sum = 0;
388     XdelX = 3*delta*delta*delta;
389     for(i=17; i<21; ++i){
390 jpye 1832 sum += atdl->a * pow(tau, atdl->t - 1) * ipow(delta, atdl->d - 1) * atdl->t *(atdl->d - XdelX);
391     ++atdl;
392 jpye 1826 }
393     phir += exp(-delta*delta*delta) * sum;
394    
395     return phir;
396     }
397    
398     #define SQ(X) ((X)*(X))
399    
400     /**
401     Second derivative of helmholtz residual function with respect to
402     delta (twice).
403     */
404     double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
405    
406     double sum;
407     double phir = 0;
408     unsigned i;
409     unsigned X;
410     double XdelX;
411    
412 jpye 1832 const HelmholtzATDL *atdl = &(data->atdl[0]);
413 jpye 1826
414     for(i=0; i<5; ++i){
415 jpye 1832 phir += atdl->a * pow(tau, atdl->t) * ipow(delta, atdl->d - 2) * (SQ(atdl->d) - X);
416     ++atdl;
417 jpye 1826 }
418    
419     sum = 0;
420     X = 1;
421     XdelX = delta;
422     for(i=5; i<10; ++i){
423 jpye 1832 sum += atdl->a * pow(tau, atdl->t) * ipow(delta, atdl->d - 2) * (SQ(XdelX) - X*XdelX - 2*atdl->d*XdelX + XdelX + SQ(atdl->d) - atdl->d);
424     ++atdl;
425 jpye 1826 }
426     phir += exp(-delta) * sum;
427    
428     sum = 0;
429     X = 2;
430     XdelX = 2*delta*delta;
431     for(i=10; i<17; ++i){
432 jpye 1832 sum += atdl->a * pow(tau, atdl->t) * ipow(delta, atdl->d - 2) * (SQ(XdelX) - X*XdelX - 2*atdl->d*XdelX + XdelX + SQ(atdl->d) - atdl->d);
433     ++atdl;
434 jpye 1826 }
435     phir += exp(-delta*delta) * sum;
436    
437     sum = 0;
438     X = 3;
439     XdelX = 3*delta*delta*delta;
440     for(i=17; i<21; ++i){
441 jpye 1832 sum += atdl->a * pow(tau, atdl->t) * ipow(delta, atdl->d - 2) * (SQ(XdelX) - X*XdelX - 2*atdl->d*XdelX + XdelX + SQ(atdl->d) - atdl->d);
442     ++atdl;
443 jpye 1826 }
444     phir += exp(-delta*delta*delta) * sum;
445    
446     return phir;
447     }
448    

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