Changing FreeStream to create inflow on all patches of type patch. Implemented Bird eqn 4.22 for the number flux and eqn 12.5 for the velocity distibution. Drawing T and U for the FreeStream from the boundaryT and boundaryU fields.

This commit is contained in:
graham
2009-04-16 19:18:15 +01:00
parent a8b37fc4a7
commit 239b954ad3
6 changed files with 335 additions and 160 deletions

View File

@ -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;
}

View File

@ -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>
@ -625,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();
}
}
@ -813,8 +821,12 @@ void Foam::DsmcCloud<ParcelType>::info() const
Info<< "Cloud name: " << this->name() << nl
<< " Number of dsmc particles = "
<< nDsmcParticles << nl
<< " Number of molecules = "
<< nDsmcParticles
<< endl;
if (nDsmcParticles)
{
Info<< " Number of molecules = "
<< nMol << nl
<< " Mass in system = "
<< returnReduce(massInSystem(), sumOp<scalar>()) << nl
@ -827,8 +839,9 @@ void Foam::DsmcCloud<ParcelType>::info() const
<< " Average internal energy = "
<< internalEnergy/nMol << nl
<< " Average total energy = "
<< (internalEnergy + linearKineticEnergy)/nMol << nl
<< (internalEnergy + linearKineticEnergy)/nMol
<< endl;
}
}
@ -874,7 +887,7 @@ Foam::scalar Foam::DsmcCloud<ParcelType>::equipartitionInternalEnergy
scalar energyRatio;
scalar P;
scalar P = -1;
do
{

View File

@ -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

View File

@ -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

View File

@ -36,31 +36,26 @@ 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
(
@ -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,17 +130,65 @@ void Foam::FreeStream<CloudType>::inflow()
Random& rndGen(cloud.rndGen());
const polyPatch& patch = mesh.boundaryMesh()[patchIndex_];
scalar sqrtPi = sqrt(mathematicalConstant::pi);
label particlesInserted = 0;
const volScalarField::GeometricBoundaryField& boundaryT
(
cloud.T().boundaryField()
);
const volVectorField::GeometricBoundaryField& boundaryU
(
cloud.U().boundaryField()
);
forAll(patches_, p)
{
label patchI = patches_[p];
const polyPatch& patch = mesh.boundaryMesh()[patchI];
// Add mass to the accumulators. negative face area dotted with the
// velocity to point flux into the domain.
forAll(particleFluxAccumulators_, i)
// Take a reference to the particleFluxAccumulator for this patch
List<Field<scalar> >& pFA = particleFluxAccumulators_[p];
forAll(pFA, i)
{
particleFluxAccumulators_[i] +=
-patch.faceAreas() & (velocity_*numberDensities_[i]*deltaT);
label typeId = moleculeTypeIds_[i];
scalar mass = cloud.constProps(typeId).mass();
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)
@ -177,22 +228,26 @@ void Foam::FreeStream<CloudType>::inflow()
// Normal unit vector *negative* so normal is pointing into the
// domain
vector nw = patch.faceAreas()[f];
nw /= -mag(nw);
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 tw1 = fC - (mesh.points()[faceVertices[0]]);
tw1 /= mag(tw1);
vector t1 = fC - (mesh.points()[faceVertices[0]]);
t1 /= mag(t1);
// 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);
// flat and n and t1 aren't perfectly orthogonal
vector t2 = n^t1;
t2 /= mag(t2);
forAll(particleFluxAccumulators_, i)
scalar faceTemperature = boundaryT[patchI][f];
const vector& faceVelocity = boundaryU[patchI][f];
forAll(pFA, i)
{
scalar& faceAccumulator = particleFluxAccumulators_[i][f];
scalar& faceAccumulator = pFA[i][f];
// Number of particles to insert
label nI = max(label(faceAccumulator), 0);
@ -203,9 +258,10 @@ void Foam::FreeStream<CloudType>::inflow()
scalar mass = cloud.constProps(typeId).mass();
for (label n = 0; n < nI; n++)
for (label i = 0; i < nI; i++)
{
// Choose a triangle to insert on, based on their relative area
// Choose a triangle to insert on, based on their relative
// area
scalar triSelection = rndGen.scalar01();
@ -239,19 +295,74 @@ void Foam::FreeStream<CloudType>::inflow()
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;
// Coefficients required for Bird eqn 12.5
scalar uNormProbCoeffA =
sCosTheta + sqrt(sqr(sCosTheta) + 2.0);
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)
{
continue;
}
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
(
temperature_,
faceTemperature,
cloud.constProps(typeId).internalDegreesOfFreedom()
);
@ -268,6 +379,7 @@ void Foam::FreeStream<CloudType>::inflow()
}
}
}
}
reduce(particlesInserted, sumOp<label>());

View File

@ -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: