zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source 
to support any number of phases
This commit is contained in:
Henry Weller
2015-06-26 15:15:10 +01:00
parent 0f94c22141
commit f67c9ab9bd
11 changed files with 99 additions and 48 deletions
Download Patch File Download Diff File Expand all files Collapse all files

23
applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H
View File

@ -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);
}
}

5
applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H
View File

@ -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

2
applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C
View File

@ -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_;

2
applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H
View File

@ -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);

25
applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C
View File

@ -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

14
applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H
View File

@ -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

19
applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
View File

@ -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;
}
// ************************************************************************* //

15
applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
View File

@ -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

30
applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C
View File

@ -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 =
surfaceScalarField DbyA(phase1_.DbyA() + phase2_.DbyA());
alphaDbyA =
fvc::interpolate(max(alpha1, scalar(0)))
*fvc::interpolate(max(alpha2, scalar(0)));
*fvc::interpolate(max(alpha2, scalar(0)))
*DbyA;
surfaceScalarField phiP
(
pPrimeByA_()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf()
);
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();

6
applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.H
View File

@ -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();
};

6
applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystemI.H
View File

@ -65,10 +65,4 @@ inline const Foam::phaseModel& Foam::twoPhaseSystem::otherPhase
}
inline Foam::tmp<Foam::surfaceScalarField>& Foam::twoPhaseSystem::pPrimeByA()
{
return pPrimeByA_;
}
// ************************************************************************* //