气固流化床启动阶段挡板内构件受力特性的CFD-DEM模拟

作者：李铁男赵碧丹赵鹏张永民王军武来源：《化工学报》日期：2022-09-05人气：1127

流化床反应器具有良好的传质传热性能，且能够较好输送固体颗粒物料，因此广泛用于工业生产过程[1]。反应器内设置合适的内构件可强化传质与传热并提高反应器性能[2-3]。基于以往研究经验，采用倾斜叶片的流化床内构件破碎气泡和改善气固接触的效果更好，而且还可以通过调整结构参数有效调节气固相的停留时间分布，具有比垂直构件和具有竖直叶片的网状格栅更好的反应强化效果，因此此类内构件在工业中应用最广泛[3]，典型的例子如百叶窗格栅、脊型和塔型内构件等。由于工业流化床内存在复杂的动态两相流动以及不同操作状态之间的切换，内构件会受到不同形式和强度的作用力，如果内构件结构设计不合理，很可能会致使内构件发生损伤甚至被破坏，进而导致工业装置出现故障和经济损失。因此，为了保障流化床内构件的长周期可靠性，系统掌握内构件在流化床内的受力特性尤为重要[4]。

早期针对流化床中内构件受力特性的研究主要采用实验手段，大多数关注的是正常流化状态下B 类颗粒流化床内圆形水平换热管的受力特性[5-9]，其中Grace等[7-8] 的研究为揭示正常流化状态下圆管受力的机理奠定了良好基础。近年来，Zhang等[4, 10-13]系统研究了化工领域A 类颗粒流化床反应器中经常使用的斜片挡板内构件的受力特性，发现启动阶段流化床层内部构件会受到一个很大的向上载荷脉冲，其峰值可达正常流化状态下内构件受到的平均载荷值的数倍甚至一个数量级以上[11]，这对工业装置内构件可靠性是一个巨大的潜在威胁，很可能是导致工业装置内构件损坏的一种新机理。

对于冷态实验室条件下一些小型简单结构内构件，实验方法可以获得较为准确的受力特性数据，但对于工业装置高温高压的苛刻条件下的大型复杂结构内构件，其在流化床中的受力特性预测则难以使用实验方法，因此，使用近年来日益成熟的计算流体力学模拟方法研究流化床中内构件的受力是一种有益的尝试。Higashida等[14]采用虚拟颗粒方法（FPM方法）模拟研究了三维鼓泡床中悬浮球体受到的动态垂直力，并将模拟计算结果与带有拉格朗日传感器系统的悬浮球体所测得的实验值相比较。研究发现，球体在装置内床层流化过程中，其动态垂直力的波动与气泡的运动密切相关，球体受到的合力由流体作用力与颗粒碰撞作用力两部分组成。Yan等[15]采用双流体模型模拟了气固两相体系内不同形状的大颗粒受力演化规律。模拟结果显示：随着表观气速的增加，大颗粒受到的作用力也相应增大。Nagahashi等[16]采用FPM模拟方法，对二维流化床中单个气泡经过水平圆管的受力机理进行了模拟研究，并将模拟结果与实验结果相比较。结果显示：模拟与实验中气泡的运动轨迹吻合较好，实验拍摄气泡经过圆管过程中，气泡尾流携带的颗粒撞击圆管会使圆管受到一个脉冲峰值。

本研究以启动阶段流化床层内部构件的非正常受力现象为研究对象，首先基于CFD-DEM数值模拟方法，建立流化床挡板内构件表面载荷强度的统计方法，并通过与前期流化床启动阶段斜片挡板内构件受力特性的实验数据[11]进行对照，以验证该载荷统计计算方法的正确性及合理性。进一步，以该方法研究不同参数（碰撞模型参数、表观气速、颗粒粒径）对启动阶段斜片挡板内构件受力特性的影响规律，旨在为下一步揭示流化床启动阶段内构件非正常受力现象的内在机理奠定良好基础，也为开发工业流化床内构件受力特性的预测工具奠定理论基础。

1 挡板受力统计计算方法

气固挡板流化床启动阶段的本质是密相颗粒料层由未流化的固定床状态逐渐转变为流化床状态的一个过渡过程，相比其他多相流模拟方法，CFD-DEM方法更适宜于准确描述这一过渡过程。因此，本研究将基于CFD-DEM方法，建立内构件表面载荷分布的统计计算方法，分析内构件表面受到的动态载荷信号。从统计力学角度分析，气固流化床密相床层中的内构件受力特性是气相分子与离散颗粒对内构件壁面作用特性的动态表征。CFD-DEM模拟方法可以获得颗粒速度、位置和相互作用力(颗粒-颗粒、颗粒-壁面)等微观颗粒相尺度信息，但是无法直接获得工程上关注的内构件表面载荷强度及分布等宏观物理量信息。

气体压力的物理本质为气固两相流体系内构件受力求解提供一定参考。内构件气固两相流体系内不同尺度下，从统计物理学可知固体颗粒与内构件表面作用的动态变化过程和气体分子与壁面作用的动态过程相类似。此外，宏观连续场信息通常都是借助颗粒流系统中每个离散颗粒和分子动力学中单个原子的作用力、位置和速度等微观信息而获得，故内构件表面受到的颗粒碰撞作用力特性可参考气体对壁面作用形成的压力过程及其动态作用力特性。在标准分子动力学模拟中，基于维里定理求解计算表面压力。对于非均匀系统，压力应为张量形式，且是与位置相关的函数[17]。其中，压力张量的公式包含两部分：第一部分称为动压力，代表颗粒的运动对压力的贡献量；第二部分称为位形压力，代表粒子间的相互作用。

