reactingEulerFoam: Support compressibility and mass-transfer independently

Now combinations of incompressible, compressible phases with or without
mass-transfer are supported efficiently.

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C

@@ -151,6 +151,16 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+template<class BasePhaseSystem>
+bool Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::transfersMass
+(
+    const phaseModel& phase
+) const
+{
+    return true;
+}
+
+
 template<class BasePhaseSystem>
 Foam::tmp<Foam::volScalarField>
 Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::dmdt

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H

@@ -120,6 +120,9 @@ public:
 
     // Member Functions
 
+        //- Return true if there is mass transfer for phase
+        virtual bool transfersMass(const phaseModel& phase) const;
+
         //- Return the interfacial mass flow rate
         virtual tmp<volScalarField> dmdt(const phasePairKey& key) const;
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C

@@ -92,20 +92,7 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::dmdt
     const Foam::phaseModel& phase
 ) const
 {
-    return tmp<volScalarField>
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                IOobject::groupName("dmdt", phase.name()),
-                this->mesh().time().timeName(),
-                this->mesh().time()
-            ),
-            this->mesh(),
-            dimensionedScalar("zero", dimDensity/dimTime, 0)
-        )
-    );
+    return tmp<volScalarField>(NULL);
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C

@@ -25,6 +25,7 @@ License
 
 #include "ThermalPhaseChangePhaseSystem.H"
 #include "fvCFD.H"
+#include "alphatPhaseChangeWallFunctionFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -39,7 +40,61 @@ ThermalPhaseChangePhaseSystem
     volatile_(this->lookup("volatile")),
     saturationModel_(saturationModel::New(this->subDict("saturationModel"))),
     massTransfer_(this->lookup("massTransfer"))
