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 "NitrogenFluidProperties.h"
11 : #include "Conversion.h"
12 : #include "MathUtils.h"
13 : #include "libmesh/utility.h"
14 :
15 : registerMooseObject("FluidPropertiesApp", NitrogenFluidProperties);
16 :
17 : InputParameters
18 22 : NitrogenFluidProperties::validParams()
19 : {
20 22 : InputParameters params = HelmholtzFluidProperties::validParams();
21 22 : params.addClassDescription("Fluid properties for Nitrogen (N2)");
22 22 : return params;
23 0 : }
24 :
25 11 : NitrogenFluidProperties::NitrogenFluidProperties(const InputParameters & parameters)
26 : : HelmholtzFluidProperties(parameters),
27 11 : _Mn2(28.01348e-3),
28 11 : _p_critical(3.3958e6),
29 11 : _T_critical(126.192),
30 11 : _rho_molar_critical(11.1839),
31 11 : _rho_critical(313.3),
32 11 : _p_triple(12.523e3),
33 11 : _T_triple(63.151)
34 : {
35 11 : }
36 :
37 : std::string
38 1 : NitrogenFluidProperties::fluidName() const
39 : {
40 1 : return "nitrogen";
41 : }
42 :
43 : Real
44 801 : NitrogenFluidProperties::molarMass() const
45 : {
46 801 : return _Mn2;
47 : }
48 :
49 : Real
50 1 : NitrogenFluidProperties::criticalPressure() const
51 : {
52 1 : return _p_critical;
53 : }
54 :
55 : Real
56 724 : NitrogenFluidProperties::criticalTemperature() const
57 : {
58 724 : return _T_critical;
59 : }
60 :
61 : Real
62 724 : NitrogenFluidProperties::criticalDensity() const
63 : {
64 724 : return _rho_critical;
65 : }
66 :
67 : Real
68 1 : NitrogenFluidProperties::triplePointPressure() const
69 : {
70 1 : return _p_triple;
71 : }
72 :
73 : Real
74 1 : NitrogenFluidProperties::triplePointTemperature() const
75 : {
76 1 : return _T_triple;
77 : }
78 :
79 : Real
80 31 : NitrogenFluidProperties::mu_from_rho_T(Real density, Real temperature) const
81 : {
82 : // Scale the input density and temperature
83 31 : const Real delta = density / _rho_critical;
84 31 : const Real tau = _T_critical / temperature;
85 31 : const Real logTstar = std::log(temperature / 98.94);
86 :
87 : // The dilute gas component
88 : Real logOmega = 0.0;
89 186 : for (std::size_t i = 0; i < _bmu.size(); ++i)
90 155 : logOmega += _bmu[i] * MathUtils::pow(logTstar, i);
91 :
92 : const Real mu0 =
93 31 : 0.0266958 * std::sqrt(1000.0 * _Mn2 * temperature) / (0.3656 * 0.3656 * std::exp(logOmega));
94 :
95 : // The residual component
96 : Real mur = 0.0;
97 186 : for (std::size_t i = 0; i < _Nmu.size(); ++i)
98 155 : mur += _Nmu[i] * std::pow(tau, _tmu[i]) * MathUtils::pow(delta, _dmu[i]) *
99 155 : std::exp(-_gammamu[i] * MathUtils::pow(delta, _lmu[i]));
100 :
101 : // The viscosity in Pa.s
102 31 : return (mu0 + mur) * 1.0e-6;
103 : }
104 :
105 : void
106 7 : NitrogenFluidProperties::mu_from_rho_T(Real density,
107 : Real temperature,
108 : Real ddensity_dT,
109 : Real & mu,
110 : Real & dmu_drho,
111 : Real & dmu_dT) const
112 : {
113 : // Scale the input density and temperature
114 7 : const Real delta = density / _rho_critical;
115 7 : const Real tau = _T_critical / temperature;
116 7 : const Real logTstar = std::log(temperature / 98.94);
117 :
118 : // The dilute gas component
119 : Real logOmega = 0.0, dlogOmega_dT = 0.0;
120 42 : for (std::size_t i = 0; i < _bmu.size(); ++i)
121 : {
122 35 : logOmega += _bmu[i] * MathUtils::pow(logTstar, i);
123 35 : dlogOmega_dT += i * _bmu[i] * MathUtils::pow(logTstar, i) / (temperature * logTstar);
124 : }
125 :
126 : const Real mu0 =
127 7 : 0.0266958 * std::sqrt(1000.0 * _Mn2 * temperature) / (0.3656 * 0.3656 * std::exp(logOmega));
128 7 : const Real dmu0_dT = 26.6958 * _Mn2 * (1.0 - 2.0 * temperature * dlogOmega_dT) *
129 7 : std::exp(-logOmega) /
130 7 : (2.0 * std::sqrt(1000.0 * _Mn2 * temperature) * 0.3656 * 0.3656);
131 :
132 : // The residual component
133 : Real mur = 0.0, dmur_drho = 0.0, dmur_dT = 0.0;
134 : Real term;
135 42 : for (std::size_t i = 0; i < _Nmu.size(); ++i)
136 : {
137 35 : term = _Nmu[i] * std::pow(tau, _tmu[i]) * MathUtils::pow(delta, _dmu[i]) *
138 35 : std::exp(-_gammamu[i] * MathUtils::pow(delta, _lmu[i]));
139 35 : mur += term;
140 35 : dmur_drho += (_dmu[i] - _lmu[i] * _gammamu[i] * MathUtils::pow(delta, _lmu[i])) * term / delta /
141 35 : _rho_molar_critical / (1000.0 * _Mn2);
142 35 : dmur_dT += -_tmu[i] * term / temperature;
143 : }
144 :
145 : // The viscosity in Pa.s
146 7 : mu = (mu0 + mur) * 1.0e-6;
147 7 : dmu_drho = dmur_drho * 1.0e-6;
148 7 : dmu_dT = (dmu0_dT + dmur_dT) * 1.0e-6 + dmu_drho * ddensity_dT;
149 7 : }
150 :
151 : Real
152 11 : NitrogenFluidProperties::mu_from_p_T(Real pressure, Real temperature) const
153 : {
154 : // Require density first
155 11 : const Real density = rho_from_p_T(pressure, temperature);
156 11 : return mu_from_rho_T(density, temperature);
157 : }
158 :
159 : void
160 4 : NitrogenFluidProperties::mu_from_p_T(
161 : Real pressure, Real temperature, Real & mu, Real & dmu_dp, Real & dmu_dT) const
162 : {
163 : Real rho, drho_dp, drho_dT;
164 4 : rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT);
165 :
166 : Real dmu_drho;
167 4 : mu_from_rho_T(rho, temperature, drho_dT, mu, dmu_drho, dmu_dT);
168 4 : dmu_dp = dmu_drho * drho_dp;
169 4 : }
170 :
171 : void
172 1 : NitrogenFluidProperties::rho_mu_from_p_T(Real pressure,
173 : Real temperature,
174 : Real & rho,
175 : Real & mu) const
176 : {
177 1 : rho = rho_from_p_T(pressure, temperature);
178 1 : mu = mu_from_rho_T(rho, temperature);
179 1 : }
180 :
181 : void
182 1 : NitrogenFluidProperties::rho_mu_from_p_T(Real pressure,
183 : Real temperature,
184 : Real & rho,
185 : Real & drho_dp,
186 : Real & drho_dT,
187 : Real & mu,
188 : Real & dmu_dp,
189 : Real & dmu_dT) const
190 : {
191 1 : rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT);
192 : Real dmu_drho;
193 1 : mu_from_rho_T(rho, temperature, drho_dT, mu, dmu_drho, dmu_dT);
194 1 : dmu_dp = dmu_drho * drho_dp;
195 1 : }
196 :
197 : Real
198 11 : NitrogenFluidProperties::k_from_rho_T(Real density, Real temperature) const
199 : {
200 : // Scale the input density and temperature
201 11 : const Real delta = density / _rho_critical;
202 11 : const Real tau = _T_critical / temperature;
203 :
204 : // The dilute gas component
205 : const Real lambda0 =
206 11 : 1.511 * mu_from_rho_T(0.0, temperature) * 1.0e6 + 2.117 / tau - 3.332 * std::pow(tau, -0.7);
207 :
208 : // The residual component
209 : Real lambdar = 0.0;
210 77 : for (std::size_t i = 0; i < _Nk.size(); ++i)
211 66 : lambdar += _Nk[i] * std::pow(tau, _tk[i]) * MathUtils::pow(delta, _dk[i]) *
212 66 : std::exp(-_gammak[i] * MathUtils::pow(delta, _lk[i]));
213 :
214 : // The thermal conductivity (note: critical enhancement not implemented)
215 11 : return (lambda0 + lambdar) * 1.0e-3;
216 : }
217 :
218 : Real
219 9 : NitrogenFluidProperties::k_from_p_T(Real pressure, Real temperature) const
220 : {
221 : // Require density first
222 9 : const Real density = rho_from_p_T(pressure, temperature);
223 9 : return k_from_rho_T(density, temperature);
224 : }
225 :
226 : void
227 1 : NitrogenFluidProperties::k_from_p_T(
228 : Real pressure, Real temperature, Real & k, Real & dk_dp, Real & dk_dT) const
229 : {
230 1 : k = this->k_from_p_T(pressure, temperature);
231 : // Calculate derivatives using finite differences
232 : const Real eps = 1.0e-6;
233 1 : const Real peps = pressure * eps;
234 1 : const Real Teps = temperature * eps;
235 :
236 1 : dk_dp = (this->k_from_p_T(pressure + peps, temperature) - k) / peps;
237 1 : dk_dT = (this->k_from_p_T(pressure, temperature + Teps) - k) / Teps;
238 1 : }
239 :
240 : std::vector<Real>
241 1 : NitrogenFluidProperties::henryCoefficients() const
242 : {
243 1 : return {-9.67578, 4.72162, 11.70585};
244 : }
245 :
246 : Real
247 3 : NitrogenFluidProperties::vaporPressure(Real temperature) const
248 : {
249 3 : if (temperature < _T_triple || temperature > _T_critical)
250 0 : throw MooseException("Temperature is out of range in " + name() + ": vaporPressure()");
251 :
252 3 : const Real Tr = temperature / _T_critical;
253 3 : const Real theta = 1.0 - Tr;
254 :
255 : const Real logpressure =
256 3 : (-6.12445284 * theta + 1.2632722 * std::pow(theta, 1.5) - 0.765910082 * std::pow(theta, 2.5) -
257 3 : 1.77570564 * Utility::pow<5>(theta)) /
258 3 : Tr;
259 :
260 3 : return _p_critical * std::exp(logpressure);
261 : }
262 :
263 : void
264 0 : NitrogenFluidProperties::vaporPressure(Real, Real &, Real &) const
265 : {
266 0 : mooseError("vaporPressure() is not implemented");
267 : }
268 :
269 : Real
270 3 : NitrogenFluidProperties::saturatedLiquidDensity(Real temperature) const
271 : {
272 3 : if (temperature < _T_triple || temperature > _T_critical)
273 0 : throw MooseException("Temperature is out of range in " + name() + ": vaporPressure()");
274 :
275 3 : const Real Tr = temperature / _T_critical;
276 3 : const Real theta = 1.0 - Tr;
277 :
278 : const Real logpressure =
279 3 : 1.48654237 * std::pow(theta, 0.3294) - 0.280476066 * std::pow(theta, 2.0 / 3.0) +
280 3 : 0.0894143085 * std::pow(theta, 8.0 / 3.0) - 0.119879866 * std::pow(theta, 35.0 / 6.0);
281 :
282 3 : return _rho_critical * std::exp(logpressure);
283 : }
284 :
285 : Real
286 3 : NitrogenFluidProperties::saturatedVaporDensity(Real temperature) const
287 : {
288 3 : if (temperature < _T_triple || temperature > _T_critical)
289 0 : throw MooseException("Temperature is out of range in " + name() + ": vaporPressure()");
290 :
291 3 : const Real Tr = temperature / _T_critical;
292 3 : const Real theta = 1.0 - Tr;
293 :
294 : const Real logpressure =
295 3 : (-1.70127164 * std::pow(theta, 0.34) - 3.70402649 * std::pow(theta, 5.0 / 6.0) +
296 3 : 1.29859383 * std::pow(theta, 7.0 / 6.0) - 0.561424977 * std::pow(theta, 13.0 / 6.0) -
297 3 : 2.68505381 * std::pow(theta, 14.0 / 3.0)) /
298 3 : Tr;
299 :
300 3 : return _rho_critical * std::exp(logpressure);
301 : }
302 :
303 : Real
304 3 : NitrogenFluidProperties::alpha(Real delta, Real tau) const
305 : {
306 : // Ideal gas component of the Helmholtz free energy
307 3 : const Real alpha0 = std::log(delta) + _a[0] * std::log(tau) + _a[1] + _a[2] * tau + _a[3] / tau +
308 3 : _a[4] / Utility::pow<2>(tau) + _a[5] / Utility::pow<3>(tau) +
309 3 : _a[6] * std::log(1.0 - std::exp(-_a[7] * tau));
310 :
311 : // Residual component of the Helmholtz free energy
312 : Real alphar = 0.0;
313 :
314 21 : for (std::size_t i = 0; i < _N1.size(); ++i)
315 18 : alphar += _N1[i] * MathUtils::pow(delta, _i1[i]) * std::pow(tau, _j1[i]);
316 :
317 81 : for (std::size_t i = 0; i < _N2.size(); ++i)
318 78 : alphar += _N2[i] * MathUtils::pow(delta, _i2[i]) * std::pow(tau, _j2[i]) *
319 78 : std::exp(-MathUtils::pow(delta, _l2[i]));
320 :
321 15 : for (std::size_t i = 0; i < _N3.size(); ++i)
322 12 : alphar += _N3[i] * MathUtils::pow(delta, _i3[i]) * std::pow(tau, _j3[i]) *
323 12 : std::exp(-_phi3[i] * Utility::pow<2>(delta - 1.0) -
324 12 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i]));
325 :
326 : // The Helmholtz free energy is the sum of these two
327 3 : return alpha0 + alphar;
328 : }
329 :
330 : Real
331 706 : NitrogenFluidProperties::dalpha_ddelta(Real delta, Real tau) const
332 : {
333 : // Ideal gas component of the Helmholtz free energy
334 706 : const Real dalpha0 = 1.0 / delta;
335 :
336 : // Residual component of the Helmholtz free energy
337 : Real dalphar = 0.0;
338 :
339 4942 : for (std::size_t i = 0; i < _N1.size(); ++i)
340 4236 : dalphar += _N1[i] * _i1[i] * MathUtils::pow(delta, _i1[i]) * std::pow(tau, _j1[i]);
341 :
342 19062 : for (std::size_t i = 0; i < _N2.size(); ++i)
343 18356 : dalphar += _N2[i] * MathUtils::pow(delta, _i2[i]) * std::pow(tau, _j2[i]) *
344 18356 : std::exp(-MathUtils::pow(delta, _l2[i])) *
345 18356 : (_i2[i] - _l2[i] * MathUtils::pow(delta, _l2[i]));
346 :
347 3530 : for (std::size_t i = 0; i < _N3.size(); ++i)
348 2824 : dalphar += _N3[i] * MathUtils::pow(delta, _i3[i]) * std::pow(tau, _j3[i]) *
349 2824 : std::exp(-_phi3[i] * Utility::pow<2>(delta - 1.0) -
350 2824 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
351 2824 : (_i3[i] - 2.0 * delta * _phi3[i] * (delta - 1.0));
352 :
353 : // The Helmholtz free energy is the sum of these two
354 706 : return dalpha0 + dalphar / delta;
355 : }
356 :
357 : Real
358 23 : NitrogenFluidProperties::dalpha_dtau(Real delta, Real tau) const
359 : {
360 : // Ideal gas component of the Helmholtz free energy
361 23 : const Real dalpha0 = _a[0] + _a[2] * tau - _a[3] / tau - 2.0 * _a[4] / Utility::pow<2>(tau) -
362 23 : 3.0 * _a[5] / Utility::pow<3>(tau) +
363 23 : _a[6] * _a[7] * tau / (std::exp(_a[7] * tau) - 1.0);
364 :
365 : // Residual component of the Helmholtz free energy
366 : Real dalphar = 0.0;
367 :
368 161 : for (std::size_t i = 0; i < _N1.size(); ++i)
369 138 : dalphar += _N1[i] * _j1[i] * MathUtils::pow(delta, _i1[i]) * std::pow(tau, _j1[i]);
370 :
371 621 : for (std::size_t i = 0; i < _N2.size(); ++i)
372 598 : dalphar += _N2[i] * _j2[i] * MathUtils::pow(delta, _i2[i]) * std::pow(tau, _j2[i]) *
373 598 : std::exp(-MathUtils::pow(delta, _l2[i]));
374 :
375 115 : for (std::size_t i = 0; i < _N3.size(); ++i)
376 92 : dalphar += _N3[i] * MathUtils::pow(delta, _i3[i]) * std::pow(tau, _j3[i]) *
377 92 : std::exp(-_phi3[i] * Utility::pow<2>(delta - 1.0) -
378 92 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
379 92 : (_j3[i] - 2.0 * tau * _beta3[i] * (tau - _gamma3[i]));
380 :
381 : // The Helmholtz free energy is the sum of these two
382 23 : return (dalpha0 + dalphar) / tau;
383 : }
384 :
385 : Real
386 20 : NitrogenFluidProperties::d2alpha_ddelta2(Real delta, Real tau) const
387 : {
388 : // Ideal gas component of the Helmholtz free energy
389 20 : const Real dalpha0 = -1.0 / delta / delta;
390 :
391 : // Residual component of the Helmholtz free energy
392 : Real dalphar = 0.0;
393 :
394 140 : for (std::size_t i = 0; i < _N1.size(); ++i)
395 120 : dalphar +=
396 120 : _N1[i] * _i1[i] * (_i1[i] - 1.0) * MathUtils::pow(delta, _i1[i]) * std::pow(tau, _j1[i]);
397 :
398 540 : for (std::size_t i = 0; i < _N2.size(); ++i)
399 520 : dalphar += _N2[i] * MathUtils::pow(delta, _i2[i]) * std::pow(tau, _j2[i]) *
400 520 : std::exp(-MathUtils::pow(delta, _l2[i])) *
401 520 : ((_i2[i] - _l2[i] * MathUtils::pow(delta, _l2[i])) *
402 520 : (_i2[i] - 1.0 - _l2[i] * MathUtils::pow(delta, _l2[i])) -
403 520 : _l2[i] * _l2[i] * MathUtils::pow(delta, _l2[i]));
404 :
405 100 : for (std::size_t i = 0; i < _N3.size(); ++i)
406 80 : dalphar += _N3[i] * MathUtils::pow(delta, _i3[i]) * std::pow(tau, _j3[i]) *
407 80 : std::exp(-_phi3[i] * Utility::pow<2>(delta - 1.0) -
408 80 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
409 80 : (Utility::pow<2>(_i3[i] - 2.0 * delta * _phi3[i] * (delta - 1.0)) - _i3[i] -
410 80 : 2.0 * delta * delta * _phi3[i]);
411 :
412 : // The Helmholtz free energy
413 20 : return dalpha0 + dalphar / delta / delta;
414 : }
415 :
416 : Real
417 13 : NitrogenFluidProperties::d2alpha_dtau2(Real delta, Real tau) const
418 : {
419 : // Ideal gas component of the Helmholtz free energy
420 13 : const Real dalpha0 = -_a[0] + 2.0 * _a[3] / tau + 6.0 * _a[4] / Utility::pow<2>(tau) +
421 13 : 12.0 * _a[5] / Utility::pow<3>(tau) -
422 13 : _a[6] * _a[7] * _a[7] * tau * tau * std::exp(_a[7] * tau) /
423 13 : Utility::pow<2>(std::exp(_a[7] * tau) - 1.0);
424 :
425 : // Residual component of the Helmholtz free energy
426 : Real dalphar = 0.0;
427 :
428 91 : for (std::size_t i = 0; i < _N1.size(); ++i)
429 78 : dalphar +=
430 78 : _N1[i] * _j1[i] * (_j1[i] - 1.0) * MathUtils::pow(delta, _i1[i]) * std::pow(tau, _j1[i]);
431 :
432 351 : for (std::size_t i = 0; i < _N2.size(); ++i)
433 338 : dalphar += _N2[i] * _j2[i] * (_j2[i] - 1.0) * MathUtils::pow(delta, _i2[i]) *
434 338 : std::pow(tau, _j2[i]) * std::exp(-MathUtils::pow(delta, _l2[i]));
435 :
436 65 : for (std::size_t i = 0; i < _N3.size(); ++i)
437 52 : dalphar += _N3[i] * MathUtils::pow(delta, _i3[i]) * std::pow(tau, _j3[i]) *
438 52 : std::exp(-_phi3[i] * Utility::pow<2>(delta - 1.0) -
439 52 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
440 52 : (Utility::pow<2>(_j3[i] - 2.0 * tau * _beta3[i] * (tau - _gamma3[i])) - _j3[i] -
441 52 : 2.0 * tau * tau * _beta3[i]);
442 :
443 : // The Helmholtz free energy is the sum of these two
444 13 : return (dalpha0 + dalphar) / tau / tau;
445 : }
446 :
447 : Real
448 20 : NitrogenFluidProperties::d2alpha_ddeltatau(Real delta, Real tau) const
449 : {
450 : // Residual component of the Helmholtz free energy (second derivative of ideal
451 : // component wrt delta and tau is 0)
452 : Real dalphar = 0.0;
453 :
454 140 : for (std::size_t i = 0; i < _N1.size(); ++i)
455 120 : dalphar += _N1[i] * _i1[i] * _j1[i] * MathUtils::pow(delta, _i1[i]) * std::pow(tau, _j1[i]);
456 :
457 540 : for (std::size_t i = 0; i < _N2.size(); ++i)
458 520 : dalphar += _N2[i] * _j2[i] * MathUtils::pow(delta, _i2[i]) * std::pow(tau, _j2[i]) *
459 520 : std::exp(-MathUtils::pow(delta, _l2[i])) *
460 520 : (_i2[i] - _l2[i] * MathUtils::pow(delta, _l2[i]));
461 :
462 100 : for (std::size_t i = 0; i < _N3.size(); ++i)
463 80 : dalphar += _N3[i] * MathUtils::pow(delta, _i3[i]) * std::pow(tau, _j3[i]) *
464 80 : std::exp(-_phi3[i] * Utility::pow<2>(delta - 1.0) -
465 80 : _beta3[i] * Utility::pow<2>(tau - _gamma3[i])) *
466 80 : (_i3[i] - 2.0 * delta * _phi3[i] * (delta - 1.0)) *
467 80 : (_j3[i] - 2.0 * tau * _beta3[i] * (tau - _gamma3[i]));
468 :
469 : // The Helmholtz free energy
470 20 : return dalphar / delta / tau;
471 : }
|
