埋管流化床内湿颗粒流动及混合特性的CFD-DEM数值模拟

刘道银，宋诚骁，王铮，马吉亮，陈晓平



埋管流化床内湿颗粒流动及混合特性的CFD-DEM数值模拟

刘道银，宋诚骁，王铮，马吉亮，陈晓平

（东南大学能源与环境学院，能源热转换及其过程测控教育部重点实验室，江苏南京 210096）

掺液滴的气固流化床，越来越多地应用于颗粒材料合成和改性工程，但是液滴使颗粒具有一定的湿度。湿颗粒具有黏性，与干燥颗粒流化特性有很大差异。本文应用迟滞碰撞模型，分别在颗粒碰撞的加载和卸载阶段考虑液体黏弹阻力（lubrication force）和液桥力引起的动能耗损，建立湿颗粒流化床CFD-DEM模型。对含有埋管的湿颗粒流化床进行数值模拟研究，比较了不同黏性条件下颗粒的流动及混合特性。发现随着颗粒黏性的增大，颗粒聚团加剧，气泡边界变得粗糙并逐渐转变为气体沟流，最终导致流化失效，埋管在床中起骨架作用，促进了沟流的形成。随着颗粒黏性的增大，颗粒的相对运动被阻碍，颗粒速度降低，床层压降逐渐减小，颗粒混合速度减慢。

湿颗粒流化；迟滞碰撞模型；颗粒聚团；颗粒混合

在众多的应用气固流化床的物理和化学工艺中，经常伴有喷雾过程，例如催化剂颗粒制备、药品颗粒制备、食品颗粒加工、废液焚烧、烯烃聚合反应等。在此类掺液滴气固流化床中，液滴使颗粒具有一定的湿度，湿颗粒具有黏性，可能导致颗粒团聚、流化不良甚至失流化现象。

湿颗粒流化特性已经引起了较多的关注。在实验研究方面，ZHOU等[1]结合图像法和声频测量技术研究了不同含液量气固流化床的流体力学特性。XU等[2]利应用图像法研究了湿颗粒条件下喷动流化床的流型、最小喷动速度、最大喷射高度等。WANG等[3]运用PIV技术研究了鼓泡流化床中湿颗粒流化形态，并与模型预测结果进行对比。MA 等[4]研究了黏性颗粒流化床中气泡特性（气泡形状、数量、直径等），研究气泡向气体沟流转变的规律。

在数值模拟研究方面，应用较多的是在DEM模型（discrete element method）[5]基础上考虑颗粒黏性，主要有三类思路：①黏性力与接触力线性叠加；②黏性力归结到碰撞恢复系数；③应用迟滞碰撞模型考虑黏性作用。

第一类思路将液桥力和碰撞力叠加。MIKAMI等[6]归纳出静态液桥力的回归形式，在弹簧阻尼碰撞模型的基础上，沿颗粒间法向叠加静态液桥力，模型可以预测聚团流化现象。在MIKAMI湿颗粒碰撞模型基础上，许多学者对不同对象湿颗粒体系，如滚筒、混合器、料仓、流化床等，进行了数值模拟研究。LIM等[7]研究了不同流化风速下湿颗粒流化床内的颗粒混合特性。WANG等[8-9]模拟流化床提升管中干湿颗粒的流化形态、聚团特性。HE等[10]对流化床中干湿颗粒的流化特性进行了研究。此类模型不能考虑液体黏度对颗粒碰撞黏结特性，因为静态液桥力由液体表面张力系数和液体含量决定，与液体黏性无关，动态液桥力虽然与液体黏性关系较大，但是在流化床条件下通常可以忽略。实验和理论研究均表明，颗粒相互靠近，排挤颗粒间的液体时产生黏弹阻力（lubrication force），对于高黏性液体而言尤为如此，因此直接忽略液体黏性的影响存在一定的局限性[11]。

第二类思路将湿颗粒碰撞动能耗损归结到碰撞恢复系数，然后用湿恢复系数（wet restitution coefficient）替代弹簧阻尼碰撞模型中的碰撞恢复系数[12]。VAN BUIJTENEN等[13]研究了湿恢复系数对喷雾流化床造粒过程中颗粒相互作用的影响。DARABI等[14]提出了一种湿恢复系数模型并对流化床进行数值模拟研究。该模型的优点在于同时考虑了液体黏性和液体表面张力对碰撞的影响，但是湿恢复系数很小时，弹簧阻尼碰撞模型中阻尼力过大，导致在加载阶段接触力可能减小的现象，偏离了实际物理过程；而当恢复系数接近0时，模型的求解存在奇点，导致该模型无法复现颗粒聚团 现象。

第三类思路运用迟滞碰撞模型[15]分阶段考虑黏弹阻力和液桥力对碰撞的影响。DAVIS等[16]研究了黏弹阻力对颗粒碰撞的影响，并提出了基于Stokes量纲数为1的碰撞恢复系数表达式。MA等[17]研究了湿颗粒的碰撞过程，区分了颗粒靠近阶段的黏弹阻力耗损和分离阶段的液桥力耗损。应用迟滞碰撞模型描述湿颗粒碰撞过程，加载过程考虑黏弹阻力耗损，卸载过程考虑液桥力耗损，既可以预测碰撞恢复系数为0时的情况（即颗粒发生黏结），同时又能考虑液体黏性和表面张力对动能耗损的影响。

常规颗粒在含埋管流化床中的运动特性已有很多研究[18-21]，但是针对湿颗粒的含埋管流化床的研究目前较为少见。本文作者将迟滞碰撞模型与CFD-DEM模型耦合，对含埋管湿颗粒流化床进行数值模拟研究，研究不同黏性和流化风速条件下的湿颗粒流动及混合特性（流动形态、颗粒速度、床层压降、混合行为）。限于计算量，本文作者针对含约20万个颗粒的二维流化床开展CFD-DEM数值模拟研究，主要目的在于从定性的角度考察分析湿颗粒流化特性以及模型适用性，同时为将来开展湿颗粒三维流化床CFD-DEM定量研究和分析奠定基础。

1 模型介绍

1.1 数学模型

在CFD-DEM模型[22-24]中，气相用Navier-Stokes方程描述，颗粒相用牛顿第二定律跟踪体系内每个颗粒的受力和运动。对于每个颗粒，考虑颗粒自身重力、周围颗粒对它的接触力以及气体对颗粒的曳力和压力梯度力。气固两相之间采用双向耦合，曳力应用Eugun和Wen&Yu关联式。CFD-DEM模型的接触力通常采用线性弹簧阻尼碰撞模型描述，而本文作者采用的接触力模型为迟滞碰撞模型，下面详细介绍。

1.2 湿颗粒碰撞模型

结合一对湿颗粒碰撞过程的接触力与形变量关系曲线（图1），对模型介绍如下[25]。图1所示中轴表示接触力，大于0表示颗粒间作用力为排斥力，小于0为吸引力。整个碰撞过程分为靠近、加载、卸载和分离。

阶段1：加载过程从图1中的点开始。当颗粒接触并开始加载时，颗粒持续减速直到速度为0（点）。段的斜率即加载弹性模量（kload）。BC段颗粒的受力情况见式(1)。

F= kload+p(1)

阶段2：在颗粒速度减小至0以后，颗粒在弹性力作用下开始反弹，卸载过程开始。卸载过程进行至点时，塑性形变量为0。段的斜率即卸载弹性模量（kunload）。段颗粒的受力情况见 式(2)。

F= kunload(–0)p(2)

k,unload与k,load的关系用黏弹恢复系数（wet）表示，wet=(k,load/k,unload)1/2。k,unload和wet是输 入参数。wet由数（）的函数表示[26]，本文取用常见形式wet=exp(–35/)。=pp/(6π2)中，p为颗粒质量；为液体黏度；为颗粒 半径[27]。

阶段3：卸载过程沿线继续进行，直到点。点以后弹性力消失，只有液桥力继续作用于颗粒之间。液桥力可以由液体表面张力系数和液体含量函数确定[7]。在本文中，液桥力的最大值由数（）给定。数定义为液桥力与颗粒重力的比值，如式(3)。

式中，为液体表面张力系数；p为颗粒 密度。

阶段4：分离阶段从点到点，在此过程中，颗粒没有直接接触，但是液桥力仍然存在，直到颗粒间距离大于极限液桥距离c（点）[7]。

除了上述两种由液体产生的能量耗损外，部分能量耗损以波等形式，耗散于空气，采用干碰撞恢复系数（dry）表示。将颗粒发生碰撞黏结的临界点记为临界碰撞速度（p,cr）。p,cr越大，表示颗粒黏性越大，则越容易发生黏结。

2 结果和讨论

2.1 模拟条件及工况介绍

模拟对象为准二维流化床，床体尺寸、颗粒物性等参数按照MA等[28]的实验条件设置。床体尺寸30cm×90cm，床体深度为1倍颗粒直径，颗粒数236148个，参数详见表1。流化床含18根埋管，埋管半径1cm，埋管布置高度分别为10cm、16cm、22cm、28cm和34cm，具体位置如图2所示。模拟物理时间12s，为了避免颗粒初始堆积状态对流化的影响，3s至12s的数据用于统计和分析。模型输入参数包括量纲为1液桥力数、液体黏度和弹性模量kunload。数由液体表面张力和含液量确定。在本文模拟中，固定和kunload，改变数和流化风速g，共12个计算工况，见表2。

2.2 流动形态及压力特征

图2为g=2.5mf时3种黏性条件下的瞬时流动形态模拟结果。随数增大，颗粒间黏性力作用越来越明显，颗粒团聚尺寸逐渐增大，气泡由较规则的圆形和帽形逐渐发展为具有分叉、边界粗糙的不规则形状，高黏性时气泡消失，形成稳定的气体沟流。MIKAMI[7]、HE[10]等的模拟结果也反映类似的气泡形态特性。此外，由颗粒颜色看到，黏性增加导致颗粒整体运动速度降低。气体沟流的出现和颗粒运动速度的降低，均表明流化床内气固流动受到抑制。埋管起到了骨架的作用，使得高黏性条件下颗粒团聚可能阻塞于埋管间、埋管与壁面间，从而促进了流化失效。在MA等[28]的实验中，黏性力的引入由多聚物涂层法实现，黏性力的大小与床温呈近似线性关系，无多聚物涂层即无黏性，有多聚物涂层、床温36℃即低黏性，有多聚物涂层、床温44℃即高黏性。将本文的模拟结果与MA等[28]的实验结果（图3）定性对照，黏性颗粒引起的聚团、气体沟流等特征基本与实验相符。

表1 数值模拟物理参数表

表2 数值模拟工况表

图4为g=2.5mf时3种颗粒黏性条件床层压降动态曲线。按照床料重量计算的床层压降理论值为3635Pa，模型计算得到的各颗粒黏性条件的床层压降与理论值接近，说明了模型合理性。图5为时均床层压降图，其中图5(a)对比不同黏性条件，图5(b)对比不同流化风速条件。随着颗粒黏性的逐渐增大，床层压降逐渐减小，并且颗粒黏性导致床层压降减小这一现象，在不同流化风速条件下都可以观察到。这是由于黏性增大导致颗粒聚团，并引起流化不良和气体沟流现象，导致床层压降下降。

图3 Ug=2.5Umf时MA[28]等实验拍摄不同颗粒黏性的瞬时流动形态图

2.3 颗粒及气体速度

图6为g=2.5mf时3种颗粒黏性条件下颗粒时均速度（p）矢量图，直观反映颗粒在床内的循环情况。干颗粒的循环状态清晰明显，随着黏性的增加，颗粒速度降低，颗粒循环逐渐变弱，=56时接近“死床”。由于埋管的存在，颗粒在床内大循环的同时有绕埋管小循环趋势，低黏性时较明显。

图7为g=2.5mf时3种颗粒黏性条件下时均气体速度（gl）矢量图。由于气泡集中在密相区床层中部，干颗粒条件下表现为床层中部gl较大，壁面附近较小。随着黏性增大，床层中部的时均气速减小，表明气固流动受到一定抑制，高黏性条件下沟流稳定存在。

图8为床层20cm高度处，时均颗粒水平分速度（px）和轴向分速度（py）的分布情况，其中图8(a)和8(b)对比不同颗粒黏性情况。随数增大，环核流动趋势越来越弱，在低黏性条件下（=6）环核流动尚可分辨，随着黏性加强（=17以上工况）环核流动基本消失。此外，由于采样点位置接近高度22cm处的埋管，可观察到该高度的4处埋管（=6cm，12cm，18cm，24cm）对颗粒整体环核流动的分解作用，其中=12cm和=18cm的埋管对颗粒速度方向改变效果最为明显。图8(c)和8(d)对比不同流化风速情况，可以看到提高风速可一定程度改善黏性颗粒流化，促进床内颗粒循环运动。

图6g=2.5mf时不同颗粒黏性条件下时均颗粒速度矢量分布（为了便于观察，每4个网格显示颗粒速度）

图7 为g=2.5mf时不同颗粒黏性条件下时均气体速度矢量分布（为了便于观察，每3个网格显示颗粒速度）

2.4 颗粒混合

图9为g=2.5mf时不同黏性条件的颗粒混合时序图，选取某一时刻为起始时刻0，时序间隔为2s，其中蓝色颗粒为示踪颗粒。为研究颗粒在水平方向的混合情况，将流化床从左到右等分10个区域，依次命名为L1～L10，示踪颗粒初始布置在L4和L5区域。同理，为研究轴向混合情况，将床体从下至上每3cm共划分10个区域，依次命名为V1～V10，示踪颗粒初始布置在V6和V7 区域。

图10为示踪颗粒在床内各区域的浓度分布曲线图，(a)、(b)为干颗粒，(c)、(d)为=17的湿颗粒。结合混合时序图与示踪颗粒浓度分布图可以发现，干颗粒在两个方向的混合明显快于湿颗粒，且随黏度增加，混合速率进一步降低。当颗粒黏度较高时，颗粒的混合进度几乎停滞，表明流化严重不良。此外，干颗粒水平方向达到混合均匀，大约需要5.5s，而竖直方向达到混合均匀，需要约2.5s，所以颗粒水平混合慢于轴向混合，该规律在黏性颗粒流化时也能观察到。这是由于流化床颗粒竖直方向运动比竖直方向运动更剧烈。

3 结论

应用CFD-DEM数值模拟耦合迟滞碰撞模型，研究埋管流化床中湿颗粒的流动及混合特性，主要结论如下。

（1）颗粒间的黏性力阻碍颗粒的相对运动，随着黏性力增大，更多颗粒形成聚团，并最终导致气泡转变为气体沟流，与实验结果对照定性合理。埋管起骨架作用，加速流化失效。时均床层压降随着颗粒黏性增大而减小。

（2）干颗粒床层中部的时均气体速度较壁面处大，高黏性条件下时均气速分布可反映稳定沟流的产生。埋管导致颗粒运动发生较大幅度的转向，并将颗粒环核流动分解为若干小环。颗粒速度随着颗粒黏性增大逐渐降低，环核流动逐渐消失，床内颗粒循环被破坏。

（3）颗粒混合速度随颗粒黏性增大而减小，高黏性时混合难以均匀。颗粒沿水平混合的速率通常大于竖直方向。

[1] ZHOU Y F，SHI Q，HUANG Z L，et al. Effects of liquid action mechanisms on hydrodynamics in liquid-containing gas-solid fluidized bed reactor[J]. Chemical Engineering Journal，2016，285：121-127.

[2] XU H B，ZHONG W Q，YU A B，et al. Spouting characteristics of wet particles in a conical-cylindrical spouted bed[J]. Industrial & Engineering Chemistry Research，2015，54（40）：9894-9902.

[3] WANG T Y，HE Y R，TANG T Q，et al. Experimental and numerical study on a bubbling fluidized bed with wet particles[J]. AIChE Journal，2016，62：1970-1985.

[4] MA J L，LIU D Y，CHEN X P. Bubble behaviors of large cohesive particles in a 2D fluidized bed[J]. Industrial & Engineering Chemistry Research，2016，55（3）：624-634.

[5] TSUJI Y，KAWAGUCHI T，TANAKA T. Discrete particles simulation of 2-dimensional fluidized bed[J]. Powder Technology，1993，77（1）：79-87.

[6] MIKAMI T，KAMIYA H，HORIO M. Numerical simulation of cohesive powder behavior in a fluidized bed[J]. Chemical Engineering Science，1998，53（10）：1927-1940.

[7] LIM C，WEE E，BENG R，et al. Mixing behaviors of wet granular materials in gas fluidized bed systems[J]. AIChE Journal，2013，59（11）：4058-4067.

[8] WANG M，ZHU W B，SUN Q Q，et al. A DEM simulation of dry and wet particle flow behaviors in riser[J]. Powder Technology，2014，267：221-233.

[9] 王猛. 循环流化床提升管内气固两相流的DEM数值模拟[D]. 哈尔滨：哈尔滨工程大学，2013.

WANG M. A DEM numerical simulation of gas-solid two-phase flow in circulating fluidized bed riser[D]. Harbin：Harbin Engineering University，2013.

[10] HE Y R，PENG W G，WANG T Y，et al. DEM study of wet cohesive particles in the presence of liquid bridges in a gas fluidized bed[J]. Mathematical Problems in Engineering，2014，1：1-14.

[11] SEVILLE J P K，CLIFT R. The effect of thin liquid layers on fluidisation characteristics[J]. Powder Technology，1984，37（1）：117-129.

[12] SUTKAR V S，DEEN N G，PATIL A V，et al. CFD-DEM model for coupled heat and mass transfer in a spout fluidized bed with liquid injection[J]. Chemical Engineering Journal，2016，288：185-197.

[13] VAN BUIJTENEN M S，DEEN N G，HEINRICH S，et al. A discrete element study of wet particle-particle interaction during granulation in a spout fluidized bed[J]. Canadian Journal of Chemical Engineering，2009，87（2）：308-317.

[14] DARABI P，POUGATCH K，SALCUDEAN M，et al. DEM investigations of fluidized beds in the presence of liquid coating[J]. Powder Technology，2011，214（3）：365-374.

[15] LUDING S. Cohesive，frictional powders: contact models for tension[J]Granular Matter，2008，10（4）：235-246.

[16] DAVIS R H，RAGER D A，GOOD B T. Elastohydrodynamic rebound of spheres from coated surfaces[J]. Journal of Fluid Mechanics，2002，468：107-119.

[17] MA J L，LIU D Y，CHEN X P. Experimental study of oblique impact between dry spheres and liquid layers[J]. Phys. Rev. E. Stat. Nolin Soft Matter Phys.，2013，88（3）：033018.

[18] YATES J G，RUIZ-MARTINEZ R S，CHEESMAN D J. Prediction of bubble size in a fluidized bed containing horizontal tubes[J]. Chemical Engineering Science，1990，45（4）：1105-1111.

