OpenFOAM 拉格朗日库深度解析

拉格朗日库类解析

类的结构

喷雾子模型

KinematicParcel

Kinematic parcel class with rotational motion (as spherical particles only) and one/two-way coupling with the continuous phase.

  • InjectionModel
    • CellZoneInjection
    • ConeInjection
    • ConeNozzleInjection
template<class CloudType>
class InjectionModelList
:
public PtrList<InjectionModel<CloudType>>

ConeNozzleInjection 构造函数:

tanVec1_ = normalised(perpendicular(direction_));
tanVec2_ = normalised(direction_ ^ tanVec1_);

ConeNozzleInjection::setPositionAndCell

scalar beta = mathematical::twoPi*rndGen.globalScalar01();
normal_ = tanVec1_*cos(beta) + tanVec2_*sin(beta);

ConeNozzleInjection::setProperties

scalar alpha = sin(coneAngle);
scalar dcorr = cos(coneAngle);

vector normal = alpha*normal_;
vector dirVec = dcorr*direction_;
dirVec += normal;
dirVec /= mag(dirVec);//dirVec是速度方向向量

ConeNozzleInjection::setProperties 函数的作用是设置出射粒子的速度和直径:

parcel.d() = sizeDistribution_->sample();

ConeNozzleInjection::setProperties 函数在 InjectionModel::inject 中被调用。
粒子直径在这里给定,然后在哪里更新呢?我知道的有破碎模型里边。在 SprayParcel::calc 中调用了破碎模型的 updatethis->d() 作为引用参数传入 update 函数,被破碎模型修改。

InjectionModel::prepareForNextTimeStep

// Number of parcels to inject
newParcels = this->parcelsToInject(t0, t1);
template<class CloudType>
Foam::label Foam::ConeNozzleInjection<CloudType>::parcelsToInject
//计算该时间步长内喷射的 parcels 数
(
const scalar time0,
const scalar time1
)
{
if ((time0 >= 0.0) && (time0 < duration_))
{
return floor((time1 - time0)*parcelsPerSecond_);
}
else
{
return 0;
}
}

KinematicCloud 中:

//- Injector models
InjectionModelList<KinematicCloud<CloudType>> injectors_;

//- construct
injectors_
(
subModelProperties_.subOrEmptyDict("injectionModels"),
*this
)

template<class CloudType>
inline const Foam::InjectionModelList<Foam::KinematicCloud<CloudType>>&
Foam::KinematicCloud<CloudType>::injectors() const
{
return injectors_;
}

injectors_.inject(cloud, td);

injectors_.updateMesh();

injectors_.info(Info);

最重要的 inject 函数在 KinematicCloud::evolveCloud 函数被调用。

distributionModel

哪里被用到?
答:InjectionModel 的派生类中,比如 ConeNozzleInjection.H 中:
类外申明:

// Forward declaration of classes
class distributionModel;

类内申明:

const autoPtr<distributionModel> sizeDistribution_;

ConeNozzleInjection.C 中,先 #include "distributionModel.H",然后在构造函数中:

sizeDistribution_
(
distributionModel::New
(
this->coeffDict().subDict("sizeDistribution"),
owner.rndGen()
)
)

这其实是很好的一个 RTS 的简单使用典范,可以参考这个,制作自己的 RTS。

  • particleForces
    • SphereDrag
  • patchInteractionModel
    • localInteraction
    • rebound
    • standardWallInteraction
    • none

ThermoParcel

Thermodynamic parcel class with one/two-way coupling with the continuous phase. Includes Kinematic parcel sub-models, plus:

  • heatTransferModel
    • RanzMarshall
  • surfaceFilmModel
    • thermoSurfaceFilm

ReactingParcel

Reacting parcel class with one/two-way coupling with the continuous phase.

  • compositionModel
    • SinglePhaseMixture
      算例文件 sprayCloudProperties:
compositionModel SinglePhaseMixture;

singlePhaseMixtureCoeffs
{
phases
(
liquid
{
C7H16 0.39;
C8H18 0.61;
}
);
}
```
注意这里给定的混合燃料的质量分数。
`CompositionModel` 中:
```C++
phasePropertiesList phaseProps_;

template<class CloudType>
Foam::CompositionModel<CloudType>::CompositionModel
(
const dictionary& dict,
CloudType& owner,
const word& type
)
:
CloudSubModelBase<CloudType>(owner, dict, typeName, type),
thermo_(owner.thermo()),
phaseProps_
(
this->coeffDict().lookup("phases"),
thermo_.carrier().species(),
thermo_.liquids().components(),
thermo_.solids().components()
)
{}

phasePropertiesList

List<phaseProperties> props_;

