含垂直列管束内构件鼓泡塔的CFD模拟

李兆奇，赵远方，管小平，王丽军，成有为，李希



含垂直列管束内构件鼓泡塔的CFD模拟

李兆奇，赵远方，管小平，王丽军，成有为，李希

（浙江大学化学工程与生物工程学系，浙江杭州310027）

针对含有密集垂直列管束内构件的气液鼓泡塔，在两相Euler二维轴对称-模型中，分别考虑气相和液相受到的列管阻力。通过引入相应的动量源、湍动源以及耗散源建立带列管内构件的鼓泡塔二维CFD模型。模型能清晰、准确地描述带列管束鼓泡塔中气液流动的特征：“烟囱效应”以及分布器影响区延长。计算得到的气含率以及液速的二维分布在宽泛的表观气速（0.12～0.62 m·s-1）范围内与实验值相符。

鼓泡塔；密集列管束；内构件；计算流体力学；模型

引 言

鼓泡塔反应器是一种工业中广泛应用的气液两相和气液固三相接触与反应装置，对于强放热反应如费托合成，工业上通常安装密集的换热列管来移除反应热。文献中关于列管内构件对鼓泡塔流动影响的实验研究已有很多[1-5]。但关于含列管内构件鼓泡塔的CFD模拟报道较少。Bernemman[1]利用一维动量方程，给定抛物形气含率和湍流黏度分布，将列管束的影响视作均匀施加的摩擦阻力考虑，引入压力梯度项，模拟得到的液速分布与实验结果符合较好。Forret等[4]认为带内构件的鼓泡塔不适合用一维轴向扩散模型（ADM），通过引入轴向液速分布、轴/径两向扩散系数，建立了2D模型来预测内构件对液相混合的影响。Larachi等[6]采用三维CFD模拟考察直径1 m的鼓泡塔内，4种不同列管布置的流场分布，结果表明，加入列管后，湍流动能降低至空塔的1/3左右。张煜等[7]将列管束对流动的影响平均化为单位体积的动量源、湍动源和耗散源，摒弃了将列管壁面考虑为几何边界的方法，使得计算量明显降低，且模拟得到的充分发展段气含率和液速分布与实验值符合较好，能够较好地解释列管束产生的“烟囱效应”。

Bernamann[1]、Forret等[4,8]、Chen等[9]、张煜 等[3]曾测量了直径0.44～1.0 m列管鼓泡塔内的液速和气含率分布，结果却不一致，对于“烟囱效应”是否存在仍有争议，赵远方等[10]在500 mm×5000 mm的鼓泡塔安装40根垂直列管束并详细测定了不同轴向高度的液体速度和气含率分布，其实验研究表明，与空塔相比，列管塔中分布器影响区的范围显著增大，当流动发展至约4倍塔径以上才能产生“烟囱效应”。含列管束内构件鼓泡塔的流体力学模型必须反映这一客观事实，而张煜等[7]的一维模型仅能描述充分发展区的“烟囱”现象，对于列管塔中不可忽略的分布器影响区的流动情况缺乏考虑，Larachi等[6]和Laborde-Boutet等[11]将列管外壁考虑为几何边界进行全塔三维模拟，但其网格密度很大，计算量相当庞大，且计算结果与文献中相关实验报道相差很大。因此，有必要发展合适的流体力学模型，从列管内构件的力学机理出发，准确高效地模拟列管塔的整体流动规律。本文借鉴张煜等[7]将列管束的影响平均化为相应源项的方式，在文献[12]基于径向力平衡的二维轴对称空塔模型的基础上，全面地考虑列管束对气泡相以及液相的力学作用，在Euler-模型中引入相应的动量源、湍动源以及耗散源，模拟列管束存在条件下鼓泡塔内流动的发展规律，在阐释“烟囱效应”的同时揭示列管束造成分布器影响区增大的流体力学机理。

1 数学模型

