Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev

applications/solvers/lagrangian/sprayFoam/Make/files

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

applications/solvers/lagrangian/sprayFoam/Make/options

@@ -0,0 +1,50 @@
+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/spray/lnInclude \
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(FOAM_SOLVERS)/lagrangian/reactingParcelFoam
+
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lcompressibleTurbulenceModel \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -llagrangianSpray \
+    -lspecie \
+    -lbasicThermophysicalModels \
+    -lliquidProperties \
+    -lliquidMixtureProperties \
+    -lsolidProperties \
+    -lsolidMixtureProperties \
+    -lthermophysicalFunctions \
+    -lreactionThermophysicalModels \
+    -lSLGThermo \
+    -lchemistryModel \
+    -lradiationModels \
+    -lODE \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -lsampling

applications/solvers/lagrangian/sprayFoam/chemistry.H

@@ -0,0 +1,44 @@
+if (chemistry.chemistry())
+{
+    Info<< "Solving chemistry" << endl;
+
+    chemistry.solve
+    (
+        runTime.value() - runTime.deltaTValue(),
+        runTime.deltaTValue()
+    );
+
+    // turbulent time scale
+    if (turbulentReaction)
+    {
+        tmp<volScalarField> tepsilon(turbulence->epsilon());
+        const volScalarField& epsilon = tepsilon();
+        tmp<volScalarField> tmuEff(turbulence->muEff());
+        const volScalarField& muEff = tmuEff();
+        tmp<volScalarField> ttc(chemistry.tc());
+        const volScalarField& tc = ttc();
+
+        forAll(epsilon, i)
+        {
+            if (epsilon[i] > 0)
+            {
+                // Chalmers PaSR model
+                scalar tk = Cmix.value()*Foam::sqrt(muEff[i]/rho[i]/epsilon[i]);
+                kappa[i] =
+                    (runTime.deltaTValue() + tc[i])
+                   /(runTime.deltaTValue() + tc[i] + tk);
+            }
+            else
+            {
+                // Return to laminar combustion
+                kappa[i] = 1.0;
+            }
+        }
+    }
+    else
+    {
+        kappa = 1.0;
+    }
+
+    chemistrySh = kappa*chemistry.Sh()();
+}

applications/solvers/lagrangian/sprayFoam/createClouds.H

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

applications/solvers/lagrangian/sprayFoam/sprayFoam.C

@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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
+    sprayFoam
+
+Description
+    Transient PIMPLE solver for compressible, laminar or turbulent flow with
+    spray parcels.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "hCombustionThermo.H"
+#include "turbulenceModel.H"
+#include "basicSprayCloud.H"
+#include "psiChemistryModel.H"
+#include "chemistrySolver.H"
+#include "radiationModel.H"
+#include "SLGThermo.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 "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    pimpleControl pimple(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        parcels.evolve();
+
+        #include "chemistry.H"
+        #include "rhoEqn.H"
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (pimple.start(); pimple.loop(); pimple++)
+        {
+            #include "UEqn.H"
+            #include "YEqn.H"
+            #include "hsEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<pimple.nCorr(); corr++)
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        rho = thermo.rho();
+
+        if (runTime.write())
+        {
+            chemistry.dQ()().write();
+        }
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //

src/lagrangian/Allwmake

@@ -7,6 +7,7 @@ wmake $makeType distributionModels
 wmake $makeType basic
 wmake $makeType solidParticle
 wmake $makeType intermediate
+wmake $makeType spray
 wmake $makeType dieselSpray
 wmake $makeType dsmc
 wmake $makeType coalCombustion

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C

@@ -126,17 +126,6 @@ void Foam::KinematicCloud<CloudType>::solve(TrackData& td)
 }
 
 
-template<class CloudType>
-void Foam::KinematicCloud<CloudType>::preEvolve()
-{
-    Info<< "\nSolving cloud " << this->name() << endl;
-
-    this->dispersion().cacheFields(true);
-    forces_.cacheFields(true);
-    updateCellOccupancy();
-}
-
-
 template<class CloudType>
 void Foam::KinematicCloud<CloudType>::buildCellOccupancy()
 {
@@ -310,6 +299,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
     U_(U),
     mu_(mu),
     g_(g),
+    pAmbient_(0.0),
     forces_
     (
         *this,
@@ -403,6 +393,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
     U_(c.U_),
     mu_(c.mu_),
     g_(c.g_),
+    pAmbient_(c.pAmbient_),
     forces_(c.forces_),
     functions_(c.functions_),
     dispersionModel_(c.dispersionModel_->clone()),
@@ -478,6 +469,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
     U_(c.U_),
     mu_(c.mu_),
     g_(c.g_),
+    pAmbient_(c.pAmbient_),
     forces_(*this, mesh),
     functions_(*this),
     dispersionModel_(NULL),
@@ -507,7 +499,7 @@ bool Foam::KinematicCloud<CloudType>::hasWallImpactDistance() const
 
 
 template<class CloudType>
-void Foam::KinematicCloud<CloudType>::checkParcelProperties
+void Foam::KinematicCloud<CloudType>::setParcelThermoProperties
 (
     parcelType& parcel,
     const scalar lagrangianDt,
@@ -518,7 +510,17 @@ void Foam::KinematicCloud<CloudType>::checkParcelProperties
     {
         parcel.rho() = constProps_.rho0();
     }
+}
 
+
+template<class CloudType>
+void Foam::KinematicCloud<CloudType>::checkParcelProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt,
+    const bool fullyDescribed
+)
+{
     const scalar carrierDt = mesh_.time().deltaTValue();
     parcel.stepFraction() = (carrierDt - lagrangianDt)/carrierDt;
     parcel.typeId() = constProps_.parcelTypeId();
@@ -600,6 +602,20 @@ void Foam::KinematicCloud<CloudType>::scaleSources()
 }
 
 
+template<class CloudType>
+void Foam::KinematicCloud<CloudType>::preEvolve()
+{
+    Info<< "\nSolving cloud " << this->name() << endl;
+
+    this->dispersion().cacheFields(true);
+    forces_.cacheFields(true);
+    updateCellOccupancy();
+
+    pAmbient_ = constProps_.dict().template
+        lookupOrDefault<scalar>("pAmbient", pAmbient_);
+}
+
+
 template<class CloudType>
 void Foam::KinematicCloud<CloudType>::evolve()
 {

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H

@@ -176,6 +176,9 @@ protected:
             //- Gravity
             const dimensionedVector& g_;
 
+            //- Averaged ambient domain pressure
+            scalar pAmbient_;
+
 
         //- Optional particle forces
         forceType forces_;
@@ -230,9 +233,6 @@ protected:
             template<class TrackData>
             void solve(TrackData& td);
 
-            //- Pre-evolve
-            void preEvolve();
-
             //- Build the cellOccupancy
             void buildCellOccupancy();
 
@@ -367,6 +367,12 @@ public:
                 //- Gravity
                 inline const dimensionedVector& g() const;
 
+                //- Return const-access to the ambient pressure
+                inline scalar pAmbient() const;
+
+                //- Return reference to the ambient pressure
+                inline scalar& pAmbient();
+
 
             //- Optional particle forces
 //            inline const typename parcelType::forceType& forces() const;
@@ -456,6 +462,12 @@ public:
             //- Total rotational kinetic energy in the system
             inline scalar rotationalKineticEnergyOfSystem() const;
 
+            //- Penetration for percentage of the current total mass
+            inline scalar penetration(const scalar& prc) const;
+
+            //- Mean diameter Dij
+            inline scalar Dij(const label i, const label j) const;
+
 
             // Fields
 
@@ -474,6 +486,14 @@ public:
 
         // Cloud evolution functions
 
+            //- Set parcel thermo properties
+            void setParcelThermoProperties
+            (
+                parcelType& parcel,
+                const scalar lagrangianDt,
+                const bool fullyDescribed
+            );
+
             //- Check parcel properties
             void checkParcelProperties
             (
@@ -514,6 +534,9 @@ public:
             //- Apply scaling to (transient) cloud sources
             void scaleSources();
 
+            //- Pre-evolve
+            void preEvolve();
+
             //- Evolve the cloud
             void evolve();
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H

@@ -109,6 +109,20 @@ inline const Foam::dimensionedVector& Foam::KinematicCloud<CloudType>::g() const
 }
 
 
+template<class CloudType>
+inline Foam::scalar Foam::KinematicCloud<CloudType>::pAmbient() const
+{
+    return pAmbient_;
+}
+
+
+template<class CloudType>
+inline Foam::scalar& Foam::KinematicCloud<CloudType>::pAmbient()
+{
+    return pAmbient_;
+}
+
+
 template<class CloudType>
 //inline const typename CloudType::parcelType::forceType&
 inline const typename Foam::KinematicCloud<CloudType>::forceType&
@@ -264,6 +278,131 @@ Foam::KinematicCloud<CloudType>::rotationalKineticEnergyOfSystem() const
 }
 
 
+template<class CloudType>
+inline Foam::scalar Foam::KinematicCloud<CloudType>::Dij
+(
+    const label i,
+    const label j
+) const
+{
+    scalar si = 0.0;
+    scalar sj = 0.0;
+    forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
+    {
+        const parcelType& p = iter();
+        si += p.nParticle()*pow(p.d(), i);
+        sj += p.nParticle()*pow(p.d(), j);
+    }
+
+    reduce(si, sumOp<scalar>());
+    reduce(sj, sumOp<scalar>());
+    sj = max(sj, VSMALL);
+
+    return si/sj;
+}
+
+
+template<class CloudType>
+inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration
+(
+    const scalar& prc
+) const
+{
+    scalar distance = 0.0;
+    scalar mTot = 0.0;
+
+    label np = this->size();
+
+    // arrays containing the parcels mass and
+    // distance from injector in ascending order
+    scalarField mass(np);
+    scalarField dist(np);
+
+    if (np > 0)
+    {
+        label n = 0;
+
+        // first arrange the parcels in ascending order
+        // the first parcel is closest to its injection position
+        // and the last one is most far away.
+        forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
+        {
+            const parcelType& p = iter();
+            scalar mi = p.nParticle()*p.mass();
+            scalar di = mag(p.position() - p.position0());
+            mTot += mi;
+
+            // insert at the last place
+            mass[n] = mi;
+            dist[n] = di;
+
+            label i = 0;
+            bool found = false;
+
+            // insert the parcel in the correct place
+            // and move the others
+            while ((i < n) && (!found))
+            {
+                if (di < dist[i])
+                {
+                    found = true;
+                    for (label j=n; j>i; j--)
+                    {
+                        mass[j] = mass[j-1];
+                        dist[j] = dist[j-1];
+                    }
+                    mass[i] = mi;
+                    dist[i] = di;
+                }
+                i++;
+            }
+            n++;
+        }
+    }
+
+    reduce(mTot, sumOp<scalar>());
+
+    if (np > 0)
+    {
+        scalar mLimit = prc*mTot;
+        scalar mOff = (1.0 - prc)*mTot;
+
+        if (np > 1)
+        {
+            // 'prc' is large enough that the parcel most far
+            // away will be used, no need to loop...
+            if (mLimit > mTot - mass[np-1])
+            {
+                distance = dist[np-1];
+            }
+            else
+            {
+                scalar mOffSum = 0.0;
+                label i = np;
+
+                while ((mOffSum < mOff) && (i>0))
+                {
+                    i--;
+                    mOffSum += mass[i];
+                }
+                distance =
+                    dist[i+1]
+                  + (dist[i] - dist[i+1])*(mOffSum - mOff)
+                   /mass[i+1] ;
+            }
+        }
+        else
+        {
+            distance = dist[0];
+        }
+    }
+
+    reduce(distance, maxOp<scalar>());
+
+    return distance;
+}
+
+
 template<class CloudType>
 inline Foam::cachedRandom& Foam::KinematicCloud<CloudType>::rndGen()
 {

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.H

@@ -174,6 +174,9 @@ public:
             //- Return const access to the coupled flag
             inline const Switch coupled() const;
 
+            //- Return non-const access to the coupled flag
+            inline Switch& coupled();
+
             //- Return const access to the cell value correction flag
             inline const Switch cellValueSourceCorrection() const;
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolutionI.H

@@ -101,6 +101,12 @@ inline Foam::scalar Foam::cloudSolution::trackTime() const
 }
 
 
+inline Foam::Switch& Foam::cloudSolution::coupled()
+{
+    return coupled_;
+}
+
+
 inline const Foam::Switch Foam::cloudSolution::coupled() const
 {
     return coupled_;

src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C

@@ -224,22 +224,18 @@ Foam::ReactingCloud<CloudType>::~ReactingCloud()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
-void Foam::ReactingCloud<CloudType>::checkParcelProperties
+void Foam::ReactingCloud<CloudType>::setParcelThermoProperties
 (
     parcelType& parcel,
     const scalar lagrangianDt,
     const bool fullyDescribed
 )
 {
-    CloudType::checkParcelProperties
+    CloudType::setParcelThermoProperties(parcel, lagrangianDt, fullyDescribed);
-    (
-        parcel,
-        lagrangianDt,
-        fullyDescribed
-    );
 
     if (!fullyDescribed)
     {
+        parcel.pc() = this->thermo().thermo().p()[parcel.cell()];
         parcel.Y() = composition().YMixture0();
     }
     else
@@ -257,6 +253,18 @@ void Foam::ReactingCloud<CloudType>::checkParcelProperties
 }
 
 
+template<class CloudType>
+void Foam::ReactingCloud<CloudType>::checkParcelProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt,
+    const bool fullyDescribed
+)
+{
+    CloudType::checkParcelProperties(parcel, lagrangianDt, fullyDescribed);
+}
+
+
 template<class CloudType>
 void Foam::ReactingCloud<CloudType>::storeState()
 {

src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H

@@ -267,6 +267,14 @@ public:
 
         // Cloud evolution functions
 
+            //- Set parcel thermo properties
+            void setParcelThermoProperties
+            (
+                parcelType& parcel,
+                const scalar lagrangianDt,
+                const bool fullyDescribed
+            );
+
             //- Check parcel properties
             void checkParcelProperties
             (

src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C

@@ -155,19 +155,14 @@ Foam::ReactingMultiphaseCloud<CloudType>::~ReactingMultiphaseCloud()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
-void Foam::ReactingMultiphaseCloud<CloudType>::checkParcelProperties
+void Foam::ReactingMultiphaseCloud<CloudType>::setParcelThermoProperties
 (
     parcelType& parcel,
     const scalar lagrangianDt,
     const bool fullyDescribed
 )
 {
-    CloudType::checkParcelProperties
+    CloudType::setParcelThermoProperties(parcel, lagrangianDt, fullyDescribed);
-    (
-        parcel,
-        lagrangianDt,
-        fullyDescribed
-    );
 
     label idGas = this->composition().idGas();
     label idLiquid = this->composition().idLiquid();
@@ -203,6 +198,18 @@ void Foam::ReactingMultiphaseCloud<CloudType>::checkParcelProperties
 }
 
 
+template<class CloudType>
+void Foam::ReactingMultiphaseCloud<CloudType>::checkParcelProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt,
+    const bool fullyDescribed
+)
+{
+    CloudType::checkParcelProperties(parcel, lagrangianDt, fullyDescribed);
+}
+
+
 template<class CloudType>
 void Foam::ReactingMultiphaseCloud<CloudType>::storeState()
 {

src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.H

@@ -41,7 +41,6 @@ SourceFiles
 #ifndef ReactingMultiphaseCloud_H
 #define ReactingMultiphaseCloud_H
 
-#include "ReactingCloud.H"
 #include "reactingMultiphaseCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -241,6 +240,14 @@ public:
 
         // Cloud evolution functions
 
+            //- Set parcel thermo properties
+            void setParcelThermoProperties
+            (
+                parcelType& parcel,
+                const scalar lagrangianDt,
+                const bool fullyDescribed
+            );
+
             //- Check parcel properties
             void checkParcelProperties
             (

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C

@@ -233,6 +233,24 @@ Foam::ThermoCloud<CloudType>::~ThermoCloud()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+void Foam::ThermoCloud<CloudType>::setParcelThermoProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt,
+    const bool fullyDescribed
+)
+{
+    CloudType::setParcelThermoProperties(parcel, lagrangianDt, fullyDescribed);
+
+    if (!fullyDescribed)
+    {
+        parcel.T() = constProps_.T0();
+        parcel.Cp() = constProps_.Cp0();
+    }
+}
+
+
 template<class CloudType>
 void Foam::ThermoCloud<CloudType>::checkParcelProperties
 (
@@ -242,12 +260,6 @@ void Foam::ThermoCloud<CloudType>::checkParcelProperties
 )
 {
     CloudType::checkParcelProperties(parcel, lagrangianDt, fullyDescribed);
-
-    if (!fullyDescribed)
-    {
-        parcel.T() = constProps_.T0();
-        parcel.Cp() = constProps_.Cp0();
-    }
 }
 
 
@@ -304,6 +316,15 @@ void Foam::ThermoCloud<CloudType>::scaleSources()
 }
 
 
+template<class CloudType>
+void Foam::ThermoCloud<CloudType>::preEvolve()
+{
+    CloudType::preEvolve();
+
+    this->pAmbient() = thermo_.thermo().p().average().value();
+}
+
+
 template<class CloudType>
 void Foam::ThermoCloud<CloudType>::evolve()
 {

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H

@@ -281,6 +281,14 @@ public:
 
         // Cloud evolution functions
 
+            //- Set parcel thermo properties
+            void setParcelThermoProperties
+            (
+                parcelType& parcel,
+                const scalar lagrangianDt,
+                const bool fullyDescribed
+            );
+
             //- Check parcel properties
             void checkParcelProperties
             (
@@ -304,6 +312,9 @@ public:
             //- Apply scaling to (transient) cloud sources
             void scaleSources();
 
+            //- Pre-evolve
+            void preEvolve();
+
             //- Evolve the cloud
             void evolve();
 

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H

@@ -124,6 +124,17 @@ public:
                 const bool readFields = true
             );
 
+            //- Construct from components
+            constantProperties
+            (
+                const label parcelTypeId,
+                const scalar rhoMin,
+                const scalar rho0,
+                const scalar minParticleMass,
+                const scalar youngsModulus,
+                const scalar poissonsRatio
+            );
+
 
         // Member functions
 
@@ -536,6 +547,15 @@ public:
                 const scalar muc        // carrier dynamic viscosity
             ) const;
 
+            //- Weber number
+            inline scalar We
+            (
+                const vector& U,        // particle velocity
+                const scalar d,         // particle diameter
+                const scalar rhoc,      // carrier density
+                const scalar sigma      // particle surface tension
+            ) const;
+
 
         // Main calculation loop
 

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H

@@ -85,6 +85,26 @@ inline Foam::KinematicParcel<ParcelType>::constantProperties::constantProperties
     }
 }
 
+template<class ParcelType>
+inline Foam::KinematicParcel<ParcelType>::constantProperties::constantProperties
+(
+    const label parcelTypeId,
+    const scalar rhoMin,
+    const scalar rho0,
+    const scalar minParticleMass,
+    const scalar youngsModulus,
+    const scalar poissonsRatio
+)
+:
+    dict_(dictionary::null),
+    parcelTypeId_(parcelTypeId),
+    rhoMin_(rhoMin),
+    rho0_(rho0),
+    minParticleMass_(minParticleMass),
+    youngsModulus_(youngsModulus),
+    poissonsRatio_(poissonsRatio)
+{}
+
 
 template<class ParcelType>
 inline Foam::KinematicParcel<ParcelType>::KinematicParcel
@@ -522,4 +542,17 @@ inline Foam::scalar Foam::KinematicParcel<ParcelType>::Re
 }
 
 
+template<class ParcelType>
+inline Foam::scalar Foam::KinematicParcel<ParcelType>::We
+(
+    const vector& U,
+    const scalar d,
+    const scalar rhoc,
+    const scalar sigma
+) const
+{
+    return rhoc*magSqr(U - Uc_)*d/(sigma + ROOTVSMALL);
+}
+
+
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H

@@ -38,8 +38,8 @@ SourceFiles
 #ifndef ReactingMultiphaseParcel_H
 #define ReactingMultiphaseParcel_H
 
-#include "ReactingParcel.H"
+#include "particle.H"
-#include "ReactingMultiphaseCloud.H"
+#include "SLGThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H

@@ -26,8 +26,7 @@ License
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class ParcelType>
-inline
+inline Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::
-Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::
 constantProperties()
 :
     ParcelType::constantProperties(),
@@ -37,8 +36,7 @@ constantProperties()
 
 
 template<class ParcelType>
-inline
+inline Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::
-Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::
 constantProperties
 (
     const constantProperties& cp
@@ -227,4 +225,5 @@ inline bool& Foam::ReactingMultiphaseParcel<ParcelType>::canCombust()
     return canCombust_;
 }
 
+
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C

@@ -24,8 +24,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "ReactingParcel.H"
-#include "mathematicalConstants.H"
 #include "specie.H"
+#include "CompositionModel.H"
+#include "mathematicalConstants.H"
 
 using namespace Foam::constant::mathematical;
 

src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.H

@@ -38,8 +38,8 @@ SourceFiles
 #ifndef ReactingParcel_H
 #define ReactingParcel_H
 
-#include "ThermoParcel.H"
+#include "particle.H"
-#include "ReactingCloud.H"
+#include "SLGThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -104,6 +104,27 @@ public:
                 const bool readFields = true
             );
 
+            //- Construct from components
+            constantProperties
+            (
+                const label parcelTypeId,
+                const scalar rhoMin,
+                const scalar rho0,
+                const scalar minParticleMass,
+                const scalar youngsModulus,
+                const scalar poissonsRatio,
+                const scalar T0,
+                const scalar TMin,
+                const scalar Cp0,
+                const scalar epsilon0,
+                const scalar f0,
+                const scalar Pr,
+                const scalar pMin,
+                const Switch& constantVolume,
+                const scalar Tvap,
+                const scalar Tbp
+            );
+
 
         // Access
 
@@ -136,7 +157,6 @@ public:
                 autoPtr<interpolation<scalar> > pInterp_;
 
 
-
     public:
 
         typedef typename ParcelType::template TrackingData<CloudType>::trackPart
@@ -327,6 +347,9 @@ public:
             //- Return the owner cell pressure
             inline scalar pc() const;
 
+            //- Return reference to the owner cell pressure
+            inline scalar& pc();
+
 
         // Edit
 
@@ -403,7 +426,7 @@ public:
                 const CompositionType& compModel
             );
 
-            //- Read - composition supplied
+            //- Write - composition supplied
             template<class CloudType>
             static void writeFields(const CloudType& c);
 

src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H

@@ -74,6 +74,49 @@ inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
 }
 
 
+template<class ParcelType>
+inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
+(
+    const label parcelTypeId,
+    const scalar rhoMin,
+    const scalar rho0,
+    const scalar minParticleMass,
+    const scalar youngsModulus,
+    const scalar poissonsRatio,
+    const scalar T0,
+    const scalar TMin,
+    const scalar Cp0,
+    const scalar epsilon0,
+    const scalar f0,
+    const scalar Pr,
+    const scalar pMin,
+    const Switch& constantVolume,
+    const scalar Tvap,
+    const scalar Tbp
+)
+:
+    ParcelType::constantProperties
+    (
+        parcelTypeId,
+        rhoMin,
+        rho0,
+        minParticleMass,
+        youngsModulus,
+        poissonsRatio,
+        T0,
+        TMin,
+        Cp0,
+        epsilon0,
+        f0,
+        Pr
+    ),
+    pMin_(pMin),
+    constantVolume_(constantVolume),
+    Tvap_(Tvap),
+    Tbp_(Tbp)
+{}
+
+
 template<class ParcelType>
 inline Foam::ReactingParcel<ParcelType>::ReactingParcel
 (
@@ -194,6 +237,13 @@ inline Foam::scalar Foam::ReactingParcel<ParcelType>::pc() const
 }
 
 
+template<class ParcelType>
+inline Foam::scalar& Foam::ReactingParcel<ParcelType>::pc()
+{
+    return pc_;
+}
+
+
 template<class ParcelType>
 inline Foam::scalar& Foam::ReactingParcel<ParcelType>::mass0()
 {

src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H

@@ -39,8 +39,8 @@ SourceFiles
 #ifndef ThermoParcel_H
 #define ThermoParcel_H
 
-#include "KinematicParcel.H"
+#include "particle.H"
-#include "ThermoCloud.H"
+#include "SLGThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -112,6 +112,23 @@ public:
                 const bool readFields = true
             );
 
+            //- Construct from components
+            constantProperties
+            (
+                const label parcelTypeId,
+                const scalar rhoMin,
+                const scalar rho0,
+                const scalar minParticleMass,
+                const scalar youngsModulus,
+                const scalar poissonsRatio,
+                const scalar T0,
+                const scalar TMin,
+                const scalar Cp0,
+                const scalar epsilon0,
+                const scalar f0,
+                const scalar Pr
+            );
+
 
         // Member functions
 

src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelI.H

@@ -81,6 +81,41 @@ inline Foam::ThermoParcel<ParcelType>::constantProperties::constantProperties
 }
 
 
