/[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 1985 - (show annotations) (download) (as text)
Tue Feb 3 01:15:17 2009 UTC (13 years, 4 months ago) by jpye
File MIME type: text/x-csrc
File size: 7536 byte(s)
Added critical term calculation in helm_resid.
Added gaussian and critical term values for water from IAPWS95.
Added some test data from IAPWS95 to water.c tests.
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
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 #ifdef TEST
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 int main(void){
185 double rho, T;
186 const HelmholtzData *d;
187
188 d = &helmholtz_data_water;
189 double maxerr = 0;
190
191 /* a simple macro to actually do the testing */
192 #define ASSERT_TOL(FN,PARAM1,PARAM2,PARAM3,VAL,TOL) {\
193 double cval; cval = FN(PARAM1,PARAM2,PARAM3);\
194 double err; err = cval - (double)(VAL);\
195 double relerrpc = (cval-(VAL))/(VAL)*100;\
196 if(fabs(relerrpc)>maxerr)maxerr=fabs(relerrpc);\
197 if(fabs(err)>fabs(TOL)){\
198 fprintf(stderr,"ERROR in line %d: value of '%s(%f,%f,%s)' = %0.8f,"\
199 " should be %f, error is %.10e (%.2f%%)!\n"\
200 , __LINE__, #FN,PARAM1,PARAM2,#PARAM3, cval, VAL,cval-(VAL)\
201 ,relerrpc\
202 );\
203 exit(1);\
204 }else{\
205 fprintf(stderr," OK, %s(%f,%f,%s) = %8.2e with %.8f%% err.\n"\
206 ,#FN,PARAM1,PARAM2,#PARAM3,VAL,relerrpc\
207 );\
208 }\
209 }
210
211
212 fprintf(stderr,"phi0 TESTS\n");
213 for(T = 300; T <= 900; T+= 100){
214 for(rho = 900; rho >= 0.9; rho*=0.5){
215 double delta = rho / d->rho_star;
216 double tau = d->T_star / T;
217 double p0 = phi0(delta,tau);
218
219 ASSERT_TOL(helm_ideal, tau, delta, d->ideal, p0, p0*1e-5);
220 }
221 }
222
223 fprintf(stderr,"IAPWS95 TABLE 6 TESTS\n");
224 T = 500.; /* K */
225 rho = 838.025; /* kg/m鲁 */
226 double tau = d->T_star / T;
227 double delta = rho / d->rho_star;
228
229 ASSERT_TOL(helm_ideal, tau, delta, d->ideal, 0.204797733E1, 1e-6);
230 ASSERT_TOL(helm_ideal_tau, tau, delta, d->ideal, 0.904611106E1, 1e-6);
231 ASSERT_TOL(HELM_IDEAL_DELTAU, tau, delta, d->ideal, 0., 1e-6);
232
233 double phitt = helm_ideal_tautau(tau, d->ideal);
234 double val = (-0.193249185E1);
235 double err = phitt - val;
236 if(fabs(err) > 1e-6){
237 fprintf(stderr,"ERROR in helm_ideal_tautau near line %d\n",__LINE__);
238 exit(1);
239 }else{
240 fprintf(stderr," OK, helm_ideal_tautau(%f) = %8.2e with %.6f%% err.\n"
241 ,tau,val,err/val*100.
242 );
243 }
244
245 /* FIXME still need to implement helm_ideal_del, helm_ideal_deldel, helm_ideal_deltau */
246
247 ASSERT_TOL(helm_resid, tau, delta, d, -0.342693206E1, 1e-8);
248 ASSERT_TOL(helm_resid_del, tau, delta, d, -0.364366650, 1e-8);
249 ASSERT_TOL(helm_resid_deldel, tau, delta, d, 0.856063701, 1e-8);
250 ASSERT_TOL(helm_resid_tau, tau, delta, d, -0.581403435E1, 1e-8);
251 ASSERT_TOL(helm_resid_tautau, tau, delta, d, -0.223440737E1, 1e-8);
252 ASSERT_TOL(helm_resid_deltau, tau, delta, d, -0.112176915e1, 1e-8);
253
254 fprintf(stderr,"Tests completed OK (maximum error = %0.2f%%)\n",maxerr);
255 exit(0);
256 }
257
258 #endif
259

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