From 2b74f9486f2a9e41f49f2286d0f2e486d2fb319a Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Fri, 7 Apr 2023 14:16:21 +0100 Subject: [PATCH] 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. --- .../lagrangian/denseParticleFoam/UcEqn.H | 23 +--- .../denseParticleFoam/denseParticleFoam.C | 114 +++++------------- .../lagrangian/denseParticleFoam/pEqn.H | 56 ++++++--- 3 files changed, 72 insertions(+), 121 deletions(-) diff --git a/applications/solvers/lagrangian/denseParticleFoam/UcEqn.H b/applications/solvers/lagrangian/denseParticleFoam/UcEqn.H index 11556c8188..d1bdcfe0c1 100644 --- a/applications/solvers/lagrangian/denseParticleFoam/UcEqn.H +++ b/applications/solvers/lagrangian/denseParticleFoam/UcEqn.H @@ -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); diff --git a/applications/solvers/lagrangian/denseParticleFoam/denseParticleFoam.C b/applications/solvers/lagrangian/denseParticleFoam/denseParticleFoam.C index 883a0d11ee..00f796597e 100644 --- a/applications/solvers/lagrangian/denseParticleFoam/denseParticleFoam.C +++ b/applications/solvers/lagrangian/denseParticleFoam/denseParticleFoam.C @@ -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::names[] = - { - "faceExplicitCellImplicit", - "faceExplicitCellLagged", - "faceImplicit" - }; - - const NamedEnum 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()) { diff --git a/applications/solvers/lagrangian/denseParticleFoam/pEqn.H b/applications/solvers/lagrangian/denseParticleFoam/pEqn.H index 8ecee4039b..2e7dbf9c94 100644 --- a/applications/solvers/lagrangian/denseParticleFoam/pEqn.H +++ b/applications/solvers/lagrangian/denseParticleFoam/pEqn.H @@ -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);