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 268 : CO2FluidProperties::validParams()
20 : {
21 268 : InputParameters params = HelmholtzFluidProperties::validParams();
22 268 : params.addClassDescription(
23 : "Fluid properties for carbon dioxide (CO2) using the Span & Wagner EOS");
24 268 : return params;
25 0 : }
26 :
27 146 : CO2FluidProperties::CO2FluidProperties(const InputParameters & parameters)
28 146 : : HelmholtzFluidProperties(parameters)
29 : {
30 146 : }
31 :
32 256 : CO2FluidProperties::~CO2FluidProperties() {}
33 :
34 : std::string
35 10 : CO2FluidProperties::fluidName() const
36 : {
37 10 : return "co2";
38 : }
39 :
40 : Real
41 95275501 : CO2FluidProperties::molarMass() const
42 : {
43 95275501 : return _Mco2;
44 : }
45 :
46 : Real
47 3 : CO2FluidProperties::criticalPressure() const
48 : {
49 3 : return _critical_pressure;
50 : }
51 :
52 : Real
53 95275502 : CO2FluidProperties::criticalTemperature() const
54 : {
55 95275502 : return _critical_temperature;
56 : }
57 :
58 : Real
59 95275502 : CO2FluidProperties::criticalDensity() const
60 : {
61 95275502 : 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 8142410 : CO2FluidProperties::meltingPressure(Real temperature) const
78 : {
79 8142410 : if (temperature < _triple_point_temperature)
80 0 : throw MooseException("Temperature is below the triple point temperature in " + name() +
81 0 : ": meltingPressure()");
82 :
83 8142410 : Real Tstar = temperature / _triple_point_temperature;
84 :
85 8142410 : return _triple_point_pressure *
86 8142410 : (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 11061 : CO2FluidProperties::vaporPressure(Real temperature) const
108 : {
109 11061 : if (temperature < _triple_point_temperature || temperature > _critical_temperature)
110 0 : throw MooseException("Temperature is out of range in " + name() + ": vaporPressure()");
111 :
112 11061 : Real Tstar = temperature / _critical_temperature;
113 :
114 : Real logpressure =
115 11061 : (-7.0602087 * (1.0 - Tstar) + 1.9391218 * std::pow(1.0 - Tstar, 1.5) -
116 11061 : 1.6463597 * Utility::pow<2>(1.0 - Tstar) - 3.2995634 * Utility::pow<4>(1.0 - Tstar)) /
117 11061 : Tstar;
118 :
119 11061 : 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 7491478 : CO2FluidProperties::saturatedLiquidDensity(Real temperature) const
130 : {
131 7491478 : if (temperature < _triple_point_temperature || temperature > _critical_temperature)
132 0 : throw MooseException("Temperature is out of range in " + name() + ": saturatedLiquiDensity()");
133 :
134 7491478 : Real Tstar = temperature / _critical_temperature;
135 :
136 7491478 : Real logdensity = 1.9245108 * std::pow(1.0 - Tstar, 0.34) -
137 7491478 : 0.62385555 * std::pow(1.0 - Tstar, 0.5) -
138 7491478 : 0.32731127 * std::pow(1.0 - Tstar, 10.0 / 6.0) +
139 7491478 : 0.39245142 * std::pow(1.0 - Tstar, 11.0 / 6.0);
140 :
141 7491478 : return _critical_density * std::exp(logdensity);
142 : }
143 :
144 : Real
145 7491478 : CO2FluidProperties::saturatedVaporDensity(Real temperature) const
146 : {
147 7491478 : if (temperature < _triple_point_temperature || temperature > _critical_temperature)
148 0 : throw MooseException("Temperature is out of range in " + name() + ": saturatedVaporDensity()");
149 :
150 7491478 : Real Tstar = temperature / _critical_temperature;
151 :
152 : Real logdensity =
153 7491478 : (-1.7074879 * std::pow(1.0 - Tstar, 0.34) - 0.82274670 * std::pow(1.0 - Tstar, 0.5) -
154 7491478 : 4.6008549 * (1.0 - Tstar) - 10.111178 * std::pow(1.0 - Tstar, 7.0 / 3.0) -
155 7491478 : 29.742252 * std::pow(1.0 - Tstar, 14.0 / 3.0));
156 :
157 7491478 : return _critical_density * std::exp(logdensity);
158 : }
159 :
160 : Real
161 900264 : CO2FluidProperties::alpha(Real delta, Real tau) const
162 : {
163 : // Ideal gas component of the Helmholtz free energy
164 : Real sum0 = 0.0;
165 5401584 : for (std::size_t i = 0; i < _a0.size(); ++i)
166 4501320 : sum0 += _a0[i] * std::log(1.0 - std::exp(-_theta0[i] * tau));
167 :
168 900264 : 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 7202112 : for (std::size_t i = 0; i < _n1.size(); ++i)
174 6301848 : phir += _n1[i] * MathUtils::pow(delta, _d1[i]) * std::pow(tau, _t1[i]);
175 :
176 25207392 : for (std::size_t i = 0; i < _n2.size(); ++i)
177 24307128 : phir += _n2[i] * MathUtils::pow(delta, _d2[i]) * std::pow(tau, _t2[i]) *
178 24307128 : std::exp(-MathUtils::pow(delta, _c2[i]));
179 :
180 5401584 : for (std::size_t i = 0; i < _n3.size(); ++i)
181 4501320 : phir += _n3[i] * MathUtils::pow(delta, _d3[i]) * MathUtils::pow(tau, _t3[i]) *
182 4501320 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _eps3[i]) -
183 4501320 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i]));
184 :
185 3601056 : for (std::size_t i = 0; i < _n4.size(); ++i)
186 : {
187 2700792 : theta = 1.0 - tau + _A4[i] * std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]));
188 2700792 : Delta = Utility::pow<2>(theta) + _B4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i]);
189 2700792 : Psi = std::exp(-_C4[i] * Utility::pow<2>(delta - 1.0) - _D4[i] * Utility::pow<2>(tau - 1.0));
190 2700792 : phir += _n4[i] * std::pow(Delta, _b4[i]) * delta * Psi;
191 : }
192 :
193 : // The Helmholtz free energy is the sum of these components
194 900264 : return phi0 + phir;
195 : }
196 :
197 : Real
198 92564697 : CO2FluidProperties::dalpha_ddelta(Real delta, Real tau) const
199 : {
200 : // Derivative of the ideal gas component wrt gamma
201 92564697 : 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 740517576 : for (std::size_t i = 0; i < _n1.size(); ++i)
208 647952879 : dphirdd += _n1[i] * _d1[i] * MathUtils::pow(delta, _d1[i] - 1.0) * std::pow(tau, _t1[i]);
209 :
210 2591811516 : for (std::size_t i = 0; i < _n2.size(); ++i)
211 2499246819 : dphirdd += _n2[i] * std::exp(-MathUtils::pow(delta, _c2[i])) *
212 2499246819 : (MathUtils::pow(delta, _d2[i] - 1.0) * std::pow(tau, _t2[i]) *
213 2499246819 : (_d2[i] - _c2[i] * MathUtils::pow(delta, _c2[i])));
214 :
215 555388182 : for (std::size_t i = 0; i < _n3.size(); ++i)
216 462823485 : dphirdd += _n3[i] * MathUtils::pow(delta, _d3[i]) * MathUtils::pow(tau, _t3[i]) *
217 462823485 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _eps3[i]) -
218 462823485 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
219 462823485 : (_d3[i] / delta - 2.0 * _alpha3[i] * (delta - _eps3[i]));
220 :
221 370258788 : for (std::size_t i = 0; i < _n4.size(); ++i)
222 : {
223 277694091 : theta = 1.0 - tau + _A4[i] * std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]));
224 277694091 : Delta = Utility::pow<2>(theta) + _B4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i]);
225 277694091 : Psi = std::exp(-_C4[i] * Utility::pow<2>(delta - 1.0) - _D4[i] * Utility::pow<2>(tau - 1.0));
226 277694091 : dPsi_dd = -2.0 * _C4[i] * (delta - 1.0) * Psi;
227 277694091 : dDelta_dd = (delta - 1.0) *
228 277694091 : (_A4[i] * theta * 2.0 / _beta4[i] *
229 277694091 : std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]) - 1.0) +
230 277694091 : 2.0 * _B4[i] * _a4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i] - 1.0));
231 :
232 277694091 : dphirdd += _n4[i] * (std::pow(Delta, _b4[i]) * (Psi + delta * dPsi_dd) +
233 277694091 : _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 92564697 : return dphi0dd + dphirdd;
238 : }
239 :
240 : Real
241 2720813 : CO2FluidProperties::dalpha_dtau(Real delta, Real tau) const
242 : {
243 : // Derivative of the ideal gas component wrt tau
244 : Real sum0 = 0.0;
245 16324878 : for (std::size_t i = 0; i < _a0.size(); ++i)
246 13604065 : sum0 += _a0[i] * _theta0[i] * (1.0 / (1.0 - std::exp(-_theta0[i] * tau)) - 1.0);
247 :
248 2720813 : 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 21766504 : for (std::size_t i = 0; i < _n1.size(); ++i)
254 19045691 : dphirdt += _n1[i] * _t1[i] * MathUtils::pow(delta, _d1[i]) * std::pow(tau, _t1[i] - 1.0);
255 :
256 76182764 : for (std::size_t i = 0; i < _n2.size(); ++i)
257 73461951 : dphirdt += _n2[i] * _t2[i] * MathUtils::pow(delta, _d2[i]) * std::pow(tau, _t2[i] - 1.0) *
258 73461951 : std::exp(-MathUtils::pow(delta, _c2[i]));
259 :
260 16324878 : for (std::size_t i = 0; i < _n3.size(); ++i)
261 13604065 : dphirdt += _n3[i] * MathUtils::pow(delta, _d3[i]) * MathUtils::pow(tau, _t3[i]) *
262 13604065 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _eps3[i]) -
263 13604065 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
264 13604065 : (_t3[i] / tau - 2.0 * _beta3[i] * (tau - _gamma3[i]));
265 :
266 10883252 : for (std::size_t i = 0; i < _n4.size(); ++i)
267 : {
268 8162439 : theta = 1.0 - tau + _A4[i] * std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]));
269 8162439 : Delta = Utility::pow<2>(theta) + _B4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i]);
270 8162439 : Psi = std::exp(-_C4[i] * Utility::pow<2>(delta - 1.0) - _D4[i] * Utility::pow<2>(tau - 1.0));
271 8162439 : dDelta_dt = -2.0 * theta * _b4[i] * std::pow(Delta, _b4[i] - 1.0);
272 8162439 : dPsi_dt = -2.0 * _D4[i] * (tau - 1.0) * Psi;
273 :
274 8162439 : 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 2720813 : return dphi0dt + dphirdt;
279 : }
280 :
281 : Real
282 1800512 : CO2FluidProperties::d2alpha_ddelta2(Real delta, Real tau) const
283 : {
284 : // Second derivative of the ideal gas component wrt gamma
285 1800512 : 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 14404096 : for (std::size_t i = 0; i < _n1.size(); ++i)
292 12603584 : d2phirdd2 += _n1[i] * _d1[i] * (_d1[i] - 1.0) * MathUtils::pow(delta, _d1[i] - 2) *
293 12603584 : std::pow(tau, _t1[i]);
294 :
295 50414336 : for (std::size_t i = 0; i < _n2.size(); ++i)
296 48613824 : d2phirdd2 += _n2[i] * std::exp(-MathUtils::pow(delta, _c2[i])) *
297 48613824 : MathUtils::pow(delta, _d2[i] - 2) * std::pow(tau, _t2[i]) *
298 48613824 : ((_d2[i] - _c2[i] * MathUtils::pow(delta, _c2[i])) *
299 48613824 : (_d2[i] - 1.0 - _c2[i] * MathUtils::pow(delta, _c2[i])) -
300 48613824 : _c2[i] * _c2[i] * MathUtils::pow(delta, _c2[i]));
301 :
302 10803072 : for (std::size_t i = 0; i < _n3.size(); ++i)
303 9002560 : d2phirdd2 +=
304 9002560 : _n3[i] * MathUtils::pow(tau, _t3[i]) *
305 9002560 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _eps3[i]) -
306 9002560 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
307 9002560 : (-2.0 * _alpha3[i] * MathUtils::pow(delta, _d3[i]) +
308 9002560 : 4.0 * _alpha3[i] * _alpha3[i] * MathUtils::pow(delta, _d3[i]) *
309 9002560 : Utility::pow<2>(delta - _eps3[i]) -
310 9002560 : 4.0 * _d3[i] * _alpha3[i] * MathUtils::pow(delta, _d3[i] - 1.0) * (delta - _eps3[i]) +
311 9002560 : _d3[i] * (_d3[i] - 1.0) * MathUtils::pow(delta, _d3[i] - 2.0));
312 :
313 7202048 : for (std::size_t i = 0; i < _n4.size(); ++i)
314 : {
315 5401536 : theta = 1.0 - tau + _A4[i] * std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]));
316 5401536 : Delta = Utility::pow<2>(theta) + _B4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i]);
317 5401536 : Psi = std::exp(-_C4[i] * Utility::pow<2>(delta - 1.0) - _D4[i] * Utility::pow<2>(tau - 1.0));
318 5401536 : dPsi_dd = -2.0 * _C4[i] * (delta - 1.0) * Psi;
319 5401536 : dDelta_dd = (delta - 1.0) *
320 5401536 : (_A4[i] * theta * 2.0 / _beta4[i] *
321 5401536 : std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]) - 1.0) +
322 5401536 : 2.0 * _B4[i] * _a4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i] - 1.0));
323 5401536 : d2Psi_dd2 = 3.0 * _D4[i] * Psi * (2.0 * _C4[i] * Utility::pow<2>(delta - 1.0) - 1.0);
324 5401536 : d2Delta_dd2 = 1.0 / (delta - 1.0) * dDelta_dd +
325 5401536 : (delta - 1.0) * (delta - 1.0) *
326 5401536 : (4.0 * _B4[i] * _a4[i] * (_a4[i] - 1.0) *
327 5401536 : std::pow(Utility::pow<2>(delta - 1.0), _a4[i] - 2.0) +
328 5401536 : 2.0 * _A4[i] * _A4[i] *
329 5401536 : Utility::pow<2>(std::pow(Utility::pow<2>(delta - 1.0),
330 5401536 : 1.0 / (2.0 * _beta4[i]) - 1.0)) /
331 5401536 : _beta4[i] / _beta4[i] +
332 5401536 : (4.0 / _beta4[i]) * _A4[i] * theta * (1.0 / (2.0 * _beta4[i]) - 1.0) *
333 5401536 : std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]) - 2.0));
334 5401536 : d2phirdd2 +=
335 5401536 : _n4[i] *
336 5401536 : (std::pow(Delta, _b4[i]) * (2.0 * dPsi_dd + delta * d2Psi_dd2) +
337 5401536 : 2.0 * _b4[i] * std::pow(Delta, _b4[i] - 1.0) * dDelta_dd * (Psi + delta * dPsi_dd) +
338 5401536 : _b4[i] *
339 5401536 : (std::pow(Delta, _b4[i] - 1.0) * d2Delta_dd2 +
340 5401536 : (_b4[i] - 1.0) * std::pow(Delta, _b4[i] - 2.0) * Utility::pow<2>(dDelta_dd)) *
341 5401536 : delta * Psi);
342 : }
343 : // The second derivative of the free energy wrt delta is the sum of these components
344 1800512 : return d2phi0dd2 + d2phirdd2;
345 : }
346 :
347 : Real
348 2700759 : 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 16204554 : for (std::size_t i = 0; i < _a0.size(); ++i)
353 13503795 : sum0 += _a0[i] * _theta0[i] * _theta0[i] * std::exp(-_theta0[i] * tau) /
354 13503795 : Utility::pow<2>(1.0 - std::exp(-_theta0[i] * tau));
355 :
356 2700759 : 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 21606072 : for (std::size_t i = 0; i < _n1.size(); ++i)
363 18905313 : d2phirdt2 += _n1[i] * _t1[i] * (_t1[i] - 1.0) * MathUtils::pow(delta, _d1[i]) *
364 18905313 : std::pow(tau, _t1[i] - 2.0);
365 :
366 75621252 : for (std::size_t i = 0; i < _n2.size(); ++i)
367 72920493 : d2phirdt2 += _n2[i] * _t2[i] * (_t2[i] - 1.0) * MathUtils::pow(delta, _d2[i]) *
368 72920493 : std::exp(-MathUtils::pow(delta, _c2[i])) * std::pow(tau, _t2[i] - 2.0);
369 :
370 16204554 : for (std::size_t i = 0; i < _n3.size(); ++i)
371 13503795 : d2phirdt2 += _n3[i] * MathUtils::pow(delta, _d3[i]) * MathUtils::pow(tau, _t3[i]) *
372 13503795 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _eps3[i]) -
373 13503795 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
374 13503795 : (Utility::pow<2>(_t3[i] / tau - 2.0 * _beta3[i] * (tau - _gamma3[i])) -
375 13503795 : _t3[i] / tau / tau - 2.0 * _beta3[i]);
376 :
377 10803036 : for (std::size_t i = 0; i < _n4.size(); ++i)
378 : {
379 8102277 : theta = 1.0 - tau + _A4[i] * std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]));
380 8102277 : Delta = Utility::pow<2>(theta) + _B4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i]);
381 8102277 : Psi = std::exp(-_C4[i] * Utility::pow<2>(delta - 1.0) - _D4[i] * Utility::pow<2>(tau - 1.0));
382 8102277 : dDelta_dt = -2.0 * theta * _b4[i] * std::pow(Delta, _b4[i] - 1.0);
383 8102277 : d2Delta_dt2 = 2.0 * _b4[i] * std::pow(Delta, _b4[i] - 1.0) +
384 8102277 : 4.0 * theta * theta * _b4[i] * (_b4[i] - 1.0) * std::pow(Delta, _b4[i] - 2.0);
385 8102277 : dPsi_dt = -2.0 * _D4[i] * (tau - 1.0) * Psi;
386 8102277 : d2Psi_dt2 = 2.0 * _D4[i] * (2.0 * _D4[i] * (tau - 1.0) * (tau - 1.0) - 1.0) * Psi;
387 8102277 : d2phirdt2 +=
388 8102277 : _n4[i] * delta *
389 8102277 : (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 2700759 : return d2phi0dt2 + d2phirdt2;
394 : }
395 :
396 : Real
397 1800512 : 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 14404096 : for (std::size_t i = 0; i < _n1.size(); ++i)
404 12603584 : d2phirddt += _n1[i] * _d1[i] * _t1[i] * MathUtils::pow(delta, _d1[i] - 1.0) *
405 12603584 : std::pow(tau, _t1[i] - 1.0);
406 :
407 50414336 : for (std::size_t i = 0; i < _n2.size(); ++i)
408 48613824 : d2phirddt += _n2[i] * std::exp(-MathUtils::pow(delta, _c2[i])) *
409 48613824 : (MathUtils::pow(delta, _d2[i] - 1.0) * _t2[i] * std::pow(tau, _t2[i] - 1.0) *
410 48613824 : (_d2[i] - _c2[i] * MathUtils::pow(delta, _c2[i])));
411 :
412 10803072 : for (std::size_t i = 0; i < _n3.size(); ++i)
413 9002560 : d2phirddt += _n3[i] * MathUtils::pow(delta, _d3[i]) * MathUtils::pow(tau, _t3[i]) *
414 9002560 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _eps3[i]) -
415 9002560 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
416 9002560 : (_d3[i] / delta - 2.0 * _alpha3[i] * (delta - _eps3[i])) *
417 9002560 : (_t3[i] / tau - 2.0 * _beta3[i] * (tau - _gamma3[i]));
418 :
419 7202048 : for (std::size_t i = 0; i < _n4.size(); ++i)
420 : {
421 5401536 : theta = 1.0 - tau + _A4[i] * std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]));
422 5401536 : Delta = Utility::pow<2>(theta) + _B4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i]);
423 5401536 : Psi = std::exp(-_C4[i] * Utility::pow<2>(delta - 1.0) - _D4[i] * Utility::pow<2>(tau - 1.0));
424 5401536 : dPsi_dd = -2.0 * _C4[i] * (delta - 1.0) * Psi;
425 5401536 : dPsi_dt = -2.0 * _D4[i] * (tau - 1.0) * Psi;
426 5401536 : d2Psi_ddt = 4.0 * _C4[i] * _D4[i] * (delta - 1.0) * (tau - 1.0) * Psi;
427 5401536 : dDelta_dd = (delta - 1.0) *
428 5401536 : (_A4[i] * theta * 2.0 / _beta4[i] *
429 5401536 : std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]) - 1.0) +
430 5401536 : 2.0 * _B4[i] * _a4[i] * std::pow(Utility::pow<2>(delta - 1.0), _a4[i] - 1.0));
431 5401536 : dDelta_dt = -2.0 * theta * _b4[i] * std::pow(Delta, _b4[i] - 1.0);
432 5401536 : d2Delta_ddt = -2.0 * _A4[i] * _b4[i] / _beta4[i] * std::pow(Delta, _b4[i] - 1.0) *
433 5401536 : (delta - 1.0) *
434 5401536 : std::pow(Utility::pow<2>(delta - 1.0), 1.0 / (2.0 * _beta4[i]) - 1.0) -
435 5401536 : 2.0 * theta * _b4[i] * (_b4[i] - 1.0) * std::pow(Delta, _b4[i] - 2.0) * dDelta_dd;
436 :
437 5401536 : d2phirddt += _n4[i] * (std::pow(Delta, _b4[i]) * (dPsi_dt + delta * d2Psi_ddt) +
438 5401536 : delta * _b4[i] * std::pow(Delta, _b4[i] - 1.0) * dDelta_dd * dPsi_dt +
439 5401536 : dDelta_dt * (Psi + delta * dPsi_dd) + d2Delta_ddt * delta * Psi);
440 : }
441 :
442 1800512 : return d2phirddt;
443 : }
444 :
445 : Real
446 89864966 : CO2FluidProperties::p_from_rho_T(Real density, Real temperature) const
447 : {
448 : // Check that the input parameters are within the region of validity
449 89864966 : 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 89864966 : if (temperature > _triple_point_temperature && temperature < _critical_temperature)
455 : {
456 7491475 : Real gas_density = saturatedVaporDensity(temperature);
457 7491475 : Real liquid_density = saturatedLiquidDensity(temperature);
458 :
459 7491475 : if (density < gas_density || density > liquid_density)
460 7480419 : pressure = HelmholtzFluidProperties::p_from_rho_T(density, temperature);
461 : else
462 11056 : pressure = vaporPressure(temperature);
463 : }
464 : else
465 82373491 : pressure = HelmholtzFluidProperties::p_from_rho_T(density, temperature);
466 :
467 89864966 : return pressure;
468 : }
469 :
470 : Real
471 8142407 : CO2FluidProperties::rho_from_p_T(Real pressure, Real temperature) const
472 : {
473 : // Check that the input parameters are within the region of validity
474 8142407 : 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 8142407 : if (((temperature > _triple_point_temperature) && (pressure > meltingPressure(temperature))) ||
479 8142407 : ((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 8142407 : Real lower_density = 100.0;
486 8142407 : 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 89864966 : auto pressure_diff = [&pressure, &temperature, this](Real x)
491 89864966 : { return p_from_rho_T(x, temperature) - pressure; };
492 :
493 8142407 : BrentsMethod::bracket(pressure_diff, lower_density, upper_density);
494 8142407 : density = BrentsMethod::root(pressure_diff, lower_density, upper_density);
495 :
496 8142407 : 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 910274 : CO2FluidProperties::mu_from_p_T(Real pressure, Real temperature) const
508 : {
509 910274 : Real rho = rho_from_p_T(pressure, temperature);
510 910274 : 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 910283 : CO2FluidProperties::mu_from_rho_T(Real density, Real temperature) const
527 : {
528 : // Check that the input parameters are within the region of validity
529 910283 : 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 910283 : Real Tstar = temperature / 251.196;
533 :
534 : // Viscosity in the zero-density limit
535 : Real sum = 0.0;
536 :
537 5461698 : for (std::size_t i = 0; i < _mu_a.size(); ++i)
538 4551415 : sum += _mu_a[i] * MathUtils::pow(std::log(Tstar), i);
539 :
540 910283 : Real mu0 = 1.00697 * std::sqrt(temperature) / std::exp(sum);
541 :
542 : // Excess viscosity due to finite density
543 910283 : Real mue = _mu_d[0] * density + _mu_d[1] * Utility::pow<2>(density) +
544 910283 : _mu_d[2] * Utility::pow<6>(density) / Utility::pow<3>(Tstar) +
545 910283 : _mu_d[3] * Utility::pow<8>(density) + _mu_d[4] * Utility::pow<8>(density) / Tstar;
546 :
547 910283 : 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 900272 : CO2FluidProperties::k_from_p_T(Real pressure, Real temperature) const
642 : {
643 : // Require density first
644 900272 : Real density = rho_from_p_T(pressure, temperature);
645 900272 : return k_from_rho_T(density, temperature);
646 : }
647 :
648 : void
649 1 : CO2FluidProperties::k_from_p_T(
650 : Real pressure, Real temperature, Real & k, Real & dk_dp, Real & dk_dT) const
651 : {
652 1 : 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 1 : const Real peps = pressure * eps;
657 1 : const Real Teps = temperature * eps;
658 :
659 1 : dk_dp = (this->k_from_p_T(pressure + peps, temperature) - k) / peps;
660 1 : dk_dT = (this->k_from_p_T(pressure, temperature + Teps) - k) / Teps;
661 1 : }
662 :
663 : Real
664 900275 : 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 900275 : 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 900275 : Real Tr = temperature / _critical_temperature;
673 900275 : Real rhor = density / _critical_density;
674 :
675 : Real sum1 = 0.0;
676 3601100 : for (std::size_t i = 0; i < _k_n1.size(); ++i)
677 2700825 : sum1 += _k_n1[i] * std::pow(Tr, _k_g1[i]) * MathUtils::pow(rhor, _k_h1[i]);
678 :
679 : Real sum2 = 0.0;
680 7202200 : for (std::size_t i = 0; i < _k_n2.size(); ++i)
681 6301925 : 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 900275 : 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 900275 : rhor *
688 900275 : std::exp(-std::pow(rhor, _k_a[0]) / _k_a[0] - Utility::pow<2>(_k_a[1] * (Tr - 1.0)) -
689 900275 : Utility::pow<2>(_k_a[2] * (rhor - 1.0))) /
690 900275 : std::pow(std::pow(Utility::pow<2>(
691 900275 : 1.0 - 1.0 / Tr +
692 900275 : _k_a[3] * std::pow(Utility::pow<2>(rhor - 1.0), 1.0 / (2.0 * _k_a[4]))),
693 : _k_a[5]) +
694 900275 : std::pow(Utility::pow<2>(_k_a[6] * (rhor - alpha)), _k_a[7]),
695 900275 : _k_a[8]);
696 :
697 900275 : return 4.81384 * (sum1 + std::exp(-5.0 * rhor * rhor) * sum2 + 0.775547504 * lambdac) / 1000.0;
698 : }
|
