摘 要：流形元是新出现的一种被较为广泛地应用于材料破坏模拟的数值分析方法。为了验证该方法在材料破坏模拟中的有效性，分别利用流形元与有限元两种不同的数值方法对岩石冲击破坏过程进行了模拟分析。模拟结果表明，流形元法对材料破坏的模拟结果更符合实际情况，这主要是由于该方法采用了两种不同的覆盖系统——数学覆盖和物理覆盖，并引入能够客观反映材料裂纹的产生与扩展准则，以及相应的块体运动理论。流形元的出现有望对材料破坏模拟开创出一条新的途径。
关键词：岩石冲击破坏；数值模拟；流形元；有限元；对比分析
 
1   引 言
自1960年Clough最先在进行飞机结构分析时提出有限元以来[1]，数值计算方法就开始在工程结构分析、流体分析及场分析等众多领域中得到广泛的应用，同时数值方法本身也得到不断的发展与完善。目前从种类上来说，主要包括有限元法、离散元法和流形元三类。这三类数值方法在实际应用中各有一定的优缺点，如有限元法在求解连续介质发生连续变形时非常有效，而离散元法在求解非连续材料发生非连续变形时非常有效。然而随着生产实践的发展，又提出了一些更新的问题如初始完整性很好的材料在外力作用下产生破坏的过程模拟等，由于有限元法是以连续性为基础，它在模拟材料破坏方面有着不可克服的困难，于是就采用一种近似的方法来模拟材料的破坏，即以某一物理量作为材料的破坏准则，如应力、应变或损伤等，然后确定一个阈值，当某单元的该物理量的值超过阈值时就认为该单元已经破坏，即把该单元从整个结构中删除，实际上是把单元的质量、刚度等乘以一个很小的数值，使该单元不参加以后的计算。这种方法在目前的计算中也得到了一定的应用[2-4]，但其显然是一种近似的方法，也可以说是有限元为弥补其在材料破坏模拟方面的不足而采用的一个补救措施，与实际情况存在着一定的差距。同样尽管离散元法在模拟材料的破坏方面有着巨大的优势，但是它却不能模拟新裂纹的产生，尽管目前也有不少学者采用虚拟裂纹的方法来模拟完整材料或部分完整材料破坏过程，并取得了良好的效果[5]，但显然这也是一种近似的方法。
于是流形元即数值流形方法(NMM)就应运而生，该方法是由美籍华人石根华博士继提出不连续变形分析方法(DDA)之后，又提出的一种新的数值方法[6]。虽然该方法是在20世纪90年代初才被提出，但由于其在模拟材料破坏方面具有明显的优越性，一开始就吸引了众多学者的关注。到目前为止，该方法已在材料破坏模拟分析中得到了成功应用[7,8]；同时该方法在块体运动模拟方面，完全吸收了DDA中关于块体运动的理论，能够很好地模拟块体破坏后的飞散过程，这一重大进展对以连续介质为基础的有限元法来说是一个重大的突破，克服了在利用有限元法计算时仅能给出待分析域内应力分布，而不能模拟其破碎及破碎后块体运动等情况的不足。目前关于流形元模拟材料破坏、运动等方面已有不少的相关研究成果发表。但是至今还没有看到关于有限元与流形元的相关对比研究成果发表，如果能够采用这两种方法对同一问题进行分析研究，则更能够很好地反映出它们之间的相互关系及各自的优缺点，因此本文就采用这两种方法对岩石冲击破坏过程进行模拟研究，并对其模拟结果进行了对比分析。
 
2  两种数值方法的动力学求解格式
岩石在冲击载荷作用下的破坏问题属于明显的断裂动力学问题。断裂动力学就是研究那些惯性效应不能忽略的断裂力学问题[9]，其求解方法明显不同于断裂静力学问题，它们的一个最重要的区别就是材料的惯性效应不能忽略，考虑了惯性效应的断裂力学就是断裂动力学或动态断裂力学。在动态加载时，试件除产生弹塑性变形外，内部各质点的自由振动要获得一定的加速度，从而产生惯性力，这就是所谓的动态加载时的惯性效应。
在动力问题的分析中，数值流形方法引入了惯性矩阵，它相当于有限元方法中的质量矩阵，以充分考虑动力学问题中的惯性效应。数值流形方法在处理动力学问题时，与处理静力学问题相比的一个最大差别就是在当前步的计算中，各单元继承了前一时间步的速度，而不是像处理静力学问题一样置当前时间步的速度为零。而这种处理方法对一般动力学问题来说过于简单，并且这种处理方法也仅仅考虑了单元的惯性效应，而没有考虑单元的阻尼，所以还不够精确。因此这里采用的数值流形方法计算程序是借用动力有限元法的求解思想，利用动力有限元中的Newmark方法改进后的数值流形方法计算程序，改进的原理及改进后的程序优越性已在相关文献中进行了论述[7]。有限元的动力计算方法同样是采用瞬态动力分析Newmark方法。
 
3  数值流形方法的计算原理
由于目前数值流形方法在实际中的应用还不是太广泛，所以有必要首先对其计算原理进行简单介绍。
3.1 动力问题的求解方程
流形方法在求解动力学问题中所使用的方程为结构动力学方程：
 式中：M为质量矩阵；C为阻尼矩阵；△δ为位移增量；分别是位移速度和加速度；K=Ke+Kcn+Kcs+Kf，其中Ke是刚度矩阵，Kcn、Kcs分别为块体及不连续面之间的接触矩阵，Kf是约束矩阵。△F为总载荷增量矩阵，△F=Fp+Fb+Ff-F0+Fcn+Fcs+Ffr，其中Fp是外载荷向量，Fb是体积力向量，Ff是已知约束位移引起的等效载荷向量，F0是初应力向量；Fcn、Fcs分别为法向和切向接触引起的等效载荷向量；Ffr，为接触面之间的摩擦力引起的等效载荷向量。
3.2 裂纹的产生及扩展机理
在外力作用下，材料的破坏过程一方面是由于材料内部已存在的裂纹被激活扩展，而另一方面就是在外载作用下，材料内部的应力超过了材料的抗拉或抗剪强度而导致材料被拉断或剪坏。而对于岩石等抗压强度远远大于其抗拉和抗剪强度的脆性材料来说，其破坏的主要形式为拉伸破坏和剪切破坏。
材料的破坏过程是裂纹在外力作用下不断产生、聚集和汇合的过程。在流形程序判断新裂纹的起裂和已存在裂纹的扩展问题中，采用了不同的标准。对于已存在的裂纹是否扩展的问题采用断裂力学的应力强度因子作为判断标准，根据线弹性断裂力学来判断裂纹的起裂与扩展[7]。对于新裂纹的产生，一般考虑采用以应力为判据。这里采用摩尔-库仑准则作为裂纹起裂的判据，其判断方法如下：①当第1主应力大于材料的抗拉强度时，材料被拉伸破坏，产生新裂纹；②当某点的最大剪应力大于材料的抗剪强度时，材料被剪切破坏，产生新裂纹。
3.3 大变形及接触问题的处理
在外力作用下，岩体内的初始裂纹将发生扩展并最终把岩体切割成相应的块体，这些块体在自重或其他外力作用下产生运动。数值流形方法对块体的运动模拟中，一定要保证块体之间无嵌入、无拉伸，它完全借鉴了DDA中块体运动学的理论。而这一过程是通过反复的开一合迭代，即在可能的接触边界之间施加和移去弹簧来实现的。在边界处块体之间有张开、锁定和滑动三种状态，当边界从张开发展到锁定时施加切向和法向弹簧；当边界从滑动到锁定时加切向弹簧；当边界之间滑动时在边界两侧加摩擦力，从而对刚度矩阵作出修正，计算接触力并叠加于本次计算的平衡力列阵中，参加下一次的迭代计算。
 
4   算例分析
4.1 爆炸冲击破坏的力学特征及荷载简化
在本算例中，取爆炸载荷作为施加的冲击载荷。相应的实际问题是对应于二维平面问题的圆形装药在有限域内爆炸后形成爆破漏斗的过程模拟。炸药在炮孔中起爆后，岩石将发生如下的破碎过程：
(1)强大的冲击波压应力使炮孔周围岩石受压破碎，在瞬时形成压缩破碎和初始裂隙。
(2)环向拉应力及应力波反射拉应力使岩石中的裂隙扩展，引起岩石进一步破裂，包括初始裂隙的扩展和二次裂隙的形成。
(3)爆生气体膨胀作用使岩石中的裂隙贯穿形成碎块，岩块运动形成爆破漏斗。
由于实际的爆破过程十分复杂，因此在满足工程要求的条件下，可以对岩石爆破破坏的过程进行简化处理，通过在炮孔内壁上施加均布的冲击三角波载荷来模拟爆炸荷载，冲击载荷的升压时间为80ms，降压时间为220ms，整个荷载作用时间为300ms，载荷的最大峰值压力为5GPa。
4.2 计算模型及计算结果分析
所建立的计算模型是二维分析域内的圆形装药在有限域岩体内爆炸后，形成爆破漏斗的过程。取计算模型的尺寸为4m×1.5m，圆形炮孔直径为O.5m，炮孔中心距上边界的距离为0.5m，距左右边界的距离均为2m。假定岩石为各向同性均质岩体，其中模型的上边界为自由边界，其余三个边界均为固定边界，计算模型示意图如图1所示。
在本算例中，所分析问题的类型为动态类型，因此所采用的材料计算参数均应为动态参数。取大理岩的参数作为本算例中的岩石参数：动弹性模量80GPa，动泊松比0.20，密度2400kg／m3，动态断裂韧性0.5MN／m3/2，节理面的摩擦角45°，黏结力10MPa，动抗拉强度30MPa。
4.2.1流形元计算结果
取计算过程中的3个瞬时状态如图2(a)～(c)所示。
4.2.2有限元计算结果
采用目前应用较为广泛的有限元数值计算软件ANSYS对上述模型进行计算分析，模拟结果如图3（略）所示。由于有限元软件还不能很好地考虑裂纹的产生及块体的运动，所以采用了有限元法中对材料破坏模拟的近似方法——单元生死法，即采用某一准则对计算结果中的单元进行删除，即删除失效的单元。这里采用Mises有效应力失效准则，当某单元的Mises应力超过1MPa时，就认为该单元已经失效，并把该单元从计算模型中删除，用于表示岩石的破坏；同时还给出了模型中的应力分布情况。
4.2.3两种计算结果的对比分析
由上述流形元和有限元模拟的结果可以看出：
(1)由于流形元采用的是材料的拉、剪破坏准则和线弹性断裂力学理论，因此它可以很容易地模拟完整岩石在冲击荷载作用下新裂纹的起裂及扩展；对于块体的运动采用的足DDA中的块体运动理论，所以它也能够很好地模拟块体的运动过程。而相比之下，由于有限元法是以连续介质为基础，它不具备考虑岩石中新裂纹的产生及扩展的功能，所以也很难得出符合岩石实际破坏过程的计算结果。但有限元为了克服其在模拟材料破坏方面的不足，引入了单元删除的方法，而图3(B)所给出的结果就是采用单元删除的方法得到的岩石破坏情况计算结果，这在一定程度上弥补了有限元在这方面的不足。
(2)从对岩石破坏过程模拟的准确性来看，岩石的破坏是一种渐进的过程，也就是说岩石的破坏是从无裂纹到一条裂纹再到多条裂纹，以至于最后的完全破坏。图2的流形元模拟结果非常客观地反映了这种破坏过程，即材料首先产生一条裂纹，然后再到多条裂纹，最后切割成块体。相比之下，有限元的计算结果则无法反映这种客观的破坏过程，在O.7ms以前岩石还处于少数几条裂纹产生与扩展阶段时，有限元的计算结果则几乎无法反映这种少量、局部的破坏情况；而当破坏的范围逐渐增大，即在3ms以后有限元则以应力的形式体现出了岩石中的应力不同，进而可以根据应力失效准则得出材料可能的破坏区域，但有限元仍无法模拟出岩石破坏后的块体运动情况。因此，从对岩石破坏过程模拟的准确性来看，流形元法则显得更为优越。
(3)从对爆炸空腔的模拟结果来看，流形元法得到的爆炸空腔是一个倒三角形的形状，这与岩石在近地表处的炸药爆炸破坏的实际结果是符合的，这是因为在近地表处由于自由面的影响，裂纹很容易扩展达到地表，形成爆破空腔。有限元模拟的爆破空腔则是在炮孔周围较大，而在地表则没有形成一个倒三角形的爆破漏斗。造成两者差别的根本原因是两种数值方法采用的算法不同，由于流形元采用的是断裂力学方法，当岩石中出现初始裂纹以后，由外力作用产生的能量将很容易使已有裂纹继续扩展，这样就很容易形成图2(c)中的两条优势裂纹。当这两条优势裂纹达到地表以后，就形成了所谓的爆破漏斗，即使这两条优势裂纹中间的岩石块体没有得到破碎也丝毫不影响漏斗区的形成。而有限元方法则不同，由于其没有裂纹的形成，所以也不存在能量的优先分配，所有的能量几乎都是按照单元距离或位置对称的原则进行分配，所以其形成的爆破空腔也是非常对称的，显然这种数学上的对称是不符合岩石破坏的实际情况的。
(4)从裂纹的发展情况来看，流形元在模拟岩石的冲击破坏过程中裂纹的发展，存在着明显的分岔现象。即裂纹并不是沿着最初出现的裂纹一直向前发展，而是在发展过程中会偏离原来的路径而分为两条或多条分支裂纹，甚至在原裂纹没有到达的地方出现了新的小裂纹，这是裂纹动态扩展中一个最常见的现象之一[10]。而有限元的模拟中却没有发现这一点，这也是由于有限元的局限性所致。
 
5   结 语
本文对目前应用最普遍的有限元和最新出现的流形元两种数值方法，在计算岩石冲击破坏中的优缺点进行了对比分析，说明了流形元在材料破坏模拟中的优势，即由于流形元是根据材料的实际破坏特点采用相应的裂纹起裂及扩展准则进行计算的，因此对材料破坏过程的模拟与实际情况是比较吻合的。而有限元法则是以连续介质为基础，没有考虑到材料的实际破坏过程，而是以应力或应变等失效准则作为判断材料破坏的判据，因而计算出的材料破坏情况与实际有较大的差距。因此，以连续介质为基础发展起来的有限元法在计算连续材料在外力作用下产生连续变形方面是具有很大优势的，而在计算不连续破坏方面则具有一定的局限性。由于流形元法是以物理覆盖和数学覆盖双重覆盖为基础建立的，因而能够方便地考虑材料在外力作用下的破坏过程，破坏后单元的重新构建也十分方便，同时由于其在块体运动方面完全借鉴了DDA中的块体运动理论，因此对材料破坏后的运动模拟也是十分有效，有望对材料的破坏模拟研究开创出一条新途径。
摘自《工程爆破》总第61期
参考文献：
[1]Clough R W．The Finite Element in Plane Stress Analysis [A]．In：Proc．2ndASCE Conf．on Electronic Computation [C]．1960．
[2]Zhi-liang Wang，Yong-chi Li，R．F．Shen．Numerical Simulation of Tensile Damage and Blast Crater in Brittle Rock Due to underground Explosion[J]．International
Journal of Rock Mechanics & Mining Sciences。2007，44：730—738．
[3]尹潘．岩石破坏机理的LCA模型数值模拟[J]．辽宁工程技术大学学报，2006，25(3)：364—366．
[4]王来贵，赵娜，周永发，等．岩石受拉破坏的数值模拟方法[J]．返宁工程技术大学学报，2007，26(2)：198—200．
[5]张国新，武晓峰．裂隙渗流对岩石边坡稳定的影响——渗流、变形耦合作用的DDA法[J]．岩石力学与工程学报，2003，22(8)：1269—1275．
[6]Shi Genhua．Manifold Method of Material Analysis[A]．In：Transactions of the Ninth Army Conference on Applied Mathematics and Computing[C]．Minneapolish，Minncsoda，USA：1992，51—76．
[7]刘红岩，秦四清，杨军．爆炸荷载下岩石破坏的数值流形方法模拟[J]．爆炸与冲击，2007，27(1)：50一56．
[8]张国新，赵妍，石根华，等．模拟岩石边坡倾倒破坏的数值流形方法[J]．岩土工程学报，2007，29(6)：800—805．
[9]范天佑．断裂理论基础[M]．北京：科学出版社，2003．
[10]刘再华，解德，王元汉，等．工程断裂动力学[M]．武汉：华中理工大学出版社，1996．