Added equipartition internal energy distribution sampling function and calling...

Added equipartition internal energy distribution sampling function and calling on particle intialisation, Maxwellian wall collision and free stream injection.  Interface with WallInteractionModel modified to pass a scalar& for the internal energy and a label to specify the typeId.

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

@@ -172,9 +172,11 @@ void Foam::DsmcCloud<ParcelType>::initialise
                         cP.mass()
                     );
 
-                    scalar Ei =
-                        0.5*cP.internalDegreesOfFreedom()
-                       *kb*temperature;
+                    scalar Ei = equipartitionInternalEnergy
+                    (
+                        temperature,
+                        cP.internalDegreesOfFreedom()
+                    );
 
                     U += velocity;
 
@@ -690,6 +692,47 @@ Foam::vector Foam::DsmcCloud<ParcelType>::equipartitionLinearVelocity
 }
 
 
+template<class ParcelType>
+Foam::scalar Foam::DsmcCloud<ParcelType>::equipartitionInternalEnergy
+(
+    scalar temperature,
+    scalar iDof
+)
+{
+    scalar Ei = 0.0;
+
+    if
+    (
+        iDof < 2.0 + SMALL
+     && iDof > 2.0 - SMALL
+    )
+    {
+        // Special case for iDof = 2, i.e. diatomics;
+        Ei = -log(rndGen_.scalar01())*kb*temperature;
+    }
+    else
+    {
+        scalar a = 0.5*iDof - 1;
+
+        scalar energyRatio;
+
+        scalar P;
+
+        do
+        {
+            energyRatio = 10*rndGen_.scalar01();
+
+            P = pow((energyRatio/a), a)*exp(a - energyRatio);
+
+        } while (P < rndGen_.scalar01());
+
+        Ei = energyRatio*kb*temperature;
+    }
+
+    return Ei;
+}
+
+
 template<class ParcelType>
 void Foam::DsmcCloud<ParcelType>::dumpParticlePositions() const
 {

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

@@ -281,8 +281,17 @@ public:
                     scalar mass
                 );
 
-                //- From the Maxwellian distribution:
+                //- Generate a random internal energy, sampled from the
+                // equilibrium distribution (Bird eqn 11.22 and 11.23 and
+                // adapting code from DSMC3.FOR)
+                scalar equipartitionInternalEnergy
+                (
+                    scalar temperature,
+                    scalar internalDegreesOfFreedom
+                );
+
 
+                // From the Maxwellian distribution:
                 //- Average particle speed
                 inline scalar maxwellianAverageSpeed
                 (

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

@@ -109,25 +109,28 @@ void Foam::DsmcParcel<ParcelType>::hitWallPatch
 {
     const constantProperties& constProps(td.cloud().constProps(typeId_));
 
+    scalar m = constProps.mass();
+
     // pre-interaction energy
-    scalar preIE = 0.5*constProps.mass()*(U_ & U_) + Ei_;
+    scalar preIE = 0.5*m*(U_ & U_) + Ei_;
 
     // pre-interaction momentum
-    vector preIMom = constProps.mass()*U_;
+    vector preIMom = m*U_;
 
     td.cloud().wallInteraction().correct
     (
         wpp,
         this->face(),
         U_,
-        constProps.mass()
+        Ei_,
+        typeId_
     );
 
     // post-interaction energy
-    scalar postIE = 0.5*constProps.mass()*(U_ & U_) + Ei_;
+    scalar postIE = 0.5*m*(U_ & U_) + Ei_;
 
     // post-interaction momentum
-    vector postIMom = constProps.mass()*U_;
+    vector postIMom = m*U_;
 
     label wppIndex = wpp.index();
 

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

@@ -81,9 +81,7 @@ Foam::scalar Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::energyRatio
                     ChiBMinusOne
                 );
         }
-    }
-
-    while (P < rndGen.scalar01());
+    } while (P < rndGen.scalar01());
 
     return energyRatio;
 }

src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C

@@ -237,9 +237,11 @@ void Foam::FreeStream<CloudType>::inflow()
 
                 U += velocity_;
 
-                scalar Ei =
-                    0.5*cloud.constProps(typeId).internalDegreesOfFreedom()
-                   *CloudType::kb*temperature_;
+                scalar Ei = cloud.equipartitionInternalEnergy
+                (
+                    temperature_,
+                    cloud.constProps(typeId).internalDegreesOfFreedom()
+                );
 
                 cloud.addNewParcel
                 (

src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C

@@ -54,7 +54,8 @@ void Foam::MaxwellianThermal<CloudType>::correct
     const wallPolyPatch& wpp,
     const label faceId,
     vector& U,
-    scalar mass
+    scalar& Ei,
+    label typeId
 )
 {
     label wppIndex = wpp.index();
@@ -102,6 +103,10 @@ void Foam::MaxwellianThermal<CloudType>::correct
 
     scalar T = cloud.T().boundaryField()[wppIndex][wppLocalFace];
 
+    scalar mass = cloud.constProps(typeId).mass();
+
+    scalar iDof = cloud.constProps(typeId).internalDegreesOfFreedom();
+
     U =
         sqrt(CloudType::kb*T/mass)
        *(
@@ -111,6 +116,8 @@ void Foam::MaxwellianThermal<CloudType>::correct
         );
 
     U += cloud.U().boundaryField()[wppIndex][wppLocalFace];
+
+    Ei = cloud.equipartitionInternalEnergy(T, iDof);
 }
 
 

src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.H

@@ -78,7 +78,8 @@ public:
             const wallPolyPatch& wpp,
             const label faceId,
             vector& U,
-            scalar mass
+            scalar& Ei,
+            label typeId
         );
 };
 

src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C

@@ -56,7 +56,8 @@ void Foam::SpecularReflection<CloudType>::correct
     const wallPolyPatch& wpp,
     const label faceId,
     vector& U,
-    scalar mass
+    scalar& Ei,
+    label typeId
 )
 {
     vector nw = wpp.faceAreas()[wpp.whichFace(faceId)];

src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.H

@@ -76,7 +76,8 @@ public:
             const wallPolyPatch& wpp,
             const label faceId,
             vector& U,
-            scalar mass
+            scalar& Ei,
+            label typeId
         );
 };
 

src/lagrangian/dsmc/submodels/WallInteractionModel/WallInteractionModel/WallInteractionModel.H

@@ -130,7 +130,8 @@ public:
             const wallPolyPatch& wpp,
             const label faceId,
             vector& U,
-            scalar mass
+            scalar& Ei,
+            label typeId
         ) = 0;
 };