/[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 1994 - (show annotations) (download) (as text)
Wed Feb 4 03:52:36 2009 UTC (11 years, 6 months ago) by jpye
File MIME type: text/x-csrc
File size: 10089 byte(s)
Still some error near critical point for calculation of speed of sound.
However this error is not reconcilable with both IAPWS95 and REFPROP8, because those
two sources give different values. The value returned by FPROPS is midway between
the two.
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 */ 0.00000000000000 /* kg/kmol -- need to look up value cited by IAPWS */
57 , /* rho_star */322. /* kg/m鲁 */
58 , /* T_star */ WATER_TSTAR /* K */
59 , &ideal_data_water
60 , 51 /* np */
61 , (const HelmholtzPowTerm[]){
62 /* a_i, t_i, d_i, l_i */
63 {0.12533547935523E-1, -0.5, 1, 0}
64 ,{0.78957634722828E1, 0.875, 1, 0}
65 ,{-0.87803203303561E1, 1, 1, 0}
66 ,{0.31802509345418, 0.5, 2, 0}
67 ,{-0.26145533859358, 0.75, 2, 0}
68 ,{-0.78199751687981E-2, 0.375, 3, 0}
69 ,{0.88089493102134E-2, 1, 4, 0}
70 ,{-0.66856572307965, 4, 1, 1}
71 ,{0.20433810950965, 6, 1, 1}
72 ,{-0.66212605039687E-4, 12, 1, 1}
73 ,{-0.19232721156002, 1, 2, 1}
74 ,{-0.25709043003438, 5, 2, 1}
75 ,{0.16074868486251, 4, 3, 1}
76 ,{-0.40092828925807E-1, 2, 4, 1}
77 ,{0.39343422603254E-6, 13, 4, 1}
78 ,{-0.75941377088144E-5, 9, 5, 1}
79 ,{0.56250979351888E-3, 3, 7, 1}
80 ,{-0.15608652257135E-4, 4, 9, 1}
81 ,{0.11537996422951E-8, 11, 10, 1}
82 ,{0.36582165144204E-6, 4, 11, 1}
83 ,{-0.13251180074668E-11, 13, 13, 1}
84 ,{-0.62639586912454E-9, 1, 15, 1}
85 ,{-0.10793600908932, 7, 1, 2}
86 ,{0.17611491008752E-1, 1, 2, 2}
87 ,{0.22132295167546, 9, 2, 2}
88 ,{-0.40247669763528, 10, 2, 2}
89 ,{0.58083399985759, 10, 3, 2}
90 ,{0.49969146990806E-2, 3, 4, 2}
91 ,{-0.31358700712549E-1, 7, 4, 2}
92 ,{-0.74315929710341, 10, 4, 2}
93 ,{0.47807329915480, 10, 5, 2}
94 ,{0.20527940895948E-1, 6, 6, 2}
95 ,{-0.13636435110343, 10, 6, 2}
96 ,{0.14180634400617E-1, 10, 7, 2}
97 ,{0.83326504880713E-2, 1, 9, 2}
98 ,{-0.29052336009585E-1, 2, 9, 2}
99 ,{0.38615085574206E-1, 3, 9, 2}
100 ,{-0.20393486513704E-1, 4, 9, 2}
101 ,{-0.16554050063734E-2, 8, 9, 2}
102 ,{0.19955571979541E-2, 6, 10, 2}
103 ,{0.15870308324157E-3, 9, 10, 2}
104 ,{-0.16388568342530E-4, 8, 12, 2}
105 ,{0.43613615723811E-1, 16, 3, 3}
106 ,{0.34994005463765E-1, 22, 4, 3}
107 ,{-0.76788197844621E-1, 23, 4, 3}
108 ,{0.22446277332006E-1, 23, 5, 3}
109 ,{-0.62689710414685E-4, 10, 14, 4}
110 ,{-0.55711118565645E-9, 50, 3, 6}
111 ,{-0.19905718354408, 44, 6, 6}
112 ,{0.31777497330738, 46, 6, 6}
113 ,{-0.11841182425981, 50, 6, 6}
114 }
115 , 3 /* gaussian terms */
116 , (const HelmholtzGausTerm[]){
117 /* n, t, d, alpha, beta, gamma, epsilon */
118 {-0.31306260323435e2, 0, 3, 20, 150, 1.21, 1}
119 ,{0.31546140237781e2, 1, 3, 20, 150, 1.21, 1}
120 ,{-0.25213154341695e4,4, 3, 20, 250, 1.25, 1}
121 }
122 , 2 /* critical terms */
123 , (const HelmholtzCritTerm[]){
124 /* n, a, b, beta, A, B, C, D */
125 {-0.14874640856724, 3.5, 0.85, 0.3, 0.32, 0.2, 28, 700}
126 ,{0.31806110878444, 3.5, 0.95, 0.3, 0.32, 0.2, 32, 800}
127 }
128 };
129
130 #ifdef TEST
131 /*
132 Test suite. These tests attempt to validate the current code using
133 a few sample figures output by REFPROP 7.0.
134
135 To run the test, compile and run as follows:
136
137 ./test.py water
138 */
139
140 /*
141 some code from freesteam, http://freesteam.sf.net/, which has been thoroughly
142 validated already.
143 */
144
145 const double n0[] = {
146 0.0 /* placeholder */,
147 -8.32044648201, 6.6832105268, 3.00632 /* const, linear, ln(tau) coeffs */
148 , 0.012436, 0.97315, 1.27950, 0.96956, 0.24873 /* exponential coeffs */
149 };
150
151 const double gamma0[] = {
152 0.0, 0.0, 0.0, 0.0,
153 1.28728967,
154 3.53734222,
155 7.74073708,
156 9.24437796,
157 27.5075105
158 };
159
160 enum Limits{
161 eGamma1 = 4,
162 eGamma2 = 9,
163 };
164
165 #include <assert.h>
166 #include <stdlib.h>
167 #include <stdio.h>
168 #include <math.h>
169 #include "ideal_impl.h"
170 #include "helmholtz_impl.h"
171
172 double phi0(const double delta, const double tau){
173 int i;
174 double sum = 0;
175 for (i = eGamma1; i < eGamma2; i++)
176 {
177 sum += n0[i]*log(1-exp(-tau*gamma0[i]));
178 }
179 sum += log(delta) + n0[1] + n0[2]*tau + n0[3]*log(tau);
180 return sum;
181 }
182
183 typedef struct{double T, rho, p, cv, w, s;} TestDataIAPWS95;
184 const TestDataIAPWS95 td[]; const unsigned ntd;
185
186 int main(void){
187 double rho, T;
188 const HelmholtzData *d;
189
190 d = &helmholtz_data_water;
191 double maxerr = 0;
192 unsigned i;
193
194 /* a simple macro to actually do the testing */
195 #define ASSERT_TOL(FN,PARAM1,PARAM2,PARAM3,VAL,TOL) {\
196 double cval; cval = FN(PARAM1,PARAM2,PARAM3);\
197 double err; err = cval - (double)(VAL);\
198 double relerrpc = (cval-(VAL))/(VAL)*100;\
199 if(fabs(relerrpc)>maxerr)maxerr=fabs(relerrpc);\
200 if(fabs(err)>fabs(TOL)){\
201 fprintf(stderr,"ERROR in line %d: value of '%s(%f,%f,%s)' = %0.8f,"\
202 " should be %f, error is %.10e (%.2f%%)!\n"\
203 , __LINE__, #FN,PARAM1,PARAM2,#PARAM3, cval, VAL,cval-(VAL)\
204 ,relerrpc\
205 );\
206 /*exit(1);*/\
207 }else{\
208 fprintf(stderr," OK, %s(%f,%f,%s) = %8.2e with %.8f%% err.\n"\
209 ,#FN,PARAM1,PARAM2,#PARAM3,VAL,relerrpc\
210 );\
211 }\
212 }
213
214
215 #if 0
216 /* these tests pass, but don't prove much */
217 fprintf(stderr,"COMPARISON OF phi0 VALUES WITH THOSE FROM FREESTEAM\n");
218 for(T = 300; T <= 900; T+= 100){
219 for(rho = 900; rho >= 0.9; rho*=0.5){
220 double delta = rho / d->rho_star;
221 double tau = d->T_star / T;
222 double p0 = phi0(delta,tau);
223
224 ASSERT_TOL(helm_ideal, tau, delta, d->ideal, p0, p0*1e-5);
225 }
226 }
227 #endif
228
229 /* LOW-LEVEL TEST DATA PROVIDED IN IAPWS95 */
230
231 fprintf(stderr,"\nIAPWS95 TABLE 6 TESTS\n");
232 T = 500.; /* K */
233 rho = 838.025; /* kg/m鲁 */
234 double tau = d->T_star / T;
235 double delta = rho / d->rho_star;
236
237 ASSERT_TOL(helm_ideal, tau, delta, d->ideal, 0.204797733E1, 1e-8);
238 ASSERT_TOL(helm_ideal_tau, tau, delta, d->ideal, 0.904611106E1, 1e-8);
239 ASSERT_TOL(HELM_IDEAL_DELTAU, tau, delta, d->ideal, 0., 1e-8);
240
241 double phitt = helm_ideal_tautau(tau, d->ideal);
242 double val = (-0.193249185E1);
243 double err = phitt - val;
244 if(fabs(err) > 1e-8){
245 fprintf(stderr,"ERROR in helm_ideal_tautau near line %d\n",__LINE__);
246 exit(1);
247 }else{
248 fprintf(stderr," OK, helm_ideal_tautau(%f) = %8.2e with %.6f%% err.\n"
249 ,tau,val,err/val*100.
250 );
251 }
252
253 /* FIXME still need to implement helm_ideal_del, helm_ideal_deldel, helm_ideal_deltau */
254
255 ASSERT_TOL(helm_resid, tau, delta, d, -0.342693206E1, 1e-8);
256 ASSERT_TOL(helm_resid_del, tau, delta, d, -0.364366650, 1e-8);
257 ASSERT_TOL(helm_resid_deldel, tau, delta, d, 0.856063701, 1e-8);
258 ASSERT_TOL(helm_resid_tau, tau, delta, d, -0.581403435E1, 1e-8);
259 ASSERT_TOL(helm_resid_tautau, tau, delta, d, -0.223440737E1, 1e-8);
260 ASSERT_TOL(helm_resid_deltau, tau, delta, d, -0.112176915e1, 1e-8);
261
262 fprintf(stderr,"\nADDITIONAL LOW-LEVEL TESTS NEAR CRITICAL POINT\n");
263
264 T = 647.; /* K */
265 rho = 358.; /* kg/m鲁 */
266 tau = d->T_star / T;
267 delta = rho / d->rho_star;
268
269 /* this test value calculated from pressure using REFPROP 8 */
270 ASSERT_TOL(helmholtz_a, T, rho, d, -8.286875181e5, 1e-4);
271 ASSERT_TOL(helm_resid_del, tau, delta, d, -7.14012024e-1, 1e-8);
272 ASSERT_TOL(helmholtz_s, T, rho, d, 4.320923066e3, 5e-8);
273 ASSERT_TOL(helmholtz_cv, T, rho, d, 6.183157277e3, 5e-7);
274 ASSERT_TOL(helmholtz_p, T, rho, d, 2.203847557e7, 7e-4);
275 ASSERT_TOL(helmholtz_cp, T, rho, d, 3.531798573e6, 1e-8);
276 ASSERT_TOL(helmholtz_w, T, rho, d, 2.52140783e2, 1e-8);
277
278 fprintf(stderr,"\nIAPWS95 TABLE 7 (SINGLE-PHASE) TESTS\n");
279 for(i=0; i<ntd; ++i){
280 double T = td[i].T;
281 double rho = td[i].rho;
282 double p = td[i].p * 1e6; /* Pa */
283 double cv = td[i].cv * 1e3; /* J/kgK */
284 double w = td[i].w; /* m/s */
285 double s = td[i].s * 1e3; /* J/kgK */
286 //fprintf(stderr,"T = %f, rho = %f, p = %f, w = %f, wcalc = %f\n",T,rho,p,w, helmholtz_w(T,rho,d));
287 ASSERT_TOL(helmholtz_s, T, rho, d, s, s*1e-8);
288 ASSERT_TOL(helmholtz_p, T, rho, d, p, p*1e-8);
289 ASSERT_TOL(helmholtz_cv, T, rho, d, cv, cv*1e-8);
290 ASSERT_TOL(helmholtz_w, T, rho, d, w, w*2e-5);
291 }
292
293
294 fprintf(stderr,"Tests completed OK (maximum error = %0.8f%%)\n",maxerr);
295 exit(0);
296 }
297
298 /* HIGHER-LEVEL TEST-DATA PROVIDED IN IAPWS95 */
299
300 const TestDataIAPWS95 td[] = {
301 {300, 0.9965560e3, 0.992418352e-1, 0.413018112e1, 0.150151914e4, 0.393062643}
302 ,{300, 0.1005308e4, 0.200022515e2, 0.406798347e1, 0.153492501e4, 0.387405401}
303 ,{300, 0.1188202e4, 0.700004704e3, 0.346135580e1, 0.244357992e4, 0.132609616}
304 ,{500, 0.4350000, 0.999679423e-1, 0.150817541e1, 0.548314253e3, 0.794488271e1}
305 ,{500, 0.4532000e1, 0.999938125, 0.166991025e1, 0.535739001e3, 0.682502725e1}
306 ,{500, 0.8380250e3, 0.100003858e2, 0.322106219e1, 0.127128441e4, 0.256690919e1}
307 ,{500, 0.1084564e4, 0.700000405e3, 0.307437693e1, 0.241200877e4, 0.203237509e1}
308 ,{647, 0.3580000e3, 0.220384756e2, 0.618315728e1, 0.252145078e3, 0.432092307e1}
309 ,{900, 0.2410000, 0.100062559, 0.175890657e1, 0.724027147e3, 0.916653194e1}
310 ,{900, 0.5261500e2, 0.200000690e2, 0.193510526e1, 0.698445674e3, 0.659070225e1}
311 ,{900, 0.8707690e3, 0.700000006e3, 0.266422350e1, 0.201933608e4, 0.417223802e1}
312 };
313
314 const unsigned ntd = sizeof(td)/sizeof(TestDataIAPWS95);
315
316 #endif
317

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