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