-{}
+{
+
+    forAllConstIter
+    (
+        phaseSystem::phasePairTable,
+        this->phasePairs_,
+        phasePairIter
+    )
+    {
+        const phasePair& pair(phasePairIter());
+
+        if (pair.ordered())
+        {
+            continue;
+        }
+
+        // Initially assume no mass transfer
+        iDmdt_.insert
+        (
+            pair,
+            new volScalarField
+            (
+                IOobject
+                (
+                    IOobject::groupName("iDmdt", pair.name()),
+                    this->mesh().time().timeName(),
+                    this->mesh(),
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                this->mesh(),
+                dimensionedScalar("zero", dimDensity/dimTime, 0)
+            )
+        );
+
+        // Initially assume no mass transfer
+        wDmdt_.insert
+        (
+            pair,
+            new volScalarField
+            (
+                IOobject
+                (
+                    IOobject::groupName("wDmdt", pair.name()),
+                    this->mesh().time().timeName(),
+                    this->mesh(),
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                this->mesh(),
+                dimensionedScalar("zero", dimDensity/dimTime, 0)
+            )
+        );
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -52,6 +107,14 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+template<class BasePhaseSystem>
+const Foam::saturationModel&
+Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::saturation() const
+{
+    return saturationModel_();
+}
+
+
 template<class BasePhaseSystem>
 Foam::autoPtr<Foam::phaseSystem::massTransferTable>
 Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
@@ -121,6 +184,9 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
 template<class BasePhaseSystem>
 void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
 {
+    typedef compressible::alphatPhaseChangeWallFunctionFvPatchScalarField
+        alphatPhaseChangeWallFunction;
+
     BasePhaseSystem::correctThermo();
 
     forAllConstIter
@@ -147,6 +213,8 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
         const volScalarField& he2(phase2.thermo().he());
 
         volScalarField& dmdt(*this->dmdt_[pair]);
+        volScalarField& iDmdt(*this->iDmdt_[pair]);
+        volScalarField& wDmdt(*this->wDmdt_[pair]);
 
         volScalarField& Tf = *this->Tf_[pair];
 
@@ -167,25 +235,28 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
 
             Tf = saturationModel_->Tsat(phase1.thermo().p());
 
-            dmdt =
-                (H1*(Tf - T1) + H2*(Tf - T2))
+            scalar iDmdtRelax(this->mesh().fieldRelaxationFactor("iDmdt"));
+
+            iDmdt =
+                (1 - iDmdtRelax)*iDmdt
+              + iDmdtRelax*(H1*(Tf - T1) + H2*(Tf - T2))
                /min
                 (
-                    (pos(dmdt)*he2 + neg(dmdt)*hef2)
-                  - (neg(dmdt)*he1 + pos(dmdt)*hef1),
+                    (pos(iDmdt)*he2 + neg(iDmdt)*hef2)
+                  - (neg(iDmdt)*he1 + pos(iDmdt)*hef1),
                     0.3*mag(hef2 - hef1)
                 );
 
-            Info<< "dmdt." << pair.name()
-                << ": min = " << min(dmdt.internalField())
-                << ", mean = " << average(dmdt.internalField())
-                << ", max = " << max(dmdt.internalField())
-                << ", integral = " << fvc::domainIntegrate(dmdt).value()
+            Info<< "iDmdt." << pair.name()
+                << ": min = " << min(iDmdt.internalField())
+                << ", mean = " << average(iDmdt.internalField())
+                << ", max = " << max(iDmdt.internalField())
+                << ", integral = " << fvc::domainIntegrate(iDmdt).value()
                 << endl;
         }
         else
         {
-            dmdt == dimensionedScalar("0", dmdt.dimensions(), 0);
+            iDmdt == dimensionedScalar("0", dmdt.dimensions(), 0);
         }
 
         volScalarField H1(this->heatTransferModels_[pair][pair.first()]->K());
@@ -200,12 +271,13 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
 
         volScalarField mDotL
         (
-            dmdt*
+            iDmdt*
             (
-                (pos(dmdt)*he2 + neg(dmdt)*hef2)
-              - (neg(dmdt)*he1 + pos(dmdt)*hef1)
+                (pos(iDmdt)*he2 + neg(iDmdt)*hef2)
+              - (neg(iDmdt)*he1 + pos(iDmdt)*hef1)
             )
         );
+
         Tf = (H1*T1 + H2*T2 + mDotL)/(H1 + H2);
 
         Info<< "Tf." << pair.name()
@@ -213,6 +285,68 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
             << ", mean = " << average(Tf.internalField())
             << ", max = " << max(Tf.internalField())
             << endl;
+
+        // Accumulate dmdt contributions from boundaries
+        if
+        (
+            phase2.mesh().foundObject<volScalarField>
+            (
+                "alphat." +  phase2.name()
+            )
+        )
+        {
+            scalar wDmdtRelax(this->mesh().fieldRelaxationFactor("wDmdt"));
+            wDmdt *= (1 - wDmdtRelax);
+
+            const volScalarField& alphat =
+                phase2.mesh().lookupObject<volScalarField>
+                (
+                    "alphat." +  phase2.name()
+                );
+
+            const fvPatchList& patches = this->mesh().boundary();
+            forAll(patches, patchi)
+            {
+                const fvPatch& currPatch = patches[patchi];
+
+                if
+                (
+                    isA<alphatPhaseChangeWallFunction>
+                    (
+                        alphat.boundaryField()[patchi]
+                    )
+                )
+                {
+                    const scalarField& patchDmdt =
+                        refCast<const alphatPhaseChangeWallFunction>
+                        (
+                            alphat.boundaryField()[patchi]
+                        ).dmdt();
+
+                    forAll(patchDmdt,facei)
+                    {
+                        label faceCelli = currPatch.faceCells()[facei];
+                        wDmdt[faceCelli] += wDmdtRelax*patchDmdt[facei];
+                    }
+                }
+            }
+
+            Info<< "wDmdt." << pair.name()
+                << ": min = " << min(wDmdt.internalField())
+                << ", mean = " << average(wDmdt.internalField())
+                << ", max = " << max(wDmdt.internalField())
+                << ", integral = " << fvc::domainIntegrate(wDmdt).value()
+                << endl;
+        }
+
+        dmdt = wDmdt + iDmdt;
+
+        Info<< "dmdt." << pair.name()
+            << ": min = " << min(dmdt.internalField())
+            << ", mean = " << average(dmdt.internalField())
+            << ", max = " << max(dmdt.internalField())
+            << ", integral = " << fvc::domainIntegrate(dmdt).value()
+            << endl;
     }
 }
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H

@@ -73,6 +73,14 @@ protected:
         // Mass transfer enabled
         Switch massTransfer_;
 
+        //- Interfacial Mass transfer rate
+        HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
+            iDmdt_;
+
+        //- Wall Mass transfer rate
+        HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
+            wDmdt_;
+
 
 public:
 
@@ -88,6 +96,9 @@ public:
 
     // Member Functions
 
+        //- Return the saturationModel
+        const saturationModel& saturation() const;
+
         //- Return the mass transfer matrices
         virtual autoPtr<phaseSystem::massTransferTable> massTransfer() const;
 

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

@@ -37,17 +37,7 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::AnisothermalPhaseModel
 )
 :
     BasePhaseModel(fluid, phaseName, index),
-    divU_
-    (
-        IOobject
-        (
-            IOobject::groupName("divU", this->name()),
-            fluid.mesh().time().timeName(),
-            fluid.mesh()
-        ),
-        fluid.mesh(),
-        dimensionedScalar("divU", dimless/dimTime, 0)
-    ),
+    divU_(NULL),
     K_
     (
         IOobject
@@ -105,7 +95,7 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
 
     volScalarField& he = this->thermo_->he();
 
-    return
+    tmp<fvScalarMatrix> tEEqn
     (
         fvm::ddt(alpha, this->rho(), he) + fvm::div(alphaRhoPhi, he)
       - fvm::Sp(contErr, he)
@@ -119,16 +109,23 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
            *fvc::interpolate(alphaEff),
             he
         )
-
-      + (
-            he.name() == this->thermo_->phasePropertyName("e")
-          ? fvc::ddt(alpha)*this->thermo().p()
-          + fvc::div(alphaPhi, this->thermo().p())
-          : -alpha*this->fluid().dpdt()
-        )
      ==
         this->Sh()
     );
+
+    // Add the appropriate pressure-work term
+    if (he.name() == this->thermo_->phasePropertyName("e"))
+    {
+        tEEqn() +=
+            fvc::ddt(alpha)*this->thermo().p()
+          + fvc::div(alphaPhi, this->thermo().p());
+    }
+    else if (this->thermo_->dpdt())
+    {
+        tEEqn() -= alpha*this->fluid().dpdt();
+    }
+
+    return tEEqn;
 }
 
 
