豎直彎管中含固體顆粒分散泡狀流沖蝕預(yù)測(cè)

彭文山，馬力，侯健，曹學(xué)文，韓明一

（1.中國(guó)船舶重工集團(tuán)公司第七二五研究所 海洋腐蝕與防護(hù)重點(diǎn)實(shí)驗(yàn)室，山東 青島 266237；2.中國(guó)石油大學(xué)（華東）儲(chǔ)運(yùn)與建筑工程學(xué)院，山東 青島 266580；3.廊坊中油朗威工程項(xiàng)目管理有限公司，河北 廊坊 065000）

隨著我國(guó)油氣開(kāi)采由陸地向海洋轉(zhuǎn)移，海底管道投用量持續(xù)增加[1]，海洋混輸管網(wǎng)得到不斷應(yīng)用?；燧敼艿乐辛黧w流動(dòng)復(fù)雜，流體對(duì)管壁沖刷嚴(yán)重，油井出砂量及原油含水量隨著油氣田開(kāi)采時(shí)間的增長(zhǎng)而增加。油氣開(kāi)采及運(yùn)輸過(guò)程中，管內(nèi)固體顆粒隨油氣流動(dòng)，極大地加劇了管內(nèi)介質(zhì)對(duì)管壁的沖擊，極易造成管道沖蝕破壞，導(dǎo)致危險(xiǎn)事故[2-3]。因此，提出準(zhǔn)確的多相流沖蝕計(jì)算方法對(duì)于預(yù)測(cè)管道可能出現(xiàn)的沖蝕最嚴(yán)重位置，進(jìn)而保障管輸安全、減少經(jīng)濟(jì)損失意義重大。

泡狀流在垂直上升氣液兩相流中是常見(jiàn)的流型[4]，例如海底垂直上升管道以及海洋立管中。含固體顆粒泡狀流中，顆粒受氣液流動(dòng)影響，空間分布隨機(jī)性大，而管道沖蝕最嚴(yán)重位置與顆粒的分布直接相關(guān)，要準(zhǔn)確仿真多相流沖蝕，需考慮氣體-液體-顆粒-管道壁面之間的多向耦合，研究具有挑戰(zhàn)性。國(guó)外部分學(xué)者開(kāi)展了含固體顆粒多相流沖蝕實(shí)驗(yàn)，在實(shí)驗(yàn)結(jié)果基礎(chǔ)上提出了一些經(jīng)驗(yàn)和半經(jīng)驗(yàn)沖蝕計(jì)算公式[5-7]，但這些公式僅可計(jì)算最大沖蝕速率，局限性較大。Chen 等人[8]提出采用將氣液兩相簡(jiǎn)化成混合流體相來(lái)計(jì)算管道內(nèi)顆粒的沖蝕，計(jì)算獲得的沖蝕結(jié)果偏小，且無(wú)法獲得管道內(nèi)部氣液分布及顆粒運(yùn)動(dòng)情況。因此以上方法僅適應(yīng)于估算管道沖蝕速率，并不能正確解釋管道沖蝕機(jī)理。因此，提出了一種基于流體體積模型（VOF）模型和離散相模型（DPM）模型的多相流瞬態(tài)沖蝕仿真計(jì)算方法，研究管內(nèi)氣液分布、顆粒運(yùn)動(dòng)對(duì)沖蝕的影響，并與經(jīng)驗(yàn)?zāi)Ｐ陀?jì)算結(jié)果、簡(jiǎn)化計(jì)算流體動(dòng)力學(xué)（CFD）計(jì)算結(jié)果及實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，驗(yàn)證該模型準(zhǔn)確性，為含固體顆粒管道的分散泡狀流沖蝕預(yù)測(cè)提供參考。

1 分散泡狀流沖蝕求解

管道內(nèi)為氣液固三相流，氣液兩相為連續(xù)相，固體顆粒為離散相，采用DPM 模型完成固體顆粒的計(jì)算。分散泡狀流沖蝕求解分為以下三步：流場(chǎng)計(jì)算、顆粒追蹤及沖蝕速率求解。采用Eulerian-Lagrangian方法[9]，求解氣液兩相連續(xù)相流場(chǎng)和離散相顆粒受力方程，進(jìn)行顆粒軌跡計(jì)算，沖蝕計(jì)算模型提取前兩步計(jì)算獲得的信息，計(jì)算得到?jīng)_蝕速率。

1.1 連續(xù)相流場(chǎng)計(jì)算

標(biāo)準(zhǔn)k-ε 模型在計(jì)算管道內(nèi)部不同截面處流體速度的誤差較小[10]，選用該模型計(jì)算管道內(nèi)部流場(chǎng)。使用VOF 模型實(shí)現(xiàn)管內(nèi)分散泡狀流的求解。另外，考慮表面張力的作用，采用Brackbill 等人[11]提出的連續(xù)表面力模型（CSF）。利用這個(gè)模型，VOF 模型中的附加表面張力通過(guò)源項(xiàng)的方式添加到動(dòng)量方程中。

1.2 離散相控制方程

顆粒軌跡追蹤過(guò)程中不考慮顆粒之間的相互作用，單位質(zhì)量固體顆粒在分散泡狀流氣液兩相流中的受力方程為：

1.3 沖蝕計(jì)算模型

目前已有許多沖蝕模型被應(yīng)用于沖蝕計(jì)算[13]，Liu等人[14]對(duì)沖蝕計(jì)算研究發(fā)現(xiàn)，廣島大學(xué)Oka 等人[15-16]提出的沖蝕模型能較為準(zhǔn)確地預(yù)測(cè)多相流條件下管道的沖蝕速率。另外，模型中考慮了較多的影響因素，適應(yīng)范圍較廣。使用該模型進(jìn)行沖蝕計(jì)算：

式中：ρw是管材密度，kg/m3；α 為顆粒沖擊角度，rad；ρw為管材密度，kg/m3；HV 是管材維氏硬度；uC為參考顆粒速度，m/s；dC為參考直徑，m；k0、k1、k2、k3、n1、n2、dC、uC為常數(shù)，見(jiàn)表1[15-16]。

表1 Oka 等沖蝕模型參數(shù)Tab.1 Parameters of erosion model proposed by Oka et al.

2 仿真模型建立

2.1 管道基本參數(shù)

Birchenough 等[17]進(jìn)行沖蝕實(shí)驗(yàn)，獲得一組豎直管道中空氣-水系統(tǒng)分散泡狀流沖蝕的實(shí)驗(yàn)結(jié)果，在進(jìn)行數(shù)值分析時(shí)，與該組實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比驗(yàn)證。Birchenough 等人的實(shí)驗(yàn)管道管徑為49 mm，工作壓力為0.2 MPa，管道材料為碳鋼，仿真參數(shù)值見(jiàn)表2，D 為管道直徑，R/D 為彎徑比，VSG、VSL分別為氣相和液相的表觀流速，ρl、ρg、ρp分別為液體、氣體、砂粒密度。彎管模型由進(jìn)口直管段、彎頭和出口直管段組成。為使得管內(nèi)流動(dòng)充分發(fā)展，并考慮模型數(shù)值計(jì)算消耗，入口直管段長(zhǎng)度取為50 倍管徑（50D），出口直管段長(zhǎng)度取為25 倍管徑（25D），如圖1 所示。水的表面張力為0.073 N/m。

表2 分散泡狀流仿真參數(shù)Tab.2 Simulation parameters of dispersed bubble flow

圖1 豎直彎管幾何模型Fig.1 Geometry model of the vertical pipe bend

2.2 網(wǎng)格劃分

為提高VOF 模型與DPM 模型耦合計(jì)算的穩(wěn)定性，計(jì)算模型采用規(guī)則四邊形面網(wǎng)格，采用六面體彎管體網(wǎng)格，網(wǎng)格劃分如圖2 所示，豎直彎管三維模型網(wǎng)格總數(shù)為524 000。

圖2 豎直彎管網(wǎng)格劃分Fig.2 Computational mesh of the vertical pipe bend

2.3 邊界條件

模型入口采用速度入口邊界，速度大小設(shè)置為氣液混合速度，出口采用壓力出口。入口水力直徑為49 mm，湍流強(qiáng)度經(jīng)計(jì)算為3.2%。DPM 模型中進(jìn)口和出口采用逃逸（Escape）條件，壁面采用反彈（Reflect）條件，采用面射流源噴射固體顆粒，固體顆粒初始速度與流體入口速度相同。管壁壁面粗糙度常數(shù)設(shè)置為0.5，并設(shè)定為“靜止壁面”和“無(wú)滑移壁面”。壁面設(shè)置恢復(fù)系數(shù)來(lái)計(jì)算固體顆粒碰撞壁面前后的速度變化，計(jì)算過(guò)程中使用Grant 和Tabakoff[18]提出的隨機(jī)顆粒-壁面碰撞反彈模型的恢復(fù)系數(shù)，形式如下：

式中：en為法向恢復(fù)系數(shù)；et為切向恢復(fù)系數(shù)。

2.4 數(shù)值求解

使用FLUENT 軟件完成沖蝕計(jì)算，基于VOF 模型和DPM 模型的CFD 沖蝕計(jì)算中，壓力-速度耦合采用PISO 方式，壓力離散采用PRESTO!方式，動(dòng)量方程、湍動(dòng)能方程、湍流擴(kuò)散率方程均使用QUICK格式，體積分?jǐn)?shù)離散采用Geo-Reconstruct 方式，來(lái)提高計(jì)算精度。為了使收斂比較平穩(wěn)，并保證模擬結(jié)果的可靠性，時(shí)間步長(zhǎng)設(shè)置為0.001 s。顆粒同樣采用瞬態(tài)計(jì)算，計(jì)算開(kāi)始前打開(kāi)離散相模型，加入離散相粒子，追蹤時(shí)間步長(zhǎng)與連續(xù)相計(jì)算步長(zhǎng)相同。在瞬態(tài)計(jì)算過(guò)程中，每計(jì)算一步，固體顆粒在入口處“噴射”一次，入口處有500 個(gè)網(wǎng)格單元，則入口處每一時(shí)間步要入射500 個(gè)固體顆粒?？偟挠?jì)算時(shí)間為8 s，時(shí)間步長(zhǎng)為0.001 s，則數(shù)值計(jì)算共噴射進(jìn)入管道的顆粒數(shù)目為4 000 000 個(gè)。顆粒軌跡追蹤時(shí)，假定顆粒不影響氣液兩相流動(dòng)，僅考慮氣液兩相對(duì)顆粒運(yùn)動(dòng)的影響，即單向耦合計(jì)算。為了便于分析管道內(nèi)部氣液固多相的運(yùn)動(dòng)情況，分別在上、下游直管段中間位置以及彎頭45°截面處設(shè)置監(jiān)控面，對(duì)瞬態(tài)計(jì)算結(jié)果進(jìn)行監(jiān)控保存如圖1 中紅色虛線位置。

3 結(jié)果分析

3.1 沖蝕結(jié)果對(duì)比

基于VOF 模型和DPM 模型的CFD 沖蝕計(jì)算中，沖蝕量隨時(shí)間增長(zhǎng)是持續(xù)增大的，沖蝕速率的大小為沖蝕量隨時(shí)間變化的快慢，選取計(jì)算過(guò)程中沖蝕量隨時(shí)間變化最快的計(jì)算結(jié)果為最大沖蝕速率。FLUENT 軟件數(shù)值模擬直接計(jì)算獲得的沖蝕單位為沖蝕速率（kg/(m2·s)），為了便于同Birchenough 等人的實(shí)驗(yàn)結(jié)果對(duì)比，沖蝕速率的單位轉(zhuǎn)換方式如下：

ER（m·kg-1）=ER（kg·m-2·s-1）/ρw（kg·m-3）/mp（kg·s-1）

式中：mp為固體顆粒質(zhì)量流量。

將瞬態(tài)模型計(jì)算結(jié)果與 Salama 經(jīng)驗(yàn)?zāi)Ｐ蚚5]、Bourgoyne 經(jīng)驗(yàn)?zāi)Ｐ蚚6]以及簡(jiǎn)化CFD 模型計(jì)算結(jié)果做對(duì)比。

Salama 經(jīng)驗(yàn)?zāi)Ｐ停?/p>

式中：Sm為管道形狀系數(shù)，m2/s2，取2000；Vm為氣液混合流速，Vm=VSG+VSL，m/s；ρm為混合流體密度，kg/m3。

Bourgoyne 經(jīng)驗(yàn)?zāi)Ｐ停?/p>

式中：Fe為砂粒侵蝕系數(shù)，取為0.01g/kg；ρp為固體顆粒密度，kg/m3；ρw為管壁材料密度，kg/m3；Wp為固體顆粒體積流量，m3/s；Apipe為管道橫截面積，m2；HL為液體體積分?jǐn)?shù)。

簡(jiǎn)化CFD 沖蝕計(jì)算思路為：沖蝕簡(jiǎn)化計(jì)算方法假設(shè)固體顆粒均勻地分布在氣液兩相中，將分散泡狀流假定為混合均勻的單相流體，采用混合密度、混合黏度以及混合速度完成穩(wěn)態(tài)的沖蝕計(jì)算，見(jiàn)式（7）、（8）。該方法與Chen 等人[8]提出的分散泡狀流沖蝕求解思路一致。

式中：μg、μl、μm分別為氣體、液體、混合流體的黏度，Pa·s。

Birchenough 等人實(shí)驗(yàn)結(jié)果為4.6×10-9m/kg，采用沖蝕預(yù)測(cè)值/實(shí)驗(yàn)值來(lái)直觀地反映不同模型的計(jì)算精度，如圖3 所示。由圖3 可知，基于VOF 模型和DPM 模型的CFD 沖蝕計(jì)算結(jié)果最接近實(shí)驗(yàn)值。經(jīng)驗(yàn)?zāi)Ｐ涂芍庇^快速地確定沖蝕速率大小，但是不能確定管道內(nèi)部的流動(dòng)過(guò)程及沖蝕細(xì)節(jié)。簡(jiǎn)化CFD 計(jì)算模型由于將流體簡(jiǎn)化，其管道最終沖蝕分布并不能反映實(shí)際沖蝕情況。

圖3 不同沖蝕模型計(jì)算結(jié)果對(duì)比Fig.3 Comparison between calculation results by different erosion models

3.2 顆粒運(yùn)動(dòng)軌跡與氣液分布關(guān)系

利用截面含液率來(lái)表示某時(shí)刻管道截面氣液分布情況，不同時(shí)刻管內(nèi)固體顆粒瞬態(tài)分布與氣液兩相分布之間的關(guān)系如圖4 所示。由圖4 可知，管道內(nèi)部顆粒稀疏的區(qū)域（箭頭指向區(qū)域），液相含量也相對(duì)較少，氣相含量高，說(shuō)明顆粒與液相分布具有密切相關(guān)性。為深入研究顆粒分布與氣液兩相之間的關(guān)系，進(jìn)一步分析了管道三個(gè)監(jiān)控面上的截面含液率與相對(duì)顆粒含量之間的關(guān)系，如圖5 所示。截面含液率的峰值變化規(guī)律與相對(duì)顆粒含量的峰值變化趨勢(shì)比較接近，截面處含液率較高時(shí)（通過(guò)液體較多），截面的顆粒含量也相對(duì)較高，說(shuō)明管內(nèi)的固體顆粒大部分分散在液相之中。截面含液率和相對(duì)顆粒含量在不同截面的波動(dòng)范圍是不同的，如圖6、圖7 所示。

圖4 彎管不同時(shí)刻顆粒運(yùn)動(dòng)軌跡與管內(nèi)氣液分布Fig.4 Particle trajectories vs gas and liquid distributions in pipe bend at different time

