| 1 |
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 |
z_tau: tau * T = Tc; |
| 79 |
z_delta: delta * rhoc = rho; |
| 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 |
z_r3_b1[i]: 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 |
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 |
|
| 288 |
z_dDELTA_ddelta[i]: 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 |
z_dPSI_ddelta[i]: dPSI_ddelta[i] = -2*C[i]*d1*PSI[i]; |
| 293 |
|
| 294 |
z_dDELTAbi_ddelta[i]: dDELTAbi_ddelta[i] = b[i] * DELTA[i]^(b[i]-1) * dDELTA_ddelta[i]; |
| 295 |
|
| 296 |
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 |
|
| 299 |
z_d2PSI_ddelta2[i]: d2PSI_ddelta2[i] = (2 * C[i] * d1^2 - 1) * 2 * C[i] * PSI[i]; |
| 300 |
|
| 301 |
z_d2DELTA_ddelta2[i]: 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 |
z_d2DELTAbi_ddelta2[i]: 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 |
z_d2DELTAbi_dtau2[i]: 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 |
z_d2PSI_dtau2[i]: d2PSI_dtau2[i] = (2 * D[i] * t1^2 - 1) * 2 * D[i] * PSI[i]; |
| 315 |
|
| 316 |
z_d2PSI_ddeltadtau[i]: d2PSI_ddeltadtau[i] = 4 * C[i] * D[i] * d1 * t1 * PSI[i]; |
| 317 |
|
| 318 |
z_d2DELTAbi_ddeltadtau[i]: 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 |
z_phir_r2[i]: 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 |
z_phirdelta_r2[i]: 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 |
z_phirtau_r2[i]: 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 |
z_phirdeltadelta_r2[i]: 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 |
z_phirtautau_r2[i]: 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 |
z_phirdeltatau_r2[i]: 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 |
za_pressure: p |
| 445 |
= rho * R * T * (1 + delta*phirdelta); |
| 446 |
|
| 447 |
za_internal_energy: u |
| 448 |
= R * T * tau * (phi0tau + phirtau); |
| 449 |
|
| 450 |
za_enthalpy: h |
| 451 |
= R * T * (1 + tau*(phi0tau + phirtau) + delta*phirdelta); |
| 452 |
|
| 453 |
za_entropy: s |
| 454 |
= R * (tau*(phi0tau + phirtau) - phi0 - phir); |
| 455 |
|
| 456 |
za_c_isochoric: cv |
| 457 |
= - R * tau^2 * (phi0tautau + phirtautau); |
| 458 |
|
| 459 |
za_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 |
|
| 473 |
METHODS |
| 474 |
METHOD default_self; |
| 475 |
RUN reset; |
| 476 |
RUN values; |
| 477 |
END default_self; |
| 478 |
|
| 479 |
METHOD specify; |
| 480 |
FIX T,rho; |
| 481 |
END specify; |
| 482 |
|
| 483 |
METHOD values; |
| 484 |
(* 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 |
END values; |
| 488 |
|
| 489 |
METHOD self_test; |
| 490 |
(* 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 |
ASSERT abs(phi0delta - (0.384236747)) < 1e-9; |
| 494 |
ASSERT abs(phirdelta - (-0.364366650)) < 1e-9; |
| 495 |
ASSERT abs(phi0deltadelta - (-0.147637878)) < 1e-9; |
| 496 |
ASSERT abs(phirdeltadelta - (0.856063701)) < 1e-9; |
| 497 |
ASSERT abs(phi0tau - (0.904611106e1)) < 1e-8; |
| 498 |
ASSERT abs(phirtau - (-0.581403435e1)) < 1e-8; |
| 499 |
ASSERT abs(phi0tautau - (-0.193249185e1)) < 1e-8; |
| 500 |
ASSERT abs(phirtautau - (-0.223440737e1)) < 1e-8; |
| 501 |
ASSERT abs(phi0deltatau - 0) < 1e-20; |
| 502 |
ASSERT abs(phirdeltatau - (-0.112176915e1)) < 1e-8; |
| 503 |
END self_test; |
| 504 |
|
| 505 |
END iapws95_1phase; |
| 506 |
|
| 507 |
ATOM density_slack REFINES mass_density; |
| 508 |
lower_bound := 0.0 {kg/m^3}; |
| 509 |
END density_slack; |
| 510 |
|
| 511 |
MODEL iapws95_2phase REFINES thermo_state; |
| 512 |
(* phase components *) |
| 513 |
Sl IS_A iapws95_1phase; |
| 514 |
Sv IS_A iapws95_1phase; |
| 515 |
|
| 516 |
(* saturation conditions *) |
| 517 |
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 |
rhof ALIASES sat.rhof; |
| 523 |
rhog ALIASES sat.rhog; |
| 524 |
rhov ALIASES Sv.rho; |
| 525 |
rhol ALIASES Sl.rho; |
| 526 |
|
| 527 |
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 |
|
| 531 |
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 |
x ALIASES xg; |
| 537 |
|
| 538 |
(* 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 |
(* |
| 546 |
z10: w = xf * Sl.w + xg * Sv.w; * check this * |
| 547 |
*) |
| 548 |
METHODS |
| 549 |
|
| 550 |
METHOD default_self; |
| 551 |
RUN reset; |
| 552 |
RUN values; |
| 553 |
Sl.rho := 800 {kg/m^3}; |
| 554 |
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 |
END default_self; |
| 561 |
|
| 562 |
METHOD specify; |
| 563 |
FIX T,x; |
| 564 |
END specify; |
| 565 |
|
| 566 |
METHOD values; |
| 567 |
T := 390 {K}; |
| 568 |
x := 0.5; |
| 569 |
END values; |
| 570 |
|
| 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 |
|
| 580 |
Sl.d1.lower_bound := -1; |
| 581 |
Sv.d1.lower_bound := -1; |
| 582 |
END bound_self; |
| 583 |
|
| 584 |
END iapws95_2phase; |
| 585 |
|
| 586 |
(* |
| 587 |
the following models are me trying to work out a robust way of calculating |
| 588 |
properties in terms of (p,h). |
| 589 |
|
| 590 |
Art suggested that I must first solve for a related point in terms (T,rho) |
| 591 |
(the 'natural parameters' then reorganise and solve for (p,h) at the |
| 592 |
point I need. |
| 593 |
*) |
| 594 |
|
| 595 |
(* |
| 596 |
Superheaded region in terms of (p,h). |
| 597 |
First load and run 'default_self' then solve. Then run 'default_ph' and |
| 598 |
try to solve again -- fails. |
| 599 |
*) |
| 600 |
MODEL iapws95_2phase_ph_test REFINES iapws95_2phase; |
| 601 |
|
| 602 |
METHODS |
| 603 |
METHOD specify; |
| 604 |
FIX T, rho; |
| 605 |
END specify; |
| 606 |
|
| 607 |
METHOD values; |
| 608 |
T := 452 {K}; |
| 609 |
x := 1; |
| 610 |
rho := 4 {kg/m^3}; |
| 611 |
END values; |
| 612 |
|
| 613 |
METHOD specify_ph; |
| 614 |
FIX p,h; |
| 615 |
END specify_ph; |
| 616 |
METHOD values_ph; |
| 617 |
p := 10 {bar}; |
| 618 |
h := 2998 {kJ/kg}; |
| 619 |
END values_ph; |
| 620 |
|
| 621 |
METHOD default_ph; |
| 622 |
RUN ClearAll; RUN specify_ph; RUN values_ph; |
| 623 |
END default_ph; |
| 624 |
|
| 625 |
METHOD self_test; |
| 626 |
ASSERT abs(T - 548.15 {K}) < 1 {K}; |
| 627 |
ASSERT abs(rho - 8.43 {kg/m^3}) < 0.2 {kg/m^3}; |
| 628 |
END self_test; |
| 629 |
|
| 630 |
END iapws95_2phase_ph_test; |
| 631 |
|
| 632 |
(* |
| 633 |
Subcooled. This one works: run the 'default_self' then solve. Then |
| 634 |
run 'default_ph' and solve: it converges fine. |
| 635 |
*) |
| 636 |
MODEL iapws95_2phase_ph_test2 REFINES iapws95_2phase; |
| 637 |
|
| 638 |
METHODS |
| 639 |
METHOD specify; |
| 640 |
FIX T,rho; |
| 641 |
END specify; |
| 642 |
METHOD values; |
| 643 |
T := 290 {K}; |
| 644 |
rho := 999 {kg/m^3}; |
| 645 |
END values; |
| 646 |
METHOD default_ph; |
| 647 |
RUN ClearAll; RUN specify_ph; RUN values_ph; |
| 648 |
END default_ph; |
| 649 |
METHOD specify_ph; |
| 650 |
FIX p,h; |
| 651 |
END specify_ph; |
| 652 |
METHOD values_ph; |
| 653 |
p := 10 {bar}; |
| 654 |
h := 418.7 {kJ/kg}; |
| 655 |
END values_ph; |
| 656 |
METHOD self_test; |
| 657 |
ASSERT abs(T - 373.15 {K}) < 1 {K}; |
| 658 |
ASSERT abs(rho - 959 {kg/m^3}) < 3 {kg/m^3}; |
| 659 |
END self_test; |
| 660 |
|
| 661 |
END iapws95_2phase_ph_test2; |
| 662 |
|
| 663 |
(* |
| 664 |
This last model *will* converge with careful coaxing |
| 665 |
You need to alternate between ROW2NORM scaling |
| 666 |
and RELNOM scaling using the QRSlv engine. |
| 667 |
Some tweaking of delta/d1 values may be necessary as |
| 668 |
well. |
| 669 |
*) |
| 670 |
MODEL iapws95_2phase_ph_test3 REFINES iapws95_1phase; |
| 671 |
|
| 672 |
METHODS |
| 673 |
METHOD specify; |
| 674 |
FIX p,h; |
| 675 |
END specify; |
| 676 |
METHOD values; |
| 677 |
p := 1 {bar}; |
| 678 |
h := 104.9 {kJ/kg}; |
| 679 |
(* free vars *) |
| 680 |
delta := 3; |
| 681 |
d1 := 2; |
| 682 |
rho := 900; |
| 683 |
END values; |
| 684 |
METHOD self_test; |
| 685 |
ASSERT abs(T - 298.15 {K}) < 1 {K}; |
| 686 |
ASSERT abs(rho - 997 {kg/m^3}) < 3 {kg/m^3}; |
| 687 |
END self_test; |
| 688 |
|
| 689 |
END iapws95_2phase_ph_test3; |
Parent Directory
| 
Revision Log