From 808ea384f0998c2619749f781aa898fdebedb6be Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Sun, 11 Aug 2013 10:57:43 +0100
Subject: [PATCH] 
 TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon:
 preliminary version

---
 .../RAS/mixtureKEpsilon/mixtureKEpsilon.C     | 698 ++++++++++++++++++
 .../RAS/mixtureKEpsilon/mixtureKEpsilon.H     | 252 +++++++
 2 files changed, 950 insertions(+)
 create mode 100644 src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C
 create mode 100644 src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.H

diff --git a/src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C b/src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C
new file mode 100644
index 00000000000..4cd773db4be
--- /dev/null
+++ b/src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C
@@ -0,0 +1,698 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "mixtureKEpsilon.H"
+#include "bound.H"
+#include "twoPhaseSystem.H"
+#include "fixedValueFvPatchFields.H"
+#include "inletOutletFvPatchFields.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+mixtureKEpsilon<BasicTurbulenceModel>::mixtureKEpsilon
+(
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName,
+    const word& type
+)
+:
+    eddyViscosity<RASModel<BasicTurbulenceModel> >
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaPhi,
+        phi,
+        transport,
+        propertiesName
+    ),
+
+    liquidTurbulencePtr_(NULL),
+
+    Cmu_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cmu",
+            this->coeffDict_,
+            0.09
+        )
+    ),
+    C1_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "C1",
+            this->coeffDict_,
+            1.44
+        )
+    ),
+    C2_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "C2",
+            this->coeffDict_,
+            1.92
+        )
+    ),
+    C3_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "C3",
+            this->coeffDict_,
+            C2_.value()
+        )
+    ),
+    Cp_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cp",
+            this->coeffDict_,
+            0.25
+        )
+    ),
+    sigmak_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "sigmak",
+            this->coeffDict_,
+            1.0
+        )
+    ),
+    sigmaEps_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "sigmaEps",
+            this->coeffDict_,
+            1.3
+        )
+    ),
+
+    k_
+    (
+        IOobject
+        (
+            IOobject::groupName("k", U.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    ),
+    epsilon_
+    (
+        IOobject
+        (
+            IOobject::groupName("epsilon", U.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    )
+{
+    bound(k_, this->kMin_);
+    bound(epsilon_, this->epsilonMin_);
+
+    if (type == typeName)
+    {
+        this->printCoeffs(type);
+    }
+}
+
+
+template<class BasicTurbulenceModel>
+wordList mixtureKEpsilon<BasicTurbulenceModel>::epsilonBoundaryTypes
+(
+    const volScalarField& epsilon
+) const
+{
+    const volScalarField::GeometricBoundaryField& ebf = epsilon.boundaryField();
+
+    wordList ebt = ebf.types();
+
+    forAll(ebf, patchi)
+    {
+        if (isA<fixedValueFvPatchScalarField>(ebf[patchi]))
+        {
+            ebt[patchi] = fixedValueFvPatchScalarField::typeName;
+        }
+    }
+
+    return ebt;
+}
+
+
+template<class BasicTurbulenceModel>
+void mixtureKEpsilon<BasicTurbulenceModel>::correctInletOutlet
+(
+    volScalarField& vsf,
+    const volScalarField& refVsf
+) const
+{
+    volScalarField::GeometricBoundaryField& bf = vsf.boundaryField();
+    const volScalarField::GeometricBoundaryField& refBf =
+        refVsf.boundaryField();
+
+    forAll(bf, patchi)
+    {
+        if
+        (
+            isA<inletOutletFvPatchScalarField>(bf[patchi])
+         && isA<inletOutletFvPatchScalarField>(refBf[patchi])
+        )
+        {
+            refCast<inletOutletFvPatchScalarField>
+            (bf[patchi]).refValue() =
+            refCast<const inletOutletFvPatchScalarField>
+            (refBf[patchi]).refValue();
+        }
+    }
+}
+
+
+template<class BasicTurbulenceModel>
+void mixtureKEpsilon<BasicTurbulenceModel>::initMixtureFields()
+{
+    if (rhom_.valid()) return;
+
+    // Local references to gas-phase properties
+    const volScalarField& kg = this->k_;
+    const volScalarField& epsilong = this->epsilon_;
+
+    // Local references to liquid-phase properties
+    mixtureKEpsilon<BasicTurbulenceModel>& turbc = this->liquidTurbulence();
+    const volScalarField& kl = turbc.k_;
+    const volScalarField& epsilonl = turbc.epsilon_;
+
+    Ct2_.set
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Ct2",
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            Ct2()
+        )
+    );
+
+    rhom_.set
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "rhom",
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            rhom()
+        )
+    );
+
+    km_.set
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "km",
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            mix(kl, kg),
+            kl.boundaryField().types()
+        )
+    );
+    correctInletOutlet(km_(), kl);
+
+    epsilonm_.set
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "epsilonm",
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            mix(epsilonl, epsilong),
+            epsilonBoundaryTypes(epsilonl)
+        )
+    );
+    correctInletOutlet(epsilonm_(), epsilonl);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+bool mixtureKEpsilon<BasicTurbulenceModel>::read()
+{
+    if (eddyViscosity<RASModel<BasicTurbulenceModel> >::read())
+    {
+        Cmu_.readIfPresent(this->coeffDict());
+        C1_.readIfPresent(this->coeffDict());
+        C2_.readIfPresent(this->coeffDict());
+        C3_.readIfPresent(this->coeffDict());
+        Cp_.readIfPresent(this->coeffDict());
+        sigmak_.readIfPresent(this->coeffDict());
+        sigmaEps_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class BasicTurbulenceModel>
+void mixtureKEpsilon<BasicTurbulenceModel>::correctNut()
+{
+    this->nut_ = Cmu_*sqr(k_)/epsilon_;
+    this->nut_.correctBoundaryConditions();
+}
+
+
+template<class BasicTurbulenceModel>
+mixtureKEpsilon<BasicTurbulenceModel>&
+mixtureKEpsilon<BasicTurbulenceModel>::liquidTurbulence() const
+{
+    if (!liquidTurbulencePtr_)
+    {
+        const volVectorField& U = this->U_;
+
+        const transportModel& gas = this->transport();
+        const twoPhaseSystem& fluid = gas.fluid();
+        const transportModel& liquid = fluid.otherPhase(gas);
+
+        liquidTurbulencePtr_ =
+           &const_cast<mixtureKEpsilon<BasicTurbulenceModel>&>
+            (
+                U.db().lookupObject<mixtureKEpsilon<BasicTurbulenceModel> >
+                (
+                    IOobject::groupName
+                    (
+                        turbulenceModel::propertiesName,
+                        liquid.name()
+                    )
+                )
+            );
+    }
+
+    return *liquidTurbulencePtr_;
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::Ct2() const
+{
+    const mixtureKEpsilon<BasicTurbulenceModel>& liquidTurbulence =
+        this->liquidTurbulence();
+
+    const transportModel& gas = this->transport();
+    const twoPhaseSystem& fluid = gas.fluid();
+    const transportModel& liquid = fluid.otherPhase(gas);
+
+    const volScalarField& alphag = this->alpha_;
+
+    volScalarField magUr(mag(liquidTurbulence.U() - this->U()));
+
+    volScalarField beta
+    (
+        (6*this->Cmu_/(4*sqrt(3.0/2.0)))
+       *alphag*fluid.drag(gas).K(magUr)/liquid.rho()
+       *(liquidTurbulence.k_/liquidTurbulence.epsilon_)
+    );
+    volScalarField Ct0((3 + beta)/(1 + beta + 2*gas.rho()/liquid.rho()));
+    volScalarField fAlphad((180 + (-4.71e3 + 4.26e4*alphag)*alphag)*alphag);
+
+    return sqr(1 + (Ct0 - 1)*exp(-fAlphad));
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::rholEff() const
+{
+    const transportModel& gas = this->transport();
+    const twoPhaseSystem& fluid = gas.fluid();
+    return fluid.otherPhase(gas).rho();
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::rhogEff() const
+{
+    const transportModel& gas = this->transport();
+    const twoPhaseSystem& fluid = gas.fluid();
+    return
+        gas.rho()
+      + fluid.Cvm()*fluid.otherPhase(gas).rho();
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::rhom() const
+{
+    const volScalarField& alphag = this->alpha_;
+    const volScalarField& alphal = this->liquidTurbulence().alpha_;
+
+    return alphal*rholEff() + alphag*rhogEff();
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::mix
+(
+    const volScalarField& fc,
+    const volScalarField& fd
+) const
+{
+    const volScalarField& alphag = this->alpha_;
+    const volScalarField& alphal = this->liquidTurbulence().alpha_;
+
+    return (alphal*rholEff()*fc + alphag*rhogEff()*fd)/rhom_();
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::mixU
+(
+    const volScalarField& fc,
+    const volScalarField& fd
+) const
+{
+    const volScalarField& alphag = this->alpha_;
+    const volScalarField& alphal = this->liquidTurbulence().alpha_;
+
+    return
+        (alphal*rholEff()*fc + alphag*rhogEff()*Ct2_()*fd)
+       /(alphal*rholEff() + alphag*rhogEff()*Ct2_());
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<surfaceScalarField> mixtureKEpsilon<BasicTurbulenceModel>::mixFlux
+(
+    const surfaceScalarField& fc,
+    const surfaceScalarField& fd
+) const
+{
+    const volScalarField& alphag = this->alpha_;
+    const volScalarField& alphal = this->liquidTurbulence().alpha_;
+
+    surfaceScalarField alphalf(fvc::interpolate(alphal));
+    surfaceScalarField alphagf(fvc::interpolate(alphag));
+
+    surfaceScalarField rholEfff(fvc::interpolate(rholEff()));
+    surfaceScalarField rhogEfff(fvc::interpolate(rhogEff()));
+
+    return
+       fvc::interpolate(rhom_()/(alphal*rholEff() + alphag*rhogEff()*Ct2_()))
+      *(alphalf*rholEfff*fc + alphagf*rhogEfff*fvc::interpolate(Ct2_())*fd);
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::bubbleG() const
+{
+    const mixtureKEpsilon<BasicTurbulenceModel>& liquidTurbulence =
+        this->liquidTurbulence();
+
+    const transportModel& gas = this->transport();
+    const twoPhaseSystem& fluid = gas.fluid();
+    const transportModel& liquid = fluid.otherPhase(gas);
+
+    volScalarField magUr(mag(liquidTurbulence.U() - this->U()));
+
+    // Lahey model
+    tmp<volScalarField> bubbleG
+    (
+        Cp_
+       *sqr(liquid)*liquid.rho()
+       *(
+            pow3(magUr)
+          + pow(fluid.drag(gas).K(magUr)*gas.d()/liquid.rho(), 3.0/4.0)
+           *pow(magUr, 9.0/4.0)
+        )
+       *gas
+       /gas.d()
+    );
+
+    // Simple model
+    // tmp<volScalarField> bubbleG
+    // (
+    //     Cp_*sqr(liquid)*gas*fluid.drag(gas).K(magUr)*sqr(magUr)
+    // );
+
+    return bubbleG;
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<fvScalarMatrix> mixtureKEpsilon<BasicTurbulenceModel>::kSource() const
+{
+    return fvm::Su(bubbleG(), km_());
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<fvScalarMatrix> mixtureKEpsilon<BasicTurbulenceModel>::epsilonSource() const
+{
+    return fvm::Su(C3_*epsilonm_()*bubbleG()/km_(), epsilonm_());
+}
+
+
+template<class BasicTurbulenceModel>
+void mixtureKEpsilon<BasicTurbulenceModel>::correct()
+{
+    const transportModel& gas = this->transport();
+    const twoPhaseSystem& fluid = gas.fluid();
+
+    // Only solve the mixture turbulence for the gas-phase
+    if (&gas != &fluid.phase1())
+    {
+        // This is the liquid phase but check the model for the gas-phase
+        // is consistent
+        this->liquidTurbulence();
+
+        return;
+    }
+
+    if (!this->turbulence_)
+    {
+        return;
+    }
+
+    // Initialise the mixture fields if they have not yet been constructed
+    initMixtureFields();
+
+    // Local references to gas-phase properties
+    const surfaceScalarField& phig = this->phi_;
+    const volVectorField& Ug = this->U_;
+    const volScalarField& alphag = this->alpha_;
+    volScalarField& kg = this->k_;
+    volScalarField& epsilong = this->epsilon_;
+    volScalarField& nutg = this->nut_;
+
+    // Local references to liquid-phase properties
+    mixtureKEpsilon<BasicTurbulenceModel>& liquidTurbulence =
+        this->liquidTurbulence();
+    const surfaceScalarField& phil = liquidTurbulence.phi_;
+    const volVectorField& Ul = liquidTurbulence.U_;
+    const volScalarField& alphal = liquidTurbulence.alpha_;
+    volScalarField& kl = liquidTurbulence.k_;
+    volScalarField& epsilonl = liquidTurbulence.epsilon_;
+    volScalarField& nutl = liquidTurbulence.nut_;
+
+    // Local references to mixture properties
+    volScalarField& rhom = rhom_();
+    volScalarField& km = km_();
+    volScalarField& epsilonm = epsilonm_();
+
+    eddyViscosity<RASModel<BasicTurbulenceModel> >::correct();
+
+    // Update the effective mixture density
+    rhom = this->rhom();
+
+    // Mixture flux
+    surfaceScalarField phim("phim", mixFlux(phil, phig));
+
+    // Mixture velocity divergence
+    volScalarField divUm
+    (
+        mixU
+        (
+            fvc::div(fvc::absolute(phil, Ul)),
+            fvc::div(fvc::absolute(phig, Ug))
+        )
+    );
+
+    tmp<volScalarField> Gc;
+    {
+        tmp<volTensorField> tgradUl = fvc::grad(Ul);
+        Gc = tmp<volScalarField>
+        (
+            new volScalarField
+            (
+                this->GName(),
+                nutl*(tgradUl() && dev(twoSymm(tgradUl())))
+            )
+        );
+        tgradUl.clear();
+
+        // Update k, epsilon and G at the wall
+        kl.boundaryField().updateCoeffs();
+        epsilonl.boundaryField().updateCoeffs();
+    }
+
+    tmp<volScalarField> Gd;
+    {
+        tmp<volTensorField> tgradUg = fvc::grad(Ug);
+        Gd = tmp<volScalarField>
+        (
+            new volScalarField
+            (
+                this->GName(),
+                nutg*(tgradUg() && dev(twoSymm(tgradUg())))
+            )
+        );
+        tgradUg.clear();
+
+        // Update k, epsilon and G at the wall
+        kg.boundaryField().updateCoeffs();
+        epsilong.boundaryField().updateCoeffs();
+    }
+
+    // Mixture turbulence generation
+    volScalarField Gm(mix(Gc, Gd));
+
+    // Mixture turbulence viscosity
+    volScalarField nutm(mixU(nutl, nutg));
+
+    // Update the mixture k and epsilon boundary conditions
+    km == mix(kl, kg);
+    bound(km, this->kMin_);
+    epsilonm == mix(epsilonl, epsilong);
+    bound(epsilonm, this->epsilonMin_);
+
+    // Dissipation equation
+    tmp<fvScalarMatrix> epsEqn
+    (
+        fvm::ddt(rhom, epsilonm)
+      + fvm::div(phim, epsilonm)
+      - fvm::Sp(fvc::ddt(rhom) + fvc::div(phim), epsilonm)
+      - fvm::laplacian(DepsilonEff(rhom*nutm), epsilonm)
+     ==
+        C1_*rhom*Gm*epsilonm/km
+      - fvm::SuSp(((2.0/3.0)*C1_)*rhom*divUm, epsilonm)
+      - fvm::Sp(C2_*rhom*epsilonm/km, epsilonm)
+      + epsilonSource()
+    );
+
+    epsEqn().relax();
+
+    epsEqn().boundaryManipulate(epsilonm.boundaryField());
+
+    solve(epsEqn);
+    bound(epsilonm, this->epsilonMin_);
+
+
+    // Turbulent kinetic energy equation
+    tmp<fvScalarMatrix> kmEqn
+    (
+        fvm::ddt(rhom, km)
+      + fvm::div(phim, km)
+      - fvm::Sp(fvc::ddt(rhom) + fvc::div(phim), km)
+      - fvm::laplacian(DkEff(rhom*nutm), km)
+     ==
+        rhom*Gm
+      - fvm::SuSp((2.0/3.0)*rhom*divUm, km)
+      - fvm::Sp(rhom*epsilonm/km, km)
+      + kSource()
+    );
+
+    kmEqn().relax();
+    solve(kmEqn);
+    bound(km, this->kMin_);
+    km.correctBoundaryConditions();
+
+    volScalarField Cc2(rhom/(alphal*rholEff() + alphag*rhogEff()*Ct2_()));
+    kl = Cc2*km;
+    kl.correctBoundaryConditions();
+    epsilonl = Cc2*epsilonm;
+    epsilonl.correctBoundaryConditions();
+    liquidTurbulence.correctNut();
+
+    Ct2_() = Ct2();
+    kg = Ct2_()*kl;
+    kg.correctBoundaryConditions();
+    epsilong = Ct2_()*epsilonl;
+    epsilong.correctBoundaryConditions();
+    nutg = Ct2_()*(liquidTurbulence.nu()/this->nu())*nutl;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.H b/src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.H
new file mode 100644
index 00000000000..ccc13b8035b
--- /dev/null
+++ b/src/TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.H
@@ -0,0 +1,252 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+Class
+    Foam::RASModels::mixtureKEpsilon
+
+Group
+    grpRASTurbulence
+
+Description
+    Standard k-epsilon turbulence model
+
+    The default model coefficients correspond to the following:
+    \verbatim
+        mixtureKEpsilonCoeffs
+        {
+            Cmu         0.09;
+            C1          1.44;
+            C2          1.92;
+            C3          -0.33;
+            sigmak      1.0;
+            sigmaEps    1.3;
+        }
+    \endverbatim
+
+SourceFiles
+    mixtureKEpsilon.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mixtureKEpsilon_H
+#define mixtureKEpsilon_H
+
+#include "RASModel.H"
+#include "eddyViscosity.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class mixtureKEpsilon Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicTurbulenceModel>
+class mixtureKEpsilon
+:
+    public eddyViscosity<RASModel<BasicTurbulenceModel> >
+{
+    // Private data
+
+        mutable mixtureKEpsilon<BasicTurbulenceModel> *liquidTurbulencePtr_;
+
+
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        mixtureKEpsilon(const mixtureKEpsilon&);
+        mixtureKEpsilon& operator=(const mixtureKEpsilon&);
+
+        //- Return the turbulence model for the other phase
+        mixtureKEpsilon<BasicTurbulenceModel>& liquidTurbulence() const;
+
+
+protected:
+
+    // Protected data
+
+        // Model coefficients
+
+            dimensionedScalar Cmu_;
+            dimensionedScalar C1_;
+            dimensionedScalar C2_;
+            dimensionedScalar C3_;
+            dimensionedScalar Cp_;
+            dimensionedScalar sigmak_;
+            dimensionedScalar sigmaEps_;
+
+        // Fields
+
+            volScalarField k_;
+            volScalarField epsilon_;
+
+        // Mixture fields
+
+            autoPtr<volScalarField> Ct2_;
+            autoPtr<volScalarField> rhom_;
+            autoPtr<volScalarField> km_;
+            autoPtr<volScalarField> epsilonm_;
+
+
+    // Protected Member Functions
+
+        wordList epsilonBoundaryTypes(const volScalarField& epsilon) const;
+
+        void correctInletOutlet
+        (
+            volScalarField& vsf,
+            const volScalarField& refVsf
+        ) const;
+
+        void initMixtureFields();
+
+        virtual void correctNut();
+
+        tmp<volScalarField> Ct2() const;
+
+        tmp<volScalarField> rholEff() const;
+        tmp<volScalarField> rhogEff() const;
+        tmp<volScalarField> rhom() const;
+
+        tmp<volScalarField> mix
+        (
+            const volScalarField& fc,
+            const volScalarField& fd
+        ) const;
+
+        tmp<volScalarField> mixU
+        (
+            const volScalarField& fc,
+            const volScalarField& fd
+        ) const;
+
+        tmp<surfaceScalarField> mixFlux
+        (
+            const surfaceScalarField& fc,
+            const surfaceScalarField& fd
+        ) const;
+
+        tmp<volScalarField> bubbleG() const;
+        virtual tmp<fvScalarMatrix> kSource() const;
+        virtual tmp<fvScalarMatrix> epsilonSource() const;
+
+        //- Return the effective diffusivity for k
+        tmp<volScalarField> DkEff(const volScalarField& nutm) const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    "DkEff",
+                    nutm/sigmak_
+                )
+            );
+        }
+
+        //- Return the effective diffusivity for epsilon
+        tmp<volScalarField> DepsilonEff(const volScalarField& nutm) const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    "DepsilonEff",
+                    nutm/sigmaEps_
+                )
+            );
+        }
+
+public:
+
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
+    typedef typename BasicTurbulenceModel::rhoField rhoField;
+    typedef typename BasicTurbulenceModel::transportModel transportModel;
+
+
+    //- Runtime type information
+    TypeName("mixtureKEpsilon");
+
+
+    // Constructors
+
+        //- Construct from components
+        mixtureKEpsilon
+        (
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+    //- Destructor
+    virtual ~mixtureKEpsilon()
+    {}
+
+
+    // Member Functions
+
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+
+        //- Return the turbulence kinetic energy
+        virtual tmp<volScalarField> k() const
+        {
+            return k_;
+        }
+
+        //- Return the turbulence kinetic energy dissipation rate
+        virtual tmp<volScalarField> epsilon() const
+        {
+            return epsilon_;
+        }
+
+        //- Solve the turbulence equations and correct the turbulence viscosity
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "mixtureKEpsilon.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
-- 
GitLab