高壓管匯沖蝕速率數值模擬新方法研究*

李建亭 曾 云 李 寧

(1.中石化中原石油工程有限公司工程服務管理中心 2.長江大學機械工程學院)

0 引 言

目前在石油工程領域，致密油和頁巖氣等非常規油氣資源的勘探與開發正愈發受到關注[1-2]。隨著非常規油氣田開發技術的日益發展，對鉆完井壓裂領域所需求的裝備性能要求也越來越高。壓裂液的輸送需要高壓、長壽命的高壓管匯，現有的高壓管匯受工作環境的制約[3]，壽命不長成為影響油氣開發經濟效益的關鍵因素。影響高壓管匯服役壽命的主要因素是沖蝕磨損。國內外學者對管匯沖蝕失效因素進行了多方面分析，壓裂液固相物理特性和運動特性是主要因素，管壁材料、內部流場結構、溫度及壓力等是次要因素[4-5]。

在考慮不同因素影響的情況下，目前關于沖蝕速率的計算還是以半經驗計算模型為主，各模型預測的結果都不是很精確[6-11]。數值模擬方法給沖蝕速率的計算提供了新的思路，尤其是計算流體力學離散單元法(CFD-DEM)的發展。以有限元和離散單元法為基礎，A.UZI等[12]提出了一維侵蝕模型(ODEM)的概念，描述了粒子壁面碰撞特性的一維流動和統計分布。使用單向耦合DPM模型進行侵蝕預測，忽略了粒子間碰撞的影響[13-15]。對顆粒體積分數較小的流體非常有效，但隨著粒子體積分數的增大，其準確性受到了質疑。離散元法(DEM)可以作為歐拉-拉格朗日模型的替代，能夠模擬粒子間的相互作用，但計算成本較高。在石油天然氣領域中，由于實際流場模擬尺度過大，采用離散元法(DEM)求解時計算工作量巨大，不太適合對流場進行模擬。文獻[16]基于Eulerian-Eulerian方法提出了混合顆粒流CFD模型，在侵蝕發生的位置使用拉格朗日模型來準確地跟蹤粒子并避免Eulerian-Eulerian方法的不確定性。這種模擬顆粒流的方法是歐拉-歐拉-拉格朗日方法，即DDPM模型，該模型介于DEM和歐拉-歐拉方法之間。目前，隨著超深井的開采，對壓裂液的要求逐漸上升，普通水基壓裂液的含砂體積分數一般在10%左右，而瓜膠壓裂液的攜砂性能好，含砂體積分數達到20%左右。文獻[17]已證實DDPM方法考慮四相耦合，其計算精確度優于DPM模型。

本文基于DDPM模型，考慮顆粒的不同物理特性，對高壓管匯不同結構彎頭進行沖蝕速率數值模擬分析，以評估高壓管匯的服役壽命。所得結論對高壓管匯的結構改進以及剩余壽命評估具有理論和工程應用意義。

1 數值模型

1.1 稠密離散相(DDPM)模型

基于標準歐拉多相流模型，本文研究的固相和液相沒有溫度傳遞，為非牛頓流體且流體不可壓縮流。連續性方程和動量方程如下[18]：

(1)

(2)

式中：α、ρ、v和p分別表示體積分數、密度、速度和壓力，下標f表示流體相。

體積分數和速度場根據拉格朗日方法求取，該方法通過將力平衡積分到每個顆粒上來計算每個顆粒的軌跡。顆粒運動軌跡控制方程如下：

(3)

式中：m為質量，u為速度，Cd為阻力系數，d為直徑，FKTGF為顆粒間碰撞和平移產生的力，下標s表示固體顆粒。

式(3)右側的前3項分別代表阻力、浮力和壓力梯度力。傳統的DPM方法也考慮了這些力。與DPM模型不同的是，該方法使用顆粒流動力學理論(KTGF)考慮粒子間碰撞和平移的影響。

