ABAQUS中不同斷層帶的煤層地應(yīng)力場反演

孔 鵬，任廣聰，薄海江，馬新仿，劉羽欣，高海濱

(1.中聯(lián)煤層氣有限責(zé)任公司 煤層氣研發(fā)中心，山西 太原 030000； 2.中國石油大學(xué)(北京) 石油工程學(xué)院，北京 102249； 3.中海油田服務(wù)股份有限公司 油田生產(chǎn)事業(yè)部，天津 300000； 4.中聯(lián)煤層氣有限責(zé)任公司 晉城分公司，山西 晉城 048000)

引 言

近年來，我國煤層氣開發(fā)取得了一定的成果，柿莊區(qū)塊是我國煤層氣勘探程度較高的煤層氣區(qū)塊之一，由于煤層自身性質(zhì)的限制，水力壓裂改造技術(shù)是煤層氣高效開發(fā)的必須手段[1-2]。水力壓裂形成的裂縫沿最大水平主應(yīng)力擴(kuò)展[3-4]，因此，全面認(rèn)識煤層初始地應(yīng)力場分布對于預(yù)測水力裂縫走向、壓裂設(shè)計等具有重要意義[5]。通過Kaiser聲發(fā)射實驗[6]和古地磁實驗[7-8]可以準(zhǔn)確地獲取地應(yīng)力的大小和方向，但是室內(nèi)實驗只能測試有限位置的數(shù)據(jù)，可以根據(jù)這些有限數(shù)據(jù)點(diǎn)進(jìn)行全區(qū)地應(yīng)力反演，得到整個區(qū)塊的地應(yīng)力大小和方向分布。地層中普遍有一定規(guī)模的斷層發(fā)育，斷層的存在會影響地應(yīng)力的連續(xù)分布狀態(tài)[9-10]，因此，模型中對于斷層處理方式的不同會對計算結(jié)果產(chǎn)生較大影響。

郭懷志[11]最早提出利用有限元數(shù)學(xué)模型回歸分析方法求解初始地應(yīng)力場，奠定了多元回歸計算地應(yīng)力場的基礎(chǔ)，后來發(fā)展出偏最小二乘法、位移反分析方法、神經(jīng)網(wǎng)絡(luò)法和灰色關(guān)聯(lián)法等。近年來，許多學(xué)者使用不同分析軟件進(jìn)行地應(yīng)力場研究，包括ABAQUS、ANSYS、FLAC3D等，回歸方法從多元線性回歸發(fā)展到支持向量回歸、多約束方法、p值法[12-14]。但是這些方法注重對算法的改進(jìn)，忽略了儲層自身的地質(zhì)狀況和巖石力學(xué)性質(zhì)等的影響。

使用有限元方法建立的地質(zhì)模型中，一般以斷層帶的形式處理斷層，即在斷層附近加密網(wǎng)格，賦予不同屬性，目前使用的斷層帶處理方法主要包括單一性質(zhì)和過渡斷層帶[12]。斷層本身是一種較大規(guī)模的天然裂隙[15-16]，因此，可以考慮通過損傷力學(xué)模型來反映斷層，利用ABAQUS中的粘結(jié)單元(cohesive)來模擬斷層，考慮巖石的力學(xué)強(qiáng)度等性質(zhì)，使結(jié)果更真實可靠。

本文使用大型有限元分析軟件ABAQUS，建立地質(zhì)模型，使用3種方法處理斷層，分別為單一性質(zhì)斷層帶、過渡斷層帶和粘結(jié)單元斷層帶，通過對3種方法結(jié)果的對比，分析各種方法的優(yōu)劣性。

1 地應(yīng)力反演求解原理

利用邊界載荷反分析方法求解地應(yīng)力場，如圖1所示，將工區(qū)兩邊鉸接固定，另外兩邊共分為m份，每次以相同大小的載荷進(jìn)行加載，則一共加載m個載荷，在目標(biāo)井位置，如果所有載荷的疊加效果與實際測得地應(yīng)力大小和方向相同，則達(dá)到反演目的，從而獲取整個工區(qū)的應(yīng)力分布場。

每個目標(biāo)井的應(yīng)力約束條件包括最大水平主應(yīng)力、最小水平主應(yīng)力和最大水平主應(yīng)力方向，但是在ABAQUS中無法直接將地應(yīng)力方向作為約束條件輸入，為方便計算，需要將主應(yīng)力大小和方向轉(zhuǎn)化為正應(yīng)力S11、S22、S12，且拉力為負(fù)，壓力為正[17-18]，具體關(guān)系如下：

(1)

(2)

(3)

其中：S11為X軸的正應(yīng)力，MPa；S22為Y軸的正應(yīng)力，MPa；S12為XY平面上的切應(yīng)力，MPa；σH為最大水平主應(yīng)力，MPa；σh為最小水平主應(yīng)力，MPa；θ為X軸與最大水平主應(yīng)力的夾角，(°)。

因此，假設(shè)共有n口目標(biāo)井，則加載m個載荷時，每個載荷都會對目標(biāo)井產(chǎn)生影響，m個載荷對1#井的應(yīng)力影響為

(4)

同理，對第n口目標(biāo)井有，

(5)

根據(jù)上述過程，m個載荷對n口目標(biāo)井加載后，產(chǎn)生系數(shù)矩陣

(6)

可以看出，當(dāng)m=3n時，該系數(shù)矩陣有唯一解，可以直接求解系數(shù)；當(dāng)m<3n時，該系數(shù)矩陣有多解，可以通過多元線性回歸分析，將常數(shù)項設(shè)為0，求取系數(shù)的大小。

2 模型建立

柿莊南區(qū)塊位于沁水盆地，地層整體為單斜構(gòu)造，局部發(fā)育有褶曲、斷層和陷落柱，其中，斷層以正斷層為主。穩(wěn)定發(fā)育的主要煤層為二疊系下統(tǒng)山西組的3號和石炭系上統(tǒng)太原組的15號2套煤層，研究工區(qū)為3號煤層。3號煤層埋深在600～750 m，呈南淺北深，中部埋深大于兩側(cè)的趨勢，局部大于900 m，平均厚度6 m。具有儲層壓力低、滲透率低、飽和度低、非均質(zhì)性強(qiáng)的特點(diǎn)。平均孔隙度為4.35%，滲透率介于(0.04～1.10)×10-3μm2，儲層壓力2.4～3.01 MPa，儲層溫度21.7～27.5 ℃。煤層(原煤)含氣量一般在13～20 m3/t，含氣量相對較高。

基于柿莊南區(qū)塊3號煤層物性、巖石力學(xué)數(shù)據(jù)及煤層氣壓裂井分布，在ABAQUS中建立工區(qū)，大小為4 km×4 km，工區(qū)外圍增加1 km寬度，消除邊界效應(yīng)的影響。選擇well 1、well 2、well 3、well 4、well 5、well 6共6口井作為約束井，井眼位置如圖2所示。工區(qū)中存在2個斷層，分別為斷層Ⅰ和斷層Ⅱ。其中，well 1與斷層Ⅰ相鄰，well 2與斷層Ⅱ距離較近。

約束井地應(yīng)力大小和方向的準(zhǔn)確性對反演結(jié)果至關(guān)重要。地應(yīng)力與泊松比、楊氏模量等有關(guān)，又因本區(qū)塊受構(gòu)造運(yùn)動的影響，水平主應(yīng)力很大部分來源于地質(zhì)構(gòu)造運(yùn)動產(chǎn)生的構(gòu)造應(yīng)力，因此，采用組合彈簧模型計算地應(yīng)力剖面較為合適[19]。水力裂縫沿最大水平主應(yīng)力方向擴(kuò)展，因此，地應(yīng)力方向根據(jù)鄰井水力壓裂微地震監(jiān)測確定。約束井主應(yīng)力大小、方向以及正應(yīng)力、切應(yīng)力如下表所示。

組合彈簧模型計算方法如下：

(7)

(8)

(9)

式中，σH為最大水平主應(yīng)力，MPa；σh為最小水平主應(yīng)力，MPa；β1為最大構(gòu)造應(yīng)力系數(shù)；β2為最小構(gòu)造應(yīng)力系數(shù)；μs為測井縱波時差，μs/m；pp為孔隙壓力，MPa；η為有效應(yīng)力系數(shù)；σv為垂向應(yīng)力，MPa；ρb為巖石密度，g/cm3；ρ0為處理研究井段頂界至井口的地層密度平均值，g/cm3；H0為處理井段頂界深度，m；H為處理井段底界深度，m。

約束井主應(yīng)力大小、方向及式(1)—式(3)計算的正應(yīng)力、切應(yīng)力如表1所示。

圖3(a)為單一性質(zhì)斷層帶，即在斷層位置加密網(wǎng)格，并且賦予相同屬性，整個模型中存在斷層帶和基質(zhì)巖石2種材料；圖3(b)為過渡斷層帶，即將斷層區(qū)域分為過渡區(qū)和中心區(qū)，過渡區(qū)和中心區(qū)賦予不同屬性，整個模型中存在過渡帶、斷層帶和基質(zhì)巖石3種材料；圖3(c)為粘結(jié)單元斷層帶，即將斷層區(qū)域分為過渡區(qū)和cohesive粘結(jié)單元區(qū)，利用粘結(jié)單元模擬斷層行為，中間線條表示cohesive粘結(jié)單元。圖3(a)和(b)是目前研究中使用的斷層帶處理方法，這兩種方法雖然將斷層和基質(zhì)的性質(zhì)區(qū)分開來，但是并沒有考慮斷層帶的力學(xué)行為，因此，建立圖3(c)的cohesive粘結(jié)單元斷層帶。

如圖4所示，將模型劃分網(wǎng)格，巖石基質(zhì)和斷層過渡區(qū)域網(wǎng)格類型均為滲流/應(yīng)力單元(CPE4P)，cohesive粘結(jié)單元網(wǎng)格類型為孔隙應(yīng)力粘結(jié)單元(COH2D4P)。使用過渡網(wǎng)格，工區(qū)采用小尺寸網(wǎng)格，斷層區(qū)域局部加密，提高計算精度；邊界區(qū)域采用大尺寸網(wǎng)格，減少計算時間。由于斷層帶處理方式不同，整體網(wǎng)格數(shù)量有微弱差異。

使用滲流/應(yīng)力單元可以考慮地層孔隙和地層流體的影響，粘結(jié)單元從力學(xué)損傷角度考慮更加符合斷層自身性質(zhì)。為了控制巖石力學(xué)之外的變量對結(jié)果的影響，孔隙度取5%，滲透率取0.1×10-3μm2。其中，巖石基質(zhì)楊氏模量、泊松比數(shù)據(jù)來源于室內(nèi)三軸壓縮實驗結(jié)果，過渡區(qū)域取值實驗結(jié)果的0.8倍，中心區(qū)域和粘結(jié)單元取值實驗結(jié)果的0.6倍。根據(jù)室內(nèi)三軸實驗結(jié)果，柿莊南煤巖楊氏模量為4.13～8.03 GPa，泊松比0.35～0.44，因此，模型基質(zhì)取值如表2所示。

表2 材料屬性Tab.2 Properties of rock

3 模型求解

(1)計算過程

以單一性質(zhì)斷層帶為例展示求解過程，如圖5所示，將工區(qū)兩邊鉸接固定，另外兩邊共劃分為8等份，每次施加20 MPa的壓力，并分別記錄8次施加應(yīng)力后目標(biāo)井位置的S11、S22、S12值(表3)。

表3 單次載荷計算結(jié)果Tab.3 Calculation results of loads

根據(jù)式6，使用多元線性回歸方法計算每個載荷對應(yīng)的系數(shù)，計算結(jié)果見表4。將校正后的8個載荷同時加載，計算結(jié)果即為初始地應(yīng)力場分布。

表4 線性回歸結(jié)果Tab.4 Linear regression results

同理，可以求得過渡斷層帶、粘結(jié)單元斷層帶2種處理方式下的載荷系數(shù)。

(2)計算結(jié)果

圖6分別為單一性質(zhì)斷層帶、過渡斷層帶和粘結(jié)單元斷層帶3種模型的最大水平主應(yīng)力計算結(jié)果?？梢钥闯?，最大水平主應(yīng)力大小在斷層帶區(qū)域及附近存在差異，尤其是斷層端部，應(yīng)力大小與附近區(qū)域存在突變，表現(xiàn)出非連續(xù)性的分布。斷層位置應(yīng)力值小于周圍基質(zhì)區(qū)域，斷層Ⅰ處最大水平主應(yīng)力一般小于10 MPa，而周圍基質(zhì)區(qū)域大于10 MPa。另外，過渡斷層帶模型中，雖然應(yīng)力分布在斷層帶附近出現(xiàn)突變，但是不會出現(xiàn)錯位，即應(yīng)力大小連續(xù)分布，在斷層處出現(xiàn)“間斷”；單一性質(zhì)斷層帶和粘結(jié)單元斷層帶則會出現(xiàn)分布錯位。

圖7分別為單一性質(zhì)斷層帶、過渡斷層帶和粘結(jié)單元斷層帶3種模型的最小水平主應(yīng)力計算結(jié)果?？梢钥闯觯?種模型最小水平主應(yīng)力大小總體一致，在斷層Ⅱ下端，應(yīng)力分布出現(xiàn)部分差異。斷層位置最小水平主應(yīng)力偏小，一般小于6 MPa。另外，過渡斷層帶模型中，應(yīng)力分布在斷層帶附近只出現(xiàn)突變，但是不會出現(xiàn)錯位。

圖8和圖9分別為單一性質(zhì)斷層帶、過渡斷層帶和粘結(jié)單元斷層帶3種模型的最大水平主應(yīng)力、最小水平主應(yīng)力方向的計算結(jié)果。

3種模型計算的應(yīng)力方向一致，最大水平主應(yīng)力方向與最小水平主應(yīng)力垂直，斷層處應(yīng)力方向發(fā)生突變，最小水平主應(yīng)力方向與斷層走向一致。

(3)結(jié)果分析

為了確定不同斷層帶對計算結(jié)果的影響，通過后處理分別讀取3種情況下斷層Ⅰ和斷層Ⅱ附近的well 1井和well 2井水平應(yīng)力大小和方向，并與實際測試值對比，結(jié)果如表5和表6所示。a表示單一斷層帶，b表示過渡斷層帶，c表示粘結(jié)單元斷層帶。

表5 水平應(yīng)力計算結(jié)果與實際測試值對比Tab.5 Comparison between calculating and testing horizontal stress

表6 水平應(yīng)力計算值與實際測試值之差Tab.6 Difference between calculating and testing horizontal stress

可以看出，3種斷層帶處理方法計算的最大水平主應(yīng)力和最小水平主應(yīng)力絕對值比較接近，與實際值均存在誤差。

其中，斷層Ⅰ附近的well 1井水平主應(yīng)力絕對值偏小，水平應(yīng)力差值也偏小。過渡斷層帶處理方式的水平應(yīng)力差計算結(jié)果最接近實際值，偏小0.12 MPa；單一斷層帶計算結(jié)果與實際值相差最大，偏小0.34 MPa；粘結(jié)單元斷層帶計算結(jié)果次之，偏小0.25 MPa。

斷層Ⅱ附近的well 2井水平主應(yīng)力絕對值偏大，水平應(yīng)力差值也偏大。其中，單一斷層帶的水平應(yīng)力差計算結(jié)果最接近實際值，偏大0.35 MPa；過渡斷層帶與實際值偏差最大，偏大0.49 MPa；粘結(jié)單元斷層帶計算結(jié)果次之，偏大0.38 MPa。

從3種斷層帶處理方法計算出的地應(yīng)力方向可以看出，粘結(jié)單元斷層帶結(jié)果最接近實際方向。Well 1井計算結(jié)果為粘結(jié)單元斷層帶最準(zhǔn)確，單一斷層帶次之，過渡斷層帶方向偏差最大；well 2井計算結(jié)果為粘結(jié)單元最準(zhǔn)確，單一斷層帶次之，過渡斷層帶方向偏差最大。因此，可以判斷，粘結(jié)單元計算的方向最準(zhǔn)確。

4 結(jié) 論

(1)在斷層帶尤其是斷層端部，應(yīng)力大小與附近區(qū)域存在突變，表現(xiàn)非連續(xù)性分布，地應(yīng)力方向在斷層帶發(fā)生突變。

(2)計算地應(yīng)力大小時，粘結(jié)單元斷層帶計算結(jié)果的誤差均為中等，單一斷層帶和過渡斷層帶計算結(jié)果不穩(wěn)定；計算地應(yīng)力方向時，粘結(jié)單元計算的水平應(yīng)力方向最準(zhǔn)確，單一斷層帶次之，過渡斷層帶計算的水平應(yīng)力方向偏差最大。

(3)同時考慮地應(yīng)力大小和方向誤差，粘結(jié)單元斷層帶處理方法計算方向最準(zhǔn)確，其原因為粘結(jié)單元將斷層視為天然裂隙處理，考慮損傷演化。模擬計算煤層地應(yīng)力場分布時，應(yīng)關(guān)注斷層等地質(zhì)構(gòu)造的巖石力學(xué)性質(zhì)，并將之與巖石基質(zhì)性質(zhì)區(qū)分開來，使結(jié)果更符合實際情況。