基于CFD-DEM方法的迂曲裂縫中支撐劑運移關(guān)鍵影響因素分析

王雪飛，王素玲，侯 峰，王 明，李雪梅，孫丹丹

(1.東北石油大學，黑龍江 大慶 163318；2.中國石油新疆油田分公司，新疆 克拉瑪依 834002；3.中科院高能物理研究所，北京 100000；4.散裂中子源科學中心，廣東 東莞 523000；5.齊齊哈爾大學，黑龍江 齊齊哈爾 161006；6.中國石油大慶油田有限責任公司，黑龍江 大慶 163514)

0 引 言

在壓裂中，壓力梯度驅(qū)動的流動以及支撐劑和周圍攜砂流體的密度差導致的對流，二者的共同作用是支撐劑在裂縫中運移的主要原因[1-8]。國內(nèi)外學者通過實驗和數(shù)值模擬方法，對支撐劑運移規(guī)律進行了大量研究。侯磊等[9]將復雜裂縫系統(tǒng)抽象成單元化物理模型，開展了分流量和轉(zhuǎn)向條件研究。喬繼彤等[10]考慮支撐劑沉降因素，建立了支撐劑二維輸送方程以及含砂液的二維運動方程，應用有限元方法模擬了水力壓裂中的支撐劑輸送過程。Kern等[11]基于室內(nèi)實驗方法，對砂漿混合物中砂粒運移機理進行了研究。Shokir等[12]采用小玻璃模型模擬水力裂縫，研究注入速率、支撐劑用量和聚合物類型等參數(shù)對支撐劑運移的影響。Chun等[13]研究了裂縫傾角、支撐劑粒徑和注入點位置對支撐劑運移的影響。張濤等[14]建立Euler-Euler流固耦合雙流體模型，模擬單縫內(nèi)水輸送支撐劑的行為。徐加祥等[15]利用COMSOL軟件中的自由與多孔介質(zhì)流動耦合模型，對壓裂液在微裂縫中的流動進行模擬。Zeng等[16]采用代表性顆粒模型(RPM)對支撐劑顆粒進行放大，研究支撐劑在裂縫中的運移機理。Siddhamshetty等[17-18]采用三維多相顆粒單元(MP-PIC)模型，利用歐拉-拉格朗日方法，模擬多尺寸支撐劑的運移。Rivas等[19]利用有限體積法建立支撐劑通過水力裂縫運移的數(shù)值模型，模擬KGD裂縫中支撐劑的運移。Kou等[20]使用CFD-DEM方法研究支撐劑在原生裂縫和次生裂縫中運移。Zeng等[21]基于浸沒邊界法和解析CFD-DEM方法模擬水-顆粒耦合機制的支撐劑運移過程。Wu等[22]利用CFD-DEM方法，模擬水平井射孔套管中支撐劑的運移過程。Wang等[23]采用CFD-DEM方法研究影響支撐劑在裂縫及裂縫網(wǎng)絡中分布規(guī)律的因素，研究裂縫閉合程度、注入速度和入口開度對單個裂縫中支撐劑分布的影響。

上述對支撐劑運移和鋪置的研究中均將裂縫簡化為平直裂縫，忽略了迂曲裂縫對支撐劑運移和鋪置的影響?，F(xiàn)場微地震監(jiān)測和實驗室模擬壓裂裂縫均顯示壓裂裂縫具有迂曲特點[24-27]。因此，研究支撐劑鋪置規(guī)律時考慮裂縫迂曲情況十分必要?；贑FD-DEM方法，研究了裂縫迂曲度、支撐劑注入速度、壓裂液砂比、支撐劑粒徑對迂曲裂縫中支撐劑運移和鋪置的影響，為現(xiàn)場壓裂作業(yè)提供理論依據(jù)。

1 顆粒與流體耦合計算方法

目前顆粒與流體的耦合計算方法主要有歐拉-歐拉法[28-29]和歐拉-拉格朗日法[30]。此次研究采用歐拉-拉格朗日法，流體采用Navier-Stokes方程求解，而顆粒采用離散元方法，利用牛頓運動定律預測顆粒運動軌跡，并判斷顆粒相互碰撞情況。歐拉-拉格朗日法結(jié)合了流體相的連續(xù)描述和基于牛頓運動方程的顆粒相拉格朗日表征，這種模擬兩相流的組合方法可提供兩相耦合微觀行為的詳細信息。

