multiphaseEuler: nAlphaSubCycles is now a Function1 of the alpha Courant number

Henry Weller · Henry Weller · commit 1706ec99bbf8 · 2025-02-16T11:46:43.000Z
The explicit MULES sub-cycling control nAlphaSubCycles is now a Function1 of the
alpha Courant number which allows a very convenient way to specify the number of
sub-cycles independent of the time-step the case is run with while maintaining
stability of the explicit solution of the phase-fraction.  For example with the
following entry in fvSolution of the pipeBend case:

solvers
{
    "alpha.*"
    {
        nAlphaSubCycles polynomial (0 3);
    }
.
.
.

nAlphaSubCycles is evaluated as three times the current maximum alpha Courant
number rounded upwards to an integer, i.e. given the maximum Courant number is
set to 1 for this case 3 sub-cycles will be executed.
diff --git a/applications/modules/multiphaseEuler/multiphaseEuler.C b/applications/modules/multiphaseEuler/multiphaseEuler.C
@@ -61,6 +61,8 @@ bool Foam::solvers::multiphaseEuler::read()
     nEnergyCorrectors =
         pimple.dict().lookupOrDefault<int>("nEnergyCorrectors", 1);
 
+    alphaControls.read(mesh.solution().solverDict("alpha"));
+
     return true;
 }
 
