ENH: PatchflowRateInjection model - refactored and improved

This commit is contained in:
andy
2013-07-12 13:14:55 +01:00
parent f67b047dbb
commit 611959fae1
2 changed files with 81 additions and 160 deletions

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -27,6 +27,7 @@ License
#include "TimeDataEntry.H"
#include "distributionModel.H"
#include "mathematicalConstants.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -39,19 +40,23 @@ Foam::PatchFlowRateInjection<CloudType>::PatchFlowRateInjection
)
:
InjectionModel<CloudType>(dict, owner, modelName,typeName),
patchName_(this->coeffDict().lookup("patchName")),
patchId_(owner.mesh().boundaryMesh().findPatchID(patchName_)),
patchArea_(0.0),
patchNormal_(vector::zero),
patchInjectionBase(owner.mesh(), this->coeffDict().lookup("patchName")),
phiName_(this->coeffDict().template lookupOrDefault<word>("phi", "phi")),
rhoName_(this->coeffDict().template lookupOrDefault<word>("rho", "rho")),
duration_(readScalar(this->coeffDict().lookup("duration"))),
concentration_(readScalar(this->coeffDict().lookup("concentration"))),
parcelsPerSecond_
concentration_
(
readScalar(this->coeffDict().lookup("parcelsPerSecond"))
TimeDataEntry<scalar>
(
owner.db().time(),
"concentration",
this->coeffDict()
)
),
parcelConcentration_
(
readScalar(this->coeffDict().lookup("parcelConcentration"))
),
U0_(vector::zero),
sizeDistribution_
(
distributionModels::distributionModel::New
@ -59,45 +64,11 @@ Foam::PatchFlowRateInjection<CloudType>::PatchFlowRateInjection
this->coeffDict().subDict("sizeDistribution"),
owner.rndGen()
)
),
cellOwners_(),
fraction_(1.0),
pMeanVolume_(0.0)
)
{
if (patchId_ < 0)
{
FatalErrorIn
(
"PatchFlowRateInjection<CloudType>::PatchFlowRateInjection"
"("
"const dictionary&, "
"CloudType&"
")"
) << "Requested patch " << patchName_ << " not found" << nl
<< "Available patches are: " << owner.mesh().boundaryMesh().names()
<< nl << exit(FatalError);
}
const polyPatch& patch = owner.mesh().boundaryMesh()[patchId_];
duration_ = owner.db().time().userTimeToTime(duration_);
updateMesh();
// TODO: retrieve mean diameter from distrution model
scalar pMeanDiameter =
readScalar(this->coeffDict().lookup("meanParticleDiameter"));
pMeanVolume_ = constant::mathematical::pi*pow3(pMeanDiameter)/6.0;
// patch geometry
label patchSize = cellOwners_.size();
label totalPatchSize = patchSize;
reduce(totalPatchSize, sumOp<label>());
fraction_ = scalar(patchSize)/totalPatchSize;
patchArea_ = gSum(mag(patch.faceAreas()));
patchNormal_ = gSum(patch.faceNormals())/totalPatchSize;
patchNormal_ /= mag(patchNormal_);
patchInjectionBase::updateMesh(owner.mesh());
// Re-initialise total mass/volume to inject to zero
// - will be reset during each injection
@ -113,20 +84,13 @@ Foam::PatchFlowRateInjection<CloudType>::PatchFlowRateInjection
)
:
InjectionModel<CloudType>(im),
patchName_(im.patchName_),
patchId_(im.patchId_),
patchArea_(im.patchArea_),
patchNormal_(im.patchNormal_),
patchInjectionBase(im),
phiName_(im.phiName_),
rhoName_(im.rhoName_),
duration_(im.duration_),
concentration_(im.concentration_),
parcelsPerSecond_(im.parcelsPerSecond_),
U0_(im.U0_),
sizeDistribution_(im.sizeDistribution_().clone().ptr()),
cellOwners_(im.cellOwners_),
fraction_(im.fraction_),
pMeanVolume_(im.pMeanVolume_)
parcelConcentration_(im.parcelConcentration_),
sizeDistribution_(im.sizeDistribution_().clone().ptr())
{}
@ -142,9 +106,7 @@ Foam::PatchFlowRateInjection<CloudType>::~PatchFlowRateInjection()
template<class CloudType>
void Foam::PatchFlowRateInjection<CloudType>::updateMesh()
{
// Set/cache the injector cells
const polyPatch& patch = this->owner().mesh().boundaryMesh()[patchId_];
cellOwners_ = patch.faceCells();
patchInjectionBase::updateMesh(this->owner().mesh());
}
@ -155,6 +117,36 @@ Foam::scalar Foam::PatchFlowRateInjection<CloudType>::timeEnd() const
}
template<class CloudType>
Foam::scalar Foam::PatchFlowRateInjection<CloudType>::flowRate() const
{
const polyMesh& mesh = this->owner().mesh();
const surfaceScalarField& phi =
mesh.lookupObject<surfaceScalarField>(phiName_);
const scalarField& phip = phi.boundaryField()[patchId_];
scalar flowRateIn = 0.0;
if (phi.dimensions() == dimVelocity*dimArea)
{
flowRateIn = max(0.0, -sum(phip));
}
else
{
const volScalarField& rho =
mesh.lookupObject<volScalarField>(rhoName_);
const scalarField& rhop = rho.boundaryField()[patchId_];
flowRateIn = max(0.0, -sum(phip/rhop));
}
reduce(flowRateIn, sumOp<scalar>());
return flowRateIn;
}
template<class CloudType>
Foam::label Foam::PatchFlowRateInjection<CloudType>::parcelsToInject
(
@ -164,10 +156,11 @@ Foam::label Foam::PatchFlowRateInjection<CloudType>::parcelsToInject
{
if ((time0 >= 0.0) && (time0 < duration_))
{
scalar nParcels = fraction_*(time1 - time0)*parcelsPerSecond_;
scalar dt = time1 - time0;
cachedRandom& rnd = this->owner().rndGen();
scalar c = concentration_.value(0.5*(time0 + time1));
scalar nParcels = fraction_*parcelConcentration_*c*flowRate()*dt;
label nParcelsToInject = floor(nParcels);
// Inject an additional parcel with a probability based on the
@ -177,7 +170,7 @@ Foam::label Foam::PatchFlowRateInjection<CloudType>::parcelsToInject
nParcelsToInject > 0
&& (
nParcels - scalar(nParcelsToInject)
> rnd.position(scalar(0), scalar(1))
> this->owner().rndGen().position(scalar(0), scalar(1))
)
)
{
@ -204,38 +197,11 @@ Foam::scalar Foam::PatchFlowRateInjection<CloudType>::volumeToInject
if ((time0 >= 0.0) && (time0 < duration_))
{
const polyMesh& mesh = this->owner().mesh();
scalar c = concentration_.value(0.5*(time0 + time1));
const surfaceScalarField& phi =
mesh.lookupObject<surfaceScalarField>(phiName_);
const scalarField& phip = phi.boundaryField()[patchId_];
scalar carrierVolume = 0.0;
if (phi.dimensions() == dimVelocity*dimArea)
{
const scalar flowRateIn = max(0.0, -sum(phip));
U0_ = -patchNormal_*flowRateIn/patchArea_;
carrierVolume = (time1 - time0)*flowRateIn;
}
else
{
const volScalarField& rho =
mesh.lookupObject<volScalarField>(rhoName_);
const scalarField& rhop = rho.boundaryField()[patchId_];
const scalar flowRateIn = max(0.0, -sum(phip/rhop));
U0_ = -patchNormal_*flowRateIn/patchArea_;
carrierVolume = (time1 - time0)*flowRateIn;
}
const scalar newParticles = concentration_*carrierVolume;
volume = pMeanVolume_*newParticles;
volume = c*(time1 - time0)*flowRate();
}
reduce(volume, sumOp<scalar>());
this->volumeTotal_ = volume;
this->massTotal_ = volume*this->owner().constProps().rho0();
@ -255,39 +221,15 @@ void Foam::PatchFlowRateInjection<CloudType>::setPositionAndCell
label& tetPtI
)
{
if (cellOwners_.size() > 0)
{
cachedRandom& rnd = this->owner().rndGen();
label cellI = rnd.position<label>(0, cellOwners_.size() - 1);
cellOwner = cellOwners_[cellI];
// The position is between the face and cell centre, which could be
// in any tet of the decomposed cell, so arbitrarily choose the first
// face of the cell as the tetFace and the first point after the base
// point on the face as the tetPt. The tracking will pick the cell
// consistent with the motion in the firsttracking step.
tetFaceI = this->owner().mesh().cells()[cellOwner][0];
tetPtI = 1;
// position perturbed between cell and patch face centres
const vector& pc = this->owner().mesh().C()[cellOwner];
const vector& pf =
this->owner().mesh().Cf().boundaryField()[patchId_][cellI];
const scalar a = rnd.sample01<scalar>();
const vector d = pf - pc;
position = pc + 0.5*a*d;
}
else
{
cellOwner = -1;
tetFaceI = -1;
tetPtI = -1;
// dummy position
position = pTraits<vector>::max;
}
patchInjectionBase::setPositionAndCell
(
this->owner().mesh(),
this->owner().rndGen(),
position,
cellOwner,
tetFaceI,
tetPtI
);
}
@ -300,8 +242,8 @@ void Foam::PatchFlowRateInjection<CloudType>::setProperties
typename CloudType::parcelType& parcel
)
{
// set particle velocity
parcel.U() = U0_;
// set particle velocity to carrier velocity
parcel.U() = this->owner().U()[parcel.cell()];
// set particle diameter
parcel.d() = sizeDistribution_->sample();

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -32,10 +32,10 @@ Description
- Total mass to inject
- Name of patch
- Injection duration
- Initial parcel velocity
- Injection target concentration/carrier volume flow rate
- Initial parcel velocity given by local flow velocity
- Parcel diameters obtained by distribution model
- Parcels injected at cell centres adjacent to patch
- Parcels injected randomly across the patch
SourceFiles
PatchFlowRateInjection.C
@ -46,15 +46,14 @@ SourceFiles
#define PatchFlowRateInjection_H
#include "InjectionModel.H"
#include "patchInjectionBase.H"
#include "TimeDataEntry.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
template<class Type>
class TimeDataEntry;
class distributionModel;
/*---------------------------------------------------------------------------*\
@ -64,22 +63,11 @@ class distributionModel;
template<class CloudType>
class PatchFlowRateInjection
:
public InjectionModel<CloudType>
public InjectionModel<CloudType>,
public patchInjectionBase
{
// Private data
//- Name of patch
const word patchName_;
//- Id of patch
const label patchId_;
//- Patch area
scalar patchArea_;
//- Patch normal direction
vector patchNormal_;
//- Name of carrier (mass or volume) flux field
const word phiName_;
@ -89,27 +77,15 @@ class PatchFlowRateInjection
//- Injection duration [s]
scalar duration_;
//- Concentration of particles to carrier [] (particles/m3)
const scalar concentration_;
//- Concentration profile of particle volume to carrier volume [-]
const TimeDataEntry<scalar> concentration_;
//- Number of parcels to introduce per second []
const label parcelsPerSecond_;
//- Initial parcel velocity [m/s]
vector U0_;
//- Parcels to introduce per unit volume flow rate m3 [n/m3]
const scalar parcelConcentration_;
//- Parcel size distribution model
const autoPtr<distributionModels::distributionModel> sizeDistribution_;
//- List of cell labels corresponding to injector positions
labelList cellOwners_;
//- Fraction of injection controlled by this processor
scalar fraction_;
//- Mean particle volume TODO: temporary measure - return from PDF
scalar pMeanVolume_;
public:
@ -152,6 +128,9 @@ public:
//- Return the end-of-injection time
scalar timeEnd() const;
//- Return the total volumetric flow rate across the patch [m3/s]
virtual scalar flowRate() const;
//- Number of parcels to introduce relative to SOI
virtual label parcelsToInject(const scalar time0, const scalar time1);