diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
index 67f549fc578e4fc18952861e64b815003da2ffcf..c84b0b7205b6156daf943742083bd19f4bfaab26 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
@@ -33,6 +33,8 @@ License
template<class ParcelType>
Foam::scalar Foam::DsmcCloud<ParcelType>::kb = 1.380650277e-23;
+template<class ParcelType>
+Foam::scalar Foam::DsmcCloud<ParcelType>::Tref = 273;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@@ -185,13 +187,113 @@ void Foam::DsmcCloud<ParcelType>::initialise
}
}
}
+
+ // Initialise the sigmaTcRMax_ field to the product of the cross section of
+ // the most abundant species and the most probable thermal speed (Bird,
+ // p222-223)
+
+ label mostAbundantType(findMax(numberDensities));
+
+ const typename ParcelType::constantProperties& cP = constProps
+ (
+ mostAbundantType
+ );
+
+ sigmaTcRMax_.internalField() = cP.sigmaT()*maxwellianMostProbableSpeed
+ (
+ temperature,
+ cP.mass()
+ );
}
template<class ParcelType>
void Foam::DsmcCloud<ParcelType>::collisions()
{
- Info<< "DsmcCloud collisions() - PLACEHOLDER" << endl;
+ scalar deltaT = mesh_.time().deltaT().value();
+
+ label collisionCandidates = 0;
+
+ label collisions = 0;
+
+ forAll(cellOccupancy_, celli)
+ {
+ const DynamicList<ParcelType*>& cellParcels(cellOccupancy_[celli]);
+
+ label nC(cellParcels.size());
+
+ if (nC > 1)
+ {
+ scalar sigmaTcRMax = sigmaTcRMax_[celli];
+
+ scalar selectedPairs = collisionSelectionRemainder_[celli]
+ + 0.5*nC*(nC-1)*nParticle_*sigmaTcRMax*deltaT
+ /mesh_.cellVolumes()[celli];
+
+ label nCandidates(selectedPairs);
+
+ collisionSelectionRemainder_[celli] = selectedPairs - nCandidates;
+
+ collisionCandidates += nCandidates;
+
+ for(label c = 0; c < nCandidates; c++)
+ {
+ // Select the first collision candidate
+
+ label candidateP = rndGen_.integer(0, nC-1);
+
+ label candidateQ = rndGen_.integer(0, nC-1);
+
+ // If the same candidate is chosen, choose again
+ while(candidateP == candidateQ)
+ {
+ candidateQ = rndGen_.integer(0, nC-1);
+ }
+
+ ParcelType* parcelP = cellParcels[candidateP];
+ ParcelType* parcelQ = cellParcels[candidateQ];
+
+ label typeIdP = parcelP->typeId();
+ label typeIdQ = parcelQ->typeId();
+
+ scalar sigmaTcR = binaryCollision().sigmaTcR
+ (
+ typeIdP,
+ typeIdQ,
+ parcelP->U(),
+ parcelQ->U()
+ );
+
+ // Update the maximum value of sigmaTcR stored, but use the
+ // initial value in the acceptance-rejection criteria because
+ // the number of collision candidates selected was based on this
+
+ if (sigmaTcR > sigmaTcRMax_[celli])
+ {
+ sigmaTcRMax_[celli] = sigmaTcR;
+ }
+
+ if ((sigmaTcR/sigmaTcRMax) > rndGen_.scalar01())
+ {
+ binaryCollision().collide
+ (
+ typeIdP,
+ typeIdQ,
+ parcelP->U(),
+ parcelQ->U()
+ );
+
+ collisions++;
+ }
+ }
+ }
+ }
+
+ Info<< " Collisions = "
+ << collisions << nl
+ << " Acceptance rate = "
+ << scalar(collisions)/scalar(collisionCandidates) << nl
+ << endl;
}
@@ -246,6 +348,19 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
typeIdList_(particleProperties_.lookup("typeIdList")),
nParticle_(readScalar(particleProperties_.lookup("nEquivalentParticles"))),
cellOccupancy_(mesh_.nCells()),
+ sigmaTcRMax_
+ (
+ IOobject
+ (
+ this->name() + "SigmaTcRMax",
+ mesh_.time().timeName(),
+ mesh_,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh_
+ ),
+ collisionSelectionRemainder_(mesh_.nCells(), 0),
constProps_(),
rndGen_(label(971501)),
binaryCollisionModel_
@@ -297,6 +412,20 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
typeIdList_(particleProperties_.lookup("typeIdList")),
nParticle_(readScalar(particleProperties_.lookup("nEquivalentParticles"))),
cellOccupancy_(),
+ sigmaTcRMax_
+ (
+ IOobject
+ (
+ this->name() + "SigmaTcRMax",
+ mesh_.time().timeName(),
+ mesh_,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh_,
+ dimensionedScalar("zero", dimensionSet(0, 3, -1, 0, 0), 0.0)
+ ),
+ collisionSelectionRemainder_(),
constProps_(),
rndGen_(label(971501)),
binaryCollisionModel_(),
@@ -322,13 +451,15 @@ Foam::DsmcCloud<ParcelType>::~DsmcCloud()
template<class ParcelType>
void Foam::DsmcCloud<ParcelType>::evolve()
{
+ buildCellOccupancy();
+
typename ParcelType::trackData td(*this);
//this->injection().inject(td);
if (debug)
{
- this->dumpParticlePositions();
+ this->dumpParticlePositions();
}
// Move the particles ballistically with their current velocities
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
index 364d8ee518f0d5ef344683e0476655c985bbf702..0cbd4aecba89c303a918cfad49a3cd4f8b770289 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
@@ -44,6 +44,7 @@ SourceFiles
#include "Random.H"
#include "fvMesh.H"
#include "volFields.H"
+#include "scalarIOField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -95,6 +96,14 @@ private:
//- A data structure holding which particles are in which cell
List<DynamicList<ParcelType*> > cellOccupancy_;
+ //- An IOField holding the value of (sigmaT * cR)max for each cell (see
+ // Bird p220). Initialised with the parcelsm updated as required, and
+ // read in on start/restart.
+ volScalarField sigmaTcRMax_;
+
+ //- A field holding the remainder from the previous collision selections
+ scalarField collisionSelectionRemainder_;
+
//- Parcel constant properties - one for each type
List<typename ParcelType::constantProperties> constProps_;
@@ -137,8 +146,12 @@ public:
// Static data members
+ //- Boltzmann constant
static scalar kb;
+ //- Reference temperature for all models
+ static scalar Tref;
+
// Constructors
@@ -181,13 +194,21 @@ public:
inline const List<word>& typeIdList() const;
//- Return the number of real particles represented by one
- //- parcel
+ // parcel
inline label nParticle() const;
//- Return the cell occupancy addressing
inline const List<DynamicList<ParcelType*> >&
cellOccupancy() const;
+ //- Return the sigmaTcRMax field. non-const access to allow
+ // updating.
+ inline volScalarField& sigmaTcRMax();
+
+ //- Return the collision selection remainder field. non-const
+ // access to allow updating.
+ inline scalarField& collisionSelectionRemainder();
+
//- Return all of the constant properties
inline const List<typename ParcelType::constantProperties>&
constProps() const;
@@ -199,17 +220,52 @@ public:
//- Return refernce to the random object
inline Random& rndGen();
+ // Kinetic theory helper functions
+
+ //- Generate a random velocity sampled from the Maxwellian speed
+ // distribution
+ vector equipartitionLinearVelocity
+ (
+ scalar temperature,
+ scalar mass
+ );
+
+ //- From the Maxwellian distribution:
+
+ //- Average particle speed
+ inline scalar maxwellianAverageSpeed
+ (
+ scalar temperature,
+ scalar mass
+ ) const;
+
+ //- RMS particle speed
+ inline scalar maxwellianRMSSpeed
+ (
+ scalar temperature,
+ scalar mass
+ ) const;
+
+ //- Most probable speed
+ inline scalar maxwellianMostProbableSpeed
+ (
+ scalar temperature,
+ scalar mass
+ ) const;
+
//- Clear the Cloud
inline void clear();
-
// Sub-models
//- Return reference to binary elastic collision model
- inline const
- BinaryCollisionModel<DsmcCloud<ParcelType> >&
+ inline const BinaryCollisionModel<DsmcCloud<ParcelType> >&
binaryCollision() const;
+ //- Return non-const reference to binary elastic collision model
+ inline BinaryCollisionModel<DsmcCloud<ParcelType> >&
+ binaryCollision();
+
//- Return reference to wall interaction model
inline const WallInteractionModel<DsmcCloud<ParcelType> >&
wallInteraction() const;
@@ -229,14 +285,6 @@ public:
//- Print cloud information
void info() const;
- //- Generate a random velocity sampled from the Maxwellian speed
- // distribution
- vector equipartitionLinearVelocity
- (
- scalar temperature,
- scalar mass
- );
-
//- Dump particle positions to .obj file
void dumpParticlePositions() const;
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
index 46320f79d5a5f9e2deb1749b9c25c943332a22ed..91d8f9488b3f19c95aa9d88e5b69e2c07ee90bec 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
@@ -72,6 +72,21 @@ Foam::DsmcCloud<ParcelType>::cellOccupancy() const
}
+template<class ParcelType>
+inline Foam::volScalarField& Foam::DsmcCloud<ParcelType>::sigmaTcRMax()
+{
+ return sigmaTcRMax_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalarField&
+Foam::DsmcCloud<ParcelType>::collisionSelectionRemainder()
+{
+ return collisionSelectionRemainder_;
+}
+
+
template<class ParcelType>
inline const Foam::List<typename ParcelType::constantProperties>&
Foam::DsmcCloud<ParcelType>::constProps() const
@@ -114,6 +129,14 @@ Foam::DsmcCloud<ParcelType>::binaryCollision() const
}
+template<class ParcelType>
+inline Foam::BinaryCollisionModel<Foam::DsmcCloud<ParcelType> >&
+Foam::DsmcCloud<ParcelType>::binaryCollision()
+{
+ return binaryCollisionModel_();
+}
+
+
template<class ParcelType>
inline const Foam::WallInteractionModel<Foam::DsmcCloud<ParcelType> >&
Foam::DsmcCloud<ParcelType>::wallInteraction() const
@@ -135,6 +158,8 @@ inline Foam::scalar Foam::DsmcCloud<ParcelType>::massInSystem() const
{
scalar sysMass = 0.0;
+ Info<< "massInSystem() not implemented" << endl;
+
return sysMass;
}
@@ -146,6 +171,40 @@ inline Foam::Random& Foam::DsmcCloud<ParcelType>::rndGen()
}
+template<class ParcelType>
+inline Foam::scalar Foam::DsmcCloud<ParcelType>::maxwellianAverageSpeed
+(
+ scalar temperature,
+ scalar mass
+) const
+{
+ return 2.0*sqrt(2.0*kb*temperature/(mathematicalConstant::pi*mass));
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::DsmcCloud<ParcelType>::maxwellianRMSSpeed
+(
+ scalar temperature,
+ scalar mass
+) const
+{
+ return sqrt(3.0*kb*temperature/mass);
+}
+
+
+template<class ParcelType>
+inline Foam::scalar
+Foam::DsmcCloud<ParcelType>::maxwellianMostProbableSpeed
+(
+ scalar temperature,
+ scalar mass
+) const
+{
+ return sqrt(2.0*kb*temperature/mass);
+}
+
+
template<class ParcelType>
inline const Foam::tmp<Foam::volScalarField>
Foam::DsmcCloud<ParcelType>::rhoN() const
diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H
index 599592f055092b75122858725e171f9da5065df2..5ea8d57624c1c7dfbcc7bcab116f734d0f7c542a 100644
--- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H
+++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H
@@ -87,6 +87,11 @@ public:
//- Particle hard sphere diameter [m] (constant)
scalar d_;
+ //- Internal degrees of freedom
+ scalar internalDegreesOfFreedom_;
+
+ //- Viscosity index
+ scalar omega_;
public:
@@ -104,6 +109,16 @@ public:
//- Return const access to the minimum particle mass
inline scalar d() const;
+
+ //- Return the reference total collision cross section
+ inline scalar sigmaT() const;
+
+ //- Return the internalDegreesOfFreedom
+ inline scalar internalDegreesOfFreedom() const;
+
+ //- Return the viscosity index
+ inline scalar omega() const;
+
};
@@ -195,9 +210,6 @@ public:
// Edit
- //- Return access to diameter
- inline scalar& d();
-
//- Return access to velocity
inline vector& U();
diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelI.H b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelI.H
index 352ea27ca1830bdadf06badc55096b8da65128ad..e4c117b2c6335c4adb4944ae20f8924867df0663 100644
--- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelI.H
+++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelI.H
@@ -41,7 +41,12 @@ inline Foam::DsmcParcel<ParcelType>::constantProperties::constantProperties
)
:
mass_(readScalar(dict.lookup("mass"))),
- d_(readScalar(dict.lookup("diameter")))
+ d_(readScalar(dict.lookup("diameter"))),
+ internalDegreesOfFreedom_
+ (
+ readScalar(dict.lookup("internalDegreesOfFreedom"))
+ ),
+ omega_(readScalar(dict.lookup("omega")))
{}
@@ -90,6 +95,31 @@ Foam::DsmcParcel<ParcelType>::constantProperties::d() const
}
+template <class ParcelType>
+inline Foam::scalar
+Foam::DsmcParcel<ParcelType>::constantProperties::sigmaT() const
+{
+ return mathematicalConstant::pi*d_*d_;
+}
+
+
+template <class ParcelType>
+inline Foam::scalar
+Foam::DsmcParcel<ParcelType>::constantProperties::internalDegreesOfFreedom()
+const
+{
+ return internalDegreesOfFreedom_;
+}
+
+
+template <class ParcelType>
+inline Foam::scalar
+Foam::DsmcParcel<ParcelType>::constantProperties::omega() const
+{
+ return omega_;
+}
+
+
// * * * * * * * * * * * trackData Member Functions * * * * * * * * * * * * //
template <class ParcelType>
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.C b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.C
index 9c4c306e941d96d1c624af4aaac00151bfd87896..04ba113e3f3a6cd55e99c9d8f92975c484a8abde 100644
--- a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.C
+++ b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.C
@@ -58,6 +58,14 @@ Foam::BinaryCollisionModel<CloudType>::owner() const
}
+template<class CloudType>
+CloudType&
+Foam::BinaryCollisionModel<CloudType>::owner()
+{
+ return owner_;
+}
+
+
template<class CloudType>
const Foam::dictionary&
Foam::BinaryCollisionModel<CloudType>::dict() const
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.H b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.H
index f8ffbe93ae085b99f7baf7b37170d233f96cf4af..5ae8ace6797ed66ff61ecc3be6ef4da5b84f3d43 100644
--- a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.H
+++ b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.H
@@ -115,6 +115,9 @@ public:
//- Return the owner cloud object
const CloudType& owner() const;
+ //- Return non-const access to the owner cloud object
+ CloudType& owner();
+
//- Return the dictionary
const dictionary& dict() const;
@@ -124,14 +127,23 @@ public:
// Member Functions
+ //- Return the collision cross section * relative velocity product
+ virtual scalar sigmaTcR
+ (
+ label typeIdP,
+ label typeIdQ,
+ const vector& UP,
+ const vector& UQ
+ ) const = 0;
+
//- Apply collision
virtual void collide
(
- label idA,
- label idB,
- vector& UA,
- vector& UB
- ) const = 0;
+ label typeIdP,
+ label typeIdQ,
+ vector& UP,
+ vector& UQ
+ ) = 0;
};
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C
index 62bbeba0865480f6db2f32805d0bb7b0d49df795..58cabefccbcfe1e71b564a37db772ab27cdf7e55 100644
--- a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C
+++ b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C
@@ -50,15 +50,80 @@ Foam::VariableHardSphere<CloudType>::~VariableHardSphere()
template <class CloudType>
-void Foam::VariableHardSphere<CloudType>::collide
+Foam::scalar Foam::VariableHardSphere<CloudType>::sigmaTcR
(
- label idA,
- label idB,
- vector& UA,
- vector& UB
+ 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::VariableHardSphere<CloudType>::collide
+(
+ label typeIdP,
+ label typeIdQ,
+ vector& UP,
+ vector& UQ
+)
+{
+ CloudType& cloud(this->owner());
+
+ Random& rndGen(cloud.rndGen());
+
+ scalar mP = cloud.constProps(typeIdP).mass();
+
+ scalar mQ = cloud.constProps(typeIdQ).mass();
+
+ vector Ucm = (mP*UP + mQ*UQ)/(mP + mQ);
+
+ scalar cR = mag(UP - UQ);
+
+ 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);
}
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.H b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.H
index bacdda6c55ad0d3674da6abd751aad2f34901461..00b7ec9527ee18ed683a7c455c1ce98966a0ea45 100644
--- a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.H
+++ b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.H
@@ -75,14 +75,23 @@ public:
// 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 idA,
- label idB,
- vector& UA,
- vector& UB
- ) const;
+ label typeIdP,
+ label typeIdQ,
+ vector& UP,
+ vector& UQ
+ );
};
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
index 050397d9370b376ceec33b94a21d2a3c439d319a..ab6939b0457723644e2bcb11010114c7e5573d41 100644
--- a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
@@ -78,7 +78,7 @@ void Foam::MaxwellianThermal<CloudType>::correct
// Wall tangential velocity (flow direction)
vector Ut = U - magUn*nw;
- CloudType& cloud(WallInteractionModel<CloudType>::owner());
+ CloudType& cloud(this->owner());
Random& rndGen(cloud.rndGen());