DSMC: Corrected LarsenBorgnakkeVariableHardSphere model

Also changed internalDegreesOfFreedom to be an integer type

src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C

@@ -120,9 +120,7 @@ void Foam::DSMCCloud<ParcelType>::initialise
         forAll(cellTets, tetI)
         {
             const tetIndices& cellTetIs = cellTets[tetI];
-
             tetPointRef tet = cellTetIs.tet(mesh_);
-
             scalar tetVolume = tet.mag();
 
             forAll(molecules, i)
@@ -239,7 +237,6 @@ void Foam::DSMCCloud<ParcelType>::collisions()
 
         if (nC > 1)
         {
-
             // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
             // Assign particles to one of 8 Cartesian subCells
 
@@ -257,14 +254,12 @@ void Foam::DSMCCloud<ParcelType>::collisions()
             forAll(cellParcels, i)
             {
                 const ParcelType& p = *cellParcels[i];
-
                 vector relPos = p.position() - cC;
 
                 label subCell =
                     pos(relPos.x()) + 2*pos(relPos.y()) + 4*pos(relPos.z());
 
                 subCells[subCell].append(i);
-
                 whichSubCell[i] = subCell;
             }
 
@@ -278,9 +273,7 @@ void Foam::DSMCCloud<ParcelType>::collisions()
                /mesh_.cellVolumes()[cellI];
 
             label nCandidates(selectedPairs);
-
             collisionSelectionRemainder_[cellI] = selectedPairs - nCandidates;
-
             collisionCandidates += nCandidates;
 
             for (label c = 0; c < nCandidates; c++)
@@ -295,7 +288,6 @@ void Foam::DSMCCloud<ParcelType>::collisions()
                 label candidateQ = -1;
 
                 List<label> subCellPs = subCells[whichSubCell[candidateP]];
-
                 label nSC = subCellPs.size();
 
                 if (nSC > 1)
@@ -307,7 +299,6 @@ void Foam::DSMCCloud<ParcelType>::collisions()
                     do
                     {
                         candidateQ = subCellPs[rndGen_.integer(0, nSC - 1)];
-
                     } while (candidateP == candidateQ);
                 }
                 else
@@ -319,7 +310,6 @@ void Foam::DSMCCloud<ParcelType>::collisions()
                     do
                     {
                         candidateQ = rndGen_.integer(0, nC - 1);
-
                     } while (candidateP == candidateQ);
                 }
 
@@ -406,15 +396,10 @@ void Foam::DSMCCloud<ParcelType>::resetFields()
     );
 
     rhoN_ = dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL);
-
     rhoM_ =  dimensionedScalar("zero",  dimensionSet(1, -3, 0, 0, 0), VSMALL);
-
     dsmcRhoN_ = dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), 0.0);
-
     linearKE_ = dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0);
-
     internalE_ = dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0);
-
     iDof_ = dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL);
 
     momentum_ = dimensionedVector
@@ -430,17 +415,11 @@ template<class ParcelType>
 void Foam::DSMCCloud<ParcelType>::calculateFields()
 {
     scalarField& rhoN = rhoN_.internalField();
-
     scalarField& rhoM = rhoM_.internalField();
-
     scalarField& dsmcRhoN = dsmcRhoN_.internalField();
-
     scalarField& linearKE = linearKE_.internalField();
-
     scalarField& internalE = internalE_.internalField();
-
     scalarField& iDof = iDof_.internalField();
-
     vectorField& momentum = momentum_.internalField();
 
     forAllConstIter(typename DSMCCloud<ParcelType>, *this, iter)
@@ -449,17 +428,11 @@ void Foam::DSMCCloud<ParcelType>::calculateFields()
         const label cellI = p.cell();
 
         rhoN[cellI]++;
-
         rhoM[cellI] += constProps(p.typeId()).mass();
-
         dsmcRhoN[cellI]++;
-
         linearKE[cellI] += 0.5*constProps(p.typeId()).mass()*(p.U() & p.U());
-
         internalE[cellI] += p.Ei();
-
         iDof[cellI] += constProps(p.typeId()).internalDegreesOfFreedom();
-
         momentum[cellI] += constProps(p.typeId()).mass()*p.U();
     }
 
@@ -732,7 +705,6 @@ Foam::DSMCCloud<ParcelType>::DSMCCloud
     )
 {
     buildConstProps();
-
     buildCellOccupancy();
 
     // Initialise the collision selection remainder to a random value between 0
@@ -971,9 +943,7 @@ Foam::DSMCCloud<ParcelType>::DSMCCloud
     inflowBoundaryModel_()
 {
     clear();
-
     buildConstProps();
-
     initialise(dsmcInitialiseDict);
 }
 
@@ -1082,7 +1052,7 @@ template<class ParcelType>
 Foam::scalar Foam::DSMCCloud<ParcelType>::equipartitionInternalEnergy
 (
     scalar temperature,
-    scalar iDof
+    direction iDof
 )
 {
     scalar Ei = 0.0;
@@ -1099,17 +1069,13 @@ Foam::scalar Foam::DSMCCloud<ParcelType>::equipartitionInternalEnergy
     else
     {
         scalar a = 0.5*iDof - 1;
-
         scalar energyRatio;
-
         scalar P = -1;
 
         do
         {
             energyRatio = 10*rndGen_.scalar01();
-
             P = pow((energyRatio/a), a)*exp(a - energyRatio);
-
         } while (P < rndGen_.scalar01());
 
         Ei = energyRatio*physicoChemical::k.value()*temperature;

src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.H

@@ -350,7 +350,7 @@ public:
                 scalar equipartitionInternalEnergy
                 (
                     scalar temperature,
-                    scalar internalDegreesOfFreedom
+                    direction internalDegreesOfFreedom
                 );
 
 

src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H

@@ -89,7 +89,7 @@ public:
             scalar d_;
 
             //- Internal degrees of freedom
-            scalar internalDegreesOfFreedom_;
+            direction internalDegreesOfFreedom_;
 
             //- Viscosity index
             scalar omega_;
@@ -119,7 +119,7 @@ public:
             inline scalar sigmaT() const;
 
             //- Return the internalDegreesOfFreedom
-            inline scalar internalDegreesOfFreedom() const;
+            inline direction internalDegreesOfFreedom() const;
 
             //- Return the viscosity index
             inline scalar omega() const;

src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelI.H

@@ -45,7 +45,7 @@ inline Foam::DSMCParcel<ParcelType>::constantProperties::constantProperties
     d_(readScalar(dict.lookup("diameter"))),
     internalDegreesOfFreedom_
     (
-        readScalar(dict.lookup("internalDegreesOfFreedom"))
+        readInt(dict.lookup("internalDegreesOfFreedom"))
     ),
     omega_(readScalar(dict.lookup("omega")))
 {}
@@ -97,7 +97,7 @@ Foam::DSMCParcel<ParcelType>::constantProperties::sigmaT() const
 
 
 template<class ParcelType>
-inline Foam::scalar
+inline Foam::direction
 Foam::DSMCParcel<ParcelType>::constantProperties::internalDegreesOfFreedom()
 const
 {

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

@@ -42,7 +42,6 @@ Foam::scalar Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::energyRatio
     Random& rndGen(cloud.rndGen());
 
     scalar ChiAMinusOne = ChiA - 1;
-
     scalar ChiBMinusOne = ChiB - 1;
 
     if (ChiAMinusOne < SMALL && ChiBMinusOne < SMALL)
@@ -51,7 +50,6 @@ Foam::scalar Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::energyRatio
     }
 
     scalar energyRatio;
-
     scalar P;
 
     do
@@ -159,9 +157,7 @@ Foam::scalar Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::sigmaTcR
     }
 
     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
@@ -199,12 +195,10 @@ void Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::collide
     // DSMC0R.FOR
 
     scalar preCollisionEiP = EiP;
-
     scalar preCollisionEiQ = EiQ;
 
-    scalar iDofP = cloud.constProps(typeIdP).internalDegreesOfFreedom();
+    direction iDofP = cloud.constProps(typeIdP).internalDegreesOfFreedom();
-
+    direction iDofQ = cloud.constProps(typeIdQ).internalDegreesOfFreedom();
-    scalar iDofQ = cloud.constProps(typeIdQ).internalDegreesOfFreedom();
 
     scalar omegaPQ =
         0.5
@@ -214,17 +208,11 @@ void Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::collide
         );
 
     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)
@@ -233,9 +221,16 @@ void Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::collide
         {
             availableEnergy += preCollisionEiP;
 
+            if (iDofP == 2)
+            {
+                scalar energyRatio = 1.0 - pow(rndGen.scalar01(), (1.0/ChiB));
+                EiP = energyRatio*availableEnergy;
+            }
+            else
+            {
                 scalar ChiA = 0.5*iDofP;
-
                 EiP = energyRatio(ChiA, ChiB)*availableEnergy;
+            }
 
             availableEnergy -= EiP;
         }
@@ -247,10 +242,16 @@ void Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::collide
         {
             availableEnergy += preCollisionEiQ;
 
-            // Change to general LB ratio calculation
+            if (iDofQ == 2)
-            scalar energyRatio = 1.0 - pow(rndGen.scalar01(),(1.0/ChiB));
+            {
-
+                scalar energyRatio = 1.0 - pow(rndGen.scalar01(), (1.0/ChiB));
                 EiQ = energyRatio*availableEnergy;
+            }
+            else
+            {
+                scalar ChiA = 0.5*iDofQ;
+                EiQ = energyRatio(ChiA, ChiB)*availableEnergy;
+            }
 
             availableEnergy -= EiQ;
         }
@@ -260,11 +261,8 @@ void Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::collide
     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 = twoPi*rndGen.scalar01();
 
     vector postCollisionRelU =
@@ -277,7 +275,6 @@ void Foam::LarsenBorgnakkeVariableHardSphere<CloudType>::collide
         );
 
     UP = Ucm + postCollisionRelU*mQ/(mP + mQ);
-
     UQ = Ucm - postCollisionRelU*mP/(mP + mQ);
 }
 

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

@@ -108,7 +108,7 @@ void Foam::MaxwellianThermal<CloudType>::correct
 
     scalar mass = cloud.constProps(typeId).mass();
 
-    scalar iDof = cloud.constProps(typeId).internalDegreesOfFreedom();
+    direction iDof = cloud.constProps(typeId).internalDegreesOfFreedom();
 
     U =
         sqrt(physicoChemical::k.value()*T/mass)

src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C

@@ -110,7 +110,7 @@ void Foam::MixedDiffuseSpecular<CloudType>::correct
 
         scalar mass = cloud.constProps(typeId).mass();
 
-        scalar iDof = cloud.constProps(typeId).internalDegreesOfFreedom();
+        direction iDof = cloud.constProps(typeId).internalDegreesOfFreedom();
 
         U =
             sqrt(physicoChemical::k.value()*T/mass)