紀 健，傅曉寧，王 毅，李俊飛，紀寶強

(1.東華理工大學土木與建筑工程學院，南昌 330000；2.中國石油新疆油田分公司，克拉瑪依 834000)

目前關于氣液兩相流聲場的研究和模擬幾乎空白，氣液兩相流管道流動的復雜性導致模擬存在很大誤差。國內外關于氣液兩相流場與聲場研究主要基于實驗分析，但實驗條件下能測量的參數(shù)很有限。通過數(shù)值模擬分析泄漏過程各個參數(shù)變化，補充實驗在描述各流場參數(shù)方面的不足。水平管道的氣液兩相流的流場噪聲產生原因是兩相流中的氣液持液率變化、湍流產生的壓力脈動和速度脈動，因此分析水平管道氣液兩相流壓力和速度變化有助于對氣液兩相流聲場特性的掌握。

通過泄漏開啟過程參數(shù)變化分析三種流型工況下泄漏過程的聲波特性，即數(shù)值模擬的目的就是分析不同流型下氣液兩相流管道泄漏工況下的流場參數(shù)，主要包括管道發(fā)生泄漏前后，泄漏孔處壓力、流速、湍流強度、渦量等流體參數(shù)變化，進而分析各流型工況下的泄漏發(fā)生時的流體發(fā)聲機理。

1 模型選擇、計算方法和邊界條件

多相流流體控制方程和數(shù)學模型：計算流體力學軟件通常采用納維-斯托克斯方程作為控制方程對氣液兩相流管道進行數(shù)值模擬[1-5]，其連續(xù)性方程、動量守恒方程和能量守恒方程為

(1)

(2)

(3)

式中：ρ為氣體密度，kg/m3；u為氣體速度，m/s；p為氣體靜壓，kPa；τij為黏性應力張量；E為總能，J；qj為熱通量，J/s；t為時間，s。

數(shù)值模擬軟件常用湍流求解的方法很多，相關研究成熟。文獻調研對比分析發(fā)現(xiàn)[6-11]，Realizablek-ε模型可以在強流線彎曲、旋渦和旋轉方面表現(xiàn)很好，且在預測流動分離和復雜二次流、圓形射流等方面模擬效果極好，因此選擇其作為湍流模型[12]，進行湍動能、耗散率輸運和湍流黏度計算。

計算流體力學提供的三種最基本多相流模型是Mixture模型、Eulerian模型和VOF(volume of fulent model)模型[13-16]，各個模型技術成熟，均得到廣泛應用和深入研究。模擬的主要目的是研究氣液兩相管內泄漏前后流動形態(tài)和氣液兩相界面不斷變化。VOF模型是基于求解體積分數(shù)方程來實現(xiàn)準確追蹤氣液兩相流的界面變化情況的，主要適用于有明顯氣液相間界面的流動模擬。因此選擇VOF模型瞬態(tài)求解，氣液兩相界面相互作用方法采用幾何重構格式，其表面張力采用連續(xù)表面力模型[17-18]。

水平氣液兩相流動過程中絕熱，氣體可壓縮，液體不可壓縮。設空氣為主相，水為次相。模擬和實驗條件均為低壓系統(tǒng)，氣體符合理想氣體狀態(tài)方程。水平管道中的流動狀態(tài)具有強瞬態(tài)特性，計算過程用壓力的隱式算子分割法(pressure implicit with splitting of operators，PISO)算法，時間步長選取1×10-6s。

計算模擬過程，監(jiān)控球閥模擬泄漏過程中的泄漏孔處的升力系數(shù)隨時間變化，以及泄漏孔處的動態(tài)壓力變化。穩(wěn)定條件需要滿足保障檢測點參數(shù)波動變化曲線穩(wěn)定以及基本的氣液兩相流動模擬收斂條件。設計入口邊界條件利用實際流動參數(shù)計算獲得。氣液入口設為與實驗過程相對應的質量流量入口，出口設置為壓力出口。

2 幾何模型、網格劃分、泄漏方式

