diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
index 8950f3d958bbeb302c0c611f903e7713f7a286c7..ff1e0b0998e785a62760e51db2c1fe92d0d56fde 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,8 +24,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "ThermalPhaseChangePhaseSystem.H"
-#include "fvCFD.H"
 #include "alphatPhaseChangeWallFunctionFvPatchScalarField.H"
+#include "fvcVolumeIntegrate.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/Make/files b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/Make/files
index 300e1902a96327b6e6b6ef3316e0a0cf82f39f2d..fb9a8d1042cbf4fcc0050bd2a4733a8e2799290f 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/Make/files
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/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
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/derivedFvPatchFields/JohnsonJacksonParticleSlip/JohnsonJacksonParticleSlipFvPatchVectorField.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/derivedFvPatchFields/JohnsonJacksonParticleSlip/JohnsonJacksonParticleSlipFvPatchVectorField.C
index 4aa484cb5f93d1adc416c780f976f12b706b8f9f..04715ecbe749c493e02d8baba65d520927a21fa1 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/derivedFvPatchFields/JohnsonJacksonParticleSlip/JohnsonJacksonParticleSlipFvPatchVectorField.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/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());
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.C
index ab5fb531ea5daa3bf81bcc7ff283b4bdaa9c4a9f..ab6d3d354a1bc7e7e98d36e54b2cee3e6ac4b569 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/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_)
-       /pow(max(alphaMax - alpha1, alphaDeltaMin_), p_);
+        Fr_*pow(max(alpha - alphaMinFriction, scalar(0)), eta_)
+       /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)
-       *(alphaMax-alpha1)
-      + p_*pow(max(alpha1 - alphaMinFriction, scalar(0)), eta_)
-    )/pow(max(alphaMax - alpha1, alphaDeltaMin_), p_ + 1.0);
+        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::JohnsonJackson::nu
 (
-    const volScalarField& alpha1,
+    const phaseModel& phase,
     const dimensionedScalar& alphaMinFriction,
     const dimensionedScalar& alphaMax,
     const volScalarField& pf,
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.H
index 1fe92e7d0ea972fbaa7943ca7586d5ee30dba32d..3013fc25efe989cae381eef0f2a0a6bc5adca659 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/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,
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C
new file mode 100644
index 0000000000000000000000000000000000000000..4f05634d6f0ce37ff6aafc1356d4a13ff047f141
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C
@@ -0,0 +1,216 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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;
+    const scalar I2Dsmall = 1.0e-15;
+
+    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/6.0
+                       *(
+                            sqr(D[celli].xx() - D[celli].yy())
+                          + sqr(D[celli].yy() - D[celli].zz())
+                          + 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();
+
+    forAll(patches, patchi)
+    {
+        if (!patches[patchi].coupled())
+        {
+            nuf.boundaryField()[patchi] =
+                (
+                    pf.boundaryField()[patchi]*sin(phi_.value())
+                   /(
+                        mag(U.boundaryField()[patchi].snGrad())
+                      + I2Dsmall
+                    )
+                );
+        }
+    }
+
+    // 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;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.H
new file mode 100644
index 0000000000000000000000000000000000000000..2f8c971894e177d333066c0b62826154eec82e74
--- /dev/null
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/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
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
index b88ec0923abf187218bce9d510032cb9673c2961..8c072802d17df4132f48723f608d1ed8279bf223 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/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(),
-                alpha1.mesh(),
+                phase.mesh().time().timeName(),
+                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,31 @@ 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())
-                  + sqr(D[celli].yy() - D[celli].zz())
-                  + sqr(D[celli].zz() - D[celli].xx()))
-                  + sqr(D[celli].xy()) + sqr(D[celli].xz())
-                  + sqr(D[celli].yz())) + I2Dsmall
+                    sqrt((1.0/3.0)*sqr(tr(D[celli])) - invariantII(D[celli]))
+                  + SMALL
+                );
+        }
+    }
+
+    const fvPatchList& patches = phase.mesh().boundary();
+    const volVectorField& U = phase.U();
+
+    forAll(patches, patchi)
+    {
+        if (!patches[patchi].coupled())
+        {
+            nuf.boundaryField()[patchi] =
+                (
+                    pf.boundaryField()[patchi]*sin(phi_.value())
+                   /(
+                        mag(U.boundaryField()[patchi].snGrad())
+                      + SMALL
+                    )
                 );
         }
     }
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.H
index d4a304eb40c1a1444ae5c2ecdc428de44dcd5471..2b388df27a0d09556884cd3f417b8a8f18eead8a 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/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,
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/frictionalStressModel/frictionalStressModel.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/frictionalStressModel/frictionalStressModel.H
index 923a1fce3f71a273916c41946b48c3b3dff1b3b5..20feced23449f4517bdd0794889fba6077773166 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/frictionalStressModel/frictionalStressModel.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/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,
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
index f3b7d5da14091749b9961211d1b2820386a45dae..ff3be20ab10647ee6b6bc32793e374f06518b9e9 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
@@ -111,6 +111,13 @@ Foam::RASModels::kineticTheoryModel::kineticTheoryModel
         coeffDict_
     ),
 
+    maxNut_
+    (
+        "maxNut",
+        dimensionSet(0,2,-1,0,0),
+        coeffDict_.lookupOrDefault<scalar>("maxNut",1000)
+    ),
+
     Theta_
     (
         IOobject
@@ -164,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)
@@ -267,7 +288,7 @@ Foam::RASModels::kineticTheoryModel::pPrime() const
         )
      +  frictionalStressModel_->frictionalPressurePrime
         (
-            alpha_,
+            phase_,
             alphaMinFriction_,
             alphaMax_
         )
@@ -519,24 +540,25 @@ void Foam::RASModels::kineticTheoryModel::correct()
         (
             frictionalStressModel_->frictionalPressure
             (
-                alpha,
+                phase_,
                 alphaMinFriction_,
                 alphaMax_
             )
         );
 
-        // Add frictional shear viscosity, Eq. 3.30, p. 52
-        nut_ += frictionalStressModel_->nu
+        nuFric_ = frictionalStressModel_->nu
         (
-            alpha,
+            phase_,
             alphaMinFriction_,
             alphaMax_,
             pf/rho,
             D
         );
 
-        // Limit viscosity
-        nut_.min(100);
+        // Limit viscosity and add frictional viscosity
+        nut_.min(maxNut_);
+        nuFric_ = min(nuFric_, maxNut_ - nut_);
+        nut_ += nuFric_;
     }
 
     if (debug)
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H
index d648de7bc1104f1b216018dc855d84dade1a433a..03f323e519e6adbdd56657e03ef554ee7034f1f1 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H
@@ -120,6 +120,9 @@ class kineticTheoryModel
             //- Residual phase fraction
             dimensionedScalar residualAlpha_;
 
+            //- Maximum turbulent viscosity
+            dimensionedScalar maxNut_;
+
 
         // Kinetic Theory Model Fields
 
@@ -135,6 +138,9 @@ class kineticTheoryModel
             //- The granular "thermal" conductivity
             volScalarField kappa_;
 