由圖6 可知，上游直管段截面處的截面含液率波動(dòng)幅度小于彎頭處和下游直管段截面處。這主要是由于相對(duì)于彎頭以及下游直管段，上游直管段中氣液兩相流動(dòng)比較均勻，固體顆粒分布也較均勻。當(dāng)顆粒在彎頭處發(fā)生碰撞后，顆粒的運(yùn)動(dòng)軌跡發(fā)生改變，并且由于氣液兩相流經(jīng)彎頭后，分散泡狀流流型發(fā)生明顯變化，在彎頭頂部及下游直管段頂部液相含量明顯增大，而在彎頭及下游直管段底部，氣相含量明顯增大，氣泡與液體的摻混更加劇烈，導(dǎo)致截面顆粒含量的波動(dòng)較大。截面顆粒含量平均水平同樣也是下游直管段和彎頭處大于上游直管段，如圖7 所示。這主要是由于部分通過(guò)彎頭的顆粒在進(jìn)入下游直管段過(guò)程中會(huì)穿越彎頭頂部和下游直管段頂部的流體，而該區(qū)域的液相比例較高，固體顆粒進(jìn)入該區(qū)域后，與大量液體進(jìn)行交互計(jì)算，液相密度及黏度比氣相大得多。由公式（1）可知，固體顆粒進(jìn)入該區(qū)域后速度減小，形成滯止區(qū)，導(dǎo)致顆粒運(yùn)動(dòng)減緩，顆粒含量變大。

3.3 沖蝕分布與氣液分布關(guān)系

基于VOF 模型和DPM 模型的瞬態(tài)沖蝕計(jì)算獲得的沖蝕分布如圖8 所示。CFD 計(jì)算8 s 時(shí)，沖蝕量為2.24× 10-4kg/m2。瞬態(tài)仿真計(jì)算過(guò)程中，管道沖蝕量隨時(shí)間變化最快時(shí)的結(jié)果為最大沖蝕速率，求得最大沖蝕速率為4.41×10-9m/kg，沖蝕最嚴(yán)重位置出現(xiàn)在彎頭出口位置處。Birchenough 等人的實(shí)驗(yàn)結(jié)果為4.6× 10-9m/kg，差別很小。氣液兩相在彎頭處分布不均，在彎頭彎頭出口處（90°處）頂部（A 區(qū)域）液相含量較多，擁有較大速度的顆粒在碰撞管壁之前與液體相互作用，使得顆粒碰撞管壁速率減小，液體的存在對(duì)彎頭起到一定的保護(hù)作用，在彎頭小于90°角度處，部分位置氣相含量較高（B 區(qū)域，如60°、45°位置），顆粒與壁面碰撞時(shí)能量依然較大。雖然90°附近（A 區(qū)域）的顆粒碰撞能較大，但是損失也多，而小于90°區(qū)域（B 區(qū)域，如60°、45°位置）碰撞能略小，但是損失也較小，綜合以上因素，導(dǎo)致沖蝕嚴(yán)重區(qū)域較大。

圖5 不同截面處含液率與顆粒含量Fig.5 Liquid hold-up vs particle mass concentration at different sections: a) Section of upstream straight pipe in bend; b) 45° section of elbow; c) Section of downstream straight pipe in bend

4 結(jié)論

1）基于VOF 模型和DPM 模型的CFD 沖蝕模型計(jì)算結(jié)果與實(shí)驗(yàn)值接近，且與經(jīng)驗(yàn)?zāi)Ｐ图昂?jiǎn)化CFD模型相比，可直觀反映管道內(nèi)部流場(chǎng)及固體顆粒運(yùn)動(dòng)與壁面沖蝕之間的關(guān)系。

圖6 不同截面處截面含液率Fig.6 Liquid hold-up at different sections

圖7 不同截面處截面顆粒含量Fig.7 Particle mass concentration at different pipe sections

圖8 基于VOF 模型和DPM 模型耦合計(jì)算彎管沖蝕分布及氣液分布Fig.8 Erosion profile and gas-liquid distribution of the pipe bend based on the VOF and DPM coupled models

2）分散泡狀流中，管內(nèi)固體顆粒大部分分散在液相中，管道不同截面處的含液率與顆粒含量相關(guān)性較大，下游直管段和彎頭處，截面含液率較大，受到彎頭處流速滯止的影響，固體顆粒的含量也較大。

3）分散泡狀流中，固體顆粒在彎頭部位受到液相的緩沖作用以及氣液分布等影響，固體顆粒對(duì)管道沖蝕區(qū)域較大，沖蝕最嚴(yán)重位置出現(xiàn)在彎頭出口處附近。