programmer's 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-2015 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_turbulence_model.h"
42 
43 /*----------------------------------------------------------------------------*/
44 
46 
47 /*=============================================================================
48  * Local Macro definitions
49  *============================================================================*/
50 
51 /*============================================================================
52  * Type definition
53  *============================================================================*/
54 
55 /* Wall function type */
56 /*--------------------*/
57 
58 typedef enum {
59 
67 
69 
70 typedef enum {
71 
74 
76 
77 /* Wall functions descriptor */
78 /*---------------------------*/
79 
80 typedef struct {
81 
82  cs_wall_f_type_t iwallf; /* wall function type */
83 
84  cs_wall_f_s_type_t iwalfs; /* wall function type for scalars */
85 
86  int iwallt; /* exchange coefficient correlation
87  - 0: not used by default
88  - 1: exchange coefficient computed with a
89  correlation */
90 
91  double ypluli; /* limit value of y+ for the viscous
92  sublayer */
93 
95 
96 /*============================================================================
97  * Global variables
98  *============================================================================*/
99 
100 /* Pointer to wall functions descriptor structure */
101 
103 
104 /*============================================================================
105  * Private function definitions
106  *============================================================================*/
107 
108 /*----------------------------------------------------------------------------*/
125 /*----------------------------------------------------------------------------*/
126 
127 inline static void
129  cs_real_t vel,
130  cs_real_t y,
131  int *iuntur,
132  cs_lnum_t *nsubla,
133  cs_lnum_t *nlogla,
134  cs_real_t *ustar,
135  cs_real_t *uk,
136  cs_real_t *yplus,
137  cs_real_t *ypup,
138  cs_real_t *cofimp)
139 {
140  const double ypluli = cs_glob_wall_functions->ypluli;
141 
142  const double ydvisc = y / l_visc;
143 
144  /* Compute the friction velocity ustar */
145 
146  *ustar = pow((vel/(cs_turb_apow * pow(ydvisc, cs_turb_bpow))), cs_turb_dpow);
147  *uk = *ustar;
148  *yplus = *ustar * ydvisc;
149 
150  /* In the viscous sub-layer: U+ = y+ */
151  if (*yplus <= ypluli) {
152 
153  *ustar = sqrt(vel / ydvisc);
154  *yplus = *ustar * ydvisc;
155  *uk = *ustar;
156  *ypup = 1.;
157  *cofimp = 0.;
158 
159  /* Disable the wall funcion count the cell in the viscous sub-layer */
160  *iuntur = 0;
161  *nsubla += 1;
162 
163  /* In the log layer */
164  } else {
165  *ypup = pow(vel, 2. * cs_turb_dpow-1.)
166  / pow(cs_turb_apow, 2. * cs_turb_dpow);
167  *cofimp = 1. + cs_turb_bpow
168  * pow(*ustar, cs_turb_bpow + 1. - 1./cs_turb_dpow)
169  * (pow(2., cs_turb_bpow - 1.) - 2.);
170 
171  /* Count the cell in the log layer */
172  *nlogla += 1;
173 
174  }
175 }
176 
177 /*----------------------------------------------------------------------------*/
196 /*----------------------------------------------------------------------------*/
197 
198 inline static void
200  cs_real_t l_visc,
201  cs_real_t vel,
202  cs_real_t y,
203  int *iuntur,
204  cs_lnum_t *nsubla,
205  cs_lnum_t *nlogla,
206  cs_real_t *ustar,
207  cs_real_t *uk,
208  cs_real_t *yplus,
209  cs_real_t *ypup,
210  cs_real_t *cofimp)
211 {
212  const double ypluli = cs_glob_wall_functions->ypluli;
213 
214  double ustarwer, ustarmin, ustaro, ydvisc;
215  double eps = 0.001;
216  int niter_max = 100;
217  int iter = 0;
218  double reynolds;
219 
220  /* Compute the local Reynolds number */
221 
222  ydvisc = y / l_visc;
223  reynolds = vel * ydvisc;
224 
225  /*
226  * Compute the friction velocity ustar
227  */
228 
229  /* In the viscous sub-layer: U+ = y+ */
230  if (reynolds <= ypluli * ypluli) {
231 
232  *ustar = sqrt(vel / ydvisc);
233  *yplus = *ustar * ydvisc;
234  *uk = *ustar;
235  *ypup = 1.;
236  *cofimp = 0.;
237 
238  /* Disable the wall funcion count the cell in the viscous sub-layer */
239  *iuntur = 0;
240  *nsubla += 1;
241 
242  /* In the log layer */
243  } else {
244 
245  /* The initial value is Wener or the minimun ustar to ensure convergence */
246  ustarwer = pow(fabs(vel) / cs_turb_apow / pow(ydvisc, cs_turb_bpow),
247  cs_turb_dpow);
248  ustarmin = exp(-cs_turb_cstlog * cs_turb_xkappa)/ydvisc;
249  ustaro = CS_MAX(ustarwer, ustarmin);
250  *ustar = (cs_turb_xkappa * vel + ustaro)
251  / (log(ydvisc * ustaro) + cs_turb_xkappa * cs_turb_cstlog + 1.);
252 
253  /* Iterative solving */
254  for (iter = 0; iter < niter_max
255  && fabs(*ustar - ustaro) >= eps * ustaro; iter++) {
256  ustaro = *ustar;
257  *ustar = (cs_turb_xkappa * vel + ustaro)
258  / (log(ydvisc * ustaro) + cs_turb_xkappa * cs_turb_cstlog + 1.);
259  }
260 
261  if (iter >= niter_max) {
262  bft_printf(_("WARNING: non-convergence in the computation\n"
263  "******** of the friction velocity\n\n"
264  "face number: %d \n"
265  "friction vel: %f \n" ), ifac, *ustar);
266  }
267 
268  *uk = *ustar;
269  *yplus = *ustar * ydvisc;
270  *ypup = *yplus / (log(*yplus) / cs_turb_xkappa + cs_turb_cstlog);
271  *cofimp = 1. - *ypup / cs_turb_xkappa * 1.5 / *yplus;
272 
273  /* Count the cell in the log layer */
274  *nlogla += 1;
275 
276  }
277 
278 }
279 
280 /*----------------------------------------------------------------------------*/
300 /*----------------------------------------------------------------------------*/
301 
302 inline static void
304  cs_real_t t_visc,
305  cs_real_t vel,
306  cs_real_t y,
307  cs_real_t kinetic_en,
308  int *iuntur,
309  cs_lnum_t *nsubla,
310  cs_lnum_t *nlogla,
311  cs_real_t *ustar,
312  cs_real_t *uk,
313  cs_real_t *yplus,
314  cs_real_t *ypup,
315  cs_real_t *cofimp)
316 {
317  const double ypluli = cs_glob_wall_functions->ypluli;
318 
319  double rcprod, ml_visc, Re, g;
320 
321  /* Compute the friction velocity ustar */
322 
323  /* Blending for very low values of k */
324  Re = sqrt(kinetic_en) * y / l_visc;
325  g = exp(-Re/11.);
326 
327  *uk = sqrt( (1.-g) * cs_turb_cmu025 * cs_turb_cmu025 * kinetic_en
328  + g * l_visc * vel / y);
329 
330  *yplus = *uk * y / l_visc;
331 
332  /* log layer */
333  if (*yplus > ypluli) {
334 
335  *ustar = vel / (log(*yplus) / cs_turb_xkappa + cs_turb_cstlog);
336  *ypup = *yplus / (log(*yplus) / cs_turb_xkappa + cs_turb_cstlog);
337  /* Mixing length viscosity */
338  ml_visc = cs_turb_xkappa * l_visc * *yplus;
339  rcprod = CS_MIN(cs_turb_xkappa, CS_MAX(1., sqrt(ml_visc / t_visc)) / *yplus);
340  *cofimp = 1. - *ypup / cs_turb_xkappa * ( 2. * rcprod - 1. / (2. * *yplus));
341 
342  *nlogla += 1;
343 
344  /* viscous sub-layer */
345  } else {
346 
347  if (*yplus > 1.e-12) {
348  *ustar = fabs(vel / *yplus); /* FIXME remove that: its is here only to
349  be fully equivalent to the former code. */
350  } else {
351  *ustar = 0.;
352  }
353  *ypup = 1.;
354  *cofimp = 0.;
355 
356  *iuntur = 0;
357  *nsubla += 1;
358 
359  }
360 }
361 
362 /*----------------------------------------------------------------------------*/
383 /*----------------------------------------------------------------------------*/
384 
385 inline static void
387  cs_real_t t_visc,
388  cs_real_t vel,
389  cs_real_t y,
390  cs_real_t kinetic_en,
391  int *iuntur,
392  cs_lnum_t *nsubla,
393  cs_lnum_t *nlogla,
394  cs_real_t *ustar,
395  cs_real_t *uk,
396  cs_real_t *yplus,
397  cs_real_t *dplus,
398  cs_real_t *ypup,
399  cs_real_t *cofimp)
400 {
401  const double ypluli = cs_glob_wall_functions->ypluli;
402 
403  double rcprod, ml_visc, Re, g;
404  /* Compute the friction velocity ustar */
405 
406  /* Blending for very low values of k */
407  Re = sqrt(kinetic_en) * y / l_visc;
408  g = exp(-Re/11.);
409 
410  *uk = sqrt( (1.-g) * cs_turb_cmu025 * cs_turb_cmu025 * kinetic_en
411  + g * l_visc * vel / y);
412 
413  *yplus = *uk * y / l_visc;
414 
415  /* Compute the friction velocity ustar */
416  *uk = cs_turb_cmu025 * sqrt(kinetic_en);
417  *yplus = *uk * y / l_visc;
418 
419  /* Log layer */
420  if (*yplus > ypluli) {
421 
422  *dplus = 0.;
423 
424  *nlogla += 1;
425 
426  /* Viscous sub-layer and therefore shift */
427  } else {
428 
429  *dplus = ypluli - *yplus;
430  *yplus = ypluli;
431 
432  /* Count the cell as if it was in the viscous sub-layer */
433  *nsubla += 1;
434 
435  }
436 
437  /* Mixing length viscosity */
438  ml_visc = cs_turb_xkappa * l_visc * *yplus;
439  rcprod = CS_MIN(cs_turb_xkappa, CS_MAX(1., sqrt(ml_visc / t_visc)) / *yplus);
440 
441  *ustar = vel / (log(*yplus) / cs_turb_xkappa + cs_turb_cstlog);
442  *ypup = (*yplus - *dplus) / (log(*yplus) / cs_turb_xkappa + cs_turb_cstlog);
443  *cofimp = 1. - *ypup
444  / cs_turb_xkappa * (2. * rcprod - 1. / (2. * *yplus - *dplus));
445 }
446 
447 /*----------------------------------------------------------------------------
448  * Compute u+ for a given yk+ between 0.1 and 200 according to the two
449  * scales wall functions using Van Driest mixing length.
450  * This function holds the coefficients of the polynome fitting log(u+).
451  *
452  * parameters:
453  * yplus <-- dimensionless distance
454  *
455  * returns:
456  * the resulting dimensionless velocity.
457  *----------------------------------------------------------------------------*/
458 
459 inline static cs_real_t
461 {
462  /* Coefficients of the polynome fitting log(u+) for yk < 200 */
463  static double aa[11] = {-0.0091921, 3.9577, 0.031578,
464  -0.51013, -2.3254, -0.72665,
465  2.969, 0.48506, -1.5944,
466  0.087309, 0.1987 };
467 
468  cs_real_t y1,y2,y3,y4,y5,y6,y7,y8,y9,y10, uplus;
469 
470  y1 = 0.25 * log(yplus);
471  y2 = y1 * y1;
472  y3 = y2 * y1;
473  y4 = y3 * y1;
474  y5 = y4 * y1;
475  y6 = y5 * y1;
476  y7 = y6 * y1;
477  y8 = y7 * y1;
478  y9 = y8 * y1;
479  y10 = y9 * y1;
480 
481  uplus = aa[0]
482  + aa[1] * y1
483  + aa[2] * y2
484  + aa[3] * y3
485  + aa[4] * y4
486  + aa[5] * y5
487  + aa[6] * y6
488  + aa[7] * y7
489  + aa[8] * y8
490  + aa[9] * y9
491  + aa[10] * y10;
492 
493  return exp(uplus);
494 }
495 
496 /*----------------------------------------------------------------------------*/
531 /*----------------------------------------------------------------------------*/
532 
533 inline static void
535  cs_real_t l_visc,
536  cs_real_t vel,
537  cs_real_t y,
538  cs_real_t kinetic_en,
539  int *iuntur,
540  cs_lnum_t *nsubla,
541  cs_lnum_t *nlogla,
542  cs_real_t *ustar,
543  cs_real_t *uk,
544  cs_real_t *yplus,
545  cs_real_t *ypup,
546  cs_real_t *cofimp,
547  cs_real_t *lmk,
548  cs_real_t kr,
549  bool wf)
550 {
551  double urplus, dup, lmk15;
552 
553  if (wf)
554  *uk = sqrt(sqrt((1.-cs_turb_crij2)/cs_turb_crij1 * rnnb * kinetic_en));
555 
556  /* Set a low threshold value in case tangential velocity is zero */
557  *yplus = CS_MAX(*uk * y / l_visc, 1.e-4);
558 
559  /* Dimensionless roughness */
560  cs_real_t krp = *uk * kr / l_visc;
561 
562  /* Extension of Van Driest mixing length according to Rotta (1962) with
563  Cebeci & Chang (1978) correlation */
564  cs_real_t dyrp = 0.9 * (sqrt(krp) - krp * exp(-krp / 6.));
565  cs_real_t yrplus = *yplus + dyrp;
566 
567  if (dyrp <= 1.e-1)
568  dup = dyrp;
569  else if (dyrp <= 200.)
570  dup = _vdriest_dupdyp_integral(dyrp);
571  else
572  dup = 16.088739022054590 + log(dyrp/200.) / cs_turb_xkappa;
573 
574  if (yrplus <= 1.e-1) {
575 
576  urplus = yrplus;
577 
578  if (wf) {
579  *iuntur = 0;
580  *nsubla += 1;
581 
582  *lmk = 0.;
583 
584  *ypup = 1.;
585 
586  *cofimp = 0.;
587  }
588 
589  } else if (yrplus <= 200.) {
590 
591  urplus = _vdriest_dupdyp_integral(yrplus);
592 
593  if (wf) {
594  *nlogla += 1;
595 
596  *ypup = *yplus / (urplus-dup);
597 
598  /* Mixing length in y+ */
599  *lmk = cs_turb_xkappa * (*yplus) *(1-exp(- (*yplus) / cs_turb_vdriest));
600 
601  /* Mixing length in 3/2*y+ */
602  lmk15 = cs_turb_xkappa * 1.5 * (*yplus) *(1-exp(- 1.5 * (*yplus)
603  / cs_turb_vdriest));
604 
605  *cofimp = 1. - (2. / (1. + *lmk) - 1. / (1. + lmk15)) * *ypup;
606  }
607 
608  } else {
609 
610  urplus = 16.088739022054590 + log(yrplus/200) / cs_turb_xkappa;
611 
612  if (wf) {
613  *nlogla += 1;
614 
615  *ypup = *yplus / (urplus-dup);
616 
617  /* Mixing length in y+ */
618  *lmk = cs_turb_xkappa * (*yplus) *(1-exp(- (*yplus) / cs_turb_vdriest));
619 
620  /* Mixing length in 3/2*y+ */
621  lmk15 = cs_turb_xkappa * 1.5 * (*yplus) *(1-exp(- 1.5 * (*yplus)
622  / cs_turb_vdriest));
623 
624  *cofimp = 1. - (2. / *lmk - 1. / lmk15) * *ypup;
625  }
626 
627  }
628 
629  *ustar = vel / (urplus-dup);
630 }
631 
632 /*----------------------------------------------------------------------------*/
661 /*----------------------------------------------------------------------------*/
662 
663 inline static void
665  cs_real_t t_visc,
666  cs_real_t vel,
667  cs_real_t y,
668  cs_real_t roughness,
669  cs_real_t kinetic_en,
670  int *iuntur,
671  cs_lnum_t *nsubla,
672  cs_lnum_t *nlogla,
673  cs_real_t *ustar,
674  cs_real_t *uk,
675  cs_real_t *yplus,
676  cs_real_t *dplus,
677  cs_real_t *ypup,
678  cs_real_t *cofimp)
679 {
680  const double ypluli = cs_glob_wall_functions->ypluli;
681 
682  double rcprod, ml_visc, Re, g;
683 
684  /* Compute the friction velocity ustar */
685 
686  /* Shifting of the wall distance to be consistant with
687  * the fully rough wall function
688  *
689  * ln((y+y0)/y0) = ln((y+y0)/alpha xi) + kappa * 5.2
690  *
691  * y0 = roughness * exp(-kappa * 8.5)
692  * */
693  double y0 = roughness*exp(-cs_turb_xkappa*cs_turb_cstlog_rough);
694 
695  /* Blending for very low values of k */
696  Re = sqrt(kinetic_en) * (y + y0) / l_visc;
697  g = exp(-Re/11.);
698 
699  *uk = sqrt( (1.-g) * cs_turb_cmu025 * cs_turb_cmu025 * kinetic_en
700  + g * l_visc * vel / (y + y0));
701 
702  double effective_visc = (l_visc + cs_turb_cstlog_alpha * roughness * *uk);
703 
704  *yplus = *uk * (y + y0) / effective_visc;
705 
706  /* NB: tends to "y/xi" in rough regime, to "y.uk/nu" in smooth regime */
707  double yk = *uk * y / l_visc;
708 
709  /* As for scalable wall functions, "y*uk/effective_visc+ dplus = yplus" */
710  *dplus = *uk * y0 / effective_visc;
711 
712  /* Log layer and shifted with the roughness */
713  if (yk > ypluli) {
714 
715  *nlogla += 1;
716 
717  /* Viscous sub-layer and therefore shift again */
718  } else {
719 
720  *dplus = ypluli - *yplus;
721  *yplus = ypluli;
722  /* Count the cell as if it was in the viscous sub-layer */
723  *nsubla += 1;
724 
725  }
726 
727  double uplus = log(*yplus) / cs_turb_xkappa + cs_turb_cstlog;
728  *ustar = vel / uplus;
729  *ypup = yk / uplus;
730 
731  /* Mixing length viscosity, compatible with both regimes */
732  ml_visc = cs_turb_xkappa * *uk * (y + y0);
733  rcprod = CS_MIN(cs_turb_xkappa, CS_MAX(1., sqrt(ml_visc / t_visc)) / *yplus);
734  *cofimp = 1. - (*yplus - *dplus) / uplus
735  / cs_turb_xkappa * ( 2. * rcprod - 1. / (2. * *yplus - *dplus));
736 
737 }
738 
739 /*----------------------------------------------------------------------------*/
759 /*----------------------------------------------------------------------------*/
760 
761 inline static void
763  cs_real_t t_visc,
764  cs_real_t vel,
765  cs_real_t y,
766  int *iuntur,
767  cs_lnum_t *nsubla,
768  cs_lnum_t *nlogla,
769  cs_real_t *ustar,
770  cs_real_t *uk,
771  cs_real_t *yplus,
772  cs_real_t *dplus,
773  cs_real_t *ypup,
774  cs_real_t *cofimp)
775 {
776  const double ypluli = cs_glob_wall_functions->ypluli;
777 
778  /* Compute the friction velocity ustar */
779 
780  *ustar = sqrt(vel * l_visc / y);
781  *yplus = *ustar * y / l_visc;
782  *uk = *ustar;
783  *ypup = 1.;
784  *cofimp = 0.;
785  *iuntur = 0;
786 
787  if (*yplus <= ypluli) {
788 
789  /* Disable the wall funcion count the cell in the viscous sub-layer */
790  *nsubla += 1;
791 
792  } else {
793 
794  /* Count the cell as if it was in the viscous sub-layer */
795  *nsubla += 1;
796 
797  }
798 }
799 
800 /*----------------------------------------------------------------------------*/
826 /*----------------------------------------------------------------------------*/
827 
828 inline static void
830  cs_real_t prt,
831  cs_real_t yplus,
832  cs_real_t dplus,
833  cs_real_t *htur,
834  cs_real_t *yplim)
835 {
836  /* Local variables */
837  double tplus;
838  double beta2,a2;
839  double yp2;
840  double prlm1;
841 
842  const double epzero = 1.e-12;
843 
844  /*==========================================================================*/
845 
846  /*==========================================================================
847  1. Initializations
848  ==========================================================================*/
849 
850  /*==========================================================================*/
851 
852  (*htur) = CS_MAX(yplus-dplus,epzero)/CS_MAX(yplus,epzero);
853 
854  prlm1 = 0.1;
855 
856  /*==========================================================================
857  2. Compute htur for small Prandtl numbers
858  ==========================================================================*/
859 
860  if (prl <= prlm1) {
861  (*yplim) = prt/(prl*cs_turb_xkappa);
862  if (yplus > (*yplim)) {
863  tplus = prl*(*yplim) + prt/cs_turb_xkappa * log(yplus/(*yplim));
864  (*htur) = prl*(yplus-dplus)/tplus;
865  }
866 
867  /*========================================================================
868  3. Compute htur for the model with three sub-layers
869  ========================================================================*/
870 
871  } else {
872  yp2 = cs_turb_xkappa*1000./prt;
873  yp2 = sqrt(yp2);
874  (*yplim) = pow(1000./prl,1./3.);
875 
876  a2 = 15.*pow(prl,2./3.);
877  beta2 = a2 - 500./ pow(yp2,2);
878 
879  if (yplus >= (*yplim) && yplus < yp2) {
880  tplus = a2 - 500./(yplus*yplus);
881  (*htur) = prl*(yplus-dplus)/tplus;
882  }
883 
884  if (yplus >= yp2) {
885  tplus = beta2 + prt/cs_turb_xkappa*log(yplus/yp2);
886  (*htur) = prl*(yplus-dplus)/tplus;
887  }
888 
889  }
890 }
891 
892 /*----------------------------------------------------------------------------*/
915 /*----------------------------------------------------------------------------*/
916 
917 inline static void
919  cs_real_t prt,
920  cs_real_t yplus,
921  cs_real_t *htur)
922 {
923  cs_real_t prlrat = prl / prt;
924 
925  /* Parameters of the numerical quadrature */
926  const int ninter_max = 100;
927  const cs_real_t ypmax = 1.e2;
928 
929  /* No correction for very small yplus */
930  if (yplus <= 0.1)
931  *htur = 1.;
932  else {
933  cs_real_t ypint = CS_MIN(yplus, ypmax);
934 
935  /* The number of sub-intervals is taken proportional to yplus and equal to
936  * ninter_max if yplus=ypmax */
937 
938  int npeff = CS_MAX((int)(ypint / ypmax * (double)(ninter_max)), 1);
939 
940  double dy = ypint / (double)(npeff);
941  cs_real_t stplus = 0.;
942  cs_real_t nut1 = 0.;
943  cs_real_t nut2 = 0.;
944 
945  for (int ip = 1; ip <= npeff; ip++) {
946  double yp = ypint * (double)(ip) / (double)(npeff);
947  nut2 = cs_turb_xkappa * yp * (1. - exp(-yp / cs_turb_vdriest));
948  stplus += dy / (1. + prlrat * 0.5 * (nut1 + nut2));
949  nut1 = nut2;
950  }
951 
952  if (yplus > ypint) {
953  cs_real_t r = prlrat * cs_turb_xkappa;
954  stplus += log( (1. + r*yplus) / (1. + r*ypint)) / r;
955  }
956 
957  if (stplus >= 1.e-6)
958  *htur = yplus / stplus;
959  else
960  *htur = 1.;
961  }
962 }
963 
964 /*============================================================================
965  * Public function definitions for Fortran API
966  *============================================================================*/
967 
968 /*----------------------------------------------------------------------------
969  * Wrapper to cs_wall_functions_velocity.
970  *----------------------------------------------------------------------------*/
971 
972 void CS_PROCF (wallfunctions, WALLFUNCTIONS)
973 (
974  const cs_int_t *const iwallf,
975  const cs_lnum_t *const ifac,
976  const cs_real_t *const viscosity,
977  const cs_real_t *const t_visc,
978  const cs_real_t *const vel,
979  const cs_real_t *const y,
980  const cs_real_t *const roughness,
981  const cs_real_t *const rnnb,
982  const cs_real_t *const kinetic_en,
983  cs_int_t *iuntur,
984  cs_lnum_t *nsubla,
985  cs_lnum_t *nlogla,
986  cs_real_t *ustar,
987  cs_real_t *uk,
988  cs_real_t *yplus,
989  cs_real_t *ypup,
990  cs_real_t *cofimp,
991  cs_real_t *dplus
992 );
993 
994 /*----------------------------------------------------------------------------
995  * Wrapper to cs_wall_functions_scalar.
996  *----------------------------------------------------------------------------*/
997 
998 void CS_PROCF (hturbp, HTURBP)
999 (
1000  const cs_int_t *const iwalfs,
1001  const cs_real_t *const prl,
1002  const cs_real_t *const prt,
1003  const cs_real_t *const yplus,
1004  const cs_real_t *const dplus,
1005  cs_real_t *htur,
1006  cs_real_t *yplim
1007 );
1008 
1009 /*=============================================================================
1010  * Public function prototypes
1011  *============================================================================*/
1012 
1013 /*----------------------------------------------------------------------------
1014  *! \brief Provide acces to cs_glob_wall_functions
1015  *
1016  * needed to initialize structure with GUI
1017  *----------------------------------------------------------------------------*/
1018 
1021 
1025 /*-------------------------------------------------------------------------------
1026  Arguments
1027  ______________________________________________________________________________.
1028  mode name role !
1029  ______________________________________________________________________________*/
1055 /*-------------------------------------------------------------------------------*/
1056 
1057 void
1059  cs_lnum_t ifac,
1060  cs_real_t l_visc,
1061  cs_real_t t_visc,
1062  cs_real_t vel,
1063  cs_real_t y,
1064  cs_real_t roughness,
1065  cs_real_t rnnb,
1066  cs_real_t kinetic_en,
1067  int *iuntur,
1068  cs_lnum_t *nsubla,
1069  cs_lnum_t *nlogla,
1070  cs_real_t *ustar,
1071  cs_real_t *uk,
1072  cs_real_t *yplus,
1073  cs_real_t *ypup,
1074  cs_real_t *cofimp,
1075  cs_real_t *dplus);
1076 
1077 /*-------------------------------------------------------------------------------*/
1078 
1092 /*-------------------------------------------------------------------------------
1093  Arguments
1094  ______________________________________________________________________________.
1095  mode name role !
1096  ______________________________________________________________________________*/
1107 /*-------------------------------------------------------------------------------*/
1108 
1109 void
1111  cs_real_t prl,
1112  cs_real_t prt,
1113  cs_real_t yplus,
1114  cs_real_t dplus,
1115  cs_real_t *htur,
1116  cs_real_t *yplim);
1117 
1118 /*----------------------------------------------------------------------------*/
1119 
1121 
1122 #endif /* __CS_WALL_FUNCTIONS_H__ */
const double cs_turb_cstlog_rough
Definition: cs_turbulence_model.c:326
const double cs_turb_xkappa
Definition: cs_turbulence_model.c:295
double precision epzero
epsilon
Definition: cstnum.f90:40
cs_wall_f_type_t iwallf
Definition: cs_wall_functions.h:82
void wallfunctions(const cs_int_t *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 roughness, const cs_real_t *const rnnb, const cs_real_t *const kinetic_en, cs_int_t *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:185
static void cs_wall_functions_2scales_scalable(cs_real_t l_visc, cs_real_t t_visc, cs_real_t vel, cs_real_t y, cs_real_t kinetic_en, int *iuntur, cs_lnum_t *nsubla, cs_lnum_t *nlogla, cs_real_t *ustar, cs_real_t *uk, cs_real_t *yplus, cs_real_t *dplus, cs_real_t *ypup, cs_real_t *cofimp)
Scalable wall function: shift the wall if .
Definition: cs_wall_functions.h:386
#define BEGIN_C_DECLS
Definition: cs_defs.h:429
integer(c_int), pointer, save iwalfs
Wall functions for scalar.
Definition: optcal.f90:552
int cs_int_t
Fortran-compatible integer.
Definition: cs_defs.h:295
double cs_turb_cmu025
Definition: cs_turbulence_model.c:356
static cs_real_t _vdriest_dupdyp_integral(cs_real_t yplus)
Definition: cs_wall_functions.h:460
cs_wall_functions_t * cs_get_glob_wall_functions(void)
Definition: cs_wall_functions.c:263
Definition: cs_wall_functions.h:72
Definition: cs_wall_functions.h:61
double precision, dimension(ncharm), save a2
Definition: cpincl.f90:226
void hturbp(const cs_int_t *const iwalfs, const cs_real_t *const prl, const cs_real_t *const prt, 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:233
Definition: cs_wall_functions.h:60
double cs_turb_dpow
Definition: cs_turbulence_model.c:345
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 roughness, 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:303
const double cs_turb_bpow
Definition: cs_turbulence_model.c:342
Definition: cs_wall_functions.h:62
static void cs_wall_functions_s_arpaci_larsen(cs_real_t prl, cs_real_t prt, 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:829
cs_wall_f_s_type_t
Definition: cs_wall_functions.h:70
Definition: cs_wall_functions.h:64
const double cs_turb_vdriest
Definition: cs_turbulence_model.c:304
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: with ...
Definition: cs_wall_functions.h:918
double cs_turb_cstlog_alpha
Definition: cs_turbulence_model.c:336
#define CS_MIN(a, b)
Definition: cs_defs.h:402
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:265
integer(c_int), pointer, save iwallf
Wall functions.
Definition: optcal.f90:547
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 roughness, cs_real_t kinetic_en, int *iuntur, cs_lnum_t *nsubla, cs_lnum_t *nlogla, cs_real_t *ustar, cs_real_t *uk, cs_real_t *yplus, cs_real_t *dplus, cs_real_t *ypup, cs_real_t *cofimp)
Two velocity scales wall function with automatic switch from rough to smooth.
Definition: cs_wall_functions.h:664
Definition: cs_wall_functions.h:63
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:292
static void cs_wall_functions_1scale_power(cs_real_t l_visc, cs_real_t vel, cs_real_t y, int *iuntur, cs_lnum_t *nsubla, cs_lnum_t *nlogla, cs_real_t *ustar, cs_real_t *uk, cs_real_t *yplus, cs_real_t *ypup, cs_real_t *cofimp)
Power law: Werner & Wengle.
Definition: cs_wall_functions.h:128
static void cs_wall_functions_2scales_vdriest(cs_real_t rnnb, cs_real_t l_visc, cs_real_t vel, cs_real_t y, cs_real_t kinetic_en, int *iuntur, cs_lnum_t *nsubla, cs_lnum_t *nlogla, cs_real_t *ustar, cs_real_t *uk, cs_real_t *yplus, cs_real_t *ypup, cs_real_t *cofimp, cs_real_t *lmk, cs_real_t kr, bool wf)
Two velocity scales wall function using Van Driest mixing length.
Definition: cs_wall_functions.h:534
#define END_C_DECLS
Definition: cs_defs.h:430
wall functions descriptor.
Definition: cs_wall_functions.h:80
#define _(String)
Definition: cs_defs.h:52
double cs_real_t
Definition: cs_defs.h:296
Definition: cs_wall_functions.h:66
const double cs_turb_cstlog
Definition: cs_turbulence_model.c:315
#define CS_PROCF(x, y)
Definition: cs_defs.h:453
#define CS_MAX(a, b)
Definition: cs_defs.h:403
int iwallt
Definition: cs_wall_functions.h:86
static void cs_wall_functions_2scales_log(cs_real_t l_visc, cs_real_t t_visc, cs_real_t vel, cs_real_t y, cs_real_t kinetic_en, int *iuntur, cs_lnum_t *nsubla, cs_lnum_t *nlogla, cs_real_t *ustar, cs_real_t *uk, cs_real_t *yplus, cs_real_t *ypup, cs_real_t *cofimp)
Log law: piecewise linear and log, with two velocity scales based on the friction and the turbulent k...
Definition: cs_wall_functions.h:303
static void cs_wall_functions_disabled(cs_real_t l_visc, cs_real_t t_visc, cs_real_t vel, cs_real_t y, int *iuntur, cs_lnum_t *nsubla, cs_lnum_t *nlogla, cs_real_t *ustar, cs_real_t *uk, cs_real_t *yplus, cs_real_t *dplus, cs_real_t *ypup, cs_real_t *cofimp)
No wall function.
Definition: cs_wall_functions.h:762
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:84
Definition: cs_field_pointer.h:70
double ypluli
Definition: cs_wall_functions.h:91
cs_wall_f_type_t
Definition: cs_wall_functions.h:58
Definition: cs_wall_functions.h:73
const double cs_turb_crij1
Definition: cs_turbulence_model.c:400
const double cs_turb_crij2
Definition: cs_turbulence_model.c:406
static void cs_wall_functions_1scale_log(cs_lnum_t ifac, cs_real_t l_visc, cs_real_t vel, cs_real_t y, int *iuntur, cs_lnum_t *nsubla, cs_lnum_t *nlogla, cs_real_t *ustar, cs_real_t *uk, cs_real_t *yplus, cs_real_t *ypup, cs_real_t *cofimp)
Log law: piecewise linear and log, with one velocity scale based on the friction. ...
Definition: cs_wall_functions.h:199
const cs_wall_functions_t * cs_glob_wall_functions
void cs_wall_functions_scalar(cs_wall_f_s_type_t iwalfs, cs_real_t prl, cs_real_t prt, 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:478
Definition: cs_wall_functions.h:65
const double cs_turb_apow
Definition: cs_turbulence_model.c:339