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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 120 - (show annotations) (download) (as text)
Mon Dec 19 06:07:52 2005 UTC (14 years, 7 months ago) by johnpye
File MIME type: text/x-ascend
File size: 18030 byte(s)
Adding demonstration to bug #165.
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 steam_state;
58
59 delta IS_A factor;
60 tau IS_A factor;
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[i] =
396 n[i]* tau^t[i] * exp( r3_b1[i] ) * (
397 - 2 * alpha[i] * delta^d[i]
398 (* >>>>>>>>> *) + 4 * alpha[i]^2 * delta^d[i] * d1^2 (* <<<<<<<<<<<<<<< BUG IS HERE! *)
399 - 4 * d[i] * alpha[i] * delta^(d[i]-1) * d1
400 + d[i] * (d[i]-1) * delta^(d[i]-2)
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
494 = sqrt(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 steam_state;

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