reactingEulerFoam: Updated phase-change support

Patches provided by Juho Peltola

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

@@ -67,26 +67,7 @@ ThermalPhaseChangePhaseSystem
                     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::READ_IF_PRESENT,
                     IOobject::AUTO_WRITE
                 ),
                 this->mesh(),
@@ -115,6 +96,101 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::saturation() const
 }
 
 
+template<class BasePhaseSystem>
+Foam::autoPtr<Foam::phaseSystem::heatTransferTable>
+Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::heatTransfer() const
+{
+    typedef compressible::alphatPhaseChangeWallFunctionFvPatchScalarField
+        alphatPhaseChangeWallFunction;
+
+    autoPtr<phaseSystem::heatTransferTable> eqnsPtr =
+        Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::heatTransfer();
+
+    phaseSystem::heatTransferTable& eqns = eqnsPtr();
+
+    // Accumulate mDotL contributions from boundaries
+    forAllConstIter
+    (
+        phaseSystem::phasePairTable,
+        this->phasePairs_,
+        phasePairIter
+    )
+    {
+        const phasePair& pair(phasePairIter());
+
+        if (pair.ordered())
+        {
+            continue;
+        }
+
+        const phaseModel& phase = pair.phase1();
+        const phaseModel& otherPhase = pair.phase2();
+
+        volScalarField mDotL
+        (
+            IOobject
+            (
+                "mDotL",
+                phase.mesh().time().timeName(),
+                phase.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            phase.mesh(),
+            dimensionedScalar("",dimensionSet(1,-1,-3,0,0),0.0)
+        );
+
+        if
+        (
+            otherPhase.mesh().foundObject<volScalarField>
+            (
+                "alphat." +  otherPhase.name()
+            )
+        )
+        {
+            const volScalarField& alphat =
+                otherPhase.mesh().lookupObject<volScalarField>
+                (
+                    "alphat." +  otherPhase.name()
+                );
+
+            const fvPatchList& patches = this->mesh().boundary();
+            forAll(patches, patchi)
+            {
+                const fvPatch& currPatch = patches[patchi];
+
+                if
+                (
+                    isA<alphatPhaseChangeWallFunction>
+                    (
+                        alphat.boundaryField()[patchi]
+                    )
+                )
+                {
+                    const scalarField& patchMDotL =
+                        refCast<const alphatPhaseChangeWallFunction>
+                        (
+                            alphat.boundaryField()[patchi]
+                        ).mDotL();
+
+                    forAll(patchMDotL,facei)
+                    {
+                        label faceCelli = currPatch.faceCells()[facei];
+                        mDotL[faceCelli] = patchMDotL[facei];
+                    }
+                }
+            }
+        }
+
+        *eqns[otherPhase.name()] -= mDotL;
+
+    }
+
+    return eqnsPtr;
+}
+
+
 template<class BasePhaseSystem>
 Foam::autoPtr<Foam::phaseSystem::massTransferTable>
 Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
@@ -189,6 +265,8 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
 
     BasePhaseSystem::correctThermo();
 
+
+
     forAllConstIter
     (
         phaseSystem::phasePairTable,
@@ -206,6 +284,18 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
         const phaseModel& phase1 = pair.phase1();
         const phaseModel& phase2 = pair.phase2();
 
+        Info<< phase1.name() << " min/max T "
+            << min(phase1.thermo().T()).value()
+            << " - "
+            << max(phase1.thermo().T()).value()
+            << endl;
+
+        Info<< phase2.name() << " min/max T "
+            << min(phase2.thermo().T()).value()
+            << " - "
+            << max(phase2.thermo().T()).value()
+            << endl;
+
         const volScalarField& T1(phase1.thermo().T());
         const volScalarField& T2(phase2.thermo().T());
 
@@ -214,7 +304,6 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
 
         volScalarField& dmdt(*this->dmdt_[pair]);
         volScalarField& iDmdt(*this->iDmdt_[pair]);
-        volScalarField& wDmdt(*this->wDmdt_[pair]);
 
         volScalarField& Tf = *this->Tf_[pair];
 
@@ -287,6 +376,21 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
             << endl;
 
         // Accumulate dmdt contributions from boundaries
+        volScalarField wDmdt
+        (
+            IOobject
+            (
+                IOobject::groupName("wDmdt", pair.name()),
+                this->mesh().time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            this->mesh(),
+            dimensionedScalar("zero", dimDensity/dimTime, 0)
+        );
+
         if
         (
             phase2.mesh().foundObject<volScalarField>
@@ -295,9 +399,6 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
             )
         )
         {
-            scalar wDmdtRelax(this->mesh().fieldRelaxationFactor("wDmdt"));
-            wDmdt *= (1 - wDmdtRelax);
-
             const volScalarField& alphat =
                 phase2.mesh().lookupObject<volScalarField>
                 (
@@ -326,7 +427,7 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
                     forAll(patchDmdt,facei)
                     {
                         label faceCelli = currPatch.faceCells()[facei];
-                        wDmdt[faceCelli] += wDmdtRelax*patchDmdt[facei];
+                        wDmdt[faceCelli] += patchDmdt[facei];
                     }
                 }
             }

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

@@ -77,10 +77,6 @@ protected:
         HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
             iDmdt_;
 
-        //- Wall Mass transfer rate
-        HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
-            wDmdt_;
-
 
 public:
 
@@ -99,6 +95,9 @@ public:
         //- Return the saturationModel
         const saturationModel& saturation() const;
 
+        //- Return the heat transfer matrices
+        virtual autoPtr<phaseSystem::heatTransferTable> heatTransfer() const;
+
         //- Return the mass transfer matrices
         virtual autoPtr<phaseSystem::massTransferTable> massTransfer() const;
 

src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatPhaseChangeWallFunction/alphatPhaseChangeWallFunctionFvPatchScalarField.C

@@ -39,7 +39,6 @@ namespace compressible
 
 defineTypeNameAndDebug(alphatPhaseChangeWallFunctionFvPatchScalarField,0);
 
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 alphatPhaseChangeWallFunctionFvPatchScalarField::
@@ -50,7 +49,8 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(p, iF),
-    dmdt_(p.size(), 0.0)
+    dmdt_(p.size(), 0.0),
+    mDotL_(p.size(), 0.0)
 {}
 
 
@@ -63,12 +63,18 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(p, iF, dict),
-    dmdt_(p.size(), 0.0)
+    dmdt_(p.size(), 0.0),
+    mDotL_(p.size(), 0.0)
 {
     if (dict.found("dmdt"))
     {
         dmdt_ = scalarField("dmdt", dict, p.size());
     }
+
+    if (dict.found("mDotL"))
+    {
+        dmdt_ = scalarField("mDotL", dict, p.size());
+    }
 }
 
 
@@ -82,7 +88,8 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(ptf, p, iF, mapper),
-    dmdt_(ptf.dmdt_, mapper)
+    dmdt_(ptf.dmdt_, mapper),
+    mDotL_(ptf.mDotL_, mapper)
 {}
 
 
@@ -93,7 +100,8 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(awfpsf),
-    dmdt_(awfpsf.dmdt_)
+    dmdt_(awfpsf.dmdt_),
+    mDotL_(awfpsf.mDotL_)
 {}
 
 
@@ -105,7 +113,8 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(awfpsf, iF),
-    dmdt_(awfpsf.dmdt_)
+    dmdt_(awfpsf.dmdt_),
+    mDotL_(awfpsf.mDotL_)
 {}
 
 
@@ -115,6 +124,7 @@ void alphatPhaseChangeWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchField<scalar>::write(os);
     dmdt_.writeEntry("dmdt", os);
+    mDotL_.writeEntry("mDotL", os);
     writeEntry("value", os);
 }
 

src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatPhaseChangeWallFunction/alphatPhaseChangeWallFunctionFvPatchScalarField.H

@@ -31,7 +31,8 @@ Description
     Abstract base-class for all alphatWallFunctions supporting phase-change.
 
 SeeAlso
-    Foam::alphatWallFunction
+    Foam::fixedValueFvPatchScalarField
+    Foam::alphatWallFunctionFvPatchScalarField
 
 SourceFiles
     alphatPhaseChangeWallFunctionFvPatchScalarField.C
@@ -51,7 +52,7 @@ namespace compressible
 {
 
 /*---------------------------------------------------------------------------*\
-            Class alphatPhaseChangeWallFunctionFvPatchScalarField Declaration
+       Class alphatPhaseChangeWallFunctionFvPatchScalarField Declaration
 \*---------------------------------------------------------------------------*/
 
 class alphatPhaseChangeWallFunctionFvPatchScalarField
@@ -65,6 +66,9 @@ protected:
         //- Rate of phase-change
         scalarField dmdt_;
 
+        //- Latent heat of the phase-change
+        scalarField mDotL_;
+
 
 public:
 
@@ -122,6 +126,12 @@ public:
             return dmdt_;
         }
 
+        //- Return the enthelpy source due to phase-change
+        const scalarField& mDotL() const
+        {
+            return mDotL_;
+        }
+
         // Evaluation functions
 
             //- Update the coefficients associated with the patch field