為了進一步對阻力系數進行定量研究，本文設壓裂液支撐劑中存在球形和非球形顆粒，分別采用兩種阻力系數經驗模型進行分析。針對球形顆粒采用D.GIDASPOW[19]提出的模型，阻力系數Cd的表達式如下：

(4)

式中：Res為顆粒的雷諾數。

(5)

式中：μ為流體的黏度。

對于非球形粒子，采用A.HAIDER等[20]提出的模型，表達式如下：

(6)

式(6)中，各類系數如下：

b1=exp(2.328 8-6.458 1φ+2.448 6φ2)

(7)

b2=0.096 4+0.556 5φ

(8)

b3=exp(4.905-13.894 4φ+

18.422 2φ2-10.259 9φ3)

(9)

b4=exp(1.4681-12.258 4φ-

20.732 2φ2+15.885 5φ3)

(10)

形狀系數φ由等效球形顆粒表面積s與非球形顆粒表面積S之比進行求解。與DPM模型不同，DDPM模型考慮了顆粒之間的碰撞力。式(3)中的KTGF模型用于預測顆粒與顆粒間碰撞或者球形顆粒平移產生的應力，表達式如下：

(11)

(12)

μs=μs,col+μs,kin+μs,fr

(13)

(14)

(15)

式中：μs,col、μs,kin及μs,fr分別表示碰撞、運動和摩擦而產生的剪切黏度，ess為粒子碰撞產生的恢復系數，θs為顆粒溫度，g0為徑向分布函數。

本文研究對象中高壓管匯中的顆粒體積分數小于堆積極限，因此忽略摩擦黏度，且計算非球形顆粒時直接等效為顆粒粒徑ds。體積黏度λs表達式如下[21]：

(16)

1.2 顆粒之間碰撞受力計算模型

根據文獻[22]，顆粒受壓力計算的表達式如下：

(17)

顆粒溫度的概念基于顆粒的隨機平移運動產生的動能，恢復系數默認為0.9。本文主要研究顆粒之間碰撞對沖蝕的影響，其徑向分布函數是修正因子，當顆粒變得稠密也就是占比升高時，顆粒之間的碰撞會發生變化，其徑向分布函數與顆粒體積占比之間的關系如圖1所示。圖1中αs,max為顆粒體積分數極限值，取0.63。

圖1 徑向分布函數與顆粒體積分數之間的關系曲線

2 沖蝕計算模型

2.1 General模型

在沖蝕模型中，本文采用最為普遍的沖蝕速率模型，表達式如下：

(18)

2.2 DDPM沖蝕模型

稠密顆粒流動中，由于近壁面固體顆粒層的屏蔽作用，顆粒間的相互作用可能會影響沖蝕速率，且離散相顆粒即便是幾乎平行壁面也同樣會導致沖蝕磨損，這樣更接近工程實際。本文研究中總沖蝕率Etotal計算式如下：

Etotal=Eabrasive+Eimp

(19)

式中：Eabrasive為固相顆粒對壁面造成的沖蝕，Eimp為考慮屏蔽效應下的顆粒對壁面的沖蝕。

在稠密顆粒流動中，接近壁面的顆?？赡軙p速，或者被附近與壁面幾乎平行運動的其他固體顆粒反彈。這種屏蔽作用會降低顆粒撞擊固體壁面造成的沖蝕率，為考慮這種影響，引入屏蔽因子fshield,計算式如下：

(20)

當局部固相體積分數接近堆積極限(αs,max=0.63)時，壁面將受到保護，不受任何顆粒的沖擊，即沖蝕率為0。

3 數值模擬結果分析與討論

3.1 幾何模型和網格劃分

結合壓裂現場實際工況，高壓管匯的內徑D=76 mm，進口長度L1=7D、出口長度L2=10D。彎徑比一共設置4組(R/D=1.0、1.5、2.0及2.5)，目的是研究不同結構下顆粒軌跡。邊界條件選擇速度進口和壓力出口。