依據(jù)實驗建立泄漏仿真模型依據(jù)前期實驗管道完成，由主管道和泄漏模擬分支管道組成，原理是利用一個T形支管安裝球閥，通過開啟球閥并更換泄漏孔徑的方法模擬各種工況泄漏，主管道參數(shù)與實驗相同，管徑37 mm、管長1 200 mm；泄漏模擬分支管道管徑10 mm、管長120 mm。主管道和分支管道中心線交點定義為原點。泄漏管道球閥的球芯中心位置為直角坐標(0，80)處，孔板壁厚2 mm，泄漏孔徑選取6 mm，氣體從主干管道左端流入，從主干管道右端流出，具體尺寸方位如圖1所示。

圖1 簡化的管道泄漏幾何模型Fig.1 Simplified geometric model of pipeline leakage

網格生成直接關系到計算問題的成敗，對數(shù)值模擬計算至關重要。網格主要分為結構網格和非結構網格兩大類。結構化網格是只包含四邊形或者六面體，非結構化網格是三角形和四面體。對于二維模擬，選用四面體網。針對水平管道流動特性和幾何形狀的復雜程度，同時考慮計算精度，時間耗費等問題，采用以下策略進行網格劃分：整體采用結構網格，氣液入口處兩相相互混合，速度梯度較大，泄漏支管管道流線發(fā)生大角度彎曲，且管道結構本身具有細長特點，為實現(xiàn)泄漏流場和聲場的轉換的準確，對氣液入口和泄漏支管上、主分支管道的相交處網格進行加密。依據(jù)以上策略得到兩種幾何模型的兩種網格數(shù)量分別為網格1和網格2，分別為25 368和64 000。對兩種不同密度的網格的計算結果進行對比，選取泄漏開啟過程壓力參數(shù)作為敏感性分析的特征參數(shù)，結果圖2所示，兩種網格質量的計算結果非常相近。因此充分考慮計算精度要求和效率，選定網格1進行數(shù)值模擬。

圖2 兩種不同網格相應的壓力變化曲線Fig.2 Pressure change cuver of two different meshing

利用動網格計算需要在用Gambit中進行網格劃分時，對閥芯與泄漏管道邊界設置如圖3所示。

圖3 幾何模型的網格劃分和邊界設置Fig.3 Structured mesh of geometries and boundary setup

圖4 泄漏模擬球閥閥門轉動規(guī)律Fig.4 Rule of valve rotation of leakage simulated ball valve

在實際工程中管道發(fā)生泄漏是在瞬間完成的。采用Fluent中動網格技術模擬泄漏并定義閥門開啟時間設置[19]。假設閥門位移與時間關系如圖4所示，時間t時刻閥門開啟，閥門開啟時間為0.1 s，保持開啟狀態(tài)直至流動穩(wěn)定。

3 數(shù)值模擬結果分析

重點分析泄漏導致水平管道內氣液兩相流場內的總壓力、流速、湍流強度和渦量這四個流動參數(shù)量的變化情況。為了更好地對比分析不同流型下泄漏對流動的影響，選擇實驗過程中三種流型對應氣液折算速度下的氣液質量流量數(shù)據(jù)，如表1所示。

表1 各流型模擬數(shù)據(jù)選取Table 1 Selection of flow pattern simulation data

以實驗模型為研究對象，進行氣液兩相流水平管道上部泄漏流場模擬。以表1中分層流對應氣液質量流量設置氣液入口邊界，泄漏孔徑6 mm為例，進行球閥轉動過程流場的數(shù)值模擬。氣液混合速度矢量變化如圖5所示，該云圖代表氣液混合速度的大小，各圖依次為球閥轉動間隔0.01 s的氣液速度變化，即從球閥開始轉動到完成轉動的各個0.01 s時間的球閥開啟不同程度的速度變化。由圖5可以得到以下結論。

圖5 不同時刻流體速度矢量圖Fig.5 Vector diagram of fluid velocity at different times

(1)當管道內流動狀態(tài)為分層流時，球閥的轉動從開啟到完成整個過程，上部泄漏支管道內流體為氣體，幾乎沒有液體存在。泄漏氣體瞬時充滿整個分支管道。氣體噴射出管道通常發(fā)生在球閥旋轉的第0.01秒稍后，球閥在第0.01秒才能實現(xiàn)真正意義的開啟。

(2)泄漏閥門真正開啟時，泄漏流體在閥門出口處形成高速流動，空氣突然膨脹、減速、擴張或碰撞產生雷諾應力或剪切應力形成湍流導致空氣動力發(fā)生。存在某一時間段，支管道中的泄漏流體流動速度大于主管道中的流體流動速度。主管內的液體需要克服重力作用才能被氣體攜帶出泄漏管道，通常情況需要較高的氣體流速。

球閥開啟前，管道內處于穩(wěn)定流動，主管道和分支管內的氣液均處于相對穩(wěn)定的狀態(tài)。球閥開啟瞬間，由于球閥內的壓力遠小于支管道內球閥前的壓力，在內外在壓差的作用下，氣體迅速流入球閥，泄漏支管中流體的速度瞬間增大。隨著球閥繼續(xù)轉動，氣體流入支管道，并隨著球閥的不斷開啟，氣速不斷增大。此時實驗過程流體動力學特性和數(shù)值模擬分析結論相同。數(shù)值模擬泄漏時，在球閥的開啟過程進行泄漏孔口的升力系數(shù)Cl監(jiān)測。Cl是由監(jiān)測點所受升力除以參考值計算的動壓得到，由升力系數(shù)定義得知其為無量綱量。

穩(wěn)定流動時，分支管道泄漏孔口升力系數(shù)為零。模擬泄漏過程，球閥開啟動作定義發(fā)生在0.5～0.6 s內，即圖4中t=0.5 s。球閥球芯開始轉動，但閥門并沒有實現(xiàn)實際性的開啟，即圖3中交界面還不存在，此時升力系數(shù)仍然為零。當隨著球閥的繼續(xù)轉動，球閥實現(xiàn)真正開啟。得到各個流型球閥模擬泄漏過程中的泄漏孔口處升力系數(shù)Cl隨時間變化曲線如圖6所示。

圖6 泄漏孔升力系數(shù)隨時間變化Fig.6 The lift coefficient of leakage hole changes with time

分層流、波浪流、段塞流泄漏模擬在球閥開啟過程，即球芯發(fā)生轉動過程中，各個流型工況下升力系數(shù)均發(fā)生了極大變化，區(qū)別于泄漏發(fā)生前均趨于水平為零值，極容易分辨識別，與理論分析相一致，即泄漏孔口處由無任何流體通過時不受任何作用力；而閥芯轉動，當有流體流入分支管道，泄漏孔口開始受到力的作用，此時升力系數(shù)發(fā)生變化。對比分析圖6升力系數(shù)Cl隨著流型從分層流、波浪流、段塞流依次變化，其各自流型下的升力系數(shù)最大值不斷增大。與理論分析相同，即升力系數(shù)與流體的動壓有關，動壓越大其絕對值越大。通過球閥的閥芯轉動實現(xiàn)開啟在模擬過程只發(fā)生一次，與實際工況相符，為非周期性。

4 數(shù)值模擬與實驗驗證

結合模擬設定，升力系數(shù)即為泄漏孔處受到的升力除以參考值計算得到的泄漏孔口處的動壓。此時最小升力系數(shù)與實驗過程泄漏發(fā)生時刻產生的最小動態(tài)壓力幅值相對應。在各個流型下分別對應的最小升力系數(shù)Cl變化曲線和最小動態(tài)壓力信號p變化曲線之間存在一個關系為:Cl/p=常數(shù)。若兩者間的關系成立可以證明泄漏過程數(shù)值模擬符合實驗。通過圖6數(shù)值模擬結果數(shù)據(jù)分析，得到最小升力系數(shù)值。

