1 #ifndef __CS_WALL_FUNCTIONS_H__
2 #define __CS_WALL_FUNCTIONS_H__
3
4 /*============================================================================
5 * Wall functions
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 * BFT library headers
32 *----------------------------------------------------------------------------*/
33
34 #include <bft_printf.h>
35
36 /*----------------------------------------------------------------------------
37 * Local headers
38 *----------------------------------------------------------------------------*/
39
40 #include "cs_base.h"
41 #include "cs_math.h"
42 #include "cs_turbulence_model.h"
43
44 /*----------------------------------------------------------------------------*/
45
46 BEGIN_C_DECLS
47
48 /*=============================================================================
49 * Local Macro definitions
50 *============================================================================*/
51
52 /*============================================================================
53 * Type definition
54 *============================================================================*/
55
56 /* Wall function type */
57 /*--------------------*/
58
59 typedef enum {
60
61 CS_WALL_F_UNSET = -999,
62 CS_WALL_F_DISABLED = 0,
63 CS_WALL_F_1SCALE_POWER = 1,
64 CS_WALL_F_1SCALE_LOG = 2,
65 CS_WALL_F_2SCALES_LOG = 3,
66 CS_WALL_F_SCALABLE_2SCALES_LOG = 4,
67 CS_WALL_F_2SCALES_VDRIEST = 5,
68 CS_WALL_F_2SCALES_SMOOTH_ROUGH = 6,
69 CS_WALL_F_2SCALES_CONTINUOUS = 7,
70
71 } cs_wall_f_type_t;
72
73 typedef enum {
74
75 CS_WALL_F_S_UNSET = -999,
76 CS_WALL_F_S_ARPACI_LARSEN = 0,
77 CS_WALL_F_S_VDRIEST = 1,
78 CS_WALL_F_S_LOUIS = 2,
79 CS_WALL_F_S_MONIN_OBUKHOV = 3,
80 CS_WALL_F_S_SMOOTH_ROUGH = 4,//TODO merge with MO ?
81
82 } cs_wall_f_s_type_t;
83
84 /* Wall functions descriptor */
85 /*---------------------------*/
86
87 typedef struct {
88
89 cs_wall_f_type_t iwallf; /* wall function type */
90
91 cs_wall_f_s_type_t iwalfs; /* wall function type for scalars */
92
93 double ypluli; /* limit value of y+ for the viscous
94 sublayer */
95
96 } cs_wall_functions_t;
97
98 /*============================================================================
99 * Global variables
100 *============================================================================*/
101
102 /* Pointer to wall functions descriptor structure */
103
104 extern const cs_wall_functions_t *cs_glob_wall_functions;
105
106 /*============================================================================
107 * Private function definitions
108 *============================================================================*/
109
110 /*----------------------------------------------------------------------------*/
111 /*!
112 * \brief Power law: Werner & Wengle.
113 *
114 * \param[in] l_visc kinematic viscosity
115 * \param[in] vel wall projected cell center velocity
116 * \param[in] y wall distance
117 * \param[out] iuntur indicator: 0 in the viscous sublayer
118 * \param[out] nsubla counter of cell in the viscous sublayer
119 * \param[out] nlogla counter of cell in the log-layer
120 * \param[out] ustar friction velocity
121 * \param[out] uk friction velocity
122 * \param[out] yplus dimensionless distance to the wall
123 * \param[out] ypup yplus projected vel ratio
124 * \param[out] cofimp \f$\frac{|U_F|}{|U_I^p|}\f$ to ensure a good
125 * turbulence production
126 */
127 /*----------------------------------------------------------------------------*/
128
129 inline static void
cs_wall_functions_1scale_power(cs_real_t l_visc,cs_real_t vel,cs_real_t y,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)130 cs_wall_functions_1scale_power(cs_real_t l_visc,
131 cs_real_t vel,
132 cs_real_t y,
133 int *iuntur,
134 cs_lnum_t *nsubla,
135 cs_lnum_t *nlogla,
136 cs_real_t *ustar,
137 cs_real_t *uk,
138 cs_real_t *yplus,
139 cs_real_t *ypup,
140 cs_real_t *cofimp)
141 {
142 const double ypluli = cs_glob_wall_functions->ypluli;
143
144 const double ydvisc = y / l_visc;
145
146 /* Compute the friction velocity ustar */
147
148 *ustar = pow((vel/(cs_turb_apow * pow(ydvisc, cs_turb_bpow))), cs_turb_dpow);
149 *uk = *ustar;
150 *yplus = *ustar * ydvisc;
151
152 /* In the viscous sub-layer: U+ = y+ */
153 if (*yplus <= ypluli) {
154
155 *ustar = sqrt(vel / ydvisc);
156 *yplus = *ustar * ydvisc;
157 *uk = *ustar;
158 *ypup = 1.;
159 *cofimp = 0.;
160
161 /* Disable the wall funcion count the cell in the viscous sub-layer */
162 *iuntur = 0;
163 *nsubla += 1;
164
165 /* In the log layer */
166 } else {
167 *ypup = pow(vel, 2. * cs_turb_dpow-1.)
168 / pow(cs_turb_apow, 2. * cs_turb_dpow);
169 *cofimp = 1. + cs_turb_bpow
170 * pow(*ustar, cs_turb_bpow + 1. - 1./cs_turb_dpow)
171 * (pow(2., cs_turb_bpow - 1.) - 2.);
172
173 /* Count the cell in the log layer */
174 *nlogla += 1;
175
176 }
177 }
178
179 /*----------------------------------------------------------------------------*/
180 /*!
181 * \brief Log law: piecewise linear and log, with one velocity scale based on
182 * the friction.
183 *
184 * \param[in] ifac face number
185 * \param[in] l_visc kinematic viscosity
186 * \param[in] vel wall projected cell center velocity
187 * \param[in] y wall distance
188 * \param[out] iuntur indicator: 0 in the viscous sublayer
189 * \param[out] nsubla counter of cell in the viscous sublayer
190 * \param[out] nlogla counter of cell in the log-layer
191 * \param[out] ustar friction velocity
192 * \param[out] uk friction velocity
193 * \param[out] yplus dimensionless distance to the wall
194 * \param[out] ypup yplus projected vel ratio
195 * \param[out] cofimp \f$\frac{|U_F|}{|U_I^p|}\f$ to ensure a good
196 * turbulence production
197 */
198 /*----------------------------------------------------------------------------*/
199
200 inline static void
cs_wall_functions_1scale_log(cs_lnum_t ifac,cs_real_t l_visc,cs_real_t vel,cs_real_t y,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)201 cs_wall_functions_1scale_log(cs_lnum_t ifac,
202 cs_real_t l_visc,
203 cs_real_t vel,
204 cs_real_t y,
205 int *iuntur,
206 cs_lnum_t *nsubla,
207 cs_lnum_t *nlogla,
208 cs_real_t *ustar,
209 cs_real_t *uk,
210 cs_real_t *yplus,
211 cs_real_t *ypup,
212 cs_real_t *cofimp)
213 {
214 const double ypluli = cs_glob_wall_functions->ypluli;
215
216 double ustarwer, ustarmin, ustaro, ydvisc;
217 double eps = 0.001;
218 int niter_max = 100;
219 int iter = 0;
220 double reynolds;
221
222 /* Compute the local Reynolds number */
223
224 ydvisc = y / l_visc;
225 reynolds = vel * ydvisc;
226
227 /*
228 * Compute the friction velocity ustar
229 */
230
231 /* In the viscous sub-layer: U+ = y+ */
232 if (reynolds <= ypluli * ypluli) {
233
234 *ustar = sqrt(vel / ydvisc);
235 *yplus = *ustar * ydvisc;
236 *uk = *ustar;
237 *ypup = 1.;
238 *cofimp = 0.;
239
240 /* Disable the wall funcion count the cell in the viscous sub-layer */
241 *iuntur = 0;
242 *nsubla += 1;
243
244 /* In the log layer */
245 } else {
246
247 /* The initial value is Wener or the minimun ustar to ensure convergence */
248 ustarwer = pow(fabs(vel) / cs_turb_apow / pow(ydvisc, cs_turb_bpow),
249 cs_turb_dpow);
250 ustarmin = exp(-cs_turb_cstlog * cs_turb_xkappa)/ydvisc;
251 ustaro = CS_MAX(ustarwer, ustarmin);
252 *ustar = (cs_turb_xkappa * vel + ustaro)
253 / (log(ydvisc * ustaro) + cs_turb_xkappa * cs_turb_cstlog + 1.);
254
255 /* Iterative solving */
256 for (iter = 0; iter < niter_max
257 && fabs(*ustar - ustaro) >= eps * ustaro; iter++) {
258 ustaro = *ustar;
259 *ustar = (cs_turb_xkappa * vel + ustaro)
260 / (log(ydvisc * ustaro) + cs_turb_xkappa * cs_turb_cstlog + 1.);
261 }
262
263 if (iter >= niter_max) {
264 bft_printf(_("WARNING: non-convergence in the computation\n"
265 "******** of the friction velocity\n\n"
266 "face id: %ld \n"
267 "friction vel: %f \n" ), (long)ifac, *ustar);
268 }
269
270 *uk = *ustar;
271 *yplus = *ustar * ydvisc;
272 *ypup = *yplus / (log(*yplus) / cs_turb_xkappa + cs_turb_cstlog);
273 *cofimp = 1. - *ypup / cs_turb_xkappa * 1.5 / *yplus;
274
275 /* Count the cell in the log layer */
276 *nlogla += 1;
277
278 }
279
280 }
281
282 /*----------------------------------------------------------------------------
283 * Compute du+/dy+ for a given yk+.
284 *
285 * parameters:
286 * yplus <-- dimensionless distance
287 *
288 * returns:
289 * the resulting dimensionless velocity.
290 *----------------------------------------------------------------------------*/
291
292 inline static cs_real_t
_uplus(cs_real_t yp,cs_real_t ka,cs_real_t B,cs_real_t cuv,cs_real_t y0,cs_real_t n)293 _uplus(cs_real_t yp,
294 cs_real_t ka,
295 cs_real_t B,
296 cs_real_t cuv,
297 cs_real_t y0,
298 cs_real_t n)
299 {
300 cs_real_t uplus, f_blend;
301
302 f_blend = exp(-0.25*cuv*pow(yp,3));
303 uplus = f_blend*yp + (log(yp)/ka +B)*(1.-exp(-pow(yp/y0,n)))*(1-f_blend);
304
305 return uplus;
306 }
307
308 /*----------------------------------------------------------------------------
309 * Compute du+/dy+ for a given yk+.
310 * parameters:
311 * yplus <-- dimensionless distance
312 * returns:
313 * the resulting dimensionless velocity gradient.
314 *----------------------------------------------------------------------------*/
315
316 inline static cs_real_t
_dupdyp(cs_real_t yp,cs_real_t ka,cs_real_t B,cs_real_t cuv,cs_real_t y0,cs_real_t n)317 _dupdyp(cs_real_t yp,
318 cs_real_t ka,
319 cs_real_t B,
320 cs_real_t cuv,
321 cs_real_t y0,
322 cs_real_t n)
323 {
324 cs_real_t dupdyp;
325
326 dupdyp = exp(-0.25*cuv*pow(yp,3))
327 - 0.75*cuv*pow(yp,3.)*exp(-0.25*cuv*pow(yp,3.))
328 + n*(1.-exp(-0.25*cuv*pow(yp,3.)))*(pow(yp,n-1.)/pow(y0,n))
329 *exp(-pow(yp/y0,n))*((1./ka)*log(yp)+B)
330 + 0.75*cuv*pow(yp,2.)*exp(-0.25*cuv*pow(yp,3.))
331 *(1.-exp(-pow(yp/y0,n)))*((1./ka)*log(yp)+B)
332 + (1./ka/yp)*(1.-exp(-pow(yp/y0,n)))*(1-exp(-0.25*cuv*pow(yp,3.)));
333
334 return dupdyp;
335 }
336
337 /*----------------------------------------------------------------------------*/
338 /*!
339 * \brief Continuous law of the wall between the linear and log law, with two
340 * velocity scales based on the friction and the turbulent kinetic energy. Can
341 * be used with RANS model either in high Reynolds or in low Reynolds
342 * (if the underlying RANS model is compatible).
343 *
344 * \param[in] rnnb \f$\vec{n}.(\tens{R}\vec{n})\f$
345 * \param[in] l_visc kinematic viscosity
346 * \param[in] t_visc turbulent kinematic viscosity
347 * \param[in] vel wall projected cell center velocity
348 * \param[in] y wall distance
349 * \param[in] kinetic_en turbulent kinetic energy
350 * \param[out] iuntur indicator: 0 in the viscous sublayer
351 * \param[out] nsubla counter of cell in the viscous sublayer
352 * \param[out] nlogla counter of cell in the log-layer
353 * \param[out] ustar friction velocity
354 * \param[out] uk friction velocity
355 * \param[out] yplus dimensionless distance to the wall
356 * \param[out] ypup yplus projected vel ratio
357 * \param[out] cofimp \f$\frac{|U_F|}{|U_I^p|}\f$ to ensure a good
358 * turbulence production
359 */
360 /*----------------------------------------------------------------------------*/
361
362 inline static void
cs_wall_functions_2scales_continuous(cs_real_t rnnb,cs_real_t l_visc,cs_real_t t_visc,cs_real_t vel,cs_real_t y,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)363 cs_wall_functions_2scales_continuous(cs_real_t rnnb,
364 cs_real_t l_visc,
365 cs_real_t t_visc,
366 cs_real_t vel,
367 cs_real_t y,
368 cs_real_t kinetic_en,
369 int *iuntur,
370 cs_lnum_t *nsubla,
371 cs_lnum_t *nlogla,
372 cs_real_t *ustar,
373 cs_real_t *uk,
374 cs_real_t *yplus,
375 cs_real_t *ypup,
376 cs_real_t *cofimp)
377 {
378 const double ypluli = cs_glob_wall_functions->ypluli;
379 double Re, g, t_visc_durb;
380 cs_real_t cstcuv, csty0, cstN;
381 cs_real_t dup1, dup2, uplus;
382
383 /* Local constants */
384 cstcuv = 1.0674e-3;
385 csty0 = 14.5e0;
386 cstN = 2.25e0;
387
388 /* Iterative process to determine uk through TKE law */
389 Re = sqrt(kinetic_en) * y / l_visc;
390 g = exp(-Re/11.);
391
392 /* Comutation of uk*/
393 *uk = sqrt( (1.-g) * cs_turb_cmu025 * cs_turb_cmu025 * kinetic_en
394 + g * l_visc * vel/y);
395
396 /* Local value of y+, estimated U+ */
397 *yplus = *uk * y / l_visc;
398 uplus = _uplus( *yplus, cs_turb_xkappa, cs_turb_cstlog, cstcuv, csty0, cstN);
399 /* Deduced velocity sclale uet*/
400 *ustar = vel / uplus;
401
402 if (*yplus < 1.e-1) {
403
404 *ypup = 1.0;
405 *cofimp = 0.0;
406
407 *iuntur = 0;
408 *nsubla += 1;
409
410 }
411 else {
412
413 /* Dimensionless velocity gradient in y+ */
414 dup1 = _dupdyp(*yplus, cs_turb_xkappa, cs_turb_cstlog,
415 cstcuv, csty0, cstN);
416 /* Dimensionless velocity gradient in 2 x y+ */
417 dup2 = _dupdyp(2.0 * *yplus, cs_turb_xkappa,
418 cs_turb_cstlog, cstcuv, csty0, cstN);
419
420 *ypup = *yplus / uplus;
421
422 /* ------------------------------------------------------------
423 * Cofimp = U,F/U,I is built so that the theoretical expression
424 * of the production P_theo = dup1 * (1.0 - dup1) is equal to
425 * P_calc = mu_t,I * ((U,I - U,F + IF*dup2)/(2IF) )^2
426 * This is a generalization of the process implemented in the 2
427 * scales wall function (iwallf = 3).
428 * ------------------------------------------------------------*/
429
430 /* Turbulent viscocity is modified for RSM so that its expression
431 * remain valid down to the wall, according to Durbin :
432 * nu_t = 0.22 * v'2 * k / eps */
433 const cs_turb_model_t *turb_model = cs_get_glob_turb_model();
434 assert(turb_model != NULL);
435 if (turb_model->itytur == 3)
436 t_visc_durb = t_visc / (kinetic_en * cs_turb_cmu ) * rnnb * 0.22;
437 else
438 t_visc_durb = t_visc;
439
440 *cofimp
441 = 1. - *ypup * (2. * sqrt(l_visc / t_visc_durb * dup1 * (1. - dup1))
442 - dup2);
443
444 /* log layer */
445 if (*yplus > ypluli) {
446 *nlogla += 1;
447 /* viscous sub-layer or buffer layer*/
448 } else {
449 *iuntur = 0;
450 *nsubla += 1;
451 }
452 }
453 }
454
455 /*----------------------------------------------------------------------------*/
456 /*!
457 * \brief Log law: piecewise linear and log, with two velocity scales based on
458 * the friction and the turbulent kinetic energy.
459 *
460 * \param[in] l_visc kinematic viscosity
461 * \param[in] t_visc turbulent kinematic viscosity
462 * \param[in] vel wall projected cell center velocity
463 * \param[in] y wall distance
464 * \param[in] kinetic_en turbulent kinetic energy
465 * \param[out] iuntur indicator: 0 in the viscous sublayer
466 * \param[out] nsubla counter of cell in the viscous sublayer
467 * \param[out] nlogla counter of cell in the log-layer
468 * \param[out] ustar friction velocity
469 * \param[out] uk friction velocity
470 * \param[out] yplus dimensionless distance to the wall
471 * \param[out] ypup yplus projected vel ratio
472 * \param[out] cofimp \f$\frac{|U_F|}{|U_I^p|}\f$ to ensure a good
473 * turbulence production
474 */
475 /*----------------------------------------------------------------------------*/
476
477 inline static void
cs_wall_functions_2scales_log(cs_real_t l_visc,cs_real_t t_visc,cs_real_t vel,cs_real_t y,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)478 cs_wall_functions_2scales_log(cs_real_t l_visc,
479 cs_real_t t_visc,
480 cs_real_t vel,
481 cs_real_t y,
482 cs_real_t kinetic_en,
483 int *iuntur,
484 cs_lnum_t *nsubla,
485 cs_lnum_t *nlogla,
486 cs_real_t *ustar,
487 cs_real_t *uk,
488 cs_real_t *yplus,
489 cs_real_t *ypup,
490 cs_real_t *cofimp)
491 {
492 const double ypluli = cs_glob_wall_functions->ypluli;
493
494 double rcprod, ml_visc, Re, g;
495
496 /* Compute the friction velocity ustar */
497
498 /* Blending for very low values of k */
499 Re = sqrt(kinetic_en) * y / l_visc;
500 g = exp(-Re/11.);
501
502 *uk = sqrt( (1.-g) * cs_turb_cmu025 * cs_turb_cmu025 * kinetic_en
503 + g * l_visc * vel / y);
504
505 *yplus = *uk * y / l_visc;
506
507 /* log layer */
508 if (*yplus > ypluli) {
509
510 *ustar = vel / (log(*yplus) / cs_turb_xkappa + cs_turb_cstlog);
511 *ypup = *yplus / (log(*yplus) / cs_turb_xkappa + cs_turb_cstlog);
512 /* Mixing length viscosity */
513 ml_visc = cs_turb_xkappa * l_visc * *yplus;
514 rcprod = CS_MIN(cs_turb_xkappa, CS_MAX(1., sqrt(ml_visc / t_visc)) / *yplus);
515 *cofimp = 1. - *ypup / cs_turb_xkappa * ( 2. * rcprod - 1. / (2. * *yplus));
516
517 *nlogla += 1;
518
519 /* viscous sub-layer */
520 } else {
521
522 if (*yplus > 1.e-12) {
523 *ustar = fabs(vel / *yplus); /* FIXME remove that: its is here only to
524 be fully equivalent to the former code. */
525 } else {
526 *ustar = 0.;
527 }
528 *ypup = 1.;
529 *cofimp = 0.;
530
531 *iuntur = 0;
532 *nsubla += 1;
533
534 }
535 }
536
537 /*----------------------------------------------------------------------------*/
538 /*!
539 * \brief Scalable wall function: shift the wall if \f$ y^+ < y^+_{lim} \f$.
540 *
541 * \param[in] l_visc kinematic viscosity
542 * \param[in] t_visc turbulent kinematic viscosity
543 * \param[in] vel wall projected cell center velocity
544 * \param[in] y wall distance
545 * \param[in] kinetic_en turbulent kinetic energy
546 * \param[out] iuntur indicator: 0 in the viscous sublayer
547 * \param[out] nsubla counter of cell in the viscous sublayer
548 * \param[out] nlogla counter of cell in the log-layer
549 * \param[out] ustar friction velocity
550 * \param[out] uk friction velocity
551 * \param[out] yplus dimensionless distance to the wall
552 * \param[out] dplus dimensionless shift to the wall for scalable
553 * wall functions
554 * \param[out] ypup yplus projected vel ratio
555 * \param[out] cofimp \f$\frac{|U_F|}{|U_I^p|}\f$ to ensure a good
556 * turbulence production
557 */
558 /*----------------------------------------------------------------------------*/
559
560 inline static void
cs_wall_functions_2scales_scalable(cs_real_t l_visc,cs_real_t t_visc,cs_real_t vel,cs_real_t y,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 * dplus,cs_real_t * ypup,cs_real_t * cofimp)561 cs_wall_functions_2scales_scalable(cs_real_t l_visc,
562 cs_real_t t_visc,
563 cs_real_t vel,
564 cs_real_t y,
565 cs_real_t kinetic_en,
566 int *iuntur,
567 cs_lnum_t *nsubla,
568 cs_lnum_t *nlogla,
569 cs_real_t *ustar,
570 cs_real_t *uk,
571 cs_real_t *yplus,
572 cs_real_t *dplus,
573 cs_real_t *ypup,
574 cs_real_t *cofimp)
575 {
576 CS_UNUSED(iuntur);
577
578 const double ypluli = cs_glob_wall_functions->ypluli;
579
580 double rcprod, ml_visc, Re, g;
581 /* Compute the friction velocity ustar */
582
583 /* Blending for very low values of k */
584 Re = sqrt(kinetic_en) * y / l_visc;
585 g = exp(-Re/11.);
586
587 *uk = sqrt( (1.-g) * cs_turb_cmu025 * cs_turb_cmu025 * kinetic_en
588 + g * l_visc * vel / y);
589
590 *yplus = *uk * y / l_visc;
591
592 /* Compute the friction velocity ustar */
593 *uk = cs_turb_cmu025 * sqrt(kinetic_en);//FIXME
594 *yplus = *uk * y / l_visc;//FIXME
595
596 /* Log layer */
597 if (*yplus > ypluli) {
598
599 *dplus = 0.;
600
601 *nlogla += 1;
602
603 /* Viscous sub-layer and therefore shift */
604 } else {
605
606 *dplus = ypluli - *yplus;
607
608 /* Count the cell as if it was in the viscous sub-layer */
609 *nsubla += 1;
610
611 }
612
613 /* Mixing length viscosity */
614 ml_visc = cs_turb_xkappa * l_visc * (*yplus + *dplus);
615 rcprod = CS_MIN(cs_turb_xkappa, CS_MAX(1., sqrt(ml_visc / t_visc)) / (*yplus + *dplus));
616
617 *ustar = vel / (log(*yplus + *dplus) / cs_turb_xkappa + cs_turb_cstlog);
618 *ypup = *yplus / (log(*yplus + *dplus) / cs_turb_xkappa + cs_turb_cstlog);
619 *cofimp = 1. - *ypup
620 / cs_turb_xkappa * (2. * rcprod - 1. / (2. * *yplus + *dplus));
621 }
622
623 /*----------------------------------------------------------------------------
624 * Compute u+ for a given yk+ between 0.1 and 200 according to the two
625 * scales wall functions using Van Driest mixing length.
626 * This function holds the coefficients of the polynome fitting log(u+).
627 *
628 * parameters:
629 * yplus <-- dimensionless distance
630 *
631 * returns:
632 * the resulting dimensionless velocity.
633 *----------------------------------------------------------------------------*/
634
635 inline static cs_real_t
_vdriest_dupdyp_integral(cs_real_t yplus)636 _vdriest_dupdyp_integral(cs_real_t yplus)
637 {
638 /* Coefficients of the polynome fitting log(u+) for yk < 200 */
639 static double aa[11] = {-0.0091921, 3.9577, 0.031578,
640 -0.51013, -2.3254, -0.72665,
641 2.969, 0.48506, -1.5944,
642 0.087309, 0.1987 };
643
644 cs_real_t y1,y2,y3,y4,y5,y6,y7,y8,y9,y10, uplus;
645
646 y1 = 0.25 * log(yplus);
647 y2 = y1 * y1;
648 y3 = y2 * y1;
649 y4 = y3 * y1;
650 y5 = y4 * y1;
651 y6 = y5 * y1;
652 y7 = y6 * y1;
653 y8 = y7 * y1;
654 y9 = y8 * y1;
655 y10 = y9 * y1;
656
657 uplus = aa[0]
658 + aa[1] * y1
659 + aa[2] * y2
660 + aa[3] * y3
661 + aa[4] * y4
662 + aa[5] * y5
663 + aa[6] * y6
664 + aa[7] * y7
665 + aa[8] * y8
666 + aa[9] * y9
667 + aa[10] * y10;
668
669 return exp(uplus);
670 }
671
672 /*----------------------------------------------------------------------------*/
673 /*!
674 * \brief Two velocity scales wall function using Van Driest mixing length.
675 *
676 * \f$ u^+ \f$ is computed as follows:
677 * \f[ u^+ = \int_0^{y_k^+} \dfrac{dy_k^+}{1+L_m^k} \f]
678 * with \f$ L_m^k \f$ standing for Van Driest mixing length:
679 * \f[ L_m^k = \kappa y_k^+ (1 - exp(\dfrac{-y_k^+}{A})) \f].
680 *
681 * A polynome fitting the integral is used for \f$ y_k^+ < 200 \f$,
682 * and a log law is used for \f$ y_k^+ >= 200 \f$.
683 *
684 * A wall roughness can be taken into account in the mixing length as
685 * proposed by Rotta (1962) with Cebeci & Chang (1978) correlation.
686 *
687 * \param[in] rnnb \f$\vec{n}.(\tens{R}\vec{n})\f$
688 * \param[in] l_visc kinematic viscosity
689 * \param[in] vel wall projected cell center velocity
690 * \param[in] y wall distance
691 * \param[in] kinetic_en turbulent kinetic energy
692 * \param[out] iuntur indicator: 0 in the viscous sublayer
693 * \param[out] nsubla counter of cell in the viscous sublayer
694 * \param[out] nlogla counter of cell in the log-layer
695 * \param[out] ustar friction velocity
696 * \param[out] uk friction velocity
697 * \param[out] yplus dimensionless distance to the wall
698 * \param[out] ypup yplus projected vel ratio
699 * \param[out] cofimp \f$\frac{|U_F|}{|U_I^p|}\f$ to ensure a good
700 * turbulence production
701 * \param[in] lmk dimensionless mixing length
702 * \param[in] kr wall roughness
703 * \param[in] wf enable full wall function computation,
704 * if false, uk is not recomputed and uplus is the
705 * only output
706 */
707 /*----------------------------------------------------------------------------*/
708
709 inline static void
cs_wall_functions_2scales_vdriest(cs_real_t rnnb,cs_real_t l_visc,cs_real_t vel,cs_real_t y,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 * lmk,cs_real_t kr,bool wf)710 cs_wall_functions_2scales_vdriest(cs_real_t rnnb,
711 cs_real_t l_visc,
712 cs_real_t vel,
713 cs_real_t y,
714 cs_real_t kinetic_en,
715 int *iuntur,
716 cs_lnum_t *nsubla,
717 cs_lnum_t *nlogla,
718 cs_real_t *ustar,
719 cs_real_t *uk,
720 cs_real_t *yplus,
721 cs_real_t *ypup,
722 cs_real_t *cofimp,
723 cs_real_t *lmk,
724 cs_real_t kr,
725 bool wf)
726 {
727 double urplus, d_up, lmk15;
728
729 if (wf)
730 *uk = sqrt(sqrt((1.-cs_turb_crij2)/cs_turb_crij1 * rnnb * kinetic_en));
731
732 /* Set a low threshold value in case tangential velocity is zero */
733 *yplus = CS_MAX(*uk * y / l_visc, 1.e-4);
734
735 /* Dimensionless roughness */
736 cs_real_t krp = *uk * kr / l_visc;
737
738 /* Extension of Van Driest mixing length according to Rotta (1962) with
739 Cebeci & Chang (1978) correlation */
740 cs_real_t dyrp = 0.9 * (sqrt(krp) - krp * exp(-krp / 6.));
741 cs_real_t yrplus = *yplus + dyrp;
742
743 if (dyrp <= 1.e-1)
744 d_up = dyrp;
745 else if (dyrp <= 200.)
746 d_up = _vdriest_dupdyp_integral(dyrp);
747 else
748 d_up = 16.088739022054590 + log(dyrp/200.) / cs_turb_xkappa;
749
750 if (yrplus <= 1.e-1) {
751
752 urplus = yrplus;
753
754 if (wf) {
755 *iuntur = 0;
756 *nsubla += 1;
757
758 *lmk = 0.;
759
760 *ypup = 1.;
761
762 *cofimp = 0.;
763 }
764
765 } else if (yrplus <= 200.) {
766
767 urplus = _vdriest_dupdyp_integral(yrplus);
768
769 if (wf) {
770 *nlogla += 1;
771
772 *ypup = *yplus / (urplus-d_up);
773
774 /* Mixing length in y+ */
775 *lmk = cs_turb_xkappa * (*yplus) *(1-exp(- (*yplus) / cs_turb_vdriest));
776
777 /* Mixing length in 3/2*y+ */
778 lmk15 = cs_turb_xkappa * 1.5 * (*yplus) *(1-exp(- 1.5 * (*yplus)
779 / cs_turb_vdriest));
780
781 *cofimp = 1. - (2. / (1. + *lmk) - 1. / (1. + lmk15)) * *ypup;
782 }
783
784 } else {
785
786 urplus = 16.088739022054590 + log(yrplus/200) / cs_turb_xkappa;
787
788 if (wf) {
789 *nlogla += 1;
790
791 *ypup = *yplus / (urplus-d_up);
792
793 /* Mixing length in y+ */
794 *lmk = cs_turb_xkappa * (*yplus) *(1-exp(- (*yplus) / cs_turb_vdriest));
795
796 /* Mixing length in 3/2*y+ */
797 lmk15 = cs_turb_xkappa * 1.5 * (*yplus) *(1-exp(- 1.5 * (*yplus)
798 / cs_turb_vdriest));
799
800 *cofimp = 1. - (2. / *lmk - 1. / lmk15) * *ypup;
801 }
802
803 }
804
805 *ustar = vel / (urplus-d_up);
806 }
807
808 /*----------------------------------------------------------------------------*/
809 /*!
810 * \brief Two velocity scales wall function with automatic switch
811 * from rough to smooth.
812 *
813 * \f$ u^+ \f$ is computed as follows:
814 * \f[ u^+ = \dfrac{1}{\kappa}
815 * \ln \left(\dfrac{(y+y_0) u_k}{\nu + \alpha \xi u_k} \right)
816 * + Cst_{smooth} \f]
817 * with \f$ \alpha = \exp \left(- \kappa(Cst_{rough}-Cst_{smooth})\right)
818 * \simeq 0.26 \f$
819 * and \f$ y_0 = \alpha \xi \exp \left(-\kappa Cst_{smooth} \right)
820 * = \xi \exp \left(-\kappa Cst_{rough} \right)
821 * \simeq \dfrac{\xi}{33}\f$.
822 *
823 * \param[in] l_visc kinematic viscosity
824 * \param[in] t_visc turbulent kinematic viscosity
825 * \param[in] vel wall projected cell center velocity
826 * \param[in] y wall distance
827 * \param[in] rough_d roughness length scale (not sand grain)
828 * \param[in] kinetic_en turbulent kinetic energy
829 * \param[out] iuntur indicator: 0 in the viscous sublayer
830 * \param[out] nsubla counter of cell in the viscous sublayer
831 * \param[out] nlogla counter of cell in the log-layer
832 * \param[out] ustar friction velocity
833 * \param[out] uk friction velocity
834 * \param[out] yplus dimensionless distance to the wall
835 * \param[out] dplus dimensionless shift to the wall for scalable
836 * wall functions
837 * \param[out] ypup yplus projected vel ratio
838 * \param[out] cofimp \f$\frac{|U_F|}{|U_I^p|}\f$ to ensure a good
839 * turbulence production
840 */
841 /*----------------------------------------------------------------------------*/
842
843 inline static void
cs_wall_functions_2scales_smooth_rough(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 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 * dplus,cs_real_t * ypup,cs_real_t * cofimp)844 cs_wall_functions_2scales_smooth_rough(cs_real_t l_visc,
845 cs_real_t t_visc,
846 cs_real_t vel,
847 cs_real_t y,
848 cs_real_t rough_d,
849 cs_real_t kinetic_en,
850 int *iuntur,
851 cs_lnum_t *nsubla,
852 cs_lnum_t *nlogla,
853 cs_real_t *ustar,
854 cs_real_t *uk,
855 cs_real_t *yplus,
856 cs_real_t *dplus,
857 cs_real_t *ypup,
858 cs_real_t *cofimp)
859 {
860 CS_UNUSED(iuntur);
861
862 const double ypluli = cs_glob_wall_functions->ypluli;
863
864 double rcprod, ml_visc, Re, g;
865
866 /* Compute the friction velocity ustar */
867
868 /* Shifting of the wall distance to be consistant with
869 * the fully rough wall function
870 *
871 * ln((y+y0)/y0) = ln((y+y0)/alpha xi) + kappa * 5.2
872 *
873 * y0 = xi * exp(-kappa * 8.5)
874 * where xi is the sand grain roughness here
875 * y0 = alpha * xi * exp(-kappa * 5.2)
876 *
877 * so:
878 * alpha = exp(-kappa * (8.5 - 5.2)) = 0.25
879 *
880 */
881 cs_real_t y0 = rough_d;
882 /* Note : Sand grain roughness given by:
883 cs_real_t sg_rough = rough_d * exp(cs_turb_xkappa*cs_turb_cstlog_rough);
884 */
885
886 /* Blending for very low values of k */
887 Re = sqrt(kinetic_en) * (y + y0) / l_visc;
888 g = exp(-Re/11.);
889
890 *uk = sqrt( (1.-g) * cs_turb_cmu025 * cs_turb_cmu025 * kinetic_en
891 + g * l_visc * vel / (y + y0));
892
893
894 *yplus = *uk * y / l_visc;
895
896 /* As for scalable wall functions, yplus is shifted of "dplus" */
897 *dplus = *uk * y0 / l_visc;
898
899 /* Shift of the velocity profile due to roughness */
900 cs_real_t shift_vel = -log(1. + y0 * exp(cs_turb_xkappa * cs_turb_cstlog) * *uk/l_visc)
901 / cs_turb_xkappa;
902
903 /* Log layer and shifted with the roughness */
904 if (*yplus > ypluli) {
905
906 *nlogla += 1;
907
908 }
909
910 /* Viscous sub-layer and therefore shift again */
911 else {
912
913 *dplus = ypluli - *yplus;
914 /* Count the cell as if it was in the viscous sub-layer */
915 *nsubla += 1;
916
917 }
918
919 cs_real_t uplus = log(*yplus + *dplus) / cs_turb_xkappa + cs_turb_cstlog + shift_vel;
920 *ustar = vel / uplus;
921 #if 0
922 bft_printf("uet=%f, u=%f, uplus=%f, yk=%f, duplus=%f\n", *ustar, vel, uplus, *yplus, 1./uplus);
923 #endif
924 *ypup = *yplus / uplus;
925
926 /* Mixing length viscosity, compatible with both regimes */
927 ml_visc = cs_turb_xkappa * l_visc * (*yplus + *dplus);
928 rcprod = CS_MIN(cs_turb_xkappa, CS_MAX(1., sqrt(ml_visc / t_visc)) / (*yplus + *dplus));
929 *cofimp = 1. - *yplus / (cs_turb_xkappa * uplus)
930 * ( 2. * rcprod - 1. / (2. * *yplus + *dplus));
931
932 }
933
934 /*----------------------------------------------------------------------------*/
935 /*!
936 * \brief No wall function.
937 *
938 * \param[in] l_visc kinematic viscosity
939 * \param[in] t_visc turbulent kinematic viscosity
940 * \param[in] vel wall projected cell center velocity
941 * \param[in] y wall distance
942 * \param[out] iuntur indicator: 0 in the viscous sublayer
943 * \param[out] nsubla counter of cell in the viscous sublayer
944 * \param[out] nlogla counter of cell in the log-layer
945 * \param[out] ustar friction velocity
946 * \param[out] uk friction velocity
947 * \param[out] yplus dimensionless distance to the wall
948 * \param[out] dplus dimensionless shift to the wall for scalable
949 * wall functions
950 * \param[out] ypup yplus projected vel ratio
951 * \param[out] cofimp \f$\frac{|U_F|}{|U_I^p|}\f$ to ensure a good
952 * turbulence production
953 */
954 /*----------------------------------------------------------------------------*/
955
956 inline static void
cs_wall_functions_disabled(cs_real_t l_visc,cs_real_t t_visc,cs_real_t vel,cs_real_t y,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 * dplus,cs_real_t * ypup,cs_real_t * cofimp)957 cs_wall_functions_disabled(cs_real_t l_visc,
958 cs_real_t t_visc,
959 cs_real_t vel,
960 cs_real_t y,
961 int *iuntur,
962 cs_lnum_t *nsubla,
963 cs_lnum_t *nlogla,
964 cs_real_t *ustar,
965 cs_real_t *uk,
966 cs_real_t *yplus,
967 cs_real_t *dplus,
968 cs_real_t *ypup,
969 cs_real_t *cofimp)
970 {
971 CS_UNUSED(t_visc);
972 CS_UNUSED(nlogla);
973 CS_UNUSED(dplus);
974
975 const double ypluli = cs_glob_wall_functions->ypluli;
976
977 /* Compute the friction velocity ustar */
978
979 *ustar = sqrt(vel * l_visc / y);
980 *yplus = *ustar * y / l_visc;
981 *uk = *ustar;
982 *ypup = 1.;
983 *cofimp = 0.;
984 *iuntur = 0;
985
986 if (*yplus <= ypluli) {
987
988 /* Disable the wall funcion count the cell in the viscous sub-layer */
989 *nsubla += 1;
990
991 } else {
992
993 /* Count the cell as if it was in the viscous sub-layer */
994 *nsubla += 1;
995
996 }
997 }
998
999 /*----------------------------------------------------------------------------*/
1000 /*!
1001 * \brief The correction of the exchange coefficient is computed thanks to a
1002 * similarity model between dynamic viscous sub-layer and themal sub-layer.
1003 *
1004 * \f$ T^+ \f$ is computed as follows:
1005 *
1006 * - For a laminar Prandtl number smaller than 0.1 (such as liquid metals),
1007 * the standard model with two sub-layers (Prandtl-Taylor) is used.
1008 *
1009 * - For a laminar Prandtl number larger than 0.1 (such as liquids and gaz),
1010 * a model with three sub-layers (Arpaci-Larsen) is used.
1011 *
1012 * The final exchange coefficient is:
1013 * \f[
1014 * h = \dfrac{K}{\centip \centf} h_{tur}
1015 * \f]
1016 *
1017 * \param[in] l_visc kinetic viscosity
1018 * \param[in] prl laminar Prandtl number
1019 * \param[in] prt turbulent Prandtl number
1020 * \param[in] rough_t scalar roughness length scale
1021 * \param[in] uk velocity scale based on TKE
1022 * \param[in] yplus dimensionless distance to the wall
1023 * \param[out] dplus dimensionless shift to the wall for scalable
1024 * wall functions
1025 * \param[out] htur correction for the exchange coefficient
1026 * (\f$ Pr y^+/T^+ \f$)
1027 * \param[out] yplim value of the limit for \f$ y^+ \f$
1028 */
1029 /*----------------------------------------------------------------------------*/
1030
1031 inline static void
cs_wall_functions_s_arpaci_larsen(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)1032 cs_wall_functions_s_arpaci_larsen(cs_real_t l_visc,
1033 cs_real_t prl,
1034 cs_real_t prt,
1035 cs_real_t rough_t,
1036 cs_real_t uk,
1037 cs_real_t yplus,
1038 cs_real_t dplus,
1039 cs_real_t *htur,
1040 cs_real_t *yplim)
1041 {
1042 /* Local variables */
1043 double tplus;
1044 double beta2,a2;
1045 double yp2;
1046 double prlm1;
1047
1048 const double epzero = cs_math_epzero;
1049
1050 /*==========================================================================
1051 1. Initializations
1052 ==========================================================================*/
1053
1054 (*htur) = prl * CS_MAX(yplus,epzero)/CS_MAX(yplus+dplus,epzero);
1055
1056 prlm1 = 0.1;
1057
1058 /* Sand grain roughness is:
1059 * zeta = z0 * exp(kappa 8.5)
1060 * Then:
1061 * hp = zeta uk / nu * exp( -kappa(8.5 - 5.2))
1062 * = z0 * uk / nu * exp(kappa * 5.2)
1063 * where 5.2 is the smooth log constant, and 8.5 the rough one
1064 *
1065 * FIXME check if we should use a molecular Schmidt number
1066 */
1067 cs_real_t hp = rough_t *uk / l_visc * exp(cs_turb_xkappa*cs_turb_cstlog);
1068
1069 /* Shift of the temperature profile due to roughness */
1070 cs_real_t shift_temp = -log(1. + hp);
1071
1072 /*==========================================================================
1073 2. Compute htur for small Prandtl numbers
1074 ==========================================================================*/
1075
1076 if (prl <= prlm1) {
1077 (*yplim) = prt/(prl*cs_turb_xkappa);
1078 if (yplus > (*yplim)) {
1079 tplus = prl*(*yplim) + prt/cs_turb_xkappa * (log((yplus+dplus)/(*yplim)) + shift_temp);
1080 (*htur) = prl * yplus / tplus;
1081 }
1082
1083 /*========================================================================
1084 3. Compute htur for the model with three sub-layers
1085 ========================================================================*/
1086
1087 } else {
1088 yp2 = sqrt(cs_turb_xkappa*1000./prt);
1089 (*yplim) = pow(1000./prl, 1./3.);
1090
1091 a2 = 15.*pow(prl, 2./3.);
1092
1093 if (yplus >= (*yplim) && yplus < yp2) {
1094 tplus = a2 - 500./((yplus+dplus)*(yplus+dplus));
1095 (*htur) = prl * yplus / tplus;
1096 }
1097
1098 if (yplus >= yp2) {
1099 beta2 = a2 - 0.5 * prt /cs_turb_xkappa;
1100 tplus = beta2 + prt/cs_turb_xkappa*log((yplus+dplus)/yp2);
1101 (*htur) = prl * yplus / tplus;
1102 }
1103
1104 }
1105 }
1106
1107 /*----------------------------------------------------------------------------*/
1108 /*!
1109 * \brief The correction of the exchange coefficient
1110 * \f$ h_{tur} = \sigma \dfrac{y^+}{t^+} \f$ is computed thanks to a
1111 * numerical integration of:
1112 * \f[
1113 * \dfrac{T^+}{\sigma}
1114 * = \int_0^{y_k^+} \dfrac{dy_k^+}{1+\dfrac{\sigma}{\sigma_t}\nu_t^+}
1115 * \f]
1116 * with \f$ \nu_t^+ = L_m^k \f$ as assumed in the derivation of the two scales
1117 * wall function using Van Driest mixing length.
1118 * Therefore \f$ L_m^k = \kappa y_k^+(1 - exp(\dfrac{-y_k^+}{A})) \f$ is taken.
1119 *
1120 * Notice that we integrate up to \f$ y^+=100 \f$ (YP100), beyond that value
1121 * the profile is prolonged by a logarithm relying on the fact that
1122 * \f$ \nu_t^+=\kappa y^+ \f$ beyond \f$ y^+=100 \f$.
1123 *
1124 * \param[in] prl molecular Prandtl number
1125 * ( \f$ Pr=\sigma=\frac{\mu C_p}{\lambda_f} \f$ )
1126 * \param[in] prt turbulent Prandtl number ( \f$ \sigma_t \f$ )
1127 * \param[in] yplus dimensionless distance to the wall
1128 * \param[out] htur correction for the exchange coefficient
1129 * (\f$ Pr y^+/T^+ \f$)
1130 */
1131 /*----------------------------------------------------------------------------*/
1132
1133 inline static void
cs_wall_functions_s_vdriest(cs_real_t prl,cs_real_t prt,cs_real_t yplus,cs_real_t * htur)1134 cs_wall_functions_s_vdriest(cs_real_t prl,
1135 cs_real_t prt,
1136 cs_real_t yplus,
1137 cs_real_t *htur)
1138 {
1139 cs_real_t prlrat = prl / prt;
1140
1141 /* Parameters of the numerical quadrature */
1142 const int ninter_max = 100;
1143 const cs_real_t ypmax = 1.e2;
1144
1145 /* No correction for very small yplus */
1146 if (yplus <= 0.1)
1147 *htur = 1.;
1148 else {
1149 cs_real_t ypint = CS_MIN(yplus, ypmax);
1150
1151 /* The number of sub-intervals is taken proportional to yplus and equal to
1152 * ninter_max if yplus=ypmax */
1153
1154 int npeff = CS_MAX((int)(ypint / ypmax * (double)(ninter_max)), 1);
1155
1156 double dy = ypint / (double)(npeff);
1157 cs_real_t stplus = 0.;
1158 cs_real_t nut1 = 0.;
1159 cs_real_t nut2 = 0.;
1160
1161 for (int ip = 1; ip <= npeff; ip++) {
1162 double yp = ypint * (double)(ip) / (double)(npeff);
1163 nut2 = cs_turb_xkappa * yp * (1. - exp(-yp / cs_turb_vdriest));
1164 stplus += dy / (1. + prlrat * 0.5 * (nut1 + nut2));
1165 nut1 = nut2;
1166 }
1167
1168 if (yplus > ypint) {
1169 cs_real_t r = prlrat * cs_turb_xkappa;
1170 stplus += log( (1. + r*yplus) / (1. + r*ypint)) / r;
1171 }
1172
1173 if (stplus >= 1.e-6)
1174 *htur = yplus / stplus;
1175 else
1176 *htur = 1.;
1177 }
1178 }
1179
1180 /*----------------------------------------------------------------------------*/
1181 /*!
1182 * \brief Rough Smooth Thermal Wall Function - Prototype
1183 *
1184 * \param[in] l_visc kinetic viscosity
1185 * \param[in] prl molecular Prandtl number
1186 * ( \f$ Pr=\sigma=\frac{\mu C_p}{\lambda_f} \f$ )
1187 * \param[in] prt turbulent Prandtl number ( \f$ \sigma_t \f$ )
1188 * \param[in] rough_t scalar roughness length scale
1189 * \param[in] uk velocity scale based on TKE
1190 * \param[in] yplus dimensionless distance to the wall
1191 * \param[in] dplus dimensionless shift to the wall for scalable
1192 * wall function
1193 * \param[out] htur correction for the exchange coefficient
1194 * (\f$ Pr y^+/T^+ \f$)
1195 */
1196 /*----------------------------------------------------------------------------*/
1197
1198 inline static void
cs_wall_functions_s_smooth_rough(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)1199 cs_wall_functions_s_smooth_rough(cs_real_t l_visc,
1200 cs_real_t prl,
1201 cs_real_t prt,
1202 cs_real_t rough_t,
1203 cs_real_t uk,
1204 cs_real_t yplus,
1205 cs_real_t dplus,
1206 cs_real_t *htur)
1207 {
1208 CS_UNUSED(prt);
1209
1210 /* Sand grain roughness is:
1211 * zeta = z0 * exp(kappa 8.5)
1212 * Then:
1213 * hp = zeta uk / nu * exp( -kappa(8.5 - 5.2))
1214 * = z0 * uk / nu * exp(kappa * 5.2)
1215 * where 5.2 is the smooth log constant, and 8.5 the rough one
1216 *
1217 * FIXME check if we should use a molecular Schmidt number
1218 */
1219 cs_real_t hp = rough_t *uk / l_visc * exp(cs_turb_xkappa*cs_turb_cstlog);
1220 const double ypluli = cs_glob_wall_functions->ypluli;
1221 const double epzero = cs_math_epzero;
1222
1223 (*htur) = CS_MAX(yplus,epzero)/CS_MAX(yplus+dplus,epzero);
1224
1225 /* Shift of the temperature profile due to roughness */
1226 cs_real_t shift_temp = -log(1. + hp);
1227
1228 if (yplus > ypluli) {
1229 cs_real_t tplus = prt * ((log(yplus+dplus) + shift_temp)/cs_turb_xkappa + cs_turb_cstlog);
1230 (*htur) = prl * yplus / tplus;
1231 }
1232 }
1233
1234 /*============================================================================
1235 * Public function definitions for Fortran API
1236 *============================================================================*/
1237
1238 /*----------------------------------------------------------------------------
1239 * Wrapper to cs_wall_functions_velocity.
1240 *----------------------------------------------------------------------------*/
1241
1242 void CS_PROCF (wallfunctions, WALLFUNCTIONS)
1243 (
1244 const int *const iwallf,
1245 const cs_lnum_t *const ifac,
1246 const cs_real_t *const viscosity,
1247 const cs_real_t *const t_visc,
1248 const cs_real_t *const vel,
1249 const cs_real_t *const y,
1250 const cs_real_t *const rough_d,
1251 const cs_real_t *const rnnb,
1252 const cs_real_t *const kinetic_en,
1253 int *iuntur,
1254 cs_lnum_t *nsubla,
1255 cs_lnum_t *nlogla,
1256 cs_real_t *ustar,
1257 cs_real_t *uk,
1258 cs_real_t *yplus,
1259 cs_real_t *ypup,
1260 cs_real_t *cofimp,
1261 cs_real_t *dplus
1262 );
1263
1264 /*----------------------------------------------------------------------------
1265 * Wrapper to cs_wall_functions_scalar.
1266 *----------------------------------------------------------------------------*/
1267
1268 void CS_PROCF (hturbp, HTURBP)
1269 (
1270 const int *const iwalfs,
1271 const cs_real_t *const l_visc,
1272 const cs_real_t *const prl,
1273 const cs_real_t *const prt,
1274 const cs_real_t *const rough_t,
1275 const cs_real_t *const uk,
1276 const cs_real_t *const yplus,
1277 const cs_real_t *const dplus,
1278 cs_real_t *htur,
1279 cs_real_t *yplim
1280 );
1281
1282 /*=============================================================================
1283 * Public function prototypes
1284 *============================================================================*/
1285
1286 /*----------------------------------------------------------------------------
1287 *! \brief Provide access to cs_glob_wall_functions
1288 *----------------------------------------------------------------------------*/
1289
1290 cs_wall_functions_t *
1291 cs_get_glob_wall_functions(void);
1292
1293 /*----------------------------------------------------------------------------*/
1294 /*!
1295 * \brief Compute the friction velocity and \f$y^+\f$ / \f$u^+\f$.
1296 *
1297 * \param[in] iwallf wall function type
1298 * \param[in] ifac face number
1299 * \param[in] l_visc kinematic viscosity
1300 * \param[in] t_visc turbulent kinematic viscosity
1301 * \param[in] vel wall projected
1302 * cell center velocity
1303 * \param[in] y wall distance
1304 * \param[in] rough_d roughness lenghth scale
1305 * \param[in] rnnb \f$\vec{n}.(\tens{R}\vec{n})\f$
1306 * \param[in] kinetic_en turbulent kinetic energy
1307 * \param[in] iuntur indicator:
1308 * 0 in the viscous sublayer
1309 * \param[in] nsubla counter of cell in the viscous
1310 * sublayer
1311 * \param[in] nlogla counter of cell in the log-layer
1312 * \param[out] ustar friction velocity
1313 * \param[out] uk friction velocity
1314 * \param[out] yplus non-dimension wall distance
1315 * \param[out] ypup yplus projected vel ratio
1316 * \param[out] cofimp \f$\frac{|U_F|}{|U_I^p|}\f$ to ensure a good
1317 * turbulence production
1318 * \param[out] dplus dimensionless shift to the wall
1319 * for scalable wall functions
1320 */
1321 /*----------------------------------------------------------------------------*/
1322
1323 void
1324 cs_wall_functions_velocity(cs_wall_f_type_t iwallf,
1325 cs_lnum_t ifac,
1326 cs_real_t l_visc,
1327 cs_real_t t_visc,
1328 cs_real_t vel,
1329 cs_real_t y,
1330 cs_real_t rough_d,
1331 cs_real_t rnnb,
1332 cs_real_t kinetic_en,
1333 int *iuntur,
1334 cs_lnum_t *nsubla,
1335 cs_lnum_t *nlogla,
1336 cs_real_t *ustar,
1337 cs_real_t *uk,
1338 cs_real_t *yplus,
1339 cs_real_t *ypup,
1340 cs_real_t *cofimp,
1341 cs_real_t *dplus);
1342
1343 /*----------------------------------------------------------------------------*/
1344 /*!
1345 * \brief Compute the correction of the exchange coefficient between the
1346 * fluid and the wall for a turbulent flow.
1347 *
1348 * This is function of the dimensionless
1349 * distance to the wall \f$ y^+ = \dfrac{\centip \centf u_\star}{\nu}\f$.
1350 *
1351 * Then the return coefficient reads:
1352 * \f[
1353 * h_{tur} = Pr \dfrac{y^+}{T^+}
1354 * \f]
1355 *
1356 * \param[in] iwalfs type of wall functions for scalar
1357 * \param[in] l_visc kinematic viscosity
1358 * \param[in] prl laminar Prandtl number
1359 * \param[in] prt turbulent Prandtl number
1360 * \param[in] rough_t scalar roughness lenghth scale
1361 * \param[in] uk velocity scale based on TKE
1362 * \param[in] yplus dimensionless distance to the wall
1363 * \param[in] dplus dimensionless distance for scalable
1364 * wall functions
1365 * \param[out] htur corrected exchange coefficient
1366 * \param[out] yplim value of the limit for \f$ y^+ \f$
1367 */
1368 /*----------------------------------------------------------------------------*/
1369
1370 void
1371 cs_wall_functions_scalar(cs_wall_f_s_type_t iwalfs,
1372 cs_real_t l_visc,
1373 cs_real_t prl,
1374 cs_real_t prt,
1375 cs_real_t rough_t,
1376 cs_real_t uk,
1377 cs_real_t yplus,
1378 cs_real_t dplus,
1379 cs_real_t *htur,
1380 cs_real_t *yplim);
1381
1382 /*----------------------------------------------------------------------------*/
1383
1384 END_C_DECLS
1385
1386 #endif /* __CS_WALL_FUNCTIONS_H__ */
1387