zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
mattijs
2009-07-23 17:17:45 +01:00
parent b01fc1d7e3 be65512449
commit 76ca0da080
29 changed files with 623 additions and 603 deletions
Download Patch File Download Diff File Expand all files Collapse all files

46
src/thermophysicalModels/pdfs/RosinRammler/RosinRammler.C
View File

@ -24,33 +24,21 @@ License
\*---------------------------------------------------------------------------*/
#include "RosinRammler.H"
#include "addToRunTimeSelectionTable.H"
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(RosinRammler, 0);
namespace Foam
{
defineTypeNameAndDebug(RosinRammler, 0);
addToRunTimeSelectionTable
(
pdf,
RosinRammler,
dictionary
);
addToRunTimeSelectionTable(pdf, RosinRammler, dictionary);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
RosinRammler::RosinRammler
(
const dictionary& dict,
Random& rndGen
)
Foam::RosinRammler::RosinRammler(const dictionary& dict, Random& rndGen)
:
pdf(dict, rndGen),
pdfDict_(dict.subDict(typeName + "PDF")),
@ -80,10 +68,10 @@ RosinRammler::RosinRammler
// find max value so that it can be normalized to 1.0
scalar sMax = 0;
label n = d_.size();
for(label i=0; i<n; i++)
for (label i=0; i<n; i++)
{
scalar s = exp(-1.0);
for(label j=0; j<n; j++)
for (label j=0; j<n; j++)
{
if (i!=j)
{
@ -96,11 +84,10 @@ RosinRammler::RosinRammler
sMax = max(sMax, s);
}
for(label i=0; i<n; i++)
for (label i=0; i<n; i++)
{
ls_[i] /= sMax;
}
}
@ -112,20 +99,20 @@ Foam::RosinRammler::~RosinRammler()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
scalar RosinRammler::sample() const
Foam::scalar Foam::RosinRammler::sample() const
{
scalar y = 0;
scalar x = 0;
label n = d_.size();
bool success = false;
while(!success)
while (!success)
{
x = minValue_ + range_*rndGen_.scalar01();
y = rndGen_.scalar01();
scalar p = 0.0;
for(label i=0; i<n; i++)
for (label i=0; i<n; i++)
{
scalar xx = pow(x/d_[i], n_[i]);
p += ls_[i]*xx*exp(-xx);
@ -140,16 +127,17 @@ scalar RosinRammler::sample() const
return x;
}
scalar RosinRammler::minValue() const
Foam::scalar Foam::RosinRammler::minValue() const
{
return minValue_;
}
scalar RosinRammler::maxValue() const
Foam::scalar Foam::RosinRammler::maxValue() const
{
return maxValue_;
}
// ************************************************************************* //
}
// ************************************************************************* //

17
src/thermophysicalModels/pdfs/RosinRammler/RosinRammler.H
View File

@ -72,10 +72,11 @@ class RosinRammler
scalar range_;
public:
//- Runtime type information
TypeName("RosinRammler");
TypeName("RosinRammler");
// Constructors
@ -88,18 +89,20 @@ public:
);
// Destructor
~RosinRammler();
//- Destructor
virtual ~RosinRammler();
// Member Functions
scalar sample() const;
//- Sample the pdf
virtual scalar sample() const;
scalar minValue() const;
scalar maxValue() const;
//- Return the minimum value
virtual scalar minValue() const;
//- Return the maximum value
virtual scalar maxValue() const;
};

46
src/thermophysicalModels/pdfs/exponential/exponential.C
View File

@ -28,29 +28,17 @@ License
#include "exponential.H"
#include "addToRunTimeSelectionTable.H"
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(exponential, 0);
addToRunTimeSelectionTable
(
pdf,
exponential,
dictionary
);
namespace Foam
{
defineTypeNameAndDebug(exponential, 0);
addToRunTimeSelectionTable(pdf, exponential, dictionary);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
exponential::exponential
(
const dictionary& dict,
Random& rndGen
)
Foam::exponential::exponential(const dictionary& dict, Random& rndGen)
:
pdf(dict, rndGen),
pdfDict_(dict.subDict(typeName + "PDF")),
@ -70,10 +58,10 @@ exponential::exponential
scalar sMax = 0;
label n = lambda_.size();
for(label i=0; i<n; i++)
for (label i=0; i<n; i++)
{
scalar s = lambda_[i]*exp(-lambda_[i]*minValue_);
for(label j=0; j<n; j++)
for (label j=0; j<n; j++)
{
if (i!=j)
{
@ -89,7 +77,6 @@ exponential::exponential
{
ls_[i] /= sMax;
}
}
@ -101,14 +88,14 @@ Foam::exponential::~exponential()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
scalar exponential::sample() const
Foam::scalar Foam::exponential::sample() const
{
scalar y = 0;
scalar x = 0;
label n = lambda_.size();
bool success = false;
while(!success)
while (!success)
{
x = minValue_ + range_*rndGen_.scalar01();
y = rndGen_.scalar01();
@ -119,25 +106,26 @@ scalar exponential::sample() const
p += ls_[i]*exp(-lambda_[i]*x);
}
if (y<p)
if (y<p)
{
success = true;
}
}
return x;
}
scalar exponential::minValue() const
Foam::scalar Foam::exponential::minValue() const
{
return minValue_;
}
scalar exponential::maxValue() const
Foam::scalar Foam::exponential::maxValue() const
{
return maxValue_;
}
// ************************************************************************* //
}
// ************************************************************************* //

18
src/thermophysicalModels/pdfs/exponential/exponential.H
View File

@ -46,7 +46,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class exponential Declaration
Class exponential Declaration
\*---------------------------------------------------------------------------*/
class exponential
@ -69,7 +69,7 @@ class exponential
public:
//- Runtime type information
TypeName("exponential");
TypeName("exponential");
// Constructors
@ -82,18 +82,20 @@ public:
);
// Destructor
~exponential();
//- Destructor
virtual ~exponential();
// Member Functions
scalar sample() const;
//- Sample the pdf
virtual scalar sample() const;
scalar minValue() const;
scalar maxValue() const;
//- Return the minimum value
virtual scalar minValue() const;
//- Return the maximum value
virtual scalar maxValue() const;
};

