MPPICFOAM 源码剖析
往期相关文章
模型简介
固相应力
数值实现
1
求解器主程序MPPICFOAM
MPPIC和DPM共用一个求解器文件，仅在定义cloud时存在差异，主程序会根据宏#ifdef MPPIC来选择使用basicKinematicTypeCloud的具体定义。根据此宏可以将粒子云定义为basicKinematicMPPICCloud或者basicKinematicCollidingCloud。这两种粒子云的差别在于计算颗粒云内粒子作用力模型的实现，简单来说，MPPIC会使用颗粒应力表示颗粒间作用力，而collidingCloud会用碰撞模型来计算作用力。
两个求解器MPPICFoam和DPMFoan仅在粒子云上存在差别，而整个算法的实现流程都是一样的，主要关注两个粒子云的定义，仅根据粒子的实现就能实现两个风格差异很大的求解器，说明了类继承和模板的强大，进一步说明了代码可复用的重要性。
basicKinematicMPPICCloud的定义
typedef MPPICCloud<KinematicCloud<Cloud<basicKinematicMPPICParcel>>> basicKinematicMPPICCloud;
KinematicCloud类模板主管颗粒云的运动，而Cloud主管颗粒云的整体演变。
从MPPICCloud开始解析看类定义template<classCloudType> class MPPICCloud:
public CloudType MPPICCloud是一个类模板，并且其模板参数为其父类。在这里也就是KinematicCloud。
Protected成员包含三个量，分别是packingModel_，dampingModel_和isotropyModel_三个模型。
之后是一些构造和用于访问成员的内联函数，成员函数中重要的主要是evolve函数和motion函数。
2
MPPICFOAM颗粒相求解过程
颗粒相求解过程：
01
02
03
04
接下来，主要讲解使用MPPICFOAM如何实现以上4个过程：
求解器程序初始化后进入主循环中，第一步就是颗粒相演变，kinematicCloud.evolve()。kinematicCloud是之前定义的basicKinematicMPPICCloud类变量，负责粒子区信息的管理，从此函数开始进行主要内容的解析。
在evolve函数中先定义了一个TrackingData的变量，此变量是建立Cloud和parcel之间信息传递的桥梁。这里的parcel就是basicKinematicMPPICParcel，一个MPPICParcel的实例。之后调用this->solve(td)，这个solve函数在MPPICCloud中是没有定义的，只能是继承自父类kinematicCloud。
进入solve，主要有三个函数：
td.cloud().preEvolve()
通过调用kinematicCloud，主要更新一些力和粒子-网格信息，cellOccupancy[iter().cell()].append(&iter())将粒子归入到每个网格内。
evolveCloud(td)
evolveCloud主要处理的粒子移动，调用的是MPPICCloud的motion函数，在这之前如果有inject，还会更新粒子的占据cell信息。
td.cloud().scaleSources()
通过传入的td参数访问MPPICCloud
在MPPICCloud.motion中定义td的track类型，可以影响之后具体的parcel运动。
td.part() = TrackData::tpLinearTrack;
CloudType::move(td,this->db().time().deltaTValue());
将td的track类型定为tpLinearTrack，随后调用了CloudType::move，这里CloudType为kinematicCloud。查找kinematicCloud的move函数，发现kinematicCloud没有重载move函数，因此调用的是Cloud的move。在此函数内有并行处理。
假设一切条件合适后，会遇到一个遍历所有parcel的循环。
forAllIter(typenameCloud<ParticleType>, *this, pIter)
{
ParticleType& p = pIter();
// Move the particle
bool keepParticle = p.move(td,trackTime);
在此循环中，根据每个具体颗粒调用对应的颗粒类的move函数（在这里为MPPICParcel）。也就是说，这些Cloud不负责具体的力和速度、位置计算，而是从全局的角度考虑整个粒子群的力学模型和管理所有粒子而已。
在MPPICParcel的move中可以找到
case TrackData::tpLinearTrack:
{
ParcelType::move(td, trackTime);
break;
注意到之前传入td.part()的正好是tpLinearTrack，因此执行这一分支语句，即颗粒正常移动位置（ParcelType对应KinematicParcel），同时在p.move中计算p.stepFraction()，即此次移动耗时百分比。
执行完p.move后，根据粒子是否运动出计算域来决定颗粒去留，之后仍是并行计算，即所有的颗粒并行计算是在Cloud内进行的。
接下来是MPPIC的三个独有模型，计算它们对颗粒速度的修正。
首先是damping模型：
//Damping
if (dampingModel_->active())
{
td.updateAverages(*this);
dampingModel_->cacheFields(true);
td.part() =TrackData::tpDampingNoTrack;
CloudType::move(td,this->db().time().deltaTValue());
td.part() = TrackData::tpCorrectTrack;
CloudType::move(td,this->db().time().deltaTValue());
dampingModel_->cacheFields(false);
}
inlineTrackingData
(
CloudType& cloud,
trackPart part = tpLinearTrack
);
trackpart采用默认参数，其构造在MPPICParcelTrackingDataI.H文件中，主要是用AveragingMethod类对象进行初始化volumeAverage_，radiusAverage_，rhoAverage_，uAverage_，uSqrAverage_，frequencyAverage_和massAverage_，即颗粒的体积，半径，密度，速度，速度方差，频率和质量。AveragingMethod是一个抽象基类，可以通过RTS选用不同的平均方法。
进入MPPICParcel的td. updateAverages，根据具体的AveragingMethod，其对应的权重因子是不同的，因此需要计算权重因子，这是一个场量。先是遍历cloud内的所有粒子，直接在tet四面体内对体积，质量和密度，速度进行累加，用到了AveragingMethod的add方法，同样需要计算权重因子。在通过将计算的速度平均（cell）插值到粒子所在位置，计算速度方差。计算质量平均和sauter直径。为了计算碰撞频率，需要计算颗粒位置的平均值插值，如体积分数，半径和速度，并还要计算权重以计算平均值。
接下来，通过dampingModel_->cacheFields(true)将欧拉场读取到个场变量内。
通过改变td的part标示为tpDampingNoTrack来改变MPPICParcel的对应状态以便在MPPICParcel中具体处理parcel的移动。再次进入kinematicCloud也就是Cloud的move函数，继续再调用MPPICParcel的move对单个parcel移动，但注意此时的td标示已经改变。从字面来看tpDampingNoTrack标示只计算damping速度修正而不移动粒子。
caseTrackData::tpDampingNoTrack:
{
p.UCorrect() =
td.cloud().dampingModel().velocityCorrection(p, trackTime);
td.keepParticle = true;
通过调用MPPICCloud的damping模型对单个粒子进行修正，再次注意，Cloud负责模型，而parcel负责具体的移动，如速度和位置，damping模型可以进行RTS。
以relax模型为例，其修正速度代码为：
const scalar x =
deltaT*oneByTimeScaleAverage_->interpolate(p.position(), tetIs);
const vector u= uAverage_->interpolate(p.position(), tetIs);
return (u - p.U())*x/(x + 2.0);
其中oneByTimeScaleAverage_表示弛豫时间倒数
，这个时间是基于cell的，因此需要插值到parcel位置。于是x就为
，而u表示cell平均颗粒速度的颗粒位置插值
，最终返回的速度修正为：
这其中涉及到弛豫时间，也就是时间尺度的计算，而时间尺度在dampingModel_->cacheFields中得到计算，以及其他欧拉场量得到读取。
弛豫时间计算如下：
oneByTimeScaleAverage_() =
(
this->timeScaleModel_->oneByTau
(
volumeAverage,
radiusAverage,
uSqrAverage,
frequencyAverage
)
)();
其中用到了timeScaleModel_->oneByTau。timeScaleModel_也是基类，需要用RTS确定具体的子类。
假设为nonEquilibrium模型。此模型的初始化是在damping的构造函数中：
timeScaleModel_
(
TimeScaleModel::New
(
this->coeffDict().subDict(TimeScaleModel::typeName)
可以看到，是通过读取damping的dict的子字典TimeScaleModel::typeName也就是“timeScaleModel”得到。为了计算弛豫时间，进入到nonEquilibrium的oneByTau函数。
staticconst scalar a =
8.0*sqrt(2.0)/(3.0*constant::mathematical::pi)
*0.25*(1.0 - e_*e_);
return a*f*alphaPacked_/max(alphaPacked_ -alpha, SMALL);
其中：
最后返回：
回到MPPICCloud.motion的damping块，这里再次将td的标示设为tpCorrectTrack，通过move修正parcel的速度和位置。进入MPPICParcel：
case TrackData::tpCorrectTrack:
{
vector U = p.U();
scalar f = p.stepFraction();
p.U() = (1.0 - f)*p.UCorrect();
ParcelType::move(td, trackTime);
p.U() = U + (p.stepFraction() -f)*p.UCorrect();
break;
}
在上一步中，p.UCorrect()已经由damping模型得到更新（p是对当前粒子的引用），并且已经减去了最初的颗粒速度p.U()，因此通过(1.0 - f)*p.UCorrect()计算得到的p.U()实际为修正速度。
执行完damping模型计算后，在MPPICCloud.motion中进入packing模型：
// Packing
if (packingModel_->active())
{
// same procedure as for damping
td.updateAverages(*this);
packingModel_->cacheFields(true);
td.part() =TrackData::tpPackingNoTrack;
CloudType::move(td,this->db().time().deltaTValue());
td.part() = TrackData::tpCorrectTrack;
CloudType::move(td,this->db().time().deltaTValue());
packingModel_->cacheFields(false);
}
其执行流程与damping模型基本一样，不同的是packing模型的调用：同样，将td.part()设定为tpPackingNoTrack，并执行move来调用packing模型，走固定程序，直接进入MPPICParcel的move程序中的packing块：
caseTrackData::tpPackingNoTrack:
{
p.UCorrect() =
td.cloud().packingModel().velocityCorrection(p, trackTime);
可以看到在MPPICParcel通过td建立的联系调用了与cloud相关联的packingModel对速度进行修正。packing模型也可以进行RTS。
以Explicit形式为例，在velocityCorrection中先计算固相应力：
constvector tauGrad =
stressAverage_->interpolateGrad(p.position(), tetIs);
固相应力stressAverage_在cacheFields中计算：
stressAverage_()=
this->particleStressModel_->tau
(
*volumeAverage_,
rhoAverage,
uSqrAverage
)();
可以看到，其调用了应力模型，particleStress模型同样也有RTS，以HarrisCrighton为例。
的计算过程为：
return
(
pSolid_
* pow(alpha, beta_)
/ denominator(alpha)
);
其中denominator(alpha)为
max
(
alphaPacked_ - alpha,
max(eps_*(1.0 - alpha), SMALL)
);
为： 
计算完颗粒应力之后，接着就需要对计算的速度进行修正，于是有：
if ((uRelative & alphaGrad)> 0)
{
dU = - deltaT*tauGrad/(p.rho()*alpha);
}
// apply the velocity limiters
return
correctionLimiting_->limitedVelocity
(
p.U(),
dU,
uMean
);
其中，dU为从连续的颗粒固相应力梯度计算的颗粒速度
。之后，通过加速度方向以及颗粒的当前运动速度以及恢复系数等参数就能够对颗粒速度进行修正。上一篇文章中提到根据计算修正速度的值时使用的是相对速度
，颗粒速度
，CorrectionLimitingMethod同样可以有RTS。
以relative为例，其返回的修正速度为
const vector uRelative = uP -uMean;
return minMod
(
dU,
- (1.0 + this->e_)*uRelative
);
即对应
if (isotropyModel_->active())
{
td.updateAverages(*this);
// apply isotropy model
isotropyModel_->calculate();
}
Isotropy块不同于前两个块，其仅仅只调用isotropyModel的calculate()对颗粒速度进行计算，更新的速度将用于下一时间步。同样isotropyModel也有RTS，如图：
exponentAverage=
exp
(
- deltaT
*this->timeScaleModel_->oneByTau
(
volumeAverage,
radiusAverage,
uSqrAverage,
frequencyAverage
)
)();
只是这里调用的时间尺度模型为isotropic，其oneByTau函数返回值为
staticconst scalar a =
8.0*sqrt(2.0)/(5.0*constant::mathematical::pi)
*0.25*(3.0 - e_)*(1.0 + e_);
return a*f*alphaPacked_/max(alphaPacked_ - alpha, SMALL);
最后返回：
接着第二步为速度取样步，对颗粒云中的每个粒子生成一个随机数
// random sampling
forAllIter(typenameCloudType, this->owner(), iter)
{
typenameCloudType::parcelType& p = iter();
const tetIndicestetIs(p.cell(), p.tetFace(), p.tetPt(), mesh);
const scalar x =exponentAverage.interpolate(p.position(), tetIs);
if (x <rndGen.sample01<scalar>())
{
const vectorr(sampleGauss(), sampleGauss(), sampleGauss());
const vector u =uAverage.interpolate(p.position(), tetIs);
const scalar uRms=
sqrt(max(uSqrAverage.interpolate(p.position(), tetIs), 0.0));
p.U() = u +r*uRms*oneBySqrtThree;
}
}
即X = rndGen.sample01<scalar>()，
当通过尺度模型计算的exponentAverage等于：
满足
时，
新的颗粒速度为：
，
否则
其中
通过r(sampleGauss(), sampleGauss(), sampleGauss())计算，sampleGauss()为Stochastic私有函数。
//conservation correction
forAllIter(typename CloudType,this->owner(), iter)
{
typename CloudType::parcelType& p =iter();
const tetIndices tetIs(p.cell(),p.tetFace(), p.tetPt(), mesh);
const vector u =uAverage.interpolate(p.position(), tetIs);
const scalar uRms =
sqrt(max(uSqrAverage.interpolate(p.position(), tetIs), 0.0));
const vector uTilde = uTildeAverage.interpolate(p.position(),tetIs);
const scalar uTildeRms =
sqrt(max(uTildeSqrAverage.interpolate(p.position(), tetIs), 0.0));
p.U() = u + (p.U() -uTilde)*uRms/max(uTildeRms, SMALL);
}
其中，
= uTildeAverage为颗粒整体平均速度，
= uRms颗粒速度均方根，
= uTildeRms为颗粒整体平均速度均方根。最终，得到速度
。
