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 "MethaneFluidProperties.h"
11 : #include "MathUtils.h"
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("FluidPropertiesApp", MethaneFluidProperties);
15 :
16 : InputParameters
17 37 : MethaneFluidProperties::validParams()
18 : {
19 37 : InputParameters params = HelmholtzFluidProperties::validParams();
20 37 : params.addClassDescription("Fluid properties for methane (CH4)");
21 37 : return params;
22 0 : }
23 :
24 19 : MethaneFluidProperties::MethaneFluidProperties(const InputParameters & parameters)
25 : : HelmholtzFluidProperties(parameters),
26 19 : _Mch4(16.0425e-3),
27 19 : _p_critical(4.5992e6),
28 19 : _T_critical(190.564),
29 19 : _rho_critical(162.66),
30 19 : _p_triple(1.169e4),
31 19 : _T_triple(90.6941)
32 : {
33 19 : }
34 :
35 19 : MethaneFluidProperties::~MethaneFluidProperties() {}
36 :
37 : std::string
38 1 : MethaneFluidProperties::fluidName() const
39 : {
40 1 : return "methane";
41 : }
42 :
43 : Real
44 3664 : MethaneFluidProperties::molarMass() const
45 : {
46 3664 : return _Mch4;
47 : }
48 :
49 : Real
50 1 : MethaneFluidProperties::criticalPressure() const
51 : {
52 1 : return _p_critical;
53 : }
54 :
55 : Real
56 3664 : MethaneFluidProperties::criticalTemperature() const
57 : {
58 3664 : return _T_critical;
59 : }
60 :
61 : Real
62 3664 : MethaneFluidProperties::criticalDensity() const
63 : {
64 3664 : return _rho_critical;
65 : }
66 :
67 : Real
68 1 : MethaneFluidProperties::triplePointPressure() const
69 : {
70 1 : return _p_triple;
71 : }
72 :
73 : Real
74 1 : MethaneFluidProperties::triplePointTemperature() const
75 : {
76 1 : return _T_triple;
77 : }
78 :
79 : Real
80 58 : MethaneFluidProperties::mu_from_p_T(Real /*pressure*/, Real temperature) const
81 : {
82 : // Check the temperature is in the range of validity (200 K <= T <= 1000 K)
83 58 : if (temperature <= 200.0 || temperature >= 1000.0)
84 0 : mooseException(
85 : "Temperature ", temperature, "K out of range (200K, 1000K) in ", name(), ": mu_from_p_T()");
86 :
87 : Real viscosity = 0.0;
88 406 : for (std::size_t i = 0; i < _a.size(); ++i)
89 348 : viscosity += _a[i] * MathUtils::pow(temperature, i);
90 :
91 58 : return viscosity * 1.e-6;
92 : }
93 :
94 : void
95 3 : MethaneFluidProperties::mu_from_p_T(
96 : Real pressure, Real temperature, Real & mu, Real & dmu_dp, Real & dmu_dT) const
97 : {
98 :
99 3 : mu = this->mu_from_p_T(pressure, temperature);
100 3 : dmu_dp = 0.0;
101 :
102 : Real dmudt = 0.0;
103 21 : for (std::size_t i = 0; i < _a.size(); ++i)
104 18 : dmudt += i * _a[i] * MathUtils::pow(temperature, i) / temperature;
105 3 : dmu_dT = dmudt * 1.e-6;
106 3 : }
107 :
108 : Real
109 54 : MethaneFluidProperties::k_from_p_T(Real /*pressure*/, Real temperature) const
110 : {
111 : // Check the temperature is in the range of validity (200 K <= T <= 1000 K)
112 54 : if (temperature <= 200.0 || temperature >= 1000.0)
113 0 : mooseException(
114 : "Temperature ", temperature, "K out of range (200K, 1000K) in ", name(), ": k()");
115 :
116 : Real kt = 0.0;
117 432 : for (std::size_t i = 0; i < _b.size(); ++i)
118 378 : kt += _b[i] * MathUtils::pow(temperature, i);
119 :
120 54 : return kt;
121 : }
122 :
123 : void
124 1 : MethaneFluidProperties::k_from_p_T(
125 : Real /*pressure*/, Real temperature, Real & k, Real & dk_dp, Real & dk_dT) const
126 : {
127 : // Check the temperature is in the range of validity (200 K <= T <= 1000 K)
128 1 : if (temperature <= 200.0 || temperature >= 1000.0)
129 0 : mooseException(
130 : "Temperature ", temperature, "K out of range (200K, 1000K) in ", name(), ": k()");
131 :
132 : Real kt = 0.0, dkt_dT = 0.0;
133 :
134 8 : for (std::size_t i = 0; i < _b.size(); ++i)
135 7 : kt += _b[i] * MathUtils::pow(temperature, i);
136 :
137 7 : for (std::size_t i = 1; i < _b.size(); ++i)
138 6 : dkt_dT += i * _b[i] * MathUtils::pow(temperature, i) / temperature;
139 :
140 1 : k = kt;
141 1 : dk_dp = 0.0;
142 1 : dk_dT = dkt_dT;
143 1 : }
144 :
145 : Real
146 3 : MethaneFluidProperties::vaporPressure(Real temperature) const
147 : {
148 3 : if (temperature < _T_triple || temperature > _T_critical)
149 0 : mooseException("Temperature is out of range in ", name(), ": vaporPressure()");
150 :
151 3 : const Real Tr = temperature / _T_critical;
152 3 : const Real theta = 1.0 - Tr;
153 :
154 3 : const Real logpressure = (-6.036219 * theta + 1.409353 * std::pow(theta, 1.5) -
155 3 : 0.4945199 * Utility::pow<2>(theta) - 1.443048 * std::pow(theta, 4.5)) /
156 3 : Tr;
157 :
158 3 : return _p_critical * std::exp(logpressure);
159 : }
160 :
161 : void
162 0 : MethaneFluidProperties::vaporPressure(Real, Real &, Real &) const
163 : {
164 0 : mooseError("vaporPressure() is not implemented");
165 : }
166 :
167 : std::vector<Real>
168 1 : MethaneFluidProperties::henryCoefficients() const
169 : {
170 1 : return {-10.44708, 4.66491, 12.12986};
171 : }
172 :
173 : Real
174 3 : MethaneFluidProperties::saturatedLiquidDensity(Real temperature) const
175 : {
176 3 : if (temperature < _T_triple || temperature > _T_critical)
177 0 : mooseException(
178 : "Temperature ", temperature, " is out of range in ", name(), ": saturatedLiquidDensity()");
179 :
180 3 : const Real Tr = temperature / _T_critical;
181 3 : const Real theta = 1.0 - Tr;
182 :
183 3 : const Real logdensity = 1.9906389 * std::pow(theta, 0.354) - 0.78756197 * std::sqrt(theta) +
184 3 : 0.036976723 * std::pow(theta, 2.5);
185 :
186 3 : return _rho_critical * std::exp(logdensity);
187 : }
188 :
189 : Real
190 3 : MethaneFluidProperties::saturatedVaporDensity(Real temperature) const
191 : {
192 3 : if (temperature < _T_triple || temperature > _T_critical)
193 0 : mooseException(
194 : "Temperature ", temperature, " is out of range in ", name(), ": saturatedVaporDensity()");
195 :
196 3 : const Real Tr = temperature / _T_critical;
197 3 : const Real theta = 1.0 - Tr;
198 :
199 : const Real logdensity =
200 3 : -1.880284 * std::pow(theta, 0.354) - 2.8526531 * std::pow(theta, 5.0 / 6.0) -
201 3 : 3.000648 * std::pow(theta, 1.5) - 5.251169 * std::pow(theta, 2.5) -
202 3 : 13.191859 * std::pow(theta, 25.0 / 6.0) - 37.553961 * std::pow(theta, 47.0 / 6.0);
203 :
204 3 : return _rho_critical * std::exp(logdensity);
205 : }
206 :
207 : Real
208 49 : MethaneFluidProperties::alpha(Real delta, Real tau) const
209 : {
210 : // Ideal gas component of the Helmholtz free energy
211 49 : Real alpha0 = std::log(delta) + 9.91243972 - 6.33270087 * tau + 3.0016 * std::log(tau);
212 :
213 294 : for (std::size_t i = 0; i < _a0.size(); ++i)
214 245 : alpha0 += _a0[i] * std::log(1.0 - std::exp(-_b0[i] * tau));
215 :
216 : // Residual component of the Helmholtz free energy
217 : Real alphar = 0.0;
218 :
219 686 : for (std::size_t i = 0; i < _t1.size(); ++i)
220 637 : alphar += _N1[i] * MathUtils::pow(delta, _d1[i]) * std::pow(tau, _t1[i]);
221 :
222 1176 : for (std::size_t i = 0; i < _t2.size(); ++i)
223 1127 : alphar += _N2[i] * MathUtils::pow(delta, _d2[i]) * std::pow(tau, _t2[i]) *
224 1127 : std::exp(-MathUtils::pow(delta, _c2[i]));
225 :
226 245 : for (std::size_t i = 0; i < _t3.size(); ++i)
227 196 : alphar += _N3[i] * MathUtils::pow(delta, _d3[i]) * std::pow(tau, _t3[i]) *
228 196 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _D3[i]) -
229 196 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i]));
230 :
231 : // The Helmholtz free energy is the sum of these two
232 49 : return alpha0 + alphar;
233 : }
234 :
235 : Real
236 3508 : MethaneFluidProperties::dalpha_ddelta(Real delta, Real tau) const
237 : {
238 : // Ideal gas component of the Helmholtz free energy
239 3508 : Real dalpha0 = 1.0 / delta;
240 :
241 : // Residual component of the Helmholtz free energy
242 : Real dalphar = 0.0;
243 :
244 49112 : for (std::size_t i = 0; i < _t1.size(); ++i)
245 45604 : dalphar += _N1[i] * _d1[i] * MathUtils::pow(delta, _d1[i]) * std::pow(tau, _t1[i]);
246 :
247 84192 : for (std::size_t i = 0; i < _t2.size(); ++i)
248 80684 : dalphar += _N2[i] * MathUtils::pow(delta, _d2[i]) * std::pow(tau, _t2[i]) *
249 80684 : std::exp(-MathUtils::pow(delta, _c2[i])) *
250 80684 : (_d2[i] - _c2[i] * MathUtils::pow(delta, _c2[i]));
251 :
252 17540 : for (std::size_t i = 0; i < _t3.size(); ++i)
253 14032 : dalphar += _N3[i] * MathUtils::pow(delta, _d3[i]) * std::pow(tau, _t3[i]) *
254 14032 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _D3[i]) -
255 14032 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
256 14032 : (_d3[i] - delta * (2.0 * _alpha3[i] * (delta - _D3[i])));
257 :
258 : // The Helmholtz free energy is the sum of these two
259 3508 : return dalpha0 + dalphar / delta;
260 : }
261 :
262 : Real
263 161 : MethaneFluidProperties::dalpha_dtau(Real delta, Real tau) const
264 : {
265 : // Ideal gas component of the Helmholtz free energy
266 161 : Real dalpha0 = -6.33270087 + 3.0016 / tau;
267 :
268 966 : for (std::size_t i = 0; i < _a0.size(); ++i)
269 805 : dalpha0 += _a0[i] * _b0[i] * (1.0 / (1.0 - std::exp(-_b0[i] * tau)) - 1.0);
270 :
271 : // Residual component of the Helmholtz free energy
272 : Real dalphar = 0.0;
273 :
274 2254 : for (std::size_t i = 0; i < _t1.size(); ++i)
275 2093 : dalphar += _N1[i] * _t1[i] * MathUtils::pow(delta, _d1[i]) * std::pow(tau, _t1[i]);
276 :
277 3864 : for (std::size_t i = 0; i < _t2.size(); ++i)
278 3703 : dalphar += _N2[i] * _t2[i] * MathUtils::pow(delta, _d2[i]) * std::pow(tau, _t2[i]) *
279 3703 : std::exp(-MathUtils::pow(delta, _c2[i]));
280 :
281 805 : for (std::size_t i = 0; i < _t3.size(); ++i)
282 644 : dalphar += _N3[i] * MathUtils::pow(delta, _d3[i]) * std::pow(tau, _t3[i]) *
283 644 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _D3[i]) -
284 644 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
285 644 : (_t3[i] + tau * (2.0 * _beta3[i] * (tau - _gamma3[i])));
286 :
287 : // The Helmholtz free energy is the sum of these two
288 161 : return dalpha0 + dalphar / tau;
289 : }
290 :
291 : Real
292 106 : MethaneFluidProperties::d2alpha_ddelta2(Real delta, Real tau) const
293 : {
294 : // Ideal gas component of the Helmholtz free energy
295 106 : Real dalpha0 = -1.0 / delta / delta;
296 :
297 : // Residual component of the Helmholtz free energy
298 : Real dalphar = 0.0;
299 :
300 1484 : for (std::size_t i = 0; i < _t1.size(); ++i)
301 1378 : dalphar +=
302 1378 : _N1[i] * _d1[i] * (_d1[i] - 1.0) * MathUtils::pow(delta, _d1[i]) * std::pow(tau, _t1[i]);
303 :
304 2544 : for (std::size_t i = 0; i < _t2.size(); ++i)
305 2438 : dalphar += _N2[i] * MathUtils::pow(delta, _d2[i]) * std::pow(tau, _t2[i]) *
306 2438 : std::exp(-MathUtils::pow(delta, _c2[i])) *
307 2438 : ((_d2[i] - _c2[i] * MathUtils::pow(delta, _c2[i])) *
308 2438 : (_d2[i] - 1.0 - _c2[i] * MathUtils::pow(delta, _c2[i])) -
309 2438 : _c2[i] * _c2[i] * MathUtils::pow(delta, _c2[i]));
310 :
311 530 : for (std::size_t i = 0; i < _t3.size(); ++i)
312 424 : dalphar += _N3[i] * std::pow(tau, _t3[i]) *
313 424 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _D3[i]) -
314 424 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
315 424 : (-2.0 * _alpha3[i] * MathUtils::pow(delta, _d3[i] + 2) +
316 424 : 4.0 * _alpha3[i] * _alpha3[i] * MathUtils::pow(delta, _d3[i] + 2) *
317 424 : MathUtils::pow(delta - _D3[i], 2) -
318 424 : 4.0 * _d3[i] * _alpha3[i] * MathUtils::pow(delta, _d3[i] + 1) * (delta - _D3[i]) +
319 424 : _d3[i] * (_d3[i] - 1) * MathUtils::pow(delta, _d3[i]));
320 :
321 : // The Helmholtz free energy is the sum of these two
322 106 : return dalpha0 + dalphar / delta / delta;
323 : }
324 :
325 : Real
326 151 : MethaneFluidProperties::d2alpha_dtau2(Real delta, Real tau) const
327 : {
328 : // Ideal gas component of the Helmholtz free energy
329 151 : Real dalpha0 = -3.0016 / tau / tau;
330 :
331 906 : for (std::size_t i = 0; i < _a0.size(); ++i)
332 : {
333 755 : Real exptau = std::exp(-_b0[i] * tau);
334 755 : dalpha0 -= _a0[i] * (_b0[i] * _b0[i] * exptau / (1.0 - exptau) / (1.0 - exptau));
335 : }
336 :
337 : // Residual component of the Helmholtz free energy
338 : Real dalphar = 0.0;
339 :
340 2114 : for (std::size_t i = 0; i < _t1.size(); ++i)
341 1963 : dalphar +=
342 1963 : _N1[i] * _t1[i] * (_t1[i] - 1.0) * MathUtils::pow(delta, _d1[i]) * std::pow(tau, _t1[i]);
343 :
344 3624 : for (std::size_t i = 0; i < _t2.size(); ++i)
345 3473 : dalphar += _N2[i] * _t2[i] * (_t2[i] - 1.0) * MathUtils::pow(delta, _d2[i]) *
346 3473 : std::pow(tau, _t2[i]) * std::exp(-MathUtils::pow(delta, _c2[i]));
347 :
348 755 : for (std::size_t i = 0; i < _t3.size(); ++i)
349 604 : dalphar += _N3[i] * MathUtils::pow(delta, _d3[i]) * std::pow(tau, _t3[i]) *
350 604 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _D3[i]) -
351 604 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
352 604 : (tau * _t3[i] - _beta3[i] * tau * tau * MathUtils::pow(tau - _gamma3[i], 2) -
353 604 : _t3[i] - 2.0 * tau * tau * _beta3[i]);
354 :
355 : // The Helmholtz free energy is the sum of these two
356 151 : return dalpha0 + dalphar / tau / tau;
357 : }
358 :
359 : Real
360 106 : MethaneFluidProperties::d2alpha_ddeltatau(Real delta, Real tau) const
361 : {
362 : // Residual component of the Helmholtz free energy
363 : Real dalphar = 0.0;
364 :
365 1484 : for (std::size_t i = 0; i < _t1.size(); ++i)
366 1378 : dalphar += _N1[i] * _d1[i] * _t1[i] * std::pow(delta, _d1[i]) * std::pow(tau, _t1[i]);
367 :
368 2544 : for (std::size_t i = 0; i < _t2.size(); ++i)
369 2438 : dalphar += _N2[i] * _t2[i] * std::pow(delta, _d2[i]) * std::pow(tau, _t2[i]) *
370 2438 : std::exp(-MathUtils::pow(delta, _c2[i])) *
371 2438 : (_d2[i] - _c2[i] * MathUtils::pow(delta, _c2[i]));
372 :
373 530 : for (std::size_t i = 0; i < _t3.size(); ++i)
374 424 : dalphar += _N3[i] * std::pow(delta, _d3[i]) * std::pow(tau, _t3[i]) *
375 424 : std::exp(-_alpha3[i] * Utility::pow<2>(delta - _D3[i]) -
376 424 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
377 424 : (_d3[i] - 2.0 * _alpha3[i] * delta * (delta - _D3[i]) *
378 424 : (_t3[i] - 2.0 * _beta3[i] * tau * (tau - _gamma3[i])));
379 :
380 : // The Helmholtz free energy is the sum of these two
381 106 : return dalphar / delta / tau;
382 : }
|
