郭興，孫曉，穆景福，喬紅軍，羅攀，李恒

（1.陜西延長石油(集團)有限責任公司研究院,西安 710065；2.陜西省二氧化碳封存與提高采收率重點實驗室,西安 710065；3.北京大學工程學院,北京 100871）

0 引言

壓裂作為一種有效的儲層改造增產(chǎn)技術，在油氣藏勘探開發(fā)過程中應用廣泛，但其浪費大量的水資源、對環(huán)境及儲層造成較大的污染等弊端和局限性也逐漸暴露出來[1–2]。如何在當前“雙碳”背景下，“清潔、低碳、環(huán)?！钡倪M行油氣田勘探開發(fā)，保障我國能源安全，是目前亟需解決的問題。

超臨界CO2由于其高密度、低黏度、高擴散性等特殊性質(zhì)[3]，非常適合頁巖油氣等非常規(guī)油資源的勘探開發(fā)[4–5]。國內(nèi)外進行了諸多理論研究和現(xiàn)場試驗以及規(guī)模化CO2壓裂技術現(xiàn)場應用[6–7]，其可以降低破裂壓力，實現(xiàn)低滲致密儲層造復雜縫網(wǎng)，不僅節(jié)省大量水資源，對環(huán)境和儲層無污染，還能在油氣井增產(chǎn)提效的同時實現(xiàn)二氧化碳埋存[8–9]。超臨界CO2壓裂過程中支撐劑在裂縫中的有效運移和沉降是壓裂施工中極其重要的環(huán)節(jié)之一。

目前對壓裂過程中支撐劑運移規(guī)律研究主要采用室內(nèi)實驗和數(shù)值模擬2種方式[10]。通過實驗研究可以很直觀地看到支撐架在運移過程中的運動規(guī)律，分析得到小尺度（相比于儲層裂縫尺度）下描述顆粒群在裂縫中的沉降規(guī)律[11]。而對于超臨界CO2壓裂，在實驗室很難通過物理手段模擬其壓裂過程中支撐劑在CO2流體流動過程中的運移情況，具有極大的限制性，難以真正反映支撐劑在真實裂縫中的運移情況，如研究的時間和空間尺度與實際生產(chǎn)情況相差太大，邊界條件設置與實際情況相差甚遠，監(jiān)測統(tǒng)計流體、顆粒兩相速度的手段不夠完善等[12]。通過利用數(shù)值計算手段則可以彌補實驗研究的不足，結合實驗研究總結得到的校正模型，能夠精確地描述大尺度下支撐劑的運移情況[13]。因此筆者通過數(shù)值計算研究不同參數(shù)對CO2壓裂液攜砂效果的影響，為超臨界CO2壓裂設計優(yōu)化及現(xiàn)場施工提供重要的技術支撐。

1 數(shù)學模型

由于歐拉-歐拉型方法將顆粒相整體視作單一連續(xù)相的假設與實際情況相差較大，而歐拉-拉格朗日型粒子類方法中的多相質(zhì)點網(wǎng)格法（MP-PIC）則是通過顆粒相壓力梯度力替代顆粒之間的碰撞，從而表征顆粒間的相互作用力；將具有相同動力學特性的顆粒群作為計算單元，以便達到減小計算量的目的，對于大規(guī)模粒子運動模擬計算具有獨特的優(yōu)勢[14]。MP-PIC方法中顆粒相的運動不存在因差分格式導致在空間上產(chǎn)生的非物理擴散，因此其運動描述更為準確[15]?？紤]到CO2黏度低，裂縫中流體運動復雜，存在湍流效應，以及由于CO2攜砂能力差，導致顆粒之間的作用以及顆粒-固壁作用不可忽略，因此采用多相質(zhì)點網(wǎng)格模型（MPPIC）方法模擬顆粒團在CO2壓裂液中的運移沉降過程。

1.1 流體控制方程

基于Patankar and Joseph[16]擴展的二維MPPIC方法，考慮大量支撐劑顆粒在裂縫擴展中的運移，流體相的控制方程可以通過Favre平均法推導得出，鑒于篇幅限制，具體過程不再贅述，具體可參考文獻[17]。流體控制方程為體積平均納維-斯托克斯方程（VANS），包括連續(xù)性方程和動量方程[15,18]：

