/[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 1840 - (show annotations) (download) (as text)
Fri Aug 29 03:52:04 2008 UTC (13 years, 10 months ago) by jpye
File MIME type: text/x-csrc
File size: 12520 byte(s)
Fixed test data for ammonia p,h,s: agreement to 0.09% for all values.
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 #ifdef TEST
61 assert(data->rho_star!=0);
62 assert(T!=0);
63 assert(!isnan(tau));
64 assert(!isnan(delta));
65 assert(!isnan(data->R));
66
67 //fprintf(stderr,"p calc: T = %f\n",T);
68 //fprintf(stderr,"p calc: tau = %f\n",tau);
69 //fprintf(stderr,"p calc: rho = %f\n",rho);
70 //fprintf(stderr,"p calc: delta = %f\n",delta);
71 //fprintf(stderr,"p calc: R*T*rho = %f\n",data->R * T * rho);
72 #endif
73
74 return data->R * T * rho * (1. + delta * helm_resid_del(tau,delta,data));
75 }
76
77 /**
78 Function to calculate internal energy from Helmholtz free energy EOS, given
79 temperature and mass density.
80
81 @param T temperature in K
82 @param rho mass density in kg/m続
83 @return internal energy in ???
84 */
85 double helmholtz_u(double T, double rho, const HelmholtzData *data){
86
87 double tau = data->T_star / T;
88 double delta = rho / data->rho_star;
89
90 #ifdef TEST
91 assert(data->rho_star!=0);
92 assert(T!=0);
93 assert(!isnan(tau));
94 assert(!isnan(delta));
95 assert(!isnan(data->R));
96 #endif
97
98 #ifdef TEST
99 fprintf(stderr,"ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
100 fprintf(stderr,"resid_tau = %f\n",helm_resid_tau(tau,delta,data));
101 fprintf(stderr,"R T = %f\n",data->R * data->T_star);
102 #endif
103
104 return data->R * data->T_star * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data));
105 }
106
107 /**
108 Function to calculate enthalpy from Helmholtz free energy EOS, given
109 temperature and mass density.
110
111 @param T temperature in K
112 @param rho mass density in kg/m続
113 @return enthalpy in J/kg
114 */
115 double helmholtz_h(double T, double rho, const HelmholtzData *data){
116
117 double tau = data->T_star / T;
118 double delta = rho / data->rho_star;
119
120 #ifdef TEST
121 assert(data->rho_star!=0);
122 assert(T!=0);
123 assert(!isnan(tau));
124 assert(!isnan(delta));
125 assert(!isnan(data->R));
126 #endif
127
128 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));
129 }
130
131 /**
132 Function to calculate entropy from Helmholtz free energy EOS, given
133 temperature and mass density.
134
135 @param T temperature in K
136 @param rho mass density in kg/m続
137 @return entropy in J/kgK
138 */
139 double helmholtz_s(double T, double rho, const HelmholtzData *data){
140
141 double tau = data->T_star / T;
142 double delta = rho / data->rho_star;
143
144 #ifdef TEST
145 assert(data->rho_star!=0);
146 assert(T!=0);
147 assert(!isnan(tau));
148 assert(!isnan(delta));
149 assert(!isnan(data->R));
150 #endif
151
152 return data->R * (
153 tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data))
154 - helm_ideal(tau,delta,data->ideal) - helm_resid(tau,delta,data)
155 );
156 }
157
158 /*---------------------------------------------
159 UTILITY FUNCTION(S)
160 */
161
162 /* ipow: public domain by Mark Stephen with suggestions by Keiichi Nakasato */
163 static double ipow(double x, int n){
164 double t = 1.0;
165
166 if(!n)return 1.0; /* At the top. x^0 = 1 */
167
168 if (n < 0){
169 n = -n;
170 x = 1.0/x; /* error if x == 0. Good */
171 } /* ZTC/SC returns inf, which is even better */
172
173 if (x == 0.0)return 0.0;
174
175 do{
176 if(n & 1)t *= x;
177 n /= 2; /* KN prefers if (n/=2) x*=x; This avoids an */
178 x *= x; /* unnecessary but benign multiplication on */
179 }while(n); /* the last pass, but the comparison is always
180 true _except_ on the last pass. */
181
182 return t;
183 }
184
185 /*---------------------------------------------
186 IDEAL COMPONENT RELATIONS
187 */
188
189 /**
190 Ideal component of helmholtz function
191 */
192 double helm_ideal(double tau, double delta, const IdealData *data){
193
194 const IdealPowTerm *pt;
195 const IdealExpTerm *et;
196
197 unsigned i;
198 double sum = log(delta) + data->c + data->m * tau - log(tau);
199 double term;
200
201 //fprintf(stderr,"constant = %f, linear = %f", data->c, data->m);
202 //fprintf(stderr,"initial terms = %f\n",sum);
203 pt = &(data->pt[0]);
204 for(i = 0; i<data->np; ++i, ++pt){
205 term = pt->a0 * pow(tau, pt->t0);
206 //fprintf(stderr,"i = %d: a0 = %f, t0 = %f, term = %f\n",i,pt->a0, pt->t0, term);
207 sum += pt->a0 * pow(tau, pt->t0);
208 }
209
210 #if 0
211 et = &(data->et[0]);
212 for(i=0; i<data->ne; ++i, ++et){
213 sum += et->b * log( 1 - exp(- et->B * tau));
214 }
215 #endif
216
217 return sum;
218 }
219
220 /**
221 Partial dervivative of ideal component of helmholtz residual function with
222 respect to tau.
223 */
224 double helm_ideal_tau(double tau, double delta, const IdealData *data){
225 const IdealPowTerm *pt;
226 const IdealExpTerm *et;
227
228 unsigned i;
229 double sum = -1./tau + data->m;
230
231 pt = &(data->pt[0]);
232 for(i = 0; i<data->np; ++i, ++pt){
233 sum += pt->a0 * pt->t0 * pow(tau, pt->t0 - 1);
234 }
235
236 #if 0
237 et = &(data->et[0]);
238 for(i=0; i<data->ne; ++i, ++et){
239 sum += et->b * log( 1 - exp(- et->B * tau));
240 }
241 #endif
242 return sum;
243 }
244
245 /**
246 Residual part of helmholtz function. Note: we have NOT prematurely
247 optimised here ;-)
248 */
249 double helm_resid(double tau, double delta, const HelmholtzData *data){
250 double sum;
251 double res = 0;
252 double delX;
253 unsigned l;
254 unsigned np, i;
255 const HelmholtzPowTerm *pt;
256
257 np = data->np;
258 pt = &(data->pt[0]);
259
260 delX = 1;
261
262 l = 0;
263 sum = 0;
264 for(i=0; i<np; ++i){
265 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
266 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
267 ++pt;
268 //fprintf(stderr,"l = %d\n",l);
269 if(i+1==np || l != pt->l){
270 if(l==0){
271 //fprintf(stderr,"Adding non-exp term\n");
272 res += sum;
273 }else{
274 //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
275 res += sum * exp(-delX);
276 }
277 /* set l to new value */
278 if(i+1!=np){
279 l = pt->l;
280 //fprintf(stderr,"New l = %d\n",l);
281 delX = ipow(delta,l);
282 sum = 0;
283 }
284 }
285 }
286
287 return res;
288 }
289
290 /**
291 Derivative of the helmholtz residual function with respect to
292 delta.
293
294 NOTE: POSSIBLY STILL AN ERROR IN THIS FUNCTION.
295 */
296 double helm_resid_del(double tau,double delta, const HelmholtzData *data){
297 double sum;
298 double res = 0;
299 double delX, XdelX;
300 unsigned l;
301 unsigned n, i;
302 const HelmholtzPowTerm *pt;
303 const HelmholtzExpTerm *et;
304
305 n = data->np;
306 pt = &(data->pt[0]);
307
308 delX = 1;
309
310 l = 0;
311 sum = 0;
312 XdelX = 0;
313 for(i=0; i<n; ++i){
314 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
315 sum += pt->a * pow(tau, pt->t) * pow(delta, pt->d - 1) * (pt->d - XdelX);
316 ++pt;
317 //fprintf(stderr,"l = %d\n",l);
318 if(i+1==n || l != pt->l){
319 if(l==0){
320 //fprintf(stderr,"Adding non-exp term\n");
321 //fprintf(stderr,"sum = %f\n",sum);
322 res += sum;
323 }else{
324 //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
325 //fprintf(stderr,"sum = %f\n",sum);
326 res += sum * exp(-delX);
327 }
328 /* set l to new value */
329 if(i+1!=n){
330 l = pt->l;
331 delX = ipow(delta,l);
332 XdelX = l * delX;
333 //fprintf(stderr,"New l = %d, XdelX = %f\n",l,XdelX);
334 sum = 0;
335 }
336 }
337 }
338
339 #if 0
340 /* now the exponential terms */
341 n = data->ne;
342 et = &(data->et[0]);
343 for(i=0; i< n; ++i){
344 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);
345
346 double del2 = delta*delta;
347 double tau2 = tau*tau;
348 double gam2 = et->gamma * et->gamma;
349 sum = -et->a * pow(tau,et->t) * ipow(delta,et->d-1)
350 * (2 * et->phi * del2 - 2 * et->phi * delta - et->d)
351 * exp(-et->phi * del2
352 + 2 * et->phi * delta
353 - et->beta * tau2
354 + 2 * et->beta * et->gamma * tau
355 - et->phi
356 - et->beta * gam2
357 );
358 fprintf(stderr,"sum = %f\n",sum);
359 res += sum;
360 ++et;
361 }
362 #endif
363
364 return res;
365 }
366
367 /**
368 Derivative of the helmholtz residual function with respect to
369 tau.
370 */
371 double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
372
373 double sum;
374 double res = 0;
375 double delX;
376 unsigned l;
377 unsigned np, i;
378 const HelmholtzPowTerm *pt;
379
380 np = data->np;
381 pt = &(data->pt[0]);
382
383 delX = 1;
384
385 l = 0;
386 sum = 0;
387 for(i=0; i<np; ++i){
388 if(pt->t){
389 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
390 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
391 }
392 ++pt;
393 //fprintf(stderr,"l = %d\n",l);
394 if(i+1==np || l != pt->l){
395 if(l==0){
396 //fprintf(stderr,"Adding non-exp term\n");
397 res += sum;
398 }else{
399 //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
400 res += sum * exp(-delX);
401 }
402 /* set l to new value */
403 if(i+1!=np){
404 l = pt->l;
405 //fprintf(stderr,"New l = %d\n",l);
406 delX = ipow(delta,l);
407 sum = 0;
408 }
409 }
410 }
411
412 return res;
413 }
414
415
416
417 /**
418 Mixed derivative of the helmholtz residual function with respect to
419 delta and tau
420 */
421 double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
422
423 double sum;
424 double phir = 0;
425 unsigned i;
426 double XdelX;
427
428 const HelmholtzPowTerm *pt = &(data->pt[0]);
429
430 for(i=0; i<5; ++i){
431 phir += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->d * pt->t;
432 ++pt;
433 }
434
435 sum = 0;
436 XdelX = delta;
437 for(i=5; i<10; ++i){
438 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->t *(pt->d - XdelX);
439 ++pt;
440 }
441 phir += exp(-delta) * sum;
442
443 sum = 0;
444 XdelX = 2*delta*delta;
445 for(i=10; i<17; ++i){
446 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->t *(pt->d - XdelX);
447 ++pt;
448 }
449 phir += exp(-delta*delta) * sum;
450
451 sum = 0;
452 XdelX = 3*delta*delta*delta;
453 for(i=17; i<21; ++i){
454 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->t *(pt->d - XdelX);
455 ++pt;
456 }
457 phir += exp(-delta*delta*delta) * sum;
458
459 return phir;
460 }
461
462 #define SQ(X) ((X)*(X))
463
464 /**
465 Second derivative of helmholtz residual function with respect to
466 delta (twice).
467 */
468 double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
469
470 double sum;
471 double phir = 0;
472 unsigned i;
473 unsigned X;
474 double XdelX;
475
476 const HelmholtzPowTerm *pt = &(data->pt[0]);
477
478 for(i=0; i<5; ++i){
479 phir += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (SQ(pt->d) - X);
480 ++pt;
481 }
482
483 sum = 0;
484 X = 1;
485 XdelX = delta;
486 for(i=5; i<10; ++i){
487 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);
488 ++pt;
489 }
490 phir += exp(-delta) * sum;
491
492 sum = 0;
493 X = 2;
494 XdelX = 2*delta*delta;
495 for(i=10; i<17; ++i){
496 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);
497 ++pt;
498 }
499 phir += exp(-delta*delta) * sum;
500
501 sum = 0;
502 X = 3;
503 XdelX = 3*delta*delta*delta;
504 for(i=17; i<21; ++i){
505 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);
506 ++pt;
507 }
508 phir += exp(-delta*delta*delta) * sum;
509
510 return phir;
511 }
512

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