建立考慮顆粒-顆粒、顆粒-壁面和流體-顆粒相互作用的固液兩相流耦合模型。對模型做以下假設：流體是不可壓縮的，流體流變性和溫度在模擬過程中不發(fā)生變化，流體與壁面間不滑移，支撐劑是規(guī)則的球形顆粒。由于歐拉-拉格朗日模型計算耗時較長，因此，采用實驗室規(guī)模的裂縫模型。

1.1 流體控制方程

流體質(zhì)量守恒方程：

(1)

流體動量守恒方程：

(2)

兩相流與單向流的區(qū)別在于顆粒在計算單元內(nèi)占據(jù)的流體體積。改變流體體積分數(shù) 將影響流體運動狀態(tài)，流體運動狀態(tài)影響顆粒受到的流體作用力，從而改變流體與顆粒動量交換力源項 。

(3)

(4)

1.2 顆粒受力及動力學方程

顆粒-流體的流動系統(tǒng)中，顆粒進行的平移和旋轉(zhuǎn)運動符合牛頓第二定律[31]，動力學方程為：

(5)

(6)

顆粒受到的流體作用力的方程為：

Ff=Fd+Fvm+Fpg+Fsa+Fb+Fm

(7)

式中：Fd為顆粒受到的流體阻力(曳力)，N；Fvm為虛擬質(zhì)量力，N；Fpg為壓力梯度力，N；Fsa為Saffman升力，N；Fb為Basset力，N；Fm為Magnus力，N。

已有研究[32]認為，作用于顆粒上的力主要是流體阻力，其他力比較小可忽略不計。因此，顆粒受到流體作用力只考慮流體阻力。

當黏性流體繞流顆粒時，流體對顆粒產(chǎn)生與流動方向一致的作用力，該作用力與流體和顆粒相對速度(uf-up)有關(guān)。當流體的流速小于顆粒的運動速度時，作用于顆粒上的力表現(xiàn)為阻力；當流體的流速大于顆粒運動速度時，作用在顆粒上的力表現(xiàn)為曳力，計算公式為：

(8)

(9)

(10)

式中：Cd為曳力系數(shù)；Rep為顆粒雷諾數(shù)；Dp為顆粒直徑，m；μf為流體表觀黏度，Pa·s。

1.3 顆粒碰撞模型

如圖1所示，在外力作用下顆粒i與顆粒j在c點碰撞，虛線表示開始碰撞時顆粒j的位置，隨著兩顆粒相對運動，產(chǎn)生法向重疊量γ和切向位移δ(單位均為m)。通過引入彈簧、阻尼器、滑動器和耦合器等計算顆粒間產(chǎn)生碰撞力，如圖2所示。

圖1 兩顆粒碰撞示意圖Fig.1 The schematic diagram of collision between two particles

圖2 彈簧-阻尼系統(tǒng)接觸力模型Fig.2 The contact force model of spring-damping system

對于三維坐標中的球體顆粒，F(xiàn)σij為作用于顆粒i的碰撞力的法向分量，可等效為彈簧和阻尼器產(chǎn)生的力之和：

(11)

式中：uij為顆粒i與顆粒j的相對速度矢量，uij=ui-uj；a為顆粒i的中心至顆粒j的中心直線方向上的單位矢量；kσ、cσ分別為顆粒i的法向彈性系數(shù)和法向阻尼系數(shù)。

碰撞力的切向分量Fτij為：

Fτij=-kτδ-cτucτ

(12)

ucτ=uij-(uij·a)a+riωi×a+rjωj×a

(13)

式中：kτ、cτ分別為切向彈性系數(shù)和切向阻尼系數(shù)；ucτ為接觸點的滑移速度，m/s；ri、rj分別為顆粒i和顆粒j的半徑，mm；ωi、ωj分別為顆粒i和顆粒j的角速度，rad/s；

綜上所述，作用在顆粒i上的合力和合力矩分別為：

(14)

式中：Fij為作用在顆粒i上的合力，N；Tij為作用在顆粒i上的合力矩，N·m；。

當多個顆粒同時與顆粒i接觸時，作用在顆粒i上的總力和總力矩為：

