含水層儲氣庫注入階段蓋層力學完整性數(shù)值模擬分析

楊軍偉， 仲家銳， 賈善坡， 滕振超， 楊典森

(1. 東北石油大學 土木建筑工程學院，黑龍江 大慶 163318； 2. 哈爾濱工程大學 航天與建筑工程學院，黑龍江 哈爾濱 150001； 3. 東北石油大學 非常規(guī)油氣研究院，黑龍江 大慶 163318； 4. 中國科學院武漢巖土力學研究所 巖土力學與工程國家重點實驗室，湖北 武漢 430071 )

0 引言

天然氣以清潔、高效等特點成為現(xiàn)代化城市燃氣發(fā)展的主要方向[1-3]。地下儲氣庫的加速建設可彌補中國天然氣儲存能力不足，保障國家能源安全并完善天然氣產業(yè)鏈。根據(jù)儲層、地質等特點，地下儲氣庫主要分為油氣藏型、含水層型、鹽穴型和廢棄礦坑型。油氣藏型儲氣庫改建是最常用、最經濟的一種儲氣方法，但在大型工業(yè)城市中心和大城市附近，并非都有適合建設地下儲氣庫的枯竭油氣田，且長距離管線輸送對安全供氣的要求很高。建造含水層型地下儲氣庫所需地層結構分布較廣泛，一般在大城市周圍可以找到合適的建庫區(qū)域[4-6]。大量天然氣注入深層地質構造產生巨大的超壓和應力狀態(tài)變化，可能危及蓋層的力學完整性，導致天然氣泄漏。因此，有必要開展對含水層儲氣庫注入階段蓋層力學完整性的研究，以確保注入的天然氣能在預期壽命內安全儲存在地下含水層儲氣庫中。

人們研究含水層儲氣庫多集中于CO2封存時蓋層的穩(wěn)定性。對于阿爾及利亞的Salah天然氣項目，RUTQVIST J等[7]建立彈性三維模型， 探討CO2注入后蓋層隆起的原因和機理，但模型存在很多不確定因素，不適合實際應用，同時忽略天然氣在注入過程中產生的應力變化對蓋層剪切及拉伸破壞的影響；VILARRASA V等[8]通過實驗測試具有代表性的原位條件下CO2儲存的巖石類型，將測得的參數(shù)作為數(shù)值模擬的輸入數(shù)據(jù)，模擬CO2注入對蓋層剪切破壞的影響，但忽略注入過程對蓋層變形的影響；在恒應力條件下，F(xiàn)AVERO V等[9]研究CO2注入對蓋層流體力學性能的影響，但忽略注入過程引起的應力變化對蓋層地質力學完整性的影響；HADIAN P等[10]應用4種實驗技術，評估超臨界CO2注入對蓋層巖石物性的影響，但忽略注入過程對蓋層力學性能的影響；在彈性狀態(tài)下，VILARRASA V等[11]研究CO2注入位置和斷層滲透率對斷層穩(wěn)定性等地質力學的影響，以及注入井周圍蓋層穩(wěn)定性的影響，但未考慮蓋層參數(shù)對蓋層力學完整性的影響；KIM K等[12]應用數(shù)值模擬，研究CO2注入對蓋層巖石地質力學性質的影響，但未考慮注入速率對蓋層地質力學性質的影響；LI Chao等[13]選擇熱彈性力學本構模型，采用熱—流—固耦合模擬方法，研究深層含水層注入大量CO2后蓋層的變形和地質力學的不穩(wěn)定，但未考慮地層參數(shù)對蓋層變形和地質力學不穩(wěn)定影響的敏感性；賀凱[14]應用實驗和數(shù)值模擬相結合的方法，在CO2地質封存過程中，對儲蓋層體系的地質力學作用影響進行研究，但未考慮對蓋層力學作用影響的敏感性。對于天然氣注入過程中蓋層的地質力學變化的研究相對較少，而應用數(shù)值模擬技術模擬含水層儲氣庫天然氣注入對蓋層力學完整性的研究更少。

筆者利用ABAQUS軟件，建立含水層儲氣庫天然氣注入階段的數(shù)值模型，探討含水層儲氣庫在天然氣注入過程中蓋層的力學行為特征，分析蓋層參數(shù)、天然氣注入速率和地應力因數(shù)對蓋層力學完整性的影響，為我國含水層地下儲氣庫蓋層穩(wěn)定性的研究提供技術支持和借鑒，也為工程實際提供參考。

1 地質力學模型

1.1 計算條件

以某含水層目標為研究對象，地層為相對完整型背斜構造，儲層為砂巖，蓋層為泥巖，前期篩選結果表明在該區(qū)塊建設含水層儲氣庫成功率較高。取99塊儲層巖心樣品進行物性測試，孔隙度最小為2.5%，最大為15.4%，平均為10.0%；滲透系數(shù)最小為4.00×10-10m/s，最大為1.29×10-6m/s，平均為1.00×10-7m/s。取16塊蓋層巖心樣品進行物性測試，孔隙度最小為6.4%，最大為26.3%，平均為10.0%；滲透系數(shù)最小為1.22×10-13m/s，最大為3.75×10-9m/s，平均為1.00×10-12m/s；突破壓力大于10.00 MPa，滿足靜態(tài)密封性要求。

對儲層和蓋層巖心樣品進行三軸壓縮實驗，不同圍壓下儲層和蓋層巖石的應力—應變曲線分別見圖1和圖2。由圖1-2可知，儲層巖石的彈性模量為4.69～24.29 GPa，平均為11.31 GPa，泊松比為0.21～0.39，平均為0.29；蓋層巖石的彈性模量為5.59～9.13 GPa，平均為7.18 GPa，泊松比為0.10～0.14，平均為0.11。對巖石應力—應變數(shù)據(jù)進行處理，得到儲層和蓋層巖石對應的強度參數(shù)：內摩擦角分別為28.2°和25.0°；黏聚力分別為17.23和20.53 MPa。

圖1 不同圍壓下儲層巖石應力—應變曲線(80 ℃)Fig.1 Stress-strain curves of reservoir rocks under different confining pressures(80 ℃)

圖2 不同圍壓下蓋層巖石應力—應變曲線(常溫)Fig.2 Stress-strain curves of caprock rocks under different confining pressures(room temperature)

對儲層和蓋層部分巖心樣品進行間接拉伸實驗。儲層和蓋層試樣的間接拉伸破壞曲線見圖3。由圖3可知，儲層巖石抗拉強度為7.44～13.13 MPa，平均為10.31 MPa；蓋層巖石抗拉強度為4.15～7.44 MPa，平均為5.52 MPa。

由于地層巖心測試資料較少，參照儲層、蓋層巖石力學實驗數(shù)據(jù)測試結果及相關文獻資料[15-16]確定模型計算參數(shù)(見表1)。

圖3 儲層和蓋層試樣間接拉伸破壞曲線Fig.3 Indirect tensile failure curves of reservoir and caprock samples

表1 模型計算參數(shù)

1.2 模型建立

研究區(qū)構造相對平緩，屬于微幅構造且儲層厚度較大，為水平構造，數(shù)值模擬采用二維軸對稱模型[17]。根據(jù)地質資料建立二維模型(見圖4)，模型為矩形，長度為5 000 m，由蓋層、儲層和底層組成，蓋層、儲層、底層厚度分別為63、276、193 m(見圖4(a))。蓋層上部有深度為3 248.0 m的地質層，以等效覆蓋層表示。天然氣沿儲層中一點以0.035 m3/s的恒定注入速率向儲層注入6個月，注入點深度為3 487.0 m。

由于研究目標是蓋層地質力學變化而不是孔隙壓力本身，模擬采用單相流體流動，流體相為非必要限制條件[18]。

以ABAQUS軟件中的流—固耦合計算模塊為求解器，求解流—固耦合問題，通過黏土Biot一維固結過程的數(shù)值模擬與理論解對比，對數(shù)值模型進行驗證[19]。模型網(wǎng)格劃分采用CAX4RP單元(四結點軸對稱四邊形單元)、雙線性位移、雙線性孔隙壓力。在儲層和蓋層中，越靠近注入點網(wǎng)格劃分越細，越遠離注入點網(wǎng)格劃分越粗，二維模型有限元網(wǎng)格劃分見圖4(b)。

假設3種相互正交的主應力中有一種具有垂直方向，垂向應力的大小由上覆巖體和流體的質量決定[20]。不同深度的垂向應力可表示為

σv=ρgz,

(1)

式中：σv為垂向應力；ρ為覆蓋層平均密度；g為重力加速度；z為深度?？紤]覆蓋層巖石的自重，由式(1)可得蓋層頂部的垂向應力為83.1 MPa(ρ=2 608 kg/m3，g=9.81 m/s2，z=3 248 m)。初始水平應力為各向異性，模型的地應力因數(shù)為0.76。初始孔隙壓力梯度為10.50 MPa/km，在天然氣注入開始時，注入點的初始孔隙壓力為36.61 MPa。

數(shù)值模擬采用總孔隙壓力的方法，力學邊界：上邊界為上覆壓力；下邊界為法向位移約束；左邊界為軸對稱邊界；右邊界為法向約束。滲透邊界：上邊界為滲透邊界；下邊界不排水；左邊界不排水；右邊界為封閉邊界。

圖4 二維模型的物理模型及網(wǎng)格劃分Fig.4 Physical model and mesh generation of two dimensional model

1.3 蓋層破壞準則

巖石破壞準則為給定巖石破壞時應力狀態(tài)需要滿足的條件，在條件外巖石不能承受施加的載荷[21]。為確定蓋層在天然氣注入過程中是否穩(wěn)定，需要定義蓋層的破壞準則。如果應力狀態(tài)低于破壞準則，則蓋層是穩(wěn)定的。采用莫爾—庫倫破壞準則[22]模擬巖石破壞的表達式為

τ=c+(σn-pf)tanφ,

(2)

式中：τ為剪切面上的剪應力；c為黏聚力；σn為剪切面上的正應力；pf為孔隙壓力；φ為內摩擦角。巖石的黏聚力和內摩擦角見表1。巖石剪切面上的剪應力和正應力分別為

(3)

(4)

式(3-4)中：σ1為最大主應力；σ3為最小主應力；θ為巖石破壞面與σ3方向的夾角[23-24]。當τ≥c+(σn-pf)tanφ時，蓋層發(fā)生剪切破壞。

巖石拉伸破壞發(fā)生在拉伸應力超過巖石抗拉強度時。材料的抗拉強度是材料破裂之前能承受的最大拉應力。根據(jù)材料第一強度理論，拉伸破壞準則[25-26]為

(5)

2 數(shù)值模擬結果

為使初始應力場和模型邊界條件之間達到力學平衡，應進行地應力平衡，其計算結果的準確性決定模擬結果的可靠性。采用ODB文件導入法對模型進行地應力平衡，得到初始地應力條件下垂向位移數(shù)量級為10-7m；垂向應力為壓應力，自上而下逐漸增大；孔隙壓力自上而下隨深度線性增加，初始地應力平衡效果較好。

2.1 孔隙壓力

離注入點水平距離10 m處取一條垂向測線CE，天然氣注入后，不同注入時間的測線CE孔隙壓力變化曲線見圖5。由圖5可知，蓋層的孔隙壓力隨深度的增加而增加，最大孔隙壓力位于注入6個月時的蓋層底部(深度為3 311.0 m)，為40.92 MPa，比初始孔隙壓力增加6.15 MPa。模型最大孔隙壓力位于天然氣注入6個月時的注入點(深度為3 487.0 m)，為59.86 MPa，是初始孔隙壓力的1.635倍。

當天然氣注入6個月時，離注入點不同水平距離取4條垂向測線，得到測線AB、CD、FG、HI的孔隙壓力變化曲線(見圖6)。由圖6可知，離注入點越近的垂向測線孔隙壓力變化越大，離注入點較遠的垂向測線孔隙壓力變化較小。

圖5 不同注入時間的測線CE孔隙壓力變化曲線Fig.5 Curves of pore pressure variation of the CE under different injection time

圖6 不同水平距離的垂向測線孔隙壓力變化曲線Fig.6 Curves of pore pressure variation of the vertical measuring lines under different horizontal distances

不同注入時間的儲、蓋層界面孔隙壓力增量變化曲線見圖7。由圖7可知，最大孔隙壓力增量發(fā)生在天然氣注入6個月時的注入點附近，為6.118 MPa，隨離注入點水平距離增加，孔隙壓力增量變小。

蓋層底部測點J孔隙壓力隨注入時間變化曲線見圖8。由圖8可知，隨注入時間增加，測點J孔隙壓力增大，且孔隙壓力在天然氣注入1個月內迅速上升，變化速率較大；隨注入時間繼續(xù)增加，孔隙壓力的增加速率逐漸變小。

圖7 不同注入時間的儲、蓋層界面孔隙壓力增量變化曲線Fig.7 Curves of pore pressure variation at the reservoir-caprock interface under different injection time

圖8 蓋層底部測點J孔隙壓力隨注入時間變化曲線Fig.8 Curve of pore pressure variation of the measuring point J at the bottom of the caprock with the injection time

2.2 蓋層位移

孔隙壓力增大引起的應力變化導致儲層孔隙膨脹，儲、蓋層界面上移。天然氣注入1個月和6個月時，儲、蓋層界面位移變化曲線見圖9。由圖9(a)可知，天然氣注入6個月時，注入點附近出現(xiàn)最大垂向位移，為0.073 m，且隨離注入點距離增加，垂向位移逐漸變??；由圖9(b)可知，天然氣注入6個月時，在離注入點水平距離1.0 km處產生最大水平位移，為0.032 m，且水平位移在離注入點一定水平距離內隨距離增加而變大，超過一定水平距離后隨距離增加而變小。

圖9 不同注入時間的儲、蓋層界面位移變化曲線Fig.9 Curves of displacement at the reservoir-caprock interface under the different injection time

儲、蓋層界面最大位移隨注入時間變化曲線見圖10。由圖10可知，隨注入時間增加，儲、蓋層界面最大垂向和最大水平位移逐漸增大，且位移增加的速率在注入1個月內變化顯著，之后位移增加的速率逐漸減小。

圖10 儲、蓋層界面最大位移隨注入時間變化曲線Fig.10 Curves of maximum displacement at the reservoir-caprock interface with injection time

2.3 測線CE應力

不同注入時間的測線CE應力變化曲線見圖11。由圖11可知，蓋層在同一注入時間內垂向應力和水平應力隨深度增加先減小后增大；相同深度時，蓋層的垂向應力和水平應力隨注入時間的增加而增大，且垂向應力和水平應力未出現(xiàn)正值，即該工況下蓋層未發(fā)生拉伸破壞。

2.4 蓋層剪切破壞

采用莫爾—庫侖失效準則[22]，得到測點J在天然氣注入過程中應力狀態(tài)變化(見圖12)。由圖12可知，天然氣注入過程中，莫爾圓沒有達到蓋層強度包絡線。蓋層在天然氣注入6個月時沒有發(fā)生破壞，并保持穩(wěn)定。蓋層強度包絡線的截距較大且斜率較陡，預計莫爾圓不會觸及蓋層破壞線。如果采用更高的注入速率，則最小主應力可能變?yōu)樨撝?。因此，需要進行參數(shù)敏感性分析，確定天然氣注入速率對蓋層剪切破壞的影響。

圖11 不同注入時間的測線CE應力變化曲線Fig.11 Curves of stress at measuring line CE under different injection time

圖12 測點J應力莫爾圓狀態(tài)變化Fig.12 State change disgram of stress Mohr circle at a measuring point J

3 參數(shù)敏感性

對蓋層參數(shù)、天然氣注入速率和地應力因數(shù)[28-29]進行敏感性分析，評價參數(shù)對蓋層變形及應力變化的敏感性。參數(shù)敏感性分析輸入數(shù)據(jù)見表2，其中A0、B0、C0、D0和E0為原工況數(shù)據(jù)。取蓋層底部為水平測線，研究參數(shù)變化對測線最大位移的影響，分析各參數(shù)對蓋層變形的敏感性；在蓋層底部取測點J，研究參數(shù)變化對測點J應力的影響，討論各參數(shù)對蓋層拉伸和剪切破壞的敏感性。

表2 蓋層參數(shù)、天然氣注入速率及地應力因數(shù)敏感性分析輸入數(shù)據(jù)

3.1 蓋層彈性模量

天然氣注入6個月時，蓋層底部最大位移隨蓋層彈性模量變化曲線見圖13。蓋層底部測點J的正應力和剪應力隨彈性模量變化曲線見圖14。由圖13-14可知，蓋層彈性模量從3.59 GPa增加到10.77 GPa，蓋層底部最大水平位移減小3.760×10-3m，最大垂向位移減小2.280×10-3m；測點J的正應力增加0.41 MPa，剪應力降低0.27 MPa。因此，蓋層彈性模量對蓋層的最大水平位移和最大垂向位移影響較小，且對最大水平位移的影響比最大垂向位移的大；對蓋層拉伸和剪切破壞趨勢影響較小。

圖13 蓋層底部最大位移隨蓋層彈性模量變化曲線Fig.13 Curves of the maximum displacement at the bottom of the caprock with Youngs modulus of caprock

圖14 測點J應力隨蓋層彈性模量變化曲線Fig.14 Curves of stress at measuring point J varying with Youngs modulus of caprock

3.2 蓋層泊松比

天然氣注入6個月時，蓋層底部最大位移隨蓋層泊松比變化曲線見圖15。蓋層底部測點J正應力和剪應力隨泊松比變化曲線見圖16。由圖15-16可知，蓋層泊松比從0.10增加到0.35，蓋層底部最大水平位移增大0.640×10-3m，最大垂向位移增大0.660×10-3m；測點J正應力降低0.08 MPa，剪應力增加0.04 MPa。因此，蓋層泊松比對蓋層的最大水平位移和最大垂向位移影響較小，可以忽略不計；對蓋層拉伸破壞和剪切破壞趨勢幾乎沒有影響。

3.3 蓋層滲透系數(shù)

天然氣注入6個月時，蓋層底部最大位移隨蓋層滲透系數(shù)變化曲線見圖17。蓋層底部測點J正應力和剪應力隨滲透系數(shù)變化曲線見圖18。由圖17-18可知，蓋層滲透系數(shù)從10-14m/s增加到10-8m/s，蓋層底部的最大水平位移增大7.080×10-3m，最大垂向位移減小6.940×10-3m；測點J的正應力增加0.01 MPa，剪應力增加0.59 MPa。因此，蓋層滲透系數(shù)對蓋層的最大水平位移和最大垂向位移影響較大；對蓋層拉伸破壞趨勢幾乎沒有影響，對剪切破壞趨勢影響較小。

圖15 蓋層底部最大位移隨蓋層泊松比變化曲線Fig.15 Curves of the maximum displacement at the bottom of the caprock with Poisson ratio of caprock

圖16 測點J應力隨蓋層泊松比變化曲線Fig.16 Curves of stress at measuring point J varying with Poisson ratio of caprock

圖17 蓋層底部最大位移隨蓋層滲透系數(shù)變化曲線Fig.17 Curves of the maximum displacement at the bottom of the caprock with permeability coefficient of caprock

3.4 天然氣注入速率

天然氣注入6個月時，蓋層底部最大位移隨注入速率變化曲線見圖19。蓋層底部測點J正應力和剪應力隨注入速率變化曲線見圖20。由圖19-20可知，天然氣注入速率從0.025 m3/s增加到0.055 m3/s，蓋層底部最大水平位移增大27.010×10-3m，最大垂向位移增大62.490×10-3m；測點J的正應力降低4.30 MPa，剪應力降低0.66 MPa。因此，天然氣注入速率對蓋層最大位移、拉伸破壞和剪切破壞趨勢影響較大，且對拉伸破壞的影響比剪切破壞的大。

圖18 測點J應力隨蓋層滲透系數(shù)變化曲線Fig.18 Curves of stress at measuring point J varying with permeability coefficient of caprock

圖19 蓋層底部最大位移隨天然氣注入速率變化曲線Fig.19 Curves of the maximum displacement at the bottom of the caprock with injection rate of gas

圖20 測點J應力隨天然氣注入速率變化曲線Fig.20 Curves of stress at measuring point J varying with injection rate of gas

3.5 地應力因數(shù)

天然氣注入6個月時，蓋層底部最大位移隨地應力因數(shù)變化曲線見圖21。蓋層底部測點J正應力和剪應力隨地應力因數(shù)變化曲線見圖22。由圖21-22可知，地應力因數(shù)由0.66增加到1.26，蓋層底部的最大水平位移和最大垂向位移幾乎沒有改變；測點J正應力增加15.02 MPa，剪應力先降低7.51 MPa，后增加7.01 MPa。因此，地應力因數(shù)對蓋層最大位移幾乎沒有影響，對蓋層拉伸破壞和剪切破壞趨勢影響很大。

圖21 蓋層底部最大位移隨地應力因數(shù)變化曲線Fig.21 Curves of the maximum displacement at the bottom of the caprock with coefficient of stress

圖22 測點J應力隨地應力因數(shù)變化曲線Fig.22 Curves of stress at measuring point J varying with coeffcient of stress

蓋層滲透系數(shù)和天然氣注入速率對蓋層變形的影響顯著，天然氣注入速率和地應力因數(shù)對蓋層應力的影響顯著。因此，在實際工程實踐中，規(guī)劃深層地質構造天然氣注入方案時，應首先考慮蓋層滲透系數(shù)、天然氣注入速率和地應力因數(shù)等因素。

4 結論

(1)以擬建含水層儲氣庫為依托工程，建立含水層儲氣庫地質力學模型，采用ABAQUS軟件對天然氣注入階段的蓋層力學完整性、蓋層變形及應力敏感性進行模擬分析。

(2)注入早期的注入點附近引起的孔隙壓力變化和地質力學變化最大，是蓋層中巖石破壞風險最大的關鍵部位。

(3)對蓋層變形最為敏感的參數(shù)為天然氣注入速率，其次為蓋層的滲透系數(shù)；對蓋層應力變化最為敏感的參數(shù)為地應力因數(shù)，其次為天然氣注入速率；蓋層泊松比和模型右側邊界類型對蓋層變形和應力變化的影響不顯著。