diff --git a/applications/utilities/mesh/generation/blockMesh/blockPoints.C b/applications/utilities/mesh/generation/blockMesh/blockPoints.C
index 0cb024215dbcf57f76a754fb2f7c9900fe802bb8..2fecb21bf54b2340e3e56afd9b5acc408d1d56a2 100644
--- a/applications/utilities/mesh/generation/blockMesh/blockPoints.C
+++ b/applications/utilities/mesh/generation/blockMesh/blockPoints.C
@@ -41,23 +41,23 @@ namespace Foam
void block::blockPoints()
{
// set local variables for mesh specification
- label ni = blockDef_.n().x();
- label nj = blockDef_.n().y();
- label nk = blockDef_.n().z();
+ const label ni = blockDef_.n().x();
+ const label nj = blockDef_.n().y();
+ const label nk = blockDef_.n().z();
- point start = blockDef_.points()[blockDef_.blockShape()[0]];
- point xEnd = blockDef_.points()[blockDef_.blockShape()[1]];
- point xyEnd = blockDef_.points()[blockDef_.blockShape()[2]];
- point yEnd = blockDef_.points()[blockDef_.blockShape()[3]];
+ const point p000 = blockDef_.points()[blockDef_.blockShape()[0]];
+ const point p100 = blockDef_.points()[blockDef_.blockShape()[1]];
+ const point p110 = blockDef_.points()[blockDef_.blockShape()[2]];
+ const point p010 = blockDef_.points()[blockDef_.blockShape()[3]];
- point zEnd = blockDef_.points()[blockDef_.blockShape()[4]];
- point xzEnd = blockDef_.points()[blockDef_.blockShape()[5]];
- point xyzEnd = blockDef_.points()[blockDef_.blockShape()[6]];
- point yzEnd = blockDef_.points()[blockDef_.blockShape()[7]];
+ const point p001 = blockDef_.points()[blockDef_.blockShape()[4]];
+ const point p101 = blockDef_.points()[blockDef_.blockShape()[5]];
+ const point p111 = blockDef_.points()[blockDef_.blockShape()[6]];
+ const point p011 = blockDef_.points()[blockDef_.blockShape()[7]];
- // set reference to the list of edge point and weighting factors
- const List<List<point> >& edgePoints = blockDef_.blockEdgePoints();
- const scalarListList& edgeWeights = blockDef_.blockEdgeWeights();
+ // list of edge point and weighting factors
+ const List<List<point> >& p = blockDef_.blockEdgePoints();
+ const scalarListList& w = blockDef_.blockEdgeWeights();
// generate vertices
@@ -69,193 +69,138 @@ void block::blockPoints()
{
label vertexNo = vtxLabel(i, j, k);
- vector edgex1 = start*(1.0 - edgeWeights[0][i])
- + xEnd*edgeWeights[0][i];
+ // points on edges
+ vector edgex1 = p000 + (p100 - p000)*w[0][i];
+ vector edgex2 = p010 + (p110 - p010)*w[1][i];
+ vector edgex3 = p011 + (p111 - p011)*w[2][i];
+ vector edgex4 = p001 + (p101 - p001)*w[3][i];
- vector edgex2 = yEnd*(1.0 - edgeWeights[1][i])
- + xyEnd*edgeWeights[1][i];
+ vector edgey1 = p000 + (p010 - p000)*w[4][j];
+ vector edgey2 = p100 + (p110 - p100)*w[5][j];
+ vector edgey3 = p101 + (p111 - p101)*w[6][j];
+ vector edgey4 = p001 + (p011 - p001)*w[7][j];
- vector edgex3 = yzEnd*(1.0 - edgeWeights[2][i])
- + xyzEnd*edgeWeights[2][i];
-
- vector edgex4 = zEnd*(1.0 - edgeWeights[3][i])
- + xzEnd*edgeWeights[3][i];
-
-
-
- vector edgey1 = start*(1.0 - edgeWeights[4][j])
- + yEnd*edgeWeights[4][j];
-
- vector edgey2 = xEnd*(1.0 - edgeWeights[5][j])
- + xyEnd*edgeWeights[5][j];
-
- vector edgey3 = xzEnd*(1.0 - edgeWeights[6][j])
- + xyzEnd*edgeWeights[6][j];
-
- vector edgey4 = zEnd*(1.0 - edgeWeights[7][j])
- + yzEnd*edgeWeights[7][j];
-
-
-
- vector edgez1 = start*(1.0 - edgeWeights[8][k])
- + zEnd*edgeWeights[8][k];
-
- vector edgez2 = xEnd*(1.0 - edgeWeights[9][k])
- + xzEnd*edgeWeights[9][k];
-
- vector edgez3 = xyEnd*(1.0 - edgeWeights[10][k])
- + xyzEnd*edgeWeights[10][k];
-
- vector edgez4 = yEnd*(1.0 - edgeWeights[11][k])
- + yzEnd*edgeWeights[11][k];
+ vector edgez1 = p000 + (p001 - p000)*w[8][k];
+ vector edgez2 = p100 + (p101 - p100)*w[9][k];
+ vector edgez3 = p110 + (p111 - p110)*w[10][k];
+ vector edgez4 = p010 + (p011 - p010)*w[11][k];
// calculate the importance factors for all edges
// x - direction
scalar impx1 =
(
- (1.0 - edgeWeights[0][i])
- *(1.0 - edgeWeights[4][j])
- *(1.0 - edgeWeights[8][k])
- + edgeWeights[0][i]
- *(1.0 - edgeWeights[5][j])
- *(1.0 - edgeWeights[9][k])
+ (1.0 - w[0][i])*(1.0 - w[4][j])*(1.0 - w[8][k])
+ + w[0][i]*(1.0 - w[5][j])*(1.0 - w[9][k])
);
scalar impx2 =
(
- (1.0 - edgeWeights[1][i])
- *edgeWeights[4][j]
- *(1.0 - edgeWeights[11][k])
- + edgeWeights[1][i]
- *edgeWeights[5][j]
- *(1.0 - edgeWeights[10][k])
+ (1.0 - w[1][i])*w[4][j]*(1.0 - w[11][k])
+ + w[1][i]*w[5][j]*(1.0 - w[10][k])
);
- scalar impx3 =
- (
- (1.0 - edgeWeights[2][i])
- *edgeWeights[7][j]
- *edgeWeights[11][k]
- + edgeWeights[2][i]
- *edgeWeights[6][j]
- *edgeWeights[10][k]
- );
+ scalar impx3 =
+ (
+ (1.0 - w[2][i])*w[7][j]*w[11][k]
+ + w[2][i]*w[6][j]*w[10][k]
+ );
- scalar impx4 =
- (
- (1.0 - edgeWeights[3][i])
- *(1.0 - edgeWeights[7][j])
- *edgeWeights[8][k]
- + edgeWeights[3][i]
- *(1.0 - edgeWeights[6][j])
- *edgeWeights[9][k]
- );
+ scalar impx4 =
+ (
+ (1.0 - w[3][i])*(1.0 - w[7][j])*w[8][k]
+ + w[3][i]*(1.0 - w[6][j])*w[9][k]
+ );
+ scalar magImpx = impx1 + impx2 + impx3 + impx4;
+ impx1 /= magImpx;
+ impx2 /= magImpx;
+ impx3 /= magImpx;
+ impx4 /= magImpx;
// y - direction
scalar impy1 =
(
- (1.0 - edgeWeights[4][j])
- *(1.0 - edgeWeights[0][i])
- *(1.0 - edgeWeights[8][k])
- + edgeWeights[4][j]
- *(1.0 - edgeWeights[1][i])
- *(1.0 - edgeWeights[11][k])
+ (1.0 - w[4][j])*(1.0 - w[0][i])*(1.0 - w[8][k])
+ + w[4][j]*(1.0 - w[1][i])*(1.0 - w[11][k])
);
scalar impy2 =
(
- (1.0 - edgeWeights[5][j])
- *edgeWeights[0][i]
- *(1.0 - edgeWeights[9][k])
- + edgeWeights[5][j]
- *edgeWeights[1][i]
- *(1.0 - edgeWeights[10][k])
+ (1.0 - w[5][j])*w[0][i]*(1.0 - w[9][k])
+ + w[5][j]*w[1][i]*(1.0 - w[10][k])
);
scalar impy3 =
(
- (1.0 - edgeWeights[6][j])
- *edgeWeights[3][i]
- *edgeWeights[9][k]
- + edgeWeights[6][j]
- *edgeWeights[2][i]
- *edgeWeights[10][k]
+ (1.0 - w[6][j])*w[3][i]*w[9][k]
+ + w[6][j]*w[2][i]*w[10][k]
);
scalar impy4 =
(
- (1.0 - edgeWeights[7][j])
- *(1.0 - edgeWeights[3][i])
- *edgeWeights[8][k]
- + edgeWeights[7][j]
- *(1.0 - edgeWeights[2][i])
- *edgeWeights[11][k]
+ (1.0 - w[7][j])*(1.0 - w[3][i])*w[8][k]
+ + w[7][j]*(1.0 - w[2][i])*w[11][k]
);
+ scalar magImpy = impy1 + impy2 + impy3 + impy4;
+ impy1 /= magImpy;
+ impy2 /= magImpy;
+ impy3 /= magImpy;
+ impy4 /= magImpy;
// z - direction
scalar impz1 =
(
- (1.0 - edgeWeights[8][k])
- *(1.0 - edgeWeights[0][i])
- *(1.0 - edgeWeights[4][j])
- + edgeWeights[8][k]
- *(1.0 - edgeWeights[3][i])
- *(1.0 - edgeWeights[7][j])
+ (1.0 - w[8][k])*(1.0 - w[0][i])*(1.0 - w[4][j])
+ + w[8][k]*(1.0 - w[3][i])*(1.0 - w[7][j])
);
scalar impz2 =
(
- (1.0 - edgeWeights[9][k])
- *edgeWeights[0][i]
- *(1.0 - edgeWeights[5][j])
- + edgeWeights[9][k]
- *edgeWeights[3][i]
- *(1.0 - edgeWeights[6][j])
+ (1.0 - w[9][k])*w[0][i]*(1.0 - w[5][j])
+ + w[9][k]*w[3][i]*(1.0 - w[6][j])
);
scalar impz3 =
(
- (1.0 - edgeWeights[10][k])
- *edgeWeights[1][i]
- *edgeWeights[5][j]
- + edgeWeights[10][k]
- *edgeWeights[2][i]
- *edgeWeights[6][j]
+ (1.0 - w[10][k])*w[1][i]*w[5][j]
+ + w[10][k]*w[2][i]*w[6][j]
);
scalar impz4 =
(
- (1.0 - edgeWeights[11][k])
- *(1.0 - edgeWeights[1][i])
- *edgeWeights[4][j]
- + edgeWeights[11][k]
- *(1.0 - edgeWeights[2][i])
- *edgeWeights[7][j]
+ (1.0 - w[11][k])*(1.0 - w[1][i])*w[4][j]
+ + w[11][k]*(1.0 - w[2][i])*w[7][j]
);
+ scalar magImpz = impz1 + impz2 + impz3 + impz4;
+ impz1 /= magImpz;
+ impz2 /= magImpz;
+ impz3 /= magImpz;
+ impz4 /= magImpz;
+
// calculate the correction vectors
- vector corx1 = impx1*(edgePoints[0][i] - edgex1);
- vector corx2 = impx2*(edgePoints[1][i] - edgex2);
- vector corx3 = impx3*(edgePoints[2][i] - edgex3);
- vector corx4 = impx4*(edgePoints[3][i] - edgex4);
+ vector corx1 = impx1*(p[0][i] - edgex1);
+ vector corx2 = impx2*(p[1][i] - edgex2);
+ vector corx3 = impx3*(p[2][i] - edgex3);
+ vector corx4 = impx4*(p[3][i] - edgex4);
- vector cory1 = impy1*(edgePoints[4][j] - edgey1);
- vector cory2 = impy2*(edgePoints[5][j] - edgey2);
- vector cory3 = impy3*(edgePoints[6][j] - edgey3);
- vector cory4 = impy4*(edgePoints[7][j] - edgey4);
+ vector cory1 = impy1*(p[4][j] - edgey1);
+ vector cory2 = impy2*(p[5][j] - edgey2);
+ vector cory3 = impy3*(p[6][j] - edgey3);
+ vector cory4 = impy4*(p[7][j] - edgey4);
- vector corz1 = impz1*(edgePoints[8][k] - edgez1);
- vector corz2 = impz2*(edgePoints[9][k] - edgez2);
- vector corz3 = impz3*(edgePoints[10][k] - edgez3);
- vector corz4 = impz4*(edgePoints[11][k] - edgez4);
+ vector corz1 = impz1*(p[8][k] - edgez1);
+ vector corz2 = impz2*(p[9][k] - edgez2);
+ vector corz3 = impz3*(p[10][k] - edgez3);
+ vector corz4 = impz4*(p[11][k] - edgez4);
// multiply by the importance factor
+
// x - direction
edgex1 *= impx1;
edgex2 *= impx2;
@@ -285,7 +230,6 @@ void block::blockPoints()
vertices_[vertexNo] += corx1 + corx2 + corx3 + corx4;
vertices_[vertexNo] += cory1 + cory2 + cory3 + cory4;
vertices_[vertexNo] += corz1 + corz2 + corz3 + corz4;
-
}
}
}
diff --git a/applications/utilities/postProcessing/patch/patchAverage/patchAverage.C b/applications/utilities/postProcessing/patch/patchAverage/patchAverage.C
index 41140a33fea2bf551f1d4de65a276e5ce11dbd13..f66a5cf2f3664e25125cd391795e12a657de7857 100644
--- a/applications/utilities/postProcessing/patch/patchAverage/patchAverage.C
+++ b/applications/utilities/postProcessing/patch/patchAverage/patchAverage.C
@@ -79,12 +79,12 @@ int main(int argc, char *argv[])
Info<< " Reading volScalarField " << fieldName << endl;
volScalarField field(fieldHeader, mesh);
- scalar area = sum(mesh.magSf().boundaryField()[patchi]);
+ scalar area = gSum(mesh.magSf().boundaryField()[patchi]);
scalar sumField = 0;
if (area > 0)
{
- sumField = sum
+ sumField = gSum
(
mesh.magSf().boundaryField()[patchi]
* field.boundaryField()[patchi]
diff --git a/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C b/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C
index 3c0e18a9cf37e841cd8ca6a79b9b4f6aeeeef11a..c4eb708ff71453935bc8813cdc6059f557fd6f58 100644
--- a/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C
+++ b/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C
@@ -75,14 +75,14 @@ int main(int argc, char *argv[])
}
// Give patch area
- Info<< " Patch area = " << sum(mesh.Sf().boundaryField()[patchi]) << endl;
+ Info<< " Patch area = " << gSum(mesh.Sf().boundaryField()[patchi]) << endl;
if (fieldHeader.headerClassName() == "volScalarField")
{
Info<< " Reading volScalarField " << fieldName << endl;
volScalarField field(fieldHeader, mesh);
- vector sumField = sum
+ vector sumField = gSum
(
mesh.Sf().boundaryField()[patchi]
* field.boundaryField()[patchi]
@@ -97,7 +97,7 @@ int main(int argc, char *argv[])
Info<< " Reading surfaceScalarField " << fieldName << endl;
surfaceScalarField field(fieldHeader, mesh);
- scalar sumField = sum(field.boundaryField()[patchi]);
+ scalar sumField = gSum(field.boundaryField()[patchi]);
Info<< " Integral of " << fieldName << " over patch "
<< patchName << '[' << patchi << ']' << " = "
diff --git a/etc/apps/cint/bashrc b/etc/apps/cint/bashrc
index 74827dfeed23e180a716d4bda555d112a0cbd4bf..11909e0487036cc94cf8dc613568c1b3bab0c01c 100644
--- a/etc/apps/cint/bashrc
+++ b/etc/apps/cint/bashrc
@@ -31,7 +31,7 @@
#
#------------------------------------------------------------------------------
-export CINTSYSDIR=~/pub/CINT/cint
+export CINTSYSDIR=~/pub/CINT/cint7
export PATH=$PATH:$CINTSYSDIR
export MANPATH=$MANPATH:$CINTSYSDIR/doc
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:$CINTSYSDIR
diff --git a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
index 00ce11355a3ceebabe00dffd59f1605453f76ee2..321eeba6a87be806089b4cafd53ea2195286ca56 100644
--- a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
+++ b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
@@ -33,15 +33,68 @@ void Foam::setRefCell
const volScalarField& field,
const dictionary& dict,
label& refCelli,
- scalar& refValue
+ scalar& refValue,
+ bool forceReference
)
{
- if (field.needReference())
+ if (field.needReference() || forceReference)
{
word refCellName = field.name() + "RefCell";
+ word refPointName = field.name() + "RefPoint";
+
word refValueName = field.name() + "RefValue";
- refCelli = readLabel(dict.lookup(refCellName));
+ if (dict.found(refCellName))
+ {
+ if (Pstream::master())
+ {
+ refCelli = readLabel(dict.lookup(refCellName));
+ }
+ else
+ {
+ refCelli = -1;
+ }
+ }
+ else if (dict.found(refPointName))
+ {
+ point refPointi(dict.lookup(refPointName));
+ refCelli = field.mesh().findCell(refPointi);
+ label hasRef = (refCelli >= 0 ? 1 : 0);
+ label sumHasRef = returnReduce<label>(hasRef, sumOp<label>());
+ if (sumHasRef != 1)
+ {
+ FatalErrorIn
+ (
+ "void Foam::setRefCell"
+ "("
+ " const volScalarField&,"
+ " const dictionary&,"
+ " label& scalar&,"
+ " bool"
+ ")"
+ )
+ << "Unable to set reference cell for field " << field.name()
+ << nl << " Reference point " << refPointName
+ << " found on multiple domains" << nl << abort(FatalError);
+ }
+ }
+ else
+ {
+ FatalErrorIn
+ (
+ "void Foam::setRefCell"
+ "("
+ " const volScalarField&,"
+ " const dictionary&,"
+ " label& scalar&,"
+ " bool"
+ ")"
+ )
+ << "Unable to set reference cell for field" << field.name() << nl
+ << " Please supply either " << refCellName
+ << " or " << refPointName << nl << abort(FatalError);
+ }
+
refValue = readScalar(dict.lookup(refValueName));
}
}
diff --git a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H
index 724f44ed175186e6c765dc4b3fd7472cace9b41c..db9858ce27aed5b35cd2518281721d71c6c524ad 100644
--- a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H
+++ b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H
@@ -52,7 +52,8 @@ void setRefCell
const volScalarField& field,
const dictionary& dict,
label& refCelli,
- scalar& refValue
+ scalar& refValue,
+ bool forceReference = false
);
}
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C
index 279e72a0c55620bc270bd21ba3d65c022529434c..ab3af1717e006ce47754d5affaa9b7d93de5bc3a 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C
@@ -234,6 +234,16 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
(
mpp.samplePatch()
);
+ if (patchID < 0)
+ {
+ FatalErrorIn
+ (
+ "void directMappedFixedValueFvPatchField<Type>::"
+ "updateCoeffs()"
+ )<< "Unable to find sample patch " << mpp.samplePatch()
+ << " for patch " << this->patch().name() << nl
+ << abort(FatalError);
+ }
typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
const word& fieldName = this->dimensionedInternalField().name();
const fieldType& sendField =
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index 9dab83d168a4e57c7041331ca349b95dd1fdcf44..b8698eb17738d2e5e39adfd6ac0eaa3568356a15 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -483,7 +483,7 @@ void Foam::fvMatrix<Type>::setReference
{
if (psi_.needReference() || forceReference)
{
- if (Pstream::master())
+ if (cell >= 0)
{
source()[cell] += diag()[cell]*value;
diag()[cell] += diag()[cell];
diff --git a/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H b/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
index 754996f27c6c6cc8797596629309270fd2e0fbc3..7881a307987c9f98d2e8df05ea5ae76a61bc3107 100644
--- a/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
@@ -164,7 +164,7 @@ inline Foam::scalar Foam::janafThermo<equationOfState>::hc() const
(
(((a[4]/5.0*Tstd + a[3]/4.0)*Tstd + a[2]/3.0)*Tstd + a[1]/2.0)*Tstd
+ a[0]
- )*Tstd
+ )*Tstd + a[5]
);
}
diff --git a/src/turbulenceModels/LES/incompressible/Make/files b/src/turbulenceModels/LES/incompressible/Make/files
index b6c245c866bb81e3c6767c2f6fbde1823bd0c430..feda56ca15f600f00a5efd2f24da1f3ab9c2cdc0 100644
--- a/src/turbulenceModels/LES/incompressible/Make/files
+++ b/src/turbulenceModels/LES/incompressible/Make/files
@@ -27,6 +27,7 @@ dynMixedSmagorinsky/dynMixedSmagorinsky.C
/*Smagorinsky2/Smagorinsky2.C*/
+kOmegaSSTSAS/kOmegaSSTSAS.C
/* Wall functions */
wallFunctions=derivedFvPatchFields/wallFunctions
diff --git a/src/turbulenceModels/LES/incompressible/kOmegaSSTSAS/kOmegaSSTSAS.C b/src/turbulenceModels/LES/incompressible/kOmegaSSTSAS/kOmegaSSTSAS.C
new file mode 100644
index 0000000000000000000000000000000000000000..a06f4298c64bdb8f8f202e781682ead49d8b8773
--- /dev/null
+++ b/src/turbulenceModels/LES/incompressible/kOmegaSSTSAS/kOmegaSSTSAS.C
@@ -0,0 +1,447 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software; you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by the
+ Free Software Foundation; either version 2 of the License, or (at your
+ option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM; if not, write to the Free Software Foundation,
+ Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "kOmegaSSTSAS.H"
+#include "addToRunTimeSelectionTable.H"
+#include "wallDist.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace LESModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(kOmegaSSTSAS, 0);
+addToRunTimeSelectionTable(LESModel, kOmegaSSTSAS, dictionary);
+
+// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+
+tmp<volScalarField> kOmegaSSTSAS::F1(const volScalarField& CDkOmega) const
+{
+ volScalarField CDkOmegaPlus = max
+ (
+ CDkOmega,
+ dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
+ );
+
+ volScalarField arg1 = min
+ (
+ min
+ (
+ max
+ (
+ (scalar(1)/betaStar_)*sqrt(k_)/(omega_*y_),
+ scalar(500)*nu()/(sqr(y_)*omega_)
+ ),
+ (4*alphaOmega2_)*k_/(CDkOmegaPlus*sqr(y_))
+ ),
+ scalar(10)
+ );
+
+ return tanh(pow4(arg1));
+}
+
+
+tmp<volScalarField> kOmegaSSTSAS::F2() const
+{
+ volScalarField arg2 = min
+ (
+ max
+ (
+ (scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_),
+ scalar(500)*nu()/(sqr(y_)*omega_)
+ ),
+ scalar(100)
+ );
+
+ return tanh(sqr(arg2));
+}
+
+
+tmp<volScalarField> kOmegaSSTSAS::Lvk2
+(
+ const volScalarField& S2
+) const
+{
+ return kappa_*sqrt(S2)
+ /(
+ mag(fvc::laplacian(U()))
+ + dimensionedScalar
+ (
+ "ROOTVSMALL",
+ dimensionSet(0, -1 , -1, 0, 0, 0, 0),
+ ROOTVSMALL
+ )
+ );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+kOmegaSSTSAS::kOmegaSSTSAS
+(
+ const volVectorField& U,
+ const surfaceScalarField& phi,
+ transportModel& transport,
+ const word& modelName
+)
+:
+ LESModel(modelName, U, phi, transport),
+
+ alphaK1_
+ (
+ dimensioned<scalar>::lookupOrAddToDict
+ (
+ "alphaK1",
+ coeffDict(),
+ 0.85034
+ )
+ ),
+ alphaK2_
+ (
+ dimensioned<scalar>::lookupOrAddToDict
+ (
+ "alphaK2",
+ coeffDict(),
+ 1.0
+ )
+ ),
+ alphaOmega1_
+ (
+ dimensioned<scalar>::lookupOrAddToDict
+ (
+ "alphaOmega1",
+ coeffDict(),
+ 0.5
+ )
+ ),
+ alphaOmega2_
+ (
+ dimensioned<scalar>::lookupOrAddToDict
+ (
+ "alphaOmega2",
+ coeffDict(),
+ 0.85616
+ )
+ ),
+ gamma1_
+ (
+ dimensioned<scalar>::lookupOrAddToDict
+ (
+ "gamma1",
+ coeffDict(),
+ 0.5532
+ )
+ ),
+ gamma2_
+ (
+ dimensioned<scalar>::lookupOrAddToDict
+ (
+ "gamma2",
+ coeffDict(),
+ 0.4403
+ )
+ ),
+ beta1_
+ (
+ dimensioned<scalar>::lookupOrAddToDict
+ (
+ "beta1",
+ coeffDict(),
+ 0.075
+ )
+ ),
+ beta2_
+ (
+ dimensioned<scalar>::lookupOrAddToDict
+ (
+ "beta2",
+ coeffDict(),
+ 0.0828
+ )
+ ),
+ betaStar_
+ (
+ dimensioned<scalar>::lookupOrAddToDict
+ (
+ "betaStar",
+ coeffDict(),
+ 0.09
+ )
+ ),
+ a1_
+ (
+ dimensioned<scalar>::lookupOrAddToDict
+ (
+ "a1",
+ coeffDict(),
+ 0.31
+ )
+ ),
+ c1_
+ (
+ dimensioned<scalar>::lookupOrAddToDict
+ (
+ "c1",
+ coeffDict(),
+ 10.0
+ )
+ ),
+ alphaPhi_
+ (
+ dimensioned<scalar>::lookupOrAddToDict
+ (
+ "alphaPhi",
+ coeffDict(),
+ 0.666667
+ )
+ ),
+ zetaTilda2_
+ (
+ dimensioned<scalar>::lookupOrAddToDict
+ (
+ "zetaTilda2",
+ coeffDict(),
+ 1.755
+ )
+ ),
+ FSAS_
+ (
+ dimensioned<scalar>::lookupOrAddToDict
+ (
+ "FSAS",
+ coeffDict(),
+ 1.25
+ )
+ ),
+
+ omega0_("omega0", dimless/dimTime, SMALL),
+ omegaSmall_("omegaSmall", dimless/dimTime, SMALL),
+ y_(mesh_),
+ Cmu_
+ (
+ dimensioned<scalar>::lookupOrAddToDict
+ (
+ "Cmu",
+ coeffDict(),
+ 0.09
+ )
+ ),
+ kappa_
+ (
+ dimensioned<scalar>::lookupOrAddToDict
+ (
+ "kappa",
+ *this,
+ 0.4187
+ )
+ ),
+
+ k_
+ (
+ IOobject
+ (
+ "k",
+ runTime_.timeName(),
+ mesh_,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh_
+ ),
+
+ omega_
+ (
+ IOobject
+ (
+ "omega",
+ runTime_.timeName(),
+ mesh_,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh_
+ ),
+
+ nuSgs_
+ (
+ IOobject
+ (
+ "nuSgs",
+ runTime_.timeName(),
+ mesh_,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh_
+ )
+{
+ printCoeffs();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void kOmegaSSTSAS::correct(const tmp<volTensorField>& gradU)
+{
+ LESModel::correct(gradU);
+
+ if (mesh_.changing())
+ {
+ y_.correct();
+ }
+
+ volScalarField S2 = magSqr(symm(gradU()));
+ gradU.clear();
+
+ volVectorField gradK = fvc::grad(k_);
+ volVectorField gradOmega = fvc::grad(omega_);
+ volScalarField L = sqrt(k_)/(pow(Cmu_, 0.25)*(omega_ + omegaSmall_));
+ volScalarField CDkOmega = (2.0*alphaOmega2_)*(gradK & gradOmega)/(omega_ + omegaSmall_);
+ volScalarField F1 = this->F1(CDkOmega);
+ volScalarField G = nuSgs_*2.0*S2;
+
+ // Turbulent kinetic energy equation
+ solve
+ (
+ fvm::ddt(k_)
+ + fvm::div(phi(), k_)
+ - fvm::Sp(fvc::div(phi()), k_)
+ - fvm::laplacian(DkEff(F1), k_)
+ ==
+ min(G, c1_*betaStar_*k_*omega_)
+ - fvm::Sp(betaStar_*omega_, k_)
+ );
+
+
+ bound(k_, k0());
+
+ volScalarField grad_omega_k = max
+ (
+ magSqr(gradOmega)/
+ sqr(omega_ + omegaSmall_),
+ magSqr(gradK)/
+ sqr(k_ + k0())
+ );
+
+ // Turbulent frequency equation
+ solve
+ (
+ fvm::ddt(omega_)
+ + fvm::div(phi(), omega_)
+ - fvm::Sp(fvc::div(phi()), omega_)
+ - fvm::laplacian(DomegaEff(F1), omega_)
+ ==
+ gamma(F1)*2.0*S2
+ - fvm::Sp(beta(F1)*omega_, omega_)
+ - fvm::SuSp // cross diffusion term
+ (
+ (F1 - scalar(1))*CDkOmega/omega_,
+ omega_
+ )
+ + FSAS_
+ *max
+ (
+ dimensionedScalar("zero",dimensionSet(0, 0 , -2, 0, 0),0. ),
+ zetaTilda2_*kappa_*S2*(L/Lvk2(S2))- 2.0/alphaPhi_*k_*grad_omega_k
+ )
+
+ );
+ bound(omega_, omega0_);
+
+ // Re-calculate viscosity
+ nuSgs_ == a1_*k_/max(a1_*omega_, F2()*sqrt(S2));
+ nuSgs_.correctBoundaryConditions();
+}
+
+
+tmp<volScalarField> kOmegaSSTSAS::epsilon() const
+{
+ return 2.0*nuEff()*magSqr(symm(fvc::grad(U())));
+}
+
+
+tmp<volSymmTensorField> kOmegaSSTSAS::B() const
+{
+ return ((2.0/3.0)*I)*k() - nuSgs()*twoSymm(fvc::grad(U()));
+}
+
+
+tmp<volSymmTensorField> kOmegaSSTSAS::devBeff() const
+{
+ return -nuEff()*dev(twoSymm(fvc::grad(U())));
+}
+
+
+tmp<fvVectorMatrix> kOmegaSSTSAS::divDevBeff(volVectorField& U) const
+{
+ return
+ (
+ - fvm::laplacian(nuEff(), U) - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
+ );
+}
+
+
+bool kOmegaSSTSAS::read()
+{
+ if (LESModel::read())
+ {
+ alphaK1_.readIfPresent(coeffDict());
+ alphaK2_.readIfPresent(coeffDict());
+ alphaOmega1_.readIfPresent(coeffDict());
+ alphaOmega2_.readIfPresent(coeffDict());
+ gamma1_.readIfPresent(coeffDict());
+ gamma2_.readIfPresent(coeffDict());
+ beta1_.readIfPresent(coeffDict());
+ beta2_.readIfPresent(coeffDict());
+ betaStar_.readIfPresent(coeffDict());
+ a1_.readIfPresent(coeffDict());
+ c1_.readIfPresent(coeffDict());
+ alphaPhi_.readIfPresent(coeffDict());
+ zetaTilda2_.readIfPresent(coeffDict());
+ FSAS_.readIfPresent(coeffDict());
+
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/turbulenceModels/LES/incompressible/kOmegaSSTSAS/kOmegaSSTSAS.H b/src/turbulenceModels/LES/incompressible/kOmegaSSTSAS/kOmegaSSTSAS.H
new file mode 100644
index 0000000000000000000000000000000000000000..8048be0adf729fce3a06c0166537ed279d6dee8e
--- /dev/null
+++ b/src/turbulenceModels/LES/incompressible/kOmegaSSTSAS/kOmegaSSTSAS.H
@@ -0,0 +1,258 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software; you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by the
+ Free Software Foundation; either version 2 of the License, or (at your
+ option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM; if not, write to the Free Software Foundation,
+ Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+ Foam::LESmodels::kOmegaSSTSAS
+
+Description
+ kOmegaSSTSAS LES turbulence model for incompressible flows
+
+ Note: does not have an explicit dependency on spatial discretisation
+ i.e. LESdelta not used.
+
+SourceFiles
+ kOmegaSSTSAS.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kOmegaSSTSAS_H
+#define kOmegaSSTSAS_H
+
+#include "LESModel.H"
+#include "volFields.H"
+#include "wallDist.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace LESModels
+{
+
+/*---------------------------------------------------------------------------*\
+ Class kOmegaSSTSAS Declaration
+\*---------------------------------------------------------------------------*/
+
+class kOmegaSSTSAS
+:
+ public LESModel
+{
+ // Private member functions
+
+ // Disallow default bitwise copy construct and assignment
+ kOmegaSSTSAS(const kOmegaSSTSAS&);
+ kOmegaSSTSAS& operator=(const kOmegaSSTSAS&);
+
+
+protected:
+
+ // Protected data
+
+ // Model constants
+
+ dimensionedScalar alphaK1_;
+ dimensionedScalar alphaK2_;
+
+ dimensionedScalar alphaOmega1_;
+ dimensionedScalar alphaOmega2_;
+
+ dimensionedScalar gamma1_;
+ dimensionedScalar gamma2_;
+
+ dimensionedScalar beta1_;
+ dimensionedScalar beta2_;
+
+ dimensionedScalar betaStar_;
+
+ dimensionedScalar a1_;
+ dimensionedScalar c1_;
+
+ dimensionedScalar alphaPhi_;
+ dimensionedScalar zetaTilda2_;
+ dimensionedScalar FSAS_;
+
+ dimensionedScalar omega0_;
+ dimensionedScalar omegaSmall_;
+
+ wallDist y_;
+ dimensionedScalar Cmu_;
+ dimensionedScalar kappa_;
+
+
+ // Fields
+
+ volScalarField k_;
+ volScalarField omega_;
+ volScalarField nuSgs_;
+
+
+ // Protected member functions
+
+ tmp<volScalarField> Lvk2
+ (
+ const volScalarField& S2
+ ) const;
+
+ tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
+ tmp<volScalarField> F2() const;
+
+ tmp<volScalarField> blend
+ (
+ const volScalarField& F1,
+ const dimensionedScalar& psi1,
+ const dimensionedScalar& psi2
+ ) const
+ {
+ return F1*(psi1 - psi2) + psi2;
+ }
+
+ tmp<volScalarField> alphaK
+ (
+ const volScalarField& F1
+ ) const
+ {
+ return blend(F1, alphaK1_, alphaK2_);
+ }
+
+ tmp<volScalarField> alphaOmega
+ (
+ const volScalarField& F1
+ ) const
+ {
+ return blend(F1, alphaOmega1_, alphaOmega2_);
+ }
+
+ tmp<volScalarField> beta
+ (
+ const volScalarField& F1
+ ) const
+ {
+ return blend(F1, beta1_, beta2_);
+ }
+
+ tmp<volScalarField> gamma
+ (
+ const volScalarField& F1
+ ) const
+ {
+ return blend(F1, gamma1_, gamma2_);
+ }
+
+
+public:
+
+ //- Runtime type information
+ TypeName("kOmegaSSTSAS");
+
+
+ // Constructors
+
+ //- Constructor from components
+ kOmegaSSTSAS
+ (
+ const volVectorField& U,
+ const surfaceScalarField& phi,
+ transportModel& transport,
+ const word& modelName = typeName
+ );
+
+
+ //- Destructor
+ virtual ~kOmegaSSTSAS()
+ {}
+
+
+ // Member Functions
+
+ //- Return SGS kinetic energy
+ tmp<volScalarField> k() const
+ {
+ return k_;
+ }
+
+ //- Return omega
+ tmp<volScalarField> omega() const
+ {
+ return omega_;
+ }
+
+ //- Return the effective diffusivity for k
+ tmp<volScalarField> DkEff(const volScalarField& F1) const
+ {
+ return tmp<volScalarField>
+ (
+ new volScalarField("DkEff", alphaK(F1)*nuSgs_ + nu())
+ );
+ }
+
+ //- Return the effective diffusivity for omega
+ tmp<volScalarField> DomegaEff(const volScalarField& F1) const
+ {
+ return tmp<volScalarField>
+ (
+ new volScalarField("DomegaEff", alphaOmega(F1)*nuSgs_ + nu())
+ );
+ }
+
+ //- Return sub-grid disipation rate
+ tmp<volScalarField> epsilon() const;
+
+ //- Return SGS viscosity
+ tmp<volScalarField> nuSgs() const
+ {
+ return nuSgs_;
+ }
+
+ //- Return the sub-grid stress tensor.
+ tmp<volSymmTensorField> B() const;
+
+ //- Return the effective sub-grid turbulence stress tensor
+ // including the laminar stress
+ tmp<volSymmTensorField> devBeff() const;
+
+ //- Return the deviatoric part of the divergence of Beff
+ // i.e. the additional term in the filtered NSE.
+ tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+
+ //- Solve the turbulence equations (k-w) and correct the turbulence viscosity
+ virtual void correct(const tmp<volTensorField>& gradU);
+
+ //- Read turbulenceProperties dictionary
+ virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //