/[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 1836 - (show annotations) (download) (as text)
Mon Aug 25 08:32:38 2008 UTC (11 years, 11 months ago) by jpye
File MIME type: text/x-csrc
File size: 11207 byte(s)
Fixed helm_resid and helm_resid_tau.
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
31 #ifdef TEST
32 #include <assert.h>
33 #include <stdlib.h>
34 #include <stdio.h>
35 #endif
36
37 /* forward decls */
38
39 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 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
47 /**
48 Function to calculate pressure from Helmholtz free energy EOS, given temperature
49 and mass density.
50
51 @param T temperature in K
52 @param rho mass density in kg/m続
53 @return pressure in Pa???
54 */
55 double helmholtz_p(double T, double rho, const HelmholtzData *data){
56
57 double tau = data->T_star / T;
58 double delta = rho / data->rho_star;
59
60 return data->R * T * rho * (1. + delta * helm_resid_del(tau,delta,data));
61 }
62
63 /**
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 double helmholtz_u(double T, double rho, const HelmholtzData *data){
72
73 double tau = data->T_star / T;
74 double delta = rho / data->rho_star;
75
76 #ifdef TEST
77 fprintf(stderr,"ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
78 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 #endif
81
82 return data->R * data->T_star * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data));
83 }
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 @return enthalpy in J/kg
92 */
93 double helmholtz_h(double T, double rho, const HelmholtzData *data){
94
95 double tau = data->T_star / T;
96 double delta = rho / data->rho_star;
97
98 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 }
100
101 /**
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 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 );
118 }
119
120 /*---------------------------------------------
121 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 if(!n)return 1.0; /* At the top. x^0 = 1 */
129
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 IDEAL COMPONENT RELATIONS
149 */
150
151 /**
152 Ideal component of helmholtz function
153 */
154 double helm_ideal(double tau, double delta, const IdealData *data){
155
156 const IdealPowTerm *pt;
157 const IdealExpTerm *et;
158
159 unsigned i;
160 double sum = log(delta) + data->c + data->m * tau - log(tau);
161 double term;
162
163 //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 }
181
182 /**
183 Partial dervivative of ideal component of helmholtz residual function with
184 respect to tau.
185 */
186 double helm_ideal_tau(double tau, double delta, const IdealData *data){
187 const IdealPowTerm *pt;
188 const IdealExpTerm *et;
189
190 unsigned i;
191 double sum = -1./tau + data->m;
192
193 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
198 #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 }
206
207 /**
208 Residual part of helmholtz function. Note: we have NOT prematurely
209 optimised here ;-)
210 */
211 double helm_resid(double tau, double delta, const HelmholtzData *data){
212 double sum;
213 double res = 0;
214 double delX;
215 unsigned l;
216 unsigned nr, i;
217 const HelmholtzATDL *atdl;
218
219 nr = data->nr;
220 atdl = &(data->atdl[0]);
221
222 delX = 1;
223
224 l = 0;
225 sum = 0;
226 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 sum += atdl->a * pow(tau, atdl->t) * ipow(delta, atdl->d);
229 ++atdl;
230 //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 }
248
249 return res;
250 }
251
252 /**
253 Derivative of the helmholtz residual function with respect to
254 delta.
255
256 THERE APPEARS TO BE AN ERROR IN THIS FUNCTION.
257 */
258 double helm_resid_del(double tau,double delta, const HelmholtzData *data){
259 double sum;
260 double res = 0;
261 double delX, XdelX;
262 unsigned l;
263 unsigned nr, i;
264 const HelmholtzATDL *atdl;
265
266 nr = data->nr;
267 atdl = &(data->atdl[0]);
268
269 delX = 1;
270
271 l = 0;
272 sum = 0;
273 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 ++atdl;
278 //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 }
299
300 return res;
301 }
302
303 /**
304 Derivative of the helmholtz residual function with respect to
305 tau.
306 */
307 double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
308
309 double sum;
310 double res = 0;
311 double delX;
312 unsigned l;
313 unsigned nr, i;
314 const HelmholtzATDL *atdl;
315
316 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 return res;
349 }
350
351
352
353 /**
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
364 const HelmholtzATDL *atdl = &(data->atdl[0]);
365
366 for(i=0; i<5; ++i){
367 phir += atdl->a * pow(tau, atdl->t - 1) * ipow(delta, atdl->d - 1) * atdl->d * atdl->t;
368 ++atdl;
369 }
370
371 sum = 0;
372 XdelX = delta;
373 for(i=5; i<10; ++i){
374 sum += atdl->a * pow(tau, atdl->t - 1) * ipow(delta, atdl->d - 1) * atdl->t *(atdl->d - XdelX);
375 ++atdl;
376 }
377 phir += exp(-delta) * sum;
378
379 sum = 0;
380 XdelX = 2*delta*delta;
381 for(i=10; i<17; ++i){
382 sum += atdl->a * pow(tau, atdl->t - 1) * ipow(delta, atdl->d - 1) * atdl->t *(atdl->d - XdelX);
383 ++atdl;
384 }
385 phir += exp(-delta*delta) * sum;
386
387 sum = 0;
388 XdelX = 3*delta*delta*delta;
389 for(i=17; i<21; ++i){
390 sum += atdl->a * pow(tau, atdl->t - 1) * ipow(delta, atdl->d - 1) * atdl->t *(atdl->d - XdelX);
391 ++atdl;
392 }
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 const HelmholtzATDL *atdl = &(data->atdl[0]);
413
414 for(i=0; i<5; ++i){
415 phir += atdl->a * pow(tau, atdl->t) * ipow(delta, atdl->d - 2) * (SQ(atdl->d) - X);
416 ++atdl;
417 }
418
419 sum = 0;
420 X = 1;
421 XdelX = delta;
422 for(i=5; i<10; ++i){
423 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 }
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 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 }
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 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 }
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