From da429d77f5bf74b6367c8383179b9ee3f00dd2ac Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Mon, 11 Nov 2019 14:41:35 +0000 Subject: [PATCH] reactingTwoPhaseEulerFoam: Significantly improved handling of the particle pressure In order to improve stability and robustness of fluidised bed cases the semi-implicit treatment of the particle pressure (pPrime) is now applied within the time-step sub-cycling along with the phase differential flux update. This allows the simulations to be performed reliably at a significantly increased maximum Courant number (up to 5 for some cases) without introducing chequerboarding patterns in regions of low particle phase fraction which occurred with the previous algorithm. The fluidisedBed tutorial has been updated to be more representative of real bubbling bed cases and to demonstrate the new pPrime functionality. Developed in collaboration with Timo Niemi, VTT. --- .../multiphaseSystem/multiphaseSystem.C | 1 - .../twoPhaseSystem/twoPhaseSystem.C | 151 ++++++++---------- .../RAS/fluidisedBed/0/Theta.particles | 5 +- .../RAS/fluidisedBed/0/U.air | 16 +- .../RAS/fluidisedBed/0/U.particles | 15 +- .../RAS/fluidisedBed/Allrun | 4 +- .../constant/turbulenceProperties.particles | 6 +- .../RAS/fluidisedBed/system/controlDict | 12 +- .../RAS/fluidisedBed/system/decomposeParDict | 41 +++++ .../RAS/fluidisedBed/system/fvSolution | 22 +-- 10 files changed, 155 insertions(+), 118 deletions(-) create mode 100644 tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/system/decomposeParDict diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C index f9a2fecae9..4a9be34270 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C @@ -328,7 +328,6 @@ void Foam::multiphaseSystem::solveAlphas - fvm::laplacian(alphaDbyAs[phase.index()], alpha, "bounded") ); - alphaEqn.relax(); alphaEqn.solve(); phase.alphaPhiRef() += alphaEqn.flux(); diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C index 66f915ed38..9da3a82126 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C @@ -174,6 +174,14 @@ void Foam::twoPhaseSystem::solve surfaceScalarField phir("phir", phi1 - phi2); + tmp alphaPhiDbyA0; + if (implicitPhasePressure() && (rAUs.size() || rAUfs.size())) + { + alphaPhiDbyA0 = + this->DByAfs(rAUs, rAUfs)[phase1_.index()] + *fvc::snGrad(alpha1, "bounded")*mesh_.magSf(); + } + for (int acorr=0; acorr trSubDeltaT; - tmp alphaDbyA; - if (implicitPhasePressure() && (rAUs.size() || rAUfs.size())) + if (LTS && nAlphaSubCycles > 1) { - const surfaceScalarField DbyA + trSubDeltaT = + fv::localEulerDdt::localRSubDeltaT(mesh_, nAlphaSubCycles); + } + + for + ( + subCycle alphaSubCycle(alpha1, nAlphaSubCycles); + !(++alphaSubCycle).end(); + ) + { + surfaceScalarField alphaPhi1 ( - this->DByAfs(rAUs, rAUfs)[phase1_.index()] + fvc::flux + ( + phi_, + alpha1, + alphaScheme + ) + + fvc::flux + ( + -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme), + alpha1, + alpharScheme + ) ); - alphaDbyA = - fvc::interpolate(max(alpha1, scalar(0))) - *fvc::interpolate(max(alpha2, scalar(0))) - *DbyA; + phase1_.correctInflowOutflow(alphaPhi1); - alphaPhi1 += - alphaDbyA()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf(); - } - - phase1_.correctInflowOutflow(alphaPhi1); - - if (nAlphaSubCycles > 1) - { - tmp trSubDeltaT; - - if (LTS) + if (alphaPhiDbyA0.valid()) { - trSubDeltaT = - fv::localEulerDdt::localRSubDeltaT(mesh_, nAlphaSubCycles); + alphaPhi1 += + fvc::interpolate(max(alpha1, scalar(0))) + *fvc::interpolate(max(scalar(1) - alpha1, scalar(0))) + *alphaPhiDbyA0(); } - for - ( - subCycle alphaSubCycle(alpha1, nAlphaSubCycles); - !(++alphaSubCycle).end(); - ) - { - surfaceScalarField alphaPhi10(alphaPhi1); - - MULES::explicitSolve - ( - geometricOneField(), - alpha1, - phi_, - alphaPhi10, - (alphaSubCycle.index()*Sp)(), - (Su - (alphaSubCycle.index() - 1)*Sp*alpha1)(), - UniformField(phase1_.alphaMax()), - zeroField() - ); - - if (alphaSubCycle.index() == 1) - { - phase1_.alphaPhiRef() = alphaPhi10; - } - else - { - phase1_.alphaPhiRef() += alphaPhi10; - } - } - - phase1_.alphaPhiRef() /= nAlphaSubCycles; - } - else - { MULES::explicitSolve ( geometricOneField(), @@ -310,21 +279,39 @@ void Foam::twoPhaseSystem::solve zeroField() ); - phase1_.alphaPhiRef() = alphaPhi1; + if (alphaSubCycle.index() == 1) + { + phase1_.alphaPhiRef() = alphaPhi1; + } + else + { + phase1_.alphaPhiRef() += alphaPhi1; + } + + if (alphaPhiDbyA0.valid()) + { + const surfaceScalarField alphaDbyA + ( + fvc::interpolate(max(alpha1, scalar(0))) + *fvc::interpolate(max(scalar(1) - alpha1, scalar(0))) + *this->DByAfs(rAUs, rAUfs)[phase1_.index()] + ); + + fvScalarMatrix alpha1Eqn + ( + fvm::ddt(alpha1) - fvc::ddt(alpha1) + - fvm::laplacian(alphaDbyA, alpha1, "bounded") + ); + + alpha1Eqn.solve(); + + phase1_.alphaPhiRef() += alpha1Eqn.flux(); + } } - if (alphaDbyA.valid()) + if (nAlphaSubCycles > 1) { - fvScalarMatrix alpha1Eqn - ( - fvm::ddt(alpha1) - fvc::ddt(alpha1) - - fvm::laplacian(alphaDbyA(), alpha1, "bounded") - ); - - alpha1Eqn.relax(); - alpha1Eqn.solve(); - - phase1_.alphaPhiRef() += alpha1Eqn.flux(); + phase1_.alphaPhiRef() /= nAlphaSubCycles; } phase1_.alphaRhoPhiRef() = diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/Theta.particles b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/Theta.particles index a429765cbf..d077c5a98e 100644 --- a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/Theta.particles +++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/Theta.particles @@ -35,7 +35,10 @@ boundaryField walls { - type zeroGradient; + type JohnsonJacksonParticleTheta; + restitutionCoefficient 0.8; + specularityCoefficient 0.01; + value uniform 1e-4; } frontAndBackPlanes diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/U.air b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/U.air index 00251d0b98..1d6b2ddadb 100644 --- a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/U.air +++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/U.air @@ -22,17 +22,17 @@ boundaryField { inlet { - type interstitialInletVelocity; - inletVelocity uniform (0 0.25 0); - alpha alpha.air; - value $internalField; + type interstitialInletVelocity; + inletVelocity uniform (0 0.25 0); + alpha alpha.air; + value $internalField; } outlet { - type pressureInletOutletVelocity; - phi phi.air; - value $internalField; + type pressureInletOutletVelocity; + phi phi.air; + value $internalField; } walls @@ -42,7 +42,7 @@ boundaryField frontAndBackPlanes { - type empty; + type empty; } } diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/U.particles b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/U.particles index e09cf09d95..cc8efd3774 100644 --- a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/U.particles +++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/U.particles @@ -22,25 +22,26 @@ boundaryField { inlet { - type fixedValue; - value uniform (0 0 0); + type fixedValue; + value uniform (0 0 0); } outlet { - type fixedValue; - value uniform (0 0 0); + type fixedValue; + value uniform (0 0 0); } walls { - type fixedValue; - value uniform (0 0 0); + type JohnsonJacksonParticleSlip; + specularityCoefficient 0.01; + value uniform (0 0 0); } frontAndBackPlanes { - type empty; + type empty; } } diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/Allrun b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/Allrun index ed61ac5435..c6c14fb6a9 100755 --- a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/Allrun +++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/Allrun @@ -6,6 +6,8 @@ cd ${0%/*} || exit 1 # Run from this directory runApplication blockMesh runApplication setFields -runApplication $(getApplication) +runApplication decomposePar +runParallel $(getApplication) +runApplication reconstructPar #------------------------------------------------------------------------------ diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/constant/turbulenceProperties.particles b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/constant/turbulenceProperties.particles index 941085d3f0..97057ae611 100644 --- a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/constant/turbulenceProperties.particles +++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/constant/turbulenceProperties.particles @@ -36,16 +36,16 @@ RAS viscosityModel Gidaspow; conductivityModel Gidaspow; granularPressureModel Lun; - frictionalStressModel JohnsonJackson; + frictionalStressModel JohnsonJacksonSchaeffer; radialModel SinclairJackson; - JohnsonJacksonCoeffs + JohnsonJacksonSchaefferCoeffs { Fr 0.05; eta 2; p 5; phi 28.5; - alphaDeltaMin 0.05; + alphaDeltaMin 0.01; } } diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/system/controlDict b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/system/controlDict index 655b50b9eb..e2b8f9e4f1 100644 --- a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/system/controlDict +++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/system/controlDict @@ -25,15 +25,15 @@ stopAt endTime; endTime 2; -deltaT 0.0002; +deltaT 2e-4; -writeControl runTime; +writeControl adjustableRunTime; writeInterval 0.01; purgeWrite 0; -writeFormat ascii; +writeFormat binary; writePrecision 6; @@ -45,11 +45,11 @@ timePrecision 6; runTimeModifiable on; -adjustTimeStep no; +adjustTimeStep yes; -maxCo 0.9; +maxCo 2; -maxDeltaT 1e-05; +maxDeltaT 0.01; functions { diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/system/decomposeParDict b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/system/decomposeParDict new file mode 100644 index 0000000000..161c141c41 --- /dev/null +++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/system/decomposeParDict @@ -0,0 +1,41 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object decomposeParDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +numberOfSubdomains 4; + +/* + Main methods are: + 1) Geometric: "simple"; "hierarchical", with ordered sorting, e.g. xyz, yxz + 2) Scotch: "scotch", when running in serial; "ptscotch", running in parallel +*/ + +method hierarchical; + +simpleCoeffs +{ + n (1 4 1); // total must match numberOfSubdomains + delta 0.001; +} + +hierarchicalCoeffs +{ + n (1 4 1); // total must match numberOfSubdomains + delta 0.001; + order xyz; +} + + +// ************************************************************************* // diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/system/fvSolution b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/system/fvSolution index be83c38bc8..bc22f7d337 100644 --- a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/system/fvSolution +++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/system/fvSolution @@ -20,11 +20,14 @@ solvers "alpha.*" { nAlphaCorr 1; - nAlphaSubCycles 2; + nAlphaSubCycles 3; implicitPhasePressure yes; - solver smoothSolver; - smoother symGaussSeidel; + extremaCoeff 1; + + solver PBiCGStab; + preconditioner DIC; + tolerance 1e-9; relTol 0; minIter 1; @@ -60,12 +63,13 @@ solvers tolerance 1e-6; relTol 0; minIter 1; + maxIter 10; } "Theta.*" { - solver smoothSolver; - smoother symGaussSeidel; + solver PBiCGStab; + preconditioner DILU; tolerance 1e-6; relTol 0; minIter 1; @@ -73,8 +77,8 @@ solvers "(k|epsilon).*" { - solver smoothSolver; - smoother symGaussSeidel; + solver PBiCGStab; + preconditioner DILU; tolerance 1e-5; relTol 0; minIter 1; @@ -84,9 +88,9 @@ solvers PIMPLE { nOuterCorrectors 3; - nCorrectors 1; + nCorrectors 2; nNonOrthogonalCorrectors 0; - faceMomentum yes; + faceMomentum no; } relaxationFactors