基于DEM-CFD的粮堆气流阻力与传热特性模拟研究

陈桂香，刘文磊，刘超赛，郑德乾，岳龙飞

(河南工业大学土木工程学院，郑州 450001)

固定床是一种充满随机填充颗粒的固定容器，广泛应用于农业干燥、化学反应以及核反应堆之中[1-3]。粮食籽粒作为常见的颗粒物质，在粮堆固定床中存在随机堆积过程，使床层内的孔隙分布十分复杂，对床层内气体流动、温度传递过程产生影响，因此研究固定床内部流体的流动阻力与传热特性，对提高粮食干燥效率与固定床的设计具有指导作用[4, 5]。

国内外学者对粮堆气流阻力与传热特性进行大量的实验与模拟研究。陈竹筠[6]利用设计的实验台，研究了不同流速、孔隙率以及床层高度对玉米堆通风阻力的影响，并对Ergun方程进行了修正。朱英开等[7]采用均匀的多孔介质假定对粮层阻力进行模拟，得到了3种粮食的压降分布规律。Kashaninejad等[8]对散装开心果进行气流阻力测定，研究了填充方式等对结果的影响，对比3种气流阻力模型，Ergun方程具有最佳的拟合效果。Zare等[9]设计了深床干燥实验装置，通过改变入风时气流温度与流量，测定了粮食温度变化规律，给出了稻谷传热效率的计算方式。Ranjbaran等[10]采用计算流体力学(CFD)方法对深床干燥器中的稻谷进行数值模拟，利用编写的自定义函数，对深床稻谷传热过程及传热性能进行了研究。

多数学者采用宏观水平的多孔介质假定模型对粮堆阻力进行模拟研究，没有考虑粮食颗粒在随机堆积过程中孔隙分布不均匀对气流特性的影响。黄凯[11]利用EDEM软件构建玉米颗粒进行随机堆积，根据孔隙结构参数信息，建立了孔道网络通风模型。Rusinek等[12]采用离散单元法(DEM)预测了油菜籽在筒仓中的温度分布，描述了油菜籽的自发热过程。刘立意等[13]采用CFD-DEM耦合方法，对2个品种稻谷的通风阻力进行数值模拟分析，得到了粮层阻力与粮层深度、表观风速的关系，但没有直接将孔隙结构导入到流体计算过程中。离散元(DEM)方法已广泛地应用于颗粒介质的堆积、压缩等研究中，基于DEM-CFD耦合的方法，能更加真实的模拟堆积床中气流流动与传热过程。

本研究通过离散单元法(DEM)构建所需孔隙率的球床堆积结构，将球床堆积信息导入前处理软件中并对颗粒间的接触进行处理，利用CFD模拟大豆固定床中的气体流动与传热过程，研究不同入口风速对粮堆气流阻力与传热特性的影响，为后续粮堆非均质孔隙结构中的多场耦合研究提供参考。

1 数值模拟模型及方法

1.1 离散元堆积模型

本研究选用PFC3D离散元软件模拟了颗粒的随机生成过程，与连续多孔介质模型相比，离散单元法(DEM)将颗粒材料视为不同粒子的集合，球形颗粒由于重力向下运动，每个粒子均遵循基本的物理和力学定律，考虑了粒子物理特性与相互作用特性，当颗粒达到稳定时，模拟堆积过程停止[14]。本研究设定颗粒的运动速度小于10-7m/s时达到稳定，共生成1 418个颗粒，表1为离散元模拟所需的参数值[15]。

表1 大豆颗粒堆积床模拟的离散元参数

在固定床层的实验与模拟研究中，需要研究壁面效应对结果的影响，Mueller等[16]研究表明，固定床直径D/dp≤10时，壁面效应对堆积填充的孔隙率影响较大。图1为本研究的几何模型，由三部分组成：空气入口段，颗粒堆积段、空气出口段。出口与入口段的长度为20 mm，堆积段长度为100 mm，固定床层的直径为70 mm，颗粒直径为6.5 mm(D/dp=10.7)。

