三維現今地應力場有限元數值計算

劉洪濤，謝 兵，楊向同，巴 旦，周鵬遙，黎莉莉

(1.塔里木油田油氣工程研究院，新疆維吾爾自治區 庫爾勒 841000；2.百通思達(北京)石油科技發展有限公司，北京 100107)

三維現今地應力場有限元數值計算

劉洪濤1，謝 兵2，楊向同1，巴 旦1，周鵬遙1，黎莉莉1

(1.塔里木油田油氣工程研究院，新疆維吾爾自治區 庫爾勒 841000；2.百通思達(北京)石油科技發展有限公司，北京 100107)

研究現今地應力場的分布規律對油氣藏的勘探、開發有重要的指導意義。依據Petrel地質模型網格節點數據，建立有限元三維精細幾何模型，并對巖體和斷層進行不同尺寸的網格剖分。綜合利用地球物理測井與巖芯資料，采用VB.net編程求得三維巖石力學參數，根據單元中心坐標分別賦予相應的參數。運用大型有限元軟件ANSYS建立了塔里木盆地##2地區三維有限元模型，以測井應力解釋剖面為井點應力約束條件，經過多次反演試算得到模型的合理邊界條件。與現場資料數據比較發現：數值模擬結果吻合性較好。可為三維現今地應力場的模擬提供數值計算方法。

有限元；三維精細模型；巖石力學參數；地應力場

0 引言

地應力場的分布規律可為油田開發中的井網布置、水力壓裂和優化開采提供重要依據。因此，根據現場勘探開發獲取的資料，進行地應力數值模擬具有重要工程價值。由于三維地質建模復雜煩瑣和三維巖石參數的求取難度大等原因，目前地應力場的模擬計算基本上從以下兩方面進行：1)忽略地層幾何形狀不規則性和斷層的影響，進行簡單的二維模擬分析[1-2]；2)以計算參考井的巖石力學參數代替一定區塊的巖石力學參數[3-5]來進行三維模擬。地下儲層通常是復雜的地質體，地層的幾何形狀和斷層等細部構造與地層的應力場分布密切相關，第1種研究方法通常不能準確分析構造應力場和斷層附近的局部應力場；第2種研究方法忽略了地層的不均質性，不能研究透鏡體、砂泥巖互層等復雜地質體的應力場。

本文依據Petrel地質模型數據，建立了與油藏目的層地質條件相一致的三維精細有限元模型。并綜合利用聲波測井、地震速度場和室內巖芯試驗等數據資料，采用VB.net編程計算了整個研究區的三維巖石力學參數。通過反演試算法修改模型邊界條件，采用ANSYS計算了儲層三維現今地應力場。

1 幾何建模與網格剖分

1.1 塔里木油田##2區塊地質概況

##2區塊位于新疆維吾爾自治區阿克蘇地區境內，構造上位于塔里木盆地庫車坳陷北部克拉蘇構造帶上，西鄰大北氣田，東接克拉2氣田，南部為拜城凹陷，含氣層系主要為下白堊統巴什基奇克組。##2構造為繼承性山前沖斷帶，定型于第四紀喜馬拉雅運動晚期，早前屬于前展式疊瓦斷序列，晚期—現今屬于后展式斷展褶皺序列，以地層褶皺為主。垂向構造上由巖石強度的差異劃分為三套力學結構層，自上而下依次為：膏巖層、碎屑巖層、剛性層。上覆膏巖層通過發生塑性變形與剪切變形吸收掉大部分擠壓能量，從而造成層間的強烈剪切變形作用。

如圖1，##1-2區塊發育5條較大的區域性斷層，延伸距離均在30 km以上，最長可達55 km，切穿整個區塊，將區塊分成6個主要的斷塊。區塊北部斷塊多發育背斜構造，形成大面積的油氣圈閉，次級斷層發育較少；南部斷塊構造相對平緩，僅有少數小面積油氣圈閉零星分布，次級斷層較為發育。區塊中部發育各種類型的斷層相關褶皺，整體反映出南北方向和東西方向的擠壓應力場。

1.2 Petrel地質模型與ANSYS幾何模型互連

為了消除邊界效應對計算模型的影響，本次所選取的計算巖體區域略大于實際模擬的##2區塊。東西方向長約35 km，南北方向長約6.6 km，平均厚度約300 m，深度方向可分為3個結構層。

圖1 ##2區塊關鍵井區地質構造(1∶100 000)

根據地質力學理論，儲層的幾何形態控制著整個研究區的現今構造應力，細部構造(褶皺和斷層等)影響著研究區的局部應力場，因此，在采用有限元分析建模時，必須考慮這兩個影響因素。由于儲層結構構造的復雜性，目前很難用有限元軟件(如ANSYS)直接進行幾何建模，而三維地質建模軟件Petrel可以建立復雜的地質模型，進行Petrel地質模型和ANSYS幾何模型的互連，可以建立接近油氣藏實際地質特征的全三維精心有限元幾何模型。采用Petrel對地質模型劃分網格，并提取各個層面的節點坐標，借用Solidworks對層面節點進行貝塞爾曲面擬合，建立各個儲層層面的幾何模型并保存為ANSYS可以識別的iges格式文件。在ANSYS中對所有層面進行實體生成從而完成幾何建模(如圖3)。提取Petrel模型中斷層三維邊界線導入ANSYS，生成斷層實體，從而保證斷層的形態(傾向、傾角和斷距)與地質構造的一致性。

圖2 Petrel地質模型(1∶400 0000)

圖3 ANSYS幾何模型

1.3 網格劃分

斷層幾何尺寸數量級一般在幾m或幾十m，而儲層巖體區域數量級在幾km。斷層的規模、幾何形態及力學性質對周圍地應力場的分布有著重要的影響[6]，在有限元分析中不能忽視。為了提高計算效率并保證斷層附近區域的計算精度，采用巖體區域和斷層獨立剖分網格。如圖4～5，巖體區域采用較粗網格(一般100 m)，而斷層采用精細網格(一般10 m)。鑒于Solid186單元是ANSYS中的一個高階三維20節點固體結構單元，具有二次位移模式，可以較好地模擬不規則的網格模型，所以巖體區域和斷層均采用Solid186單元。

圖4 巖體區域粗網格模型(1∶400 0000)

圖5 斷層細網格模型

2 三維巖石力學參數的計算

目前求取巖石力學參數的方法主要有兩種：一是在實驗室內對巖芯樣品進行實測；二是根據地球物理測井資料推求巖石力學參數，并根據室內試驗結果進行校正。室內三軸試驗有較高的可參考性，然而由于巖芯樣品的獲取有限，不能得到整個目的層的物性參數。聲波測井具有測量的連續性和較高分辨率的優點，可以精確地確定探測范圍內的物性參數和地層參數，但只能得出井口附近的巖體力學參數，不能獲得整個區域尤其是斷層的力學參數。地震速度場與聲波測井相比，準確性較差并存在概率誤差，但所獲的速度場是三維數據體，能夠計算整個區域的力學參數。

2.1 動態參數的求取

根據彈性理論，動態彈性參數與巖石縱橫波波速及密度存在式(1)的關系[7]。其中巖石縱波波速可通過聲波測井資料得到，橫波波速可通過縱橫波測井或經驗公式、統計關系求得。

(1)

式中：Vp為縱波波速，m/s；Vs為橫波波速，m/s；ρ為巖石密度，g/cm3；μd為動態泊松比；Ed為動態彈性模量，MPa。

2.2 靜態參數的求取

巖石力學特性參數的靜態值和動態值存在著一定的差值，動態力學參數不能直接用于有限元模擬分析中。Montmayeur[8]和Yale[9]指出靜態力學參數更能代表石油工程中巖石的受載條件。林英松[10]等在三軸應力下對砂、泥巖等巖芯進行巖石力學參數的動靜態同步測試，并進行動靜態參數的線形回歸，結果表明：巖石的動靜態楊氏模量之間有較好的相關性，而動靜態泊松比之間的關系不明顯。根據某地區室內三軸試驗求得的巖芯樣品靜態參數，進行動靜態參數的數據回歸，得出動靜態巖石力學參數的統計關系，如式(2)。

(2)

式中：μs為靜態泊松比；μd為動態泊松比；Es為靜態彈性模量，MPa；Ed為動態彈性模量，MPa。

2.3 編程計算

地震速度場數據體包含了整個研究工區的三維波速信息，在此基礎上利用動靜態參數計算公式，可以建立整個研究工區的三維巖石力學參數。為了便于參數調整與多次重復計算，進行了VB.net相關程序的開發。程序的計算流程如圖6所示。

對于每個有限元網格單元，根據單元中心坐標，在巖石力學參數數據體中循環查找距離最近的參數，找到后賦予本單元。依次循環，給每個單元賦予相應的巖石參數。

圖6 巖石力學參數程序計算流程圖

3 計算與結果分析

3.1 邊界施加

邊界條件主要包括模型的邊界應力方向和大小，應力方向主要是依據該區域地質構造演化和生產井的鉆井誘導縫和井壁坍塌方位確定，如圖7所示，主應力主方向為NNE方向。應力大小主要依據測井資料、室內試驗和水力壓裂等資料確定。由于邊界應力場分布的復雜性，初始邊界條件難以準確設定，因而，應力場的數值模擬實際上是一個正演與反演的多次反復計算修正的過程。

本次計算中，初始邊界條件設定為水平最大主應力180 MPa，水平最小主應力為150 MPa，垂向主應力為160～170 MPa。

圖7 井點裂縫電子成像圖

定義A是實測井點應力場矩陣，B為計算井點應力場矩陣，則定義誤差Δ：

Δ=‖A-B‖2。

(3)

以生產井井點上實測的應力場為約束條件，校核ANSYS計算的井點應力場，用線性迭加原理來反復修改邊界條件，直到誤差Δ滿足需求，停止循環。為方便結果的校正與邊界條件的修改，本文開發了用戶VB.net校核程序。程序的計算流程如圖8所示。

3.2 結果分析

提取出最終計算結果中生產井井點最大水平主應力，與測井實測值進行對比分析，表1展示了##2-2-14等6口井的對比結果，應力差值均控制在15 MPa以內(誤差在10%以內)，滿足工程需要。

據據##2地質構造演化和生產井的鉆井誘導縫和井壁坍塌方位，確定了##205、##2-2-14、##2-2-12等20口井井點在巴一底最大水平主應力方位。如圖9所示，數值模擬結果與實測結果相比：大部分井水平最大主應力方向能夠較好地吻合，##205、##201、##203、##2-1-7 4口井的方向相差較大。根據地質模型判斷這4口井沒位于斷層上，主應力方向不可能發生較大程度的跳動，從而判斷這4口井的主應力方向實測解釋可能存在一定誤差。

表1 井點水平最大主應力模擬值與實測值比較

圖8 邊界條件校正程序流程圖

圖9 20口井井點水平最大主應力方向模擬值(箭頭表示)與實測值(橢圓表示)對比(1∶100 000)

##2區塊的最大水平主應力云圖如圖10所示，最大水平主應力方向如圖11所示。模型中部地區是背斜構造，最大地應力主要是南北方向的擠壓力；從中間區域向西部，應力方向由原來的近南北向逐漸逆時針旋轉成近北東—南西向；從中間區域向東部，應力方向則由原來的近南北向逐漸順時針旋轉成近北西—南東向。斷層與周圍地層的地應力相比，由于斷層的彈性模量較小、泊松比較大，斷裂(斷層)帶應力顯著較小，成為應力釋放區。

4 結語

通過計算實例的分析，可以得出以下結論：

1)利用Petrel建立地質模型并采集網格節點數據的方法，可以建立全三維精細ANSYS幾何模型。

2)對巖體區域和斷層進行獨立劃分網格，不但提高了軟件計算效率，而且可以保證斷層局部區域的計算精度。

3)綜合利用地球物理測井資料，采用VB.net計算可以有效地獲取三維巖石力學參數；由于地震速度場數值的精度較差，對速度場數據進行更好的校正是提高巖石力學參數的計算精確的關鍵。

4)本文提供的計算地應力場的方法能夠更準確地模擬地應力場，可為油田開發中的井網布置、水力壓裂和優化開采提供一定的參考。

圖10 ##2區塊水平最大主應力大小云圖(1∶200 000)

圖11 ##2區塊水平最大主應力方向(1∶200 000)

[1] 唐治，潘一山，閻海鵬，等.基于ANSYS的二維地應力場的分析[J].科學技術與工程，2010，10(28): 6926-6929.

[2] 劉顯太，戴俊生，徐建春，等.純41斷塊沙四段現今地應力場有限元模擬[J].石油勘探與開發，2003，30(3):126-128.

[3] 馮立，張學婧.海拉爾盆地三維地應力場數值模擬[J].大慶油田地質與開發，2011，30(2):115-119.

[4] 張廣明，熊春明，劉合，等.復雜斷塊地應力數值模擬方法研究[J].斷塊油田，2011，18(6):710-713.

[5] 曾聯波，肖淑蓉，羅安湘.陜甘寧盆地中部靖安地區現今應力場三維有限元數值模擬及其在油田開發中的意義[J].地質力學學報，1998，4(3):58-62.

[6] 沈海超，程遠方，王京印，等.斷層對地應力場影響的有限元研究[J].大慶石油地質與開發，2007，2(26):34-37.

[7] 路保平，鮑洪志.巖石力學參數求取方法進展[J].石油鉆探技術，2005，33(5):44-47.

[8] Montmayeur H，Graves R M.Prediction of static elastic /mechanical properties of consolidated and unconsolidated sands from acoustic measurements: correlatious[C]// SPE Annual Technical Conference and Exhibition.Louisiana:Society of Petroleum Engineers，1986:1-16.

[9] Yale D P，Jamieson W H.Static and dynamic rock mechanical properties in the Hugo ton and Panoma field，Kansas[C]// SPE Mid-Continent Gas Symposium.Texas: Society of Petroleum Engineers，1994:209-219.

[10] 林英松，葛洪奎，王順昌.巖石動靜力學參數的試驗研究[J].巖石力學與工程學報，1998，17(2):216-222.

The FEA Numerical Algorithm for 3D Modern Crustal Stress Field

LIU Hong-tao，et al.

(ExplorationandDevelopmentResearchInstituteofTarimOilfield，KorlaXinjiangUygur841000，China)

The study on modern Crustal stress field distribution is with guiding significance to prospect and exploit oil and gas reservoirs.3D refined finite element model is established in terms of Petrel geology model grid point data，rock and faults are meshed in different sizes grid.Comprehensively utilizing geophysical well logging information and Core data，3D rock mechanical parameters are assigned to unit according to center coordinates，are programmed by VB.NET program.3D modern of Tarim basin ##2 is established by large ANSYS software to test that the well stress interpretation profile is well point stress constraint condition，and boundary conditions are got after several inversion test.The result is full of much accuracy by comprising with spot materials data.This paper can provide a numerical algorithm for 3D modern crustal stress field.

finite element；3D fine model；rock mechanical parameters；crustal stress field

10.3969/j.issn.1009-8984.2016.04.019

2016-08-01

劉洪濤(1983-)，男(漢)，黑龍江肇州，工程師 主要研究試油完井工程。

TE32

A

1009-8984(2016)04-0071-05