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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 186 - (hide annotations) (download) (as text)
Fri Jan 13 03:36:44 2006 UTC (14 years, 8 months ago) by johnpye
File MIME type: text/x-ascend
File size: 19668 byte(s)
Modularising the 'loop' components in my steam system test case.
Changing some models to use the new 'FIX varlist' syntax (bug 204)
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     tau = Tc / T;
79     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 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     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     PSI[i] = exp(-C[i]*d1^2 - D[i]*t1^2);
285 johnpye 84 (*theta_expr:*) theta[i] = 1 - tau + A[i] * (d1^2)^( 1/(2 * beta_r4[i]) );
286     (*DELTA_expr:*) DELTA[i] = theta[i]^2 + B[i] * (d1^2)^a[i];
287 johnpye 78
288 johnpye 84 (*dDELTA_ddelta_expr:*) 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 84 (*dPSI_ddelta_expr:*) dPSI_ddelta[i] = -2*C[i]*d1*PSI[i];
293 johnpye 78
294 johnpye 84 (*dDELTAbi_ddelta_expr:*) dDELTAbi_ddelta[i] = b[i] * DELTA[i]^(b[i]-1) * dDELTA_ddelta[i];
295 johnpye 78
296 johnpye 84 (*dDELTAbi_dtau_expr:*) dDELTAbi_dtau[i] = -2 * theta[i] * b[i] * DELTA[i]^(b[i]-1);
297     (*dPSI_dtau_expr:*) dPSI_dtau[i] = -2 * D[i] * t1 * PSI[i];
298 johnpye 78
299 johnpye 84 (*d2PSI_ddelta2_expr:*) d2PSI_ddelta2[i] = (2 * C[i] * d1^2 - 1) * 2 * C[i] * PSI[i];
300 johnpye 78
301 johnpye 84 (*d2DELTA_ddelta2_expr:*) 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 96 (*d2DELTAbi_ddelta2_expr:*) 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 johnpye 78
311 johnpye 96 (*d2DELTAbi_dtau2_expr:*) 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 johnpye 78
314 johnpye 84 (*d2PSI_dtau2_expr:*) d2PSI_dtau2[i] = (2 * D[i] * t1^2 - 1) * 2 * D[i] * PSI[i];
315 johnpye 78
316 johnpye 84 (*d2PSI_ddeltadtau_expr:*) d2PSI_ddeltadtau[i] = 4 * C[i] * D[i] * d1 * t1 * PSI[i];
317 johnpye 78
318 johnpye 96 (*d2DELTAbi_ddeltadtau_expr:*) 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 johnpye 78
322     END FOR;
323    
324 johnpye 85 phir_r2[range_r2] IS_A factor;
325     FOR i IN range_r2 CREATE
326     phir_r2[i] = n[i] * delta^d[i] * tau^t[i] *
327     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     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 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     phirtau_r2[i] = n[i]*t[i]*delta^d[i]*tau^(t[i]-1) * exp(-delta^c[i]);
361     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     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 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     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 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     phirdeltatau_r2[i] =
425     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     pressure: p
445 johnpye 84 = rho * R * T * (1 + delta*phirdelta);
446 johnpye 78
447     internal_energy: u
448 johnpye 84 = R * T * tau * (phi0tau + phirtau);
449 johnpye 78
450     enthalpy: h
451 johnpye 84 = R * T * (1 + tau*(phi0tau + phirtau) + delta*phirdelta);
452 johnpye 78
453     entropy: s
454 johnpye 84 = R * (tau*(phi0tau + phirtau) - phi0 - phir);
455 johnpye 78
456     c_isochoric: cv
457 johnpye 84 = - R * tau^2 * (phi0tautau + phirtautau);
458 johnpye 78
459     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 105 spd_sound: w^2
467     = 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 78
472     METHODS
473     METHOD default_self;
474 johnpye 83 RUN ClearAll;
475     RUN specify;
476     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     z10: w = xf * Sl.w + xg * Sv.w; (* check this *)
546 johnpye 106 METHODS
547    
548     METHOD default_self;
549     RUN ClearAll;
550     RUN specify;
551     RUN values;
552     END default_self;
553    
554     METHOD specify;
555 johnpye 186 FIX T,x;
556 johnpye 106 END specify;
557    
558     METHOD values;
559     T := 390 {K};
560     x := 0.5;
561     END values;
562    
563 johnpye 108 END iapws95_2phase;
564 johnpye 106

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