豎直彎管中含固體顆粒分散泡狀流沖蝕預測

彭文山，馬力，侯健，曹學文，韓明一

（1.中國船舶重工集團公司第七二五研究所 海洋腐蝕與防護重點實驗室，山東 青島 266237；2.中國石油大學（華東）儲運與建筑工程學院，山東 青島 266580；3.廊坊中油朗威工程項目管理有限公司，河北 廊坊 065000）

隨著我國油氣開采由陸地向海洋轉移，海底管道投用量持續增加[1]，海洋混輸管網得到不斷應用。混輸管道中流體流動復雜，流體對管壁沖刷嚴重，油井出砂量及原油含水量隨著油氣田開采時間的增長而增加。油氣開采及運輸過程中，管內固體顆粒隨油氣流動，極大地加劇了管內介質對管壁的沖擊，極易造成管道沖蝕破壞，導致危險事故[2-3]。因此，提出準確的多相流沖蝕計算方法對于預測管道可能出現的沖蝕最嚴重位置，進而保障管輸安全、減少經濟損失意義重大。

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

1 分散泡狀流沖蝕求解

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

1.1 連續相流場計算

標準k-ε 模型在計算管道內部不同截面處流體速度的誤差較小[10]，選用該模型計算管道內部流場。使用VOF 模型實現管內分散泡狀流的求解。另外，考慮表面張力的作用，采用Brackbill 等人[11]提出的連續表面力模型（CSF）。利用這個模型，VOF 模型中的附加表面張力通過源項的方式添加到動量方程中。

1.2 離散相控制方程

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

1.3 沖蝕計算模型

目前已有許多沖蝕模型被應用于沖蝕計算[13]，Liu等人[14]對沖蝕計算研究發現，廣島大學Oka 等人[15-16]提出的沖蝕模型能較為準確地預測多相流條件下管道的沖蝕速率。另外，模型中考慮了較多的影響因素，適應范圍較廣。使用該模型進行沖蝕計算：

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

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

2 仿真模型建立

2.1 管道基本參數

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

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

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

2.2 網格劃分

為提高VOF 模型與DPM 模型耦合計算的穩定性，計算模型采用規則四邊形面網格，采用六面體彎管體網格，網格劃分如圖2 所示，豎直彎管三維模型網格總數為524 000。

圖2 豎直彎管網格劃分Fig.2 Computational mesh of the vertical pipe bend

2.3 邊界條件

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

式中：en為法向恢復系數；et為切向恢復系數。

2.4 數值求解

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

3 結果分析

3.1 沖蝕結果對比

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

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

式中：mp為固體顆粒質量流量。

將瞬態模型計算結果與 Salama 經驗模型[5]、Bourgoyne 經驗模型[6]以及簡化CFD 模型計算結果做對比。

Salama 經驗模型：

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

Bourgoyne 經驗模型：

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

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

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

Birchenough 等人實驗結果為4.6×10-9m/kg，采用沖蝕預測值/實驗值來直觀地反映不同模型的計算精度，如圖3 所示。由圖3 可知，基于VOF 模型和DPM 模型的CFD 沖蝕計算結果最接近實驗值。經驗模型可直觀快速地確定沖蝕速率大小，但是不能確定管道內部的流動過程及沖蝕細節。簡化CFD 計算模型由于將流體簡化，其管道最終沖蝕分布并不能反映實際沖蝕情況。

圖3 不同沖蝕模型計算結果對比Fig.3 Comparison between calculation results by different erosion models

3.2 顆粒運動軌跡與氣液分布關系

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

圖4 彎管不同時刻顆粒運動軌跡與管內氣液分布Fig.4 Particle trajectories vs gas and liquid distributions in pipe bend at different time

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

3.3 沖蝕分布與氣液分布關系

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

圖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 結論

1）基于VOF 模型和DPM 模型的CFD 沖蝕模型計算結果與實驗值接近，且與經驗模型及簡化CFD模型相比，可直觀反映管道內部流場及固體顆粒運動與壁面沖蝕之間的關系。

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

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

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

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

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