60
src/thermophysicalModels/pdfs/general/general.C
View File

@ -24,33 +24,20 @@ License
\*---------------------------------------------------------------------------*/
#include "general.H"
#include "addToRunTimeSelectionTable.H"
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(general, 0);
addToRunTimeSelectionTable
(
pdf,
general,
dictionary
);
namespace Foam
{
defineTypeNameAndDebug(general, 0);
addToRunTimeSelectionTable(pdf, general, dictionary);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
general::general
(
const dictionary& dict,
Random& rndGen
)
Foam::general::general(const dictionary& dict, Random& rndGen)
:
pdf(dict, rndGen),
pdfDict_(dict.subDict(typeName + "PDF")),
@ -63,16 +50,15 @@ general::general
// normalize the pdf
scalar yMax = 0;
for(label i=0; i<nEntries_; i++)
for (label i=0; i<nEntries_; i++)
{
yMax = max(yMax, xy_[i][1]);
}
for(label i=0; i<nEntries_; i++)
for (label i=0; i<nEntries_; i++)
{
xy_[i][1] /= yMax;
}
}
@ -84,21 +70,21 @@ Foam::general::~general()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
scalar general::sample() const
Foam::scalar Foam::general::sample() const
{
scalar y = 0;
scalar x = 0;
bool success = false;
while(!success)
while (!success)
{
x = minValue_ + range_*rndGen_.scalar01();
y = rndGen_.scalar01();
bool intervalFound = false;
label i = -1;
while(!intervalFound)
while (!intervalFound)
{
i++;
if ( (x>xy_[i][0]) && (x<xy_[i+1][0]) )
@ -106,27 +92,33 @@ scalar general::sample() const
intervalFound = true;
}
}
scalar p = xy_[i][1] + (x-xy_[i][0])*(xy_[i+1][1]-xy_[i][1])/(xy_[i+1][0]-xy_[i][0]);
if (y<p)
scalar p =
xy_[i][1]
+ (x-xy_[i][0])
*(xy_[i+1][1]-xy_[i][1])
/(xy_[i+1][0]-xy_[i][0]);
if (y<p)
{
success = true;
}
}
return x;
}
scalar general::minValue() const
Foam::scalar Foam::general::minValue() const
{
return minValue_;
}
scalar general::maxValue() const
Foam::scalar Foam::general::maxValue() const
{
return maxValue_;
}
// ************************************************************************* //
}
// ************************************************************************* //

20
src/thermophysicalModels/pdfs/general/general.H
View File

