天然氣水合物數值模擬方法及其應用

郭朝斌，張可霓，凌璐璐

（1. 北京師范大學水科學研究院，北京 100875；2. 地下水污染控制與修復教育部工程研究中心，北京 100875）

天然氣水合物數值模擬方法及其應用

郭朝斌1,2，張可霓1,2，凌璐璐1,2

（1. 北京師范大學水科學研究院，北京 100875；2. 地下水污染控制與修復教育部工程研究中心，北京 100875）

在開采天然氣水合物前有必要用數值模擬方法對其形成或分解的規律進行研究以及對開采方案進行評估。數值模擬可根據野外和室內實驗數據，定量計算壓力、溫度和天然氣水合物的動力學轉換關系，評價天然氣水合物的開采潛力，從而預測復雜系統的熱物理動力過程。TOUGH+HYDRATE因其能夠模擬水合物分解的多組分多相流過程，較其他軟件相對成熟而得到廣泛應用。概述其模擬的基本原理和方法，總結利用該軟件模擬包括降壓、注熱以及降壓和注熱相結合等方法的應用實例，結合我國天然氣水合物調查與研究現狀提出數值模擬方面的工作建議。

天然氣水合物；多組分多相流；數值模擬；TOUGH+HYDRATE軟件

天然氣水合物是由天然氣與水在高壓低溫條件下形成的類冰狀結晶物質，因其外觀似冰且遇火即可燃燒，故又稱為“可燃冰”、“氣冰”、“固體瓦斯”。在標準狀態下，一單位體積天然氣水合物的分解最多可產生164單位體積的甲烷氣體，因而其是一種重要的未來潛在能源。天然氣水合物是20世紀科學考察中發現的一種新的礦產資源，主要分布在大陸永久凍土、島嶼的斜坡地帶、活動和被動大陸邊緣的隆起處、極地大陸架以及海洋和一些內陸湖的深水環境。

有關學者估計，全球天然氣水合物中的碳含量是當前已探明儲量的所有化石燃料（包括石油、天然氣和煤）中碳總和的兩倍。作為高效、儲量巨大的新型能源，其廣闊的開發前景引起世界上很多國家的關注和重視，相繼開展勘查與研究。美、日等國已在各自海域發現并開采出天然氣水合物。自1999年開始，中國先后在南海、東海、陸上凍土區和國際海底區域開展天然氣水合物的調查研究，陸續發現了一系列地質、地球物理和地球化學異常標志，并在南海北部神狐地區、祁連山凍土區等地成功獲取水合物實物樣品，取得了找礦工作的重大突破[1]。據測算，中國南海天然氣水合物的資源量為700億噸油當量，約相當于中國陸上石油與天然氣資源量總數的二分之一。而在我國祁連山凍土區鉆探獲取天然氣水合物，是世界中低緯度高山凍土區首次發現，具有重大的科研、經濟與戰略意義[2]。

目前開采天然氣水合物的方法主要有降壓法、注熱法、化學試劑法以及這些方法的綜合。數值模擬可以充分利用野外和室內實驗數據，定量計算壓力、溫度和天然氣水合物的動力學轉換關系，評價天然氣水合物的開采潛力，從而預測復雜系統的熱—物理動力過程[3,4]。數值模擬有助于全面認識開發過程中各因素的敏感性，對水合物的產氣量進行定量估算，可為開采方案的確定提供科學依據，是對天然氣水合物開采方案進行評價的一種重要手段。

本文概要闡述天然氣水合物數值模擬的原理與方法，并結合我國在南海等地區開展的相關調查研究，介紹其實際應用。

1 數值模擬的研發及其特點

隨著世界各國對天然氣水合物日益重視，天然氣水合物開采模擬研究得到了快速發展。目前國際上比較知名的天然氣水合物模擬器包括：MH21-HYDRES、STOMP-HYD、CMG-STARS、HydrResSim、TOUGH+HYDRATE等。

MH21-HYDRES是日本企業、政府和學術界聯合對甲烷水合物的研究項目MH21第一階段研發的大型模擬器，可以模擬熱激發法、降壓法以及各種方法相互結合的方法開采甲烷的過程。該模擬器可適用于三維笛卡兒坐標系或者二維徑向坐標系下動態網格加密、4種組分（甲烷、水、甲醇和鹽分）、5種相態（氣相、液相、冰相、水合物相、鹽）的狀況。為減少計算量，該模擬器的網格剖分一般比較粗，由此也容易產生數值誤差，這是該模擬器存在的不足。一般通過在局部重點區域應用動態網格加密方法，來減少計算量和保證計算精度。