由表1實驗數(shù)據(jù)選取對應各個流型工況下氣液流量模擬泄漏，泄漏孔徑為6 mm采集到相應各個流型下的動態(tài)壓力信號，得最小升力系數(shù)和最小動態(tài)壓力值如表2所示。分別對表2中最小升力系數(shù)和最小動態(tài)壓力幅值離散點進行相關性驗證，得到升力系數(shù)和動態(tài)壓力變化擬合曲線和相關系數(shù)，如圖7所示。

表2 球閥開啟模擬與實驗過程參數(shù)對比Table 2 Comparison of simulation and experimental process parameters of ball valve opening

圖7 數(shù)值模擬升力系數(shù)和實驗動壓擬合關系曲線Fig.7 Fitting curve of lift coefficient of numerical simulation and experimental dynamic pressure

擬合得到相關系數(shù)為R2=0.999 7，說明兩者之間具有很好的相關性，且由此可得Cl/p=0.000 1，即泄漏模擬開啟過程得到的數(shù)值模擬結果與泄漏實驗結果相符合。

5 泄漏開啟過程參數(shù)分析

5.1 分層流泄漏開啟模擬參數(shù)分析

分析分層流泄漏開啟流場，首先數(shù)值模擬4 mm和6 mm泄漏孔徑時泄漏開啟過程，流場分析得總壓力、流速、渦量、湍流強度變化如圖8所示。

圖8 分層流泄漏模擬過程每0.01 s各參數(shù)變化曲線Fig.8 Simulates the variation curve of each parameter every 0.01 s in the leakage process of stratified flow

其中，球閥轉動0.1 s內的每間隔0.01 s采樣保存一次，共采集11個點的參數(shù)變化。采樣點1為球閥球芯即將發(fā)生轉動，采樣點11處為完成球芯轉動，球閥開啟過程完成。由圖8分析發(fā)現(xiàn)各參數(shù)變化具有以下特點。

由圖8(a)發(fā)現(xiàn)，球閥真正開啟過程瞬間為圖中第2個采樣點，即第0.01 s處，此時總壓力曲線出現(xiàn)最大值，說明球閥真正轉動開始瞬間管道內的總壓最大，此時泄漏瞬時聲波信號幅值最大，即管道內的氣體從閥芯連通支管部分向外泄漏，泄漏孔徑越大，瞬時聲波信號越強烈，波動幅值越大。

由圖8(b)發(fā)現(xiàn)，在第0.01 s處為開啟過程流體流速變化曲線的最小值點。流速在開啟瞬間迅速減小到最低，然后隨著球閥的轉動時間增長再逐漸增加。出現(xiàn)這種現(xiàn)象的原因，同理第一個采樣點沒有泄漏發(fā)生，此時刻泄漏前最后穩(wěn)定流動時刻，此時的管道內流體流動最大值仍然出現(xiàn)在主管道，該值為主管道內最大氣液流速。泄漏發(fā)生在泄漏孔發(fā)生相分離，主管道和分支管道出現(xiàn)氣液重新分配，主管道減小而分支管的流速增加，壓力勢能還未完全轉換為流體動能。隨著球閥球芯轉動，越來越多氣體流入分支管道，球閥完全開啟后泄漏流體流速到達最大。隨泄漏孔徑增大，泄漏流體流速增大。

由圖8(c)發(fā)現(xiàn)，渦量在球閥開啟過程呈波浪式變化。球閥開始轉動時先增大后減小，然后又增大再減小直到開啟完成，變小趨于穩(wěn)定。整個過程渦量具有兩個極大值，說明在球閥的開啟過程中，球閥先造成流體擾動引起渦量發(fā)生很大變化，而流體又與管壁相互作用，造成渦量不穩(wěn)定。隨著球閥開啟的完成，分層流渦量大小恢復到泄漏前大小，此時泄漏孔徑對渦量影響不大。由理論分析知渦量主要產生偶極子聲源，球閥開始轉動前幾個采樣點造成管道內流場產生渦量很大。

