微极键基近场动力学的隐式求解方法及应用

周倍合， 朱其志， 李惟简

（1.河海大学岩石力学与堤坝工程教育部重点实验室，南京 210098；2.河海大学江苏省岩土工程技术工程研究中心，南京 210098）

固体材料破裂问题是固体力学研究的经典难题之一［1］，根据研究手段不同现有研究方法大致分为两类：物理试验和数值模拟. 相比于试验方法，数值模拟方法费用低、效率高，对计算模型的形状和尺寸没有特殊要求，边界条件与荷载的设置更加自由，突破了试验设备条件和试件制备方法的局限性. 目前，已出现多种数值方法，其中扩展有限元法（XFEM）和离散元法（DEM）是发展相对成熟并广受关注的方法. XFEM方法在有限元方法（FEM）的基础上引入额外的节点自由度和局部强化函数，可用于解决裂纹、孔洞、夹杂等间断问题，无须重新划分网格［3］，但不适用于模拟复杂交叉、分叉等多裂纹扩展贯通问题. 基于非连续性假设的DEM［6］方法可模拟裂纹的扩展贯通现象，但在反映材料连续区域的力学行为时存在明显缺陷.

Silling 等［8］提出基于非局部理论的PD 理论，使用空间积分方程重构固体力学运动方程，使方程在不连续处也有定义，实现在统一框架下自然地描述裂纹扩展过程. 早期的PD理论被称为“键基”模型，存在泊松比固定的缺陷［10］. 为了克服泊松比限制问题并使PD理论与经典连续介质力学的联系更加紧密，Silling等［9］提出了“态”的概念并发展了“态基”模型. 除态基PD模型外，学者们还提出多种改进的键基PD模型以解决泊松比限制问题，Han等［11］对这些模型进行了详细介绍. 其中，Zhu等［12］提出了考虑键切向变形的双参数键基PD模型（XPDM），在均匀变形条件下该模型可重现任意泊松比材料的变形，但在模拟非均匀变形问题存在不足. Diana等［13］提出的MPPD模型在XPDM模型的基础上增加了局部转动自由度，可考虑键的刚体转动对自身切向变形的影响，并使用该方法成功模拟了非均匀变形问题.

使用显式差分技术对运动方程进行时域积分是近场动力学的常规数值求解方法，对于静态和准静态问题，需要使用动态松弛技术对系统引入人工阻尼加快计算收敛速度［14］. 选择合适的人工阻尼和时间步长对收敛速度和计算精度具有很大影响. 为解决这一问题，本文从能量角度出发，依据最小势能原理推导了MPPD 模型物质点系统的控制方程，构造了系统的非局部刚度矩阵，实现用隐式方法求解静态和准静态问题. 使用该隐式方法对悬臂梁进行了弹性变形分析，并对双悬臂梁试件中的Ⅰ型裂纹和L-形板中的Ⅰ-Ⅱ复合型裂纹扩展过程进行了模拟.

1 微极键基近场动力学模型

根据近场动力学方法的非局部特点，材料区域R 中的某个物质点i 将与其近场域Hx（在二维和三维下分别为具有特征半径δ 的圆和球体）中的其他质点产生相互作用，这种相互作用通过“键”进行传递. 与只考虑轴向力的BB-PD模型不同，MPPD模型中考虑了键的切向变形和物质点的旋转，解除泊松比限制问题.可将MPPD模型中的键理想化为三个弹簧的组合，以键的初始状态为参考，三个弹簧分为：约束物质点沿着键的轴向和切向发生相对运动的弹簧以及约束物质点发生相对旋转的弹簧［13］，如图1所示.

图1 微极键基近场动力学模型示意图Fig.1 Diagrammatic sketch of MPPD model

根据牛顿第二定律，t时刻物质点i 的运动方程为：

其中：ρ 为物质密度；J 为单位体积该物质的转动惯量；ü和θ̈分别表示加速度和角加速度；f 和m 分别表示物质点i 和j 之间的相互作用力和相互作用力矩；b 和c 分别表示物质点i 处的体力密度和力矩密度；Vj为物质点j 的体积. 当式（3）中两个方程等号的左端项都为零时，物质点i 处于平衡状态.

在小变形条件下，根据式（2）可得键的轴向伸长率l 可表示为：

式中：n∥=ξ ||ξ =(cos α,sin α)是初始状态下键的单位方向向量；α 为初始状下键的方位角. 规定逆时针旋转方向为正，两个质点在t时刻的转角分别为θi和θj，则小变形情况下键的切向变形可表示为：

式中：n⊥=(-sin α,cos α)，取n⊥的方向为切向变形的正方向. 以n∥和n⊥为正交基，可以建立键的局部坐标系统，如图1（b）. 两质点的相对旋转角为：

对于均质各向同性线弹性材料，假设组成键的三个弹簧的劲度系数分别为kn、kt和kϑ，在t 时刻物质点i 所受的弹簧力和力矩分别为：

式中：向量外积n∥×n⊥表示力矩的正方向. 根据上述条件，物质点i 的应变能密度可表示为：