文献[3, 10]中的实验研究表明，与空塔相比，列管塔中流型主要表现为两个特征：一是“烟囱效应”；二是分布器影响区延长。文献[12]曾针对无列管条件下的气液鼓泡塔建立了基于径向力平衡的两相Euler二维轴对称-模型，模型计算稳定性较好，易于收敛，本文在该模型的基础上综合考虑列管束对流动所产生上述两方面影响的力学机制，其一表现为列管束对液相的作用，主要是抑制液相湍动，湍流黏度下降导致充分发展段的液速分布比空塔更为陡峭，本文将文献[7]的一维模型中列管束产生的动量源、湍动源及耗散源拓展至二维轴对称模型；其二表现为密集列管束限制气泡向塔中心聚集的径向运动，迫使气泡群趋向于沿直线式浮升，导致流动发展变缓，而以往文献中关于列管束对气泡作用的流体力学研究很少，本文将参照气泡群壁面润滑力的概念，首次提出密集列管束对气泡群作用的力学源项。通过综合考虑列管束对液相以及气泡相的力学作用，建立一个较为完善的列管塔流体力学模型来考察带列管条件下鼓泡塔中的流动情况。

1.1 列管束源项的流体力学机制

在流体力学模型中将列管束考虑成源项基于以下两点假设：①鼓泡塔内布置的竖直列管足够密集，从而列管对流动的阻碍和耗散作用可以平均化表示成单位体积的动量源、湍动源和耗散源；②忽略列管壁面附近的局部分布，只考虑远离壁面的鼓泡塔内速度和气含率的“宏观”分布。列管对鼓泡塔气液两相流的力学作用主要表现为两个方面：一是对液相平均速度以及脉动速度的轴向摩擦和径向形体阻力，文献[7]的一维模型中阐述了列管束对液相脉动的阻碍是产生“烟囱”效应的主要原因；另一方面，列管对气泡径向运动的阻力是导致流动发展变慢，分布器影响区延长的原因。基于上述考虑，本文将列管内构件作用归纳如下。

1.1.1 动量源 列管对液体的轴向摩擦阻力及径向形体阻力如文献[7]中使用的表达式

式中，S为摩擦阻力系数，与列管塔Reynolds数有关，在本文所考察的湍动鼓泡塔操作范围内，S0.005；S为单位体积流体中列管的比表面积，m-1；K是von kármán系数，K的取值可由流体横向通过圆柱的绕流阻力系数公式[13]估算，在本文考察的流速范围内，可取近似值，K1.7；L是时均速度的径向分量，m·s-1；F是单位体积流体内列管迎风面积，m-1。

文献中关于列管对气泡群作用力的研究较少，而关于塔壁边界对气泡的作用，Antal等[14]及Tomiyama[15]曾经提出壁面润滑力的概念，在两相管道流中，气泡受到管壁附近强烈剪切的影响，会受到一个垂直于流动方向，驱使气泡远离壁面的作用力。赵远方等[10]的实验研究表明，列管内构件的存在限制了气泡的径向和切向运动，本文参考壁面润滑力的表达形式[16]，提出列管壁对气泡群的力学作用源项，表达式如下

式中，列管对气泡群运动阻力与气相运动速度和列管密度呈正比，IG为阻力系数，取值为0.002；G为气相径向速度，m·s-1；S为气液滑移速度，m·s-1；为列管所占截面积比；B为气泡直径，m；BBSL，为气泡Reynolds数。

1.1.2 湍流动能和耗散率源项 选用标准-模型来封闭动量方程，分别在和方程中增加列管束产生的湍动源和耗散源项。鼓泡塔充分发展段内竖直列管束对流体湍动主要有两方面作用：一是列管对液体径向脉动的形状阻力，导致湍流动能降低，用体积源kr（m2·s-3）来表示；二是列管束对流体轴向流动产生的摩擦阻力，作用在管壁附近形成湍流边界层剪切区，把平均流动转化为湍动，使得湍流动能上升，用体积源sz（m2·s-3）表示[7]。

式中，c为模型参数，取值为3.0。

体积源sz为列管在近壁面处流体剪切产成的湍动动能，由流动阻力关系可知其与轴向平均速度3次方以及列管比表面积呈正比，表述如下

1.2 守恒方程与本构关系

在Euler-Euler双流体模型中考虑列管束产生的动量源、湍动源和耗散源，基本运动控制方程如下：

连续性方程

动量方程

式中，DLTD分别代表曳力、升力和湍流扩散力，N·m-3；internals表示列管束对液相的摩擦阻力kf,L、形体阻力kd,L以及对气泡群的阻力IG，N·m-3。其中，曳力公式采用Schiller & Naumann曳力模型

气泡假设为均匀气泡，气泡直径（B）通过轴向浮力和曳力相互平衡确定

除曳力外，气液相间考虑升力[17-18]与湍动扩散力[19]的作用

混合相的湍流动能、湍流耗散以及湍流黏度由标准-模型推导得到

式中，sz与kr分别表示由于列管束存在导致湍流动能的增加与减少，其中参数C1、C2Cσσ的取值分别为1.44、1.92、0.09、1.0、1.3，式(13)最后一项表示与列管产生的湍动源项相对应的耗散相，参数C3通过实验确定，C31.4。

1.3 模型对象及求解方法

1.3.1 模拟对象 本文模拟对象为均匀布置垂直列管束的500 mm（476 mm id）气液鼓泡塔，其中列管直径为25 mm，管长4 m，管间距分别取60 mm（列管数40）、75 mm（31）与90 mm（18），分别考察了0.12、0.31、0.62 m·s-13个表观气速下列管塔（40，31，18）与空塔（0）中的流动情况，采用二维轴对称模型。网格在Gambit 2.2.30里生成，导入商业软件Fluent 6.3求解。网格采用结构化网格划分，并对边壁边界层格网加密，共计30×800个网格，如图1所示。

图1 网格示意图

1.3.2 数值解法 模型通过非稳态求解，离散格式时间导数采用一阶隐式格式，空间导数采用一阶迎风格式，压力-速度耦合采用phase-coupled simple方法。初始条件设置：塔内净液位高0（03m），0以下G0，0以上G1。边界条件为：进 口采用速度进口边界条件；出口采用压力出口；中心轴采用轴对称边界条件；壁面采用标准壁面函数，气液相均为速度无滑移条件。数值模拟的初 始时间步长为0.0001 s，待计算稳定后调整为0.0005 s。

2 结果与讨论

本文的实验数据主要来源于文献[3]及文献[9]，文献[3]做了关于500 mm鼓泡塔中充分发展段的气含率及液速径向分布的实验研究，阐述了列管内构件产生“烟囱效应”的力学机理，文献[9]在500 mm鼓泡塔中补充测定了表观气速0.12、0.31、0.62 m·s-1时含40根均匀布置垂直列管及空塔中不同轴向高度的气含率及轴向液速分布数据，实验表明除“烟囱效应”外，列管内构件对鼓泡塔中流动发展存在阻碍作用，主要表现为列管束阻碍发展区的气泡运动，导致分布器影响区延长。

2.1 气含率模拟结果

图2为不同表观气速下列管塔与空塔中气含率分布的模拟结果。模型能清晰反映加入列管后鼓泡中流动的变化规律。在均匀进气的条件下，空塔中（0）气泡迅速向塔中心聚集，形成气含率中心高边壁低的稳定分布。而在均匀布置了40根垂直列管的塔中（40），可以观察到气泡缓慢向塔中心运动的过程，分布器影响区约为4倍塔径左右，在分布器影响区，气含率初始分布较为平坦，随着气泡不断上升，塔中心气含率逐渐升高，在较高的轴向位置气含率分布才趋于稳定，流体充分发展。此外，该模型也清晰地反映了列管塔的另一特性“烟囱效应”，列管塔中充分发展段中心气含率高出空塔20%～30%。从图中可以看出，在列管塔中，分布器影响区约占整个鼓泡塔50%左右的区域，对于列管塔的流体力学模拟，采用文献[7]一维模型仅能计算充分发展区的流动分布，是不合理的。本文全面地考虑了列管束的作用机制，能获得较准确的二维分布图形。

图2 不同表观气速下列管塔与空塔中气含率分布云图

(a)G0.12 m·s-1,40; (b)G0.31 m·s-1,40; (c)G0.62 m·s-1,40; (d)G0.12 m·s-1,0; (e)G0.31 m·s-1,0;(f)G0.62 m·s-1,0

图3列举了4个轴向位置不同表观气速下气含率径向分布的模拟结果与实验值对比，从图中可以看出，在1.23 m的轴向位置，气含率沿径向均呈双峰分布，在1.73 m处，中心气含率较为平坦。而当＞2.23 m时，气含率分布趋于稳定，呈中心高边壁低的抛物形分布。该模型对不同表观气速下的气含率计算均与实验结果符合较好，在宽泛表观气速范围内（G0.12～0.62 m·s-1），计算结果与实验值相对误差在10%以内。

图3 不同表观气速下列管塔气含率分布模拟结果与实验值比较

2.2 轴向液速模拟结果

图4为列管塔与空塔轴向液速分布云图。从图中可以看出，鼓泡塔中液相呈现出中心向上、近壁面向下的流动循环。而空塔与列管塔的区别在于，在空塔中，液相速度在近分布器处即达到最大值，液相循环较为稳定，而列管塔中的液速呈逐渐发展的趋势，这与气含率的分布相似，证明了鼓泡塔中液相流动为气含率不均匀分布所驱动这一机理。列管的存在阻碍了气泡的径向运动，从而减缓了液相循环的发展，而在较高的轴向位置，气含率分布趋于稳定，当气相的宏观径向速度为零时，塔中心的轴向液速达到最大值，并不再随高度变化。在有列管存在时，充分发展后的轴向液速明显高于空塔，这是由于在充分发展区，列管的力学作用主要表现为对液相径向脉动的阻碍，导致湍流黏度下降，引发“烟囱效应”。

图4 不同表观气速下列管塔与空塔中轴向液速分布云图

(a)G0.12 m·s-1,40; (b)G0.31 m·s-1,40; (c)G0.62 m·s-1,40; (d)G0.12 m·s-1,0; (e)G0.31 m·s-1,0;(f)G0.62 m·s-1,0

图5为4个轴向位置不同表观气速下液速径向分布的模拟结果与实验值对比，模拟得到的规律与实验结果一致。随着轴向高度的增加，中心液速不断增大，液速分布逐渐趋于陡峭，当＞2.23 m时充分发展。模型对于各个条件下的液速模拟均与实验结果相符，最大相对误差不超过5%，表明本文提出的关于列管塔的流体力学模型适用条件较广，模型中考虑列管的力学机理准确，对含内构件鼓泡塔流型预测具有一定意义。

图5 不同表观气速下列管塔轴向液速分布模拟结果与实验值比较

图6为在500 mm鼓泡塔中安装不同数目列管时充分发展段的液速分布，从图中可以看出，该模型对不同管间距条件下的模拟结果也与实验值相符，表明该模型对于均匀布置列管束的鼓泡塔流型预测具有普遍适用性。模型中通过实验数据确定的参数C3适用于均匀布置换热列管束，列管数＜40 , 表观气速G为0.12～0.62 m·s-1湍流鼓泡区的鼓泡塔流体力学计算。

图6 不同列管密度条件下液速分布模拟结果与实验值比较

通过模型计算，可以更好地理解分布器影响区延长以及“烟囱效应”产生的原因。列管的加入导致流体的径向和切向运动受到阻碍，使得气泡更趋向于轴向直线形浮升，初始相对均匀的气含率分布能够沿列管方向扩展到很大范围，流动难以达到充分发展，导致分布器影响区延长；列管对液体湍动的阻碍作用使得液相湍流黏度下降，这一点得到了相关文献实验测量的证实[1,3-4]，由于鼓泡塔内液体的运动是由密度差引起的剪切流动，在气含率分布一定的情况下，液体湍流黏度的下降将导致速度径向梯度的增加，表现为在列管塔的充分发展段，液速分布比空塔更为陡峭，形成所谓的“烟囱效应”[3]。

2.3 源项影响及参数敏感性分析

本文在空塔模型的基础上，引入了3个动量源项以及两个湍动源项，分别为IG、kf,Lkd,Lkrsz其中IG作用于气泡相，主要描述列管塔中分布器影响区延长的现象，kf,Lkd,Lkrsz作用于液相，主要描述充分发展段的“烟囱效应”。为了解各个源项对流动影响的强弱程度，取表观气速G0.31 m·s-1、列管数目40的条件为计算示例，在模型中摒弃单个源项，考察流场的变化情况。图7为气泡阻尼项IG0时，塔内气含率及轴向液速分布，从图中可以看出，从分布器向上，流动迅速进入充分发展区，这与实验情况不符，但是气含率及液速在充分发展区的分布与图2 (b)及图4 (b)一致，表明该源项仅作用于分布器影响区，其表达式(3)中的阻力系数IG为经验参数，通过实验确定，适用于均匀布置列管束的鼓泡塔，列管数＜40, 表观气速0.12 m·s-1＜G＜0.62 m·s-1。

图7 源项FIG0时气含率及液速的分布云图

图8为在模型中分别不考虑源项kf,Lkd,Lkrsz时，充分发展区液速、湍流动能以及湍流黏度的变化趋势，从图8 (a)中可以看出，仅有kr对液速分布有显著影响，当kr0时，中心液速下降50%左右，而图8 (b)中5组湍流动能的分布曲线并没有明显差异，这表明，kr项的物理意义是列管表面与脉动的相互作用导致耗散项增大，从而使湍流黏度降低，从图8 (c)中可知kr项的加入会使湍流黏度减小20%左右。

图8 不考虑单个源项时液速、湍流动能及湍流黏度的分布

由以上分析可知，对列管塔充分发展区的流动参数分布起决定性作用的是列管束对脉动速度的形体阻力，因此需要着重分析与形体阻力相关参数的敏感性，而从式(4)可知，参数c与充分发展区的液速分布无关，所以仅需要考察K及C3的影响，取表观气速G0.31 m·s-1、列管数目40为条件进行计算分析，结果如图9所示，当K增大或减小50%时，中心液速分别增大31%或减小20%，而当C3增大或减小50%时，中心液速分别减小61%或增大117%，因此，C3是本文列管塔流体力学模型中最敏感的参数。

图9 CK及Cε3对液速分布的影响

3 结 论

（1）本文将列管束的作用考虑为对气液两相的作用力以及湍动源项，在描述“烟囱效应”的同时，阐释了带列管内构件鼓泡塔分布器影响区延长的力学机理。

（2）该模型在宽泛的表观气速（0.12～0.62 m·s-1）条件下，计算给出的气含率以及液速的二维分布与实验值吻合良好，能准确地描述含有列管内构件鼓泡塔的流动特性。

科技型中小企业技术创新投资二元平衡度与创新绩效关系研究…………………………………………………………… 张旭军，杨 勇，张 楠，段习贤，刘伦斌（69）

符 号 说 明

CIG——阻力系数 CL——升力系数 CTD——湍流扩散力系数 ct——模型参数 Fkd,L——列管对液体的形体阻力，N·m-3 Fkf ,L——列管对液体的摩擦阻力，N·m-3 g——重力加速度，m·s-2 k——湍流动能，m·s-2 P——压力，Pa R——塔半径，m r——径向位置，m uL——液体轴向平均速度，m·s-1 VG——表观气速，m·s-1 VS——气泡滑移速度，m·s-1 Z——轴向高度，m αG——局部气含率 ε——湍流耗散率，m2·s-3 ρL——液体密度，kg·m-3 ρG——气体密度，kg·m-3 ρm——气液混合物密度，kg·m-3 μt——湍流黏度，Pa·s 下角标 G——气相 L——液相 m——混合相

References

[1] Bernemann K. Zur fluiddynamik und zum vermischungsverhalten der fluessigen phase in blasensaeulen mit laengsangestroemten rohrbuendeln [D]. Dortmund: University Dortmund, 1989

[2] Chen J, Li F, Degaleesan S, Gupta Pt, Al-Dahhan M H, Dudukovic M P, Toseland B A. Fluid dynamic parameters in bubble columns with internals[J]...., 1999, 54: 2187-2197

[3] Zhang Yu (张煜), Lu Jia (卢佳), Wang Lijun (王丽军), Li Xi (李希). Studies on hydrodynamics of turbulent slurry bubble column (Ⅲ): Effect of vertical pipe bundles [J].(化工学报), 2009, 60 (5): 1135-1140

[4] Forret A, Schweitzer J M, Gauthier T, Krishna R, Schweich D. Liquid dispersion in large diameter bubble columns, with and without internals [J]....., 2003, 81: 360-366

[5] Youssef A A, Hamed M E, Grimes J T, Al-Dahhan M H, Dudukovic M P. Hydrodynamics of pilot-scale bubble columns: effect of internals [J]....., 2013, 52 (1): 43-55

[7] Zhang Yu (张煜), Li Hongbo (李红波), Li Zhaoqi (李兆奇), Wang Lijun (王丽军), Li Xi (李希). Hydrodynamics of turbulent slurry bubble column (Ⅳ): Modeling and simulation of bubble column with vertical pipe bundle [J].(化工学报), 2011, 62 (12): 3373-3380

[8] Forret A, Schweitzer J M, Gauthier R,. Influence of scale on the hydrodynamics of bubble column reactors: an experimental study in columns of 0.1,0.4 and 1m diameters [J]., 2003, 58: 719-724

[9] Chen J, Li F, Degaleesan S,. Fluid dynamic parameters in bubble columns with internals [J]., 1999, 54: 2187-2197

[10] Zhao Yuanfang (赵远方), Li Zhaoqi (李兆奇), Guan Xiaoping (管小平), Wang Lijun (王丽军), Cheng Youwei (成有为), Li Xi (李希). The effect of pipe bundles on the sparger zone development in bubble columns [J].(高校化学工程学报), DOI：10.3969/j.issn.1003-9015.2015.01.30.01

[11] Laborde-Boutet C, Larachi F, Dromard N, Delsart O. CFD simulations of hydrodynamics/thermal coupling phenomena in a bubble column with internals [J].., 2010, 56: 2397-2411

[12] Li Zhaoqi (李兆奇), Wang Lijun (王丽军), Guan Xiaoping (管小平), Cheng Youwei (成有为), Li Xi (李希). A 2-D hydrodynamic model of bubble column based on lateral-force balance [J].(化工学报), 2014, 65 (11): 4222-4230

[13] Robert H Perry, Don W Green. Perry’s Chemical Engineer’s Handbook [M]. 7th ed. Beijing: Science Press, 2001

[14] Antal S P, Lahey R T, Flaherty J E. Analysis of phase distribution in fully developed laminar bubbly two-phase flow [J]., 1991 (7): 635-652

[15] Tomiyama A. Struggle with computational bubble dynamics//Third International Conference on Multiphase Flow, ICM'F’98 [C]. Lyon, France, 1998

[16] Lucas D, Krepper E, Prasser H M. Use of models for lift, wall and turbulent dispersion forces acting on bubbles for poly-disperse flows [J]., 2007, 62 (15): 4146-4157

[17] Auton T R. The lift force on a spherical body in a rotational flow [J]..., 1987, 183: 199-218

[18] Thomas N H, Auton T R, Sene K, Hunt J C R. Entrapment and transport of bubbles by transient large eddies in multiphase turbulent shear flows//Stephens H S, Stapleton C A, eds. International Conference on Physical Modelling of Multi-Phase Flows [C]. Cranfield, BHRA Fluid Engineering, 1983:169-184

[19] Lopez de Bertodano M. Turbulent bubbly two-phase flow in a triangular duct [D]. New York: Rensselaer Rolytechnic Institute, 1992

CFD simulation of bubble column with vertical bundle internals

LI Zhaoqi, ZHAO Yuanfang, GUAN Xiaoping, WANG Lijun, CHENG Youwei, LI Xi

Department of Chemical EngineeringZhejiang UniversityHangzhouZhejiangChina

For gas-liquid bubble column with dense vertical bundle internals, the resistance force of bundles on gas and liquid phases was taken into account in the two-phase Euler axisymmetric-model. By introduction of corresponding momentum source, turbulence source and dissipation source, a 2-D CFD model of bubble column with bundle internals was established. The model could clearly and accurately demonstrate the features of gas-liquid flow in bubble column with bundles: “funnel effect” and extension of sparger zone. The calculated two-dimensional distributions of hold-up and liquid velocity agreed with experimental data within a wide range of superficial gas velocity (0.12—0.62 m·s-1).

bubble column; dense vertical bundles; internals; CFD; model

date: 2014-06-12.

WANG Lijun, wang_lijun@zju.edu.cn

10.11949/j.issn.0438-1157.20140884

TQ 021.1

A

0438—1157（2015）03—0932—10

国家高技术研究发展计划项目（2011AA05A205）；国家自然科学基金项目（U1162125，U1361112）。

2014-06-12收到初稿，2014-12-16收到修改稿。

联系人：王丽军。第一作者：李兆奇（1985—），男，博士研究生。

supported by the National High Technology Research and Development Program of China (2011AA05A205) and the National Natural Science Foundation of China (U1162125, U1361112).