1引言
　　自1960年Clough最先在进行飞机结构分析时提出有限元以来[1]，数值计算方法就开始在工程结构分析、流体分析及场分析等众多领域中得到广泛的应用，同时数值方法本身也得到不断的发展与完善。目前从种类上来说，主要包括有限元法、离散元法和流形元三类。这三类数值方法在实际应用中各有一定的优缺点，如有限元法在求解连续介质发生连续变形时非常有效，而离散元法在求解非连续材料发生非连续变形时非常有效。然而随着生产实践的发展，又提出了一些更新的问题如初始完整性很好的材料在外力作用下产生破坏的过程模拟等，由于有限元法是以连续性为基础，它在模拟材料破坏方面有着不可克服的困难，于是就采用一种近似的方法来模拟材料的破坏，即以某一物理量作为材料的破坏准则，如应力、应变或损伤等，然后确定一个阈值，当某单元的该物理量的值超过阈值时就认为该单元已经破坏，即把该单元从整个结构中删除，实际上是把单元的质量、刚度等乘以一个很小的数值，使该单元不参加以后的计算。这种方法在目前的计算中也得到了一定的应用[2-4]，但其显然是一种近似的方法，也可以说是有限元为弥补其在材料破坏模拟方面的不足而采用的一个补救措施，与实际情况存在着一定的差距。同样尽管离散元法在模拟材料的破坏方面有着巨大的优势，但是它却不能模拟新裂纹的产生，尽管目前也有不少学者采用虚拟裂纹的方法来模拟完整材料或部分完整材料破坏过程，并取得了良好的效果[5]，但显然这也是一种近似的方法。
　　于是流形元即数值流形方法(NMM)就应运而生，该方法是由美籍华人石根华博士继提出不连续变形分析方法(DDA)之后，又提出的一种新的数值方法[6]。虽然该方法是在20世纪90年代初才被提出，但由于其在模拟材料破坏方面具有明显的优越性，一开始就吸引了众多学者的关注。到目前为止，该方法已在材料破坏模拟分析中得到了成功应用[7,8]；同时该方法在块体运动模拟方面，完全吸收了DDA中关于块体运动的理论，能够很好地模拟块体破坏后的飞散过程，这一重大进展对以连续介质为基础的有限元法来说是一个重大的突破，克服了在利用有限元法计算时仅能给出待分析域内应力分布，而不能模拟其破碎及破碎后块体运动等情况的不足。目前关于流形元模拟材料破坏、运动等方面已有不少的相关研究成果发表。但是至今还没有看到关于有限元与流形元的相关对比研究成果发表，如果能够采用这两种方法对同一问题进行分析研究，则更能够很好地反映出它们之间的相互关系及各自的优缺点，因此本文就采用这两种方法对岩石冲击破坏过程进行模拟研究，并对其模拟结果进行了对比分析。
　　2两种数值方法的动力学求解格式
　　岩石在冲击载荷作用下的破坏问题属于明显的断裂动力学问题。断裂动力学就是研究那些惯性效应不能忽略的断裂力学问题[9]，其求解方法明显不同于断裂静力学问题，它们的一个最重要的区别就是材料的惯性效应不能忽略，考虑了惯性效应的断裂力学就是断裂动力学或动态断裂力学。在动态加载时，试件除产生弹塑性变形外，内部各质点的自由振动要获得一定的加速度，从而产生惯性力，这就是所谓的动态加载时的惯性效应。
　　在动力问题的分析中，数值流形方法引入了惯性矩阵，它相当于有限元方法中的质量矩阵，以充分考虑动力学问题中的惯性效应。数值流形方法在处理动力学问题时，与处理静力学问题相比的一个最大差别就是在当前步的计算中，各单元继承了前一时间步的速度，而不是像处理静力学问题一样置当前时间步的速度为零。而这种处理方法对一般动力学问题来说过于简单，并且这种处理方法也仅仅考虑了单元的惯性效应，而没有考虑单元的阻尼，所以还不够精确。因此这里采用的数值流形方法计算程序是借用动力有限元法的求解思想，利用动力有限元中的Newmark方法改进后的数值流形方法计算程序，改进的原理及改进后的程序优越性已在相关文献中进行了论述[7]。有限元的动力计算方法同样是采用瞬态动力分析Newmark方法。
　　3数值流形方法的计算原理
　　由于目前数值流形方法在实际中的应用还不是太广泛，所以有必要首先对其计算原理进行简单介绍。
　　3.1动力问题的求解方程
　　流形方法在求解动力学问题中所使用的方程为结构动力学方程：

　　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所示。