基于非局部微分算子的近场动力学及固体材料应力数值模拟

李树忱，马鹏飞，王修伟，刘祥坤

(山东大学岩土与结构工程研究中心，山东，济南 250061)

固体材料的力学行为一直是工程结构领域关注的热点，探究其失稳破坏机理对工程防灾减灾具有重要意义。多年来，国内外学者针对此问题已取得了丰硕研究成果[1-4]，采用理论方法探究材料破坏机制时通常需要一定简化，无法很好地反映真实的情况，利用试验手段存在成本较高、时间周期长等困难。随着计算机技术发展，数值模拟手段越来越多的被用来研究材料的力学特性。

有限元法[5-7]、离散元法[8-9]、无网格方法[10]等数值方法的相继出现为探究材料力学行为提供了成本低且较快捷的方式，在模拟裂纹扩展时微分控制方程会遇到尖端奇异性等难题。扩展有限元理论运用预先设定裂纹扩展方向及长度来模拟材料破坏现象[11-12]。分子动力学方法可较好的体现结构微观信息，但是在模拟裂纹扩展等尺度时会遇到计算效率低、耗费时间长等不足[13]。

近场动力学方法是采用非局部作用特性来描述材料力学行为的理论[14]。该方法核心为求解空间积分方程，因此可较好的避免尖端奇异性等难题。近场动力学理论最早由SILLING 等[15]提出，黄丹等[16-17]首次将该理论引入国内，并针对混凝土等固体材料变形破坏开展研究。ZHOU 等[18]和谷新保等[19]利用近场动力学方法模拟了岩石类材料的裂纹扩展过程。WANG 等[20-21]建立了共轭键基近场动力学模型以突破泊松比的限制。ZHU等[22]提出了具有键旋转效应的近场动力学公式并对脆性材料变形行了模拟。王超等[23]在模拟潜艇破冰上浮的现象中运用了近场动力学方法。黄小华等[24]提出了双参数模型并对冲击荷载对破坏的影响进行了分析。HAN 等[25]利用连续介质力学理论与近场动力学耦合的方式对材料变形进行模拟。

近场动力学中以键为基础的微弹性模型具有广泛的应用[26]，微弹性模型控制方程为积分方程，在计算时不涉及应力及应变，但经典固体材料力学特性分析时应力与应变更为直观，并且大多破坏准则为基于应力的准则，在近场动力学理论中无法直接应用，文献[27]中提出利用力通量的几何等价关系计算应力，但其计算过程复杂且计算量相对较大。

本文结合前人研究基础，在近场动力学中引入非局部微分算子求解理论，建立微弹性应力分析模型。为验证有效性，利用扩展后的方法对几种二维平面问题进行模拟，同时将计算值与经典理论解对比，并且对粒子离散间距、泰勒项数及权函数的数值收敛性进行研究。结果表明：本文提出的方法可较好的计算结构处的应力值并进行相应的分析。

1 基本理论

1.1 近场动力学理论

近场动力学理论与传统方法相比最大的不同为非局部作用。如图1 所示非局部作用是指在一定范围内存在相互作用力，x′为x作 用域 δ内具有物理意义的物质点，在初始参考坐标系中两点的相对位置为x′-x=ξ ，并且x在t时刻的控制方程为：

式中：u(x,t)与u(x′,t) 分 别为x与x′变形时的位移，u(x′,t)-u(x,t) =η为变形后的相对位移，因此，此时x与x′相对位置可表示 ξ+η ； ρ(x) 与u¨(x,t)分别为密度与加速度；f为x与x′之间相互作用力；b(x,t)为t时刻体积力密度。

1.2 微弹性模型

在微弹性模型中相互作用力为相对位置 ξ与位移 η的函数，同时x与x′之间f大小相等方向相反，因此，式(1)可重新表示为：

式中，1/2 为每个相互作用键中的能量一半属于物质点x，其中的微势能ω可写为：

2 非局部微分算子

2.1 积分方程