圖2為高壓管匯幾何模型及網格劃分圖。圖2中所有網格采用六面體結構網格，網格尺寸已經過無關性驗證。

圖2 高壓管匯幾何模型及網格劃分圖

本文結合壓裂工況不同排量計算，設置支撐劑顆粒密度為2 750 kg/m3，連續相流體選擇不可壓縮流體，并結合不同工況選用DPM模型和DDPM模型分別進行數值模擬，以便對比分析。表1為不同工況下的參數取值。

表1 不同工況下的參數取值

3.2 不同模型下流場特性分析

在沖蝕磨損研究中，影響因素主要有以下5個方面：①粒子的物理特性和軌跡，②局部的流場和湍流，③固體壁面條件，④多相流的影響，⑤材料磨損造成局部的空穴。在本文的研究中，因素③和因素④在沖蝕模型中已經進行了定義。因此，結合不同工況下定義的參數，主要針對高壓管匯的流場特性和顆粒軌跡進行分析。

流場分析中發現，高壓管匯的沖蝕失效主要發生在彎頭處，因此本研究針對彎頭不同偏轉角(0°、15°、30°、45°、60°、90°)下彎頭截面速度展開分析。在流速12 m/s、質量流量5.14 kg/s、R/D=1.5的工況下，高壓管匯彎頭處截面云圖如圖3所示。

由圖3可知，在彎頭偏轉角0°～60°區域可以明顯看到流動分層現象，結合伯努利方程可知，彎頭內側壓力小于彎頭外側壓力，從而形成高速流動區域。在彎頭偏轉角60°～75°區域可以發現一個現象，即彎頭靠近內側出現了局部流速較小區域，這是一種二次流動現象，隨著角度繼續增大接近90°，明顯出現了低速區域，然后在90°截面后出口直管段，速度在內壁與外壁發生了明顯變化。出口直管端靠近外側的流速明顯大于內側流速，這也說明從彎頭區域開始，整個高壓管匯以外側出現沖蝕為主。從湍動能和湍流耗散率的角度來看，在不同的偏轉角截面中，靠近彎頭內壁的湍動能逐漸增強，湍流耗散率也逐漸增大。這也充分反映在這個區域流動出現了拐點，形成了二次流動。

圖3 高壓管匯彎頭處截面云圖

3.3 不同模型下的沖蝕速率對比分析

采用DPM模型和DDPM模型對高壓管匯4種工況、4種結構設計開展正交試驗，共16組，其數值仿真結果如圖4所示。

圖4 不同模型下高壓管匯沖蝕速率模擬曲線

由圖4a可知：隨著彎徑比的增大，沖蝕速率呈遞減趨勢，彎徑比在1.5～2.0時沖蝕速率出現一個相差數量級的下跌，這可能是因為DPM模型計算顆粒運動時主要圍繞流體的運動軌跡；隨著彎徑比的增大，在彎頭區域流動比較充分，顆粒在一定程度上接近彎頭外側壁面區域時不會由于速度過快而匯聚集中沖擊壁面，這也表明DPM計算方法存在一定缺陷，并未考慮顆粒與顆粒之間的碰撞。由圖4b可知，采用DDPM模型時，彎徑比在1.0～2.5時，沖蝕速率先減小后上升，R/D=1.5時最小，與DPM模型稍有不同，結合引言中DPM模型的劣勢，說明DPM模型在模擬固相體積分數較高的液固兩相流時存在一定的偏差。

