Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
reactingMultiphaseEulerFoam: New Euler-Euler multiphase solver
Supporting any number of phases with heat and mass transfer, phase-change and reactions
  • Loading branch information
Henry Weller committed Sep 11, 2015
1 parent 7c87973 commit cc5f67a
Show file tree
Hide file tree
Showing 172 changed files with 13,507 additions and 148 deletions.
Expand Up @@ -18,7 +18,7 @@ forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
new fvVectorMatrix
(
fvm::ddt(alpha, U)
+ fvm::div(phase.phiAlpha(), U)
+ fvm::div(phase.alphaPhi(), U)

+ (alpha/phase.rho())*fluid.Cvm(phase)*
(
Expand Down
Expand Up @@ -61,18 +61,18 @@ void Foam::multiphaseSystem::calcAlphas()

void Foam::multiphaseSystem::solveAlphas()
{
PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size());
PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
int phasei = 0;

forAllIter(PtrDictionary<phaseModel>, phases_, iter)
{
phaseModel& phase1 = iter();
volScalarField& alpha1 = phase1;

phase1.phiAlpha() =
phase1.alphaPhi() =
dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0);

phiAlphaCorrs.set
alphaPhiCorrs.set
(
phasei,
new surfaceScalarField
Expand All @@ -87,7 +87,7 @@ void Foam::multiphaseSystem::solveAlphas()
)
);

surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];
surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];

forAllIter(PtrDictionary<phaseModel>, phases_, iter2)
{
Expand Down Expand Up @@ -118,7 +118,7 @@ void Foam::multiphaseSystem::solveAlphas()
"div(phir," + alpha2.name() + ',' + alpha1.name() + ')'
);

phiAlphaCorr += fvc::flux
alphaPhiCorr += fvc::flux
(
-fvc::flux(-phir, phase2, phirScheme),
phase1,
Expand All @@ -127,21 +127,21 @@ void Foam::multiphaseSystem::solveAlphas()
}

// Ensure that the flux at inflow BCs is preserved
forAll(phiAlphaCorr.boundaryField(), patchi)
forAll(alphaPhiCorr.boundaryField(), patchi)
{
fvsPatchScalarField& phiAlphaCorrp =
phiAlphaCorr.boundaryField()[patchi];
fvsPatchScalarField& alphaPhiCorrp =
alphaPhiCorr.boundaryField()[patchi];

if (!phiAlphaCorrp.coupled())
if (!alphaPhiCorrp.coupled())
{
const scalarField& phi1p = phase1.phi().boundaryField()[patchi];
const scalarField& alpha1p = alpha1.boundaryField()[patchi];

forAll(phiAlphaCorrp, facei)
forAll(alphaPhiCorrp, facei)
{
if (phi1p[facei] < 0)
{
phiAlphaCorrp[facei] = alpha1p[facei]*phi1p[facei];
alphaPhiCorrp[facei] = alpha1p[facei]*phi1p[facei];
}
}
}
Expand All @@ -153,7 +153,7 @@ void Foam::multiphaseSystem::solveAlphas()
geometricOneField(),
phase1,
phi_,
phiAlphaCorr,
alphaPhiCorr,
zeroField(),
zeroField(),
1,
Expand All @@ -164,7 +164,7 @@ void Foam::multiphaseSystem::solveAlphas()
phasei++;
}

MULES::limitSum(phiAlphaCorrs);
MULES::limitSum(alphaPhiCorrs);

volScalarField sumAlpha
(
Expand All @@ -184,19 +184,19 @@ void Foam::multiphaseSystem::solveAlphas()
{
phaseModel& phase1 = iter();

surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei];
phiAlpha += upwind<scalar>(mesh_, phi_).flux(phase1);
surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase1);

MULES::explicitSolve
(
geometricOneField(),
phase1,
phiAlpha,
alphaPhi,
zeroField(),
zeroField()
);

phase1.phiAlpha() += phiAlpha;
phase1.alphaPhi() += alphaPhi;

Info<< phase1.name() << " volume fraction, min, max = "
<< phase1.weightedAverage(mesh_.V()).value()
Expand Down
Expand Up @@ -100,11 +100,11 @@ Foam::phaseModel::phaseModel
mesh,
dimensionedVector("0", dimVelocity/dimTime, vector::zero)
),
phiAlpha_
alphaPhi_
(
IOobject
(
IOobject::groupName("phiAlpha", phaseName),
IOobject::groupName("alphaPhi", phaseName),
mesh.time().timeName(),
mesh
),
Expand Down
Expand Up @@ -80,7 +80,7 @@ class phaseModel
volVectorField DDtU_;

//- Volumetric flux of the phase
surfaceScalarField phiAlpha_;
surfaceScalarField alphaPhi_;

//- Volumetric flux for the phase
autoPtr<surfaceScalarField> phiPtr_;
Expand Down Expand Up @@ -198,14 +198,14 @@ public:
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
Expand Down
Expand Up @@ -6,5 +6,6 @@ wclean libso phaseSystems
wclean libso interfacialModels
wclean libso interfacialCompositionModels
reactingTwoPhaseEulerFoam/Allwclean
reactingMultiphaseEulerFoam/Allwclean

# ----------------------------------------------------------------- end-of-file
1 change: 1 addition & 0 deletions applications/solvers/multiphase/reactingEulerFoam/Allwmake
Expand Up @@ -8,5 +8,6 @@ wmake libso phaseSystems
wmake libso interfacialModels
wmake libso interfacialCompositionModels
reactingTwoPhaseEulerFoam/Allwmake
reactingMultiphaseEulerFoam/Allwmake

# ----------------------------------------------------------------- end-of-file
Expand Up @@ -21,7 +21,6 @@ EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude

LIB_LIBS = \
-lreactingTwoPhaseSystem \
-lfluidThermophysicalModels \
-lreactionThermophysicalModels \
-lspecie
Expand Up @@ -10,7 +10,6 @@ EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude

LIB_LIBS = \
-lreactingTwoPhaseSystem \
-lcompressibleTransportModels \
-lfluidThermophysicalModels \
-lspecie
Expand Up @@ -132,7 +132,10 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const
volScalarField muAlphaI
(
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
Expand Down
Expand Up @@ -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>
Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
Expand Down
Expand Up @@ -121,10 +121,10 @@ public:
// Member Functions

//- Return the interfacial mass flow rate
virtual tmp<volScalarField> dmdt
(
const phasePairKey& key
) const;
virtual tmp<volScalarField> dmdt(const phasePairKey& key) const;

//- Return the total interfacial mass transfer rate for phase
virtual tmp<volScalarField> dmdt(const phaseModel& phase) const;

//- Return the momentum transfer matrices