twoPhaseEulerFoam: Corrected kineticTheory frictional stress at walls

Patch contributed by Juho Peltola, VTT

The new JohnsonJacksonSchaefferFrictionalStress model is included.

Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=2058

applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/Make/files

@@ -31,6 +31,7 @@ kineticTheoryModels/frictionalStressModel/frictionalStressModel/frictionalStress
 kineticTheoryModels/frictionalStressModel/frictionalStressModel/newFrictionalStressModel.C
 kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.C
 kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
+kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C
 
 kineticTheoryModels/derivedFvPatchFields/JohnsonJacksonParticleTheta/JohnsonJacksonParticleThetaFvPatchScalarField.C
 kineticTheoryModels/derivedFvPatchFields/JohnsonJacksonParticleSlip/JohnsonJacksonParticleSlipFvPatchVectorField.C

applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/derivedFvPatchFields/JohnsonJacksonParticleSlip/JohnsonJacksonParticleSlipFvPatchVectorField.C

@@ -190,6 +190,14 @@ void Foam::JohnsonJacksonParticleSlipFvPatchVectorField::updateCoeffs()
         )
     );
 
+    const scalarField nuFric
+    (
+        patch().lookupPatchField<volScalarField, scalar>
+        (
+            IOobject::groupName("nuFric", phased.name())
+        )
+    );
+
     word ThetaName(IOobject::groupName("Theta", phased.name()));
 
     const fvPatchScalarField& Theta
@@ -222,7 +230,7 @@ void Foam::JohnsonJacksonParticleSlipFvPatchVectorField::updateCoeffs()
        *gs0
        *specularityCoefficient_.value()
        *sqrt(3.0*Theta)
-       /max(6.0*nu*alphaMax.value(), SMALL)
+       /max(6.0*(nu - nuFric)*alphaMax.value(), SMALL)
     );
 
     this->valueFraction() = c/(c + patch().deltaCoeffs());

applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/derivedFvPatchFields/JohnsonJacksonParticleTheta/JohnsonJacksonParticleThetaFvPatchScalarField.C

@@ -65,7 +65,8 @@ JohnsonJacksonParticleThetaFvPatchScalarField
     mixedFvPatchScalarField(ptf, p, iF, mapper),
     restitutionCoefficient_(ptf.restitutionCoefficient_),
     specularityCoefficient_(ptf.specularityCoefficient_)
-{}
+{
+}
 
 
 Foam::JohnsonJacksonParticleThetaFvPatchScalarField::

applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,15 +81,16 @@ Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::frictionalStressModels::JohnsonJackson::
 frictionalPressure
 (
-    const volScalarField& alpha1,
+    const phaseModel& phase,
     const dimensionedScalar& alphaMinFriction,
     const dimensionedScalar& alphaMax
 ) const
 {
+    const volScalarField& alpha = phase;
 
     return
-        Fr_*pow(max(alpha1 - alphaMinFriction, scalar(0)), eta_)
+        Fr_*pow(max(alpha - alphaMinFriction, scalar(0)), eta_)
-       /pow(max(alphaMax - alpha1, alphaDeltaMin_), p_);
+       /pow(max(alphaMax - alpha, alphaDeltaMin_), p_);
 }
 
 
@@ -97,24 +98,26 @@ Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::frictionalStressModels::JohnsonJackson::
 frictionalPressurePrime
 (
-    const volScalarField& alpha1,
+    const phaseModel& phase,
     const dimensionedScalar& alphaMinFriction,
     const dimensionedScalar& alphaMax
 ) const
 {
+    const volScalarField& alpha = phase;
+
     return Fr_*
     (
-        eta_*pow(max(alpha1 - alphaMinFriction, scalar(0)), eta_ - 1.0)
+        eta_*pow(max(alpha - alphaMinFriction, scalar(0)), eta_ - 1.0)
-       *(alphaMax-alpha1)
+       *(alphaMax-alpha)
-      + p_*pow(max(alpha1 - alphaMinFriction, scalar(0)), eta_)
+      + p_*pow(max(alpha - alphaMinFriction, scalar(0)), eta_)
-    )/pow(max(alphaMax - alpha1, alphaDeltaMin_), p_ + 1.0);
+    )/pow(max(alphaMax - alpha, alphaDeltaMin_), p_ + 1.0);
 }
 
 
 Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::frictionalStressModels::JohnsonJackson::nu
 (
-    const volScalarField& alpha1,
+    const phaseModel& phase,
     const dimensionedScalar& alphaMinFriction,
     const dimensionedScalar& alphaMax,
     const volScalarField& pf,

applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -93,21 +93,21 @@ public:
 
         virtual tmp<volScalarField> frictionalPressure
         (
-            const volScalarField& alpha1,
+            const phaseModel& phase,
             const dimensionedScalar& alphaMinFriction,
             const dimensionedScalar& alphaMax
         ) const;
 
         virtual tmp<volScalarField> frictionalPressurePrime
         (
-            const volScalarField& alpha1,
+            const phaseModel& phase,
             const dimensionedScalar& alphaMinFriction,
             const dimensionedScalar& alphaMax
         ) const;
 
         virtual tmp<volScalarField> nu
         (
-            const volScalarField& alpha1,
+            const phaseModel& phase,
             const dimensionedScalar& alphaMinFriction,
             const dimensionedScalar& alphaMax,
             const volScalarField& pf,

applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C

@@ -0,0 +1,208 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 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 "JohnsonJacksonSchaefferFrictionalStress.H"
+#include "addToRunTimeSelectionTable.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace kineticTheoryModels
+{
+namespace frictionalStressModels
+{
+    defineTypeNameAndDebug(JohnsonJacksonSchaeffer, 0);
+
+    addToRunTimeSelectionTable
+    (
+        frictionalStressModel,
+        JohnsonJacksonSchaeffer,
+        dictionary
+    );
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::kineticTheoryModels::frictionalStressModels::
+JohnsonJacksonSchaeffer::JohnsonJacksonSchaeffer
+(
+    const dictionary& dict
+)
+:
+    frictionalStressModel(dict),
+    coeffDict_(dict.subDict(typeName + "Coeffs")),
+    Fr_("Fr", dimensionSet(1, -1, -2, 0, 0), coeffDict_),
+    eta_("eta", dimless, coeffDict_),
+    p_("p", dimless, coeffDict_),
+    phi_("phi", dimless, coeffDict_),
+    alphaDeltaMin_("alphaDeltaMin", dimless, coeffDict_)
+{
+    phi_ *= constant::mathematical::pi/180.0;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::kineticTheoryModels::frictionalStressModels::
+JohnsonJacksonSchaeffer::~JohnsonJacksonSchaeffer()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::kineticTheoryModels::frictionalStressModels::
+JohnsonJacksonSchaeffer::frictionalPressure
+(
+    const phaseModel& phase,
+    const dimensionedScalar& alphaMinFriction,
+    const dimensionedScalar& alphaMax
+) const
+{
+    const volScalarField& alpha = phase;
+
+    return
+        Fr_*pow(max(alpha - alphaMinFriction, scalar(0)), eta_)
+       /pow(max(alphaMax - alpha, alphaDeltaMin_), p_);
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::kineticTheoryModels::frictionalStressModels::
+JohnsonJacksonSchaeffer::frictionalPressurePrime
+(
+    const phaseModel& phase,
+    const dimensionedScalar& alphaMinFriction,
+    const dimensionedScalar& alphaMax
+) const
+{
+    const volScalarField& alpha = phase;
+
+    return Fr_*
+    (
+        eta_*pow(max(alpha - alphaMinFriction, scalar(0)), eta_ - 1.0)
+       *(alphaMax-alpha)
+      + p_*pow(max(alpha - alphaMinFriction, scalar(0)), eta_)
+    )/pow(max(alphaMax - alpha, alphaDeltaMin_), p_ + 1.0);
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::kineticTheoryModels::frictionalStressModels::
+JohnsonJacksonSchaeffer::nu
+(
+    const phaseModel& phase,
+    const dimensionedScalar& alphaMinFriction,
+    const dimensionedScalar& alphaMax,
+    const volScalarField& pf,
+    const volSymmTensorField& D
+) const
+{
+    const volScalarField& alpha = phase;
+
+    tmp<volScalarField> tnu
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "JohnsonJacksonSchaeffer:nu",
+                phase.mesh().time().timeName(),
+                phase.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            phase.mesh(),
+            dimensionedScalar("nu", dimensionSet(0, 2, -1, 0, 0), 0.0)
+        )
+    );
+
+    volScalarField& nuf = tnu.ref();
+
+    forAll(D, celli)
+    {
+        if (alpha[celli] > alphaMinFriction.value())
+        {
+            nuf[celli] =
+                0.5*pf[celli]*sin(phi_.value())
+               /(
+                    sqrt((1.0/3.0)*sqr(tr(D[celli])) - invariantII(D[celli]))
+                  + SMALL
+                );
+        }
+    }
+
+    const fvPatchList& patches = phase.mesh().boundary();
+    const volVectorField& U = phase.U();
+
+    volScalarField::Boundary& nufBf = nuf.boundaryFieldRef();
+
+    forAll(patches, patchi)
+    {
+        if (!patches[patchi].coupled())
+        {
+            nufBf[patchi] =
+                (
+                    pf.boundaryField()[patchi]*sin(phi_.value())
+                   /(
+                        mag(U.boundaryField()[patchi].snGrad())
+                      + SMALL
+                    )
+                );
+        }
+    }
+
+    // Correct coupled BCs
+    nuf.correctBoundaryConditions();
+
+    return tnu;
+}
+
+
+bool Foam::kineticTheoryModels::frictionalStressModels::
+JohnsonJacksonSchaeffer::read()
+{
+    coeffDict_ <<= dict_.subDict(typeName + "Coeffs");
+
+    Fr_.read(coeffDict_);
+    eta_.read(coeffDict_);
+    p_.read(coeffDict_);
+
+    phi_.read(coeffDict_);
+    phi_ *= constant::mathematical::pi/180.0;
+
+    alphaDeltaMin_.read(coeffDict_);
+
+    return true;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.H

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 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::kineticTheoryModels::frictionalStressModels::JohnsonJacksonSchaeffer
+
+Description
+
+SourceFiles
+    JohnsonJacksonSchaefferFrictionalStress.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef JohnsonJacksonSchaeffer_H
+#define JohnsonJacksonSchaeffer_H
+
+#include "frictionalStressModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace kineticTheoryModels
+{
+namespace frictionalStressModels
+{
+
+/*---------------------------------------------------------------------------*\
+               Class JohnsonJacksonSchaeffer Declaration
+\*---------------------------------------------------------------------------*/
+
+class JohnsonJacksonSchaeffer
+:
+    public frictionalStressModel
+{
+    // Private data
+
+        dictionary coeffDict_;
+
+        //- Material constant for frictional normal stress
+        dimensionedScalar Fr_;
+
+        //- Material constant for frictional normal stress
+        dimensionedScalar eta_;
+
+        //- Material constant for frictional normal stress
+        dimensionedScalar p_;
+
+        //- Angle of internal friction
+        dimensionedScalar phi_;
+
+        //- Lower limit for (alphaMax - alpha1)
+        dimensionedScalar alphaDeltaMin_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("JohnsonJacksonSchaeffer");
+
+
+    // Constructors
+
+        //- Construct from components
+        JohnsonJacksonSchaeffer(const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~JohnsonJacksonSchaeffer();
+
+
+    // Member functions
+
+        virtual tmp<volScalarField> frictionalPressure
+        (
+            const phaseModel& phase,
+            const dimensionedScalar& alphaMinFriction,
+            const dimensionedScalar& alphaMax
+        ) const;
+
+        virtual tmp<volScalarField> frictionalPressurePrime
+        (
+            const phaseModel& phase,
+            const dimensionedScalar& alphaMinFriction,
+            const dimensionedScalar& alphaMax
+        ) const;
+
+        virtual tmp<volScalarField> nu
+        (
+            const phaseModel& phase,
+            const dimensionedScalar& alphaMinFriction,
+            const dimensionedScalar& alphaMax,
+            const volScalarField& pf,
+            const volSymmTensorField& D
+        ) const;
+
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace frictionalStressModels
+} // End namespace kineticTheoryModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C

@@ -74,14 +74,16 @@ Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::
 frictionalPressure
 (
-    const volScalarField& alpha1,
+    const phaseModel& phase,
     const dimensionedScalar& alphaMinFriction,
     const dimensionedScalar& alphaMax
 ) const
 {
+    const volScalarField& alpha = phase;
+
     return
         dimensionedScalar("1e24", dimensionSet(1, -1, -2, 0, 0), 1e24)
-       *pow(Foam::max(alpha1 - alphaMinFriction, scalar(0)), 10.0);
+       *pow(Foam::max(alpha - alphaMinFriction, scalar(0)), 10.0);
 }
 
 
@@ -89,31 +91,31 @@ Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::
 frictionalPressurePrime
 (
-    const volScalarField& alpha1,
+    const phaseModel& phase,
     const dimensionedScalar& alphaMinFriction,
     const dimensionedScalar& alphaMax
 ) const
 {
+    const volScalarField& alpha = phase;
+
     return
         dimensionedScalar("1e25", dimensionSet(1, -1, -2, 0, 0), 1e25)
-       *pow(Foam::max(alpha1 - alphaMinFriction, scalar(0)), 9.0);
+       *pow(Foam::max(alpha - alphaMinFriction, scalar(0)), 9.0);
 }
 
 
 Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::nu
 (
-    const volScalarField& alpha1,
+    const phaseModel& phase,
     const dimensionedScalar& alphaMinFriction,
     const dimensionedScalar& alphaMax,
     const volScalarField& pf,
     const volSymmTensorField& D
 ) const
 {
-    const scalar I2Dsmall = 1.0e-15;
+    const volScalarField& alpha = phase;
 
-    // Creating nu assuming it should be 0 on the boundary which may not be
-    // true
     tmp<volScalarField> tnu
     (
         new volScalarField
@@ -121,13 +123,13 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::nu
             IOobject
             (
                 "Schaeffer:nu",
-                alpha1.mesh().time().timeName(),
+                phase.mesh().time().timeName(),
-                alpha1.mesh(),
+                phase.mesh(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
-            alpha1.mesh(),
+            phase.mesh(),
             dimensionedScalar("nu", dimensionSet(0, 2, -1, 0, 0), 0.0)
         )
     );
@@ -136,16 +138,33 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::nu
 
     forAll(D, celli)
     {
-        if (alpha1[celli] > alphaMinFriction.value())
+        if (alpha[celli] > alphaMinFriction.value())
         {
             nuf[celli] =
                 0.5*pf[celli]*sin(phi_.value())
                /(
-                    sqrt(1.0/6.0*(sqr(D[celli].xx() - D[celli].yy())
+                    sqrt((1.0/3.0)*sqr(tr(D[celli])) - invariantII(D[celli]))
-                  + sqr(D[celli].yy() - D[celli].zz())
+                  + SMALL
-                  + sqr(D[celli].zz() - D[celli].xx()))
+                );
-                  + sqr(D[celli].xy()) + sqr(D[celli].xz())
+        }
-                  + sqr(D[celli].yz())) + I2Dsmall
+    }
+
+    const fvPatchList& patches = phase.mesh().boundary();
+    const volVectorField& U = phase.U();
+
+    volScalarField::Boundary& nufBf = nuf.boundaryFieldRef();
+
+    forAll(patches, patchi)
+    {
+        if (!patches[patchi].coupled())
+        {
+            nufBf[patchi] =
+                (
+                    pf.boundaryField()[patchi]*sin(phi_.value())
+                   /(
+                        mag(U.boundaryField()[patchi].snGrad())
+                      + SMALL
+                    )
                 );
         }
     }

applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,21 +81,21 @@ public:
 
         virtual tmp<volScalarField> frictionalPressure
         (
-            const volScalarField& alpha1,
+            const phaseModel& phase,
             const dimensionedScalar& alphaMinFriction,
             const dimensionedScalar& alphaMax
         ) const;
 
         virtual tmp<volScalarField> frictionalPressurePrime
         (
-            const volScalarField& alpha1,
+            const phaseModel& phase,
             const dimensionedScalar& alphaMinFriction,
             const dimensionedScalar& alphaMax
         ) const;
 
         virtual tmp<volScalarField> nu
         (
-            const volScalarField& alpha1,
+            const phaseModel& phase,
             const dimensionedScalar& alphaMinFriction,
             const dimensionedScalar& alphaMax,
             const volScalarField& pf,

applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/frictionalStressModel/frictionalStressModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,6 +36,7 @@ SourceFiles
 #include "volFields.H"
 #include "dimensionedTypes.H"
 #include "runTimeSelectionTables.H"
+#include "phaseModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -107,21 +108,21 @@ public:
 
         virtual tmp<volScalarField> frictionalPressure
         (
-            const volScalarField& alpha1,
+            const phaseModel& phase,
             const dimensionedScalar& alphaMinFriction,
             const dimensionedScalar& alphaMax
         ) const = 0;
 
         virtual tmp<volScalarField> frictionalPressurePrime
         (
-            const volScalarField& alpha1f,
+            const phaseModel& phase,
             const dimensionedScalar& alphaMinFriction,
             const dimensionedScalar& alphaMax
         ) const = 0;
 
         virtual tmp<volScalarField> nu
         (
-            const volScalarField& alpha1,
+            const phaseModel& phase,
             const dimensionedScalar& alphaMinFriction,
             const dimensionedScalar& alphaMax,
             const volScalarField& pf,

applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C

@@ -26,6 +26,7 @@ License
 #include "kineticTheoryModel.H"
 #include "mathematicalConstants.H"
 #include "twoPhaseSystem.H"
+#include "fvOptions.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -43,10 +44,7 @@ Foam::RASModels::kineticTheoryModel::kineticTheoryModel
 :
     eddyViscosity
     <
-        RASModel<EddyDiffusivity<ThermalDiffusivity
+        RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>>
-        <
-            PhaseCompressibleTurbulenceModel<phaseModel>
-        >>>
     >
     (
         type,
@@ -113,6 +111,13 @@ Foam::RASModels::kineticTheoryModel::kineticTheoryModel
         coeffDict_
     ),
 
+    maxNut_
+    (
+        "maxNut",
+        dimensionSet(0,2,-1,0,0),
+        coeffDict_.lookupOrDefault<scalar>("maxNut",1000)
+    ),
+
     Theta_
     (
         IOobject
@@ -166,6 +171,20 @@ Foam::RASModels::kineticTheoryModel::kineticTheoryModel
         ),
         U.mesh(),
         dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0)
+    ),
+
+    nuFric_
+    (
+        IOobject
+        (
+            IOobject::groupName("nuFric", phase.name()),
+            U.time().timeName(),
+            U.mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        U.mesh(),
+        dimensionedScalar("zero", dimensionSet(0, 2, -1, 0, 0), 0.0)
     )
 {
     if (type == typeName)
@@ -189,10 +208,7 @@ bool Foam::RASModels::kineticTheoryModel::read()
     (
         eddyViscosity
         <
-            RASModel<EddyDiffusivity<ThermalDiffusivity
+            RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>>
-            <
-                PhaseCompressibleTurbulenceModel<phaseModel>
-            >>>
         >::read()
     )
     {
@@ -272,7 +288,7 @@ Foam::RASModels::kineticTheoryModel::pPrime() const
         )
      +  frictionalStressModel_->frictionalPressurePrime
         (
-            alpha_,
+            phase_,
             alphaMinFriction_,
             alphaMax_
         )
@@ -424,6 +440,8 @@ void Foam::RASModels::kineticTheoryModel::correct()
         // 'thermal' conductivity (Table 3.3, p. 49)
         kappa_ = conductivityModel_->kappa(alpha, Theta_, gs0_, rho, da, e_);
 
+        fv::options& fvOptions(fv::options::New(mesh_));
+
         // Construct the granular temperature equation (Eq. 3.20, p. 44)
         // NB. note that there are two typos in Eq. 3.20:
         //     Ps should be without grad
@@ -443,10 +461,13 @@ void Foam::RASModels::kineticTheoryModel::correct()
           + fvm::Sp(-gammaCoeff, Theta_)
           + fvm::Sp(-J1, Theta_)
           + fvm::Sp(J2/(Theta_ + ThetaSmall), Theta_)
+          + fvOptions(alpha, rho, Theta_)
         );
 
         ThetaEqn.relax();
+        fvOptions.constrain(ThetaEqn);
         ThetaEqn.solve();
+        fvOptions.correct(Theta_);
     }
     else
     {
@@ -519,24 +540,25 @@ void Foam::RASModels::kineticTheoryModel::correct()
         (
             frictionalStressModel_->frictionalPressure
             (
-                alpha,
+                phase_,
                 alphaMinFriction_,
                 alphaMax_
             )
         );
 
-        // Add frictional shear viscosity, Eq. 3.30, p. 52
+        nuFric_ = frictionalStressModel_->nu
-        nut_ += frictionalStressModel_->nu
         (
-            alpha,
+            phase_,
             alphaMinFriction_,
             alphaMax_,
             pf/rho,
             D
         );
 
-        // Limit viscosity
+        // Limit viscosity and add frictional viscosity
-        nut_.min(100);
+        nut_.min(maxNut_);
+        nuFric_ = min(nuFric_, maxNut_ - nut_);
+        nut_ += nuFric_;
     }
 
     if (debug)

applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H

@@ -47,8 +47,7 @@ SourceFiles
 
 #include "RASModel.H"
 #include "eddyViscosity.H"
-#include "PhaseCompressibleTurbulenceModel.H"
+#include "phaseCompressibleTurbulenceModel.H"
-#include "ThermalDiffusivity.H"
 #include "EddyDiffusivity.H"
 #include "phaseModel.H"
 #include "dragModel.H"
@@ -74,10 +73,7 @@ class kineticTheoryModel
 :
     public eddyViscosity
     <
-        RASModel<EddyDiffusivity<ThermalDiffusivity
+        RASModel<EddyDiffusivity<phaseCompressibleTurbulenceModel>>
-        <
-            PhaseCompressibleTurbulenceModel<phaseModel>
-        >>>
     >
 {
     // Private data
@@ -124,6 +120,9 @@ class kineticTheoryModel
             //- Residual phase fraction
             dimensionedScalar residualAlpha_;
 
+            //- Maximum turbulent viscosity
+            dimensionedScalar maxNut_;
+
 
         // Kinetic Theory Model Fields
 
@@ -139,6 +138,9 @@ class kineticTheoryModel
             //- The granular "thermal" conductivity
             volScalarField kappa_;
 
+            //- The frictional viscosity
+            volScalarField nuFric_;
+
 
     // Private Member Functions
 

applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phaseCompressibleTurbulenceModel.H

@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 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/>.
+
+Typedef
+    Foam::phaseCompressibleTurbulenceModel
+
+Description
+    Typedef for phaseCompressibleTurbulenceModel
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef phaseCompressibleTurbulenceModel_H
+#define phaseCompressibleTurbulenceModel_H
+
+#include "phaseCompressibleTurbulenceModelFwd.H"
+#include "PhaseCompressibleTurbulenceModel.H"
+#include "ThermalDiffusivity.H"
+#include "phaseModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+      typedef ThermalDiffusivity<PhaseCompressibleTurbulenceModel<phaseModel>>
+          phaseCompressibleTurbulenceModel;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phaseCompressibleTurbulenceModelFwd.H

@@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 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/>.
+
+Typedef
+    Foam::phaseCompressibleTurbulenceModel
+
+Description
+    Forward declaration of typedef for phaseCompressibleTurbulenceModel
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef phaseCompressibleTurbulenceModelFwd_H
+#define phaseCompressibleTurbulenceModelFwd_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    class phaseModel;
+
+    template<class TransportModel>
+    class PhaseCompressibleTurbulenceModel;
+
+    template<class BasicTurbulenceModel>
+    class ThermalDiffusivity;
+
+      typedef ThermalDiffusivity<PhaseCompressibleTurbulenceModel<phaseModel>>
+          phaseCompressibleTurbulenceModel;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //