基于 ANSYS Workbench 的高精度自动化三维地质建模方法——以天然气水合物相关的Slipstream海底滑坡为例

隆松伯1 何涛1,† 梁前勇2 蓝坤1 林进清2 董一飞2 何健2

1.造山带与地壳演化教育部重点实验室, 北京大学地球与空间科学学院, 北京 100871; 2.中国地质调查局广州海洋地质调查局, 广州 510760; †通信作者, E-mail: taohe@pku.edu.cn

摘要 在采用有限元软件 ANSYS 分析水合物相关的 Slipstream 海底滑坡体时, 针对多波束测深获得的复杂滑塌区海底地貌, 首先通过滑坡陡壁的形态恢复滑塌面被掩埋的部分, 然后通过与周围海脊的形态相似性对比, 重建滑塌发生前的海脊地貌。在此基础上, 利用 ANSYS Workbench 的 Jscript 脚本功能, 实现三维地质模型的高精度自动构建, 极大地提高了复杂坡体结构的建模效率, 可以为后续数值分析的准确性提供关键保障。对 Slipstream 海底滑塌坡体进行自重条件下稳定性的初步模拟, 结果显示, 以海底之下约 100m 的浅层水合物富集层为界, 其上的沉积层中最大剪应力呈现楔形分布的高值条带, 与该处滑坡后的阶梯状地貌吻合, 验证了地质模型的准确性和建模方法的有效性。

关键词 三维地质建模; ANSYS Workbench; 天然气水合物; 海底滑坡

天然气水合物是一种储量巨大、前景广阔的清洁燃料能源[12], 主要分布于大陆边缘陆坡区海水深度为几百米以上的沉积物中, 少部分赋存于大陆永久性冻土带[3]。天然气水合物与海底沉积层的稳定性有着密切的联系。与水合物相关的海底滑坡是海洋地质研究领域的一个热点问题。水合物层与上覆和下伏沉积物之间的物理性质差异会形成应力集中, 并且, 天然气水合物很容易因外界扰动而分解, 导致的沉积物液化和孔隙超压会大大地降低沉积层的力学稳定性。例如, Bouriak 等[4]和 Bugge 等[5]证实, 挪威外海 Storrega 区域 3 次主要滑塌事件中的第二次是由天然气水合物分解触发的。Booth 等[6]发现海底大多数滑坡分布在天然气水合物稳定带上限附近。但是, 水合物控制滑塌的具体机制目前尚不明确, 前人的研究要么限于现象描述, 要么受限于滑坡区三维地质建模的难度, 应力分布的定量分析模型比较简单和粗糙[7]

在地质过程数值模拟分析中, 地质建模是关键步骤之一, 而复杂地质构造的高精度三维建模则是难点。从20世纪80年代一些商业采矿软件包含基本的三维建模功能[8], 到1993年Houlding[9]提出三维地质建模概念, 再到今天各种建模软件和建模方法层出不穷, 地质建模技术经过三十余年的发展, 已基本上满足研究和工程的需求。但是, 地质空间关系的复杂性、软件功能的局限性和应用需求的独特性依然困扰着广大研究者, 限制了数值模拟方法的效果[10]。如何完善地质建模方法, 尤其是提高建模的速度、精度和自动化程度, 仍是地质模拟研究中的一个重要课题。

目前, 大型通用商业化有限元软件 ANSYS 在边坡稳定性分析中得到广泛应用。虽然 ANSYS 具有较强的规则形状建模功能, 但是面对复杂地貌的情况, 在创建高精度层面时, 其经典界面的人机交互建模方式效率非常低。由于应力应变分布与地质体形状密切相关, 因此许多研究者探索利用第三方的专业建模工具(如 AutoCAD)来创建复杂层面或实体模型[1112], 但在导入 ANSYS 的过程中存在兼容性风险和繁琐的二次修饰问题。针对这种情况, 在水合物相关的 Slipstream 海底滑坡研究中, 我们依托 ANSYS Workbench 内建的 Jscript 脚本功能, 开发一套自动化的复杂层面高精度建模方法, 高效地构建滑塌区复杂海底地貌的三维地质模型, 初步得到自重条件下的坡体内部应力分布, 加深了对水合物控制滑塌机制的理解。

1 Slipstream海底滑塌地貌重建

1.1 研究区地质背景

Slipstream 海底滑坡位于加拿大温哥华岛外卡斯卡迪(Cascadia)陆缘坡脚的变形前缘区。在此处, 胡安德富卡(Juan de Fuca)板块以约45 mm/a的速率俯冲于北美板块之下, 数百米厚的第四纪洋盆沉积物几乎完全从俯冲洋壳上刮削下来, 形成巨厚的增生楔[13]。这些增生楔相互挤压、缩短和重叠, 从而在陆坡底部形成众多与板块边界近平行的长条状背斜山脊, Slipstream 山脊就是其中之一(图 1)。该区域还拥有极为丰富的天然气水合物资源, 一直是海洋地质研究的焦点区域, 因此拥有丰富的地质、地球物理勘探资料。本研究采用的多波速测深、地震成像、地球物理测井等数据来自 SeaJade (Seafloor Earthquake Array — Japan Canada CascaDia Experi-ment)项目以及综合大洋钻探计划(Integrated Ocean Drilling Program, IODP)。

多波束测深得到的 Slipstream 滑坡现今海底地貌(图 1 和图 2(a))显示: 滑坡所在的变形前缘背斜山脊大部分坡体保持完整; 靠近背斜脊线的滑坡面暴露而陡直, 但滑坡面下部被大量滑坡残积物掩埋; 滑坡舌处地貌凌乱, 大量垮塌的坡体沉积物在此散落和堆积。

1.2 Slipstream滑坡前的山脊地貌重建

根据 Slipstream 滑坡区的地貌信息, 结合一般情况下滑坡面下部沿滑坡方向呈弧形、垂直于滑坡方向呈下凹谷状的特征[14]以及周围无滑坡山脊的形态, 可以大致推测滑坡面被掩埋部分的深度, 完成滑坡前坡体原始形态的地貌重建。多波束测深原始数据为经纬度坐标, 采样间隔为 0.00014°, 由于ANSYS 建模为直角坐标系, 因此转为通用的横墨卡托格网系统(universal transverse mercator, UTM), 得到大约10 m × 15 m的地貌网格。

我们截取的工作区域是以 Slipstream 山脊为中心的长 8km (X方向)、宽 4km (Y方向)的矩形。预测滑坡面下部被掩埋部分的形态时, 在沿滑坡方向(Y 方向)的截面上进行处理。以 X = −2000 m (图 2 (a)中红色曲线)的截面为例, 先提取此处水深数据, 画出滑坡后海底地貌的剖面图(图 2(b)中红色曲线), 然后依据滑坡面上部陡壁、下部凹面的普遍形态关系, 估计被掩埋的滑坡面位置(图 2(b)中蓝色曲线)。具体操作时, 在需要移除滑坡残积物的部分添加数个控制点(图 2(b)中紫色空心圆点), 通过这些控制点进行样条插值, 得到剖面曲线。最后预测得到的滑坡面如图 2(c)所示, 用于与周围山脊形态进行比较, 恢复原始 Slipstream 山脊形貌, 并与后续模拟的应力集中层位进行位置比较。

说明: width=425.15,height=464.85

据文献[14]修改。红色五角星为U1326测井位置; 红色三角形为04, 30, 31, 32, 33, 34和35号海底地震仪(Ocean Bottom Seismometer, OBS)的位置; 绿色实线为OBS测线1和4; 黄色实线为单道地震测线1, 4, 9, 11, 15, 16, 18和20; 叠加的等深线间隔25 m; 红色虚线框为三维地质建模区域; 红色箭头指示板块运动方向; 右上图中红色方框指示左下图的范围

图1 研究区构造背景、Slipstream滑坡区多波束海底地貌和地震勘探位置

Fig. 1 Tectonic background of study area, multibeam bathymetry of Slipstream submarine slide and seismic survey locations

我们的目标是模拟和分析滑坡发生前的临界应力状态, 研究水合物相关的滑坡发生机制, 因此需要恢复滑坡前山脊坡体的初始地貌形态。参考周围无滑坡山脊的近同心椭圆状等深线形态, 以预测得到的无掩埋滑坡面(图 2(c))为基础, 对其等深线数据(图 3(a))进行处理, 在坡体滑塌区域沿长轴方向(X 方向)设立控制点, 通过样条插值方法构建滑坡前的等深线, 最终的恢复结果如图 3(b)所示。图 2(d)展示重建的滑坡前 Slipstream 山脊地貌, 可以看到恢复的滑坡区与周围未滑坡坡体自然地融合, 地貌形态一致性较好, 可以保证后续三维地质建模的 质量。

说明: width=476.25,height=340.2

(a) 多波束测深获得的当前滑坡地貌; (b) 沿(a)中红色剖面进行的滑坡面预测; (c) 无残积体掩埋的滑坡面; (d) 恢复的滑坡前山脊形态

图2 滑坡面预测和原始地貌重建

Fig. 2 Glide plane estimation and original ridge morphology reconstruction

说明: width=444.7,height=154.55

图3 预测的滑坡面等深线(a)与恢复的滑坡前山脊等深线(b)对比

Fig. 3 Comparison of bathymetry contour maps for estimated glide plane (a) and pre-sliding ridge (b)

2 ANSYS三维地质建模

2.1 海底地貌建模

传统的 ANSYS 坡体稳定性分析方法采用经典界面, 曲面绘制能力有限, 无法自动地将海底地貌的 X, Y, Z 离散坐标点拟合成一个完整曲面, 在构建不规则起伏层面时, 需要采用人机交互方式, 先连点成样条曲线、再通过蒙皮(skin)操作绘制通过这些曲线的空间曲面。即使采用 ANSYS 的 APDL命令流进行参数化建模, 操作也很繁琐, 并且, 为了避免因曲面生成时间过长和生成的面太过破碎导致后面的布尔操作失败, 样条曲线的控制点数不宜超过 6 个, 生成面的曲线每次不宜超过 9 条。可见,在高精度的大量空间数据条件下, 这种方法非常费时费力, 以至难以完成建模。如果用专业建模软件生成曲面, 再导入 ANSYS, 则往往存在兼容性风险和二次修饰问题, 并且无法修改, 导致模型后期维护僵化。

ANSYS Workbench里的Model Designer模块具有类似 AutoCAD 的三维建模能力, 并可以通过Jscript 脚本文件, 自动执行建模指令。我们以此为基础, 发展了一套用 Matlab 编制建模脚本文件的方法, 实现复杂地质层面高精度建模过程的自动化和参数化, 从而脱离庞杂繁琐的人机交互建模工作。下面, 以 Slipstream 山脊的高精度三维地质建模为例进行说明。模型的顶面为恢复的滑坡前海底地形(图 2(d)), 空间网格大小为 10m × 10m, 共有 321201个坐标点。考虑到模型顶面(海底)的高程在−2106~ −2565 m之间, 水合物埋藏深度小于−3000 m, 将模型底面定为深度固定在−3000 m的水平面。

与通常的专业建模软件相似, Model Designer的建模方式如下: 在活动平面(Active Plane)上由线(Line)构成线条草图(Sketch), 再将线条草图包围的面积创建为剖面(Plane), 最后对一系列剖面进行蒙皮操作形成实体模型(Body)。因此, 建模所用的 Jscript 脚本命令首先是在 planeiSketchesOnly(p)函数中获得活动平面的坐标系位置, 创建草图 p.Ski(p 代表构建草图的剖面, i代表草图编号), 命名为Sketchi。如图 4(a)所示, 模型剖面的线条草图包含3 条直线(底边 p.Lni2, 两个侧边 p.Lni1 和 p.Lni3)和 1 条顶部曲线 p.Spi (其中 p.Spi是调用海底地貌深度数据, 用样条曲线方式生成), 再通过尺寸约束, 将这 4 条边围成模型在对应坐标位置的线条草图 Sketchi。接下来, 调用 planeiSketchesOnly (new Object())函数命令, 将创建的剖面 planei赋值给变量 psi。创建线条草图 Sketch1, 并最终生成剖面plane1的代码举例如下。

function plane1SketchesOnly (p)

{

//获得活动平面的坐标系位置

p.Plane = agb.GetActivePlane();

p.Origin = p.Plane.GetOrigin();

p.XAxis = p.Plane.GetXAxis();

p.YAxis = p.Plane.GetYAxis();

//新建一个草图并命名为Sketch1

p.Sk1 = p.Plane.newSketch();

p.Sk1.Name = "Sketch1";

//调用对应坐标位置的地貌数据并生成四条边

with (p.Sk1)

{

p.Ln11 = Line(0, −2544.1743, 0, −3000);

p.Ln12 = Line(0, −3000, 4000, −3000);

p.Ln13 = Line(4000, −3000, 4000, −2564.6686);

p.Sp1 = SplineBegin();

//用样条曲线生成顶部曲线

with(p.Sp1)

{

SplineFlexibility = agc.Yes;

说明: width=340.2,height=113.4

(a)中字母p代表剖面(plane), i为剖面的编号, Ln代表直线, Sp代表样条曲线; (b)中红色箭头表示将沿与剖面垂直方向创建的一系列剖面, 通过蒙皮处理形成三维实体模型

图4 由一系列草图剖面(a)构成三维模型(b)

Fig. 4 3D modeling using sequential sketch planes

SplineXY(950, −2544.1743);

SplineXY(940, −2544.8855);

......

SplineXY(0, −2564.6686);

SplineFitPtEnd();

}

}

//通过尺寸约束将四条边围成Sketch1

with (p.Plane)

{

HorizontalCon(p.Ln12);

VerticalCon(p.Ln11);

VerticalCon(p.Ln13);

CoincidentCon(p.Ln11.Base, 0, −2544.1743,

p.Sp1.Base, 0, −2544.1743);

CoincidentCon(p.Ln13.End, 4000, −2564.6686,

p.Sp1.End, 4000, −2564.6686);

CoincidentCon(p.Ln11.End, 0, −3000,

p.Ln12.Base, 0, −3000);

CoincidentCon(p.Ln12.End, 4000, −3000,

p.Ln13.Base, 4000, −3000);

}

p.Plane.EvalDimCons();

//返回由线条草图Sketch1生成的剖面p

return p;

}

//创建剖面plane1并赋值给变量ps1

var plane1 = agb.GetXYPlane();

agb.SetActivePlane (plane1);

var ps1 = plane1SketchesOnly (new Object())。

然后, 在前一个剖面(例如 planei)的位置基础上, 沿着坐标 Z方向移动一个网格距离 10m(图4(b)中红色箭头方向), 重复上述方法, 创建下一个剖面 planei+1, 并赋值给对应的变量 psi+1。例如, plane2的代码如下。

//在剖面plane1的基础上创建plane2

var plane2 = agb.PlaneFromPlane(plane1);

//沿坐标Z方向移动10 m

plane2.AddTransform(agc.XformZOffset, 10);

agb.regen();

agb.SetActivePlane (plane2);

//赋值给变量ps2

var ps2 = plane2SketchesOnly (new Object())。

在逐次完成所有剖面的创建后, 就可以用蒙皮操作创建三维实体。为了提高蒙皮的成功率, 对于地形起伏太大的高精度建模, 可以将模型分割为数行或者数列, 建立较小的局部剖面, 经过蒙皮操作之后, 会自动地融合成一个完整的实体模型。例如, 将 101 个剖面经蒙皮操作成为一个实体 Skin1 的代码如下。

var Skin1 = agb.Skin(agc.Add, agc.No, 0.0, 0.0);

Skin1.Name = "Skin1"

Skin1.AddBaseObject(ps1.Sk1);

Skin1.AddBaseObject(ps2.Sk2);

……

Skin1.AddBaseObject(ps101.Sk101);

agb.Regen()。

可见, 上述建模过程主要是循环地调用地貌数据, 因此我们编写了 Matlab 工具, 可以按照研究需求自动产生包含建模命令的 Jscript 脚本文件, 然后在 Model Designer 模块中运行该脚本文件即可自动完成 Slipstream 山脊的高精度三维实体模型, 无需人工干预, 大大减少了工作量。

经蒙皮操作后, 相邻剖面线条间最后会生成一个面。为方便后期施加模型载荷和约束, 可以进一步将这些面融合(merge)成一个面。因此, 最终生成的实体模型只含6个面(顶、底和4个侧面)。

2.2 坡体内部速度结构建模

沿着 Slipstream 山脊长轴中线的地震成像剖面(图 1 中测线 4)显示, 坡体内部存在两个主要的强 地震反射, 分别对应浅部深度约为 100m(meters below seafloor)的水合物富集砂体 GHS (gas hydrate sand)和深度为 265~275m 的水合物似海底反射层BSR(bottom simulating reflector), 因而可以大致将模型分为 5 层(图 5)。Yelisetti等[15]根据 OBS 走时数据, 反演出该处的纵波速度结构, 从而得到各层的厚度信息(表 1)。进一步反演出该处的横波速度结构, 并估算层 2 (GHS)和层4 (BSR)的水合物饱和度约为40% (何涛等未发表资料)。根据反演得出的纵横波速度结构和邻近的 IODP U1326(图 1 左上角红色五角星处)钻井获得的测井数据, 可以确定后续应力分布模拟计算所需的各层沉积物平均密度和泊松比(表1)。

当数据质量高并且信息充分时, 对坡体内部的沉积层, 可以采取与上述地貌建模相同的方法进行高精度的层面建模。由于本研究区缺乏高精度的三维地震数据, 仅有 5 条测线穿过 Slipstream 山脊(图1), 因此采用简化方式对坡体内部 5 个沉积层建模, 即利用 GHS 和 BSR 这两个水合物层大致与海底地貌平行的关系, 通过提取并且平移坡体模型的上表面到指定深度的位置, 再沿深度方向拉伸(thin)到相应厚度的方法来构建层 2 和层 4。最后, 使用布尔(boolean)操作将层 2 和层 4 插入坡体模型中, 分割出层1、层3和层5。建好的模型如图6所示。

说明: width=391.2,height=255.1

4 条彩色粗线条为地震剖面上拾取的强反射事件界面, 分别指示浅部高速层(层 2)和 BSR 层(层 4)的顶面和底面, 从而将浅部约200 m的沉积层按照速度结构从海底开始向下划分为5层

图5 沿Slipstream山脊长轴中线的地震测线4成像结果

Fig. 5 Migrated seismic section of single channel line-4 along the long axis of Slipstream

表1 地层物性特征参数

Table 1 Physical parameters of each layer

层序号厚度/m水合物饱和度/%平均密度/ (g∙cm−3)泊松比纵波速度/ (m∙s−1) 1100 02.020.4871530 2 30402.040.4851580+3.9d 3135 01.930.4601260+3.3d 4 15402.090.4302090 5120 02.210.4201100+2.0d

说明: d为沉积层相对于海底的深度。

2.3 有限元网格剖分

模型的网格剖分(mesh)是有限元数值模拟成功与否的关键步骤之一。本文的三维模型使用 ANSYS Workbench 的多域扫掠网格剖分方式, 即将模型自动地分成多个规则区域, 然后对每一个区域进行扫掠网格剖分。该方法以具有 20 个节点(8 个顶点加上 12 条边的中点)的高次六面体单元为主进行自由剖分, 每个节点都可以进行求解, 并且具有在地形复杂区域自动加密单元和退化为四面体、棱锥、棱柱等单元的能力。经过收敛测试和比较, 本文模型最后剖分为 44055 个单元(总节点数为 242084), 其中最小的单元边长约为 15m, 尺度与地貌数据的UTM网格相当, 达到精度要求。

说明: width=221.15,height=93.6

红色数字为与图5对应的模型层号

图6 Slipstream坡体三维地质模型

Fig. 6 Slipstream slope 3D geological model

说明: width=340.2,height=189.95

(a)沿滑塌区中线的最大剪应力分布; (b)为(a)中红色虚线框的放大部分, 红色实线代表现在海底地貌, 红色虚线代表预测的滑坡面位置, 黄线划分楔形的应力高值带并指示潜在破裂面

图7 最大剪应力分布

Fig. 7 Maximum shear stress distribution

3 自重条件下的应力数值模拟

本文对 Slipstream 山脊初步开展自重条件下坡体内部应力分布的模型分析。由于滑坡前坡体内部各层之间没有相互滑动, 因此用 ANSYS 的接触分析方法将各沉积层接触面绑在一起, 使分离的 5 个沉积层组合为一个整体。考虑到 Slipstream 山脊的大部分坡体现今仍然保持完好, 因此对模型底面施加固定约束(fixed supports), 对四周侧面施加用于模拟对称边界的无摩擦约束(frictionless support)。最后, 添加标准重力加速度(standard earth gravity)载荷进行数值求解。

说明: width=221.15,height=221.15

(a) 浅部水合物富集层(橙色虚线)之上的最大剪应力高值条带将沉积层分割成类似楔形的块体, 坡脚处的块体在外界扰动时最容易失稳(红色箭头指示楔形块体滑动方向); (b) 坡脚沉积层滑走后, 后面紧邻的块体失去支撑而失稳滑动; (c) 沉积层块体向后逐次垮塌, 直到接近山脊处, 达到新的平衡; (d) 滑塌的沉积物一部分冲出坡脚散布在深海盆, 另外一部分残留在坡脚上方, 掩埋了下部的滑坡面

图8 浅部水合物富集层相关的低角度海底滑坡模型

Fig. 8 Sketch cartoon of low angle submarine slide along the shallow high hydrate concentration layer

X=−2000m(图 2(a)中红色曲线)处的滑坡区中线截面为例, 自由重力载荷下的最大剪应力分布如图 7 所示。最大剪应力分布的第一个特征是应力在浅部水合物富集层的上表面发生应力集中, 并且与推测的滑坡面深度大致相同(图 7(b)中红色虚线)。Yelisetti 等[15]已经发现滑坡面与浅水合物层具有空间一致性, 本文的结果进一步表明, 该处的高砂质含量和高水合物饱和度形成一个力学性能突变层, 导致应力集中而使滑坡面位于此深度。

最大剪应力分布的第二个特征是浅部水合物富集层之上的沉积层中应力高值条带呈楔形或三角形(图 7(b)中黄线)。Kvalstad 等[7]在研究挪威西海岸的 Storegga 滑塌时, 观察到滑坡后海底存在明显的阶梯状地貌, 与 Slipstream 滑塌区(图 2(a))相似。因此, 我们推测 Slipstream 滑塌的原因和过程大致如图 8 所示: 胡安德富卡板块在变形前缘俯冲时产生孔隙超压、大地震等破坏坡体稳定性的扰动因素, 诱发坡脚处沉积层最先失稳, 然后, 浅部水合物富集层造成的楔形剪应力集中条带导致坡体以楔形或三角形块体逐渐垮塌, 滑坡壁逐步后退, 直至垮塌接近山脊处后形成新的平衡。

4 结论

1)在原始地貌未被完全破坏的情况下, 可以通过地貌形态的对比分析, 逐步获得完整滑坡面和滑坡前的原始山脊坡体。

2)本研究编制的 Matlab 程序可以自动地将地形数据转换为 Jscript 建模脚本, 从而可在 ANSYS Workbench 中实现三维地质模型的自动构建, 使人机交互次数大幅度地减少, 极大地提高对复杂层面进行高精度建模的效率。结合布尔运算等操作, 可以构建如 Slipstream 滑坡区一样具复杂沉积层结构的精细三维地质模型, 为后续数值分析的准确性提供关键保障。

3)对 Slipstream 山脊模型初步进行自重条件下的稳定性模拟, 结果显示滑坡面与海底之下约 100m 的浅部水合物富集层大致处于同一深度, 最大剪应力在该层的上表面集中, 在该层之上的沉积层中呈现楔形或三角形的高值条带, 与滑坡后海底的阶梯状地貌相吻合。这些结果表明 Slipstream 滑坡与浅部水合物富集层之间存在密切联系, 也验证了本文地质建模方法的有效性。

参考文献

[1]Milkov A V. Global estimates of hydrate-bound gas in marine sediments: how much is really out there?. Earth-Science Reviews, 2004, 66(3/4): 183‒197

[2]Collett T S. Energy resource potential of natural gas hydrates. Aapg Bulletin, 2002, 86(11): 1971‒1992

[3]蒋向明. 天然气水合物的形成条件及成因分析. 中国煤炭地质, 2009, 21(12): 7‒11

[4]Bouriak S, Vanneste M, Saoutkine A. Inferred gas hydrates and clay diapirs near the Storegga Slide on the southern edge of the Vøring Plateau, offshore Norway. Marine Geology, 2000, 163: 125‒148

[5]Bugge T, Befring S, Belderson R H, et al. A giant three-stage submarine slide off Norway. Geo-Marine Letters, 1987, 7(4): 191‒198

[6]Booth J S, Winters W J, Dillon W P. Circumstantial evidence of gas hydrate and slope failure associations on the United States, Atlantic continental margin. Annals of the New York Academy of Sciences, 1994, 715: 487‒489

[7]Kvalstad T J, Andresen L, Forsberg C F, et al. The Storegga slide: evaluation of triggering sources and slide mechanics. Marine & Petroleum Geology, 2005, 22(1/2): 245‒256

[8]王明华, 白云. 三维地质建模研究现状与发展趋势. 土工基础, 2006, 20(4): 68‒70

[9]Houlding S W. 3D geoscience modeling: computer techniques for geological characterization. Berlin: Springer-Verlag, 1994

[10]武强, 徐华. 三维地质建模与可视化方法研究. 中国科学(D辑), 2004, 34(1): 54‒60

[11]郝钟雄. ANSYS 与 CAD 软件的接口问题研究. 机械设计与制造, 2007(7): 75‒76

[12]张晓东, 于擘, 王袖和. Ansys workbench参数传递辅助建立连续颌骨质量变化的种植体骨块三维有限元模型. 中国厂矿医学, 2008, 21(2): 174‒175

[13]Riedel M, Collett T S, Malone M J, et al. Proceedings of the integrated ocean drilling program expedition 311 [R]. Washington, DC: Integrated Ocean Drill Pro-gram, 2006

[14]姜德义. 边坡稳定性分析与滑坡防治. 重庆: 重庆大学出版社, 2005

[15]Yelisetti S, Spence G D, Riedel M. Role of gas hy-drates in slope failure on frontal ridge of northern Cascadia margin. Geophysical Journal International, 2014, 199(1): 441‒458

A High Precision Automatic 3D Geological Modeling Method Based on ANSYS Workbench: A Case Study of Gas Hydrate- related Slipstream Submarine Slide

LONG Songbo1, HE Tao1,†, LIANG Qianyong2, LAN Kun1, LIN Jinqing2, DONG Yifei2, HE Jian2

1.The Key Laboratory of Orogenic Belts and Crustal Evolution (MOE), School of Earth and Space Sciences, Peking University, Beijing 100871; 2. Guangzhou Marine Geological Survey, China Geological Survey, Guangzhou 510760; † Corresponding author, E-mail: taohe@pku.edu.cn

Abstract In the study of gas hydrate-related Slipstream submarine slide, the finite element analysis software ANSYS is used to construct the 3D model with complicated submarine slump topography acquired by multibeam sounding system. The lower part of sliding surface buried by slump accumulation is estimated from the main scarp geometry, and the original ground surface before slump is reconstructed according to the morphological similarity of surrounding ridges. Then, the high precision 3D geological model is automatically completed by running Jscript file in ANSYS Workbench, which greatly improves the efficiency of complex geometric modeling and thus provides a key guarantee for the accuracy of subsequent finite element numerical analysis. The stability simulation of Slipstream Ridge under its self-weight condition showes that the maximum shear stress in sediments above a shallow gas hydrate concentration layer at about 100 meters below seafloor is distributed as a series of high value bands in wedge shape, which matches well with the stepped topography observed on the current slump surface and verifies the accuracy of the 3D geological model and the validity of theproposed modeling method.

Key words 3D geological modeling; ANSYS Workbench; gas hydrate; submarine slide

中图分类号 P628

doi: 10.13209/j.0479-8023.2018.051

国家自然科学基金(40904029, 41676032)和中国地质调查局天然气水合物专项(DD20160217)资助

收稿日期: 2017−07−27;

修回日期: 2017−09−24;

网络出版日期: 2018−07−04