岩石巴西圆盘动态破坏过程数值模拟邱流潮（中国农业大学应用力学系, 北京 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 uM -=+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 ss 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 verticalstresses along vertical diameter4 结 语本文应用有限元法与离散元法的耦合模型对脆性岩石的劈裂破坏过程进行数值模拟，分析结果验证了这种数值方法分析脆性岩石变形与破坏过程的可行性。
需要说明的是，本文只是这方面的一种尝试,在很多方面（包括断裂模型以及网格重划分）都还需要改进，随着研究的深入,将不断趋于完善。
参考文献：[1] W.C. Zhu, C.A. Tang. Numerical simulation ofBrazilian disk rock failure under static and dynamic loading [J]. International Journal of Rock Mechanics & Mining Sciences 2006, 43: 236–252. [2] P. A. Cundall, A computer model for simulatingprogressive large scale movements in blocky rock systems[M]．Proceedings of the International Symposium Rock Fracture，ISRM ，Nancy，vo1．1，1971.[3] 王泳嘉,邢纪波. 离散单元法及其在岩土力学中的应用[M]. 沈阳:东北工学院出版社,1991. (Wang YJ ,Xing JB. The Discrete Element Method and It s Applicationin Soil and Rock Mechanics [M]. Shenyang: Northeast Engineering College ,1991. (in Chinese))[4] 张楚汉, 金峰, 周元德. 岩石和混凝土离散-接触-断裂分析[M].清华大学出版社,2008. (Zhang CH, Jin F, Zhou YD. Discrete-Contact-Fracture Analysis in Soil and Rock Mechanics [M]. Beijing: Tsinghua University, 2008.(in Chinese))[5] 李世海, 高波, 燕琳. 三峡永久闸高边坡开挖三维离散元数值模拟[J]. 岩土力学, 2002, 23(3): 272-277.(3-D simulation of the excavation of high steep slope of Three-Gorges permanent lock by distinct element method. Rock and soil mechanics 2002, 23(3): 272-277. (in Chinese))[6] C. A. Tang, Micromechanical Model for Simulatingthe Fracture Process of Rock, Rock Mech. Rock Engng.37 (2004), 25–56.[7] A. MUNJIZA, D.R.J.OWEN, N. BICANIC. A combinedfinite-discrete element method in transient dynamics of fracturing solids. Engineering Computations 1995, 12:145-174.[8] 唐志平. 激光辐照下充压柱壳失效的三维离散元模拟[J ]. 爆炸与冲击,2001 ,21 (l) :127.(Z. P. Tang. 3D Discrete Meco-Element Simulation of Failure Processes for Cylindric Shells Subjected to Inner Pressure and Surface Laser Irradiation [J ] .Explosion and Shock Waves ,2001 ,21 (l) :127. (in Chinese))[9] 胥建龙,唐志平. 离散元与有限元结合的多尺度方法及其应用[J ]. 计算物理,2003 ,20 (6) :477. (Xu J L ,Tang Z P. Combined Discrete/ Finite Element Multi2Scale Numerical Method and It s Application [J].Chinese Journal of Computational Physics ,2003 ,20 (6) :477. (in Chinese))[10] 邱流潮, 基于联合有限-离散元法的混凝土重力坝地震破坏过程仿真[J]，水力发电，2009, 35(5): 36-38.(Qiu Liu-chao, Seismic Rupture Process Simulation of Concrete Gravity Dam based on the Combined Finite-discrete Element Method [J]. Water Power, 2009, 35(5): 36-38.)[11] A. Munjiza, K. R. F. Andrews And J. K. White,Combined Single And Smeared Crack Model In Combined Finite-Discrete Element Analysis [J]. Int. J. Numer.Meth. Engng. 44 (1999) 41-57.[12] W.F. Chen, Plasticity in Reinforced Concrete [M],McGrawHill, New York, 1982.Numerical modeling of Failure of Brittle Rock using Combined Finite-DiscreteElements MethodQIU Liu-chao(Department of Applied Mechanics,China Agricultural University, Beijing 100083, China)Abstract: The theory of combined Finite-Discrete Elements Method is briefly summarized and used to simulate the rock deformation and failure process. Simulations of both the intact rock behaviour and the development of fractures in Brazilian tests are performed. The simulation results identify the failure mechanisms and failure patterns and validate the nature of tensile crack initiation and propagation in brittle rock. The numerical results also show that the combined Finite-Discrete Elements Method is suitable for modeling the rock deformation and failure process.Key words: rock failure; Brazilian tests; numerical simulation; finite element method; discrete element method.。