選擇4組工況中流速和顆粒質量流量最小的一組沖蝕云圖進行分析，如圖5所示。由圖5可知：不論是DPM模型還是DDPM模型，在彎徑比為1.5和2.0時沖蝕集中位置沒有另外幾種彎徑比明顯，且造成的沖蝕面積沒有彎徑比為1.0和2.5時大，這是因為彎徑比較小時，壓裂液流經彎頭處，其流動方向由于結構變化急促，粒子更容易沖出液相，沖擊內壁，造成更為嚴重的沖蝕磨損；彎徑比變大時，在彎頭處運動方向改變較為平緩，慣性相對較小，對壁面的沖擊較輕，但由于DDPM模型考慮顆粒之間的碰撞，在彎徑比過大的區域，匯聚的顆粒增多，顆粒之間的碰撞導致沖蝕速率反而會上升，但隨著彎徑比繼續增大，沖蝕速率呈遞減趨勢。

圖5 不同模型下管匯沖蝕云圖對比

3.4 顆粒形狀系數對沖蝕速率和沖蝕區域的影響

在高壓管匯的彎頭處容易形成沖蝕的區域主要為彎頭處的外壁和出口直管端的側壁。圖6為4種不同結構下顆粒軌跡圖。圖6中的流速為12 m/s，質量流量為5.14 kg/s。4種不同結構的顆粒軌跡圖皆表明兩處區域(彎頭外壁和出口直管端側壁)是顆粒軌跡匯集之處。產生這種現象的原因是：流體從直管進入彎頭時存在局部阻力，導致水力動能減小，且3.2節分析過彎頭處外壁壓力大于內壁壓力，內側壁形成高速流動區域，外側壁流速相對較低，但外壁區域固相顆粒所受離心力大，導致沖蝕磨損明顯增大；在彎頭出口直管端由于固相顆粒撞擊壁面反彈增多，再加上二次流的牽引作用，導致出口直管側壁也容易形成沖蝕嚴重區域。

圖6 4種不同結構下顆粒軌跡圖

為了進一步研究固相顆粒的物理特性，分別設置0.9、0.7、0.5和0.3共4種不同的形狀系數來分析顆粒形狀系數對沖蝕速率和沖蝕區域的影響，結果如圖7所示。由圖7可知，隨著顆粒形狀系數的減小(顆粒逐漸變得尖銳)，彎頭部分的最大沖蝕速率呈遞增趨勢，本工況速度10 m/s，質量流量4.17 kg/s，彎徑比1.5。顆粒的形狀系數從1.0減小至0.3時，最大沖蝕速率從5.726×10-5kg/(m2·s)增大到7.79×10-5kg/(m2·s)。由圖7可以發現，在彎頭截面偏轉角為86°～88°時，也就是靠近彎頭和出口端直管的交匯處沖蝕速率最大。由圖7可以發現一個現象，當固相顆粒形狀系數大于0.5時，沖蝕速率隨著彎頭截面偏轉角的增大而增大，但當固相顆粒形狀系數小于0.5時，在偏轉角15°～45°之間，沖蝕速率呈現遞減的趨勢。其原因可能是偏轉角15°～45°這段區域的外壁內壓大，且由于形狀系數小固相顆粒的拖曳力系數大，顆粒與流體的跟隨性很好，顆粒被流體帶著繞過障礙物，而不是撞擊壁面，顆粒與壁面的碰撞概率減小，產生了一定的撞擊遲緩效應。

圖7 顆粒形狀系數對沖蝕速率和沖蝕區域的影響

3.5 顆粒斯托克斯系數對沖蝕速率的影響

高壓管匯輸送的壓裂液具有高密度和高黏度(曳力系數考慮黏度的影響)等特性，可對固相顆粒的運動以及沖蝕產生一定影響。與顆粒氣流相比，當顆粒體積分數較高時，顆粒載荷對流場的影響通常變得很重要，為此引入一個無因次量——顆粒斯托克斯系數(St)來研究固相顆粒對高壓管匯的沖蝕行為，以探索流場與固相顆粒特性之間的復雜耦合效應。St定義為顆粒響應時間與流體流動時間的比值，表示如下：

(21)

