键相关近场动力学及其与经典连续介质力学的形貌耦合断裂分析综述

一、引言

工程结构服役期间，金属、地质材料与混凝土等材料中极易发生裂纹的萌生、扩展与分叉，在数值模拟中准确捕捉裂纹的完整演化过程具有重要的工程意义。经典断裂力学以应力强度因子、能量释放率等少量参数表征裂纹场，建立了成熟的断裂判据体系，但通常需要预先已知裂纹几何与扩展路径，且依赖小范围屈服假设，当断裂过程区扩大、开裂模式高度复杂时，这种少参数描述便难以胜任。内聚力模型（CZM）通过在潜在裂纹路径上预置牵引—分离关系，能够以统一方式描述不同裂纹模式与界面脱粘，然而常规有限元实现中界面单元须沿预设路径预插入，迫使裂纹沿单元边界扩展，导致结果具有网格依赖性。相场方法用连续相场变量将尖锐裂纹正则化为有限宽度损伤带，无需显式追踪裂纹面即可自然捕捉复杂裂纹，但其能量泛函高度非线性且含内长度参数，断裂区需要足够精细的网格才能解析相场梯度，计算开销大且对长度参数敏感。

近场动力学（PD）是Silling于2000年提出的非局部固体力学理论，其内力以有限邻域内物质点间相互作用的积分形式表述，用非局部作用替代空间导数，因此即使位移场不连续时公式依然有效，能够自然表征裂纹。PD已广泛应用于冲击爆破、热—力耦合、功能材料电—化学—力学效应等复杂断裂问题的数值模拟，也可与相场、非局部损伤理论等方法结合形成互补能力。然而，适合工程应用的PD本构建模仍面临挑战。基于键的PD（BBPD）模型中键的本构响应仅取决于该键自身的变形，其等效泊松比为固定值（如二维平面应力条件下为1/3），难以表征典型工程材料；非普通态基PD（NOSBPD）模型通过定义与点关联的非局部变形梯度，可直接调用经典连续介质力学（CCM）本构关系，但普遍存在零能模式问题，会引起位移场与应变场的数值失稳，已有的稳定化方案往往引入额外的非物理稳定项或调节参数。

为此，Chen提出了键相关近场动力学（BAPD）模型：为每条键引入键相关变形梯度，使键的变形被直接评估而非由点关联梯度推断，从而抑制NOSBPD中常见的伪模式。后续研究系统考察了其稳定性与收敛性，导出了弱形式，证明键相关策略在避免零能模式的同时保留了直接采用经典本构的能力。但与NOSBPD对每个邻域只计算一个点关联变形梯度不同，BAPD须为每条键计算相应的键相关变形梯度，计算成本显著增加。提高计算效率的一条有效途径是PD–CCM耦合建模：仅在裂纹区域使用PD，其余区域由CCM（有限元）描述。已有工作包括重叠域方法、PD点与有限元的界面耦合算法、局部自适应细化与可变视域策略等。基于应变能密度等价，Lubineau等构造了空间变化的“Morphing（形貌）"函数，使BBPD与CCM在统一框架内实现材料属性的平滑过渡，该方法随后被推广到态基PD与CCM的耦合，并发展出可随裂纹演化自适应调整PD域的形貌耦合框架，且已在商业软件中实现。

韩非团队在“Computational Engineering, Finance, and Science"期刊发表了文章“An improved bond-associated peridynamic model and its adaptive coupling with CCM for fracture analysis"，该文将现有形貌类PD–CCM耦合框架推广到BAPD模型，用于准静态与动态断裂分析。在均匀小变形与线弹性假设下，基于BAPD与CCM在物质点处的应变能密度等价，重新表述了BAPD力态的修正因子；在此基础上建立了以“Morphing"函数实现两模型刚度平滑过渡的BAPD–CCM耦合模型：裂纹区域由BAPD描述裂纹的萌生与扩展，其键相关变形梯度允许直接使用CCM本构且不易出现零能模式，其余区域由CCM建模以降低计算成本。通过二维与三维的准静态及动态算例，从求解精度、断裂模拟能力与计算效率三方面对所提模型进行了系统评估。

二、CCM与PD基本理论

CCM理论在拉格朗日框架下由线动量守恒原理导出运动方程：物质点加速度由其一阶Piola–Kirchhoff应力张量的散度与体力共同决定，应力张量与变形梯度（单位张量与位移梯度之和）功共轭。PD理论则将运动方程表示为对邻域内力态的体积分形式：某点的加速度由其视域（horizon）δ邻域内所有键的力态差与体力决定，其中键向量定义为两物质点的参考位置差，力态的具体定义随PD模型而异。

