0.   引言

在大型建筑构件生产过程中，由于施工工艺等原因，会产生大量初始缺陷。当荷载超过承载力范围时，这些初始缺陷可能导致构件整体失效，从而引发工程事故。随着计算机技术的发展，数值试验成为一种重要的研究手段，因此可通过数值方法对缺陷的裂纹扩展规律进行深入研究。

有限元法是目前广泛应用于固体力学和结构分析等的连续型数值方法。对于动态破裂问题，包括裂纹萌生、裂纹扩展、裂纹分叉、复杂裂纹相互作用等非连续问题，目前已有长足发展。Belytschko等^[1]基于扩展有限元(the extended finite element method，XFEM)方法提出了一种能够处理裂纹尖端位于单元内部的方法，并用于动态裂纹扩展的模拟。另一种处理非连续问题的方法是在有限单元边上设置黏结面，裂纹则只能沿着单元边界产生。Xu等^[2]采用这一方法模拟了裂纹扩展问题，并分析了裂纹扩展速度过快的原因。另外，其他连续方法在模拟非连续问题上得到了广泛发展，例如岩石破裂过程分析系统(rock failure process analysis system, RFPA)通过单元弱化，实现了模拟岩石从弹性阶段到完全破坏全过程。梁正召等^[3]采用RFPA模拟了岩石三维裂纹扩展规律。马鹏飞等^[4]采用近场动力学方法进行了裂纹扩展模拟。

与连续方法相比，离散型数值方法优势在于能够处理非连续问题。目前，离散单元法(discrete element method, DEM)^[5]已经广泛应用于研究脆性材料的破坏。其中，PFC(particle flow code)^[6]属于采用颗粒离散的一种DEM。Kou等^[7]基于PFC的黏结颗粒模型研究了不同工况下裂纹分叉问题，并与其他数值方法模拟的裂纹扩展速度进行了对比。另外，还有其他离散型方法也被广范应用与研究裂纹扩展问题。例如，Zhao等^[8]采用离散弹簧模型(the distinct lattice spring model, DLSM)模拟了自相似裂纹问题，从能量角度解释了在没有考虑损伤本构模型的方法中模拟结果与试验结果相吻合的原因。

数值流形法(numerical manifold method, NMM)^[9-10]既能够处理连续问题，又能够处理非连续问题，广泛应用于裂纹扩展问题的研究。然而NMM难以用于动态裂纹扩展研究。

非连续变形分析(discontinuous deformation analysis，DDA)最初由Shi^[11]提出，广泛应用于解决大尺度的岩石工程问题。一些学者采用DDA方法模拟材料破坏，实现从弹性阶段到破坏阶段全过程的模拟。王士民等^[12-13]为了研究节理岩体的失稳破坏问题，提出采用非连续子母块体理论，并采用增广拉格朗日乘子处理子块体接触。Ning等^[14-15]基于子块体DDA模拟了波传播、爆破等岩石动态破坏问题，并探讨了减小网格依赖性的解决方法。焦玉勇等^[16-17]通过设置虚节理将模型离散为DDA块体，并发展了新的接触本构模型模拟岩石破坏。与修正接触模型不同，徐栋栋等^[18]通过在子块体间插入黏结单元作为总体刚度矩阵的荷载项，实现应变软化的模拟。Xia等^[19]引入新的拉伸破坏准则，采用DDA模拟了巴西劈裂试验验证了算法的有效性，并分析模拟了地震滑坡算例。以上研究都采用了子块体DDA的方法来模拟岩石破坏，并以三角形离散的方法为主，从改进本构模型或破坏准则的角度进行裂纹扩展的定性分析。

与三角形离散相比，Voronoi多边形离散的DDA方法计算效率更高，一定程度上能够降低网格依赖性，其离散形式在模拟岩石材料时更加接近于晶粒结构^[20]。本文采用罚弹簧处理块体间的接触，由于罚弹簧取值太大会导致刚度矩阵病态，而罚弹簧取值较小，就必须考虑弹簧的变形量，因此为进行定量分析，需要对这一方法进行参数标定。本文首先基于Voronoi多边形离散的DDA方法提出均布弹簧刚度的算法使模拟材料破坏更加合理。然后，通过模拟单轴压缩和带切槽的巴西半圆盘试验进行参数标定，进而尝试使用Voronoi多边形离散的DDA量化研究脆性材料的动态裂纹扩展问题。

1.   Voronoi多边形离散的DDA方法改进

DDA方法通过Newmark时间积分求解块体系统的运动方程，其在解决动力学问题具有一定的优势。方法最初针对工程问题提出，将问题域按照岩体真实节理划分块体，节理参数按照实际参数取值。但是，由于原始DDA中块体本身无法发生破坏，因此不能很好地解决材料尺度的问题。为了使DDA能够模拟岩石等材料从弹性阶段到破坏的全过程，需要设置虚节理对材料进行离散，并设置虚节理强度，保证计算的连续性。当虚节理间的应力水平满足破坏准则时，虚节理变为真实节理，裂纹产生。本文采用最大拉应力准则和莫尔-库仑准则：

式中，$\sigma$和$\tau$是虚节理处的拉应力和切应力，${\sigma _{\text{t}}}$，c，φ分别为虚节理的抗拉强度、黏聚力和内摩擦角。

本文采用Voronoi多边形离散，其具有离散形式与岩石矿物晶粒结构相似，计算效率更高，在一定程度上能够降低其网格依赖性等优点^[20]。在模拟单轴压缩过程中发现，随机Voronoi离散试样在失稳前会优先在多边形短边处产生拉伸破坏，影响试样初期裂纹形态的发展。如图 1(a)为两个边长不等的块体1，2和块体3黏结。块体3上部施加速度荷载，模型底部固定。如图 1(b)，采用原始DDA计算，由于块体2的黏结边长l较小，则块体2会先于块体1发生破坏。这是由于原始DDA中的破坏判别与黏结长度(即边长的一半)直接关联：

图  1  不同边长块体的破坏模型及块体1，2的y方向位移–时间关系

Figure  1.  Failure model with different size blocks and y-direction displacement-time curves of blocks 1 and 2

下载: 全尺寸图片 幻灯片

式中${k_{\text{n}}}$，${k_{\text{s}}}$分别为法向、切向弹簧刚度；${d_{\text{n}}}$，${d_{\text{s}}}$分别为接触对法向位移、切向位移；φ为内摩擦角；c为黏聚力；l为黏结长度；${\sigma _{\text{t}}}$为弹簧抗拉强度。显然，受拉时${d_{\text{n}}}$相同，则块体2由于黏结长度小而先发生拉坏。剪切也有类似问题。这会造成Voronoi离散试样的网格依赖性增强。为减小Voronoi离散的DDA方法在模拟岩石破坏时的网格依赖性，本文对其进行了改进。

首先，确定模型的平均虚节理长度，通过计算均布弹簧刚度，从而进行破坏判断：

式中l_aver为模型内共边接触的平均黏结长度，k^*为均布弹簧刚度，l为当前接触的黏结长度。如图 1(c)所示，改进后，块体1和块体2的位移时间曲线完全重合，说明两者同时达到破坏条件，解决了上述问题。由于算法只是对破坏强度进行了修正，因此改进前后对算法的接触收敛性并没有太大影响。为了验证改进后的DDA方法在模拟岩石破坏中的有效性，进行了单轴压缩和巴西劈裂模拟，微观参数选取与先前研究相同，采用DDA动态计算，时间步长为10 μs，并统计了裂纹微破裂方向角。单轴压缩结果如图 2(a)所示，改进后的DDA在加载过程中，拉伸破坏更加集中在拉应力最大的方向，微破坏分布角度与虚节理长度的依赖性减小。如图 2(b)~(d)所示，随着加载时间增加，施加荷载提高，改进前后的单轴压缩微破裂角度分布趋于相同，但改进后DDA微裂纹分布更加合理。如图 2(e)，(f)所示，改进后的单轴压缩剪切破坏出现的更晚，但试样完全破坏后剪切破坏的数量和分布角度两者差异不大。图 3显示改进前后试样的破坏模式都具有明显竖向剥落的宏观拉伸裂纹，改进后的破坏试样宏观剪切破裂呈更加明显的X形。图 4，5显示改进后的DDA方法试样的裂纹分布比原始DDA更加集中于加载轴方向，但其破坏模式差异较小。通过单轴压缩和巴西劈裂模拟，进一步证明了这一改进在模拟脆性材料破坏时，能够一定程度上降低网格依赖性。

图  2  单轴压缩试样微破裂角度分布

Figure  2.  Distribution of microcrack angle of uniaxial compression specimens

下载: 全尺寸图片 幻灯片

图  3  单轴压缩试样破坏模式

Figure  3.  Failure patterns of uniaxial compression specimens

下载: 全尺寸图片 幻灯片

图  4  巴西劈裂试样微破裂角度分布

Figure  4.  Distribution of microcrack angle of Brazilian split specimens

下载: 全尺寸图片 幻灯片

图  5  巴西劈裂试样破坏模式

Figure  5.  Failure patterns of Brazilian split specimens

下载: 全尺寸图片 幻灯片

2.   细观参数标定

目前，采用DDA方法模拟材料破坏的研究逐渐丰富了起来，但很少有文献系统地研究过DDA方法的参数标定问题，大多数研究仍然停留在定性分析。本文针对Voronoi多边形离散的DDA方法进行参数标定，为后续定量分析提供基础。本文的参数表定工作分为两部分：通过模拟单轴压缩试验进行弹性参数(弹性模量、泊松比)的标定；通过模拟带切槽巴西半圆盘试验进行强度参数(虚节理的拉伸强度、黏聚力和内摩擦角)的标定。

2.1   弹性参数标定

由于离散数值方法的计算结果与离散程度有关，因此在标定前需要确定合适的块体尺寸。通过模拟单轴压缩的弹性变形阶段，确定块体数量对连续体的弹性模量、泊松比这类弹性参数的影响。单轴压缩试样尺寸为0.05 m×0.1 m，不考虑端部摩擦，在上下两个加载端设置4个位移点，分别施加0.01 m/s的速度荷载。试验微观参数选择如表 1所示，本小节只对弹性参数进行标定，而表 1中的强度微观参数则是通过3.2节的带切槽巴西半圆盘试验得到。由于岩石材料在细观尺度下矿物颗粒尺寸通常只有几毫米，因此本文综合考虑计算效率和精度，子块体尺寸分别选择1.0，1.5，2.0，2.5，3.0，4.0 mm离散试样进行单轴压缩模拟。如表 2所示，子块体尺寸在小于4.0 mm范围内，其宏观弹性参数相差并不大。因此后续数值试验可选择小于4.0 mm的块体尺寸进行离散。因此，在进行PMMA材料的弹性参数标定时，选择块体尺寸为2.0 mm的单轴压缩模拟试验。如表 3所示模型的宏观弹性参数与实际接近，作为后续试验的参数取值依据。

表  1  DDA微观参数取值

Table  1.  Values of DDA microscopic parameters

密度/(kg·m-3)	弹性模量/GPa	泊松比	内摩擦角/(°)	端部摩擦角	黏聚力/MPa	抗拉强度/MPa	法向弹簧刚度/GPa	切向弹簧刚度/GPa	最大位移比	时间步长
1190	4.4	0.38	30	0	70	7	13.2	5.28	1×10-5	0(自动计算时步)

下载: 导出CSV 

| 显示表格

表  2  块体数量对岩石材料弹性参数的影响

Table  2.  The influence of the number of blocks on the elastic parameters

块体尺寸/mm	块体数量/个	弹性模量/GPa	泊松比
1.0	6129	59.30	0.268
1.5	2672	59.05	0.264
2.0	1593	58.62	0.264
2.5	1051	59.10	0.264
3.0	701	58.61	0.269
4.0	430	58.82	0.273
注：Voronoi离散(R=0.3)。