[19] HULL A S，CHEN Z，FRITZ J W，et al. Influence of horizontal tube banks on the behavior of bubbling fluidized beds：1. Bubble hydrodynamics[J]. Powder Technology，1999，103（3）：230-242.

[20] ASEGEHEGN T W，SCHREIBER M，KRAUTZ H J. Investigation of bubble behavior in fluidized beds with and without immersed horizontal tubes using a digital image analysis technique[J]. Powder Technology，2011，210（3）：248-260.

[21] ASEGEHEGN T W，SCHREIBER M，KRAUTZ H J. Numerical simulation and experimental validation of bubble behavior in 2D gas–solid fluidized beds with immersed horizontal tubes[J]. Chemical Engineering Science，2011，66（21）：5410-5427.

[22] VAN WACHEM B G M，VAN DER SCHAAF J，SCHOUTEN J C，et al. Experimental validation of Lagrangian–Eulerian simulations of fluidized beds[J]. Powder Technology，2001，116（2）：155-165.

[23] LIU D Y，BU C S，CHEN X P. Development and test of CFD-DEM model for complex geometry：a coupling algorithm for fluent and DEM[J]. Computers & Chemical Engineering，2013，58：260-268.

[24] LUO K，YANG S L，FANG M M，et al. LES-DEM investigation of the solid transportation mechanism in a 3-D bubbling fluidized bed. part I: Hydrodynamics[J]. Powder Technology，2014，256：385-394.

[25] LIU D Y，VAN WACHEM B G M，MUDDE R F，et al. An adhesive CFD-DEM model for simulating nanoparticle agglomerate fluidization[J]. AIChE Journal，2016，62（7）：2259-2270.

[26] LEGENDRE D，ZENIT R，DANIEL C，et al. A note on the modelling of the bouncing of spherical drops or solid spheres on a wall in viscous fluid[J]. Chemical Engineering Science，2006，61（11）：3543-3549.

[27] DAVIS R H，SERAVSSOL J M，HINCH E J. The elastohydrodynamic collision of 2 spheres[J]. Journal of Fluid Mechanics，1986，163：479-497.

[28] MA J L，LIU D Y，CHEN X P. Bubbling behaviors of cohesive particles in a two-dimensional fluidized bed with immersed tubes[J]. Particuology，2017，31：152-160. http:// dx.doi.org/10.1016/ j.partic. 20106.005.

CFD-DEM simulation of wet particles flow and mixing behavior in fluidized bed with immersed tubes

LIU Daoyin，SONG Chengxiao，WANG Zheng，MA Jiliang，CHEN Xiaoping

（Key Laboratory of Energy Thermal Conversion and Control of Ministry of Education，School of Energy and Environment，Southeast University，Nanjing 210096，Jiangsu，China）

Gas-solid fluidized beds with the presence of some liquid droplets are applied widely in industry. Due to the existence of droplets，particles become wet and viscous，which makes fluidization significantly different with dry particles. In this paper，a hysteresis contact model is applied to consider energy dissipation of lubrication force and liquid bridge force during wet particle collision，which is coupled with CFD-DEM model. The fluidization and mixing characteristics of wet particles under different conditions in a 2D fluidized bed with immersed tubes are compared. Generally，it is found that with the presence of liquid，the particles begin to agglomerate and the bubbles become gas channels. The immersed tubes play the part of frameworks，promoting the formation of gas channels. With increasing cohesiveness，the particle velocity decreases generally and the internal recirculation of wet particles is restrained. It is also found that with an increase in the cohesiveness，the pressure drop of the fluidized bed and the mixing speed of particles decrease.

wet particle；hysteresis contact model；agglomerate；solid mixing

TQ021

A

1000–6613（2017）06–2070–08

10.16085/j.issn.1000-6613.2017.06.015

2016-09-20；

2016-12-08。

国家自然科学基金项目（51306035，51676042）。

刘道银（1982—），男，副研究员，研究方向为多相流及数值模拟。E-mail：dyliu@seu.edu.cn。