diff --git a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
index ff6166ffd020eba945bd4da3f1af8de345ce82ec..aafa52ca59affe2779e4a0acaa3063895261c23b 100644
--- a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
+++ b/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;
diff --git a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.H b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.H
index eadbbd85ef315d7ba9ef57568584ee728cabef9b..6b57ae1ed7decbc068273b7e9a84ef1ae3c85cab 100644
--- a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.H
+++ b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.H
@@ -350,7 +350,7 @@ public:
scalar equipartitionInternalEnergy
(
scalar temperature,
- scalar internalDegreesOfFreedom
+ direction internalDegreesOfFreedom
);
diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H
index 3d9106c7ef194de2e145d972c71ed654b0e8435a..10d93b50ddbb48a4ab15ad5145dceef44982acc7 100644
--- a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.H
+++ b/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;
diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelI.H b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelI.H
index 50a23c908934fa422744f349ab9fdb31d9e76313..54de74c96de031c42210bca43e06207a98f367fd 100644
--- a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcelI.H
+++ b/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
{
diff --git a/src/lagrangian/DSMC/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C b/src/lagrangian/DSMC/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C
index 70979c901aa27d4e0606899c46b9369dcd3fa39e..9139e09be3b0f15f21f5ccb914a4d6e04bb4a92e 100644
--- a/src/lagrangian/DSMC/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C
+++ b/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();
-
- scalar iDofQ = cloud.constProps(typeIdQ).internalDegreesOfFreedom();
+ direction iDofP = cloud.constProps(typeIdP).internalDegreesOfFreedom();
+ direction 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;
- scalar ChiA = 0.5*iDofP;
-
- EiP = energyRatio(ChiA, ChiB)*availableEnergy;
+ 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
- scalar energyRatio = 1.0 - pow(rndGen.scalar01(),(1.0/ChiB));
-
- EiQ = energyRatio*availableEnergy;
+ if (iDofQ == 2)
+ {
+ 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);
}
diff --git a/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C b/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
index 9c976fc50b806482629960b42973f36cc8829390..3bc25401edb2fa2faf9ff4dc538c0bd30b46c284 100644
--- a/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
+++ b/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)
diff --git a/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C b/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C
index c4ef48c0c9fabc2f2ba9de74910f38f4e3398f9c..06e5a7079c87b7f6f811e526515de1005b0e6ede 100644
--- a/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C
+++ b/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)