diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C index 172cfb0db5..351f4ca80d 100644 --- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C +++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C @@ -39,6 +39,137 @@ namespace compressible namespace RASModels { +tmp +mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus +( + const scalarField& magUp +) const +{ + const label patchI = patch().index(); + + const RASModel& rasModel = db().lookupObject("RASProperties"); + const scalarField& y = rasModel.y()[patchI]; + const scalarField& muw = rasModel.mu().boundaryField()[patchI]; + const fvPatchScalarField& rho = rasModel.rho().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 = rho[facei]*magUpara*y[facei]/muw[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] = max(0.0, yp); + } + } + } + else + { + // Smooth Walls + forAll(yPlus, facei) + { + const scalar magUpara = magUp[facei]; + const scalar Re = rho[facei]*magUpara*y[facei]/muw[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; + yp = (kappaRe + yp)/(1.0 + log(E_*yp)); + + } while + ( + mag(ryPlusLam*(yp - yPlusLast)) > 0.0001 + && ++iter < 10 + ); + + yPlus[facei] = max(0.0, yp); + } + } + + return tyPlus; +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField:: @@ -123,140 +254,24 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcMut() 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& muw = rasModel.mu().boundaryField()[patchI]; + const fvPatchScalarField& rho = rasModel.rho().boundaryField()[patchI]; scalarField magUp = mag(Uw.patchInternalField() - Uw); - const fvPatchScalarField& rho = rasModel.rho().boundaryField()[patchI]; + tmp tyPlus = calcYPlus(magUp); + scalarField& yPlus = tyPlus(); tmp tmutw(new scalarField(patch().size(), 0.0)); scalarField& mutw = tmutw(); - 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(mutw, facei) - { - const scalar magUpara = magUp[facei]; - const scalar Re = rho[facei]*magUpara*y[facei]/muw[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; - if( yPlus < 0 ) - { - yPlus = 0; - } - } - else - { - // Ensure immediate end and mutw = 0 - yPlus = 0; - } - - } while - ( - mag(ryPlusLam*(yPlus - yPlusLast)) > 0.0001 - && ++iter < 10 - && yPlus > VSMALL - ); - - if (yPlus > yPlusLam) - { - mutw[facei] = muw[facei]*(yPlus*yPlus/Re - 1); - } - } - } - } - else - { - // Smooth Walls - forAll(mutw, facei) - { - const scalar magUpara = magUp[facei]; - const scalar Re = rho[facei]*magUpara*y[facei]/muw[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) - { - mutw[facei] = muw[facei]*(yPlus*yPlus/Re - 1); - } + const scalar Re = rho[facei]*magUp[facei]*y[facei]/muw[facei]; + mutw[facei] = muw[facei]*(sqr(yPlus[facei])/Re - 1); } } @@ -267,13 +282,13 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcMut() const tmp mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::yPlus() const { - notImplemented - ( - "mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::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); } @@ -283,6 +298,8 @@ void mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::write ) const { fixedValueFvPatchScalarField::write(os); + os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl; + os.writeKeyword("E") << E_ << token::END_STATEMENT << nl; os.writeKeyword("roughnessHeight") << roughnessHeight_ << token::END_STATEMENT << nl; os.writeKeyword("roughnessConstant") diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H index e2434c8183..2414825d21 100644 --- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H +++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H @@ -74,6 +74,9 @@ protected: // Protected member functions + //- Calculate yPLus + virtual tmp calcYPlus(const scalarField& magUp) const; + //- Calculate the turbulence viscosity virtual tmp calcMut() const; diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C index f2337dab47..493cf1598d 100644 --- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C +++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C @@ -39,6 +39,50 @@ namespace compressible namespace RASModels { +// * * * * * * * * * * * * Protected 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 fvPatchScalarField& rhow = rasModel.rho().boundaryField()[patchI]; + const fvPatchScalarField& muw = rasModel.mu().boundaryField()[patchI]; + + tmp tyPlus(new scalarField(patch().size(), 0.0)); + scalarField& yPlus = tyPlus(); + + forAll(yPlus, faceI) + { + scalar kappaRe = kappa_*magUp[facei]*y[faceI]/(muw[faceI]/rhow[faceI]); + + scalar yp = yPlusLam; + scalar ryPlusLam = 1.0/yp; + + int iter = 0; + scalar yPlusLast = 0.0; + + do + { + yPlusLast = yp; + yp = (kappaRe + yp)/(1.0 + log(E_*yp)); + + } while (mag(ryPlusLam*(yp - yPlusLast)) > 0.01 && ++iter < 10); + + yPlus[facei] = max(0.0, yp); + } + + return tyPlus; + +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // mutSpalartAllmarasStandardWallFunctionFvPatchScalarField:: @@ -108,40 +152,23 @@ mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcMut() 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 fvPatchScalarField& rhow = rasModel.rho().boundaryField()[patchI]; - const fvPatchScalarField& muw = rasModel.mu().boundaryField()[patchI]; + tmp tyPlus = calcYPlus(magUp); + scalarField& yPlus = tyPlus(); + tmp tmutw(new scalarField(patch().size(), 0.0)); scalarField& mutw = tmutw(); - forAll(mutw, faceI) + forAll(yPlus, faceI) { - scalar magUpara = magUp[faceI]; - - scalar kappaRe = kappa_*magUpara*y[faceI]/(muw[faceI]/rhow[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) - { - mutw[faceI] = muw[faceI]*(yPlus*kappa_/log(E_*yPlus) - 1.0); + mutw[faceI] = + muw[faceI]*(yPlus[facei]*kappa_/log(E_*yPlus[facei]) - 1.0); } } @@ -152,13 +179,12 @@ mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcMut() const tmp mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::yPlus() const { - notImplemented - ( - "mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::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/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.H b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.H index 6813e59667..09ba112ce4 100644 --- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.H +++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.H @@ -60,6 +60,9 @@ protected: // Protected member functions + //- Calculate yPLus + virtual tmp calcYPlus(const scalarField& magUp) const; + //- Calculate the turbulence viscosity virtual tmp calcMut() const; 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 d6b96ed5b4..f44da462b9 100644 --- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C +++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C @@ -155,7 +155,7 @@ nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus do { yPlusLast = yp; - yPlus = (kappaRe + yp)/(1.0 + log(E_*yp)); + yp = (kappaRe + yp)/(1.0 + log(E_*yp)); } while(mag(ryPlusLam*(yp - yPlusLast)) > 0.0001 && ++iter < 10); 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 25c9b87d98..450c857a75 100644 --- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C +++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C @@ -39,7 +39,7 @@ namespace incompressible namespace RASModels { -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // tmp nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcYPlus @@ -70,11 +70,11 @@ nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcYPlus do { yPlusLast = yp; - yPlus = (kappaRe + yp)/(1.0 + log(E_*yp)); + yp = (kappaRe + yp)/(1.0 + log(E_*yp)); } while(mag(ryPlusLam*(yp - yPlusLast)) > 0.01 && ++iter < 10 ); - yPlus[facei] = yp; + yPlus[facei] = max(0.0, yp); } return tyPlus;