基于多場耦合效應(yīng)的水平井壓裂應(yīng)力場分析

官 兵，李士斌，張立剛，陳雙慶，王業(yè)強.

(1.東北石油大學(xué)石油工程學(xué)院，黑龍江大慶 163318；2.大慶油田有限責(zé)任公司第三采油廠，黑龍江大慶 163000)

?

基于多場耦合效應(yīng)的水平井壓裂應(yīng)力場分析

官 兵1，李士斌1，張立剛1，陳雙慶1，王業(yè)強2.

(1.東北石油大學(xué)石油工程學(xué)院，黑龍江大慶 163318；2.大慶油田有限責(zé)任公司第三采油廠，黑龍江大慶 163000)

為了有效分析水平井壓裂過程中的應(yīng)力場分布情況，考慮應(yīng)力-滲流-溫度多物理場因素的耦合作用，利用有限元數(shù)值模擬手段，開展了對水平井壓裂裂縫近井應(yīng)力場的模擬分析。首先通過巖石物理力學(xué)參數(shù)測定試驗，得到了有限元分析的基本數(shù)據(jù)；之后引入應(yīng)力-滲流-溫度耦合理論，建立了水平井?dāng)_動應(yīng)力場的有限元分析模型，模擬計算得到研究區(qū)的壓力、溫度和應(yīng)力分布情況；最后，綜合多場耦合效應(yīng)得到近水平井筒裂縫周圍的應(yīng)力場變化規(guī)律，確定裂縫在平行井筒方向?qū)ζ渲車鷳?yīng)力場的擾動范圍主要集中在1倍的裂縫長度以內(nèi)，裂縫產(chǎn)生的誘導(dǎo)應(yīng)力降低了原井筒附近的壓應(yīng)力大小。將數(shù)值解結(jié)果與已有擾動應(yīng)力數(shù)學(xué)模型計算的相應(yīng)解析解進行對比，驗證了有限元模型求解的合理性，為水平井壓裂應(yīng)力場的分析研究提供了一定的基礎(chǔ)和參考。

水平井壓裂；多場耦合；應(yīng)力場；有限元；合理性

水平井壓裂應(yīng)力場是影響裂縫起裂和擴展行為的關(guān)鍵因素之一[1-2]，而壓裂過程中應(yīng)力場又受控于多個物理場的共同作用。水平井壓裂改造過程中，井筒進入儲層，打破原始地層的力學(xué)平衡，在原地應(yīng)力場的約束和井筒內(nèi)液柱壓力的共同作用下，井筒周圍產(chǎn)生應(yīng)力集中效應(yīng)，壓裂裂縫產(chǎn)生誘導(dǎo)應(yīng)力，二者疊加使地層巖石發(fā)生變形，從而擾動原始地應(yīng)力場。同時由于注入壓力與地層孔隙壓力的差異造成壓裂液滲濾引起附加應(yīng)力，加之注入液體溫度與儲層溫度的差異也會產(chǎn)生熱應(yīng)力作用。因而可知，水平井壓裂是一個應(yīng)力-滲流-溫度多場耦合的過程。

早在1995年，陳勉[3]等人就基于多孔彈性理論、疊加原理，建立了考慮內(nèi)壓、地應(yīng)力及壓裂液滲流效應(yīng)的大斜度井筒周圍應(yīng)力分布計算模型；2004年，徐嚴(yán)波[4]在此基礎(chǔ)上引入地層溫度產(chǎn)生的熱應(yīng)力，建立了水平井筒周圍的應(yīng)力分布模型；之后劉洪[5]、尹建[6]、翁定為[7]等都對人工裂縫的誘導(dǎo)應(yīng)力場進行了研究，推導(dǎo)出了它的計算數(shù)學(xué)模型，其中尹建[8]、翁定為等人還利用有限元數(shù)值模擬方法分析了壓裂裂縫周圍的應(yīng)力場變化?；谶@些研究，已經(jīng)建立了比較成熟的水平井井筒及壓裂裂縫所產(chǎn)生的誘導(dǎo)應(yīng)力場計算模型，并且提出了多場耦合作用的影響，但是對于同時考慮流-固-熱耦合理論和水平井與裂縫共同形成的擾動應(yīng)力場研究卻較少。本文在前人研究的基礎(chǔ)上，針對考慮應(yīng)力-滲流-溫度多場作用各種因素影響下的水平井壓裂擾動應(yīng)力場計算模型開展了進一步的研究，提出有限元數(shù)值分析方法，模擬計算了壓裂過程中多物理場耦合作用下的應(yīng)力場變化，并將其與解析解對比進行了合理性驗證。

1 試驗數(shù)據(jù)

選取大慶油田外圍的扶楊油層致密砂巖組作為本文的研究對象，遵循中華人民共和國地質(zhì)礦產(chǎn)行業(yè)標(biāo)準(zhǔn)(DZ/T0276.4—2015)《巖石物理力學(xué)性質(zhì)試驗規(guī)程》，挑選完整且表面沒有明顯裂紋的巖石制作成直徑×高度規(guī)格為25 mm×50 mm的圓柱體巖心若干，記錄好取樣深度和對應(yīng)標(biāo)號。通過開展孔隙度、滲透率測試及三軸應(yīng)力試驗等獲得了研究區(qū)儲層巖石的基本物性參數(shù)和力學(xué)參數(shù)(表1)，為儲層裂縫數(shù)值模擬分析提供了基礎(chǔ)。

表1 取樣巖心試驗數(shù)據(jù)結(jié)果

2 多場耦合理論

多場耦合作用是指一個系統(tǒng)中所存在的物理場相互影響，區(qū)別于單一的物理作用，不同物理場中每一個都在其他物理場中扮演著不同的角色，始終存在耦合關(guān)系，包括其起始狀態(tài)和整個過程。分析一個復(fù)雜多物理場的前提是引入耦合關(guān)系?；趪鴥?nèi)外大量學(xué)者[9-11]對不同物理場耦合關(guān)系數(shù)學(xué)模型的分析，推導(dǎo)出了一組考慮多物理場完全耦合而又具有可求解性的方程組。其中，假設(shè)以拉應(yīng)力為正，孔隙水壓力為負，得到應(yīng)力-滲流-溫度耦合模型的完整方程組如下：

(1)

式中G——巖石剪切模量，Pa；ν——地層泊松比；p——地層孔隙水壓力，Pa；T——地層溫度，℃；T0——無應(yīng)力狀態(tài)下的絕對溫度，℃；t——時間，s；K′——排水體積模量，且K′=2G(1+ν)/3(1-2ν)，Pa；

Ks——巖石骨架的有效體積模量，Pa；

αp——地層Biot系數(shù)，αp=1-K′/Ks，αp≤1，其值取決于地層巖石壓縮性；

αT——巖體熱膨脹系數(shù)，℃-1；

al——流體熱膨脹系數(shù)，℃-1；

as——巖石骨架的體積膨脹系數(shù)，℃-1；

Fi——i方向的體積力，N；

ui(i=x,y,z)——i方向的位移，m；

εv——體積應(yīng)變，εv=εxx+εyy+εzz；

Φ——孔隙度，%；

κ——巖石固有滲透系數(shù)，m2；

Vl——多孔介質(zhì)中流體速度矢量，s-1；

βl——流體的體積模量，Pa；

ρl——流體密度，kg·m-3；

γ——γ=(2μ+3λ)β，μ和λ為拉梅常數(shù)，J/m3，β為各向同性固體的線性熱膨脹系數(shù)，℃-1；

g——重力加速度，m·s-2；

z——垂直坐標(biāo)，m；

(ρcp)t——充滿了流體的多孔介質(zhì)的比熱容，J/(m3·℃)；

Kt——充滿了流體的多孔介質(zhì)的熱傳導(dǎo)系數(shù)，W/(m3·℃)；

qt——充滿了流體的多孔介質(zhì)的熱源強度，W/m3。

由于應(yīng)力-滲流-溫度多場耦合方程組既考慮了不同物理場的相互作用又是時間的物理過程，有典型的非線性，除少數(shù)特殊情況可求出解析解外，一般需要用數(shù)值方法進行求解?；趹?yīng)力-滲流-溫度多場耦合理論模型的建立，本文選用專門針對偏微分方程組求解的有限元分析軟件COMSOLMultiphysics對其進行求解。

3 有限元模型求解

3.1 模型建立

利用有限元數(shù)值分析方法分析水平井壓裂裂縫擾動應(yīng)力場，首先需要建立擾動應(yīng)力場的三維幾何模型。設(shè)定地層模型的長寬高尺寸分別為2000m×1000m×1000m；井筒半徑為5m，水平井筒長1000m，假設(shè)水平井筒沿最小水平主應(yīng)力方向延伸，裂縫與水平井筒夾角為90°；壓裂人工裂縫形態(tài)設(shè)定為橢球型，其半長為100m，半高為50m，開度為5m。繪制其幾何模型如圖1所示。其中，直角坐標(biāo)系(x,y,z)與地層巖石所受的3個主應(yīng)力(σh,σH,σv)方向相對應(yīng)。幾何模型建立以后，需要將其轉(zhuǎn)化為有限元數(shù)學(xué)模型。有限元模擬數(shù)學(xué)模型的建立首先要按照有限元數(shù)值分析所要求的數(shù)學(xué)、力學(xué)規(guī)則進行單元劃分。本模型采用自由四面體單元對模型進行剖分，為了提高計算精度，在剖分井筒和裂縫單元時需要進行網(wǎng)格加密處理，最后得到研究區(qū)塊的擾動應(yīng)力場有限元模型如圖2所示。

圖1 擾動應(yīng)力場的三維幾何模型示意圖Fig.1 3D geometry model of disturbed stress field

通過COMSOL有限元數(shù)值分析軟件，在建立有限元模型之后，只要確定儲層巖石的物理力學(xué)參數(shù)，并給予模型合理的施加載荷和邊界約束，就能夠?qū)λ骄畨毫褦_動應(yīng)力場的分布進行模擬計算。

圖2 擾動應(yīng)力場的有限元模型示意圖Fig.2 Finite element model of disturbed stress field

3.2 基本數(shù)據(jù)

假設(shè)地層巖石為多孔彈性介質(zhì)，其中的流體流動滿足Darcy定律，井筒與裂縫中充滿單相液流，井筒與裂縫完全連通，井壁與裂縫面所受壓力和溫度等同于其中液體的壓力和溫度，注入液體比地層溫度低24℃，其他給定分析水平井壓裂擾動應(yīng)力場分布的模型計算基本數(shù)據(jù)見表2。

3.3 求解結(jié)果

根據(jù)多場耦合理論模型及地層巖石的物理力學(xué)參數(shù)，對構(gòu)建的水平井壓裂裂縫有限元分析模型進行計算，得到了水平井壓裂裂縫的應(yīng)力-滲流-溫度多場耦合的模擬計算結(jié)果。

表2 模型計算基本數(shù)據(jù)

圖3是水平井壓裂過程中的流場壓力作用云圖。如圖3a所示，地層整體的壓力分布在20 MPa左右，井筒和裂縫內(nèi)的流體壓力較高為30 MPa，而地層多孔介質(zhì)具有滲透作用，井筒和裂縫周圍會發(fā)生一定的液體滲濾作用(其流動方向如圖中流線方向所示)，因而也會產(chǎn)生壓力傳遞，從而出現(xiàn)壓力過渡階段(其作用方向如圖中紅色箭頭方向，呈發(fā)散狀)。截取水平井筒軸線所在的二維水平面，分析水平井筒和壓裂裂縫周圍的局部壓力變化，如圖3b所示。井壁和裂縫面處的壓力最高為30 MPa，而隨著沿箭頭方向距離的增加，壓力逐漸降低，直至接近原地層孔隙壓力20 MPa。

圖3 流場壓力計算云圖及等值線圖Fig.3 Flow field pressure result and contour map

圖4是水平井壓裂過程中的溫度變化云圖。如圖4a所示，地層整體溫度隨著深度的增加而逐漸升高，其分布范圍為20～59.903℃。由于井筒內(nèi)注入初始液體溫度較低，假設(shè)井筒和裂縫壁面處的溫度與地面溫度一致為20℃。同樣截取二維水平面，如圖4b所示。水平井筒所在深度的原地層溫度分布在44 MPa左右，地層溫度高于井筒內(nèi)的液體溫度，從而產(chǎn)生傳熱過程，地層溫度向井筒方向傳遞，使井筒和裂縫周圍的溫度逐漸升高。

圖4 溫度場的計算云圖及等值線圖Fig.4 Temperature field result and contour map

水平井筒壓裂的流場壓力和溫度場變化都會影響其周圍的應(yīng)力分布，同時考慮原始地層的應(yīng)力場作用，可以得到多場耦合作用下的水平井壓裂過程中的地層水平主應(yīng)力的分布變化情況，如圖5和圖6所示。其中壓應(yīng)力為負，拉應(yīng)力為正。

圖5 最小水平主應(yīng)力計算云圖及等值線圖Fig.5 Minimum horizontal principal stress result and contour map

如圖5a所示，地層整體的最小水平主應(yīng)力方向沒有發(fā)生明顯偏轉(zhuǎn)，依舊平行于x軸方向；地層整體的最小水平主應(yīng)力大小集中在30 MPa左右，與原地應(yīng)力場設(shè)置一致，只在裂縫周圍發(fā)生了顯著變化。截取水平井筒軸線所在的二維水平面，選取壓裂裂縫周圍的局部應(yīng)力場進行分析，如圖5b所示。從圖中可以看出水平井筒和裂縫的出現(xiàn)使原地應(yīng)力場的應(yīng)力發(fā)生了擾動，裂縫周圍應(yīng)力場發(fā)生了一定的偏轉(zhuǎn)，其擾動范圍為-37.536～6.5907 MPa；而只有在裂縫壁面附近出現(xiàn)了正的拉應(yīng)力，其應(yīng)力值為6.6 MPa左右；隨著距離裂縫壁面水平距離x軸方向的增加，其周圍應(yīng)力由拉應(yīng)力變?yōu)閴簯?yīng)力而逐漸增大，裂縫尖端出現(xiàn)較大的應(yīng)力集中，為-37 MPa左右。

圖6 最大水平主應(yīng)力計算云圖及等值線圖Fig.6 Maximum horizontal principal stress result and contour map

如圖6a所示，地層整體的最大水平主應(yīng)力方向沒有發(fā)生明顯偏轉(zhuǎn)，依舊平行于y軸方向；地層整體的最大水平主應(yīng)力大小集中在50 MPa左右，與原地應(yīng)力場設(shè)置一致，與最小水平主應(yīng)力不同的是，水平井井筒和裂縫周圍都發(fā)生了較為明顯的變化。同樣截取水平井筒軸線所在的二維水平面，分析壓裂裂縫周圍的局部應(yīng)力場，如圖6b所示。裂縫周圍應(yīng)力場發(fā)生了一定的偏轉(zhuǎn)，而井筒和裂縫同時作用，使最大水平主應(yīng)力發(fā)生擾動，其擾動范圍為-97.062～-29.388 MPa；裂縫壁面附近出現(xiàn)了正的拉應(yīng)力，與原地應(yīng)力場的負壓應(yīng)力疊加，其值為-29 MPa左右；隨著距離裂縫壁面x軸方向水平距離的增加，其周圍壓應(yīng)力逐漸增大，裂縫尖端出現(xiàn)較大的應(yīng)力集中，為-97 MPa左右，井筒周圍的應(yīng)力變化主要分布在y軸方向，其變化規(guī)律同裂縫。

通過以上對水平井壓裂裂縫應(yīng)力場、滲流場和溫度場的分析可以看出，水平井壓裂裂縫會使井壁附近的原始地層條件發(fā)生擾動，它是一個多場耦合作用的過程，因而在分析其壓裂應(yīng)力場擾動規(guī)律的同時，必須考慮其他物理場的作用效果。

3.4 結(jié)果分析

有限元數(shù)值分析方法能夠快速高效地對多場耦合方程組進行求解，但是針對水平井分段壓裂擾動應(yīng)力場的計算分析，有限元分析是否正確合理還需要展開進一步研究。以大慶外圍扶楊油層某水平井分段壓裂的水平井作為研究對象，提取水平井筒井壁所在直線(x軸方向)上的計算應(yīng)力值，繪制出擾動應(yīng)力變化規(guī)律的曲線。通過將有限元數(shù)值模擬計算的數(shù)值結(jié)果與已有的成熟擾動應(yīng)力數(shù)學(xué)模型[5,12-13]計算的解析結(jié)果進行對比(圖7)，驗證有限元模型求解的合理性。

圖7 井壁處擾動應(yīng)力數(shù)值解與解析解對比Fig.7 Comparison between the numerical and analytical solutions of disturbance stress in well wall

從圖7中可以看出，針對水平井壓裂擾動應(yīng)力場分析這一問題，利用有限元數(shù)值模擬軟件得到的結(jié)果曲線在局部范圍內(nèi)會出現(xiàn)一定的波動，這是由于在進行三維建模求解過程中網(wǎng)格的精度劃分所造成的。而受限于數(shù)值模擬軟件和計算機硬件配置，由網(wǎng)格精度所引起的計算誤差無法避免，但不影響其對應(yīng)力大小隨距離變化規(guī)律的分析。在平行井筒方向，裂縫對其周圍的應(yīng)力場擾動范圍主要分布在1倍的裂縫長度以內(nèi)，擾動后的最小水平主應(yīng)力下降速率先緩后快，當(dāng)其與裂縫面x軸方向的距離接近1倍裂縫長度后，裂縫產(chǎn)生的擾動逐漸消失；而擾動后的最大主應(yīng)力先快速降至最低，在其0.8倍裂縫半長處出現(xiàn)最低值，之后又緩慢上升，直至其影響程度趨于0。

通過與已有數(shù)學(xué)模型所計算的解析解結(jié)果進行對比，兩種方法所得到的應(yīng)力變化規(guī)律基本一致，其結(jié)果基本相符。因此，本文利用有限元數(shù)值分析方法展開對水平井壓裂擾動應(yīng)力分析是合理的。

4 結(jié)論

(1)通過對研究區(qū)致密砂巖進行的密度、孔隙度、滲透率基本物性測試，以及三軸應(yīng)力試驗，掌握了該儲層巖石的物性特征：孔隙度平均值為4.70%；滲透率平均值小于0.5 mD，分布在納達西級別；巖石楊氏模量平均值為20387.1 MPa，泊松比平均值為0.130，抗壓強度平均值為213.54 MPa。

(2)基于應(yīng)力-滲流-溫度多場耦合方程，利用有限元數(shù)值分析軟件COMSOL Multiphysics建立了水平井分段壓裂擾動應(yīng)力場的有限元分析模型，提出了求解多物理場的有限元數(shù)值方法，并計算了裂縫周圍的壓力、溫度和應(yīng)力分布情況。

(3)針對水平井分段壓裂問題，選取相同研究對象，將有限元分析方法模擬計算得到的數(shù)值解與已有水平井壓裂擾動應(yīng)力數(shù)學(xué)模型計算的解析解進行對比，得到的求解結(jié)果基本一致，驗證了有限元數(shù)值法研究水平井壓裂擾動應(yīng)力場的合理性。

[1] MANCHANDA R, SHARMA M M. Impact of completion design on fracture complexity in horizontal wells[C]. SPE159899, 2012.

[2] ROUSSEL N P, MANCHANDA R, SHARMA M M. Implications of fracturing pressure data recorded during a horizontal completion on stage spacing design[C]. SPE152631, 2012.

[3] 陳勉,陳治喜,黃榮樽.大斜度井水壓裂縫起裂研究[J].石油大學(xué)學(xué)報(自然科學(xué)版),1995,19(2):30-35.

[4] 徐嚴(yán)波.水平井水力壓裂基礎(chǔ)理論研究[D].成都:西南石油大學(xué),2004:1-128.

[5] 劉洪,胡永全,趙金洲,等.重復(fù)壓裂氣井誘導(dǎo)應(yīng)力場模擬研究[J].巖石力學(xué)與工程學(xué)報,2004,23(23):4022-4027.

[6] 尹建,郭建春,曾凡輝.水平井分段壓裂射孔間距優(yōu)化方法[J].石油鉆探技術(shù),2012,40(5):67-71.

[7] 翁定為,嚴(yán)星明,盧擁軍,等.考慮應(yīng)力干擾的致密油氣壓裂設(shè)計及實現(xiàn)方法[J].中國礦業(yè)大學(xué)學(xué)報,2014,43(4):639-645.

[8] 尹建.水平井分段壓裂誘導(dǎo)應(yīng)力場研究與應(yīng)用[D].成都:西南石油大學(xué),2014:1-75..

[9] 王自明,江同文,宋文杰,等.彈性氣藏中氣體滲流的氣-固-熱耦合數(shù)學(xué)模型[J].天然氣工業(yè),2003,23(3):85-87.

[10] 賈善坡,鄒臣頌,王越之,等.基于熱-流-固耦合模型的石油鉆井施工過程數(shù)值分析[J].巖土力學(xué),2012,33(S2):321-328.

[11] 曹偉,康志勤,呂義清,等.基于COMSOL Multiphysics的煤層氣對流熱采分析[J].地下空間與工程學(xué)報,2014,10(5):1139-1145.

[12] 邵尚奇,田守嶒,李根生,等.水平井縫網(wǎng)壓裂裂縫間距的優(yōu)化[J].石油鉆探技術(shù),2014,42(1):86-90.

[13] 曾順鵬,張國強,韓家新,等.多裂縫應(yīng)力陰影效應(yīng)模型及水平井分段壓裂優(yōu)化設(shè)計[J].天然氣工業(yè),2015,35(3):55-59.

Analysis of Stress Field of Horizontal Well FracturingBased on Multi-Field Coupling Effect

Guan Bing1, Li Shibin1, Zhang Ligang1, Chen Shuangqing1, Wang Yeqiang2

(1.School of Petroleum Engineering, Northeast Petroleum University, Daqing, Heilongjiang 163318, China;2.No.3 Oil Production Plant of Daqing Oilfield Co., Ltd., Daqing, Heilongjiang 163000, China)

In order to analyze the distribution of stress field in the process of horizontal wells fracturing effectively, the multi-field coupling effect was considered and the finite element numerical simulation analysis was applied. Firstly, through the determination test of rock physical and mechanical parameters, the basic data were obtained. Then introducing the theory of stress - seepage - temperature coupling, the finite element analysis model of the disturbance stress field in horizontal well fracturing was established. The results of pressure, temperature and stress distribution in the study area were simulated and calculated. Finally, combined with the multi-field coupling effect, the change law of the stress field around the fractures was achieved. It was found that in the parallel wellbore direction the stress disturbance range of one fracture was mainly concentrated within 1 time of its length and its induced stress made the original compressive stress around the wellbore reduced. Comparing the numerical solution results with the corresponding analytic solutions of the existing disturbed stress mathematic model, the rationality of the finite element model was verified. It can provide a certain basis and reference to study the stress field in horizontal wells fracturing.

horizontal well fracturing; multi-field coupling; stress field; finite element; rationality

東北石油大學(xué)研究生創(chuàng)新科研項目“致密油氣藏體積壓裂局部應(yīng)力場改造控制機理研究”(YJSCX2015-017NEPU)、國家自然科學(xué)基金項目“致密儲層井工廠壓裂應(yīng)力場演化機理及調(diào)控方法”(51504067)聯(lián)合資助。

官兵(1990—)，女，黑龍江省雙城市人，2013年畢業(yè)于東北石油大學(xué)石油工程專業(yè)，在讀博士研究生，主要從事非常規(guī)壓裂裂縫應(yīng)力場數(shù)值模擬分析的研究工作。郵箱：sygcgb1990@163.com.

TE357

A