updates + added porous zones

applications/solvers/Lagrangian/trackedReactingParcelFoam/UEqn.H

@@ -10,6 +10,37 @@
 
     UEqn.relax();
 
+    tmp<volScalarField> trAU;
+    tmp<volTensorField> trTU;
+
+    if (pressureImplicitPorosity)
+    {
+        tmp<volTensorField> tTU = tensor(I)*UEqn.A();
+        pZones.addResistance(UEqn, tTU());
+        trTU = inv(tTU());
+        trTU().rename("rAU");
+
+        volVectorField gradp = fvc::grad(p);
+
+        for (int UCorr=0; UCorr<nUCorr; UCorr++)
+        {
+            U = trTU() & (UEqn.H() - gradp);
+        }
+        U.correctBoundaryConditions();
+    }
+    else
+    {
+        pZones.addResistance(UEqn);
+
+        solve
+        (
+            UEqn == -fvc::grad(p)
+        );
+
+        trAU = 1.0/UEqn.A();
+        trAU().rename("rAU");
+    }
+
     if (momentumPredictor)
     {
         solve(UEqn == -fvc::grad(p));

applications/solvers/Lagrangian/trackedReactingParcelFoam/createFields.H

-    Info<< "Reading thermophysical properties\n" << endl;
+    Info<< "Reading thermophysical properties" << nl << endl;
 
     autoPtr<hCombustionThermo> thermo
     (
@@ -28,52 +28,7 @@
         thermo->rho()
     );
 
-// lagrangian coal density field
-/*    volScalarField rholc
-    (
-        IOobject
-        (
-            "rholc",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0, 0, 0), 0.0)
-    );
-
-// lagrangian limestone density field
-    volScalarField rhols
-    (
-        IOobject
-        (
-            "rhols",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0, 0, 0), 0.0)
-    );
-
-// lagrangian total density field
-    volScalarField rhol
-    (
-        IOobject
-        (
-            "rhol",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0, 0, 0), 0.0)
-    );*/
-
-    Info<< "\nReading field U\n" << endl;
+    Info<< "Reading field U" << nl << endl;
     volVectorField U
     (
         IOobject
@@ -103,7 +58,7 @@
         dimensionedScalar("zero", dimless, 0.0)
     );
 
-    Info<< "Creating turbulence model\n" << endl;
+    Info<< "Creating turbulence model" << nl << endl;
     autoPtr<compressible::turbulenceModel> turbulence
     (
         compressible::turbulenceModel::New
@@ -115,11 +70,11 @@
         )
     );
 
-    Info<< "Creating field DpDt\n" << endl;
+    Info<< "Creating field DpDt" << nl << endl;
     volScalarField DpDt =
         fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
-    Info << "Constructing chemical mechanism" << endl;
+    Info << "Constructing chemical mechanism" << nl << endl;
     chemistryModel chemistry
     (
         thermo(),
@@ -134,8 +89,23 @@
     }
     fields.add(h);
 
-    Info<< "Creating radiation model\n" << endl;
-    autoPtr<radiation::radiationModel> radiation
-    (
-        radiation::radiationModel::New(T)
-    );
+    Info<< "Creating porous zones" << nl << endl;
+    porousZones pZones(mesh);
+    Switch pressureImplicitPorosity(false);
+
+    label nUCorr = 0;
+    if (pZones.size())
+    {
+        // nUCorrectors for pressureImplicitPorosity
+        if (mesh.solutionDict().subDict("PISO").found("nUCorrectors"))
+        {
+            mesh.solutionDict().subDict("PISO").lookup("nUCorrectors")
+                >> nUCorr;
+        }
+
+        if (nUCorr > 0)
+        {
+            pressureImplicitPorosity = true;
+        }
+    }
+

applications/solvers/Lagrangian/trackedReactingParcelFoam/pEqn.H

 rho = thermo->rho();
 
-volScalarField rUA = 1.0/UEqn.A();
-U = rUA*UEqn.H();
+if (pressureImplicitPorosity)
+{
+    U = trTU()&UEqn.H();
+}
+else
+{
+    U = trAU()*UEqn.H();
+}
 
 if (transonic)
 {
@@ -11,17 +17,28 @@ if (transonic)
         fvc::interpolate(thermo->psi())
        *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+//          + fvc::ddtPhiCorr(rUA, rho, U, phi)
         )
     );
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
+        tmp<fvScalarMatrix> lapTerm;
+
+        if (pressureImplicitPorosity)
+        {
+            lapTerm = fvm::laplacian(rho*trTU(), p);
+        }
+        else
+        {
+            lapTerm = fvm::laplacian(rho*trAU(), p);
+        }
+
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rUA, p)
+          - lapTerm()
          ==
             reactingParcels.Srho()
         );
@@ -37,19 +54,30 @@ if (transonic)
 else
 {
     phi =
-        fvc::interpolate(rho)*
-        (
+        fvc::interpolate(rho)
+       *(
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+          + fvc::ddtPhiCorr(trAU(), rho, U, phi)
         );
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
+        tmp<fvScalarMatrix> lapTerm;
+
+        if (pressureImplicitPorosity)
+        {
+            lapTerm = fvm::laplacian(rho*trTU(), p);
+        }
+        else
+        {
+            lapTerm = fvm::laplacian(rho*trAU(), p);
+        }
+
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
           + fvc::div(phi)
-          - fvm::laplacian(rho*rUA, p)
+          - lapTerm()
          ==
             reactingParcels.Srho()
         );
@@ -66,7 +94,15 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rUA*fvc::grad(p);
+if (pressureImplicitPorosity)
+{
+    U -= trTU()&fvc::grad(p);
+}
+else
+{
+    U -= trAU()*fvc::grad(p);
+}
 U.correctBoundaryConditions();
 
 DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+

applications/solvers/Lagrangian/trackedReactingParcelFoam/trackedReactingParcelFoam.C

@@ -36,60 +36,60 @@ Description
 #include "chemistrySolver.H"
 #include "ReactingCloudThermoTypes.H"
 #include "radiationModel.H"
+#include "porousZones.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
-#   include "setRootCase.H"
-
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readChemistryProperties.H"
-#   include "readEnvironmentalProperties.H"
-#   include "createFields.H"
-#   include "createClouds.H"
-#   include "readPISOControls.H"
-#   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
-#   include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readChemistryProperties.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "createRadiationModel.H"
+    #include "createClouds.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readTimeControls.H"
-#       include "readPISOControls.H"
-#       include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "readTimeControls.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        Info << "Evolving reacting cloud" << endl;
-
         reactingParcels.evolve();
 
         reactingParcels.info();
 
-#       include "chemistry.H"
-#       include "rhoEqn.H"
+        #include "chemistry.H"
+        #include "rhoEqn.H"
 
         // --- PIMPLE loop
         for (int ocorr=1; ocorr<=nOuterCorr; ocorr++)
         {
-#           include "UEqn.H"
-#           include "YEqn.H"
+            #include "UEqn.H"
+            #include "YEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
-#               include "hEqn.H"
-#               include "pEqn.H"
+                #include "hEqn.H"
+                #include "pEqn.H"
             }
 
             Info<< "T gas min/max   = " << min(T).value() << ", "
@@ -102,7 +102,7 @@ int main(int argc, char *argv[])
 
         if (runTime.write())
         {
-#           include "additionalOutput.H"
+            #include "additionalOutput.H"
         }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"