应力-应变在经典Navier 平衡方程中为位移的微分函数，而非局部微分算子试图寻找其积分表达以适用于近场动力学模型，本节以二维为例推导微分算子的形式。如图1 所示，x与x′相互用点并且相对位置x′-x=ξ ，可考虑在x处进行Taylor级数展开：

2.2 应力-应变表达

由近场动力学与经典Navier 平衡方程对应关系[28]，这里将近场动力学中某点应力分量表示为：

本节对完整平板、含圆孔平板及含裂隙平板在外荷载作用下的变形进行模拟，并将应力分布结果与经典理论解对比以验证本文方法的有效性。现在考虑一个各向同性的平板在两端受到均匀荷载其尺寸及荷载分布如图2 所示。

平板尺寸1 m×0.5 m，采用平面应力条件计算，弹性模量E=18 GPa ， 泊松比ν=1/3 ，密度ρ=2730 kg/m3，两侧荷载为p0=15 MPa，以中心为坐标原点应力分布的理论解为：

3 数值验证

3.1 平板拉伸

计算时将模型离散为 200×100节点，间距取Δx=0.005 m，相互作用半径 δ=3.015Δx同时采用稳态求解方案，计算结果对比如图3 所示。

由图3 可知利用本文方法所计算的应力结果与理论解相比在模型中间区域吻合良好，但在边界及角落处存在较大误差，主要原因是边界附近的节点邻域不完整导致模量值失真从而求解位移场出现偏差，因此，基于位移求解的应力场在边界区域同样会出现误差。

3.2 含圆孔平板

为了进一步验证本文方法的有效性，3.2 节对具有高次精确解的含圆孔平板在荷载下变形进行模拟并将计算结果与经典理论解对比。现在考虑一个中间含有孔的方板在两端受到均匀荷载，尺寸及荷载分布如图4 所示。

平板尺寸0.5 m×0.5 m 并且中间含R=0.025 m的圆孔，采用平面应力计算，弹性模量E=15 GPa，泊松比ν=1/3 ，密度 ρ=2650 kg/m3，两侧荷载为p0=10 MPa，以中心为极坐标原点应力分布的经典理论解为：

由图5 可知平板的孔口区域表现出明显的应力集中现象，孔口应力远大于距孔口较远处的应力。同时最大和最小的应力均出现在孔边上，由开孔所引发的引力扰动也主要发生在1.5 倍圆孔直径的范围之内，从定性的方面验证本文方法的有效性。由直角坐标系与极坐标系的对应关系，可分别选取竖直方向、水平方向与环绕孔口区域的计算值与理论值比较，从定量的角度验证有效性，对比结果如图6 所示。

如图6 所示 σxx和 σyy沿竖直与水平方向对比， σxy沿半径 1.4R的环向进行对比。利用本文方法所得到的结果与理论解的整体偏差很小，误差仅出现在靠近圆孔边缘附近。理论解圆孔附近应力集中系数为3，本文方法所预测的应力集中系数为2.7 左右，逐渐远离孔口区域的计算结果收敛于理论解，靠近模型边界时会发生略微的失真但整体吻合较好，因此，可以验证本文方法的有效性。

3.3 含裂隙平板

实际结构破坏时往往会伴随裂纹的萌生与扩展，探究裂纹周围应力场分布对裂纹的发展具有指导意义，因此，本节利用该方法计算裂纹尖端应力场并将结果与理论解对比验证方法的有效性。

现在考虑一个中间含有裂纹的方板在四周受到均匀荷载，尺寸及荷载分布如图7 所示。平板尺寸0.5 m×0.5 m 并且中间含L=0.1 m的裂纹，采用平面应力进行计算，弹性模量E=15 GPa，泊松比v=1/3 ， 密度 ρ=2600 kg/m3，四周施加的荷载为p0=1 MPa，该问题的经典断裂力学解为：

式中：a为预制裂隙一半的长度；KI为应力强度因子，计算得KI=0.396 MPa·m1/2； θ与r为裂纹尖端附近点的角度与半径。计算时将模型离散为150×150 节 点，相互作用半径 δ=3.015Δx，其余计算参数与3.2 节相同应力计算结果如图8 所示。

由图8 可知在应力作用下平板裂隙处存在明显的应力集中现象，而且随着与裂隙距离的增加逐渐收敛于边界值，从定性方面验证方法有效性。

由图9 为理论解的局部对比，本文方法计算的结果与理论解相比在裂纹尖端处应力分布与最大值基本一致，而远离尖端处时计算解则会与理论解产生偏差并逐渐收敛于边界条件值。偏差产生的原因是经典断裂力学解的适用范围为裂纹尖端处，并且有模型无限大的前提假定，模型尺寸相比裂纹不够无限大导致理论解与计算结果不同，并且由于在力学分析中往往关注裂纹尖端处的应力状况，因此，可以验证本文方法的有效性。

3.4 裂纹扩展过程

前面几节所计算应力均为稳定状态下的结果，为了验证本文方法同样适用于动态裂纹扩展，本节将探究含预制裂纹材料在压缩破坏下的应力变化过程。计算模型在图7 的基础上进行修改，考虑一个中间含有裂纹板在竖向受持续位移边界条件，平板尺寸0.5 m×1.0 m，裂隙长度 0.28 m，裂隙与竖直方向的为夹角 30°，弹性模量取E=4.71 GPa，密度取 ρ=2600 kg/m3，临界伸长率取s0=0.0318，在持续位移加载下发生逐渐破坏，损伤及最大主应力分布过程如图10 所示。

由图10 可知，预制裂隙在竖向持续加载的条件下发生渐进破坏，裂纹扩展的方向与加载方向相同，其中最大主应力出现正负交替分布的情况(拉为正压为负)。在损伤初期阶段拉应力出现在裂纹端部靠近加载端的部位，随着荷载的不断增加拉应力不断变大，裂纹出现萌生扩展，拉应力持续出现在裂纹尖端处促使裂纹向最大主应力方向扩展，这过程与经典分析一致，可验证方法有效性。

4 数值收敛性分析

4.1 离散间距影响

本节开展非局部微分算子离散间距、泰勒项数及权函数收敛性研究，探究不同因素对结果的影响。以3.2 节含圆孔平板为例，将模型分别离散为 50×50 、 80×80 及 100×100节点，其余计算参数与3.2 节相同，x=0.0方向结果对比如图11 所示。

由图11 可知，随着离散节点的增加粒子离散间距随之减小，节点更加密集，在x=0.0方向上σxx与 σyy的最终计算结果更趋近于理论解，这一结论同样适用于y=0.0方向，在此不再赘述。由于增加节点数量会成倍的提高计算效率，因此，在实际计算时应该合理选择离散间距兼顾效率与准确度。

4.2 泰勒项数影响

本节计算条件与3.2 节一致，其中泰勒项数N分别取2、3 及5，x=0.0方向对比如图12 所示。

由图12 可知，随着泰勒项数的增加，对x=0.0 方 向上 σxx与 σyy最终计算结果的影响并不显著。泰勒项数的增加会导致计算效率降低，因此，在实际计算时应该同样应该合理选择确保效率。

4.3 权函数影响

本节计算条件与3.2 节一致，其中权函数除了指数分布外，同时引入如下所示的Gauss 函数、三次样条及四次样条权函数，其中d=|ξ|/δ，探究权函数对结果的影响，结果对比如图13 所示。

由图13 可知，Gauss 权函数计算结果与理论解偏差最小，该结论在 σyy边界区域较为明显，计算效果比样条权函数好，因此，在计算时可考虑采用Gauss 权函数。

5 结论

本文建立了基于非局部微分算子的近场动力学应力分析模型并对其有效性及收敛性进行了研究，取得以下结论：

(1) 基于非局部微分算子的近场动力学应力分析模型可以较好计算固体材料的应力分布状况，适用于完整、非完整材料，为经典近场动力学方法提供了应力分析思路。

(2) 离散间距及权函数对数值收敛结果具有显著影响，其中Gauss 权函数具有较好地收敛效果，同时计算精度随着离散间距的减小而提高但计算效率会降低，实际计算时应合理选择离散间距。