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

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