www.mooseframework.org
CO2FluidProperties.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
8 #include "CO2FluidProperties.h"
9 #include "BrentsMethod.h"
10 
11 template <>
12 InputParameters
14 {
15  InputParameters params = validParams<SinglePhaseFluidPropertiesPT>();
16  params.addClassDescription(
17  "Fluid properties for carbon dioxide (CO2) using the Span & Wagner EOS");
18  return params;
19 }
20 
21 CO2FluidProperties::CO2FluidProperties(const InputParameters & parameters)
22  : SinglePhaseFluidPropertiesPT(parameters)
23 {
24 }
25 
27 
28 std::string
30 {
31  return "co2";
32 }
33 
34 Real
36 {
37  return _Mco2;
38 }
39 
40 Real
42 {
43  return _critical_pressure;
44 }
45 
46 Real
48 {
49  return _critical_temperature;
50 }
51 
52 Real
54 {
55  return _critical_density;
56 }
57 
58 Real
60 {
62 }
63 
64 Real
66 {
68 }
69 
70 Real
72 {
73  if (temperature < _triple_point_temperature)
74  mooseError(
75  "Temperature is below the triple point temperature in ", name(), ": meltingPressure()");
76 
77  Real Tstar = temperature / _triple_point_temperature;
78 
79  return _triple_point_pressure *
80  (1.0 + 1955.539 * (Tstar - 1.0) + 2055.4593 * std::pow(Tstar - 1.0, 2.0));
81 }
82 
83 Real
85 {
86  if (temperature > _triple_point_temperature)
87  mooseError(
88  "Temperature is above the triple point temperature in ", name(), ": sublimationPressure()");
89 
90  Real Tstar = temperature / _triple_point_temperature;
91 
93  std::exp((-14.740846 * (1.0 - Tstar) + 2.4327015 * std::pow(1.0 - Tstar, 1.9) -
94  5.3061778 * std::pow(1.0 - Tstar, 2.9)) /
95  Tstar);
96 
97  return pressure;
98 }
99 
100 Real
102 {
103  if (temperature < _triple_point_temperature || temperature > _critical_temperature)
104  mooseError("Temperature is out of range in ", name(), ": vaporPressure()");
105 
106  Real Tstar = temperature / _critical_temperature;
107 
108  Real logpressure =
109  (-7.0602087 * (1.0 - Tstar) + 1.9391218 * std::pow(1.0 - Tstar, 1.5) -
110  1.6463597 * std::pow(1.0 - Tstar, 2.0) - 3.2995634 * std::pow(1.0 - Tstar, 4.0)) /
111  Tstar;
112 
113  return _critical_pressure * std::exp(logpressure);
114 }
115 
116 Real
118 {
119  if (temperature < _triple_point_temperature || temperature > _critical_temperature)
120  mooseError("Temperature is out of range in ", name(), ": saturatedLiquiDensity()");
121 
122  Real Tstar = temperature / _critical_temperature;
123 
124  Real logdensity = 1.9245108 * std::pow(1.0 - Tstar, 0.34) -
125  0.62385555 * std::pow(1.0 - Tstar, 0.5) -
126  0.32731127 * std::pow(1.0 - Tstar, 10.0 / 6.0) +
127  0.39245142 * std::pow(1.0 - Tstar, 11.0 / 6.0);
128 
129  return _critical_density * std::exp(logdensity);
130 }
131 
132 Real
134 {
135  if (temperature < _triple_point_temperature || temperature > _critical_temperature)
136  mooseError("Temperature is out of range in ", name(), ": saturatedVaporDensity()");
137 
138  Real Tstar = temperature / _critical_temperature;
139 
140  Real logdensity =
141  (-1.7074879 * std::pow(1.0 - Tstar, 0.34) - 0.82274670 * std::pow(1.0 - Tstar, 0.5) -
142  4.6008549 * (1.0 - Tstar) - 10.111178 * std::pow(1.0 - Tstar, 7.0 / 3.0) -
143  29.742252 * std::pow(1.0 - Tstar, 14.0 / 3.0));
144 
145  return _critical_density * std::exp(logdensity);
146 }
147 
148 Real
149 CO2FluidProperties::phiSW(Real delta, Real tau) const
150 {
151  // Ideal gas component of the Helmholtz free energy
152  Real sum0 = 0.0;
153  for (std::size_t i = 0; i < _a0.size(); ++i)
154  sum0 += _a0[i] * std::log(1.0 - std::exp(-_theta0[i] * tau));
155 
156  Real phi0 = std::log(delta) + 8.37304456 - 3.70454304 * tau + 2.5 * std::log(tau) + sum0;
157 
158  // Residual component of the Helmholtz free energy
159  Real theta, Delta, Psi;
160  Real phir = 0.0;
161  for (std::size_t i = 0; i < _n1.size(); ++i)
162  phir += _n1[i] * std::pow(delta, _d1[i]) * std::pow(tau, _t1[i]);
163 
164  for (std::size_t i = 0; i < _n2.size(); ++i)
165  phir += _n2[i] * std::pow(delta, _d2[i]) * std::pow(tau, _t2[i]) *
166  std::exp(-std::pow(delta, _c2[i]));
167 
168  for (std::size_t i = 0; i < _n3.size(); ++i)
169  phir += _n3[i] * std::pow(delta, _d3[i]) * std::pow(tau, _t3[i]) *
170  std::exp(-_alpha3[i] * std::pow(delta - _eps3[i], 2.0) -
171  _beta3[i] * std::pow(tau - _gamma3[i], 2.0));
172 
173  for (std::size_t i = 0; i < _n4.size(); ++i)
174  {
175  theta = 1.0 - tau + _A4[i] * std::pow(std::pow(delta - 1.0, 2.0), 1.0 / (2.0 * _beta4[i]));
176  Delta = std::pow(theta, 2) + _B4[i] * std::pow(std::pow(delta - 1.0, 2.0), _a4[i]);
177  Psi = std::exp(-_C4[i] * std::pow(delta - 1.0, 2.0) - _D4[i] * std::pow(tau - 1.0, 2.0));
178  phir += _n4[i] * std::pow(Delta, _b4[i]) * delta * Psi;
179  }
180 
181  // The Helmholtz free energy is the sum of these components
182  return phi0 + phir;
183 }
184 
185 Real
186 CO2FluidProperties::dphiSW_dd(Real delta, Real tau) const
187 {
188  // Derivative of the ideal gas component wrt gamma
189  Real dphi0dd = 1.0 / delta;
190 
191  // Derivative of the residual component wrt gamma
192  Real theta, Delta, Psi, dDelta_dd, dPsi_dd;
193  Real dphirdd = 0.0;
194 
195  for (std::size_t i = 0; i < _n1.size(); ++i)
196  dphirdd += _n1[i] * _d1[i] * std::pow(delta, _d1[i] - 1.0) * std::pow(tau, _t1[i]);
197 
198  for (std::size_t i = 0; i < _n2.size(); ++i)
199  dphirdd += _n2[i] * std::exp(-std::pow(delta, _c2[i])) *
200  (std::pow(delta, _d2[i] - 1.0) * std::pow(tau, _t2[i]) *
201  (_d2[i] - _c2[i] * std::pow(delta, _c2[i])));
202 
203  for (std::size_t i = 0; i < _n3.size(); ++i)
204  dphirdd += _n3[i] * std::pow(delta, _d3[i]) * std::pow(tau, _t3[i]) *
205  std::exp(-_alpha3[i] * std::pow(delta - _eps3[i], 2.0) -
206  _beta3[i] * std::pow(tau - _gamma3[i], 2.0)) *
207  (_d3[i] / delta - 2.0 * _alpha3[i] * (delta - _eps3[i]));
208 
209  for (std::size_t i = 0; i < _n4.size(); ++i)
210  {
211  theta = 1.0 - tau + _A4[i] * std::pow(std::pow(delta - 1.0, 2.0), 1.0 / (2.0 * _beta4[i]));
212  Delta = std::pow(theta, 2.0) + _B4[i] * std::pow(std::pow(delta - 1.0, 2.0), _a4[i]);
213  Psi = std::exp(-_C4[i] * std::pow(delta - 1.0, 2.0) - _D4[i] * std::pow(tau - 1.0, 2.0));
214  dPsi_dd = -2.0 * _C4[i] * (delta - 1.0) * Psi;
215  dDelta_dd = (delta - 1.0) *
216  (_A4[i] * theta * 2.0 / _beta4[i] *
217  std::pow(std::pow(delta - 1.0, 2.0), 1.0 / (2.0 * _beta4[i]) - 1.0) +
218  2.0 * _B4[i] * _a4[i] * std::pow(std::pow(delta - 1.0, 2.0), _a4[i] - 1.0));
219 
220  dphirdd += _n4[i] * (std::pow(Delta, _b4[i]) * (Psi + delta * dPsi_dd) +
221  _b4[i] * std::pow(Delta, _b4[i] - 1.0) * dDelta_dd * delta * Psi);
222  }
223 
224  // The derivative of the free energy wrt delta is the sum of these components
225  return dphi0dd + dphirdd;
226 }
227 
228 Real
229 CO2FluidProperties::dphiSW_dt(Real delta, Real tau) const
230 {
231  // Derivative of the ideal gas component wrt tau
232  Real sum0 = 0.0;
233  for (std::size_t i = 0; i < _a0.size(); ++i)
234  sum0 += _a0[i] * _theta0[i] * (1.0 / (1.0 - std::exp(-_theta0[i] * tau)) - 1.0);
235 
236  Real dphi0dt = -3.70454304 + 2.5 / tau + sum0;
237 
238  // Derivative of the residual component wrt tau
239  Real theta, Delta, Psi, dDelta_dt, dPsi_dt;
240  Real dphirdt = 0.0;
241  for (std::size_t i = 0; i < _n1.size(); ++i)
242  dphirdt += _n1[i] * _t1[i] * std::pow(delta, _d1[i]) * std::pow(tau, _t1[i] - 1.0);
243 
244  for (std::size_t i = 0; i < _n2.size(); ++i)
245  dphirdt += _n2[i] * _t2[i] * std::pow(delta, _d2[i]) * std::pow(tau, _t2[i] - 1.0) *
246  std::exp(-std::pow(delta, _c2[i]));
247 
248  for (std::size_t i = 0; i < _n3.size(); ++i)
249  dphirdt += _n3[i] * std::pow(delta, _d3[i]) * std::pow(tau, _t3[i]) *
250  std::exp(-_alpha3[i] * std::pow(delta - _eps3[i], 2.0) -
251  _beta3[i] * std::pow(tau - _gamma3[i], 2.0)) *
252  (_t3[i] / tau - 2.0 * _beta3[i] * (tau - _gamma3[i]));
253 
254  for (std::size_t i = 0; i < _n4.size(); ++i)
255  {
256  theta = 1.0 - tau + _A4[i] * std::pow(std::pow(delta - 1.0, 2.0), 1.0 / (2.0 * _beta4[i]));
257  Delta = std::pow(theta, 2.0) + _B4[i] * std::pow(std::pow(delta - 1.0, 2.0), _a4[i]);
258  Psi = std::exp(-_C4[i] * std::pow(delta - 1.0, 2.0) - _D4[i] * std::pow(tau - 1.0, 2.0));
259  dDelta_dt = -2.0 * theta * _b4[i] * std::pow(Delta, _b4[i] - 1.0);
260  dPsi_dt = -2.0 * _D4[i] * (tau - 1.0) * Psi;
261 
262  dphirdt += _n4[i] * delta * (Psi * dDelta_dt + std::pow(Delta, _b4[i]) * dPsi_dt);
263  }
264 
265  // The derivative of the free energy wrt tau is the sum of these components
266  return dphi0dt + dphirdt;
267 }
268 
269 Real
270 CO2FluidProperties::d2phiSW_dd2(Real delta, Real tau) const
271 {
272  // Second derivative of the ideal gas component wrt gamma
273  Real d2phi0dd2 = -1.0 / delta / delta;
274 
275  // Second derivative of the residual component wrt gamma
276  Real d2phirdd2 = 0.0;
277  Real theta, Delta, Psi, dDelta_dd, dPsi_dd, d2Delta_dd2, d2Psi_dd2;
278 
279  for (std::size_t i = 0; i < _n1.size(); ++i)
280  d2phirdd2 +=
281  _n1[i] * _d1[i] * (_d1[i] - 1.0) * std::pow(delta, _d1[i] - 2.0) * std::pow(tau, _t1[i]);
282 
283  for (std::size_t i = 0; i < _n2.size(); ++i)
284  d2phirdd2 += _n2[i] * std::exp(-std::pow(delta, _c2[i])) * std::pow(delta, _d2[i] - 2.0) *
285  std::pow(tau, _t2[i]) * ((_d2[i] - _c2[i] * std::pow(delta, _c2[i])) *
286  (_d2[i] - 1.0 - _c2[i] * std::pow(delta, _c2[i])) -
287  _c2[i] * _c2[i] * std::pow(delta, _c2[i]));
288 
289  for (std::size_t i = 0; i < _n3.size(); ++i)
290  d2phirdd2 +=
291  _n3[i] * std::pow(tau, _t3[i]) * std::exp(-_alpha3[i] * std::pow(delta - _eps3[i], 2.0) -
292  _beta3[i] * std::pow(tau - _gamma3[i], 2.0)) *
293  (-2.0 * _alpha3[i] * std::pow(delta, _d3[i]) +
294  4.0 * _alpha3[i] * _alpha3[i] * std::pow(delta, _d3[i]) * std::pow(delta - _eps3[i], 2.0) -
295  4.0 * _d3[i] * _alpha3[i] * std::pow(delta, _d3[i] - 1.0) * (delta - _eps3[i]) +
296  _d3[i] * (_d3[i] - 1.0) * std::pow(delta, _d3[i] - 2.0));
297 
298  for (std::size_t i = 0; i < _n4.size(); ++i)
299  {
300  theta = 1.0 - tau + _A4[i] * std::pow(std::pow(delta - 1.0, 2.0), 1.0 / (2.0 * _beta4[i]));
301  Delta = std::pow(theta, 2.0) + _B4[i] * std::pow(std::pow(delta - 1.0, 2.0), _a4[i]);
302  Psi = std::exp(-_C4[i] * std::pow(delta - 1.0, 2.0) - _D4[i] * std::pow(tau - 1.0, 2.0));
303  dPsi_dd = -2.0 * _C4[i] * (delta - 1.0) * Psi;
304  dDelta_dd = (delta - 1.0) *
305  (_A4[i] * theta * 2.0 / _beta4[i] *
306  std::pow(std::pow(delta - 1.0, 2.0), 1.0 / (2.0 * _beta4[i]) - 1.0) +
307  2.0 * _B4[i] * _a4[i] * std::pow(std::pow(delta - 1.0, 2.0), _a4[i] - 1.0));
308  d2Psi_dd2 = 3.0 * _D4[i] * Psi * (2.0 * _C4[i] * std::pow(delta - 1.0, 2.0) - 1.0);
309  d2Delta_dd2 =
310  1.0 / (delta - 1.0) * dDelta_dd +
311  (delta - 1.0) * (delta - 1.0) *
312  (4.0 * _B4[i] * _a4[i] * (_a4[i] - 1.0) *
313  std::pow(std::pow(delta - 1.0, 2.0), _a4[i] - 2.0) +
314  2.0 * _A4[i] * _A4[i] *
315  std::pow(std::pow(std::pow(delta - 1.0, 2.0), 1.0 / (2.0 * _beta4[i]) - 1.0),
316  2.0) /
317  _beta4[i] / _beta4[i] +
318  (4.0 / _beta4[i]) * _A4[i] * theta * (1.0 / (2.0 * _beta4[i]) - 1.0) *
319  std::pow(std::pow(delta - 1.0, 2.0), 1.0 / (2.0 * _beta4[i]) - 2.0));
320  d2phirdd2 +=
321  _n4[i] *
322  (std::pow(Delta, _b4[i]) * (2.0 * dPsi_dd + delta * d2Psi_dd2) +
323  2.0 * _b4[i] * std::pow(Delta, _b4[i] - 1.0) * dDelta_dd * (Psi + delta * dPsi_dd) +
324  _b4[i] * (std::pow(Delta, _b4[i] - 1.0) * d2Delta_dd2 +
325  (_b4[i] - 1.0) * std::pow(Delta, _b4[i] - 2.0) * std::pow(dDelta_dd, 2.0)) *
326  delta * Psi);
327  }
328  // The second derivative of the free energy wrt delta is the sum of these components
329  return d2phi0dd2 + d2phirdd2;
330 }
331 
332 Real
333 CO2FluidProperties::d2phiSW_dt2(Real delta, Real tau) const
334 {
335  // Second derivative of the ideal gas component wrt tau
336  Real sum0 = 0.0;
337  for (std::size_t i = 0; i < _a0.size(); ++i)
338  sum0 += _a0[i] * _theta0[i] * _theta0[i] * std::exp(-_theta0[i] * tau) *
339  std::pow(1.0 - std::exp(-_theta0[i] * tau), -2.0);
340 
341  Real d2phi0dt2 = -2.5 / tau / tau - sum0;
342 
343  // Second derivative of the residual component wrt tau
344  Real d2phirdt2 = 0.0;
345  Real theta, Delta, Psi, dPsi_dt, dDelta_dt, d2Delta_dt2, d2Psi_dt2;
346 
347  for (std::size_t i = 0; i < _n1.size(); ++i)
348  d2phirdt2 +=
349  _n1[i] * _t1[i] * (_t1[i] - 1.0) * std::pow(delta, _d1[i]) * std::pow(tau, _t1[i] - 2.0);
350 
351  for (std::size_t i = 0; i < _n2.size(); ++i)
352  d2phirdt2 += _n2[i] * _t2[i] * (_t2[i] - 1.0) * std::pow(delta, _d2[i]) *
353  std::exp(-std::pow(delta, _c2[i])) * std::pow(tau, _t2[i] - 2.0);
354 
355  for (std::size_t i = 0; i < _n3.size(); ++i)
356  d2phirdt2 += _n3[i] * std::pow(delta, _d3[i]) * std::pow(tau, _t3[i]) *
357  std::exp(-_alpha3[i] * std::pow(delta - _eps3[i], 2.0) -
358  _beta3[i] * std::pow(tau - _gamma3[i], 2.0)) *
359  (std::pow(_t3[i] / tau - 2.0 * _beta3[i] * (tau - _gamma3[i]), 2.0) -
360  _t3[i] / tau / tau - 2.0 * _beta3[i]);
361 
362  for (std::size_t i = 0; i < _n4.size(); ++i)
363  {
364  theta = 1.0 - tau + _A4[i] * std::pow(std::pow(delta - 1.0, 2.0), 1.0 / (2.0 * _beta4[i]));
365  Delta = std::pow(theta, 2.0) + _B4[i] * std::pow(std::pow(delta - 1.0, 2.0), _a4[i]);
366  Psi = std::exp(-_C4[i] * std::pow(delta - 1.0, 2.0) - _D4[i] * std::pow(tau - 1.0, 2.0));
367  dDelta_dt = -2.0 * theta * _b4[i] * std::pow(Delta, _b4[i] - 1.0);
368  d2Delta_dt2 = 2.0 * _b4[i] * std::pow(Delta, _b4[i] - 1.0) +
369  4.0 * theta * theta * _b4[i] * (_b4[i] - 1.0) * std::pow(Delta, _b4[i] - 2.0);
370  dPsi_dt = -2.0 * _D4[i] * (tau - 1.0) * Psi;
371  d2Psi_dt2 = 2.0 * _D4[i] * (2.0 * _D4[i] * (tau - 1.0) * (tau - 1.0) - 1.0) * Psi;
372  d2phirdt2 += _n4[i] * delta * (Psi * d2Delta_dt2 + 2.0 * dDelta_dt * dPsi_dt +
373  std::pow(Delta, _b4[i]) * d2Psi_dt2);
374  }
375 
376  // The second derivative of the free energy wrt tau is the sum of these components
377  return d2phi0dt2 + d2phirdt2;
378 }
379 
380 Real
381 CO2FluidProperties::d2phiSW_ddt(Real delta, Real tau) const
382 {
383  // Note: second derivative of the ideal gas component wrt delta and tau is 0
384  // Derivative of the residual component wrt gamma
385  Real theta, Delta, Psi, dDelta_dd, dPsi_dd, dDelta_dt, dPsi_dt, d2Delta_ddt, d2Psi_ddt;
386  Real d2phirddt = 0.0;
387  for (std::size_t i = 0; i < _n1.size(); ++i)
388  d2phirddt +=
389  _n1[i] * _d1[i] * _t1[i] * std::pow(delta, _d1[i] - 1.0) * std::pow(tau, _t1[i] - 1.0);
390 
391  for (std::size_t i = 0; i < _n2.size(); ++i)
392  d2phirddt += _n2[i] * std::exp(-std::pow(delta, _c2[i])) *
393  (std::pow(delta, _d2[i] - 1.0) * _t2[i] * std::pow(tau, _t2[i] - 1.0) *
394  (_d2[i] - _c2[i] * std::pow(delta, _c2[i])));
395 
396  for (std::size_t i = 0; i < _n3.size(); ++i)
397  d2phirddt += _n3[i] * std::pow(delta, _d3[i]) * std::pow(tau, _t3[i]) *
398  std::exp(-_alpha3[i] * std::pow(delta - _eps3[i], 2.0) -
399  _beta3[i] * std::pow(tau - _gamma3[i], 2.0)) *
400  (_d3[i] / delta - 2.0 * _alpha3[i] * (delta - _eps3[i])) *
401  (_t3[i] / tau - 2.0 * _beta3[i] * (tau - _gamma3[i]));
402 
403  for (std::size_t i = 0; i < _n4.size(); ++i)
404  {
405  theta = 1.0 - tau + _A4[i] * std::pow(std::pow(delta - 1.0, 2.0), 1.0 / (2.0 * _beta4[i]));
406  Delta = std::pow(theta, 2.0) + _B4[i] * std::pow(std::pow(delta - 1.0, 2.0), _a4[i]);
407  Psi = std::exp(-_C4[i] * std::pow(delta - 1.0, 2.0) - _D4[i] * std::pow(tau - 1.0, 2.0));
408  dPsi_dd = -2.0 * _C4[i] * (delta - 1.0) * Psi;
409  dPsi_dt = -2.0 * _D4[i] * (tau - 1.0) * Psi;
410  d2Psi_ddt = 4.0 * _C4[i] * _D4[i] * (delta - 1.0) * (tau - 1.0) * Psi;
411  dDelta_dd = (delta - 1.0) *
412  (_A4[i] * theta * 2.0 / _beta4[i] *
413  std::pow(std::pow(delta - 1.0, 2.0), 1.0 / (2.0 * _beta4[i]) - 1.0) +
414  2.0 * _B4[i] * _a4[i] * std::pow(std::pow(delta - 1.0, 2.0), _a4[i] - 1.0));
415  dDelta_dt = -2.0 * theta * _b4[i] * std::pow(Delta, _b4[i] - 1.0);
416  d2Delta_ddt = -2.0 * _A4[i] * _b4[i] / _beta4[i] * std::pow(Delta, _b4[i] - 1.0) *
417  (delta - 1.0) *
418  std::pow(std::pow(delta - 1.0, 2.0), 1.0 / (2.0 * _beta4[i]) - 1.0) -
419  2.0 * theta * _b4[i] * (_b4[i] - 1.0) * std::pow(Delta, _b4[i] - 2.0) * dDelta_dd;
420 
421  d2phirddt += _n4[i] * (std::pow(Delta, _b4[i]) * (dPsi_dt + delta * d2Psi_ddt) +
422  delta * _b4[i] * std::pow(Delta, _b4[i] - 1.0) * dDelta_dd * dPsi_dt +
423  dDelta_dt * (Psi + delta * dPsi_dd) + d2Delta_ddt * delta * Psi);
424  }
425 
426  return d2phirddt;
427 }
428 
429 Real
431 {
432  // Scale the input density and temperature
433  Real delta = density / _critical_density;
434  Real tau = _critical_temperature / temperature;
435 
436  return _Rco2 * temperature * density * delta * dphiSW_dd(delta, tau);
437 }
438 
439 Real
441 {
442  // Check that the input parameters are within the region of validity
443  if (temperature < 216.0 || temperature > 1100.0 || density <= 0.0)
444  mooseError("Parameters out of range in ", name(), ": pressure()");
445 
446  Real pressure = 0.0;
447 
448  if (temperature > _triple_point_temperature && temperature < _critical_temperature)
449  {
450  Real gas_density = saturatedVaporDensity(temperature);
451  Real liquid_density = saturatedLiquidDensity(temperature);
452 
453  if (density < gas_density || density > liquid_density)
454  pressure = pressureSW(density, temperature);
455  else
456  pressure = vaporPressure(temperature);
457  }
458  else
459  pressure = pressureSW(density, temperature);
460 
461  return pressure;
462 }
463 
464 Real
466 {
467  // Check that the input parameters are within the region of validity
468  if (temperature < 216.0 || temperature > 1100.0 || pressure <= 0.0)
469  mooseError("Parameters out of range in ", name(), ": rho()");
470 
471  // Also check that the pressure and temperature are not in the solid phase region
472  if (((temperature > _triple_point_temperature) && (pressure > meltingPressure(temperature))) ||
473  ((temperature < _triple_point_temperature) && (pressure > sublimationPressure(temperature))))
474  mooseError(
475  "Input pressure and temperature in ", name(), ": rho() correspond to solid CO2 phase");
476 
477  Real density;
478  // Initial estimate of a bracketing interval for the density
479  Real lower_density = 100.0;
480  Real upper_density = 1000.0;
481 
482  // The density is found by finding the zero of the pressure calculated using the
483  // Span and Wagner EOS minus the input pressure
484  auto pressure_diff = [&pressure, &temperature, this](Real x) {
485  return this->pressure(x, temperature) - pressure;
486  };
487 
488  BrentsMethod::bracket(pressure_diff, lower_density, upper_density);
489  density = BrentsMethod::root(pressure_diff, lower_density, upper_density);
490 
491  return density;
492 }
493 
494 void
496  Real pressure, Real temperature, Real & rho, Real & drho_dp, Real & drho_dT) const
497 {
498  rho = this->rho(pressure, temperature);
499  // Scale the density and temperature
500  Real delta = rho / _critical_density;
501  Real tau = _critical_temperature / temperature;
502  Real dpdd = dphiSW_dd(delta, tau);
503  Real d2pdd2 = d2phiSW_dd2(delta, tau);
504 
505  drho_dp = 1.0 / (_Rco2 * temperature * delta * (2.0 * dpdd + delta * d2pdd2));
506  drho_dT =
507  rho * (tau * d2phiSW_ddt(delta, tau) - dpdd) / temperature / (2.0 * dpdd + delta * d2pdd2);
508 }
509 
510 Real
512 {
513  Real rho = this->rho(pressure, temperature);
514  return this->mu_from_rho_T(rho, temperature);
515 }
516 
517 void
519  Real pressure, Real temperature, Real & mu, Real & dmu_dp, Real & dmu_dT) const
520 {
521  Real rho, drho_dp, drho_dT;
522  this->rho_dpT(pressure, temperature, rho, drho_dp, drho_dT);
523 
524  Real dmu_drho;
525  mu_drhoT_from_rho_T(rho, temperature, drho_dT, mu, dmu_drho, dmu_dT);
526  dmu_dp = dmu_drho * drho_dp;
527 }
528 
529 Real
531 {
532  // Check that the input parameters are within the region of validity
533  if (temperature < 216.0 || temperature > 1000.0 || density > 1400.0)
534  mooseError("Parameters out of range in ", name(), ": mu()");
535 
536  Real Tstar = temperature / 251.196;
537  const std::vector<Real> a{0.235156, -0.491266, 5.211155e-2, 5.347906e-2, -1.537102e-2};
538  const std::vector<Real> d{
539  0.4071119e-2, 0.7198037e-4, 0.2411697e-16, 0.2971072e-22, -0.1627888e-22};
540 
541  // Viscosity in the zero-density limit
542  Real sum = 0.0;
543 
544  for (std::size_t i = 0; i < a.size(); ++i)
545  sum += a[i] * std::pow(std::log(Tstar), i);
546 
547  Real mu0 = 1.00697 * std::sqrt(temperature) / std::exp(sum);
548 
549  // Excess viscosity due to finite density
550  Real mue = d[0] * density + d[1] * std::pow(density, 2) +
551  d[2] * std::pow(density, 6) / std::pow(Tstar, 3) + d[3] * std::pow(density, 8) +
552  d[4] * std::pow(density, 8) / Tstar;
553 
554  return (mu0 + mue) * 1.0e-6; // convert to Pa.s
555 }
556 
557 void
559  Real temperature,
560  Real ddensity_dT,
561  Real & mu,
562  Real & dmu_drho,
563  Real & dmu_dT) const
564 {
565  // Check that the input parameters are within the region of validity
566  if (temperature < 216.0 || temperature > 1000.0 || density > 1400.0)
567  mooseError("Parameters out of range in ", name(), ": mu_drhoT()");
568 
569  Real Tstar = temperature / 251.196;
570  Real dTstar_dT = 1.0 / 251.196;
571  const std::vector<Real> a{0.235156, -0.491266, 5.211155e-2, 5.347906e-2, -1.537102e-2};
572  const std::vector<Real> d{
573  0.4071119e-2, 0.7198037e-4, 0.2411697e-16, 0.2971072e-22, -0.1627888e-22};
574 
575  // Viscosity in the zero-density limit. Note this is only a function of T.
576  // Start the sum at i = 1 so the derivative is defined
577  Real sum0 = a[0], dsum0_dTstar = 0.0;
578 
579  for (std::size_t i = 1; i < a.size(); ++i)
580  {
581  sum0 += a[i] * std::pow(std::log(Tstar), i);
582  dsum0_dTstar += i * a[i] * std::pow(std::log(Tstar), i - 1) / Tstar;
583  }
584 
585  Real mu0 = 1.00697 * std::sqrt(temperature) / std::exp(sum0);
586  Real dmu0_dT = (0.5 * 1.00697 / std::sqrt(temperature) -
587  1.00697 * std::sqrt(temperature) * dsum0_dTstar * dTstar_dT) /
588  std::exp(sum0);
589 
590  // Excess viscosity due to finite density
591  Real mue = d[0] * density + d[1] * std::pow(density, 2) +
592  d[2] * std::pow(density, 6) / std::pow(Tstar, 3) + d[3] * std::pow(density, 8) +
593  d[4] * std::pow(density, 8) / Tstar;
594 
595  Real dmue_drho = d[0] + 2.0 * d[1] * density +
596  6.0 * d[2] * std::pow(density, 5) / std::pow(Tstar, 3) +
597  8.0 * d[3] * std::pow(density, 7) + 8.0 * d[4] * std::pow(density, 7) / Tstar;
598 
599  Real dmue_dT = (-3.0 * d[2] * std::pow(density, 6) / std::pow(Tstar, 4) -
600  d[4] * std::pow(density, 8) / Tstar / Tstar) *
601  dTstar_dT;
602 
603  // Viscosity in Pa.s is
604  mu = (mu0 + mue) * 1.0e-6;
605  dmu_drho = dmue_drho * 1.0e-6;
606  dmu_dT = (dmu0_dT + dmue_dT) * 1.0e-6 + dmu_drho * ddensity_dT;
607 }
608 
609 Real
611 {
612  // This correlation uses temperature in C
613  Real Tc = temperature - _T_c2k;
614  // The parial volume
615  Real V = 37.51 - 9.585e-2 * Tc + 8.74e-4 * Tc * Tc - 5.044e-7 * Tc * Tc * Tc;
616 
617  return 1.0e6 * _Mco2 / V;
618 }
619 
620 Real
622 {
623  return henryConstantIAPWS(temperature, -8.55445, 4.01195, 9.52345);
624 }
625 
626 void
627 CO2FluidProperties::henryConstant_dT(Real temperature, Real & Kh, Real & dKh_dT) const
628 {
629  henryConstantIAPWS_dT(temperature, Kh, dKh_dT, -8.55445, 4.01195, 9.52345);
630 }
631 
632 Real
634 {
635  // Require density first
636  Real density = rho(pressure, temperature);
637  // Scale the input density and temperature
638  Real delta = density / _critical_density;
639  Real tau = _critical_temperature / temperature;
640 
641  return _Rco2 * temperature * tau * dphiSW_dt(delta, tau);
642 }
643 
644 void
646  Real pressure, Real temperature, Real & e, Real & de_dp, Real & de_dT) const
647 {
648  // Require density first
649  Real density = rho(pressure, temperature);
650  // Scale the input density and temperature
651  Real delta = density / _critical_density;
652  Real tau = _critical_temperature / temperature;
653 
654  e = _Rco2 * temperature * tau * dphiSW_dt(delta, tau);
655  Real dpdd = dphiSW_dd(delta, tau);
656  Real d2pdd2 = d2phiSW_dd2(delta, tau);
657  Real d2pddt = d2phiSW_ddt(delta, tau);
658 
659  de_dp = tau * d2pddt / (density * (2.0 * dpdd + delta * d2pdd2));
660  de_dT = -_Rco2 * (delta * tau * d2pddt * (dpdd - tau * d2pddt) / (2.0 * dpdd + delta * d2pdd2) +
661  tau * tau * d2phiSW_dt2(delta, tau));
662 }
663 
664 void
666  Real temperature,
667  Real & rho,
668  Real & drho_dp,
669  Real & drho_dT,
670  Real & e,
671  Real & de_dp,
672  Real & de_dT) const
673 {
674  rho_dpT(pressure, temperature, rho, drho_dp, drho_dT);
675  e_dpT(pressure, temperature, e, de_dp, de_dT);
676 }
677 
678 Real
680 {
681  // Require density first
682  Real density = rho(pressure, temperature);
683  // Scale the input density and temperature
684  Real delta = density / _critical_density;
685  Real tau = _critical_temperature / temperature;
686 
687  Real speed2 =
688  2.0 * delta * dphiSW_dd(delta, tau) + delta * delta * d2phiSW_dd2(delta, tau) -
689  std::pow(delta * dphiSW_dd(delta, tau) - delta * tau * d2phiSW_ddt(delta, tau), 2.0) /
690  (tau * tau * d2phiSW_dt2(delta, tau));
691 
692  return std::sqrt(_Rco2 * temperature * speed2);
693 }
694 
695 Real
697 {
698  // Require density first
699  Real density = rho(pressure, temperature);
700  // Scale the input density and temperature
701  Real delta = density / _critical_density;
702  Real tau = _critical_temperature / temperature;
703 
704  Real heat_capacity =
705  -tau * tau * d2phiSW_dt2(delta, tau) +
706  std::pow(delta * dphiSW_dd(delta, tau) - delta * tau * d2phiSW_ddt(delta, tau), 2.0) /
707  (2.0 * delta * dphiSW_dd(delta, tau) + delta * delta * d2phiSW_dd2(delta, tau));
708 
709  return _Rco2 * heat_capacity;
710 }
711 
712 Real
714 {
715  // Require density first
716  Real density = rho(pressure, temperature);
717  // Scale the input density and temperature
718  Real delta = density / _critical_density;
719  Real tau = _critical_temperature / temperature;
720 
721  return -_Rco2 * tau * tau * d2phiSW_dt2(delta, tau);
722 }
723 
724 Real
726 {
727  Real rho = this->rho(pressure, temperature);
728  return this->k_from_rho_T(rho, temperature);
729 }
730 
731 void
733  Real /*pressure*/, Real /*temperature*/, Real & /*k*/, Real & /*dk_dp*/, Real & /*dk_dT*/) const
734 {
735  mooseError(name(), ": k_dpT() is not implemented");
736 }
737 
738 Real
740 {
741  // Check the temperature is in the range of validity (216.592 K <= T <= 1000 K)
742  if (temperature <= _triple_point_temperature || temperature >= 1000.0)
743  mooseError("Temperature ", temperature, "K out of range (200K, 1000K) in ", name(), ": k()");
744 
745  const std::vector<Real> g1{0.0, 0.0, 1.5};
746  const std::vector<Real> g2{0.0, 1.0, 1.5, 1.5, 1.5, 3.5, 5.5};
747  const std::vector<Real> h1{1.0, 5.0, 1.0};
748  const std::vector<Real> h2{1.0, 2.0, 0.0, 5.0, 9.0, 0.0, 0.0};
749  const std::vector<Real> n1{7.69857587, 0.159885811, 1.56918621};
750  const std::vector<Real> n2{
751  -6.73400790, 16.3890156, 3.69415242, 22.3205514, 66.1420950, -0.171779133, 0.00433043347};
752  const std::vector<Real> a{
753  3.0, 6.70697, 0.94604, 0.3, 0.3, 0.39751, 0.33791, 0.77963, 0.79857, 0.9, 0.02, 0.2};
754 
755  // Scaled variables
756  Real Tr = temperature / _critical_temperature;
757  Real rhor = density / _critical_density;
758 
759  Real sum1 = 0.0;
760  for (std::size_t i = 0; i < n1.size(); ++i)
761  sum1 += n1[i] * std::pow(Tr, g1[i]) * std::pow(rhor, h1[i]);
762 
763  Real sum2 = 0.0;
764  for (std::size_t i = 0; i < n2.size(); ++i)
765  sum2 += n2[i] * std::pow(Tr, g2[i]) * std::pow(rhor, h2[i]);
766 
767  // Near-critical enhancement
768  Real alpha = 1.0 - a[9] * std::acosh(1.0 + a[10] * std::pow(std::pow(1.0 - Tr, 2.0), a[11]));
769  Real lambdac =
770  rhor * std::exp(-std::pow(rhor, a[0]) / a[0] - std::pow(a[1] * (Tr - 1.0), 2.0) -
771  std::pow(a[2] * (rhor - 1.0), 2.0)) /
772  std::pow(std::pow(std::pow(1.0 - 1.0 / Tr +
773  a[3] * std::pow(std::pow(rhor - 1.0, 2.0), 1.0 / (2.0 * a[4])),
774  2.0),
775  a[5]) +
776  std::pow(std::pow(a[6] * (rhor - alpha), 2.0), a[7]),
777  a[8]);
778 
779  return 4.81384 * (sum1 + std::exp(-5.0 * rhor * rhor) * sum2 + 0.775547504 * lambdac) / 1000.0;
780 }
781 
782 Real
784 {
785  // Require density first
786  Real density = rho(pressure, temperature);
787  // Scale the input density and temperature
788  Real delta = density / _critical_density;
789  Real tau = _critical_temperature / temperature;
790 
791  return _Rco2 * (tau * dphiSW_dt(delta, tau) - phiSW(delta, tau));
792 }
793 
794 Real
796 {
797  // Require density first
798  Real density = rho(pressure, temperature);
799  // Scale the input density and temperature
800  Real delta = density / _critical_density;
801  Real tau = _critical_temperature / temperature;
802 
803  return _Rco2 * temperature * (tau * dphiSW_dt(delta, tau) + delta * dphiSW_dd(delta, tau));
804 }
805 
806 void
808  Real pressure, Real temperature, Real & h, Real & dh_dp, Real & dh_dT) const
809 {
810  // Require density first
811  Real density = rho(pressure, temperature);
812  // Scale the input density and temperature
813  Real delta = density / _critical_density;
814  Real tau = _critical_temperature / temperature;
815  Real dpdd = dphiSW_dd(delta, tau);
816  Real d2pdd2 = d2phiSW_dd2(delta, tau);
817  Real d2pddt = d2phiSW_ddt(delta, tau);
818 
819  h = _Rco2 * temperature * (tau * dphiSW_dt(delta, tau) + delta * dpdd);
820  dh_dp = (dpdd + delta * d2pdd2 + tau * d2pddt) / (density * (2.0 * dpdd + delta * d2pdd2));
821  dh_dT = _Rco2 * delta * dpdd * (1.0 - tau * d2pddt / dpdd) * (1.0 - tau * d2pddt / dpdd) /
822  (2.0 + delta * d2pdd2 / dpdd) -
823  _Rco2 * tau * tau * d2phiSW_dt2(delta, tau);
824 }
825 
826 Real CO2FluidProperties::beta(Real /*pressure*/, Real /*temperature*/) const
827 {
828  mooseError(name(), ": beta() not implemented yet");
829 }
std::vector< Real > _n4
std::vector< Real > _t1
virtual Real beta(Real pressure, Real temperature) const override
Thermal expansion coefficient.
std::vector< Real > _d3
std::vector< Real > _n2
virtual Real k_from_rho_T(Real density, Real temperature) const override
Thermal conductivity as a function of density and temperature.
Real criticalDensity() const
CO2 critical density.
virtual Real cp(Real pressure, Real temperature) const override
Isobaric specific heat capacity.
std::vector< Real > _a4
Real criticalTemperature() const
CO2 critical temperature.
std::vector< Real > _c2
std::vector< Real > _beta4
std::vector< Real > _a0
Coefficients for the ideal gas component of the Helmholtz free energy.
std::vector< Real > _D4
virtual void henryConstant_dT(Real temperature, Real &Kh, Real &dKh_dT) const override
Henry&#39;s law constant for dissolution in water and derivative wrt temperature.
const Real _critical_density
Critical density (kg/m^3)
Real molarMass() const override
Molar mass.
virtual Real cv(Real pressure, Real temperature) const override
Isochoric specific heat.
const Real _critical_pressure
Critical pressure (Pa)
const std::string density
Definition: NS.h:15
virtual Real mu(Real pressure, Real temperature) const override
virtual void h_dpT(Real pressure, Real temperature, Real &h, Real &dh_dp, Real &dh_dT) const override
Enthalpy and its derivatives wrt pressure and temperature.
virtual Real c(Real pressure, Real temperature) const override
Speed of sound.
std::vector< Real > _gamma3
const std::string temperature
Definition: NS.h:25
std::vector< Real > _t3
std::vector< Real > _n3
virtual Real rho(Real pressure, Real temperature) const override
Density.
Real d2phiSW_dt2(Real delta, Real tau) const
Second derivative of Helmholtz free energy wrt tau.
Real d2phiSW_dd2(Real delta, Real tau) const
Second derivative of Helmholtz free energy wrt delta.
std::vector< Real > _b4
Common class for single phase fluid properties using a pressure and temperature formulation.
virtual void mu_dpT(Real pressure, Real temperature, Real &mu, Real &dmu_dp, Real &dmu_dT) const override
InputParameters validParams< CO2FluidProperties >()
virtual Real mu_from_rho_T(Real density, Real temperature) const override
CO2FluidProperties(const InputParameters &parameters)
virtual std::string fluidName() const override
Fluid name.
InputParameters validParams< SinglePhaseFluidPropertiesPT >()
std::vector< Real > _theta0
Real phiSW(Real delta, Real tau) const
Helmholtz free energy for CO2 From Span and Wagner (reference above)
std::vector< Real > _A4
virtual Real h(Real p, Real T) const override
Specific enthalpy.
Real criticalPressure() const
CO2 critical pressure.
virtual Real k(Real pressure, Real temperature) const override
Thermal conductivity.
std::vector< Real > _alpha3
Real triplePointTemperature() const
CO2 triple point temperature.
const Real _T_c2k
Conversion of temperature from Celcius to Kelvin.
std::vector< Real > _n1
Coefficients for the residual component of the Helmholtz free energy.
Real sublimationPressure(Real temperature) const
Sublimation pressure.
const Real _Rco2
Specific gas constant (J/mol/K)
Real root(std::function< Real(Real)> const &f, Real x1, Real x2, Real tol=1.0e-12)
Finds the root of a function using Brent&#39;s method.
Definition: BrentsMethod.C:58
std::vector< Real > _C4
Real dphiSW_dt(Real delta, Real tau) const
Derivative of Helmholtz free energy wrt tau.
virtual void mu_drhoT_from_rho_T(Real density, Real temperature, Real ddensity_dT, Real &mu, Real &dmu_drho, Real &dmu_dT) const override
Dynamic viscosity and its derivatives wrt density and temperature.
Real dphiSW_dd(Real delta, Real tau) const
Derivative of Helmholtz free energy wrt delta.
std::vector< Real > _d1
Real d2phiSW_ddt(Real delta, Real tau) const
Second derivative of Helmholtz free energy wrt delta and tau.
std::vector< Real > _t2
virtual void rho_dpT(Real pressure, Real temperature, Real &rho, Real &drho_dp, Real &drho_dT) const override
Density and its derivatives wrt pressure and temperature.
const Real _triple_point_temperature
Triple point temperature (K)
Real saturatedLiquidDensity(Real temperature) const
Saturated liquid density of CO2 Valid for temperatures between the triple point temperature and criti...
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
Real pressure(Real density, Real temperature) const
Pressure as a function of density and temperature From Span and Wagner (reference above) ...
virtual Real e(Real pressure, Real temperature) const override
Internal energy.
virtual void henryConstantIAPWS_dT(Real temperature, Real &Kh, Real &dKh_dT, Real A, Real B, Real C) const
IAPWS formulation of Henry&#39;s law constant for dissolution in water and derivative wrt temperature...
const Real _critical_temperature
Critical temperature (K)
Real pressureSW(Real density, Real temperature) const
Internal function to calculate pressure as a function of density and temperature using the Span and W...
virtual void k_dpT(Real pressure, Real temperature, Real &k, Real &dk_dp, Real &dk_dT) const override
Thermal conductivity and its derivatives wrt pressure and temperature.
Real vaporPressure(Real temperature) const
Vapor pressure.
std::vector< Real > _d2
std::vector< Real > _beta3
virtual Real henryConstantIAPWS(Real temperature, Real A, Real B, Real C) const
IAPWS formulation of Henry&#39;s law constant for dissolution in water From Guidelines on the Henry&#39;s con...
Real triplePointPressure() const
CO2 triple point pressure.
virtual Real henryConstant(Real temperature) const override
Henry&#39;s law constant for dissolution in water.
const Real _triple_point_pressure
Triple point pressure (Pa)
virtual void e_dpT(Real pressure, Real temperature, Real &e, Real &de_dp, Real &de_dT) const override
Internal energy and its derivatives wrt pressure and temperature.
Real saturatedVaporDensity(Real temperature) const
Saturated vapor density of CO2 Valid for temperatures between the triple point temperature and critic...
std::vector< Real > _B4
Real meltingPressure(Real temperature) const
Melting pressure.
virtual void rho_e_dpT(Real pressure, Real temperature, Real &rho, Real &drho_dp, Real &drho_dT, Real &e, Real &de_dp, Real &de_dT) const override
Density and internal energy and their derivatives wrt pressure and temperature.
std::vector< Real > _eps3
virtual Real s(Real pressure, Real temperature) const override
Specific entropy.
const Real _Mco2
Molar mass of CO2 (kg/mol)
Real partialDensity(Real temperature) const
Partial density of dissolved CO2 From Garcia, Density of aqueous solutions of CO2, LBNL-49023 (2001)
void bracket(std::function< Real(Real)> const &f, Real &x1, Real &x2)
Function to bracket a root of a given function.
Definition: BrentsMethod.C:14