, 姜霖松1, 吴丹2, 解茂昭11. 大连理工大学 能源与动力学院, 辽宁 大连 116024;
2. 中国船舶重工集团公司第760研究所, 辽宁 大连 116013
收稿日期：2017-03-31
基金项目：国家自然科学基金资助项目(51576029，51606176)；中央高校基本科研业务费专项资金资助项目(DUT16QY19)。
作者简介：刘宏升(1978-)，男，辽宁大连人，大连理工大学副教授。
摘要：采用离散元软件LIGGGHTS重现球形堆积床的重力堆积过程, 建立了球形堆积床的三维随机结构模型, 同时构建了顺排结构和插排结构两种有序堆积模型与之对比.采用CFD软件Fluent14.0对三种结构堆积床内的湍流流动进行了数值模拟研究.通过宏观阻力系数和压强降的分析, 验证了随机模型的有效性, 并讨论了球形堆积床随机结构内速度场特性和湍流参数的变化规律.结果表明:与经典的Ergun公式相比, 随机模型的宏观阻力系数和压强降的整体误差低于6%;与有序堆积床相比, 随机堆积模型具有流体速度场均匀性好、尾旋尺寸大、扰动强烈及瞬时耗散率极值高等特点.
关键词：多孔介质堆积床随机结构速度场湍流特性
Simulation on Turbulence Characteristics in Three-Dimensional Sphere Packed Bed
LIU Hong-sheng1
, JIANG Lin-song1, WU Dan2, XIE Mao-zhao11. School of Energy and Power Engineering, Dalian University of Technology, Dalian 116024, China;
2. The 760th Research Institute of China Ship Building Industry Corporation, Dalian 116013, China
Corresponding author: LIU Hong-sheng, E-mail: lhsh@dlut.edu.cn
Abstract: Discrete element software LIGGGHTS is introduced to reproduce the forming process of sphere packed bed by gravity, and a 3D random geometric model similar to the actual accumulation structure is presented. Two ordered models, the in-line packed model and the staggered packed model, are also established to contrast with the random packed model. CFD software of Fluent 14.0 is employed to study the turbulence flows in the packed bed. The effectiveness of the random model is validated by analysis of the macroscopic pressure drop and the resistance coefficient, and the velocity field and the turbulence parameters of the three structures are predicted and discussed. Simulation results show that the macroscopic resistance coefficient and the macroscopic pressure drop of the random packed bed model are consistent with the classical formula of Ergun, with maximum relative error less than 6%. Compared with the orderly packed bed, the random packed model is characterized with more homogeneous velocity field, larger tail vortex size, stronger disturbance and extremely larger turbulence dissipation.
Key words: porous mediumpacked bedrandom structurevelocity fieldturbulence characteristic
多孔介质内的燃烧又称为过滤燃烧, 是指采用多孔介质材料取代自由空间, 使可燃混合气或燃油蒸汽经固体多孔介质或颗粒堆积床, 在类似于过滤情况下发生的燃烧, 具有燃烧稳定、传热特性好、燃烧效率高、排放水平低等优势[1].
近二十多年来, 多孔介质燃烧技术逐渐在国际燃烧界形成研究热点[2], 相关研究工作大部分是针对气体燃料预混燃烧机理进行, 研究内容包括燃烧特性与热效率[3]、火焰稳定性[4]、贫燃极限[5]等方面.
多孔介质是一种具有多尺度结构特性的复杂几何系统, 其内部孔隙的形状、大小及其空间分布通常都是随机的.建立具有实际多孔介质几何特征的三维随机模型, 是深入理解过滤燃烧机理的重要前提.迄今为止, 各国研究者采用的多孔介质结构重构方法[6-8]包括:物理方法、数值方法、分形理论等.由于多孔介质复杂的结构性, 生成的网格数量巨大, 对计算要求较高.传统的理论与数值模拟研究一般采用宏观描述的简化方法, 通过对基本方程进行宏观体积平均并选取相应经验公式[9], 或将三维复杂结构简化为二维结构[10-12], 建立具有周期性的阵列方块、圆球或椭球结构; 这样简化了问题, 但也与实际情况有很大差别.
本文采用离散元软件LIGGGHTS重现小球重力堆积过程, 充分考虑堆积过程中重力、球间碰撞及壁面碰撞引起的轨迹变化, 得到与实际堆积结构相近的三维球堆积床随机结构模型, 通过分析随机结构堆积床内的宏观阻力系数与压强降变化规律, 确定模型的有效性.同时建立小球尺寸相同的顺排、插排有序堆积模型, 与随机结构对比, 分析三种结构堆积床内流体流动的速度场、湍动能及耗散率等参数的发展规律, 旨在为多孔介质内过滤燃烧的深入研究奠定基础.
1 物理模型与计算方法本文建立的球堆积结构包括对称结构与随机结构, 其中对称结构又按照顺排和插排两种方式构建.在小球顺排单元中(见图 1a), 8个小球分别置于正方体的8个顶点, 相连两球之间球心距相等; 将多个单元有序排列, 可得到顺排型小球堆积床.在插排单元中(见图 1b), 5个小球分别置于正四棱锥的5个顶点处, 相连两球之间球心距相等; 将单元阵列按规律有序排列后, 得到插排型小球堆积床, 其中插排堆积床侧壁处的小球设为半球结构.
图 1(Fig. 1)
 | 图 1 三种堆积方式的单元体和堆积体Fig.1 Cell body and packed bed of different structures (a)—顺排单元；(b)—插排单元；(c)—随机堆积. |
小球随机排列结构采用国际开源离散元软件LIGGGHTS, 通过重新编程设定重建小球的重力堆积过程, 得到与实际堆积结构更为接近的球堆积床三维随机模型.在仿真堆积过程中, 首先设置小球的堆积参数, 在小球堆积区域上方设定投掷高度与投掷区域面积, 在区域内小球可以从任一随机位置开始自由落体运动, 如图 1c所示; 然后建立与小球的材料硬度、表面粗糙度相关的碰撞模型, 通过软件求解随机统计方程, 进而确定小球碰撞过程的运动轨迹及落点位置, 模拟小球的碰撞、堆积过程, 最终得到随机成形堆积体, 实现对小球堆积床真实的重力堆积过程的模拟.
2 数学模型根据雷诺数判断多孔介质区内的流动为湍流, 故采用标准雷诺数k－ε湍流模型; 对于气体流动问题可忽略重力的影响.具体控制方程如下[13].
连续性方程:
 | (1) |
 | (2) |
 | (3) |
 | (4) |
边界条件:
入口z=0: u=v=0, w=w0;
出口
壁面条件:u=v=w=k=ε=0.
采用通用CFD软件Fluent14.0求解上述方程.以空气作为计算工质, 因不涉及温度变化, 物性参数均认为是常数, 取ρ=1.205 kg/m3, μ=1.818 9×10-5 kg/(m·s).利用GAMBIT软件, 进行非结构化网格划分, 并对网格进行无关性验证, 得到局部加密后的网格.最终得到顺排、插排和随机堆积情况下网格数分别为1 583 960, 1 143 933和1 755 443, 其中随机堆积结构的小球位置数据由LIGGGHTS软件的计算结果导出.由于随机结构中小球分布没有规律, 为了达到无关性标准, 在小球交界面附近的局部空间内需要对网格进一步加密, 因此随机堆积结构网格数量最多.
在三种小球堆积床结构中, 小球直径为5 mm, 顺排、插排和随机堆积结构对应的孔隙率分别是0.472, 0.290和0.389.由于小球间距很小, 流动截面剧烈变化引起的“缩-扩”现象并不明显, 大量计算表明, “缩-扩”现象并没有造成明显的压力梯度突变, 因此在进行压力计算时采用SIMPLE格式, 而动量、湍动能及耗散率均采用具有三阶精度的QUICK格式.经过验证, 残差收敛判断标准为10-6时, 其对计算结果的影响可以忽略不计.
3 模型验证为了验证随机结构模型在流动模拟方面的有效性, 以文献[14]中球堆积床多孔介质的几何参数作为建模和流动计算的参考数据, 建立了随机结构的多孔介质堆积区, 将其填充至界面边长为17.6 mm, 总长为450 mm的方形管道内, 其中多孔介质堆积区长度为150 mm, 布置在管道的中心位置, 入口段和出口段长度均为150 mm, 这样可以避免入口效应及出口效应对多孔介质区域内流动模拟的影响.
在5种不同雷诺数
图 2(Fig. 2)
 | 图 2 宏观阻力系数和压强降随雷诺数Re的变化Fig.2 Variations of the macroscopic pressure drops and resistance coefficient with Re |
由图 2可以看出, 随着雷诺数的增加, 压强降整体呈增加趋势; 低雷诺数(103以下)时, 计算值与Ergun公式相对误差在3%以下, 高雷诺数时压强降的相对误差保持在6%左右.阻力系数f随雷诺数的增加呈下降趋势, 低雷诺数时的计算值与Ergun公式基本重合; 当雷诺数超过8×104后, 相对误差范围稳定在5%~7%之间, 表明球堆积床随机结构模型在预测多孔介质压强降和摩擦阻力系数方面是有效的.
大量堆积结果表明, 由于小球的下落及碰撞过程具有很强的随机性, 即使在相同的堆积参数条件下, LIGGGHTS软件建立的随机模型在其局部微观结构上仍存在较大差异; 但对大量堆积数据进行统计后发现, 在小球尺寸相同、堆积区域相同的条件下, 当堆积的小球数量达到一定值以后, 随机性对由积分平均法得到的宏观阻力系数和宏观压强降并没有影响.
为验证随机结构模型在孔隙率预测上的有效性, 本文按照文献[15]中的实验工况, 计算了随机多孔介质结构的径向孔隙率[15], 将多孔介质区域沿径向(x方向)均分为多个条形区域, 通过积分计算各剖面上的孔隙面积, 并计算出各剖面的局部孔隙率, 得到不同x位置处的孔隙率值.从图 3中可以看出, 计算值与实验结果保持了相似的变化规律, 整个计算区域内最大相对误差为9.1%, 出现在x为24 mm处; 该位置是靠近壁面处的第一个孔隙率波谷位置.出现偏差的主要原因是, 文献[15]实验结果中的孔隙率是对三维模型取环形区域进行面积分得到的, 而本文仅取三维随机模型中的一个二维切面进行计算, 孔隙率存在一定的不均匀性.此外模拟结果中相邻波峰或波谷之间的距离均小于实验值; 考虑到文献[15]中指出该间距与堆积球体直径是基本相当的, 而本文模型二维切面中的圆面直径都是小于或等于三维堆积模型中的球体直径, 因此模拟值低于实验值也是合理的.
图 3(Fig. 3)
 | 图 3 局部孔隙率模拟值与实验值对比Fig.3 Local porosity contrast between simulation and experiment |
4 计算结果与分析图 4给出了在入口速度10 m/s条件下, 三种不同堆积结构内的流线分布图.在顺排堆积结构(图 4a)内, 各横截面上小球彼此平行, 小球之间的空隙形成明显的z向直通道, 该通道内的气流速度明显高于流域内的其他位置; 该结构下气流瞬时最大速度为62.502 m/s.在堆积区尾流附近处受到多孔介质扰动而形成局部的尾流涡旋, 但旋涡尺度很小, 对尾流区流场的影响并不明显, 尾流区的流线整体上呈平行趋势.
图 4(Fig. 4)
 | 图 4 多孔介质结构内的流线图Fig.4 Streamline in packed bed structures (a)—顺排结构；(b)—插排结构；(c)—随机结构. |
