chtMultiRegionFoam: Added support for reactions

chtMultiRegionFoam now supports reaction/combustion modelling in fluid regions in the same way as reactingFoam.
2025-12-08 06:57:46 +00:00 · 2017-12-13 08:37:25 +00:00
parent 94d05421d3
commit 7c237a59d0
11 changed files with 160 additions and 22 deletions
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
@ -1,4 +1,3 @@
-
 EXE_INC = \
    -I. \
    -I./fluid \
@ -12,7 +11,11 @@ EXE_INC = \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
@ -22,8 +25,12 @@ EXE_INC = \
 EXE_LIBS = \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
-    -lsolidThermo \
    -lspecie \
+    -lreactionThermophysicalModels \
+    -lsolidThermo \
+    -lchemistryModel \
+    -lODE \
+    -lcombustionModels \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
    -lmeshTools \
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
@ -34,8 +34,9 @@ Description
 \*---------------------------------------------------------------------------*/

 #include "fvCFD.H"
-#include "rhoThermo.H"
 #include "turbulentFluidThermoModel.H"
+#include "rhoReactionThermo.H"
+#include "CombustionModel.H"
 #include "fixedGradientFvPatchFields.H"
 #include "regionProperties.H"
 #include "compressibleCourantNo.H"
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/EEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/EEqn.H
@ -15,10 +15,11 @@
            )
          : -dpdt
        )
-      - fvm::laplacian(turb.alphaEff(), he)
+      - fvm::laplacian(turbulence.alphaEff(), he)
     ==
        rho*(U&g)
      + rad.Sh(thermo, he)
+      + Qdot
      + fvOptions(rho, he)
    );

--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
@ -6,7 +6,7 @@
    (
        fvm::ddt(rho, U) + fvm::div(phi, U)
      + MRF.DDt(rho, U)
-      + turb.divDevRhoReff(U)
+      + turbulence.divDevRhoReff(U)
     ==
        fvOptions(rho, U)
    );
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/YEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/YEqn.H
@ -0,0 +1,61 @@
+tmp<fv::convectionScheme<scalar>> mvConvection(nullptr);
+
+if (Y.size())
+{
+    mvConvection = tmp<fv::convectionScheme<scalar>>
+    (
+        fv::convectionScheme<scalar>::New
+        (
+            mesh,
+            fields,
+            phi,
+            mesh.divScheme("div(phi,Yi_h)")
+        )
+    );
+}
+
+{
+    reaction.correct();
+    Qdot = reaction.Qdot();
+    volScalarField Yt
+    (
+        IOobject("Yt", runTime.timeName(), mesh),
+        mesh,
+        dimensionedScalar("Yt", dimless, 0)
+    );
+
+    forAll(Y, i)
+    {
+        if (i != inertIndex && composition.active(i))
+        {
+            volScalarField& Yi = Y[i];
+
+            fvScalarMatrix YiEqn
+            (
+                fvm::ddt(rho, Yi)
+              + mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence.muEff(), Yi)
+             ==
+                reaction.R(Yi)
+              + fvOptions(rho, Yi)
+            );
+
+            YiEqn.relax();
+
+            fvOptions.constrain(YiEqn);
+
+            YiEqn.solve(mesh.solver("Yi"));
+
+            fvOptions.correct(Yi);
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+    }
+
+    if (Y.size())
+    {
+        Y[inertIndex] = scalar(1) - Yt;
+        Y[inertIndex].max(0.0);
+    }
+}
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@ -1,5 +1,5 @@
 // Initialise fluid field pointer lists
-PtrList<rhoThermo> thermoFluid(fluidRegions.size());
+PtrList<rhoReactionThermo> thermoFluid(fluidRegions.size());
 PtrList<volScalarField> rhoFluid(fluidRegions.size());
 PtrList<volVectorField> UFluid(fluidRegions.size());
 PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
@ -7,11 +7,15 @@ PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
 PtrList<uniformDimensionedScalarField> hRefFluid(fluidRegions.size());
 PtrList<volScalarField> ghFluid(fluidRegions.size());
 PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
-PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
+PtrList<compressible::turbulenceModel> turbulenceFluid(fluidRegions.size());
+PtrList<CombustionModel<rhoReactionThermo>> reactionFluid(fluidRegions.size());
 PtrList<volScalarField> p_rghFluid(fluidRegions.size());
 PtrList<radiation::radiationModel> radiation(fluidRegions.size());
 PtrList<volScalarField> KFluid(fluidRegions.size());
 PtrList<volScalarField> dpdtFluid(fluidRegions.size());
+PtrList<multivariateSurfaceInterpolationScheme<scalar>::fieldTable>
+    fieldsFluid(fluidRegions.size());
+PtrList<volScalarField> QdotFluid(fluidRegions.size());

 List<scalar> initialMassFluid(fluidRegions.size());
 List<bool> frozenFlowFluid(fluidRegions.size(), false);
@ -29,12 +33,7 @@ forAll(fluidRegions, i)
        << fluidRegions[i].name() << nl << endl;

    Info<< "    Adding to thermoFluid\n" << endl;
-
-    thermoFluid.set
-    (
-        i,
-        rhoThermo::New(fluidRegions[i]).ptr()
-    );
+    thermoFluid.set(i, rhoReactionThermo::New(fluidRegions[i]).ptr());

    Info<< "    Adding to rhoFluid\n" << endl;
    rhoFluid.set