(15)

式中：Fi為作用在顆粒i上的總力，N；Ti為作用在顆粒i上的總力矩，N·m；。

1.4 CFD-DEM耦合計算方法

CFD-DEM耦合求解流程：首先建立連續(xù)相與離散相數(shù)值模型，設置邊界條件，計算流體體積分數(shù)αf和動量交換力源項Sf，啟動CFD求解器，計算首次迭代的動量離散方程中的系數(shù)和常數(shù)項，求解動量離散方程，對速度和壓力進行修正，并判斷是否滿足收斂準則。若不滿足，則重新求解動量離散方程的系數(shù)和常數(shù)，再次修正速度和壓力，直至滿足收斂準則；若滿足收斂準則，流體域推進至下一個計算時間步。啟動DEM求解器，根據(jù)插值函數(shù)獲得顆粒中心位置流體速度uf和顆粒受到的流體作用力Ff，使用質(zhì)點坐標方法追蹤顆粒并更新顆粒位置和顆粒運動速度up，判斷顆粒碰撞情況并計算顆粒碰撞力，進行顆粒動力學方程求解，更新顆粒運動速度up，再次計算體積分數(shù)αf和動量交換力源項Sf，開始下一次流體求解，直至耦合求解結(jié)束。

2 邊界條件及驗證

2.1 邊界條件建立

在壓裂過程中，先將前置液注入地層，打開裂縫，然后注入均勻混合的壓裂液。為了準確模擬這一過程，先在CFD中計算穩(wěn)態(tài)流場直至收斂，將穩(wěn)態(tài)流場作為初始條件，建立裂縫物理模型，對流體域進行六面體網(wǎng)格劃分，并對邊界層進行網(wǎng)格加密處理。采用SIMPLEC方法和瞬態(tài)解算器求解流體控制方程。設置速度入口邊界和壓力出口邊界，顆粒與流體初始速度相同。

顆粒選用石英砂材質(zhì)，密度為2 650 kg/m3，彈性模量為5.5×1010Pa，泊松比為0.13。裂縫邊界為巖石，密度為3 000 kg/m3，彈性模量為6.5×1010Pa，泊松比為0.20。顆粒采用Hertz-Mindlin接觸模型，顆粒接觸參數(shù)采用虛擬實驗進行標定，恢復系數(shù)為0.40，靜摩擦系數(shù)為0.68，動摩擦系數(shù)為0.01。

2.2 數(shù)值模型驗證

基于建立的CFD-DEM耦合數(shù)值模型，對文獻[33]中“一”型裂縫砂丘鋪置形態(tài)實驗過程進行數(shù)值模擬，對比結(jié)果如圖3所示。由圖3可知：“一”型裂縫中砂丘鋪置高度和鋪置的距離、鋪置形態(tài)基本一致，驗證了建立的數(shù)值模型的準確性。

圖3 “一”型裂縫中支撐劑鋪置砂丘形態(tài)對比Fig.3 The comparison of proppant displacement shape in linear fracture

3 迂曲裂縫中支撐劑運移影響因素分析

裂縫迂曲度是裂縫中心線長度與裂縫首尾直線距離的比值。數(shù)值模擬中，流體和支撐劑從同一入口流入，初始速度相同。

3.1 迂曲度對支撐劑運移的影響

模擬條件：裂縫迂曲度為1.0、1.2、1.5、2.0，如圖4所示；流體注入速度為1.0 m/s，模擬時間為0.5 s。壓裂液在不同迂曲裂縫中的流速如圖5所示，流速損失率隨進入裂縫距離變化如圖6所示，支撐劑在裂縫內(nèi)分布如圖7所示。

圖4 裂縫迂曲度示意圖Fig.4 The schematic diagram of crack tortuosity

圖5 裂縫中流速分布Fig.5 The flow rate distribution in fracture

圖6 流速損失率Fig.6 The loss ratio of flow rate

圖7 迂曲度對裂縫中支撐劑分布的影響Fig.7 The effect of tortuosity on proppant distribution in fractures

由圖5可知：壓裂液沿裂縫流動時，流速隨著進入裂縫的距離增大而減小。其原因為壓裂液具有一定的黏度，在流動過程中與裂縫壁面間的摩阻力和裂縫形狀引起的局部阻力造成流體能量損失。迂曲度為1.0的裂縫內(nèi)流速呈線性下降，迂曲度為1.2、1.5、2.0的裂縫內(nèi)流速呈非線性下降。由圖6可知，隨著流體進入裂縫距離的增大，流體流速損失率增大。迂曲度越大流速損失越快，進入裂縫20 mm處，與迂曲度為1.0的裂縫中的流速損失率相比，迂曲度為1.2、1.5、2.0的裂縫的流速損失率分別增大2.2、3.7、4.5倍。其原因為流體進入裂縫距離相同時，迂曲度越大實際流過距離越長，流體在迂曲裂縫內(nèi)還會形成渦旋，流體產(chǎn)生的沿程和局部壓力損失增大。進入裂縫初期流速損失率增長最快，隨著進入裂縫距離的增大，流速損失率逐漸減緩。

由圖7a可知：迂曲度為1.0時，支撐劑向前平穩(wěn)運移，形成的砂丘坡度較緩，鋪置距離較長，支撐裂縫面積大。砂丘底層顆粒受上部顆粒的壓實作用和周圍顆粒的阻力停止運動，砂丘中部顆粒受周圍顆粒的阻力和壓裂液的攜帶，以較小速度向前滾動，大部分上層顆粒以較大速度和流體一起向前運移。

由圖7b、c、d可知：迂曲度由1.2增至2.0時，迂曲裂縫阻力較大，流體流動狀態(tài)復雜，流體流動速度不均勻；顆粒與壁面以及顆粒與顆粒之間碰撞加劇，導致流體和顆粒劇烈減速，顆粒在裂縫內(nèi)發(fā)生沉降，形成堵塞，造成大面積支撐劑沉降，只有少部分上層顆粒以較小速度向前滾動；砂丘鋪置距離較短，坡度較大，顆粒在裂縫入口堆積，無法有效支撐裂縫。

3.2 注入速度對支撐劑運移的影響

模擬條件：支撐劑注入速度為0.1、0.5、1.0、1.5 m/s，裂縫迂曲度為1.0、1.2、1.5、2.0，流體和支撐劑的初始速度相同，壓裂液砂比為10%。支撐劑在裂縫中分布如圖8所示(以迂曲度為1.2的裂縫為例，下同)；各迂曲度裂縫中注入速度對鋪置距離的影響如圖9所示。

圖8 注入速度對裂縫中支撐劑分布的影響Fig.8 The effect of injection rate on proppant distribution in fractures

圖9 注入速度對鋪置距離的影響Fig.9 The analysis of the effect of injection rate on displacement distance

由圖8a可知：注入速度為0.1 m/s時，砂丘堆積在裂縫入口處，繼續(xù)注入支撐劑砂丘高度增長迅速，最終堵塞裂縫入口，后續(xù)支撐劑將無法進入裂縫，而砂丘長度增長甚微，說明后續(xù)支撐劑起不到支撐裂縫作用。

由圖8b可知：注入速度為0.5 m/s時，迂曲裂縫內(nèi)的顆粒與顆粒、顆粒與壁面劇烈碰撞，流體局部壓力損失增大，大部分支撐劑在裂縫迂曲處沉降堆積成砂丘，繼續(xù)注入支撐劑，在裂縫迂曲狹窄處繼續(xù)堆積形成堵塞，后續(xù)注入支撐劑不能進入裂縫深處，導致裂縫不能得到有效支撐。

由圖8c可知：注入速度為1.0 m/s時，支撐劑在迂曲裂縫中減速明顯，但仍以一定速度邊向前運移邊沉降形成砂丘，相比注入速度低的情況，高注入速度形成的砂丘形狀較平緩，砂丘鋪置面積增大，繼續(xù)注入支撐劑能夠延長裂縫支撐長度。

由圖8d可知，注入速度為1.5 m/s時，支撐劑克服裂縫迂曲帶來的阻力，在裂縫中鋪置距離更長，更均勻。

由圖9可知：支撐劑的鋪置距離隨注入速度的增加而增大，但不同迂曲度的裂縫中支撐劑對注入速度的敏感度不同，注入速度從0.1 m/s增至1.5 m/s，迂曲度為1.0、1.2、1.5、2.0的裂縫中支撐劑鋪置距離分別延長1.85、1.72、1.52、1.30倍。

