張井志，李蔚

?

毛細管內(nèi)氣液Taylor流動的氣泡及阻力特性

張井志1,2，李蔚1

（1浙江大學能源工程學院，浙江杭州310027；2浙江大學能源工程學院，先進航空發(fā)動機協(xié)同創(chuàng)新中心，浙江杭州 310027）

采用相對坐標系方法，研究毛細管（2mm）內(nèi)充分發(fā)展垂直上升氣液Taylor流動，分析兩種工作介質(zhì)下Taylor氣泡的形狀、上升速度、液膜厚度以及壓降特性。結(jié)果表明：隨著兩相表觀速度（tp）增大，Taylor氣泡長度增大，氣泡尾部曲率半徑增大。氣泡長度及內(nèi)部回流區(qū)隨著氣泡體積分數(shù)（g）增大而增大，量綱1液膜厚度與氣泡上升速度與毛細數(shù)（）正相關，模擬結(jié)果與經(jīng)驗公式吻合較好。摩擦阻力因子（c）隨tp與g的增大而降低，N2/乙二醇為工質(zhì)的Taylor流動c低于單相情況，而N2/水為工質(zhì)的Taylor流動c高于單相情況。Kreutzer等的流型依賴公式以及Lockhart等的分離模型可較好預測本文的兩相壓降，模擬結(jié)果與預測值的誤差在±10%以內(nèi)，常規(guī)通道所推薦5仍然適用于本文毛細管情況。

Taylor氣泡；氣液兩相流；數(shù)值模擬；毛細管；相對坐標系

引 言

Taylor流動是一種基本的氣液兩相流流型，廣泛應用于化工以及電子芯片散熱行業(yè)。隨著制造技術的發(fā)展以及電子器件、化工機械的微小化，微通道內(nèi)的兩相流動引起各國學者的關注。微通道具有較小當量直徑，流動形式主要為層流，很難采用常規(guī)通道中利用湍流擾動的方式進一步強化換熱。而微通道內(nèi)兩相流動，由于氣液兩相界面的存在，使液側(cè)流體產(chǎn)生擾動強化傳熱傳質(zhì)過程，微通道內(nèi)穩(wěn)定周期性Taylor流動可以滿足化工以及芯片散熱要求。

早期Taylor流動的研究主要采用實驗與理論方法，主要針對通道壓降、兩相流流型圖以及液膜厚度等方面分析。Kreutzer等[1]測定了小通道內(nèi)Taylor流動的壓降情況，結(jié)果表明壓降主要與/有關，并提出了一個依賴流型的微通道Taylor流動的壓降公式。Liu等[2]實驗研究了豎直通道內(nèi)的Taylor流動情況，并提出了氣泡上升速度以及壓降實驗關系式。Han等[3]利用激光位移計測量了液膜厚度，指出液膜厚度主要依賴毛細數(shù)（）Reynolds數(shù)（）。

隨著計算性能的提升以及兩相數(shù)值研究方法的發(fā)展，數(shù)值模擬逐漸成為研究毛細管內(nèi)Taylor流動的重要手段。相比實驗方法，數(shù)值計算可以獲得精確的內(nèi)部流動情況。黨敏輝等[4]的模擬結(jié)果表明入口形式不同，氣泡長度受入口速度的影響不同。Qian等[5]研究了T型微通道內(nèi)的Taylor流動，模擬結(jié)果與實驗結(jié)果吻合較好，給出了氣柱與液柱長度的公式。Shao等[6]研究了入口段對于氣泡大小的影響，模擬結(jié)果表明，氣泡大小隨噴嘴直徑以及壁厚增大而增大。賀瀟等[7]對垂直以及傾斜上升管內(nèi)氣液兩相彈狀流進行研究，指出壁面切應力隨著Froude數(shù)增大而增大。Taha等[8]、Zheng等[9]、Araújo等[10]均采用兩相相對坐標系的方法研究豎直上升通道內(nèi)的Taylor流動。Asadolahi等[11]利用動網(wǎng)格技術，追蹤單個Taylor氣泡的變化，分析了充分發(fā)展狀態(tài)下氣泡的換熱和壓降情況。

對于微通道內(nèi)的Taylor流動已有較多的實驗模擬研究，但充分發(fā)展狀態(tài)下流動分析工作仍然不多。本文主要利用相對坐標系的方法，研究管徑2 mm毛細管內(nèi)，工作介質(zhì)為N2/水和N2/乙二醇充分發(fā)展狀態(tài)氣液Taylor流動，分析了影響氣泡上升的速度、液膜厚度及壓降的因素，并與經(jīng)驗公式進行對比。

1 數(shù)值計算模型

氣液兩相的數(shù)值模擬采用Ansys Fluent 12.0進行，采用VOF模型捕捉氣液界面。VOF模型基于兩相無法互相滲透的假設，在連續(xù)性、動量方程的基礎上，增加體積分數(shù)方程，求解計算域內(nèi)部的氣液分布情況。VOF模型質(zhì)量守恒特性較好，并具有較高的界面精度，廣泛應用于氣液兩相流數(shù)值模擬。

1.1 控制方程

計算區(qū)域流體控制方程如下：

連續(xù)性方程

動量方程

(2)

體積分數(shù)方程

其中

1.2 幾何模型及計算方法

本文采用二維軸對稱模型分析毛細管內(nèi)Taylor流動，計算域如圖1所示，長度設定為20 mm（10）。入口采用速度入口邊界條件，速度初始分布為管內(nèi)層流充分發(fā)展的速度分布，即拋物線形分布。出口采用壓力出口邊界條件，出口壓力設為標準大氣壓。上壁面設為無滑移的壁面，壁面接觸角設置為90°，采用連續(xù)表面張力模型（CSF）考慮表面張力作用。對稱軸采用軸邊界，重力方向沿軸負方向。

圖1 計算模型

將全場的軸向速度設定為入口速度值，利用用戶自定義函數(shù)(UDF)將出口的速度分量賦予入口。假設Taylor氣泡初始為圓柱形，放置在計算域的中心，其長度主要由初始氣體體積分數(shù)（g）決定。計算結(jié)果表明，最終結(jié)果與初始氣泡的形狀無關，初始形狀僅影響計算時間。采用相對坐標系方法，整個計算域的運動坐標系固定于Taylor氣泡上，當氣泡的移動速度趨于穩(wěn)定值時得到穩(wěn)定的流場以及氣液分界面。氣泡速度的定義如式(7)所示

1.3 網(wǎng)格無關性

計算域采用結(jié)構(gòu)化網(wǎng)格進行劃分，近壁面附近漸進加密，網(wǎng)格的樣式見圖1。對網(wǎng)格進行無關性驗證，如圖2所示，3種網(wǎng)格下Taylor流動氣泡形狀基本一致，其中粗網(wǎng)格的氣泡長度稍高于加密后的情況，兩種加密方法氣泡形狀及長度基本一致。由于本文采用變時間步長，為節(jié)約計算時間，本文模擬采用的網(wǎng)格數(shù)目為95×1540（徑向網(wǎng)格數(shù)目×軸向網(wǎng)格數(shù)目）。

圖2 無關性檢測

2 結(jié)果分析與討論

2.1 氣泡形狀與流線

圖3為不同流速下，兩種工質(zhì)Taylor氣泡的形狀以及流線圖，其中藍色代表氣泡區(qū)域，紅色代表液體區(qū)域。Taylor氣泡主要由近似球形的頭部、尾部以及圓柱形的氣柱組成，在氣泡與壁面之間存在一層薄液膜，此部分區(qū)域的液體相對氣泡的速度向下。工質(zhì)為N2/水的Taylor流動，在氣泡內(nèi)部存在3個比較明顯的旋渦區(qū)域，而為N2/乙二醇時僅能觀測到1個較大的旋渦區(qū)域，且液膜厚度與氣泡長度均大于前者。Taylor氣泡外部存在兩個較大的旋渦，增強核心區(qū)域與壁面區(qū)域間傳熱、傳質(zhì)過程，強化傳熱、傳質(zhì)效果。隨著入口速度（tp）增大，兩種工質(zhì)的Taylor氣泡的頭部均拉長，尾部曲率半徑增大，趨于扁平，而后者的變化趨勢更加劇烈。

