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

This commit is contained in:
mattijs
2013-05-31 17:46:44 +01:00
20 changed files with 2365 additions and 609 deletions

View File

@ -284,6 +284,16 @@ Foam::tmp<Foam::Field<Type> > Foam::cyclicACMIFvPatchField<Type>::snGrad
}
template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::updateCoeffs()
{
const scalarField& mask = cyclicACMIPatch_.cyclicACMIPatch().mask();
const fvPatchField<Type>& npf = nonOverlapPatchField();
const_cast<fvPatchField<Type>&>(npf).updateCoeffs(mask);
}
template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::evaluate
(
@ -376,6 +386,20 @@ Foam::cyclicACMIFvPatchField<Type>::gradientBoundaryCoeffs() const
}
template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::manipulateMatrix
(
fvMatrix<Type>& matrix
)
{
// blend contrubutions from the coupled and non-overlap patches
const fvPatchField<Type>& npf = nonOverlapPatchField();
const scalarField& mask = cyclicACMIPatch_.cyclicACMIPatch().mask();
const_cast<fvPatchField<Type>&>(npf).manipulateMatrix(matrix, mask);
}
template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::write(Ostream& os) const
{

View File

@ -178,12 +178,34 @@ public:
//- Return reference to non-overlapping patchField
const fvPatchField<Type>& nonOverlapPatchField() const;
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
scalarField& result,
const scalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
) const;
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
Field<Type>&,
const Field<Type>&,
const scalarField&,
const Pstream::commsTypes commsType
) const;
//- Return patch-normal gradient
virtual tmp<Field<Type> > snGrad
(
const scalarField& deltaCoeffs
) const;
//- Update the coefficients associated with the patch field
void updateCoeffs();
//- Evaluate the patch field
virtual void evaluate
(
@ -226,24 +248,8 @@ public:
// evaluation of the gradient of this patchField
virtual tmp<Field<Type> > gradientBoundaryCoeffs() const;
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
scalarField& result,
const scalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
) const;
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
Field<Type>&,
const Field<Type>&,
const scalarField&,
const Pstream::commsTypes commsType
) const;
//- Manipulate matrix
virtual void manipulateMatrix(fvMatrix<Type>& matrix);
// Cyclic AMI coupled interface functions

View File

@ -42,6 +42,7 @@ Foam::fvPatchField<Type>::fvPatchField
patch_(p),
internalField_(iF),
updated_(false),
manipulatedMatrix_(false),
patchType_(word::null)
{}
@ -58,6 +59,7 @@ Foam::fvPatchField<Type>::fvPatchField
patch_(p),
internalField_(iF),
updated_(false),
manipulatedMatrix_(false),
patchType_(word::null)
{}
@ -75,6 +77,7 @@ Foam::fvPatchField<Type>::fvPatchField
patch_(p),
internalField_(iF),
updated_(false),
manipulatedMatrix_(false),
patchType_(ptf.patchType_)
{}
@ -92,6 +95,7 @@ Foam::fvPatchField<Type>::fvPatchField
patch_(p),
internalField_(iF),
updated_(false),
manipulatedMatrix_(false),
patchType_(dict.lookupOrDefault<word>("patchType", word::null))
{
if (dict.found("value"))
@ -133,6 +137,7 @@ Foam::fvPatchField<Type>::fvPatchField
patch_(ptf.patch_),
internalField_(ptf.internalField_),
updated_(false),
manipulatedMatrix_(false),
patchType_(ptf.patchType_)
{}
@ -148,6 +153,7 @@ Foam::fvPatchField<Type>::fvPatchField
patch_(ptf.patch_),
internalField_(iF),
updated_(false),
manipulatedMatrix_(false),
patchType_(ptf.patchType_)
{}
@ -267,6 +273,28 @@ void Foam::fvPatchField<Type>::rmap
}
template<class Type>
void Foam::fvPatchField<Type>::updateCoeffs()
{
updated_ = true;
}
template<class Type>
void Foam::fvPatchField<Type>::updateCoeffs(const scalarField& weights)
{
if (!updated_)
{
updateCoeffs();
Field<Type>& fld = *this;
fld *= weights;
updated_ = true;
}
}
template<class Type>
void Foam::fvPatchField<Type>::evaluate(const Pstream::commsTypes)
{
@ -276,13 +304,25 @@ void Foam::fvPatchField<Type>::evaluate(const Pstream::commsTypes)
}
updated_ = false;
manipulatedMatrix_ = false;
}
template<class Type>
void Foam::fvPatchField<Type>::manipulateMatrix(fvMatrix<Type>& matrix)
{
// do nothing
manipulatedMatrix_ = true;
}
template<class Type>
void Foam::fvPatchField<Type>::manipulateMatrix
(
fvMatrix<Type>& matrix,
const scalarField& weights
)
{
manipulatedMatrix_ = true;
}

View File

@ -93,6 +93,10 @@ class fvPatchField
// the construction of the matrix
bool updated_;
//- Update index used so that manipulateMatrix is called only once
// during the construction of the matrix
bool manipulatedMatrix_;
//- Optional patch type, used to allow specified boundary conditions
// to be applied to constraint patches by providing the constraint
// patch type as 'patchType'
@ -327,6 +331,12 @@ public:
return updated_;
}
//- Return true if the matrix has already been manipulated
bool manipulatedMatrix() const
{
return manipulatedMatrix_;
}
// Mapping functions
@ -365,10 +375,12 @@ public:
//- Update the coefficients associated with the patch field
// Sets Updated to true
virtual void updateCoeffs()
{
updated_ = true;
}
virtual void updateCoeffs();
//- Update the coefficients associated with the patch field
// and apply weight field
// Sets Updated to true
virtual void updateCoeffs(const scalarField& weights);
//- Return internal field next to patch as patch field
virtual tmp<Field<Type> > patchInternalField() const;
@ -479,6 +491,13 @@ public:
//- Manipulate matrix
virtual void manipulateMatrix(fvMatrix<Type>& matrix);
//- Manipulate matrix with given weights
virtual void manipulateMatrix
(
fvMatrix<Type>& matrix,
const scalarField& weights
);
// I-O

View File

@ -342,6 +342,18 @@ Foam::label Foam::cyclicACMIPolyPatch::nonOverlapPatchID() const
<< exit(FatalError);
}
if (nonOverlapPatchID_ < index())
{
FatalErrorIn("cyclicPolyAMIPatch::neighbPatchID() const")
<< "Boundary ordering error: " << type()
<< " patch must be defined prior to its non-overlapping patch"
<< nl
<< type() << " patch: " << name() << ", ID:" << index() << nl
<< "Non-overlap patch: " << nonOverlapPatchName_
<< ", ID:" << nonOverlapPatchID_ << nl
<< exit(FatalError);
}
const polyPatch& noPp = this->boundaryMesh()[nonOverlapPatchID_];
bool ok = true;

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,6 +25,7 @@ License
#include "backwardsCompatibilityWallFunctions.H"
#include "volFields.H"
#include "calculatedFvPatchField.H"
#include "alphatWallFunctionFvPatchScalarField.H"
#include "mutkWallFunctionFvPatchScalarField.H"

View File

@ -55,6 +55,67 @@ scalar epsilonLowReWallFunctionFvPatchScalarField::yPlusLam
}
void epsilonLowReWallFunctionFvPatchScalarField::calculate
(
const turbulenceModel& turbulence,
const List<scalar>& cornerWeights,
const fvPatch& patch,
scalarField& G,
scalarField& epsilon
)
{
const label patchI = patch.index();
const scalarField& y = turbulence.y()[patchI];
const scalar Cmu25 = pow025(Cmu_);
const scalar Cmu75 = pow(Cmu_, 0.75);
const tmp<volScalarField> tk = turbulence.k();
const volScalarField& k = tk();
const tmp<volScalarField> tmu = turbulence.mu();
const scalarField& muw = tmu().boundaryField()[patchI];
const tmp<volScalarField> tmut = turbulence.mut();
const volScalarField& mut = tmut();
const scalarField& mutw = mut.boundaryField()[patchI];
const scalarField& rhow = turbulence.rho().boundaryField()[patchI];
const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
const scalarField magGradUw(mag(Uw.snGrad()));
// Set epsilon and G
forAll(mutw, faceI)
{
label cellI = patch.faceCells()[faceI];
scalar yPlus = Cmu25*sqrt(k[cellI])*y[faceI]/muw[faceI]/rhow[faceI];
scalar w = cornerWeights[faceI];
if (yPlus > yPlusLam_)
{
epsilon[cellI] = w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
}
else
{
epsilon[cellI] =
w*2.0*k[cellI]*muw[faceI]/rhow[faceI]/sqr(y[faceI]);
}
G[cellI] =
w
*(mutw[faceI] + muw[faceI])
*magGradUw[faceI]
*Cmu25*sqrt(k[cellI])
/(kappa_*y[faceI]);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
epsilonLowReWallFunctionFvPatchScalarField::
@ -119,84 +180,6 @@ epsilonLowReWallFunctionFvPatchScalarField
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void epsilonLowReWallFunctionFvPatchScalarField::updateCoeffs()
{
if (updated())
{
return;
}
const label patchI = patch().index();
const turbulenceModel& turbulence =
db().lookupObject<turbulenceModel>("turbulenceModel");
const scalarField& y = turbulence.y()[patchI];
volScalarField& G =
const_cast<volScalarField&>
(
db().lookupObject<volScalarField>
(
turbulence.GName()
)
);
DimensionedField<scalar, volMesh>& epsilon =
const_cast<DimensionedField<scalar, volMesh>&>
(
dimensionedInternalField()
);
const tmp<volScalarField> tk = turbulence.k();
const volScalarField& k = tk();
const tmp<volScalarField> tmu = turbulence.mu();
const scalarField& muw = tmu().boundaryField()[patchI];
const tmp<volScalarField> tmut = turbulence.mut();
const volScalarField& mut = tmut();
const scalarField& mutw = mut.boundaryField()[patchI];
const scalarField& rhow = turbulence.rho().boundaryField()[patchI];
const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
const scalarField magGradUw(mag(Uw.snGrad()));
const scalar Cmu25 = pow025(Cmu_);
const scalar Cmu75 = pow(Cmu_, 0.75);
// Set epsilon and G
forAll(mutw, faceI)
{
label faceCellI = patch().faceCells()[faceI];
scalar yPlus = Cmu25*sqrt(k[faceCellI])*y[faceI]/muw[faceI]/rhow[faceI];
if (yPlus > yPlusLam_)
{
epsilon[faceCellI] = Cmu75*pow(k[faceCellI], 1.5)/(kappa_*y[faceI]);
}
else
{
epsilon[faceCellI] =
2.0*k[faceCellI]*muw[faceI]/rhow[faceI]/sqr(y[faceI]);
}
G[faceCellI] =
(mutw[faceI] + muw[faceI])
*magGradUw[faceI]
*Cmu25*sqrt(k[faceCellI])
/(kappa_*y[faceI]);
}
fixedInternalValueFvPatchField<scalar>::updateCoeffs();
// TODO: perform averaging for cells sharing more than one boundary face
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makePatchTypeField

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -97,6 +97,16 @@ protected:
//- Calculate the Y+ at the edge of the laminar sublayer
scalar yPlusLam(const scalar kappa, const scalar E);
//- Calculate the epsilon and G
virtual void calculate
(
const turbulenceModel& turbulence,
const List<scalar>& cornerWeights,
const fvPatch& patch,
scalarField& G,
scalarField& epsilon
);
public:
@ -165,14 +175,6 @@ public:
new epsilonLowReWallFunctionFvPatchScalarField(*this, iF)
);
}
// Member functions
// Evaluation functions
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
};

