ENH: Modified chtMultiRegionFoam including porous media and

secondary regions loops

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files

 fluid/compressibleCourantNo.C
 solid/solidRegionDiffNo.C
-
 chtMultiRegionFoam.C
 
-EXE = $(FOAM_APPBIN)/chtMultiRegionFoam
+EXE = $(FOAM_USER_APPBIN)/chtMultiRegionFoam
 

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options

+
 EXE_INC = \
-    -Ifluid \
-    -Isolid \
-    -Iinclude \
+    -I./fluid \
+    -I./solid \
+    -I./porousFluid \
+    -I./porousSolid \
+    -I./include \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basicSolidThermo/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude
+
 
 EXE_LIBS = \
     -lbasicThermophysicalModels \
@@ -16,7 +25,9 @@ EXE_LIBS = \
     -lspecie \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
-    -lcompressibleLESModels \
     -lmeshTools \
     -lfiniteVolume \
-    -lradiationModels
+    -lradiationModels \
+    -lfieldSources \
+    -lregionModels \
+    -lsampling

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

@@ -26,7 +26,13 @@ Application
 
 Description
     Combination of heatConductionFoam and buoyantFoam for conjugate heat
-    transfer between a solid region and fluid region
+    transfer between a solid region and fluid region. It includes
+    porous media in the primary fluid region treated explicitly.
+
+    It handles secondary fluid or solid circuits which can be coupled
+    thermally with the main fluid region. i.e radiators, etc.
+
+    The secondary fluid region is
 
 \*---------------------------------------------------------------------------*/
 
@@ -39,6 +45,9 @@ Description
 #include "solidRegionDiffNo.H"
 #include "basicSolidThermo.H"
 #include "radiationModel.H"
+#include "porousZones.H"
+#include "IObasicSourceList.H"
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,12 +60,16 @@ int main(int argc, char *argv[])
 
     #include "createFluidMeshes.H"
     #include "createSolidMeshes.H"
+    #include "createPorousFluidRegions.H"
+    #include "createPorousSolidMeshes.H"
 
     #include "createFluidFields.H"
     #include "createSolidFields.H"
+    #include "createPorousFluidFields.H"
+    #include "createPorousSolidFields.H"
 
-    #include "initContinuityErrs.H"
 
+    #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "readSolidTimeControls.H"
 
@@ -65,14 +78,12 @@ int main(int argc, char *argv[])
     #include "solidRegionDiffusionNo.H"
     #include "setInitialMultiRegionDeltaT.H"
 
-
     while (runTime.run())
     {
         #include "readTimeControls.H"
         #include "readSolidTimeControls.H"
         #include "readPIMPLEControls.H"
 
-
         #include "compressibleMultiRegionCourantNo.H"
         #include "solidRegionDiffusionNo.H"
         #include "setMultiRegionDeltaT.H"
@@ -105,6 +116,24 @@ int main(int argc, char *argv[])
                 #include "solveFluid.H"
             }
 
+            forAll(porousFluidRegions, i)
+            {
+                Info<< "\nSolving for fluid porous region "
+                    << porousFluidRegions[i].name() << endl;
+                #include "setPorousFluidFields.H"
+                #include "readPorousFluidRegionPIMPLEControls.H"
+                #include "solvePorousFluid.H"
+            }
+
+            forAll(porousSolidRegions, i)
+            {
+                Info<< "\nSolving for porous solid region "
+                    << porousSolidRegions[i].name() << endl;
+                #include "setPorousRegionSolidFields.H"
+                #include "readPorousSolidMultiRegionPIMPLEControls.H"
+                #include "solvePorousSolid.H"
+            }
+
             forAll(solidRegions, i)
             {
                 Info<< "\nSolving for solid region "
@@ -113,6 +142,7 @@ int main(int argc, char *argv[])
                 #include "readSolidMultiRegionPIMPLEControls.H"
                 #include "solveSolid.H"
             }
+
         }
 
         runTime.write();

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/Make/options

 EXE_INC = \
-    -I.. \
     -Ifluid \
     -Isolid \
+    -I./porousFluid \
+    -I./porousSolid \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basicSolidThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude
+    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
@@ -20,4 +25,7 @@ EXE_LIBS = \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
-    -lradiationModels
+    -lradiationModels \
+    -lfieldSources \
+    -lregionModels \
+    -lsampling

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C

@@ -36,6 +36,8 @@ Description
 #include "regionProperties.H"
 #include "basicSolidThermo.H"
 #include "radiationModel.H"
+#include "porousZones.H"
+#include "IObasicSourceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,9 +50,13 @@ int main(int argc, char *argv[])
 
     #include "createFluidMeshes.H"
     #include "createSolidMeshes.H"
+    #include "createPorousFluidRegions.H"
+    #include "createPorousSolidMeshes.H"
 
     #include "createFluidFields.H"
     #include "createSolidFields.H"
+    #include "createPorousFluidFields.H"
+    #include "createPorousSolidFields.H"
 
     #include "initContinuityErrs.H"
 
@@ -68,6 +74,24 @@ int main(int argc, char *argv[])
             #include "solveFluid.H"
         }
 