迂曲度為1.0的裂縫中，支撐劑運移時所受到阻力較小，隨著注入速度增加鋪置距離增長較快。隨著裂縫迂曲度的增大，支撐劑向前運移阻力越大，受力復雜，隨著注入速度增加鋪置距離增長變慢。當注入速度小于1.0 m/s時，支撐劑和流體的動能較小，支撐劑在裂縫內(nèi)運移受阻，發(fā)生快速沉降，鋪置距離較短。而注入速度大于1.0 m/s時，支撐劑和流體具有足夠的動能抵抗迂曲裂縫的形狀阻力，使得支撐劑越過迂曲地帶后仍然具有一定動能，運移一段距離才發(fā)生沉降。因此，迂曲度越大的裂縫，需要更大的注入速度才能獲得較理想的鋪置效果，注入速度較小時，支撐劑會阻塞裂縫入口，使后續(xù)支撐劑無法運移至裂縫深處，支撐裂縫效果較差。

3.3 壓裂液砂比對支撐劑運移的影響

模擬條件：壓裂液砂比為10%、20%、30%、40%，注入速度為1.0 m/s，支撐劑粒徑為0.6 mm。支撐劑在裂縫中分布如圖10所示。

由圖10可知：壓裂液砂比由10%增至40%時，支撐劑鋪置距離縮短。其原因為裂縫中注入高砂比壓裂液時，顆粒與顆粒、顆粒與壁面碰撞增多，顆粒動能損失增大，顆粒發(fā)生快速沉降，在裂縫入口附近發(fā)生堵塞，阻礙后續(xù)顆粒的進入，無法獲得長距離支撐效果。而低砂比壓裂液運移時，顆粒碰撞次數(shù)減少，顆粒能夠跟隨壓裂液進入裂縫深處，鋪置距離增大。由此可見，砂比越小支撐劑運移過程阻力損失越小，鋪置越均勻，但減小砂比會增加泵注時間。

圖10 壓裂液砂比對裂縫中支撐劑分布的影響Fig.10 The effect of fracturing fluid-proppant ratio on proppant distribution in fractures

3.4 支撐劑粒徑對支撐劑運移的影響

模擬條件：支撐劑粒徑為0.2、0.4、0.6、0.8 mm，注入速度為1.0 m/s，壓裂液砂比為10%。撐劑在裂縫中分布如圖11所示。

由圖11可知：支撐劑粒徑越大，壓裂液對支撐劑的攜帶性越弱，在迂曲裂縫內(nèi)通過性越差，與壁面發(fā)生多次碰撞后發(fā)生快速沉降，支撐劑鋪置距離短，且不均勻。支撐劑粒徑越小，壓裂液對其攜帶能力越強，支撐劑通過迂曲裂縫數(shù)量增多，能夠形成鋪置距離長且均勻的裂縫。

圖11 支撐劑粒徑比對裂縫中支撐劑分布的影響Fig.11 The effect of proppant particle size ratio on proppant distribution in fractures

4 結(jié) 論

(1) 壓裂液在裂縫中流動，入口流速相同條件下，隨著裂縫迂曲度的增大，壓裂液速度損失越大，進入裂縫初期流速衰減越快，隨著進入裂縫距離的增大，流速衰減變緩。裂縫迂曲度越小，支撐劑鋪置越均勻，鋪置距離越長，迂曲度過大會導致支撐劑在裂縫入口處堵塞，支撐劑運移距離縮短。

(2) 隨著支撐劑入口速度的增大，不同迂曲裂縫中支撐劑運移距離均有所增加，迂曲度較大的裂縫，支撐劑運移距離隨著速度的增加緩慢增加，迂曲度小的裂縫支撐劑運移距離隨著速度增加明顯增大。支撐劑注入速度大于1.0 m/s時，不同迂曲裂縫內(nèi)支撐劑均能鋪置一定距離，對裂縫形成有效支撐，速度越大支撐劑鋪置距離越長分布越均勻，對裂縫的支撐效果越好。

(3) 支撐劑在迂曲裂縫中運移距離隨壓裂液砂比的增加而減小。支撐劑粒徑越大，流體對其攜帶能力越差，支撐劑鋪置距離越短，且鋪置不均勻。