Commit 85b27f67 authored by Henry Weller's avatar Henry Weller
Browse files

reactingMultiphaseEulerFoam: New Euler-Euler multiphase solver

Supporting any number of phases with heat and mass transfer, phase-change and reactions
parent 34f060cf
...@@ -18,7 +18,7 @@ forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter) ...@@ -18,7 +18,7 @@ forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
new fvVectorMatrix new fvVectorMatrix
( (
fvm::ddt(alpha, U) fvm::ddt(alpha, U)
+ fvm::div(phase.phiAlpha(), U) + fvm::div(phase.alphaPhi(), U)
+ (alpha/phase.rho())*fluid.Cvm(phase)* + (alpha/phase.rho())*fluid.Cvm(phase)*
( (
......
...@@ -61,7 +61,7 @@ void Foam::multiphaseSystem::calcAlphas() ...@@ -61,7 +61,7 @@ void Foam::multiphaseSystem::calcAlphas()
void Foam::multiphaseSystem::solveAlphas() void Foam::multiphaseSystem::solveAlphas()
{ {
PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size()); PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
int phasei = 0; int phasei = 0;
forAllIter(PtrDictionary<phaseModel>, phases_, iter) forAllIter(PtrDictionary<phaseModel>, phases_, iter)
...@@ -69,10 +69,10 @@ void Foam::multiphaseSystem::solveAlphas() ...@@ -69,10 +69,10 @@ void Foam::multiphaseSystem::solveAlphas()
phaseModel& phase1 = iter(); phaseModel& phase1 = iter();
volScalarField& alpha1 = phase1; volScalarField& alpha1 = phase1;
phase1.phiAlpha() = phase1.alphaPhi() =
dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0); dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0);
phiAlphaCorrs.set alphaPhiCorrs.set
( (
phasei, phasei,
new surfaceScalarField new surfaceScalarField
...@@ -87,7 +87,7 @@ void Foam::multiphaseSystem::solveAlphas() ...@@ -87,7 +87,7 @@ void Foam::multiphaseSystem::solveAlphas()
) )
); );
surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei]; surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
forAllIter(PtrDictionary<phaseModel>, phases_, iter2) forAllIter(PtrDictionary<phaseModel>, phases_, iter2)
{ {
...@@ -118,7 +118,7 @@ void Foam::multiphaseSystem::solveAlphas() ...@@ -118,7 +118,7 @@ void Foam::multiphaseSystem::solveAlphas()
"div(phir," + alpha2.name() + ',' + alpha1.name() + ')' "div(phir," + alpha2.name() + ',' + alpha1.name() + ')'
); );
phiAlphaCorr += fvc::flux alphaPhiCorr += fvc::flux
( (
-fvc::flux(-phir, phase2, phirScheme), -fvc::flux(-phir, phase2, phirScheme),
phase1, phase1,
...@@ -127,21 +127,21 @@ void Foam::multiphaseSystem::solveAlphas() ...@@ -127,21 +127,21 @@ void Foam::multiphaseSystem::solveAlphas()
} }
// Ensure that the flux at inflow BCs is preserved // Ensure that the flux at inflow BCs is preserved
forAll(phiAlphaCorr.boundaryField(), patchi) forAll(alphaPhiCorr.boundaryField(), patchi)
{ {
fvsPatchScalarField& phiAlphaCorrp = fvsPatchScalarField& alphaPhiCorrp =
phiAlphaCorr.boundaryField()[patchi]; alphaPhiCorr.boundaryField()[patchi];
if (!phiAlphaCorrp.coupled()) if (!alphaPhiCorrp.coupled())
{ {
const scalarField& phi1p = phase1.phi().boundaryField()[patchi]; const scalarField& phi1p = phase1.phi().boundaryField()[patchi];
const scalarField& alpha1p = alpha1.boundaryField()[patchi]; const scalarField& alpha1p = alpha1.boundaryField()[patchi];
forAll(phiAlphaCorrp, facei) forAll(alphaPhiCorrp, facei)
{ {
if (phi1p[facei] < 0) if (phi1p[facei] < 0)
{ {
phiAlphaCorrp[facei] = alpha1p[facei]*phi1p[facei]; alphaPhiCorrp[facei] = alpha1p[facei]*phi1p[facei];
} }
} }
} }
...@@ -153,7 +153,7 @@ void Foam::multiphaseSystem::solveAlphas() ...@@ -153,7 +153,7 @@ void Foam::multiphaseSystem::solveAlphas()
geometricOneField(), geometricOneField(),
phase1, phase1,
phi_, phi_,
phiAlphaCorr, alphaPhiCorr,
zeroField(), zeroField(),
zeroField(), zeroField(),
1, 1,
...@@ -164,7 +164,7 @@ void Foam::multiphaseSystem::solveAlphas() ...@@ -164,7 +164,7 @@ void Foam::multiphaseSystem::solveAlphas()
phasei++; phasei++;
} }
MULES::limitSum(phiAlphaCorrs); MULES::limitSum(alphaPhiCorrs);
volScalarField sumAlpha volScalarField sumAlpha
( (
...@@ -184,19 +184,19 @@ void Foam::multiphaseSystem::solveAlphas() ...@@ -184,19 +184,19 @@ void Foam::multiphaseSystem::solveAlphas()
{ {
phaseModel& phase1 = iter(); phaseModel& phase1 = iter();
surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei]; surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
phiAlpha += upwind<scalar>(mesh_, phi_).flux(phase1); alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase1);
MULES::explicitSolve MULES::explicitSolve
( (
geometricOneField(), geometricOneField(),
phase1, phase1,
phiAlpha, alphaPhi,
zeroField(), zeroField(),
zeroField() zeroField()
); );
phase1.phiAlpha() += phiAlpha; phase1.alphaPhi() += alphaPhi;
Info<< phase1.name() << " volume fraction, min, max = " Info<< phase1.name() << " volume fraction, min, max = "
<< phase1.weightedAverage(mesh_.V()).value() << phase1.weightedAverage(mesh_.V()).value()
......
...@@ -100,11 +100,11 @@ Foam::phaseModel::phaseModel ...@@ -100,11 +100,11 @@ Foam::phaseModel::phaseModel
mesh, mesh,
dimensionedVector("0", dimVelocity/dimTime, vector::zero) dimensionedVector("0", dimVelocity/dimTime, vector::zero)
), ),
phiAlpha_ alphaPhi_
( (
IOobject IOobject
( (
IOobject::groupName("phiAlpha", phaseName), IOobject::groupName("alphaPhi", phaseName),
mesh.time().timeName(), mesh.time().timeName(),
mesh mesh
), ),
......
...@@ -80,7 +80,7 @@ class phaseModel ...@@ -80,7 +80,7 @@ class phaseModel
volVectorField DDtU_; volVectorField DDtU_;
//- Volumetric flux of the phase //- Volumetric flux of the phase
surfaceScalarField phiAlpha_; surfaceScalarField alphaPhi_;
//- Volumetric flux for the phase //- Volumetric flux for the phase
autoPtr<surfaceScalarField> phiPtr_; autoPtr<surfaceScalarField> phiPtr_;
...@@ -198,14 +198,14 @@ public: ...@@ -198,14 +198,14 @@ public:
return phiPtr_(); return phiPtr_();
} }
const surfaceScalarField& phiAlpha() const const surfaceScalarField& alphaPhi() const
{ {
return phiAlpha_; return alphaPhi_;
} }
surfaceScalarField& phiAlpha() surfaceScalarField& alphaPhi()
{ {
return phiAlpha_; return alphaPhi_;
} }
//- Correct the phase properties //- Correct the phase properties
......
...@@ -6,5 +6,6 @@ wclean libso phaseSystems ...@@ -6,5 +6,6 @@ wclean libso phaseSystems
wclean libso interfacialModels wclean libso interfacialModels
wclean libso interfacialCompositionModels wclean libso interfacialCompositionModels
reactingTwoPhaseEulerFoam/Allwclean reactingTwoPhaseEulerFoam/Allwclean
reactingMultiphaseEulerFoam/Allwclean
# ----------------------------------------------------------------- end-of-file # ----------------------------------------------------------------- end-of-file
...@@ -8,5 +8,6 @@ wmake libso phaseSystems ...@@ -8,5 +8,6 @@ wmake libso phaseSystems
wmake libso interfacialModels wmake libso interfacialModels
wmake libso interfacialCompositionModels wmake libso interfacialCompositionModels
reactingTwoPhaseEulerFoam/Allwmake reactingTwoPhaseEulerFoam/Allwmake
reactingMultiphaseEulerFoam/Allwmake
# ----------------------------------------------------------------- end-of-file # ----------------------------------------------------------------- end-of-file
...@@ -21,7 +21,6 @@ EXE_INC = \ ...@@ -21,7 +21,6 @@ EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude -I$(LIB_SRC)/meshTools/lnInclude
LIB_LIBS = \ LIB_LIBS = \
-lreactingTwoPhaseSystem \
-lfluidThermophysicalModels \ -lfluidThermophysicalModels \
-lreactionThermophysicalModels \ -lreactionThermophysicalModels \
-lspecie -lspecie
...@@ -10,7 +10,6 @@ EXE_INC = \ ...@@ -10,7 +10,6 @@ EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude -I$(LIB_SRC)/meshTools/lnInclude
LIB_LIBS = \ LIB_LIBS = \
-lreactingTwoPhaseSystem \
-lcompressibleTransportModels \ -lcompressibleTransportModels \
-lfluidThermophysicalModels \ -lfluidThermophysicalModels \
-lspecie -lspecie
...@@ -132,7 +132,10 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const ...@@ -132,7 +132,10 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const
volScalarField muAlphaI volScalarField muAlphaI
( (
alpha1*rho1*nu1*alpha2*rho2*nu2 alpha1*rho1*nu1*alpha2*rho2*nu2
/(alpha1*rho1*nu1 + alpha2*rho2*nu2) /(
max(alpha1, pair_.phase1().residualAlpha())*rho1*nu1
+ max(alpha2, pair_.phase2().residualAlpha())*rho2*nu2
)
); );
volScalarField ReI volScalarField ReI
......
...@@ -164,6 +164,60 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::dmdt ...@@ -164,6 +164,60 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::dmdt
} }
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::dmdt
(
const Foam::phaseModel& phase
) const
{
tmp<volScalarField> tdmdt
(
new volScalarField
(
IOobject
(
IOobject::groupName("dmdt", phase.name()),
this->mesh_.time().timeName(),
this->mesh_
),
this->mesh_,
dimensionedScalar("zero", dimDensity/dimTime, 0)
)
);
forAllConstIter
(
phaseSystem::phasePairTable,
this->phasePairs_,
phasePairIter
)
{
const phasePair& pair(phasePairIter());
if (pair.ordered())
{
continue;
}
const phaseModel* phase1 = &pair.phase1();
const phaseModel* phase2 = &pair.phase2();
forAllConstIter(phasePair, pair, iter)
{
if (phase1 == &phase)
{
tdmdt() += this->dmdt(pair);
}
Swap(phase1, phase2);
}
}
return tdmdt;
}
template<class BasePhaseSystem> template<class BasePhaseSystem>
Foam::autoPtr<Foam::phaseSystem::momentumTransferTable> Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
......
...@@ -121,10 +121,10 @@ public: ...@@ -121,10 +121,10 @@ public:
// Member Functions // Member Functions
//- Return the interfacial mass flow rate //- Return the interfacial mass flow rate
virtual tmp<volScalarField> dmdt virtual tmp<volScalarField> dmdt(const phasePairKey& key) const;
(
const phasePairKey& key //- Return the total interfacial mass transfer rate for phase
) const; virtual tmp<volScalarField> dmdt(const phaseModel& phase) const;
//- Return the momentum transfer matrices //- Return the momentum transfer matrices
virtual autoPtr<phaseSystem::momentumTransferTable> virtual autoPtr<phaseSystem::momentumTransferTable>
......
...@@ -85,6 +85,30 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::dmdt ...@@ -85,6 +85,30 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::dmdt
} }
template<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
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)
)
);
}
template<class BasePhaseSystem> template<class BasePhaseSystem>
Foam::autoPtr<Foam::phaseSystem::heatTransferTable> Foam::autoPtr<Foam::phaseSystem::heatTransferTable>
Foam::HeatTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const Foam::HeatTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
......
...@@ -95,6 +95,9 @@ public: ...@@ -95,6 +95,9 @@ public:
const phasePairKey& key const phasePairKey& key
) const; ) const;
//- Return the total interfacial mass transfer rate for phase
virtual tmp<volScalarField> dmdt(const phaseModel& phase) const;
//- Return the heat transfer matrices //- Return the heat transfer matrices
virtual autoPtr<phaseSystem::heatTransferTable> virtual autoPtr<phaseSystem::heatTransferTable>
heatTransfer() const; heatTransfer() const;
......
...@@ -182,7 +182,6 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Kd ...@@ -182,7 +182,6 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Kd
) )
); );
// Add the implicit part of the drag force
forAllConstIter forAllConstIter
( (
phaseSystem::KdTable, phaseSystem::KdTable,
...@@ -442,7 +441,6 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const ...@@ -442,7 +441,6 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
} }
// Add the implicit part of the drag force // Add the implicit part of the drag force
/* ***HGW Currently this is handled in the pEqn
forAllConstIter forAllConstIter
( (
phaseSystem::KdTable, phaseSystem::KdTable,
...@@ -466,7 +464,6 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const ...@@ -466,7 +464,6 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
Swap(phase, otherPhase); Swap(phase, otherPhase);
} }
} }
*/
// Update the virtual mass coefficients // Update the virtual mass coefficients
forAllConstIter forAllConstIter
...@@ -499,20 +496,57 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const ...@@ -499,20 +496,57 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
const volVectorField& U = phase->U(); const volVectorField& U = phase->U();
const surfaceScalarField& phi = phase->phi(); const surfaceScalarField& phi = phase->phi();
*eqns[phase->name()] += *eqns[phase->name()] -=
- Vm Vm
*( *(
fvm::ddt(U) fvm::ddt(U)
+ fvm::div(phi, U) + fvm::div(phi, U)
- fvm::Sp(fvc::div(phi), U) - fvm::Sp(fvc::div(phi), U)
- otherPhase->DUDt() - otherPhase->DUDt()
) )
- this->MRF_.DDt(Vm, U - otherPhase->U()); + this->MRF_.DDt(Vm, U - otherPhase->U());
Swap(phase, otherPhase); Swap(phase, otherPhase);
} }
} }
return eqnsPtr;
}
template<class BasePhaseSystem>
Foam::autoPtr<Foam::PtrList<Foam::volVectorField> >
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
{
autoPtr<PtrList<volVectorField> > tFs
(
new PtrList<volVectorField>(this->phases().size())
);
PtrList<volVectorField>& Fs = tFs();
forAll(Fs, phasei)
{
Fs.set
(
phasei,
new volVectorField
(
IOobject
(
liftModel::typeName + ":F",
this->mesh_.time().timeName(),
this->mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
this->mesh_,
dimensionedVector("zero", liftModel::dimF, vector::zero)
)
);
}
// Add the lift force // Add the lift force
forAllConstIter forAllConstIter
( (
...@@ -525,8 +559,8 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const ...@@ -525,8 +559,8 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
const phasePair& pair(this->phasePairs_[liftModelIter.key()]); const phasePair& pair(this->phasePairs_[liftModelIter.key()]);