+        forAll(porousFluidRegions, i)
+        {
+            Info<< "\nSolving for fluid porous region "
+                << porousFluidRegions[i].name() << endl;
+            #include "setPorousFluidFields.H"
+            #include "readPorousFluidRegionSIMPLEControls.H"
+            #include "solvePorousFluid.H"
+        }
+
+        forAll(porousSolidRegions, i)
+        {
+            Info<< "\nSolving for porous solid region "
+                << porousSolidRegions[i].name() << endl;
+            #include "setPorousRegionSolidFields.H"
+            #include "readPorousSolidMultiRegionSIMPLEControls.H"
+            #include "solvePorousSolid.H"
+        }
+
         forAll(solidRegions, i)
         {
             Info<< "\nSolving for solid region "

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/UEqn.H

@@ -3,10 +3,14 @@
     (
         fvm::div(phi, U)
       + turb.divDevRhoReff(U)
+      + sources(rho, U)
     );
 
     UEqn().relax();
 
+    // Add porous resistance
+    pZones.addResistance(UEqn());
+
     solve
     (
         UEqn()

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H

@@ -18,6 +18,9 @@
     PtrList<dimensionedScalar> rhoMax(fluidRegions.size());
     PtrList<dimensionedScalar> rhoMin(fluidRegions.size());
 
+    PtrList<IObasicSourceList> heatSources(fluidRegions.size());
+    PtrList<porousZones> porousZonesFluid(fluidRegions.size());
+
     // Populate fluid field pointer lists
     forAll(fluidRegions, i)
     {
@@ -209,6 +212,20 @@
                 )
             )
         );
+
+        Info<< "    Adding sources\n" << endl;
+        heatSources.set
+        (
+            i,
+            new IObasicSourceList(fluidRegions[i])
+        );
+
+        Info<< "    Adding porous zones\n" << endl;
+        porousZonesFluid.set
+        (
+            i,
+            new porousZones(fluidRegions[i])
+        );
     }
 
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidMeshes.H

-    PtrList<fvMesh> fluidRegions(rp.fluidRegionNames().size());
+    const wordList fluidNames(rp["fluid"]);
 
-    forAll(rp.fluidRegionNames(), i)
+    PtrList<fvMesh> fluidRegions(fluidNames.size());
+
+    forAll(fluidNames, i)
     {
-        Info<< "Create fluid mesh for region " << rp.fluidRegionNames()[i]
+        Info<< "Create fluid mesh for region " << fluidNames[i]
             << " for time = " << runTime.timeName() << nl << endl;
 
         fluidRegions.set
@@ -12,7 +14,7 @@
             (
                 IOobject
                 (
-                    rp.fluidRegionNames()[i],
+                    fluidNames[i],
                     runTime.timeName(),
                     runTime,
                     IOobject::MUST_READ

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H

@@ -7,6 +7,7 @@
      ==
       - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")
       + rad.Sh(thermo)
+      + sources(rho, h)
     );
 
     hEqn.relax();

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/readFluidMultiRegionSIMPLEControls.H

@@ -3,9 +3,3 @@
     const int nNonOrthCorr =
         simple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
 
-    const bool momentumPredictor =
-        simple.lookupOrDefault("momentumPredictor", true);
-
-    const bool transonic =
-        simple.lookupOrDefault("transonic", false);
-

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H

@@ -12,6 +12,10 @@
     const volScalarField& psi = thermo.psi();
     volScalarField& h = thermo.h();
 
+    IObasicSourceList& sources = heatSources[i];
+
+    const porousZones& pZones = porousZonesFluid[i];
+
     const dimensionedScalar initialMass
     (
         "initialMass",

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/UPorousFluidEqn.H

+    // Solve the Momentum equation
+    tmp<fvVectorMatrix> porousUEqn
+    (
+        fvm::div(porousPhi, porousU)
+      + turbPorous.divDevRhoReff(porousU)
+      + porousSources(porousRho, porousU)
+    );
+
+    porousUEqn().relax();
+
+    solve(porousUEqn() == -fvc::grad(porousP));

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/createPorousFluidFields.H

+    // Initialise porous field pointer lists
+    PtrList<basicRhoThermo> thermoPorous(porousFluidRegions.size());
+    PtrList<volScalarField> rhoPorous(porousFluidRegions.size());
+    PtrList<volScalarField> kappaPorous(porousFluidRegions.size());
+    PtrList<volVectorField> UPorous(porousFluidRegions.size());
+    PtrList<surfaceScalarField> phiPorous(porousFluidRegions.size());
+    PtrList<compressible::turbulenceModel> turbulencePorous
+    (
+        porousFluidRegions.size()
+    );
+    PtrList<volScalarField> pPorous(porousFluidRegions.size());
+
+    List<scalar> initialMassFluidPorous(porousFluidRegions.size());
+    List<label> pRefCellFluidPorous(porousFluidRegions.size(),0);
+    List<scalar> pRefValueFluidPorous(porousFluidRegions.size(),0.0);
+
+    PtrList<dimensionedScalar> rhoMaxPorous(fluidRegions.size());
+    PtrList<dimensionedScalar> rhoMinPorous(fluidRegions.size());
+
+    PtrList<IObasicSourceList> heatPorousSources
+    (
+        porousFluidRegions.size()
+    );
+
+    forAll(porousFluidRegions, i)
+    {
+        Info<< "Reading fluid mesh thermophysical properties for porous "
+            << porousFluidRegions[i].name() << nl << endl;
+
+        Info<< "    Adding to thermoFluid porous\n" << endl;
+
+        thermoPorous.set
+        (
+            i,
+            basicRhoThermo::New(porousFluidRegions[i]).ptr()
+        );
+
+        Info<< "    Adding to rhoPorous\n" << endl;
+        rhoPorous.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "rho",
+                    runTime.timeName(),
+                    porousFluidRegions[i],
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                thermoPorous[i].rho()
+            )
+        );
+
+        Info<< "    Adding to UPorous\n" << endl;
+        UPorous.set
+        (
+            i,
+            new volVectorField
+            (
+                IOobject
+                (
+                    "U",
+                    runTime.timeName(),
+                    porousFluidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                porousFluidRegions[i]
+            )
+        );
+
+        Info<< "    Adding to phiPorous\n" << endl;
+        phiPorous.set
+        (
+            i,
+            new surfaceScalarField
+            (
+                IOobject
+                (
+                    "phi",
+                    runTime.timeName(),
+                    porousFluidRegions[i],
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                linearInterpolate(rhoPorous[i]*UPorous[i])
+                    & porousFluidRegions[i].Sf()
+            )
+        );
+
+        Info<< "    Adding turbulence to porous\n" << endl;
+        turbulencePorous.set
+        (
+            i,
+            compressible::turbulenceModel::New
+            (
+                rhoPorous[i],
+                UPorous[i],
+                phiPorous[i],
+                thermoPorous[i]
+            ).ptr()
+        );
+
+        Info<< "    Adding to kappaFluid\n" << endl;
+        kappaPorous.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "kappaPorous",
+                    runTime.timeName(),
+                    porousFluidRegions[i],
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                thermoPorous[i].Cp()*thermoPorous[i].alpha()
+            )
+        );
+
+        pPorous.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "p",
+                    runTime.timeName(),
+                    porousFluidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                porousFluidRegions[i]
+            )
+        );
+
+        setRefCell
+        (
+            thermoPorous[i].p(),
+            pPorous[i],
+            porousFluidRegions[i].solutionDict().subDict("SIMPLE"),
+            pRefCellFluidPorous[i],
+            pRefValueFluidPorous[i]
+        );
+
+        rhoMaxPorous.set
+        (
+            i,
+            new dimensionedScalar
+            (
+                porousFluidRegions[i].solutionDict().subDict("SIMPLE").lookup
+                (
+                    "rhoMax"
+                )
+            )
+        );
+
+        rhoMinPorous.set
+        (
+            i,
+            new dimensionedScalar
+            (
+                porousFluidRegions[i].solutionDict().subDict("SIMPLE").lookup
+                (
+                    "rhoMin"
+                )
+            )
+        );
+
+        heatPorousSources.set
+        (
+            i,
+            new IObasicSourceList(porousFluidRegions[i])
+        );
+    }
+
+

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/createPorousFluidRegions.H

+    const wordList porousFluidNames(rp["porousFluid"]);
+
+    PtrList<fvMesh> porousFluidRegions(porousFluidNames.size());
+
+    forAll (porousFluidNames, iPorous)
+    {
+        const word porousFluidName = porousFluidNames[iPorous];
+        Info<< "Create porous fluid region " << porousFluidName
+            << nl << endl;
+
+        porousFluidRegions.set
+        (
+            iPorous,
+            new fvMesh
+            (
+                IOobject
+                (
+                    porousFluidName,
+                    runTime.timeName(),
+                    runTime,
+                    IOobject::MUST_READ
+                )
+            )
+        );
+    }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/hPorousFluidEqn.H

+{
+    fvScalarMatrix hPorousEqn
+    (
+        fvm::div(porousPhi, porousH)
+      - fvm::Sp(fvc::div(porousPhi), porousH)
+      - fvm::laplacian(turbPorous.alphaEff(), porousH)
+     ==
+      - fvc::div(porousPhi, 0.5*magSqr(porousU), "div(phi,K)")
+      + porousSources(porousRho, porousH)
+    );
+
+    hPorousEqn.relax();
+    hPorousEqn.solve();
+
+    porousThermo.correct();
+
+    Info<< "Min/max in the porous T:"
+        << min(porousThermo.T()).value() << ' '
+        << max(porousThermo.T()).value() << endl;
+}

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/pPorousFluidEqn.H

+porousRho = porousThermo.rho();
+porousRho = max(porousRho, rhoMin);
+porousRho = min(porousRho, rhoMax);
+porousRho.relax();
+
+volScalarField rAUPorous(1.0/porousUEqn().A());
+
+porousU = rAUPorous*porousUEqn().H();
+porousUEqn.clear();
+
+bool closedVolume = false;
+
+porousPhi =
+    fvc::interpolate(porousRho)