翟 誠,孫可明
(1.山西工程技術學院,山西 陽泉 045000;2.遼寧工業大學,遼寧 錦州 121001;3.青島理工大學,山東 青島 266520)
天然氣水合物因具有儲量豐富、清潔、高效等優點,被認為是極具開發價值和良好發展前景的新型替代能源之一。勘探成果表明,中國的天然氣水合物資源主要蘊藏在南海北部的神狐海域、青藏高原陸地凍土區以及漠河盆地陸地凍土區[1-3]。降壓法因具有操作簡單、經濟、高效等優點,被認為是未來工業化開采水合物的重要方法之一。但如果防控措施不當,水合物降壓分解引起儲層失穩破壞后,會導致包括海底邊坡失穩等一系列災害效應的發生[4-5]。目前,很多國內外學者運用數值模擬方法分析天然氣水合物分解引起的儲層變形破壞問題,并取得了一些具有借鑒意義的成果[6-14]。但在降壓分解條件下,缺少影響水合物沉積層近井儲層變形破壞因素的敏感性分析。鑒于此,以中國南海神狐海域天然氣水合物降壓開采為工程背景,采用數值模擬與正交實驗相結合的方法,基于ABAQUS有限元軟件的USDFLD子程序進行二次開發,探討初始水合物飽和度、井底壓力和水合物藏所處地應力狀態對水合物沉積層近井儲層變形破壞影響的敏感程度。
建立模型時,做出如下假設:①將水合物沉積層看作含有固、液、氣三相的多孔介質,其中,水合物與沉積物骨架膠結構成多孔介質的固體骨架,甲烷氣體和水則充填在固體骨架圍成的孔隙中;②固體骨架為各向同性的彈塑性材料,其變形滿足小變形假設,且在變形過程中沉積物骨架的密度和水合物的密度不發生變化;③不考慮水合物降壓分解過程中儲層溫度的變化,即溫度對水合物沉積層有效應力和孔隙流體性質的影響忽略不計;④甲烷氣體為理想氣體,水合物沉積層中水、氣的滲流滿足Darcy定律。
水合物分解的產氣速率采用Kim等[15]建立的分解動力學方程來計算,即:
(1)
As=φShAhs
(2)
lgpe=0.0342(T-273.15)+0.0005(T-273.15)2+6.4804
(3)

根據式(1)和化學反應方程式CH4·5.75H2O=CH4+5.75H2O,可得單位體積水合物的分解速率和產水速率分別為:
(4)

變形場方程包括固體骨架的靜力平衡方程、幾何方程和彈塑性本構方程3個部分。根據有效應力原理和彈性力學理論,水合物沉積層固體骨架的靜力平衡方程可用張量形式表示為:
(5)
(6)

孔隙中氣相壓力與水相壓力之間的關系滿足如下的方程:
pg-pw=pc
(7)
式中:pc為毛管力,Pa。
采用VG模型計算毛管力,其表達式為:

(8)
式中:Swr為殘余水飽和度;Swmax為最大含水飽和度;a、m、n為VG模型參數,且m=1-1/n。
幾何方程為:
(9)
式中:εij為應變張量;u為位移,m;ui,j、uj,i為采用張量表示的不同方向位移的偏微分形式。
固體骨架增量形式的彈塑性本構方程為:
[dσ′]=[De][dεe]=
[De]{[dε]-[dεP]-[dεpl]}
(10)
式中:[dσ′]為有效應力增量矩陣;[De]為彈性模量矩陣;[dεe]為彈性應變增量矩陣;[dε]為總應變增量矩陣;[dεP]為孔隙壓力產生的壓應變增量矩陣;[dεpl]為塑性應變增量矩陣。
孔隙壓力變化產生的單向應變增量為:
(11)
式中:Ksh為固體骨架的體積模量,Pa;vsh為固體骨架的泊松比;Esh為固體骨架的彈性模量,Pa。
塑性應變增量為:
(12)
(13)

考慮水合物分解效應和固體骨架變形對孔隙流體滲流的影響,根據質量守恒方程和達西定律,可推導得到氣、水兩相的流固耦合滲流場方程。
氣相:

(14)
孔隙中氣相相對滲透率為[18]:
(15)
水相:
(16)
孔隙中水相的相對滲透率為[19]:
(17)
水合物相的質量守恒方程:
(18)

將上述數學模型與水合物沉積層彈性模量和滲透率隨有效應力及水合物飽和度變化的關系模型結合[20],再給定具體的初始條件和邊界條件,即可構成水合物降壓分解誘發儲層變形破壞的流固耦合彈塑性模型。
建立如圖1a所示的有限元模型進行正交數值模擬實驗,并對井眼附近區域進行網格細化,坐標如圖1b所示。整體模型的幾何尺寸為:邊長L=10 m,井眼半徑r=0.2 m。地質條件參考中國南海神狐海域水合物藏進行選取:海水深度為1 235 m,水合物藏距離海床的深度為182 m,儲層沉積物骨架主要由粉砂巖和黏土構成。模擬所用基本參數如表1所示[7,10,13]。
變形場邊界條件如下。
(1)σBC=σhmax,σCD=σhmin,(σBC、σCD分別為作用在BC、CD邊的應力,Pa;σhmax、σhmin分別為最大和最小水平地應力,Pa)。
(2)UAB=Uy=0,UDE=Ux=0,AE為固定端約束(UAB、UDE為AB、DE邊位移,m;Ux、Uy為x、y方向位移,m)。
滲流場邊界條件:

(2)AB、DE為滲流封閉邊界。

圖1 有限元模型示意圖

表1 數值模擬基本參數


表2 正交數值模擬實驗因素水平

表3 正交數值模擬實驗方案及結果
根據表3的實驗結果可分別計算得到以PEEQmax和塑性區范圍為評價指標的各因素水平的均值(表4、5),并進一步采用直觀分析法對表4、5進行極差分析,從而確定每個因素對近井儲層變形破壞影響的敏感程度。由表6、7的極差數據分析可知:各因素對近井儲層 的影響程度由大到小依次為井底壓力、有效主應力差、初始水合物飽和度,說明井底壓力對近井儲層變形破壞的影響最大,有效主應力差的影響次之,而初始水合物飽和度的影響最小;各因素對近井儲層塑性區范圍的影響程度由大到小依次為有效主應力差、井底壓力、初始水合物飽和度,說明有效主應力差對近井儲層塑性區范圍的影響最大,井底壓力的影響次之,而初始水合物飽和度的影響最小。

表4 以PEEQmax為評價指標的各因素水平的實驗均值

表5 以塑性區范圍為評價指標的各因素水平的實驗均值

表6 PEEQmax的直觀分析

表7 塑性區范圍的直觀分析
為了進一步說明各因素對近井儲層變形破壞和塑性區范圍影響的變化趨勢,根據表4、5的分析結果,分別作出各影響因素對近井儲層影響的直觀分析圖(圖2、3)。

圖2 各因素對近井儲層PEEQmax的影響規律

圖3 各因素對近井儲層塑性區范圍的影響規律
由表6、7和圖2、3可得到如下結論。
(1) 近井儲層發生變形破壞的程度主要受井底壓力和有效主應力差的影響。其中,井底壓力對PEEQmax的影響最為顯著,而塑性區范圍則受有效主應力差的影響最為明顯。相比較于前2個參數,初始水合物飽和度對近井儲層變形破壞的影響相對較小。
(2) 隨著初始水合物飽和度的增大,近井儲層PEEQmax和塑性區范圍逐漸變小,即發生變形破壞的程度逐漸減弱。這是因為,初始水合物的含量越多,其與沉積物骨架膠結后形成的儲層固體骨架的力學性質越好,抵抗變形破壞的能力也就越強。
(3) 井底壓力越低,近井儲層PEEQmax和塑性區的范圍越大。這是因為,井底壓力越低,水合物降壓分解的效果越明顯,近井儲層水合物分解區的范圍以及分解區儲層的彈性模量、黏聚力以及摩擦角等力學參數下降的幅度越大,其發生變形破壞的程度也就越嚴重。雖然井底壓力的降低,會使近井儲層的有效應力增大,在一定程度上抑制和削弱水合物飽和度降低引起的儲層力學參數劣化的程度,但這與水合物分解帶來的負面效應相比,后者占據主導地位。因此,近井儲層的抗剪強度總體仍表現出下降趨勢,并最終發生剪切破壞。
(4) 隨著有效主應力差的增大,近井儲層PEEQmax和塑性區的范圍越來越大。根據Mohr-Coulomb準則,在其他條件相同的前提下,近井儲層的有效主應力差越大,發生剪切破壞的程度越嚴重。另外,受水平地應力非均勻性的影響,近井儲層塑性區范圍的分布表現出明顯的各向異性,主要沿最大水平地應力方向擴展。圖4為選取的其中3組實驗的近井儲層塑性區分布等值線圖以及沿儲層AB方向PEEQ分布的對比分析圖,由圖4可明顯看出塑性區沿最大水平地應力方向的擴展更為顯著。

圖4 近井儲層塑性區分布等值線及對比分析
比較實驗因素水平變化和誤差波動對實驗結果的差異,將表3中的“空白列”作為誤差列對數值模擬的實驗結果進行方差分析,結果如表8、9所示。由表8、9可知:影響近井儲層PEEQmax和塑性區范圍的3個因素中,井底壓力和有效主應力差2個影響因素對應的P值皆小于0.05,說明它們是影響近井儲層變形破壞的顯著因素。根據F值的大小進行排序,3個因素對近井儲層PEEQmax的影響程度由大到小依次為井底壓力、有效主應力差、初始水合物飽和度;對塑性區范圍的影響程度由大到小依次為有效主應力差、井底壓力、初始水合物飽和度。這與極差分析所得結論是一致的。

表8 PEEQmax的方差分析

表9 塑性區范圍的方差分析
(1) 對水合物沉積層近井儲層變形破壞的影響程度由大到小依次為井底壓力、有效主 應力差、初始水合物飽和度,對塑性區范圍的影響程度由大到小依次為有效主應差、 井底壓力、初始水合物飽和度。
(2) 井底壓力和有效主應力差是影響近井儲層變形破壞的2個顯著因素,其中前者的影響在于水合物降壓分解引起的近井儲層力學性質的劣化,而后者對近井儲層變形破壞的影響則是因為初始地應力狀態的非均勻性。
(3) 使用降壓法進行天然氣水合物開采時,為確保安全、可控開采,水合物沉積層所處的地應力狀態和降壓時選用的井底壓力是2個需要重點考慮的因素,再結合水合物藏的地質條件、測井資料、取心資料以及水合物降壓分解效應對儲層力學參數的影響規律,對井底壓力進行優化設計,這是下一步的工作重點。