1 /*============================================================================
2 * Wall functions
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_config.h"
28 #include "cs_defs.h"
29
30 /*----------------------------------------------------------------------------
31 * Standard C library headers
32 *----------------------------------------------------------------------------*/
33
34 #include <assert.h>
35 #include <errno.h>
36 #include <stdio.h>
37 #include <stdarg.h>
38 #include <string.h>
39 #include <math.h>
40 #include <stdlib.h>
41 #include <float.h>
42
43 #if defined(HAVE_MPI)
44 #include <mpi.h>
45 #endif
46
47 /*----------------------------------------------------------------------------
48 * BFT library headers
49 *----------------------------------------------------------------------------*/
50
51 #include <bft_error.h>
52 #include <bft_mem.h>
53 #include <bft_printf.h>
54
55 /*----------------------------------------------------------------------------
56 * Local headers
57 *----------------------------------------------------------------------------*/
58
59 #include "cs_log.h"
60 #include "cs_mesh.h"
61 #include "cs_mesh_quantities.h"
62 #include "cs_turbulence_model.h"
63
64 /*----------------------------------------------------------------------------
65 * Header for the current file
66 *----------------------------------------------------------------------------*/
67
68 #include "cs_wall_functions.h"
69
70 /*----------------------------------------------------------------------------*/
71
72 BEGIN_C_DECLS
73
74 /*=============================================================================
75 * Additional doxygen documentation
76 *============================================================================*/
77
78 /*!
79 \file cs_wall_functions.c
80 Wall functions descriptor and computation.
81 */
82 /*----------------------------------------------------------------------------*/
83
84 /*! \struct cs_wall_functions_t
85
86 \brief wall functions descriptor.
87
88 Members of this wall functions descriptor are publicly accessible, to allow
89 for concise syntax, as it is expected to be used in many places.
90
91 \var cs_wall_functions_t::iwallf
92 Indicates the type of wall function used for the velocity
93 boundary conditions on a frictional wall.\n
94 - CS_WALL_F_DISABLED: no wall functions
95 - CS_WALL_F_1SCALE_POWER: one scale of friction velocities (power law)
96 - CS_WALL_F_1SCALE_LOG: one scale of friction velocities (log law)
97 - CS_WALL_F_2SCALES_LOG: two scales of friction velocities (log law)
98 - CS_WALL_F_SCALABLE_2SCALES_LOG: two scales of friction velocities
99 (log law) (scalable wall functions)
100 - CS_WALL_F_2SCALES_VDRIEST: two scales of friction velocities
101 (mixing length based on V. Driest analysis)
102 - CS_WALL_F_2SCALES_SMOOTH_ROUGH: wall function unifying rough and smooth
103 friction regimes
104 - CS_WALL_F_2SCALES_CONTINUOUS: All \f$ y^+ \f$ for low Reynolds models\n
105 \ref iwallf is initialised to CS_WALL_F_1SCALE_LOG for \ref iturb = 10,
106 40, 41 or 70
107 (mixing length, LES and Spalart Allmaras).\n
108 \ref iwallf is initialised to CS_WALL_F_DISABLED for \ref iturb = 0, 32,
109 50 or 51\n
110 \ref iwallf is initialised to CS_WALL_F_2SCALES_LOG for \ref iturb = 20,
111 21, 30, 31 or 60
112 (\f$k-\epsilon\f$, \f$R_{ij}-\epsilon\f$ LRR, \f$R_{ij}-\epsilon\f$ SSG
113 and \f$k-\omega\f$ SST models).\n
114 The v2f model (\ref iturb=50) is not designed to use wall functions
115 (the mesh must be low Reynolds).\n
116 The value \ref iwallf = CS_WALL_F_2SCALES_LOG is not compatible with
117 \ref iturb=0, 10, 40
118 or 41 (laminar, mixing length and LES).\n
119 Concerning the \f$k-\epsilon\f$ and \f$R_{ij}-\epsilon\f$ models, the
120 two-scales model is usually at least as satisfactory as the one-scale
121 model.\n
122 The scalable wall function allows to virtually shift the wall when
123 necessary in order to be always in a logarithmic layer. It is used to make up for
124 the problems related to the use of High-Reynolds models on very refined
125 meshes.\n
126 Useful if \ref iturb is different from 50.
127 \var cs_wall_functions_t::iwalfs
128 wall functions for scalar
129 - CS_WALL_F_S_ARPACI_LARSEN: three layers (Arpaci and Larsen) or two layers (Prandtl-Taylor) for
130 Prandtl number smaller than 0.1
131 - CS_WALL_F_S_VDRIEST: consistant with the 2 scales wall function for velocity using Van
132 Driest mixing length
133 - CS_WALL_F_S_LOUIS: default wall function for atmospheric flows for
134 potential temperature. This has influence on the dynamic.
135 - CS_WALL_F_S_MONIN_OBUKHOV: Monin Obukhov wall function for atmospheric flows for
136 potential temperature. This has influence on the dynamic.
137
138 \var cs_wall_functions_t::ypluli
139 limit value of \f$y^+\f$ for the viscous sublayer
140
141 \ref ypluli depends on the chosen wall function: it is initialized to
142 10.88 for the scalable wall function (\ref iwallf=CS_WALL_F_SCALABLE_2SCALES_LOG), otherwise it is
143 initialized to \f$1/\kappa\approx 2,38\f$. In LES, \ref ypluli is taken
144 by default to be 10.88. Always useful.
145
146 */
147 /*----------------------------------------------------------------------------*/
148
149 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
150
151 /*============================================================================
152 * Static global variables
153 *============================================================================*/
154
155 /* wall functions structure and associated pointer */
156
157 static cs_wall_functions_t _wall_functions =
158 {
159 .iwallf = -999,
160 .iwalfs = -999,
161 .ypluli = -1e13
162 };
163
164 const cs_wall_functions_t * cs_glob_wall_functions = &_wall_functions;
165
166 /*============================================================================
167 * Prototypes for functions intended for use only by Fortran wrappers.
168 * (descriptions follow, with function bodies).
169 *============================================================================*/
170
171 void
172 cs_f_wall_functions_get_pointers(int **iwallf,
173 int **iwalfs,
174 double **ypluli);
175
176 /*============================================================================
177 * Fortran wrapper function definitions
178 *============================================================================*/
179
180 /*----------------------------------------------------------------------------
181 * Get pointers to members of the wall functions structure.
182 *
183 * This function is intended for use by Fortran wrappers, and
184 * enables mapping to Fortran global pointers.
185 *
186 * parameters:
187 * iwallf --> pointer to cs_glob_wall_functions->iwallf
188 * iwalfs --> pointer to cs_glob_wall_functions->iwalfs
189 * ypluli --> pointer to cs_glob_wall_functions->ypluli
190 *----------------------------------------------------------------------------*/
191
192 void
cs_f_wall_functions_get_pointers(int ** iwallf,int ** iwalfs,double ** ypluli)193 cs_f_wall_functions_get_pointers(int **iwallf,
194 int **iwalfs,
195 double **ypluli)
196 {
197 *iwallf = (int *)&(_wall_functions.iwallf);
198 *iwalfs = (int *)&(_wall_functions.iwalfs);
199 *ypluli = &(_wall_functions.ypluli);
200 }
201
202 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
203
204 /*============================================================================
205 * Public function definitions for Fortran API
206 *============================================================================*/
207
208 /*----------------------------------------------------------------------------
209 * Wrapper to cs_wall_functions_velocity
210 *----------------------------------------------------------------------------*/
211
CS_PROCF(wallfunctions,WALLFUNCTIONS)212 void CS_PROCF (wallfunctions, WALLFUNCTIONS)
213 (
214 const int *const iwallf,
215 const cs_lnum_t *const ifac,
216 const cs_real_t *const l_visc,
217 const cs_real_t *const t_visc,
218 const cs_real_t *const vel,
219 const cs_real_t *const y,
220 const cs_real_t *const rough_d,
221 const cs_real_t *const rnnb,
222 const cs_real_t *const kinetic_en,
223 int *iuntur,
224 cs_lnum_t *nsubla,
225 cs_lnum_t *nlogla,
226 cs_real_t *ustar,
227 cs_real_t *uk,
228 cs_real_t *yplus,
229 cs_real_t *ypup,
230 cs_real_t *cofimp,
231 cs_real_t *dplus
232 )
233 {
234 assert(*iwallf >= 0 && *iwallf <= 7);
235
236 cs_wall_functions_velocity((cs_wall_f_type_t)*iwallf,
237 *ifac,
238 *l_visc,
239 *t_visc,
240 *vel,
241 *y,
242 *rough_d,
243 *rnnb,
244 *kinetic_en,
245 iuntur,
246 nsubla,
247 nlogla,
248 ustar,
249 uk,
250 yplus,
251 ypup,
252 cofimp,
253 dplus);
254 }
255
256 /*----------------------------------------------------------------------------
257 * Wrapper to cs_wall_functions_scalar
258 *----------------------------------------------------------------------------*/
259
CS_PROCF(hturbp,HTURBP)260 void CS_PROCF (hturbp, HTURBP)
261 (
262 const int *const iwalfs,
263 const cs_real_t *const l_visc,
264 const cs_real_t *const prl,
265 const cs_real_t *const prt,
266 const cs_real_t *const rough_t,
267 const cs_real_t *const uk,
268 const cs_real_t *const yplus,
269 const cs_real_t *const dplus,
270 cs_real_t *htur,
271 cs_real_t *yplim
272 )
273 {
274 cs_wall_functions_scalar((cs_wall_f_s_type_t)*iwalfs,
275 *l_visc,
276 *prl,
277 *prt,
278 *rough_t,
279 *uk,
280 *yplus,
281 *dplus,
282 htur,
283 yplim);
284 }
285
286 /*=============================================================================
287 * Public function definitions
288 *============================================================================*/
289
290 /*----------------------------------------------------------------------------
291 *! \brief Provide access to cs_glob_wall_functions
292 *
293 * needed to initialize structure with GUI
294 *----------------------------------------------------------------------------*/
295
296 cs_wall_functions_t *
cs_get_glob_wall_functions(void)297 cs_get_glob_wall_functions(void)
298 {
299 return &_wall_functions;
300 }
301
302 /*----------------------------------------------------------------------------*/
303 /*!
304 * \brief Compute the friction velocity and \f$y^+\f$ / \f$u^+\f$.
305 *
306 * \param[in] iwallf wall function type
307 * \param[in] ifac face number
308 * \param[in] l_visc kinematic viscosity
309 * \param[in] t_visc turbulent kinematic viscosity
310 * \param[in] vel wall projected cell center velocity
311 * \param[in] y wall distance
312 * \param[in] rough_d roughness length scale
313 * (not sand grain roughness)
314 * \param[in] rnnb \f$\vec{n}.(\tens{R}\vec{n})\f$
315 * \param[in] kinetic_en turbulent kinetic energy (cell center)
316 * \param[in,out] iuntur indicator: 0 in the viscous sublayer
317 * \param[in,out] nsubla counter of cell in the viscous sublayer
318 * \param[in,out] nlogla counter of cell in the log-layer
319 * \param[out] ustar friction velocity
320 * \param[out] uk friction velocity
321 * \param[out] yplus dimensionless distance to the wall
322 * \param[out] ypup yplus projected vel ratio
323 * \param[out] cofimp \f$\frac{|U_F|}{|U_I^p|}\f$ to ensure a good
324 * turbulence production
325 * \param[out] dplus dimensionless shift to the wall for scalable
326 * wall functions
327 */
328 /*----------------------------------------------------------------------------*/
329
330 void
cs_wall_functions_velocity(cs_wall_f_type_t iwallf,cs_lnum_t ifac,cs_real_t l_visc,cs_real_t t_visc,cs_real_t vel,cs_real_t y,cs_real_t rough_d,cs_real_t rnnb,cs_real_t kinetic_en,int * iuntur,cs_lnum_t * nsubla,cs_lnum_t * nlogla,cs_real_t * ustar,cs_real_t * uk,cs_real_t * yplus,cs_real_t * ypup,cs_real_t * cofimp,cs_real_t * dplus)331 cs_wall_functions_velocity(cs_wall_f_type_t iwallf,
332 cs_lnum_t ifac,
333 cs_real_t l_visc,
334 cs_real_t t_visc,
335 cs_real_t vel,
336 cs_real_t y,
337 cs_real_t rough_d,
338 cs_real_t rnnb,
339 cs_real_t kinetic_en,
340 int *iuntur,
341 cs_lnum_t *nsubla,
342 cs_lnum_t *nlogla,
343 cs_real_t *ustar,
344 cs_real_t *uk,
345 cs_real_t *yplus,
346 cs_real_t *ypup,
347 cs_real_t *cofimp,
348 cs_real_t *dplus)
349 {
350 cs_real_t lmk;
351 bool wf = true;
352
353 /* Pseudo shift of the wall, 0 by default */
354 *dplus = 0.;
355
356 /* Activation of wall function by default */
357 *iuntur = 1;
358
359 /* Get the adjacent border cell and its fluid/solid tag */
360 cs_lnum_t cell_id = cs_glob_mesh->b_face_cells[ifac-1] ;
361
362 /* If the cell is a solid cell, disable wall functions */
363
364 const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
365 if (mq->has_disable_flag) {
366 if (mq->c_disable_flag[cell_id])
367 iwallf = CS_WALL_F_DISABLED;
368 }
369
370 /* Sand Grain roughness */
371 cs_real_t sg_rough = rough_d * exp(cs_turb_xkappa*cs_turb_cstlog_rough);
372
373 switch (iwallf) {
374 case CS_WALL_F_DISABLED:
375 cs_wall_functions_disabled(l_visc,
376 t_visc,
377 vel,
378 y,
379 iuntur,
380 nsubla,
381 nlogla,
382 ustar,
383 uk,
384 yplus,
385 dplus,
386 ypup,
387 cofimp);
388 break;
389 case CS_WALL_F_1SCALE_POWER:
390 cs_wall_functions_1scale_power(l_visc,
391 vel,
392 y,
393 iuntur,
394 nsubla,
395 nlogla,
396 ustar,
397 uk,
398 yplus,
399 ypup,
400 cofimp);
401 break;
402 case CS_WALL_F_1SCALE_LOG:
403 cs_wall_functions_1scale_log(ifac,
404 l_visc,
405 vel,
406 y,
407 iuntur,
408 nsubla,
409 nlogla,
410 ustar,
411 uk,
412 yplus,
413 ypup,
414 cofimp);
415 break;
416 case CS_WALL_F_2SCALES_LOG:
417 cs_wall_functions_2scales_log(l_visc,
418 t_visc,
419 vel,
420 y,
421 kinetic_en,
422 iuntur,
423 nsubla,
424 nlogla,
425 ustar,
426 uk,
427 yplus,
428 ypup,
429 cofimp);
430 break;
431 case CS_WALL_F_SCALABLE_2SCALES_LOG:
432 cs_wall_functions_2scales_scalable(l_visc,
433 t_visc,
434 vel,
435 y,
436 kinetic_en,
437 iuntur,
438 nsubla,
439 nlogla,
440 ustar,
441 uk,
442 yplus,
443 dplus,
444 ypup,
445 cofimp);
446 break;
447 case CS_WALL_F_2SCALES_VDRIEST:
448 cs_wall_functions_2scales_vdriest(rnnb,
449 l_visc,
450 vel,
451 y,
452 kinetic_en,
453 iuntur,
454 nsubla,
455 nlogla,
456 ustar,
457 uk,
458 yplus,
459 ypup,
460 cofimp,
461 &lmk,
462 sg_rough,
463 wf);
464 break;
465 case CS_WALL_F_2SCALES_SMOOTH_ROUGH:
466 cs_wall_functions_2scales_smooth_rough(l_visc,
467 t_visc,
468 vel,
469 y,
470 rough_d,
471 kinetic_en,
472 iuntur,
473 nsubla,
474 nlogla,
475 ustar,
476 uk,
477 yplus,
478 dplus,
479 ypup,
480 cofimp);
481 break;
482 case CS_WALL_F_2SCALES_CONTINUOUS:
483 cs_wall_functions_2scales_continuous(rnnb,
484 l_visc,
485 t_visc,
486 vel,
487 y,
488 kinetic_en,
489 iuntur,
490 nsubla,
491 nlogla,
492 ustar,
493 uk,
494 yplus,
495 ypup,
496 cofimp);
497 break;
498 default:
499 break;
500 }
501 }
502
503 /*----------------------------------------------------------------------------*/
504 /*!
505 * \brief Compute the correction of the exchange coefficient between the
506 * fluid and the wall for a turbulent flow.
507 *
508 * This is function of the dimensionless
509 * distance to the wall \f$ y^+ = \dfrac{\centip \centf u_\star}{\nu}\f$.
510 *
511 * Then the return coefficient reads:
512 * \f[
513 * h_{tur} = Pr \dfrac{y^+}{T^+}
514 * \f]
515 *
516 * \param[in] iwalfs type of wall functions for scalar
517 * \param[in] l_visc kinematic viscosity
518 * \param[in] prl laminar Prandtl number
519 * \param[in] prt turbulent Prandtl number
520 * \param[in] rough_t scalar roughness lenghth scale
521 * \param[in] uk velocity scale based on TKE
522 * \param[in] yplus dimensionless distance to the wall
523 * \param[in] dplus dimensionless distance for scalable
524 * wall functions
525 * \param[out] htur corrected exchange coefficient
526 * \param[out] yplim value of the limit for \f$ y^+ \f$
527 */
528 /*----------------------------------------------------------------------------*/
529
530 void
cs_wall_functions_scalar(cs_wall_f_s_type_t iwalfs,cs_real_t l_visc,cs_real_t prl,cs_real_t prt,cs_real_t rough_t,cs_real_t uk,cs_real_t yplus,cs_real_t dplus,cs_real_t * htur,cs_real_t * yplim)531 cs_wall_functions_scalar(cs_wall_f_s_type_t iwalfs,
532 cs_real_t l_visc,
533 cs_real_t prl,
534 cs_real_t prt,
535 cs_real_t rough_t,
536 cs_real_t uk,
537 cs_real_t yplus,
538 cs_real_t dplus,
539 cs_real_t *htur,
540 cs_real_t *yplim)
541 {
542 switch (iwalfs) {
543 case CS_WALL_F_S_ARPACI_LARSEN:
544 cs_wall_functions_s_arpaci_larsen(l_visc,
545 prl,
546 prt,
547 rough_t,
548 uk,
549 yplus,
550 dplus,
551 htur,
552 yplim);
553 break;
554 case CS_WALL_F_S_VDRIEST:
555 cs_wall_functions_s_vdriest(prl,
556 prt,
557 yplus,
558 htur);
559 break;
560 case CS_WALL_F_S_SMOOTH_ROUGH:
561 cs_wall_functions_s_smooth_rough(l_visc,
562 prl,
563 prt,
564 rough_t,
565 uk,
566 yplus,
567 dplus,
568 htur);
569 break;
570 default:
571 /* TODO Monin Obukhov or Louis atmospheric wall function
572 * must be adapted to smooth wall functions.
573 * Arpaci and Larsen wall functions are put as in previous versions of
574 * Code_Saturne.
575 * */
576 cs_wall_functions_s_arpaci_larsen(l_visc,
577 prl,
578 prt,
579 rough_t,
580 uk,
581 yplus,
582 dplus,
583 htur,
584 yplim);
585 break;
586 }
587 }
588
589 /*----------------------------------------------------------------------------*/
590
591 END_C_DECLS
592