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