在均匀各向同性膨胀条件下γ=ϑ=0，此时MPPD模型的应变能密度只与l 有关. 将此时MPPD模型的应变能密度与经典连续介质力学的应变能密度比较，可得到键参数kn与材料宏观弹性参数之间的关系. 当kn已知，对于均匀的纯剪切变形情况ϑ=0，同样可通过比较应变能密度可得到kt的表达式［13］.

平面应力情况下：

平面应变情况下：

在（9）和（10）式中：E 为杨氏模量；ν 为泊松比；h 为模型的厚度. 引入物质点旋转自由度可提高对非均匀变形问题模拟的准确性，平面问题中kϑ的表达式为［15-16］：

2 隐式求解方法

2.1 系统控制方程

数值求解时需要将模型离散化，若离散后任意质点i 的近场域中有Mi个物质点，式（8）的离散形式为：

式中：ξij= ||xj-xi为初始状态时物质点i 和j 之间的距离；ϑij=θj-θi表示变形后物质点i 和j 的相对转角.υc(j)为物质点j 的体积修正系数，表达式为［17］：

其中，r=Δx 2，Δx 为物质点间距. 因为参数kn、kt和kϑ是基于具有完整近场范围物质点的应变能密度确定的，所以对于靠近自由表面或者材料的界面处的物质点需要进行键参数修正. 本文使用能量法［17］计算表面效应修正因子Gij.

假设离散后模型中共有N 个物质点，则物体的总势能Etot(X)可表示为：

式中：X=(x′1,x′2,…,x′N)T表示物质点系统的当前构型；x′i(i=1,2,…,N)为变形后质点i 的位置矢量；Vi表示质点i 的体积；bi和ui分别表示质点i 的体力密度矢量和位移矢量.

因为近场动力学中物质点之间存在的非局部作用与原子尺度有限元方法［18］中的多体相互作用类似，根据最小势能原理，当离散后的物质点系统总势能最小时物体处于平衡状态. 假定X=X∗时物质点系统的总势能最小，则

将式（16）代入式（15）得到整个系统的控制方程：

其中，位移u=X-X0，K 为系统的非局部刚度矩阵：

P 为系统的不平衡力，可表示为：

对于非线性系统，需要对式（17）进行迭代计算，每个迭代步中更新刚度矩阵和不平衡力，直到系统不平衡力P 等于零.

2.2 单元刚度矩阵和不平衡力

将物质点看作有限元节点时，式（17）与经典有限元方法的控制方程一致. 有限元方法常用的单元形式只能描述邻近节点之间的局部作用，不适用于考虑非局部作用PD模型. 为了描述这种非局部作用，可将中心质点i 与其近场域中所有质点组成的集合作为一个单元，如图2所示. 单元中，中心点与近场域中所有非中心点间存在相互作用，但非中心质点之间互不影响. 因此，在单元刚度矩阵和不平衡力中，只有与中心点有关联的位置上元素不为零，其余的元素均为零. 物质点i 处对应单元的刚度矩阵和不平衡力分别为：

图2 MPPD-FEM模型中的一个单元Fig.2 Illustration of the element used in MPPD-FEM model

其中，j=1,2,…,Mi，Mi表示物质点i 近场域中包含其他物质点的个数. 在MPPD模型中每个物质点具有三个自由度，物质点i 对应单元的刚度矩阵为3(Mi+1)阶的方阵.

2.3 损伤的描述

近场动力学中，材料损伤可通过消除物质点之间的相互作用来描述. 在微观弹脆性（PMB）模型中，物质点之间键的伸长率超过临界值lc时，键失效［19］. 引入标量特征函数λ(η,ξ)来描述键的失效情况：

由于MPPD模型解除了泊松比限制，不同泊松比时的临界伸长率可按以下方法计算［20］：

其中，Gc为断裂力学中的能量释放率；β=3δ 4π；β′=0.238 73δ；μ 和κ 分别为剪切模量和体积模量：

物质点的局部损伤定义为该物质点近场域中失效键的数量与初始状态时键总数的加权比值：

使用微观弹脆性（PMB）模型能够有效模拟岩石、混凝土等脆性材料的破坏问题. 为简化问题，本文基于键基模型的微观脆性材料进行研究，使用最大伸长率准则判断MPPD模型中键的轴向弹簧kn是否发生破坏，并假设轴向弹簧破坏时切向弹簧和旋转弹簧也一同破坏，对应物质点之间的相互作用失效.

2.4 隐式方法的实现流程

微观弹脆性材料的试件在出现损伤前为线性系统，荷载作用下发生线弹性变形. 随着荷载增加，物质点系统中有达到临界伸长率的键失效时，材料开始出现损伤. 继续加载过程，失效键的数量不断增加，材料中损伤累积并导致力学性能逐渐衰减. 因此，当计算过程中有新的键失效时，需要及时更新整体刚度矩阵，将失效键对应的贡献值从中减去. 为了准确模拟材料的损伤变化过程并且兼顾计算效率，参考文献［21］提出的C方案，每个迭代步中只断开达到临界伸长率且伸长率最大的Nb个键，然后更新刚度矩阵. 当某个迭代步中没有新的键失效时，增加荷载然后再进行修正刚度矩阵的迭代过程，直到荷载大小达到停止条件或试件完全破坏时结束程序. 隐式算法的程序实现过程如图3所示，本文在计算过程中取Nb=4 .

