diff --git a/applications/solvers/Lagrangian/trackedReactingParcelFoam/Make/files b/applications/solvers/Lagrangian/trackedReactingParcelFoam/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..4bba4a1cc67a87dc9faa94e0b8734795709b1094
--- /dev/null
+++ b/applications/solvers/Lagrangian/trackedReactingParcelFoam/Make/files
@@ -0,0 +1,3 @@
+trackedReactingParcelFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/trackedReactingParcelFoam
diff --git a/applications/solvers/Lagrangian/trackedReactingParcelFoam/Make/options b/applications/solvers/Lagrangian/trackedReactingParcelFoam/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..b0e474c308f5b7208b8f099789ca3247bae306fa
--- /dev/null
+++ b/applications/solvers/Lagrangian/trackedReactingParcelFoam/Make/options
@@ -0,0 +1,39 @@
+EXE_INC = \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I${LIB_SRC}/meshTools/lnInclude \
+ -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+ -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+ -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+ -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/pdfs/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/liquids/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/liquidMixture/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/solids/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/solidMixture/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+ -I$(LIB_SRC)/ODE/lnInclude
+
+EXE_LIBS = \
+ -L$(FOAM_USER_LIBBIN) \
+ -lfiniteVolume \
+ -lmeshTools \
+ -lcompressibleRASModels \
+ -lcompressibleLESModels \
+ -llagrangian \
+ -llagrangianIntermediate \
+ -lspecie \
+ -lbasicThermophysicalModels \
+ -lliquids \
+ -lliquidMixture \
+ -lsolids \
+ -lsolidMixture \
+ -lthermophysicalFunctions \
+ -lcombustionThermophysicalModels \
+ -lchemistryModel \
+ -lradiation \
+ -lODE
diff --git a/applications/solvers/Lagrangian/trackedReactingParcelFoam/UEqn.H b/applications/solvers/Lagrangian/trackedReactingParcelFoam/UEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..b34cd840019803bf5a04c431e4472b089b511ee0
--- /dev/null
+++ b/applications/solvers/Lagrangian/trackedReactingParcelFoam/UEqn.H
@@ -0,0 +1,16 @@
+ fvVectorMatrix UEqn
+ (
+ fvm::ddt(rho, U)
+ + fvm::div(phi, U)
+ + turbulence->divDevRhoReff(U)
+ ==
+ rho.dimensionedInternalField()*g
+ + reactingParcels.SU()
+ );
+
+ UEqn.relax();
+
+ if (momentumPredictor)
+ {
+ solve(UEqn == -fvc::grad(p));
+ }
diff --git a/applications/solvers/Lagrangian/trackedReactingParcelFoam/YEqn.H b/applications/solvers/Lagrangian/trackedReactingParcelFoam/YEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..687acd4e70a4ec429839065c244d6c636a20e9df
--- /dev/null
+++ b/applications/solvers/Lagrangian/trackedReactingParcelFoam/YEqn.H
@@ -0,0 +1,45 @@
+
+tmp<fv::convectionScheme<scalar> > mvConvection
+(
+ fv::convectionScheme<scalar>::New
+ (
+ mesh,
+ fields,
+ phi,
+ mesh.divScheme("div(phi,Yi_h)")
+ )
+);
+
+
+{
+ label inertIndex = -1;
+ volScalarField Yt = 0.0*Y[0];
+
+ for (label i=0; i<Y.size(); i++)
+ {
+ if (Y[i].name() != inertSpecie)
+ {
+ volScalarField& Yi = Y[i];
+ solve
+ (
+ fvm::ddt(rho, Yi)
+ + mvConvection->fvmDiv(phi, Yi)
+ - fvm::laplacian(turbulence->muEff(), Yi)
+ ==
+ reactingParcels.Srho(i)
+ + kappa*chemistry.RR(i)().dimensionedInternalField()
+ );
+
+ Yi.max(0.0);
+ Yt += Yi;
+ }
+ else
+ {
+ inertIndex = i;
+ }
+ }
+
+ Y[inertIndex] = scalar(1) - Yt;
+ Y[inertIndex].max(0.0);
+
+}
diff --git a/applications/solvers/Lagrangian/trackedReactingParcelFoam/additionalOutput.H b/applications/solvers/Lagrangian/trackedReactingParcelFoam/additionalOutput.H
new file mode 100644
index 0000000000000000000000000000000000000000..9edd35eb7cab9acd294ba6a9e73b6ee8087b1c0d
--- /dev/null
+++ b/applications/solvers/Lagrangian/trackedReactingParcelFoam/additionalOutput.H
@@ -0,0 +1,48 @@
+{
+ tmp<volScalarField> tdQ
+ (
+ new volScalarField
+ (
+ IOobject
+ (
+ "dQ",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar
+ (
+ "zero",
+ dimensionSet(1, -3, -1, 0, 0, 0, 0),
+ 0.0
+ )
+ )
+ );
+
+ scalarField& dQ = tdQ();
+
+ scalarField cp(dQ.size(), 0.0);
+
+ forAll(Y, i)
+ {
+ volScalarField RRi = chemistry.RR(i);
+
+ forAll(h, celli)
+ {
+ scalar Ti = T[celli];
+ cp[celli] += Y[i][celli]*chemistry.specieThermo()[i].Cp(Ti);
+ scalar hi = chemistry.specieThermo()[i].h(Ti);
+ scalar RR = RRi[celli];
+ dQ[celli] -= hi*RR;
+ }
+ }
+
+ forAll(dQ, celli)
+ {
+ dQ[celli] /= cp[celli];
+ }
+
+ tdQ().write();
+}
diff --git a/applications/solvers/Lagrangian/trackedReactingParcelFoam/chemistry.H b/applications/solvers/Lagrangian/trackedReactingParcelFoam/chemistry.H
new file mode 100644
index 0000000000000000000000000000000000000000..07b1e9953b0db867186f6c668d27a9415a26c265
--- /dev/null
+++ b/applications/solvers/Lagrangian/trackedReactingParcelFoam/chemistry.H
@@ -0,0 +1,25 @@
+{
+ Info << "Solving chemistry" << endl;
+
+ chemistry.solve
+ (
+ runTime.value() - runTime.deltaT().value(),
+ runTime.deltaT().value()
+ );
+
+ // turbulent time scale
+ if (turbulentReaction)
+ {
+ DimensionedField<scalar, volMesh> tk =
+ Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
+ DimensionedField<scalar, volMesh> tc =
+ chemistry.tc()().dimensionedInternalField();
+
+ // Chalmers PaSR model
+ kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
+ }
+ else
+ {
+ kappa = 1.0;
+ }
+}
diff --git a/applications/solvers/Lagrangian/trackedReactingParcelFoam/createClouds.H b/applications/solvers/Lagrangian/trackedReactingParcelFoam/createClouds.H
new file mode 100644
index 0000000000000000000000000000000000000000..e8a02c79ab8a9bfd2f7e377dd3aec8eb324d9b18
--- /dev/null
+++ b/applications/solvers/Lagrangian/trackedReactingParcelFoam/createClouds.H
@@ -0,0 +1,42 @@
+Info<< "\nConstructing interpolation" << endl;
+
+Info << "\nConstructing gas properties" << endl;
+/*
+PtrList<specieConstProperties> gasProperties(Y.size());
+forAll(gasProperties, i)
+{
+ gasProperties.set
+ (
+ i,
+ new specieConstProperties
+ (
+ dynamic_cast<const multiComponentMixture<constTransport<
+ specieThermo<hConstThermo<perfectGas> > > >&>
+ (thermo()).speciesData()[i]
+ )
+ );
+}
+*/
+PtrList<specieReactingProperties> gasProperties(Y.size());
+forAll(gasProperties, i)
+{
+ gasProperties.set
+ (
+ i,
+ new specieReactingProperties
+ (
+ dynamic_cast<const reactingMixture&>(thermo()).speciesData()[i]
+ )
+ );
+}
+
+Info<< "\nConstructing reacting cloud" << endl;
+trackedReactingCloud reactingParcels
+(
+ "reactingCloud1",
+ rho,
+ U,
+ g,
+ thermo(),
+ gasProperties
+);
diff --git a/applications/solvers/Lagrangian/trackedReactingParcelFoam/createFields.H b/applications/solvers/Lagrangian/trackedReactingParcelFoam/createFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..079a4f2e1d1ba7b48486909490f1abcd7a7f0352
--- /dev/null
+++ b/applications/solvers/Lagrangian/trackedReactingParcelFoam/createFields.H
@@ -0,0 +1,141 @@
+ Info<< "Reading thermophysical properties\n" << endl;
+
+ autoPtr<hCombustionThermo> thermo
+ (
+ hCombustionThermo::New(mesh)
+ );
+
+ combustionMixture& composition = thermo->composition();
+ PtrList<volScalarField>& Y = composition.Y();
+
+ word inertSpecie(thermo->lookup("inertSpecie"));
+
+ volScalarField& p = thermo->p();
+ volScalarField& h = thermo->h();
+ const volScalarField& T = thermo->T();
+ const volScalarField& psi = thermo->psi();
+
+ volScalarField rho
+ (
+ IOobject
+ (
+ "rho",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ 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;
+ volVectorField U
+ (
+ IOobject
+ (
+ "U",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+ );
+
+# include "compressibleCreatePhi.H"
+
+ DimensionedField<scalar, volMesh> kappa
+ (
+ IOobject
+ (
+ "kappa",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh,
+ dimensionedScalar("zero", dimless, 0.0)
+ );
+
+ Info<< "Creating turbulence model\n" << endl;
+ autoPtr<compressible::turbulenceModel> turbulence
+ (
+ compressible::turbulenceModel::New
+ (
+ rho,
+ U,
+ phi,
+ thermo()
+ )
+ );
+
+ Info<< "Creating field DpDt\n" << endl;
+ volScalarField DpDt =
+ fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+ Info << "Constructing chemical mechanism" << endl;
+ chemistryModel chemistry
+ (
+ thermo(),
+ rho
+ );
+
+ multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+ forAll (Y, i)
+ {
+ fields.add(Y[i]);
+ }
+ fields.add(h);
+
+ Info<< "Creating radiation model\n" << endl;
+ autoPtr<radiation::radiationModel> radiation
+ (
+ radiation::radiationModel::New(T)
+ );
diff --git a/applications/solvers/Lagrangian/trackedReactingParcelFoam/hEqn.H b/applications/solvers/Lagrangian/trackedReactingParcelFoam/hEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..32e71555369290be7d019af0f4d6ad15473ee064
--- /dev/null
+++ b/applications/solvers/Lagrangian/trackedReactingParcelFoam/hEqn.H
@@ -0,0 +1,20 @@
+{
+ fvScalarMatrix hEqn
+ (
+ fvm::ddt(rho, h)
+ + fvm::div(phi, h)
+ - fvm::laplacian(turbulence->alphaEff(), h)
+ ==
+ DpDt
+ + reactingParcels.Sh()
+ + radiation->Sh(thermo())
+ );
+
+ hEqn.relax();
+
+ hEqn.solve();
+
+ thermo->correct();
+
+ radiation->correct();
+}
diff --git a/applications/solvers/Lagrangian/trackedReactingParcelFoam/pEqn.H b/applications/solvers/Lagrangian/trackedReactingParcelFoam/pEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..3dee469329ccf12db6b1cfc74a55c4b919af1ce4
--- /dev/null
+++ b/applications/solvers/Lagrangian/trackedReactingParcelFoam/pEqn.H
@@ -0,0 +1,72 @@
+rho = thermo->rho();
+
+volScalarField rUA = 1.0/UEqn.A();
+U = rUA*UEqn.H();
+
+if (transonic)
+{
+ surfaceScalarField phid
+ (
+ "phid",
+ fvc::interpolate(thermo->psi())
+ *(
+ (fvc::interpolate(U) & mesh.Sf())
+ + fvc::ddtPhiCorr(rUA, rho, U, phi)
+ )
+ );
+
+ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+ {
+ fvScalarMatrix pEqn
+ (
+ fvm::ddt(psi, p)
+ + fvm::div(phid, p)
+ - fvm::laplacian(rho*rUA, p)
+ ==
+ reactingParcels.Srho()
+ );
+
+ pEqn.solve();
+
+ if (nonOrth == nNonOrthCorr)
+ {
+ phi == pEqn.flux();
+ }
+ }
+}
+else
+{
+ phi =
+ fvc::interpolate(rho)*
+ (
+ (fvc::interpolate(U) & mesh.Sf())
+ + fvc::ddtPhiCorr(rUA, rho, U, phi)
+ );
+
+ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+ {
+ fvScalarMatrix pEqn
+ (
+ fvm::ddt(psi, p)
+ + fvc::div(phi)
+ - fvm::laplacian(rho*rUA, p)
+ ==
+ reactingParcels.Srho()
+ );
+
+ pEqn.solve();
+
+ if (nonOrth == nNonOrthCorr)
+ {
+ phi += pEqn.flux();
+ }
+ }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+U -= rUA*fvc::grad(p);
+U.correctBoundaryConditions();
+
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/applications/solvers/Lagrangian/trackedReactingParcelFoam/readChemistryProperties.H b/applications/solvers/Lagrangian/trackedReactingParcelFoam/readChemistryProperties.H
new file mode 100644
index 0000000000000000000000000000000000000000..1a60e6fb34645a004fd39321f7a54d3bd5b45381
--- /dev/null
+++ b/applications/solvers/Lagrangian/trackedReactingParcelFoam/readChemistryProperties.H
@@ -0,0 +1,22 @@
+Info<< "Reading chemistry properties\n" << endl;
+
+IOdictionary chemistryProperties
+(
+ IOobject
+ (
+ "chemistryProperties",
+ runTime.constant(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ )
+);
+
+Switch turbulentReaction(chemistryProperties.lookup("turbulentReaction"));
+
+dimensionedScalar Cmix("Cmix", dimless, 1.0);
+
+if (turbulentReaction)
+{
+ chemistryProperties.lookup("Cmix") >> Cmix;
+}
diff --git a/applications/solvers/Lagrangian/trackedReactingParcelFoam/rhoEqn.H b/applications/solvers/Lagrangian/trackedReactingParcelFoam/rhoEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..03eed2dc8219a3e6e836c5b48089ba6e97a90c3f
--- /dev/null
+++ b/applications/solvers/Lagrangian/trackedReactingParcelFoam/rhoEqn.H
@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+ Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Global
+ rhoEqn
+
+Description
+ Solve the continuity for density.
+
+\*---------------------------------------------------------------------------*/
+
+{
+ solve
+ (
+ fvm::ddt(rho)
+ + fvc::div(phi)
+ ==
+ reactingParcels.Srho()
+ );
+}
+
+// ************************************************************************* //
diff --git a/applications/solvers/Lagrangian/trackedReactingParcelFoam/trackedReactingParcelFoam.C b/applications/solvers/Lagrangian/trackedReactingParcelFoam/trackedReactingParcelFoam.C
new file mode 100644
index 0000000000000000000000000000000000000000..07838ea813e41a6888893cf9d25bc0fca169bfa6
--- /dev/null
+++ b/applications/solvers/Lagrangian/trackedReactingParcelFoam/trackedReactingParcelFoam.C
@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+ Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "hCombustionThermo.H"
+#include "turbulenceModel.H"
+#include "trackedReactingCloud.H"
+#include "chemistryModel.H"
+#include "chemistrySolver.H"
+#include "ReactingCloudThermoTypes.H"
+#include "radiationModel.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"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ Info<< "\nStarting time loop\n" << endl;
+
+ while (runTime.run())
+ {
+# 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"
+
+ // --- PIMPLE loop
+ for (int ocorr=1; ocorr<=nOuterCorr; ocorr++)
+ {
+# include "UEqn.H"
+# include "YEqn.H"
+
+ // --- PISO loop
+ for (int corr=1; corr<=nCorr; corr++)
+ {
+# include "hEqn.H"
+# include "pEqn.H"
+ }
+
+ Info<< "T gas min/max = " << min(T).value() << ", "
+ << max(T).value() << endl;
+ }
+
+ turbulence->correct();
+
+ rho = thermo->rho();
+
+ if (runTime.write())
+ {
+# include "additionalOutput.H"
+ }
+
+ Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+ << " ClockTime = " << runTime.elapsedClockTime() << " s"
+ << nl << endl;
+ }
+
+ Info<< "End\n" << endl;
+
+ return(0);
+}
+
+
+// ************************************************************************* //