我在几年前做了个一拉格朗日粒子直径变化的代码。我现在已经记不起来太多了。能想起来的是：
粒子直径跟流场的湍流动能耗散率有关系
植入的是一个跟PBM有关的粒径变化方程，考虑了coalescense和破碎
代码可以编译并且测试过可以体现粒径变化
可以把下面的MomentumCloud.C：替换到OpenFOAM原来的代码然后编译。搞明白之后看看思路然后适配到你们的算法上。
MomentumCloud.C：
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "MomentumCloud.H"
#include "integrationScheme.H"
#include "interpolation.H"
#include "subCycleTime.H"
#include "InjectionModelList.H"
#include "DispersionModel.H"
#include "PatchInteractionModel.H"
#include "StochasticCollisionModel.H"
#include "SurfaceFilmModel.H"
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
template<class CloudType>
void Foam::MomentumCloud<CloudType>::setModels()
{
dispersionModel_.reset
(
DispersionModel<MomentumCloud<CloudType>>::New
(
subModelProperties_,
*this
).ptr()
);
patchInteractionModel_.reset
(
PatchInteractionModel<MomentumCloud<CloudType>>::New
(
subModelProperties_,
*this
).ptr()
);
stochasticCollisionModel_.reset
(
StochasticCollisionModel<MomentumCloud<CloudType>>::New
(
subModelProperties_,
*this
).ptr()
);
surfaceFilmModel_.reset
(
SurfaceFilmModel<MomentumCloud<CloudType>>::New
(
subModelProperties_,
*this
).ptr()
);
UIntegrator_.reset
(
integrationScheme::New
(
"U",
solution_.integrationSchemes()
).ptr()
);
}
template<class CloudType>
template<class TrackCloudType>
void Foam::MomentumCloud<CloudType>::solve
(
TrackCloudType& cloud,
typename parcelType::trackingData& td
)
{
if (solution_.steadyState())
{
cloud.storeState();
cloud.preEvolve();
evolveCloud(cloud, td);
if (solution_.coupled())
{
cloud.relaxSources(cloud.cloudCopy());
}
}
else
{
/////////////////////////////////
Info<< "Dyfluid: Evolve function in MomentumCloud.C" << nl;
List<DynamicList<label>> PartLabelPre(U_.size());
List<DynamicList<label>> PartLabelPost(U_.size());
DynamicList<label> cellWithPartPre;
DynamicList<label> cellWithPartPost;
pNumber_.primitiveFieldRef() = 0.0;
forAllIter(typename MomentumCloud<CloudType>, *this, iter)
{
parcelType& p = iter();
pNumber_[p.cell()] += p.nParticle();
PartLabelPre[p.cell()].append(p.origId());
}
forAll(U_, cell)
{
if (PartLabelPre[cell].size() != 0)
{
cellWithPartPre.append(cell);
}
}
//Info<< cellWithPartPre << nl;
/////////////////////////////////
cloud.preEvolve();
evolveCloud(cloud, td);
if (solution_.coupled())
{
cloud.scaleSources();
}
/////////////////////////////////
forAllIter(typename MomentumCloud<CloudType>, *this, iter)
{
parcelType& p = iter();
PartLabelPost[p.cell()].append(p.origId());
}
const volScalarField& epsi = U_.mesh().lookupObject<volScalarField>("epsilon.water");
scalar C1 = 0.00481;
scalar C2 = 0.08;
scalar sigma = 0.072;
scalar rhoc = 998.0;
scalar D1 = 0.88;
scalar D2 = 9e6;
scalar muc = 1e-3;
const scalar dt = this->mesh().time().deltaTValue();
forAll(U_, cell)
{
// if there are particles in some cell
if (PartLabelPost[cell].size() != 0)
{
cellWithPartPost.append(cell);
scalar allVolume = 0.0;
scalar allnParticle = 0.0;
scalar allnParticleNew = 0.0;
scalar diamAve = 0.0;
//Info<< "Cell [" << cell << "] has " << PartLabelPost[cell].size()
// << " particles" << nl;
// loop all over particles
forAllIter(typename MomentumCloud<CloudType>, *this, iter)
{
parcelType& p = iter();
// loop particles label in this cell
forAll(PartLabelPost[cell], i)
{
// if this particle belongs to this cell
if (p.origId() == PartLabelPost[cell](i))
{
allVolume += p.nParticle()*M_PI/6.0*pow3(p.d());
allnParticle += p.nParticle();
diamAve = pow(allVolume/allnParticle*6.0/M_PI, 1.0/3.0);
//Info<< " Labels are " << PartLabelPost[cell](i)
// << ", its nParticle is " << p.nParticle()
// << ", volume is " << allVolume << nl;
}
}
}
const scalar d = diamAve;
// cell manipulation, calculate average diameter
if (d > SMALL)
{
const scalar epsilon = epsi[cell];
scalar gN =
allnParticle*C1*pow(epsilon, 1.0/3.0)/pow(d, 2.0/3.0)
*exp(-C2*sigma/(rhoc*pow(epsilon, 2.0/3.0)*pow(d, 5.0/3.0)));
scalar aSqrN =
sqr(allnParticle)*D1*pow(epsilon, 1.0/3.0)*4.0*sqrt(2.0)*pow(d, 7.0/6.0)
*exp(-D2*muc*rhoc*epsilon/sqr(sigma)*pow4(d)/16.0);
allnParticleNew = allnParticle + (gN - aSqrN)*dt;
diamAve = diamAve*pow(allnParticle/allnParticleNew, 1.0/3.0);
//Info<< "gN = " << gN << nl;
//Info<< "aSqrN = " << aSqrN << nl;
//Info<< "Particle increased by " << allnParticleNew - allnParticle << nl;
//Info<< "diamAve is " << diamAve << nl;
}
scalar newVolume = 0.0;
forAllIter(typename MomentumCloud<CloudType>, *this, iter)
{
parcelType& p = iter();
forAll(PartLabelPost[cell], i)
{
if (p.origId() == PartLabelPost[cell](i))
{
p.nParticle() = allnParticleNew/PartLabelPost[cell].size();
p.d() = diamAve;
newVolume += p.nParticle()*M_PI/6.0*pow3(p.d());
//Info<< " after brecoa, volume is " << newVolume
// << "p.nParticle is " << p.nParticle() << nl;
}
}
}
}
}
/////////////////////////////////
}
cloud.info();
cloud.postEvolve();
if (solution_.steadyState())
{
cloud.restoreState();
}
}
template<class CloudType>
void Foam::MomentumCloud<CloudType>::buildCellOccupancy()
{
if (cellOccupancyPtr_.empty())
{
cellOccupancyPtr_.reset
(
new List<DynamicList<parcelType*>>(this->mesh().nCells())
);
}
else if (cellOccupancyPtr_().size() != this->mesh().nCells())
{
// If the size of the mesh has changed, reset the
// cellOccupancy size
cellOccupancyPtr_().setSize(this->mesh().nCells());
}
List<DynamicList<parcelType*>>& cellOccupancy = cellOccupancyPtr_();
forAll(cellOccupancy, cO)
{
cellOccupancy[cO].clear();
}
forAllIter(typename MomentumCloud<CloudType>, *this, iter)
{
cellOccupancy[iter().cell()].append(&iter());
}
}
template<class CloudType>
void Foam::MomentumCloud<CloudType>::updateCellOccupancy()
{
// Only build the cellOccupancy if the pointer is set, i.e. it has
// been requested before.
if (cellOccupancyPtr_.valid())
{
buildCellOccupancy();
}
}
template<class CloudType>
template<class TrackCloudType>
void Foam::MomentumCloud<CloudType>::evolveCloud
(
TrackCloudType& cloud,
typename parcelType::trackingData& td
)
{
if (solution_.coupled())
{
cloud.resetSourceTerms();
}
if (solution_.transient())
{
label preInjectionSize = this->size();
this->surfaceFilm().inject(cloud);
// Update the cellOccupancy if the size of the cloud has changed
// during the injection.
if (preInjectionSize != this->size())
{
updateCellOccupancy();
preInjectionSize = this->size();
}
injectors_.inject(cloud, td);
// Assume that motion will update the cellOccupancy as necessary
// before it is required.
cloud.motion(cloud, td);
stochasticCollision().update(td, solution_.trackTime());
}
else
{
// this->surfaceFilm().injectSteadyState(cloud);
injectors_.injectSteadyState(cloud, td, solution_.trackTime());
CloudType::move(cloud, td, solution_.trackTime());
}
}
template<class CloudType>
void Foam::MomentumCloud<CloudType>::postEvolve()
{
Info<< endl;
if (debug)
{
this->writePositions();
}
this->dispersion().cacheFields(false);
forces_.cacheFields(false);
functions_.postEvolve();
solution_.nextIter();
if (this->db().time().writeTime())
{
outputProperties_.writeObject
(
IOstream::ASCII,
IOstream::currentVersion,
this->db().time().writeCompression(),
true
);
}
}
template<class CloudType>
void Foam::MomentumCloud<CloudType>::cloudReset(MomentumCloud<CloudType>& c)
{
CloudType::cloudReset(c);
rndGen_ = c.rndGen_;
forces_.transfer(c.forces_);
functions_.transfer(c.functions_);
injectors_.transfer(c.injectors_);
dispersionModel_.reset(c.dispersionModel_.ptr());
patchInteractionModel_.reset(c.patchInteractionModel_.ptr());
stochasticCollisionModel_.reset(c.stochasticCollisionModel_.ptr());
surfaceFilmModel_.reset(c.surfaceFilmModel_.ptr());
UIntegrator_.reset(c.UIntegrator_.ptr());
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class CloudType>
Foam::MomentumCloud<CloudType>::MomentumCloud
(
const word& cloudName,
const volScalarField& rho,
const volVectorField& U,
const volScalarField& mu,
const dimensionedVector& g,
const bool readFields
)
:
CloudType(cloudName, rho, U, mu, g, false),
cloudCopyPtr_(nullptr),
particleProperties_
(
IOobject
(
cloudName + "Properties",
this->mesh().time().constant(),
this->mesh(),
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
),
outputProperties_
(
IOobject
(
cloudName + "OutputProperties",
this->mesh().time().timeName(),
"uniform"/cloud::prefix/cloudName,
this->mesh(),
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
)
),
solution_(this->mesh(), particleProperties_.subDict("solution")),
constProps_(particleProperties_),
subModelProperties_
(
particleProperties_.subOrEmptyDict("subModels", true)
),
rndGen_(0),
cellOccupancyPtr_(),
cellLengthScale_(mag(cbrt(this->mesh().V()))),
pNumber_
(
IOobject
(
"pNumbers",
this->mesh().time().timeName(),
this->mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
this->mesh(),
dimensionedScalar(dimless, 0.0)
),
rho_(rho),
U_(U),
mu_(mu),
g_(g),
pAmbient_(0.0),
forces_
(
*this,
this->mesh(),
subModelProperties_.subOrEmptyDict
(
"particleForces",
true
),
true
),
functions_
(
*this,
particleProperties_.subOrEmptyDict("cloudFunctions"),
true
),
injectors_
(
subModelProperties_.subOrEmptyDict("injectionModels"),
*this
),
dispersionModel_(nullptr),
patchInteractionModel_(nullptr),
stochasticCollisionModel_(nullptr),
surfaceFilmModel_(nullptr),
UIntegrator_(nullptr),
UTrans_
(
new volVectorField::Internal
(
IOobject
(
this->name() + ":UTrans",
this->db().time().timeName(),
this->db(),
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
this->mesh(),
dimensionedVector(dimMass*dimVelocity, Zero)
)
),
UCoeff_
(
new volScalarField::Internal
(
IOobject
(
this->name() + ":UCoeff",
this->db().time().timeName(),
this->db(),
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
this->mesh(),
dimensionedScalar( dimMass, 0)
)
)
{
setModels();
if (readFields)
{
parcelType::readFields(*this);
this->deleteLostParticles();
}
if (solution_.resetSourcesOnStartup())
{
resetSourceTerms();
}
}
template<class CloudType>
Foam::MomentumCloud<CloudType>::MomentumCloud
(
const word& cloudName,
const volScalarField& rho,
const volVectorField& U,
const dimensionedVector& g,
const fluidThermo& carrierThermo,
const bool readFields
)
:
MomentumCloud(cloudName, rho, U, carrierThermo.mu(), g, readFields)
{}
template<class CloudType>
Foam::MomentumCloud<CloudType>::MomentumCloud
(
MomentumCloud<CloudType>& c,
const word& name
)
:
CloudType(c, name),
cloudCopyPtr_(nullptr),
particleProperties_(c.particleProperties_),
outputProperties_(c.outputProperties_),
solution_(c.solution_),
constProps_(c.constProps_),
subModelProperties_(c.subModelProperties_),
rndGen_(c.rndGen_),
cellOccupancyPtr_(nullptr),
cellLengthScale_(c.cellLengthScale_),
pNumber_(c.pNumber_),
rho_(c.rho_),
U_(c.U_),
mu_(c.mu_),
g_(c.g_),
pAmbient_(c.pAmbient_),
forces_(c.forces_),
functions_(c.functions_),
injectors_(c.injectors_),
dispersionModel_(c.dispersionModel_->clone()),
patchInteractionModel_(c.patchInteractionModel_->clone()),
stochasticCollisionModel_(c.stochasticCollisionModel_->clone()),
surfaceFilmModel_(c.surfaceFilmModel_->clone()),
UIntegrator_(c.UIntegrator_->clone()),
UTrans_
(
new volVectorField::Internal
(
IOobject
(
this->name() + ":UTrans",
this->db().time().timeName(),
this->db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
c.UTrans_()
)
),
UCoeff_
(
new volScalarField::Internal
(
IOobject
(
name + ":UCoeff",
this->db().time().timeName(),
this->db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
c.UCoeff_()
)
)
{}
template<class CloudType>
Foam::MomentumCloud<CloudType>::MomentumCloud
(
const fvMesh& mesh,
const word& name,
const MomentumCloud<CloudType>& c
)
:
CloudType(mesh, name, c),
cloudCopyPtr_(nullptr),
particleProperties_
(
IOobject
(
name + "Properties",
mesh.time().constant(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
)
),
outputProperties_
(
IOobject
(
name + "OutputProperties",
this->mesh().time().timeName(),
"uniform"/cloud::prefix/name,
this->mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
)
),
solution_(mesh),
constProps_(),
subModelProperties_(dictionary::null),
rndGen_(0),
cellOccupancyPtr_(nullptr),
cellLengthScale_(c.cellLengthScale_),
pNumber_(c.pNumber_),
rho_(c.rho_),
U_(c.U_),
mu_(c.mu_),
g_(c.g_),
pAmbient_(c.pAmbient_),
forces_(*this, mesh),
functions_(*this),
injectors_(*this),
dispersionModel_(nullptr),
patchInteractionModel_(nullptr),
stochasticCollisionModel_(nullptr),
surfaceFilmModel_(nullptr),
UIntegrator_(nullptr),
UTrans_(nullptr),
UCoeff_(nullptr)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class CloudType>
Foam::MomentumCloud<CloudType>::~MomentumCloud()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class CloudType>
void Foam::MomentumCloud<CloudType>::setParcelThermoProperties
(
parcelType& parcel,
const scalar lagrangianDt
)
{
parcel.rho() = constProps_.rho0();
}
template<class CloudType>
void Foam::MomentumCloud<CloudType>::checkParcelProperties
(
parcelType& parcel,
const scalar lagrangianDt,
const bool fullyDescribed
)
{
const scalar carrierDt = this->mesh().time().deltaTValue();
parcel.stepFraction() = (carrierDt - lagrangianDt)/carrierDt;
if (parcel.typeId() == -1)
{
parcel.typeId() = constProps_.parcelTypeId();
}
}
template<class CloudType>
void Foam::MomentumCloud<CloudType>::storeState()
{
cloudCopyPtr_.reset
(
static_cast<MomentumCloud<CloudType>*>
(
clone(this->name() + "Copy").ptr()
)
);
}
template<class CloudType>
void Foam::MomentumCloud<CloudType>::restoreState()
{
cloudReset(cloudCopyPtr_());
cloudCopyPtr_.clear();
}
template<class CloudType>
void Foam::MomentumCloud<CloudType>::resetSourceTerms()
{
UTransRef().field() = Zero;
UCoeffRef().field() = 0.0;
}
template<class CloudType>
template<class Type>
void Foam::MomentumCloud<CloudType>::relax
(
DimensionedField<Type, volMesh>& field,
const DimensionedField<Type, volMesh>& field0,
const word& name
) const
{
const scalar coeff = solution_.relaxCoeff(name);
field = field0 + coeff*(field - field0);
}
template<class CloudType>
template<class Type>
void Foam::MomentumCloud<CloudType>::scale
(
DimensionedField<Type, volMesh>& field,
const word& name
) const
{
const scalar coeff = solution_.relaxCoeff(name);
field *= coeff;
}
template<class CloudType>
void Foam::MomentumCloud<CloudType>::relaxSources
(
const MomentumCloud<CloudType>& cloudOldTime
)
{
this->relax(UTrans_(), cloudOldTime.UTrans_(), "U");
this->relax(UCoeff_(), cloudOldTime.UCoeff_(), "U");
}
template<class CloudType>
void Foam::MomentumCloud<CloudType>::scaleSources()
{
this->scale(UTrans_(), "U");
this->scale(UCoeff_(), "U");
}
template<class CloudType>
void Foam::MomentumCloud<CloudType>::preEvolve()
{
// force calculation of mesh dimensions - needed for parallel runs
// with topology change due to lazy evaluation of valid mesh dimensions
label nGeometricD = this->mesh().nGeometricD();
Info<< "\nSolving " << nGeometricD << "-D cloud " << this->name() << endl;
this->dispersion().cacheFields(true);
forces_.cacheFields(true);
updateCellOccupancy();
pAmbient_ = constProps_.dict().template
lookupOrDefault<scalar>("pAmbient", pAmbient_);
functions_.preEvolve();
}
template<class CloudType>
void Foam::MomentumCloud<CloudType>::evolve()
{
if (solution_.canEvolve())
{
typename parcelType::trackingData td(*this);
solve(*this, td);
}
}
template<class CloudType>
template<class TrackCloudType>
void Foam::MomentumCloud<CloudType>::motion
(
TrackCloudType& cloud,
typename parcelType::trackingData& td
)
{
CloudType::move(cloud, td, solution_.trackTime());
updateCellOccupancy();
}
template<class CloudType>
void Foam::MomentumCloud<CloudType>::patchData
(
const parcelType& p,
const polyPatch& pp,
vector& nw,
vector& Up
) const
{
p.patchData(nw, Up);
Up /= p.mesh().time().deltaTValue();
// If this is a wall patch, then there may be a non-zero tangential velocity
// component; the lid velocity in a lid-driven cavity case, for example. We
// want the particle to interact with this velocity, so we look it up in the
// velocity field and use it to set the wall-tangential component.
if (!this->mesh().moving() && isA<wallPolyPatch>(pp))
{
const label patchi = pp.index();
const label patchFacei = pp.whichFace(p.face());
// We only want to use the boundary condition value only if it is set
// by the boundary condition. If the boundary values are extrapolated
// (e.g., slip conditions) then they represent the motion of the fluid
// just inside the domain rather than that of the wall itself.
if (U_.boundaryField()[patchi].fixesValue())
{
const vector Uw1 = U_.boundaryField()[patchi][patchFacei];
const vector& Uw0 =
U_.oldTime().boundaryField()[patchi][patchFacei];
const scalar f = p.currentTimeFraction();
const vector Uw = Uw0 + f*(Uw1 - Uw0);
const tensor nnw = nw*nw;
Up = (nnw & Up) + Uw - (nnw & Uw);
}
}
}
template<class CloudType>
void Foam::MomentumCloud<CloudType>::updateMesh()
{
updateCellOccupancy();
injectors_.updateMesh();
cellLengthScale_ = mag(cbrt(this->mesh().V()));
}
template<class CloudType>
void Foam::MomentumCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
{
Cloud<parcelType>::autoMap(mapper);
updateMesh();
}
template<class CloudType>
void Foam::MomentumCloud<CloudType>::info()
{
vector linearMomentum = linearMomentumOfSystem();
reduce(linearMomentum, sumOp<vector>());
scalar linearKineticEnergy = linearKineticEnergyOfSystem();
reduce(linearKineticEnergy, sumOp<scalar>());
Info<< "Cloud: " << this->name() << nl
<< " Current number of parcels = "
<< returnReduce(this->size(), sumOp<label>()) << nl
<< " Current mass in system = "
<< returnReduce(massInSystem(), sumOp<scalar>()) << nl
<< " Linear momentum = "
<< linearMomentum << nl
<< " |Linear momentum| = "
<< mag(linearMomentum) << nl
<< " Linear kinetic energy = "
<< linearKineticEnergy << nl;
injectors_.info(Info);
this->surfaceFilm().info(Info);
this->patchInteraction().info(Info);
}
// ************************************************************************* //
