采用由中国工程物理研究院流体物理研究所自主编写的格点-弹簧模型(Lattice-spring model)^[18-19]。该方法与物质点法、近场动力学以及格子模型等相似，主要解决大变形、裂纹、复杂界面等^[20-21]问题，现已在脆性介质断裂研究中发挥了重要作用。其算法原理：将连续介质离散分化为由格点/颗粒和弹簧相互连接的网络，通过计算格点与弹簧之间的演化过程和响应特征来获得连续介质整体系统的演化机制与规律。格点携带材料中物质微元的质量、位置、速度等信息，弹簧在两格点之间施加相互作用，并存储应变能和损伤断裂信息。图1显示了两颗粒之间的相互作用模型，邻近颗粒i和j之间分别设定沿法向和切向的弹簧。两颗粒中心连接的法向弹簧主要起着抵抗拉伸和压缩的作用，刚度系数为k_n；垂直两颗粒中心线且作用于颗粒接触边缘的切向弹簧主要起着抵抗剪切和转动的作用，刚度系数为k_τ。法向和切向作用力表示为

式中：$ {\delta }_{ij}^{n} $和$ {\delta }_{ij}^{\tau } $分别为颗粒i和j沿法向和切向的相对位移，$ {f}_{ij}^{n} $和$ {f}_{ij}^{\tau } $分别为作用于法向和切向上的力。

计算模型中的核心问题是准确地设定弹簧的刚度系数。Gusev^[22]提出了一种有限元-离散元参数映射方法，其基本思想是：先将目标材料的弹性常数转换为有限元模型中的相互作用参数，再进行映射，将有限元模型参数转换为弹簧刚度系数。Gusev验证了模拟不同周期性复合介质格点-弹簧模型的准确性。Yu等^[23]对比了通过格点-弹簧模型设定的孔洞模型的损伤演化结果与回收实验材料的微观形貌，进一步验证了格点-弹簧模型的准确性。

在本计算模型中，统一采用了Griffith能量平衡原理的断裂判据^[24]，即模型中相邻格点微元之间存在一根法向拉压弹簧和一根切向剪切弹簧。一对弹簧中储存的总应变能是法向弹簧中储存的拉伸应变能与切向弹簧中储存的剪切应变能之和；当弹簧的总应变能等于或大于形成新微裂纹表面所对应的表面能时，即发生断裂。一般认为，静水压压缩不会导致断裂，故模型中没有计入法向弹簧储存的压缩应变能对断裂的影响。具体而言，计算模型的断裂判据如下

式中：$ {U}_{ij} $为一对弹簧中储存的总应变能，$ {U}_{ij}^{n} $和$ {U}_{ij}^{\tau } $分别为法向和切向应变能，$ {c}_{0} $为两格点间微裂纹的长度，$ \gamma $为断裂表面能。当弹簧满足以上条件时发生断裂，两个格点之间原本存在的相互作用力（抗拉伸和抗剪切作用力）消失，在目标样品中形成一小段裂纹；而在靠近裂纹两侧的格点仍存在其他相互作用，其周围区域应力分布不均匀，当某一时刻或某一区域产生应力集中时，就有可能导致下一对弹簧的断裂，如此演化，弹簧断裂判据依次进行。

本研究所建立的基于格点-弹簧模型的冲击压缩构型如图2所示：左端红色部分代表铜（Cu）冲击加载平板，选择铜材料的原因在于相同的冲击条件下铜冲击平板所携带的能量高，获得的冲击响应信号更显著；冲击平板以恒定速度v_p沿x轴正方向冲击目标样品，如黄色箭头所示，使得冲击波在目标样品中自左向右传播；中间蓝色部分表示目标样品，即PMMA多孔材料，白色圆形部分表示孔洞；右端灰色部分为氟化锂（LiF）窗口，用于防止冲击波在目标样品后界面突然卸载而无法获得正常信号；模型上、下边界都设定为周期性边界条件。在本研究中，将x轴称为横向，y轴称为纵向。

