Added Larsen Borgnakke internal energy redistribution and all supporting variables and function calls. Added energy and momentum monitoring functions. Added U and T fields to be used for boundary conditions and field measurement. Modified constructors accordingly. Now reading dsmcInitialiseDict in the Cloud, not the application. Initialisation dict now reads a subdict of <species keyword> <number density> entries.

applications/solvers/dsmc/dsmcFoam/createFields.H

@@ -1,4 +1,60 @@
 
+    // volScalarField rhoN
+    // (
+    //     IOobject
+    //     (
+    //         "rhoN",
+    //         runTime.timeName(),
+    //         mesh,
+    //         IOobject::MUST_READ,
+    //         IOobject::AUTO_WRITE
+    //     ),
+    //     mesh
+    // );
+
+
+    // volScalarField rhoM
+    // (
+    //     IOobject
+    //     (
+    //         "rhoM",
+    //         runTime.timeName(),
+    //         mesh,
+    //         IOobject::MUST_READ,
+    //         IOobject::AUTO_WRITE
+    //     ),
+    //     mesh
+    // );
+
+
+    Info<< "\nReading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "\nReading field T\n" << endl;
+    volScalarField T
+    (
+        IOobject
+        (
+            "T",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
     Info<< "Constructing dsmcCloud " << endl;
 
-    dsmcCloud dsmc("dsmc", mesh);
+    dsmcCloud dsmc("dsmc", T, U);

applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C

@@ -26,8 +26,8 @@ Application
     dsmcFoam
 
 Description
-    Initialise a case for dsmcFoam by reading the initialise subdictionary in
+    Initialise a case for dsmcFoam by reading the initialisation dictionary
-    dsmcProperties
+    system/dsmcInitialise
 
 \*---------------------------------------------------------------------------*/
 
@@ -42,21 +42,11 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
 
-    IOdictionary dsmcInitialiseDict
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-    (
-        IOobject
-        (
-            "dsmcInitialiseDict",
-            mesh.time().system(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
 
     Info<< "Initialising dsmc for Time = " << runTime.timeName() << nl << endl;
 
-    dsmcCloud dsmc("dsmc", mesh, dsmcInitialiseDict);
+    dsmcCloud dsmc("dsmc", mesh);
 
     label totalMolecules = dsmc.size();
 

src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C

@@ -41,8 +41,7 @@ Foam::scalar Foam::DsmcCloud<ParcelType>::Tref = 273;
 template<class ParcelType>
 void Foam::DsmcCloud<ParcelType>::buildConstProps()
 {
-    Info<< nl << "typeIds found " << typeIdList_ << endl;
+    Info<< nl << "Constructing constant properties for" << endl;
-
     constProps_.setSize(typeIdList_.size());
 
     dictionary moleculeProperties
@@ -54,6 +53,8 @@ void Foam::DsmcCloud<ParcelType>::buildConstProps()
     {
         const word& id(typeIdList_[i]);
 
+        Info<< "    " << id << endl;
+
         const dictionary& molDict(moleculeProperties.subDict(id));
 
         constProps_[i] =
@@ -65,7 +66,6 @@ void Foam::DsmcCloud<ParcelType>::buildConstProps()
 template<class ParcelType>
 void Foam::DsmcCloud<ParcelType>::buildCellOccupancy()
 {
-
     forAll(cellOccupancy_, cO)
     {
         cellOccupancy_[cO].clear();
@@ -91,21 +91,21 @@ void Foam::DsmcCloud<ParcelType>::initialise
 
     const vector velocity(dsmcInitialiseDict.lookup("velocity"));
 
-    List<word> molecules
+    const dictionary& numberDensitiesDict
     (
-        dsmcInitialiseDict.lookup("molecules")
+        dsmcInitialiseDict.subDict("numberDensities")
     );
 
-    Field<scalar> numberDensities
+    List<word> molecules(numberDensitiesDict.toc());
-    (
-        dsmcInitialiseDict.lookup("numberDensities")
-    );
 
-    if(molecules.size() != numberDensities.size())
+    Field<scalar> numberDensities(molecules.size());
+
+    forAll(molecules, i)
     {
-        FatalErrorIn("Foam::Foam::DsmcCloud<ParcelType>::initialise")
+        numberDensities[i] = readScalar
-            << "molecules and numberDensities must be the same size."
+        (
-            << nl << abort(FatalError);
+            numberDensitiesDict.lookup(molecules[i])
+        );
     }
 
     numberDensities /= nParticle_;
@@ -171,6 +171,10 @@ void Foam::DsmcCloud<ParcelType>::initialise
                         cP.mass()
                     );
 
+                    scalar Ei =
+                        0.5*cP.internalDegreesOfFreedom()
+                       *kb*temperature;
+
                     U += velocity;
 
                     if (cell >= 0)
@@ -179,6 +183,7 @@ void Foam::DsmcCloud<ParcelType>::initialise
                         (
                             p,
                             U,
+                            Ei,
                             cell,
                             typeId
                         );
@@ -282,7 +287,9 @@ void Foam::DsmcCloud<ParcelType>::collisions()
                         typeIdP,
                         typeIdQ,
                         parcelP->U(),
-                        parcelQ->U()
+                        parcelQ->U(),
+                        parcelP->Ei(),
+                        parcelQ->Ei()
                     );
 
                     collisions++;
@@ -295,11 +302,18 @@ void Foam::DsmcCloud<ParcelType>::collisions()
 
     reduce(collisionCandidates, sumOp<label>());
 
-    Info<< "    Collisions                      = "
+    if (collisionCandidates)
-        << collisions << nl
+    {
-        << "    Acceptance rate                 = "
+        Info<< "    Collisions                      = "
-        << scalar(collisions)/scalar(collisionCandidates) << nl
+            << collisions << nl
-        << endl;
+            << "    Acceptance rate                 = "
+            << scalar(collisions)/scalar(collisionCandidates) << nl
+            << endl;
+    }
+    else
+    {
+        Info<< "    No collisions" << endl;
+    }
 }
 
 
@@ -310,6 +324,7 @@ void Foam::DsmcCloud<ParcelType>::addNewParcel
 (
     const vector& position,
     const vector& U,
+    const scalar Ei,
     const label cellId,
     const label typeId
 )
@@ -319,6 +334,7 @@ void Foam::DsmcCloud<ParcelType>::addNewParcel
         *this,
         position,
         U,
+        Ei,
         cellId,
         typeId
     );
@@ -333,20 +349,21 @@ template<class ParcelType>
 Foam::DsmcCloud<ParcelType>::DsmcCloud
 (
     const word& cloudType,
-    const fvMesh& mesh
+    const volScalarField& T,
+    const volVectorField& U
 )
 :
-    Cloud<ParcelType>(mesh, cloudType, false),
+    Cloud<ParcelType>(T.mesh(), cloudType, false),
     DsmcBaseCloud(),
     cloudType_(cloudType),
-    mesh_(mesh),
+    mesh_(T.mesh()),
     particleProperties_
     (
         IOobject
         (
             cloudType + "Properties",
-            mesh.time().constant(),
+            mesh_.time().constant(),
-            mesh,
+            mesh_,
             IOobject::MUST_READ,
             IOobject::NO_WRITE
         )
@@ -369,6 +386,8 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
     collisionSelectionRemainder_(mesh_.nCells(), 0),
     constProps_(),
     rndGen_(label(971501)),
+    T_(T),
+    U_(U),
     binaryCollisionModel_
     (
         BinaryCollisionModel<DsmcCloud<ParcelType> >::New
@@ -396,8 +415,7 @@ template<class ParcelType>
 Foam::DsmcCloud<ParcelType>::DsmcCloud
 (
     const word& cloudType,
-    const fvMesh& mesh,
+    const fvMesh& mesh
-    const IOdictionary& dsmcInitialiseDict
 )
     :
     Cloud<ParcelType>(mesh, cloudType, false),
@@ -409,8 +427,8 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
         IOobject
         (
             cloudType + "Properties",
-            mesh.time().constant(),
+            mesh_.time().constant(),
-            mesh,
+            mesh_,
             IOobject::MUST_READ,
             IOobject::NO_WRITE
         )
@@ -434,6 +452,43 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
     collisionSelectionRemainder_(),
     constProps_(),
     rndGen_(label(971501)),
+    T_
+    (
+        volScalarField
+        (
+            IOobject
+            (
+                "T",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero",  dimensionSet(0, 0, 0, 1, 0), 0.0)
+        )
+    ),
+    U_
+    (
+        volVectorField
+        (
+            IOobject
+            (
+                "U",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedVector
+            (
+                "zero",
+                dimensionSet(0, 1, -1, 0, 0),
+                vector::zero
+            )
+        )
+    ),
     binaryCollisionModel_(),
     wallInteractionModel_()
 {
@@ -441,6 +496,18 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
 
     buildConstProps();
 
+    IOdictionary dsmcInitialiseDict
+    (
+        IOobject
+        (
+            "dsmcInitialiseDict",
+            mesh_.time().system(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
     initialise(dsmcInitialiseDict);
 }
 
@@ -478,9 +545,14 @@ template<class ParcelType>
 void Foam::DsmcCloud<ParcelType>::info() const
 {
     vector linearMomentum = linearMomentumOfSystem();
-
     reduce(linearMomentum, sumOp<vector>());
 
+    scalar linearKineticEnergy = linearKineticEnergyOfSystem();
+    reduce(linearKineticEnergy, sumOp<scalar>());
+
+    scalar internalEnergy = internalEnergyOfSystem();
+    reduce(internalEnergy, sumOp<scalar>());
+
     Info<< "Cloud name: " << this->name() << nl
         << "    Current number of parcels       = "
         << returnReduce(this->size(), sumOp<label>()) << nl
@@ -491,7 +563,11 @@ void Foam::DsmcCloud<ParcelType>::info() const
         << "    Linear momentum magnitude       = "
         << mag(linearMomentum) << nl
         << "    Linear kinetic energy           = "
-        << returnReduce(linearKineticEnergyOfSystem(), sumOp<scalar>()) << nl
+        << linearKineticEnergy << nl
+        << "    Internal energy                 = "
+        << internalEnergy << nl
+        << "    Total energy                    = "
+        << internalEnergy + linearKineticEnergy << nl
         << endl;
 }
 

src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H

@@ -107,9 +107,19 @@ class DsmcCloud
         //- Random number generator
         Random rndGen_;
 
+
+        // References to the macroscopic fields
+
+            //- Temperature
+            const volScalarField& T_;
+
+            //- Velocity
+            const volVectorField& U_;
+
+
         // References to the cloud sub-models
 
-            //- Binary  Collision model
+            //- Binary collision model
             autoPtr<BinaryCollisionModel<DsmcCloud<ParcelType> > >
                 binaryCollisionModel_;
 
@@ -156,16 +166,15 @@ public:
         DsmcCloud
         (
             const word& cloudType,
-            const fvMesh& mesh
+            const volScalarField& T,
+            const volVectorField& U
         );
 
-        //- Construct given name, mesh, and typeIdListDict.  Used by
+        //- Construct given name and mesh.  Used to initialise.
-        // dsmcInitialse.
         DsmcCloud
         (
             const word& cloudType,
-            const fvMesh& mesh,
+            const fvMesh& mesh
-            const IOdictionary& dsmcInitialiseDict
         );
 
 
@@ -185,6 +194,9 @@ public:
                 //- Return refernce to the mesh
                 inline const fvMesh& mesh() const;
 
+
+            // References to the dsmc specific data
+
                 //- Return particle properties dictionary
                 inline const IOdictionary& particleProperties() const;
 
@@ -219,6 +231,15 @@ public:
                 inline Random& rndGen();
 
 
+            // References to the macroscopic fields
+
+                //- Return macroscopic temperature
+                inline const volScalarField& T() const;
+
+                //- Return macroscopic velocity
+                inline const volVectorField& U() const;
+
+
             // Kinetic theory helper functions
 
                 //- Generate a random velocity sampled from the Maxwellian speed
@@ -252,9 +273,6 @@ public:
                     scalar mass
                 ) const;
 
-            //- Clear the Cloud
-            inline void clear();
-
 
             // Sub-models
 
@@ -289,6 +307,9 @@ public:
             //- Total linear kinetic energy in the system
             inline scalar linearKineticEnergyOfSystem() const;
 
+            //- Total internal energy in the system
+            inline scalar internalEnergyOfSystem() const;
+
             //- Print cloud information
             void info() const;
 
@@ -307,12 +328,6 @@ public:
                 //- Return the field of number of DSMC particles
                 inline const tmp<volScalarField> rhoNP() const;
 
-                //- Return the velocity field
-                inline const tmp<volVectorField> U() const;
-
-                //- Return the temperature field
-                inline const tmp<volScalarField> totalKE() const;
-
 
         // Cloud evolution functions
 
@@ -321,12 +336,19 @@ public:
             (
                 const vector& position,
                 const vector& U,
+                const scalar Ei,
                 const label cellId,
                 const label typeId
             );
 
             //- Evolve the cloud (move, collide)
             void evolve();
+
+
+            //- Clear the Cloud
+            inline void clear();
+
+
 };
 
 

src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H

@@ -48,7 +48,6 @@ Foam::DsmcCloud<ParcelType>::particleProperties() const
 }
 
 
-
 template<class ParcelType>
 inline const Foam::List<Foam::word>&
 Foam::DsmcCloud<ParcelType>::typeIdList() const
@@ -115,9 +114,23 @@ Foam::DsmcCloud<ParcelType>::constProps
 
 
 template<class ParcelType>
-inline void Foam::DsmcCloud<ParcelType>::clear()
+inline Foam::Random& Foam::DsmcCloud<ParcelType>::rndGen()
 {
-    return IDLList<ParcelType>::clear();
+    return rndGen_;
+}
+
+
+template<class ParcelType>
+inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::T() const
+{
+    return T_;
+}
+
+
+template<class ParcelType>
+inline const Foam::volVectorField& Foam::DsmcCloud<ParcelType>::U() const
+{
+    return U_;
 }
 
 
@@ -158,9 +171,19 @@ inline Foam::scalar Foam::DsmcCloud<ParcelType>::massInSystem() const
 {
     scalar sysMass = 0.0;
 
-    Info<< "massInSystem() not implemented" << endl;
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
 
-    return sysMass;
+        const typename ParcelType::constantProperties& cP = constProps
+        (
+            p.typeId()
+        );
+
+        sysMass += cP.mass();
+    }
+
+    return nParticle_*sysMass;
 }
 
 
@@ -208,9 +231,19 @@ Foam::DsmcCloud<ParcelType>::linearKineticEnergyOfSystem() const
 
 
 template<class ParcelType>
-inline Foam::Random& Foam::DsmcCloud<ParcelType>::rndGen()
+inline Foam::scalar
+Foam::DsmcCloud<ParcelType>::internalEnergyOfSystem() const
 {
-    return rndGen_;
+    scalar internalEnergy = 0.0;
+
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+
+        internalEnergy += p.Ei();
+    }
+
+    return nParticle_*internalEnergy;
 }
 
 
@@ -327,96 +360,10 @@ Foam::DsmcCloud<ParcelType>::rhoNP() const
 
 
 template<class ParcelType>
-inline const Foam::tmp<Foam::volVectorField>
+inline void Foam::DsmcCloud<ParcelType>::clear()
-Foam::DsmcCloud<ParcelType>::U() const
 {
-    tmp<volScalarField> tU
+    return IDLList<ParcelType>::clear();
-    (
-        new volVectorField
-        (
-            IOobject
-            (
-                this->name() + "U",
-                this->db().time().timeName(),
-                this->db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh_,
-            dimensionedVector
-            (
-                "zero",
-                dimensionSet(0, 1, -1, 0, 0),
-                vector::zero
-            )
-        )
-    );
-
-    return tU;
 }
 
 
-template<class ParcelType>
-inline const Foam::tmp<Foam::volScalarField>
-Foam::DsmcCloud<ParcelType>::totalKE() const
-{
-    tmp<volScalarField> ttotalKE
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                this->name() + "totalKE",
-                this->db().time().timeName(),
-                this->db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh_,
-            dimensionedScalar("zero", dimensionSet(0, 0, 0, 1, 0), 0.0)
-        )
-    );
-
-    return ttotalKE;
-}
-
-// template<class ParcelType>
-// inline const Foam::tmp<Foam::volScalarField>
-// Foam::DsmcCloud<ParcelType>::alpha() const
-// {
-//     tmp<volScalarField> talpha
-//     (
-//         new volScalarField
-//         (
-//             IOobject
-//             (
-//                 this->name() + "Alpha",
-//                 this->db().time().timeName(),
-//                 this->db(),
-//                 IOobject::NO_READ,
-//                 IOobject::NO_WRITE,
-//                 false
-//             ),
-//             mesh_,
-//             dimensionedScalar("zero", dimensionSet(0, 0, 0, 0, 0), 0.0)
-//         )
-//     );
-
-//     scalarField& alpha = talpha().internalField();
-//     forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
-//     {
-//         const ParcelType& p = iter();
-//         const label cellI = p.cell();
-
-//         alpha[cellI] += p.nParticle()*p.mass();
-//     }
-
-//     alpha /= (mesh().cellVolumes()*rho_);
-
-//     return talpha;
-// }
-
-
 // ************************************************************************* //

src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C

@@ -39,10 +39,11 @@ namespace Foam
 Foam::dsmcCloud::dsmcCloud
 (
     const word& cloudType,
-    const fvMesh& mesh
+    const volScalarField& T,
+    const volVectorField& U
 )
 :
-    DsmcCloud<dsmcParcel>(cloudType, mesh)
+    DsmcCloud<dsmcParcel>(cloudType, T, U)
 {
     dsmcParcel::readFields(*this);
 }
@@ -51,11 +52,10 @@ Foam::dsmcCloud::dsmcCloud
 Foam::dsmcCloud::dsmcCloud
 (
     const word& cloudType,
-    const fvMesh& mesh,
+    const fvMesh& mesh
-    const IOdictionary& dsmcInitialiseDict
 )
 :
-    DsmcCloud<dsmcParcel>(cloudType, mesh, dsmcInitialiseDict)
+    DsmcCloud<dsmcParcel>(cloudType, mesh)
 {}
 
 

src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H

@@ -73,15 +73,15 @@ public:
         dsmcCloud
         (
             const word& cloudType,
-            const fvMesh& mesh
+            const volScalarField& T,
+            const volVectorField& U
         );
 
-        //- Construct from with the typeIdListDict, used to initialise
+        //- Construct from name and mesh, used to initialise.
         dsmcCloud
         (
             const word& cloudType,
-            const fvMesh& mesh,
+            const fvMesh& mesh
-            const IOdictionary& dsmcInitialiseDict
         );
 
     //- Destructor

src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C

@@ -65,8 +65,7 @@ bool Foam::DsmcParcel<ParcelType>::move
             {
                 if
                 (
-                    isType<processorPolyPatch>
+                    isType<processorPolyPatch>(pbMesh[p.patch(p.face())])
-                        (pbMesh[p.patch(p.face())])
                 )
                 {
                     td.switchProcessor = true;

src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H

@@ -162,6 +162,10 @@ protected:
             //- Velocity of Parcel [m/s]
             vector U_;
 
+            //- Internal energy of the Parcel, covering all non-translational
+            // degrees of freedom [J]
+            scalar Ei_;
+
             //- Parcel type id
             label typeId_;
 
@@ -182,6 +186,7 @@ public:
             DsmcCloud<ParcelType>& owner,
             const vector& position,
             const vector& U,
+            const scalar Ei,
             const label celli,
             const label typeId
         );
@@ -211,11 +216,17 @@ public:
             //- Return const access to velocity
             inline const vector& U() const;
 
+            //- Return const access to internal energy
+            inline scalar Ei() const;
+
         // Edit
 
             //- Return access to velocity
             inline vector& U();
 
+            //- Return access to internal energy
+            inline scalar& Ei();
+
 
         // Main calculation loop
 

src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelI.H

@@ -67,12 +67,14 @@ inline Foam::DsmcParcel<ParcelType>::DsmcParcel
     DsmcCloud<ParcelType>& owner,
     const vector& position,
     const vector& U,
+    const scalar Ei,
     const label celli,
     const label typeId
 )
 :
     Particle<ParcelType>(owner, position, celli),
     U_(U),
+    Ei_(Ei),
     typeId_(typeId)
 {}
 
@@ -146,6 +148,13 @@ inline const Foam::vector& Foam::DsmcParcel<ParcelType>::U() const
 }
 
 
+template <class ParcelType>
+inline Foam::scalar Foam::DsmcParcel<ParcelType>::Ei() const
+{
+    return Ei_;
+}
+
+
 template <class ParcelType>
 inline Foam::vector& Foam::DsmcParcel<ParcelType>::U()
 {
@@ -153,4 +162,11 @@ inline Foam::vector& Foam::DsmcParcel<ParcelType>::U()
 }
 
 
+template <class ParcelType>
+inline Foam::scalar& Foam::DsmcParcel<ParcelType>::Ei()
+{
+    return Ei_;
+}
+
+
 // ************************************************************************* //

src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelIO.C

@@ -41,6 +41,7 @@ Foam::DsmcParcel<ParcelType>::DsmcParcel
 :
     Particle<ParcelType>(cloud, is, readFields),
     U_(vector::zero),
+    Ei_(0.0),
     typeId_(-1)
 {
     if (readFields)
@@ -48,6 +49,7 @@ Foam::DsmcParcel<ParcelType>::DsmcParcel
         if (is.format() == IOstream::ASCII)
         {
             is >> U_;
+            Ei_ = readScalar(is);
             typeId_ = readLabel(is);
         }
         else
@@ -56,6 +58,7 @@ Foam::DsmcParcel<ParcelType>::DsmcParcel
             (
                 reinterpret_cast<char*>(&U_),
                 sizeof(U_)
+              + sizeof(Ei_)
               + sizeof(typeId_)
             );
         }
@@ -84,6 +87,9 @@ void Foam::DsmcParcel<ParcelType>::readFields
     IOField<vector> U(c.fieldIOobject("U", IOobject::MUST_READ));
     c.checkFieldIOobject(c, U);
 
+    IOField<scalar> Ei(c.fieldIOobject("Ei", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, Ei);
+
     IOField<label> typeId(c.fieldIOobject("typeId", IOobject::MUST_READ));
     c.checkFieldIOobject(c, typeId);
 
@@ -93,6 +99,7 @@ void Foam::DsmcParcel<ParcelType>::readFields
         ParcelType& p = iter();
 
         p.U_ = U[i];
+        p.Ei_ = Ei[i];
         p.typeId_ = typeId[i];
         i++;
     }
@@ -110,6 +117,7 @@ void Foam::DsmcParcel<ParcelType>::writeFields
     label np =  c.size();
 
     IOField<vector> U(c.fieldIOobject("U", IOobject::NO_READ), np);
+    IOField<scalar> Ei(c.fieldIOobject("Ei", IOobject::NO_READ), np);
     IOField<label> typeId(c.fieldIOobject("typeId", IOobject::NO_READ), np);
 
     label i = 0;
@@ -118,11 +126,13 @@ void Foam::DsmcParcel<ParcelType>::writeFields
         const DsmcParcel<ParcelType>& p = iter();
 
         U[i] = p.U();
+        Ei[i] = p.Ei();
         typeId[i] = p.typeId();
         i++;
     }
 
     U.write();
+    Ei.write();
     typeId.write();
 }
 
@@ -140,6 +150,7 @@ Foam::Ostream& Foam::operator<<
     {
         os  << static_cast<const Particle<ParcelType>& >(p)
             << token::SPACE << p.U()
+            << token::SPACE << p.Ei()
             << token::SPACE << p.typeId();
     }
     else
@@ -149,6 +160,7 @@ Foam::Ostream& Foam::operator<<
         (
             reinterpret_cast<const char*>(&p.U_),
             sizeof(p.U())
+          + sizeof(p.Ei())
           + sizeof(p.typeId())
         );
     }

src/lagrangian/dsmc/parcels/derived/dsmcParcel/dsmcParcel.C

@@ -43,6 +43,7 @@ Foam::dsmcParcel::dsmcParcel
     DsmcCloud<dsmcParcel>& owner,
     const vector& position,
     const vector& U,
+    const scalar Ei,
     const label celli,
     const label typeId
 )
@@ -52,6 +53,7 @@ Foam::dsmcParcel::dsmcParcel
         owner,
         position,
         U,
+        Ei,
         celli,
         typeId
     )

src/lagrangian/dsmc/parcels/derived/dsmcParcel/dsmcParcel.H

@@ -66,6 +66,7 @@ public:
             DsmcCloud<dsmcParcel>& owner,
             const vector& position,
             const vector& U,
+            const scalar Ei,
             const label celli,
             const label typeId
         );

src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelBinaryCollisionModels.C

@@ -27,6 +27,7 @@ License
 #include "dsmcParcel.H"
 #include "DsmcCloud.H"
 #include "VariableHardSphere.H"
+#include "LarsenBorgnakkeVariableHardSphere.H"
 
 namespace Foam
 {
@@ -39,6 +40,12 @@ namespace Foam
         DsmcCloud,
         dsmcParcel
     );
+    makeBinaryCollisionModelType
+    (
+        LarsenBorgnakkeVariableHardSphere,
+        DsmcCloud,
+        dsmcParcel
+    );
 };
 
 

src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.H

@@ -139,7 +139,9 @@ public:
             label typeIdP,
             label typeIdQ,
             vector& UP,
-            vector& UQ
+            vector& UQ,
+            scalar& EiP,
+            scalar& EiQ
         ) = 0;
 };
 

src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C

@@ -0,0 +1,265 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 "LarsenBorgnakkeVariableHardSphere.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::scalar Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::energyRatio
+(
+    scalar ChiA,
+    scalar ChiB
+)
+{
+    CloudType& cloud(this->owner());
+
+    Random& rndGen(cloud.rndGen());
+
+    scalar ChiAMinusOne = ChiA - 1;
+
+    scalar ChiBMinusOne = ChiB - 1;
+
+    if (ChiAMinusOne < SMALL && ChiBMinusOne < SMALL)
+    {
+        return rndGen.scalar01();
+    }
+
+    scalar energyRatio;
+
+    scalar P;
+
+    do
+    {
+        P = 0;
+
+        energyRatio = rndGen.scalar01();
+
+        if (ChiAMinusOne < SMALL)
+        {
+            P = 1.0 - pow(energyRatio, ChiB);
+        }
+        else if (ChiBMinusOne < SMALL)
+        {
+            P = 1.0 - pow(energyRatio, ChiA);
+        }
+        else
+        {
+            P = pow
+            (
+                (ChiAMinusOne + ChiBMinusOne)*energyRatio/ChiAMinusOne,
+                ChiAMinusOne
+            )
+           *pow
+            (
+                (ChiAMinusOne + ChiBMinusOne)*(1 - energyRatio)/ChiBMinusOne,
+                ChiBMinusOne
+            );
+        }
+    }
+
+    while (P < rndGen.scalar01());
+
+    return energyRatio;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::LarsenBorgnakkeVariableHardSphere
+(
+    const dictionary& dict,
+    CloudType& cloud
+)
+:
+    BinaryCollisionModel<CloudType>(dict, cloud, typeName),
+    relaxationCollisionNumber_
+    (
+        readScalar(this->coeffDict().lookup("relaxationCollisionNumber"))
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::~LarsenBorgnakkeVariableHardSphere()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+template <class CloudType>
+Foam::scalar Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::sigmaTcR
+(
+    label typeIdP,
+    label typeIdQ,
+    const vector& UP,
+    const vector& UQ
+) const
+{
+    const CloudType& cloud(this->owner());
+
+    scalar dPQ =
+        0.5
+       *(
+            cloud.constProps(typeIdP).d()
+          + cloud.constProps(typeIdQ).d()
+        );
+
+    scalar omegaPQ =
+        0.5
+       *(
+            cloud.constProps(typeIdP).omega()
+          + cloud.constProps(typeIdQ).omega()
+        );
+
+    scalar cR = mag(UP - UQ);
+
+    scalar mP = cloud.constProps(typeIdP).mass();
+
+    scalar mQ = cloud.constProps(typeIdQ).mass();
+
+    scalar mR = mP*mQ/(mP + mQ);
+
+    // calculating cross section = pi*dPQ^2, where dPQ is from Bird, eq. 4.79
+    scalar sigmaTPQ =
+        mathematicalConstant::pi*dPQ*dPQ
+       *pow(2.0*CloudType::kb*CloudType::Tref/(mR*cR*cR), omegaPQ - 0.5)
+       /exp(Foam::lgamma(2.5 - omegaPQ));
+
+    return sigmaTPQ*cR;
+}
+
+
+template <class CloudType>
+void Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::collide
+(
+    label typeIdP,
+    label typeIdQ,
+    vector& UP,
+    vector& UQ,
+    scalar& EiP,
+    scalar& EiQ
+)
+{
+    CloudType& cloud(this->owner());
+
+    Random& rndGen(cloud.rndGen());
+
+    scalar inverseCollisionNumber = 1/relaxationCollisionNumber_;
+
+    // Larsen Borgnakke internal energy redistribution part.  Using the serial
+    // application of the LB method, as per the INELRS subroutine in Bird's
+    // DSMC0R.FOR
+
+    scalar preCollisionEiP = EiP;
+
+    scalar preCollisionEiQ = EiQ;
+
+    scalar iDofP = cloud.constProps(typeIdP).internalDegreesOfFreedom();
+
+    scalar iDofQ = cloud.constProps(typeIdQ).internalDegreesOfFreedom();
+
+    scalar omegaPQ =
+        0.5
+       *(
+            cloud.constProps(typeIdP).omega()
+          + cloud.constProps(typeIdQ).omega()
+        );
+
+    scalar mP = cloud.constProps(typeIdP).mass();
+
+    scalar mQ = cloud.constProps(typeIdQ).mass();
+
+    scalar mR = mP*mQ/(mP + mQ);
+
+    vector Ucm = (mP*UP + mQ*UQ)/(mP + mQ);
+
+    scalar cRsqr = magSqr(UP - UQ);
+
+    scalar availableEnergy = 0.5*mR*cRsqr;
+
+    scalar ChiB = 2.5 - omegaPQ;
+
+    if (iDofP > 0)
+    {
+        if (inverseCollisionNumber > rndGen.scalar01())
+        {
+            availableEnergy += preCollisionEiP;
+
+            scalar ChiA = 0.5*iDofP;
+
+            EiP = energyRatio(ChiA, ChiB)*availableEnergy;
+
+            availableEnergy -= EiP;
+        }
+    }
+
+    if (iDofQ > 0)
+    {
+        if (inverseCollisionNumber > rndGen.scalar01())
+        {
+            availableEnergy += preCollisionEiQ;
+
+            // Change to general LB ratio calculation
+            scalar energyRatio = 1.0 - pow(rndGen.scalar01(),(1.0/ChiB));
+
+            EiQ = energyRatio*availableEnergy;
+
+            availableEnergy -= EiQ;
+        }
+    }
+
+    // Rescale the translational energy
+    scalar cR = sqrt(2.0*availableEnergy/mR);
+
+    // Variable Hard Sphere collision part
+
+    scalar cosTheta = 2.0*rndGen.scalar01() - 1.0;
+
+    scalar sinTheta = sqrt(1.0 - cosTheta*cosTheta);
+
+    scalar phi = 2.0*mathematicalConstant::pi*rndGen.scalar01();
+
+    vector postCollisionRelU =
+        cR
+       *vector
+        (
+            cosTheta,
+            sinTheta*cos(phi),
+            sinTheta*sin(phi)
+        );
+
+    UP = Ucm + postCollisionRelU*mQ/(mP + mQ);
+
+    UQ = Ucm - postCollisionRelU*mP/(mP + mQ);
+}
+
+
+// ************************************************************************* //

src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.H

@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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::LarsenBorgnakkeVariableHardSphere
+
+Description
+    Variable Hard Sphere BinaryCollision Model with Larsen Borgnakke internal
+    energy redistribution.  Based on the INELRS subroutine in Bird's DSMC0R.FOR
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LarsenBorgnakkeVariableHardSphere_H
+#define LarsenBorgnakkeVariableHardSphere_H
+
+#include "BinaryCollisionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                      Class LarsenBorgnakkeVariableHardSphere Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class LarsenBorgnakkeVariableHardSphere
+:
+    public BinaryCollisionModel<CloudType>
+{
+    // Private data
+
+        //- Temperature
+        const scalar relaxationCollisionNumber_;
+
+
+    // Private Member Functions
+
+        //- Calculate the energy ratio for distribution to internal degrees of
+        // freedom
+        scalar energyRatio
+        (
+            scalar ChiA,
+            scalar ChiB
+        );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("LarsenBorgnakkeVariableHardSphere");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        LarsenBorgnakkeVariableHardSphere
+        (
+            const dictionary& dict,
+            CloudType& cloud
+        );
+
+
+    // Destructor
+    virtual ~LarsenBorgnakkeVariableHardSphere();
+
+
+    // Member Functions
+
+        //- Return the collision cross section * relative velocity product
+        virtual scalar sigmaTcR
+        (
+            label typeIdP,
+            label typeIdQ,
+            const vector& UP,
+            const vector& UQ
+        ) const;
+
+        //- Apply collision
+        virtual void collide
+        (
+            label typeIdP,
+            label typeIdQ,
+            vector& UP,
+            vector& UQ,
+            scalar& EiP,
+            scalar& EiQ
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "LarsenBorgnakkeVariableHardSphere.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C

@@ -98,7 +98,9 @@ void Foam::VariableHardSphere<CloudType>::collide
     label typeIdP,
     label typeIdQ,
     vector& UP,
-    vector& UQ
+    vector& UQ,
+    scalar& EiP,
+    scalar& EiQ
 )
 {
     CloudType& cloud(this->owner());

src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.H

@@ -84,7 +84,9 @@ public:
             label typeIdP,
             label typeIdQ,
             vector& UP,
-            vector& UQ
+            vector& UQ,
+            scalar& EiP,
+            scalar& EiQ
         );
 };