鄧金根，陳崢嶸，，耿亞楠，劉書杰，朱海燕

(1.中國石油大學(xué)油氣資源與探測國家重點(diǎn)實(shí)驗(yàn)室，北京102249;2.中海油研究總院，北京100027)

目前，常規(guī)能源的開發(fā)和利用已逐漸不能滿足現(xiàn)代工業(yè)的發(fā)展需要，勘探開發(fā)的重點(diǎn)投向了難以開采的頁巖氣。頁巖儲層是一種低孔隙度和特低滲透率的致密儲層，其生產(chǎn)完全依賴于壓裂效果。對頁巖進(jìn)行壓裂施工必須對頁巖起裂的規(guī)律有足夠的認(rèn)識，其中地應(yīng)力模型是指導(dǎo)頁巖壓裂施工的極為重要的理論基礎(chǔ)。目前地應(yīng)力預(yù)測模型大部分建立在地層均質(zhì)各向同性的基礎(chǔ)上，而頁巖儲層具有較強(qiáng)的非均質(zhì)性和各向異性。Thiercelin［1］討論了線彈性的橫觀各向同性地層地應(yīng)力的計(jì)算模型，Higgins［2］在研究 Baxter頁巖地層的地應(yīng)力時分別應(yīng)用了各向同性和橫觀各向同性的地應(yīng)力計(jì)算模型，并對二者進(jìn)行了計(jì)算和比較。目前，多數(shù)學(xué)者的理論都只考慮上覆地層壓力和構(gòu)造應(yīng)力的影響，對地應(yīng)力的描述并不全面，因此并不能完全描述橫觀各向同性地層地應(yīng)力的特性。筆者在前人理論研究的基礎(chǔ)上，結(jié)合滲流和孔隙壓力的影響建立考慮頁巖儲層橫觀各向同性的地應(yīng)力模型。

1 頁巖本構(gòu)特性

頁巖由于其沉積歷史的特殊性［3］通常具有良好的分層結(jié)構(gòu)，在構(gòu)造和組成上呈連續(xù)變化，因此巖石的性質(zhì)在平面和垂向上會有區(qū)別，頁巖巖石呈現(xiàn)出橫觀各向同性的特性［4-5］，其巖石力學(xué)性質(zhì)也表現(xiàn)出橫觀各向同性［6］，即頁巖儲層水平方向是各向同性的而垂直方向則表現(xiàn)為各向異性(圖1)。因此，對頁巖儲層可通過橫觀各向同性模型進(jìn)行研究。

圖1 頁巖橫觀各向同性單元體Fig.1 Element of shale representing transversely isotropy body

對于橫觀各向同性的頁巖，應(yīng)力和應(yīng)變之間的本構(gòu)方程可通過修正后的胡克定律［7］表示為

式中，σij為二階應(yīng)力張量的分量形式，MPa;εkl為二階應(yīng)變張量的分量形式;Cijkl為四階彈性剛度張量的分量形式，MPa;i，j，k，l=1，2，3。

二階應(yīng)力張量和二階應(yīng)變張量實(shí)際分別都只有6個獨(dú)立的分量，因此四階彈性剛度張量可表示為6×6矩陣。根據(jù)地層的對稱性，剛度張量的獨(dú)立參數(shù)將減少，因此剛度張量C可表示為

對于橫觀各向同性地層，矩陣各參數(shù)有以下關(guān)系 C11=C22，C12=C21，C13=C31=C23=C32，C44=，其中 C66可以表示為因此可以用5個獨(dú)立的變量表示矩陣C［9-10］。要完全描述矩陣C需要求得以上5個獨(dú)立的參數(shù)，該5個參數(shù)［11］可以表示為

式中，ρ為體積密度，kg/m3;vPH和vPV分別為水平方向和垂直方向的P波波速，m/s;vSH和vSV分別為水平方向和垂直方向的S波波速，m/s;v45為與地層對稱軸成45°的P波(m=－1)或S波(m=1)波速，m/s。

經(jīng)過聲波測井的測量，可求得橫觀各向同性的5個獨(dú)立參數(shù)。則垂直和水平方向上的彈性模量和泊松比［12-13］分別為

由彈性模量和泊松比表示的頁巖本構(gòu)模型為

其中

式中，σ和τ分別為作用在單元體上的正應(yīng)力和剪應(yīng)力，MPa;ε和γ分別為單元體上的正應(yīng)變和剪應(yīng)變;Eh和Ev分別為水平方向和垂直方向上的彈性模量，MPa;νh和νv分別為施加水平和垂直應(yīng)變時水平應(yīng)變的泊松比;Gv為垂直平面的剪切模量，MPa。

地層力學(xué)參數(shù)的測量通常通過應(yīng)力應(yīng)變的關(guān)系或聲波速度測量得到，通過應(yīng)力應(yīng)變關(guān)系施加載荷測試，一般測定的是標(biāo)準(zhǔn)的力學(xué)彈性參數(shù)(如彈性模量、剪切模量和泊松比等)，聲波速度測試則是測量聲波的強(qiáng)度參數(shù)(如C11、C33和C12等)，而聲波實(shí)際上也是施加在巖石上的載荷，聲學(xué)性質(zhì)與聲波傳播過程中巖石的應(yīng)力應(yīng)變有關(guān)。通過上述分析，這兩種力學(xué)參數(shù)的測量方法等價，通過以上兩種方法都可以解決巖石力學(xué)參數(shù)的測量問題，在應(yīng)用過程中可根據(jù)實(shí)際情況擇優(yōu)使用。

2 頁巖地應(yīng)力預(yù)測模型的建立

地應(yīng)力場中垂直地應(yīng)力由上覆地層的重力產(chǎn)生，可通過地層的密度和厚度的積分計(jì)算得到。水平地應(yīng)力的形成由地層構(gòu)造運(yùn)動產(chǎn)生的應(yīng)變、不同地層之間地層性質(zhì)的各向異性和非均質(zhì)性及其他地質(zhì)因素共同作用產(chǎn)生。為了預(yù)測水平地應(yīng)力并構(gòu)建理想的地應(yīng)力模型，首先必須對所研究地層的物理力學(xué)參數(shù)的各向異性情況有足夠的了解。建立頁巖儲層地應(yīng)力模型必須考慮頁巖的橫觀各向同性的特性。

Thiercelin討論了線彈性各向同性和橫觀各向同性的地層地應(yīng)力計(jì)算模型，對于各向同性地層，在平面應(yīng)變的假設(shè)下，考慮上覆地層壓力和構(gòu)造應(yīng)力的影響，水平地應(yīng)力模型［1］可表示為

式中，σv為垂直地應(yīng)力，MPa;σH和σh分別為水平最大和最小地應(yīng)力，MPa;pp為地層孔隙壓力，MPa;α為有效應(yīng)力系數(shù)，通常取1;εH和εh分別為水平最大和最小構(gòu)造應(yīng)變。

對于橫觀各向同性的地層，Thiercelin在各向同性地應(yīng)力模型的基礎(chǔ)上考慮了地層垂直和水平方向上力學(xué)性質(zhì)的差別，建立了橫觀各向同性地層地應(yīng)力模型，即

式中，ξ為滲流系數(shù)，通常取0。

其他地質(zhì)條件也可能產(chǎn)生水平地應(yīng)力，如孔隙壓力的變化。結(jié)合上述因素，基于線彈性假設(shè)建立考慮地層橫觀各向同性、孔隙壓力變化、構(gòu)造載荷和滲流的有效水平地應(yīng)力模型，即

式中，σh，g，σh，p和 σh，t分別為上覆地層壓力、孔隙壓力變化和構(gòu)造應(yīng)力產(chǎn)生的水平地應(yīng)力，MPa。

圖2所示為頁巖儲層中影響水平地應(yīng)力的各種因素。

圖2 頁巖儲層中影響水平地應(yīng)力的各種因素Fig.2 Factors of influencing horizontal in-situ stress in shale formation

如果忽略構(gòu)造應(yīng)力的影響，所產(chǎn)生的水平地應(yīng)力將相同(σH=σh)，如果某一方向上存在構(gòu)造作用產(chǎn)生的變形，則將導(dǎo)致水平地應(yīng)力在不同方向上的應(yīng)力差(σH＞σh)。相應(yīng)的地應(yīng)力為

其中

孔隙壓力的變化通過系數(shù)(1+K0)改變水平地應(yīng)力。地質(zhì)上的構(gòu)造應(yīng)變則通過壓縮和拉伸使地應(yīng)力相應(yīng)的增加和減少，二者之間的比例關(guān)系可用地層的各向異性參數(shù)K1和K2表示。對于橫觀各向同性地層，如果不存在構(gòu)造應(yīng)力產(chǎn)生的變形，則兩個水平地應(yīng)力相等，而一旦考慮了構(gòu)造應(yīng)力，兩個水平地應(yīng)力將產(chǎn)生差別。除以上因素外，還有一些影響地應(yīng)力預(yù)測的地質(zhì)力學(xué)因素。該模型相對簡單且所需的數(shù)據(jù)可通過測井解釋方法獲取，便于實(shí)際應(yīng)用。

3 算 例

圖3 Eh/Ev對最大和最小水平地應(yīng)力的影響Fig.3 Effect of Eh/Evon the maximum and minimum horizontal in-situ stress

根據(jù)上述理論分析，編寫計(jì)算機(jī)程序模擬Eh/Ev、νv/νh、孔隙壓力變化對地應(yīng)力的影響。模擬參數(shù):地層平均密度為2.61 g/cm3，有效應(yīng)力系數(shù)為1，滲流系數(shù)為0，εH和 εh分別為0.3和0.1。

3.1 Eh/Ev對頁巖儲層地應(yīng)力的影響

模擬參數(shù):垂直和水平方向的泊松比均為0.25，孔隙壓力變化為0.5 MPa。圖3為Eh/Ev對最大和最小水平地應(yīng)力的影響。隨著Eh/Ev的增大，地層的各向異性更強(qiáng)，水平地應(yīng)力將整體變大。當(dāng)Eh/Ev=4時，最大和最小水平地應(yīng)力分別是各向同性(Eh/Ev=1)情況下的2.30和2.57倍，平均增加144%，因此忽略橫觀各向同性對地層地應(yīng)力的影響，預(yù)測結(jié)果會產(chǎn)生較大誤差，各向同性的預(yù)測結(jié)果明顯偏低。

3.2 νv/νh對頁巖儲層地應(yīng)力的影響

模擬參數(shù):垂直和水平方向的彈性模量均為16 GPa，孔隙壓力變化為0.5 MPa。圖4 為 νv/νh對最大和最小水平地應(yīng)力的影響。

圖4 vv/vh對最大和最小水平地應(yīng)力的影響Fig.4 Effect of vv/vhon the maximum and minimum horizontal in-situ stress

從圖4中可以看出，最大和最小水平地應(yīng)力隨νv/νh值增大而增大。其中，νv/νh=2.2 的最大和最小水平地應(yīng)力分別是νv/νh=1時的1.15和1.18倍，平均增加16.5%。因此垂直和水平方向上的泊松比對地應(yīng)力的預(yù)測亦會產(chǎn)生較大的影響。

3.3 孔隙壓力變化對頁巖儲層地應(yīng)力的影響

圖5 孔隙壓力變化對最大和最小水平地應(yīng)力的影響Fig.5 Effect of pore pressure on the maximum and minimum horizontal in-situ stress

模擬參數(shù):垂直和水平方向的彈性模量均為16 GPa，垂直和水平方向的泊松比均為0.2。圖5為孔隙壓力變化對最大和最小水平地應(yīng)力的影響。預(yù)測地應(yīng)力考慮孔隙壓力的變化時，最大和最小水平地應(yīng)力隨著孔隙應(yīng)力變化的增大而增大。其中，孔隙壓力變化值為2.4 MPa時的最大和最小水平地應(yīng)力分別是孔隙壓力不變化時的1.13和1.15倍，平均增加13.8%。因此，孔隙壓力的變化對地應(yīng)力的預(yù)測也是一個不可忽視的因素。

3.4 地層升沉對頁巖儲層地應(yīng)力的影響

一些學(xué)者在建立地應(yīng)力預(yù)測模型時還考慮了地層的升沉作用(地層上升或者沉降)對地應(yīng)力預(yù)測的影響［14-15］。該模型表達(dá)式為

對應(yīng)的最大和最小水平地應(yīng)力分別為

式中，σh，u為升沉作用影響產(chǎn)生的水平地應(yīng)力，MPa;h為當(dāng)前地層深度，m，Δh為當(dāng)前地層深度與上升或者沉降之前地層深度的差值，即升沉高度，m。

模擬參數(shù):垂直和水平方向的彈性模量均為16 GPa，垂直和水平方向的泊松比均為0.2，孔隙壓力變化為0.5 MPa，其他參數(shù)同上。圖6為升沉高度對最大和最小水平地應(yīng)力的影響。最大和最小水平地應(yīng)力隨著升沉高度的增大而增大，但升沉作用對地應(yīng)力的預(yù)測影響不大。

圖6 升沉高度對最大和最小水平地應(yīng)力的影響Fig.6 Effect of subsidence or uplift on the maximum and minimum horizontal in-situ stress

4 模型驗(yàn)證

選取Baxter頁巖的地層參數(shù)進(jìn)行模擬計(jì)算說明式(2)的準(zhǔn)確性，將計(jì)算結(jié)果與地層瞬間閉合壓力(ISIP)進(jìn)行比較(圖7)，可以看出本文模型和Thiercelin模型的結(jié)果都與ISIP基本一致，本文模型結(jié)果誤差最大僅為1.2%，所得地應(yīng)力符合Baxter頁巖地應(yīng)力的計(jì)算情況［2］，驗(yàn)證了所建立的地應(yīng)力模型的準(zhǔn)確性。相同條件下，各向同性模型計(jì)算的地應(yīng)力結(jié)果與實(shí)際誤差較大，并不準(zhǔn)確，因此采用橫觀各向同性地應(yīng)力模型更符合現(xiàn)場實(shí)際情況。

圖7 各向同性模型、Thiercelin模型和本文模型計(jì)算的最小水平地應(yīng)力與ISIP數(shù)據(jù)的比較Fig.7 Comparison of ISIP and the minimum horizontal in-situ stress computed by isotropy，Thiercelin and transversely isotropy models

5 結(jié)論

(1)水平與垂直彈性模量之比、垂直與水平泊松比的比值和孔隙壓力的變化對地應(yīng)力的預(yù)測影響較大，且以上3種因素增大，地應(yīng)力的預(yù)測結(jié)果也將增大。升沉作用對地應(yīng)力的影響較小。

(2)對于頁巖儲層，橫觀各向同性的地應(yīng)力模型要比各向同性的地應(yīng)力模型更能表示地層地應(yīng)力的真實(shí)情況。準(zhǔn)確計(jì)算地應(yīng)力的分布對頁巖儲層現(xiàn)場的完井壓裂增產(chǎn)設(shè)計(jì)有重要的指導(dǎo)作用。

［1］THIERCELIN M J，PLUMB R A.A core-based prediction of lithologic stress contrasts in east Texas formations［J］.SPE Formation Evaluation，1994，9(4):251-258.

［2］HIGGINS S，GOODWIN S，DONALD A，et al.Anisotropic stress models improve completion design in the Baxter shale［R］.SPE 115736，2008.

［3］RIVERA S R，HANDWERGER D，KIESCHNICK J，et al.Accounting for heterogeneity provides a new perspective for completions in tight gas shales［R］.American Rock Mechanics Association 05-758，2005.

［4］LINGS M L，PENNIGNTON D S，DASH D F T.Anisotropic stiffness parameters and their measurement in a stiff natural clay［J］.Geotechnique，2000，50(2):109-125.

［5］GRAHAM J，CROOKS J H A，BELL A L.Time effects on the stress-strain behavior of natural soft clays ［J］.Geotechnique，1983，33(3):327-340.

［6］SAYERS C M.Seismic anisotropy of shales［J］.Geophysical Prospecting，2005，53(5):667-676.

［7］ GAUTAM R.Anisotropy in deformations and hydraulic properties of Colorado shale［D］.Calgary:Civil Engineering，University of Calgary，2004.

［8］NYE J F.Physical properties of crystals［M］.Oxford:Oxford University Press，1985.

［9］NORRIS A N，SINHA B K.Weak elastic anisotropy and the tube wave［J］.Geophysics，1993，58(8):1091-1098.

［10］WALSH J，SINHA B，DONALD A.Formation anisotropy parameters using borehole sonic data［R］.Society of Petrophysicists and Well-Log Analysis 2006-TT，2006.

［11］HORNBY B E，HOWIE J M，INCE D W.Anisotropy correction for deviated-well sonic logs:application to seismic well tie［J］.Geophysics，2003，68(2):464-471.

［12］SAROUT J，MOLEZ L，GUEGUEN Y，et al.Shale dynamic properties and anisotropy under triaxial loading:experimental and theoretical investigations［J］.Physical and Chemistry of the Earth，2007，32(8/14):896-906.

［13］RIVERA S R，DEENDAYALU C，YANG Y K.Unlocking the unconventional oil and gas reservoirs:the effect of laminated heterogeneity in wellbore stability and completion of light gas shale reservois［R］.Offshore Technology Conference 20269，2009.

［14］SAFDAR K，SAJJAD A，HAN H X，et al.Importance of shale anisotropy in estimating in-situ stresses and wellbore stability analysis in Horn river basin［R］.SPE 149433，2011.

［15］PRIOUL R，KARPFINGER F，DEENADAYALU C，et al.Improving fracture initiation predictions on arbitrarily oriented wells in anisotropic shales［R］.SPE 147462，2011.

［16］RIVERA S R，DEENADAYALU C，CHERTOV M，et al.Improving horizontal completions on heterogeneous tight shales［R］.SPE 146998，2011.