This repository has been archived on 2023-11-20. You can view files and clone it, but cannot push or open issues or pull requests.
foam-extend4.1-coherent-io/applications/solvers/multiphase/twoPhaseEulerFoam/turbulenceModel/wallFunctions.H

86 lines
2.3 KiB
C

{
labelList cellBoundaryFaceCount(epsilon.size(), 0);
scalar Cmu25 = ::pow(Cmu, 0.25);
scalar Cmu75 = ::pow(Cmu, 0.75);
const fvPatchList& patches = mesh.boundary();
//- Initialise the near-wall P field to zero
forAll(patches, patchi)
{
const fvPatch& currPatch = patches[patchi];
if (isType<wallFvPatch>(currPatch))
{
forAll(currPatch, facei)
{
label faceCelli = currPatch.faceCells()[facei];
epsilon[faceCelli] = 0.0;
G[faceCelli] = 0.0;
}
}
}
//- Accumulate the wall face contributions to epsilon and G
// Increment cellBoundaryFaceCount for each face for averaging
forAll(patches, patchi)
{
const fvPatch& currPatch = patches[patchi];
if (isType<wallFvPatch>(currPatch))
{
const scalarField& nuw = nutb.boundaryField()[patchi];
scalarField magFaceGradU = mag(U.boundaryField()[patchi].snGrad());
forAll(currPatch, facei)
{
label faceCelli = currPatch.faceCells()[facei];
scalar yPlus =
Cmu25*y[patchi][facei]
*::sqrt(k[faceCelli])
/nub.value();
// For corner cells (with two boundary or more faces),
// epsilon and G in the near-wall cell are calculated
// as an average
cellBoundaryFaceCount[faceCelli]++;
epsilon[faceCelli] +=
Cmu75*::pow(k[faceCelli], 1.5)
/(kappa*y[patchi][facei]);
if (yPlus > 11.6)
{
G[faceCelli] +=
nuw[facei]*magFaceGradU[facei]
*Cmu25*::sqrt(k[faceCelli])
/(kappa*y[patchi][facei]);
}
}
}
}
// perform the averaging
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
if (isType<wallFvPatch>(curPatch))
{
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
epsilon[faceCelli] /= cellBoundaryFaceCount[faceCelli];
G[faceCelli] /= cellBoundaryFaceCount[faceCelli];
}
}
}
}