圖3 N2/H2O、N2/(CH2OH)2為工質(zhì)時Taylor流動流線圖及氣泡形狀

圖4 兩種工作介質(zhì)Taylor氣泡輪廓圖

由于具有較高的黏性系數(shù)，工質(zhì)為N2/乙二醇的Taylor氣泡相對穩(wěn)定，尾部無不穩(wěn)定區(qū)域。隨著流速的增大，Taylor氣泡尾部逐漸趨于平坦，氣泡形狀更趨近子彈形。工質(zhì)為N2/水時，表面張力要高于工質(zhì)為N2/乙二醇，Taylor氣泡更容易保持半球形的尾部。

tp0.3662 m·s-1，不同g時Taylor氣泡形狀如圖4 (b)所示，氣泡長度隨著g的增大基本呈線性增大。在同樣的tp與g下工質(zhì)為N2/乙二醇的Taylor氣泡長度要高于工質(zhì)為N2/水的情況。整個Taylor氣泡形狀基本不隨g改變而變化，氣泡頭部、尾部形狀，尾部不穩(wěn)定區(qū)域，以及頭部到穩(wěn)定液膜區(qū)域的長度保持一致，決定氣泡形狀的主要因素為入口的兩相表觀流速與工質(zhì)物性。

2.2 氣泡的移動速度與液膜厚度

Taylor氣泡的上升速度b高于tp，氣泡相對與液體向上運動，導致了薄液膜區(qū)域的形成。Liu等[2]實驗研究了豎直上升毛細管內(nèi)的Taylor流動，指出b/tp主要由決定，并提出了b/tp的表達[式(8)]，適用范圍0.0002＜＜0.39。Liu等[2]的實驗結(jié)果以及本文的模擬結(jié)果如圖5所示，由圖可得，式(8)可以較好地預測本文模擬結(jié)果，模擬結(jié)果與公式預測值的誤差在-5%～+10%。b/tp隨著增大而增大，增大趨勢逐漸降低。在同樣的tp、g下，Taylor氣泡隨著增大變得細長，更多氣泡處于核心區(qū)域。由充分發(fā)展的層流管內(nèi)流動速度分布可得，速度由核心區(qū)域至壁面呈現(xiàn)拋物線形遞減關系。高下，Taylor氣泡占據(jù)更多的高流速區(qū)域是導致tp/b上升的一個原因。

Taylor氣泡頭部為近似圓球或橢球形，液膜厚度沿著頭部到尾部方向先降低后趨于穩(wěn)定，在尾部出現(xiàn)不穩(wěn)定的波動，以穩(wěn)定的液膜處的厚度作為Taylor流動的液膜厚度。Irandoust等[12]實驗測定了Taylor流動的液膜厚度，并提出了液膜厚度的預測公式[式(9)]。Aussillous等[13]利用實驗數(shù)據(jù)分析了液膜厚度的影響因素，指出在慣性力可以忽略的條件下，式(10)可以較好地預測液膜厚度。

圖6為本文模擬結(jié)果，Leung等[14]以及Taylor[15]的實驗結(jié)果給出了量綱1液膜隨的變化規(guī)律，模擬結(jié)果與文獻實驗值吻合較好。式(9)可以較好地預測N2/水的Taylor流動的液膜厚度，模擬值與預測值的誤差在±20%以內(nèi)。N2/乙二醇的Taylor流動，慣性力的作用相對較小，其液膜厚度更符合式(10)，當＞0.15后，模擬值與預測值的誤差在10%以內(nèi)。Leung等[14]實驗結(jié)果同樣表明，式(9)可以更好地預測N2/水為工質(zhì)的Taylor液膜厚度。在固定的g條件下，氣泡長度隨著的增大而增大，氣泡變得瘦長，液膜區(qū)域流過的液體增多，增大。

