updated lagrangian and radiation libraries

src/lagrangian/Allwmake

@@ -3,6 +3,5 @@ set -x
 
 wmake libso basic
 wmake libso solidParticle
-wmake libso basic
 wmake libso intermediate
 wmake libso dieselSpray

src/lagrangian/basic/Cloud/Cloud.C

@@ -43,7 +43,6 @@ Foam::Cloud<ParticleType>::Cloud
     cloud(pMesh),
     IDLList<ParticleType>(particles),
     polyMesh_(pMesh),
-    cloudName_("defaultCloud"),
     allFaces_(pMesh.faces()),
     points_(pMesh.points()),
     cellFaces_(pMesh.cells()),
@@ -65,7 +64,6 @@ Foam::Cloud<ParticleType>::Cloud
     cloud(pMesh, cloudName),
     IDLList<ParticleType>(particles),
     polyMesh_(pMesh),
-    cloudName_(cloudName),
     allFaces_(pMesh.faces()),
     points_(pMesh.points()),
     cellFaces_(pMesh.cells()),

src/lagrangian/basic/Cloud/Cloud.H

@@ -26,7 +26,6 @@ Class
     Foam::Cloud
 
 Description
-    Foam::Cloud
 
 SourceFiles
     Cloud.C
@@ -74,7 +73,6 @@ class Cloud
     // Private data
 
         const polyMesh& polyMesh_;
-        const word cloudName_;
         const faceList& allFaces_;
         const vectorField& points_;
         const cellList& cellFaces_;
@@ -152,12 +150,6 @@ public:
                 return polyMesh_;
             }
 
-            //- Return the cloud name reference
-            const word& cloudName() const
-            {
-                return cloudName_;
-            }
-
             //- Is this global face an internal face?
             bool internalFace(const label facei) const
             {
@@ -213,6 +205,11 @@ public:
                 return IDLList<ParticleType>::end();
             };
 
+            void clear()
+            {
+                return IDLList<ParticleType>::clear();
+            };
+
 
         // Edit
 

src/lagrangian/basic/Cloud/CloudIO.C

@@ -60,7 +60,6 @@ Foam::Cloud<ParticleType>::Cloud
 :
     cloud(pMesh),
     polyMesh_(pMesh),
-    cloudName_("defaultCloud"),
     allFaces_(pMesh.faces()),
     points_(pMesh.points()),
     cellFaces_(pMesh.cells()),
@@ -83,7 +82,6 @@ Foam::Cloud<ParticleType>::Cloud
 :
     cloud(pMesh, cloudName),
     polyMesh_(pMesh),
-    cloudName_(cloudName),
     allFaces_(pMesh.faces()),
     points_(pMesh.points()),
     cellFaces_(pMesh.cells()),

src/lagrangian/basic/IOPosition/IOPosition.C

@@ -22,6 +22,8 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
+Description
+
 \*---------------------------------------------------------------------------*/
 
 #include "IOPosition.H"

src/lagrangian/basic/Particle/Particle.C

@@ -39,7 +39,7 @@ template<class ParticleType>
 Foam::labelList Foam::Particle<ParticleType>::findFaces
 (
     const vector& position
-)
+) const
 {
     const polyMesh& mesh = cloud_.polyMesh_;
     const labelList& faces = mesh.cells()[celli_];
@@ -69,7 +69,7 @@ Foam::labelList Foam::Particle<ParticleType>::findFaces
     const vector& position,
     const label celli,
     const scalar stepFraction
-)
+) const
 {
     const polyMesh& mesh = cloud_.polyMesh_;
     const labelList& faces = mesh.cells()[celli];
@@ -94,11 +94,11 @@ Foam::labelList Foam::Particle<ParticleType>::findFaces
 
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::Particle<ParticleType>::prepareForParallelTransfer
 (
     const label patchi,
-    TrackingData& td
+    TrackData& td
 )
 {
     // Convert the face index to be local to the processor patch
@@ -107,11 +107,11 @@ void Foam::Particle<ParticleType>::prepareForParallelTransfer
 
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::Particle<ParticleType>::correctAfterParallelTransfer
 (
     const label patchi,
-    TrackingData& td
+    TrackData& td
 )
 {
     const processorPolyPatch& ppp =
@@ -156,7 +156,15 @@ void Foam::Particle<ParticleType>::correctAfterParallelTransfer
     }
 
     // Reset the face index for the next tracking operation
-    facei_ = -1;
+    if (stepFraction_ > (1.0 - SMALL))
+    {
+        stepFraction_ = 1.0;
+        facei_ = -1;
+    }
+    else
+    {
+        facei_ += ppp.start();
+    }
 }
 
 
@@ -181,11 +189,11 @@ Foam::Particle<ParticleType>::Particle
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 Foam::label Foam::Particle<ParticleType>::track
 (
     const vector& endPosition,
-    TrackingData& td
+    TrackData& td
 )
 {
     facei_ = -1;
@@ -200,6 +208,7 @@ Foam::label Foam::Particle<ParticleType>::track
 }
 
 
+
 template<class ParticleType>
 Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
 {
@@ -207,13 +216,12 @@ Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
     return track(endPosition, dummyTd);
 }
 
-
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 Foam::scalar Foam::Particle<ParticleType>::trackToFace
 (
     const vector& endPosition,
-    TrackingData& td
+    TrackData& td
 )
 {
     const polyMesh& mesh = cloud_.polyMesh_;
@@ -301,8 +309,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
             {
                 FatalErrorIn
                 (
-                    "Particle::trackToFace"
+                    "Particle::trackToFace(const vector&, TrackData&)"
-                    "(const vector&, TrackingData&)"
                 )<< "addressing failure" << nl
                  << abort(FatalError);
             }
@@ -389,7 +396,6 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     return trackFraction;
 }
 
-
 template<class ParticleType>
 Foam::scalar Foam::Particle<ParticleType>::trackToFace
 (
@@ -400,7 +406,6 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     return trackToFace(endPosition, dummyTd);
 }
 
-
 template<class ParticleType>
 void Foam::Particle<ParticleType>::transformPosition(const tensor& T)
 {
@@ -419,11 +424,11 @@ void Foam::Particle<ParticleType>::transformProperties(const vector&)
 
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::Particle<ParticleType>::hitWedgePatch
 (
     const wedgePolyPatch& wpp,
-    TrackingData&
+    TrackData&
 )
 {
     vector nf = wpp.faceAreas()[wpp.whichFace(facei_)];
@@ -434,11 +439,11 @@ void Foam::Particle<ParticleType>::hitWedgePatch
 
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::Particle<ParticleType>::hitSymmetryPatch
 (
     const symmetryPolyPatch& spp,
-    TrackingData&
+    TrackData&
 )
 {
     vector nf = spp.faceAreas()[spp.whichFace(facei_)];
@@ -449,11 +454,11 @@ void Foam::Particle<ParticleType>::hitSymmetryPatch
 
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::Particle<ParticleType>::hitCyclicPatch
 (
     const cyclicPolyPatch& cpp,
-    TrackingData&
+    TrackData&
 )
 {
     label patchFacei_ = cpp.whichFace(facei_);
@@ -481,31 +486,31 @@ void Foam::Particle<ParticleType>::hitCyclicPatch
 
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::Particle<ParticleType>::hitProcessorPatch
 (
     const processorPolyPatch& spp,
-    TrackingData& td
+    TrackData& td
 )
 {}
 
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::Particle<ParticleType>::hitWallPatch
 (
     const wallPolyPatch& spp,
-    TrackingData&
+    TrackData&
 )
 {}
 
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::Particle<ParticleType>::hitPatch
 (
     const polyPatch& spp,
-    TrackingData&
+    TrackData&
 )
 {}
 

src/lagrangian/basic/Particle/Particle.H

@@ -26,7 +26,6 @@ Class
     Foam::Particle
 
 Description
-    Foam::Particle
 
 \*---------------------------------------------------------------------------*/
 
@@ -79,6 +78,41 @@ class Particle
     public IDLList<ParticleType>::link
 {
 
+public:
+
+    //- Class used to pass tracking data to the trackToFace function
+    class trackData
+    {
+
+        // Private data
+
+            //- Reference to the cloud containing this particle
+            Cloud<ParticleType>& cloud_;
+
+
+    public:
+
+        bool switchProcessor;
+        bool keepParticle;
+
+
+        // Constructors
+
+            inline trackData
+            (
+                Cloud<ParticleType>& cloud
+            );
+
+
+        // Member functions
+
+            //- Return a reference to the cloud
+            inline Cloud<ParticleType>& cloud();
+    };
+
+
+private:
+
     // Private data
 
         //- Reference to the particle cloud
@@ -108,7 +142,7 @@ class Particle
             const vector& to,
             const label facei,
             const scalar stepFraction
-        );
+        ) const;
 
         //- Return the 'lambda' value for the position, p, on the face,
         // where, p = from + lamda*(to - from)
@@ -118,13 +152,13 @@ class Particle
             const vector& from,
             const vector& to,
             const label facei
-        );
+        ) const;
 
         //- Return the faces between position and cell centre
         labelList findFaces
         (
             const vector& position
-        );
+        ) const;
 
         //- Return the faces between position and cell centre
         labelList findFaces
@@ -132,80 +166,91 @@ class Particle
             const vector& position,
             const label celli,
             const scalar stepFraction
-        );
+        ) const;
 
 
 protected:
 
+    // Patch interactions
+
         //- Overridable function to handle the particle hitting a wedgePatch
-        template<class TrackingData>
+        template<class TrackData>
         void hitWedgePatch
         (
             const wedgePolyPatch&,
-            TrackingData& td
+            TrackData& td
         );
 
-        //- Overridable function to handle the particle hitting a symmetryPatch
+        //- Overridable function to handle the particle hitting a
-        template<class TrackingData>
+        //  symmetryPatch
+        template<class TrackData>
         void hitSymmetryPatch
         (
             const symmetryPolyPatch&,
-            TrackingData& td
+            TrackData& td
         );
 
         //- Overridable function to handle the particle hitting a cyclicPatch
-        template<class TrackingData>
+        template<class TrackData>
         void hitCyclicPatch
         (
             const cyclicPolyPatch&,
-            TrackingData& td
+            TrackData& td
         );
 
-        //- Overridable function to handle the particle hitting a processorPatch
+        //- Overridable function to handle the particle hitting a
-        template<class TrackingData>
+        //  processorPatch
+        template<class TrackData>
         void hitProcessorPatch
         (
             const processorPolyPatch&,
-            TrackingData& td
+            TrackData& td
         );
 
         //- Overridable function to handle the particle hitting a wallPatch
-        template<class TrackingData>
+        template<class TrackData>
         void hitWallPatch
         (
             const wallPolyPatch&,
-            TrackingData& td
+            TrackData& td
         );
 
-        //- Overridable function to handle the particle hitting a general patch
+        //- Overridable function to handle the particle hitting a
-        template<class TrackingData>
+        //  general patch
+        template<class TrackData>
         void hitPatch
         (
             const polyPatch&,
-            TrackingData& td
+            TrackData& td
         );
 
+
+    // Transformations
+
         //- Transform the position the particle
         //  according to the given transformation tensor
-        void transformPosition(const tensor& T);
+        virtual void transformPosition(const tensor& T);
 
         //- Transform the physical properties of the particle
         //  according to the given transformation tensor
-        void transformProperties(const tensor& T);
+        virtual void transformProperties(const tensor& T);
 
         //- Transform the physical properties of the particle
         //  according to the given separation vector
-        void transformProperties(const vector& separation);
+        virtual void transformProperties(const vector& separation);
+
+
+    // Parallel transfer
 
         //- Convert global addressing to the processor patch
         //  local equivalents
-        template<class TrackingData>
+        template<class TrackData>
-        void prepareForParallelTransfer(const label patchi, TrackingData& td);
+        void prepareForParallelTransfer(const label patchi, TrackData& td);
 
         //- Convert processor patch addressing to the global equivalents
         //  and set the celli to the face-neighbour
-        template<class TrackingData>
+        template<class TrackData>
-        void correctAfterParallelTransfer(const label patchi, TrackingData& td);
+        void correctAfterParallelTransfer(const label patchi, TrackData& td);
 
 
 public:
@@ -217,36 +262,6 @@ public:
     TypeName("Particle");
 
 
-    //- Class used to pass tracking data to the trackToFace function
-    class trackData
-    {
-
-        // Private data
-
-            //- Reference to the cloud containing this particle
-            Cloud<ParticleType>& cloud_;
-
-
-    public:
-
-        bool switchProcessor;
-        bool keepParticle;
-
-
-        // Constructors
-
-            inline trackData
-            (
-                Cloud<ParticleType>& cloud
-            );
-
-
-        // Member functions
-
-            inline Cloud<ParticleType>& cloud();
-    };
-
-
     // Constructors
 
         //- Construct from components
@@ -265,7 +280,8 @@ public:
             bool readFields = true
         );
 
-        //- Factory class to read-construct particles used for parallel transfer
+        //- Factory class to read-construct particles used for
+        //  parallel transfer
         class iNew
         {
 
@@ -288,7 +304,7 @@ public:
         };
 
 
-    // Destructor
+    //- Destructor
 
         virtual ~Particle()
         {}
@@ -299,7 +315,7 @@ public:
         // Access
 
             //- Return true if particle is in cell
-            inline bool inCell();
+            inline bool inCell() const;
 
             //- Return true if position is in cell i
             inline bool inCell
@@ -307,7 +323,7 @@ public:
                 const vector& position,
                 const label celli,
                 const scalar stepFraction
-            );
+            ) const;
 
             //- Return current particle position
             inline const vector& position() const;
@@ -337,7 +353,11 @@ public:
             inline label patch(const label facei) const;
 
             //- Which face of this patch is this particle on
-            inline label patchFace(const label patchi, const label facei) const;
+            inline label patchFace
+            (
+                const label patchi,
+                const label facei
+            ) const;
 
             //- The nearest distance to a wall that
             //  the particle can be in the n direction
@@ -359,14 +379,14 @@ public:
             //  the fraction of the time-step completed.
             //  Returns the boundary face index if the track stops at the
             //  boundary, -1 otherwise.
-            template<class TrackingData>
+            template<class TrackData>
             label track
             (
                 const vector& endPosition,
-                TrackingData& td
+                TrackData& td
             );
 
-            //- Calls the templated track with dummy TrackingData
+            //- Calls the templated track with dummy TrackData
             label track(const vector& endPosition);
 
             //- Track particle to a given position and returns 1.0 if the
@@ -375,17 +395,18 @@ public:
             //  completed.
             //  on entry 'stepFraction()' should be set to the fraction of the
             //  time-step at which the tracking starts.
-            template<class TrackingData>
+            template<class TrackData>
             scalar trackToFace
             (
                 const vector& endPosition,
-                TrackingData& td
+                TrackData& td
             );
 
-            //- Calls the templated trackToFace with dummy TrackingData
+            //- Calls the templated trackToFace with dummy TrackData
             scalar trackToFace(const vector& endPosition);
 
-            //- Return the index of the face to be used in the interpolation routine
+            //- Return the index of the face to be used in the interpolation
+            //  routine
             inline label faceInterpolation() const;
 
 

src/lagrangian/basic/Particle/ParticleI.H

@@ -40,7 +40,7 @@ inline scalar Particle<ParticleType>::lambda
     const vector& to,
     const label facei,
     const scalar stepFraction
