李明川，樊栓獅，徐賦海，嚴(yán)科，黃愛先

（1中國石油大學(xué)（華東）石油工程學(xué)院，山東青島266580；2華南理工大學(xué)化學(xué)與化工學(xué)院，廣東廣州510641；3中國石化勝利油田分公司東辛采油廠，山東東營257094；4中國石化勝利油田分公司勝利采油廠，山東東營257051）

引 言

天然氣水合物是由碳?xì)錃怏w和水形成具有籠型結(jié)構(gòu)的、類似冰狀的、可燃燒的晶體化合物，俗稱“可燃冰”[1-3]。天然氣水合物為高度壓縮的天然氣資源，1 m3天然氣水合物能分解出160～180 m3（標(biāo)準(zhǔn)）的天然氣[4]。天然氣水合物廣泛存在于海洋沉積物和永凍土中，據(jù)估計全球天然氣水合物中的甲烷碳含量達(dá)1016kg或含有20×1015m3的甲烷氣，因此天然氣水合物被公認(rèn)為是21世紀(jì)具有良好前景的后續(xù)能源[5]。

天然氣水合物的熱分解過程類似于冰的消融過程，是一個考慮移動邊界的熱傳導(dǎo)相變問題，也稱為兩相Stefan相變問題[6]。Stefan邊界問題的研究可以追溯到19世紀(jì)，Lamé等[7]在假設(shè)地球為半無限液體的條件下，研究了地球冷卻過程中產(chǎn)生的固體地殼的厚度與時間的平方根呈正比的類似問題。Neumann假定初始溫度高于熔點時，得出了半無限大Stefan問題的特解，簡稱Neumann解，他的成果直到20世紀(jì)初才由Weber出版發(fā)表出來[8]。最后解決這問題的是Stefan，他將問題轉(zhuǎn)化為兩個相的互變問題，利用Lamé等的假設(shè)條件，得到了更為完全的解。Stefan第一次從實驗到理論對移動邊界的相變問題進(jìn)行了詳細(xì)的研究，因此自由或移動邊界問題也稱為Stefan邊界問題[9]。Stefan移動邊界問題逐漸被廣泛應(yīng)用于固液相變、化學(xué)反應(yīng)、溶質(zhì)擴(kuò)散、液氣相變等諸多領(lǐng)域[10]。

天然氣水合物熱分解Stefan理論是一個新興的跨學(xué)科應(yīng)用領(lǐng)域，部分學(xué)者在這方面進(jìn)行了大量研究工作。Makogon[11]首次使用類似經(jīng)典Stefan的固體融化問題，描述了水合物的降壓分解過程，獲得了壓力分布的自相似解。Selim等[12-13]認(rèn)為水合物分解過程是一個移動邊界消融問題，利用Von-karman積分得到了溫度分布的近似解。Yousif等[14]建立了水合物分解界面的移動邊界模型，并用該模型重新評估了水合物分解的氣體體積。Tsypkin[15]建立了包含水合物分解前緣和冰融化前緣的兩個未知移動邊界的Stefan廣義數(shù)學(xué)模型。唐良廣等[16]將水合物分解過程看成界面運(yùn)動過程，獲得了熱傳導(dǎo)模型的解析解。李明川等[17]建立了水合物熱分解過程Stefan徑向模型，得到了Paterson指數(shù)積分函數(shù)精確解。Li等[18-19]提出了水合物熱分解準(zhǔn)穩(wěn)態(tài)Stefan相變模型，確定了水合物熱分解前緣位置數(shù)學(xué)模型。

綜上，前人雖給出了水合物熱分解Stefan相變模型，但未能詳細(xì)給出建模、求解及檢驗過程，這有礙于學(xué)科、行業(yè)的應(yīng)用性研究進(jìn)展。作者對水合物熱分解Stefan相變過程進(jìn)行了系統(tǒng)性的研究，包括連續(xù)熱傳導(dǎo)模型、不連續(xù)界面Rankine-Hugoniot跳躍條件和Stefn耦合界面條件等理論。本文在單相連續(xù)水合物控制體能量守恒積分模型基礎(chǔ)上，建立了考慮具有尖銳移動邊界分隔的水合物熱分解控制體的分解界面耦合Stefan能量守恒條件。給出了一半無限大天然氣水合物熱分解Stefan相變模型及定解條件，利用Boltzmann相似變量求得了誤差函數(shù)形式的Neumann解，并對界面耦合Stefan超越方程進(jìn)行了單調(diào)性證明，確定了Stefan模型解的唯一性，結(jié)合實例進(jìn)行了分析計算。通過本文的理論研究工作，可促進(jìn)水合物熱分解相變規(guī)律、相變機(jī)理的揭示，為推動天然氣水合物工業(yè)化開采提供理論依據(jù)。

1 控制體熱傳導(dǎo)Stefan模型

有一體積一定的連續(xù)的單相水合物控制體V，外表面積S，任取一微元dV和任一表面積dS，單位法向量n由控制體指向外部（圖1）[20]。若控制體無內(nèi)部熱源和外來能量做功，能量隨時間的變化率與通過表面積S傳導(dǎo)到控制體V中的熱量凈速率相等（無對流），其數(shù)學(xué)形式為：

圖1 連續(xù)水合物控制體Fig.1 Continuous hydrate control volume

式中，(ρε)為水合物單位體積的總能量，kJ/m3；q為通過單位表面積的熱傳導(dǎo)量，kJ/m2。總能量為內(nèi)能、動能及勢能之和，忽略勢能，則內(nèi)能為：

體積一定的水合物控制體V與時間t無關(guān)，熱量(ρh)連續(xù)地通過微元體dV，式（4）左端的體積分和對時間的微分運(yùn)算可互換，移動到積分式中，d/dt變?yōu)??t；q=k·gradT連續(xù)地通過微元表面dS，式（4）右端可根據(jù)散度定理[21]，

式（7）即是單相連續(xù)區(qū)域的熱傳導(dǎo)方程。

接下來將式（4）應(yīng)用于具有跳躍移動界面的水合物熱分解控制體（圖2）[22]。t時刻，界面Σ將控制體V分為液相分解區(qū)Vl和固相水合物區(qū)Vs，表面積分為Sl和Ss，單位法向量n l和ns分別指向各區(qū)域外方向。經(jīng)過一小段時間δt后，界面移動到Σ'位置，變化的體積為δV，意味著在這一小段時間內(nèi)有δV體積的固相水合物參與到熱分解過程，相應(yīng)的液相分解區(qū)增加δV的體積（忽略分解氣體的影響），單位法向量n*由液相分解區(qū)Vl指向固相水合物區(qū)Vs。在這個過程中，由界面Σ到界面Σ'通過的焓h和溫度梯度k·gradT是不連續(xù)的，但在各自的區(qū)域（液相分解區(qū)Vl和固相水合物區(qū)Vs）是連續(xù)的。為了推導(dǎo)界面的能量守恒條件，將式（4）左端分解到每個區(qū)域上分別積分，研究它們隨時間的變化。

圖2 界面不連續(xù)水合物控制體Fig.2 Interface discontinuous hydrate control volume

在時間t時刻:

注意到式（10）右端最后一極限項，表示在δt時間界面Σ→Σ'所增加的δV體積內(nèi)的熱焓值。當(dāng)δt→0時，δV/δt=v n*dΣ，其中v n*是界面微元dΣ運(yùn)動速度，方向垂直于界面Σ指向固相水合物；在δt→0時δV所在的空間趨近于表面Σ，水合物分解潛熱L=hl-hs，則式（10）的δV積分變?yōu)棣卜e分:

由式（4）分別代替式（11）中Vl、Vs上的積分，相應(yīng)液相分解區(qū)和固相水合物在δt→0時增加了分解界面面積Σ，則:

這里單位法向量n*在界面Σ上由液相分解區(qū)指向固相水合物區(qū)，兩相表面上的法向量有n l=n*，將式（12）右端按照兩相表面Sl、Ss和界面Σ分別積分，且S=Sl+Ss，但積分式中Sl與Ss方向相反，則有：

將式（13）重新整理:

對任意的界面Σ，被積函數(shù)在界面上的所有點均消失，即環(huán)繞界面Σ積分等于0[23]。對比式（14）和式（4）可得：

式(15)即為水合物熱分解的界面耦合Stefan能量守恒跳躍條件。

2 半無限大天然氣水合物熱分解Stefan模型及Neumann解

水合物分解界面條件：

定義水合物分解界面移動速度v n*=dxfdt，則液相分解區(qū)Vl和固相水合物區(qū)Vs的耦合Stefan條件由式（15）可得：

式（19）～式（24）構(gòu)成了半無限大水合物藏一維Stefan模型，該模型的解是關(guān)于液相分解區(qū)和固相水合物區(qū)溫度分布及未知前緣位置的移動邊界的函數(shù)，而且溫度與移動邊界有關(guān)。由此可得方程的解（見附錄）。

移動界面位置：

液相分解區(qū)溫度：

固相水合物區(qū)溫度：

式（26）和式（27）的溫度通過界面方程式（24）耦合一起，其中λ是超越方程式（28）的解。

對該超越方程進(jìn)行單調(diào)性證明，以確定其解具有唯一性[24-25]。

3 結(jié)果與討論

已知一半無限大一維純甲烷水合物儲藏，水合物物性如表1所示。

表1 水合物熱分解一維模型基礎(chǔ)數(shù)據(jù)［26-27］Table 1 One dimensional model basic data for thermal dissociation hydrate［26-27］

水合物在分解過程中保持恒定壓力，忽略注入熱水的影響，可通過Deaton相平衡關(guān)聯(lián)方程[28]計算水合物的熱分解溫度Tm：

通過計算在壓力5.55 MPa下，水合物的熱分解溫度為5.305089℃。理論計算表明，1個單位體積的可燃冰可以分解為164個單位體積的天然氣及0.8個單位體積的水[29-30]，根據(jù)Clausius-Clapeyron方程[29-30]可計算水合物摩爾分解熱L為14.83763 W/mol。

式中，ng為分解甲烷氣體的量，mol；z為甲烷氣體偏差因子；R為通用氣體參數(shù)，R=0.008314 MPa·m3/(kmol·K)。

從圖3可以看出，隨著λ增加超越方程急劇下降，隨后平緩變化，但總體上保持單調(diào)遞減，驗證了對式（28）的證明。在保持注入溫度Tinj=100℃時，超越方程唯一解為λ=1.5021[圖3(a)]；隨著注入溫度的增加（60，100，200℃），超越方程的解為λ=1.3515，1.5021，1.6888，對超越方程的影響逐漸增大，但增加幅度減小[圖3(b)]。

圖3 超越方程單調(diào)性及解關(guān)系曲線Fig.3 Relationship curves of the monotonicity and the solution of the transcendental equation

從圖4看出，液相水合物分解區(qū)溫度隨著距離的增加急劇下降至前緣xf的分解溫度Tm=5.305089℃；經(jīng)歷不連續(xù)的躍變后，溫度穿透進(jìn)入固相水合物區(qū)，下降幅度逐漸減小直至水合物藏溫度TR=3.50℃，并保持不變。保持注入溫度Tinj=100℃時，隨著時間增加水合物分解前緣位置逐漸增大，溫度穿透深度xd（定義：溫度穿越分解前緣所能到達(dá)的距離，即xd=x TR-xTm）逐漸增大；5、10和20 d的水合物分解前緣xf分別到達(dá)3.3568、4.7472和6.7135 m，穿透深度分別為4.8622-3.3568=1.5054 m、6.8933-4.7472=2.1461 m 和 9.7180-6.7135=3.0045 m，穿透深度增加幅度逐漸減小[圖4(a)]。隨著注入溫度的增加（60、100、200℃），液相水合物分解區(qū)溫度下降幅度增大，水合物分解前緣位置增加；在10 d時間，分解前緣分別到達(dá)4.2712、4.7472和5.3372 m；穿透深度逐漸降低，分別為6.5634-4.2712=2.2922 m、6.8933-4.7472=2.1461 m 和7.3266-5.3372=1.9894 m，這是因為固相水合物區(qū)溫度下降梯度增大，在相同時間10 d的穿透深度降低了[圖4(b)]。

圖4 水合物體溫度隨距離關(guān)系曲線Fig.4 Relationship curves of temperature of hydrate control body vs.distance

從圖5可以看出，隨著時間增加，固相水合物區(qū)溫度從水合物藏溫度TR=3.50℃緩慢上升，至分解溫度Tm=5.305089℃后溫度迅速增大，但增加幅度逐漸降低。保持注入溫度100℃時，達(dá)到分解前緣xf=6、8和10 m的時間分別為15.9745、28.3993和44.3738 d，溫度穿透水合物體后，從水合物藏溫度TR=3.50℃到分解溫度Tm=5.305089℃所需穿透時間td（定義：溫度穿越指定前緣所需要的時間，即td=tTR-tTm）分別為15.9745-7.6=8.3745 d、28.3993-13.6=14.7993 d和44.3738-21.2=23.1738 d，穿透時間隨前緣位置增大，增加幅度變大[圖5(a)]。保持注入溫度100℃到1000 d時，在分解前緣xf=6、8和10 m的溫度分別為79.2572、72.5967和66.1462℃。隨著注入溫度的增加（60，100，200℃），液相水合物分解區(qū)溫度逐漸增大，隨著時間增加上升幅度減??；不同注入溫度到達(dá)指定前緣位置xf=8 m所需時間分別為35.0814、28.3993和22.4674 d，穿透時間為35.0814-11.9=23.1814 d、28.3993-11=17.3993 d和22.4674-10=12.4674 d，穿透時間減小幅度降低[圖5(b)]。注入恒定60、100和200℃溫度在指定分解前緣xf=8 m到1000 d所能達(dá)到的溫度分別為43.7979、72.5967和144.6165℃，且增加幅度在擴(kuò)大。

圖5 水合物體溫度隨時間關(guān)系曲線Fig.5 Relationship curvesof temperature of hydrate control body vs.time

從圖6可以看出，隨著注入時間增加，水合物分解前緣逐漸增大，但增加幅度在減小。隨著注入溫度的增加（60、100、200℃），水合物分解前緣增加，10 d時分解前緣達(dá)到xf=4.2712、4.7412和5.3372 m，前緣增加幅度降低；隨著注入溫度的增加（60、100、200℃），到達(dá)分解前緣xf=8 m所需時間減小，分別為35.0814、28.3993和22.4674 d，減小幅度降低。

圖6 水合物分解前緣隨時間關(guān)系Fig.6 Relationship curves of the dissociation frontal brim of hydrate vs.time

接下來對溫度敏感性參數(shù)進(jìn)行分析。從圖7可以看出，溫度對λ和xf的影響逐漸增大，但增加幅度減??；溫度對xd和td的影響逐漸減小，減小幅度降低。擬合溫度與超越方程的解λ具有較好的二項式關(guān)系λ=-1.356×10-5Tinj2+0.005934Tinj+1.044，擬合精度SSE（殘差平方和）：1.827×10-31，R2：1，修正R2：NaN，RMSE（標(biāo)準(zhǔn)差）：NaN；在10 d時擬合溫度與分解前緣xf之間具有較好二項式關(guān)系xf=-4.286×10-5+0.01876Tinj+3.3，擬合精度SSE（殘差平方和）：7.889×10-31，R2：1，修正R2：NaN，RMSE（標(biāo)準(zhǔn)差）：NaN[圖7(a)]。擬合溫度與穿透深度xd具有較好的冪函數(shù)關(guān)系xd=4.515Tinj-0.09629-0.7516,擬合精度SSE（殘差平方和）：2.02×10-12，R2：1，修正R2：NaN，RMSE（標(biāo)準(zhǔn)差）：NaN；分解前緣8 m處擬合溫度與穿透時間td具有較好的二項式關(guān)系td=6.802×10-4T-2.534Tinj+35.94，擬合精度SSE（殘差平方和）：9.624×10-28，R2：1，修 正R2：NaN，RMSE（標(biāo) 準(zhǔn) 差）：NaN[圖7(b)]。

圖7 溫度敏感性參數(shù)關(guān)系Fig.7 Relationship curves of the sensitivity analysis of temperature

4 結(jié) 論

天然氣水合物是21世紀(jì)公認(rèn)的良好后續(xù)能源，其熱分解過程是一個考慮移動邊界的熱傳導(dǎo)Stefan相變問題。本文建立了單相無熱源連續(xù)水合物微元控制體的守恒積分熱傳導(dǎo)模型，在移動邊界分隔的水合物熱分解控制體各連續(xù)相（液相和固相）能量守恒方程基礎(chǔ)上，建立了尖銳界面上不連續(xù)熱焓和溫度梯度的耦合Stefan跳躍條件。以此為基礎(chǔ)，給出了一維半無限大純水合物藏?zé)岱纸鈹?shù)學(xué)模型及定解條件，利用Boltzmann相似變量求得了水合物熱分解Stefan相變模型的溫度分布和分解前緣的Neumann解，并對超越方程進(jìn)行了單調(diào)性證明，確定了Stefan模型解的唯一性。通過Deaton相平衡關(guān)聯(lián)方程和Clausius-Clapeyron方程計算了該純水合物藏的熱分解溫度和摩爾分解熱，利用MATLAB程序驗證了超越方程的單調(diào)性和Stefan模型解的唯一性，分析了溫度對超越方程的影響和解唯一性的影響，并開展了熱分解水合物體相變過程的溫度分布、分解前緣規(guī)律性研究和溫度對超越方程解和分解前緣的敏感性擬合研究。

符號說明

c——熱容，J/K

h——單位體積的熱焓，kJ/cm3

kl,ks——分別為液相和固相熱導(dǎo)率，W/(m·K)

L——水合物摩爾分解熱潛熱，W/mol

ng——甲烷氣體的量，mol

p——壓力，MPa

q——單位表面的傳熱量，kJ/m2

R——通用氣體常數(shù)，R=0.008314 MPa·m3/(kmol·K)

S,Sl,Ss——分別為控制體、液相分解區(qū)和固相水合物區(qū)的外表面積，m2

TR,Tinj,Tl,Tm,Ts——分別為初始、注入、液相區(qū)、分解和固相區(qū)的溫度，K

t,td——分別為分解時間和穿透時間，d

u——內(nèi)能，kJ

V,Vl,Vs——分別為控制體、液相分解區(qū)和固相水合物區(qū)的體積，m3

v n*——界面移動速度，m/s

xf,xd——分別為分解前緣位置和穿透深度，m

z——?dú)怏w壓縮因子

α——熱擴(kuò)散系數(shù)，m2/s

λ——超越方程變量

ρ——密度，kg/m3

下角標(biāo)

f——前緣

inj——注入

l——液相

m——分解

s——固相

引入Boltzmann相似變量ξ=xt[31]，則解的形式為：

F(ξ)是一個未知函數(shù)，將式(A2)代入熱傳導(dǎo)方程：

將上述熱傳導(dǎo)方程轉(zhuǎn)為二階線性齊次微分方程：

該方程具有積分因子解：

其中C1是積分常數(shù)。式（A3）兩端同時乘以M(ξ)，根據(jù)微分乘積法則：

積分可得：

C2是積分常數(shù)。由微積分基本定理可知，式（A6）的解為：

A為常數(shù)。所以方程的解為：