天然氣水合物熱分解Stefan相變數(shù)學(xué)模擬研究

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

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

引 言

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

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

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

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

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

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

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

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

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

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

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

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

在時(shí)間t時(shí)刻:

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

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

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

將式（13）重新整理:

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

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

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

水合物分解界面條件：

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

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

移動(dòng)界面位置：

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

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

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

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

3 結(jié)果與討論

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

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

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

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

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

從圖3可以看出，隨著λ增加超越方程急劇下降，隨后平緩變化，但總體上保持單調(diào)遞減，驗(yàn)證了對(duì)式（28）的證明。在保持注入溫度Tinj=100℃時(shí)，超越方程唯一解為λ=1.5021[圖3(a)]；隨著注入溫度的增加（60，100，200℃），超越方程的解為λ=1.3515，1.5021，1.6888，對(duì)超越方程的影響逐漸增大，但增加幅度減小[圖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℃時(shí)，隨著時(shí)間增加水合物分解前緣位置逐漸增大，溫度穿透深度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時(shí)間，分解前緣分別到達(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，這是因?yàn)楣滔嗨衔飬^(qū)溫度下降梯度增大，在相同時(shí)間10 d的穿透深度降低了[圖4(b)]。

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

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

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

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

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

接下來(lái)對(duì)溫度敏感性參數(shù)進(jìn)行分析。從圖7可以看出，溫度對(duì)λ和xf的影響逐漸增大，但增加幅度減小；溫度對(duì)xd和td的影響逐漸減小，減小幅度降低。擬合溫度與超越方程的解λ具有較好的二項(xiàng)式關(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時(shí)擬合溫度與分解前緣xf之間具有較好二項(xiàng)式關(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處擬合溫度與穿透時(shí)間td具有較好的二項(xiàng)式關(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ù)能源，其熱分解過(guò)程是一個(gè)考慮移動(dòng)邊界的熱傳導(dǎo)Stefan相變問(wèn)題。本文建立了單相無(wú)熱源連續(xù)水合物微元控制體的守恒積分熱傳導(dǎo)模型，在移動(dòng)邊界分隔的水合物熱分解控制體各連續(xù)相（液相和固相）能量守恒方程基礎(chǔ)上，建立了尖銳界面上不連續(xù)熱焓和溫度梯度的耦合Stefan跳躍條件。以此為基礎(chǔ)，給出了一維半無(wú)限大純水合物藏?zé)岱纸鈹?shù)學(xué)模型及定解條件，利用Boltzmann相似變量求得了水合物熱分解Stefan相變模型的溫度分布和分解前緣的Neumann解，并對(duì)超越方程進(jìn)行了單調(diào)性證明，確定了Stefan模型解的唯一性。通過(guò)Deaton相平衡關(guān)聯(lián)方程和Clausius-Clapeyron方程計(jì)算了該純水合物藏的熱分解溫度和摩爾分解熱，利用MATLAB程序驗(yàn)證了超越方程的單調(diào)性和Stefan模型解的唯一性，分析了溫度對(duì)超越方程的影響和解唯一性的影響，并開(kāi)展了熱分解水合物體相變過(guò)程的溫度分布、分解前緣規(guī)律性研究和溫度對(duì)超越方程解和分解前緣的敏感性擬合研究。

符號(hào)說(shuō)明

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——分別為分解時(shí)間和穿透時(shí)間，d

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

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

v n*——界面移動(dòng)速度，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(ξ)是一個(gè)未知函數(shù)，將式(A2)代入熱傳導(dǎo)方程：

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

該方程具有積分因子解：

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

積分可得：

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

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