XiFluid: Extended to support rich inhomogeneous mixtures and exhaust gas recirculation

The implementation is based on the work undertaken by Weller in 2002 to extend
the Weller b-Xi combustion model to non-premixed combustion:

    Weller, H. G. (2002, August).
    The Application of the Weller Combustion Models to
    Non-Premixed Combustion.
    (Technical Report TR/HGW/03)

Class
    Foam::solvers::XiFluid

Description
    Solver module for compressible premixed/partially-premixed combustion with
    turbulence modelling.

    Combusting RANS code using the Weller b-Xi two-equation combustion model.
    Xi may be obtained by either the solution of the Xi transport equation or
    from an algebraic expression.

    Reference:
    \verbatim
        Weller, H. G. (1993).
        The development of a new flame area combustion model
        using conditional averaging.
        Thermo-fluids section report TF 9307.
    \endverbatim

    Both approaches are based on Gulder's flame speed correlation which has been
    shown to be appropriate by comparison with the results from the spectral
    model.

    Reference:
    \verbatim
        Weller, H. G., Marooney, C. J., & Gosman, A. D. (1991, January).
        A new spectral method for calculation of the time-varying area
        of a laminar flame in homogeneous turbulence.
        In Symposium (International) on Combustion
        (Vol. 23, No. 1, pp. 629-636). Elsevier.
    \endverbatim

    Strain effects are incorporated directly into the Xi equation
    but not in the algebraic approximation.  Further work need to be
    done on this issue, particularly regarding the enhanced removal rate
    caused by flame compression.  Analysis using results of the spectral
    model will be required.

    For cases involving very lean Propane flames or other flames which are
    very strain-sensitive, a transport equation for the laminar flame
    speed is present.  This equation is derived using heuristic arguments
    involving the strain time scale and the strain-rate at extinction.
    the transport velocity is the same as that for the Xi equation.

    Reference:
    \verbatim
        Weller, H. G., Tabor, G., Gosman, A. D., & Fureby, C. (1998, January).
        Application of a flame-wrinkling LES combustion model
        to a turbulent mixing layer.
        In Symposium (International) on combustion
        (Vol. 27, No. 1, pp. 899-907). Elsevier.
    \endverbatim

    For inhomogeneous mixtures, in addition to the regress variable \c b, it is
    necessary to solve for the mixture-fraction \c ft provided by the \c
    leanInhomogeneousMixture and also the fuel mass-fraction \fu if there are
    rich regions in the mixture, provided by the \c inhomogeneousMixture.
    Details of the extention of the Weller b-Xi combustion model to non-premixed
    combustion can be found in the Technical Report TR/HGW/03.

    Reference:
    \verbatim
        Weller, H. G. (2002, August).
        The Application of the Weller Combustion Models to
        Non-Premixed Combustion.
        (Technical Report TR/HGW/03)
    \endverbatim

    For inhomogeneous mixtures with exhaust gas re-circulation it is necessary
    to additionally solve for the recirculated exhaust gas mass-fraction \c
    egr which is provided by the \c inhomogeneousEGRMixture mixture.

    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
    pseudo-transient and steady simulations.

    Reference:
    \verbatim
        Greenshields, C. J., & Weller, H. G. (2022).
        Notes on Computational Fluid Dynamics: General Principles.
        CFD Direct Ltd.: Reading, UK.
    \endverbatim

    Optional fvModels and fvConstraints are provided to enhance the simulation
    in many ways including adding various sources, chemical reactions,
    combustion, Lagrangian particles, radiation, surface film etc. and
    constraining or limiting the solution.

Author: Henry Weller

applications/modules/XiFluid/XiFluid.C

@@ -105,6 +105,16 @@ Foam::solvers::XiFluid::XiFluid(fvMesh& mesh)
         fields.add(thermo_.Y("ft"));
     }
 
+    if (thermo_.containsSpecie("fu"))
+    {
+        fields.add(thermo_.Y("fu"));
+    }
+
+    if (thermo_.containsSpecie("egr"))
+    {
+        fields.add(thermo_.Y("egr"));
+    }
+
     fields.add(b);
     fields.add(thermo.he());
     fields.add(thermo.heu());

