ACTA Scientiarum Naturalium Universitatis Pekinensis
A High Precision Automatic 3D Geological Modeling Method Based on ANSYS Workbench: A Case Study of Gas Hydraterelated Slipstream Submarine Slide
LONG Songbo, HE Tao, LIANG Qianyong, et al
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 the proposed modeling method. Key words 3D geological modeling; ANSYS Workbench; gas hydrate; submarine slide
天然气水合物是一种储量巨大、前景广阔的清洁燃料能源[1–2], 主要分布于大陆边缘陆坡区海水深度为几百米以上的沉积物中, 少部分赋存于大陆永久性冻土带[3]。天然气水合物与海底沉积层的稳定性有着密切的联系。与水合物相关的海底滑坡是海洋地质研究领域的一个热点问题。水合物层与上覆和下伏沉积物之间的物理性质差异会形成应力集中, 并且, 天然气水合物很容易因外界扰动而分解,导致的沉积物液化和孔隙超压会大大地降低沉积层的力学稳定性。例如, Bouriak 等[4]和 Bugge 等[5]证实, 挪威外海 Storrega 区域3次主要滑塌事件中的
[6]第二次是由天然气水合物分解触发的。Booth等发现海底大多数滑坡分布在天然气水合物稳定带上限附近。但是, 水合物控制滑塌的具体机制目前尚不明确, 前人的研究要么限于现象描述, 要么受限于滑坡区三维地质建模的难度, 应力分布的定量分析模型比较简单和粗糙[7]。
在地质过程数值模拟分析中, 地质建模是关键步骤之一, 而复杂地质构造的高精度三维建模则是难点。从 20 世纪 80 年代一些商业采矿软件包含基本的三维建模功能[8], 到 1993 年 Houlding[9]提出三维地质建模概念, 再到今天各种建模软件和建模方法层出不穷, 地质建模技术经过三十余年的发展,已基本上满足研究和工程的需求。但是, 地质空间关系的复杂性、软件功能的局限性和应用需求的独特性依然困扰着广大研究者, 限制了数值模拟方法的效果[10]。如何完善地质建模方法, 尤其是提高建模的速度、精度和自动化程度, 仍是地质模拟研究中的一个重要课题。
目前, 大型通用商业化有限元软件ANSYS在边坡稳定性分析中得到广泛应用。虽然ANSYS具有较强的规则形状建模功能, 但是面对复杂地貌的情况, 在创建高精度层面时, 其经典界面的人机交互建模方式效率非常低。由于应力应变分布与地质体形状密切相关, 因此许多研究者探索利用第三方的专业建模工具(如 AUTOCAD)来创建复杂层面或实体模型[11–12],但在导入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 Experiment)项目以及综合大洋钻探计划(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 山脊为中心的长8 km (X 方向)、宽4 km (Y 方向)的矩形。预测滑坡面下部被掩埋部分的形态时, 在沿滑坡方向(Y 方向)的截面上进行处理。以 X = −2000 m (图 2 (a)中红色曲线)的截面为例, 先提取此处水深数据,画出滑坡后海底地貌的剖面图(图 2(b)中红色曲线),然后依据滑坡面上部陡壁、下部凹面的普遍形态关系, 估计被掩埋的滑坡面位置(图 2(b)中蓝色曲线)。
具体操作时, 在需要移除滑坡残积物的部分添加数个控制点(图2(b)中紫色空心圆点), 通过这些控制点进行样条插值, 得到剖面曲线。最后预测得到的滑坡面如图2(c)所示, 用于与周围山脊形态进行比较, 恢复原始 Slipstream 山脊形貌, 并与后续模拟的应力集中层位进行位置比较。我们的目标是模拟和分析滑坡发生前的临界应 力状态, 研究水合物相关的滑坡发生机制, 因此需要恢复滑坡前山脊坡体的初始地貌形态。参考周围无滑坡山脊的近同心椭圆状等深线形态, 以预测得到的无掩埋滑坡面(图 2(c))为基础, 对其等深线数据(图 3(a))进行处理, 在坡体滑塌区域沿长轴方向(X方向)设立控制点, 通过样条插值方法构建滑坡前的等深线, 最终的恢复结果如图 3(b)所示。图 2
(d)展示重建的滑坡前 Slipstream 山脊地貌, 可以看到恢复的滑坡区与周围未滑坡坡体自然地融合, 地貌形态一致性较好, 可以保证后续三维地质建模的质量。
2 ANSYS 三维地质建模2.1 海底地貌建模
传统的ANSYS坡体稳定性分析方法采用经典界面, 曲面绘制能力有限, 无法自动地将海底地貌
的 X, Y, Z 离散坐标点拟合成一个完整曲面, 在构建不规则起伏层面时, 需要采用人机交互方式, 先连点成样条曲线、再通过蒙皮(skin)操作绘制通过这些曲线的空间曲面。即使采用ANSYS 的 APDL命令流进行参数化建模, 操作也很繁琐, 并且, 为了避免因曲面生成时间过长和生成的面太过破碎导致后面的布尔操作失败, 样条曲线的控制点数不宜超过6 个, 生成面的曲线每次不宜超过9条。可见,在高精度的大量空间数据条件下, 这种方法非常费时费力, 以至难以完成建模。如果用专业建模软件生成曲面, 再导入ANSYS, 则往往存在兼容性风险和二次修饰问题, 并且无法修改, 导致模型后期维护僵化。
ANSYS Workbench 里的 Model Designer 模块具有类似 AUTOCAD的三维建模能力, 并可以通过Jscript 脚本文件, 自动执行建模指令。我们以此为基础, 发展了一套用Matlab编制建模脚本文件的方法, 实现复杂地质层面高精度建模过程的自动化和参数化, 从而脱离庞杂繁琐的人机交互建模工作。下面, 以 Slipstream山脊的高精度三维地质建模为例进行说明。模型的顶面为恢复的滑坡前海底地形(图 2(d)), 空间网格大小为10 m × 10 m, 共有 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;
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方向移动一个网格距离10 m (图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)显示, 坡体内部存在两个主要的强地震反射, 分别对应浅部深度约为100 m (meters below seafloor)的水合物富集砂体 GHS (gas hydrate sand)和深度为 265~275 m的水合物似海底反射层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 所示。
2.3 有限元网格剖分
模型的网格剖分(mesh)是有限元数值模拟成功
与否的关键步骤之一。本文的三维模型使用ANSYS Workbench的多域扫掠网格剖分方式, 即将模型自动地分成多个规则区域, 然后对每一个区域进行扫掠网格剖分。该方法以具有20个节点(8个顶点加上 12条边的中点)的高次六面体单元为主进行自由剖分, 每个节点都可以进行求解, 并且具有在地形复杂区域自动加密单元和退化为四面体、棱锥、棱柱等单元的能力。经过收敛测试和比较, 本文模型最后剖分为44055个单元(总节点数为 242084), 其中最小的单元边长约为15 m, 尺度与地貌数据的UTM 网格相当, 达到精度要求。
3 自重条件下的应力数值模拟
本文对 Slipstream山脊初步开展自重条件下坡体内部应力分布的模型分析。由于滑坡前坡体内部各层之间没有相互滑动, 因此用ANSYS的接触分析方法将各沉积层接触面绑在一起, 使分离的 5 个沉积层组合为一个整体。考虑到Slipstream 山脊的大部分坡体现今仍然保持完好, 因此对模型底面施加固定约束(fixed supports), 对四周侧面施加用于模拟对称边界的无摩擦约束(frictionless support)。最后, 添加标准重力加速度(standard earth gravity)载荷进行数值求解。
以 X = −2000 m (图 2(a)中红色曲线)处的滑坡区中线截面为例, 自由重力载荷下的最大剪应力分布如图7所示。最大剪应力分布的第一个特征是应力在浅部水合物富集层的上表面发生应力集中, 并且与推测的滑坡面深度大致相同(图 7(b)中红色虚
[15]线)。Yelisetti 等 已经发现滑坡面与浅水合物层具有空间一致性, 本文的结果进一步表明, 该处的高砂质含量和高水合物饱和度形成一个力学性能突变层, 导致应力集中而使滑坡面位于此深度。最大剪应力分布的第二个特征是浅部水合物富集层之上的沉积层中应力高值条带呈楔形或三角形
[7] (图 7(b)中黄线)。Kvalstad 等 在研究挪威西海岸的 Storegga 滑塌时, 观察到滑坡后海底存在明显的阶梯状地貌, 与 Slipstream 滑塌区(图2(a))相似。因 此, 我们推测 Slipstream 滑塌的原因和过程大致如图 8 所示: 胡安德富卡板块在变形前缘俯冲时产生孔隙超压、大地震等破坏坡体稳定性的扰动因素,
诱发坡脚处沉积层最先失稳, 然后, 浅部水合物富集层造成的楔形剪应力集中条带导致坡体以楔形或三角形块体逐渐垮塌, 滑坡壁逐步后退, 直至垮塌接近山脊处后形成新的平衡。
4 结论
1) 在原始地貌未被完全破坏的情况下, 可以通过地貌形态的对比分析, 逐步获得完整滑坡面和滑坡前的原始山脊坡体。
2) 本研究编制的 Matlab程序可以自动地将地形数据转换为 Jscript 建模脚本, 从而可在ANSYS Workbench中实现三维地质模型的自动构建, 使人机交互次数大幅度地减少, 极大地提高对复杂层面进行高精度建模的效率。结合布尔运算等操作, 可以构建如 Slipstream滑坡区一样具复杂沉积层结构的精细三维地质模型, 为后续数值分析的准确性提供关键保障。
3) 对 Slipstream山脊模型初步进行自重条件下的稳定性模拟, 结果显示滑坡面与海底之下约100 m的浅部水合物富集层大致处于同一深度, 最大剪应力在该层的上表面集中, 在该层之上的沉积层中呈现楔形或三角形的高值条带, 与滑坡后海底的阶梯状地貌相吻合。这些结果表明Slipstream 滑坡与浅部水合物富集层之间存在密切联系, 也验证了本文地质建模方法的有效性。 [1] Milkov A V. Global estimates of hydrate-bound gas
183‒197 in marine sediments: how much is really out there?. Earth-science Reviews, 2004, 66(3/4):
1971‒1992 [2] Collett T S. Energy resource potential of natural gas hydrates. Aapg Bulletin, 2002, 86(11):
7‒11 [3] 蒋向明. 天然气水合物的形成条件及成因分析. 中国煤炭地质, 2009, 21(12): [4] Bouriak S, Vanneste M, Saoutkine A. Inferred gas hydrates and clay diapirs near the Storegga Slide on
125‒148 the southern edge of the Vøring Plateau, offshore Norway. Marine Geology, 2000, 163: [5] Bugge T, Befring S, Belderson R H, et al. A giant
191‒198 three-stage submarine slide off Norway. Geo-marine Letters, 1987, 7(4): [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.
487‒489 Annals of the New York Academy of Sciences, 1994, 715: [7] Kvalstad T J, Andresen L, Forsberg C F, et al. The Storegga slide: evaluation of triggering sources and
245‒256 slide mechanics. Marine & Petroleum Geology, 2005, 22(1/2):
68‒70 [8] 王明华, 白云. 三维地质建模研究现状与发展趋势.土工基础, 2006, 20(4): [9] Houlding S W. 3D geoscience modeling: computer techniques for geological characterization. Berlin: Springer-verlag, 1994
54‒60 [10] 武强, 徐华. 三维地质建模与可视化方法研究. 中国科学(D 辑), 2004, 34(1):
75‒76 [11] 郝钟雄. ANSYS 与 CAD 软件的接口问题研究. 机械设计与制造, 2007(7): [12] 张晓东, 于擘, 王袖和. Ansys workbench 参数传递
174‒175辅助建立连续颌骨质量变化的种植体骨块三维有限元模型. 中国厂矿医学, 2008, 21(2): [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 Program, 2006 [14] 姜德义. 边坡稳定性分析与滑坡防治. 重庆: 重庆大学出版社, 2005 [15] Yelisetti S, Spence G D, Riedel M. Role of gas hydrates in slope failure on frontal ridge of northern
441‒458 Cascadia margin. Geophysical Journal International, 2014, 199(1):