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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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