ENH: Added film+pyrolysis solver

applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/Make/files

@@ -0,0 +1,3 @@
+reactingParcelFilmPyrolysisFoam.C
+
+EXE = $(FOAM_APPBIN)/reactingParcelFilmPyrolysisFoam

applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/Make/options

@@ -0,0 +1,51 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solid/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basicSolidThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solidChemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/combustionModels/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/regionModels/pyrolysisModels/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -lspecie \
+    -lbasicThermophysicalModels \
+    -lsolidProperties \
+    -lsolidMixtureProperties \
+    -lthermophysicalFunctions \
+    -lreactionThermophysicalModels \
+    -lSLGThermo \
+    -lchemistryModel \
+    -lsolidChemistryModel \
+    -lcombustionModels \
+    -lregionModels \
+    -lradiationModels \
+    -lsurfaceFilmModels \
+    -lpyrolysisModels \
+    -llagrangianIntermediate \
+    -lODE \
+    -lsampling

applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/UEqn.H

@@ -0,0 +1,26 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        parcels.SU(U)
+    );
+
+    UEqn.relax();
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+          ==
+            fvc::reconstruct
+            (
+                (
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                )*mesh.magSf()
+            )
+        );
+    }

applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/YhsEqn.H

@@ -0,0 +1,69 @@
+tmp<fv::convectionScheme<scalar> > mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_hs)")
+    )
+);
+{
+    combustion->correct();
+    dQ = combustion->dQ();
+    label inertIndex = -1;
+    volScalarField Yt = 0.0*Y[0];
+
+    forAll(Y, i)
+    {
+        if (Y[i].name() != inertSpecie)
+        {
+            volScalarField& Yi = Y[i];
+            fvScalarMatrix R = combustion->R(Yi);
+
+            solve
+            (
+                fvm::ddt(rho, Yi)
+              + mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence->alphaEff(), Yi)
+              ==
+                parcels.SYi(i, Yi)
+              + surfaceFilm.Srho(i)
+              + R,
+                mesh.solver("Yi")
+            );
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+        else
+        {
+            inertIndex = i;
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+
+    fvScalarMatrix hsEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+        DpDt
+      + dQ
+      + radiation->Shs(thermo)
+      + parcels.Sh(hs)
+      + surfaceFilm.Sh()
+    );
+
+    hsEqn.relax();
+    hsEqn.solve();
+
+    thermo.correct();
+
+    radiation->correct();
+
+    Info<< "min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
+}

applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/createClouds.H

@@ -0,0 +1,9 @@
+Info<< "\nConstructing reacting cloud" << endl;
+basicReactingCloud parcels
+(
+    "reactingCloud1",
+    rho,
+    U,
+    g,
+    slgThermo
+);

applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/createFields.H

@@ -0,0 +1,154 @@
+    Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<hsCombustionThermo> pThermo
+    (
+        hsCombustionThermo::New(mesh)
+    );
+    hsCombustionThermo& thermo = pThermo();
+
+    SLGThermo slgThermo(mesh, thermo);
+
+    basicMultiComponentMixture& composition = thermo.composition();
+    PtrList<volScalarField>& Y = composition.Y();
+
+    const word inertSpecie(thermo.lookup("inertSpecie"));
+
+    Info<< "Creating field rho\n" << endl;
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+    volScalarField& p = thermo.p();
+    volScalarField& hs = thermo.hs();
+    const volScalarField& T = thermo.T();
+    const volScalarField& psi = thermo.psi();
+
+    Info<< "\nReading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "compressibleCreatePhi.H"
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    IOdictionary combustionProperties
+    (
+        IOobject
+        (
+            "combustionProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    );
+
+    Info<< "Creating combustion model\n" << endl;
+    autoPtr<combustionModel> combustion
+    (
+        combustionModel::combustionModel::New
+        (
+            combustionProperties,
+            thermo,
+            turbulence(),
+            phi,
+            rho
+        )
+    );
+
+    volScalarField dQ
+    (
+        IOobject
+        (
+            "dQ",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("dQ", dimMass/pow3(dimTime)/dimLength, 0.0)
+    );
+
+    Info<< "Creating field DpDt\n" << endl;
+    volScalarField DpDt
+    (
+        "DpDt",
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
+    );
+
+
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+
+    surfaceScalarField ghf("gh", g & mesh.Cf());
+
+    volScalarField p_rgh
+    (
+        IOobject
+        (
+            "p_rgh",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // Force p_rgh to be consistent with p
+    p_rgh = p - rho*gh;
+
+    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+    forAll(Y, i)
+    {
+        fields.add(Y[i]);
+    }
+    fields.add(hs);
+
+    IOdictionary additionalControlsDict
+    (
+        IOobject
+        (
+            "additionalControls",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    );
+
+    Switch solvePrimaryRegion
+    (
+        additionalControlsDict.lookup("solvePrimaryRegion")
+    );

applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/createPyrolysisModel.H

@@ -0,0 +1,6 @@
+Info<< "Creating pyrolysis model" << endl;
+
+autoPtr<regionModels::pyrolysisModels::pyrolysisModel> pyrolysis
+(
+    regionModels::pyrolysisModels::pyrolysisModel::New(mesh)
+);

applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/createSurfaceFilmModel.H

@@ -0,0 +1,7 @@
+Info<< "\nConstructing surface film model" << endl;
+
+typedef regionModels::surfaceFilmModels::surfaceFilmModel filmModelType;
+
+autoPtr<filmModelType> tsurfaceFilm(filmModelType::New(mesh, g));
+filmModelType& surfaceFilm = tsurfaceFilm();
+

applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/pEqn.H

@@ -0,0 +1,50 @@
+rho = thermo.rho();
+
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf(rAU.name(), fvc::interpolate(rho*rAU));
+U = rAU*UEqn.H();
+
+surfaceScalarField phiU
+(
+    fvc::interpolate(rho)
+   *(
+        (fvc::interpolate(U) & mesh.Sf())
+      + fvc::ddtPhiCorr(rAU, rho, U, phi)
+    )
+);
+
+phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
+
+for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
+{
+    fvScalarMatrix p_rghEqn
+    (
+        fvc::ddt(psi, rho)*gh
+      + fvc::div(phi)
+      + fvm::ddt(psi, p_rgh)
+      - fvm::laplacian(rhorAUf, p_rgh)
+     ==
+        parcels.Srho()
+      + surfaceFilm.Srho()
+    );
+
+    p_rghEqn.solve
+    (
+        mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
+    );
+
+    if (nonOrth == pimple.nNonOrthCorr())
+    {
+        phi += p_rghEqn.flux();
+    }
+}
+
+p = p_rgh + rho*gh;
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
+U.correctBoundaryConditions();
+
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/reactingParcelFilmPyrolysisFoam.C

@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    reactingParcelFilmPyrolysisFoam
+
+Description
+    Transient PISO solver for compressible, laminar or turbulent flow with
+    reacting Lagrangian parcels, surface film and pyrolysis modelling.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "turbulenceModel.H"
+#include "basicReactingCloud.H"
+#include "surfaceFilmModel.H"
+#include "pyrolysisModel.H"
+#include "radiationModel.H"
+#include "SLGThermo.H"
+#include "hsCombustionThermo.H"
+#include "solidChemistryModel.H"
+#include "combustionModel.H"
+#include "pimpleControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readChemistryProperties.H"
+    #include "readGravitationalAcceleration.H"
+    #include "createFields.H"
+    #include "createClouds.H"
+    #include "createRadiationModel.H"
+    #include "createSurfaceFilmModel.H"
+    #include "createPyrolysisModel.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+    #include "readPyrolysisTimeControls.H"
+
+    pimpleControl pimple(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "solidRegionDiffusionNo.H"
+        #include "setMultiRegionDeltaT.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        parcels.evolve();
+
+        surfaceFilm.evolve();
+
+        pyrolysis->evolve();
+
+        if (solvePrimaryRegion)
+        {
+            #include "rhoEqn.H"
+
+            // --- PIMPLE loop
+            for (pimple.start(); pimple.loop(); pimple++)
+            {
+                #include "UEqn.H"
+                #include "YhsEqn.H"
+
+                // --- PISO loop
+                for (int corr=1; corr<=pimple.nCorr(); corr++)
+                {
+                    #include "pEqn.H"
+                }
+
+                if (pimple.turbCorr())
+                {
+                    turbulence->correct();
+                }
+            }
+
+            rho = thermo.rho();
+        }
+        else
+        {
+            runTime.write();
+        }
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/readChemistryProperties.H

@@ -0,0 +1,23 @@
+Info<< "Reading chemistry properties\n" << endl;
+
+IOdictionary chemistryProperties
+(
+    IOobject
+    (
+        "chemistryProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE,
+        false
+    )
+);
+
+Switch turbulentReaction(chemistryProperties.lookup("turbulentReaction"));
+
+dimensionedScalar Cmix("Cmix", dimless, 1.0);
+
+if (turbulentReaction)
+{
+    chemistryProperties.lookup("Cmix") >> Cmix;
+}

applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/readPyrolysisTimeControls.H

@@ -0,0 +1,34 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    readPyrolysisTimeControls
+
+Description
+
+
+\*---------------------------------------------------------------------------*/
+
+scalar maxDi = pyrolysis->maxDiff();
+
+// ************************************************************************* //

applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/rhoEqn.H

@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    rhoEqn
+
+Description
+    Solve the continuity for density.
+
+\*---------------------------------------------------------------------------*/
+
+{
+    solve
+    (
+        fvm::ddt(rho)
+      + fvc::div(phi)
+      ==
+        parcels.Srho(rho)
+      + surfaceFilm.Srho()
+    );
+}
+
+// ************************************************************************* //

applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/setMultiRegionDeltaT.H

@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    setMultiRegionDeltaT
+
+Description
+    Reset the timestep to maintain a constant maximum Courant numbers.
+    Reduction of time-step is immediate, but increase is damped to avoid
+    unstable oscillations.
+
+\*---------------------------------------------------------------------------*/
+
+if (adjustTimeStep)
+{
+    if (CoNum == -GREAT)
+    {
+        CoNum = SMALL;
+    }
+
+    if (DiNum == -GREAT)
+    {
+        DiNum = SMALL;
+    }
+
+
+    const scalar TFactorFluid = maxCo/(CoNum + SMALL);
+    const scalar TFactorSolid = maxDi/(DiNum + SMALL);
+    const scalar TFactorFilm = maxCo/(surfaceFilm.CourantNumber() + SMALL);
+
+    const scalar dt0 = runTime.deltaTValue();
+
+    runTime.setDeltaT
+    (
+        min
+        (
+            dt0*min(min(TFactorFluid, min(TFactorFilm, TFactorSolid)), 1.2),
+            maxDeltaT
+        )
+    );
+}
+
+
+// ************************************************************************* //

applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/solidRegionDiffusionNo.H

@@ -0,0 +1,2 @@
+scalar DiNum =  pyrolysis->solidRegionDiffNo();
+