图1 几何模型示意图

1.2 网格划分及独立性验证

1.2.1 网格划分

采用自开发的脚本程序，生成上述构建的随机堆积颗粒与几何模型，重建出的球形颗粒之间会形成点接触或狭小的间隙，不利于网格的生成，直接进行网格划分时，接触点附近的网格单元会有较大的偏斜率，导致模拟难以收敛，影响计算精度，因此需要对重建颗粒间的接触进行处理。现有的主要处理方法包括缩径法、扩径法、搭桥法和截面法[17]。缩径与扩径法改变了球体的半径，使整体的孔隙率发生变化，截面法与搭桥法只对接触点局部进行处理，保证了球床几何模型的完整性。本研究选用搭桥法进行接触处理，能更准确的模拟动与传热的过程，圆柱形桥的直径为0.1dp，高度为0.05dp，网格划分与接触处理如图2所示。

图2 网格划分与接触处理图

1.2.2 网格独立性验证

为保证模拟结果的准确性，进行网格独立性研究，并选择合适的网格尺寸。采用非结构四面体网格对模型进行划分，设置4种网格尺寸分别为1/5、1/10、1/15dp和1/20dp。在同一模拟设置条件下，网格独立性计算结果如图3所示，当网格尺寸为1/15dp时，固定床中定点处流体压降和温度计算结果无明显差别，满足无关性要求。因此，本次模拟选用最大网格尺寸1/15dp进行研究，网格数量为612万。

图3 网格独立性验证

1.3 CFD计算模型

本研究采用ANSYS Fluent对流动与传热进行计算求解，模型计算区域的进口和出口分别设置为速度入口边界和压力出口边界，侧壁与粮食颗粒边壁设置为固定边界。Dixon等[18]对固定床层中的流动进行了细致的划分，研究表明当颗粒雷诺(Reynolds)数Rem大于300时，床层内的流动属于湍流，在本次模拟中入口风速设置在0.5～2 m/s范围内，Rem为358～1 432，由于固定床内部通道复杂，故选用RNG-k-ε湍流模型，壁面进行标准函数处理。控制方程采用基于压力-速度方法进行稳态求解，压力-速度耦合求解器采用SIMPLE算法，压力项采用PRESTO!格式离散，其余方程采用二阶迎风格式离散，为了得到更精确的结果，所有变量的收敛残差设置为10-5，表2为模拟初始参数条件。

表2 CFD数值模拟参数

2 计算结果分析

2.1 径向孔隙率分布

为了得到径向孔隙率分布，将床层沿径向划分为多个同轴圆柱形，如图4所示。

图4 球床径向剖面图

孔隙率为圆柱面截取的流体面积与圆环面积的比值，计算方法如式(1)所示。

(1)

式中：As为径向半径r处圆柱面上固体区域的面积/m2；H为颗粒堆积区域的高度/m。

通过Fluent中的切面工具沿固定床径向设置45个同心圆柱面，利用Report工具计算出柱面上的孔隙率，离散元(DEM)堆积结构径向孔隙分布如图5所示。De Klery[19]提出了一种广泛用于预测圆柱形填充床内孔隙率分布的经验公式，如式(2)所示。

(2)

由图5可知，根据离散元(DEM)堆积结构计算出的径向孔隙率分布与De Klery[19]提出的经验公式基本一致，径向孔隙率向床层中心方向振荡，离中心距离越近振幅越小，并趋于恒定值。当无量纲距离n=0(r=R)时，径向孔隙率最大为1，这是因为最外侧圆柱面与球颗粒只存在接触点，随着无量纲距离增大，孔隙率逐渐变小，最终趋近于球床平均孔隙率。

图5 球床径向孔隙率分布

2.2 压力场分析

粮堆压力降是固定床中流动特性的重要特征参数，它受颗粒雷诺数与壁面阻力系数的影响。修正后的颗粒雷诺数计算公式为：

(3)

