diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
index b6d20bd630f09482206239cc63e20e4a63f443f5..c588bb0cb76e30cafb874674d2f8508c6dbd5de0 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
@@ -555,6 +555,11 @@ void Foam::DsmcCloud<ParcelType>::evolve()
template<class ParcelType>
void Foam::DsmcCloud<ParcelType>::info() const
{
+ label nDsmcParticles = this->size();
+ reduce(nDsmcParticles, sumOp<label>());
+
+ scalar nMol = nDsmcParticles*nParticle_;
+
vector linearMomentum = linearMomentumOfSystem();
reduce(linearMomentum, sumOp<vector>());
@@ -565,20 +570,22 @@ void Foam::DsmcCloud<ParcelType>::info() const
reduce(internalEnergy, sumOp<scalar>());
Info<< "Cloud name: " << this->name() << nl
- << " Current number of parcels = "
- << returnReduce(this->size(), sumOp<label>()) << nl
- << " Current mass in system = "
+ << " Number of dsmc particles = "
+ << nDsmcParticles << nl
+ << " Number of molecules = "
+ << nMol << nl
+ << " Mass in system = "
<< returnReduce(massInSystem(), sumOp<scalar>()) << nl
- << " Linear momentum = "
- << linearMomentum << nl
- << " Linear momentum magnitude = "
- << mag(linearMomentum) << nl
- << " Linear kinetic energy = "
- << linearKineticEnergy << nl
- << " Internal energy = "
- << internalEnergy << nl
- << " Total energy = "
- << internalEnergy + linearKineticEnergy << nl
+ << " Average linear momentum = "
+ << linearMomentum/nMol << nl
+ << " |Average linear momentum| = "
+ << mag(linearMomentum)/nMol << nl
+ << " Average linear kinetic energy = "
+ << linearKineticEnergy/nMol << nl
+ << " Average internal energy = "
+ << internalEnergy/nMol << nl
+ << " Average total energy = "
+ << (internalEnergy + linearKineticEnergy)/nMol << nl
<< endl;
}
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
index 5d9edc9c06f2ce2448cb53bc1ce6f3377fa42f8b..d0bb26bc9ca81fe468784c833bcda165e1714ba8 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
@@ -97,7 +97,7 @@ class DsmcCloud
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
+ // Bird p220). Initialised with the parcels, updated as required, and
// read in on start/restart.
volScalarField sigmaTcRMax_;
diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
index ed55544077b937e09ac54c7a38bdc44a9ff91fd0..a69f74aa994d88459e395570adecbe65be66700e 100644
--- a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
+++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
@@ -126,8 +126,59 @@ void Foam::FreeStream<CloudType>::inflow()
{
particleFluxAccumulators_[i] +=
-patch.faceAreas() & (velocity_*numberDensities_[i]*deltaT);
+ }
+
+ forAll(patch, f)
+ {
+ // Loop over all faces as the outer loop to avoid calculating
+ // geometrical properties multiple times.
+
+ labelList faceVertices = patch[f];
+
+ label nVertices = faceVertices.size();
+
+ label globalFaceIndex = f + patch.start();
- forAll(particleFluxAccumulators_[i], f)
+ label cell = mesh.faceOwner()[globalFaceIndex];
+
+ const vector& fC = patch.faceCentres()[f];
+
+ scalar fA = mag(patch.faceAreas()[f]);
+
+ // Cummulative triangle area fractions
+ List<scalar> cTriAFracs(nVertices);
+
+ for (label v = 0; v < nVertices - 1; v++)
+ {
+ const point& vA = mesh.points()[faceVertices[v]];
+
+ const point& vB = mesh.points()[faceVertices[(v + 1)]];
+
+ cTriAFracs[v] =
+ 0.5*mag((vA - fC)^(vB - fC))/fA
+ + cTriAFracs[max((v - 1), 0)];
+ }
+
+ // Force the last area fraction value to 1.0 to avoid any
+ // rounding/non-flat face errors giving a value < 1.0
+ cTriAFracs[nVertices - 1] = 1.0;
+
+ // Normal unit vector *negative* so normal is pointing into the
+ // domain
+ vector nw = patch.faceAreas()[f];
+ nw /= -mag(nw);
+
+ // Wall tangential unit vector. Use the direction between the
+ // face centre and the first vertex in the list
+ vector tw1 = fC - (mesh.points()[faceVertices[0]]);
+ tw1 /= mag(tw1);
+
+ // Other tangential unit vector. Rescaling in case face is not
+ // flat and nw and tw1 aren't perfectly orthogonal
+ vector tw2 = nw ^ tw1;
+ tw2 /= mag(tw2);
+
+ forAll(particleFluxAccumulators_, i)
{
scalar& faceAccumulator = particleFluxAccumulators_[i][f];
@@ -140,39 +191,44 @@ void Foam::FreeStream<CloudType>::inflow()
scalar mass = cloud.constProps(typeId).mass();
- labelList faceVertices = patch[f];
+ for (label n = 0; n < nI; n++)
+ {
+ // Choose a triangle to insert on, based on their relative area
- label globalFaceIndex = f + patch.start();
+ scalar triSelection = rndGen.scalar01();
- label cell = mesh.faceOwner()[globalFaceIndex];
+ // Selected triangle
+ label sTri = -1;
- const vector& fC = patch.faceCentres()[f];
+ forAll(cTriAFracs, tri)
+ {
+ sTri = tri;
- for (label n = 0; n < nI; n++)
- {
- // Temporarily insert particles half way between the face and
- // cell centres
- vector p = 0.5*(fC + mesh.cellCentres()[cell]);
+ if (cTriAFracs[tri] >= triSelection)
+ {
+ break;
+ }
+ }
- // Normal unit vector *negative* so normal is pointing into the
- // domain
- vector nw = patch.faceAreas()[f];
- nw /= -mag(nw);
+ // Randomly distribute the points on the triangle, using the
+ // method from:
+ // Generating Random Points in Triangles
+ // by Greg Turk
+ // from "Graphics Gems", Academic Press, 1990
+ // http://tog.acm.org/GraphicsGems/gems/TriPoints.c
- // Wall tangential unit vector. Use the direction between the
- // face centre and the first vertex in the list
- vector tw1 = fC - (mesh.points()[faceVertices[0]]);
- tw1 /= mag(tw1);
+ const point& A = fC;
+ const point& B = mesh.points()[faceVertices[sTri]];
+ const point& C =
+ mesh.points()[faceVertices[(sTri + 1) % nVertices]];
- // Other tangential unit vector. Rescaling in case face is not
- // flat and nw and tw1 aren't perfectly orthogonal
- vector tw2 = nw ^ tw1;
- tw2 /= mag(tw2);
+ scalar s = rndGen.scalar01();
+ scalar t = sqrt(rndGen.scalar01());
- scalar C = sqrt(CloudType::kb*temperature_/mass);
+ point p = (1 - t)*A + (1 - s)*t*B + s*t*C;
vector U =
- C
+ sqrt(CloudType::kb*temperature_/mass)
*(
rndGen.GaussNormal()*tw1
+ rndGen.GaussNormal()*tw2
@@ -204,14 +260,6 @@ void Foam::FreeStream<CloudType>::inflow()
Info<< " Particles inserted = "
<< particlesInserted << endl;
- // Info<< "insert particles now! " << nl
- // << temperature_ << nl
- // << velocity_ << nl
- // << moleculeTypeIds_ << nl
- // << numberDensities_ << nl
- // << particleFluxAccumulators_ << nl
- // << patch
- // << endl;
}
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
index 23cd0bed34613953eca7b34f55296559bc100af0..a6ed8e5718381c05491a1cf7c1aa859fa1472753 100644
--- a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
@@ -102,10 +102,8 @@ void Foam::MaxwellianThermal<CloudType>::correct
scalar T = cloud.T().boundaryField()[wppIndex][wppLocalFace];
- scalar C = sqrt(CloudType::kb*T/mass);
-
U =
- C
+ sqrt(CloudType::kb*T/mass)
*(
rndGen.GaussNormal()*tw1
+ rndGen.GaussNormal()*tw2