基于双流体模型的气液两相流三维瞬态数值模拟

基于双流体模型的气液两相流三维瞬态数值模拟

气液两相流模拟鼓泡塔检测器具有结构简单、操作方便、相距面大、传热传递质量高、传热传递质量好等特点。广泛应用于土木工程、生物工程、环境和能源工业领域。过去20年,有关鼓泡塔已经做了大量的研究工作,包括数值模拟和实验测量,但该类反应器的设计与放大仍然很困难。由于气体分散相,鼓泡塔内的气/液相流动十分复杂,到目前为止,文献报道的研究结果大多以宏观数据为主,目前的实验技术很难测量鼓泡塔内的流场细节。当气含率较高,流型为非均相流时,流场细节的测量则更为困难,因此,高表观气速下的局部气含率及气液相速度的实验数据报道很少。但仅凭整体宏观特性数据进行鼓泡塔的设计与放大很不可靠,整体属性不能敏锐地反映鼓泡塔的尺寸和结构的影响,只能通过局部属性来表征,因此,有必要获取鼓泡塔内的局部流场信息。如前所述,靠实验手段测量局部流场相当困难,随着计算流体力学的发展、计算机性能和数值模型可靠性的提高,依靠数值模拟获取更多的局部流场信息已经成为可能。本文将采用数值模拟对鼓泡塔流动进行三维瞬态分析,以获取更多的流场信息。鼓泡塔数值模拟属于气液两相流模拟。两相流模拟主要有两种方法,即双流体模型和欧拉-拉格朗日法。其中欧拉-拉格朗日法需跟踪每个气泡,对计算资源的要求很高,目前只能模拟低气含率的两相流;双流体模型对各相分别求解动量方程和组分方程,降低了计算量。本文采用双流体模型。目前双欧拉法进行气液两相流模拟存在两大困难,一是正确描述两相间的相互作用力;二是采用恰当的湍流模型模拟气液相湍流。关于相间作用力的模拟存在很大的争议,但通常认为曳力作用最为主要,Dudukovic等、Krishna等仅考虑曳力就获得了与实验相一致的时均结果。在曳力的描述中,气泡尺寸最为关键,然而大部分文献都采用单一气泡尺寸模型,这显然与实际情况不符。尽管通过选择合适的气泡尺寸可以获得与实验结果相吻合的时均结果,但是需要通过反复的尝试、筛选,明显降低了数值模拟的预测性。为了考虑气泡的尺寸变化,Ishii等[5,6,7,8,9,10,11,12]已经做了一定的工作,采用了不同的气泡尺寸模型。其中,Ishii等采用气泡界面传输模型来模拟气泡界面密度,而Sanyal等采用气泡簇平衡模型(PBM模型)预测气泡尺寸。气泡界面传输模型只需解一组气泡界面密度标量方程,而PBM模型需对每一组气泡求解组分标量方程,显然,后者比前者需要更多的计算资源与计算时间,特别是PBM模型的各气泡组分方程相互耦合,增加了收敛的难度。尽管最新的FLUENT商业软件中已经耦合了该PBM模型,但考虑到三维模拟的计算量,本文仍采用Ishii的单方程气泡界面密度模型。计算表明,该模型能大大降低对计算量的需求,并具有很好的收敛性,且计算结果与实验值比较吻合。液相湍流的模拟是多相流模拟的难点,不同的学者采用了不同的湍流模型。Ranade等采用单相标准k-ε模型,气相假定为层流。Sanyal等采用了考虑气泡影响的k-ε模型模拟液相,气相湍动描述则基于Tchen的理论。Zhou等采用二阶矩模型模拟气液相湍流,并考虑了两相之间的脉动相关。Zhang等采用大涡模拟对液相湍流进行模拟。可见,并没有普遍适用的湍流模型可供使用。鉴于此,本文仍采用修正k-ε模型,气相湍流模拟则基于Tchen的理论,重点考察单方程气泡模型能否预测鼓泡塔内的气液界面浓度。1数学模型1.1湍流黏度模型双流体模型控制方程可在单相流控制方程的基础上采用集平均方法推导。本文中假定液体为连续相,气体为分散相不可压缩,且不考虑相间质量和热量传递,控制方程如下:质量守恒方程∂αi∂t+∇⋅(αiUi)=0(1)∂αi∂t+∇⋅(αiUi)=0(1)动量守恒方程∂αiρiUi∂t+∇⋅(αiρiUiUi)=-αi⋅∇pi+∇⋅(αiτi)+Fi+αiρig(2)∂αiρiUi∂t+∇⋅(αiρiUiUi)=−αi⋅∇pi+∇⋅(αiτi)+Fi+αiρig(2)τi=(μt+μl)[∇Ui+(∇Ui)Τ-23Ι∇⋅Ui](3)τi=(μt+μl)[∇Ui+(∇Ui)T−23I∇⋅Ui](3)式中Fi为相间作用力,N;τi表示黏性应力张量,Pa;μl表示分子黏度,Pa·s;μt为湍流黏度,Pa·s。液相湍流黏度采用k-ε模型来计算。模型方程为μt‚l=Cμρlk2ε(4)μt‚l=Cμρlk2ε(4)∂αlρlkl∂t+∇⋅(αlρlklUl)=∇⋅(αlμt‚lσk∇kl)+αlGk‚l-αlρlεl+αlρlΠk‚l(5)∂αlρlkl∂t+∇⋅(αlρlklUl)=∇⋅(αlμt‚lσk∇kl)+αlGk‚l−αlρlεl+αlρlΠk‚l(5)∂αlρlεl∂t+∇⋅(αlρlεlUl)=∇⋅(αlμt‚lσε∇εl)+αlεlkl(Cε1Gk‚l-Cε2ρlεl)+αlρlΠε‚l(6)∂αlρlεl∂t+∇⋅(αlρlεlUl)=∇⋅(αlμt‚lσε∇εl)+αlεlkl(Cε1Gk‚l−Cε2ρlεl)+αlρlΠε‚l(6)式中Πk,l、Πε,l表示分散相引起的湍流附加作用项;参数Cε1、Cε2、Cμ、σε及σk为常数,分别取1.44、1.92、0.09、1.0、1.3。1.2气泡破碎过程气泡尺寸是模拟两相间动量传递的关键参数,气泡在鼓泡塔内的上升过程中将经历多次聚并与破碎,尺寸将发生显著变化。当气泡所受到的外界应力大于表面张力时,气泡破碎,而当两个或多个气泡相互碰撞且气泡间液膜小于一定厚度时将聚并。根据Ishii等推导的气泡数密度传输方程,气泡的界面密度传输标量方程可表示为∂ρa∂t+∇⋅(ρaUg)=Sbk-Sco(7)∂ρa∂t+∇⋅(ρaUg)=Sbk−Sco(7)d=6α/a(8)式中Sbk表示气泡破碎导致的气液界面增加,根据气泡破碎理论,气泡破碎的主要机理是液相湍流涡碰撞,但只有尺寸小于气泡的湍流涡的碰撞才有效。气泡破碎速率取决于气泡与湍流涡的碰撞频率,并且气泡Weber数必须大于某一临界值气泡才会破碎,Wu等依据气体分子动理论,推导出气泡碰撞速率表达式Sbk=118Cbkut(a2α)(1-WecrWe)1/2exp(WecrWe)(9)Sbk=118Cbkut(a2α)(1−WecrWe)1/2exp(WecrWe)(9)We=ρutdσ(10)ut=√2(εd)1/3(11)式中ut表示处于湍流能谱惯性子区的湍流涡速度。式(7)中Sco表示气泡聚并所引起的相界面减少。Wu等认为气泡聚并主要因气泡之间的碰撞产生,碰撞频率与湍流脉动有关,若只考虑气泡之间的二元碰撞,聚并速率可以表示为Sco=-Cco(uta2)3π×[1α1/3max(α1/3max-α1/3)]{1-exp[-α1/3max(αmax/α)1/3-1]}(12)Sbk与Cco为可调节参数,本文中取0.3;Wecr数根据文献取0.18;αmax为最大气含率,取0.8。式(11)中耗散率可通过求解k-ε方程获取。本文仅考虑两相之间的曳力作用,曳力系数采用Schiller&Naumann模型FD‚g=-0.75αg(1-αg)ρlCDdB|Ug-Ul|(Ug-Ul)CD={24(1+0.15Re0.687)/ReRe≤10000.44Re>1000(13)2速度进口边界本文模拟一圆柱形鼓泡塔,塔径D=190mm,净液位高度H=1000mm,分布器采用多孔分布结构,孔径1mm,网格结构如图1所示。模拟中,鼓泡塔顶部采用压力出口边界条件,底部采用速度进口边界,壁面为非滑移边界。初始时间步长Δt=0.001s,待计算达到稳定后取Δt=0.005s。离散格式时间导数采用一阶迎风格式,空间导数采用Quick格式,压力-速度耦合采用phase-coupledsimple方法,压力松弛因子为0.7,动量方程松弛因子取0.3,其余采用默认值。气泡界面密度传输模型通过在Fluent中添加一组标量传输方程来实现,标量方程的聚并项和破碎项可通过用户自定义程序来实现。待流动充分发展并达到拟稳态后,即可对特定位置上的瞬时物理量进行记录,并保存在外部文件中,物理量的实时记录可通过UDF程序来实现。3结果与讨论3.1气含率分布模拟结果本文模拟了0.02、0.05、0.10、0.12、0.15、0.20、0.25、0.30m·s-18个表观气速,涉及均相流型与非均相流型。图2描述了表观气速0.12m·s-1时中心轴截面的液相速度流型,箭头长度表示速度大小,颜色表示气含率。由图可知,本文采用的模型有效地揭示出鼓泡塔内气液两相流的瞬态行为,特别是大涡的发展,这种大尺度涡的发展呈现明显摆动。整体气含率是表征鼓泡塔性能的重要参数,也是衡量数值模拟结果是否可靠的重要依据。图3虚线表示整体气含率模拟值,实线为实验结果,可以看出模拟结果与实验结果基本吻合。模拟曲线反映了气含率随表观气速增加而增大的趋势,但随着表观气速增加,气含率的增加幅度减小,若将气含率表示为表观气速的函数,则指数因子应小于1。气含率分布使鼓泡塔在径向方向形成密度差,推动了整体循环流动的形成,因此,预测局部气含率分布十分重要。图4描述了表观气速为0.12m·s-1时气含率的径向分布,由图可知,随着高径比的增大,塔中心气含率逐步提高,分布曲线斜率增大,这是因为随着高径比的增大,流型逐渐发展,多孔分布器形成的小气泡经过不断的聚并,形成越来越多的大气泡,并逐步向塔心聚集所致。但高径比为2、3处的气含率分布差别已很小,说明流动已发展充分。另外,图4对气含率分布模拟值与文献实验结果进行了比较。尽管文献取于H/D=2.8处,而模拟结果取于H/D=3处,但可以认为这两个位置上的气含率基本相同。可以看出,模拟结果与实验结果吻合较好,但塔中心模拟值比实验值稍高,随着远离中心轴,模拟值反而低于实验值,但误差都小于10%,进一步靠近壁面,模拟结果与实验结果很接近。另外,模拟曲线相对光滑,而实验曲线则显得过于突兀。图5为不同表观气速下的气含率径向分布,可以看出,局部气含率随表观气速的增大而增大,特别是塔中心,局部气含率增加明显。当表观气速较高时(>0.15m·s-1),壁面附近的气含率增加不明显,甚至有所降低,这是因为壁面处的气含率主要是由回流气泡组成,当表观气速增大时,大气泡增加,小气泡减少,导致壁面附近的回流气泡减少,气含率降低。另外,表观气速很低时(Ug=0.02m·s-1),气含率分布相当均匀,说明此时正处于均相流型,随着气速的增大,气含率分布曲线越来越陡,说明流型逐步由均相流转变为非均相流。3.2表观气速对液相速度的影响图6为中心线轴向液相速度Riquarts关系式预测值与模拟结果的比较。从图中可以看出,模拟值与关系式预测值吻合较好。当气速低于0.15m·s-1,模拟值比计算值要高10%左右,而当气速高于0.15m·s-1,模拟结果稍低于计算值,但都处于工程问题所能接受的误差范围内。图7给出了表观气速为0.12m·s-1时的轴向液相速度径向分布。由图可以看出,塔中心处轴向液相速度随高径比的增大而增大,但H/D>2以后,速度增加值减小,曲线斜率变化也减小,说明H/D>2以后流动逐渐充分发展。另外,该图还对模拟结果与文献的实验值进行了定量比较,可以看出,模拟结果定性地反映了速度径向分布规律,模拟较合理地预测了液相速度拐点。但鼓泡塔的中心模拟值比实验值略低,进一步靠近壁面,误差有所增大。目前,有关鼓泡塔数值模拟的文献中仅有少部分工作对气相速度进行了预测,并与实验结果进行比较,大多仅预报液相速度,主要是由于气相的速度不易测量,特别是表观气速较大时(>0.1m·s-1),气相速度测量难度非常大。目前对非均匀两相流气相速度数据报道的甚少,因此,本文仅给出表观气速为0.12m·s-1时轴向气相速度模拟结果(图8)。图8可以根据气相对液相的滑移来解释,对于一定尺寸的气泡,气/液相间滑移速度Δu与气泡在静止液体中的上升速度有相同的量级,若在图7的液相速度曲线上叠加一滑移速度,则可获得图8。从图中可以看出,随着高径比增大,气相速度逐渐增大,特别是塔中心(r/R<0.6)区域气相速度增加明显,这是因为气泡在上升的过程中聚并形成了大气泡,上升速度加快,而在壁面附近气泡尺寸变化不大,以小气泡为主(大气泡向塔中心运动),因此壁面附近气相速度沿轴向没有大的变化。图8中气相轴向速度都大于零,说明从长时均结果来看,气相在各个空间位置上均呈向上运动,但并不表示气泡没有回流,事实上壁面附近空间点轴向气相速度的时间监控曲线中仍存在负值,说明仍有气泡被液体卷吸回流。图9显示了不同表观气速下时均轴向液相速度径向分布的模拟结果。所有的时均曲线均呈抛物线形,中心点的液相速度随表观气速的提高而增加。从图中可以看出,液相速度的拐点发生在0.6<r/R<0.8范围内,这与文献结果相符,但拐点位置随表观气速增加而越来越靠近壁面,这是因为气速增大时,气泡的摆动更为剧烈。另外,液相速度最低点随表观气速变化并不明显,这可能是因壁面影响造成的。图10(a)、(b)给出了液相与气相法向雷诺应力分量¯u′zu′z的径向分布结果。从图可以看出,气/液相的轴向脉动速度分量最大值不是在塔中心,而是在r/R=0.6附近,即液相速度分布的拐点位置,说明该处速度变化剧烈,湍动强度大。比较图10(a)、(b)可以看出,气相的平均脉动速度要略大于液相,这与周力行的结论一致。由于缺乏可用实验数据,本文只对模拟结果进行定性分析。Degaleesan通过实验获得162mm鼓泡塔内表观气速Ug=0.3m·s-1时的液相轴向雷诺应力分量,其最大值约为3500cm2·s-2,比本文模拟结果大1倍(忽略塔径的影响),特别是塔中心,模拟结果要比实验值低一个数量级。Degaleesan等还测得了140mm鼓泡塔内表观气速Ug=0.12m·s-1时的液体轴向雷诺应力分量,若不考虑塔径效应,本文模拟结果仅为实验结果的1/3,但曲线的形状与实验曲线比较接近。总体说来,模拟的雷诺应力结果与实验值相差较大。Chen等采用相同湍流模型得到的雷诺应力模拟结果也同样偏低,而Zhou等采用二阶矩代数应力模型则较好地预测了气/液相雷诺应力。可见本文所采用的湍流模型不能有效地模拟气液相湍流。由于k-ε湍流模型为雷诺时均模型,时均化将使许多瞬态的高频流动特性被过滤掉,模拟得到的仅为拟瞬态结果,因此,采用拟瞬态结果来计算雷诺应力将会偏低。尽管本文所采用的湍流模型低估了湍流脉动速度,但该模型却较为合理地预测了液相湍流动能,图11给出了表观气速0.12m·s-1时的湍流动能模拟值,并与实验值进行了比较,预测结果与实验值较为吻合。3.3表观气速对模拟结果的影响体积传质系数是评估反应器性能最为重要的参数,对于快反应,体积传质系数(Kla)决定了反应器的生产能力及生产效率。Kla是传质系数和相间比表面积的乘积。关于传质系数Kl有不同的理论,如双膜理论、湍流涡理论等,但各理论之间存在相互矛盾的地方,因此,本文将着重考察气液界面表面积密度与表观气速的关系。关于界面密度已有报道,本文的模拟结果将与文献进行定量对比(图12),由图12可以看出,表观气速较低时(<0.10m·s-1)时,模拟结果与文献比较接近,而