脆性材料具有“硬且脆”的特点，可以承受很高的外加载荷作用而不发生屈服和塑性变形，但是一旦外加载荷超过其弹性极限，则通常会发生迅速的裂纹扩展贯通和灾难性断裂。典型的脆性材料包括岩石、陶瓷、玻璃、固体炸药等，在强动载荷的爆炸与冲击加载下易发生断裂和破碎，造成功能失效和事故危害。举例来说：在利用爆破方式进行矿山掘进时，爆炸波与复杂地质条件的综合作用易导致岩爆等矿难事故；铁电陶瓷可以在冲击波压缩下瞬时释放脉冲大电流，但是冲击破坏也易导致其介电击穿和放电失效；防弹玻璃广泛应用于各类军事车辆，但遭受大口径弹药的打击时，其对车内人员的防护能力仍显不足；固体炸药填装于各类武器弹药之中，但是当经受意外的跌落、撞击时，炸药易受到损伤并发生爆炸。深入开展脆性材料动态断裂机制的理解和断裂行为的调控研究，既属于力学研究领域的前沿挑战，也对军事、工业等领域中提升脆性材料的可靠性和安全性具有重要意义。

断裂问题是力学界的百年难题，众多学者通过实验、理论和数值模拟对裂纹扩展、介质破碎等现象开展了研究。在脆性材料内单裂纹的拉伸扩展研究中，裂纹扩展路径的振荡和分叉问题得到了持续关注^[1–3]：从能量转换角度看，周围介质注入到裂纹尖端的能量越多，裂纹就以越快的速度扩展从而消耗能量；如果直线开裂不足以完全消耗注入的能量时，裂纹就会发生振荡^[4]；当注入裂尖的能量进一步增加，裂纹将从振荡演化到分叉^[5]，针对一些特定材料，利用裂纹扩展的两个临界速度可以较好地表征发生路径振荡与裂尖分叉的起始条件。另一方面，超高速碰撞条件下弹丸和靶板的破碎以及碎片云的形成过程也得到了较深入的研究^[6–8]：在相关的实验与数值模拟中，通常无法追踪单个碎片的断裂行为，为此仅对大量碎片整体的破碎尺寸、速度分布、飞散过程、对后靶的毁伤效应等进行评估。爆炸与冲击加载注入的能量通常不足以使脆性材料彻底破碎，但足以驱动大量裂纹同时扩展，形成裂纹网络。相比于过程清晰的单裂纹以及统计描述的碎片云而言，这类裂纹网络的动态演化问题目前仍缺少深入的探索和系统的研究。裂纹网络研究的难点在于：实验测量中现有的高速摄影、X射线照相、质子照相等原位诊断手段均不具备足够高的空间分辨能力以追踪每条裂纹的扩展；但理论分析上裂纹的密度通常又没有达到适于作统计描述的程度；数值模拟中应用最为广泛的有限元、有限差分等计算方法在处理爆炸与冲击加载导致的剪切裂纹随机萌生、裂纹面挤压摩擦等问题时仍面临较大的算法困难。

实验、理论、数值模拟多方面的困难共同制约了对爆炸与冲击加载下岩石、陶瓷、玻璃、固体炸药等各类脆性材料动态断裂行为的深入研究。为此，相关领域的研究者们建设新装置、提出新理论、发展新方法，以期获得对裂纹网络动态演化机理的系统认识，进而实现对脆性材料可靠性和安全性的设计调控。本文将聚焦数值模拟方面：首先介绍多种发展或改进用于脆性断裂模拟的介观尺度计算方法，对比分析其优缺点；然后，系统介绍在描述脆性材料的断裂和破碎行为方面具有天然优势的一种粒子类方法—格子模型；接着列举了将格子模型应用于脆性断裂研究的一系列代表性成果；进而探讨现有格子模型存在的不足，指出将来的发展和改进方向；最后，给出总结与展望。

1.   典型的脆性断裂介观模拟方法

介观尺度的模拟计算方法既可以描述材料的微结构和缺陷对每条裂纹扩展过程的调制，又可以统计获得大量裂纹演化对材料宏观尺度力学响应的影响，因而非常适合于爆炸与冲击加载下脆性材料的动态断裂模拟研究。本节分为两部分，分别介绍两大类介观模拟方法，即依赖网格的计算方法（有限元方法等）和无网格/粒子方法（物质点法、近场动力学方法、格子模型等）在脆性材料动态断裂研究领域的发展现状。

1.1   依赖网格的计算方法

1.1.1   有限元方法

有限元方法（Finite Element Method，FEM）是目前最成熟且应用最广泛的力学响应计算方法，其基本原理是利用网格将连续的材料划分成有限数量的多边形小单元，并假设每个单元内部介质的变形状态均可以通过对单元顶角上“节点”变形状态的插值计算得到^[9]。这样整个材料中每个点上位移、应力、应变状态的求解就简化为网格上有限数量节点的位移、应力、应变状态求解。基于材料的本构关系（应力与弹性系数、应变、应变率、温度等之间的函数关系），结合外界施加的始边值条件（飞片碰撞速度、周期性边界或自由边界等条件），可以迭代求解出材料中各个节点上的粒子速度、应力、应变等。在有限元模型中还可以考虑温度场、电磁场等与材料的相互作用，实现多物理耦合计算^[10]。虽然有限元方法的优点众多，但是依赖网格的计算方法普遍存在不适合模拟大变形和损伤断裂现象的缺点，导致其在脆性断裂研究中的应用较少。

Espinosa等^[11–13]基于有限元方法建立了多晶陶瓷的冲击波压缩模型。在一个长条形样品中，利用“Voronoi Tessellations”技术划分出大量形貌随机的多边形，每个多边形代表一个晶粒，在陶瓷晶粒之间，还设定气孔、微裂纹和玻璃相等微结构/缺陷。每个晶粒内部都被三角形网格划分成多个有限元单元，两两晶粒之间采用无厚度的四边形“内聚力”单元进行连接。四边形单元在断裂判据满足后将发生断裂，然而三角形单元却无法断裂，因此该模型只能模拟陶瓷的沿晶断裂，不能表现穿晶断裂。在冲击波从上至下扫过样品的压缩过程中，当剪应力导致的晶粒间剪切变形超过断裂判据后，晶界上开始出现裂纹。

1.1.2   流体动力学方法

流体动力学方法原本应用于流体物质的计算，但是在极高的冲击应力下固体材料抵抗剪切变形的能力（屈服强度）与冲击波施加的剪切应力相比可以忽略不计，此时固体也像流体一样无法保持自身形状，因此强冲击条件下流体动力学方法也常用于固体冲击响应的计算。流体动力学方法中考虑的控制方程与固体力学中不同，但通常也用有限元或有限差分进行求解。Bourne等^[14]利用多介质欧拉流体动力学程序研究了AD995氧化铝多晶陶瓷的冲击响应。模型样品的晶粒依照实验样品冲击前的显微形貌照片进行划分，实际样品中晶粒内部和晶界上的气孔均被呈现在模型样品之中。为简化影响因素，模型中采用最简单的Murnaghan状态方程和弹性-理想塑性本构关系。模拟中发现，随着冲击波向前传播，原本平整的波阵面由于样品的介观不均匀性而变得参差不齐。在Hugoniot状态中大部分晶粒都已经进入塑性变形阶段，但是也有少量晶粒正好取向到弹性极限最高的晶向，仍处于弹性变形阶段，导致宏观波剖面上出现上、下两个弹塑性转变点。该模型只能根据应力场的分布推测缺陷周期可能形成沿晶和穿晶裂纹的位置，并不能真正模拟出裂纹扩展过程；在孔洞周围应该存在的应力集中效应在模型中也没有被明显地表现出来；由于没有裂纹扩展和应力集中，孔洞塌缩压实所需要的应力被大大提高，以至于模型中没有观察到孔洞塌缩现象。

有限元方法等基于连续介质力学框架而建立起的计算方法具有理论严密、数值精度高、发展成熟和应用广泛等优点，但其依赖于网格的求解方法，在处理大变形、断裂和破碎等复杂问题时会遇到诸如网格畸变扭曲、控制方程中存在对流项、断裂面上方程无法求导、物质界面不易确定等各种数值困难^[15]。因此，相关领域研究人员的努力目标是发展新的不依赖于网格的计算方法，以克服有限元法的缺点，但同时又在新方法中保留有限元方法的诸多优点。

1.2   无网格/粒子方法

无网格/粒子方法中也存在网格的概念和定义，但其算法的设定回避了网格畸变、断裂所带来的各种数值困难，近年来得到充分重视和快速发展。目前，已经发展出多种理论框架和计算方法，广泛地应用于动态裂纹扩展、金属加工成形、高速和超高速碰撞、流体动力学等问题研究^[15]。然而，当前适用于研究脆性材料中动态裂纹扩展，特别是爆炸与冲击加载下裂纹网络动态演化的方法仍较少报道，本节将介绍其中具有代表性的模型。

1.2.1   物质点法

物质点法（Material Point Method）最初于1994年由Sulsky等^[16–17]提出并应用于塑性大变形问题研究。其基本思想是将要研究的样品离散成为模型中的一组质点，质点的相对运动就代表了样品的变形。基于密度集中的概念，每个质点都集中携带其周围一定区域内介质的质量，模型中所有需要被记录的状态量都定义在质点之上。在模拟中的每个时间步上都做一次映射：首先将质点所携带的信息（如质量、材质、位置等）映射到一张新建的网格上；然后在网格中计算出当前时间步的应力和应变状态等，这里的网格计算直接调用成熟的有限元算法；接着再做一次反向映射，将网格中算出的信息转换为质点速度等信息；最后根据质点速度和时间步长（两个时间步之间的间隔）算出下个时间步中质点的位置。进入下一步计算时，重新映射生成新的没有发生畸变和扭曲的网格。这样一套迭代计算，既充分地发挥了有限元求解的优点，又完全回避了有限元的各种数值困难。物质点法被广泛应用于超高速碰撞、侵彻/穿甲等塑性大变形问题的模拟^[18–19]。但是，大部分物质点法对断裂问题的模拟较为定性，力学中经典的拉伸型、剪切型裂纹扩展模式及动态拉伸所产生的层裂现象等均需要进行适当的设定后才能被模拟^[15]。

Li等^[20]将物质点法应用于脆性圆盘经受碰撞之后断裂特征的模拟。他们对模型的改进之处在于引入Weibull的断裂强度可变性理论作为断裂判据。在每个时间步中计算所有质点上的最大正应力$\sigma$，并利用（1）式计算在一定体积V_ref之中介质的“生存概率”${P_{\rm{s}}}\left( {{V_{{\rm{ref}}}}} \right)$

式中：${\sigma _0}$为参考应力，m为Weibull系数。受破坏的质点被称为“幽灵质点”（Ghost Points），它们不再承载任何应力和应变，但是携带质量以保证质量矩阵不发生变化。“生存概率”的断裂判据是将物质点法应用于准确描述断裂现象的初步尝试之一。生存概率${P_{\rm{s}}}\left( {{V_{{\rm{ref}}}}} \right)$的设定带有很强的经验性；“幽灵质点”相当于一个空洞，不适合描述爆炸与冲击加载下剪切断裂过程中的两侧介质相互挤压。目前，将物质点法应用于断裂问题研究是学者们的努力方向之一^[21–22]。

1.2.2   近场动力学方法

近场动力学（Peridynamics）方法最初于2000年由Silling^[23]提出并致力于描述材料断裂问题的新方法。在传统的连续介质力学框架下，裂纹作为一个间断面，会导致运动方程的不连续，因而无法通过求导获得应变、应力等信息。此时裂纹界面必须被处理成为一种边界条件；但材料在复杂加载下产生裂纹的位置又是不可预知的，因此在模拟过程中实时地确认“内部边界”（裂纹）的位置和状态就成为困难的数值问题^[24]。近场动力学的核心思路就是在求解作用于质点（节点）的力时，将求导过程反过来变成积分过程；积分并不要求函数的连续性，这就从根本上克服了裂纹间断面带来的数值困难。基于“键”的近场动力学在计算方面类似于分子动力学，模型中利用“对势”计算某质点受到的周围一定“范围”（Horizon）内其他质点的作用力^[23]。与分子动力学的不同在于，近场动力学本质上是一种宏观连续介质力学方法，而非微观方法，模型中对势和范围的选择均不需要微观上的具体物理意义；分子动力学中的原子并不会“记忆”自身的初始位置和变形过程（即迭代过程中这些信息并不参与计算），但近场动力学中的质点会记忆参考构型。

Ha等^[25–26]利用键基的近场动力学研究了脆性材料中一条张开型（I型）裂纹在向前传播过程中随速度不同而出现的分叉、不对称分叉、连续分叉、次生裂纹（与主裂纹成90°角）等复杂现象，成功地再现了实验中观察到的几种典型裂纹分叉现象，并能根据模型中的初始和边界条件、应力场分布和演化过程等解释这些现象的形成机制和条件。研究发现，控制这些丰富现象的主要机制在于应变怎样从加载边界传递到裂纹的尖端区域。最初提出的键基近场动力学还存在一些不足之处。为克服这些问题，Silling等^[27]提出一种经过推广的态基近场动力学。目前，近场动力学的应用正从均匀的各向同性材料拓展到非均匀的各向异性材料之中^[28]。由于计算效率不及有限元，在出现裂纹的间断区域采用近场动力学，而在无裂纹的连续区域采用有限元方法是将来的发展方向之一^[29]。

1.2.3   格子模型

格子模型（Lattice Model）的起源可以追溯到早期的几类典型模型，如框架网络方法（Frame Network Method）^[30]、凝聚态物理中模拟原子间谐振势能的方法^[31–33]、1979年Cundall^[34]提出的离散元方法（Discrete Element Method）等。格子模型并非一开始就具备完整的理论框架和规范，而是在固体材料的弹性和断裂性质研究中逐渐发展成熟。其核心优势在于清晰的物理图像而非严密的数学框架，当其被不同的研究者应用于不同领域时，模型的具体设定都会被相应地修改，以便突出该领域中最关键的影响因素^[35–36]。这就导致格子模型存在大量的变体，且文献中对其称呼也多种多样，如格点-弹簧模型（Lattice-Spring Model）、离散元方法、弹簧网络模型（Spring Network Model）等。

格子模型与近年来逐渐成熟的近场动力学方法非常类似，在处理断裂和破碎问题时具有天然优势。与近场动力学方法的不同在于：（1）格子模型中格点间的相互作用通常是局域的，即一个格点只与其直接相邻的格点相互作用（但也存在相互作用非局域的格子模型）；（2）对于格点间相互作用的计算，研究者们提出了多种处理方法，以表现不同材料的力学特性；（3）格子模型中格点速度和位移更新与分子动力学方法相同，仅取决于当前状态，而不需要用到初始的参考构型。格子模型的物理图像清晰，算法简单易行，已广泛应用于研究具有复杂微结构的非均质材料中发生的断裂和破碎现象。

2.   脆性材料动态断裂的格子模型

2.1   基本原理

如图1所示，格子模型及其各种变体都是由格点/颗粒和弹簧/键/梁组成的网格结构，其中：格点携带介质的质量、位置、速度等信息；弹簧则决定相互作用力、加速度、阻尼、应力和损伤断裂等；网格描述格点与弹簧的拓扑结构与排布规律，但网格本身并不参与计算。格子模型的计算过程与分子动力学方法^[37]非常相似：在当前时间步下，根据格点的相对位置和弹簧的设定（作用方式和作用力参数），计算出每个格点的加速度；根据当前时间步的格点速度和加速度，通过运动学积分得出下一时刻格点的位置和速度；根据格点的位移算出应变，根据弹簧中作用力的方向和大小算出应力，再由应力和应变等判断弹簧是否断裂并形成微裂纹。模型在外加载荷下发生的变形、断裂等演化代表和反映了被模拟的真实材料在外加载荷下将会出现的变形和断裂演化。

图  1  格子模型示意图

Figure  1.  Schematic of the lattice model

下载: 全尺寸图片 幻灯片

2.2   相互作用的设定

在格子模型中，格点间的相互作用通常被形象地理解为由于连接两格点的弹簧发生了拉、压变形而导致的张力和斥力。弹簧相互作用力设定最简单的格子模型，被称为“Hookean模型”^[38–40]。它是一个中心力场模型，其中的格点只受到中心力（只与两格点间相对距离有关的力）的作用，其基本形式可以用下式表示

式中：f_ij、k_ij和$\Delta {l_{ij}}$分别为格点j对于格点i的作用力、两个格点间连接弹簧的刚度系数及相对于初始平衡距离的弹簧伸长量。Hookean模型中两两格点之间并不会表现出抵抗剪切变形的效果，但模型整体上是抗剪的。

中心力场的设定导致Hookean模型只具有固定的泊松比，而不能模拟各种泊松比的材料。为克服这一局限，多种改进模型被提出，其中最常用的是在两两格点之间增加剪切弹簧的“Born模型”^[41–43]，其形式可以表示为

式中：上标n和$\tau$分别表示法向和切向分量，${\delta _{ij}}$表示格点j与格点i之间的相对位移。这里法向定义为穿过两个格点，并从格点i指向格点j的方向；切向则垂直于法向。（3）式中的法向相对位移等同于（2）式中的弹簧伸长量。在加入切向作用力以阻碍两两节点间的剪切变形之后，Born模型就可以通过调节切向的弹簧刚度系数$k_{ij}^\tau$控制模型整体的泊松比。其他相互作用力设定还包括在一个格点的两个相邻法向弹簧之间增加一根“环向弹簧”，用于阻碍两根法向弹簧之间的夹角变化^[44]，从而起到与增加切向弹簧一样的效果，获得泊松比可调的模型。在两两格点之间还可以设定弹性力学概念中的“梁”，以同时具备抵抗拉伸、压缩、剪切和弯曲的作用^[45]。总之，相互作用的设定，可以根据所研究问题的具体特点进行相应的改进。目前既能表现不同的泊松比，形式又最为简单的Born模型得到了最多的应用。

Born模型的一大缺点是普遍不具备应变能的转动守恒性质^[39]。这一缺点在研究小变形量问题时并不突出，但却对本文所关注的动态断裂、裂纹网络扩展、介质破碎等现象的模拟产生影响。2011年Zhao等^[46]分析了导致Born模型转动不守恒的原理，并提出了“基于局域应变的方法”以改进切向相对位移的计算，避免引入转动相关项，成功地恢复了Born模型的转动守恒性质。喻寅等^[47]则通过在普通的Born模型中引入一个额外的转动自由度，使得Born模型自动获得应变能的转动守恒性质。

2.3   作用力参数的确定

相互作用方式的设定定性地决定了模型的力学特性，而作用力参数的选取则进一步定量地决定模型的响应。这是格子模型中最重要和最困难的问题，以往的很多工作都只通过经验或半经验方式设定作用力参数，限制了模型定量计算的能力。理论上严密的参数设定方法一直是格子模型研究中关注的中心问题，为此研究者们提出了很多有效的方法^{[39, 46, 48–50]}。其中，应用最广泛的是基于Hertz接触理论的参数设定方式^[34]，其法向作用力F_n的计算为

式中：$\delta$、E^*和R^*分别为两个颗粒之间的法向重叠量、等效杨氏模量和等效半径。E^*和R^*分别定义为

式中：${E_i}$、${\nu _i}$、${R_i}$ 和${E_j}$、${\nu _j}$、${R_j}$分别为颗粒i和颗粒j的杨氏模量、泊松比和颗粒半径。这一算法在描述散体/粉体等非连续介质的挤压作用时最为有效。

针对连续介质的相互作用，1996年Ostoja-Starzewski等^{[44, 51]}报道了一套基于“应变能相等”原则而建立的参数设定方法。模型中每个格点通过法向弹簧与周围6个最近邻格点相连，假设周围6个格点相对于中心格点的位移已知，则根据连续介质力学计算得到的应变能为

式中： ${\sigma }$和${\varepsilon }$分别为应力和应变矢量，在体积V内对其积分就得到应变能；C表示材料的刚度张量，${\sigma } = {C} \cdot {\varepsilon }$；由于讨论的是线弹性和均匀变形状态，故可以直接完成积分。将所有弹簧中存储的弹性势能相加也能得到这一微元中的总弹性势能

式中： ${E_b} = {\left( {k{u} \cdot {u}} \right)^b}$表示第b根弹簧上的势能，$k$和u分别表示弹簧的刚度系数和位移差，${N_b}$为弹簧的总数。如果模型能够正确地反映真实材料的变形状态，那么必然有

进一步地，应变${\varepsilon}$和位移差u之间可以相互表示。例如可以将位移差变换为应变，此时模型微元的弹性势能为

式中：新矩阵k含有弹簧刚度系数的信息，可以利用k与材料的刚度张量C的对应关系求出弹簧的刚度系数。类似的原理也被其他研究者应用于计算和设计各种不同相互作用方式的格子模型^{[39, 46]}。

2004年Gusev^[52]提出了以有限元为媒介进行格子模型参数映射的方法，其基本思想是：首先建立一个格子模型和有限元共用的模型网格，然后将待模拟材料的弹性常数转换为有限元网格中的相互作用参数，再借助共用的网格，映射出格子模型中的相互作用参数（弹簧刚度系数）。在有限元中，各个单元之间的位移和作用力是相互关联的，为了得到整个网格的变形状态，所有单元的作用力–位移关系需要被组装在一起并整体求解^[9]。设节点总数为N，则有

或

式中：F和${\delta}$分别表示整个网格中每个节点受到的外力和位移，K被称为“总体刚度矩阵”（总刚）。有限元方法通过求解（11）式获得模型的整体变形，但是在格子模型的计算中，却需要单独计算两两格点之间的相互作用。所以在参数映射过程中需要先分解总刚，获得节点j对于节点i的作用力，如（12）式所示

式中：${\delta _{ij}} = {\delta _i} - {\delta _j}$。在格子模型中，颗粒i与颗粒j之间的相互作用正好由（12）式进行描述，从而保证有限元和格子模型计算得出相同的作用力和应变能。有限元方法是公认的具有严谨数学基础和良好计算精度的算法，意味着经过参数映射的格子模型具有定量表现材料弹性变形响应的能力。

2.4   断裂判据的选取

在格子模型中最常用的断裂判据设定为：当弹簧的相对位移量、主应力、剪应力等指标超过某个临界值之后发生弹簧断裂，该临界值可参考实验数据并结合经验进行设定。另一类常用的断裂判据则基于Griffith能量平衡原理^[53]，即当模型中某一对弹簧（法向和切向）中存储的弹性势能（机械能）大于或等于形成一段长度为$c$的微裂纹所消耗的能量${R_0}c$时，就判定弹簧满足断裂条件。其中，“附着功”R₀表征各种阻碍裂纹扩展机制所消耗的总能量，可由材料的断裂韧性导出。在喻寅等^[54]开发的格子模型中，将法向弹簧处于拉伸状态时的弹性势能与切向弹簧处于剪切状态时的弹性势能相加，再将总弹性势能与${R_0}c$进行比较，判断是否发生断裂。他们认为，纯压缩（不带剪应力，例如均匀介质的静水压压缩）并不会引起材料破坏，所以在计算总弹性势能时并不计入法向弹簧处于压缩状态时的能量。当断裂判据满足后，弹簧将永久性地断开，两个颗粒之间失去抗拉伸和抗剪切的作用力。但是受破坏的颗粒之间仍然有其他相互作用，当它们的间距小于平衡距离时，将产生沿法向的排斥力和黏性力以及沿切向的干摩擦力。

能量判据的优点在于可以方便地统一判定拉、压、剪不同加载状态下法向和切向弹簧的断裂行为。但是在对线弹性弹簧使用能量判据时，模型的断裂应变等存在网格尺寸依赖性：在某一固定的网格尺寸下，从断裂表面能可以导出相应的断裂应变；但是当固定断裂表面能不变时，随着网格尺寸的减小，模型的断裂应变将不协调地增大。这一问题在均分网格中可以通过经验性的设定规避，但在使用局部区域加密网格的模型时会导致样品各区域的断裂特性不一致。其根本原因在于，脆性材料虽然在宏观上表现出从线弹性直接进入脆断的响应，但是微观尺寸上原子之间从拉伸到断裂的响应却是非线性的；如果强行用线弹性弹簧统一描述从宏观到微观的响应，必然造成模型失真。针对线弹性弹簧采用能量判据时面临的缺陷，张振南等^[55]提出一种新的两段线性弹簧以描述拉伸断裂。这一设定综合了分子动力学作用势和有限元内聚力模型的特征及优点，可以通过一套参数校验方法，确保弹簧拉伸断裂所对应的断裂能、断裂应变、断裂强度与实验测量结果相符。

2.5   微结构与缺陷的设置

材料的宏观响应由基质的力学性质和微结构共同决定。本节所介绍的格子模型的相互作用力形式、弹簧刚度系数设定、断裂判据选取等都确保了模型可以准确地反映材料基质的力学性质。再依据材料的微介观特征，采用适当的网格，引入适当的介观结构，就可以模拟出真实脆性材料的动态断裂过程。如图2左侧所示，格子模型中可使用规则（结构化）与随机（非规则）两大类网格定义格点与弹簧的相对位置和排布规律：规则网格适用于陶瓷等晶体材料，有利于表现材料的各向异性，而随机网格适用于玻璃等非晶材料，有利于表现材料的无序特性；规则网格对裂纹的扩展方向产生显著的调控效应，有利于表现晶体的解理断裂、剪切滑移等特征，而随机网格不限制裂纹的扩展方向，有利于表现非晶的翅裂纹、扩展路径不稳定性等特征。

图  2  格子模型中网格、微结构、缺陷设定示意图

Figure  2.  Representatives of spring networks, microstructures and defects in the lattice model

下载: 全尺寸图片 幻灯片

脆性材料中常见的微结构和缺陷包括晶界、孔洞、第二相颗粒、微裂纹等。图2右侧给出了几种典型陶瓷材料的微结构和缺陷建模示意。例如对于致密的透明陶瓷，建模时主要考虑其晶粒随机取向及两两晶粒之间的晶界；建模时常利用“Voronoi Tessellations”技术实现晶粒多面体的划分^[11]，并对每个晶粒内的规则网格进行随机旋转表现出取向的随机性；如此建模后可观察到剪切裂纹穿过晶界后的偏折和分叉。对于疏松的多孔陶瓷，建模时还需要在其中设定孔洞；孔洞通常是在样品中随机选择圆心点，然后删除圆心周围给定半径内的一簇格点；持续选点、删除，直到满足预设的气孔率。对于颗粒增强复合陶瓷，建模时需要在基体材料中加入第二相颗粒材料，两种材料将分别设定不同的变形响应和断裂判据；两种材料之间的弹簧通常取两组参数的平均值。对于玻璃、塑料等非晶材料，则应该在随机网格之上进行微结构/缺陷设定。对于岩石、混凝土等由大量晶体颗粒组成的材料，如果注重研究其大尺度力学响应，则既可采用随机网格进行建模，也可采用规则网格建模，并为弹簧设定随机的参数。

3.   格子模型研究脆性断裂的代表性成果

3.1   国外研究的代表性成果

Zapperi等^[56]运用一种最简单的格子模型—“随机电熔丝模型”，研究了非均质脆性材料的拉伸断裂演化过程。该模型以弹簧参数的随机性表现石膏、混凝土等非均质材料内部的无序特性；基于标量弹性理论中胡克定律（$\sigma = E\varepsilon$）与电学中电流、电导率、电压关系（$I = \sigma V$）的数学形式相同的特征，利用电学方程的求解算法计算格子模型中每根弹簧上的应力负载。当一根弹簧上的负载达到最大值时，不是将弹簧断开，而是使其弹性系数按一定比例减小。这种情况下总应变的增加与部分弹簧弹性系数的减小相互平衡，使得总应力稳定在一个平台上，相当于出现了宏观塑性。这一设定的考量在于当材料的一个小区域中出现微裂纹时，其影响可能不是使周围介质完全失效，而是使小区域中的力学性质受到减损，即相当于使弹簧弹性系数按一定比例减小。这套模型设定表现出了脆性材料断裂过程中的声发射、自组织演化和雪崩行为。

Kale等^[57–58]利用格子模型研究了非均质脆性材料中的弹性–塑性–脆性转变与雪崩行为。其格子模型中的弹簧被设定为弹性–塑性硬化响应：弹性段的刚度系数$k = 2E/3\left( {1 + \nu } \right)$，塑性段的刚度系数${k_{\rm{T}}} = 2{E_{\rm{T}}}/3\left( {1 + \nu } \right)$，其中E、E_T、$\nu$分别为弹性模量、硬化模量和泊松比；格点间作用力F与伸长量 $\delta$ 的关系为

式中：${\delta _{\rm{y}}}$、${\delta _{\rm{f}}}$、${\delta _{\rm{u}}}$分别为弹簧发生弹塑性转变、拉伸断裂与准弹性卸载时的伸长量。为表现材料内部的无序性，弹簧参数进行了随机设定。研究发现，弹性参数的变化显著地调控脆性材料的宏观应力–应变响应，可以模拟出从弹性–塑性硬化到弹性–脆性断裂的多种响应特征；对无序程度的设定可以控制理想塑性与塑性硬化材料中的塑性应变雪崩行为。

Ostoja-Starzewski等^[59]利用格点正三角形排布，但弹簧的杨氏模量和拉伸强度设有随机起伏的格子模型，研究了预制孔洞的环氧树脂板在准静态拉伸下的裂纹随机扩展行为。模拟中设定了不同网格尺寸的数值试样，并参照对应的实验样品在试样中设定位置无序排列的31个直径为1/4英寸的圆孔，再以2.5～40.0 cm/s的不同恒定速度对试样施加准静态拉伸。模拟发现，在格点间距为0.1～0.01 cm的4种不同网格中，基本不存在裂纹扩展路径的网格依赖性，均成功重现了在实验观察中占主导的裂纹特征。针对格点间距为0.02 cm的网格，研究了环氧树脂的无序特征对裂纹扩展路径的影响，对比发现当随机起伏的变化范围设定在[–0.5×10^–4, 0.5×10^–4]时，杨氏模量随机性的影响比拉伸强度随机性的影响更强；但当随机起伏的变化范围设定在[–0.5×10^–8, 0.5×10^–8]时，杨氏模型与拉伸强度随机性的影响均可忽略不计。

Krajcinovic等^[60]利用随机三角形网格的格子模型研究了脆性板被高速物体侵彻时的空腔膨胀特性。他们将弹簧设定为拉伸时以线弹性响应，而压缩时则以非线性响应，即

式中：${\lambda _{0ij}}$和${\lambda _{ij}}$分别为初始与当前时刻格点i与格点j之间的距离，${k_{ij}}$为弹簧刚度，B定义了排斥墙的斜率。非线性的相互作用使得格点在相互靠近时的排斥力能够更迅速地增长，导致冲击波速度和介质抗压缩能力增加。针对刚性弹体侵彻脆性材料所导致的空腔膨胀过程，研究了中心空腔以不同速率向四周膨胀时样品中损伤的演化特征，如空腔周围（包括压碎区、损伤区、弹性区）的损伤度、损伤分布、扩展速率等。

Buxton等^[39]利用正六面体网格的三维格子模型研究了带有夹杂物的材料中弹性和塑性变形的特征。在模型中考虑了格点与其最近邻格点（棱边顶点）和次近邻格点（面对角线顶点）的相互作用，并基于系统应变能与弹簧势能相等的要求推导获得了杨氏模量、泊松比、体弹模量等与弹簧法向及切向系数的关系。为描述塑性变形，假设材料为各向同性硬化；对于已进入塑性屈服阶段的弹簧，减小其弹性模量，但保证应力变化的连续性。模拟表明：在弹性变形条件下，模型可以准确表现各向同性材料的变形行为；针对材料中存在球形夹杂物的情况，模型计算的弹性应力场分布与理论预测符合较好；在考虑基体材料的塑性变形后，模型计算的球形夹杂物周围应力场分布与Wilner等^[61]的预测符合良好。

Horie等^[48–49]利用考虑晶界结构与晶粒取向的格子模型，研究了冲击波加载下多晶铜的冲击塑性响应和拉伸层裂行为。为将格子模型应用于延性金属材料模拟，Horie等在法向弹簧中引入Lennard-Jones势（（15）式），以描述冲击加载下材料的非线性压缩行为。

式中：${ {f}}_{\rm{p}}^{ij}$为格点i与格点j之间的排斥力，r为两格点间当前时刻的距离，r₀为格点间的平衡距离；其余参数a、m、n根据材料的Hugoniot冲击压缩数据进行拟合；在拉伸阶段，引入了温度相关的屈服条件和理想塑性响应。对于切向弹簧，设定为弹性–理想塑性响应。模型中由一组正三角形排布的格点构成一个晶粒，大量晶粒紧密排列构成整个数值试样；每个晶粒中的格点均设定一个排布方向，以表现多晶材料中的晶向随机性。模拟中观察到了冲击加载导致的晶粒转动变形，分析了弹性波、塑性波的展宽行为，探讨了层裂破坏及准弹性再加载现象的内在机制。

Wang等^[62–63]利用格子模型的一种变体—离散元模型，研究了准静态单轴加载下脆性材料样品中的翅裂纹演化和压碎破裂过程。该模型中：考虑了格点间的法向、切向作用力及滚动与扭转的力矩；利用四元数描述格点在三维空间中的旋转；通过将格点间的相对旋转分解为两次顺序无关的转动，唯一地确定相对旋转产生的力矩；基于有限变形的方式，而非增量方式计算力和力矩，获得了更好的数值稳定性。利用二维情况下格点粒径随机分布的网格模拟了翅裂纹萌生、扩展的过程；发现只有当法向、切向与转动作用都存在时，模拟计算才能重现实验中观察到的翅裂纹特征。利用三维情况下格点以面心立方排布的网格模拟了岩石的压缩破裂行为，成功再现了可在岩石单轴压缩破坏实验中观察到的主裂纹快速贯通过程，及崩裂飞散后两个加载端面上各自残留的未受破坏的岩石尖锥特征。

Wang等^[64]将离散元与有限元相耦合，研究了地下岩石开采中的爆破过程。模拟中利用有限元商业软件计算爆炸及气体产物膨胀过程，将爆炸载荷及气体压强作为边界条件施加到利用离散元模拟的岩石之上，研究了岩石中裂纹的扩展与失效过程及预先存在的岩石解理断层对于裂纹的影响。模拟发现：在起爆点附近，爆炸波的应力幅值超过岩石的屈服极限，导致岩石压缩粉碎；在更远区域，失效模式则转变为拉伸加载导致的径向裂纹扩展。当岩石中预先存在解理断层时，爆炸波将在断层上发生反射和耗散，显著影响岩石破碎特征；一个约45°沉降的断层可造成岩石中最严重的拉伸和压缩失效；但当断层的刚度和摩擦增加时，岩石失效区的尺寸将明显地减小。

3.2   国内研究的代表性成果

张振南等^{[55, 65]}针对岩石材料动态断裂的准确预测问题，发展了一种新的格子模型—离散虚内键模型。该模型由多边形/多面体元胞密堆积组成；元胞多边形/多面体的每个顶点均为一个格点；元胞内两两格点之间均设有弹簧，弹簧作用力函数由虚内键超弹性势函数导出；为使模型具有可调的泊松比，在元胞内两根相邻弹簧的夹角上采用广义的弹脆性Stillinger-Weber势函数设定变形抗力；为了统一描述岩石材料的拉伸强度、断裂应变、断裂表面能等特征，提出在弹簧压缩阶段及拉伸未达到强度极限时采用超弹性势函数，而在弹簧拉伸超过强度极限后改用双线性作用力。研究表明，通过一组实验对双线性拉伸断裂的分段函数进行参数校准后，模型可以精确地重现准静态和动力加载下岩石的裂纹形貌与载荷-位移函数。进一步加入JWL爆轰产物状态方程后，该模型可有效地模拟采矿爆破过程中冲击波、爆生气体、地应力等对岩石断裂特征的影响。

刘晓星等^[66–68]利用可描述粉体材料烧结致密化的离散元模型，研究了多孔电极材料的烧结动力学过程，及其达到部分烧结状态后的断裂行为和强度特征。在温度和压力的作用下，烧结致密化模型允许格点从细小离散的粉体状态开始相互吞并和长大，从而逐渐形成尺寸大但数量少、格点之间结合紧密的陶瓷块体。他们还模拟研究了粉体（初始格点）尺寸等因素对烧结动力学过程和烧结结果的影响；进一步针对烧结制备的脆性电极材料，改进了格点间大接触面积条件下的作用力模型，模拟研究了压缩和拉伸加载下的应力-应变关系、损伤成核、裂纹扩展、失效模式与强度阈值等。

吕文银^[69]对氧化铝陶瓷、玻璃等脆性材料在中/高应变率加载下的冲击破坏过程进行了模拟研究。在对陶瓷的建模中，将晶界设定为弱界面，晶界上的断裂能根据两晶粒的晶向角之差变化，对裂纹的扩展产生影响。模拟带孔隙的氧化铝陶瓷，所获得的压缩强度与实验结果基本吻合。对比发现：多孔陶瓷表现出一定的应变率敏感性，而致密陶瓷则未表现出率敏感性；且多孔陶瓷的裂纹密度时程曲线具有明显的“台阶”，而致密陶瓷中没有。进一步的分析表明，这些“台阶”对应多孔陶瓷中两种止裂机制—应力松弛和晶界转向。将数值模拟与高速摄像相结合，研究了圆柱玻璃杆撞击刚性壁过程，发现其破坏模式包括撞击端的破坏波传播和自由端的层裂。

吴建奎^[70]研究了脆性晶体材料在高应力率拉伸加载下的裂纹孕育起裂、高速扩展与振荡分叉过程。格子模型采用正三角形网格，预设一段半无限裂纹；初始时刻对裂纹两侧的格点分别施加V与–V的初始速度，从而在裂纹面上产生高应力率加载。模拟结果表明，高速裂纹的扩展呈现出单裂纹传播、微分叉、宏观分叉的演化过程。一旦达到起裂的临界条件，裂纹速度会瞬间起跳达到0.5倍瑞利波速，且随后的速度增长和扩展路径具有明显的相关性：单裂纹阶段对应裂纹速度突然起跳后的缓慢上升过程，微分叉阶段对应裂纹速度的微小振荡，而宏观分叉会造成裂纹速度的剧烈振荡。

王文强等^[71]利用正三角形网格的格子模型探讨了橡胶薄膜中裂纹扩展路径正弦形振荡的奇特现象。通常认为，裂纹路径的不稳定性只有在裂尖速度超过某个较高临界值后才会出现，但是在橡胶薄膜中低速裂纹即会导致正弦形振荡。王文强等分析认为，超弹性、黏弹性、非局域弹性是决定橡胶材料动态断裂行为的3个关键因素；如图3(a)所示，其模型中考虑了格点与其最近邻和次近邻的相互作用，从而表现出弹性变形的非局域效应；再通过弹簧作用力函数的设定，表现出超弹性和黏弹性。模拟结果表明：非局域效应尤为重要；在不考虑非局域效应（不计算次近邻格点间的相互作用）时，如图3(b)所示，裂纹在低速下只能直线扩展，而加入非局域效应后，如图3(c)所示，低速裂纹发生正弦形振荡。此外，当黏弹性过小时，裂纹也沿直线传播；超弹性对于表现橡胶的大变形很关键。通过调节超弹性、黏弹性和非局域弹性的参数，裂纹的路径、传播速度以及正弦形振荡的波长和振幅均可以调控，大多数实验结果可以被定量地重现。

图  3  橡胶薄膜裂纹扩展的格子模型：（a）格点间的最近邻与次近邻相互作用，（b）不考虑次近邻作用时裂纹直线传播，（c）考虑次近邻作用以表现非局域效应后裂纹扩展路径出现正弦形振荡

Figure  3.  Lattice model for the crack propagation in a rubber film: (a) the interaction among the nearest lattices and the next nearest neighbors; (b) the crack propagates linearly when the interaction with the next nearest neighbor were ignored; (c) the crack propagates oscillatorily when the nonlocal effect contributed by the next nearest neighbors was modeled

下载: 全尺寸图片 幻灯片

傅华^[72]利用格子模型与有限元相结合的方法，研究了高聚物黏结炸药在动态加载下的压剪变形、热点生成、反应点火等过程。他对炸药晶体颗粒采用有限元建模，以表现晶粒具有的黏弹塑性响应；对黏结剂采用格子模型建模，以表现黏结剂在拉伸加载下的脆性；炸药晶体边界与黏结剂边界的接触面采用有限元与格子模型相衔接的过渡方法。如图4（a）所示，衔接的关键是在两模型的边界上加入一层过渡微元，它们既属于有限元网格，也属于格子模型的正三角形网格，将同时参与两个模型的计算，并实现两模型间边界条件的传递。图4（b）展示了炸药模型中有限元晶粒和格子模型黏结剂的显微形貌。如图4（c）所示，模拟表明，应力波扫过后压剪摩擦形成的热点区域多集中在晶粒之间的区域，多边形晶粒间的相互作用是热点生成的重要原因。

图  4  格子模型和有限元网格结合示意图（a），炸药模型图（b）（其中蓝色基体为黏结剂，红色颗粒为炸药晶体），应力波扫过后黏结剂与炸药晶粒的摩擦升温（c）

Figure  4.  (a) Schematic of a model combined by lattice model and finite element method; (b) the model of polymer-bonded explosives (Blue matrix represent binder, red particles represent HMX crystals); (c) the temperature rise induced by the friction between explosive particles and binders under dynamic loading

下载: 全尺寸图片 幻灯片

于继东等^[73]建立了高聚物黏结炸药挤压点火与裂缝燃烧的力-热-化学耦合计算模型，可以同时描述裂纹萌生扩展、产物气体扩散和断面点火燃烧等复杂过程。他们利用格子模型描述炸药颗粒与黏结剂的损伤断裂，利用气相离散元模型描述格点热分解之后气体产物的膨胀效应，利用WSB两步法化学反应模型计算固体格点向气体格点转变时的反应速率和温度。模拟获得的典型结果如图5所示：初始样品中存在与加载方向呈45°角的微裂纹。当冲击波从下至上扫过样品后，微裂纹的两个尖端上萌生出两条翅裂纹；微裂纹上下两平面出现严重的摩擦滑移，其尖端还扩展出剪切型裂纹，均导致显著的温度上升，使得炸药晶体开始燃烧并生成气体。模拟结果表明，在冲击加载达到压力平衡之后，裂纹面的燃烧过程进一步驱动翅裂纹扩展，直到翅裂纹贯穿模型样品。

图  5  裂纹扩展、气体扩散和燃烧反应耦合的格子模型计算结果

Figure  5.  Simulations of crack extension, gas diffusion and combustion by lattice model

下载: 全尺寸图片 幻灯片

喻寅等^{[47, 54, 74]}利用正三角形网格的二维格子模型研究了多孔陶瓷、透明陶瓷、复合陶瓷等材料中的介观损伤演化机制与宏观冲击响应规律。其模型的主要特点是考虑了陶瓷材料中晶粒取向、晶界强度、孔洞裂纹、第二相颗粒等微结构特征（见图2）对宏观响应的影响。例如，对于多孔陶瓷，模拟了冲击波扫过后孔洞的变形和塌缩过程（图6（a）），且通过回收实验观察到与模拟结果相同的孔洞破坏特征（图6（b））；模拟发现了由于剪切裂纹网络扩展和破碎颗粒转动所导致的滑移和转动变形新机制（图6（c）），分析了脆性材料中孔洞、裂纹等微介观演化对于宏观冲击平衡态的影响（图6（d）），揭示了多孔脆性材料中宏观冲击波剖面传播与介观损伤演化之间的对应和关联，解释了冲击波逐渐演变为双波结构的原因（图6（e）～图6（f）），探讨了孔洞随机排布以及以三角形（图6（g））、四边形（图6（h））和六边形（图6（i））点阵规则排布对受冲击脆性材料的介观损伤演化和宏观冲击响应的影响。

图  6  多孔陶瓷冲击破坏的介观机制和宏观特征

Figure  6.  Mesoscopic deformation and macroscopic response of shocked porous ceramics

下载: 全尺寸图片 幻灯片

4.   格子模型的主要问题与改进方向

（1）非线性响应

格子模型在对弹脆性材料的建模与模拟研究中取得了巨大的成功，但仍缺乏对金属等延展性材料的表现能力；且在不同的应力、温度条件下，岩石、陶瓷、玻璃、固体炸药等材料也可能展现出不可忽略的黏性和塑性响应；在极端压缩条件下，材料模量随强度逐渐升高的现象也不适合用线弹性的格子模型进行描述。因此，需要系统地发展可精确描述材料非线性响应的格子模型。Ostoja-Starzewski等^[57–58]及Horie等^[48–49]已尝试为弹簧设定弹性–塑性硬化响应；Krajcinovic等^[60]使用非线性增长的排斥力描述材料的冲击响应；张振南等^[55]提出弹簧断裂的双线性响应。在他们的研究基础上，未来的弹簧作用力可能会吸收多种设定的优点，精细地表现出不同压强、应变率、温度等条件下的非线性压缩、“位错”滑移（弹簧的断开和愈合）、韧脆转变等行为。

（2）多物理耦合

格子模型目前以处理准静态和动态的损伤断裂演化问题为主，但在诸多与脆性断裂密切相关的科学问题中，材料的电、光、热、化学演化同样起着至关重要的作用。例如，铁电陶瓷在冲击波压缩下的放电失效现象是动态断裂与电学击穿相耦合的问题；透明陶瓷的冲击消光和冲击发光现象是动态断裂与裂纹表面光学散射、裂纹内部摩擦升温发光相耦合的问题；矿山开采与岩爆灾害等是动态断裂与爆生气体膨胀相耦合的问题。因此，需要实现格子模型的功能拓展及其与其他物理模型的有机结合。刘晓星等^[66–68]发展了格点的烧结动力学模型，描述了粉体通过烧结致密化转变为块体的过程；于继东等^[73]开发了力-热-化学耦合的炸药挤压点火与裂缝燃烧模型。未来的格子将会充分借鉴他们的多物理建模方案，实现对不同材料中力、热、光、电、化学演化的耦合模拟。

（3）跨尺度衔接

格子模型目前的研究尺度主要在微米级到厘米级之间。然而，在军事、工业等应用中通常需要回答微介观的损伤断裂如何影响宏观大尺度上的装置失效、建筑物破坏等问题。例如，钻地弹侵彻建筑物过程中，微介观尺度的炸药颗粒的挤压摩擦与宏观尺度的弹体、装药压剪变形均非常敏感地影响着弹药的安定性；在对其进行准确预测时，宏观整弹响应的模拟与炸药局部热点区域微介观演化的模拟有着同等重要的地位。因此，需要推动介观尺度的格子模型与更大和更小尺度上模拟计算的衔接。张振南等^[65]发展了网格自适应的格子模型，在岩石模型中远离裂纹的区域使用宏观大网格，在靠近裂纹尖端区域则不断加密形成介观小网格；傅华^[72]发展了格子模型与有限元相结合的算法，原理上讲其有限元网格可在非关键区域自适应地放大到宏观尺度，以节省计算开销。未来的格子模型将继承和创新，在跨尺度的动态断裂问题中追求物理真实与计算效率的兼顾。

5.   总结与展望

爆炸与冲击加载下的脆性材料动态断裂研究有重要的军事、工业需求。本文首先对比了有限元方法、流体动力学方法、物质点法、近场动力学方法、格子模型等多种脆性断裂模拟中常用的介观尺度网格方法与无网格/粒子方法的特点；进而讨论了格子模型的基本原理、相互作用设定、作用力参数确定、断裂判据选取、微结构/缺陷设置等建模过程；再分别介绍了国内外在使用格子模型开展脆性断裂研究中取得的一部分代表性成果；最后分析了格子模型在非线性响应、多物理耦合、跨尺度衔接等方面存在的不足以及改进方向。目前，格子模型已在岩石压缩破裂、陶瓷力学失效、炸药挤压点火等问题的研究中发挥了重要作用，将来持续改进与完善的格子必将在地质灾害预测、陶瓷可靠性提升、弹药安全性设计，以及其他更广泛的应用领域中发挥新的更大的作用。