李興雨 李琳 譚義海

摘 要：為了研究虹吸管中氣泡運動行為與變形特征，對虹吸管中氣泡運動過程進行數(shù)值模擬，結(jié)果表明：上行管內(nèi)氣泡縱橫比越小，氣泡速度越大，流量的變化對氣泡縱橫比影響很小，且氣泡形狀與靜水中氣泡相圖相吻合，氣泡在上升過程中會發(fā)生橫向振蕩，振蕩幅度隨著氣泡直徑的增大先增大后減小;下行管內(nèi)流量不同，氣泡運動方向也不同，氣泡變形情況與上行管內(nèi)一致;中行管內(nèi)氣泡為貼壁狀態(tài)流動，隨著氣泡直徑增大，氣泡速度減小，當氣泡直徑增加到一定數(shù)值時，氣泡會發(fā)生破裂，出現(xiàn)較小的氣泡，氣泡速度又出現(xiàn)增大的現(xiàn)象。

關(guān)鍵詞：虹吸管;數(shù)值模擬;氣泡速度;縱橫比;振蕩幅度

中圖分類號：TV672.5 文獻標志碼：A

doi：10.3969/j.issn.1000-1379.2020.12.032

Abstract：The aim of this study was to detect the movement behavior and deformation characteristics of the air bubbles in siphon pipeline. It was essential to use numerical methods to figure out the movement of the air bubbles in the siphon pipeline. The results show that the bubble aspect ratio in the ascending tube is smaller and the bubble velocity is bigger. The change of flow has little effect on the aspect ratio and the bubble shape is consistent with its phase diagram in still water. During the progress of rising， the bubble oscillates laterally and the amplitude of oscillation first increases and then decreases with the increase of bubble diameter. The flow shows a wide difference in the downpipe， the movement direction of the bubble is also different and the bubble deformation is in agreement with the ascending tube. The bubbles in the middle pipe flow in the way with attaching the sidewall. Meanwhile， the bubble diameter increases with the decreasing of the bubble velocity. The bubble breaks at a certain diameter and smaller bubbles appear， then bubble velocity increases again. The results can provide a reference for the study of gas-liquid two-phase flow in siphon pipeline.

Key words： siphon; numerical methods; bubble velocity; aspect ratio; oscillation amplitude

水平管段距離長、真空度大的虹吸管輸水管道是坎兒井式地下水庫的重要組成部分。李琳等[1]、許史等[2-3]在新疆臺蘭河某地下水庫長距離虹吸管道水力學模型試驗中發(fā)現(xiàn)，安裝高度小于7 m時，虹吸管內(nèi)出現(xiàn)偽空化現(xiàn)象，流動介質(zhì)由單一液相流轉(zhuǎn)變?yōu)闅馀萘鳌⑦^渡流和氣團流。王夢婷[4]、張小瑩等[5]的試驗結(jié)果表明，虹吸管安裝高度為6～8 m時，各水頭運行時上行管均以氣泡流為主，進入水平管路后，低水頭運行時氣泡在運動過程中不斷聚合形成氣團和長度為30～60 cm的大氣囊，氣囊體積隨運行時間增加而增大，當其貫穿整個斷面時管路斷流。為了保證虹吸管正常運行，應盡量避免氣泡在一定條件下聚合成氣囊，而這一問題與氣泡在動水和負壓運行的管道內(nèi)的動力學行為密切相關(guān)。