-)
+) const
 {
     const polyMesh& mesh = cloud_.polyMesh_;
     bool movingMesh = mesh.moving();
@@ -71,7 +71,6 @@ inline scalar Particle<ParticleType>::lambda
 
         // for a moving mesh we need to reconstruct the old
         // Sf and Cf from oldPoints (they aren't stored)
-        // NN.
 
         const vectorField& oldPoints = mesh.oldPoints();
 
@@ -100,7 +99,7 @@ inline scalar Particle<ParticleType>::lambda
             // find center of rotation
             vector omega = Sf0 ^ Sf;
             scalar magOmega = mag(omega);
-            omega /= magOmega+SMALL;
+            omega /= magOmega + SMALL;
             vector n0 = omega ^ Sf0;
             scalar lam = ((Cf - Cf0) & Sf)/(n0 & Sf);
             vector r0 = Cf0 + lam*n0;
@@ -188,7 +187,7 @@ inline scalar Particle<ParticleType>::lambda
     const vector& from,
     const vector& to,
     const label facei
-)
+) const
 {
     const polyMesh& mesh = cloud_.polyMesh_;
 
@@ -235,7 +234,7 @@ inline scalar Particle<ParticleType>::lambda
 
 
 template<class ParticleType>
-inline bool Particle<ParticleType>::inCell()
+inline bool Particle<ParticleType>::inCell() const
 {
     labelList faces = findFaces(position_);
 
@@ -249,7 +248,7 @@ inline bool Particle<ParticleType>::inCell
     const vector& position,
     const label celli,
     const scalar stepFraction
-)
+) const
 {
     labelList faces = findFaces(position, celli, stepFraction);
 
@@ -277,7 +276,6 @@ inline Cloud<ParticleType>& Particle<ParticleType>::trackData::cloud()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-//- Return reference to the particle cloud
 template<class ParticleType>
 inline const Cloud<ParticleType>& Particle<ParticleType>::cloud() const
 {
@@ -313,7 +311,6 @@ inline label Particle<ParticleType>::face() const
 }
 
 
-//- Is the particle on a boundary face?
 template<class ParticleType>
 inline bool Particle<ParticleType>::onBoundary() const
 {

src/lagrangian/basic/indexedParticle/indexedParticle.H

@@ -26,7 +26,6 @@ Class
     Foam::indexedParticle
 
 Description
-    An indexed Particle
 
 SourceFiles
     indexedParticleI.H

src/lagrangian/basic/indexedParticle/indexedParticleCloud.C

@@ -22,6 +22,8 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
+Description
+
 \*---------------------------------------------------------------------------*/
 
 #include "indexedParticle.H"

src/lagrangian/basic/lagrangianField/LagrangianField.H

@@ -26,7 +26,6 @@ Class
     Foam::LagrangianField
 
 Description
-    Lagrangian Field
 
 SourceFiles
     LagrangianFieldI.H

src/lagrangian/basic/passiveParticle/passiveParticle.H

@@ -26,7 +26,6 @@ Class
     Foam::passiveParticle
 
 Description
-    A passive Particle
 
 SourceFiles
     passiveParticleI.H

src/lagrangian/basic/passiveParticle/passiveParticleCloud.C

@@ -22,6 +22,8 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
+Description
+
 \*---------------------------------------------------------------------------*/
 
 #include "passiveParticle.H"

src/lagrangian/basic/polyMeshInfo/polyMeshInfo.H

@@ -26,7 +26,7 @@ Class
     Foam::polyMeshInfo
 
 Description
-    Foam::polyMeshInfo
+
 
 SourceFiles
     polyMeshInfo.C

src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "Analytical.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Analytical<Type>::Analytical
+(
+    const word& phiName,
+    const dictionary& dict
+)
+:
+    IntegrationScheme<Type>(phiName, dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Analytical<Type>::~Analytical()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Type Foam::Analytical<Type>::integrate
+(
+    const Type phi,
+    const scalar dt,
+    const Type alpha,
+    const scalar beta
+) const
+{
+    return alpha + (phi - alpha)*exp(-beta*dt);
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H

@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::Analytical
+
+Description
+    Analytical integration
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Analytical_H
+#define Analytical_H
+
+#include "IntegrationScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class Analytical Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class Analytical
+:
+    public IntegrationScheme<Type>
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("Analytical");
+
+
+    // Constructors
+
+        //- Construct from components
+        Analytical
+        (
+            const word& phiName,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+
+        virtual ~Analytical();
+
+
+    // Member Functions
+
+        //- Perform the integration
+        virtual Type integrate
+        (
+            const Type phi,
+            const scalar dt,
+            const Type alpha,
+            const scalar beta
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "Analytical.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "Euler.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Euler<Type>::Euler
+(
+    const word& phiName,
+    const dictionary& dict
+)
+:
+    IntegrationScheme<Type>(phiName, dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Euler<Type>::~Euler()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Type Foam::Euler<Type>::integrate
+(
+    const Type phi,
+    const scalar dt,
+    const Type alpha,
+    const scalar beta
+) const
+{
+    return (phi + dt*alpha)/(1.0 + dt/beta);
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H

@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::Euler
+
+Description
+    Euler-implicit integration
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Euler_H
+#define Euler_H
+
+#include "IntegrationScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Euler Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class Euler
+:
+    public IntegrationScheme<Type>
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("Euler");
+
+
+    // Constructors
+
+        //- Construct from components
+        Euler
+        (
+            const word& phiName,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+
+        virtual ~Euler();
+
+
+    // Member Functions
+
+        //- Perform the integration
+        virtual Type integrate
+        (
+            const Type phi,
+            const scalar dt,
+            const Type alpha,
+            const scalar beta
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "Euler.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.C

@@ -0,0 +1,54 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "IntegrationScheme.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::IntegrationScheme<Type>::IntegrationScheme
+(
+    const word& phiName,
+    const dictionary& dict
+)
+:
+   phiName_(phiName),
+   dict_(dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor    * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::IntegrationScheme<Type>::~IntegrationScheme()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+#include "newIntegrationScheme.C"
+
+// ************************************************************************* //

src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H

@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+
+Class
+    Foam::IntegrationScheme
+
+Description
+    Top level model for Integration schemes
+
+SourceFiles
+    IntegrationScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IntegrationScheme_H
+#define IntegrationScheme_H
+
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "dictionary.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class IntegrationScheme Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class IntegrationScheme
+{
+    // Private data
+
+        //- Name of the Integration variable
+        const word& phiName_;
+
+        //- Reference to the dictionary
+        const dictionary& dict_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        IntegrationScheme(const IntegrationScheme&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const IntegrationScheme&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("IntegrationScheme");
+
+
+    //- Declare runtime constructor selection table
+
+         declareRunTimeSelectionTable
+         (
+             autoPtr,
+             IntegrationScheme,
+             dictionary,
+             (
+                 const word& phiName,
+                 const dictionary& dict
+             ),
+             (phiName, dict)
+         );
+
+
+    // Constructors
+
+        //- Construct from components
+        IntegrationScheme
+        (
+            const word& phiName,
+            const dictionary& dict
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected radiation model
+        static autoPtr<IntegrationScheme> New
+        (
+            const word& phiName,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+
+        virtual ~IntegrationScheme();
+
+
+    // Member Functions
+
+        //- Perform the Integration
+        virtual Type integrate
+        (
+            const Type phi,
+            const scalar dt,
+            const Type alpha,
+            const scalar beta
+        ) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeIntegrationScheme(Type)                                           \
+                                                                              \
+    defineNamedTemplateTypeNameAndDebug(IntegrationScheme<Type>, 0);          \
+                                                                              \
+    defineTemplateRunTimeSelectionTable                                       \
+    (                                                                         \
+        IntegrationScheme<Type>,                                               \
+        dictionary                                                            \
+    );
+
+
+#define makeIntegrationSchemeType(SS, Type)                                   \
+                                                                              \
+    defineNamedTemplateTypeNameAndDebug(SS<Type>, 0);                         \
+                                                                              \
+    IntegrationScheme<Type>::adddictionaryConstructorToTable<SS<Type> >          \
+        add##SS##Type##ConstructorToTable_;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "IntegrationScheme.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+

src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationSchemesFwd.H

@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+InClass
+    Foam::IntegrationScheme
+
+Description
+
+SourceFiles
+    IntegrationScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IntegrationSchemesFwd_H
+#define IntegrationSchemesFwd_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+class IntegrationScheme;
+
+typedef IntegrationScheme<scalar> scalarIntegrationScheme;
+typedef IntegrationScheme<vector> vectorIntegrationScheme;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/newIntegrationScheme.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "IntegrationScheme.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::autoPtr<Foam::IntegrationScheme<Type> >
+Foam::IntegrationScheme<Type>::New
+(
+    const word& phiName,
+    const dictionary& dict
+)
+{
+    word IntegrationSchemeTypeName;
+
+    dict.lookup(phiName) >> IntegrationSchemeTypeName;
+
+    Info<< "Selecting " << phiName << " IntegrationScheme "
+        << IntegrationSchemeTypeName << endl;
+
+    typename dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(IntegrationSchemeTypeName);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "IntegrationScheme::New(const dictionary&)"
+        )   << "Unknown IntegrationScheme type "
+            << IntegrationSchemeTypeName << nl << nl
+            << "Valid IntegrationScheme types are:" << nl
+            << dictionaryConstructorTablePtr_->toc() << nl
+            << exit(FatalError);
+    }
+
+    return autoPtr<IntegrationScheme<Type> >(cstrIter()(phiName, dict));
+}
+
+// ************************************************************************* //

src/lagrangian/intermediate/IntegrationScheme/makeIntegrationSchemes.C

@@ -0,0 +1,46 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "IntegrationScheme.H"
+#include "Euler.H"
+#include "Analytical.H"
+
+#include "scalar.H"
+#include "vector.H"
+
+namespace Foam
+{
+    makeIntegrationScheme(scalar);
+    makeIntegrationSchemeType(Euler, scalar);
+    makeIntegrationSchemeType(Analytical, scalar);
+
+    makeIntegrationScheme(vector);
+    makeIntegrationSchemeType(Euler, vector);
+    makeIntegrationSchemeType(Analytical, vector);
+};
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/Make/files

@@ -44,5 +44,8 @@ parcels/derived/reactingParcel/makeReactingParcelWallInteractionModels.C
 submodels/addOns/radiation/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C
 submodels/addOns/radiation/scatter/cloudScatter/cloudScatter.C
 
+/* integration schemes */
+IntegrationScheme/makeIntegrationSchemes.C
+
 
 LIB = $(FOAM_LIBBIN)/liblagrangianIntermediate

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

@@ -30,6 +30,8 @@ License
 #include "InjectionModel.H"
 #include "WallInteractionModel.H"
 
+#include "IntegrationScheme.H"
+
 #include "interpolationCellPoint.H"
 
 // * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
@@ -153,7 +155,6 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
     kinematicCloud(),
     cloudType_(cloudType),
     mesh_(rho.mesh()),
-    runTime_(rho.time()),
     vpi_(vpi),
     particleProperties_
     (
@@ -170,7 +171,7 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
     parcelTypeId_(readLabel(particleProperties_.lookup("parcelTypeId"))),
     coupled_(particleProperties_.lookup("coupled")),
     rndGen_(label(0)),
-    time0_(runTime_.value()),
+    time0_(this->db().time().value()),
     parcelBasisType_(particleProperties_.lookup("parcelBasisType")),
     parcelBasis_(pbNumber),
     massTotal_
@@ -182,6 +183,7 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
     U_(U),
     mu_(mu),
     g_(g),
+    interpolationSchemes_(particleProperties_.subDict("interpolationSchemes")),
     dispersionModel_
     (
         DispersionModel<KinematicCloud<ParcelType> >::New
@@ -214,6 +216,14 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
             *this
         )
     ),
+    UIntegrator_
+    (
+        vectorIntegrationScheme::New
+        (
+            "U",
+            particleProperties_.subDict("integrationSchemes")
+        )
+    ),
     nInjections_(0),
     nParcelsAdded_(0),
     nParcelsAddedTotal_(0),
@@ -221,8 +231,8 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
     (
         IOobject
         (
-            this->cloudName() + "UTrans",
+            this->name() + "UTrans",
-            runTime_.timeName(),
+            this->db().time().timeName(),
             this->db(),
             IOobject::NO_READ,
             IOobject::NO_WRITE,
@@ -235,8 +245,8 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
     (
         IOobject
         (
-            this->cloudName() + "UCoeff",
+            this->name() + "UCoeff",
-            runTime_.timeName(),
+            this->db().time().timeName(),
             this->db(),
             IOobject::NO_READ,
             IOobject::NO_WRITE,
@@ -288,17 +298,34 @@ void Foam::KinematicCloud<ParcelType>::resetSourceTerms()
 template<class ParcelType>
 void Foam::KinematicCloud<ParcelType>::evolve()
 {
-    interpolationCellPoint<scalar> rhoInterp(vpi_, rho_);
+    autoPtr<interpolation<scalar> > rhoInterpolator = interpolation<scalar>::New
-    interpolationCellPoint<vector> UInterp(vpi_, U_);
+    (
-    interpolationCellPoint<scalar> muInterp(vpi_, mu_);
+        interpolationSchemes_,
+        vpi_,
+        rho_
+    );
+
+    autoPtr<interpolation<vector> > UInterpolator = interpolation<vector>::New
+    (
+        interpolationSchemes_,
+        vpi_,
+        U_
+    );
+
+    autoPtr<interpolation<scalar> > muInterpolator = interpolation<scalar>::New
+    (
+        interpolationSchemes_,
+        vpi_,
+        mu_
+    );
 
     typename ParcelType::trackData td
     (
         *this,
         constProps_,
-        rhoInterp,
+        rhoInterpolator(),
-        UInterp,
+        UInterpolator(),
-        muInterp,
+        muInterpolator(),
         g_.value()
     );
 
@@ -330,7 +357,7 @@ void Foam::KinematicCloud<ParcelType>::inject
     TrackingData& td
 )
 {
-    scalar time = runTime_.value();
+    scalar time = this->db().time().value();
 
     scalar pRho = td.constProps().rho0();
 
@@ -363,7 +390,7 @@ void Foam::KinematicCloud<ParcelType>::inject
     // Duration of injection period during this timestep
     scalar deltaT = min
     (
-        runTime().deltaT().value(),
+        this->db().time().deltaT().value(),
         min
         (
             time - this->injection().timeStart(),
@@ -430,8 +457,8 @@ void Foam::KinematicCloud<ParcelType>::inject
 
             scalar dt = time - timeInj;
 
-            pPtr->stepFraction() = (runTime_.deltaT().value() - dt)
+            pPtr->stepFraction() = (this->db().time().deltaT().value() - dt)
-                /runTime_.deltaT().value();
+                /this->time().deltaT().value();
 
             this->injectParcel(td, pPtr);
          }
@@ -466,7 +493,7 @@ void Foam::KinematicCloud<ParcelType>::postInjectCheck()
 {
     if (nParcelsAdded_)
     {
-        Pout<< "\n--> Cloud: " << this->cloudName() << nl <<
+        Pout<< "\n--> Cloud: " << this->name() << nl <<
                "    Added " << nParcelsAdded_ <<  " new parcels" << nl << endl;
     }
 
@@ -474,7 +501,7 @@ void Foam::KinematicCloud<ParcelType>::postInjectCheck()
     nParcelsAdded_ = 0;
 
     // Set time for start of next injection
-    time0_ = runTime_.value();
+    time0_ = this->db().time().value();
 
     // Increment number of injections
     nInjections_++;
@@ -484,7 +511,7 @@ void Foam::KinematicCloud<ParcelType>::postInjectCheck()
 template<class ParcelType>
 void Foam::KinematicCloud<ParcelType>::info() const
 {
-    Info<< "Cloud name: " << this->cloudName() << nl
+    Info<< "Cloud name: " << this->name() << nl
         << "    Parcels added during this run   = "
         << returnReduce(nParcelsAddedTotal_, sumOp<label>()) << nl
         << "    Mass introduced during this run = "
@@ -502,8 +529,8 @@ void Foam::KinematicCloud<ParcelType>::dumpParticlePositions() const
 {
     OFstream pObj
     (
-        this->runTime().path()/"parcelPositions_"
+        this->db().time().path()/"parcelPositions_"
-      + this->cloudName() + "_"
+      + this->name() + "_"
       + name(this->nInjections_) + ".obj"
     );
 

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

@@ -27,7 +27,6 @@ Class
 
 Description
     Templated base class for kinematic cloud
-
     - Kinematic only
       - Dispersion model
       - Drag model
@@ -53,7 +52,8 @@ SourceFiles
 #include "volPointInterpolation.H"
 #include "fvMatrices.H"
 #include "fvm.H"
-//#include "DimensionedField.H"
+
+#include "IntegrationSchemesFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -108,7 +108,7 @@ private:
 
         //- References to the mesh and time databases
         const fvMesh& mesh_;
-        const Time& runTime_;
+//        const Time& runTime_;
 
         //- Reference to the interpolation for the carrier phase to the parcels
         const volPointInterpolation& vpi_;
@@ -165,6 +165,11 @@ private:
             const dimensionedVector& g_;
 
 
+        // Interpolation
+
+            dictionary interpolationSchemes_;
+
+
         // References to the cloud sub-models
 
             //- Dispersion model
@@ -183,6 +188,12 @@ private:
                 wallInteractionModel_;
 
 
+        // Reference to the particle integration schemes
+
+            //- Velocity integration
+            autoPtr<vectorIntegrationScheme> UIntegrator_;
+
+
         // Counters
 
             //- Number of injections counter
@@ -273,9 +284,6 @@ public:
                 //- Return refernce to the mesh
                 inline const fvMesh& mesh() const;
 
-                //- Return referece to the time database
-                inline const Time& runTime() const;
-
                 //- Retuen reference to the interpolation
                 inline const volPointInterpolation& vpi() const;
 
@@ -313,6 +321,12 @@ public:
                 inline const dimensionedVector& g() const;
 
 
+            // Interpolations
+
+                //- Return reference to the interpolation dictionary
+                inline const dictionary& interpolationSchemes() const;
+
+
             // Sub-models
 
                 //- Return reference to dispersion model
@@ -335,6 +349,12 @@ public:
                     wallInteraction() const;
 
 
+            // Integration schemes
+
+                //-Return reference to velocity integration
+                inline const vectorIntegrationScheme& UIntegrator() const;
+
+
             // Sources
 
                 // Momentum
@@ -345,8 +365,11 @@ public:
                     //- Coefficient for carrier phase U equation
                     inline DimensionedField<scalar, volMesh>& UCoeff();
 
-                    //- Momentum source term
+                    //- Return tmp momentum source term - fully explicit
-                    inline tmp<fvVectorMatrix> SU(volVectorField& U) const;
+                    inline tmp<DimensionedField<vector, volMesh> > SU1() const;
+
+                    //- Return tmp momentum source term - semi-implicit
+                    inline tmp<fvVectorMatrix> SU2(volVectorField& U) const;
 
 
 

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

@@ -47,13 +47,6 @@ inline const Foam::fvMesh& Foam::KinematicCloud<ParcelType>::mesh() const
 }
 
 
-template<class ParcelType>
-inline const Foam::Time& Foam::KinematicCloud<ParcelType>::runTime() const
-{
-    return runTime_;
-}
-
-
 template<class ParcelType>
 inline const Foam::volPointInterpolation&
 Foam::KinematicCloud<ParcelType>::vpi() const
@@ -106,6 +99,13 @@ Foam::KinematicCloud<ParcelType>::g() const
     return g_;
 }
 
+template<class ParcelType>
+inline const Foam::dictionary&
+Foam::KinematicCloud<ParcelType>::interpolationSchemes() const
+{
+    return interpolationSchemes_;
+}
+
 
 template<class ParcelType>
 inline const Foam::DispersionModel<Foam::KinematicCloud<ParcelType> >&
@@ -147,6 +147,14 @@ Foam::KinematicCloud<ParcelType>::wallInteraction() const
 }
 
 
+template<class ParcelType>
+inline const Foam::vectorIntegrationScheme&
+Foam::KinematicCloud<ParcelType>::UIntegrator() const
+{
+    return UIntegrator_;
+}
+
+
 template<class ParcelType>
 inline const Foam::label Foam::KinematicCloud<ParcelType>::nInjections() const
 {
@@ -220,15 +228,51 @@ Foam::KinematicCloud<ParcelType>::UCoeff()
 
 
 template<class ParcelType>
-inline Foam::tmp<Foam::fvVectorMatrix>
+inline Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh> >
-Foam::KinematicCloud<ParcelType>::SU(volVectorField& U) const
+Foam::KinematicCloud<ParcelType>::SU1() const
 {
-    Info<< "UTrans min/max = "
+    tmp<DimensionedField<vector, volMesh> > tSU1
-        << min(UTrans_) << ", " << max(UTrans_) << endl;
+    (
-    Info<< "UCoeff min/max = "
+        new DimensionedField<vector, volMesh>
-        << min(UCoeff_) << ", " << max(UCoeff_) << endl;
+        (
+            IOobject
+            (
+                this->name() + "SU1",
+                this->db().time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            this->mesh(),
+            dimensionedVector
+            (
+                 "zero",
+                 dimensionSet(1, -2, -2, 0, 0),
+                 vector::zero
+            )
+        )
+    );
 
-    return UTrans_/(mesh_.V()*runTime_.deltaT())
+    vectorField& SU1 = tSU1().field();
+    SU1 = UTrans_/(mesh_.V()*this->db().time().deltaT());
+
+    return tSU1;
+}
+
+
+template<class ParcelType>
+inline Foam::tmp<Foam::fvVectorMatrix>
+Foam::KinematicCloud<ParcelType>::SU2(volVectorField& U) const
+{
+    if (debug)
+    {
+        Info<< "UTrans min/max = "
+            << min(UTrans_) << ", " << max(UTrans_) << endl;
+        Info<< "UCoeff min/max = "
+            << min(UCoeff_) << ", " << max(UCoeff_) << endl;
+    }
+
+    return UTrans_/(mesh_.V()*this->db().time().deltaT())
          - fvm::Sp(UCoeff_/mesh_.V(), U)
          + UCoeff_/mesh_.V()*U;
 }
@@ -244,8 +288,8 @@ Foam::KinematicCloud<ParcelType>::theta() const
         (
             IOobject
             (
-                this->cloudName() + "Theta",
+                this->name() + "Theta",
-                runTime_.timeName(),
+                this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
@@ -281,8 +325,8 @@ Foam::KinematicCloud<ParcelType>::alpha() const
         (
             IOobject
             (
-                this->cloudName() + "Alpha",
+                this->name() + "Alpha",
-                runTime_.timeName(),
+                this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,

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

@@ -86,15 +86,15 @@ Foam::ReactingCloud<ParcelType>::ReactingCloud
             (
                 IOobject
                 (
-                     this->cloudName() + "rhoTrans" + name(i),
+                     this->name() + "rhoTrans" + name(i),
-                     this->runTime().timeName(),
+                     this->db().time().timeName(),
                      this->db(),
                      IOobject::NO_READ,
                      IOobject::NO_WRITE,
                      false
                 ),
                 this->mesh(),
-                dimensionedScalar("zero", dimensionSet(1, 0, 0, 0, 0), 0.0)
+                dimensionedScalar("zero", dimMass, 0.0)
             )
         );
     }
@@ -128,23 +128,58 @@ void Foam::ReactingCloud<ParcelType>::evolve()
     const volScalarField cp = carrierThermo_.Cp();
     const volScalarField& p = carrierThermo_.p();
 
-    interpolationCellPoint<scalar> rhoInterp(this->vpi(), this->rho());
+    autoPtr<interpolation<scalar> > rhoInterpolator = interpolation<scalar>::New
-    interpolationCellPoint<vector> UInterp(this->vpi(), this->U());
+    (
-    interpolationCellPoint<scalar> muInterp(this->vpi(), this->mu());
+        this->interpolationSchemes(),
-    interpolationCellPoint<scalar> TInterp(this->vpi(), T);
+        this->vpi(),
-    interpolationCellPoint<scalar> cpInterp(this->vpi(), cp);
+        this->rho()
-    interpolationCellPoint<scalar> pInterp(this->vpi(), p);
+    );
+
+    autoPtr<interpolation<vector> > UInterpolator = interpolation<vector>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        this->U()
+    );
+
+    autoPtr<interpolation<scalar> > muInterpolator = interpolation<scalar>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        this->mu()
+    );
+
+    autoPtr<interpolation<scalar> > TInterpolator = interpolation<scalar>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        T
+    );
+
+    autoPtr<interpolation<scalar> > cpInterpolator = interpolation<scalar>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        cp
+    );
+
+    autoPtr<interpolation<scalar> > pInterpolator = interpolation<scalar>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        p
+    );
 
     typename ParcelType::trackData td
     (
         *this,
         constProps_,
-        rhoInterp,
+        rhoInterpolator(),
-        UInterp,
+        UInterpolator(),
-        muInterp,
+        muInterpolator(),
-        TInterp,
+        TInterpolator(),
-        cpInterp,
+        cpInterpolator(),
-        pInterp,
+        pInterpolator(),
         this->g().value()
     );
 
@@ -176,7 +211,7 @@ void Foam::ReactingCloud<ParcelType>::inject
     TrackingData& td
 )
 {
-    scalar time = this->runTime().value();
+    scalar time = this->db().time().value();
 
     scalar pRho = td.constProps().rho0();
 
@@ -213,7 +248,7 @@ void Foam::ReactingCloud<ParcelType>::inject
     // Duration of injection period during this timestep
     scalar deltaT = min
     (
-        this->runTime().deltaT().value(),
+        this->db().time().deltaT().value(),
         min
         (
             time - this->injection().timeStart(),
@@ -284,8 +319,8 @@ void Foam::ReactingCloud<ParcelType>::inject
 
             scalar dt = time - timeInj;
 
-            pPtr->stepFraction() = (this->runTime().deltaT().value() - dt)
+            pPtr->stepFraction() = (this->db().time().deltaT().value() - dt)
-                /this->runTime().deltaT().value();
+                /this->db().time().deltaT().value();
 
             injectParcel(td, pPtr);
          }

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

@@ -27,7 +27,6 @@ Class
 
 Description
     Templated base class for reactive cloud
-
     - Adds to kinematic cloud
       - Heat transfer
 
@@ -182,13 +181,15 @@ public:
                     inline PtrList<DimensionedField<scalar, volMesh> >&
                         rhoTrans();
 
-                    //- Retun tmp mass source for field i
+                    //- Return tmp mass source for field i
+                    //  Fully explicit
                     inline tmp<DimensionedField<scalar, volMesh> >
-                        Srho(const label i) const;
+                        Srho1(const label i) const;
 
-                    //- Retun tmp total mass source for carrier phase
+                    //- Return tmp total mass source for carrier phase
+                    //  Fully explicit
                     inline tmp<DimensionedField<scalar, volMesh> >
-                        Srho() const;
+                        Srho1() const;
 
 
         // Check

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

@@ -92,15 +92,15 @@ Foam::ReactingCloud<ParcelType>::rhoTrans()
 
 template<class ParcelType>
 inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
-Foam::ReactingCloud<ParcelType>::Srho(const label i) const
+Foam::ReactingCloud<ParcelType>::Srho1(const label i) const
 {
-    return rhoTrans_[i]/(this->runTime().deltaT()*this->mesh().V());
+    return rhoTrans_[i]/(this->db().time().deltaT()*this->mesh().V());
 }
 
 
 template<class ParcelType>
 inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
-Foam::ReactingCloud<ParcelType>::Srho() const
+Foam::ReactingCloud<ParcelType>::Srho1() const
 {
     tmp<DimensionedField<scalar, volMesh> > trhoTrans
     (
@@ -108,8 +108,8 @@ Foam::ReactingCloud<ParcelType>::Srho() const
         (
             IOobject
             (
-                this->cloudName() + "rhoTrans",
+                this->name() + "rhoTrans",
-                this->runTime().timeName(),
+                this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
@@ -120,13 +120,13 @@ Foam::ReactingCloud<ParcelType>::Srho() const
         )
     );
 
-    DimensionedField<scalar, volMesh>& sourceField = trhoTrans();
+    scalarField& sourceField = trhoTrans().field();
     forAll (rhoTrans_, i)
     {
         sourceField += rhoTrans_[i];
     }
 
-    return trhoTrans/(this->runTime().deltaT()*this->mesh().V());
+    return trhoTrans/(this->db().time().deltaT()*this->mesh().V());
 }
 
 // ************************************************************************* //

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

@@ -22,7 +22,7 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-InClass
+Class
     Foam::cloudThermoTypes
 
 Description

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

@@ -63,13 +63,21 @@ Foam::ThermoCloud<ParcelType>::ThermoCloud
             *this
         )
     ),
+    TIntegrator_
+    (
+        scalarIntegrationScheme::New
+        (
+            "T",
+            this->particleProperties().subDict("integrationSchemes")
+        )
+    ),
     radiation_(this->particleProperties().lookup("radiation")),
     hTrans_
     (
         IOobject
         (
-            this->cloudName() + "hTrans",
+            this->name() + "hTrans",
-            this->runTime().timeName(),
+            this->db().time().timeName(),
             this->db(),
             IOobject::NO_READ,
             IOobject::NO_WRITE,
@@ -82,8 +90,8 @@ Foam::ThermoCloud<ParcelType>::ThermoCloud
     (
         IOobject
         (
-            this->cloudName() + "hCoeff",
+            this->name() + "hCoeff",
-            this->runTime().timeName(),
+            this->db().time().timeName(),
             this->db(),
             IOobject::NO_READ,
             IOobject::NO_WRITE,
@@ -119,21 +127,50 @@ void Foam::ThermoCloud<ParcelType>::evolve()
     const volScalarField& T = carrierThermo_.T();
     const volScalarField cp = carrierThermo_.Cp();
 
-    interpolationCellPoint<scalar> rhoInterp(this->vpi(), this->rho());
+    autoPtr<interpolation<scalar> > rhoInterpolator = interpolation<scalar>::New
-    interpolationCellPoint<vector> UInterp(this->vpi(), this->U());
+    (
-    interpolationCellPoint<scalar> muInterp(this->vpi(), this->mu());
+        this->interpolationSchemes(),
-    interpolationCellPoint<scalar> TInterp(this->vpi(), T);
+        this->vpi(),
-    interpolationCellPoint<scalar> cpInterp(this->vpi(), cp);
+        this->rho()
+    );
+
+    autoPtr<interpolation<vector> > UInterpolator = interpolation<vector>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        this->U()
+    );
+
+    autoPtr<interpolation<scalar> > muInterpolator = interpolation<scalar>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        this->mu()
+    );
+
+    autoPtr<interpolation<scalar> > TInterpolator = interpolation<scalar>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        T
+    );
+
+    autoPtr<interpolation<scalar> > cpInterpolator = interpolation<scalar>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        cp
+    );
 
     typename ParcelType::trackData td
     (
         *this,
         constProps_,
-        rhoInterp,
+        rhoInterpolator(),
-        UInterp,
+        UInterpolator(),
-        muInterp,
+        muInterpolator(),
-        TInterp,
+        TInterpolator(),
-        cpInterp,
+        cpInterpolator(),
         this->g().value()
     );
 
@@ -165,123 +202,8 @@ void Foam::ThermoCloud<ParcelType>::inject
     TrackingData& td
 )
 {
+    // Injection is same as for KinematicCloud<ParcelType>
     KinematicCloud<ParcelType>::inject(td);
-/*
-    scalar time = this->runTime().value();
-
-    scalar pRho = td.constProps().rho0();
-
-    // Number of parcels to introduce during this timestep
-    const label nParcels = this->injection().nParcelsToInject
-    (
-        this->nInjections(),
-        this->time0(),
-        time
-    );
-
-    // Return if no parcels are required
-    if (!nParcels)
-    {
-        this->postInjectCheck();
-        return;
-    }
-
-    // Volume of particles to introduce during this timestep
-    scalar pVolume = this->injection().volume
-    (
-         this->time0(),
-         time,
-         this->meshInfo()
-    );
-
-    // Volume fraction to introduce during this timestep
-    scalar pVolumeFraction =
-        this->injection().volumeFraction(this->time0(), time);
-
-    // Duration of injection period during this timestep
-    scalar deltaT = min
-    (
-        this->runTime().deltaT().value(),
-        min
-        (
-            time - this->injection().timeStart(),
-            this->injection().timeEnd() - this->time0()
-        )
-    );
-
-    // Pad injection time if injection starts during this timestep
-    scalar padTime = max
-    (
-        0.0,
-        this->injection().timeStart() - this->time0()
-    );
-
-    // Introduce new parcels linearly with time
-    for (label iParcel=0; iParcel<nParcels; iParcel++)
-    {
-        // Calculate the pseudo time of injection for parcel 'iParcel'
-        scalar timeInj = this->time0() + padTime + deltaT*iParcel/nParcels;
-
-        // Determine injected parcel properties
-        vector pPosition = this->injection().position
-        (
-            iParcel,
-            timeInj,
-            this->meshInfo(),
-            this->rndGen()
-        );
-
-        // Diameter of parcels
-        scalar pDiameter = this->injection().d0(iParcel, timeInj);
-
-        // Number of particles per parcel
-        scalar pNumberOfParticles = this->setNumberOfParticles
-        (
-            nParcels,
-            pDiameter,
-            pVolumeFraction,
-            pRho,
-            pVolume
-        );
-
-        // Velocity of parcels
-        vector pU = this->injection().velocity(iParcel, timeInj);
-
-        // Determine the injection cell
-        label pCell = -1;
-        this->setInjectorCellAndPosition(pCell, pPosition);
-
-        if (pCell >= 0)
-        {
-            // construct the parcel that is to be injected
-            ParcelType* pPtr = new ParcelType
-            (
-                td.cloud(),
-                this->parcelTypeId(),
-                pPosition,
-                pCell,
-                pDiameter,
-                pU,
-                pNumberOfParticles,
-                td.constProps()
-            );
-
-            scalar dt = time - timeInj;
-
-            pPtr->stepFraction() = (this->runTime().deltaT().value() - dt)
-                /this->runTime().deltaT().value();
-
-            this->injectParcel(td, pPtr);
-         }
-    }
-
-    this->postInjectCheck();
-
-    if (debug)
-    {
-        this->dumpParticlePositions();
-    }
-*/
 }
 
 

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

@@ -27,7 +27,6 @@ Class
 
 Description
     Templated base class for thermodynamic cloud
-
     - Adds to kinematic cloud
       - Heat transfer
 
@@ -87,6 +86,12 @@ class ThermoCloud
                 heatTransferModel_;
 
 
+        // Reference to the particle integration schemes
+
+            //- Temperature integration
+            autoPtr<scalarIntegrationScheme> TIntegrator_;
+
+
         // Modelling options
 
             //- Include radiation
@@ -154,6 +159,12 @@ public:
                     heatTransfer() const;
 
 
+            // Integration schemes
+
+                //-Return reference to velocity integration
+                inline const scalarIntegrationScheme& TIntegrator() const;
+
+
             // Modelling options
 
                 //- Radiation flag
@@ -170,8 +181,11 @@ public:
                     //- Coefficient for carrier phase h equation
                     inline DimensionedField<scalar, volMesh>& hCoeff();
 
-                    //- Enthalpy source term
+                    //- return tmp enthalpy source term - fully explicit
-                    inline tmp<fvScalarMatrix> Sh(volScalarField& h) const;
+                    inline tmp<DimensionedField<scalar, volMesh> > Sh1() const;
+
+                    //- Return tmp enthalpy source term - semi-implicit
+                    inline tmp<fvScalarMatrix> Sh2(volScalarField& h) const;
 
 
                 // Radiation - overrides thermoCloud virtual abstract members

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

@@ -52,6 +52,14 @@ Foam::ThermoCloud<ParcelType>::heatTransfer() const
 }
 
 
+template<class ParcelType>
+inline const Foam::scalarIntegrationScheme&
+Foam::ThermoCloud<ParcelType>::TIntegrator() const
+{
+    return TIntegrator_;
+}
+
+
 template<class ParcelType>
 inline const bool Foam::ThermoCloud<ParcelType>::radiation() const
 {
@@ -75,18 +83,55 @@ Foam::ThermoCloud<ParcelType>::hCoeff()
 }
 
 
+template<class ParcelType>
+inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::ThermoCloud<ParcelType>::Sh1() const
+{
+    tmp<DimensionedField<scalar, volMesh> > tSh1
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                this->name() + "Sh1",
+                this->db().time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            this->mesh(),
+            dimensionedScalar
+            (
+                "zero",
+                dimMass/dimLength/pow3(dimTime),
+                0.0
+            )
+        )
+    );
+
+    scalarField& Sh1 = tSh1().field();
+    Sh1 = hTrans_/(this->mesh().V()*this->db().time().deltaT());
+
+    return tSh1;
+}
+
+
 template<class ParcelType>
 inline Foam::tmp<Foam::fvScalarMatrix>
-Foam::ThermoCloud<ParcelType>::Sh(volScalarField& h) const
+Foam::ThermoCloud<ParcelType>::Sh2(volScalarField& h) const
 {
     const volScalarField cp = carrierThermo_.Cp();
 
-    Info<< "hTrans min/max = "
+    if (debug)
-        << min(hTrans_) << ", " << max(hTrans_) << endl;
+    {
-    Info<< "hCoeff min/max = "
+        Info<< "hTrans min/max = "
-        << min(hCoeff_) << ", " << max(hCoeff_) << endl;
+            << min(hTrans_) << ", " << max(hTrans_) << endl;
+        Info<< "hCoeff min/max = "
+            << min(hCoeff_) << ", " << max(hCoeff_) << endl;
+    }
 
-    return hTrans_/(this->mesh().V()*this->runTime().deltaT())
+    return hTrans_/(this->mesh().V()*this->db().time().deltaT())
          - fvm::Sp(hCoeff_/(cp*this->mesh().V()), h)
          + hCoeff_/(cp*this->mesh().V())*h;
 }
@@ -102,30 +147,34 @@ Foam::ThermoCloud<ParcelType>::Ep() const
         (
             IOobject
             (
-                this->cloudName() + "radiationEp",
+                this->name() + "radiationEp",
-                this->runTime().timeName(),
+                this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
             this->mesh(),
-            dimensionedScalar("zero", dimensionSet(1, -1, -3, 0, 0), 0.0)
+            dimensionedScalar("zero", dimMass/dimLength/pow3(dimTime), 0.0)
         )
     );
 
-    scalarField& Ep = tEp().internalField();
+    // Need to check if coupled as field is created on-the-fly
-    const scalarField& V = this->pMesh().cellVolumes();
+    if (radiation_ && this->coupled())
-    const scalar epsilon = constProps_.epsilon0();
-
-    forAllConstIter(typename ThermoCloud<ParcelType>, *this, iter)
     {
-        const ParcelType& p = iter();
+        scalarField& Ep = tEp().internalField();
-        const label cellI = p.cell();
+        const scalarField& V = this->mesh().V();
-        Ep[cellI] += p.nParticle()*p.areaP()*pow4(p.T());
+        const scalar epsilon = constProps_.epsilon0();
-    }
 
-    Ep *= epsilon*radiation::sigmaSB.value()/(V*mathematicalConstant::pi);
+        forAllConstIter(typename ThermoCloud<ParcelType>, *this, iter)
+        {
+            const ParcelType& p = iter();
+            const label cellI = p.cell();
+            Ep[cellI] += p.nParticle()*p.areaP()*pow4(p.T());
+        }
+
+        Ep *= epsilon*radiation::sigmaSB.value()/V;
+    }
 
     return tEp;
 }
@@ -141,30 +190,34 @@ Foam::ThermoCloud<ParcelType>::ap() const
         (
             IOobject
             (
-                this->cloudName() + "radiationAp",
+                this->name() + "radiationAp",
-                this->runTime().timeName(),
+                this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
             this->mesh(),
-            dimensionedScalar("zero", dimensionSet(0, -1, 0, 0, 0), 0.0)
+            dimensionedScalar("zero", dimless/dimLength, 0.0)
         )
     );
 
-    scalarField& ap = tap().internalField();
+    // Need to check if coupled as field is created on-the-fly
-    const scalarField& V = this->pMesh().cellVolumes();
+    if (radiation_ && this->coupled())
-    const scalar epsilon = constProps_.epsilon0();
-
-    forAllConstIter(typename ThermoCloud<ParcelType>, *this, iter)
     {
-        const ParcelType& p = iter();
+        scalarField& ap = tap().internalField();
-        const label cellI = p.cell();
+        const scalarField& V = this->mesh().V();
-        ap[cellI] += p.nParticle()*p.areaP();
+        const scalar epsilon = constProps_.epsilon0();
-    }
 
-    ap *= epsilon/V;
+        forAllConstIter(typename ThermoCloud<ParcelType>, *this, iter)
+        {
+            const ParcelType& p = iter();
+            const label cellI = p.cell();
+            ap[cellI] += p.nParticle()*p.areaP();
+        }
+
+        ap *= epsilon/V;
+    }
 
     return tap;
 }
@@ -180,32 +233,36 @@ Foam::ThermoCloud<ParcelType>::sigmap() const
         (
             IOobject
             (
-                this->cloudName() + "radiationSigmap",
+                this->name() + "radiationSigmap",
-                this->runTime().timeName(),
+                this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
             this->mesh(),
-            dimensionedScalar("zero", dimensionSet(0, -1, 0, 0, 0), 0.0)
+            dimensionedScalar("zero", dimless/dimLength, 0.0)
         )
     );
 
-    scalarField& sigmap = tsigmap().internalField();
+    // Need to check if coupled as field is created on-the-fly
-
+    if (radiation_ && this->coupled())
-    const scalarField& V = this->pMesh().cellVolumes();
-    const scalar epsilon = constProps_.epsilon0();
-    const scalar f = constProps_.f0();
-
-    forAllConstIter(typename ThermoCloud<ParcelType>, *this, iter)
     {
-        const ParcelType& p = iter();
+        scalarField& sigmap = tsigmap().internalField();
-        const label cellI = p.cell();
-        sigmap[cellI] += p.nParticle()*p.areaP();
-    }
 
-    sigmap *= (1.0 - f)*(1.0 - epsilon)/V;
+        const scalarField& V = this->mesh().V();
+        const scalar epsilon = constProps_.epsilon0();
+        const scalar f = constProps_.f0();
+
+        forAllConstIter(typename ThermoCloud<ParcelType>, *this, iter)
+        {
+            const ParcelType& p = iter();
+            const label cellI = p.cell();
+            sigmap[cellI] += p.nParticle()*p.areaP();
+        }
+
+        sigmap *= (1.0 - f)*(1.0 - epsilon)/V;
+    }
 
     return tsigmap;
 }

src/lagrangian/intermediate/clouds/derived/basicReactingCloud/basicReactingCloud.H

@@ -26,7 +26,6 @@ Class
     Foam::basicReactingCloud
 
 Description
-    The basicReactingCloud
 
 SourceFiles
     basicReactingCloud.C

src/lagrangian/intermediate/clouds/derived/basicThermoCloud/basicThermoCloud.H

@@ -26,7 +26,6 @@ Class
     Foam::basicThermoCloud
 
 Description
-    The basicThermoCloud
 
 SourceFiles
     basicThermoCloud.C

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

@@ -30,29 +30,58 @@ License
 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
 
 template<class ParcelType>
-template<class TrackingData>
+void Foam::KinematicParcel<ParcelType>::updateCellQuantities
+(
+    trackData& td,
+    const scalar dt,
+    const label celli
+)
+{
+    rhoc_ = td.rhoInterp().interpolate(this->position(), celli);
+    Uc_ = td.UInterp().interpolate(this->position(), celli);
+    muc_ = td.muInterp().interpolate(this->position(), celli);
+
+    // Apply dispersion components to carrier phase velocity
+    Uc_ = td.cloud().dispersion().update
+    (
+        dt,
+        celli,
+        U_,
+        Uc_,
+        UTurb_,
+        tTurb_
+    );
+}
+
+
+template<class ParcelType>
 void Foam::KinematicParcel<ParcelType>::calcCoupled
 (
-    TrackingData& td,
+    trackData& td,
-    const label celli,
     const scalar dt,
-    const scalar rhoc,
+    const label celli
-    vector& Uc,
-    const scalar muc
 )
 {
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Define local properties at beginning of timestep
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    const scalar mass0 = mass();
+//    const scalar mass0 = mass();
-    const vector U0 = U_;
+//    const vector U0 = U_;
+
+
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Initialise transfer terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Momentum transfer from the particle to the carrier phase
+    vector dUTrans = vector::zero;
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate velocity - update U
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar Cud = 0.0;
-    const vector U1 = calcVelocity(td, dt, rhoc, Uc, muc, Cud);
+    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -60,10 +89,10 @@ void Foam::KinematicParcel<ParcelType>::calcCoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~
 
     // Update momentum transfer
-    td.cloud().UTrans()[celli] += nParticle_*mass0*(U0 - U1);
+    td.cloud().UTrans()[celli] += nParticle_*dUTrans;
 
     // Accumulate coefficient to be applied in carrier phase momentum coupling
-    td.cloud().UCoeff()[celli] += nParticle_*mass0*Cud;
+    td.cloud().UCoeff()[celli] += nParticle_*mass()*Cud;
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -74,72 +103,62 @@ void Foam::KinematicParcel<ParcelType>::calcCoupled
 
 
 template<class ParcelType>
-template<class TrackingData>
 void Foam::KinematicParcel<ParcelType>::calcUncoupled
 (
-    TrackingData& td,
+    trackData& td,
-    const scalar dt,
+    const scalar dt
-    const scalar rhoc,
-    vector& Uc,
-    const scalar muc
 )
 {
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Initialise transfer terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Momentum transfer from the particle to the carrier phase
+    vector dUTrans = vector::zero;
+
+
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate velocity - update U
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar Cud = 0.0;
-    this->U() = calcVelocity(td, dt, rhoc, Uc, muc, Cud);
+    this->U() = calcVelocity(td, dt, Cud, dUTrans);
 }
 
 
 template<class ParcelType>
-template<class TrackingData>
 Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
 (
-    TrackingData& td,
+    trackData& td,
     const scalar dt,
-    const scalar rhoc,
+    scalar& Cud,
-    vector& Uc,
+    vector& dUTrans
-    const scalar muc,
-    scalar& Cud
 )
 {
     // Correct carrier phase velocity for 2-D slab cases
     const polyMeshInfo& meshInfo = td.cloud().meshInfo();
     if (meshInfo.caseIs2dSlab())
     {
-        Uc.component(meshInfo.emptyComponent()) = 0.0;
+        Uc_.component(meshInfo.emptyComponent()) = 0.0;
     }
 
-    // Update relative velocity
-    Ur_ = U_ - Uc;
-
     // Return linearised term from drag model
-//    const scalar Cud = td.cloud().drag().Cu
+    Cud = td.cloud().drag().Cu(U_ - Uc_, d_, rhoc_, rho_, muc_);
-    Cud = td.cloud().drag().Cu
-    (
-        Ur_,
-        d_,
-        rhoc,
-        rho_,
-        muc
-    );
-
-    // Update velocity - treat as 3-D
-    const vector ap = (1.0 - rhoc/rho_)*td.g();
-    const scalar bp = 1.0/Cud;
 
 
     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Set new particle velocity
     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Analytical
+    // Update velocity - treat as 3-D
-//    const scalar expTerm = exp(-dt/bp);
+    const scalar bp = 1.0/(Cud + VSMALL);
-//    vector Unew = Uc + (U_ - Uc)*expTerm + ap*bp*(1.0 - expTerm);
+    const vector ap = Uc_/bp + rhoc_/rho_*td.g();
 
-    // Euler-implicit
+    vector Unew = td.cloud().UIntegrator().integrate(U_, dt, ap, bp);
-    vector Unew = (U_ + dt*(ap + Uc/bp))/(1.0 + dt/bp);
+
+//    Info<< "U_, Unew = " << U_ << ", " << Unew << endl;
+
+    // Calculate the momentum transfer to the continuous phase
+    dUTrans = -mass()*(Unew - U_);
 
     // Make corrections for 2-D cases
     if (meshInfo.caseIs2d())
@@ -148,6 +167,7 @@ Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
         {
             // Remove the slab normal parcel velocity component
             Unew.component(meshInfo.emptyComponent()) = 0.0;
+            dUTrans.component(meshInfo.emptyComponent()) = 0.0;
 
             // Snap parcels to central plane
             this->position().component(meshInfo.emptyComponent()) =
@@ -173,10 +193,9 @@ Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template <class ParcelType>
-template <class TrackingData>
 bool Foam::KinematicParcel<ParcelType>::move
 (
-    TrackingData& td
+    trackData& td
 )
 {
     td.switchProcessor = false;
@@ -198,39 +217,21 @@ bool Foam::KinematicParcel<ParcelType>::move
         // since this will change if a face is hit
         label celli = this->cell();
 
-        dt *= trackToFace(this->position() + dt*U_, td);
+        dt *= this->trackToFace(this->position() + dt*U_, td);
 
         tEnd -= dt;
         this->stepFraction() = 1.0 - tEnd/deltaT;
 
-        cellPointWeight cpw
+        // Update cell based properties
-        (
+        updateCellQuantities(td, dt, celli);
-            mesh,
-            this->position(),
-            celli,
-            faceInterpolation()
-        );
-        scalar rhoc = td.rhoInterp().interpolate(cpw);
-        vector Uc = td.UInterp().interpolate(cpw);
-        scalar muc = td.muInterp().interpolate(cpw);
-
-        Uc = td.cloud().dispersion().update
-        (
-            dt,
-            celli,
-            U_,
-            Uc,
-            UTurb_,
-            tTurb_
-        );
 
         if (td.cloud().coupled())
         {
-            calcCoupled(td, celli, dt, rhoc, Uc, muc);
+            calcCoupled(td, dt, celli);
         }
         else
         {
-            calcUncoupled(td, dt, rhoc, Uc, muc);
+            calcUncoupled(td, dt);
         }
 
         if (this->onBoundary() && td.keepParticle)
@@ -254,11 +255,11 @@ bool Foam::KinematicParcel<ParcelType>::move
 
 
 template <class ParcelType>
-template <class TrackingData>
+template <class TrackData>
 void Foam::KinematicParcel<ParcelType>::hitProcessorPatch
 (
     const processorPolyPatch&,
-    TrackingData& td
+    TrackData& td
 )
 {
     td.switchProcessor = true;
@@ -275,11 +276,11 @@ void Foam::KinematicParcel<ParcelType>::hitProcessorPatch
 
 
 template <class ParcelType>
-template <class TrackingData>
+template <class TrackData>
 void Foam::KinematicParcel<ParcelType>::hitWallPatch
 (
     const wallPolyPatch& wpp,
-    TrackingData& td
+    TrackData& td
 )
 {
     td.cloud().wallInteraction().correct(wpp, this->face(), U_);
@@ -296,11 +297,11 @@ void Foam::KinematicParcel<ParcelType>::hitWallPatch
 
 
 template <class ParcelType>
-template <class TrackingData>
+template <class TrackData>
 void Foam::KinematicParcel<ParcelType>::hitPatch
 (
     const polyPatch&,
-    TrackingData& td
+    TrackData& td
 )
 {
     td.keepParticle = false;

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

@@ -26,10 +26,8 @@ Class
     Foam::KinematicParcel
 
 Description
-    Kinematic parcel class with one-way coupling with the continuous
+    Kinematic parcel class with one/two-way coupling with the continuous
-    phase.
+    phase. Sub-models include:
-
-    Sub-models include:
     - drag
     - break-up
     - wall interactions
@@ -61,6 +59,8 @@ namespace Foam
 template<class ParcelType>
 class KinematicParcel;
 
+// Forward declaration of friend functions
+
 template<class ParcelType>
 Ostream& operator<<
 (
@@ -78,80 +78,9 @@ class KinematicParcel
     public Particle<ParcelType>
 {
 
-    // Private member data
-
-        //- Parcel type id
-        label typeId_;
-
-        //- Diameter [m]
-        scalar d_;
-
-        //- Velocity of Parcel [m/s]
-        vector U_;
-
-        //- Relative velocity of Parcel [m/s]
-        vector Ur_;
-
-        //- Number of particles in Parcel
-        scalar nParticle_;
-
-        //- Density [kg/m3]
-        scalar rho_;
-
-        //- Time spent in turbulent eddy
-        scalar tTurb_;
-
-        //- Turbulent velocity fluctuation
-        vector UTurb_;
-
-
-protected:
-
-    // Protected member functions
-
-        template<class TrackingData>
-        void calcCoupled
-        (
-            TrackingData& td,
-            const label celli,
-            const scalar dt,
-            const scalar rhoc,
-            vector& Uc,
-            const scalar muc
-        );
-
-        template<class TrackingData>
-        void calcUncoupled
-        (
-            TrackingData& td,
-            const scalar dt,
-            const scalar rhoc,
-            vector& Uc,
-            const scalar muc
-        );
-
-        //- Calculate new particle velocity
-        template<class TrackingData>
-        vector calcVelocity
-        (
-            TrackingData& td,
-            const scalar dt,
-            const scalar rhoc,
-            vector& Uc,
-            const scalar muc,
-            scalar& Cud
-        );
-
-
 public:
 
-    //- Runtime type information
+    //- Class to hold kinematic particle constant properties
-    TypeName("KinematicParcel");
-
-    friend class Cloud<ParcelType>;
-
-
-    //- Class to hold particle constant properties
     class constantProperties
     {
 
@@ -182,7 +111,7 @@ public:
     };
 
 
-    //- Class used to pass tracking data to the trackToFace function
+    //- Class used to pass kinematic tracking data to the trackToFace function
     class trackData
     :
         public Particle<ParcelType>::trackData
@@ -198,9 +127,9 @@ public:
 
             //- Interpolators for continuous phase fields
 
-                const interpolationCellPoint<scalar>& rhoInterp_;
+                const interpolation<scalar>& rhoInterp_;
-                const interpolationCellPoint<vector>& UInterp_;
+                const interpolation<vector>& UInterp_;
-                const interpolationCellPoint<scalar>& muInterp_;
+                const interpolation<scalar>& muInterp_;
 
             //- Local gravitational or other body-force acceleration
             const vector& g_;
@@ -218,9 +147,9 @@ public:
             (
                 KinematicCloud<ParcelType>& cloud,
                 const constantProperties& constProps,
-                const interpolationCellPoint<scalar>& rhoInterp,
+                const interpolation<scalar>& rhoInterp,
-                const interpolationCellPoint<vector>& UInterp,
+                const interpolation<vector>& UInterp,
-                const interpolationCellPoint<scalar>& muInterp,
+                const interpolation<scalar>& muInterp,
                 const vector& g
             );
 
@@ -231,16 +160,96 @@ public:
 
             inline const constantProperties& constProps() const;
 
-            inline const interpolationCellPoint<scalar>& rhoInterp() const;
+            inline const interpolation<scalar>& rhoInterp() const;
 
-            inline const interpolationCellPoint<vector>& UInterp() const;
+            inline const interpolation<vector>& UInterp() const;
 
-            inline const interpolationCellPoint<scalar>& muInterp() const;
+            inline const interpolation<scalar>& muInterp() const;
 
             inline const vector& g() const;
     };
 
 
+protected:
+
+    // Protected member data
+
+        // Parcel properties
+
+            //- Parcel type id
+            label typeId_;
+
+            //- Diameter [m]
+            scalar d_;
+
+            //- Velocity of Parcel [m/s]
+            vector U_;
+
+            //- Number of particles in Parcel
+            scalar nParticle_;
+
+            //- Density [kg/m3]
+            scalar rho_;
+
+            //- Time spent in turbulent eddy [s]
+            scalar tTurb_;
+
+            //- Turbulent velocity fluctuation [m/s]
+            vector UTurb_;
+
+
+        // Cell-based quantities
+
+            // - Density
+            scalar rhoc_;
+
+            // - Velocity
+            vector Uc_;
+
+            // - Viscosity
+            scalar muc_;
+
+
+    // Protected member functions
+
+        virtual void updateCellQuantities
+        (
+            trackData& td,
+            const scalar dt,
+            const label celli
+        );
+
+        virtual void calcCoupled
+        (
+            trackData& td,
+            const scalar dt,
+            const label celli
+        );
+
+        virtual void calcUncoupled
+        (
+            trackData& td,
+            const scalar dt
+        );
+
+        //- Calculate new particle velocity
+        virtual vector calcVelocity
+        (
+            trackData& td,
+            const scalar dt,
+            scalar& Cud,
+            vector& dUTrans
+        );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("KinematicParcel");
+
+    friend class Cloud<ParcelType>;
+
+
     // Constructors
 
         //- Construct from components
@@ -276,10 +285,10 @@ public:
         // Access
 
             //- Return type id
-            inline const label& typeId() const;
+            inline const label typeId() const;
 
             //- Return diameter
-            inline const scalar& d() const;
+            inline const scalar d() const;
             inline scalar& d();
 
             //- Return velocity
@@ -291,15 +300,15 @@ public:
             inline vector& Ur();
 
             //- Return number of particles
-            inline const scalar& nParticle() const;
+            inline const scalar nParticle() const;
             inline scalar& nParticle();
 
             //- Return density
-            inline const scalar& rho() const;
+            inline const scalar rho() const;
             inline scalar& rho();
 
             //- Return time spent in turbulent eddy
-            inline const scalar& tTurb() const;
+            inline const scalar tTurb() const;
             inline scalar& tTurb();
 
             //- Return turbulent velocity fluctuation
@@ -310,7 +319,8 @@ public:
             //  the particle can be in the n direction
             inline scalar wallImpactDistance(const vector& n) const;
 
-            //- Return the index of the face to be used in the interpolation routine
+            //- Return the index of the face to be used in the interpolation
+            //  routine
             inline label faceInterpolation() const;
 
             //- Particle volume
@@ -329,51 +339,53 @@ public:
         // Tracking
 
             //- Move the parcel
-            template<class TrackingData>
+            virtual bool move
-            bool move
             (
-                TrackingData& td
+                trackData& td
             );
 
+
+        // Patch interactions
+
             //- Overridable function to handle the particle hitting a
             //  processorPatch
-            template<class TrackingData>
+            template<class TrackData>
             void hitProcessorPatch
             (
                 const processorPolyPatch&,
-                TrackingData& td
+                TrackData& td
             );
 
             //- Overridable function to handle the particle hitting a
             //  processorPatch without trackData
-            virtual void hitProcessorPatch
+            void hitProcessorPatch
             (
                 const processorPolyPatch&,
                 int&
             );
 
             //- Overridable function to handle the particle hitting a wallPatch
-            template<class TrackingData>
+            template<class TrackData>
             void hitWallPatch
             (
                 const wallPolyPatch&,
-                TrackingData& td
+                TrackData& td
             );
 
             //- Overridable function to handle the particle hitting a wallPatch
             //  without trackData
-            virtual void hitWallPatch
+            void hitWallPatch
             (
                 const wallPolyPatch&,
                 int&
             );
 
             //- Overridable function to handle the particle hitting a polyPatch
-            template<class TrackingData>
+            template<class TrackData>
             void hitPatch
             (
                 const polyPatch&,
-                TrackingData& td
+                TrackData& td
             );
 
             //- Overridable function to handle the particle hitting a polyPatch
@@ -386,11 +398,11 @@ public:
 
             //- Transform the physical properties of the particle
             //  according to the given transformation tensor
-            void transformProperties(const tensor& T);
+            virtual void transformProperties(const tensor& T);
 
             //- Transform the physical properties of the particle
             //  according to the given separation vector
-            void transformProperties(const vector& separation);
+            virtual void transformProperties(const vector& separation);
 
 
         // I-O
@@ -405,6 +417,7 @@ public:
                 const KinematicCloud<ParcelType>& c
             );
 
+
     // Ostream Operator
 
         friend Ostream& operator<< <ParcelType>

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

@@ -42,9 +42,9 @@ inline Foam::KinematicParcel<ParcelType>::trackData::trackData
 (
     KinematicCloud<ParcelType>& cloud,
     const constantProperties& constProps,
-    const interpolationCellPoint<scalar>& rhoInterp,
+    const interpolation<scalar>& rhoInterp,
-    const interpolationCellPoint<vector>& UInterp,
+    const interpolation<vector>& UInterp,
-    const interpolationCellPoint<scalar>& muInterp,
+    const interpolation<scalar>& muInterp,
     const vector& g
 )
 :
@@ -75,11 +75,13 @@ inline Foam::KinematicParcel<ParcelType>::KinematicParcel
     typeId_(typeId),
     d_(d0),
     U_(U0),
-    Ur_(vector::zero),
     nParticle_(nParticle0),
     rho_(constProps.rho0()),
     tTurb_(0.0),
-    UTurb_(vector::zero)
+    UTurb_(vector::zero),
+    rhoc_(0.0),
+    Uc_(vector::zero),
+    muc_(0.0)
 {}
 
 
@@ -120,7 +122,7 @@ Foam::KinematicParcel<ParcelType>::trackData::constProps() const
 
 
 template<class ParcelType>
-inline const Foam::interpolationCellPoint<Foam::scalar>&
+inline const Foam::interpolation<Foam::scalar>&
 Foam::KinematicParcel<ParcelType>::trackData::rhoInterp() const
 {
     return rhoInterp_;
@@ -128,7 +130,7 @@ Foam::KinematicParcel<ParcelType>::trackData::rhoInterp() const
 
 
 template <class ParcelType>
-inline const Foam::interpolationCellPoint<Foam::vector>&
+inline const Foam::interpolation<Foam::vector>&
 Foam::KinematicParcel<ParcelType>::trackData::UInterp() const
 {
     return UInterp_;
@@ -136,7 +138,7 @@ Foam::KinematicParcel<ParcelType>::trackData::UInterp() const
 
 
 template<class ParcelType>
-inline const Foam::interpolationCellPoint<Foam::scalar>&
+inline const Foam::interpolation<Foam::scalar>&
 Foam::KinematicParcel<ParcelType>::trackData::muInterp() const
 {
     return muInterp_;
@@ -154,14 +156,14 @@ Foam::KinematicParcel<ParcelType>::trackData::g() const
 // * * * * * * * * * * KinematicParcel Member Functions  * * * * * * * * * * //
 
 template <class ParcelType>
-inline const Foam::label& Foam::KinematicParcel<ParcelType>::typeId() const
+inline const Foam::label Foam::KinematicParcel<ParcelType>::typeId() const
 {
     return typeId_;
 }
 
 
 template <class ParcelType>
-inline const Foam::scalar& Foam::KinematicParcel<ParcelType>::d() const
+inline const Foam::scalar Foam::KinematicParcel<ParcelType>::d() const
 {
     return d_;
 }
@@ -214,21 +216,7 @@ inline Foam::vector& Foam::KinematicParcel<ParcelType>::U()
 
 
 template <class ParcelType>
-inline const Foam::vector& Foam::KinematicParcel<ParcelType>::Ur() const
+inline const Foam::scalar Foam::KinematicParcel<ParcelType>::nParticle() const
-{
-    return Ur_;
-}
-
-
-template <class ParcelType>
-inline Foam::vector& Foam::KinematicParcel<ParcelType>::Ur()
-{
-    return Ur_;
-}
-
-
-template <class ParcelType>
-inline const Foam::scalar& Foam::KinematicParcel<ParcelType>::nParticle() const
 {
     return nParticle_;
 }
@@ -242,7 +230,7 @@ inline Foam::scalar& Foam::KinematicParcel<ParcelType>::nParticle()
 
 
 template <class ParcelType>
-inline const Foam::scalar& Foam::KinematicParcel<ParcelType>::rho() const
+inline const Foam::scalar Foam::KinematicParcel<ParcelType>::rho() const
 {
     return rho_;
 }
@@ -256,7 +244,7 @@ inline Foam::scalar& Foam::KinematicParcel<ParcelType>::rho()
 
 
 template <class ParcelType>
-inline const Foam::scalar& Foam::KinematicParcel<ParcelType>::tTurb() const
+inline const Foam::scalar Foam::KinematicParcel<ParcelType>::tTurb() const
 {
     return tTurb_;
 }
@@ -286,7 +274,7 @@ inline Foam::vector& Foam::KinematicParcel<ParcelType>::UTurb()
 template <class ParcelType>
 inline Foam::scalar Foam::KinematicParcel<ParcelType>::volume() const
 {
-    return mathematicalConstant::pi/6*pow(d_, 3);
+    return mathematicalConstant::pi/6.0*pow3(d_);
 }
 
 
@@ -300,7 +288,7 @@ inline Foam::scalar Foam::KinematicParcel<ParcelType>::mass() const
 template <class ParcelType>
 inline Foam::scalar Foam::KinematicParcel<ParcelType>::areaP() const
 {
-    return 0.25*mathematicalConstant::pi*d_*d_;
+    return 0.25*areaS();
 }
 
 

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C

@@ -41,11 +41,13 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
     typeId_(0),
     d_(0.0),
     U_(vector::zero),
-    Ur_(vector::zero),
+    nParticle_(0.0),
-    nParticle_(0),
     rho_(0.0),
     tTurb_(0.0),
-    UTurb_(vector::zero)
+    UTurb_(vector::zero),
+    rhoc_(0.0),
+    Uc_(vector::zero),
+    muc_(0.0)
 {
     if (readFields)
     {
@@ -54,8 +56,7 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
             typeId_ = readLabel(is);
             d_ = readScalar(is);
             is >> U_;
-            is >> Ur_;
+            nParticle_ = readScalar(is);
-            nParticle_ = readLabel(is);
             rho_ = readScalar(is);
             tTurb_ = readScalar(is);
             is >> UTurb_;
@@ -68,7 +69,6 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
                 sizeof(typeId_)
               + sizeof(d_)
               + sizeof(U_)
-              + sizeof(Ur_)
               + sizeof(nParticle_)
               + sizeof(rho_)
               + sizeof(tTurb_)
@@ -106,9 +106,6 @@ void Foam::KinematicParcel<ParcelType>::readFields
     IOField<vector> U(c.fieldIOobject("U"));
     c.checkFieldIOobject(c, U);
 
-    IOField<vector> Ur(c.fieldIOobject("Ur"));
-    c.checkFieldIOobject(c, Ur);
-
     IOField<scalar> nParticle(c.fieldIOobject("nParticle"));
     c.checkFieldIOobject(c, nParticle);
 
@@ -129,7 +126,6 @@ void Foam::KinematicParcel<ParcelType>::readFields
         p.typeId_ = typeId[i];
         p.d_ = d[i];
         p.U_ = U[i];
-        p.Ur_ = Ur[i];
         p.nParticle_ = nParticle[i];
         p.rho_ = rho[i];
         p.tTurb_ = tTurb[i];
@@ -152,7 +148,6 @@ void Foam::KinematicParcel<ParcelType>::writeFields
     IOField<label> typeId(c.fieldIOobject("typeId"), np);
     IOField<scalar> d(c.fieldIOobject("d"), np);
     IOField<vector> U(c.fieldIOobject("U"), np);
-    IOField<vector> Ur(c.fieldIOobject("Ur"), np);
     IOField<scalar> nParticle(c.fieldIOobject("nParticle"), np);
     IOField<scalar> rho(c.fieldIOobject("rho"), np);
     IOField<scalar> tTurb(c.fieldIOobject("tTurb"), np);
@@ -166,7 +161,6 @@ void Foam::KinematicParcel<ParcelType>::writeFields
         typeId[i] = p.typeId();
         d[i] = p.d();
         U[i] = p.U();
-        Ur[i] = p.Ur();
         nParticle[i] = p.nParticle();
         rho[i] = p.rho();
         tTurb[i] = p.tTurb();
@@ -177,7 +171,6 @@ void Foam::KinematicParcel<ParcelType>::writeFields
     typeId.write();
     d.write();
     U.write();
-    Ur.write();
     nParticle.write();
     rho.write();
     tTurb.write();
@@ -200,7 +193,6 @@ Foam::Ostream& Foam::operator<<
             << token::SPACE << p.typeId()
             << token::SPACE << p.d()
             << token::SPACE << p.U()
-            << token::SPACE << p.Ur()
             << token::SPACE << p.nParticle()
             << token::SPACE << p.rho()
             << token::SPACE << p.tTurb()
@@ -211,11 +203,10 @@ Foam::Ostream& Foam::operator<<
         os  << static_cast<const Particle<ParcelType>& >(p);
         os.write
         (
-            reinterpret_cast<const char*>(&p.typeId()),
+            reinterpret_cast<const char*>(p.typeId()),
             sizeof(p.typeId())
           + sizeof(p.d())
           + sizeof(p.U())
-          + sizeof(p.Ur())
           + sizeof(p.nParticle())
           + sizeof(p.rho())
           + sizeof(p.tTurb())

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

@@ -29,28 +29,45 @@ License
 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
 
 template<class ParcelType>
-template<class TrackingData>
+void Foam::ReactingParcel<ParcelType>::updateCellQuantities
+(
+    trackData& td,
+    const scalar dt,
+    const label celli
+)
+{
+    ThermoParcel<ParcelType>::updateCellQuantities(td, dt, celli);
+
+    pc_ = td.pInterp().interpolate(this->position(), celli);
+}
+
+
+template<class ParcelType>
 void Foam::ReactingParcel<ParcelType>::calcCoupled
 (
-    TrackingData& td,
+    trackData& td,
-    const label celli,
     const scalar dt,
-    const scalar rhoc,
+    const label celli
-    vector& Uc,
-    const scalar muc,
-    const scalar Tc,
-    const scalar cpc,
-    const scalar pc
 )
 {
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Define local properties at beginning of timestep
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    const vector U0 = this->U();
+    const vector U0 = this->U_;
-    const scalar T0 = this->T();
     const scalar mass0 = this->mass();
-    const scalar cp0 = this->cp();
+    const scalar np0 = this->nParticle_;
-    const scalar np0 = this->nParticle();
+    const scalar T0 = this->T_;
+    const scalar cp0 = this->cp_;
+
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Initialise transfer terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Momentum transfer from the particle to the carrier phase
+    vector dUTrans = vector::zero;
+
+    // Enthalpy transfer from the particle to the carrier phase
+    scalar dhTrans = 0.0;
 
     // Mass transfer from particle to carrier phase
     // - components exist in particle already
@@ -69,14 +86,14 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     // Calculate heat transfer - update T
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar htc = 0.0;
-    scalar T1 = calcHeatTransfer(td, celli, dt, rhoc, Uc, muc, Tc, cpc, htc);
+    scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate velocity - update U
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar Cud = 0.0;
-    const vector U1 = calcVelocity(td, dt, rhoc, Uc, muc, Cud);
+    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -88,18 +105,7 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate surface reactions
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    calcSurfaceReactions
+    calcSurfaceReactions(td, dt, celli, T0, T1, dMassMTSR, dMassSR);
-    (
-        td,
-        dt,
-        celli,
-        rhoc,
-        Tc,
-        T0,
-        T1,
-        dMassMTSR,
-        dMassSR
-    );
 
     // New total mass
     const scalar mass1 = mass0 - sum(dMassMT) - dMassMTSR;
@@ -112,7 +118,7 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     // Specific heat capacity of non-volatile components
     const scalar cpNonVolatile =
         (
-            YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc, Tc)
+            YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, this->Tc_)
           + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_)
         )/(YMixture_[1] + YMixture_[2]);
 
@@ -129,18 +135,18 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     // Transfer mass lost from particle to carrier mass source
     forAll(dMassMT, i)
     {
-        td.cloud().rhoTrans(i)[celli] +=
+        td.cloud().rhoTrans(i)[celli] += np0*(dMassMT[i] + dMassSR[i]);
-            np0*(dMassMT[i] + dMassSR[i]);
     }
 
     // Update momentum transfer
-    td.cloud().UTrans()[celli] += np0*(mass0*U0 - mass1*U1);
+    td.cloud().UTrans()[celli] += np0*dUTrans;
 
     // Accumulate coefficient to be applied in carrier phase momentum coupling
     td.cloud().UCoeff()[celli] += np0*mass0*Cud;
 
     // Update enthalpy transfer
-    td.cloud().hTrans()[celli] += np0*(mass0*cp0*T0 - mass1*cp1*T1);
+//    td.cloud().hTrans()[celli] += np0*(mass0*cp0*T0 - mass1*cp1*T1);
+    td.cloud().hTrans()[celli] += np0*((mass0*cp0 - mass1*cp1)*T0 + dhTrans);
 
     // Accumulate coefficient to be applied in carrier phase enthalpy coupling
     td.cloud().hCoeff()[celli] += np0*htc*this->areaS();
@@ -166,45 +172,48 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     else
     {
-        this->U() = U1;
+        this->U_ = U1;
-        this->T() = T1;
+        this->T_ = T1;
-        this->cp() = cp1;
+        this->cp_ = cp1;
 
         // Update particle density or diameter
         if (td.cloud().massTransfer().changesVolume())
         {
-            this->d() = cbrt(mass1/this->rho()*6.0/mathematicalConstant::pi);
+            this->d_ = cbrt(mass1/this->rho_*6.0/mathematicalConstant::pi);
         }
         else
         {
-            this->rho() = mass1/this->volume();
+            this->rho_ = mass1/this->volume();
         }
     }
 }
 
 
 template<class ParcelType>
-template<class TrackingData>
 void Foam::ReactingParcel<ParcelType>::calcUncoupled
 (
-    TrackingData& td,
+    trackData& td,
-    const label celli,
     const scalar dt,
-    const scalar rhoc,
+    const label celli
-    vector& Uc,
-    const scalar muc,
-    const scalar Tc,
-    const scalar cpc,
-    const scalar pc
 )
 {
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Define local properties at beginning of timestep
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    const scalar T0 = this->T();
+    const scalar T0 = this->T_;
     const scalar mass0 = this->mass();
 //    const scalar cp0 = this->cp();
 
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Initialise transfer terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Momentum transfer from the particle to the carrier phase
+    vector dUTrans = vector::zero;
+
+    // Enthalpy transfer from the particle to the carrier phase
+    scalar dhTrans = 0.0;
+
     // Mass transfer from particle to carrier phase
     // - components exist in particle already
     scalarList dMassMT(td.cloud().gases().size(), 0.0);
@@ -222,7 +231,7 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     // Calculate heat transfer - update T
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar htc = 0.0;
-    scalar T1 = calcHeatTransfer(td, celli, dt, rhoc, Uc, muc, Tc, cpc, htc);
+    scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
 
     // Limit new temp max by vapourisarion temperature
     T1 = min(td.constProps().Tvap(), T1);
@@ -231,8 +240,8 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate velocity - update U
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar Fd = 0.0;
+    scalar Cud = 0.0;
-    const vector U1 = calcVelocity(td, dt, rhoc, Uc, muc, Fd);
+    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -249,8 +258,6 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
         td,
         dt,
         celli,
-        rhoc,
-        Tc,
         T0,
         T1,
         dMassMTSR,
@@ -268,7 +275,7 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     // Specific heat capacity of non-volatile components
     const scalar cpNonVolatile =
         (
-            YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc, Tc)
+            YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, this->Tc_)
           + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_)
         )/(YMixture_[1] + YMixture_[2]);
 
@@ -290,28 +297,27 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     else
     {
-        this->U() = U1;
+        this->U_ = U1;
-        this->T() = T1;
+        this->T_ = T1;
-        this->cp() = cp1;
+        this->cp_ = cp1;
 
         // Update particle density or diameter
         if (td.cloud().massTransfer().changesVolume())
         {
-            this->d() = cbrt(mass1/this->rho()*6.0/mathematicalConstant::pi);
+            this->d_ = cbrt(mass1/this->rho_*6.0/mathematicalConstant::pi);
         }
         else
         {
-            this->rho() = mass1/this->volume();
+            this->rho_ = mass1/this->volume();
         }
     }
 }
 
 
 template<class ParcelType>
-template<class TrackingData>
 void Foam::ReactingParcel<ParcelType>::calcMassTransfer
 (
-    TrackingData& td,
+    trackData& td,
     const scalar dt,
     const scalar T0,
     const scalar T1,
@@ -333,8 +339,8 @@ void Foam::ReactingParcel<ParcelType>::calcMassTransfer
     if
     (
         !td.cloud().massTransfer().active()
-     || this->T()<td.constProps().Tvap()
+     || this->T_<td.constProps().Tvap()
-     || this->T()<td.constProps().Tbp()
+     || this->T_<td.constProps().Tbp()
     )
     {
         return;
@@ -371,14 +377,11 @@ void Foam::ReactingParcel<ParcelType>::calcMassTransfer
 
 
 template<class ParcelType>
-template<class TrackingData>
 void Foam::ReactingParcel<ParcelType>::calcSurfaceReactions
 (
-    TrackingData& td,
+    trackData& td,
     const scalar dt,
     const label celli,
-    const scalar rhoc,
-    const scalar Tc,
     const scalar T0,
     const scalar T1,
     scalar& dMassMTSR,
@@ -397,11 +400,11 @@ void Foam::ReactingParcel<ParcelType>::calcSurfaceReactions
     (
         dt,
         celli,
-        this->d(),
+        this->d_,
         T0,
         T1,
-        Tc,
+        this->Tc_,
-        rhoc,
+        this->rhoc_,
         this->mass(),
         YGas_,
         YLiquid_,
@@ -415,92 +418,6 @@ void Foam::ReactingParcel<ParcelType>::calcSurfaceReactions
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class ParcelType>
-template<class TrackingData>
-bool Foam::ReactingParcel<ParcelType>::move
-(
-    TrackingData& td
-)
-{
-    td.switchProcessor = false;
-    td.keepParticle = true;
-
-    const polyMesh& mesh = td.cloud().pMesh();
-    const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
-
-    const scalar deltaT = mesh.time().deltaT().value();
-    scalar tEnd = (1.0 - this->stepFraction())*deltaT;
-    const scalar dtMax = tEnd;
-
-    while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
-    {
-        // Set the Lagrangian time-step
-        scalar dt = min(dtMax, tEnd);
-
-        // Remember which cell the parcel is in
-        // since this will change if a face is hit
-        label celli = this->cell();
-
-        dt *= trackToFace(this->position() + dt*this->U(), td);
-
-        tEnd -= dt;
-
-        this->stepFraction() = 1.0 - tEnd/deltaT;
-
-        // Avoid div0 in reacting sub-models
-        if (dt < SMALL)
-        {
-            break;
-        }
-
-        cellPointWeight cpw
-        (
-            mesh,
-            this->position(),
-            celli,
-            this->faceInterpolation()
-        );
-        scalar rhoc = td.rhoInterp().interpolate(cpw);
-        vector Uc = td.UInterp().interpolate(cpw);
-        scalar muc = td.muInterp().interpolate(cpw);
-        scalar Tc = td.TInterp().interpolate(cpw);
-        scalar cpc = td.cpInterp().interpolate(cpw);
-        scalar pc = td.pInterp().interpolate(cpw);
-
-        Uc = td.cloud().dispersion().update
-        (
-            dt,
-            celli,
-            this->U(),
-            Uc,
-            this->UTurb(),
-            this->tTurb()
-        );
-
-        if (td.cloud().coupled())
-        {
-            calcCoupled(td, celli, dt, rhoc, Uc, muc, Tc, cpc, pc);
-        }
-        else
-        {
-            calcUncoupled(td, celli, dt, rhoc, Uc, muc, Tc, cpc, pc);
-        }
-
-        if (this->onBoundary() && td.keepParticle)
-        {
-            if (this->face() > -1)
-            {
-                if (isType<processorPolyPatch>(pbMesh[this->patch(this->face())]))
-                {
-                    td.switchProcessor = true;
-                }
-            }
-        }
-    }
-
-    return td.keepParticle;
-}
-
 
 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 

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

@@ -26,16 +26,14 @@ Class
     Foam::ReactingParcel
 
 Description
-    Reacting parcel class with one-way coupling with the continuous
+    Reactinf parcel class with one/two-way coupling with the continuous
-    phase.
+    phase. Includes thermo parcel sub-models, plus:
-
-    Includes thermo parcel sub-models, plus:
     - combustion
 
 SourceFiles
-    reactingParcelI.H
+    ReactingParcelI.H
     ReactingParcel.C
-    reactingParcelIO.C
+    ReactingParcelIO.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -75,96 +73,9 @@ class ReactingParcel
     public ThermoParcel<ParcelType>
 {
 
-    // Private data
-
-        //- Initial particle mass
-        scalar mass0_;
-
-        //- Mass fractions of mixture
-        scalarField YMixture_;
-
-        //- Mass fractions of gases
-        scalarField YGas_;
-
-        //- Mass fractions of liquids
-        scalarField YLiquid_;
-
-        //- Mass fractions of solids
-        scalarField YSolid_;
-
-        //- Flag to say that the particle is allowed to combust
-        //  Only true one ALL volatiles have been evolved from the particle
-        bool canCombust_;
-
-
-protected:
-
-    // Protected member functions
-
-        template<class TrackingData>
-        void calcCoupled
-        (
-            TrackingData& td,
-            const label celli,
-            const scalar dt,
-            const scalar rhoc,
-            vector& Uc,
-            const scalar muc,
-            const scalar Tc,
-            const scalar Cpc,
-            const scalar pc
-        );
-
-        template<class TrackingData>
-        void calcUncoupled
-        (
-            TrackingData& td,
-            const label celli,
-            const scalar dt,
-            const scalar rhoc,
-            vector& Uc,
-            const scalar muc,
-            const scalar Tc,
-            const scalar Cpc,
-            const scalar pc
-        );
-
-        //- Calculate mass transfer
-        template<class TrackingData>
-        void calcMassTransfer
-        (
-            TrackingData& td,
-            const scalar dt,
-            const scalar T0,
-            const scalar T1,
-            scalarList& dMassMT
-        );
-
-        //- Calculate surface reactions
-        template<class TrackingData>
-        void calcSurfaceReactions
-        (
-            TrackingData& td,
-            const scalar dt,
-            const label celli,
-            const scalar rhoc,
-            const scalar Tc,
-            const scalar T0,
-            const scalar T1,
-            scalar& dMassMTSR,
-            scalarList& dMassMT
-        );
-
-
 public:
 
-    //- Runtime type information
+    //- Class to hold reacting particle constant properties
-    TypeName("ReactingParcel");
-
-    friend class Cloud<ParcelType>;
-
-
-    //- Class to hold particle constant properties
     class constantProperties
     :
         public ThermoParcel<ParcelType>::constantProperties
@@ -196,7 +107,7 @@ public:
     };
 
 
-    //- Class used to pass tracking data to the trackToFace function
+    //- Class used to pass reacting tracking data to the trackToFace function
     class trackData
     :
         public ThermoParcel<ParcelType>::trackData
@@ -212,7 +123,7 @@ public:
 
             //- Interpolators for continuous phase fields
 
-                const interpolationCellPoint<scalar>& pInterp_;
+                const interpolation<scalar>& pInterp_;
 
     public:
 
@@ -222,12 +133,12 @@ public:
             (
                 ReactingCloud<ParcelType>& cloud,
                 const constantProperties& constProps,
-                const interpolationCellPoint<scalar>& rhoInterp,
+                const interpolation<scalar>& rhoInterp,
-                const interpolationCellPoint<vector>& UInterp,
+                const interpolation<vector>& UInterp,
-                const interpolationCellPoint<scalar>& muInterp,
+                const interpolation<scalar>& muInterp,
-                const interpolationCellPoint<scalar>& TInterp,
+                const interpolation<scalar>& TInterp,
-                const interpolationCellPoint<scalar>& CpInterp,
+                const interpolation<scalar>& CpInterp,
-                const interpolationCellPoint<scalar>& pInterp,
+                const interpolation<scalar>& pInterp,
                 const vector& g
             );
 
@@ -238,10 +149,96 @@ public:
 
             inline const constantProperties& constProps() const;
 
-            inline const interpolationCellPoint<scalar>& pInterp() const;
+            inline const interpolation<scalar>& pInterp() const;
     };
 
 
+protected:
+
+    // Protected data
+
+        // Parcel properties
+
+            //- Initial particle mass
+            scalar mass0_;
+
+            //- Mass fractions of mixture
+            scalarField YMixture_;
+
+            //- Mass fractions of gases
+            scalarField YGas_;
+
+            //- Mass fractions of liquids
+            scalarField YLiquid_;
+
+            //- Mass fractions of solids
+            scalarField YSolid_;
+
+            //- Flag to say that the particle is allowed to combust
+            //  Only true after volatile content falls below threshold value
+            bool canCombust_;
+
+
+        // Cell-based quantities
+
+            //- Pressure
+            scalar pc_;
+
+
+    // Protected member functions
+
+        virtual void updateCellQuantities
+        (
+            trackData& td,
+            const scalar dt,
+            const label celli
+        );
+
+        virtual void calcCoupled
+        (
+            trackData& td,
+            const scalar dt,
+            const label celli
+        );
+
+        virtual void calcUncoupled
+        (
+            trackData& td,
+            const scalar dt,
+            const label celli
+        );
+
+        //- Calculate mass transfer
+        virtual void calcMassTransfer
+        (
+            trackData& td,
+            const scalar dt,
+            const scalar T0,
+            const scalar T1,
+            scalarList& dMassMT
+        );
+
+        //- Calculate surface reactions
+        virtual void calcSurfaceReactions
+        (
+            trackData& td,
+            const scalar dt,
+            const label celli,
+            const scalar T0,
+            const scalar T1,
+            scalar& dMassMTSR,
+            scalarList& dMassMT
+        );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("ReactingParcel");
+
+    friend class Cloud<ParcelType>;
+
+
     // Constructors
 
         //- Construct from components
@@ -301,16 +298,6 @@ public:
             inline scalar& mass0();
 
 
-        // Tracking
-
-            //- Move parcel
-            template<class TrackingData>
-            bool move
-            (
-                TrackingData&
-            );
-
-
         // I-O
 
             static void readFields

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

@@ -43,12 +43,12 @@ inline Foam::ReactingParcel<ParcelType>::trackData::trackData
 (
     ReactingCloud<ParcelType>& cloud,
     const constantProperties& constProps,
-    const interpolationCellPoint<scalar>& rhoInterp,
+    const interpolation<scalar>& rhoInterp,
-    const interpolationCellPoint<vector>& UInterp,
+    const interpolation<vector>& UInterp,
-    const interpolationCellPoint<scalar>& muInterp,
+    const interpolation<scalar>& muInterp,
-    const interpolationCellPoint<scalar>& TInterp,
+    const interpolation<scalar>& TInterp,
-    const interpolationCellPoint<scalar>& CpInterp,
+    const interpolation<scalar>& CpInterp,
-    const interpolationCellPoint<scalar>& pInterp,
+    const interpolation<scalar>& pInterp,
     const vector& g
 )
 :
@@ -101,7 +101,8 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
     YMixture_(YMixture0),
     YGas_(YGas0),
     YLiquid_(YLiquid0),
-    YSolid_(YSolid0)
+    YSolid_(YSolid0),
+    pc_(0.0)
 {
     // Set initial parcel mass
     mass0_ = this->mass();
@@ -145,7 +146,7 @@ Foam::ReactingParcel<ParcelType>::trackData::constProps() const
 
 
 template<class ParcelType>
-inline const Foam::interpolationCellPoint<Foam::scalar>&
+inline const Foam::interpolation<Foam::scalar>&
 Foam::ReactingParcel<ParcelType>::trackData::pInterp() const
 {
     return pInterp_;

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

@@ -42,7 +42,8 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
     YMixture_(0),
     YGas_(0),
     YLiquid_(0),
-    YSolid_(0)
+    YSolid_(0),
+    pc_(0.0)
 {
     if (readFields)
     {

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

@@ -29,42 +29,61 @@ License
 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
 
 template<class ParcelType>
-template<class TrackingData>
+void Foam::ThermoParcel<ParcelType>::updateCellQuantities
+(
+    trackData& td,
+    const scalar dt,
+    const label celli
+)
+{
+    KinematicParcel<ParcelType>::updateCellQuantities(td, dt, celli);
+
+    Tc_ = td.TInterp().interpolate(this->position(), celli);
+    cpc_ = td.cpInterp().interpolate(this->position(), celli);
+}
+
+
+template<class ParcelType>
 void Foam::ThermoParcel<ParcelType>::calcCoupled
 (
-    TrackingData& td,
+    trackData& td,
-    const label celli,
     const scalar dt,
-    const scalar rhoc,
+    const label celli
-    vector& Uc,
-    const scalar muc,
-    const scalar Tc,
-    const scalar cpc
 )
 {
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Define local properties at beginning of timestep
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    const vector U0 = this->U();
+    const vector U0 = this->U_;
-    const scalar T0 = this->T();
     const scalar mass0 = this->mass();
-    const scalar cp0 = this->cp();
+    const scalar np0 = this->nParticle_;
-    const scalar np0 = this->nParticle();
+//    const scalar T0 = T_;
+//    const scalar cp0 = cp_;
+
+
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Initialise transfer terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Momentum transfer from the particle to the carrier phase
+    vector dUTrans = vector::zero;
+
+    // Enthalpy transfer from the particle to the carrier phase
+    scalar dhTrans = 0.0;
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate velocity - update U
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar Cud = 0.0;
-    const vector U1 = calcVelocity(td, dt, rhoc, Uc, muc, Cud);
+    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate heat transfer - update T
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar htc = 0.0;
-    const scalar T1 =
+    const scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
-        calcHeatTransfer(td, celli, dt, rhoc, Uc, muc, Tc, cpc, htc);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -72,13 +91,13 @@ void Foam::ThermoParcel<ParcelType>::calcCoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~
 
     // Update momentum transfer
-    td.cloud().UTrans()[celli] += np0*(mass0*(U0 - U1));
+    td.cloud().UTrans()[celli] += np0*dUTrans;
 
     // Accumulate coefficient to be applied in carrier phase momentum coupling
     td.cloud().UCoeff()[celli] += np0*mass0*Cud;
 
     // Update enthalpy transfer
-    td.cloud().hTrans()[celli] += np0*mass0*cp0*(T0 - T1);
+    td.cloud().hTrans()[celli] += np0*dhTrans;
 
     // Accumulate coefficient to be applied in carrier phase enthalpy coupling
     td.cloud().hCoeff()[celli] += np0*htc*this->areaS();
@@ -93,74 +112,76 @@ void Foam::ThermoParcel<ParcelType>::calcCoupled
 
 
 template<class ParcelType>
-template<class TrackingData>
 void Foam::ThermoParcel<ParcelType>::calcUncoupled
 (
-    TrackingData& td,
+    trackData& td,
-    const label celli,
     const scalar dt,
-    const scalar rhoc,
+    const label celli
-    vector& Uc,
-    const scalar muc,
-    const scalar Tc,
-    const scalar cpc
 )
 {
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Initialise transfer terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Momentum transfer from the particle to the carrier phase
+    vector dUTrans = vector::zero;
+
+    // Enthalpy transfer from the particle to the carrier phase
+    scalar dhTrans = 0.0;
+
+
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate velocity - update U
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar Cud = 0.0;
-    this->U() = calcVelocity(td, dt, rhoc, Uc, muc, Cud);
+    this->U_ = calcVelocity(td, dt, Cud, dUTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate heat transfer - update T
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar htc = 0.0;
-    this->T() = calcHeatTransfer(td, celli, dt, rhoc, Uc, muc, Tc, cpc, htc);
+    T_ = calcHeatTransfer(td, dt, celli, htc, dhTrans);
 }
 
 
 template<class ParcelType>
-template<class TrackingData>
 Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 (
-    TrackingData& td,
+    trackData& td,
-    const label celli,
     const scalar dt,
-    const scalar rhoc,
+    const label celli,
-    const vector& Uc,
+    scalar& htc,
-    const scalar muc,
+    scalar& dhTrans
-    const scalar Tc,
-    const scalar cpc,
-    scalar& htc
 )
 {
     if (!td.cloud().heatTransfer().active())
     {
+        htc = 0.0;
+        dhTrans = 0.0;
         return T_;
     }
 
     // Calc heat transfer coefficient
     htc = td.cloud().heatTransfer().h
     (
-        this->d(),
+        this->d_,
-        this->Ur(),
+        this->U_ - this->Uc_,
-        rhoc,
+        this->rhoc_,
-        this->rho(),
+        this->rho_,
-        cpc,
+        cpc_,
         cp_,
-        muc
+        this->muc_
     );
 
     // Determine ap and bp coefficients
-    scalar ap = Tc;
+    scalar ap = Tc_;
     scalar bp = htc;
     if (td.cloud().radiation())
     {
         // Carrier phase incident radiation field
-        // Currently the G field is not interpolated to the parcel position
+        // - The G field is not interpolated to the parcel position
-        // - instead, the cell centre value is applied directly
+        //   Instead, the cell centre value is applied directly
         const scalarField& G = td.cloud().mesh().objectRegistry
             ::lookupObject<volScalarField>("G");
 
@@ -168,21 +189,21 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
         const scalar sigma = radiation::sigmaSB.value();
         const scalar epsilon = td.constProps().epsilon0();
         const scalar epsilonSigmaT3 = epsilon*sigma*pow3(T_);
-        ap = (htc*Tc + 0.25*epsilon*G[celli])/(htc + epsilonSigmaT3);
+        ap = (htc*Tc_ + 0.25*epsilon*G[celli])/(htc + epsilonSigmaT3);
         bp += epsilonSigmaT3;
     }
-    bp *= 6.0/(this->rho()*this->d()*cp_);
+    bp *= 6.0/(this->rho_*this->d_*cp_);
 
 
     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Set new particle temperature
     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Analytical
+    scalar Tnew = td.cloud().TIntegrator().integrate(T_, dt, ap, bp);
-    const scalar Tnew = ap + (T_ - ap)*exp(-bp*dt);
 
-    // Euler-implicit
+    Info<< "T_, Tnew = " << T_ << ", " << Tnew << endl;
-//    const scalar Tnew = (T_ + dt*ap*bp)/(1.0 + dt*bp);
+
+    dhTrans = -this->mass()*cp_*(Tnew - T_);
 
     return Tnew;
 }
@@ -190,110 +211,6 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class ParcelType>
-template<class TrackingData>
-bool Foam::ThermoParcel<ParcelType>::move
-(
-    TrackingData& td
-)
-{
-    td.switchProcessor = false;
-    td.keepParticle = true;
-
-    const polyMesh& mesh = td.cloud().pMesh();
-    const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
-
-    const scalar deltaT = mesh.time().deltaT().value();
-    scalar tEnd = (1.0 - this->stepFraction())*deltaT;
-    const scalar dtMax = tEnd;
-
-    while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
-    {
-        // Set the Lagrangian time-step
-        scalar dt = min(dtMax, tEnd);
-
-        // Remember which cell the parcel is in
-        // since this will change if a face is hit
-        label celli = this->cell();
-
-        dt *= trackToFace(this->position() + dt*this->U(), td);
-
-        tEnd -= dt;
-        this->stepFraction() = 1.0 - tEnd/deltaT;
-
-        cellPointWeight cpw
-        (
-            mesh,
-            this->position(),
-            celli,
-            this->faceInterpolation()
-        );
-        scalar rhoc = td.rhoInterp().interpolate(cpw);
-        vector Uc = td.UInterp().interpolate(cpw);
-        scalar muc = td.muInterp().interpolate(cpw);
-        scalar Tc = td.TInterp().interpolate(cpw);
-        scalar cpc = td.cpInterp().interpolate(cpw);
-
-        Uc = td.cloud().dispersion().update
-        (
-            dt,
-            celli,
-            this->U(),
-            Uc,
-            this->UTurb(),
-            this->tTurb()
-        );
-
-        if (td.cloud().coupled())
-        {
-            calcCoupled(td, celli, dt, rhoc, Uc, muc, Tc, cpc);
-        }
-        else
-        {
-            calcUncoupled(td, celli, dt, rhoc, Uc, muc, Tc, cpc);
-        }
-
-        if (this->onBoundary() && td.keepParticle)
-        {
-            if (this->face() > -1)
-            {
-                if
-                (
-                    isType<processorPolyPatch>
-                        (pbMesh[this->patch(this->face())])
-                )
-                {
-                    td.switchProcessor = true;
-                }
-            }
-        }
-    }
-
-    return td.keepParticle;
-}
-
-
-template<class ParcelType>
-template<class TrackingData>
-void Foam::ThermoParcel<ParcelType>::hitWallPatch
-(
-    const wallPolyPatch& wpp,
-    TrackingData& td
-)
-{
-    td.cloud().wallInteraction().correct(wpp, this->face(), this->U());
-}
-
-
-
-template<class ParcelType>
-void Foam::ThermoParcel<ParcelType>::hitWallPatch
-(
-    const wallPolyPatch& wpp,
-    int&
-)
-{}
-
 
 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 

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

@@ -26,16 +26,14 @@ Class
     Foam::ThermoParcel
 
 Description
-    Thermodynamic parcel class with one-way coupling with the continuous
+    Thermodynamic parcel class with one/two-way coupling with the continuous
-    phase.
+    phase. Includes Kinematic parcel sub-models, plus:
-
-    Includes Kinematic parcel sub-models, plus:
     - heat transfer
 
 SourceFiles
-    thermoParcelI.H
+    ThermoParcelI.H
     ThermoParcel.C
-    thermoParcelIO.C
+    ThermoParcelIO.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -75,72 +73,9 @@ class ThermoParcel
     public KinematicParcel<ParcelType>
 {
 
-    // Private data
-
-        //- Parcel properties
-
-            //- Temperature [K]
-            scalar T_;
-
-            //- Specific heat capacity [J/(kg.K)]
-            scalar cp_;
-
-
-protected:
-
-    // Protected member functions
-
-        template<class TrackingData>
-        void calcCoupled
-        (
-            TrackingData& td,
-            const label celli,
-            const scalar dt,
-            const scalar rhoc,
-            vector& Uc,
-            const scalar muc,
-            const scalar Tc,
-            const scalar cpc
-        );
-
-        template<class TrackingData>
-        void calcUncoupled
-        (
-            TrackingData& td,
-            const label celli,
-            const scalar dt,
-            const scalar rhoc,
-            vector& Uc,
-            const scalar muc,
-            const scalar Tc,
-            const scalar cpc
-        );
-
-        //- Calculate new particle temperature
-        template<class TrackingData>
-        scalar calcHeatTransfer
-        (
-            TrackingData& td,
-            const label celli,
-            const scalar dt,
-            const scalar rhoc,
-            const vector& Uc,
-            const scalar muc,
-            const scalar Tc,
-            const scalar cpc,
-            scalar& htc
-        );
-
-
 public:
 
-    //- Runtime type information
+    //- Class to hold thermo particle constant properties
-    TypeName("ThermoParcel");
-
-    friend class Cloud<ParcelType>;
-
-
-    //- Class to hold particle constant properties
     class constantProperties
     :
         public KinematicParcel<ParcelType>::constantProperties
@@ -180,7 +115,7 @@ public:
     };
 
 
-    //- Class used to pass tracking data to the trackToFace function
+    //- Class used to pass thermo tracking data to the trackToFace function
     class trackData
     :
         public KinematicParcel<ParcelType>::trackData
@@ -196,8 +131,8 @@ public:
 
             //- Interpolators for continuous phase fields
 
-                const interpolationCellPoint<scalar>& TInterp_;
+                const interpolation<scalar>& TInterp_;
-                const interpolationCellPoint<scalar>& cpInterp_;
+                const interpolation<scalar>& cpInterp_;
 
 
     public:
@@ -208,11 +143,11 @@ public:
             (
                 ThermoCloud<ParcelType>& cloud,
                 const constantProperties& constProps,
-                const interpolationCellPoint<scalar>& rhoInterp,
+                const interpolation<scalar>& rhoInterp,
-                const interpolationCellPoint<vector>& UInterp,
+                const interpolation<vector>& UInterp,
-                const interpolationCellPoint<scalar>& muInterp,
+                const interpolation<scalar>& muInterp,
-                const interpolationCellPoint<scalar>& TInterp,
+                const interpolation<scalar>& TInterp,
-                const interpolationCellPoint<scalar>& cpInterp,
+                const interpolation<scalar>& cpInterp,
                 const vector& g
             );
 
@@ -223,12 +158,76 @@ public:
 
             inline const constantProperties& constProps() const;
 
-            inline const interpolationCellPoint<scalar>& TInterp() const;
+            inline const interpolation<scalar>& TInterp() const;
 
-            inline const interpolationCellPoint<scalar>& cpInterp() const;
+            inline const interpolation<scalar>& cpInterp() const;
     };
 
 
+protected:
+
+    // Protected data
+
+        // Parcel properties
+
+            //- Temperature [K]
+            scalar T_;
+
+            //- Specific heat capacity [J/(kg.K)]
+            scalar cp_;
+
+
+        // Call-based quantities
+
+            //- Temperature
+            scalar Tc_;
+
+            //- Specific heat capacity
+            scalar cpc_;
+
+
+    // Protected member functions
+
+        virtual void updateCellQuantities
+        (
+            trackData& td,
+            const scalar dt,
+            const label celli
+        );
+
+        virtual void calcCoupled
+        (
+            trackData& td,
+            const scalar dt,
+            const label celli
+        );
+
+        virtual void calcUncoupled
+        (
+            trackData& td,
+            const scalar dt,
+            const label celli
+        );
+
+        //- Calculate new particle temperature
+        virtual scalar calcHeatTransfer
+        (
+            trackData& td,
+            const scalar dt,
+            const label celli,
+            scalar& htc,
+            scalar& dhTrans
+        );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("ThermoParcel");
+
+    friend class Cloud<ParcelType>;
+
+
     // Constructors
 
         //- Construct from components
@@ -271,30 +270,6 @@ public:
             inline const scalar cp() const;
             inline scalar& cp();
 
-            //- Move parcel
-            template<class TrackingData>
-            bool move
-            (
-                 TrackingData&
-            );
-
-
-        //- Patch interaction
-
-            //- Overridable function to handle the particle hitting a wallPatch
-            template<class TrackingData>
-            void hitWallPatch
-            (
-                 const wallPolyPatch&,
-                 TrackingData& td
-            );
-
-            void hitWallPatch
-            (
-                 const wallPolyPatch&,
-                 int&
-            );
-
 
         // I-O
 

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

@@ -45,11 +45,11 @@ inline Foam::ThermoParcel<ParcelType>::trackData::trackData
 (
     ThermoCloud<ParcelType>& cloud,
     const constantProperties& constProps,
-    const interpolationCellPoint<scalar>& rhoInterp,
+    const interpolation<scalar>& rhoInterp,
-    const interpolationCellPoint<vector>& UInterp,
+    const interpolation<vector>& UInterp,
-    const interpolationCellPoint<scalar>& muInterp,
+    const interpolation<scalar>& muInterp,
-    const interpolationCellPoint<scalar>& TInterp,
+    const interpolation<scalar>& TInterp,
-    const interpolationCellPoint<scalar>& cpInterp,
+    const interpolation<scalar>& cpInterp,
     const vector& g
 )
 :
@@ -94,7 +94,9 @@ inline Foam::ThermoParcel<ParcelType>::ThermoParcel
         constProps
     ),
     T_(constProps.T0()),
-    cp_(constProps.cp0())
+    cp_(constProps.cp0()),
+    Tc_(0.0),
+    cpc_(0.0)
 {}
 
 
@@ -151,7 +153,7 @@ Foam::ThermoParcel<ParcelType>::trackData::constProps() const
 
 
 template<class ParcelType>
-inline const Foam::interpolationCellPoint<Foam::scalar>&
+inline const Foam::interpolation<Foam::scalar>&
 Foam::ThermoParcel<ParcelType>::trackData::TInterp() const
 {
     return TInterp_;
@@ -159,7 +161,7 @@ Foam::ThermoParcel<ParcelType>::trackData::TInterp() const
 
 
 template<class ParcelType>
-inline const Foam::interpolationCellPoint<Foam::scalar>&
+inline const Foam::interpolation<Foam::scalar>&
 Foam::ThermoParcel<ParcelType>::trackData::cpInterp() const
 {
     return cpInterp_;

src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C

@@ -39,7 +39,9 @@ Foam::ThermoParcel<ParcelType>::ThermoParcel
 :
     KinematicParcel<ParcelType>(cloud, is, readFields),
     T_(0.0),
-    cp_(0.0)
+    cp_(0.0),
+    Tc_(0.0),
+    cpc_(0.0)
 {
     if (readFields)
     {

src/lagrangian/intermediate/parcels/derived/kinematicParcel/kinematicParcel.H

@@ -26,7 +26,7 @@ Class
     Foam::kinematicParcel
 
 Description
-    Foam::kinematicParcel
+
 
 SourceFiles
     kinematicParcel.C

src/lagrangian/intermediate/parcels/derived/reactingParcel/reactingParcel.H

@@ -26,7 +26,7 @@ Class
     Foam::reactingParcel
 
 Description
-    Foam::reactingParcel
+
 
 SourceFiles
     reactingParcel.C

src/lagrangian/intermediate/parcels/derived/thermoParcel/thermoParcel.H

@@ -26,7 +26,7 @@ Class
     Foam::thermoParcel
 
 Description
-    Foam::thermoParcel
+
 
 SourceFiles
     thermoParcel.C

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

@@ -26,7 +26,6 @@ Class
     Foam::DispersionModel
 
 Description
-    Templated dispersion model class
 
 \*---------------------------------------------------------------------------*/
 

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/DispersionRASModel/DispersionRASModel.H

@@ -26,10 +26,6 @@ Class
     Foam::DispersionRASModel
 
 Description
-    Templated dispersion RAS model
-
-SourceFiles
-    DispersionRASModel.C
 
 \*---------------------------------------------------------------------------*/
 

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.H

@@ -27,7 +27,7 @@ Class
 
 Description
     The velocity is perturbed in the direction of -grad(k), with a
-    Gaussian random number distribution with variance sigma,
+    Gaussian random number distribution with variance sigma.
     where sigma is defined below
 
 \*---------------------------------------------------------------------------*/

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C

@@ -22,6 +22,8 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
+Description
+
 \*---------------------------------------------------------------------------*/
 
 #include "StochasticDispersionRAS.H"

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.H

@@ -26,10 +26,9 @@ Class
     Foam::StochasticDispersionRAS
 
 Description
-    The velocity is perturbed in random direction.
+    The velocity is perturbed in random direction, with a
- 
+    Gaussian random number distribution with variance sigma.
-    A Gaussian random number distribution is used with variance sigma,
+    where sigma is defined below 
-    where sigma is defined below
 
 \*---------------------------------------------------------------------------*/
 

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

@@ -25,6 +25,7 @@ License
 Class
     Foam::DragModel
 
+
 Description
     Templated drag model class
 
@@ -121,10 +122,7 @@ public:
         virtual bool active() const = 0;
 
         //- Return drag coefficient
-        virtual scalar Cd
+        virtual scalar Cd(const scalar Re) const = 0;
-        (
-            const scalar Re
-        ) const = 0;
 
         //- Return linearised coefficient for velocity equation
         //  Drag force per unit particle mass = Cu(U - Up)

src/lagrangian/intermediate/submodels/Kinematic/DragModel/NoDrag/NoDrag.H

@@ -74,10 +74,7 @@ public:
 
         bool active() const;
 
-        scalar Cd
+        scalar Cd(const scalar) const;
-        (
-            const scalar
-        ) const;
 };
 
 

src/lagrangian/intermediate/submodels/Kinematic/DragModel/SphereDrag/SphereDrag.H

@@ -74,10 +74,7 @@ public:
 
         bool active() const;
 
-        scalar Cd
+        scalar Cd(const scalar Re) const;
-        (
-            const scalar Re
-        ) const;
 };
 
 

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

@@ -25,6 +25,7 @@ License
 Class
     Foam::InjectionModel
 
+
 Description
     Templated injection model class
 

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

@@ -46,7 +46,7 @@ Foam::ManualInjection<CloudType>::ManualInjection
         IOobject
         (
             positionsFile_,
-            owner.runTime().constant(),
+            owner.db().time().constant(),
             owner.mesh(),
             IOobject::MUST_READ,
             IOobject::NO_WRITE

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

@@ -27,12 +27,10 @@ Class
 
 Description
     Manual injection
-
-    For manual injection,
     - User specifies
-      -# Total mass to inject
+      - Total mass to inject
-      -# Parcel positions in file \<positionsFile\>
+      - Parcel positions in file <positionsFile>
-      -# Initial parcel velocity
+      - Initial parcel velocity
     - Parcel diameters obtained by PDF model
     - All parcels introduced at the start of the calculation
 

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

@@ -25,6 +25,7 @@ License
 Class
     Foam::WallInteractionModel
 
+
 Description
     Templated wall interaction model class
 

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

@@ -25,9 +25,9 @@ License
 Class
     Foam::CompositionModel
 
-Description
-    Templated reacting parcel composition model class.
 
+Description
+    Templated reacting parcel composition model class
     Consists of gases (via thermo package), liquids and solids
 
 SourceFiles

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

@@ -28,7 +28,6 @@ Class
 
 Description
     Templated parcel single mixture fraction class
-
     - Each phase sums to a mass fraction of 1
 
 SourceFiles

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

@@ -25,6 +25,7 @@ License
 Class
     Foam::MassTransferModel
 
+
 Description
     Templated mass transfer model class
 

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

@@ -25,6 +25,7 @@ License
 Class
     Foam::SurfaceReactionModel
 
+
 Description
     Templated surface reaction model class
 

src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/HeatTransferModel/HeatTransferModel.H

@@ -25,6 +25,7 @@ License
 Class
     Foam::HeatTransferModel
 
+
 Description
     Templated heat transfer model class
 

src/lagrangian/intermediate/submodels/addOns/radiation/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C

@@ -72,7 +72,7 @@ Foam::radiation::cloudAbsorptionEmission::~cloudAbsorptionEmission()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::cloudAbsorptionEmission::a() const
+Foam::radiation::cloudAbsorptionEmission::aDisp() const
 {
     tmp<volScalarField> ta
     (
@@ -107,7 +107,7 @@ Foam::radiation::cloudAbsorptionEmission::a() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::cloudAbsorptionEmission::e() const
+Foam::radiation::cloudAbsorptionEmission::eDisp() const
 {
     tmp<volScalarField> te
     (
@@ -132,7 +132,7 @@ Foam::radiation::cloudAbsorptionEmission::e() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::cloudAbsorptionEmission::E() const
+Foam::radiation::cloudAbsorptionEmission::EDisp() const
 {
     tmp<volScalarField> tE
     (

src/lagrangian/intermediate/submodels/addOns/radiation/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.H

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::radiation::cloudAbsorptionEmission
+    Foam::cloudAbsorptionEmission
 
 Description
     Retrieves absorption/emission data from a cloud object
@@ -86,18 +86,24 @@ public:
 
     // Member Operators
 
-        // Edit
-
         // Access
 
-            //- Return absorption coefficient
+            // Absorption coefficient
-            tmp<volScalarField> a() const;
 
-            //- Return emission coefficient
+                //- Absorption coefficient for dispersed phase
-            tmp<volScalarField> e() const;
+                tmp<volScalarField> aDisp() const;
 
-            //- Return emission contribution
+
-            tmp<volScalarField> E() const;
+            // Emission coefficient
+
+                //- Emission coefficient for dispersed phase
+                tmp<volScalarField> eDisp() const;
+
+
+            // Emission contribution
+
+                //- Return emission contribution for dispersed phase
+                tmp<volScalarField> EDisp() const;
 };
 
 

src/lagrangian/intermediate/submodels/addOns/radiation/scatter/cloudScatter/cloudScatter.H

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::radiation::cloudScatter
+    Foam::scatterModel
 
 Description
     Cloud radiation scatter model

src/thermophysicalModels/radiation/Make/files

@@ -24,8 +24,8 @@ submodels/absorptionEmissionModel/binaryAbsorptionEmission/binaryAbsorptionEmiss
 
 
 /* Boundary conditions */
-boundaryConditions/MarshakRadiation/MarshakRadiationMixedFvPatchScalarField.C
+derivedFvPatchFields/MarshakRadiation/MarshakRadiationMixedFvPatchScalarField.C
-boundaryConditions/MarshakRadiationFixedT/MarshakRadiationFixedTMixedFvPatchScalarField.C
+derivedFvPatchFields/MarshakRadiationFixedT/MarshakRadiationFixedTMixedFvPatchScalarField.C
 
 
 LIB = $(FOAM_LIBBIN)/libradiation

src/thermophysicalModels/radiation/radiationModel/P1/P1.C

@@ -53,13 +53,9 @@ namespace Foam
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
-Foam::radiation::P1::P1
+Foam::radiation::P1::P1(const volScalarField& T)
-(
-    const fvMesh& mesh,
-    const basicThermo& thermo
-)
 :
-    radiationModel(typeName, mesh, thermo),
+    radiationModel(typeName, T),
     G_
     (
         IOobject
@@ -143,12 +139,10 @@ void Foam::radiation::P1::correct()
     {
         return;
     }
-
     a_ = absorptionEmission_->a();
     e_ = absorptionEmission_->e();
     E_ = absorptionEmission_->E();
     const volScalarField sigmaEff = scatter_->sigmaEff();
-    const volScalarField& T = thermo_.T();
 
     // Construct diffusion
     const volScalarField gamma
@@ -170,7 +164,7 @@ void Foam::radiation::P1::correct()
         fvm::laplacian(gamma, G_)
       - fvm::Sp(a_, G_)
      ==
-      - 4.0*(e_*radiation::sigmaSB*pow4(T) + mathematicalConstant::pi*E_)
+      - 4.0*(e_*radiation::sigmaSB*pow4(T_) + E_)
     );
 }
 
@@ -187,17 +181,26 @@ Foam::tmp<Foam::volScalarField> Foam::radiation::P1::Rp() const
                 mesh_.time().timeName(),
                 mesh_,
                 IOobject::NO_READ,
-                IOobject::NO_WRITE
+                IOobject::NO_WRITE,
+                false
             ),
-            4.0*e_*radiation::sigmaSB
+            4.0*absorptionEmission_->eCont()*radiation::sigmaSB
         )
     );
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::radiation::P1::Ru() const
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::radiation::P1::Ru() const
 {
-    return a_*G_ - 4.0*mathematicalConstant::pi*E_;
+    const DimensionedField<scalar, volMesh>& G =
+        G_.dimensionedInternalField();
+    const DimensionedField<scalar, volMesh> E =
+        absorptionEmission_->ECont()().dimensionedInternalField();
+    const DimensionedField<scalar, volMesh> a =
+        absorptionEmission_->aCont()().dimensionedInternalField();
+
+    return  a*G - 4.0*E;
 }
 
 

src/thermophysicalModels/radiation/radiationModel/P1/P1.H

@@ -93,7 +93,7 @@ public:
     // Constructors
 
         //- Construct from components
-        P1(const fvMesh& mesh, const basicThermo& thermo);
+        P1(const volScalarField& T);
 
 
     // Destructor
@@ -118,7 +118,7 @@ public:
             virtual tmp<volScalarField> Rp() const;
 
             //- Source term component (constant)
-            virtual tmp<volScalarField> Ru() const;
+            virtual tmp<DimensionedField<scalar, volMesh> > Ru() const;
 };
 
 

src/thermophysicalModels/radiation/radiationModel/noRadiation/noRadiation.C

@@ -51,13 +51,9 @@ namespace Foam
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
-Foam::radiation::noRadiation::noRadiation
+Foam::radiation::noRadiation::noRadiation(const volScalarField& T)
-(
-    const fvMesh& mesh,
-    const basicThermo& thermo
-)
 :
-    radiationModel(typeName, mesh, thermo)
+    radiationModel(typeName, T)
 {}
 
 
@@ -99,8 +95,7 @@ Foam::tmp<Foam::volScalarField> Foam::radiation::noRadiation::Rp() const
             dimensionedScalar
             (
                 "Rp",
-                radiation::sigmaSB.dimensions()/dimLength
+                radiation::sigmaSB.dimensions()/dimLength,
-                   *pow3(dimTemperature),
                 0.0
             )
         )
@@ -108,11 +103,12 @@ Foam::tmp<Foam::volScalarField> Foam::radiation::noRadiation::Rp() const
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::radiation::noRadiation::Ru() const
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::radiation::noRadiation::Ru() const
 {
-    return tmp<volScalarField>
+    return tmp<DimensionedField<scalar, volMesh> >
     (
-        new volScalarField
+        new DimensionedField<scalar, volMesh>
         (
             IOobject
             (

src/thermophysicalModels/radiation/radiationModel/noRadiation/noRadiation.H

@@ -73,7 +73,7 @@ public:
     // Constructors
 
         //- Construct from components
-        noRadiation(const fvMesh& mesh, const basicThermo& thermo);
+        noRadiation(const volScalarField& T);
 
 
     // Destructor
@@ -95,7 +95,7 @@ public:
             tmp<volScalarField> Rp() const;
 
             //- Source term component (constant)
-            tmp<volScalarField> Ru() const;
+            tmp<DimensionedField<scalar, volMesh> > Ru() const;
 };
 
 

src/thermophysicalModels/radiation/radiationModel/radiationModel/newRadiationModel.C

@@ -38,8 +38,7 @@ namespace radiation
 
 autoPtr<radiationModel> radiationModel::New
 (
-    const fvMesh& mesh,
+    const volScalarField& T
-    const basicThermo& thermo
 )
 {
     word radiationModelTypeName;
@@ -53,8 +52,8 @@ autoPtr<radiationModel> radiationModel::New
             IOobject
             (
                 "radiationProperties",
-                mesh.time().constant(),
+                T.mesh().time().constant(),
-                mesh.db(),
+                T.mesh().db(),
                 IOobject::MUST_READ,
                 IOobject::NO_WRITE
             )
@@ -73,7 +72,7 @@ autoPtr<radiationModel> radiationModel::New
     {
         FatalErrorIn
         (
-            "radiationModel::New(const fvMesh&, const basicThermo&)"
+            "radiationModel::New(const volScalarField&)"
         )   << "Unknown radiationModel type " << radiationModelTypeName
             << nl << nl
             << "Valid radiationModel types are :" << nl
@@ -81,7 +80,7 @@ autoPtr<radiationModel> radiationModel::New
             << exit(FatalError);
     }
 
-    return autoPtr<radiationModel>(cstrIter()(mesh, thermo));
+    return autoPtr<radiationModel>(cstrIter()(T));
 }
 
 

src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.C

@@ -46,8 +46,7 @@ namespace Foam
 Foam::radiation::radiationModel::radiationModel
 (
     const word& type,
-    const fvMesh& mesh,
+    const volScalarField& T
-    const basicThermo& thermo
 )
 :
     IOdictionary
@@ -55,18 +54,18 @@ Foam::radiation::radiationModel::radiationModel
         IOobject
         (
             "radiationProperties",
-            mesh.time().constant(),
+            T.mesh().time().constant(),
-            mesh.db(),
+            T.mesh().db(),
             IOobject::MUST_READ,
             IOobject::NO_WRITE
         )
     ),
-    thermo_(thermo),
+    T_(T),
-    mesh_(mesh),
+    mesh_(T.mesh()),
     radiation_(lookup("radiation")),
     radiationModelCoeffs_(subDict(type + "Coeffs")),
     absorptionEmission_(absorptionEmissionModel::New(*this, mesh_)),
-    scatter_(scatterModel::New(*this, mesh))
+    scatter_(scatterModel::New(*this, mesh_))
 {}
 
 
@@ -96,18 +95,18 @@ bool Foam::radiation::radiationModel::read()
 
 Foam::tmp<Foam::fvScalarMatrix> Foam::radiation::radiationModel::Sh
 (
-    volScalarField& h
+    basicThermo& thermo
 ) const
 {
-    const volScalarField cp = thermo_.Cp();
+    volScalarField& h = thermo.h();
-    const volScalarField& T = thermo_.T();
+    const volScalarField cp = thermo.Cp();
-    const volScalarField T3 = pow3(T);
+    const volScalarField T3 = pow3(T_);
 
     return
     (
         Ru()
-      - fvm::Sp(Rp()*T3/cp, h)
+      - fvm::Sp(4.0*Rp()*T3/cp, h)
-      - Rp()*T3*(T - h/cp)
+      - Rp()*T3*(T_ - 4.0*h/cp)
     );
 }
 

src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.H

@@ -73,8 +73,8 @@ protected:
 
     // Protected data
 
-        //- Reference to the thermo package
+        //- Reference to the temperature field
-        const basicThermo& thermo_;
+        const volScalarField& T_;
 
         //- Reference to the mesh
         const fvMesh& mesh_;
@@ -120,10 +120,9 @@ public:
              radiationModel,
              dictionary,
              (
-                 const fvMesh& mesh,
+                 const volScalarField& T
-                 const basicThermo& thermo
              ),
-             (mesh, thermo)
+             (T)
          );
 
 
@@ -133,8 +132,7 @@ public:
         radiationModel
         (
             const word& type,
-            const fvMesh& mesh,
+            const volScalarField& T
-            const basicThermo& thermo
         );
 
 
@@ -143,8 +141,7 @@ public:
          //- Return a reference to the selected radiation model
          static autoPtr<radiationModel> New
          (
-             const fvMesh& mesh,
+             const volScalarField& T
-             const basicThermo& thermo
          );
 
 
@@ -170,12 +167,12 @@ public:
             virtual tmp<volScalarField> Rp() const = 0;
 
             //- Source term component (constant)
-            virtual tmp<volScalarField> Ru() const = 0;
+            virtual tmp<DimensionedField<scalar, volMesh> > Ru() const = 0;
 
-            //- Source term
+            //- Enthalpy source term
             virtual tmp<fvScalarMatrix> Sh
             (
-                volScalarField& h
+                basicThermo& thermo
             ) const;
 };
 

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/absorptionEmissionModel/absorptionEmissionModel.C

@@ -56,4 +56,165 @@ Foam::radiation::absorptionEmissionModel::~absorptionEmissionModel()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::absorptionEmissionModel::a() const
+{
+    return aDisp() + aCont();
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::absorptionEmissionModel::aCont() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "aCont",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimless/dimLength, 0.0)
+        )
+    );
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::absorptionEmissionModel::aDisp() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "aDisp",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimless/dimLength, 0.0)
+        )
+    );
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::absorptionEmissionModel::e() const
+{
+    return eDisp() + eCont();
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::absorptionEmissionModel::eCont() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "eCont",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimless/dimLength, 0.0)
+        )
+    );
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::absorptionEmissionModel::eDisp() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "eDisp",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimless/dimLength, 0.0)
+        )
+    );
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::absorptionEmissionModel::E() const
+{
+    return EDisp() + ECont();
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::absorptionEmissionModel::ECont() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "ECont",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimMass/dimLength/pow3(dimTime), 0.0)
+        )
+    );
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::absorptionEmissionModel::EDisp() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "EDisp",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimMass/dimLength/pow3(dimTime), 0.0)
+        )
+    );
+}
+
+
 // ************************************************************************* //

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/absorptionEmissionModel/absorptionEmissionModel.H

@@ -69,7 +69,7 @@ public:
     //- Runtime type information
     TypeName("absorptionEmissionModel");
 
-    // Declare runtime constructor selection table
+    //- Declare runtime constructor selection table
 
         declareRunTimeSelectionTable
         (
@@ -94,7 +94,7 @@ public:
         );
 
 
-    // Selector
+    //- Selector
 
         static autoPtr<absorptionEmissionModel> New
         (
@@ -103,7 +103,7 @@ public:
         );
 
 
-    // Destructor
+    //- Destructor
 
         virtual ~absorptionEmissionModel();
 
@@ -112,14 +112,40 @@ public:
 
         // Access
 
-            //- Return absorption coefficient
+            // Absorption coefficient
-            virtual tmp<volScalarField> a() const = 0;
 
-            //- Return emission coefficient
+                //- Absorption coefficient (net)
-            virtual tmp<volScalarField> e() const = 0;
+                virtual tmp<volScalarField> a() const;
 
-            //- Return emission contribution
+                //- Absorption coefficient for continuous phase
-            virtual tmp<volScalarField> E() const = 0;
+                virtual tmp<volScalarField> aCont() const;
+
+                //- Absorption coefficient for dispersed phase
+                virtual tmp<volScalarField> aDisp() const;
+
+
+            // Emission coefficient
+
+                //- Emission coefficient (net)
+                virtual tmp<volScalarField> e() const;
+
+                //- Return emission coefficient for continuous phase
+                virtual tmp<volScalarField> eCont() const;
+
+                //- Return emission coefficient for dispersed phase
+                virtual tmp<volScalarField> eDisp() const;
+
+
+            // Emission contribution
+
+                //- Emission contribution (net)
+                virtual tmp<volScalarField> E() const;
+
+                //- Emission contribution for continuous phase
+                virtual tmp<volScalarField> ECont() const;
+
+                //- Emission contribution for dispersed phase
+                virtual tmp<volScalarField> EDisp() const;
 };
 
 

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/binaryAbsorptionEmission/binaryAbsorptionEmission.C

@@ -74,21 +74,45 @@ Foam::radiation::binaryAbsorptionEmission::~binaryAbsorptionEmission()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volScalarField> Foam::radiation::binaryAbsorptionEmission::a() const
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::binaryAbsorptionEmission::aCont() const
 {
-    return model1_->a() + model2_->a();
+    return model1_->aCont() + model2_->aCont();
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::radiation::binaryAbsorptionEmission::e() const
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::binaryAbsorptionEmission::aDisp() const
 {
-    return model1_->e() + model2_->e();
+    return model1_->aDisp() + model2_->aDisp();
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::radiation::binaryAbsorptionEmission::E() const
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::binaryAbsorptionEmission::eCont() const
 {
-    return model1_->E() + model2_->E();
+    return model1_->eCont() + model2_->eCont();
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::binaryAbsorptionEmission::eDisp() const
+{
+    return model1_->eDisp() + model2_->eDisp();
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::binaryAbsorptionEmission::ECont() const
+{
+    return model1_->ECont() + model2_->ECont();
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::binaryAbsorptionEmission::EDisp() const
+{
+    return model1_->EDisp() + model2_->EDisp();
 }
 
 

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/binaryAbsorptionEmission/binaryAbsorptionEmission.H

@@ -89,18 +89,33 @@ public:
 
     // Member Operators
 
-        // Edit
-
         // Access
 
-            //- Return absorption coefficient
+            // Absorption coefficient
-            tmp<volScalarField> a() const;
 
-            //- Return emission coefficient
+                //- Absorption coefficient for continuous phase
-            tmp<volScalarField> e() const;
+                tmp<volScalarField> aCont() const;
 
-            //- Return emission contribution
+                //- Absorption coefficient for dispersed phase
-            tmp<volScalarField> E() const;
+                tmp<volScalarField> aDisp() const;
+
+
+            // Emission coefficient
+
+                //- Emission coefficient for continuous phase
+                tmp<volScalarField> eCont() const;
+
+                //- Emission coefficient for dispersed phase
+                tmp<volScalarField> eDisp() const;
+
+
+            // Emission contribution
+
+                //- Emission contribution for continuous phase
+                tmp<volScalarField> ECont() const;
+
+                //- Emission contribution for continuous phase
+                tmp<volScalarField> EDisp() const;
 };
 
 

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/constantAbsorptionEmission/constantAbsorptionEmission.C

@@ -70,7 +70,7 @@ Foam::radiation::constantAbsorptionEmission::~constantAbsorptionEmission()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::constantAbsorptionEmission::a() const
+Foam::radiation::constantAbsorptionEmission::aCont() const
 {
     tmp<volScalarField> ta
     (
@@ -95,7 +95,7 @@ Foam::radiation::constantAbsorptionEmission::a() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::constantAbsorptionEmission::e() const
+Foam::radiation::constantAbsorptionEmission::eCont() const
 {
     tmp<volScalarField> te
     (
@@ -119,7 +119,8 @@ Foam::radiation::constantAbsorptionEmission::e() const
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::radiation::constantAbsorptionEmission::E() const
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::constantAbsorptionEmission::ECont() const
 {
     tmp<volScalarField> tE
     (

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/constantAbsorptionEmission/constantAbsorptionEmission.H

@@ -26,7 +26,8 @@ Class
     Foam::radiation::constantAbsorptionEmission
 
 Description
-    Constant radiation absorption and emission coefficients
+    Constant radiation absorption and emission coefficients for continuous
+    phase
 
 SourceFiles
     constantAbsorptionEmission.C
@@ -92,18 +93,24 @@ public:
 
     // Member Operators
 
-        // Edit
-
         // Access
 
-            //- Return absorption coefficient
+            // Absorption coefficient
-            tmp<volScalarField> a() const;
 
-            //- Return emission coefficient
+                //- Absorption coefficient for continuous phase
-            tmp<volScalarField> e() const;
+                tmp<volScalarField> aCont() const;
 
-            //- Return emission contribution
+
-            tmp<volScalarField> E() const;
+            // Emission coefficient
+
+                //- Emission coefficient for continuous phase
+                tmp<volScalarField> eCont() const;
+
+
+            // Emission contribution
+
+                //- Emission contribution for continuous phase
+                tmp<volScalarField> ECont() const;
 };
 
 

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/noAbsorptionEmission/noAbsorptionEmission.C

@@ -63,75 +63,4 @@ Foam::radiation::noAbsorptionEmission::~noAbsorptionEmission()
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField>
-Foam::radiation::noAbsorptionEmission::a() const
-{
-    return tmp<volScalarField>
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "a",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh_,
-            dimensionedScalar("a", dimless/dimLength, 0.0)
-        )
-    );
-}
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::radiation::noAbsorptionEmission::e() const
-{
-    return tmp<volScalarField>
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "e",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh_,
-            dimensionedScalar("e", dimless/dimLength, 0.0)
-        )
-    );
-}
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::radiation::noAbsorptionEmission::E() const
-{
-    return tmp<volScalarField>
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "E",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh_,
-            dimensionedScalar("E", dimMass/dimLength/pow3(dimTime), 0.0)
-        )
-    );
-}
-
-
 // ************************************************************************* //

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/noAbsorptionEmission/noAbsorptionEmission.H

@@ -54,12 +54,6 @@ class noAbsorptionEmission
     public absorptionEmissionModel
 {
 
-    // Private data
-
-        //- Coefficients dictionary
-        dictionary coeffsDict_;
-
-
 public:
 
     //- Runtime type information
@@ -79,22 +73,6 @@ public:
     // Destructor
 
         ~noAbsorptionEmission();
-
-
-    // Member Operators
-
-        // Edit
-
-        // Access
-
-            //- Return absorption coefficient
-            tmp<volScalarField> a() const;
-
-            //- Return emission coefficient
-            tmp<volScalarField> e() const;
-
-            //- Return emission contribution
-            tmp<volScalarField> E() const;
 };