From c8135cee57b2d5fd534d85dcc78714b91ec6ed38 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Fri, 25 Sep 2015 19:00:07 +0100
Subject: [PATCH] reactingEulerFoam: Further improvements to the handling of
mass-transfer between incompressible and compressible phases
---
.../HeatTransferPhaseSystem.C | 10 +++++++
.../HeatTransferPhaseSystem.H | 3 ++
.../reactionThermo/hRefConstThermos.C | 10 +++++++
.../multiphaseSystem/multiphaseSystem.C | 6 ----
.../multiphaseSystem/multiphaseSystem.H | 2 +-
.../reactingMultiphaseEulerFoam/pU/pEqn.H | 8 +++--
.../reactingTwoPhaseEulerFoam/pU/pEqn.H | 29 ++++++++++++++++---
.../reactingTwoPhaseEulerFoam/pUf/pEqn.H | 29 ++++++++++++++++---
.../twoPhaseSystem/twoPhaseSystem.C | 6 ++++
.../twoPhaseSystem/twoPhaseSystem.H | 6 ++++
10 files changed, 92 insertions(+), 17 deletions(-)
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C
index 95636c37497..07b41eb7710 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C
@@ -57,6 +57,16 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::~HeatTransferPhaseSystem()
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+template<class BasePhaseSystem>
+bool Foam::HeatTransferPhaseSystem<BasePhaseSystem>::transfersMass
+(
+ const phaseModel& phase
+) const
+{
+ return false;
+}
+
+
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::HeatTransferPhaseSystem<BasePhaseSystem>::dmdt
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.H
index cef3b696bd6..0c908239fdb 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.H
@@ -89,6 +89,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
(
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C
index 76c2351109c..f43303a2b7e 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C
@@ -254,6 +254,16 @@ makeReactionMixtureThermo
);
+makeReactionMixtureThermo
+(
+ rhoThermo,
+ rhoReactionThermo,
+ heRhoThermo,
+ multiComponentMixture,
+ constRefGasHThermoPhysics
+);
+
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index a79e9393b5a..5512fd65599 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -515,12 +515,6 @@ 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
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
index a9390f2589e..e33bf8c4878 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
@@ -176,7 +176,7 @@ public:
) const = 0;
//- Return true if there is mass transfer for phase
- virtual bool transfersMass(const phaseModel& phase) const;
+ virtual bool transfersMass(const phaseModel& phase) const = 0;
//- Return the total interfacial mass transfer rate for phase
virtual tmp<volScalarField> dmdt(const phaseModel& phase) const = 0;
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
index 230b6501065..2703d75ef1f 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
@@ -357,11 +357,15 @@ while (pimple.correct())
{
if (pEqnComps.set(phasei))
{
- pEqnComps[phasei] -= fluid.dmdt(phase);
+ pEqnComps[phasei] -= fluid.dmdt(phase)/rho;
}
else
{
- pEqnComps.set(phasei, fvm::Su(-fluid.dmdt(phase), rho));
+ pEqnComps.set
+ (
+ phasei,
+ fvm::Su(-fluid.dmdt(phase)/rho, p_rgh)
+ );
}
}
}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
index c8a9e282889..1af167cc408 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
@@ -247,7 +247,7 @@ while (pimple.correct())
pEqnComp1 =
(
- phase1.continuityError() - fluid.dmdt()
+ phase1.continuityError()
- fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
)/rho1
+ (alpha1/rho1)*correction
@@ -270,7 +270,7 @@ while (pimple.correct())
pEqnComp2 =
(
- phase2.continuityError() + fluid.dmdt()
+ phase2.continuityError()
- fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
)/rho2
+ (alpha2/rho2)*correction
@@ -288,7 +288,7 @@ while (pimple.correct())
{
pEqnComp1 =
(
- phase1.continuityError() - fluid.dmdt()
+ phase1.continuityError()
- fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
)/rho1
+ (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
@@ -298,13 +298,34 @@ while (pimple.correct())
{
pEqnComp2 =
(
- phase2.continuityError() + fluid.dmdt()
+ phase2.continuityError()
- fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
)/rho2
+ (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
}
}
+ if (fluid.transfersMass())
+ {
+ if (pEqnComp1.valid())
+ {
+ pEqnComp1() -= fluid.dmdt()/rho1;
+ }
+ else
+ {
+ pEqnComp1 = fvm::Su(-fluid.dmdt()/rho1, p_rgh);
+ }
+
+ if (pEqnComp2.valid())
+ {
+ pEqnComp2() += fluid.dmdt()/rho2;
+ }
+ else
+ {
+ pEqnComp2 = fvm::Su(fluid.dmdt()/rho2, p_rgh);
+ }
+ }
+
// Cache p prior to solve for density update
volScalarField p_rgh_0(p_rgh);
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
index c3f7d472fa8..e2edf14bed0 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
@@ -231,7 +231,7 @@ while (pimple.correct())
{
pEqnComp1 =
(
- phase1.continuityError() - fluid.dmdt()
+ phase1.continuityError()
- fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
)/rho1
+ (alpha1/rho1)*correction
@@ -247,7 +247,7 @@ while (pimple.correct())
{
pEqnComp2 =
(
- phase2.continuityError() + fluid.dmdt()
+ phase2.continuityError()
- fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
)/rho2
+ (alpha2/rho2)*correction
@@ -265,7 +265,7 @@ while (pimple.correct())
{
pEqnComp1 =
(
- phase1.continuityError() - fluid.dmdt()
+ phase1.continuityError()
- fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
)/rho1
+ (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
@@ -275,13 +275,34 @@ while (pimple.correct())
{
pEqnComp2 =
(
- phase2.continuityError() + fluid.dmdt()
+ phase2.continuityError()
- fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
)/rho2
+ (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
}
}
+ if (fluid.transfersMass())
+ {
+ if (pEqnComp1.valid())
+ {
+ pEqnComp1() -= fluid.dmdt()/rho1;
+ }
+ else
+ {
+ pEqnComp1 = fvm::Su(-fluid.dmdt()/rho1, p_rgh);
+ }
+
+ if (pEqnComp2.valid())
+ {
+ pEqnComp2() += fluid.dmdt()/rho2;
+ }
+ else
+ {
+ pEqnComp2 = fvm::Su(fluid.dmdt()/rho2, p_rgh);
+ }
+ }
+
// Cache p prior to solve for density update
volScalarField p_rgh_0("p_rgh_0", p_rgh);
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index 187ca4c2388..7b7ddee6f49 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -168,6 +168,12 @@ Foam::twoPhaseSystem::D() const
}
+bool Foam::twoPhaseSystem::transfersMass() const
+{
+ return transfersMass(phase1());
+}
+
+
Foam::tmp<Foam::volScalarField>
Foam::twoPhaseSystem::dmdt() const
{
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
index deca55fa6e5..3c4c98f5892 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
@@ -78,6 +78,9 @@ class twoPhaseSystem
// Multiplies the phase-fraction gradient
virtual tmp<volScalarField> D(const phasePairKey& key) const = 0;
+ //- Return true if there is mass transfer for phase
+ virtual bool transfersMass(const phaseModel& phase) const = 0;
+
//- Return the interfacial mass flow rate for phase pair
virtual tmp<volScalarField> dmdt(const phasePairKey& key) const = 0;
@@ -192,6 +195,9 @@ public:
// Multiplies the phase-fraction gradient
tmp<volScalarField> D() const;
+ //- Return true if there is mass transfer
+ bool transfersMass() const;
+
//- Return the interfacial mass flow rate
tmp<volScalarField> dmdt() const;
--
GitLab