8.3
general documentation
cs_wall_functions.h
Go to the documentation of this file.
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-2024 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 "bft_printf.h"
35#include "cs_base.h"
36#include "cs_math.h"
37#include "cs_turbulence_model.h"
38
39/*----------------------------------------------------------------------------*/
40
42
43/*=============================================================================
44 * Local Macro definitions
45 *============================================================================*/
46
47/*============================================================================
48 * Type definition
49 *============================================================================*/
50
51/* Wall function type */
52/*--------------------*/
53
54typedef enum {
55
65
67
68typedef enum {
69
75 CS_WALL_F_S_SMOOTH_ROUGH = 4,//TODO merge with MO ?
76
78
79/* Wall functions descriptor */
80/*---------------------------*/
81
82typedef struct {
83
84 cs_wall_f_type_t iwallf; /* wall function type */
85
86 cs_wall_f_s_type_t iwalfs; /* wall function type for scalars */
87
88 double ypluli; /* limit value of y+ for the viscous
89 sublayer */
90
92
93/*============================================================================
94 * Global variables
95 *============================================================================*/
96
97/* Pointer to wall functions descriptor structure */
98
100
101/*============================================================================
102 * Private function definitions
103 *============================================================================*/
104
105/*----------------------------------------------------------------------------*/
122/*----------------------------------------------------------------------------*/
123
124inline static void
127 cs_real_t y,
128 int *iuntur,
129 cs_gnum_t *nsubla,
130 cs_gnum_t *nlogla,
131 cs_real_t *ustar,
132 cs_real_t *uk,
134 cs_real_t *ypup,
135 cs_real_t *cofimp)
136{
137 const double ypluli = cs_glob_wall_functions->ypluli;
138
139 const double ydvisc = y / l_visc;
140
141 /* Compute the friction velocity ustar */
142
143 *ustar = pow((vel/(cs_turb_apow * pow(ydvisc, cs_turb_bpow))), cs_turb_dpow);
144 *uk = *ustar;
145 *yplus = *ustar * ydvisc;
146
147 /* In the viscous sub-layer: U+ = y+ */
148 if (*yplus <= ypluli) {
149
150 *ustar = sqrt(vel / ydvisc);
151 *yplus = *ustar * ydvisc;
152 *uk = *ustar;
153 *ypup = 1.;
154 *cofimp = 0.;
155
156 /* Disable the wall funcion count the cell in the viscous sub-layer */
157 *iuntur = 0;
158 *nsubla += 1;
159
160 /* In the log layer */
161 } else {
162 *ypup = pow(vel, 2. * cs_turb_dpow-1.)
163 / pow(cs_turb_apow, 2. * cs_turb_dpow);
164 *cofimp = 1. + cs_turb_bpow
165 * pow(*ustar, cs_turb_bpow + 1. - 1./cs_turb_dpow)
166 * (pow(2., cs_turb_bpow - 1.) - 2.);
167
168 /* Count the cell in the log layer */
169 *nlogla += 1;
170
171 }
172}
173
174/*----------------------------------------------------------------------------*/
193/*----------------------------------------------------------------------------*/
194
195inline static void
198 cs_real_t y,
199 int *iuntur,
200 cs_gnum_t *nsubla,
201 cs_gnum_t *nlogla,
202 cs_real_t *ustar,
203 cs_real_t *uk,
205 cs_real_t *ypup,
206 cs_real_t *cofimp)
207{
208 const double ypluli = cs_glob_wall_functions->ypluli;
209
210 double ustarwer, ustarmin, ustaro, ydvisc;
211 double eps = 0.001;
212 int niter_max = 100;
213 int iter = 0;
214 double reynolds;
215
216 /* Compute the local Reynolds number */
217
218 ydvisc = y / l_visc;
219 reynolds = vel * ydvisc;
220
221 /*
222 * Compute the friction velocity ustar
223 */
224
225 /* In the viscous sub-layer: U+ = y+ */
226 if (reynolds <= ypluli * ypluli) {
227
228 *ustar = sqrt(vel / ydvisc);
229 *yplus = *ustar * ydvisc;
230 *uk = *ustar;
231 *ypup = 1.;
232 *cofimp = 0.;
233
234 /* Disable the wall funcion count the cell in the viscous sub-layer */
235 *iuntur = 0;
236 *nsubla += 1;
237
238 /* In the log layer */
239 } else {
240
241 /* The initial value is Wener or the minimun ustar to ensure convergence */
242 ustarwer = pow(fabs(vel) / cs_turb_apow / pow(ydvisc, cs_turb_bpow),
244 ustarmin = exp(-cs_turb_cstlog * cs_turb_xkappa)/ydvisc;
245 ustaro = cs_math_fmax(ustarwer, ustarmin);
246 *ustar = (cs_turb_xkappa * vel + ustaro)
247 / (log(ydvisc * ustaro) + cs_turb_xkappa * cs_turb_cstlog + 1.);
248
249 /* Iterative solving */
250 for (iter = 0; iter < niter_max
251 && fabs(*ustar - ustaro) >= eps * ustaro; iter++) {
252 ustaro = *ustar;
253 *ustar = (cs_turb_xkappa * vel + ustaro)
254 / (log(ydvisc * ustaro) + cs_turb_xkappa * cs_turb_cstlog + 1.);
255 }
256
257 if (iter >= niter_max) {
258 bft_printf(_("WARNING: non-convergence in the computation\n"
259 "******** of the friction velocity\n\n"
260 "friction vel: %f \n" ), *ustar);
261 }
262
263 *uk = *ustar;
264 *yplus = *ustar * ydvisc;
265 *ypup = *yplus / (log(*yplus) / cs_turb_xkappa + cs_turb_cstlog);
266 *cofimp = 1. - *ypup / cs_turb_xkappa * 1.5 / *yplus;
267
268 /* Count the cell in the log layer */
269 *nlogla += 1;
270
271 }
272
273}
274
275/*----------------------------------------------------------------------------
276 * Compute du+/dy+ for a given yk+.
277 *
278 * parameters:
279 * yplus <-- dimensionless distance
280 *
281 * returns:
282 * the resulting dimensionless velocity.
283 *----------------------------------------------------------------------------*/
284
285inline static cs_real_t
287 cs_real_t ka,
288 cs_real_t B,
289 cs_real_t cuv,
290 cs_real_t y0,
291 cs_real_t n)
292{
293 cs_real_t uplus, f_blend;
294
295 f_blend = exp(-0.25*cuv*pow(yp,3));
296 uplus = f_blend*yp + (log(yp)/ka +B)*(1.-exp(-pow(yp/y0,n)))*(1-f_blend);
297
298 return uplus;
299}
300
301/*----------------------------------------------------------------------------
302 * Compute du+/dy+ for a given yk+.
303 * parameters:
304 * yplus <-- dimensionless distance
305 * returns:
306 * the resulting dimensionless velocity gradient.
307 *----------------------------------------------------------------------------*/
308
309inline static cs_real_t
311 cs_real_t ka,
312 cs_real_t B,
313 cs_real_t cuv,
314 cs_real_t y0,
315 cs_real_t n)
316{
317 cs_real_t dupdyp;
318
319 dupdyp = exp(-0.25*cuv*pow(yp,3))
320 - 0.75*cuv*pow(yp,3.)*exp(-0.25*cuv*pow(yp,3.))
321 + n*(1.-exp(-0.25*cuv*pow(yp,3.)))*(pow(yp,n-1.)/pow(y0,n))
322 *exp(-pow(yp/y0,n))*((1./ka)*log(yp)+B)
323 + 0.75*cuv*pow(yp,2.)*exp(-0.25*cuv*pow(yp,3.))
324 *(1.-exp(-pow(yp/y0,n)))*((1./ka)*log(yp)+B)
325 + (1./ka/yp)*(1.-exp(-pow(yp/y0,n)))*(1-exp(-0.25*cuv*pow(yp,3.)));
326
327 return dupdyp;
328}
329
330/*----------------------------------------------------------------------------*/
353/*----------------------------------------------------------------------------*/
354
355inline static void
357 cs_real_t l_visc,
358 cs_real_t t_visc,
360 cs_real_t y,
361 cs_real_t kinetic_en,
362 int *iuntur,
363 cs_gnum_t *nsubla,
364 cs_gnum_t *nlogla,
365 cs_real_t *ustar,
366 cs_real_t *uk,
368 cs_real_t *ypup,
369 cs_real_t *cofimp)
370{
371 const double ypluli = cs_glob_wall_functions->ypluli;
372 double Re, g, t_visc_durb;
373 cs_real_t cstcuv, csty0, cstN;
374 cs_real_t dup1, dup2, uplus;
375
376 /* Local constants */
377 cstcuv = 1.0674e-3;
378 csty0 = 14.5e0;
379 cstN = 2.25e0;
380
381 /* Iterative process to determine uk through TKE law */
382 Re = sqrt(kinetic_en) * y / l_visc;
383 g = exp(-Re/11.);
384
385 /* Comutation of uk*/
386 *uk = sqrt( (1.-g) * sqrt(cs_turb_cmu) * kinetic_en
387 + g * l_visc * vel/y);
388
389 /* Local value of y+, estimated U+ */
390 *yplus = *uk * y / l_visc;
391 uplus = _uplus( *yplus, cs_turb_xkappa, cs_turb_cstlog, cstcuv, csty0, cstN);
392 /* Deduced velocity sclale uet*/
393 *ustar = vel / uplus;
394
395 if (*yplus < 1.e-1) {
396
397 *ypup = 1.0;
398 *cofimp = 0.0;
399
400 *iuntur = 0;
401 *nsubla += 1;
402
403 }
404 else {
405
406 /* Dimensionless velocity gradient in y+ */
408 cstcuv, csty0, cstN);
409 /* Dimensionless velocity gradient in 2 x y+ */
410 dup2 = _dupdyp(2.0 * *yplus, cs_turb_xkappa,
411 cs_turb_cstlog, cstcuv, csty0, cstN);
412
413 *ypup = *yplus / uplus;
414
415 /* ------------------------------------------------------------
416 * Cofimp = U,F/U,I is built so that the theoretical expression
417 * of the production P_theo = dup1 * (1.0 - dup1) is equal to
418 * P_calc = mu_t,I * ((U,I - U,F + IF*dup2)/(2IF) )^2
419 * This is a generalization of the process implemented in the 2
420 * scales wall function (iwallf = 3).
421 * ------------------------------------------------------------*/
422
423 /* Turbulent viscocity is modified for RSM so that its expression
424 * remain valid down to the wall, according to Durbin :
425 * nu_t = 0.22 * v'2 * k / eps */
426 const cs_turb_model_t *turb_model = cs_get_glob_turb_model();
427 assert(turb_model != NULL);
428 if (turb_model->itytur == 3)
429 t_visc_durb = t_visc / (kinetic_en * cs_turb_cmu ) * rnnb * 0.22;
430 else
431 t_visc_durb = t_visc;
432
433 *cofimp
434 = 1. - *ypup * (2. * sqrt(l_visc / t_visc_durb * dup1 * (1. - dup1))
435 - dup2);
436
437 /* log layer */
438 if (*yplus > ypluli) {
439 *nlogla += 1;
440 /* viscous sub-layer or buffer layer*/
441 } else {
442 *iuntur = 0;
443 *nsubla += 1;
444 }
445 }
446}
447
448/*----------------------------------------------------------------------------*/
468/*----------------------------------------------------------------------------*/
469
470inline static void
472 cs_real_t t_visc,
474 cs_real_t y,
475 cs_real_t kinetic_en,
476 int *iuntur,
477 cs_gnum_t *nsubla,
478 cs_gnum_t *nlogla,
479 cs_real_t *ustar,
480 cs_real_t *uk,
482 cs_real_t *ypup,
483 cs_real_t *cofimp)
484{
485 const double ypluli = cs_glob_wall_functions->ypluli;
486
487 double rcprod, ml_visc, Re, g;
488
489 /* Compute the friction velocity ustar */
490
491 /* Blending for very low values of k */
492 Re = sqrt(kinetic_en) * y / l_visc;
493 g = exp(-Re/11.);
494
495 *uk = sqrt( (1.-g) * sqrt(cs_turb_cmu) * kinetic_en
496 + g * l_visc * vel / y);
497
498 *yplus = *uk * y / l_visc;
499
500 /* log layer */
501 if (*yplus > ypluli) {
502
503 *ustar = vel / (log(*yplus) / cs_turb_xkappa + cs_turb_cstlog);
504 *ypup = *yplus / (log(*yplus) / cs_turb_xkappa + cs_turb_cstlog);
505 /* Mixing length viscosity */
506 ml_visc = cs_turb_xkappa * l_visc * *yplus;
508 cs_math_fmax(1., sqrt(ml_visc / t_visc)) / *yplus);
509 *cofimp = 1. - *ypup / cs_turb_xkappa * ( 2. * rcprod - 1. / (2. * *yplus));
510
511 *nlogla += 1;
512
513 /* viscous sub-layer */
514 }
515 else {
516
517 if (*yplus > 1.e-12) {
518 *ustar = fabs(vel / *yplus); /* FIXME remove that: its is here only to
519 be fully equivalent to the former code. */
520 }
521 else {
522 *ustar = 0.;
523 }
524 *ypup = 1.;
525 *cofimp = 0.;
526
527 *iuntur = 0;
528 *nsubla += 1;
529
530 }
531}
532
533/*----------------------------------------------------------------------------*/
554/*----------------------------------------------------------------------------*/
555
556inline static void
558 cs_real_t t_visc,
560 cs_real_t y,
561 cs_real_t kinetic_en,
562 int *iuntur,
563 cs_gnum_t *nsubla,
564 cs_gnum_t *nlogla,
565 cs_real_t *ustar,
566 cs_real_t *uk,
568 cs_real_t *dplus,
569 cs_real_t *ypup,
570 cs_real_t *cofimp)
571{
572 CS_UNUSED(iuntur);
573
574 const double ypluli = cs_glob_wall_functions->ypluli;
575
576 double rcprod, ml_visc, Re, g;
577 /* Compute the friction velocity ustar */
578
579 /* Blending for very low values of k */
580 Re = sqrt(kinetic_en) * y / l_visc;
581 g = exp(-Re/11.);
582
583 *uk = sqrt( (1.-g) * sqrt(cs_turb_cmu) * kinetic_en
584 + g * l_visc * vel / y);
585
586 *yplus = *uk * y / l_visc;
587
588 /* Compute the friction velocity ustar */
589 *uk = sqrt(sqrt(cs_turb_cmu) *kinetic_en);//FIXME
590 *yplus = *uk * y / l_visc;//FIXME
591
592 /* Log layer */
593 if (*yplus > ypluli) {
594
595 *dplus = 0.;
596
597 *nlogla += 1;
598
599 /* Viscous sub-layer and therefore shift */
600 } else {
601
602 *dplus = ypluli - *yplus;
603
604 /* Count the cell as if it was in the viscous sub-layer */
605 *nsubla += 1;
606
607 }
608
609 /* Mixing length viscosity */
610 ml_visc = cs_turb_xkappa * l_visc * (*yplus + *dplus);
612 cs_math_fmax(1., sqrt(ml_visc / t_visc))
613 / (*yplus + *dplus));
614
615 *ustar = vel / (log(*yplus + *dplus) / cs_turb_xkappa + cs_turb_cstlog);
616 *ypup = *yplus / (log(*yplus + *dplus) / cs_turb_xkappa + cs_turb_cstlog);
617 *cofimp = 1. - *ypup
618 / cs_turb_xkappa * (2. * rcprod - 1. / (2. * *yplus + *dplus));
619}
620
621/*----------------------------------------------------------------------------
622 * Compute u+ for a given yk+ between 0.1 and 200 according to the two
623 * scales wall functions using Van Driest mixing length.
624 * This function holds the coefficients of the polynome fitting log(u+).
625 *
626 * parameters:
627 * yplus <-- dimensionless distance
628 *
629 * returns:
630 * the resulting dimensionless velocity.
631 *----------------------------------------------------------------------------*/
632
633inline static cs_real_t
635{
636 /* Coefficients of the polynome fitting log(u+) for yk < 200 */
637 static double aa[11] = {-0.0091921, 3.9577, 0.031578,
638 -0.51013, -2.3254, -0.72665,
639 2.969, 0.48506, -1.5944,
640 0.087309, 0.1987 };
641
642 cs_real_t y1,y2,y3,y4,y5,y6,y7,y8,y9,y10, uplus;
643
644 y1 = 0.25 * log(yplus);
645 y2 = y1 * y1;
646 y3 = y2 * y1;
647 y4 = y3 * y1;
648 y5 = y4 * y1;
649 y6 = y5 * y1;
650 y7 = y6 * y1;
651 y8 = y7 * y1;
652 y9 = y8 * y1;
653 y10 = y9 * y1;
654
655 uplus = aa[0]
656 + aa[1] * y1
657 + aa[2] * y2
658 + aa[3] * y3
659 + aa[4] * y4
660 + aa[5] * y5
661 + aa[6] * y6
662 + aa[7] * y7
663 + aa[8] * y8
664 + aa[9] * y9
665 + aa[10] * y10;
666
667 return exp(uplus);
668}
669
670/*----------------------------------------------------------------------------*/
705/*----------------------------------------------------------------------------*/
706
707inline static void
709 cs_real_t l_visc,
711 cs_real_t y,
712 cs_real_t kinetic_en,
713 int *iuntur,
714 cs_gnum_t *nsubla,
715 cs_gnum_t *nlogla,
716 cs_real_t *ustar,
717 cs_real_t *uk,
719 cs_real_t *ypup,
720 cs_real_t *cofimp,
721 cs_real_t *lmk,
722 cs_real_t kr,
723 bool wf)
724{
725 double urplus, d_up, lmk15;
726
727 if (wf)
728 *uk = sqrt(sqrt((1.-cs_turb_crij2)/cs_turb_crij1 * rnnb * kinetic_en));
729
730 /* Set a low threshold value in case tangential velocity is zero */
731 *yplus = cs_math_fmax(*uk * y / l_visc, 1.e-4);
732
733 /* Dimensionless roughness */
734 cs_real_t krp = *uk * kr / l_visc;
735
736 /* Extension of Van Driest mixing length according to Rotta (1962) with
737 Cebeci & Chang (1978) correlation */
738 cs_real_t dyrp = 0.9 * (sqrt(krp) - krp * exp(-krp / 6.));
739 cs_real_t yrplus = *yplus + dyrp;
740
741 if (dyrp <= 1.e-1)
742 d_up = dyrp;
743 else if (dyrp <= 200.)
744 d_up = _vdriest_dupdyp_integral(dyrp);
745 else
746 d_up = 16.088739022054590 + log(dyrp/200.) / cs_turb_xkappa;
747
748 if (yrplus <= 1.e-1) {
749
750 urplus = yrplus;
751
752 if (wf) {
753 *iuntur = 0;
754 *nsubla += 1;
755
756 *lmk = 0.;
757
758 *ypup = 1.;
759
760 *cofimp = 0.;
761 }
762
763 } else if (yrplus <= 200.) {
764
765 urplus = _vdriest_dupdyp_integral(yrplus);
766
767 if (wf) {
768 *nlogla += 1;
769
770 *ypup = *yplus / (urplus-d_up);
771
772 /* Mixing length in y+ */
773 *lmk = cs_turb_xkappa * (*yplus) *(1-exp(- (*yplus) / cs_turb_vdriest));
774
775 /* Mixing length in 3/2*y+ */
776 lmk15 = cs_turb_xkappa * 1.5 * (*yplus) *(1-exp(- 1.5 * (*yplus)
777 / cs_turb_vdriest));
778
779 *cofimp = 1. - (2. / (1. + *lmk) - 1. / (1. + lmk15)) * *ypup;
780 }
781
782 } else {
783
784 urplus = 16.088739022054590 + log(yrplus/200) / cs_turb_xkappa;
785
786 if (wf) {
787 *nlogla += 1;
788
789 *ypup = *yplus / (urplus-d_up);
790
791 /* Mixing length in y+ */
792 *lmk = cs_turb_xkappa * (*yplus) *(1-exp(- (*yplus) / cs_turb_vdriest));
793
794 /* Mixing length in 3/2*y+ */
795 lmk15 = cs_turb_xkappa * 1.5 * (*yplus) *(1-exp(- 1.5 * (*yplus)
796 / cs_turb_vdriest));
797
798 *cofimp = 1. - (2. / *lmk - 1. / lmk15) * *ypup;
799 }
800
801 }
802
803 *ustar = vel / (urplus-d_up);
804}
805
806/*----------------------------------------------------------------------------*/
839/*----------------------------------------------------------------------------*/
840
841inline static void
843 cs_real_t t_visc,
845 cs_real_t y,
846 cs_real_t rough_d,
847 cs_real_t kinetic_en,
848 int *iuntur,
849 cs_gnum_t *nsubla,
850 cs_gnum_t *nlogla,
851 cs_real_t *ustar,
852 cs_real_t *uk,
854 cs_real_t *dplus,
855 cs_real_t *ypup,
856 cs_real_t *cofimp)
857{
858 CS_UNUSED(iuntur);
859
860 const double ypluli = cs_glob_wall_functions->ypluli;
861
862 double rcprod, ml_visc, Re, g;
863
864 /* Compute the friction velocity ustar */
865
866 /* Shifting of the wall distance to be consistant with
867 * the fully rough wall function
868 *
869 * ln((y+y0)/y0) = ln((y+y0)/alpha xi) + kappa * 5.2
870 *
871 * y0 = xi * exp(-kappa * 8.5)
872 * where xi is the sand grain roughness here
873 * y0 = alpha * xi * exp(-kappa * 5.2)
874 *
875 * so:
876 * alpha = exp(-kappa * (8.5 - 5.2)) = 0.25
877 *
878 */
879 cs_real_t y0 = rough_d;
880 /* Note : Sand grain roughness given by:
881 cs_real_t sg_rough = rough_d * exp(cs_turb_xkappa*cs_turb_cstlog_rough);
882 */
883
884 /* Blending for very low values of k */
885 Re = sqrt(kinetic_en) * (y + y0) / l_visc;
886 g = exp(-Re/11.);
887
888 *uk = sqrt( (1.-g) * sqrt(cs_turb_cmu) * kinetic_en
889 + g * l_visc * vel / (y + y0));
890
891
892 *yplus = *uk * y / l_visc;
893
894 /* As for scalable wall functions, yplus is shifted of "dplus" */
895 *dplus = *uk * y0 / l_visc;
896
897 /* Shift of the velocity profile due to roughness */
898 cs_real_t shift_vel = -log(1. + y0 * exp(cs_turb_xkappa * cs_turb_cstlog) * *uk/l_visc)
900
901 /* Log layer and shifted with the roughness */
902 if (*yplus > ypluli) {
903
904 *nlogla += 1;
905
906 }
907
908 /* Viscous sub-layer and therefore shift again */
909 else {
910
911 *dplus = ypluli - *yplus;
912 /* Count the cell as if it was in the viscous sub-layer */
913 *nsubla += 1;
914
915 }
916
917 cs_real_t uplus = log(*yplus + *dplus)
918 / cs_turb_xkappa + cs_turb_cstlog + shift_vel;
919 *ustar = vel / uplus;
920#if 0
921 bft_printf("uet=%f, u=%f, uplus=%f, yk=%f, duplus=%f\n",
922 *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);
929 cs_math_fmax(1., sqrt(ml_visc / t_visc))
930 / (*yplus + *dplus));
931 *cofimp = 1. - *yplus / (cs_turb_xkappa * uplus)
932 * ( 2. * rcprod - 1. / (2. * *yplus + *dplus));
933
934}
935
936/*----------------------------------------------------------------------------*/
956/*----------------------------------------------------------------------------*/
957
958inline static void
960 cs_real_t t_visc,
962 cs_real_t y,
963 int *iuntur,
964 cs_gnum_t *nsubla,
965 cs_gnum_t *nlogla,
966 cs_real_t *ustar,
967 cs_real_t *uk,
969 cs_real_t *dplus,
970 cs_real_t *ypup,
971 cs_real_t *cofimp)
972{
973 CS_UNUSED(t_visc);
974 CS_UNUSED(nlogla);
975 CS_UNUSED(dplus);
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 /* Count the cell as if it was in the viscous sub-layer */
987 *nsubla += 1;
988}
989
990/*----------------------------------------------------------------------------*/
1020/*----------------------------------------------------------------------------*/
1021
1022inline static void
1024 cs_real_t prl,
1025 cs_real_t prt,
1026 cs_real_t rough_t,
1027 cs_real_t uk,
1029 cs_real_t dplus,
1030 cs_real_t *htur,
1031 cs_real_t *yplim)
1032{
1033 /* Local variables */
1034 double tplus;
1035 double beta2,a2;
1036 double yp2;
1037 double prlm1;
1038
1039 const double epzero = cs_math_epzero;
1040
1041 /*==========================================================================
1042 1. Initializations
1043 ==========================================================================*/
1044
1045 (*htur) = cs_math_fmax(yplus, epzero) / cs_math_fmax(yplus+dplus, epzero);
1046
1047 prlm1 = 0.1;
1048
1049 /* Sand grain roughness is:
1050 * zeta = z0 * exp(kappa 8.5)
1051 * Then:
1052 * hp = zeta uk / nu * exp( -kappa(8.5 - 5.2))
1053 * = z0 * uk / nu * exp(kappa * 5.2)
1054 * where 5.2 is the smooth log constant, and 8.5 the rough one
1055 *
1056 * FIXME check if we should use a molecular Schmidt number
1057 */
1058 cs_real_t hp = rough_t *uk / l_visc * exp(cs_turb_xkappa*cs_turb_cstlog);
1059
1060 /* Shift of the temperature profile due to roughness */
1061 cs_real_t shift_temp = -log(1. + hp);
1062
1063 /*==========================================================================
1064 2. Compute htur for small Prandtl numbers
1065 ==========================================================================*/
1066
1067 if (prl <= prlm1) {
1068 (*yplim) = prt/(prl*cs_turb_xkappa);
1069 if (yplus > (*yplim)) {
1070 tplus = prl*(*yplim) + prt/cs_turb_xkappa
1071 * (log((yplus+dplus)/(*yplim)) + shift_temp);
1072 (*htur) = prl * yplus / tplus;
1073 }
1074
1075 /*========================================================================
1076 3. Compute htur for the model with three sub-layers
1077 ========================================================================*/
1078
1079 } else {
1080 yp2 = sqrt(cs_turb_xkappa*1000./prt);
1081 (*yplim) = pow(1000./prl, 1./3.);
1082
1083 a2 = 15.*pow(prl, 2./3.);
1084
1085 if (yplus >= (*yplim) && yplus < yp2) {
1086 tplus = a2 - 500./((yplus+dplus)*(yplus+dplus));
1087 (*htur) = prl * yplus / tplus;
1088 }
1089
1090 if (yplus >= yp2) {
1091 beta2 = a2 - 0.5 * prt /cs_turb_xkappa;
1092 tplus = beta2 + prt/cs_turb_xkappa*log((yplus+dplus)/yp2);
1093 (*htur) = prl * yplus / tplus;
1094 }
1095
1096 }
1097}
1098
1099/*----------------------------------------------------------------------------*/
1123/*----------------------------------------------------------------------------*/
1124
1125inline static void
1127 cs_real_t prt,
1129 cs_real_t *htur)
1130{
1131 cs_real_t prlrat = prl / prt;
1132
1133 /* Parameters of the numerical quadrature */
1134 const int ninter_max = 100;
1135 const cs_real_t ypmax = 1.e2;
1136
1137 /* No correction for very small yplus */
1138 if (yplus <= 0.1)
1139 *htur = 1.;
1140 else {
1141 cs_real_t ypint = cs_math_fmin(yplus, ypmax);
1142
1143 /* The number of sub-intervals is taken proportional to yplus and equal to
1144 * ninter_max if yplus=ypmax */
1145
1146 int npeff = cs_math_fmax((int)(ypint / ypmax * (double)(ninter_max)), 1);
1147
1148 double dy = ypint / (double)(npeff);
1149 cs_real_t stplus = 0.;
1150 cs_real_t nut1 = 0.;
1151 cs_real_t nut2 = 0.;
1152
1153 for (int ip = 1; ip <= npeff; ip++) {
1154 double yp = ypint * (double)(ip) / (double)(npeff);
1155 nut2 = cs_turb_xkappa * yp * (1. - exp(-yp / cs_turb_vdriest));
1156 stplus += dy / (1. + prlrat * 0.5 * (nut1 + nut2));
1157 nut1 = nut2;
1158 }
1159
1160 if (yplus > ypint) {
1161 cs_real_t r = prlrat * cs_turb_xkappa;
1162 stplus += log( (1. + r*yplus) / (1. + r*ypint)) / r;
1163 }
1164
1165 if (stplus >= 1.e-6)
1166 *htur = yplus / stplus;
1167 else
1168 *htur = 1.;
1169 }
1170}
1171
1172/*----------------------------------------------------------------------------*/
1188/*----------------------------------------------------------------------------*/
1189
1190inline static void
1192 cs_real_t prl,
1193 cs_real_t prt,
1194 cs_real_t rough_t,
1195 cs_real_t uk,
1197 cs_real_t dplus,
1198 cs_real_t *htur)
1199{
1200 CS_UNUSED(prt);
1201
1202 /* Sand grain roughness is:
1203 * zeta = z0 * exp(kappa 8.5)
1204 * Then:
1205 * hp = zeta uk / nu * exp( -kappa(8.5 - 5.2))
1206 * = z0 * uk / nu * exp(kappa * 5.2)
1207 * where 5.2 is the smooth log constant, and 8.5 the rough one
1208 *
1209 * FIXME check if we should use a molecular Schmidt number
1210 */
1211 cs_real_t hp = rough_t *uk / l_visc * exp(cs_turb_xkappa*cs_turb_cstlog);
1212 const double ypluli = cs_glob_wall_functions->ypluli;
1213 const double epzero = cs_math_epzero;
1214
1216
1217 /* Shift of the temperature profile due to roughness */
1218 cs_real_t shift_temp = -log(1. + hp);
1219
1220 if (yplus > ypluli) {
1221 cs_real_t tplus = prt * ( (log(yplus+dplus) + shift_temp)/cs_turb_xkappa
1222 + cs_turb_cstlog);
1223 (*htur) = prl * yplus / tplus;
1224 }
1225}
1226
1227/*============================================================================
1228 * Public function definitions for Fortran API
1229 *============================================================================*/
1230
1231/*=============================================================================
1232 * Public function prototypes
1233 *============================================================================*/
1234
1235/*----------------------------------------------------------------------------
1236 *! \brief Provide access to cs_glob_wall_functions
1237 *----------------------------------------------------------------------------*/
1238
1241
1242/*----------------------------------------------------------------------------*/
1269/*----------------------------------------------------------------------------*/
1270
1271void
1273 cs_real_t l_visc,
1274 cs_real_t t_visc,
1275 cs_real_t vel,
1276 cs_real_t y,
1277 cs_real_t rough_d,
1278 cs_real_t rnnb,
1279 cs_real_t kinetic_en,
1280 int *iuntur,
1281 cs_gnum_t *nsubla,
1282 cs_gnum_t *nlogla,
1283 cs_real_t *ustar,
1284 cs_real_t *uk,
1286 cs_real_t *ypup,
1287 cs_real_t *cofimp,
1288 cs_real_t *dplus);
1289
1290/*----------------------------------------------------------------------------*/
1315/*----------------------------------------------------------------------------*/
1316
1317void
1319 cs_real_t l_visc,
1320 cs_real_t prl,
1321 cs_real_t prt,
1322 cs_real_t rough_t,
1323 cs_real_t uk,
1325 cs_real_t dplus,
1326 cs_real_t *htur,
1327 cs_real_t *yplim);
1328
1329/*----------------------------------------------------------------------------*/
1337/*----------------------------------------------------------------------------*/
1338
1339void
1341 cs_real_t *st_exp,
1342 cs_real_t *st_imp);
1343
1344/*----------------------------------------------------------------------------*/
1345
1347
1348#endif /* __CS_WALL_FUNCTIONS_H__ */
int bft_printf(const char *const format,...)
Replacement for printf() with modifiable behavior.
Definition: bft_printf.cpp:140
#define BEGIN_C_DECLS
Definition: cs_defs.h:542
double cs_real_t
Floating-point value.
Definition: cs_defs.h:342
#define _(String)
Definition: cs_defs.h:67
uint64_t cs_gnum_t
global mesh entity number
Definition: cs_defs.h:325
#define CS_UNUSED(x)
Definition: cs_defs.h:528
#define END_C_DECLS
Definition: cs_defs.h:543
@ eps
Definition: cs_field_pointer.h:73
@ vel
Definition: cs_field_pointer.h:70
@ yplus
Definition: cs_field_pointer.h:242
const cs_real_t cs_math_epzero
static CS_F_HOST_DEVICE cs_real_t cs_math_fmin(cs_real_t x, cs_real_t y)
Compute the min value of two real values.
Definition: cs_math.h:484
static CS_F_HOST_DEVICE cs_real_t cs_math_fmax(cs_real_t x, cs_real_t y)
Compute the max value of two real values.
Definition: cs_math.h:503
double cs_turb_vdriest
Definition: cs_turbulence_model.cpp:446
double cs_turb_crij2
Definition: cs_turbulence_model.cpp:551
double cs_turb_crij1
Definition: cs_turbulence_model.cpp:545
cs_turb_model_t * cs_get_glob_turb_model(void)
Provide write access to turbulence model structure.
Definition: cs_turbulence_model.cpp:1425
double cs_turb_cstlog
Definition: cs_turbulence_model.cpp:458
double cs_turb_dpow
Definition: cs_turbulence_model.cpp:491
double cs_turb_apow
Definition: cs_turbulence_model.cpp:485
double cs_turb_cmu
Definition: cs_turbulence_model.cpp:499
double cs_turb_bpow
Definition: cs_turbulence_model.cpp:488
double cs_turb_xkappa
Definition: cs_turbulence_model.cpp:437
cs_wall_f_s_type_t
Definition: cs_wall_functions.h:68
@ CS_WALL_F_S_SMOOTH_ROUGH
Definition: cs_wall_functions.h:75
@ CS_WALL_F_S_MONIN_OBUKHOV
Definition: cs_wall_functions.h:74
@ CS_WALL_F_S_VDRIEST
Definition: cs_wall_functions.h:72
@ CS_WALL_F_S_ARPACI_LARSEN
Definition: cs_wall_functions.h:71
@ CS_WALL_F_S_UNSET
Definition: cs_wall_functions.h:70
@ CS_WALL_F_S_LOUIS
Definition: cs_wall_functions.h:73
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)
Definition: cs_wall_functions.h:310
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)
Compute the correction of the exchange coefficient between the fluid and the wall for a turbulent flo...
Definition: cs_wall_functions.cpp:438
void cs_immersed_boundary_wall_functions(int f_id, cs_real_t *st_exp, cs_real_t *st_imp)
Compute boundary contributions for all immersed boundaries.
Definition: cs_wall_functions.cpp:509
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_gnum_t *nsubla, cs_gnum_t *nlogla, cs_real_t *ustar, cs_real_t *uk, cs_real_t *yplus, cs_real_t *ypup, cs_real_t *cofimp)
Continuous law of the wall between the linear and log law, with two velocity scales based on the fric...
Definition: cs_wall_functions.h:356
static void cs_wall_functions_1scale_log(cs_real_t l_visc, cs_real_t vel, cs_real_t y, int *iuntur, cs_gnum_t *nsubla, cs_gnum_t *nlogla, cs_real_t *ustar, cs_real_t *uk, cs_real_t *yplus, cs_real_t *ypup, cs_real_t *cofimp)
Log law: piecewise linear and log, with one velocity scale based on the friction.
Definition: cs_wall_functions.h:196
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_gnum_t *nsubla, cs_gnum_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)
Two velocity scales wall function with automatic switch from rough to smooth.
Definition: cs_wall_functions.h:842
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)
Rough Smooth Thermal Wall Function - Prototype.
Definition: cs_wall_functions.h:1191
static void cs_wall_functions_1scale_power(cs_real_t l_visc, cs_real_t vel, cs_real_t y, int *iuntur, cs_gnum_t *nsubla, cs_gnum_t *nlogla, cs_real_t *ustar, cs_real_t *uk, cs_real_t *yplus, cs_real_t *ypup, cs_real_t *cofimp)
Power law: Werner & Wengle.
Definition: cs_wall_functions.h:125
const cs_wall_functions_t * cs_glob_wall_functions
static cs_real_t _vdriest_dupdyp_integral(cs_real_t yplus)
Definition: cs_wall_functions.h:634
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)
The correction of the exchange coefficient is computed thanks to a similarity model between dynamic v...
Definition: cs_wall_functions.h:1023
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_gnum_t *nsubla, cs_gnum_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)
No wall function.
Definition: cs_wall_functions.h:959
cs_wall_functions_t * cs_get_glob_wall_functions(void)
Definition: cs_wall_functions.cpp:218
void cs_wall_functions_velocity(cs_wall_f_type_t iwallf, 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_gnum_t *nsubla, cs_gnum_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)
Compute the friction velocity and / .
Definition: cs_wall_functions.cpp:251
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_gnum_t *nsubla, cs_gnum_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)
Two velocity scales wall function using Van Driest mixing length.
Definition: cs_wall_functions.h:708
static void cs_wall_functions_s_vdriest(cs_real_t prl, cs_real_t prt, cs_real_t yplus, cs_real_t *htur)
The correction of the exchange coefficient is computed thanks to a numerical integration of:
Definition: cs_wall_functions.h:1126
cs_wall_f_type_t
Definition: cs_wall_functions.h:54
@ CS_WALL_F_1SCALE_LOG
Definition: cs_wall_functions.h:59
@ CS_WALL_F_1SCALE_POWER
Definition: cs_wall_functions.h:58
@ CS_WALL_F_2SCALES_SMOOTH_ROUGH
Definition: cs_wall_functions.h:63
@ CS_WALL_F_2SCALES_LOG
Definition: cs_wall_functions.h:60
@ CS_WALL_F_2SCALES_VDRIEST
Definition: cs_wall_functions.h:62
@ CS_WALL_F_DISABLED
Definition: cs_wall_functions.h:57
@ CS_WALL_F_UNSET
Definition: cs_wall_functions.h:56
@ CS_WALL_F_2SCALES_CONTINUOUS
Definition: cs_wall_functions.h:64
@ CS_WALL_F_SCALABLE_2SCALES_LOG
Definition: cs_wall_functions.h:61
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_gnum_t *nsubla, cs_gnum_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)
Scalable wall function: shift the wall if .
Definition: cs_wall_functions.h:557
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_gnum_t *nsubla, cs_gnum_t *nlogla, cs_real_t *ustar, cs_real_t *uk, cs_real_t *yplus, cs_real_t *ypup, cs_real_t *cofimp)
Log law: piecewise linear and log, with two velocity scales based on the friction and the turbulent k...
Definition: cs_wall_functions.h:471
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)
Definition: cs_wall_functions.h:286
double precision epzero
epsilon
Definition: cstnum.f90:40
real(c_double), pointer, save ypluli
limit value of for the viscous sublayer. ypluli depends on the chosen wall function: it is initializ...
Definition: cstphy.f90:313
integer(c_int), pointer, save iwallf
Wall functions Indicates the type of wall function used for the velocity boundary conditions on a fri...
Definition: optcal.f90:340
integer(c_int), pointer, save iwalfs
Wall functions for scalar.
Definition: optcal.f90:345
real(c_double), dimension(:), pointer, save a2
Definition: cpincl.f90:109
Turbulence model general options descriptor.
Definition: cs_turbulence_model.h:127
int itytur
Definition: cs_turbulence_model.h:154
wall functions descriptor.
Definition: cs_wall_functions.h:82
double ypluli
Definition: cs_wall_functions.h:88
cs_wall_f_s_type_t iwalfs
Definition: cs_wall_functions.h:86
cs_wall_f_type_t iwallf
Definition: cs_wall_functions.h:84