基于流固耦合加载技术的DDA方法.pdf

第 33 卷 第 1 期 2016 年 3 月 爆 破 BLASTING Vol 33 No 1 Mar 2016 doi 10 3963/ j issn 1001 -487X 2016 01 002 基于流固耦合加载技术的 DDA 方法* 余德运 1， 刘殿书1， 王庆斌2， 何成龙3 （1 中国矿业大学 （北京）力学与建筑工程学院， 北京 100083； 2 保利民爆哈密有限公司， 哈密 839200； 3 北京理工大学 爆炸科学与技术国家重点实验室， 北京 100081） 摘 要 针对目前 DDA 方法尚没有实现对爆炸荷载耦合加载， 而常用的两种爆炸荷载等效加载方法又存 在不足的现状， 在原 DDA 程序中引入基于空间网格覆盖判断的耦合加载子程序， 并用来模拟集中药包爆破。 结果表明 基于网格覆盖判断的耦合技术的 DDA 方法， 能体现岩石爆破中爆生气体和爆炸冲击波共同作用； 能实现爆生产物气体对周围岩石介质的准确、 持续加载； 能体现最小抵抗线对爆生气体作用时间的影响； 能 够较好地对模拟岩石爆破过程中， 爆炸产物气体对孔壁的作用， 产物气体在裂缝中泄漏时对裂缝的扩张作 用， 以及抛掷过程中产物气体对岩块的驱动作用。 关键词 非连续变形分析；流固耦合；数值模拟 中图分类号 O342 文献标识码 A 文章编号 1001 -487X （2016） 01 -0006 -06 An Improved DDA based on Fluid-Structure Coupling Loading Technology YU De-yun1， LIU Dian-shu1， WANG Qing-bin2， HE Cheng-long3 （1 School of Mechanics Civil Engineering， China University of Mining and Technology（Beijing） ， Beijing 100083， China； 2 Poly Civil explosive Hami Co Ltd， Hami 839200， China； 3 State Key Laboratory of Explosion Science and Technology， Beijing Institute of Technology， Beijing 100081， China） Abstract The present DDA could not per the blast load by fluid-structure coupling， and there two e- quivalent s have deficiency The coupling subroutine based on spatial grid coverage judgment was added to the original DDA program， and used to simulate the charge blasting The result showed that the extended DDA program could achieve accurate and sustained blasting load， it could reflect the effect of minimum burdens on the detonation gas loading time and blasting scope， it could simulate the blowing， the cavity expanding， the explosive gassed escaping process Key words DDA；fluid-structure interaction；numerical simulation 收稿日期 2015 -10 -14 作者简介 余德运 （1981 - ） ， 男， 博士后， 湖北武汉人， 主要从事爆破 安全技术及数值仿真计算研究，（E-mail） yudeyun2000 163 com。 基金项目 国家自然科学基金项目 （No 51374038） 重复爆破荷载作 用下岩石动态疲劳损伤特性研究 DDA （Discontinuous Deation Analysis） 方法 是 1984 年由石根华提出的一种对诸如岩石等不连 续块体系统的力学响应问题进行数值计算的方 法 ［1］。DDA （非连续变形分析） 方法自其问世以来 已经在岩土工程中得到了广泛的研究和应用。岩石 爆破过程中， 爆炸产物气体对孔壁的作用涉及流固 耦合问题， 产物气体在裂缝中泄漏时对裂缝的扩张 作用和抛掷过程中产物气体对岩块的驱动作用等都 属于流固耦合问题。目前的 DDA 方法， 尚没有实现 对爆炸荷载的耦合加载， 现在主要采用两种等效的 方式实现爆炸荷载的加载 第一种方法是直接利用 原经典 DDA 程序中的点载荷加载方式， 将炮孔壁的 爆炸压力等效为半理论半经验的爆炸荷载压力曲线 （如三角型脉冲荷载） ， 直接以点载荷的形式加 载 ［2］； 第二种方法是甯尤军提出的根据爆生气体理 想状态方程计算爆腔即时压力进行加载 ［3］， 该方法 通过跟踪炮孔的扩张和炮孔周边裂隙的发展贯通、 将爆生产物气体视为理想气体， 采用理想气体状态 方程计算爆腔的即时压力， 并将这一压力作用到主 炮孔内壁和贯通裂隙面上。 第一种爆炸载荷加载方法中， 采用半理论半经 验的爆炸荷载压力曲线， 点载荷的坐标位置、 加载点 数量和加载时间在建模时确定， 在块体运动过程中 点载荷与块体相对位置保持不变、 点载荷数量保持 不变。爆炸载荷作用时间和压力大小由载荷曲线确 定， 与爆腔发展不相关。该方法不能体现爆生气体 产物泄漏对压力的影响， 不能体现爆生产物气体对 裂缝的扩张作用。 第二种爆炸载荷加载方法中， 不考虑爆轰产物 气体由瞬时爆轰压力膨胀到某一临界压力这一过 程， 即不考虑炸药爆炸这一冲击载荷的升压过程， 只 考虑爆生气体从某一临界压力膨胀至大气压这一卸 载过程。该方法能够体现爆生产物气体产物泄漏对 压力的影响， 能够体现产物气体对裂缝的扩张作用， 也能够对 DDA 块体进行持续加载。但是， 该方法在 计算过程中， 将爆生产物气体视为理想气体， 并且还 认为是完全气体 （比热比 γ 为常数） ， 爆腔内各点的 压力都相等， 因此， 该方法得到的爆腔压力是近似 值， 难以实现爆生产物气体对周围介质的准确加载。 为了准确地反映岩石爆破过程中， 炸药爆炸产 生的冲击波及爆轰产物气体与周围结构的相互作 用， 得到比较准确的爆炸流场发展过程和较为准确 的爆炸加载， 体现爆炸产物气体对岩石破裂的作用， 以及实现对块体进行持续加载， 采用一种基于重叠 网格空间覆盖判断的流固耦合计算技术 ［4， 5］， 在原 DDA 程序中， 添加采用流体动力学计算炸药爆炸压 力和基于空间覆盖判断处理耦合作用的耦合加载子 程序， 实现岩石爆破中爆生产物气体对周围岩石介 质的准确、 持续加载。 1 基于网格覆盖的流固耦合技术的 DDA 方法 1 1 计算流程 在对爆炸流场与结构相互作用这个流固耦合问 题进行数值模拟时， 可以通过 Eulerian 流体弹塑性 模型实现对流场和固体结构的统一描述， 流体和固 体之间的相互作用体现为同一计算区域中不同物质 界面间的识别及处理 ［6］。其中， 岩石变形及运动过 程由 DDA 方法进行计算， 炮孔内炸药爆炸产生的气 体压力由多物质流体动力学方法进行计算， 流场网 格与 DDA 块体在空间上可以自由重叠。程序计算 过程中， 以流体程序为主线进行求解， 在流体程序的 时间步间进行耦合处理和 DDA 求解。两种程序按 各自的时间步长计算规则计算时间步长， 选择两种 程序的时间步中较小的时间步作为整个程序的时间 步长。程序总的计算流程如图 1 所示。 图 1 总的计算流程 Fig 1 The whole calculation process 1 2 流体与 DDA 块体空间覆盖判断 在对 DDA 单元划分和对流体系进行网格划分 时， 为了提高求解效率， 对单元尺寸提出要求 （a） DDA 块体单元根据实际需要， 正常划分块 体单元。 （b） 流体网格为规则的矩形网格， 且流体 单元网格最大边长小于 DDA 块体单元最小边长。 图 2 为流体网格与 DDA 块体在空间上出现重 叠， 即两者的网格出现了空间覆盖情况， 对两者空间 相对位置进行判断。对于流场单元， 根据其被 DDA 块体单元重叠情况， 可把流场区域分为 4 类 未被覆 盖区、 全部被覆盖区、 部分被覆盖区域和边界与 DDA 边界重合区域， 在图 2 中， 分别用 “Ⅰ、 Ⅱ、 Ⅲ、 Ⅳ” 标记。 对于图 2 中的四种类型的流场区域， 按以下原 则进行处理。 （1） 对于未被 DDA 块体覆盖的流场区域Ⅰ， 流 体单元的四个节点都位于结构区域之外， 且没有结 构单元的节点位于这个流体单元内， 即这个流体单 元没有被 DDA 块体覆盖。此类流体单元， 在流场计 算中不需要进行任何特殊处理。 7第 33 卷 第 1 期 余德运， 刘殿书， 王庆斌， 等 基于流固耦合加载技术的 DDA 方法 图 2 流体网格与 DDA 块体空间覆盖 Fig 2 The fluid grid covered with DDA block space （2） 对于全部被 DDA 块体覆盖的流场区域Ⅱ， 流体单元的四个节点都位于结构区域内， 即这个流 体单元完全被 DDA 块体覆盖。此类流体单元， 在流 场计算中将这些网格设置成虚网格， 在流体程序的 时间循环中不需要进行计算， 通向此类网格的物质 输运也被截断。 （3）对于部分被 DDA 块体覆盖的流场区域Ⅲ， 流场单元有 1 至 3 个节点位于结构区域内， 即这个 流体单元部分被 DDA 块体覆盖。流固耦合相互作 用就是通过对这类网格的处理得以实现的。对这类 网格的处理包含两部分 对网格量的修正和对输运 的修正。 （4）对于边界与 DDA 块体边界重合的流场区 域Ⅳ， 流体单元的边与 DDA 块体的边重合， 其又可 分为两种情况， 如图 3 所示。结构单元边界与流体 单元的两条边相交和结构单元边界与流体单元的三 条边相交。对于此类流体单元， 在流场计算中也需 对网格量和输运进行修正。 图 3 流体网格与 DDA 块体边界重合情况 Fig 3 The fluid grid coincide with DDA block boundary 1 3 流体网格量和输运的修正 根据流固耦合作用物理过程特点， DDA 块体在 流体单元体积应力载荷下产生运动， 反过来流体单 元也受到 DDA 块体的反作用。DDA 块体对流体的 作用体现在两个方面 对网格中流体物质的压缩 （反之则是稀疏， 在这里作为一致过程进行处理） 和 对流体物质在网格间输运的阻挡， 其中 DDA 块体对 流体单元的压缩作用表现在导致网格内物质体积变 化和物质速度的变化， DDA 块体对流体物质在网格 间输运的阻挡表现在输运量的变化， 所以必需对流 体网格量和输运进行修正。 1 3 1 流体网格量的修正 DDA 块体对流体网格中物质的压缩导致网格 中流体物质体积和速度这两个方面发生变化， 需要 分别进行处理。 如图 4 所示， 在两个时间步之间， 流体网格中的 DDA 块体边界发生了运动， 流体网格的未被覆盖部 分体积变小， 网格内物质被压缩， 导致密度和压力增 大。通过几何判断， 可以得到流体网格被 DDA 块体 覆盖的面积， 将剩余面积与初始面积的比率作为修 正系数， 以此来对网格中物质密度进行修正， 然后使 用状态方程， 重新计算流体网格压力。 图 4 流体网格的压缩 Fig 4 The compression of fluid grid 对于被 DDA 块体部分覆盖的流体网格， 其被覆 盖面积计算， 采用的是将任意多边形分解为多个三 角形 （如图 5 所示） ， 各三角形的面积利用单纯形积 分方法计算 ［7］， 然后求所有三角形面积之和， 即为 流体单元被 DDA 块体覆盖的多边形的面积。 1 3 2 流体网格输运的修正 在流体动力学计算程序中， 一般由应力效应步 计算网格速度得到输运速度。在本文的多物质流体 动力学程序中， 采用贡献-接受网格法进行网格间的 物质输运 ［8］， 流体物质通过网格边界输运到相邻网 格中， 但对于被 DDA 块体覆盖的流体网格， 由于 8爆 破 2016 年 3 月 DDA 块体的进入， 流体单元与 DDA 块体在空间上 重叠， 某些流体网格边界被 DDA 块体全部或部分覆 盖， 如图 6 （a） 中的流体网格边 1 和边 2 被部分覆 盖， 图 6 （b） 中的流体网格边 1 和边 3 被部分覆盖边 2 被全部覆盖， 需要对通过这几条边的物质输运进 行输运体积修正 首先根据 DDA 块体侵入网格过程 中的动量守恒来求得这些流体网格速度， 速度的方 向由 DDA 块体边界段的法向确定； 再根据速度的变 化得到流体网格能量的变化。 图 5 流体网格覆盖面积计算 Fig 5 The fluid grid coverage calculation 图 6 DDA 块体对流场中物质的压缩和输运阻挡 Fig 6 The compression and transportation to the material in fluid grid 程序计算时， 先对被部分覆盖的流体网格边进 行标记， 求网格边的覆盖率， 然后修正这条边上的输 运体积。 1 4 流场对 DDA 块体压力加载 流场对 DDA 块体的作用， 是通过在 DDA 块体 边界单元的外边上施加压力的方式实现。因此， 只 需要处理单元边为固体边界边情况的单元即可。某 个 DDA 块体边界有可能位于一个流体网格， 也有可 能属于多个流体网格。对于 DDA 块体边界边位于 一个流体网格的情况， 使用这个流体网格中的压力 作为这个单元边的加载条件； 对于边界边位于多个 流体网格的情况， 如图 7 所示。对于此种情况， 需要 根据此 DDA 块体单元的边被不同流体网格所截片 段的长度， 对这个几流体网格中的压力通过加权平 均计算这条边上的压力载荷 P L1P1 L2P2 L3P3 L （1） 式中 L1、 L2、 L3分别为此 DDA 块体的某边在三个流 体网格中的长度， 且有 L1 L2 L3 L； P1、 P2、 P3分 别为三个流体网格的压力。 计算得到了这个 DDA 块体单元边上的压力后， 把这个压力值作为边界压力载荷加载作用于 DDA 的 边上。流体对 DDA 块体边界的载荷形式为线载荷。 图 7 流场对 DDA 块体边界的压力加载 Fig 7 The pressure loading on DDA block boundary 2 数值算例 在现实条件下， 炸药在岩石中爆炸时， 爆炸产生 的压力高达 104MPa 量级， 而最坚硬的岩石的抗压 强度也只有 102MPa 量级， 所以在高温高压的爆炸 产物作用下， 邻近装药的岩石受到强烈压缩， 结构完 全破坏， 颗料被压碎， 形成一个 “粉碎区” ［9］。而在 DDA 程序中， DDA 单元本身不会消失， 本文中的 DDA 程序还没有实现 DDA 子块体的内部开裂， 为 了合理解决这一矛盾， 在建模时， 采用于不耦合装药 方法， 以降低直接作用在 DDA 块体上的爆生产物气 体压力。 2 1 模型 为了提高计算效率， 在确定流场区域时， 只需将 可能发生流固耦合作用的范围建立流场区域即可。 当爆腔扩大至自由面， 即腔内爆生产体产物逸出至 空气中时， 爆生气体产物体积迅速扩大， 其压力也迅 速降至大气压。因此， 流场计算范围以爆腔中心到 自由面的距离来确定。 图 8 为集中药包在岩体中爆破的模型。药包半 径为 0 30 m， 岩体平面尺寸为 长 8 5 m， 厚度为 5 0 m， 药包埋深分别为 1 0 m、 1 5 m 和 2 0 m， 即 集中药包爆破的最小抵抗线分别为 W 1 0 m、 W 1 5 m 和 W 2 0 m。岩石材料及节理参数见 表 1， 炸药参数见表 2。 2 2 模拟结果分析 图 9 图11 分别是最小抵抗线为 1 0 m、 1 5 m 和 2 0 m 时集中药包爆破的抛掷过程。 9第 33 卷 第 1 期 余德运， 刘殿书， 王庆斌， 等 基于流固耦合加载技术的 DDA 方法 图 8 爆破漏斗计算模型 Fig 8 A model for calculating the blasting funnel 表 1 爆破岩石块体材料及节理参数 Table 1 The parameters of block and joint 岩石材料参数 节理参数 密度 ρ/ （kgm -3） 弹性模量 E/ GPa 泊松比 ν 抗拉强度 σft/ MPa 抗拉强度 σft/ MPa 摩擦角 φ/ 粘结力 c/ MPa 抗拉强度 σt/ MPa 2650450 24603 525205 0 表 2 炸药材料参数 Table 2 The explosive material parameters 材料 密度 ρ/ （kgm -3） 爆速 v/ （ms -1） C-j 压力/ GPa A/ GPaB/ GPa 材料常数 C/ GPaR1R2ω 弹性模量/ GPa 2岩石炸药120035001 84214 40 1820 5164 20 90 15492 图 9 最小抵抗线 W 1 0 m 的抛掷过程 Fig 9 Casting process with the minimum burden W 1 0 m 图 10 最小抵抗线 W 1 5 m 的抛掷过程 Fig 10 Casting process with the minimum burden W 1 5 m 图 11 最小抵抗线 W 2 0 m 的抛掷过程 Fig 11 Casting process with the minimum burden W 2 0 mm 01爆 破 2016 年 3 月 从图 9、 图 10 和图 11 中可以看出， 药包上方块 体受爆生气体膨胀和爆炸冲击波的共同作用， 爆生 气体的作用体现在对爆腔的膨胀鼓包； 爆炸冲击波 的作用可从图9 （b） 、 10 （b） 和11 （b） 中可以看到 爆 炸冲击波传播至自由面时发生反射， 压缩波反射产 生拉伸波， 使自由面附近块体被抛起。因此， 爆炸耦 合加载方式的 DDA 程序， 能够体现爆生气体和冲击 波的共同作用。 图 12 为抵抗线分别为 W 1 0 m、 W 1 5 m 和 W 2 0 m 时， 集中药包在岩体中爆炸时的爆腔 中心压力-时间曲线。 图 12 不同抵抗线条件下集中药包爆腔中心压力-时间曲线 Fig 12 The pressure-time curve of concentrated charge blasting under different W 从图 12 可以看出， 爆腔中心的压力可明显分为 两个阶段 快速的加载阶段和动态卸载阶段。对于 加载阶段， 当最小抵抗线分别为 1 0 m、 1 5 cm 和 2 0 m 时， 压力很快上升到最大值， 其加载阶段压力 曲线基本重合， 这说明加载过程基本不受最小抵抗 线的影响， 这是因为在此阶段， 主要是炸药爆炸过 程， 爆炸压力峰值及爆轰过程由炸药本身的特性决 定， 所以在此阶段， 加载曲线重合； 对于动态卸载阶 段， 爆腔中心的压力衰减速度、 压力大小及作用时间 会随着最小抵抗线的不同而不同， 爆腔压力衰减速 度随最小抵抗线的增大而减小， 爆生气体压力作用 时间随最小抵抗线的增大而增大， 即最小抵抗线的 增大， 有利于提高爆生气体的能量利用率。因此， 爆 炸耦合加载方式的 DDA 程序， 能够得到较准确且持 续的爆炸载荷， 且能够体现最小抵抗线对爆生产物 气体作用时间的影响。 3 结论 （1） 基于网格覆盖的流固耦合技术的 DDA 方 法， 与直接加载半理论半经验的爆炸荷载压力曲线 相比， 能够实现爆生气体对周围介质的持续加载； 与 根据爆生气体理想状态方程计算爆腔即时压力进行 加载相比， 充分考虑了炸药爆轰至爆生产物气体膨 胀这一全过程， 可以得到更准确的爆炸载荷。 （2） 基于网格覆盖的流固耦合技术的 DDA 方 法， 能够体现爆生气体和冲击波的共同作用， 能体现 最小抵抗线对爆生气体作用时间及爆破效果的影 响， 能够较好地模拟岩石爆破过程中， 爆炸产物气体 对孔壁的作用， 产物气体在裂缝中泄漏时对裂缝的 扩张作用， 抛掷过程中产物气体对岩块的驱动作用。 参考文献 （References） ［1］ SHI G H， GOODMAN R E Discontinuous deation a- nalysis ［C］ ∥Proceedings of the 25th US Symposium of Rock Mechanics Evanston， USA， 1984 269-277 ［2］ A Mortazavi， P D Katsabanis Modelling burden size and strata dip effects on the surface blasting process ［j］ In- ternational journal of Rock Mechanics and Mining Sci- ences， 2001（38） 481-498 ［3］ 甯尤军， 杨 军， 陈鹏万 节理岩体爆破的 DDA 方法模 拟 ［j］ 岩土力学， 2010， 31 （7） 2259-2263 ［3］ NING You-jun， YANG jun， CHEN Peng-wan Numerical simulation of rock blasting in jointed rock mass by DDA ［ j］ Rock and Soil Mechanics， 2012， 31 （7） 2256-2263 （in Chinese） ［4］ 浦锡锋， 王仲琦， 白春华， 等 多物质流体动力学方法 与结构动力学方法结合的流固耦合计算技术 ［j］ 计 算物理， 2010， 27 （6） 833-839 ［4］ PU Xi-fen， WANG Zhong-qi， BAI Chun-hua， et al Cou- pling eulerian multi-material hydrodynamic and lagrangian structure dynamic for numerical simulation of flu- id-structure interaction［j］ Chinese journal of Computa- tional Physics， 2010， 27 （6） 833-839 （in Chinese） ［5］ 浦锡锋， 王仲琦， 余德运 多物质流体动力学方法与 DDA 方法结合的流固耦合计算技术 ［j］ 兵工学报， 2012， 33 （S1） 221-226 ［5］ PU Xi-fen， WANG Zhong-qi， YU De-yun Coupling multi- material hydrodynamic and dda mtehod for numerical sim- ulation of explosion in rock［ j］ Acta Armamentarii， 2012， 33 （S1） 221-226 （in Chinese） ［6］ 恽寿榕， 涂侯杰， 梁德寿， 等 爆炸力学计算方法 ［M］ 北京 北京理工大学出版社， 1995 ［7］ 石根华 数值流行方法与非连续变形分析 ［M］ 裴觉 民， 译 北京 清华大学出版社， 1997 ［8］ 顾尔祚 流体力学有限差分法基础 ［M］ 上海 上海交 通大学出版社， 1988 ［9］ 陈宝心， 杨勤荣 爆破动力学基础 ［M］ 武汉 湖北科 学技术出版社， 2005 11第 33 卷 第 1 期 余德运， 刘殿书， 王庆斌， 等 基于流固耦合加载技术的 DDA 方法