@@ -140,7 +137,7 @@ bool Foam::AnisothermalPhaseModel<BasePhaseModel>::compressible() const
 
 
 template<class BasePhaseModel>
-const Foam::volScalarField&
+const Foam::tmp<Foam::volScalarField>&
 Foam::AnisothermalPhaseModel<BasePhaseModel>::divU() const
 {
     return divU_;
@@ -149,7 +146,10 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::divU() const
 
 template<class BasePhaseModel>
 void
-Foam::AnisothermalPhaseModel<BasePhaseModel>::divU(const volScalarField& divU)
+Foam::AnisothermalPhaseModel<BasePhaseModel>::divU
+(
+    const tmp<volScalarField>& divU
+)
 {
     divU_ = divU;
 }

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

@@ -55,7 +55,7 @@ class AnisothermalPhaseModel
     // Private data
 
         //- Dilatation rate
-        volScalarField divU_;
+        tmp<volScalarField> divU_;
 
         //- Kinetic energy
         volScalarField K_;
@@ -95,10 +95,10 @@ public:
             virtual bool compressible() const;
 
             //- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
-            virtual const volScalarField& divU() const;
+            virtual const tmp<volScalarField>& divU() const;
 
             //- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
-            virtual void divU(const volScalarField& divU);
+            virtual void divU(const tmp<volScalarField>& divU);
 
             //- Return the phase kinetic energy
             virtual const volScalarField& K() const;

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

@@ -165,16 +165,17 @@ bool Foam::phaseModel::compressible() const
 }
 
 
-const Foam::volScalarField& Foam::phaseModel::divU() const
+const Foam::tmp<Foam::volScalarField>& Foam::phaseModel::divU() const
 {
     notImplemented("Foam::phaseModel::divU()");
-    return volScalarField::null();
+    static tmp<Foam::volScalarField> divU_(NULL);
+    return divU_;
 }
 
 
-void Foam::phaseModel::divU(const volScalarField& divU)
+void Foam::phaseModel::divU(const tmp<volScalarField>& divU)
 {
-    WarningIn("phaseModel::divU(const volScalarField& divU)")
+    WarningIn("phaseModel::divU(const tmp<volScalarField>& divU)")
         << "Attempt to set the dilatation rate of an incompressible phase"
         << endl;
 }

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

@@ -216,10 +216,10 @@ public:
             virtual bool compressible() const;
 
             //- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
-            virtual const volScalarField& divU() const;
+            virtual const tmp<volScalarField>& divU() const;
 
             //- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
-            virtual void divU(const volScalarField& divU);
+            virtual void divU(const tmp<volScalarField>& divU);
 
             //- Return the phase kinetic energy
             virtual const volScalarField& K() const;

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C

@@ -120,6 +120,19 @@ constTransport
     >
 > constRefRhoConstEThermoPhysics;
 
+typedef
+constTransport
+<
+    species::thermo
+    <
+        hRefConstThermo
+        <
+            rhoConst<specie>
+        >,
+        sensibleEnthalpy
+    >
+> constRefRhoConstHThermoPhysics;
+
 
 // pureMixture, sensibleEnthalpy:
 
@@ -229,6 +242,18 @@ makeReactionMixtureThermo
 );
 
 
+// multiComponentMixture, sensibleEnthalpy:
+
+makeReactionMixtureThermo
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    multiComponentMixture,
+    constRefRhoConstHThermoPhysics
+);
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C

@@ -204,7 +204,7 @@ void Foam::multiphaseSystem::solveAlphas()
                 mesh_
             ),
             mesh_,
-            dimensionedScalar("Sp", phase.divU().dimensions(), 0.0)
+            dimensionedScalar("Sp", divU.dimensions(), 0.0)
         );
 
         volScalarField::DimensionedInternalField Su
@@ -220,9 +220,9 @@ void Foam::multiphaseSystem::solveAlphas()
             divU*min(alpha, scalar(1))
         );
 
-        if (phase.compressible())
+        if (phase.divU().valid())
         {
-            const scalarField& dgdt = phase.divU();
+            const scalarField& dgdt = phase.divU()();
 
             forAll(dgdt, celli)
             {
@@ -247,9 +247,9 @@ void Foam::multiphaseSystem::solveAlphas()
 
             if (&phase2 == &phase) continue;
 
-            if (phase2.compressible())
+            if (phase2.divU().valid())
             {
-                const scalarField& dgdt2 = phase2.divU();
+                const scalarField& dgdt2 = phase2.divU()();
 
                 forAll(dgdt2, celli)
                 {
@@ -515,6 +515,12 @@ Foam::multiphaseSystem::~multiphaseSystem()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+bool Foam::multiphaseSystem::transfersMass(const phaseModel& phase) const
+{
+    return false;
+}
+
+
 Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::surfaceTension
 (
     const phaseModel& phase1

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H

@@ -175,6 +175,9 @@ public:
             const PtrList<volScalarField>& rAUs
         ) const = 0;
 
+        //- Return true if there is mass transfer for phase
+        virtual bool transfersMass(const phaseModel& phase) const;
+
         //- Return the total interfacial mass transfer rate for phase
         virtual tmp<volScalarField> dmdt(const phaseModel& phase) const = 0;
 

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H

@@ -311,7 +311,7 @@ while (pimple.correct())
                     phasei,
                     (
                         (
-                            phase.continuityError() - fluid.dmdt(phase)
+                            phase.continuityError()
                           - fvc::Sp
                             (
                                 fvc::ddt(alpha) + fvc::div(phase.alphaPhi()),
@@ -340,7 +340,7 @@ while (pimple.correct())
                     phasei,
                     (
                         (