在插排堆积区(图 4b)中, z向直通道被打断, 流场与小球交替出现, 流线也随之密集、稀疏地周期性变化, 在流动过程中流体不断地掺混并分离, 导致速度方向不断变化; 该结构下气流瞬时最大速度为201.945 m/s.插排结构尾流区内流线受到的扰动要强于顺排结构, 靠近流道壁面附近出现少量旋涡, 但流线整体上仍接近直线分布, 这是因为插排分布较顺排分布更为紧密, 将空间分割成更小的区域则孔隙尺度更小, 湍流强度随之减小.
在随机堆积区(图 4c)中, 流线呈随机且无序状态, 气流瞬时最大速度达到132.038 m/s, 介于顺排、插排之间.与对称结构不同的是, 在随机堆积区尾流处出现了非常明显的旋涡, 尺寸在15 mm左右, 旋涡几乎充满了整个流道, 这是由于随机堆积区域的尾端处小球排列同样是随机的, 形成的流道方向各异, 并且会有凸起的小球对流体造成扰动; 因此随机堆积床对脱离堆积区域流体的扰动最为强烈, 表明随机结构更容易在尾流区形成湍流, 从而方便不同流体的相互掺混.
由流线局部放大图可以看出, 在顺排结构(图 5a)中, 小球接触面的迎流面流线密集, 背流面上速度矢量则较稀疏且速度方向杂乱, 形成了明显的涡旋, 这与速度矢量图是一致的.这除了是因为球面绕流产生的扰动外, 还由于流体绕过接触面后流动空间突然变大, 导致雷诺数骤然变化, 从而为涡旋形成提供了有利条件.在插排结构(图 5b)中, 小球背流面虽然速度方向有所改变, 但并没有形成涡旋, 这主要是因为插排堆积结构相对紧密, 分割后的区域孔隙尺度过小, 并不足以产生涡旋.
图 5(Fig. 5)
 | 图 5 多孔介质结构内流线局部放大图Fig.5 Local larger image of streamline in packed bed structure (a)—顺排结构；(b)—插排结构；(c)—随机结构. |
随机结构内的流线分布存在两种情况, 在图 5c左图中的小球表面没有涡旋, 而右图中小球的背流面产生了明显涡旋, 这在一定程度上说明随机堆积床内流线的随机性与无序性, 而且随机结构多孔介质内的速度分布最为均匀, 反映了随机多孔介质内均质流场的特性, 这必将成为多孔介质内均质燃烧的有利条件.
图 6给出了沿流动方向(z向)不同截面上速度积分平均值的分布规律.可以看出, 由于结构的对称性, 插排与顺排多孔结构中的速度分布呈现明显的周期性.这里为了避免进出口效应的影响, z值范围取在20~120 mm内.在顺排堆积结构中, z向直通道区域内的平均速度明显大于其他位置, 从而使速度出现周期性波动, 每个周期跨域的尺度为6.7 mm左右, 振幅中心处平均速度为12.954 m/s.在插排结构中, 虽然未形成明显的z向直通道, 但该结构的孔隙率最低, 流道上孔隙与小球交替出现使流体不断地掺混并分离, 致使速度方向不断变化, 因此插排结构横截面上的平均速度最高, 振幅中心处平均速度为52.448 15 m/s, 而且其速度波动振幅最大.在随机结构下, 速度分布相对稳定, 不再具有明显的周期性, 振幅中心处平均速度在13.771 764 m/s, 介于顺排、插排结构之间.
图 6(Fig. 6)
 | 图 6 流动方向上截面平均速度分布Fig.6 Distribution of average velocity in flow direction on the cross section |
