趙振偉，尚新春

(1.北京科技大學(xué)土木與環(huán)境工程學(xué)院，北京 100083；2.北京科技大學(xué)應(yīng)用科學(xué)學(xué)院，北京 100083)

天然氣水合物是由天然氣和水在一定溫度和壓力條件下形成的一種冰狀晶體，其廣泛存在于永久凍土帶和海底沉積物孔隙當(dāng)中。天然氣水合物的含碳量大約是地球上其他所有常規(guī)能源含碳量的兩倍，被認(rèn)為是未來一種理想能源[1-2]。天然氣水合的分解，主要受壓力和溫度控制，研究水合物開采過程中溫度場(chǎng)和壓力場(chǎng)的分布，能夠在實(shí)際生產(chǎn)中確定水合物的分解范圍，預(yù)測(cè)天然氣的產(chǎn)量。通過對(duì)天然氣水合物分解面位置的分析，可以找出影響分解面移動(dòng)速度的因素，為提高天然氣的產(chǎn)出效率提供可行的方法。對(duì)天然氣水合物開采過程中溫度場(chǎng)和壓力場(chǎng)的分析，是近20年多來的研究熱點(diǎn)。

作者最近研究了水合物開采過程中的壓力場(chǎng)和溫度場(chǎng)，在求解壓力場(chǎng)時(shí)只考慮了水的滲流，于是壓力場(chǎng)控制方程為線性方程[9]。本文在文獻(xiàn)[1]的基礎(chǔ)上，進(jìn)一步考慮了水合物相變潛熱的影響，在分解面上加入能量守恒方程，對(duì)水合物開采過程中溫壓場(chǎng)的軸對(duì)稱模型進(jìn)行了改進(jìn)，并通過波爾茲曼變換方法，得到溫度場(chǎng)壓力場(chǎng)耦合模型的半解析-半數(shù)值解。在確定分解面位置時(shí)，將分解面上的質(zhì)量守恒方程轉(zhuǎn)化為一常微分方程來求解。與文獻(xiàn)[2]不同的是，本文考慮了氣體的流動(dòng)，壓力場(chǎng)為非線性方程，并在分解面上考慮了氣體密度隨壓力和溫度的變化。

1 數(shù)學(xué)模型

考慮水合物區(qū)一豎直井，通過降壓法開采水合物。假設(shè)壓力場(chǎng)和溫度場(chǎng)以井筒垂直軸線為中心成軸對(duì)稱分布且與深度無關(guān)。引入柱坐標(biāo)系(r,θ,z)。本文將水合物的開采過程看作是一個(gè)移動(dòng)邊界問題，在水合物開采過程中，礦區(qū)被分為兩個(gè)區(qū)域：分解區(qū)r0≤r≤R和未分解區(qū)R≤r≤+∞，r=R(t)表示移動(dòng)邊界的位置。記Pn(r,t)為壓力場(chǎng)和Tn(r,t)為溫度場(chǎng)。本文各個(gè)力學(xué)量的下標(biāo)n=1、2，分別代表分解區(qū)和未分解區(qū)。

Goodrz Ahmadi(2007)給出了如下壓力場(chǎng)的非線性控制方程[3]：

(1)

(2)

因沉積物的導(dǎo)熱率較低，這里忽略了熱傳導(dǎo)效應(yīng)，僅考慮對(duì)流效應(yīng)，其溫度場(chǎng)控制方程為[3]：

(3)

式中，常數(shù)δ和η分別為孔隙中氣體的節(jié)流系數(shù)和絕熱系數(shù)；cv和c分別為氣體和沉積物的熱容。

假設(shè)在井口以及無窮遠(yuǎn)處壓力和溫度均為常數(shù)，則邊界條件如下：

(4)

(5)

并且初始條件如下：

(6)

(7)

在分解面上水合物因分解而產(chǎn)生氣體，其質(zhì)量守恒方程為[5]：

(8)

其中氣體體密度滿足如下關(guān)系[8]：

(9)

T0=273.15K；P0=1.01×105Pa；ρ0是氣體在T0、P0條件下的密度；z=0.88，是氣體的壓縮系數(shù)；TD、PD是水合物相變時(shí)的溫度和壓力。

假設(shè)氣體流速滿足滲流的達(dá)西定律：

(10)

將式(9)和式(10)帶入質(zhì)量守恒方程可得到：

(11)

其次，在分解面上，由于水合物分解會(huì)吸收很大部分熱量，所以在分解面上的能量平衡方程時(shí)，應(yīng)考慮水合物相變潛熱的影響。而沉積物的導(dǎo)熱率相對(duì)較低，在此忽略熱傳導(dǎo)的影響，可以推導(dǎo)出考慮水合物相變影響的的能量守恒方程如下：

(12)

式中：L是水合物的相變潛熱；ρh是水合物的密度。

將式(9)和式(10)代入式(12)可得：

(13)

最后，在分解面上，水合物的溫度壓力滿足水合物的相平衡條件[2]：

(14)

式中：a=0.0342K-1；b=0.0005K-2；c=6.4804。

在分解面上需滿足壓力和溫度的連續(xù)性條件：

(15)

(16)

至此，確定水合物分解過程的壓力場(chǎng)Pn(r,t)和溫度場(chǎng)Tn(r,t)(n=1,2)在數(shù)學(xué)上被歸結(jié)為：在分解區(qū)r0≤r≤R(t)和未分解區(qū)R(t)≤r≤+∞上, 求解偏微分方程組(2)和(3)的解, 使其滿足初-邊值條件(4)～(7), 同時(shí)滿足移動(dòng)邊界R=Rc(t)上的質(zhì)量和能量守恒方程(11)和(13) 以及連續(xù)性條件(15)～(16)。

2 問題的求解

引入波爾茲曼變換[4]

(17)

將式(17)代入方程(2)可得到壓力場(chǎng)的常微分方程表達(dá)形式：

(18)

(19)

利用式(17)，邊界條件(4)、初始條件(6)和界面連續(xù)性條件(15)可轉(zhuǎn)換為：

(20)

(21)

(22)

對(duì)方程(18)、(19)進(jìn)行兩次積分，再利用邊界條件(20)～(22)，可得到壓力場(chǎng)的解析表達(dá)式：

(23)

(24)

同樣，將式(17)以及壓力場(chǎng)解析表達(dá)式(23)、(24)代入方程(3)，可得溫度場(chǎng)的常微分方程表達(dá)形式：

(25)

(26)

利用式(17)，邊界條件(5)初始條件(7)及界面連續(xù)性條件(16)可轉(zhuǎn)換為：

(27)

(28)

(29)

結(jié)合邊界條件(27)～(29)，求解常微分方程(25)、(26)，可得到溫度場(chǎng)的解析表達(dá)式：

(30)

(31)

上述壓力場(chǎng)和溫度場(chǎng)的表達(dá)式中，PD、TD及與分解面位置有關(guān)的uc為不確定量，可通過界面質(zhì)量守恒方程(11)和能量守恒方程(13)及水合物的相平衡方程(14)來確定。在確定上述三個(gè)不定量時(shí)，首先將方程(11)和(13)相減可得到：

(32)

將方程(14)代入到式(32)，可得關(guān)于TD的非線性方程，利用數(shù)值方法求解出TD，再由方程(14)確定出PD。將壓力場(chǎng)Pn的表達(dá)式(23)、(24)代入能量守恒方程(13)，得到如下關(guān)于Rc(t)常微分方程的初值問題：

(33)

Rc(t0)=r0

通過數(shù)值方法求解上述常微分方程的初值問題(33)，即可確定t時(shí)刻的分解面位置Rc(t)。

3 數(shù)值算例

下面給出溫度場(chǎng)和壓力場(chǎng)分布以及分解面位置的數(shù)值計(jì)算結(jié)果，計(jì)算參數(shù)如下：

c=2000 J/kg,cv=3000 J/kg,r0=0.1 m,

α=0.08,β=0.1,L=43500 J/kg,T0=273.15K,ε=0.129,

P0=1.01×105Pa,ρ0=0.706 kg/m3,Ρh=910 kg/m3,δ=8.0×10-7K/Pa,η=3.2×10-6K/Pa,K1=5.2md,K2=0.4md,Φ= 0.2。

圖1給出了井筒壓力為4MPa、分解時(shí)間分別為30d、60d和90d時(shí)，壓力和溫度沿水平徑向的分布。由圖1中可以看出，由于水合物開采時(shí)井筒壓力急劇降低，井筒處的壓力變化梯度較大；在未分解區(qū)，分解面附近的壓力最小，向無窮遠(yuǎn)處逐漸增大到地層的初始?jí)毫?，并且壓力梯度越來越小。在同一水平位置，壓力隨著時(shí)間逐漸降低，但是對(duì)時(shí)間變化率卻在減小。由圖2可以看出，由于這里考慮了水合物相變潛熱的影響，水合物分解會(huì)吸收部分熱量，溫度的最低值在分解面附近，向兩邊逐漸過渡到地層的初始溫度，分解面附近溫度梯度較大。

圖1 不同時(shí)間壓力分布

圖2 不同時(shí)間溫度分布

圖3給出了不同的井筒壓力條件下，水合物的分解面位置距井筒中心的水平距離與時(shí)間的關(guān)系。從圖中可以看出，在井筒壓力一定的條件下，分解距離隨時(shí)間逐步增大，但對(duì)時(shí)間的變化率卻在不斷減小。在開始時(shí)刻，分解前緣移動(dòng)速度較快，隨著分解時(shí)間的進(jìn)行，分解前緣的移動(dòng)速度會(huì)逐漸變慢。通過對(duì)不同井筒壓力下分解距離的比較可以看出，井筒壓力越低，分解速度越快。當(dāng)時(shí)間為90d時(shí),井筒壓力從5MPa降低到3MPa，分解距離從7m增加到24m。由此可見，水合物的分解速度對(duì)井筒壓力比較敏感，所以在開采過程中，應(yīng)盡量降低井筒壓力以加快水合物的分解。

圖3 分解面距井筒中心的水平距離與時(shí)間的關(guān)系

通過壓力場(chǎng)的解析表達(dá)式，單位高度上，天然氣水合物開采時(shí)的產(chǎn)氣速率，可通過井筒處氣體的流速確定：

(34)

通過氣體產(chǎn)生速率的表達(dá)式(34)，可以得到降壓法開采水合物時(shí)產(chǎn)氣速率與時(shí)間的關(guān)系。圖4給出了井筒壓力為3MPa、4MPa、5MPa時(shí)，氣體產(chǎn)生速率隨時(shí)間的變化曲線。從圖4中可以看出，產(chǎn)氣速率隨時(shí)間的增加在減小，開始時(shí)減小幅度較大，但是經(jīng)過一段時(shí)間，產(chǎn)氣速率逐步趨于一穩(wěn)定值。隨著井口壓力的降低，氣體產(chǎn)氣速率會(huì)明顯提高。

圖4 產(chǎn)氣速率與時(shí)間的關(guān)系

4 結(jié) 論

(1)通過水合物分解過程中的壓力場(chǎng)和溫度場(chǎng)分布情況看，在井口處壓力梯度最大，壓力值最??；隨著徑向坐標(biāo)的增大，壓力逐漸變大并趨于地層的初始?jí)毫?；受水合物分解的影響，分解面附近溫度較低，溫度梯度較大。

(2)水合物開始分解時(shí)，天然氣的產(chǎn)氣速率隨時(shí)間在逐步遞減，當(dāng)經(jīng)過一段時(shí)間后，產(chǎn)氣率趨于一穩(wěn)定值。

(3)天然氣的產(chǎn)氣速率受井口壓力的影響明顯，通過降低井口的壓力，可以明顯提高天然氣的產(chǎn)出速率。在天然氣水合物的實(shí)際開采中，可以通過降低井口壓力來提高生產(chǎn)效率。

[1] Ahmadi G, Ji C, Smith DH. Production of natural gas from methane hydrate by a constant down-hole pressure well[J]. Energy Conversion and Management，2007，48(7)： 2053-2068.

[2] Zhenwei Zhao，Xinchun Shang. Analysis for temperature and pressure fields in process of hydrate dissociation by depressurization[J]. International Journal for Numerical and Analytical Methods in Geomechanics，2010(3).

[3] Ji C, Ahmadi G, Smith DH. Natural gas production from hydrate decomposition by depressurization [J]. Chemical Engineering Science，2001，56(20)：5801-5814.

[4] M.Janz. Moisture diffusivities evaluated at high moisture levels from a series of water absorption tests [J]. Materials and Structures，2007(35)：141-148.