氣泡以分散相的形式在水中的運動過程是一種復雜、不穩(wěn)定和非線性的水動力學現(xiàn)象。Sanada等[6]、Wang等[7]通過試驗研究和數(shù)值模擬對不同黏度液體中的氣泡運動進行了研究，得到了不同直徑下氣泡的形狀。程軍明等[8]對氣泡在靜水中上升和破裂過程進行了數(shù)值模擬，獲得了氣泡上升速度和氣泡形狀之間的關(guān)系。陳如艷等[9]通過試驗研究得到水平管內(nèi)氣液兩相流中液相流速、氣相流速和氣泡尺寸的相互影響規(guī)律。Bhaga等[10]研究了浮力氣泡在黏性液體中的上升運動過程，得到了氣泡形狀及終端上升速度與雷諾數(shù)Re、莫頓數(shù)Mo和埃奧特沃斯數(shù)Eo之間關(guān)系的氣泡相圖。Ellingsen等[11]使用高速攝像機對氣泡在各種局部流動條件及不同水質(zhì)情況下的動力學特性進行了研究，并得到了在不同條件下氣泡的運動及形變過程。Ohta等[12-13]先后運用VOF數(shù)值模型，揭示了當流動條件為低Mo和Eo時，氣泡上升行為與其初始形狀密切相關(guān)。前人的研究主要集中于靜水中氣泡的運動特性，而對于動水和負壓運行下虹吸管內(nèi)的氣泡運動特性尚未見報道。本文應用標準k-ε模型，同時結(jié)合VOF方法對虹吸管上、中、下行管路中單個氣泡在不同流量條件下的運動行為進行模擬，為探明虹吸管內(nèi)氣泡聚合形成氣囊和氣泡獨立運動產(chǎn)生條件奠定基礎，同時為實際工程中合理布置虹吸管路排氣裝置、保證管路正常運行提供設計依據(jù)。

1 湍流數(shù)學模型

1.1 控制方程

1.2 標準k-ε模型和VOF法

在整個虹吸管內(nèi)，除了彎頭處以外，其他地方流線基本相互平行，因此選擇計算量小、計算精度和收斂性均較好的標準k-ε模型進行計算。湍動能k和耗散率ε方程為式中：σk、σε分別為湍動能k和湍動能耗散率ε對應的Prandtl數(shù)，均為1.39;t為時間;Gk為湍動能k的生成項，Gk=μt（uixj+ujxj）uixj，其中μt=ρCμk2ε，Cμ=0.084 5;C1ε=1.44;C2ε=1.92;C3ε=0.09;Gb是因浮力影響引起的湍動能而產(chǎn)生的，Gb=βgiμtσtuixi，β=0.012，σt為紊動普朗特數(shù)，其他參數(shù)取值見文獻[15];YM為有壓縮流脈動膨脹對總耗散率的影響值，在本文中模擬的液相認為是不可壓縮的，因此不考慮該項。

精確描述氣液兩相的運動界面是研究氣泡行為的關(guān)鍵。VOF方法追蹤的是網(wǎng)格中流體體積，具有容易實現(xiàn)、計算量小和精度高等優(yōu)點，因此本文采用VOF法追蹤氣液交界面。當αw=0時，表示管道內(nèi)沒有水，被氣充滿;當水的體積分數(shù)αw=1時，表示管內(nèi)被水充滿;當0≤αw≤1時，表示管內(nèi)被水和氣充滿。αw的控制方程為

2 計算區(qū)域的離散及邊界條件

由于本文模擬的管道為圓形截面，氣泡為球形，管道與氣泡都具有軸對稱性，因此建立二維模型進行數(shù)值模擬。模型由上行管、中行管和下行管三部分構(gòu)成，其中上行管和下行管高度均為6 m，中行管長度為18 m，管徑為2 cm。分別對單個氣泡在上行管、中行管、下行管內(nèi)運動過程進行數(shù)值模擬研究。虹吸管模型邊界采用速度進口邊界條件，液相初始速度分別為0.06、0.15、0.2 m/s且垂直于管路的橫斷面，氣相的初始速度為0。管道出口為壓力邊界，相對壓強為0，管道壁面設為無滑移固壁，采用標準壁面函數(shù)進行修正。當t=0時，氣泡位于上行管距離管口50 cm中心處。計算區(qū)域網(wǎng)格劃分采用四邊形網(wǎng)格和三角形網(wǎng)格，由于整體劃分網(wǎng)格數(shù)目較多，計算時間過長，因此對上行管、中行管和下行管網(wǎng)格局部加密。氣泡在不同管段運動時，需運動一段距離后速度才能達到穩(wěn)定。在上行管中，對氣泡起點以上3 m范圍內(nèi)的管道進行加密;當氣泡在中行管運動時，對距離氣泡起點0.5 m處之后的3 m范圍內(nèi)的管道進行加密;下行管中，由于流量不同，氣泡運動方向也不同，因此針對不同流量分別做了加密方案，流量較小時，對氣泡起點以上3 m范圍內(nèi)的管道進行加密，流量較大時，對氣泡起點以下3 m范圍內(nèi)的管道進行加密;網(wǎng)格單元最小尺寸為4×10-4m。控制方程的離散采用有限體積法，從穩(wěn)定性、精確性、適用性方面考慮，控制方程中的對流項和擴散項均采用二階迎風格式進行離散。離散后的線性代數(shù)方程組采用SIMPLEC算法迭代求解，計算時間步長為1×10-4 s。

使用標準的k-ε模型和VOF法，通過FLUENT軟件進行模擬，設置表面張力系數(shù)σ=0.072 8 N/m，模擬虹吸管輸水流量分別為0.018 8、0.047 1 L/s，氣泡直徑分別為1、4、5、6、7 mm的管路內(nèi)氣泡運動行為特征。

3 計算結(jié)果與分析

3.1 數(shù)學模型的驗證

采用文獻[16]和文獻[17]的試驗結(jié)果對本文的數(shù)學模型及經(jīng)驗參數(shù)進行驗證。

3.1.1 物理模型試驗裝置

文獻[16]試驗裝置見圖1。試驗裝置高度為1 m，長、寬分別為20.32、7.56 cm，液面高度為0.8 m。使用空氣注射器注入所需大小的氣泡，通過水箱與管道的壓力差使得氣泡進入試驗裝置。氣泡出口位置為半球形倒置的杯罩，手動旋轉(zhuǎn)杯罩釋放氣泡，使用高速攝像機在測量區(qū)域進行測量。

3.1.2 數(shù)學模型的建立與對比

文獻[18]中指出計算區(qū)域的寬度大于4倍氣泡直徑時，可忽略邊界對氣泡運動特性的影響。為了驗證擬采用的數(shù)學模型，先建立高度為1 m、直徑為50 mm的一端封閉一端開口的豎直管道模型，出口邊界設為壓力出口（相對壓強為0），網(wǎng)格尺寸分別為6×10-4、5×10-4、4×10-4、3×10-4 m，分別模擬了直徑為4、6、8 mm的氣泡在靜水中的上升和變形過程，氣泡直徑為6 mm時不同網(wǎng)格尺寸下氣泡速度隨時間變化情況見圖2。由圖2可知，網(wǎng)格尺寸越小，氣泡速度越快，當網(wǎng)格尺寸小于4×10-4 m時，氣泡速度不再隨著網(wǎng)格尺寸的變化而發(fā)生變化。網(wǎng)格尺寸越小，需要的網(wǎng)格數(shù)量越大，計算耗時就越長，為了減少計算時間，選擇網(wǎng)格尺寸為4×10-4m進行計算。由圖2還可得到，網(wǎng)格尺寸為4×10-4m時模擬的氣泡速度與實測的結(jié)果吻合較好。當t=0.15 s時，氣泡速度模擬結(jié)果與試驗結(jié)果相差最小，相對誤差為1.7%;當t=0.25 s時，氣泡速度模擬結(jié)果與試驗結(jié)果相差最大，相對誤差為4.9%。

將氣泡直徑為6 mm時的計算結(jié)果與文獻[17]中實測的氣泡相圖及文獻[16]的試驗結(jié)果進行對比，結(jié)果如表1和圖2所示。表1為不同直徑氣泡在上升過程中不同時刻所對應的Re和Eo，并將氣泡形狀模擬結(jié)果與文獻[17]結(jié)果進行對比。從表1中可以看出3種直徑的氣泡運動過程與文獻[17]中氣泡相圖的結(jié)果一致。文獻[19]中提到，氣泡在流體中上升時，氣泡形狀主要取決于氣泡雷諾數(shù)Re、埃奧特沃斯數(shù)Eo和莫頓數(shù)Mo三個無因次準數(shù)，它們的表達式分別為

3.2 氣泡豎直運動和水平運動過程及形狀變化

在負壓條件下，不同流量時虹吸管上、中、下行管內(nèi)直徑d分別為5、6、7 mm的氣泡運動過程及形狀變化見圖3～圖5。氣泡形變過程用不同時刻的縱橫比（即氣泡短軸與長軸的比值）變化來表示，t=0時，縱橫比為1，氣泡縱橫比越接近1表明氣泡變形量越小。

圖3為上行管內(nèi)氣泡運動變形過程，t=0時氣泡為球形。當流量為0.018 8 L/s時，隨著氣泡上升，氣泡變成橢球狀，隨著直徑的增大，變形越來越明顯，如d=5 mm的氣泡經(jīng)過4 s后，縱橫比由t=0時的1變?yōu)?.64，而d=7 mm的氣泡經(jīng)過相同時間縱橫比由t=0時的1變?yōu)?.55。當流量由0.018 8 L/s增大至0.047 1 L/s時，氣泡仍然保持橢球狀上升，縱橫比隨氣泡直徑增大而減小，但是流量的增大對氣泡形狀的變化影響很小，如d=5 mm和d=7 mm的氣泡經(jīng)過6 s后縱橫比分別由1變?yōu)?.63和0.56，與流量為0.018 8 L/s時的縱橫比基本相同。

圖4為球形氣泡在下行管隨水流運動時氣泡變形過程。氣泡在下行或上行過程中主要經(jīng)歷了球形和橢球形的變形過程。計算不同時刻的Re、Eo和Mo，與文獻[17]的氣泡相圖進行對比，結(jié)果表明：動水條件下虹吸管上、下行管內(nèi)氣泡在不同時刻的形狀變化過程模擬結(jié)果與文獻[17]的靜水條件下氣泡運動變形過程一致。當流量為0.018 8 L/s時，因流速較小，流體曳力（物體在流體中有相對運動時，會受到流體的阻力，阻力大小由相對速度差決定）較小，在浮力作用下d=5 mm的氣泡在下行管內(nèi)向上運動，氣泡運動速度小于上行管中氣泡運動速度，此時下行管氣泡縱橫比較大，氣泡運動過程中變形量略小于上行管內(nèi)的。如t=4 s時，d=5 mm的氣泡在下行管內(nèi)氣泡縱橫比比上行管內(nèi)增加2%。當流量為0.047 1 L/s時，流體曳力大于浮力，氣泡隨水流一起向下運動，隨著氣泡速度的增大，相對速度差減小，流體曳力減小，下行管氣泡運動速度小于上行管內(nèi)的，因此氣泡變形量仍小于上行管內(nèi)的氣泡變形量，經(jīng)過4 s后，d=7 mm的氣泡縱橫比比上行管中增大了9%。

圖5為氣泡在中行管運動變形過程，不同流量下氣泡變形過程基本一致。當t=0時，氣泡為球形。當t>0時，在浮力和壁面黏附力作用下氣泡變?yōu)榘霗E球狀，自管道中心向上管壁運動，保持貼壁流運動。當Q=0.018 8 L/s時，經(jīng)過相同時間，氣泡縱橫比隨氣泡直徑的增大而減小。如d=5 mm和7 mm的氣泡經(jīng)過4 s后縱橫比分別由1變?yōu)?.24和0.17，對同一直徑氣泡而言，經(jīng)過相同時間，氣泡縱橫比隨流量增大而減小。如Q=0.018 8 L/s時，d=5 mm的氣泡經(jīng)過4 s后縱橫比較Q=0.047 1 L/s時增大了11%。

3.3 虹吸管內(nèi)不同直徑氣泡速度與縱橫比

圖6為Q=0.018 8 L/s時不同直徑氣泡的速度、縱橫比隨時間的變化情況。其中：v、v1、v2分別為氣泡直徑為5、6、7 mm時的氣泡運動速度，e、e1、e2分別為氣泡直徑為5、6、7 mm時的氣泡縱橫比。

氣泡直徑一定時，縱橫比越大，氣泡速度越小，這一規(guī)律與文獻[16]中提到的靜水條件下氣泡運動過程中縱橫比與其速度成反比的規(guī)律一致。流量增大時，氣泡縱橫比與速度的關(guān)系和小流量時一致（限于篇幅，文中未給出其他流量時的圖形）。從圖6（a）可以看出，氣泡直徑越大，速度越小，但不同直徑氣泡的速度相差不大，如d=5 mm和d=7 mm時，最大差為t=1.0 s時的12%，最小差為t=2 s時的2%。其原因是氣泡直徑增大，所受浮力增大，在浮力作用下，氣泡與液體之間的相對速度差增大，導致流體曳力增大，抑制氣泡上升，隨著氣泡直徑的增大，管壁黏附力將會阻礙氣泡運動，使得大、小氣泡的運動速度差別不明顯。

