天然气水合物热分解Stefan相变数学模拟研究

李明川，樊栓狮，徐赋海，严科，黄爱先

（1中国石油大学（华东）石油工程学院，山东青岛266580；2华南理工大学化学与化工学院，广东广州510641；3中国石化胜利油田分公司东辛采油厂，山东东营257094；4中国石化胜利油田分公司胜利采油厂，山东东营257051）

引 言

天然气水合物是由碳氢气体和水形成具有笼型结构的、类似冰状的、可燃烧的晶体化合物，俗称“可燃冰”[1-3]。天然气水合物为高度压缩的天然气资源，1 m3天然气水合物能分解出160～180 m3（标准）的天然气[4]。天然气水合物广泛存在于海洋沉积物和永冻土中，据估计全球天然气水合物中的甲烷碳含量达1016kg或含有20×1015m3的甲烷气，因此天然气水合物被公认为是21世纪具有良好前景的后续能源[5]。

天然气水合物的热分解过程类似于冰的消融过程，是一个考虑移动边界的热传导相变问题，也称为两相Stefan相变问题[6]。Stefan边界问题的研究可以追溯到19世纪，Lamé等[7]在假设地球为半无限液体的条件下，研究了地球冷却过程中产生的固体地壳的厚度与时间的平方根呈正比的类似问题。Neumann假定初始温度高于熔点时，得出了半无限大Stefan问题的特解，简称Neumann解，他的成果直到20世纪初才由Weber出版发表出来[8]。最后解决这问题的是Stefan，他将问题转化为两个相的互变问题，利用Lamé等的假设条件，得到了更为完全的解。Stefan第一次从实验到理论对移动边界的相变问题进行了详细的研究，因此自由或移动边界问题也称为Stefan边界问题[9]。Stefan移动边界问题逐渐被广泛应用于固液相变、化学反应、溶质扩散、液气相变等诸多领域[10]。

天然气水合物热分解Stefan理论是一个新兴的跨学科应用领域，部分学者在这方面进行了大量研究工作。Makogon[11]首次使用类似经典Stefan的固体融化问题，描述了水合物的降压分解过程，获得了压力分布的自相似解。Selim等[12-13]认为水合物分解过程是一个移动边界消融问题，利用Von-karman积分得到了温度分布的近似解。Yousif等[14]建立了水合物分解界面的移动边界模型，并用该模型重新评估了水合物分解的气体体积。Tsypkin[15]建立了包含水合物分解前缘和冰融化前缘的两个未知移动边界的Stefan广义数学模型。唐良广等[16]将水合物分解过程看成界面运动过程，获得了热传导模型的解析解。李明川等[17]建立了水合物热分解过程Stefan径向模型，得到了Paterson指数积分函数精确解。Li等[18-19]提出了水合物热分解准稳态Stefan相变模型，确定了水合物热分解前缘位置数学模型。

综上，前人虽给出了水合物热分解Stefan相变模型，但未能详细给出建模、求解及检验过程，这有碍于学科、行业的应用性研究进展。作者对水合物热分解Stefan相变过程进行了系统性的研究，包括连续热传导模型、不连续界面Rankine-Hugoniot跳跃条件和Stefn耦合界面条件等理论。本文在单相连续水合物控制体能量守恒积分模型基础上，建立了考虑具有尖锐移动边界分隔的水合物热分解控制体的分解界面耦合Stefan能量守恒条件。给出了一半无限大天然气水合物热分解Stefan相变模型及定解条件，利用Boltzmann相似变量求得了误差函数形式的Neumann解，并对界面耦合Stefan超越方程进行了单调性证明，确定了Stefan模型解的唯一性，结合实例进行了分析计算。通过本文的理论研究工作，可促进水合物热分解相变规律、相变机理的揭示，为推动天然气水合物工业化开采提供理论依据。

1 控制体热传导Stefan模型

有一体积一定的连续的单相水合物控制体V，外表面积S，任取一微元dV和任一表面积dS，单位法向量n由控制体指向外部（图1）[20]。若控制体无内部热源和外来能量做功，能量随时间的变化率与通过表面积S传导到控制体V中的热量净速率相等（无对流），其数学形式为：

图1 连续水合物控制体Fig.1 Continuous hydrate control volume

式中，(ρε)为水合物单位体积的总能量，kJ/m3；q为通过单位表面积的热传导量，kJ/m2。总能量为内能、动能及势能之和，忽略势能，则内能为：

体积一定的水合物控制体V与时间t无关，热量(ρh)连续地通过微元体dV，式（4）左端的体积分和对时间的微分运算可互换，移动到积分式中，d/dt变为∂∂t；q=k·gradT连续地通过微元表面dS，式（4）右端可根据散度定理[21]，

式（7）即是单相连续区域的热传导方程。

接下来将式（4）应用于具有跳跃移动界面的水合物热分解控制体（图2）[22]。t时刻，界面Σ将控制体V分为液相分解区Vl和固相水合物区Vs，表面积分为Sl和Ss，单位法向量n l和ns分别指向各区域外方向。经过一小段时间δt后，界面移动到Σ'位置，变化的体积为δV，意味着在这一小段时间内有δV体积的固相水合物参与到热分解过程，相应的液相分解区增加δV的体积（忽略分解气体的影响），单位法向量n*由液相分解区Vl指向固相水合物区Vs。在这个过程中，由界面Σ到界面Σ'通过的焓h和温度梯度k·gradT是不连续的，但在各自的区域（液相分解区Vl和固相水合物区Vs）是连续的。为了推导界面的能量守恒条件，将式（4）左端分解到每个区域上分别积分，研究它们随时间的变化。

图2 界面不连续水合物控制体Fig.2 Interface discontinuous hydrate control volume

在时间t时刻:

注意到式（10）右端最后一极限项，表示在δt时间界面Σ→Σ'所增加的δV体积内的热焓值。当δt→0时，δV/δt=v n*dΣ，其中v n*是界面微元dΣ运动速度，方向垂直于界面Σ指向固相水合物；在δt→0时δV所在的空间趋近于表面Σ，水合物分解潜热L=hl-hs，则式（10）的δV积分变为Σ积分:

由式（4）分别代替式（11）中Vl、Vs上的积分，相应液相分解区和固相水合物在δt→0时增加了分解界面面积Σ，则:

这里单位法向量n*在界面Σ上由液相分解区指向固相水合物区，两相表面上的法向量有n l=n*，将式（12）右端按照两相表面Sl、Ss和界面Σ分别积分，且S=Sl+Ss，但积分式中Sl与Ss方向相反，则有：

将式（13）重新整理:

对任意的界面Σ，被积函数在界面上的所有点均消失，即环绕界面Σ积分等于0[23]。对比式（14）和式（4）可得：

式(15)即为水合物热分解的界面耦合Stefan能量守恒跳跃条件。

2 半无限大天然气水合物热分解Stefan模型及Neumann解

水合物分解界面条件：

定义水合物分解界面移动速度v n*=dxfdt，则液相分解区Vl和固相水合物区Vs的耦合Stefan条件由式（15）可得：

式（19）～式（24）构成了半无限大水合物藏一维Stefan模型，该模型的解是关于液相分解区和固相水合物区温度分布及未知前缘位置的移动边界的函数，而且温度与移动边界有关。由此可得方程的解（见附录）。

移动界面位置：

液相分解区温度：

固相水合物区温度：

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

对该超越方程进行单调性证明，以确定其解具有唯一性[24-25]。

3 结果与讨论

已知一半无限大一维纯甲烷水合物储藏，水合物物性如表1所示。

表1 水合物热分解一维模型基础数据［26-27］Table 1 One dimensional model basic data for thermal dissociation hydrate［26-27］

水合物在分解过程中保持恒定压力，忽略注入热水的影响，可通过Deaton相平衡关联方程[28]计算水合物的热分解温度Tm：

通过计算在压力5.55 MPa下，水合物的热分解温度为5.305089℃。理论计算表明，1个单位体积的可燃冰可以分解为164个单位体积的天然气及0.8个单位体积的水[29-30]，根据Clausius-Clapeyron方程[29-30]可计算水合物摩尔分解热L为14.83763 W/mol。

式中，ng为分解甲烷气体的量，mol；z为甲烷气体偏差因子；R为通用气体参数，R=0.008314 MPa·m3/(kmol·K)。

从图3可以看出，随着λ增加超越方程急剧下降，随后平缓变化，但总体上保持单调递减，验证了对式（28）的证明。在保持注入温度Tinj=100℃时，超越方程唯一解为λ=1.5021[图3(a)]；随着注入温度的增加（60，100，200℃），超越方程的解为λ=1.3515，1.5021，1.6888，对超越方程的影响逐渐增大，但增加幅度减小[图3(b)]。

图3 超越方程单调性及解关系曲线Fig.3 Relationship curves of the monotonicity and the solution of the transcendental equation

从图4看出，液相水合物分解区温度随着距离的增加急剧下降至前缘xf的分解温度Tm=5.305089℃；经历不连续的跃变后，温度穿透进入固相水合物区，下降幅度逐渐减小直至水合物藏温度TR=3.50℃，并保持不变。保持注入温度Tinj=100℃时，随着时间增加水合物分解前缘位置逐渐增大，温度穿透深度xd（定义：温度穿越分解前缘所能到达的距离，即xd=x TR-xTm）逐渐增大；5、10和20 d的水合物分解前缘xf分别到达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℃），液相水合物分解区温度下降幅度增大，水合物分解前缘位置增加；在10 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，这是因为固相水合物区温度下降梯度增大，在相同时间10 d的穿透深度降低了[图4(b)]。

图4 水合物体温度随距离关系曲线Fig.4 Relationship curves of temperature of hydrate control body vs.distance

从图5可以看出，随着时间增加，固相水合物区温度从水合物藏温度TR=3.50℃缓慢上升，至分解温度Tm=5.305089℃后温度迅速增大，但增加幅度逐渐降低。保持注入温度100℃时，达到分解前缘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℃），液相水合物分解区温度逐渐增大，随着时间增加上升幅度减小；不同注入温度到达指定前缘位置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所能达到的温度分别为43.7979、72.5967和144.6165℃，且增加幅度在扩大。

图5 水合物体温度随时间关系曲线Fig.5 Relationship curvesof temperature of hydrate control body vs.time

从图6可以看出，随着注入时间增加，水合物分解前缘逐渐增大，但增加幅度在减小。随着注入温度的增加（60、100、200℃），水合物分解前缘增加，10 d时分解前缘达到xf=4.2712、4.7412和5.3372 m，前缘增加幅度降低；随着注入温度的增加（60、100、200℃），到达分解前缘xf=8 m所需时间减小，分别为35.0814、28.3993和22.4674 d，减小幅度降低。

图6 水合物分解前缘随时间关系Fig.6 Relationship curves of the dissociation frontal brim of hydrate vs.time

接下来对温度敏感性参数进行分析。从图7可以看出，温度对λ和xf的影响逐渐增大，但增加幅度减小；温度对xd和td的影响逐渐减小，减小幅度降低。拟合温度与超越方程的解λ具有较好的二项式关系λ=-1.356×10-5Tinj2+0.005934Tinj+1.044，拟合精度SSE（残差平方和）：1.827×10-31，R2：1，修正R2：NaN，RMSE（标准差）：NaN；在10 d时拟合温度与分解前缘xf之间具有较好二项式关系xf=-4.286×10-5+0.01876Tinj+3.3，拟合精度SSE（残差平方和）：7.889×10-31，R2：1，修正R2：NaN，RMSE（标准差）：NaN[图7(a)]。拟合温度与穿透深度xd具有较好的幂函数关系xd=4.515Tinj-0.09629-0.7516,拟合精度SSE（残差平方和）：2.02×10-12，R2：1，修正R2：NaN，RMSE（标准差）：NaN；分解前缘8 m处拟合温度与穿透时间td具有较好的二项式关系td=6.802×10-4T-2.534Tinj+35.94，拟合精度SSE（残差平方和）：9.624×10-28，R2：1，修 正R2：NaN，RMSE（标 准 差）：NaN[图7(b)]。

图7 温度敏感性参数关系Fig.7 Relationship curves of the sensitivity analysis of temperature

4 结 论

天然气水合物是21世纪公认的良好后续能源，其热分解过程是一个考虑移动边界的热传导Stefan相变问题。本文建立了单相无热源连续水合物微元控制体的守恒积分热传导模型，在移动边界分隔的水合物热分解控制体各连续相（液相和固相）能量守恒方程基础上，建立了尖锐界面上不连续热焓和温度梯度的耦合Stefan跳跃条件。以此为基础，给出了一维半无限大纯水合物藏热分解数学模型及定解条件，利用Boltzmann相似变量求得了水合物热分解Stefan相变模型的温度分布和分解前缘的Neumann解，并对超越方程进行了单调性证明，确定了Stefan模型解的唯一性。通过Deaton相平衡关联方程和Clausius-Clapeyron方程计算了该纯水合物藏的热分解温度和摩尔分解热，利用MATLAB程序验证了超越方程的单调性和Stefan模型解的唯一性，分析了温度对超越方程的影响和解唯一性的影响，并开展了热分解水合物体相变过程的温度分布、分解前缘规律性研究和温度对超越方程解和分解前缘的敏感性拟合研究。

符号说明

c——热容，J/K

h——单位体积的热焓，kJ/cm3

kl,ks——分别为液相和固相热导率，W/(m·K)

L——水合物摩尔分解热潜热，W/mol

ng——甲烷气体的量，mol

p——压力，MPa

q——单位表面的传热量，kJ/m2

R——通用气体常数，R=0.008314 MPa·m3/(kmol·K)

S,Sl,Ss——分别为控制体、液相分解区和固相水合物区的外表面积，m2

TR,Tinj,Tl,Tm,Ts——分别为初始、注入、液相区、分解和固相区的温度，K

t,td——分别为分解时间和穿透时间，d

u——内能，kJ

V,Vl,Vs——分别为控制体、液相分解区和固相水合物区的体积，m3

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

xf,xd——分别为分解前缘位置和穿透深度，m

z——气体压缩因子

α——热扩散系数，m2/s

λ——超越方程变量

ρ——密度，kg/m3

下角标

f——前缘

inj——注入

l——液相

m——分解

s——固相

引入Boltzmann相似变量ξ=xt[31]，则解的形式为：

F(ξ)是一个未知函数，将式(A2)代入热传导方程：

将上述热传导方程转为二阶线性齐次微分方程：

该方程具有积分因子解：

其中C1是积分常数。式（A3）两端同时乘以M(ξ)，根据微分乘积法则：

积分可得：

C2是积分常数。由微积分基本定理可知，式（A6）的解为：

A为常数。所以方程的解为：