在NOSBPD模型中，力态由影响函数、一阶Piola–Kirchhoff应力、点关联形状张量的逆与键参考位置向量构造，其中点关联变形梯度通过邻域内加权的当前构型与参考构型位置向量态的张量积积分获得。由于本质上是加权平均，点关联变形梯度部分忽略了变形的局部一致性、方向性与峰值特征：同一条键由两端点各自的点关联变形梯度预测的变形一般不一致，即不满足运动学约束。

图1  NOSBPD（a）与BAPD（b）模型变形梯度构型的对比

BAPD模型将变形梯度与单条键相关联（图1）：力态等于修正因子与键相关力态的乘积，键相关力态由影响函数、键相关应力、键相关形状张量逆与键参考位置向量构造。键相关变形梯度在两物质点邻域的交叠域上积分定义，形状张量同样定义于该交叠域上，键相关应力则由交叠域内应变能密度对键相关变形梯度的偏导给出。由于点关联变形梯度是邻域平均的结果，族内多条键共享同一变形梯度，单键的变形扰动在平均过程中被削弱甚至抹平，模型丧失分辨单键变形的能力并可能产生零能模式；而键相关变形梯度与每条键一一对应，即使某条键发生变形扰动也能在模型中得到反映，从而缓解乃至避免零能模式。

在小变形假设下对两模型线性化：CCM中一阶Piola–Kirchhoff应力可由Cauchy应力近似，后者由广义胡克定律通过弹性刚度张量与无穷小应变张量联系；BAPD中键相关应力与键相关无穷小应变的关系按CCM类比建立，键相关位移梯度由键相关变形梯度表达式在无穷小条件下退化和对导获得。当键相关邻域与点邻域等距同构时，交叠域关于键对称，键相关形状张量与变形梯度满足对键取反的对称性，运动学约束严格成立，相应的运动学约束BAPD模型避免了NOSBPD中两端正负键力态不一致的问题。

三、基于应变能密度等价的改进BAPD模型

针对线弹性材料承受均匀无穷小变形的情形，本文基于BAPD与CCM在物质点处应变能密度的等价重新表述BAPD力态中的修正因子。这种在均匀变形下建立模型等价性的处理沿用了Silling等关联NOSBPD与CCM时所引入的对应性定义，并不限制模型只能应用于均匀变形场。将键能量密度沿其邻域积分得到BAPD应变能密度；在邻域内位移梯度近似均匀的假设下，键相关位移梯度退化为均匀梯度，键相关应变与应力分别等于均匀应变与均匀应力，且转动部分对能量无贡献。

由能量等价条件可以推出：为保证任意均匀应变下BAPD应变能密度与CCM一致，键相关形状张量逆与含修正因子的键参考位置张量积的影响函数加权积分必须恰好等于四阶单位张量。传统修正因子定义为键相关族体积与邻域体积之比，仅是一种体积度量；本文在球形邻域与径向影响函数假设下，要求键相关形状张量与点关联形状张量沿任意键方向具有相同的二次型值，推导出新的解析修正因子χa，即键方向上点关联形状张量逆与键相关形状张量逆二次型之比，并在附录中严格证明了该修正因子使能量等价条件精确成立。

四、BAPD-CCM耦合模型

图2  BAPD-CCM模型中BAPD域、CCM域与过渡域的分布

BAPD-CCM模型通过应变能密度等价将两模型耦合，并采用“Morphing"策略实现过渡域内刚度贡献的平滑过渡。计算域被划分为三部分（图2）：CCM域Ω1由CCM控制，BAPD域Ω2由BAPD控制，二者之间的过渡域Ωm实现平滑衔接。控制方程包括运动学许可与兼容性条件（无穷小应变定义、相对位移、给定位移边界）、动力学许可条件（含应力散度与力态积分之和的运动方程、给定面力边界）以及本构方程（CCM应力由其刚度张量与应变联系，BAPD力态由修正因子、影响函数、刚度张量、键相关应变、键相关形状张量逆与键参考位置构造）。

为实现两模型刚度张量的平滑过渡，引入取值于[0,1]的“Morphing"函数α：BAPD刚度张量取α倍的全刚度张量，其余刚度贡献则分配给CCM。当α=0且其邻域内各点α均为0时，该点位于CCM域，总刚度由CCM提供；当α=1且邻域内均为1时，该点位于BAPD域，总刚度由BAPD提供，其应变能密度遵循BAPD表达式；当α介于0与1之间时，该点位于过渡域，总刚度由CCM与BAPD共同贡献。由过渡域应变能密度与全CCM应变能密度的等价条件，可解析导出CCM域刚度张量的表达式，其中BAPD刚度贡献以两端点α平均的方式进入积分。

