Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "CO2FluidProperties.h"
11 : #include "BrentsMethod.h"
12 : #include "Conversion.h"
13 : #include "MathUtils.h"
14 : #include "libmesh/utility.h"
15 :
16 : registerMooseObject("FluidPropertiesApp", CO2FluidProperties);
17 :
18 : InputParameters
19 201 : CO2FluidProperties::validParams()
20 : {
21 201 : InputParameters params = HelmholtzFluidProperties::validParams();
22 201 : params.addClassDescription(
23 : "Fluid properties for carbon dioxide (CO2) using the Span & Wagner EOS");
24 201 : return params;
25 0 : }
26 :
27 107 : CO2FluidProperties::CO2FluidProperties(const InputParameters & parameters)
28 107 : : HelmholtzFluidProperties(parameters)
29 : {
30 107 : }
31 :
32 89 : CO2FluidProperties::~CO2FluidProperties() {}
33 :
34 : std::string
35 8 : CO2FluidProperties::fluidName() const
36 : {
37 8 : return "co2";
38 : }
39 :
40 : Real
41 57449729 : CO2FluidProperties::molarMass() const
42 : {
43 57449729 : return _Mco2;
44 : }
45 :
46 : Real
47 3 : CO2FluidProperties::criticalPressure() const
48 : {
49 3 : return _critical_pressure;
50 : }
51 :
52 : Real
53 57449730 : CO2FluidProperties::criticalTemperature() const
54 : {
55 57449730 : return _critical_temperature;
56 : }
57 :
58 : Real
59 57449730 : CO2FluidProperties::criticalDensity() const
60 : {
61 57449730 : return _critical_density;
62 : }
63 :
64 : Real
65 3 : CO2FluidProperties::triplePointPressure() const
66 : {
67 3 : return _triple_point_pressure;
68 : }
69 :
70 : Real
71 3 : CO2FluidProperties::triplePointTemperature() const
72 : {
73 3 : return _triple_point_temperature;
74 : }
75 :
76 : Real
77 4908974 : CO2FluidProperties::meltingPressure(Real temperature) const
78 : {
79 4908974 : if (temperature < _triple_point_temperature)
80 0 : throw MooseException("Temperature is below the triple point temperature in " + name() +
81 0 : ": meltingPressure()");
82 :
83 4908974 : Real Tstar = temperature / _triple_point_temperature;
84 :
85 4908974 : return _triple_point_pressure *
86 4908974 : (1.0 + 1955.539 * (Tstar - 1.0) + 2055.4593 * Utility::pow<2>(Tstar - 1.0));
87 : }
88 :
89 : Real
90 1 : CO2FluidProperties::sublimationPressure(Real temperature) const
91 : {
92 1 : if (temperature > _triple_point_temperature)
93 0 : throw MooseException("Temperature is above the triple point temperature in " + name() +
94 0 : ": sublimationPressure()");
95 :
96 1 : Real Tstar = temperature / _triple_point_temperature;
97 :
98 1 : Real pressure = _triple_point_pressure *
99 1 : std::exp((-14.740846 * (1.0 - Tstar) + 2.4327015 * std::pow(1.0 - Tstar, 1.9) -
100 1 : 5.3061778 * std::pow(1.0 - Tstar, 2.9)) /
101 1 : Tstar);
102 :
103 1 : return pressure;
104 : }
105 :
106 : Real
107 6993 : CO2FluidProperties::vaporPressure(Real temperature) const
108 : {
109 6993 : if (temperature < _triple_point_temperature || temperature > _critical_temperature)
110 0 : throw MooseException("Temperature is out of range in " + name() + ": vaporPressure()");
111 :
112 6993 : Real Tstar = temperature / _critical_temperature;
113 :
114 : Real logpressure =
115 6993 : (-7.0602087 * (1.0 - Tstar) + 1.9391218 * std::pow(1.0 - Tstar, 1.5) -
116 6993 : 1.6463597 * Utility::pow<2>(1.0 - Tstar) - 3.2995634 * Utility::pow<4>(1.0 - Tstar)) /
117 6993 : Tstar;
118 :
119 6993 : return _critical_pressure * std::exp(logpressure);
120 : }
121 :
122 : void
123 0 : CO2FluidProperties::vaporPressure(Real, Real &, Real &) const
124 : {
125 0 : mooseError("vaporPressure() is not implemented");
126 : }
127 :
128 : Real
129 4510570 : CO2FluidProperties::saturatedLiquidDensity(Real temperature) const
130 : {
131 4510570 : if (temperature < _triple_point_temperature || temperature > _critical_temperature)
132 0 : throw MooseException("Temperature is out of range in " + name() + ": saturatedLiquiDensity()");
133 :
134 4510570 : Real Tstar = temperature / _critical_temperature;
135 :
136 4510570 : Real logdensity = 1.9245108 * std::pow(1.0 - Tstar, 0.34) -
137 4510570 : 0.62385555 * std::pow(1.0 - Tstar, 0.5) -
138 4510570 : 0.32731127 * std::pow(1.0 - Tstar, 10.0 / 6.0) +
139 4510570 : 0.39245142 * std::pow(1.0 - Tstar, 11.0 / 6.0);
140 :
141 4510570 : return _critical_density * std::exp(logdensity);
142 : }
143 :
144 : Real
145 4510570 : CO2FluidProperties::saturatedVaporDensity(Real temperature) const
146 : {
147 4510570 : if (temperature < _triple_point_temperature || temperature > _critical_temperature)
148 0 : throw MooseException("Temperature is out of range in " + name() + ": saturatedVaporDensity()");
149 :
150 4510570 : Real Tstar = temperature / _critical_temperature;
151 :
152 : Real logdensity =
153 4510570 : (-1.7074879 * std::pow(1.0 - Tstar, 0.34) - 0.82274670 * std::pow(1.0 - Tstar, 0.5) -
154 4510570 : 4.6008549 * (1.0 - Tstar) - 10.111178 * std::pow(1.0 - Tstar, 7.0 / 3.0) -
155 4510570 : 29.742252 * std::pow(1.0 - Tstar, 14.0 / 3.0));
156 :
157 4510570 : return _critical_density * std::exp(logdensity);
158 : }
159 :
160 : Real
161 540996 : CO2FluidProperties::alpha(Real delta, Real tau) const
162 : {
163 : // Ideal gas component of the Helmholtz free energy
164 : Real sum0 = 0.0;
165 3245976 : for (std::size_t i = 0; i < _a0.size(); ++i)
166 2704980 : sum0 += _a0[i] * std::log(1.0 - std::exp(-_theta0[i] * tau));
167 :
168 540996 : Real phi0 = std::log(delta) + 8.37304456 - 3.70454304 * tau + 2.5 * std::log(tau) + sum0;
169 :
170 : // Residual component of the Helmholtz free energy
171 : Real theta, Delta, Psi;
172 : Real phir = 0.0;
173 4327968 : for (std::size_t i = 0; i < _n1.size(); ++i)
174 3786972 : phir += _n1[i] * MathUtils::pow(delta, _d1[i]) * std::pow(tau, _t1[i]);
175 :
176 15147888 : for (std::size_t i = 0; i < _n2.size(); ++i)
177 14606892 : phir += _n2[i] * MathUtils::pow(delta, _d2[i]) * std::pow(tau, _t2[i]) *
178 14606892 : std::exp(-MathUtils::pow(delta, _c2[i]));
179 :
180 3245976 : for (std::size_t i = 0; i < _n3.size(); ++i)
181 2704980 : phir += _n3[i] * MathUtils::pow(delta, _d3[i]) * MathUtils::pow(tau, _t3[i]) *
182 2704980 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _eps3[i]) -
183 2704980 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i]));
184 :
185 2163984 : for (std::size_t i = 0; i < _n4.size(); ++i)
186 : {
187 1622988 : theta = 1.0 - tau + _A4[i] * std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]));
188 1622988 : Delta = Utility::pow<2>(theta) + _B4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i]);
189 1622988 : Psi = std::exp(-_C4[i] * Utility::pow<2>(delta - 1.0) - _D4[i] * Utility::pow<2>(tau - 1.0));
190 1622988 : phir += _n4[i] * std::pow(Delta, _b4[i]) * delta * Psi;
191 : }
192 :
193 : // The Helmholtz free energy is the sum of these components
194 540996 : return phi0 + phir;
195 : }
196 :
197 : Real
198 55816739 : CO2FluidProperties::dalpha_ddelta(Real delta, Real tau) const
199 : {
200 : // Derivative of the ideal gas component wrt gamma
201 55816739 : Real dphi0dd = 1.0 / delta;
202 :
203 : // Derivative of the residual component wrt gamma
204 : Real theta, Delta, Psi, dDelta_dd, dPsi_dd;
205 : Real dphirdd = 0.0;
206 :
207 446533912 : for (std::size_t i = 0; i < _n1.size(); ++i)
208 390717173 : dphirdd += _n1[i] * _d1[i] * MathUtils::pow(delta, _d1[i] - 1.0) * std::pow(tau, _t1[i]);
209 :
210 1562868692 : for (std::size_t i = 0; i < _n2.size(); ++i)
211 1507051953 : dphirdd += _n2[i] * std::exp(-MathUtils::pow(delta, _c2[i])) *
212 1507051953 : (MathUtils::pow(delta, _d2[i] - 1.0) * std::pow(tau, _t2[i]) *
213 1507051953 : (_d2[i] - _c2[i] * MathUtils::pow(delta, _c2[i])));
214 :
215 334900434 : for (std::size_t i = 0; i < _n3.size(); ++i)
216 279083695 : dphirdd += _n3[i] * MathUtils::pow(delta, _d3[i]) * MathUtils::pow(tau, _t3[i]) *
217 279083695 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _eps3[i]) -
218 279083695 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
219 279083695 : (_d3[i] / delta - 2.0 * _alpha3[i] * (delta - _eps3[i]));
220 :
221 223266956 : for (std::size_t i = 0; i < _n4.size(); ++i)
222 : {
223 167450217 : theta = 1.0 - tau + _A4[i] * std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]));
224 167450217 : Delta = Utility::pow<2>(theta) + _B4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i]);
225 167450217 : Psi = std::exp(-_C4[i] * Utility::pow<2>(delta - 1.0) - _D4[i] * Utility::pow<2>(tau - 1.0));
226 167450217 : dPsi_dd = -2.0 * _C4[i] * (delta - 1.0) * Psi;
227 167450217 : dDelta_dd = (delta - 1.0) *
228 167450217 : (_A4[i] * theta * 2.0 / _beta4[i] *
229 167450217 : std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]) - 1.0) +
230 167450217 : 2.0 * _B4[i] * _a4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i] - 1.0));
231 :
232 167450217 : dphirdd += _n4[i] * (std::pow(Delta, _b4[i]) * (Psi + delta * dPsi_dd) +
233 167450217 : _b4[i] * std::pow(Delta, _b4[i] - 1.0) * dDelta_dd * delta * Psi);
234 : }
235 :
236 : // The derivative of the free energy wrt delta is the sum of these components
237 55816739 : return dphi0dd + dphirdd;
238 : }
239 :
240 : Real
241 1643001 : CO2FluidProperties::dalpha_dtau(Real delta, Real tau) const
242 : {
243 : // Derivative of the ideal gas component wrt tau
244 : Real sum0 = 0.0;
245 9858006 : for (std::size_t i = 0; i < _a0.size(); ++i)
246 8215005 : sum0 += _a0[i] * _theta0[i] * (1.0 / (1.0 - std::exp(-_theta0[i] * tau)) - 1.0);
247 :
248 1643001 : Real dphi0dt = -3.70454304 + 2.5 / tau + sum0;
249 :
250 : // Derivative of the residual component wrt tau
251 : Real theta, Delta, Psi, dDelta_dt, dPsi_dt;
252 : Real dphirdt = 0.0;
253 13144008 : for (std::size_t i = 0; i < _n1.size(); ++i)
254 11501007 : dphirdt += _n1[i] * _t1[i] * MathUtils::pow(delta, _d1[i]) * std::pow(tau, _t1[i] - 1.0);
255 :
256 46004028 : for (std::size_t i = 0; i < _n2.size(); ++i)
257 44361027 : dphirdt += _n2[i] * _t2[i] * MathUtils::pow(delta, _d2[i]) * std::pow(tau, _t2[i] - 1.0) *
258 44361027 : std::exp(-MathUtils::pow(delta, _c2[i]));
259 :
260 9858006 : for (std::size_t i = 0; i < _n3.size(); ++i)
261 8215005 : dphirdt += _n3[i] * MathUtils::pow(delta, _d3[i]) * MathUtils::pow(tau, _t3[i]) *
262 8215005 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _eps3[i]) -
263 8215005 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
264 8215005 : (_t3[i] / tau - 2.0 * _beta3[i] * (tau - _gamma3[i]));
265 :
266 6572004 : for (std::size_t i = 0; i < _n4.size(); ++i)
267 : {
268 4929003 : theta = 1.0 - tau + _A4[i] * std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]));
269 4929003 : Delta = Utility::pow<2>(theta) + _B4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i]);
270 4929003 : Psi = std::exp(-_C4[i] * Utility::pow<2>(delta - 1.0) - _D4[i] * Utility::pow<2>(tau - 1.0));
271 4929003 : dDelta_dt = -2.0 * theta * _b4[i] * std::pow(Delta, _b4[i] - 1.0);
272 4929003 : dPsi_dt = -2.0 * _D4[i] * (tau - 1.0) * Psi;
273 :
274 4929003 : dphirdt += _n4[i] * delta * (Psi * dDelta_dt + std::pow(Delta, _b4[i]) * dPsi_dt);
275 : }
276 :
277 : // The derivative of the free energy wrt tau is the sum of these components
278 1643001 : return dphi0dt + dphirdt;
279 : }
280 :
281 : Real
282 1081970 : CO2FluidProperties::d2alpha_ddelta2(Real delta, Real tau) const
283 : {
284 : // Second derivative of the ideal gas component wrt gamma
285 1081970 : Real d2phi0dd2 = -1.0 / delta / delta;
286 :
287 : // Second derivative of the residual component wrt gamma
288 : Real d2phirdd2 = 0.0;
289 : Real theta, Delta, Psi, dDelta_dd, dPsi_dd, d2Delta_dd2, d2Psi_dd2;
290 :
291 8655760 : for (std::size_t i = 0; i < _n1.size(); ++i)
292 7573790 : d2phirdd2 += _n1[i] * _d1[i] * (_d1[i] - 1.0) * MathUtils::pow(delta, _d1[i] - 2) *
293 7573790 : std::pow(tau, _t1[i]);
294 :
295 30295160 : for (std::size_t i = 0; i < _n2.size(); ++i)
296 29213190 : d2phirdd2 += _n2[i] * std::exp(-MathUtils::pow(delta, _c2[i])) *
297 29213190 : MathUtils::pow(delta, _d2[i] - 2) * std::pow(tau, _t2[i]) *
298 29213190 : ((_d2[i] - _c2[i] * MathUtils::pow(delta, _c2[i])) *
299 29213190 : (_d2[i] - 1.0 - _c2[i] * MathUtils::pow(delta, _c2[i])) -
300 29213190 : _c2[i] * _c2[i] * MathUtils::pow(delta, _c2[i]));
301 :
302 6491820 : for (std::size_t i = 0; i < _n3.size(); ++i)
303 5409850 : d2phirdd2 +=
304 5409850 : _n3[i] * MathUtils::pow(tau, _t3[i]) *
305 5409850 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _eps3[i]) -
306 5409850 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
307 5409850 : (-2.0 * _alpha3[i] * MathUtils::pow(delta, _d3[i]) +
308 5409850 : 4.0 * _alpha3[i] * _alpha3[i] * MathUtils::pow(delta, _d3[i]) *
309 5409850 : Utility::pow<2>(delta - _eps3[i]) -
310 5409850 : 4.0 * _d3[i] * _alpha3[i] * MathUtils::pow(delta, _d3[i] - 1.0) * (delta - _eps3[i]) +
311 5409850 : _d3[i] * (_d3[i] - 1.0) * MathUtils::pow(delta, _d3[i] - 2.0));
312 :
313 4327880 : for (std::size_t i = 0; i < _n4.size(); ++i)
314 : {
315 3245910 : theta = 1.0 - tau + _A4[i] * std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]));
316 3245910 : Delta = Utility::pow<2>(theta) + _B4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i]);
317 3245910 : Psi = std::exp(-_C4[i] * Utility::pow<2>(delta - 1.0) - _D4[i] * Utility::pow<2>(tau - 1.0));
318 3245910 : dPsi_dd = -2.0 * _C4[i] * (delta - 1.0) * Psi;
319 3245910 : dDelta_dd = (delta - 1.0) *
320 3245910 : (_A4[i] * theta * 2.0 / _beta4[i] *
321 3245910 : std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]) - 1.0) +
322 3245910 : 2.0 * _B4[i] * _a4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i] - 1.0));
323 3245910 : d2Psi_dd2 = 3.0 * _D4[i] * Psi * (2.0 * _C4[i] * Utility::pow<2>(delta - 1.0) - 1.0);
324 3245910 : d2Delta_dd2 = 1.0 / (delta - 1.0) * dDelta_dd +
325 3245910 : (delta - 1.0) * (delta - 1.0) *
326 3245910 : (4.0 * _B4[i] * _a4[i] * (_a4[i] - 1.0) *
327 3245910 : std::pow(Utility::pow<2>(delta - 1.0), _a4[i] - 2.0) +
328 3245910 : 2.0 * _A4[i] * _A4[i] *
329 3245910 : Utility::pow<2>(std::pow(Utility::pow<2>(delta - 1.0),
330 3245910 : 1.0 / (2.0 * _beta4[i]) - 1.0)) /
331 3245910 : _beta4[i] / _beta4[i] +
332 3245910 : (4.0 / _beta4[i]) * _A4[i] * theta * (1.0 / (2.0 * _beta4[i]) - 1.0) *
333 3245910 : std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]) - 2.0));
334 3245910 : d2phirdd2 +=
335 3245910 : _n4[i] *
336 3245910 : (std::pow(Delta, _b4[i]) * (2.0 * dPsi_dd + delta * d2Psi_dd2) +
337 3245910 : 2.0 * _b4[i] * std::pow(Delta, _b4[i] - 1.0) * dDelta_dd * (Psi + delta * dPsi_dd) +
338 3245910 : _b4[i] *
339 3245910 : (std::pow(Delta, _b4[i] - 1.0) * d2Delta_dd2 +
340 3245910 : (_b4[i] - 1.0) * std::pow(Delta, _b4[i] - 2.0) * Utility::pow<2>(dDelta_dd)) *
341 3245910 : delta * Psi);
342 : }
343 : // The second derivative of the free energy wrt delta is the sum of these components
344 1081970 : return d2phi0dd2 + d2phirdd2;
345 : }
346 :
347 : Real
348 1622945 : CO2FluidProperties::d2alpha_dtau2(Real delta, Real tau) const
349 : {
350 : // Second derivative of the ideal gas component wrt tau
351 : Real sum0 = 0.0;
352 9737670 : for (std::size_t i = 0; i < _a0.size(); ++i)
353 8114725 : sum0 += _a0[i] * _theta0[i] * _theta0[i] * std::exp(-_theta0[i] * tau) /
354 8114725 : Utility::pow<2>(1.0 - std::exp(-_theta0[i] * tau));
355 :
356 1622945 : Real d2phi0dt2 = -2.5 / tau / tau - sum0;
357 :
358 : // Second derivative of the residual component wrt tau
359 : Real d2phirdt2 = 0.0;
360 : Real theta, Delta, Psi, dPsi_dt, dDelta_dt, d2Delta_dt2, d2Psi_dt2;
361 :
362 12983560 : for (std::size_t i = 0; i < _n1.size(); ++i)
363 11360615 : d2phirdt2 += _n1[i] * _t1[i] * (_t1[i] - 1.0) * MathUtils::pow(delta, _d1[i]) *
364 11360615 : std::pow(tau, _t1[i] - 2.0);
365 :
366 45442460 : for (std::size_t i = 0; i < _n2.size(); ++i)
367 43819515 : d2phirdt2 += _n2[i] * _t2[i] * (_t2[i] - 1.0) * MathUtils::pow(delta, _d2[i]) *
368 43819515 : std::exp(-MathUtils::pow(delta, _c2[i])) * std::pow(tau, _t2[i] - 2.0);
369 :
370 9737670 : for (std::size_t i = 0; i < _n3.size(); ++i)
371 8114725 : d2phirdt2 += _n3[i] * MathUtils::pow(delta, _d3[i]) * MathUtils::pow(tau, _t3[i]) *
372 8114725 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _eps3[i]) -
373 8114725 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
374 8114725 : (Utility::pow<2>(_t3[i] / tau - 2.0 * _beta3[i] * (tau - _gamma3[i])) -
375 8114725 : _t3[i] / tau / tau - 2.0 * _beta3[i]);
376 :
377 6491780 : for (std::size_t i = 0; i < _n4.size(); ++i)
378 : {
379 4868835 : theta = 1.0 - tau + _A4[i] * std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]));
380 4868835 : Delta = Utility::pow<2>(theta) + _B4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i]);
381 4868835 : Psi = std::exp(-_C4[i] * Utility::pow<2>(delta - 1.0) - _D4[i] * Utility::pow<2>(tau - 1.0));
382 4868835 : dDelta_dt = -2.0 * theta * _b4[i] * std::pow(Delta, _b4[i] - 1.0);
383 4868835 : d2Delta_dt2 = 2.0 * _b4[i] * std::pow(Delta, _b4[i] - 1.0) +
384 4868835 : 4.0 * theta * theta * _b4[i] * (_b4[i] - 1.0) * std::pow(Delta, _b4[i] - 2.0);
385 4868835 : dPsi_dt = -2.0 * _D4[i] * (tau - 1.0) * Psi;
386 4868835 : d2Psi_dt2 = 2.0 * _D4[i] * (2.0 * _D4[i] * (tau - 1.0) * (tau - 1.0) - 1.0) * Psi;
387 4868835 : d2phirdt2 +=
388 4868835 : _n4[i] * delta *
389 4868835 : (Psi * d2Delta_dt2 + 2.0 * dDelta_dt * dPsi_dt + std::pow(Delta, _b4[i]) * d2Psi_dt2);
390 : }
391 :
392 : // The second derivative of the free energy wrt tau is the sum of these components
393 1622945 : return d2phi0dt2 + d2phirdt2;
394 : }
395 :
396 : Real
397 1081970 : CO2FluidProperties::d2alpha_ddeltatau(Real delta, Real tau) const
398 : {
399 : // Note: second derivative of the ideal gas component wrt delta and tau is 0
400 : // Derivative of the residual component wrt gamma
401 : Real theta, Delta, Psi, dDelta_dd, dPsi_dd, dDelta_dt, dPsi_dt, d2Delta_ddt, d2Psi_ddt;
402 : Real d2phirddt = 0.0;
403 8655760 : for (std::size_t i = 0; i < _n1.size(); ++i)
404 7573790 : d2phirddt += _n1[i] * _d1[i] * _t1[i] * MathUtils::pow(delta, _d1[i] - 1.0) *
405 7573790 : std::pow(tau, _t1[i] - 1.0);
406 :
407 30295160 : for (std::size_t i = 0; i < _n2.size(); ++i)
408 29213190 : d2phirddt += _n2[i] * std::exp(-MathUtils::pow(delta, _c2[i])) *
409 29213190 : (MathUtils::pow(delta, _d2[i] - 1.0) * _t2[i] * std::pow(tau, _t2[i] - 1.0) *
410 29213190 : (_d2[i] - _c2[i] * MathUtils::pow(delta, _c2[i])));
411 :
412 6491820 : for (std::size_t i = 0; i < _n3.size(); ++i)
413 5409850 : d2phirddt += _n3[i] * MathUtils::pow(delta, _d3[i]) * MathUtils::pow(tau, _t3[i]) *
414 5409850 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _eps3[i]) -
415 5409850 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
416 5409850 : (_d3[i] / delta - 2.0 * _alpha3[i] * (delta - _eps3[i])) *
417 5409850 : (_t3[i] / tau - 2.0 * _beta3[i] * (tau - _gamma3[i]));
418 :
419 4327880 : for (std::size_t i = 0; i < _n4.size(); ++i)
420 : {
421 3245910 : theta = 1.0 - tau + _A4[i] * std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]));
422 3245910 : Delta = Utility::pow<2>(theta) + _B4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i]);
423 3245910 : Psi = std::exp(-_C4[i] * Utility::pow<2>(delta - 1.0) - _D4[i] * Utility::pow<2>(tau - 1.0));
424 3245910 : dPsi_dd = -2.0 * _C4[i] * (delta - 1.0) * Psi;
425 3245910 : dPsi_dt = -2.0 * _D4[i] * (tau - 1.0) * Psi;
426 3245910 : d2Psi_ddt = 4.0 * _C4[i] * _D4[i] * (delta - 1.0) * (tau - 1.0) * Psi;
427 3245910 : dDelta_dd = (delta - 1.0) *
428 3245910 : (_A4[i] * theta * 2.0 / _beta4[i] *
429 3245910 : std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]) - 1.0) +
430 3245910 : 2.0 * _B4[i] * _a4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i] - 1.0));
431 3245910 : dDelta_dt = -2.0 * theta * _b4[i] * std::pow(Delta, _b4[i] - 1.0);
432 3245910 : d2Delta_ddt = -2.0 * _A4[i] * _b4[i] / _beta4[i] * std::pow(Delta, _b4[i] - 1.0) *
433 3245910 : (delta - 1.0) *
434 3245910 : std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]) - 1.0) -
435 3245910 : 2.0 * theta * _b4[i] * (_b4[i] - 1.0) * std::pow(Delta, _b4[i] - 2.0) * dDelta_dd;
436 :
437 3245910 : d2phirddt += _n4[i] * (std::pow(Delta, _b4[i]) * (dPsi_dt + delta * d2Psi_ddt) +
438 3245910 : delta * _b4[i] * std::pow(Delta, _b4[i] - 1.0) * dDelta_dd * dPsi_dt +
439 3245910 : dDelta_dt * (Psi + delta * dPsi_dd) + d2Delta_ddt * delta * Psi);
440 : }
441 :
442 1081970 : return d2phirddt;
443 : }
444 :
445 : Real
446 54190754 : CO2FluidProperties::p_from_rho_T(Real density, Real temperature) const
447 : {
448 : // Check that the input parameters are within the region of validity
449 54190754 : if (temperature < 216.0 || temperature > 1100.0 || density <= 0.0)
450 0 : throw MooseException("Parameters out of range in " + name() + ": pressure()");
451 :
452 : Real pressure = 0.0;
453 :
454 54190754 : if (temperature > _triple_point_temperature && temperature < _critical_temperature)
455 : {
456 4510567 : Real gas_density = saturatedVaporDensity(temperature);
457 4510567 : Real liquid_density = saturatedLiquidDensity(temperature);
458 :
459 4510567 : if (density < gas_density || density > liquid_density)
460 4503579 : pressure = HelmholtzFluidProperties::p_from_rho_T(density, temperature);
461 : else
462 6988 : pressure = vaporPressure(temperature);
463 : }
464 : else
465 49680187 : pressure = HelmholtzFluidProperties::p_from_rho_T(density, temperature);
466 :
467 54190754 : return pressure;
468 : }
469 :
470 : Real
471 4908971 : CO2FluidProperties::rho_from_p_T(Real pressure, Real temperature) const
472 : {
473 : // Check that the input parameters are within the region of validity
474 4908971 : if (temperature < 216.0 || temperature > 1100.0 || pressure <= 0.0)
475 0 : throw MooseException("Parameters out of range in " + name() + ": rho_from_p_T()");
476 :
477 : // Also check that the pressure and temperature are not in the solid phase region
478 4908971 : if (((temperature > _triple_point_temperature) && (pressure > meltingPressure(temperature))) ||
479 4908971 : ((temperature < _triple_point_temperature) && (pressure > sublimationPressure(temperature))))
480 0 : throw MooseException("Input pressure and temperature in " + name() +
481 0 : ": rho_from_p_T() correspond to solid CO2 phase");
482 :
483 : Real density;
484 : // Initial estimate of a bracketing interval for the density
485 4908971 : Real lower_density = 100.0;
486 4908971 : Real upper_density = 1000.0;
487 :
488 : // The density is found by finding the zero of the pressure calculated using the
489 : // Span and Wagner EOS minus the input pressure
490 54190754 : auto pressure_diff = [&pressure, &temperature, this](Real x)
491 54190754 : { return p_from_rho_T(x, temperature) - pressure; };
492 :
493 4908971 : BrentsMethod::bracket(pressure_diff, lower_density, upper_density);
494 4908971 : density = BrentsMethod::root(pressure_diff, lower_density, upper_density);
495 :
496 4908971 : return density;
497 : }
498 :
499 : void
500 10 : CO2FluidProperties::rho_from_p_T(
501 : Real pressure, Real temperature, Real & rho, Real & drho_dp, Real & drho_dT) const
502 : {
503 10 : HelmholtzFluidProperties::rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT);
504 10 : }
505 :
506 : Real
507 551002 : CO2FluidProperties::mu_from_p_T(Real pressure, Real temperature) const
508 : {
509 551002 : Real rho = rho_from_p_T(pressure, temperature);
510 551002 : return mu_from_rho_T(rho, temperature);
511 : }
512 :
513 : void
514 7 : CO2FluidProperties::mu_from_p_T(
515 : Real pressure, Real temperature, Real & mu, Real & dmu_dp, Real & dmu_dT) const
516 : {
517 : Real rho, drho_dp, drho_dT;
518 7 : HelmholtzFluidProperties::rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT);
519 :
520 : Real dmu_drho;
521 7 : mu_from_rho_T(rho, temperature, drho_dT, mu, dmu_drho, dmu_dT);
522 7 : dmu_dp = dmu_drho * drho_dp;
523 7 : }
524 :
525 : Real
526 551011 : CO2FluidProperties::mu_from_rho_T(Real density, Real temperature) const
527 : {
528 : // Check that the input parameters are within the region of validity
529 551011 : if (temperature < 216.0 || temperature > 1000.0 || density > 1400.0)
530 0 : throw MooseException("Parameters out of range in " + name() + ": mu_from_rho_T()");
531 :
532 551011 : Real Tstar = temperature / 251.196;
533 :
534 : // Viscosity in the zero-density limit
535 : Real sum = 0.0;
536 :
537 3306066 : for (std::size_t i = 0; i < _mu_a.size(); ++i)
538 2755055 : sum += _mu_a[i] * MathUtils::pow(std::log(Tstar), i);
539 :
540 551011 : Real mu0 = 1.00697 * std::sqrt(temperature) / std::exp(sum);
541 :
542 : // Excess viscosity due to finite density
543 551011 : Real mue = _mu_d[0] * density + _mu_d[1] * Utility::pow<2>(density) +
544 551011 : _mu_d[2] * Utility::pow<6>(density) / Utility::pow<3>(Tstar) +
545 551011 : _mu_d[3] * Utility::pow<8>(density) + _mu_d[4] * Utility::pow<8>(density) / Tstar;
546 :
547 551011 : return (mu0 + mue) * 1.0e-6; // convert to Pa.s
548 : }
549 :
550 : void
551 10 : CO2FluidProperties::mu_from_rho_T(Real density,
552 : Real temperature,
553 : Real ddensity_dT,
554 : Real & mu,
555 : Real & dmu_drho,
556 : Real & dmu_dT) const
557 : {
558 : // Check that the input parameters are within the region of validity
559 10 : if (temperature < 216.0 || temperature > 1000.0 || density > 1400.0)
560 0 : throw MooseException("Parameters out of range in " + name() + ": mu_drhoT()");
561 :
562 10 : Real Tstar = temperature / 251.196;
563 : Real dTstar_dT = 1.0 / 251.196;
564 :
565 : // Viscosity in the zero-density limit. Note this is only a function of T.
566 : // Start the sum at i = 1 so the derivative is defined
567 10 : Real sum0 = _mu_a[0], dsum0_dTstar = 0.0;
568 :
569 50 : for (std::size_t i = 1; i < _mu_a.size(); ++i)
570 : {
571 40 : sum0 += _mu_a[i] * MathUtils::pow(std::log(Tstar), i);
572 40 : dsum0_dTstar += i * _mu_a[i] * MathUtils::pow(std::log(Tstar), i - 1) / Tstar;
573 : }
574 :
575 10 : Real mu0 = 1.00697 * std::sqrt(temperature) / std::exp(sum0);
576 10 : Real dmu0_dT = (0.5 * 1.00697 / std::sqrt(temperature) -
577 10 : 1.00697 * std::sqrt(temperature) * dsum0_dTstar * dTstar_dT) /
578 10 : std::exp(sum0);
579 :
580 : // Excess viscosity due to finite density
581 10 : Real mue = _mu_d[0] * density + _mu_d[1] * Utility::pow<2>(density) +
582 10 : _mu_d[2] * Utility::pow<6>(density) / Utility::pow<3>(Tstar) +
583 10 : _mu_d[3] * Utility::pow<8>(density) + _mu_d[4] * Utility::pow<8>(density) / Tstar;
584 :
585 10 : Real dmue_drho = _mu_d[0] + 2.0 * _mu_d[1] * density +
586 10 : 6.0 * _mu_d[2] * Utility::pow<5>(density) / Utility::pow<3>(Tstar) +
587 10 : 8.0 * _mu_d[3] * Utility::pow<7>(density) +
588 10 : 8.0 * _mu_d[4] * Utility::pow<7>(density) / Tstar;
589 :
590 10 : Real dmue_dT = (-3.0 * _mu_d[2] * Utility::pow<6>(density) / Utility::pow<4>(Tstar) -
591 10 : _mu_d[4] * Utility::pow<8>(density) / Tstar / Tstar) *
592 10 : dTstar_dT;
593 :
594 : // Viscosity in Pa.s is
595 10 : mu = (mu0 + mue) * 1.0e-6;
596 10 : dmu_drho = dmue_drho * 1.0e-6;
597 10 : dmu_dT = (dmu0_dT + dmue_dT) * 1.0e-6 + dmu_drho * ddensity_dT;
598 10 : }
599 :
600 : void
601 1 : CO2FluidProperties::rho_mu_from_p_T(Real pressure, Real temperature, Real & rho, Real & mu) const
602 : {
603 1 : rho = rho_from_p_T(pressure, temperature);
604 1 : mu = mu_from_rho_T(rho, temperature);
605 1 : }
606 :
607 : void
608 1 : CO2FluidProperties::rho_mu_from_p_T(Real pressure,
609 : Real temperature,
610 : Real & rho,
611 : Real & drho_dp,
612 : Real & drho_dT,
613 : Real & mu,
614 : Real & dmu_dp,
615 : Real & dmu_dT) const
616 : {
617 1 : rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT);
618 : Real dmu_drho;
619 1 : mu_from_rho_T(rho, temperature, drho_dT, mu, dmu_drho, dmu_dT);
620 1 : dmu_dp = dmu_drho * drho_dp;
621 1 : }
622 :
623 : Real
624 3 : CO2FluidProperties::partialDensity(Real temperature) const
625 : {
626 : // This correlation uses temperature in C
627 3 : Real Tc = temperature - _T_c2k;
628 : // The partial volume
629 3 : Real V = 37.51 - 9.585e-2 * Tc + 8.74e-4 * Tc * Tc - 5.044e-7 * Tc * Tc * Tc;
630 :
631 3 : return 1.0e6 * _Mco2 / V;
632 : }
633 :
634 : std::vector<Real>
635 7 : CO2FluidProperties::henryCoefficients() const
636 : {
637 7 : return {-8.55445, 4.01195, 9.52345};
638 : }
639 :
640 : Real
641 541006 : CO2FluidProperties::k_from_p_T(Real pressure, Real temperature) const
642 : {
643 : // Require density first
644 541006 : Real density = rho_from_p_T(pressure, temperature);
645 541006 : return k_from_rho_T(density, temperature);
646 : }
647 :
648 : void
649 3 : CO2FluidProperties::k_from_p_T(
650 : Real pressure, Real temperature, Real & k, Real & dk_dp, Real & dk_dT) const
651 : {
652 3 : k = this->k_from_p_T(pressure, temperature);
653 : // Calculate derivatives using finite differences. Note: this will be slow as
654 : // multiple calculations of density are required
655 : const Real eps = 1.0e-6;
656 3 : const Real peps = pressure * eps;
657 3 : const Real Teps = temperature * eps;
658 :
659 3 : dk_dp = (this->k_from_p_T(pressure + peps, temperature) - k) / peps;
660 3 : dk_dT = (this->k_from_p_T(pressure, temperature + Teps) - k) / Teps;
661 3 : }
662 :
663 : Real
664 541009 : CO2FluidProperties::k_from_rho_T(Real density, Real temperature) const
665 : {
666 : // Check the temperature is in the range of validity (216.592 K <= T <= 1000 K)
667 541009 : if (temperature <= _triple_point_temperature || temperature >= 1000.0)
668 0 : throw MooseException("Temperature " + Moose::stringify(temperature) +
669 0 : "K out of range (200K, 1000K) in " + name() + ": k()");
670 :
671 : // Scaled variables
672 541009 : Real Tr = temperature / _critical_temperature;
673 541009 : Real rhor = density / _critical_density;
674 :
675 : Real sum1 = 0.0;
676 2164036 : for (std::size_t i = 0; i < _k_n1.size(); ++i)
677 1623027 : sum1 += _k_n1[i] * std::pow(Tr, _k_g1[i]) * MathUtils::pow(rhor, _k_h1[i]);
678 :
679 : Real sum2 = 0.0;
680 4328072 : for (std::size_t i = 0; i < _k_n2.size(); ++i)
681 3787063 : sum2 += _k_n2[i] * std::pow(Tr, _k_g2[i]) * MathUtils::pow(rhor, _k_h2[i]);
682 :
683 : // Near-critical enhancement
684 : Real alpha =
685 541009 : 1.0 - _k_a[9] * std::acosh(1.0 + _k_a[10] * std::pow(Utility::pow<2>(1.0 - Tr), _k_a[11]));
686 : Real lambdac =
687 541009 : rhor *
688 541009 : std::exp(-std::pow(rhor, _k_a[0]) / _k_a[0] - Utility::pow<2>(_k_a[1] * (Tr - 1.0)) -
689 541009 : Utility::pow<2>(_k_a[2] * (rhor - 1.0))) /
690 541009 : std::pow(std::pow(Utility::pow<2>(
691 541009 : 1.0 - 1.0 / Tr +
692 541009 : _k_a[3] * std::pow(Utility::pow<2>(rhor - 1.0), 1.0 / (2.0 * _k_a[4]))),
693 : _k_a[5]) +
694 541009 : std::pow(Utility::pow<2>(_k_a[6] * (rhor - alpha)), _k_a[7]),
695 541009 : _k_a[8]);
696 :
697 541009 : return 4.81384 * (sum1 + std::exp(-5.0 * rhor * rhor) * sum2 + 0.775547504 * lambdac) / 1000.0;
698 : }
|
