王明,劉巨保,王雪飛,姚利明,張西偉,張志龍,楊明
(1.東北石油大學 機械科學與工程學院,黑龍江 大慶 163318;2.新加坡南洋理工大學 機械與航天工程學院,新加坡 639798;3.海洋石油工程股份有限公司 安裝事業部,天津 300451)
沖蝕磨損是移動的固體顆粒撞擊表面,將材料從固體表面除去的過程[1]。壓裂施工作業中,顆粒從流體中獲得動能,隨流體運動,沖擊連接管匯及壓裂工具,導致管壁材料損耗,壁厚變薄,管匯強度降低,甚至破裂穿孔[2-3]。早期進行了一些典型的試驗來研究侵蝕機理[4-7],但沖蝕是一個十分復雜的過程,流動條件、顆粒材料屬性、幾何形狀等因素均可對沖蝕結果產生影響[8]。單純依靠試驗,難以實現復雜因素的沖蝕影響研究。隨著計算機科學的發展,數值模擬逐漸成為復雜因素對沖蝕影響研究的重要手段。為了準確預測沖蝕過程,建立了兩種有效的顆粒尺度模型——計算流體動力學-離散相模型(CFD-DPM)[9]和計算流體動力學-離散元方法(CFD-DEM)[10]。
CFD-DPM 方法可以考慮顆粒-流體相互作用的雙向耦合,也可以不考慮顆粒對流體作用的單相耦合,較適合用于不考慮顆粒間相互作用的稀疏流[9]。Wang[11]建立了長半徑彎頭沖蝕模型。李昭乾[12]研究了管道外徑、彎曲半徑、顆粒直徑與顆粒流速等參數對90°彎管沖蝕率的影響。Edwards[13]和張孟昀[14]對90°彎管和盲通管的沖蝕特性進行了研究。Mazumder[15]研究了U 型彎道內液體和氣體速度對最大沖蝕位置的影響。楊德成[16]研究了固體顆粒參數與彎管結構對沖蝕位置和沖蝕速率的影響。
CFD-DEM 方法可以跟蹤單個顆粒運動,考慮顆粒間相互作用,能夠更準確地預測顆粒運動、反彈及二次流對沖蝕的影響,稠密流和稀疏流均適用。Uuemis 和Kleis[17]研究發現,隨著顆粒濃度的增加,絕對沖蝕率(總侵蝕質量)隨著濃度的增加而增加,但比沖蝕率(每單位質量的顆粒撞擊壁面時從表面除去的物質質量)減少。Chen[18]研究發現,在高濃度下,顆粒從壁面反彈時,它們會撞擊向壁面移動的顆粒,使其減速,出現“屏蔽”或“緩沖”現象,并不總是較高的顆粒濃度就會出現較高的沖蝕速率。上述研究表明,在進行沖蝕研究時,不能忽略顆粒-顆粒的相互作用。Markus[19]采用CFD-DEM 方法,對顆粒沖蝕進料管進行了數值模擬,局部高沖蝕率的計算結果與試驗結果一致。姚利明[20]采用CFD-DEM 方法對縮徑管進行了顆粒的碰撞和沖蝕特性研究。Jukai[21]采用CFD-DEM 方法,對直徑為40 mm 且彎曲角度為60°、45°的彎管沖蝕進行了預測。研究結果表明,CFD-DEM 方法可以有效地對沖蝕磨損和顆粒流體流動進行數值模擬。
上述研究中,在沖蝕磨損影響因素方面,多采用特定因素對沖蝕結果的影響進行研究(控制變量方法只改變單一參數),未對沖蝕磨損影響因素間的交互作用進行討論,且單因素分析計算工作量巨大。因此,本文采用CFD-DEM 計算方法,結合正交試驗設計,分析各因素對壓裂管匯及井下噴砂器等工具中常見的變徑T 型流道沖蝕的影響,力求更加全面、準確、高效地分析變徑T 型流道在壓裂施工中的沖蝕磨損規律,為后續工程應用提供理論支撐。
計算流體力學與離散單元法耦合的 CFD-DEM方法,其連續流體的流動是基于CFD 中局部平均概念求解Navier-Stokes 方程獲得,離散顆粒的運動是基于DEM 中求解牛頓第二運動定律獲得[22-24]。CFDDEM 方法能夠準確地描述每個計算時間步長,顆粒-流體、顆粒-顆粒、顆粒-壁面間的動量交換及碰撞時能量的傳遞[25]。顆粒在流體內部受到流體的作用力和顆粒與顆粒、顆粒與壁面的碰撞力,顆粒的存在使流體的流動狀態發生改變。
取管道內流體與顆粒群為研究對象,建立圓管內流體與顆粒群耦合的兩相流模型,如圖1 所示。該模型共分2 個域,顆粒群域Ωp和流體域Ωf。在圓管內兩相流體系中任取一控制體dVc,dVf、dVp為流體和顆粒在控制體內占據的體積。

圖1 流體與顆粒群耦合模型示意圖Fig.1 Schematic diagram of fluid and particle group coupling model
1.1.1 流體控制方程
流體控制方程為[26]:

式中:ρf為流體密度;uf為流體速度矢量;αf為流體的體積分數;Sf為流體與顆粒動量交換源項。
固液兩相流體與單向流體的區別在于,顆粒的加入占據了流體計算單元的體積,相應流體的體積分數αf發生改變,見式(3)。

動量交換源項是作用在流體計算單元內的流體阻力(曳力)總和的體積平均項。阻力(曳力)是由于兩相之間的相對運動產生的,通過計算阻力的動量交換源項可以實現兩相之間的耦合,源項Sf的計算公式為:

式中:Fd為顆粒受到的阻力(曳力)。
1.1.2 顆粒動力學方程
顆粒的平移和旋轉運動,根據牛頓第二定律,其計算公式為:

式中:mi為顆粒i質量;Fcn,ij和Fct,ij為顆粒i與顆粒j的法向和切向碰撞力;Ii為顆粒i的慣性矩;ω i為顆粒的角速度;Tc為顆粒碰撞引起的力矩;Fg為顆粒重力;Fvm為虛擬質量力;Fpg為壓力梯度力;Fsl為Saffiman 升力;Fb為Basset 力;Fm為Magnus 力。
文獻[27]認為流體作用在顆粒上的力主要是阻力,因此本文只關注重力和阻力。采用廣泛應用的Di Felice阻力模型[28]來描述流體與顆粒之間的相互作用。模型公式為:

顆粒在流體內部受到流體的作用力,顆粒的存在改變流場的流動狀態,使作用在顆粒上的力實時更新。顆粒對流體的影響通過將動量交換源項Sf和體積分數αf引入流體控制方程(1)和方程(2)中,流體對顆粒的影響通過將阻力Fd引入到顆粒動力學方程(5)中,實現顆粒-流體雙向耦合。
運動顆粒與壁面碰撞是造成壁面沖蝕磨損的主要因素,因此提高顆粒碰撞搜索精度是準確預測沖蝕磨損的關鍵。顆粒-壁面碰撞簡化為顆粒與一個大直徑顆粒發生碰撞。顆粒碰撞搜索過程如下:
1)計算域內選取搜索碰撞的顆粒為目標顆粒,以目標顆粒一個計算時間步長Δt的運動位移為立方體對角線,建立空間搜索網格體,如圖2a 所示。由于計算域內不同顆粒在Δt時間步長內運動的位移不同,所構建的空間網格也會出現不同尺寸,即為變網格。
2)對在計算時間步長內運動軌跡與空間網格體存在交集的顆粒進行篩選,若運動軌跡與空間網格體存在交集,則可能與目標顆粒發生碰撞,若不存在交集,則不可能與目標顆粒發生碰撞,如圖2a 所示。

圖2 顆粒碰撞搜索示意圖Fig.2 Search diagram of particle collision: (a) schematic diagram of variable grid establishment(c) schematic diagram of particle spacing screening


顆粒碰撞采用軟球模型,如圖3 所示。顆粒間接觸過程簡化為彈簧振子的阻尼振動,法向碰撞力Fcn,ij是彈簧和法向阻尼作用在顆粒i上的彈性力和阻尼作用的合力。根據Hertz 接觸理論,顆粒碰撞力表示為:

圖3 顆粒碰撞示意圖Fig.3 Schematic diagram of particle collision


式中:μ為靜摩擦因數;k為切向單位向量。
顆粒位置變化直接導致流體體積分數和動量交換力源項的變化,因此根據流體單元內流體體積分數和動量交換力源項的變化,建立流體與顆粒群耦合的收斂條件,如式(13)—(14)。