@ -62,16 +62,18 @@ class general
List<pair> xy_;
label nEntries_;
//- min and max values of the distribution
scalar minValue_;
scalar maxValue_;
scalar range_;
public:
//- Runtime type information
TypeName("general");
TypeName("general");
// Constructors
@ -85,17 +87,19 @@ public:
// Destructor
~general();
virtual ~general();
// Member Functions
scalar sample() const;
//- Sample the pdf
virtual scalar sample() const;
scalar minValue() const;
scalar maxValue() const;
//- Return the minimum value
virtual scalar minValue() const;
//- Return the maximum value
virtual scalar maxValue() const;
};
@ -105,10 +109,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//#include "generalI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

51
src/thermophysicalModels/pdfs/normal/normal.C
View File

@ -24,33 +24,20 @@ License
\*---------------------------------------------------------------------------*/
#include "normal.H"
#include "addToRunTimeSelectionTable.H"
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(normal, 0);
addToRunTimeSelectionTable
(
pdf,
normal,
dictionary
);
namespace Foam
{
defineTypeNameAndDebug(normal, 0);
addToRunTimeSelectionTable(pdf, normal, dictionary);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
normal::normal
(
const dictionary& dict,
Random& rndGen
)
Foam::normal::normal(const dictionary& dict, Random& rndGen)
:
pdf(dict, rndGen),
pdfDict_(dict.subDict(typeName + "PDF")),
@ -63,11 +50,11 @@ normal::normal
{
scalar sMax = 0;
label n = strength_.size();
for(label i=0; i<n; i++)
for (label i=0; i<n; i++)
{
scalar x = expectation_[i];
scalar s = strength_[i];
for(label j=0; j<n; j++)
for (label j=0; j<n; j++)
{
if (i!=j)
{
@ -80,11 +67,10 @@ normal::normal
sMax = max(sMax, s);
}
for(label i=0; i<n; i++)
for (label i=0; i<n; i++)
{
strength_[i] /= sMax;
}
}
@ -96,20 +82,20 @@ Foam::normal::~normal()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
scalar normal::sample() const
Foam::scalar Foam::normal::sample() const
{
scalar y = 0;
scalar x = 0;
label n = expectation_.size();
bool success = false;
while(!success)
while (!success)
{
x = minValue_ + range_*rndGen_.scalar01();
y = rndGen_.scalar01();
scalar p = 0.0;
for(label i=0; i<n; i++)
for (label i=0; i<n; i++)
{
scalar nu = expectation_[i];
scalar sigma = variance_[i];
@ -118,25 +104,26 @@ scalar normal::sample() const
p += s*exp(-0.5*v*v);
}
if (y<p)
if (y<p)
{
success = true;
}
}
return x;
}
scalar normal::minValue() const
Foam::scalar Foam::normal::minValue() const
{
return minValue_;
}
scalar normal::maxValue() const
Foam::scalar Foam::normal::maxValue() const
{
return maxValue_;
}
// ************************************************************************* //
}
// ************************************************************************* //

21
src/thermophysicalModels/pdfs/normal/normal.H
View File

@ -73,10 +73,11 @@ class normal
scalar range_;
public:
//- Runtime type information
TypeName("normal");
TypeName("normal");
// Constructors
@ -89,18 +90,20 @@ public:
);
// Destructor
~normal();
//- Destructor
virtual ~normal();
// Member Functions
scalar sample() const;
//- Sample the pdf
virtual scalar sample() const;
scalar minValue() const;
scalar maxValue() const;
//- Return the minimum value
virtual scalar minValue() const;
//- Return the maximum value
virtual scalar maxValue() const;
};
@ -110,10 +113,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//#include "normalI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

12
src/thermophysicalModels/pdfs/pdf/newPdf.C
View File

@ -28,10 +28,7 @@ License
// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
namespace Foam
{
autoPtr<Foam::pdf> Foam::pdf::New
Foam::autoPtr<Foam::pdf> Foam::pdf::New
(
const dictionary& dict,
Random& rndGen
@ -47,8 +44,8 @@ autoPtr<Foam::pdf> Foam::pdf::New
if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalErrorIn("pdf::New(const dictionary&, Random&)")
<< "unknown pdf type " << pdfType << endl << endl
<< "Valid pdf types are :" << endl
<< "unknown pdf type " << pdfType << nl << nl
<< "Valid pdf types are:" << nl
<< dictionaryConstructorTablePtr_->toc()
<< exit(FatalError);
}
@ -56,6 +53,5 @@ autoPtr<Foam::pdf> Foam::pdf::New
return autoPtr<pdf>(cstrIter()(dict, rndGen));
}
// ************************************************************************* //
}
// ************************************************************************* //