由圖8(d)發(fā)現(xiàn)，隨著球閥開啟，管道內湍流強度不斷減少，該最大湍流強度值通過分析發(fā)現(xiàn)通常出現(xiàn)在管道出口處，在分層流上部泄漏時泄漏流體只有氣體噴射，說明在球閥開啟的瞬間湍流強度最大，此時湍流噪聲作為產生四極子聲源是分層流泄漏瞬時聲波產生的主要原因之一。

分層流管道閥門球芯開始轉動，連通分支管瞬間，管道內的壓力最大，內外壓差最大，此時氣液泄漏流體流速最小。球閥開始轉動，附近產生的流體總壓、渦量、湍流強度變化出現(xiàn)最大值，說明分層流泄漏聲波信號主要產生于球閥開啟瞬間。分層流泄漏聲源組成主要為偶極子聲源和四極子聲源。隨著泄漏孔徑增大，球閥開啟過程的最大壓力和流速值增大，渦量和湍流強度的變化劇烈程度增加。

5.2 波浪流泄漏開啟模擬參數(shù)分析

同理，分析波浪流泄漏流場，各參數(shù)最大值變化如圖9所示，分析發(fā)現(xiàn)各參數(shù)變化具有以下特點。

圖9 波浪流泄漏模擬過程每0.01 s各參數(shù)變化曲線Fig.9 Simulates the variation curve of each parameter every 0.01 s in the leakage process of wave flow

由圖9(a)發(fā)現(xiàn)，波浪流與分層流總壓變化規(guī)律相同，球閥轉動開始瞬間即第2個采樣點處，管道內的總壓最大，即此時管道內外壓差最大產生瞬時泄漏聲波幅值最大。前文理論分析得知波浪流的管道流動背景較大，這種泄漏聲波信號與分層流相比容易被淹沒在管道的流動背景噪聲中。

由圖9(b)發(fā)現(xiàn)，球閥在開啟瞬間，最大流體流速曲線在球閥轉動第0.01 s處流速大大減小，在第0.02 s到達最小值。球閥均隨著開啟時間延長，氣液流速逐漸增大但是增大趨勢變緩。產生這一現(xiàn)象的原因：管道內流體波浪流動特性，氣液相界面為高速流動的氣體攜帶液體向前波動。球閥開啟，氣體迅速擴散出管道，管道內氣體流速大大減小。隨著球閥繼續(xù)開啟，支管內氣量越來越多，流速越來越大，同時主管道內的總流速降低。

由圖9(c)發(fā)現(xiàn)，球閥真正開啟時刻、渦量幅值急劇減小。波浪流穩(wěn)定流動具有氣液間相互作用大、渦量大的特點。球閥瞬間開啟造成管道內氣體流速減低，減小了氣液兩相間的相互作用，從而管道內流動背景噪聲減小。整個球閥轉動過程造成分支管道內渦量增加，相對于其開啟大大減緩主管道內的渦量相比作用仍較小。

由圖9(d)發(fā)現(xiàn)，整個球閥開啟過程湍流強度變化呈現(xiàn)波浪式變化。因為波浪流管道內穩(wěn)定流動時，湍流變化同樣呈現(xiàn)波浪式變化，球閥開啟在一定程度上減緩了主管道內的流體湍流強度，分支管道內的湍流強度變化與主管道變化相比很小。

總體上說，波浪流發(fā)生泄漏，球閥開啟瞬間，管道內外壓差最大，產生泄漏瞬時聲波信號幅值最大。球閥在真正開啟時刻附近造成氣液流速大大降低且管道內渦量急劇減少，一定程度上減緩了管道內的湍流強度。說明上部泄漏大大削弱管道內的流動噪聲信號。因此波浪流發(fā)生上部泄漏時，可通過分析泄漏發(fā)生后管道內的流動噪聲的明顯減小來進行分析判斷。

5.3 段塞流泄漏開啟模擬參數(shù)分析

同理，分析段塞流泄漏流場，各參數(shù)最大值的變化如圖10所示。分析發(fā)現(xiàn)，主要有以下結論。

