朱維耀，張啟濤，岳 明，張燎源

1) 北京科技大學土木與資源工程學院，北京 100083 2) 中國石化勝利油田分公司石油工程技術(shù)研究院，東營 257000

水力壓裂技術(shù)已經(jīng)廣泛應(yīng)用于非常規(guī)油氣資源的開發(fā)中[1-5]. 在壓裂過程中需要注入大量的固體支撐劑顆粒，目的在于防止人造裂縫在閉合壓力的作用下完全封閉，從而為地下油氣資源的采出提供高導(dǎo)流能力的通道，提升油氣田的開發(fā)效果[6-7]. 支撐劑的分布狀態(tài)將直接影響生產(chǎn)井的生產(chǎn)效果，針對直井單裂縫中支撐劑在注入過程后的運移與沉降過程，前人已經(jīng)進行了較為詳細的研究[8-14]. 然而在實際開發(fā)過程中，地下的裂縫主要以相互交錯的縫網(wǎng)形式存在[15-16]，針對裂縫網(wǎng)絡(luò)內(nèi)的支撐劑分布，一些學者進行了物理實驗的研究[17-19]，但是相關(guān)數(shù)值模擬及理論研究相對缺乏. 同時，目前支撐劑的研究仍然主要集中在支撐劑的運移過程上，針對支撐劑分布狀態(tài)對生產(chǎn)動態(tài)規(guī)律影響的模擬研究相對較少，相關(guān)問題具有較大的研究價值.

為了準確地模擬支撐劑在裂縫的運移分布規(guī)律以及鋪砂過程結(jié)束后裂縫在實際儲層中的生產(chǎn)表現(xiàn)，本文在之前已經(jīng)建立的裂縫鋪砂數(shù)學模型[12]、裂縫導(dǎo)流能力數(shù)學模型[20]及含啟動壓力梯度致密儲層開發(fā)數(shù)學模型[21]的基礎(chǔ)上，建立了致密油壓裂開發(fā)全耦合計算模型，研究支撐劑在裂縫網(wǎng)絡(luò)的分布狀態(tài)對致密儲層開采動態(tài)規(guī)律的影響，以期為現(xiàn)場開發(fā)提供指導(dǎo).

1 數(shù)學模型選擇

1.1 裂縫網(wǎng)絡(luò)攜砂液流動數(shù)學模型

攜砂液主要由固體支撐劑和壓裂液組成，因此裂縫中的攜砂液流動是一種典型的固液兩相流問題. 采用混合物模型對攜砂液流動進行模擬計算. 混合物模型是一種簡化的兩相模型，只求解一組動量方程，即混合物的連續(xù)性方程和固相輸運方程，具有較高的計算效率和收斂性，適合裂縫等復(fù)雜幾何的模擬計算[22-24]. 本次模擬不考慮溫度對兩相流動的影響.

模型動量方程可以表示為[12]：

其中，u是混合物的質(zhì)量平均速度，m·s-1；p是壓力，Pa；g是重力加速度，m·s-2；cs是顆粒質(zhì)量分數(shù)；uslip是固相和液相間的相對速度，m·s-1；ρ是混合物的密度，kg·m-3，表示為：

其中，ρf和ρs分別是純液體和固體的密度，kg·m-3；φs是固相的體積分數(shù). 根據(jù)Krieger表達式[25]，混合物黏度η表示為：

其中，ηf是純液體的動力黏度，Pa·s；φmax是最大的顆粒濃度.

混合物模型中的連續(xù)性方程為[12]：

固相的輸運方程為[12]：

其中，us是固相速度，m·s-1，它可以表示為[12]：

將式（6）代入式（5）中[12]：

基于以上方程，使用COMSOL多物理場模擬軟件對攜砂液流動進行模擬.

1.2 裂縫網(wǎng)絡(luò)幾何參數(shù)和致密儲油層開發(fā)數(shù)學模型

鋪砂過程結(jié)束后，裂縫寬度和裂縫滲透率都會因受壓發(fā)生改變. 朱維耀等[20]根據(jù)支撐劑的最密排列以及Carman-Kozeny公式，結(jié)合彈性力學的相關(guān)知識，建立了考慮鋪砂濃度、閉合壓力、支撐劑粒徑、支撐劑材料以及支撐劑嵌入的裂縫滲透率與縫寬（即導(dǎo)流能力）計算模型. 通過該模型，鋪砂計算出的支撐劑顆粒濃度可以轉(zhuǎn)換為此時的裂縫寬度和滲透率，進而可以用于之后的致密儲層黑油模型計算中.

受壓后的裂縫寬度計算公式如下[20]：

其中，wfh為受壓后裂縫的寬度，m；r為支撐劑直徑，m；C為鋪砂濃度，由上文計算出的支撐劑顆粒濃度換算得到，kg·m-2；mo為單個支撐劑質(zhì)量，kg；h1為支撐劑嵌入裂縫內(nèi)的距離，m，其計算公式如下[20]：

其中，F(xiàn)為一個支撐劑所受的壓力，N；pc為閉合壓力，Pa；μ為支撐劑的彈性模量，Pa；E為支撐劑的泊松比；μ1為壁面材料的彈性模量，Pa；E1為壁面材料的泊松比.

式（8）中的Bt為破碎率，與支撐劑材料及閉合壓力大小有關(guān)：

其中，pm為支撐劑開始發(fā)生破裂的臨界閉合壓力，Pa；a和b分別為實驗擬合系數(shù)，量綱為一.

裂縫滲透率計算公式如下[20]：

其中，no為單位面積支撐劑的個數(shù)；τ為迂曲度.

獲得縫網(wǎng)裂縫寬度和滲透率之后，需要置入儲層模型中進行模擬生產(chǎn)，從而評價壓裂改造效果，獲得生產(chǎn)動態(tài)規(guī)律. 選取含啟動壓力梯度致密油儲層進行模擬. 針對含啟動壓力梯度的致密油儲層開發(fā)模擬，朱維耀等[21]提出了體積力近似法，可以較為準確地對含啟動壓力梯度儲層流體流動進行計算模擬. 本文模型假設(shè)儲層中的流體只能先流入裂縫，之后才能通過裂縫流入井筒.

基于以上方程，使用COMSOL多物理場模擬軟件及Livelink for MATLAB研究裂縫網(wǎng)絡(luò)支撐劑非均勻分布對致密儲層開采動態(tài)規(guī)律的影響.

2 裂縫網(wǎng)絡(luò)支撐劑運移分布模擬

2.1 裂縫網(wǎng)絡(luò)幾何建模與參數(shù)選取

基于Warren-Root正交裂縫網(wǎng)絡(luò)模型[26]以及Mayerhofer等[27]的縫網(wǎng)模型，進行了裂縫網(wǎng)絡(luò)幾何建模（圖1）. 裂縫網(wǎng)絡(luò)由兩條人工主裂縫和兩條橫向次級裂縫組成，主裂縫與次級裂縫長度均為100 m，主裂縫初始縫寬為0.02 m，次級裂縫初始縫寬為0.005 m，井筒半徑為0.1 m. 攜砂液由井筒首先流入主裂縫后再流入次級裂縫.

圖1 裂縫網(wǎng)絡(luò)幾何模型示意圖Fig.1 Illustration of fracture network model

為了簡化幾何模型，同時減少模型計算量，在實際計算時對裂縫網(wǎng)絡(luò)模型進行了對稱簡化處理（圖2）. 編號1、2和3為攜砂液出口，設(shè)置為壓力出口邊界且無固體通量；4、5和6a為對稱邊界；6b為入口，設(shè)置為速度入口邊界.

圖2 裂縫網(wǎng)絡(luò)對稱簡化示意圖Fig.2 Symmetric fracture network model

對上述幾何模型進行網(wǎng)格劃分，使用有限元法進行求解. 為了更精確地捕捉固相參數(shù)，在靠近入口和下邊界處進行網(wǎng)格加密處理. 計算中使用的基本參數(shù)如表1中所示.

表1 基本計算參數(shù)Table 1 Basic calculation parameters

2.2 裂縫網(wǎng)絡(luò)鋪砂計算結(jié)果分析

（1）裂縫網(wǎng)絡(luò)支撐劑運移與分布狀態(tài)規(guī)律.

根據(jù)上述模型與計算參數(shù)，首先對裂縫網(wǎng)絡(luò)內(nèi)的攜砂液流動進行模擬，計算時間為1000 s，計算結(jié)果如圖3所示.

圖3 裂縫網(wǎng)絡(luò)鋪砂過程示意圖Fig.3 Illustration of proppant-laden fluid injection process

從圖3中可以看出，裂縫網(wǎng)絡(luò)中支撐劑濃度的分布狀態(tài)要明顯不同于單裂縫. 支撐劑在進入裂縫網(wǎng)絡(luò)之后逐漸以“砂包”狀緩慢向前推進，在到達主裂縫與分支縫交匯點后，逐漸形成支撐劑堆積，裂縫交匯點處的支撐劑沉降堆積高度要明顯高于其他部分. 繼續(xù)注入攜砂液，支撐劑開始進入次級裂縫，從圖3（c）中可以看出，鋪砂過程結(jié)束后次級裂縫外側(cè)部分的砂堤堆積高度要比主裂縫中小25%～50%，其支撐劑懸浮區(qū)則明顯更大. 相比較而言，次級裂縫內(nèi)側(cè)部分由于溝通了兩條主裂縫，其支撐劑堆積高度也要明顯高于外側(cè)部分，改造效果更加明顯.

（2）鋪砂參數(shù)對裂縫網(wǎng)絡(luò)平均導(dǎo)流能力影響.

為了研究支撐劑物性參數(shù)以及鋪砂施工參數(shù)對鋪砂效果的影響，選取不同支撐劑粒徑、材料和攜砂液注入速度進行鋪砂模擬. 通過模擬得到縫網(wǎng)各處的支撐劑濃度后，通過式（8）和式（11）計算縫網(wǎng)內(nèi)各部分導(dǎo)流能力，再對縫網(wǎng)進行整體積分，取得裂縫網(wǎng)絡(luò)的平均導(dǎo)流能力，并將該導(dǎo)流能力參數(shù)作為衡量鋪砂效果進行參數(shù)優(yōu)選的定量指標.

為了研究支撐劑粒徑對裂縫網(wǎng)絡(luò)整體導(dǎo)流能力的影響，選取直徑為200、300和400 μm的支撐劑顆粒進行計算，在其他因素相同的情況下，研究不同閉合壓力下，不同粒徑支撐劑所形成裂縫平均導(dǎo)流能力. 計算結(jié)果如圖4所示，縫網(wǎng)平均導(dǎo)流能力與支撐劑粒徑大小呈正相關(guān)關(guān)系：在閉合壓力為10 MPa的條件下，粒徑400 μm計算出的導(dǎo)流能力要比200 μm情況大10 D·cm左右；當閉合壓力壓力增大到60 MPa時，二者則相差3 D·cm. 整體而言，大粒徑支撐劑的對裂縫導(dǎo)流能力的正作用較為明顯.

圖4 支撐劑粒徑對縫網(wǎng)整體導(dǎo)流能力的影響Fig.4 Effect of proppant diameter on fracture network conductivity

為了研究支撐劑材料對裂縫網(wǎng)絡(luò)整體導(dǎo)流能力的影響，選取石英和陶粒作為模擬對象. 陶粒彈性模量為100 GPa，泊松比為0.25；石英彈性模量為30 GPa，泊松比為0.15. 計算結(jié)果如圖5所示，在相同鋪砂條件下，石英支撐劑計算出的縫網(wǎng)平均導(dǎo)流能力要小于陶粒. 二者之間的差異大小與閉合壓力呈正相關(guān)關(guān)系. 當閉合壓力為0時，二者的導(dǎo)流能力基本相等，而當閉合壓力上升到60 MPa時，石英的破碎率遠高于陶粒，其導(dǎo)流能力要比陶粒小50%. 因此，在實際選取支撐劑種類時，支撐劑的材料也是必須要考慮的因素. 變形能力小，材料強度大，不易破碎的材料對于水力壓裂改造效果具有較大的正影響.

圖5 支撐劑材料對縫網(wǎng)整體導(dǎo)流能力的影響Fig.5 Effect of proppant materials on fracture network conductivity

為了研究攜砂液注入速度對裂縫網(wǎng)絡(luò)整體導(dǎo)流能力的影響，選取攜砂液的法向入口流速分別為0.3，0.4，0.5和0.6 m·s-1，計算裂縫網(wǎng)絡(luò)平均導(dǎo)流能力，計算結(jié)果如圖6所示.

圖6 攜砂液注入速度對縫網(wǎng)整體導(dǎo)流能力的影響Fig.6 Effect of proppant-laden fluid injection velocity on fracture network conductivity

從圖6可以發(fā)現(xiàn)，注入攜砂液速度越大，其裂縫網(wǎng)絡(luò)的整體導(dǎo)流能力越大，二者呈現(xiàn)正相關(guān)關(guān)系. 但是隨著流入速度逐漸增加，裂縫網(wǎng)絡(luò)整體導(dǎo)流能力的增長逐漸變緩，當注入速度大于0.6 m·s-1，時，再繼續(xù)增加注入速度，導(dǎo)流能力增長比較有限. 通過計算注入砂量大小以及裂縫平均導(dǎo)流能力的比值，最佳的攜砂液注入速度為0.4 m·s-1左右.

3 裂縫網(wǎng)絡(luò)支撐劑非均勻分布對致密儲層開采動態(tài)規(guī)律影響

傳統(tǒng)油藏數(shù)值模擬模型一般假定壓裂裂縫內(nèi)部滲透率相等，縫寬相等. 而實際情況下，由于支撐劑在裂縫中的運移分布差異，裂縫各處的滲透率與縫寬存在差異，這種差異往往會導(dǎo)致計算誤差，影響模擬計算及結(jié)果的準確性. 通過對比支撐劑理想分布條件和支撐劑非均勻分布條件下的致密儲層生產(chǎn)表現(xiàn)，定量研究支撐劑非均勻分布對油藏數(shù)值模擬計算結(jié)果的影響，模型計算參數(shù)如表2所示.

表 2 基本計算參數(shù)Table 2 Basic calculation parameters

圖7 裂縫網(wǎng)絡(luò)內(nèi)支撐劑理想均勻分布與不均勻分布條件下的致密儲層開發(fā)300 d壓力場對比. （a）支撐劑均勻分布狀態(tài)；（b）基于支撐劑均勻分布的儲層壓力場分布；（c）支撐劑不均勻分布；（d）基于支撐劑不均勻分布的儲層壓力場分布Fig.7 Comparison between pressure distribution based on proppant idealized and uneven distribution at 300 days in a tight oil reservoir: (a) idealized proppant distribution; (b) pressure distribution with even proppant distribution; (c) uneven proppant distribution; (d) pressure distribution with uneven proppant distribution

支撐劑理想均勻分布如圖7（a）所示，支撐劑非均勻分布如圖7（c）所示. 兩種分布狀態(tài)衍生計算出的致密油儲層生產(chǎn)300 d后壓力場分布分別如圖7（b）和圖7（d）所示. 可以看出考慮支撐劑非均勻分布的條件下，裂縫網(wǎng)絡(luò)的有效控制面積相對較小，壓降主要集中在主裂縫周邊；而支撐劑理想分布的條件下，裂縫網(wǎng)絡(luò)的有效控制面積更大，壓力波及范圍更廣，壓降集中在主裂縫以及次級縫周邊. 因此在油藏數(shù)值模擬的過程中，將所有裂縫都假定為理想等寬裂縫的做法并不合理，需要針對儲層特性以及壓裂施工方式進行具體分析.

研究同時發(fā)現(xiàn)，基質(zhì)滲透率大小會對考慮與不考慮支撐劑不均勻分布兩種情況之間的產(chǎn)量差異產(chǎn)生影響. 在其他參數(shù)相同的條件下，改變儲層的基質(zhì)滲透率，計算兩種情況下的產(chǎn)量增長系數(shù)[27]進行對比，計算結(jié)果如圖8所示.

圖8 基質(zhì)滲透率對計算產(chǎn)量差異的影響Fig.8 Effect of matrix permeability on difference in calculated productivity

從圖8中可以看出，兩種情況下的產(chǎn)量差異與基質(zhì)滲透率呈反比關(guān)系，該結(jié)論與Kong等[28]的計算結(jié)果基本一致. 當模型基質(zhì)滲透率為0.05 mD時，支撐劑理想均勻分布條件下的產(chǎn)量要比考慮支撐劑不均勻分布的情況高41.7%，在這種情況下不能忽略支撐劑非均勻分布的影響. 而當基質(zhì)滲透率為5 mD時，二者的產(chǎn)量差異在5%以內(nèi)，此時如果用理想化等寬等滲裂縫進行油藏數(shù)值模擬則較為合理.

4 實例應(yīng)用

以勝利油田某低滲透區(qū)塊一水平井為例，完鉆井深為3542 m，油層厚度為20 m. 根據(jù)壓裂層測井解釋情況，平均滲透率為5.8 mD，地面原油平均密度為0.858 g·cm-3，平均原油黏度為16.68 mPa·s，水平井壓裂14段，平均段間距為67.2 m，平均裂縫半長為152.1 m. 依照現(xiàn)場提供模擬參數(shù)，首先對目標儲層進行鋪砂模擬，而后轉(zhuǎn)為生產(chǎn)模擬，計算產(chǎn)量與實際油田產(chǎn)量如圖9所示.

從對比圖中可以看出，鋪砂-生產(chǎn)耦合模擬計算得出的產(chǎn)量曲線與實際水平井產(chǎn)量曲線較為貼合，驗證了模型的準確性，同時為模型今后進一步應(yīng)用與實際生產(chǎn)奠定了基礎(chǔ). 該模型具有較大的應(yīng)用與拓展空間.

圖9 耦合計算模型結(jié)果與實際生產(chǎn)數(shù)據(jù)對比Fig.9 Comparison of oil rate between coupling model and oil field data

5 結(jié)論

為了研究裂縫網(wǎng)絡(luò)內(nèi)支撐劑運移分布規(guī)律以及支撐劑不均勻分布對致密儲層的影響，基于作者提出的多個數(shù)學模型，構(gòu)建了致密儲層水力壓裂生產(chǎn)耦合計算模型，對致密油儲層壓裂開發(fā)生產(chǎn)全過程進行模擬，得到了以下幾點結(jié)論：

（1）水平井裂縫網(wǎng)絡(luò)中的支撐劑運移與分布狀態(tài)與垂直井單裂縫存在明顯差異. 在裂縫網(wǎng)絡(luò)中，支撐劑會在裂縫交匯處產(chǎn)生堆積，堆積高度高于縫網(wǎng)其他部分. 次級裂縫中的支撐劑更多的處于懸浮狀態(tài)，且支撐劑堆積高度相較于主裂縫小25%～50%. 溝通次級縫具有更高的支撐劑沉降堆積程度，有利于油氣資源的開發(fā).

（2）裂縫網(wǎng)絡(luò)平均導(dǎo)流能力與支撐劑粒徑大小呈正比. 縫網(wǎng)導(dǎo)流能力還與法向注入速度呈正比，考慮注砂經(jīng)濟性，該模型條件下最佳的攜砂液注入速度為0.4 mμs-1左右.

（3）在進行油藏數(shù)值模擬時應(yīng)考慮支撐劑非均勻分布對計算結(jié)果的潛在影響. 計算結(jié)果顯示，當基質(zhì)滲透率為0.05 mD時，支撐劑理想均勻分布條件下計算出的產(chǎn)量要比考慮支撐劑非均勻分布的情況高41.7%；當基質(zhì)滲透率逐漸增大到5 mD時，產(chǎn)量差異在5%以內(nèi)，此時忽略支撐劑非均勻分布相對合理.