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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 166 - (show annotations) (download) (as text)
Thu Jan 5 09:51:50 2006 UTC (16 years, 8 months ago) by johnpye
File MIME type: text/x-ascend
File size: 20935 byte(s)
Working on re-formulating the IAPWS-95 two-phase version with 'slack variables'
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 (* phase components *)
509 Sf IS_A iapws95_1phase;
510 Sg IS_A iapws95_1phase;
511 xf IS_A fraction;
512 xg IS_A fraction;
513 xfslack IS_A fraction;
514 xgslack IS_A fraction;
515
516 z01: xf + xg = 1.0;
517
518 (* complementarity... ??? *)
519 z02: xf * xfslack = 0.0; (* xfslack must be zero, else xf must be zero (ie no liquid) *)
520 z03: xg * xgslack = 0.0; (* xgslack must be zero, else xg must be zero (ie no gas) *)
521 (* so, if we have a non-zero slack variable, that phase will have disappeared and we don't need to worry! *)
522
523 (* saturation conditions *)
524 sat IS_A iapws_sat_curves;
525 sat.rhof, Sf.rho ARE_THE_SAME;
526 sat.rhog, Sg.rho ARE_THE_SAME;
527 Sf.T, Sg.T, sat.T ARE_THE_SAME;
528 rhof ALIASES Sf.rho;
529 rhog ALIASES Sg.rho;
530 T, Sf.T ARE_THE_SAME;
531 p, Sf.p ARE_THE_SAME;
532
533 x ALIASES xg;
534
535 z04: rhof * rhog = rho * ( xf * rhog + xg * rhof);
536 z05: u = xf * Sf.u + xg * Sg.u;
537 z06: h = xf * Sf.h + xg * Sg.h;
538 z07: s = xf * Sf.s + xg * Sg.s;
539 z08: cp = xf * Sf.cp + xg * Sg.cp;
540 z09: cv = xf * Sf.cv + xg * Sg.cv;
541 z10: w = xf * Sf.w + xg * Sg.w; (* check this *)
542 METHODS
543
544 METHOD default_self;
545 RUN ClearAll;
546 RUN specify;
547 RUN values;
548 END default_self;
549
550 METHOD specify;
551 T.fixed := TRUE;
552 x.fixed := TRUE;
553 END specify;
554
555 METHOD values;
556 T := 390 {K};
557 x := 0.5;
558 END values;
559
560 END iapws95_2phase;
561
562 (*
563 MODEL iapws95 REFINES thermo_state;
564 S1 IS_A iapws95_1phase;
565 S2 IS_A iapws95_2phase;
566 sat IS_A iapws_sat_curves;
567
568 T, sat.T, S1.T, S2.T ARE_THE_SAME;
569 rho, S1.rho, S2.rho ARE_THE_SAME;
570 h, S1.h, S2.h ARE_THE_SAME;
571 u, S1.u, S2.u ARE_THE_SAME;
572 p, S1.p, S2.p ARE_THE_SAME;
573 s, S1.s, S2.s ARE_THE_SAME;
574 cp, S1.cp, S2.cp ARE_THE_SAME;
575 cv, S1.cv, S2.cv ARE_THE_SAME;
576 w, S1.w, S2.w ARE_THE_SAME;
577
578 CONDITIONAL
579 sub_test: rho >= sat.rhof;
580 super_test: rho <= sat.rhog;
581 sat1_test: rho < sat.rhof;
582 sat2_test: rho > sat.rhog;
583 END CONDITIONAL;
584
585 issat1 IS_A boolean_var;
586 sat1_expr: issat1 == SATISFIED(sat1_test);
587 issat2 IS_A boolean_var;
588 sat2_expr: issat1 == SATISFIED(sat2_test);
589 issat IS_A boolean_var;
590 saturation_expr: issat == issat1 OR issat2;
591 issub IS_A boolean_var;
592 subcooled_expr: issub == SATISFIED(sub_test);
593 issuper IS_A boolean_var;
594 superheated_expr: issuper == SATISFIED(super_test);
595
596 saturation_when: WHEN (issat)
597 CASE TRUE:
598 USE S2;
599 CASE FALSE:
600 USE S1;
601 END WHEN;
602
603 METHODS
604 METHOD default_self;
605 RUN reset;
606 RUN values;
607 END default_self;
608
609 METHOD specify;
610 T.fixed := TRUE;
611 rho.fixed := TRUE;
612 END specify;
613
614 METHOD values;
615 (* fixed vars *)
616 T := 300 {K};
617 rho := 100 {kg/m^3};
618 (* initial values *)
619 issat := TRUE;
620 issuper := FALSE;
621 issub := FALSE;
622 p := 1 {atm};
623 h := 400 {kJ/kg};
624 END values;
625
626 END iapws95;

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