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 (isA<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 (isA<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             (
  42                 mag(U.boundaryField()[patchi].snGrad())
  43             );
  44
  45             forAll(curPatch, facei)
  46             {
  47                 label faceCelli = curPatch.faceCells()[facei];
  48
  49                 scalar yPlus =
  50                     Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])
  51                    /(muw[facei]/rhow[facei]);
  52
  53                 // For corner cells (with two boundary or more faces),
  54                 // epsilon and G in the near-wall cell are calculated
  55                 // as an average
  56
  57                 cellBoundaryFaceCount[faceCelli]++;
  58
  59                 epsilon[faceCelli] +=
  60                     Cmu75*rho[faceCelli]*::pow(k[faceCelli], 1.5)
  61                    /(kappa_*y[patchi][facei]);
  62
  63                 if (yPlus > 11.6)
  64                 {
  65                     G[faceCelli] +=
  66                         mutw[facei]*magFaceGradU[facei]
  67                        *Cmu25*::sqrt(k[faceCelli])
  68                        /(kappa_*y[patchi][facei]);
  69                 }
  70             }
  71         }
  72     }
  73
  74
  75     // perform the averaging
  76
  77     forAll(patches, patchi)
  78     {
  79         const fvPatch& curPatch = patches[patchi];
  80
  81         if (isA<wallFvPatch>(curPatch))
  82         {
  83             forAll(curPatch, facei)
  84             {
  85                 label faceCelli = curPatch.faceCells()[facei];
  86
  87                 epsilon[faceCelli] /= cellBoundaryFaceCount[faceCelli];
  88                 G[faceCelli] /= cellBoundaryFaceCount[faceCelli];
  89             }
  90         }
  91     }
  92 }