reactingTwoPhaseEulerFoam: Generalize the handling of the dilatation rate

to support any number of phases

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H

@@ -76,7 +76,7 @@ tmp<surfaceScalarField> phiF2;
       + (fvc::interpolate(rAU1*F) & mesh.Sf())
     );
 
-    // Phase-1 dispersion, lift and wall-lubrication flux
+    // Phase-2 dispersion, lift and wall-lubrication flux
     phiF2 =
     (
       - Df2*snGradAlpha1
@@ -339,11 +339,14 @@ while (pimple.correct())
             }
 
             // Compressibility correction for phase-fraction equations
-            fluid.dgdt() =
-            (
-                alpha1*(pEqnComp2 & p_rgh)
-              - alpha2*(pEqnComp1 & p_rgh)
-            );
+            if (phase1.compressible())
+            {
+                phase1.D(pEqnComp1 & p_rgh);
+            }
+            if (phase2.compressible())
+            {
+                phase2.D(pEqnComp2 & p_rgh);
+            }
 
             // Optionally relax pressure for velocity correction
             p_rgh.relax();

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H

@@ -323,11 +323,15 @@ while (pimple.correct())
             U2.correctBoundaryConditions();
             fvOptions.correct(U2);
 
-            fluid.dgdt() =
-            (
-                alpha1*(pEqnComp2 & p_rgh)
-              - alpha2*(pEqnComp1 & p_rgh)
-            );
+            // Compressibility correction for phase-fraction equations
+            if (phase1.compressible())
+            {
+                phase1.D(pEqnComp1 & p_rgh);
+            }
+            if (phase2.compressible())
+            {
+                phase2.D(pEqnComp2 & p_rgh);
+            }
         }
     }
 

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C

@@ -36,6 +36,17 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::AnisothermalPhaseModel
 )
 :
     BasePhaseModel(fluid, phaseName),
+    D_
+    (
+        IOobject
+        (
+            IOobject::groupName("D", this->name()),
+            fluid.mesh().time().timeName(),
+            fluid.mesh()
+        ),
+        fluid.mesh(),
+        dimensionedScalar("D", dimless/dimTime, 0)
+    ),
     K_
     (
         IOobject
@@ -120,4 +131,27 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
 }
 
 
+template<class BasePhaseModel>
+bool Foam::AnisothermalPhaseModel<BasePhaseModel>::compressible() const
+{
+    return true;
+}
+
+
+template<class BasePhaseModel>
+Foam::tmp<Foam::volScalarField>
+Foam::AnisothermalPhaseModel<BasePhaseModel>::D() const
+{
+    return D_;
+}
+
+
+template<class BasePhaseModel>
+void
+Foam::AnisothermalPhaseModel<BasePhaseModel>::D(const volScalarField& D)
+{
+    D_ = D;
+}
+
+
 // ************************************************************************* //

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H

@@ -54,6 +54,9 @@ class AnisothermalPhaseModel
 {
     // Private data
 
+        //- Dilatation
+        volScalarField D_;
+
         //- Kinetic energy
         volScalarField K_;
 
@@ -79,6 +82,18 @@ public:
 
         //- Return the enthalpy equation
         virtual tmp<fvScalarMatrix> heEqn();
+
+
+        // Compressibility (variable density)
+
+            //- Return true if the phase is compressible otherwise false
+            virtual bool compressible() const;
+
+            //- Phase dilatation rate ((alpha/rho)*Drho/Dt)
+            virtual tmp<volScalarField> D() const;
+
+            //- Set phase dilatation rate ((alpha/rho)*Drho/Dt)
+            virtual void D(const volScalarField& D);
 };
 
 

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C

@@ -134,4 +134,24 @@ bool Foam::phaseModel::read()
 }
 
 
+bool Foam::phaseModel::compressible() const
+{
+    return false;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::phaseModel::D() const
+{
+    return tmp<volScalarField>();
+}
+
+
+void Foam::phaseModel::D(const volScalarField& D)
+{
+    WarningIn("phaseModel::D(const volScalarField& D)")
+        << "Attempt to set the dilatation rate of an incompressible phase"
+        << endl;
+}
+
+
 // ************************************************************************* //

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H

@@ -83,6 +83,7 @@ public:
     //- Runtime type information
     ClassName("phaseModel");
 
+
     // Declare runtime construction
 
         declareRunTimeSelectionTable
@@ -163,6 +164,18 @@ public:
         virtual bool read();
 
 
+        // Compressibility (variable density)
+
+            //- Return true if the phase is compressible otherwise false
+            virtual bool compressible() const;
+
+            //- Phase dilatation rate ((alpha/rho)*Drho/Dt)
+            virtual tmp<volScalarField> D() const;
+
+            //- Set phase dilatation rate ((alpha/rho)*Drho/Dt)
+            virtual void D(const volScalarField& D);
+
+
         // Thermo
 
             //- Return const access to the thermophysical model

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C

@@ -183,20 +183,6 @@ Foam::phaseSystem::phaseSystem
         ),
         mesh,
         dimensionedScalar("dpdt", dimPressure/dimTime, 0)
