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

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