郭旭桓, 张子木, 赵秋月, 张廷安
东北大学 多金属共生矿生态化冶金教育部重点实验室, 辽宁 沈阳 110819
收稿日期：2016-10-31
基金项目：国家自然科学基金资助项目(51204040, U1202274)；中央高校基本科研业务费专项资金资助项目(N130607001, N140204017)；辽宁省教育厅科学技术研究项目(L2014096)；辽宁省自然科学基金资助项目(2015020591)。
作者简介：郭旭桓(1986-), 男, 山西洪洞人, 东北大学博士研究生;
张廷安(1960-), 男, 河南周口人, 东北大学教授，博士生导师。

摘要：采用数值模拟方法研究了压力能驱动的自搅拌反应器内液相流场的分布规律.首先对标准κ-ε模型、RNG κ-ε模型、realizable κ-ε模型下的模拟结果与PIV实验结果进行对比, 确定了合适的湍流模型; 然后研究了不同入口压力、液位高度对反应器内流场分布的影响.结果表明:标准κ-ε模型能更精确地模拟该反应器内的流场分布, 反应器内流场稳定时间为12 s以上.增加入口压力和液位高度均有利于反应器内流体速度均匀分布, 入口压力为3 MPa时几乎不存在速度死区.
关键词：CFD(计算流体力学)速度分布湍动能自搅拌管式反应器
Numerical Simulation of Liquid Phase Flow in a Self-Stirring Reactor
GUO Xu-huan, ZHANG Zi-mu, ZHAO Qiu-yue, ZHANG Ting-an
Key Laboratory of Ecological Metallurgy of Multi-metal Intergrown Ores Ministry of Education, Northeastern University, Shenyang 110819, China
Corresponding author: ZHAO Qiu-yue, E-mail: zhaoqy@smm.neu.edu.cn
Abstract: The distribution of liquid phase flow field in a self-stirring reactor driven by pressure energy was studied by numerical simulation method. Firstly, the appropriate turbulence model was determined by comparing the simulation results of the standard k-ε model, RNG k-ε model, and realizable k-ε model with the PIV experimental results.Then the influence of inlet pressure and liquid level height on the fluid field distribution in the reactor was studied. The results showed that the standard k-ε turbulence model can accurately simulate fluid field distribution of the reactor, and the stable time is more than 12 s. Increasing inlet pressure and liquid height is favorable to the even distribution of fluid velocity in the reactor, and there is almost no dead zone of velocity when inlet pressure is 3 MPa.
Key Words: CFD(computational fluid dynamics)velocity distributionturbulent kinetic energyself-stirringtubular reactor
搅拌反应器广泛应用在湿法冶金的配料、浆化、浸出、结晶、溶解、还原、分解和萃取等工序中.北京有色冶金设计研究总院对全国有色系统冶炼厂搅拌设备的调查和统计结果表明, 许多湿法冶炼车间50%以上的功率是消耗在搅拌作业上^[1].高效搅拌反应器的研究一直是学术界和工业界关注的热点.
2008年Kasat等^[2]采用Fluent软件, 模拟研究了带挡板的泥浆搅拌槽内固-液两相流的固体浓度分布、搅拌功率、轴向速度分布及混合时间, 与实验数据吻合良好.2012年Wadnerkar等^[3]研究了带挡板的搅拌槽内Gidaspow, Wen and Yu, Brucato和改进的Brucato四种不同曳力模型下的云高度及湍动能, 结果表明改进的Brucato模型与实验结果更加吻合.2013年Qi等^[4]研究了带挡板的搅拌槽内固-液两相流动规律, 模拟结果表明:功率准数随着固体浓度的增加而增加, 固体颗粒平均浓度的增加和粒径的增加均不利于其在液相中的均匀分布, 黏度的增加有利于固体颗粒的均匀分布, 搅拌桨正下方底部存在一定的死区.2014年Zhao等^[5]研究了改进的Intermig桨搅拌槽内的流场分布, 结果表明改进的Intermig桨改善了搅拌槽内流场的流型, 降低了功耗.2016年Yang等^[6]研究了圆柱形椭圆底搅拌槽内气液固三相流动行为, 结果表明:改进的气泡尺寸模型和曳力模型更好地预测了气泡在反应器内的分布规律.2015年Qiao等^[7]采用CFD方法, 研究了圆柱形平底搅拌槽内固-液悬浮临界转速, 得到了在接近临界悬浮转速时功率准数、局部固含率、标准偏差、云高度、无因次混合时间和固体颗粒轴向速度6个主要参数的变化规律.
CFD数值模拟技术结合物理实验验证成为研究搅拌反应器的主要方法.目前的研究大多集中在搅拌槽反应器上, 对管式搅拌反应器的研究较少.东北大学张廷安^[8]课题组围绕节能, 从传热传质、化学反应和抑制结疤三个主题原创性地开发了一种新型的加压反应器-压力能驱动的自搅拌反应器, 已获得发明专利授权(授权公告号为:CN103495374B).张子木等^[9-10]初步研究了液位高度对流场的影响, 对搅拌转速与原动力-压力能的关系没有作细致的研究.本文对该反应器的从动桨桨型作了进一步改进, 并借助FLUENT15.0软件中Eulerian两相流模型、标准κ-ε湍流模型详细研究了不同入口压力、液位高度等工艺条件对反应器内流场的影响规律, 为该反应器的进一步开发奠定基础.
1 几何模型与数值模拟方法1.1 几何模型实验用反应器长L=954 mm, 内径D=200 mm, 进出口管直径d₁=16 mm, 框式从动桨叶片顶端离管壁d₂=5 mm, 桨叶高度h₁=95 mm, 宽度h₂=142 mm, 间距W=112.5 mm, 涡轮搅拌桨为16叶圆盘桨, 涡轮直径d₃=145 mm, 涡轮叶片宽h₃=40 mm.采用Workbench15.0软件按实验用反应器1:1绘制几何模型，如图 1所示.
图 1(Fig. 1)

图 1 几何模型示意图(mm)Fig.1 Schematic diagram of geometric model

1.2 网格划分用ICEM15.0划分四面体、六面体混合网格, 对搅拌桨区域和进出口处的网格进行加密处理, 根据不同模型, 网格数量在30万至120万不等.图 2所示为反应器和搅拌桨网格模型.
图 2(Fig. 2)

图 2 网格示意图Fig.2 Schematic diagram of grid (a)—反应器网格划分；(b)—涡轮桨的网格划分；(c)—从动桨的网格划分.

1.3 边界条件及数值解法用FLUENT15.0有限体积法离散并求解控制方程, 选用速度压力耦合的SIMPLE算法, 二阶迎风格式, 所有项的残差收敛均设为10^-5.流体流动设为非定常流动, 用滑移网格法处理搅拌桨叶与静止区域的相对运动.用Eulerian模型模拟不同液位高度的流场特性, 用标准壁面函数法处理近壁区域, 涡轮桨和搅拌轴及桨叶设为移动壁面边界条件, 筒体设为静止壁面边界条件, 具体设置如下:入口为压力入口, 压力值为0.5~3 MPa; 出口为自由出口.根据转速(单位rad/s)倒数的1/10的时间步长选取方法, 本研究中考察的搅拌转速介于50~150 r/min, 推算出时间步长设为0.005 s较好.
2 模拟结果与分析2.1 模型验证为了能够同时满足模拟结果的准确性和计算时间的合理性, 首先对30万~120万的不同网格密度进行网格无关性检查, 结果发现95万网格模拟结果和计算时间均较为合理, 最终采用95万网格用于本文的模拟计算.
分别用标准κ-ε模型、RNG κ-ε模型、realizable κ-ε模型三种湍流模型, 在相同操作参数下进行计算模拟, 与同条件下PIV实验所测结果进行对比.取搅拌轴正前方2 cm截面上中间2个桨叶区域为PIV测量区域, 所测流场矢量图与模拟结果对比见图 3.
图 3(Fig. 3)

图 3 不同模型模拟结果与PIV实验结果对比图Fig.3 Comparison charts of PIV results and the simulation results of different models (a)—标准κ-ε模型；(b)—PIV实验结果; (c)—RNG κ-ε；(d)—realizable κ-ε.

从图 3可以看出, 标准κ-ε模型的模拟结果更接近实验结果.本文后边的模拟均选用标准κ-ε湍流模型.
2.2 同一转速不同时刻的速度分布为了方便分析, 首先将流场划分为四个区域, 反应器左端壁面到涡轮搅拌桨之间为区域1, 0 < Z/L < 0.11;涡轮搅拌桨到框式从动桨之间为区域2, 0.11 < Z/L < 0.21;框式从动桨搅拌区域为区域3, 0.21 < Z/L < 0.73;其余区域为区域4, 0.73 < Z/L < 1.为了对反应器内的流动情况进行定量分析, 分别选取了3条轴向监测线和7个径向监测面, 并进行了无因次处理，具体如图 4所示.
图 4(Fig. 4)

图 4 监测线及监测面示意图Fig.4 Schematic diagram of the different traverses and sections

首先对液面高度为200 mm(满水状态)、入口压力为2 MPa情况下反应器内流体的流动情况进行了数值计算, 得到不同时刻轴向中心监测面上的速度分布如图 5所示.
图 5(Fig. 5)

图 5 不同时刻轴向中心监测面上的速度分布Fig.5 The velocity distribution in the axial center section at different time (a)—t=0.5 s；(b)—t=1 s; (c)—t=3 s；(d)—t=9 s; (e)—t=15 s；(f)—t=30 s.

由图 5a~图 5c可以看出, 由于时间不同, 搅拌桨叶的位置不同, 虽然整体上流体都是由入口流向出口, 但是不同区域的局部流场有些区别, 在区域1, 不同时间均保持几乎相同的速度分布和流场分布, 这可能是由于涡轮桨为16叶对称结构, 可以认为该区域内流体流动行为不随时间变化.在区域1内, 由于入口高压流体的冲击和涡轮桨转动的双重作用, 流体在沿轴向流动的同时, 在搅拌轴上、下侧各形成一个轴向环流, 有效抑制了左端角部区域的死区形成.在区域2, 流体均沿着轴向流动.在区域3, 由于桨叶的转动和排水作用, 在4个桨叶周围均形成一个轴向环流, 环流区域的大小随着桨叶的转动周期性地变化, 桨叶附近的环流区域较大.在区域4, 搅拌轴上方存在一个大的轴向环流区, 保证在没有搅拌桨叶的情况下, 该区域也不会形成死区.从图 5d~图 5f可以看出, 在不同时刻桨叶位置相同的情况下, 流场的速度分布及流型几乎完全相同, 该反应器内流动行为呈周期性变化.除入口与出口外, 在搅拌桨叶的带动作用下, 其附近的流体速度较大, 搅拌轴附近的流体速度较小, 区域4的轴向环流中心存在速度死区.由图中可以看出9 s以后的流场基本已经稳定, 故后文中计算模拟时间均为15 s.
2.3 不同入口压力下的速度分布2.3.1 轴向速度分布为了准确反映反应器内流体的轴向速度分布, 分别取径向高度为Y/D=0.2, 0.75和0.9三条轴向监测线上的速度分布, 见图 6.
图 6(Fig. 6)

图 6 不同径向高度的轴向流体速度变化曲线Fig.6 Profiles of the axial liquid velocity at different radial heights (a)—Y/D=0.2；(b)—Y/D=0.75；(c)—Y/D=0.9.

从图 6中可以看出, 3条轴向监测线上的速度分布表现出相似的变化规律:区域1, 流体在壁面边界层内速度为零, 离开边界层后流体轴向流动到涡轮桨壁面后回流, 加之涡轮桨的排水作用, 在该区域形成了两边速度大、中间速度小的流动行为, 在Z/L为0.025处速度达到最小, 涡轮搅拌桨处在入口正下方, 在入口处高压流体的驱动下该处为反应器内速度最大点; 区域2, 在搅拌桨叶轴向环流的影响下, 流体速度逐渐减小, 3条监测线均在Z/L为0.2左右达到最小, 与径向距离无关, 从动桨的左端边缘为Z/L=0.21;区域3, 搅拌轴上方流体受到桨叶环流的影响, 速度逐渐减小, 在到达下个桨叶边端时达到最小, 搅拌轴下方区域速度基本稳定; 区域4, 流体速度基本稳定, 径向距离Y/D=0.75的轴向监测线因为经过该区域的轴向环流中心, 在Z/L为0.85处出现速度极小值, 搅拌轴下方监测线受出口影响, 在Z/L为0.97处出现速度极低点.
比较不同搅拌转速的轴向监测线可以看出, 搅拌转速的增加不仅可以提高各个区域的轴向速度, 还可以明显改善区域4内搅拌轴上侧环流中心的速度死区.各径向监测面的平均速度的轴向分布也表现出了该规律, 见图 7.
图 7(Fig. 7)

图 7 不同转速下的轴向平均速度Fig.7 Axial velocity at different rotation speeds

搅拌转速的增加明显改善了反应器内的湍流状态, 使速度分布更加均匀, 不断减小直至消除速度死区的范围.
2.3.2 径向速度分布为了研究反应器内流体的径向流动规律，取轴向距离Z/L=0.55位置处径向监测面不同转速下的速度分布，如图 8所示.从图 8可以看出, 该监测面上流体主要做绕轴环流, 在搅拌桨叶与轴之间区域形成一个小漩涡.搅拌转速的增加对流场的流型影响较小, 但使得监测面内的速度分布更加均匀, 且明显减小了速度死区.
图 8(Fig. 8)

图 8 不同转速下Z/L=0.55处径向监测面的速度分布Fig.8 Velocity distribution of the section Z/L=0.55 at different rotation speeds (a)—65 r/min；(b)—100 r/min；(c)—150 r/min.

图 9给出了该径向监测面上相应转速下的湍动能分布.从图中可以看出低转速时监测面上桨叶附近流体的湍流强度较大, 其余区域均较小.随着搅拌转速的增大, 桨叶附近和其前一时刻转动区域的流体湍流强度增加明显, 这能够有效抑制反应器边端死区的形成.
图 9(Fig. 9)

图 9 不同转速下Z/L=0.55处径向监测面的湍动能分布Fig.9 Turbulent kinetic energy distribution of the section Z/L=0.55 at different rotation speeds (a)—65 r/min；(b)—100 r/min；(c)—150 r/min.

2.4 不同液位高度下的速度分布本节主要研究在固定的搅拌转速下, 液位高度Y分别为100, 133, 167 mm(Y/D=0.5, 0.67, 0.83)时反应器内流体的流动行为.数值模拟所需的入口压力边界条件由实验获得.
图 10为径向高度Y/D分别为0.04和0.2的2组轴向监测线上的速度分布.从图 10a中可以看出, 三种液位高度情况下监测线上流体的速度分布规律大体相同, 由于靠近反应器底部, 除靠近出口处流体流速较大外, 流速基本保持稳定没有太大波动.从图 10b可以看出, 三种不同条件表现出了不同的流动规律, 在区域1, 由于液位较低, 100, 133 mm两种条件下的流体受到入口高压流体的影响较小, 流体速度较为稳定; 167 mm影响较大, 速度逐渐减小.由于三者转速相同, 在所以在涡轮搅拌桨处速度相同.区域2三者均呈现出速度不断增加的趋势, 在Z/L=0.18左右达到极大值.在区域3, 受搅拌桨叶的驱动作用, 在桨叶附近流体速度均匀稳定, 在两桨叶之间速度先减小后增大.在区域4, 133, 167 mm两种条件下流体受出口的影响, 在出口处速度较大, 出口两端速度逐渐减小, 100 mm条件下的流体由于液位高度较低, 在反应器搅拌的湍动影响下, 出口附近流体液面不平, 导致出现一个速度极低点.
图 10(Fig. 10)

图 10 不同径向高度的轴向液速变化曲线图Fig.10 Profiles of the axial liquid velocity at different radial heights (a)—Y/D=0.04；(b)—Y/D=0.2.

图 11为不同液位高度下7个径向监测面的轴向平均速度分布, 从图中可以看出, 液位高度越高, 轴向平均速度分布越平稳.
图 11(Fig. 11)

图 11 不同液位高度下的轴向平均速度Fig.11 Mean axial velocity at different liquid heights

3 结论1) 对比不同湍流模型模拟结果与PIV实验结果的速度矢量图, 确定了标准κ-ε模型能更精确地模拟该反应器内的流动规律.
2) 分析了搅拌转速为100 r/min条件下不同时间的流场分布, 确定反应器内流场稳定时间为12 s以上.
3) 模拟了不同入口压力下流体的流动规律, 发现反应器内整体的速度分布较为平稳, 入口压力的增加有利于流体速度的均匀分布, 并能减小直至消除流动死区的范围.
4) 液位高度的增加有利于反应器内流体速度的均匀分布.
5) 通过对该反应器内单相流的宏观流动行为进行模拟, 对该反应器的流体动力学特性有了科学的了解, 同时也为该反应器内固-液两相流动行为的研究打下了基础.
参考文献

[1]	厉衡隆, 顾松青, 李金鹏, 等. 铝冶炼生产技术手册[M]. 北京: 冶金工业出版社, 2011: 238-400. ( Li Heng-long, Gu Song-qing, Li Jin-peng, et al. Technical manual of aluminum smelting[M]. Beijing: Metallurgical Industry Press, 2011: 238-400.)
[2]	Kasat G R, Khopkar A R, Ranade V V, et al. CFD simulation of liquid-phase mixing in solid-liquid stirred reactor[J].Chemical Engineering Science, 2008, 63(15): 3877–3885.DOI:10.1016/j.ces.2008.04.018
[3]	Wadnerkar D, Utikar R P, Tade M O, et al. CFD simulation of solid-liquid stirred tanks[J].Advanced Powder Technology, 2012, 23(4): 445–453.DOI:10.1016/j.apt.2012.03.007
[4]	Qi N N, Zhang H, Zhang K, et al. CFD simulation of particle suspension in a stirred tank[J].Particuology, 2013, 11(3): 317–326.DOI:10.1016/j.partic.2012.03.003
[5]	Zhao H L, Zhang Z M, Zhang T A, et al. Experimental and CFD studies of solid-liquid slurry tank stirred with an improved intermig impeller[J].Transactions of Nonferrous Metals Society of China, 2014, 24(8): 2650–2659.DOI:10.1016/S1003-6326(14)63395-1
[6]	Yang S F, Li X Y, Yang C, et al. Computational fluid dynamics simulation and experimental measurement of gas and solid holdup distributions in a gas-liquid-solid stirred reactor[J].Industrial and Engineering Chemistry Research, 2016, 55(12): 3276–3286.DOI:10.1021/acs.iecr.5b03163
[7]	Qiao S C, Wang R J, Yang X X, et al. CFD simulation of complete drawdown of floating solids in a stirred tank[J].Canadian Journal of Chemical Engineering, 2015, 93(1): 141–149.DOI:10.1002/cjce.v93.1
[8]	张廷安. 一种自搅拌管式溶出反应器: CN103495374A[P]. 2014-01-08. ( Zhang Ting-an. A self-stirring tube type dissolution reactor: CN103495374A[P]. 2014-01-08. ))
[9]	张子木, 吕超, 赵秋月, 等. 压力能驱动的自搅拌反应器流动特性数值模拟[J].东北大学学报(自然科学版), 2015, 36(7): 962–965. ( Zhang Zi-mu, Lyu Chao, Zhao Qiu-yue, et al. Numerical simulation on flow characteristics of self-stirring reactor driven by pressure energy[J].Journal of Northeastern University(Natural Science), 2015, 36(7): 962–965.)
[10]	Zhang Z M, Zhao Q Y, Lyu C, et al. Characteristics of fluid flows in a novel self-stirring reactor[J].The Journal of the Minerals, Metals & Materials Society, 2014, 66(7): 1197–1201.