圖10 段塞流泄漏模擬過程每0.01 s各參數(shù)變化曲線Fig.10 Simulates the variation curve of each parameter every 0.01 s in the leakage process of slug flow

(1)球閥在開啟過程，管道內總壓波動性不斷加強，且泄漏孔徑越大波動性越強。流速在開啟過程中出現(xiàn)最大峰值，流速變化具有波動性。湍流強度在泄漏過程呈先減小后增加趨勢，但是球閥未開啟與完成球閥開啟后的差值不大。

(2)渦量在球閥球芯開啟過程中變化不明顯，在球閥完全開啟后渦量增加，遠大于未開啟動狀態(tài)，說明此時管道內產生的泄漏聲源最大。結合文中段塞流泄漏流動特征以及泄漏持續(xù)聲源發(fā)生的時間分析得知，渦量可能是段塞流泄漏持續(xù)聲源的重要組成部分。通過分析管道泄漏發(fā)生前后渦量變化可實現(xiàn)泄漏判斷。

5.4 不同流型泄漏開啟過程模擬對比

對分層流、波浪流和段塞流不同工況下，泄漏開啟過程中各參數(shù)即氣液流體的壓力、流速、湍流強度和渦量最大值變化進行分析對比，得到不同流型泄漏流場的泄漏開啟時各參數(shù)特征如下。

5.4.1 分層流

球閥開啟瞬間，管道內壓力最大、流速最小，內外壓差最大，產生泄漏聲波信號最大幅值。開啟附近的渦量和湍流強度最大，泄漏聲源主要由偶極子和四極子聲源組成。

5.4.2 波浪流

球閥開啟瞬間，管道內壓力最大、流速最小，內外壓差最大，產生泄漏瞬時聲波信號，渦量最小，湍動強度呈減小趨勢。渦量和湍流強度均減小可知流動噪聲減小，泄漏前后流動噪聲變化可以進行泄漏判斷。

5.4.3 段塞流

球閥開啟過程，管道內總壓波動性不斷加強，且出現(xiàn)流速最大值；球閥完全開啟時渦量最大。泄漏發(fā)生前后渦量變化可以進行上部泄漏判斷。

由表3知，分層流通過泄漏發(fā)生時的瞬時泄漏聲波信號就可以很好地進行泄漏判斷。而波浪流和段塞流的泄漏判斷必須依賴于對瞬時泄漏聲波信號和泄漏后管道內流體流動狀態(tài)，要得到較為準確的泄漏判斷，需要綜合分析這兩種信號的泄漏特性，這與前期獲取的實驗數(shù)據(jù)分析相一致。

6 結論

采用簡化的幾何模型，基于Fluent動網格技術通過UDF編程實現(xiàn)氣液兩相流管道泄漏開啟過程和完成的流場模擬，從泄漏流場各參數(shù)變化分析泄漏聲源特性與理論相驗證，是為了補充前期的實驗在采集流動參數(shù)方面的不足，主要得到如下結論。

(1)分層流在球閥真正開啟的瞬間，管道內總壓、渦量、湍流強度均為最大值，流速為最小值。泄漏聲波信號主要產生于球閥開啟瞬間，瞬時聲源能夠很好地進行泄漏判斷。

(2)波浪流球閥開啟瞬間，管道內總壓、渦量均為最大值，流速、湍流強度均為最小值。上部泄漏后的管道內壓力和流速較穩(wěn)定。

(3)段塞流在球閥開啟過程，總壓、氣液流速和湍流強度變化規(guī)律性不強；開啟過程會出現(xiàn)流速最大值；球閥完全開啟時渦量最大。需要主要靠泄漏完成流場參數(shù)變化來進行判斷。通過分析泄漏后流場參數(shù)變化發(fā)現(xiàn)，上部泄漏主要造成泄漏孔后主干道內流體流動趨于流動狀態(tài)不斷變化。

需要注意的是，為了彌補聲波泄漏檢測實驗過程中無法測量的參數(shù)問題而進行該理論研究，不可代替實驗。