CFD-DEM 耦合計算流程如圖4 所示,CFD 計算采用SIMPLE 算法求解流體域連續性方程和動量方程,根據顆粒和流體的相對速度,計算顆粒受到的流體作用力Ff,傳遞給DEM 求解器,并啟動DEM 計算。DEM 采用牛頓第二運動定律求解顆粒運動方程,根據顆粒在Δtp的運動軌跡,搜索和判斷顆粒碰撞,若顆粒產生碰撞,需修正顆粒計算時間步長。根據顆粒在計算域內位置和受力的變化,計算得到流體的體積分數αf和動量交換力源項Sf,若不滿足式(13)和式(14),傳遞給CFD 求解器,自動更新流體計算時間步長,并推進CFD 計算。

圖4 CFD-DEM 耦合計算流程圖Fig.4 Flow chart of CFD-DEM coupling calculation


隨著對沖蝕規律的不斷研究,得到了許多預測侵蝕的經驗公式。本文以液固兩相流體動力學為基礎,將T 型流道的幾何特征和材料特性考慮在內,沖蝕模型采用Tulsa 大學提出的半經驗沖蝕模型[29]:

基于CFD 和DEM 平臺,采用Delphi 編程二次開發了流體與顆粒群耦合動力學分析程序。
采用本文建立的耦合計算模型與文獻[20]的試驗結果進行對比。試驗采用大徑為38 mm、小徑為12 mm 的縮徑管為研究對象。模擬條件完全來源于試驗,流體密度為1000 kg/m3,黏度為180 MPa·s,速度為11.8 m/s,顆粒直徑為0.65 mm,砂比為35%,管道材料彈性模量為2.06×1011Pa,密度為7850 kg/m3,泊松比為0.3。試驗沖蝕結果與數值模擬沖蝕結果對比圖如圖5 所示。試驗沖蝕深度與模擬沖蝕深度對比曲線如圖6 所示。由圖5 和圖6 可知,模擬與試驗結果吻合,驗證了本文所建立的數值方法的可靠性。

圖5 試驗沖蝕結果與數值模擬結果對比圖Fig.5 Comparison of experimental erosion results and numerical simulation results

圖6 試驗沖蝕深度與模擬沖蝕深度對比曲線Fig.6 Comparison curve of experimental erosion depth and simulated erosion depth
變徑T 型流道支流道直徑D=38 mm,主流道為變徑段,變徑比s取0.4、0.6、0.8、1.0。變徑T 型流道幾何模型如圖7 所示。

圖7 變徑T 型流道結構示意圖Fig.7 Schematic diagram of variable diameter T-shaped runner structure
根據葡扶179-斜89 井壓裂施工方案,顆粒材質選用石英砂,壓裂施工使用100 目、40—70 目、20—40 目3 種規格支撐劑,支撐劑粒徑取0.125、0.3、0.6、0.85 mm,彈性模量為2.1×108Pa,泊松比為0.2,密度為2650 kg/m3。根據石油與天然氣行業標準ST/Y 6376—2008《壓裂液通用技術條件》,確定壓裂液黏度為180 MPa·s。流道壁面材料為鋼,材料密度為7870 kg/m3,泊松比為0.28,彈性模量為2.12×1011Pa。
流體采用SIMPLEC 方法和瞬態解算器求解連續相流的控制方程,采用標準k-?湍流模型和近壁區處理的標準壁函數。速度入口邊界,流體與顆粒速度相同,出口采用壓力出口。
顆粒-顆粒或顆粒-壁面接觸采用 Hertz-Mindlin(無滑移)接觸模型,其揭示了單個顆粒的運動特性,并跟蹤了顆粒與顆粒和顆粒與壁面的碰撞。顆粒與壁面碰撞的接觸參數,通過虛擬試驗方法進行標定,恢復系數為0.4,靜摩擦因數為0.68,動摩擦因數為0.01。
采用ANSYS ICEM CFD 軟件對變徑T 型流道進行六面體網格劃分,局部進行網格加密。以流道內平均壓差值作為指標,進行網格無關性驗證。以變徑比s=0.4 的T 型流道為例,如圖8 所示,網格大于40萬以上,平均壓差值變化平穩,可以保證計算精度。

圖8 網格無關性驗證Fig.8 Grid independence verification
本文選用入口流速、砂比、支撐劑粒徑、變徑比4 個參數作為影響因素,對變徑T 型流道的沖蝕進行數值仿真。正交試驗因素與水平劃分見表1。仿真試驗結果如表2 所示。

表1 沖蝕影響因素-水平劃分Tab.1 Erosion influencing factors-level division

表2 沖蝕磨損多參數分析正交試驗表Tab.2 Orthogonal test table for multi-parameter analysis of erosion and wear
對表2 進行分析,得到徑向流速、切向流速、管壁壓力、最大沖蝕速率均值主效應,如圖9—12 所示。其中因素變化對均值主效應的響應優先排序見表3。

表3 因素變化各均值主效應優先排序表Tab.3 Priority ranking table of main effects of each mean value of factor changes

圖9 徑向流速主效應圖Fig.9 Radial velocity main effect diagram
以徑向流速作為衡量變徑T 型流道沖蝕嚴重度的指標時,各因素的影響由大到小依次為入口流速、支撐劑粒徑、變徑比、砂比,其值分別為25 m、0.6 mm、1、15%時,徑向流速均值主效應最大,沖蝕最嚴重。
以切向流速作為衡量變徑T 型流道沖蝕嚴重度的指標時,各因素的影響由大到小依次為入口流速、變徑比、砂比及支撐劑粒徑,其值分別為25 m/s、1、5%、0.125 mm,切向流速均值主效應最大,沖蝕效果最弱。
以管壁壓力作為衡量變徑T 型流道沖蝕嚴重度的指標時,各因素的影響由大到小依次為入口流速、支撐劑粒徑、變徑比、砂比,其值分別為25 m/s、0.85 mm、0.6、25%,管壁壓力均值主效應最大,沖蝕最嚴重。

圖10 切向流速主效應圖Fig.10 Main effect diagram of tangential velocity

圖11 管壁壓力主效應圖Fig.11 Main effect diagram of pipe wall pressure

圖12 最大沖蝕速率主效應圖Fig.12 Main effect diagram of maximum erosion rate
以最大沖蝕速率作為衡量變徑T 型流道沖蝕嚴重度的指標時,各因素的影響由大到小依次為入口流速、支撐劑粒徑、變徑比、砂比,其值分別為25 m/s、0.6 mm、0.4、35%,最大沖蝕速率均值主效應最大,沖蝕最嚴重。
綜上所述,以徑向流速作為衡量沖蝕程度的指標時,入口流速影響最大,其他因素影響不顯著。以切向流速作為衡量沖蝕程度的指標時,入口流速產生的影響最大,其次是變徑比,砂比和支撐劑粒徑影響不顯著。以管壁壓力和最大沖蝕速率作為指標時,入口流速影響最大,其他因素影響不顯著。以最大沖蝕速率作為衡量沖蝕程度的指標時,入口流速影響最大,其次是支撐劑粒徑、變徑比、砂比。
由圖13—16 可知,入口流速與砂比、支撐劑粒徑、變徑比各因素交互作用,對四項試驗指標產生的影響最為顯著,各水平整體變化趨勢基本一致。砂比、支撐劑粒徑、變徑比三因素兩兩交互作用對試驗指標存在一定影響,但效果沒有入口流速明顯。下文采用單因素分析,進一步驗證上述結論。

圖13 徑向流速交互作用圖Fig.13 Radial velocity interaction diagram

圖14 切向流速交互作用圖Fig.14 Interaction diagram of tangential flow velocity

圖15 管壁壓力交互作用圖Fig.15 Diagram of the interaction of pipe wall pressure

圖16 最大沖蝕速率交互作用圖Fig.16 Interaction diagram of maximum erosion rate
通過控制單一變量的方式分析各因素對變徑T型流道的沖蝕影響,模擬參數取第4 水平。
砂比35%、支撐劑粒徑0.85 mm、變徑比1 不變,入口流速為10、15、20、25 m/s 時,變徑T 型流道沖蝕速率對比如圖17 所示,內部顆粒堆積狀態對比如圖18 所示,最大沖蝕速率隨入口流速的變化曲線如圖19 所示,沖蝕速率隨位置的變化規律如圖20 所示。

圖17 入口流速不同時變徑T 型流道沖蝕速率對比圖Fig.17 Comparison of erosion rates of variable diameter T-shaped runners with different inlet velocities

圖18 入口流速不同時變徑T 型流道內部顆粒堆積狀態對比圖Fig.18 Comparison of particle accumulation conditions in the variable diameter T-shaped flow channel with different inlet flow rates

圖19 最大沖蝕速率隨入口流速的變化Fig.19 Curve of maximum erosion rate versus inlet flow velocity

圖20 不同流速最大沖蝕速率隨位置的變化規律Fig.20 The curve of the maximum erosion rate with different flow velocity changes with position
由圖17 可知,流道最大沖蝕速率和沖蝕面積隨入口流速的增加而增大,最大沖蝕速率位置在主流道與支流道交界處,形成沖蝕圓環,隨著入口流速的增加,圓環面積逐漸擴大。由圖17 和圖18 可知,兩相流體通過T 型流道的主流道后,流體向支流道出口處流動,形成分叉,在迎流面形成低流速區域,顆粒在此處形成堆積,堆積顆粒對底部壁面形成保護。隨入口流速的增大,流體對顆粒的攜帶能力增強,堆積顆粒寬度從35 mm 減小到23 mm,沖蝕圓環中心面積逐漸減小。兩相流體在堆積顆粒周邊與支流道壁面碰撞,形成沖蝕圓環。入口流速越大,支撐劑具有的動能越大,對壁面碰撞的沖蝕速率越大,反彈后的顆粒還具有足夠的動能對壁面造成沖蝕,所以沖蝕面積隨入口流速的增大而增大。沖蝕速率變化趨勢與主效應圖中,入口流速的變化趨勢相同。由圖19 可知,最大沖蝕速率隨入口流速的增大而增大。
由圖20 可知,隨入口流速的增大,最大沖蝕速率位置基本不變,在距離主流道中心30 mm 附近。同一入口流速下,沿支流道出口方向,沖蝕速率逐漸減小,入口流速越大,沖蝕速率的變化速率越快,沖蝕區域逐漸向支流道兩側出口處擴大。
入口流速25 m/s、支撐劑粒徑0.85 mm、變徑比1 不變,砂比為5%、15%、25%、35%時,變徑T 型流道沖蝕速率對比如圖21 所示,最大沖蝕速率隨砂比的變化曲線如圖22 所示。

圖21 砂比不同變徑T 型流道沖蝕速率對比圖Fig.21 Comparison of erosion rates of T-shaped runners with different sand ratios

圖22 最大沖蝕速率隨砂比變化曲線Fig.22 Curve of maximum erosion rate with sand ratio
砂比為5%~35%時,隨著砂比的增大,變徑T 型流道的最大沖蝕速率和沖蝕面積均增大。主要原因是顆粒直徑一定,砂比增大,顆粒數量增多,單位時間內撞擊流道次數增多,同時砂比增大后,增大了顆粒與顆粒、顆粒與壁面的碰撞次數和碰撞范圍。因此,砂比增大會使流道壁面的平均沖蝕率增大。最大沖蝕速率位置基本不變。
在本文中,顆粒砂比范圍為5%~35%,流道壁面最大沖蝕速度從2.81×10?4mm/s 增加到5.33×10?4mm/s,可以看出,砂比對最大沖蝕速率的影響沒有入口流速那樣顯著。最大沖蝕速率的變化趨勢與砂比變化正相關。本文中砂比取值過小,未能全面反映最大沖蝕速率與砂比之間的關系。
入口流速25 m/s、砂比35%、變徑比1 不變,支撐劑粒徑為0.125、0.3、0.6、0.85 mm,得到最大沖蝕速率隨支撐劑粒徑的變化曲線,如圖23 所示。由圖23 可知,支撐劑粒徑在0.125~0.85 mm 時,最大沖蝕速率先變大后變小。在0.125~0.6 mm 時,最大沖蝕速率隨著支撐劑粒徑的增大而增大;在 0.6~0.85 mm 時,最大沖蝕速率隨著支撐劑粒徑的增大而減小。其主要原因是,粒徑在0.125~0.6 mm 時,粒徑增加,顆粒具有的動能逐漸增大,顆粒撞擊壁面的磨損程度也逐漸增大。當支撐劑粒徑增大到 0.6~0.850 mm 時,流體對顆粒的攜帶能力減弱,顆粒具有的動能下降,顆粒與壁面碰撞接觸面積也增大,相同砂比條件下,顆粒直徑增大,顆粒數量降低,同一位置顆粒對壁面的碰撞頻率降低,最終導致顆粒直徑增大到一定程度后,對壁面造成的最大沖蝕速率逐漸降低。

