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
+
+// ************************************************************************* //