-    ),
-
-    dgdt_
-    (
-        IOobject
-        (
-            "dgdt",
-            mesh.time().timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("dgdt", dimless/dimTime, 0)
     )
 {
     // Blending methods

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H

@@ -174,9 +174,6 @@ protected:
         //- Rate of change of pressure
         volScalarField dpdt_;
 
-        //- Dilatation
-        volScalarField dgdt_;
-
         //- Blending methods
         blendingMethodTable blendingMethods_;
 
@@ -366,12 +363,6 @@ public:
             //- Access the rate of change of the pressure
             inline volScalarField& dpdt();
 
-            //- Constant access the dilatation parameter
-            inline const volScalarField& dgdt() const;
-
-            //- Access the dilatation parameter
-            inline volScalarField& dgdt();
-
             //- Access a sub model between a phase pair
             template <class modelType>
             const modelType& lookupSubModel(const phasePair& key) const;

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H

@@ -54,17 +54,5 @@ inline Foam::volScalarField& Foam::phaseSystem::dpdt()
     return dpdt_;
 }
 
-inline const Foam::volScalarField& Foam::phaseSystem::dgdt() const
-{
-    return dgdt_;
-}
-
-
-inline Foam::volScalarField& Foam::phaseSystem::dgdt()
-{
-    return dgdt_;
-}
-
-
 
 // ************************************************************************* //

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C

@@ -170,12 +170,6 @@ void Foam::twoPhaseSystem::solve()
     volScalarField& alpha1 = phase1_;
     volScalarField& alpha2 = phase2_;
 
-    const surfaceScalarField& phi = this->phi();
-    const surfaceScalarField& phi1 = phase1_.phi();
-    const surfaceScalarField& phi2 = phase2_.phi();
-
-    const volScalarField& dgdt = this->dgdt();
-
     const dictionary& alphaControls = mesh.solverDict(alpha1.name());
 
     label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
@@ -184,22 +178,59 @@ void Foam::twoPhaseSystem::solve()
     word alphaScheme("div(phi," + alpha1.name() + ')');
     word alpharScheme("div(phir," + alpha1.name() + ')');
 
-    alpha1.correctBoundaryConditions();
+    const surfaceScalarField& phi = this->phi();
+    const surfaceScalarField& phi1 = phase1_.phi();
+    const surfaceScalarField& phi2 = phase2_.phi();
 
+    // Construct the dilatation rate source term
+    tmp<volScalarField::DimensionedInternalField> tdgdt;
+
+    if (phase1_.compressible() && phase2_.compressible())
+    {
+        tdgdt =
+        (
+            alpha1.dimensionedInternalField()
+           *phase2_.D()().dimensionedInternalField()
+          - alpha2.dimensionedInternalField()
+           *phase1_.D()().dimensionedInternalField()
+        );
+    }
+    else if (phase1_.compressible())
+    {
+        tdgdt =
+        (
+          - alpha2.dimensionedInternalField()
+           *phase1_.D()().dimensionedInternalField()
+        );
+    }
+    else if (phase2_.compressible())
+    {
+        tdgdt =
+        (
+            alpha1.dimensionedInternalField()
+           *phase2_.D()().dimensionedInternalField()
+        );
+    }
+
+    alpha1.correctBoundaryConditions();
+    surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0))));
 
     surfaceScalarField phic("phic", phi);
     surfaceScalarField phir("phir", phi1 - phi2);
 
-    surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0))));
+    tmp<surfaceScalarField> alpha1alpha2f;
 
     if (pPrimeByA_.valid())
     {
+        alpha1alpha2f =
+            fvc::interpolate(max(alpha1, scalar(0)))
+           *fvc::interpolate(max(alpha2, scalar(0)));
+
         surfaceScalarField phiP
         (
             pPrimeByA_()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf()
         );
 
-        phic += alpha1f*phiP;
         phir += phiP;
     }
 
@@ -214,7 +245,7 @@ void Foam::twoPhaseSystem::solve()
                 mesh
             ),
             mesh,
-            dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
+            dimensionedScalar("Sp", dimless/dimTime, 0.0)
         );
 
         volScalarField::DimensionedInternalField Su
@@ -230,6 +261,10 @@ void Foam::twoPhaseSystem::solve()
             fvc::div(phi)*min(alpha1, scalar(1))
         );
 
+        if (tdgdt.valid())
+        {
+            scalarField& dgdt = tdgdt();
+
             forAll(dgdt, celli)
             {
                 if (dgdt[celli] > 0.0)
@@ -242,6 +277,7 @@ void Foam::twoPhaseSystem::solve()
                     Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
                 }
             }
+        }
 
         surfaceScalarField alphaPhic1
         (
@@ -336,7 +372,7 @@ void Foam::twoPhaseSystem::solve()
             fvScalarMatrix alpha1Eqn
             (
                 fvm::ddt(alpha1) - fvc::ddt(alpha1)
-              - fvm::laplacian(alpha1f*pPrimeByA_(), alpha1, "bounded")
+              - fvm::laplacian(alpha1alpha2f*pPrimeByA_(), alpha1, "bounded")
             );
 
             alpha1Eqn.relax();