Multiphase flow computational model for extraction of gas hydrates in marine soft soils
-
摘要: 天然气水合物是一种储量巨大的非常规清洁能源,其安全、高效开采对于保障我国能源安全意义重大。降压法是目前水合物开采中使用最广、最具前景的一种开采方式。本文针对水合物降压开采中动态多相渗流过程开展理论及数值模拟研究。考虑水合物分解与多相渗流之间的动态耦合,建立了耦合水合物动力学分解与多孔介质非饱和渗流的数学模型,并基于COMSOL多物理场耦合有限元软件进行数值实现。模型通过与Masuda水合物分解试验进行对比验证了有效性。基于数学模型对水合物开采的关键影响因素开展敏感性分析,分析结果表明:水合物饱和度降低速率随储层渗透率、降压幅度的增加而加快。Abstract: The gas hydrate is a type of unconventional clean energy with substantial reserves. The successful and safely extraction of gas hydrates is of great significance in guaranteeing the energy safety of China. At present the depressurization method is a kind of the most widely used and promising extraction means for the gas hydrates. The problems of multiphase flow in the reservoir caused by hydrate phase transition are studied theoretically and numerically. Firstly, a mathematical model for coupling the kinetic decomposition of the gas hydrates and the multiphase flow in porous media is established. The proposed model is numerically implemented by the software COMSOL Multiphysics, and its effectiveness is validated through the Masuda's experiments. Then, the sensitivity analysis for the key factors which have impact on the extraction of the gas hydrates is conducted. The results show that the gas production rate increases with the increase of the permeability and pressure drop amplitude.
-
Keywords:
- gas hydrate /
- depressurization method /
- multiphase flow /
- phase transition /
- energy geotechnics
-
0. 引言
土地荒漠化是当今世界面临的一个严重的环境和社会经济问题。荒漠化土地面积的迅速扩展,造成大范围土地退化和巨大经济损失[1]。中国土地荒漠化和沙化情况严重,根据国家林业局发布的全国荒漠化和沙化状况公报[2],近年来尽管荒漠化、沙化土地面积逐年减少,但总体情况仍较为严重:根据第五次全国荒漠化和沙化监测工作结果,中国荒漠化土地面积达261.16万km2,占国土面积27.20%;沙化土地面积达172.12万km2,占国土面积17.93%。土地荒漠化的加剧,导致沙尘暴等恶劣天气频发,对经济发展、生态文明、居民健康都造成了严重影响。因此,沙漠治理迫在眉睫。
近年来,生物诱导碳酸钙沉积技术作为一种新型土体固化和防渗技术,被广泛应用于土木工程和岩土工程领域[3]。生物诱导碳酸钙沉积技术目前主要分为两类:微生物诱导碳酸钙沉积[4](microbial induced calcite precipitation, MICP)和酶诱导碳酸钙沉积[5](enzyme induced calcite precipitation, EICP)。两者作用机理类似,利用脲酶的尿素水解作用诱导产生碳酸钙晶体,在岩土体材料中起到填充和胶结作用,从而实现土体加固的目的[6-8]。二者区别在于,MICP方法利用的是含脲酶的活性微生物,EICP所需脲酶事先从含脲酶的微生物和植物细胞中提取获得。对于MICP方法,其脲酶来源主要为细菌(如巴氏生孢八叠球菌)[9],加固粗粒土效果显著,能有效提高处理土体强度和防渗性能,但细菌培养过程复杂,且缺乏有效的监测难以评定其对环境的影响。对于EICP方法,其固化细粒土效果较好,且脲酶来源可以为植物,更加容易获得。高纯度的商品化脲酶成本高,难以大规模使用。豆类植物脲酶含量最为丰富[10],从豆类植物中快速提取粗酶液,可有效降低成本[11]。
考虑将大豆脲酶诱导碳酸钙技术与黄原胶相结合,对前者进行优化。黄原胶是一种生物胞外多糖,是水凝胶的一种,通常在食品和医学中,用作增稠剂、稳定剂、黏合剂等[12]。近年来,有研究将水凝胶作为一种新型环境友好型土体改良材料,进行土体加固[13]、土体防渗[14]和防扬尘[15-16]等研究,并取得了显著的研究成果。黄原胶对土体进行处理,可以起到保水和固化等作用,且黄原胶的固化作用具有延性,这和SICP固化中的碳酸钙脆性固化不同,二者共同可能会起到更好的效果。
本文将大豆脲酶诱导碳酸钙沉积技术与黄原胶相结合,通过测定不同浓度混合溶液处理试样的表面强度、风沙吹蚀强度、钙离子含量,并结合微观观测,评价大豆脲酶诱导碳酸钙沉积与黄原胶结合的防风固沙效果,并探究其机理。
1. 试验材料与试验方法
1.1 试验材料
(1)试验用砂
试验用砂取自宁夏中卫沙坡头区域,采用筛分法测得粒径集中分布在区间0.075~0.25 mm,属于细砂,颗粒形状近于球形,级配不良,粒径分布曲线见图1。砂土土粒比重为2.66,自然状态下含水率为0.31%。干燥状态下测得砂样最小干密度为1.44 g/cm3,最大干密度为1.72 g/cm3。对砂样进行X射线衍射分析(XRD),其主要成分为二氧化硅。
(2)大豆脲酶胶结液
试验所用大豆脲酶提取步骤如下[11]:
a)将大豆烘干置于粉碎机中磨成豆粉;
b)按照固液比为1∶25取豆粉溶于去离子水中,充分搅拌后静置1 h,待豆渣沉淀;
c)将豆粉溶液放入离心机离心15 min,离心机参数设置为4℃、3000 r/min;
d)取离心管中上清液,即得大豆脲酶溶液。
用电导率法测量脲酶活性[17](电导率仪为雷磁DDB-303A型),经测定,本文试验所用大豆脲酶水解尿素的活性在25°C室温环境下约为6.7 mmol/(L·min)。
将大豆脲酶溶液和盐溶液(氯化钙-尿素溶液)1∶1混合,即为大豆脲酶胶结液。
(3)黄原胶
试验用黄原胶为食品级黄原胶,该品类黄原胶价格经济,且已有试验证明黄原胶可用于土体固化[13]。将黄原胶溶于去离子水配成黄原胶溶液,与大豆脲酶胶结液1∶1混合,即为处理液。
1.2 试验方法
(1)试样制备与处理
试验土样均用金属浅盘制做,金属盘尺寸为:27 cm×20 cm×4.8 cm(长×宽×高)。将取自宁夏沙坡头的风积砂过0.5 mm筛,筛除里面的树枝、枯草等,控制各试样密度大约为1.5 g/cm3。按照4 L/m2的用量,使用加压式喷壶将处理液均匀喷洒在砂盘表面。为了让碳酸钙充分沉积,处理完毕的试样置于室温下反应24 h,待反应充分完成后,将试样置于60℃环境烘干,进行后续试验。采用这种试验方法是为了模拟防风固沙的实际环境条件。在表层土湿润的状态下,一般不存在土体吹蚀问题,因而需要在接近环境条件(地表温度)下进行干燥,测试其防风性能。
考虑处理液浓度对试验的影响,通过控制变量设置不同处理液成分浓度进行对比,每种配方制作3个平行试样,处理液配方见表1。
表 1 不同处理液配方Table 1. Composition of different treatment liquids配方编号 大豆脲酶胶结液浓度(混合后)/(mol·L-1) 黄原胶溶液浓度(混合后)/(g·L-1) U-0-0 纯去离子水处理空白样 T-0-0.1 — 0.1 T-0-0.3 — 0.3 T-0-1 — 1 T-0.1-0 0.1 — T-0.3-0 0.3 — T-0.1-0.1 0.1 0.1 T-0.1-0.3 0.1 0.3 T-0.1-1 0.1 1 T-0.3-0.1 0.3 0.1 T-0.3-0.3 0.3 0.3 T-0.3-1 0.3 1 注: 试样中编号U表示未处理试样(untreated soil),T表示处理试样(treated soil),后面的数字表示大豆脲酶胶结液浓度(0.1,0.3 mol/L),最后的数字表示黄原胶溶液浓度(0.1,0.3,1 g/L)。(2)表面强度试验
使用数显式推拉力计对烘干试样进行表面强度测量,测量探头为平头,截面为圆形,面积约19.6 mm2。测量时,将推拉力计置于水平台面上,转动手柄使探头以3 mm/s恒定速度匀速贯入试样表面,读取峰值贯入力。考虑到喷洒处理无法完全均匀、土体本身的性质存在差异以及为减小试样处理和反应过程中的随机误差,每个试样取10个测量点,取平均峰值贯入力记为该试样表面强度表征值。
(3)碳酸钙含量测定
由于碳酸钙晶体和黄原胶的填充、胶结作用,在处理试样的表面会形成一层硬壳层,对硬壳层取样进行碳酸钙含量测定,以评价黄原胶对处理试样碳酸钙生成的影响。按国家标准《水质钙的测定EDTA滴定法》[18]进行,每个试样进行3次平行测定,取平均值换算后得到试样表面生成的碳酸钙含量,并计算加入钙盐转换为碳酸钙的转化率。
(4)风沙吹蚀试验
评定处理试样防风固沙性能的直观试验即风蚀试验。考虑沙漠实际风蚀情况,设置不含沙风蚀和含沙风蚀两种情况,使用两组平行试样对照,试验前测量每盆试样的烘干质量。根据风速与风力的表达关系式V=0.84×F2/3,F为蒲氏风级,以沙坡头地区为例统计宁夏风力情况,如表2所示,在近5年统计时间内,大部分时间沙坡头风力为3级以下,对应风速为5 m/s以下,仅有1 d,风力达7级。以此为依据设计在5,10,15 m/s(分别对应风力等级为3级、5级、7级)风速条件下吹蚀试样,同时,根据水土保持学说[19],随着风速的增加,输沙量增加,以此为参考,考虑实际试验条件,以每分钟吹过土样单位表面积的沙量为指标,设定含沙风蚀试验组中,5,10,15 m/s风速条件下的输沙量分别为18.5,37,55.5 kg/ (m2·min)。风蚀试验持风时间为15 min,测试样的质量损失率,并观察表面硬壳层破损情况,进行比较。为评价恶劣条件下处理试样的防风固沙性能,将15 m/s风蚀时间持续至60 min,计算试样的质量损失率,并观察表面硬壳层破损情况。风沙吹蚀试验装置如图2所示。
表 2 宁夏沙坡头风力统计(2014-07—2019-11)Table 2. Wind power statistics in Shapotou, Ningxia风力 无风 1级(0.3~1.5 m/s) 2级(1.6~3.3 m/s) 3级(3.4~5.4 m/s) 天数/d 530 292 524 275 风力 4级(5.5~7.9 m/s) 5级(8.0~10.7 m/s) 6级(10.8~13.8 m/s) 7级(13.9~17.1 m/s) 天数/d 68 41 5 1 (5)SEM扫描电镜试验
取处理土样的硬壳层,用去离子水轻轻冲洗去硬壳层表面未反应完的无机盐与杂质后,于60°C环境烘干,随后用镊子小心夹取表面中心的薄片黏于金属板上进行喷金处理(如图3所示),喷金处理完成后,利用扫描电子显微镜(飞泉200环境扫描电镜),从微观角度直观观测单纯黄原胶作用、单纯大豆脲酶胶结液作用以及大豆脲酶胶结液结合黄原胶共同作用下,土颗粒联接形式的不同,揭示优化机理。
2. 结果与分析
2.1 碳酸钙生成量
如图4所示,SICP和黄原胶共同处理的土体,碳酸钙含量随着SICP浓度升高而升高,而钙离子转化率均在50%~60%波动,几乎不受SICP及黄原胶浓度影响,可以认为在本文试验参数范围内,黄原胶不会抑制大豆脲酶诱导碳酸钙生成过程。与此同时,黄原胶的抗酶解反应特性[12]也使得其在大豆脲酶诱导碳酸钙沉积过程中能保持自身性质。
2.2 表面强度
如图5所示,与纯去离子水处理试样U-0-0相比,黄原胶溶液和大豆脲酶诱导碳酸钙沉积均能提高处理土体的表面强度,且随着浓度增加,处理土体表面强度增加。总体上,碳酸钙对土体的胶结作用强于黄原胶对土体的胶结作用,当单纯黄原胶用量达到1 g/L时,试样峰值贯入力为23.9 N,而单纯SICP处理浓度达0.3 mol/L时,峰值贯入力为73 N。将大豆脲酶胶结液与黄原胶结合后,处理土体的表面强度提高明显,其联合胶结强度优于二者胶结强度简单叠加的效果,峰值贯入力可达138 N。
黄原胶通过共价键作用以及黏结土颗粒固化土体,固化作用具有延性,大豆脲酶胶结液通过诱导生成碳酸钙晶体胶结土颗粒,固化作用具有脆性,两者结合起到了更好的固化效果,表明黄原胶对大豆脲酶诱导碳酸钙沉积处理土体有改良优化作用。
2.3 风沙吹蚀试验
如图6所示,可以看到只用去离子水喷洒的空白样U-0-0在10 m/s 的风速下被吹蚀,只用0.1 g/L黄原胶溶液处理的土样,在10 m/s风速下开始出现质量损失,并随着风速增大、吹风时间增长,质量损失增大。其余试样在不含沙风蚀试验中均未出现破坏情况,质量损失微小或完全无损失。
含沙风蚀试验结果如图7~9所示,试样风蚀质量损失率均随着风速增加、吹蚀时间增长而增长。并且,试样的抗风能力随着SCIP处理液和黄原胶的增加呈现增强趋势。其中,只喷洒黄原胶溶液的土样,同等条件下的质量损失率也较大(图7),纯大豆脲酶胶结液处理的试样次之,将大豆脲酶胶结液和黄原胶共同作用后,质量损失率显著降低(图8,9)。
单纯黄原胶处理的试样,实际固化物含量较低,土颗粒之间仅存在表层固结,无法承受高风速、长时间吹蚀,胶结作用弱于碳酸钙。而当黄原胶与碳酸钙晶体共同作用,土颗粒之间黏结的更加牢固,处理土样的抗拉强度和胶结程度提高,对抗风剪和风夹带沙颗粒冲击磨蚀作用增强,整体抗风沙吹蚀性能得到了提升。
对比不含沙风蚀与含沙风蚀的试验结果,在15 m/s不含沙风吹蚀1 h工况下,仅处理所得强度最低的T-0-0.1试样质量损失率为3.8%,其余试样均无质量损失,而在15 m/s含沙风吹蚀1 h工况下,所有试样均有不同程度质量损失,损失率在0.7%~66.7%不等,其中,处理所得强度最低的T-0-0.1试样质量损失率达66.7%。
沙漠现场的风往往夹杂着沙子形成风沙流吹蚀,风沙流不仅包含风的剪切作用,更有沙颗粒对土体表面的冲击和磨蚀作用,对土体的破坏要强于不含沙风蚀,且随着风速提高,风夹带沙粒量增多,对地剪切力和沙粒冲击磨蚀作用都显著增大。本文在生物诱导碳酸钙沉积相关研究领域首次考虑风夹带沙颗粒侵蚀作用。试验中,当风速大于5 m/s(临界起沙速度)时,沙颗粒克服重力作用,在风力夹带下形成风沙流,由于本文所用沙漠现场沙粒径集中分布在75~250 μm,风夹带沙粒运动方式以跃移为主,跃移颗粒在运动过程中不断从气流中吸收能量,冲击地面,这种冲击效能往往引起土体颗粒的分离[20]。同时,由于处理土样表面固结形成屏障,含沙粒的气流绕流障碍物时,到达障碍物正前缘的气流受到阻挡,往障碍物的两侧流动。结果使通过两侧的气流增多,流速增加。这些沙粒将对障碍物前缘和侧边的沙面产生冲击效果,后续气流将把飞起沙粒带走,一旦形成小的侵蚀沟坑,气流会在这些位置脱离和回流,加重对沟槽的冲击,使沟槽加深 [21]。因此,在测定处理试样防风固沙性能时,用含有沙子的风进行吹蚀,更贴合沙漠实际情况。
统观试样风蚀量与表面强度的关系,如图10所示,基本呈现出,随着表面强度增加,试样风蚀量降低、抗风沙侵蚀能力增加的趋势,当峰值贯入力达到70 N时,风蚀量在所有工况下均低于5%,表明试样抗风蚀强度与表面强度存在正相关关系。随着胶结方案的优化,土体整体的强度提高,薄弱区域也得到了加固,从而能达到较好的防风固沙效果。
观察土样表面硬壳层在风沙吹蚀过程中的破坏发展情况(以图11,12为例),同等吹蚀条件下,加入黄原胶的SICP处理的土样表面硬壳层保留程度较好,磨蚀处亦未贯穿硬壳层,且当黄原胶浓度提高至1 g/L后,在15 m/s,输沙量为55.5 kg/(m2﹒min)的长时间吹蚀情况下,硬壳层仍保留完好,虽有磨蚀,但未出现破裂,最大磨蚀深度仅3 mm,远小于硬壳层厚度2 cm。表明黄原胶的加入,优化了大豆脲酶诱导碳酸钙的单一胶结模式,提高了土样表面硬壳层整体性,延缓了土样破坏的进程,改善了处理土样防风固沙性能。
考虑到宁夏沙坡头地区风力情况(如表2),认为15 m/s吹蚀1 h工况可以表征实际至少一年的风蚀破坏情况,即大豆脲酶胶结液与黄原胶联合防风固沙技术理论固沙时长可达至少一年,这有利于将该技术与植物固沙相结合,为植物生长提供了充足时间条件。因此,用大豆脲酶胶结液与黄原胶结合,对前者进行改进优化来处理沙漠,防风固沙,具有实际工程应用前景。
2.4 SEM扫描电镜试验
图13展示了黄原胶处理土样表面硬壳层的胶结形态,呈现出胶体固化的形态。图14展示了大豆脲酶诱导碳酸钙沉积处理土样表面硬壳层的胶结形态,呈现出结晶的形态。在图15中,在土颗粒与土颗粒接触的位置,同时存在胶体黏结与矿物晶体胶结作用,这可能是大豆脲酶诱导碳酸钙沉积结合黄原胶处理的土样,表面强度和抗风沙吹蚀强度比单一处理方式得到的强度高的原因:由于胶体具有黏性和流动性,在土体中,可以在土颗粒之间流动、附着、包裹,并且在重力作用下,更易于流动聚集至低势能点(即颗粒间空隙处)。
在这样的流动过程中,沙子颗粒以及晶体之间的接触面积增加,在胶体的包裹、黏结、填充下,处理土体的整体性提高。且胶体残余水分的脱水会导致胶体的持续硬化,因此,由胶体包裹的土颗粒与生成晶体之间的相互作用会更紧密,结构整体性及强度上升,抗风沙吹蚀能力增强。宏观的强度和力学性能和微观结构的关系仍需进一步的研究探索。
3. 结论
本文研究了利用黄原胶优化大豆脲酶诱导碳酸钙沉积防风固沙的办法,通过试验结果,可以得到以下结论:
(1)黄原胶能提高SICP固土所得表面强度,且表面强度随着大豆脲酶胶结液浓度和黄原胶浓度上升而上升,将0.3 mol/L的SICP胶结液与1 g/L黄原胶溶液混合,峰值贯入力可达138 N。黄原胶的加入对大豆脲酶催化水解尿素结合外加钙离子形成碳酸钙过程未表现出明显的抑制作用,可以认为在本研究试验参数范围内黄原胶不会抑制碳酸钙晶体的生成。
(2)试样的抗风能力随着SICP处理液和黄原胶的增加呈现增强趋势。只喷洒黄原胶溶液的土样,同等条件下的质量损失率也较大,将大豆脲酶胶结液和黄原胶共同作用后,质量损失率显著降低。试样风蚀质量损失率均随着风速增加、吹蚀时间增长而增长。同时发现,风蚀试验中,如果风夹带沙粒等固体颗粒,其侵蚀破坏能力明显更强。
(3)本文模拟沙漠地区风沙流吹蚀作用,发现黄原胶能增强SICP技术抗风沙侵蚀能力,随着大豆脲酶胶结液和黄原胶浓度提高,风蚀量降低,且与表面强度成正相关,当峰值贯入力达到70 N时,风蚀量在所有工况下均低于5%。
(4)由于胶体包裹、填充与黏结在晶体与土颗粒的空隙间,晶体与土颗粒之间的胶结更加紧密,土样的致密程度提高,表现为表面强度和抗风蚀强度的提高,这可能是黄原胶的加入优化碳酸钙与土颗粒之间胶结模式的原因。
-
表 1 数值模型参数表
Table 1 Parameters of numerical model
物性参数 取值 单位 初始水合物饱和度 0.443 — 初始水饱和度 0.351 — 储层孔隙度 0.182 — 储层绝对渗透率 96.7 mD 初始孔隙压力 3.75 MPa 降压边界压力 2.84 MPa 水合物地层温度 275.45 K 水的密度ρw 1000 kg/m3 甲烷天然气密度ρg 0.684 kg/m3 水合物密度ρh 917 kg/m3 水合物摩尔质量Mh 0.124 kg/mol 甲烷天然气摩尔质量Mg 0.016 kg/mol 水摩尔质量Mw 0.018 kg/mol 气相动力黏度μg 1.84×10 -5 Pa·s 液相动力黏度μw 1.01×10-5 Pa·s 液相残余饱和度Swr 0.1 — 气相残余饱和度Sgr 0.05 — VG模型参数m 0.45 — 初始毛细管力 1 kPa 表 2 敏感性分析参数表
Table 2 Parameters for sensitivity analysis
编号 影响因素 取值 单位 1 初始渗透率 40,60,80,100 mD 2 井筒降压大小 1.5,2.0,2.5,2.84 MPa -
[1] SUM A K, KOH C A, SLOAN E D. Clathrate hydrates: from laboratory science to engineering practice[J]. Industrial & Engineering Chemistry Research, 2009, 48(16): 7457–7465.
[2] 田慧会, 韦昌富, 颜荣涛, 等. 粉土中二氧化碳水合物分解过程的核磁试验研究[J]. 中国科学: 物理学力学天文学, 2019, 49(3): 173–180. doi: 10.3969/j.issn.0253-2778.2019.03.001 TIAN Hui-hui, WEI Chang-fu, YAN Rong-tao, et al. A NMR-based analysis of carbon dioxide hydrate dissociation process in silt[J]. Scientia Sinica (Physica, Mechanica & Astronomica), 2019, 49(3): 173–180. (in Chinese) doi: 10.3969/j.issn.0253-2778.2019.03.001
[3] DAI S, SEOL Y. Water permeability in hydrate-bearing sediments: A pore-scale study[J]. Geophysical Research Letters, 2014, 41(12): 4176–4184. doi: 10.1002/2014GL060535
[4] 刘乐乐, 张准, 宁伏龙, 等. 含水合物沉积物渗透率分形模型[J]. 中国科学: 物理学力学天文学, 2019, 49(3): 165–172. https://www.cnki.com.cn/Article/CJFDTOTAL-JGXK201903013.htm LIU Le-le, ZHANG Zhun, NING Fu-long, et al. A fractal model for the relative permeability prediction of hydrate-bearing sediments[J]. Scientia Sinica (Physica, Mechanica & Astronomica), 2019, 49(3): 165–172. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-JGXK201903013.htm
[5] SINGH H, MYSHAKIN E M, SEOL Y. A nonempirical relative permeability model for hydrate-bearing sediments[J]. Society of Petroleum Engineers Journal, 2019, 24(2): 547–562.
[6] 蔡建超, 夏宇轩, 徐赛, 等. 含水合物沉积物多相渗流特性研究进展[J]. 力学学报, 2020, 52(1): 208–223. https://www.cnki.com.cn/Article/CJFDTOTAL-LXXB202001019.htm CAI Jian-chao, XIA Yu-xuan, XU Sai, et al. Advances in multiphase seepage characteristics of natural gas hydrate sediments[J]. Chinese Journal of Theoretical and Applied Mechanics, 2020, 52(1): 208–223. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-LXXB202001019.htm
[7] KIM H C, BISHNOI P R, HEIDEMANN R A, et al. Kinetics of methane hydrate decomposition[J]. Chemical Engineering Science, 1987, 42(7): 1645–1653. doi: 10.1016/0009-2509(87)80169-0
[8] SUN X, NANCHARY N, MOHANTY K K. 1-D modeling of hydrate depressurization in porous media[J]. Transport in Porous Media, 2005, 58(3): 315–338. doi: 10.1007/s11242-004-1410-x
[9] KAMATH V A. A perspective on gas production from hydrates[C]// JNOC's Methane Hydrate International Symposium, 1998: 20–22.
[10] SUN X, LUO H, SOGA K. A coupled thermal-hydraulic- mechanical-chemical (THMC) model for methane hydrate bearing sediments using COMSOL Multiphysics[J]. Journal of Zhejiang University-SCIENCE A, 2018, 19(8): 600–623. doi: 10.1631/jzus.A1700464
[11] VAN GENUCHTEN M T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980, 44(5): 892–898. doi: 10.2136/sssaj1980.03615995004400050002x
[12] MASUDA Y. Numerical calculation of gas production performance from reservoirs containing natural gas hydrates[C]// Annual Technical Conference, San Antonio, Texas, 1997.
[13] HARDWICK J S, MATHIAS S A. Masuda's sandstone core hydrate dissociation experiment revisited[J]. Chemical Engineering Science, 2018, 175: 98–109. doi: 10.1016/j.ces.2017.09.003