第７卷第３期２００７年６月
过程工程学报
ＴｈｅＣｈｉｎｅｓｅＪｏ啪ａ１ｏｆＰｒｏｃｅｓｓＥｎｇｉｎｅｅｒｉｎｇ
、，０１．７Ｎｏ．３
Ｊ１】ｎｅ２００７气固流化床内鼓泡行为的离散颗粒数值模拟
吴春亮１’２，詹杰民２
（１．广东海洋大学工程学院，广东湛江５２４０８８；２．中山大学应用力学与工程系，广东广州５１０２７５）
摘要：采用离散颗粒数值模型探讨了气固流化床内多个气泡的运动规律．数值计算考虑了颗粒与气体的相互作用以及颗粒间的相互作用，颗粒与气体的相互作用采用Ｇｉｄａｓｐｏｗ阻力修正公式，颗粒间的相互作用则采用硬球碰撞模式描述．给出了一种在二维非结构网格下精确计算颗粒体积分数及两相全隐耦合求解颗粒速度的方法，同时，在硬球模型中对颗粒碰撞事件采用散列表的处理方法，大大提高了计算效率．数值计算结果表明，由于多个气泡的相互作用而导致气泡的横向运动显著，气泡运动轨迹与气泡形状依赖于流化床内静压力的分布，气泡总是选择压力小、压差梯度大的路径上升，这一现象与气液分离流中气泡的运动规律类似．
关键词：气固流化床；离散颗粒模型；硬球模式
中图分类号：ＴＱ０１３．２文献标识码：Ａ文章编号：１００９—６０６ｘ（２００７）０３—０４３２—０７
１前言
稠密颗粒两相流在工业上应用非常广泛，如鼓泡式流化床和循环流化床反应装置、循环流化床燃烧炉以及鼓泡塔气液反应器、沉降室和旋流分离器等．另外，部分管道输运中也采用稠密颗粒两相流操作．
大部分的数学模型对连续相采用Ｅｕｌｅｒ方法描述，而基于对颗粒的不同描述方式则分以下两大类：（１）基于Ｅｕｌｅｒ体系的双流体模型；（２）基于Ｌａｇｒａｌｌｇｅ描述体系的离散颗粒模型．另外，一类对颗粒相采用Ｅｕｌｅｒ和Ｌａ铲ａ１１９ｅ组合的模型也得到了发展．有别于稀疏的两相流，稠密的颗粒两相流中颗粒体积含量一般大于１０％，颗粒碰撞的时间尺度小于颗粒驰豫时间尺度，因此数值模型中除了全面考虑流体与颗粒之间的质量、动量和能量的交换，即双向耦合（１’ｗｏ．ｗａｙｃｏｕｐｌｉｎｇ）之外，还需要考虑颗粒与颗粒间的相互作用．
当采用基于Ｌａｇ啪ｇｅ描述方法的离散颗粒模型时，颗粒间的作用模式分为两种：软球模式（ＳｏＲｓｐｈｅｒｅｍｏｄｅｌ）和硬球模式（Ｈａｒｄｓｐｈｅｒｅｍｏｄｅｌ），二者的区别在于颗粒的碰撞方程是采用微分形式还是积分形式．软球模式通过引入弹簧、缓冲器、滑块等概念来描述颗粒之间的接触作用力，允许颗粒之间有小量的重叠，采用刚度系数、阻尼系数和摩擦系数来表达接触位移与颗粒受力的关系．由于采用了微分形式的颗粒碰撞方程，软球模式不但可以得到接触前后颗粒速度的关系，而且还可以计算出接触过程中颗粒的具体受力，正因如此，应用该模式时需要非常小的时间步长（一般在１．０×１０“ｓ，如果颗粒刚度太大则需更小）．在每一个时间步，颗粒的位置
一旦确定，颗粒之间的接触关系就得以确定，于是每个颗粒的受力也确定，这样通过运动方程可以计算出该时刻颗粒的速度以及下一时刻颗粒的位置，该模型的这一特点称之为“时间驱动”．
硬球模式采用积分形式的颗粒运动方程，忽略碰撞的具体过程，通过引入法向恢复系数、切向恢复系数及摩擦系数来描述碰撞前后的颗粒速度关系，假设所有的碰撞是瞬时的二体碰撞，这样在每一次碰撞发生之前更新颗粒位置，同时检测即将发生的碰撞，即颗粒的运动是由一系列的碰撞事件来推动，这一特点称为“事件驱动”．硬球模型中，颗粒位置严格地互不重叠．Ｈｏｏｍａｎｓ掣ｌ】首先实现了基于硬球离散颗粒模型的气固流化床数值模拟，欧阳洁等［２，３】采用离散颗粒模型计算了鼓泡式、循环式流化床内气固流态，数值模型捕捉到了鼓泡、塞状流以及颗粒团聚等现象．Ｂｏｋｋｅｒｓ等［４】采用三维硬球离散颗粒模型，对比了两种气固阻力公式下流化床内单个鼓泡形状的差异，并计算了双组分颗粒的分离度，与实验数据相当吻合．Ｌｕ等［５】结合气相亚网格尺度模型和硬球离散颗粒模式，计算了流化床内颗粒的径向分布函数以及颗粒脉动能，对颗粒动力理论进行了初步验证．刘阳等［６】采用该模型模拟了流化床内不同密度的双组分颗粒分离现象．
在有多个射流口或是自由鼓泡的情况下，如果床层较深，流化床中容易同时形成多个气泡．气泡的运动轨迹、运动方式直接影响到流化床的操作性能．目前，对多个气泡相互作用的数值研究甚少，本工作将采用离散颗粒模型对此展开初步研究．
收稿日期：２００６一ｌＯ一０９，修回日期：２００７—０３—２０
作者简介：吴春亮（１９７７一），男，湖北省咸宁市人，博士，讲师，工程力学专业，Ｔｅｌ：０７５９－２３８３７２０，Ｅ－ｍｉｌ：ｃｈ吼１ｉ粕ｇｖｌｒｌｌ＠ｇｎｌａｉｌ．ｃｏｍ．
第３期吴春亮等：气固流化床内鼓泡行为的离散颗粒数值模拟４３３
２气体一颗粒控制方程
离散颗粒模型对于气相采用Ｅｕｌｅｒ描述方法，其连续性方程为
善（成气）＋Ｖ?（几龟％）＝ｏ，
动量方程为
暑（岛％％）＋Ｖ‘（肪％％％）＝一％即＋Ｖ‘（％元）＋＆＋％肪ｇ．式中，庙为气相密度，龟为空隙率即气相体积含量，乏为剪切应力张量．在Ｌａｇｒａｎｇｅ体系中，单颗粒的运动方程可以根据牛顿第二运动定律表述为
饰导一咋即＋篙（％刊Ｗ．
这里只考虑了两相之间的滞止阻力、重力和浮力，下标Ｇ和Ｐ分别表示气相和颗粒物理量，即Ｇ为流体在颗粒位置的速度，∥为阻力系数，本工作中采用Ｇｉｄａｓｐｏｗ修正阻力系数吐
当旬≤ｏ．８时，∥＝１５０哗＋１．７５（１一龟）争Ｉ％一脚Ｉ；
ｓ＠；、咚、。
当◇ｏ．８时，∥＝三ｃＤ鱼与导风１％一引莳６５，
其中ＣＤ：』若（，删５硝８７）嬲圳ｏｏ，
【ｏ．４４Ｒ邰＞１０００
月如：鱼＆剑旦二型．
∥
其中，绵为颗粒直径，∥为气体动力粘性系数．在气相动量方程中，＆为单位体积内气体和颗粒之间的动量交换，可以对控制体积内的所有颗粒进行积分求得：
＆一吉膳等（％刊鼬喃肌
其中，硒数表明颗粒在积分的控制体内，矿为积分的控制体体积，＾垆为颗粒数目．
颗粒间的作用采用硬球模式描述．该模型假设碰
撞是瞬时的刚性二体碰撞，碰撞的滑动过程中，摩擦力遵从ＣｏｕｌｏＩｎｂ定律．如图１所示，考虑质量和半径分别为聊１，ｒ１以及聊２，圪的２个刚性球１与２的碰撞，速度分别为＇，１，屹，角速度为劬，娩，在碰撞的接触点处，颗粒１相对颗粒２的相对速度为
ｌ，１２＝（Ｖ１一Ｖ２）一（，ｉｑ＋吃％）×一．
定义两颗粒正碰的质心连线为法线方向，与之垂直
的方向为切线方向：
肛商，
Ｉ五一而ｌ
扣静精，
这里上标０表示碰撞前的物理量，Ｇ，：为两颗粒质心的相对速度．设颗粒１受到的冲量为Ｉ，，应用动量定理和动量矩定理，有
％（＇，１一ｌ，？）＝．，，
聊２（１，２一．Ｉ，：）＝一．，，
ｉ（劬一钟）＝Ｉ，×刀，
垒（吃一硝）＝Ｉ，×栉，
其中，为转动惯量．现在需要引入几个参数来关联碰撞前后物理量之间的关系，从而表达由于非完全弹性碰撞引起的能量耗散．首先引入法向恢复系数％，得到法向碰撞冲量为
小吲?珊：（１坞）焉?
再引入摩擦系数体摩擦系数能够表示在接触的滑动过程中摩擦力与正压力的关系，在碰撞过程中，这一关系就是切向冲量和法向冲量之间的关系，服从Ｃｏｕｌｏｍｂ定律：
Ｉ－，×一ｌ＝一肼（．，?以）．
如果切向的相对速度足够大，两颗粒在碰撞过程中一直保持滑移状态（Ｓｌｉｄｉｎｇ）；否则，经过一段滑移之后，两颗粒保持粘附状态（Ｓｔｉｃｋｉｌｌｇ）．对于这种情况，类似法向恢复系数，引入切向恢复系数白来表示碰撞前后两颗粒切向相对速度的关系：
图１两个颗粒的碰撞示意图
Ｆｉｇ．１１１ｌｕｓ信ａｔｉｏｎｏｆｃｏｌｌｉｓｉｏｎｂｅ锕ｅｅｎｔ、）ｌｒｏｐａｎｉｃｌｅｓ对于上述两种情况，分别求得切向的冲量为
４３４
过程工程学报第７卷
当肼＜至等笔胖时，‘舢ｍ。：一所以；
～
７（巩＋埘，）＾
。…““
～“
当舻群时，‘蛐刊Ⅷ孺．
３数值方法简述
本工作采用非结构化网格下的有限体积法离散气相输运方程，单相流的ＳＩＭＰＬＥ算法用来求解气相速度与压力的耦合［８，９】．离散后的气相速度方程、压力校正方程采用代数多重网格法进行求解．由于离散颗粒模型采用Ｌａ铲ａｎｇｅ方法追踪具体的每一个颗粒运动，于是网格单元上的空隙率可以根据颗粒占据该网格的份额来计算．一般地，颗粒尺寸要远远小于Ｅｕｌｅｒ网格尺寸，但是从我们的数值经验来看，仍需精确计算颗粒的体积份额．如果简单将颗粒体积累加到颗粒中心位置所在的网格，则导致气相连续性方程计算难以收敛．例如，图２（ａ）所示正方形网格计算到单元１的颗粒体积分数约为０．４，图２（ｂ）则约为０．２，图２（ａ）的计算会带来５０％的相对误差，并且有可能使某个单元的颗粒体积含量超过紧致含量．因此当颗粒边界跨越多个网格时，必须准确计算落入各个网格的颗粒面积．
图２粗略计算颗粒体积分数可能带来５０％的相对误差
Ｆｉｇ．２
５０％ｅｒｒｏｒｓｃａｕｓｅｄｂｙｏｍｉｔｔｉｎｇｐａｆｔｉｃｌｅｓｈａｐｅ
文献［１０］给出了一种在结构网格下准确计算空隙率的方法．在非结构网格下，如图３所示，一个颗粒与网格之间的关系只可能是图中所示的３种之一，方便起见，记Ｃｆ为颗粒ｐ中心所在的网格，记Ｇ（胛＝１，２，…）为颗粒边界跨越的网格，所有的坐标位置都转换为以颗粒中心为原点的局部坐标．颗粒对其跨越的网格的面积贡献
分别记为彳ｊ，么。（肛１，２，…）．
图３颗粒与网格的３种位置关系
Ｆｉｇ．３Ｔｈｒｅｅｋｉｎｄｓｏｆｐａｎｉｃｌｅｐｏｓｉｔｉｏｎｓ
ｏｎ
１１ｎｓ饥ｌｃｔｕｒｅｄ鲥ｄ
对于情况１，４＝４＝石带，４，＝ｏ：对于情况２，彳。＝
口＂。。卅’＂。。）一口州—爿’。，曲，４＝万名一∑４，彳＂。。，彳’＂。分
别为２条边心和ＰＱ对应的较大的扇形面积和较小的
扇形面积，而彳枷，爿’删则为对应的２个三角形＾幻Ⅳ，
ＰＤＱ的面积［如图４（ａ），图３的Ｃａｓｅ２是特例，即么’＂。。＝０，彳’。．们＝０］．
对于情况３，４＝刀带一泓，４＝４触＋４—４，此时，
颗粒与某个单元ｅ的２个面（边）有交点，其中颗粒对单元ｅ的三角形面积爿＋的正负贡献根据下列原则来确定：如图４①）所示，如果颗粒的径线（指颗粒中心与交点的连线，如０ＩⅣ或Ｄ＾力与该单元的另外一条边（指不与
交点对应的边，如ａＩⅣ对应为Ｍ矿或Ｇ啊对应为Ⅳｎ相
交，则该径线与圆心组成的三角形对该单元是正的面积贡献，否则为负的贡献．
在计算实现上，对每个颗粒中心所在的网格单元执行一遍界面（边）循环，以区分不同的颗粒位置情况．虽然计算过程稍微复杂，但是该方法计算出的颗粒体积分数是充分准确的，且适合于二维情况下所有的结构化网格、非结构化网格、结构化和非结构化混合网格．由文献［１，８］描述的ＤＰＭ计算流程可以看出，在每个时间步，颗粒体积分数只计算１次，因此稍微付出的计算量可以接受．从计算实验来看，采用上述方法计算颗粒体积分数的时间与计算颗粒碰撞或求解连续相离散方程组的时间相比小得多．另外，颗粒占重叠网格的体积份额可以用来加权计算颗粒对气相的阻力作用，这样处理对源
项的计算起到了光滑化的作用，有利于气相方程组的迭代收敛．
（ａ）
（ｂ）
图４两种不同颗粒位置下颗粒体积分数的计算
Ｆｉｇ．４Ｃａｌｃｕｌａｔｉｏｎｏｆｐａｎｉｃｌｅｖｏｌｕｌｎｅ丘ａｃｔｉｏｎｗｈｅｎ
ａ
ｐａｎｉｃｌｅ
ｏｖｅｒｌ印ｓ埘ｔｈｔ１１ｅＥｕｌｅｒｉａｌｌｃｅｌｌ
硬球离散颗粒模型采用颗粒运动一碰撞解耦的方法，分两步求解颗粒的运动【８】．第一步计算两相作用，颗粒固定在空间位置，直接求解颗粒动量方程，此时考虑颗粒所受到的所有非碰撞外力，如重力、气相阻力、
升力等．求解颗粒运动的第二步主要是计算颗粒间的作
用．硬球模式采用碰撞动力学描述颗粒问的作用，一旦