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 "GeneralFunctorFluidProps.h"
11 : #include "NS.h" // Variable Term Names
12 : #include "HeatTransferUtils.h"
13 : #include "NavierStokesMethods.h"
14 : #include "SinglePhaseFluidProperties.h"
15 :
16 : registerMooseObject("NavierStokesApp", GeneralFunctorFluidProps);
17 : registerMooseObject("NavierStokesApp", NonADGeneralFunctorFluidProps);
18 :
19 : template <bool is_ad>
20 : InputParameters
21 818 : GeneralFunctorFluidPropsTempl<is_ad>::validParams()
22 : {
23 818 : auto params = FunctorMaterial::validParams();
24 818 : params.addRequiredParam<UserObjectName>(NS::fluid, "Fluid properties object");
25 818 : params.addClassDescription("Creates functor fluid properties using a (P, T) formulation");
26 :
27 818 : params.addRequiredParam<MooseFunctorName>(NS::pressure, "Pressure");
28 818 : params.addRequiredParam<MooseFunctorName>(NS::T_fluid, "Fluid temperature");
29 818 : params.addRequiredParam<MooseFunctorName>(NS::speed, "Velocity norm");
30 818 : params.addParam<MooseFunctorName>(NS::density, "Density");
31 1636 : params.addParam<bool>(
32 : "force_define_density",
33 1636 : false,
34 : "Whether to force the definition of a density functor from the fluid properties");
35 1636 : params.addParam<bool>("neglect_derivatives_of_density_time_derivative",
36 1636 : true,
37 : "Whether to neglect the derivatives with regards to nonlinear variables "
38 : "of the density time derivatives");
39 :
40 1636 : params.addParam<FunctionName>(
41 1636 : "mu_rampdown", 1, "A function describing a ramp down of viscosity over time");
42 818 : params.addRequiredParam<MooseFunctorName>(NS::porosity, "porosity");
43 1636 : params.addRequiredParam<MooseFunctorName>(
44 : "characteristic_length", "characteristic length for Reynolds number calculation");
45 :
46 : // Dynamic pressure
47 1636 : params.addParam<bool>("solving_for_dynamic_pressure",
48 1636 : false,
49 : "Whether to solve for the dynamic pressure instead of the total pressure");
50 818 : params.addParam<Point>("reference_pressure_point",
51 818 : Point(0, 0, 0),
52 : "Point at which the gravity term for the static pressure is zero");
53 1636 : params.addParam<Real>("reference_pressure", 1e5, "Total pressure at the reference point");
54 1636 : params.addParam<Point>("gravity", Point(0, 0, -9.81), "Gravity vector");
55 :
56 1636 : params.addParamNamesToGroup(
57 : "solving_for_dynamic_pressure reference_pressure_point reference_pressure",
58 : "Dynamic pressure");
59 :
60 : // Property names
61 1636 : params.addParam<MooseFunctorName>(
62 : "density_name", NS::density, "Name to give to the density functor");
63 1636 : params.addParam<MooseFunctorName>(
64 : "dynamic_viscosity_name", NS::mu, "Name to give to the dynamic viscosity functor");
65 1636 : params.addParam<MooseFunctorName>(
66 : "specific_heat_name", NS::cp, "Name to give to the specific heat (cp) functor");
67 1636 : params.addParam<MooseFunctorName>(
68 : "thermal_conductivity_name", NS::k, "Name to give to the thermal conductivity functor");
69 1636 : params.addParamNamesToGroup(
70 : "density_name dynamic_viscosity_name specific_heat_name thermal_conductivity_name",
71 : "Functor property names");
72 818 : return params;
73 0 : }
74 :
75 : template <bool is_ad>
76 442 : GeneralFunctorFluidPropsTempl<is_ad>::GeneralFunctorFluidPropsTempl(
77 : const InputParameters & parameters)
78 : : FunctorMaterial(parameters),
79 884 : _fluid(UserObjectInterface::getUserObject<SinglePhaseFluidProperties>(NS::fluid)),
80 442 : _eps(getFunctor<GenericReal<is_ad>>(NS::porosity)),
81 884 : _d(getFunctor<Real>("characteristic_length")),
82 :
83 884 : _pressure_is_dynamic(getParam<bool>("solving_for_dynamic_pressure")),
84 884 : _reference_pressure_point(getParam<Point>("reference_pressure_point")),
85 884 : _reference_pressure_value(getParam<Real>("reference_pressure")),
86 884 : _gravity_vec(getParam<Point>("gravity")),
87 :
88 442 : _pressure(getFunctor<GenericReal<is_ad>>(NS::pressure)),
89 442 : _T_fluid(getFunctor<GenericReal<is_ad>>(NS::T_fluid)),
90 442 : _speed(getFunctor<GenericReal<is_ad>>(NS::speed)),
91 884 : _force_define_density(getParam<bool>("force_define_density")),
92 442 : _rho(getFunctor<GenericReal<is_ad>>(NS::density)),
93 442 : _mu_rampdown(getFunction("mu_rampdown")),
94 442 : _neglect_derivatives_of_density_time_derivative(
95 442 : getParam<bool>("neglect_derivatives_of_density_time_derivative")),
96 :
97 884 : _density_name(getParam<MooseFunctorName>("density_name")),
98 884 : _dynamic_viscosity_name(getParam<MooseFunctorName>("dynamic_viscosity_name")),
99 884 : _specific_heat_name(getParam<MooseFunctorName>("specific_heat_name")),
100 1326 : _thermal_conductivity_name(getParam<MooseFunctorName>("thermal_conductivity_name"))
101 : {
102 : // Check parameters
103 872 : if (!_pressure_is_dynamic &&
104 2162 : (isParamSetByUser("reference_pressure_point") || isParamSetByUser("reference_pressure")))
105 2 : paramError("solving_for_dynamic_pressure",
106 : "'reference_pressure_point' and 'reference_pressure' should not be set unless "
107 : "solving for the dynamic pressure");
108 :
109 : //
110 : // Set material properties functors
111 : //
112 :
113 : // If pressure is dynamic (else case), we must obtain the total pressure
114 : // We duplicate all this code. An alternative would be to redefine the pressure functor.
115 : // The issue with that is that you need the density functor to define the pressure,
116 : // and the pressure functor to define the density.
117 : // This could be solved by keeping a pointer to the pressure functor as an attribute and set
118 : // the pressure functor after the density functor has been defined.
119 440 : if (!_pressure_is_dynamic)
120 : {
121 428 : if (!isParamValid(NS::density) || _force_define_density)
122 1218 : addFunctorProperty<GenericReal<is_ad>>(
123 406 : _density_name,
124 41041755 : [this](const auto & r, const auto & t) -> GenericReal<is_ad>
125 41041755 : { return _fluid.rho_from_p_T(_pressure(r, t), _T_fluid(r, t)); });
126 :
127 1284 : addFunctorProperty<GenericReal<is_ad>>(
128 : NS::cv,
129 350 : [this](const auto & r, const auto & t) -> GenericReal<is_ad>
130 350 : { return _fluid.cv_from_p_T(_pressure(r, t), _T_fluid(r, t)); });
131 :
132 1284 : const auto & cp = addFunctorProperty<GenericReal<is_ad>>(
133 428 : _specific_heat_name,
134 11058863 : [this](const auto & r, const auto & t) -> GenericReal<is_ad>
135 11058863 : { return _fluid.cp_from_p_T(_pressure(r, t), _T_fluid(r, t)); });
136 :
137 1284 : const auto & mu = addFunctorProperty<GenericReal<is_ad>>(
138 428 : _dynamic_viscosity_name,
139 19894078 : [this](const auto & r, const auto & t) -> GenericReal<is_ad>
140 19894078 : { return _mu_rampdown(r, t) * _fluid.mu_from_p_T(_pressure(r, t), _T_fluid(r, t)); });
141 :
142 1284 : const auto & k = addFunctorProperty<GenericReal<is_ad>>(
143 428 : _thermal_conductivity_name,
144 11676793 : [this](const auto & r, const auto & t) -> GenericReal<is_ad>
145 11676793 : { return _fluid.k_from_p_T(_pressure(r, t), _T_fluid(r, t)); });
146 :
147 : //
148 : // Time derivatives of fluid properties
149 : //
150 428 : if (_neglect_derivatives_of_density_time_derivative)
151 : {
152 1624 : addFunctorProperty<GenericReal<is_ad>>(
153 406 : NS::time_deriv(_density_name),
154 2046990 : [this](const auto & r, const auto & t) -> GenericReal<is_ad>
155 : {
156 : Real rho, drho_dp, drho_dT;
157 2046990 : _fluid.rho_from_p_T(MetaPhysicL::raw_value(_pressure(r, t)),
158 2046990 : MetaPhysicL::raw_value(_T_fluid(r, t)),
159 : rho,
160 : drho_dp,
161 : drho_dT);
162 6140970 : return drho_dp * _pressure.dot(r, t) + drho_dT * _T_fluid.dot(r, t);
163 : });
164 : }
165 : else
166 : {
167 88 : addFunctorProperty<GenericReal<is_ad>>(
168 22 : NS::time_deriv(_density_name),
169 25600 : [this](const auto & r, const auto & t) -> GenericReal<is_ad>
170 : {
171 : GenericReal<is_ad> rho, drho_dp, drho_dT;
172 25600 : _fluid.rho_from_p_T(_pressure(r, t), _T_fluid(r, t), rho, drho_dp, drho_dT);
173 76800 : return drho_dp * _pressure.dot(r, t) + drho_dT * _T_fluid.dot(r, t);
174 : });
175 : }
176 :
177 1712 : addFunctorProperty<GenericReal<is_ad>>(
178 428 : NS::time_deriv(_specific_heat_name),
179 350 : [this](const auto & r, const auto & t) -> GenericReal<is_ad>
180 : {
181 : Real dcp_dp, dcp_dT, dummy;
182 350 : const auto raw_pressure = MetaPhysicL::raw_value(_pressure(r, t));
183 350 : const auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
184 350 : _fluid.cp_from_p_T(raw_pressure, raw_T_fluid, dummy, dcp_dp, dcp_dT);
185 :
186 700 : return dcp_dp * _pressure.dot(r, t) + dcp_dT * _T_fluid.dot(r, t);
187 : });
188 :
189 : //
190 : // Temperature and pressure derivatives, to help with computing time derivatives
191 : //
192 :
193 1712 : const auto & drho_dp = addFunctorProperty<Real>(
194 428 : derivativePropertyNameFirst(_density_name, NS::pressure),
195 700 : [this](const auto & r, const auto & t) -> Real
196 : {
197 : Real drho_dp, drho_dT, dummy;
198 700 : auto raw_pressure = MetaPhysicL::raw_value(_pressure(r, t));
199 700 : auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
200 :
201 700 : _fluid.rho_from_p_T(raw_pressure, raw_T_fluid, dummy, drho_dp, drho_dT);
202 700 : return drho_dp;
203 : });
204 :
205 1712 : const auto & drho_dT = addFunctorProperty<Real>(
206 428 : derivativePropertyNameFirst(_density_name, NS::T_fluid),
207 700 : [this](const auto & r, const auto & t) -> Real
208 : {
209 : Real drho_dp, drho_dT, dummy;
210 700 : auto raw_pressure = MetaPhysicL::raw_value(_pressure(r, t));
211 700 : auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
212 :
213 700 : _fluid.rho_from_p_T(raw_pressure, raw_T_fluid, dummy, drho_dp, drho_dT);
214 700 : return drho_dT;
215 : });
216 :
217 1712 : const auto & dcp_dp = addFunctorProperty<Real>(
218 428 : derivativePropertyNameFirst(_specific_heat_name, NS::pressure),
219 700 : [this](const auto & r, const auto & t) -> Real
220 : {
221 : Real dcp_dp, dcp_dT, dummy;
222 700 : auto raw_pressure = MetaPhysicL::raw_value(_pressure(r, t));
223 700 : auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
224 :
225 700 : _fluid.cp_from_p_T(raw_pressure, raw_T_fluid, dummy, dcp_dp, dcp_dT);
226 700 : return dcp_dp;
227 : });
228 :
229 1712 : const auto & dcp_dT = addFunctorProperty<Real>(
230 428 : derivativePropertyNameFirst(_specific_heat_name, NS::T_fluid),
231 700 : [this](const auto & r, const auto & t) -> Real
232 : {
233 : Real dcp_dp, dcp_dT, dummy;
234 700 : auto raw_pressure = MetaPhysicL::raw_value(_pressure(r, t));
235 700 : auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
236 :
237 700 : _fluid.cp_from_p_T(raw_pressure, raw_T_fluid, dummy, dcp_dp, dcp_dT);
238 700 : return dcp_dT;
239 : });
240 :
241 1712 : const auto & dmu_dp = addFunctorProperty<Real>(
242 428 : derivativePropertyNameFirst(_dynamic_viscosity_name, NS::pressure),
243 1050 : [this](const auto & r, const auto & t) -> Real
244 : {
245 : Real dmu_dp, dmu_dT, dummy;
246 1050 : auto raw_pressure = MetaPhysicL::raw_value(_pressure(r, t));
247 1050 : auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
248 :
249 1050 : _fluid.mu_from_p_T(raw_pressure, raw_T_fluid, dummy, dmu_dp, dmu_dT);
250 1050 : return _mu_rampdown(r, t) * dmu_dp;
251 : });
252 :
253 1712 : const auto & dmu_dT = addFunctorProperty<Real>(
254 428 : derivativePropertyNameFirst(_dynamic_viscosity_name, NS::T_fluid),
255 1050 : [this](const auto & r, const auto & t) -> Real
256 : {
257 : Real dmu_dp, dmu_dT, dummy;
258 1050 : auto raw_pressure = MetaPhysicL::raw_value(_pressure(r, t));
259 1050 : auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
260 :
261 1050 : _fluid.mu_from_p_T(raw_pressure, raw_T_fluid, dummy, dmu_dp, dmu_dT);
262 1050 : return _mu_rampdown(r, t) * dmu_dT;
263 : });
264 :
265 1712 : const auto & dk_dp = addFunctorProperty<Real>(
266 428 : derivativePropertyNameFirst(_thermal_conductivity_name, NS::pressure),
267 700 : [this](const auto & r, const auto & t) -> Real
268 : {
269 : Real dk_dp, dk_dT, dummy;
270 700 : auto raw_pressure = MetaPhysicL::raw_value(_pressure(r, t));
271 700 : auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
272 :
273 700 : _fluid.k_from_p_T(raw_pressure, raw_T_fluid, dummy, dk_dp, dk_dT);
274 700 : return dk_dp;
275 : });
276 :
277 1712 : const auto & dk_dT = addFunctorProperty<Real>(
278 428 : derivativePropertyNameFirst(_thermal_conductivity_name, NS::T_fluid),
279 700 : [this](const auto & r, const auto & t) -> Real
280 : {
281 : Real dk_dp, dk_dT, dummy;
282 700 : auto raw_pressure = MetaPhysicL::raw_value(_pressure(r, t));
283 700 : auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
284 :
285 700 : _fluid.k_from_p_T(raw_pressure, raw_T_fluid, dummy, dk_dp, dk_dT);
286 700 : return dk_dT;
287 : });
288 :
289 : //
290 : // Fluid adimensional quantities, used in numerous correlations
291 : //
292 :
293 1284 : addFunctorProperty<GenericReal<is_ad>>(
294 : NS::Prandtl,
295 5300 : [&cp, &mu, &k](const auto & r, const auto & t) -> GenericReal<is_ad>
296 : {
297 : static constexpr Real small_number = 1e-8;
298 :
299 : using std::max;
300 5300 : return HeatTransferUtils::prandtl(cp(r, t), mu(r, t), max(k(r, t), small_number));
301 : });
302 :
303 1284 : addFunctorProperty<Real>(
304 : derivativePropertyNameFirst(NS::Prandtl, NS::pressure),
305 350 : [&mu, &cp, &k, &dmu_dp, &dcp_dp, &dk_dp](const auto & r, const auto & t) -> Real
306 : {
307 350 : return NS::prandtlPropertyDerivative(MetaPhysicL::raw_value(mu(r, t)),
308 350 : MetaPhysicL::raw_value(cp(r, t)),
309 350 : MetaPhysicL::raw_value(k(r, t)),
310 350 : dmu_dp(r, t),
311 350 : dcp_dp(r, t),
312 700 : dk_dp(r, t));
313 : });
314 :
315 1284 : addFunctorProperty<Real>(
316 : derivativePropertyNameFirst(NS::Prandtl, NS::T_fluid),
317 350 : [&mu, &cp, &k, &dmu_dT, &dcp_dT, &dk_dT](const auto & r, const auto & t) -> Real
318 : {
319 350 : return NS::prandtlPropertyDerivative(MetaPhysicL::raw_value(mu(r, t)),
320 350 : MetaPhysicL::raw_value(cp(r, t)),
321 350 : MetaPhysicL::raw_value(k(r, t)),
322 350 : dmu_dT(r, t),
323 350 : dcp_dT(r, t),
324 700 : dk_dT(r, t));
325 : });
326 :
327 : //
328 : // (pore / particle) Reynolds number based on superficial velocity and
329 : // characteristic length. Only call Reynolds() one time to compute all three so that
330 : // we don't redundantly check that viscosity is not too close to zero.
331 : //
332 :
333 1284 : const auto & Re = addFunctorProperty<GenericReal<is_ad>>(
334 : NS::Reynolds,
335 6700 : [this, &mu](const auto & r, const auto & t) -> GenericReal<is_ad>
336 : {
337 : static constexpr Real small_number = 1e-8;
338 : using std::max;
339 : return max(
340 4950 : HeatTransferUtils::reynolds(
341 6700 : _rho(r, t), _eps(r, t) * _speed(r, t), _d(r, t), max(mu(r, t), small_number)),
342 5250 : small_number);
343 : });
344 :
345 1284 : addFunctorProperty<Real>(
346 : derivativePropertyNameFirst(NS::Reynolds, NS::pressure),
347 350 : [this, &Re, &mu, &drho_dp, &dmu_dp](const auto & r, const auto & t) -> Real
348 : {
349 350 : return NS::reynoldsPropertyDerivative(MetaPhysicL::raw_value(Re(r, t)),
350 350 : MetaPhysicL::raw_value(_rho(r, t)),
351 350 : MetaPhysicL::raw_value(mu(r, t)),
352 350 : drho_dp(r, t),
353 700 : dmu_dp(r, t));
354 : });
355 :
356 1284 : addFunctorProperty<Real>(
357 : derivativePropertyNameFirst(NS::Reynolds, NS::T_fluid),
358 350 : [this, &Re, &mu, &drho_dT, &dmu_dT](const auto & r, const auto & t) -> Real
359 : {
360 350 : return NS::reynoldsPropertyDerivative(MetaPhysicL::raw_value(Re(r, t)),
361 350 : MetaPhysicL::raw_value(_rho(r, t)),
362 350 : MetaPhysicL::raw_value(mu(r, t)),
363 350 : drho_dT(r, t),
364 700 : dmu_dT(r, t));
365 : });
366 :
367 : // (hydraulic) Reynolds number
368 1284 : addFunctorProperty<GenericReal<is_ad>>(
369 : NS::Reynolds_hydraulic,
370 350 : [this, &Re](const auto & r, const auto & t) -> GenericReal<is_ad>
371 : {
372 : static constexpr Real small_number = 1e-8;
373 : using std::max;
374 1050 : return Re(r, t) / max(1 - _eps(r, t), small_number);
375 : });
376 :
377 : // (interstitial) Reynolds number
378 1284 : addFunctorProperty<GenericReal<is_ad>>(
379 : NS::Reynolds_interstitial,
380 350 : [this, &Re](const auto & r, const auto & t) -> GenericReal<is_ad>
381 350 : { return Re(r, t) / _eps(r, t); });
382 : }
383 : else
384 : {
385 :
386 : const auto & rho =
387 0 : (!isParamValid(NS::density) || _force_define_density)
388 72 : ? addFunctorProperty<GenericReal<is_ad>>(
389 12 : _density_name,
390 8400 : [this](const auto & r, const auto & t) -> GenericReal<is_ad>
391 : {
392 8400 : auto total_pressure = _pressure(r, t) + _reference_pressure_value;
393 : // TODO: we should be integrating this term
394 8400 : const auto rho_approx = _fluid.rho_from_p_T(total_pressure, _T_fluid(r, t));
395 8400 : total_pressure +=
396 8400 : rho_approx * _gravity_vec * (r.getPoint() - _reference_pressure_point);
397 8400 : return _fluid.rho_from_p_T(total_pressure, _T_fluid(r, t));
398 : })
399 0 : : getFunctor<GenericReal<is_ad>>(_density_name);
400 :
401 36 : addFunctorProperty<GenericReal<is_ad>>(
402 : NS::cv,
403 175 : [this, &rho](const auto & r, const auto & t) -> GenericReal<is_ad>
404 : {
405 525 : const auto total_pressure =
406 175 : _pressure(r, t) + _reference_pressure_value +
407 175 : rho(r, t) * _gravity_vec * (r.getPoint() - _reference_pressure_point);
408 175 : return _fluid.cv_from_p_T(total_pressure, _T_fluid(r, t));
409 : });
410 :
411 36 : const auto & cp = addFunctorProperty<GenericReal<is_ad>>(
412 12 : _specific_heat_name,
413 700 : [this, &rho](const auto & r, const auto & t) -> GenericReal<is_ad>
414 : {
415 2100 : const auto total_pressure =
416 700 : _pressure(r, t) + _reference_pressure_value +
417 700 : rho(r, t) * _gravity_vec * (r.getPoint() - _reference_pressure_point);
418 700 : return _fluid.cp_from_p_T(total_pressure, _T_fluid(r, t));
419 : });
420 :
421 36 : const auto & mu = addFunctorProperty<GenericReal<is_ad>>(
422 12 : _dynamic_viscosity_name,
423 1925 : [this, &rho](const auto & r, const auto & t) -> GenericReal<is_ad>
424 : {
425 5775 : const auto total_pressure =
426 1925 : _pressure(r, t) + _reference_pressure_value +
427 1925 : rho(r, t) * _gravity_vec * (r.getPoint() - _reference_pressure_point);
428 1925 : return _mu_rampdown(r, t) * _fluid.mu_from_p_T(total_pressure, _T_fluid(r, t));
429 : });
430 :
431 36 : const auto & k = addFunctorProperty<GenericReal<is_ad>>(
432 12 : _thermal_conductivity_name,
433 700 : [this, &rho](const auto & r, const auto & t) -> GenericReal<is_ad>
434 : {
435 2100 : const auto total_pressure =
436 700 : _pressure(r, t) + _reference_pressure_value +
437 700 : rho(r, t) * _gravity_vec * (r.getPoint() - _reference_pressure_point);
438 700 : return _fluid.k_from_p_T(total_pressure, _T_fluid(r, t));
439 : });
440 :
441 : //
442 : // Time derivatives of fluid properties
443 : //
444 12 : if (_neglect_derivatives_of_density_time_derivative)
445 : {
446 48 : addFunctorProperty<GenericReal<is_ad>>(
447 12 : NS::time_deriv(_density_name),
448 175 : [this, &rho](const auto & r, const auto & t) -> GenericReal<is_ad>
449 : {
450 525 : const auto total_pressure =
451 175 : _pressure(r, t) + _reference_pressure_value +
452 175 : rho(r, t) * _gravity_vec * (r.getPoint() - _reference_pressure_point);
453 : Real rho, drho_dp, drho_dT;
454 175 : _fluid.rho_from_p_T(MetaPhysicL::raw_value(total_pressure),
455 175 : MetaPhysicL::raw_value(_T_fluid(r, t)),
456 : rho,
457 : drho_dp,
458 : drho_dT);
459 350 : return drho_dp * _pressure.dot(r, t) + drho_dT * _T_fluid.dot(r, t);
460 : });
461 : }
462 : else
463 : {
464 0 : addFunctorProperty<GenericReal<is_ad>>(
465 0 : NS::time_deriv(_density_name),
466 0 : [this, &rho](const auto & r, const auto & t) -> GenericReal<is_ad>
467 : {
468 0 : const auto total_pressure =
469 0 : _pressure(r, t) + _reference_pressure_value +
470 0 : rho(r, t) * _gravity_vec * (r.getPoint() - _reference_pressure_point);
471 : GenericReal<is_ad> rho, drho_dp, drho_dT;
472 0 : _fluid.rho_from_p_T(total_pressure, _T_fluid(r, t), rho, drho_dp, drho_dT);
473 0 : return drho_dp * _pressure.dot(r, t) + drho_dT * _T_fluid.dot(r, t);
474 : });
475 : }
476 :
477 48 : addFunctorProperty<GenericReal<is_ad>>(
478 12 : NS::time_deriv(_specific_heat_name),
479 175 : [this, &rho](const auto & r, const auto & t) -> GenericReal<is_ad>
480 : {
481 525 : const auto total_pressure =
482 175 : _pressure(r, t) + _reference_pressure_value +
483 350 : rho(r, t) * _gravity_vec * (r.getPoint() - _reference_pressure_point);
484 : Real dcp_dp, dcp_dT, dummy;
485 : auto raw_pressure = MetaPhysicL::raw_value(total_pressure);
486 175 : auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
487 175 : _fluid.cp_from_p_T(raw_pressure, raw_T_fluid, dummy, dcp_dp, dcp_dT);
488 :
489 350 : return dcp_dp * _pressure.dot(r, t) + dcp_dT * _T_fluid.dot(r, t);
490 : });
491 :
492 : //
493 : // Temperature and pressure derivatives, to help with computing time derivatives
494 : //
495 :
496 48 : const auto & drho_dp = addFunctorProperty<Real>(
497 12 : derivativePropertyNameFirst(_density_name, NS::pressure),
498 350 : [this, &rho](const auto & r, const auto & t) -> Real
499 : {
500 1050 : const auto total_pressure =
501 350 : _pressure(r, t) + _reference_pressure_value +
502 700 : rho(r, t) * _gravity_vec * (r.getPoint() - _reference_pressure_point);
503 : Real drho_dp, drho_dT, dummy;
504 : auto raw_pressure = MetaPhysicL::raw_value(total_pressure);
505 350 : auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
506 :
507 350 : _fluid.rho_from_p_T(raw_pressure, raw_T_fluid, dummy, drho_dp, drho_dT);
508 350 : return drho_dp;
509 : });
510 :
511 48 : const auto & drho_dT = addFunctorProperty<Real>(
512 12 : derivativePropertyNameFirst(_density_name, NS::T_fluid),
513 350 : [this, &rho](const auto & r, const auto & t) -> Real
514 : {
515 1050 : const auto total_pressure =
516 350 : _pressure(r, t) + _reference_pressure_value +
517 700 : rho(r, t) * _gravity_vec * (r.getPoint() - _reference_pressure_point);
518 : Real drho_dp, drho_dT, dummy;
519 : auto raw_pressure = MetaPhysicL::raw_value(total_pressure);
520 350 : auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
521 :
522 350 : _fluid.rho_from_p_T(raw_pressure, raw_T_fluid, dummy, drho_dp, drho_dT);
523 350 : return drho_dT;
524 : });
525 :
526 48 : const auto & dcp_dp = addFunctorProperty<Real>(
527 12 : derivativePropertyNameFirst(_specific_heat_name, NS::pressure),
528 350 : [this, &rho](const auto & r, const auto & t) -> Real
529 : {
530 1050 : const auto total_pressure =
531 350 : _pressure(r, t) + _reference_pressure_value +
532 700 : rho(r, t) * _gravity_vec * (r.getPoint() - _reference_pressure_point);
533 : Real dcp_dp, dcp_dT, dummy;
534 : auto raw_pressure = MetaPhysicL::raw_value(total_pressure);
535 350 : auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
536 :
537 350 : _fluid.cp_from_p_T(raw_pressure, raw_T_fluid, dummy, dcp_dp, dcp_dT);
538 350 : return dcp_dp;
539 : });
540 :
541 48 : const auto & dcp_dT = addFunctorProperty<Real>(
542 12 : derivativePropertyNameFirst(_specific_heat_name, NS::T_fluid),
543 350 : [this, &rho](const auto & r, const auto & t) -> Real
544 : {
545 1050 : const auto total_pressure =
546 350 : _pressure(r, t) + _reference_pressure_value +
547 700 : rho(r, t) * _gravity_vec * (r.getPoint() - _reference_pressure_point);
548 : Real dcp_dp, dcp_dT, dummy;
549 : auto raw_pressure = MetaPhysicL::raw_value(total_pressure);
550 350 : auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
551 :
552 350 : _fluid.cp_from_p_T(raw_pressure, raw_T_fluid, dummy, dcp_dp, dcp_dT);
553 350 : return dcp_dT;
554 : });
555 :
556 48 : const auto & dmu_dp = addFunctorProperty<Real>(
557 12 : derivativePropertyNameFirst(_dynamic_viscosity_name, NS::pressure),
558 525 : [this, &rho](const auto & r, const auto & t) -> Real
559 : {
560 1575 : const auto total_pressure =
561 525 : _pressure(r, t) + _reference_pressure_value +
562 1050 : rho(r, t) * _gravity_vec * (r.getPoint() - _reference_pressure_point);
563 : Real dmu_dp, dmu_dT, dummy;
564 : auto raw_pressure = MetaPhysicL::raw_value(total_pressure);
565 525 : auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
566 :
567 525 : _fluid.mu_from_p_T(raw_pressure, raw_T_fluid, dummy, dmu_dp, dmu_dT);
568 525 : return _mu_rampdown(r, t) * dmu_dp;
569 : });
570 :
571 48 : const auto & dmu_dT = addFunctorProperty<Real>(
572 12 : derivativePropertyNameFirst(_dynamic_viscosity_name, NS::T_fluid),
573 525 : [this, &rho](const auto & r, const auto & t) -> Real
574 : {
575 1575 : const auto total_pressure =
576 525 : _pressure(r, t) + _reference_pressure_value +
577 1050 : rho(r, t) * _gravity_vec * (r.getPoint() - _reference_pressure_point);
578 : Real dmu_dp, dmu_dT, dummy;
579 : auto raw_pressure = MetaPhysicL::raw_value(total_pressure);
580 525 : auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
581 :
582 525 : _fluid.mu_from_p_T(raw_pressure, raw_T_fluid, dummy, dmu_dp, dmu_dT);
583 525 : return _mu_rampdown(r, t) * dmu_dT;
584 : });
585 :
586 48 : const auto & dk_dp = addFunctorProperty<Real>(
587 12 : derivativePropertyNameFirst(_thermal_conductivity_name, NS::pressure),
588 350 : [this, &rho](const auto & r, const auto & t) -> Real
589 : {
590 1050 : const auto total_pressure =
591 350 : _pressure(r, t) + _reference_pressure_value +
592 700 : rho(r, t) * _gravity_vec * (r.getPoint() - _reference_pressure_point);
593 : Real dk_dp, dk_dT, dummy;
594 : auto raw_pressure = MetaPhysicL::raw_value(total_pressure);
595 350 : auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
596 :
597 350 : _fluid.k_from_p_T(raw_pressure, raw_T_fluid, dummy, dk_dp, dk_dT);
598 350 : return dk_dp;
599 : });
600 :
601 48 : const auto & dk_dT = addFunctorProperty<Real>(
602 12 : derivativePropertyNameFirst(_thermal_conductivity_name, NS::T_fluid),
603 350 : [this, &rho](const auto & r, const auto & t) -> Real
604 : {
605 1050 : const auto total_pressure =
606 350 : _pressure(r, t) + _reference_pressure_value +
607 700 : rho(r, t) * _gravity_vec * (r.getPoint() - _reference_pressure_point);
608 : Real dk_dp, dk_dT, dummy;
609 : auto raw_pressure = MetaPhysicL::raw_value(total_pressure);
610 350 : auto raw_T_fluid = MetaPhysicL::raw_value(_T_fluid(r, t));
611 :
612 350 : _fluid.k_from_p_T(raw_pressure, raw_T_fluid, dummy, dk_dp, dk_dT);
613 350 : return dk_dT;
614 : });
615 :
616 : //
617 : // Fluid adimensional quantities, used in numerous correlations
618 : //
619 :
620 36 : addFunctorProperty<GenericReal<is_ad>>(
621 : NS::Prandtl,
622 175 : [&cp, &mu, &k](const auto & r, const auto & t) -> GenericReal<is_ad>
623 : {
624 : static constexpr Real small_number = 1e-8;
625 : using std::max;
626 :
627 175 : return HeatTransferUtils::prandtl(cp(r, t), mu(r, t), max(k(r, t), small_number));
628 : });
629 :
630 36 : addFunctorProperty<Real>(
631 : derivativePropertyNameFirst(NS::Prandtl, NS::pressure),
632 175 : [&mu, &cp, &k, &dmu_dp, &dcp_dp, &dk_dp](const auto & r, const auto & t) -> Real
633 : {
634 175 : return NS::prandtlPropertyDerivative(MetaPhysicL::raw_value(mu(r, t)),
635 175 : MetaPhysicL::raw_value(cp(r, t)),
636 175 : MetaPhysicL::raw_value(k(r, t)),
637 175 : dmu_dp(r, t),
638 175 : dcp_dp(r, t),
639 350 : dk_dp(r, t));
640 : });
641 :
642 36 : addFunctorProperty<Real>(
643 : derivativePropertyNameFirst(NS::Prandtl, NS::T_fluid),
644 175 : [&mu, &cp, &k, &dmu_dT, &dcp_dT, &dk_dT](const auto & r, const auto & t) -> Real
645 : {
646 175 : return NS::prandtlPropertyDerivative(MetaPhysicL::raw_value(mu(r, t)),
647 175 : MetaPhysicL::raw_value(cp(r, t)),
648 175 : MetaPhysicL::raw_value(k(r, t)),
649 175 : dmu_dT(r, t),
650 175 : dcp_dT(r, t),
651 350 : dk_dT(r, t));
652 : });
653 :
654 : //
655 : // (pore / particle) Reynolds number based on superficial velocity and
656 : // characteristic length. Only call Reynolds() one time to compute all three so that
657 : // we don't redundantly check that viscosity is not too close to zero.
658 : //
659 :
660 36 : const auto & Re = addFunctorProperty<GenericReal<is_ad>>(
661 : NS::Reynolds,
662 875 : [this, &mu](const auto & r, const auto & t) -> GenericReal<is_ad>
663 : {
664 : static constexpr Real small_number = 1e-8;
665 : using std::max;
666 :
667 : return max(
668 0 : HeatTransferUtils::reynolds(
669 875 : _rho(r, t), _eps(r, t) * _speed(r, t), _d(r, t), max(mu(r, t), small_number)),
670 2625 : small_number);
671 : });
672 :
673 36 : addFunctorProperty<Real>(
674 : derivativePropertyNameFirst(NS::Reynolds, NS::pressure),
675 175 : [this, &Re, &mu, &drho_dp, &dmu_dp](const auto & r, const auto & t) -> Real
676 : {
677 175 : return NS::reynoldsPropertyDerivative(MetaPhysicL::raw_value(Re(r, t)),
678 175 : MetaPhysicL::raw_value(_rho(r, t)),
679 175 : MetaPhysicL::raw_value(mu(r, t)),
680 175 : drho_dp(r, t),
681 350 : dmu_dp(r, t));
682 : });
683 :
684 36 : addFunctorProperty<Real>(
685 : derivativePropertyNameFirst(NS::Reynolds, NS::T_fluid),
686 175 : [this, &Re, &mu, &drho_dT, &dmu_dT](const auto & r, const auto & t) -> Real
687 : {
688 175 : return NS::reynoldsPropertyDerivative(MetaPhysicL::raw_value(Re(r, t)),
689 175 : MetaPhysicL::raw_value(_rho(r, t)),
690 175 : MetaPhysicL::raw_value(mu(r, t)),
691 175 : drho_dT(r, t),
692 350 : dmu_dT(r, t));
693 : });
694 :
695 : // (hydraulic) Reynolds number
696 36 : addFunctorProperty<GenericReal<is_ad>>(
697 : NS::Reynolds_hydraulic,
698 175 : [this, &Re](const auto & r, const auto & t) -> GenericReal<is_ad>
699 : {
700 : static constexpr Real small_number = 1e-8;
701 : using std::max;
702 525 : return Re(r, t) / max(1 - _eps(r, t), small_number);
703 : });
704 :
705 : // (interstitial) Reynolds number
706 36 : addFunctorProperty<GenericReal<is_ad>>(
707 : NS::Reynolds_interstitial,
708 175 : [this, &Re](const auto & r, const auto & t) -> GenericReal<is_ad>
709 175 : { return Re(r, t) / _eps(r, t); });
710 : }
711 10538 : }
712 :
713 : template class GeneralFunctorFluidPropsTempl<false>;
714 : template class GeneralFunctorFluidPropsTempl<true>;
|
