diff --git a/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C b/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C
index 1ea52cc6aae883bd4f5fa85e93b0294ca799d638..cb068960e8e459dc3c3fc7db1d45e458f649c552 100644
--- a/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C
+++ b/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C
@@ -91,7 +91,7 @@ int main(int argc, char *argv[])
runTime.write();
- Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+ Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
diff --git a/applications/utilities/postProcessing/sampling/sample/sample.C b/applications/utilities/postProcessing/sampling/sample/sample.C
index 2451522ffe5d2a9f08e72e3fe882eebaade5bb82..12864fdceb7daf69c129b1e4b9eba68d04f46005 100644
--- a/applications/utilities/postProcessing/sampling/sample/sample.C
+++ b/applications/utilities/postProcessing/sampling/sample/sample.C
@@ -97,10 +97,11 @@ using namespace Foam;
int main(int argc, char *argv[])
{
timeSelector::addOptions();
+# include "addRegionOption.H"
# include "setRootCase.H"
# include "createTime.H"
instantList timeDirs = timeSelector::select0(runTime, args);
-# include "createMesh.H"
+# include "createNamedMesh.H"
IOsampledSets sSets
(
diff --git a/applications/utilities/postProcessing/sampling/sample/sampleDict b/applications/utilities/postProcessing/sampling/sample/sampleDict
index 559eb2843678da0a0ea06eec30c154c6a82957b7..4284b2878c3ae308b563648976ce56fb28eed41c 100644
--- a/applications/utilities/postProcessing/sampling/sample/sampleDict
+++ b/applications/utilities/postProcessing/sampling/sample/sampleDict
@@ -105,6 +105,14 @@ sets
end (2 0.51 0.005);
nPoints 10;
}
+
+ somePoints
+ {
+ type cloud;
+ axis xyz;
+ points ((0.049 0.049 0.005)(0.051 0.049 0.005));
+ }
+
);
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C
index 0c3bf4340838d01e7dbbf264a9fb6fdb948ed6ce..d85e393fd8d430687f8c0d1eaaabac5d5c586ab6 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C
@@ -296,6 +296,16 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
this->operator==(newValues);
+ if (debug)
+ {
+ Info<< "directMapped on field:" << fldName
+ << " patch:" << this->patch().name()
+ << " avg:" << gAverage(*this)
+ << " min:" << gMin(*this)
+ << " max:" << gMax(*this)
+ << endl;
+ }
+
fixedValueFvPatchField<Type>::updateCoeffs();
}
diff --git a/src/lagrangian/basic/Particle/Particle.C b/src/lagrangian/basic/Particle/Particle.C
index ec1d09a00f94244a5183af5a09fe41aae70f216c..6c8cbe4aa106f7181d0f93a84a8e6c5c5cd2996c 100644
--- a/src/lagrangian/basic/Particle/Particle.C
+++ b/src/lagrangian/basic/Particle/Particle.C
@@ -390,7 +390,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
// slightly towards the cell-centre.
if (trackFraction < SMALL)
{
- position_ += 1.0e-6*(mesh.cellCentres()[celli_] - position_);
+ position_ += 1.0e-3*(mesh.cellCentres()[celli_] - position_);
}
return trackFraction;
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
index 4c4bcd0a43dc965a3bd1daae9b06801b0d359c77..53177e9afebcc5a263917230543768ee48c5ebe1 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
@@ -37,6 +37,7 @@ Foam::scalar Foam::DsmcCloud<ParcelType>::kb = 1.380650277e-23;
template<class ParcelType>
Foam::scalar Foam::DsmcCloud<ParcelType>::Tref = 273;
+
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class ParcelType>
@@ -113,87 +114,161 @@ void Foam::DsmcCloud<ParcelType>::initialise
numberDensities /= nParticle_;
- scalar x0 = mesh_.bounds().min().x();
- scalar xR = mesh_.bounds().max().x() - x0;
+ forAll(mesh_.cells(), cell)
+ {
+ const vector& cC = mesh_.cellCentres()[cell];
+ const labelList& cellFaces = mesh_.cells()[cell];
+ const scalar cV = mesh_.cellVolumes()[cell];
+
+ label nTets = 0;
- scalar y0 = mesh_.bounds().min().y();
- scalar yR = mesh_.bounds().max().y() - y0;
+ // Each face is split into nEdges (or nVertices) - 2 tets.
+ forAll(cellFaces, face)
+ {
+ nTets += mesh_.faces()[cellFaces[face]].size() - 2;
+ }
- scalar z0 = mesh_.bounds().min().z();
- scalar zR = mesh_.bounds().max().z() - z0;
+ // Calculate the cumulative tet volumes circulating around the cell and
+ // record the vertex labels of each.
+ scalarList cTetVFracs(nTets, 0.0);
- forAll(molecules, i)
- {
- const word& moleculeName(molecules[i]);
+ List<labelList> tetPtIs(nTets, labelList(3,-1));
- label typeId(findIndex(typeIdList_, moleculeName));
+ // Keep track of which tet this is.
+ label tet = 0;
- if (typeId == -1)
+ forAll(cellFaces, face)
{
- FatalErrorIn("Foam::DsmcCloud<ParcelType>::initialise")
- << "typeId " << moleculeName << "not defined." << nl
- << abort(FatalError);
- }
+ const labelList& facePoints = mesh_.faces()[cellFaces[face]];
- const typename ParcelType::constantProperties& cP = constProps(typeId);
+ label pointI = 1;
+ while ((pointI + 1) < facePoints.size())
+ {
- scalar numberDensity = numberDensities[i];
+ const vector& pA = mesh_.points()[facePoints[0]];
+ const vector& pB = mesh_.points()[facePoints[pointI]];
+ const vector& pC = mesh_.points()[facePoints[pointI + 1]];
- scalar spacing = pow(numberDensity,-(1.0/3.0));
+ cTetVFracs[tet] =
+ mag(((pA - cC) ^ (pB - cC)) & (pC - cC))/(cV*6.0)
+ + cTetVFracs[max((tet - 1),0)];
- int ni = label(xR/spacing) + 1;
- int nj = label(yR/spacing) + 1;
- int nk = label(zR/spacing) + 1;
+ tetPtIs[tet][0] = facePoints[0];
+ tetPtIs[tet][1] = facePoints[pointI];
+ tetPtIs[tet][2] = facePoints[pointI + 1];
- vector delta(xR/ni, yR/nj, zR/nk);
+ pointI++;
+ tet++;
+ }
+ }
- scalar pert = spacing;
+ // Force the last volume fraction value to 1.0 to avoid any
+ // rounding/non-flat face errors giving a value < 1.0
+ cTetVFracs[nTets - 1] = 1.0;
- for (int i = 0; i < ni; i++)
+ forAll(molecules, i)
{
- for (int j = 0; j < nj; j++)
+ const word& moleculeName(molecules[i]);
+
+ label typeId(findIndex(typeIdList_, moleculeName));
+
+ if (typeId == -1)
{
- for (int k = 0; k < nk; k++)
- {
- point p
- (
- x0 + (i + 0.5)*delta.x(),
- y0 + (j + 0.5)*delta.y(),
- z0 + (k + 0.5)*delta.z()
- );
+ FatalErrorIn("Foam::DsmcCloud<ParcelType>::initialise")
+ << "typeId " << moleculeName << "not defined." << nl
+ << abort(FatalError);
+ }
- p.x() += pert*(rndGen_.scalar01() - 0.5);
- p.y() += pert*(rndGen_.scalar01() - 0.5);
- p.z() += pert*(rndGen_.scalar01() - 0.5);
+ const typename ParcelType::constantProperties& cP =
+ constProps(typeId);
- label cell = mesh_.findCell(p);
+ scalar numberDensity = numberDensities[i];
- vector U = equipartitionLinearVelocity
- (
- temperature,
- cP.mass()
- );
+ // Calculate the number of particles required
+ scalar particlesRequired = numberDensity*mesh_.cellVolumes()[cell];
- scalar Ei = equipartitionInternalEnergy
- (
- temperature,
- cP.internalDegreesOfFreedom()
- );
+ // Only integer numbers of particles can be inserted
+ label nParticlesToInsert = label(particlesRequired);
+
+ // Add another particle with a probability proportional to the
+ // remainder of taking the integer part of particlesRequired
+ if ((particlesRequired - nParticlesToInsert) > rndGen_.scalar01())
+ {
+ nParticlesToInsert++;
+ }
+
+ for (label pI = 0; pI < nParticlesToInsert; pI++)
+ {
+ // Choose a random point in a generic tetrahedron
+
+ scalar s = rndGen_.scalar01();
+ scalar t = rndGen_.scalar01();
+ scalar u = rndGen_.scalar01();
+
+ if (s + t > 1.0)
+ {
+ s = 1.0 - s;
+ t = 1.0 - t;
+ }
+
+ if (t + u > 1.0)
+ {
+ scalar tmp = u;
+ u = 1.0 - s - t;
+ t = 1.0 - tmp;
+ }
+ else if (s + t + u > 1.0)
+ {
+ scalar tmp = u;
+ u = s + t + u - 1.0;
+ s = 1.0 - t - tmp;
+ }
+
+ // Choose a tetrahedron to insert in, based on their relative
+ // volumes
+ scalar tetSelection = rndGen_.scalar01();
- U += velocity;
+ // Selected tetrahedron
+ label sTet = -1;
- if (cell >= 0)
+ forAll(cTetVFracs, tet)
+ {
+ sTet = tet;
+
+ if (cTetVFracs[tet] >= tetSelection)
{
- addNewParcel
- (
- p,
- U,
- Ei,
- cell,
- typeId
- );
+ break;
}
}
+
+ vector p =
+ (1 - s - t - u)*cC
+ + s*mesh_.points()[tetPtIs[sTet][0]]
+ + t*mesh_.points()[tetPtIs[sTet][1]]
+ + u*mesh_.points()[tetPtIs[sTet][2]];
+
+ vector U = equipartitionLinearVelocity
+ (
+ temperature,
+ cP.mass()
+ );
+
+ scalar Ei = equipartitionInternalEnergy
+ (
+ temperature,
+ cP.internalDegreesOfFreedom()
+ );
+
+ U += velocity;
+
+ addNewParcel
+ (
+ p,
+ U,
+ Ei,
+ cell,
+ typeId
+ );
}
}
}
@@ -276,7 +351,7 @@ void Foam::DsmcCloud<ParcelType>::collisions()
scalar selectedPairs =
collisionSelectionRemainder_[celli]
- + 0.5*nC*(nC-1)*nParticle_*sigmaTcRMax*deltaT
+ + 0.5*nC*(nC - 1)*nParticle_*sigmaTcRMax*deltaT
/mesh_.cellVolumes()[celli];
label nCandidates(selectedPairs);
@@ -551,6 +626,13 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
buildConstProps();
buildCellOccupancy();
+
+ // Initialise the collision selection remainder to a random value between 0
+ // and 1.
+ forAll(collisionSelectionRemainder_, i)
+ {
+ collisionSelectionRemainder_[i] = rndGen_.scalar01();
+ }
}
@@ -739,22 +821,27 @@ void Foam::DsmcCloud<ParcelType>::info() const
Info<< "Cloud name: " << this->name() << nl
<< " Number of dsmc particles = "
- << nDsmcParticles << nl
- << " Number of molecules = "
- << nMol << nl
- << " Mass in system = "
- << returnReduce(massInSystem(), sumOp<scalar>()) << 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
+ << nDsmcParticles
<< endl;
+
+ if (nDsmcParticles)
+ {
+ Info<< " Number of molecules = "
+ << nMol << nl
+ << " Mass in system = "
+ << returnReduce(massInSystem(), sumOp<scalar>()) << 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
+ << endl;
+ }
}
@@ -785,7 +872,11 @@ Foam::scalar Foam::DsmcCloud<ParcelType>::equipartitionInternalEnergy
{
scalar Ei = 0.0;
- if (iDof < 2.0 + SMALL && iDof > 2.0 - SMALL)
+ if (iDof < SMALL)
+ {
+ return Ei;
+ }
+ else if (iDof < 2.0 + SMALL && iDof > 2.0 - SMALL)
{
// Special case for iDof = 2, i.e. diatomics;
Ei = -log(rndGen_.scalar01())*kb*temperature;
@@ -796,7 +887,7 @@ Foam::scalar Foam::DsmcCloud<ParcelType>::equipartitionInternalEnergy
scalar energyRatio;
- scalar P;
+ scalar P = -1;
do
{
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
index b246b0499040fa15a68ab38ceb2543b8249e43f1..20bed1ee7b14429c66fcc90b2dcb52023ff73241 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
@@ -299,6 +299,12 @@ public:
scalar mass
) const;
+ inline scalarField maxwellianAverageSpeed
+ (
+ scalarField temperature,
+ scalar mass
+ ) const;
+
//- RMS particle speed
inline scalar maxwellianRMSSpeed
(
@@ -306,6 +312,12 @@ public:
scalar mass
) const;
+ inline scalarField maxwellianRMSSpeed
+ (
+ scalarField temperature,
+ scalar mass
+ ) const;
+
//- Most probable speed
inline scalar maxwellianMostProbableSpeed
(
@@ -313,6 +325,11 @@ public:
scalar mass
) const;
+ inline scalarField maxwellianMostProbableSpeed
+ (
+ scalarField temperature,
+ scalar mass
+ ) const;
// Sub-models
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
index 73899636e8f6cf2a7e47754a15c8cf8aeb711712..d39de4954faa1a062ab15da0d878e736d542b42e 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
@@ -302,6 +302,17 @@ inline Foam::scalar Foam::DsmcCloud<ParcelType>::maxwellianAverageSpeed
}
+template<class ParcelType>
+inline Foam::scalarField Foam::DsmcCloud<ParcelType>::maxwellianAverageSpeed
+(
+ scalarField 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
(
@@ -313,6 +324,17 @@ inline Foam::scalar Foam::DsmcCloud<ParcelType>::maxwellianRMSSpeed
}
+template<class ParcelType>
+inline Foam::scalarField Foam::DsmcCloud<ParcelType>::maxwellianRMSSpeed
+(
+ scalarField temperature,
+ scalar mass
+) const
+{
+ return sqrt(3.0*kb*temperature/mass);
+}
+
+
template<class ParcelType>
inline Foam::scalar
Foam::DsmcCloud<ParcelType>::maxwellianMostProbableSpeed
@@ -325,6 +347,18 @@ Foam::DsmcCloud<ParcelType>::maxwellianMostProbableSpeed
}
+template<class ParcelType>
+inline Foam::scalarField
+Foam::DsmcCloud<ParcelType>::maxwellianMostProbableSpeed
+(
+ scalarField 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
@@ -343,7 +377,7 @@ Foam::DsmcCloud<ParcelType>::rhoN() const
false
),
mesh_,
- dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0)
+ dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL)
)
);
@@ -380,7 +414,7 @@ Foam::DsmcCloud<ParcelType>::rhoM() const
false
),
mesh_,
- dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), 0.0)
+ dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), VSMALL)
)
);
@@ -568,7 +602,7 @@ Foam::DsmcCloud<ParcelType>::iDof() const
false
),
mesh_,
- dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0)
+ dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL)
)
);
diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
index f4203f4d66f7e85eaf335eee7c5c3b80521dd4e0..cd18d9064d3dfbe6d5126f6621c96d8af32a5df7 100644
--- a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
+++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
@@ -36,32 +36,27 @@ Foam::FreeStream<CloudType>::FreeStream
)
:
InflowBoundaryModel<CloudType>(dict, cloud, typeName),
- patchIndex_(),
- temperature_(readScalar(this->coeffDict().lookup("temperature"))),
- velocity_(this->coeffDict().lookup("velocity")),
+ patches_(),
moleculeTypeIds_(),
numberDensities_(),
particleFluxAccumulators_()
{
- word patchName = this->coeffDict().lookup("patch");
+ // Identify which patches to use
- patchIndex_ = cloud.mesh().boundaryMesh().findPatchID(patchName);
+ DynamicList<label> patches;
- const polyPatch& patch = cloud.mesh().boundaryMesh()[patchIndex_];
-
- if (patchIndex_ == -1)
+ forAll(cloud.mesh().boundaryMesh(), p)
{
- FatalErrorIn
- (
- "Foam::FreeStream<CloudType>::FreeStream"
- "("
- "const dictionary&, "
- "CloudType&"
- ")"
- ) << "patch " << patchName << " not found." << nl
- << abort(FatalError);
+ const polyPatch& patch = cloud.mesh().boundaryMesh()[p];
+
+ if (patch.type() == polyPatch::typeName)
+ {
+ patches.append(p);
+ }
}
+ patches_.transfer(patches);
+
const dictionary& numberDensitiesDict
(
this->coeffDict().subDict("numberDensities")
@@ -69,10 +64,24 @@ Foam::FreeStream<CloudType>::FreeStream
List<word> molecules(numberDensitiesDict.toc());
- numberDensities_.setSize(molecules.size());
+ // Initialise the particleFluxAccumulators_
+ particleFluxAccumulators_.setSize(patches_.size());
+
+ forAll(patches_, p)
+ {
+ const polyPatch& patch = cloud.mesh().boundaryMesh()[patches_[p]];
+
+ particleFluxAccumulators_[p] = List<Field<scalar> >
+ (
+ molecules.size(),
+ Field<scalar>(patch.size(), 0.0)
+ );
+ }
moleculeTypeIds_.setSize(molecules.size());
+ numberDensities_.setSize(molecules.size());
+
forAll(molecules, i)
{
numberDensities_[i] = readScalar
@@ -97,12 +106,6 @@ Foam::FreeStream<CloudType>::FreeStream
}
numberDensities_ /= cloud.nParticle();
-
- particleFluxAccumulators_.setSize
- (
- molecules.size(),
- Field<scalar>(patch.size(), 0)
- );
}
@@ -127,144 +130,262 @@ void Foam::FreeStream<CloudType>::inflow()
Random& rndGen(cloud.rndGen());
- const polyPatch& patch = mesh.boundaryMesh()[patchIndex_];
+ scalar sqrtPi = sqrt(mathematicalConstant::pi);
label particlesInserted = 0;
- // Add mass to the accumulators. negative face area dotted with the
- // velocity to point flux into the domain.
+ const volScalarField::GeometricBoundaryField& boundaryT
+ (
+ cloud.T().boundaryField()
+ );
- forAll(particleFluxAccumulators_, i)
- {
- particleFluxAccumulators_[i] +=
- -patch.faceAreas() & (velocity_*numberDensities_[i]*deltaT);
- }
+ const volVectorField::GeometricBoundaryField& boundaryU
+ (
+ cloud.U().boundaryField()
+ );
- forAll(patch, f)
+ forAll(patches_, p)
{
- // Loop over all faces as the outer loop to avoid calculating
- // geometrical properties multiple times.
+ label patchI = patches_[p];
- labelList faceVertices = patch[f];
+ const polyPatch& patch = mesh.boundaryMesh()[patchI];
- label nVertices = faceVertices.size();
+ // Add mass to the accumulators. negative face area dotted with the
+ // velocity to point flux into the domain.
- label globalFaceIndex = f + patch.start();
+ // Take a reference to the particleFluxAccumulator for this patch
+ List<Field<scalar> >& pFA = particleFluxAccumulators_[p];
- label cell = mesh.faceOwner()[globalFaceIndex];
-
- const vector& fC = patch.faceCentres()[f];
+ forAll(pFA, i)
+ {
+ label typeId = moleculeTypeIds_[i];
- scalar fA = mag(patch.faceAreas()[f]);
+ scalar mass = cloud.constProps(typeId).mass();
- // Cummulative triangle area fractions
- List<scalar> cTriAFracs(nVertices);
+ scalarField mostProbableSpeed
+ (
+ cloud.maxwellianMostProbableSpeed
+ (
+ boundaryT[patchI],
+ mass
+ )
+ );
+
+ // Dotting boundary velocity with the face unit normal (which points
+ // out of the domain, so it must be negated), dividing by the most
+ // probable speed to form molecularSpeedRatio * cosTheta
+
+ scalarField sCosTheta =
+ boundaryU[patchI]
+ & -patch.faceAreas()/mag(patch.faceAreas())
+ /mostProbableSpeed;
+
+ // From Bird eqn 4.22
+
+ pFA[i] +=
+ mag(patch.faceAreas())*numberDensities_[i]*deltaT
+ *mostProbableSpeed
+ *(
+ exp(-sqr(sCosTheta)) + sqrtPi*sCosTheta*(1 + erf(sCosTheta))
+ )
+ /(2.0*sqrtPi);
+ }
- for (label v = 0; v < nVertices - 1; v++)
+ forAll(patch, f)
{
- const point& vA = mesh.points()[faceVertices[v]];
+ // Loop over all faces as the outer loop to avoid calculating
+ // geometrical properties multiple times.
- const point& vB = mesh.points()[faceVertices[(v + 1)]];
+ labelList faceVertices = patch[f];
- cTriAFracs[v] =
- 0.5*mag((vA - fC)^(vB - fC))/fA
- + cTriAFracs[max((v - 1), 0)];
- }
+ label nVertices = faceVertices.size();
- // 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;
+ label globalFaceIndex = f + patch.start();
- // Normal unit vector *negative* so normal is pointing into the
- // domain
- vector nw = patch.faceAreas()[f];
- nw /= -mag(nw);
+ label cell = mesh.faceOwner()[globalFaceIndex];
- // 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 vector& fC = patch.faceCentres()[f];
- // 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 fA = mag(patch.faceAreas()[f]);
- forAll(particleFluxAccumulators_, i)
- {
- scalar& faceAccumulator = particleFluxAccumulators_[i][f];
+ // Cumulative triangle area fractions
+ List<scalar> cTriAFracs(nVertices);
- // Number of particles to insert
- label nI = max(label(faceAccumulator), 0);
+ for (label v = 0; v < nVertices - 1; v++)
+ {
+ const point& vA = mesh.points()[faceVertices[v]];
- faceAccumulator -= nI;
+ const point& vB = mesh.points()[faceVertices[(v + 1)]];
- label typeId = moleculeTypeIds_[i];
+ cTriAFracs[v] =
+ 0.5*mag((vA - fC)^(vB - fC))/fA
+ + cTriAFracs[max((v - 1), 0)];
+ }
- scalar mass = cloud.constProps(typeId).mass();
+ // 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 n = patch.faceAreas()[f];
+ n /= -mag(n);
+
+ // Wall tangential unit vector. Use the direction between the
+ // face centre and the first vertex in the list
+ vector t1 = fC - (mesh.points()[faceVertices[0]]);
+ t1 /= mag(t1);
- for (label n = 0; n < nI; n++)
+ // Other tangential unit vector. Rescaling in case face is not
+ // flat and n and t1 aren't perfectly orthogonal
+ vector t2 = n^t1;
+ t2 /= mag(t2);
+
+ scalar faceTemperature = boundaryT[patchI][f];
+
+ const vector& faceVelocity = boundaryU[patchI][f];
+
+ forAll(pFA, i)
{
- // Choose a triangle to insert on, based on their relative area
+ scalar& faceAccumulator = pFA[i][f];
- scalar triSelection = rndGen.scalar01();
+ // Number of whole particles to insert
+ label nI = max(label(faceAccumulator), 0);
- // Selected triangle
- label sTri = -1;
+ // Add another particle with a probability proportional to the
+ // remainder of taking the integer part of faceAccumulator
+ if ((faceAccumulator - nI) > rndGen.scalar01())
+ {
+ nI++;
+ }
+
+ faceAccumulator -= nI;
+
+ label typeId = moleculeTypeIds_[i];
+
+ scalar mass = cloud.constProps(typeId).mass();
- forAll(cTriAFracs, tri)
+ for (label i = 0; i < nI; i++)
{
- sTri = tri;
+ // Choose a triangle to insert on, based on their relative
+ // area
- if (cTriAFracs[tri] >= triSelection)
+ scalar triSelection = rndGen.scalar01();
+
+ // Selected triangle
+ label sTri = -1;
+
+ forAll(cTriAFracs, tri)
{
- break;
+ sTri = tri;
+
+ if (cTriAFracs[tri] >= triSelection)
+ {
+ break;
+ }
}
- }
- // 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
+ // 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
- const point& A = fC;
- const point& B = mesh.points()[faceVertices[sTri]];
- const point& C =
+ const point& A = fC;
+ const point& B = mesh.points()[faceVertices[sTri]];
+ const point& C =
mesh.points()[faceVertices[(sTri + 1) % nVertices]];
- scalar s = rndGen.scalar01();
- scalar t = sqrt(rndGen.scalar01());
+ scalar s = rndGen.scalar01();
+ scalar t = sqrt(rndGen.scalar01());
- point p = (1 - t)*A + (1 - s)*t*B + s*t*C;
+ point p = (1 - t)*A + (1 - s)*t*B + s*t*C;
- vector U =
- sqrt(CloudType::kb*temperature_/mass)
- *(
- rndGen.GaussNormal()*tw1
- + rndGen.GaussNormal()*tw2
- - sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw
+ // Velocity generation
+
+ scalar mostProbableSpeed
+ (
+ cloud.maxwellianMostProbableSpeed
+ (
+ faceTemperature,
+ mass
+ )
);
- U += velocity_;
+ scalar sCosTheta = (faceVelocity & n)/mostProbableSpeed;
- scalar Ei = cloud.equipartitionInternalEnergy
- (
- temperature_,
- cloud.constProps(typeId).internalDegreesOfFreedom()
- );
+ // Coefficients required for Bird eqn 12.5
+ scalar uNormProbCoeffA =
+ sCosTheta + sqrt(sqr(sCosTheta) + 2.0);
- cloud.addNewParcel
- (
- p,
- U,
- Ei,
- cell,
- typeId
- );
-
- particlesInserted++;
+ scalar uNormProbCoeffB =
+ 0.5*
+ (
+ 1.0
+ + sCosTheta*(sCosTheta - sqrt(sqr(sCosTheta) + 2.0))
+ );
+
+ // Equivalent to the QA value in Bird's DSMC3.FOR
+ scalar randomScaling = 3.0;
+
+ if (sCosTheta < -3)
+ {
+ randomScaling = mag(sCosTheta) + 1;
+ }
+
+ scalar P = -1;
+
+ // Normalised candidates for the normal direction velocity
+ // component
+ scalar uNormal;
+ scalar uNormalThermal;
+
+ // Select a velocity using Bird eqn 12.5
+ do
+ {
+ uNormalThermal =
+ randomScaling*(2.0*rndGen.scalar01() - 1);
+
+ uNormal = uNormalThermal + sCosTheta;
+
+ if (uNormal < 0.0)
+ {
+ P = -1;
+ }
+ else
+ {
+ P = 2.0*uNormal/uNormProbCoeffA
+ *exp(uNormProbCoeffB - sqr(uNormalThermal));
+ }
+
+ } while (P < rndGen.scalar01());
+
+ vector U =
+ sqrt(CloudType::kb*faceTemperature/mass)
+ *(
+ rndGen.GaussNormal()*t1
+ + rndGen.GaussNormal()*t2
+ )
+ + mostProbableSpeed*uNormal*n;
+
+ scalar Ei = cloud.equipartitionInternalEnergy
+ (
+ faceTemperature,
+ cloud.constProps(typeId).internalDegreesOfFreedom()
+ );
+
+ cloud.addNewParcel
+ (
+ p,
+ U,
+ Ei,
+ cell,
+ typeId
+ );
+
+ particlesInserted++;
+ }
}
}
}
diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.H b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.H
index f1bd4cc32133490ea027136a6de2f36798898574..03c2fe2abb72fe657ef9905162585099c445693d 100644
--- a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.H
+++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.H
@@ -26,8 +26,10 @@ Class
Foam::FreeStream
Description
- Inserting new particles across the faces of a specified patch for a free
- stream. Uniform values of temperature, velocity and number densities
+ Inserting new particles across the faces of a all patched of type
+ "patch" for a free stream. Uniform values number density, temperature
+ and velocity sourced face-by-face from the boundaryT and boundaryU fields
+ of the cloud.
\*---------------------------------------------------------------------------*/
@@ -52,14 +54,8 @@ class FreeStream
{
// Private data
- //- Index of patch to introduce particles across
- label patchIndex_;
-
- //- Temperature of the free stream
- scalar temperature_;
-
- //- Velocity of the free stream
- vector velocity_;
+ //- The indices of patches to introduce molecules across
+ labelList patches_;
//- The molecule types to be introduced
List<label> moleculeTypeIds_;
@@ -67,10 +63,13 @@ class FreeStream
//- The number density of the species in the inflow
Field<scalar> numberDensities_;
- //- A List of Fields, one Field for every species to be introduced, each
- // field entry corresponding to a face on the patch to be injected
- // across.
- List<Field<scalar> > particleFluxAccumulators_;
+ //- A List of Lists of Fields specifying carry-over of mass flux from
+ // one timestep to the next
+ // + Outer List - one inner List for each patch
+ // + Inner List - one Field for every species to be introduced
+ // + Each field entry corresponding to a face to be injected across
+ // with a particular species
+ List<List<Field<scalar> > > particleFluxAccumulators_;
public:
diff --git a/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C
index 3aaa9c514a462f7d2850d0035ab34478d70e7817..273c88f3699822572bce150ad303b3224f5574ce 100644
--- a/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C
+++ b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C
@@ -137,7 +137,7 @@ void Foam::dsmcFields::write()
iDofMeanName
);
- if (min(rhoNMean).value() > VSMALL)
+ if (min(mag(rhoNMean)).value() > VSMALL)
{
Info<< "Calculating dsmcFields." << endl;
@@ -223,10 +223,9 @@ void Foam::dsmcFields::write()
}
else
{
- Info<< "Small or negative value (" << min(rhoNMean)
+ Info<< "Small value (" << min(mag(rhoNMean))
<< ") found in rhoNMean field. "
- << "Not calculating dsmcFields to avoid division by zero "
- << "or invalid results."
+ << "Not calculating dsmcFields to avoid division by zero."
<< endl;
}
}
diff --git a/src/sampling/sampledSet/sampledSets/sampledSets.C b/src/sampling/sampledSet/sampledSets/sampledSets.C
index 3f7b7a3323ede9cbc888a87d1929da1ba887e405..4d71403277a790c84179c0f5612cbff226f6754e 100644
--- a/src/sampling/sampledSet/sampledSets/sampledSets.C
+++ b/src/sampling/sampledSet/sampledSets/sampledSets.C
@@ -236,8 +236,6 @@ Foam::sampledSets::sampledSets
loadFromFiles_(loadFromFiles),
outputPath_(fileName::null),
searchEngine_(mesh_, true),
-// pMeshPtr_(NULL),
-// pInterpPtr_(NULL),
fieldNames_(),
interpolationScheme_(word::null),
writeFormat_(word::null)
@@ -250,6 +248,10 @@ Foam::sampledSets::sampledSets
{
outputPath_ = mesh_.time().path()/name_;
}
+ if (mesh_.name() != fvMesh::defaultRegion)
+ {
+ outputPath_ = outputPath_/mesh_.name();
+ }
read(dict);
}