寧伏龍 ，竇曉峰 ，孫嘉鑫 ，劉志超 ，李彥龍 ，李曉東 ，趙穎杰 ，劉昌嶺，陸紅鋒，于彥江 ，李芷 ，羅強 ，曹鑫鑫

1 中國地質(zhì)大學(xué)(武漢)工程學(xué)院，武漢430074

2 科技部地球深部鉆探與深地資源開發(fā)國際聯(lián)合研究中心，武漢 430074

3 中國地質(zhì)調(diào)查局青島海洋地質(zhì)研究所，青島 266071

4 中國地質(zhì)調(diào)查局廣州海洋地質(zhì)調(diào)查局，廣州 510075

0 前言

天然氣水合物是由水分子形成的氫鍵網(wǎng)絡(luò)(籠狀結(jié)構(gòu))包裹客體天然氣分子在高壓低溫環(huán)境下形成的似冰狀結(jié)晶化合物[1-2]。自然界最常見的天然氣水合物客體分子為甲烷，因此也常被稱為甲烷水合物。天然氣水合物由于其資源儲量大、分布廣泛、清潔低碳等優(yōu)點[3-4]，被公認為是潛在的替代能源[5-6]，其中90%以上都分布在海洋區(qū)域[7]。因此，海域水合物是資源開發(fā)的重點。我國南海也蘊藏有豐富的水合物資源，經(jīng)過近20 年的不懈努力，我國在南海北部陸坡圈定了11 個水合物成礦遠景區(qū)、19 個成礦區(qū)帶、25 個有利區(qū)塊、24 個鉆探目標(biāo)區(qū)，取得了一系列重大找礦成果，極大豐富了我國的能源儲備[8-9]。尤其是2017 年在南海神狐海域成功實施了首次水合物試采[10]，令世界矚目。從5 月10 日成功產(chǎn)氣點火，至7 月9 日第一口試采井順利關(guān)井，總產(chǎn)氣量達30.9×104m3，連續(xù)試采時間世界第一[11-12]。同年，國務(wù)院正式批準(zhǔn)將天然氣水合物列為我國第173 個礦種[13]。2020 年3 月，我國海域天然氣水合物第二輪試采又迎來歷史性突破，總產(chǎn)氣量高達86.14×104m3，日均產(chǎn)氣量為2.87×104m3，實現(xiàn)了“探索性試采”向“試驗性試采”的重大跨越[4]。試采成功為推動水合物的產(chǎn)業(yè)化發(fā)展邁出了關(guān)鍵一步。然而要想實現(xiàn)產(chǎn)業(yè)化生產(chǎn)(如產(chǎn)氣量＞40 000 m3/d，生產(chǎn)10 年以上)，則必須解決制約水合物安全高效可控開采的瓶頸難題。

由于海域水合物儲層通常埋深淺，固結(jié)程度弱，鉆采過程中容易發(fā)生井壁失穩(wěn)、固井失效、出砂、海底沉降、滑坡等一系列工程與地質(zhì)問題[7,14-15]。其中開采導(dǎo)致的泥砂產(chǎn)出已成為水合物產(chǎn)業(yè)化亟需突破的瓶頸之一。所謂出砂就是指儲層開采過程中，固體顆粒主要是砂顆粒從儲層骨架脫離并隨氣、水流體運移至井筒內(nèi)產(chǎn)出的現(xiàn)象。2007 年加拿大Mallik水合物試開采井、2012 年美國阿拉斯加CO2置換開采甲烷水合物試驗井以及2013 和2017 年日本南海海槽水合物試采井都出現(xiàn)了嚴重的出砂問題[16]，極大地影響了產(chǎn)氣效率，甚至導(dǎo)致試采不得不提前終止。盡管我國南海粉砂質(zhì)水合物儲層第一次試采出砂量少，但防砂后的出水不足也導(dǎo)致降壓效果受到影響進而制約了產(chǎn)能提升。產(chǎn)生上述生產(chǎn)狀況的原因，歸根結(jié)底還是水合物儲層出砂機理與規(guī)律未全面厘清與掌握，從而導(dǎo)致防砂方法不科學(xué)，防砂技術(shù)不過關(guān)。只有揭示了出砂機理并系統(tǒng)掌握水合物儲層出砂規(guī)律，才能優(yōu)化設(shè)計防砂方案，形成科學(xué)合理的完井生產(chǎn)制度，最終實現(xiàn)水合物安全高效可控開采。

由于水合物開采現(xiàn)場試驗成本高、風(fēng)險大，室內(nèi)實驗和數(shù)值模擬成為水合物開采出砂機理、規(guī)律預(yù)測及其控制的主要研究手段。其中，數(shù)值模擬可為室內(nèi)實驗研究和現(xiàn)場試采提供重要理論依據(jù)與技術(shù)支撐，具有提高實驗效率、節(jié)約時間和材料成本等優(yōu)點，因此成為水合物出砂預(yù)測和控制優(yōu)化設(shè)計的首選方法。目前，數(shù)值模擬方法已經(jīng)在水合物開采產(chǎn)能預(yù)測方面得到了廣泛運用，比如目標(biāo)靶區(qū)產(chǎn)能評價[17-22]、不同開采方法或不同井型的產(chǎn)能對比[23-29]以及儲層增產(chǎn)改造下的產(chǎn)能分析等[30-32]。然而，針對水合物儲層出砂問題的數(shù)值模擬研究公開發(fā)表的文獻相對較少，相關(guān)研究還處于初期探索階段，急需梳理出后續(xù)需要解決的關(guān)鍵問題。

因此，本文在簡要介紹水合物開采儲層出砂行為特征的基礎(chǔ)上，從連續(xù)介質(zhì)和離散介質(zhì)兩個方面系統(tǒng)總結(jié)了水合物開采儲層出砂數(shù)值模擬研究現(xiàn)狀，并對今后水合物出砂數(shù)值模擬分析與預(yù)測研究提出了幾點建議，旨在提高出砂數(shù)值模擬準(zhǔn)確性和適用性，進而為后續(xù)水合物防砂和生產(chǎn)制度優(yōu)化設(shè)計奠定良好的技術(shù)基礎(chǔ)。

1 水合物儲層開采出砂行為分析

基于井筒的水合物開采過程實際上是井周儲層中固態(tài)水合物分解成氣態(tài)天然氣和液態(tài)水并伴隨泥砂等固體產(chǎn)出的過程。因此，與傳統(tǒng)油氣資源開發(fā)不同的是，水合物開采儲層出砂問題涉及水合物相態(tài)變化[14]，是一個復(fù)雜的熱—流—固—化多場耦合過程[33]，具有多尺度、強耦合的時空演化特征，導(dǎo)致對儲層特別是井周儲層響應(yīng)行為，如出砂預(yù)測和控制變得異常困難，成為水合物安全高效開采的瓶頸之一，也是其中的難點和堵點所在，亟待突破。

圖1 描述了水合物開采以及儲層出砂過程中的溫度—滲流—應(yīng)力—化學(xué)多場耦合作用關(guān)系。以降壓開采為例，降壓開采過程中，通過在井底抽汲地層流體使得儲層孔隙壓力降低，根據(jù)有效應(yīng)力原理，受降壓影響區(qū)域儲層的有效應(yīng)力增大，從而導(dǎo)致該區(qū)域一定范圍內(nèi)的應(yīng)力場改變(過程1、6)；反之，應(yīng)力場的改變會影響沉積物骨架變形，進而對儲層孔隙度和滲透率以及毛細管壓力產(chǎn)生影響，最終影響地層流體遷移(過程6)；當(dāng)儲層孔隙壓力低于水合物相平衡壓力時，水合物分解，固態(tài)水合物轉(zhuǎn)變?yōu)橐簯B(tài)水和氣態(tài)甲烷(過程3、7)；水合物分解釋放的大量氣體又會導(dǎo)致孔隙壓力的急劇增高，進而影響應(yīng)力場(過程1)；水合物作為儲層一部分，其分解將導(dǎo)致儲層孔隙或裂隙空間增大并提高儲層滲透率，利于壓力傳遞，促進水合物的進一步分解，極大增加了水合物儲層破壞出砂的可能性(過程3、7)；隨著水氣產(chǎn)出，孔隙壓力降低將導(dǎo)致有效應(yīng)力增加，從而擠壓沉積物骨架顆粒，進而降低孔隙度，部分抵消由于水合物分解而提高的地層滲透率，從而使水氣運移阻力增大，反過來阻礙水氣運移產(chǎn)出(過程1、6)；水合物分解是一個吸熱反應(yīng)，使得分解區(qū)溫度降低，分解區(qū)溫度降低反過來又會抑制水合物的分解，甚至引起局部水合物二次形成(過程4)；水合物分解會導(dǎo)致儲層強度降低，當(dāng)有效應(yīng)力或流體作用力超過儲層強度將導(dǎo)致儲層破壞，從而引起出砂，進而導(dǎo)致應(yīng)力松弛[35](過程9)；開采初期，井眼溫度與儲層原始溫度并不一致，二者的溫度差導(dǎo)致井周區(qū)域溫度場的改變。此外水合物相變和多相滲流傳熱同樣會引起溫度場的變化(過程4、8、10)，而溫度場的改變又會產(chǎn)生兩方面的影響：(1)改變流體的密度、黏度，影響流體的遷移(過程8、10)；(2)作為一種溫度荷載(熱應(yīng)力)直接改變儲層的應(yīng)力狀態(tài)(過程5)。

上述水合物開采儲層出砂多場耦合機理復(fù)雜，非線性強，難以在數(shù)學(xué)上獲得解析解。因而，基于連續(xù)介質(zhì)和離散單元法(簡稱離散元)的數(shù)值求解方法成為分析多場耦合作用下水合物儲層出砂機理、規(guī)律及其控制的有效手段。

圖1 溫度—滲流—應(yīng)力—化學(xué)多場耦合作用示意圖(改自Sun[18]和Fuente[34])Fig. 1 Schematic diagram of thermal-hydraulic-mechanical-chemical multi-field coupling (modified from Sun[18] and Fuente[34])

2 基于連續(xù)介質(zhì)理論的出砂數(shù)值模擬

連續(xù)介質(zhì)理論將含水合物巖土體視為連續(xù)介質(zhì)材料，基于常規(guī)多場耦合理論建立出砂分析模型，所求解方程組包括水、氣、水合物三相質(zhì)量守恒方程、能量守恒方程、力學(xué)平衡方程、動量守恒方程及一些輔助方程如力學(xué)本構(gòu)方程、滲透率控制方程、水合物聚散動力學(xué)方程、水合物相平衡條件、飽和度和孔隙度方程等，最后根據(jù)初始和邊界條件量化求解并利用出砂準(zhǔn)則進行出砂判斷?？梢哉f，水合物開采出砂研究是多場耦合模型應(yīng)用的經(jīng)典案例之一，其中，出砂準(zhǔn)則是出砂模擬的關(guān)鍵性控制條件。

2.1 摩爾—庫倫強度準(zhǔn)則

摩爾—庫倫(Mohr-Coulomb)本構(gòu)模型采用張拉剪切組合破壞準(zhǔn)則，亦稱摩爾—庫倫強度準(zhǔn)則(見式(1))。對于水合物而言，當(dāng)水合物開采井周儲層所受的拉應(yīng)力或剪應(yīng)力超過儲層固有的抗拉或抗剪強度時，儲層就會發(fā)生張拉破壞或剪切破壞，視為出砂。

目前較為成熟的水合物儲層數(shù)值模型和相應(yīng) 的 模 擬 軟 件 有MH21-HYDRES[36]、STOMPHYD[37]、CMG-STARS[21,38]、UMSICHT HyRes[39]和TOUGH+HYDRATE[40]，它們在水合物相變和流體運移方面均有很好的適用性。在此基礎(chǔ)上，耦合力學(xué)響應(yīng)分析模塊開發(fā)了溫度—滲流—力學(xué)—相變四場耦合的水合物開采綜合響應(yīng)模型，如TOUGU+HYDRATE+FLAC3D[41]和TOUGH+HYDRATE+Biot[42]， 對應(yīng)的數(shù)值模擬結(jié)果更貼近工程背景。如孫嘉鑫等人[43]借 鑒Moridis[40]、Rutqvist[41]等研究工作，構(gòu)建了TOUGU+HYDRATE+FLAC3D耦合模型并采用摩爾-庫倫理想彈塑性本構(gòu)模型及對應(yīng)強度準(zhǔn)則(式(1))對南海神狐海域W19 站位水合物藏鉆采過程中的力學(xué)響應(yīng)特性進行了分析，并對出砂風(fēng)險進行了初步評估。模擬結(jié)果表明：在儲層尺度和短期模擬條件下，采用逐步降壓開采方案，儲層有效主應(yīng)力演化路徑始終不會超過峰值強度曲線，見圖2，即南海神狐海域W19 站位出現(xiàn)地層屈服破壞出砂的可能性不大。式(1)中：σ1和σ3分別為最大和最小主應(yīng)力，單位MPa；Nθ為內(nèi)摩擦角θ的函數(shù)，Nθ=(1+sinθ)/(1-sinθ);σmax=c/tanθ；c為黏聚力是水合物飽和度SH(下標(biāo)H表示水合物)的函數(shù)，單位MPa，c(SH)=c(SH=0)+(1-sinθ)/(2cosθ)αSHβ,其中c(SH= 0)是地層中無水合物時的黏聚力，α和β分別是甲烷水合物的恒定材料參數(shù)，α控制強度增加的數(shù)值，無量綱，β是水合物飽和度的增加程度，無量綱；σt為抗拉強度，單位MPa，材料抗拉強度不超過σmax。

圖2 降壓開采過程中，井周3 處不同深度對應(yīng)的有效主應(yīng)力路徑演化情況(mbsf—meter below seafloor)[43]Fig. 2 Evolution of effective principal stress path corresponding to different depths of the well during the depressurization(mbsf—meter below seafloor)[43]

圖3 降壓開采過程中，距開采井水平距離為0.6 m和9.3 m處對應(yīng)的有效主應(yīng)力路徑演化情況[42]Fig. 3 Evolution of effective principal stress path at the position of 0.6 m and 9.3 m horizontally away from well during the depressurization[42]

同樣，Jin等人[42]通過構(gòu)建TOUGH+HYDRATE+Biot耦合模型對南海神狐海域SH2 站位水合物藏在水平井鉆采過程中的力學(xué)響應(yīng)特性進行了分析，發(fā)現(xiàn)儲層有效主應(yīng)力路徑與Zhang等[44]提出的Mohr-Coulomb強度曲線非常靠近(圖3)，理論上開采過程中存在出砂風(fēng)險。Zhang等[44]提出的強度準(zhǔn)則是根據(jù)粉質(zhì)黏土三軸測試數(shù)據(jù)建立的，與我國南海泥質(zhì)粉砂儲層物性相近。但是有效主應(yīng)力路徑與Rutqvist等[41]采用的根據(jù)日本南海海槽砂質(zhì)儲層性質(zhì)建立的Mohr-Coulomb強度準(zhǔn)則曲線則相去甚遠(圖3)?？梢?，儲層屈服判斷以及出砂風(fēng)險評估與沉積物物性息息相關(guān)。

2.2 侵蝕準(zhǔn)則

基于侵蝕力學(xué)發(fā)展而來的侵蝕準(zhǔn)則在常規(guī)砂巖油藏開采出砂模擬研究方面應(yīng)用較為普遍[45]?；谇治g準(zhǔn)則建立的出砂模型一般假設(shè)地層砂在流體的動力作用下發(fā)生侵蝕，固體砂轉(zhuǎn)變?yōu)榱骰斑M而運移產(chǎn)出，砂粒運移過程中可能出現(xiàn)沉積和再次舉升，如圖4 所示。

針對水合物儲層，Uchida[46]等在Papamichos與Malmanger所建出砂模型的基礎(chǔ)上考慮水合物影響，首次提出了適用于水合物地層的砂脫離數(shù)學(xué)模型(式(2))：

式(2)中，mssi(下標(biāo)ssi表示固體砂狀態(tài))是固體砂的單位體積質(zhì)量，單位kg/m3；ω2是砂脫離速度的控制參數(shù)，單位h-1；Mdtc(上標(biāo)dtc表示脫離)是砂脫離勢能，無量綱；sw、sg分別是水、氣流體的飽和度 ，無量綱；iw、ig分別是水、氣流體的水力梯度，無量綱，和分別是水、氣引起砂脫離的臨界水力梯度，無量綱；H(x)是單位階躍函數(shù)(即當(dāng)x≤0 時，H(x)=0；當(dāng)x＞0 時，H(x)=1)，該方程大括號內(nèi)第一部分代表了水砂混合物中砂粒質(zhì)量比例，第二部分代表了氣砂混合物中砂粒質(zhì)量比例。

水合物儲層砂脫離模型(式(2))的亮點在于：

(1)使用臨界水力梯度表征出砂驅(qū)動力進行砂脫離判斷，并將臨界水力梯度與水合物飽和度相關(guān)聯(lián)，假設(shè)臨界水力梯度隨水合物飽和度的增大而指數(shù)上升。關(guān)系式(以水相為例)如下：

圖4 地層侵蝕出砂示意圖Fig. 4 Schematic diagram of sand production by erosion

式(3)中，沒有水合物影響的臨界水力梯度，無量綱；ω3是與水合物相關(guān)的臨界水力梯度增長因子，無量綱；SH是水合物飽和度，無量綱。

(2)引入脫砂勢能概念，并假設(shè)砂脫離勢能與剪切變形呈線性相關(guān)，關(guān)系式如下：

式(4)中，dMdtc+是砂脫離勢能增量，無量綱；dMdtc-是砂脫離勢能減小量，無量綱；ω4是與剪切變形相關(guān)的脫砂勢能增量因子，無量綱；dεd是偏應(yīng)變增量，無量綱；mssi0是初始的原位砂質(zhì)量(下標(biāo)ssi0 表示初始固體砂狀態(tài))，單位g；是初始的砂脫離勢能，取值0。此外，該模型還考慮了出砂過程中砂顆粒運動狀態(tài)的轉(zhuǎn)變，包括砂脫離運移后的沉降和再次舉升。模擬結(jié)果表明(圖5(a))，控制降壓速率能夠有效降低開采出砂量。模型的不足之處在于假設(shè)的臨界水力梯度隨水合物飽和度的增大而指數(shù)上升可能不完全適用于分解時水合物飽和度降低過程，同時模型還引入了6 個未經(jīng)數(shù)學(xué)推導(dǎo)和實驗支撐的參數(shù)[35]。最近，Uchida等人[47]將出砂模型中引入的6 個參數(shù)對出砂量的相對影響進行了評價，結(jié)果表明出砂量主要由ω3和ω4主導(dǎo)。此外，Uchida等人[48]還對均質(zhì)和泥砂互層兩種類型的含水合物沉積物出砂情況進行了模擬，發(fā)現(xiàn)與均質(zhì)情況相比，存在泥砂互層時，泥砂界面處的剪切變形以及水力梯度均較大，相應(yīng)的出砂量也較大且出砂持續(xù)時間更久(圖5(b))。

Akaki和Kimoto[49]基于多相混合物理論，建立了一個預(yù)測水合物分解過程中包括侵蝕出砂在內(nèi)的化學(xué)—熱—力多場耦合行為的數(shù)值模型。該模型將水相和固相間的相互作用力作為侵蝕出砂驅(qū)動力(式(5)和(6))，并考慮到侵蝕行為受顆粒粒徑控制，建立了顆粒粒徑與起蝕條件以及侵蝕質(zhì)量轉(zhuǎn)化速率間的本構(gòu)關(guān)系。針對甲烷水合物的降壓開采模擬結(jié)果表明，侵蝕出砂情況下的產(chǎn)氣量比不考慮侵蝕出砂情況下的產(chǎn)氣量更大，這是因為侵蝕出砂會增大儲層滲透率，有利于壓力傳遞。侵蝕出砂還會促進土體骨架應(yīng)力和體積應(yīng)變的增長。此外，模擬結(jié)果也很好地展現(xiàn)了顆粒粒徑對出砂的影響。

式(5)和(6)中：W(SH)是水相與固相間的作用力，單位kN；是水相與砂相間的作用力，單位kN；是水相與水合物相間的作用力，單位kN；nW是水相體積分數(shù)，無量綱；kW是水相滲透率系數(shù)，單位m/s；vW是水相速度矢量，單位m/s；vSH是固相速度矢量，單位m/s；WW是水相和固相間的相對速度單位m/s；R是侵蝕驅(qū)動力，單位kN/m3；nflow是水流方向單位矢量，無量綱；Rcr是臨界侵蝕驅(qū)動力，當(dāng)R＞Rcr時，發(fā)生侵蝕，單位kN/m3；Rcr0是參考臨界驅(qū)動力，單位kN/m3；h(S) 是表征抗侵蝕能力的函數(shù)，與應(yīng)力水平、儲層土體結(jié)構(gòu)性質(zhì)以及水合物飽和度等因素相關(guān)，需要根據(jù)實驗數(shù)據(jù)確定，無量綱。

圖5 (a)不同降壓速率下的出砂量和產(chǎn)氣量對比(P8t2 表示2 d內(nèi)降壓8 MPa，降壓速率為4 MPa/d；H5t2P8t4 表示井內(nèi)注熱升溫5 K并在第2～4 d內(nèi)降壓8 MPa，降壓速率為4 MPa/天；P8t4、P8t6、P8t8、H5t2P8t6、H5t2P8t8、H5t2P8t10 同理。)[46](b)均質(zhì)和泥砂互層型水合物沉積物在30 d內(nèi)的單位深度出砂量對比[48]Fig. 5 (a) Comparison of sand production and gas production with different depressurization rates (P8t2,the depressurization rate of 4 MPa/d,that is 8 MPa drawdown in 2 days; H5t2P8t4,well temperature raised 5 K and throughout the operation,pressure reduction 8 MPa in 2-4 days with the depressurization rate of 4 MPa/d; P8t4、P8t6、P8t8、H5t2P8t6、H5t2P8t8、H5t2P8t10 have the same meanings)[46] (b) Comparison of sand production of homogeneous and interbedded hydrate-bearing sediments over a 30 day period[48]

圖6 (a)不同控砂標(biāo)準(zhǔn)下的產(chǎn)能對比和(b)出砂量對比[38]Fig. 6 (a) Comparison of gas production and (b) sand production applying the different sand control criteria[38]

Yu等[38]借助CMG-STARS軟件建立了泥質(zhì)粉砂水合物儲層開采過程中泥砂—水合物—氣—水多相流動數(shù)值模型，其中泥砂侵蝕脫離速率與臨界水流速相關(guān)(式(7))。之后對不同防砂設(shè)計方案下的產(chǎn)能與出砂量進行了模擬評估，見圖6。發(fā)現(xiàn)防砂設(shè)計對于水合物藏的持續(xù)性開采具有至關(guān)重要的影響，不恰當(dāng)?shù)姆郎霸O(shè)計會導(dǎo)致近井壁區(qū)域地層大量出砂，進而使得黏土堵塞產(chǎn)氣通道，降低產(chǎn)能。模擬獲得的最優(yōu)控砂標(biāo)準(zhǔn)為Dp=3～4 μm(Dp表示控砂標(biāo)準(zhǔn)參數(shù)，只有尺寸小于Dp的砂?？梢酝ㄟ^防砂裝置流入井眼)。

式(7)中：msi(下標(biāo)si表示砂)是每個控制體積內(nèi)固體砂的質(zhì)量，單位kg/m3；ω是泥砂脫離速度的控制參數(shù)，單位d-1；U是水流速，單位m/d；Ucri是引起泥砂脫離的臨界水流速，單位m/d。

2.3 等效塑性應(yīng)變準(zhǔn)則

油氣開采過程中，近井壁儲層在地應(yīng)力和孔隙壓力的聯(lián)合作用下易發(fā)生屈服破壞，這主要是因為儲層應(yīng)力水平超過彈性范圍后，儲層產(chǎn)生一定塑性變形所致[50]。儲層達到塑性狀態(tài)時仍能抵抗一定荷載，地層不出砂，只有當(dāng)儲層塑性應(yīng)變超過極限值即臨界塑性應(yīng)變后，儲層才會破壞[33]。

Yan等人[51]將等效塑性應(yīng)變作為出砂判斷準(zhǔn)則(式(8)和(9))，并利用TOUGH+HYDRATE軟件結(jié)合有限元軟件ABAQUS建立二維裸眼完井模型并求解了自定義水合物開采出砂多場耦合方程組，系統(tǒng)模擬分析了水合物分解效應(yīng)、開采壓差、井底溫度、儲層滲透率、初始水合物飽和度以及各向異性原位應(yīng)力等敏感性參數(shù)對出砂的影響規(guī)律。模擬結(jié)果表明，井眼效應(yīng)導(dǎo)致的應(yīng)力集中和水合物分解效應(yīng)引起的地層強度降低是儲層出砂的根本原因。

式(8)和(9)中，εp是等效塑性應(yīng)變，無量綱是3 個方向的塑性主應(yīng)變，無量綱。

上述構(gòu)建的出砂數(shù)值模型在表征近井壁儲層由于出砂而發(fā)生漸進破壞方面存在局限性(圖7(a))。采用任意拉格朗日—歐拉自適應(yīng)網(wǎng)格劃分技術(shù)(ALE)能夠有效解決這一問題。該方法能夠使網(wǎng)格獨立于材料運動，即使出砂過程中破壞或被侵蝕的儲層介質(zhì)被移除導(dǎo)致儲層發(fā)生嚴重變形(圖7(b))，出砂數(shù)值模型依舊能夠通過自適應(yīng)調(diào)整而保持高質(zhì)量的網(wǎng)格[52]。隨著出砂的發(fā)生，侵蝕面的邊界條件也隨之更新。張懷文等人[53]將侵蝕準(zhǔn)則與等效塑性應(yīng)變準(zhǔn)則相結(jié)合(式(10))并采用ALE方法進行了水合物開采出砂模擬研究，表明水合物地層的最大等效塑性應(yīng)變以及破壞出砂風(fēng)險均隨開采時間的增加而增大。此外還將出砂現(xiàn)象分為輕微出砂、暫穩(wěn)定出砂、快速出砂和后穩(wěn)定出砂4 個階段。

式(10)中，V為地層出砂速度，單位kg/(m2·s)；φ為地層孔隙度，無量綱；C1為流體中顆粒密度，單位kg/m3；vw為孔隙流體速度，單位m/s；η為地層出砂系數(shù)，無量綱，取決于地層的等效塑性應(yīng)變，同時η1最大不超過η2；εp′為地層最大可承受塑性應(yīng)變，無量綱；η1、η2為實驗系數(shù)，無量綱。

圖7 (a)未采用ALE自適應(yīng)網(wǎng)格的水合物開采出砂模型井周等效塑性應(yīng)變分布情況[51](b)采用ALE自適應(yīng)網(wǎng)格的出砂模型井周等效塑性應(yīng)變分布情況[53]Fig. 7 (a) Distribution of equivalent plastic strain around well in sand production model of hydrate production without ALE adaptive meshing[51] (b) with ALE adaptive meshing[53]

2.4 力平衡準(zhǔn)則

作用于物體的力或力距之和等于零的狀態(tài)，稱為力平衡。

針對水合物儲層開采出砂問題，劉浩伽等[54]考慮水合物的孔隙填充賦存模式，基于水合物分解區(qū)地層松散沉積物球型顆粒堆積假設(shè)，通過分析水合物分解前緣和分解區(qū)內(nèi)部地層砂粒的受力情況(圖8)，并根據(jù)力矩平衡條件建立了松散沉積物中地層砂粒啟動運移臨界流速的計算模型，砂粒啟動運移的臨界條件就是力矩平衡狀態(tài)被打破。關(guān)系如下：

(1)水合物分解前緣地層砂粒啟動運移條件

砂粒啟動運移臨界流速

(2)分解區(qū)內(nèi)部地層砂粒啟動運移條件

砂粒啟動運移臨界流速

其中

式(11)～(17)中，F(xiàn)g為微粒的浮重，單位N；Fn表示砂粒間的相互擠壓力，是上覆地層壓力沿顆粒接觸面法線方向的分量，單位N；Fh是液流滲流產(chǎn)生的推動力，單位N；Fd為水流的牽引力，單位N；Fe為內(nèi)聚力，單位N；Ff是砂粒間的摩擦力，單位N；Fl為流體的舉升力，單位N；rs為地層砂顆粒的半徑，單位m；μ為流體動力黏度，單位Pa·s；ρ為流體中地層砂微粒密度，單位kg/m3；ρ液為驅(qū)替流體的密度，單位kg/m3；h為儲層埋藏深度，單位m；H為海底之上的水深，單位m；φ為水合物分解區(qū)孔隙度，無量綱；ω為流體牽引力常數(shù)，無量綱；χ為舉升系數(shù)，無量綱；k為儲層滲透率，單位m2。

敏感性因素分析結(jié)果表明：水合物分解前緣和分解區(qū)內(nèi)部微粒的臨界流速均隨著粒徑、排布角度、顆粒摩擦系數(shù)的增大而增大。膠結(jié)物性質(zhì)、流體中的電解質(zhì)類型以及濃度等參數(shù)也均對臨界出砂流速產(chǎn)生一定的影響。

董長銀等人[55]考慮水合物的膠結(jié)賦存模式，構(gòu)建了水合物儲層砂顆粒結(jié)構(gòu)模型。該模型利用水合物強度比例系數(shù)和影響系數(shù)表征水合物分解對膠結(jié)強度的影響，并以力平衡準(zhǔn)則作為出砂判斷條件，即當(dāng)砂脫離動力大于阻力時，骨架砂顆粒剝離。模擬得到了前端類蚯蚓洞與后端連續(xù)垮塌的復(fù)合出砂形態(tài)。

圖8 分解區(qū)松散沉積物中砂粒受力示意圖(a)分解前緣(b)分解區(qū)內(nèi)部[54]Fig. 8 Force diagram of sands in unconsolidated sediments of hydrate decomposed zone (a) at dissociation front (b) inside dissociation front[54]

3 基于離散介質(zhì)理論的出砂數(shù)值模擬

連續(xù)介質(zhì)模型雖然考慮了多場耦合作用，但往往忽略了泥砂產(chǎn)出的離散特征以及顆粒間的作用力對出砂的影響。常規(guī)砂巖油藏出砂模擬研究結(jié)果表明，基于離散介質(zhì)理論的離散單元法建立的出砂預(yù)測模型能夠不受變形量限制直觀地模擬砂粒從地層骨架脫離、運移的全過程，從細觀尺度揭示儲層出砂機制與規(guī)律。然而，目前采用離散元模擬預(yù)測水合物儲層開采出砂的文獻報道還不多。

與常規(guī)砂巖油藏出砂離散元模型相比，水合物儲層出砂離散元模型必須綜合考慮以下3 個核心要素：一是水合物賦存模式仿真；二是水合物分解效應(yīng)；三是流固耦合作用。

3.1 水合物賦存模式仿真

自然界中的水合物多以孔隙填充、膠結(jié)、包裹和骨架支撐4 種形式賦存于海洋或大陸凍土沉積物中[56-61]，見圖9。

(1)孔隙填充型：水合物在沉積物骨架顆粒邊界上成核，并在孔隙空間內(nèi)自由生長。但是，孔隙內(nèi)水合物不會將周圍骨架顆粒膠結(jié)起來。

(2)膠結(jié)型：水合物在沉積物骨架顆粒間成核生長，并將骨架顆粒膠結(jié)在一起。

(3)包裹型：水合物在沉積物骨架顆粒表面及彎液處成核生長，最終將骨架顆粒包覆[62]。

(4)骨架支撐型：水合物橋接鄰近的顆粒，并通過成為承重骨架的一部分而為整體骨架提供力學(xué)穩(wěn)定性。當(dāng)水合物飽和度超過SH= 25%～40%時，孔隙填充型水合物通常會轉(zhuǎn)化為骨架支撐型水合物[63-65]。

不同的水合物賦存模式所對應(yīng)的沉積物宏觀物理力學(xué)性質(zhì)不同[57,66]，是影響水合物儲層開采出砂的重要因素之一。因此有必要根據(jù)實際儲層細觀結(jié)構(gòu)，建立不同分布模式下的含水合物沉積物開采出砂模型，系統(tǒng)性地對水合物開采過程中的出砂問題進行定性和定量的描述。雖然基于離散單元法的水合物開采出砂研究很少報道，但通過離散單元法構(gòu)建數(shù)值力學(xué)實驗如雙軸壓縮實驗(二維)和三軸壓縮實驗(三維)研究水合物對沉積物力學(xué)性質(zhì)的影響取得了長足的發(fā)展，其中不同分布模式的含水合物沉積物試樣離散元制備方法是水合物儲層出砂模型構(gòu)建的基礎(chǔ)。

圖9 水合物賦存模式Fig. 9 Pore habits of hydrate in sediments

圖10 不同賦存模式下的含水合物沉積物數(shù)值試樣示意圖(a)沉積物骨架(b)包裹型(c)孔隙填充型—水合物顆粒(d)孔隙填充型—水合物團簇(e)膠結(jié)型—水合物顆粒(f)混合型(g)膠結(jié)型—膠結(jié)接觸(改自Katagiri[73]、Yu[70]、賀潔[68])Fig. 10 Schematic diagram of numerical samples of hydrate-bearing sediments under different pore habits (a) sediment skeleton (b) grain-coating (c) pore-filling,hydrate particle (d) pore-filling,hydrate cluster (e) cementation,hydrate particle (f)combination (g) cementation,cementing contact (modified from Katagiri[73]、Yu[70]、HeJie[68])

3.1.1 孔隙填充型

Brugada等[65]提出了“Soil-＞Hydrate”和“Soil+Hydrate”兩種不同的孔隙填充型含水合物沉積物數(shù)值成樣方法。具體步驟如下：

(1)“Soil-＞Hydrate”方法，首先采用“半徑膨脹法”在一定區(qū)域內(nèi)生成指定粒徑級配、數(shù)目和初始孔隙度φ1的沉積物骨架顆粒(圖10(a))，然后通過伺服機制將沉積物骨架顆粒在一定應(yīng)力條件下固結(jié)壓實，達到目標(biāo)孔隙度φ2(φ2＜φ1)，接著，在上述區(qū)域內(nèi)隨機生成水合物顆粒(水合物顆粒粒徑小于骨架顆粒粒徑，一般差一個量級)，并要求水合物顆粒與已有顆粒不產(chǎn)生重疊，如此循環(huán)，直到指定的水合物顆粒數(shù)目滿足所需的水合物飽和度(圖10(c))。

(2)“Soil+Hydrate”方法的不同之處在于，在生成孔隙度φ1的沉積物骨架顆粒之后，接著按“Soil-＞Hydrate”方法生成水合物顆粒，最后再通過伺服機制將沉積物骨架顆粒和水合物顆粒一起固結(jié)壓實到目標(biāo)孔隙度φ2。上述兩種方法均可通過改變水合物顆粒生成數(shù)目來制備不同飽和度的含水合物沉積物(圖10(c))?！癝oil+Hydrate”方法對應(yīng)含水合物沉積物的室內(nèi)制備方法——預(yù)拌合法。

賀潔和蔣明鏡[67-68]在上述方法的基礎(chǔ)上，提出了采用水合物團簇制樣的新方法。首先采用“分層欠壓法[69]”在一定區(qū)域內(nèi)生成指定粒徑級配、數(shù)目和初始孔隙度φ1的沉積物骨架顆粒。然后，在上述沉積物骨架中隨機生成指定數(shù)目和粒徑Ra的圓球顆粒，接著，逐一刪除圓球顆粒，并在半徑為Ra的球體范圍內(nèi)生成單一顆粒半徑為r的水合物小顆粒，最后通過平行膠結(jié)模型施加強膠結(jié)作用將這些水合物顆粒凝聚成為整體，形成水合物團簇，即以水合物團簇替換圓球顆粒(圖10(d))。該方法的優(yōu)點是能夠展現(xiàn)水合物自身破壞對沉積物力學(xué)性質(zhì)的影響。

如前所述，當(dāng)模型中水合物顆粒足夠多，即水合物飽和度增大到某一特定值時，水合物小顆粒就有可能承擔(dān)骨架應(yīng)力，也就是說孔隙填充型水合物轉(zhuǎn)化為骨架支撐型水合物。

3.1.2 膠結(jié)型

Yu等[70]通過膠結(jié)模型將顆粒膠結(jié)以制備膠結(jié)型含水合物沉積物。首先生成沉積物骨架顆粒，然后，在給定飽和度下將所需水合物顆粒平均分為6 組，分階段施加6 個不同方向的重力(類似的，也可以給顆粒施加不同方向的速度[71])以實現(xiàn)水合物顆粒在膠結(jié)過程中的均勻分布。每一階段水合物顆粒都會沿某一特定重力方向沉降，直至與骨架顆粒接觸而停止并逐步累積。需要滿足的條件是保證水合物顆粒與骨架顆?；蚺c其他水合物顆粒至少有兩次接觸。一旦模型體系達到平衡狀態(tài)，就通過平行膠結(jié)模型在所有的骨架顆粒與水合物顆粒接觸點以及水合物與水合物顆粒接觸點形成膠結(jié)(圖10(e))。

Vinod等[72]則在模型中引入了一種“虛擬水合物”方法，即將儲層砂粒以顆粒形式直觀地表示，而水合物則由顆粒間的“鍵合”來表示。待沉積物骨架顆粒生成后，指定骨架顆粒間接觸模型為平行膠結(jié)模型即可(圖10(g))。水合物飽和度的差異通過平行膠結(jié)模型中的膠結(jié)強度參數(shù)來反映。

3.1.3 包裹型

Katagiri等[73]在已生成的沉積物骨架顆粒上將其半徑增大來制備包裹型含水合物沉積物，半徑增量所形成的的殼層代表水合物包裹層(圖10(b))。利用幾何學(xué)方法計算半徑增量所形成的殼層總體積與孔隙體積的比值即為水合物飽和度。

蔣明鏡等[62,74]同樣考慮將水合物簡化為骨架顆粒表面的均勻包覆層，但是在數(shù)值試樣中并沒有將水合物層顯示刻畫，而是通過改進已有離散元膠結(jié)模型并施加于骨架顆粒上以等效顆粒接觸處水合物的膠結(jié)作用，能夠?qū)崿F(xiàn)水合物的傳遞作用以及增強顆粒的抗轉(zhuǎn)動能力。

3.1.4 混合型

楊期君和趙春風(fēng)[75]提出了孔隙填充和膠結(jié)混合型的含水合物沉積物制備方法，該方法與“Soil-＞Hydrate”方法基本相同，差別在于后者在生成水合物顆粒時首先將骨架顆粒位置固定并縮小，然后在試樣孔隙中生成指定數(shù)目的水合物小顆粒，接著將骨架顆粒直徑逐漸膨脹到原有直徑并取消對其的固定，將模型計算至初始平衡。最后，在骨架顆粒與水合物顆粒接觸點以及水合物與水合物顆粒接觸點分配平行膠結(jié)模型(圖10(f))。

周博等[76]采用了更加簡化的成樣方法。首先通過計算將滿足飽和度要求的水合物小顆粒轉(zhuǎn)化為砂粒粒徑級配曲線中的一部分，然后根據(jù)初始孔隙度和粒徑級配曲線同時生成試樣所需的骨架大顆粒和水合物小顆粒，接著，平衡顆粒消除顆粒間重疊量，最后施加固結(jié)壓力，并在模型計算平衡后在骨架顆粒之間分配線性接觸模型，在骨架顆粒與水合物顆粒接觸點以及水合物與水合物顆粒接觸點分配平行膠結(jié)模型(圖10(f))。

當(dāng)然，室內(nèi)實驗結(jié)果表明，不同的水合物飽和度以及水合物生成方式會產(chǎn)生特定賦存類型的含水合物沉積物。因此，在進行數(shù)值模擬時也應(yīng)該遵循室內(nèi)實驗規(guī)律合理選擇建模方式。

3.2 水合物分解效應(yīng)

開采過程中，水合物分解是出砂的重要誘因之一。借鑒常規(guī)油氣出砂觀點，在膠結(jié)型水合物儲層開采過程中，隨著水合物的分解，地層的膠結(jié)物質(zhì)逐漸減少，膠結(jié)性能不斷下降，砂粒受到的運動學(xué)約束隨之減弱[33]。因此，膠結(jié)型水合物分解可視為儲層顆粒間黏結(jié)強度的退化過程?？紤]到這一點，蔣明鏡等[77-78]針對膠結(jié)型水合物儲層，建立了考慮溫度—水壓—力學(xué)特性綜合影響的含水合物沉積物二維微觀膠結(jié)模型，見式(18)～(23)。該膠結(jié)模型能夠反映溫度和水壓力對于水合物膠結(jié)剛度和強度的耦合影響：當(dāng)位于水合物生成區(qū)域內(nèi)的溫壓點向相平衡曲線靠近，即L減小時，水合物膠結(jié)弱化，膠結(jié)剛度和強度隨之減??；當(dāng)溫壓點剛好位于相平衡曲線上或脫離相平衡區(qū)域進入水合物分解區(qū)域時，水合物則即刻分解導(dǎo)致膠結(jié)強度和剛度為0。

具體公式如下：

圖11 TOUGH+HYDRATE+FLAC3D+PFC3D耦合分析流程Fig. 11 Coupling analysis progress of TOUGH+HYDRATE+FLAC3D+PFC3D

式(18)～(23)中：Rcb是水合物膠結(jié)抗壓強度，單位N/m；Rtb是水合物膠結(jié)抗拉強度，單位N/m；B是水合物膠結(jié)寬度，單位m，與水合物飽和度存在定量關(guān)系；1 表示垂直于紙面方向的單位長度；qmax,c是水壓σw時水合物的壓縮峰值偏應(yīng)力，單位MPa；qmax,t是水壓σw時水合物的拉伸峰值偏應(yīng)力，單位MPa；σw是水壓，單位Pa；T是溫度，單位K；是水合物接觸法向剛度，單位N；E是水合物彈性模量，單位Pa；h是膠結(jié)水合物的平均厚度，單位m；pa是一個大氣壓值，取值1.01×105Pa；L是引入的反映水合物賦存溫壓點到其相平衡曲線距離的無量綱參數(shù)；α1和α2分別是qmax,c隨L線性擬合關(guān)系線的斜率和縱軸截距，無量綱；α3和α4分別是qmax,t隨L線性擬合關(guān)系線的斜率和縱軸截距，無量綱；β1和β2分別是E/pa隨L線性擬合關(guān)系線的斜率和縱軸截距，無量綱。

寧 伏 龍 等[79]通 過TOUGH+HYDRATE+FLAC-3D+PFC3D多軟件耦合構(gòu)建了基于連續(xù)—離散介質(zhì)理論的水合物開采儲層響應(yīng)與出砂數(shù)值模型。其分析流程如圖11 所示。首先借鑒Rutqvist 等人[41]工作，利用勞倫斯伯克利國家實驗室開發(fā)的水合物—熱—流耦合開采模擬軟件TOUGH+HYDRATE與ITASCA公司開發(fā)的巖土體力學(xué)穩(wěn)定分析軟件FLAC3D進行耦合，建立熱—流—固(THM)多場耦合水合物開采數(shù)值模型并進行儲層響應(yīng)分析。然后，將模擬獲取的應(yīng)力、流速和水合物飽和度等數(shù)據(jù)作為邊界和初始條件傳遞給離散元模擬軟件PFC3D建立的細觀出砂模型進行耦合計算，最后分析出砂機理、規(guī)律并預(yù)測出砂量。

為反映水合物分解效應(yīng)，基于上述連續(xù)—離散介質(zhì)分析模型引入了黏結(jié)腐蝕概念，即在細觀出砂模型中通過編寫FISH子程序使平行膠結(jié)模型中的黏結(jié)強度參數(shù)(黏結(jié)法向強度pb_ten和黏結(jié)細觀黏聚力pb_coh隨水合物飽和度降低而遞減，從而等效模擬水合物的分解過程及其對儲層力學(xué)性質(zhì)的影響。其假定沉積物強度隨水合物飽和度下降而線性降低，公式如下：

式(24)和(25)中：pb_ten(SH)是水合物分解過程中沉積物的細觀黏結(jié)法向強度，單位MPa；SH0是沉積物的初始水合物飽和度，無量綱；SH是分解過程中水合物飽和度，無量綱；pb_ten(SH0)是水合物分解前的沉積物細觀黏結(jié)法向強度，單位MPa；pb_coh(SH)是水合物分解過程中沉積物的細觀黏聚力，單位MPa；pb_coh(SH0)是水合物分解前的沉積物細觀黏聚力，單位MPa。

圖12 流固耦合作用下的水合物開采出砂示意圖Fig. 12 Schematic diagram of sand production during gas recovery from hydrate-bearing sediments under fluid-solid coupling

3.3 流固耦合作用

如前所述，水合物開采降壓分解會產(chǎn)生氣體和水，水氣運移產(chǎn)生的流體作用力會影響出砂。因此，出砂主要是水合物分解效應(yīng)、地應(yīng)力和流體拖曳共同作用的結(jié)果(圖12)。有效應(yīng)力擠壓和滲流侵蝕會導(dǎo)致砂粒從地層骨架脫離并運移，而顆粒的運動則會改變地層孔滲特性從而間接影響流體運動，這種流固耦合作用是控制水合物儲層尤其是泥質(zhì)粉砂儲層出砂的重要因素[80]。

目前有若干將離散單元法與流體流動模型相結(jié)合的流固耦合建模方法，包括DEM-流體作用力、DEM-Darcy(達西定律)、DEM-CFD(計算流體動力學(xué))、DEM-LBM(格子玻爾茲曼方法)、DEM-SPH(光滑粒子流)、DEM-DNS(直接數(shù)值模擬)等[81-82]。其中DEM-LBM和DEM-CFD已被廣泛用于常規(guī)油氣儲層出砂機理及規(guī)律的研究。

Boutt等[83]、Han和Cundall[84]以及Ghassemi和Pak[85]均建立了二維DEM-LBM耦合數(shù)值模型進行出砂模擬研究。其中Boutt等定量分析了顆粒形狀和摩擦系數(shù)對于出砂的影響，發(fā)現(xiàn)累計出砂量會先增長然后趨于平穩(wěn)，較高的摩擦系數(shù)會抑制初期的出砂速率，但是穩(wěn)定階段的出砂量可能受顆粒形狀影響更大。砂拱的穩(wěn)定性問題一直是出砂機理研究的重點，砂拱對于阻止地層砂粒進入井筒進而抑制出砂具有重要作用。Han和Cundall以砂拱效應(yīng)為研究點，通過改變流體壓力梯度值再現(xiàn)了射孔區(qū)域砂拱形成、坍塌和重建的整個過程。Ghassemi和Pak的研究表明，地層出砂量取決于流體速率和地應(yīng)力。在給定地應(yīng)力下，當(dāng)流體速率超過臨界流速就會造成大量出砂，而增大地應(yīng)力能夠提高臨界流速。由于二維模型的幾何局限性，儲層的三維特征(比如孔隙度)無法得到真實反映，因此學(xué)者們致力于三維出砂數(shù)值模型的開發(fā)和應(yīng)用。

Zhou等[86]利用DEM-CFD耦合三維模擬方法，綜合考慮流體壓力梯度和拖曳力，研究了流體流動侵蝕弱膠結(jié)砂巖出砂以及形成空腔的過程，模擬結(jié)果還發(fā)現(xiàn)了和物理實驗一致的出砂團簇現(xiàn)象。但是Zhou等建立的三維出砂模型中只考慮了徑向流體流動而并未控制徑向邊界應(yīng)力，與實際情況略為不符。Climent等[87]同樣建立了DEM-CFD耦合的三維出砂模型，研究了遠場應(yīng)力與外邊界孔隙壓力(即流體壓力)對出砂的影響。但是該模型只考慮近井壁出砂引起的砂巖孔隙度的變化而并未關(guān)聯(lián)滲透率。Cui等[88]以及劉先珊和許明[89-90]在上述DEM-CFD耦合模型的基礎(chǔ)上，運用幾何學(xué)的方法，并考慮砂巖變形和出砂引起的砂巖滲透率的改變，研究了流固耦合作用下的出砂細觀力學(xué)機制。模擬結(jié)果均表明，流體運動對砂巖力學(xué)特性的影響不可忽略，它關(guān)系著砂巖屈服模式以及出砂量的大小。此外Cui等還詳細研究了邊界應(yīng)力對于出砂的影響，結(jié)果表明邊界應(yīng)力主導(dǎo)出砂響應(yīng)。與各向同性邊界應(yīng)力作用情況相比較，在不等邊界應(yīng)力作用下的地層出砂量更大。更多有關(guān)出砂流固耦合數(shù)值模擬研究見表1。

以上流固耦合作用研究的都是常規(guī)油氣儲層出砂問題，但適用于水合物儲層的流固耦合離散元模型還未見報道，原因或許是含相變的流固耦合作用過于復(fù)雜，建模較為困難?？紤]到氣體黏度低，對砂粒的拖曳效果弱[49,54]，后續(xù)可忽略氣體作用，以簡化基于PFC或者LIGGGHTS的水合物儲層流固耦合離散元數(shù)值建模和計算。

表1 出砂流固耦合數(shù)值模擬研究Table 1 Summary of the numerical simulation of fluid-soild coupling in sand production

4 討論與建議

目前水合物開采出砂數(shù)值模擬研究尚處起步階段，數(shù)值建模過程往往對一些因素進行假設(shè)和理想化設(shè)計，模擬結(jié)果也缺乏充分的實驗和現(xiàn)場驗證，導(dǎo)致數(shù)值模擬和分析的準(zhǔn)確性和適用性還不夠。如已有研究表明在實驗室模擬條件下水合物形成會排擠周圍顆粒產(chǎn)生“脹聚”效應(yīng)，但是在現(xiàn)場儲層條件下由于孔隙流體的存在，水合物分解產(chǎn)生的“縮散”效應(yīng)并不會讓周圍顆粒簡單復(fù)原或者復(fù)填。因此，如何準(zhǔn)確仿真孔隙尺度下水合物分解中和分解后骨架砂顆粒運移和相互作用是今后數(shù)值模擬需要解決的一個關(guān)鍵問題。未來建議在如下方面繼續(xù)加大研究力度：

(1)三維連續(xù)介質(zhì)模型構(gòu)建

連續(xù)介質(zhì)模型需要在井眼周圍劃分足夠精細的網(wǎng)格才能較準(zhǔn)確模擬水合物開采儲層響應(yīng)過程，這對三維模型的實現(xiàn)提出了很高的要求，巨大的網(wǎng)格數(shù)量必然消耗大量的計算資源和時間。因此，當(dāng)前連續(xù)介質(zhì)出砂模擬工作大多基于平面應(yīng)變假設(shè)或軸對稱假設(shè)從而構(gòu)建二維模型來完成。然而，二維模擬結(jié)果的可靠性較低。隨著后續(xù)數(shù)值計算能力的提升，依托測井解釋數(shù)據(jù)建立能夠刻畫儲層結(jié)構(gòu)特征如地層起伏、夾層狀態(tài)以及物性參數(shù)空間非均質(zhì)特性的三維連續(xù)介質(zhì)模型會是水合物儲層出砂模擬走向精細化的必經(jīng)之路。同時弱固結(jié)的水合物儲層變形和出砂也容易導(dǎo)致模型網(wǎng)格產(chǎn)生畸變，導(dǎo)致不收斂發(fā)生，后續(xù)或許可以考慮物質(zhì)點法等無網(wǎng)格方法進行水合物出砂數(shù)值模擬研究。

(2)滲透率模型修正

應(yīng)用離散元模型對水合物出砂流固耦合問題進行分析時，需要對流體流動方案做進一步修正。例如，在目前的DEM-CFD模型中，滲透率通常只與顆粒受力運移產(chǎn)出引起的孔隙度變化相關(guān)聯(lián)(參見下面的Kozeny-Carman方程，式(26))。然而，顆粒間膠結(jié)質(zhì)的退化，比如膠結(jié)型水合物逐步分解，同樣會影響儲層孔隙網(wǎng)絡(luò)結(jié)構(gòu)，進而改變儲層滲透率。也就是說，水合物降壓開采過程中儲層滲透率響應(yīng)行為，很大程度上取決于沉積物有效應(yīng)力與水合物飽和度共同控制的儲層孔隙結(jié)構(gòu)時空演化規(guī)律。因此，在離散元難以顯式刻畫水合物分解進程引起孔隙度變化的前提下，就需要對Kozeny-Carman方程進行修正，將式(26)中孔隙度φ與水合物飽和度相關(guān)聯(lián)以進一步描述水合物儲層出砂孔滲變化特征，建立綜合考慮水合物分解、出砂以及地應(yīng)力敏感的低滲水合物儲層滲透率模型。

圖13 水合物儲層顆粒間微力測試系統(tǒng)Fig. 13 Micro-force testing system between particles in hydrate reservoir

式(26)中：k為滲透率，單位m2，c為常數(shù)，無量綱；dp為顆粒的平均直徑，單位m；φ為孔隙度，無量綱。

(3)準(zhǔn)確高效的細觀參數(shù)標(biāo)定

有別于連續(xù)介質(zhì)模型對于復(fù)雜宏觀本構(gòu)關(guān)系的依賴性，離散單元法著重于細觀尺度研究，一組精確的細觀參數(shù)是該方法準(zhǔn)確揭示儲層出砂細觀機理的前提。泥砂顆粒間以及泥砂顆粒與水合物顆粒間的細觀作用力是厘清水合物開采出砂行為的關(guān)鍵。目前，獲取細觀參數(shù)的方法主要有兩種：直接測量法和試錯法[106]。直接測量法，顧名思義就是通過物理實驗直接測得顆粒尺度材料屬性，并且測量難度隨顆粒尺寸的減小而增大。試錯法就是借助數(shù)值力學(xué)實驗進行細觀參數(shù)標(biāo)定。細觀參數(shù)標(biāo)定是離散元模擬的主要難點之一，以顆粒流軟件PFC中內(nèi)置的平行膠結(jié)模型為例，其相關(guān)參數(shù)多達9 個[75]，一般的參數(shù)標(biāo)定過程就是通過不斷修改各參數(shù)值即試錯，直到獲得一組能夠使得數(shù)值實驗結(jié)果與材料宏觀性質(zhì)大體一致的細觀參數(shù)為止[107]。在各參數(shù)無明確量化關(guān)系的情況下，刪繁就簡，減少不確定參數(shù)的數(shù)量是簡化標(biāo)定的有效手段。中國地質(zhì)大學(xué)(武漢)聯(lián)合廣州海洋地質(zhì)調(diào)查局研發(fā)了一套水合物儲層顆粒間微力測試系統(tǒng)，見圖13。可以直接測得細觀顆粒尺度水合物抗拉強度和抗剪強度以及水合物分解后顆粒間的摩擦力，然后再通過試錯法標(biāo)定獲得其余參數(shù)值，這種直接測量法結(jié)合試錯法的綜合方法能夠更加準(zhǔn)確高效地為離散單元法的出砂數(shù)值模擬提供細觀參數(shù)。

(4)細粒的精細刻畫表征

我國南海水合物儲層主要為黏土質(zhì)粉砂，中值粒徑為10～20 μm[108]。較之日本砂質(zhì)儲層(中值粒徑約150 μm[109])，我國水合物儲層的細粒特點突出。并且細粒黏土常呈薄片狀，其表面帶負電，有雙電層結(jié)構(gòu)，這使得黏土顆粒間的相互作用與無黏性的砂顆粒相比要復(fù)雜的多。黏土的宏觀力學(xué)性質(zhì)與其微觀結(jié)構(gòu)密切相關(guān)，而離散單元法是研究巖土材料行為宏微觀聯(lián)系的重要工具[110]。此外，儲層中還富含有孔蟲殼體顆粒，因此利用離散單元法進行泥質(zhì)粉砂型水合物儲層出砂模擬，不但要考慮有孔蟲、黏土、水合物、砂粒的結(jié)構(gòu)特征差異和細觀力學(xué)參數(shù)差異，還要充分考慮有孔蟲、水合物、砂粒和黏土三者間的相互接觸關(guān)系。其中，黏土顆粒間的雙電層斥力和范德華力能夠使其在不發(fā)生機械接觸的情況下就產(chǎn)生碰撞響應(yīng)，需要發(fā)展新的接觸算法來表征這一現(xiàn)象。宏觀連續(xù)介質(zhì)模擬方面，細粒的影響需要通過本構(gòu)模型和強度準(zhǔn)則進行刻畫表征。目前，諸多基于Duncan-Chang本構(gòu)模型[111-112]、臨界狀態(tài)模型[113-114]和損傷本構(gòu)模型[115-122]的非線性彈性以及彈塑性含水合物沉積物本構(gòu)模型被提出，張峰等[123]對此進行了綜述，闡述了各模型的優(yōu)缺點，然而對于最優(yōu)本構(gòu)模型的選擇卻莫衷一是。此外，2.1 節(jié)中Jin等[42]研究結(jié)果也表明采用由不同儲層物性參數(shù)對應(yīng)建立的強度準(zhǔn)則進行水合物藏開采出砂風(fēng)險評估會得到差異較大的結(jié)果。水合物開采出砂問題歸根結(jié)底是力學(xué)失穩(wěn)所致，因此與研究對象相符的、能夠全面反映多因素條件下含水合物沉積物物理力學(xué)性質(zhì)的、經(jīng)過試采現(xiàn)場實測數(shù)據(jù)校核優(yōu)化的本構(gòu)模型和強度準(zhǔn)則亟待發(fā)展，這對于出砂問題的量化分析尤為重要。

(5)考慮顆粒破碎現(xiàn)象

水合物儲層降壓開采會導(dǎo)致地層有效應(yīng)力增大，從而可能引起有孔蟲或泥砂顆粒破碎[124]，見圖14。實驗結(jié)果[124]也表明，水合物的存在能夠阻礙砂粒的破壞，而水合物分解則會加劇砂粒破碎。砂粒破碎會使得地層收縮壓實、滲透性降低，并產(chǎn)生更細的顆粒，從而加劇出砂[125]。Zhang等[126-128]基于斷裂力學(xué)提出了一種關(guān)聯(lián)顆粒破碎程度和能量的破碎本構(gòu)模型，該模型可以反映可壓碎顆粒集料的力學(xué)響應(yīng)。連續(xù)介質(zhì)模擬方面可以借鑒Zhang等的研究成果，在上述破碎本構(gòu)模型的基礎(chǔ)上再添加水合物的影響，從而發(fā)展適用于水合物儲層的顆粒(包括有孔蟲顆粒、砂顆粒、水合物顆粒)破碎本構(gòu)模型。離散介質(zhì)方面，建立顆粒破碎模型一般有3 種方法：替換法、團簇法和組合法[129-130]。替換法思路為當(dāng)顆粒受力情況滿足預(yù)先定義的破碎準(zhǔn)則時，采用一組無黏結(jié)的小顆粒替換發(fā)生破碎的原始顆粒。團簇法通過黏結(jié)一定數(shù)量的較小顆粒形成團聚體，顆粒的破碎與否通過黏結(jié)強度來判斷。賀潔等[68]構(gòu)建水合物團簇的思想就屬于團簇法。復(fù)合法則是將置換法和團簇法相結(jié)合發(fā)展起來的，初始模型為單顆粒集合，模擬過程中篩選出高應(yīng)力顆粒并被團聚體所取代，而其余顆粒仍然是單顆粒。通過這種方法，可以用團聚體的破碎來精確地模擬顆粒的破碎，而團聚體的破碎是由黏結(jié)破壞決定的。該方法比置換法更準(zhǔn)確，比團簇法更高效[129]。因此可選用組合法來描述水合物儲層的顆粒破碎過程。

圖14 砂粒在應(yīng)力作用下發(fā)生破碎(a)砂粒破碎前(b)砂粒破碎后[125]Fig. 14 Particle crushing in hydrate-bearing sands under stress (a) before sand crushing (b) after sand crushing[125]

(6)分階段建模

根據(jù)水合物的分解狀態(tài)，水合物降壓開采過程中井周儲層分別處于未分解、分解中、完全分解3 個階段。降壓初期，開采壓差不高，水合物未分解，但是存在應(yīng)力變化和流體流動，此時可能僅有少量砂粒從井周脫離，不構(gòu)成出砂風(fēng)險?？刹捎秒x散單元法預(yù)測此時的出砂情況，重在揭示不同水合物分布模式下的儲層砂粒的啟動條件(臨界流速和有效應(yīng)力)；當(dāng)開采壓差增大，水合物脫離相平衡條件開始分解。這一階段多場耦合作用最為復(fù)雜，可采用連續(xù)介質(zhì)和離散介質(zhì)耦合方法進行研究(圖15)，離散介質(zhì)區(qū)域用于厘清水合物在出砂機理中的角色并獲得瞬時出砂量。此外，離散介質(zhì)區(qū)域?qū)崟r獲得的儲層孔滲參數(shù)可直接傳遞給連續(xù)介質(zhì)模型，然后連續(xù)介質(zhì)區(qū)域結(jié)合出砂準(zhǔn)則預(yù)測整體出砂量；第3 階段是近井壁區(qū)域內(nèi)水合物完全分解后，此時與第一階段的區(qū)別主要是儲層物性大大改變，原先被水合物占據(jù)的空間被重新分配，導(dǎo)致滲透率、孔隙度和力學(xué)強度(殘余強度)變化很大。這一階段的出砂行為極有可能引發(fā)井周儲層掏空進而破壞開采井井身結(jié)構(gòu)甚至導(dǎo)致開采井報廢，不容忽視。發(fā)展新的本構(gòu)關(guān)系刻畫這一過程并結(jié)合連續(xù)介質(zhì)模擬方法進行研究是第三階段的重點。此外，上述3 個階段后續(xù)研究還要重點關(guān)注出砂與防砂工藝的整體建模與分析，如基于CT掃描的篩網(wǎng)建模和基于LIGGGHTS的離散元擋砂模擬。

(7)機器學(xué)習(xí)

水合物開采出砂是復(fù)雜的強耦合非線性問題，使得數(shù)值模型的建立以及數(shù)值求解尤為復(fù)雜且耗時。而時下蓬勃發(fā)展的機器學(xué)習(xí)方法兼具解決復(fù)雜非線性問題、擅長挖掘數(shù)據(jù)間內(nèi)在規(guī)律以及快速計算等優(yōu)點[131-132]，因此很適合用來解決水合物開采出砂問題。大數(shù)據(jù)是機器學(xué)習(xí)的基礎(chǔ)，然而當(dāng)前水合物實際場地試采數(shù)據(jù)較少，因此可利用已發(fā)展的水合物開采出砂數(shù)值模型進行不同參數(shù)及工況下的數(shù)值模擬以獲取建模所需的訓(xùn)練樣本并訓(xùn)練得到對應(yīng)的初始機器學(xué)習(xí)模型，再利用少量的試采數(shù)據(jù)樣本和較多的實驗數(shù)據(jù)樣本對初始機器學(xué)習(xí)模型進行修正。反過來，訓(xùn)練好的的機器學(xué)習(xí)模型還能用于數(shù)值模型的優(yōu)化。數(shù)據(jù)驅(qū)動的機器學(xué)習(xí)與模型驅(qū)動的數(shù)值模擬可以相互促進，相輔相成，將加快水合物出砂問題的攻克進程。

圖15 連續(xù)—離散介質(zhì)耦合的水合物開采出砂數(shù)值模型示意圖Fig. 15 Schematic diagram of numerical model of sand production during gas recovery from hydrate-bearing sediments based on the coupling of continuous and discrete medium

5 總結(jié)

出砂問題是掣肘水合物產(chǎn)業(yè)化進程的難題之一。要厘清水合物儲層出砂規(guī)律、揭示其機理并對出砂進行防控離不開準(zhǔn)確的出砂分析和預(yù)測數(shù)值模型?；谶B續(xù)介質(zhì)和離散單元法的數(shù)值模擬技術(shù)是當(dāng)前出砂研究的主要手段之一。在水合物多場耦合模型研究的基礎(chǔ)上引入Mohr-Coulomb強度準(zhǔn)則、侵蝕準(zhǔn)則、等效塑性應(yīng)變準(zhǔn)則和力平衡準(zhǔn)則進行出砂判斷，能夠較系統(tǒng)厘清出砂規(guī)律并從工程上預(yù)測出砂量。不同于連續(xù)介質(zhì)方法所遵循的宏觀唯象研究思路，離散單元法從顆粒材料間內(nèi)部相互作用力出發(fā)，可以捕捉每個砂粒在出砂進程中的運動和協(xié)作以及破壞機制，利于從細觀尺度闡明出砂機理。但由于水合物賦存模式多樣、分解效應(yīng)復(fù)雜、流固耦合作用交錯導(dǎo)致相應(yīng)的離散元出砂模型鮮有，相關(guān)研究有待推進。相信隨著現(xiàn)場試采與物理實驗樣本的增多以及對水合物儲層出砂現(xiàn)象認識的深入和數(shù)值技術(shù)的進步，水合物開采儲層出砂數(shù)值模擬會更加精確更加系統(tǒng)。建議今后從三維連續(xù)介質(zhì)模型構(gòu)建、滲透率模型修正、準(zhǔn)確高效細觀參數(shù)標(biāo)定、細粒的精細刻畫表征、考慮顆粒破碎與抗轉(zhuǎn)動作用、分階段建模以及機器學(xué)習(xí)等多個方面進一步深化當(dāng)前的水合物儲層出砂數(shù)值模擬研究。此外，儲層出砂是儲層工程地質(zhì)特性與人工鉆采活動共同作用的結(jié)果，出砂數(shù)值模擬研究除了關(guān)注儲層本身外，還應(yīng)充分結(jié)合防砂方式、生產(chǎn)規(guī)程等進行整體研究，實現(xiàn)出砂—產(chǎn)能—儲層穩(wěn)定一體化協(xié)同評價和預(yù)測，以更加科學(xué)的理論分析指導(dǎo)我國海域天然氣水合物開采防砂控泥完井設(shè)計，推動水合物產(chǎn)業(yè)化進程。