/[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 1857 - (show annotations) (download) (as text)
Thu Sep 11 14:46:51 2008 UTC (13 years, 9 months ago) by jpye
File MIME type: text/x-csrc
File size: 13553 byte(s)
p, u, h working for nitrogen, still problems with s.
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 TEST
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 Calculation zero-pressure specific heat capacity
166 */
167 double helmholtz_cp0(double T, const HelmholtzData *data){
168 double val = helm_cp0(T,data->ideal);
169 #if 0
170 fprintf(stderr,"val = %f\n",val);
171 #endif
172 return val;
173 }
174
175 /*---------------------------------------------
176 UTILITY FUNCTION(S)
177 */
178
179 /* ipow: public domain by Mark Stephen with suggestions by Keiichi Nakasato */
180 static double ipow(double x, int n){
181 double t = 1.0;
182
183 if(!n)return 1.0; /* At the top. x^0 = 1 */
184
185 if (n < 0){
186 n = -n;
187 x = 1.0/x; /* error if x == 0. Good */
188 } /* ZTC/SC returns inf, which is even better */
189
190 if (x == 0.0)return 0.0;
191
192 do{
193 if(n & 1)t *= x;
194 n /= 2; /* KN prefers if (n/=2) x*=x; This avoids an */
195 x *= x; /* unnecessary but benign multiplication on */
196 }while(n); /* the last pass, but the comparison is always
197 true _except_ on the last pass. */
198
199 return t;
200 }
201
202 /**
203 Residual part of helmholtz function.
204 */
205 double helm_resid(double tau, double delta, const HelmholtzData *data){
206 #if 0
207 double sum, res = 0;
208 unsigned n, i;
209 const HelmholtzPowTerm *pt;
210 const HelmholtzExpTerm *et;
211
212 n = data->np;
213 pt = &(data->pt[0]);
214
215 /* power terms */
216 sum = 0;
217 unsigned oldl;
218 for(i=0; i<n; ++i){
219 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
220 oldl = pt->l;
221 ++pt;
222 if(i+1==n || oldl != pt->l){
223 if(oldl == 0){
224 res += sum;
225 }else{
226 res += sum * exp(-ipow(delta,pt->l));
227 }
228 sum = 0;
229 }
230 }
231
232 /* now the exponential terms */
233 n = data->ne;
234 et = &(data->et[0]);
235 for(i=0; i< n; ++i){
236 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);
237
238 double e1 = -et->phi * delta*delta
239 + 2 * et->phi * delta
240 - et->beta * tau * tau
241 + 2 * et->beta * et->gamma * tau
242 - et->phi
243 - et->beta * et->gamma * et->gamma;
244 sum = et->a * pow(tau,et->t) * ipow(delta,et->d) * exp(e1);
245 //fprintf(stderr,"sum = %f\n",sum);
246 res += sum;
247 ++et;
248 }
249
250 return res;
251 }
252
253 #else
254 double sum;
255 double res = 0;
256 double delX;
257 unsigned l;
258 unsigned n, i;
259 const HelmholtzPowTerm *pt;
260 const HelmholtzExpTerm *et;
261
262 n = data->np;
263 pt = &(data->pt[0]);
264
265 delX = 1;
266
267 l = 0;
268 sum = 0;
269 for(i=0; i<n; ++i){
270 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
271 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
272 ++pt;
273 //fprintf(stderr,"l = %d\n",l);
274 if(i+1==n || l != pt->l){
275 if(l==0){
276 //fprintf(stderr,"Adding non-exp term\n");
277 res += sum;
278 }else{
279 //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
280 res += sum * exp(-delX);
281 }
282 /* set l to new value */
283 if(i+1!=n){
284 l = pt->l;
285 //fprintf(stderr,"New l = %d\n",l);
286 delX = ipow(delta,l);
287 sum = 0;
288 }
289 }
290 }
291
292 /* now the exponential terms */
293 n = data->ne;
294 et = &(data->et[0]);
295 for(i=0; i< n; ++i){
296 #ifdef TEST
297 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);
298 #endif
299
300 double e1 = -et->phi * delta*delta
301 + 2 * et->phi * delta
302 - et->beta * tau * tau
303 + 2 * et->beta * et->gamma * tau
304 - et->phi
305 - et->beta * et->gamma * et->gamma;
306 sum = et->a * pow(tau,et->t) * ipow(delta,et->d) * exp(e1);
307 //fprintf(stderr,"sum = %f\n",sum);
308 res += sum;
309 ++et;
310 }
311
312 return res;
313 }
314 #endif
315
316 /**
317 Derivative of the helmholtz residual function with respect to
318 delta.
319 */
320 double helm_resid_del(double tau,double delta, const HelmholtzData *data){
321 double sum, res = 0;
322 double dell, ldell;
323 unsigned n, i;
324 const HelmholtzPowTerm *pt;
325 const HelmholtzExpTerm *et;
326
327 n = data->np;
328 pt = &(data->pt[0]);
329
330 sum = 0;
331 dell = ipow(delta,pt->l);
332 ldell = pt->l * dell;
333 unsigned oldl;
334 for(i=0; i<n; ++i){
335 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
336 oldl = pt->l;
337 ++pt;
338 if(i+1==n || oldl != pt->l){
339 if(oldl == 0){
340 res += sum;
341 }else{
342 res += sum * exp(-dell);
343 }
344 sum = 0;
345 dell = ipow(delta,pt->l);
346 ldell = pt->l*dell;
347 }
348 }
349
350 #if 1
351 /* now the exponential terms */
352 n = data->ne;
353 et = &(data->et[0]);
354 for(i=0; i< n; ++i){
355 //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);
356
357 double del2 = delta*delta;
358 double tau2 = tau*tau;
359 double gam2 = et->gamma * et->gamma;
360 double e1 = -et->phi * del2
361 + 2 * et->phi * delta
362 - et->beta * tau2
363 + 2 * et->beta * et->gamma * tau
364 - et->phi
365 - et->beta * gam2;
366 sum = -et->a * pow(tau,et->t) * ipow(delta,et->d-1)
367 * (2 * et->phi * del2 - 2 * et->phi * delta - et->d)
368 * exp(e1);
369 //fprintf(stderr,"sum = %f\n",sum);
370 res += sum;
371 ++et;
372 }
373 #endif
374
375 return res;
376 }
377
378 /**
379 Derivative of the helmholtz residual function with respect to
380 tau.
381 */
382 double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
383
384 double sum;
385 double res = 0;
386 double delX;
387 unsigned l;
388 unsigned n, i;
389 const HelmholtzPowTerm *pt;
390 const HelmholtzExpTerm *et;
391
392 n = data->np;
393 pt = &(data->pt[0]);
394
395 delX = 1;
396
397 l = 0;
398 sum = 0;
399 for(i=0; i<n; ++i){
400 if(pt->t){
401 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
402 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
403 }
404 ++pt;
405 //fprintf(stderr,"l = %d\n",l);
406 if(i+1==n || l != pt->l){
407 if(l==0){
408 //fprintf(stderr,"Adding non-exp term\n");
409 res += sum;
410 }else{
411 //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
412 res += sum * exp(-delX);
413 }
414 /* set l to new value */
415 if(i+1!=n){
416 l = pt->l;
417 //fprintf(stderr,"New l = %d\n",l);
418 delX = ipow(delta,l);
419 sum = 0;
420 }
421 }
422 }
423
424 #if 1
425 /* now the exponential terms */
426 n = data->ne;
427 et = &(data->et[0]);
428 for(i=0; i< n; ++i){
429 //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);
430
431 double tau2 = tau*tau;
432 double del2 = delta*delta;
433 double gam2 = et->gamma * et->gamma;
434 double e1 = -et->phi * del2
435 + 2 * et->phi * delta
436 - et->beta * tau2
437 + 2 * et->beta * et->gamma * tau
438 - et->phi
439 - et->beta * gam2;
440 sum = -et->a * pow(tau,et->t - 1) * ipow(delta,et->d)
441 * (2 * et->beta * tau2 - 2 * et->beta * et->gamma * tau - et->t)
442 * exp(e1);
443 //fprintf(stderr,"sum = %f\n",sum);
444 res += sum;
445 ++et;
446 }
447 #endif
448
449 return res;
450 }
451
452
453
454 /**
455 Mixed derivative of the helmholtz residual function with respect to
456 delta and tau
457 */
458 double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
459
460 double sum;
461 double phir = 0;
462 unsigned i;
463 double XdelX;
464
465 const HelmholtzPowTerm *pt = &(data->pt[0]);
466
467 for(i=0; i<5; ++i){
468 phir += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->d * pt->t;
469 ++pt;
470 }
471
472 sum = 0;
473 XdelX = delta;
474 for(i=5; i<10; ++i){
475 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->t *(pt->d - XdelX);
476 ++pt;
477 }
478 phir += exp(-delta) * sum;
479
480 sum = 0;
481 XdelX = 2*delta*delta;
482 for(i=10; i<17; ++i){
483 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->t *(pt->d - XdelX);
484 ++pt;
485 }
486 phir += exp(-delta*delta) * sum;
487
488 sum = 0;
489 XdelX = 3*delta*delta*delta;
490 for(i=17; i<21; ++i){
491 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->t *(pt->d - XdelX);
492 ++pt;
493 }
494 phir += exp(-delta*delta*delta) * sum;
495
496 return phir;
497 }
498
499 #define SQ(X) ((X)*(X))
500
501 /**
502 Second derivative of helmholtz residual function with respect to
503 delta (twice).
504 */
505 double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
506
507 double sum;
508 double phir = 0;
509 unsigned i;
510 unsigned X;
511 double XdelX;
512
513 const HelmholtzPowTerm *pt = &(data->pt[0]);
514
515 for(i=0; i<5; ++i){
516 phir += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (SQ(pt->d) - X);
517 ++pt;
518 }
519
520 sum = 0;
521 X = 1;
522 XdelX = delta;
523 for(i=5; i<10; ++i){
524 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);
525 ++pt;
526 }
527 phir += exp(-delta) * sum;
528
529 sum = 0;
530 X = 2;
531 XdelX = 2*delta*delta;
532 for(i=10; i<17; ++i){
533 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);
534 ++pt;
535 }
536 phir += exp(-delta*delta) * sum;
537
538 sum = 0;
539 X = 3;
540 XdelX = 3*delta*delta*delta;
541 for(i=17; i<21; ++i){
542 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);
543 ++pt;
544 }
545 phir += exp(-delta*delta*delta) * sum;
546
547 return phir;
548 }
549

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