气固两相流体系内构件受力载荷的求解计算需要了解颗粒流体系内应力张量的组成形式。连续介质场通常需要由离散粒子数据构造，这些离散数据是每个原子或颗粒的位置、速度和相互作用力。在颗粒系统中，颗粒动理学理论是在Chapman等[18]的气体动理学理论的基础上发展起来的。一般而言，封闭固相应力的方法主要有两类：一类是颗粒动理学理论(KTGF)；另一类是基于CFD-DEM方法的数值统计分析。KTGF的推导方法主要分为以下几种。(1) 考虑非弹性碰撞及颗粒体积影响。Gidaspow[19]基于Chapman等[18]建立的气体动理论，进而考虑非弹性碰撞影响和颗粒体积影响等不可忽略因素，建立起经典的颗粒动理论。(2) 考虑固相处于近平衡稳态情况。Rao等[20]进一步考虑固相处于近平衡稳态时，通过完善颗粒速度分布函数，导出更合理的固相本构关系。(3) 考虑体系内介尺度结构广泛存在情况。Zhao等[21]考虑气固两相流系统中团聚物等介尺度结构广泛存在，建立了适用性更广且更准确的多尺度动理论。而对于另一类基于CFD-DEM方法的数值计算，可通过捕捉颗粒系统中单颗粒的运动特性，选择合理的统计平均方式，获得合理的固相本构关系。颗粒系统中应力张量的求解方法主要分为三种。（1）体积平均法。Babic[22]最早提出了体积平均法，而后Zhu等[23]在此基础上进一步发展，将颗粒质量、速度和相互作用力等微观性质与体系密度、应力和偶应力等宏观特性联系起来。（2）平面法。Mehrabadi等[24]采用虚功原理应用平面法推导出应力张量表达式，并证明平面法推导出的应力张量形式与体积平均法等效。（3）粗粒化方法[25]。基于假设颗粒间为二元碰撞，每一对颗粒都只有单一的一个接触点，碰撞接触区域被一个接触点取代，且碰撞不是瞬时的，将离散数据体积平均化并用连续性描述。此方法推导应力张量表达式不需要假设颗粒是刚性的或球形的。综上所述，颗粒体系内的固相应力张量的组成形式主要包含两个部分：(1) 动应力(kinetic stress)，由颗粒体系内局部体积中颗粒质心的速度相对于局部体速度的脉动产生；(2) 接触应力(contact stress)，由体系内颗粒之间的直接接触产生。

当颗粒体系内有挡板内构件存在时，内构件边壁与颗粒间的碰撞也会对应力张量有贡献作用。为了建立本研究中挡板内构件表面载荷的统计计算公式，需要了解内构件的存在对应力张量的贡献形式。Weinhart等[26]采用Goldhirsch[25]的体积平均方法导出了边界存在所贡献的边界应力形式。研究发现，当体系内有壁面边界存在下，应力张量除了由颗粒速度波动产生的动应力和流动颗粒间碰撞产生的碰撞应力外，还有边壁与颗粒碰撞的碰撞贡献。根据牛顿第三定律可知，壁面边界对固相颗粒的碰撞作用与颗粒施加于壁面的碰撞作用数值上相等。因此，统计挡板表面受到颗粒的碰撞载荷形式与边界碰撞应力的统计形式类似。

统计计算密相床层中的挡板内构件表面受到颗粒的碰撞力，旨在将微观的颗粒壁面碰撞信息与宏观的挡板表面受力信息联系起来，并将信息在时间与空间上统计平均得到内构件表面受到颗粒碰撞应力及其分布特性。本研究的挡板内构件形式为斜片挡板，其受到颗粒的碰撞作用主要表现在与气体分布板通入气体的流向相垂直的上、下两表面所受的颗粒碰撞力。

根据柯西基本定律，挡板内构件表面上单位面积受到的颗粒碰撞力与挡板表面面积平均的碰撞应力张量关系如式（1）所示。

$T = σ \cdot v$ (1)

式中， v 为单位向量； σ 为碰撞应力张量； T 为碰撞力。

基于挡板内构件表面受到的碰撞载荷与固相应力中内构件边壁的贡献数值上相等，故需要找到一个合理的平均化方法将颗粒的微观碰撞信息转化为挡板表面受到的宏观应力信息。对体系内边壁存在的影响处理，Ries等[27]提出了镜像系统(the mirrored system)的方法。此方法的思想是通过类似将体系内边壁一侧的颗粒关于边壁镜像对称，来代替边壁存在对固相应力的影响。因此，基于已知挡板内构件的表面面积，为确定挡板表面面积平均的碰撞应力张量 σ，可参考此方法处理内构件表面。

由于本研究所采用的颗粒物料为单一粒径的颗粒(带粒径分布的也可以以同样的方法统计分析，但是此时两颗粒质心间位移连线的长度l不是常数)，为了更合理有效地统计计算挡板表面受到的碰撞应力，将挡板表面紧贴着的一层颗粒关于挡板表面镜像对称，此时与颗粒相互作用的挡板表面等效为一层粒径与床内颗粒相同且平铺面积与挡板表面面积相同的平面颗粒。因此，挡板表面受到颗粒的碰撞作用等效为单层平面颗粒与虚拟颗粒之间的碰撞作用，可认为单层平面颗粒与碰撞颗粒组成的颗粒群间存在一个法向量与挡板表面法向量相一致的“虚拟平面”，如图1所示。

图1

图1 密相床层中挡板表面与颗粒接触示意图(a)；等效为“虚拟平面”后颗粒碰撞接触图示(b)

Fig.1 Schematic diagram of contact between the baffle and particles in dense bed(a) and contact between an imaginary plane and particles in a dense bed(b)

水平放置的挡板主要以上、下两个表面为碰撞受力面，故在此以挡板上表面为例（下表面处理方法一致）。单位面积的“虚拟平面”上存在足够多数量的颗粒间碰撞，定义 T 为单位“虚拟平面”面积上碰撞颗粒对单层平面颗粒施加的碰撞力，则 T 如式（2）所示。

$T = \sum_{i = 1}^{N_{c}} F^{i}$ (2)

式中， Fi 为单位面积“虚拟平面”的上、下第i对颗粒的碰撞接触力；Nc为单位面积上的碰撞接触点数。

