基于非局部LSM 优化的近场动力学及脆性材料变形模拟

马鹏飞，李树忱，王修伟，王曼灵

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

固体材料的力学特性始终是工程结构领域关注的热点，探究其失稳机制对工程防灾减灾领域意义重大[1-6]，采用理论方法探究裂纹扩展机制时通常需要大量不能真实反映工程问题，而利用实验方式时存在成本高、时间长等难题。随着计算机技术的发展，数值模拟方法越来越多的被用来探究结构变形及破坏问题。

有限元法[7-9]、离散元法[10-11]、无网格方法[12-13]等数值方法的出现为探究结构变形提供了低廉、快捷的途径，但在模拟材料裂纹扩展时会遇到尖端奇异性及裂纹成核等难题。扩展有限元理论采用预先设定裂纹的扩展方向及长度来解决结构裂纹扩展问题[14-15]。分子动力学方法可以较好的反应结构微观尺度，但是在处理宏观尺度问题时遇到计算效率低、计算时间长等难题[16]。

近场动力学是利用非局部作用思想来描述物质力学行为的方法[17]。该理论求解空间积分型运动方程而非微分方程可以很好的避免裂纹尖端奇异性等问题。近场动力学最早由Silling 等[18]提出，黄丹等[19-20]将其引入国内，针对混凝土变形破坏、巴西劈裂等问题进行了模拟。Zhou 等[21]利用近场动力学理论模拟了岩石类材料的裂纹扩展问题。王超等[22]在冰桨接触及模拟潜艇破冰上浮问题中应用了近场动力学。Zhu 等[23]运用该方法对不含不同预制倾角的岩石试件破坏过程进行了模拟。王超聪等[24]采用近场动力学理论模拟了热防护材料烧蚀过程。黄小华等[25]提出了单键双参数近场动力学模型并对冲击荷载对破坏的影响进行了分析。

以键为基础的微弹性模型(PMB)思路清晰应用广泛[26]。由于该模型过度简化，使得利用该模型时泊松比为定值，此外边界粒子作用区域不完整，计算时存在由“表面效应”引发的误差，同时计算时没有考虑非局部作用程度与作用距离之间的关系导致定量计算时精度受到影响。

本文结合前人的基础，将近场动力学中的微模量与弹性常数建立联系，通过刚度等效的方式建立线性方程组，并且利用求解不定方程组最小二乘最小范数的方式来矫正微模量，在此基础上引入可以反映非局部作用程度的核函数来提高计算精度，利用优化后的方法对二维脆性材料在荷载作用下的变形及裂纹扩展过程进行了模拟。结果表明，本文提出的方法可以有效地降低边界处的误差并减少误差产生的范围，在模拟裂纹扩展时会有更好的收敛效果，较好的应用于固体材料的变形与破坏。

1 理论与假定

1.1 近场动力学理论

图1 物质点相互作用Fig. 1 Interaction of material points

1.2 微弹性模型

1.3 小变形假定

图2 小变形假设Fig. 2 Small deformation assumption

可发现模量与文献[26]中相同，这可以证明本文小变形假定与键基近场动力学理论的一致性。c仅由体积模量K决定，导致该模型适用具有局限，另外在边界区域由于积分区域不完整会产生较大的计算误差。

2 非局部最小二乘优化

2.1 非局部优化

已证明模量c在传统模型中为常量，通常非局部作用程度会随距离而改变，使用常量c会导致定量计算时精度受到影响，Huang 等[28]针对该问题引入反映非局部作用程度的核函数从而提高定量计算的精度。本节在前人基础上，在小变形假设推导的理论基础上引入非局部核函数，以提高后续LSM 优化计算精度，式(9)作用力中的可重新表示为：

2.2 最小二乘方法优化

本节在小变形假定及非局部优化的基础上，将近场动力学模量与连续介质力学中弹性常数建立联系，并利用最小二乘方法优化键常数以此减少边界误差甚至突破适用局限性。式(12)可以表示为：

在经典连续介质力学二维问题中，应力-应变及刚度张量的一般形式为：

利用式(42)可建立传统连续介质力学刚度矩阵与近场动力学键模量之间的关系从而描述材料的变形特性，为后续裂纹扩展模拟奠定基础，类似式(43)可以将所有模量之间的关系表达为：

最小范数最小二乘解可采用Moore-Penrose 伪逆矩阵求解：

上述问题可采用二次规划(QR)处理，二次规划是非线性规划中的特殊规划问题，在约束最小二乘问题求解方面效果较好。以拉格朗日方法为例求解等式约束二次规划问题，考虑如下问题：

3 数值验证

3.1 单轴荷载

为验证本文方法的有效性，本节对二维平板在单轴荷载条件下的变形进行模拟，并将结果与理论及经典矫正结果做出对比。矩形板的尺寸及所受荷载如图3 所示。

图3 单轴加载模型 /mm Fig. 3 Uniaxial loading model

试件尺寸1 m×0.5 m，采用平面应力方式计算，试件弹性模量取E=20 GPa ，泊松比取ν=1/3，密度 取 ρ=2800 kg/m3，拉 伸 荷 载 取p0=20 MPa，以中心为坐标原点，试件在荷载作用下x及y方向位移的理论解为：