基于准等熵加载技术研究背景，本研究设计的孔洞模型及其参数如下：整体模型大小设定为600 μm × 300 μm，Cu、PMMA多孔材料和LiF窗口材料的厚度（$\delta $）均为200 μm，单元格点半径为1 μm，计算单元格点数为413 816。设计孔洞排列模型的孔隙率（$ \varPhi$）均控制在10%，孔洞半径R设定为10 μm，整体模型所包含的计算单元格点数为396 735。材料的相关参数及断裂表面能（$\gamma $）^[25]如表1所示，其中：$\;\rho $为密度；E为杨氏模量；$\;\mu$为泊松比；PMMA的断裂表面能$\gamma $没有采用准静态下的实验参数，而是先基于密实PMMA的冲击加载实验结果反推高应变率条件下PMMA的断裂参数，然后将其应用于本研究中的PMMA多孔材料计算。计算中，孔洞直接作为微结构加入，无孔洞的地方仍为密实PMMA。根据低速下准等熵加载技术研究的需要，结合本程序的实际计算能力，平板冲击加载速度设定为v_p = 100 m/s。

如前所述，大部分研究人员都在关注多孔材料的孔隙形状、排列等微观结构对多孔材料冲击响应的综合作用^[8-14]，很少关注孔洞排列单一因素对多孔材料冲击响应的影响。为了研究孔洞排列对PMMA多孔材料冲击响应的影响，设计了随机排列（Random lattice）、四角点阵（Square lattice）、三角点阵（Triangular lattice）以及锥形的递减排列（Decreasing lattice）和递增排列（Increasing lattice），如图3所示。

图4展示了5种孔洞排列模型PMMA多孔材料在早期单次冲击加载时孔洞裂纹的扩展过程，为了观察孔洞的裂纹扩展演化，每种排列模型的记录时刻为50、100、130、180 ns。当冲击波传播50 ns时，靠近冲击端孔洞体积收缩，孔洞外壁开始破碎，纵向区域萌发裂纹。100 ns时，靠近冲击端孔洞已经开始被破碎颗粒填充，裂纹已经扩展并沿45°倾角方向延伸。130 ns时，各个模型的冲击波正好传播至PMMA多孔材料与LiF窗口界面位置，即冲击波已通过整个孔洞区域，说明孔洞排列方式不影响冲击波的传播速度；另外，也可以观察到冲击波传播完整个孔洞区域时，所有孔洞未发生完全塌缩变形，见130 ns时各个模型中孔洞的形貌。当冲击加载至180 ns时，孔洞塌缩程度开始呈现不同的现象。随机排列和三角点阵模型中，孔洞附近区域均已萌发裂纹，并产生破碎颗粒；在四角点阵模型中，虽然孔洞附近区域已经萌发裂纹，但有一列孔洞并未展开裂纹，孔洞外壁没有被完全压缩破坏，如图4（b）中180 ns时刻裂纹扩展的虚线框所示；同样，在递减、递增排列模型中部分孔洞外壁未被压缩损伤，没有产生破碎颗粒，与喻寅等^[26]的研究结论相似，规则排列中孔洞与孔洞之间的损伤区域可以产生应力屏蔽响应，使得一定范围内孔洞附近区域的损伤破坏受到抑制。由此可见，四角点阵、递增排列和递减排列方式的孔洞结构能够有效减缓PMMA多孔材料内部孔洞区域的应力集中和孔洞塌缩变形速率，其中四角点阵方式的孔洞结构影响效果最显著。

从图4中也可以观察到，每种排列模型的孔洞纵向附近区域优先出现裂纹，这是因为冲击波接触孔洞界面时反射稀疏波^[27]，稀疏波主要作用于孔洞横向附近区域，从而降低横向区域的压应力，导致孔洞纵向附近区域应力集中并产生裂纹。从理论公式上也可以分析出，孔洞附近区域出现应力集中现象，容易导致材料被压缩破坏。然而在冲击加载过程中应力集中并不是导致孔洞周围区域塌缩破坏的直接原因，通常还与材料的剪应力有关。二维平面内的最大分解剪应力定义为^[23]

式中：$ {\sigma }_{1} $和$ {\sigma }_{2} $为材料主应力。在理想的平面冲击状态下，目标样品处于一维应变状态，而且沿x轴和y轴的正应力$ {\sigma }_{x} $、$ {\sigma }_{y} $分别对应材料主应力$ {\sigma }_{1} $和$ {\sigma }_{2} $。最大分解剪应力$ {\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}} $的方向与x轴的夹角约为45°。可以观察到，孔洞附近区域应力偏离一维应变状态，主应力方向会偏离x轴和y轴，因此需要根据$ {\sigma }_{x} $、$ {\sigma }_{y} $计算出两个主应力^[28]，即

从图5(a)中的孔洞附近区域可以看到，远离孔洞区域的剪应力很小，且沿着远离孔洞区域扩展，孔洞附近区域容易产生剪应力集中，导致该区域的弹簧先达到断裂判据要求，剪切裂纹从孔洞纵向区域向远离孔洞内部扩展。由2.2节中各个模型的应力分布可知，部分孔洞周围区域应力达到300 MPa，而在图4对应的部分孔洞裂纹扩展图中并没有出现裂纹。这是因为单轴压缩下PMMA材料处于压剪复合状态，且静水压力较低而剪切应力分量较大，此刻PMMA材料发生断裂不是因为静水压力，而是因为剪切应力达到了阈值。在本计算的冲击加载条件下，PMMA多孔材料处于一维应变状态，该状态下静水压力高，而剪切应力分量小，导致部分孔洞周围区域应力虽然达到300 MPa，但其对应的剪切应力分量（即切向弹簧中的剪切应变能）并未达到断裂阈值，因此PMMA多孔材料中有部分孔洞并未出现裂纹。

随着冲击波的传播，裂纹逐渐向远离孔洞内部方向扩展，形成剪切裂纹，孔洞周围区域受到剪切破坏；且相邻孔洞之间的裂纹相互贯通，加快了孔洞的塌缩变形，破碎颗粒逐渐填满孔洞内部而达到材料密实化，如图5(b)所示。

图6记录了5种孔洞排列模型中孔洞附近区域早期被压缩损伤时沿冲击波方向（横向）的应力分布变化，即沿x轴的正应力$ {\sigma }_{x} $的分布。孔洞周围虽然存在应力集中和应力松弛区域，但是都在冲击波平均应力幅值上下浮动，冲击波应力很高，本研究将平均冲击应力幅值设为正值。在50和100 ns时刻，冲击波在PMMA多孔材料内传播，可以清楚地观察到当冲击波扫过孔洞时，每种模型中孔洞周围的应力分布不均匀，横、纵向区域应力变化显著，即孔洞横向区域应力低，纵向区域应力高，且应力集中于孔洞纵向区域。这是因为PMMA多孔材料在横向受到动态冲击加载时，冲击波沿着PMMA多孔材料横向传播，当冲击波传播到孔洞自由面时，会反射稀疏波，根据动量守恒定律，稀疏波起卸载作用，导致碰撞瞬间横向自由面区域应力迅速下降^[27]。这也解释了前面所述的裂纹萌发于孔洞纵向区域的现象。130 ns时刻，可以观察到四角点阵模型中同一层孔洞的横向区域是低应力区域且相互连通，降低了PMMA多孔材料与窗口界面处应力，影响了冲击波阵面的平整性。同样，在该时刻，三角点阵、递减排列、递增排列模型中，靠近PMMA多孔材料与窗口界面处的孔洞引起波阵面不平整，导致界面处应力不均匀且得到降低效果，影响了PMMA多孔材料与窗口界面处粒子速度的大小。待到冲击加载经历180 ns后，根据波阻抗匹配原理，当冲击波从低阻抗材料传播到高阻抗材料界面时，会反射比原冲击波速度更快的压缩波^[27]，LiF窗口材料的波阻抗远大于PMMA的波阻抗，因此界面处会反射压缩波，当压缩波在PMMA多孔材料中传播时，PMMA多孔材料右端孔洞周围应力迅速升高，进一步加快孔洞的压缩破坏。