损伤描述采用键临界伸长准则：脆性材料中一旦某键的伸长达到临界值，该键发生不可逆断裂，键的状态用标量变量表征。为将微观断键与宏观断裂联系起来，对每个物质点定义标量损伤变量（取值0～1），等于邻域内断键体积占比。在BAPD-CCM模型中，断键仅限于BAPD域内发生，域外损伤恒为零，保证了损伤演化的局部性。

五、数值算法

图3  BAPD-CCM模型的离散方式及刚度从CCM到BAPD的过渡

数值求解方面，CCM模型由有限元求解，BAPD模型采用无网格节点积分；时间离散对准静态与动态问题均采用显式格式，准静态问题采用自适应动力松弛（ADR）方法，动态问题采用中心差分法。此外引入强度诱导与断键诱导两种准则，在分析过程中自适应扩展BAPD域。计算模型离散为单元与节点（图3）：当节点位于CCM域时，保留单元连接关系，内力由有限元格式的CCM模型计算；当节点位于BAPD域时，解除该节点所有连接关系，单元节点转化为独立的PD粒子，内力由BAPD模型提供；当节点位于过渡域时，单元贡献与粒子贡献并存，两种模型各提供部分内力。由于过渡域内单元节点与PD粒子共位，耦合界面处的位移连续性得以自然保证。

六、模型验证

6.1 改进BAPD模型的评估

首先数值检验不同修正因子下BAPD与CCM的能量等价性，以等价条件中积分张量对单位张量的符合程度度量误差。由于修正因子为标量且影响函数仅依赖键长，积分项具有径向对称结构：在球形邻域上所有非对角项由对称性抵消，对角项相同，因此只需考察单一主值λ——理论上其值应为1。结果表明：在视域δ与网格尺寸之比取3、4、5时，传统体积比修正因子χv对应的主值分别为1.090、1.049、1.041，均存在正偏差；而新修正因子χa在所有视域下均给出λ=1.000，与附录中的理论证明一致。值得注意的是，由于χv常使λ大于1，耦合公式中的张量项在α趋于1时不能退化为单位张量，可能导致BAPD域内出现负的CCM刚度张量，这也从侧面说明χa更适合耦合模型。

图4  矩形板的几何与边界条件

随后考察含中心水平预裂纹（长约10 mm）的40 mm×100 mm二维矩形板（图4）：下边界固定，上边界施加5 MPa均匀牵引力，网格尺寸0.5 mm，视域取3.03倍网格尺寸，材料为E=6 GPa、ν=0.25的各向同性线弹性体。以FEM解为参考，在x=15 mm处选取观测线，比较两种修正因子下BAPD模型的位移剖面（图5）。

图5  观测线上位移uy与ux的对比

对于水平位移ux，三种解在观测范围内总体接近，但放大视图显示χa的结果比χv更接近FEM参考解；对于竖向位移uy，两BAPD解在裂纹前方与FEM解较为接近，穿过预裂纹后差异明显增大，χa结果始终更贴近参考解，而χv偏差较大。为量化整体偏差，评估相对L2范数、相对H1半范数与相对能量范数三种误差指标：χa对应误差分别为1.5574%、1.8124%与1.8135%，而χv对应误差分别高达5.7974%、6.2821%与6.2494%，表明在该算例中χa产生的整体位移与能量误差显著小于χv。

6.2 准静态载荷下的BAPD-CCM模型

图6  几何与边界条件（a）、CCM/过渡/BAPD域的划分（b）及“Morphing"函数α的分布（c）

采用承受准静态单轴拉伸的100 mm方板评估BAPD-CCM模型的精度（图6）：上下边缘施加uy=±1 mm的竖向位移，两边缘中点处固定水平位移；网格尺寸0.5 mm，共40,000个四边形单元与40,401个节点；视域取3.03倍网格尺寸并采用常数影响函数；材料参数为E=72 GPa、ν=1/3。板中央嵌入“兔形"BAPD域，过渡域宽6倍网格尺寸，图6（c）给出了“Morphing"函数α的分布。

图7  准静态加载下的位移云图：FEM（a）、BAPD-CCM（b）与BAPD（c）的ux；FEM（d）、BAPD-CCM（e）与BAPD（f）的uy

图7比较了FEM、BAPD-CCM与纯BAPD三种解的位移云图：三者的位移场均大体均匀，但纯BAPD解在板角附近出现明显弯曲的等值线（尤其图7（c）），表现出边界效应，而BAPD-CCM模型则不存在此类问题。

图8  准静态加载下BAPD与BAPD-CCM模型的对比：相对误差范数（a）与每载荷步计算时间（b）

以FEM解为参考量化整体误差（图8（a））：纯BAPD的相对L2误差、H1半范数误差与能量范数误差分别为0.108%、1.786%与1.882%，而BAPD-CCM相应误差仅为0.049%、0.026%与0.016%，明显更小。效率方面（图8（b）），单载荷步BAPD耗时38.307 s，BAPD-CCM耗时21.857 s，提速约43%。

6.3 动态载荷下的BAPD-CCM模型

动态算例沿用上一小节的几何与BAPD域设置：材料参数为E=72 GPa、ν=1/3、ρ=2450 kg/m³；下边缘竖向固定并在其中点约束水平位移以抑制刚体运动，上边缘施加1 GPa牵引力；网格尺寸1 mm（10,000个单元、10,201个节点），时间步长5×10⁻⁹ s，共计算50,000步（总时长0.25 ms）。

图9  不同时刻FEM（左）与BAPD-CCM（右）竖向位移uy云图的对比

为清晰观察位移波的传播，图9给出了前20 µs内FEM与BAPD-CCM的uy云图对比：两种位移场在空间分布与幅值上均高度吻合，但BAPD-CCM模型在耦合界面处存在轻微偏差。

图10  BAPD-CCM与FEM的动态响应对比：(a)~(c)三个代表性节点处的竖向位移uy时程；(d)应变能与动能时程

进一步监测三个代表性节点处的竖向位移时程（图10（a）~（c）），BAPD-CCM与FEM位移曲线的最大偏差低于3%；对所有时间步比较应变能与动能（图10（d）），二者的变化趋势与量值均相似，验证了耦合模型在动力学问题中的求解精度。

七、数值算例

7.1 准静态加载下的I型断裂

图11  I型断裂的几何与边界条件（a）及BAPD-CCM模型中的域划分（b）

考虑边长3.2 m的方板，上边缘中部含长1.625 m的中心预裂纹（图11）：材料为各向同性线弹性，E=28.3 GPa、ν=0.3，断裂能Gc=196 J/m²；视域取3.03倍网格尺寸，全域离散为26,002个节点与25,600个单元（网格尺寸20 mm）；采用断键诱导准则，在裂纹周围布置BAPD域，位移加载分100个增量步施加。

图12  不同增量步下α值（a）、单元与粒子（b）及损伤云图（c）的演化

图12展示了观测域内若干载荷步下α值、单元/PD粒子拓扑与损伤域的演化。第35步时几乎没有损伤发生；第36步时裂纹附近若干键开始断裂；随着断键累积，α>0的区域逐渐向外扩展，该域内单元节点迅速转化为PD粒子，损伤域沿裂纹路径快速向下生长；此后直到第99步，损伤前沿仅略有推进，表明裂纹一旦萌生，主裂纹扩展即在短时间内完成，与脆性断裂的特征一致。

图13  BAPD-CCM模型（a）与BAPD模型（b）的损伤云图

图13比较了BAPD-CCM与纯BAPD模型的最终损伤云图：两者裂纹路径总体一致，但损伤带宽度不同。可能的原因是，断键诱导算法中BAPD域仅在新断键附近激活并随裂纹生长渐进扩展，将损伤演化限制在裂纹路径附近的狭窄区域内，从而得到更纤细的裂纹形态。

图14  力—CMOD曲线与计算时间的对比

力—CMOD曲线（图14（a））表明BAPD与BAPD-CCM均能再现实验的总体趋势：峰值载荷对应大规模断裂的起始，其后载荷显著下降。计算时间对比（图14（b））表明，BAPD-CCM显著降低了成本，总模拟时间减少72.144%。

7.2 脆性材料中的动态裂纹分叉

图15  裂纹分叉算例的几何与边界条件

本小节重访经典的动态裂纹分叉基准算例：试件为100 mm×40 mm的狭长矩形板，中线含50 mm预裂纹（图15）；全域离散为64,000个均匀四边形单元（网格尺寸0.25 mm），视域取4.04倍网格尺寸；材料为钠钙玻璃，假设为脆性、均匀、各向同性线弹性，ρ=2.44×10³ kg/m³、E=72 GPa、ν=1/3，断裂能Gc=135 J/m²；上下边缘施加12 MPa均匀法向牵引力并保持恒定，总加载时长30 µs，由30,000个时间步求解；采用强度诱导准则激活BAPD域。

