applications/solvers/multiphase/settlingFoam/wallFunctions.H

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