@ -156,8 +155,8 @@ forAll(fluidRegions, i)
        )
    );

-    Info<< "    Adding to turbulence\n" << endl;
-    turbulence.set
+    Info<< "    Adding to turbulenceFluid\n" << endl;
+    turbulenceFluid.set
    (
        i,
        compressible::turbulenceModel::New
@ -169,6 +168,17 @@ forAll(fluidRegions, i)
        ).ptr()
    );

+    Info<< "    Adding to reactionFluid\n" << endl;
+    reactionFluid.set
+    (
+        i,
+        CombustionModel<rhoReactionThermo>::New
+        (
+            thermoFluid[i],
+            turbulenceFluid[i]
+        )
+    );
+
    p_rghFluid.set
    (
        i,
@ -191,6 +201,7 @@ forAll(fluidRegions, i)

    fluidRegions[i].setFluxRequired(p_rghFluid[i].name());

+    Info<< "    Adding to radiationFluid\n" << endl;
    radiation.set
    (
        i,
@ -232,6 +243,37 @@ forAll(fluidRegions, i)
        )
    );

+    Info<< "    Adding to fieldsFluid\n" << endl;
+    fieldsFluid.set
+    (
+        i,
+        new multivariateSurfaceInterpolationScheme<scalar>::fieldTable
+    );
+    forAll(thermoFluid[i].composition().Y(), j)
+    {
+        fieldsFluid[i].add(thermoFluid[i].composition().Y()[j]);
+    }
+    fieldsFluid[i].add(thermoFluid[i].he());
+
+    Info<< "    Adding to QdotFluid\n" << endl;
+    QdotFluid.set
+    (
+        i,
+        new volScalarField
+        (
+            IOobject
+            (
+                "Qdot",
+                runTime.timeName(),
+                fluidRegions[i],
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            fluidRegions[i],
+            dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
+        )
+    );
+
    const dictionary& pimpleDict =
        fluidRegions[i].solutionDict().subDict("PIMPLE");
    pimpleDict.readIfPresent("frozenFlow", frozenFlowFluid[i]);
@ -250,7 +292,7 @@ forAll(fluidRegions, i)
        new fv::options(fluidRegions[i])
    );

-    turbulence[i].validate();
+    turbulenceFluid[i].validate();

    if (pimpleDict.isDict("residualControl"))
    {
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@ -1,13 +1,33 @@
    fvMesh& mesh = fluidRegions[i];

-    rhoThermo& thermo = thermoFluid[i];
+    CombustionModel<rhoReactionThermo>& reaction = reactionFluid[i];
+
+    rhoReactionThermo& thermo = reaction.thermo();
    thermo.validate(args.executable(), "h", "e");

+    basicMultiComponentMixture& composition = thermo.composition();
+    PtrList<volScalarField>& Y = composition.Y();
+
+    label inertIndex = -1;
+    if (Y.size())
+    {
+        const word inertSpecie(thermo.lookup("inertSpecie"));
+        if (!composition.species().found(inertSpecie))
+        {
+            FatalIOErrorIn(args.executable().c_str(), thermo)
+                << "Inert specie " << inertSpecie
+                << " not found in available species "
+                << composition.species()
+                << exit(FatalIOError);
+        }
+        inertIndex = composition.species()[inertSpecie];
+    }
+
    volScalarField& rho = rhoFluid[i];
    volVectorField& U = UFluid[i];
    surfaceScalarField& phi = phiFluid[i];

-    compressible::turbulenceModel& turb = turbulence[i];
+    compressible::turbulenceModel& turbulence = turbulenceFluid[i];
    volScalarField& K = KFluid[i];
    volScalarField& dpdt = dpdtFluid[i];

@ -20,6 +40,11 @@
    const volScalarField& gh = ghFluid[i];
    const surfaceScalarField& ghf = ghfFluid[i];

+    multivariateSurfaceInterpolationScheme<scalar>::fieldTable& fields =
+        fieldsFluid[i];
+
+    volScalarField& Qdot = QdotFluid[i];
+
    radiation::radiationModel& rad = radiation[i];

    IOMRFZoneList& MRF = MRFfluid[i];
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
@ -15,6 +15,7 @@ else
    }

    #include "UEqn.H"
+    #include "YEqn.H"
    #include "EEqn.H"

    // --- PISO loop
@ -23,7 +24,7 @@ else
        #include "pEqn.H"
    }

-    turb.correct();
+    turbulence.correct();

    rho = thermo.rho();
 }
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/constant/bottomWater/thermophysicalProperties
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/constant/bottomWater/thermophysicalProperties
@ -17,7 +17,7 @@ FoamFile
 thermoType
 {
    type            heRhoThermo;
-    mixture         pureMixture;
+    mixture         singleComponentMixture;
    transport       const;
    thermo          hConst;
    equationOfState rhoConst;
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/constant/topAir/thermophysicalProperties
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/constant/topAir/thermophysicalProperties
@ -18,7 +18,7 @@ FoamFile
 thermoType
 {
    type            heRhoThermo;
-    mixture         pureMixture;
+    mixture         singleComponentMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/constant/bottomAir/thermophysicalProperties
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/constant/bottomAir/thermophysicalProperties
@ -18,7 +18,7 @@ FoamFile
 thermoType
 {
    type            heRhoThermo;
-    mixture         pureMixture;
+    mixture         singleComponentMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;