8.2
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 
54 typedef enum {
55 
65 
67 
68 typedef enum {
69 
75  CS_WALL_F_S_SMOOTH_ROUGH = 4,//TODO merge with MO ?
76 
78 
79 /* Wall functions descriptor */
80 /*---------------------------*/
81 
82 typedef 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 
124 inline static void
126  cs_real_t vel,
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,
133  cs_real_t *yplus,
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 
195 inline static void
197  cs_real_t vel,
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,
204  cs_real_t *yplus,
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),
243  cs_turb_dpow);
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 
285 inline 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 
309 inline 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 
355 inline static void
357  cs_real_t l_visc,
358  cs_real_t t_visc,
359  cs_real_t vel,
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,
367  cs_real_t *yplus,
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) * cs_turb_cmu025 * cs_turb_cmu025 * 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 
470 inline static void
472  cs_real_t t_visc,
473  cs_real_t vel,
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,
481  cs_real_t *yplus,
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) * cs_turb_cmu025 * cs_turb_cmu025 * 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;
507  rcprod = cs_math_fmin(cs_turb_xkappa,
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 
556 inline static void
558  cs_real_t t_visc,
559  cs_real_t vel,
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,
567  cs_real_t *yplus,
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) * cs_turb_cmu025 * cs_turb_cmu025 * kinetic_en
584  + g * l_visc * vel / y);
585 
586  *yplus = *uk * y / l_visc;
587 
588  /* Compute the friction velocity ustar */
589  *uk = cs_turb_cmu025 * sqrt(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);
611  rcprod = cs_math_fmin(cs_turb_xkappa,
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 
633 inline 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 
707 inline static void
709  cs_real_t l_visc,
710  cs_real_t vel,
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,
718  cs_real_t *yplus,
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 
841 inline static void
843  cs_real_t t_visc,
844  cs_real_t vel,
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,
853  cs_real_t *yplus,
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) * cs_turb_cmu025 * cs_turb_cmu025 * 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)
899  / cs_turb_xkappa;
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);
928  rcprod = cs_math_fmin(cs_turb_xkappa,
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 
958 inline static void
960  cs_real_t t_visc,
961  cs_real_t vel,
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,
968  cs_real_t *yplus,
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 
1022 inline static void
1024  cs_real_t prl,
1025  cs_real_t prt,
1026  cs_real_t rough_t,
1027  cs_real_t uk,
1028  cs_real_t yplus,
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 
1125 inline static void
1127  cs_real_t prt,
1128  cs_real_t yplus,
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 
1190 inline static void
1192  cs_real_t prl,
1193  cs_real_t prt,
1194  cs_real_t rough_t,
1195  cs_real_t uk,
1196  cs_real_t yplus,
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 
1215  (*htur) = cs_math_fmax(yplus,epzero)/cs_math_fmax(yplus+dplus,epzero);
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  * Wrapper to cs_wall_functions_velocity.
1233  *----------------------------------------------------------------------------*/
1234 
1235 void CS_PROCF (wallfunctions, WALLFUNCTIONS)
1236 (
1237  const int *const iwallf,
1238  const cs_lnum_t *const ifac,
1239  const cs_real_t *const viscosity,
1240  const cs_real_t *const t_visc,
1241  const cs_real_t *const vel,
1242  const cs_real_t *const y,
1243  const cs_real_t *const rough_d,
1244  const cs_real_t *const rnnb,
1245  const cs_real_t *const kinetic_en,
1246  int *iuntur,
1247  cs_lnum_t *nsubla,
1248  cs_lnum_t *nlogla,
1249  cs_real_t *ustar,
1250  cs_real_t *uk,
1251  cs_real_t *yplus,
1252  cs_real_t *ypup,
1253  cs_real_t *cofimp,
1254  cs_real_t *dplus
1255 );
1256 
1257 /*----------------------------------------------------------------------------
1258  * Wrapper to cs_wall_functions_scalar.
1259  *----------------------------------------------------------------------------*/
1260 
1261 void CS_PROCF (hturbp, HTURBP)
1262 (
1263  const int *const iwalfs,
1264  const cs_real_t *const l_visc,
1265  const cs_real_t *const prl,
1266  const cs_real_t *const prt,
1267  const cs_real_t *const rough_t,
1268  const cs_real_t *const uk,
1269  const cs_real_t *const yplus,
1270  const cs_real_t *const dplus,
1271  cs_real_t *htur,
1272  cs_real_t *yplim
1273 );
1274 
1275 /*=============================================================================
1276  * Public function prototypes
1277  *============================================================================*/
1278 
1279 /*----------------------------------------------------------------------------
1280  *! \brief Provide access to cs_glob_wall_functions
1281  *----------------------------------------------------------------------------*/
1282 
1285 
1286 /*----------------------------------------------------------------------------*/
1313 /*----------------------------------------------------------------------------*/
1314 
1315 void
1317  cs_real_t l_visc,
1318  cs_real_t t_visc,
1319  cs_real_t vel,
1320  cs_real_t y,
1321  cs_real_t rough_d,
1322  cs_real_t rnnb,
1323  cs_real_t kinetic_en,
1324  int *iuntur,
1325  cs_gnum_t *nsubla,
1326  cs_gnum_t *nlogla,
1327  cs_real_t *ustar,
1328  cs_real_t *uk,
1329  cs_real_t *yplus,
1330  cs_real_t *ypup,
1331  cs_real_t *cofimp,
1332  cs_real_t *dplus);
1333 
1334 /*----------------------------------------------------------------------------*/
1359 /*----------------------------------------------------------------------------*/
1360 
1361 void
1363  cs_real_t l_visc,
1364  cs_real_t prl,
1365  cs_real_t prt,
1366  cs_real_t rough_t,
1367  cs_real_t uk,
1368  cs_real_t yplus,
1369  cs_real_t dplus,
1370  cs_real_t *htur,
1371  cs_real_t *yplim);
1372 
1373 /*----------------------------------------------------------------------------*/
1381 /*----------------------------------------------------------------------------*/
1382 
1383 void
1385  cs_real_t *st_exp,
1386  cs_real_t *st_imp);
1387 
1388 /*----------------------------------------------------------------------------*/
1389 
1391 
1392 #endif /* __CS_WALL_FUNCTIONS_H__ */
int bft_printf(const char *const format,...)
Replacement for printf() with modifiable behavior.
Definition: bft_printf.c:140
#define BEGIN_C_DECLS
Definition: cs_defs.h:528
double cs_real_t
Floating-point value.
Definition: cs_defs.h:332
#define _(String)
Definition: cs_defs.h:63
#define CS_PROCF(x, y)
Definition: cs_defs.h:560
unsigned long cs_gnum_t
global mesh entity number
Definition: cs_defs.h:310
#define CS_UNUSED(x)
Definition: cs_defs.h:514
#define END_C_DECLS
Definition: cs_defs.h:529
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:325
@ eps
Definition: cs_field_pointer.h:71
@ vel
Definition: cs_field_pointer.h:68
@ yplus
Definition: cs_field_pointer.h:240
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:181
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:200
double cs_turb_vdriest
Definition: cs_turbulence_model.c:455
double cs_turb_crij2
Definition: cs_turbulence_model.c:545
double cs_turb_crij1
Definition: cs_turbulence_model.c:539
double cs_turb_cmu025
Definition: cs_turbulence_model.c:508
double cs_turb_cstlog
Definition: cs_turbulence_model.c:466
cs_turb_model_t * cs_get_glob_turb_model(void)
Provide write access to turbulence model structure.
Definition: cs_turbulence_model.c:1496
double cs_turb_dpow
Definition: cs_turbulence_model.c:497
double cs_turb_apow
Definition: cs_turbulence_model.c:491
double cs_turb_cmu
Definition: cs_turbulence_model.c:505
double cs_turb_bpow
Definition: cs_turbulence_model.c:494
double cs_turb_xkappa
Definition: cs_turbulence_model.c:446
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.c: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.c: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
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.c: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
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)
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
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)
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
cs_wall_functions_t * cs_get_glob_wall_functions(void)
Definition: cs_wall_functions.c:218
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:434
integer(c_int), pointer, save iwalfs
Wall functions for scalar.
Definition: optcal.f90:439
real(c_double), dimension(:), pointer, save a2
Definition: cpincl.f90:109
Turbulence model general options descriptor.
Definition: cs_turbulence_model.h:115
int itytur
Definition: cs_turbulence_model.h:139
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