气固两相流化床在能源和冶金等一些重要领域已被技术人员广泛应用^[1-2],对其进行研究具有重要意义,而且对流化过程操作条件的优化也具有指导意义^[3-4]。目前,流化床主要工作参数多按照实验结果或经验设定,而对流化过程的影响因素无法精确掌握,易造成流化效果不佳,影响最终结果。利用计算机对流化床进行数值模拟已经慢慢成为人们研究流化床内气固两相流动特性的主要技术手段。该方法不仅能够对流化床内两相流动特征进行精确地计算、预测与研究,而且对于一些用传统实验方法无法获取和难以测定的细节问题,其能以直观的方式进行展现^[5-9]。Prashant等^[10]模拟了小尺度致密气固流化床,研究了壁面附近和中部区域上、下床高度的固体速度差异。Maysam等^[11]分析流化床中的脉冲流,指出利用脉冲流可以提高流化的质量与效率。Wang等^[12]针对单个孤立颗粒,采用CFD-DEM方法对传热进行分析,指出流体流速和颗粒直径对对流换热的影响比传导换热的影响更显著。李九如等^[13]利用欧拉双流体模型构建了合理的气固两相流模型,模拟了鼓泡沸腾床的气固两相流动过程,分析了其流动特性。王坤等^[14]对气固流化床进行了数值模拟,使用欧拉-拉格朗日模型,分析了流化床内颗粒速度、粒度分布、压力与固含率。

国内外专家针对流化床流动状态模拟和实验做了诸多研究,但在数值模拟方面多是二维的,一些学者也开展过三维模拟,主要是基于欧拉-欧拉模型(TFM),将流化床中固相颗粒视作另一流体类型,改变了颗粒实际受力情况,由于欧拉双流体模型在分析颗粒运动特性中有一定局限性。笔者采用欧拉-拉格朗日模型,基于计算流体力学(CFD)和离散单元法(EDEM)对流化床气固两相瞬态流场进行三维数值模拟,研究不同进风速度下流化床内流场特点和颗粒运动,为探究流化床不同工况下最佳进风速度提供理论依据。

气相的质量和动量守恒方程^[15]为

式中:ε_g——孔隙率;

ρ_g——气相的密度;

u_g——气相速度矢量;

p——流化床内静压;

τ_g——气相黏性应力张量;

g——重力加速度;

F_DEM——气相与颗粒相对运动时运动的颗粒产生的作用在气相上的相互作用力。

F_DEM=∑β(u_p-u_g)+F_o,(3)

式中:β——气相与颗粒的能量传递系数;

F_o——其他作用力。

此质量与动量守恒方程也可用于计算颗粒相。

对于颗粒相的相关运动,将在DEM求解器中进行计算,主要包括颗粒的平移和旋转,主要通过计算牛顿第二定律^[16],即

式中:k——颗粒i相接触的颗粒总数量;

F_n,ij——i、j两颗粒碰撞时的法向接触力;

F_t,ij——i、j两颗粒碰撞时的切向接触力;

F_p,i、F_d,i、F_τ,i——流化风所施加的压力梯度力、曳力以及黏性力;

I_i——单颗粒的转动惯量。

软球模型作为目前应用最广的颗粒接触碰撞模型为

式中:k_n、k_t——法向和切向的等效弹簧刚度;

c_n——法向阻尼系数;

Δu_n、Δu_t——法向和切向偏移量。

等效的弹簧刚度和阻尼系数基于赫兹-梅林颗粒碰撞模型^[17]得到。

以FL-5型流化床平台作为研究对象,结构如图1所示。其主要由流化床床体、进气管和出口管三个部分组成。将流化床床体分为三个区域,上层为抖袋及其他部件,中间部分为颗粒流化区域,下方为物料颗粒。为了减少计算量省略掉物料下方孔板、上层抖袋及其他一些零部件,利用三维建模软件对流化床流体域进行建模。

图1 FL-5流化床结构示意
Fig.1 FL-5 schematic of fluidized bed structure fig

采用Fluent Meshing对流体计算域进行网格划分^[18],为了确保模拟精度,内部采用结构化的六面体网格,外部设置三层多面体网格作为边界层。设立流化床的边界条件包括底部速度入口、顶部压力出口及壁面边界。

在EDEM中考虑颗粒与颗粒、颗粒与壁面间的相互作用力,假设固相颗粒均为单个均质球体,采用软球模型,接触模型选取Hertz-Mindlin无滑移模型,在模拟中固体颗粒不产生任何变化。在CFD中进行瞬态模拟,湍流模型使用SST k-ω模型,流化床进气管设置6种不同进风速度,出气管设置压力出口为0 Pa,容器壁面为无滑移壁面,采用有限体积法对气相进行求解,压力-速度耦合采用Simple算法实现。CFD-DEM耦合具体模型参数设置分别为:颗粒泊松比0.23; 密度1 000 kg/m³; 剪切模量1×10^-7 Pa; 恢复系数0.5,静摩擦系数0.5,动摩擦系数0.01; 颗粒半径为3 mm; 颗粒数量为20 000个; 时间步长3×10^-5 s,进风速度分别为33、35、37、39、41、43 m/s; 曳力模型为Gidaspow模型; 迭代次数600; 时间步长0.01 s。

