/[ascend]/trunk/models/johnpye/fprops/water.c
ViewVC logotype

Contents of /trunk/models/johnpye/fprops/water.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1996 - (show annotations) (download) (as text)
Thu Feb 5 09:42:42 2009 UTC (13 years, 4 months ago) by jpye
File MIME type: text/x-csrc
File size: 11364 byte(s)
Added acentric factor to HelmholtzData.
Working on adding calculation of saturation curve using Maxwell phase-equilibrium condition (per IAPWS95).
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 */
19
20 #include "water.h"
21
22 #define WATER_R 461.51805 /* J/kg路K */
23 #define WATER_TSTAR 647.096 /* K */
24
25 /**
26 Ideal gas data for Water/Steam, from IAPWS-95.
27 http://www.iapws.org/relguide/IAPWS95.pdf
28 */
29 const IdealData ideal_data_water = {
30 -8.32044648201 /* constant */
31 , 6.6832105268 /* linear */
32 , WATER_TSTAR /* Tstar */
33 , WATER_R /* cpstar J/kgK */
34 , 1 /* power terms */
35 , (const IdealPowTerm[]){
36 {1. + 3.00632, 0}
37 }
38 , 5 /* exponential terms */
39 , (const IdealExpTerm []){
40 /* b, beta */
41 {0.012436, 1.28728967*WATER_TSTAR}
42 ,{0.97315, 3.53734222*WATER_TSTAR}
43 ,{1.27950, 7.74073708*WATER_TSTAR}
44 ,{0.96956, 9.24437796*WATER_TSTAR}
45 ,{0.24873, 27.5075105*WATER_TSTAR}
46 }
47 };
48
49
50 /**
51 Residual (non-ideal) property data for Water/Steam, from IAPWS-95.
52 http://www.iapws.org/relguide/IAPWS95.pdf
53 */
54 const HelmholtzData helmholtz_data_water = {
55 /* R */ WATER_R /* J/kg/K */
56 , /* M */ 18.015242 /* kg/kmol -- G S Kell, J Phys Chem Ref Data (6) 1109 (1977) */
57 , /* rho_star */322. /* kg/m鲁 */
58 , /* T_star */ WATER_TSTAR /* K */
59 , 0.344 /* acentric factor, source: Reid, Prausnitz & Polling */
60 , &ideal_data_water
61 , 51 /* np */
62 , (const HelmholtzPowTerm[]){
63 /* a_i, t_i, d_i, l_i */
64 {0.12533547935523E-1, -0.5, 1, 0}
65 ,{0.78957634722828E1, 0.875, 1, 0}
66 ,{-0.87803203303561E1, 1, 1, 0}
67 ,{0.31802509345418, 0.5, 2, 0}
68 ,{-0.26145533859358, 0.75, 2, 0}
69 ,{-0.78199751687981E-2, 0.375, 3, 0}
70 ,{0.88089493102134E-2, 1, 4, 0}
71 ,{-0.66856572307965, 4, 1, 1}
72 ,{0.20433810950965, 6, 1, 1}
73 ,{-0.66212605039687E-4, 12, 1, 1}
74 ,{-0.19232721156002, 1, 2, 1}
75 ,{-0.25709043003438, 5, 2, 1}
76 ,{0.16074868486251, 4, 3, 1}
77 ,{-0.40092828925807E-1, 2, 4, 1}
78 ,{0.39343422603254E-6, 13, 4, 1}
79 ,{-0.75941377088144E-5, 9, 5, 1}
80 ,{0.56250979351888E-3, 3, 7, 1}
81 ,{-0.15608652257135E-4, 4, 9, 1}
82 ,{0.11537996422951E-8, 11, 10, 1}
83 ,{0.36582165144204E-6, 4, 11, 1}
84 ,{-0.13251180074668E-11, 13, 13, 1}
85 ,{-0.62639586912454E-9, 1, 15, 1}
86 ,{-0.10793600908932, 7, 1, 2}
87 ,{0.17611491008752E-1, 1, 2, 2}
88 ,{0.22132295167546, 9, 2, 2}
89 ,{-0.40247669763528, 10, 2, 2}
90 ,{0.58083399985759, 10, 3, 2}
91 ,{0.49969146990806E-2, 3, 4, 2}
92 ,{-0.31358700712549E-1, 7, 4, 2}
93 ,{-0.74315929710341, 10, 4, 2}
94 ,{0.47807329915480, 10, 5, 2}
95 ,{0.20527940895948E-1, 6, 6, 2}
96 ,{-0.13636435110343, 10, 6, 2}
97 ,{0.14180634400617E-1, 10, 7, 2}
98 ,{0.83326504880713E-2, 1, 9, 2}
99 ,{-0.29052336009585E-1, 2, 9, 2}
100 ,{0.38615085574206E-1, 3, 9, 2}
101 ,{-0.20393486513704E-1, 4, 9, 2}
102 ,{-0.16554050063734E-2, 8, 9, 2}
103 ,{0.19955571979541E-2, 6, 10, 2}
104 ,{0.15870308324157E-3, 9, 10, 2}
105 ,{-0.16388568342530E-4, 8, 12, 2}
106 ,{0.43613615723811E-1, 16, 3, 3}
107 ,{0.34994005463765E-1, 22, 4, 3}
108 ,{-0.76788197844621E-1, 23, 4, 3}
109 ,{0.22446277332006E-1, 23, 5, 3}
110 ,{-0.62689710414685E-4, 10, 14, 4}
111 ,{-0.55711118565645E-9, 50, 3, 6}
112 ,{-0.19905718354408, 44, 6, 6}
113 ,{0.31777497330738, 46, 6, 6}
114 ,{-0.11841182425981, 50, 6, 6}
115 }
116 , 3 /* gaussian terms */
117 , (const HelmholtzGausTerm[]){
118 /* n, t, d, alpha, beta, gamma, epsilon */
119 {-0.31306260323435e2, 0, 3, 20, 150, 1.21, 1}
120 ,{0.31546140237781e2, 1, 3, 20, 150, 1.21, 1}
121 ,{-0.25213154341695e4,4, 3, 20, 250, 1.25, 1}
122 }
123 , 2 /* critical terms */
124 , (const HelmholtzCritTerm[]){
125 /* n, a, b, beta, A, B, C, D */
126 {-0.14874640856724, 3.5, 0.85, 0.3, 0.32, 0.2, 28, 700}
127 ,{0.31806110878444, 3.5, 0.95, 0.3, 0.32, 0.2, 32, 800}
128 }
129 };
130
131 #ifdef TEST
132 /*
133 Test suite. These tests attempt to validate the current code using
134 a few sample figures output by REFPROP 7.0.
135
136 To run the test, compile and run as follows:
137
138 ./test.py water
139 */
140
141 /*
142 some code from freesteam, http://freesteam.sf.net/, which has been thoroughly
143 validated already.
144 */
145
146 const double n0[] = {
147 0.0 /* placeholder */,
148 -8.32044648201, 6.6832105268, 3.00632 /* const, linear, ln(tau) coeffs */
149 , 0.012436, 0.97315, 1.27950, 0.96956, 0.24873 /* exponential coeffs */
150 };
151
152 const double gamma0[] = {
153 0.0, 0.0, 0.0, 0.0,
154 1.28728967,
155 3.53734222,
156 7.74073708,
157 9.24437796,
158 27.5075105
159 };
160
161 enum Limits{
162 eGamma1 = 4,
163 eGamma2 = 9,
164 };
165
166 #include <assert.h>
167 #include <stdlib.h>
168 #include <stdio.h>
169 #include <math.h>
170 #include "ideal_impl.h"
171 #include "helmholtz_impl.h"
172
173 double phi0(const double delta, const double tau){
174 int i;
175 double sum = 0;
176 for (i = eGamma1; i < eGamma2; i++)
177 {
178 sum += n0[i]*log(1-exp(-tau*gamma0[i]));
179 }
180 sum += log(delta) + n0[1] + n0[2]*tau + n0[3]*log(tau);
181 return sum;
182 }
183
184 typedef struct{double T, rho, p, cv, w, s;} TestDataIAPWS95;
185 const TestDataIAPWS95 td[]; const unsigned ntd;
186
187 typedef struct{double T, p, rho_f, rho_g, h_f, h_g, s_f, s_g;} TestDataSat;
188 const TestDataSat tds[]; const unsigned ntds;
189
190 int main(void){
191 double rho, T;
192 const HelmholtzData *d;
193
194 d = &helmholtz_data_water;
195 double maxerr = 0;
196 unsigned i;
197
198 /* a simple macro to actually do the testing */
199 #define ASSERT_TOL(FN,PARAM1,PARAM2,PARAM3,VAL,TOL) {\
200 double cval; cval = FN(PARAM1,PARAM2,PARAM3);\
201 double err; err = cval - (double)(VAL);\
202 double relerrpc = (cval-(VAL))/(VAL)*100;\
203 if(fabs(relerrpc)>maxerr)maxerr=fabs(relerrpc);\
204 if(fabs(err)>fabs(TOL)){\
205 fprintf(stderr,"ERROR in line %d: value of '%s(%f,%f,%s)' = %0.8f,"\
206 " should be %f, error is %.10e (%.7f%%)!\n"\
207 , __LINE__, #FN,PARAM1,PARAM2,#PARAM3, cval, VAL,cval-(VAL)\
208 ,relerrpc\
209 );\
210 /*exit(1);*/\
211 }else{\
212 fprintf(stderr," OK, %s(%f,%f,%s) = %8.2e with %.8f%% err.\n"\
213 ,#FN,PARAM1,PARAM2,#PARAM3,VAL,relerrpc\
214 );\
215 }\
216 }
217
218 #if 0
219 /* these tests pass, but don't prove much */
220 fprintf(stderr,"COMPARISON OF phi0 VALUES WITH THOSE FROM FREESTEAM\n");
221 for(T = 300; T <= 900; T+= 100){
222 for(rho = 900; rho >= 0.9; rho*=0.5){
223 double delta = rho / d->rho_star;
224 double tau = d->T_star / T;
225 double p0 = phi0(delta,tau);
226
227 ASSERT_TOL(helm_ideal, tau, delta, d->ideal, p0, p0*1e-5);
228 }
229 }
230 #endif
231
232 /* LOW-LEVEL TEST DATA PROVIDED IN IAPWS95 */
233
234 fprintf(stderr,"\nIAPWS95 TABLE 6 TESTS\n");
235 T = 500.; /* K */
236 rho = 838.025; /* kg/m鲁 */
237 double tau = d->T_star / T;
238 double delta = rho / d->rho_star;
239
240 ASSERT_TOL(helm_ideal, tau, delta, d->ideal, 0.204797733E1, 1e-8);
241 ASSERT_TOL(helm_ideal_tau, tau, delta, d->ideal, 0.904611106E1, 1e-8);
242 ASSERT_TOL(HELM_IDEAL_DELTAU, tau, delta, d->ideal, 0., 1e-8);
243
244 double phitt = helm_ideal_tautau(tau, d->ideal);
245 double val = (-0.193249185E1);
246 double err = phitt - val;
247 if(fabs(err) > 1e-8){
248 fprintf(stderr,"ERROR in helm_ideal_tautau near line %d\n",__LINE__);
249 exit(1);
250 }else{
251 fprintf(stderr," OK, helm_ideal_tautau(%f) = %8.2e with %.6f%% err.\n"
252 ,tau,val,err/val*100.
253 );
254 }
255
256 /* FIXME still need to implement helm_ideal_del, helm_ideal_deldel, helm_ideal_deltau */
257
258 ASSERT_TOL(helm_resid, tau, delta, d, -0.342693206E1, 1e-8);
259 ASSERT_TOL(helm_resid_del, tau, delta, d, -0.364366650, 1e-8);
260 ASSERT_TOL(helm_resid_deldel, tau, delta, d, 0.856063701, 1e-8);
261 ASSERT_TOL(helm_resid_tau, tau, delta, d, -0.581403435E1, 1e-8);
262 ASSERT_TOL(helm_resid_tautau, tau, delta, d, -0.223440737E1, 1e-8);
263 ASSERT_TOL(helm_resid_deltau, tau, delta, d, -0.112176915e1, 1e-8);
264
265 #if 0
266 fprintf(stderr,"\nADDITIONAL LOW-LEVEL TESTS NEAR CRITICAL POINT\n");
267
268 T = 647.; /* K */
269 rho = 358.; /* kg/m鲁 */
270 tau = d->T_star / T;
271 delta = rho / d->rho_star;
272
273 /* this test value calculated from pressure using REFPROP 8 */
274 ASSERT_TOL(helmholtz_a, T, rho, d, -8.286875181e5, 1e-4);
275 ASSERT_TOL(helm_resid_del, tau, delta, d, -7.14012024e-1, 1e-8);
276 ASSERT_TOL(helmholtz_s, T, rho, d, 4.320923066e3, 5e-8);
277 ASSERT_TOL(helmholtz_cv, T, rho, d, 6.183157277e3, 5e-7);
278 ASSERT_TOL(helmholtz_p, T, rho, d, 2.203847557e7, 7e-4);
279 ASSERT_TOL(helmholtz_cp, T, rho, d, 3.531798573e6, 1e-8);
280 ASSERT_TOL(helmholtz_w, T, rho, d, 2.52140783e2, 1e-8);
281 #endif
282
283 fprintf(stderr,"\nIAPWS95 TABLE 7 (SINGLE-PHASE) TESTS\n");
284 for(i=0; i<ntd; ++i){
285 double T = td[i].T;
286 double rho = td[i].rho;
287 double p = td[i].p * 1e6; /* Pa */
288 double cv = td[i].cv * 1e3; /* J/kgK */
289 double w = td[i].w; /* m/s */
290 double s = td[i].s * 1e3; /* J/kgK */
291 //fprintf(stderr,"T = %f, rho = %f, p = %f, w = %f, wcalc = %f\n",T,rho,p,w, helmholtz_w(T,rho,d));
292 ASSERT_TOL(helmholtz_s, T, rho, d, s, s*1e-8);
293 ASSERT_TOL(helmholtz_p, T, rho, d, p, p*1e-8);
294 ASSERT_TOL(helmholtz_cv, T, rho, d, cv, cv*1e-8);
295 ASSERT_TOL(helmholtz_w, T, rho, d, w, w*2e-5);
296 }
297
298 #ifdef NOT_YET_FDA_APPROVED
299 fprintf(stderr,"\nIAPWS95 TABLE 8 (SATURATION) TESTS\n");
300 for(i=0; i<ntds; ++i){
301 double T = tds[i].T;
302 double p = tds[i].p * 1e6; /* Pa */
303 double rho_f = tds[i].rho_f;
304 double rho_g = tds[i].rho_g;
305 double h_f = tds[i].h_f * 1e3;
306 double h_g = tds[i].h_g * 1e3;
307 double s_f = tds[i].s_f * 1e3;
308 double s_g = tds[i].s_g * 1e3;
309 fprintf(stderr,"T = %f, p = %f, rho_f = %f, rho_g = %f\n",T,p,rho_f, rho_g);
310 double rho_f_eval, rho_g_eval, p_eval;
311
312 int res;
313 res = helmholtz_sat_t(T, &p, &rho_f, &rho_g, d);
314 }
315 #endif
316
317 fprintf(stderr,"Tests completed OK (maximum error = %0.8f%%)\n",maxerr);
318 exit(0);
319 }
320
321 /* HIGHER-LEVEL TEST-DATA PROVIDED IN IAPWS95 */
322
323 const TestDataIAPWS95 td[] = {
324 {300, 0.9965560e3, 0.992418352e-1, 0.413018112e1, 0.150151914e4, 0.393062643}
325 ,{300, 0.1005308e4, 0.200022515e2, 0.406798347e1, 0.153492501e4, 0.387405401}
326 ,{300, 0.1188202e4, 0.700004704e3, 0.346135580e1, 0.244357992e4, 0.132609616}
327 ,{500, 0.4350000, 0.999679423e-1, 0.150817541e1, 0.548314253e3, 0.794488271e1}
328 ,{500, 0.4532000e1, 0.999938125, 0.166991025e1, 0.535739001e3, 0.682502725e1}
329 ,{500, 0.8380250e3, 0.100003858e2, 0.322106219e1, 0.127128441e4, 0.256690919e1}
330 ,{500, 0.1084564e4, 0.700000405e3, 0.307437693e1, 0.241200877e4, 0.203237509e1}
331 ,{647, 0.3580000e3, 0.220384756e2, 0.618315728e1, 0.252145078e3, 0.432092307e1}
332 ,{900, 0.2410000, 0.100062559, 0.175890657e1, 0.724027147e3, 0.916653194e1}
333 ,{900, 0.5261500e2, 0.200000690e2, 0.193510526e1, 0.698445674e3, 0.659070225e1}
334 ,{900, 0.8707690e3, 0.700000006e3, 0.266422350e1, 0.201933608e4, 0.417223802e1}
335 };
336
337 const unsigned ntd = sizeof(td)/sizeof(TestDataIAPWS95);
338
339 const TestDataSat tds[] = {
340 /* T, p (MPa), rho_f, rho_g, h_f (kJ/kg), h_g (kJ/kg), s_f (kJ/kgK), s_g (kJ/kgK) */
341 {450, 0.932203564, 0.890341250e3, 0.481200360e1, 0.749161585e3, 0.277441078e4, 0.210865845e1, 0.660921221e1}
342 ,{275, 0.698451167e-3, 0.999887406e3, 0.550664919e-2, 0.775972202e1, 0.250428995e4, 0.283094670e-1, 0.910660121e1}
343 ,{625, 0.169082693e2, 0.567090385e3, 0.118290280e3, 0.168626976e4, 0.255071625e4, 0.380194683e1, 0.518506121e1}
344 };
345
346 const unsigned ntds = sizeof(tds)/sizeof(TestDataSat);
347
348 #endif
349

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