Mehrabadi等[24,28]通过虚功原理证明了单位“虚拟平面”面积上颗粒碰撞力的统计求和可写成统计平均的颗粒簇体积V上的求和，此采样颗粒簇由与挡板上表面接触的真实颗粒和镜像虚拟挡板颗粒组成。“虚拟平面”上面积平均应力张量 σ 与采样体积内颗粒间相互作用产生的碰撞应力张量 σ ′相等[29-30]，即

$σ = σ^{'}$ (3)

因此，由边壁表面受到的颗粒碰撞作用贡献形式可知，统计平均的颗粒簇体积中颗粒间碰撞应力张量 σ 的统计计算公式为

$σ_{α β}^{w} = \frac{1}{V} \sum l_{α} \otimes F_{β} = \frac{1}{V} \sum_{i = 1}^{N_{c}} l_{i α} F_{i β}$ (4)

式中，上角标w代表壁面;下角标α和β表示方向分量。

本研究中，挡板上、下表面的表面积均为S；颗粒粒径为dp；一对颗粒的每次碰撞接触中，两颗粒质心间位移连线的长度用l表示，即l的大小为颗粒质心到碰撞接触点间位移长度的2倍； F 为两颗粒间的碰撞接触力；Vw为挡板表面相对应的均匀统计颗粒簇体积，如图1(b)中虚线框所示，其值大小为以挡板表面积为底，颗粒粒径为高的立方体体积，即 $V_{w} = S d_{p}$ ；N为某一时刻挡板表面的总碰撞接触点数。因此，挡板表面受到颗粒间碰撞应力张量 σ 统计计算公式为

$σ_{α β}^{w} = \frac{1}{V_{w}} \sum l_{α} \otimes F_{β} = \frac{1}{V_{w}} \sum_{i = 1}^{N} l_{i α} F_{i β} = \frac{1}{S d_{p}} \sum_{i = 1}^{N} l_{i α} F_{i β}$ (5)

从式（5）可知，在数值模拟计算程序中，每一时间步求得挡板表面受到颗粒碰撞的应力张量表达式为

$σ = (\begin{matrix} σ_{11} & σ_{12} & σ_{13} \\ σ_{21} & σ_{22} & σ_{23} \\ σ_{31} & σ_{32} & σ_{33} \end{matrix}) = \frac{1}{S d_{p}} \sum_{i = 1}^{N} (\begin{matrix} l_{i 1} F_{i 1} & l_{i 1} F_{i 2} & l_{i 1} F_{i 3} \\ l_{i 2} F_{i 1} & l_{i 2} F_{i 2} & l_{i 2} F_{i 3} \\ l_{i 3} F_{i 1} & l_{i 3} F_{i 2} & l_{i 3} F_{i 3} \end{matrix})$ (6)

式中，每个分量表示挡板表面受到的均布载荷； $σ_{33}$ 为垂直挡板表面方向的碰撞均布载荷，为主要关注参数。运用此公式统计分析CFD-DEM的模拟结果，即可分析出颗粒与壁面作用的应力分布特性。

2 模拟的几何模型及物理参数设置

本研究选取Liu等[11]的三维方形冷模流化床实验装置（图2）作为模拟研究对象。装置的主体由方形流化床床体、预分布器、气体分布板、单个斜片挡板和旋风分离器等组成。方形流化床的截面面积为0.3 m×0.3 m，总高度为5 m。单个斜片由单个板条和两个固定底座组成，固定底座通过螺栓将板条的两端固定在流化床的床壁上，使得板条可稳定置于密相床层中。其中，挡板板条的长为0.3 m、宽为0.05 m、厚度为0.008 m，且板条的长度可刚好横跨床层截面，与流化床的截面边长一致。单个斜片挡板水平安置在流化床边壁的中心位置处，安装高度距流化床底部分布板0.5 m。

图2

图2 三维方形冷模流化床实验装置示意图[11]

Fig.2 Schematic diagram of the cold three-dimensional fluidized bed with a square cross-section[11]

为了便于计算，本研究将上述实验装置进行简化处理，省略了流化床顶部的旋风分离器和底部的气体分配室。方形流化床几何模型的长度和宽度与实验等同，即长和宽均为0.3 m。实际实验中流化床内床层膨胀高度最高膨胀至2 m左右，考虑到节省计算资源和提高计算效率，将流化床的床层稀相空间的高度缩减，故方形流化床的模型总高度定为2.2 m。模型中，挡板内构件的安装位置与实验相一致，安装在距床层底部分布板0.5 m的流化床壁面中心位置处，且挡板长度保持不变，仍为0.3 m。为简化网格划分并成功生成均匀六面体网格(六面体网格尺寸为0.01 m×0.01 m×0.01 m)，故将单个斜片挡板内构件的宽度调整为0.06 m，厚度定为0.01 m。而对于本研究中的无挡板自由床而言，除略去挡板流化床内斜片挡板外，其余几何结构均与挡板流化床一致。最终得到模拟采用的方形流化床几何结构如图3所示。

图3

图3 模拟采用的流化床几何结构

Fig.3 The geometry of fluidized bed for simulation

本文选用刘对平[4] 实验中的粒径为595 μm的B类非球形石英砂颗粒作为模拟体系内的固体介质，密度等相关物料基本性质与实验测量值一致。实验前期未曾测量石英砂颗粒的球形度信息，故根据测量颗粒的起始流化气速等参数信息，采用Hua等[31]提出的方法求解估算颗粒的球形度，最终确定模拟采用的颗粒球形度为0.86。流化床中床层物料的初始自由堆积高度为1 m。模拟体系内流体采用常温常压空气（25℃、1.2 kg/m3、1.8×10-5 Pa·s）。

3 数学模型及模拟设置

3.1 数学模型

本研究所采用的粗粒化CFD-DEM方法是基于Lu等[32] 提出的EMMS-DPM (energy-minimization multi-scale-discrete particle method，基于能量最小多尺度离散颗粒)方法。故本研究中的粗颗粒质量 $m_{C G P} = k^{3} m_{p}$ ，其中，k为粗粒化率（ $k = d_{C G P} / d_{p}$ ），mp为真实颗粒的质量。颗粒-颗粒碰撞作用与颗粒-壁面碰撞作用的计算使用Peng等[33]提出的方法，采用不考虑历史碰撞影响的Hooke（胡克线弹性碰撞）模型求解计算。计算颗粒-颗粒碰撞与颗粒-壁面碰撞之间的相互作用力时，尽管非球形颗粒与颗粒和壁面之间碰撞接触力的计算，可以通过组合颗粒球元构建非球形颗粒几何模型的方法实现[34-36]，但将大幅增加计算量。因此，本文在计算颗粒-颗粒以及颗粒-壁面间相互作用力时仍然将颗粒视为球形颗粒。显然，这是一种提高计算精度与保证合理计算量间采取的折中方法，研究表明这种方法可以获得较为合理的模拟结果[37-39]。气固相间曳力的计算采用Hua等[31]提出的非球形颗粒的曳力模型求解。在Ergun曳力关联式中引入球形度，并利用Ganser非球形颗粒曳力模型[40]来计算Wen和Yu曳力关联式中的单颗粒曳力系数Cd。

气固系统DEM的控制方程组分别为气相、固相的质量守恒和动量守恒方程：

$\frac{\partial}{\partial t} (ε_{g} ρ_{g}) + \nabla \cdot (ε_{g} ρ_{g} u_{g}) = 0$ (7) $\begin{array}{l} \frac{\partial}{\partial t} (ε_{g} ρ_{g} u_{g}) + \nabla \cdot (ε_{g} ρ_{g} u_{g} u_{g}) = - ε_{g} \nabla p + ε_{g} ρ_{g} g + \\ \nabla \cdot (ε τ_{g}) - \frac{\sum_{i = 1}^{N_{C G P, c e l l}} F_{d r a g, i}}{V_{c e l l}} \end{array}$ (8) $τ_{g} = μ_{g} (\nabla u_{g} + \nabla u_{g}^{T}) - \frac{2}{3} u_{g} (\nabla \cdot u_{g}) I$ (9) $m_{C G P, i} \frac{d u_{C G P, i}}{d t} = - \frac{π}{6} d_{C G P, i}^{3} \nabla p + m_{C G P, i} g + F_{d r a g, i} + \sum_{j = 1}^{N} F_{c, i j}$ (10) $I_{C G P, i} \frac{d ω_{C G P, i}}{d t} = R \times F_{i j}^{t} - μ_{t} R_{C G P, i} |F_{i j}^{n}| ω_{C G P, i}$ (11)

式中，ρg、 ug和p分别为流体密度、速度和压强；t为时间； g 为重力加速度；εg为网格空隙率，即计算网格内流体的体积分数；NCGP,cell为计算网格中粗颗粒数目；Vcell为当前计算流体网格的体积； Fdrag,i 为计算网格内颗粒i与气流之间的曳力交换； τg为牛顿流体黏度应力张量；μg为流体的剪切黏度； I 为单位张量； $u_{C G P, i}$ 为粗颗粒平动速度； $d C G P, i$ 为粗颗粒粒径；N为某一时刻与颗粒i相互作用的颗粒总数; $I_{C G P, i}$ 为粗颗粒的转动惯量； $ω_{C G P, i}$ 为粗颗粒角速度； $F_{i j}^{n}$ 、 $F_{i j}^{t}$ 分别为法向和切向碰撞接触力； R 为颗粒质心到碰撞接触点的位移矢量；μt为滚动摩擦因数； $R_{C G P, i}$ 为粗颗粒的半径。

鉴于气固流化床内颗粒之间的碰撞接触较为频繁，本研究采用弹簧-阻尼模型对粗颗粒之间的碰撞过程进行求解计算[33]，颗粒所受的碰撞接触力等于法向接触力与切向接触力之和。

$F_{c, i j} = F_{i j}^{n} + F_{i j}^{t}$ (12) $F_{i j}^{n} = - k_{n} δ_{n} n_{i j} - γ_{n} u_{i j}^{n}$ (13) $F_{i j}^{t} = \{\begin{array}{l} - k_{t} δ_{t} t_{i j} - γ_{t} u_{i j}^{t} \\ - μ_{t} |F_{i j}^{n}| t_{i j} \end{array}$ $\begin{array}{l} f o r \\ f o r \end{array}$ $\begin{array}{l} |F_{i j}^{t}| \leq μ_{f} |F_{i j}^{n}| \\ |F_{i j}^{t}| > μ_{f} |F_{i j}^{n}| \end{array}$ (14) $k_{n} = \frac{16}{15} \sqrt[]{R_{C G P}^{*}} Y^{*} {(\frac{15 m_{C G P}^{*} u_{c}^{2}}{16 \sqrt[]{R_{C G P}^{*}} Y^{*}})}^{1 / 5}$ (15) $\frac{1}{R_{C G P}^{*}} = \frac{1}{R_{C G P, i}} + \frac{1}{R_{C G P, j}}$ (16) $\frac{1}{Y^{*}} = \frac{1 - v_{i}^{2}}{Y_{i}} + \frac{1 - v_{j}^{2}}{Y_{j}}$ (17) $\frac{1}{m_{C G P}^{*}} = \frac{1}{m_{C G P, i}} + \frac{1}{m_{C G P, j}}$ (18) $δ_{n} = (R_{i} + R_{j}) - |r_{j} - r_{i}|$ (19) $γ_{n} = \sqrt[]{\frac{4 m_{C G P}^{*} k_{n}}{1 + {(π / l n \sqrt[]{e_{C G P}})}^{2}}}$ (20) $e_{C G P} = \sqrt[]{1 + k (e_{p}^{2} - 1)}$ (21) $u_{i j}^{n} = (u_{i j} \cdot n_{i j}) n_{i j}$ (22) $u_{i j}^{t} = u_{i j} - u_{i j}^{n}$ (23) $m a x (|F_{i j}^{t}|) = μ_{f} |F_{i j}^{n}|$ (24)

式中，kn为法向弹性系数；kt为切向弹性系数，且kt=kn；γn为法向阻尼系数；γt为切向阻尼系数，且γt=γn；δn为两碰撞颗粒间的法向重叠量；δt为两颗粒间的切向重叠量； nij 代表颗粒i质心指向颗粒j的法向单位向量； tij 代表颗粒间的切向单位向量； $u_{i j}^{n}$ 是颗粒i与颗粒j的法向相对速度； $u_{i j}^{t}$ 代表两颗粒间的切向相对速度；μf为滑动摩擦因数； $R_{C G P}^{*}$ 为粗颗粒的有效半径； $Y^{*}$ 为杨氏模量； $m_{C G P}^{*}$ 为粗颗粒的有效质量；uc为颗粒的特征速度； $R_{C G P, i}$ 和 $R_{C G P, j}$ 分别为粗颗粒i和粗颗粒j的半径；vi 和vj 分别为颗粒i和颗粒j的泊松比；Yi 和Yj 分别为颗粒i和颗粒j的杨氏模量； $m_{C G P, i}$ 和 $m_{C G P, j}$ 分别为粗颗粒i和粗颗粒j的质量；Ri 和Rj 分别为颗粒i和颗粒j的半径； ri 和 rj 分别为颗粒i和颗粒j的位置矢量；ep为真实颗粒的恢复系数。

本研究中固相颗粒为非球形的石英砂颗粒，故采用Hua等[31]提出的非球形颗粒的曳力模型来确定本研究系统内气固相间曳力。其中具体的曳力系数可表示为[41]

$β_{E r g u n - W e n - Y u - G a n s e r} = φ_{g p} β_{E r g u n} + (1 - φ_{g p}) β_{W e n - Y u - G a n s e r}$ (25) $β_{E r g u n} = \frac{180 μ_{g} ε_{p}^{2}}{ψ^{2} d_{a v e}^{2} ε_{g}} + \frac{2 ρ_{g} ε_{p} |u_{g} - u_{a v e}|}{ψ d_{a v e}}$ (26) $β_{W e n - Y u - G a n s e r} = \frac{3}{4} C_{D} \frac{ρ_{g} ε_{p} ε_{g} |u_{g} - u_{a v e}|}{d_{a v e}} ε_{g}^{- 2.65}$ (27) $φ_{g p} = \frac{a r c t a n [180 \times 2.0 (ε_{p} - ε_{s m f})]}{π} + 0.5$ (28) $u_{a v e} = \frac{\sum_{i = 1}^{N_{C G P, c e l l}} m_{i} u_{i}}{\sum_{i = 1}^{N_{C G P, c e l l}} m_{i}}$ (29) $d_{a v e} = {(\frac{\sum_{i = 1}^{N_{C G P, c e l l}} d_{C G P, i}^{3}}{N_{C G P, c e l l}})}^{1 / 3} / k$ (30) $C_{D} = \frac{24}{R e_{p} K_{1}} [1.0 + 0.1118 {(R e_{p} K_{1} K_{2})}^{0.6567}] + \frac{0.4305 K_{1}}{1 + \frac{3350}{R e_{p} K_{1} K_{2}}}$ (31) $R e_{p} = \frac{ε_{g} ρ_{g} d_{a v e} |u_{g} - u_{a v e}|}{μ_{g}}$ (32) $K_{1} = {(\frac{1}{3} \frac{d_{A}}{d_{V}} + \frac{2}{3} ψ^{- 0.5})}^{- 1} - 2.25 \frac{d_{V}}{D}$ (33) $K_{2} = 10^{1.8148 {(- l g ψ)}^{0.5743}}$ (34) $β_{i} = \frac{ε_{i} d_{a v e}}{ε_{p} d_{i}} β_{E r g u n - W e n - Y u - G a n s e r}$ (35) $β_{i} = β_{E r g u n - W e n - Y u - G a n s e r}$ (36) $F_{d r a g, i} = \frac{β_{i}}{ε_{i}} \times \frac{π}{6} d_{C G P, i}^{3} (u_{g} - u_{i})$ (37)

式中，εp为网格中颗粒的体积分数；Ψ为颗粒的球形度；mi 为粗颗粒质量； ui 为粗颗粒速度；k为粗粒化率；Rep为颗粒Reynolds数；K1和K2分别为颗粒斯托克斯形状因子与牛顿形状因子；dA为颗粒等投影面积球当量直径；dV为颗粒等体积球当量直径，且dA和dV两参数难以通过实验测得，故假设dA=dV=dsauter为颗粒Sauter平均直径；D为流化床的水力直径；εi 为网格内i类颗粒的体积分数；di 为颗粒的粒径。由于本研究体系内颗粒为单一粒径颗粒，故di = dave， εi =εp。

3.2 模拟设置

本研究所需方形气固挡板流化床与无挡板自由床的三维几何结构由ANSYS ICEM软件进行建立，模拟参数如表1所示。采用均匀的六面体结构化网格进行划分绘制完成对整个计算床体的网格表征，而后将网格导入开源计算流体力学软件OpenFOAM进行求解计算。流体计算网格用于求解CFD中计算流体的空隙率及其他物理量，对保证数值计算的精确性和稳定性具有重要意义。流体计算网格的尺寸需足够精细方可获得充足的尺度信息[41-42]，且为保证计算资源与时间成本的合理性，故取边长为10 mm的六面体网格作为后续模拟研究的基础，网格尺寸约为真实颗粒粒径的15倍，亦是粗粒化后粗颗粒的2~5倍，故可保证用颗粒中心算法可合理统计计算网格中的空隙率及固相速度等物理量。分别在几何模型的X、Y和Z三个方向上绘制30×30×220个均匀流体计算网格。此外，本文所采用的固相颗粒物料石英砂颗粒为Geldart B类颗粒，划分的流体计算网格足以获得网格无关性的模拟结果[43]。

表1 挡板床与自由床数值模拟参数

Table 1 Simulation parameters for the baffled and free fluidized bed

Parameter		Value
bed size	Lx ×Ly ×Lz /mm	300×300×2200
baffle size	lx ×ly ×lz /mm	60×300×10
particle	mean Sauter particle diameter/μm	595
	density/(kg/m3)	2906
	voidage at the minimum fluidization condition	0.47
	minimum fluidization velocity/(m/s)	0.33
	sphericity	0.86
	coarse-graining ratio	5
	coarse-grained particle number	3459863
	restitution coefficient	0.90
	friction coefficient	0.30
	rolling friction coefficient	0.01
	characteristic velocity/(m/s)	0.5
	Young’s modulus/Pa	1×108
	Poisson’s ratio	0.3
	time step/s	1×10-5
gas	density/(kg/m3)	1.2
	viscosity/(Pa·s)	1.8×10-5
	superficial gas velocity/(m/s)	0.4, 0.6, 0.9
	operating pressure/Pa	101325
	gas grid size/mm	10
	time step/s	1×10-4

新窗口打开| 下载CSV

图3显示了本文数值模拟的边界条件设置情况。流化床床层底部作为气相流体的入口，将其设置为速度进口边界条件，底部采用均匀进气，且速度大小依据操作气速设定。气相出口位于反应器几何模型的顶部，并设置为压力出口边界条件，且压力值与大气压相同。为模拟前期实验中二级旋风分离器及滤袋的作用，在反应器几何模型顶部设定不允许颗粒逸出，以保证模拟体系内的颗粒存量恒定。挡板内构件设置为与流化床壁面一致，设定为无滑移的边壁。为确保计算资源与时间成本的合理性，本研究选取模拟离散颗粒的粗粒化率k=5，且流体计算网格的尺寸是粗粒化后粗颗粒的2~5倍，可保证模拟计算的准确性[42, 44]。根据流化床内真实颗粒数目确定的粗粒化颗粒数为3459863个。

4 模拟结果与讨论

4.1 碰撞受力统计计算公式验证

挡板表面法向应力分量 σ33在颗粒与壁面的碰撞作用中占绝对主导地位，其统计结果的数量级远远大于其余8个分量。因此，挡板表面法向应力分量 σ33可作为评价挡板表面受颗粒碰撞作用大小的参数，其余8个分量可忽略。此外，由前面分析可知，挡板表面法向应力分量 σ33为挡板上、下表面受到颗粒碰撞的均布载荷差，即挡板受到的颗粒碰撞载荷为挡板下表面受到的碰撞均布载荷与上表面受到的碰撞载荷的差值，即

$q = σ_{33} = σ_{33}^{d o w n} - σ_{33}^{u p}$ (38)

式中， q 为挡板受到的颗粒碰撞均布载荷，方向竖直向上； $σ_{33}^{d o w n}$ 为挡板下表面受到的颗粒碰撞均布载荷； $σ_{33}^{u p}$ 为挡板上表面受到的颗粒碰撞均布载荷。

从图4中可以看出，流化床启动阶段内构件所受碰撞载荷的CFD-DEM模拟结果与实验结果随时间演变的趋势大体上一致，数值模拟统计的挡板受均布载荷变化曲线显示在流化气体刚通入的启动阶段会产生一个较大的脉冲峰值，而后挡板内构件所受的载荷减小并进一步波动变化。与实验值相比，启动阶段的第一个波峰及波谷值出现的时间更早，这说明应用本文所建立的载荷统计计算公式可以定性地成功模拟统计出床层由固定床向流化床转变过程中内构件受到非正常受力特性这一现象。此外，图中CFD-DEM数值模拟结果的脉冲峰值与实验值大小基本一致，说明模拟可以(半)定量复现实验中观测到的气固流化床反应器中密相床层流态转变时的受力特性。最后，在后续稳定阶段载荷波峰的出现与气泡尾流内颗粒与壁面作用紧密相关。综上所述，挡板内构件受颗粒碰撞载荷统计计算公式可以较为准确地预测气固流化床启动阶段内构件的受力特性。

图4

图4 流化床启动阶段内构件所受载荷实验结果与CFD-DEM模拟结果比较

Fig.4 Comparison of the experimental and simulation results of the stress exerted on the baffle during start-up of the fluidized bed

4.2 颗粒碰撞模型参数对受力特性的影响

考虑到内构件表面载荷主要由固相颗粒与内构件表面的碰撞作用贡献，在CFD-DEM数值模拟方法中，离散颗粒与挡板壁面的碰撞作用主要由DEM统计这些微观物理量场信息。DEM方法中，颗粒碰撞模型参数为碰撞作用的主要影响因素，颗粒杨氏模量、碰撞恢复系数、滑动摩擦因数和滚动摩擦因数这些碰撞参数对内构件表面与颗粒间的碰撞作用密切相关，并对挡板内构件表面载荷的影响尤为显著。因此，考察固相颗粒的杨氏模量、摩擦因数和恢复系数对挡板表面受碰撞力的影响有重要意义。

颗粒杨氏模量是固相颗粒物性的重要参数，表征了颗粒表面的硬度。杨氏模量的大小取决于固相颗粒自身的性质，其数值越大颗粒越不容易发生变形，反映了固相颗粒材料的刚性[45]。此外，在计算机数值模拟方法中，离散颗粒的杨氏模量作为模型碰撞参数是CFD-DEM数值模拟计算的一个重要参数。颗粒杨氏模量很大程度上影响数值模拟的计算量及计算效率，其决定了DEM中计算的时间步长。在研究气固两相间的时均流动特性时，颗粒的杨氏模量对统计计算的流动参数无显著影响，可以将杨氏模量数值设置小些以减少计算量，提高计算效率[46]。

对于密相床层的启动阶段，体系内颗粒-挡板壁面间的碰撞载荷统计则与杨氏模量的大小密切相关。基于本研究所选用的Hooke模型，法向弹性系数kn与杨氏模量直接相关。因此，为了探究内构件的受力特性，需要研究杨氏模量的大小对气固流化床启动阶段密相床层中挡板内构件表面受到的碰撞载荷的影响。固相颗粒物料仍为相同的B类石英砂颗粒，其杨氏模量的取值范围介于0.1 GPa和10 GPa量级之间，其中前者为研究气固流动特性统计时均参数常选用的经验值，而后者是真实石英砂颗粒的杨氏模量量级，而挡板内构件的杨氏模量与颗粒的杨氏模量设置相同。本小节选用的颗粒杨氏模量具体数值及对应DEM与CFD中的时间步如表2所示，模拟计算体系及其余模拟参数设置与前面小节模拟工况的设置一致。

表2 挡板床数值模拟参数

Table 2 Simulation parameters for the baffled fluidized bed

Parameter	Value
Young’s modulus/Pa	1×108	5×108	1×109	2×109	3×109	5×109	1×1010
DEM-time step/s	1×10-5	1×10-5	1×10-6	1×10-6	1×10-6	1×10-6	1×10-6
CFD-time step/s	1×10-4	1×10-4	1×10-5	1×10-5	1×10-5	1×10-5	1×10-5

新窗口打开| 下载CSV

图5显示了固相颗粒选取不同杨氏模量下密相床层中挡板内构件表面碰撞载荷随时间的演化情况。发现流化床自启动开始至出现第一个完整峰值脉冲载荷阶段，因颗粒杨氏模量的不同，峰值载荷的大小与峰值出现的时间均受到影响发生相应变化，而后载荷的变化不再受颗粒杨氏模量的影响。对图中启动阶段出现的第一个载荷峰值位置局部放大，如图5内附图所示。（1）在颗粒初始化自由堆积成密相床层的阶段(0 s之前)，颗粒堆积过程中体系逐渐密实，新力链不断生成。当床层高度稳定，此时体系足够密实，力链充分发展几乎不发生变化，挡板内构件表面载荷几乎稳定不变[47]。颗粒堆积的随机性与堆积时间的差异使得未通入气体前挡板上、下表面受到的颗粒的挤压碰撞作用不同，致使挡板内构件表面载荷存在一定的差异性。但是大体上，较大的颗粒杨氏模量使得颗粒-壁面间的法向接触力增大，致使固定床状态时挡板受到的碰撞作用载荷更大。（2）当体系内通入流化空气，启动阶段开始至第一个峰值脉冲出现(0~0.12 s)，峰值出现的时间随杨氏模量的增加而缩减，峰值载荷的大小随颗粒杨氏模量的增大而减小，且当杨氏模量增大至10 GPa量级时，此影响不再显著，峰值大小不再发生明显变化。说明颗粒杨氏模量的增加，加快了由固定床转化为流化床的速度，杨氏模量大的颗粒床层内挡板表面载荷率先达到峰值。但由于杨氏模量增大，颗粒的表面硬度更大、变形程度减小，致使颗粒的碰撞重叠量减小，碰撞接触时间明显缩短，使得内构件表面的碰撞作用载荷减小。此外，杨氏模量大的颗粒床层因颗粒硬度大、变形小，在流体通入时会更容易打破“自锁”状态，固相颗粒较快发生松动而流化，可能也是使得挡板表面受碰撞载荷峰值较小的原因。（3）当挡板表面受到的第一个峰值脉冲结束后(0.12 s之后)，内构件表面载荷变化基本不受颗粒杨氏模量的影响。体系进入正常流化阶段后，颗粒杨氏模量大的床层中颗粒-壁面的碰撞作用更剧烈，发生碰撞更为频繁，以及在弹性系数增加的共同影响下，在一定程度上抵消了因碰撞重叠量减小而使得内构件表面载荷减小的趋势。

图5

图5 不同颗粒杨氏模量下颗粒床层内挡板受到载荷随时间的变化

Fig.5 Effect of particle Young's modulus on the stress exerted on the baffle immersed in the fluidized bed

流化床启动阶段，未流化密相床层会被气节推起一定高度，使得挡板内构件受到较大的碰撞作用。而在这一短暂过程中，颗粒间的碰撞作用较为强烈，碰撞重叠量会更大，为保证模拟计算的精度，统计床层启动阶段颗粒最大重叠量与颗粒粒径的比值随时间的变化情况，如图6所示。从图中可以看出，每个工况下颗粒最大重叠量与粒径的比值都小于2%，说明启动阶段模拟计算的准确性。此外，随着颗粒杨氏模量的增大，最大重叠量与粒径的比值减小，也进一步说明了颗粒硬度增大，难以发生变形的特性。

图6

图6 不同杨氏模量下颗粒最大重叠量与粒径比值随时间的变化

Fig.6 Variation of the ratio of the maximum particle-particle overlap to particle size with time under different Young's modulus

颗粒碰撞恢复系数是颗粒的重要物性参数之一，可以反映颗粒碰撞后恢复到变形前初始状态的能力，即可以反映颗粒在接触碰撞过程中能量耗散的强弱[48]。颗粒摩擦因数是固相颗粒物料的一个基本物性参数，与接触面的粗糙程度有关。颗粒-颗粒与颗粒-壁面之间的摩擦会消耗颗粒的动能，进而影响体系内气固相间的交换量以及固相颗粒在切向方向上的运动情况[49]。颗粒滚动摩擦因数作为固相颗粒物料的重要物理力学参数之一，也会影响颗粒体系内形成的微观力学结构。颗粒间的滚动摩擦因数主要影响颗粒系统的堆积形态[50]，而在DEM碰撞模型中，颗粒滚动摩擦因数的取值一般远小于固相颗粒的直径[51]。因此，依据单一变量原则，分别选取0.80、0.85、0.90这三个值作为模拟所需的颗粒碰撞恢复系数；选取颗粒的滑动摩擦因数为0.30、0.35和0.40；取颗粒的滚动摩擦因数为0.01、0.03和0.05。其余参数的设置与之前保持一致，模拟结果如图7所示，可以看出碰撞恢复系数、滑动摩擦因数和滚动摩擦因数对挡板内构件的受力影响并不显著。

图7

图7 不同颗粒碰撞恢复系数(a)、滑动摩擦因数(b)、滚动摩擦因数(c)下颗粒床层内挡板受到载荷随时间的变化

Fig.7 Effect of particle restitution coefficient (a), particle sliding friction coefficient (b), and particle rolling friction coefficient (c) on the stress exerted on the baffle immersed in the fluidized bed

4.3 表观气速对受力特性的影响

流化床反应器的流化风速是控制流化床运行及气固两相流动的重要操作条件。因此，研究表观气速对挡板内构件受力特性的影响具有重要意义。基于4.2节关于碰撞模型参数对受力特性的影响研究，现选定一组关于颗粒属性的最佳模型参数，即颗粒的杨氏模量为1×1010 Pa，DEM中时间步长为1×10-6 s，CFD中时间步长为1×10-5 s，碰撞恢复系数为0.90，滑动摩擦因数为0.30，滚动摩擦因数为0.01。根据固相颗粒的属性及最小流化气速，现选取用于模拟研究的表观气速的一组值0.5、0.6、0.9 m/s。

图8为不同表观气速下，流化床启动阶段密相床层中挡板内构件受到的碰撞载荷瞬态变化情况。模拟结果显示，对于初始堆积床层，表观气速增大，相当于给予流化床的气速增量增加，挡板内构件的受力载荷强度也相应增大。当表观气速增加的幅度越大，内构件表面载荷的数值增长的幅度也越大。出现这一现象的原因在于表观气速增大，气速增量越大，流化床底部形成的气节高度更高，床层中固相颗粒的动量越大，对挡板内构件表面的撞击更为剧烈，因此载荷强度在数值上表现出更大的情况。说明在流化床床层启动阶段，表观气速的大小是影响挡板内构件受力的主要因素之一。

图8

图8 不同表观气速下颗粒床层内挡板受到载荷随时间的变化

Fig.8 Effect of superficial gas velocity on the stress exerted on the baffle immersed in the fluidized bed

上述模拟结果说明在实际流化床启动操作过程中，设置较小的表观气速可明显减小内构件受到的碰撞力。因此，要避免在较大的气速通入下完成床层启动操作，防止挡板内构件因瞬时受到较大的作用力而发生强度破坏，导致反应器装置无法正常运行。

4.4 颗粒粒径对受力特性的影响

流化床体系内，固相颗粒粒径的改变可使系统的最小流化速度发生变化，从而影响体系内的流态转变及颗粒碰撞动力学[49]。因此，采用不同粒径的颗粒物料作为床料进行数值模拟，基于单一变量原则，在表观气速分别按统一设定的工况和以最小流化气速倍数设定的工况下，研究颗粒粒径变化对密相床层中挡板内构件受力特性的影响。模拟研究所选用的固相颗粒仅粒径不同，其余的颗粒物性参数均一致（如颗粒数目）。

基于单一变量原则，将流化床内的流化气速采用统一设定，并将其设置为0.6 m/s，以此来研究相同表观气速下颗粒尺寸对内构件受力特性的影响。本节关于颗粒属性的模型参数设置与4.3节一致，其余参数参考表1设置。

图9显示了表观气速统一设定为0.6 m/s工况下，不同粒径的颗粒床层启动阶段密相床层中挡板内构件表面载荷强度的瞬态演变情况。从图中可看出随着颗粒粒径的增大，内构件表面受的碰撞载荷强度整体有减小趋势，尤其在第一个完整峰值载荷出现的区间，载荷强度有明显减小的趋势，在0.06 s位置表现更为明显。出现这一现象的原因可能在于，随着颗粒尺寸的增加，体系内的最小流化气速相应增加，而在相同表观气速下粒径较大的固相颗粒床层易出现流化不良的现象，床层中气泡的尺寸相对较小，气固两相的运动及相间作用程度减弱，挡板内构件受到的载荷作用也相应减小，这与前期实验中的现象基本一致。

图9

图9 不同颗粒粒径床层启动阶段挡板受到载荷随时间的变化

Fig.9 The variation of stress on the baffle immersed in the fluidized bed with time using different particle sizes

5 结论

本研究基于CFD-DEM数值模拟方法，基于刘对平[4]前期实验所用的三维方形冷模流化床装置作为模拟体系，以工业上广泛应用的单个斜片挡板内构件作为研究对象，建立了挡板内构件表面受力载荷强度的统计计算方法，并借助此统计计算公式分析研究了气固物性对挡板床启动阶段密相床层中单个水平斜片挡板内构件的受力特性的影响规律。主要结论如下。

（1）基于CFD-DEM模拟方法建立了气固流化床中浸没在密相床层的挡板内构件表面受颗粒碰撞载荷分布的统计计算公式，可用于定量分析内构件表面受到的动态碰撞载荷信号。通过与课题组前期实验测量数据对比，建立的内构件表面载荷统计计算公式可以定性和(半)定量复现实验中流化床启动阶段内构件受力特性的变化，为研究内构件在床层启动阶段流型转变过程中出现这一非正常受力现象的内在机理奠定基础。

（2）颗粒碰撞模型参数中，杨氏模量是影响内构件受力特性的关键性因素，而碰撞恢复系数、滑动摩擦因数和滚动摩擦因数对挡板内构件的受力影响并不显著。研究发现流化床自启动开始至出现第一个完整峰值脉冲载荷阶段，不同的颗粒杨氏模量下，内构件受力的峰值载荷大小与峰值出现的时间均受到影响发生相应变化（随着杨氏模量的增大，峰值载荷先逐渐变小后趋于稳定），而后阶段载荷的变化几乎不再受颗粒杨氏模量的影响。此外，在受杨氏模量影响阶段，随着杨氏模量的增大，第一个完整峰值载荷相应减小，当杨氏模量增大至颗粒物理真实值，峰值载荷不再发生变化。

（3）气固流化床体系内表观气速是影响内构件受力的主要因素之一。表观气速越大，密相床层中的挡板内构件受力载荷强度越大，与实验结论一致。为保障工业上反应器装置的正常运行和内部构件及其支撑结构的稳定性，要在相对较小的气速增幅下完成流化床的启动操作过程。

关键字：优秀论文

上一篇：磁场对湿颗粒流化床系统中介尺度结构影响机制研究
下一篇：光纤内窥图像法测量MTO催化剂表观形貌及其积炭量的实验研究

栏目分类

热门排行

推荐信息

期刊知识