MP-PIC固相运动的数值实现
在此基础上，本期阐述MI-PIC模型的数值实现，重点介绍颗粒相运动求解过程，对固相应力模型进行调整的固相速度修正模型以及气固间插值和平均方法。
欧拉-拉格朗日方法中，颗粒相和连续相在计算尺度和描述方式的不同导致不能用与单一流体相同的方式进行求解，而是在各自的相空间分别求解。
MP-PIC方法气固场求解过程的四个阶段（图1）
a | 颗粒相信息集中 - 平均过程； |
b | 流体场和颗粒相连续场求解； |
c | 流体信息和颗粒连续信息插值到离散颗粒位置； |
d | 颗粒速度和位置更新。 |
图1 颗粒相求解的四个阶段
从图1可以看到：由于流体相的信息存储在计算网格体心，而颗粒相的计算粒子随时间变化在空间不断更换位置，为了在拉格朗日相空间和欧拉网格间传递颗粒特性来实现相间耦合，必须使用插值方法来完成此项操作。
通常插值权重函数为三个方向上插值算子的乘积，且这三个算子需要满足相互正交和有界。如x方向的算子[1]：
方程1
式中，xξ表示网格节点位置，xp表示计算颗粒位置。
一般来说，任何满足此条件的函数都可以用来计算权重，当使用笛卡尔网格-切割体网格，为了减少计算复杂度，通常采用三线性插值函数。例如，网格上的固相体积分数可以通过如下插值操作得到，
方程2
其中，Vp,i 是颗粒体积，np,i 是数值粒子中包含的颗粒数，Np是当前网格ξ内的数值粒子数量，Vξ为网格体积。此外，在更新颗粒速度时需要计算颗粒应力的梯度，对于梯度的计算如下式所示：
方程3
式中，Nξ 表示数值粒子周围的网格，即粒子位置的场梯度为周围网格场梯度的插值的和。插值函数的乘积和梯度的详细定义和运算可以参考文献[2]。
在MP-PIC中，将颗粒相信息分配到每个颗粒，让它们单独携带并在一个时间步内以彼此独立的特定属性运动。当进入下一个时间步后，先将网格内属于此网格或影响此网格的所有粒子信息收集，在欧拉网格上重新构建出连续的颗粒相场以便计算颗粒间作用力。因而可以将颗粒相既当做连续相又可以当做离散相。
图2为MP-PIC求解流程图，其中（a）为气固相求解，（b）为固相求解。
图2 MP-PIC求解流程图
新时间层的颗粒位置可以通过对方程积分得到更新：
方程4
颗粒速度则通过对加速度方程积分得到，对所有项都采用新时间层的量后有隐式形式
方程5
其中：
为颗粒位置的插值隐式速度；
为颗粒位置的插值隐式气相压力梯度；
为颗粒位置的插值隐式固相压力梯度。
、
和
为从欧拉网格中心点上计算得到的平均场量在数值粒子所处位置的插值。
图3为考虑了颗粒阻尼耗散和各项同性松弛对颗粒运动影响的示意图。在这其中颗粒应力模型通常是一个分段函数，在颗粒体积分数接近最大堆积浓度时，会得到间断解，容易导致数值不稳定。
图3 数值颗粒运动过程
由于固相应力的计算相对于其他几项对颗粒浓度的变化极为敏感，为了维持求解过程的稳定性，速度方程中的颗粒应力项的计算需要特殊处理。此外，当颗粒应力占主导时，由其计算的修正速度非常大而与实际不符，因而不能直接在方程5中使用。通常方程5中颗粒速度的计算被分割成两个部分，即速度为固相应力和其他所有作用力时间积分的总和，这样能为计算带来更大的稳定性。
对于某一特定粒子，
方程6
因此，排除了固相应力后计算的中间速度
为：
方程7
从连续的颗粒固相应力梯度计算的颗粒速度为：
方程8
为保障计算的稳定性，需要忽略分母中的隐式项，因而
重新写为：
方程9
1
固相速度修正
通过对固相应力模型的分析（详见：MP-PIC颗粒间作用力：固相应力）可以知道，颗粒应力模型计算出的速度将远大于颗粒运动速度，如果直接用于颗粒运动，将会违背守恒原理。而实际过程中，颗粒接近大量堆积颗粒时是逐渐减速后反射的。因此，应力模型的作用应该是防止颗粒进入已经处于堆积状态的颗粒群而出现过度堆积。同时，应力模型给出了颗粒所受应力的方向，即加速度方向。通过加速度方向以及颗粒的当前运动速度以及恢复系数等参数就能够对颗粒速度进行修正。
根据对浓度计算方式的不同，速度修正方法分为显式和隐式方法。
在更新颗粒速度的计算中，为了保持计算的稳定性，减小颗粒应力项对颗粒体积分数的敏感性，将颗粒速度的计算分割成了两部分。通过对前一部分的计算，可以对速度进行更新，并通过此中间速度得到颗粒云的中间位置
方程10
根据得到的中间颗粒位置，可以计算出中间颗粒浓度，并用此浓度近似下一时刻浓度从而计算颗粒应力和移动颗粒，并最终得到新的颗粒浓度分布
方程11
由前面分析可知，由于颗粒应力仅在颗粒接近堆积浓度时才变得重要，因而对于绝大部分区域，中间颗粒浓度基本上等于下一时间步的颗粒浓度。并在此实现中将显式地用于颗粒应力的计算进而计算修正速度。为了更为直观的介绍速度修正过程，假设颗粒仅在垂直方向运动，如图4所示。
图4 显式速度修正
图4展示了单个粒子相对于堆积区的三种典型的运动状态：
A为颗粒处于堆积区内部，运动受到周围粒子或壁面的限制，运动速度基本与堆积区颗粒的整体速度一致；
B为颗粒朝堆积区运动，将受到颗粒应力作用而逐渐减速，甚至会反弹；
C为颗粒远离堆积区，不再受颗粒应力作用。
方程12表示当应力模型计算的速度朝上时（仅在垂直方向），用于限制颗粒速度的速度修正方程。图4右侧的四种情况则更为详细地描述了如何确定修正速度，具体的实现过程可以参考Garg[3]。
方程12
隐式形式中不再用中间浓度来计算颗粒应力，而是直接计算出下一时间步的颗粒浓度。类似于双流体模型，对固相建立连续方程或者用颗粒控制方程得到：
方程13
将第二项用中间时间层的量近似，并分裂成两项：
方程14
之后认为第二项满足中间时间层的连续方程，并将第三项中的速度用速度应力方程代替：
方程15
最后对最后一项用链式法则，得到关于固相浓度的方程：
方程16
求得新时间层的固相浓度后就可以计算出单元颗粒速度。
相比于通常的欧拉-拉格朗日方法，MP-PIC方法颗粒相与流体网格结合得更为紧密，需要频繁地在两者之间交换信息。精确而高效的插值和平均算法是保证MP-PIC计算鲁棒性和一致性的前提，特别是在非结构网格中，由于丢失了与坐标轴线相关的拓扑信息，插值和平均算法需要更具有通用性。存储流体信息的欧拉网格是按一定规则组织的离散的空间点集，插值计算即意味着基于这些空间点的坐标和场值，以及假定的场变化来构建连续的场分布。对于结构化且网格线与坐标轴线共线的网格，有许多高阶、有效的插值格式，这些高精度插值格式通过在三个维度上进行重复操作即可实现将数据从欧拉网格传递到空间中任意点。在实际的气固流动模拟中，绝大多数的设备装置都具有复杂的几何结构。虽然也可以通过台阶网格（Step Grid）来近似边界，但对计算域仍会引入较大的误差。
相对来说，非结构网格对流场几何边界的保真程度相当高，能很好地用于流体相的数值模拟。然而，当涉及到气-固两相间插值操作时则会对计算过程以及程序开发带来了很大的困难，很多在结构网格上充分发展，相当成熟的高精度方法不能在非结构网格上使用。此外，在非结构网格上选择插值方法时，有几点要求需要确保，如全面覆盖流场、插值量在空间中连续变化、具有能恢复构造离散点值的能力、仅与几何相关以及强连通性。由于涉及到空间插值，就会伴随有空间分解，将复杂的多面体网格单元分解为基础的四面体能很方便的进行插值计算。虽然也可以在六面体等复杂的单元内插值，但无疑会极大地增加计算难度。本文涉及的非结构网格主要针对多面体网格（Polyhedronmesh）。在接下来的内容中将介绍单元分解方法以及相关的四面体内插值算法。
Cell-Vertex Method
Cell-Vertex方法主要基于将每个多面体网格单元分解成多个四面体单元，并在四面体单元内进行插值计算。
分解策略如下：
1
对于给定的多面体网格单元，初步分解成多个小多面体，每个小多面体由单元中心（Cell Centre）和每个面的所有顶点组成；
2
对于每个小多面体，基于原多面体的面再次进行分解，根据使多面体质量最优的原则，在面的多个顶点中选择一个点作为分解基点，将面上与基点不相邻点间连线，从而分割原来的面，最后用组成分割成的三角面的三个顶点与单元中心点组成四面体。
如图5所示，分解后得到的一个单元体由面基础点Base，构成一条边的A、B点以及网格中心点C构成。
图5 计算网格分解为四面体单元
将计算网格分解成由四面体组成的子网格系统后，此网格系统由单元中心点和各个单元顶点组成，其中单元中心点存储了所有的流体信息（同位网格）。由于插值过程在每个四面体内进行，因此在插值前需要知道单元顶点的流场值。这一步由单元中心到单元顶点的平均过程完成。单元顶点处的值通过对环绕此顶点的所有单元值进行平均，由于假设流场量在空间中线性变化，而各单元中心距此顶点距离不同，因而其影响也不尽相同，需要通过权重来表示。在这里使用中心距顶点的距离的倒数作为权重。计算表达式为：
方程17
方程18
其中Φv 和Φi 分别表示单元顶点和第个单元中心的场量，wi为归一化的权重值，ri 为从单元中心到顶点的距离。采用距离的倒数作为权重计算是众多计算方法中一种应用较广且比较简单的方式，此外，还可以通过拉格朗日乘子法（Lagrangemultipliers）来计算权重。
权重计算
在将整个计算域用上述任意一种方法分解后，会得到众多的四面体单元，这些单元相互不重叠，且完全充满整个空间。计算域中任意一点必然会位于其中某一四面体内或者单元中心、单元顶点甚至单元面上。这样，由于四面体结构的固定性，任意一个特定点上的流场值都可以通过四面体单元四个顶点上的值插值得到。这其中涉及到四面体内权重的计算，在这里使用基于重心的权重算法。
假设流场在空间连续线性地变化，因此就可以根据空间点的位置计算权重，这些权重值在颗粒位置的下一次更新前可以一直使用。考虑一个由四个顶点r1，r2，r3，r4组成的四面体，其内部任意点r都可以用四个顶点的组合表示，即
方程19
通过将r用分量形式表示，加上λ的和为1的条件，可以建立关于λ的代数方程组：
方程20
其中T 为3×3的矩阵，其具体的表达式为[4]：
方程21
计算得到的作为当前点r的权重值。
通过积分牛顿运动方程后，颗粒移动到新位置，因而使颗粒云的空间分布得到更新。从连续颗粒应力模型计算速度的方程可知，计算颗粒中间速度需要使用颗粒群整体信息的平均值，即拉格朗日场的平均表述形式。其位于气相流场的网格节点上，包括颗粒浓度、颗粒应力梯度和颗粒整体的平均速度以及这些场的梯度场。而粒子是离散的，因此需要将这些离散的信息收集起来建立连续的场。平均方法主要有单元平均，双重网格平均以及矩平均[5]，这里主要介绍双重网格平均。
双重网格平均在单元分解的基础上判定颗粒归属的四面体，对于单元中心点，颗粒将单元内所有颗粒信息都映射到单元中心上，而对于单元顶点则需要建立双重网格，在整个映射过程中，仍然使用基于四面体重心计算的权重系数。
方程22
方程23
其中Φ为需要重构的平均颗粒场，N为与中心或顶点相关的颗粒数，Vtet为颗粒所在四面体的体积，除以4是因为双重网格的原因。
由于假设构建的场在四面体内线性变化，因此，在计算过程中可以认为梯度在四面体内为常量，梯度计算采用有限元中的一阶梯度估算，即线性梯度重构。假设一个偏离四面体顶点x1的任意点x1+h，其值可以通过泰勒展开式近似为[6]：
方程24
因而，对其他三个顶点都用此式近似，得到如下3×3系统：
方程25
从而可以计算出当前四面体的颗粒平均场梯度▽Φ。
下 期 预 告
本篇介绍了颗粒相运动求解过程，对固相应力模型进行调整的固相速度修正模型以及气固间插值和平均方法。下期将对MPPIC模型在OpenFoam中的实现代码进行解析，主要包括其三大模型：
damping模型
packing模型
isotropy模型
参考文献
[1] O'Rourke, P.J., J.U. Brackbill, and B. Larrouturou, On Particle-Grid Interpolation and Calculating Chemistry in Particle-in-Cell Methods. Journal of Computational Physics, 1993. 109(1): p. 37-52.
[2] Snider, D.M., An Incompressible Three-Dimensional Multiphase Particle-in-Cell Model for Dense Particle Flows. Journal of Computational Physics, 2001. 170(2): p. 523-549.
[3] R. Garg, J.F.D., Documentation of open-source MFIX–PIC software for gas-solids flows. https://mfix.netl.doe.gov/documentation/mfix_pic_doc.pdf, 2013.
[4] Goldschmidt, M.J.V., R. Beetstra, and J.A.M. Kuipers, Hydrodynamic modelling of dense gas-fluidised beds: comparison and validation of 3D discrete particle and continuum models. Powder Technology, 2004. 142(1): p. 23-47.
[5] Kuipers, J.A.M., A two-fluid micro balance model for fluidized beds. 1990, University of Twente.
[6] Correa, C.D., R. Hero, and K.-L. Ma, Acomparison of gradient estimation methods for volume rendering on unstructuredmeshes. IEEE Transactions on Visualization and Computer Graphics, 2011. 17(3):p. 305-319.
• end •