2.3 壓降特性

計算域的摩擦阻力系數(shù)定義如下

Taylor流動的計算域阻力壓降主要由液柱壓降與維持氣泡的拉普拉斯壓降組成。由于Taylor氣泡的存在，液柱所占據(jù)的區(qū)域要小于計算域的長度。由圖7 (a)可得，與單相流動類似，c隨著tp增大而降低，降低趨勢逐漸減小。由于乙二醇的黏度系數(shù)高于水，以N2/乙二醇為工質(zhì)Taylor流動的c遠大于N2/水的情況。c（N2/乙二醇）＜16，表明以N2/乙二醇為工質(zhì)Taylor流動的壓降要低于單相流動，氣泡的拉普拉斯壓降的作用不足以抵消液柱長度縮短帶來的阻力壓降降低。c（N2/水）＞16，Taylor流動壓降高于單相情況，在同樣的流速情況下，Taylor氣泡對計算域壓降的影響要高于N2/乙二醇的情況。

由圖7 (b)可得，c隨著g的增大而降低，基本呈線性變化的規(guī)律，主要是由于彈狀流氣泡的長度隨著g的增大而增大，液柱長度逐漸降低。由于氣泡的頭部及尾部形狀不隨g的變化而改變，拉普拉斯壓降保持不變，整場的阻力系數(shù)隨著液柱長度的降低而降低。

圖7 fcRe隨Vtp 和ξg 的變化規(guī)律

Kreutzer等[1]利用數(shù)值模擬與實驗，分析了慣性力以及表面張力對Taylor流動的影響，認為slug主要與slug/(/)0.33有關，壓降關系如式(13)所示。由于實驗中工質(zhì)不純而引起的Marangoni效應，對于模擬結(jié)果，對于實驗結(jié)果。Warnier等[16]指出，當＞150時，更能精確預測實驗結(jié)果；當＜150時，可以更好地吻合實驗結(jié)果。

由圖8可得，式(13)可以很好地預測本文的模擬結(jié)果以及Walsh等[17]、Kreutzer等[1]的結(jié)果。在低時，預測精度較高，相反具有較好的預測效果。由于Taylor氣泡的作用，slug＞16，當液柱長度趨于無窮大時，slug16，與單相情況一致。對于本文模擬情況，slug隨著的增大而降低，slug隨著的降低而降低。

圖8 fslugRe隨Lslug/d(Ca/Re)0.33的變化規(guī)律

由式(13)可以推導出全場的壓強梯度表達式

其中，對于工質(zhì)為N2/水，，工質(zhì)為N2/乙二醇，。由圖9可得，式(14)可以很好地預測本文的模擬壓降，誤差均在±10%以內(nèi)，式(14)的正確應用依賴于液柱長度的測量精度。在工程應用中，由于條件的限制很難測量液柱的長度。兩相壓力梯度同樣可以由Lockhart等[18]提出的分離模型（LMC模型）進行計算，此種方法的優(yōu)點在于方便使用，各個參數(shù)均由入口條件決定。LMC模型的描述如式(15)～式(17)所示，對于氣、液兩相均為層流，5。

由圖9可以看出，LMC模型同樣具有較高的精度，模擬值與預測值的誤差在±10%以內(nèi)。從本質(zhì)上講，兩種壓降公式均是在單相流體壓降梯度的基礎上考慮第二相的作用，進行修正得到全場的壓降。Kreutzer等[1]的公式更依賴于流型，可以從本質(zhì)上解釋壓降變化的原因，而LMC模型的適用性更廣，不再局限于周期性Taylor流動。5為常規(guī)通道的推薦值，但從本文模擬結(jié)果可以看出，其也可較好地預測微小通道內(nèi)的壓強梯度。