为了探究流化床不同进风速度对气固流化床两相流化流动特性的影响,分别用33、35、37、39、41、43 m/s,6种进风速度对气固流化床进行了6 s的数值模拟。

进风速度为37 m/s时x-z平面与y-z平面上0～6 s不同时间的颗粒瞬时运动状态如图2所示。

图2 v=37 m/s时不同时刻颗粒运动状态
Fig.2 Particle motion state at different time for v=37 m/s

从图2可以看出,初始状态时,空气通过进气管以一定速度进入流化床内部,扰动气固两相流场,突然增加的气体轴向流动速度,使颗粒整体随气流开始向上不断翻腾运动; 在t=0.1 s时,底部料仓沉积的颗粒受到进气管瞬间进入气流的冲击作用,几乎全部被吹起,运动非常激烈,且中间上部的颗粒运动速度最快; 同样在0.1～1 s时间段内,可以发现颗粒运动状态也发生了极大的变化,较多颗粒回落至流化床底部; 与t=1 s时相比,t=2 s时颗粒运动平稳了许多,开始时被气流冲击到上方的颗粒都已平稳下落,流化状态相对稳定。t=3 s之后,流化状态达到稳定水平。

从图2b可以看出,受风力影响,颗粒大多从左下方靠近壁面以较大速度向上运动,而越向上流场的湍动能提供的升力就会越小,所以在向上运动过程中颗粒速度逐渐减小。当流场湍动能提供的升力无法支撑颗粒重力时,颗粒就会四散开来向下运动,而且在料仓的右下角有一部分颗粒会静止不动,始终没有参与流化过程。造成这种状态的原因主要有两个,一是由于流化床自身几何形结构决定,流化床料仓呈下小上大倒锥形,通过料仓底部进入的风由于流化区域结构横截面不断增大,致使风力不断减小,当风力不足以支撑物料重力时物料颗粒就会四散下落; 二是因为流化床进气管的位置。由图2b可见,进气管位于流化床底部右侧,那么高速气流通过进气管进入流化床内部到达底仓位置时左侧的气流流速就会比右侧气流流速大,图2b中1 s时左侧颗粒被高高吹起,而右侧仅有一部分处在下落过程,甚至右下角一部分颗粒在刚开始被吹起又下落至底部后,会静止不动; 而这恰巧又加剧了这一现象的产生,堆积在右下角的一小部分颗粒会一定程度阻挡气体流动,所以气体由于压力原因会进一步从左侧进入料仓,对颗粒进行流化,直至达到相对平稳的状态。

由图2可见,当流化状态达到稳定后,速度较大的区域有两个,一个是进气管对侧左下角区域,由于进气管位置原因此处风速最大,所以能够带动附近颗粒以较快速度向上运动; 另一个是进气管同侧未进入流化状态颗粒的上方区域,同样由于结构问题,首先右下角一部分颗粒在流化过程中几乎处于静止状态,会阻挡气体的进入,无法为该区域提供较大升力,所以该区域颗粒只能下落,其次流场在右边靠近壁面进入流化床后,由于出气管尺寸有限,并不能及时流出,所以流场会向右下方流动,在流化床内形成环形流场,会带给颗粒一个向下的力,最终致使此区域颗粒会快速向下运动。当颗粒流动达到相对稳定后,颗粒的速度在运动过程中会有一条明显的分界线,线的左侧颗粒向上运动,线的右侧颗粒向下运动,而在线附近区域,气流带给颗粒的升力与颗粒的重力相近,所以速度较小。

不同进风速度条件下,物料颗粒平均速度随时间变化曲线如图3所示。从图3可以看出,风机开始工作的瞬间,颗粒猛然间受到强烈气流的影响速度产生极大变化,在料仓底部由静止状态加速到喷动速度,高速气流与颗粒接触时,产生喷动轴线方向力和颗粒切线方向力,引起颗粒旋转上升运动,当颗粒到达一定高度后,出现喷涌,随着喷涌高度的增加,床层压强降低,在重力和减弱的驱动气流的合力作用下,颗粒速度降低,颗粒开始在压强小的区域回落,整体速度减小。随后率先回落到底部的颗粒由于高速气流的持续存在,会再次加速向上运动,与此时正在下落的部分颗粒发生碰撞,最终大部分颗粒会一同再次向上喷涌,造成了颗粒平均速度的第二次攀升。随着时间的推移,经过两次剧烈波动后,颗粒运动速度波动越来越小,逐渐达到平稳状态。

