基于TCK损伤本构的岩石爆破效应数值模拟

Numerical simulation of effect of rock blasting based on TCK damage constitutive model
WANG Zhi-liang1, ZHENG Ming-xin2
(1. College of Engineering Science, University of Science and Technology of China, Hefei 230027, China; 2. Institute of Bridge & Road and Geotechnical Engineering, East China Jiaotong University, Nanchang 330013, China)
析软件，尽管其材料库比较丰富，但它还为用户自 定义本构提供了方便的接口功能，上述岩石损伤本 构的嵌入流程如图 1 所示。
ε v = ε v + dε v
Ls-dyna接口
ε v > 0?
Y
Cd = 5 ⎛ P ⎞ ⎛ K IC ⎞ −2 v max k⎜ ⎟ ε ⎟ ⎜ 2 ⎝ 3K ⎠ ⎝ ρ C ⎠
1 引 言
岩石在爆破载荷下相继损伤和破碎，正确地评 估岩石的损伤积累和准确地预测岩石破裂的范围一 直是工程爆破中两个重要任务，前人已对此做出不 少努力
[1-3]
展的。该模型认为，岩体预先存在着许多随机分布 的裂纹。在拉载的作用下，这些裂纹将被激活、成 核甚至贯通；在受压状态下材料服从理想弹塑性本 构，忽略了材料的硬化、软化等特性。 另一方面，半无限岩体介质中工程爆破是生产 实践中最常见的， 根据介质爆破理论和应力波知识， 高温、高压爆轰波在离炮孔稍远处演化为应力波， 当应力波遇到临空自由面时，其将反射为拉伸波， 拉伸波的峰值强度一旦超过介质的抗拉强度，就会 在介质内产生拉裂破坏区，这个拉裂破坏范围的合 理确定，对临空面附近建（构）筑物的保护以及整 个施工的安全都具有重要参考价值。
p 有效塑性应变 ε eff 和当前屈服极限 σ y 表示为 p ε eff =∫ t 0
⎛ 16 ⎞ ν = ν ⎜1 − Cd ⎟ 9 ⎝ ⎠
式中：ν 为岩体的初始泊松比。
（2）
2 p ij dt ε 3
（11） （12）
p σ y = σ 0 + β Epε eff
在某一压力水平 P 下，单位体积岩体的裂纹激 活数目的统计，Kipp 和 Grady[7]建议采用 Weibull 分布函数： ⎛ P ⎞ N =ξ⎜ ⎟ ⎝ 3K ⎠
第 29 卷第 1 期 2008 年 1 月
文章编号：1000－7598－(2008) 01－0230－05
岩 土 力 学 Rock and Soil Mechanics
Vol.29 No.1 Jan. 2008
基于 TCK 损伤本构的岩石爆破效应数值模拟
王志亮 1，郑明新 2
（1.中国科学技术大学 工程科学学院，合肥 230027；2.华东交通大学 道桥与岩土工程研究所，南昌 330013）
m
式中： β 为硬化参数； Ep 为塑性硬化模量，缺省时 p e ij ij − ε ij 取为弹性模量 E 的 0.1 倍[8]；而 ε =ε 。 式 （9） 中， 当 φ ≤ 0 时为弹性或中性加载； φ >0 为塑性硬化。 2.3 炸药的状态方程 对于炸药， 采用 JWL 状态方程[5]来模拟其压力
摘
要：脆性岩石在爆破载荷下的动态破碎行为是工程技术人员所关心的，采用数值法预估岩体破裂范围和损伤大小对施工
安全具有重大意义，其中合理的岩石损伤模型是关键。假设岩体破坏近似服从 Mises 屈服准则，且考虑到岩体的应变硬化特 性，把 Taylor-Chen-Kuszmaul（TCK）模型中的拉裂损伤演化方程和材料双线性、弹塑性本构耦合到一起，并简洁地嵌入到 有限元分析软件 LS-DYNA 中，通过半无限岩体中柱形药包爆破实例验证了其合理性与准确性。对球形药包爆破问题也进行 了数值分析。该损伤模型及其数值模拟在工程中具有一定的参考价值。 关 键 词：岩石；爆破效应；损伤模型；用户子程序；拉裂；数值模拟 文献标识码：A 中图分类号：TU 452
232
岩
土
力
学
2008 年
3.2
岩石本构嵌入流程
。由 实际计算域尺寸为 1 250 cm × 1 250 cm（图 2） 于炸药的半径相对来说很小，故而网格划分较密， 共有 62 500 个单元。 边界条件如下：左边界为对称边界，上表面为 临空自由边界，下边界和右边界均为非反射边界。
LS-DYNA 作为一个瞬时的、大变形有限元分
。为了描述岩石断裂和破碎的主要特征，
加强人们对岩石爆破物理过程的理解，实现智能化 爆破，当前在材料本构理论基础上改进或发展新型 的、合理的岩石爆破损伤模型仍具有重要意义。 TCK 模型是一个较早的岩石损伤模型，它由 Taylor等人在Kipp和Grady模型 （1980年） 基础上， 利用了含裂纹体的宏观等效弹性模量表达式 而发
m 2
N
各向同性硬化 的弹塑性本构
炸药
2.5 m 2.5 m 12.5 m 7.5 m 12.5 m
1 (ν ) = ∂f1 ∂ν C d f ∂ν ∂Cd
1 (ν )] / 9 d + Cd f = 16[ f1 (ν )C D
∆t D=D+D N
∆σ ij = K dδ ij ∆ε kk + 2Gd ∆eij
[4] [1] [2]
收稿日期：2006-05-11 作者简介：王志亮，男，1969 年生，博士，副教授，主要从事于岩土工程、工程力学等方面教研工作。E-mail: zliangw@ustc.edu.cn
第1期
王志亮等：基于 TCK 损伤本构的岩石爆破效应数值模拟
231
m 2
本文拟基于大型有限元分析软件 LS-DYNA 的 用户自定义材料本构的功能[5, 6]，把 TCK 模型中的 拉裂损伤演化规则和双线形弹塑性材料本构耦合起 来，并嵌入到该动力分析软件中，对半无限岩体中 爆破诱发损伤的发展规律及其分布范围进行深入探 讨，其方法、思路可为类似工程提供一定的参考。
2 φ = σ i2 − σ y
（9）
（1）
式中： σ i =
3 sij sij ；应力偏量 sij 表达式为 2 1 sij = σ ij − σ kk δ ij 3
（10）
式中： f1 (ν ) = (1 − ν 2 ) /(1 − 2ν ) ； K 为岩体的初始体 积模量； K ，ν 分别为损伤材料的体积模量和泊松 比； Cd 为岩体中的裂纹密度； D 为损伤标量。 为了简化计算，对于ν 采用如下的表达式[2]：
Abstract: Engineers and technicians are concerned with the dynamic fracture behaviors of brittle rock under blast loading. It is of very importance for construction safety to predict the fragmentation scope and damage value by use of numerical approach. However, a reasonable damage constitutive model is the key thereinto. By assuming rock observes the Mises yield criterion and considering the strain-hardening behavior of rock, the tension damage evolution equation of Taylor-Chen-Kuszmaul model is coupled with the bilinear elastoplastic constitutive model. Subsequently, the developed constitutive model is concisely implemented into the FEM software, LS-DYNA. The rationality and accuracy are tested via a practical example of cylindrical charge explosion in a semi-infinite rock mass. At last, one spherical explosion is numerically analyzed. This damage model and related numerical simulations can be of certain value in a practical rock blasting. Key words: rock; blasting effect; damage model; user-defined subroutine; tensile fracturing; numerical simulation
Cd = N (ω r 3 )
3 岩石损伤本构数值嵌入流程
3.1 数值分析工具 本文采用著名的显式动力分析有限元软件LS-
（5）
DYNA来完成计算的。该软件为爆炸、冲击乃至物
体跌落等提供了分析工具，后来经过多次的功能扩 充和改进，发展成为当今强有力的非线性动力有限 元软件。
式中： ω 为未知的比例系数。 把式（3） 、式（4）代入式（5） ，并令 k = ξω ， 则
K d = (1 − D) K Gd = (1 − D)G
Leabharlann Baidu
图 2 柱形药包爆破的计算域尺寸 Fig.2 Geometry of explosive column
D > 1.0 ?
Y
D = 1.0
图 3 给出了炸药起爆后 2 ms 的损伤分布云图。 可看出损伤由两部分组成：柱形炮孔周围是由较高 的径向、切向拉应力形成的破碎区域；而从药卷上 方到地表间的损伤区域是由应力波在地表自由面反 射后拉裂形成的，整个外形象一顶“帽子” ，而帽子 中、上部便构成了地表附近的“爆破漏斗”[10]。由 于有限元计算中要求介质是连续的，因此，爆破造 成的介质抛掷现象不能模拟。实际上，这个“爆破 漏斗”中部分岩石碎块要抛射出去的[11]。合理地确 定这个“爆破漏斗”范围具有重要意义，它关系到 该范围内的建筑物安全和土石方的开挖费用等。
K 16 =1− f1 (ν )Cd = 1 − D 9 K
依据广义 Hooke 定律，应力-应变关系表达为
∆σ ij = K dδ ij ∆ε kk + 2Gd ∆eij
（8）
式中：Kd 和 Gd 均为损伤材料的模量；∆ε kk 和 ∆eij 为 体应变和偏应变增量； δ ij 代表 Kronecker 记号。 2.2 双线形弹塑性本构 对于岩体介质处于压缩状态下时，设其行为服 从各向同性硬化的弹塑性本构[5]，其屈服函数为
Cd =
5 ⎛ P ⎞ ⎛ K IC ⎞ −2 v max k⎜ ⎟ ε ⎟ ⎜ 2 ⎝ 3K ⎠ ⎝ ρ C ⎠
（6）
由式（1）可得出损伤演化方程为
D=
16 f1 (ν )Cd 9
（7）
2 材料的本构关系
2.1 TCK 拉裂损伤演化方程 岩石中缺陷和裂纹是预先存在的，它们使得岩 体存在初始损伤。 这些裂纹在外载下会进一步发展。 Budiansky 和O’Connell[4]曾利用“自洽法”确定出 布满裂纹固体的弹性模量。TCK模型认为，岩体内 随机分布的裂纹在拉载下将被激活并相互作用，对 于存在预裂纹的各向同性岩体介质，其有效体积模 量为
式中：A，B， R1 ， R2 和 ω 均为材料参数；E 为比 内能，其初始值为 E0 ；相对体积 V 的初始值为 V0 。
1 ⎛ 20 K IC r= ⎜ v max 2⎜ ⎝ ρ Cε
⎞ ⎟ ⎟ ⎠
2/3
（4）
式中： ρ 为材料密度；C 为纵波波速； K IC 为材料 v max 为断裂时的体应变率。 的断裂韧性； ε 而裂纹密度 Cd 和平均碎块半径 r 间关系为
（3）
式中： ξ ，m 分别为 Weibull 分布参数。 针对脆性材料在动态断裂中的平均碎片块度 半径，Grady 曾提出了一个[2]计算式：
p 和相对体积 V 之间关系：
ω ⎞ − R1V ω ⎞ − R2V ω E ⎛ ⎛ （13） + B ⎜1 − + p = A ⎜1 − ⎟e ⎟e V ⎝ R1V ⎠ ⎝ R2V ⎠