STOMP-HYD是由西北太平洋國家實驗室（PNNL）在1996年由美國能源部支持研發而成，用來模擬地下多相流體的運動和變化規律。其可模擬平衡模型和動力學模型下的水合物的形成和分解，主要用到4個質量守恒方程和1個能量守恒方程，水合物、冰、析出鹽和客體介質等被當作非移動相[5]。

CMG-STARS由加拿大CMG公司(Computer Modelling Group Ltd)研發，可用于三組流動、多組分流體的流動模擬。

TOUGH+HYDRATE是由美國伯克利國家實驗室（LBNL）研發的TOUGH+家族中模擬天然氣水合物的軟件，其能模擬多組分（包括添加劑）以及多種相態，并不斷改進和完善，使其在天然氣水合物模擬方面具有較強優勢，應用也更廣。

1998年，美國勞倫斯伯克利國家實驗室的Moridis在通用多相滲流模擬計算軟件TOUGH2的基礎上，發展了模擬水合物的模塊EOSHYDR，可以對凍土層和海洋地層中平衡狀態下的游離的甲烷水合物的進行簡單模擬。

2003年，Moridis對EOSHYDR進行改善和增強，發布了同樣基于TOUGH2的EOSHYDR2模塊，可模擬平衡和動力學狀態下的二元水合物。

2005年，Moridis進一步改進模型，發展為TOUGHFx/HYDRATE軟件，并提供了水合物形成和分解的平衡和動力學模型。模型中考慮四相（氣相、液相、冰相、水合物相）、9組分（水合物、水、天然甲烷、水合物分解出來的甲烷、第二種碳氫組分、鹽、水溶性抑制劑與組分熱焓等），各組分存在于各相中。該模型可以剖析在降壓、注熱、加入抑制劑等條件下水合物分解的機理，可描述在非絕熱條件下水合物分解和形成的過程。其能夠比較準確地刻畫溫度、壓力等多重參數及其變化下，水合物在地層中的運動過程。

2008年發展成TOUGH+HYDRATE，其是繼TOUGH2之后的新一代通用模擬器，用于對含有水合物的地質環境系統的模擬。在天然氣水合物沉積層中，如永久凍土層和深海沉積層，只要滿足達西定律，通過對質能平衡方程耦合的求解，TOUGH+HYDRATE能模擬非等溫條件下水合物分解、相態的變化、流體和熱流的流動情況[6]。

TOUGH+HYDRATE包含了兩個模型，分別是對水合物的形成分解的平衡模型和動力學模型。模型可對水、天然氣、天然氣水合物、可溶性的抑制劑等四個質量組分，以及氣相、液相、冰相和水合物相等四種可能的相態的分析與模擬。其采用標準FORTRAN95編寫而成，新的代碼采用面向對象的編程技術，并利用動態內存分布，有效地減少了空間的需求。

由張可霓等開發的pT+H軟件是TOUGH+HYDRATE的并行版本，可以模擬包括水合物、水、甲烷、鹽、添加劑、熱等涉及26種可能相態的多組分運移[7]。pT+H軟件既可以模擬平衡和動態分解過程，又可以模擬可能的降壓、熱、添加劑以及混合情況的分解過程，可在超大規模計算機系統中進行并行計算，是模擬器的理想選擇。在pT+H中利用域分解技術和MPI(message passing interface)實現并行計算。利用METIS軟件中的劃分算法，將模擬區域劃分成子區域。在多核或多CPU電腦中，每一個核/CPU處理一個劃分出來的子區域，將一個大區域劃分成子區域并由多個核/CPU進行熱動力學變量的計算，從而極大地提高了計算效率。另外，pT+H中用Aztec并行線性求解器求解局部線性方程。并行計算的高效性和靈活性，使pT+H能夠模擬大規模的實際場地。

下文將主要基于TOUGH+HYDRATE，介紹數值模擬的基本原理、方法與應用實踐。

2 數值模擬基本原理與方法

2.1 組分和相態

數值模擬可分析模擬的組分包括水合物(h)、水(w)、甲烷(m)、水溶性抑制劑(i)、熱( )等，所用的天然氣水合物的相態圖如圖1所示[6]。在相圖中，當溫度為273.2K、壓力達到2.3MPa時，會出現四相點（臨界點），即圖中的Q1點。

圖1 天然氣水合物相態圖Fig.1 Phase diagram of natural gas hydrate

在數值模擬中，通過幾個主變量來描述系統的熱動力學狀態，其它的次變量通過與主變量的關系進行求解得出。T+H中共有26種相及相組合，其中13種為水合物平衡反應的相態，13種和動力學反應相關。對于有抑制劑的系統，則增加關于抑制劑的主變量即抑制劑在液相中的質量分數。無抑制劑系統在各相態組合情況下的主變量如表1所示。

表1 平衡狀態下天然氣水合物模擬中的主變量Table 1 Primary variables in equilibrium hydrate simulations without inhibitor

2.2 控制方程和輔助方程

將模擬區域通過積分有限差分法進行離散，在通過積分有限差分法分割的每個網格中都滿足質量和能量守恒方程。

其中，積分區域Vn為所考察區域的任意的一個子區域，該子區域以為邊界。式左中的Mk表示單位體積的質量，其上標k視不同對象而定，對水合物成分k即h，水則為w，甲烷m，水溶性抑制劑i，熱組分θ。F為通量項，q為源匯項。

數值模擬中包含平衡模型和動力學模型，A、G、H、I分別代表液相、氣相、水合物相和冰相。在平衡模型中，

式中：?為孔隙度，ρβ為相β的密度，Sβ為相β的飽和度，Xβ

k為相β中各組分質量分數。

在動力學模型中下，

通量項；

源匯項；

除滿足主要控制方程外，還需要滿足一些輔助方程：

流體飽和度總和為1：

與β相中的k成分質量分數之和為1：

在多相流體的模擬中，由于流體的物理性質隨著熱力學參數（溫度、壓力、組分質量分數）而改變，相對滲透系數與毛細壓力是飽和度的函數。在數值模擬中，常用到9種毛細壓力—飽和度的關系曲線，包含線性關系曲線、Pickens模型、Narasimhan模型、Milly模型、Udell and Fitch模型、Van Genuchten模型和Brooks-Corey模型等。經常采用的Van Genuchten模型為：

式中，Pcap表示液相毛細壓力，滿足-Pmax≤Pcap≤0的限制條件；λ為模型擬合參數。

式中，S*表示有效水飽和度、S1r表示殘余水飽和度、S1s表示飽和水飽和度。

對于相對滲透率和飽和度之間的關系，也有9種關系模型，如線性模型、Pickens模型、Corey模型、Grant模型、Van Genuchten-Mualem模型等。其中Van Genuchten-Mualem模型為：

式中，kr1為液相相對滲透率、krg為氣相相對滲透率。

式中，Sgr為殘余氣飽和度。

2.3 空間和時間離散

數值模擬中需要對模擬區域進行離散，通常采用積分有限差分法進行空間離散，其可將模擬區域離散成任意形狀的多面體。在計算過程中，只需獲得單元的體積、相鄰兩網格的接觸面積與距離以及它們之間連線與垂向的交角等相關參數。在處理任意形狀的單元時不必考慮坐標系統，同時也不受單元塊鄰近單元數限制。

圖2中，Anm為表面，Fnm為組分F通過Anm表面的平均值。

圖2 積分差分法中空間離散和流動項計算示意Fig.2 Space discretization and geometry data in the integral finite difference method

時間采用一階有限差分離散。時間和空間的離散結果用殘差方程的形式表示為：

其中，xt表示t時刻的主變量，表示在區域n范圍內的k組分的殘差，Vn表示區域n的體積，M表示組分單位體積的質量或者熱能。求解非線性方程，采用Newton-Raphson迭代求解。

3 數值模擬技術應用實踐

數值模擬技術已在天然氣水合物的調查研究方面得到廣泛應用，發揮了重要作用。我國在此領域也有長足進步，成果豐碩。

有關學者通過數值模擬方法對南海神狐區域的產氣量潛力進行了評估[8,9]。吳能友等根據神狐海域有鉆孔確定的含水合物區的各種參數，認為在50%概率條件下，該區水合物資源量約為160×108m3[10]。

李小森等針對深海條件下天然氣水合物礦藏的單井降壓法開采，對開采過程和甲烷氣體擴散進行了數值模擬[11,12]。結果表明，開采中井口產氣速率是一個升高—降低—波動升高的過程，水合物分解產生的氣體有一部分通過上蓋層溢出，在一定程度上將增加大氣中溫室氣體含量；開采初期水合物分解速率降低的主要原因是由于水合物分解產生的甲烷氣體在地層中大量累積，后期分解速率產生波動則是因為發生了“氣穴現象”；在井口附近，因壓力變化較快，致使水合物分解最為劇烈。

胡立堂等針對南海神狐海域天然氣水合物特征，利用pT+H軟件和水平井技術對SH2站位天然氣水合物開采進行了模擬研究，討論了單純注熱和注熱與降壓結合兩種模式下的水合物開采效果[13]。結果認為，單純提高注水溫度對于CH4產氣量沒有明顯改善，因此開采時選擇適宜溫度即可；溫度太低時可能形成二次水合物，不利于后期開采，而溫度太高也不經濟；注熱降壓比單純注熱開采的效果有較大改善，可作為神狐海域天然氣水合物開發利用的優選方案。

張可霓等根據南海神狐海域SH2鉆探數據，利用多個水平井進行注熱和降壓相結合開采方案進行模擬研究[14]。網格剖分如圖3所示，圖4為注熱溫度50℃的方案中不同時間水合物飽和度的分布模擬結果。注熱對注入井附近水合物飽和度分布影響較大，早期的分解主要集中在注入井附近，逐漸形成圓柱狀飽和度分布。在經過近一年的產氣后，分解前緣接近產氣井。在產氣井附近的壓降同樣引起產氣井附近小范圍水合物的分解。在兩種方法引起水合物分解邊緣相遇時，產氣率達到峰值。由模擬結果可知，產氣井附近的小尺度壓降對水合物分解的影響范圍有限。

圖3 雙井方案二維網格剖分示意圖Fig.3 2D model mesh for the two-well gas production system

圖4 注熱開采與降壓開采綜合方案水合物飽和度分布模擬熱量通過下部井注入50℃的水Fig.4 Evaluation of hydrate saturation distribution during gas production with heat stimulation from the hydrate accumulations. Heat was injected by applying a constant temperature of 50°C at the lower well

數值模擬技術還可為鉆井施工方案提供支持。寧伏龍等利用數值模擬，分析了水基鉆井液侵入海洋含水合物地層過程及其對井周地層的影響[15]。據此初步掌握了在過平衡鉆井條件下，當鉆井液溫度高于地層水合物穩定溫度時，鉆井液侵入含水合物地層的基本特征。

4 結語

數值模擬是對天然氣水合物開采方案進行比選與評價的重要手段，能夠針對多組分多相流的模擬技術得到廣泛應用。隨著我國在天然氣水合物勘查方面取得重大突破，有必要在借鑒國外先進模擬技術的同時，研發具有自主知識產權的實用軟件和分析系統，以促進天然氣水合物的開發利用。

References)

[1] 張洪濤,張海啟,祝有海. 中國天然氣水合物調查研究現狀及其進展[J]. 中國地質,2007,34(6):953-961.

Zhang H T, Zhang H Q, Zhu Y H. Gas hydrate investigation and research in China: Present status and progress[J].Geology in China,2007,34(6):963-961.

[2] 祝有海,張永勤,文懷軍,等. 祁連山凍土區天然氣水合物及其基本特征[J]. 地球學報,2010,31(1):7-16.

Zhu Y H, Zhang Y Q, Wen H J, et al. Gas hydrates in the Qilian mountain permafrost and their basic characteristics [J].Acta Geoscientica Sinica,2010,31(1):7-16.

[3] 張郁,李小森,李剛,等. 天然氣水合物分解和開采的機理及數學模型研究綜述[J]. 現代地質,2010,24(5):979-985.

Zhang Y, Li X S, Li G. Review: Mathematic modeling of natural gas hydrate exploitation[J].Geoscience, 2010,24(5):979-985.

[4] 李淑霞,陳月明,杜慶軍,等. 天然氣水合物開采方法及數值模擬研究評述[J]. 中國石油大學學報(自然科學版),2006,(3):146-150.

Li S X, Chen Y M, Du Q J, et al. Commentary of production method and numerical simulation of natural gas hydrates[J].Journal of China University of Petroleum (Edition of Natural Science),2006,(3):146-150.

[5] White M D, Oostrom M. Subsurface transport over multiple phases application guide[R]. PNNL-11218. Pacific Northwest National Laboratory, Washington US, 1997.

[6] Moridis G J, Kowalsky M B, Pruess K. TOUGH+HYDRATE user’s manual: A code for the simulation of system behavior in hydrate-bearing geologic media[R]. Lawrence Berkeley National Laboratory, Berkeley US, 2008.

[7] Zhang K N, Moridis G J, Wu Y S, et al. A domain decomposition approach forlarge-scale simulations of flow processes in hydratebearing geologic media[A]. Proceedings of the 6th International Conference on Gas Hydrates(ICGH 2008)[C]. Vancouver Canada, 2008.

[8] Li G, Moridis G J, Zhang K N, et al. Evaluation of gas production potential from marine gas hydrate deposits in Shenhu area of South China Sea[J].Energy & Fuels, 2010, 24 (11):6018-6033.

[9] Su Z, Cao Y C, Wu N Y, et al. Numerical investigation on methane hydrate accumulation in Shenhu area, northern continental slope of South China Sea[J].Marine and Petroleum Geology, 2012,38(1):158-165.

[10] Wu N Y, Yang S X, Zhang H Q, et al. Gas hydrate system of Shenhu area, northern South China Sea: Wire-line logging, geochemical results and preliminary resources estimates[A]. Proceedings of the Offshore Technology Conference[C]. Texas US, 2010.

[11] 李小森,陳琦,李剛,等. 海底水合物礦藏降壓開采與甲烷氣體擴散過程的數值模擬[J]. 現代地質,2010, 24(3):598-606.

Li X S, Chen Q, Li G., et al. The simulation of hydrate production by depressurization in deep ocean[J].Geoscience, 2010,24(3):598-606.

[12] 李剛,李小森,陳琦,等. 南海神狐海域天然氣水合物開采數值模擬[J]. 化學學報,2010,(11):1083-1092.

Li G, Li X S, Chen Q, et al. Numerical simulation of gas production from gas hydrate zone in Shenhu area, South China Sea[J].Acta Chimica Sinica,2010,(11):1083-1092.

[13] 胡立堂,張可霓,高童. 南海神狐海域天然氣水合物注熱降壓開采數值模擬研究[J]. 現代地質,2011,25(4): 675-681.

Hu L T, Zhang K N, Gao T. Numerical studies of gas production from gas hydrate zone using heat injection and depressurization in Shenhu area, the South China Sea[J].Geoscience, 2011,25(4):675-681.

[14] Zhang K N, Moridis G J, Wu N Y, et al. Evaluation of alternative horizontal well designs for gas production from hydrate deposits in the Shenhu area, South China Sea[A]. International Oil and Gas Conference and Exhibition in China[C]. Beijing China, 2010.

[15] 寧伏龍,張可霓,吳能友,等. 鉆井液侵入海洋含水合物地層的一維數值模擬研究[J]. 地球物理學報, 2013,56(1):204-218.

Ning F L, Zhang K N, Wu N Y, et al. The invasion of waterbased drilling fluid into oceanic gas hydrate-bearing sediment: One-dimension numerical simulations[J].Chinese Journal of Geophysics, 2013,56(1):204-218.

Numerical Simulation Methods and their Application to Natural Gas Hydrate Exploration

GUO Chao-Bin1,2, ZHANG Ke-Ni1,2, LING Lu-Lu1,2
(1.College of Water Science, Beijing Normal University ,Beijing 100875, China;
2. Engineering Research Center for Groundwater Pollution Control and Remediation, Ministry of Education, Beijing 100875, China)

Before drilling for natural gas hydrate (NGH), it is necessary to use a numerical simulation method to study its formation and/or decomposition, and to identify a suitable scheme for its exploration. According to both field and laboratory experimental data, numerical simulations can be used to quantitatively calculate the dynamic transformation relationships among pressure, temperature, and NGH, and to evaluate the potential for NGH production by predicting the thermal physical dynamic processes involved in its occurrence. The TOUGH+HYDRATE software is used to simulate the decomposition process of NGH with respect to multicomponent and multiphase flow. The basic principles and methods used in the software are described, and an example is given of the application of simulations using the software, including depressurization, thermal stimulation, and the combination of both. The paper discusses current research and numerical simulation work being conducted with regard to investigations of NGH in China.

natural gas hydrate; multicomponent and multiphase flow; numerical simulation; TOUGH+HYDRATE software

P744.4

A

2095-1329(2013)02-0071-05

10.3969/j.issn.2095-1329.2013.02.017

2013-03-19

2013-05-04

郭朝斌(1989-),男,碩士生,主要從事多相流體數值模擬研究.

電子郵箱：cugbgcb@163.com

聯系電話：18901289515

國家能源應用技術及工程示范項目(NY20111102-1)