22
src/thermophysicalModels/pdfs/pdf/pdf.C
View File

@ -27,24 +27,16 @@ License
#include "pdf.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(pdf, 0);
defineRunTimeSelectionTable(pdf, dictionary);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
defineTypeNameAndDebug(pdf, 0);
defineRunTimeSelectionTable(pdf, dictionary);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
pdf::pdf
(
const dictionary& dict,
Random& rndGen
)
Foam::pdf::pdf(const dictionary& dict, Random& rndGen)
:
dict_(dict),
rndGen_(rndGen)
@ -53,10 +45,8 @@ pdf::pdf
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
pdf::~pdf()
Foam::pdf::~pdf()
{}
// ************************************************************************* //
} // end namespace Foam

54
src/thermophysicalModels/pdfs/pdf/pdf.H
View File

@ -91,44 +91,35 @@ protected:
public:
//-Runtime type information
TypeName("pdf");
TypeName("pdf");
// Declare runtime constructor selection table
declareRunTimeSelectionTable
//- Declare runtime constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
pdf,
dictionary,
(
autoPtr,
pdf,
dictionary,
(
const dictionary& dict,
Random& rndGen
),
(dict, rndGen)
);
const dictionary& dict,
Random& rndGen
),
(dict, rndGen)
);
// Constructors
//- Construct from dictionary
pdf
(
const dictionary& dict,
Random& rndGen
);
// Selectors
static autoPtr<pdf> New
(
const dictionary& dict,
Random& rndGen
);
pdf(const dictionary& dict, Random& rndGen);
// Destructor
//- Selector
static autoPtr<pdf> New(const dictionary& dict, Random& rndGen);
virtual ~pdf();
//- Destructor
virtual ~pdf();
// Member Functions
@ -136,7 +127,10 @@ public:
//- Sample the pdf
virtual scalar sample() const = 0;
//- Return the minimum value
virtual scalar minValue() const = 0;
//- Return the maximum value
virtual scalar maxValue() const = 0;
};
@ -147,10 +141,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//#include "pdfI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

35
src/thermophysicalModels/pdfs/uniform/uniform.C
View File

@ -27,29 +27,17 @@ License
#include "uniform.H"
#include "addToRunTimeSelectionTable.H"
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(uniform, 0);
addToRunTimeSelectionTable
(
pdf,
uniform,
dictionary
);
namespace Foam
{
defineTypeNameAndDebug(uniform, 0);
addToRunTimeSelectionTable(pdf, uniform, dictionary);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
uniform::uniform
(
const dictionary& dict,
Random& rndGen
)
Foam::uniform::uniform(const dictionary& dict, Random& rndGen)
:
pdf(dict, rndGen),
pdfDict_(dict.subDict(typeName + "PDF")),
@ -67,21 +55,22 @@ Foam::uniform::~uniform()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
scalar uniform::sample() const
Foam::scalar Foam::uniform::sample() const
{
return (minValue_ + rndGen_.scalar01()*range_);
}
scalar uniform::minValue() const
Foam::scalar Foam::uniform::minValue() const
{
return minValue_;
}
scalar uniform::maxValue() const
Foam::scalar Foam::uniform::maxValue() const
{
return maxValue_;
}
// ************************************************************************* //
}
// ************************************************************************* //

21
src/thermophysicalModels/pdfs/uniform/uniform.H
View File

@ -56,6 +56,7 @@ class uniform
// Private data
dictionary pdfDict_;
//- min and max values of the distribution
scalar minValue_;
scalar maxValue_;
@ -67,7 +68,7 @@ class uniform
public:
//- Runtime type information
TypeName("uniform");
TypeName("uniform");
// Constructors
@ -80,18 +81,20 @@ public:
);
// Destructor
~uniform();
//- Destructor
virtual ~uniform();
// Member Functions
scalar sample() const;
//- Sample the pdf
virtual scalar sample() const;
scalar minValue() const;
scalar maxValue() const;
//- Return the minimum value
virtual scalar minValue() const;
//- Return the maximum value
virtual scalar maxValue() const;
};
@ -101,10 +104,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//#include "uniformI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