图 7给出了入口速度10 m/s条件下, 三种堆积床结构中心截面湍动能分布图.可以看出, 顺排和插排堆积区域内湍动能呈有序化分布, 高湍动能区出现在y向(垂直于主流方向)接触的两小球接触点前后位置, 低湍动能区出现在每个小球接触面背流面.对比观察图 7b, 可以看出, 高湍动能区出现在三小球挤压形成的喉道处, 之后因流体在接触面周围充分掺混, 湍动能水平降低, 在接触面附近的背流区则出现低湍动能区; 随着流体进入喉道, 湍动能再次升高, 如此往复.由图 7c可以看出, 随机结构湍动能分布与对称结构存在很大差异, 堆积区域内湍动能为无序分布.
图 7(Fig. 7)
 | 图 7 不同结构中心截面湍动能分布Fig.7 Turbulent kinetic energy distribution on central section in different structures (a)—顺排结构；(b)—插排结构；(c)—随机结构. |
图 8给出了在流动方向(z向)上截面湍动能积分平均值的分布规律.可以看出, 在顺排结构中湍动能分布呈现明显的周期性, 湍动能数值在38.56~49.77 m2/s2范围波动, 局部瞬时湍动能的最大值为111.257 9 m2/s2.在插排结构中, 湍动能在129.27~170.34 m2/s2范围波动, 虽然湍动能出现了增减交替的变化规律, 而且波动幅度比顺排更明显, 但这种变化的周期性并不明显, 这表明插排的交错结构使流体不断地掺混, 速度梯度很大, 从而破坏了湍动能周期性变化规律; 其局部湍动能的最大值达到440.476 7 m2/s2, 远高于顺排结构.在随机结构中湍动能呈无序性分布, 最高湍动能为230.487 2 m2/s2, 介于顺排与插排结构之间; 沿流动方向整体上看, 湍动能略有增大的趋势, 平均湍动能数值从11.78 m2/s2增至39.72 m2/s2, 这是由于随机结构对流动具有强烈的扰动效果, 增强了局部孔隙内的速度脉动, 从而使湍动能有所增加.
图 8(Fig. 8)
 | 图 8 流动方向截面平均湍动能分布Fig.8 Distribution of average turbulent kinetic energy in flow direction on the cross section |
图 9给出了在流动方向(z向)上截面湍流耗散率的积分平均值分布规律.可以看出, 在顺排结构中湍流耗散率呈现明显的周期性, 湍流耗散率在3.324×105 ~7.582×105 m2/s3范围波动, 振幅中心处湍流耗散率为5.741×106 m2/s3, 湍流耗散率局部最大瞬时值为7.410 99×106 m2/s3; 在插排结构中湍流耗散率在3.922×106~6.576×106 m2/s3范围波动, 变化幅度最大, 局部最大湍流耗散率为4.950 323×107 m2/s3; 在随机结构中湍流耗散率最低且呈无序性分布, 湍流耗散率从9.826×104 m2/s3增至5.469×105 m2/s3, 整体上要低于顺排与插排结构, 但是其最高耗散率高达7.739 875×108 m2/s3, 远远高于两种对称结构, 这是因为随机堆积的排列空间扭曲度较大, 造成局部区域耗散率值极高的现象, 因此在进行随机结构耗散率研究时, 有必要对局部区域进行特殊分析, 以考虑局部极端现象对燃烧过程造成的影响.
图 9(Fig. 9)
 | 图 9 多孔介质结构的湍流耗散率分布Fig.9 Distribution of average turbulent dissipation rate in flow direction on the cross section |