從圖6（b）可以看出，不同直徑的氣泡在下行管內(nèi)的流速均小于上行管內(nèi)的，且同一時刻氣泡縱橫比越大，速度越小。Q=0.018 8 L/s時，在下行管中，氣泡向上運動，由于水流方向與上行管內(nèi)水流方向相反，因此速度比上行管內(nèi)小。氣泡速度隨著氣泡直徑的增大略有減小，如氣泡直徑由5 mm增加到7 mm時，氣泡速度由0.16 m/s減小19%，其原因是較大的氣泡所受管壁黏附作用影響較大，使得阻力增大。當流量增加到0.047 1 L/s時，流體曳力增大，氣泡向下運動，氣泡的縱橫比與速度的關(guān)系也與上述一致;氣泡直徑對氣泡速度的影響很小。

從圖6（c）可以看出，在中行管內(nèi)不同直徑的氣泡先進行短暫的減速運動后達到終速度，與上行管相比速度有所減小，其原因是在中行管內(nèi)氣泡先向上運動，直到與管壁貼合，壁面黏附力增大，浮力不再做功，氣泡速度減小。當直徑d=5～6 mm時，氣泡速度會隨著氣泡直徑的增大而減小。當氣泡直徑進一步增大，如d=7 mm時，氣泡速度增大，其原因是氣泡直徑過大，運動過程中發(fā)生了破裂，形成了較小的氣泡，使得速度又出現(xiàn)增大的現(xiàn)象。由于中行管內(nèi)是貼壁流，因此氣泡縱橫比與上行管相比相差很大，如Q=0.018 8 L/s時，d=5 mm的氣泡經(jīng)過4 s后，中行管內(nèi)氣泡縱橫比較上行管內(nèi)減小65%。

3.4 氣泡運動軌跡

氣泡在上升過程中會發(fā)生橫向位移，其原因是氣泡在上升的過程中會發(fā)生形變使得受力不平衡，導致了橫向振蕩。將氣泡中心點的位置坐標（0，0.6）作為t=0時刻氣泡的運動位置，流量為0.047 1 L/s時不同直徑的氣泡在上升過程的運動軌跡見圖7，圖中x/d表示氣泡的中心點位置沿管徑方向的相對位置，y/hs表示氣泡的中心點位置與虹吸管高度的相對距離。由圖7可知，氣泡在上行管中運動時，隨著氣泡直徑的增大，振蕩幅度發(fā)生變化。當上升氣泡直徑較小時，如d=1 mm時，氣泡上升中的振蕩幅度較小，基本沿直線上升。隨著氣泡直徑的增大，如d=4 mm時，氣泡呈“之”字形上升。當氣泡直徑增大到7 mm時，氣泡的振蕩幅度卻開始減小，其原因是直徑較大的氣泡在上升過程中受到壁面黏滯作用的影響增大，對氣泡的振蕩起阻礙作用，使得振蕩幅度變小。在不同流量的情況下，氣泡的運動軌跡呈現(xiàn)相似的變化。

4 結(jié) 論

本研究模擬了單個氣泡在靜水中的上升過程，并與文獻[17]的氣泡相圖結(jié)果進行比較，驗證了模擬結(jié)果的正確性，在此基礎上對虹吸管道水流進行數(shù)值模擬研究，并對模擬結(jié)果進行分析。

（1）氣泡在上行管上升過程中，由球狀變?yōu)闄E球狀，縱橫比逐漸減小，氣泡直徑越大，縱橫比變化越劇烈，而流量的變化對氣泡縱橫比的影響很小。氣泡速度與縱橫比有關(guān)，縱橫比越大，氣泡速度越小。氣泡直徑較小時，氣泡沿直線上升，隨著氣泡直徑增大，氣泡發(fā)生振蕩，呈“之”字形上升。當氣泡直徑增大到一定程度時，氣泡受到壁面的黏滯作用影響增大，對氣泡的振蕩起阻礙作用，使得振蕩幅度變小。