16
src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.C
View File

@ -67,8 +67,7 @@ scalar mutRoughWallFunctionFvPatchScalarField::fnRough
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
mutRoughWallFunctionFvPatchScalarField::
mutRoughWallFunctionFvPatchScalarField
mutRoughWallFunctionFvPatchScalarField::mutRoughWallFunctionFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF
@ -80,8 +79,7 @@ mutRoughWallFunctionFvPatchScalarField
{}
mutRoughWallFunctionFvPatchScalarField::
mutRoughWallFunctionFvPatchScalarField
mutRoughWallFunctionFvPatchScalarField::mutRoughWallFunctionFvPatchScalarField
(
const mutRoughWallFunctionFvPatchScalarField& ptf,
const fvPatch& p,
@ -95,8 +93,7 @@ mutRoughWallFunctionFvPatchScalarField
{}
mutRoughWallFunctionFvPatchScalarField::
mutRoughWallFunctionFvPatchScalarField
mutRoughWallFunctionFvPatchScalarField::mutRoughWallFunctionFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
@ -109,8 +106,7 @@ mutRoughWallFunctionFvPatchScalarField
{}
mutRoughWallFunctionFvPatchScalarField::
mutRoughWallFunctionFvPatchScalarField
mutRoughWallFunctionFvPatchScalarField::mutRoughWallFunctionFvPatchScalarField
(
const mutRoughWallFunctionFvPatchScalarField& rwfpsf
)
@ -121,8 +117,7 @@ mutRoughWallFunctionFvPatchScalarField
{}
mutRoughWallFunctionFvPatchScalarField::
mutRoughWallFunctionFvPatchScalarField
mutRoughWallFunctionFvPatchScalarField::mutRoughWallFunctionFvPatchScalarField
(
const mutRoughWallFunctionFvPatchScalarField& rwfpsf,
const DimensionedField<scalar, volMesh>& iF
@ -221,6 +216,7 @@ tmp<scalarField> mutRoughWallFunctionFvPatchScalarField::calcMut() const
void mutRoughWallFunctionFvPatchScalarField::write(Ostream& os) const
{
fvPatchField<scalar>::write(os);
writeLocalEntries(os);
Cs_.writeEntry("Cs", os);
Ks_.writeEntry("Ks", os);
writeEntry("value", os);

275
src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C
View File

@ -39,6 +39,138 @@ namespace compressible
namespace RASModels
{
tmp<scalarField>
mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
(
const scalarField& magUp
) const
{
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
const scalarField& y = rasModel.y()[patchI];
const scalarField& muw = rasModel.mu().boundaryField()[patchI];
const fvPatchScalarField& rho = rasModel.rho().boundaryField()[patchI];
tmp<scalarField> 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 +255,24 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcMut() const
const RASModel& rasModel = db().lookupObject<RASModel>("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<scalarField> tyPlus = calcYPlus(magUp);
scalarField& yPlus = tyPlus();
tmp<scalarField> 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 +283,13 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcMut() const
tmp<scalarField>
mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::yPlus() const
{
notImplemented
(
"mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::yPlus()"
"const"
);
const label patchI = patch().index();
return tmp<scalarField>(NULL);
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
const scalarField magUp = mag(Uw.patchInternalField() - Uw);
return calcYPlus(magUp);
}
@ -283,6 +299,7 @@ void mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::write
) const
{
fixedValueFvPatchScalarField::write(os);
writeLocalEntries(os);
os.writeKeyword("roughnessHeight")
<< roughnessHeight_ << token::END_STATEMENT << nl;
os.writeKeyword("roughnessConstant")

3
src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H
View File

@ -74,6 +74,9 @@ protected:
// Protected member functions
//- Calculate yPLus
virtual tmp<scalarField> calcYPlus(const scalarField& magUp) const;
//- Calculate the turbulence viscosity
virtual tmp<scalarField> calcMut() const;

90
src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C
View File

@ -39,6 +39,49 @@ namespace compressible
namespace RASModels
{
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
tmp<scalarField>
mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcYPlus
(
const scalarField& magUp
) const
{
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("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<scalarField> 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::
@ -107,41 +150,22 @@ mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcMut() const
const RASModel& rasModel = db().lookupObject<RASModel>("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<scalarField> tyPlus = calcYPlus(magUp);
scalarField& yPlus = tyPlus();
tmp<scalarField> 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 +176,12 @@ mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcMut() const
tmp<scalarField>
mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::yPlus() const
{
notImplemented
(
"mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::yPlus() "
"const"
);
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
const scalarField magUp = mag(Uw.patchInternalField() - Uw);
return tmp<scalarField>(NULL);
return calcYPlus(magUp);
}
@ -168,8 +191,7 @@ void mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::write
) const
{
fvPatchField<scalar>::write(os);
os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
writeLocalEntries(os);
writeEntry("value", os);
}

3
src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.H
View File

@ -60,6 +60,9 @@ protected:
// Protected member functions
//- Calculate yPLus
virtual tmp<scalarField> calcYPlus(const scalarField& magUp) const;
//- Calculate the turbulence viscosity
virtual tmp<scalarField> calcMut() const;

11
src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasWallFunction/mutSpalartAllmarasWallFunctionFvPatchScalarField.C
View File

@ -39,7 +39,7 @@ namespace compressible
namespace RASModels
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
tmp<scalarField> mutSpalartAllmarasWallFunctionFvPatchScalarField::calcUTau
(
@ -173,11 +173,8 @@ mutSpalartAllmarasWallFunctionFvPatchScalarField::calcMut() const
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
const scalarField magGradU = mag(Uw.snGrad());
const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
const scalarField& muw = rasModel.mu().boundaryField()[patchI];
return max(0.0, rhow*sqr(calcUTau(magGradU))/magGradU - muw);
@ -191,11 +188,8 @@ mutSpalartAllmarasWallFunctionFvPatchScalarField::yPlus() const
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const scalarField& y = rasModel.y()[patchI];
const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
const scalarField& muw = rasModel.mu().boundaryField()[patchI];
return y*calcUTau(mag(Uw.snGrad()))/(muw/rhow);
@ -208,8 +202,7 @@ void mutSpalartAllmarasWallFunctionFvPatchScalarField::write
) const
{
fvPatchField<scalar>::write(os);
os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
writeLocalEntries(os);
writeEntry("value", os);
}

8
src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutWallFunction/mutWallFunctionFvPatchScalarField.C
View File

@ -189,10 +189,16 @@ tmp<scalarField> mutWallFunctionFvPatchScalarField::yPlus() const
void mutWallFunctionFvPatchScalarField::write(Ostream& os) const
{
fvPatchField<scalar>::write(os);
writeLocalEntries(os);
writeEntry("value", os);
}
void mutWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
{
os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
writeEntry("value", os);
}

3
src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutWallFunction/mutWallFunctionFvPatchScalarField.H
View File

@ -79,6 +79,9 @@ protected:
//- Calculate the turbulence viscosity
virtual tmp<scalarField> calcMut() const;
//- Write local wall function variables
virtual void writeLocalEntries(Ostream&) const;
public:

4
src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutRoughWallFunction/nutRoughWallFunctionFvPatchScalarField.C
View File

@ -223,9 +223,7 @@ tmp<scalarField> nutRoughWallFunctionFvPatchScalarField::calcNut() const
void nutRoughWallFunctionFvPatchScalarField::write(Ostream& os) const
{
fvPatchField<scalar>::write(os);
os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
writeLocalEntries(os);
Cs_.writeEntry("Cs", os);
Ks_.writeEntry("Ks", os);
writeEntry("value", os);

269
src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C
View File

@ -39,6 +39,134 @@ namespace incompressible
namespace RASModels
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
tmp<scalarField>
nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
(
const scalarField& magUp
) const
{
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
const scalarField& y = rasModel.y()[patchI];
const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
tmp<scalarField> 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;
yp = (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<RASModel>("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<scalarField> tyPlus = calcYPlus(magUp);
scalarField& yPlus = tyPlus();
tmp<scalarField> 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<scalarField>
nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::yPlus() const
{
notImplemented
(
"nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::yPlus()"
"const"
);
const label patchI = patch().index();
return tmp<scalarField>(NULL);
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
const scalarField magUp = mag(Uw.patchInternalField() - Uw);
return calcYPlus(magUp);
}
@ -273,15 +294,15 @@ void nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::write
Ostream& os
) const
{
fixedValueFvPatchScalarField::write(os);
os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
fvPatchField<scalar>::write(os);
writeLocalEntries(os);
os.writeKeyword("roughnessHeight")
<< roughnessHeight_ << token::END_STATEMENT << nl;
os.writeKeyword("roughnessConstant")
<< roughnessConstant_ << token::END_STATEMENT << nl;
os.writeKeyword("roughnessFudgeFactor")
<< roughnessFudgeFactor_ << token::END_STATEMENT << nl;
writeEntry("value", os);
}

3
src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H
View File

@ -72,6 +72,9 @@ class nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
// Protected member functions
//- Calculate yPLus
virtual tmp<scalarField> calcYPlus(const scalarField& magUp) const;
//- Calculate the turbulence viscosity
virtual tmp<scalarField> calcNut() const;

88
src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C
View File

@ -39,6 +39,48 @@ namespace incompressible
namespace RASModels
{
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
tmp<scalarField>
nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcYPlus
(
const scalarField& magUp
) const
{
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
const scalarField& y = rasModel.y()[patchI];
const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
tmp<scalarField> 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;
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 * * * * * * * * * * * * * * //
nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::
@ -107,37 +149,22 @@ nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcNut() const
const RASModel& rasModel = db().lookupObject<RASModel>("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<scalarField> tyPlus = calcYPlus(magUp);
scalarField& yPlus = tyPlus();
tmp<scalarField> 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<scalarField>
nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::yPlus() const
{
notImplemented
(
"nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::yPlus() "
"const"
);
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
const scalarField magUp = mag(Uw.patchInternalField() - Uw);
return tmp<scalarField>(NULL);
return calcYPlus(magUp);
}
@ -163,9 +189,9 @@ void nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::write
Ostream& os
) const
{
fixedValueFvPatchScalarField::write(os);
os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
fvPatchField<scalar>::write(os);
writeLocalEntries(os);
writeEntry("value", os);
}

3
src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.H
View File

@ -60,6 +60,9 @@ protected:
// Protected member functions
//- Calculate yPLus
virtual tmp<scalarField> calcYPlus(const scalarField& magUp) const;
//- Calculate the turbulence viscosity
virtual tmp<scalarField> calcNut() const;

13
src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasWallFunction/nutSpalartAllmarasWallFunctionFvPatchScalarField.C
View File

@ -51,8 +51,7 @@ tmp<scalarField> nutSpalartAllmarasWallFunctionFvPatchScalarField::calcUTau
const fvPatchVectorField& Uw =
rasModel.U().boundaryField()[patch().index()];
scalarField magUp = mag(Uw.patchInternalField() - Uw);
const scalarField magUp = mag(Uw.patchInternalField() - Uw);
const scalarField& nuw = rasModel.nu().boundaryField()[patch().index()];
const scalarField& nutw = *this;
@ -167,9 +166,7 @@ nutSpalartAllmarasWallFunctionFvPatchScalarField::calcNut() const
const RASModel& rasModel = db().lookupObject<RASModel>("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 +180,7 @@ nutSpalartAllmarasWallFunctionFvPatchScalarField::yPlus() const
const RASModel& rasModel = db().lookupObject<RASModel>("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;
@ -197,9 +192,9 @@ void nutSpalartAllmarasWallFunctionFvPatchScalarField::write
Ostream& os
) const
{
fixedValueFvPatchScalarField::write(os);
os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
fvPatchField<scalar>::write(os);
writeLocalEntries(os);
writeEntry("value", os);
}

11
src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.C
View File

@ -40,7 +40,7 @@ namespace incompressible
namespace RASModels
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void nutWallFunctionFvPatchScalarField::checkType()
{
@ -158,7 +158,6 @@ tmp<scalarField> nutWallFunctionFvPatchScalarField::calcNut() const
const scalar Cmu25 = pow(Cmu_, 0.25);
tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
scalarField& nutw = tnutw();
@ -197,10 +196,16 @@ tmp<scalarField> nutWallFunctionFvPatchScalarField::yPlus() const
void nutWallFunctionFvPatchScalarField::write(Ostream& os) const
{
fvPatchField<scalar>::write(os);
writeLocalEntries(os);
writeEntry("value", os);
}
void nutWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
{
os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
writeEntry("value", os);
}

3
src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.H
View File

@ -79,6 +79,9 @@ protected:
//- Calculate the turbulence viscosity
virtual tmp<scalarField> calcNut() const;
//- Write local wall function variables
virtual void writeLocalEntries(Ostream&) const;
public: