李 鐵，宋濟(jì)洋，呂昌堯，楊志瑞

(東北電力大學(xué) 能源與動力工程學(xué)院，吉林 吉林 132012)

冷卻管內(nèi)氣液界面流動特性數(shù)值模擬研究

李 鐵，宋濟(jì)洋，呂昌堯，楊志瑞

(東北電力大學(xué) 能源與動力工程學(xué)院，吉林 吉林 132012)

采用數(shù)值模擬的方法，并應(yīng)用VOF多相流模型，對管內(nèi)液膜的速度分布特性、管內(nèi)動壓力分布特性進(jìn)行了分析，探討了氣液兩相速度差與氣液兩相湍動能大小的關(guān)系及不同液膜入口厚度下的液膜流型。結(jié)果表明：液膜平均速度在軸向位置先增大，后趨于穩(wěn)定；液體入口流量對液膜軸向速度分布影響顯著；液膜沿管徑方向速度梯度較大；管內(nèi)動壓力受液體入口速度的影響明顯；氣液兩相速度差越大，湍動能越大；液膜入口厚度越厚，液膜越穩(wěn)定。

液膜；氣液界面；VOF模型；數(shù)值模擬

冷卻管是德士古水煤漿氣化爐激冷室的重要組成部分，是氣液并流的通道。高溫合成氣與液膜接觸時發(fā)生劇烈熱質(zhì)傳遞。若水膜斷裂，一方面降低了液膜的冷卻作用；另一方面使得冷卻管內(nèi)壁直接與高溫氣體直接接觸而燒壞，因此對冷卻管內(nèi)壁液膜流動特性與斷裂的研究具有重要意義[1]。引起液膜斷裂的因素包括影響流動的因素和影響熱質(zhì)傳遞的因素。

許杰等[2]人通過實驗對旋流式降膜流動規(guī)律進(jìn)行了總結(jié)：在相同入口流量下，管的上半部分液膜厚度的波動較大，且在豎直高度500 mm處不再變化。李鐵等[3]人研究表明，液體入口流量越大，在壁面形成的液膜均勻性越好，越不容易斷裂。顏留成等[4]人對洗滌冷卻環(huán)流體冷態(tài)流動建立數(shù)學(xué)模型，通過對冷卻環(huán)出口處的流動狀態(tài)進(jìn)行試驗研究驗證了數(shù)學(xué)模型的可行性。結(jié)果表明：入水口周向位置的射流孔出口平均速度最小，兩入水口之間周向位置的射流孔出口平均流速最大。

因此，本文主要分析了影響流動的因素及液膜流動過程中管內(nèi)主要參數(shù)的變化規(guī)律。

1 物理模型

冷卻管內(nèi)實際工作環(huán)境很復(fù)雜，創(chuàng)建物理模型時需進(jìn)行一定簡化。由于合成氣中灰份濃度比較低，所以忽略氣相中灰份對合成氣流場的影響，并將合成氣體看做空氣。本文依據(jù)華東理工大學(xué)開發(fā)的新型洗滌冷卻室的實驗裝置，對其洗滌冷卻管進(jìn)行1∶1的比例建模[5]。圖1給出了氣液兩相管內(nèi)并行流動示意圖，在冷卻管內(nèi)氣液兩相并行向下流動，合成氣在管中心流動，激冷水沿管內(nèi)壁呈膜狀流動。表1給出了數(shù)值計算參數(shù)，氣液兩相的入口速度方向均與重力方向一致，并規(guī)定軸向1.5 m處為兩相的入口端。

圖1 管內(nèi)氣液并行流動示意圖

表1 冷卻管尺寸及計算參數(shù)

參數(shù)數(shù)值參數(shù)數(shù)值冷卻管直徑d/m0.15冷卻管長度h/m1.5液膜入口厚度/m0.005操作壓力p/kPa101.325氣體密度ρ/(kg·m-3)1.225氣體入口速度vg/(m·s-1)1-6液體入口流量Q/(m3·h-1)4-9.8

2 數(shù)學(xué)模型

2.1 VOF模型

因為兩相流動過程中沒有發(fā)生相互滲透，所以本文采用VOF模型[6]。Fluent軟件應(yīng)用有限體積法來離散輸運方程，氣液界面運動的追蹤是基于某相的分布，在每個控制容積內(nèi)，所有相體積分率的和為1，所有變量及其屬性的區(qū)域被各相共享并且代表了容積平均值，因此每一相的容積比率在每一位置是可知的。這樣，在任何給定單元內(nèi)的變量及其屬性或者純粹代表了一相，或代表了相的混合，這取決于容積比率值，即：

αg=0(純氣相) ，(1)

αg=1(純液相) ，(2)

0<αg<1(氣液混合相) ，(3)

在VOF方法中出現(xiàn)了許多界面重構(gòu)的方法，本文采用的是幾何重構(gòu)技術(shù)。由已知的VOF函數(shù)在單個網(wǎng)格上構(gòu)造有斜率的直線來逼近運動界面，然后再由流體的輸運特性，構(gòu)造下一時刻的VOF函數(shù)值[7]。

2.2 動量方程

基于歐拉法，建立連續(xù)介質(zhì)流動的控制方程，動量方程為

式中： ρ為密度；u為速度；p為壓力；g為重力加速度；Fs為外力；μ為粘性常數(shù)。

采用不可壓縮湍流的標(biāo)準(zhǔn)k-ε模型

其中，Gk、G1ε和G2ε為常數(shù)。

2.3 表面張力與壁面粘附處理方法

采用連續(xù)表面張力模型(CFS)描述自由表面張力[8]。在VOF模型中，表面張力作為動量方程的源項。每個網(wǎng)格內(nèi)僅存兩相介質(zhì)，方程如下：

設(shè)流體與壁面的接觸角，即壁面處界面的切線和界面的夾角為θw，壁面的單位法向量和單位切向量分別為nw和Tw，則近壁面單元的表面法向量為

nb=nwcosθw+Twsinθw.(9)

3 網(wǎng)格的劃分與數(shù)值算法

網(wǎng)格劃分的重點在于液膜區(qū)域網(wǎng)格的處理。液膜厚度與管徑相差2個數(shù)量級，而且在后續(xù)的計算中發(fā)現(xiàn)，液膜厚度變化范圍很大，所以液膜區(qū)域的網(wǎng)格要足夠密才能保證計算的精度。采用結(jié)構(gòu)化網(wǎng)格，氣體主體部分網(wǎng)格邊長為2 mm。對液膜區(qū)域進(jìn)行了局部加密處理，網(wǎng)格邊長為1 mm。采用有限體積隱式迭代，對體積分?jǐn)?shù)方程、動量方程分離求解。對壓力項離散應(yīng)用Body Force Weighted方法；對壓力速度耦合采用PISO算法，該算法以SIMPLE算法為基礎(chǔ)，做了2個附加修正：臨近修正[9]和偏斜修正[10]。為了計算簡便，本文沿冷卻管對稱中心選取1/2對稱面進(jìn)行計算[11]。

4 數(shù)值計算結(jié)果與分析

4.1 液膜流動過程中的形態(tài)變化

圖2給出了常壓下vg=2.5 m/s，Q=4 m3·h-1時，不同時刻不同高度區(qū)間內(nèi)的液膜流動形態(tài)。

圖2 vg=2.5 m/s，Q=4 m3·h-1時不同時刻不同高度區(qū)間內(nèi)液膜流動形態(tài)示意圖

圖3給出了常壓下vg=2.5 m/s，Q=8.2 m3·h-1時，不同時刻不同高度區(qū)間內(nèi)的液膜流動形態(tài)。

圖3 vg=2.5 m/s，Q=8.2 m3·h-1時不同時刻不同高度區(qū)間內(nèi)液膜流動形態(tài)示意圖

由于液膜非常薄，需截取不同高度上的局部位置進(jìn)行觀察。液膜沿管壁流動時流動形態(tài)發(fā)生多次變化。t=0.4 s，液膜出現(xiàn)斷裂帶。t=0.8 s，有部分液滴脫離管壁飛濺到冷卻管中。t=1.0 s，粘附在壁面的液滴重新聚集成連續(xù)不間斷的液膜向下流動。在氣液并行流動過程中，液膜是一種破碎、聚集、再破碎、再聚集的循環(huán)流動過程。對比圖2和圖3可知，液體入口流量越大，液膜越厚，發(fā)生斷裂的位置減少。

