diff --git a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C
index ec0340f23ecc97b4c52e6a381d702c29a3d0e3ed..1e6809288d83e4730978c3f1d437ee37ac168905 100644
--- a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C
+++ b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C
@@ -47,26 +47,30 @@ addToRunTimeSelectionTable(RASModel, ShihQuadraticKE, dictionary);
void ShihQuadraticKE::correctNut()
{
- nut_ = Cmu_*sqr(k_)/epsilon_;
- #include "wallNonlinearViscosityI.H"
+ correctNonlinearStress(fvc::grad(U_));
}
void ShihQuadraticKE::correctNonlinearStress(const volTensorField& gradU)
{
- nonlinearStress_ = symm
- (
- pow(k_, 3.0)/sqr(epsilon_)
+ volSymmTensorField S(symm(gradU));
+ volTensorField W(skew(gradU));
+
+ volScalarField sBar((k_/epsilon_)*sqrt(2.0)*mag(symm(gradU)));
+ volScalarField wBar((k_/epsilon_)*sqrt(2.0)*mag(skew(gradU)));
+
+ volScalarField Cmu(2.0/(3.0*(Cmu1_ + sBar + Cmu2_*wBar)));
+
+ nut_ = Cmu*sqr(k_)/epsilon_;
+ nut_.correctBoundaryConditions();
+
+ nonlinearStress_ =
+ pow3(k_)/((A1_ + pow3(sBar))*sqr(epsilon_))
*(
- Ctau1_/fEta_
- *(
- (gradU & gradU)
- + (gradU & gradU)().T()
- )
- + Ctau2_/fEta_*(gradU & T(gradU))
- + Ctau3_/fEta_*(T(gradU) & gradU)
- )
- );
+ beta1_*dev(symm(S&S))
+ + beta2_*symm((W&S) - (S&W))
+ + beta3_*dev(symm(W&W))
+ );
}
@@ -132,7 +136,7 @@ ShihQuadraticKE::ShihQuadraticKE
1.3
)
),
- A1_
+ Cmu1_
(
dimensioned<scalar>::lookupOrAddToDict
(
@@ -141,68 +145,49 @@ ShihQuadraticKE::ShihQuadraticKE
1.25
)
),
- A2_
- (
- dimensioned<scalar>::lookupOrAddToDict
- (
- "A2",
- coeffDict_,
- 1000.0
- )
- ),
- Ctau1_
- (
- dimensioned<scalar>::lookupOrAddToDict
- (
- "Ctau1",
- coeffDict_,
- -4.0
- )
- ),
- Ctau2_
+ Cmu2_
(
dimensioned<scalar>::lookupOrAddToDict
(
- "Ctau2",
+ "alphaKsi",
coeffDict_,
- 13.0
+ 0.9
)
),
- Ctau3_
+ A1_
(
dimensioned<scalar>::lookupOrAddToDict
(
- "Ctau3",
+ "A2",
coeffDict_,
- -2.0
+ 1000.0
)
),
- alphaKsi_
+ beta1_
(
dimensioned<scalar>::lookupOrAddToDict
(
- "alphaKsi",
+ "beta1",
coeffDict_,
- 0.9
+ 3.0
)
),
-
- kappa_
+ beta2_
(
dimensioned<scalar>::lookupOrAddToDict
(
- "kappa_",
+ "beta2",
coeffDict_,
- 0.41
+ 15.0
)
),
- E_
+ beta3_
(
dimensioned<scalar>::lookupOrAddToDict
(
- "E",
+ "beta3",
coeffDict_,
- 9.8
+ -19.0
)
),
@@ -230,28 +215,14 @@ ShihQuadraticKE::ShihQuadraticKE
IOobject::AUTO_WRITE
),
mesh_
- ),
-
- eta_
- (
- k_/bound(epsilon_, epsilonMin_)
- *sqrt(2.0*magSqr(0.5*(fvc::grad(U) + T(fvc::grad(U)))))
- ),
- ksi_
- (
- k_/epsilon_
- *sqrt(2.0*magSqr(0.5*(fvc::grad(U) - T(fvc::grad(U)))))
- ),
- Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
- fEta_(A2_ + pow(eta_, 3.0))
+ )
{
bound(k_, kMin_);
- // already bounded: bound(epsilon_, epsilonMin_);
+ bound(epsilon_, epsilonMin_);
if (type == typeName)
{
correctNut();
- correctNonlinearStress(fvc::grad(U));
printCoeffs(type);
}
}
@@ -267,15 +238,12 @@ bool ShihQuadraticKE::read()
C2_.readIfPresent(coeffDict());
sigmak_.readIfPresent(coeffDict());
sigmaEps_.readIfPresent(coeffDict());
+ Cmu1_.readIfPresent(coeffDict());
+ Cmu2_.readIfPresent(coeffDict());
A1_.readIfPresent(coeffDict());
- A2_.readIfPresent(coeffDict());
- Ctau1_.readIfPresent(coeffDict());
- Ctau2_.readIfPresent(coeffDict());
- Ctau3_.readIfPresent(coeffDict());
- alphaKsi_.readIfPresent(coeffDict());
-
- kappa_.readIfPresent(coeffDict());
- E_.readIfPresent(coeffDict());
+ beta1_.readIfPresent(coeffDict());
+ beta2_.readIfPresent(coeffDict());
+ beta3_.readIfPresent(coeffDict());
return true;
}
@@ -298,16 +266,14 @@ void ShihQuadraticKE::correct()
tmp<volTensorField> tgradU = fvc::grad(U_);
const volTensorField& gradU = tgradU();
- // generation term
- tmp<volScalarField> S2 = symm(gradU) && gradU;
-
volScalarField G
(
GName(),
- Cmu_*sqr(k_)/epsilon_*S2 - (nonlinearStress_ && gradU)
+ (nut_*twoSymm(gradU) - nonlinearStress_) && gradU
);
- #include "nonLinearWallFunctionsI.H"
+ // Update epsilon and G at the wall
+ epsilon_.boundaryField().updateCoeffs();
// Dissipation equation
tmp<fvScalarMatrix> epsEqn
@@ -322,14 +288,13 @@ void ShihQuadraticKE::correct()
epsEqn().relax();
- #include "wallDissipationI.H"
+ epsEqn().boundaryManipulate(epsilon_.boundaryField());
solve(epsEqn);
bound(epsilon_, epsilonMin_);
// Turbulent kinetic energy equation
-
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(k_)
@@ -345,14 +310,7 @@ void ShihQuadraticKE::correct()
bound(k_, kMin_);
- // Re-calculate viscosity
-
- eta_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU + T(gradU))));
- ksi_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU - T(gradU))));
- Cmu_ = 2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_));
- fEta_ = A2_ + pow(eta_, 3.0);
-
- correctNut();
+ // Re-calculate viscosity and non-linear stress
correctNonlinearStress(gradU);
}
diff --git a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.H b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.H
index 01c6fe83c9d961b04ea50140059020f4a4100032..2348c638917343bf93017b7f7503a68f1f0b6ae0 100644
--- a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.H
+++ b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.H
@@ -77,15 +77,12 @@ protected:
dimensionedScalar C2_;
dimensionedScalar sigmak_;
dimensionedScalar sigmaEps_;
+ dimensionedScalar Cmu1_;
+ dimensionedScalar Cmu2_;
dimensionedScalar A1_;
- dimensionedScalar A2_;
- dimensionedScalar Ctau1_;
- dimensionedScalar Ctau2_;
- dimensionedScalar Ctau3_;
- dimensionedScalar alphaKsi_;
-
- dimensionedScalar kappa_;
- dimensionedScalar E_;
+ dimensionedScalar beta1_;
+ dimensionedScalar beta2_;
+ dimensionedScalar beta3_;
// Fields
@@ -93,11 +90,6 @@ protected:
volScalarField k_;
volScalarField epsilon_;
- volScalarField eta_;
- volScalarField ksi_;
- volScalarField Cmu_;
- volScalarField fEta_;
-
// Protected Member Functions
diff --git a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/include/nonLinearWallFunctionsI.H b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/include/nonLinearWallFunctionsI.H
deleted file mode 100644
index c9a4feb16b66c244f88f3cb162e3574e1bb500b8..0000000000000000000000000000000000000000
--- a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/include/nonLinearWallFunctionsI.H
+++ /dev/null
@@ -1,137 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
- \\/ 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 3 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, see <http://www.gnu.org/licenses/>.
-
-Global
- nonLinearwallFunctions
-
-Description
- Calculate wall generation and dissipation from wall-functions
- for non-linear models.
-
-\*---------------------------------------------------------------------------*/
-
-{
- labelList cellBoundaryFaceCount(epsilon_.size(), 0);
-
- scalar yPlusLam = nutkWallFunctionFvPatchScalarField::yPlusLam
- (
- kappa_.value(),
- E_.value()
- );
-
- const fvPatchList& patches = mesh_.boundary();
-
- //- Initialise the near-wall G and epsilon fields to zero
- forAll(patches, patchi)
- {
- const fvPatch& curPatch = patches[patchi];
-
- if (isA<wallFvPatch>(curPatch))
- {
- forAll(curPatch, facei)
- {
- label faceCelli = curPatch.faceCells()[facei];
-
- epsilon_[faceCelli] = 0.0;
- G[faceCelli] = 0.0;
- }
- }
- }
-
- const volScalarField nuLam(this->nu());
-
- //- Accumulate the wall face contributions to epsilon and G
- // Increment cellBoundaryFaceCount for each face for averaging
- forAll(patches, patchi)
- {
- const fvPatch& curPatch = patches[patchi];
-
- if (isA<wallFvPatch>(curPatch))
- {
- #include "checkPatchFieldTypes.H"
-
- const scalarField& nuw = nuLam.boundaryField()[patchi];
- const scalarField& nutw = nut_.boundaryField()[patchi];
-
- const scalarField magFaceGradU
- (
- mag(U_.boundaryField()[patchi].snGrad())
- );
-
- forAll(curPatch, facei)
- {
- label faceCelli = curPatch.faceCells()[facei];
-
- //- Using local Cmu !
- scalar Cmu25 = pow025(Cmu_[faceCelli]);
- scalar Cmu75 = pow(Cmu_[faceCelli], 0.75);
-
- scalar yPlus =
- Cmu25*y_[patchi][facei]
- *sqrt(k_[faceCelli])
- /nuw[facei];
-
- // For corner cells (with two boundary or more faces),
- // epsilon and G in the near-wall cell are calculated
- // as an average
-
- cellBoundaryFaceCount[faceCelli]++;
-
- epsilon_[faceCelli] +=
- Cmu75*pow(k_[faceCelli], 1.5)
- /(kappa_.value()*y_[patchi][facei]);
-
- if (yPlus > yPlusLam)
- {
- G[faceCelli] +=
- (nutw[facei] + nuw[facei])
- *magFaceGradU[facei]
- *Cmu25*sqrt(k_[faceCelli])
- /(kappa_.value()*y_[patchi][facei])
- - (nonlinearStress_[faceCelli] && gradU[faceCelli]);
- }
- }
- }
- }
-
- // Perform the averaging
-
- forAll(patches, patchi)
- {
- const fvPatch& curPatch = patches[patchi];
-
- if (isA<wallFvPatch>(curPatch))
- {
- forAll(curPatch, facei)
- {
- label faceCelli = curPatch.faceCells()[facei];
-
- epsilon_[faceCelli] /= cellBoundaryFaceCount[faceCelli];
- G[faceCelli] /= cellBoundaryFaceCount[faceCelli];
- }
- }
- }
-}
-
-
-// ************************************************************************* //
diff --git a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/include/wallNonlinearViscosityI.H b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/include/wallNonlinearViscosityI.H
deleted file mode 100644
index b6936b49c2e234cdd634bcb3a63a79f4648b92c1..0000000000000000000000000000000000000000
--- a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/include/wallNonlinearViscosityI.H
+++ /dev/null
@@ -1,78 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
- \\/ 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 3 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, see <http://www.gnu.org/licenses/>.
-
-Global
- wallNonlinearViscosity
-
-Description
- Calculate wall viscosity for non-linear models
-
-\*---------------------------------------------------------------------------*/
-
-{
- const fvPatchList& patches = mesh_.boundary();
-
- const scalar yPlusLam = nutkWallFunctionFvPatchScalarField::yPlusLam
- (
- kappa_.value(),
- E_.value()
- );
-
- const volScalarField nuLam(this->nu());
-
- forAll(patches, patchi)
- {
- const fvPatch& curPatch = patches[patchi];
-
- if (isA<wallFvPatch>(curPatch))
- {
- const scalarField& nuw = nuLam.boundaryField()[patchi];
- scalarField& nutw = nut_.boundaryField()[patchi];
-
- forAll(curPatch, facei)
- {
- label faceCelli = curPatch.faceCells()[facei];
-
- //- Using local Cmu
- scalar Cmu25 = pow025(Cmu_[faceCelli]);
-
- scalar yPlus =
- Cmu25*y_[patchi][facei]*sqrt(k_[faceCelli])/nuw[facei];
-
- if (yPlus > yPlusLam)
- {
- nutw[facei] =
- nuw[facei]
- *(yPlus*kappa_.value()/log(E_.value()*yPlus) - 1);
- }
- else
- {
- nutw[facei] = 0.0;
- }
- }
- }
- }
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/pitzDaily/constant/turbulenceProperties b/tutorials/incompressible/simpleFoam/pitzDaily/constant/turbulenceProperties
index 1c97bd5f78af9408c84a4fad9110acc20f731980..520cb4e3d2d7ba867cb2097e3b40e013993be1c4 100644
--- a/tutorials/incompressible/simpleFoam/pitzDaily/constant/turbulenceProperties
+++ b/tutorials/incompressible/simpleFoam/pitzDaily/constant/turbulenceProperties
@@ -19,7 +19,9 @@ simulationType RAS;
RAS
{
- RASModel v2f; //kEpsilon; //realizableKE; //kOmega; //kOmegaSST;
+ // Tested with ShihQuadraticKE, v2f, kEpsilon, realizableKE,
+ // kOmega, kOmegaSST;
+ RASModel v2f;
turbulence on;
diff --git a/tutorials/incompressible/simpleFoam/pitzDaily/system/fvSchemes b/tutorials/incompressible/simpleFoam/pitzDaily/system/fvSchemes
index 82b90f408ce987249a5bcf642ce259857fba212c..c91475d4264df7231f765e21f8cebbead8a18f44 100644
--- a/tutorials/incompressible/simpleFoam/pitzDaily/system/fvSchemes
+++ b/tutorials/incompressible/simpleFoam/pitzDaily/system/fvSchemes
@@ -34,6 +34,7 @@ divSchemes
div(phi,omega) bounded Gauss limitedLinear 1;
div(phi,v2) bounded Gauss limitedLinear 1;
div((nuEff*dev2(T(grad(U))))) Gauss linear;
+ div(nonlinearStress) Gauss linear;
}
laplacianSchemes