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 29832 : 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 2815441 : Reynolds(const T1 & volume_fraction, const T2 & rho, const T3 & vel, const T4 & D_h, const T5 & mu)
119 : {
120 2815450 : 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 927161 : 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 17177 : 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 17177 : const auto diffusivity = HeatTransferUtils::thermalDiffusivity(k, rho, cp);
161 17177 : 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 1 : return gravity_magnitude * beta * dT * std::pow(D_h, 3) * (rho_liquid * rho_liquid) /
186 1 : (mu_liquid * mu_liquid);
187 : }
188 :
189 : /**
190 : * Compute Laplace number (or coefficient)
191 : *
192 : * @param surf_tension Surface tension
193 : * @param delta_rho Difference in density of phases
194 : * @param gravity_magnitude Gravitational acceleration magnitude
195 : *
196 : * @return Laplace number
197 : */
198 : template <typename T1, typename T2>
199 : auto
200 : Laplace(const T1 & surf_tension, const T2 & delta_rho, const Real & gravity_magnitude)
201 : {
202 : return std::sqrt(surf_tension / (gravity_magnitude * delta_rho));
203 : }
204 :
205 : /**
206 : * Compute viscosity number (or coefficient)
207 : *
208 : * @param viscosity Viscosity
209 : * @param surf_tension Surface tension
210 : * @param rho_k Density of k-th phase of interest
211 : * @param delta_rho Density difference
212 : * @param gravity_magnitude Gravitational acceleration magnitude
213 : *
214 : * @return viscosity number
215 : */
216 : template <typename T1, typename T2, typename T3, typename T4>
217 : auto
218 1 : viscosityNumber(const T1 & viscosity,
219 : const T2 & surf_tension,
220 : const T3 & rho_k,
221 : const T4 & delta_rho,
222 : const Real & gravity_magnitude)
223 : {
224 1 : return viscosity /
225 1 : std::sqrt(rho_k * surf_tension * std::sqrt(surf_tension / gravity_magnitude / delta_rho));
226 : }
227 :
228 : using NS::wallHeatTransferCoefficient;
229 :
230 : /**
231 : * Compute Dean number
232 : * @param Re Reynolds number
233 : * @param doD tube diameter to coil diameter ratio
234 : *
235 : * @return Dean number
236 : */
237 : template <typename T1, typename T2>
238 : auto
239 : Dean(const T1 & Re, const T2 & doD)
240 : {
241 : return Re * std::sqrt(doD);
242 : }
243 :
244 : /**
245 : * Computes velocity and its derivatives from alpha*rho*A and alpha*rho*u*A
246 : *
247 : * @param[in] arhoA alpha*rho*A
248 : * @param[in] arhouA alpha*rho*u*A
249 : * @param[out] vel velocity
250 : * @param[out] dvel_darhoA derivative of velocity w.r.t. alpha*rho*A
251 : * @param[out] dvel_darhouA derivative of velocity w.r.t. alpha*rho*u*A
252 : */
253 : void
254 : vel_from_arhoA_arhouA(Real arhoA, Real arhouA, Real & vel, Real & dvel_darhoA, Real & dvel_darhouA);
255 :
256 : /**
257 : * Computes velocity from alpha*rho*A and alpha*rho*u*A
258 : *
259 : * @param arhoA alpha*rho*A
260 : * @param arhouA alpha*rho*u*A
261 : * @return velocity
262 : */
263 : ADReal vel_from_arhoA_arhouA(ADReal arhoA, ADReal arhouA);
264 :
265 : /**
266 : * Derivative of velocity w.r.t. alpha*rho*A
267 : *
268 : * @param[in] arhoA alpha*rho*A
269 : * @param[in] arhouA alpha*rho*u*A
270 : * @returns derivative of velocity w.r.t. alpha*rho*A
271 : */
272 : Real dvel_darhoA(Real arhoA, Real arhouA);
273 :
274 : /**
275 : * Derivative of velocity w.r.t. alpha*rho*u*A
276 : *
277 : * @param[in] arhoA alpha*rho*A
278 : * @returns derivative of velocity w.r.t. alpha*rho*u*A
279 : */
280 : Real dvel_darhouA(Real arhoA);
281 :
282 : /**
283 : * Computes density and its derivatives from alpha*rho*A, alpha, and area.
284 : *
285 : * @param[in] arhoA alpha*rho*A
286 : * @param[in] alpha volume fraction
287 : * @param[in] A area
288 : * @param[out] rho density
289 : * @param[out] drho_darhoA derivative of density w.r.t. alpha*rho*A
290 : * @param[out] drho_dalpha derivative of density w.r.t. alpha
291 : */
292 : void rho_from_arhoA_alpha_A(
293 : Real arhoA, Real alpha, Real A, Real & rho, Real & drho_darhoA, Real & drho_dalpha);
294 :
295 : /**
296 : * Computes density from alpha*rho*A, alpha, and area.
297 : *
298 : * @param[in] arhoA alpha*rho*A
299 : * @param[in] alpha volume fraction
300 : * @param[in] A area
301 : * @returns density
302 : */
303 : ADReal rho_from_arhoA_alpha_A(ADReal arhoA, ADReal alpha, ADReal A);
304 :
305 : /**
306 : * Computes specific volume and its derivatives from rho*A, and area.
307 : *
308 : * @param[in] rhoA rho*A
309 : * @param[in] A area
310 : * @param[out] dv_drhoA derivative of specific volume w.r.t. rho*A
311 : */
312 : void v_from_rhoA_A(Real rhoA, Real A, Real & v, Real & dv_drhoA);
313 :
314 : /**
315 : * Computes specific volume and its derivatives from rho*A, and area.
316 : *
317 : * @param[in] rhoA rho*A
318 : * @param[in] A area
319 : * @returns specific volume
320 : */
321 : ADReal v_from_rhoA_A(ADReal rhoA, ADReal A);
322 :
323 : /**
324 : * Computes specific volume and its derivatives from alpha*rho*A, volume fraction, and area.
325 : *
326 : * @param[in] arhoA alpha*rho*A
327 : * @param[in] alpha volume fraction
328 : * @param[in] A area
329 : * @param[out] dv_darhoA derivative of specific volume w.r.t. alpha*rho*A
330 : * @param[out] dv_dalpha derivative of specific volume w.r.t. volume fraction
331 : */
332 : void
333 : v_from_arhoA_alpha_A(Real arhoA, Real alpha, Real A, Real & v, Real & dv_darhoA, Real & dv_dalpha);
334 :
335 : /**
336 : * Computes specific volume and its derivatives from alpha*rho*A, volume fraction, and area.
337 : *
338 : * @param[in] arhoA alpha*rho*A
339 : * @param[in] alpha volume fraction
340 : * @param[in] A area
341 : * @returns specific volume
342 : */
343 : ADReal v_from_arhoA_alpha_A(ADReal arhoA, ADReal alpha, ADReal A);
344 :
345 : /**
346 : * Computes specific volume and its derivative with respect to density
347 : *
348 : * @param[in] rho density
349 : * @param[in] v specific volume
350 : * @param[in] dv_drho derivative of specific volume w.r.t. density
351 : */
352 : void v_from_rho(Real rho, Real & v, Real & dv_drho);
353 :
354 : /**
355 : * Derivative of specific volume wrt alpha_liquid
356 : *
357 : * Makes sense only when using 7-equation model
358 : * @param area - The cross-sectional area
359 : * @param arhoA - density equation solution variable: alpha*rho*A
360 : * @param is_liquid - True if the specific volume corresponds to liquid phase
361 : */
362 : Real dv_dalpha_liquid(Real area, Real arhoA, bool is_liquid);
363 :
364 : /**
365 : * Derivative of specific volume wrt density equation solution variable
366 : *
367 : * @param area - Cross-sectional area
368 : * @param arhoA - density equation solution variable: alpha*rho*A
369 : */
370 : Real dv_darhoA(Real area, Real arhoA);
371 :
372 : /**
373 : * Computes specific internal energy and its derivatives from alpha*rho*A, alpha*rho*u*A, and
374 : * alpha*rho*E*A
375 : *
376 : * @param[in] arhoA alpha*rho*A
377 : * @param[in] arhouA alpha*rho*u*A
378 : * @param[in] arhoEA alpha*rho*E*A
379 : * @param[out] e specific internal energy
380 : * @param[out] de_darhoA derivative of specific internal energy w.r.t. alpha*rho*A
381 : * @param[out] de_darhouA derivative of specific internal energy w.r.t. alpha*rho*u*A
382 : * @param[out] de_darhoEA derivative of specific internal energy w.r.t. alpha*rho*E*A
383 : */
384 : void e_from_arhoA_arhouA_arhoEA(Real arhoA,
385 : Real arhouA,
386 : Real arhoEA,
387 : Real & e,
388 : Real & de_darhoA,
389 : Real & de_darhouA,
390 : Real & de_darhoEA);
391 :
392 : ADReal e_from_arhoA_arhouA_arhoEA(ADReal arhoA, ADReal arhouA, ADReal arhoEA);
393 :
394 : /**
395 : * Computes specific internal energy and its derivatives from specific total energy and velocity
396 : *
397 : * @param[in] E specific total energy
398 : * @param[in] vel velocity
399 : * @param[out] e specific internal energy
400 : * @param[out] de_dE derivative of specific internal energy w.r.t. specific total energy
401 : * @param[out] de_dvel derivative of specific internal energy w.r.t. velocity
402 : */
403 : void e_from_E_vel(Real E, Real vel, Real & e, Real & de_dE, Real & de_dvel);
404 :
405 : /**
406 : * Computes specific internal energy from specific total energy and velocity
407 : *
408 : * @param[in] E specific total energy
409 : * @param[in] vel velocity
410 : * @returns specific internal energy
411 : */
412 : ADReal e_from_E_vel(ADReal E, ADReal vel);
413 :
414 : /**
415 : * Derivative of specific internal energy wrt density of the phase (rhoA or arhoA)
416 : *
417 : * @param arhoA - density equation solution variable: alpha*rho*A
418 : * @param arhouA - momentum equation solution variable: alpha*rho*u*A
419 : * @param arhoEA - energy equation solution variable: alpha*rho*E*A
420 : */
421 : Real de_darhoA(Real arhoA, Real arhouA, Real arhoEA);
422 :
423 : /**
424 : * Derivative of specific internal energy wrt momentum of the phase (rhouA or arhouA)
425 : *
426 : * @param arhoA - density equation solution variable: alpha*rho*A
427 : * @param arhouA - momentum equation solution variable: alpha*rho*u*A
428 : */
429 : Real de_darhouA(Real arhoA, Real arhouA);
430 :
431 : /**
432 : * Derivative of specific internal energy wrt total energy of the phase (rhoEA or arhoEA)
433 : *
434 : * @param arhoA - density equation solution variable: alpha*rho*A
435 : */
436 : Real de_darhoEA(Real arhoA);
437 :
438 : /**
439 : * Computes specific total energy and its derivatives from alpha*rho*A and alpha*rho*E*A
440 : *
441 : * @param[in] arhoA alpha*rho*A
442 : * @param[in] arhoEA alpha*rho*E*A
443 : * @param[out] E specific total energy
444 : * @param[out] dE_darhoA derivative of specific total energy w.r.t. alpha*rho*A
445 : * @param[out] dE_darhoEA derivative of specific total energy w.r.t. alpha*rho*E*A
446 : */
447 : void E_from_arhoA_arhoEA(Real arhoA, Real arhoEA, Real & E, Real & dE_darhoA, Real & dE_darhoEA);
448 :
449 : /**
450 : * Computes specific total energy from alpha*rho*A and alpha*rho*E*A
451 : *
452 : * @param[in] arhoA alpha*rho*A
453 : * @param[in] arhoEA alpha*rho*E*A
454 : * @returns specific total energy
455 : */
456 : ADReal E_from_arhoA_arhoEA(ADReal arhoA, ADReal arhoEA);
457 :
458 : /**
459 : * Computes specific total energy and its derivatives from specific internal energy and velocity
460 : *
461 : * @param[in] e specific internal energy
462 : * @param[in] vel velocity
463 : * @param[out] E specific total energy
464 : * @param[out] dE_de derivative of specific total energy w.r.t. specific internal energy
465 : * @param[out] dE_dvel derivative of specific total energy w.r.t. velocity
466 : */
467 : void E_from_e_vel(Real e, Real vel, Real & E, Real & dE_de, Real & dE_dvel);
468 :
469 : /**
470 : * Computes specific enthalpy and its derivatives from specific internal energy, pressure, and
471 : * density
472 : *
473 : * @param[in] e specific internal energy
474 : * @param[in] p pressure
475 : * @param[in] rho density
476 : * @param[out] h specific enthalpy
477 : * @param[out] dh_de derivative of specific enthalpy w.r.t. specific internal energy
478 : * @param[out] dh_dp derivative of specific enthalpy w.r.t. pressure
479 : * @param[out] dh_drho derivative of specific enthalpy w.r.t. density
480 : */
481 : void h_from_e_p_rho(Real e, Real p, Real rho, Real & h, Real & dh_de, Real & dh_dp, Real & dh_drho);
482 :
483 : ADReal h_from_e_p_rho(ADReal e, ADReal p, ADReal rho);
484 :
485 : /**
486 : * Determine if inlet boundary condition should be applied
487 : *
488 : * @return true if the flow conditions are inlet, false otherwise
489 : * @param vel Velocity of the phase
490 : * @param normal Outward normal vector
491 : */
492 : bool isInlet(Real vel, Real normal);
493 : bool isInlet(ADReal vel, Real normal);
494 :
495 : /**
496 : * Determine if outlet boundary condition should be applied
497 : *
498 : * @return true if the flow conditions are outlet, false otherwise
499 : * @param vel Velocity of the phase
500 : * @param normal Outward normal vector
501 : */
502 : bool isOutlet(Real vel, Real normal);
503 : bool isOutlet(ADReal vel, Real normal);
504 : }
|