将模型离散为 200×100节点，节点间距为Δx=0.005 m ， 半径取 δ=3.015Δx，不采用任何矫正措施的结果与理论解的误差百分比如图4 所示。

图4 无矫正计算结果相对误差Fig. 4 Relative error of uncorrected calculation results

由图4 可知不利用任何矫正措施的模型存在较大的误差，特别是在角落及边界区域，在加载方向及非加载方向最大误差分别为21.33%与 32.05%，误差主要集中在角落处及加载端处。

边界产生大量误差的主要原因是边界附近的节点邻域不完整。目前较为常用的方法是利用应变能密度等价引入修正系数，利用该方法计算结果与误差百分比分别如图5 与图6 所示。

图5 应变能密度计算结果Fig. 5 Calculation results by strain energy density

图6 应变能密度相对误差Fig. 6 Relative error by strain energy density

计算解与理论解的误差百分比为了方便对比误差百分比低于0.1%时不考虑。由图6 可知，修正之后的计算结果较大改善，在角落与加载端部的误差显著降低，中部误差基本消失，加载方向的误差最大为1.72%，在非加载方向的误差最大为12.73%。

而利用本文提出的方法计算结果与误差百分比分别如图7 与图8 所示。不同方法计算结果的相对误差对比如表1 所示。

表1 最大相对误差对比表Table 1 Strain softening stress-strain curve of rock

图7 LSM 单轴加载结果Fig. 7 Uniaxial calculation results by LSM

图8 LSM 单轴相对误差Fig. 8 Uniaxial relative error by LSM

本文方法计算的相对误差与经典方法相比，在加载方向上最大相对误差没有较大的改善，但是误差产生的范围所减小。此外在非加载方向相比应变能密度方法有较大的改善，从12.73%降至4.35%并且产生误差的范围也有所降低，说明本文提出方法在模拟结构变形方面的可行性。

3.2 双轴荷载

为进一步验证方法的有效性，双轴荷载条件下的变形进行模拟，尺寸及所受荷载如图9 所示。

图9 双轴加载模型 /mmFig. 9 Biaxial loading model

两侧拉伸荷载取px=20 MPa，上下拉伸荷载取py=15 MPa，位移理论解可用叠加法得：

模型离散及节点间距等参数与3.1 节一致，经典应变能密度等价计算误差如图10 所示。

图10 双轴应变能密度相对误差Fig. 10 Biaxial relative error by strain energy density

经典方法在两个加载方向的误差最大分别为1.57%与11.71%，对比图6 与图10，双轴加载条件下边界的最大误差相比单轴加载有所降低产生误差的范围有较小的增加。利用本文的方法模拟结果与误差分布如图11 与图12 所示。

图11 LSM 双轴加载结果Fig. 11 Biaxial calculation results by LSM

本文方法在两个加载方向的误差最大分别为1.53%与3.29%。对比图8 与图12，双轴最大相对误差相比单轴有所降低，并且对比图10 与图12，最大误差比经典方法具有改善，尤其在最小主应力方向，且产生误差的范围也有所减少，进一步验证了本文方法的有效性。

图12 LSM 双轴相对误差Fig. 12 Biaxial relative error by LSM

3.3 裂纹扩展

为验证本文方法研究材料断裂的有效性，本节对含裂纹的脆性材料在荷载作用下的动态破裂过程进行模拟并将结果与文献[29]对比，尺寸及所受荷载如图13 所示。

图13 裂纹扩展模型[29] /mmFig. 13 Crack growth model[29]

参数与文献[29]相同取弹性模量E=72 GPa，泊松比ν=0.22 ，密度 ρ =2440 kg/m3，两侧拉伸荷载p0=14 MPa ，断裂能G=135 J/m2。

为探究相同计算成本下的收敛结果，模型离散选取与一组参考文献类似的节点数 200×80，间距为 Δx=0.5 mm ， 半径取 δ=3.015Δx，不同时步nt下裂纹扩展情况以及裂纹扩展速度随时间的变化对比分别如图14 与图15 所示。

图14 不同时步损伤分布Fig. 14 Damage distribution at different time steps

图15 裂纹扩展速度对比Fig. 15 Crack growth speed comparison

对比本文模拟结果与文献[29]结果可知本文优化的方法可以较好的反映脆性材料的裂纹扩展过程，捕捉到裂纹分叉的现象。此外由图15 可知，在大致相等的计算成本下本文方法得到的裂纹扩展速度收敛效果相比文献[29]有所提升，这其中表现在收敛速度与收敛结果上，利用本文方法在同等计算成本下可更接近实验测得值。

4 结论

本文建立了基于非局部最小二乘优化的近场动力学数值模型并对脆性材料的变形及破坏进行了模拟，取得以下结论：

(1) 基于最小二乘最小范数最优解的近场动力学方法可以较好模拟固体材料的变形与破坏。与经典的近场动力学方法相比，本文方法可以较好的降低边界处的最大相对误差及产生误差的范围，同时可以较为准确的反映脆性材料的裂纹扩展过程。

(2) 在荷载条件下，利用近场动力学方法所模拟的结果与理论结果所产生的误差在两个方向分别存在于角落与最大主应力边界，并且双轴加载条件下产生的最大误差相比单轴条件有所降低。