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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 327 - (show annotations) (download) (as text)
Fri Feb 24 02:04:49 2006 UTC (14 years, 8 months ago) by johnpye
File MIME type: text/x-ascend
File size: 21395 byte(s)
Adding some labels to equations in iapws95.a4c file
Adding some scaling commands, but still to no avail.
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 http://freesteam.sourceforge.net/
39
40 This program is free software; you can redistribute it
41 and/or modify it under the terms of the GNU General Public
42 License as published by the Free Software Foundation; either
43 version 2 of the License, or (at your option) any later
44 version.
45
46 This program is distributed in the hope that it will be
47 useful, but WITHOUT ANY WARRANTY; without even the implied
48 warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
49 PURPOSE. See the GNU General Public License for more
50 details.
51
52 You should have received a copy of the GNU General Public
53 License along with this program; if not, write to the Free
54 Software Foundation, Inc., 59 Temple Place, Suite 330,
55 Boston, MA 02111-1307 USA
56
57 *)
58
59 MODEL iapws95_1phase REFINES thermo_state;
60
61 delta IS_A positive_variable;
62 tau IS_A positive_variable;
63
64 (*-------------- CONSTANTS ---------------*)
65 rhoc IS_A mass_density_constant;
66 Tc IS_A temperature_constant;
67
68 rhoc "density of water at the critical point"
69 :== 322 {kg/m^3};
70
71 Tc "temperature of water at the critical point"
72 :== 647.096 {K};
73
74 R IS_A specific_gas_constant;
75 R "specific gas constant for water"
76 :== 0.46151805 {kJ/kg/K};
77
78 z_tau: tau = Tc / T;
79 z_delta: delta = rho / rhoc;
80
81
82 range_0 IS_A set OF integer_constant;
83 range_0 :== [1..8];
84
85 range_01 IS_A set OF integer_constant;
86 range_01 :== [4..8];
87
88 range_r1 IS_A set OF integer_constant;
89 range_r1 :== [1..7];
90
91 range_r2 IS_A set OF integer_constant;
92 range_r2 :== [8..51];
93
94 range_r3 IS_A set OF integer_constant;
95 range_r3 :== [52..54];
96
97 range_r4 IS_A set OF integer_constant;
98 range_r4 :== [55..56];
99
100 n0[range_0] IS_A real_constant;
101 n0[1] :== -8.32044648201;
102
103 n0[2] :== 6.6832105268;
104 n0[3] :== 3.00632;
105 n0[4] :== 0.012436;
106
107 n0[5] :== 0.97315;
108 n0[6] :== 1.27950;
109 n0[7] :== 0.96956;
110 n0[8] :== 0.24873;
111
112 gamma0[range_01] IS_A real_constant;
113 gamma0[4] :== 1.28728967;
114 gamma0[5] :== 3.53734222;
115 gamma0[6] :== 7.74073708;
116 gamma0[7] :== 9.24437796;
117 gamma0[8] :== 27.5075105;
118
119 n[1..56] IS_A real_constant;
120 n[1] :== 0.12533547935523e-1; n[2] :== 0.78957634722828e+1; n[3] :== -0.87803203303561e+1; n[4] :== 0.31802509345418e+0;
121 n[5] :== -0.26145533859358e+0; n[6] :== -0.78199751687981e-2; n[7] :== 0.88089493102134e-2;
122 n[8] :== -0.66856572307965e+0; n[9] :== 0.20433810950965e+0; n[10] :== -0.66212605039687e-4;
123 n[11] :== -0.19232721156002e+0; n[12] :== -0.25709043003438e+0; n[13] :== 0.16074868486251e+0; n[14] :== -0.40092828925807e-1;
124 n[15] :== 0.39343422603254e-6; n[16] :== -0.75941377088144e-5; n[17] :== 0.56250979351888e-3;
125 n[18] :== -0.15608652257135e-4; n[19] :== 0.11537996422951e-8; n[20] :== 0.36582165144204e-6;
126 n[21] :== -0.13251180074668e-11; n[22] :== -0.62639586912454e-9; n[23] :== -0.10793600908932e+0; n[24] :== 0.17611491008752e-1;
127 n[25] :== 0.22132295167546e+0; n[26] :== -0.40247669763528e+0; n[27] :== 0.58083399985759e+0;
128 n[28] :== 0.49969146990806e-2; n[29] :== -0.31358700712549e-1; n[30] :== -0.74315929710341e+0;
129 n[31] :== 0.47807329915480e+0; n[32] :== 0.20527940895948e-1; n[33] :== -0.13636435110343e+0; n[34] :== 0.14180634400617e-1;
130 n[35] :== 0.83326504880713e-2; n[36] :== -0.29052336009585e-1; n[37] :== 0.38615085574206e-1;
131 n[38] :== -0.20393486513704e-1; n[39] :== -0.16554050063734e-2; n[40] :== 0.19955571979541e-2;
132 n[41] :== 0.15870308324157e-3; n[42] :== -0.16388568342530e-4; n[43] :== 0.43613615723811e-1; n[44] :== 0.34994005463765e-1;
133 n[45] :== -0.76788197844621e-1; n[46] :== 0.22446277332006e-1; n[47] :== -0.62689710414685e-4;
134 n[48] :== -0.55711118565645e-9; n[49] :== -0.19905718354408e+0; n[50] :== 0.31777497330738e+0;
135 n[51] :== -0.11841182425981e+0; n[52] :== -0.31306260323435e+2; n[53] :== 0.31546140237781e+2; n[54] :== -0.25213154341695e+4;
136 n[55] :== -0.14874640856724e+0; n[56] :== 0.31806110878444e+0;
137
138 c[1..51] IS_A integer_constant;
139 c[1] :== 0; c[2] :== 0; c[3] :== 0; c[4] :== 0; c[5] :== 0;
140 c[6] :== 0; c[7] :== 0; c[8] :== 1; c[9] :== 1; c[10] :== 1;
141 c[11] :== 1; c[12] :== 1; c[13] :== 1; c[14] :== 1; c[15] :== 1;
142 c[16] :== 1; c[17] :== 1; c[18] :== 1; c[19] :== 1; c[20] :== 1;
143 c[21] :== 1; c[22] :== 1; c[23] :== 2; c[24] :== 2; c[25] :== 2;
144 c[26] :== 2; c[27] :== 2; c[28] :== 2; c[29] :== 2; c[30] :== 2;
145 c[31] :== 2; c[32] :== 2; c[33] :== 2; c[34] :== 2; c[35] :== 2;
146 c[36] :== 2; c[37] :== 2; c[38] :== 2; c[39] :== 2; c[40] :== 2;
147 c[41] :== 2; c[42] :== 2; c[43] :== 3; c[44] :== 3; c[45] :== 3;
148 c[46] :== 3; c[47] :== 4; c[48] :== 6; c[49] :== 6; c[50] :== 6;
149 c[51] :== 6;
150
151 d[1..54] IS_A integer_constant;
152 d[1] :== 1; d[2] :== 1; d[3] :== 1; d[4] :== 2; d[5] :== 2;
153 d[6] :== 3; d[7] :== 4; d[8] :== 1; d[9] :== 1; d[10] :== 1;
154 d[11] :== 2; d[12] :== 2; d[13] :== 3; d[14] :== 4; d[15] :== 4;
155 d[16] :== 5; d[17] :== 7; d[18] :== 9; d[19] :== 10; d[20] :== 11;
156 d[21] :== 13; d[22] :== 15; d[23] :== 1; d[24] :== 2; d[25] :== 2;
157 d[26] :== 2; d[27] :== 3; d[28] :== 4; d[29] :== 4; d[30] :== 4;
158 d[31] :== 5; d[32] :== 6; d[33] :== 6; d[34] :== 7; d[35] :== 9;
159 d[36] :== 9; d[37] :== 9; d[38] :== 9; d[39] :== 9; d[40] :== 10;
160 d[41] :== 10; d[42] :== 12; d[43] :== 3; d[44] :== 4; d[45] :== 4;
161 d[46] :== 5; d[47] :== 14; d[48] :== 3; d[49] :== 6; d[50] :== 6;
162 d[51] :== 6; d[52] :== 3; d[53] :== 3; d[54] :== 3;
163
164 t[1..54] IS_A real_constant;
165 t[1] :== -0.5; t[2] :== 0.875; t[3] :== 1; t[4] :== 0.5; t[5] :== 0.75;
166 t[6] :== 0.375; t[7] :== 1; t[8] :== 4; t[9] :== 6; t[10] :== 12;
167 t[11] :== 1; t[12] :== 5; t[13] :== 4; t[14] :== 2; t[15] :== 13;
168 t[16] :== 9; t[17] :== 3; t[18] :== 4; t[19] :== 11; t[20] :== 4;
169 t[21] :== 13; t[22] :== 1; t[23] :== 7; t[24] :== 1; t[25] :== 9;
170 t[26] :== 10; t[27] :== 10; t[28] :== 3; t[29] :== 7; t[30] :== 10;
171 t[31] :== 10; t[32] :== 6; t[33] :== 10; t[34] :== 10; t[35] :== 1;
172 t[36] :== 2; t[37] :== 3; t[38] :== 4; t[39] :== 8; t[40] :== 6;
173 t[41] :== 9; t[42] :== 8; t[43] :== 16; t[44] :== 22; t[45] :== 23;
174 t[46] :== 23; t[47] :== 10; t[48] :== 50; t[49] :== 44; t[50] :== 46;
175 t[51] :== 50; t[52] :== 0; t[53] :== 1; t[54] :== 4;
176
177 (* Correlation parameters *)
178
179 (* TODO convert from C to ASCEND arrays? Note trickiness with 0- and 1-based array indices. *)
180
181 a[range_r4] IS_A real_constant;
182 a[55]:==3.5;
183 a[56]:==3.5;
184
185 b[range_r4] IS_A real_constant;
186 b[55]:== 0.85;
187 b[56]:== 0.95;
188
189 A[range_r4] IS_A real_constant;
190 A[55]:==0.32;
191 A[56]:==0.32;
192
193 B[range_r4] IS_A real_constant;
194 B[55]:==0.2;
195 B[56]:==0.2;
196
197 C[range_r4] IS_A real_constant;
198 C[55]:==28;
199 C[56]:==32;
200
201 D[range_r4] IS_A real_constant;
202 D[55]:==700;
203 D[56]:==800;
204
205 beta_r4[range_r4] IS_A real_constant;
206 beta_r4[55]:==0.3;
207 beta_r4[56]:==0.3;
208
209 alpha[range_r3] IS_A integer_constant;
210 alpha[52]:==20;
211 alpha[53]:==20;
212 alpha[54]:==20;
213
214 beta[range_r3] IS_A real_constant;
215 beta[52]:==150;
216 beta[53]:==150;
217 beta[54]:==250;
218
219 gamma[range_r3] IS_A real_constant;
220 gamma[52]:==1.21;
221 gamma[53]:==1.21;
222 gamma[54]:==1.25;
223
224 (*--------------- DIMENSIONLESS PARTIAL DERIVATIVES ---------------- *)
225
226 (*------------ IDEAL PARTS ------------*)
227
228 phi0 IS_A factor;
229 phi0_expr: phi0 =
230 SUM[ n0[i]*ln(1-exp(-tau*gamma0[i])) | i IN [range_01] ]
231 + ln(delta) + n0[1] + n0[2]*tau + n0[3]*ln(tau);
232
233 phi0delta IS_A factor;
234 phi0delta_expr: phi0delta = 1.0/delta;
235
236 phi0deltadelta IS_A factor;
237 phi0deltadelta_expr: phi0deltadelta =
238 -1.0/(delta*delta);
239
240
241 phi0tau IS_A factor;
242 phi0tau_expr: phi0tau =
243 n0[2] + n0[3]/tau
244 + SUM[ n0[i]*gamma0[i]*(1/(1-exp(-tau*gamma0[i])) - 1) | i IN [range_01] ];
245
246 phi0deltatau IS_A real_constant;
247 phi0deltatau :== 0.0;
248
249 phi0tautau IS_A factor;
250 phi0tautau_expr: phi0tautau
251 = -n0[3] / tau^2
252 - SUM [ n0[i] * gamma0[i]^2 * exp(-gamma0[i] * tau) / ( 1 - exp(-gamma0[i] * tau) )^2 | i IN range_01 ];
253
254 (*----------- 'REAL' PARTS -- CLOSE YOUR EYES -----------*)
255
256 d1 IS_A factor;
257 d1_expr: d1 = delta - 1;
258
259 t1 IS_A factor;
260 t1_expr: t1 = tau -1;
261
262 r3_b1[range_r3] IS_A factor;
263 FOR i IN range_r3 CREATE
264 z_r3_b1[i]: r3_b1[i] = -alpha[i]*d1^2 - beta[i]* (tau - gamma[i])^2;
265 END FOR;
266
267 PSI[range_r4] IS_A factor;
268 theta[range_r4] IS_A factor;
269 DELTA[range_r4] IS_A factor;
270 dDELTA_ddelta[range_r4] IS_A factor;
271 dPSI_ddelta[range_r4] IS_A factor;
272 dDELTAbi_ddelta[range_r4] IS_A factor;
273 dDELTAbi_dtau[range_r4] IS_A factor;
274 dPSI_dtau[range_r4] IS_A factor;
275 d2DELTA_ddelta2[range_r4] IS_A factor;
276 d2DELTAbi_ddelta2[range_r4] IS_A factor;
277 d2PSI_ddelta2[range_r4] IS_A factor;
278 d2DELTAbi_dtau2[range_r4] IS_A factor;
279 d2PSI_dtau2[range_r4] IS_A factor;
280 d2PSI_ddeltadtau[range_r4] IS_A factor;
281 d2DELTAbi_ddeltadtau[range_r4] IS_A factor;
282
283 FOR i IN range_r4 CREATE
284 z_PSI[i]: PSI[i] = exp(-C[i]*d1^2 - D[i]*t1^2);
285 z_theta[i]: theta[i] = 1 - tau + A[i] * (d1^2)^( 1/(2 * beta_r4[i]) );
286 z_DELTA[i]: DELTA[i] = theta[i]^2 + B[i] * (d1^2)^a[i];
287
288 z_dDELTA_ddelta[i]: dDELTA_ddelta[i] = d1*(A[i]*theta[i]*2/beta_r4[i]*
289 (d1^2)^(1/(2*beta_r4[i]) - 1) + 2*B[i]*a[i]*
290 (d1^2)^(a[i] - 1) );
291
292 z_dPSI_ddelta[i]: dPSI_ddelta[i] = -2*C[i]*d1*PSI[i];
293
294 z_dDELTAbi_ddelta[i]: dDELTAbi_ddelta[i] = b[i] * DELTA[i]^(b[i]-1) * dDELTA_ddelta[i];
295
296 z_dDELTAbi_dtau[i]: dDELTAbi_dtau[i] = -2 * theta[i] * b[i] * DELTA[i]^(b[i]-1);
297 z_dPSI_dtau[i]: dPSI_dtau[i] = -2 * D[i] * t1 * PSI[i];
298
299 z_d2PSI_ddelta2[i]: d2PSI_ddelta2[i] = (2 * C[i] * d1^2 - 1) * 2 * C[i] * PSI[i];
300
301 z_d2DELTA_ddelta2[i]: d2DELTA_ddelta2[i] = 1/d1*dDELTA_ddelta[i] + d1^2*(4*B[i]*a[i]*
302 (a[i]-1)*(d1^2)^(a[i]-2) + 2*A[i]^2*
303 (1/(beta_r4[i]^2))*((d1^2) ^ (1/(2*beta_r4[i])-1))^2) +
304 A[i]*theta[i]*4/beta_r4[i]*(1/(2*beta_r4[i]) - 1)*
305 (d1^2)^(1/(2*beta_r4[i]) - 2);
306
307 z_d2DELTAbi_ddelta2[i]: d2DELTAbi_ddelta2[i] =b[i] * (
308 DELTA[i]^(b[i]-1) *d2DELTA_ddelta2[i]
309 + (b[i]-1) * DELTA[i]^(b[i]-2) * dDELTA_ddelta[i]^2 );
310
311 z_d2DELTAbi_dtau2[i]: d2DELTAbi_dtau2[i] = 2 * b[i] * DELTA[i]^(b[i]-1) +
312 4 * theta[i]^2 * b[i]*(b[i]-1) * DELTA[i]^(b[i]-2);
313
314 z_d2PSI_dtau2[i]: d2PSI_dtau2[i] = (2 * D[i] * t1^2 - 1) * 2 * D[i] * PSI[i];
315
316 z_d2PSI_ddeltadtau[i]: d2PSI_ddeltadtau[i] = 4 * C[i] * D[i] * d1 * t1 * PSI[i];
317
318 z_d2DELTAbi_ddeltadtau[i]: d2DELTAbi_ddeltadtau[i] =
319 - A[i] * b[i]* (2/beta_r4[i]) * DELTA[i]^(b[i]-1) * d1 * (d1^2)^(1/(2*beta_r4[i])-1)
320 - 2 * theta[i] * b[i] *(b[i]-1) * DELTA[i]^(b[i]-2) * dDELTA_ddelta[i];
321
322 END FOR;
323
324 phir_r2[range_r2] IS_A factor;
325 FOR i IN range_r2 CREATE
326 z_phir_r2[i]: phir_r2[i] = n[i] * delta^d[i] * tau^t[i] *
327 exp(-delta^c[i]);
328 END FOR;
329
330 phir IS_A factor;
331 phir_expr: phir
332 =
333 SUM[ n[i] * delta^d[i] * tau^t[i] | i IN [range_r1] ]
334 + SUM[ phir_r2[i] | i IN [range_r2] ]
335 + SUM[ n[i] * delta^d[i] * tau^t[i] * exp(
336 -alpha[i]*d1^2 - beta[i]*
337 (tau - gamma[i])*(tau - gamma[i])
338 ) | i IN [range_r3] ]
339 + SUM[ n[i] * DELTA[i]^b[i] * delta * PSI[i] | i IN [range_r4] ];
340
341
342 phirdelta_r2[range_r2] IS_A factor;
343 FOR i IN range_r2 CREATE
344 z_phirdelta_r2[i]: phirdelta_r2[i] = n[i] * exp(-delta^c[i]) * (delta^(d[i]-1) *
345 tau^t[i] * (d[i] - c[i] * delta^c[i]) );
346 END FOR;
347
348 phirdelta IS_A factor;
349 phirdelta_expr: phirdelta =
350 SUM[ n[i] * d[i] * delta^(d[i] - 1) * tau^t[i] | i IN range_r1 ]
351 + SUM[ phirdelta_r2[i] | i IN range_r2 ]
352 + SUM[ n[i]*delta^d[i] * tau^t[i] * exp( r3_b1[i] ) *
353 (d[i]/delta - 2*alpha[i]*d1) | i IN range_r3 ]
354 + SUM[ n[i] * (
355 DELTA[i]^b[i] * (PSI[i] + delta * dPSI_ddelta[i] )
356 + dDELTAbi_ddelta[i] * delta * PSI[i] ) | i IN range_r4 ];
357
358 phirtau_r2[range_r2] IS_A factor;
359 FOR i IN range_r2 CREATE
360 z_phirtau_r2[i]: phirtau_r2[i] = n[i]*t[i]*delta^d[i]*tau^(t[i]-1) * exp(-delta^c[i]);
361 END FOR;
362
363 phirtau IS_A factor;
364 phirtau_expr: phirtau =
365 SUM[ n[i] * t[i] * delta^d[i] * tau^(t[i]-1) | i IN range_r1 ]
366 + SUM[ phirtau_r2[i] | i IN range_r2 ]
367 + SUM[ n[i] * delta^d[i] * tau^t[i] * exp( r3_b1[i] )*
368 (t[i]/tau - 2*beta[i]*(tau - gamma[i])) | i IN range_r3 ]
369 + SUM[ n[i]*delta*(dDELTAbi_dtau[i] * PSI[i] +
370 DELTA[i]^b[i]*dPSI_dtau[i]) | i IN range_r4 ];
371
372 phirdeltadelta_r2[range_r2] IS_A factor;
373 FOR i IN range_r2 CREATE
374 z_phirdeltadelta_r2[i]: phirdeltadelta_r2[i] =
375 n[i] * exp(-delta^c[i]) * (
376 delta^(d[i]-2) *
377 tau^t[i] * (
378 ( d[i] - c[i]*delta^c[i] )*(d[i]- 1 - c[i] * delta^c[i] )
379 - c[i]^2 * delta^c[i]
380 )
381 );
382 END FOR;
383
384 phirdeltadelta IS_A factor;
385 phirdeltadelta_expr: phirdeltadelta =
386 SUM[
387 n[i] * d[i] * (d[i] - 1) * delta^(d[i]-2) * tau^t[i] | i IN range_r1 ]
388 + SUM [phirdeltadelta_r2[i] | i IN range_r2 ]
389 + SUM [
390 n[i]* tau^t[i] * exp( r3_b1[i] ) * (
391 - 2 * alpha[i] * delta^d[i]
392 + 4 * alpha[i]*alpha[i] * delta^d[i] * d1^2 (* d1 = d-1 = d-eps *)
393 - 4 * d[i] * alpha[i] * delta^(d[i]-1) * d1
394 + d[i] * (d[i]-1) * delta^(d[i]-1)
395 ) | i IN range_r3 ]
396 + SUM [n[i]*(
397 DELTA[i]^b[i] * (2*dPSI_ddelta[i] + delta*d2PSI_ddelta2[i] )
398 + 2*dDELTAbi_ddelta[i] * (PSI[i] + delta*dPSI_ddelta[i])
399 + d2DELTAbi_ddelta2[i] * delta * PSI[i] ) | i IN range_r4 ];
400
401
402
403 phirtautau_r2[range_r2] IS_A factor;
404 FOR i IN range_r2 CREATE
405 z_phirtautau_r2[i]: phirtautau_r2[i] = n[i]*t[i]*(t[i]-1)* delta^d[i] *
406 tau^(t[i]-2) * exp(- delta^c[i] );
407 END FOR;
408
409 phirtautau IS_A factor;
410 phirtautau_expr: phirtautau =
411 SUM[ n[i] * t[i] * (t[i]-1) * delta^d[i] *
412 tau^(t[i]-2) | i IN range_r1 ]
413 + SUM[ phirtautau_r2[i] | i IN range_r2 ]
414 + SUM[ n[i] * delta^d[i] * tau^t[i] * exp( r3_b1[i] ) *
415 ( (t[i]/tau - 2*beta[i]* (tau - gamma[i]) )^2 -
416 t[i]/tau^2 - 2*beta[i] ) | i IN range_r3 ]
417 + SUM[ n[i] * delta * (d2DELTAbi_dtau2[i] * PSI[i] +
418 2 * dDELTAbi_dtau[i] * dPSI_dtau[i]
419 + DELTA[i]^b[i] * d2PSI_dtau2[i] ) | i IN range_r4 ];
420
421
422 phirdeltatau_r2[range_r2] IS_A factor;
423 FOR i IN range_r2 CREATE
424 z_phirdeltatau_r2[i]: phirdeltatau_r2[i] =
425 n[i]*t[i] * delta^(d[i]-1) * tau^(t[i]-1)
426 * (d[i] - c[i] * delta^c[i]) * exp(-delta^c[i]);
427 END FOR;
428
429 phirdeltatau IS_A factor;
430 phirdeltatau_expr: phirdeltatau =
431 SUM[ n[i]*d[i]*t[i] * delta^(d[i]-1) * tau^(t[i]-1) | i IN range_r1 ]
432
433 + SUM[ phirdeltatau_r2[i] | i IN range_r2 ]
434 + SUM[ n[i] * delta^d[i] * tau^t[i] * exp( r3_b1[i] ) *
435 (d[i]/delta - 2*alpha[i]*d1)
436 * (t[i]/tau - 2*beta[i]* (tau - gamma[i]) ) | i IN range_r3 ]
437 + SUM[ n[i]*( DELTA[i]^b[i] * (dPSI_dtau[i] + delta *
438 d2PSI_ddeltadtau[i] ) + delta * dDELTAbi_ddelta[i] * dPSI_dtau[i] +
439 dDELTAbi_dtau[i] * (PSI[i] + delta * dPSI_ddelta[i] ) +
440 d2DELTAbi_ddeltadtau[i] * delta * PSI[i] ) | i IN range_r4 ];
441
442 (*--------- THERMO PROPERTY RELATIONS ----------- *)
443
444 za_pressure: p
445 = rho * R * T * (1 + delta*phirdelta);
446
447 za_internal_energy: u
448 = R * T * tau * (phi0tau + phirtau);
449
450 za_enthalpy: h
451 = R * T * (1 + tau*(phi0tau + phirtau) + delta*phirdelta);
452
453 za_entropy: s
454 = R * (tau*(phi0tau + phirtau) - phi0 - phir);
455
456 za_c_isochoric: cv
457 = - R * tau^2 * (phi0tautau + phirtautau);
458
459 za_c_isobaric: cp
460 = - R * (
461 tau^2 * (phi0tautau + phirtautau)
462 + ( ( 1 + delta*phirdelta - delta*tau*phirdeltatau )^2 )
463 / ( 1 + 2*delta*phirdelta + delta^2 * phirdeltadelta )
464 );
465
466 (* spd_sound: w^2
467 = R * T * ( 1 + 2*delta*phirdelta + delta^2 * phirdeltadelta -
468 ( 1 + delta*phirdelta - delta*tau*phirdeltatau )^2
469 / ( tau^2 * (phi0tautau + phirtautau) )
470 );
471 *)
472
473 METHODS
474 METHOD default_self;
475 RUN reset;
476 RUN values;
477 END default_self;
478
479 METHOD specify;
480 FIX T,rho;
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 ATOM density_slack REFINES mass_density;
508 lower_bound := 0.0 {kg/m^3};
509 END density_slack;
510
511 MODEL iapws95_2phase REFINES thermo_state;
512 (* phase components *)
513 Sl IS_A iapws95_1phase;
514 Sv IS_A iapws95_1phase;
515
516 (* saturation conditions *)
517 sat IS_A iapws_sat_density;
518
519 T, Sl.T, Sv.T, sat.T ARE_THE_SAME;
520 p, Sl.p ARE_THE_SAME;
521
522 rhof ALIASES sat.rhof;
523 rhog ALIASES sat.rhog;
524 rhov ALIASES Sv.rho;
525 rhol ALIASES Sl.rho;
526
527 rhofslack, rhogslack IS_A density_slack;
528 z02a: rhogslack = rhog - rhov; (* rhog >= rhov since gas is most dense at saturation*)
529 z03a: rhofslack = rhol - rhof; (* rhol >= rhof since it's compressed liquid. but what about the density anomaly..??? *)
530
531 xf IS_A fraction;
532 xg IS_A fraction;
533 z01: xf + xg = 1.0;
534 z02: xf * rhogslack = 0.0 {kg/m^3}; (* xf==0 or rhog==rhov *)
535 z03: xg * rhofslack = 0.0 {kg/m^3}; (* xg==0 or rhof==rhol *)
536 x ALIASES xg;
537
538 (* two-phase properties *)
539 z04: Sl.rho * Sv.rho = rho * ( xf * Sv.rho + xg * Sl.rho);
540 z05: u = xf * Sl.u + xg * Sv.u;
541 z06: h = xf * Sl.h + xg * Sv.h;
542 z07: s = xf * Sl.s + xg * Sv.s;
543 z08: cp = xf * Sl.cp + xg * Sv.cp;
544 z09: cv = xf * Sl.cv + xg * Sv.cv;
545 (*
546 z10: w = xf * Sl.w + xg * Sv.w; * check this *
547 *)
548 METHODS
549
550 METHOD default_self;
551 RUN reset;
552 RUN values;
553 Sl.rho := 800 {kg/m^3};
554 Sl.delta := 2.48;
555 Sl.d1 := 1.48;
556 Sv.delta := 0.012;
557 Sv.d1 := -0.987;
558 RUN scale_self;
559 RUN bound_self;
560 END default_self;
561
562 METHOD specify;
563 FIX T,x;
564 END specify;
565
566 METHOD values;
567 T := 390 {K};
568 x := 0.5;
569 END values;
570
571 METHOD scale_self;
572 Sl.rho.nominal := 800 {kg/m^3};
573 Sv.rho.nominal := 1 {kg/m^3};
574 END scale_self;
575
576 METHOD bound_self;
577 (*Sl.tau.lower_bound := 1; *)(* T = T_c *)
578 (*Sl.tau.upper_bound := 2.3690133; *)(* T = 273.15 *)
579
580 Sl.d1.lower_bound := -1;
581 Sv.d1.lower_bound := -1;
582 END bound_self;
583
584 END iapws95_2phase;
585
586
587 MODEL iapws95_2phase_ph_test REFINES iapws95_2phase;
588
589 METHODS
590 METHOD specify;
591 FIX p,h;
592 END specify;
593 METHOD values;
594 p := 10 {bar};
595 h := 2998 {kJ/kg};
596 END values;
597
598 METHOD self_test;
599 ASSERT abs(T - 548.15 {K}) < 1 {K};
600 ASSERT abs(rho - 8.43 {kg/m^3}) < 0.2 {kg/m^3};
601 END self_test;
602
603 END iapws95_2phase_ph_test;
604
605 MODEL iapws95_2phase_ph_test2 REFINES iapws95_2phase;
606
607 METHODS
608 METHOD specify;
609 FIX p,h;
610 END specify;
611 METHOD values;
612 p := 10 {bar};
613 h := 418.7 {kJ/kg};
614 END values;
615 METHOD self_test;
616 ASSERT abs(T - 373.15 {K}) < 1 {K};
617 ASSERT abs(rho - 959 {kg/m^3}) < 3 {kg/m^3};
618 END self_test;
619
620 END iapws95_2phase_ph_test2;
621
622 (*
623 This last model *will* converge with careful coaxing
624 You need to alternate between ROW2NORM scaling
625 and RELNOM scaling using the QRSlv engine.
626 Some tweaking of delta/d1 values may be necessary as
627 well.
628 *)
629 MODEL iapws95_2phase_ph_test3 REFINES iapws95_1phase;
630
631 METHODS
632 METHOD specify;
633 FIX p,h;
634 END specify;
635 METHOD values;
636 p := 1 {bar};
637 h := 104.9 {kJ/kg};
638 (* free vars *)
639 delta := 3;
640 d1 := 2;
641 rho := 900;
642 END values;
643 METHOD self_test;
644 ASSERT abs(T - 298.15 {K}) < 1 {K};
645 ASSERT abs(rho - 997 {kg/m^3}) < 3 {kg/m^3};
646 END self_test;
647
648 END iapws95_2phase_ph_test3;

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