MP-PIC颗粒间作用力：固相应力
上期文章简要介绍了MP-PIC模型
为进一步解决堆积颗粒的极限体积分数问题，本期文章将主要介绍MP-PIC颗粒间作用力：固相应力。
在MP-PIC方法的早期文献中，仅有法向应力模型，且将其表示为非对角线元素为零的固相应力的梯度，后续出现的MP-PIC文献中，增加了粘性应力模型来修正之前颗粒束碰撞后颗粒无散射的问题。
图1演示了两股稀相颗粒射流在（a）仅有固相应力，（b）增加阻尼模型以及（c）增加阻尼模型和各向同性模型后颗粒散射效果。
图1 两股稀相颗粒射流不同模型下的散射效果[2]
对于浓相颗粒流，在这两类固相应力模型中，法向应力模型对计算过程的稳定性起到了决定性作用[1]，这是因为当固相处于浓相并接近堆积浓度时，如果法向应力模型不能阻止颗粒过度堆积，颗粒浓度就有可能进一步增加而超出理论上限从而导致模拟很快变得不稳定甚至发散。
目前，基于MP-PIC方法的软件或求解器主要有Barracuda、Fluent DDPM、MFIX-PIC和OpenFOAM中的MPPICFoam。其中Barracuda、Fluent DDPM为商业软件，MFIX-PIC和MPPICFoam属于开源软件，用户可进行二次开发加入自定义的应力模型。Fluent DDPM，MFIX-PIC和MPPICFoam采用了各具特色的数值算法和应力模型，但这些软件或求解器运行于不同软件平台。
鉴于法向应力模型对MP-PIC模拟的重要性，本篇将对比几种不同颗粒法向应力模型。为避免混淆，在之后内容中将法相应力模型称为固相应力模型。
模型1
Harris & Crighton模型
在大多数MP-PIC文献中，固相应力的表达式采用的是基于Harris & Crighton模型[3]的扩展形式，其基于流体网格的连续表达式为：
方程1
式中Ps 为压力常数，ζ为10-7量级的小数，γ为2到5之间的常数。分母中第二项的加入是为了消除当Ɛp 接近堆积浓度时分母为零的奇异问题，同时其物理意义是因为当颗粒在接近极限浓度时会发生转移和重排轻微超出物理极限。
可以看到，由于方程的高度非线性，只有当颗粒浓度接近堆积极限时此作用力才会起主导作用，而颗粒处于稀相时，颗粒应力对颗粒的作用很小甚至可以忽略。
模型2
Lun模型
在Fluent DDPM模型[3]中，对颗粒动量方程采用了特殊处理来防止颗粒体积分数超出用户给定的极限值以避免出现非物理值，从而允许用户模拟各类高浓度颗粒流。其采用的颗粒应力模型与TFM中的摩擦应力模型一致，由颗粒动理论给出。Harris 和 Crighton的颗粒应力模型仅与颗粒浓度相关而没有考虑颗粒速度的波动，而Lun[4]提出的基于浓相气体动理论的连续固相应力模型则考虑了颗粒温度影响，
方程2
其中：e为恢复系数；
是平均颗粒密度。
颗粒温度的计算式：
方程3
其中C为瞬时速度减去基于网格计算的颗粒平均速度。
径向分布函数：
方程4
模型3
指数形式固相应力模型
基于相压力模型的指数形式固相应力（Exponential）模型源自常颗粒粘度模型[5]（Constant Particle Viscosity - CPV）。这里只用其中的颗粒压力模型，其颗粒间压力的计算基于弹性概念，表示为
方程5
其中c1 和c2 均为经验常量，Ɛg,min为能够达到的最小气相体积分数。相比于之前的两种颗粒应力模型，指数模型在颗粒堆积浓度处没有间断的跃升，而是平滑过渡。但相比于其他应力模型，由于在堆积浓度计算的应力值较小，不能有效阻止颗粒过度堆积，而是在超出堆积浓度一定量后才急剧增加。
Harris & Crighton模型 VS Lun模型 VS 指数形式（Exponential）模型
假定：
颗粒堆积浓度为0.62；
Harris 和 Crighton模型中Ps为1Pa，β为3；
Lun模型中为0.9；
指数模型中c1 为1000，c2 为500。
计算的相应的颗粒应力值如图2所示：
图2 不同颗粒应力模型计算的颗粒应力值
从图中可知，Harris & Crighton模型以及Lun模型在堆积浓度处产生了间断，而指数模型则连续变化，且得到的颗粒应力值很小。
模型4
Coarse Grain 模型
由于应力模型的唯一作用是为了防止颗粒过度堆积，因此基于离散元模型直接计算颗粒间作用力的方式也同样可用于MP-PIC模拟。但由于需要具体计算颗粒对间作用力，而不是用模型的方式，其计算量将是所有模型中最大的。
Patankar[6, 7]在MP-PIC中引入颗粒-颗粒和颗粒-壁面碰撞力来确保颗粒-颗粒以及颗粒-壁面间没有任何的重叠区从而避免在所有计算域中颗粒体积分数超出堆积极限。颗粒-颗粒碰撞和颗粒-壁面作用使用DEM模型计算：
方程6
方程7
表示由于与q粒子碰撞而作用于p颗粒上的作用力kc，ηc为颗粒刚度，Rep为阻尼系数，表示数值颗粒p的有效半径。
同样，颗粒-壁面作用力可以通过下式计算：
方程8
方程9
其中dpw为数值粒子中心到壁面的距离。
因此最终数值粒子p的受力为所有碰撞作用力的和为：方程10
其中N表示所有与p有接触的数值粒子。
模型5
着色函数
从图2可以看到，应力模型计算的固相应力在颗粒浓度较低时基本可以忽略，而在接近或超过堆积浓度时又极其的大，通过其计算的修正速度将会违背物理事实（反射速度远大于碰撞前颗粒速度）。
由于颗粒应力模型的唯一作用是阻止颗粒进入到已经达到堆积浓度的流体计算网格内。因此，是否精确计算颗粒应力对颗粒速度的计算是无关紧要的，于是通过使用着色函数[8]标记出那些超过颗粒堆积浓度的流体单元，也同样能实现颗粒应力模型的功能。这样，颗粒的修正速度不再需要通过应力模型计算，而仅需要通过应力模型得到修正速度的方向即可，颗粒速度则可以碰撞前的颗粒和堆积颗粒的平均速度获得。
着色函数定义为
方程11
通过此函数代替颗粒法相应力，可以计算修正速度的方向
方程12
即修正速度的方向为颗粒浓度梯度的反方向。由于直接使用浓度梯度作为修正速度，因此，将不能显式或隐式速度修正方法，而需要进行特殊的速度分支判断[8]。
Barracuda软件使用Harris & Crighton模型，Fluent DDPM使用Lun模型，MFIX-PIC使用着色函数，基于OpenFOAM的开源和容易扩展的特性，OpenFOAM的MPPICFoam可实现不同的应力模型，方便对模型进行对比。
图3为在MPPICFoam不同应力模型下Goldschmidt流化床的瞬时颗粒运动模拟结果对比。
图3 不同应力模型下颗粒运动模拟结果对比
从图3可以看到：
前面4种应力模型都能预测出清晰的气泡，后2种应力模型的模拟结果中颗粒则比较均匀的分布于床体空间内，没有气泡产生；
Harris和Lun应力模型计算的气泡无论在大小和生成频率上都比较一致，且都为典型的“肾脏”形状气泡；
Exponential应力模型产生的气泡会相对中心轴线有一定偏移，且由于当颗粒浓度接近于堆积浓度时，计算的应力较小，在一些区域出现了固相浓度大于堆积浓度的现象；
Color的应力模型虽然也能产生气泡，但是同一时刻存在着多个气泡，而且可以看到，在壁面附近也出现了气泡。
关于MI-PIC模型的数值实现将在下一篇文章MP-PIC固相运动数值算法中进一步阐述。
参考文献
[1] Enwald, H. and A.E. Almstedt, Fluid dynamics of a pressurized fluidized bed: comparison between numerical solutions from two-fluid models and experimental results. Chemical Engineering Science, 1999. 54(3): p. 329-342.
[2] O'Rourke, P.J. and D.M. Snider, Inclusion of collisional return-to-isotropy in the MP-PIC method. Chemical Engineering Science, 2012. 80: p. 39-54.
[3] ANSYS Fluent Theory Guide. Release 15.0. 2013: ANSYS, Inc.
[4] Lun, C.K.K., S.B. Savage, and D.J. Jeffrey, Kinetic theories for granular flow: inelastic particles in Couette flow and slightly inelastic particles in a general flowfield. Journal of fluid mechanics, 1984. 140: p. 223-256.
[5] Johansson, K., B.G.M. van Wachem, and A.E. Almstedt, Experimental validation of CFD models for fluidized beds: Influence of particle stress models, gas phase compressibility and air inflow models. Chemical Engineering Science, 2006. 61(5): p. 1705-1717.
[6] N.A.Patankar, D.D.J., modeling and numerical simulation of particulate flows by the Eulerian-Langrangian approach. International Journal of Multiphase Flow, 2001. 27: p. 1659-1684.
[7] Patankar, N.A. and D.D. Joseph, Lagrangian numerical simulation of particulate flows. International Journal of Multiphase Flow, 2001. 27(10): p. 1685-1706.
[8] 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.