下载: 导出CSV 

| 显示表格

表  3  PMMA和模型的宏观参数值

Table  3.  Value of macroscopic parameters of PMMA and model

类型	密度/(kg·m-3)	弹性模量/GPa	泊松比
PMMA材料	1190	3.24	0.350
模型	1190	3.28	0.357

下载: 导出CSV 

| 显示表格

2.2   强度参数标定

为了使用Voronoi多边形离散的DDA方法研究动态断裂问题，需要对DDA强度参数进行标定。本文采用带切槽的半圆盘(semi-circular bend，SCB)试验进行PMMA的混合模式断裂测试。如图 6(a)所示，模型尺寸与试验尺寸保持一致，半径R为50 mm，预制切槽长度a为15 mm，选择a / R = 0.3和S/R = 0.43进行纯I型和I-II型混合断裂试验。通过设置切槽角度β为0°，10°，20°，30°，40°，43°，47°，50°，模拟断裂模式从纯Ⅰ型断裂向纯Ⅱ型断裂渐变。模型采用尺寸渐变的Voronoi离散，最小块体尺寸为0.5 mm，集中在预制切槽两侧，从而降低裂纹扩展路径的网格依赖性，最大块体尺寸为3 mm，保证模型宏观弹性参数受到块体尺寸影响较小。

图  6  SCB的Voronoi多边形离散模型及正则化临界应力强度因子对比

Figure  6.  Voronoi discretization model of SCB and comparison of critical normalized stress intensity factors from experimental and numerical tests

下载: 全尺寸图片 幻灯片

对于SCB试样，其临界应力强度因子计算表达式为

式中Y_I，Y_II为SCB试样的几何参数，与a / R，S / R，β相关，Ayatollahi等^[21]给出了SCB在不同角度下Y_I，Y_II的取值。如图 6(b)，通过数值标定得到的正则化临界应力强度因子K_If/K_Ic和K_IIf/K_Ic与已有的试验结果对比，数值结果基本处于试验得到的范围内，且随着角度变化，Ⅰ型和Ⅱ型应力强度因子变化趋势与试验结果也较为吻合。其中，试验测得PMMA的Ⅰ型断裂韧度K_Ic为2.13 MPa·${{\text{m}}^{\frac{1}{2}}}$。表 4为数值标定的荷载值与试验荷载均值的对比，除30°时误差较大外，其他角度荷载误差均不超过10%。由于DDA方法本身属于离散型数值方法，结果受多个因素影响，这种误差可以接受。

表  4  数值标定荷载与试验荷载对比

Table  4.  Comparison of numerical and experimental load

裂纹角度/(°)	数值荷载/kN	试验荷载均值/kN	误差/%
0	2.28	2.38	4.20
10	2.47	2.53	2.37
20	2.28	2.45	6.94
30	2.53	3.03	16.50
40	3.54	3.73	5.09
43	3.58	3.63	1.38
47	4.14	4.13	0.24
50	3.97	4.24	6.37

下载: 导出CSV 

| 显示表格

图 7(a)中SCB试样的切槽倾斜角为0˚，在荷载作用下产生纯Ⅰ型裂纹，裂纹路径沿从初始预制裂纹沿着加载轴方向垂直延伸。图 7(b)~(h)随着SCB试样的切槽倾斜角β增大，裂纹均从预制裂纹尖端向载荷点扩展。因此，采用Voronoi多边形离散的DDA方法能够很好地模拟出裂纹扩展路径，并与已有的试验结果保持一致。

图  7  SCB模拟的裂纹扩展路径与试验对比

Figure  7.  Comparison of numerical and experimental crack propagation paths

下载: 全尺寸图片 幻灯片

综上所述，Voronoi多边形离散的DDA方法能够较为准确地模拟PMMA这类脆性材料的断裂行为。这部分标定工作是后续DDA方法模拟动态裂纹扩展的基础。

3.   数值试验

本节将基于标定后的参数(表 2)，用Voronoi离散的DDA方法模拟自相似裂纹扩展、裂纹分叉和紧凑拉伸试验，尝试从能量耗散、裂纹扩展形态、裂纹扩展速度等方面，定性和定量地分析动态裂纹扩展的模拟结果。

3.1   自相似裂纹动态扩展

建筑结构部件的初始缺陷与整个结构相比非常微小，在准静态拉伸应力荷载下，易发生整个结构失效。断裂力学中提出的自发裂纹扩展，即是这一概念在工程中的体现，其中Ⅰ型自相似裂纹属于自发裂纹的一种。自相似裂纹的扩展需要满足以下两个条件：稳定的力场和初始裂纹长度为0。Xia等^[22]通过试验对自相似裂纹进行了详细研究。与稳态裂纹不同，自相似裂纹的断裂能随着裂纹长度增加而增加，这与静态裂纹扩展相同。

图 8所示为0.3 m×0.3 m的Voronoi多边形模型，模型两端施加拉伸位移荷载U_y=0.05 mm，加载过程中保持准静态施加，计算采用DDA动态计算，时间步长为1 μs。在模型中部设置长为0.02 m的预设节理，预设节理的强度与其他虚节理强度相同，从而保持力场达到稳定前模型的连续性。位移荷载施加过程中通过设置动力系数为0.95来消耗多余动能，以保证模型内部应力平衡。在模拟过程中，当力场稳定后，动力系数恢复为1，取消阻尼作用，并将预设裂纹的强度设置为零，实现裂纹的突然释放，以此模拟自发裂纹产生时引发弹性波的传播。图 9为不同时刻下，裂纹突然产生引起弹性波的速度云图。由于施加的荷载较小，因此无法引发裂纹的继续扩展。在100 μs时，裂纹释放引发的弹性波已经到达上下两端边界。在200 μs时，裂纹引发的弹性波已经到达预制裂纹面。这一结果与Zhao等^[8]基于DLSM模拟自发裂纹的结果一致。

图  8  自相似裂纹动态扩展计算模型

Figure  8.  Dynamic crack propagation model of self-similar crack

下载: 全尺寸图片 幻灯片

图  9  预加载试样自发裂纹产生引发弹性波的速度云图

Figure  9.  Elastic waves induced by spontaneous crack in a pre-loaded specimen

下载: 全尺寸图片 幻灯片

为研究自相似裂纹的扩展规律，增大模型两端的位移荷载至U_y=0.18 mm，与前一算例相同，待模型中力场稳定后突然释放预制裂纹。当裂纹尖端的应力状态达到破坏强度时，裂纹将扩展，模拟结果如图 10所示。随着裂纹的扩展，主裂纹两侧的微裂纹数量增多，最后甚至产生明显的次级宏观裂纹(图 10(d))。

图  10  DDA模拟的自相似裂纹扩展

Figure  10.  Self-similar crack propagation predicted by DDA

下载: 全尺寸图片 幻灯片

为了分析不同加载位移下的裂纹扩展规律及网格规则性对结果的影响，采用了规则和不规则网格进行模拟，加载位移U_y=0.1，0.14，0.18 mm。其中规则Voronoi多边形与正六边形相同，只是将裂纹扩展路径预先设置为直线。如图 11为时间在200 μs时规则网格和不规则网格在不同位移荷载下的裂纹形态。对于规则Voronoi多边形(正六边形)离散模型，荷载增加其破坏形式均与图 11(a)相似：主裂纹沿着预制裂纹面进行扩展，并同样伴随有少量微裂纹产生。对于非规则网格，裂纹形态如图 11(c)，(d)所示，随着位移荷载增大，随机Voronoi多边形离散模型的主裂纹周围产生的微裂纹及次级宏观裂纹逐渐增多。这是由于块体和弹簧中存储的弹性能释放，导致主裂纹及其周围的微裂纹产生，并伴随着模型内块体动能的产生和模型的应力重分布。因此对于非规则网格，随着位移荷载增大，模型储存的弹性能越多，转化为模型的动能和微裂纹产生的数量越多，从而保证了自相似裂纹以恒定速度向两侧传播。当荷载足够大时，由于主裂纹两侧的微裂纹逐渐增多，会逐渐形成较明显的次级宏观裂纹，但主裂纹并不会消失。