4.2 不同液膜入口厚度下的流型

圖4給出了vg=2.5 m/s，Q=8.2 m3·h-1時，不同液膜入口厚度下管內(nèi)入口區(qū)域的液膜流型。

圖4 不同液膜入口厚度下流型

圖5給出了常壓下vg=2.5 m/s，Q=8.2 m3·h-1時，不同液膜入口厚度下液膜波動曲線。

結(jié)合液膜流型圖和流動過程中液膜厚度波動曲線可知，不同液膜入口厚度對液膜流型影響顯著。整體上看，液膜在沿管壁向下流動的過程中，逐漸變薄，最后趨于穩(wěn)定，不同液膜入口厚度下，液膜在流動過程中的變化趨勢基本一致，但液膜入口厚度的增大使得液膜流動過程中的整體厚度變厚。從波動曲線圖可以看出，當(dāng)液膜入口厚度較小時，流動過程中液膜波動頻率較大，一直持續(xù)到管的出口，液膜厚度不穩(wěn)定。隨著液膜入口厚度的增加，波動減弱?？梢?，液膜入口厚度越大，液膜越穩(wěn)定。

圖5 不同液膜入口厚度下的液膜波動

4.3 液膜速度分布

圖6給出了氣體入口速度vg=2.5 m/s，液體流量Q=8.2 m3·h-1不變時，兩相速度分布。

圖6 兩相速度分布

圖7給出了常壓下氣體入口速度vg=2.5 m/s時，不同液體入口流量下液膜軸向平均速度分布；圖8給出了相同工況下的液膜徑向速度分布。其中,0 mm處為液膜與壁面交界處。

圖7 軸向平均速度分布

圖8 徑向速度分布

由圖7可知，液膜平均速度大小與液體入口流量成正比。液體流量從4 m3·h-1增大到8.2 m3·h-1時，液膜速度增幅程度較大。繼續(xù)增大液體流量到Q=9.8 m3·h-1，液膜速度增量變小。液膜在重力和壁面剪切力的同時作用下加速運動一段時間后，流速趨于穩(wěn)定。在液膜入口區(qū)域(1.4 m-1.5 m)，液體入口流量越小，曲線斜率越大，即加速度越大。這是因為液體入口流量較小時，液膜的初速度小，壁面對液膜的剪切力小，所以液膜加速度較大。當(dāng)流量變大，液膜的初速度變大，壁面剪切力變大，加速度變小。液膜加速度越大，增速越快，液膜受力到達(dá)平衡時間越短。故液體流量較小時，液膜速度率先達(dá)到穩(wěn)定。

由圖8可知，由于壁面存在粘性底層，所以與壁面接觸處的液膜速度為0。液膜靠近壁面處有較大的徑向速度梯度，速度增大幅度明顯，隨著膜厚的增加，速度增大趨于平緩。液體入口流量越大，速度梯度越大，液膜整體的徑向速度越大。液體入口流量不同情況下，液膜徑向速度有著相似的分布規(guī)律。這是由于液膜內(nèi)部各質(zhì)點粘性力相互作用的結(jié)果。

4.4 管內(nèi)動壓力分布

圖9給出了氣體入口速度vg=2.5 m/s，液體流量Q=8.2 m3·h-1時，兩相動壓力分布。

圖9 管內(nèi)兩相動壓力分布

圖10(a)、圖10(b)分別給出了氣體入口速度vg=2.5 m/s，液體流量Q=8.2 m3·h-1時，兩相動壓力分布。

由圖10(a)、圖10(b)可知，管內(nèi)氣體的動壓力沿軸向方向逐漸變小，液膜的動壓力沿軸向方向逐漸增大，最后兩者趨于穩(wěn)定。由動壓力的計算公式P=1/2ρv2可知，管內(nèi)動壓力大小取決于兩相的密度和速度。流動過程中，兩者的密度大小不變，所以氣液動壓力的大小完全取決于兩者的速度變化規(guī)律。氣體在流動時速度略有減小，所以動壓力沿軸向只降低了0.35 Pa，減小幅度不大。液膜速度沿軸向先增加后不變，動壓力增加了1 381 Pa，增幅程度較為明顯。液膜動壓力分布曲線與液膜軸向速度分布曲線呈相同的變化趨勢。

圖10 管內(nèi)氣液兩相動壓力沿軸向分布曲線

圖11(a)給出了氣體入口速度vg=2.5 m/s不變時，不同液速下管內(nèi)動壓力軸向分布。圖11(b)給出了液體入口流量Q=8.2 m3·h-1不變時，不同氣速下管內(nèi)動壓力軸向分布。

圖11 氣液兩相速度變化對管內(nèi)動壓力分布的影響

對比圖11(a)、圖11(b)兩種情況，當(dāng)分別增大氣、液兩相入口速度時，管內(nèi)的動壓力都有所增加，但改變液相入口速度時管內(nèi)動壓力的增量顯著，這是因為液體的密度要遠(yuǎn)遠(yuǎn)大于氣體密度。管內(nèi)動壓力的分布曲線與液相動壓力的分布曲線基本一致，這進(jìn)一步說明了液相動壓力的改變對管內(nèi)動壓力的分布起主導(dǎo)作用。

4.5 液膜流動過程中管內(nèi)湍動能分布特性

湍動能是衡量湍流發(fā)展和衰退的指標(biāo)，氣液間湍流強度的大小是影響氣液間熱質(zhì)傳遞的重要因素之一。氣液間發(fā)生湍流流動時會產(chǎn)生大量的漩渦，以漩渦擴(kuò)散為主，同時會增大傳質(zhì)質(zhì)量。圖12給出了常壓下vg=2.5 m/s，Q=8.2 m3·h-1不變時，液膜向下流動過程中管內(nèi)氣液兩相湍動能分布。

氣液相湍動能較大的部分主要分布在氣液交界面處，而氣體中心的湍動能很小。在進(jìn)口端，液膜在速度較大的氣體擾動下，發(fā)生湍流流動，氣液間發(fā)生激烈的能量傳遞，此時的湍動能較大。隨著液膜流動的進(jìn)行，湍動區(qū)域變長，當(dāng)t=0.3 s時，液膜發(fā)生堆積，液膜厚度變厚，氣體流動方向與降膜方向不再平行，氣體的流動受到液膜的阻擋，擾動加大，導(dǎo)致液膜局部湍動能變大。當(dāng)t=0.5 s時，液膜發(fā)生斷流在壁面上形成液滴，氣體掠過液滴而產(chǎn)生的局部漩渦使局部的湍動能增大幅度更加明顯。

圖12 不同時刻管內(nèi)湍動能分布

圖13給出了常壓下vg=2.5 m/s，Q=8.2 m3·h-1不變時，不同軸向位置截面上平均湍動能分布。圖14給出了相同工況下，不同軸向位置上兩相的速度差分布，其中液相速度取液膜的平均速度。

圖13 不同軸向位置截面上平均湍動能分布

圖14 不同軸向位置氣液兩相速度差分布

對比兩圖可知，不同軸向位置截面上的平均湍動能與氣液兩相的速度差有著相似的變化規(guī)律。當(dāng)兩相速度差較大時，管內(nèi)湍動能較大，隨著兩相速度差逐漸縮小到穩(wěn)定不變，管內(nèi)不同截面上的平均湍動能沿軸向位置同時逐漸變小，最后趨于穩(wěn)定。由分析結(jié)果可知，兩相速度差對管內(nèi)的湍動能分布有著重要的影響。

5 結(jié) 論

本文采用數(shù)值模擬的方法，對管內(nèi)液膜的速度分布特性、管內(nèi)動壓力分布特性進(jìn)行了分析，探討了氣液兩相速度差與氣液兩相湍動能大小的關(guān)系，不同液膜入口厚度下的液膜流型。得出以下結(jié)論：

(1)液體入口厚度越大，液膜越穩(wěn)定；

(2)液膜軸向平均速度在軸向位置先逐漸增大，后趨于穩(wěn)定；

(3)初始液體入口流量越小，加速度越大。液膜靠近壁面處有較大的徑向速度梯度，速度增大幅度明顯，隨著膜厚的增加，速度增大趨于平緩；

(4)分別增大氣、液兩相入口速度時，管內(nèi)的動壓力都有所增加，但改變液相入口速度時管內(nèi)動壓力的增量顯著。管內(nèi)動壓力的分布曲線與液相動壓力的分布曲線基本一致，這進(jìn)一步說明了液相動壓力的改變對管內(nèi)動壓力的分布起主導(dǎo)作用；

(5)不同軸向位置截面上的平均湍動能與氣液兩相的速度差有著相似的變化規(guī)律，兩相速度差越大，湍流流動產(chǎn)生的湍動能越大，兩相速度差對管內(nèi)的湍動能分布有著重要的影響。

[1] 王輔臣，于廣鎖，龔欣，等.大型煤氣化技術(shù)的研究與發(fā)展[J].化工進(jìn)展.2009，28(2)：173-180.

[2] 許杰.洗滌冷卻室內(nèi)流動及熱質(zhì)傳遞過程研究[D].上海：華東理工大學(xué)，2003.

[3] 李鐵.下降管內(nèi)壁激冷水降膜流動特性[J].東南大學(xué)學(xué)報，2006，36(6)：962-966.

[4] 顏留成，王亦飛，王靈萍，等.洗滌冷卻環(huán)內(nèi)液體冷態(tài)流動行為數(shù)值模擬[J].化工學(xué)報，2014，65(6)：2928-2933.

[5] 王靈萍，王亦飛，郭強強，等.洗滌冷卻管內(nèi)垂直降膜流動特性[J].化工學(xué)報，2013，64(6)：1959-1969.

[6] 劉儒勛，王志峰.數(shù)值模擬方法和運動界面追蹤[M].合肥：中國科學(xué)技術(shù)大學(xué)出版社，2001.

[7] D.L.Youngs.Time-dependent multi-material flow with large fluid distortion[M].New York：Academic Press，1982.

[8] J.U.ABrackbill.Continuum method formodeling surfacetension[J].Journal of Computational Physics，1992，100(3)：335-354.

[9] R.I.Issa.Solution of implicitly discretized fluid flow equations by operatorsplitting [J].Jowrnal of Computional Physics，1986，62(2)：40-65.

[10] J.Ferziger,L.PericM.Computationalmethods for fluiddynamics[M].Heidelberg：Springer-Verlag，1996.

[11] 李鐵，袁竹林.德士古氣化爐激冷室下降管內(nèi)氣液傳熱傳質(zhì)過程研究[J].東北電力大學(xué)學(xué)報，2013，33(3)：48-51.

Numerical Simulation of Gas-liquid Interface Flow Characferistics in Scrubbing-cooling Pipe

LI Tie，SONG Ji-yang，LV Chang-yao，YANG Zhi-rui

(Energy Resource and Power Engineering College，Northeast Electric Power University，Jilin Jilin 132012)

The distribution of the liquid film velocity and the dynamic pressure in the vertical pipe was simulated by using VOF multiphase flow model.The relationship between the velocity of gas-liquid two-phase and turbulent kinetic energy was discussed in the same method.Flow pattern under different inlet thickness of liquid film was discussed on the same time.The results showed that：average speed of the liquid film was increased in axial position and then stabilized；Axial velocity distribution was influenced significantly by liquid inlet flux；Liquid film velocity gradient along the direction of the diameter was larger；Dynamic pressure in the vertical pipe was influenced by liquid inlet velocity obvious；when the different value of gas-liquid two-phase velocity was bigger，the turbulent kinetic energy was larger.The thicker of the liquid film inlet thickness，the more stable of the liquid film.

Liquid film；Gas-liquid interface；Volume of fluid；Numerical simulation

2016-11-12

李 鐵(1979-)，女，博士，副教授，主要研究方向：多相流及數(shù)值模擬.

1005-2992(2017)01-0087-08

TQ021

A

電子郵箱： litie0622@163.com(李鐵)；707559438@qq.com(宋濟(jì)洋)；493566668@qq.com(呂昌堯)；445353083@qq.com(楊志瑞)