Skip to content
Permalink
Browse files

Rationalised wall function implementation to avoid complex and incons…

…istent coefficients

All wall functions now operate collaboratively, obtaining the Cmu, kappa and E
coefficients and yPlusLam from the nutWallFunction base class.  Now these
optional inputs need only be specified in the nut boundary condition with the k,
epsilon, omega, v2 and f wall functions obtaining these values from there.  This
is much simpler to specify and avoids inconsistencies in the operation of the
wall functions for the different turbulence fields.

The code has also been rationalised and simplified avoiding unnecessary code
and duplication.
  • Loading branch information...
Henry Weller
Henry Weller committed May 31, 2019
1 parent fc4d7b9 commit 1e2550b6cd5eca51b7ffbe89dd72ab7d8a4426dc
Showing with 284 additions and 895 deletions.
  1. +1 −3 ...alphatFixedDmdtWallBoilingWallFunction/alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField.C
  2. +39 −56 ...atPhaseChangeJayatillekeWallFunction/alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField.C
  3. +2 −12 ...atPhaseChangeJayatillekeWallFunction/alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField.H
  4. +25 −13 ...ivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
  5. +36 −61 ...phatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C
  6. +2 −12 ...phatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.H
  7. +29 −66 ...phatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C
  8. +2 −3 ...phatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.H
  9. +19 −80 ...ds/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
  10. +2 −32 ...ds/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H
  11. +20 −102 ...derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.C
  12. +2 −45 ...derivedFvPatchFields/wallFunctions/fWallFunctions/fWallFunction/fWallFunctionFvPatchScalarField.H
  13. +21 −91 ...chFields/wallFunctions/kqRWallFunctions/kLowReWallFunction/kLowReWallFunctionFvPatchScalarField.C
  14. +2 −28 ...chFields/wallFunctions/kqRWallFunctions/kLowReWallFunction/kLowReWallFunctionFvPatchScalarField.H
  15. +5 −35 ...derivedFvPatchFields/wallFunctions/kqRWallFunctions/kqRWallFunction/kqRWallFunctionFvPatchField.C
  16. +0 −6 ...derivedFvPatchFields/wallFunctions/kqRWallFunctions/kqRWallFunction/kqRWallFunctionFvPatchField.H
  17. +1 −0 ...elds/wallFunctions/nutWallFunctions/nutLowReWallFunction/nutLowReWallFunctionFvPatchScalarField.H
  18. +1 −0 ...ds/wallFunctions/nutWallFunctions/nutURoughWallFunction/nutURoughWallFunctionFvPatchScalarField.H
  19. +1 −0 ...lFunctions/nutWallFunctions/nutUSpaldingWallFunction/nutUSpaldingWallFunctionFvPatchScalarField.H
  20. +1 −0 ...unctions/nutWallFunctions/nutUTabulatedWallFunction/nutUTabulatedWallFunctionFvPatchScalarField.H
  21. +1 −0 ...vPatchFields/wallFunctions/nutWallFunctions/nutUWallFunction/nutUWallFunctionFvPatchScalarField.H
  22. +19 −1 ...dFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.H
  23. +1 −0 ...ds/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.H
  24. +1 −0 ...vPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.H
  25. +21 −74 ...chFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
  26. +6 −26 ...chFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H
  27. +21 −103 ...ivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.C
  28. +2 −45 ...ivedFvPatchFields/wallFunctions/v2WallFunctions/v2WallFunction/v2WallFunctionFvPatchScalarField.H
  29. +1 −1 tutorials/incompressible/simpleFoam/T3A/0/nut
@@ -48,9 +48,7 @@ alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
relax_(1.0),
fixedDmdt_(0.0),
L_(0.0)
{
checkType();
}
{}


alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::
@@ -24,14 +24,11 @@ License
\*---------------------------------------------------------------------------*/

#include "alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField.H"
#include "fvPatchFieldMapper.H"
#include "addToRunTimeSelectionTable.H"

#include "phaseSystem.H"
#include "compressibleTurbulenceModel.H"
#include "ThermalDiffusivity.H"
#include "PhaseCompressibleTurbulenceModel.H"
#include "wallFvPatch.H"
#include "addToRunTimeSelectionTable.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@@ -51,18 +48,6 @@ label alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::maxIters_

// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //

void alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::checkType()
{
if (!isA<wallFvPatch>(patch()))
{
FatalErrorInFunction
<< "Patch type for patch " << patch().name() << " must be wall\n"
<< "Current patch type is " << patch().type() << nl
<< exit(FatalError);
}
}


tmp<scalarField>
alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::Psmooth
(
@@ -76,6 +61,7 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::Psmooth
tmp<scalarField>
alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::yPlusTherm
(
const nutWallFunctionFvPatchScalarField& nutw,
const scalarField& P,
const scalarField& Prat
) const
@@ -89,9 +75,10 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::yPlusTherm

for (int i=0; i<maxIters_; i++)
{
scalar f = ypt - (log(E_*ypt)/kappa_ + P[facei])/Prat[facei];
scalar df = 1 - 1.0/(ypt*kappa_*Prat[facei]);
scalar yptNew = ypt - f/df;
const scalar f =
ypt - (log(nutw.E()*ypt)/nutw.kappa() + P[facei])/Prat[facei];
const scalar df = 1 - 1.0/(ypt*nutw.kappa()*Prat[facei]);
const scalar yptNew = ypt - f/df;

if (yptNew < vSmall)
{
@@ -138,7 +125,13 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::calcAlphat
IOobject::groupName(turbulenceModel::propertiesName, phase.name())
);

const scalar Cmu25 = pow025(Cmu_);
const nutWallFunctionFvPatchScalarField& nutw =
refCast<const nutWallFunctionFvPatchScalarField>
(
turbModel.nut()().boundaryField()[patchi]
);

const scalar Cmu25 = pow025(nutw.Cmu());

const scalarField& y = turbModel.y()[patchi];

@@ -183,7 +176,7 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::calcAlphat
// Thermal sublayer thickness
scalarField P(this->Psmooth(Prat));

scalarField yPlusTherm(this->yPlusTherm(P, Prat));
scalarField yPlusTherm(this->yPlusTherm(nutw, P, Prat));

tmp<scalarField> talphatConv(new scalarField(this->size()));
scalarField& alphatConv = talphatConv.ref();
@@ -195,21 +188,31 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::calcAlphat
scalar alphaEff = 0.0;
if (yPlus[facei] < yPlusTherm[facei])
{
scalar A = qDot[facei]*rhow[facei]*uTau[facei]*y[facei];
scalar B = qDot[facei]*Pr[facei]*yPlus[facei];
scalar C = Pr[facei]*0.5*rhow[facei]*uTau[facei]*sqr(magUp[facei]);
const scalar A = qDot[facei]*rhow[facei]*uTau[facei]*y[facei];

const scalar B = qDot[facei]*Pr[facei]*yPlus[facei];

const scalar C =
Pr[facei]*0.5*rhow[facei]*uTau[facei]*sqr(magUp[facei]);

alphaEff = A/(B + C + vSmall);
}
else
{
scalar A = qDot[facei]*rhow[facei]*uTau[facei]*y[facei];
scalar B =
qDot[facei]*Prt_*(1.0/kappa_*log(E_*yPlus[facei]) + P[facei]);
scalar magUc =
uTau[facei]/kappa_*log(E_*yPlusTherm[facei]) - mag(Uw[facei]);
scalar C =
const scalar A = qDot[facei]*rhow[facei]*uTau[facei]*y[facei];

const scalar B =
qDot[facei]*Prt_
*(1.0/nutw.kappa()*log(nutw.E()*yPlus[facei]) + P[facei]);

const scalar magUc =
uTau[facei]/nutw.kappa()
*log(nutw.E()*yPlusTherm[facei]) - mag(Uw[facei]);

const scalar C =
0.5*rhow[facei]*uTau[facei]
*(Prt_*sqr(magUp[facei]) + (Pr[facei] - Prt_)*sqr(magUc));

alphaEff = A/(B + C + vSmall);
}

@@ -231,13 +234,8 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
)
:
alphatPhaseChangeWallFunctionFvPatchScalarField(p, iF),
Prt_(0.85),
Cmu_(0.09),
kappa_(0.41),
E_(9.8)
{
checkType();
}
Prt_(0.85)
{}


alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::
@@ -249,10 +247,7 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
)
:
alphatPhaseChangeWallFunctionFvPatchScalarField(p, iF, dict),
Prt_(dict.lookupOrDefault<scalar>("Prt", 0.85)),
Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
E_(dict.lookupOrDefault<scalar>("E", 9.8))
Prt_(dict.lookupOrDefault<scalar>("Prt", 0.85))
{}


@@ -266,10 +261,7 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
)
:
alphatPhaseChangeWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
Prt_(ptf.Prt_),
Cmu_(ptf.Cmu_),
kappa_(ptf.kappa_),
E_(ptf.E_)
Prt_(ptf.Prt_)
{}


