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

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