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 : #pragma once
11 :
12 : #include "libmesh/libmesh_common.h"
13 : #include "libmesh/vector_value.h"
14 : #include "libmesh/dense_vector.h"
15 :
16 : #include "ADReal.h"
17 : #include "NavierStokesMethods.h"
18 : #include "HeatTransferUtils.h"
19 :
20 : using namespace libMesh;
21 :
22 : namespace THM
23 : {
24 :
25 : // Default value for magnitude of acceleration due to gravity
26 : static const Real gravity_const = 9.81;
27 :
28 : // Default value for gravitational acceleration vector
29 : static VectorValue<Real> default_gravity_vector = VectorValue<Real>(0.0, 0.0, -gravity_const);
30 :
31 : // Stefan-Boltzman constant, in [W/m^2-K]
32 : static const Real Stefan_Boltzman_const = 5.670e-8;
33 :
34 : /**
35 : * The sign function
36 : * @param val The argument of the sign function
37 : * @return -1 for negative values, 0 for zero and 1 for positive values
38 : */
39 : template <typename T>
40 : int
41 : sgn(T val)
42 : {
43 : return (T(0) < val) - (val < T(0));
44 : }
45 :
46 : /**
47 : * Tests if two real-valued vectors are equal within some absolute tolerance
48 : *
49 : * @param[in] a First vector
50 : * @param[in] b Second vector
51 : * @param[in] tol Absolute tolerance
52 : */
53 : bool absoluteFuzzyEqualVectors(const RealVectorValue & a,
54 : const RealVectorValue & b,
55 2 : const Real & tol = libMesh::TOLERANCE * libMesh::TOLERANCE);
56 :
57 : /**
58 : * Tests if two real-valued vectors are parallel within some absolute tolerance
59 : *
60 : * @param[in] a First vector
61 : * @param[in] b Second vector
62 : * @param[in] tol Absolute tolerance
63 : */
64 : bool areParallelVectors(const RealVectorValue & a,
65 : const RealVectorValue & b,
66 11680 : const Real & tol = libMesh::TOLERANCE * libMesh::TOLERANCE);
67 :
68 : /**
69 : * Tests if two real-valued vectors are in the same direction
70 : *
71 : * @param[in] a First vector
72 : * @param[in] b Second vector
73 : * @param[in] tol Absolute tolerance
74 : */
75 : bool haveSameDirection(const RealVectorValue & a,
76 : const RealVectorValue & b,
77 2 : const Real & tol = libMesh::TOLERANCE * libMesh::TOLERANCE);
78 :
79 : /**
80 : * Computes a derivative of a fraction using quotient rule for a derivative
81 : * w.r.t. a scalar quantity
82 : *
83 : * @param[in] num numerator value
84 : * @param[in] den denominator value
85 : * @param[in] dnum_dy derivative of numerator value
86 : * @param[in] dden_dy derivative of denominator value
87 : */
88 : Real
89 : applyQuotientRule(const Real & num, const Real & den, const Real & dnum_dy, const Real & dden_dy);
90 :
91 : /**
92 : * Computes a derivative of a fraction using quotient rule for a derivative
93 : * w.r.t. a vector quantity
94 : *
95 : * @param[in] num numerator value
96 : * @param[in] den denominator value
97 : * @param[in] dnum_dy derivative of numerator value
98 : * @param[in] dden_dy derivative of denominator value
99 : */
100 : DenseVector<Real> applyQuotientRule(const Real & num,
101 : const Real & den,
102 : const DenseVector<Real> & dnum_dy,
103 : const DenseVector<Real> & dden_dy);
104 :
105 : /**
106 : * Compute Reynolds number
107 : *
108 : * @param volume_fraction The volume fraction of the phase
109 : * @param rho The density of the phase
110 : * @param vel The velocity of the phase
111 : * @param D_h The hydraulic diameter
112 : * @param mu The viscosity of the phase
113 : *
114 : * @return Reynolds number
115 : */
116 : template <typename T1, typename T2, typename T3, typename T4, typename T5>
117 : auto
118 40617 : Reynolds(const T1 & volume_fraction, const T2 & rho, const T3 & vel, const T4 & D_h, const T5 & mu)
119 : {
120 40623 : return volume_fraction * HeatTransferUtils::reynolds(rho, vel, D_h, mu);
121 : }
122 :
123 : /**
124 : * Compute Prandtl number
125 : * @param cp Specific heat
126 : * @param mu Dynamic viscosity
127 : * @param k Thermal conductivity
128 : *
129 : * @return Prandtl number
130 : */
131 : template <typename T1, typename T2, typename T3>
132 : auto
133 : Prandtl(const T1 & cp, const T2 & mu, const T3 & k)
134 : {
135 15530 : return HeatTransferUtils::prandtl(cp, mu, k);
136 : }
137 :
138 : /**
139 : * Compute Peclet number
140 : *
141 : * @param volume_fraction The volume fraction of the phase
142 : * @param rho The density of the phase
143 : * @param vel The velocity of the phase
144 : * @param D_h The hydraulic diameter
145 : * @param k Thermal conductivity
146 : * @param cp Specific heat
147 : * @param k Thermal conductivity
148 : *
149 : * @return Peclet number
150 : */
151 : template <typename T1, typename T2, typename T3, typename T4, typename T5, typename T6>
152 : auto
153 12362 : Peclet(const T1 & volume_fraction,
154 : const T2 & cp,
155 : const T3 & rho,
156 : const T4 & vel,
157 : const T5 & D_h,
158 : const T6 & k)
159 : {
160 12362 : const auto diffusivity = HeatTransferUtils::thermalDiffusivity(k, rho, cp);
161 12362 : return volume_fraction * HeatTransferUtils::peclet(vel, D_h, diffusivity);
162 : }
163 :
164 : /**
165 : * Compute Grashof number
166 : *
167 : * @param beta Thermal expansion coefficient
168 : * @param dT |T_w - T|
169 : * @param D_h Hydraulic diameter
170 : * @param rho_liquid Density of liquid
171 : * @param mu_liquid Viscosity of liquid
172 : * @param gravity_magnitude Gravitational acceleration magnitude
173 : *
174 : * @return Grashof number
175 : */
176 : template <typename T1, typename T2, typename T3, typename T4, typename T5>
177 : auto
178 1 : Grashof(const T1 & beta,
179 : const T2 & dT,
180 : const T3 & D_h,
181 : const T4 & rho_liquid,
182 : const T5 & mu_liquid,
183 : const Real & gravity_magnitude)
184 : {
185 : using std::pow;
186 1 : return gravity_magnitude * beta * dT * pow(D_h, 3) * (rho_liquid * rho_liquid) /
187 1 : (mu_liquid * mu_liquid);
188 : }
189 :
190 : /**
191 : * Compute Laplace number (or coefficient)
192 : *
193 : * @param surf_tension Surface tension
194 : * @param delta_rho Difference in density of phases
195 : * @param gravity_magnitude Gravitational acceleration magnitude
196 : *
197 : * @return Laplace number
198 : */
199 : template <typename T1, typename T2>
200 : auto
201 : Laplace(const T1 & surf_tension, const T2 & delta_rho, const Real & gravity_magnitude)
202 : {
203 : using std::sqrt;
204 : return sqrt(surf_tension / (gravity_magnitude * delta_rho));
205 : }
206 :
207 : /**
208 : * Compute viscosity number (or coefficient)
209 : *
210 : * @param viscosity Viscosity
211 : * @param surf_tension Surface tension
212 : * @param rho_k Density of k-th phase of interest
213 : * @param delta_rho Density difference
214 : * @param gravity_magnitude Gravitational acceleration magnitude
215 : *
216 : * @return viscosity number
217 : */
218 : template <typename T1, typename T2, typename T3, typename T4>
219 : auto
220 1 : viscosityNumber(const T1 & viscosity,
221 : const T2 & surf_tension,
222 : const T3 & rho_k,
223 : const T4 & delta_rho,
224 : const Real & gravity_magnitude)
225 : {
226 : using std::sqrt;
227 1 : return viscosity /
228 1 : sqrt(rho_k * surf_tension * sqrt(surf_tension / gravity_magnitude / delta_rho));
229 : }
230 :
231 : using NS::wallHeatTransferCoefficient;
232 :
233 : /**
234 : * Compute Dean number
235 : * @param Re Reynolds number
236 : * @param doD tube diameter to coil diameter ratio
237 : *
238 : * @return Dean number
239 : */
240 : template <typename T1, typename T2>
241 : auto
242 : Dean(const T1 & Re, const T2 & doD)
243 : {
244 : using std::sqrt;
245 : return Re * sqrt(doD);
246 : }
247 :
248 : /**
249 : * Computes velocity and its derivatives from alpha*rho*A and alpha*rho*u*A
250 : *
251 : * @param[in] arhoA alpha*rho*A
252 : * @param[in] arhouA alpha*rho*u*A
253 : * @param[out] vel velocity
254 : * @param[out] dvel_darhoA derivative of velocity w.r.t. alpha*rho*A
255 : * @param[out] dvel_darhouA derivative of velocity w.r.t. alpha*rho*u*A
256 : */
257 : void
258 : vel_from_arhoA_arhouA(Real arhoA, Real arhouA, Real & vel, Real & dvel_darhoA, Real & dvel_darhouA);
259 :
260 : /**
261 : * Computes velocity from alpha*rho*A and alpha*rho*u*A
262 : *
263 : * @param arhoA alpha*rho*A
264 : * @param arhouA alpha*rho*u*A
265 : * @return velocity
266 : */
267 : ADReal vel_from_arhoA_arhouA(ADReal arhoA, ADReal arhouA);
268 :
269 : /**
270 : * Derivative of velocity w.r.t. alpha*rho*A
271 : *
272 : * @param[in] arhoA alpha*rho*A
273 : * @param[in] arhouA alpha*rho*u*A
274 : * @returns derivative of velocity w.r.t. alpha*rho*A
275 : */
276 : Real dvel_darhoA(Real arhoA, Real arhouA);
277 :
278 : /**
279 : * Derivative of velocity w.r.t. alpha*rho*u*A
280 : *
281 : * @param[in] arhoA alpha*rho*A
282 : * @returns derivative of velocity w.r.t. alpha*rho*u*A
283 : */
284 : Real dvel_darhouA(Real arhoA);
285 :
286 : /**
287 : * Computes density and its derivatives from alpha*rho*A, alpha, and area.
288 : *
289 : * @param[in] arhoA alpha*rho*A
290 : * @param[in] alpha volume fraction
291 : * @param[in] A area
292 : * @param[out] rho density
293 : * @param[out] drho_darhoA derivative of density w.r.t. alpha*rho*A
294 : * @param[out] drho_dalpha derivative of density w.r.t. alpha
295 : */
296 : void rho_from_arhoA_alpha_A(
297 : Real arhoA, Real alpha, Real A, Real & rho, Real & drho_darhoA, Real & drho_dalpha);
298 :
299 : /**
300 : * Computes density from alpha*rho*A, alpha, and area.
301 : *
302 : * @param[in] arhoA alpha*rho*A
303 : * @param[in] alpha volume fraction
304 : * @param[in] A area
305 : * @returns density
306 : */
307 : ADReal rho_from_arhoA_alpha_A(ADReal arhoA, ADReal alpha, ADReal A);
308 :
309 : /**
310 : * Computes specific volume and its derivatives from rho*A, and area.
311 : *
312 : * @param[in] rhoA rho*A
313 : * @param[in] A area
314 : * @param[out] dv_drhoA derivative of specific volume w.r.t. rho*A
315 : */
316 : void v_from_rhoA_A(Real rhoA, Real A, Real & v, Real & dv_drhoA);
317 :
318 : /**
319 : * Computes specific volume and its derivatives from rho*A, and area.
320 : *
321 : * @param[in] rhoA rho*A
322 : * @param[in] A area
323 : * @returns specific volume
324 : */
325 : ADReal v_from_rhoA_A(ADReal rhoA, ADReal A);
326 :
327 : /**
328 : * Computes specific volume and its derivatives from alpha*rho*A, volume fraction, and area.
329 : *
330 : * @param[in] arhoA alpha*rho*A
331 : * @param[in] alpha volume fraction
332 : * @param[in] A area
333 : * @param[out] dv_darhoA derivative of specific volume w.r.t. alpha*rho*A
334 : * @param[out] dv_dalpha derivative of specific volume w.r.t. volume fraction
335 : */
336 : void
337 : v_from_arhoA_alpha_A(Real arhoA, Real alpha, Real A, Real & v, Real & dv_darhoA, Real & dv_dalpha);
338 :
339 : /**
340 : * Computes specific volume and its derivatives from alpha*rho*A, volume fraction, and area.
341 : *
342 : * @param[in] arhoA alpha*rho*A
343 : * @param[in] alpha volume fraction
344 : * @param[in] A area
345 : * @returns specific volume
346 : */
347 : ADReal v_from_arhoA_alpha_A(ADReal arhoA, ADReal alpha, ADReal A);
348 :
349 : /**
350 : * Computes specific volume and its derivative with respect to density
351 : *
352 : * @param[in] rho density
353 : * @param[in] v specific volume
354 : * @param[in] dv_drho derivative of specific volume w.r.t. density
355 : */
356 : void v_from_rho(Real rho, Real & v, Real & dv_drho);
357 :
358 : /**
359 : * Derivative of specific volume wrt alpha_liquid
360 : *
361 : * Makes sense only when using 7-equation model
362 : * @param area - The cross-sectional area
363 : * @param arhoA - density equation solution variable: alpha*rho*A
364 : * @param is_liquid - True if the specific volume corresponds to liquid phase
365 : */
366 : Real dv_dalpha_liquid(Real area, Real arhoA, bool is_liquid);
367 :
368 : /**
369 : * Derivative of specific volume wrt density equation solution variable
370 : *
371 : * @param area - Cross-sectional area
372 : * @param arhoA - density equation solution variable: alpha*rho*A
373 : */
374 : Real dv_darhoA(Real area, Real arhoA);
375 :
376 : /**
377 : * Computes specific internal energy and its derivatives from alpha*rho*A, alpha*rho*u*A, and
378 : * alpha*rho*E*A
379 : *
380 : * @param[in] arhoA alpha*rho*A
381 : * @param[in] arhouA alpha*rho*u*A
382 : * @param[in] arhoEA alpha*rho*E*A
383 : * @param[out] e specific internal energy
384 : * @param[out] de_darhoA derivative of specific internal energy w.r.t. alpha*rho*A
385 : * @param[out] de_darhouA derivative of specific internal energy w.r.t. alpha*rho*u*A
386 : * @param[out] de_darhoEA derivative of specific internal energy w.r.t. alpha*rho*E*A
387 : */
388 : void e_from_arhoA_arhouA_arhoEA(Real arhoA,
389 : Real arhouA,
390 : Real arhoEA,
391 : Real & e,
392 : Real & de_darhoA,
393 : Real & de_darhouA,
394 : Real & de_darhoEA);
395 :
396 : ADReal e_from_arhoA_arhouA_arhoEA(ADReal arhoA, ADReal arhouA, ADReal arhoEA);
397 :
398 : /**
399 : * Computes specific internal energy and its derivatives from specific total energy and velocity
400 : *
401 : * @param[in] E specific total energy
402 : * @param[in] vel velocity
403 : * @param[out] e specific internal energy
404 : * @param[out] de_dE derivative of specific internal energy w.r.t. specific total energy
405 : * @param[out] de_dvel derivative of specific internal energy w.r.t. velocity
406 : */
407 : void e_from_E_vel(Real E, Real vel, Real & e, Real & de_dE, Real & de_dvel);
408 :
409 : /**
410 : * Computes specific internal energy from specific total energy and velocity
411 : *
412 : * @param[in] E specific total energy
413 : * @param[in] vel velocity
414 : * @returns specific internal energy
415 : */
416 : ADReal e_from_E_vel(ADReal E, ADReal vel);
417 :
418 : /**
419 : * Derivative of specific internal energy wrt density of the phase (rhoA or arhoA)
420 : *
421 : * @param arhoA - density equation solution variable: alpha*rho*A
422 : * @param arhouA - momentum equation solution variable: alpha*rho*u*A
423 : * @param arhoEA - energy equation solution variable: alpha*rho*E*A
424 : */
425 : Real de_darhoA(Real arhoA, Real arhouA, Real arhoEA);
426 :
427 : /**
428 : * Derivative of specific internal energy wrt momentum of the phase (rhouA or arhouA)
429 : *
430 : * @param arhoA - density equation solution variable: alpha*rho*A
431 : * @param arhouA - momentum equation solution variable: alpha*rho*u*A
432 : */
433 : Real de_darhouA(Real arhoA, Real arhouA);
434 :
435 : /**
436 : * Derivative of specific internal energy wrt total energy of the phase (rhoEA or arhoEA)
437 : *
438 : * @param arhoA - density equation solution variable: alpha*rho*A
439 : */
440 : Real de_darhoEA(Real arhoA);
441 :
442 : /**
443 : * Computes specific total energy and its derivatives from alpha*rho*A and alpha*rho*E*A
444 : *
445 : * @param[in] arhoA alpha*rho*A
446 : * @param[in] arhoEA alpha*rho*E*A
447 : * @param[out] E specific total energy
448 : * @param[out] dE_darhoA derivative of specific total energy w.r.t. alpha*rho*A
449 : * @param[out] dE_darhoEA derivative of specific total energy w.r.t. alpha*rho*E*A
450 : */
451 : void E_from_arhoA_arhoEA(Real arhoA, Real arhoEA, Real & E, Real & dE_darhoA, Real & dE_darhoEA);
452 :
453 : /**
454 : * Computes specific total energy from alpha*rho*A and alpha*rho*E*A
455 : *
456 : * @param[in] arhoA alpha*rho*A
457 : * @param[in] arhoEA alpha*rho*E*A
458 : * @returns specific total energy
459 : */
460 : ADReal E_from_arhoA_arhoEA(ADReal arhoA, ADReal arhoEA);
461 :
462 : /**
463 : * Computes specific total energy and its derivatives from specific internal energy and velocity
464 : *
465 : * @param[in] e specific internal energy
466 : * @param[in] vel velocity
467 : * @param[out] E specific total energy
468 : * @param[out] dE_de derivative of specific total energy w.r.t. specific internal energy
469 : * @param[out] dE_dvel derivative of specific total energy w.r.t. velocity
470 : */
471 : void E_from_e_vel(Real e, Real vel, Real & E, Real & dE_de, Real & dE_dvel);
472 :
473 : /**
474 : * Computes specific enthalpy and its derivatives from specific internal energy, pressure, and
475 : * density
476 : *
477 : * @param[in] e specific internal energy
478 : * @param[in] p pressure
479 : * @param[in] rho density
480 : * @param[out] h specific enthalpy
481 : * @param[out] dh_de derivative of specific enthalpy w.r.t. specific internal energy
482 : * @param[out] dh_dp derivative of specific enthalpy w.r.t. pressure
483 : * @param[out] dh_drho derivative of specific enthalpy w.r.t. density
484 : */
485 : void h_from_e_p_rho(Real e, Real p, Real rho, Real & h, Real & dh_de, Real & dh_dp, Real & dh_drho);
486 :
487 : ADReal h_from_e_p_rho(ADReal e, ADReal p, ADReal rho);
488 :
489 : /**
490 : * Determine if inlet boundary condition should be applied
491 : *
492 : * @return true if the flow conditions are inlet, false otherwise
493 : * @param vel Velocity of the phase
494 : * @param normal Outward normal vector
495 : */
496 : bool isInlet(Real vel, Real normal);
497 : bool isInlet(ADReal vel, Real normal);
498 :
499 : /**
500 : * Determine if outlet boundary condition should be applied
501 : *
502 : * @return true if the flow conditions are outlet, false otherwise
503 : * @param vel Velocity of the phase
504 : * @param normal Outward normal vector
505 : */
506 : bool isOutlet(Real vel, Real normal);
507 : bool isOutlet(ADReal vel, Real normal);
508 : }
|
