基于CFD-DEM算法的壓裂支撐劑縫內鋪砂規律*

郭慧娟 高暢 韓金井 范應璞 Guodong Zhan 王宏偉 張金亞

(1.中國石油集團工程技術研究院有限公司 2. 中國石油大學(北京)機械與儲運工程學院 3. 中國石油集團川慶鉆探公司井下作業公司 4. Saudi Aramco Drilling Technology Division EXPEC Advanced Research Center)

0 引 言

近年來，隨著頁巖氣、致密氣的大規模開發，水力壓裂技術成為油氣資源增產穩產的關鍵技術[1-3]。該技術在儲層壓裂后井下壓力小于閉合壓力時，支撐劑保持裂縫張開形成高滲流通道，進而實現生產，是目前低滲透油氣藏、頁巖氣藏實現開發的有效手段和利器。但由于當前現場施工條件以及成本限制，采用的傳統支撐劑與壓裂液使得支撐劑顆粒在輸送至目標距離前就大量沉積在近井地帶[4]，裂縫有效支撐距離短，開發效果不理想。因此，研究支撐劑在壓裂縫內鋪砂規律的影響具有重要意義。

目前，很多學者對支撐劑在裂縫內的沉降運移規律進行研究[5-8]。L.R.KERN等[9]對單一流速下支撐劑在水中流動首次進行了試驗，發現了攜砂液速度、支撐劑濃度以及支撐劑粒徑大小具有函數關系。R.E.BABCOCK等[10]在L.R.KERN研究的基礎上，針對牛頓流體與非牛頓流體流體進行物理試驗，結果表明支撐劑濃度是裂縫高度的連續函數。張濤等[11]采用CFD數值模擬的方法，建立雙歐拉兩相流模型，研究了清水壓裂中不同進口速度及位置，砂礫密度的支撐劑沉降行為，考慮了流動的湍流效應、高質量分數下顆粒間的摩擦應力。任嵐等[12]采用有限體積法求解，分析了支撐劑密度為1 300與2 600 kg/m3在復雜裂縫中沉降運移行為，并分析了在分支裂縫阻力效應下多種因素對于支撐劑運移的影響。梁瑩等[13]采用線性擬合的方法對低密度支撐劑在3種不同黏度壓裂液下，不同排量和不同砂比砂堤的鋪置形態進行擬合，應用于現場試驗并具有良好的指導作用。ZHANG G.D.等[14-15]采用耦合EDEM與Fluent的方法，對比了多尺寸顆粒相對于均尺寸顆粒對支撐劑放置的影響，得出不同粒徑的支撐劑形成的不同沙丘形態流場產生影響。WANG X.Y.等[16]研究裂縫閉合、注入速度和入口不均勻打開對單個裂縫中支撐劑分布的影響，以及支撐劑直徑、密度和射孔高度差對支撐劑分選的影響。ZENG J.S.等[17-18]用CFD-DEM耦合研究支撐劑顆粒在交叉裂縫中的運移機理，考慮不同裂縫交角以及顆粒聚團效應，發現支撐劑進入支縫的比例與裂縫交角呈負角相關，在狹縫中會出現支撐劑非均勻聚團現象。支撐劑運移的本質問題是流體與固態顆粒的復雜耦合傳動問題。基于計算流體力學建立的支撐劑運移模型大都采用歐拉框架，將流體-顆粒系統假設為互相摻混的2個流體組成的系統，要求固液兩相均為連續相[19-20]，適用于顆粒質量分數低于10%的情況，且忽略了顆粒體積對流場的影響，對支撐劑在壓裂縫中動力學表征較為有限。

為此，筆者采用計算流體力學-離散元(CFD-DEM)耦合方法對矩形壓裂縫中顆粒的運移進行數值計算，采用單一因素控制變量法探究支撐劑密度、攜砂液流速、質量分數以及壓裂液黏度等因素對壓裂縫支撐劑運移的影響規律，以期為選擇合適的攜砂液及泵注排量提供理論依據。

1 模型與網格

本文將輸送模型簡化為矩形，建立中等矩形裂縫的三維模型，采用計算流體力學-離散元(CFD-DEM)耦合方法對支撐劑的運移展布進行數值模擬。建立的三維幾何模型如圖1所示。圖1中幾何模型為長3 000 mm、寬6 mm、高600 mm的單一矩形平板裂縫。數值模型參考了可視化沉降試驗的平板裂縫，左側為3個均勻分布的入口，大小為60 mm×6 mm (頂部入口中心距離上壁面50 mm，底部入口中心距離下壁面50 mm，中部入口位于壁面正中間)，右側出口大小為40 mm×6 mm，模型采用Gambit進行結構性網格劃分，長度方向上網格數目為600個，高度方向上網格數為120個，寬度方向上網格數為2個，網格數量總計144 000個。假設流體為不可壓縮牛頓流體。顆粒為規則的圓形顆粒，球體不發生形變且與流體在入口處充分均勻混合。

圖1 三維幾何模型示意圖Fig.1 Schematic diagram of the 3D geometric model

2 數值模型與求解方法

在支撐劑運移的數值計算中，不僅要考慮顆粒運動對流體的影響，也要考慮流體對顆粒的反作用力。壓裂縫中流體流動控制方程為質量守恒及動量守恒方程，采用k-ε湍流模型使連續相計算方程組封閉，顆粒相采用EDEM中的單顆粒平動方程、轉動方程及固液相間作用力模型。

2.1 連續相控制方程

2.1.1 質量守恒方程

流體相的運動由Navier-Stokes方程控制，方程如下：

(1)

(2)

式中：ρ為流體的密度，kg/m3；εl為流體體積分數；ul為流動速度，m/s；k為流體相湍動能，m2/s2；ε為湍流耗散率，m2/s3；t為時間，s；Gk為湍動能產生相，kg/(m·s3)；μt為湍流黏度，Pa·s；μ為液相動力黏度，Pa·s；σk為湍動輪k對應的Prandtl數，σk=1.0；ui為X方向速度分量。

2.1.2 動量方程

-α1?p1+?τ1+α1ρ1g+β(vs-v1)

(3)

τl=μ(?u+?uT)

(4)

式中：p1為流體壓力，Pa；g為重力加速度，g=9.81 m/s2；τ1為剪切應力張量，Pa；β為相間動量交換系數，kg/(m3·s)；α1為動量修正系數；vs為出口流速，m/s；v1為入口流速，m/s；u為時均速度，m/s。

2.1.3 湍流方程

本文采用Realizablek-ε模型，在應變率大的情況下，會產生正應力小于0的情況，該模型對正應力進行了數學修正。

Gk+Gb-ρε-YM

(5)

湍流擴散方程：

(6)

式中：σε為湍動能ε對應的Prandtl數，σε=1.2；Gk、Gb、YM參數與標準化k-ε模型相同，Pa/s；C1ε、C3ε、C1、C2、A0為經驗常數；ak為修正系數；Meff為液體動力黏度，Pa·s；E為時均應變率，s-1；v為液體流速，m/s。

2.2 顆粒相控制方程

2.2.1 單顆粒受力方程

Fk=FD+FVM+FL+Fg+FF

(7)

式中：FF為浮力，N；Fg為重力，N；FL為升力，N；FVM為虛擬質量力，N；FD為拖曳力，N。

2.2.2 單顆粒平動方程和轉動方程

EDEM通過跟蹤單個顆粒軌跡來預測固相行為，顆粒的運動遵循牛頓第二定律。

平動方程：

(8)

轉動方程：

(9)

式中：mp，i為顆粒i的質量，kg；up，i為顆粒i的線速度，m/s；Fpc，ij為顆粒i與其他顆粒接觸產生的接觸力，N；Flp，i為流體對顆粒i的作用力，N；Ipc，ij為顆粒i的轉動慣量，kg·m2；ωp，i為顆粒i的角速度，rad/s；Tpc，ij為顆粒i與顆粒j接觸產生的接觸力矩，N·m。

(10)

2.3 兩相相間耦合模型

在固液兩相流動中，顆粒受到液體的多種力，這里主要考慮浮力和拖曳力。

2.3.1 固液相間作用力模型

浮力為：

(11)

拖曳力為：

(12)

(13)

虛擬質量力：

(14)

升力為：

FL=CLαsρl(ur)×(?×ul)

(15)

2.3.2 液相對固相的曳力模型

液相對固相的曳力計算采用Gidaspow模型，該模型是基于Ergum模型和Wen&Yu模型得到的經驗模型，適用于稠密的流化床模擬。其中εl代表液體體積分數，當液體體積分數小于0.8時，顆粒雷諾數對曳力無影響；當液體體積分數大于0.8時，顆粒雷諾數對曳力產生影響。

(16)

(17)

(18)

式中：Rep為顆粒雷諾數，無因次；εs為固相體積分數。

2.4 邊界條件與求解方法

上述控制方程采用有限體積法求解，湍流模型采用Realizablek-ε模型，強化壁面湍流，邊界設置為速度入口，壓力出口0.1 MPa，支撐劑泊松比0.5，靜摩擦因數0.5，動摩擦因數0.08，顆粒直徑1 mm，無滑移壁面條件，采用SIMPLE算法求解速度與壓力的耦合，設置DEM的時間步長2×10-5s，模擬時長20 s。

3 模擬方法驗證

為了驗證流固耦合模型的準確性，對試驗流速為0.64 m/s時支撐劑運移的模擬結果與韓琦[21]的平板裂縫沉降試驗結果進行對比(見圖2)。試驗裝置為大型可視化單一平板裂縫模型，對比參數如表1所示。

圖2 試驗流速0.64 m/s試驗模擬對比圖Fig.2 Experimental and simulation results at the flow rate of 0.64 m/s

表1 試驗模擬對照參數表Table 1 Experiment and simulation parameters

由圖2可以看出，砂堤整體呈現出“兩段狀”(由于計算資源受限，模擬計算時間為30 s)，兩段砂堤的中間部分沉降床高度與計算時間呈正相關，逐漸呈現出從“2個砂堤→1個砂堤”的趨勢。試驗主砂堤的位置約為2.15 m處，模擬時主砂堤的位置約為1.75 m處。分析砂丘位置產生差異的主要原因是仿真時采用顆粒直徑為1 mm進行計算，而試驗時采用1～2 mm混合直徑支撐劑，造成初始堆積位置的差異。可以看出支撐劑的運移沉降行為試驗與模擬情況較為相符，證明所建立流固耦合模型的準確性。

4 模擬結果分析

4.1 支撐劑密度對運移規律的影響

在攜砂液黏度為1 mPa·s，射流速度為1.6 m/s，顆粒質量分數4%工況下，分別對3種密度的支撐劑進行模擬，結果如圖3所示。

本文用攜砂液傳輸末端質量分數梯度角α來評價支撐劑的運移情況，傳輸末端質量分數梯度角α表示傳輸末端攜砂液質量分數沿流動方向的變化速率。由于入口的射流作用，攜砂液在入口處呈現“彈狀流”狀態，流場擾動較大。隨時間推移攜砂液裂縫中流動，流場逐漸平緩，此時支撐劑的受力情況較為穩定，傳輸末端質量分數梯度角α越大，表明沉降床鋪置更均勻且至裂縫深處，輸砂行為整體更接近于活塞流。由圖3a可知：高密度支撐劑(2 850 kg/m3)由于重力的對流作用在入口處迅速發生沉降，砂堤整體長度為1.50 m，傳輸末端質量分數梯度角α急劇變小；中密度支撐劑(1 950 kg/m3)在距離入口0.15 m處發生沉降，砂堤整體長度為2.15 m，傳輸末端質量分數梯度角α接近150°；低密度支撐劑(1 190 kg/ m3)在距離入口0.30 m處發生沉降，運移至裂縫最末端，在裂縫中鋪置更均勻，傳輸能力更強。隨著注入時間的延長，沉降床的高度穩定提升，傳輸末端質量分數梯度角α為平角。

圖4為不同密度支撐劑的液相體積分數分布圖(圖4a、圖4b、圖4c體積分數分別是24.6%、27.3%、32.8%)。較之高體積分數支撐劑，中低密度支撐劑的流化區域更大，且在沉降床表層流化度更高，表明支撐劑的密度越小，支撐劑堆積越松散，運移性能越好。

圖3 不同密度支撐劑運移圖Fig.3 Migration of proppants with different densities

圖4 不同密度支撐劑液相體積分數分布Fig.4 Liquid phase volume fraction distribution of proppants with different densities

圖5為3種密度支撐劑不同時刻的運移過程圖。高密度支撐劑(見圖5a)在入口處迅速沉降，隨著混砂液向前推進，液體的動能逐漸耗盡，支撐劑曳力小于重力，在注入時間為5 s時，在入口附近支撐劑沉降和堆積為3個小沙丘。隨著注入時間延長至10 s時，整體表現為沉降床高度逐漸增高，砂堤形態并無明顯變化。注入時間15 s時，砂堤高度增長速度變快，在越過第一個沙丘后，流過沉降床上方的攜砂液流速變大，注入的顆粒逐漸向沙丘背側堆積，因此沙丘形態逐步由3個沙丘變為2個沙丘，之后延長注入時間只沿沉降床高度方向變化。中密度支撐劑(見圖5b)在距離入口0.15 m處發生沉降，隨著攜砂液向裂縫深處運移，沉降床逐漸由2個沙丘合并為1個“駝峰式”沙丘。低密度支撐劑(見圖5c)由于湍流作用，支撐劑在出口被“吸出”，支撐劑沉降床無明顯沙丘。

圖5 3種密度支撐劑的不同時刻的運移過程圖Fig.5 Migration process of proppants with three densities at different time

4.2 攜砂液黏度對運移規律的影響

