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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 327 - (hide 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 aw0a 179 REQUIRE "system.a4l";
2     REQUIRE "atoms.a4l";
3     REQUIRE "johnpye/thermo_types.a4c";
4 johnpye 112 REQUIRE "johnpye/iapws_sat_curves.a4c";
5 johnpye 78
6     (*************************************************************************
7    
8     The thermodynamic properties of water calculated with the
9 johnpye 90 IAPWS95 equations. Variables and (example/possible) units:
10 johnpye 78
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 johnpye 167 http://freesteam.sourceforge.net/
39 johnpye 78
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 johnpye 107
59 johnpye 108 MODEL iapws95_1phase REFINES thermo_state;
60 johnpye 78
61 johnpye 105 delta IS_A positive_variable;
62     tau IS_A positive_variable;
63 johnpye 78
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 johnpye 84
74     R IS_A specific_gas_constant;
75     R "specific gas constant for water"
76     :== 0.46151805 {kJ/kg/K};
77 johnpye 78
78 johnpye 327 z_tau: tau = Tc / T;
79     z_delta: delta = rho / rhoc;
80 johnpye 78
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 johnpye 84 phi0_expr: phi0 =
230 johnpye 96 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 johnpye 78
233     phi0delta IS_A factor;
234 johnpye 84 phi0delta_expr: phi0delta = 1.0/delta;
235 johnpye 78
236     phi0deltadelta IS_A factor;
237 johnpye 84 phi0deltadelta_expr: phi0deltadelta =
238 johnpye 78 -1.0/(delta*delta);
239    
240    
241     phi0tau IS_A factor;
242 johnpye 84 phi0tau_expr: phi0tau =
243 johnpye 78 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 johnpye 84 phi0tautau_expr: phi0tautau
251 johnpye 96 = -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 johnpye 78
254     (*----------- 'REAL' PARTS -- CLOSE YOUR EYES -----------*)
255    
256     d1 IS_A factor;
257 johnpye 84 d1_expr: d1 = delta - 1;
258 johnpye 78
259     t1 IS_A factor;
260 johnpye 84 t1_expr: t1 = tau -1;
261 johnpye 78
262     r3_b1[range_r3] IS_A factor;
263     FOR i IN range_r3 CREATE
264 johnpye 327 z_r3_b1[i]: r3_b1[i] = -alpha[i]*d1^2 - beta[i]* (tau - gamma[i])^2;
265 johnpye 78 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 johnpye 327 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 johnpye 78
288 johnpye 327 z_dDELTA_ddelta[i]: dDELTA_ddelta[i] = d1*(A[i]*theta[i]*2/beta_r4[i]*
289 johnpye 78 (d1^2)^(1/(2*beta_r4[i]) - 1) + 2*B[i]*a[i]*
290     (d1^2)^(a[i] - 1) );
291    
292 johnpye 327 z_dPSI_ddelta[i]: dPSI_ddelta[i] = -2*C[i]*d1*PSI[i];
293 johnpye 78
294 johnpye 327 z_dDELTAbi_ddelta[i]: dDELTAbi_ddelta[i] = b[i] * DELTA[i]^(b[i]-1) * dDELTA_ddelta[i];
295 johnpye 78
296 johnpye 327 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 johnpye 78
299 johnpye 327 z_d2PSI_ddelta2[i]: d2PSI_ddelta2[i] = (2 * C[i] * d1^2 - 1) * 2 * C[i] * PSI[i];
300 johnpye 78
301 johnpye 327 z_d2DELTA_ddelta2[i]: d2DELTA_ddelta2[i] = 1/d1*dDELTA_ddelta[i] + d1^2*(4*B[i]*a[i]*
302 johnpye 78 (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 johnpye 327 z_d2DELTAbi_ddelta2[i]: d2DELTAbi_ddelta2[i] =b[i] * (
308 johnpye 96 DELTA[i]^(b[i]-1) *d2DELTA_ddelta2[i]
309     + (b[i]-1) * DELTA[i]^(b[i]-2) * dDELTA_ddelta[i]^2 );
310 johnpye 78
311 johnpye 327 z_d2DELTAbi_dtau2[i]: d2DELTAbi_dtau2[i] = 2 * b[i] * DELTA[i]^(b[i]-1) +
312 johnpye 96 4 * theta[i]^2 * b[i]*(b[i]-1) * DELTA[i]^(b[i]-2);
313 johnpye 78
314 johnpye 327 z_d2PSI_dtau2[i]: d2PSI_dtau2[i] = (2 * D[i] * t1^2 - 1) * 2 * D[i] * PSI[i];
315 johnpye 78
316 johnpye 327 z_d2PSI_ddeltadtau[i]: d2PSI_ddeltadtau[i] = 4 * C[i] * D[i] * d1 * t1 * PSI[i];
317 johnpye 78
318 johnpye 327 z_d2DELTAbi_ddeltadtau[i]: d2DELTAbi_ddeltadtau[i] =
319 johnpye 96 - 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 johnpye 78
322     END FOR;
323    
324 johnpye 85 phir_r2[range_r2] IS_A factor;
325     FOR i IN range_r2 CREATE
326 johnpye 327 z_phir_r2[i]: phir_r2[i] = n[i] * delta^d[i] * tau^t[i] *
327 johnpye 85 exp(-delta^c[i]);
328     END FOR;
329    
330 johnpye 78 phir IS_A factor;
331 johnpye 84 phir_expr: phir
332 johnpye 78 =
333 johnpye 84 SUM[ n[i] * delta^d[i] * tau^t[i] | i IN [range_r1] ]
334 johnpye 85 + SUM[ phir_r2[i] | i IN [range_r2] ]
335 johnpye 84 + SUM[ n[i] * delta^d[i] * tau^t[i] * exp(
336 johnpye 78 -alpha[i]*d1^2 - beta[i]*
337     (tau - gamma[i])*(tau - gamma[i])
338 johnpye 84 ) | i IN [range_r3] ]
339     + SUM[ n[i] * DELTA[i]^b[i] * delta * PSI[i] | i IN [range_r4] ];
340 johnpye 78
341    
342 johnpye 85 phirdelta_r2[range_r2] IS_A factor;
343     FOR i IN range_r2 CREATE
344 johnpye 327 z_phirdelta_r2[i]: phirdelta_r2[i] = n[i] * exp(-delta^c[i]) * (delta^(d[i]-1) *
345 johnpye 85 tau^t[i] * (d[i] - c[i] * delta^c[i]) );
346     END FOR;
347    
348 johnpye 78 phirdelta IS_A factor;
349 johnpye 84 phirdelta_expr: phirdelta =
350 johnpye 78 SUM[ n[i] * d[i] * delta^(d[i] - 1) * tau^t[i] | i IN range_r1 ]
351 johnpye 85 + SUM[ phirdelta_r2[i] | i IN range_r2 ]
352 johnpye 78 + 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 johnpye 85 phirtau_r2[range_r2] IS_A factor;
359     FOR i IN range_r2 CREATE
360 johnpye 327 z_phirtau_r2[i]: phirtau_r2[i] = n[i]*t[i]*delta^d[i]*tau^(t[i]-1) * exp(-delta^c[i]);
361 johnpye 85 END FOR;
362    
363 johnpye 78 phirtau IS_A factor;
364 johnpye 84 phirtau_expr: phirtau =
365 johnpye 78 SUM[ n[i] * t[i] * delta^d[i] * tau^(t[i]-1) | i IN range_r1 ]
366 johnpye 85 + SUM[ phirtau_r2[i] | i IN range_r2 ]
367 johnpye 78 + 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 johnpye 84 DELTA[i]^b[i]*dPSI_dtau[i]) | i IN range_r4 ];
371 johnpye 78
372 johnpye 85 phirdeltadelta_r2[range_r2] IS_A factor;
373     FOR i IN range_r2 CREATE
374 johnpye 327 z_phirdeltadelta_r2[i]: phirdeltadelta_r2[i] =
375 johnpye 85 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 johnpye 84
384 johnpye 107 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 johnpye 96 n[i]* tau^t[i] * exp( r3_b1[i] ) * (
391     - 2 * alpha[i] * delta^d[i]
392 johnpye 104 + 4 * alpha[i]*alpha[i] * delta^d[i] * d1^2 (* d1 = d-1 = d-eps *)
393 johnpye 96 - 4 * d[i] * alpha[i] * delta^(d[i]-1) * d1
394 johnpye 105 + d[i] * (d[i]-1) * delta^(d[i]-1)
395 johnpye 107 ) | i IN range_r3 ]
396     + SUM [n[i]*(
397 johnpye 78 DELTA[i]^b[i] * (2*dPSI_ddelta[i] + delta*d2PSI_ddelta2[i] )
398     + 2*dDELTAbi_ddelta[i] * (PSI[i] + delta*dPSI_ddelta[i])
399 johnpye 107 + d2DELTAbi_ddelta2[i] * delta * PSI[i] ) | i IN range_r4 ];
400 johnpye 78
401 johnpye 84
402 johnpye 96
403 johnpye 85 phirtautau_r2[range_r2] IS_A factor;
404     FOR i IN range_r2 CREATE
405 johnpye 327 z_phirtautau_r2[i]: phirtautau_r2[i] = n[i]*t[i]*(t[i]-1)* delta^d[i] *
406 johnpye 85 tau^(t[i]-2) * exp(- delta^c[i] );
407     END FOR;
408    
409 johnpye 78 phirtautau IS_A factor;
410 johnpye 84 phirtautau_expr: phirtautau =
411 johnpye 78 SUM[ n[i] * t[i] * (t[i]-1) * delta^d[i] *
412     tau^(t[i]-2) | i IN range_r1 ]
413 johnpye 85 + SUM[ phirtautau_r2[i] | i IN range_r2 ]
414 johnpye 78 + 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 johnpye 84 2 * dDELTAbi_dtau[i] * dPSI_dtau[i]
419     + DELTA[i]^b[i] * d2PSI_dtau2[i] ) | i IN range_r4 ];
420 johnpye 78
421 johnpye 85
422     phirdeltatau_r2[range_r2] IS_A factor;
423     FOR i IN range_r2 CREATE
424 johnpye 327 z_phirdeltatau_r2[i]: phirdeltatau_r2[i] =
425 johnpye 85 n[i]*t[i] * delta^(d[i]-1) * tau^(t[i]-1)
426 johnpye 104 * (d[i] - c[i] * delta^c[i]) * exp(-delta^c[i]);
427 johnpye 85 END FOR;
428    
429 johnpye 78 phirdeltatau IS_A factor;
430 johnpye 84 phirdeltatau_expr: phirdeltatau =
431     SUM[ n[i]*d[i]*t[i] * delta^(d[i]-1) * tau^(t[i]-1) | i IN range_r1 ]
432 johnpye 78
433 johnpye 85 + SUM[ phirdeltatau_r2[i] | i IN range_r2 ]
434 johnpye 78 + SUM[ n[i] * delta^d[i] * tau^t[i] * exp( r3_b1[i] ) *
435 johnpye 84 (d[i]/delta - 2*alpha[i]*d1)
436     * (t[i]/tau - 2*beta[i]* (tau - gamma[i]) ) | i IN range_r3 ]
437 johnpye 78 + 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 johnpye 84 d2DELTAbi_ddeltadtau[i] * delta * PSI[i] ) | i IN range_r4 ];
441 johnpye 78
442     (*--------- THERMO PROPERTY RELATIONS ----------- *)
443    
444 johnpye 327 za_pressure: p
445 johnpye 84 = rho * R * T * (1 + delta*phirdelta);
446 johnpye 78
447 johnpye 327 za_internal_energy: u
448 johnpye 84 = R * T * tau * (phi0tau + phirtau);
449 johnpye 78
450 johnpye 327 za_enthalpy: h
451 johnpye 84 = R * T * (1 + tau*(phi0tau + phirtau) + delta*phirdelta);
452 johnpye 78
453 johnpye 327 za_entropy: s
454 johnpye 84 = R * (tau*(phi0tau + phirtau) - phi0 - phir);
455 johnpye 78
456 johnpye 327 za_c_isochoric: cv
457 johnpye 84 = - R * tau^2 * (phi0tautau + phirtautau);
458 johnpye 78
459 johnpye 327 za_c_isobaric: cp
460 johnpye 96 = - R * (
461     tau^2 * (phi0tautau + phirtautau)
462     + ( ( 1 + delta*phirdelta - delta*tau*phirdeltatau )^2 )
463     / ( 1 + 2*delta*phirdelta + delta^2 * phirdeltadelta )
464     );
465 johnpye 78
466 johnpye 296 (* spd_sound: w^2
467 johnpye 105 = R * T * ( 1 + 2*delta*phirdelta + delta^2 * phirdeltadelta -
468 johnpye 78 ( 1 + delta*phirdelta - delta*tau*phirdeltatau )^2
469 johnpye 84 / ( tau^2 * (phi0tautau + phirtautau) )
470 johnpye 105 );
471 johnpye 296 *)
472 johnpye 78
473     METHODS
474     METHOD default_self;
475 johnpye 327 RUN reset;
476 johnpye 83 RUN values;
477 johnpye 78 END default_self;
478    
479     METHOD specify;
480 johnpye 186 FIX T,rho;
481 johnpye 78 END specify;
482    
483     METHOD values;
484 johnpye 96 (* 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 johnpye 78 END values;
488    
489 johnpye 154 METHOD self_test;
490 johnpye 155 (* 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 johnpye 154 ASSERT abs(phi0delta - (0.384236747)) < 1e-9;
494 johnpye 155 ASSERT abs(phirdelta - (-0.364366650)) < 1e-9;
495     ASSERT abs(phi0deltadelta - (-0.147637878)) < 1e-9;
496 johnpye 154 ASSERT abs(phirdeltadelta - (0.856063701)) < 1e-9;
497     ASSERT abs(phi0tau - (0.904611106e1)) < 1e-8;
498 johnpye 155 ASSERT abs(phirtau - (-0.581403435e1)) < 1e-8;
499     ASSERT abs(phi0tautau - (-0.193249185e1)) < 1e-8;
500     ASSERT abs(phirtautau - (-0.223440737e1)) < 1e-8;
501 johnpye 154 ASSERT abs(phi0deltatau - 0) < 1e-20;
502 johnpye 155 ASSERT abs(phirdeltatau - (-0.112176915e1)) < 1e-8;
503 johnpye 154 END self_test;
504    
505 johnpye 108 END iapws95_1phase;
506 johnpye 78
507 johnpye 186 ATOM density_slack REFINES mass_density;
508     lower_bound := 0.0 {kg/m^3};
509     END density_slack;
510    
511 johnpye 108 MODEL iapws95_2phase REFINES thermo_state;
512 johnpye 166 (* phase components *)
513 johnpye 186 Sl IS_A iapws95_1phase;
514     Sv IS_A iapws95_1phase;
515 johnpye 166
516     (* saturation conditions *)
517 johnpye 186 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 johnpye 167 rhof ALIASES sat.rhof;
523     rhog ALIASES sat.rhog;
524 johnpye 186 rhov ALIASES Sv.rho;
525     rhol ALIASES Sl.rho;
526 johnpye 106
527 johnpye 186 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 johnpye 167
531 johnpye 186 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 johnpye 166 x ALIASES xg;
537 johnpye 106
538 johnpye 186 (* 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 johnpye 296 (*
546     z10: w = xf * Sl.w + xg * Sv.w; * check this *
547     *)
548 johnpye 106 METHODS
549    
550     METHOD default_self;
551 johnpye 324 RUN reset;
552 johnpye 106 RUN values;
553 johnpye 320 Sl.rho := 800 {kg/m^3};
554 johnpye 327 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 johnpye 106 END default_self;
561    
562     METHOD specify;
563 johnpye 186 FIX T,x;
564 johnpye 106 END specify;
565    
566     METHOD values;
567     T := 390 {K};
568     x := 0.5;
569     END values;
570 johnpye 327
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 johnpye 106
580 johnpye 327 Sl.d1.lower_bound := -1;
581     Sv.d1.lower_bound := -1;
582     END bound_self;
583    
584 johnpye 108 END iapws95_2phase;
585 johnpye 106
586 johnpye 324
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 johnpye 327
598 johnpye 324 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 johnpye 327
603 johnpye 324 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 johnpye 327 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;
649    

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