reactingTwoPhaseEulerFoam: Improved support for boiling/condensation

Includes many contributions from Juho Peltola

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/EEqns.H

@@ -47,3 +47,9 @@
 }
 
 fluid.correctThermo();
+
+Info<< "  phase1.thermo().T(): " << min(phase1.thermo().T()).value()
+    << " - " << max(phase1.thermo().T()).value() << endl;
+
+Info<< "  phase2.thermo().T(): " << min(phase2.thermo().T()).value()
+    << " - " << max(phase2.thermo().T()).value() << endl;

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C

@@ -60,13 +60,13 @@ Foam::heatTransferModels::RanzMarshall::~RanzMarshall()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::heatTransferModels::RanzMarshall::K() const
+Foam::heatTransferModels::RanzMarshall::K(const scalar residualAlpha) const
 {
     volScalarField Nu(scalar(2) + 0.6*sqrt(pair_.Re())*cbrt(pair_.Pr()));
 
     return
         6.0
-       *max(pair_.dispersed(), residualAlpha_)
+       *max(pair_.dispersed(), residualAlpha)
        *pair_.continuous().kappa()
        *Nu
        /sqr(pair_.dispersed().d());

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.H

@@ -79,7 +79,7 @@ public:
     // Member Functions
 
         //- The heat transfer function K used in the enthalpy equation
-        tmp<volScalarField> K() const;
+        tmp<volScalarField> K(const scalar residualAlpha) const;
 };
 
 

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/heatTransferModel/heatTransferModel.C

@@ -65,4 +65,13 @@ Foam::heatTransferModel::~heatTransferModel()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::heatTransferModel::K() const
+{
+    return K(residualAlpha_.value());
+}
+
+
 // ************************************************************************* //

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/heatTransferModel/heatTransferModel.H

@@ -117,7 +117,13 @@ public:
         //- The heat transfer function K used in the enthalpy equation
         //    ddt(alpha1*rho1*ha) + ... = ... K*(Ta - Tb)
         //    ddt(alpha2*rho2*hb) + ... = ... K*(Tb - Ta)
-        virtual tmp<volScalarField> K() const = 0;
+        tmp<volScalarField> K() const;
+
+        //- The heat transfer function K used in the enthalpy equation
+        //    ddt(alpha1*rho1*ha) + ... = ... K*(Ta - Tb)
+        //    ddt(alpha2*rho2*hb) + ... = ... K*(Tb - Ta)
+        //  with a specified residual volume fraction
+        virtual tmp<volScalarField> K(const scalar residualAlpha) const = 0;
 };
 
 

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/sphericalHeatTransfer/sphericalHeatTransfer.C

@@ -65,11 +65,14 @@ Foam::heatTransferModels::sphericalHeatTransfer::~sphericalHeatTransfer()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::heatTransferModels::sphericalHeatTransfer::K() const
+Foam::heatTransferModels::sphericalHeatTransfer::K
+(
+    const scalar residualAlpha
+) const
 {
     return
         60.0
-       *max(pair_.dispersed(), residualAlpha_)
+       *max(pair_.dispersed(), residualAlpha)
        *pair_.continuous().kappa()
        /sqr(pair_.dispersed().d());
 }

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/interfacialModels/heatTransferModels/sphericalHeatTransfer/sphericalHeatTransfer.H

@@ -79,7 +79,7 @@ public:
     // Member Functions
 
         //- The heat transfer function K used in the enthalpy equation
-        tmp<volScalarField> K() const;
+        tmp<volScalarField> K(const scalar residualAlpha) const;
 };
 
 

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

@@ -223,6 +223,66 @@ Foam::BlendedInterfacialModel<ModelType>::K() const
 }
 
 
+template<class ModelType>
+Foam::tmp<Foam::volScalarField>
+Foam::BlendedInterfacialModel<ModelType>::K(const scalar residualAlpha) const
+{
+    tmp<volScalarField> f1, f2;
+
+    if (model_.valid() || model1In2_.valid())
+    {
+        f1 = blending_.f1(phase1_, phase2_);
+    }
+
+    if (model_.valid() || model2In1_.valid())
+    {
+        f2 = blending_.f2(phase1_, phase2_);
+    }
+
+    tmp<volScalarField> x
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                ModelType::typeName + ":K",
+                phase1_.mesh().time().timeName(),
+                phase1_.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            phase1_.mesh(),
+            dimensionedScalar("zero", ModelType::dimK, 0)
+        )
+    );
+
+    if (model_.valid())
+    {
+        x() += model_->K(residualAlpha)*(scalar(1) - f1() - f2());
+    }
+    if (model1In2_.valid())
+    {
+        x() += model1In2_->K(residualAlpha)*f1;
+    }
+    if (model2In1_.valid())
+    {
+        x() += model2In1_->K(residualAlpha)*f2;
+    }
+
+    if
+    (
+        correctFixedFluxBCs_
+     && (model_.valid() || model1In2_.valid() || model2In1_.valid())
+    )
+    {
+        correctFixedFluxBCs(x());
+    }
+
+    return x;
+}
+
+
 template<class ModelType>
 Foam::tmp<Foam::surfaceScalarField>
 Foam::BlendedInterfacialModel<ModelType>::Kf() const

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

@@ -136,6 +136,9 @@ public:
         //- Return the blended force coefficient
         tmp<volScalarField> K() const;
 
+        //- Return the blended force coefficient with a specified residual alpha
+        tmp<volScalarField> K(const scalar residualAlpha) const;
+
         //- Return the face blended force coefficient
         tmp<surfaceScalarField> Kf() const;
 

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/Make/files

@@ -29,4 +29,6 @@ twoPhaseSystem/twoPhaseSystem.C
 twoPhaseSystem/newTwoPhaseSystem.C
 twoPhaseSystem/twoPhaseSystems.C
 
+reactionThermo/hRefConstThermos.C
+
 LIB = $(FOAM_LIBBIN)/libreactingTwoPhaseSystem

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

@@ -192,11 +192,11 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
         const volVectorField& U2(pair.phase2().U());
 
         const volScalarField dmdt(this->dmdt(pair));
-        const volScalarField dmdt12(posPart(dmdt));
+        const volScalarField dmdt21(posPart(dmdt));
-        const volScalarField dmdt21(negPart(dmdt));
+        const volScalarField dmdt12(negPart(dmdt));
 
-        *eqns[pair.phase1().name()] += fvm::Sp(dmdt21, U1) - dmdt21*U2;
+        *eqns[pair.phase1().name()] += dmdt21*U2 - fvm::Sp(dmdt21, U1);
-        *eqns[pair.phase2().name()] += dmdt12*U1 - fvm::Sp(dmdt12, U2);
+        *eqns[pair.phase2().name()] -= dmdt12*U1 - fvm::Sp(dmdt12, U2);
     }
 
     return eqnsPtr;
@@ -283,7 +283,7 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
         }
     }
 
-    // Source term due to mass trasfer
+    // Source term due to mass transfer
     forAllConstIter
     (
         phaseSystem::phasePairTable,
@@ -308,17 +308,19 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
         const volScalarField& K2(phase2.K());
 
         const volScalarField dmdt(this->dmdt(pair));
-        const volScalarField dmdt12(posPart(dmdt));
+        const volScalarField dmdt21(posPart(dmdt));
-        const volScalarField dmdt21(negPart(dmdt));
+        const volScalarField dmdt12(negPart(dmdt));
         const volScalarField& Tf(*Tf_[pair]);
 
         *eqns[phase1.name()] +=
-            fvm::Sp(dmdt21, he1) + dmdt21*K1
+            dmdt21*(phase1.thermo().he(phase1.thermo().p(), Tf))
-          - dmdt21*(phase2.thermo().he(phase2.thermo().p(), Tf) + K2);
+          - fvm::Sp(dmdt21, he1)
+          + dmdt21*(K2 - K1);
 
-        *eqns[phase2.name()] +=
+        *eqns[phase2.name()] -=
-            dmdt12*(phase1.thermo().he(phase1.thermo().p(), Tf) + K1)
+            dmdt12*(phase2.thermo().he(phase2.thermo().p(), Tf))
-          - fvm::Sp(dmdt12, he2) - dmdt12*K2;
+          - fvm::Sp(dmdt12, he2)
+          + dmdt12*(K1 - K2);
     }
 
     return eqnsPtr;

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

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "ThermalPhaseChangePhaseSystem.H"
+#include "fvCFD.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -36,7 +37,8 @@ ThermalPhaseChangePhaseSystem
 :
     HeatAndMassTransferPhaseSystem<BasePhaseSystem>(mesh),
     volatile_(this->lookup("volatile")),
-    saturationModel_(saturationModel::New(this->subDict("saturationModel")))
+    saturationModel_(saturationModel::New(this->subDict("saturationModel"))),
+    massTransfer_(this->lookup("massTransfer"))
 {}
 
 
@@ -143,37 +145,79 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
         const phaseModel& phase1 = pair.phase1();
         const phaseModel& phase2 = pair.phase2();
 
-        volScalarField& Tf = *this->Tf_[pair];
-        Tf = saturationModel_->Tsat(phase1.thermo().p());
-
-        Info<< "Tf." << pair.name()
-            << ": min = " << min(Tf.internalField())
-            << ", mean = " << average(Tf.internalField())
-            << ", max = " << max(Tf.internalField())
-            << endl;
-
-        volScalarField& dmdt(*this->dmdt_[pair]);
-
-        volScalarField H1(this->heatTransferModels_[pair][pair.first()]->K());
-        volScalarField H2(this->heatTransferModels_[pair][pair.second()]->K());
-
         const volScalarField& T1(phase1.thermo().T());
         const volScalarField& T2(phase2.thermo().T());
 
         const volScalarField& he1(phase1.thermo().he());
         const volScalarField& he2(phase2.thermo().he());
 
-        volScalarField hef2(phase2.thermo().he(phase2.thermo().p(), Tf));
+        volScalarField& dmdt(*this->dmdt_[pair]);
-        volScalarField hef1(phase1.thermo().he(phase1.thermo().p(), Tf));
 
-        dmdt =
+        volScalarField& Tf = *this->Tf_[pair];
-            (H2*(Tf - T1) + H1*(Tf - T2))
+
-           /min
+        volScalarField hef1(phase1.thermo().he(phase1.thermo().p(), Tf));
+        volScalarField hef2(phase2.thermo().he(phase2.thermo().p(), Tf));
+
+        if (massTransfer_ )
+        {
+            volScalarField H1
+            (
+                this->heatTransferModels_[pair][pair.first()]->K(0)
+            );
+
+            volScalarField H2
+            (
+                this->heatTransferModels_[pair][pair.second()]->K(0)
+            );
+
+            Tf = saturationModel_->Tsat(phase1.thermo().p());
+
+            dmdt =
+                (H1*(Tf - T1) + H2*(Tf - T2))
+               /min
+                (
+                    (pos(dmdt)*he2 + neg(dmdt)*hef2)
+                  - (neg(dmdt)*he1 + pos(dmdt)*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()
+                << endl;
+        }
+        else
+        {
+            dmdt == dimensionedScalar("0", dmdt.dimensions(), 0);
+        }
+
+        volScalarField H1(this->heatTransferModels_[pair][pair.first()]->K());
+        volScalarField H2(this->heatTransferModels_[pair][pair.second()]->K());
+
+        // Limit the H[12] boundary field to avoid /0
+        const scalar HLimit = 1e-4;
+        H1.boundaryField() =
+            max(H1.boundaryField(), phase1.boundaryField()*HLimit);
+        H2.boundaryField() =
+            max(H2.boundaryField(), phase2.boundaryField()*HLimit);
+
+        volScalarField mDotL
+        (
+            dmdt*
             (
                 (pos(dmdt)*he2 + neg(dmdt)*hef2)
-              - (neg(dmdt)*he1 + pos(dmdt)*hef1),
+              - (neg(dmdt)*he1 + pos(dmdt)*hef1)
-                0.3*(hef1 - hef2)
+            )
-            );
+        );
+        Tf = (H1*T1 + H2*T2 + mDotL)/(H1 + H2);
+
+        Info<< "Tf." << pair.name()
+            << ": min = " << min(Tf.internalField())
+            << ", mean = " << average(Tf.internalField())
+            << ", max = " << max(Tf.internalField())
+            << endl;
     }
 }
 

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

@@ -43,6 +43,7 @@ SourceFiles
 
 #include "HeatAndMassTransferPhaseSystem.H"
 #include "saturationModel.H"
+#include "Switch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -69,6 +70,9 @@ protected:
         //- The saturation model used to evaluate Tsat = Tf
         autoPtr<saturationModel> saturationModel_;
 
+        // Mass transfer enabled
+        Switch massTransfer_;
+
 
 public:
 

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

@@ -0,0 +1,226 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 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 "makeReactionThermo.H"
+#include "makeThermo.H"
+
+#include "rhoReactionThermo.H"
+#include "heRhoThermo.H"
+
+#include "specie.H"
+#include "perfectGas.H"
+#include "perfectFluid.H"
+#include "rhoConst.H"
+#include "incompressiblePerfectGas.H"
+#include "hConstThermo.H"
+#include "janafThermo.H"
+#include "sensibleEnthalpy.H"
+#include "absoluteEnthalpy.H"
+#include "absoluteInternalEnergy.H"
+#include "thermo.H"
+
+#include "constTransport.H"
+#include "sutherlandTransport.H"
+
+#include "pureMixture.H"
+#include "homogeneousMixture.H"
+#include "inhomogeneousMixture.H"
+#include "veryInhomogeneousMixture.H"
+#include "multiComponentMixture.H"
+#include "reactingMixture.H"
+#include "singleStepReactingMixture.H"
+
+#include "thermoPhysicsTypes.H"
+
+#include "hRefConstThermo.H"
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Thermo type typedefs:
+
+typedef
+constTransport
+<
+    species::thermo
+    <
+        hRefConstThermo
+        <
+            perfectGas<specie>
+        >,
+        sensibleEnthalpy
+    >
+> constRefGasHThermoPhysics;
+
+typedef
+constTransport
+<
+    species::thermo
+    <
+        hRefConstThermo
+        <
+            perfectFluid<specie>
+        >,
+        sensibleEnthalpy
+    >
+> constRefFluidHThermoPhysics;
+
+typedef
+constTransport
+<
+    species::thermo
+    <
+        hRefConstThermo
+        <
+            perfectGas<specie>
+        >,
+        sensibleInternalEnergy
+    >
+> constRefGasEThermoPhysics;
+
+typedef
+constTransport
+<
+    species::thermo
+    <
+        hRefConstThermo
+        <
+            perfectFluid<specie>
+        >,
+        sensibleInternalEnergy
+    >
+> constRefFluidEThermoPhysics;
+
+
+// pureMixture, sensibleEnthalpy:
+
+makeThermo
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleEnthalpy,
+    hRefConstThermo,
+    perfectGas,
+    specie
+);
+
+makeThermo
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleEnthalpy,
+    hRefConstThermo,
+    perfectFluid,
+    specie
+);
+
+makeThermo
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleEnthalpy,
+    hRefConstThermo,
+    rhoConst,
+    specie
+);
+
+
+// pureMixture, sensibleInternalEnergy:
+
+makeThermo
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleInternalEnergy,
+    hRefConstThermo,
+    perfectGas,
+    specie
+);
+
+makeThermo
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleInternalEnergy,
+    hRefConstThermo,
+    perfectFluid,
+    specie
+);
+
+makeThermo
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleInternalEnergy,
+    hRefConstThermo,
+    rhoConst,
+    specie
+);
+
+
+// multiComponentMixture, sensibleInternalEnergy:
+
+makeReactionMixtureThermo
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    multiComponentMixture,
+    constRefGasEThermoPhysics
+);
+
+makeReactionMixtureThermo
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    multiComponentMixture,
+    constRefFluidEThermoPhysics
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //