賈石川，胡立強(qiáng)，肖智勇，王剛,2，姜楓

(1.山東科技大學(xué) 山東省土木工程防災(zāi)減災(zāi)重點(diǎn)實(shí)驗(yàn)室，山東 青島 266590；2.山東科技大學(xué) 礦山災(zāi)害預(yù)防控制省部共建國家重點(diǎn)實(shí)驗(yàn)室培育基地，山東 青島 266590)

0 引 言

自從北美成功開發(fā)頁巖氣以來，頁巖氣等非常規(guī)能源受到了能源界的廣泛關(guān)注[1]。中國頁巖氣可采資源量為1.5×1013～2.5×1013m3,與美國資源量相當(dāng)，現(xiàn)采用大規(guī)模水力壓裂和微地震監(jiān)測(cè)技術(shù)進(jìn)行開采[2]。中國在頁巖氣開發(fā)探索中形成了新的開采機(jī)制，提高了單井產(chǎn)能[3]。非常規(guī)能源包括煤層氣、致密砂巖氣和頁巖氣,為了減少對(duì)常規(guī)能源的依賴，非常規(guī)能源的商業(yè)化開發(fā)變得越來越迫切，滲透率通常被認(rèn)為是開采的關(guān)鍵參數(shù)之一，是決定非常規(guī)能源商業(yè)開采可行性的重要因素[4]。實(shí)驗(yàn)室測(cè)量滲透率分為穩(wěn)態(tài)法和非穩(wěn)態(tài)法[5]，滲流介質(zhì)有煤油、水、氣體等多種[6]。穩(wěn)態(tài)法在不同驅(qū)替壓力下監(jiān)測(cè)流量，通過滑脫效應(yīng)改進(jìn)，得到絕對(duì)滲透率；非穩(wěn)態(tài)法中的脈沖衰減法根據(jù)脈沖壓力與時(shí)間的關(guān)系得到滲透率。在測(cè)量滲透率方法的發(fā)展過程中，由于脈沖法測(cè)量致密巖石滲透率(10-6～10 mD)[7]花費(fèi)時(shí)間少，精度高，一直被廣大學(xué)者采用和發(fā)展[8-13]。

脈沖法和穩(wěn)態(tài)法的區(qū)別在于脈沖法不需要測(cè)量流速，只需要記錄壓力隨時(shí)間的變化關(guān)系就可以測(cè)出滲透率[14]，同時(shí)脈沖法在高圍壓條件下進(jìn)行，更加符合現(xiàn)場(chǎng)情況。然而脈沖壓力響應(yīng)受到多種因素影響，包括巖樣滲透率、孔隙度、橫截面積、巖石體壓縮性、基質(zhì)壓縮性以及測(cè)量所用的上下游儲(chǔ)層體積等。頁巖滲透率極低，即便采用脈沖法測(cè)量也需要很長的時(shí)間，且脈沖法對(duì)氣密性要求很高。此外，孔隙度測(cè)量以及反映壓力變化的壓力傳感器精度也會(huì)影響壓力響應(yīng)，實(shí)驗(yàn)裝置的局限性也使得實(shí)驗(yàn)室研究各種因素對(duì)壓力響應(yīng)的影響變得困難。為了更好地分析壓力響應(yīng)，數(shù)值模擬方法更加實(shí)際可行。

LIN W[15]采用解析解和數(shù)值解求解了滲透率巖石的一維流體流動(dòng)微分方程，然而在解析解中并未考慮下游儲(chǔ)層的影響，其數(shù)值解直接采用TRUMP代碼，未給出詳細(xì)的推導(dǎo)；F.Civan[16]提出了用實(shí)用有限分析法求解微分方程的空間/時(shí)間解的數(shù)值方法，與基于分離變量和拉普拉斯變換的原始有限解析求解方法所產(chǎn)生的繁瑣離散化格式相比，可以得到更加簡便的離散化格式；王自明等[17]利用有限元與有限差分法求解流固耦合方程，交替迭代求解得到參數(shù)的解；劉文超等[18]通過有限差分法研究了不穩(wěn)定滲流邊界數(shù)值解，算法簡單，但局限于低滲非達(dá)西流等。

綜上所述，在多孔介質(zhì)的滲流計(jì)算問題中，有限差分法求解偏微分方程研究較少。本文在F.Civan[16]提出的實(shí)用有限分析法求解微分方程的基礎(chǔ)上，將其應(yīng)用于多孔可壓縮介質(zhì)的一維瞬態(tài)流動(dòng)的偏微分控制方程，提出用有限差分法求解非線性偏微分方程。此外，在MATLAB環(huán)境下編制計(jì)算程序代碼，模擬求解脈沖壓力隨時(shí)間和空間的變化趨勢(shì)。通過改變模型中的參數(shù)值，重點(diǎn)研究滲透率、孔隙度和上下游儲(chǔ)層體積對(duì)脈沖壓力響應(yīng)的影響，以期為非常規(guī)儲(chǔ)層能源商業(yè)化開采提供指導(dǎo)意見。

1 模型建立

在測(cè)量花崗巖滲透率時(shí)，W.F.Brace等[9]首先提出可壓縮流體通過多孔可壓縮介質(zhì)的一維瞬態(tài)流動(dòng)控制偏微分方程，將壓力簡化為距離x和時(shí)間t的函數(shù)，即

(1)

其中，L為試樣長度，c與地質(zhì)材料的儲(chǔ)存性能有關(guān)，是一個(gè)與儲(chǔ)層內(nèi)部連接的孔隙度、基質(zhì)和巖體壓縮性、流體性能有關(guān)的函數(shù)[19]。

利用達(dá)西定律和質(zhì)量守恒原理，試樣在上下游儲(chǔ)層的初始值可以簡化為

p(x,0)=p0,0

(2)

p(0,t)=pu(t),t≥0，

(3)

p(L,t)=pd(t),t≥0，

(4)

式中，pu(t),pd(t)分別為t時(shí)刻上下游儲(chǔ)層的壓力。

式(2)表示t為0時(shí)巖心內(nèi)部初始?jí)毫閜0，式(3)表示t時(shí)刻x=0處的壓力為pu(t)，式(4)表示t時(shí)刻x=L處的壓力為pd(t)，其中pd(0)=p0，pu(0)=p0+Δp。

老齡化為養(yǎng)老保障帶來巨大壓力。隨著人口平均預(yù)期壽命的上升，越來越多的人開始享受養(yǎng)老保障的福利，但上交老年保險(xiǎn)的比例不是很高，導(dǎo)致保險(xiǎn)行業(yè)盈利受阻，甚至負(fù)盈利[5]。人口老齡化會(huì)對(duì)醫(yī)療保險(xiǎn)和老年保險(xiǎn)形成沖擊，不能有效解決這個(gè)問題將給我國經(jīng)濟(jì)社會(huì)帶來負(fù)面影響。

上下游儲(chǔ)層的邊界值可以簡化為

(5)

(6)

(7)

(8)

(9)

式中：A,k,φ分別為橫截面面積，滲透率和孔隙度；Vu,Vd分別為上下游儲(chǔ)層體積；μ和β分別為測(cè)試流體的黏度和壓縮系數(shù)；βb為試樣體壓縮系數(shù)；βs為固體基質(zhì)的壓縮系數(shù)。

由于式(1)無法直接求解，且給出的解析解形式復(fù)雜，所以人們很難正確理解流體隨時(shí)間和空間變化的流動(dòng)行為。此外，邊界條件的復(fù)雜性也增加了解析解的難度。為了正確理解可壓縮流體通過多孔可壓縮介質(zhì)的一維瞬態(tài)流動(dòng)時(shí)空特性，采用有限差分法分析偏微分控制方程，并通過編制MATLAB程序?qū)⒎蔷€性偏微分方程代碼化，加入邊界條件和初始條件進(jìn)行數(shù)值計(jì)算。

壓力p(x,t)對(duì)x進(jìn)行泰勒多項(xiàng)式展開，

(11)

(12)

其中o(Δx3)為x的三階等價(jià)無窮小。將式(11)與(12)相加，

(13)

忽略x的高階等價(jià)無窮小，式(13)可改寫為

(14)

壓力p(x,t)對(duì)t進(jìn)行泰勒多項(xiàng)式展開，

(15)

忽略t的二階等價(jià)無窮小，將式(14)和(15)代入式(10)，則得到有限差分形式為

(16)

利用有限差分法計(jì)算巖心上下游邊界條件的壓力梯度[20]，

(17)

(18)

故上游邊界條件式(5)有限差分形式為

(19)

下游邊界條件式(6)有限差分形式為

(20)

對(duì)式(16)、式(19)～(20)通過MATLAB編程進(jìn)行時(shí)間步求解，可得到壓力隨時(shí)間和空間的變化曲線。

2 模型參數(shù)分析

非常規(guī)儲(chǔ)層中的一些參數(shù)如滲透率、孔隙度等，對(duì)于不同巖石其數(shù)量級(jí)相差很大，此外在實(shí)驗(yàn)室測(cè)量過程中，試樣尺寸、固體壓縮性、基質(zhì)壓縮性和上下游儲(chǔ)層體積大小均對(duì)壓力響應(yīng)具有一定影響。本文研究對(duì)象主要為煤巖，采用氦氣作為滲流介質(zhì)，數(shù)值模擬基本參數(shù)如表1所示。為了研究這些參數(shù)對(duì)壓力響應(yīng)的影響，基于建立的數(shù)值模型，在MATLAB環(huán)境下編程進(jìn)行計(jì)算，計(jì)算過程如圖1所示：首先根據(jù)表1輸入初始計(jì)算參數(shù)，選擇合適的時(shí)間步長dt與距離步長dx，計(jì)算中選取的最大時(shí)長為10 000 s，時(shí)間步長為1 s，L=10 cm，距離步長為0.05 cm；其次輸入上下游的邊界條件差分形式，即式(19)～(20)，最后輸入式(16)，得到時(shí)間域1：t_len-1與空間域1：x_len-1上的壓力值。

圖1 模型計(jì)算流程圖Fig.1 Flow chart of model calculation

表1 數(shù)值模擬基本參數(shù)Tab.1 Basic parameters used in numerical simulation

試樣βb、βs數(shù)值很小，計(jì)算過程中對(duì)c的影響可以忽略不計(jì)。假設(shè)氦氣壓縮系數(shù)和黏度系數(shù)保持不變，本部分僅研究滲透率、孔隙度和儲(chǔ)層體積3個(gè)因素對(duì)壓力響應(yīng)的影響。

2.1 滲透率對(duì)脈沖壓力響應(yīng)的影響

數(shù)值模擬的介質(zhì)為煤巖，滲透率為0.000 1～100 mD[22]，分別選0.1，0.5，1，2 mD進(jìn)行計(jì)算,結(jié)果如圖2所示。由圖2(a)可知，隨著滲透率降低，壓力減小的速度變緩，表明到達(dá)最終平衡的時(shí)間增長。由圖2(b)可以看出，不同滲透率下的壓力響應(yīng)均包含線性區(qū)域和非線性區(qū)域，且非線性區(qū)域出現(xiàn)的時(shí)間較短。非線性區(qū)域代表初始脈沖時(shí)，試樣中的氣流處于非穩(wěn)流狀態(tài)，隨著脈沖進(jìn)行，試樣內(nèi)的氣流逐漸變成穩(wěn)流狀態(tài)，即圖2(b)中的線性區(qū)域。隨著滲透率增大，線性擬合曲線的斜率絕對(duì)值變大，表明所需平衡時(shí)間減小。此外，線性區(qū)域出現(xiàn)的時(shí)間隨著滲透率降低而滯后，表明在其他因素相同的情況下，平衡壓力在低滲透率的巖石中達(dá)到穩(wěn)定梯度下降所需時(shí)間較長。不同滲透率下模擬所用各參數(shù)的值見表2，可以看出滲透率同時(shí)影響c與λu，λd，并與這3個(gè)數(shù)值呈負(fù)相關(guān)關(guān)系。

圖2 滲透率對(duì)脈沖壓力響應(yīng)的影響曲線Fig.2 Effect curves of permeability on pulse pressure response

表2 不同滲透率下模型中各參數(shù)的值Tab.2 Values of various parameters of the model under different permeability

2.2 孔隙度對(duì)脈沖壓力響應(yīng)影響

由式(7)可知，當(dāng)滲透率保持不變時(shí)，隨著孔隙度增大，c也逐漸增大，數(shù)值模擬中選擇0.5%，2%，5% 3組孔隙度進(jìn)行計(jì)算，模擬所用參數(shù)如表3所示。圖3(a)顯示了不同孔隙度下的脈沖壓力響應(yīng)，φ=0.5%時(shí)，脈沖壓力下降速度最慢，但其平衡時(shí)間明顯小于其他兩組數(shù)據(jù)。隨著孔隙度增大，脈沖壓力在前期下降速度逐漸增大，當(dāng)上下游儲(chǔ)層脈沖壓力接近時(shí)，孔隙度較大的巖樣脈沖壓力平衡所需時(shí)間明顯變長。

表3 不同孔隙度下模型中各參數(shù)的值Tab.3 Values of various parameters in the model under different porosity

一般認(rèn)為孔隙度對(duì)滲透率有一定影響，為了簡化分析，本文僅考慮單變量造成的影響，即不考慮各因素間的相互影響。當(dāng)滲透率保持不變時(shí)，孔隙度增加勢(shì)必會(huì)引起其他改變。分析認(rèn)為，巖樣內(nèi)的孔隙主要分為貫通的孔隙和封閉的孔隙。結(jié)合圖3(b)可知，孔隙度較大時(shí)，氣體滲透巖心過程達(dá)到穩(wěn)流狀態(tài)所需的時(shí)間更長，表明相同滲透率下，孔隙度大的巖心其內(nèi)部孔隙構(gòu)造更復(fù)雜，上下游壓差較小時(shí)，平衡所需時(shí)間更長?？紫抖容^大使得氣體在最初滲透過程中脈沖壓力變化更為明顯。

圖3 孔隙度對(duì)脈沖壓力響應(yīng)的影響曲線Fig.3 Effect of porosity on pulse pressure response

2.3 儲(chǔ)層體積對(duì)脈沖壓力響應(yīng)影響

相對(duì)于不同的上下游儲(chǔ)層體積，相同的儲(chǔ)層體積不僅可以簡化計(jì)算方法，還能夠減少實(shí)驗(yàn)時(shí)間[23]。常規(guī)的壓力脈沖裝置上下游儲(chǔ)層體積相同，因此，本部分主要研究上下游儲(chǔ)層體積相同時(shí)，不同的儲(chǔ)層體積對(duì)脈沖壓力響應(yīng)的影響。模擬選用實(shí)驗(yàn)室常用儲(chǔ)層體積5,50，200 cm3進(jìn)行計(jì)算，其他參數(shù)如表4所示。圖4(a)為儲(chǔ)層體積對(duì)脈沖壓力響應(yīng)的影響，當(dāng)上下游儲(chǔ)層體積為5 cm3時(shí)，脈沖壓力下降速度最快，且平衡所需時(shí)間最短，隨著儲(chǔ)層體積增大，壓力變化幅度逐漸減少，最終平衡時(shí)間也隨著儲(chǔ)層體積增加而增加，表明較小的儲(chǔ)層體積可以縮短平衡時(shí)間。儲(chǔ)層體積過小也會(huì)存在一些弊端，如圖4(b)所示，儲(chǔ)層體積越小，達(dá)到穩(wěn)流狀態(tài)所需時(shí)間越長，并且較小的儲(chǔ)層體積對(duì)于試樣的非均質(zhì)性更加敏感，當(dāng)采用小的儲(chǔ)層體積時(shí)，脈沖壓力變化引起的溫度變化以及外界擾動(dòng)等因素對(duì)測(cè)量結(jié)果造成的影響將更加顯著。此外，當(dāng)孔隙體積保持不變時(shí)，儲(chǔ)層體積較小時(shí)相同的孔隙體積誤差引起的滲透率計(jì)算誤差更大[24]，因此，在實(shí)驗(yàn)室進(jìn)行致密巖石滲透率實(shí)驗(yàn)時(shí)應(yīng)選擇合適的儲(chǔ)層體積進(jìn)行測(cè)量。

表4 不同儲(chǔ)層體積下模型中各參數(shù)的值Tab.4 Values of parameters in the model under different reservoir volumes

圖4 儲(chǔ)層體積對(duì)脈沖壓力響應(yīng)的影響曲線Fig.4 Effect curves of reservoir volume on pulse pressure response

3 結(jié) 論

(1)提出了采用有限差分法計(jì)算可壓縮流體通過多孔可壓縮介質(zhì)的一維瞬態(tài)流動(dòng)非線性偏微分方程，并在MATLAB環(huán)境下研究了滲透率、孔隙度、儲(chǔ)層體積對(duì)脈沖壓力響應(yīng)的影響。

(2)隨著滲透率增加，壓力減小速度變快，所需平衡時(shí)間變短。在脈沖過程中，氣體達(dá)到穩(wěn)流狀態(tài)的時(shí)間隨著滲透率增加而減少。

(3)脈沖壓力在孔隙度較大的試樣中前期變化更大，上下游儲(chǔ)層壓力差較小時(shí)，平衡時(shí)間與孔隙度成正比，表明滲透率、橫截面面積、上下游儲(chǔ)層體積保持不變時(shí)，孔隙度越大，其內(nèi)部孔隙越復(fù)雜，脈沖壓力主要經(jīng)歷快速減少階段和緩慢減少階段。

(4)上下游儲(chǔ)層體積相同時(shí)，隨著儲(chǔ)層體積增加，平衡所需時(shí)間變長，最終平衡壓力增大，氣流達(dá)到穩(wěn)定狀態(tài)時(shí)間變短。