岩石巴西圆盘动态破坏过程数值模拟_邱流潮

岩石巴西圆盘动态破坏过程数值模拟

邱流潮

（中国农业大学应用力学系, 北京 100083, email:qiuliuchao@ ）

摘 要：岩石变形与破坏过程的研究对于岩石破碎等与工程密切相关的问题具有重要的作用。本文采用有限元与离散元耦合数值方法模拟岩石的变形和破坏过程。通过对巴西圆盘实验的数值计算，再现了脆性岩石的破坏机理和破坏模式以及岩石裂缝萌生、裂缝扩展以及完全破坏的全过程。

关键词：岩石破坏；巴西圆盘实验；数值模拟；有限元法；离散元法 中图分类号：TU452 文献标识码：A 文章编号：

0 引 言

岩石变形与破坏过程的研究对于岩石破碎等与工程密切相关的问题具有重要的作用。与其它研究方法相比较而言，数值方法可以考虑岩石的非均匀特性和所处的复杂受力条件，具有较强的适用性[1]

。目前，有限元和离散元法已被广泛应用于工程

问题的数值分析中。Cundall [2]，王泳嘉与邢纪波[3]

，

张楚汉与金峰[4],李世海[5]

等成功应用离散元技术

分析岩石工程问题。唐春安[6]

采用有限元法并结合细观材料本构模型分析岩石渐进破坏过程。一般而言，有限元方法适合于解决连续介质问题,而离散单元法适合于非连续介质问题。因此,将离散元法与有限单元法有机地结合起来,可以充分发挥各自的长处,特别是在模拟细观多相材料动态特征以及涉及从连续介质到非连续介质转化的研究方面（譬如岩石的渐进破坏过程等）,其优势尤为突出。近

年来，Munjiza 与Owen [7]

利用离散元与有限元结合的方法分析多重断裂固体和离散系统的模拟。唐志

平等人[8-9]

用离散元有限元相结合计算受拉铝板在激光辐照下的破坏模拟以及多尺度结合应用研究。

邱流潮[10]

应用有限-离散元耦合法对混凝土重力坝地震破坏过程进行了仿真。本文简要介绍了有限元与离散元耦合方法的基本原理，并用该方法研究岩石试样的劈裂破坏过程。数值模拟再现了岩石从起裂，裂纹扩展，到完全破碎的全过程。

1 运动方

可变形体动力平衡方程的有限元离散形式为：

)(int u f f u C u

M -=+ext （1）

式中，M ,C 分别表示质量矩阵和阻尼矩阵。ext f 表

示包括接触力在内的所有外力向量，)(int u f 表示变

形引起的内力向量。u

，u ，u 分别表示加速度，速度和位移向量。对于显式时间积分方法，质量矩阵M 需要对角化处理。如果粘性阻尼采用质量比例

基金项目：中央高校基本科研业务费专项资金(2009140)和国家自然科学基金（10872214）资助

阻尼，则阻尼矩阵C 定义为C =αM , α常数比例因子。当(1)式求解准静力问题时,则采用动力松弛法求解，选取较高的阻尼值（譬如，常数比例因子α取0.3）可使计算在足够的时间后收敛到静力解。

刚体运动由平动和转动两部分组成，刚体i 的运动方程为：

i i i i i i c m G F v a -=+ （2）

i i dt d M L =/

（3）

其中，

i a ，i v 分别为刚体i 形心的加速度和速度；i m 为刚体i 的质量；阻尼i c =αi m ，α为质量阻尼系数；i F 为刚体受到的除重力以外的所有力的合力，包括接触力和外荷载；i G 为刚体的重力；i L 为刚体的角动量，i M 为刚体受到的合力矩。

2 断裂模型

断裂模型是有限元与离散元耦合方法中的重要

组成部分。目前，有限元和离散元法中采用的断裂模型有很多中，包括全局方法和局部方法。而局部方法包括虚裂缝法（基于软化或损伤模型）以及单

裂缝法。在本文的计算中，采用Munjiza等人[11]

建议的虚裂缝和单裂缝的组合模型。其基本思路是：在到达材料极限抗拉强度之前，采用常规的有限元技术结合材料本够关系进行计算。而进入软化阶段后，则应用单裂缝模型进行模拟。为了简单起见，假设裂缝沿单元的边展开。当某个单元边上的黏结应力达到材料的抗拉强度时，裂缝便沿该单元边萌生。而当裂缝开度超过临界缝宽时，裂缝扩展，此时应力减少到零。这种方法的主要好处是：它能模拟多重裂缝的萌生和扩展而不需要额外增加太多的计算时间。

3 数值计算

本文以巴西圆盘试样作为研究对象，该试样的几何尺寸以及有限元网格如图1所示。岩石试样的直径为150 mm，在试样的上、下部分别有一个加载板，同时，在每个加载板与试样之间有一个垫层，

数值计算中，加载板设为刚体。垫层为均匀弹性材料，弹性模量为60GPa，泊松比均为0.2。该岩石试样的弹性模量为30GPa，泊松比均为0.2，抗拉

强度为3.2 MPa，密度为2700 kg/m 3

，断裂能为100 N/m。加载板以垂直速度0.003 m/s 运动（上板向下，下板向上）来模拟加载过程。计算模型采用常应变三角形单元离散（假设为平面应力状态），整个模型包含1114个单元和3342个节点。

图2给出了巴西圆盘试样中心点处水平应力的响应时程。从图中可以看出，在达到峰值载荷之前应力－时间曲线基本上呈线性趋势。而超过极限载荷后，试样承载能力急剧下降。

h o r i z o n t a l s t r e s

s e (M P a )

time (ms)

图 2. 试样中心点处水平应力的响应时程

Fig. 2. Stress-time curve at the central of rock disk

图3给出了几个典型时刻（对应图2中应力－时间曲线上的圆点）的开裂模式和水平应力分布。从图中可以看出，大约在 t = 38 ms 时裂缝开始在试样中心点处萌生，因为此处的水平应力达到了最大拉应力。随后裂缝沿垂直径向扩展，在大约在 t = 55 ms 时整个圆盘完全断开。

contour at selected times

图4给出了几个典型时刻（对应图2中应力－时间曲线上的圆点）沿圆盘式样垂直径向的水平应力和垂直应力分布。从图中可以看出, 圆盘式样中应力分布为拉－压类型。在达到峰值加载之前，圆盘式样垂直径向大约四分之三的范围的水平应力均为拉应力，而沿整个垂直径向的垂直应力都是压应力。紧靠加载垫片下的区域，水平应力和垂直应力幅值增长剧烈，且这些区域基本处于三向压力状态。以上这些数值预测结论和文献[12]中的线弹性理论分析结果是吻合的。随着加载的进行，大约在 t = 38 ms 时裂缝开始在试样中心点处萌生，随后裂缝沿垂直径向扩展，在大约在 t = 55 ms 时沿整个垂直径向的水平应力和垂直应力都减小为零，说明此时整个圆盘沿垂直径已完全断开。

y (m m )

stress (MPa)

y (m m )

stress (MPa)

y (m m )

stress (MPa)

y (m m )

stress (MPa)

图 4. 圆盘垂直径向的水平应力和垂直应力分布 Fig. 4. Distribution of horizontal and vertical

stresses along vertical diameter