和 
分别是颗粒体积和直径，
是颗粒滑动摩擦系数，
和 
分别是法向和切向重叠量，
是颗粒的杨氏模量，
是泊松比。
和 
分别是流体粘度和密度，
是流体相的体积分数，
是流体速度。
基于欧拉模型对流体相进行模拟。其质量和动量守恒的控制方程为： 其中， 
和 
分别是流体密度、平均速度和脉动速度、压力、流体剪切应力张量和重力加速度。 
是流体相的体积分数。在模拟中，雷诺应力张量项通过k-ε模型求解。 
是在给定位置和时间上，所有颗粒对流体施加的每单位体积的合力，由下式给出： 其中， 
是流体单元的体积，N是流体单元中的总颗粒数。
2.2. 模型实现
在这项工作中，模拟是基于内部DEM代码和ANSYS-Fluent的框架进行的，其中颗粒运动通过内部DEM代码在GPU上计算，流体流动则由CPU上的ANSYS Fluent解决。CFD与DEM之间的耦合计算（如确定颗粒的 CFD 单元位置、考虑其体积分数以及将相应的源项合并到动量方程）基于DDPM方案中的函数。DDPM通过将离散相的体积分数引入质量和动量守恒方程中，解决了标准DPM的局限性，允许在稠密颗粒系统中考虑颗粒与连续相之间的相互作用。两个求解器之间的通信通过Windows Winsock编程实现。 1）CFD在CPU上的实现
CFD求解器利用MPI在CPU上进行并行化，将计算域划分为多个子域，并将每个子域分配给不同的线程。每个核心都与其他线程同时在其各自的数据集上执行相同的程序。具体来说，在 ANSYS Fluent 中，进程由Cortex、主机进程和一组n个线程的进程组成。主机进程接收来自 Cortex 的命令，通过 process-0 向所有进程发送命令，并从 process-0 收集消息并执行写入文件、数据传输和在所有数据上显示消息等操作。CFD中的计算流程如图2所示。在模拟中，颗粒位置/属性和套接字初始化首先由DEFINE_INIT的UDF（用户定义函数）在主机进程中读取并运行，并且流体流动更新在线程进程中执行。之后，根据从 DEM 获得的颗粒信息，在主机进程中准备DPM写入并将其分发到所有线程。随后，在线程进程中执行DDPM方案，包括确定颗粒的CFD单元位置，考虑其体积分数，并将相应的源项合并到动量方程。之后，存储在不同线程进程中的颗粒阻力通过node-0传输到主机进程，然后将颗粒信息发送到主机进程中的DEM。最后，在所有进程中清除DPM写入，以避免CFD计算时的中重复计算，为下一次循环迭代做好准备。在这一部分中，使用了UDF DLLEXP的两个函数，“fl_dpm_pre_solve_update”函数用于清除DPM写入，而“fl_dpm_post_solve_update”函数用于创建写入、发送/接收数据和计算阻力。需要注意的是，通过利用ANSYS Fluent中内置的DDPM（密集离散相模型）方案，模拟可以有效利用结构化和非结构化网格。与通过UDF DEFINE_SOURCE直接添加动量源项相比，该方法在CFD模拟中提供了更好的收敛性和稳定性，特别是对于涉及复杂几何形状和非结构化网格的情况。主要是由于DDPM方案中使用了基于节点的平滑算法，可以通过平滑加载将颗粒的影响均匀分布到相邻网格节点和相邻单元上，从而减少了网格依赖性和计算不稳定性。 图2. CFD-DEM耦合计算流程，方框内的步骤在GPU上执行 2）DEM在GPU上的实现
图2还展示了DEM在GPU上的计算流程。计算从数据初始化开始，从输入文件读取颗粒属性，然后进行套接字初始化并建立与 CFD 中套接字的连接。随后，根据输入数据在CPU和GPU中创建阵列的内存，并根据颗粒大小生成颗粒搜索网格。然后，DEM中使用的边界壁面被离散化并分配给颗粒搜索网格。需要注意的是，上述模拟处理在模拟期间仅执行一次。之后，一旦从 CFD 接收到颗粒属性，就会将其复 制到 GPU 内存，然后在 GPU 上执行 DEM 计算，其中包括基于颗粒搜索网格创建颗粒邻居列表，计算颗粒之间以及颗粒与壁面之间的相互作用，以及更新颗粒的位置和速度。DEM 计算完成后，将 DEM 时间与 CFD 时间进行比较。如果二者同步，颗粒信息从GPU复 制到CPU并发送到CFD。否则，将进入新的DEM循环，直到与CFD时间步同步。 在GPU DEM中，一系列三角面元用于表示具有复杂几何的边界壁面以进行高效计算，每个三角面元都被分配一个唯一的编号并与DEM单元相关联。如图3所示，如果三角面元的离散点落在相应的DEM单元内，则将该三角面元添加到该单元的数组中。离散点之间的最小距离设置为DEM单元大小的1/10和最大面长度的最小值，以确保每个三角面元都被正确分配给DEM单元。需要注意的是，该算法仅在主模拟循环之前执行一次，并将排序后的数组存储在GPU内存中以便后续使用。三角面元分配的伪代码可以在补充材料的表S1中看到。 图3. DEM 单元中的边界壁面和壁面的离散化，用于识别所占用的壁面的 DEM 单元。 颗粒与三角面元所组成的边界壁面之间力的计算遵循以下步骤。首先，根据颗粒的位置确定颗粒所在DEM单元的索引。然后，在 26 个相邻 DEM 单元以及该单元内包含的所有三角面元上执行循环。在循环过程中，检查颗粒与三角面元之间的距离，并确定颗粒与壁面之间的接触类型。最后，根据接触类型计算颗粒-壁面之间的作用力。关于颗粒与壁面之间接触检测的详细信息可以在Su等的工作中看到。颗粒-壁面作用力计算的伪代码可以在补充材料的表S2中看到。
3.1.流动动力学验证
根据实验室规模二维填充床实验测量中的空腔形态进行验证。填充床的尺寸为长260毫米、高600毫米、宽17.7毫米。用于向填充床供应气体的风口直径和深度分别为10毫米和30毫米。床中的颗粒直径为3-4毫米，密度为1,000 kg/m³。填充床的高度为350毫米。实验中注入床中的气体密度为1.23 kg/m³，粘度为1.78 × 10^(-5) kg/(m s)，流量为16 m³/h。在模拟中，使用相似原理对几何尺寸和颗粒直径进行了缩放。缩放后，几何长、高、宽分别为2,600毫米、3,500毫米和88.6毫米。模拟中注入床中的气体密度为0.716 kg/m³，速度为232 m/s，颗粒直径为40毫米，填充床中的颗粒总数为11,000。颗粒的杨氏模量为2 × 10^7 Pa，颗粒之间的滚动摩擦系数为0.01。恢复系数设置为0.6，颗粒之间的摩擦系数设置为0.63，而颗粒与壁面之间的摩擦系数设置为0.56。颗粒和气体的其他属性与实验保持一致。 DEM网格尺寸设置为颗粒直径的2.5倍，相邻颗粒搜索间隙为颗粒直径的0.8倍，并在后续研究中保持一致。CFD网格尺寸为44.3毫米，CFD 生成的网格数为 7,640。CFD的时间步长为2 × 10-4秒，DEM的时间步长为1 × 10-5秒。使用4个CPU核心进行CFD和1个GPU进行DEM的模拟需要20分钟来计算1秒的物理时间。需要注意的是，在以下工作中，模拟在搭载Intel Core i9-13900 K CPU和NVIDIA GeForce RTX 3090 Ti GPU的计算机上执行。所有模拟都在两个求解器中都使用了双精度浮点运算。 图4显示了模拟和实验测量之间的空腔形态的对比。结果表明，实验和模拟中的空腔形状具有可比性。具体来说，在实验中，空腔的深度为105.1毫米，高度为120.0毫米，高度与深度的比值为1.14。相比之下，在模拟中，空腔的深度为0.97米，高度为1.196米，高度与深度的比值为1.23。实验与模拟之间的差异约为8.1%，表明模型在空腔形态方面得到了验证。
(a) 实验（Lu等，2020），(b) 模拟和(c) 模拟中心面的体积分数 3.2. GPU CFD-DEM效率验证
为了验证使用基于GPU的DEM-CFD模拟的加速效果，将模拟结果与之前使用CPU进行的模拟进行了比较。该模拟基于实验室规模的鼓泡流化床实验，模拟条件的详细信息如表2和图5所示。模拟和实验都使用了总质量为1.9千克的尼龙珠。 图6显示了在气体速度为2.19 m/s时颗粒的平均水平/垂直速度。测量区域高度为0.0672 m，横向均匀分布5个监测点。每个区域的高度、长度和宽度分别为0.0457 m、0.0457 m和0.003 m。平均速度采用物理时间5s后的平均值，以减少启动效应的影响。图6 a-b显示，在相应区域，无论是水平速度还是垂直速度，都与实验数据相当。此外，模拟中床层在0.0413 m至0.3461 m高度之间的压降为0.70 kPa，这也与实验结果的0.69 kPa相当。
图6. 模拟结果与实验数据（Gopalan等，2016）的比较 气体速度为2.19 m/s时的(a)时间平均颗粒水平速度和(b)时间平均颗粒垂直速度 为了验证基于GPU的CFD-DEM模拟的效率，与之前使用基于 CPU 并行化的粗粒度方法的模拟进行了比较进行了比较。该模拟是在NCI的8核超级计算机Gadi上使用MFIX-DEM软件进行的。图7展示了不同方法在物理时间为10秒时的计算时间。结果表明，基于GPU的模拟在计算时间上大约比基于CPU的模拟低8.5倍。此外，尽管粗粒度方法使用了11,618个颗粒来模拟包含92,948个颗粒的实验室规模的鼓泡流化床，但基于GPU的模拟在计算时间上甚至超过了粗粒度比为2的粗粒度方法。需要注意的是，目前的模拟使用的是基于游戏的GPU，该 GPU 并未针对双精度计算进行优化。如果使用专业显卡，计算效率将大大提高。
基于CPU的CFD-DEM和粗粒度方法（CGM）的计算时间，数据来自参考文献Wang和Shen（2022）。 本节通过模拟高炉回旋区来展示模型在处理具有大量颗粒和复杂几何形状的风口的有效性和能力。分析了风口角度对回旋区动力学的影响，并考虑了风口的详细情况，预测了基于颗粒碰撞的风口磨损。 4.1. 模拟条件
图8显示了模拟中使用的几何形状和网格。整个高炉（下部）由40个风口组成。为了降低计算成本，仅模拟了高炉的一部分，代表整体的1/40，并且包含单个风口。此外，模拟排除了受气流影响较小的死区，只关注活性焦炭区域。考虑了三个具有不同倾斜角度的风口。向下倾斜5°（即-5°风口角度）的风口几何细节与工业中使用的相似。水平风口和10°风口的长度与-5°风口相同，而其他尺寸和几何特征则根据-5°风口进行调整。模拟中使用的尺寸和几何细节的详细信息如图8b所示。
(b) 几何细节（槽型高炉），(c) CFD网格和(d)模拟中使用的初始颗粒 图8c展示了模拟中使用的网格。总共使用了9,100个多面体网格，其中风口附近网格尺寸为0.05米，其他区域为0.08米。模拟中共使用了550,000个直径为20毫米、密度为850 kg/m^3的焦炭颗粒。在计算之前，所有颗粒都在高炉中随机生成，然后在重力的影响下降落，形成了高度为1,900毫米的填充床，直到系统的动能耗散至接近零，如图8d所示。最后，气体通过风口入口引入高炉，形成回旋区。高炉模拟的参数细节如表3所示。
4.2. 模型的有效性
图9展示了总物理时间为5秒的高炉回旋区CFD-DEM模拟中计算时间的百分比。模拟在 Intel Core i9-13900 K 的 8 核 CPU 和 NVIDIA GeForce RTX 3090 Ti 的 GPU 上进行。5秒模拟的总计算时间为12.3小时。值得注意的是，CFD求解器占据了大部分计算时间，占总计算时间的84%，而DEM求解器仅占16%，表明CFD求解器是CFD-DEM模拟中的主要瓶颈。此外，在DEM求解器中，创建颗粒邻居列表和计算颗粒-颗粒相互作用占据了大部分计算时间。相比之下，颗粒-壁面相互作用所花费的计算时间不到20%，显示出基于GPU的颗粒-壁面接触检测算法的高效性。 图9. 不同求解器和关键DEM子程序的计算时间百分比。 4.3. 回旋区动力学
图10显示了不同风口角度的回旋区动力学，其中空隙率等值面为0.55，已证明空隙率为0.55可以准确捕捉回旋区空腔。结果表明，不同风口角度下回旋区的形成是相似的。最初，回旋区沿着垂直于风口的方向扩展，随后逐渐深入高炉内部，导致回旋区高度降低。最终，焦炭颗粒和气体之间实现动态平衡，导致所有风口角度的回旋区在 0.5 秒后保持稳定。有趣的是，回旋区表现出与风口尖端相对应的更大对齐度，并具有与风口尖端相同的倾斜角度。此外，回旋区高度随着风口角度的增加而减小。
图10. 不同风口角度(a) 水平，(b) -5°和(c)-10°下回旋区的形成。回旋区形状是使用等值面方法生成的，空隙率为 0.55。 已知回旋区的尺寸对高炉的效率和生产率有显著影响。为了评估不同风口角度下回旋区的尺寸，图11比较了不同风口角度下回旋区尺寸随时间的变化。图11a展示了不同风口角度下回旋区深度随时间的变化，回旋区深度最初迅速增加，然后在0.5秒后趋于稳定，观察到最大回旋区深度出现在-5°风口，但仅略大于-10°的风口得到的结果。图11 b-e显示了不同风口角度下高炉内回旋区高度、宽度、表面积和体积的演变，最初都表现出增加，随后随时间减少。此后，略有增加，最终在 0.5 秒后达到稳定状态。然而，不同风口角度下回旋区尺寸的变化趋势存在显著差异。具体来说，随着风口角度的增加，回旋区高度减小，而回旋区宽度在不同风口角度下趋于相似。此外，在-5°风口角度下观察到最大的表面积，但仅略大于水平风口下所得到的表面积。类似地，回旋区在-5°的风口处达到其最大体积，但回旋区体积明显大于水平风口。此外，-5°风口的回旋区尺寸在0.5秒后的稳定阶段显示出较小的波动，表明与其他两个风口角度相比，回旋区形成更为稳定。上述结果表明，高炉中-5°风口相比其他风口配置具有更好的性能。 图11. 在空气体积分数为0.55的等值面下，不同风口角度下回旋区特性的演变：(a) 深度，(b) 高度，(c) 宽度，(d) 表面积和(e)体积 图12显示了在1秒到5秒的回旋区稳定期间的平均尺寸。结果表明，无论是回旋区深度还是宽度，都在-5°风口处达到最大值，并且在-10°风口也观察到相似的值。同时，随着风口角度的增加，回旋区高度逐渐减小。此外，与图11中的结果一致，-5°风口处获得了最大的回旋区体积，而进一步增加风口角度会导致回旋区体积减小。
4.4. 流场和颗粒动力学
图13显示了不同风口角度下的中心截面和水平截面的风速。结果表明，随着风口角度的增加，高风速区域扩大。此外，在水平风口顶部存在一个停滞区域，定义为低速区域（在模拟中小于1 m/s）。该停滞区域随着风口角度的增加而减小，并且在风口下方的区域也观察到类似的趋势，这表明较大的风口角度增强了高炉内的流动性。
图13. （a）水平风口、（b）-5°风口和（c）-10°风口的中心截面（上）和水平截面（下）的风速分布。 图14显示了具有三个颗粒厚度的中心平面上的颗粒速度分布。结果表明，对于不同的风口角度，回旋区中心的颗粒速度较高，而回旋区边界附近的速度相对较低。该现象可能是由于颗粒频繁碰撞而导致的在向前移动时速度下降。然而，随着风口角度的增加，高速颗粒可以达到更低的位置。图14还突出显示了风口附近的颗粒速度矢量，表明回旋区有两个循环——一个在风口的顶部，另一个在风口的底部—导致颗粒不断被送入回旋区。
图14. （a）水平风口、（b）-5°风口和（c）-10°风口的中心截面的颗粒速度分布（上）和风口附近的颗粒局部速度矢量（下）。
4.5. 颗粒活跃区
如图14所示，尽管入口速度为173 m/s，但大多数颗粒几乎是静止的。为了定义颗粒的活跃区，基于颗粒参与回旋区所需的最小速度，设定了一个0.05 m/s的颗粒速度阈值（图14）。图15a展示了在物理时间从1秒到5秒范围内，不同风口角度下的颗粒活跃区。结果表明，颗粒活跃区是一个近球形的形状，其高度比回旋区更高，对于不同的风口角度，都近似于一个椭球体。此外，风口尖端附近的活跃区对于不同的风口角度都倾向于与风口尖端垂直，与回旋区相似。该结果表明风口角度与BF内颗粒活跃区的形状/位置之间存在相关性。随着风口角度的增加，颗粒活跃区的高度降低，说明由于风口倾斜，入口空气提供的能量被消耗在风口尖端底部的颗粒碰撞中。图15b展示了颗粒活跃区的体积和表面积。与回旋区类似，颗粒活跃区在-5°风口处具有最大的体积和表面积，而水平风口的颗粒活跃区体积和表面积最小，这表明风口倾斜增强了BF中的颗粒运动。需要注意的是，如图S1和图S2所示的补充材料中，无论颗粒速度阈值如何变化，颗粒活跃区随风口角度变化的趋势都保持一致。尽管由于阈值定义的不同，形状、体积和表面积会有所变化，但观察到的趋势与图15中所示的趋势相似。 图15. 在5.0秒时，不同风口角度下(a)颗粒速度为0.05 m/s的等值面，(b) 颗粒活跃区的平均体积和表面积。 4.6. 磨损预测
已知风口的磨损是高炉风口失效的主要原因之一。为了预测风口的磨损状况，这里采用了广泛使用的Archard磨损模型来测量风口的磨粒磨损，其表达式为： 其中，dw 和A是几何形状的磨损深度和单元面积，Fn 是法向接触力，ΔU 是接触时颗粒滑动距离，K 是磨损系数，Hv 是维氏硬度。通常，风口由铜制成，铜的磨损系数通常为 6.0 × 10^-4，维氏硬度为 3.5 × 10^8 Pa。 图16a显示了风口上的局部磨损深度；可以看出，对于不同的风口角度，磨损发生在风口末端（风口尖端）附近，但具体的磨损模式因风口角度而异。在水平风口的情况下，最严重的磨损位于风口尖端的大部分区域，而风口尖端底部的磨损相对较低。对于-5°风口，严重的磨损出现在风口尖端的底部和顶部的小部分区域。模拟结果与行业报告中风口失效模式的描述相类似，该报告指出风口尖端顶部和底部的损坏最为严重，如图16b和图16c所示。类似的趋势也在-10°的风口中发现。结果表明，随着风口角度的增加，磨损减少。水平风口在风口尖端顶部和底部区域的显著磨损归因于这些区域存在循环区和风口尖端与颗粒之间的碰撞角度，如图14所示。风口尖端顶部和底部区域较高的碰撞角度导致这些区域磨损严重。通过采用倾斜风口，可以减少颗粒与风口壁之间的碰撞角度，从而随着风口角度的增加而减少磨损。
图16. (a) 5秒时不同风口角度下风口的局部磨损深度预测，(b) 工业中风口尖端的局部缺陷（Rittenschober et al., 2015），其中颜色越深表示风口损坏的频率越高，(c) 工业中风口尖端不同位置的损坏百分比（Rittenschober et al., 2015），以及(d) 风口末端随时间变化的平均磨损深度，线条是对应数据的线性拟合。平均磨损深度是风口尖端区域的平均值。 图16d显示了平均磨损深度的定量分析。它显示了在0.5秒之前磨损深度迅速增加，之后对于不同的风口角度，磨损深度随时间线性增加。水平风口、-5°风口和-10°风口的斜率分别为3.88 × 10-9、3.77 × 10-9和3.62 × 10-9，所有风口角度的截距均为2.4 × 10-9。随着风口角度的增加，斜率减小，表明在相同入口速度下，随着风口角度的增加，磨损减少。值得注意的是，尽管在5秒时磨损深度看似最小，但考虑到风口的寿命，可能会变得相对显著。通常，风口在高炉中的平均寿命为281天，风口尖端处水平、-5°和-10°风口的平均磨损深度分别为11.4毫米、11.0毫米和10.6毫米。此外，通过计算在5秒内的最大磨损并假设这种最大磨损随着时间线性增加，水平、-5°和-10°风口在281天后的最大磨损深度分别高于68.4毫米、66.0毫米和63.5毫米，表明风口损坏主要是由颗粒对风口的磨损造成的。 在本研究中，开发了一种基于CPU-GPU跨平台的CFD-DEM耦合方法，用于模拟具有复杂几何形状的大规模气固流动。颗粒运动由基于GPU的DEM求解器模拟，而流场则在CPU上使用ANSYS Fluent进行计算；CFD和DEM之间的耦合通过Winsock数据传输实现。然后应用该模型在颗粒尺度上模拟高炉回旋区动力学。可以得出以下结论： （1）该模型在空腔形态方面得到了验证。该模型的有效性在小型流化床中得到了验证，并且计算时间与CGM方法使用的先前模拟进行了比较，显示出与实验结果良好的一致性以及与CGM方法相比优越的计算速度。 （2）CFD计算时间占CFD-DEM模拟的84%，而DEM计算仅占16%。在GPU上的DEM模拟中，大部分时间用于计算颗粒-颗粒力和创建邻居列表。颗粒-壁面力的计算仅占DEM计算时间的16%。 （3）在回旋区形成过程中，回旋区的形状与风口尖端对齐，并且具有与风口尖端相同的倾斜角度。此外，当风口向下倾斜5°时，回旋区显示出其最大体积和颗粒活跃区。此外，使用倾斜风口的高炉中，气体和颗粒的流动性都得到了增强。 （4）对于水平风口，磨损累积在顶部风口尖端；当风口倾斜时，磨损移动到底部。值得注意的是，最严重的磨损发生在倾斜风口的风口尖端底部，与工业风口中观察到的故障位置一致。随着风口角度的增加，磨损减少，表明较大的风口角度可以更好地防止颗粒碰撞。 总体而言，开发的CFD-DEM模型证明了其在模拟复杂几何形状的稠密气固流动的有效性。该模型在研究高炉回旋区动力学方面的成功应用凸显了其在以大量颗粒和复杂几何形状为特征的其他工业过程中的应用潜力。 翻译转载自：《Chemical Engineering Science》“GPU-accelerated CFD-DEM modeling of gas-solid flow with complex geometry and an application to raceway dynamics in industry-scale blast furnaces” 
和 
分别是颗粒i的平移速度、角速度、惯性矩和质量。 
是从颗粒i的中心指向接触点的半径向量。 
是滚动摩擦系数。 
是颗粒i的单位角速度。 
和 
分别是颗粒i与其邻近颗粒j之间的法向接触力和切向接触力。这些方程的详细信息见表1（Zhu等，2022）。其中： 
