mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev
This commit is contained in:
mattijs
@ -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;
|
||||
}
|
||||
|
||||
@ -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;
|
||||
|
||||
@ -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;
|
||||
|
||||
scalar y0 = mesh_.bounds().min().y();
|
||||
scalar yR = mesh_.bounds().max().y() - y0;
|
||||
|
||||
scalar z0 = mesh_.bounds().min().z();
|
||||
scalar zR = mesh_.bounds().max().z() - z0;
|
||||
|
||||
forAll(molecules, i)
|
||||
forAll(mesh_.cells(), cell)
|
||||
{
|
||||
const word& moleculeName(molecules[i]);
|
||||
const vector& cC = mesh_.cellCentres()[cell];
|
||||
const labelList& cellFaces = mesh_.cells()[cell];
|
||||
const scalar cV = mesh_.cellVolumes()[cell];
|
||||
|
||||
label typeId(findIndex(typeIdList_, moleculeName));
|
||||
label nTets = 0;
|
||||
|
||||
if (typeId == -1)
|
||||
// Each face is split into nEdges (or nVertices) - 2 tets.
|
||||
forAll(cellFaces, face)
|
||||
{
|
||||
FatalErrorIn("Foam::DsmcCloud<ParcelType>::initialise")
|
||||
<< "typeId " << moleculeName << "not defined." << nl
|
||||
<< abort(FatalError);
|
||||
nTets += mesh_.faces()[cellFaces[face]].size() - 2;
|
||||
}
|
||||
|
||||
const typename ParcelType::constantProperties& cP = constProps(typeId);
|
||||
// Calculate the cumulative tet volumes circulating around the cell and
|
||||
// record the vertex labels of each.
|
||||
scalarList cTetVFracs(nTets, 0.0);
|
||||
|
||||
scalar numberDensity = numberDensities[i];
|
||||
List<labelList> tetPtIs(nTets, labelList(3,-1));
|
||||
|
||||
scalar spacing = pow(numberDensity,-(1.0/3.0));
|
||||
// Keep track of which tet this is.
|
||||
label tet = 0;
|
||||
|
||||
int ni = label(xR/spacing) + 1;
|
||||
int nj = label(yR/spacing) + 1;
|
||||
int nk = label(zR/spacing) + 1;
|
||||
|
||||
vector delta(xR/ni, yR/nj, zR/nk);
|
||||
|
||||
scalar pert = spacing;
|
||||
|
||||
for (int i = 0; i < ni; i++)
|
||||
forAll(cellFaces, face)
|
||||
{
|
||||
for (int j = 0; j < nj; j++)
|
||||
const labelList& facePoints = mesh_.faces()[cellFaces[face]];
|
||||
|
||||
label pointI = 1;
|
||||
while ((pointI + 1) < facePoints.size())
|
||||
{
|
||||
for (int k = 0; k < nk; k++)
|
||||
|
||||
const vector& pA = mesh_.points()[facePoints[0]];
|
||||
const vector& pB = mesh_.points()[facePoints[pointI]];
|
||||
const vector& pC = mesh_.points()[facePoints[pointI + 1]];
|
||||
|
||||
cTetVFracs[tet] =
|
||||
mag(((pA - cC) ^ (pB - cC)) & (pC - cC))/(cV*6.0)
|
||||
+ cTetVFracs[max((tet - 1),0)];
|
||||
|
||||
tetPtIs[tet][0] = facePoints[0];
|
||||
tetPtIs[tet][1] = facePoints[pointI];
|
||||
tetPtIs[tet][2] = facePoints[pointI + 1];
|
||||
|
||||
pointI++;
|
||||
tet++;
|
||||
}
|
||||
}
|
||||
|
||||
// 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;
|
||||
|
||||
forAll(molecules, i)
|
||||
{
|
||||
const word& moleculeName(molecules[i]);
|
||||
|
||||
label typeId(findIndex(typeIdList_, moleculeName));
|
||||
|
||||
if (typeId == -1)
|
||||
{
|
||||
FatalErrorIn("Foam::DsmcCloud<ParcelType>::initialise")
|
||||
<< "typeId " << moleculeName << "not defined." << nl
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
const typename ParcelType::constantProperties& cP =
|
||||
constProps(typeId);
|
||||
|
||||
scalar numberDensity = numberDensities[i];
|
||||
|
||||
// Calculate the number of particles required
|
||||
scalar particlesRequired = numberDensity*mesh_.cellVolumes()[cell];
|
||||
|
||||
// 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)
|
||||
{
|
||||
point p
|
||||
(
|
||||
x0 + (i + 0.5)*delta.x(),
|
||||
y0 + (j + 0.5)*delta.y(),
|
||||
z0 + (k + 0.5)*delta.z()
|
||||
);
|
||||
s = 1.0 - s;
|
||||
t = 1.0 - t;
|
||||
}
|
||||
|
||||
p.x() += pert*(rndGen_.scalar01() - 0.5);
|
||||
p.y() += pert*(rndGen_.scalar01() - 0.5);
|
||||
p.z() += pert*(rndGen_.scalar01() - 0.5);
|
||||
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;
|
||||
}
|
||||
|
||||
label cell = mesh_.findCell(p);
|
||||
// Choose a tetrahedron to insert in, based on their relative
|
||||
// volumes
|
||||
scalar tetSelection = rndGen_.scalar01();
|
||||
|
||||
vector U = equipartitionLinearVelocity
|
||||
(
|
||||
temperature,
|
||||
cP.mass()
|
||||
);
|
||||
// Selected tetrahedron
|
||||
label sTet = -1;
|
||||
|
||||
scalar Ei = equipartitionInternalEnergy
|
||||
(
|
||||
temperature,
|
||||
cP.internalDegreesOfFreedom()
|
||||
);
|
||||
forAll(cTetVFracs, tet)
|
||||
{
|
||||
sTet = tet;
|
||||
|
||||
U += velocity;
|
||||
|
||||
if (cell >= 0)
|
||||
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
|
||||
{
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
@ -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)
|
||||
)
|
||||
);
|
||||
|
||||
|
||||
@ -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)
|
||||
const volVectorField::GeometricBoundaryField& boundaryU
|
||||
(
|
||||
cloud.U().boundaryField()
|
||||
);
|
||||
|
||||
forAll(patches_, p)
|
||||
{
|
||||
particleFluxAccumulators_[i] +=
|
||||
-patch.faceAreas() & (velocity_*numberDensities_[i]*deltaT);
|
||||
}
|
||||
label patchI = patches_[p];
|
||||
|
||||
forAll(patch, f)
|
||||
{
|
||||
// Loop over all faces as the outer loop to avoid calculating
|
||||
// geometrical properties multiple times.
|
||||
const polyPatch& patch = mesh.boundaryMesh()[patchI];
|
||||
|
||||
labelList faceVertices = patch[f];
|
||||
// Add mass to the accumulators. negative face area dotted with the
|
||||
// velocity to point flux into the domain.
|
||||
|
||||
label nVertices = faceVertices.size();
|
||||
// Take a reference to the particleFluxAccumulator for this patch
|
||||
List<Field<scalar> >& pFA = particleFluxAccumulators_[p];
|
||||
|
||||
label globalFaceIndex = f + patch.start();
|
||||
|
||||
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++)
|
||||
forAll(pFA, i)
|
||||
{
|
||||
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];
|
||||
|
||||
// Number of particles to insert
|
||||
label nI = max(label(faceAccumulator), 0);
|
||||
|
||||
faceAccumulator -= nI;
|
||||
|
||||
label typeId = moleculeTypeIds_[i];
|
||||
|
||||
scalar mass = cloud.constProps(typeId).mass();
|
||||
|
||||
for (label n = 0; n < nI; n++)
|
||||
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);
|
||||
}
|
||||
|
||||
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();
|
||||
|
||||
label cell = mesh.faceOwner()[globalFaceIndex];
|
||||
|
||||
const vector& fC = patch.faceCentres()[f];
|
||||
|
||||
scalar fA = mag(patch.faceAreas()[f]);
|
||||
|
||||
// Cumulative triangle area fractions
|
||||
List<scalar> cTriAFracs(nVertices);
|
||||
|
||||
for (label v = 0; v < nVertices - 1; v++)
|
||||
{
|
||||
// Choose a triangle to insert on, based on their relative area
|
||||
const point& vA = mesh.points()[faceVertices[v]];
|
||||
|
||||
scalar triSelection = rndGen.scalar01();
|
||||
const point& vB = mesh.points()[faceVertices[(v + 1)]];
|
||||
|
||||
// Selected triangle
|
||||
label sTri = -1;
|
||||
cTriAFracs[v] =
|
||||
0.5*mag((vA - fC)^(vB - fC))/fA
|
||||
+ cTriAFracs[max((v - 1), 0)];
|
||||
}
|
||||
|
||||
forAll(cTriAFracs, tri)
|
||||
// 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);
|
||||
|
||||
// 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)
|
||||
{
|
||||
scalar& faceAccumulator = pFA[i][f];
|
||||
|
||||
// Number of whole particles to insert
|
||||
label nI = max(label(faceAccumulator), 0);
|
||||
|
||||
// Add another particle with a probability proportional to the
|
||||
// remainder of taking the integer part of faceAccumulator
|
||||
if ((faceAccumulator - nI) > rndGen.scalar01())
|
||||
{
|
||||
sTri = tri;
|
||||
|
||||
if (cTriAFracs[tri] >= triSelection)
|
||||
{
|
||||
break;
|
||||
}
|
||||
nI++;
|
||||
}
|
||||
|
||||
// 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
|
||||
faceAccumulator -= nI;
|
||||
|
||||
const point& A = fC;
|
||||
const point& B = mesh.points()[faceVertices[sTri]];
|
||||
const point& C =
|
||||
label typeId = moleculeTypeIds_[i];
|
||||
|
||||
scalar mass = cloud.constProps(typeId).mass();
|
||||
|
||||
for (label i = 0; i < nI; i++)
|
||||
{
|
||||
// Choose a triangle to insert on, based on their relative
|
||||
// area
|
||||
|
||||
scalar triSelection = rndGen.scalar01();
|
||||
|
||||
// Selected triangle
|
||||
label sTri = -1;
|
||||
|
||||
forAll(cTriAFracs, tri)
|
||||
{
|
||||
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
|
||||
|
||||
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
|
||||
);
|
||||
scalar uNormProbCoeffB =
|
||||
0.5*
|
||||
(
|
||||
1.0
|
||||
+ sCosTheta*(sCosTheta - sqrt(sqr(sCosTheta) + 2.0))
|
||||
);
|
||||
|
||||
particlesInserted++;
|
||||
// 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++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -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:
|
||||
|
||||
@ -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;
|
||||
}
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user