+template<class ParcelType>
+inline Foam::ThermoParcel<ParcelType>::constantProperties::constantProperties
+(
+    const label parcelTypeId,
+    const scalar rhoMin,
+    const scalar rho0,
+    const scalar minParticleMass,
+    const scalar youngsModulus,
+    const scalar poissonsRatio,
+    const scalar T0,
+    const scalar TMin,
+    const scalar Cp0,
+    const scalar epsilon0,
+    const scalar f0,
+    const scalar Pr
+)
+:
+    ParcelType::constantProperties
+    (
+        parcelTypeId,
+        rhoMin,
+        rho0,
+        minParticleMass,
+        youngsModulus,
+        poissonsRatio
+    ),
+    T0_(T0),
+    TMin_(TMin),
+    Cp0_(Cp0),
+    epsilon0_(epsilon0),
+    f0_(f0),
+    Pr_(Pr)
+{}
+
+
 template<class ParcelType>
 inline Foam::ThermoParcel<ParcelType>::ThermoParcel
 (

src/lagrangian/intermediate/parcels/derived/basicReactingParcel/makeBasicReactingParcelSubmodels.C

@@ -45,9 +45,6 @@ License
 
 namespace Foam
 {
-    typedef basicReactingCloud::thermoCloudType thermoCloudType2;
-    typedef basicReactingCloud::kinematicCloudType kinematicCloudType2;
-
     makeParcelCloudFunctionObjects(basicReactingCloud);
 
     // Kinematic sub-models

src/lagrangian/intermediate/parcels/include/makeParcelInjectionModels.H

@@ -30,6 +30,7 @@ License
 
 #include "CellZoneInjection.H"
 #include "ConeInjection.H"
+#include "ConeNozzleInjection.H"
 #include "FieldActivatedInjection.H"
 #include "InflationInjection.H"
 #include "KinematicLookupTableInjection.H"

src/lagrangian/intermediate/parcels/include/makeReactingMultiphaseParcelInjectionModels.H

@@ -30,6 +30,7 @@ License
 
 #include "CellZoneInjection.H"
 #include "ConeInjection.H"
+#include "ConeNozzleInjection.H"
 #include "FieldActivatedInjection.H"
 #include "ManualInjection.H"
 #include "NoInjection.H"

src/lagrangian/intermediate/parcels/include/makeReactingParcelInjectionModels.H

@@ -30,6 +30,7 @@ License
 
 #include "CellZoneInjection.H"
 #include "ConeInjection.H"
+#include "ConeNozzleInjection.H"
 #include "FieldActivatedInjection.H"
 #include "ManualInjection.H"
 #include "NoInjection.H"

src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.C

@@ -27,6 +27,7 @@ License
 #include "Pstream.H"
 #include "ListListOps.H"
 #include "surfaceWriter.H"
+#include "globalIndex.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/NoCollision/NoCollision.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -87,7 +87,7 @@ public:
         //  this model this will always be 1.
         virtual label nSubCycles() const;
 
-        //- Flag to indicate whether model activates injection model
+        //- Flag to indicate whether model activates collision model
         virtual bool active() const;
 
         //- Indicates whether model determines wall collisions or not,

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/CellZoneInjection/CellZoneInjection.C

@@ -158,45 +158,6 @@ void Foam::CellZoneInjection<CloudType>::setPositions
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::CellZoneInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((0.0 >= time0) && (0.0 < time1))
-    {
-        return positions_.size();
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::CellZoneInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    // All parcels introduced at SOI
-    if ((0.0 >= time0) && (0.0 < time1))
-    {
-        return this->volumeTotal_;
-    }
-    else
-    {
-        return 0.0;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -316,6 +277,43 @@ Foam::scalar Foam::CellZoneInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::CellZoneInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((0.0 >= time0) && (0.0 < time1))
+    {
+        return positions_.size();
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::CellZoneInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    // All parcels introduced at SOI
+    if ((0.0 >= time0) && (0.0 < time1))
+    {
+        return this->volumeTotal_;
+    }
+    else
+    {
+        return 0.0;
+    }
+}
+
+
 template<class CloudType>
 void Foam::CellZoneInjection<CloudType>::setPositionAndCell
 (

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/CellZoneInjection/CellZoneInjection.H

@@ -95,25 +95,6 @@ class CellZoneInjection
         void setPositions(const labelList& cellZoneCells);
 
 
-protected:
-
-    // Protected member functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        label parcelsToInject
-        (
-            const scalar time0,
-            const scalar time1
-        );
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        scalar volumeToInject
-        (
-            const scalar time0,
-            const scalar time1
-        );
-
-
 public:
 
     //- Runtime type information
@@ -147,6 +128,12 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        scalar volumeToInject(const scalar time0, const scalar time1);
+
 
         // Injection geometry
 

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjection/ConeInjection.C

@@ -30,53 +30,6 @@ License
 
 using namespace Foam::constant::mathematical;
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::ConeInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        const scalar targetVolume = flowRateProfile_().integrate(0, time1);
-
-        const label targetParcels =
-            parcelsPerInjector_*targetVolume/this->volumeTotal_;
-
-        const label nToInject = targetParcels - nInjected_;
-
-        nInjected_ += nToInject;
-
-        return positionAxis_.size()*nToInject;
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::ConeInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        return flowRateProfile_().integrate(time0, time1);
-    }
-    else
-    {
-        return 0.0;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -195,6 +148,51 @@ Foam::scalar Foam::ConeInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::ConeInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        const scalar targetVolume = flowRateProfile_().integrate(0, time1);
+
+        const label targetParcels =
+            parcelsPerInjector_*targetVolume/this->volumeTotal_;
+
+        const label nToInject = targetParcels - nInjected_;
+
+        nInjected_ += nToInject;
+
+        return positionAxis_.size()*nToInject;
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::ConeInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return flowRateProfile_().integrate(time0, time1);
+    }
+    else
+    {
+        return 0.0;
+    }
+}
+
+
 template<class CloudType>
 void Foam::ConeInjection<CloudType>::setPositionAndCell
 (

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjection/ConeInjection.H

@@ -113,17 +113,6 @@ class ConeInjection
             vectorList tanVec2_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -157,6 +146,13 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
+
 
         // Injection geometry
 

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.C

@@ -0,0 +1,374 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "ConeNozzleInjection.H"
+#include "DataEntry.H"
+#include "mathematicalConstants.H"
+#include "distributionModel.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    InjectionModel<CloudType>(dict, owner, typeName),
+    injectionMethod_(imPoint),
+    outerNozzleDiameter_
+    (
+        readScalar(this->coeffDict().lookup("outerNozzleDiameter"))
+    ),
+    innerNozzleDiameter_
+    (
+        readScalar(this->coeffDict().lookup("innerNozzleDiameter"))
+    ),
+    duration_(readScalar(this->coeffDict().lookup("duration"))),
+    position_(this->coeffDict().lookup("position")),
+    injectorCell_(-1),
+    tetFaceI_(-1),
+    tetPtI_(-1),
+    direction_(this->coeffDict().lookup("direction")),
+    parcelsPerSecond_
+    (
+        readScalar(this->coeffDict().lookup("parcelsPerSecond"))
+    ),
+    volumeFlowRate_
+    (
+        DataEntry<scalar>::New
+        (
+            "volumeFlowRate",
+            this->coeffDict()
+        )
+    ),
+    Cd_
+    (
+        DataEntry<scalar>::New
+        (
+            "Cd",
+            this->coeffDict()
+        )
+    ),
+    thetaInner_
+    (
+        DataEntry<scalar>::New
+        (
+            "thetaInner",
+            this->coeffDict()
+        )
+    ),
+    thetaOuter_
+    (
+        DataEntry<scalar>::New
+        (
+            "thetaOuter",
+            this->coeffDict()
+        )
+    ),
+    sizeDistribution_
+    (
+        distributionModels::distributionModel::New
+        (
+            this->coeffDict().subDict("sizeDistribution"),
+            owner.rndGen()
+        )
+    ),
+    tanVec1_(vector::zero),
+    tanVec2_(vector::zero),
+    normal_(vector::zero)
+{
+    if (innerNozzleDiameter_ >= outerNozzleDiameter_)
+    {
+        FatalErrorIn
+        (
+            "Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection"
+            "("
+                "const dictionary&, "
+                "CloudType&"
+            ")"
+        )<< "innerNozzleDiameter >= outerNozzleDiameter" << nl
+         << exit(FatalError);
+    }
+
+    word injectionMethodType = this->coeffDict().lookup("injectionMethod");
+
+    if (injectionMethodType == "disc")
+    {
+        injectionMethod_ = imDisc;
+    }
+    else if (injectionMethodType == "point")
+    {
+        injectionMethod_ = imPoint;
+
+        // Set/cache the injector cell
+        this->findCellAtPosition
+        (
+            injectorCell_,
+            tetFaceI_,
+            tetPtI_,
+            position_,
+            false
+        );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "Foam::InjectionModel<CloudType>::InjectionModel"
+            "("
+                "const dictionary&, "
+                "CloudType&"
+            ")"
+        )<< "injectionMethod must be either 'point' or 'disc'" << nl
+         << exit(FatalError);
+    }
+
+    cachedRandom& rndGen = this->owner().rndGen();
+
+    // Normalise direction vector
+    direction_ /= mag(direction_);
+
+    // Determine direction vectors tangential to direction
+    vector tangent = vector::zero;
+    scalar magTangent = 0.0;
+
+    while(magTangent < SMALL)
+    {
+        vector v = rndGen.sample01<vector>();
+
+        tangent = v - (v & direction_)*direction_;
+        magTangent = mag(tangent);
+    }
+
+    tanVec1_ = tangent/magTangent;
+    tanVec2_ = direction_^tanVec1_;
+
+    // Set total volume to inject
+    this->volumeTotal_ = volumeFlowRate_().integrate(0.0, duration_);
+}
+
+
+template<class CloudType>
+Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection
+(
+    const ConeNozzleInjection<CloudType>& im
+)
+:
+    InjectionModel<CloudType>(im),
+    injectionMethod_(im.injectionMethod_),
+    outerNozzleDiameter_(im.outerNozzleDiameter_),
+    innerNozzleDiameter_(im.innerNozzleDiameter_),
+    duration_(im.duration_),
+    position_(im.position_),
+    injectorCell_(im.injectorCell_),
+    direction_(im.direction_),
+    parcelsPerSecond_(im.parcelsPerSecond_),
+    volumeFlowRate_(im.volumeFlowRate_().clone().ptr()),
+    Cd_(im.Cd_().clone().ptr()),
+    thetaInner_(im.thetaInner_().clone().ptr()),
+    thetaOuter_(im.thetaOuter_().clone().ptr()),
+    sizeDistribution_(im.sizeDistribution_().clone().ptr()),
+    tanVec1_(im.tanVec1_),
+    tanVec2_(im.tanVec1_),
+    normal_(im.normal_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::ConeNozzleInjection<CloudType>::~ConeNozzleInjection()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::ConeNozzleInjection<CloudType>::timeEnd() const
+{
+    return this->SOI_ + duration_;
+}
+
+
+template<class CloudType>
+Foam::label Foam::ConeNozzleInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return floor((time1 - time0)*parcelsPerSecond_);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::ConeNozzleInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return volumeFlowRate_().integrate(time0, time1);
+    }
+    else
+    {
+        return 0.0;
+    }
+}
+
+
+template<class CloudType>
+void Foam::ConeNozzleInjection<CloudType>::setPositionAndCell
+(
+    const label,
+    const label,
+    const scalar,
+    vector& position,
+    label& cellOwner,
+    label& tetFaceI,
+    label& tetPtI
+)
+{
+    cachedRandom& rndGen = this->owner().rndGen();
+
+    scalar beta = mathematical::twoPi*rndGen.sample01<scalar>();
+    normal_ = tanVec1_*cos(beta) + tanVec2_*sin(beta);
+
+    switch (injectionMethod_)
+    {
+        case imPoint:
+        {
+            position = position_;
+            cellOwner = injectorCell_;
+            tetFaceI = tetFaceI_;
+            tetPtI = tetPtI_;
+
+            break;
+        }
+        case imDisc:
+        {
+            scalar frac = rndGen.sample01<scalar>();
+            scalar dr = outerNozzleDiameter_ - innerNozzleDiameter_;
+            scalar r = 0.5*(innerNozzleDiameter_ + frac*dr);
+            position = position_ + r*normal_;
+
+            this->findCellAtPosition
+            (
+                cellOwner,
+                tetFaceI,
+                tetPtI,
+                position,
+                false
+            );
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "void Foam::ConeNozzleInjection<CloudType>::setPositionAndCell"
+                "("
+                    "const label, "
+                    "const label, "
+                    "const scalar, "
+                    "vector&, "
+                    "label&"
+                ")"
+            )<< "Unknown injectionMethod type" << nl
+             << exit(FatalError);
+        }
+    }
+}
+
+
+template<class CloudType>
+void Foam::ConeNozzleInjection<CloudType>::setProperties
+(
+    const label parcelI,
+    const label,
+    const scalar time,
+    typename CloudType::parcelType& parcel
+)
+{
+    cachedRandom& rndGen = this->owner().rndGen();
+
+    // set particle velocity
+    const scalar deg2Rad = mathematical::pi/180.0;
+
+    scalar t = time - this->SOI_;
+    scalar ti = thetaInner_().value(t);
+    scalar to = thetaOuter_().value(t);
+    scalar coneAngle = rndGen.sample01<scalar>()*(to - ti) + ti;
+
+    coneAngle *= deg2Rad;
+    scalar alpha = sin(coneAngle);
+    scalar dcorr = cos(coneAngle);
+
+    vector normal = alpha*normal_;
+    vector dirVec = dcorr*direction_;
+    dirVec += normal;
+    dirVec /= mag(dirVec);
+
+    scalar Ao = 0.25*mathematical::pi*outerNozzleDiameter_*outerNozzleDiameter_;
+    scalar Ai = 0.25*mathematical::pi*innerNozzleDiameter_*innerNozzleDiameter_;
+    scalar massFlowRate =
+        this->massTotal()*volumeFlowRate_().value(t)/this->volumeTotal();
+
+    scalar Umag = massFlowRate/(parcel.rho()*Cd_().value(t)*(Ao - Ai));
+    parcel.U() = Umag*dirVec;
+
+    // set particle diameter
+    parcel.d() = sizeDistribution_->sample();
+}
+
+
+template<class CloudType>
+bool Foam::ConeNozzleInjection<CloudType>::fullyDescribed() const
+{
+    return false;
+}
+
+
+template<class CloudType>
+bool Foam::ConeNozzleInjection<CloudType>::validInjection(const label)
+{
+    return true;
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.H

@@ -0,0 +1,235 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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/>.
+
+Class
+    Foam::ConeNozzleInjection
+
+Description
+    Cone injection
+
+    - User specifies
+      - time of start of injection
+      - injector position
+      - direction (along injection axis)
+      - parcel flow rate
+      - discharge coefficient, Cd
+      - inner and outer cone angles
+
+    - Parcel diameters obtained by size distribution model
+
+    - Parcel velocity is calculated as:
+
+        U = V_dot/(A * Cd), where V_dot is the volume flow rate
+
+    Based on the old 'unitInjection' model
+
+SourceFiles
+    ConeNozzleInjection.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ConeNozzleInjection_H
+#define ConeNozzleInjection_H
+
+#include "InjectionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+template<class Type>
+class DataEntry;
+
+class distributionModel;
+
+/*---------------------------------------------------------------------------*\
+                     Class ConeNozzleInjection Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class ConeNozzleInjection
+:
+    public InjectionModel<CloudType>
+{
+public:
+
+    //- Injection method enumeration
+    enum injectionMethod
+    {
+        imPoint,
+        imDisc
+    };
+
+
+private:
+
+    // Private data
+
+        //- point/disc injection method
+        injectionMethod injectionMethod_;
+
+        //- Outer nozzle diameter [m]
+        const scalar outerNozzleDiameter_;
+
+        //- Inner nozzle diameter [m]
+        const scalar innerNozzleDiameter_;
+
+        //- Injection duration [s]
+        const scalar duration_;
+
+        //- Injector position [m]
+        vector position_;
+
+        //- Cell containing injector position []
+        label injectorCell_;
+
+        //- Index of tet face for injector cell
+        label tetFaceI_;
+
+        //- Index of tet point for injector cell
+        label tetPtI_;
+
+        //- Injector direction []
+        vector direction_;
+
+        //- Number of parcels to introduce per second []
+        const label parcelsPerSecond_;
+
+        //- Volume flow rate of parcels to introduce relative to SOI [m^3/s]
+        const autoPtr<DataEntry<scalar> > volumeFlowRate_;
+
+        //- Discharge coefficient, relative to SOI [m/s]
+        const autoPtr<DataEntry<scalar> > Cd_;
+
+        //- Inner cone angle relative to SOI [deg]
+        const autoPtr<DataEntry<scalar> > thetaInner_;
+
+        //- Outer cone angle relative to SOI [deg]
+        const autoPtr<DataEntry<scalar> > thetaOuter_;
+
+        //- Parcel size PDF model
+        const autoPtr<distributionModels::distributionModel> sizeDistribution_;
+
+
+        // Tangential vectors to the direction vector
+
+            //- First tangential vector
+            vector tanVec1_;
+
+            //- Second tangential vector
+            vector tanVec2_;
+
+            //- injection vector orthogonal to direction
+            vector normal_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("coneNozzleInjection");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        ConeNozzleInjection(const dictionary& dict, CloudType& owner);
+
+        //- Construct copy
+        ConeNozzleInjection(const ConeNozzleInjection<CloudType>& im);
+
+        //- Construct and return a clone
+        virtual autoPtr<InjectionModel<CloudType> > clone() const
+        {
+            return autoPtr<InjectionModel<CloudType> >
+            (
+                new ConeNozzleInjection<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~ConeNozzleInjection();
+
+
+    // Member Functions
+
+        //- Return the end-of-injection time
+        scalar timeEnd() const;
+
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
+
+        // Injection geometry
+
+            //- Set the injection position and owner cell
+            virtual void setPositionAndCell
+            (
+                const label parcelI,
+                const label nParcels,
+                const scalar time,
+                vector& position,
+                label& cellOwner,
+                label& tetFaceI,
+                label& tetPtI
+            );
+
+            //- Set the parcel properties
+            virtual void setProperties
+            (
+                const label parcelI,
+                const label nParcels,
+                const scalar time,
+                typename CloudType::parcelType& parcel
+            );
+
+            //- Flag to identify whether model fully describes the parcel
+            virtual bool fullyDescribed() const;
+
+            //- Return flag to identify whether or not injection of parcelI is
+            //  permitted
+            virtual bool validInjection(const label parcelI);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "ConeNozzleInjection.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/FieldActivatedInjection/FieldActivatedInjection.C

@@ -29,44 +29,6 @@ License
 
 using namespace Foam::constant::mathematical;
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::FieldActivatedInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if (sum(nParcelsInjected_) < nParcelsPerInjector_*positions_.size())
-    {
-        return positions_.size();
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::FieldActivatedInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if (sum(nParcelsInjected_) < nParcelsPerInjector_*positions_.size())
-    {
-        return this->volumeTotal_/nParcelsPerInjector_;
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -180,16 +142,45 @@ Foam::FieldActivatedInjection<CloudType>::~FieldActivatedInjection()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
-bool Foam::FieldActivatedInjection<CloudType>::active() const
+Foam::scalar Foam::FieldActivatedInjection<CloudType>::timeEnd() const
 {
-    return true;
+    return GREAT;
 }
 
 
 template<class CloudType>
-Foam::scalar Foam::FieldActivatedInjection<CloudType>::timeEnd() const
+Foam::label Foam::FieldActivatedInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
 {
-    return GREAT;
+    if (sum(nParcelsInjected_) < nParcelsPerInjector_*positions_.size())
+    {
+        return positions_.size();
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::FieldActivatedInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if (sum(nParcelsInjected_) < nParcelsPerInjector_*positions_.size())
+    {
+        return this->volumeTotal_/nParcelsPerInjector_;
+    }
+    else
+    {
+        return 0;
+    }
 }
 
 

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/FieldActivatedInjection/FieldActivatedInjection.H

@@ -115,17 +115,6 @@ class FieldActivatedInjection
                 sizeDistribution_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -156,12 +145,15 @@ public:
 
     // Member Functions
 
-        //- Flag to indicate whether model activates injection model
-        bool active() const;
-
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
 
         // Injection geometry
 

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InflationInjection/InflationInjection.C

@@ -31,7 +31,129 @@ License
 
 using namespace Foam::constant::mathematical;
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::InflationInjection<CloudType>::InflationInjection
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    InjectionModel<CloudType>(dict, owner, typeName),
+    generationSetName_(this->coeffDict().lookup("generationCellSet")),
+    inflationSetName_(this->coeffDict().lookup("inflationCellSet")),
+    generationCells_(),
+    inflationCells_(),
+    duration_(readScalar(this->coeffDict().lookup("duration"))),
+    flowRateProfile_
+    (
+        DataEntry<scalar>::New
+        (
+            "flowRateProfile",
+            this->coeffDict()
+        )
+    ),
+    growthRate_
+    (
+        DataEntry<scalar>::New
+        (
+            "growthRate",
+            this->coeffDict()
+        )
+    ),
+    newParticles_(),
+    volumeAccumulator_(0.0),
+    fraction_(1.0),
+    selfSeed_(this->coeffDict().lookupOrDefault("selfSeed", false)),
+    dSeed_(SMALL),
+    sizeDistribution_
+    (
+        distributionModels::distributionModel::New
+        (
+            this->coeffDict().subDict("sizeDistribution"),
+            owner.rndGen()
+        )
+    )
+{
+    if (selfSeed_)
+    {
+        dSeed_ = readScalar(this->coeffDict().lookup("dSeed"));
+    }
+
+    cellSet generationCells(this->owner().mesh(), generationSetName_);
+
+    generationCells_ = generationCells.toc();
+
+    cellSet inflationCells(this->owner().mesh(), inflationSetName_);
+
+    // Union of cellSets
+    inflationCells |= generationCells;
+
+    inflationCells_ = inflationCells.toc();
+
+    if (Pstream::parRun())
+    {
+        scalar generationVolume = 0.0;
+
+        forAll(generationCells_, gCI)
+        {
+            label cI = generationCells_[gCI];
+
+            generationVolume += this->owner().mesh().cellVolumes()[cI];
+        }
+
+        scalar totalGenerationVolume = generationVolume;
+
+        reduce(totalGenerationVolume, sumOp<scalar>());
+
+        fraction_ = generationVolume/totalGenerationVolume;
+    }
+
+    // Set total volume/mass to inject
+    this->volumeTotal_ = fraction_*flowRateProfile_().integrate(0.0, duration_);
+    this->massTotal_ *= fraction_;
+}
+
+
+template<class CloudType>
+Foam::InflationInjection<CloudType>::InflationInjection
+(
+    const Foam::InflationInjection<CloudType>& im
+)
+:
+    InjectionModel<CloudType>(im),
+    generationSetName_(im.generationSetName_),
+    inflationSetName_(im.inflationSetName_),
+    generationCells_(im.generationCells_),
+    inflationCells_(im.inflationCells_),
+    duration_(im.duration_),
+    flowRateProfile_(im.flowRateProfile_().clone().ptr()),
+    growthRate_(im.growthRate_().clone().ptr()),
+    newParticles_(im.newParticles_),
+    volumeAccumulator_(im.volumeAccumulator_),
+    fraction_(im.fraction_),
+    selfSeed_(im.selfSeed_),
+    dSeed_(im.dSeed_),
+    sizeDistribution_(im.sizeDistribution_().clone().ptr())
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::InflationInjection<CloudType>::~InflationInjection()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::InflationInjection<CloudType>::timeEnd() const
+{
+    return this->SOI_ + duration_;
+}
+
 
 template<class CloudType>
 Foam::label Foam::InflationInjection<CloudType>::parcelsToInject
@@ -305,130 +427,6 @@ Foam::scalar Foam::InflationInjection<CloudType>::volumeToInject
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::InflationInjection<CloudType>::InflationInjection
-(
-    const dictionary& dict,
-    CloudType& owner
-)
-:
-    InjectionModel<CloudType>(dict, owner, typeName),
-    generationSetName_(this->coeffDict().lookup("generationCellSet")),
-    inflationSetName_(this->coeffDict().lookup("inflationCellSet")),
-    generationCells_(),
-    inflationCells_(),
-    duration_(readScalar(this->coeffDict().lookup("duration"))),
-    flowRateProfile_
-    (
-        DataEntry<scalar>::New
-        (
-            "flowRateProfile",
-            this->coeffDict()
-        )
-    ),
-    growthRate_
-    (
-        DataEntry<scalar>::New
-        (
-            "growthRate",
-            this->coeffDict()
-        )
-    ),
-    newParticles_(),
-    volumeAccumulator_(0.0),
-    fraction_(1.0),
-    selfSeed_(this->coeffDict().lookupOrDefault("selfSeed", false)),
-    dSeed_(SMALL),
-    sizeDistribution_
-    (
-        distributionModels::distributionModel::New
-        (
-            this->coeffDict().subDict("sizeDistribution"),
-            owner.rndGen()
-        )
-    )
-{
-    if (selfSeed_)
-    {
-        dSeed_ = readScalar(this->coeffDict().lookup("dSeed"));
-    }
-
-    cellSet generationCells(this->owner().mesh(), generationSetName_);
-
-    generationCells_ = generationCells.toc();
-
-    cellSet inflationCells(this->owner().mesh(), inflationSetName_);
-
-    // Union of cellSets
-    inflationCells |= generationCells;
-
-    inflationCells_ = inflationCells.toc();
-
-    if (Pstream::parRun())
-    {
-        scalar generationVolume = 0.0;
-
-        forAll(generationCells_, gCI)
-        {
-            label cI = generationCells_[gCI];
-
-            generationVolume += this->owner().mesh().cellVolumes()[cI];
-        }
-
-        scalar totalGenerationVolume = generationVolume;
-
-        reduce(totalGenerationVolume, sumOp<scalar>());
-
-        fraction_ = generationVolume/totalGenerationVolume;
-    }
-
-    // Set total volume/mass to inject
-    this->volumeTotal_ = fraction_*flowRateProfile_().integrate(0.0, duration_);
-    this->massTotal_ *= fraction_;
-}
-
-
-template<class CloudType>
-Foam::InflationInjection<CloudType>::InflationInjection
-(
-    const Foam::InflationInjection<CloudType>& im
-)
-:
-    InjectionModel<CloudType>(im),
-    generationSetName_(im.generationSetName_),
-    inflationSetName_(im.inflationSetName_),
-    generationCells_(im.generationCells_),
-    inflationCells_(im.inflationCells_),
-    duration_(im.duration_),
-    flowRateProfile_(im.flowRateProfile_().clone().ptr()),
-    growthRate_(im.growthRate_().clone().ptr()),
-    newParticles_(im.newParticles_),
-    volumeAccumulator_(im.volumeAccumulator_),
-    fraction_(im.fraction_),
-    selfSeed_(im.selfSeed_),
-    dSeed_(im.dSeed_),
-    sizeDistribution_(im.sizeDistribution_().clone().ptr())
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::InflationInjection<CloudType>::~InflationInjection()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::scalar Foam::InflationInjection<CloudType>::timeEnd() const
-{
-    return this->SOI_ + duration_;
-}
-
-
 template<class CloudType>
 void Foam::InflationInjection<CloudType>::setPositionAndCell
 (

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InflationInjection/InflationInjection.H

@@ -111,17 +111,6 @@ class InflationInjection
         const autoPtr<distributionModels::distributionModel> sizeDistribution_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -155,6 +144,12 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
 
         // Injection geometry
 

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C

@@ -101,46 +101,6 @@ void Foam::InjectionModel<CloudType>::writeProps()
 }
 
 
-template<class CloudType>
-Foam::label Foam::InjectionModel<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    notImplemented
-    (
-        "Foam::label Foam::InjectionModel<CloudType>::parcelsToInject"
-        "("
-            "const scalar, "
-            "const scalar"
-        ") const"
-    );
-
-    return 0;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::InjectionModel<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    notImplemented
-    (
-        "Foam::scalar Foam::InjectionModel<CloudType>::volumeToInject"
-        "("
-            "const scalar, "
-            "const scalar"
-        ") const"
-    );
-
-    return 0.0;
-}
-
-
 template<class CloudType>
 bool Foam::InjectionModel<CloudType>::validInjection(const label parcelI)
 {
@@ -510,6 +470,54 @@ Foam::scalar Foam::InjectionModel<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::InjectionModel<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    notImplemented
+    (
+        "Foam::label Foam::InjectionModel<CloudType>::parcelsToInject"
+        "("
+            "const scalar, "
+            "const scalar"
+        ") const"
+    );
+
+    return 0;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::InjectionModel<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    notImplemented
+    (
+        "Foam::scalar Foam::InjectionModel<CloudType>::volumeToInject"
+        "("
+            "const scalar, "
+            "const scalar"
+        ") const"
+    );
+
+    return 0.0;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::InjectionModel<CloudType>::averageParcelMass()
+{
+    label nTotal = parcelsToInject(0.0, timeEnd() - timeStart());
+    return massTotal_/nTotal;
+}
+
+
 template<class CloudType>
 template<class TrackData>
 void Foam::InjectionModel<CloudType>::inject(TrackData& td)
@@ -522,6 +530,7 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
     const scalar time = this->owner().db().time().value();
     const scalar trackTime = this->owner().solution().trackTime();
     const polyMesh& mesh = this->owner().mesh();
+    typename TrackData::cloudType& cloud = td.cloud();
 
     // Prepare for next time step
     label parcelsAdded = 0;
@@ -583,11 +592,14 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
                     tetPtI
                 );
 
+                // Check/set new parcel thermo properties
+                cloud.setParcelThermoProperties(*pPtr, dt, fullyDescribed());
+
                 // Assign new parcel properties in injection model
                 setProperties(parcelI, newParcels, timeInj, *pPtr);
 
-                // Check new parcel properties
+                // Check/set new parcel injection properties
-                td.cloud().checkParcelProperties(*pPtr, dt, fullyDescribed());
+                cloud.checkParcelProperties(*pPtr, dt, fullyDescribed());
 
                 // Apply correction to velocity for 2-D cases
                 meshTools::constrainDirection
@@ -639,6 +651,7 @@ void Foam::InjectionModel<CloudType>::injectSteadyState
     }
 
     const polyMesh& mesh = this->owner().mesh();
+    typename TrackData::cloudType& cloud = td.cloud();
 
     // Reset counters
     time0_ = 0.0;
@@ -688,11 +701,14 @@ void Foam::InjectionModel<CloudType>::injectSteadyState
                 tetPtI
             );
 
+            // Check/set new parcel thermo properties
+            cloud.setParcelThermoProperties(*pPtr, 0.0, fullyDescribed());
+
             // Assign new parcel properties in injection model
             setProperties(parcelI, newParcels, 0.0, *pPtr);
 
-            // Check new parcel properties
+            // Check/set new parcel injection properties
-            td.cloud().checkParcelProperties(*pPtr, 0.0, fullyDescribed());
+            cloud.checkParcelProperties(*pPtr, 0.0, fullyDescribed());
 
             // Apply correction to velocity for 2-D cases
             meshTools::constrainDirection

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H

@@ -143,12 +143,6 @@ protected:
 
     // Protected Member Functions
 
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
         //- Additional flag to identify whether or not injection of parcelI is
         //  permitted
         virtual bool validInjection(const label parcelI);
@@ -266,6 +260,24 @@ public:
             //- Return the end-of-injection time
             virtual scalar timeEnd() const;
 
+            //- Number of parcels to introduce relative to SOI
+            virtual label parcelsToInject
+            (
+                const scalar time0,
+                const scalar time1
+            );
+
+            //- Volume of parcels to introduce relative to SOI
+            virtual scalar volumeToInject
+            (
+                const scalar time0,
+                const scalar time1
+            );
+
+            //- Return the average parcel mass over the injection period
+            scalar averageParcelMass();
+
+
             // Counters
 
                 //- Return the number of injections

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/KinematicLookupTableInjection.C

@@ -26,46 +26,6 @@ License
 #include "KinematicLookupTableInjection.H"
 #include "scalarIOList.H"
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::KinematicLookupTableInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        return floor(injectorCells_.size()*(time1 - time0)*parcelsPerSecond_);
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::KinematicLookupTableInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    scalar volume = 0.0;
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        forAll(injectors_, i)
-        {
-            volume += injectors_[i].mDot()/injectors_[i].rho()*(time1 - time0);
-        }
-    }
-
-    return volume;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -156,6 +116,44 @@ Foam::scalar Foam::KinematicLookupTableInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::KinematicLookupTableInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return floor(injectorCells_.size()*(time1 - time0)*parcelsPerSecond_);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::KinematicLookupTableInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    scalar volume = 0.0;
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        forAll(injectors_, i)
+        {
+            volume += injectors_[i].mDot()/injectors_[i].rho()*(time1 - time0);
+        }
+    }
+
+    return volume;
+}
+
+
 template<class CloudType>
 void Foam::KinematicLookupTableInjection<CloudType>::setPositionAndCell
 (

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/KinematicLookupTableInjection.H

@@ -91,17 +91,6 @@ class KinematicLookupTableInjection
         labelList injectorTetPts_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -138,6 +127,12 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
 
         // Injection geometry
 

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ManualInjection/ManualInjection.C

@@ -30,45 +30,6 @@ License
 
 using namespace Foam::constant::mathematical;
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::ManualInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((0.0 >= time0) && (0.0 < time1))
-    {
-        return positions_.size();
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::ManualInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    // All parcels introduced at SOI
-    if ((0.0 >= time0) && (0.0 < time1))
-    {
-        return this->volumeTotal_;
-    }
-    else
-    {
-        return 0.0;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -192,6 +153,43 @@ Foam::scalar Foam::ManualInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::ManualInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((0.0 >= time0) && (0.0 < time1))
+    {
+        return positions_.size();
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::ManualInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    // All parcels introduced at SOI
+    if ((0.0 >= time0) && (0.0 < time1))
+    {
+        return this->volumeTotal_;
+    }
+    else
+    {
+        return 0.0;
+    }
+}
+
+
 template<class CloudType>
 void Foam::ManualInjection<CloudType>::setPositionAndCell
 (

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ManualInjection/ManualInjection.H

@@ -85,17 +85,6 @@ class ManualInjection
         const autoPtr<distributionModels::distributionModel> sizeDistribution_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -129,6 +118,12 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
 
         // Injection geometry
 

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/NoInjection/NoInjection.C

@@ -26,30 +26,6 @@ License
 #include "NoInjection.H"
 #include "addToRunTimeSelectionTable.H"
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::NoInjection<CloudType>::parcelsToInject
-(
-    const scalar,
-    const scalar
-)
-{
-    return 0;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::NoInjection<CloudType>::volumeToInject
-(
-    const scalar,
-    const scalar
-)
-{
-    return 0.0;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -89,6 +65,28 @@ Foam::scalar Foam::NoInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::NoInjection<CloudType>::parcelsToInject
+(
+    const scalar,
+    const scalar
+)
+{
+    return 0;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::NoInjection<CloudType>::volumeToInject
+(
+    const scalar,
+    const scalar
+)
+{
+    return 0.0;
+}
+
+
 template<class CloudType>
 void Foam::NoInjection<CloudType>::setPositionAndCell
 (

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/NoInjection/NoInjection.H

@@ -51,16 +51,6 @@ class NoInjection
 :
     public InjectionModel<CloudType>
 {
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
 
 public:
 
@@ -98,6 +88,12 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
 
         // Injection geometry
 

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/PatchInjection.C

@@ -27,64 +27,6 @@ License
 #include "DataEntry.H"
 #include "distributionModel.H"
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::PatchInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        scalar nParcels =fraction_*(time1 - time0)*parcelsPerSecond_;
-
-        cachedRandom& rnd = this->owner().rndGen();
-
-        label nParcelsToInject = floor(nParcels);
-
-        // Inject an additional parcel with a probability based on the
-        // remainder after the floor function
-        if
-        (
-            nParcelsToInject > 0
-         && (
-               nParcels - scalar(nParcelsToInject)
-             > rnd.position(scalar(0), scalar(1))
-            )
-        )
-        {
-            ++nParcelsToInject;
-        }
-
-        return nParcelsToInject;
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::PatchInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        return fraction_*flowRateProfile_().integrate(time0, time1);
-    }
-    else
-    {
-        return 0.0;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -182,6 +124,62 @@ Foam::scalar Foam::PatchInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::PatchInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        scalar nParcels =fraction_*(time1 - time0)*parcelsPerSecond_;
+
+        cachedRandom& rnd = this->owner().rndGen();
+
+        label nParcelsToInject = floor(nParcels);
+
+        // Inject an additional parcel with a probability based on the
+        // remainder after the floor function
+        if
+        (
+            nParcelsToInject > 0
+         && (
+               nParcels - scalar(nParcelsToInject)
+             > rnd.position(scalar(0), scalar(1))
+            )
+        )
+        {
+            ++nParcelsToInject;
+        }
+
+        return nParcelsToInject;
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::PatchInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return fraction_*flowRateProfile_().integrate(time0, time1);
+    }
+    else
+    {
+        return 0.0;
+    }
+}
+
+
 template<class CloudType>
 void Foam::PatchInjection<CloudType>::setPositionAndCell
 (

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/PatchInjection.H

@@ -94,17 +94,6 @@ class PatchInjection
         scalar fraction_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -138,6 +127,12 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
 
         // Injection geometry
 

src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H

@@ -38,6 +38,7 @@ SourceFiles
 #ifndef CompositionModel_H
 #define CompositionModel_H
 
+#include "SubModelBase.H"
 #include "IOdictionary.H"
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"

src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.C

@@ -25,46 +25,6 @@ License
 
 #include "ReactingLookupTableInjection.H"
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::ReactingLookupTableInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        return floor(injectorCells_.size()*(time1 - time0)*parcelsPerSecond_);
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::ReactingLookupTableInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    scalar volume = 0.0;
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        forAll(injectors_, i)
-        {
-            volume += injectors_[i].mDot()/injectors_[i].rho()*(time1 - time0);
-        }
-    }
-
-    return volume;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -155,6 +115,44 @@ Foam::scalar Foam::ReactingLookupTableInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::ReactingLookupTableInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return floor(injectorCells_.size()*(time1 - time0)*parcelsPerSecond_);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::ReactingLookupTableInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    scalar volume = 0.0;
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        forAll(injectors_, i)
+        {
+            volume += injectors_[i].mDot()/injectors_[i].rho()*(time1 - time0);
+        }
+    }
+
+    return volume;
+}
+
+
 template<class CloudType>
 void Foam::ReactingLookupTableInjection<CloudType>::setPositionAndCell
 (

src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.H

@@ -94,17 +94,6 @@ class ReactingLookupTableInjection
         labelList injectorTetPts_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -145,6 +134,13 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
+
 
         // Injection geometry
 

src/lagrangian/intermediate/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.C

@@ -25,48 +25,6 @@ License
 
 #include "ReactingMultiphaseLookupTableInjection.H"
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label
-Foam::ReactingMultiphaseLookupTableInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        return floor(injectorCells_.size()*(time1 - time0)*parcelsPerSecond_);
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar
-Foam::ReactingMultiphaseLookupTableInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    scalar volume = 0.0;
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        forAll(injectors_, i)
-        {
-            volume += injectors_[i].mDot()/injectors_[i].rho()*(time1 - time0);
-        }
-    }
-
-    return volume;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -161,6 +119,46 @@ Foam::ReactingMultiphaseLookupTableInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label
+Foam::ReactingMultiphaseLookupTableInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return floor(injectorCells_.size()*(time1 - time0)*parcelsPerSecond_);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar
+Foam::ReactingMultiphaseLookupTableInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    scalar volume = 0.0;
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        forAll(injectors_, i)
+        {
+            volume += injectors_[i].mDot()/injectors_[i].rho()*(time1 - time0);
+        }
+    }
+
+    return volume;
+}
+
+
 template<class CloudType>
 void Foam::ReactingMultiphaseLookupTableInjection<CloudType>::setPositionAndCell
 (

src/lagrangian/intermediate/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.H

@@ -97,17 +97,6 @@ class ReactingMultiphaseLookupTableInjection
         labelList injectorTetPts_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -148,6 +137,12 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
 
         // Injection geometry
 

src/lagrangian/intermediate/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.C

@@ -26,46 +26,6 @@ License
 #include "ThermoLookupTableInjection.H"
 #include "scalarIOList.H"
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::ThermoLookupTableInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        return floor(injectorCells_.size()*(time1 - time0)*parcelsPerSecond_);
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::ThermoLookupTableInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    scalar volume = 0.0;
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        forAll(injectors_, i)
-        {
-            volume += injectors_[i].mDot()/injectors_[i].rho()*(time1 - time0);
-        }
-    }
-
-    return volume;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -156,6 +116,44 @@ Foam::scalar Foam::ThermoLookupTableInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::ThermoLookupTableInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return floor(injectorCells_.size()*(time1 - time0)*parcelsPerSecond_);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::ThermoLookupTableInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    scalar volume = 0.0;
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        forAll(injectors_, i)
+        {
+            volume += injectors_[i].mDot()/injectors_[i].rho()*(time1 - time0);
+        }
+    }
+
+    return volume;
+}
+
+
 template<class CloudType>
 void Foam::ThermoLookupTableInjection<CloudType>::setPositionAndCell
 (

src/lagrangian/intermediate/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.H

@@ -93,17 +93,6 @@ class ThermoLookupTableInjection
         labelList injectorTetPts_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -141,6 +130,13 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
+
 
         // Injection geometry
 

src/lagrangian/spray/Make/files

@@ -0,0 +1,7 @@
+clouds/baseClasses/sprayCloud/sprayCloud.C
+
+SPRAYPARCEL=parcels/derived/basicSprayParcel
+$(SPRAYPARCEL)/defineBasicSprayParcel.C
+$(SPRAYPARCEL)/makeBasicSprayParcelSubmodels.C
+
+LIB = $(FOAM_LIBBIN)/liblagrangianSpray

src/lagrangian/spray/Make/options

@@ -0,0 +1,46 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/LES/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+LIB_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -ldistributionModels \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -lliquidProperties \
+    -lliquidMixtureProperties \
+    -lsolidProperties \
+    -lsolidMixtureProperties \
+    -lspecie \
+    -lbasicThermophysicalModels \
+    -lreactionThermophysicalModels \
+    -lSLGThermo \
+    -lchemistryModel \
+    -lradiationModels \
+    -lODE \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -ldynamicFvMesh \
+    -lsurfaceFilmModels

src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.C

@@ -0,0 +1,371 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "SprayCloud.H"
+#include "AtomizationModel.H"
+#include "BreakupModel.H"
+#include "StochasticCollisionModel.H"
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::SprayCloud<CloudType>::setModels()
+{
+    atomizationModel_.reset
+    (
+        AtomizationModel<SprayCloud<CloudType> >::New
+        (
+            this->subModelProperties(),
+            *this
+        ).ptr()
+    );
+
+    breakupModel_.reset
+    (
+        BreakupModel<SprayCloud<CloudType> >::New
+        (
+            this->subModelProperties(),
+            *this
+        ).ptr()
+    );
+
+    stochasticCollisionModel_.reset
+    (
+        StochasticCollisionModel<SprayCloud<CloudType> >::New
+        (
+            this->subModelProperties(),
+            *this
+        ).ptr()
+    );
+}
+
+
+template<class CloudType>
+void Foam::SprayCloud<CloudType>::cloudReset
+(
+    SprayCloud<CloudType>& c
+)
+{
+    CloudType::cloudReset(c);
+
+    atomizationModel_.reset(c.atomizationModel_.ptr());
+    breakupModel_.reset(c.breakupModel_.ptr());
+    stochasticCollisionModel_.reset(c.stochasticCollisionModel_.ptr());
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::SprayCloud<CloudType>::SprayCloud
+(
+    const word& cloudName,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const dimensionedVector& g,
+    const SLGThermo& thermo,
+    bool readFields
+)
+:
+    CloudType(cloudName, rho, U, g, thermo, false),
+    sprayCloud(),
+    cloudCopyPtr_(NULL),
+    averageParcelMass_(0.0),
+    atomizationModel_(NULL),
+    breakupModel_(NULL),
+    stochasticCollisionModel_(NULL)
+{
+    if (this->solution().active())
+    {
+        setModels();
+
+        averageParcelMass_ = this->injection().averageParcelMass();
+
+        if (readFields)
+        {
+            parcelType::readFields(*this, this->composition());
+        }
+
+        Info << "Average parcel mass: " << averageParcelMass_ << endl;
+    }
+}
+
+
+template<class CloudType>
+Foam::SprayCloud<CloudType>::SprayCloud
+(
+    SprayCloud<CloudType>& c,
+    const word& name
+)
+:
+    CloudType(c, name),
+    sprayCloud(),
+    cloudCopyPtr_(NULL),
+    averageParcelMass_(c.averageParcelMass_),
+    atomizationModel_(c.atomizationModel_->clone()),
+    breakupModel_(c.breakupModel_->clone()),
+    stochasticCollisionModel_(c.stochasticCollisionModel_->clone())
+{}
+
+
+template<class CloudType>
+Foam::SprayCloud<CloudType>::SprayCloud
+(
+    const fvMesh& mesh,
+    const word& name,
+    const SprayCloud<CloudType>& c
+)
+:
+    CloudType(mesh, name, c),
+    sprayCloud(),
+    cloudCopyPtr_(NULL),
+    averageParcelMass_(0.0),
+    atomizationModel_(NULL),
+    breakupModel_(NULL),
+    stochasticCollisionModel_(NULL)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::SprayCloud<CloudType>::~SprayCloud()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::SprayCloud<CloudType>::setParcelThermoProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt,
+    const bool fullyDescribed
+)
+{
+    CloudType::setParcelThermoProperties(parcel, lagrangianDt, fullyDescribed);
+
+    if (!fullyDescribed)
+    {
+        const liquidMixtureProperties& liqMix = this->composition().liquids();
+
+        const scalarField& Y(parcel.Y());
+        scalarField X(liqMix.X(Y));
+
+        // override rho and Cp from constantProperties
+        parcel.Cp() = liqMix.Cp(parcel.pc(), parcel.T(), X);
+        parcel.rho() = liqMix.rho(parcel.pc(), parcel.T(), X);
+    }
+}
+
+
+template<class CloudType>
+void Foam::SprayCloud<CloudType>::checkParcelProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt,
+    const bool fullyDescribed
+)
+{
+    CloudType::checkParcelProperties(parcel, lagrangianDt, fullyDescribed);
+
+    // store the injection position and initial drop size
+    parcel.position0() = parcel.position();
+    parcel.d0() = parcel.d();
+
+    parcel.y() = breakup().y0();
+    parcel.yDot() = breakup().yDot0();
+
+    parcel.liquidCore() = atomization().initLiquidCore();
+}
+
+
+template<class CloudType>
+void Foam::SprayCloud<CloudType>::storeState()
+{
+    cloudCopyPtr_.reset
+    (
+        static_cast<SprayCloud<CloudType>*>
+        (
+            clone(this->name() + "Copy").ptr()
+        )
+    );
+}
+
+
+template<class CloudType>
+void Foam::SprayCloud<CloudType>::restoreState()
+{
+    cloudReset(cloudCopyPtr_());
+    cloudCopyPtr_.clear();
+}
+
+
+template<class CloudType>
+void Foam::SprayCloud<CloudType>::evolve()
+{
+    if (this->solution().canEvolve())
+    {
+        typename parcelType::template
+            TrackingData<SprayCloud<CloudType> > td(*this);
+
+        this->solve(td);
+    }
+}
+
+
+template<class CloudType>
+template<class TrackData>
+void Foam::SprayCloud<CloudType>::motion(TrackData& td)
+{
+    const scalar dt = this->solution().trackTime();
+
+    td.part() = TrackData::tpLinearTrack;
+    CloudType::move(td, dt);
+
+    this->updateCellOccupancy();
+
+    if (stochasticCollision().active())
+    {
+        const liquidMixtureProperties& liqMix = this->composition().liquids();
+
+        label i = 0;
+        forAllIter(typename SprayCloud<CloudType>, *this, iter)
+        {
+            label j = 0;
+            forAllIter(typename SprayCloud<CloudType>, *this, jter)
+            {
+                if (j > i)
+                {
+                    parcelType& p = iter();
+                    scalar Vi = this->mesh().V()[p.cell()];
+                    scalarField X1(liqMix.X(p.Y()));
+                    scalar sigma1 = liqMix.sigma(p.pc(), p.T(), X1);
+                    scalar mp = p.mass()*p.nParticle();
+
+                    parcelType& q = jter();
+                    scalar Vj = this->mesh().V()[q.cell()];
+                    scalarField X2(liqMix.X(q.Y()));
+                    scalar sigma2 = liqMix.sigma(q.pc(), q.T(), X2);
+                    scalar mq = q.mass()*q.nParticle();
+
+                    bool updateProperties = stochasticCollision().update
+                    (
+                        dt,
+                        this->rndGen(),
+                        p.position(),
+                        mp,
+                        p.d(),
+                        p.nParticle(),
+                        p.U(),
+                        p.rho(),
+                        p.T(),
+                        p.Y(),
+                        sigma1,
+                        p.cell(),
+                        Vi,
+                        q.position(),
+                        mq,
+                        q.d(),
+                        q.nParticle(),
+                        q.U(),
+                        q.rho(),
+                        q.T(),
+                        q.Y(),
+                        sigma2,
+                        q.cell(),
+                        Vj
+                    );
+
+                    // for coalescence we need to update the density and
+                    // the diameter cause of the temp/conc/mass-change
+                    if (updateProperties)
+                    {
+                        if (mp > VSMALL)
+                        {
+                            scalarField Xp(liqMix.X(p.Y()));
+                            p.rho() = liqMix.rho(p.pc(), p.T(), Xp);
+                            p.Cp() = liqMix.Cp(p.pc(), p.T(), Xp);
+                            p.d() =
+                                cbrt
+                                (
+                                    6.0*mp
+                                   /(
+                                        p.nParticle()
+                                       *p.rho()
+                                       *constant::mathematical::pi
+                                    )
+                                );
+                        }
+
+                        if (mq > VSMALL)
+                        {
+                            scalarField Xq(liqMix.X(q.Y()));
+                            q.rho() = liqMix.rho(q.pc(), q.T(), Xq);
+                            q.Cp() = liqMix.Cp(q.pc(), q.T(), Xq);
+                            q.d() =
+                                cbrt
+                                (
+                                    6.0*mq
+                                   /(
+                                        q.nParticle()
+                                       *q.rho()
+                                       *constant::mathematical::pi
+                                    )
+                                );
+                        }
+                    }
+                }
+                j++;
+            }
+            i++;
+        }
+
+        // remove coalesced particles (diameter set to 0)
+        forAllIter(typename SprayCloud<CloudType>, *this, iter)
+        {
+            parcelType& p = iter();
+            if (p.mass() < VSMALL)
+            {
+                deleteParticle(p);
+            }
+        }
+    }
+}
+
+
+template<class CloudType>
+void Foam::SprayCloud<CloudType>::info() const
+{
+    CloudType::info();
+    scalar d32 = 1.0e+6*this->Dij(3, 2);
+    scalar pen = this->penetration(0.95);
+
+    Info << "    D32 (mu)                        = " << d32 << endl;
+    Info << "    Liquid penetration 95% mass (m) = " << pen << endl;
+}
+
+
+// ************************************************************************* //

src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.H

@@ -0,0 +1,262 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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/>.
+
+Class
+    Foam::SprayCloud
+
+Description
+    Templated base class for spray cloud
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SprayCloud_H
+#define SprayCloud_H
+
+#include "sprayCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+template<class CloudType>
+class AtomizationModel;
+
+template<class CloudType>
+class BreakupModel;
+
+template<class CloudType>
+class StochasticCollisionModel;
+
+/*---------------------------------------------------------------------------*\
+                      Class SprayCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class SprayCloud
+:
+    public CloudType,
+    public sprayCloud
+{
+public:
+
+    // Public typedefs
+
+        //- Type of cloud this cloud was instantiated for
+        typedef CloudType cloudType;
+
+        //- Type of parcel the cloud was instantiated for
+        typedef typename CloudType::particleType parcelType;
+
+        //- Convenience typedef for this cloud type
+        typedef SprayCloud<CloudType> sprayCloudType;
+
+
+private:
+
+    // Private data
+
+        //- Cloud copy pointer
+        autoPtr<SprayCloud<CloudType> > cloudCopyPtr_;
+
+        //- Average parcel mass
+        scalar averageParcelMass_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        SprayCloud(const SprayCloud&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const SprayCloud&);
+
+
+protected:
+
+    // Protected data
+
+        // References to the cloud sub-models
+
+            //- Atomization model
+            autoPtr<AtomizationModel<SprayCloud<CloudType> > >
+                atomizationModel_;
+
+            //- Break-up model
+            autoPtr<BreakupModel<SprayCloud<CloudType> > > breakupModel_;
+
+            //- Collision model
+            autoPtr<StochasticCollisionModel<SprayCloud<CloudType> > >
+                stochasticCollisionModel_;
+
+
+    // Protected Member Functions
+
+        // Initialisation
+
+            //- Set cloud sub-models
+            void setModels();
+
+
+        // Cloud evolution functions
+
+            //- Reset state of cloud
+            void cloudReset(SprayCloud<CloudType>& c);
+
+
+public:
+
+    // Constructors
+
+        //- Construct given carrier gas fields
+        SprayCloud
+        (
+            const word& cloudName,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const dimensionedVector& g,
+            const SLGThermo& thermo,
+            bool readFields = true
+        );
+
+        //- Copy constructor with new name
+        SprayCloud(SprayCloud<CloudType>& c, const word& name);
+
+        //- Copy constructor with new name - creates bare cloud
+        SprayCloud
+        (
+            const fvMesh& mesh,
+            const word& name,
+            const SprayCloud<CloudType>& c
+        );
+
+
+        //- Construct and return clone based on (this) with new name
+        virtual autoPtr<Cloud<parcelType> > clone(const word& name)
+        {
+            return autoPtr<Cloud<parcelType> >
+            (
+                new SprayCloud(*this, name)
+            );
+        }
+
+        //- Construct and return bare clone based on (this) with new name
+        virtual autoPtr<Cloud<parcelType> > cloneBare(const word& name) const
+        {
+            return autoPtr<Cloud<parcelType> >
+            (
+                new SprayCloud(this->mesh(), name, *this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~SprayCloud();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return a reference to the cloud copy
+            inline const SprayCloud& cloudCopy() const;
+
+            //- Return const-access to the average parcel mass
+            inline scalar averageParcelMass() const;
+
+
+            // Sub-models
+
+                //- Return const-access to the atomization model
+                inline const AtomizationModel<SprayCloud<CloudType> >&
+                    atomization() const;
+
+                //- Return const-access to the breakup model
+                inline const BreakupModel<SprayCloud<CloudType> >&
+                    breakup() const;
+
+                //- Return const-access to the breakup model
+                inline const StochasticCollisionModel<SprayCloud<CloudType> >&
+                    stochasticCollision() const;
+
+
+        // Check
+
+            //- Print cloud information
+            void info() const;
+
+
+        // Cloud evolution functions
+
+            //- Set parcel thermo properties
+            void setParcelThermoProperties
+            (
+                parcelType& parcel,
+                const scalar lagrangianDt,
+                const bool fullyDescribed
+            );
+
+            //- Check parcel properties
+            void checkParcelProperties
+            (
+                parcelType& parcel,
+                const scalar lagrangianDt,
+                const bool fullyDescribed
+            );
+
+            //- Store the current cloud state
+            void storeState();
+
+            //- Reset the current cloud to the previously stored state
+            void restoreState();
+
+            //- Evolve the spray (inject, move)
+            void evolve();
+            //- Particle motion
+
+            template<class TrackData>
+            void motion(TrackData& td);
+};
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "SprayCloudI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "SprayCloud.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloudI.H

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+inline const Foam::SprayCloud<CloudType>&
+Foam::SprayCloud<CloudType>::cloudCopy() const
+{
+    return cloudCopyPtr_();
+}
+
+
+template<class CloudType>
+inline const Foam::AtomizationModel<Foam::SprayCloud<CloudType> >&
+Foam::SprayCloud<CloudType>::atomization() const
+{
+    return atomizationModel_;
+}
+
+
+template<class CloudType>
+inline const Foam::BreakupModel<Foam::SprayCloud<CloudType> >&
+Foam::SprayCloud<CloudType>::breakup() const
+{
+    return breakupModel_;
+}
+
+
+template<class CloudType>
+inline const Foam::StochasticCollisionModel<Foam::SprayCloud<CloudType> >&
+Foam::SprayCloud<CloudType>::stochasticCollision() const
+{
+    return stochasticCollisionModel_;
+}
+
+
+template<class CloudType>
+inline Foam::scalar Foam::SprayCloud<CloudType>::averageParcelMass() const
+{
+    return averageParcelMass_;
+}
+

src/lagrangian/spray/clouds/baseClasses/sprayCloud/sprayCloud.C

@@ -0,0 +1,48 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "sprayCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(sprayCloud, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sprayCloud::sprayCloud()
+{}
+
+
+// * * * * * * * * * * * * * * * *  Destructors  * * * * * * * * * * * * * * //
+
+Foam::sprayCloud::~sprayCloud()
+{}
+
+
+// ************************************************************************* //

src/lagrangian/spray/clouds/baseClasses/sprayCloud/sprayCloud.H

@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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/>.
+
+Class
+    Foam::sprayCloud
+
+Description
+    Virtual abstract base class for templated SprayCloud
+
+SourceFiles
+    sprayCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sprayCloud_H
+#define sprayCloud_H
+
+#include "typeInfo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class sprayCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+class sprayCloud
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        sprayCloud(const sprayCloud&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const sprayCloud&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("sprayCloud");
+
+    // Constructors
+
+        //- Null constructor
+        sprayCloud();
+
+
+    //- Destructor
+    virtual ~sprayCloud();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/clouds/derived/basicSprayCloud/basicSprayCloud.H

@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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/>.
+
+Class
+    Foam::basicSprayCloud
+
+Description
+    Cloud class to introduce reacting spray parcels
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef basicSprayCloud_H
+#define basicSprayCloud_H
+
+#include "Cloud.H"
+#include "KinematicCloud.H"
+#include "ThermoCloud.H"
+#include "ReactingCloud.H"
+#include "SprayCloud.H"
+#include "basicSprayParcel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef SprayCloud
+    <
+        ReactingCloud
+        <
+            ThermoCloud
+            <
+                KinematicCloud
+                <
+                    Cloud
+                    <
+                        basicSprayParcel
+                    >
+                >
+            >
+        >
+    > basicSprayCloud;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+

src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C

@@ -0,0 +1,490 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "SprayParcel.H"
+#include "CompositionModel.H"
+#include "AtomizationModel.H"
+
+// * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
+
+template<class ParcelType>
+template<class TrackData>
+void Foam::SprayParcel<ParcelType>::setCellValues
+(
+    TrackData& td,
+    const scalar dt,
+    const label cellI
+)
+{
+    ParcelType::setCellValues(td, dt, cellI);
+}
+
+
+template<class ParcelType>
+template<class TrackData>
+void Foam::SprayParcel<ParcelType>::cellValueSourceCorrection
+(
+    TrackData& td,
+    const scalar dt,
+    const label cellI
+)
+{
+    ParcelType::cellValueSourceCorrection(td, dt, cellI);
+}
+
+
+template<class ParcelType>
+template<class TrackData>
+void Foam::SprayParcel<ParcelType>::calc
+(
+    TrackData& td,
+    const scalar dt,
+    const label cellI
+)
+{
+    bool coupled = td.cloud().solution().coupled();
+
+    // check if parcel belongs to liquid core
+    if (liquidCore() > 0.5)
+    {
+        // liquid core parcels should not interact with the gas
+        if (td.cloud().solution().coupled())
+        {
+            td.cloud().solution().coupled() = false;
+        }
+    }
+
+    // store the parcel properties
+    const scalarField& Y(this->Y());
+    scalarField X(td.cloud().composition().liquids().X(Y));
+
+    scalar T0 = this->T();
+    this->Cp() = td.cloud().composition().liquids().Cp(this->pc_, T0, X);
+    scalar rho0 = td.cloud().composition().liquids().rho(this->pc_, T0, X);
+    this->rho() = rho0;
+
+    ParcelType::calc(td, dt, cellI);
+
+    if (td.keepParticle)
+    {
+        // update Cp, diameter and density due to change in temperature
+        // and/or composition
+        scalar T1 = this->T();
+        const scalarField& Y1(this->Y());
+        scalarField X1(td.cloud().composition().liquids().X(Y1));
+
+        this->Cp() = td.cloud().composition().liquids().Cp(this->pc_, T1, X1);
+
+        scalar rho1 = td.cloud().composition().liquids().rho(this->pc_, T1, X1);
+        this->rho() = rho1;
+        scalar d1 = this->d()*cbrt(rho0/rho1);
+        this->d() = d1;
+
+        if (liquidCore() > 0.5)
+        {
+            calcAtomization(td, dt, cellI);
+
+            // preserve the total mass/volume by increasing the number of
+            // particles in parcels due to breakup
+            scalar d2 = this->d();
+            this->nParticle() *= pow3(d1/d2);
+        }
+        else
+        {
+            calcBreakup(td, dt, cellI);
+        }
+    }
+
+    // restore coupled
+    td.cloud().solution().coupled() = coupled;
+}
+
+
+template<class ParcelType>
+template<class TrackData>
+void Foam::SprayParcel<ParcelType>::calcAtomization
+(
+    TrackData& td,
+    const scalar dt,
+    const label cellI
+)
+{
+    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
+    const CompositionModel<reactingCloudType>& composition =
+        td.cloud().composition();
+
+    typedef typename TrackData::cloudType::sprayCloudType sprayCloudType;
+    const AtomizationModel<sprayCloudType>& atomization =
+        td.cloud().atomization();
+
+
+    // cell state info is updated in ReactingParcel calc
+    const scalarField& Y(this->Y());
+    scalarField X(composition.liquids().X(Y));
+
+    scalar rho = composition.liquids().rho(this->pc(), this->T(), X);
+    scalar mu = composition.liquids().mu(this->pc(), this->T(), X);
+    scalar sigma = composition.liquids().sigma(this->pc(), this->T(), X);
+
+    // Average molecular weight of carrier mix - assumes perfect gas
+    scalar Wc = this->rhoc_*specie::RR*this->Tc()/this->pc();
+    scalar R = specie::RR/Wc;
+    scalar Tav = atomization.Taverage(this->T(), this->Tc());
+
+    // calculate average gas density based on average temperature
+    scalar rhoAv = this->pc()/(R*Tav);
+
+    scalar soi = td.cloud().injection().timeStart();
+    scalar currentTime = td.cloud().db().time().value();
+    const vector& pos = this->position();
+    const vector& injectionPos = this->position0();
+
+    // disregard the continous phase when calculating the relative velocity
+    // (in line with the deactivated coupled assumption)
+    scalar Urel = mag(this->U());
+
+    scalar t0 = max(0.0, currentTime - this->age() - soi);
+    scalar t1 = min(t0 + dt, td.cloud().injection().timeEnd() - soi);
+    // this should be the vol flow rate from when the parcel was injected
+    scalar volFlowRate = td.cloud().injection().volumeToInject(t0, t1)/dt;
+
+    scalar chi = 0.0;
+    if (atomization.calcChi())
+    {
+        chi = this->chi(td, X);
+    }
+
+    atomization.update
+    (
+        dt,
+        this->d(),
+        this->liquidCore(),
+        this->tc(),
+        rho,
+        mu,
+        sigma,
+        volFlowRate,
+        rhoAv,
+        Urel,
+        pos,
+        injectionPos,
+        td.cloud().pAmbient(),
+        chi,
+        td.cloud().rndGen()
+    );
+}
+
+
+template<class ParcelType>
+template<class TrackData>
+void Foam::SprayParcel<ParcelType>::calcBreakup
+(
+    TrackData& td,
+    const scalar dt,
+    const label cellI
+)
+{
+    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
+    const CompositionModel<reactingCloudType>& composition =
+        td.cloud().composition();
+
+    typedef typename TrackData::cloudType cloudType;
+    typedef typename cloudType::parcelType parcelType;
+    typedef typename cloudType::forceType forceType;
+
+    const parcelType& p = static_cast<const parcelType&>(*this);
+    const forceType& forces = td.cloud().forces();
+
+    if (td.cloud().breakup().solveOscillationEq())
+    {
+        solveTABEq(td, dt);
+    }
+
+    // cell state info is updated in ReactingParcel calc
+    const scalarField& Y(this->Y());
+    scalarField X(composition.liquids().X(Y));
+
+    scalar rho = composition.liquids().rho(this->pc(), this->T(), X);
+    scalar mu = composition.liquids().mu(this->pc(), this->T(), X);
+    scalar sigma = composition.liquids().sigma(this->pc(), this->T(), X);
+
+    // Average molecular weight of carrier mix - assumes perfect gas
+    scalar Wc = this->rhoc()*specie::RR*this->Tc()/this->pc();
+    scalar R = specie::RR/Wc;
+    scalar Tav = td.cloud().atomization().Taverage(this->T(), this->Tc());
+
+    // calculate average gas density based on average temperature
+    scalar rhoAv = this->pc()/(R*Tav);
+    scalar muAv = this->muc();
+    vector Urel = this->U() - this->Uc();
+    scalar Urmag = mag(Urel);
+    scalar Re = this->Re(this->U(), this->d(), rhoAv, muAv);
+
+    const scalar mass = p.mass();
+    const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, muAv);
+    const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, muAv);
+    scalar tMom = 1.0/(Fcp.Sp() + Fncp.Sp());
+
+    const vector g = td.cloud().g().value();
+
+    scalar massChild = 0.0;
+    scalar dChild = 0.0;
+    if
+    (
+        td.cloud().breakup().update
+        (
+            dt,
+            g,
+            this->d(),
+            this->tc(),
+            this->ms(),
+            this->nParticle(),
+            this->KHindex(),
+            this->y(),
+            this->yDot(),
+            this->d0(),
+            rho,
+            mu,
+            sigma,
+            this->U(),
+            rhoAv,
+            muAv,
+            Urel,
+            Urmag,
+            tMom,
+            td.cloud().averageParcelMass(),
+            dChild,
+            massChild,
+            td.cloud().rndGen()
+        )
+    )
+    {
+        scalar Re = rhoAv*Urmag*dChild/muAv;
+        this->mass0() -= massChild;
+
+        // Add child parcel as copy of parent
+        SprayParcel<ParcelType>* child = new SprayParcel<ParcelType>(*this);
+        child->mass0() = massChild;
+        child->d() = dChild;
+        child->nParticle() = massChild/rho*this->volume(dChild);
+
+        const forceSuSp Fcp =
+            forces.calcCoupled(*child, dt, massChild, Re, muAv);
+        const forceSuSp Fncp =
+            forces.calcNonCoupled(*child, dt, massChild, Re, muAv);
+
+        child->liquidCore() = 0.0;
+        child->KHindex() = 1.0;
+        child->y() = td.cloud().breakup().y0();
+        child->yDot() = td.cloud().breakup().yDot0();
+        child->tc() = -GREAT;
+        child->ms() = 0.0;
+        child->injector() = this->injector();
+        child->tMom() = 1.0/(Fcp.Sp() + Fncp.Sp());
+        child->user() = 0.0;
+        child->setCellValues(td, dt, cellI);
+
+        td.cloud().addParticle(child);
+    }
+}
+
+
+template<class ParcelType>
+template<class TrackData>
+Foam::scalar Foam::SprayParcel<ParcelType>::chi
+(
+    TrackData& td,
+    const scalarField& X
+) const
+{
+    // modifications to take account of the flash boiling on primary break-up
+
+    static label nIter = 200;
+
+    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
+    const CompositionModel<reactingCloudType>& composition =
+        td.cloud().composition();
+
+    scalar chi = 0.0;
+    scalar T0 = this->T();
+    scalar Tc0 = this->Tc();
+    scalar p0 = this->pc();
+    scalar pAmb = td.cloud().pAmbient();
+
+    scalar pv = composition.liquids().pv(p0, T0, X);
+
+    forAll(composition.liquids(), i)
+    {
+        if (pv >= 0.999*pAmb)
+        {
+            // liquid is boiling - calc boiling temperature
+
+            const liquidProperties& liq = composition.liquids().properties()[i];
+            scalar TBoil = T0;
+
+            for (label k=0; k<nIter; k++)
+            {
+                scalar pBoil = liq.pv(p0, TBoil);
+
+                if (pBoil > p0)
+                {
+                    TBoil = TBoil - (T0 - Tc0)/nIter;
+                }
+                else
+                {
+                    break;
+                }
+            }
+
+            scalar hl = liq.hl(pAmb, TBoil);
+            scalar iTp = liq.h(pAmb, T0) - liq.rho(pAmb, T0);
+            scalar iTb = liq.h(pAmb, TBoil) - pAmb/liq.rho(pAmb, TBoil);
+
+            chi += X[i]*(iTp - iTb)/hl;
+        }
+    }
+
+    chi = min(1.0, max(chi, 0.0));
+
+    return chi;
+}
+
+
+template<class ParcelType>
+template<class TrackData>
+void Foam::SprayParcel<ParcelType>::solveTABEq
+(
+    TrackData& td,
+    const scalar dt
+)
+{
+    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
+    const CompositionModel<reactingCloudType>& composition =
+        td.cloud().composition();
+
+    const scalar& TABCmu = td.cloud().breakup().TABCmu();
+    const scalar& TABWeCrit = td.cloud().breakup().TABWeCrit();
+    const scalar& TABComega = td.cloud().breakup().TABComega();
+
+    scalar r = 0.5*this->d();
+    scalar r2 = r*r;
+    scalar r3 = r*r2;
+
+    const scalarField& Y(this->Y());
+    scalarField X(composition.liquids().X(Y));
+
+    scalar rho = composition.liquids().rho(this->pc(), this->T(), X);
+    scalar mu = composition.liquids().mu(this->pc(), this->T(), X);
+    scalar sigma = composition.liquids().sigma(this->pc(), this->T(), X);
+
+    // inverse of characteristic viscous damping time
+    scalar rtd = 0.5*TABCmu*mu/(rho*r2);
+
+    // oscillation frequency (squared)
+    scalar omega2 = TABComega*sigma/(rho*r3) - rtd*rtd;
+
+    if(omega2 > 0)
+    {
+        scalar omega = sqrt(omega2);
+        scalar rhoc = this->rhoc();
+        scalar Wetmp = this->We(this->U(), r, rhoc, sigma)/TABWeCrit;
+
+        scalar y1 = this->y() - Wetmp;
+        scalar y2 = this->yDot()/omega;
+
+        // update distortion parameters
+        scalar c = cos(omega*dt);
+        scalar s = sin(omega*dt);
+        scalar e = exp(-rtd*dt);
+        y2 = (this->yDot() + y1*rtd)/omega;
+
+        this->y() = Wetmp + e*(y1*c + y2*s);
+        if (this->y() < 0)
+        {
+            this->y() = 0.0;
+            this->yDot() = 0.0;
+        }
+        else
+        {
+            this->yDot() = (Wetmp - this->y())*rtd + e*omega*(y2*c - y1*s);
+        }
+    }
+    else
+    {
+        // reset distortion parameters
+        this->y() = 0;
+        this->yDot() = 0;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class ParcelType>
+Foam::SprayParcel<ParcelType>::SprayParcel(const SprayParcel<ParcelType>& p)
+:
+    ParcelType(p),
+    d0_(p.d0_),
+    position0_(p.position0_),
+    liquidCore_(p.liquidCore_),
+    KHindex_(p.KHindex_),
+    y_(p.y_),
+    yDot_(p.yDot_),
+    tc_(p.tc_),
+    ms_(p.ms_),
+    injector_(p.injector_),
+    tMom_(p.tMom_),
+    user_(p.user_)
+{}
+
+
+template <class ParcelType>
+Foam::SprayParcel<ParcelType>::SprayParcel
+(
+    const SprayParcel<ParcelType>& p,
+    const polyMesh& mesh
+)
+:
+    ParcelType(p, mesh),
+    d0_(p.d0_),
+    position0_(p.position0_),
+    liquidCore_(p.liquidCore_),
+    KHindex_(p.KHindex_),
+    y_(p.y_),
+    yDot_(p.yDot_),
+    tc_(p.tc_),
+    ms_(p.ms_),
+    injector_(p.injector_),
+    tMom_(p.tMom_),
+    user_(p.user_)
+{}
+
+
+// * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //
+
+#include "SprayParcelIO.C"
+
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H

@@ -0,0 +1,421 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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/>.
+
+Class
+    Foam::SprayParcel
+
+Description
+    Reacing spray parcel, with added functionality for atomization and breakup
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SprayParcel_H
+#define SprayParcel_H
+
+#include "particle.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class ParcelType>
+class SprayParcel;
+
+template<class ParcelType>
+Ostream& operator<<
+(
+    Ostream&,
+    const SprayParcel<ParcelType>&
+);
+
+/*---------------------------------------------------------------------------*\
+                          Class SprayParcel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ParcelType>
+class SprayParcel
+:
+    public ParcelType
+{
+
+protected:
+
+    // Protected data
+
+        // Spray parcel properties
+
+            //- Initial droplet diameter
+            scalar d0_;
+
+            //- Injection position
+            vector position0_;
+
+            //- Part of liquid core ( >0.5=liquid, <0.5=droplet )
+            scalar liquidCore_;
+
+            //- Index for KH Breakup
+            scalar KHindex_;
+
+            //- Spherical deviation
+            scalar y_;
+
+            //- Rate of change of spherical deviation
+            scalar yDot_;
+
+            //- Characteristic time (used in atomization and/or breakup model)
+            scalar tc_;
+
+            //- Stripped parcel mass due to breakup
+            scalar ms_;
+
+            //- Injected from injector (needed e.g. for calculating distance
+            //  from injector)
+            scalar injector_;
+
+            //- Momentum relaxation time (needed for calculating parcel acc.)
+            scalar tMom_;
+
+            //- Passive scalar (extra variable to be defined by user)
+            scalar user_;
+
+
+public:
+
+    // Static data members
+
+        //- String representation of properties
+        static string propHeader;
+
+        //- Runtime type information
+        TypeName("SprayParcel");
+
+
+    // Constructors
+
+        //- Construct from owner, position, and cloud owner
+        //  Other properties initialised as null
+        inline SprayParcel
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI
+        );
+
+        //- Construct from components
+        inline SprayParcel
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI,
+            const label typeId,
+            const scalar nParticle0,
+            const scalar d0,
+            const scalar dTarget0,
+            const vector& U0,
+            const vector& f0,
+            const vector& angularMomentum0,
+            const vector& torque0,
+            const scalarField& Y0,
+            const scalar liquidCore,
+            const scalar KHindex,
+            const scalar y,
+            const scalar yDot,
+            const scalar tc,
+            const scalar ms,
+            const scalar injector,
+            const scalar tMom,
+            const scalar user,
+            const typename ParcelType::constantProperties& constProps
+        );
+
+        //- Construct from Istream
+        SprayParcel
+        (
+            const polyMesh& mesh,
+            Istream& is,
+            bool readFields = true
+        );
+
+        //- Construct as a copy
+        SprayParcel
+        (
+            const SprayParcel& p,
+            const polyMesh& mesh
+        );
+
+        //- Construct as a copy
+        SprayParcel(const SprayParcel& p);
+
+        //- Construct and return a (basic particle) clone
+        virtual autoPtr<particle> clone() const
+        {
+            return autoPtr<particle>(new SprayParcel<ParcelType>(*this));
+        }
+
+        //- Construct and return a (basic particle) clone
+        virtual autoPtr<particle> clone(const polyMesh& mesh) const
+        {
+            return autoPtr<particle>
+            (
+                new SprayParcel<ParcelType>(*this, mesh)
+            );
+        }
+
+        //- Factory class to read-construct particles used for
+        //  parallel transfer
+        class iNew
+        {
+            const polyMesh& mesh_;
+
+        public:
+
+            iNew(const polyMesh& mesh)
+            :
+                mesh_(mesh)
+            {}
+
+            autoPtr<SprayParcel<ParcelType> > operator()(Istream& is) const
+            {
+                return autoPtr<SprayParcel<ParcelType> >
+                (
+                    new SprayParcel<ParcelType>(mesh_, is, true)
+                );
+            }
+        };
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return const access to initial droplet diameter
+            inline scalar d0() const;
+
+            //- Return const access to initial droplet position
+            inline const vector& position0() const;
+
+            //- Return const access to liquid core
+            inline scalar liquidCore() const;
+
+            //- Return const access to Kelvin-Helmholtz breakup index
+            inline scalar KHindex() const;
+
+            //- Return const access to spherical deviation
+            inline scalar y() const;
+
+            //- Return const access to rate of change of spherical deviation
+            inline scalar yDot() const;
+
+            //- Return const access to atomization characteristic time
+            inline scalar tc() const;
+
+            //- Return const access to stripped parcel mass
+            inline scalar ms() const;
+
+            //- Return const access to injector id
+            inline scalar injector() const;
+
+            //- Return const access to momentum relaxation time
+            inline scalar tMom() const;
+
+            //- Return const access to passive user scalar
+            inline scalar user() const;
+
+
+        // Edit
+
+            //- Return access to initial droplet diameter
+            inline scalar& d0();
+
+            //- Return access to initial droplet position
+            inline vector& position0();
+
+            //- Return access to liquid core
+            inline scalar& liquidCore();
+
+            //- Return access to Kelvin-Helmholtz breakup index
+            inline scalar& KHindex();
+
+            //- Return access to spherical deviation
+            inline scalar& y();
+
+            //- Return access to rate of change of spherical deviation
+            inline scalar& yDot();
+
+            //- Return access to atomization characteristic time
+            inline scalar& tc();
+
+            //- Return access to stripped parcel mass
+            inline scalar& ms();
+
+            //- Return access to injector id
+            inline scalar& injector();
+
+            //- Return access to momentum relaxation time
+            inline scalar& tMom();
+
+            //- Return access to passive user scalar
+            inline scalar& user();
+
+
+        // Main calculation loop
+
+            //- Set cell values
+            template<class TrackData>
+            void setCellValues
+            (
+                TrackData& td,
+                const scalar dt,
+                const label cellI
+            );
+
+            //- Correct parcel properties according to atomization model
+            template<class TrackData>
+            void calcAtomization
+            (
+                TrackData& td,
+                const scalar dt,
+                const label cellI
+            );
+
+            //- Correct parcel properties according to breakup model
+            template<class TrackData>
+            void calcBreakup
+            (
+                TrackData& td,
+                const scalar dt,
+                const label cellI
+            );
+
+            //- Correct cell values using latest transfer information
+            template<class TrackData>
+            void cellValueSourceCorrection
+            (
+                TrackData& td,
+                const scalar dt,
+                const label cellI
+            );
+
+            //- Correct surface values due to emitted species
+            template<class TrackData>
+            void correctSurfaceValues
+            (
+                TrackData& td,
+                const label cellI,
+                const scalar T,
+                const scalarField& Cs,
+                scalar& rhos,
+                scalar& mus,
+                scalar& Pr,
+                scalar& kappa
+            );
+
+            //- Update parcel properties over the time interval
+            template<class TrackData>
+            void calc
+            (
+                TrackData& td,
+                const scalar dt,
+                const label cellI
+            );
+
+            //- Calculate the chi-factor for flash-boiling for the
+            //  atomization model
+            template<class TrackData>
+            scalar chi
+            (
+                TrackData& td,
+                const scalarField& X
+            ) const;
+
+            //- Solve the TAB equation
+            template<class TrackData>
+            void solveTABEq
+            (
+                TrackData& td,
+                const scalar dt
+            );
+
+
+        // I-O
+
+            //- Read
+            template<class CloudType, class CompositionType>
+            static void readFields
+            (
+                CloudType& c,
+                const CompositionType& compModel
+            );
+
+            //- Read - no composition
+            template<class CloudType>
+            static void readFields(CloudType& c);
+
+            //- Write
+            template<class CloudType, class CompositionType>
+            static void writeFields
+            (
+                const CloudType& c,
+                const CompositionType& compModel
+            );
+
+            //- Write - composition supplied
+            template<class CloudType>
+            static void writeFields(const CloudType& c);
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<< <ParcelType>
+        (
+            Ostream&,
+            const SprayParcel<ParcelType>&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "SprayParcelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "SprayParcel.C"
+#endif
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelI.H

@@ -0,0 +1,270 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParcelType>
+inline Foam::SprayParcel<ParcelType>::SprayParcel
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI
+)
+:
+    ParcelType(mesh, position, cellI, tetFaceI, tetPtI),
+    d0_(this->d()),
+    position0_(position),
+    liquidCore_(0.0),
+    KHindex_(0.0),
+    y_(0.0),
+    yDot_(0.0),
+    tc_(0.0),
+    ms_(0.0),
+    injector_(1.0),
+    tMom_(GREAT),
+    user_(0.0)
+{}
+
+
+template<class ParcelType>
+inline Foam::SprayParcel<ParcelType>::SprayParcel
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI,
+    const label typeId,
+    const scalar nParticle0,
+    const scalar d0,
+    const scalar dTarget0,
+    const vector& U0,
+    const vector& f0,
+    const vector& angularMomentum0,
+    const vector& torque0,
+    const scalarField& Y0,
+    const scalar liquidCore,
+    const scalar KHindex,
+    const scalar y,
+    const scalar yDot,
+    const scalar tc,
+    const scalar ms,
+    const scalar injector,
+    const scalar tMom,
+    const scalar user,
+    const typename ParcelType::constantProperties& constProps
+)
+:
+    ParcelType
+    (
+        mesh,
+        position,
+        cellI,
+        tetFaceI,
+        tetPtI,
+        typeId,
+        nParticle0,
+        d0,
+        dTarget0,
+        U0,
+        f0,
+        angularMomentum0,
+        torque0,
+        Y0,
+        constProps
+    ),
+    d0_(d0),
+    position0_(position),
+    liquidCore_(liquidCore),
+    KHindex_(KHindex),
+    y_(y),
+    yDot_(yDot),
+    tc_(tc),
+    ms_(ms),
+    injector_(injector),
+    tMom_(tMom),
+    user_(user)
+{}
+
+
+// * * * * * * * * * * SprayParcel Member Functions * * * * * * * * * * * * //
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::d0() const
+{
+    return d0_;
+}
+
+
+template<class ParcelType>
+inline const Foam::vector& Foam::SprayParcel<ParcelType>::position0() const
+{
+    return position0_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::liquidCore() const
+{
+    return liquidCore_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::KHindex() const
+{
+    return KHindex_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::y() const
+{
+    return y_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::yDot() const
+{
+    return yDot_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::tc() const
+{
+    return tc_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::ms() const
+{
+    return ms_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::injector() const
+{
+    return injector_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::tMom() const
+{
+    return tMom_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::user() const
+{
+    return user_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::d0()
+{
+    return d0_;
+}
+
+
+template<class ParcelType>
+inline Foam::vector& Foam::SprayParcel<ParcelType>::position0()
+{
+    return position0_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::liquidCore()
+{
+    return liquidCore_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::KHindex()
+{
+    return KHindex_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::y()
+{
+    return y_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::yDot()
+{
+    return yDot_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::tc()
+{
+    return tc_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::ms()
+{
+    return ms_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::injector()
+{
+    return injector_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::tMom()
+{
+    return tMom_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::user()
+{
+    return user_;
+}
+
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelIO.C

@@ -0,0 +1,337 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "SprayParcel.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template <class ParcelType>
+Foam::string Foam::SprayParcel<ParcelType>::propHeader =
+    ParcelType::propHeader
+  + " d0"
+  + " position0"
+  + " liquidCore"
+  + " KHindex"
+  + " y"
+  + " yDot"
+  + " tc"
+  + " ms"
+  + " injector"
+  + " tMom"
+  + " user";
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::SprayParcel<ParcelType>::SprayParcel
+(
+    const polyMesh& mesh,
+    Istream& is,
+    bool readFields
+)
+:
+    ParcelType(mesh, is, readFields),
+    d0_(0.0),
+    position0_(vector::zero),
+    liquidCore_(0.0),
+    KHindex_(0.0),
+    y_(0.0),
+    yDot_(0.0),
+    tc_(0.0),
+    ms_(0.0),
+    injector_(1.0),
+    tMom_(GREAT),
+    user_(0.0)
+{
+    if (readFields)
+    {
+
+        if (is.format() == IOstream::ASCII)
+        {
+            d0_ = readScalar(is);
+            is >> position0_;
+            liquidCore_ = readScalar(is);
+            KHindex_ = readScalar(is);
+            y_ = readScalar(is);
+            yDot_ = readScalar(is);
+            tc_ = readScalar(is);
+            ms_ = readScalar(is);
+            injector_ = readScalar(is);
+            tMom_ = readScalar(is);
+            user_ = readScalar(is);
+        }
+        else
+        {
+            is.read
+            (
+                reinterpret_cast<char*>(&d0_),
+                sizeof(d0_)
+              + sizeof(position0_)
+              + sizeof(liquidCore_)
+              + sizeof(KHindex_)
+              + sizeof(y_)
+              + sizeof(yDot_)
+              + sizeof(tc_)
+              + sizeof(ms_)
+              + sizeof(injector_)
+              + sizeof(tMom_)
+              + sizeof(user_)
+            );
+        }
+    }
+
+    // Check state of Istream
+    is.check
+    (
+        "SprayParcel<ParcelType>::SprayParcel"
+        "("
+            "const polyMesh, "
+            "Istream&, "
+            "bool"
+        ")"
+    );
+}
+
+
+template<class ParcelType>
+template<class CloudType>
+void Foam::SprayParcel<ParcelType>::readFields(CloudType& c)
+{
+    if (!c.size())
+    {
+        return;
+    }
+
+    ParcelType::readFields(c);
+}
+
+
+template<class ParcelType>
+template<class CloudType, class CompositionType>
+void Foam::SprayParcel<ParcelType>::readFields
+(
+    CloudType& c,
+    const CompositionType& compModel
+)
+{
+    if (!c.size())
+    {
+        return;
+    }
+
+    ParcelType::readFields(c, compModel);
+
+    IOField<scalar> d0(c.fieldIOobject("d0", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, d0);
+
+    IOField<vector> position0
+    (
+        c.fieldIOobject("position0", IOobject::MUST_READ)
+    );
+    c.checkFieldIOobject(c, position0);
+
+    IOField<scalar> liquidCore(c.fieldIOobject
+    (
+        "liquidCore", IOobject::MUST_READ)
+    );
+    c.checkFieldIOobject(c, liquidCore);
+
+    IOField<scalar> KHindex(c.fieldIOobject("KHindex", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, KHindex);
+
+    IOField<scalar> y(c.fieldIOobject("y", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, y);
+
+    IOField<scalar> yDot(c.fieldIOobject("yDot", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, yDot);
+
+    IOField<scalar> tc(c.fieldIOobject("tc", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, tc);
+
+    IOField<scalar> ms(c.fieldIOobject("ms", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, ms);
+
+    IOField<scalar> injector(c.fieldIOobject("injector", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, injector);
+
+    IOField<scalar> tMom(c.fieldIOobject("tMom", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, tMom);
+
+    IOField<scalar> user(c.fieldIOobject("user", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, user);
+
+    label i = 0;
+    forAllIter(typename Cloud<SprayParcel<ParcelType> >, c, iter)
+    {
+        SprayParcel<ParcelType>& p = iter();
+        p.d0_ = d0[i];
+        p.position0_ = position0[i];
+        p.liquidCore_ = liquidCore[i];
+        p.KHindex_ = KHindex[i];
+        p.y_ = y[i];
+        p.yDot_ = yDot[i];
+        p.tc_ = tc[i];
+        p.ms_ = ms[i];
+        p.injector_ = injector[i];
+        p.tMom_ = tMom[i];
+        p.user_ = user[i];
+        i++;
+    }
+}
+
+
+template<class ParcelType>
+template<class CloudType>
+void Foam::SprayParcel<ParcelType>::writeFields(const CloudType& c)
+{
+    ParcelType::writeFields(c);
+}
+
+
+template<class ParcelType>
+template<class CloudType, class CompositionType>
+void Foam::SprayParcel<ParcelType>::writeFields
+(
+    const CloudType& c,
+    const CompositionType& compModel
+)
+{
+    ParcelType::writeFields(c, compModel);
+
+    label np = c.size();
+
+    IOField<scalar> d0(c.fieldIOobject("d0", IOobject::NO_READ), np);
+    IOField<vector> position0
+    (
+        c.fieldIOobject("position0", IOobject::NO_READ),
+        np
+    );
+    IOField<scalar> liquidCore
+    (
+        c.fieldIOobject("liquidCore", IOobject::NO_READ),
+        np
+    );
+    IOField<scalar> KHindex(c.fieldIOobject("KHindex", IOobject::NO_READ), np);
+    IOField<scalar> y(c.fieldIOobject("y", IOobject::NO_READ), np);
+    IOField<scalar> yDot(c.fieldIOobject("yDot", IOobject::NO_READ), np);
+    IOField<scalar> tc(c.fieldIOobject("tc", IOobject::NO_READ), np);
+    IOField<scalar> ms(c.fieldIOobject("ms", IOobject::NO_READ), np);
+    IOField<scalar> injector
+    (
+        c.fieldIOobject("injector", IOobject::NO_READ),
+        np
+    );
+    IOField<scalar> tMom(c.fieldIOobject("tMom", IOobject::NO_READ), np);
+    IOField<scalar> user(c.fieldIOobject("user", IOobject::NO_READ), np);
+
+    label i = 0;
+    forAllConstIter(typename Cloud<SprayParcel<ParcelType> >, c, iter)
+    {
+        const SprayParcel<ParcelType>& p = iter();
+        d0[i] = p.d0_;
+        position0[i] = p.position0_;
+        liquidCore[i] = p.liquidCore_;
+        KHindex[i] = p.KHindex_;
+        y[i] = p.y_;
+        yDot[i] = p.yDot_;
+        tc[i] = p.tc_;
+        ms[i] = p.ms_;
+        injector[i] = p.injector_;
+        tMom[i] = p.tMom_;
+        user[i] = p.user_;
+        i++;
+    }
+
+    d0.write();
+    position0.write();
+    liquidCore.write();
+    KHindex.write();
+    y.write();
+    yDot.write();
+    tc.write();
+    ms.write();
+    injector.write();
+    tMom.write();
+    user.write();
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const SprayParcel<ParcelType>& p
+)
+{
+    if (os.format() == IOstream::ASCII)
+    {
+        os  << static_cast<const ParcelType&>(p)
+        << token::SPACE << p.d0()
+        << token::SPACE << p.position0()
+        << token::SPACE << p.liquidCore()
+        << token::SPACE << p.KHindex()
+        << token::SPACE << p.y()
+        << token::SPACE << p.yDot()
+        << token::SPACE << p.tc()
+        << token::SPACE << p.ms()
+        << token::SPACE << p.injector()
+        << token::SPACE << p.tMom()
+        << token::SPACE << p.user();
+    }
+    else
+    {
+        os  << static_cast<const ParcelType&>(p);
+        os.write
+        (
+            reinterpret_cast<const char*>(&p.d0_),
+            sizeof(p.d0())
+          + sizeof(p.position0())
+          + sizeof(p.liquidCore())
+          + sizeof(p.KHindex())
+          + sizeof(p.y())
+          + sizeof(p.yDot())
+          + sizeof(p.tc())
+          + sizeof(p.ms())
+          + sizeof(p.injector())
+          + sizeof(p.tMom())
+          + sizeof(p.user())
+        );
+    }
+
+    // Check state of Ostream
+    os.check
+    (
+        "Ostream& operator<<(Ostream&, const SprayParcel<ParcelType>&)"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/derived/basicSprayParcel/basicSprayParcel.H

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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/>.
+
+Class
+    Foam::basicSprayParcel
+
+Description
+    Definition of spray parcel
+
+SourceFiles
+    basicSprayParcel.C
+    basicSprayParcelIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef basicSprayParcel_H
+#define basicSprayParcel_H
+
+#include "contiguous.H"
+#include "particle.H"
+#include "KinematicParcel.H"
+#include "ThermoParcel.H"
+#include "ReactingParcel.H"
+#include "SprayParcel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef SprayParcel
+    <
+        ReactingParcel
+        <
+            ThermoParcel
+            <
+                KinematicParcel
+                <
+                    particle
+                >
+            >
+        >
+    > basicSprayParcel;
+
+    template<>
+    inline bool contiguous<basicSprayParcel>()
+    {
+        return false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/derived/basicSprayParcel/defineBasicSprayParcel.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "basicSprayParcel.H"
+#include "Cloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebug(basicSprayParcel, 0);
+    defineTemplateTypeNameAndDebug(Cloud<basicSprayParcel>, 0);
+};
+
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/derived/basicSprayParcel/makeBasicSprayParcelSubmodels.C

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "basicSprayCloud.H"
+
+#include "makeParcelCloudFunctionObjects.H"
+
+// Kinematic
+#include "makeThermoParcelForces.H" // thermo variant
+#include "makeParcelDispersionModels.H"
+#include "makeSprayParcelInjectionModels.H" // Spray variant
+#include "makeParcelPatchInteractionModels.H"
+
+// Thermodynamic
+#include "makeParcelHeatTransferModels.H"
+
+// Reacting
+#include "makeReactingParcelCompositionModels.H"
+#include "makeReactingParcelPhaseChangeModels.H"
+#include "makeReactingParcelSurfaceFilmModels.H"
+
+// Spray
+#include "makeSprayParcelAtomizationModels.H"
+#include "makeSprayParcelBreakupModels.H"
+#include "makeSprayParcelCollisionModels.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makeParcelCloudFunctionObjects(basicSprayCloud);
+
+    // Kinematic sub-models
+    makeThermoParcelForces(basicSprayCloud);
+    makeParcelDispersionModels(basicSprayCloud);
+    makeSprayParcelInjectionModels(basicSprayCloud);
+    makeParcelPatchInteractionModels(basicSprayCloud);
+
+    // Thermo sub-models
+    makeParcelHeatTransferModels(basicSprayCloud);
+
+    // Reacting sub-models
+    makeReactingParcelCompositionModels(basicSprayCloud);
+    makeReactingParcelPhaseChangeModels(basicSprayCloud);
+    makeReactingParcelSurfaceFilmModels(basicSprayCloud);
+
+    // Spray sub-models
+    makeSprayParcelAtomizationModels(basicSprayCloud);
+    makeSprayParcelBreakupModels(basicSprayCloud);
+    makeSprayParcelCollisionModels(basicSprayCloud);
+};
+
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/include/makeSprayParcelAtomizationModels.H

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef makeSprayParcelAtomizationModels_H
+#define makeSprayParcelAtomizationModels_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "NoAtomization.H"
+#include "BlobsSheetAtomization.H"
+#include "LISAAtomization.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeSprayParcelAtomizationModels(CloudType)                           \
+                                                                              \
+    makeAtomizationModel(CloudType);                                          \
+    makeAtomizationModelType(NoAtomization, CloudType);                       \
+    makeAtomizationModelType(BlobsSheetAtomization, CloudType);               \
+    makeAtomizationModelType(LISAAtomization, CloudType);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/include/makeSprayParcelBreakupModels.H

@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef makeSprayParcelBreakupModels_H
+#define makeSprayParcelBreakupModels_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "NoBreakup.H"
+#include "PilchErdman.H"
+#include "ReitzDiwakar.H"
+#include "ReitzKHRT.H"
+#include "TAB.H"
+#include "ETAB.H"
+#include "SHF.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeSprayParcelBreakupModels(CloudType)                               \
+                                                                              \
+    makeBreakupModel(CloudType)                                               \
+    makeBreakupModelType(NoBreakup, CloudType);                               \
+    makeBreakupModelType(PilchErdman, CloudType);                             \
+    makeBreakupModelType(ReitzDiwakar, CloudType);                            \
+    makeBreakupModelType(ReitzKHRT, CloudType);                               \
+    makeBreakupModelType(TAB, CloudType);                                     \
+    makeBreakupModelType(ETAB, CloudType);                                    \
+    makeBreakupModelType(SHF, CloudType);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/include/makeSprayParcelCollisionModels.H

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef makeSprayParcelCollisionModels_H
+#define makeSprayParcelCollisionModels_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "NoStochasticCollision.H"
+#include "ORourkeCollision.H"
+#include "TrajectoryCollision.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeSprayParcelCollisionModels(CloudType)                             \
+                                                                              \
+    makeStochasticCollisionModel(CloudType);                                  \
+    makeStochasticCollisionModelType(NoStochasticCollision, CloudType);       \
+    makeStochasticCollisionModelType(ORourkeCollision, CloudType);            \
+    makeStochasticCollisionModelType(TrajectoryCollision, CloudType);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/include/makeSprayParcelInjectionModels.H

@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef makeSprayParcelInjectionModels_H
+#define makeSprayParcelInjectionModels_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "CellZoneInjection.H"
+#include "ConeInjection.H"
+#include "ConeNozzleInjection.H"
+#include "FieldActivatedInjection.H"
+#include "InflationInjection.H"
+#include "ManualInjection.H"
+#include "NoInjection.H"
+#include "PatchInjection.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeSprayParcelInjectionModels(CloudType)                             \
+                                                                              \
+    makeInjectionModel(CloudType);                                            \
+                                                                              \
+    makeInjectionModelType(CellZoneInjection, CloudType);                     \
+    makeInjectionModelType(ConeInjection, CloudType);                         \
+    makeInjectionModelType(ConeNozzleInjection, CloudType);                   \
+    makeInjectionModelType(FieldActivatedInjection, CloudType);               \
+    makeInjectionModelType(InflationInjection, CloudType);                    \
+    makeInjectionModelType(ManualInjection, CloudType);                       \
+    makeInjectionModelType(NoInjection, CloudType);                           \
+    makeInjectionModelType(PatchInjection, CloudType);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/AtomizationModel/AtomizationModel/AtomizationModel.C

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "AtomizationModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::AtomizationModel<CloudType>::AtomizationModel
+(
+    CloudType& owner
+)
+:
+    SubModelBase<CloudType>(owner)
+{}
+
+
+template<class CloudType>
+Foam::AtomizationModel<CloudType>::AtomizationModel
+(
+    const AtomizationModel<CloudType>& am
+)
+:
+    SubModelBase<CloudType>(am)
+{}
+
+
+template<class CloudType>
+Foam::AtomizationModel<CloudType>::AtomizationModel
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& type
+)
+:
+    SubModelBase<CloudType>(owner, dict, type)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::AtomizationModel<CloudType>::~AtomizationModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::AtomizationModel<CloudType>::initLiquidCore() const
+{
+    notImplemented
+    (
+        "Foam::scalar "
+        "Foam::AtomizationModel<CloudType>::initLiquidCore() const"
+    );
+
+    return 0.0;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::AtomizationModel<CloudType>::Taverage
+(
+    const scalar& Tl,
+    const scalar& Tc
+) const
+{
+    return (2.0*Tl + Tc)/3.0;
+}
+
+
+template<class CloudType>
+bool Foam::AtomizationModel<CloudType>::calcChi() const
+{
+    notImplemented("bool Foam::AtomizationModel<CloudType>::calcChi()");
+
+    return false;
+}
+
+
+template<class CloudType>
+void Foam::AtomizationModel<CloudType>::update
+(
+    const scalar dt,
+    scalar& d,
+    scalar& liquidCore,
+    scalar& tc,
+    const scalar rho,
+    const scalar mu,
+    const scalar sigma,
+    const scalar volFlowRate,
+    const scalar rhoAv,
+    const scalar Urel,
+    const vector& pos,
+    const vector& injectionPos,
+    const scalar pAmbient,
+    const scalar chi,
+    cachedRandom& rndGen
+) const
+{
+    notImplemented
+    (
+        "void Foam::AtomizationModel<CloudType>::update"
+        "("
+            "const scalar, "
+            "scalar&, "
+            "scalar&, "
+            "scalar&, "
+            "const scalar, "
+            "const scalar, "
+            "const scalar, "
+            "const scalar, "
+            "const scalar, "
+            "const scalar, "
+            "const vector&, "
+            "const vector&, "
+            "const scalar, "
+            "const scalar, "
+            "cachedRandom&"
+        ") const"
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "AtomizationModelNew.C"
+
+// ************************************************************************* //
+

src/lagrangian/spray/submodels/AtomizationModel/AtomizationModel/AtomizationModel.H

@@ -0,0 +1,188 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::AtomizationModel
+
+Description
+    Templated atomization model class
+
+SourceFiles
+    AtomizationModel.C
+    AtomizationModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef AtomizationModel_H
+#define AtomizationModel_H
+
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "SubModelBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class AtomizationModel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class AtomizationModel
+:
+    public SubModelBase<CloudType>
+{
+public:
+
+    //- Runtime type information
+    TypeName("atomizationModel");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        AtomizationModel,
+        dictionary,
+        (
+            const dictionary& dict,
+            CloudType& owner
+        ),
+        (dict, owner)
+    );
+
+
+    // Constructors
+
+        //- Construct null from owner
+        AtomizationModel(CloudType& owner);
+
+        //- Construct from dictionary
+        AtomizationModel
+        (
+            const dictionary& dict,
+            CloudType& owner,
+            const word& type
+        );
+
+        //- Construct copy
+        AtomizationModel(const AtomizationModel<CloudType>& am);
+
+        //- Construct and return a clone
+        virtual autoPtr<AtomizationModel<CloudType> > clone() const
+        {
+            return autoPtr<AtomizationModel<CloudType> >
+            (
+                new AtomizationModel<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~AtomizationModel();
+
+
+    //- Selector
+    static autoPtr<AtomizationModel<CloudType> > New
+    (
+        const dictionary& dict,
+        CloudType& owner
+    );
+
+
+    // Member Functions
+
+        //- initial value of liquidCore
+        virtual scalar initLiquidCore() const;
+
+        //- Average temperature calculation
+        scalar Taverage(const scalar& Tliq, const scalar& Tc) const;
+
+        //- flag to indicate if chi needs to be calculated
+        virtual bool calcChi() const;
+
+        virtual void update
+        (
+            const scalar dt,
+            scalar& d,
+            scalar& liquidCore,
+            scalar& tc,
+            const scalar rho,
+            const scalar mu,
+            const scalar sigma,
+            const scalar volFlowRate,
+            const scalar rhoAv,
+            const scalar Urel,
+            const vector& pos,
+            const vector& injectionPos,
+            const scalar pAmbient,
+            const scalar chi,
+            cachedRandom& rndGen
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeAtomizationModel(CloudType)                                       \
+                                                                              \
+    typedef CloudType::sprayCloudType sprayCloudType;                         \
+    defineNamedTemplateTypeNameAndDebug                                       \
+    (                                                                         \
+        AtomizationModel<sprayCloudType>,                                     \
+        0                                                                     \
+    );                                                                        \
+    defineTemplateRunTimeSelectionTable                                       \
+    (                                                                         \
+        AtomizationModel<sprayCloudType>,                                     \
+        dictionary                                                            \
+    );
+
+
+#define makeAtomizationModelType(SS, CloudType)                               \
+                                                                              \
+    typedef CloudType::sprayCloudType sprayCloudType;                         \
+    defineNamedTemplateTypeNameAndDebug(SS<sprayCloudType>, 0);               \
+                                                                              \
+    AtomizationModel<sprayCloudType>::                                        \
+        adddictionaryConstructorToTable<SS<sprayCloudType> >                  \
+            add##SS##CloudType##sprayCloudType##ConstructorToTable_;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "AtomizationModel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/AtomizationModel/AtomizationModel/AtomizationModelNew.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "AtomizationModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::autoPtr<Foam::AtomizationModel<CloudType> >
+Foam::AtomizationModel<CloudType>::New
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+{
+    word AtomizationModelType(dict.lookup("atomizationModel"));
+
+    Info<< "Selecting AtomizationModel " << AtomizationModelType << endl;
+
+    typename dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(AtomizationModelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "AtomizationModel<CloudType>::New"
+            "("
+                "const dictionary&, "
+                "CloudType&"
+            ")"
+        )   << "Unknown AtomizationModelType type "
+            << AtomizationModelType
+            << ", constructor not in hash table" << nl << nl
+            << "    Valid AtomizationModel types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc() << exit(FatalError);
+    }
+
+    return autoPtr<AtomizationModel<CloudType> >(cstrIter()(dict, owner));
+}
+
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.C

@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "BlobsSheetAtomization.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::BlobsSheetAtomization<CloudType>::BlobsSheetAtomization
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    AtomizationModel<CloudType>(dict, owner, typeName),
+    B_(readScalar(this->coeffDict().lookup("B"))),
+    angle_(readScalar(this->coeffDict().lookup("angle")))
+{}
+
+
+template <class CloudType>
+Foam::BlobsSheetAtomization<CloudType>::BlobsSheetAtomization
+(
+    const BlobsSheetAtomization<CloudType>& am
+)
+:
+    AtomizationModel<CloudType>(am),
+    B_(am.B_),
+    angle_(am.B_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::BlobsSheetAtomization<CloudType>::~BlobsSheetAtomization()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::BlobsSheetAtomization<CloudType>::initLiquidCore() const
+{
+    return 1.0;
+}
+
+
+template<class CloudType>
+bool Foam::BlobsSheetAtomization<CloudType>::calcChi() const
+{
+    return false;
+}
+
+
+template<class CloudType>
+void Foam::BlobsSheetAtomization<CloudType>::update
+(
+    const scalar dt,
+    scalar& d,
+    scalar& liquidCore,
+    scalar& tc,
+    const scalar rho,
+    const scalar mu,
+    const scalar sigma,
+    const scalar volFlowRate,
+    const scalar rhoAv,
+    const scalar Urel,
+    const vector& pos,
+    const vector& injectionPos,
+    const scalar pAmbient,
+    const scalar chi,
+    cachedRandom& rndGen
+) const
+{
+    scalar lBU =
+        B_
+      * sqrt
+        (
+            rho*sigma*d*cos(angle_*constant::mathematical::pi/360.0)
+          / sqr(rhoAv*Urel)
+        );
+
+    scalar pWalk = mag(pos - injectionPos);
+
+    if (pWalk > lBU)
+    {
+        liquidCore = 0.0;
+    }
+}
+
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.H

@@ -0,0 +1,141 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::BlobsSheetAtomization
+
+Description
+    Primary Breakup Model for pressure swirl atomizers.
+
+    Accurate description in
+    @verbatim
+    Z. Han, S. Parrish, P.V. Farrell, R.D. Reitz
+    "Modeling Atomization Processes Of Pressure Swirl Hollow-Cone Fuel Sprays"
+    Atomization and Sprays, vol. 7, pp. 663-684, 1997
+
+    and
+
+    L. Allocca, G. Bella, A. De Vita, L. Di Angelo
+    "Experimental Validation of a GDI Spray Model"
+    SAE Technical Paper Series, 2002-01-1137
+    @endverbatim
+
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BlobsSheetAtomization_H
+#define BlobsSheetAtomization_H
+
+#include "AtomizationModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                   Class BlobsSheetAtomization Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class BlobsSheetAtomization
+:
+    public AtomizationModel<CloudType>
+{
+private:
+
+    scalar B_;
+    scalar angle_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("blobsSheetAtomization");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        BlobsSheetAtomization(const dictionary& dict, CloudType& cloud);
+
+        //- Construct copy
+        BlobsSheetAtomization(const BlobsSheetAtomization<CloudType>& am);
+
+        //- Construct and return a clone
+        virtual autoPtr<AtomizationModel<CloudType> > clone() const
+        {
+            return autoPtr<AtomizationModel<CloudType> >
+            (
+                new BlobsSheetAtomization<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~BlobsSheetAtomization();
+
+
+    // Member Functions
+
+        //- initial value of liquidCore
+        virtual scalar initLiquidCore() const;
+
+        //- flag to indicate if chi needs to be calculated
+        virtual bool calcChi() const;
+
+        virtual void update
+        (
+            const scalar dt,
+            scalar& d,
+            scalar& liquidCore,
+            scalar& tc,
+            const scalar rho,
+            const scalar mu,
+            const scalar sigma,
+            const scalar volHlowRate,
+            const scalar rhoAv,
+            const scalar Urel,
+            const vector& pos,
+            const vector& injectionPos,
+            const scalar pAmbient,
+            const scalar chi,
+            cachedRandom& rndGen
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "BlobsSheetAtomization.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.C

@@ -0,0 +1,279 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "LISAAtomization.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::LISAAtomization<CloudType>::LISAAtomization
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    AtomizationModel<CloudType>(dict, owner, typeName),
+    Cl_(readScalar(this->coeffDict().lookup("Cl"))),
+    cTau_(readScalar(this->coeffDict().lookup("cTau"))),
+    Q_(readScalar(this->coeffDict().lookup("Q"))),
+    lisaExp_(readScalar(this->coeffDict().lookup("lisaExp"))),
+    injectorDirection_(this->coeffDict().lookup("injectorDirection")),
+    SMDCalcMethod_(this->coeffDict().lookup("SMDCalculationMethod"))
+{
+    // Note: Would be good if this could be picked up from the injector
+    injectorDirection_ /= mag(injectorDirection_);
+
+    if (SMDCalcMethod_ == "method1")
+    {
+        SMDMethod_ = method1;
+    }
+    else if (SMDCalcMethod_ == "method2")
+    {
+        SMDMethod_ = method2;
+    }
+    else
+    {
+        SMDMethod_ = method2;
+        Info<< "Warning: SMDCalculationMethod " << SMDCalcMethod_
+            << " unknown. Options are (method1 | method2). Using method2"
+            << endl;
+    }
+}
+
+template <class CloudType>
+Foam::LISAAtomization<CloudType>::LISAAtomization
+(
+    const LISAAtomization<CloudType>& am
+)
+:
+    AtomizationModel<CloudType>(am),
+    Cl_(am.Cl_),
+    cTau_(am.cTau_),
+    Q_(am.Q_),
+    lisaExp_(am.lisaExp_),
+    injectorDirection_(am.injectorDirection_),
+    SMDCalcMethod_(am.SMDCalcMethod_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::LISAAtomization<CloudType>::~LISAAtomization()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::LISAAtomization<CloudType>::initLiquidCore() const
+{
+    return 1.0;
+}
+
+
+template<class CloudType>
+bool Foam::LISAAtomization<CloudType>::calcChi() const
+{
+    return true;
+}
+
+
+template<class CloudType>
+void Foam::LISAAtomization<CloudType>::update
+(
+    const scalar dt,
+    scalar& d,
+    scalar& liquidCore,
+    scalar& tc,
+    const scalar rho,
+    const scalar mu,
+    const scalar sigma,
+    const scalar volFlowRate,
+    const scalar rhoAv,
+    const scalar Urel,
+    const vector& pos,
+    const vector& injectionPos,
+    const scalar pAmbient,
+    const scalar chi,
+    cachedRandom& rndGen
+) const
+{
+    if (volFlowRate < SMALL)
+    {
+        return;
+    }
+
+    scalar tau = 0.0;
+    scalar dL = 0.0;
+    scalar k = 0.0;
+
+    // update atomization characteristic time
+    tc += dt;
+
+    scalar We = 0.5*rhoAv*sqr(Urel)*d/sigma;
+    scalar nu = mu/rho;
+
+    scalar Q = rhoAv/rho;
+
+    vector diff = pos - injectionPos;
+    scalar pWalk = mag(diff);
+    scalar traveledTime = pWalk/Urel;
+
+    scalar h = diff & injectorDirection_;
+    scalar delta = sqrt(sqr(pWalk) - sqr(h));
+
+    scalar hSheet = volFlowRate/(constant::mathematical::pi*delta*Urel);
+
+    // update drop diameter
+    d = min(d, hSheet);
+
+    if (We > 27.0/16.0)
+    {
+        scalar kPos = 0.0;
+        scalar kNeg = Q*sqr(Urel)*rho/sigma;
+
+        scalar derivPos = sqrt(Q*sqr(Urel));
+
+        scalar derivNeg =
+            (
+                8.0*sqr(nu)*pow3(kNeg)
+              + Q*sqr(Urel)*kNeg
+              - 3.0*sigma/2.0/rho*sqr(kNeg)
+            )
+          / sqrt
+            (
+                4.0*sqr(nu)*pow4(kNeg)
+              + Q*sqr(Urel)*sqr(kNeg)
+              - sigma*pow3(kNeg)/rho
+            )
+          - 4.0*nu*kNeg;
+
+        scalar kOld = 0.0;
+
+        for (label i=0; i<40; i++)
+        {
+            k = kPos - (derivPos/((derivNeg - derivPos)/(kNeg - kPos)));
+
+            scalar derivk =
+                (
+                    8.0*sqr(nu)*pow3(k)
+                  + Q*sqr(Urel)*k
+                  - 3.0*sigma/2.0/rho*sqr(k)
+                )
+              / sqrt
+                (
+                    4.0*sqr(nu)*pow4(k)
+                  + Q*sqr(Urel)*sqr(k)
+                  - sigma*pow3(k)/rho
+                )
+              - 4.0*nu*k;
+
+            if (derivk > 0)
+            {
+                derivPos = derivk;
+                kPos = k;
+            }
+            else
+            {
+                derivNeg = derivk;
+                kNeg = k;
+            }
+
+            if (mag(k - kOld)/k < 1e-4)
+            {
+                break;
+            }
+
+            kOld = k;
+        }
+
+        scalar omegaS =
+          - 2.0*nu*sqr(k)
+          + sqrt
+            (
+                4.0*sqr(nu)*pow4(k)
+              + Q*sqr(Urel)*sqr(k)
+              - sigma*pow3(k)/rho
+            );
+
+        tau = cTau_/omegaS;
+
+        dL = sqrt(8.0*d/k);
+    }
+    else
+    {
+        k = rhoAv*sqr(Urel)/2.0*sigma;
+
+        scalar J = 0.5*traveledTime*hSheet;
+
+        tau = pow(3.0*cTau_, 2.0/3.0)*cbrt(J*sigma/(sqr(Q)*pow4(Urel)*rho));
+
+        dL = sqrt(4.0*d/k);
+    }
+
+    scalar kL = 1.0/(dL*sqrt(0.5 + 1.5 * mu/sqrt(rho*sigma*dL)));
+
+    scalar dD = cbrt(3.0*constant::mathematical::pi*sqr(dL)/kL);
+
+    scalar atmPressure = 1.0e+5;
+
+    scalar pRatio = pAmbient/atmPressure;
+
+    dD = dD*pow(pRatio, lisaExp_);
+
+    scalar pExp = 0.135;
+
+    //  modifing dD to take account of flash boiling
+    dD = dD*(1.0 - chi*pow(pRatio, -pExp));
+    scalar lBU = Cl_ * mag(Urel)*tau;
+
+    if (pWalk > lBU)
+    {
+        scalar x = 0;
+
+        switch (SMDMethod_)
+        {
+            case method1:
+            {
+                #include "LISASMDCalcMethod1.H"
+                break;
+            }
+            case method2:
+            {
+                #include "LISASMDCalcMethod2.H"
+                break;
+            }
+        }
+
+        //  New droplet properties
+        liquidCore = 0.0;
+        d = x;
+        tc = 0.0;
+    }
+}
+
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.H

@@ -0,0 +1,161 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::LISAAtomization
+
+Description
+    Primary Breakup Model for pressure swirl atomizers.
+
+    Accurate description in
+    @verbatim
+    P.K. Senecal, D.P. Schmidt, I. Nouar, C.J. Rutland, R.D. Reitz, M. Corradini
+    "Modeling high-speed viscous liquid sheet atomization"
+    International Journal of Multiphase Flow 25 (1999) pags. 1073-1097
+    @endverbatim
+
+    and
+
+    @verbatim
+    D.P. Schmidt, I. Nouar, P.K. Senecal, C.J. Rutland, J.K. Martin, R.D. Reitz
+    "Pressure-Swirl Atomization in the Near Field"
+    SAE Techical Paper Series 1999-01-0496
+    @endverbatim
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LISAAtomization_H
+#define LISAAtomization_H
+
+#include "AtomizationModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                       Class LISAAtomization Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class LISAAtomization
+:
+    public AtomizationModel<CloudType>
+{
+
+public:
+
+    //- Enumeration for SMD calculations
+    enum SMDMethods
+    {
+        method1,
+        method2
+    };
+
+
+private:
+
+    // Private data
+
+        scalar Cl_;
+        scalar cTau_;
+        scalar Q_;
+        scalar lisaExp_;
+        vector injectorDirection_;
+        word SMDCalcMethod_;
+
+        SMDMethods SMDMethod_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("LISA");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        LISAAtomization(const dictionary&, CloudType&);
+
+        //- Construct copy
+        LISAAtomization(const LISAAtomization<CloudType>& am);
+
+        //- Construct and return a clone
+        virtual autoPtr<AtomizationModel<CloudType> > clone() const
+        {
+            return autoPtr<AtomizationModel<CloudType> >
+            (
+                new LISAAtomization<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~LISAAtomization();
+
+
+    // Member Functions
+
+        //- initial value of liquidCore
+        virtual scalar initLiquidCore() const;
+
+        //- flag to indicate if chi needs to be calculated
+        virtual bool calcChi() const;
+
+        virtual void update
+        (
+            const scalar dt,
+            scalar& d,
+            scalar& liquidCore,
+            scalar& tc,
+            const scalar rho,
+            const scalar mu,
+            const scalar sigma,
+            const scalar volFlowRate,
+            const scalar rhoAv,
+            const scalar Urel,
+            const vector& pos,
+            const vector& injectionPos,
+            const scalar pAmbient,
+            const scalar chi,
+            cachedRandom& rndGen
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "LISAAtomization.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISASMDCalcMethod1.H

@@ -0,0 +1,33 @@
+{
+    // calculate the new diameter with a Rosin Rammler distribution
+
+    scalar minValue = min(d, dD/10.0);
+    scalar maxValue = dD;
+
+    if(maxValue - minValue < SMALL)
+    {
+        minValue = d/10.0;
+    }
+
+    scalar range = maxValue - minValue;
+
+    scalar y = 0;
+
+    bool success = false;
+
+    while(!success)
+    {
+
+        x = minValue + range*rndGen.sample01<scalar>();
+        y = rndGen.sample01<scalar>();
+        scalar p = 0.0;
+        scalar nExp = 1;
+        scalar xx = pow(x/dD, nExp);
+
+        p = xx*exp(-xx);
+        if (y<p)
+        {
+            success = true;
+        }
+    }
+}

src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISASMDCalcMethod2.H

@@ -0,0 +1,53 @@
+{
+    // calculate the new diameter with the standard 1D Rosin Rammler
+    // distribution. Calculation of the mean radius based on SMR rs.
+    // Coefficient factorGamma depends on nExp. Note that Reitz either used
+    // (Schmidt et al., 1999-01-0496) or skipped (Senecal et al.) this factor!
+    // scalar factorGamma = 0.75*sqrt(mathematicalConstant::pi);   //nExp=2
+    scalar factorGamma = 1.;
+    scalar delta = dD/factorGamma;
+
+    // dD is the SMD, and the delta is calculated using gamma
+    // function. Here we assume nExp = 2
+    scalar minValue = dD/10.0;
+    scalar maxValue = dD;
+
+    // The pdf value for 4.0*delta is already very small.
+    // scalar maxValue = delta*4.0;
+
+    if (maxValue - minValue < SMALL)
+    {
+        minValue = maxValue/20.0;
+    }
+
+    scalar range = maxValue - minValue;
+
+    scalar nExp = 3;
+    FixedList<scalar, 500> rrd;
+    scalar probFactorMin = exp(-pow(minValue/delta, nExp));
+    scalar probFactorMax = exp(-pow(maxValue/delta, nExp));
+    scalar probFactor = 1.0/(probFactorMin - probFactorMax);
+
+    forAll(rrd, n)
+    {
+        scalar xx = minValue + range*n/500;
+        rrd[n] = (probFactorMin - exp(-pow(xx/delta, nExp)))*probFactor;
+    }
+
+
+    bool success = false;
+
+    scalar y = rndGen.sample01<scalar>();
+    label k = 0;
+
+    while(!success && (k<500))
+    {
+        if (rrd[k] > y)
+        {
+            success = true;
+        }
+        k++;
+    }
+
+    x = minValue + range*(k - 0.5)/500.0;
+}

src/lagrangian/spray/submodels/AtomizationModel/NoAtomization/NoAtomization.C

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "NoAtomization.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::NoAtomization<CloudType>::NoAtomization
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    AtomizationModel<CloudType>(owner)
+{}
+
+
+template <class CloudType>
+Foam::NoAtomization<CloudType>::NoAtomization
+(
+    const NoAtomization<CloudType>& am
+)
+:
+    AtomizationModel<CloudType>(am)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::NoAtomization<CloudType>::~NoAtomization()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+bool Foam::NoAtomization<CloudType>::active() const
+{
+    return false;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::NoAtomization<CloudType>::initLiquidCore() const
+{
+    return 0.0;
+}
+
+
+template<class CloudType>
+bool Foam::NoAtomization<CloudType>::calcChi() const
+{
+    return false;
+}
+
+
+template<class CloudType>
+void Foam::NoAtomization<CloudType>::update
+(
+    const scalar dt,
+    scalar& d,
+    scalar& liquidCore,
+    scalar& tc,
+    const scalar rho,
+    const scalar mu,
+    const scalar sigma,
+    const scalar volFlowRate,
+    const scalar rhoAv,
+    const scalar Urel,
+    const vector& pos,
+    const vector& injectionPos,
+    const scalar pAmbient,
+    const scalar chi,
+    cachedRandom& rndGen
+) const
+{}
+
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/AtomizationModel/NoAtomization/NoAtomization.H

@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::NoAtomization
+
+Description
+    Dummy phase change model for 'none'
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef NoAtomization_H
+#define NoAtomization_H
+
+#include "AtomizationModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                       Class NoAtomization Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class NoAtomization
+:
+    public AtomizationModel<CloudType>
+{
+public:
+
+    //- Runtime type information
+    TypeName("none");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        NoAtomization(const dictionary& dict, CloudType& cloud);
+
+        //- Construct copy
+        NoAtomization(const NoAtomization<CloudType>& am);
+
+        //- Construct and return a clone
+        virtual autoPtr<AtomizationModel<CloudType> > clone() const
+        {
+            return autoPtr<AtomizationModel<CloudType> >
+            (
+                new NoAtomization<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~NoAtomization();
+
+
+    // Member Functions
+
+        //- Flag to indicate whether model activates atomization model
+        virtual bool active() const;
+
+        //- Initial value of liquidCore
+        virtual scalar initLiquidCore() const;
+
+        //- Flag to indicate if chi needs to be calculated
+        virtual bool calcChi() const;
+
+        virtual void update
+        (
+            const scalar dt,
+            scalar& d,
+            scalar& liquidCore,
+            scalar& tc,
+            const scalar rho,
+            const scalar mu,
+            const scalar sigma,
+            const scalar volFlowRate,
+            const scalar rhoAv,
+            const scalar Urel,
+            const vector& pos,
+            const vector& injectionPos,
+            const scalar pAmbient,
+            const scalar chi,
+            cachedRandom& rndGen
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "NoAtomization.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C

@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "BreakupModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::BreakupModel<CloudType>::BreakupModel
+(
+    CloudType& owner
+)
+:
+    SubModelBase<CloudType>(owner),
+    solveOscillationEq_(false),
+    y0_(0.0),
+    yDot0_(0.0),
+    TABComega_(0.0),
+    TABCmu_(0.0),
+    TABWeCrit_(0.0)
+{}
+
+
+template<class CloudType>
+Foam::BreakupModel<CloudType>::BreakupModel
+(
+    const BreakupModel<CloudType>& bum
+)
+:
+    SubModelBase<CloudType>(bum),
+    solveOscillationEq_(bum.solveOscillationEq_),
+    y0_(bum.y0_),
+    yDot0_(bum.yDot0_),
+    TABComega_(bum.TABComega_),
+    TABCmu_(bum.TABCmu_),
+    TABWeCrit_(bum.TABWeCrit_)
+{}
+
+
+template<class CloudType>
+Foam::BreakupModel<CloudType>::BreakupModel
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& type
+)
+:
+    SubModelBase<CloudType>(owner, dict, type),
+    solveOscillationEq_(this->coeffDict().lookup("solveOscillationEq")),
+    y0_(0.0),
+    yDot0_(0.0),
+    TABComega_(0.0),
+    TABCmu_(0.0),
+    TABWeCrit_(0.0)
+{
+    if (solveOscillationEq_)
+    {
+        const dictionary TABcoeffsDict(dict.subDict("TABCoeffs"));
+        y0_ = TABcoeffsDict.template lookupOrDefault<scalar>("y0", 0.0);
+        yDot0_ = TABcoeffsDict.template lookupOrDefault<scalar>("yDot0", 0.0);
+        TABComega_ =
+            TABcoeffsDict.template lookupOrDefault<scalar>("Comega", 8.0);
+        TABCmu_ = TABcoeffsDict.template lookupOrDefault<scalar>("Cmu", 10.0);
+        TABWeCrit_ =
+            TABcoeffsDict.template lookupOrDefault<scalar>("WeCrit", 12.0);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::BreakupModel<CloudType>::~BreakupModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+bool Foam::BreakupModel<CloudType>::update
+(
+    const scalar dt,
+    const vector& g,
+    scalar& d,
+    scalar& tc,
+    scalar& ms,
+    scalar& nParticle,
+    scalar& KHindex,
+    scalar& y,
+    scalar& yDot,
+    const scalar d0,
+    const scalar rho,
+    const scalar mu,
+    const scalar sigma,
+    const vector& U,
+    const scalar rhoc,
+    const scalar muc,
+    const vector& Urel,
+    const scalar Urmag,
+    const scalar tMom,
+    const scalar averageParcelMass,
+    scalar& dChild,
+    scalar& massChild,
+    cachedRandom& rndGen
+) const
+{
+    notImplemented
+    (
+        "bool Foam::BreakupModel<CloudType>::update"
+        "("
+            "const scalar, "
+            "const vector&, "
+            "scalar&, "
+            "scalar&, "
+            "scalar&, "
+            "scalar&, "
+            "scalar&, "
+            "scalar&, "
+            "scalar&, "
+            "const scalar, "
+            "const scalar, "
+            "const scalar, "
+            "const scalar, "
+            "const vector&, "
+            "const scalar, "
+            "const scalar, "
+            "const vector&, "
+            "const scalar, "
+            "const scalar, "
+            "const scalar, "
+            "scalar&, "
+            "scalar&, "
+            "cachedRandom&"
+        ") const;"
+    );
+
+    return false;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "BreakupModelNew.C"
+
+// ************************************************************************* //
+

src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H

@@ -0,0 +1,235 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::BreakupModel
+
+Description
+    Templated break-up model class
+
+SourceFiles
+    BreakupModel.C
+    BreakupModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BreakupModel_H
+#define BreakupModel_H
+
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class BreakupModel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class BreakupModel
+:
+    public SubModelBase<CloudType>
+{
+
+protected:
+
+    // Protected data
+
+        Switch solveOscillationEq_;
+
+        scalar y0_;
+        scalar yDot0_;
+        scalar TABComega_;
+        scalar TABCmu_;
+        scalar TABWeCrit_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("breakupModel");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        BreakupModel,
+        dictionary,
+        (
+            const dictionary& dict,
+            CloudType& owner
+        ),
+        (dict, owner)
+    );
+
+
+    // Constructors
+
+        //- Construct null from owner
+        BreakupModel(CloudType& owner);
+
+        //- Construct from dictionary
+        BreakupModel
+        (
+            const dictionary& dict,
+            CloudType& owner,
+            const word& type
+        );
+
+        //- Construct copy
+        BreakupModel(const BreakupModel<CloudType>& bum);
+
+        //- Construct and return a clone
+        virtual autoPtr<BreakupModel<CloudType> > clone() const
+        {
+            return autoPtr<BreakupModel<CloudType> >
+            (
+                new BreakupModel<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~BreakupModel();
+
+
+    //- Selector
+    static autoPtr<BreakupModel<CloudType> > New
+    (
+        const dictionary& dict,
+        CloudType& owner
+    );
+
+
+    // Access
+
+        inline const Switch& solveOscillationEq() const
+        {
+            return solveOscillationEq_;
+        }
+
+        inline const scalar& y0() const
+        {
+            return y0_;
+        }
+
+        inline const scalar& yDot0() const
+        {
+            return yDot0_;
+        }
+
+        inline const scalar& TABComega() const
+        {
+            return TABComega_;
+        }
+
+        inline const scalar& TABCmu() const
+        {
+            return TABCmu_;
+        }
+
+        inline const scalar& TABWeCrit() const
+        {
+            return TABWeCrit_;
+        }
+
+
+    // Member Functions
+
+        //- Update the parcel properties and return true if a child parcel
+        //  should be added
+        virtual bool update
+        (
+            const scalar dt,
+            const vector& g,
+            scalar& d,
+            scalar& tc,
+            scalar& ms,
+            scalar& nParticle,
+            scalar& KHindex,
+            scalar& y,
+            scalar& yDot,
+            const scalar d0,
+            const scalar rho,
+            const scalar mu,
+            const scalar sigma,
+            const vector& U,
+            const scalar rhoc,
+            const scalar muc,
+            const vector& Urel,
+            const scalar Urmag,
+            const scalar tMom,
+            const scalar averageParcelMass,
+            scalar& dChild,
+            scalar& massChild,
+            cachedRandom& rndGen
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeBreakupModel(CloudType)                                           \
+                                                                              \
+    typedef CloudType::sprayCloudType sprayCloudType;                         \
+    defineNamedTemplateTypeNameAndDebug                                       \
+    (                                                                         \
+        BreakupModel<sprayCloudType>,                                         \
+        0                                                                     \
+    );                                                                        \
+    defineTemplateRunTimeSelectionTable                                       \
+    (                                                                         \
+        BreakupModel<sprayCloudType>,                                         \
+        dictionary                                                            \
+    );
+
+
+#define makeBreakupModelType(SS, CloudType)                                   \
+                                                                              \
+    typedef CloudType::sprayCloudType sprayCloudType;                         \
+    defineNamedTemplateTypeNameAndDebug(SS<sprayCloudType>, 0);               \
+                                                                              \
+    BreakupModel<sprayCloudType>::                                            \
+        adddictionaryConstructorToTable<SS<sprayCloudType> >                  \
+            add##SS##CloudType##sprayCloudType##ConstructorToTable_;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "BreakupModel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModelNew.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "BreakupModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::autoPtr<Foam::BreakupModel<CloudType> >
+Foam::BreakupModel<CloudType>::New
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+{
+    word BreakupModelType(dict.lookup("breakupModel"));
+
+    Info<< "Selecting BreakupModel " << BreakupModelType << endl;
+
+    typename dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(BreakupModelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "BreakupModel<CloudType>::New"
+            "("
+                "const dictionary&, "
+                "CloudType&"
+            ")"
+        )   << "Unknown BreakupModelType type "
+            << BreakupModelType
+            << ", constructor not in hash table" << nl << nl
+            << "    Valid BreakupModel types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc() << exit(FatalError);
+    }
+
+    return autoPtr<BreakupModel<CloudType> >(cstrIter()(dict, owner));
+}
+
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C

@@ -0,0 +1,228 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "ETAB.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::ETAB<CloudType>::ETAB
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    BreakupModel<CloudType>(dict, owner, typeName),
+    Cmu_(this->coeffDict().template lookupOrDefault<scalar>("Cmu", 10.0)),
+    Comega_(this->coeffDict().template lookupOrDefault<scalar>("Comega", 8.0)),
+    k1_(this->coeffDict().template lookupOrDefault<scalar>("k1", 0.2)),
+    k2_(this->coeffDict().template lookupOrDefault<scalar>("k2", 0.2)),
+    WeCrit_
+(
+    this->coeffDict().template lookupOrDefault<scalar>("WeCrit", 12.0)
+),
+    WeTransition_
+    (
+        this->coeffDict().template lookupOrDefault<scalar>
+        (
+            "WeTransition",
+            100.0
+        )
+    ),
+    AWe_(0.0)
+{
+    scalar k21 = k2_/k1_;
+    AWe_ = (k21*sqrt(WeTransition_) - 1.0)/pow4(WeTransition_);
+
+    if (!BreakupModel<CloudType>::solveOscillationEq_)
+    {
+        Info<< "Warning: solveOscillationEq is set to "
+            << BreakupModel<CloudType>::solveOscillationEq_ << nl
+            << " Setting it to true in order for the ETAB model to work."
+            << endl;
+        BreakupModel<CloudType>::solveOscillationEq_ = true;
+    }
+}
+
+
+template<class CloudType>
+Foam::ETAB<CloudType>::ETAB(const ETAB<CloudType>& bum)
+:
+    BreakupModel<CloudType>(bum),
+    Cmu_(bum.Cmu_),
+    Comega_(bum.Comega_),
+    k1_(bum.k1_),
+    k2_(bum.k2_),
+    WeCrit_(bum.WeCrit_),
+    WeTransition_(bum.WeTransition_),
+    AWe_(bum.AWe_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::ETAB<CloudType>::~ETAB()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+bool Foam::ETAB<CloudType>::update
+(
+    const scalar dt,
+    const vector& g,
+    scalar& d,
+    scalar& tc,
+    scalar& ms,
+    scalar& nParticle,
+    scalar& KHindex,
+    scalar& y,
+    scalar& yDot,
+    const scalar d0,
+    const scalar rho,
+    const scalar mu,
+    const scalar sigma,
+    const vector& U,
+    const scalar rhoc,
+    const scalar muc,
+    const vector& Urel,
+    const scalar Urmag,
+    const scalar tMom,
+    const scalar averageParcelMass,
+    scalar& dChild,
+    scalar& massChild,
+    cachedRandom& rndGen
+) const
+{
+    scalar r  = 0.5*d;
+    scalar r2 = r*r;
+    scalar r3 = r*r2;
+
+    scalar semiMass = nParticle*pow3(d);
+
+    // inverse of characteristic viscous damping time
+    scalar rtd = 0.5*Cmu_*mu/(rho*r2);
+
+    // oscillation frequency (squared)
+    scalar omega2 = Comega_*sigma/(rho*r3) - rtd*rtd;
+
+    if (omega2 > 0)
+    {
+        scalar omega = sqrt(omega2);
+        scalar romega = 1.0/omega;
+
+        scalar We = rhoc*sqr(Urmag)*r/sigma;
+        scalar Wetmp = We/WeCrit_;
+
+        scalar y1 = y - Wetmp;
+        scalar y2 = yDot*romega;
+
+        scalar a = sqrt(y1*y1 + y2*y2);
+
+        // scotty we may have break-up
+        if (a + Wetmp > 1.0)
+        {
+            scalar phic = y1/a;
+
+            // constrain phic within -1 to 1
+            phic = max(min(phic, 1), -1);
+
+            scalar phit = acos(phic);
+            scalar phi = phit;
+            scalar quad = -y2/a;
+            if (quad < 0)
+            {
+                phi = 2*constant::mathematical::pi - phit;
+            }
+
+            scalar tb = 0;
+
+            if (mag(y) < 1.0)
+            {
+                scalar theta = acos((1.0 - Wetmp)/a);
+
+                if (theta < phi)
+                {
+                    if (2*constant::mathematical::pi - theta >= phi)
+                    {
+                        theta = -theta;
+                    }
+                    theta += 2*constant::mathematical::pi;
+                }
+                tb = (theta - phi)*romega;
+
+                // breakup occurs
+                if (dt > tb)
+                {
+                    y = 1.0;
+                    yDot = -a*omega*sin(omega*tb + phi);
+                }
+            }
+
+            // update droplet size
+            if (dt > tb)
+            {
+                scalar sqrtWe = AWe_*pow4(We) + 1.0;
+                scalar Kbr = k1_*omega*sqrtWe;
+
+                if (We > WeTransition_)
+                {
+                    sqrtWe = sqrt(We);
+                    Kbr =k2_*omega*sqrtWe;
+                }
+
+                scalar rWetmp = 1.0/Wetmp;
+                scalar cosdtbu = max(-1.0, min(1.0, 1.0 - rWetmp));
+                scalar dtbu = romega*acos(cosdtbu);
+                scalar decay = exp(-Kbr*dtbu);
+
+                scalar rNew = decay*r;
+                if (rNew < r)
+                {
+                    d = 2.0*rNew;
+                    y = 0.0;
+                    yDot = 0.0;
+                }
+            }
+        }
+    }
+    else
+    {
+        // reset droplet distortion parameters
+        y = 0;
+        yDot = 0;
+    }
+
+    // update the nParticle count to conserve mass
+    nParticle = semiMass/pow3(d);
+
+    // Do not add child parcel
+    return false;
+}
+
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.H

@@ -0,0 +1,159 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::ETAB
+
+Description
+    The Enhanced TAB model.
+
+    Described in the papers below.
+    @verbatim
+    F.X. Tanner
+        "Liquid Jet Atomization and Droplet Breakup Modeling of
+        Non-Evaporating Diesel Fuel Sprays"
+        SAE 970050,
+        SAE Transactions: Journal of Engines, Vol 106, Sec 3 pp 127-140
+
+    F.X. Tanner and G. Weisser
+        "Simulation of Liquid Jet Atomization for
+        Fuel Sprays by Means of Cascade Drop Breakup Model"
+        SAE 980808
+        SAE Technical Paper Series
+    @endverbatim
+
+See Also
+    The TAB model
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ETAB_H
+#define ETAB_H
+
+#include "BreakupModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                       Class ETAB Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class ETAB
+:
+    public BreakupModel<CloudType>
+{
+private:
+
+    // Private data
+
+        // Model constants
+
+            // Cmu_ and Comega_ are the same as in the TAB model
+            scalar Cmu_;
+            scalar Comega_;
+
+            scalar k1_;
+            scalar k2_;
+            scalar WeCrit_;
+            scalar WeTransition_;
+            scalar AWe_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("ETAB");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        ETAB(const dictionary&, CloudType&);
+
+        //- Construct copy
+        ETAB(const ETAB<CloudType>& bum);
+
+        //- Construct and return a clone
+        virtual autoPtr<BreakupModel<CloudType> > clone() const
+        {
+            return autoPtr<BreakupModel<CloudType> >
+            (
+                new ETAB<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~ETAB();
+
+
+    // Member Functions
+
+        //- update the parcel properties
+        virtual bool update
+        (
+            const scalar dt,
+            const vector& g,
+            scalar& d,
+            scalar& tc,
+            scalar& ms,
+            scalar& nParticle,
+            scalar& KHindex,
+            scalar& y,
+            scalar& yDot,
+            const scalar d0,
+            const scalar rho,
+            const scalar mu,
+            const scalar sigma,
+            const vector& U,
+            const scalar rhoc,
+            const scalar muc,
+            const vector& Urel,
+            const scalar Urmag,
+            const scalar tMom,
+            const scalar averageParcelMass,
+            scalar& dChild,
+            scalar& massChild,
+            cachedRandom& rndGen
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "ETAB.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/BreakupModel/NoBreakup/NoBreakup.C

@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "NoBreakup.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::NoBreakup<CloudType>::NoBreakup
+(
+    const dictionary&,
+    CloudType& owner
+)
+:
+    BreakupModel<CloudType>(owner)
+{}
+
+
+
+template <class CloudType>
+Foam::NoBreakup<CloudType>::NoBreakup(const NoBreakup<CloudType>& bum)
+:
+    BreakupModel<CloudType>(bum)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::NoBreakup<CloudType>::~NoBreakup()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+bool Foam::NoBreakup<CloudType>::active() const
+{
+    return false;
+}
+
+
+template<class CloudType>
+bool Foam::NoBreakup<CloudType>::update
+(
+    const scalar dt,
+    const vector& g,
+    scalar& d,
+    scalar& tc,
+    scalar& ms,
+    scalar& nParticle,
+    scalar& KHindex,
+    scalar& y,
+    scalar& yDot,
+    const scalar d0,
+    const scalar rho,
+    const scalar mu,
+    const scalar sigma,
+    const vector& U,
+    const scalar rhoc,
+    const scalar muc,
+    const vector& Urel,
+    const scalar Urmag,
+    const scalar tMom,
+    const scalar averageParcelMass,
+    scalar& dChild,
+    scalar& massChild,
+    cachedRandom& rndGen
+) const
+{
+    // Do nothing
+    return false;
+}
+
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/BreakupModel/NoBreakup/NoBreakup.H

@@ -0,0 +1,127 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::NoBreakup
+
+Description
+    Dummy breakup model for 'none'
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef NoBreakup_H
+#define NoBreakup_H
+
+#include "BreakupModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                       Class NoBreakup Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class NoBreakup
+:
+    public BreakupModel<CloudType>
+{
+public:
+
+    //- Runtime type information
+    TypeName("none");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        NoBreakup(const dictionary&, CloudType&);
+
+        //- Construct copy
+        NoBreakup(const NoBreakup<CloudType>& bum);
+
+        //- Construct and return a clone
+        virtual autoPtr<BreakupModel<CloudType> > clone() const
+        {
+            return autoPtr<BreakupModel<CloudType> >
+            (
+                new NoBreakup<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~NoBreakup();
+
+
+    // Member Functions
+
+        //- Flag to indicate whether model activates break-up model
+        virtual bool active() const;
+
+        //- update the parcel properties
+        virtual bool update
+        (
+            const scalar dt,
+            const vector& g,
+            scalar& d,
+            scalar& tc,
+            scalar& ms,
+            scalar& nParticle,
+            scalar& KHindex,
+            scalar& y,
+            scalar& yDot,
+            const scalar d0,
+            const scalar rho,
+            const scalar mu,
+            const scalar sigma,
+            const vector& U,
+            const scalar rhoc,
+            const scalar muc,
+            const vector& Urel,
+            const scalar Urmag,
+            const scalar tMom,
+            const scalar averageParcelMass,
+            scalar& dChild,
+            scalar& massChild,
+            cachedRandom& rndGen
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "NoBreakup.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.C

@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PilchErdman.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::PilchErdman<CloudType>::PilchErdman
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    BreakupModel<CloudType>(dict, owner, typeName),
+    B1_(this->coeffDict().template lookupOrDefault<scalar>("B1", 0.375)),
+    B2_(this->coeffDict().template lookupOrDefault<scalar>("B2", 0.236))
+{}
+
+
+template <class CloudType>
+Foam::PilchErdman<CloudType>::PilchErdman(const PilchErdman<CloudType>& bum)
+:
+    BreakupModel<CloudType>(bum),
+    B1_(bum.B1_),
+    B2_(bum.B2_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::PilchErdman<CloudType>::~PilchErdman()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+bool Foam::PilchErdman<CloudType>::update
+(
+    const scalar dt,
+    const vector& g,
+    scalar& d,
+    scalar& tc,
+    scalar& ms,
+    scalar& nParticle,
+    scalar& KHindex,
+    scalar& y,
+    scalar& yDot,
+    const scalar d0,
+    const scalar rho,
+    const scalar mu,
+    const scalar sigma,
+    const vector& U,
+    const scalar rhoc,
+    const scalar muc,
+    const vector& Urel,
+    const scalar Urmag,
+    const scalar tMom,
+    const scalar averageParcelMass,
+    scalar& dChild,
+    scalar& massChild,
+    cachedRandom& rndGen
+) const
+{
+    scalar semiMass = nParticle*pow3(d);
+    scalar We = 0.5*rhoc*sqr(Urmag)*d/sigma;
+    scalar Oh = mu/pow(rho*d*sigma, 0.5);
+
+    scalar Wec = 6.0*(1.0 + 1.077*pow(Oh, 1.6));
+
+    if (We > Wec)
+    {
+        // We > 1335, wave crest stripping
+        scalar taubBar = 5.5;
+
+        if (We > 175.0)
+        {
+            // sheet stripping
+            taubBar = 0.766*pow(2.0*We - 12.0, 0.25);
+        }
+        else if (We > 22.0)
+        {
+            // Bag-and-stamen breakup
+            taubBar = 14.1*pow(2.0*We - 12.0, -0.25);
+        }
+        else if (We > 9.0)
+        {
+            // Bag breakup
+            taubBar = 2.45*pow(2.0*We - 12.0, 0.25);
+        }
+        else if (We > 6.0)
+        {
+            // Vibrational breakup
+            taubBar = 6.0*pow(2.0*We - 12.0, -0.25);
+        }
+
+        scalar rho12 = pow(rhoc/rho, 0.5);
+
+        scalar Vd = Urmag*rho12*(B1_*taubBar * B2_*taubBar*taubBar);
+        scalar Vd1 = sqr(1.0 - Vd/Urmag);
+        Vd1 = max(Vd1, SMALL);
+        scalar Ds = 2.0*Wec*sigma*Vd1/(Vd1*rhoc*sqr(Urmag));
+        scalar A = Urmag*rho12/d;
+
+        scalar taub = taubBar/A;
+
+        scalar frac = dt/taub;
+
+        // update the droplet diameter according to the rate eq. (implicitly)
+        d = (d + frac*Ds)/(1.0 + frac);
+
+        // correct the number of particles to conserve mass
+        nParticle = semiMass/pow3(d);
+    }
+
+    return false;
+}
+
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.H

@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::PilchErdman
+
+Description
+    secondary breakup model
+
+    @verbatim
+    Pilch, M. and Erdman, C.A.
+    "Use of breakup time data and velocity history data
+     to predict the maximum size of stable fragments for acceleration
+     induced breakup of a liquid drop."
+    Int. J. Multiphase Flows 13 (1987), 741-757
+    @endverbatim
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PilchErdman_H
+#define PilchErdman_H
+
+#include "BreakupModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                       Class PilchErdman Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class PilchErdman
+:
+    public BreakupModel<CloudType>
+{
+private:
+
+    // Private data
+
+        scalar B1_;
+        scalar B2_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("PilchErdman");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        PilchErdman(const dictionary&, CloudType&);
+
+        //- Construct copy
+        PilchErdman(const PilchErdman<CloudType>& bum);
+
+        //- Construct and return a clone
+        virtual autoPtr<BreakupModel<CloudType> > clone() const
+        {
+            return autoPtr<BreakupModel<CloudType> >
+            (
+                new PilchErdman<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~PilchErdman();
+
+
+    // Member Functions
+
+        //- update the parcel properties
+        virtual bool update
+        (
+            const scalar dt,
+            const vector& g,
+            scalar& d,
+            scalar& tc,
+            scalar& ms,
+            scalar& nParticle,
+            scalar& KHindex,
+            scalar& y,
+            scalar& yDot,
+            const scalar d0,
+            const scalar rho,
+            const scalar mu,
+            const scalar sigma,
+            const vector& U,
+            const scalar rhoc,
+            const scalar muc,
+            const vector& Urel,
+            const scalar Urmag,
+            const scalar tMom,
+            const scalar averageParcelMass,
+            scalar& dChild,
+            scalar& massChild,
+            cachedRandom& rndGen
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "PilchErdman.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/BreakupModel/ReitzDiwakar/ReitzDiwakar.C

@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "ReitzDiwakar.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::ReitzDiwakar<CloudType>::ReitzDiwakar
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    BreakupModel<CloudType>(dict, owner, typeName),
+    Cbag_(this->coeffDict().template lookupOrDefault<scalar>("Cbag", 6.0)),
+    Cb_(this->coeffDict().template lookupOrDefault<scalar>("Cb", 0.785)),
+    Cstrip_(this->coeffDict().template lookupOrDefault<scalar>("Cstrip", 0.5)),
+    Cs_(this->coeffDict().template lookupOrDefault<scalar>("Cs", 10.0))
+{}
+
+
+template <class CloudType>
+Foam::ReitzDiwakar<CloudType>::ReitzDiwakar(const ReitzDiwakar<CloudType>& bum)
+:
+    BreakupModel<CloudType>(bum),
+    Cbag_(bum.Cbag_),
+    Cb_(bum.Cb_),
+    Cstrip_(bum.Cstrip_),
+    Cs_(bum.Cs_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::ReitzDiwakar<CloudType>::~ReitzDiwakar()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+bool Foam::ReitzDiwakar<CloudType>::update
+(
+    const scalar dt,
+    const vector& g,
+    scalar& d,
+    scalar& tc,
+    scalar& ms,
+    scalar& nParticle,
+    scalar& KHindex,
+    scalar& y,
+    scalar& yDot,
+    const scalar d0,
+    const scalar rho,
+    const scalar mu,
+    const scalar sigma,
+    const vector& U,
+    const scalar rhoc,
+    const scalar muc,
+    const vector& Urel,
+    const scalar Urmag,
+    const scalar tMom,
+    const scalar averageParcelMass,
+    scalar& dChild,
+    scalar& massChild,
+    cachedRandom& rndGen
+) const
+{
+    scalar d1 = d;
+    scalar nuc = muc/rhoc;
+    scalar We = 0.5*rhoc*sqr(Urmag)*d/sigma;
+    scalar Re = Urmag*d/nuc;
+
+    scalar sqRey = sqrt(Re);
+
+    if (We > Cbag_)
+    {
+        if (We > Cstrip_*sqRey)
+        {
+            scalar dStrip = sqr(2.0*Cstrip_*sigma)/(rhoc*pow3(Urmag)*muc);
+            scalar tauStrip = Cs_*d*sqrt(rho/rhoc)/Urmag;
+            scalar fraction = dt/tauStrip;
+
+            // new droplet diameter, implicit calculation
+            d = (fraction*dStrip + d)/(1.0 + fraction);
+        }
+        else
+        {
+            scalar dBag = 2.0*Cbag_*sigma/(rhoc*sqr(Urmag));
+
+            scalar tauBag = Cb_*d*sqrt(rho*d/sigma);
+
+            scalar fraction = dt/tauBag;
+
+            // new droplet diameter, implicit calculation
+            d = (fraction*dBag + d)/(1.0 + fraction);
+        }
+
+        // preserve the total mass/volume by updating the number of
+        // particles in parcels due to breakup
+        nParticle *= pow(d1/d, 3.0);
+    }
+
+    return false;
+}
+
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/BreakupModel/ReitzDiwakar/ReitzDiwakar.H

@@ -0,0 +1,152 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::ReitzDiwakar
+
+Description
+    secondary breakup model
+
+    @verbatim
+    Reitz, R.D.
+    "Modelling atomization processes in highpressure vaporizing sprays"
+    Atomization and Spray Technology 3 (1987), 309-337
+    @endverbatim
+
+    @verbatim
+    Reitz, R.D. and Diwakar, R.
+    "Effect of drop breakup on fuel sprays"
+    SAE Tech. paper series, 860469 (1986)
+    @endverbatim
+
+    @verbatim
+    Reitz, R.D. and Diwakar, R.
+    "Structure of high-pressure fuel sprays"
+    SAE Tech. paper series, 870598 (1987)
+    @endverbatim
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ReitzDiwakar_H
+#define ReitzDiwakar_H
+
+#include "BreakupModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                       Class ReitzDiwakar Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class ReitzDiwakar
+:
+    public BreakupModel<CloudType>
+{
+private:
+
+    // Private data
+
+        scalar Cbag_;
+        scalar Cb_;
+        scalar Cstrip_;
+        scalar Cs_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("ReitzDiwakar");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        ReitzDiwakar(const dictionary&, CloudType&);
+
+        //- Construct copy
+        ReitzDiwakar(const ReitzDiwakar<CloudType>& bum);
+
+        //- Construct and return a clone
+        virtual autoPtr<BreakupModel<CloudType> > clone() const
+        {
+            return autoPtr<BreakupModel<CloudType> >
+            (
+                new ReitzDiwakar<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~ReitzDiwakar();
+
+
+    // Member Functions
+
+        //- update the parcel properties
+        virtual bool update
+        (
+            const scalar dt,
+            const vector& g,
+            scalar& d,
+            scalar& tc,
+            scalar& ms,
+            scalar& nParticle,
+            scalar& KHindex,
+            scalar& y,
+            scalar& yDot,
+            const scalar d0,
+            const scalar rho,
+            const scalar mu,
+            const scalar sigma,
+            const vector& U,
+            const scalar rhoc,
+            const scalar muc,
+            const vector& Urel,
+            const scalar Urmag,
+            const scalar tMom,
+            const scalar averageParcelMass,
+            scalar& dChild,
+            scalar& massChild,
+            cachedRandom& rndGen
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "ReitzDiwakar.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //