diff --git a/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/Make/files b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..7910a0296578f689ba0cb09fe8212f09a2eae6ae
--- /dev/null
+++ b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/Make/files
@@ -0,0 +1,3 @@
+reactingParcelFilmPyrolysisFoam.C
+
+EXE = $(FOAM_APPBIN)/reactingParcelFilmPyrolysisFoam
diff --git a/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/Make/options b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..7782cc45ccb80be0afa7df6874b97d62551a200a
--- /dev/null
+++ b/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
diff --git a/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/UEqn.H b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/UEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..a64e50a2d24d19b3f1dd25e4bc73acbcfc27a356
--- /dev/null
+++ b/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()
+ )
+ );
+ }
diff --git a/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/YhsEqn.H b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/YhsEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..ae70f37d8cd3293765fd439594e0c3bee7697f09
--- /dev/null
+++ b/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;
+}
diff --git a/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/createClouds.H b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/createClouds.H
new file mode 100644
index 0000000000000000000000000000000000000000..c568be12a16faa253b69d064499038de506a36b5
--- /dev/null
+++ b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/createClouds.H
@@ -0,0 +1,9 @@
+Info<< "\nConstructing reacting cloud" << endl;
+basicReactingCloud parcels
+(
+ "reactingCloud1",
+ rho,
+ U,
+ g,
+ slgThermo
+);
diff --git a/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/createFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..96a3457b8378d863a4d44024cd9fc4d5e3730d09
--- /dev/null
+++ b/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")
+ );
diff --git a/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/createPyrolysisModel.H b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/createPyrolysisModel.H
new file mode 100644
index 0000000000000000000000000000000000000000..e5b8d4e45793215a334c71b2c42d5939272e7f64
--- /dev/null
+++ b/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)
+);
diff --git a/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/createSurfaceFilmModel.H b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/createSurfaceFilmModel.H
new file mode 100644
index 0000000000000000000000000000000000000000..1db870f0aa76f02a89cc98748782f97fd0f1242b
--- /dev/null
+++ b/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();
+
diff --git a/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/pEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..bbfb5be051c06d5378b0e15d7f33cff76a19f1a3
--- /dev/null
+++ b/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);
diff --git a/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/reactingParcelFilmPyrolysisFoam.C b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/reactingParcelFilmPyrolysisFoam.C
new file mode 100644
index 0000000000000000000000000000000000000000..6a46b2fbfb157574c01fabda9f7a8c406c0c8304
--- /dev/null
+++ b/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);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/readChemistryProperties.H b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/readChemistryProperties.H
new file mode 100644
index 0000000000000000000000000000000000000000..f0bcf7597fcf71f1e9b8ee2dbc879200a85fa2cc
--- /dev/null
+++ b/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;
+}
diff --git a/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/readPyrolysisTimeControls.H b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/readPyrolysisTimeControls.H
new file mode 100644
index 0000000000000000000000000000000000000000..43a4e9e36335b5aeb3821f0641922510e0de1e48
--- /dev/null
+++ b/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();
+
+// ************************************************************************* //
diff --git a/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/rhoEqn.H b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/rhoEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..617f1df8a64a0a2dd875e401948334968037fd3b
--- /dev/null
+++ b/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()
+ );
+}
+
+// ************************************************************************* //
diff --git a/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/setMultiRegionDeltaT.H b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/setMultiRegionDeltaT.H
new file mode 100644
index 0000000000000000000000000000000000000000..6c5a6a7215c09a4b77ea90abec664099f7c39e35
--- /dev/null
+++ b/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
+ )
+ );
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/solidRegionDiffusionNo.H b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/solidRegionDiffusionNo.H
new file mode 100644
index 0000000000000000000000000000000000000000..a6ab0eb14a6741700fa0379190cbc7b52b21426d
--- /dev/null
+++ b/applications/solvers/lagrangian/reactingParcelFilmPyrolysisFoam/solidRegionDiffusionNo.H
@@ -0,0 +1,2 @@
+scalar DiNum = pyrolysis->solidRegionDiffNo();
+