View File

@ -26,10 +26,10 @@ License
#include "epsilonWallFunctionFvPatchScalarField.H"
#include "compressible/turbulenceModel/turbulenceModel.H"
#include "fvPatchFieldMapper.H"
#include "fvMatrix.H"
#include "volFields.H"
#include "addToRunTimeSelectionTable.H"
#include "mutWallFunctionFvPatchScalarField.H"
#include "wallFvPatch.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -62,6 +62,193 @@ void epsilonWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
}
void epsilonWallFunctionFvPatchScalarField::setMaster()
{
if (master_ != -1)
{
return;
}
const volScalarField& epsilon =
static_cast<const volScalarField&>(this->dimensionedInternalField());
const volScalarField::GeometricBoundaryField& bf = epsilon.boundaryField();
label master = -1;
forAll(bf, patchI)
{
if (isA<epsilonWallFunctionFvPatchScalarField>(bf[patchI]))
{
epsilonWallFunctionFvPatchScalarField& epf = epsilonPatch(patchI);
if (master == -1)
{
master = patchI;
}
epf.master() = master;
}
}
}
void epsilonWallFunctionFvPatchScalarField::createAveragingWeights()
{
if (initialised_)
{
return;
}
const volScalarField& epsilon =
static_cast<const volScalarField&>(this->dimensionedInternalField());
const volScalarField::GeometricBoundaryField& bf = epsilon.boundaryField();
const fvMesh& mesh = epsilon.mesh();
volScalarField weights
(
IOobject
(
"weights",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false // do not register
),
mesh,
dimensionedScalar("zero", dimless, 0.0)
);
DynamicList<label> epsilonPatches(bf.size());
forAll(bf, patchI)
{
if (isA<epsilonWallFunctionFvPatchScalarField>(bf[patchI]))
{
epsilonPatches.append(patchI);
const labelUList& faceCells = bf[patchI].patch().faceCells();
forAll(faceCells, i)
{
label cellI = faceCells[i];
weights[cellI]++;
}
}
}
cornerWeights_.setSize(bf.size());
forAll(epsilonPatches, i)
{
label patchI = epsilonPatches[i];
const fvPatchField& wf = weights.boundaryField()[patchI];
cornerWeights_[patchI] = 1.0/wf.patchInternalField();
}
G_.setSize(dimensionedInternalField().size(), 0.0);
epsilon_.setSize(dimensionedInternalField().size(), 0.0);
initialised_ = true;
}
epsilonWallFunctionFvPatchScalarField&
epsilonWallFunctionFvPatchScalarField::epsilonPatch(const label patchI)
{
const volScalarField& epsilon =
static_cast<const volScalarField&>(this->dimensionedInternalField());
const volScalarField::GeometricBoundaryField& bf = epsilon.boundaryField();
const epsilonWallFunctionFvPatchScalarField& epf =
refCast<const epsilonWallFunctionFvPatchScalarField>(bf[patchI]);
return const_cast<epsilonWallFunctionFvPatchScalarField&>(epf);
}
void epsilonWallFunctionFvPatchScalarField::calculateTurbulenceFields
(
const turbulenceModel& turbulence,
scalarField& G0,
scalarField& epsilon0
)
{
// accumulate all of the G and epsilon contributions
forAll(cornerWeights_, patchI)
{
if (!cornerWeights_[patchI].empty())
{
epsilonWallFunctionFvPatchScalarField& epf = epsilonPatch(patchI);
const List<scalar>& w = cornerWeights_[patchI];
epf.calculate(turbulence, w, epf.patch(), G0, epsilon0);
}
}
// apply zero-gradient condition for epsilon
forAll(cornerWeights_, patchI)
{
if (!cornerWeights_[patchI].empty())
{
epsilonWallFunctionFvPatchScalarField& epf = epsilonPatch(patchI);
epf == scalarField(epsilon0, epf.patch().faceCells());
}
}
}
void epsilonWallFunctionFvPatchScalarField::calculate
(
const turbulenceModel& turbulence,
const List<scalar>& cornerWeights,
const fvPatch& patch,
scalarField& G,
scalarField& epsilon
)
{
const label patchI = patch.index();
const scalarField& y = turbulence.y()[patchI];
const scalar Cmu25 = pow025(Cmu_);
const scalar Cmu75 = pow(Cmu_, 0.75);
const tmp<volScalarField> tk = turbulence.k();
const volScalarField& k = tk();
const tmp<volScalarField> tmu = turbulence.mu();
const scalarField& muw = tmu().boundaryField()[patchI];
const tmp<volScalarField> tmut = turbulence.mut();
const volScalarField& mut = tmut();
const scalarField& mutw = mut.boundaryField()[patchI];
const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
const scalarField magGradUw(mag(Uw.snGrad()));
// Set epsilon and G
forAll(mutw, faceI)
{
label cellI = patch.faceCells()[faceI];
scalar w = cornerWeights[faceI];
epsilon[cellI] = w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
G[cellI] =
w
*(mutw[faceI] + muw[faceI])
*magGradUw[faceI]
*Cmu25*sqrt(k[cellI])
/(kappa_*y[faceI]);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
@ -70,10 +257,15 @@ epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
const DimensionedField<scalar, volMesh>& iF
)
:
fixedInternalValueFvPatchField<scalar>(p, iF),
fixedValueFvPatchField<scalar>(p, iF),
Cmu_(0.09),
kappa_(0.41),
E_(9.8)
E_(9.8),
G_(),
epsilon_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
@ -87,10 +279,15 @@ epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
const fvPatchFieldMapper& mapper
)
:
fixedInternalValueFvPatchField<scalar>(ptf, p, iF, mapper),
fixedValueFvPatchField<scalar>(ptf, p, iF, mapper),
Cmu_(ptf.Cmu_),
kappa_(ptf.kappa_),
E_(ptf.E_)
E_(ptf.E_),
G_(),
epsilon_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
@ -103,10 +300,15 @@ epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
const dictionary& dict
)
:
fixedInternalValueFvPatchField<scalar>(p, iF, dict),
fixedValueFvPatchField<scalar>(p, iF, dict),
Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
E_(dict.lookupOrDefault<scalar>("E", 9.8))
E_(dict.lookupOrDefault<scalar>("E", 9.8)),
G_(),
epsilon_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
@ -117,10 +319,15 @@ epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
const epsilonWallFunctionFvPatchScalarField& ewfpsf
)
:
fixedInternalValueFvPatchField<scalar>(ewfpsf),
fixedValueFvPatchField<scalar>(ewfpsf),
Cmu_(ewfpsf.Cmu_),
kappa_(ewfpsf.kappa_),
E_(ewfpsf.E_)
E_(ewfpsf.E_),
G_(),
epsilon_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
@ -132,10 +339,15 @@ epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
const DimensionedField<scalar, volMesh>& iF
)
:
fixedInternalValueFvPatchField<scalar>(ewfpsf, iF),
fixedValueFvPatchField<scalar>(ewfpsf, iF),
Cmu_(ewfpsf.Cmu_),
kappa_(ewfpsf.kappa_),
E_(ewfpsf.E_)
E_(ewfpsf.E_),
G_(),
epsilon_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
@ -143,6 +355,38 @@ epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
scalarField& epsilonWallFunctionFvPatchScalarField::G(bool init)
{
if (patch().index() == master_)
{
if (init)
{
G_ = 0.0;
}
return G_;
}
return epsilonPatch(master_).G();
}
scalarField& epsilonWallFunctionFvPatchScalarField::epsilon(bool init)
{
if (patch().index() == master_)
{
if (init)
{
epsilon_ = 0.0;
}
return epsilon_;
}
return epsilonPatch(master_).epsilon(init);
}
void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
{
if (updated())
@ -150,76 +394,168 @@ void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
return;
}
const label patchI = patch().index();
const turbulenceModel& turbulence =
db().lookupObject<turbulenceModel>("turbulenceModel");
db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
const scalar Cmu25 = pow025(Cmu_);
const scalar Cmu75 = pow(Cmu_, 0.75);
setMaster();
const scalarField& y = turbulence.y()[patchI];
volScalarField& G =
const_cast<volScalarField&>
(
db().lookupObject<volScalarField>
(
turbulence.GName()
)
);
DimensionedField<scalar, volMesh>& epsilon =
const_cast<DimensionedField<scalar, volMesh>&>
(
dimensionedInternalField()
);
const tmp<volScalarField> tk = turbulence.k();
const volScalarField& k = tk();
const scalarField& muw = turbulence.mu().boundaryField()[patchI];
const tmp<volScalarField> tmut = turbulence.mut();
const volScalarField& mut = tmut();
const scalarField& mutw = mut.boundaryField()[patchI];
const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
const scalarField magGradUw(mag(Uw.snGrad()));
// Set epsilon and G
forAll(mutw, faceI)
if (patch().index() == master_)
{
label faceCellI = patch().faceCells()[faceI];
epsilon[faceCellI] = Cmu75*pow(k[faceCellI], 1.5)/(kappa_*y[faceI]);
G[faceCellI] =
(mutw[faceI] + muw[faceI])
*magGradUw[faceI]
*Cmu25*sqrt(k[faceCellI])
/(kappa_*y[faceI]);
createAveragingWeights();
calculateTurbulenceFields(turbulence, G(true), epsilon(true));
}
fixedInternalValueFvPatchField<scalar>::updateCoeffs();
const scalarField& G0 = this->G();
const scalarField& epsilon0 = this->epsilon();
// TODO: perform averaging for cells sharing more than one boundary face
typedef DimensionedField<scalar, volMesh> FieldType;
FieldType& G =
const_cast<FieldType&>
(
db().lookupObject<FieldType>(turbulence.GName())
);
FieldType& epsilon = const_cast<FieldType&>(dimensionedInternalField());
forAll(*this, faceI)
{
label cellI = patch().faceCells()[faceI];
G[cellI] = G0[cellI];
epsilon[cellI] = epsilon0[cellI];
}
fvPatchField<scalar>::updateCoeffs();
}
void epsilonWallFunctionFvPatchScalarField::evaluate
void epsilonWallFunctionFvPatchScalarField::updateCoeffs
(
const Pstream::commsTypes commsType
const scalarField& weights
)
{
fixedInternalValueFvPatchField<scalar>::evaluate(commsType);
if (updated())
{
return;
}
const turbulenceModel& turbulence =
db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
setMaster();
if (patch().index() == master_)
{
createAveragingWeights();
calculateTurbulenceFields(turbulence, G(true), epsilon(true));
}
const scalarField& G0 = this->G();
const scalarField& epsilon0 = this->epsilon();
typedef DimensionedField<scalar, volMesh> FieldType;
FieldType& G =
const_cast<FieldType&>
(
db().lookupObject<FieldType>(turbulence.GName())
);
FieldType& epsilon = const_cast<FieldType&>(dimensionedInternalField());
// only set the values if the weights are < 1 - tolerance
forAll(weights, faceI)
{
scalar w = weights[faceI];
if (w < 1.0 - 1e-6)
{
label cellI = patch().faceCells()[faceI];
G[cellI] = w*G[cellI] + (1.0 - w)*G0[cellI];
epsilon[cellI] = w*epsilon[cellI] + (1.0 - w)*epsilon0[cellI];
}
}
fvPatchField<scalar>::updateCoeffs();
}
void epsilonWallFunctionFvPatchScalarField::manipulateMatrix
(
fvMatrix<scalar>& matrix
)
{
if (manipulatedMatrix())
{
return;
}
matrix.setValues(patch().faceCells(), patchInternalField());
fvPatchField<scalar>::manipulateMatrix(matrix);
}
void epsilonWallFunctionFvPatchScalarField::manipulateMatrix
(
fvMatrix<scalar>& matrix,
const Field<scalar>& weights
)
{
if (manipulatedMatrix())
{
return;
}
// filter weights so that we only apply the constraint where the
// weight > SMALL
DynamicList<label> constraintCells(weights.size());
DynamicList<scalar> constraintEpsilon(weights.size());
const labelUList& faceCells = patch().faceCells();
const DimensionedField<scalar, volMesh>& epsilon
= dimensionedInternalField();
label nConstrainedCells = 0;
forAll(weights, faceI)
{
// only set the values if the weights are < 1 - tolerance
if (weights[faceI] < (1.0 - 1e-6))
{
nConstrainedCells++;
label cellI = faceCells[faceI];
constraintCells.append(cellI);
constraintEpsilon.append(epsilon[cellI]);
}
}
if (debug)
{
Pout<< "Patch: " << patch().name()
<< ": number of constrained cells = " << nConstrainedCells
<< " out of " << patch().size()
<< endl;
}
matrix.setValues
(
constraintCells,
scalarField(constraintEpsilon.xfer())
);
fvPatchField<scalar>::manipulateMatrix(matrix);
}
void epsilonWallFunctionFvPatchScalarField::write(Ostream& os) const
{
fixedInternalValueFvPatchField<scalar>::write(os);
fixedValueFvPatchField<scalar>::write(os);
writeLocalEntries(os);
writeEntry("value", os);
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -71,7 +71,7 @@ SourceFiles
#ifndef compressibleEpsilonWallFunctionFvPatchScalarField_H
#define compressibleEpsilonWallFunctionFvPatchScalarField_H
#include "fixedInternalValueFvPatchField.H"
#include "fixedValueFvPatchField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -80,13 +80,15 @@ namespace Foam
namespace compressible
{
class turbulenceModel;
/*---------------------------------------------------------------------------*\
Class epsilonWallFunctionFvPatchScalarField Declaration
\*---------------------------------------------------------------------------*/
class epsilonWallFunctionFvPatchScalarField
:
public fixedInternalValueFvPatchField<scalar>
public fixedValueFvPatchField<scalar>
{
protected:
@ -101,6 +103,21 @@ protected:
//- E coefficient
scalar E_;
//- Local copy of turbulence G field
scalarField G_;
//- Local copy of turbulence epsilon field
scalarField epsilon_;
//- Initialised flag
bool initialised_;
//- Master patch ID
label master_;
//- List of averaging corner weights
List<List<scalar> > cornerWeights_;
// Protected Member Functions
@ -110,6 +127,44 @@ protected:
//- Write local wall function variables
virtual void writeLocalEntries(Ostream&) const;
//- Set the master patch - master is responsible for updating all
// wall function patches
virtual void setMaster();
//- Create the averaging weights for cells which are bounded by
// multiple wall function faces
virtual void createAveragingWeights();
//- Helper function to return non-const access to an epsilon patch
virtual epsilonWallFunctionFvPatchScalarField& epsilonPatch
(
const label patchI
);
//- Main driver to calculate the turbulence fields
virtual void calculateTurbulenceFields
(
const turbulenceModel& turbulence,
scalarField& G0,
scalarField& epsilon0
);
//- Calculate the epsilon and G
virtual void calculate
(
const turbulenceModel& turbulence,
const List<scalar>& cornerWeights,
const fvPatch& patch,
scalarField& G,
scalarField& epsilon
);
//- Return non-const access to the master patch ID
virtual label& master()
{
return master_;
}
public:
@ -180,15 +235,39 @@ public:
}
//- Destructor
virtual ~epsilonWallFunctionFvPatchScalarField()
{}
// Member functions
// Access
//- Return non-const access to the master's G field
scalarField& G(bool init = false);
//- Return non-const access to the master's epsilon field
scalarField& epsilon(bool init = false);
// Evaluation functions
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Evaluate the patchField
virtual void evaluate(const Pstream::commsTypes);
//- Update the coefficients associated with the patch field
virtual void updateCoeffs(const scalarField& weights);
//- Manipulate matrix
virtual void manipulateMatrix(fvMatrix<scalar>& matrix);
//- Manipulate matrix with given weights
virtual void manipulateMatrix
(
fvMatrix<scalar>& matrix,
const scalarField& weights
);
// I-O

View File

@ -26,10 +26,11 @@ License
#include "omegaWallFunctionFvPatchScalarField.H"
#include "compressible/turbulenceModel/turbulenceModel.H"
#include "fvPatchFieldMapper.H"
#include "fvMatrix.H"
#include "volFields.H"
#include "addToRunTimeSelectionTable.H"
#include "mutWallFunctionFvPatchScalarField.H"
#include "wallFvPatch.H"
#include "mutWallFunctionFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -63,6 +64,198 @@ void omegaWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
}
void omegaWallFunctionFvPatchScalarField::setMaster()
{
if (master_ != -1)
{
return;
}
const volScalarField& omega =
static_cast<const volScalarField&>(this->dimensionedInternalField());
const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
label master = -1;
forAll(bf, patchI)
{
if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchI]))
{
omegaWallFunctionFvPatchScalarField& epf = omegaPatch(patchI);
if (master == -1)
{
master = patchI;
}
epf.master() = master;
}
}
}
void omegaWallFunctionFvPatchScalarField::createAveragingWeights()
{
if (initialised_)
{
return;
}
const volScalarField& omega =
static_cast<const volScalarField&>(this->dimensionedInternalField());
const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
const fvMesh& mesh = omega.mesh();
volScalarField weights
(
IOobject
(
"weights",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false // do not register
),
mesh,
dimensionedScalar("zero", dimless, 0.0)
);
DynamicList<label> omegaPatches(bf.size());
forAll(bf, patchI)
{
if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchI]))
{
omegaPatches.append(patchI);
const labelUList& faceCells = bf[patchI].patch().faceCells();
forAll(faceCells, i)
{
label cellI = faceCells[i];
weights[cellI]++;
}
}
}
cornerWeights_.setSize(bf.size());
forAll(omegaPatches, i)
{
label patchI = omegaPatches[i];
const fvPatchField& wf = weights.boundaryField()[patchI];
cornerWeights_[patchI] = 1.0/wf.patchInternalField();
}
G_.setSize(dimensionedInternalField().size(), 0.0);
omega_.setSize(dimensionedInternalField().size(), 0.0);
initialised_ = true;
}
omegaWallFunctionFvPatchScalarField&
omegaWallFunctionFvPatchScalarField::omegaPatch(const label patchI)
{
const volScalarField& omega =
static_cast<const volScalarField&>(this->dimensionedInternalField());
const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
const omegaWallFunctionFvPatchScalarField& epf =
refCast<const omegaWallFunctionFvPatchScalarField>(bf[patchI]);
return const_cast<omegaWallFunctionFvPatchScalarField&>(epf);
}
void omegaWallFunctionFvPatchScalarField::calculateTurbulenceFields
(
const turbulenceModel& turbulence,
scalarField& G0,
scalarField& omega0
)
{
// accumulate all of the G and omega contributions
forAll(cornerWeights_, patchI)
{
if (!cornerWeights_[patchI].empty())
{
omegaWallFunctionFvPatchScalarField& epf = omegaPatch(patchI);
const List<scalar>& w = cornerWeights_[patchI];
epf.calculate(turbulence, w, epf.patch(), G0, omega0);
}
}
// apply zero-gradient condition for omega
forAll(cornerWeights_, patchI)
{
if (!cornerWeights_[patchI].empty())
{
omegaWallFunctionFvPatchScalarField& epf = omegaPatch(patchI);
epf == scalarField(omega0, epf.patch().faceCells());
}
}
}
void omegaWallFunctionFvPatchScalarField::calculate
(
const turbulenceModel& turbulence,
const List<scalar>& cornerWeights,
const fvPatch& patch,
scalarField& G,
scalarField& omega
)
{
const label patchI = patch.index();
const scalarField& y = turbulence.y()[patchI];
const scalar Cmu25 = pow025(Cmu_);
const tmp<volScalarField> tk = turbulence.k();
const volScalarField& k = tk();
const scalarField& rhow = turbulence.rho().boundaryField()[patchI];
const tmp<volScalarField> tmu = turbulence.mu();
const scalarField& muw = tmu().boundaryField()[patchI];
const tmp<volScalarField> tmut = turbulence.mut();
const volScalarField& mut = tmut();
const scalarField& mutw = mut.boundaryField()[patchI];
const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
const scalarField magGradUw(mag(Uw.snGrad()));
// Set omega and G
forAll(mutw, faceI)
{
label cellI = patch.faceCells()[faceI];
scalar w = cornerWeights[faceI];
scalar omegaVis = 6.0*muw[faceI]/(rhow[faceI]*beta1_*sqr(y[faceI]));
scalar omegaLog = sqrt(k[cellI])/(Cmu25*kappa_*y[faceI]);
omega[cellI] = w*sqrt(sqr(omegaVis) + sqr(omegaLog));
G[cellI] =
w
*(mutw[faceI] + muw[faceI])
*magGradUw[faceI]
*Cmu25*sqrt(k[cellI])
/(kappa_*y[faceI]);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
@ -71,13 +264,17 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
const DimensionedField<scalar, volMesh>& iF
)
:
fixedInternalValueFvPatchField<scalar>(p, iF),
fixedValueFvPatchField<scalar>(p, iF),
Cmu_(0.09),
kappa_(0.41),
E_(9.8),
beta1_(0.075),
yPlusLam_(mutWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_))
yPlusLam_(mutWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
G_(),
omega_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
@ -91,12 +288,17 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
const fvPatchFieldMapper& mapper
)
:
fixedInternalValueFvPatchField<scalar>(ptf, p, iF, mapper),
fixedValueFvPatchField<scalar>(ptf, p, iF, mapper),
Cmu_(ptf.Cmu_),
kappa_(ptf.kappa_),
E_(ptf.E_),
beta1_(ptf.beta1_),
yPlusLam_(ptf.yPlusLam_)
yPlusLam_(ptf.yPlusLam_),
G_(),
omega_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
@ -109,12 +311,17 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
const dictionary& dict
)
:
fixedInternalValueFvPatchField<scalar>(p, iF, dict),
fixedValueFvPatchField<scalar>(p, iF, dict),
Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
E_(dict.lookupOrDefault<scalar>("E", 9.8)),
beta1_(dict.lookupOrDefault<scalar>("beta1", 0.075)),
yPlusLam_(mutWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_))
yPlusLam_(mutWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
G_(),
omega_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
@ -125,12 +332,17 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
const omegaWallFunctionFvPatchScalarField& owfpsf
)
:
fixedInternalValueFvPatchField<scalar>(owfpsf),
fixedValueFvPatchField<scalar>(owfpsf),
Cmu_(owfpsf.Cmu_),
kappa_(owfpsf.kappa_),
E_(owfpsf.E_),
beta1_(owfpsf.beta1_),
yPlusLam_(owfpsf.yPlusLam_)
yPlusLam_(owfpsf.yPlusLam_),
G_(),
omega_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
@ -142,13 +354,17 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
const DimensionedField<scalar, volMesh>& iF
)
:
fixedInternalValueFvPatchField<scalar>(owfpsf, iF),
fixedValueFvPatchField<scalar>(owfpsf, iF),
Cmu_(owfpsf.Cmu_),
kappa_(owfpsf.kappa_),
E_(owfpsf.E_),
beta1_(owfpsf.beta1_),
yPlusLam_(owfpsf.yPlusLam_)
yPlusLam_(owfpsf.yPlusLam_),
G_(),
omega_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
@ -156,6 +372,38 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
scalarField& omegaWallFunctionFvPatchScalarField::G(bool init)
{
if (patch().index() == master_)
{
if (init)
{
G_ = 0.0;
}
return G_;
}
return omegaPatch(master_).G();
}
scalarField& omegaWallFunctionFvPatchScalarField::omega(bool init)
{
if (patch().index() == master_)
{
if (init)
{
omega_ = 0.0;
}
return omega_;
}
return omegaPatch(master_).omega(init);
}
void omegaWallFunctionFvPatchScalarField::updateCoeffs()
{
if (updated())
@ -163,68 +411,168 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs()
return;
}
const label patchI = patch().index();
const turbulenceModel& turbulence =
db().lookupObject<turbulenceModel>("turbulenceModel");
const scalarField& y = turbulence.y()[patch().index()];
db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
const scalar Cmu25 = pow025(Cmu_);
setMaster();
volScalarField& G = const_cast<volScalarField&>
(
db().lookupObject<volScalarField>
(
turbulence.GName()
)
);
DimensionedField<scalar, volMesh>& omega =
const_cast<DimensionedField<scalar, volMesh>&>
(
dimensionedInternalField()
);
const tmp<volScalarField> tk = turbulence.k();
const volScalarField& k = tk();
const scalarField& rhow = turbulence.rho().boundaryField()[patchI];
const scalarField& muw = turbulence.mu().boundaryField()[patchI];
const tmp<volScalarField> tmut = turbulence.mut();
const volScalarField& mut = tmut();
const scalarField& mutw = mut.boundaryField()[patchI];
const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
const scalarField magGradUw(mag(Uw.snGrad()));
// Set omega and G
forAll(mutw, faceI)
if (patch().index() == master_)
{
label faceCellI = patch().faceCells()[faceI];
scalar omegaVis = 6.0*muw[faceI]/(rhow[faceI]*beta1_*sqr(y[faceI]));
scalar omegaLog = sqrt(k[faceCellI])/(Cmu25*kappa_*y[faceI]);
omega[faceCellI] = sqrt(sqr(omegaVis) + sqr(omegaLog));
G[faceCellI] =
(mutw[faceI] + muw[faceI])
*magGradUw[faceI]
*Cmu25*sqrt(k[faceCellI])
/(kappa_*y[faceI]);
createAveragingWeights();
calculateTurbulenceFields(turbulence, G(true), omega(true));
}
fixedInternalValueFvPatchField<scalar>::updateCoeffs();
const scalarField& G0 = this->G();
const scalarField& omega0 = this->omega();
// TODO: perform averaging for cells sharing more than one boundary face
typedef DimensionedField<scalar, volMesh> FieldType;
FieldType& G =
const_cast<FieldType&>
(
db().lookupObject<FieldType>(turbulence.GName())
);
FieldType& omega = const_cast<FieldType&>(dimensionedInternalField());
forAll(*this, faceI)
{
label cellI = patch().faceCells()[faceI];
G[cellI] = G0[cellI];
omega[cellI] = omega0[cellI];
}
fvPatchField<scalar>::updateCoeffs();
}
void omegaWallFunctionFvPatchScalarField::updateCoeffs
(
const scalarField& weights
)
{
if (updated())
{
return;
}
const turbulenceModel& turbulence =
db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
setMaster();
if (patch().index() == master_)
{
createAveragingWeights();
calculateTurbulenceFields(turbulence, G(true), omega(true));
}
const scalarField& G0 = this->G();
const scalarField& omega0 = this->omega();
typedef DimensionedField<scalar, volMesh> FieldType;
FieldType& G =
const_cast<FieldType&>
(
db().lookupObject<FieldType>(turbulence.GName())
);
FieldType& omega = const_cast<FieldType&>(dimensionedInternalField());
// only set the values if the weights are < 1 - tolerance
forAll(weights, faceI)
{
scalar w = weights[faceI];
if (w < 1.0 - 1e-6)
{
label cellI = patch().faceCells()[faceI];
G[cellI] = w*G[cellI] + (1.0 - w)*G0[cellI];
omega[cellI] = w*omega[cellI] + (1.0 - w)*omega0[cellI];
}
}
fvPatchField<scalar>::updateCoeffs();
}
void omegaWallFunctionFvPatchScalarField::manipulateMatrix
(
fvMatrix<scalar>& matrix
)
{
if (manipulatedMatrix())
{
return;
}
matrix.setValues(patch().faceCells(), patchInternalField());
fvPatchField<scalar>::manipulateMatrix(matrix);
}
void omegaWallFunctionFvPatchScalarField::manipulateMatrix
(
fvMatrix<scalar>& matrix,
const Field<scalar>& weights
)
{
if (manipulatedMatrix())
{
return;
}
// filter weights so that we only apply the constraint where the
// weight > SMALL
DynamicList<label> constraintCells(weights.size());
DynamicList<scalar> constraintomega(weights.size());
const labelUList& faceCells = patch().faceCells();
const DimensionedField<scalar, volMesh>& omega
= dimensionedInternalField();
label nConstrainedCells = 0;
forAll(weights, faceI)
{
// only set the values if the weights are < 1 - tolerance
if (weights[faceI] < (1.0 - 1e-6))
{
nConstrainedCells++;
label cellI = faceCells[faceI];
constraintCells.append(cellI);
constraintomega.append(omega[cellI]);
}
}
if (debug)
{
Pout<< "Patch: " << patch().name()
<< ": number of constrained cells = " << nConstrainedCells
<< " out of " << patch().size()
<< endl;
}
matrix.setValues
(
constraintCells,
scalarField(constraintomega.xfer())
);
fvPatchField<scalar>::manipulateMatrix(matrix);
}
void omegaWallFunctionFvPatchScalarField::write(Ostream& os) const
{
fixedInternalValueFvPatchField<scalar>::write(os);
fixedValueFvPatchField<scalar>::write(os);
writeLocalEntries(os);
writeEntry("value", os);
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -76,7 +76,7 @@ SourceFiles
#ifndef compressibleOmegaWallFunctionFvPatchScalarField_H
#define compressibleOmegaWallFunctionFvPatchScalarField_H
#include "fixedInternalValueFvPatchField.H"
#include "fixedValueFvPatchField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -85,13 +85,15 @@ namespace Foam
namespace compressible
{
class turbulenceModel;
/*---------------------------------------------------------------------------*\
Class omegaWallFunctionFvPatchScalarField Declaration
Class omegaWallFunctionFvPatchScalarField Declaration
\*---------------------------------------------------------------------------*/
class omegaWallFunctionFvPatchScalarField
:
public fixedInternalValueFvPatchField<scalar>
public fixedValueFvPatchField<scalar>
{
protected:
@ -112,6 +114,21 @@ protected:
//- Y+ at the edge of the laminar sublayer
scalar yPlusLam_;
//- Local copy of turbulence G field
scalarField G_;
//- Local copy of turbulence omega field
scalarField omega_;
//- Initialised flag
bool initialised_;
//- Master patch ID
label master_;
//- List of averaging corner weights
List<List<scalar> > cornerWeights_;
// Protected Member Functions
@ -121,6 +138,44 @@ protected:
//- Write local wall function variables
virtual void writeLocalEntries(Ostream&) const;
//- Set the master patch - master is responsible for updating all
// wall function patches
virtual void setMaster();
//- Create the averaging weights for cells which are bounded by
// multiple wall function faces
virtual void createAveragingWeights();
//- Helper function to return non-const access to an omega patch
virtual omegaWallFunctionFvPatchScalarField& omegaPatch
(
const label patchI
);
//- Main driver to calculate the turbulence fields
virtual void calculateTurbulenceFields
(
const turbulenceModel& turbulence,
scalarField& G0,
scalarField& omega0
);
//- Calculate the omega and G
virtual void calculate
(
const turbulenceModel& turbulence,
const List<scalar>& cornerWeights,
const fvPatch& patch,
scalarField& G,
scalarField& omega
);
//- Return non-const access to the master patch ID
virtual label& master()
{
return master_;
}
public:
@ -193,11 +248,33 @@ public:
// Member functions
// Access
//- Return non-const access to the master's G field
scalarField& G(bool init = false);
//- Return non-const access to the master's omega field
scalarField& omega(bool init = false);
// Evaluation functions
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Update the coefficients associated with the patch field
virtual void updateCoeffs(const scalarField& weights);
//- Manipulate matrix
virtual void manipulateMatrix(fvMatrix<scalar>& matrix);
//- Manipulate matrix with given weights
virtual void manipulateMatrix
(
fvMatrix<scalar>& matrix,
const scalarField& weights
);
// I-O

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,6 +25,7 @@ License
#include "backwardsCompatibilityWallFunctions.H"
#include "volFields.H"
#include "calculatedFvPatchField.H"
#include "nutkWallFunctionFvPatchScalarField.H"
#include "nutLowReWallFunctionFvPatchScalarField.H"

View File

@ -55,6 +55,64 @@ scalar epsilonLowReWallFunctionFvPatchScalarField::yPlusLam
}
void epsilonLowReWallFunctionFvPatchScalarField::calculate
(
const turbulenceModel& turbulence,
const List<scalar>& cornerWeights,
const fvPatch& patch,
scalarField& G,
scalarField& epsilon
)
{
const label patchI = patch.index();
const scalarField& y = turbulence.y()[patchI];
const scalar Cmu25 = pow025(Cmu_);
const scalar Cmu75 = pow(Cmu_, 0.75);
const tmp<volScalarField> tk = turbulence.k();
const volScalarField& k = tk();
const tmp<volScalarField> tnu = turbulence.nu();
const scalarField& nuw = tnu().boundaryField()[patchI];
const tmp<volScalarField> tnut = turbulence.nut();
const volScalarField& nut = tnut();
const scalarField& nutw = nut.boundaryField()[patchI];
const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
const scalarField magGradUw(mag(Uw.snGrad()));
// Set epsilon and G
forAll(nutw, faceI)
{
label cellI = patch.faceCells()[faceI];
scalar yPlus = Cmu25*sqrt(k[cellI])*y[faceI]/nuw[faceI];
scalar w = cornerWeights[faceI];
if (yPlus > yPlusLam_)
{
epsilon[cellI] = w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
}
else
{
epsilon[cellI] = w*2.0*k[cellI]*nuw[faceI]/sqr(y[faceI]);
}
G[cellI] =
w
*(nutw[faceI] + nuw[faceI])
*magGradUw[faceI]
*Cmu25*sqrt(k[cellI])
/(kappa_*y[faceI]);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
epsilonLowReWallFunctionFvPatchScalarField::
@ -119,81 +177,6 @@ epsilonLowReWallFunctionFvPatchScalarField
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void epsilonLowReWallFunctionFvPatchScalarField::updateCoeffs()
{
if (updated())
{
return;
}
const label patchI = patch().index();
const turbulenceModel& turbulence =
db().lookupObject<turbulenceModel>("turbulenceModel");
const scalarField& y = turbulence.y()[patchI];
volScalarField& G =
const_cast<volScalarField&>
(
db().lookupObject<volScalarField>
(
turbulence.GName()
)
);
DimensionedField<scalar, volMesh>& epsilon =
const_cast<DimensionedField<scalar, volMesh>&>
(
dimensionedInternalField()
);
const tmp<volScalarField> tk = turbulence.k();
const volScalarField& k = tk();
const tmp<volScalarField> tnu = turbulence.nu();
const scalarField& nuw = tnu().boundaryField()[patchI];
const tmp<volScalarField> tnut = turbulence.nut();
const volScalarField& nut = tnut();
const scalarField& nutw = nut.boundaryField()[patchI];
const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
const scalarField magGradUw(mag(Uw.snGrad()));
const scalar Cmu25 = pow025(Cmu_);
const scalar Cmu75 = pow(Cmu_, 0.75);
// Set epsilon and G
forAll(nutw, faceI)
{
label faceCellI = patch().faceCells()[faceI];
scalar yPlus = Cmu25*sqrt(k[faceCellI])*y[faceI]/nuw[faceI];
if (yPlus > yPlusLam_)
{
epsilon[faceCellI] = Cmu75*pow(k[faceCellI], 1.5)/(kappa_*y[faceI]);
}
else
{
epsilon[faceCellI] = 2.0*k[faceCellI]*nuw[faceI]/sqr(y[faceI]);
}
G[faceCellI] =
(nutw[faceI] + nuw[faceI])
*magGradUw[faceI]
*Cmu25*sqrt(k[faceCellI])
/(kappa_*y[faceI]);
}
fixedInternalValueFvPatchField<scalar>::updateCoeffs();
// TODO: perform averaging for cells sharing more than one boundary face
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makePatchTypeField

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -97,6 +97,16 @@ protected:
//- Calculate the Y+ at the edge of the laminar sublayer
scalar yPlusLam(const scalar kappa, const scalar E);
//- Calculate the epsilon and G
virtual void calculate
(
const turbulenceModel& turbulence,
const List<scalar>& cornerWeights,
const fvPatch& patch,
scalarField& G,
scalarField& epsilon
);
public:
@ -166,13 +176,9 @@ public:
);
}
// Member functions
// Evaluation functions
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Destructor
virtual ~epsilonLowReWallFunctionFvPatchScalarField()
{}
};

View File

@ -26,11 +26,10 @@ License
#include "epsilonWallFunctionFvPatchScalarField.H"
#include "incompressible/turbulenceModel/turbulenceModel.H"
#include "fvPatchFieldMapper.H"
#include "fvMatrix.H"
#include "volFields.H"
#include "addToRunTimeSelectionTable.H"
#include "wallFvPatch.H"
#include "nutkWallFunctionFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -39,7 +38,7 @@ namespace Foam
namespace incompressible
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
void epsilonWallFunctionFvPatchScalarField::checkType()
{
@ -63,118 +62,160 @@ void epsilonWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF
)
:
fixedInternalValueFvPatchField<scalar>(p, iF),
Cmu_(0.09),
kappa_(0.41),
E_(9.8)
void epsilonWallFunctionFvPatchScalarField::setMaster()
{
checkType();
}
epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
(
const epsilonWallFunctionFvPatchScalarField& ptf,
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
fixedInternalValueFvPatchField<scalar>(ptf, p, iF, mapper),
Cmu_(ptf.Cmu_),
kappa_(ptf.kappa_),
E_(ptf.E_)
{
checkType();
}
epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const dictionary& dict
)
:
fixedInternalValueFvPatchField<scalar>(p, iF, dict),
Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
E_(dict.lookupOrDefault<scalar>("E", 9.8))
{
checkType();
}
epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
(
const epsilonWallFunctionFvPatchScalarField& ewfpsf
)
:
fixedInternalValueFvPatchField<scalar>(ewfpsf),
Cmu_(ewfpsf.Cmu_),
kappa_(ewfpsf.kappa_),
E_(ewfpsf.E_)
{
checkType();
}
epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
(
const epsilonWallFunctionFvPatchScalarField& ewfpsf,
const DimensionedField<scalar, volMesh>& iF
)
:
fixedInternalValueFvPatchField<scalar>(ewfpsf, iF),
Cmu_(ewfpsf.Cmu_),
kappa_(ewfpsf.kappa_),
E_(ewfpsf.E_)
{
checkType();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
{
if (updated())
if (master_ != -1)
{
return;
}
const label patchI = patch().index();
const volScalarField& epsilon =
static_cast<const volScalarField&>(this->dimensionedInternalField());
const volScalarField::GeometricBoundaryField& bf = epsilon.boundaryField();
label master = -1;
forAll(bf, patchI)
{
if (isA<epsilonWallFunctionFvPatchScalarField>(bf[patchI]))
{
epsilonWallFunctionFvPatchScalarField& epf = epsilonPatch(patchI);
if (master == -1)
{
master = patchI;
}
epf.master() = master;
}
}
}
void epsilonWallFunctionFvPatchScalarField::createAveragingWeights()
{
if (initialised_)
{
return;
}
const volScalarField& epsilon =
static_cast<const volScalarField&>(this->dimensionedInternalField());
const volScalarField::GeometricBoundaryField& bf = epsilon.boundaryField();
const fvMesh& mesh = epsilon.mesh();
volScalarField weights
(
IOobject
(
"weights",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false // do not register
),
mesh,
dimensionedScalar("zero", dimless, 0.0)
);
DynamicList<label> epsilonPatches(bf.size());
forAll(bf, patchI)
{
if (isA<epsilonWallFunctionFvPatchScalarField>(bf[patchI]))
{
epsilonPatches.append(patchI);
const labelUList& faceCells = bf[patchI].patch().faceCells();
forAll(faceCells, i)
{
label cellI = faceCells[i];
weights[cellI]++;
}
}
}
cornerWeights_.setSize(bf.size());
forAll(epsilonPatches, i)
{
label patchI = epsilonPatches[i];
const fvPatchField& wf = weights.boundaryField()[patchI];
cornerWeights_[patchI] = 1.0/wf.patchInternalField();
}
G_.setSize(dimensionedInternalField().size(), 0.0);
epsilon_.setSize(dimensionedInternalField().size(), 0.0);
initialised_ = true;
}
epsilonWallFunctionFvPatchScalarField&
epsilonWallFunctionFvPatchScalarField::epsilonPatch(const label patchI)
{
const volScalarField& epsilon =
static_cast<const volScalarField&>(this->dimensionedInternalField());
const volScalarField::GeometricBoundaryField& bf = epsilon.boundaryField();
const epsilonWallFunctionFvPatchScalarField& epf =
refCast<const epsilonWallFunctionFvPatchScalarField>(bf[patchI]);
return const_cast<epsilonWallFunctionFvPatchScalarField&>(epf);
}
void epsilonWallFunctionFvPatchScalarField::calculateTurbulenceFields
(
const turbulenceModel& turbulence,
scalarField& G0,
scalarField& epsilon0
)
{
// accumulate all of the G and epsilon contributions
forAll(cornerWeights_, patchI)
{
if (!cornerWeights_[patchI].empty())
{
epsilonWallFunctionFvPatchScalarField& epf = epsilonPatch(patchI);
const List<scalar>& w = cornerWeights_[patchI];
epf.calculate(turbulence, w, epf.patch(), G0, epsilon0);
}
}
// apply zero-gradient condition for epsilon
forAll(cornerWeights_, patchI)
{
if (!cornerWeights_[patchI].empty())
{
epsilonWallFunctionFvPatchScalarField& epf = epsilonPatch(patchI);
epf == scalarField(epsilon0, epf.patch().faceCells());
}
}
}
void epsilonWallFunctionFvPatchScalarField::calculate
(
const turbulenceModel& turbulence,
const List<scalar>& cornerWeights,
const fvPatch& patch,
scalarField& G,
scalarField& epsilon
)
{
const label patchI = patch.index();
const turbulenceModel& turbulence =
db().lookupObject<turbulenceModel>("turbulenceModel");
const scalarField& y = turbulence.y()[patchI];
const scalar Cmu25 = pow025(Cmu_);
const scalar Cmu75 = pow(Cmu_, 0.75);
volScalarField& G =
const_cast<volScalarField&>
(
db().lookupObject<volScalarField>
(
turbulence.GName()
)
);
DimensionedField<scalar, volMesh>& epsilon =
const_cast<DimensionedField<scalar, volMesh>&>
(
dimensionedInternalField()
);
const tmp<volScalarField> tk = turbulence.k();
const volScalarField& k = tk();
@ -192,35 +233,329 @@ void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
// Set epsilon and G
forAll(nutw, faceI)
{
label faceCellI = patch().faceCells()[faceI];
label cellI = patch.faceCells()[faceI];
epsilon[faceCellI] = Cmu75*pow(k[faceCellI], 1.5)/(kappa_*y[faceI]);
scalar w = cornerWeights[faceI];
G[faceCellI] =
(nutw[faceI] + nuw[faceI])
epsilon[cellI] += w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
G[cellI] +=
w
*(nutw[faceI] + nuw[faceI])
*magGradUw[faceI]
*Cmu25*sqrt(k[faceCellI])
*Cmu25*sqrt(k[cellI])
/(kappa_*y[faceI]);
}
fixedInternalValueFvPatchField<scalar>::updateCoeffs();
// TODO: perform averaging for cells sharing more than one boundary face
}
void epsilonWallFunctionFvPatchScalarField::evaluate
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
(
const Pstream::commsTypes commsType
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchField<scalar>(p, iF),
Cmu_(0.09),
kappa_(0.41),
E_(9.8),
G_(),
epsilon_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
(
const epsilonWallFunctionFvPatchScalarField& ptf,
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
fixedValueFvPatchField<scalar>(ptf, p, iF, mapper),
Cmu_(ptf.Cmu_),
kappa_(ptf.kappa_),
E_(ptf.E_),
G_(),
epsilon_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const dictionary& dict
)
:
fixedValueFvPatchField<scalar>(p, iF, dict),
Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
E_(dict.lookupOrDefault<scalar>("E", 9.8)),
G_(),
epsilon_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
(
const epsilonWallFunctionFvPatchScalarField& ewfpsf
)
:
fixedValueFvPatchField<scalar>(ewfpsf),
Cmu_(ewfpsf.Cmu_),
kappa_(ewfpsf.kappa_),
E_(ewfpsf.E_),
G_(),
epsilon_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
(
const epsilonWallFunctionFvPatchScalarField& ewfpsf,
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchField<scalar>(ewfpsf, iF),
Cmu_(ewfpsf.Cmu_),
kappa_(ewfpsf.kappa_),
E_(ewfpsf.E_),
G_(),
epsilon_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
scalarField& epsilonWallFunctionFvPatchScalarField::G(bool init)
{
if (patch().index() == master_)
{
if (init)
{
G_ = 0.0;
}
return G_;
}
return epsilonPatch(master_).G();
}
scalarField& epsilonWallFunctionFvPatchScalarField::epsilon(bool init)
{
if (patch().index() == master_)
{
if (init)
{
epsilon_ = 0.0;
}
return epsilon_;
}
return epsilonPatch(master_).epsilon(init);
}
void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
{
if (updated())
{
return;
}
const turbulenceModel& turbulence =
db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
setMaster();
if (patch().index() == master_)
{
createAveragingWeights();
calculateTurbulenceFields(turbulence, G(true), epsilon(true));
}
const scalarField& G0 = this->G();
const scalarField& epsilon0 = this->epsilon();
typedef DimensionedField<scalar, volMesh> FieldType;
FieldType& G =
const_cast<FieldType&>
(
db().lookupObject<FieldType>(turbulence.GName())
);
FieldType& epsilon = const_cast<FieldType&>(dimensionedInternalField());
forAll(*this, faceI)
{
label cellI = patch().faceCells()[faceI];
G[cellI] = G0[cellI];
epsilon[cellI] = epsilon0[cellI];
}
fvPatchField<scalar>::updateCoeffs();
}
void epsilonWallFunctionFvPatchScalarField::updateCoeffs
(
const scalarField& weights
)
{
fixedInternalValueFvPatchField<scalar>::evaluate(commsType);
if (updated())
{
return;
}
const turbulenceModel& turbulence =
db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
setMaster();
if (patch().index() == master_)
{
createAveragingWeights();
calculateTurbulenceFields(turbulence, G(true), epsilon(true));
}
const scalarField& G0 = this->G();
const scalarField& epsilon0 = this->epsilon();
typedef DimensionedField<scalar, volMesh> FieldType;
FieldType& G =
const_cast<FieldType&>
(
db().lookupObject<FieldType>(turbulence.GName())
);
FieldType& epsilon = const_cast<FieldType&>(dimensionedInternalField());
// only set the values if the weights are < 1 - tolerance
forAll(weights, faceI)
{
scalar w = weights[faceI];
if (w < 1.0 - 1e-6)
{
label cellI = patch().faceCells()[faceI];
G[cellI] = w*G[cellI] + (1.0 - w)*G0[cellI];
epsilon[cellI] = w*epsilon[cellI] + (1.0 - w)*epsilon0[cellI];
}
}
fvPatchField<scalar>::updateCoeffs();
}
void epsilonWallFunctionFvPatchScalarField::manipulateMatrix
(
fvMatrix<scalar>& matrix
)
{
if (manipulatedMatrix())
{
return;
}
matrix.setValues(patch().faceCells(), patchInternalField());
fvPatchField<scalar>::manipulateMatrix(matrix);
}
void epsilonWallFunctionFvPatchScalarField::manipulateMatrix
(
fvMatrix<scalar>& matrix,
const Field<scalar>& weights
)
{
if (manipulatedMatrix())
{
return;
}
// filter weights so that we only apply the constraint where the
// weight > SMALL
DynamicList<label> constraintCells(weights.size());
DynamicList<scalar> constraintEpsilon(weights.size());
const labelUList& faceCells = patch().faceCells();
const DimensionedField<scalar, volMesh>& epsilon
= dimensionedInternalField();
label nConstrainedCells = 0;
forAll(weights, faceI)
{
// only set the values if the weights are < 1 - tolerance
if (weights[faceI] < (1.0 - 1e-6))
{
nConstrainedCells++;
label cellI = faceCells[faceI];
constraintCells.append(cellI);
constraintEpsilon.append(epsilon[cellI]);
}
}
if (debug)
{
Pout<< "Patch: " << patch().name()
<< ": number of constrained cells = " << nConstrainedCells
<< " out of " << patch().size()
<< endl;
}
matrix.setValues
(
constraintCells,
scalarField(constraintEpsilon.xfer())
);
fvPatchField<scalar>::manipulateMatrix(matrix);
}
void epsilonWallFunctionFvPatchScalarField::write(Ostream& os) const
{
fixedInternalValueFvPatchField<scalar>::write(os);
fixedValueFvPatchField<scalar>::write(os);
writeLocalEntries(os);
writeEntry("value", os);
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -71,7 +71,7 @@ SourceFiles
#ifndef epsilonWallFunctionFvPatchScalarField_H
#define epsilonWallFunctionFvPatchScalarField_H
#include "fixedInternalValueFvPatchField.H"
#include "fixedValueFvPatchField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -80,13 +80,15 @@ namespace Foam
namespace incompressible
{
class turbulenceModel;
/*---------------------------------------------------------------------------*\
Class epsilonWallFunctionFvPatchScalarField Declaration
\*---------------------------------------------------------------------------*/
class epsilonWallFunctionFvPatchScalarField
:
public fixedInternalValueFvPatchField<scalar>
public fixedValueFvPatchField<scalar>
{
protected:
@ -101,6 +103,21 @@ protected:
//- E coefficient
scalar E_;
//- Local copy of turbulence G field
scalarField G_;
//- Local copy of turbulence epsilon field
scalarField epsilon_;
//- Initialised flag
bool initialised_;
//- Master patch ID
label master_;
//- List of averaging corner weights
List<List<scalar> > cornerWeights_;
// Protected Member Functions
@ -110,6 +127,44 @@ protected:
//- Write local wall function variables
virtual void writeLocalEntries(Ostream&) const;
//- Set the master patch - master is responsible for updating all
// wall function patches
virtual void setMaster();
//- Create the averaging weights for cells which are bounded by
// multiple wall function faces
virtual void createAveragingWeights();
//- Helper function to return non-const access to an epsilon patch
virtual epsilonWallFunctionFvPatchScalarField& epsilonPatch
(
const label patchI
);
//- Main driver to calculate the turbulence fields
virtual void calculateTurbulenceFields
(
const turbulenceModel& turbulence,
scalarField& G0,
scalarField& epsilon0
);
//- Calculate the epsilon and G
virtual void calculate
(
const turbulenceModel& turbulence,
const List<scalar>& cornerWeights,
const fvPatch& patch,
scalarField& G,
scalarField& epsilon
);
//- Return non-const access to the master patch ID
virtual label& master()
{
return master_;
}
public:
@ -179,16 +234,39 @@ public:
);
}
//- Destructor
virtual ~epsilonWallFunctionFvPatchScalarField()
{}
// Member functions
// Access
//- Return non-const access to the master's G field
scalarField& G(bool init = false);
//- Return non-const access to the master's epsilon field
scalarField& epsilon(bool init = false);
// Evaluation functions
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Evaluate the patchField
virtual void evaluate(const Pstream::commsTypes);
//- Update the coefficients associated with the patch field
virtual void updateCoeffs(const scalarField& weights);
//- Manipulate matrix
virtual void manipulateMatrix(fvMatrix<scalar>& matrix);
//- Manipulate matrix with given weights
virtual void manipulateMatrix
(
fvMatrix<scalar>& matrix,
const scalarField& weights
);
// I-O

View File

@ -26,11 +26,11 @@ License
#include "omegaWallFunctionFvPatchScalarField.H"
#include "incompressible/turbulenceModel/turbulenceModel.H"
#include "fvPatchFieldMapper.H"
#include "fvMatrix.H"
#include "volFields.H"
#include "addToRunTimeSelectionTable.H"
#include "wallFvPatch.H"
#include "nutkWallFunctionFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -64,127 +64,159 @@ void omegaWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF
)
:
fixedInternalValueFvPatchField<scalar>(p, iF),
Cmu_(0.09),
kappa_(0.41),
E_(9.8),
beta1_(0.075),
yPlusLam_(nutkWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_))
void omegaWallFunctionFvPatchScalarField::setMaster()
{
checkType();
}
omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
(
const omegaWallFunctionFvPatchScalarField& ptf,
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
fixedInternalValueFvPatchField<scalar>(ptf, p, iF, mapper),
Cmu_(ptf.Cmu_),
kappa_(ptf.kappa_),
E_(ptf.E_),
beta1_(ptf.beta1_),
yPlusLam_(ptf.yPlusLam_)
{
checkType();
}
omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const dictionary& dict
)
:
fixedInternalValueFvPatchField<scalar>(p, iF, dict),
Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
E_(dict.lookupOrDefault<scalar>("E", 9.8)),
beta1_(dict.lookupOrDefault<scalar>("beta1", 0.075)),
yPlusLam_(nutkWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_))
{
checkType();
}
omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
(
const omegaWallFunctionFvPatchScalarField& owfpsf
)
:
fixedInternalValueFvPatchField<scalar>(owfpsf),
Cmu_(owfpsf.Cmu_),
kappa_(owfpsf.kappa_),
E_(owfpsf.E_),
beta1_(owfpsf.beta1_),
yPlusLam_(owfpsf.yPlusLam_)
{
checkType();
}
omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
(
const omegaWallFunctionFvPatchScalarField& owfpsf,
const DimensionedField<scalar, volMesh>& iF
)
:
fixedInternalValueFvPatchField<scalar>(owfpsf, iF),
Cmu_(owfpsf.Cmu_),
kappa_(owfpsf.kappa_),
E_(owfpsf.E_),
beta1_(owfpsf.beta1_),
yPlusLam_(owfpsf.yPlusLam_)
{
checkType();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void omegaWallFunctionFvPatchScalarField::updateCoeffs()
{
if (updated())
if (master_ != -1)
{
return;
}
const label patchI = patch().index();
const volScalarField& omega =
static_cast<const volScalarField&>(this->dimensionedInternalField());
const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
label master = -1;
forAll(bf, patchI)
{
if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchI]))
{
omegaWallFunctionFvPatchScalarField& epf = omegaPatch(patchI);
if (master == -1)
{
master = patchI;
}
epf.master() = master;
}
}
}
void omegaWallFunctionFvPatchScalarField::createAveragingWeights()
{
if (initialised_)
{
return;
}
const volScalarField& omega =
static_cast<const volScalarField&>(this->dimensionedInternalField());
const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
const fvMesh& mesh = omega.mesh();
volScalarField weights
(
IOobject
(
"weights",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false // do not register
),
mesh,
dimensionedScalar("zero", dimless, 0.0)
);
DynamicList<label> omegaPatches(bf.size());
forAll(bf, patchI)
{
if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchI]))
{
omegaPatches.append(patchI);
const labelUList& faceCells = bf[patchI].patch().faceCells();
forAll(faceCells, i)
{
label cellI = faceCells[i];
weights[cellI]++;
}
}
}
cornerWeights_.setSize(bf.size());
forAll(omegaPatches, i)
{
label patchI = omegaPatches[i];
const fvPatchField& wf = weights.boundaryField()[patchI];
cornerWeights_[patchI] = 1.0/wf.patchInternalField();
}
G_.setSize(dimensionedInternalField().size(), 0.0);
omega_.setSize(dimensionedInternalField().size(), 0.0);
initialised_ = true;
}
omegaWallFunctionFvPatchScalarField&
omegaWallFunctionFvPatchScalarField::omegaPatch(const label patchI)
{
const volScalarField& omega =
static_cast<const volScalarField&>(this->dimensionedInternalField());
const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
const omegaWallFunctionFvPatchScalarField& epf =
refCast<const omegaWallFunctionFvPatchScalarField>(bf[patchI]);
return const_cast<omegaWallFunctionFvPatchScalarField&>(epf);
}
void omegaWallFunctionFvPatchScalarField::calculateTurbulenceFields
(
const turbulenceModel& turbulence,
scalarField& G0,
scalarField& omega0
)
{
// accumulate all of the G and omega contributions
forAll(cornerWeights_, patchI)
{
if (!cornerWeights_[patchI].empty())
{
omegaWallFunctionFvPatchScalarField& epf = omegaPatch(patchI);
const List<scalar>& w = cornerWeights_[patchI];
epf.calculate(turbulence, w, epf.patch(), G0, omega0);
}
}
// apply zero-gradient condition for omega
forAll(cornerWeights_, patchI)
{
if (!cornerWeights_[patchI].empty())
{
omegaWallFunctionFvPatchScalarField& epf = omegaPatch(patchI);
epf == scalarField(omega0, epf.patch().faceCells());
}
}
}
void omegaWallFunctionFvPatchScalarField::calculate
(
const turbulenceModel& turbulence,
const List<scalar>& cornerWeights,
const fvPatch& patch,
scalarField& G,
scalarField& omega
)
{
const label patchI = patch.index();
const turbulenceModel& turbulence =
db().lookupObject<turbulenceModel>("turbulenceModel");
const scalarField& y = turbulence.y()[patchI];
const scalar Cmu25 = pow025(Cmu_);
volScalarField& G =
const_cast<volScalarField&>
(
db().lookupObject<volScalarField>
(
turbulence.GName()
)
);
DimensionedField<scalar, volMesh>& omega =
const_cast<DimensionedField<scalar, volMesh>&>
(
dimensionedInternalField()
);
const tmp<volScalarField> tk = turbulence.k();
const volScalarField& k = tk();
@ -202,30 +234,343 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs()
// Set omega and G
forAll(nutw, faceI)
{
label faceCellI = patch().faceCells()[faceI];
label cellI = patch.faceCells()[faceI];
scalar w = cornerWeights[faceI];
scalar omegaVis = 6.0*nuw[faceI]/(beta1_*sqr(y[faceI]));
scalar omegaLog = sqrt(k[faceCellI])/(Cmu25*kappa_*y[faceI]);
scalar omegaLog = sqrt(k[cellI])/(Cmu25*kappa_*y[faceI]);
omega[faceCellI] = sqrt(sqr(omegaVis) + sqr(omegaLog));
omega[cellI] = w*sqrt(sqr(omegaVis) + sqr(omegaLog));
G[faceCellI] =
(nutw[faceI] + nuw[faceI])
G[cellI] =
w
*(nutw[faceI] + nuw[faceI])
*magGradUw[faceI]
*Cmu25*sqrt(k[faceCellI])
*Cmu25*sqrt(k[cellI])
/(kappa_*y[faceI]);
}
}
fixedInternalValueFvPatchField<scalar>::updateCoeffs();
// TODO: perform averaging for cells sharing more than one boundary face
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchField<scalar>(p, iF),
Cmu_(0.09),
kappa_(0.41),
E_(9.8),
beta1_(0.075),
yPlusLam_(nutkWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
G_(),
omega_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
(
const omegaWallFunctionFvPatchScalarField& ptf,
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
fixedValueFvPatchField<scalar>(ptf, p, iF, mapper),
Cmu_(ptf.Cmu_),
kappa_(ptf.kappa_),
E_(ptf.E_),
beta1_(ptf.beta1_),
yPlusLam_(ptf.yPlusLam_),
G_(),
omega_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const dictionary& dict
)
:
fixedValueFvPatchField<scalar>(p, iF, dict),
Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
E_(dict.lookupOrDefault<scalar>("E", 9.8)),
beta1_(dict.lookupOrDefault<scalar>("beta1", 0.075)),
yPlusLam_(nutkWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
G_(),
omega_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
(
const omegaWallFunctionFvPatchScalarField& owfpsf
)
:
fixedValueFvPatchField<scalar>(owfpsf),
Cmu_(owfpsf.Cmu_),
kappa_(owfpsf.kappa_),
E_(owfpsf.E_),
beta1_(owfpsf.beta1_),
yPlusLam_(owfpsf.yPlusLam_),
G_(),
omega_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
(
const omegaWallFunctionFvPatchScalarField& owfpsf,
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchField<scalar>(owfpsf, iF),
Cmu_(owfpsf.Cmu_),
kappa_(owfpsf.kappa_),
E_(owfpsf.E_),
beta1_(owfpsf.beta1_),
yPlusLam_(owfpsf.yPlusLam_),
G_(),
omega_(),
initialised_(false),
master_(-1),
cornerWeights_()
{
checkType();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
scalarField& omegaWallFunctionFvPatchScalarField::G(bool init)
{
if (patch().index() == master_)
{
if (init)
{
G_ = 0.0;
}
return G_;
}
return omegaPatch(master_).G();
}
scalarField& omegaWallFunctionFvPatchScalarField::omega(bool init)
{
if (patch().index() == master_)
{
if (init)
{
omega_ = 0.0;
}
return omega_;
}
return omegaPatch(master_).omega(init);
}
void omegaWallFunctionFvPatchScalarField::updateCoeffs()
{
if (updated())
{
return;
}
const turbulenceModel& turbulence =
db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
setMaster();
if (patch().index() == master_)
{
createAveragingWeights();
calculateTurbulenceFields(turbulence, G(true), omega(true));
}
const scalarField& G0 = this->G();
const scalarField& omega0 = this->omega();
typedef DimensionedField<scalar, volMesh> FieldType;
FieldType& G =
const_cast<FieldType&>
(
db().lookupObject<FieldType>(turbulence.GName())
);
FieldType& omega = const_cast<FieldType&>(dimensionedInternalField());
forAll(*this, faceI)
{
label cellI = patch().faceCells()[faceI];
G[cellI] = G0[cellI];
omega[cellI] = omega0[cellI];
}
fvPatchField<scalar>::updateCoeffs();
}
void omegaWallFunctionFvPatchScalarField::updateCoeffs
(
const scalarField& weights
)
{
if (updated())
{
return;
}
const turbulenceModel& turbulence =
db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
setMaster();
if (patch().index() == master_)
{
createAveragingWeights();
calculateTurbulenceFields(turbulence, G(true), omega(true));
}
const scalarField& G0 = this->G();
const scalarField& omega0 = this->omega();
typedef DimensionedField<scalar, volMesh> FieldType;
FieldType& G =
const_cast<FieldType&>
(
db().lookupObject<FieldType>(turbulence.GName())
);
FieldType& omega = const_cast<FieldType&>(dimensionedInternalField());
// only set the values if the weights are < 1 - tolerance
forAll(weights, faceI)
{
scalar w = weights[faceI];
if (w < 1.0 - 1e-6)
{
label cellI = patch().faceCells()[faceI];
G[cellI] = w*G[cellI] + (1.0 - w)*G0[cellI];
omega[cellI] = w*omega[cellI] + (1.0 - w)*omega0[cellI];
}
}
fvPatchField<scalar>::updateCoeffs();
}
void omegaWallFunctionFvPatchScalarField::manipulateMatrix
(
fvMatrix<scalar>& matrix
)
{
if (manipulatedMatrix())
{
return;
}
matrix.setValues(patch().faceCells(), patchInternalField());
fvPatchField<scalar>::manipulateMatrix(matrix);
}
void omegaWallFunctionFvPatchScalarField::manipulateMatrix
(
fvMatrix<scalar>& matrix,
const Field<scalar>& weights
)
{
if (manipulatedMatrix())
{
return;
}
// filter weights so that we only apply the constraint where the
// weight > SMALL
DynamicList<label> constraintCells(weights.size());
DynamicList<scalar> constraintomega(weights.size());
const labelUList& faceCells = patch().faceCells();
const DimensionedField<scalar, volMesh>& omega
= dimensionedInternalField();
label nConstrainedCells = 0;
forAll(weights, faceI)
{
// only set the values if the weights are < 1 - tolerance
if (weights[faceI] < (1.0 - 1e-6))
{
nConstrainedCells++;
label cellI = faceCells[faceI];
constraintCells.append(cellI);
constraintomega.append(omega[cellI]);
}
}
if (debug)
{
Pout<< "Patch: " << patch().name()
<< ": number of constrained cells = " << nConstrainedCells
<< " out of " << patch().size()
<< endl;
}
matrix.setValues
(
constraintCells,
scalarField(constraintomega.xfer())
);
fvPatchField<scalar>::manipulateMatrix(matrix);
}
void omegaWallFunctionFvPatchScalarField::write(Ostream& os) const
{
fixedInternalValueFvPatchField<scalar>::write(os);
fixedValueFvPatchField<scalar>::write(os);
writeLocalEntries(os);
writeEntry("value", os);
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -76,7 +76,7 @@ SourceFiles
#ifndef omegaWallFunctionFvPatchScalarField_H
#define omegaWallFunctionFvPatchScalarField_H
#include "fixedInternalValueFvPatchField.H"
#include "fixedValueFvPatchField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -85,13 +85,15 @@ namespace Foam
namespace incompressible
{
class turbulenceModel;
/*---------------------------------------------------------------------------*\
Class omegaWallFunctionFvPatchScalarField Declaration
Class omegaWallFunctionFvPatchScalarField Declaration
\*---------------------------------------------------------------------------*/
class omegaWallFunctionFvPatchScalarField
:
public fixedInternalValueFvPatchField<scalar>
public fixedValueFvPatchField<scalar>
{
protected:
@ -112,6 +114,21 @@ protected:
//- Y+ at the edge of the laminar sublayer
scalar yPlusLam_;
//- Local copy of turbulence G field
scalarField G_;
//- Local copy of turbulence omega field
scalarField omega_;
//- Initialised flag
bool initialised_;
//- Master patch ID
label master_;
//- List of averaging corner weights
List<List<scalar> > cornerWeights_;
// Protected Member Functions
@ -121,6 +138,44 @@ protected:
//- Write local wall function variables
virtual void writeLocalEntries(Ostream&) const;
//- Set the master patch - master is responsible for updating all
// wall function patches
virtual void setMaster();
//- Create the averaging weights for cells which are bounded by
// multiple wall function faces
virtual void createAveragingWeights();
//- Helper function to return non-const access to an omega patch
virtual omegaWallFunctionFvPatchScalarField& omegaPatch
(
const label patchI
);
//- Main driver to calculate the turbulence fields
virtual void calculateTurbulenceFields
(
const turbulenceModel& turbulence,
scalarField& G0,
scalarField& omega0
);
//- Calculate the omega and G
virtual void calculate
(
const turbulenceModel& turbulence,
const List<scalar>& cornerWeights,
const fvPatch& patch,
scalarField& G,
scalarField& omega
);
//- Return non-const access to the master patch ID
virtual label& master()
{
return master_;
}
public:
@ -193,11 +248,33 @@ public:
// Member functions
// Access
//- Return non-const access to the master's G field
scalarField& G(bool init = false);
//- Return non-const access to the master's omega field
scalarField& omega(bool init = false);
// Evaluation functions
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Update the coefficients associated with the patch field
virtual void updateCoeffs(const scalarField& weights);
//- Manipulate matrix
virtual void manipulateMatrix(fvMatrix<scalar>& matrix);
//- Manipulate matrix with given weights
virtual void manipulateMatrix
(
fvMatrix<scalar>& matrix,
const scalarField& weights
);
// I-O