From c7836f9d83bc12a6cf6d2f8377ac4de016aefaaa Mon Sep 17 00:00:00 2001 From: Kutalmis Bercin Date: Tue, 3 May 2022 13:57:21 +0100 Subject: [PATCH] BUG: nutWallFunction: prevent casting errors (fixes #1457, #1966) Previously, a nutWallFunctionFvPatchScalarField ref should be created in epsilon, k, and omega wall functions to fetch various common wall-function coefficients necessary to carry out and complete local operations inside these wall functions. However, this arrangement required the use of a nut wall function, even when unnecessary, when any of non-nut wall functions are being used. Therefore, some users had been redundantly restrained and obstructed with rather obscure casting-error messages. Also, the wall-function coefficients Cmu, kappa and E have been obtained from the specified nutWallFunction in order to ensure that each patch possesses the same set of values for these coefficients. Although the motivation sounds reasonable, it has also been putting redundant restraints on users and disregarding the specifics of each wall-function. For example, the variation of epsilon in near-wall regions is usually very steep and non-monotonic specific - an expert user may therefore want to use an epsilon-specific coefficient, and this was not allowed by the previous arrangement. This commit introduces a new class (i.e. wallFunctionCoefficients) comprising all common wall-function coefficients and yPlus calculations. --- .../turbulenceModels/Make/files | 2 + .../epsilonWallFunctionFvPatchScalarField.C | 25 ++-- .../epsilonWallFunctionFvPatchScalarField.H | 4 + .../kLowReWallFunctionFvPatchScalarField.C | 27 ++-- .../kLowReWallFunctionFvPatchScalarField.H | 4 + ...utUBlendedWallFunctionFvPatchScalarField.C | 5 +- .../nutURoughWallFunctionFvPatchScalarField.C | 20 ++- ...tUSpaldingWallFunctionFvPatchScalarField.C | 9 +- .../nutUWallFunctionFvPatchScalarField.C | 19 ++- .../nutWallFunctionFvPatchScalarField.C | 59 ++------ .../nutWallFunctionFvPatchScalarField.H | 39 +----- .../nutkRoughWallFunctionFvPatchScalarField.C | 8 +- .../nutkWallFunctionFvPatchScalarField.C | 11 +- .../omegaWallFunctionFvPatchScalarField.C | 23 +++- .../omegaWallFunctionFvPatchScalarField.H | 4 + .../wallFunctionCoefficients.C | 96 +++++++++++++ .../wallFunctionCoefficients.H | 129 ++++++++++++++++++ ...atmEpsilonWallFunctionFvPatchScalarField.C | 18 ++- .../atmNutUWallFunctionFvPatchScalarField.C | 4 +- .../atmNutWallFunctionFvPatchScalarField.C | 5 +- .../atmNutkWallFunctionFvPatchScalarField.C | 5 +- .../atmOmegaWallFunctionFvPatchScalarField.C | 13 +- .../nutkFilmWallFunctionFvPatchScalarField.C | 9 +- 23 files changed, 380 insertions(+), 158 deletions(-) create mode 100644 src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/wallFunction/wallFunctionCoefficients/wallFunctionCoefficients.C create mode 100644 src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/wallFunction/wallFunctionCoefficients/wallFunctionCoefficients.H diff --git a/src/TurbulenceModels/turbulenceModels/Make/files b/src/TurbulenceModels/turbulenceModels/Make/files index f542cc7395..52b7043056 100644 --- a/src/TurbulenceModels/turbulenceModels/Make/files +++ b/src/TurbulenceModels/turbulenceModels/Make/files @@ -32,6 +32,8 @@ derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchField.C /* Wall function BCs */ wallFunctions = derivedFvPatchFields/wallFunctions +$(wallFunctions)/wallFunction/wallFunctionCoefficients/wallFunctionCoefficients.C + nutWallFunctions = $(wallFunctions)/nutWallFunctions $(nutWallFunctions)/nutWallFunction/nutWallFunctionFvPatchScalarField.C diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C index d076876144..34af105367 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C @@ -202,8 +202,8 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculate { const label patchi = patch.index(); - const nutWallFunctionFvPatchScalarField& nutw = - nutWallFunctionFvPatchScalarField::nutw(turbModel, patchi); + const tmp tnutw = turbModel.nut(patchi); + const scalarField& nutw = tnutw(); const scalarField& y = turbModel.y()[patchi]; @@ -217,8 +217,10 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculate const scalarField magGradUw(mag(Uw.snGrad())); - const scalar Cmu25 = pow025(nutw.Cmu()); - const scalar Cmu75 = pow(nutw.Cmu(), 0.75); + const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); + const scalar Cmu75 = pow(wallCoeffs_.Cmu(), 0.75); + const scalar kappa = wallCoeffs_.kappa(); + const scalar yPlusLam = wallCoeffs_.yPlusLam(); // Set epsilon and G forAll(nutw, facei) @@ -236,13 +238,13 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculate // Contribution from the inertial sublayer const scalar epsilonLog = - w*Cmu75*pow(k[celli], 1.5)/(nutw.kappa()*y[facei]); + w*Cmu75*pow(k[celli], 1.5)/(kappa*y[facei]); switch (blending_) { case blendingType::STEPWISE: { - if (lowReCorrection_ && yPlus < nutw.yPlusLam()) + if (lowReCorrection_ && yPlus < yPlusLam) { epsilonBlended = epsilonVis; } @@ -285,14 +287,14 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculate epsilon0[celli] += epsilonBlended; - if (!(lowReCorrection_ && yPlus < nutw.yPlusLam())) + if (!(lowReCorrection_ && yPlus < yPlusLam)) { G0[celli] += w *(nutw[facei] + nuw[facei]) *magGradUw[facei] *Cmu25*sqrt(k[celli]) - /(nutw.kappa()*y[facei]); + /(kappa*y[facei]); } } } @@ -310,6 +312,8 @@ void Foam::epsilonWallFunctionFvPatchScalarField::writeLocalEntries { os.writeEntry("n", n_); } + + wallCoeffs_.writeEntries(os); } @@ -328,6 +332,7 @@ epsilonWallFunctionFvPatchScalarField lowReCorrection_(false), initialised_(false), master_(-1), + wallCoeffs_(), G_(), epsilon_(), cornerWeights_() @@ -349,6 +354,7 @@ epsilonWallFunctionFvPatchScalarField lowReCorrection_(ptf.lowReCorrection_), initialised_(false), master_(-1), + wallCoeffs_(ptf.wallCoeffs_), G_(), epsilon_(), cornerWeights_() @@ -385,6 +391,7 @@ epsilonWallFunctionFvPatchScalarField lowReCorrection_(dict.getOrDefault("lowReCorrection", false)), initialised_(false), master_(-1), + wallCoeffs_(dict), G_(), epsilon_(), cornerWeights_() @@ -406,6 +413,7 @@ epsilonWallFunctionFvPatchScalarField lowReCorrection_(ewfpsf.lowReCorrection_), initialised_(false), master_(-1), + wallCoeffs_(ewfpsf.wallCoeffs_), G_(), epsilon_(), cornerWeights_() @@ -425,6 +433,7 @@ epsilonWallFunctionFvPatchScalarField lowReCorrection_(ewfpsf.lowReCorrection_), initialised_(false), master_(-1), + wallCoeffs_(ewfpsf.wallCoeffs_), G_(), epsilon_(), cornerWeights_() diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H index 8baec73e06..cd09146965 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H @@ -173,6 +173,7 @@ SourceFiles #define epsilonWallFunctionFvPatchScalarField_H #include "fixedValueFvPatchField.H" +#include "wallFunctionCoefficients.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -230,6 +231,9 @@ protected: //- Master patch ID label master_; + //- Wall-function coefficients + wallFunctionCoefficients wallCoeffs_; + //- Local copy of turbulence G field scalarField G_; diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/kqRWallFunctions/kLowReWallFunction/kLowReWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/kqRWallFunctions/kLowReWallFunction/kLowReWallFunctionFvPatchScalarField.C index 51886bc3aa..56eca2ac2b 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/kqRWallFunctions/kLowReWallFunction/kLowReWallFunctionFvPatchScalarField.C +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/kqRWallFunctions/kLowReWallFunction/kLowReWallFunctionFvPatchScalarField.C @@ -42,6 +42,7 @@ void Foam::kLowReWallFunctionFvPatchScalarField::writeLocalEntries os.writeEntryIfDifferent("Ck", -0.416, Ck_); os.writeEntryIfDifferent("Bk", 8.366, Bk_); os.writeEntryIfDifferent("C", 11.0, C_); + wallCoeffs_.writeEntries(os); } @@ -57,7 +58,8 @@ Foam::kLowReWallFunctionFvPatchScalarField::kLowReWallFunctionFvPatchScalarField Ceps2_(1.9), Ck_(-0.416), Bk_(8.366), - C_(11.0) + C_(11.0), + wallCoeffs_() {} @@ -73,7 +75,8 @@ Foam::kLowReWallFunctionFvPatchScalarField::kLowReWallFunctionFvPatchScalarField Ceps2_(ptf.Ceps2_), Ck_(ptf.Ck_), Bk_(ptf.Bk_), - C_(ptf.C_) + C_(ptf.C_), + wallCoeffs_(ptf.wallCoeffs_) {} @@ -96,7 +99,8 @@ Foam::kLowReWallFunctionFvPatchScalarField::kLowReWallFunctionFvPatchScalarField ), Ck_(dict.getOrDefault("Ck", -0.416)), Bk_(dict.getOrDefault("Bk", 8.366)), - C_(dict.getOrDefault("C", 11.0)) + C_(dict.getOrDefault("C", 11.0)), + wallCoeffs_(dict) {} @@ -109,7 +113,8 @@ Foam::kLowReWallFunctionFvPatchScalarField::kLowReWallFunctionFvPatchScalarField Ceps2_(kwfpsf.Ceps2_), Ck_(kwfpsf.Ck_), Bk_(kwfpsf.Bk_), - C_(kwfpsf.C_) + C_(kwfpsf.C_), + wallCoeffs_(kwfpsf.wallCoeffs_) {} @@ -123,7 +128,8 @@ Foam::kLowReWallFunctionFvPatchScalarField::kLowReWallFunctionFvPatchScalarField Ceps2_(kwfpsf.Ceps2_), Ck_(kwfpsf.Ck_), Bk_(kwfpsf.Bk_), - C_(kwfpsf.C_) + C_(kwfpsf.C_), + wallCoeffs_(kwfpsf.wallCoeffs_) {} @@ -147,9 +153,6 @@ void Foam::kLowReWallFunctionFvPatchScalarField::updateCoeffs() ) ); - const nutWallFunctionFvPatchScalarField& nutw = - nutWallFunctionFvPatchScalarField::nutw(turbModel, patchi); - const scalarField& y = turbModel.y()[patchi]; const tmp tnuw = turbModel.nu(patchi); @@ -158,7 +161,9 @@ void Foam::kLowReWallFunctionFvPatchScalarField::updateCoeffs() const tmp tk = turbModel.k(); const volScalarField& k = tk(); - const scalar Cmu25 = pow025(nutw.Cmu()); + const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); + const scalar kappa = wallCoeffs_.kappa(); + const scalar yPlusLam = wallCoeffs_.yPlusLam(); scalarField& kw = *this; @@ -169,9 +174,9 @@ void Foam::kLowReWallFunctionFvPatchScalarField::updateCoeffs() const scalar uTau = Cmu25*sqrt(k[celli]); const scalar yPlus = uTau*y[facei]/nuw[facei]; - if (yPlus > nutw.yPlusLam()) + if (yPlus > yPlusLam) { - kw[facei] = Ck_/nutw.kappa()*log(yPlus) + Bk_; + kw[facei] = Ck_/kappa*log(yPlus) + Bk_; } else { diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/kqRWallFunctions/kLowReWallFunction/kLowReWallFunctionFvPatchScalarField.H b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/kqRWallFunctions/kLowReWallFunction/kLowReWallFunctionFvPatchScalarField.H index 951da88423..5b74b11dd1 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/kqRWallFunctions/kLowReWallFunction/kLowReWallFunctionFvPatchScalarField.H +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/kqRWallFunctions/kLowReWallFunction/kLowReWallFunctionFvPatchScalarField.H @@ -101,6 +101,7 @@ SourceFiles #define kLowReWallFunctionFvPatchScalarField_H #include "fixedValueFvPatchField.H" +#include "wallFunctionCoefficients.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -131,6 +132,9 @@ protected: //- C coefficient scalar C_; + //- Wall-function coefficients + wallFunctionCoefficients wallCoeffs_; + // Protected Member Functions diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUBlendedWallFunction/nutUBlendedWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUBlendedWallFunction/nutUBlendedWallFunctionFvPatchScalarField.C index c68c8cffac..405032a8e4 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUBlendedWallFunction/nutUBlendedWallFunctionFvPatchScalarField.C +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUBlendedWallFunction/nutUBlendedWallFunctionFvPatchScalarField.C @@ -92,6 +92,9 @@ Foam::nutUBlendedWallFunctionFvPatchScalarField::calcUTau const scalarField& nutw = *this; + const scalar kappa = wallCoeffs_.kappa(); + const scalar E = wallCoeffs_.E(); + forAll(uTaup, facei) { scalar ut = sqrt((nutw[facei] + nuw[facei])*magGradU[facei]); @@ -104,7 +107,7 @@ Foam::nutUBlendedWallFunctionFvPatchScalarField::calcUTau const scalar yPlus = y[facei]*ut/nuw[facei]; const scalar uTauVis = magUp[facei]/yPlus; const scalar uTauLog = - kappa_*magUp[facei]/log(max(E_*yPlus, 1 + 1e-4)); + kappa*magUp[facei]/log(max(E*yPlus, 1 + 1e-4)); const scalar utNew = pow(pow(uTauVis, n_) + pow(uTauLog, n_), 1.0/n_); diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutURoughWallFunction/nutURoughWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutURoughWallFunction/nutURoughWallFunctionFvPatchScalarField.C index a9c0c3f8b2..5810b81ca5 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutURoughWallFunction/nutURoughWallFunctionFvPatchScalarField.C +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutURoughWallFunction/nutURoughWallFunctionFvPatchScalarField.C @@ -56,6 +56,8 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcNut() const // The flow velocity at the adjacent cell centre const scalarField magUp(mag(Uw.patchInternalField() - Uw)); + const scalar yPlusLam = wallCoeffs_.yPlusLam(); + tmp tyPlus = calcYPlus(magUp); scalarField& yPlus = tyPlus.ref(); @@ -64,7 +66,7 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcNut() const forAll(yPlus, facei) { - if (yPlusLam_ < yPlus[facei]) + if (yPlusLam < yPlus[facei]) { const scalar Re = magUp[facei]*y[facei]/nuw[facei] + ROOTVSMALL; nutw[facei] = nuw[facei]*(sqr(yPlus[facei])/Re - 1); @@ -95,6 +97,10 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus const tmp tnuw = turbModel.nu(patchi); const scalarField& nuw = tnuw(); + const scalar kappa = wallCoeffs_.kappa(); + const scalar E = wallCoeffs_.E(); + const scalar yPlusLam = wallCoeffs_.yPlusLam(); + tmp tyPlus(new scalarField(patch().size(), Zero)); scalarField& yPlus = tyPlus.ref(); @@ -114,9 +120,9 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus { const scalar magUpara = magUp[facei]; const scalar Re = magUpara*y[facei]/nuw[facei]; - const scalar kappaRe = kappa_*Re; + const scalar kappaRe = kappa*Re; - scalar yp = yPlusLam_; + scalar yp = yPlusLam; const scalar ryPlusLam = 1.0/yp; int iter = 0; @@ -155,7 +161,7 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus (c_1*sint_2*KsPlus/t_1) + (c_3*logt_1*cos(t_2)); } - scalar denom = 1.0 + log(E_*yp) - G - yPlusGPrime; + scalar denom = 1.0 + log(E*yp) - G - yPlusGPrime; if (mag(denom) > VSMALL) { yp = (kappaRe + yp*(1 - yPlusGPrime))/denom; @@ -178,9 +184,9 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus { const scalar magUpara = magUp[facei]; const scalar Re = magUpara*y[facei]/nuw[facei]; - const scalar kappaRe = kappa_*Re; + const scalar kappaRe = kappa*Re; - scalar yp = yPlusLam_; + scalar yp = yPlusLam; const scalar ryPlusLam = 1.0/yp; int iter = 0; @@ -189,7 +195,7 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus do { yPlusLast = yp; - yp = (kappaRe + yp)/(1.0 + log(E_*yp)); + yp = (kappaRe + yp)/(1.0 + log(E*yp)); } while diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUSpaldingWallFunction/nutUSpaldingWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUSpaldingWallFunction/nutUSpaldingWallFunctionFvPatchScalarField.C index fe6ae4103a..33fac7aab0 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUSpaldingWallFunction/nutUSpaldingWallFunctionFvPatchScalarField.C +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUSpaldingWallFunction/nutUSpaldingWallFunctionFvPatchScalarField.C @@ -126,6 +126,9 @@ Foam::nutUSpaldingWallFunctionFvPatchScalarField::calcUTau const tmp tnuw = turbModel.nu(patchi); const scalarField& nuw = tnuw(); + const scalar kappa = wallCoeffs_.kappa(); + const scalar E = wallCoeffs_.E(); + const scalarField& nutw = *this; tmp tuTau(new scalarField(patch().size(), Zero)); @@ -146,18 +149,18 @@ Foam::nutUSpaldingWallFunctionFvPatchScalarField::calcUTau do { - scalar kUu = min(kappa_*magUp[facei]/ut, 50); + scalar kUu = min(kappa*magUp[facei]/ut, 50); scalar fkUu = exp(kUu) - 1 - kUu*(1 + 0.5*kUu); scalar f = - ut*y[facei]/nuw[facei] + magUp[facei]/ut - + 1/E_*(fkUu - 1.0/6.0*kUu*sqr(kUu)); + + 1/E*(fkUu - 1.0/6.0*kUu*sqr(kUu)); scalar df = y[facei]/nuw[facei] + magUp[facei]/sqr(ut) - + 1/E_*kUu*fkUu/ut; + + 1/E*kUu*fkUu/ut; scalar uTauNew = ut + f/df; err[facei] = mag((ut - uTauNew)/ut); diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUWallFunction/nutUWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUWallFunction/nutUWallFunctionFvPatchScalarField.C index 4c1923bb1c..c18d5df262 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUWallFunction/nutUWallFunctionFvPatchScalarField.C +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUWallFunction/nutUWallFunctionFvPatchScalarField.C @@ -53,6 +53,9 @@ Foam::nutUWallFunctionFvPatchScalarField::calcNut() const const tmp tnuw = turbModel.nu(patchi); const scalarField& nuw = tnuw(); + const scalar kappa = wallCoeffs_.kappa(); + const scalar E = wallCoeffs_.E(); + tmp tyPlus = calcYPlus(magUp); const scalarField& yPlus = tyPlus(); @@ -67,7 +70,7 @@ Foam::nutUWallFunctionFvPatchScalarField::calcNut() const // Inertial sublayer contribution const scalar nutLog = nuw[facei] - *(yPlus[facei]*kappa_/log(max(E_*yPlus[facei], 1 + 1e-4)) - 1.0); + *(yPlus[facei]*kappa/log(max(E*yPlus[facei], 1 + 1e-4)) - 1.0); nutw[facei] = blend(nutVis, nutLog, yPlus[facei]); } @@ -96,14 +99,18 @@ Foam::nutUWallFunctionFvPatchScalarField::calcYPlus const tmp tnuw = turbModel.nu(patchi); const scalarField& nuw = tnuw(); + const scalar kappa = wallCoeffs_.kappa(); + const scalar E = wallCoeffs_.E(); + const scalar yPlusLam = wallCoeffs_.yPlusLam(); + tmp tyPlus(new scalarField(patch().size(), Zero)); scalarField& yPlus = tyPlus.ref(); forAll(yPlus, facei) { - const scalar kappaRe = kappa_*magUp[facei]*y[facei]/nuw[facei]; + const scalar kappaRe = kappa*magUp[facei]*y[facei]/nuw[facei]; - scalar yp = yPlusLam_; + scalar yp = yPlusLam; const scalar ryPlusLam = 1.0/yp; int iter = 0; @@ -112,7 +119,7 @@ Foam::nutUWallFunctionFvPatchScalarField::calcYPlus do { yPlusLast = yp; - yp = (kappaRe + yp)/(1.0 + log(E_*yp)); + yp = (kappaRe + yp)/(1.0 + log(E*yp)); } while (mag(ryPlusLam*(yp - yPlusLast)) > 0.01 && ++iter < 10 ); @@ -219,12 +226,14 @@ Foam::nutUWallFunctionFvPatchScalarField::yPlus() const const scalarField magUp(mag(Uw.patchInternalField() - Uw)); const scalarField magGradUw(mag(Uw.snGrad())); + const scalar yPlusLam = wallCoeffs_.yPlusLam(); + tmp tyPlus = calcYPlus(magUp); scalarField& yPlus = tyPlus.ref(); forAll(yPlus, facei) { - if (yPlusLam_ > yPlus[facei]) + if (yPlusLam > yPlus[facei]) { // viscous sublayer yPlus[facei] = diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.C index 7013d1c191..b2678c966c 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.C +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016, 2019 OpenFOAM Foundation - Copyright (C) 2019-2020 OpenCFD Ltd. + Copyright (C) 2019-2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -92,10 +92,8 @@ void Foam::nutWallFunctionFvPatchScalarField::writeLocalEntries Ostream& os ) const { - os.writeEntry("Cmu", Cmu_); - os.writeEntry("kappa", kappa_); - os.writeEntry("E", E_); os.writeEntryIfDifferent("U", word::null, UName_); + wallCoeffs_.writeEntries(os); } @@ -111,10 +109,7 @@ Foam::nutWallFunctionFvPatchScalarField::nutWallFunctionFvPatchScalarField blending_(blendingType::STEPWISE), n_(4.0), UName_(word::null), - Cmu_(0.09), - kappa_(0.41), - E_(9.8), - yPlusLam_(yPlusLam(kappa_, E_)) + wallCoeffs_() { checkType(); } @@ -132,10 +127,7 @@ Foam::nutWallFunctionFvPatchScalarField::nutWallFunctionFvPatchScalarField blending_(ptf.blending_), n_(ptf.n_), UName_(ptf.UName_), - Cmu_(ptf.Cmu_), - kappa_(ptf.kappa_), - E_(ptf.E_), - yPlusLam_(ptf.yPlusLam_) + wallCoeffs_(ptf.wallCoeffs_) { checkType(); } @@ -168,13 +160,7 @@ Foam::nutWallFunctionFvPatchScalarField::nutWallFunctionFvPatchScalarField ) ), UName_(dict.getOrDefault("U", word::null)), - Cmu_(dict.getOrDefault("Cmu", 0.09)), - kappa_ - ( - dict.getCheckOrDefault("kappa", 0.41, scalarMinMax::ge(SMALL)) - ), - E_(dict.getCheckOrDefault("E", 9.8, scalarMinMax::ge(SMALL))), - yPlusLam_(yPlusLam(kappa_, E_)) + wallCoeffs_(dict) { checkType(); } @@ -189,10 +175,7 @@ Foam::nutWallFunctionFvPatchScalarField::nutWallFunctionFvPatchScalarField blending_(wfpsf.blending_), n_(wfpsf.n_), UName_(wfpsf.UName_), - Cmu_(wfpsf.Cmu_), - kappa_(wfpsf.kappa_), - E_(wfpsf.E_), - yPlusLam_(wfpsf.yPlusLam_) + wallCoeffs_(wfpsf.wallCoeffs_) { checkType(); } @@ -208,10 +191,7 @@ Foam::nutWallFunctionFvPatchScalarField::nutWallFunctionFvPatchScalarField blending_(wfpsf.blending_), n_(wfpsf.n_), UName_(wfpsf.UName_), - Cmu_(wfpsf.Cmu_), - kappa_(wfpsf.kappa_), - E_(wfpsf.E_), - yPlusLam_(wfpsf.yPlusLam_) + wallCoeffs_(wfpsf.wallCoeffs_) { checkType(); } @@ -235,23 +215,6 @@ Foam::nutWallFunctionFvPatchScalarField::nutw } -Foam::scalar Foam::nutWallFunctionFvPatchScalarField::yPlusLam -( - const scalar kappa, - const scalar E -) -{ - scalar ypl = 11.0; - - for (label i = 0; i < 10; ++i) - { - ypl = log(max(E*ypl, 1.0))/kappa; - } - - return ypl; -} - - Foam::scalar Foam::nutWallFunctionFvPatchScalarField::blend ( const scalar nutVis, @@ -265,7 +228,7 @@ Foam::scalar Foam::nutWallFunctionFvPatchScalarField::blend { case blendingType::STEPWISE: { - if (yPlus > yPlusLam_) + if (yPlus > wallCoeffs_.yPlusLam()) { nutw = nutLog; } @@ -310,12 +273,6 @@ Foam::scalar Foam::nutWallFunctionFvPatchScalarField::blend } -Foam::scalar Foam::nutWallFunctionFvPatchScalarField::yPlusLam() const -{ - return yPlusLam_; -} - - void Foam::nutWallFunctionFvPatchScalarField::updateCoeffs() { if (updated()) diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.H b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.H index 2c9713049b..4c697461e8 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.H +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.H @@ -167,6 +167,7 @@ SourceFiles #define nutWallFunctionFvPatchScalarField_H #include "fixedValueFvPatchFields.H" +#include "wallFunctionCoefficients.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -214,18 +215,8 @@ protected: // retrieved from the turbulence model word UName_; - //- Empirical model coefficient - scalar Cmu_; - - //- von Kármán constant - scalar kappa_; - - //- Wall roughness parameter - scalar E_; - - //- Estimated y+ value at the intersection - //- of the viscous and inertial sublayers - scalar yPlusLam_; + //- Wall-function coefficients + wallFunctionCoefficients wallCoeffs_; // Protected Member Functions @@ -296,24 +287,6 @@ public: // Member Functions - //- Return Cmu - scalar Cmu() const - { - return Cmu_; - } - - //- Return kappa - scalar kappa() const - { - return kappa_; - } - - //- Return E - scalar E() const - { - return E_; - } - //- Return the nut patchField for the given wall patch static const nutWallFunctionFvPatchScalarField& nutw ( @@ -321,12 +294,6 @@ public: const label patchi ); - //- Estimate the y+ at the intersection of the two sublayers - static scalar yPlusLam(const scalar kappa, const scalar E); - - //- Return the estimated y+ at the two-sublayer intersection - scalar yPlusLam() const; - //- Return the blended nut according to the chosen blending treatment scalar blend ( diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.C index c51e798168..85089d2d5b 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.C +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.C @@ -76,7 +76,9 @@ calcNut() const const tmp tnuw = turbModel.nu(patchi); const scalarField& nuw = tnuw(); - const scalar Cmu25 = pow025(Cmu_); + const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); + const scalar kappa = wallCoeffs_.kappa(); + const scalar E = wallCoeffs_.E(); tmp tnutw(new scalarField(*this)); scalarField& nutw = tnutw.ref(); @@ -89,7 +91,7 @@ calcNut() const const scalar yPlus = uStar*y[facei]/nuw[facei]; const scalar KsPlus = uStar*Ks_[facei]/nuw[facei]; - scalar Edash = E_; + scalar Edash = E; if (2.25 < KsPlus) { Edash /= fnRough(KsPlus, Cs_[facei]); @@ -105,7 +107,7 @@ calcNut() const min ( nuw[facei] - *(yPlus*kappa_/log(max(Edash*yPlus, 1+1e-4)) - 1), + *(yPlus*kappa/log(max(Edash*yPlus, 1+1e-4)) - 1), 2*limitingNutw ), 0.5*limitingNutw ); diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C index 5f5311adb5..5e0c841f6e 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C @@ -55,7 +55,9 @@ calcNut() const const tmp tnuw = turbModel.nu(patchi); const scalarField& nuw = tnuw(); - const scalar Cmu25 = pow025(Cmu_); + const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); + const scalar kappa = wallCoeffs_.kappa(); + const scalar E = wallCoeffs_.E(); tmp tnutw(new scalarField(patch().size(), Zero)); scalarField& nutw = tnutw.ref(); @@ -71,7 +73,7 @@ calcNut() const // Inertial sublayer contribution const scalar nutLog = - nuw[facei]*(yPlus*kappa_/log(max(E_*yPlus, 1 + 1e-4)) - 1.0); + nuw[facei]*(yPlus*kappa/log(max(E*yPlus, 1 + 1e-4)) - 1.0); nutw[facei] = blend(nutVis, nutLog, yPlus); } @@ -180,7 +182,8 @@ yPlus() const const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi]; const scalarField magGradUw(mag(Uw.snGrad())); - const scalar Cmu25 = pow025(Cmu_); + const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); + const scalar yPlusLam = wallCoeffs_.yPlusLam(); auto tyPlus = tmp::New(patch().size(), Zero); auto& yPlus = tyPlus.ref(); @@ -190,7 +193,7 @@ yPlus() const // inertial sublayer yPlus[facei] = Cmu25*y[facei]*sqrt(kwc[facei])/nuw[facei]; - if (yPlusLam_ > yPlus[facei]) + if (yPlusLam > yPlus[facei]) { // viscous sublayer yPlus[facei] = diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C index 827bd54efc..a692abfd79 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C @@ -202,8 +202,8 @@ void Foam::omegaWallFunctionFvPatchScalarField::calculate { const label patchi = patch.index(); - const nutWallFunctionFvPatchScalarField& nutw = - nutWallFunctionFvPatchScalarField::nutw(turbModel, patchi); + const tmp tnutw = turbModel.nut(patchi); + const scalarField& nutw = tnutw(); const scalarField& y = turbModel.y()[patchi]; @@ -217,7 +217,9 @@ void Foam::omegaWallFunctionFvPatchScalarField::calculate const scalarField magGradUw(mag(Uw.snGrad())); - const scalar Cmu25 = pow025(nutw.Cmu()); + const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); + const scalar kappa = wallCoeffs_.kappa(); + const scalar yPlusLam = wallCoeffs_.yPlusLam(); // Set omega and G forAll(nutw, facei) @@ -230,13 +232,13 @@ void Foam::omegaWallFunctionFvPatchScalarField::calculate const scalar omegaVis = 6.0*nuw[facei]/(beta1_*sqr(y[facei])); // Contribution from the inertial sublayer - const scalar omegaLog = sqrt(k[celli])/(Cmu25*nutw.kappa()*y[facei]); + const scalar omegaLog = sqrt(k[celli])/(Cmu25*kappa*y[facei]); switch (blending_) { case blendingType::STEPWISE: { - if (yPlus > nutw.yPlusLam()) + if (yPlus > yPlusLam) { omega0[celli] += w*omegaLog; } @@ -296,14 +298,14 @@ void Foam::omegaWallFunctionFvPatchScalarField::calculate } } - if (!(blending_ == blendingType::STEPWISE) || yPlus > nutw.yPlusLam()) + if (!(blending_ == blendingType::STEPWISE) || yPlus > yPlusLam) { G0[celli] += w *(nutw[facei] + nuw[facei]) *magGradUw[facei] *Cmu25*sqrt(k[celli]) - /(nutw.kappa()*y[facei]); + /(kappa*y[facei]); } } } @@ -321,6 +323,8 @@ void Foam::omegaWallFunctionFvPatchScalarField::writeLocalEntries { os.writeEntry("n", n_); } + + wallCoeffs_.writeEntries(os); } @@ -338,6 +342,7 @@ Foam::omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField initialised_(false), master_(-1), beta1_(0.075), + wallCoeffs_(), G_(), omega_(), cornerWeights_() @@ -358,6 +363,7 @@ Foam::omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField initialised_(false), master_(-1), beta1_(ptf.beta1_), + wallCoeffs_(ptf.wallCoeffs_), G_(), omega_(), cornerWeights_() @@ -393,6 +399,7 @@ Foam::omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField initialised_(false), master_(-1), beta1_(dict.getOrDefault("beta1", 0.075)), + wallCoeffs_(dict), G_(), omega_(), cornerWeights_() @@ -436,6 +443,7 @@ Foam::omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField initialised_(false), master_(-1), beta1_(owfpsf.beta1_), + wallCoeffs_(owfpsf.wallCoeffs_), G_(), omega_(), cornerWeights_() @@ -454,6 +462,7 @@ Foam::omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField initialised_(false), master_(-1), beta1_(owfpsf.beta1_), + wallCoeffs_(owfpsf.wallCoeffs_), G_(), omega_(), cornerWeights_() diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H index 9703832ee2..f448fd2943 100644 --- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H @@ -201,6 +201,7 @@ SourceFiles #define omegaWallFunctionFvPatchScalarField_H #include "fixedValueFvPatchField.H" +#include "wallFunctionCoefficients.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -264,6 +265,9 @@ protected: //- beta1 coefficient scalar beta1_; + //- Wall-function coefficients + wallFunctionCoefficients wallCoeffs_; + //- Local copy of turbulence G field scalarField G_; diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/wallFunction/wallFunctionCoefficients/wallFunctionCoefficients.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/wallFunction/wallFunctionCoefficients/wallFunctionCoefficients.C new file mode 100644 index 0000000000..ecf132e86a --- /dev/null +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/wallFunction/wallFunctionCoefficients/wallFunctionCoefficients.C @@ -0,0 +1,96 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2022 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "wallFunctionCoefficients.H" +#include "dictionary.H" +#include "MinMax.H" + +// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // + +namespace Foam +{ + +//- Estimate the y+ at the intersection of the two sublayers +static scalar calcYPlusLam +( + const scalar kappa, + const scalar E +) +{ + scalar ypl = 11; + + for (int iter = 0; iter < 10; ++iter) + { + ypl = log(max(E*ypl, scalar(1)))/kappa; + } + + return ypl; +} + +} // End namespace Foam + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::wallFunctionCoefficients::wallFunctionCoefficients() +: + Cmu_(0.09), + kappa_(0.41), + E_(9.8), + yPlusLam_(calcYPlusLam(kappa_, E_)) +{} + + +Foam::wallFunctionCoefficients::wallFunctionCoefficients +( + const dictionary& dict +) +: + Cmu_(dict.getOrDefault("Cmu", 0.09)), + kappa_ + ( + dict.getCheckOrDefault("kappa", 0.41, scalarMinMax::ge(SMALL)) + ), + E_(dict.getCheckOrDefault("E", 9.8, scalarMinMax::ge(SMALL))), + yPlusLam_(calcYPlusLam(kappa_, E_)) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::wallFunctionCoefficients::writeEntries +( + Ostream& os +) const +{ + os.writeEntryIfDifferent("Cmu", 0.09, Cmu_); + os.writeEntryIfDifferent("kappa", 0.41, kappa_); + os.writeEntryIfDifferent("E", 9.8, E_); +} + + +// ************************************************************************* // diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/wallFunction/wallFunctionCoefficients/wallFunctionCoefficients.H b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/wallFunction/wallFunctionCoefficients/wallFunctionCoefficients.H new file mode 100644 index 0000000000..d994da37e1 --- /dev/null +++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/wallFunction/wallFunctionCoefficients/wallFunctionCoefficients.H @@ -0,0 +1,129 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2022 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::wallFunctionCoefficients + +Description + Class to describe the wall-function coefficients being + used in the wall function boundary conditions. + +SourceFiles + wallFunctionCoefficients.C + +\*---------------------------------------------------------------------------*/ + +#ifndef Foam_wallFunctionCoefficients_H +#define Foam_wallFunctionCoefficients_H + +#include "scalarFwd.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward Declarations +class dictionary; +class Ostream; + +/*---------------------------------------------------------------------------*\ + Class wallFunctionCoefficients Declaration +\*---------------------------------------------------------------------------*/ + +class wallFunctionCoefficients +{ + // Private Data + + //- Empirical model coefficient + scalar Cmu_; + + //- von Karman constant + scalar kappa_; + + //- Wall roughness parameter + scalar E_; + + //- Estimated y+ value at the intersection + //- of the viscous and inertial sublayers + scalar yPlusLam_; + + +public: + + // Constructors + + //- Construct with default coefficients + wallFunctionCoefficients(); + + //- Construct from dictionary + explicit wallFunctionCoefficients(const dictionary& dict); + + + // Member Functions + + // Access + + //- Return the object: Cmu + scalar Cmu() const noexcept + { + return Cmu_; + } + + //- Return the object: kappa + scalar kappa() const noexcept + { + return kappa_; + } + + //- Return the object: E + scalar E() const noexcept + { + return E_; + } + + //- Return the object: yPlusLam + scalar yPlusLam() const noexcept + { + return yPlusLam_; + } + + + // I-O + + //- Write wall-function coefficients as dictionary entries + void writeEntries(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmEpsilonWallFunction/atmEpsilonWallFunctionFvPatchScalarField.C b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmEpsilonWallFunction/atmEpsilonWallFunctionFvPatchScalarField.C index 76486e01d2..adb64daef1 100644 --- a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmEpsilonWallFunction/atmEpsilonWallFunctionFvPatchScalarField.C +++ b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmEpsilonWallFunction/atmEpsilonWallFunctionFvPatchScalarField.C @@ -45,8 +45,8 @@ void Foam::atmEpsilonWallFunctionFvPatchScalarField::calculate { const label patchi = patch.index(); - const nutWallFunctionFvPatchScalarField& nutw = - nutWallFunctionFvPatchScalarField::nutw(turbModel, patchi); + const tmp tnutw = turbModel.nut(patchi); + const scalarField& nutw = tnutw(); const scalarField& y = turbModel.y()[patchi]; @@ -60,8 +60,10 @@ void Foam::atmEpsilonWallFunctionFvPatchScalarField::calculate const scalarField magGradUw(mag(Uw.snGrad())); - const scalar Cmu25 = pow025(nutw.Cmu()); - const scalar Cmu75 = pow(nutw.Cmu(), 0.75); + const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); + const scalar Cmu75 = pow(wallCoeffs_.Cmu(), 0.75); + const scalar kappa = wallCoeffs_.kappa(); + const scalar yPlusLam = wallCoeffs_.yPlusLam(); const scalar t = db().time().timeOutputValue(); const scalarField z0(z0_->value(t)); @@ -92,16 +94,16 @@ void Foam::atmEpsilonWallFunctionFvPatchScalarField::calculate // (PGVB:Eq. 7, RH:Eq. 8) scalar epsilonc = - w*Cmu75*pow(k[celli], 1.5)/(nutw.kappa()*(y[facei] + z0[facei])); + w*Cmu75*pow(k[celli], 1.5)/(kappa*(y[facei] + z0[facei])); scalar Gc = w *(nutw[facei] + nuw[facei]) *magGradUw[facei] *Cmu25*sqrt(k[celli]) - /(nutw.kappa()*(y[facei] + z0[facei])); + /(kappa*(y[facei] + z0[facei])); - if (lowReCorrection_ && yPlus < nutw.yPlusLam()) + if (lowReCorrection_ && yPlus < yPlusLam) { epsilonc = w*2.0*k[celli]*nuw[facei]/sqr(y[facei] + z0[facei]); Gc = 0; @@ -125,6 +127,8 @@ void Foam::atmEpsilonWallFunctionFvPatchScalarField::writeLocalEntries { z0_->writeData(os); } + + wallCoeffs_.writeEntries(os); } diff --git a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutUWallFunction/atmNutUWallFunctionFvPatchScalarField.C b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutUWallFunction/atmNutUWallFunctionFvPatchScalarField.C index 3095929c17..ffc129b90b 100644 --- a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutUWallFunction/atmNutUWallFunctionFvPatchScalarField.C +++ b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutUWallFunction/atmNutUWallFunctionFvPatchScalarField.C @@ -62,6 +62,8 @@ tmp atmNutUWallFunctionFvPatchScalarField::calcNut() const const tmp tnuw = turbModel.nu(patchi); const scalarField& nuw = tnuw(); + const scalar kappa = wallCoeffs_.kappa(); + const scalar t = db().time().timeOutputValue(); const scalarField z0(z0_->value(t)); @@ -84,7 +86,7 @@ tmp atmNutUWallFunctionFvPatchScalarField::calcNut() const forAll(nutw, facei) { const scalar Edash = (y[facei] + z0[facei])/(z0[facei] + 1e-4); - const scalar uStar = magUpn[facei]*kappa_/log(max(Edash, 1.0 + 1e-4)); + const scalar uStar = magUpn[facei]*kappa/log(max(Edash, 1.0 + 1e-4)); nutw[facei] = sqr(uStar)/max(magUpn[facei], 1e-6)*y[facei] - nuw[facei]; } diff --git a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutWallFunction/atmNutWallFunctionFvPatchScalarField.C b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutWallFunction/atmNutWallFunctionFvPatchScalarField.C index 92a80115ab..0a1d6b3dc7 100644 --- a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutWallFunction/atmNutWallFunctionFvPatchScalarField.C +++ b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutWallFunction/atmNutWallFunctionFvPatchScalarField.C @@ -64,7 +64,8 @@ tmp atmNutWallFunctionFvPatchScalarField::calcNut() const auto tnutw = tmp::New(*this); auto& nutw = tnutw.ref(); - const scalar Cmu25 = pow025(Cmu_); + const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); + const scalar kappa = wallCoeffs_.kappa(); const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi]; const scalarField magUp(mag(Uw.patchInternalField() - Uw)); @@ -95,7 +96,7 @@ tmp atmNutWallFunctionFvPatchScalarField::calcNut() const const scalar Edash = (y[facei] + z0[facei])/(z0[facei] + z0Min_); // (RH:Eq. 6) - const scalar uStarU = magUp[facei]*kappa_/log(max(Edash, 1 + SMALL)); + const scalar uStarU = magUp[facei]*kappa/log(max(Edash, 1 + SMALL)); // (RH:Eq. 7) const scalar uStarK = Cmu25*Foam::sqrt(k[celli]); diff --git a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutkWallFunction/atmNutkWallFunctionFvPatchScalarField.C b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutkWallFunction/atmNutkWallFunctionFvPatchScalarField.C index 3d1cfe6c14..79e7ef7e85 100644 --- a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutkWallFunction/atmNutkWallFunctionFvPatchScalarField.C +++ b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutkWallFunction/atmNutkWallFunctionFvPatchScalarField.C @@ -65,7 +65,8 @@ tmp atmNutkWallFunctionFvPatchScalarField::calcNut() const auto tnutw = tmp::New(*this); auto& nutw = tnutw.ref(); - const scalar Cmu25 = pow025(Cmu_); + const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); + const scalar kappa = wallCoeffs_.kappa(); const scalar t = db().time().timeOutputValue(); const scalarField z0(z0_->value(t)); @@ -94,7 +95,7 @@ tmp atmNutkWallFunctionFvPatchScalarField::calcNut() const const scalar yPlus = uStar*y[facei]/nuw[facei]; const scalar Edash = (y[facei] + z0[facei])/z0[facei]; - nutw[facei] = nuw[facei]*(yPlus*kappa_/log(max(Edash, 1 + 1e-4)) - 1); + nutw[facei] = nuw[facei]*(yPlus*kappa/log(max(Edash, 1 + 1e-4)) - 1); } if (boundNut_) diff --git a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmOmegaWallFunction/atmOmegaWallFunctionFvPatchScalarField.C b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmOmegaWallFunction/atmOmegaWallFunctionFvPatchScalarField.C index 9c20bfa3a0..3c8b247742 100644 --- a/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmOmegaWallFunction/atmOmegaWallFunctionFvPatchScalarField.C +++ b/src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmOmegaWallFunction/atmOmegaWallFunctionFvPatchScalarField.C @@ -45,8 +45,8 @@ void Foam::atmOmegaWallFunctionFvPatchScalarField::calculate { const label patchi = patch.index(); - const nutWallFunctionFvPatchScalarField& nutw = - nutWallFunctionFvPatchScalarField::nutw(turbModel, patchi); + const tmp tnutw = turbModel.nut(patchi); + const scalarField& nutw = tnutw(); const scalarField& y = turbModel.y()[patchi]; @@ -60,7 +60,8 @@ void Foam::atmOmegaWallFunctionFvPatchScalarField::calculate const scalarField magGradUw(mag(Uw.snGrad())); - const scalar Cmu25 = pow025(nutw.Cmu()); + const scalar Cmu25 = pow025(wallCoeffs_.Cmu()); + const scalar kappa = wallCoeffs_.kappa(); const scalar t = db().time().timeOutputValue(); const scalarField z0(z0_->value(t)); @@ -88,14 +89,14 @@ void Foam::atmOmegaWallFunctionFvPatchScalarField::calculate const scalar w = cornerWeights[facei]; omega0[celli] += - w*sqrt(k[celli])/(Cmu25*nutw.kappa()*(y[facei] + z0[facei])); + w*sqrt(k[celli])/(Cmu25*kappa*(y[facei] + z0[facei])); G0[celli] += w *(nutw[facei] + nuw[facei]) *magGradUw[facei] *Cmu25*sqrt(k[celli]) - /(nutw.kappa()*(y[facei] + z0[facei])); + /(kappa*(y[facei] + z0[facei])); } } @@ -109,6 +110,8 @@ void Foam::atmOmegaWallFunctionFvPatchScalarField::writeLocalEntries { z0_->writeData(os); } + + wallCoeffs_.writeEntries(os); } diff --git a/src/regionModels/surfaceFilmModels/derivedFvPatchFields/wallFunctions/nutkFilmWallFunction/nutkFilmWallFunctionFvPatchScalarField.C b/src/regionModels/surfaceFilmModels/derivedFvPatchFields/wallFunctions/nutkFilmWallFunction/nutkFilmWallFunctionFvPatchScalarField.C index d39e106e09..72718f1841 100644 --- a/src/regionModels/surfaceFilmModels/derivedFvPatchFields/wallFunctions/nutkFilmWallFunction/nutkFilmWallFunctionFvPatchScalarField.C +++ b/src/regionModels/surfaceFilmModels/derivedFvPatchFields/wallFunctions/nutkFilmWallFunction/nutkFilmWallFunctionFvPatchScalarField.C @@ -92,7 +92,8 @@ tmp nutkFilmWallFunctionFvPatchScalarField::calcUTau const tmp tnuw = turbModel.nu(patchi); const scalarField& nuw = tnuw(); - const scalar Cmu25 = pow(Cmu_, 0.25); + const scalar Cmu25 = pow(wallCoeffs_.Cmu(), 0.25); + const scalar kappa = wallCoeffs_.kappa(); forAll(uTau, facei) { @@ -108,7 +109,7 @@ tmp nutkFilmWallFunctionFvPatchScalarField::calcUTau if (yPlus > yPlusCrit_) { scalar expTerm = exp(min(50.0, B_*mStar)); - scalar powTerm = pow(yPlus, mStar/kappa_); + scalar powTerm = pow(yPlus, mStar/kappa); factor = mStar/(expTerm*powTerm - 1.0 + ROOTVSMALL); } else @@ -163,8 +164,6 @@ void nutkFilmWallFunctionFvPatchScalarField::writeLocalEntries ); os.writeEntryIfDifferent("B", 5.5, B_); os.writeEntryIfDifferent("yPlusCrit", 11.05, yPlusCrit_); - os.writeEntryIfDifferent("Cmu", 0.09, Cmu_); - os.writeEntryIfDifferent("kappa", 0.41, kappa_); } @@ -266,7 +265,7 @@ tmp nutkFilmWallFunctionFvPatchScalarField::yPlus() const void nutkFilmWallFunctionFvPatchScalarField::write(Ostream& os) const { - fvPatchField::write(os); + nutWallFunctionFvPatchScalarField::write(os); writeLocalEntries(os); writeEntry("value", os); }