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 commit 'OpenCFD/master' into olesenm

Browse Source
This commit is contained in:
Mark Olesen
2009-07-24 09:19:46 +02:00
parent b095d78496 2ee961c502
commit 63dfefa5e8
34 changed files with 691 additions and 611 deletions
Download Patch File Download Diff File Expand all files Collapse all files

58
applications/utilities/mesh/manipulation/createBaffles/createBaffles.C
View File

@ -26,8 +26,8 @@ Description
Makes internal faces into boundary faces. Does not duplicate points, unlike
mergeOrSplitBaffles.
Note: if any coupled patch face is selected for baffling automatically
the opposite member is selected for baffling as well. Note that this
Note: if any coupled patch face is selected for baffling the opposite
member has to be selected for baffling as well. Note that this
is the same as repatching. This was added only for convenience so
you don't have to filter coupled boundary out of your set.
@ -128,6 +128,7 @@ int main(int argc, char *argv[])
argList::validArgs.append("faceZone");
argList::validArgs.append("patch");
argList::validOptions.insert("additionalPatches", "(patch2 .. patchN)");
argList::validOptions.insert("internalFacesOnly", "");
argList::validOptions.insert("overwrite", "");
# include "setRootCase.H"
@ -183,6 +184,12 @@ int main(int argc, char *argv[])
bool overwrite = args.optionFound("overwrite");
bool internalFacesOnly = args.optionFound("internalFacesOnly");
if (internalFacesOnly)
{
Info<< "Not converting faces on non-coupled patches." << nl << endl;
}
// Read objects in time directory
@ -234,7 +241,21 @@ int main(int argc, char *argv[])
// guarantees that when e.g. creating a cyclic all faces from one
// side come first and faces from the other side next.
// Whether first use of face (modify) or consecutive (add)
PackedBoolList modifiedFace(mesh.nFaces());
// Never modify coupled faces
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
forAll(pp, i)
{
modifiedFace[pp.start()+i] = 1;
}
}
}
label nModified = 0;
forAll(newPatches, i)
{
@ -281,6 +302,8 @@ int main(int argc, char *argv[])
modifiedFace // modify or add status
);
}
nModified++;
}
}
@ -333,16 +356,27 @@ int main(int argc, char *argv[])
// Modify any boundary faces
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Normal boundary:
// - move to new patch. Might already be back-to-back baffle
// you want to add cyclic to. Do warn though.
//
// Processor boundary:
// - do not move to cyclic
// - add normal patches though.
// For warning once per patch.
labelHashSet patchWarned;
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (patches[newPatchI].coupled() && pp.coupled())
if (pp.coupled() && patches[newPatchI].coupled())
{
// Do not allow coupled faces to be moved to different coupled
// patches.
}
else
else if (pp.coupled() || !internalFacesOnly)
{
forAll(pp, i)
{
@ -352,6 +386,19 @@ int main(int argc, char *argv[])
if (zoneFaceI != -1)
{
if (patchWarned.insert(patchI))
{
WarningIn(args.executable())
<< "Found boundary face (in patch " << pp.name()
<< ") in faceZone " << fZone.name()
<< " to convert to baffle patch "
<< patches[newPatchI].name()
<< endl
<< " Run with -internalFacesOnly option"
<< " if you don't wish to convert"
<< " boundary faces." << endl;
}
modifyOrAddFace
(
meshMod,
@ -364,6 +411,7 @@ int main(int argc, char *argv[])
fZone.flipMap()[zoneFaceI], // face flip in zone
modifiedFace // modify or add status
);
nModified++;
}
}
}
@ -371,7 +419,7 @@ int main(int argc, char *argv[])
}
Info<< "Converted " << returnReduce(modifiedFace.count(), sumOp<label>())
Info<< "Converted " << returnReduce(nModified, sumOp<label>())
<< " faces into boundary faces on patch " << patchName << nl << endl;
if (!overwrite)

4
applications/utilities/postProcessing/dataConversion/foamToEnsightParts/foamToEnsightParts.C
View File

@ -36,8 +36,8 @@ Usage
@param -ascii \n
Write Ensight data in ASCII format instead of "C Binary"
@param -zeroTime \n
Include the often incomplete initial conditions.
@param -noZero \n
Exclude the often incomplete initial conditions.
@param -index \<start\>\n
Ignore the time index contained in the time file and use a

7
src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
View File

@ -1112,6 +1112,13 @@ bool Foam::cyclicPolyPatch::order
return false;
}
if (pp.size()&1)
{
FatalErrorIn("cyclicPolyPatch::order(..)")
<< "Size of cyclic " << name() << " should be a multiple of 2"
<< ". It is " << pp.size() << abort(FatalError);
}
label halfSize = pp.size()/2;
// Supplied primitivePatch already with new points.

5
src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
View File

@ -191,7 +191,10 @@ void Foam::polyTopoChange::countMap
}
Foam::labelHashSet Foam::polyTopoChange::getSetIndices(const PackedBoolList& lst)
Foam::labelHashSet Foam::polyTopoChange::getSetIndices
(
const PackedBoolList& lst
)
{
labelHashSet values(lst.count());
forAll(lst, i)

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
// ************************************************************************* //

2
src/topoChangerFvMesh/Make/files
View File

@ -5,6 +5,8 @@ linearValveFvMesh/linearValveFvMesh.C
linearValveLayersFvMesh/linearValveLayersFvMesh.C
*/
movingConeTopoFvMesh/movingConeTopoFvMesh.C
/*
mixerFvMesh/mixerFvMesh.C
*/
LIB = $(FOAM_LIBBIN)/libtopoChangerFvMesh

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: