水射流破碎南海含水合物沉积物数值模拟研究

潘栋彬， 陈 晨,3， 杨 林， 靳成才， 李曦桐， 刘书源， 陈 勇， 朱 江

(1.吉林大学建设工程学院，吉林 长春 130026； 2.自然资源部复杂条件钻采技术重点实验室，吉林 长春 130026； 3.超硬材料国家重点实验室，吉林 长春 130012)

0 引言

天然气水合物(Natural Gas Hydrates，简称NGH)是在一定条件下由轻烃、二氧化碳及硫化氢等小分子气体与水相互作用形成的白色固态结晶物质[1]。天然气水合物具有能量密度大、一定的形成条件、储量大和分布广四个主要特点[2]。天然气水合物一般形成于高压低温环境，单位体积的天然气水合物分解后能够产生164 m3的天然气与0.8 m3的水，能量密度远大于石油、煤炭、天然气等常规能源；全球天然气水合物总量大约达到7.6×1018m3，其主要分布于海底与陆地永久冻土区，全球的水合物中含有机碳总量相当于已探明的常规化石有机碳储量的2倍[3-4]。由于天然气水合物的资源潜力巨大，近几十年在世界范围内兴起了天然气水合物勘查、试采、开发的研究热潮，美国、日本、加拿大、印度与韩国等国已经制定了针对天然气水合物的长期研究计划，如何安全高效以及环保地开采天然气水合物资源已经成为世界各国关注的前沿与焦点[5-6]。

继2007年与2013年在南海钻获天然气水合物实物样品以来[7-8]，我国于2017年在南海神狐海域成功试采天然气水合物[9]，这展示了我国南海巨大的天然气水合物资源开发潜力。我国南海水合物有着赋存温度低、压力高、埋藏深度浅、没有致密盖层、储层胶结强度低、易于碎化等特点。与传统油气藏与砂岩水合物储层相比，我国南海水合物储层的骨架不具有岩石构架，水合物本身就是储层骨架的一部分，储层稳定性不足，水合物分解后容易导致储层出现失稳垮塌，进而引发海底滑坡等地质灾害。目前一般用热激发法、降压法、注化学试剂法与CO2置换法等开采方法对海洋天然气水合物进行开采[10-11]。但是单纯地采用这四种开采方法局限性较大：降压法可能引起水合物的二次生成或者冰的生成，堵塞渗透路径，于长期开采不利；热激法具有热利用效率较低的问题，且只能进行小范围加热；化学试剂法的缺点是化学试剂昂贵，对天然气水合物层的作用缓慢，存在环境污染问题；CO2置换法开采周期长，要求水合物储层具有较高的渗透性。因此，针对我国此类非成岩天然气水合物资源，采取安全有效的方法对此类水合物资源进行开采有极大的必要性。

目前，高压水射流方法在传统能源开采方面已得到广泛应用，技术较为成熟[12-16]，其具有无污染、能量利用率高、低热、反作用小等特点，且能够大幅度提高采矿速率、降低成本。高文爽、陈晨等[11]提出了采用高压热射流的方法开采天然气水合物的思路，高压热射流不仅可以将天然气水合物切削成便于水力输送的细小颗粒，还能够不断地热激水合物使其分解，热能利用率高。周守为等人提出了水合物射流破碎固态流化开采技术，同时，在2017年5月该技术成功应用于我国南海神狐海域天然气水合物的试采工程中[17]。这说明射流开采天然气水合物是完全可行的，我国南海水合物弱胶结、非成岩、渗透性低的赋存特点能够满足高压水射流破碎的条件。采用高压水射流切割、破碎南海天然气水合物储层是一种研究潜力巨大的新开采方法，具有重大的研究价值与意义。

水射流切割、破碎含水合物沉积物是高压射流开采水合物方法的核心步骤之一，因此探究射流破碎天然气水合物过程与机理则有极大必要性。在淹没状态下，水射流破碎含水合物沉积物过程属于流体和固体的非线性碰撞动力耦合问题，即流固耦合问题。LS-DYNA是一种基于显式积分的有限元软件，能够采用ALE数值计算方法解决流固耦合问题，其有效性受到国内外学者公认。因此，本文采用著名的LS-DYNA显示动力分析有限元程序，对淹没状态下，水射流破碎南海海底含水合物沉积物过程进行数值模拟研究，研究了不同射流速度对高压水射流作用下含水合物沉积物破碎效果的影响规律。

1 ALE算法理论基础

1.1 ALE算法基本思想

用于计算的网格单元不会附着于流体质点，网格的位置能够在空间坐标系中进行移动，且其运动状态能够进行自由选定。ALE算法兼具Lagrange和Euler算法两者的优点，即在物质运动边界的处理上吸收了Lagrange算法的特点，因而可以对固体结构边界的运动进行准确描述；同时也吸收了Euler算法的特点，网格单元与物质本身相互独立，网格单元之间的物质能够进行流动，且网格可以在数值计算求解时根据所需的参数进行位置调整，防止对计算网格造成严重的畸变。

在ALE方法中，已知参考构型，因此在描述各物理量时需引入一个参考坐标y，即F=F(y,t)。各个物理量的物质导数被转换为对应物理量的参考导数：

(1)

式中：Xi——拉格朗日坐标系；yi——欧拉坐标系；wi——相对速度，wi=u-v；u——物质运动速度；v——空间网格运动速度。

1.2 ALE算法的流体控制方程

将式(1)中的控制方程组进行坐标转换，转换到参考坐标系：

(1)质量守恒方程

(2)

(2)动量守恒方程

(3)

(3)能量守恒方程

(4)

式中：ρ——物质密度；fi——单位质量体力；σij——柯西应力张量；e——单位质量的内能；qi——热通量；i、j——坐标轴方向，(i=1,2,3；j=1,2,3)。

2 材料模型

2.1 水的材料模式

本文采用MAT_NULL本构模型，为了给材料补充完整的应力张量，MAT_NULL本构模型要与材料的状态方程联立，共同对流体行为进行描述，其中状态方程可定义压缩材料与膨胀材料的压力。本文水介质状态方程采用EOS_GRUNEISEN，该状态方程表达式为：

(γ0+aμ)E

(5)

式中：ρ0——流体初始密度；μ——空气的比体积；C——冲击波速度与质点速度up关系曲线的纵轴截距；γ0——Gruneisen系数；a——一阶体积修正系数；S1、S2、S3——关系曲线的斜率系数；E——单位体积初始内能。

水的参数为[18]：密度1.05 g/cm3，压力截断值-10 Pa，动力粘性系数8.50×103Pa·μs，us-up曲线截距C=0.148 cm/us，Gruneisen伽马γ0=0.5。

2.2 含水合物沉积物材料模式

*MAT_FHWA_SOIL是一种常用于模拟土壤的本构模型，属于非线性弹塑性各向同性损伤材料模型。土体可看作是离散颗粒相互粘结而形成的集合体[19]，粘聚力较小，它与摩擦强度共同决定了土体强度。本文数值模拟分析的研究对象为海底含水合物沉积物，水合物以分散状、填隙状、脉状及层状等状态赋存于砂质、粘土质粉砂和粘土沉积物中，水合物与沉积物通过胶结作用形成统一整体[20]。与原始沉积物相比，含水合物沉积物具有更大的粘聚力与内摩擦角，固结强度更高。因此，本文将含水合物沉积物模型简化为力学性能得到改善的土体，并采用*MAT_FHWA_SOIL本构模型对含水合物沉积物进行描述。含水合物沉积物材料参数为[21-23]：密度2.0 g/cm3，比重2.7，体积模量9.3 GPa，剪切模量4.2 GPa，内摩擦角0.436 rad，内聚力45 kPa，含水率50%。

3 数值仿真模型

3.1 计算模型建立

含水合物沉积物整体几何模型尺寸为4 cm×4 cm×10 cm；水域整体几何模型尺寸为5 cm×5 cm×11 cm(水域模型高度视靶距而定，保证沉积物底面与水域底面距离为0.5 cm)；射流源为圆柱状，半径为0.1 cm(大小视实验喷嘴半径而定)，高度为0.05 cm；射流破碎含水合物沉积物模型是轴对称的，为了减少计算量，提高计算精度与效率，因此只针对模型的1/4进行模拟计算分析，淹没射流破碎含水合物沉积物三维模型(1/4模型)如图1所示。

3.2 网格划分

采用映射网格法对射流破碎含水合物几何模型进行网格划分，为了在提高计算效率的同时保证计算精度，对几何模型进行网格局部加密。水射流破碎含水合物沉积物时存在明显的局部效应，因此，对水域几何模型中心处的附近区域进行网格加密。图2中(a)与(b)分别显示出了水域和含水合物沉积物网格划分情况。

3.3 边界条件

对1/4模型的2个对称面施加法线方向的位移约束，即限制模型YOZ对称面上X轴方向位移与限制XOZ对称面上Y轴方向位移；同时，限制含水合物沉积物模型底面的所有自由度，保证水射流冲击过程中含水合物沉积物模型不发生运动。同时，为模拟半无限大的含水合物沉积物边界，将1/4模型中含水合物沉积物模型的底面与两个侧面设置为无反射边界，同时也将水域模型底面与两个侧面的设置为无反射边界。

图2 水域与含水合物沉积物网格局部加密图

4 数值模拟结果和讨论

4.1 射流速度对冲蚀深度的影响规律

图3所示的为射流速度与冲蚀体积关系图，从图3中可以看出，随着射流速度的增大，冲蚀深度逐渐增大，两者呈线性递增关系，这与Lu等[24]人的射流破岩试验结果一致；同时，在同一射流速度条件下，喷距越小，冲蚀深度越大。分析认为，射流速度越大，射流所具有的能量越大，在同一条件下，射流冲蚀含水合物沉积物能够到达的深度越大；射流速度一定，喷距越小，射流能量耗散得越少，冲蚀深度越大。在分析射流速度与岩石冲蚀深度关系时，Momber将射流破岩速度分为潜伏段与侵蚀段，当射流达到阈值速度时岩石才能够开始被破碎。射流冲蚀含水合物沉积物与射流破岩类似，由图3、图4可知，当射流速度小于一定值时，冲蚀深度与冲蚀体积为零，射流无法对含水合物沉积物产生破碎。

4.2 射流速度对冲蚀体积的影响规律

图4所示的为射流速度与冲蚀体积关系图，由图4可以看出，随着射流速度的增大，含水合物沉积物冲蚀体积逐渐增大，当射流速度较小时，曲线斜率较低，冲蚀体积增大得较慢；当射流速度较大时，曲线斜率较大，冲蚀体积增长得较快，这与Momber[25]的射流破岩试验结果一致。含水合物沉积物冲蚀体积是含水合物沉积物的轴向冲蚀与径向冲蚀共同作用的结果。水射流处于淹没状态，射流在运动的过程中卷吸周围的水，射流冲蚀流量逐渐增加。当其他条件一定时，射流速度越大，水射流具有的能量越大，射流在冲蚀含水合物沉积物的同时能够卷吸更多的水，冲蚀有效面积越大，对含水合物沉积物的轴向与径向冲蚀作用增强，轴向冲蚀体积与径向冲蚀体积迅速增大，加大了含水合物沉积物冲蚀体积递增速率。在实际海底射流冲蚀开采水合物的过程中，可以在设备条件允许的情况下提高射流速度，提高开采效率。

图3 射流速度与冲蚀深度关系图

图4 射流速度与冲蚀体积关系图

5 结论

(1)随着射流速度的增大，冲蚀深度逐渐增大，两者呈线性递增关系。水射流破碎含水合物沉积物存在临界破碎压力，当水射流压力大于临界破碎压力时，才能对含水合物沉积物进行破碎。

(2)含水合物沉积物冲蚀体积是轴向冲蚀与径向冲蚀共同作用的结果，射流速度越大，对含水合物沉积物的轴向与径向冲蚀作用增强，加大了含水合物沉积物冲蚀体积递增速率。在实际海底射流冲蚀开采水合物的过程中，可以在设备条件允许的情况下提高射流速度，提高开采效率。

(3)射流开采天然气水合物是完全可行的，我国南海水合物弱胶结、非成岩、渗透性低的赋存特点能够满足高压水射流破碎的条件，采用高压水射流切割、破碎南海天然气水合物储层的是一种研究潜力巨大的新开采方法，具有重大的研究价值与意义。