基于“地質-工程”耦合平臺的順北斷溶體儲層地應力分布規律分析

趙兵 , 羅攀登, 張妹珠, 周志敏, 黃欣

(1.中國石油化工股份有限公司西北油田分公司, 烏魯木齊 830011; 2.武漢大學科學技術發展研究院, 武漢 430072;3.武漢大學土木建筑工程學院, 武漢 430072)

隨著油田勘探和開發的不斷加深和擴展,碳酸鹽巖儲層逐漸進入了人們的視野,該種類型油藏具有儲量大、分布廣等優點,已成為了人們開采和研究的重點對象[1-2]。塔里木盆地的順北油田是新開發的儲量較大的碳酸鹽巖斷溶體儲層之一,儲層中的斷裂破碎區域如天然洞穴、天然裂縫內儲藏有充足的油氣資源[3]。針對斷溶體儲層的特點(孔隙度低、滲透率低、非均質性強等)提出水力壓裂開采技術,該技術通過增加儲層中裂縫和溶洞之間的有效溝通路徑,降低分隔性,從而實現油田增產[4]。研究表明,水力壓裂過程中的水力裂縫的擴展軌跡與地應力具有很強的相關性[5]。順北油田屬于碳酸鹽巖斷控裂縫-洞穴型儲層,該油田主要位于交匯斷裂帶區域,地質條件復雜。為了利用水力裂縫溝通盡可能多的縫洞儲集體以實現增產,則必須掌握儲層三維地應力場空間分布規律。因此,如何準確獲得順北油田三維地應力場的精細空間分布規律,成為利用水力壓裂技術進行油藏開采增產的關鍵問題。

在地應力測量和分析方面,現場地應力測量方法根據測量原理可分為[6]：一是破裂法,如水壓致裂法、聲發射法和井壁崩落法;二是變形法,如孔徑變形法和孔壁應變法等;三是地球物理法,包括微地震波法和地面電位法等;四是構造法,如宏觀構造反推法、顯微構造分析法。在實際情形中,地質體中的縫洞儲集體以及破裂帶會對地應力場分布情況產生較大影響,現場實測的地應力值存在一定的離散型,僅依靠現場測量無法全面準確地反映出地質體地應力場的空間分布。因此,目前許多學者采用有限元軟件對地層進行模擬計算,以精細刻畫地應力場空間分布。王愛國等[7]利用子模型法提取目標區域邊界條件,在構造地應力場的尺度對黑山峽大柳樹壩址區進行數值模擬分析計算,描述了地震斷層處地應力場特性,對壩體破裂危險性進行了分析評價。白克新[8]針對含斷層地層巷道覆巖應力分布分析問題,根據彈性力學理論,利用數值模擬軟件建立了平面模型。他對存在斷層及正常條件下地層的應力特征及位移規律進行了比較分析,并利用物理相似模擬方法獲取開采工作面接近斷層過程中的巖層破裂、垮落特征以及地層應力演化規律,提出了回采巷道的合理支護方案。李萬才等[9]以某壩址岸坡作為研究對象,建立三維數值模型,分析了在河谷下切影響下地質體應力、位移場特征以及塑性區的演化過程。然而,對于地質結構復雜的油氣藏儲層來說,由于有限元等軟件的建模能力較弱,難以建立精細的三維儲層模型。

在石油工程中,常利用地質建模軟件來構建精細的三維儲層模型,輔助人員進行開采策略調整。在地質建模方面,Petrel軟件作為一款的地質建模軟件被廣泛應用于油氣開采行業。該軟件采用角點網格構建三維地質模型,能夠精細地刻畫出地質體中的斷層、侵入、褶皺、尖滅等地質構造情形,同時還整合地質體的地震監測數據、油藏屬性等儲層信息。在實際應用中,胡向陽等[10]為解決碳酸鹽巖縫洞型油藏三維地質建模的難題,提出了多元約束三維地質建模方法。該方法能夠表征了縫洞儲集體的空間分布,并為該區的增儲建產提供有效保障。孫先達等[11]為了實現營城組升平地區火山巖儲層的三維可視化動態表述和展示,通過Petrel軟件利用構造層面及斷層數據建立了構造模型和斷層模型。在此基礎上,他們采用確定性建模和隨機建模結合,構建了該區火山巖相三維地質模型,實現在三維空間上刻畫了典型火山巖體的巖相特征,并研究其在三度空間上的變化規律。

為了滿足地質-工程一體化的需求,目前許多學者利用地質建模軟件完成了數值仿真模型的構建。郭穎星等[12]為解決包含斷層的復雜地質構造建模問題,開發了Petrel與Adina等軟件網格數據之間的接口程序。該程序可以實現地質模型網格數據與仿真軟件網格數據的相互移植,從而實現有限元模型的精確建模。劉鈺洋等[13]開發了一系列算法,并將其集成到Petrel2ANSYS中,實現了Petrel中角點網格的3D屬性模型與ANSYS中有限元網格3D模型之間的雙向轉換。徐珂等[14]構建了高尚堡油田深層油藏南區三維巖石力學場,通過Petrel地質模型與ANSYS有限元模型聯合建模,對該地區三維地應力場空間分布規律開展了研究,并為油藏開發過程中壓裂優勢區和優勢井段的優選提供了理論依據。Park等[15]開發了PET2OGS一套3種算法來集成靜態模型(Petrel)和動態模型(OpenGeoSys),轉換有限差分法(FDM)網格為有限元法網格。

目前對油藏儲層的應力場研究主要集中于油藏區域的整體應力場分布,沒有針對縫洞型儲集體油藏地應力場的研究,缺乏對斷溶體儲層區域及內部、不同斷裂帶交匯區三維地應力場分布規律的分析。AiFrac是適用于水力壓裂、油氣開采等領域的多物理場耦合作用巖體裂縫擴展軟件[16]。現通過研發“地質-工程”耦合平臺,將地質模型轉化為AiFrac仿真軟件使用的模型,完成地質模型和數值仿真一體化,以實現三維地應力場和水力壓裂的精細化模擬,形成實時驅動的儲層地質和壓裂過程數字孿生平臺。研究成果將為油氣開采提供可靠的決策支持和優化方案。

1 “地質-工程”耦合平臺

雖然Petrel等地質建模軟件能夠構建三維地質模型,對儲層進行精細刻畫并描述儲層的空間分布和屬性分布特征,但在計算地應力場方面效率較低,不能夠勝任大型三維地應力場的模擬。相比之下,有限元等數值仿真軟件相較于地質建模軟件則具備很強的大型地應力場計算能力,但在三維精細化建模方面能力較弱。因此,將數值仿真軟件與地質建模軟件結合使用是十分有必要的,這種方法既可以精細表征儲層特征,也可以準確地模擬的地應力場分布規律。利用“地質-工程”耦合平臺,將地質建模軟件與數值仿真軟件AiFrac整合[17]。該平臺將地質模型轉換成AiFrac計算所需的幾何數據,并將儲層地質參數信息融入AiFrac網格數據中,從而實現三維地應力場的精細化仿真。

1.1 地質模型

三維地質建模軟件能夠創建三維構造框架模型,利用地球物理和地質知識庫綜合分析油藏開發的全過程,幫助人們理解復雜的地質問題。地質建模常常采用的網格為角點網格[18]。角點網格是由4個坐標線以及8個網格節點構成的網格體,該網格體內部又可以分成i×j×k個子單元(i、j、k分別為角點網格在X、Y、Z方向上的數量),且相鄰子單元的相鄰界面節點坐標可以不完全相同,以作為斷層或者裂縫。對于單體結構單元,4條坐標線中的任何一條都可以減少到零長度,即網格單元存在著退化,一般退化的網格會出現在侵蝕等復雜地質構造處。其中,網格體在平面上可以分為i×j個子平面網格,并且單元網格的長、寬大小可變[18]。在平面外方向上將連接頂底網格點的網格面分成k層,并且該網格面可以是傾斜的。每個網格節點的坐標值可以根據頂底邊界的節點坐標值進行內插計算得出。則A、E兩點的X、Y坐標,可由頂部I點與底部J點的X、Y、Z坐標值通過線性內插的方式得到。

(1)

(2)

式中：A點坐標為(XA,YA,ZA);E點的坐標為(XE,YE,ZE);I點的坐標為(XI,YI,ZI);J點坐標為(XJ,YJ,ZJ)。

在建模時,首先,利用井分層點數據插值生成構造層面,再依據斷層模型生成包含頂面、中面以及底面的空間格架網格模型,然后,將地層層面信息加入其中,最后,通過設置垂向單元的厚度、單元的個數進行三維網格的創建以及地層網格的精細化。

1.2 AiFrac數值仿真模型

采用的軟件AiFrac基于有限-無網格法(finite element-meshfree method, FEMM)原理來模擬三維水力擴展和地應力場[19-20],該方法最早由Rajendran等[21]于2007年首次提出。有限-無網格法中,根據單元與裂隙的位置關系,將單元分為3種類型,即常規的普通單元、橋單元和裂隙單元,如圖1所示。被裂隙穿過的單元定義為裂隙單元,與裂隙單元相鄰的為橋單元,其余的單元為普通單元。同時,該方法定義兩種節點,裂隙單元上所有節點為單位分解節點,其他節點為普通的有限元節點。

圖1 有限-無網格法中的3種單元和兩種節點

對于非裂縫的普通網格單元,采用單位分解法(partition of unity method, PUM)構建位移函數。對于一個四面體單元,其計算域Ω內位移函數由4個節點的權函數和4個節點的位移得

(3)

式(3)中：wi(x)為一系列非負的權函數;i為與節

對于裂隙穿過的裂隙單元,采用Shepard公式構造單元節點權函數[22],可表示為

(4)

式(4)中：φ′ι(x)為形函數構成的子權重函數。

(5)

式(5)的含義是一個被裂縫分割的單元只有節點都處在裂縫面的一側才認為節點彼此可見。因此,對于考慮裂隙面穿過的任意單元中的計算點Pi,可以用式(6)計算φ′i。

(6)

式(6)中：φi可以由任意坐標為x∈ψΩ的計算點Pi的體積坐標計算得到,可表示為

(7)

式(7)中：v(P1P2P3P4)為頂點P1、P2、P3和P4的四面體體積;v[P(x)PiPjPk]為由單元中任意點P(x)和3個頂點Pi、Pj、Pk組成的四面體體積。

對于橋單元,權函數采用有限元中的形函數形式。

1.3 地質模型和AiFrac模型耦合

“地質-工程”耦合平臺采用C++語言編寫。C++語言同時擁有高級語言和匯編語言的優點,同時具有簡潔靈活、數據結構豐富和程序執行效率高等特點,保證了平臺能夠快速高效地將地質模型中的角點網格單元轉換為可以供AiFrac計算使用的四面體網格單元。主要使用步驟如下。

步驟1使用地質建模軟件中的“Eclipse Keywords format”格式將儲層地質體模型數據導出。

步驟2利用“地質-工程”耦合平臺的讀入功能,讀取“.GRDECL”后綴名的數據文件,并獲取SPECGRID、COORD和ZCORN等網格坐標參數,以及PERMX、PORO和ACTNUM網格屬性信息,再對讀入的數據進行處理,生成AiFrac仿真模型所需的“.inp”后綴名數據文件。

步驟3再將數據文件導入到AiFrac數值仿真軟件中,實現地質模型和數值仿真的一體化。

“地質-工程”耦合平臺運行邏輯如圖2所示,具體如下：①檢查輸入模型的有效性并讀取文件,解析角點網格模型數據,記錄角點網格模型的節點、網格和屬性信息;②判斷構建的網格模型單元的有效性,將無效單元刪除,同時移除無效節點,然后更新節點單元關系;③通過嚙合重復節點,減少生成的網格中重復節點的數量,同時將角點網格轉換成有限元網格;④將四面體網格的幾何、屬性信息輸出,形成可供AiFrac仿真軟件運行的輸入文件。

圖2 “地質-工程”耦合平臺運行邏輯

2 順北斷溶體儲層精細化建模

2.1 基本地質特征

選取塔里木盆地順托果勒隆起構造帶中的順北5號斷裂帶部分區域作為研究對象,該區域的地質構造主要受斷裂帶的控制,儲層類型為裂縫-洞穴型碳酸鹽巖儲集體。此外,順北地區的主要目標層位為一間房組和鷹山組,在早期拉張環境下,該地區沉積了厚達3 000 m的寒武系-中奧陶統臺地相碳酸鹽巖地層。

2.2 儲層力學參數的確定

通過現場采樣進行室內巖石力學實驗,以及處理分析測井資料,可以獲取彈性模量、泊松比等巖石力學參數。以往的研究成果表明：通常情況下,斷裂帶的彈性模型等參數相比儲層同性地層會有一定程度降低,而泊松比略大。目前還沒有確定斷裂帶巖石力學參數的方法,一般可以將斷裂帶的彈性模量設置為儲層同性地層楊氏模量的50%～70%[23-24]。斷裂帶的彈性模量和泊松比與其復雜程度有關,復雜程度越大,其彈性模量通常會越小,而泊松比則會越大。基于前期獲得的巖石力學實驗和測井數據,合理選取巖石力學參數[25]。具體巖石力學參數如表1所示。

表1 計算模型的巖石力學參數Table 1 Rock mechanics parameters of the calculation model

2.3 邊界條件和載荷確定

水平最大主應力方向一般是利用電成像、偶極聲波資料、古地磁及波速各向異性等方法獲得。由于斷裂構造、巖體結構和性質不同,造成地應力方向發生了顯著變化,通過平均各單井的水平最大主應力方向數據,確定研究區域現今水平最大主應力方向為NE41.4°[25]。

現今地應力大小主要可以通過室內巖心實驗、測井數據和水力壓裂法等方法獲得,通過對前期完成的單井實測數據取平均值,確定研究區的邊界荷載分別為水平最大主應力177.78 MPa和水平最小主應力為137.01 MPa,垂向應力通常等于巖石的上覆巖層壓力,可以利用密度測井估算獲得[25]。

2.4 基于地質數據體的三維仿真模型建立

基于油田構造特征、三維地震解釋和測井數據等工程勘測資料,利用地質建模軟件建立斷裂帶儲集體三維地質模型。地質模型采用平面分段、縱向分層的思路進行建立。具體而言,首先,基于地震解釋等工程資料,采用地震解釋和地質研究相結合的方式,對儲集體的斷層以及斷裂組合特征進行描述;然后,基于地震數據,建立斷裂帶儲集體空間模型;最后,結合試井、巖心、測井數據,建立儲集體屬性模型。由于斷裂帶巖體的孔隙度與儲層原巖有較大差異,因此,可以通過孔隙度來表征斷溶體,并將其與儲層原巖進行區分。三維地質模型如圖3所示。圖3(a)和圖3(b)分別為地質模型中的溶洞分布圖和斷裂帶分布圖,比例尺表征的是水平方向上的比例關系。

基于前期已構建的順北5斷裂帶的地質模型,利用“地質-工程”耦合平臺,構建基于地質數據體的三維仿真模型。由于研究區域的地質數據體三維模型尺寸大、網格數量多,使用精細數值仿真軟件直接對原尺寸網格的模型進行計算會面臨耗時較長、內存占用多等問題。因此,在對地質數據體三維模型進行網格轉換前,會對模型網格進行適當的粗化。為了方便加載模擬真實的邊界條件情形,將地質模型的周邊未建模區域進行了適當填充,對其進行加框處理,從而生成一個長方體形狀的數值仿真計算模型。模型邊界方向分別對應實際研究區域水平最大主應力方向和水平最小主應力方向。采用1.3節中所介紹的方法,首先,將地質模型按照“Eclipse Keywords format”格式導出角點網格的幾何和屬性數據,然后,利用“地質-工程”耦合平臺將數據轉換生成后綴為.inp格式的AiFrac仿真輸入文件,最后,開始AiFrac仿真計算。模型加載示意如圖4所示,在水平方向上施加擠壓應力。在垂向上施加重力。本次研究所建立的AiFrac三維仿真模型包含1.787×107個三維四面體單元、4.78×106個節點。

圖4 模型加載示意圖

3 三維地應力場分布特征

通過AiFrac模擬計算,可以得到順北5號斷裂帶斷溶體儲層地應力場分布規律。圖5(a)為最大主應力大小分布,可以發現,現今最大主應力主要以擠壓應力為主,在遠離斷裂帶的地層處,地應力場均勻分布,最大主應力的數值主要范圍為170～180 MPa,而在斷裂帶的內部,最大主應力會偏小,其數值范圍為小于170 MPa。

圖5 儲層區域主應力分布

圖5(b)為中間主應力大小分布,可以發現,現今中間主應力同樣以擠壓應力為主,在遠離斷裂帶的地層處,地應力很長均勻分布,中間主應力的數值主要范圍為170～180 MPa,比最大主應力數值略微偏小,而在斷裂帶的內部,最小主應力會偏小,其數值范圍為小于170 MPa。

圖5(c)為最小主應力大小分布圖,可以發現,現今最小主應力同樣以擠壓應力為主,在遠離斷裂帶的地層處,地應力場均勻分布,最小主應力的數值主要范圍為133～137 MPa,而在斷裂帶的內部,最小主應力會偏小,其數值范圍為小于133 MPa。

圖6為三軸主應力大小剖面圖,圖6(a)、圖6(b)、圖6(c)和圖6(d)分別為剖面位置、最大主應力、中間主應力和最小主應力。可以發現,與前面的平面三軸主應力分布規律基本一致,在斷裂帶內部,三軸主應力大小偏小,而在斷裂帶端部及斷裂帶交匯,三軸主應力偏大。同時,由于垂向應力增加,三軸主應力大小整體上隨著地層深度的增加而增加。

圖6 儲層區域應力場分布剖面

斷裂帶等非連續結構面對地應力的大小和方向都有著顯著的影響,這是影響順北地區地應力場分布的主要因素之一。由于巖石力學參數弱化,斷裂帶內部應力顯著降低,產生明顯應力釋放現象。同時,斷裂帶端部附近則產生了明顯的應力集中現象。另外,交叉等類型的斷層交匯區也容易發生應力集中現象。此外,斷裂帶的規模越大,引起的地應力擾動帶范圍越大,同時,斷裂帶的傾角、走向,也會對地應力場的分布產生影響。綜上所述,順北斷溶體儲層區域及內部、不同斷裂帶交匯區三維地應力場分布差異明顯,斷裂帶的規模、傾角、走向以及形態以及等因素在不同程度上均對地應力狀態產生影響。

4 結論

針對順北油田斷溶體儲層區域及內部、不同斷層交匯區三維應力場精細空間分布刻畫的難題,通過開發“地質-工程”耦合平臺,提出了一種精細化三維地應力場仿真模擬方法,該方法實現了地質模型和數值仿真的一體化。此外,基于“地質-工程”耦合平臺,開展了順北碳酸鹽巖儲集體三維地應力場分布規律研究。得到以下主要結論。

(1)通過開發“地質-工程”耦合平臺,可以將基于角點網格的地質模型信息轉換到AiFrac數值仿真軟件中,實現了地質模型和數值仿真一體化。該平臺結合了兩種方法的優勢：一是克服了地質模型軟件無法和力學模擬的軟件銜接問題,實現了地質建模和地應力場反演的一體化;二是在三維地應力場數值模型中保留了地質模型的儲層特征。可以快速和準確地將包含斷層等復雜地質構造的地質模型轉換成三維地應力場反演數值模型,保留地質模型中的網格和屬性信息。

(2)以順北油田為工程實例,基于開發的“地質-工程”耦合平臺,將地質模型轉換為AiFrac仿真模型,同時確定儲層力學參數以及邊界條件,實現了對碳酸鹽巖儲集體三維地應力場的精細研究。研究結果表明,斷裂帶等非連續結構面對三維地應力的大小和方向都有著顯著的影響,這是影響順北地區三維地應力場分布的主要因素之一。斷裂帶處材料參數的弱化,會導致在斷裂帶內部出現應力釋放,從而引起斷裂帶內部應力降低,同時在斷裂帶端部附近則產生明顯的應力集中現象。順北斷溶體儲層區域及內部、不同斷裂帶交匯區三維地應力場分布差異也很明顯,斷裂帶的規模、傾角、走向以及形態等因素在不同程度上都會對地應力狀態產生影響。