图  11  规则网格的自相似裂纹扩展和不同预荷载下不规则网格的自相似裂纹扩展

Figure  11.  Self-similar crack propagation for both regular and irregular discretizations under different pre-loads

下载: 全尺寸图片 幻灯片

图 12(a)为随机Voronoi多边形离散模型的裂纹扩展长度与时间关系图。由于力场稳定，满足自相似裂纹的特征，其裂纹扩展长度与时间基本呈线性关系，即裂纹扩展速度是恒定的。位移荷载分别为U_y=0.1，0.14，0.18 mm时，裂纹扩展速度依次增大分别为637.06，700.94，851.55 m/s。图 12(b)为规则Voronoi多边形(正六边形)离散模型的裂纹扩展长度与时间关系图，位移荷载分别为U_y=0.1，0.14，0.18 mm时，裂纹扩展速度依次增大分别为795.35，860.65，990.77 m/s，同一位移荷载下裂纹扩展速度恒定的规律更加明显。由于规则离散模型产生的微裂纹相比不规则离散模型较少，因此更多的弹性能转化为裂纹扩展的动能，造成裂纹扩展速度较快。另外，由于规则模型的宏观参数和随机模型的宏观参数存在一定的差异，这也可能导致规则离散模型的裂纹扩展速度高于不规则离散模型。

图  12  裂纹扩展长度与时间关系

Figure  12.  Relationship between crack length and time

下载: 全尺寸图片 幻灯片

3.2   裂纹分叉

如图 13所示为0.04 m×0.1 m的矩形平板，并在轴线处设置长为${a_0}$= 0.05 m的预制切槽，在平板上下边界施加均布力荷载为${\sigma _0}$= 1 MPa。模型的微观参数依照标定的表 1进行取值。模型共有6180个多边形块体，在裂纹扩展路径方向进行了局部加密，块体尺寸为0.5~2 mm，选择DDA动态计算，时间步长为0.1 μs。

图  13  裂纹分叉数值模型

Figure  13.  Numerical model of crack branching

下载: 全尺寸图片 幻灯片

由于DDA中块体为常应变弹性块体，借鉴Zhou等^[23]计算弹性能时的简化公式，模型中每个块体中储存的弹性能为

式中，E为块体弹性模量，$\nu$为块体泊松比。通过改变施加荷载${\sigma _0}$的大小，可以改变每个块体中弹性能的大小。当块体存储的弹性能W₀全部转化为动态破裂能G_c时，Zhou等^[23]认为G_c会随着裂纹扩展速度增加而增加，当达到裂纹扩展速度极限值时，G_c接近于无限大，从而导致裂纹扩展速度达到极限值。如图 14(a)所示，PMMA材料在较高荷载作用下出现裂纹分叉的现象。在较高的荷载作用下，W₀值较高而当达到G_c时，裂纹就以相应较快的裂纹速度扩展。并且当G_c较高时，裂纹扩展速度达到极限值，此时发生裂纹分叉现象。

图  14  PMMA的裂纹扩展与分叉对比

Figure  14.  Crack propagation and branching in PMMA

下载: 全尺寸图片 幻灯片

如图 14(b)~(d)所示，试样在36 μs时发生裂纹分叉，60 μs时裂纹以两条宏观裂纹继续扩展，84 μs时两条主裂纹两侧出现大量分支裂纹和微裂纹。这是由于随着W₀增大，裂纹的G_c值也逐渐增大，裂纹扩展速度加快。当W₀无法完全通过G_c耗散时，便通过产生的分支裂纹和微裂耗散多余的能量，这与Zhou等^[23]的分析结果相吻合。图 15对比了DDA与扩展有限元(extended-FEM)^[1]、近场动力学(peridynamic)^[24]、黏结颗粒模型(the bond-particle method)^[7]、晶格法(the lattice-based approch)^[25]等其他数值方法模拟的裂纹扩展速度。对于DDA方法，当施加的荷载超过材料强度时，发生起裂现象。当裂纹扩展速度到达峰值，接近材料的Rayleigh波速938 m/s时，发生裂纹分叉。对比其他模拟方法发现，DDA的裂纹传播速度略高于其他数值方法，接近于Rayleigh波速传播。同时，检测到某些时刻的裂纹扩展速度高于Rayleigh波速，这是由于微裂纹出现可能先于主裂纹，导致监测到的波速过快，通过光滑后，可以避免出现这一问题，同样的现象也出现在Remmers等^[26]的研究中。