5 结论1) 三维随机结构模型的宏观阻力系数和宏观压强降的计算值与经典Ergun吻合较好, 相对误差低于6%, 局部孔隙率计算结果与实验值最大误差在9%以内, 证明随机结构模型可以用于湍流流场的流动分析.
2) 顺排结构中的速度、湍动能与湍流耗散率均呈周期性规律分布, 周期宽度在6.7 mm左右, 平均速度在12.954 m/s附近波动.
3) 插排结构孔隙率最低, 插排结构内的速度、湍动能与湍流耗散率值均为最大, 平均速度在52.448 15 m/s附近波动.
4) 随机堆积床内速度分布相对均匀, 平均速度介于顺排、插排之间, 但随机堆积结构对尾流区的扰动最强烈, 可形成尺寸约15 mm的涡旋.
参考文献
| [1] | Mujeebu M A, Abdullah M Z, Bakar M Z A, et al. Combustion in porous media and its applications:a comprehensive survey[J].Journal of Environmental Management, 2009, 90(8): 2287–2312.DOI:10.1016/j.jenvman.2008.10.009 |
| [2] | Mujeebu M A, Abdullah M Z, Mohamad A A, et al. Trends in modeling of porous media combustion[J].Progress in Energy and Combustion Science, 2010, 36(6): 627–650.DOI:10.1016/j.pecs.2010.02.002 |
| [3] | Yu B H, Kum S M, Lee C E. Combustion characteristics and thermal efficiency for premixed porous-media types of burners[J].Energy, 2013, 53(1): 343–350. |
| [4] | Gao H B, Qu Z G, Feng X B, et al. Methane/air premixed combustion in a two-layer porous burner with different foam materials[J].Fuel, 2014, 115: 154–161.DOI:10.1016/j.fuel.2013.06.023 |
| [5] | Stelzner B, Keramiotis C, Voss S, et al. Analysis of the flame structure for lean methane-air combustion in porous inert media by resolving the hydroxyl radical[J].Proceedings of the Combustion Institute, 2015, 35(3): 3381–3388.DOI:10.1016/j.proci.2014.06.151 |
| [6] | Song Y C, Jiang L L, Liu Y, et al. An experimental study on CO2/water displacement in porous media using high-resolution magnetic resonance imaging[J].International Journal of Greenhouse Gas Control, 2012, 10: 501–509.DOI:10.1016/j.ijggc.2012.07.017 |
| [7] | Roozbahani M M, Huat B B K, Asadi A. The effect of different random number distributions on the porosity of spherical particles[J].Advanced Power Technology, 2013, 24(1): 26–35.DOI:10.1016/j.apt.2012.01.006 |
| [8] | Miao T J, Yu B M, Duan Y G, et al. A fractal model for spherical seepage in porous media[J].International Communications in Heat and Mass Transfer, 2014, 58: 71–78.DOI:10.1016/j.icheatmasstransfer.2014.08.023 |
| [9] | Lutsenko N A. Modeling of heterogeneous combustion in porous media under free convection[J].Proceedings of the Combustion Institute, 2013, 34(2): 2289–2294.DOI:10.1016/j.proci.2012.06.147 |
| [10] | Nakayama A, Kuwahara F. A macroscopic turbulence model for flow in a porous medium[J].Journal of Fluids Engineering, 1999, 121(2): 427–433.DOI:10.1115/1.2822227 |
| [11] | Kuwahara F, Kameyama Y, Yamashita S, et al. Numerical modeling of turbulent flow in porous media using a spatially periodic array[J].Journal of Porous Media, 1998, 1(1): 47–55.DOI:10.1615/JPorMedia.v1.i1 |
| [12] | Pedras M H J, de Lemos M J S. Computation of turbulent flow in porous media using a low-Reynolds k-model and an infinite array of transversally displaced elliptic rods[J].Numerical Heat Transfer, 2003, 43(6): 585–602.DOI:10.1080/10407780307349 |
| [13] | 陶文铨. 数值传热学[M]. 西安: 西安交通大学出版社, 2001. ( Tao Wen-quan. Numerical heat transfer[M]. Xi'an: Xi'an Jiaotong University Press, 2001.) |
| [14] | Ergun S. Fluid flow through packed columns[J].Chemical Engineering Progress, 1952, 48(2): 89–94. |
| [15] | Mueller G E. Radial void fraction distributions in randomly packed fixed beds of uniformly sized spheres in cylindrical containers[J].Powder Technology, 1992, 72(3): 268–275. |