图3 不同进风风速下物料颗粒平均速度随时间变化曲线
Fig.3 Curves of average velocity of material particles with time under different inlet wind speeds

整体来看,在6种不同进风速度的作用下,物料颗粒的平均速度曲线趋势走向是相同的,随时间推移呈波浪式前进状态,且波动具有一定的周期性,随进风速度增加颗粒运动周期变长,但都逐渐趋于平稳。不仅如此随着进风速度的增加,颗粒平均运动速度也会变大,颗粒的速度波动也越来越激烈。显然进风速度越大,对整个流场的流动造成较大的扰动就会越大,流动变得复杂, 增加了流场的湍流混沌特性,气相可以带动颗粒更快进行无规则运动。

流化床需要在风力作用下使内部物料充分进行混合,但是进风速度太大会使得流化床床层高度过高,床内流化的物料与上层区域抖袋及其他结构接触,容易堵塞抖袋,导致抖袋返喷次数的增加和物料的损失,同时也会使上层压强增大,影响流化效果。而进风速度太小会使流化床内部存在涡流,使得部分床层内颗粒静止或者参与混合机会减少造成混合不充分。

不同风速下底仓、流化区域、上层区域物料质量随时间分布曲线,如图4所示。由图4可见,流化状态下颗粒主要在底仓和流化区内运动。除开始时高速气体与颗粒床层接触的瞬间产生了较大波动,造成颗粒质量流动不稳定外,随时间推移质量分布都越来越平稳,仅存在较小波动。由图4a可以看出随进风速度增加,在底仓的颗粒质量逐渐减少,底仓堆积的颗粒越少,进入流化状态的颗粒就越多,流化效果就会相对较好。由图4b可知,当达到平稳后,在流化区域的质量并非是随进风速度的增加而增大,进风速度过大或过小都会使在流化区域的颗粒质量降低,在当前模拟条件下,进风速度为39 m/s时流化区域颗粒质量最大。由图4c可知,随着进风速度不断增加,到达上层区域的颗粒会越多,说明随进风风速增大颗粒流化的床层高度会越高,但根据流化床内部实际情况,在上层区间有抖袋及一些其他装置,流化过程中并不希望床层高度过高。综合分析在当前模拟条件下,进风速度为39 m/s左右时流化效果相对较好。

图4 不同进风速度下物料质量随时间分布曲线
Fig.4 Distribution curve of material mass with time at different air inlet speeds

在EDEM颗粒仿真运动过程中,从填充的颗粒中,随机选取5个颗粒作为示踪粒子,对不同进风速度下颗粒运动轨迹进行分析,示踪颗粒运动轨迹如图5所示。

图5 不同风速下颗粒运动轨迹
Fig.5 Trajectory of particles at different wind speeds

从图5可以看出,静止颗粒在高速气流吹动下上下翻腾^[19],由于结构原因颗粒在进风速度的升力作用下由左侧快速向上运动,达到一定高度后升力减弱,颗粒就会向四周压力相对小的区域扩散,逐渐回落,形成左侧上升右侧下降的环形结构^[20],不仅如此,在速度较大的近壁面两侧由于湍流动能和重力的综合影响,也会使颗粒在壁面做回旋的环形运动。落回底仓后,会再次被强烈的风速沿左侧壁面吹起,以此往复,直至达到相对平衡状态,在此过程中颗粒在流化床内被不断混合。对比不同风速颗粒轨迹分析,进风速度较小时,颗粒大多在底仓运动,无法被吹起达到流化状态,底部及两侧面湍流较多,颗粒运动不规则,混合不够均匀。随着进风速度的不断增加,流化的环形结构逐渐清晰,形成大循环,在底仓的颗粒渐少,进入流化区域的颗粒更多,颗粒流动稳定性更强,流化状态更好,混合更加均匀; 但当进风速度继续增大时,颗粒运动速度会越大,运动状态也越复杂,颗粒喷出增多,飞溅现象明显,床层高度会随之增加,流化区域的颗粒会进入上层区域,对颗粒的流化状态造成破坏。

采用CFD-DEM耦合对流化床内气固两相流进行数值模拟,研究了不同进风速度条件下对流化床内颗粒运动状态的影响。

(1)由于FL-5型流化床结构原因,达到稳定后,物料颗粒会跟随气流方向从进气管对侧快速上升,达到流化高度后,从两侧或进气管同侧位置四散下落,形成特殊的环形流化结构。结合物料颗粒平均速度分析,流化达到稳定后颗粒运动会具有一定的周期性,且进风速度越大周期会越长。

(2)根据对流化过程中颗粒质量分布及运动轨迹分析,进风速度过大或过小都会影响颗粒的流化效果,速度过小,颗粒多在侧面及底部运动,无法对颗粒进行较好混合; 速度过大,颗粒飞溅明显,进入上层区域的颗粒会对流化状态造成破坏。综合当前模拟条件,进风速度为39 m/s左右时,流化效果相对较好。