/[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 1829 - (show annotations) (download) (as text)
Fri Aug 22 05:00:53 2008 UTC (15 years, 9 months ago) by jpye
File MIME type: text/x-csrc
File size: 19726 byte(s)
Added entropy relation and tests. All entropy tests pass to within 0.03%.
1 /* ASCEND modelling environment
2 Copyright (C) 2006 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 /* Property data for Ammonia, from Tillner-Roth, Harms-Watzenberg and
38 Baehr, 'Eine neue Fundamentalgleichung f端r Ammoniak', DKV-Tagungsbericht,
39 20:167-181, 1993. This is the ammmonia property correlation recommended
40 by NIST in its program REFPROP 7.0. */
41 const HelmholtzData helmholtz_data_ammonia = {
42 /* R */ 488.189 /* J/kg/K */
43 , /* rho_star */225. /* kg/m続 */
44 , /* T_star */ 405.40 /* K */
45
46 , {
47 /* a0_1 */ -15.815020
48 ,/* a0_2 */ 4.255726
49 ,/* a0_3 */ 11.474340
50 ,/* a0_4 */ -1.296211
51 ,/* a0_5 */ 0.5706757
52 }
53
54 , {
55 /* a_i, t_i, d_i */
56 /* 1 */{0.4554431E-1, -0.5 , 2}
57 ,{0.7238548E+0, 0.5 , 1 }
58 ,{0.1229470E-1, 1 , 4 }
59 ,{-0.1858814E+1, 1.5 , 1 }
60 /* 5 */,{0.2141882E-10, 3 , 15 }
61 ,{-0.1430020E-1, 0 , 3 }
62 ,{0.3441324E+0, 3 , 3 }
63 ,{-0.2873571E+0, 4 , 1 }
64 ,{0.2352589E-4, 4 , 8 }
65 /* 10 */,{-0.3497111E-1, 5 , 2}
66 ,{0.2397852E-1, 3 , 1}
67 ,{0.1831117E-2, 5 , 8}
68 ,{-0.4085375E-1, 6 , 1}
69 ,{0.2379275E+0, 8 , 2}
70 /* 15 */,{-0.3548972E-1, 8 , 3}
71 ,{-0.1823729E+0, 10, 2}
72 ,{0.2281556E-1, 10 , 4}
73 ,{-0.6663444E-2, 5 , 3}
74 ,{-0.8847486E-2, 7.5, 1}
75 /* 20 */,{0.2272635E-2 , 15 , 2}
76 ,{-0.5588655E-3, 30, 4}
77 }
78 };
79
80 /* forward decls */
81
82 static double helm_ideal(double tau, double delta, const HelmholtzData *data);
83 static double helm_ideal_tau(double tau, double delta, const HelmholtzData *data);
84 static double helm_resid(double tau, double delta, const HelmholtzData *data);
85 static double helm_resid_del(double tau, double delta, const HelmholtzData *data);
86 static double helm_resid_tau(double tau, double delta, const HelmholtzData *data);
87 static double helm_resid_deltau(double tau, double delta, const HelmholtzData *data);
88 static double helm_resid_deldel(double tau, double delta, const HelmholtzData *data);
89
90 /**
91 Function to calculate pressure from Helmholtz free energy EOS, given temperature
92 and mass density.
93
94 @param T temperature in K
95 @param rho mass density in kg/m続
96 @return pressure in Pa???
97 */
98 double helmholtz_p(double T, double rho, const HelmholtzData *data){
99
100 double tau = data->T_star / T;
101 double delta = rho / data->rho_star;
102
103 return data->R * T * rho * (1. + delta * helm_resid_del(tau,delta,data));
104 }
105
106 /**
107 Function to calculate internal energy from Helmholtz free energy EOS, given
108 temperature and mass density.
109
110 @param T temperature in K
111 @param rho mass density in kg/m続
112 @return internal energy in ???
113 */
114 double helmholtz_u(double T, double rho, const HelmholtzData *data){
115
116 double tau = data->T_star / T;
117 double delta = rho / data->rho_star;
118
119 fprintf(stderr,"ideal_tau = %f\n",helm_ideal_tau(tau,delta,data));
120 fprintf(stderr,"resid_tau = %f\n",helm_resid_tau(tau,delta,data));
121
122 fprintf(stderr,"R T = %f\n",data->R * data->T_star);
123
124 return data->R * data->T_star * (helm_ideal_tau(tau,delta,data) + helm_resid_tau(tau,delta,data));
125 }
126
127 /**
128 Function to calculate enthalpy from Helmholtz free energy EOS, given
129 temperature and mass density.
130
131 @param T temperature in K
132 @param rho mass density in kg/m続
133 @return enthalpy in J/kg
134 */
135 double helmholtz_h(double T, double rho, const HelmholtzData *data){
136
137 double tau = data->T_star / T;
138 double delta = rho / data->rho_star;
139
140 return data->R * T * (1 + tau * (helm_ideal_tau(tau,delta,data) + helm_resid_tau(tau,delta,data)) + delta*helm_resid_del(tau,delta,data));
141 }
142
143 /**
144 Function to calculate entropy from Helmholtz free energy EOS, given
145 temperature and mass density.
146
147 @param T temperature in K
148 @param rho mass density in kg/m続
149 @return entropy in J/kgK
150 */
151 double helmholtz_s(double T, double rho, const HelmholtzData *data){
152
153 double tau = data->T_star / T;
154 double delta = rho / data->rho_star;
155
156 return data->R * (
157 tau * (helm_ideal_tau(tau,delta,data) + helm_resid_tau(tau,delta,data))
158 - helm_ideal(tau,delta,data) - helm_resid(tau,delta,data)
159 );
160 }
161
162 /*---------------------------------------------
163 UTILITY FUNCTION(S)
164 */
165
166 /* ipow: public domain by Mark Stephen with suggestions by Keiichi Nakasato */
167 static double ipow(double x, int n){
168 double t = 1.0;
169
170 if(!n)return 1.0; /* At the top. x^0 = 1 */
171
172 if (n < 0){
173 n = -n;
174 x = 1.0/x; /* error if x == 0. Good */
175 } /* ZTC/SC returns inf, which is even better */
176
177 if (x == 0.0)return 0.0;
178
179 do{
180 if(n & 1)t *= x;
181 n /= 2; /* KN prefers if (n/=2) x*=x; This avoids an */
182 x *= x; /* unnecessary but benign multiplication on */
183 }while(n); /* the last pass, but the comparison is always
184 true _except_ on the last pass. */
185
186 return t;
187 }
188
189 /*---------------------------------------------
190 IDEAL COMPONENT RELATIONS
191 */
192
193 /**
194 Ideal component of helmholtz function
195 */
196 double helm_ideal(double tau, double delta, const HelmholtzData *data){
197 const double *a0;
198
199 double tau13 = pow(tau,1./3.);
200 double taum32 = pow(tau,-3./2.);
201 double taum74 = pow(tau,-7./4.);
202
203 a0 = &(data->a0[0]);
204
205 return log(delta) + a0[0] + a0[1] * tau - log(tau) + a0[2] * tau13 + a0[3] * taum32 + a0[4] * taum74;
206 }
207
208 /**
209 Partial dervivative of ideal component of helmholtz residual function with
210 respect to tau.
211 */
212 double helm_ideal_tau(double tau, double delta, const HelmholtzData *data){
213 const double *a0;
214
215 double taum114 = pow(tau,-11./4.);
216 double taum52 = pow(tau,-5./2.);
217 double taum23 = pow(tau,-2./3.);
218
219 //fprintf(stderr,"tau = %f, taum23 = %f\n",tau,taum23);
220 //fprintf(stderr,"tau = %f, taum52 = %f\n",tau,taum52);
221 //fprintf(stderr,"tau = %f, taum114 = %f\n",tau,taum114);
222
223 a0 = &(data->a0[0]);
224
225 //unsigned i; for(i=0;i<5;++i)fprintf(stderr,"a0[%u] = %f\n",i,a0[i]);
226
227 double res;
228 res = a0[2]/3.*taum23 - 1./tau - 3./2.*a0[3]*taum52 - 7./4.*a0[4]*taum114 + a0[1];
229 //fprintf(stderr,"res = %f\n",res);
230 return res;
231 }
232
233 /**
234 Residual part of helmholtz function. Note: we have NOT prematurely
235 optimised here ;-)
236 */
237 double helm_resid(double tau, double delta, const HelmholtzData *data){
238
239 double sum;
240 double phir = 0;
241 unsigned i;
242
243 const HelmholtzATD *atd = &(data->atd[0]);
244
245 for(i=0; i<5; ++i){
246 phir += atd->a * pow(tau, atd->t) * ipow(delta, atd->d);
247 ++atd;
248 }
249
250 sum = 0;
251 for(i=5; i<10; ++i){
252 sum += atd->a * pow(tau, atd->t) * ipow(delta, atd->d);
253 ++atd;
254 }
255 phir += exp(-delta) * sum;
256
257 sum = 0;
258 for(i=10; i<17; ++i){
259 sum += atd->a * pow(tau, atd->t) * ipow(delta, atd->d);
260 ++atd;
261 }
262 phir += exp(-delta*delta) * sum;
263
264 sum = 0;
265 for(i=17; i<21; ++i){
266 sum += atd->a * pow(tau, atd->t) * ipow(delta, atd->d);
267 ++atd;
268 }
269 phir += exp(-delta*delta*delta) * sum;
270
271 return phir;
272 }
273
274 /**
275 Derivative of the helmholtz residual function with respect to
276 delta.
277 */
278 double helm_resid_del(double tau,double delta, const HelmholtzData *data){
279
280 double sum;
281 double phir = 0;
282 unsigned i;
283 double XdelX;
284
285 const HelmholtzATD *atd = &(data->atd[0]);
286
287 for(i=0; i<5; ++i){
288 fprintf(stderr,"i = %d, a = %e, t = %f, d = %d\n",i+1, atd->a, atd->t, atd->d);
289 phir += atd->a * pow(tau, atd->t) * pow(delta, atd->d - 1) * atd->d;
290 ++atd;
291 }
292
293 sum = 0;
294 XdelX = delta;
295 for(i=5; i<10; ++i){
296 fprintf(stderr,"i = %d, a = %e, t = %f, d = %d\n",i+1, atd->a, atd->t, atd->d);
297 sum += atd->a * pow(tau, atd->t) * pow(delta, atd->d - 1) * (atd->d - delta);
298 ++atd;
299 }
300 //fprintf(stderr,"sum = %f\n",sum);
301 phir += exp(-delta) * sum;
302
303 sum = 0;
304 XdelX = 2.*delta*delta;
305 for(i=10; i<17; ++i){
306 fprintf(stderr,"i = %d, a = %e, t = %f, d = %d\n",i+1, atd->a, atd->t, atd->d);
307 sum += atd->a * pow(tau, atd->t) * pow(delta, atd->d - 1) * (atd->d - XdelX);
308 ++atd;
309 }
310 //fprintf(stderr,"sum = %f\n",sum);
311 phir += exp(-delta*delta) * sum;
312
313 sum = 0;
314 XdelX = 3.*delta*delta*delta;
315 for(i=17; i<21; ++i){
316 fprintf(stderr,"i = %d, a = %e, t = %f, d = %d\n",i+1, atd->a, atd->t, atd->d);
317 sum += atd->a * pow(tau, atd->t) * pow(delta, atd->d - 1) * (atd->d - XdelX);
318 ++atd;
319 }
320 //fprintf(stderr,"sum = %f\n",sum);
321 phir += exp(-delta*delta*delta) * sum;
322
323 return phir;
324 }
325
326 /**
327 Derivative of the helmholtz residual function with respect to
328 tau.
329 */
330 double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
331
332 double sum;
333 double res = 0;
334 unsigned i;
335
336 const HelmholtzATD *atd = &(data->atd[0]);
337
338 for(i=0; i<5; ++i){
339 if(atd->t){
340 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d\n",i+1, atd->a, atd->t, atd->d);
341 res += atd->a * pow(tau, atd->t - 1) * ipow(delta, atd->d) * atd->t;
342 }
343 ++atd;
344 }
345
346 sum = 0;
347 for(i=5; i<10; ++i){
348 if(atd->t){
349 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d\n",i+1, atd->a, atd->t, atd->d);
350 sum += atd->a * pow(tau, atd->t - 1) * ipow(delta, atd->d) * atd->t;
351 }
352 ++atd;
353 }
354 res += exp(-delta) * sum;
355
356 sum = 0;
357 for(i=10; i<17; ++i){
358 if(atd->t){
359 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d\n",i+1, atd->a, atd->t, atd->d);
360 sum += atd->a * pow(tau, atd->t - 1) * ipow(delta, atd->d) * atd->t;
361 }
362 ++atd;
363 }
364 res += exp(-delta*delta) * sum;
365
366 sum = 0;
367 for(i=17; i<21; ++i){
368 if(atd->t){
369 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d\n",i+1, atd->a, atd->t, atd->d);
370 sum += atd->a * pow(tau, atd->t - 1) * ipow(delta, atd->d) * atd->t;
371 }
372 ++atd;
373 }
374 res += exp(-delta*delta*delta) * sum;
375
376 return res;
377 }
378
379
380
381 /**
382 Mixed derivative of the helmholtz residual function with respect to
383 delta and tau
384 */
385 double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
386
387 double sum;
388 double phir = 0;
389 unsigned i;
390 double XdelX;
391
392 const HelmholtzATD *atd = &(data->atd[0]);
393
394 for(i=0; i<5; ++i){
395 phir += atd->a * pow(tau, atd->t - 1) * ipow(delta, atd->d - 1) * atd->d * atd->t;
396 ++atd;
397 }
398
399 sum = 0;
400 XdelX = delta;
401 for(i=5; i<10; ++i){
402 sum += atd->a * pow(tau, atd->t - 1) * ipow(delta, atd->d - 1) * atd->t *(atd->d - XdelX);
403 ++atd;
404 }
405 phir += exp(-delta) * sum;
406
407 sum = 0;
408 XdelX = 2*delta*delta;
409 for(i=10; i<17; ++i){
410 sum += atd->a * pow(tau, atd->t - 1) * ipow(delta, atd->d - 1) * atd->t *(atd->d - XdelX);
411 ++atd;
412 }
413 phir += exp(-delta*delta) * sum;
414
415 sum = 0;
416 XdelX = 3*delta*delta*delta;
417 for(i=17; i<21; ++i){
418 sum += atd->a * pow(tau, atd->t - 1) * ipow(delta, atd->d - 1) * atd->t *(atd->d - XdelX);
419 ++atd;
420 }
421 phir += exp(-delta*delta*delta) * sum;
422
423 return phir;
424 }
425
426 #define SQ(X) ((X)*(X))
427
428 /**
429 Second derivative of helmholtz residual function with respect to
430 delta (twice).
431 */
432 double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
433
434 double sum;
435 double phir = 0;
436 unsigned i;
437 unsigned X;
438 double XdelX;
439
440 const HelmholtzATD *atd = &(data->atd[0]);
441
442 for(i=0; i<5; ++i){
443 phir += atd->a * pow(tau, atd->t) * ipow(delta, atd->d - 2) * (SQ(atd->d) - X);
444 ++atd;
445 }
446
447 sum = 0;
448 X = 1;
449 XdelX = delta;
450 for(i=5; i<10; ++i){
451 sum += atd->a * pow(tau, atd->t) * ipow(delta, atd->d - 2) * (SQ(XdelX) - X*XdelX - 2*atd->d*XdelX + XdelX + SQ(atd->d) - atd->d);
452 ++atd;
453 }
454 phir += exp(-delta) * sum;
455
456 sum = 0;
457 X = 2;
458 XdelX = 2*delta*delta;
459 for(i=10; i<17; ++i){
460 sum += atd->a * pow(tau, atd->t) * ipow(delta, atd->d - 2) * (SQ(XdelX) - X*XdelX - 2*atd->d*XdelX + XdelX + SQ(atd->d) - atd->d);
461 ++atd;
462 }
463 phir += exp(-delta*delta) * sum;
464
465 sum = 0;
466 X = 3;
467 XdelX = 3*delta*delta*delta;
468 for(i=17; i<21; ++i){
469 sum += atd->a * pow(tau, atd->t) * ipow(delta, atd->d - 2) * (SQ(XdelX) - X*XdelX - 2*atd->d*XdelX + XdelX + SQ(atd->d) - atd->d);
470 ++atd;
471 }
472 phir += exp(-delta*delta*delta) * sum;
473
474 return phir;
475 }
476
477
478 /*
479 Test suite. These tests attempt to validate the current code using
480 a few sample figures output by REFPROP 7.0.
481
482 To run the test, compile and run as follows:
483
484 gcc helmholtz.c -DTEST -o helmholtz -lm && ./helmholtz
485 */
486 #ifdef TEST
487
488 /* a simple macro to actually do the testing */
489 #define ASSERT_TOL(EXPR,VAL,TOL) {\
490 double cval; cval = (EXPR);\
491 double err; err = cval - (double)(VAL);\
492 double relerrpc = (cval-(VAL))/(VAL)*100;\
493 if(fabs(err)>TOL){\
494 fprintf(stderr,"ERROR in line %d: value of '%s' = %f, should be %f, error is %f (%.2f%%)!\n"\
495 , __LINE__, #EXPR, cval, VAL,cval-(VAL),relerrpc);\
496 exit(1);\
497 }else{\
498 fprintf(stderr," OK, %s = %8.2e with %.2f%% err.\n",#EXPR,VAL,relerrpc);\
499 /*fprintf(stderr," (err = %8.2e, tol = %8.2e, calc = %8.2e)\n",fabs(err),TOL,cval);*/\
500 }\
501 }
502
503 int main(void){
504 double rho, T, p, h, u;
505 const HelmholtzData *d;
506
507 d = &helmholtz_data_ammonia;
508
509
510 #if 0
511 fprintf(stderr,"PRESSURE TESTS\n");
512
513 fprintf(stderr,"p(T,rho) = 0.1 MPa\n");
514 ASSERT_TOL(helmholtz_p(273.15+-70,724.75,d), 0.1E6, 7E3);
515 ASSERT_TOL(helmholtz_p(273.15+-60,713.65,d), 0.1E6, 7E3);
516 ASSERT_TOL(helmholtz_p(273.15+-50,702.11,d), 0.1E6, 7E3);
517 ASSERT_TOL(helmholtz_p(273.15+-40,690.16,d), 0.1E6, 7E3);
518 ASSERT_TOL(helmholtz_p(273.15+-33.588,682.29,d), 0.1E6, 7E3);
519 ASSERT_TOL(helmholtz_p(273.15+ 0,0.74124,d), 0.1E6, 5E3);
520 ASSERT_TOL(helmholtz_p(273.15+100,0.55135,d), 0.1E6, 5E3);
521 ASSERT_TOL(helmholtz_p(273.15+250,0.39203,d), 0.1E6, 5E3);
522 ASSERT_TOL(helmholtz_p(273.15+420,0.29562,d), 0.1E6, 5E3);
523 fprintf(stderr,"p(T,rho) = 1 MPa\n");
524 ASSERT_TOL(helmholtz_p(273.15+-70,725.06,d), 1E6, 1E3);
525 ASSERT_TOL(helmholtz_p(273.15+ 0,638.97,d), 1E6, 1E3);
526 ASSERT_TOL(helmholtz_p(273.15+ 30,7.5736,d), 1E6, 1E3);
527 ASSERT_TOL(helmholtz_p(273.15+150,4.9817,d), 1E6, 1E3);
528 ASSERT_TOL(helmholtz_p(273.15+200,4.4115,d), 1E6, 1E3);
529 ASSERT_TOL(helmholtz_p(273.15+350,3.3082,d), 1E6, 1E3);
530 ASSERT_TOL(helmholtz_p(273.15+420,2.9670,d), 1E6, 1E3);
531
532 fprintf(stderr,"p(T,rho) = 10 MPa\n");
533 //ASSERT_TOL(helmholtz_p(273.15+-40.,694.67,d), 10E6, 1E3);
534 //ASSERT_TOL(helmholtz_p(273.15+-20.,670.55,d), 10E6, 1E3);
535 //ASSERT_TOL(helmholtz_p(273.15+50,573.07,d), 10E6, 1E3);
536 //ASSERT_TOL(helmholtz_p(273.15+110,441.77,d), 10E6, 1E3);
537 ASSERT_TOL(helmholtz_p(273.15+150,74.732,d), 10E6, 1E3);
538 ASSERT_TOL(helmholtz_p(273.15+200,54.389,d), 10E6, 1E3);
539 ASSERT_TOL(helmholtz_p(273.15+350,35.072,d), 10E6, 1E3);
540 ASSERT_TOL(helmholtz_p(273.15+420,30.731,d), 10E6, 1E3);
541
542 fprintf(stderr,"p(T,rho) = 20 MPa\n");
543 ASSERT_TOL(helmholtz_p(273.15+150,359.41,d), 20E6, 1E4);
544 ASSERT_TOL(helmholtz_p(273.15+200,152.83,d), 20E6, 1E4);
545 ASSERT_TOL(helmholtz_p(273.15+350,74.590,d), 20E6, 1E4);
546 ASSERT_TOL(helmholtz_p(273.15+420,63.602,d), 20E6, 1E4);
547
548 //fprintf(stderr,"IDEAL HELMHOLTZ COMPONENT\n");
549 //ASSERT_TOL(helm_ideal(273.15, 0)
550
551 /* this offset is required to attain agreement with values from REFPROP */
552 double Z = -1635.7e3 + 1492.411e3;
553
554 fprintf(stderr,"ENTHALPY TESTS\n");
555 fprintf(stderr,"h(T,rho) at p = 0.1 MPa\n");
556 ASSERT_TOL(helmholtz_h(273.15+-60, 713.65,d), Z+75.166e3, 0.2e3);
557 ASSERT_TOL(helmholtz_h(273.15+ 0,0.76124,d), Z+1635.7e3, 0.2e3);
558 ASSERT_TOL(helmholtz_h(273.15+ 50,0.63869,d), Z+1744.0e3, 0.2e3);
559 ASSERT_TOL(helmholtz_h(273.15+200,0.43370,d), Z+2087.0e3, 0.2e3);
560 ASSERT_TOL(helmholtz_h(273.15+300,0.35769,d), Z+2340.0e3, 1e3);
561 ASSERT_TOL(helmholtz_h(273.15+420,0.29562,d), Z+2674.3e3, 1e3);
562
563 fprintf(stderr,"h(T,rho) at p = 1 MPa\n");
564 ASSERT_TOL(helmholtz_h(273.15+150,4.9817,d), Z+1949.1e3, 1e3);
565 ASSERT_TOL(helmholtz_h(273.15+200,4.4115,d), Z+2072.7e3, 1e3);
566 ASSERT_TOL(helmholtz_h(273.15+350,3.3082,d), Z+2468.2e3, 1e3);
567 ASSERT_TOL(helmholtz_h(273.15+420,2.9670,d), Z+2668.6e3, 1e3);
568
569 fprintf(stderr,"h(T,rho) at p = 10 MPa\n");
570 ASSERT_TOL(helmholtz_h(273.15+-50,706.21,d), Z+127.39e3, 2e3);
571 ASSERT_TOL(helmholtz_h(273.15+-0,645.04,d), Z+349.53e3, 2e3);
572
573 ASSERT_TOL(helmholtz_h(273.15+150,74.732,d), Z+1688.5e3, 1e3);
574 ASSERT_TOL(helmholtz_h(273.15+200,54.389,d), Z+1908.0e3, 1e3);
575 ASSERT_TOL(helmholtz_h(273.15+350,35.072,d), Z+2393.4e3, 1e3);
576 ASSERT_TOL(helmholtz_h(273.15+420,30.731,d), Z+2611.8e3, 1e3);
577
578 fprintf(stderr,"h(T,rho) at p = 20 MPa\n");
579 ASSERT_TOL(helmholtz_h(273.15+-50,710.19,d), Z+136.54e3, 0.1e3);
580 ASSERT_TOL(helmholtz_h(273.15+ 30,612.22,d), Z+493.28e3, 0.1e3);
581 ASSERT_TOL(helmholtz_h(273.15+150,359.41,d), Z+1162.5e3, 0.1e3);
582 ASSERT_TOL(helmholtz_h(273.15+200,152.83,d), Z+1662.9e3, 0.1e3);
583 ASSERT_TOL(helmholtz_h(273.15+350,74.590,d), Z+2393.4e3, 1e3);
584 ASSERT_TOL(helmholtz_h(273.15+420,63.602,d), Z+2611.8e3, 1e3);
585
586 fprintf(stderr,"h(T,rho) at p = 100 MPa\n");
587 ASSERT_TOL(helmholtz_h(273.15+ 0,690.41,d), Z+422.69e3, 0.5e3);
588 ASSERT_TOL(helmholtz_h(273.15+100,591.07,d), Z+850.44e3, 0.1e3);
589 ASSERT_TOL(helmholtz_h(273.15+250,437.69,d), Z+1506.6e3, 1e3);
590 ASSERT_TOL(helmholtz_h(273.15+420,298.79,d), Z+2252.3e3, 1e3);
591
592 #endif
593
594 fprintf(stderr,"ENTROPY TESTS\n");
595
596 /* offset required to attain agreement with REFPROP */
597 double Y = -471.596704;
598
599 fprintf(stderr,"s(T,rho) at p = 0.1 MPa\n");
600 ASSERT_TOL(helmholtz_s(273.15+-60, 713.65,d), Y+0.36737e3, 0.5);
601 ASSERT_TOL(helmholtz_s(273.15+ 0,0.76124,d), Y+6.8900e3, 0.5);
602 ASSERT_TOL(helmholtz_s(273.15+ 50,0.63869,d), Y+7.2544e3, 0.5);
603 ASSERT_TOL(helmholtz_s(273.15+200,0.43370,d), Y+8.1232e3, 0.5);
604 ASSERT_TOL(helmholtz_s(273.15+300,0.35769,d), Y+8.6084e3, 1);
605 ASSERT_TOL(helmholtz_s(273.15+420,0.29562,d), Y+9.1365e3, 1);
606
607 fprintf(stderr,"s(T,rho) at p = 1 MPa\n");
608 ASSERT_TOL(helmholtz_s(273.15+-50,702.49,d), Y+0.56381e3, 0.5);
609 ASSERT_TOL(helmholtz_s(273.15+150,4.9817,d), Y+6.7008e3, 0.5);
610 ASSERT_TOL(helmholtz_s(273.15+200,4.4115,d), Y+6.9770e3, 0.5);
611 ASSERT_TOL(helmholtz_s(273.15+350,3.3082,d), Y+7.7012e3, 0.5);
612 ASSERT_TOL(helmholtz_s(273.15+420,2.9670,d), Y+8.0059e3, 0.5);
613
614 fprintf(stderr,"s(T,rho) at p = 10 MPa\n");
615 ASSERT_TOL(helmholtz_s(273.15+-70,728.11,d), Y+0.14196e3, 1);
616 ASSERT_TOL(helmholtz_s(273.15+-50,706.21,d), Y+0.54289e3, 1);
617 ASSERT_TOL(helmholtz_s(273.15+-20,670.55,d), Y+1.0975e3, 1);
618 ASSERT_TOL(helmholtz_s(273.15+ 0,645.04,d), Y+1.4403e3, 1);
619 ASSERT_TOL(helmholtz_s(273.15+125.17,356.70,d), Y+3.5463e3, 1);
620
621 ASSERT_TOL(helmholtz_s(273.15+125.17,121.58,d), Y+4.5150e3, 1);
622 ASSERT_TOL(helmholtz_s(273.15+200,54.389,d), Y+5.5906e3, 1);
623 ASSERT_TOL(helmholtz_s(273.15+350,35.072,d), Y+6.4850e3, 1);
624 ASSERT_TOL(helmholtz_s(273.15+420,30.731,d), Y+6.8171e3, 1);
625
626 fprintf(stderr,"s(T,rho) at p = 20 MPa\n");
627 ASSERT_TOL(helmholtz_s(273.15+-50,710.19,d), Y+0.52061e3, 1);
628 ASSERT_TOL(helmholtz_s(273.15+ 30,612.22,d), Y+1.8844e3, 1);
629 ASSERT_TOL(helmholtz_s(273.15+150,359.41,d), Y+3.7164e3, 1);
630 ASSERT_TOL(helmholtz_s(273.15+200,152.83,d), Y+4.8376e3, 1);
631 ASSERT_TOL(helmholtz_s(273.15+350,74.590,d), Y+6.0407e3, 1);
632 ASSERT_TOL(helmholtz_s(273.15+420,63.602,d), Y+6.4066e3, 1);
633
634 fprintf(stderr,"s(T,rho) at p = 100 MPa\n");
635 ASSERT_TOL(helmholtz_s(273.15+ 0,690.41,d), Y+1.2158e3, 1);
636 ASSERT_TOL(helmholtz_s(273.15+100,591.07,d), Y+2.5499e3, 1);
637 ASSERT_TOL(helmholtz_s(273.15+250,437.69,d), Y+4.0264e3, 1);
638 ASSERT_TOL(helmholtz_s(273.15+420,298.79,d), Y+5.2620e3, 1);
639
640
641 fprintf(stderr,"Tests completed OK\n");
642 exit(0);
643 }
644 #endif

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