diff --git a/applications/solvers/dsmc/dsmcFoam/createFields.H b/applications/solvers/dsmc/dsmcFoam/createFields.H
index 68a44f03fd3e14c76ad7f4c3abddd098155c1c4f..5f586665bbac29f6dc2d8ce19976267fa4679faf 100644
--- a/applications/solvers/dsmc/dsmcFoam/createFields.H
+++ b/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);
diff --git a/applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C b/applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C
index 34727c0bfffdb011c9f0f86a1e6c89b0e07a5b03..04c1914f3ff4d6b7d7df4300901fc747450a29ad 100644
--- a/applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C
+++ b/applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C
@@ -26,8 +26,8 @@ Application
dsmcFoam
Description
- Initialise a case for dsmcFoam by reading the initialise subdictionary in
- dsmcProperties
+ Initialise a case for dsmcFoam by reading the initialisation dictionary
+ 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();
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
index 00164e4b5181966b31113b4eb5303a0509c8e653..1867670f70bba1d4e780d38ec8f58ab99801effc 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
+++ b/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
- (
- dsmcInitialiseDict.lookup("numberDensities")
- );
+ List<word> molecules(numberDensitiesDict.toc());
+
+ Field<scalar> numberDensities(molecules.size());
- if(molecules.size() != numberDensities.size())
+ forAll(molecules, i)
{
- FatalErrorIn("Foam::Foam::DsmcCloud<ParcelType>::initialise")
- << "molecules and numberDensities must be the same size."
- << nl << abort(FatalError);
+ numberDensities[i] = readScalar
+ (
+ 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 = "
- << collisions << nl
- << " Acceptance rate = "
- << scalar(collisions)/scalar(collisionCandidates) << nl
- << endl;
+ if (collisionCandidates)
+ {
+ Info<< " Collisions = "
+ << collisions << nl
+ << " 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,
+ mesh_.time().constant(),
+ 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 IOdictionary& dsmcInitialiseDict
+ const fvMesh& mesh
)
:
Cloud<ParcelType>(mesh, cloudType, false),
@@ -409,8 +427,8 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
IOobject
(
cloudType + "Properties",
- mesh.time().constant(),
- mesh,
+ mesh_.time().constant(),
+ 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;
}
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
index 1d071afa5bbe0609f7c34be33f3ea985bc87e56d..8ea6f31fe1207815aebcf69916fbf99efda8d928 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
+++ b/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
- // dsmcInitialse.
+ //- Construct given name and mesh. Used to initialise.
DsmcCloud
(
const word& cloudType,
- const fvMesh& mesh,
- const IOdictionary& dsmcInitialiseDict
+ const fvMesh& mesh
);
@@ -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();
+
+
};
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
index 3cbc44e98b9df2c396b3f761f55c112bc9869222..057570030e6b1c8c49ecb93b42d32c2dd0404635 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
+++ b/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>
-Foam::DsmcCloud<ParcelType>::U() const
-{
- tmp<volScalarField> tU
- (
- 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
+inline void Foam::DsmcCloud<ParcelType>::clear()
{
- 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;
+ return IDLList<ParcelType>::clear();
}
-// 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;
-// }
-
// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C
index 38a7c827f416a7c37973d43d11d973aca625e3fd..000ffcd53efb82dc4454009b3a28dc9c715cf124 100644
--- a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C
+++ b/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 IOdictionary& dsmcInitialiseDict
+ const fvMesh& mesh
)
:
- DsmcCloud<dsmcParcel>(cloudType, mesh, dsmcInitialiseDict)
+ DsmcCloud<dsmcParcel>(cloudType, mesh)
{}
diff --git a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H
index 7e0ee0cf8423e438428d191a7c1eab6db3c1b443..2f88bbaad563dbb99fef716f22256416ae98bc15 100644
--- a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H
+++ b/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 IOdictionary& dsmcInitialiseDict
+ const fvMesh& mesh
);
//- Destructor
diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C
index 941aaacdc64dca8dac4fb877ce5d56f669eba6d2..1b38143c90fc0e144b6f023a79975e06d4f6164b 100644
--- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C
+++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C
@@ -65,8 +65,7 @@ bool Foam::DsmcParcel<ParcelType>::move
{
if
(
- isType<processorPolyPatch>
- (pbMesh[p.patch(p.face())])
+ isType<processorPolyPatch>(pbMesh[p.patch(p.face())])
)
{
td.switchProcessor = true;
diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H
index 91f17669dd5cc2071cd05f863b2d7a99947f23fa..d8a6664a3afd8678710a54efef4fa64656f3a00a 100644
--- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H
+++ b/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
diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelI.H b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelI.H
index e4c117b2c6335c4adb4944ae20f8924867df0663..ff86419186fee680205fe2b70e0d0d2a6c705e61 100644
--- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelI.H
+++ b/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_;
+}
+
+
// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelIO.C b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelIO.C
index eb49da999696e4d0f6f960fce68ae8815ecca872..b8bc688c09640bfd462f754fc87175a46334c9ff 100644
--- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelIO.C
+++ b/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())
);
}
diff --git a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/dsmcParcel.C b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/dsmcParcel.C
index a5efcbca6a10c451d94882abe2bdb91172525f52..8f52fcdabf62a8dbd9e7cba06295a3ea41cbbfba 100644
--- a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/dsmcParcel.C
+++ b/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
)
diff --git a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/dsmcParcel.H b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/dsmcParcel.H
index 684af80b6a6b760b8cce3dcf14b52855ffa07100..24d9c9ada4f4654f8911c151e78b806d782807eb 100644
--- a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/dsmcParcel.H
+++ b/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
);
diff --git a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelBinaryCollisionModels.C b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelBinaryCollisionModels.C
index ea170786005337a17e8e53068e43a2e640ecffc1..de2444687b785586345a50d8eaf922bd63073b10 100644
--- a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelBinaryCollisionModels.C
+++ b/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
+ );
};
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.H b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.H
index cc279b51b70a5f25c1b831cee679393dd28ece58..771115c98428b86253f8d4920b9ad577067b7bd5 100644
--- a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.H
+++ b/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;
};
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C
new file mode 100644
index 0000000000000000000000000000000000000000..7837396e9d8f6f63ec834203f8067a5c800c9daf
--- /dev/null
+++ b/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);
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.H b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.H
new file mode 100644
index 0000000000000000000000000000000000000000..1e0f8bd9755d59781aa8466590dbd4185971f589
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C
index 433ecd24cf06a1daac4b3fa47c408afda5697cdb..e740bee537931bdc797d57b05354de59e2a1fac2 100644
--- a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.C
+++ b/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());
diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.H b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.H
index 22550288b4e95be7c8dee39058d012f38a24be2b..a1dbdb9a270c0e7c8dcb3021601fc88355e035c2 100644
--- a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/VariableHardSphere/VariableHardSphere.H
+++ b/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
);
};