图  15  不同数值方法得到的裂纹扩展速度与时间关系对比图

Figure  15.  Comparison of predicted crack velocity versus time by different numerical methods

下载: 全尺寸图片 幻灯片

3.3   紧凑拉伸试验

为了进一步验证动态裂纹扩展与加载速率的关系，进行了紧凑拉伸试验的模拟研究。紧凑拉伸试验模型尺寸为0.2 m×0.2 m，凹槽尺寸为0.064 m×0.018 m。模型共3656个子块体，块体尺寸为0.5~4 mm，采用DDA动态计算，时间步长为1 μs。凹槽面一端为固定边界，另一端施加均匀的速度荷载。模型仍然采用PMMA参数，施加速度分别为0.2，1，3 m/s。在不同的速度荷载下，试样出现了不同的裂纹扩展形式。

如图 16所示，从上往下试验和模拟的速度荷载逐渐增大。当加载速度较小时，裂纹扩展形式接近于准静态加载的结果，即裂纹沿着试样中心扩展；当加载速度逐渐增大，裂纹逐渐向加载端方向倾斜，并当加载速度足够大时，裂纹会产生分叉现象，试样的断裂模式逐渐由I型断裂向I-II混合型转化，这一结果与已有试验及模拟结果类似^[27]。由于模型凹槽一端固定另一端加载，随着加载速度的增加，达到模型强度的应力水平逐渐向加载端靠拢，因此造成裂纹会逐渐向加载端的方向倾斜的现象。从能量耗散的角度分析，随着加载速度提高，模型中储存的弹性能逐渐增大，当加载速度达到一定值时，多余的能量通过产生大量微裂纹进行耗散，因此产生了裂纹分叉现象。

图  16  不同冲击速度下裂纹扩展路径的试验^[27]和模拟结果对比

Figure  16.  Comparison of experimental and simulation results of crack propagation path under different loading velocities

下载: 全尺寸图片 幻灯片

4.   结语

本文基于Voronoi多边形离散的DDA方法，通过单轴压缩和带切槽半圆盘断裂试验进行了参数标定。在此基础上，尝试定量分析脆性材料的动态裂纹扩展问题。

通过采用均布弹簧处理块体间的接触，修正后的方法在一定程度上能够解决大量短边优先破坏的情况，进一步减小Voronoi离散的网格依赖性；通过模拟单轴压缩试验和带切槽巴西半圆盘试验分别对弹性参数和强度参数进行标定，参数标定是后续定量分析的基础；基于Voronoi多边形离散的DDA能够模拟自相似裂纹扩展速度恒定，以及随荷载增加裂纹扩展速度增加的现象；对于裂纹分叉试验，DDA能够模拟裂纹扩展速度达到极限值时，发生裂纹分叉的现象，并且之后裂纹扩展速度发生下降；通过紧凑拉伸试验，验证了Voronoi多边形离散的DDA方法能够再现不同加载速度下，脆性材料产生不同破坏模式的规律。

本文研究中存在的主要问题包括：由于采用的本构模型没有考虑率效应的影响，后续需要在本构中考虑率效应的影响；本文的本构模型采用的是弹脆性本构，没有考虑损伤弱化，因此宏观损伤主要体现在微裂纹的增加；目前采用的强度准则为最大拉应力和摩尔库伦准则准则，导致在模拟裂纹扩展时，微裂纹有可能先于主裂纹出现。