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-2021 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 "cs_defs.h" 35 36 /*----------------------------------------------------------------------------*/ 37 38 BEGIN_C_DECLS 39 40 /*============================================================================= 41 * Macro definitions 42 *============================================================================*/ 43 44 /*============================================================================ 45 * Type definitions 46 *============================================================================*/ 47 48 /*============================================================================ 49 * Global variables 50 *============================================================================*/ 51 52 /*============================================================================= 53 * Public function prototypes 54 *============================================================================*/ 55 56 /*----------------------------------------------------------------------------*/ 57 /*! 58 * Initialize turbulence model boundary condition ids. 59 */ 60 /*----------------------------------------------------------------------------*/ 61 62 void 63 cs_turbulence_model_init_bc_ids(void); 64 65 /*----------------------------------------------------------------------------*/ 66 /*! 67 * \brief Free memory allocations for turbulence boundary conditions ids. 68 */ 69 /*----------------------------------------------------------------------------*/ 70 71 void 72 cs_turbulence_model_free_bc_ids(void); 73 74 /*----------------------------------------------------------------------------*/ 75 /*! 76 * \brief Calculation of \f$ u^\star \f$, \f$ k \f$ and \f$\varepsilon \f$ 77 * from a diameter \f$ D_H \f$ and the reference velocity \f$ U_{ref} \f$ 78 * for a circular duct flow with smooth wall 79 * (use for inlet boundary conditions). 80 * 81 * Both \f$ u^\star \f$ and\f$ (k,\varepsilon )\f$ are returned, so that 82 * the user may compute other values of \f$ k \f$ and \f$ \varepsilon \f$ 83 * with \f$ u^\star \f$. 84 * 85 * We use the laws from Idel'Cik, i.e. 86 * the head loss coefficient \f$ \lambda \f$ is defined by: 87 * \f[ |\dfrac{\Delta P}{\Delta x}| = 88 * \dfrac{\lambda}{D_H} \frac{1}{2} \rho U_{ref}^2 \f] 89 * 90 * then the relation reads \f$u^\star = U_{ref} \sqrt{\dfrac{\lambda}{8}}\f$. 91 * \f$\lambda \f$ depends on the hydraulic Reynolds number 92 * \f$ Re = \dfrac{U_{ref} D_H}{ \nu} \f$ and is given by: 93 * - for \f$ Re < 2000 \f$ 94 * \f[ \lambda = \dfrac{64}{Re} \f] 95 * 96 * - for \f$ Re > 4000 \f$ 97 * \f[ \lambda = \dfrac{1}{( 1.8 \log_{10}(Re)-1.64 )^2} \f] 98 * 99 * - for \f$ 2000 < Re < 4000 \f$, we complete by a straight line 100 * \f[ \lambda = 0.021377 + 5.3115. 10^{-6} Re \f] 101 * 102 * From \f$ u^\star \f$, we can estimate \f$ k \f$ and \f$ \varepsilon\f$ 103 * from the well known formulae of developped turbulence 104 * 105 * \f[ k = \dfrac{u^{\star 2}}{\sqrt{C_\mu}} \f] 106 * \f[ \varepsilon = \dfrac{ u^{\star 3}}{(\kappa D_H /10)} \f] 107 * 108 * \param[in] uref2 square of the reference flow velocity 109 * \param[in] dh hydraulic diameter \f$ D_H \f$ 110 * \param[in] rho mass density \f$ \rho \f$ 111 * \param[in] mu dynamic viscosity \f$ \nu \f$ 112 * \param[out] ustar2 square of friction speed 113 * \param[out] k calculated turbulent intensity \f$ k \f$ 114 * \param[out] eps calculated turbulent dissipation 115 * \f$ \varepsilon \f$ 116 */ 117 /*----------------------------------------------------------------------------*/ 118 119 void 120 cs_turbulence_bc_ke_hyd_diam(double uref2, 121 double dh, 122 double rho, 123 double mu, 124 double *ustar2, 125 double *k, 126 double *eps); 127 128 /*----------------------------------------------------------------------------*/ 129 /*! 130 * \brief Calculation of \f$ k \f$ and \f$\varepsilon\f$ 131 * from a diameter \f$ D_H \f$, a turbulent intensity \f$ I \f$ 132 * and the reference velocity \f$ U_{ref} \f$ 133 * for a circular duct flow with smooth wall 134 * (for inlet boundary conditions). 135 * 136 * \f[ k = 1.5 I {U_{ref}}^2 \f] 137 * \f[ \varepsilon = 10 \dfrac{{C_\mu}^{0.75} k^{1.5}}{ \kappa D_H} \f] 138 * 139 * \param[in] uref2 square of the reference flow velocity 140 * \param[in] t_intensity turbulent intensity \f$ I \f$ 141 * \param[in] dh hydraulic diameter \f$ D_H \f$ 142 * \param[out] k calculated turbulent intensity \f$ k \f$ 143 * \param[out] eps calculated turbulent dissipation 144 * \f$ \varepsilon \f$ 145 */ 146 /*----------------------------------------------------------------------------*/ 147 148 void 149 cs_turbulence_bc_ke_turb_intensity(double uref2, 150 double t_intensity, 151 double dh, 152 double *k, 153 double *eps); 154 155 /*----------------------------------------------------------------------------*/ 156 /*! 157 * \brief Set inlet boundary condition values for turbulence variables based 158 * on a diameter \f$ D_H \f$ and the reference velocity \f$ U_{ref} \f$ 159 * for a circular duct flow with smooth wall 160 * (use for inlet boundary conditions). 161 * 162 * We use the laws from Idel'Cik, i.e. 163 * the head loss coefficient \f$ \lambda \f$ is defined by: 164 * \f[ |\dfrac{\Delta P}{\Delta x}| = 165 * \dfrac{\lambda}{D_H} \frac{1}{2} \rho U_{ref}^2 \f] 166 * 167 * then the relation reads \f$u^\star = U_{ref} \sqrt{\dfrac{\lambda}{8}}\f$. 168 * \f$\lambda \f$ depends on the hydraulic Reynolds number 169 * \f$ Re = \dfrac{U_{ref} D_H}{ \nu} \f$ and is given by: 170 * - for \f$ Re < 2000 \f$ 171 * \f[ \lambda = \dfrac{64}{Re} \f] 172 * 173 * - for \f$ Re > 4000 \f$ 174 * \f[ \lambda = \dfrac{1}{( 1.8 \log_{10}(Re)-1.64 )^2} \f] 175 * 176 * - for \f$ 2000 < Re < 4000 \f$, we complete by a straight line 177 * \f[ \lambda = 0.021377 + 5.3115. 10^{-6} Re \f] 178 * 179 * From \f$ u^\star \f$, we can estimate \f$ k \f$ and \f$ \varepsilon\f$ 180 * from the well known formulae of developped turbulence 181 * 182 * \param[in] face_id boundary face id 183 * \param[in] uref2 square of the reference flow velocity 184 * \param[in] dh hydraulic diameter \f$ D_H \f$ 185 * \param[in] rho mass density \f$ \rho \f$ 186 * \param[in] mu dynamic viscosity \f$ \nu \f$ 187 * \param[out] rcodcl boundary condition values 188 */ 189 /*----------------------------------------------------------------------------*/ 190 191 void 192 cs_turbulence_bc_inlet_hyd_diam(cs_lnum_t face_id, 193 double uref2, 194 double dh, 195 double rho, 196 double mu, 197 double *rcodcl); 198 199 /*----------------------------------------------------------------------------*/ 200 /*! 201 * \brief Set inlet boundary condition values for turbulence variables based 202 * on a diameter \f$ D_H \f$, a turbulent intensity \f$ I \f$ 203 * and the reference velocity \f$ U_{ref} \f$ 204 * for a circular duct flow with smooth wall. 205 * 206 * \param[in] face_id boundary face id 207 * \param[in] uref2 square of the reference flow velocity 208 * \param[in] t_intensity turbulent intensity \f$ I \f$ 209 * \param[in] dh hydraulic diameter \f$ D_H \f$ 210 * \param[out] rcodcl boundary condition values 211 */ 212 /*----------------------------------------------------------------------------*/ 213 214 void 215 cs_turbulence_bc_inlet_turb_intensity(cs_lnum_t face_id, 216 double uref2, 217 double t_intensity, 218 double dh, 219 double *rcodcl); 220 221 /*----------------------------------------------------------------------------*/ 222 /*! 223 * \brief Set inlet boundary condition values for turbulence variables based 224 * on given k and epsilon values. 225 * 226 * \param[in] face_id boundary face id 227 * \param[in] k turbulent kinetic energy 228 * \param[in] eps turbulent dissipation 229 * \param[out] rcodcl boundary condition values 230 */ 231 /*----------------------------------------------------------------------------*/ 232 233 void 234 cs_turbulence_bc_inlet_k_eps(cs_lnum_t face_id, 235 double k, 236 double eps, 237 double *rcodcl); 238 239 /*----------------------------------------------------------------------------*/ 240 /*! 241 * \brief Set inlet boundary condition values for turbulence variables based 242 * on given k and epsilon values only if not already initialized. 243 * 244 * \param[in] face_id boundary face id 245 * \param[in] k turbulent kinetic energy 246 * \param[in] eps turbulent dissipation 247 * \param[in] vel_dir velocity direction 248 * \param[in] shear_dir shear direction 249 * \param[out] rcodcl boundary condition values 250 */ 251 /*----------------------------------------------------------------------------*/ 252 253 void 254 cs_turbulence_bc_set_uninit_inlet_k_eps(cs_lnum_t face_id, 255 double k, 256 double eps, 257 double vel_dir[], 258 double shear_dir[], 259 double *rcodcl); 260 261 /*----------------------------------------------------------------------------*/ 262 /*! 263 * \brief Compute matrix \f$\tens{\alpha}\f$ used in the computation of the 264 * Reynolds stress tensor boundary conditions. 265 * 266 * We note \f$\tens{R}_g\f$ the Reynolds Stress tensor in the global reference 267 * frame (mesh reference frame) and \f$\tens{R}_l\f$ the Reynolds stress 268 * tensor in the local reference frame (reference frame associated to the 269 * boundary face). 270 * 271 * \f$\tens{P}_{lg}\f$ is the change of basis orthogonal matrix from local 272 * to global reference frame. 273 274 * \f$\tens{\alpha}\f$ is a 6 by 6 matrix such that: 275 * \f[ 276 * \vect{R}_{g,\fib} = \tens{\alpha} \vect{R}_{g,\centip} + \vect{R}_{g}^* 277 * \f] 278 * where symetric tensors \f$\tens{R}_g\f$ have been unfolded as follows: 279 * \f[ 280 * \vect{R}_g = \transpose{\left(R_{g,11},R_{g,22},R_{g,33}, 281 * R_{g,12},R_{g,13},R_{g,23}\right)} 282 * \f]. 283 * 284 * \f$\tens{\alpha}\f$ is defined so that \f$ \tens{R}_{g,\fib} \f$ is computed 285 * as a function of \f$\tens{R}_{g,\centip}\f$ as follows: 286 * \f[ 287 * \tens{R}_{g,\fib}=\tens{P}_{lg}\tens{R}_{l,\fib}\transpose{\tens{P}_{lg}} 288 * \f] 289 * 290 * with 291 * \f[ 292 * \tens{R}_{l,\fib} = 293 * \begin{bmatrix} 294 * R_{l,11,\centip} & u^* u_k & c R_{l,13,\centip}\\ 295 * u^* u_k & R_{l,22,\centip} & 0 \\ 296 * c R_{l,13,\centip} & 0 & R_{l,33,\centip} 297 * \end{bmatrix} + 298 * \underbrace{\begin{bmatrix} 299 * 0 & u^* u_k & 0 \\ 300 * u^* u_k & 0 & 0 \\ 301 * 0 & 0 & 0 302 * \end{bmatrix}}_{\tens{R}_l^*} 303 * \f] 304 * 305 * and 306 * \f$\tens{R}_{l,\centip}=\transpose{\tens{P}_{lg}}\tens{R}_{g,\centip} 307 * \tens{P}_{lg}\f$. 308 * 309 * Constant c is chosen depending on the type of the boundary face: 310 * \f$c = 0\f$ at a wall face, \f$c = 1\f$ at a symmetry face. 311 * 312 * \param[in] is_sym Constant c in description above 313 * (1 at a symmetry face, 0 at a wall face) 314 * \param[in] p_lg change of basis matrix (local to global) 315 * \param[out] alpha transformation matrix 316 */ 317 /*----------------------------------------------------------------------------*/ 318 319 void 320 cs_turbulence_bc_rij_transform(int is_sym, 321 cs_real_t p_lg[3][3], 322 cs_real_t alpha[6][6]); 323 324 /*----------------------------------------------------------------------------*/ 325 326 END_C_DECLS 327 328 #endif /* __CS_TURBULENCE_BC_H__ */ 329