8.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-2023 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_MAX(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_MIN(cs_turb_xkappa, CS_MAX(1., sqrt(ml_visc / t_visc)) / *yplus);
508  *cofimp = 1. - *ypup / cs_turb_xkappa * ( 2. * rcprod - 1. / (2. * *yplus));
509 
510  *nlogla += 1;
511 
512  /* viscous sub-layer */
513  } else {
514 
515  if (*yplus > 1.e-12) {
516  *ustar = fabs(vel / *yplus); /* FIXME remove that: its is here only to
517  be fully equivalent to the former code. */
518  } else {
519  *ustar = 0.;
520  }
521  *ypup = 1.;
522  *cofimp = 0.;
523 
524  *iuntur = 0;
525  *nsubla += 1;
526 
527  }
528 }
529 
530 /*----------------------------------------------------------------------------*/
551 /*----------------------------------------------------------------------------*/
552 
553 inline static void
555  cs_real_t t_visc,
556  cs_real_t vel,
557  cs_real_t y,
558  cs_real_t kinetic_en,
559  int *iuntur,
560  cs_gnum_t *nsubla,
561  cs_gnum_t *nlogla,
562  cs_real_t *ustar,
563  cs_real_t *uk,
564  cs_real_t *yplus,
565  cs_real_t *dplus,
566  cs_real_t *ypup,
567  cs_real_t *cofimp)
568 {
569  CS_UNUSED(iuntur);
570 
571  const double ypluli = cs_glob_wall_functions->ypluli;
572 
573  double rcprod, ml_visc, Re, g;
574  /* Compute the friction velocity ustar */
575 
576  /* Blending for very low values of k */
577  Re = sqrt(kinetic_en) * y / l_visc;
578  g = exp(-Re/11.);
579 
580  *uk = sqrt( (1.-g) * cs_turb_cmu025 * cs_turb_cmu025 * kinetic_en
581  + g * l_visc * vel / y);
582 
583  *yplus = *uk * y / l_visc;
584 
585  /* Compute the friction velocity ustar */
586  *uk = cs_turb_cmu025 * sqrt(kinetic_en);//FIXME
587  *yplus = *uk * y / l_visc;//FIXME
588 
589  /* Log layer */
590  if (*yplus > ypluli) {
591 
592  *dplus = 0.;
593 
594  *nlogla += 1;
595 
596  /* Viscous sub-layer and therefore shift */
597  } else {
598 
599  *dplus = ypluli - *yplus;
600 
601  /* Count the cell as if it was in the viscous sub-layer */
602  *nsubla += 1;
603 
604  }
605 
606  /* Mixing length viscosity */
607  ml_visc = cs_turb_xkappa * l_visc * (*yplus + *dplus);
608  rcprod = CS_MIN(cs_turb_xkappa, CS_MAX(1., sqrt(ml_visc / t_visc)) / (*yplus + *dplus));
609 
610  *ustar = vel / (log(*yplus + *dplus) / cs_turb_xkappa + cs_turb_cstlog);
611  *ypup = *yplus / (log(*yplus + *dplus) / cs_turb_xkappa + cs_turb_cstlog);
612  *cofimp = 1. - *ypup
613  / cs_turb_xkappa * (2. * rcprod - 1. / (2. * *yplus + *dplus));
614 }
615 
616 /*----------------------------------------------------------------------------
617  * Compute u+ for a given yk+ between 0.1 and 200 according to the two
618  * scales wall functions using Van Driest mixing length.
619  * This function holds the coefficients of the polynome fitting log(u+).
620  *
621  * parameters:
622  * yplus <-- dimensionless distance
623  *
624  * returns:
625  * the resulting dimensionless velocity.
626  *----------------------------------------------------------------------------*/
627 
628 inline static cs_real_t
630 {
631  /* Coefficients of the polynome fitting log(u+) for yk < 200 */
632  static double aa[11] = {-0.0091921, 3.9577, 0.031578,
633  -0.51013, -2.3254, -0.72665,
634  2.969, 0.48506, -1.5944,
635  0.087309, 0.1987 };
636 
637  cs_real_t y1,y2,y3,y4,y5,y6,y7,y8,y9,y10, uplus;
638 
639  y1 = 0.25 * log(yplus);
640  y2 = y1 * y1;
641  y3 = y2 * y1;
642  y4 = y3 * y1;
643  y5 = y4 * y1;
644  y6 = y5 * y1;
645  y7 = y6 * y1;
646  y8 = y7 * y1;
647  y9 = y8 * y1;
648  y10 = y9 * y1;
649 
650  uplus = aa[0]
651  + aa[1] * y1
652  + aa[2] * y2
653  + aa[3] * y3
654  + aa[4] * y4
655  + aa[5] * y5
656  + aa[6] * y6
657  + aa[7] * y7
658  + aa[8] * y8
659  + aa[9] * y9
660  + aa[10] * y10;
661 
662  return exp(uplus);
663 }
664 
665 /*----------------------------------------------------------------------------*/
700 /*----------------------------------------------------------------------------*/
701 
702 inline static void
704  cs_real_t l_visc,
705  cs_real_t vel,
706  cs_real_t y,
707  cs_real_t kinetic_en,
708  int *iuntur,
709  cs_gnum_t *nsubla,
710  cs_gnum_t *nlogla,
711  cs_real_t *ustar,
712  cs_real_t *uk,
713  cs_real_t *yplus,
714  cs_real_t *ypup,
715  cs_real_t *cofimp,
716  cs_real_t *lmk,
717  cs_real_t kr,
718  bool wf)
719 {
720  double urplus, d_up, lmk15;
721 
722  if (wf)
723  *uk = sqrt(sqrt((1.-cs_turb_crij2)/cs_turb_crij1 * rnnb * kinetic_en));
724 
725  /* Set a low threshold value in case tangential velocity is zero */
726  *yplus = CS_MAX(*uk * y / l_visc, 1.e-4);
727 
728  /* Dimensionless roughness */
729  cs_real_t krp = *uk * kr / l_visc;
730 
731  /* Extension of Van Driest mixing length according to Rotta (1962) with
732  Cebeci & Chang (1978) correlation */
733  cs_real_t dyrp = 0.9 * (sqrt(krp) - krp * exp(-krp / 6.));
734  cs_real_t yrplus = *yplus + dyrp;
735 
736  if (dyrp <= 1.e-1)
737  d_up = dyrp;
738  else if (dyrp <= 200.)
739  d_up = _vdriest_dupdyp_integral(dyrp);
740  else
741  d_up = 16.088739022054590 + log(dyrp/200.) / cs_turb_xkappa;
742 
743  if (yrplus <= 1.e-1) {
744 
745  urplus = yrplus;
746 
747  if (wf) {
748  *iuntur = 0;
749  *nsubla += 1;
750 
751  *lmk = 0.;
752 
753  *ypup = 1.;
754 
755  *cofimp = 0.;
756  }
757 
758  } else if (yrplus <= 200.) {
759 
760  urplus = _vdriest_dupdyp_integral(yrplus);
761 
762  if (wf) {
763  *nlogla += 1;
764 
765  *ypup = *yplus / (urplus-d_up);
766 
767  /* Mixing length in y+ */
768  *lmk = cs_turb_xkappa * (*yplus) *(1-exp(- (*yplus) / cs_turb_vdriest));
769 
770  /* Mixing length in 3/2*y+ */
771  lmk15 = cs_turb_xkappa * 1.5 * (*yplus) *(1-exp(- 1.5 * (*yplus)
772  / cs_turb_vdriest));
773 
774  *cofimp = 1. - (2. / (1. + *lmk) - 1. / (1. + lmk15)) * *ypup;
775  }
776 
777  } else {
778 
779  urplus = 16.088739022054590 + log(yrplus/200) / cs_turb_xkappa;
780 
781  if (wf) {
782  *nlogla += 1;
783 
784  *ypup = *yplus / (urplus-d_up);
785 
786  /* Mixing length in y+ */
787  *lmk = cs_turb_xkappa * (*yplus) *(1-exp(- (*yplus) / cs_turb_vdriest));
788 
789  /* Mixing length in 3/2*y+ */
790  lmk15 = cs_turb_xkappa * 1.5 * (*yplus) *(1-exp(- 1.5 * (*yplus)
791  / cs_turb_vdriest));
792 
793  *cofimp = 1. - (2. / *lmk - 1. / lmk15) * *ypup;
794  }
795 
796  }
797 
798  *ustar = vel / (urplus-d_up);
799 }
800 
801 /*----------------------------------------------------------------------------*/
834 /*----------------------------------------------------------------------------*/
835 
836 inline static void
838  cs_real_t t_visc,
839  cs_real_t vel,
840  cs_real_t y,
841  cs_real_t rough_d,
842  cs_real_t kinetic_en,
843  int *iuntur,
844  cs_gnum_t *nsubla,
845  cs_gnum_t *nlogla,
846  cs_real_t *ustar,
847  cs_real_t *uk,
848  cs_real_t *yplus,
849  cs_real_t *dplus,
850  cs_real_t *ypup,
851  cs_real_t *cofimp)
852 {
853  CS_UNUSED(iuntur);
854 
855  const double ypluli = cs_glob_wall_functions->ypluli;
856 
857  double rcprod, ml_visc, Re, g;
858 
859  /* Compute the friction velocity ustar */
860 
861  /* Shifting of the wall distance to be consistant with
862  * the fully rough wall function
863  *
864  * ln((y+y0)/y0) = ln((y+y0)/alpha xi) + kappa * 5.2
865  *
866  * y0 = xi * exp(-kappa * 8.5)
867  * where xi is the sand grain roughness here
868  * y0 = alpha * xi * exp(-kappa * 5.2)
869  *
870  * so:
871  * alpha = exp(-kappa * (8.5 - 5.2)) = 0.25
872  *
873  */
874  cs_real_t y0 = rough_d;
875  /* Note : Sand grain roughness given by:
876  cs_real_t sg_rough = rough_d * exp(cs_turb_xkappa*cs_turb_cstlog_rough);
877  */
878 
879  /* Blending for very low values of k */
880  Re = sqrt(kinetic_en) * (y + y0) / l_visc;
881  g = exp(-Re/11.);
882 
883  *uk = sqrt( (1.-g) * cs_turb_cmu025 * cs_turb_cmu025 * kinetic_en
884  + g * l_visc * vel / (y + y0));
885 
886 
887  *yplus = *uk * y / l_visc;
888 
889  /* As for scalable wall functions, yplus is shifted of "dplus" */
890  *dplus = *uk * y0 / l_visc;
891 
892  /* Shift of the velocity profile due to roughness */
893  cs_real_t shift_vel = -log(1. + y0 * exp(cs_turb_xkappa * cs_turb_cstlog) * *uk/l_visc)
894  / cs_turb_xkappa;
895 
896  /* Log layer and shifted with the roughness */
897  if (*yplus > ypluli) {
898 
899  *nlogla += 1;
900 
901  }
902 
903  /* Viscous sub-layer and therefore shift again */
904  else {
905 
906  *dplus = ypluli - *yplus;
907  /* Count the cell as if it was in the viscous sub-layer */
908  *nsubla += 1;
909 
910  }
911 
912  cs_real_t uplus = log(*yplus + *dplus) / cs_turb_xkappa + cs_turb_cstlog + shift_vel;
913  *ustar = vel / uplus;
914 #if 0
915  bft_printf("uet=%f, u=%f, uplus=%f, yk=%f, duplus=%f\n", *ustar, vel, uplus, *yplus, 1./uplus);
916 #endif
917  *ypup = *yplus / uplus;
918 
919  /* Mixing length viscosity, compatible with both regimes */
920  ml_visc = cs_turb_xkappa * l_visc * (*yplus + *dplus);
921  rcprod = CS_MIN(cs_turb_xkappa, CS_MAX(1., sqrt(ml_visc / t_visc)) / (*yplus + *dplus));
922  *cofimp = 1. - *yplus / (cs_turb_xkappa * uplus)
923  * ( 2. * rcprod - 1. / (2. * *yplus + *dplus));
924 
925 }
926 
927 /*----------------------------------------------------------------------------*/
947 /*----------------------------------------------------------------------------*/
948 
949 inline static void
951  cs_real_t t_visc,
952  cs_real_t vel,
953  cs_real_t y,
954  int *iuntur,
955  cs_gnum_t *nsubla,
956  cs_gnum_t *nlogla,
957  cs_real_t *ustar,
958  cs_real_t *uk,
959  cs_real_t *yplus,
960  cs_real_t *dplus,
961  cs_real_t *ypup,
962  cs_real_t *cofimp)
963 {
964  CS_UNUSED(t_visc);
965  CS_UNUSED(nlogla);
966  CS_UNUSED(dplus);
967 
968  const double ypluli = cs_glob_wall_functions->ypluli;
969 
970  /* Compute the friction velocity ustar */
971 
972  *ustar = sqrt(vel * l_visc / y);
973  *yplus = *ustar * y / l_visc;
974  *uk = *ustar;
975  *ypup = 1.;
976  *cofimp = 0.;
977  *iuntur = 0;
978 
979  if (*yplus <= ypluli) {
980 
981  /* Disable the wall funcion count the cell in the viscous sub-layer */
982  *nsubla += 1;
983 
984  } else {
985 
986  /* Count the cell as if it was in the viscous sub-layer */
987  *nsubla += 1;
988 
989  }
990 }
991 
992 /*----------------------------------------------------------------------------*/
1022 /*----------------------------------------------------------------------------*/
1023 
1024 inline static void
1026  cs_real_t prl,
1027  cs_real_t prt,
1028  cs_real_t rough_t,
1029  cs_real_t uk,
1030  cs_real_t yplus,
1031  cs_real_t dplus,
1032  cs_real_t *htur,
1033  cs_real_t *yplim)
1034 {
1035  /* Local variables */
1036  double tplus;
1037  double beta2,a2;
1038  double yp2;
1039  double prlm1;
1040 
1041  const double epzero = cs_math_epzero;
1042 
1043  /*==========================================================================
1044  1. Initializations
1045  ==========================================================================*/
1046 
1047  (*htur) = CS_MAX(yplus,epzero)/CS_MAX(yplus+dplus,epzero);
1048 
1049  prlm1 = 0.1;
1050 
1051  /* Sand grain roughness is:
1052  * zeta = z0 * exp(kappa 8.5)
1053  * Then:
1054  * hp = zeta uk / nu * exp( -kappa(8.5 - 5.2))
1055  * = z0 * uk / nu * exp(kappa * 5.2)
1056  * where 5.2 is the smooth log constant, and 8.5 the rough one
1057  *
1058  * FIXME check if we should use a molecular Schmidt number
1059  */
1060  cs_real_t hp = rough_t *uk / l_visc * exp(cs_turb_xkappa*cs_turb_cstlog);
1061 
1062  /* Shift of the temperature profile due to roughness */
1063  cs_real_t shift_temp = -log(1. + hp);
1064 
1065  /*==========================================================================
1066  2. Compute htur for small Prandtl numbers
1067  ==========================================================================*/
1068 
1069  if (prl <= prlm1) {
1070  (*yplim) = prt/(prl*cs_turb_xkappa);
1071  if (yplus > (*yplim)) {
1072  tplus = prl*(*yplim) + prt/cs_turb_xkappa * (log((yplus+dplus)/(*yplim)) + shift_temp);
1073  (*htur) = prl * yplus / tplus;
1074  }
1075 
1076  /*========================================================================
1077  3. Compute htur for the model with three sub-layers
1078  ========================================================================*/
1079 
1080  } else {
1081  yp2 = sqrt(cs_turb_xkappa*1000./prt);
1082  (*yplim) = pow(1000./prl, 1./3.);
1083 
1084  a2 = 15.*pow(prl, 2./3.);
1085 
1086  if (yplus >= (*yplim) && yplus < yp2) {
1087  tplus = a2 - 500./((yplus+dplus)*(yplus+dplus));
1088  (*htur) = prl * yplus / tplus;
1089  }
1090 
1091  if (yplus >= yp2) {
1092  beta2 = a2 - 0.5 * prt /cs_turb_xkappa;
1093  tplus = beta2 + prt/cs_turb_xkappa*log((yplus+dplus)/yp2);
1094  (*htur) = prl * yplus / tplus;
1095  }
1096 
1097  }
1098 }
1099 
1100 /*----------------------------------------------------------------------------*/
1124 /*----------------------------------------------------------------------------*/
1125 
1126 inline static void
1128  cs_real_t prt,
1129  cs_real_t yplus,
1130  cs_real_t *htur)
1131 {
1132  cs_real_t prlrat = prl / prt;
1133 
1134  /* Parameters of the numerical quadrature */
1135  const int ninter_max = 100;
1136  const cs_real_t ypmax = 1.e2;
1137 
1138  /* No correction for very small yplus */
1139  if (yplus <= 0.1)
1140  *htur = 1.;
1141  else {
1142  cs_real_t ypint = CS_MIN(yplus, ypmax);
1143 
1144  /* The number of sub-intervals is taken proportional to yplus and equal to
1145  * ninter_max if yplus=ypmax */
1146 
1147  int npeff = CS_MAX((int)(ypint / ypmax * (double)(ninter_max)), 1);
1148 
1149  double dy = ypint / (double)(npeff);
1150  cs_real_t stplus = 0.;
1151  cs_real_t nut1 = 0.;
1152  cs_real_t nut2 = 0.;
1153 
1154  for (int ip = 1; ip <= npeff; ip++) {
1155  double yp = ypint * (double)(ip) / (double)(npeff);
1156  nut2 = cs_turb_xkappa * yp * (1. - exp(-yp / cs_turb_vdriest));
1157  stplus += dy / (1. + prlrat * 0.5 * (nut1 + nut2));
1158  nut1 = nut2;
1159  }
1160 
1161  if (yplus > ypint) {
1162  cs_real_t r = prlrat * cs_turb_xkappa;
1163  stplus += log( (1. + r*yplus) / (1. + r*ypint)) / r;
1164  }
1165 
1166  if (stplus >= 1.e-6)
1167  *htur = yplus / stplus;
1168  else
1169  *htur = 1.;
1170  }
1171 }
1172 
1173 /*----------------------------------------------------------------------------*/
1189 /*----------------------------------------------------------------------------*/
1190 
1191 inline static void
1193  cs_real_t prl,
1194  cs_real_t prt,
1195  cs_real_t rough_t,
1196  cs_real_t uk,
1197  cs_real_t yplus,
1198  cs_real_t dplus,
1199  cs_real_t *htur)
1200 {
1201  CS_UNUSED(prt);
1202 
1203  /* Sand grain roughness is:
1204  * zeta = z0 * exp(kappa 8.5)
1205  * Then:
1206  * hp = zeta uk / nu * exp( -kappa(8.5 - 5.2))
1207  * = z0 * uk / nu * exp(kappa * 5.2)
1208  * where 5.2 is the smooth log constant, and 8.5 the rough one
1209  *
1210  * FIXME check if we should use a molecular Schmidt number
1211  */
1212  cs_real_t hp = rough_t *uk / l_visc * exp(cs_turb_xkappa*cs_turb_cstlog);
1213  const double ypluli = cs_glob_wall_functions->ypluli;
1214  const double epzero = cs_math_epzero;
1215 
1216  (*htur) = CS_MAX(yplus,epzero)/CS_MAX(yplus+dplus,epzero);
1217 
1218  /* Shift of the temperature profile due to roughness */
1219  cs_real_t shift_temp = -log(1. + hp);
1220 
1221  if (yplus > ypluli) {
1222  cs_real_t tplus = prt * ((log(yplus+dplus) + shift_temp)/cs_turb_xkappa + 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:514
double cs_real_t
Floating-point value.
Definition: cs_defs.h:319
#define CS_MIN(a, b)
Definition: cs_defs.h:477
#define _(String)
Definition: cs_defs.h:63
#define CS_MAX(a, b)
Definition: cs_defs.h:478
#define CS_PROCF(x, y)
Definition: cs_defs.h:528
unsigned long cs_gnum_t
global mesh entity number
Definition: cs_defs.h:298
#define CS_UNUSED(x)
Definition: cs_defs.h:500
#define END_C_DECLS
Definition: cs_defs.h:515
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:313
@ eps
Definition: cs_field_pointer.h:71
@ vel
Definition: cs_field_pointer.h:68
@ yplus
Definition: cs_field_pointer.h:239
const cs_real_t cs_math_epzero
double cs_turb_vdriest
Definition: cs_turbulence_model.c:453
double cs_turb_crij2
Definition: cs_turbulence_model.c:543
double cs_turb_crij1
Definition: cs_turbulence_model.c:537
double cs_turb_cmu025
Definition: cs_turbulence_model.c:506
double cs_turb_cstlog
Definition: cs_turbulence_model.c:464
cs_turb_model_t * cs_get_glob_turb_model(void)
Provide write access to turbulence model structure.
Definition: cs_turbulence_model.c:1494
double cs_turb_dpow
Definition: cs_turbulence_model.c:495
double cs_turb_apow
Definition: cs_turbulence_model.c:489
double cs_turb_cmu
Definition: cs_turbulence_model.c:503
double cs_turb_bpow
Definition: cs_turbulence_model.c:492
double cs_turb_xkappa
Definition: cs_turbulence_model.c:444
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:537
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:608
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:837
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:1192
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:629
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:1025
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:950
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:350
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:703
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:1127
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:281
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)
Definition: cs_wall_functions.c:216
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:554
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:317
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
double precision, dimension(ncharm), save a2
Definition: cpincl.f90:235
double precision, dimension(4, npot), save ka
Definition: ppcpfu.f90:163
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