@@ -280,10 +272,7 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
)
:
alphatPhaseChangeWallFunctionFvPatchScalarField(awfpsf),
Prt_(awfpsf.Prt_),
Cmu_(awfpsf.Cmu_),
kappa_(awfpsf.kappa_),
E_(awfpsf.E_)
Prt_(awfpsf.Prt_)
{}


@@ -295,10 +284,7 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
)
:
alphatPhaseChangeWallFunctionFvPatchScalarField(awfpsf, iF),
Prt_(awfpsf.Prt_),
Cmu_(awfpsf.Cmu_),
kappa_(awfpsf.kappa_),
E_(awfpsf.E_)
Prt_(awfpsf.Prt_)
{}


@@ -324,9 +310,6 @@ void alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::write
{
fvPatchField<scalar>::write(os);
writeEntry(os, "Prt", Prt_);
writeEntry(os, "Cmu", Cmu_);
writeEntry(os, "kappa", kappa_);
writeEntry(os, "E", E_);
writeEntry(os, "dmdt", dmdt_);
writeEntry(os, "value", *this);
}
@@ -63,6 +63,7 @@ SourceFiles
#define compressible_alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField_H

#include "alphatPhaseChangeWallFunctionFvPatchScalarField.H"
#include "nutWallFunctionFvPatchScalarField.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@@ -87,15 +88,6 @@ protected:
//- Turbulent Prandtl number
scalar Prt_;

//- Cmu coefficient
scalar Cmu_;

//- Von Karman constant
scalar kappa_;

//- E coefficient
scalar E_;

// Solution parameters

static scalar maxExp_;
@@ -105,15 +97,13 @@ protected:

// Protected Member Functions

//- Check the type of the patch
void checkType();

//- 'P' function
tmp<scalarField> Psmooth(const scalarField& Prat) const;

//- Calculate y+ at the edge of the thermal laminar sublayer
tmp<scalarField> yPlusTherm
(
const nutWallFunctionFvPatchScalarField& nutw,
const scalarField& P,
const scalarField& Prat
) const;
@@ -24,17 +24,12 @@ License
\*---------------------------------------------------------------------------*/

#include "alphatWallBoilingWallFunctionFvPatchScalarField.H"
#include "fvPatchFieldMapper.H"
#include "addToRunTimeSelectionTable.H"

#include "phaseSystem.H"
#include "compressibleTurbulenceModel.H"
#include "ThermalDiffusivity.H"
#include "PhaseCompressibleTurbulenceModel.H"
#include "saturationModel.H"
#include "wallFvPatch.H"
#include "uniformDimensionedFields.H"
#include "mathematicalConstants.H"
#include "addToRunTimeSelectionTable.H"

using namespace Foam::constant::mathematical;

@@ -428,9 +423,13 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
)
);

const tmp<scalarField> tnutw = turbModel.nut(patchi);
const nutWallFunctionFvPatchScalarField& nutw =
refCast<const nutWallFunctionFvPatchScalarField>
(
turbModel.nut()().boundaryField()[patchi]
);

const scalar Cmu25(pow025(Cmu_));
const scalar Cmu25(pow025(nutw.Cmu()));

const scalarField& y = turbModel.y()[patchi];

@@ -474,7 +473,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
// Thermal sublayer thickness
const scalarField P(this->Psmooth(Prat));

const scalarField yPlusTherm(this->yPlusTherm(P, Prat));
const scalarField yPlusTherm(this->yPlusTherm(nutw, P, Prat));

tmp<volScalarField> tCp = liquid.thermo().Cp();
const volScalarField& Cp = tCp();
@@ -517,11 +516,24 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
{
// Liquid temperature at y+=250 is estimated from logarithmic
// thermal wall function (Koncar, Krepper & Egorov, 2005)
const scalarField Tplus_y250(Prt_*(log(E_*250)/kappa_ + P));
const scalarField Tplus_y250
(
Prt_*(log(nutw.E()*250)/nutw.kappa() + P)
);

const scalarField Tplus(Prt_*(log(E_*yPlus)/kappa_ + P));
scalarField Tl(Tw - (Tplus_y250/Tplus)*(Tw - Tc));
Tl = max(Tc - 40, Tl);
const scalarField Tplus
(
Prt_*(log(nutw.E()*yPlus)/nutw.kappa() + P)
);

const scalarField Tl
(
max
(
Tc - 40,
Tw - (Tplus_y250/Tplus)*(Tw - Tc)
)
);

// Nucleation site density:
const scalarField N

0 comments on commit 1e2550b

Please sign in to comment.
You can’t perform that action at this time.