Foam::phasePropertiesList::phasePropertiesList
(
Istream& is,
const wordList& gasNames,
const wordList& liquidNames,
const wordList& solidNames
)
:
props_(is),
phaseTypeNames_(),
stateLabels_()
{
forAll(props_, i)
{
props_[i].reorder(gasNames, liquidNames, solidNames);
}

phaseTypeNames_.setSize(props_.size());
stateLabels_.setSize(props_.size());
forAll(props_, i)
{
phaseTypeNames_[i] = props_[i].phaseTypeName();
stateLabels_[i] = props_[i].stateLabel();
}
}

phaseProperties

Foam::phaseProperties::phaseProperties(Istream& is)
:
phase_(UNKNOWN),
stateLabel_("(unknown)"),
names_(0),
Y_(0),
carrierIds_(0)
{
is.check("Foam::phaseProperties::phaseProperties(Istream& is)");

dictionaryEntry phaseInfo(dictionary::null, is);

phase_ = phaseTypeNames[phaseInfo.keyword()];
stateLabel_ = phaseToStateLabel(phase_);

if (phaseInfo.size() > 0)
{
label nComponents = phaseInfo.size();
names_.setSize(nComponents, "unknownSpecie");
Y_.setSize(nComponents, 0.0);
carrierIds_.setSize(nComponents, -1);

label cmptI = 0;
forAllConstIter(IDLList<entry>, phaseInfo, iter)
{
names_[cmptI] = iter().keyword();
Y_[cmptI] = readScalar(phaseInfo.lookup(names_[cmptI]));
cmptI++;
}

checkTotalMassFraction();
}
}

这个 CompositionModel 在哪里使用呢?
首先 compositionModel_ReactingCloud 中定义,并提供了一个接口函数 composition()
ReactingParcelcalcPhaseChange 函数中,如下片段中被用到:

const CompositionModel<reactingCloudType>& composition =
cloud.composition();

scalarField X(composition.liquids().X(Y));

forAll(dMassPC, i)
{
const label cid = composition.localToCarrierId(idPhase, i);

const scalar dh = phaseChange.dh(cid, i, td.pc(), Tdash);
Sh -= dMassPC[i]*dh/dt;
}

forAll(dMassPC, i)
{
const label cid = composition.localToCarrierId(idPhase, i);

const scalar Cp = composition.carrier().Cp(cid, td.pc(), Tsdash);
const scalar W = composition.carrier().W(cid);
const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);

const scalar Dab =
composition.liquids().properties()[i].D(td.pc(), Tsdash, Wc);

// Molar flux of species coming from the particle (kmol/m^2/s)
N += Ni;

// Sum of Ni*Cpi*Wi of emission species
NCpW += Ni*Cp*W;

// Concentrations of emission species
Cs[cid] += Ni*d/(2.0*Dab);
}

SprayParcelcalc 函数中,如下片段中被用到:

const CompositionModel<reactingCloudType>& composition =
cloud.composition();

scalarField X0(composition.liquids().X(this->Y()));

scalar TMax = composition.liquids().Tc(X0);

if (composition.liquids().pv(pc0, T0, X0) >= pc0*0.999)
{
// Set TMax to boiling temperature
TMax = composition.liquids().pvInvert(pc0, X0);
}

......

SprayParcelcalcAtomization 函数中,如下片段中被用到:

scalar chi = 0.0;
if (atomization.calcChi())
{
chi = this->chi(cloud, td, composition.liquids().X(this->Y()));
}
  • phaseChangeModel
    • liquidEvaporationBoil
    • liquidEvaporation

spray

  • atomizationModel
    • BlobsSheetAtomization
    • LISAAtomization
  • breakupModel
    • ETAB
    • PilchErdman
    • ReitzDiwakar
    • ReitzKHRT
    • SHF
    • TAB
  • stochasticCollisionModel
    • ORourke
    • trajectory
    • none

turbulence

  • dispersionModel
    • stochasticDispersionRAS
    • GradientDispersionRAS

从求解器到库的调用过程

solver 中:

basicSprayCloud parcels
(
"sprayCloud",
rho,
U,
g,
slgThermo
);

parcels.evolve();

其中:

typedef SprayCloud
<
ReactingCloud
<
ThermoCloud
<
KinematicCloud
<
Cloud
<
basicSprayParcel
>
>
>
>
> basicSprayCloud;

typedef SprayParcel
<
ReactingParcel
<
ThermoParcel
<
KinematicParcel
<
particle
>
>
>
> basicSprayParcel;

我编写的函数调用过程:

SprayCloud::evolve 里边调用了 KinematicCloud::solve
KinematicCloud::solve 里边调用了:

ThermoCloud::preEvolve 里边调用了 KinematicCloud::preEvolve
KinematicCloud::evolveCloud 里边调用了:

ReactingCloud::resetSourceTerms 里边调用了:

InjectionModel::inject 里边调用了:

KinematicCloud::motion里边调用了 Cloud::move
Cloud::move 又调用了 KinematicParcel::move

KinematicParcel::move 里边调用了:

SprayParcel::setCellValues 里边调用了:

SprayParcel::calc 里边调用了:

ReactingParcel::calc 里边调用了:

SprayCloud::info 里边调用了:

经过验证,ThermoParcel::calc 函数从来没有被调用!
KinematicParcel::calc 函数也从来没有被调用!

move 定义于 KinematicParcelCloud 中。
Cloud 中的 move 一个时间步长内只调用一次。

Added 10 new parcels,对应 KinematicParcel::move 就调用 10 次,这里是在 InjectionModel::inject 中被调用的。后来如果破碎了,parcels 数变成 20 了,那么 move 函数又调用 20 次,这里是在 KinematicCloud::motion 中被调用的。

KinematicParcel::move 调用了 SprayParcel::calc,并且每个 move 函数内,calc 调用不止一次。
SprayParcel::calc中调用了喷雾的 calcAtomizationcalcBreakup,因此这个里边,粒子直径被修改了,粒子个数也可能被修改。

以下 UML 图来自openfoamwiki,基于OF2.3.X,可供参考。




喷雾信息输出到info:

Solving 3-D cloud sprayCloud

Cloud: sprayCloud injector: model1
Added 18 new parcels

Cloud: sprayCloud
Current number of parcels = 18
Current mass in system = 1.19539e-11
Linear momentum = (-1.22046e-12 2.17761e-13 -1.86459e-11)
|Linear momentum| = 1.86871e-11
Linear kinetic energy = 2.07355e-11
model1:
number of parcels added = 18
mass introduced = 1.19539e-11
Parcel fate (number, mass)
- escape = 0, 0
- stick = 0, 0
Temperature min/max = 358.01, 358.026
Mass transfer phase change = 5.38154e-18
D10, D32, Dmax (mu) = 178.438, 274.262, 337.897
Liquid penetration 95% mass (m) = 1.59326e-06

代码中的调用顺序:
KinematicCloud::solve 中:

cloud.info();

这里的 cloud 类型是模板 TrackCloudType。为了搞清楚它的真实面目,我们修改一下这里的代码,然后编译(这里必须执行 Lagrangian 库下的 Allwmake,然后编译求解器)。
后来发现好像只需要编译求解器就行了,因为那些是模板,在求解器中才实例化?
添加下面的代码到任意位置:

cloud.abc();

编译库是并不会报错,编译求解器时报错如下:

src/lagrangian/intermediate/lnInclude/KinematicCloud.C(97): error: class "Foam::SprayCloud<Foam::ReactingCloud<Foam::ThermoCloud<Foam::KinematicCloud<Foam::Cloud<Foam::basicSprayParcel>>>>>" has no member "abc"
cloud.abc();

所以这里的 cloudSprayCloud,它的 info

template<class CloudType>
void Foam::SprayCloud<CloudType>::info()
{
CloudType::info();
scalar d32 = 1.0e+6*this->Dij(3, 2);
scalar d10 = 1.0e+6*this->Dij(1, 0);
scalar dMax = 1.0e+6*this->Dmax();
scalar pen = this->penetration(0.95);

Info << " D10, D32, Dmax (mu) = " << d10 << ", " << d32
<< ", " << dMax << nl
<< " Liquid penetration 95% mass (m) = " << pen << endl;
}

其中

template<class CloudType>
inline Foam::scalar Foam::KinematicCloud<CloudType>::Dmax() const
{
scalar d = -great;
forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
{
const parcelType& p = iter();
d = max(d, p.d());
}

reduce(d, maxOp<scalar>());

return max(0.0, d);
}

SprayCloudCloudTypeReactingCloud

template<class CloudType>
void Foam::ReactingCloud<CloudType>::info()
{
CloudType::info();

this->phaseChange().info(Info);
}

ReactingCloudCloudTypeThermoCloud

template<class CloudType>
void Foam::ThermoCloud<CloudType>::info()
{
CloudType::info();

Info<< " Temperature min/max = " << Tmin() << ", " << Tmax()
<< endl;
}

ThermoCloudCloudTypeKinematicCloud

template<class CloudType>
void Foam::KinematicCloud<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);
}
文章作者: Yan Zhang
文章链接: https://openfoam.top/lagrangian/
版权声明: 本博客所有文章除特别声明外,均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 OpenFOAM 成长之路
微信打赏给博主更多动力吧~