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 "HydrogenFluidProperties.h"
11 : #include "Conversion.h"
12 : #include "MathUtils.h"
13 : #include "libmesh/utility.h"
14 :
15 : registerMooseObject("FluidPropertiesApp", HydrogenFluidProperties);
16 :
17 : InputParameters
18 22 : HydrogenFluidProperties::validParams()
19 : {
20 22 : InputParameters params = HelmholtzFluidProperties::validParams();
21 22 : params.addClassDescription("Fluid properties for Hydrogen (H2)");
22 22 : return params;
23 0 : }
24 :
25 11 : HydrogenFluidProperties::HydrogenFluidProperties(const InputParameters & parameters)
26 : : HelmholtzFluidProperties(parameters),
27 11 : _Mh2(2.01588e-3),
28 11 : _p_critical(1.315e6),
29 11 : _T_critical(33.19),
30 11 : _rho_molar_critical(15.508),
31 11 : _rho_critical(1000.0 * _rho_molar_critical * _Mh2),
32 11 : _p_triple(7.7e3),
33 11 : _T_triple(13.952)
34 : {
35 11 : }
36 :
37 : std::string
38 1 : HydrogenFluidProperties::fluidName() const
39 : {
40 1 : return "hydrogen";
41 : }
42 :
43 : Real
44 819 : HydrogenFluidProperties::molarMass() const
45 : {
46 819 : return _Mh2;
47 : }
48 :
49 : Real
50 1 : HydrogenFluidProperties::criticalPressure() const
51 : {
52 1 : return _p_critical;
53 : }
54 :
55 : Real
56 819 : HydrogenFluidProperties::criticalTemperature() const
57 : {
58 819 : return _T_critical;
59 : }
60 :
61 : Real
62 819 : HydrogenFluidProperties::criticalDensity() const
63 : {
64 819 : return _rho_critical;
65 : }
66 :
67 : Real
68 1 : HydrogenFluidProperties::triplePointPressure() const
69 : {
70 1 : return _p_triple;
71 : }
72 :
73 : Real
74 1 : HydrogenFluidProperties::triplePointTemperature() const
75 : {
76 1 : return _T_triple;
77 : }
78 :
79 : Real
80 21 : HydrogenFluidProperties::mu_from_rho_T(Real density, Real temperature) const
81 : {
82 : // Scaled variables
83 21 : const Real Tstar = temperature / 30.41;
84 21 : const Real logTstar = std::log(Tstar);
85 21 : const Real Tr = temperature / _T_critical;
86 21 : const Real rhor = density / 90.5;
87 :
88 : // Ideal gas component
89 : Real sum = 0.0;
90 126 : for (std::size_t i = 0; i < _amu.size(); ++i)
91 105 : sum += _amu[i] * MathUtils::pow(logTstar, i);
92 :
93 : const Real mu0 =
94 21 : 0.021357 * std::sqrt(1000.0 * _Mh2 * temperature) / (0.297 * 0.297 * std::exp(sum));
95 :
96 : // The excess contribution due to density
97 : Real sumr = 0.0;
98 168 : for (std::size_t i = 0; i < _bmu.size(); ++i)
99 147 : sumr += _bmu[i];
100 :
101 21 : const Real mu1 = MathUtils::pow(0.297, 3) * sumr * mu0 / Tstar;
102 :
103 : // The viscosity is then
104 : const Real mu =
105 21 : mu0 + mu1 * density +
106 21 : _cmu[0] * rhor * rhor *
107 21 : std::exp(_cmu[1] * Tr + _cmu[2] / Tr + _cmu[3] * rhor * rhor / (_cmu[4] + Tr) +
108 21 : _cmu[5] * MathUtils::pow(rhor, 6));
109 :
110 21 : return mu * 1.0e-6;
111 : }
112 :
113 : void
114 7 : HydrogenFluidProperties::mu_from_rho_T(Real density,
115 : Real temperature,
116 : Real ddensity_dT,
117 : Real & mu,
118 : Real & dmu_drho,
119 : Real & dmu_dT) const
120 : {
121 : // Scaled variables
122 7 : const Real Tstar = temperature / 30.41;
123 7 : const Real logTstar = std::log(Tstar);
124 7 : const Real Tr = temperature / _T_critical;
125 7 : const Real rhor = density / 90.5;
126 : const Real drhor_drho = 1.0 / 90.5;
127 7 : const Real dTr_dT = 1.0 / _T_critical;
128 :
129 : // The dilute gas component
130 : Real sum = 0.0, dsum_dT = 0.0;
131 42 : for (std::size_t i = 0; i < _amu.size(); ++i)
132 : {
133 35 : sum += _amu[i] * MathUtils::pow(logTstar, i);
134 35 : dsum_dT += i * _amu[i] * MathUtils::pow(logTstar, i) / (temperature * logTstar);
135 : }
136 :
137 : const Real mu0 =
138 7 : 0.021357 * std::sqrt(1000.0 * _Mh2 * temperature) / (0.297 * 0.297 * std::exp(sum));
139 7 : const Real dmu0_dT = 21.357 * _Mh2 * (1.0 - 2.0 * temperature * dsum_dT) * std::exp(-sum) /
140 7 : (2.0 * std::sqrt(1000.0 * _Mh2 * temperature) * 0.297 * 0.297);
141 :
142 : // The excess contribution due to density
143 : Real sumr = 0.0;
144 56 : for (std::size_t i = 0; i < _bmu.size(); ++i)
145 49 : sumr += _bmu[i];
146 :
147 7 : const Real mu1 = MathUtils::pow(0.297, 3) * sumr * mu0 / Tstar;
148 : const Real dmu1_dT =
149 7 : MathUtils::pow(0.297, 3) * sumr * (dmu0_dT / Tstar - mu0 / (30.41 * Tstar * Tstar));
150 :
151 : // The viscosity and derivatives are then
152 7 : const Real exponent = _cmu[1] * Tr + _cmu[2] / Tr + _cmu[3] * rhor * rhor / (_cmu[4] + Tr) +
153 7 : _cmu[5] * MathUtils::pow(rhor, 6);
154 : const Real dexponent_drho =
155 7 : (2.0 * _cmu[3] * rhor / (_cmu[4] + Tr) + 6.0 * _cmu[5] * MathUtils::pow(rhor, 5)) *
156 7 : drhor_drho;
157 : const Real dexponent_dT =
158 7 : (_cmu[1] - _cmu[2] / Tr / Tr - _cmu[3] * rhor * rhor / (_cmu[4] + Tr) / (_cmu[4] + Tr)) *
159 7 : dTr_dT;
160 :
161 7 : mu = (mu0 + mu1 * density + _cmu[0] * rhor * rhor * std::exp(exponent)) * 1.0e-6;
162 7 : dmu_drho =
163 7 : (mu1 + _cmu[0] * rhor * std::exp(exponent) * (2.0 * drhor_drho + rhor * dexponent_drho)) *
164 : 1.0e-6;
165 7 : dmu_dT =
166 7 : (dmu0_dT + density * dmu1_dT + _cmu[0] * rhor * rhor * dexponent_dT * std::exp(exponent)) *
167 7 : 1.0e-6 +
168 7 : dmu_drho * ddensity_dT;
169 7 : }
170 :
171 : Real
172 12 : HydrogenFluidProperties::mu_from_p_T(Real pressure, Real temperature) const
173 : {
174 : // Require density first
175 12 : const Real density = rho_from_p_T(pressure, temperature);
176 12 : return mu_from_rho_T(density, temperature);
177 : }
178 :
179 : void
180 4 : HydrogenFluidProperties::mu_from_p_T(
181 : Real pressure, Real temperature, Real & mu, Real & dmu_dp, Real & dmu_dT) const
182 : {
183 : Real rho, drho_dp, drho_dT;
184 4 : rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT);
185 :
186 : Real dmu_drho;
187 4 : mu_from_rho_T(rho, temperature, drho_dT, mu, dmu_drho, dmu_dT);
188 4 : dmu_dp = dmu_drho * drho_dp;
189 4 : }
190 :
191 : void
192 1 : HydrogenFluidProperties::rho_mu_from_p_T(Real pressure,
193 : Real temperature,
194 : Real & rho,
195 : Real & mu) const
196 : {
197 1 : rho = rho_from_p_T(pressure, temperature);
198 1 : mu = mu_from_rho_T(rho, temperature);
199 1 : }
200 :
201 : void
202 1 : HydrogenFluidProperties::rho_mu_from_p_T(Real pressure,
203 : Real temperature,
204 : Real & rho,
205 : Real & drho_dp,
206 : Real & drho_dT,
207 : Real & mu,
208 : Real & dmu_dp,
209 : Real & dmu_dT) const
210 : {
211 1 : rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT);
212 : Real dmu_drho;
213 1 : mu_from_rho_T(rho, temperature, drho_dT, mu, dmu_drho, dmu_dT);
214 1 : dmu_dp = dmu_drho * drho_dp;
215 1 : }
216 :
217 : Real
218 12 : HydrogenFluidProperties::k_from_rho_T(Real density, Real temperature) const
219 : {
220 : // // Scaled variables
221 12 : const Real Tr = temperature / 33.145;
222 12 : const Real rhor = density / 31.262;
223 :
224 : // The ideal gas component
225 : Real sum1 = 0.0;
226 96 : for (std::size_t i = 0; i < _a1k.size(); ++i)
227 84 : sum1 += _a1k[i] * MathUtils::pow(Tr, i);
228 :
229 : Real sum2 = 0.0;
230 60 : for (std::size_t i = 0; i < _a2k.size(); ++i)
231 48 : sum2 += _a2k[i] * MathUtils::pow(Tr, i);
232 :
233 12 : const Real lambda0 = sum1 / sum2;
234 :
235 : // The excess contribution due to density
236 : Real lambdah = 0.0;
237 72 : for (std::size_t i = 0; i < _b1k.size(); ++i)
238 60 : lambdah += (_b1k[i] + _b2k[i] * Tr) * MathUtils::pow(rhor, i + 1);
239 :
240 : // The critical enhancement
241 12 : const Real lambdac = 6.24e-4 / (-2.58e-7 + std::abs(Tr - 1.0)) *
242 12 : std::exp(-MathUtils::pow(0.837 * (rhor - 1.0), 2));
243 :
244 : // The thermal conductivity
245 12 : return lambda0 + lambdah + lambdac;
246 : }
247 :
248 : Real
249 9 : HydrogenFluidProperties::k_from_p_T(Real pressure, Real temperature) const
250 : {
251 : // Require density first
252 9 : const Real density = rho_from_p_T(pressure, temperature);
253 9 : return k_from_rho_T(density, temperature);
254 : }
255 :
256 : void
257 1 : HydrogenFluidProperties::k_from_p_T(
258 : Real pressure, Real temperature, Real & k, Real & dk_dp, Real & dk_dT) const
259 : {
260 1 : k = this->k_from_p_T(pressure, temperature);
261 : // Calculate derivatives using finite differences
262 : const Real eps = 1.0e-6;
263 1 : const Real peps = pressure * eps;
264 1 : const Real Teps = temperature * eps;
265 :
266 1 : dk_dp = (this->k_from_p_T(pressure + peps, temperature) - k) / peps;
267 1 : dk_dT = (this->k_from_p_T(pressure, temperature + Teps) - k) / Teps;
268 1 : }
269 :
270 : std::vector<Real>
271 1 : HydrogenFluidProperties::henryCoefficients() const
272 : {
273 1 : return {-4.73284, 6.08954, 6.06066};
274 : }
275 :
276 : Real
277 2 : HydrogenFluidProperties::vaporPressure(Real temperature) const
278 : {
279 2 : if (temperature < _T_triple || temperature > _T_critical)
280 0 : throw MooseException("Temperature is out of range in " + name() + ": vaporPressure()");
281 :
282 2 : const Real Tr = temperature / _T_critical;
283 2 : const Real theta = 1.0 - Tr;
284 :
285 2 : const Real logpressure = (-4.89789 * theta + 0.988588 * std::pow(theta, 1.5) +
286 2 : 0.349689 * Utility::pow<2>(theta) + 0.499356 * std::pow(theta, 2.85)) /
287 2 : Tr;
288 :
289 2 : return _p_critical * std::exp(logpressure);
290 : }
291 :
292 : void
293 0 : HydrogenFluidProperties::vaporPressure(Real, Real &, Real &) const
294 : {
295 0 : mooseError("vaporPressure() is not implemented");
296 : }
297 :
298 : Real
299 3 : HydrogenFluidProperties::alpha(Real delta, Real tau) const
300 : {
301 : // Ideal gas component of the Helmholtz free energy
302 3 : Real alpha0 = std::log(delta) + 1.5 * std::log(tau) - 1.4579856475 + 1.888076782 * tau;
303 :
304 18 : for (std::size_t i = 0; i < _a.size(); ++i)
305 15 : alpha0 += _a[i] * std::log(1.0 - std::exp(_b[i] * tau));
306 :
307 : // Residual component of the Helmholtz free energy
308 : Real alphar = 0.0;
309 :
310 24 : for (std::size_t i = 0; i < _t1.size(); ++i)
311 21 : alphar += _N1[i] * MathUtils::pow(delta, _d1[i]) * std::pow(tau, _t1[i]);
312 :
313 9 : for (std::size_t i = 0; i < _t2.size(); ++i)
314 6 : alphar += _N2[i] * MathUtils::pow(delta, _d2[i]) * std::pow(tau, _t2[i]) * std::exp(-delta);
315 :
316 18 : for (std::size_t i = 0; i < _t3.size(); ++i)
317 15 : alphar += _N3[i] * MathUtils::pow(delta, _d3[i]) * std::pow(tau, _t3[i]) *
318 15 : std::exp(_phi3[i] * Utility::pow<2>(delta - _D3[i]) +
319 15 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i]));
320 :
321 : // The Helmholtz free energy is the sum of these two
322 3 : return alpha0 + alphar;
323 : }
324 :
325 : Real
326 801 : HydrogenFluidProperties::dalpha_ddelta(Real delta, Real tau) const
327 : {
328 : // Ideal gas component of the Helmholtz free energy
329 801 : Real dalpha0 = 1.0 / delta;
330 :
331 : // Residual component of the Helmholtz free energy
332 : Real dalphar = 0.0;
333 :
334 6408 : for (std::size_t i = 0; i < _t1.size(); ++i)
335 5607 : dalphar += _N1[i] * _d1[i] * MathUtils::pow(delta, _d1[i]) * std::pow(tau, _t1[i]);
336 :
337 2403 : for (std::size_t i = 0; i < _t2.size(); ++i)
338 1602 : dalphar += _N2[i] * MathUtils::pow(delta, _d2[i]) * std::pow(tau, _t2[i]) * std::exp(-delta) *
339 1602 : (_d2[i] - delta);
340 :
341 4806 : for (std::size_t i = 0; i < _t3.size(); ++i)
342 4005 : dalphar += _N3[i] * MathUtils::pow(delta, _d3[i]) * std::pow(tau, _t3[i]) *
343 4005 : std::exp(_phi3[i] * Utility::pow<2>(delta - _D3[i]) +
344 4005 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
345 4005 : (_d3[i] + delta * (2.0 * _phi3[i] * (delta - _D3[i])));
346 :
347 : // The Helmholtz free energy is the sum of these two
348 801 : return dalpha0 + dalphar / delta;
349 : }
350 :
351 : Real
352 23 : HydrogenFluidProperties::dalpha_dtau(Real delta, Real tau) const
353 : {
354 : // Ideal gas component of the Helmholtz free energy
355 23 : Real dalpha0 = 1.5 / tau + 1.888076782;
356 :
357 138 : for (std::size_t i = 0; i < _a.size(); ++i)
358 115 : dalpha0 += _a[i] * _b[i] * (1.0 - 1.0 / (1.0 - std::exp(_b[i] * tau)));
359 :
360 : // Residual component of the Helmholtz free energy
361 : Real dalphar = 0.0;
362 :
363 184 : for (std::size_t i = 0; i < _t1.size(); ++i)
364 161 : dalphar += _N1[i] * _t1[i] * MathUtils::pow(delta, _d1[i]) * std::pow(tau, _t1[i]);
365 :
366 69 : for (std::size_t i = 0; i < _t2.size(); ++i)
367 46 : dalphar +=
368 46 : _N2[i] * _t2[i] * MathUtils::pow(delta, _d2[i]) * std::pow(tau, _t2[i]) * std::exp(-delta);
369 :
370 138 : for (std::size_t i = 0; i < _t3.size(); ++i)
371 115 : dalphar += _N3[i] * MathUtils::pow(delta, _d3[i]) * std::pow(tau, _t3[i]) *
372 115 : std::exp(_phi3[i] * Utility::pow<2>(delta - _D3[i]) +
373 115 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
374 115 : (_t3[i] + tau * (2.0 * _beta3[i] * (tau - _gamma3[i])));
375 :
376 : // The Helmholtz free energy is the sum of these two
377 23 : return dalpha0 + dalphar / tau;
378 : }
379 :
380 : Real
381 20 : HydrogenFluidProperties::d2alpha_ddelta2(Real delta, Real tau) const
382 : {
383 : // Ideal gas component of the Helmholtz free energy
384 20 : Real dalpha0 = -1.0 / delta / delta;
385 :
386 : // Residual component of the Helmholtz free energy
387 : Real dalphar = 0.0;
388 :
389 160 : for (std::size_t i = 0; i < _t1.size(); ++i)
390 140 : dalphar +=
391 140 : _N1[i] * _d1[i] * (_d1[i] - 1.0) * MathUtils::pow(delta, _d1[i]) * std::pow(tau, _t1[i]);
392 :
393 60 : for (std::size_t i = 0; i < _t2.size(); ++i)
394 40 : dalphar += _N2[i] * MathUtils::pow(delta, _d2[i]) * std::pow(tau, _t2[i]) * std::exp(-delta) *
395 40 : (delta * delta - 2.0 * _d2[i] * delta + _d2[i] * (_d2[i] - 1.0));
396 :
397 120 : for (std::size_t i = 0; i < _t3.size(); ++i)
398 100 : dalphar += _N3[i] * MathUtils::pow(delta, _d3[i]) * std::pow(tau, _t3[i]) *
399 100 : std::exp(_phi3[i] * Utility::pow<2>(delta - _D3[i]) +
400 100 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
401 100 : (_d3[i] * _d3[i] +
402 100 : 2.0 * delta * delta * _phi3[i] *
403 100 : (1.0 + 2.0 * _phi3[i] * (_D3[i] - delta) * (_D3[i] - delta)) +
404 100 : _d3[i] * (4.0 * delta * _phi3[i] * (delta - _D3[i]) - 1.0));
405 :
406 : // The Helmholtz free energy is the sum of these two
407 20 : return dalpha0 + dalphar / delta / delta;
408 : }
409 :
410 : Real
411 13 : HydrogenFluidProperties::d2alpha_dtau2(Real delta, Real tau) const
412 : {
413 : // Ideal gas component of the Helmholtz free energy
414 13 : Real dalpha0 = -1.5 / tau / tau;
415 :
416 78 : for (std::size_t i = 0; i < _a.size(); ++i)
417 : {
418 65 : Real exptau = std::exp(_b[i] * tau);
419 65 : dalpha0 -= _a[i] * (_b[i] * _b[i] * exptau / (1.0 - exptau) * (exptau / (1.0 - exptau) + 1.0));
420 : }
421 :
422 : // Residual component of the Helmholtz free energy
423 : Real dalphar = 0.0;
424 :
425 104 : for (std::size_t i = 0; i < _t1.size(); ++i)
426 91 : dalphar +=
427 91 : _N1[i] * _t1[i] * (_t1[i] - 1.0) * MathUtils::pow(delta, _d1[i]) * std::pow(tau, _t1[i]);
428 :
429 39 : for (std::size_t i = 0; i < _t2.size(); ++i)
430 26 : dalphar += _N2[i] * _t2[i] * (_t2[i] - 1.0) * MathUtils::pow(delta, _d2[i]) *
431 26 : std::pow(tau, _t2[i]) * std::exp(-delta);
432 :
433 78 : for (std::size_t i = 0; i < _t3.size(); ++i)
434 65 : dalphar += _N3[i] * MathUtils::pow(delta, _d3[i]) * std::pow(tau, _t3[i]) *
435 65 : std::exp(_phi3[i] * Utility::pow<2>(delta - _D3[i]) +
436 65 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
437 65 : (_t3[i] * _t3[i] +
438 65 : 2.0 * _beta3[i] * tau * tau *
439 65 : (1.0 + 2.0 * _beta3[i] * MathUtils::pow(tau - _gamma3[i], 2)) -
440 65 : _t3[i] * (1.0 + 4.0 * _beta3[i] * tau * (tau - _gamma3[i])));
441 :
442 : // The Helmholtz free energy is the sum of these two
443 13 : return dalpha0 + dalphar / tau / tau;
444 : }
445 :
446 : Real
447 20 : HydrogenFluidProperties::d2alpha_ddeltatau(Real delta, Real tau) const
448 : {
449 : // Residual component of the Helmholtz free energy
450 : Real dalphar = 0.0;
451 :
452 160 : for (std::size_t i = 0; i < _t1.size(); ++i)
453 140 : dalphar += _N1[i] * _d1[i] * _t1[i] * std::pow(delta, _d1[i]) * std::pow(tau, _t1[i]);
454 :
455 60 : for (std::size_t i = 0; i < _t2.size(); ++i)
456 40 : dalphar += _N2[i] * _t2[i] * std::pow(delta, _d2[i]) * std::pow(tau, _t2[i]) *
457 40 : std::exp(-delta) * (_d2[i] - delta);
458 :
459 120 : for (std::size_t i = 0; i < _t3.size(); ++i)
460 100 : dalphar += _N3[i] * std::pow(delta, _d3[i]) * std::pow(tau, _t3[i]) *
461 100 : std::exp(_phi3[i] * Utility::pow<2>(delta - _D3[i]) +
462 100 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
463 100 : (_d3[i] + delta * (2.0 * _phi3[i] * (delta - _D3[i]))) *
464 100 : (_t3[i] + 2.0 * _beta3[i] * tau * (tau - _gamma3[i]));
465 :
466 : // The Helmholtz free energy is the sum of these two
467 20 : return dalphar / delta / tau;
468 : }
|