图3 隐式求解方法的流程图Fig.3 Flow chart of the implicit solution method

3 算例与分析

3.1 弹性变形分析

如图4 所示，平面应力状态的悬臂梁，跨长为L=1 m，高度H=0.3 m，弹性模量E=200 GPa，泊松比ν=0.2. 梁的左端完全固定，右端自由并受到竖直向下的均匀分布的荷载，荷载合力F=1 kN. 用均匀分布的物质点对模型进行离散，物质点间距Δx=0.01 m，近场域半径δ=3.015 Δx .

图5（a）为使用MPPD方法模拟得到的沿x 和y 方向位移分量分布云图，FEM 法所得结果如图5（b）所示.对比可知，两种不同方法的模拟结果十分接近，位移场分布规律相同，但在数值上存在细微差别. MPPD方法计算的沿x 和y 方向的位移分量最大值都略小于有限元的结果，但相对于有限元的误差不超过1.0% .

图4 悬臂梁示意图Fig.4 Geometry and boundary conditions of a cantilever beam

3.2 裂纹扩展模拟

3.2.1 双悬臂梁中Ⅰ型裂纹扩展模拟 平面应力状态双悬臂梁的几何信息与边界条件如图6所示，厚度取h=10 mm，弹性模量E=200 GPa，泊松比ν=0.2，能量释放率Gc=962.8 J/m2. 预留初始裂纹位于试件的水平中心线上，长度为100 mm. 用均匀分布的物质点对模型进行离散，物质点间距Δx=0.5 mm，近场域半径δ=3.015 Δx. 双悬臂梁试件中裂纹先沿着水平方向进行开展，当裂纹接近试件固定端时将会发生弯折，如图7所示. 裂纹沿直线扩展阶段的位移-荷载曲线见图8，其结果与Abaqus中基于虚拟裂纹闭合技术的有限元分析结果相吻合.

图5 悬臂梁模型的位移分量分布云图（单位：m）Fig.5 Node displacements obtained by MPPD and FEM

图6 双悬臂梁几何示意图Fig.6 Geometry and boundary conditions of a double cantilever beam

图7 双悬臂梁试件中的裂纹路径Fig.7 Crack paths of double cantilever beam

图8 双悬臂梁裂纹开展的位移-荷载曲线Fig.8 Displacement-load curves of double cantilever beam crack propagation

3.2.2 L-形板中裂纹扩展模拟 普通混凝土材料的L-形板试件几何信息与边界条件如图9（a）所示，厚度为h=100 mm，弹性模量为E=25.85 GPa，泊松比ν=0.18，能量释放率Gc=65 J/m2，按平面应力问题处理.模型用均匀分布的物质点离散，物质点间距Δx=2.5 mm. 为了准确模拟L-形板中混合型裂纹的开展情况，适当增大了近场域的半径，取δ=5.015 Δx. 图10展示了不同位移条件时的模拟结果，由图可知，裂纹从L-形板的角点处出现，初始阶段与水平线成一定角度（约为26°），随着加载位移增大，裂纹开展方向逐渐趋于水平并向试件的右边界扩展. 模拟结果与图9（b）中试验结果［22］对比具有较好的一致性.

图9 L-形板示意图和试验观测的裂纹形状Fig.9 Sketch map of L-shape structure test and the experimentally observed crack path

图10 L-形板中裂纹扩展模拟结果Fig.10 Crack propagation mode in L-shape structure test

4 结论

1）本文从能量角度出发，依据最小势能原理重新推导了MPPD模型物质点系统的控制方程，提出一种适用于MPPD模型的隐式求解方法，用于提高静态和准静态问题的求解效率. 提出一种可以描述物质点之间非局部作用的单元形式，不再以单独的近场动力学键为单元，可提高组装整体刚度矩阵的效率.

2）分别使用隐式求解格式的MPPD方法和有限元方法对非均匀变形的悬臂梁进行弹性响应分析，得到物质点沿x 和y 方向位移分量的分布云图. 不同方法模拟得到的位移场分布规律相同，在数值上存在微小差别，但满足计算误差要求，验证了MPPD模型隐式求解方法的有效性.

3）使用按隐式方法求解的MPPD模型对双悬臂梁试件中的裂纹扩展现象进行模拟，裂纹沿预期方向稳定开展，并且裂纹水平直线阶段的位移-荷载曲线与Abaqus中基于虚拟裂纹闭合技术的有限元分析结果吻合良好. 此外，还对普通混凝土材料的L-形板进行了模拟，所得混合型裂纹的开展路径与试验结果十分接近. 通过以上两种算例验证了MPPD 模型隐式求解方法用于模拟静态和准静态条件下裂纹扩展问题的可靠性.