3 結(jié) 論

（1）Taylor氣泡長度隨著tp與g的增大而增大，由于具有較小的表面張力系數(shù)與較大的黏度，N2/乙二醇的Taylor流動尾部曲率變化較大，同時尾部不會出現(xiàn)不穩(wěn)定區(qū)域。

（2）b/tp與/隨著的增大而增大，Taylor氣泡隨著的增大占據(jù)更多的核心區(qū)域，導致了b上升、增大，文獻經(jīng)驗公式可以較好地預測本文的模擬結(jié)果。

（3）c隨著tp與g的增大而降低，N2/乙二醇的Taylor流動c要低于單相情況，而N2/水為工質(zhì)的Taylor流動c高于單相情況。slug主要與量綱1的液柱長度以及物性有關，而與進口流速無明顯關系。

（4）Kreutzer等[1]的流型依賴公式以及Lockhart等[18]的兩相分離模型均可以較好地預測本文的模擬結(jié)果，模擬結(jié)果與預測值的誤差在±10%以內(nèi)，常規(guī)通道所推薦5仍然適用于本文模擬的毛細管情況。

（5）從本文模擬來看，對于直徑2 mm的圓管，其阻力特性與常規(guī)通道基本一致。相對常規(guī)通道，微小通道內(nèi)表面張力的作用逐漸占據(jù)主導地位，氣泡可以維持較為穩(wěn)定的形態(tài)。對毛細管內(nèi)Taylor氣泡的進一步工作需要考慮管徑的尺度效應，重力、黏性力、慣性力及表面張力的相對作用。

符 號 說 明

a——公式系數(shù) C——Chisholm系數(shù) Ca——毛細數(shù) d——直徑，m F——動量方程源項，N·m-3 f——摩擦阻力系數(shù) g——重力加速度，m·s-2 L——長度，m P——壓強，Pa R——半徑，m Re——Reynolds數(shù) r——徑向坐標，m t——時間，s V——速度，m·s-1 v——速度矢量，m·s-1 X——Lockhart-Martinelli參數(shù) x——軸向坐標，m α——單元格體積分數(shù) β——單相與兩相表觀速度之比 δ——液膜厚度 ξ——氣泡與計算域體積之比 ρ——密度，kg·m-3 Ф——Lockhart-Martinelli參數(shù) 下角標 b——氣泡 c——計算域 G——氣相 L——液相 slug——液柱 tp——兩相

References

[1] Kreutzer M T, Kapteijn F, Moulijn J A,. Inertial and interfacial effects on pressure drop of Taylor flow in capillaries [J]., 2005, 51 (9): 2428-2440

[2] Liu H, Vandu C O, Krishna R. Hydrodynamics of Taylor flow in vertical capillaries: flow regimes, bubble rise velocity, liquid slug length, and pressure drop [J]., 2005, 44 (14): 4884-4897

[3] Han Y, Shikazono N. Measurement of liquid film thickness in micro square channel [J]., 2009, 35 (10): 896-903

[4] Dang Minhui (黨敏輝), Ren Mingyue (任明月), Chen Guangwen (陳光文). Effect of microchannel inlet configuration on Taylor bubble formation in microreactors [J].(化工學報), 2014, 65 (3): 805-812

[5] Qian D, Lawal A. Numerical study on gas and liquid slugs for Taylor flow in a T-junction microchannel [J]., 2006, 61 (23): 7609-7625

[6] Shao N, Salman W, Gavriilidis A,. CFD simulations of the effect of inlet conditions on Taylor flow formation [J]., 2008, 29 (6): 1603-1611

[7] He Xiao (賀瀟), Che Defu (車得福). CFD simulation of wall shear stress in vertical and inclined upward slug gas-liquid flow [J].() (化工學報), 2008, 59 (6): 1390-1395

[8] Taha T, Cui Z. Hydrodynamics of slug flow inside capillaries [J]., 2004, 59 (6): 1181-1190