图16  6.5 µs时刻的最大主应力云图（a），以及若干时刻的α值分布（b）与损伤云图（c）

图16给出了该问题的动态断裂演化过程：t=6.4 µs时，上下边界产生的应力波在板中部相遇并在预裂纹产生显著应力集中；此刻部分物质点满足最大主应力准则而失效，裂纹附近自适应诱导出BAPD域，损伤开始形核。随后BAPD域跟随前进的裂尖向右传播，域内损伤持续生长；从t=14.5 µs到t=20 µs，损伤开始呈现分叉形态，BAPD域也随之向上下扩展。可能的机理是，裂尖附近反射应力波与后续入射波的相互作用形成两个对称的偏轴应力集中域，进而促进裂纹分叉。

图17  本文BAPD-CCM模型预测的裂尖速度与实验数据及已有数值结果的对比

裂尖速度的评估结果（图17）表明，模型预测值与实验数据及已有数值模拟结果均吻合良好。

7.3 落锤冲击下的混凝土板

图18  混凝土落锤冲击试验示意

最后研究落锤冲击下的方形混凝土板（图18）。试验装置简化为图19（a）所示构型，外载荷直接取自试验实测数据（图19（b）），加载半径在前5 ms内由0线性增至5 mm。网格尺寸1 mm，共80,104个六面体单元与91,926个节点，视域取3.03倍网格尺寸；材料密度与杨氏模量分别取ρ=2.5×10³ kg/m³与E=30 GPa，泊松比0.2，临界伸长取1.0×10⁻⁴；总计算时长37.4 ms，共500,000个时间步；采用强度诱导准则诱导BAPD域。

图19  落锤冲击载荷的简化示意（a）及相应试验力—时间历程（b）

为便于观察，图20与图21分别给出若干时刻α>0点与损伤点（ϕ>0）的分布，其中图20为俯视图，图21为仰视图。冲击中心周围的域最先被激活并切换为BAPD域，随后随损伤演化逐渐向外扩展；损伤域在BAPD域内萌生并沿板面传播，逐渐发展为更宽的损伤带，最终形成四条主导裂纹。对比上下表面可见：背面（下表面）总体上α>0域出现更早、范围更大，损伤也更显著；正面（上表面）则主要表现为冲击点附近的局部化损伤。

图20  俯视图：14.9 ms（a）、18.7 ms（b）、22.4 ms（c）、26.2 ms（d）时刻α>0的点，以及14.9 ms（e）、18.7 ms（f）、22.4 ms（g）、26.2 ms（h）时刻损伤ϕ>0的点

图21  仰视图：14.9 ms（a）、18.7 ms（b）、22.4 ms（c）、26.2 ms（d）时刻α>0的点，以及14.9 ms（e）、18.7 ms（f）、22.4 ms（g）、26.2 ms（h）时刻损伤ϕ>0的点

图22将数值模拟得到的裂纹路径与试验结果进行对比，二者吻合，进一步证明了BAPD-CCM模型捕捉动态断裂行为的有效性。

图22  BAPD-CCM模拟（a）与试验（b）所得裂纹路径的对比

八、结论

本文基于应变能密度等价改进了BAPD模型，并进一步建立了其与CCM的自适应耦合模型。主要结论如下：（1）基于BAPD与CCM在物质点处应变能密度等价推导的新修正因子χa，在多种视域下均精确满足能量等价条件，相比传统体积比修正因子χv显著降低了位移与能量误差（L2误差由5.80%降至1.56%）；（2）以“Morphing"函数实现刚度平滑过渡的BAPD-CCM耦合模型，在准静态与动态算例中均与FEM参考解高度吻合，误差范数远小于纯BAPD模型，且单步计算提速约43%；（3）准静态I型断裂算例中，耦合模型在再现力—CMOD曲线与实验裂纹路径的同时，将总计算时间减少72%以上；（4）动态裂纹分叉算例准确捕捉了分叉形态，裂尖速度与实验数据吻合；（5）三维混凝土板落锤冲击算例准确预测了四条主导裂纹的形成与上下表面损伤差异，裂纹路径与试验结果吻合。所提模型兼顾了BAPD描述裂纹萌生与扩展的优势和CCM的高效性，为工程结构的准静态与动态断裂分析提供了一条精度与效率兼备的数值途径。