/[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 1989 - (show annotations) (download) (as text)
Tue Feb 3 07:29:29 2009 UTC (13 years, 4 months ago) by jpye
File MIME type: text/x-csrc
File size: 9687 byte(s)
Added helmholtz_w and helmholtz_cp.
Added test for helmholtz_w with water.
This exposes the fact that helm_resid_deldel and helm_resid_deltau need critical terms added.
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 fprintf(stderr,"phi0 TESTS\n");
216 for(T = 300; T <= 900; T+= 100){
217 for(rho = 900; rho >= 0.9; rho*=0.5){
218 double delta = rho / d->rho_star;
219 double tau = d->T_star / T;
220 double p0 = phi0(delta,tau);
221
222 ASSERT_TOL(helm_ideal, tau, delta, d->ideal, p0, p0*1e-5);
223 }
224 }
225
226 /* LOW-LEVEL TEST DATA PROVIDED IN IAPWS95 */
227
228 fprintf(stderr,"\nIAPWS95 TABLE 6 TESTS\n");
229 T = 500.; /* K */
230 rho = 838.025; /* kg/m鲁 */
231 double tau = d->T_star / T;
232 double delta = rho / d->rho_star;
233
234 ASSERT_TOL(helm_ideal, tau, delta, d->ideal, 0.204797733E1, 1e-8);
235 ASSERT_TOL(helm_ideal_tau, tau, delta, d->ideal, 0.904611106E1, 1e-8);
236 ASSERT_TOL(HELM_IDEAL_DELTAU, tau, delta, d->ideal, 0., 1e-8);
237
238 double phitt = helm_ideal_tautau(tau, d->ideal);
239 double val = (-0.193249185E1);
240 double err = phitt - val;
241 if(fabs(err) > 1e-8){
242 fprintf(stderr,"ERROR in helm_ideal_tautau near line %d\n",__LINE__);
243 exit(1);
244 }else{
245 fprintf(stderr," OK, helm_ideal_tautau(%f) = %8.2e with %.6f%% err.\n"
246 ,tau,val,err/val*100.
247 );
248 }
249
250 /* FIXME still need to implement helm_ideal_del, helm_ideal_deldel, helm_ideal_deltau */
251
252 ASSERT_TOL(helm_resid, tau, delta, d, -0.342693206E1, 1e-8);
253 ASSERT_TOL(helm_resid_del, tau, delta, d, -0.364366650, 1e-8);
254 ASSERT_TOL(helm_resid_deldel, tau, delta, d, 0.856063701, 1e-8);
255 ASSERT_TOL(helm_resid_tau, tau, delta, d, -0.581403435E1, 1e-8);
256 ASSERT_TOL(helm_resid_tautau, tau, delta, d, -0.223440737E1, 1e-8);
257 ASSERT_TOL(helm_resid_deltau, tau, delta, d, -0.112176915e1, 1e-8);
258
259 fprintf(stderr,"\nADDITIONAL LOW-LEVEL TESTS NEAR CRITICAL POINT\n");
260
261 T = 647.; /* K */
262 rho = 358.; /* kg/m鲁 */
263 tau = d->T_star / T;
264 delta = rho / d->rho_star;
265
266 /* this test value calculated from pressure using REFPROP 8 */
267 ASSERT_TOL(helmholtz_a, T, rho, d, -8.286875181e5, 1e-4);
268 ASSERT_TOL(helm_resid_del, tau, delta, d, -7.14012024e-1, 1e-8);
269
270 fprintf(stderr,"\nIAPWS95 SINGLE-PHASE TESTS\n");
271 for(i=0; i<ntd; ++i){
272 double T = td[i].T;
273 double rho = td[i].rho;
274 double p = td[i].p * 1e6; /* Pa */
275 double cv = td[i].cv * 1e3; /* J/kgK */
276 double w = td[i].w; /* m/s */
277 double s = td[i].s * 1e3; /* J/kgK */
278 //fprintf(stderr,"T = %f, rho = %f, p = %f, w = %f, wcalc = %f\n",T,rho,p,w, helmholtz_w(T,rho,d));
279 ASSERT_TOL(helmholtz_s, T, rho, d, s, s*1e-8);
280 ASSERT_TOL(helmholtz_p, T, rho, d, p, p*1e-8);
281 ASSERT_TOL(helmholtz_cv, T, rho, d, cv, cv*1e-8);
282 ASSERT_TOL(helmholtz_w, T, rho, d, w, w*1e-8);
283 }
284
285 fprintf(stderr,"Tests completed OK (maximum error = %0.2f%%)\n",maxerr);
286 exit(0);
287 }
288
289 /* HIGHER-LEVEL TEST-DATA PROVIDED IN IAPWS95 */
290
291 const TestDataIAPWS95 td[] = {
292 {300, 0.9965560e3, 0.992418352e-1, 0.413018112e1, 0.150151914e4, 0.393062643}
293 ,{300, 0.1005308e4, 0.200022515e2, 0.406798347e1, 0.153492501e4, 0.387405401}
294 ,{300, 0.1188202e4, 0.700004704e3, 0.346135580e1, 0.244357992e4, 0.132609616}
295 ,{500, 0.4350000, 0.999679423e-1, 0.150817541e1, 0.548314253e3, 0.794488271e1}
296 ,{500, 0.4532000e1, 0.999938125, 0.166991025e1, 0.535739001e3, 0.682502725e1}
297 ,{500, 0.8380250e3, 0.100003858e2, 0.322106219e1, 0.127128441e4, 0.256690919e1}
298 ,{500, 0.1084564e4, 0.700000405e3, 0.307437693e1, 0.241200877e4, 0.203237509e1}
299 ,{647, 0.3580000e3, 0.220384756e2, 0.618315728e1, 0.252145078e3, 0.432092307e1}
300 ,{900, 0.2410000, 0.100062559, 0.175890657e1, 0.724027147e3, 0.916653194e1}
301 ,{900, 0.5261500e2, 0.200000690e2, 0.193510526e1, 0.698445674e3, 0.659070225e1}
302 ,{900, 0.8707690e3, 0.700000006e3, 0.266422350e1, 0.201933608e4, 0.417223802e1}
303 };
304
305 const unsigned ntd = sizeof(td)/sizeof(TestDataIAPWS95);
306
307 #endif
308

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