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.

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

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C

@@ -174,6 +174,14 @@ void Foam::twoPhaseSystem::solve
 
     surfaceScalarField phir("phir", phi1 - phi2);
 
+    tmp<surfaceScalarField> 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<nAlphaCorr; acorr++)
     {
         volScalarField::Internal Sp
@@ -219,6 +227,20 @@ void Foam::twoPhaseSystem::solve
             }
         }
 
+        tmp<volScalarField> trSubDeltaT;
+
+        if (LTS && nAlphaSubCycles > 1)
+        {
+            trSubDeltaT =
+                fv::localEulerDdt::localRSubDeltaT(mesh_, nAlphaSubCycles);
+        }
+
+        for
+        (
+            subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+            !(++alphaSubCycle).end();
+        )
+        {
             surfaceScalarField alphaPhi1
             (
                 fvc::flux
@@ -235,69 +257,16 @@ void Foam::twoPhaseSystem::solve
                 )
             );
 
-        tmp<surfaceScalarField> alphaDbyA;
-        if (implicitPhasePressure() && (rAUs.size() || rAUfs.size()))
-        {
-            const surfaceScalarField DbyA
-            (
-                this->DByAfs(rAUs, rAUfs)[phase1_.index()]
-            );
-
-            alphaDbyA =
-                fvc::interpolate(max(alpha1, scalar(0)))
-               *fvc::interpolate(max(alpha2, scalar(0)))
-               *DbyA;
-
-            alphaPhi1 +=
-                alphaDbyA()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf();
-        }
-
             phase1_.correctInflowOutflow(alphaPhi1);
 
-        if (nAlphaSubCycles > 1)
+            if (alphaPhiDbyA0.valid())
             {
-            tmp<volScalarField> trSubDeltaT;
+                alphaPhi1 +=
-
+                    fvc::interpolate(max(alpha1, scalar(0)))
-            if (LTS)
+                   *fvc::interpolate(max(scalar(1) - alpha1, scalar(0)))
-            {
+                   *alphaPhiDbyA0();
-                trSubDeltaT =
-                    fv::localEulerDdt::localRSubDeltaT(mesh_, nAlphaSubCycles);
             }
 
-            for
-            (
-                subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
-                !(++alphaSubCycle).end();
-            )
-            {
-                surfaceScalarField alphaPhi10(alphaPhi1);
-
-                MULES::explicitSolve
-                (
-                    geometricOneField(),
-                    alpha1,
-                    phi_,
-                    alphaPhi10,
-                    (alphaSubCycle.index()*Sp)(),
-                    (Su - (alphaSubCycle.index() - 1)*Sp*alpha1)(),
-                    UniformField<scalar>(phase1_.alphaMax()),
-                    zeroField()
-                );
-
-                if (alphaSubCycle.index() == 1)
-                {
-                    phase1_.alphaPhiRef() = alphaPhi10;
-                }
-                else
-                {
-                    phase1_.alphaPhiRef() += alphaPhi10;
-                }
-            }
-
-            phase1_.alphaPhiRef() /= nAlphaSubCycles;
-        }
-        else
-        {
             MULES::explicitSolve
             (
                 geometricOneField(),
@@ -310,22 +279,40 @@ void Foam::twoPhaseSystem::solve
                 zeroField()
             );
 
+            if (alphaSubCycle.index() == 1)
+            {
                 phase1_.alphaPhiRef() = alphaPhi1;
             }
-
+            else
-        if (alphaDbyA.valid())
             {
+                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")
+                  - fvm::laplacian(alphaDbyA, alpha1, "bounded")
                 );
 
-            alpha1Eqn.relax();
                 alpha1Eqn.solve();
 
                 phase1_.alphaPhiRef() += alpha1Eqn.flux();
             }
+        }
+
+        if (nAlphaSubCycles > 1)
+        {
+            phase1_.alphaPhiRef() /= nAlphaSubCycles;
+        }
 
         phase1_.alphaRhoPhiRef() =
             fvc::interpolate(phase1_.rho())*phase1_.alphaPhi();

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

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/U.particles

@@ -34,7 +34,8 @@ boundaryField
 
     walls
     {
-        type               fixedValue;
+        type            JohnsonJacksonParticleSlip;
+        specularityCoefficient  0.01;
         value           uniform (0 0 0);
     }
 

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

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

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
 {

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

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/system/fvSolution

@@ -20,11 +20,14 @@ solvers
     "alpha.*"
     {
         nAlphaCorr      1;
-        nAlphaSubCycles 2;
+        nAlphaSubCycles 3;
         implicitPhasePressure yes;
 
-        solver          smoothSolver;
+        extremaCoeff    1;
-        smoother        symGaussSeidel;
+
+        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;
+        solver          PBiCGStab;
-        smoother        symGaussSeidel;
+        preconditioner  DILU;
         tolerance       1e-6;
         relTol          0;
         minIter         1;
@@ -73,8 +77,8 @@ solvers
 
     "(k|epsilon).*"
     {
-        solver          smoothSolver;
+        solver          PBiCGStab;
-        smoother        symGaussSeidel;
+        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