在支撐劑密度為1 190 kg/m3，射流速度為1.6 m/s，顆粒質量分數4%工況下，分別對攜砂液黏度為1、2、3 mPa·s進行模擬。圖6為不同黏度壓裂液的運移規律對比圖。由圖6可知，由于液體黏度升高，顆粒的隨附性變強，大量支撐劑隨壓裂液在出口流出，沉降顆粒變少，所以高黏度壓裂液較之中低黏度壓裂液的沉降床高度變低，同時液體的湍流強度變弱，具體體現為懸砂區的“波紋”情況逐漸減緩。說明支撐劑的運移性能隨著壓裂液的黏度增加而逐漸變好。在裂縫最末端支撐劑沿壁面向出口延伸，砂丘高度逐漸變高。當液體黏度為1 mPa·s時，最小體積分數為24.6%，而液體黏度為2和3 mPa·s時最小體積分數均為19%。這說明當液體黏度增加至一定程度時，支撐劑的懸浮性能達到極限，繼續增加壓裂液的黏度并不會對支撐劑的懸浮性能產生積極的影響，反而影響支撐劑的返排，對地層產生傷害。

圖6 不同黏度壓裂液支撐劑運移圖Fig.6 Proppant migration under different fracturing fluid viscosities

4.3 支撐劑流速對運移規律的影響

在攜砂液黏度為1 mPa·s，支撐劑密度為1 190 kg/m3，顆粒質量分數4%工況下，分別對現場施工排量為7.0、4.0、2.5 m3/min進行模擬。對應的模擬速度如表2所示

表2 現場施工排量與模擬速度對比表Table 2 Comparison ofon-site pump rate and simulated flow rate

圖7為不同流速支撐劑的運移規律圖。由圖7可知：速度為1.60 m/s時，支撐劑在距離入口0.3 m時就會發生沉降，顆粒主要堆積在中部靠近入口一側；速度為2.67 m/s時，支撐劑在距離入口0.9 m處發生堆積，顆粒主要堆積在中部靠近出口側，由于尾部沉降床高度增高，支撐劑的出口速度逐漸變大，輸砂區的面積較大且整體湍流強度更強；速度為4.27 m/s時，支撐劑沉降堆積在出口附近，運移狀態發生明顯改變，由于動能較大，輸砂區表現為整體橫向運移，攜砂液流動接近于活塞流，在出口被“吸出”。

圖7 不同流速支撐劑的運移規律圖Fig.7 Proppant migration at different flow rates

4.4 支撐劑顆粒質量分數對運移規律的影響

在支撐劑密度為1 190 kg/m3，攜砂液黏度為1 mPa·s，射流速度為1.6 m/s工況下，分別對顆粒質量分數為4%、6%、8%、10%的工況進行模擬，結果如圖8所示。

由圖8可以看出，顆粒質量分數與沉降床高度呈正相關性，相比圖8a與圖8b的沉降床高度增長，圖8c與圖8d的增長幅度減緩。隨著砂比的繼續增加，支撐劑擾動效應使得支撐劑增加速率下降，支撐劑沉降床的高度逐漸接近支撐劑的“平衡高度”，攜砂液的流速也逐漸增大至“平衡流速”，相當于增加了輸砂液的流速，使得沉降的支撐劑變少，懸浮態支撐劑明顯增多，支撐劑高度增加逐漸減緩。當沉降床砂堤高度達到“平衡高度”時，支撐劑的重力與浮力達到平衡，不會繼續沉降。顆粒質量分數6%、8%、10%對應流化區的液相體積分數分別為21.8%、24.6%、30.2%，證明支撐劑的沉降速度逐漸降低，高濃度支撐劑由于更為接近“平衡高度”。沉降床表面流化區域更大，一部分支撐劑由于堆積滑落至入口處產生堆積，長時間注入將堵塞入口，發生“卡泵”現象。

圖8 不同顆粒質量分數支撐劑運移圖Fig.8 Proppant migration under different particle concentrations

5 結 論

(1)同工況下，低密度支撐劑輸送距離較中密度支撐劑提升30.3%，較高密度支撐劑提高86.7%，密度較低的支撐劑輸送性能更好。在現場作業時，可以采取梯度密度支撐劑的注入方式，可以優先注入低密度支撐劑運送至裂縫深處，后期分別注入中高密度支撐劑對整個裂縫進行填充。

(2)一定程度上增加壓裂液黏度，會使攜砂液的湍流程度變弱，壓裂液的攜砂性變強，在實際壓裂作業中，可采用高黏性壓裂液將支撐劑攜帶到裂縫深處。

(3)較高的流速可以使支撐劑沉降距離入口越來越遠，實際壓裂作業時，可以考慮脈沖變速注入的方式，既可以防止入口堵塞，也可以保證一部分支撐劑在較高流速攜砂液的帶動下進入裂縫末端。

(4)適當提高支撐劑質量分數可以提前達到支撐劑沉降床高度，但質量分數過高的支撐劑會在入口出現堆積現象，混輸泵工作時，容易引起“卡泵”現象。現場壓裂作業時，要保證合理的注入濃度，也可以采取變速加砂和梯度混砂的方式避免堵塞情況發生，同時提高支撐劑的輸送效率。