Commit bf95b5c2 authored by Henry Weller's avatar Henry Weller
Browse files

reactingTwoPhaseEulerFoam: Change the implicit handling of phase-pressure and dispersion

to support any number of phases
parent 4d6823c3
......@@ -39,7 +39,7 @@ tmp<surfaceScalarField> phiF2;
volScalarField D(fluid.D());
// Phase-1 turbulent dispersion and particle-pressure flux
surfaceScalarField Df1
tmp<surfaceScalarField> DbyA1
(
fvc::interpolate
(
......@@ -48,7 +48,7 @@ tmp<surfaceScalarField> phiF2;
);
// Phase-2 turbulent dispersion and particle-pressure flux
surfaceScalarField Df2
tmp<surfaceScalarField> DbyA2
(
fvc::interpolate
(
......@@ -56,13 +56,6 @@ tmp<surfaceScalarField> phiF2;
)
);
// Cache the net diffusivity for implicit diffusion treatment in the
// phase-fraction equation
if (implicitPhasePressure)
{
fluid.pPrimeByA() = Df1 + Df2;
}
// Lift and wall-lubrication forces
volVectorField F(fluid.F());
......@@ -72,16 +65,24 @@ tmp<surfaceScalarField> phiF2;
// Phase-1 dispersion, lift and wall-lubrication flux
phiF1 =
(
Df1*snGradAlpha1
DbyA1()*snGradAlpha1
+ (fvc::interpolate(rAU1*F) & mesh.Sf())
);
// Phase-2 dispersion, lift and wall-lubrication flux
phiF2 =
(
- Df2*snGradAlpha1
- DbyA2()*snGradAlpha1
- (fvc::interpolate(rAU2*F) & mesh.Sf())
);
// Cache the phase diffusivities for implicit treatment in the
// phase-fraction equation
if (implicitPhasePressure)
{
phase1.DbyA(DbyA1);
phase2.DbyA(DbyA2);
}
}
......
......@@ -69,11 +69,12 @@ tmp<surfaceScalarField> Ff2;
fvc::interpolate(D + phase2.turbulence().pPrime())
);
// Cache the net diffusivity for implicit diffusion treatment in the
// Cache the phase diffusivities for implicit treatment in the
// phase-fraction equation
if (implicitPhasePressure)
{
fluid.pPrimeByA() = rAUf1*Df1 + rAUf2*Df2;
phase1.DbyA(rAUf1*Df1);
phase2.DbyA(rAUf2*Df2);
}
// Lift and wall-lubrication forces
......
......@@ -139,7 +139,7 @@ bool Foam::AnisothermalPhaseModel<BasePhaseModel>::compressible() const
template<class BasePhaseModel>
Foam::tmp<Foam::volScalarField>
const Foam::volScalarField&
Foam::AnisothermalPhaseModel<BasePhaseModel>::divU() const
{
return divU_;
......
......@@ -90,7 +90,7 @@ public:
virtual bool compressible() const;
//- Phase dilatation rate (d(alpha)/dt + div(alpha*phi))
virtual tmp<volScalarField> divU() const;
virtual const volScalarField& divU() const;
//- Set phase dilatation rate (d(alpha)/dt + div(alpha*phi))
virtual void divU(const volScalarField& divU);
......
......@@ -257,6 +257,31 @@ Foam::MovingPhaseModel<BasePhaseModel>::UEqn()
}
template<class BasePhaseModel>
const Foam::surfaceScalarField&
Foam::MovingPhaseModel<BasePhaseModel>::DbyA() const
{
if (DbyA_.valid())
{
return DbyA_;
}
else
{
return surfaceScalarField::null();
}
}
template<class BasePhaseModel>
void Foam::MovingPhaseModel<BasePhaseModel>::DbyA
(
const tmp<surfaceScalarField>& DbyA
)
{
DbyA_ = DbyA;
}
template<class BasePhaseModel>
Foam::tmp<Foam::volVectorField>
Foam::MovingPhaseModel<BasePhaseModel>::U() const
......
......@@ -88,6 +88,10 @@ class MovingPhaseModel
//- Continuity error
volScalarField continuityError_;
//- Phase diffusivity divided by the momentum coefficient.
// Used for implicit treatment of the phase pressure and dispersion
tmp<surfaceScalarField> DbyA_;
// Private static member functions
......@@ -124,6 +128,16 @@ public:
//- Return the momentum equation
virtual tmp<fvVectorMatrix> UEqn();
// Implicit phase pressure and dispersion support
//- Return the phase diffusivity divided by the momentum coefficient
virtual const surfaceScalarField& DbyA() const;
//- Set the phase diffusivity divided by the momentum coefficient
virtual void DbyA(const tmp<surfaceScalarField>& DbyA);
// Momentum
//- Constant access the velocity
......
......@@ -140,9 +140,10 @@ bool Foam::phaseModel::compressible() const
}
Foam::tmp<Foam::volScalarField> Foam::phaseModel::divU() const
const Foam::volScalarField& Foam::phaseModel::divU() const
{
return tmp<volScalarField>();
notImplemented("Foam::phaseModel::divU()");
return volScalarField::null();
}
......@@ -154,4 +155,18 @@ void Foam::phaseModel::divU(const volScalarField& divU)
}
const Foam::surfaceScalarField& Foam::phaseModel::DbyA() const
{
return surfaceScalarField::null();
}
void Foam::phaseModel::DbyA(const tmp<surfaceScalarField>& DbyA)
{
WarningIn("phaseModel::DbyA(const surfaceScalarField& DbyA)")
<< "Attempt to set the dilatation rate of an incompressible phase"
<< endl;
}
// ************************************************************************* //
......@@ -169,13 +169,22 @@ public:
//- Return true if the phase is compressible otherwise false
virtual bool compressible() const;
//- Phase dilatation rate (d(alpha)/dt + div(alpha*phi))
virtual tmp<volScalarField> divU() const;
//- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
virtual const volScalarField& divU() const;
//- Set phase dilatation rate (d(alpha)/dt + div(alpha*phi))
//- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
virtual void divU(const volScalarField& divU);
// Implicit phase pressure and dispersion support
//- Return the phase diffusivity divided by the momentum coefficient
virtual const surfaceScalarField& DbyA() const;
//- Set the phase diffusivity divided by the momentum coefficient
virtual void DbyA(const tmp<surfaceScalarField>& DbyA);
// Thermo
//- Return const access to the thermophysical model
......
......@@ -190,9 +190,9 @@ void Foam::twoPhaseSystem::solve()
tdgdt =
(
alpha2.dimensionedInternalField()
*phase1_.divU()().dimensionedInternalField()
*phase1_.divU().dimensionedInternalField()
- alpha1.dimensionedInternalField()
*phase2_.divU()().dimensionedInternalField()
*phase2_.divU().dimensionedInternalField()
);
}
else if (phase1_.compressible())
......@@ -200,7 +200,7 @@ void Foam::twoPhaseSystem::solve()
tdgdt =
(
alpha2.dimensionedInternalField()
*phase1_.divU()().dimensionedInternalField()
*phase1_.divU().dimensionedInternalField()
);
}
else if (phase2_.compressible())
......@@ -208,7 +208,7 @@ void Foam::twoPhaseSystem::solve()
tdgdt =
(
- alpha1.dimensionedInternalField()
*phase2_.divU()().dimensionedInternalField()
*phase2_.divU().dimensionedInternalField()
);
}
......@@ -218,20 +218,18 @@ void Foam::twoPhaseSystem::solve()
surfaceScalarField phic("phic", phi);
surfaceScalarField phir("phir", phi1 - phi2);
tmp<surfaceScalarField> alpha1alpha2f;
tmp<surfaceScalarField> alphaDbyA;
if (pPrimeByA_.valid())
if (notNull(phase1_.DbyA()) && notNull(phase2_.DbyA()))
{
alpha1alpha2f =
fvc::interpolate(max(alpha1, scalar(0)))
*fvc::interpolate(max(alpha2, scalar(0)));
surfaceScalarField DbyA(phase1_.DbyA() + phase2_.DbyA());
surfaceScalarField phiP
(
pPrimeByA_()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf()
);
alphaDbyA =
fvc::interpolate(max(alpha1, scalar(0)))
*fvc::interpolate(max(alpha2, scalar(0)))
*DbyA;
phir += phiP;
phir += DbyA*fvc::snGrad(alpha1, "bounded")*mesh_.magSf();
}
for (int acorr=0; acorr<nAlphaCorr; acorr++)
......@@ -367,12 +365,12 @@ void Foam::twoPhaseSystem::solve()
phase1_.alphaPhi() = alphaPhic1;
}
if (pPrimeByA_.valid())
if (alphaDbyA.valid())
{
fvScalarMatrix alpha1Eqn
(
fvm::ddt(alpha1) - fvc::ddt(alpha1)
- fvm::laplacian(alpha1alpha2f*pPrimeByA_(), alpha1, "bounded")
- fvm::laplacian(alphaDbyA, alpha1, "bounded")
);
alpha1Eqn.relax();
......
......@@ -63,9 +63,6 @@ protected:
//- Phase model 2
phaseModel& phase2_;
//- Optional dispersion diffusivity
tmp<surfaceScalarField> pPrimeByA_;
public:
......@@ -173,9 +170,6 @@ public:
//- Return the virtual mass model for the given phase
const virtualMassModel& virtualMass(const phaseModel& phase) const;
//- Return non-const access to the dispersion diffusivity
inline tmp<surfaceScalarField>& pPrimeByA();
};
......
......@@ -65,10 +65,4 @@ inline const Foam::phaseModel& Foam::twoPhaseSystem::otherPhase
}
inline Foam::tmp<Foam::surfaceScalarField>& Foam::twoPhaseSystem::pPrimeByA()
{
return pPrimeByA_;
}
// ************************************************************************* //
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment