Merge branch 'feature-vof-turbulence-correction' into 'develop'

New VOF multiphaseStabilizedTurbulence fvOption

See merge request Development/openfoam!316

src/TurbulenceModels/turbulenceModels/RAS/RNGkEpsilon/RNGkEpsilon.C

@@ -251,23 +251,30 @@ void RNGkEpsilon<BasicTurbulenceModel>::correct()
     const rhoField& rho = this->rho_;
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
-    volScalarField& nut = this->nut_;
+    const volScalarField::Internal unlimitedNut(Cmu_*sqr(k_())/epsilon_());
+    const volScalarField& nut = this->nut_;
     fv::options& fvOptions(fv::options::New(this->mesh_));
 
     eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
 
-    volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));
+    const volScalarField::Internal divU
+    (
+        fvc::div(fvc::absolute(this->phi(), U))().v()
+    );
 
     tmp<volTensorField> tgradU = fvc::grad(U);
-    volScalarField S2((tgradU() && dev(twoSymm(tgradU()))));
+    const volScalarField::Internal GbyNu
+    (
+        tgradU().v() && dev(twoSymm(tgradU().v()))
+    );
     tgradU.clear();
 
-    volScalarField G(this->GName(), nut*S2);
+    const volScalarField::Internal G(this->GName(), unlimitedNut*GbyNu);
 
-    volScalarField eta(sqrt(mag(S2))*k_/epsilon_);
+    const volScalarField::Internal eta(sqrt(mag(GbyNu))*k_/epsilon_);
-    volScalarField eta3(eta*sqr(eta));
+    const volScalarField::Internal eta3(eta*sqr(eta));
 
-    volScalarField R
+    const volScalarField::Internal R
     (
         ((eta*(-eta/eta0_ + scalar(1)))/(beta_*eta3 + scalar(1)))
     );
@@ -282,9 +289,9 @@ void RNGkEpsilon<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, epsilon_)
       - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
      ==
-        (C1_ - R)*alpha*rho*G*epsilon_/k_
+        (C1_ - R)*alpha()*rho()*G*epsilon_()/k_()
-      - fvm::SuSp(((2.0/3.0)*C1_ - C3_)*alpha*rho*divU, epsilon_)
+      - fvm::SuSp(((2.0/3.0)*C1_ - C3_)*alpha()*rho()*divU, epsilon_)
-      - fvm::Sp(C2_*alpha*rho*epsilon_/k_, epsilon_)
+      - fvm::Sp(C2_*alpha()*rho()*epsilon_()/k_(), epsilon_)
       + epsilonSource()
       + fvOptions(alpha, rho, epsilon_)
     );
@@ -305,9 +312,9 @@ void RNGkEpsilon<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, k_)
       - fvm::laplacian(alpha*rho*DkEff(), k_)
      ==
-        alpha*rho*G
+        alpha()*rho()*GbyNu*nut()
-      - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
+      - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
-      - fvm::Sp(alpha*rho*epsilon_/k_, k_)
+      - fvm::Sp(alpha()*rho()*epsilon_()/k_(), k_)
       + kSource()
       + fvOptions(alpha, rho, k_)
     );

src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C

@@ -231,22 +231,24 @@ void kEpsilon<BasicTurbulenceModel>::correct()
     const rhoField& rho = this->rho_;
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
-    volScalarField& nut = this->nut_;
+    const volScalarField::Internal unlimitedNut(Cmu_*sqr(k_())/epsilon_());
+    const volScalarField& nut = this->nut_;
+
     fv::options& fvOptions(fv::options::New(this->mesh_));
 
     eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
 
-    volScalarField::Internal divU
+    const volScalarField::Internal divU
     (
         fvc::div(fvc::absolute(this->phi(), U))().v()
     );
 
     tmp<volTensorField> tgradU = fvc::grad(U);
-    volScalarField::Internal G
+    const volScalarField::Internal GbyNu
     (
-        this->GName(),
+        tgradU().v() && dev(twoSymm(tgradU().v()))
-        nut.v()*(dev(twoSymm(tgradU().v())) && tgradU().v())
     );
+    const volScalarField::Internal G(this->GName(), unlimitedNut*GbyNu);
     tgradU.clear();
 
     // Update epsilon and G at the wall
@@ -280,7 +282,7 @@ void kEpsilon<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, k_)
       - fvm::laplacian(alpha*rho*DkEff(), k_)
      ==
-        alpha()*rho()*G
+        alpha()*rho()*GbyNu*nut()
       - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
       - fvm::Sp(alpha()*rho()*epsilon_()/k_(), k_)
       + kSource()

src/TurbulenceModels/turbulenceModels/RAS/kOmega/kOmega.C

@@ -191,19 +191,23 @@ void kOmega<BasicTurbulenceModel>::correct()
     const rhoField& rho = this->rho_;
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
-    volScalarField& nut = this->nut_;
+    const volScalarField::Internal unlimitedNut(k_()/omega_());
+    const volScalarField& nut = this->nut_;
     fv::options& fvOptions(fv::options::New(this->mesh_));
 
     eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
 
-    volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));
+    const volScalarField::Internal divU
+    (
+        fvc::div(fvc::absolute(this->phi(), U))().v()
+    );
 
     tmp<volTensorField> tgradU = fvc::grad(U);
-    volScalarField G
+    const volScalarField::Internal GbyNu
     (
-        this->GName(),
+        tgradU().v() && dev(twoSymm(tgradU().v()))
-        nut*(tgradU() && dev(twoSymm(tgradU())))
     );
+    const volScalarField::Internal G(this->GName(), unlimitedNut*GbyNu);
     tgradU.clear();
 
     // Update omega and G at the wall
@@ -216,9 +220,9 @@ void kOmega<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, omega_)
       - fvm::laplacian(alpha*rho*DomegaEff(), omega_)
      ==
-        gamma_*alpha*rho*G*omega_/k_
+        gamma_*alpha()*rho()*G*omega_()/k_()
-      - fvm::SuSp(((2.0/3.0)*gamma_)*alpha*rho*divU, omega_)
+      - fvm::SuSp(((2.0/3.0)*gamma_)*alpha()*rho()*divU, omega_)
-      - fvm::Sp(beta_*alpha*rho*omega_, omega_)
+      - fvm::Sp(beta_*alpha()*rho()*omega_(), omega_)
       + fvOptions(alpha, rho, omega_)
     );
 
@@ -237,9 +241,9 @@ void kOmega<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, k_)
       - fvm::laplacian(alpha*rho*DkEff(), k_)
      ==
-        alpha*rho*G
+        alpha()*rho()*GbyNu*nut()
-      - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
+      - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
-      - fvm::Sp(Cmu_*alpha*rho*omega_, k_)
+      - fvm::Sp(Cmu_*alpha()*rho()*omega_(), k_)
       + fvOptions(alpha, rho, k_)
     );
 

src/finiteVolume/cfdTools/general/fvOptions/fvOption.C

@@ -58,9 +58,10 @@ Foam::fv::option::option
     coeffs_(dict.optionalSubDict(modelType + "Coeffs")),
     active_(dict_.lookupOrDefault<Switch>("active", true)),
     fieldNames_(),
-    applied_()
+    applied_(),
+    log(true)
 {
-    Info<< incrIndent << indent << "Source: " << name_ << endl << decrIndent;
+    Log << incrIndent << indent << "Source: " << name_ << endl << decrIndent;
 }
 
 

src/finiteVolume/cfdTools/general/fvOptions/fvOption.H

@@ -6,7 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
--------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
 
@@ -102,6 +103,8 @@ public:
     //- Runtime type information
     TypeName("option");
 
+    //- Switch write log to Info
+    bool log;
 
     // Declare run-time constructor selection table
 

src/finiteVolume/cfdTools/general/fvOptions/fvOptionIO.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,6 +56,8 @@ bool Foam::fv::option::read(const dictionary& dict)
 {
     dict.readIfPresent("active", active_);
 
+    log = dict.getOrDefault("log", true);
+
     coeffs_ = dict.optionalSubDict(modelType_ + "Coeffs");
 
     return true;

src/fvOptions/Make/files

@@ -18,6 +18,7 @@ $(derivedSources)/explicitPorositySource/explicitPorositySource.C
 $(derivedSources)/jouleHeatingSource/jouleHeatingSource.C
 $(derivedSources)/meanVelocityForce/meanVelocityForce.C
 $(derivedSources)/meanVelocityForce/patchMeanVelocityForce/patchMeanVelocityForce.C
+$(derivedSources)/multiphaseStabilizedTurbulence/multiphaseStabilizedTurbulence.C
 $(derivedSources)/phaseLimitStabilization/phaseLimitStabilization.C
 $(derivedSources)/radialActuationDiskSource/radialActuationDiskSource.C
 $(derivedSources)/rotorDiskSource/rotorDiskSource.C
@@ -57,4 +58,5 @@ $(derivedConstraints)/velocityDampingConstraint/velocityDampingConstraint.C
 corrections/limitTemperature/limitTemperature.C
 corrections/limitVelocity/limitVelocity.C
 
+
 LIB = $(FOAM_LIBBIN)/libfvOptions

src/fvOptions/sources/derived/multiphaseStabilizedTurbulence/multiphaseStabilizedTurbulence.C

@@ -0,0 +1,241 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "multiphaseStabilizedTurbulence.H"
+#include "fvMatrices.H"
+#include "turbulentTransportModel.H"
+#include "gravityMeshObject.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+    defineTypeNameAndDebug(multiphaseStabilizedTurbulence, 0);
+
+    addToRunTimeSelectionTable
+    (
+        option,
+        multiphaseStabilizedTurbulence,
+        dictionary
+    );
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fv::multiphaseStabilizedTurbulence::multiphaseStabilizedTurbulence
+(
+    const word& sourceName,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    option(sourceName, modelType, dict, mesh),
+    rhoName_(coeffs_.getOrDefault<word>("rho", "rho")),
+    Cmu_
+    (
+        dimensionedScalar::lookupOrAddToDict
+        (
+            "Cmu",
+            coeffs_,
+            0.09
+        )
+    ),
+    C_
+    (
+        dimensionedScalar::lookupOrAddToDict
+        (
+            "C",
+            coeffs_,
+            1.51
+        )
+    ),
+    lambda2_
+    (
+        dimensionedScalar::lookupOrAddToDict
+        (
+            "lambda2",
+            coeffs_,
+            0.1
+        )
+    ),
+    alpha_
+    (
+        dimensionedScalar::lookupOrAddToDict
+        (
+            "alpha",
+            coeffs_,
+            1.36
+        )
+    )
+{
+    fieldNames_.setSize(2, "undefined");
+
+    // Note: incompressible only
+    const auto* turbPtr =
+        mesh_.findObject<incompressible::turbulenceModel>
+        (
+            turbulenceModel::propertiesName
+        );
+
+    if (turbPtr)
+    {
+        const tmp<volScalarField>& tk = turbPtr->k();
+        fieldNames_[0] = tk().name();
+
+        const tmp<volScalarField>& tnut = turbPtr->nut();
+        fieldNames_[1] = tnut().name();
+
+        Log << "    Applying model to " << fieldNames_[0]
+            << " and " << fieldNames_[1] << endl;
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Unable to find incompressible turbulence model"
+            << exit(FatalError);
+    }
+
+    applied_.setSize(fieldNames_.size(), false);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::fv::multiphaseStabilizedTurbulence::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<scalar>& eqn,
+    const label fieldi
+)
+{
+    // Not applicable to compressible cases
+    NotImplemented;
+}
+
+
+void Foam::fv::multiphaseStabilizedTurbulence::addSup
+(
+    fvMatrix<scalar>& eqn,
+    const label fieldi
+)
+{
+    if (fieldi != 0)
+    {
+        return;
+    }
+
+    Log << this->name() << ": applying bouyancy production term to "
+        << eqn.psi().name() << endl;
+
+    // Buoyancy production in k eqn
+
+    const auto* turbPtr =
+        mesh_.findObject<incompressible::turbulenceModel>
+        (
+            turbulenceModel::propertiesName
+        );
+
+    if (!turbPtr)
+    {
+        FatalErrorInFunction
+            << "Unable to find incompressible turbulence model"
+            << exit(FatalError);
+    }
+
+
+    tmp<volScalarField> tepsilon = turbPtr->epsilon();
+    const volScalarField& epsilon = tepsilon();
+    const volScalarField& k = eqn.psi();
+
+    // Note: using solver density field for incompressible multiphase cases
+    const auto& rho = mesh_.lookupObject<volScalarField>(rhoName_);
+
+    const auto& g = meshObjects::gravity::New(mesh_.time());
+
+    const dimensionedScalar eps0("eps0", epsilon.dimensions(), SMALL);
+
+    // Note: differing from reference by replacing nut/k by Cmu*k/epsilon
+    const volScalarField GbyK
+    (
+        "GbyK",
+        Cmu_*k/(epsilon + eps0)*alpha_*(g & fvc::grad(rho))/rho
+    );
+
+    return eqn -= fvm::SuSp(GbyK, k);
+}
+
+
+void Foam::fv::multiphaseStabilizedTurbulence::correct(volScalarField& field)
+{
+    if (field.name() != fieldNames_[1])
+    {
+        return;
+    }
+
+    Log << this->name() << ": correcting " << field.name() << endl;
+
+    const auto* turbPtr =
+        mesh_.findObject<incompressible::turbulenceModel>
+        (
+            turbulenceModel::propertiesName
+        );
+
+    // nut correction
+    const auto& U = turbPtr->U();
+    tmp<volScalarField> tepsilon = turbPtr->epsilon();
+    const auto& epsilon = tepsilon();
+    tmp<volScalarField> tk = turbPtr->k();
+    const auto& k = tk();
+
+    tmp<volTensorField> tgradU = fvc::grad(U);
+    const auto& gradU = tgradU();
+    const dimensionedScalar pSmall("pSmall", dimless/sqr(dimTime), SMALL);
+    const volScalarField pRat
+    (
+        magSqr(symm(gradU))/(magSqr(skew(gradU)) + pSmall)
+    );
+
+    const volScalarField epsilonTilde
+    (
+        max
+        (
+            epsilon,
+            lambda2_*C_*pRat*epsilon
+        )
+    );
+
+    field = Cmu_*sqr(k)/epsilonTilde;
+    field.correctBoundaryConditions();
+}
+
+
+// ************************************************************************* //

src/fvOptions/sources/derived/multiphaseStabilizedTurbulence/multiphaseStabilizedTurbulence.H

@@ -0,0 +1,196 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+    Foam::fv::multiphaseStabilizedTurbulence
+
+Group
+    grpFvOptionsSources
+
+Description
+    Applies corrections to the turbulence kinetic energy equation and turbulence
+    viscosity field for incompressible multiphase flow cases.
+
+    Turbulence kinetic energy is over-predicted in VOF solvers at the phase
+    interface and throughout the water column in nearly-potential flow regions
+    beneath surface waves.
+
+    This fvOption applies corrections based on the references:
+    \verbatim
+        Buoyancy source term in turbulence kinetic energy equation:
+            Devolder, B., Rauwoens, P., and Troch, P. (2017).
+            Application of a buoyancy-modified k-w SST turbulence model to
+            simulate wave run-up around a monopile subjected to regular waves
+            using OpenFOAM.
+            Coastal Engineering, 125, 81-94.
+
+        Correction to turbulence viscosity:
+            Larsen, B.E. and Fuhrman, D.R. (2018).
+            On the over-production of turbulence beneath surface waves in
+            Reynolds-averaged Navier-Stokes models
+            J. Fluid Mech, 853, 419-460
+    \endverbatim
+
+Usage
+    Example usage:
+
+    \verbatim
+    multiphaseStabilizedTurbulence1
+    {
+        type            multiphaseStabilizedTurbulence;
+        active          yes;
+
+        multiphaseStabilizedTurbulenceCoeffs
+        {
+            // Optional coefficients
+            lambda2         0.1;   // A value of 0 sets the nut correction to 0
+            Cmu             0.09;  // from k-epsilon model
+            C               1.51;  // from k-omega model
+            alpha           1.36;  // 1/Prt
+        }
+    }
+    \endverbatim
+
+    The model C coefficient for the k-epsilon model equates to C2/C1 = 1.33;
+    the (default) value of 1.51 comes from the k-omega model and is more
+    conservative.
+
+SourceFiles
+    multiphaseStabilizedTurbulence.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fv_multiphaseStabilizedTurbulence_H
+#define fv_multiphaseStabilizedTurbulence_H
+
+#include "fvOption.H"
+#include "dimensionedScalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+               Class multiphaseStabilizedTurbulence Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiphaseStabilizedTurbulence
+:
+    public option
+{
+    // Private data
+
+        //- Name of density field
+        const word rhoName_;
+
+
+        // Model coefficients
+
+            //- k-epsilon model Cmu coefficient
+            dimensionedScalar Cmu_;
+
+            //- Model coefficient
+            //  For k-epsilon models this equates to C2/C1 = 1.33 and for
+            //  k-omega = 1.51.  Here the default is the more conservative 1.51
+            dimensionedScalar C_;
+
+            //- lambda2 coefficient; default = 0.1
+            dimensionedScalar lambda2_;
+
+            //- Buoyancy coefficient
+            dimensionedScalar alpha_;
+
+
+     // Private Member Functions
+
+        //- No copy construct
+        multiphaseStabilizedTurbulence
+        (
+            const multiphaseStabilizedTurbulence&
+        ) = delete;
+
+        //- No copy assignment
+        void operator=(const multiphaseStabilizedTurbulence&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("multiphaseStabilizedTurbulence");
+
+
+    // Constructors
+
+        //- Construct from explicit source name and mesh
+        multiphaseStabilizedTurbulence
+        (
+            const word& sourceName,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    // Member Functions
+
+        //- Add explicit contribution to compressible k equation
+        virtual void addSup
+        (
+            const volScalarField& rho,
+            fvMatrix<scalar>& eqn,
+            const label fieldi
+        );
+
+        //- Add explicit contribution to incompressible k equation
+        virtual void addSup
+        (
+            fvMatrix<scalar>& eqn,
+            const label fieldi
+        );
+
+        //- Correct the turbulence viscosity
+        virtual void correct(volScalarField& field);
+
+        //- Read source dictionary
+        virtual bool read(const dictionary& dict)
+        {
+            return true;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //