多孔材料中甲烷水合物生成的传热数值模拟研究

郎雪梅，姚柳眉，樊栓狮，李刚，王燕鸿，2

（1 华南理工大学化学与化工学院，广东 广州 510640； 2 广东省先进绝缘涂料工程技术研究中心，广东 珠海 519175；3 广东省燃料电池技术重点实验室，广东 广州 510640）

引 言

天然气作为21 世纪最重要的能源之一被广泛应用于交通运输和生产生活[1]。随着天然气需求的不断增加，安全、经济的天然气储运技术也更为重要。目前天然气的储运技术主要有管道运输、液化储运（LNG）、压缩储运（CNG）。管道储运广泛应用于工业和居民生活，适合长期持续的远距离输送，但工程投资高，灵活性差。天然气液化储运能够提供约600 m3CH4(STP)/m3的高储存容量[2]，目前主要用于远洋输送方面，但储存温度低以及连续不断的气体蒸发对生产设备和储存装置要求较高[3]。天然气压缩储运是将天然气加压到20～25 MPa 进行储运，适合于中小规模的储运，灵活性高。但高压储运增加了运输风险[4]。水合物储运（NGH）是近几年发展起来的天然气储运技术，因其固态储存、储气密度高、经济性高、安全性高得到越来越多的关注[5-6]。

水合物储运是利用水和天然气在低温高压条件下生成固体水合物进行天然气储运的技术。在标准条件下1 m3水合物可以储存180 m3的天然气[7-9]，且可以在冰点以下和大气压条件下储存数月[10-12]。根据日本三井公司的测算，在500 km 范围内NGH 运输成本比LNG 低18%～24%[13]，适合于短距离小规模的天然气运输。但天然气水合物的生成速率慢，限制了NGH 技术的应用[14]。Fan 等[15-16]利用泡沫铝、聚氨酯泡沫等多孔材料与动力学促进剂耦合，突破了水合物形成速度的障碍，实验室内实现了20 min 内将水合储甲烷密度提升至接近理论值。

从研究角度看，水合物储运已具备实现工业化的潜力。但水合物的生长是传质传热控制的反应，因此在放大实验中存在诸多不确定因素[17]，利用Fluent 数值模拟可以将水合反应中的复杂生成传热过程数值化，清晰了解反应器内部热量分布及反应程度，从而更好地指导水合储气反应进行放大。

Fluent 模拟可用于研究水合物分解[18]和水合物生成[19]。Neto 等[20]开发了一种预测海底管道水合物形成的Fluent 模型，该模型运用双膜理论探究气液相传质过程，基于传质动力学模型搭建水合物生成数学模型，但模拟结果未与实验验证。Song 等[21]采用了渗透率模型和反应动力学模型来研究砂质沉积物中水合物的非均匀生成，结果显示初始流体分布不同对水合物生成过程有很大影响，模拟与实验结果吻合良好。但关于水合物在反应器内的模拟研究未见报道。

本文针对水合物反应器中甲烷水合物生成传热过程，建立基于化学反应动力学模型和多孔材料内传质传热模型的甲烷水合物生成传热数学模型。通过与Fan 等[15-16]的实验结果比较验证数学模型准确性。该模型有助于分析填充多孔材料的反应器内水合物的生成情况，指导水合物放大反应器的设计。

1 多孔材料内甲烷水合物生成传热数学模型

1.1 生成反应

水合物的生成需要低温高压的环境[9]，是气体和液体生成固体的一个过程。假设本文温压条件下生成的甲烷水合物为sⅠ型水合物[22-23]，水合数n为5.75[9]。甲烷水合物的摩尔平衡如式(1)所示。

1.2 基本假设

多孔材料中甲烷水合物生成过程涉及温压变化、多相流以及化学反应。在数值模拟计算的过程中，需要应用水合物生成过程中的Arrhenius 方程、质量守恒方程和能量守恒方程。参照天然气水合物有关的理论和现有的研究成果[24-25]，建立了数学模型，并做出下面的假设：

(1)忽略体系中的动力学促进剂，假设生成反应中只由液态水和甲烷气生成水合物；

(2)模型建立只考虑三相，即水相、甲烷相和水合物相；

(3)忽略甲烷在水中的溶解效应和水蒸气的存在[26]；

(4)假设多孔材料内的所有相在同一时刻都具有相同的压力；

(5)忽略辐射换热和节流效应，流体与固体之间保持局部热平衡。

1.3 控制方程

本文选择基于化学反应动力学模型和多孔材料内传质传热模型[27]，建立甲烷水合物生成传热数学模型。模型中的水相、甲烷气相和水合物相质量守恒方程如下：

水相

式中，K是多孔材料的绝对渗透率，μm2；krw是水相的相对渗透率，μm2；krg是气相的相对渗透率，μm2；μw和μg分别为水和气相的黏度，mPa·s；φ是孔隙率；ρw、ρg和ρh分别为水、气和水合物相的密度，g/cm3；pg是气相压力，MPa；mw和mg分别为单位时间单位体积生成的甲烷水合物中消耗的水质量和气质量，g/(cm3·s)；mh是单位时间单位体积生成的甲烷水合物质量，g/(cm3·s)；Sw、Sg和Sh分别为水相、甲烷气相和水合物相的饱和度。

反应釜内部填充多孔材料，因此建立多孔介质模型。假设多孔材料是均匀分布，在模拟过程增加一个源项Si，该源项由两部分组成，一个黏性损失项α和一个惯性损失项Cn。黏性损失项为水相、气相和水合物相在各个方向上的渗透率，计算式如下[28]。

式中，Si为动量方程中i方向上的源项，kg/(m2·s2)；α为相对渗透率，m-2；ui为i方向流体速度，m/s；Cn为惯性阻力系数，m-1；u为流体速度，m/s；Dp为孔径，m。

不同反应时间水相、甲烷气相和水合物相饱和度的具体数值可从Fluent 中得到，它们之间关系如下。

式(4)中的质量源项mh采用Arrhenius 方程进行建模[29-30]，并结合水合物反应界面面积进行校正。模型搭建过程参考了Kim 等[31]的水合物分解动力学模型。当p≥pe时，可作为水合物生成数学模型。因此甲烷水合物质量源项mh可用式(9)计算。

式(13)有以下假设：水合物生长基于多孔材料孔隙表面；反应界面面积取决于基底面和气相水相之间的接触可能性。式中，Ap是多孔材料内的孔隙表面积以及水合反应过程中形成的水合物表面积，m2/m3，使用Kozeny-Carman 方程进行估算[34]，用式(14)计算。

式中，φe为有效孔隙率，定义为气相和水相体积与全部孔隙体积的比值；Г是水合反应中活跃的孔隙表面积分数。绝对渗透率具有如下形式[35]，随水合物饱和度的增加而减小。

多孔材料区域中气体、水和水合物之间保持总能量平衡，流体与固体之间保持局部热平衡。总能量是内能和动能之和，能量平衡也就是控制体内各相内能、相间流动能与反应热的动态平衡，因此能量守恒方程经过简化，为

式中，C1和C2均为常数，C1=5.31×107J/kmol；C2=-15850 J/(kmol·K)。

针对上述数学模型构建三相质量源项和水合物能量源项的UDF 程序。UDF 程序利用Fluent 软件中的预定义Define 宏命令实现。使用“#DEFINE”宏命令定义源项中全局变量，DEFINE_SOURCE 宏命令定义程序中三相质量源项和能量源项的变化，因此数学模型较为复杂，于是采用编译UDF 的方式实现程序在Fluent中的装载。

2 数值模拟条件及方法

2.1 物理模型与网格划分

通过与实验数据对比验证建立的甲烷水合物生成传热数学模型准确性。根据文献[15-16]报道的多孔材料水合物储气实验进行建模。图1为实验中反应釜的模型，上方为气体填充区域，下方为多孔材料区域；多孔材料区域的高度根据实验过程填充多孔材料量和溶液量来确定，用该物理模型进行甲烷水合物生成传热模拟实验。

图1 物理模型Fig.1 Physical model

网格划分工具采用Fluent Meshing，验证模型的网格划分采用先划分表面网格再划分体网格的方法，同时划分边界层和加密流体域连接之间的网格。最小网格尺寸为2 mm，体网格划分采用了Poly-Hexcore网格，图2为模型的网格划分。

图2 模型网格划分Fig.2 Grid system of model

2.2 边界条件及初始化

选择聚氨酯泡沫（PU）和泡沫铝（AF）作为多孔填料进行水合物储甲烷实验模拟验证[15-16]。水合物设置选用冰的基本参数，再根据甲烷水合物的物理参数进行修改；在单元区域条件中根据多孔材料参数设置多孔区域条件，将多孔材料区域和气体填充区域的交界面设置为interior 边界条件。反应釜壁面无滑移且设置为温度壁面条件。采用标准初始化，两组模拟的运行参数见表1。

表1 模拟运行参数Table 1 Operating parameters of the simulation

多孔材料区域水合物生成时气、水、水合物饱和度变化对釜内介质表观热导率λa的影响如式(20)。

式中，λw、λg、λh和λp分别为水、甲烷、甲烷水合物和多孔材料的热导率，W/(m·K)。

2.3 数值模拟方法

本文采用商业计算流体力学软件Ansys Fluent进行数值模拟，使用瞬态计算，时间步长根据实际需求调整。多相流模型选择欧拉多相流，湍流模型选择Realizableκ-ε湍流模型（反应釜内的计算单元存在因水合物生成而产生压降的局部湍流，故使用湍流模型）。采用双精度求解器、Phase Coupled SIMPLE 算法，除压力项采用二阶精度来提高计算准确性，湍流、动量和能量等方程中对流项均采用一阶精度。松弛因子设定保持默认，以标准化残差合理振荡波动判断收敛。

3 模拟结果及分析

3.1 模型验证

图3 为文献[15-16]中甲烷水合物在PU 和AF 多孔材料中的生成动力学结果。用本文描述的数值模型方法进行模拟，模拟的时间根据文献[15-16]中水合物快速生成的时间段进行选择，分别模拟了33.5 min 和13.5 min。模拟初始填充气体量、水量和多孔材料量按照实验设定，模拟过程壁面温度、反应釜内压力和多孔材料孔隙度的设定均与两组实验过程保持一致。

图3 甲烷水合物生成动力学[15-16]Fig.3 Kinetic of methane hydrate formation[15-16]

图4(a)、(b)为甲烷水合物在聚氨酯泡沫PU 和泡沫铝AF 中的生成动力学模拟数据与实验数据的比较，绝对平均相对误差AARD 分别为4.38% 和5.22%，可以看到实验和模拟值生成趋势接近，模拟储气量和实验储气量基本一致，表明该数值模型建模中假设简化合理，所建甲烷水合物生成传热模型适用于该体系。这两组数据也可以明显看到水合物在泡沫铝多孔材料中比在聚氨酯泡沫多孔材料中生成更快。

图4 模拟数据与实验数据比较Fig.4 Comparison of simulation data and experimental data

3.2 多孔材料热导率对水合物生成传热过程的影响

图5 为聚氨酯泡沫内z=0 平面甲烷水合物体积分数随时间的变化。可以看到水合物生成区域在多孔介质中。从前15.0 min 云图中看到反应初期水合物生成比较均匀，在整个多孔介质填充区域均有生成。15.0 min 以后，随着生成量的增加，水合物呈现中部区域生成量低于靠近冷壁面区域，气液相交界面生成量较多。图6 为聚氨酯泡沫内z=0 平面甲烷水合物温度随时间的变化。可以看到反应釜内的温度分布与水合物生成量密切相关。10.0 min 之前，系统整体温度较低，水合物生成较快，释放出大量的热。由于聚氨酯泡沫导热差，随着水合物生成量增大，多孔材料区域温度快速升高，尤其是中心区域，超过276 K。到15.0 min时，除冷壁面附近，大部分区域温度升至276 K 以上。受到温度升高的热抑制和釜内压力降低影响，水合物生成速率开始逐渐减缓，冷量从壁面缓慢传递至反应釜中心，釜内温度开始下降。到33.0 min 左右，水合物生成基本完成，温度也恢复到环境温度。

图5 聚氨酯泡沫PU内z=0平面甲烷水合物体积分数随时间的变化Fig.5 Methane hydrate volume fraction at z=0 plane with time variation in PU foam system

图6 聚氨酯泡沫PU内z=0平面甲烷水合物温度随时间的变化Fig.6 Methane hydrate temperature at z=0 plane with time variation in PU foam system

图7 为泡沫铝内z=0 平面甲烷水合物体积分数随时间的变化，可见在泡沫铝内的水合物生成传热过程与聚氨酯泡沫中明显不同。前10.0 min，泡沫铝材料区域水合物的体积分数保持均匀分布。这是因为初期釜内温度较低，泡沫铝材料又能提供较多的成核位点[15]，水合物可以在整个区域内均匀生成。由于泡沫铝热导率为聚氨酯泡沫的6000倍，尽管水合物大量生成放热，但并未出现PU 体系中出现的中心过热现象，而是在整个AF 填充区域无较高温升，见图8 水合过程温度云图。AF 及时将中心区域的热传递到壁面而被环境吸收，最高中心温度为274 K。整个储气过程持续13.5 min 完成，较PU体系缩短20.0 min。

图7 泡沫铝AF内z=0平面甲烷水合物体积分数随时间的变化Fig.7 Methane hydrate volume fraction at z=0 plane with time variation in AF foam system

图8 泡沫铝AF内z=0平面甲烷水合物温度随时间的变化Fig.8 Methane hydrate temperature at z=0 plane with time variation in AF foam system

图9 为水合反应中心的水合物、气、水饱和度和釜内介质表观热导率随时间的变化，可从看出，随着水合物饱和度增加和气、水饱和度下降，两个多孔材料体系中介质表观热导率都随之增加。反应结束时PU体系和AF体系介质表观热导率分别为0.23 W/(m·K)和64.82 W/(m·K)。

图9 水合反应中心的水合物、气、水饱和度和介质表观热导率随时间的变化Fig.9 Plot of Sh,Sg,Sw and λa of hydration reaction center with time

为了能够更加直观比较两种多孔材料的导热能力，对不同高度的水合物温度数据进行汇总，图10 为PU 和AF 体系内水合物温度计算值取点位置，在反应釜的中线位置，z=0 mm 和x=0 mm 处，沿y轴均匀取点。均在多孔材料区域取两个测温点，一个测温点在多孔材料中心区域，一个测温点在多孔材料区域和气体填充区域交界面，并按照两点间距取气体填充区域测温点，因此PU 体系多孔材料区域每隔20 mm 取一个点，多孔材料区域2 个位置点，气体填充区域4 个位置点；AF 体系每隔18 mm取一个点，多孔材料区域2个位置点，气体填充区域5个位置点。图11(a)、(b)分别为聚氨酯泡沫PU 体系和泡沫铝AF 体系内不同高度甲烷水合物温度随时间的变化，可以明显看出两个体系最高温度相差9.0 K。甲烷水合物生成时会放出大量的热，造成周围温度上升。由于多孔材料填充区域是水合物生成区，温度波动比较大，气体填充区域温度波动较为平缓。在PU 体系中，由于甲烷水合物初期大量生成，造成温度上升，而PU 的热导率较小，生成热无法及时排出，从而出现了局部高温。过高的温度不仅减缓了水合物的生成，还造成周围水合物的分解，分解吸收热量又引起了温度的下降。故而出现了温度骤升骤降的现象。而在AF 体系中，AF 的热导率高，热量传递快，水合物生成热及时通过泡沫铝骨架导向壁面，故图中的温度波动更加平缓，温差更小。PU 体系中，从5.0 min 到30.0 min，水合物生成中心区域一直保持较高温度，最高温度甚至升高至283 K，直到水合反应基本完成温度才开始下降，而气体填充区域温度则基本保持在273.0～274.0 K 内；而AF 体系中，AF 填充的水合区域温度基本保持在273.0～274.0 K 内，该区域温度在快速反应阶段略有升高但也很快降至环境温度，反应釜整体最高温升不超过1.5 K。

图10 聚氨酯泡沫PU和泡沫铝AF体系内水合物温度计算值取点位置Fig.10 The point of hydrate temperature in PU and AF porous materials system

图11 不同高度甲烷水合物温度随时间的变化Fig.11 Methane hydrate temperature variation with time at different heights

综上，水合反应中，热量传递是影响水合物生成速率的关键因素之一。填料热导率越大，越有利于水合生成热移除以及水合物生成。导热不良时，容易在中心形成局部过热，对水合物生成造成热抑制。

4 结 论

本文建立了多孔材料内甲烷水合物生成传热数学模型，通过与实验数据对比验证了该模型的可靠性，并对不同热导率的多孔材料内水合反应过程进行数值模拟，得到以下结论。

（1）模型验证中AARD 小于6%，实验和模拟储气量趋势接近，表明该水合物生成传热模型基本准确，可靠性高。

（2）水合反应过程中，热量传递是影响水合物生成速率的关键因素之一。

（3）填料的热导率越大，越有利于水合生成热的移除以及水合物的生成。导热不良时，容易在中心部分形成局部过热，对水合物生成造成热抑制。

（4）该模型可用于水合反应器的设计和优化。在进行甲烷水合反应器设计时，应重点考虑反应器内的传热控制。