图7显示了5种孔洞排列模型在目标样品与窗口界面处的宏观粒子速度剖面。从图7中可以观察到，每种排列方式的PMMA多孔材料被完全压缩破坏后，其粒子速度曲线均存在多个速度平台，依次标记为Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ，曲线整体呈阶梯式上升趋势。当冲击波第1次扫完孔洞并传播到PMMA多孔材料与窗口界面时，因LiF窗口的阻抗高于PMMA的阻抗，将反射压缩波，使得该界面处粒子速度上升，同时因孔洞塌缩变形降低粒子速度导致第Ⅰ平台的形成。同理，随着冲击时间的增加，对应的第Ⅱ、Ⅲ、Ⅳ、Ⅴ平台都是由多种波系与孔洞结构相互作用形成，直到第Ⅴ平台时，孔洞完全塌缩，孔洞内部被破碎颗粒填充密实化达到Hugoniot状态。虽然每种排列模型的粒子速度曲线趋势相同，但是对应每个阶段的平台速度大小却存在差异，且从第Ⅱ平台开始，差异更显著，尤其是随机排列模型，其对应的粒子速度曲线都是最低的。同时不难发现，四角点阵模型对应的粒子速度曲线的第Ⅰ、Ⅱ、Ⅲ平台左端有凸起，即每当多种波系作用于PMMA多孔材料与窗口界面时，粒子速度较快地升高而后降低形成平台。如前所述，这是因为规则的四角点阵排列结构左右两端孔洞受到损伤时产生的应力场区域可延缓中间孔洞损伤，进一步延缓孔洞的塌缩变形，导致多种波系作用于界面时延缓粒子速度的降低。

从PMMA多孔材料孔洞密实化到Hugoniot状态，可以看出，孔洞排布方式不同导致粒子速度各不相同，记录第Ⅴ平台的平均粒子速度并列于表2。从对应状态的应力云图中读取PMMA多孔材料与窗口界面附近区域的平均应力，也记录在表2中。从图7可以看到，引起第Ⅰ平台升起时刻相同，即冲击波扫完各个孔洞排列模型中的孔洞后到达PMMA多孔材料与窗口界面位置的时间相同，表明孔洞排布结构不影响冲击波传播速度，但影响粒子速度。同时，PMMA多孔材料的体积收缩，密度改变，也影响波阵面后压力。由表2中的数据可知，孔洞排列方式对PMMA多孔材料的粒子速度剖面和波阵面后压力的影响程度各不相同，其中孔洞的随机排列模型对降低粒子速度最有效，四角点阵排列模型对降低波阵面后压力贡献最大。因为当冲击波经过四角点阵孔洞区域后，每层4个孔洞均处于压缩状态，对于内侧的两列孔洞而言，左右两侧都受到其他孔洞的应力屏蔽影响，而对于外侧两列孔洞，左右两侧中有一侧受其他孔洞的应力屏蔽影响，另一侧承受密实材料区域的压缩，因此在早期孔洞损伤过程中，四角点阵中外侧两列孔洞附近区域出现相对严重的损伤裂纹，且随着时间的推移，外侧孔洞附近出现较严重的横向拉伸和剪切破坏而形成低应力区，对内侧两列孔洞起到明显的损伤屏蔽作用，从而达到降低波阵面后压力的作用^[26]。这说明通过孔洞排布设计，可以改变PMMA多孔材料的冲击响应行为。

采用格点-弹簧模型，分析了5种孔洞排列方式下PMMA多孔材料的介观裂纹扩展与非弹性变形、应力变化特征以及宏观粒子速度剖面曲线，得到以下结论。

（1）孔洞破坏形式以剪切断裂为主；四角点阵排列模型能有效减缓孔洞附近区域的应力集中和孔洞非弹性变形；在每种排列模型中均观察到了孔洞与孔洞之间的扩展裂纹相互贯通，促进材料的密实化。

（2）在各种排列模型中孔洞侧向（垂直于冲击方向）区域应力较低，纵向（沿冲击方向）区域应力较高，裂纹萌发于孔洞纵向区域；四角点阵模型中同一层孔洞的侧向区域是相互连通的低应力区，降低了界面处应力；四角点阵、三角点阵、递减排列、递增排列模型都显著影响冲击波阵面的平整性。

（3）孔洞排列方式不影响冲击波的传播速度（传播速度由总的孔隙率决定），但影响粒子速度剖面和波阵面后压力；孔洞的随机排列模型对降低粒子速度最有效，四角点阵排列模型对降低波阵面后压力贡献最大。