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 "CaloricallyImperfectGas.h"
11 : #include "Conversion.h"
12 : #include "metaphysicl/raw_type.h"
13 : #include "Function.h"
14 : #include "BrentsMethod.h"
15 :
16 : registerMooseObject("FluidPropertiesApp", CaloricallyImperfectGas);
17 :
18 : InputParameters
19 49 : CaloricallyImperfectGas::validParams()
20 : {
21 49 : InputParameters params = SinglePhaseFluidProperties::validParams();
22 49 : params += NaNInterface::validParams();
23 :
24 98 : params.addRequiredParam<Real>("molar_mass", "Constant molar mass of the fluid [kg/mol]");
25 98 : params.addRequiredParam<FunctionName>("mu", "Dynamic viscosity, [Pa-s]");
26 98 : params.addRequiredParam<FunctionName>("k", "Thermal conductivity, [W/(m-K)]");
27 98 : params.addParam<Real>("T_c", 0, "Critical temperature, [K]");
28 98 : params.addParam<Real>("rho_c", 0, "Critical density, [kg/m3]");
29 98 : params.addParam<Real>("e_c", 0, "Specific internal energy at the critical point, [J/kg]");
30 98 : params.addRequiredParam<FunctionName>(
31 : "e", "Specific internal energy as a function of temperature [J/kg]");
32 98 : params.addRequiredParam<Real>("min_temperature", "Smallest temperature for lookup tables [K]");
33 98 : params.addRequiredParam<Real>("max_temperature", "Largest temperature for lookup tables [K]");
34 98 : params.addParam<Real>(
35 98 : "temperature_resolution", 1, "Step size in temperature for creating the inverse lookup T(e)");
36 98 : params.addParam<bool>("out_of_bound_error", true, "Error if out of bounds");
37 49 : params.addClassDescription("Fluid properties for an ideal gas with imperfect caloric behavior.");
38 :
39 49 : return params;
40 0 : }
41 :
42 26 : CaloricallyImperfectGas::CaloricallyImperfectGas(const InputParameters & parameters)
43 : : SinglePhaseFluidProperties(parameters),
44 : NaNInterface(this),
45 26 : _molar_mass(getParam<Real>("molar_mass")),
46 26 : _R_specific(_R / _molar_mass),
47 52 : _T_c(getParam<Real>("T_c")),
48 52 : _rho_c(getParam<Real>("rho_c")),
49 52 : _e_c(getParam<Real>("e_c")),
50 52 : _min_temperature(getParam<Real>("min_temperature")),
51 52 : _max_temperature(getParam<Real>("max_temperature")),
52 52 : _delta_T(getParam<Real>("temperature_resolution")),
53 52 : _out_of_bound_error(getParam<bool>("out_of_bound_error")),
54 52 : _tol(1e-4)
55 : {
56 26 : }
57 :
58 : void
59 26 : CaloricallyImperfectGas::initialSetup()
60 : {
61 26 : SinglePhaseFluidProperties::initialSetup();
62 26 : _e_T = &getFunction("e");
63 26 : _mu_T = &getFunction("mu");
64 26 : _k_T = &getFunction("k");
65 :
66 26 : setupLookupTables();
67 24 : }
68 :
69 : void
70 26 : CaloricallyImperfectGas::setupLookupTables()
71 : {
72 : // estimate number of points in inverse lookup & then adjust
73 : // the _delta_T so (n - 1) * _delta_T = _max_temperature - _min_temperature
74 26 : unsigned int n = std::floor((_max_temperature - _min_temperature) / _delta_T) + 1;
75 26 : _delta_T = (_max_temperature - _min_temperature) / ((Real)n - 1.0);
76 :
77 : // ensure e(T) is monotonic
78 148451 : for (unsigned int j = 0; j < n; ++j)
79 : {
80 148427 : Real temperature = _min_temperature + j * _delta_T;
81 :
82 : // Note that the function behavior at the end points
83 : // can lead to de/dT = 0 & we want to allow that
84 148427 : Real cv = _e_T->timeDerivative(temperature, Point());
85 148427 : if (cv < 0 || (cv == 0 && j > 0 && j < n - 1))
86 2 : mooseError("e(T) is not monotonically increasing with T");
87 : }
88 :
89 : // backward lookup tables
90 24 : _min_e = _e_T->value(_min_temperature, Point());
91 24 : _max_e = _e_T->value(_max_temperature, Point());
92 24 : _delta_e = (_max_e - _min_e) / ((Real)n - 1.0);
93 24 : _min_h = _min_e + _R_specific * _min_temperature;
94 24 : _max_h = _max_e + _R_specific * _max_temperature;
95 24 : _delta_h = (_max_h - _min_h) / ((Real)n - 1.0);
96 :
97 24 : _T_e_lookup.resize(n);
98 24 : _T_h_lookup.resize(n);
99 24 : _Z_T_lookup.resize(n);
100 148248 : for (unsigned int j = 0; j < n; ++j)
101 : {
102 : // internal energy
103 : {
104 148224 : Real internal_energy = _min_e + j * _delta_e;
105 :
106 : // guarding against roundoff when summing
107 : // at least once I saw roundoff cause an error here
108 148224 : internal_energy = std::min(internal_energy, _max_e);
109 :
110 : // The temperature is found by e - e(T) = 0
111 148224 : Real min = _min_temperature;
112 148224 : Real max = _max_temperature;
113 1227694 : auto e_diff = [&internal_energy, this](Real x)
114 1227694 : { return this->e_from_p_T(0.0, x) - internal_energy; };
115 148224 : BrentsMethod::bracket(e_diff, min, max);
116 148224 : Real temperature = BrentsMethod::root(e_diff, min, max);
117 148224 : _T_e_lookup[j] = temperature;
118 : }
119 :
120 : // enthalpy
121 : {
122 148224 : Real enthalpy = _min_h + j * _delta_h;
123 :
124 : // guarding against roundoff when summing
125 : // at least once I saw roundoff cause an error here
126 148224 : enthalpy = std::min(enthalpy, _max_h);
127 :
128 : // The temperature is found by h - h(T) = 0
129 148224 : Real min = _min_temperature;
130 148224 : Real max = _max_temperature;
131 1222935 : auto h_diff = [&enthalpy, this](Real x) { return this->h_from_p_T(0.0, x) - enthalpy; };
132 148224 : BrentsMethod::bracket(h_diff, min, max);
133 296448 : Real temperature = BrentsMethod::root(h_diff, min, max);
134 148224 : _T_h_lookup[j] = temperature;
135 : }
136 :
137 : // Z(T)
138 : {
139 148224 : if (j == 0)
140 24 : _Z_T_lookup[j] = 0;
141 : else
142 : {
143 148200 : Real temperature = _min_temperature + j * _delta_T;
144 148200 : Real temperature_prev = _min_temperature + (j - 1) * _delta_T;
145 148200 : Real f1 = cv_from_T(temperature) / temperature;
146 148200 : Real f2 = cv_from_T(temperature_prev) / temperature_prev;
147 148200 : _Z_T_lookup[j] = _Z_T_lookup[j - 1] + 0.5 * _delta_T * (f1 + f2);
148 : }
149 : }
150 : }
151 24 : }
152 :
153 : std::string
154 0 : CaloricallyImperfectGas::fluidName() const
155 : {
156 0 : return "caloric_imperfect_gas";
157 : }
158 :
159 : Real
160 1229328 : CaloricallyImperfectGas::e_from_T(Real T) const
161 : {
162 1229328 : if (T < _min_temperature || T > _max_temperature)
163 0 : outOfBounds("e_from_T", T, _min_temperature, _max_temperature);
164 :
165 1229328 : return _e_T->value(T, Point());
166 : }
167 :
168 : Real
169 1223784 : CaloricallyImperfectGas::h_from_T(Real T) const
170 : {
171 1223784 : if (T < _min_temperature || T > _max_temperature)
172 0 : outOfBounds("h_from_T", T, _min_temperature, _max_temperature);
173 1223784 : return _e_T->value(T, Point()) + _R_specific * T;
174 : }
175 :
176 : Real
177 4924 : CaloricallyImperfectGas::cp_from_T(Real T) const
178 : {
179 4924 : if (T < _min_temperature || T > _max_temperature)
180 2 : outOfBounds("cp_from_T", T, _min_temperature, _max_temperature);
181 :
182 4922 : return _e_T->timeDerivative(T, Point()) + _R_specific;
183 : }
184 :
185 : void
186 1597 : CaloricallyImperfectGas::cv_from_T(Real T, Real & cv, Real & dcv_dT) const
187 : {
188 1597 : if (T < _min_temperature || T > _max_temperature)
189 0 : outOfBounds("cv_from_T (3 args)", T, _min_temperature, _max_temperature);
190 : Real pert = 1.0e-7;
191 1597 : cv = cv_from_T(T);
192 1597 : Real cv_pert = cv_from_T(T * (1 + pert));
193 1597 : dcv_dT = (cv_pert - cv) / (T * pert);
194 1597 : }
195 :
196 : void
197 1399 : CaloricallyImperfectGas::cp_from_T(Real T, Real & cp, Real & dcp_dT) const
198 : {
199 1399 : if (T < _min_temperature || T > _max_temperature)
200 0 : outOfBounds("cp_from_T (3 args)", T, _min_temperature, _max_temperature);
201 : Real pert = 1.0e-7;
202 1399 : cp = cp_from_T(T);
203 1399 : Real cp_pert = cp_from_T(T * (1 + pert));
204 1399 : dcp_dT = (cp_pert - cp) / (T * pert);
205 1399 : }
206 :
207 : Real
208 11170 : CaloricallyImperfectGas::T_from_e(Real e) const
209 : {
210 11170 : if (e < _min_e || e > _max_e)
211 0 : outOfBounds("T_from_e", e, _min_e, _max_e);
212 :
213 11170 : if (e < _min_e)
214 0 : return _T_e_lookup[0];
215 11170 : else if (e > _max_e)
216 0 : return _T_e_lookup.back();
217 :
218 11170 : unsigned int index = std::floor((e - _min_e) / _delta_e);
219 11170 : Real x = (e - _min_e - index * _delta_e) / _delta_e;
220 11170 : return x * _T_e_lookup[index + 1] + (1 - x) * _T_e_lookup[index];
221 : }
222 :
223 : Real
224 436 : CaloricallyImperfectGas::T_from_h(Real h) const
225 : {
226 436 : if (h < _min_h || h > _max_h)
227 0 : outOfBounds("h_from_e", h, _min_h, _max_h);
228 :
229 436 : if (h < _min_h)
230 0 : return _T_h_lookup[0];
231 436 : else if (h > _max_h)
232 0 : return _T_h_lookup.back();
233 :
234 436 : unsigned int index = std::floor((h - _min_h) / _delta_h);
235 436 : Real x = (h - _min_h - index * _delta_h) / _delta_h;
236 436 : return x * _T_h_lookup[index + 1] + (1 - x) * _T_h_lookup[index];
237 : }
238 :
239 : Real
240 3152 : CaloricallyImperfectGas::Z_from_T(Real T) const
241 : {
242 3152 : if (T < _min_temperature || T > _max_temperature)
243 0 : outOfBounds("Z_from_T", T, _min_temperature, _max_temperature);
244 :
245 3152 : unsigned int index = std::floor((T - _min_temperature) / _delta_T);
246 3152 : Real x = (T - _min_temperature - index * _delta_T) / _delta_T;
247 3152 : return x * _Z_T_lookup[index + 1] + (1 - x) * _Z_T_lookup[index];
248 : }
249 :
250 : Real
251 2548 : CaloricallyImperfectGas::p_from_v_e(Real v, Real e) const
252 : {
253 2548 : if (v == 0.0)
254 0 : return getNaN("Invalid value of specific volume detected (v = " + Moose::stringify(v) + ").");
255 :
256 2548 : return _R_specific * T_from_e(e) / v;
257 : }
258 :
259 : ADReal
260 198 : CaloricallyImperfectGas::p_from_v_e(const ADReal & v, const ADReal & e) const
261 : {
262 198 : if (v.value() == 0.0)
263 0 : return getNaN("Invalid value of specific volume detected (v = " + Moose::stringify(v.value()) +
264 0 : ").");
265 :
266 198 : Real x = 0;
267 : Real raw1 = v.value();
268 198 : Real raw2 = e.value();
269 198 : Real dxd1 = 0;
270 198 : Real dxd2 = 0;
271 198 : p_from_v_e(raw1, raw2, x, dxd1, dxd2);
272 :
273 198 : ADReal result = x;
274 198 : result.derivatives() = v.derivatives() * dxd1 + e.derivatives() * dxd2;
275 198 : return result;
276 : }
277 :
278 : void
279 201 : CaloricallyImperfectGas::p_from_v_e(Real v, Real e, Real & p, Real & dp_dv, Real & dp_de) const
280 : {
281 201 : p = p_from_v_e(v, e);
282 201 : Real T = T_from_e(e);
283 201 : dp_dv = -_R_specific * T / v / v;
284 201 : dp_de = _R_specific / v / cv_from_T(T);
285 201 : }
286 :
287 : Real
288 8019 : CaloricallyImperfectGas::T_from_v_e(Real /*v*/, Real e) const
289 : {
290 8019 : return T_from_e(e);
291 : }
292 :
293 : ADReal
294 396 : CaloricallyImperfectGas::T_from_v_e(const ADReal & v, const ADReal & e) const
295 : {
296 396 : if (v.value() == 0.0)
297 0 : return getNaN("Invalid value of specific volume detected (v = " + Moose::stringify(v.value()) +
298 0 : ").");
299 :
300 396 : Real x = 0;
301 : Real raw1 = v.value();
302 396 : Real raw2 = e.value();
303 396 : Real dxd1 = 0;
304 396 : Real dxd2 = 0;
305 396 : T_from_v_e(raw1, raw2, x, dxd1, dxd2);
306 :
307 396 : ADReal result = x;
308 396 : result.derivatives() = v.derivatives() * dxd1 + e.derivatives() * dxd2;
309 396 : return result;
310 : }
311 :
312 : void
313 402 : CaloricallyImperfectGas::T_from_v_e(Real v, Real e, Real & T, Real & dT_dv, Real & dT_de) const
314 : {
315 402 : T = T_from_e(e);
316 402 : dT_dv = 0.0;
317 402 : dT_de = 1.0 / cv_from_v_e(v, e);
318 402 : }
319 :
320 : Real
321 440 : CaloricallyImperfectGas::c_from_v_e(Real v, Real e) const
322 : {
323 440 : Real T = T_from_v_e(v, e);
324 :
325 440 : const Real c2 = gamma_from_v_e(v, e) * _R_specific * T;
326 440 : if (c2 < 0)
327 0 : return getNaN("Sound speed squared (gamma * R * T) is negative: c2 = " + Moose::stringify(c2) +
328 : ".");
329 :
330 440 : return std::sqrt(c2);
331 : }
332 :
333 : ADReal
334 198 : CaloricallyImperfectGas::c_from_v_e(const ADReal & v, const ADReal & e) const
335 : {
336 198 : const auto T = T_from_v_e(v, e);
337 :
338 396 : const auto c2 = gamma_from_v_e(v, e) * _R_specific * T;
339 198 : if (MetaPhysicL::raw_value(c2) < 0)
340 0 : return getNaN("Sound speed squared (gamma * R * T) is negative: c2 = " +
341 0 : Moose::stringify(MetaPhysicL::raw_value(c2)) + ".");
342 :
343 198 : return std::sqrt(c2);
344 : }
345 :
346 : void
347 3 : CaloricallyImperfectGas::c_from_v_e(Real v, Real e, Real & c, Real & dc_dv, Real & dc_de) const
348 : {
349 : Real T, dT_dv, dT_de;
350 3 : T_from_v_e(v, e, T, dT_dv, dT_de);
351 :
352 : Real gamma, dgamma_dv, dgamma_dT;
353 3 : gamma_from_v_e(v, e, gamma, dgamma_dv, dgamma_dT);
354 3 : c = std::sqrt(gamma * _R_specific * T);
355 :
356 3 : const Real dc_dT = 0.5 / c * _R_specific * (gamma + T * dgamma_dT);
357 3 : dc_dv = dc_dT * dT_dv;
358 3 : dc_de = dc_dT * dT_de;
359 3 : }
360 :
361 : Real
362 207 : CaloricallyImperfectGas::c_from_p_T(Real p, Real T) const
363 : {
364 207 : return std::sqrt(gamma_from_p_T(p, T) * _R_specific * T);
365 : }
366 :
367 : ADReal
368 198 : CaloricallyImperfectGas::c_from_p_T(const ADReal & p, const ADReal & T) const
369 : {
370 396 : return std::sqrt(gamma_from_p_T(p, T) * _R_specific * T);
371 : }
372 :
373 : void
374 2 : CaloricallyImperfectGas::c_from_p_T(
375 : const Real p, const Real T, Real & c, Real & dc_dp, Real & dc_dT) const
376 : {
377 : Real gamma, dgamma_dp, dgamma_dT;
378 2 : gamma_from_p_T(p, T, gamma, dgamma_dp, dgamma_dT);
379 2 : c = std::sqrt(gamma * _R_specific * T);
380 2 : dc_dp = 0;
381 2 : dc_dT = 0.5 / c * _R_specific * (gamma + T * dgamma_dT);
382 2 : }
383 :
384 : Real
385 1686 : CaloricallyImperfectGas::cp_from_v_e(Real v, Real e) const
386 : {
387 1686 : Real T = T_from_v_e(v, e);
388 1686 : return cp_from_T(T);
389 : }
390 :
391 : void
392 403 : CaloricallyImperfectGas::cp_from_v_e(Real v, Real e, Real & cp, Real & dcp_dv, Real & dcp_de) const
393 : {
394 403 : Real T = T_from_v_e(v, e);
395 : Real dcp_dT;
396 403 : cp_from_T(T, cp, dcp_dT);
397 403 : dcp_dv = 0.0;
398 403 : dcp_de = dcp_dT / cv_from_T(T);
399 403 : }
400 :
401 : Real
402 2088 : CaloricallyImperfectGas::cv_from_v_e(Real v, Real e) const
403 : {
404 2088 : Real T = T_from_v_e(v, e);
405 2088 : return cv_from_T(T);
406 : }
407 :
408 : ADReal
409 198 : CaloricallyImperfectGas::cv_from_v_e(ADReal v, ADReal e) const
410 : {
411 198 : Real x = 0;
412 198 : Real raw1 = v.value();
413 198 : Real raw2 = e.value();
414 198 : Real dxd1 = 0;
415 198 : Real dxd2 = 0;
416 198 : cv_from_v_e(raw1, raw2, x, dxd1, dxd2);
417 :
418 198 : ADReal result = x;
419 198 : result.derivatives() = v.derivatives() * dxd1 + e.derivatives() * dxd2;
420 198 : return result;
421 : }
422 :
423 : void
424 601 : CaloricallyImperfectGas::cv_from_v_e(Real v, Real e, Real & cv, Real & dcv_dv, Real & dcv_de) const
425 : {
426 601 : Real T = T_from_v_e(v, e);
427 : Real dcv_dT;
428 601 : cv_from_T(T, cv, dcv_dT);
429 601 : dcv_dv = 0.0;
430 601 : dcv_de = dcv_dT / cv;
431 601 : }
432 :
433 : Real
434 1042 : CaloricallyImperfectGas::gamma_from_v_e(Real v, Real e) const
435 : {
436 1042 : return cp_from_v_e(v, e) / cv_from_v_e(v, e);
437 : }
438 :
439 : ADReal
440 396 : CaloricallyImperfectGas::gamma_from_v_e(ADReal v, ADReal e) const
441 : {
442 396 : Real x = 0;
443 396 : Real raw1 = v.value();
444 396 : Real raw2 = e.value();
445 396 : Real dxd1 = 0;
446 396 : Real dxd2 = 0;
447 396 : gamma_from_v_e(raw1, raw2, x, dxd1, dxd2);
448 :
449 396 : ADReal result = x;
450 396 : result.derivatives() = v.derivatives() * dxd1 + e.derivatives() * dxd2;
451 396 : return result;
452 : }
453 :
454 : void
455 401 : CaloricallyImperfectGas::gamma_from_v_e(
456 : Real v, Real e, Real & gamma, Real & dgamma_dv, Real & dgamma_de) const
457 : {
458 : Real cp, dcp_dv, dcp_de;
459 401 : cp_from_v_e(v, e, cp, dcp_dv, dcp_de);
460 : Real cv, dcv_dv, dcv_de;
461 401 : cv_from_v_e(v, e, cv, dcv_dv, dcv_de);
462 401 : gamma = cp / cv;
463 401 : dgamma_dv = 0.0;
464 401 : dgamma_de = 1.0 / cv * dcp_de - gamma / cv * dcv_de;
465 401 : }
466 :
467 : void
468 994 : CaloricallyImperfectGas::gamma_from_p_T(
469 : Real p, Real T, Real & gamma, Real & dgamma_dp, Real & dgamma_dT) const
470 : {
471 : Real cp, dcp_dp, dcp_dT;
472 994 : cp_from_p_T(p, T, cp, dcp_dp, dcp_dT);
473 : Real cv, dcv_dp, dcv_dT;
474 994 : cv_from_p_T(p, T, cv, dcv_dp, dcv_dT);
475 994 : gamma = cp / cv;
476 994 : dgamma_dp = 0.0;
477 994 : dgamma_dT = 1.0 / cv * dcp_dT - gamma / cv * dcv_dT;
478 994 : }
479 :
480 : Real
481 216 : CaloricallyImperfectGas::gamma_from_p_T(Real p, Real T) const
482 : {
483 216 : return cp_from_p_T(p, T) / cv_from_p_T(p, T);
484 : }
485 :
486 : ADReal
487 990 : CaloricallyImperfectGas::gamma_from_p_T(ADReal p, ADReal T) const
488 : {
489 990 : Real x = 0;
490 990 : Real raw1 = p.value();
491 990 : Real raw2 = T.value();
492 990 : Real dxd1 = 0;
493 990 : Real dxd2 = 0;
494 990 : gamma_from_p_T(raw1, raw2, x, dxd1, dxd2);
495 990 : ADReal result = x;
496 990 : result.derivatives() = p.derivatives() * dxd1 + T.derivatives() * dxd2;
497 990 : return result;
498 : }
499 :
500 : Real
501 637 : CaloricallyImperfectGas::mu_from_v_e(Real v, Real e) const
502 : {
503 637 : const Real T = T_from_v_e(v, e);
504 637 : const Real p = p_from_v_e(v, e);
505 637 : return mu_from_p_T(p, T);
506 : }
507 :
508 : void
509 3 : CaloricallyImperfectGas::mu_from_v_e(Real v, Real e, Real & mu, Real & dmu_dv, Real & dmu_de) const
510 : {
511 3 : const Real T = T_from_v_e(v, e);
512 3 : const Real p = p_from_v_e(v, e);
513 : Real dmu_dp, dmu_dT;
514 3 : mu_from_p_T(p, T, mu, dmu_dp, dmu_dT);
515 : // only dmu_de = dmu/dT * dT/de = 1/cv * dmu/dT is nonzero
516 : // (dmu/dv)_e = 0 because e only depends on T and e is constant
517 3 : dmu_dv = 0.0;
518 3 : dmu_de = dmu_dT / cv_from_p_T(p, T);
519 3 : }
520 :
521 : Real
522 637 : CaloricallyImperfectGas::k_from_v_e(Real v, Real e) const
523 : {
524 637 : const Real T = T_from_v_e(v, e);
525 637 : const Real p = p_from_v_e(v, e);
526 637 : return k_from_p_T(p, T);
527 : }
528 :
529 : void
530 3 : CaloricallyImperfectGas::k_from_v_e(Real v, Real e, Real & k, Real & dk_dv, Real & dk_de) const
531 : {
532 3 : const Real T = T_from_v_e(v, e);
533 3 : const Real p = p_from_v_e(v, e);
534 : Real dk_dp, dk_dT;
535 3 : k_from_p_T(p, T, k, dk_dp, dk_dT);
536 : // only dk_de = dk/dT * dT/de = 1/cv * dk/dT is nonzero
537 : // (dk/dv)_e = 0 because e only depends on T and e is constant
538 3 : dk_dv = 0.0;
539 3 : dk_de = dk_dT / cv_from_p_T(p, T);
540 3 : }
541 :
542 : Real
543 447 : CaloricallyImperfectGas::s_from_v_e(Real v, Real e) const
544 : {
545 447 : Real T = T_from_v_e(v, e);
546 447 : Real Z = Z_from_T(T);
547 447 : return Z + _R_specific * std::log(v);
548 : }
549 :
550 : void
551 4 : CaloricallyImperfectGas::s_from_v_e(Real v, Real e, Real & s, Real & ds_dv, Real & ds_de) const
552 : {
553 4 : s = s_from_v_e(v, e);
554 : // see documentation for derivation
555 4 : ds_dv = _R_specific / v;
556 4 : ds_de = 1.0 / T_from_v_e(v, e);
557 4 : }
558 :
559 : Real
560 2690 : CaloricallyImperfectGas::s_from_p_T(Real p, Real T) const
561 : {
562 2690 : Real Z = Z_from_T(T);
563 2690 : Real v = 1.0 / rho_from_p_T(p, T);
564 2690 : return Z + _R_specific * std::log(v);
565 : }
566 :
567 : void
568 6 : CaloricallyImperfectGas::s_from_p_T(Real p, Real T, Real & s, Real & ds_dp, Real & ds_dT) const
569 : {
570 6 : s = s_from_p_T(p, T);
571 6 : ds_dp = -_R_specific / p;
572 6 : ds_dT = cp_from_p_T(p, T) / T;
573 6 : }
574 :
575 : Real
576 15 : CaloricallyImperfectGas::s_from_h_p(Real h, Real p) const
577 : {
578 15 : Real T = T_from_p_h(p, h);
579 15 : Real v = 1.0 / rho_from_p_T(p, T);
580 15 : Real Z = Z_from_T(T);
581 15 : return Z + _R_specific * std::log(v);
582 : }
583 :
584 : void
585 4 : CaloricallyImperfectGas::s_from_h_p(Real h, Real p, Real & s, Real & ds_dh, Real & ds_dp) const
586 : {
587 4 : s = s_from_h_p(h, p);
588 4 : ds_dh = 1.0 / T_from_p_h(p, h);
589 4 : ds_dp = -_R_specific / p;
590 4 : }
591 :
592 : Real
593 198 : CaloricallyImperfectGas::rho_from_p_s(Real p, Real s) const
594 : {
595 2036 : auto s_diff = [&p, &s, this](Real x) { return this->s_from_p_T(p, x) - s; };
596 198 : Real min = _min_temperature;
597 198 : Real max = _max_temperature;
598 198 : BrentsMethod::bracket(s_diff, min, max);
599 198 : Real T = BrentsMethod::root(s_diff, min, max);
600 198 : return rho_from_p_T(p, T);
601 : }
602 :
603 : Real
604 212 : CaloricallyImperfectGas::e_from_v_h(Real /*v*/, Real h) const
605 : {
606 212 : Real T = T_from_h(h);
607 212 : return e_from_T(T);
608 : }
609 :
610 : void
611 3 : CaloricallyImperfectGas::e_from_v_h(Real v, Real h, Real & e, Real & de_dv, Real & de_dh) const
612 : {
613 3 : e = e_from_v_h(v, h);
614 3 : Real cv = cv_from_v_e(v, e);
615 3 : Real cp = cp_from_v_e(v, e);
616 3 : de_dv = 0.0;
617 3 : de_dh = cv / cp;
618 3 : }
619 :
620 : Real
621 3316 : CaloricallyImperfectGas::rho_from_p_T(Real p, Real T) const
622 : {
623 3316 : return p * _molar_mass / (_R * T);
624 : }
625 :
626 : ADReal
627 198 : CaloricallyImperfectGas::rho_from_p_T(const ADReal & p, const ADReal & T) const
628 : {
629 198 : return p * _molar_mass / (_R * T);
630 : }
631 :
632 : void
633 0 : CaloricallyImperfectGas::rho_from_p_T(
634 : const ADReal & p, const ADReal & T, ADReal & rho, ADReal & drho_dp, ADReal & drho_dT) const
635 : {
636 0 : rho = rho_from_p_T(p, T);
637 0 : drho_dp = _molar_mass / (_R * T);
638 0 : drho_dT = -p * _molar_mass / (_R * T * T);
639 0 : }
640 :
641 : void
642 2 : CaloricallyImperfectGas::rho_from_p_T(
643 : Real p, Real T, Real & rho, Real & drho_dp, Real & drho_dT) const
644 : {
645 2 : rho = rho_from_p_T(p, T);
646 2 : drho_dp = _molar_mass / (_R * T);
647 2 : drho_dT = -p * _molar_mass / (_R * T * T);
648 2 : }
649 :
650 : Real
651 206 : CaloricallyImperfectGas::e_from_p_rho(Real p, Real rho) const
652 : {
653 206 : return e_from_p_rho_template(p, rho);
654 : }
655 :
656 : ADReal
657 198 : CaloricallyImperfectGas::e_from_p_rho(const ADReal & p, const ADReal & rho) const
658 : {
659 198 : return e_from_p_rho_template(p, rho);
660 : }
661 :
662 : void
663 2 : CaloricallyImperfectGas::e_from_p_rho(
664 : Real p, Real rho, Real & e, Real & de_dp, Real & de_drho) const
665 : {
666 2 : e_from_p_rho_template(p, rho, e, de_dp, de_drho);
667 2 : }
668 :
669 : void
670 0 : CaloricallyImperfectGas::e_from_p_rho(
671 : const ADReal & p, const ADReal & rho, ADReal & e, ADReal & de_dp, ADReal & de_drho) const
672 : {
673 0 : e_from_p_rho_template(p, rho, e, de_dp, de_drho);
674 0 : }
675 :
676 : Real
677 406 : CaloricallyImperfectGas::e_from_T_v(Real T, Real /*v*/) const
678 : {
679 406 : return e_from_T(T);
680 : }
681 :
682 : void
683 200 : CaloricallyImperfectGas::e_from_T_v(Real T, Real v, Real & e, Real & de_dT, Real & de_dv) const
684 : {
685 200 : e = e_from_T_v(T, v);
686 200 : de_dT = cv_from_T_v(T, v);
687 200 : de_dv = 0.0;
688 200 : }
689 :
690 : ADReal
691 198 : CaloricallyImperfectGas::e_from_T_v(const ADReal & T, const ADReal & v) const
692 : {
693 198 : Real x = 0;
694 198 : Real raw1 = T.value();
695 198 : Real raw2 = v.value();
696 198 : Real dxd1 = 0;
697 198 : Real dxd2 = 0;
698 198 : e_from_T_v(raw1, raw2, x, dxd1, dxd2);
699 :
700 198 : ADReal result = x;
701 198 : result.derivatives() = T.derivatives() * dxd1 + v.derivatives() * dxd2;
702 198 : return result;
703 : }
704 :
705 : Real
706 212 : CaloricallyImperfectGas::p_from_T_v(Real T, Real v) const
707 : {
708 212 : return _R_specific * T / v;
709 : }
710 :
711 : void
712 4 : CaloricallyImperfectGas::p_from_T_v(Real T, Real v, Real & p, Real & dp_dT, Real & dp_dv) const
713 : {
714 4 : p = _R_specific * T / v;
715 4 : dp_dT = _R_specific / v;
716 4 : dp_dv = -_R_specific * T / (v * v);
717 4 : }
718 :
719 : Real
720 203 : CaloricallyImperfectGas::h_from_T_v(Real T, Real /*v*/) const
721 : {
722 203 : return h_from_T(T);
723 : }
724 :
725 : void
726 2 : CaloricallyImperfectGas::h_from_T_v(Real T, Real /*v*/, Real & h, Real & dh_dT, Real & dh_dv) const
727 : {
728 2 : h = h_from_T(T);
729 2 : dh_dT = cp_from_T(T);
730 2 : dh_dv = 0.0;
731 2 : }
732 :
733 : Real
734 6 : CaloricallyImperfectGas::s_from_T_v(Real T, Real v) const
735 : {
736 6 : Real p = p_from_T_v(T, v);
737 6 : return s_from_p_T(p, T);
738 : }
739 :
740 : void
741 2 : CaloricallyImperfectGas::s_from_T_v(Real T, Real v, Real & s, Real & ds_dT, Real & ds_dv) const
742 : {
743 : Real p, dp_dT_v, dp_dv_T;
744 : Real ds_dp_T, ds_dT_p;
745 2 : p_from_T_v(T, v, p, dp_dT_v, dp_dv_T);
746 2 : s_from_p_T(p, T, s, ds_dp_T, ds_dT_p);
747 2 : ds_dT = ds_dT_p + ds_dp_T * dp_dT_v;
748 2 : ds_dv = ds_dp_T * dp_dv_T;
749 2 : }
750 :
751 : Real
752 398 : CaloricallyImperfectGas::cv_from_T_v(Real T, Real /*v*/) const
753 : {
754 398 : return cv_from_T(T);
755 : }
756 :
757 0 : Real CaloricallyImperfectGas::e_spndl_from_v(Real /*v*/) const { return _e_c; }
758 :
759 : void
760 0 : CaloricallyImperfectGas::v_e_spndl_from_T(Real /*T*/, Real & v, Real & e) const
761 : {
762 0 : v = 1. / _rho_c;
763 0 : e = _e_c;
764 0 : }
765 :
766 : Real
767 1223579 : CaloricallyImperfectGas::h_from_p_T(Real /*p*/, Real T) const
768 : {
769 1223579 : return h_from_T(T);
770 : }
771 :
772 : void
773 2 : CaloricallyImperfectGas::h_from_p_T(Real p, Real T, Real & h, Real & dh_dp, Real & dh_dT) const
774 : {
775 2 : h = h_from_p_T(p, T);
776 2 : dh_dp = 0.0;
777 2 : dh_dT = cp_from_p_T(p, T);
778 2 : }
779 :
780 : Real
781 1228710 : CaloricallyImperfectGas::e_from_p_T(Real /*p*/, Real T) const
782 : {
783 1228710 : return e_from_T(T);
784 : }
785 :
786 : ADReal
787 396 : CaloricallyImperfectGas::e_from_p_T(ADReal p, ADReal T) const
788 : {
789 396 : Real x = 0;
790 396 : Real raw1 = p.value();
791 396 : Real raw2 = T.value();
792 396 : Real dxd1 = 0;
793 396 : Real dxd2 = 0;
794 396 : e_from_p_T(raw1, raw2, x, dxd1, dxd2);
795 :
796 396 : ADReal result = x;
797 396 : result.derivatives() = p.derivatives() * dxd1 + T.derivatives() * dxd2;
798 396 : return result;
799 : }
800 :
801 : void
802 398 : CaloricallyImperfectGas::e_from_p_T(Real p, Real T, Real & e, Real & de_dp, Real & de_dT) const
803 : {
804 398 : e = e_from_p_T(p, T);
805 398 : de_dp = 0.0;
806 398 : de_dT = cv_from_p_T(p, T);
807 398 : }
808 :
809 : Real
810 0 : CaloricallyImperfectGas::molarMass() const
811 : {
812 0 : return _molar_mass;
813 : }
814 :
815 : Real
816 0 : CaloricallyImperfectGas::criticalTemperature() const
817 : {
818 0 : return _T_c;
819 : }
820 :
821 : Real
822 0 : CaloricallyImperfectGas::criticalDensity() const
823 : {
824 0 : return _rho_c;
825 : }
826 :
827 : Real
828 0 : CaloricallyImperfectGas::criticalInternalEnergy() const
829 : {
830 0 : return _e_c;
831 : }
832 :
833 : Real
834 224 : CaloricallyImperfectGas::T_from_p_h(Real /*p*/, Real h) const
835 : {
836 224 : return T_from_h(h);
837 : }
838 :
839 : void
840 2 : CaloricallyImperfectGas::T_from_p_h(Real p, Real h, Real & T, Real & dT_dp, Real & dT_dh) const
841 : {
842 2 : T = T_from_p_h(p, h);
843 2 : dT_dp = 0;
844 2 : dT_dh = 1.0 / cp_from_p_T(p, T);
845 2 : }
846 :
847 : Real
848 1029 : CaloricallyImperfectGas::cv_from_p_T(Real /* pressure */, Real temperature) const
849 : {
850 1029 : return cv_from_T(temperature);
851 : }
852 :
853 : void
854 996 : CaloricallyImperfectGas::cv_from_p_T(
855 : Real /*p*/, Real T, Real & cv, Real & dcv_dp, Real & dcv_dT) const
856 : {
857 996 : cv_from_T(T, cv, dcv_dT);
858 996 : dcv_dp = 0.0;
859 996 : }
860 :
861 : Real
862 438 : CaloricallyImperfectGas::cp_from_p_T(Real /* pressure */, Real temperature) const
863 : {
864 438 : return cp_from_T(temperature);
865 : }
866 :
867 : void
868 996 : CaloricallyImperfectGas::cp_from_p_T(
869 : Real /*p*/, Real T, Real & cp, Real & dcp_dp, Real & dcp_dT) const
870 : {
871 996 : cp_from_T(T, cp, dcp_dT);
872 996 : dcp_dp = 0.0;
873 996 : }
874 :
875 : Real
876 853 : CaloricallyImperfectGas::mu_from_p_T(Real /* pressure */, Real temperature) const
877 : {
878 853 : return _mu_T->value(temperature, Point());
879 : }
880 :
881 : void
882 6 : CaloricallyImperfectGas::mu_from_p_T(Real p, Real T, Real & mu, Real & dmu_dp, Real & dmu_dT) const
883 : {
884 6 : mu = this->mu_from_p_T(p, T);
885 6 : dmu_dp = 0.0;
886 : Real pert = 1.0e-7;
887 6 : Real mu2 = this->mu_from_p_T(p, T * (1 + pert));
888 6 : dmu_dT = (mu2 - mu) / (T * pert);
889 6 : }
890 :
891 : Real
892 853 : CaloricallyImperfectGas::k_from_p_T(Real /* pressure */, Real temperature) const
893 : {
894 853 : return _k_T->value(temperature, Point());
895 : }
896 :
897 : void
898 6 : CaloricallyImperfectGas::k_from_p_T(Real p, Real T, Real & k, Real & dk_dp, Real & dk_dT) const
899 : {
900 6 : k = k_from_p_T(p, T);
901 6 : dk_dp = 0.0;
902 : Real pert = 1.0e-7;
903 6 : Real k2 = this->k_from_p_T(p, T * (1 + pert));
904 6 : dk_dT = (k2 - k) / (T * pert);
905 6 : }
906 :
907 : Real
908 432 : CaloricallyImperfectGas::g_from_v_e(Real v, Real e) const
909 : {
910 : // use the definition of the Gibbs free energy
911 432 : Real T = T_from_v_e(v, e);
912 432 : Real p = p_from_v_e(v, e);
913 432 : return h_from_p_T(p, T) - T * s_from_p_T(p, T);
914 : }
915 :
916 : Real
917 0 : CaloricallyImperfectGas::p_from_h_s(Real h, Real s) const
918 : {
919 0 : Real T = T_from_h(h);
920 0 : Real e = e_from_T(T);
921 0 : Real Z = Z_from_T(T);
922 0 : Real v = std::exp((s - Z) / _R_specific);
923 0 : return p_from_v_e(v, e);
924 : }
925 :
926 : void
927 0 : CaloricallyImperfectGas::p_from_h_s(Real h, Real s, Real & p, Real & dp_dh, Real & dp_ds) const
928 : {
929 0 : p = p_from_h_s(h, s);
930 : Real pert = 1e-7;
931 0 : Real p_pert = p_from_h_s(h * (1 + pert), s);
932 0 : dp_dh = (p_pert - p) / (h * pert);
933 0 : p_pert = p_from_h_s(h, s * (1 + pert));
934 0 : dp_ds = (p_pert - p) / (s * pert);
935 0 : }
936 :
937 : void
938 2 : CaloricallyImperfectGas::outOfBounds(const std::string & function,
939 : Real value,
940 : Real min,
941 : Real max) const
942 : {
943 2 : std::stringstream ss;
944 4 : ss << "Function " << function << " encountered argument value of " << value
945 4 : << " which is outside of the bounds of " << min << " .. " << max;
946 2 : if (_out_of_bound_error)
947 2 : mooseError(ss.str());
948 : else
949 0 : mooseWarning(ss.str());
950 2 : }
|