當St>1，粒子對時間的響應遠不及流體傳播對時間的響應，因此有足夠的時間來存儲粒子以響應流場的變化，并且它們可以緊隨流體的流動；當St<1時，情況恰恰相反，粒子將獨立于流體而運動。為了研究顆粒的斯托克斯系數對沖蝕的影響，設置固相顆粒為單一直徑(0.25 mm)，質量流量為0.7 kg/s，流速為3、6、9和12 m/s。根據式(21)計算得St分別為0.413 7、0.827 4、1.241 1和1.654 8。圖8為不同St下沖蝕率云圖。

圖8 不同顆粒斯托克斯系數下沖蝕率云圖

由圖8可知：當St<1時，主要沖蝕區域在彎頭靠近出口直管端兩處側壁，這是因為一部分顆粒隨二次流作用對側壁造成了沖擊磨損；當St>1，二次流對顆粒的影響作用逐漸減弱，固相顆粒主要受慣性力，顆粒動量能克服湍流渦團的影響，流場對顆粒的作用減弱，所以主要沖蝕區域集中在彎頭外壁。圖9為不同St下顆粒軌跡、流場軌跡及速度矢量對比圖。從圖9可以看到，當固相顆粒斯托克斯系數小于1時，顆粒軌跡和流場軌跡明顯不同，固相顆粒在彎頭出口靠近側壁區域隨著二次流形成了復雜的沖擊和反彈軌跡，這也是彎頭出口直管端兩處側壁區域形成主要沖蝕的根本原因。

圖9 不同顆粒斯托克斯系數下顆粒軌跡、流場軌跡線及速度矢量對比圖

無綱量用于描述顆粒與流體之間的相互作用，在實際情況中，彎管中顆粒的體積分數大于0.01，而DPM模型忽略了顆粒之間碰撞的假設不是很合理。圖10是不同St下2種模型計算的沖蝕速率對比圖。

圖10 不同固相顆粒斯托克斯系數下2種模型計算的沖蝕速率對比圖

從圖10可以看出，當斯托克斯系數接近1時，DDPM模型計算的沖蝕率出現一個陡峭的上升趨勢。DDPM模型的計算結果存在一個臨界的顆粒載荷，對于這個拐點的存在，改變管徑、顆粒粒徑和流體流速皆可改變固相的斯托克斯系數，從而間接改變沖蝕區域，這對高壓管匯易沖蝕區域的防護措施制定具有一定的指導意義。

4 結 論

本文使用CFD方法，分別采用DPM模型和DDPM模型，分析了彎徑比、壓裂工況參數、固相顆粒物理特性及顆粒斯托克斯系數對高壓管匯沖蝕速率及沖蝕區域的影響，得到以下結論：

(1)高壓彎管在60°～75°偏轉角截面區域逐漸出現流動分層現象，二次流動現象出現，湍動能逐漸增強，湍流耗散率也逐漸增大。該流場模擬結果間接解釋了彎頭出口靠近直管端兩側壁面出現沖蝕區域的現象。

(2)實際工況中，顆粒在彎頭處會匯聚，其體積分數達到20%，甚至更高。DDPM模型考慮了四相耦合，考慮了顆粒之間的碰撞和作用力，在模擬顆粒軌跡上要優于DPM模型。

(3)最主要的沖蝕區域還是彎頭出口外壁區域，顆粒形狀系數小于0.5(尖銳顆粒)時，在彎頭15°～45°偏轉截面沖蝕減輕。影響沖蝕區域的因素是彎徑比和顆粒粒徑。

(4)固相顆粒斯托克斯系數主要與顆粒粒徑、流體流速及管道直徑有關，在St>1與St<1時顆粒軌跡出現明顯差別，沖蝕區域出現一定改變。該理論對實際工況中不同沖蝕區域的防護措施制定具有一定的指導意義。

(5)DDPM模型計算結果表明存在一個臨界的顆粒載荷，考慮顆粒斯托克斯系數，顆粒體積分數對沖蝕速率具有重要影響。