/[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 1844 - (show annotations) (download) (as text)
Sat Aug 30 07:33:11 2008 UTC (13 years, 10 months ago) by jpye
File MIME type: text/x-csrc
File size: 12213 byte(s)
Cleaned up helmholtz_resid_del, seems to be fixed now.
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 double helm_resid_del(double tau,double delta, const HelmholtzData *data){
298 double sum, res = 0;
299 double dell, ldell;
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 sum = 0;
309 dell = ipow(delta,pt->l);
310 ldell = pt->l * dell;
311 unsigned oldl;
312 for(i=0; i<n; ++i){
313 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
314 oldl = pt->l;
315 ++pt;
316 if(i+1==n || oldl != pt->l){
317 if(oldl == 0){
318 res += sum;
319 }else{
320 res += sum * exp(-dell);
321 }
322 sum = 0;
323 dell = ipow(delta,pt->l);
324 ldell = pt->l*dell;
325 }
326 }
327
328 #if 1
329 /* now the exponential terms */
330 n = data->ne;
331 et = &(data->et[0]);
332 for(i=0; i< n; ++i){
333 //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);
334
335 double del2 = delta*delta;
336 double tau2 = tau*tau;
337 double gam2 = et->gamma * et->gamma;
338 sum = -et->a * pow(tau,et->t) * ipow(delta,et->d-1)
339 * (2 * et->phi * del2 - 2 * et->phi * delta - et->d)
340 * exp(-et->phi * del2
341 + 2 * et->phi * delta
342 - et->beta * tau2
343 + 2 * et->beta * et->gamma * tau
344 - et->phi
345 - et->beta * gam2
346 );
347 //fprintf(stderr,"sum = %f\n",sum);
348 res += sum;
349 ++et;
350 }
351 #endif
352
353 return res;
354 }
355
356 /**
357 Derivative of the helmholtz residual function with respect to
358 tau.
359 */
360 double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
361
362 double sum;
363 double res = 0;
364 double delX;
365 unsigned l;
366 unsigned np, i;
367 const HelmholtzPowTerm *pt;
368
369 np = data->np;
370 pt = &(data->pt[0]);
371
372 delX = 1;
373
374 l = 0;
375 sum = 0;
376 for(i=0; i<np; ++i){
377 if(pt->t){
378 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
379 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
380 }
381 ++pt;
382 //fprintf(stderr,"l = %d\n",l);
383 if(i+1==np || l != pt->l){
384 if(l==0){
385 //fprintf(stderr,"Adding non-exp term\n");
386 res += sum;
387 }else{
388 //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
389 res += sum * exp(-delX);
390 }
391 /* set l to new value */
392 if(i+1!=np){
393 l = pt->l;
394 //fprintf(stderr,"New l = %d\n",l);
395 delX = ipow(delta,l);
396 sum = 0;
397 }
398 }
399 }
400
401 return res;
402 }
403
404
405
406 /**
407 Mixed derivative of the helmholtz residual function with respect to
408 delta and tau
409 */
410 double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
411
412 double sum;
413 double phir = 0;
414 unsigned i;
415 double XdelX;
416
417 const HelmholtzPowTerm *pt = &(data->pt[0]);
418
419 for(i=0; i<5; ++i){
420 phir += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->d * pt->t;
421 ++pt;
422 }
423
424 sum = 0;
425 XdelX = delta;
426 for(i=5; i<10; ++i){
427 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->t *(pt->d - XdelX);
428 ++pt;
429 }
430 phir += exp(-delta) * sum;
431
432 sum = 0;
433 XdelX = 2*delta*delta;
434 for(i=10; i<17; ++i){
435 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->t *(pt->d - XdelX);
436 ++pt;
437 }
438 phir += exp(-delta*delta) * sum;
439
440 sum = 0;
441 XdelX = 3*delta*delta*delta;
442 for(i=17; i<21; ++i){
443 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->t *(pt->d - XdelX);
444 ++pt;
445 }
446 phir += exp(-delta*delta*delta) * sum;
447
448 return phir;
449 }
450
451 #define SQ(X) ((X)*(X))
452
453 /**
454 Second derivative of helmholtz residual function with respect to
455 delta (twice).
456 */
457 double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
458
459 double sum;
460 double phir = 0;
461 unsigned i;
462 unsigned X;
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) * ipow(delta, pt->d - 2) * (SQ(pt->d) - X);
469 ++pt;
470 }
471
472 sum = 0;
473 X = 1;
474 XdelX = delta;
475 for(i=5; i<10; ++i){
476 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);
477 ++pt;
478 }
479 phir += exp(-delta) * sum;
480
481 sum = 0;
482 X = 2;
483 XdelX = 2*delta*delta;
484 for(i=10; i<17; ++i){
485 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);
486 ++pt;
487 }
488 phir += exp(-delta*delta) * sum;
489
490 sum = 0;
491 X = 3;
492 XdelX = 3*delta*delta*delta;
493 for(i=17; i<21; ++i){
494 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);
495 ++pt;
496 }
497 phir += exp(-delta*delta*delta) * sum;
498
499 return phir;
500 }
501

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