+            //- The frictional viscosity
+            volScalarField nuFric_;
+
 
     // Private Member Functions
 
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/T.gas b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/T.gas
new file mode 100644
index 0000000000000000000000000000000000000000..b997ba8d87f6993e039af15822ff20ef1ae7fb74
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/T.gas
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T.gas;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 0 0 1 0 0 0];
+
+internalField       uniform 300;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 300;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.gas;
+        inletValue      uniform 300;
+        value           uniform 300;
+    }
+
+    "wall.*"
+    {
+        type            zeroGradient;
+    }
+
+    "frontAndBack.*"
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/T.solids b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/T.solids
new file mode 100644
index 0000000000000000000000000000000000000000..f89c81a228121fb7c9f8bf40cb654863af290e75
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/T.solids
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T.solids;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 0 0 1 0 0 0];
+
+internalField       uniform 300;
+
+boundaryField
+{
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.solids;
+        inletValue      uniform 300;
+        value           uniform 300;
+    }
+
+    "wall.*"
+    {
+        type            zeroGradient;
+    }
+
+    "frontAndBack.*"
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/Theta.solids b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/Theta.solids
new file mode 100644
index 0000000000000000000000000000000000000000..15e74863026c46ab8db16d6a4896ff6ee1a0118a
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/Theta.solids
@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      Theta.solids;
+}
+// ************************************************************************* //
+
+dimensions          [0 2 -2 0 0 0 0];
+
+internalField       uniform 0;
+
+referenceLevel      1e-4;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 1e-4;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    "wall.*"
+    {
+        type            JohnsonJacksonParticleTheta;
+        restitutionCoefficient 0.2;
+        specularityCoefficient   0.1;
+        muF             0.25;
+        sigma           2;
+        value           uniform 1e-4;
+    }
+
+    "frontAndBack.*"
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/U.gas b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/U.gas
new file mode 100644
index 0000000000000000000000000000000000000000..c2eba7f790c9c2056d4892e0b8a6e3d685862d37
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/U.gas
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U.gas;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    outlet
+    {
+        type            pressureInletOutletVelocity;
+        phi             phi.gas;
+        value           $internalField;
+    }
+
+    "wall.*"
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    "frontAndBack.*"
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/U.solids b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/U.solids
new file mode 100644
index 0000000000000000000000000000000000000000..83efb2917f8162c1447bc7e6cb18dd931c5bfe66
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/U.solids
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U.solids;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform (0 -0.2 0);
+    }
+
+    outlet
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    "wall.*"
+    {
+        type            JohnsonJacksonParticleSlip;
+        restitutionCoefficient 0.2;
+        specularityCoefficient   0.1;
+        muF             0.25;
+        sigma           2;
+        value           uniform (0 0 0);
+    }
+
+    "frontAndBack.*"
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/alpha.solids b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/alpha.solids
new file mode 100644
index 0000000000000000000000000000000000000000..7b10e13ea304e9742e3812e2cc472b45ce6106ef
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/alpha.solids
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alpha.solids;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.3;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    "wall.*"
+    {
+        type            zeroGradient;
+    }
+
+    "frontAndBack.*"
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/alphat.gas b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/alphat.gas
new file mode 100644
index 0000000000000000000000000000000000000000..a427a3d22089f842b607efd6b3a95e9efb8e1571
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/alphat.gas
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alphat.gas;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    "wall.*"
+    {
+        type            compressible::alphatWallFunction;
+        Prt             0.85;
+        value           $internalField;
+    }
+
+    "frontAndBack.*"
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/alphat.solids b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/alphat.solids
new file mode 100644
index 0000000000000000000000000000000000000000..629d9b81b2f7592eba1aede232120a7e2a938493
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/alphat.solids
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alphat.solids;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    "wall.*"
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    "frontAndBack.*"
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/epsilon.gas b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/epsilon.gas
new file mode 100644
index 0000000000000000000000000000000000000000..d71220e918cbb5cdc0e6524ec2f0fe5d4128f242
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/epsilon.gas
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      epsilon.gas;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 2 -3 0 0 0 0];
+
+internalField       uniform 1e-3;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.gas;
+        inletValue      $internalField;
+        value           $internalField;
+    }
+
+    "wall.*"
+    {
+        type            epsilonWallFunction;
+        value           $internalField;
+    }
+
+    "frontAndBack.*"
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/k.gas b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/k.gas
new file mode 100644
index 0000000000000000000000000000000000000000..7b712450dbfbc4661e4cafad3dfb1db514b467de
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/k.gas
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      k.gas;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 2 -2 0 0 0 0];
+
+internalField       uniform 1e-4;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.gas;
+        inletValue      $internalField;
+        value           $internalField;
+    }
+
+    "wall.*"
+    {
+        type            kqRWallFunction;
+        value           $internalField;
+    }
+
+    "frontAndBack.*"
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/nut.gas b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/nut.gas
new file mode 100644
index 0000000000000000000000000000000000000000..48bbb42864d6a9634bb033d1f6d11df1b00b781f
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/nut.gas
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      nut.gas;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 2 -1 0 0 0 0];
+
+internalField       uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    "wall.*"
+    {
+        type            nutkWallFunction;
+        value           $internalField;
+    }
+
+    "frontAndBack.*"
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/nut.solids b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/nut.solids
new file mode 100644
index 0000000000000000000000000000000000000000..978d7897046179ac59325dc1f87c8121a103c976
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/nut.solids
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      nut.solids;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 2 -1 0 0 0 0];
+
+internalField       uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    "wall.*"
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    "frontAndBack.*"
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/p b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/p
new file mode 100644
index 0000000000000000000000000000000000000000..2ab8e8cdcccf07a4eb0de680c83a7dbfd681efbb
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/p
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [1 -1 -2 0 0 0 0];
+
+internalField       uniform 1e5;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    "wall.*"
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    "frontAndBack.*"
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/p_rgh b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..0b84a879587303f2279405666bf5d0d61e4ff843
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/0/p_rgh
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [1 -1 -2 0 0 0 0];
+
+internalField       uniform 1e5;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedFluxPressure;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            prghTotalPressure;
+        p0              $internalField;
+        U               U.gas;
+        phi             phi.gas;
+        rho             thermo:rho.gas;
+        value           $internalField;
+    }
+
+    "wall.*"
+    {
+        type            fixedFluxPressure;
+        value           $internalField;
+    }
+
+    "frontAndBack.*"
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/g b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/g
new file mode 100644
index 0000000000000000000000000000000000000000..e0ac2653b5b370ad62f6770588121d30cac51627
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           ( 0 -9.81 0 );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/phaseProperties b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/phaseProperties
new file mode 100644
index 0000000000000000000000000000000000000000..d53b4d7196fad22aec495f92d8c5b7ab70920b79
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/phaseProperties
@@ -0,0 +1,119 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      phaseProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+type    heatAndMomentumTransferTwoPhaseSystem;
+
+phases (solids gas);
+
+solids
+{
+    type          purePhaseModel;
+
+    diameterModel constant;
+
+    constantCoeffs
+    {
+        d               462e-6;
+    }
+
+    residualAlpha   1e-5;
+}
+
+gas
+{
+    type          purePhaseModel;
+
+    diameterModel constant;
+    constantCoeffs
+    {
+        d               1;
+    }
+    residualAlpha   1e-5;
+}
+
+blending
+{
+    default
+    {
+        type            none;
+        residualAlpha   1e-6;
+        continuousPhase gas;
+    }
+}
+
+surfaceTension
+(
+    (gas and solids)
+    {
+        type            constant;
+        sigma           0;
+    }
+);
+
+aspectRatio
+(
+);
+
+drag
+(
+    (solids in gas)
+    {
+        type            GidaspowErgunWenYu;
+        residualAlpha   1e-5;
+        residualRe      1e-5;
+        swarmCorrection
+        {
+            type        none;
+        }
+    }
+);
+
+virtualMass
+(
+    (solids in gas)
+    {
+        type            constantCoefficient;
+        Cvm             0;
+    }
+);
+
+heatTransfer
+(
+    (solids in gas)
+    {
+        type            RanzMarshall;
+        residualAlpha   1e-4;
+    }
+);
+
+lift
+(
+);
+
+wallLubrication
+(
+);
+
+turbulentDispersion
+(
+);
+
+// Minimum allowable pressure
+pMin            10000;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/thermophysicalProperties.gas b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/thermophysicalProperties.gas
new file mode 100644
index 0000000000000000000000000000000000000000..a2a98f45f2c4e162d59c784de1b5e58f8a22348b
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/thermophysicalProperties.gas
@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties.gas;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectGas;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   28.9;
+    }
+    thermodynamics
+    {
+        Cp          1007;
+        Hf          0;
+    }
+    equationOfState
+    {
+        rho         1.2;
+    }
+    transport
+    {
+        mu          1.84e-05;
+        Pr          0.7;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/thermophysicalProperties.solids b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/thermophysicalProperties.solids
new file mode 100644
index 0000000000000000000000000000000000000000..84dae4b7166c241435c8fc55987ee9a24ce067d8
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/thermophysicalProperties.solids
@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties.solids;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState rhoConst;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   100;
+    }
+    equationOfState
+    {
+        rho         2480;
+    }
+    thermodynamics
+    {
+        Cp          6000;
+        Hf          0;
+    }
+    transport
+    {
+        mu          0;
+        Pr          24.47;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/turbulenceProperties.gas b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/turbulenceProperties.gas
new file mode 100644
index 0000000000000000000000000000000000000000..d92fbf9947c21aa42bebab8def094f82ace2005e
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/turbulenceProperties.gas
@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties.gas;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  RAS;
+
+RAS
+{
+    RASModel kEpsilon;
+
+    turbulence      on;
+    printCoeffs     on;
+}
+
+LES
+{
+    LESModel Smagorinsky;
+
+    turbulence      on;
+    printCoeffs     on;
+
+    delta cubeRootVol;
+
+    cubeRootVolCoeffs
+    {
+        deltaCoeff 1;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/turbulenceProperties.solids b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/turbulenceProperties.solids
new file mode 100644
index 0000000000000000000000000000000000000000..04793a45c21c537f1d65870f4372fea7a7c2224d
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/constant/turbulenceProperties.solids
@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties.solids;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  RAS;
+
+RAS
+{
+    RASModel kineticTheory;
+
+    turbulence      on;
+    printCoeffs     on;
+
+    kineticTheoryCoeffs
+    {
+        equilibrium             off;
+
+        e                       0.8;
+        alphaMax                0.65;
+        alphaMinFriction        0.5;
+        residualAlpha           1e-6;
+
+        viscosityModel          Syamlal;
+        conductivityModel       Syamlal;
+        granularPressureModel   SyamlalRogersOBrien;
+        frictionalStressModel   Schaeffer;
+        radialModel             CarnahanStarling;
+
+        SchaefferCoeffs
+        {
+            phi                     36;
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/system/blockMeshDict b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/system/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..ebf2ada3b0145848997436388a6b41acd3c4ac99
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/system/blockMeshDict
@@ -0,0 +1,128 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+vertices
+(
+    ( 0.0  0.0  -0.01)
+    ( 0.04 0.0  -0.01)
+    ( 0.04 0.15  -0.01)
+    ( 0.0  0.15  -0.01)
+    ( 0.0  0.0   0.01)
+    ( 0.04 0.0   0.01)
+    ( 0.04 0.15   0.01)
+    ( 0.0  0.15   0.01)
+
+    ( 0.04  0.0  -0.01)
+    ( 0.5 0.0  -0.01)
+    ( 0.5 0.04  -0.01)
+    ( 0.04  0.04  -0.01)
+    ( 0.04  0.0   0.01)
+    ( 0.5 0.0   0.01)
+    ( 0.5 0.04   0.01)
+    ( 0.04  0.04   0.01)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (16 60 1) simpleGrading (1 1 1)
+    hex (8 9 10 11 12 13 14 15) (180 16 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+    wall1
+    {
+        type wall;
+        faces
+        (
+            (0 4 7 3)
+            (1 5 4 0)
+        );
+    }
+    wall_merge
+    {
+        type wall;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+    frontAndBack1
+    {
+        type empty;
+        faces
+        (
+            (0 3 2 1)
+            (4 5 6 7)
+        );
+    }
+    wall2
+    {
+        type wall;
+        faces
+        (
+            (11 15 14 10)
+            (9 13 12 8)
+        );
+    }
+    wall_merge2
+    {
+        type wall;
+        faces
+        (
+            (8 12 15 11)
+        );
+    }
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (10 14 13 9)
+        );
+    }
+    frontAndBack2
+    {
+        type empty;
+        faces
+        (
+            (8 11 10 9)
+            (12 13 14 15)
+        );
+    }
+
+);
+
+mergePatchPairs
+(
+    (wall_merge wall_merge2)
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/system/controlDict b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..f929ca081e9e0ab9142f219e8a754064247647eb
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/system/controlDict
@@ -0,0 +1,55 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     reactingTwoPhaseEulerFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         2;
+
+deltaT          1e-4;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  9;
+
+writeCompression compressed;
+
+timeFormat      general;
+
+timePrecision   8;
+
+runTimeModifiable on;
+
+adjustTimeStep  yes;
+
+maxCo           0.1;
+
+maxDeltaT       0.01;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/system/fvSchemes b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..293328da6c694edaeb03b00089f1f712180b7fc6
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/system/fvSchemes
@@ -0,0 +1,74 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default     Euler;
+}
+
+gradSchemes
+{
+    default     Gauss linear;
+}
+
+divSchemes
+{
+    default                         none;
+
+    "div\(phi,alpha.*\)"            Gauss vanLeer;
+    "div\(phir,alpha.*\)"           Gauss vanLeer;
+
+    "div\(alphaRhoPhi.*,U.*\)"      Gauss limitedLinearV 1;
+    "div\(phi.*,U.*\)"              Gauss limitedLinearV 1;
+
+    "div\(alphaRhoPhi.*,(h|e).*\)"  Gauss limitedLinear 1;
+    "div\(alphaRhoPhi.*,K.*\)"      Gauss limitedLinear 1;
+    "div\(alphaPhi.*,p\)"           Gauss limitedLinear 1;
+
+    div(alphaRhoPhi.solids,Theta.solids) Gauss limitedLinear 1;
+
+    "div\(alphaRhoPhi.*,(k|epsilon).*\)"  Gauss limitedLinear 1;
+
+    div((((alpha.gas*thermo:rho.gas)*nuEff.gas)*dev2(T(grad(U.gas))))) Gauss linear;
+
+     div((((thermo:rho.solids*nut.solids)*dev2(T(grad(U.solids))))+(((thermo:rho.solids*lambda.solids)*div(phi.solids))*I)))  Gauss linear;
+}
+
+laplacianSchemes
+{
+    default     Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default     linear;
+}
+
+snGradSchemes
+{
+    default    corrected;
+}
+
+fluxRequired
+{
+    default     no;
+    p_rgh;
+    alpha.solids;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/system/fvSolution b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..10ff75d72cefbb285a9cc7cc97e12e1c24ba2b80
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/LBend/system/fvSolution
@@ -0,0 +1,111 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    "alpha.*"
+    {
+        nAlphaCorr      2;
+        nAlphaSubCycles 1;
+
+        implicitPhasePressure yes;
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-10;
+        relTol          0;
+        minIter         1;
+    }
+
+    p_rgh
+    {
+        solver          GAMG;
+        smoother        DIC;
+        nPreSweeps      0;
+        nPostSweeps     2;
+        nFinestSweeps   2;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        processorAgglomerator procFaces;
+        nAgglomeratingCells 800;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+        tolerance       1e-8;
+        relTol          0.01;
+    }
+
+    p_rghFinal
+    {
+        $p_rgh;
+        relTol          0;
+    }
+
+    "U.*"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-7;
+        relTol          0;
+        minIter         1;
+    }
+
+    "(h|e).*"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-8;
+        relTol          0;
+        minIter         1;
+        maxIter         10;
+    }
+
+    "Theta.*"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-7;
+        relTol          0;
+        minIter         1;
+    }
+
+    "(k|epsilon).*"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-5;
+        relTol          0;
+        minIter         1;
+    }
+}
+
+PIMPLE
+{
+    nOuterCorrectors 3;
+    nCorrectors      1;
+    nNonOrthogonalCorrectors 0;
+    faceMomentum     yes;
+}
+
+relaxationFactors
+{
+    equations
+    {
+        ".*"            1;
+    }
+}
+
+
+// ************************************************************************* //