From f57bfe7c0529372a07d0062fd8f0fb35e503f55f Mon Sep 17 00:00:00 2001 From: andy Date: Thu, 23 Jul 2009 14:22:05 +0100 Subject: [PATCH] adding yPlus calcs for incompressible S-A models --- ...ndardRoughWallFunctionFvPatchScalarField.C | 263 ++++++++++-------- ...ndardRoughWallFunctionFvPatchScalarField.H | 3 + ...asStandardWallFunctionFvPatchScalarField.C | 82 ++++-- ...asStandardWallFunctionFvPatchScalarField.H | 3 + ...rtAllmarasWallFunctionFvPatchScalarField.C | 4 - 5 files changed, 202 insertions(+), 153 deletions(-) diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C index 7359e6f081..d6b96ed5b4 100644 --- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C +++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C @@ -39,6 +39,134 @@ namespace incompressible namespace RASModels { +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +tmp +nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus +( + const scalarField& magUp +) const +{ + const label patchI = patch().index(); + + const RASModel& rasModel = db().lookupObject("RASProperties"); + const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_); + const scalarField& y = rasModel.y()[patchI]; + const scalarField& nuw = rasModel.nu().boundaryField()[patchI]; + + tmp tyPlus(new scalarField(patch().size(), 0.0)); + scalarField& yPlus = tyPlus(); + + if (roughnessHeight_ > 0.0) + { + // Rough Walls + const scalar c_1 = 1/(90 - 2.25) + roughnessConstant_; + static const scalar c_2 = 2.25/(90 - 2.25); + static const scalar c_3 = 2.0*atan(1.0)/log(90/2.25); + static const scalar c_4 = c_3*log(2.25); + + //if (KsPlusBasedOnYPlus_) + { + // If KsPlus is based on YPlus the extra term added to the law + // of the wall will depend on yPlus + forAll(yPlus, facei) + { + const scalar magUpara = magUp[facei]; + const scalar Re = magUpara*y[facei]/nuw[facei]; + const scalar kappaRe = kappa_*Re; + + scalar yp = yPlusLam; + const scalar ryPlusLam = 1.0/yp; + + int iter = 0; + scalar yPlusLast = 0.0; + scalar dKsPlusdYPlus = roughnessHeight_/y[facei]; + + // Enforce the roughnessHeight to be less than the distance to + // the first cell centre + if (dKsPlusdYPlus > 1) + { + dKsPlusdYPlus = 1; + } + + // Additional tuning parameter (fudge factor) - nominally = 1 + dKsPlusdYPlus *= roughnessFudgeFactor_; + + do + { + yPlusLast = yp; + + // The non-dimensional roughness height + scalar KsPlus = yp*dKsPlusdYPlus; + + // The extra term in the law-of-the-wall + scalar G = 0.0; + + scalar yPlusGPrime = 0.0; + + if (KsPlus >= 90) + { + const scalar t_1 = 1 + roughnessConstant_*KsPlus; + G = log(t_1); + yPlusGPrime = roughnessConstant_*KsPlus/t_1; + } + else if (KsPlus > 2.25) + { + const scalar t_1 = c_1*KsPlus - c_2; + const scalar t_2 = c_3*log(KsPlus) - c_4; + const scalar sint_2 = sin(t_2); + const scalar logt_1 = log(t_1); + G = logt_1*sint_2; + yPlusGPrime = + (c_1*sint_2*KsPlus/t_1) + (c_3*logt_1*cos(t_2)); + } + + scalar denom = 1.0 + log(E_*yp) - G - yPlusGPrime; + if (mag(denom) > VSMALL) + { + yp = (kappaRe + yp*(1 - yPlusGPrime))/denom; + } + } while + ( + mag(ryPlusLam*(yp - yPlusLast)) > 0.0001 + && ++iter < 10 + && yp > VSMALL + ); + + yPlus[facei] = yp; + } + } + } + else + { + // Smooth Walls + forAll(yPlus, facei) + { + const scalar magUpara = magUp[facei]; + const scalar Re = magUpara*y[facei]/nuw[facei]; + const scalar kappaRe = kappa_*Re; + + scalar yp = yPlusLam; + const scalar ryPlusLam = 1.0/yp; + + int iter = 0; + scalar yPlusLast = 0.0; + + do + { + yPlusLast = yp; + yPlus = (kappaRe + yp)/(1.0 + log(E_*yp)); + + } while(mag(ryPlusLam*(yp - yPlusLast)) > 0.0001 && ++iter < 10); + + yPlus[facei] = yp; + } + } + + return tyPlus; +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField:: @@ -123,131 +251,24 @@ nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcNut() const const RASModel& rasModel = db().lookupObject("RASProperties"); const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_); const scalarField& y = rasModel.y()[patchI]; - const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI]; const scalarField& nuw = rasModel.nu().boundaryField()[patchI]; // The flow velocity at the adjacent cell centre - scalarField magUp = mag(Uw.patchInternalField() - Uw); + const scalarField magUp = mag(Uw.patchInternalField() - Uw); + + tmp tyPlus = calcYPlus(magUp); + scalarField& yPlus = tyPlus(); tmp tnutw(new scalarField(patch().size(), 0.0)); scalarField& nutw = tnutw(); - if (roughnessHeight_ > 0.0) + forAll(yPlus, facei) { - // Rough Walls - const scalar c_1 = 1/(90 - 2.25) + roughnessConstant_; - static const scalar c_2 = 2.25/(90 - 2.25); - static const scalar c_3 = 2.0*atan(1.0)/log(90/2.25); - static const scalar c_4 = c_3*log(2.25); - - //if (KsPlusBasedOnYPlus_) + if (yPlus[facei] > yPlusLam) { - // If KsPlus is based on YPlus the extra term added to the law - // of the wall will depend on yPlus. - forAll(nutw, facei) - { - const scalar magUpara = magUp[facei]; - const scalar Re = magUpara*y[facei]/nuw[facei]; - const scalar kappaRe = kappa_*Re; - - scalar yPlus = yPlusLam; - const scalar ryPlusLam = 1.0/yPlus; - - int iter = 0; - scalar yPlusLast = 0.0; - scalar dKsPlusdYPlus = roughnessHeight_/y[facei]; - - // Enforce the roughnessHeight to be less than the distance to - // the first cell centre - if (dKsPlusdYPlus > 1) - { - dKsPlusdYPlus = 1; - } - - // Additional tuning parameter (fudge factor) - nominally = 1 - dKsPlusdYPlus *= roughnessFudgeFactor_; - - do - { - yPlusLast = yPlus; - - // The non-dimensional roughness height. - scalar KsPlus = yPlus*dKsPlusdYPlus; - - // The extra term in the law-of-the-wall. - scalar G = 0.0; - - scalar yPlusGPrime = 0.0; - - if (KsPlus >= 90) - { - const scalar t_1 = 1 + roughnessConstant_*KsPlus; - G = log(t_1); - yPlusGPrime = roughnessConstant_*KsPlus/t_1; - } - else if (KsPlus > 2.25) - { - const scalar t_1 = c_1*KsPlus - c_2; - const scalar t_2 = c_3*log(KsPlus) - c_4; - const scalar sint_2 = sin(t_2); - const scalar logt_1 = log(t_1); - G = logt_1*sint_2; - yPlusGPrime = - (c_1*sint_2*KsPlus/t_1) + (c_3*logt_1*cos(t_2)); - } - - scalar denom = 1.0 + log(E_*yPlus) - G - yPlusGPrime; - if (mag(denom) > VSMALL) - { - yPlus = (kappaRe + yPlus*(1 - yPlusGPrime))/denom; - } - else - { - // Ensure immediate end and nutw = 0. - yPlus = 0; - } - - } while - ( - mag(ryPlusLam*(yPlus - yPlusLast)) > 0.0001 - && ++iter < 10 - && yPlus > VSMALL - ); - - if (yPlus > yPlusLam) - { - nutw[facei] = nuw[facei]*(yPlus*yPlus/Re - 1); - } - } - } - } - else - { - // Smooth Walls. - forAll(nutw, facei) - { - const scalar magUpara = magUp[facei]; - const scalar Re = magUpara*y[facei]/nuw[facei]; - const scalar kappaRe = kappa_*Re; - - scalar yPlus = yPlusLam; - const scalar ryPlusLam = 1.0/yPlus; - - int iter = 0; - scalar yPlusLast = 0.0; - - do - { - yPlusLast = yPlus; - yPlus = (kappaRe + yPlus)/(1.0 + log(E_*yPlus)); - - } while(mag(ryPlusLam*(yPlus - yPlusLast)) > 0.0001 && ++iter < 10); - - if (yPlus > yPlusLam) - { - nutw[facei] = nuw[facei]*(yPlus*yPlus/Re - 1); - } + const scalar Re = magUp[facei]*y[facei]/nuw[facei]; + nutw[facei] = nuw[facei]*(sqr(yPlus[facei])/Re - 1); } } @@ -258,13 +279,13 @@ nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcNut() const tmp nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::yPlus() const { - notImplemented - ( - "nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::yPlus()" - "const" - ); + const label patchI = patch().index(); - return tmp(NULL); + const RASModel& rasModel = db().lookupObject("RASProperties"); + const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI]; + const scalarField magUp = mag(Uw.patchInternalField() - Uw); + + return calcYPlus(magUp); } diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H index 5e819756db..8d72e5ad13 100644 --- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H +++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H @@ -72,6 +72,9 @@ class nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField // Protected member functions + //- Calculate yPLus + virtual tmp calcYPlus(const scalarField& magUp) const; + //- Calculate the turbulence viscosity virtual tmp calcNut() const; diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C index aef85e9b38..25c9b87d98 100644 --- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C +++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C @@ -39,6 +39,48 @@ namespace incompressible namespace RASModels { +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +tmp +nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcYPlus +( + const scalarField& magUp +) const +{ + const label patchI = patch().index(); + + const RASModel& rasModel = db().lookupObject("RASProperties"); + const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_); + const scalarField& y = rasModel.y()[patchI]; + const scalarField& nuw = rasModel.nu().boundaryField()[patchI]; + + tmp tyPlus(new scalarField(patch().size(), 0.0)); + scalarField& yPlus = tyPlus(); + + forAll(yPlus, facei) + { + scalar kappaRe = kappa_*magUp[facei]*y[facei]/nuw[facei]; + + scalar yp = yPlusLam; + scalar ryPlusLam = 1.0/yp; + + int iter = 0; + scalar yPlusLast = 0.0; + + do + { + yPlusLast = yp; + yPlus = (kappaRe + yp)/(1.0 + log(E_*yp)); + + } while(mag(ryPlusLam*(yp - yPlusLast)) > 0.01 && ++iter < 10 ); + + yPlus[facei] = yp; + } + + return tyPlus; +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // nutSpalartAllmarasStandardWallFunctionFvPatchScalarField:: @@ -107,37 +149,22 @@ nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcNut() const const RASModel& rasModel = db().lookupObject("RASProperties"); const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_); - const scalarField& y = rasModel.y()[patchI]; - const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI]; - const scalarField magUp = mag(Uw.patchInternalField() - Uw); - const scalarField& nuw = rasModel.nu().boundaryField()[patchI]; + tmp tyPlus = calcYPlus(magUp); + scalarField& yPlus = tyPlus(); + tmp tnutw(new scalarField(patch().size(), 0.0)); scalarField& nutw = tnutw(); - forAll(nutw, facei) + forAll(yPlus, facei) { - scalar kappaRe = kappa_*magUp[facei]*y[facei]/nuw[facei]; - - scalar yPlus = yPlusLam; - scalar ryPlusLam = 1.0/yPlus; - - int iter = 0; - scalar yPlusLast = 0.0; - - do + if (yPlus[facei] > yPlusLam) { - yPlusLast = yPlus; - yPlus = (kappaRe + yPlus)/(1.0 + log(E_*yPlus)); - - } while(mag(ryPlusLam*(yPlus - yPlusLast)) > 0.01 && ++iter < 10 ); - - if (yPlus > yPlusLam) - { - nutw[facei] = nuw[facei]*(yPlus*kappa_/log(E_*yPlus) - 1.0); + nutw[facei] = + nuw[facei]*(yPlus[facei]*kappa_/log(E_*yPlus[facei]) - 1.0); } } @@ -148,13 +175,12 @@ nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcNut() const tmp nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::yPlus() const { - notImplemented - ( - "nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::yPlus() " - "const" - ); + const label patchI = patch().index(); + const RASModel& rasModel = db().lookupObject("RASProperties"); + const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI]; + const scalarField magUp = mag(Uw.patchInternalField() - Uw); - return tmp(NULL); + return calcYPlus(magUp); } diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.H index 80e525569a..f54bb57c8a 100644 --- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.H +++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.H @@ -60,6 +60,9 @@ protected: // Protected member functions + //- Calculate yPLus + virtual tmp calcYPlus(const scalarField& magUp) const; + //- Calculate the turbulence viscosity virtual tmp calcNut() const; diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasWallFunction/nutSpalartAllmarasWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasWallFunction/nutSpalartAllmarasWallFunctionFvPatchScalarField.C index 80aa3e3199..ff02c3043b 100644 --- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasWallFunction/nutSpalartAllmarasWallFunctionFvPatchScalarField.C +++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasWallFunction/nutSpalartAllmarasWallFunctionFvPatchScalarField.C @@ -167,9 +167,7 @@ nutSpalartAllmarasWallFunctionFvPatchScalarField::calcNut() const const RASModel& rasModel = db().lookupObject("RASProperties"); const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI]; - const scalarField magGradU = mag(Uw.snGrad()); - const scalarField& nuw = rasModel.nu().boundaryField()[patchI]; return max(0.0, sqr(calcUTau(magGradU))/magGradU - nuw); @@ -183,9 +181,7 @@ nutSpalartAllmarasWallFunctionFvPatchScalarField::yPlus() const const RASModel& rasModel = db().lookupObject("RASProperties"); const scalarField& y = rasModel.y()[patchI]; - const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI]; - const scalarField& nuw = rasModel.nu().boundaryField()[patchI]; return y*calcUTau(mag(Uw.snGrad()))/nuw;