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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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