/[ascend]/trunk/models/johnpye/iapws95.a4c
ViewVC logotype

Contents of /trunk/models/johnpye/iapws95.a4c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 155 - (show annotations) (download) (as text)
Fri Dec 23 00:59:32 2005 UTC (16 years, 9 months ago) by johnpye
File MIME type: text/x-ascend
File size: 20420 byte(s)
Fixing up self-testing models
1 REQUIRE "system.a4l";
2 REQUIRE "atoms.a4l";
3 REQUIRE "johnpye/thermo_types.a4c";
4 REQUIRE "johnpye/iapws_sat_curves.a4c";
5
6 (*************************************************************************
7
8 The thermodynamic properties of water calculated with the
9 IAPWS95 equations. Variables and (example/possible) units:
10
11 T Temperature, K
12 rho Density, kg/m^3
13 p Pressure, MPa
14 u Specific internal energy, kJ/kg
15 h Specific enthalpy, kJ/kg
16 s Specific entropy, kJ/(kg*K)
17 cv Isochoric specific heat, kJ/(kg*K)
18 cp Isobaric specific heat, kJ/(kg*K)
19 w Speed of sound, m/s
20
21 References:
22
23 [1] The International Association for the Properties of
24 Water and Steam, "Release on the IAPWS Formulation 1995
25 for the Thermodynamic Properties of Ordinary Water
26 Substance for General and Scientific Use", dated
27 September 1996, Fredericia, Denmark. See the "Releases"
28 section of the website http://www.iapws.org/.
29
30 [2] NIST Chemistry Webbook:
31 http://webbook.nist.gov/chemistry/fluid/
32
33 ----------------------------------------------------------------------
34
35 freesteam-ascend - IAPWS-95 steam library for ASCEND
36 Copyright (C) John Pye 2005
37 derived from work by Don Peterson for freesteam, (C) 2004.
38
39 This program is free software; you can redistribute it
40 and/or modify it under the terms of the GNU General Public
41 License as published by the Free Software Foundation; either
42 version 2 of the License, or (at your option) any later
43 version.
44
45 This program is distributed in the hope that it will be
46 useful, but WITHOUT ANY WARRANTY; without even the implied
47 warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
48 PURPOSE. See the GNU General Public License for more
49 details.
50
51 You should have received a copy of the GNU General Public
52 License along with this program; if not, write to the Free
53 Software Foundation, Inc., 59 Temple Place, Suite 330,
54 Boston, MA 02111-1307 USA
55
56 *)
57
58 MODEL iapws95_1phase REFINES thermo_state;
59
60 delta IS_A positive_variable;
61 tau IS_A positive_variable;
62
63 (*-------------- CONSTANTS ---------------*)
64 rhoc IS_A mass_density_constant;
65 Tc IS_A temperature_constant;
66
67 rhoc "density of water at the critical point"
68 :== 322 {kg/m^3};
69
70 Tc "temperature of water at the critical point"
71 :== 647.096 {K};
72
73 R IS_A specific_gas_constant;
74 R "specific gas constant for water"
75 :== 0.46151805 {kJ/kg/K};
76
77 tau = Tc / T;
78 delta = rho / rhoc;
79
80
81 range_0 IS_A set OF integer_constant;
82 range_0 :== [1..8];
83
84 range_01 IS_A set OF integer_constant;
85 range_01 :== [4..8];
86
87 range_r1 IS_A set OF integer_constant;
88 range_r1 :== [1..7];
89
90 range_r2 IS_A set OF integer_constant;
91 range_r2 :== [8..51];
92
93 range_r3 IS_A set OF integer_constant;
94 range_r3 :== [52..54];
95
96 range_r4 IS_A set OF integer_constant;
97 range_r4 :== [55..56];
98
99 n0[range_0] IS_A real_constant;
100 n0[1] :== -8.32044648201;
101
102 n0[2] :== 6.6832105268;
103 n0[3] :== 3.00632;
104 n0[4] :== 0.012436;
105
106 n0[5] :== 0.97315;
107 n0[6] :== 1.27950;
108 n0[7] :== 0.96956;
109 n0[8] :== 0.24873;
110
111 gamma0[range_01] IS_A real_constant;
112 gamma0[4] :== 1.28728967;
113 gamma0[5] :== 3.53734222;
114 gamma0[6] :== 7.74073708;
115 gamma0[7] :== 9.24437796;
116 gamma0[8] :== 27.5075105;
117
118 n[1..56] IS_A real_constant;
119 n[1] :== 0.12533547935523e-1; n[2] :== 0.78957634722828e+1; n[3] :== -0.87803203303561e+1; n[4] :== 0.31802509345418e+0;
120 n[5] :== -0.26145533859358e+0; n[6] :== -0.78199751687981e-2; n[7] :== 0.88089493102134e-2;
121 n[8] :== -0.66856572307965e+0; n[9] :== 0.20433810950965e+0; n[10] :== -0.66212605039687e-4;
122 n[11] :== -0.19232721156002e+0; n[12] :== -0.25709043003438e+0; n[13] :== 0.16074868486251e+0; n[14] :== -0.40092828925807e-1;
123 n[15] :== 0.39343422603254e-6; n[16] :== -0.75941377088144e-5; n[17] :== 0.56250979351888e-3;
124 n[18] :== -0.15608652257135e-4; n[19] :== 0.11537996422951e-8; n[20] :== 0.36582165144204e-6;
125 n[21] :== -0.13251180074668e-11; n[22] :== -0.62639586912454e-9; n[23] :== -0.10793600908932e+0; n[24] :== 0.17611491008752e-1;
126 n[25] :== 0.22132295167546e+0; n[26] :== -0.40247669763528e+0; n[27] :== 0.58083399985759e+0;
127 n[28] :== 0.49969146990806e-2; n[29] :== -0.31358700712549e-1; n[30] :== -0.74315929710341e+0;
128 n[31] :== 0.47807329915480e+0; n[32] :== 0.20527940895948e-1; n[33] :== -0.13636435110343e+0; n[34] :== 0.14180634400617e-1;
129 n[35] :== 0.83326504880713e-2; n[36] :== -0.29052336009585e-1; n[37] :== 0.38615085574206e-1;
130 n[38] :== -0.20393486513704e-1; n[39] :== -0.16554050063734e-2; n[40] :== 0.19955571979541e-2;
131 n[41] :== 0.15870308324157e-3; n[42] :== -0.16388568342530e-4; n[43] :== 0.43613615723811e-1; n[44] :== 0.34994005463765e-1;
132 n[45] :== -0.76788197844621e-1; n[46] :== 0.22446277332006e-1; n[47] :== -0.62689710414685e-4;
133 n[48] :== -0.55711118565645e-9; n[49] :== -0.19905718354408e+0; n[50] :== 0.31777497330738e+0;
134 n[51] :== -0.11841182425981e+0; n[52] :== -0.31306260323435e+2; n[53] :== 0.31546140237781e+2; n[54] :== -0.25213154341695e+4;
135 n[55] :== -0.14874640856724e+0; n[56] :== 0.31806110878444e+0;
136
137 c[1..51] IS_A integer_constant;
138 c[1] :== 0; c[2] :== 0; c[3] :== 0; c[4] :== 0; c[5] :== 0;
139 c[6] :== 0; c[7] :== 0; c[8] :== 1; c[9] :== 1; c[10] :== 1;
140 c[11] :== 1; c[12] :== 1; c[13] :== 1; c[14] :== 1; c[15] :== 1;
141 c[16] :== 1; c[17] :== 1; c[18] :== 1; c[19] :== 1; c[20] :== 1;
142 c[21] :== 1; c[22] :== 1; c[23] :== 2; c[24] :== 2; c[25] :== 2;
143 c[26] :== 2; c[27] :== 2; c[28] :== 2; c[29] :== 2; c[30] :== 2;
144 c[31] :== 2; c[32] :== 2; c[33] :== 2; c[34] :== 2; c[35] :== 2;
145 c[36] :== 2; c[37] :== 2; c[38] :== 2; c[39] :== 2; c[40] :== 2;
146 c[41] :== 2; c[42] :== 2; c[43] :== 3; c[44] :== 3; c[45] :== 3;
147 c[46] :== 3; c[47] :== 4; c[48] :== 6; c[49] :== 6; c[50] :== 6;
148 c[51] :== 6;
149
150 d[1..54] IS_A integer_constant;
151 d[1] :== 1; d[2] :== 1; d[3] :== 1; d[4] :== 2; d[5] :== 2;
152 d[6] :== 3; d[7] :== 4; d[8] :== 1; d[9] :== 1; d[10] :== 1;
153 d[11] :== 2; d[12] :== 2; d[13] :== 3; d[14] :== 4; d[15] :== 4;
154 d[16] :== 5; d[17] :== 7; d[18] :== 9; d[19] :== 10; d[20] :== 11;
155 d[21] :== 13; d[22] :== 15; d[23] :== 1; d[24] :== 2; d[25] :== 2;
156 d[26] :== 2; d[27] :== 3; d[28] :== 4; d[29] :== 4; d[30] :== 4;
157 d[31] :== 5; d[32] :== 6; d[33] :== 6; d[34] :== 7; d[35] :== 9;
158 d[36] :== 9; d[37] :== 9; d[38] :== 9; d[39] :== 9; d[40] :== 10;
159 d[41] :== 10; d[42] :== 12; d[43] :== 3; d[44] :== 4; d[45] :== 4;
160 d[46] :== 5; d[47] :== 14; d[48] :== 3; d[49] :== 6; d[50] :== 6;
161 d[51] :== 6; d[52] :== 3; d[53] :== 3; d[54] :== 3;
162
163 t[1..54] IS_A real_constant;
164 t[1] :== -0.5; t[2] :== 0.875; t[3] :== 1; t[4] :== 0.5; t[5] :== 0.75;
165 t[6] :== 0.375; t[7] :== 1; t[8] :== 4; t[9] :== 6; t[10] :== 12;
166 t[11] :== 1; t[12] :== 5; t[13] :== 4; t[14] :== 2; t[15] :== 13;
167 t[16] :== 9; t[17] :== 3; t[18] :== 4; t[19] :== 11; t[20] :== 4;
168 t[21] :== 13; t[22] :== 1; t[23] :== 7; t[24] :== 1; t[25] :== 9;
169 t[26] :== 10; t[27] :== 10; t[28] :== 3; t[29] :== 7; t[30] :== 10;
170 t[31] :== 10; t[32] :== 6; t[33] :== 10; t[34] :== 10; t[35] :== 1;
171 t[36] :== 2; t[37] :== 3; t[38] :== 4; t[39] :== 8; t[40] :== 6;
172 t[41] :== 9; t[42] :== 8; t[43] :== 16; t[44] :== 22; t[45] :== 23;
173 t[46] :== 23; t[47] :== 10; t[48] :== 50; t[49] :== 44; t[50] :== 46;
174 t[51] :== 50; t[52] :== 0; t[53] :== 1; t[54] :== 4;
175
176 (* Correlation parameters *)
177
178 (* TODO convert from C to ASCEND arrays? Note trickiness with 0- and 1-based array indices. *)
179
180 a[range_r4] IS_A real_constant;
181 a[55]:==3.5;
182 a[56]:==3.5;
183
184 b[range_r4] IS_A real_constant;
185 b[55]:== 0.85;
186 b[56]:== 0.95;
187
188 A[range_r4] IS_A real_constant;
189 A[55]:==0.32;
190 A[56]:==0.32;
191
192 B[range_r4] IS_A real_constant;
193 B[55]:==0.2;
194 B[56]:==0.2;
195
196 C[range_r4] IS_A real_constant;
197 C[55]:==28;
198 C[56]:==32;
199
200 D[range_r4] IS_A real_constant;
201 D[55]:==700;
202 D[56]:==800;
203
204 beta_r4[range_r4] IS_A real_constant;
205 beta_r4[55]:==0.3;
206 beta_r4[56]:==0.3;
207
208 alpha[range_r3] IS_A integer_constant;
209 alpha[52]:==20;
210 alpha[53]:==20;
211 alpha[54]:==20;
212
213 beta[range_r3] IS_A real_constant;
214 beta[52]:==150;
215 beta[53]:==150;
216 beta[54]:==250;
217
218 gamma[range_r3] IS_A real_constant;
219 gamma[52]:==1.21;
220 gamma[53]:==1.21;
221 gamma[54]:==1.25;
222
223 (*--------------- DIMENSIONLESS PARTIAL DERIVATIVES ---------------- *)
224
225 (*------------ IDEAL PARTS ------------*)
226
227 phi0 IS_A factor;
228 phi0_expr: phi0 =
229 SUM[ n0[i]*ln(1-exp(-tau*gamma0[i])) | i IN [range_01] ]
230 + ln(delta) + n0[1] + n0[2]*tau + n0[3]*ln(tau);
231
232 phi0delta IS_A factor;
233 phi0delta_expr: phi0delta = 1.0/delta;
234
235 phi0deltadelta IS_A factor;
236 phi0deltadelta_expr: phi0deltadelta =
237 -1.0/(delta*delta);
238
239
240 phi0tau IS_A factor;
241 phi0tau_expr: phi0tau =
242 n0[2] + n0[3]/tau
243 + SUM[ n0[i]*gamma0[i]*(1/(1-exp(-tau*gamma0[i])) - 1) | i IN [range_01] ];
244
245 phi0deltatau IS_A real_constant;
246 phi0deltatau :== 0.0;
247
248 phi0tautau IS_A factor;
249 phi0tautau_expr: phi0tautau
250 = -n0[3] / tau^2
251 - SUM [ n0[i] * gamma0[i]^2 * exp(-gamma0[i] * tau) / ( 1 - exp(-gamma0[i] * tau) )^2 | i IN range_01 ];
252
253 (*----------- 'REAL' PARTS -- CLOSE YOUR EYES -----------*)
254
255 d1 IS_A factor;
256 d1_expr: d1 = delta - 1;
257
258 t1 IS_A factor;
259 t1_expr: t1 = tau -1;
260
261 r3_b1[range_r3] IS_A factor;
262 FOR i IN range_r3 CREATE
263 r3_b1[i] = -alpha[i]*d1^2 - beta[i]* (tau - gamma[i])^2;
264 END FOR;
265
266 PSI[range_r4] IS_A factor;
267 theta[range_r4] IS_A factor;
268 DELTA[range_r4] IS_A factor;
269 dDELTA_ddelta[range_r4] IS_A factor;
270 dPSI_ddelta[range_r4] IS_A factor;
271 dDELTAbi_ddelta[range_r4] IS_A factor;
272 dDELTAbi_dtau[range_r4] IS_A factor;
273 dPSI_dtau[range_r4] IS_A factor;
274 d2DELTA_ddelta2[range_r4] IS_A factor;
275 d2DELTAbi_ddelta2[range_r4] IS_A factor;
276 d2PSI_ddelta2[range_r4] IS_A factor;
277 d2DELTAbi_dtau2[range_r4] IS_A factor;
278 d2PSI_dtau2[range_r4] IS_A factor;
279 d2PSI_ddeltadtau[range_r4] IS_A factor;
280 d2DELTAbi_ddeltadtau[range_r4] IS_A factor;
281
282 FOR i IN range_r4 CREATE
283 PSI[i] = exp(-C[i]*d1^2 - D[i]*t1^2);
284 (*theta_expr:*) theta[i] = 1 - tau + A[i] * (d1^2)^( 1/(2 * beta_r4[i]) );
285 (*DELTA_expr:*) DELTA[i] = theta[i]^2 + B[i] * (d1^2)^a[i];
286
287 (*dDELTA_ddelta_expr:*) dDELTA_ddelta[i] = d1*(A[i]*theta[i]*2/beta_r4[i]*
288 (d1^2)^(1/(2*beta_r4[i]) - 1) + 2*B[i]*a[i]*
289 (d1^2)^(a[i] - 1) );
290
291 (*dPSI_ddelta_expr:*) dPSI_ddelta[i] = -2*C[i]*d1*PSI[i];
292
293 (*dDELTAbi_ddelta_expr:*) dDELTAbi_ddelta[i] = b[i] * DELTA[i]^(b[i]-1) * dDELTA_ddelta[i];
294
295 (*dDELTAbi_dtau_expr:*) dDELTAbi_dtau[i] = -2 * theta[i] * b[i] * DELTA[i]^(b[i]-1);
296 (*dPSI_dtau_expr:*) dPSI_dtau[i] = -2 * D[i] * t1 * PSI[i];
297
298 (*d2PSI_ddelta2_expr:*) d2PSI_ddelta2[i] = (2 * C[i] * d1^2 - 1) * 2 * C[i] * PSI[i];
299
300 (*d2DELTA_ddelta2_expr:*) d2DELTA_ddelta2[i] = 1/d1*dDELTA_ddelta[i] + d1^2*(4*B[i]*a[i]*
301 (a[i]-1)*(d1^2)^(a[i]-2) + 2*A[i]^2*
302 (1/(beta_r4[i]^2))*((d1^2) ^ (1/(2*beta_r4[i])-1))^2) +
303 A[i]*theta[i]*4/beta_r4[i]*(1/(2*beta_r4[i]) - 1)*
304 (d1^2)^(1/(2*beta_r4[i]) - 2);
305
306 (*d2DELTAbi_ddelta2_expr:*) d2DELTAbi_ddelta2[i] =b[i] * (
307 DELTA[i]^(b[i]-1) *d2DELTA_ddelta2[i]
308 + (b[i]-1) * DELTA[i]^(b[i]-2) * dDELTA_ddelta[i]^2 );
309
310 (*d2DELTAbi_dtau2_expr:*) d2DELTAbi_dtau2[i] = 2 * b[i] * DELTA[i]^(b[i]-1) +
311 4 * theta[i]^2 * b[i]*(b[i]-1) * DELTA[i]^(b[i]-2);
312
313 (*d2PSI_dtau2_expr:*) d2PSI_dtau2[i] = (2 * D[i] * t1^2 - 1) * 2 * D[i] * PSI[i];
314
315 (*d2PSI_ddeltadtau_expr:*) d2PSI_ddeltadtau[i] = 4 * C[i] * D[i] * d1 * t1 * PSI[i];
316
317 (*d2DELTAbi_ddeltadtau_expr:*) d2DELTAbi_ddeltadtau[i] =
318 - A[i] * b[i]* (2/beta_r4[i]) * DELTA[i]^(b[i]-1) * d1 * (d1^2)^(1/(2*beta_r4[i])-1)
319 - 2 * theta[i] * b[i] *(b[i]-1) * DELTA[i]^(b[i]-2) * dDELTA_ddelta[i];
320
321 END FOR;
322
323 phir_r2[range_r2] IS_A factor;
324 FOR i IN range_r2 CREATE
325 phir_r2[i] = n[i] * delta^d[i] * tau^t[i] *
326 exp(-delta^c[i]);
327 END FOR;
328
329 phir IS_A factor;
330 phir_expr: phir
331 =
332 SUM[ n[i] * delta^d[i] * tau^t[i] | i IN [range_r1] ]
333 + SUM[ phir_r2[i] | i IN [range_r2] ]
334 + SUM[ n[i] * delta^d[i] * tau^t[i] * exp(
335 -alpha[i]*d1^2 - beta[i]*
336 (tau - gamma[i])*(tau - gamma[i])
337 ) | i IN [range_r3] ]
338 + SUM[ n[i] * DELTA[i]^b[i] * delta * PSI[i] | i IN [range_r4] ];
339
340
341 phirdelta_r2[range_r2] IS_A factor;
342 FOR i IN range_r2 CREATE
343 phirdelta_r2[i] = n[i] * exp(-delta^c[i]) * (delta^(d[i]-1) *
344 tau^t[i] * (d[i] - c[i] * delta^c[i]) );
345 END FOR;
346
347 phirdelta IS_A factor;
348 phirdelta_expr: phirdelta =
349 SUM[ n[i] * d[i] * delta^(d[i] - 1) * tau^t[i] | i IN range_r1 ]
350 + SUM[ phirdelta_r2[i] | i IN range_r2 ]
351 + SUM[ n[i]*delta^d[i] * tau^t[i] * exp( r3_b1[i] ) *
352 (d[i]/delta - 2*alpha[i]*d1) | i IN range_r3 ]
353 + SUM[ n[i] * (
354 DELTA[i]^b[i] * (PSI[i] + delta * dPSI_ddelta[i] )
355 + dDELTAbi_ddelta[i] * delta * PSI[i] ) | i IN range_r4 ];
356
357 phirtau_r2[range_r2] IS_A factor;
358 FOR i IN range_r2 CREATE
359 phirtau_r2[i] = n[i]*t[i]*delta^d[i]*tau^(t[i]-1) * exp(-delta^c[i]);
360 END FOR;
361
362 phirtau IS_A factor;
363 phirtau_expr: phirtau =
364 SUM[ n[i] * t[i] * delta^d[i] * tau^(t[i]-1) | i IN range_r1 ]
365 + SUM[ phirtau_r2[i] | i IN range_r2 ]
366 + SUM[ n[i] * delta^d[i] * tau^t[i] * exp( r3_b1[i] )*
367 (t[i]/tau - 2*beta[i]*(tau - gamma[i])) | i IN range_r3 ]
368 + SUM[ n[i]*delta*(dDELTAbi_dtau[i] * PSI[i] +
369 DELTA[i]^b[i]*dPSI_dtau[i]) | i IN range_r4 ];
370
371 phirdeltadelta_r2[range_r2] IS_A factor;
372 FOR i IN range_r2 CREATE
373 phirdeltadelta_r2[i] =
374 n[i] * exp(-delta^c[i]) * (
375 delta^(d[i]-2) *
376 tau^t[i] * (
377 ( d[i] - c[i]*delta^c[i] )*(d[i]- 1 - c[i] * delta^c[i] )
378 - c[i]^2 * delta^c[i]
379 )
380 );
381 END FOR;
382
383 phirdeltadelta IS_A factor;
384 phirdeltadelta_expr: phirdeltadelta =
385 SUM[
386 n[i] * d[i] * (d[i] - 1) * delta^(d[i]-2) * tau^t[i] | i IN range_r1 ]
387 + SUM [phirdeltadelta_r2[i] | i IN range_r2 ]
388 + SUM [
389 n[i]* tau^t[i] * exp( r3_b1[i] ) * (
390 - 2 * alpha[i] * delta^d[i]
391 + 4 * alpha[i]*alpha[i] * delta^d[i] * d1^2 (* d1 = d-1 = d-eps *)
392 - 4 * d[i] * alpha[i] * delta^(d[i]-1) * d1
393 + d[i] * (d[i]-1) * delta^(d[i]-1)
394 ) | i IN range_r3 ]
395 + SUM [n[i]*(
396 DELTA[i]^b[i] * (2*dPSI_ddelta[i] + delta*d2PSI_ddelta2[i] )
397 + 2*dDELTAbi_ddelta[i] * (PSI[i] + delta*dPSI_ddelta[i])
398 + d2DELTAbi_ddelta2[i] * delta * PSI[i] ) | i IN range_r4 ];
399
400
401
402 phirtautau_r2[range_r2] IS_A factor;
403 FOR i IN range_r2 CREATE
404 phirtautau_r2[i] = n[i]*t[i]*(t[i]-1)* delta^d[i] *
405 tau^(t[i]-2) * exp(- delta^c[i] );
406 END FOR;
407
408 phirtautau IS_A factor;
409 phirtautau_expr: phirtautau =
410 SUM[ n[i] * t[i] * (t[i]-1) * delta^d[i] *
411 tau^(t[i]-2) | i IN range_r1 ]
412 + SUM[ phirtautau_r2[i] | i IN range_r2 ]
413 + SUM[ n[i] * delta^d[i] * tau^t[i] * exp( r3_b1[i] ) *
414 ( (t[i]/tau - 2*beta[i]* (tau - gamma[i]) )^2 -
415 t[i]/tau^2 - 2*beta[i] ) | i IN range_r3 ]
416 + SUM[ n[i] * delta * (d2DELTAbi_dtau2[i] * PSI[i] +
417 2 * dDELTAbi_dtau[i] * dPSI_dtau[i]
418 + DELTA[i]^b[i] * d2PSI_dtau2[i] ) | i IN range_r4 ];
419
420
421 phirdeltatau_r2[range_r2] IS_A factor;
422 FOR i IN range_r2 CREATE
423 phirdeltatau_r2[i] =
424 n[i]*t[i] * delta^(d[i]-1) * tau^(t[i]-1)
425 * (d[i] - c[i] * delta^c[i]) * exp(-delta^c[i]);
426 END FOR;
427
428 phirdeltatau IS_A factor;
429 phirdeltatau_expr: phirdeltatau =
430 SUM[ n[i]*d[i]*t[i] * delta^(d[i]-1) * tau^(t[i]-1) | i IN range_r1 ]
431
432 + SUM[ phirdeltatau_r2[i] | i IN range_r2 ]
433 + SUM[ n[i] * delta^d[i] * tau^t[i] * exp( r3_b1[i] ) *
434 (d[i]/delta - 2*alpha[i]*d1)
435 * (t[i]/tau - 2*beta[i]* (tau - gamma[i]) ) | i IN range_r3 ]
436 + SUM[ n[i]*( DELTA[i]^b[i] * (dPSI_dtau[i] + delta *
437 d2PSI_ddeltadtau[i] ) + delta * dDELTAbi_ddelta[i] * dPSI_dtau[i] +
438 dDELTAbi_dtau[i] * (PSI[i] + delta * dPSI_ddelta[i] ) +
439 d2DELTAbi_ddeltadtau[i] * delta * PSI[i] ) | i IN range_r4 ];
440
441 (*--------- THERMO PROPERTY RELATIONS ----------- *)
442
443 pressure: p
444 = rho * R * T * (1 + delta*phirdelta);
445
446 internal_energy: u
447 = R * T * tau * (phi0tau + phirtau);
448
449 enthalpy: h
450 = R * T * (1 + tau*(phi0tau + phirtau) + delta*phirdelta);
451
452 entropy: s
453 = R * (tau*(phi0tau + phirtau) - phi0 - phir);
454
455 c_isochoric: cv
456 = - R * tau^2 * (phi0tautau + phirtautau);
457
458 c_isobaric: cp
459 = - R * (
460 tau^2 * (phi0tautau + phirtautau)
461 + ( ( 1 + delta*phirdelta - delta*tau*phirdeltatau )^2 )
462 / ( 1 + 2*delta*phirdelta + delta^2 * phirdeltadelta )
463 );
464
465 spd_sound: w^2
466 = R * T * ( 1 + 2*delta*phirdelta + delta^2 * phirdeltadelta -
467 ( 1 + delta*phirdelta - delta*tau*phirdeltatau )^2
468 / ( tau^2 * (phi0tautau + phirtautau) )
469 );
470
471 METHODS
472 METHOD default_self;
473 RUN ClearAll;
474 RUN specify;
475 RUN values;
476 END default_self;
477
478 METHOD specify;
479 T.fixed := TRUE;
480 rho.fixed := TRUE;
481 END specify;
482
483 METHOD values;
484 (* these are the test values from page 14 of the IAPWS-95 release *)
485 T := 500 {K};
486 rho := 838.025 {kg/m^3};
487 END values;
488
489 METHOD self_test;
490 (* test the result values from page 14 of the IAPWS-95 release *)
491 ASSERT abs(phi0 - 0.204797734E1) < 1e-8;
492 ASSERT abs(phir - (-0.342693206E1)) < 1e-8;
493 ASSERT abs(phi0delta - (0.384236747)) < 1e-9;
494 ASSERT abs(phirdelta - (-0.364366650)) < 1e-9;
495 ASSERT abs(phi0deltadelta - (-0.147637878)) < 1e-9;
496 ASSERT abs(phirdeltadelta - (0.856063701)) < 1e-9;
497 ASSERT abs(phi0tau - (0.904611106e1)) < 1e-8;
498 ASSERT abs(phirtau - (-0.581403435e1)) < 1e-8;
499 ASSERT abs(phi0tautau - (-0.193249185e1)) < 1e-8;
500 ASSERT abs(phirtautau - (-0.223440737e1)) < 1e-8;
501 ASSERT abs(phi0deltatau - 0) < 1e-20;
502 ASSERT abs(phirdeltatau - (-0.112176915e1)) < 1e-8;
503 END self_test;
504
505 END iapws95_1phase;
506
507 MODEL iapws95_2phase REFINES thermo_state;
508 Sf IS_A iapws95_1phase;
509 Sg IS_A iapws95_1phase;
510 sat IS_A iapws_sat_curves;
511 sat.rhof, Sf.rho ARE_THE_SAME;
512 sat.rhog, Sg.rho ARE_THE_SAME;
513 Sf.T, Sg.T, sat.T ARE_THE_SAME;
514 rhof ALIASES Sf.rho;
515 rhog ALIASES Sg.rho;
516 T, Sf.T ARE_THE_SAME;
517 p, Sf.p ARE_THE_SAME;
518
519 x IS_A fraction;
520
521 rhog*rhof = rho*rhog + rho*(rhof - rhog) *x;
522
523 u = Sf.u + (Sg.u -Sf.u) *x;
524
525 h = Sf.h + (Sg.h -Sf.h) *x;
526 s = Sf.s + (Sg.s -Sf.s) *x;
527 cp = Sf.cp + (Sg.cp -Sf.cp) *x;
528 cv = Sf.cv + (Sg.cv -Sf.cv) *x;
529 w = Sf.w + (Sg.w - Sf.w) *x; (* check this *)
530 METHODS
531
532 METHOD default_self;
533 RUN ClearAll;
534 RUN specify;
535 RUN values;
536 END default_self;
537
538 METHOD specify;
539 T.fixed := TRUE;
540 x.fixed := TRUE;
541 END specify;
542
543 METHOD values;
544 T := 390 {K};
545 x := 0.5;
546 END values;
547
548 END iapws95_2phase;
549
550 MODEL iapws95 REFINES thermo_state;
551 S1 IS_A iapws95_1phase;
552 S2 IS_A iapws95_2phase;
553 sat IS_A iapws_sat_curves;
554
555 T, sat.T, S1.T, S2.T ARE_THE_SAME;
556 rho, S1.rho, S2.rho ARE_THE_SAME;
557 h, S1.h, S2.h ARE_THE_SAME;
558 u, S1.u, S2.u ARE_THE_SAME;
559 p, S1.p, S2.p ARE_THE_SAME;
560 s, S1.s, S2.s ARE_THE_SAME;
561 cp, S1.cp, S2.cp ARE_THE_SAME;
562 cv, S1.cv, S2.cv ARE_THE_SAME;
563 w, S1.w, S2.w ARE_THE_SAME;
564
565 CONDITIONAL
566 sub_test: rho >= sat.rhof;
567 super_test: rho <= sat.rhog;
568 sat1_test: rho < sat.rhof;
569 sat2_test: rho > sat.rhog;
570 END CONDITIONAL;
571
572 issat1 IS_A boolean_var;
573 sat1_expr: issat1 == SATISFIED(sat1_test);
574 issat2 IS_A boolean_var;
575 sat2_expr: issat1 == SATISFIED(sat2_test);
576 issat IS_A boolean_var;
577 saturation_expr: issat == issat1 OR issat2;
578 issub IS_A boolean_var;
579 subcooled_expr: issub == SATISFIED(sub_test);
580 issuper IS_A boolean_var;
581 superheated_expr: issuper == SATISFIED(super_test);
582
583 saturation_when: WHEN (issat)
584 CASE TRUE:
585 USE S2;
586 CASE FALSE:
587 USE S1;
588 END WHEN;
589
590 METHODS
591 METHOD default_self;
592 RUN reset;
593 RUN values;
594 END default_self;
595
596 METHOD specify;
597 T.fixed := TRUE;
598 rho.fixed := TRUE;
599 END specify;
600
601 METHOD values;
602 (* fixed vars *)
603 T := 300 {K};
604 rho := 100 {kg/m^3};
605 (* initial values *)
606 issat := TRUE;
607 issuper := FALSE;
608 issub := FALSE;
609 p := 1 {atm};
610 h := 400 {kJ/kg};
611 END values;
612
613 END iapws95;

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