applications/modules/XiFluid/XiFluid.H

@@ -28,9 +28,9 @@ Description
     Solver module for compressible premixed/partially-premixed combustion with
     turbulence modelling.
 
-    Combusting RANS code using the b-Xi two-equation model.
-    Xi may be obtained by either the solution of the Xi transport
-    equation or from an algebraic expression.
+    Combusting RANS code using the Weller b-Xi two-equation combustion model.
+    Xi may be obtained by either the solution of the Xi transport equation or
+    from an algebraic expression.
 
     Reference:
     \verbatim
@@ -74,6 +74,25 @@ Description
         (Vol. 27, No. 1, pp. 899-907). Elsevier.
     \endverbatim
 
+    For inhomogeneous mixtures, in addition to the regress variable \c b, it is
+    necessary to solve for the mixture-fraction \c ft provided by the \c
+    leanInhomogeneousMixture and also the fuel mass-fraction \fu if there are
+    rich regions in the mixture, provided by the \c inhomogeneousMixture.
+    Details of the extention of the Weller b-Xi combustion model to non-premixed
+    combustion can be found in the Technical Report TR/HGW/03.
+
+    Reference:
+    \verbatim
+        Weller, H. G. (2002, August).
+        The Application of the Weller Combustion Models to
+        Non-Premixed Combustion.
+        (Technical Report TR/HGW/03)
+    \endverbatim
+
+    For inhomogeneous mixtures with exhaust gas re-circulation it is necessary
+    to additionally solve for the recirculated exhaust gas mass-fraction \c
+    egr which is provided by the \c inhomogeneousEGRMixture mixture.
+
     Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
     pseudo-transient and steady simulations.
 
@@ -176,6 +195,22 @@ protected:
             const volScalarField& Db
         );
 
+        //- Solve the fu equation for partially- and non- premixed mixtures
+        void fuSolve
+        (
+            const fv::convectionScheme<scalar>& mvConvection,
+            const volScalarField& Db,
+            const volScalarField& bSource
+        );
+
+        //- Solve the egr mass-fraction equation
+        //  for mixtures with exhaust gas recirculation
+        void egrSolve
+        (
+            const fv::convectionScheme<scalar>& mvConvection,
+            const volScalarField& Db
+        );
+
         //- Apply the early kernel growth correction to the flame-wrinkling Xi
         //  to compensate for the inaccurate flame surface area estimate
         //  when the kernel is not fully-developed,

applications/modules/XiFluid/thermophysicalPredictor.C

@@ -27,6 +27,7 @@ License
 #include "bXiIgnition.H"
 #include "fvcDdt.H"
 #include "fvmDiv.H"
+#include "fvcSup.H"
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
@@ -54,6 +55,75 @@ void Foam::solvers::XiFluid::ftSolve
 }
 
 
+void Foam::solvers::XiFluid::fuSolve
+(
+    const fv::convectionScheme<scalar>& mvConvection,
+    const volScalarField& Db,
+    const volScalarField& bSource
+)
+{
+    volScalarField& fu = thermo_.Y("fu");
+    const volScalarField& ft = thermo_.Y("ft");
+    const volScalarField& b(b_);
+
+    // Progress variable
+    const volScalarField c("c", scalar(1) - b);
+
+    // Unburnt gas density
+    const volScalarField rhou("rhou", thermo.rhou());
+
+    const volScalarField fres(thermo_.fres());
+    const volScalarField fuFres(max(fu - fres, scalar(0)));
+
+    const volScalarField fuDot =
+        bSource/(b + 0.001)
+      - 2*Db*c
+       *mag(fvc::grad(ft))/max(ft, 1e-6)
+       *mag(fvc::grad(fuFres))/max(fuFres, 1e-6);
+
+    fvScalarMatrix fuEqn
+    (
+        fvm::ddt(rho, fu)
+      + mvConvection.fvmDiv(phi, fu)
+      - fvm::laplacian(Db, fu)
+     ==
+        fvm::Sp(fuDot, fu)
+      - fuDot*fres
+      + fvModels().source(rho, fu)
+    );
+
+    fuEqn.relax();
+    fvConstraints().constrain(fuEqn);
+    fuEqn.solve();
+    fvConstraints().constrain(fu);
+}
+
+
+void Foam::solvers::XiFluid::egrSolve
+(
+    const fv::convectionScheme<scalar>& mvConvection,
+    const volScalarField& Db
+)
+{
+    volScalarField& egr = thermo_.Y("egr");
+
+    fvScalarMatrix egrEqn
+    (
+        fvm::ddt(rho, egr)
+      + mvConvection.fvmDiv(phi, egr)
+      - fvm::laplacian(Db, egr)
+     ==
+        fvModels().source(rho, egr)
+    );
+
+    egrEqn.relax();
+    fvConstraints().constrain(egrEqn);
+    egrEqn.solve();
+    fvConstraints().constrain(egr);
+}
+
+
+
 Foam::tmp<Foam::volScalarField> Foam::solvers::XiFluid::XiCorr
 (
     const volScalarField& Xi,
@@ -148,16 +218,21 @@ void Foam::solvers::XiFluid::bSolve
         fvc::interpolate(rhou*Su*XiCorr(Xi, nf, dMgb))*nf
     );
 
+    fvScalarMatrix bSourceEqn
+    (
+        fvModels().source(rho, b)
+      - fvm::div(phiSt, b)
+      + fvm::Sp(fvc::div(phiSt), b)
+    );
+
     // Create b equation
     fvScalarMatrix bEqn
     (
         fvm::ddt(rho, b)
       + mvConvection.fvmDiv(phi, b)
-      + fvm::div(phiSt, b)
-      - fvm::Sp(fvc::div(phiSt), b)
       - fvm::laplacian(Db, b)
      ==
-        fvModels().source(rho, b)
+        bSourceEqn
     );
 
     // Solve for b and constrain
@@ -171,6 +246,11 @@ void Foam::solvers::XiFluid::bSolve
 
     // Correct the laminar flame speed
     SuModel_->correct();
+
+    if (thermo_.containsSpecie("fu"))
+    {
+        fuSolve(mvConvection, Db, (bSourceEqn & b));
+    }
 }
 
 
@@ -295,6 +375,20 @@ void Foam::solvers::XiFluid::thermophysicalPredictor()
         bSolve(mvConvection(), Db);
         EauSolve(mvConvection(), Db);
     }
+    else
+    {
+        if (thermo_.containsSpecie("fu"))
+        {
+            volScalarField& fu = thermo_.Y("fu");
+            const volScalarField& ft = thermo_.Y("ft");
+            fu = ft;
+        }
+    }
+
+    if (thermo_.containsSpecie("egr"))
+    {
+        egrSolve(mvConvection(), Db);
+    }
 
     EaSolve(mvConvection(), Db);
 

src/OpenFOAM/db/functionObjects/timeControl/timeControl.C

@@ -148,7 +148,7 @@ void Foam::timeControl::read(const dictionary& dict)
         {
             const word timesName(prefix_ + "Times");
             const word frequenciesName(prefix_ + "Frequencies");
-            const bool repeat = dict.lookupOrDefault("writeRepeat", false);
+            const bool repeat = dict.lookupOrDefault(prefix_ + "Repeat", false);
 
             timeDelta_ =
                 dict.lookupOrDefault

tutorials/XiFluid/stratified/0/T

@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    left
+    {
+        type            symmetryPlane;
+    }
+
+    right
+    {
+        type            symmetryPlane;
+    }
+
+    top
+    {
+        type            symmetryPlane;
+    }
+
+    bottom
+    {
+        type            symmetryPlane;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/XiFluid/stratified/0/Tu

@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    format      ascii;
+    class       volScalarField;
+    object      Tu;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    left
+    {
+        type            symmetryPlane;
+    }
+
+    right
+    {
+        type            symmetryPlane;
+    }
+
+    top
+    {
+        type            symmetryPlane;
+    }
+
+    bottom
+    {
+        type            symmetryPlane;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/XiFluid/stratified/0/U

@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    left
+    {
+        type            symmetryPlane;
+    }
+
+    right
+    {
+        type            symmetryPlane;
+    }
+
+    top
+    {
+        type            symmetryPlane;
+    }
+
+    bottom
+    {
+        type            symmetryPlane;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //