7.1
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-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 
47 
48 /*=============================================================================
49  * Local Macro definitions
50  *============================================================================*/
51 
52 /*============================================================================
53  * Type definition
54  *============================================================================*/
55 
56 /* Wall function type */
57 /*--------------------*/
58 
59 typedef enum {
60 
70 
72 
73 typedef enum {
74 
80  CS_WALL_F_S_SMOOTH_ROUGH = 4,//TODO merge with MO ?
81 
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 
97 
98 /*============================================================================
99  * Global variables
100  *============================================================================*/
101 
102 /* Pointer to wall functions descriptor structure */
103 
105 
106 /*============================================================================
107  * Private function definitions
108  *============================================================================*/
109 
110 /*----------------------------------------------------------------------------*/
127 /*----------------------------------------------------------------------------*/
128 
129 inline static void
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 /*----------------------------------------------------------------------------*/
198 /*----------------------------------------------------------------------------*/
199 
200 inline static void
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
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
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 /*----------------------------------------------------------------------------*/
360 /*----------------------------------------------------------------------------*/
361 
362 inline static void
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 /*----------------------------------------------------------------------------*/
475 /*----------------------------------------------------------------------------*/
476 
477 inline static void
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 /*----------------------------------------------------------------------------*/
558 /*----------------------------------------------------------------------------*/
559 
560 inline static void
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
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 /*----------------------------------------------------------------------------*/
707 /*----------------------------------------------------------------------------*/
708 
709 inline static void
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 /*----------------------------------------------------------------------------*/
841 /*----------------------------------------------------------------------------*/
842 
843 inline static void
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 /*----------------------------------------------------------------------------*/
954 /*----------------------------------------------------------------------------*/
955 
956 inline static void
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 /*----------------------------------------------------------------------------*/
1029 /*----------------------------------------------------------------------------*/
1030 
1031 inline static void
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 /*----------------------------------------------------------------------------*/
1131 /*----------------------------------------------------------------------------*/
1132 
1133 inline static void
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 /*----------------------------------------------------------------------------*/
1196 /*----------------------------------------------------------------------------*/
1197 
1198 inline static void
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 
1292 
1293 /*----------------------------------------------------------------------------*/
1321 /*----------------------------------------------------------------------------*/
1322 
1323 void
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 /*----------------------------------------------------------------------------*/
1368 /*----------------------------------------------------------------------------*/
1369 
1370 void
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 
1385 
1386 #endif /* __CS_WALL_FUNCTIONS_H__ */
const double cs_turb_xkappa
Definition: cs_turbulence_model.c:387
int itytur
Definition: cs_turbulence_model.h:138
double precision epzero
epsilon
Definition: cstnum.f90:40
Definition: cs_wall_functions.h:75
cs_wall_f_type_t iwallf
Definition: cs_wall_functions.h:89
cs_turb_model_t * cs_get_glob_turb_model(void)
Provide write access to turbulence model structure.
Definition: cs_turbulence_model.c:1422
#define BEGIN_C_DECLS
Definition: cs_defs.h:510
integer(c_int), pointer, save iwalfs
Wall functions for scalar.
Definition: optcal.f90:576
const cs_real_t cs_math_epzero
Turbulence model general options descriptor.
Definition: cs_turbulence_model.h:114
#define CS_UNUSED(x)
Definition: cs_defs.h:496
double cs_turb_cmu025
Definition: cs_turbulence_model.c:448
cs_wall_functions_t * cs_get_glob_wall_functions(void)
Definition: cs_wall_functions.c:297
Definition: cs_wall_functions.h:76
Definition: cs_wall_functions.h:63
double precision, dimension(ncharm), save a2
Definition: cpincl.f90:233
double cs_turb_cmu
Definition: cs_turbulence_model.c:445
Definition: cs_wall_functions.h:79
Definition: cs_wall_functions.h:62
double cs_turb_dpow
Definition: cs_turbulence_model.c:438
double cs_real_t
Floating-point value.
Definition: cs_defs.h:322
const double cs_turb_bpow
Definition: cs_turbulence_model.c:435
Definition: cs_field_pointer.h:68
Definition: cs_wall_functions.h:64
void hturbp(const int *const iwalfs, const cs_real_t *const l_visc, const cs_real_t *const prl, const cs_real_t *const prt, const cs_real_t *const rough_t, const cs_real_t *const uk, const cs_real_t *const yplus, const cs_real_t *const dplus, cs_real_t *htur, cs_real_t *yplim)
Definition: cs_wall_functions.c:261
double cs_turb_crij2
Definition: cs_turbulence_model.c:485
cs_wall_f_s_type_t
Definition: cs_wall_functions.h:73
Definition: cs_wall_functions.h:66
const double cs_turb_vdriest
Definition: cs_turbulence_model.c:396
Definition: cs_wall_functions.h:61
Definition: cs_wall_functions.h:78
#define CS_MIN(a, b)
Definition: cs_defs.h:473
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:303
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:571
Definition: cs_wall_functions.h:65
double cs_turb_crij1
Definition: cs_turbulence_model.c:479
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:316
Definition: cs_wall_functions.h:69
#define END_C_DECLS
Definition: cs_defs.h:511
wall functions descriptor.
Definition: cs_wall_functions.h:87
#define _(String)
Definition: cs_defs.h:66
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.c:531
Definition: cs_wall_functions.h:80
Definition: cs_wall_functions.h:68
const double cs_turb_cstlog
Definition: cs_turbulence_model.c:407
#define CS_PROCF(x, y)
Definition: cs_defs.h:524
#define CS_MAX(a, b)
Definition: cs_defs.h:474
int bft_printf(const char *const format,...)
Replacement for printf() with modifiable behavior.
Definition: bft_printf.c:140
cs_wall_f_s_type_t iwalfs
Definition: cs_wall_functions.h:91
double precision, dimension(4, npot), save ka
Definition: ppcpfu.f90:159
Definition: cs_field_pointer.h:71
double ypluli
Definition: cs_wall_functions.h:93
cs_wall_f_type_t
Definition: cs_wall_functions.h:59
Definition: cs_wall_functions.h:77
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)
Compute the friction velocity and / .
Definition: cs_wall_functions.c:331
const cs_wall_functions_t * cs_glob_wall_functions
void wallfunctions(const int *const iwallf, const cs_lnum_t *const ifac, const cs_real_t *const viscosity, const cs_real_t *const t_visc, const cs_real_t *const vel, const cs_real_t *const y, const cs_real_t *const rough_d, const cs_real_t *const rnnb, const cs_real_t *const 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)
Definition: cs_wall_functions.c:213
Definition: cs_field_pointer.h:236
Definition: cs_wall_functions.h:67
const double cs_turb_apow
Definition: cs_turbulence_model.c:432