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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 169 - (hide annotations) (download) (as text)
Fri Jan 6 02:50:41 2006 UTC (14 years, 9 months ago) by johnpye
File MIME type: text/x-ascend
File size: 19344 byte(s)
Added failing version for bug #201.
1 johnpye 169 REQUIRE "system.a4l";
2     REQUIRE "atoms.a4l";
3     REQUIRE "johnpye/thermo_types.a4c";
4     REQUIRE "johnpye/iapws_sat_curves.a4c";
5    
6     (*************************************************************************
7    
8     The thermodynamic properties of water calculated with the
9     IAPWS95 equations. Variables and (example/possible) units:
10    
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     http://freesteam.sourceforge.net/
39    
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    
59     MODEL iapws95_1phase REFINES thermo_state;
60    
61     delta IS_A positive_variable;
62     tau IS_A positive_variable;
63    
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    
74     R IS_A specific_gas_constant;
75     R "specific gas constant for water"
76     :== 0.46151805 {kJ/kg/K};
77    
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     phi0_expr: phi0 =
230     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    
233     phi0delta IS_A factor;
234     phi0delta_expr: phi0delta = 1.0/delta;
235    
236     phi0deltadelta IS_A factor;
237     phi0deltadelta_expr: phi0deltadelta =
238     -1.0/(delta*delta);
239    
240    
241     phi0tau IS_A factor;
242     phi0tau_expr: phi0tau =
243     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     phi0tautau_expr: phi0tautau
251     = -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    
254     (*----------- 'REAL' PARTS -- CLOSE YOUR EYES -----------*)
255    
256     d1 IS_A factor;
257     d1_expr: d1 = delta - 1;
258    
259     t1 IS_A factor;
260     t1_expr: t1 = tau -1;
261    
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     (*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    
288     (*dDELTA_ddelta_expr:*) dDELTA_ddelta[i] = d1*(A[i]*theta[i]*2/beta_r4[i]*
289     (d1^2)^(1/(2*beta_r4[i]) - 1) + 2*B[i]*a[i]*
290     (d1^2)^(a[i] - 1) );
291    
292     (*dPSI_ddelta_expr:*) dPSI_ddelta[i] = -2*C[i]*d1*PSI[i];
293    
294     (*dDELTAbi_ddelta_expr:*) dDELTAbi_ddelta[i] = b[i] * DELTA[i]^(b[i]-1) * dDELTA_ddelta[i];
295    
296     (*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    
299     (*d2PSI_ddelta2_expr:*) d2PSI_ddelta2[i] = (2 * C[i] * d1^2 - 1) * 2 * C[i] * PSI[i];
300    
301     (*d2DELTA_ddelta2_expr:*) d2DELTA_ddelta2[i] = 1/d1*dDELTA_ddelta[i] + d1^2*(4*B[i]*a[i]*
302     (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     (*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    
311     (*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    
314     (*d2PSI_dtau2_expr:*) d2PSI_dtau2[i] = (2 * D[i] * t1^2 - 1) * 2 * D[i] * PSI[i];
315    
316     (*d2PSI_ddeltadtau_expr:*) d2PSI_ddeltadtau[i] = 4 * C[i] * D[i] * d1 * t1 * PSI[i];
317    
318     (*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    
322     END FOR;
323    
324     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     phir IS_A factor;
331     phir_expr: phir
332     =
333     SUM[ n[i] * delta^d[i] * tau^t[i] | i IN [range_r1] ]
334     + SUM[ phir_r2[i] | i IN [range_r2] ]
335     + SUM[ n[i] * delta^d[i] * tau^t[i] * exp(
336     -alpha[i]*d1^2 - beta[i]*
337     (tau - gamma[i])*(tau - gamma[i])
338     ) | i IN [range_r3] ]
339     + SUM[ n[i] * DELTA[i]^b[i] * delta * PSI[i] | i IN [range_r4] ];
340    
341    
342     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     phirdelta IS_A factor;
349     phirdelta_expr: phirdelta =
350     SUM[ n[i] * d[i] * delta^(d[i] - 1) * tau^t[i] | i IN range_r1 ]
351     + SUM[ phirdelta_r2[i] | i IN range_r2 ]
352     + 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     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     phirtau IS_A factor;
364     phirtau_expr: phirtau =
365     SUM[ n[i] * t[i] * delta^d[i] * tau^(t[i]-1) | i IN range_r1 ]
366     + SUM[ phirtau_r2[i] | i IN range_r2 ]
367     + 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     DELTA[i]^b[i]*dPSI_dtau[i]) | i IN range_r4 ];
371    
372     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    
384     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     n[i]* tau^t[i] * exp( r3_b1[i] ) * (
391     - 2 * alpha[i] * delta^d[i]
392     + 4 * alpha[i]*alpha[i] * delta^d[i] * d1^2 (* d1 = d-1 = d-eps *)
393     - 4 * d[i] * alpha[i] * delta^(d[i]-1) * d1
394     + d[i] * (d[i]-1) * delta^(d[i]-1)
395     ) | i IN range_r3 ]
396     + SUM [n[i]*(
397     DELTA[i]^b[i] * (2*dPSI_ddelta[i] + delta*d2PSI_ddelta2[i] )
398     + 2*dDELTAbi_ddelta[i] * (PSI[i] + delta*dPSI_ddelta[i])
399     + d2DELTAbi_ddelta2[i] * delta * PSI[i] ) | i IN range_r4 ];
400    
401    
402    
403     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     phirtautau IS_A factor;
410     phirtautau_expr: phirtautau =
411     SUM[ n[i] * t[i] * (t[i]-1) * delta^d[i] *
412     tau^(t[i]-2) | i IN range_r1 ]
413     + SUM[ phirtautau_r2[i] | i IN range_r2 ]
414     + 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     2 * dDELTAbi_dtau[i] * dPSI_dtau[i]
419     + DELTA[i]^b[i] * d2PSI_dtau2[i] ) | i IN range_r4 ];
420    
421    
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     * (d[i] - c[i] * delta^c[i]) * exp(-delta^c[i]);
427     END FOR;
428    
429     phirdeltatau IS_A factor;
430     phirdeltatau_expr: phirdeltatau =
431     SUM[ n[i]*d[i]*t[i] * delta^(d[i]-1) * tau^(t[i]-1) | i IN range_r1 ]
432    
433     + SUM[ phirdeltatau_r2[i] | i IN range_r2 ]
434     + SUM[ n[i] * delta^d[i] * tau^t[i] * exp( r3_b1[i] ) *
435     (d[i]/delta - 2*alpha[i]*d1)
436     * (t[i]/tau - 2*beta[i]* (tau - gamma[i]) ) | i IN range_r3 ]
437     + 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     d2DELTAbi_ddeltadtau[i] * delta * PSI[i] ) | i IN range_r4 ];
441    
442     (*--------- THERMO PROPERTY RELATIONS ----------- *)
443    
444     pressure: p
445     = rho * R * T * (1 + delta*phirdelta);
446    
447     internal_energy: u
448     = R * T * tau * (phi0tau + phirtau);
449    
450     enthalpy: h
451     = R * T * (1 + tau*(phi0tau + phirtau) + delta*phirdelta);
452    
453     entropy: s
454     = R * (tau*(phi0tau + phirtau) - phi0 - phir);
455    
456     c_isochoric: cv
457     = - R * tau^2 * (phi0tautau + phirtautau);
458    
459     c_isobaric: cp
460     = - R * (
461     tau^2 * (phi0tautau + phirtautau)
462     + ( ( 1 + delta*phirdelta - delta*tau*phirdeltatau )^2 )
463     / ( 1 + 2*delta*phirdelta + delta^2 * phirdeltadelta )
464     );
465    
466     spd_sound: w^2
467     = R * T * ( 1 + 2*delta*phirdelta + delta^2 * phirdeltadelta -
468     ( 1 + delta*phirdelta - delta*tau*phirdeltatau )^2
469     / ( tau^2 * (phi0tautau + phirtautau) )
470     );
471    
472     METHODS
473     METHOD default_self;
474     RUN ClearAll;
475     RUN specify;
476     RUN values;
477     END default_self;
478    
479     METHOD specify;
480     T.fixed := TRUE;
481     rho.fixed := TRUE;
482     END specify;
483    
484     METHOD values;
485     (* these are the test values from page 14 of the IAPWS-95 release *)
486     T := 500 {K};
487     rho := 838.025 {kg/m^3};
488     END values;
489    
490     METHOD self_test;
491     (* test the result values from page 14 of the IAPWS-95 release *)
492     ASSERT abs(phi0 - 0.204797734E1) < 1e-8;
493     ASSERT abs(phir - (-0.342693206E1)) < 1e-8;
494     ASSERT abs(phi0delta - (0.384236747)) < 1e-9;
495     ASSERT abs(phirdelta - (-0.364366650)) < 1e-9;
496     ASSERT abs(phi0deltadelta - (-0.147637878)) < 1e-9;
497     ASSERT abs(phirdeltadelta - (0.856063701)) < 1e-9;
498     ASSERT abs(phi0tau - (0.904611106e1)) < 1e-8;
499     ASSERT abs(phirtau - (-0.581403435e1)) < 1e-8;
500     ASSERT abs(phi0tautau - (-0.193249185e1)) < 1e-8;
501     ASSERT abs(phirtautau - (-0.223440737e1)) < 1e-8;
502     ASSERT abs(phi0deltatau - 0) < 1e-20;
503     ASSERT abs(phirdeltatau - (-0.112176915e1)) < 1e-8;
504     END self_test;
505    
506     END iapws95_1phase;
507    
508     MODEL iapws95_2phase REFINES thermo_state;
509     (* phase components *)
510     Sf IS_A iapws95_1phase;
511     Sg IS_A iapws95_1phase;
512    
513     (* saturation conditions *)
514     sat IS_A iapws_sat_curves;
515     Sf.T, Sg.T, sat.T ARE_THE_SAME;
516     rhof ALIASES sat.rhof;
517     rhog ALIASES sat.rhog;
518     rhov ALIASES Sg.rho;
519     rhol ALIASES Sf.rho;
520     T, Sf.T ARE_THE_SAME;
521     p, Sf.p ARE_THE_SAME;
522    
523     xf IS_A fraction;
524     xg IS_A fraction;
525     z01: xf + xg = 1.0;
526     z02: xf * (rhov-rhog) = 0.0 {kg/m^3};
527     z03: xg * (rhof-rhol) = 0.0 {kg/m^3};
528     x ALIASES xg;
529    
530     (* two-phase properties *)
531     z04: Sf.rho * Sg.rho = rho * ( xf * Sg.rho + xg * Sf.rho);
532     z05: u = xf * Sf.u + xg * Sg.u;
533     z06: h = xf * Sf.h + xg * Sg.h;
534     z07: s = xf * Sf.s + xg * Sg.s;
535     z08: cp = xf * Sf.cp + xg * Sg.cp;
536     z09: cv = xf * Sf.cv + xg * Sg.cv;
537     z10: w = xf * Sf.w + xg * Sg.w; (* check this *)
538     METHODS
539    
540     METHOD default_self;
541     RUN ClearAll;
542     RUN specify;
543     RUN values;
544     END default_self;
545    
546     METHOD specify;
547     T.fixed := TRUE;
548     x.fixed := TRUE;
549     END specify;
550    
551     METHOD values;
552     T := 390 {K};
553     x := 0.5;
554     END values;
555    
556     END iapws95_2phase;
557    

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