From e2f568cd8495c9fd41f798f60c3ef6c33ba0f763 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 16 Mar 2016 17:54:03 +0000 Subject: [PATCH] New version of block-coupled k-omega SST model. Robert Keser --- .../RAS/blockKOmegaSST/blockKOmegaSST.C | 552 ++++++++++++++++++ .../RAS/blockKOmegaSST/blockKOmegaSST.H | 274 +++++++++ 2 files changed, 826 insertions(+) create mode 100644 src/turbulenceModels/incompressible/RAS/blockKOmegaSST/blockKOmegaSST.C create mode 100644 src/turbulenceModels/incompressible/RAS/blockKOmegaSST/blockKOmegaSST.H diff --git a/src/turbulenceModels/incompressible/RAS/blockKOmegaSST/blockKOmegaSST.C b/src/turbulenceModels/incompressible/RAS/blockKOmegaSST/blockKOmegaSST.C new file mode 100644 index 000000000..baab6c628 --- /dev/null +++ b/src/turbulenceModels/incompressible/RAS/blockKOmegaSST/blockKOmegaSST.C @@ -0,0 +1,552 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "blockKOmegaSST.H" +#include "fvBlockMatrix.H" +#include "addToRunTimeSelectionTable.H" + +#include "backwardsCompatibilityWallFunctions.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace incompressible +{ +namespace RASModels +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(blockKOmegaSST, 0); +addToRunTimeSelectionTable(RASModel, blockKOmegaSST, dictionary); + +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + +tmp blockKOmegaSST::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 blockKOmegaSST::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 blockKOmegaSST::F3() const +{ + tmp arg3 = min + ( + 150*nu()/(omega_*sqr(y_)), + scalar(10) + ); + + return 1 - tanh(pow4(arg3)); +} + + +tmp blockKOmegaSST::F23() const +{ + tmp f23(F2()); + + if (F3_) + { + f23() *= F3(); + } + + return f23; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +blockKOmegaSST::blockKOmegaSST +( + const volVectorField& U, + const surfaceScalarField& phi, + transportModel& lamTransportModel +) +: + RASModel(typeName, U, phi, lamTransportModel), + + alphaK1_ + ( + dimensionedScalar::lookupOrAddToDict + ( + "alphaK1", + coeffDict_, + 0.85034 + ) + ), + alphaK2_ + ( + dimensionedScalar::lookupOrAddToDict + ( + "alphaK2", + coeffDict_, + 1.0 + ) + ), + alphaOmega1_ + ( + dimensionedScalar::lookupOrAddToDict + ( + "alphaOmega1", + coeffDict_, + 0.5 + ) + ), + alphaOmega2_ + ( + dimensionedScalar::lookupOrAddToDict + ( + "alphaOmega2", + coeffDict_, + 0.856 + ) + ), + gamma1_ + ( + dimensionedScalar::lookupOrAddToDict + ( + "gamma1", + coeffDict_, + 5.0/9.0 + ) + ), + gamma2_ + ( + dimensionedScalar::lookupOrAddToDict + ( + "gamma2", + coeffDict_, + 0.44 + ) + ), + beta1_ + ( + dimensionedScalar::lookupOrAddToDict + ( + "beta1", + coeffDict_, + 0.075 + ) + ), + beta2_ + ( + dimensionedScalar::lookupOrAddToDict + ( + "beta2", + coeffDict_, + 0.0828 + ) + ), + betaStar_ + ( + dimensionedScalar::lookupOrAddToDict + ( + "betaStar", + coeffDict_, + 0.09 + ) + ), + a1_ + ( + dimensionedScalar::lookupOrAddToDict + ( + "a1", + coeffDict_, + 0.31 + ) + ), + b1_ + ( + dimensionedScalar::lookupOrAddToDict + ( + "b1", + coeffDict_, + 1.0 + ) + ), + c1_ + ( + dimensionedScalar::lookupOrAddToDict + ( + "c1", + coeffDict_, + 10.0 + ) + ), + F3_ + ( + Switch::lookupOrAddToDict + ( + "F3", + coeffDict_, + false + ) + ), + + y_(mesh_), + + k_ + ( + IOobject + ( + "k", + runTime_.timeName(), + U_.db(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + autoCreateK("k", mesh_, U_.db()) + ), + omega_ + ( + IOobject + ( + "omega", + runTime_.timeName(), + U_.db(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + autoCreateOmega("omega", mesh_, U_.db()) + ), + nut_ + ( + IOobject + ( + "nut", + runTime_.timeName(), + U_.db(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + autoCreateNut("nut", mesh_, U_.db()) + ), + kOmega_ + ( + IOobject + ( + "kOmega", + runTime_.timeName(), + U_.db(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedVector2("zero", dimless, vector2::zero) + ) +{ + bound(k_, k0_); + bound(omega_, omega0_); + + nut_ = + ( + a1_*k_/ + max + ( + a1_*omega_, + b1_*F23()*sqrt(2.0)*mag(symm(fvc::grad(U_))) + ) + ); + nut_.correctBoundaryConditions(); + + printCoeffs(); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +tmp blockKOmegaSST::R() const +{ + return tmp + ( + new volSymmTensorField + ( + IOobject + ( + "R", + runTime_.timeName(), + U_.db(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)), + k_.boundaryField().types() + ) + ); +} + + +tmp blockKOmegaSST::devReff() const +{ + return tmp + ( + new volSymmTensorField + ( + IOobject + ( + "devRhoReff", + runTime_.timeName(), + U_.db(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + -nuEff()*dev(twoSymm(fvc::grad(U_))) + ) + ); +} + + +tmp blockKOmegaSST::divDevReff(volVectorField& U) const +{ + return + ( + - fvm::laplacian(nuEff(), U) + - fvc::div(nuEff()*dev(fvc::grad(U)().T())) + ); +} + + +bool blockKOmegaSST::read() +{ + if (RASModel::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()); + b1_.readIfPresent(coeffDict()); + c1_.readIfPresent(coeffDict()); + F3_.readIfPresent("F3", coeffDict()); + + return true; + } + else + { + return false; + } +} + + +void blockKOmegaSST::correct() +{ + // Bound in case of topological change + // HJ, 22/Aug/2007 + if (mesh_.changing()) + { + bound(k_, k0_); + bound(omega_, omega0_); + } + + RASModel::correct(); + + if (!turbulence_) + { + return; + } + + if (mesh_.changing()) + { + y_.correct(); + } + + const volScalarField S2(2*magSqr(symm(fvc::grad(U_)))); + volScalarField G("RASModel::G", nut_*S2); + + // Make coupled matrix + fvBlockMatrix kOmegaEqn(kOmega_); + + // Update omega and G at the wall + omega_.boundaryField().updateCoeffs(); + + const volScalarField CDkOmega + ( + (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_ + ); + + const volScalarField F1(this->F1(CDkOmega)); + + // Record coupling equations when omega is set. Same equations + // will be eliminated from k coupling as well due to wall functions + // HJ, 13/Nov/2015 + labelList eliminated; + + // Turbulent frequency equation + { + fvScalarMatrix omegaEqn + ( + fvm::ddt(omega_) + + fvm::div(phi_, omega_) + + fvm::SuSp(-fvc::div(phi_), omega_) + - fvm::laplacian(DomegaEff(F1), omega_) + == + gamma(F1) + *min + ( + S2, + (c1_/a1_)*betaStar_*omega_*max(a1_*omega_, b1_*F23()*sqrt(S2)) + ) + - fvm::Sp(2*beta(F1)*omega_, omega_) + + beta(F1)*sqr(omega_) + - fvm::SuSp + ( + (F1 - scalar(1))*CDkOmega/omega_, + omega_ + ) + ); + + omegaEqn.relax(); + omegaEqn.completeAssembly(); + + kOmegaEqn.insertEquation(1, omegaEqn); + + // Add coupling term + volScalarField couplingOmega + ( + "couplingOmega", + -pos((scalar(1) - F1)*CDkOmega)*(scalar(1) - F1)*CDkOmega/k_ + ); + scalarField& couplingOmegaIn = couplingOmega.internalField(); + + // Collect cell indices where wall functions are active from the + // omega wall functions + eliminated = omegaEqn.eliminatedEqns().toc(); + + // Eliminate coupling in wall function cells for omega + forAll (eliminated, cellI) + { + couplingOmegaIn[eliminated[cellI]] *= 0; + } + + // Insert coupling + kOmegaEqn.insertEquationCoupling(1, 0, couplingOmega); + } + + // Turbulent kinetic energy equation + { + fvScalarMatrix kEqn + ( + fvm::ddt(k_) + + fvm::div(phi_, k_) + + fvm::SuSp(-fvc::div(phi_), k_) + - fvm::laplacian(DkEff(F1), k_) + == + min(G, c1_*betaStar_*k_*omega_) + - fvm::Sp(betaStar_*omega_, k_) + ); + + kEqn.relax(); + kEqn.completeAssembly(); + + kOmegaEqn.insertEquation(0, kEqn); + + // Add coupling term + volScalarField couplingK + ( + "couplingK", + -pos(G - c1_*betaStar_*k_*omega_)*c1_*betaStar_*k_ + ); + scalarField& couplingKIn = couplingK.internalField(); + + // Eliminate coupling in wall function cells for k + forAll (eliminated, cellI) + { + couplingKIn[eliminated[cellI]] *= 0; + } + + // Insert coupling + kOmegaEqn.insertEquationCoupling(0, 1, couplingK); + } + + // Update source coupling: coupling terms eliminated from source + kOmegaEqn.updateSourceCoupling(); + + kOmegaEqn.solve(); + + // Retrieve solution + kOmegaEqn.retrieveSolution(0, k_.internalField()); + kOmegaEqn.retrieveSolution(1, omega_.internalField()); + + bound(k_, k0_); + bound(omega_, omega0_); + + k_.correctBoundaryConditions(); + omega_.correctBoundaryConditions(); + + // Re-calculate viscosity + // Fixed sqrt(2) error. HJ, 10/Jun/2015 + nut_ = a1_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2)); + nut_ = min(nut_, nuRatio()*nu()); + nut_.correctBoundaryConditions(); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace RASModels +} // End namespace incompressible +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/turbulenceModels/incompressible/RAS/blockKOmegaSST/blockKOmegaSST.H b/src/turbulenceModels/incompressible/RAS/blockKOmegaSST/blockKOmegaSST.H new file mode 100644 index 000000000..255970377 --- /dev/null +++ b/src/turbulenceModels/incompressible/RAS/blockKOmegaSST/blockKOmegaSST.H @@ -0,0 +1,274 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Class + Foam::incompressible::RASModels::blockKOmegaSST + +Description + Implementation of the k-omega-SST turbulence model for incompressible + flows using a block-coupled solution. + + The default model coefficients correspond to the following: + @verbatim + blockKOmegaSSTCoeffs + { + alphaK1 0.85; + alphaK2 1.0; + alphaOmega1 0.5; + alphaOmega2 0.856; + beta1 0.075; + beta2 0.0828; + betaStar 0.09; + gamma1 5/9; + gamma2 0.44; + a1 0.31; + b1 1.0; + c1 10.0; + F3 no; + } + @endverbatim + +SourceFiles + blockKOmegaSST.C + +\*---------------------------------------------------------------------------*/ + +#ifndef blockKOmegaSST_H +#define blockKOmegaSST_H + +#include "RASModel.H" +#include "wallDist.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace incompressible +{ +namespace RASModels +{ + +/*---------------------------------------------------------------------------*\ + Class kOmega Declaration +\*---------------------------------------------------------------------------*/ + +class blockKOmegaSST +: + public RASModel +{ + // Private data + + // Model coefficients + dimensionedScalar alphaK1_; + dimensionedScalar alphaK2_; + + dimensionedScalar alphaOmega1_; + dimensionedScalar alphaOmega2_; + + dimensionedScalar gamma1_; + dimensionedScalar gamma2_; + + dimensionedScalar beta1_; + dimensionedScalar beta2_; + + dimensionedScalar betaStar_; + + dimensionedScalar a1_; + dimensionedScalar b1_; + dimensionedScalar c1_; + + Switch F3_; + + + //- Wall distance field + // Note: different to wall distance in parent RASModel + wallDist y_; + + // Fields + + volScalarField k_; + volScalarField omega_; + volScalarField nut_; + + // Block vector field for k-omega + volVector2Field kOmega_; + + + // Private member functions + + tmp F1(const volScalarField& CDkOmega) const; + tmp F2() const; + tmp F3() const; + tmp F23() const; + + tmp blend + ( + const volScalarField& F1, + const dimensionedScalar& psi1, + const dimensionedScalar& psi2 + ) const + { + return F1*(psi1 - psi2) + psi2; + } + + tmp alphaK + ( + const volScalarField& F1 + ) const + { + return blend(F1, alphaK1_, alphaK2_); + } + + tmp alphaOmega + ( + const volScalarField& F1 + ) const + { + return blend(F1, alphaOmega1_, alphaOmega2_); + } + + tmp beta + ( + const volScalarField& F1 + ) const + { + return blend(F1, beta1_, beta2_); + } + + tmp gamma + ( + const volScalarField& F1 + ) const + { + return blend(F1, gamma1_, gamma2_); + } + + +public: + + //- Runtime type information + TypeName("blockKOmegaSST"); + + + // Constructors + + //- Construct from components + blockKOmegaSST + ( + const volVectorField& U, + const surfaceScalarField& phi, + transportModel& transport + ); + + + //- Destructor + + virtual ~blockKOmegaSST() + {} + + + // Member Functions + + //- Return the turbulence viscosity + virtual tmp nut() const + { + return nut_; + } + + //- Return the effective diffusivity for k + tmp DkEff(const volScalarField& F1) const + { + return tmp + ( + new volScalarField("DkEff", alphaK(F1)*nut_ + nu()) + ); + } + + //- Return the effective diffusivity for omega + tmp DomegaEff(const volScalarField& F1) const + { + return tmp + ( + new volScalarField("DomegaEff", alphaOmega(F1)*nut_ + nu()) + ); + } + + //- Return the turbulence kinetic energy + virtual tmp k() const + { + return k_; + } + + //- Return the turbulence specific dissipation rate + virtual tmp omega() const + { + return omega_; + } + + //- Return the turbulence kinetic energy dissipation rate + virtual tmp epsilon() const + { + return tmp + ( + new volScalarField + ( + IOobject + ( + "epsilon", + mesh_.time().timeName(), + mesh_ + ), + betaStar_*k_*omega_, + omega_.boundaryField().types() + ) + ); + } + + //- Return the Reynolds stress tensor + virtual tmp R() const; + + //- Return the effective stress tensor including the laminar stress + virtual tmp devReff() const; + + //- Return the source term for the momentum equation + virtual tmp divDevReff(volVectorField& U) const; + + //- Solve the turbulence equations and correct the turbulence viscosity + virtual void correct(); + + //- Read RASProperties dictionary + virtual bool read(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace RASModels +} // namespace incompressible +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* //