diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H
index 30ce7b221f3084ce391135271c292dd750056bdd..9a4174e8c13b7636e6b07ead46e30c2bee86c18d 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H
+++ b/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();
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
index c48f0b07c7ebf6c3d300492eca21709660b38088..05ed9e7d087616929a066c304209c0434d4791a1 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
+++ b/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);
+ }
}
}
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
index d62af01b53940703013ee81719186933b2c63e5c..99e16af8a0537ee0cf127498d939fd0519fd0285 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
+++ b/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;
+}
+
+
// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
index a9c3953543704b9f39cc0440435208a6d798b36d..c69156b356d17b0ac7863dd4be43865d7fd6dca4 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
+++ b/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);
};
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
index 3518650190401af7f7d8e1af4c70c2319457a016..d1cc94e6c193519f818da71daa2b04b6b14b6d76 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
+++ b/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;
+}
+
+
// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
index 0afb9f7855073a7383d0cc6850898cd75aadc447..e42a3a70b441cddf1713904a5e7f3f44bb22942e 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
+++ b/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
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
index 9144093db5a9784454e6717e077b1ee8e5ed4ff6..c89d0eec3e7b1a33403159bf57d91c6625cfda27 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
+++ b/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
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
index d7d21c430de267b01eb9fcf77394e59983212859..46b5ec20e0ffc83d272e94f4495cf27e59b8acd2 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
+++ b/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;
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
index 4095e3739f9337209ffdf78fe2630d9f00371864..adb3813bbd3a202b3435e0fae0b6e37eb1b93399 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
+++ b/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_;
-}
-
-
// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C
index 114f7416331008b040772396d48b9b5d461f9d67..00ebfcffb98faf5e534c8c93cf73c3a3942c94dc 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C
+++ b/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,16 +261,21 @@ void Foam::twoPhaseSystem::solve()
fvc::div(phi)*min(alpha1, scalar(1))
);
- forAll(dgdt, celli)
+ if (tdgdt.valid())
{
- if (dgdt[celli] > 0.0)
- {
- Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
- Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
- }
- else if (dgdt[celli] < 0.0)
+ scalarField& dgdt = tdgdt();
+
+ forAll(dgdt, celli)
{
- Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
+ if (dgdt[celli] > 0.0)
+ {
+ Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+ Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+ }
+ else if (dgdt[celli] < 0.0)
+ {
+ Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
+ }
}
}
@@ -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();