denseParticleFoam: New face-stabilised phase drag formulation
to provide consistent and stable continuous phase velocity solution without staggering patterns at the boundary with packed regions of dispersed phase.
This commit is contained in:
Henry Weller
@ -3,28 +3,12 @@ fvVectorMatrix UcEqn
|
||||
fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
|
||||
- fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
|
||||
+ continuousPhaseTurbulence->divDevTau(Uc)
|
||||
==
|
||||
(1.0/rhoc)*cloudSU
|
||||
);
|
||||
|
||||
UcEqn.relax();
|
||||
|
||||
fvConstraints.constrain(UcEqn);
|
||||
|
||||
volScalarField rAUc(1.0/UcEqn.A());
|
||||
volScalarField rASpUc(1.0/(UcEqn.A() - cloudSUp/rhoc));
|
||||
surfaceScalarField rASpUcf("Dp", fvc::interpolate(rASpUc));
|
||||
|
||||
surfaceScalarField phicSUSu
|
||||
(
|
||||
fvc::flux(rASpUc*cloudSUu/rhoc)
|
||||
+ rASpUcf*(g & mesh.Sf())
|
||||
);
|
||||
surfaceScalarField phicSUSp
|
||||
(
|
||||
fvc::interpolate(rASpUc*cloudSUp/rhoc)
|
||||
);
|
||||
|
||||
if (pimple.momentumPredictor())
|
||||
{
|
||||
solve
|
||||
@ -33,10 +17,11 @@ if (pimple.momentumPredictor())
|
||||
==
|
||||
fvc::reconstruct
|
||||
(
|
||||
(phicSUSu + phicSUSp*phic)/rASpUcf
|
||||
- fvc::snGrad(p)*mesh.magSf()
|
||||
Fgf - fvc::snGrad(p)*mesh.magSf()
|
||||
- Dcf*(phic - phid)
|
||||
)
|
||||
+ (1.0/rhoc)*(fvm::Sp(cloudSUp, Uc) - cloudSUp*Uc)
|
||||
+ Dc*fvc::reconstruct(phic - phid)
|
||||
+ Fd - fvm::Sp(Dc, Uc)
|
||||
);
|
||||
|
||||
fvConstraints.constrain(Uc);
|
||||
|
||||
@ -31,50 +31,6 @@ Description
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "NamedEnum.H"
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
enum class cloudForceSplit
|
||||
{
|
||||
faceExplicitCellImplicit, // Implicit part of the cloud force added to
|
||||
// the cell momentum equation. Explicit part
|
||||
// to the face momentum equation. This is the
|
||||
// least likely to create staggering patterns
|
||||
// in the velocity field, but it can create
|
||||
// unphysical perturbations in cell
|
||||
// velocities even when particles and flow
|
||||
// have the similar velocities.
|
||||
|
||||
faceExplicitCellLagged, // Entire cloud force evaluated explicitly
|
||||
// and added to the face momentum equation.
|
||||
// Lagged correction (i.e.,
|
||||
// fvm::Sp(cloudSU.diag(), Uc) -
|
||||
// cloudSU.diag()*Uc) added to the cell
|
||||
// momentum equation. This creates physical
|
||||
// cell velocities when particles and flow
|
||||
// have the same velocity, but can also
|
||||
// result in staggering patterns in packed
|
||||
// beds. Unsuitable for MPPIC.
|
||||
|
||||
faceImplicit // Implicit and explicit parts of the force
|
||||
// both added to the face momentum equation.
|
||||
// Behaves somewhere between the other two.
|
||||
};
|
||||
|
||||
template<>
|
||||
const char* NamedEnum<cloudForceSplit, 3>::names[] =
|
||||
{
|
||||
"faceExplicitCellImplicit",
|
||||
"faceExplicitCellLagged",
|
||||
"faceImplicit"
|
||||
};
|
||||
|
||||
const NamedEnum<cloudForceSplit, 3> cloudForceSplitNames;
|
||||
}
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#include "argList.H"
|
||||
#include "timeSelector.H"
|
||||
#include "viscosityModel.H"
|
||||
@ -165,69 +121,61 @@ int main(int argc, char *argv[])
|
||||
alphacf = fvc::interpolate(alphac);
|
||||
alphaPhic = alphacf*phic;
|
||||
|
||||
// Cloud forces
|
||||
fvVectorMatrix cloudSU(clouds.SU(Uc));
|
||||
volVectorField cloudSUu
|
||||
// Dispersed phase drag force
|
||||
volVectorField Fd
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"cloudSUu",
|
||||
"Fd",
|
||||
runTime.name(),
|
||||
mesh
|
||||
),
|
||||
mesh,
|
||||
dimensionedVector(dimForce/dimVolume, Zero),
|
||||
dimensionedVector(dimAcceleration, Zero),
|
||||
zeroGradientFvPatchVectorField::typeName
|
||||
);
|
||||
volScalarField cloudSUp
|
||||
|
||||
// continuous-dispersed phase drag coefficient
|
||||
volScalarField Dc
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"cloudSUp",
|
||||
"Dc",
|
||||
runTime.name(),
|
||||
mesh
|
||||
),
|
||||
mesh,
|
||||
dimensionedScalar(dimForce/dimVelocity/dimVolume, Zero),
|
||||
dimensionedScalar(dimless/dimTime, Zero),
|
||||
zeroGradientFvPatchVectorField::typeName
|
||||
);
|
||||
|
||||
const cloudForceSplit cloudSUSplit =
|
||||
pimple.dict().found("cloudForceSplit")
|
||||
? cloudForceSplitNames.read(pimple.dict().lookup("cloudForceSplit"))
|
||||
: cloudForceSplit::faceExplicitCellImplicit;
|
||||
|
||||
switch (cloudSUSplit)
|
||||
{
|
||||
case cloudForceSplit::faceExplicitCellImplicit:
|
||||
cloudSUu.primitiveFieldRef() = -cloudSU.source()/mesh.V();
|
||||
cloudSUu.correctBoundaryConditions();
|
||||
cloudSUp.primitiveFieldRef() = Zero;
|
||||
cloudSUp.correctBoundaryConditions();
|
||||
//cloudSU.diag() = cloudSU.diag();
|
||||
cloudSU.source() = Zero;
|
||||
break;
|
||||
const fvVectorMatrix cloudSU(clouds.SU(Uc));
|
||||
|
||||
case cloudForceSplit::faceExplicitCellLagged:
|
||||
cloudSUu.primitiveFieldRef() =
|
||||
(cloudSU.diag()*Uc() - cloudSU.source())/mesh.V();
|
||||
cloudSUu.correctBoundaryConditions();
|
||||
cloudSUp.primitiveFieldRef() = Zero;
|
||||
cloudSUp.correctBoundaryConditions();
|
||||
//cloudSU.diag() = cloudSU.diag();
|
||||
cloudSU.source() = cloudSU.diag()*Uc();
|
||||
break;
|
||||
// Dispersed phase drag force
|
||||
Fd.primitiveFieldRef() = -cloudSU.source()/mesh.V()/rhoc;
|
||||
Fd.correctBoundaryConditions();
|
||||
|
||||
case cloudForceSplit::faceImplicit:
|
||||
cloudSUu.primitiveFieldRef() = -cloudSU.source()/mesh.V();
|
||||
cloudSUu.correctBoundaryConditions();
|
||||
cloudSUp.primitiveFieldRef() = cloudSU.diag()/mesh.V();
|
||||
cloudSUp.correctBoundaryConditions();
|
||||
cloudSU.diag() = Zero;
|
||||
cloudSU.source() = Zero;
|
||||
break;
|
||||
// Continuous phase drag coefficient
|
||||
Dc.primitiveFieldRef() = -cloudSU.diag()/mesh.V()/rhoc;
|
||||
Dc.correctBoundaryConditions();
|
||||
}
|
||||
|
||||
// Face continuous-dispersed phase drag coefficient
|
||||
const surfaceScalarField Dcf(fvc::interpolate(Dc));
|
||||
|
||||
// Face dispersed phase drag force
|
||||
const surfaceScalarField Fdf(fvc::flux(Fd));
|
||||
|
||||
// Effective flux of the dispersed phase
|
||||
const surfaceScalarField phid
|
||||
(
|
||||
Fdf/(Dcf + dimensionedScalar(Dc.dimensions(), small))
|
||||
);
|
||||
|
||||
// Face buoyancy force
|
||||
const surfaceScalarField Fgf(g & mesh.Sf());
|
||||
|
||||
// --- Pressure-velocity PIMPLE corrector loop
|
||||
while (pimple.loop())
|
||||
{
|
||||
|
||||
@ -1,37 +1,46 @@
|
||||
{
|
||||
volVectorField HbyA(constrainHbyA(rAUc*UcEqn.H(), Uc, p));
|
||||
volVectorField HbyASp(rASpUc/rAUc*HbyA);
|
||||
const volScalarField rAUc(1.0/UcEqn.A());
|
||||
const volScalarField r1ADUc(1/(1 + rAUc*Dc));
|
||||
|
||||
surfaceScalarField phiHbyASp
|
||||
const surfaceScalarField rAUcf(fvc::interpolate(rAUc));
|
||||
const surfaceScalarField r1ADUcf(1/(1 + rAUcf*fvc::interpolate(Dc)));
|
||||
const surfaceScalarField rADUcf("Dp", r1ADUcf*rAUcf);
|
||||
|
||||
volVectorField HbyA(constrainHbyA(rAUc*UcEqn.H(), Uc, p));
|
||||
|
||||
surfaceScalarField phiHbyAD
|
||||
(
|
||||
"phiHbyASp",
|
||||
"phiHbyAD",
|
||||
(
|
||||
fvc::flux(HbyASp)
|
||||
+ alphacf*rASpUcf*fvc::ddtCorr(Uc, phic, Ucf)
|
||||
r1ADUcf
|
||||
*(
|
||||
fvc::flux(HbyA)
|
||||
+ alphacf*rAUcf*fvc::ddtCorr(Uc, phic, Ucf)
|
||||
)
|
||||
)
|
||||
);
|
||||
|
||||
if (p.needReference())
|
||||
{
|
||||
fvc::makeRelative(phiHbyASp, Uc);
|
||||
adjustPhi(phiHbyASp, Uc, p);
|
||||
fvc::makeAbsolute(phiHbyASp, Uc);
|
||||
fvc::makeRelative(phiHbyAD, Uc);
|
||||
adjustPhi(phiHbyAD, Uc, p);
|
||||
fvc::makeAbsolute(phiHbyAD, Uc);
|
||||
}
|
||||
|
||||
phiHbyASp += phicSUSu;
|
||||
phiHbyAD += rADUcf*(Fgf + Fdf);
|
||||
|
||||
// Update the pressure BCs to ensure flux consistency
|
||||
constrainPressure(p, Uc, phiHbyASp, rASpUcf);
|
||||
constrainPressure(p, Uc, phiHbyAD, rADUcf);
|
||||
|
||||
// Non-orthogonal pressure corrector loop
|
||||
while (pimple.correctNonOrthogonal())
|
||||
{
|
||||
fvScalarMatrix pEqn
|
||||
(
|
||||
fvm::laplacian(alphacf*rASpUcf, p)
|
||||
fvm::laplacian(alphacf*rADUcf, p)
|
||||
==
|
||||
fvc::ddt(alphac)
|
||||
+ fvc::div(alphacf*phiHbyASp)
|
||||
+ fvc::div(alphacf*phiHbyAD)
|
||||
);
|
||||
|
||||
pEqn.setReference
|
||||
@ -44,18 +53,27 @@
|
||||
|
||||
if (pimple.finalNonOrthogonalIter())
|
||||
{
|
||||
phic = phiHbyASp - pEqn.flux()/alphacf;
|
||||
phic = phiHbyAD - pEqn.flux()/alphacf;
|
||||
|
||||
// Explicitly relax pressure for momentum corrector
|
||||
p.relax();
|
||||
|
||||
Uc =
|
||||
HbyA
|
||||
+ rAUc
|
||||
*fvc::reconstruct
|
||||
(
|
||||
(phicSUSu + phicSUSp*phic - pEqn.flux()/alphacf)/rASpUcf
|
||||
r1ADUc
|
||||
*(
|
||||
HbyA
|
||||
+ rAUc
|
||||
*(
|
||||
fvc::reconstruct
|
||||
(
|
||||
Fgf - pEqn.flux()/alphacf/rADUcf
|
||||
- Dcf*(phic - phid)
|
||||
)
|
||||
+ Dc*fvc::reconstruct(phic - phid)
|
||||
+ Fd
|
||||
)
|
||||
);
|
||||
|
||||
Uc.correctBoundaryConditions();
|
||||
fvConstraints.constrain(Uc);
|
||||
|
||||
|
||||
Reference in New Issue
Block a user