一種應用于射流破巖的孔隙形狀近似法

馬小晶,周新超,肖新朋,程 璨

(新疆大學 電氣工程學院,烏魯木齊 830017)

隨著常規(guī)油氣資源的減少以及能源需求的增長,石油天然氣的開發(fā)逐漸從常規(guī)油氣資源轉(zhuǎn)移到了非常規(guī)油氣資源[1]。由于我國蘊藏著較豐富的砂巖致密油(發(fā)展致密油的開采技術(shù)是十分必要的,且致密油的存儲規(guī)律和存儲形態(tài)與巖體內(nèi)部的孔隙結(jié)構(gòu)及孔隙間的相互連通性有著密切的關系[2],因此,需要對含有孔隙結(jié)構(gòu)的巖體在常規(guī)作業(yè)中產(chǎn)生的影響展開深入研究[3]。

高壓水射流能夠大面積且有導向的碎裂儲層巖體,方便致密油等油氣資源的解吸與滲流,已被廣泛應用于鉆井開采領域[4]。眾多學者采用不同的數(shù)值模擬方法研究了各類射流破巖問題[5],但針對射流沖擊含孔隙巖體的破損過程以及孔隙結(jié)構(gòu)對破巖效果的研究相對較少。為了深入探究含孔隙巖體的破損機制,需要構(gòu)建能夠真實反映巖體中復雜孔隙結(jié)構(gòu)和分布的模型,為進一步提高巖體破損問題的模擬精度提供途徑。

目前,國內(nèi)外學者已對孔隙結(jié)構(gòu)的模型構(gòu)建技術(shù)及相關問題展開了大量研究。隋微波等[6]將真實的孔隙抽象簡化為橢球和球體,研究了三維抽象孔隙模型中孔隙結(jié)構(gòu)特征對孔隙體積壓縮系數(shù)的影響;蔡少斌等[7]為探究油氣資源在巖石多孔介質(zhì)內(nèi)部的遷移過程,建立了孔隙空間模型,并對多相流體的運動進行了模擬研究;李俊鍵等[8]提取了孔隙網(wǎng)絡結(jié)構(gòu)模型,分析了微觀孔喉結(jié)構(gòu)的非均質(zhì)性對剩余油分布形態(tài)的影響。為了構(gòu)建更加接近真實巖體結(jié)構(gòu)的數(shù)值模型,李博等[9]采用數(shù)字巖心技術(shù)(digital core technology,DCT)生成了含裂隙的巖體三維微觀結(jié)構(gòu)模型,計算了不同孔隙度下基質(zhì)的介電性質(zhì)。數(shù)字巖心技術(shù)能夠生成非常接近真實巖體的微觀孔隙結(jié)構(gòu),但需要具體的巖心樣本且生成步驟較復雜,并存在計算困難等問題[10]。

本文基于光滑粒子流體動力學(smoothed particle hydrodynamics,SPH)方法,提出一種孔隙形狀近似(pore shape approximation,PSA)法,建立了含孔隙結(jié)構(gòu)巖體的數(shù)值模型,通過模擬孔隙巖體的單軸壓縮試驗,驗證了該孔隙結(jié)構(gòu)生成方法的有效性;在此基礎上,引入描述射流和巖體力學特性的NULL型和H-J-C型本構(gòu)關系,構(gòu)建出射流沖擊巖體的動力學模型,并將模擬結(jié)果與試驗結(jié)果進行對比分析;之后模擬研究了射流沖擊含孔隙結(jié)構(gòu)巖體的動態(tài)過程,并將孔隙形狀近似法與數(shù)字巖心技術(shù)生成巖體模型的破損過程進行了對比分析;最后著重研究了含不同的孔隙形狀、孔隙尺寸及各種形狀孔隙占比的巖體在射流沖蝕下的破損規(guī)律。

1 SPH方法的基本理論

SPH方法作為一種純拉格朗日的無網(wǎng)格粒子法,是將計算域離散為一系列可運動的粒子,粒子的物理屬性可以通過求解非穩(wěn)態(tài)N-S控制方程得到,適用于模擬沖擊作用下的大變形和破損現(xiàn)象。粒子的場函數(shù)可通過周圍粒子的函數(shù)進行加權(quán)近似求解,即[11]

(1)

式中:f(x)為坐標向量x的函數(shù);Ω為x的支持域;x-x′為粒子間距離;h為光滑長度;W是核函數(shù),即:

(2)

式中,d為空間維數(shù)。用粒子近似法將連續(xù)形式積分方程寫為離散形式的方程,即:

(3)

式中:ρi和mi分別為粒子i的密度和質(zhì)量。

2 巖體內(nèi)孔隙模型的構(gòu)建方法

2.1 孔隙形狀近似法

孔隙的形狀、尺寸以及連通性等會直接影響開采效率,但巖體內(nèi)部的孔隙結(jié)構(gòu)十分復雜,數(shù)值建模時常會盡可能的簡化孔隙結(jié)構(gòu)。Silin等[12]將巖體中的孔隙和孔隙連通喉道分別近似為球體和棍狀,建立了微觀孔隙結(jié)構(gòu)的數(shù)值模型;Sadowsky等[13-14]將孔隙形狀視為圓柱體和狹長裂縫等,建模時均抽象為橢球體。

為了構(gòu)建更加貼近真實的巖體內(nèi)部孔隙結(jié)構(gòu)的模型,本文掃描了4塊砂巖樣本,采用的CT型號為 phoenix v|tome| x m,掃描精度的分辨率為30 μm,巖石樣本的直徑為30 mm,高度50 mm。每塊樣本得到175張圖片,通過處理和分析上述圖片研究了砂巖內(nèi)部的孔隙形狀。圖1給出了砂巖樣本和孔隙近似形狀的示意圖,砂巖內(nèi)部的孔隙形狀及分布都較為復雜,不可能只存在某單一形狀的孔隙,鑒于此,本文提出一種基于SPH方法的孔隙形狀近似法,在前人將孔隙簡化為球體和橢球體的基礎上,結(jié)合掃描圖增加了其他三種近似形狀,即:薄片體、正三棱柱和正三棱錐,并將孔隙喉道近似為棍狀。

(a) 砂巖樣本

2.2 離散孔隙的構(gòu)建原理

由于研究采用了適用于模擬大變形和破損現(xiàn)象的SPH方法,因此本文將射流、巖體及其內(nèi)部的孔隙結(jié)構(gòu)均離散為一系列的粒子,圖2給出了不同形狀的孔隙的模型示意圖,圖中θ代表孔隙尺寸的指標。

圖2 不同形狀的孔隙模型示意圖Fig.2 Model of pores with different shapes

研究中采用了蒙特卡洛法隨機生成孔隙[15],孔隙均勻分布在巖體內(nèi)部,并且允許孔隙之間存在部分重疊,同時在孔隙之間還構(gòu)造了多條起連通作用的孔隙喉道,使構(gòu)建的孔隙結(jié)構(gòu)更加接近天然巖體的形成過程。圖3給出了幾種形狀的孔隙及孔隙喉道在巖體內(nèi)部的分布方向。

(a) 薄片體和橢球體

LS-DYNA軟件中生成巖體模型的步驟為:

步驟1沿X軸排布一列粒子,改變Z軸位移,繼續(xù)沿X軸排布粒子,如此反復,最終構(gòu)成一個粒子面;

步驟2改變Y軸位移,按步驟1沿Y軸生成一面粒子,重復上述過程,最終生成一個由SPH粒子構(gòu)成的三維巖體,存儲在K文件中。

本研究通過K文件獲得巖體模型中粒子的坐標位置,提取部分巖體粒子并替換為孔隙粒子,即可生成不同形狀的孔隙,構(gòu)建出含孔隙結(jié)構(gòu)的巖體模型。以正方體孔隙為例,生成過程如圖4所示,步驟為:

圖4 生成正方體孔隙的示意圖Fig.4 Schematic diagram of generating cube pores

步驟1使用隨機函數(shù)抽取巖體粒子的坐標位置,確定出生成孔隙的基準粒子Pi的坐標信息;

步驟2按照粒子的間隔距離,在基準粒子的位置上分別沿X軸、Y軸和Z軸三個方向增加2個粒子,確定出構(gòu)成正方體結(jié)構(gòu)的8個粒子坐標,之后將這8個粒子的屬性替換為孔隙粒子,從而形成正方體孔隙;

步驟3采用隨機函數(shù)抽取粒子位置,重復上述步驟1和2不斷生成孔隙,直到孔隙個數(shù)達到要求。

2.3 孔隙形狀近似法的驗證

為了驗證本文提出的孔隙形狀近似法的有效性,采用SPH方法和文獻[16]的砂巖的單軸壓縮試驗數(shù)據(jù)對該試驗進行了模擬,模型尺寸3 mm×3 mm×3 mm的正方體,孔隙特性等參數(shù)參見文獻[16],巖體上表面施加位移荷載,底面建立剛體網(wǎng)格并施加固定約束。圖5為孔隙巖體Von-Mises應力分布云圖,由圖5可知,由于巖體內(nèi)部有孔隙結(jié)構(gòu),巖石內(nèi)部的Von-Mises應力分布并不均勻,并且在孔隙附近存在應力集中現(xiàn)象,這與文獻[16]中采用網(wǎng)格法得到的模擬結(jié)論一致。

文獻模型[16](t=0.2 ms)

圖6為經(jīng)過數(shù)值計算得到的孔隙巖體的應力應變曲線,與文獻[16]中的參考模型擬合曲線相比,本文在前期的上升趨勢保持一致,并且差距較小。根據(jù)該曲線可得出本文孔隙巖體模型的整體彈性模量為31.987 GPa。表1給出了孔隙巖體的整體彈性模量,由表1可知,文獻[16]采用數(shù)字巖心方法和本文提出的孔隙形狀近似法構(gòu)建的孔隙巖體模型的彈性模量均小于無孔隙砂巖的彈性模量,這說明巖體內(nèi)部的孔隙結(jié)構(gòu)在一定程度上會影響其力學特性,與文獻[17]的研究結(jié)論吻合。因此,為了分析更加接近實際的射流破巖動態(tài)過程,本文在巖體中引入孔隙結(jié)構(gòu)。與文獻[16]中數(shù)字巖心方法構(gòu)建的孔隙巖體模型相比,本研究的整體彈性模量更加接近實際值,這表明基于孔隙形狀近似法建立的孔隙巖體模型更符合實際,應用于數(shù)值模擬可獲得更加精確的研究結(jié)果。

表1 孔隙巖體的整體彈性模量Tab.1 Overall elastic modulus of porous rock

圖6 孔隙巖體的應力應變曲線Fig.6 Stress-strain curve of porous rock

3 孔隙形狀近似法在射流破巖的應用

3.1 巖體和射流的本構(gòu)關系

本文構(gòu)建了基于SPH方法的水射流破巖模型,水射流采用NULL本構(gòu)關系進行描述,并賦予Gruneisen狀態(tài)方程[18],即:

(γ0+aμ)Ea

(4)

式中:C為沖擊波速度與質(zhì)點速度變化曲線截距;S1、S2、S3為沖擊波與質(zhì)點速度變化曲線的斜率系數(shù);γ0為Gruneisen常數(shù);a為一階體積修正量。具體參數(shù)如表2所示[19]。

表2 水射流的參數(shù)Tab.2 Parameters of water jet

研究選用了存儲致密油的砂巖,材料參數(shù)如表3所示[20],滿足H-J-C本構(gòu)關系,屈服面方程為[21]

表3 砂巖的參數(shù)Tab.3 Parameters of sandstone

σ*=[A(1-D)+Bp*N](1-Clnε*)

(5)

式中:σ*為巖體沖擊載荷下等效應力與靜態(tài)屈服強度之比;P*為巖體實際壓力與靜態(tài)屈服強度之比;ε*為沖擊載荷下巖體應變率與靜態(tài)應變率之比;A、B、N和C為材料的強度參數(shù),D為損傷因子。

3.2 射流沖擊含孔隙巖體的SPH模型

針對研究的射流沖擊含孔隙巖體問題,本文采用LS-DYNA軟件進行分析。模型簡化為水射流、含孔隙結(jié)構(gòu)的巖體和孔隙物質(zhì)三種物質(zhì),暫不考慮氣相環(huán)境的影響。研究采用了1/2對稱模型,對稱面施加SPH_SYMMETRY_PLANE約束,并對巖石周圍與底面施加SPH_NON_REFLECTING約束模擬巖體的無限大邊界,同時約束巖石的豎向移動。對稱后的巖體為10 mm×5 mm×5 mm的長方體,射流直徑為1 mm,長度為5 mm。模型共由691 460個粒子構(gòu)成,其中5 460個為射流粒子,孔隙及孔隙喉道內(nèi)部填充水,其本構(gòu)關系和材料參數(shù)與水相同。圖7為射流沖擊孔隙巖體的模型示意圖。

圖7 射流沖擊孔隙巖體的模型示意圖Fig.7 Model of jet impacting porous rock

3.3 模型的驗證與分析

(1) 射流破巖問題

為了驗證本文建立的射流沖擊巖體模型的有效性,對Xue等[22]研究的水射流沖擊無孔隙煤巖問題進行了模擬,各項參數(shù)均與文獻設置保持一致,巖體部分的模型是半徑25 mm,高度50 mm的圓柱體,射流直徑為1 mm,長度為50 mm。

圖8和圖9分別給出了t=0.06 ms時煤巖的破損坑形狀與破損坑深度隨沖擊時間的變化曲線。從圖8可以看出,本模型的破損坑深度、形狀與文獻[22]的模擬結(jié)果基本相同,符合試驗結(jié)果的沖孔形態(tài)。由圖9可知,隨著沖擊時間的增加,煤巖的破損坑深度不斷增大,與文獻[22]的模擬結(jié)果基本一致,驗證了本文建立的射流破巖模型的有效性。

本文模擬結(jié)果(t=0.06 ms)

圖9 煤巖破損坑深度隨時間的變化曲線Fig.9 Variation curve of coal rock damage pit depth with time

(2) 孔隙形狀近似法與數(shù)字巖心技術(shù)的對比

在上節(jié)所構(gòu)建的水射流沖擊無孔隙巖體模型的基礎上,采用孔隙形狀近似法建立了含孔隙巖體的數(shù)值模型,模擬研究水射流沖擊破巖的動態(tài)過程,并與數(shù)字巖心技術(shù)的模擬結(jié)果進行了對比分析。

圖10為模型的構(gòu)建流程,DCT是對砂巖的CT掃描圖片進行處理,通過Avizo等軟件提取孔隙結(jié)構(gòu),而PSA模型則只需根據(jù)掃描得到的孔隙度和孔隙尺寸等核心特性參數(shù)即可生成孔隙結(jié)構(gòu)。DCT所構(gòu)建的巖體模型數(shù)據(jù)來源于Berea砂巖,孔隙度為14.9%,本文采用PSA構(gòu)建出孔隙度和孔隙平均尺寸與其相同的巖體模型,兩種模型都采用了1/2對稱模型,巖體尺寸為10 mm×5 mm×5 mm的長方體,射流直徑為1 mm,長度為5 mm。

(a) DCT

圖11為上述兩種方法所構(gòu)建的孔隙巖體模型在水射流沖擊作用下的破損演化過程,圖中藍色、綠色和紅色分別為水、孔隙和巖體。由圖11可以看出,在射流的持續(xù)沖擊作用下,兩種巖體模型的破損坑均逐漸增大,破損形狀的變化過程基本一致。

(a) DCT

本文選取了破損坑橫截面的破損寬度指標w和破損坑縱截面的破損深度指標h以及破損面積指標s來反映巖體的破損情況,其中S是初始巖體縱截面的面積,s為該截面破損坑的面積,如圖12所示。

圖12 巖體破損指標的示意圖Fig.12 Diagram of rock damage index

圖13給出了數(shù)字巖心模型和孔隙形狀近似法模型的破損指標隨沖擊時間的變化規(guī)律。

(a) 破損寬度

由圖13可知,上述兩種模型的破損寬度除在t=0.012 ms和t=0.02 ms時存在較小差異外,其他沖擊時間的變化趨勢基本相同;兩種模型的破損面積的變化曲線則幾乎完全重合。由此可知,采用DCT和PSA所建的巖體模型可得到較為一致的破損規(guī)律,這說明本文所提的孔隙形狀近似法應用于射流破巖問題是可行的。

表4給出了孔隙形狀近似法和數(shù)字巖心技術(shù)所建模型的計算總時間。兩者的模型參數(shù)完全一致,沖擊總時長為0.021 ms,分別選取3種不同孔隙平均尺寸進行對比,數(shù)值計算處理器為英特爾i5-10400F。由表4可知,在相同的孔隙平均尺寸下,采用數(shù)字巖心技術(shù)生成的孔隙巖體模型的模擬耗時大于孔隙形狀近似法,尤其是孔隙平均尺寸較大時,數(shù)字巖心技術(shù)構(gòu)建的模型還會出現(xiàn)計算崩潰的現(xiàn)象。顯然,本文提出的孔隙形狀近似法構(gòu)建的模型計算效率高且穩(wěn)定性好。

表4 不同模型計算時間Tab.4 Calculation time of different models

4 孔隙對巖體破損影響的研究

將孔隙形狀近似法應用于射流破巖的模擬研究中,分析含不同的孔隙尺寸、孔隙形狀以及各形狀孔隙占比的巖體在射流沖蝕下的破損規(guī)律。巖體模型的參數(shù)均取尺寸為10 mm×5 mm×5 mm,孔隙度為21%(19%為孔隙,2%為孔隙喉道),射流速度為500 m/s,沖擊總時間為0.21 ms,取破損面積為巖體破損指標。

4.1 孔隙的尺寸和形狀

為了分析孔隙尺寸和形狀對射流作用下巖體破損規(guī)律的影響,本文對內(nèi)部僅含單一形狀孔隙和孔隙喉道的巖體的射流破損過程進行了模擬研究。通過研究4塊砂巖樣本的實際孔隙結(jié)構(gòu),確定了巖體模型中各個形狀孔隙尺寸指標θ(參見本文2.2節(jié)的圖2)的范圍,球體和橢球體孔隙的θ取值分別為0.214 mm、0.285 mm,0.357 mm、0.428 mm、0.5 mm、0.571 mm、0.642 mm和0.714 mm,薄片體、正三棱柱和正三棱錐的θ取值分別為0.142 mm、0.285 mm,0.357 mm、0.428 mm、0.5 mm、0.571 mm和0.642 mm,其中橢球體和薄片體的縱橫比定義為長軸/短軸,縱橫比的取值范圍1.3～2.2。

不同的孔隙形狀和尺寸會直接影響巖體的破損程度,圖14給出了t=0.021 ms時含不同孔隙形狀和尺寸巖體的破損情況。由圖14可知,由于孔隙在巖體內(nèi)部隨機分布,巖體會呈現(xiàn)出不同破損形狀,孔隙分布較多的巖體區(qū)域在射流作用下更容易破損。從圖14還可知,當0.142 mm<θ<0.375 mm時,在相同的孔隙尺寸下,孔隙形狀對巖體破損程度的影響較為接近;當θ>0.375 mm時,孔隙形狀為薄片體、正三棱柱的巖體破損程度大于球體、橢球體的巖體。

圖15給出t=0.021 ms時不同尺寸孔隙的破損指標對比。由圖15(a)可知,當θ<0.571 mm時,薄片體孔隙的破損寬度最大,且在θ=0.375 mm時達到最大值。從圖15(b)可知,當θ=0.428 mm時,含橢球體、薄片體和正三棱柱孔隙巖體的破損深度較大,其中薄片體孔隙巖體的破損深度最大;當θ=0.5 mm和0.571 mm時,含薄片體、正三棱柱和正三棱錐的孔隙巖體破損深度較大。由圖15(c)可知,當0.142 mm<θ<0.375 mm時,在相同孔隙尺寸下,孔隙形狀對巖體破損面積的影響較為接近;當θ>0.375 mm時,含薄片體、正三棱柱的孔隙巖體的破損面積大于球體、橢球體的巖體;當θ=0.571 mm時,正三棱錐的孔隙巖體破損面積最大。

(a) 破損寬度

綜上,薄片體孔隙在特定的尺寸下,破損寬度、破損深度和破損面積遠大于其他形狀的孔隙;當孔隙尺寸較大時,孔隙形狀為薄片體、正三棱柱和正三棱錐的巖體在射流沖擊作用下更容易破損,不同形狀的孔隙都存在使巖體破損面積達到最大的孔隙尺寸。

4.2 不同形狀的孔隙在巖體中的占比

本文分析了不同形狀孔隙占比對射流破巖效果的影響,構(gòu)建了6種孔隙巖體模型,孔隙平均尺寸保持一致。巖體模型A、B、C和D是基于4塊砂巖樣本的孔隙特性進行構(gòu)建,巖體模型E和F則根據(jù)文獻[6]中Berea砂巖和Bentheimer砂巖進行構(gòu)建,孔隙形狀占比如表5所示。

表5 模型中的孔隙形狀占比Tab.5 Proportion of pore shape in the model

圖16為t=0.021 ms時不同巖體模型在射流沖擊作用下的破損效果。由圖16可知,不同巖體模型的破損坑形狀有明顯的差異。圖17給出了不同模型巖體的破損面積隨沖擊時間的變化曲線。由圖可知,模型C和E以及F的破損面積大于其他模型,其中模型E的破損面積最大,結(jié)合表5可以發(fā)現(xiàn),巖體模型C、E和F中薄片體、正三棱柱和正三棱錐孔隙占比相對較高,說明這類形狀孔隙的存在會使巖體更容易發(fā)生破損。

圖16 不同巖體模型的破損效果(t=0.021 ms)Fig.16 Damage effect of different models (t=0.021 ms)

圖17 不同模型下破損面積隨沖擊時間的變化曲線Fig.17 Variation curve of damage area with impact time under different models

5 結(jié) 論

(1) 針對巖體內(nèi)部的孔隙結(jié)構(gòu),提出一種用于SPH方法的孔隙形狀近似法,構(gòu)建了孔隙巖體的數(shù)值模型,通過模擬巖體的單軸壓縮試驗,驗證了含孔隙結(jié)構(gòu)巖體模型的有效性。

(2) 將孔隙近似法應用于基于SPH方法的射流破巖模擬研究中,并與數(shù)字巖心技術(shù)的模擬結(jié)果進行了對比分析,結(jié)果表明數(shù)字巖心技術(shù)和孔隙形狀近似法構(gòu)造的孔隙巖體模型具有相同的破損規(guī)律,并且孔隙形狀近似法具有較好的穩(wěn)定性和效率。

(3) 模擬研究了孔隙形狀、孔隙尺寸和不同形狀孔隙的占比對射流破巖的影響,研究結(jié)果表明,不同形狀的孔隙均存在一個使巖體破損最嚴重的孔隙尺寸;孔隙為薄片體、正三棱柱和正三棱錐,且孔隙的占比較高時,巖體則更容易被射流破損。針對目前在射流破巖領域中對含孔隙巖體的研究較少的情況,提出了一種應用于射流破巖領域的孔隙形狀近似法。