圖23 最大沖蝕速率隨支撐劑粒徑的變化曲線Fig.23 Curve of maximum erosion rate versus roppant particle size
入口流速25 m/s、砂比35%、支撐劑粒徑0.85 mm不變,變徑比為0.4、0.6、0.8、1.0,得到變徑T 型流道沖蝕速率對比如圖24 所示,內部顆粒堆積狀態對比如圖25 所示,最大沖蝕速率隨變徑比的變化曲線如圖26 所示,最大沖蝕速率隨位置的變化規律如圖27 所示。

圖24 變徑比不同下變徑T 型流道沖蝕速率對比Fig.24 Comparison of erosion rates of T-shaped runners with different reducing diameter ratios

圖25 變徑比不同下變徑T 型流道內部顆粒堆積狀態對比Fig.25 Comparison of particle accumulation conditions in T-shaped runners with different reducing diameter ratios

圖26 最大沖蝕速率隨變徑比的變化曲線Fig.26 Curve of maximum erosion rate changing with ratediameter ratio

圖27 不同變徑比下最大沖蝕速率隨位置的變化規律Fig.27 The curve of the maximum erosion rate with different variable diameter ratios as a function of position
由圖24 可知,變徑比從0.4 到1.0 時,流道壁面最大沖蝕速率從6.82×10?4mm/s 減小到5.33×10?4mm/s。由圖25 可知,顆粒在主流道迎流面上形成堆積,顆粒堆積的寬度從2 mm 增大到23 mm。形成顆粒堆積的原因是,流體出主流道向兩分支流道出口流動時,在迎流面形成低流速區域,顆粒在此處堆積。變徑比越小,形成的低流速區域面積越小,顆粒撞擊壁面起到的緩沖作用越小。入口流速相同時,變徑比為0.4、0.6、0.8、1.0,支流道出口中心處流體速度分別為4.23、7.69、9.41、11.08 m/s,證明變徑比越小,流體壓力損失越大,流體與支管路壁面碰撞更劇烈,顆粒與壁面作用面積小,形成高沖蝕率。因此,最大沖蝕速率隨著變徑比的增大而減小。
由圖24 和圖27 可知,隨變徑比的增大,最大沖蝕速率位置向支流道出口處移動,形成沖蝕圓環。變徑比為0.4 時,最大沖蝕速率在主流道中心對應的迎流壁面上。變徑比為0.6、0.8、1.0 時,最大沖蝕位置分別移動至距離主流道中心10、20、30 mm 附近,沖蝕面積向出口處擴大。
1)采用CFD-DEM 耦合模型,對變徑T 型流道四因素四水平的正交試驗進行數值計算,得出對沖蝕程度4 個試驗指標影響最大的因素是入口流速,其次是支撐劑粒徑、變徑比和砂比。各因素交互作用分析得出,入口流速與其他3 個因素交互作用的影響最大,其他因素交互影響不顯著。
2)隨著入口流速的增大,顆粒在低流速區域的堆積寬度減小12 mm,最大沖蝕速率位置不變,在距主流道中心兩側30 mm 附近,最大沖蝕速率隨流速增大4.612×10?4mm/s。支撐劑粒徑從0.125 mm到0.6 mm,最大沖蝕速率增大;從0.6 mm 到0.85 mm,最大沖蝕速率降低。變徑比為0.4~1.0 時,低流速區域堆積顆粒的寬度增大21 mm,最大沖蝕速率位置向支流道出口移動24 mm。砂比為5%~35%時,最大沖蝕速率隨砂比增大2.52×10?4mm/s,沖蝕位置不變。入口流速對沖蝕速率的影響排在第一位,支撐劑粒徑對沖蝕速率的影響排第二位,變徑比對沖蝕速率的影響排第三位,砂比對沖蝕速率的影響排第四位,與主效應分析結論一致。