（2）氣泡在下行管運動過程中，氣泡變形情況與上行管相似，但氣泡受到浮力和曳力作用，下行管內(nèi)氣泡速度和變形程度小于上行管內(nèi)的。

（3）氣泡在中行管運動過程中，氣泡經(jīng)過短暫時間運動至管道頂部，在此之后為貼壁運動，氣泡由球狀變?yōu)榘霗E球狀。隨著氣泡直徑的增大，氣泡縱橫比逐漸減小;縱橫比大小與流量有關(guān)，流量越大，縱橫比越小。氣泡速度隨氣泡直徑增大而逐漸減小，氣泡發(fā)生破裂后的氣泡速度大于原來氣泡的速度。

參考文獻：

[1] 李琳，邱秀云，許史，等.長距離虹吸管道輸水水力學模型試驗研究[J].南水北調(diào)與水利科技，2010，8（3）：106-109.

[2] 許史，李琳，邱秀云，等.長距離虹吸管輸水試驗研究初探[J].中國農(nóng)村水利水電，2010（3）：70-72.

[3] 許史.長距離虹吸管輸水試驗研究[D].烏魯木齊：新疆農(nóng)業(yè)大學，2010：13-24.

[4] 王夢婷.真空管道氣液兩相流動水力特性試驗研究[D].烏魯木齊：新疆農(nóng)業(yè)大學，2015：19-69.

[5] 張小瑩，李琳，譚義海，等.虹吸管道坡度對氣液兩相流動特性影響的試驗研究[J].農(nóng)業(yè)工程學報，2017，33（14）：122-129.

[6] SANADA T，SUGIHARA K，SHIROTA M，et al. Motion and Drag of a Single Bubble in Super-Purified Water [J]. Fluid Dynamics Research，2008，40（7/8）：534-545.

[7] WANG Huanran，LI Yang，YANG Dong，et al. On the Shape Feature of a Single Bubble Rising in Viscous Liquids [J]. Engineering Thermophysics，2009，30（9）：1492-1494.

[8] 程軍明，吳偉烽，聶娟，等.氣泡在靜水中上升破裂產(chǎn)生射流特性的數(shù)值模擬[J].江蘇大學學報（自然科學版），2014，35（5）：513-517.

[9] 陳如艷，仇性啟，王曉琦.水平管內(nèi)氣液兩相流研究[J].石油化工設備，2007，36（6）：37-38.

[10] BHAGA D，WEBER M E. Bubbles in Viscous Liquids： Shapes，Wake and Velocities [J]. Fluid Mech， 1981，105：61-85.

[11] ELLINGSEN K，RISSO F. On the Rise of an Ellipsoidal Bubble in Water：Oscillatory Paths and Liquid Induced Velocity [J]. Computer Phys， 2008，227：3358-3382.

[12] OHTA M，IMURA T， YOSHIDA Y， et al. A Computational Study of the Effect of Initial Bubble Conditions on the Motion of a Gas Bubble Rising in Viscous Liquids [J]. Multiphase Flow，2005，31（2）：223-237.

[13] BONOMETTI T，MAGNAUDET J. Transition from Spherical Cap to Toroidal Bubbles [J]. Phys Fluids，2006，18（5）：1-12.

[14] BRACKBILL J，KOTHE D. A Continuum Method for Modeling Surface Tension [J]. Journal of Computational Phvsics，1992，100：335-345.

[15] ISSA R I，OLIVEIRA. Numerical Prediction of Phase Separation in TwoPhase Flow Through TJunctions [J]. Computers Fluids，1994，23（2）：347-372.

[16] 劉柳.垂直上升管中氣泡動力學特性實驗研究[D].長沙：中南大學，2013：13-17.

[17] 李仲春，宋小明，姜勝耀，等.靜水中較大氣泡運動特性實驗研究[J].核動力工程，2015，36（1）：161-164.

[18] 何丹，李彥鵬，劉艷艷.初始形狀對浮升氣泡動力特性的影響[J].西安交通大學學報，2011，45（1）：43-47.

[19] 劉靜如.非牛頓流體中多氣泡相互作用、聚并與破裂過程的數(shù)值模擬[D].天津：天津大學，2014：1-10.

【責任編輯 許立新】