式中：Rem为颗粒雷诺数；ρ为流体密度/kg/m3；u为入口流速/m/s；dp为颗粒直径/m；ε为球床孔隙率；μ为流体的动力粘度/(Pa·s)。

Ergun方程是用于计算多孔介质中阻力系数与压力降最为广泛的方法[20]，Reichelt[21]考虑壁面效应影响，对Ergun方程的相关系数进行修正，图6为不同入口风速下颗粒雷诺数与阻力系数值。

图6 不同入口风速下颗粒雷诺数与阻力系数

本研究选用Ergun方程与Reichelt方程验证DEM-CFD数值模拟的正确性。具体如式(4)～式(8)所示。

(4)

(5)

(6)

(7)

Bw=[1.15(dp/D)2+0.87]2

(8)

图7为不同入口风速下球床压力降模拟结果与Ergun、Reichelt方程计算结果对比。DEM-CFD模拟结果与公式计算结果基本一致，粮层压力降与入口速度近似呈二次函数关系，拟合公式见式(9)。

(9)

在入口速度小于1 m/s时，模拟结果与两种方程计算结果较为吻合；当入口速度大于1 m/s时，模拟结果与Reichelt方程的计算结果平均误差为5.41%，由于Reichelt方程考虑了壁面效应对压力降的影响，与Ergun方程计算结果相差较大。

图7 不同入口风速下球床压力降

由图8可知，4种入口风速下压力分布趋于一致，沿出口方向压力逐渐降低，压力分布沿床层高度呈现明显分层现象。对比4组模拟结果，随着入口速度的增加，出口处的压力下降的越快，与压力降计算结果一致。在进口与出口段处压力分布恒定，床层内部由于颗粒的阻挡，使得球床内局部压力变化明显。

2.3 速度场分布

由图9可知，4种入口风速下形成的速度场分布相似，随着入口速度的增加，球床内部风速也越大。沿球床轴向方向上，入口处风速与填充段前段颗粒接触时，颗粒会对气流形成阻挡，使入口风速不均匀分布；在填充段内部，由于颗粒的随机分布流体通道十分复杂性，使得填充段内速度分布存在明显差异，局部平均风速增大到入口风速的2.5倍；在出口端，气流会出现停滞与回流现象。沿球床径向方向上，近壁面处孔隙率最大，向床层中心方向逐渐减小，可以发现孔隙度较高区域的气流速度明显高于其它区域。

2.4 温度场分布

由图10可知，随着入口风速的增加，球床内温度下降区域越大。由于近壁面处空气流速高于床层中心区域，流体与近壁面区域的颗粒表面热交换更为充分，近壁面处流体的平均温度低于中心区域，在低速条件下，壁面处与中心区域的温度差更为明显。球床内部孔道弯曲复杂，流体温度下降幅度较小，在出口段处，由于流体速度的停滞与回流，出口段流体温度明显高于堆积内部温度。

3 结论

本研究通过离散元法(DEM)在固定床层内生成大豆颗粒的随机堆积结构，并使用CFD方法直接模拟所创建的颗粒与流体区域。研究了不同入口速度下床层内的流动与传热特性，包括床层孔隙率分布、速度场分布、压力降以及温度分布等规律。结果表明：

图8 不同入口风速下球床压力分布云图

图9 不同入口风速下X=0截面处速度分布云图

图10 不同入口风速下X=0截面处温度分布云图

固定床近壁面处孔隙率最大并趋近于1，径向孔隙率向床层中心方向呈振荡趋势，离床层中心距离越近，孔隙率与振幅逐渐减小，最终趋近于球床平均孔隙率。随着入口风速增加，压力降增加，粮堆压力降与入口速度近似呈二次函数关系，考虑了壁面效应的Reichelt方程较Ergun方程对压力降的计算更准确。固定床中气流速度与温度场分布受孔隙率分布、入口风速影响较大，床层内部速度与温度场分布不均匀。近壁面区域气流速度是床层中心速度的2.5倍，相比于中心区域，近壁面处流体的平均温度降低更快，在出口处会出现气流的停滞区。