论文导读：：在连续体动力问题中心差分算法的基础上，引进准脆性材料的破坏准则，节点单元的破裂算法、离散子块的接触搜索及接触力计算等，对准脆性连续体在冲击载荷作用下的破裂破坏问题进行数值模拟。通过数值算例，给出结构在冲击载荷作用下裂纹产生和扩展的模拟结果图，初步验证程序的正确性和可应用性，为模拟连续体转变为非连续体这一复杂物理过程提出新方法和新思路。

　　论文关键词：有限元，离散元，准脆性，破裂，冲击

　　1 引言

　　准脆性材料结构如混凝土、岩石等在外力作用下产生的损伤断裂等破坏过程涉及到从连续体向非连续体转化的复杂物理过程，其破坏过程的研究对于结构工程、地下工程以及采矿工程等都有重要的意义。而对这一过程的模拟在数值建模、算法等诸多方面存在很多困难，迄今尚未完全成熟。

　　由于离散元法等用于散体分析的数值方法在处理结构破碎后的动力分析时有很大的优势，所以很多研究者采用离散元法与其它方法结合的算法来模拟结构多重破裂的动力学行为[1-7]。另一类很重要的方法是从连续体方法（有限元）出发进行研究[8-16]。英国Swansea大学Owen领军的研究团队[8-11]、Munjiza[12-13] 和唐春安[14]等做了大量的研究。如Klerck[11]研究的混凝土冲击问题的数值模拟；Klerck[11]和Cottrell构建了准脆性材料在拉压应力场中破碎过程的力学模型，伴随着材料的破坏同时产生了肉眼可见的裂纹。徐泳与Swansea大学合作，研究了土壤破碎问题，并与试验结果对比[15]；

　　本文方法可称为可离散有限元法，主要是结合传统有限元及离散元方法的优点，在连续体有限变形理论和动力问题中心差分算法的基础上，首先定义出连续固体，建立合理的破裂破碎准则，然后施加外来作用令其产生破裂和破碎，将有限元法与断裂-损伤力学融合在一起用于分析准脆性固体的破碎。具有固体特征的准脆性材料连续体在破裂破碎前后形态和本构特征都用动力学有限元规范和断裂-损伤力学进行描述，计算破碎后的子块碰撞作用时要引入离散元法的接触搜索和接触规律处理等。

　　2 可离散有限元法基本理论

　　2.1 动力有限元法

　　有限元动力学方程可以表达为

　　（1）

　　式中分别是位移、速度和加速度向量。M是质量矩阵，在中心差分法中一定采用集中质量矩阵，使得M成为对角矩阵。C是阻尼矩阵，一般采用Rayleigh阻尼（或广义Rayleigh阻尼）的形式，使得C可用M与K表达，例如表达为。

　　C = a1M + a2K。（2）

　　a1和a2是两个常系数，K是刚度矩阵，但是在中心差分法中并不需要组集成总刚度矩阵。F是外载荷向量。把方程改写成

　　（3）

　　方程右端是可以在单元级组集的列向量离散元，方程左端是一个对角矩阵乘以加速度向量，所以上式可以简化成

　　mii = fi ，（i=1,2,……）

　　速度的中心差分可表示为：

　　破裂（4）

　　将（4）式中的步长减小一半，就可有

　　类似有加速度的中心差分公式

　　破裂

　　当加速度求得后，就可求得速度

　　或可写成

　　破裂（5）

　　再由速度求得位移

　　当从运动方程解得节点的位移向量后，则可利用几何方程和物理方程计算所需要的应变和应力。

　　2.2 准脆性材料的应变软化失效模型

　　试验表明，一般混凝土材料在应力达到抗拉强度之前表现为线弹性特性，超过这个极限，随应变的缓慢增加拉伸应力逐渐减小。一维的拉伸应变软化表示为：

　　（6）

　　上述表达式也可由图1来表示。

　　在弹性加载过程中，应力与应变是线性关系，即图1（a）中的直线段部分，该直线的斜率是弹性模量E。当应力达到材料的抗拉强度ft时，开始出现应变软化，软化曲线是一个单调的、斜率为负的连续函数。最终，假设在某一有限的应变或时应力达到零。当时为加载；当时为卸载中国期刊全文数据库。

　　（a）全程应变（b）非弹性应变

　　（a） Fullstrain（b） Inelastic strain

　　图1 一维拉伸应变软化模型

　　Fig.1 one-dimensional tensile strain softening model

　　模拟计算时，为简单起见，对准脆性材料可以采用线性软化模型，则在第一主应力方向上的抗拉强度

　　（7）

　　其中，ft是初始抗拉强度，H可由下式给出：

　　（8）

　　其中，为单元等效尺寸，Gf为材料特定的断裂能。对于任意形状的有限元网格，应变软化会在一个单元内产生，若该单元的等效尺寸为，则断裂能可表示为：

　　（9）

　　其中，为该单元更新后的断裂能密度。

　　2.3 单元的拓扑更新

　　当拉伸应力下降为零时，在与其垂直的方向上引入裂纹。定义失效系数为断裂应变与临界断裂应变的比值。在高斯点k的局部失效系数Fk可表示为：

　　（10）

　　其中， Fk的断裂方向为qk（为垂直于失效软化方向）。将该时刻邻近单元的失效系数和断裂方向在节点上进行加权平均，来得到全局的失效映射。节点p处加权平均后的失效系数和失效方向表示为：

　　（11）

　　其中，Nadj为邻近高斯点个数，wk为权重系数，通常取为单元面积或体积。

　　当某一节点的失效系数大于1，则在该点处引入一条裂纹。在定义了失效节点的断裂面以及加权平均后产生的断裂方向后，裂纹以及相关的拓扑更新就会被引入。裂纹面代表了真实的裂纹方向且通常情况下不会恰巧与单元边界重合（如图2（a））。因此，实际的裂纹扩展方向需要进行局部的网格更新，这就可能产生新的单元，这种裂纹生成方式叫做单元内断裂（如图2（b））。如果只是单纯地将单元一分为二，该方式可能会导致几何形状狭长的不理想新单元产生，这会减小计算的时间步长。除此之外，还有另外一种裂纹生成方式，即单元间断裂方式（如图2（c）），该方式不进行局部的网格更新，裂纹只在与预计的断裂方向最接近的单元边界产生离散元，因此该方式的网格会影响裂纹的走向，对网格划分精度要求较高。

　　（a）加权平均后的结点失效方向

　　（a）Node failure direction afterweighted average

　　（b）单元内断裂

　　（b）Rupture in elements

　　（c） 单元间断裂

　　（c）Rupture between elements

　　图2 裂纹生成方式

　　Fig.2Crack formation way

　　当裂纹被引入之后，单元的损伤状态变量由于裂纹的形成而被恢复为初始值。理想的情况是只有垂直于裂纹的方向上的单元损伤被恢复为零，但是若有网格更新而生成了额外的新单元，则所有失效单元的状态变量全部初始化为零。

　　3 数值算例

　　假设冲击物为准脆性材料，被冲击物强度很高，冲击速度100m/s，计算模型的各项参数如表1所示，图3（a）给出计算用初始的网格剖分。图3（b）-（f）给出了准脆性材料冲击物裂纹开始产生及扩展的过程。由图可见，裂纹从接触处开始产生，随着时步的增加，裂纹扩展，由小的裂纹慢慢扩展成大的裂缝，接下来会有离散子块产生。

　　表1 准脆性材料冲击模型参数

　　ρ（kg/m3）

　　密度

　　E（GPa）

　　弹性模量

　　Ν

　　泊松比

　　ft（MPa）

　　抗拉强度

　　Gf（J/m2）

　　断裂能

　　p0

　　罚系数

　　冲击物

　　2340

　　26

　　0.18

　　3.15

　　3.0e+2

　　10e+9

　　被冲击物

　　3.0e+4

　　210

　　0.3

　　1.0e+14

　　3.0e+12

　　10e+9

　　the parameters of quasibrittle material impact model

　　table1

　　（a） （b）

　　（c）（d）

　　（e）（f）

　　图3 准脆性冲击问题裂纹产生及扩展图

　　Fig3 cracking and extension figures of quasi brittle material impact problem

　　（a）初始网格剖分 （b）裂纹开始产生（c）裂纹扩展 （d）- （f）开始产生较大的裂缝

　　（a） Initial mesh （b） Began to crack （c）crack propagation （d）-（f） Began to produce larger crack

　　4 结论

　　基于有限元动力问题中心差分算法，给出用于模拟冲击载荷作用下准脆性连续体多重破裂问题的可离散有限元法，针对准脆性材料的结构可以分裂、破碎的特点，建议采用不建立刚度矩阵的中心差分法求解。在准脆性材料只产生裂纹但没有达到破碎分离和不研究二次破碎问题时，线性小变形理论也还是可以采用的。但是为了描述分离子块的运动、碰撞及二次破碎问题，由于小块一般均有大位移、大转动，应该采用有限变形理论。由数值计算，我们发觉裂纹的发生、发展以及破裂现象对网格剖分是敏感的，但是只要网格剖分的足够密，裂纹的的发展趋势大致相同。

　　参考文献

　　1.OňateE, Rojek J. Combination of discrete element and finite element methods fordynamic analysis of geomechanics problems[J]. Computer methods in appliedmechanics and engineering, 2004, 193（27-29）： 3087~3128

　　2.YuanqiangTan, Dongmin Yang, Yong Sheng. Discrete element method（DEM） modeling offracture and damage in the machining process of polycrystalline SiC[J]. Journalof the European Ceramic Society, 2009, 29（6）： 1029~1037

　　3.HoonLee, Ji Hoe Kwon, Kwan Ho Kim, Hee Chan Cho. Application of DEM model tobreakage and liberation behavior of recycled aggregate from impact-breakage ofconcrete waste[J]. Minerals Engineering, 2008, 21:761~765

　　4.邢纪波，俞良群，王泳嘉。三维梁—颗粒模型与岩石材料细观力学行为模拟。岩石力学与工程学报，1999,16（6）：627~630.（XingJ B, Yu L Q, Wang Y J. 3-D Beam-Particle Model for Simulating Meso-mechanical Behaviorof Rock Material[J]. YanTu LiXue Yu GongCheng XueBao/Chinese Journal of RockMechanics and Engineering, 1999,16（6）：627~630 （in Chinese））

　　5.胥建龙，唐志平。离散元与有限元结合的多尺度计算方法及其应用。计算物理，2003（20）：477-482.（Xu K L, Tang Z P. Combined Discrete/Finite Element MultiscaleNumerical Method and Its Application[J]. JiSuan WuLi/Chinese Journal of ComputationalPhysics, 1999,16（6）：627~630（in Chinese））

　　6.傅华，刘仓理，王文强等。冲击动力学中离散元与有限元相结合的计算方法研究。高压物理学报，2006,20（4）：379~385.（FuH, Liu C L, Wang W Q. A Combined Discrete/ Finite Element Method in Shock Dynamics[J].GaoYa WuLi XueBao/ Chinese Journal of High Pressure Physics, 2006, 20（4）：379~385 （in Chinese））

　　7.侯艳丽，张冲，张楚汉等。拱坝沿建基面上滑溃决的可变形离散元仿真。岩土工程学报，2005，27（6）：657~661.（Hou Y L, Zhang C, Zhang CH. Simulationof upward-sliding failure of interface in arch dams by deformable distinctelements[J]. YanTu GongCheng XueBao/ Chinese Journal of Geotechnical Engineering, 2006, 20（4）：379~385 （inChinese））

　　8.Owen DR J, Feng Y T, Klerck P A, Yu J. Computationalstategies for discrete systems and multi-fracturing materials[J].Engineering Mechanics and Computation, 2001:135-146

　　9.Owen DR J, Feng Y T. Parallelised finite/discrete element simulation of multi-fracturingsolids and discrete systems. Engineering computations, 2001,18（3,4）： 557-576

　　10.CottrellM G, Yu J, Owen D R J. The adaptive and erosive numerical modeling of confinedboron carbide subjected to large-scale dynamic loadings with element conversionto undeformable meshless particles[J]. International Journal of ImpactEngineering, 2003,28（9）：1017-1035

　　11.KlerckP A, Sellers E J, Owen D R J. Discrete fracture in quasi-brittlematerials under compressive and tensile stress states[J]. Computer methods inapplied mechanics and engineering, 2004, 193（27-29）： 3035-3056

　　12.MunjizaA, Owen D R J, Bicanic N. A combined finite-discrete element method intransient dynamics of fracturing solids[J]. International Journal ofComputational Engineering Science,1995,12:145-174

　　13.Munjiza A. The combined finite-discrete elementmethod[M]. New York: John Wiley & Sons. 2004

　　14.唐春安，王述红，傅宇方，岩石破裂过程数值试验，[M].北京：科学出版社，2003 （Tang Chunan, WangShuhong, Fu Yufang. Numerical Test of Rock Failure. Beijing: Science Press,2003 （in Chinese））

　　15.YanjieLi, Yong Xu, Wenbing Huang,Y. Feng, D R J Owen. Experimental and numericalresearch on bulldozer working process[J]. Chinese Journal of MechanicalEngineering, 2007, 20（2）：41-46