式中，νf為流體運動黏度，其與局部顆粒濃度有關，m2/s。

由于顆粒的存在，攜砂液的有效黏度高于純流體，因此該部分使用Barree和Conway（1994）的黏度模型：

式中，v0是純流體的運動黏度，m2/s；φcp是緊密堆積顆粒體積分數(shù)，%。

1.2 顆粒動力學方程

顆粒群的運動描述采用代表質(zhì)點法的思想，即一個質(zhì)點代表多個顆粒，而這些顆粒具有相同的性質(zhì)如密度、大小等，它們具有相似的運動行為，因此可以看作一個粒子團，運移過程中質(zhì)量保持不變。質(zhì)點運動由牛頓第二定律控制，其控制方程[20–21]可表示為：

采用Snider[22]的粒子法向應力模型來解決粒子碰撞和緊密堆積問題，τp可表示為：

式中，ps和 β 為模型參數(shù)，Pa和無量綱；φcp是緊密堆積時的顆粒體積分數(shù)，%。這里ps和 β分別設置為4 Pa和2[22]。

采用基于BGK近似模型的一個粒子分布函數(shù)輸運方程來對描述碰撞阻尼機制[21]，為了考慮壁摩擦效應，這里只考慮顆粒緊密堆積區(qū)域。為粒子相的質(zhì)量平均速度，表示為可以表示為：

式中，α為顆粒速度分布的質(zhì)量加權標準差，無量綱；η為模型系數(shù)，無量綱；r32為Sauter平均半徑，m，可表示為[15]：

式中，rp是粒子半徑，m；在顆粒均勻情況下，Sauter平均半徑是顆粒相的平均半徑。

流體-顆粒之間存在相互作用，其中曳力最為重要。這里采用Wen-Yu模型[23]，則曳力系數(shù)可表示為：

式中，r是顆粒半徑，m；為顆粒與流體密度比，無量綱；Cd為無量綱的曳力系數(shù)，表示為：

1.3 求解方法

流體部分求解采用SIMPLE壓力校正迭代方法進行求解，時間項離散采用歐拉格式，對流項采用顯式一階迎風格式。顆粒部分與顆粒速度相關的項均采用隱式離散，其余項如壓力，顆粒應力等都由流場數(shù)據(jù)顯式給出。由于MP-PIC為混合型方法，即流體描述為歐拉型方法，而顆粒描述采用的是拉格朗日型方法，二者之間的數(shù)據(jù)交換需要通過插值得到，采用二維雙線性插值。

數(shù)值求解過程中裂縫擴展、顆粒團運動和流體運動3個過程是通過順序和顯示耦合實現(xiàn)的，即假設流體和顆粒運動對裂縫擴展沒有影響，首先將裂縫寬度嵌入到控制流方程中，以考慮裂縫傳播效應；然后，流體-顆粒相互作用是完全耦合的，為簡化計算，只考慮曳力和流體-顆粒相互作用力等。具體計算程序如圖1所示，其中上面的箭頭表示3種物理之間的數(shù)據(jù)傳輸。上標“n”和“n+1”分別表示第n個和（n+1）個時間步；q濾失是 流體濾失量；w是 裂縫寬度；φp是顆粒的體積分數(shù)；分別是流體和顆粒速度。

圖1 計算程序流程圖

2 實驗對比驗證

通過與實驗結果對比驗證方法的可行性，隨后利用到大規(guī)模問題當中，由簡單邊界到復雜邊界條件出發(fā)，逐步考慮在以上不同因素組合下支撐劑運移的最終分布，并作相關參數(shù)敏感性分析，找到影響支撐劑運移效果的關鍵因素。

2.1 物模實驗

物模實驗采用可視化平面裂縫模擬實驗系統(tǒng)，實驗系統(tǒng)主體為可視化平面裂縫模型，為透明鋼化玻璃，兩端設置有模擬井筒，配合流體泵和混砂罐，用于模擬支撐劑在垂直裂縫中運移及鋪置。

2.1.1 實驗方案

以滑溜水壓裂液為實驗對象，使用40/70目陶粒支撐劑（視密度為2800 kg/m3），在排量為8 m3/min，并設定砂比為6%的條件下進行物理模型實驗，實驗過程中保證均勻加砂，通過記錄某些支撐劑的運移軌跡和時間，獲得支撐劑運動速度，并對砂堤高度進行記錄。

2.1.2 實驗結果

在裂縫模型中出現(xiàn)支撐劑時開始計時，每隔可知，隨著攜砂液不斷流入裂縫，不同時間所形成砂堤的形狀大體一致，但砂堤高度隨著時間逐漸增大。通過圖2不難看出，裂縫模型總長L=4 m，裂縫各部位的沉砂量也差異較大，約95%的支撐劑均鋪置于裂縫長度L=0～3 m的裂縫區(qū)間位置處，即近縫口位置的3/4裂縫長度內(nèi)，僅有約5%的支撐劑鋪置于近裂縫尖端1/4裂縫長度內(nèi)。

圖2 實驗系統(tǒng)示意圖

2.2 數(shù)模對比驗證

為便于分析，將支撐劑運移在垂直裂縫剖面上的砂堤高度分為2類，即支撐劑沉降區(qū)和支撐劑懸浮區(qū)，分別取支撐劑體積分數(shù)為40%和10%為分界線[12]，即體積分數(shù)大于40%時支撐劑完全沉降，認為其為純固相；體積分數(shù)在10%～40%范圍內(nèi)支撐劑處于懸浮狀態(tài)，認為其為固液混合項；支撐劑體積分數(shù)小于10%時可以認為其為純液相。

數(shù)值計算中參數(shù)均與物模實驗參數(shù)一致，唯一不同的是在處理粒徑上采取40/70目的平均粒徑作為代表。時間節(jié)點分別為0.5，1，1.5，2 min的數(shù)值模擬結果如圖3所示。

從圖3可以看出，在砂堤長度上，大部分支撐劑均處于距入口3 m以內(nèi)，這與實驗結果圖3也比較吻合。圖4為不同時刻室內(nèi)實驗和數(shù)值計算的砂堤峰值高度對比，不難看出其結果相差不大。綜合比較二者結果，差異較小，說明了該模型及數(shù)值計算的正確性和可靠性。

圖3 砂比為6%時實驗砂堤形態(tài)分布圖

圖4 砂比為6%時砂堤形態(tài)數(shù)值計算結果

3 超臨界CO2壓裂縫內(nèi)支撐劑運移規(guī)律分析

二氧化碳物性參數(shù)隨溫度和壓力的變化而變化，其相變情況復雜，在CO2壓裂過程中其地面管線流動、井筒流動等溫壓變化較大的過程，研究多基于CO2物性模型進行耦合計算分析[24]。而本文研究的為CO2壓裂過程中其進入地層后在裂縫擴展過程中的支撐劑運移規(guī)律，考慮到CO2在裂縫中的平均流速較小，且可以認為二氧化碳的密度和壓力與裂縫中運移過程中的地層應力相比變化不大，為便于計算，本文計算假設CO2是弱可壓縮性，主要考慮二氧化碳的平均黏度和密度對支撐劑運移的影響。

3.1 考慮流體濾失的簡單矩形裂縫區(qū)域內(nèi)的支撐劑運移分析

為了便于對比分析，考慮CO2流體濾失，不考慮裂縫擴展和裂縫寬度變化，考察區(qū)域簡化為一個矩形區(qū)域，長160 m，高20 m，裂縫寬度設為3 mm，如圖5所示。從左側均勻速率注入，右側定壓，假設整個矩形裂縫區(qū)域內(nèi)流體濾失速率一致，最終使得右側出口流體速率由于濾失作用接近于零，計算分析支撐劑在一矩形裂縫區(qū)域內(nèi)的運移情況。

圖5 不同時刻實驗和數(shù)值計算的砂堤峰值高度對比圖

假設CO2在攜砂過程中密度為600 kg/m3，初始黏度為0.1 mPa·s，假設增黏后黏度為10 mPa·s，支撐劑密度為2650 kg/m3和1000 kg/m32種，支撐劑尺寸為20/40目和40/70目兩種，排量5 m3/min，即入口速率約為0.2 m/s；因此本算例共設計支撐劑密度、尺寸和攜砂液黏度3個參數(shù)分別各2個取值正交共8組數(shù)值實驗。計算時間共30 min，計算結果如圖6所示，支撐劑沉降區(qū)和支撐劑懸浮區(qū)分別取支撐劑體積分數(shù)為40%和10%為分界線。從圖6可以看出，當流體黏度較低時，即黏度為0.1 mPa·s時，均會在近井區(qū)域形成砂床。由于未經(jīng)增黏，液態(tài)二氧化碳的黏度比清水還要低很多，顆粒雷諾數(shù)非常大，導致顆?？焖僭诮畢^(qū)域沉降堆積形成砂床。由于黏度極低，優(yōu)化其他參數(shù)對結果影響不大。這些近井區(qū)域支撐劑堆積預示著著實際壓裂施工中的砂堵現(xiàn)象，泵壓會急劇升高，說明攜砂效果極差。

圖6 簡單矩形裂縫示意圖

將二氧化碳增黏100倍后，攜砂效果有明顯改善。但是僅增黏，效果也不是很理想。當二氧化碳黏度為10 mPa·s時，相比于支撐劑直徑為0.6 mm和密度為2650 kg/m3時，減小支撐劑尺寸，即直徑為0.2 mm時，或者減小支撐劑密度，即密度為1000 kg/m3時，均可以提高攜砂效果。但優(yōu)化支撐劑密度對攜砂效果的提高更為明顯。

不難發(fā)現(xiàn)，將二氧化碳增黏，同時減小支撐劑密度和尺寸，即二氧化碳黏度為10 mPa·s，支撐劑密度為1000 kg/m3，支撐劑直徑為0.2 mm時，其攜砂效果最理想，僅在裂縫前端沉降堆積，且支撐劑懸浮區(qū)大面積覆蓋了裂縫。

通過分析可以發(fā)現(xiàn)，顆粒密度會影響收尾速度，改變流場，一般密度越大顆粒收尾速度大沉降越快，然而在顆粒尺寸小的配置下，導致一方面沉降速度小而改變流場能力凸顯出來，低黏條件下易產(chǎn)生異重流；顆粒粒徑會影響收尾速度，粒徑越大沉降越快，易形成砂床；攜砂液黏度會影響顆粒收尾速度及流場受顆粒群擾動能力，黏度越大，顆粒沉降越慢，同時流場受到壁面黏性力增大，即使是大密度顆粒也難以影響流場。

3.2 考慮流體濾失和裂縫擴展的PKN型裂縫中的支撐劑運移分析

考慮CO2流體濾失，基于PKN型裂縫擴展模型[15,25,26]，計算分析CO2壓裂裂縫中的支撐劑的運移情況，充分考慮了實際復雜邊界條件。由于篇幅所限，本部分計算所用PKN型裂縫擴展模型不做詳述，具體參考詳見文獻[15,25,26]。在PKN型裂縫的假設中，裂縫是等高的，裂縫斷面為橢圓形，最大縫寬在裂縫中部，前緣不斷向前擴展，如圖7所示。

圖7 不同模擬條件下支撐劑運移結果

3.2.1 不同支撐劑及流體黏度影響分析

考慮一定排量下支撐劑粒徑、密度以及流體黏度影響，只考慮裂縫擴展，暫不考慮流體濾失的情況下。假設CO2在攜砂過程中密度為600 kg/m3，初始黏度為0.1 mPa·s，假設增黏后黏度為10 mPa·s，支撐劑密度為2650 kg/m3和1200 kg/m32種，支撐劑尺寸為20/40目和40/70目2種，排量5 m3/min，即入口速率約為0.2 m/s；因此本算例共設計支撐劑密度、尺寸和攜砂液黏度3個參數(shù)分別各兩個取值正交共8組數(shù)值實驗。計算時間共30 min，計算結果如圖8所示，支撐劑沉降區(qū)和支撐劑懸浮區(qū)分別取支撐劑體積分數(shù)為40%和10%為分界線。

圖8 PKN型裂縫示意圖

由圖8可以，看出支撐劑在擴展裂縫中的運移行為大致與在固定邊界的矩形裂縫中的結果類似，但支撐劑在豎直方向上的砂堤濃度分布比簡單裂縫中更為復雜，尤其支撐劑懸浮區(qū)形狀在豎直方向上縫前端部更高，也說明了考慮裂縫擴展后，在裂縫高度方向上中部的速度較大。增黏后支撐劑沉降區(qū)域面積減少，說明支撐劑沒有太多沉降，當支撐劑密度為1200 kg/m3，支撐劑直徑為0.2 mm，二氧化碳黏度為10 mPa·s時，支撐劑大部分處于懸浮區(qū)，攜砂效果是最好的，這與固定邊界的矩形裂縫中的最優(yōu)結果也是一致的。

3.2.2 不同排量影響分析

考慮到實際增黏效果達不到10 mPa·s這個量級，根據(jù)實際情況，假設二氧化碳增黏后黏度為2.5 mPa·s，顆粒尺寸為40/70目，密度取1200 kg/m3，排量分別為1、3、5 m3/min的3組算例以作比較。不同排量下支撐劑運移的數(shù)值計算結果如圖9所示。由圖9可以發(fā)現(xiàn)，增加黏度至2.5 mPa·s后，采用超輕支撐劑與細尺寸顆粒組合，攜砂效果仍然十分理想，與10 mPa·s情況相差不大，能夠滿足工程需求。此外，值得注意的是，排量為3 m3/min與5 m3/min的結果從顆粒分布上看差異不大。這種情況說明在該排量范圍下，支撐劑的平面鋪置效果類似，而當排量過小時，支撐劑的平面鋪置效果亦顯然不如高排量的效果。

圖9 PKN裂縫中的支撐劑運移結果

3.2.3 流體濾失影響分析

為了考慮流體濾失的影響，這里分別針對小尺寸低密度支撐劑在黏度分別為0.1 mPa·s和2.5 mPa·s 2種情形下進行模擬。參考壓裂方案設計模擬計算結果，當黏度為0.1 mPa·s時流體濾失率為20.53%，黏度為2.5 mPa·s時流體濾失率為6.26%。計算結果如圖10所示。

圖10 增黏至2.5 mPa·s時不同排量下支撐劑運移結果

由圖10可以發(fā)現(xiàn)，在當前濾失率條件下，支撐劑運移情況與不考慮流體濾失情況相差很小。黏度為2.5 mPa·s時，由于濾失率本身很小，為6%左右，因此運移結果十分接近；黏度為0.1 mPa·s時，盡管流體濾失20%以上，但是由于該情況下支撐劑主要在近井范圍堆積，流體濾失對顆粒運移的效果在近井周圍的影響并不明顯，因此依然與不考慮濾失情況下支撐劑運移結果相似。

4 結論與建議

考慮CO2壓裂過程中攜砂液中支撐劑顆粒之間相互作用，采用歐拉-拉格朗日方法中的多相質(zhì)點網(wǎng)格方法，建立超臨界CO2壓裂支撐劑運移數(shù)學模型，通過室內(nèi)支撐劑運移物模實驗驗證模型準確性，最后對超臨界CO2壓裂支撐劑運移規(guī)律計算和分析，并給出相應的建議。

1.CO2壓裂液的增黏是不可或缺的。未增黏CO2由于黏度低，攜砂效果極差，優(yōu)化其他參數(shù)對攜砂效果影響不大；CO2黏度增加到2.5 mPa·s可有效提高攜砂效果；采用超輕支撐劑與細尺寸顆粒組合，攜砂效果與增黏到10 mPa·s效果相差不大。

2.減小支撐劑密度和尺寸可以提高CO2壓裂液攜砂效果。優(yōu)化支撐劑密度比尺寸對攜砂效果的影響更為明顯；CO2黏度較低時，減小支撐劑密度和粒徑可以提高攜砂效果，在攜砂液充分增黏的前提下可以選擇目數(shù)稍大的以有效支撐裂縫。在施工成本允許的前提下可以利用一部分超輕支撐劑，與常規(guī)支撐劑交替注入，從而同時保證裂縫上端與低端同時支撐，提高最終的支撐面積。

3.增大排量可以提高攜砂效果，但排量繼續(xù)增大，其攜砂效果變化較??；工程上不需要只為了提高攜砂效果而過大的增加排量，導致摩阻增大，造成井口超壓和成本升高。

4.流體濾失對CO2攜砂效果影響不大；只要滿足攜砂和成本要求，不需要通過對CO2過度增黏來減少流體濾失。從CO2壓裂與地層埋存角度出發(fā)，CO2流體在裂縫擴展過程中的濾失是有益的。