@@ -247,7 +249,8 @@ void Foam::solvers::multiphaseEuler::prePredictor()
 {
     if (pimple.thermophysics() || pimple.flow())
     {
-        fluid_.solve(rAs);
+        alphaControls.correct(CoNum);
+        fluid_.solve(alphaControls, rAs);
         fluid_.correct();
         fluid_.correctContinuityError();
     }
diff --git a/applications/modules/multiphaseEuler/multiphaseEuler.H b/applications/modules/multiphaseEuler/multiphaseEuler.H
@@ -94,6 +94,8 @@ protected:
         //  Defaults to 1
         int nEnergyCorrectors;
 
+        phaseSystem::alphaControl alphaControls;
+
 
     //- Optional LTS reciprocal time-step field
     tmp<volScalarField> trDeltaT;
diff --git a/applications/modules/multiphaseEuler/phaseSystem/phaseSystem/phaseSystem.H b/applications/modules/multiphaseEuler/phaseSystem/phaseSystem/phaseSystem.H
@@ -74,6 +74,46 @@ class phaseSystem
 {
 public:
 
+    // Public Types
+
+        struct alphaControl
+        {
+            //- Function to calculate the number of explicit MULES sub-cycles
+            //  from the alpha Courant number
+            autoPtr<Function1<scalar>> nAlphaSubCyclesPtr;
+
+            label nAlphaSubCycles;
+
+            //- Number of alpha correctors
+            //  Usually used to improve the accuracy at high Courant numbers
+            //  with semi-implicit MULES, MULESCorr = true
+            label nAlphaCorr;
+
+            scalar vDotResidualAlpha;
+
+            void correct(const scalar CoNum)
+            {
+                nAlphaSubCycles = ceil(nAlphaSubCyclesPtr->value(CoNum));
+            }
+
+            void read(const dictionary& dict)
+            {
+                nAlphaSubCyclesPtr = Function1<scalar>::New
+                (
+                    "nAlphaSubCycles",
+                    dimless,
+                    dimless,
+                    dict
+                );
+
+                nAlphaCorr = dict.lookupOrDefault<label>("nAlphaCorr", 1);
+
+                vDotResidualAlpha =
+                    dict.lookupOrDefault("vDotResidualAlpha", 1e-4);
+            }
+        };
+
+
     // Public Typedefs
 
         typedef HashPtrTable<fvVectorMatrix> momentumTransferTable;
@@ -528,7 +568,11 @@ public:
         // Evolution
 
             //- Solve for the phase fractions
-            virtual void solve(const PtrList<volScalarField>& rAs);
+            virtual void solve
+            (
+                const alphaControl& alphaControls,
+                const PtrList<volScalarField>& rAs
+            );
 
             //- Correct the fluid properties other than those listed below
             virtual void correct();
diff --git a/applications/modules/multiphaseEuler/phaseSystem/phaseSystem/phaseSystemSolve.C b/applications/modules/multiphaseEuler/phaseSystem/phaseSystem/phaseSystemSolve.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2024 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2025 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -41,20 +41,15 @@ License
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::phaseSystem::solve(const PtrList<volScalarField>& rAs)
+void Foam::phaseSystem::solve
+(
+    const alphaControl& alphaControls,
+    const PtrList<volScalarField>& rAs
+)
 {
-    const dictionary& alphaControls = mesh_.solution().solverDict("alpha");
-
-    const label nAlphaSubCycles(alphaControls.lookup<label>("nAlphaSubCycles"));
-    const label nAlphaCorr(alphaControls.lookup<label>("nAlphaCorr"));
-
+    const label nAlphaSubCycles = alphaControls.nAlphaSubCycles;
     const bool LTS = fv::localEulerDdt::enabled(mesh_);
 
-    const scalar vDotResidualAlpha
-    (
-        alphaControls.lookupOrDefault("vDotResidualAlpha", 1e-4)
-    );
-
     // Optional reference phase which is not solved for
     // but obtained from the sum of the other phases
     phaseModel* referencePhasePtr = nullptr;
@@ -122,7 +117,7 @@ void Foam::phaseSystem::solve(const PtrList<volScalarField>& rAs)
         }
     }
 
-    for (int acorr=0; acorr<nAlphaCorr; acorr++)
+    for (int acorr=0; acorr<alphaControls.nAlphaCorr; acorr++)
     {
         PtrList<volScalarField::Internal> Sps(phases().size());
         PtrList<volScalarField::Internal> Sus(phases().size());
@@ -201,16 +196,28 @@ void Foam::phaseSystem::solve(const PtrList<volScalarField>& rAs)
                     {
                         Sp[celli] -=
                             vDot[celli]
-                           /max(1 - alpha[celli], vDotResidualAlpha);
+                           /max
+                            (
+                                1 - alpha[celli],
+                                alphaControls.vDotResidualAlpha
+                            );
                         Su[celli] +=
                             vDot[celli]
-                           /max(1 - alpha[celli], vDotResidualAlpha);
+                           /max
+                            (
+                                1 - alpha[celli],
+                                alphaControls.vDotResidualAlpha
+                            );
                     }
                     else if (vDot[celli] < 0)
                     {
                         Sp[celli] +=
                             vDot[celli]
-                           /max(alpha[celli], vDotResidualAlpha);
+                           /max
+                            (
+                                alpha[celli],
+                                alphaControls.vDotResidualAlpha
+                            );
                     }
                 }
             }
diff --git a/applications/modules/multiphaseEuler/phaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C b/applications/modules/multiphaseEuler/phaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2025 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -175,10 +175,11 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::specieTransfer() const
 template<class BasePhaseSystem>
 void Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::solve
 (
+    const phaseSystem::alphaControl& alphaControls,
     const PtrList<volScalarField>& rAs
 )
 {
-    BasePhaseSystem::solve(rAs);
+    BasePhaseSystem::solve(alphaControls, rAs);
 
     forAll(populationBalances_, i)
     {
diff --git a/applications/modules/multiphaseEuler/phaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H b/applications/modules/multiphaseEuler/phaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2025 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -97,7 +97,11 @@ public:
             specieTransfer() const;
 
         //- Solve all population balance equations
-        virtual void solve(const PtrList<volScalarField>& rAs);
+        virtual void solve
+        (
+            const phaseSystem::alphaControl& alphaControls,
+            const PtrList<volScalarField>& rAs
+        );
 
         //- Correct derived properties
         virtual void correct();
diff --git a/tutorials/multiphaseEuler/bubbleColumn/system/fvSolution b/tutorials/multiphaseEuler/bubbleColumn/system/fvSolution
@@ -18,7 +18,6 @@ solvers
 {
     "alpha.*"
     {
-        nAlphaCorr      1;
         nAlphaSubCycles 2;
     }
 
diff --git a/tutorials/multiphaseEuler/pipeBend/system/fvSolution b/tutorials/multiphaseEuler/pipeBend/system/fvSolution
@@ -18,7 +18,7 @@ solvers
     "alpha.*"
     {
         nAlphaCorr      1;
-        nAlphaSubCycles 3;
+        nAlphaSubCycles polynomial (0 3);
         implicitPhasePressure yes;
         solver          PBiCGStab;
         preconditioner  DIC;

Original file line number	Diff line number	Diff line change
`@@ -61,6 +61,8 @@ bool Foam::solvers::multiphaseEuler::read()`
`61`	`61`	`nEnergyCorrectors =`
`62`	`62`	`pimple.dict().lookupOrDefault<int>("nEnergyCorrectors", 1);`
`63`	`63`
	`64`	`+ alphaControls.read(mesh.solution().solverDict("alpha"));`
	`65`	`+`
`64`	`66`	`return true;`
`65`	`67`	`}`
`66`	`68`
`@@ -247,7 +249,8 @@ void Foam::solvers::multiphaseEuler::prePredictor()`
`247`	`249`	`{`
`248`	`250`	`if (pimple.thermophysics() \|\| pimple.flow())`
`249`	`251`	`{`
`250`		`- fluid_.solve(rAs);`
	`252`	`+ alphaControls.correct(CoNum);`
	`253`	`+ fluid_.solve(alphaControls, rAs);`
`251`	`254`	`fluid_.correct();`
`252`	`255`	`fluid_.correctContinuityError();`
`253`	`256`	`}`
Original file line number	Diff line number	Diff line change
`@@ -2,7 +2,7 @@`
`2`	`2`	`========= \|`
`3`	`3`	`\\ / F ield \| OpenFOAM: The Open Source CFD Toolbox`
`4`	`4`	`\\ / O peration \| Website: https://openfoam.org`
`5`		`- \\ / A nd \| Copyright (C) 2017-2023 OpenFOAM Foundation`
	`5`	`+ \\ / A nd \| Copyright (C) 2017-2025 OpenFOAM Foundation`
`6`	`6`	`\\/ M anipulation \|`
`7`	`7`	`-------------------------------------------------------------------------------`
`8`	`8`	`License`
`@@ -175,10 +175,11 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::specieTransfer() const`
`175`	`175`	`template<class BasePhaseSystem>`
`176`	`176`	`void Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::solve`
`177`	`177`	`(`
	`178`	`+ const phaseSystem::alphaControl& alphaControls,`
`178`	`179`	`const PtrList<volScalarField>& rAs`
`179`	`180`	`)`
`180`	`181`	`{`
`181`		`- BasePhaseSystem::solve(rAs);`
	`182`	`+ BasePhaseSystem::solve(alphaControls, rAs);`
`182`	`183`
`183`	`184`	`forAll(populationBalances_, i)`
`184`	`185`	`{`
Original file line number	Diff line number	Diff line change
`@@ -18,7 +18,6 @@ solvers`
`18`	`18`	`{`
`19`	`19`	`"alpha.*"`
`20`	`20`	`{`
`21`		`- nAlphaCorr 1;`
`22`	`21`	`nAlphaSubCycles 2;`
`23`	`22`	`}`
`24`	`23`
Original file line number	Diff line number	Diff line change
`@@ -18,7 +18,7 @@ solvers`
`18`	`18`	`"alpha.*"`
`19`	`19`	`{`
`20`	`20`	`nAlphaCorr 1;`
`21`		`- nAlphaSubCycles 3;`
	`21`	`+ nAlphaSubCycles polynomial (0 3);`
`22`	`22`	`implicitPhasePressure yes;`
`23`	`23`	`solver PBiCGStab;`
`24`	`24`	`preconditioner DIC;`