[9] Zheng D, He X, Che D. CFD simulations of hydrodynamic characteristics in a gas-liquid vertical upward slug flow [J]., 2007, 50 (21): 4151-4165

[10] Araújo J D P, Miranda J M, Campos J B L M. Flow of two consecutive Taylor bubbles through a vertical column of stagnant liquid—a CFD study about the influence of the leading bubble on the hydrodynamics of the trailing one [J]., 2013, 97: 16-33

[11] Asadolahi A N, Gupta R, Fletcher D F,. CFD approaches for the simulation of hydrodynamics and heat transfer in Taylor flow [J]., 2011, 66 (22): 5575-5584

[12] Irandoust S, Andersson B. Liquid-film in taylor flow through a capillary [J]., 1989, 28 (11): 1684-1688

[13] Aussillous P, Quéré D. Quick deposition of a fluid on the wall of a tube [J].(), 2000, 12 (10): 2367-2371

[14] Leung S S, Gupta R, Fletcher D F,. Effect of flow characteristics on Taylor flow heat transfer [J]., 2011, 51 (4): 2010-2020

[15] Taylor G I. Deposition of a viscous fluid on the wall of a tube [J]., 1961, 10 (2): 161-165

[16] Warnier M J F, de Croon M, Rebrov E V,. Pressure drop of gas-liquid Taylor flow in round micro-capillaries for low to intermediate Reynolds numbers [J]., 2010, 8 (1): 33-45

[17] Walsh E, Muzychka Y, Walsh P,. Pressure drop in two phase slug/bubble flows in mini scale capillaries [J]., 2009, 35 (10): 879-884

[18] Lockhart R, Martinelli R. Proposed correlation of data for isothermal two-phase, two-component flow in pipes [J]...., 1949, 45 (1): 39-48

Bubble and frictional characteristics of gas-liquid Taylor flow in capillary tube

ZHANG Jingzhi1,2, LI Wei1

(College of Energy EngineeringZhejiang UniversityHangzhouZhejiangChina;CoInnovation Center for Advanced AeroEngineCollege of Energy EngineeringZhejiang UniversityHangzhouZhejiangChina

In order to obtain the frictional characteristics of fully developed Taylor flow in the vertical capillary tube, numerical simulations of the flow in the capillary tube with diameter of 2 mm were conducted by using the moving frame reference method. The shape, rising velocity of Taylor bubble, liquid film thickness and pressure drop were obtained using two different working fluids and analyzed. Simulation results showed that the length of Taylor bubble and the radius of curvature increased with increasing two-phase superficial velocitytp. The length of Taylor bubble also increased with increasing gas voidg, while the nose and tail of Taylor bubble were independent ofg. Dimensionless thickness of liquid film and rising velocity of Taylor bubbles were proportional to capillary number. Friction factorcdecreased with increasingtpandg. Thecof Taylor flow with N2/ (CH2OH)2as workingfluid was lower than that of single phase with the sametp, while thecfor N2/H2O was higher than that of single phase. The model proposed by Lockhart and Martinelli, and the flow pattern dependent model proposed by Kreutzer. could predict the pressure drop obtained from simulation with an error of±10%. The Chisholm number5 which was recommended for conventional tube when both phases were laminar was also reasonable for the capillary tube in the simulation work.

Taylor bubble; gas-liquid flow; numerical simulation; capillary tubes; relative coordinate

2014-10-28.

Prof. LI Wei, weili96@zju.edu.cn

10.11949/j.issn.0438-1157.20141622

TK 124

A

0438—1157（2015）03—0942—07

浙江省自然科學基金項目（Z13E060001）；國家科技支撐計劃項目（2012BAA10B01）。

2014-10-28收到初稿，2014-12-11收到修改稿。

聯(lián)系人：李蔚。第一作者：張井志（1988—），男，博士研究生。

supported by the Natural Science Foundation of Zhejiang Province (Z13E060001) and the Chinese National Key Technology R&D Program (2012BAA10B01).