欧拉型颗粒模型
欧拉型颗粒的理论基础
概述
老师，欧拉型颗粒模型是什么？是把粒子当作连续体来处理吗？
没错。欧拉型颗粒模型（Eulerian Granular Model）是将粉末或颗粒集合体作为“颗粒相”这种拟连续体来处理，并与气相一起在欧拉-欧拉法的框架下求解的方法。广泛应用于流化床、气力输送、旋风分离器、粉末混合等气固两相流问题。
和DEM-CFD有什么区别呢？
DEM-CFD是离散地追踪每个粒子，但当粒子数量达到数百万至数十亿的工业规模时，计算成本会变得非常庞大。欧拉型颗粒模型将粒子群作为连续体处理，因此即使是大型系统也能在现实的计算时间内求解。不过，它无法直接处理单个粒子的接触力或形状效应。
控制方程
颗粒相的方程是怎样的？
为了描述颗粒相的运动，使用KTGF（颗粒流动力学理论）。固相的连续性方程与通常的欧拉-欧拉法相同，但其特点是固相应力张量采用由KTGF推导出的本构关系。
KTGF是气体分子运动论的粒子版本吗？
正是如此。用颗粒温度 $\Theta_s$ 来表征粒子的速度脉动。
颗粒温度的输运方程如下。
右边各项依次是：剪切产生项、扩散项（$\kappa_s$是扩散系数）、非弹性碰撞耗散项（$\gamma_s$）、气-固相间能量交换项（$\phi_{gs}$）。
固相的本构关系
固相压力和粘度是怎么确定的呢？
固相压力 $p_s$ 由颗粒温度和体积分数求得。Lun等人（1984）的模型具有代表性。
这里 $e_{ss}$ 是粒子间的恢复系数，$g_0$ 是径向分布函数。当体积分数接近最密填充率 $\alpha_{s,max}$（约0.63）时，$g_0$ 会急剧增加，这描述了粒子密集状态下的接触压力。
| 本构关系 | 模型示例 | 物理量 |
|---|---|---|
| 固相压力 | Lun, Syamlal-O'Brien | $p_s(\alpha_s, \Theta_s)$ |
| 固相粘度 | Gidaspow, Syamlal | $\mu_s(\alpha_s, \Theta_s)$ |
| 固相体积粘度 | Lun et al. | $\lambda_s$ |
| 摩擦应力 | Schaeffer, Johnson-Jackson | 密填充区域的应力 |
沙漏的物理——颗粒既非“流体”也非“固体”
为什么沙漏中的沙子看起来在“流动”，却能在沙堆的斜坡上静止？欧拉型颗粒模型为回答这个问题提供了关键。颗粒是需要流体力学连续体方程和固体力学弹性压力两者的“第三状态”。颗粒温度这一概念，是模仿分子的热运动，将粒子速度的脉动视为温度的独创性想法，由Jenkins & Richman在1980年代从气体分子运动论推导出来。如果没有这个理论，甚至连预测CFB锅炉的粒子循环速度在一个数量级内都将是困难的。
数值解法与实现
数值解法的细节
请告诉我欧拉型颗粒模型在数值上的要点。
最大的难点是固相体积分数接近最密填充率 $\alpha_{s,max}$ 时的处理。在这个区域，固相压力会急剧增大，数值上容易变得不稳定。
摩擦应力模型很重要，当 $\alpha_s > \alpha_{s,min}$（通常为0.5）时，使用Schaffer模型或Johnson & Jackson模型来添加摩擦压力和摩擦粘度。
曳力模型的选择
气固两相流的曳力模型该怎么选？
Gidaspow模型是最通用的。它在体积分数0.8处切换Ergun方程（密填充区域）和Wen-Yu方程（稀疏区域）。
| 区域 | $\alpha_g$ | 模型 | 公式 | ||
|---|---|---|---|---|---|
| 密填充 | $< 0.8$ | Ergun | $\beta = 150 \frac{\alpha_s^2 \mu_g}{\alpha_g d_s^2} + 1.75 \frac{\alpha_s \rho_g \ | \mathbf{u}_g - \mathbf{u}_s\ | }{d_s}$ |
| 稀疏 | $\geq 0.8$ | Wen-Yu | $\beta = \frac{3}{4} C_D \frac{\alpha_s \alpha_g \rho_g \ | \mathbf{u}_g - \mathbf{u}_s\ | }{d_s} \alpha_g^{-2.65}$ |
切换点不连续不会出问题吗？
确实如此，切换点的不连续性有时会引起空隙率的振荡。Huilin-Gidaspow模型引入了平滑的过渡函数进行了改进。Syamlal-O'Brien模型也使用全域连续的公式，因此稳定性更高。
OpenFOAM中的实现
在OpenFOAM中使用哪个求解器？
multiphaseEulerFoam 支持欧拉颗粒模型。KTGF的各个本构关系可以在 kineticTheoryModel 类中选择。主要设置在 constant/phaseProperties 中进行。
Fluent中的设置
在Ansys Fluent中，在欧拉多相流模型内启用颗粒相。重要的设置项目如下。
| 设置 | 推荐 | 备注 |
|---|---|---|
| Granular Viscosity | Gidaspow | 标准 |
| Granular Bulk Viscosity | Lun et al. | 体积粘度 |
| Frictional Viscosity | Schaeffer | 密填充区域 |
| Packing Limit | 0.63 | 随机填充率 |
| Restitution Coefficient | 0.9 | 玻璃珠典型值 |
KTGF收敛之壁——与“压力变为负值”现象的斗争
实现欧拉型颗粒模型时，首先遇到的难题就是固相压力变为负值导致发散的问题。这是在颗粒温度急剧收敛到0时发生的数值不稳定现象，设置底限值（最小值约1e-10 m²/s²）实际上已成为必须的对策。ANSYS的文档中给出了“作为偏微分方程求解还是作为代数近似求解”的选择，但在高密度填充（α_s > 0.4）时，不用PDE法往往精度不够，需要在现场判断计算成本与精度的权衡。