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 "TemperaturePressureFunctionFluidProperties.h"
11 : #include "NewtonInversion.h"
12 :
13 : registerMooseObject("FluidPropertiesApp", TemperaturePressureFunctionFluidProperties);
14 :
15 : InputParameters
16 102 : TemperaturePressureFunctionFluidProperties::validParams()
17 : {
18 102 : InputParameters params = SinglePhaseFluidProperties::validParams();
19 204 : params.addRequiredParam<FunctionName>(
20 : "k", "Thermal conductivity function of temperature and pressure [W/(m-K)]");
21 204 : params.addRequiredParam<FunctionName>("rho",
22 : "Density function of temperature and pressure [kg/m^3]");
23 204 : params.addRequiredParam<FunctionName>(
24 : "mu", "Dynamic viscosity function of temperature and pressure [Pa-s]");
25 :
26 204 : params.addParam<FunctionName>(
27 : "cp", "Isobaric specific heat function of temperature and pressure [J/(kg-K)]");
28 204 : params.addRangeCheckedParam<Real>(
29 : "cv", 0, "cv >= 0", "Constant isochoric specific heat [J/(kg-K)]");
30 204 : params.addParam<Real>("e_ref", 0, "Specific internal energy at the reference temperature");
31 204 : params.addParam<Real>("T_ref", 0, "Reference temperature for the specific internal energy");
32 204 : params.addParam<Real>("dT_integration_intervals",
33 204 : 10,
34 : "Size of intervals for integrating cv(T) to compute e(T) from e(T_ref)");
35 :
36 102 : params.addClassDescription(
37 : "Single-phase fluid properties that allows to provide thermal "
38 : "conductivity, density, and viscosity as functions of temperature and pressure.");
39 102 : return params;
40 0 : }
41 :
42 54 : TemperaturePressureFunctionFluidProperties::TemperaturePressureFunctionFluidProperties(
43 54 : const InputParameters & parameters)
44 : : SinglePhaseFluidProperties(parameters),
45 54 : _initialized(false),
46 54 : _cv(getParam<Real>("cv")),
47 54 : _cv_is_constant(_cv != 0),
48 108 : _e_ref(getParam<Real>("e_ref")),
49 108 : _T_ref(getParam<Real>("T_ref")),
50 162 : _integration_dT(getParam<Real>("dT_integration_intervals"))
51 : {
52 162 : if (isParamValid("cp") && _cv_is_constant)
53 0 : paramError("cp", "The parameter 'cp' may only be specified if 'cv' is unspecified or is zero.");
54 54 : }
55 :
56 : void
57 49 : TemperaturePressureFunctionFluidProperties::initialSetup()
58 : {
59 49 : _k_function = &getFunction("k");
60 49 : _rho_function = &getFunction("rho");
61 49 : _mu_function = &getFunction("mu");
62 100 : _cp_function = isParamValid("cp") ? &getFunction("cp") : nullptr;
63 49 : _initialized = true;
64 49 : }
65 :
66 : std::string
67 1 : TemperaturePressureFunctionFluidProperties::fluidName() const
68 : {
69 1 : return "TemperaturePressureFunctionFluidProperties";
70 : }
71 :
72 : Real
73 52 : TemperaturePressureFunctionFluidProperties::T_from_v_e(Real v, Real e) const
74 : {
75 52 : if (_cv_is_constant)
76 48 : return _T_ref + (e - _e_ref) / _cv;
77 : else
78 : {
79 4 : const Real p0 = _p_initial_guess;
80 4 : const Real T0 = _T_initial_guess;
81 : Real p, T;
82 4 : bool conversion_succeeded = true;
83 4 : p_T_from_v_e(v, e, p0, T0, p, T, conversion_succeeded);
84 4 : if (conversion_succeeded)
85 4 : return T;
86 : else
87 0 : mooseError("T_from_v_e calculation failed.");
88 : }
89 : }
90 :
91 : Real
92 4 : TemperaturePressureFunctionFluidProperties::T_from_p_h(Real p, Real h) const
93 : {
94 : auto lambda = [&](Real p, Real current_T, Real & new_h, Real & dh_dp, Real & dh_dT)
95 11 : { h_from_p_T(p, current_T, new_h, dh_dp, dh_dT); };
96 4 : Real T = FluidPropertiesUtils::NewtonSolve(
97 4 : p, h, _T_initial_guess, _tolerance, lambda, name() + "::T_from_p_h", _max_newton_its)
98 4 : .first;
99 : // check for nans
100 4 : if (std::isnan(T))
101 0 : mooseError("Conversion from pressure (p = ",
102 : p,
103 : ") and enthalpy (h = ",
104 : h,
105 : ") to temperature failed to converge.");
106 4 : return T;
107 : }
108 :
109 : Real
110 8 : TemperaturePressureFunctionFluidProperties::T_from_p_rho(Real p, Real rho) const
111 : {
112 : auto lambda = [&](Real p, Real current_T, Real & new_rho, Real & drho_dp, Real & drho_dT)
113 24 : { rho_from_p_T(p, current_T, new_rho, drho_dp, drho_dT); };
114 8 : Real T = FluidPropertiesUtils::NewtonSolve(
115 8 : p, rho, _T_initial_guess, _tolerance, lambda, name() + "::T_from_p_rho")
116 8 : .first;
117 : // check for nans
118 8 : if (std::isnan(T))
119 0 : mooseError("Conversion from pressure (p = ",
120 : p,
121 : ") and density (rho = ",
122 : rho,
123 : ") to temperature failed to converge.");
124 8 : return T;
125 : }
126 :
127 : Real
128 36 : TemperaturePressureFunctionFluidProperties::cp_from_v_e(Real v, Real e) const
129 : {
130 36 : if (_cv_is_constant)
131 : {
132 18 : Real p = p_from_v_e(v, e);
133 18 : Real T = T_from_v_e(v, e);
134 18 : return cp_from_p_T(p, T);
135 : }
136 : else
137 : {
138 18 : const Real p0 = _p_initial_guess;
139 18 : const Real T0 = _T_initial_guess;
140 : Real p, T;
141 18 : bool conversion_succeeded = true;
142 18 : p_T_from_v_e(v, e, p0, T0, p, T, conversion_succeeded);
143 18 : if (conversion_succeeded)
144 18 : return cp_from_p_T(p, T);
145 : else
146 0 : mooseError("cp_from_v_e calculation failed. p= ", p, " T = ", T);
147 : }
148 : }
149 :
150 : void
151 4 : TemperaturePressureFunctionFluidProperties::cp_from_v_e(
152 : Real v, Real e, Real & cp, Real & dcp_dv, Real & dcp_de) const
153 : {
154 4 : cp = cp_from_v_e(v, e);
155 : // Using finite difference to get around difficulty of implementation
156 : Real eps = 1e-10;
157 4 : Real cp_pert = cp_from_v_e(v * (1 + eps), e);
158 4 : dcp_dv = (cp_pert - cp) / eps / v;
159 4 : cp_pert = cp_from_v_e(v, e * (1 + eps));
160 4 : dcp_de = (cp_pert - cp) / eps / e;
161 4 : }
162 :
163 : Real
164 26 : TemperaturePressureFunctionFluidProperties::cv_from_v_e(Real v, Real e) const
165 : {
166 26 : if (_cv_is_constant)
167 14 : return _cv;
168 : else
169 : {
170 12 : const Real p0 = _p_initial_guess;
171 12 : const Real T0 = _T_initial_guess;
172 : Real p, T;
173 12 : bool conversion_succeeded = true;
174 12 : p_T_from_v_e(v, e, p0, T0, p, T, conversion_succeeded);
175 12 : if (conversion_succeeded)
176 12 : return cv_from_p_T(p, T);
177 : else
178 0 : mooseError("cp_from_v_e calculation failed.");
179 : }
180 : }
181 :
182 : void
183 4 : TemperaturePressureFunctionFluidProperties::cv_from_v_e(
184 : Real v, Real e, Real & cv, Real & dcv_dv, Real & dcv_de) const
185 : {
186 4 : if (_cv_is_constant)
187 : {
188 2 : cv = cv_from_v_e(v, e);
189 2 : dcv_dv = 0.0;
190 2 : dcv_de = 0.0;
191 : }
192 : else
193 : {
194 2 : const Real p0 = _p_initial_guess;
195 2 : const Real T0 = _T_initial_guess;
196 : Real p, T;
197 2 : bool conversion_succeeded = true;
198 2 : p_T_from_v_e(v, e, p0, T0, p, T, conversion_succeeded);
199 : Real dcv_dp, dcv_dT;
200 2 : cv_from_p_T(p, T, cv, dcv_dp, dcv_dT);
201 2 : if (!conversion_succeeded)
202 0 : mooseError("cp_from_v_e and derivatives calculation failed.");
203 :
204 : Real p1, T1;
205 2 : p_T_from_v_e(v * (1 + 1e-6), e, p0, T0, p1, T1, conversion_succeeded);
206 2 : Real dp_dv = (p1 - p) / (v * 1e-6);
207 2 : Real dT_dv = (T1 - T) / (v * 1e-6);
208 2 : if (!conversion_succeeded)
209 0 : mooseError("cp_from_v_e and derivatives calculation failed.");
210 :
211 : Real p2, T2;
212 2 : p_T_from_v_e(v, e * (1 + 1e-6), p0, T0, p2, T2, conversion_succeeded);
213 2 : Real dp_de = (p2 - p) / (e * 1e-6);
214 2 : Real dT_de = (T2 - T) / (e * 1e-6);
215 2 : if (!conversion_succeeded)
216 0 : mooseError("cp_from_v_e and derivatives calculation failed.");
217 :
218 2 : dcv_dv = dcv_dp * dp_dv + dcv_dT * dT_dv;
219 2 : dcv_de = dcv_dp * dp_de + dcv_dT * dT_de;
220 : }
221 4 : }
222 :
223 : Real
224 26 : TemperaturePressureFunctionFluidProperties::p_from_v_e(Real v, Real e) const
225 : {
226 26 : const Real T = T_from_v_e(v, e);
227 : // note that p and T inversion in the definition of lambda
228 : auto lambda = [&](Real T, Real current_p, Real & new_rho, Real & drho_dT, Real & drho_dp)
229 78 : { rho_from_p_T(current_p, T, new_rho, drho_dp, drho_dT); };
230 26 : Real p = FluidPropertiesUtils::NewtonSolve(
231 26 : T, 1. / v, _p_initial_guess, _tolerance, lambda, name() + "::p_from_v_e")
232 26 : .first;
233 : // check for nans
234 26 : if (std::isnan(p))
235 0 : mooseError("Conversion from specific volume (v = ",
236 : v,
237 : ") and specific energy (e = ",
238 : e,
239 : ") to pressure failed to converge.");
240 26 : return p;
241 : }
242 :
243 : Real
244 4 : TemperaturePressureFunctionFluidProperties::mu_from_v_e(Real v, Real e) const
245 : {
246 4 : if (_cv_is_constant)
247 : {
248 2 : Real temperature = T_from_v_e(v, e);
249 2 : Real pressure = p_from_v_e(v, e);
250 2 : return mu_from_p_T(pressure, temperature);
251 : }
252 : else
253 : {
254 2 : const Real p0 = _p_initial_guess;
255 2 : const Real T0 = _T_initial_guess;
256 : Real p, T;
257 2 : bool conversion_succeeded = true;
258 2 : p_T_from_v_e(v, e, p0, T0, p, T, conversion_succeeded);
259 2 : if (conversion_succeeded)
260 2 : return mu_from_p_T(p, T);
261 : else
262 0 : mooseError("mu_from_v_e calculation failed.");
263 : }
264 : }
265 :
266 : Real
267 4 : TemperaturePressureFunctionFluidProperties::k_from_v_e(Real v, Real e) const
268 : {
269 4 : if (_cv_is_constant)
270 : {
271 2 : Real temperature = T_from_v_e(v, e);
272 2 : Real pressure = p_from_v_e(v, e);
273 2 : return k_from_p_T(pressure, temperature);
274 : }
275 : else
276 : {
277 2 : const Real p0 = _p_initial_guess;
278 2 : const Real T0 = _T_initial_guess;
279 : Real p, T;
280 2 : bool conversion_succeeded = true;
281 2 : p_T_from_v_e(v, e, p0, T0, p, T, conversion_succeeded);
282 2 : if (conversion_succeeded)
283 2 : return k_from_p_T(p, T);
284 : else
285 0 : mooseError("k_from_v_e calculation failed.");
286 : }
287 : }
288 :
289 : Real
290 18046812 : TemperaturePressureFunctionFluidProperties::rho_from_p_T(Real pressure, Real temperature) const
291 : {
292 18046812 : if (!_initialized)
293 4 : const_cast<TemperaturePressureFunctionFluidProperties *>(this)->initialSetup();
294 18046812 : return _rho_function->value(0, Point(temperature, pressure, 0));
295 : }
296 :
297 : void
298 18046341 : TemperaturePressureFunctionFluidProperties::rho_from_p_T(
299 : Real pressure, Real temperature, Real & rho, Real & drho_dp, Real & drho_dT) const
300 : {
301 18046341 : rho = rho_from_p_T(pressure, temperature);
302 18046341 : const RealVectorValue grad_function = _rho_function->gradient(0, Point(temperature, pressure, 0));
303 18046341 : drho_dT = grad_function(0);
304 18046341 : drho_dp = grad_function(1);
305 18046341 : }
306 :
307 : void
308 4 : TemperaturePressureFunctionFluidProperties::rho_from_p_T(const ADReal & pressure,
309 : const ADReal & temperature,
310 : ADReal & rho,
311 : ADReal & drho_dp,
312 : ADReal & drho_dT) const
313 : {
314 4 : rho = rho_from_p_T(pressure, temperature);
315 : const ADRealVectorValue grad_function =
316 4 : _rho_function->gradient(0, Point(temperature.value(), pressure.value(), 0));
317 4 : drho_dT = grad_function(0);
318 4 : drho_dp = grad_function(1);
319 4 : }
320 :
321 : Real
322 322 : TemperaturePressureFunctionFluidProperties::v_from_p_T(Real pressure, Real temperature) const
323 : {
324 322 : return 1.0 / rho_from_p_T(pressure, temperature);
325 : }
326 :
327 : void
328 298 : TemperaturePressureFunctionFluidProperties::v_from_p_T(
329 : Real pressure, Real temperature, Real & v, Real & dv_dp, Real & dv_dT) const
330 : {
331 298 : v = v_from_p_T(pressure, temperature);
332 :
333 : Real rho, drho_dp, drho_dT;
334 298 : rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT);
335 :
336 298 : dv_dp = -v * v * drho_dp;
337 298 : dv_dT = -v * v * drho_dT;
338 298 : }
339 :
340 : Real
341 73 : TemperaturePressureFunctionFluidProperties::h_from_p_T(Real pressure, Real temperature) const
342 : {
343 73 : Real e = e_from_p_T(pressure, temperature);
344 73 : Real rho = rho_from_p_T(pressure, temperature);
345 73 : return e + pressure / rho;
346 : }
347 :
348 : void
349 15 : TemperaturePressureFunctionFluidProperties::h_from_p_T(
350 : Real pressure, Real temperature, Real & h, Real & dh_dp, Real & dh_dT) const
351 : {
352 : Real e, de_dp, de_dT;
353 15 : e_from_p_T(pressure, temperature, e, de_dp, de_dT);
354 : Real rho, drho_dp, drho_dT;
355 15 : rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT);
356 15 : h = e + pressure / rho;
357 15 : dh_dp = de_dp + 1. / rho - pressure / rho / rho * drho_dp;
358 15 : dh_dT = de_dT - pressure * drho_dT / rho / rho;
359 15 : }
360 :
361 : Real
362 763 : TemperaturePressureFunctionFluidProperties::e_from_p_T(Real pressure, Real temperature) const
363 : {
364 763 : if (_cv_is_constant)
365 127 : return _e_ref + _cv * (temperature - _T_ref);
366 : else
367 : {
368 636 : const int n_intervals = std::ceil(std::abs(temperature - _T_ref) / _integration_dT);
369 636 : const auto h = (temperature - _T_ref) / n_intervals;
370 : Real integral = 0;
371 : // Centered step integration is second-order
372 18046148 : for (const auto i : make_range(n_intervals))
373 18045512 : integral += cv_from_p_T(pressure, _T_ref + (i + 0.5) * h);
374 636 : integral *= h;
375 : // we are still missing the dV or dP term to go from V/P_ref (e_ref = e(T_ref, V/P_ref))
376 : // to current V. The dT term is the largest one though
377 636 : return _e_ref + integral;
378 : }
379 : }
380 :
381 : void
382 315 : TemperaturePressureFunctionFluidProperties::e_from_p_T(
383 : Real pressure, Real temperature, Real & e, Real & de_dp, Real & de_dT) const
384 : {
385 315 : if (_cv_is_constant)
386 : {
387 11 : e = e_from_p_T(pressure, temperature);
388 11 : de_dp = 0.0;
389 11 : de_dT = _cv;
390 : }
391 : else
392 : {
393 304 : e = e_from_p_T(pressure, temperature);
394 304 : Real ep = e_from_p_T(pressure * (1 + 1e-8), temperature);
395 304 : de_dp = (ep - e) / (pressure * 1e-8);
396 304 : de_dT = cv_from_p_T(pressure, temperature);
397 : }
398 315 : }
399 :
400 : Real
401 4 : TemperaturePressureFunctionFluidProperties::e_from_p_rho(Real p, Real rho) const
402 : {
403 4 : if (_cv_is_constant)
404 : {
405 2 : Real temperature = T_from_p_rho(p, rho);
406 2 : return _e_ref + _cv * (temperature - _T_ref);
407 : }
408 : else
409 : {
410 2 : Real temperature = T_from_p_rho(p, rho);
411 2 : return e_from_p_T(p, temperature);
412 : }
413 : }
414 :
415 : Real
416 0 : TemperaturePressureFunctionFluidProperties::beta_from_p_T(Real pressure, Real temperature) const
417 : {
418 : Real rho, drho_dp, drho_dT;
419 0 : rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT);
420 0 : return -drho_dT / rho;
421 : }
422 :
423 : Real
424 18045965 : TemperaturePressureFunctionFluidProperties::cp_from_p_T(Real p, Real T) const
425 : {
426 18045965 : if (_cv_is_constant)
427 : {
428 : Real rho, drho_dp, drho_dT;
429 81 : rho_from_p_T(p, T, rho, drho_dp, drho_dT);
430 : // Wikipedia notation for thermal expansion / compressibility coefficients
431 81 : Real alpha = -drho_dT / rho;
432 81 : Real beta = -drho_dp / rho;
433 81 : return _cv + MathUtils::pow(alpha, 2) * T / rho / beta;
434 : }
435 : else
436 : {
437 18045884 : if (!_initialized)
438 1 : const_cast<TemperaturePressureFunctionFluidProperties *>(this)->initialSetup();
439 18045884 : return _cp_function->value(0, Point(T, p, 0));
440 : }
441 : }
442 :
443 : void
444 4 : TemperaturePressureFunctionFluidProperties::cp_from_p_T(
445 : Real pressure, Real temperature, Real & cp, Real & dcp_dp, Real & dcp_dT) const
446 : {
447 4 : if (_cv_is_constant)
448 : {
449 2 : cp = cp_from_p_T(pressure, temperature);
450 : Real eps = 1e-8;
451 2 : Real cp_1p = cp_from_p_T(pressure * (1 + eps), temperature);
452 2 : Real cp_1T = cp_from_p_T(pressure, temperature * (1 + eps));
453 2 : dcp_dp = (cp_1p - cp) / (pressure * eps);
454 2 : dcp_dT = (cp_1T - cp) / (temperature * eps);
455 : }
456 : else
457 : {
458 2 : cp = cp_from_p_T(pressure, temperature);
459 : const RealVectorValue grad_function =
460 2 : _cp_function->gradient(0, Point(temperature, pressure, 0));
461 2 : dcp_dT = grad_function(0);
462 2 : dcp_dp = grad_function(1);
463 : }
464 4 : }
465 :
466 : Real
467 18045911 : TemperaturePressureFunctionFluidProperties::cv_from_p_T(Real pressure, Real temperature) const
468 : {
469 18045911 : if (_cv_is_constant)
470 59 : return _cv;
471 : else
472 : {
473 : Real rho, drho_dp, drho_dT;
474 18045852 : rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT);
475 : // Wikipedia notation for thermal expansion / compressibility coefficients
476 18045852 : Real alpha = -drho_dT / rho;
477 18045852 : Real beta = -drho_dp / rho;
478 18045852 : return cp_from_p_T(pressure, temperature) - MathUtils::pow(alpha, 2) * temperature / rho / beta;
479 : }
480 : }
481 :
482 : void
483 6 : TemperaturePressureFunctionFluidProperties::cv_from_p_T(
484 : Real pressure, Real temperature, Real & cv, Real & dcv_dp, Real & dcv_dT) const
485 : {
486 6 : if (_cv_is_constant)
487 : {
488 2 : cv = cv_from_p_T(pressure, temperature);
489 2 : dcv_dp = 0.0;
490 2 : dcv_dT = 0.0;
491 : }
492 : else
493 : {
494 4 : cv = cv_from_p_T(pressure, temperature);
495 : Real eps = 1e-10;
496 4 : Real cv_1p = cv_from_p_T(pressure * (1 + eps), temperature);
497 4 : Real cv_1T = cv_from_p_T(pressure, temperature * (1 + eps));
498 4 : dcv_dp = (cv_1p - cv) / (pressure * eps);
499 4 : dcv_dT = (cv_1T - cv) / (temperature * eps);
500 : }
501 6 : }
502 :
503 : Real
504 80 : TemperaturePressureFunctionFluidProperties::mu_from_p_T(Real pressure, Real temperature) const
505 : {
506 80 : if (!_initialized)
507 0 : const_cast<TemperaturePressureFunctionFluidProperties *>(this)->initialSetup();
508 80 : return _mu_function->value(0, Point(temperature, pressure, 0));
509 : }
510 :
511 : void
512 6 : TemperaturePressureFunctionFluidProperties::mu_from_p_T(
513 : Real pressure, Real temperature, Real & mu, Real & dmu_dp, Real & dmu_dT) const
514 : {
515 6 : mu = mu_from_p_T(pressure, temperature);
516 6 : const RealVectorValue grad_function = _mu_function->gradient(0, Point(temperature, pressure, 0));
517 6 : dmu_dT = grad_function(0);
518 6 : dmu_dp = grad_function(1);
519 6 : }
520 :
521 : Real
522 81 : TemperaturePressureFunctionFluidProperties::k_from_p_T(Real pressure, Real temperature) const
523 : {
524 81 : if (!_initialized)
525 0 : const_cast<TemperaturePressureFunctionFluidProperties *>(this)->initialSetup();
526 81 : return _k_function->value(0, Point(temperature, pressure, 0));
527 : }
528 :
529 : void
530 4 : TemperaturePressureFunctionFluidProperties::k_from_p_T(
531 : Real pressure, Real temperature, Real & k, Real & dk_dp, Real & dk_dT) const
532 : {
533 4 : k = k_from_p_T(pressure, temperature);
534 4 : const RealVectorValue grad_function = _k_function->gradient(0, Point(temperature, pressure, 0));
535 4 : dk_dT = grad_function(0);
536 4 : dk_dp = grad_function(1);
537 4 : }
|
