diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
index a78e0bdb84ed3c9d046dec62bb596c24cd4bc144..554c18542d9fd800d7d52ac0fd4892d4e8ec2e91 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
@@ -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,6 +64,196 @@ 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 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 omega and G
+ forAll(nutw, faceI)
+ {
+ label cellI = patch.faceCells()[faceI];
+
+ scalar w = cornerWeights[faceI];
+
+ scalar omegaVis = 6.0*nuw[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
+ *(nutw[faceI] + nuw[faceI])
+ *magGradUw[faceI]
+ *Cmu25*sqrt(k[cellI])
+ /(kappa_*y[faceI]);
+ }
+}
+
+
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
@@ -72,12 +262,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_(nutkWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_))
+ yPlusLam_(nutkWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
+ G_(),
+ omega_(),
+ initialised_(false),
+ master_(-1),
+ cornerWeights_()
{
checkType();
}
@@ -91,12 +286,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 +309,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_(nutkWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_))
+ yPlusLam_(nutkWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
+ G_(),
+ omega_(),
+ initialised_(false),
+ master_(-1),
+ cornerWeights_()
{
checkType();
}
@@ -125,12 +330,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,12 +352,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();
}
@@ -155,6 +370,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())
@@ -162,70 +409,168 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs()
return;
}
- const label patchI = patch().index();
-
const turbulenceModel& turbulence =
- db().lookupObject<turbulenceModel>("turbulenceModel");
- const scalarField& y = turbulence.y()[patchI];
+ db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
- const scalar Cmu25 = pow025(Cmu_);
+ 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;
- volScalarField& G =
- const_cast<volScalarField&>
+ FieldType& G =
+ const_cast<FieldType&>
(
- db().lookupObject<volScalarField>
- (
- turbulence.GName()
- )
+ db().lookupObject<FieldType>(turbulence.GName())
);
- DimensionedField<scalar, volMesh>& omega =
- const_cast<DimensionedField<scalar, volMesh>&>
+ 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&>
(
- dimensionedInternalField()
+ db().lookupObject<FieldType>(turbulence.GName())
);
- const tmp<volScalarField> tk = turbulence.k();
- const volScalarField& k = tk();
+ FieldType& omega = const_cast<FieldType&>(dimensionedInternalField());
- const tmp<volScalarField> tnu = turbulence.nu();
- const scalarField& nuw = tnu().boundaryField()[patchI];
+ // only set the values if the weights are < 1 - tolerance
+ forAll(weights, faceI)
+ {
+ scalar w = weights[faceI];
- const tmp<volScalarField> tnut = turbulence.nut();
- const volScalarField& nut = tnut();
- const scalarField& nutw = nut.boundaryField()[patchI];
+ if (w < 1.0 - 1e-6)
+ {
+ label cellI = patch().faceCells()[faceI];
- const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
+ G[cellI] = w*G[cellI] + (1.0 - w)*G0[cellI];
+ omega[cellI] = w*omega[cellI] + (1.0 - w)*omega0[cellI];
+ }
+ }
- const scalarField magGradUw(mag(Uw.snGrad()));
+ fvPatchField<scalar>::updateCoeffs();
+}
- // Set omega and G
- forAll(nutw, faceI)
+
+void omegaWallFunctionFvPatchScalarField::manipulateMatrix
+(
+ fvMatrix<scalar>& matrix
+)
+{
+ if (manipulatedMatrix())
{
- label faceCellI = patch().faceCells()[faceI];
+ return;
+ }
- scalar omegaVis = 6.0*nuw[faceI]/(beta1_*sqr(y[faceI]));
+ matrix.setValues(patch().faceCells(), patchInternalField());
- scalar omegaLog = sqrt(k[faceCellI])/(Cmu25*kappa_*y[faceI]);
+ fvPatchField<scalar>::manipulateMatrix(matrix);
+}
- omega[faceCellI] = sqrt(sqr(omegaVis) + sqr(omegaLog));
- G[faceCellI] =
- (nutw[faceI] + nuw[faceI])
- *magGradUw[faceI]
- *Cmu25*sqrt(k[faceCellI])
- /(kappa_*y[faceI]);
+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;
}
- fixedInternalValueFvPatchField<scalar>::updateCoeffs();
+ matrix.setValues
+ (
+ constraintCells,
+ scalarField(constraintomega.xfer())
+ );
- // TODO: perform averaging for cells sharing more than one boundary face
+ fvPatchField<scalar>::manipulateMatrix(matrix);
}
void omegaWallFunctionFvPatchScalarField::write(Ostream& os) const
{
- fixedInternalValueFvPatchField<scalar>::write(os);
+ fixedValueFvPatchField<scalar>::write(os);
writeLocalEntries(os);
writeEntry("value", os);
}
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H
index 9ee99bd6c6927e76a57325b96fbeb8b1be76d1a1..957d0f945eed05607aa6ba122c42714395d24ebb 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H
@@ -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