applications/solvers/multiphase/twoPhaseEulerFoam/turbulenceModel/wallFunctions.H

   1 {
   2     labelList cellBoundaryFaceCount(epsilon.size(), 0);
   3
   4     scalar Cmu25 = ::pow(Cmu, 0.25);
   5     scalar Cmu75 = ::pow(Cmu, 0.75);
   6
   7     const fvPatchList& patches = mesh.boundary();
   8
   9     //- Initialise the near-wall P field to zero
  10     forAll(patches, patchi)
  11     {
  12         const fvPatch& currPatch = patches[patchi];
  13
  14         if (isType<wallFvPatch>(currPatch))
  15         {
  16             forAll(currPatch, facei)
  17             {
  18                 label faceCelli = currPatch.faceCells()[facei];
  19
  20                 epsilon[faceCelli] = 0.0;
  21                 G[faceCelli] = 0.0;
  22             }
  23         }
  24     }
  25
  26     //- Accumulate the wall face contributions to epsilon and G
  27     //  Increment cellBoundaryFaceCount for each face for averaging
  28     forAll(patches, patchi)
  29     {
  30         const fvPatch& currPatch = patches[patchi];
  31
  32         if (isType<wallFvPatch>(currPatch))
  33         {
  34             const scalarField& nuw = nutb.boundaryField()[patchi];
  35
  36             scalarField magFaceGradU = mag(U.boundaryField()[patchi].snGrad());
  37
  38             forAll(currPatch, facei)
  39             {
  40                 label faceCelli = currPatch.faceCells()[facei];
  41
  42                 scalar yPlus =
  43                     Cmu25*y[patchi][facei]
  44                     *::sqrt(k[faceCelli])
  45                     /nub.value();
  46
  47                 // For corner cells (with two boundary or more faces),
  48                 // epsilon and G in the near-wall cell are calculated
  49                 // as an average
  50
  51                 cellBoundaryFaceCount[faceCelli]++;
  52
  53                 epsilon[faceCelli] +=
  54                      Cmu75*::pow(k[faceCelli], 1.5)
  55                     /(kappa*y[patchi][facei]);
  56
  57                 if (yPlus > 11.6)
  58                 {
  59                     G[faceCelli] +=
  60                         nuw[facei]*magFaceGradU[facei]
  61                         *Cmu25*::sqrt(k[faceCelli])
  62                         /(kappa*y[patchi][facei]);
  63                 }
  64             }
  65         }
  66     }
  67
  68
  69     // perform the averaging
  70
  71     forAll(patches, patchi)
  72     {
  73         const fvPatch& curPatch = patches[patchi];
  74
  75         if (isType<wallFvPatch>(curPatch))
  76         {
  77             forAll(curPatch, facei)
  78             {
  79                 label faceCelli = curPatch.faceCells()[facei];
  80
  81                 epsilon[faceCelli] /= cellBoundaryFaceCount[faceCelli];
  82                 G[faceCelli] /= cellBoundaryFaceCount[faceCelli];
  83             }
  84         }
  85     }
  86 }