9.2
general documentation
cs_turbulence_bc.h
Go to the documentation of this file.
1#ifndef CS_TURBULENCE_BC_H
2#define CS_TURBULENCE_BC_H
3
4/*============================================================================
5 * Base turbulence boundary conditions.
6 *============================================================================*/
7
8/*
9 This file is part of code_saturne, a general-purpose CFD tool.
10
11 Copyright (C) 1998-2026 EDF S.A.
12
13 This program is free software; you can redistribute it and/or modify it under
14 the terms of the GNU General Public License as published by the Free Software
15 Foundation; either version 2 of the License, or (at your option) any later
16 version.
17
18 This program is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
21 details.
22
23 You should have received a copy of the GNU General Public License along with
24 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
25 Street, Fifth Floor, Boston, MA 02110-1301, USA.
26*/
27
28/*----------------------------------------------------------------------------*/
29
30/*----------------------------------------------------------------------------
31 * Local headers
32 *----------------------------------------------------------------------------*/
33
34#include "base/cs_defs.h"
35
36/*=============================================================================
37 * Macro definitions
38 *============================================================================*/
39
40/*============================================================================
41 * Type definitions
42 *============================================================================*/
43
44/*============================================================================
45 * Global variables
46 *============================================================================*/
47
48/*=============================================================================
49 * Public function prototypes
50 *============================================================================*/
51
52/*----------------------------------------------------------------------------*/
53/*
54 * Initialize turbulence model boundary condition pointers.
55 */
56/*----------------------------------------------------------------------------*/
57
58extern "C" void
59cs_turbulence_bc_init_pointers(void);
60
61/*----------------------------------------------------------------------------*/
62/*
63 * \brief Free memory allocations for turbulence boundary condition pointers.
64 */
65/*----------------------------------------------------------------------------*/
66
67extern "C" void
68cs_turbulence_bc_free_pointers(void);
69
70/*----------------------------------------------------------------------------*/
71/*
72 * \brief Calculation of \f$ u^\star \f$, \f$ k \f$ and \f$\varepsilon \f$
73 * from a diameter \f$ D_H \f$ and the reference velocity \f$ U_{ref} \f$
74 * for a circular duct flow with smooth wall
75 * (use for inlet boundary conditions).
76 *
77 * Both \f$ u^\star \f$ and\f$ (k,\varepsilon )\f$ are returned, so that
78 * the user may compute other values of \f$ k \f$ and \f$ \varepsilon \f$
79 * with \f$ u^\star \f$.
80 *
81 * We use the laws from Idel'Cik, i.e.
82 * the head loss coefficient \f$ \lambda \f$ is defined by:
83 * \f[ |\dfrac{\Delta P}{\Delta x}| =
84 * \dfrac{\lambda}{D_H} \frac{1}{2} \rho U_{ref}^2 \f]
85 *
86 * then the relation reads \f$u^\star = U_{ref} \sqrt{\dfrac{\lambda}{8}}\f$.
87 * \f$\lambda \f$ depends on the hydraulic Reynolds number
88 * \f$ Re = \dfrac{U_{ref} D_H}{ \nu} \f$ and is given by:
89 * - for \f$ Re < 2000 \f$
90 * \f[ \lambda = \dfrac{64}{Re} \f]
91 *
92 * - for \f$ Re > 4000 \f$
93 * \f[ \lambda = \dfrac{1}{( 1.8 \log_{10}(Re)-1.64 )^2} \f]
94 *
95 * - for \f$ 2000 < Re < 4000 \f$, we complete by a straight line
96 * \f[ \lambda = 0.021377 + 5.3115. 10^{-6} Re \f]
97 *
98 * From \f$ u^\star \f$, we can estimate \f$ k \f$ and \f$ \varepsilon\f$
99 * from the well known formulae of developped turbulence
100 *
101 * \f[ k = \dfrac{u^{\star 2}}{\sqrt{C_\mu}} \f]
102 * \f[ \varepsilon = \dfrac{ u^{\star 3}}{(\kappa D_H /10)} \f]
103 *
104 * \param[in] uref2 square of the reference flow velocity
105 * \param[in] dh hydraulic diameter \f$ D_H \f$
106 * \param[in] rho mass density \f$ \rho \f$
107 * \param[in] mu dynamic viscosity \f$ \nu \f$
108 * \param[out] ustar2 square of friction speed
109 * \param[out] k calculated turbulent intensity \f$ k \f$
110 * \param[out] eps calculated turbulent dissipation
111 * \f$ \varepsilon \f$
112 */
113/*----------------------------------------------------------------------------*/
114
115extern "C" void
116cs_turbulence_bc_ke_hyd_diam(double uref2,
117 double dh,
118 double rho,
119 double mu,
120 double *ustar2,
121 double *k,
122 double *eps);
123
124/*----------------------------------------------------------------------------*/
125/*
126 * \brief Calculation of \f$ k \f$ and \f$\varepsilon\f$
127 * from a diameter \f$ D_H \f$, a turbulent intensity \f$ I \f$
128 * and the reference velocity \f$ U_{ref} \f$
129 * for a circular duct flow with smooth wall
130 * (for inlet boundary conditions).
131 *
132 * \f[ k = 1.5 I {U_{ref}}^2 \f]
133 * \f[ \varepsilon = 10 \dfrac{{C_\mu}^{0.75} k^{1.5}}{ \kappa D_H} \f]
134 *
135 * \param[in] uref2 square of the reference flow velocity
136 * \param[in] t_intensity turbulent intensity \f$ I \f$
137 * \param[in] dh hydraulic diameter \f$ D_H \f$
138 * \param[out] k calculated turbulent intensity \f$ k \f$
139 * \param[out] eps calculated turbulent dissipation
140 * \f$ \varepsilon \f$
141 */
142/*----------------------------------------------------------------------------*/
143
144extern "C" void
145cs_turbulence_bc_ke_turb_intensity(double uref2,
146 double t_intensity,
147 double dh,
148 double *k,
149 double *eps);
150
151/*----------------------------------------------------------------------------*/
152/*
153 * \brief Set inlet boundary condition values for turbulence variables based
154 * on a diameter \f$ D_H \f$ and the reference velocity \f$ U_{ref} \f$
155 * for a circular duct flow with smooth wall
156 * (use for inlet boundary conditions).
157 *
158 * We use the laws from Idel'Cik, i.e.
159 * the head loss coefficient \f$ \lambda \f$ is defined by:
160 * \f[ |\dfrac{\Delta P}{\Delta x}| =
161 * \dfrac{\lambda}{D_H} \frac{1}{2} \rho U_{ref}^2 \f]
162 *
163 * then the relation reads \f$u^\star = U_{ref} \sqrt{\dfrac{\lambda}{8}}\f$.
164 * \f$\lambda \f$ depends on the hydraulic Reynolds number
165 * \f$ Re = \dfrac{U_{ref} D_H}{ \nu} \f$ and is given by:
166 * - for \f$ Re < 2000 \f$
167 * \f[ \lambda = \dfrac{64}{Re} \f]
168 *
169 * - for \f$ Re > 4000 \f$
170 * \f[ \lambda = \dfrac{1}{( 1.8 \log_{10}(Re)-1.64 )^2} \f]
171 *
172 * - for \f$ 2000 < Re < 4000 \f$, we complete by a straight line
173 * \f[ \lambda = 0.021377 + 5.3115. 10^{-6} Re \f]
174 *
175 * From \f$ u^\star \f$, we can estimate \f$ k \f$ and \f$ \varepsilon\f$
176 * from the well known formulae of developped turbulence
177 *
178 * \param[in] face_id boundary face id
179 * \param[in] uref2 square of the reference flow velocity
180 * \param[in] dh hydraulic diameter \f$ D_H \f$
181 * \param[in] rho mass density \f$ \rho \f$
182 * \param[in] mu dynamic viscosity \f$ \nu \f$
183 */
184/*----------------------------------------------------------------------------*/
185
186extern "C" void
187cs_turbulence_bc_inlet_hyd_diam(cs_lnum_t face_id,
188 double uref2,
189 double dh,
190 double rho,
191 double mu);
192
193/*----------------------------------------------------------------------------*/
194/*
195 * \brief Set inlet boundary condition values for turbulence variables based
196 * on a diameter \f$ D_H \f$, a turbulent intensity \f$ I \f$
197 * and the reference velocity \f$ U_{ref} \f$
198 * for a circular duct flow with smooth wall.
199 *
200 * \param[in] face_id boundary face id
201 * \param[in] uref2 square of the reference flow velocity
202 * \param[in] t_intensity turbulent intensity \f$ I \f$
203 * \param[in] dh hydraulic diameter \f$ D_H \f$
204 */
205/*----------------------------------------------------------------------------*/
206
207extern "C" void
208cs_turbulence_bc_inlet_turb_intensity(cs_lnum_t face_id,
209 double uref2,
210 double t_intensity,
211 double dh);
212
213/*----------------------------------------------------------------------------*/
214/*
215 * \brief Set inlet boundary condition values for turbulence variables based
216 * on given k and epsilon values.
217 *
218 * \param[in] face_id boundary face id
219 * \param[in] k turbulent kinetic energy
220 * \param[in] eps turbulent dissipation
221 */
222/*----------------------------------------------------------------------------*/
223
224extern "C" void
225cs_turbulence_bc_inlet_k_eps(cs_lnum_t face_id,
226 double k,
227 double eps);
228
229/*----------------------------------------------------------------------------*/
230/*
231 * \brief Set inlet boundary condition values for turbulence variables based
232 * on given k and epsilon or Rij values.
233 *
234 * \param[in] face_id boundary face id
235 * \param[in] k turbulent kinetic energy
236 * \param[in] eps turbulent dissipation
237 */
238/*----------------------------------------------------------------------------*/
239
240extern "C" void
241cs_turbulence_bc_set_uninit_inlet(cs_lnum_t face_num,
242 double k,
243 double rij[],
244 double eps);
245
246/*----------------------------------------------------------------------------*/
247/*
248 * \brief Set inlet boundary condition values for turbulence variables based
249 * on a diameter \f$ D_H \f$ and the reference velocity \f$ U_{ref} \f$
250 * for a circular duct flow with smooth wall, only if not already set.
251 *
252 * Apart from assigning values where not already initialized, this function
253 * is similar to \ref cs_turbulence_bc_inlet_hyd_diam.
254 *
255 * \param[in] face_id boundary face id
256 * \param[in] vel_dir velocity direction (not normalized, or NULL)
257 * \param[in] shear_dir shear direction (or NULL), it also contains the
258 * level of anisotropy (Rnt = a_nt k)
259 * \param[in] uref2 square of the reference flow velocity
260 * \param[in] dh hydraulic diameter \f$ D_H \f$
261 * \param[in] rho mass density \f$ \rho \f$
262 * \param[in] mu dynamic viscosity \f$ \nu \f$
263 */
264/*----------------------------------------------------------------------------*/
265
266extern "C" void
267cs_turbulence_bc_set_uninit_inlet_hyd_diam(cs_lnum_t face_id,
268 double vel_dir[],
269 double shear_dir[],
270 double uref2,
271 double dh,
272 double rho,
273 double mu);
274
275/*----------------------------------------------------------------------------*/
276/*
277 * \brief Set inlet boundary condition values for turbulence variables based
278 * on a diameter \f$ D_H \f$, a turbulent intensity \f$ I \f$
279 * and the reference velocity \f$ U_{ref} \f$
280 * for a circular duct flow with smooth wall, only if not already set.
281 *
282 * Apart from assigning values where not already initialized, this function
283 * is similar to \ref cs_turbulence_bc_inlet_turb_intensity.
284 *
285 * \param[in] face_id boundary face id
286 * \param[in] uref2 square of the reference flow velocity
287 * \param[in] t_intensity turbulent intensity \f$ I \f$
288 * \param[in] dh hydraulic diameter \f$ D_H \f$
289 */
290/*----------------------------------------------------------------------------*/
291
292extern "C" void
293cs_turbulence_bc_set_uninit_inlet_turb_intensity(cs_lnum_t face_id,
294 double uref2,
295 double t_intensity,
296 double dh);
297
298/*----------------------------------------------------------------------------*/
299/*
300 * \brief Set inlet boundary condition values for turbulence variables based
301 * on given k and epsilon values only if not already initialized.
302 *
303 * \param[in] face_id boundary face id
304 * \param[in] k turbulent kinetic energy
305 * \param[in] eps turbulent dissipation
306 * \param[in] vel_dir velocity direction
307 * \param[in] shear_dir shear direction
308 */
309/*----------------------------------------------------------------------------*/
310
311extern "C" void
312cs_turbulence_bc_set_uninit_inlet_k_eps(cs_lnum_t face_id,
313 double k,
314 double eps,
315 double vel_dir[],
316 double shear_dir[]);
317
318/*----------------------------------------------------------------------------*/
319/*
320 * \brief Assign homogeneous Neumann turbulent boundary conditions to
321 * a given face.
322 *
323 * This is useful for outgoing flow.
324 */
325/*----------------------------------------------------------------------------*/
326
327extern "C" void
328cs_turbulence_bc_set_hmg_neumann(cs_lnum_t face_id);
329
330/*----------------------------------------------------------------------------*/
331/*
332 * \brief Compute matrix \f$\tens{\alpha}\f$ used in the computation of the
333 * Reynolds stress tensor boundary conditions.
334 *
335 * We note \f$\tens{R}_g\f$ the Reynolds Stress tensor in the global reference
336 * frame (mesh reference frame) and \f$\tens{R}_l\f$ the Reynolds stress
337 * tensor in the local reference frame (reference frame associated to the
338 * boundary face).
339 *
340 * \f$\tens{P}_{lg}\f$ is the change of basis orthogonal matrix from local
341 * to global reference frame.
342
343 * \f$\tens{\alpha}\f$ is a 6 by 6 matrix such that:
344 * \f[
345 * \vect{R}_{g,\fib} = \tens{\alpha} \vect{R}_{g,\centip} + \vect{R}_{g}^*
346 * \f]
347 * where symetric tensors \f$\tens{R}_g\f$ have been unfolded as follows:
348 * \f[
349 * \vect{R}_g = \transpose{\left(R_{g,11},R_{g,22},R_{g,33},
350 * R_{g,12},R_{g,13},R_{g,23}\right)}
351 * \f].
352 *
353 * \f$\tens{\alpha}\f$ is defined so that \f$ \tens{R}_{g,\fib} \f$ is computed
354 * as a function of \f$\tens{R}_{g,\centip}\f$ as follows:
355 * \f[
356 * \tens{R}_{g,\fib}=\tens{P}_{lg}\tens{R}_{l,\fib}\transpose{\tens{P}_{lg}}
357 * \f]
358 *
359 * with
360 * \f[
361 * \tens{R}_{l,\fib} =
362 * \begin{bmatrix}
363 * R_{l,11,\centip} & u^* u_k & c R_{l,13,\centip}\\
364 * u^* u_k & R_{l,22,\centip} & 0 \\
365 * c R_{l,13,\centip} & 0 & R_{l,33,\centip}
366 * \end{bmatrix} +
367 * \underbrace{\begin{bmatrix}
368 * 0 & u^* u_k & 0 \\
369 * u^* u_k & 0 & 0 \\
370 * 0 & 0 & 0
371 * \end{bmatrix}}_{\tens{R}_l^*}
372 * \f]
373 *
374 * and
375 * \f$\tens{R}_{l,\centip}=\transpose{\tens{P}_{lg}}\tens{R}_{g,\centip}
376 * \tens{P}_{lg}\f$.
377 *
378 * Constant c is chosen depending on the type of the boundary face:
379 * \f$c = 0\f$ at a wall face, \f$c = 1\f$ at a symmetry face.
380 *
381 * \param[in] is_sym Constant c in description above
382 * (1 at a symmetry face, 0 at a wall face)
383 * \param[in] p_lg change of basis matrix (local to global)
384 * \param[out] alpha transformation matrix
385 */
386/*----------------------------------------------------------------------------*/
387
388#ifdef __cplusplus
389
390CS_F_HOST_DEVICE inline void
391cs_turbulence_bc_rij_transform(int is_sym,
392 cs_real_t p_lg[3][3],
393 cs_real_t alpha[6][6])
394{
395 cs_real_t p_lg2[3][3];
396 for (int ii = 0; ii < 3; ii++) {
397 for (int jj = 0; jj < 3; jj++)
398 p_lg2[ii][jj] = cs::pow2(p_lg[ii][jj]);
399 }
400
401 /* alpha(i,j) for i in [1,3] and j in [1,3]: 9 terms */
402 for (int ii = 0; ii < 3; ii++) {
403 for (int jj = 0; jj < 3; jj++) {
404 alpha[jj][ii] = p_lg2[0][ii] * p_lg2[0][jj]
405 + p_lg2[1][ii] * p_lg2[1][jj]
406 + p_lg2[2][ii] * p_lg2[2][jj]
407 + 2. * is_sym * p_lg[0][ii] * p_lg[2][ii]
408 * p_lg[0][jj] * p_lg[2][jj];
409 }
410 }
411
412 /* alpha(i,j) for i in [1,3] and j in [4,6]: 9 terms */
413
414 /* Correcponding line to j */
415 const int _jj_to_kk[3] = {0, 1, 0};
416 /* Corresponding column to j*/
417 const int _jj_to_pp[3] = {1, 2, 2};
418
419 for (int ii = 0; ii < 3; ii++) {
420 for (int jj = 0; jj < 3; jj++) {
421
422 int kk = _jj_to_kk[jj];
423 int pp = _jj_to_pp[jj];
424
425 alpha[jj + 3][ii] =
426 2. * ( p_lg2[0][ii] * p_lg[0][kk] * p_lg[0][pp]
427 + p_lg2[1][ii] * p_lg[1][kk] * p_lg[1][pp]
428 + p_lg2[2][ii] * p_lg[2][kk] * p_lg[2][pp]
429 + is_sym * p_lg[2][ii] * p_lg[0][ii]
430 * ( p_lg[0][kk]*p_lg[2][pp]
431 + p_lg[2][kk]*p_lg[0][pp]));
432 }
433 }
434
435 /* alpha(i,j) for i in [4,6] and j in [1,3]: 9 terms */
436
437 for (int ii = 0; ii < 3; ii++) {
438 for (int jj = 0; jj < 3; jj++) {
439
440 int kk = _jj_to_kk[ii];
441 int pp = _jj_to_pp[ii];
442
443 /* Note: could be simplified because it is 0.5*alpha[jj+3][ii] */
444 alpha[jj][ii + 3] =
445 p_lg[0][kk] * p_lg[0][pp] * p_lg2[0][jj]
446 + p_lg[1][kk] * p_lg[1][pp] * p_lg2[1][jj]
447 + p_lg[2][kk] * p_lg[2][pp] * p_lg2[2][jj]
448 + is_sym * p_lg[2][jj] * p_lg[0][jj]
449 * ( p_lg[0][kk]*p_lg[2][pp]
450 + p_lg[2][kk]*p_lg[0][pp]);
451 }
452 }
453
454 /* alpha(i,j) for i in [4,6] and j in [4,6]: 9 terms */
455 for (int ii = 0; ii < 3; ii++) {
456 for (int jj = 0; jj < 3; jj++) {
457
458 int kk = _jj_to_kk[ii];
459 int pp = _jj_to_pp[ii];
460
461 int jj1 = _jj_to_kk[jj];
462 int jj2 = _jj_to_pp[jj];
463
464 alpha[jj + 3][ii + 3] =
465 2. * ( p_lg[0][kk] * p_lg[0][pp] * p_lg[0][jj1] * p_lg[0][jj2]
466 + p_lg[1][kk] * p_lg[1][pp] * p_lg[1][jj1] * p_lg[1][jj2]
467 + p_lg[2][kk] * p_lg[2][pp] * p_lg[2][jj1] * p_lg[2][jj2])
468 + is_sym * ( p_lg[0][kk]*p_lg[2][pp]
469 + p_lg[2][kk]*p_lg[0][pp])
470 * ( p_lg[2][jj1]*p_lg[0][jj2]
471 + p_lg[0][jj1]*p_lg[2][jj2]);
472 }
473 }
474}
475
476#endif
477
478/*----------------------------------------------------------------------------*/
479
480#endif /* CS_TURBULENCE_BC_H */
cs_defs.h
CS_F_HOST_DEVICE
#define CS_F_HOST_DEVICE
Definition: cs_defs.h:555
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:332
cs_lnum_t
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:325
k
@ k
Definition: cs_field_pointer.h:68
mu
@ mu
Definition: cs_field_pointer.h:102
eps
@ eps
Definition: cs_field_pointer.h:69
rho
@ rho
Definition: cs_field_pointer.h:96
rij
@ rij
Definition: cs_field_pointer.h:71
cs_turbulence_bc_ke_turb_intensity
void cs_turbulence_bc_ke_turb_intensity(double uref2, double t_intensity, double dh, double *k, double *eps)
Calculation of and from a diameter , a turbulent intensity and the reference velocity for a circu...
Definition: cs_turbulence_bc.cpp:701
cs_turbulence_bc_inlet_hyd_diam
void cs_turbulence_bc_inlet_hyd_diam(cs_lnum_t face_id, double uref2, double dh, double rho, double mu)
Set inlet boundary condition values for turbulence variables based on a diameter and the reference v...
Definition: cs_turbulence_bc.cpp:745
cs_turbulence_bc_set_uninit_inlet_turb_intensity
void cs_turbulence_bc_set_uninit_inlet_turb_intensity(cs_lnum_t face_id, double uref2, double t_intensity, double dh)
Set inlet boundary condition values for turbulence variables based on a diameter ,...
Definition: cs_turbulence_bc.cpp:978
cs_turbulence_bc_inlet_k_eps
void cs_turbulence_bc_inlet_k_eps(cs_lnum_t face_id, double k, double eps)
Set inlet boundary condition values for turbulence variables based on given k and epsilon values.
Definition: cs_turbulence_bc.cpp:797
cs_turbulence_bc_free_pointers
void cs_turbulence_bc_free_pointers(void)
Free memory allocations for turbulence boundary condition pointers.
Definition: cs_turbulence_bc.cpp:615
cs_turbulence_bc_set_uninit_inlet
void cs_turbulence_bc_set_uninit_inlet(cs_lnum_t face_num, double k, double rij[], double eps)
Definition: cs_turbulence_bc.cpp:952
cs_turbulence_bc_set_uninit_inlet_k_eps
void cs_turbulence_bc_set_uninit_inlet_k_eps(cs_lnum_t face_id, double k, double eps, double vel_dir[], double shear_dir[])
Set inlet boundary condition values for turbulence variables based on given k and epsilon values only...
Definition: cs_turbulence_bc.cpp:1005
cs_turbulence_bc_ke_hyd_diam
void cs_turbulence_bc_ke_hyd_diam(double uref2, double dh, double rho, double mu, double *ustar2, double *k, double *eps)
Calculation of , and from a diameter and the reference velocity for a circular duct flow with smo...
Definition: cs_turbulence_bc.cpp:669
cs_turbulence_bc_init_pointers
void cs_turbulence_bc_init_pointers(void)
Initialize turbulence model boundary condition pointers.
Definition: cs_turbulence_bc.cpp:524
cs_turbulence_bc_inlet_turb_intensity
void cs_turbulence_bc_inlet_turb_intensity(cs_lnum_t face_id, double uref2, double t_intensity, double dh)
Set inlet boundary condition values for turbulence variables based on a diameter ,...
Definition: cs_turbulence_bc.cpp:773
cs_turbulence_bc_rij_transform
CS_F_HOST_DEVICE void cs_turbulence_bc_rij_transform(int is_sym, cs_real_t p_lg[3][3], cs_real_t alpha[6][6])
Definition: cs_turbulence_bc.h:391
cs_turbulence_bc_set_hmg_neumann
void cs_turbulence_bc_set_hmg_neumann(cs_lnum_t face_id)
Assign homogeneous Neumann turbulent boundary conditions to a given face.
Definition: cs_turbulence_bc.cpp:814
cs_turbulence_bc_set_uninit_inlet_hyd_diam
void cs_turbulence_bc_set_uninit_inlet_hyd_diam(cs_lnum_t face_id, double vel_dir[], double shear_dir[], double uref2, double dh, double rho, double mu)
Set inlet boundary condition values for turbulence variables based on a diameter and the reference v...
Definition: cs_turbulence_bc.cpp:925
cs::pow2
CS_F_HOST_DEVICE float pow2(float x)
Compute the square of a real value.
Definition: cs_math.h:744