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

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

doi:10.6052/j.issn.1000-4750.2021.03.0188

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

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

基金项目: 国家自然科学基金项目(51879150，41831278)

详细信息

作者简介:
马鹏飞(1994−)，男，山东人，博士生，主要从事深部岩体破坏机理研究(E-mail: mapengfeisdu@163.com)

王修伟(1997−)，男，山东人，博士生，主要从岩土工程数值模拟方面的研究(E-mail: wxwsdu@126.com)

王曼灵(1994−)，女，河南人，博士生，主要从岩土工程数值模拟方面的研究(E-mail: manling_wang@163.com)

通讯作者:
李树忱(1973−)，男，黑龙江人，教授，博士，博导，主要从事深部岩体破坏机理研究(E-mail: shuchenli@sdu.edu.cn)

中图分类号: O346.1+1
计量
- 文章访问数: 495
- HTML全文浏览量: 188
- PDF下载量: 122
出版历程
- 收稿日期: 2021-03-13
- 修回日期: 2021-05-30
- 网络出版日期: 2021-06-17
- 刊出日期: 2022-05-26

PERIDYNAMIC METHOD BASED ON NONLOCAL LSM OPTIMIZATION AND DEFORMATION SIMULATION OF BRITTLE MATERIALS

Geotechnical and Structural Engineering Research Center, Shandong University, Jinan, Shandong 250061, China

摘要

摘要: 在经典近场动力学模型的基础上，通过小变形假定将近场动力学中的微模量与经典理论中的弹性常数建立联系，引入可以反映非局部作用特性的核函数提高计算精度，利用刚度等效的方式建立有关微模量的线性方程组，并通过寻求不定线性方程组最小二乘(LSM)最小范数解的方式对近场动力学中微模量进行优化，根据二次规划得到最优非负模量。利用优化后的方法对二维平板在单轴和双轴荷载作用下的变形及含预制裂纹脆性材料在荷载下的裂纹扩展进行了模拟并将结果与理论经典近场动力学方法结果对比。结果表明：优化后的方法可以较好的反映结构在荷载条件下的变形与破坏特性，与经典方法相比材料变形模拟在最大误差及误差范围具有良好的改善，并且模拟裂纹扩展过程在同等计算成本下具有更优的收敛速度及收敛结果，进一步验证了所提出方法的有效性，有着较为广泛的应用前景。
- 数值模拟 /
- 非局部方法 /
- 近场动力学 /
- 最小二乘法 /
- 材料变形
Abstract: Based on the peridynamic model, the relationship between the micro modulus in peridynamics and the elastic constants in classical theory is established by assuming small deformation. The kernel function which can reflect the nonlocal action characteristics is introduced to improve the calculation accuracy, and the linear equations related to the micro modulus are established by means of stiffness equivalence. The least square minimum (LSM) norm solution of uncertain linear equations is used to optimize the micro modulus in peridynamics, and the optimal nonnegative modulus is obtained through quadratic programming. The deformation of two-dimensional plate under uniaxial and biaxial loads and the crack growth of brittle material with prefabricated cracks under loads are simulated by the optimized method, and the results are compared with those by classical peridynamics method. The results show that the optimized method can better reflect the deformation and failure characteristics of the structure under load conditions, and the maximum error and error range of material deformation simulation are better than the classical method; and moreover, the simulation of crack growth process has better convergence speed and convergence results under the same calculation cost, which verifies the effectiveness and the wide application prospect of the proposed method.
- numerical simulation /
- nonlocal method /
- peridynamics /
- least square method /
- material deformation

HTML全文

固体材料的力学特性始终是工程结构领域关注的热点，探究其失稳机制对工程防灾减灾领域意义重大^[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所示，非局部作用是指在物质点在一定距离内存在相互作用力， $x'$ 为距 $x$ 半径 $\delta$ 内非局部作用点，参考坐标系中两点相对位置为 $x'{\rm{ - }}x{\rm{ = }}\xi$ ，点 $x$ 在 $t$ 时刻的平衡方程为：

$\rho ( x )\ddot u( {x,t} ) = \int_{{H_x}} {f( {x',x,u( {x',t} ),u( {x,t} ),t} )} {\rm{d}}V' + b( {x,t} )$

(1)

式中： $u\left( {x,t} \right)$ 与 $u\left( {x',t} \right)$ 分别为 $x$ 与 $x'$ 变形时的位移， $u\left( {x',t} \right){\rm{ - }}u\left( {x,t} \right){\rm{ = }}\eta$ 为变形后的相对位移，因此变形后位置为 $\xi {\rm{ + }}\eta$ ； $\rho \left( x \right)$ 与 $\ddot u\left( {x,t} \right)$ 分别 $x$ 的密度与加速度； $f$ 为 $x$ 与 $x'$ 相互作用力； $b\left( {x,t} \right)$ 为 $t$ 时刻外力荷载。

图 1 物质点相互作用

Figure 1. Interaction of material points

下载: 全尺寸图片幻灯片

1.2 微弹性模型

微弹性模型中假定作用力取决于相对位置 $\xi$ 与相对位移 $\eta$ ，而且 $x$ 与 $x'$ 之间 $f$ 大小相等方向相反，式(1)可以写为：

$\rho \left( x \right)\ddot u\left( {x,t} \right) = \int_{{H_x}} {f\left( {u' - u,x' - x} \right)} {\rm{d}}V' + b\left( {x,t} \right)$

(2)

参考文献[23]微弹性模型作用力函数f可由微势能得到：

$f\left( {\xi {\rm{,}}\eta } \right){\rm{ = }}\frac{{\partial \omega (\xi ,\eta )}}{{\partial \eta }},\;\;\;\forall \xi ,\eta$

(3)

微势能表示储存在材料中单位体积平方的能量，因此， $x$ 处的应变能密度可表示为：

$W = \frac{1}{2}\int_{{H_x}} {\omega (\xi ,\eta ){\rm{d}}V'}$

(4)

式中：1/2为每个键相对于 $x$ 储存一半的能量；PMB模型中微势能 $\omega$ 可表示为：

$\omega (\xi ,\eta ) = \frac{{c\left( \xi \right){s^2}\left| \xi \right|}}{2}$

(5)

式中： $c\left( \xi \right)$ 为键 $\xi$ 的微模量；对于经典的微弹性模型通常取常数 $c$ ； $s$ 为体现物质点间的变形伸长率：

$s = \frac{{\left| {\xi + \eta } \right| - \left| \xi \right|}}{{\left| \eta \right|}}\qquad$

(6)

式中， $\left| \xi \right|$ 与 $\left| {\xi + \eta } \right|$ 分别为变形前后的长度，当超过临界值 ${s_{\rm{0}}}$ 时键会发生不可逆的断裂。而力的方向与变形后相对位置相同，因此， $f$ 可表示为：

$f\left( {\xi {\rm{,}}\eta } \right){\rm{ = }}cs\mu \frac{{\xi + \eta }}{{\left| {\xi + \eta } \right|}}$

(7)

式(7)中， $\mu$ 为控制相互作用键断裂的函数。此外近场动力学理论中无需引入外部准则，并且由 $x$ 处的破坏情况定义损伤：

$D\left( x \right) = 1 - \frac{{\displaystyle\int\nolimits_{{H_x}} {\mu {\rm{d}}V} }}{{\displaystyle\int\nolimits_{{H_x}} {{\rm{d}}V} }}$

(8)

式中， $D$ 的取值为0~1，0代表 $x$ 处还未破坏，而1代表完全破坏。

1.3 小变形假定

如1.2节所述微弹性模型中作用力 $f$ 及键伸长率 $s$ 均为变形后相对位移的 $\eta$ 的非线性函数。而在弹性力学中，应力及应变可表示为变形的线性函数，因此，结合经典弹性力学及参考文献[27]可引入小变形假设 $\eta \ll {\rm{1}}$ ，将式(7) $f$ 表示为：

$f\left( {\xi {\rm{,}}\eta } \right){\rm{ = }}c\mu \frac{{{\eta _n}}}{\xi }\frac{\xi }{{\left| \xi \right|}}$

(9)

式中， ${\eta _n}$ 为未变形键 $\xi$ 方向上的相对位移分量，对比式(7)与式(9)可发现 ${\eta _n}{\rm{/}}\left| \xi \right|$ 为未变形键 $\xi$ 方向上的线性伸长率，与弹性力学中应变的定义类似，将式(9)代入式(2)不考虑键断裂函数时得控制方程为：

$\rho \left( x \right)\ddot u\left( {x,t} \right) = \int_{{H_x}} {c\frac{{{\eta _n}}}{\xi }\frac{\xi }{{\left| \xi \right|}}} {\rm{d}}V' + b\left( {x,t} \right)$

(10)

由于该小变形假设是基于微弹性模型，因此将式(9)代入式(3)并对 $\eta$ 积分可得微势能 $\omega$ ：

$\omega (\xi ,\eta ) = \frac{{\rm{1}}}{2}\frac{{c\eta _n^2}}{{\left| \xi \right|}}\qquad\quad$

(11)

将微势能 $\omega$ 代入式(4)可得 $x$ 处应变能密度：

${W^{PD}} = \frac{1}{2}\int_{{H_x}} {\frac{{\rm{1}}}{2}\frac{{c\eta _n^2}}{{\left| \xi \right|}}{\rm{d}}V'}$

(12)

对于二维问题，若相互作用区域为半径 $\delta$ 得圆形，则式(13)可表示为：

${W^{PD}} = \frac{c}{4}\int_0^{2\pi } {\int_0^\delta {\frac{{\eta _n^2}}{{\left| \xi \right|}}} } \left| \xi \right|{\rm{d}}\left| \xi \right|{\rm{d}}\theta$

(13)

式(13)中 ${\eta _n}$ 为唯一的未知量，假设应变无穷小，如图2所示，对于各向同性材料 ${\eta _n}$ 可写为：

${\eta _n} = \left| \xi \right|( {{\varepsilon _{11}}{{\cos }^2}\theta + {\varepsilon _{22}}{{\sin }^2}\theta } )$

(14)

图 2 小变形假设

Figure 2. Small deformation assumption

下载: 全尺寸图片幻灯片

代入式(14)化简可得：

${W^{PD}} = \frac{c}{4}\int_0^{2\pi } {\int_0^\delta {{{( {{\varepsilon _{11}}{{\cos }^2}\theta + {\varepsilon _{22}}{{\sin }^2}\theta } )}^{\rm{2}}}} } {| \xi |^{\rm{2}}}{\rm{d}}| \xi |{\rm{d}}\theta$

(15)

这里，可以假设 ${\varepsilon _{{\rm{11}}}}{\rm{ = }}{\varepsilon _{{\rm{22}}}}{\rm{ = }}\varepsilon$ 并积分式(15)得：

${W^{PD}} = \frac{{\pi c{\delta ^3}{\varepsilon ^2}}}{{\rm{6}}}$

(16)

二维平面应力条件下，应变能密度可表示为：

${W^{EM}} = \frac{E}{{2( {1 - {\nu ^2}} )}}( {\varepsilon _{11}^2 + \varepsilon _{22}^2} ) + \frac{E}{{1 - {\nu ^2}}}{\varepsilon _{{\rm{11}}}}{\varepsilon _{{\rm{22}}}}$

(17)

同样，假定 ${\varepsilon _{{\rm{11}}}}{\rm{ = }}{\varepsilon _{{\rm{22}}}}{\rm{ = }}\varepsilon$ ，与式(16)等价得：

$c = \frac{{6E}}{{\pi {\delta ^3}\left( {1 - \nu } \right)}}$

(18)

在平面应力问题中体积模量 $K$ 可表示：

$K = \frac{E}{{2\left( {1 - \nu } \right)}}\;\;\;$

(19)

代入式(22)得：

$c = \frac{{{\rm{12}}K}}{{\pi {\delta ^{\rm{3}}}}}\qquad$

(20)

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

2 非局部最小二乘优化

2.1 非局部优化

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

$f\left( {\xi {\rm{,}}\eta } \right){\rm{ = }}c\left( \xi \right)\mu \frac{{{\eta _n}}}{\xi }\frac{\xi }{{\left| \xi \right|}}$

(21)

式中， $c\left( \xi \right)$ 为随变形前相对位置 $\xi$ 变化函数：

$c\left( \xi \right){\rm{ = }}c\left( {0,\delta } \right)g\left( {\xi ,\delta } \right)$

(22)

式中， $g\left( {\xi ,\delta } \right)$ 是可反映非局部作用程度的核函数修正项，需满足条件可参考文献[25]，本文核函数修正项选用线形核函数 $g$ ：

$g\left( {\xi ,\delta } \right){\rm{ = 1 - }}\frac{{\left| \xi \right|}}{\delta }$

(23)

将式(23)代入式(22)再代入式(21)，暂不考虑键断裂的情况：

$f\left( {\xi {\rm{,}}\eta } \right){\rm{ = }}c\left( {0,\delta } \right)\left( {{\rm{1 - }}\frac{{\left| \xi \right|}}{\delta }} \right)\frac{{{\eta _n}}}{\xi }\frac{\xi }{{\left| \xi \right|}}$

(24)

代入式(3)对 $\eta$ 积分可得微势能函数 $\omega$ ：

$\omega (\xi ,\eta ) = \frac{{\rm{1}}}{2}c\left( {0,\delta } \right)\left( {{\rm{1 - }}\frac{{\left| \xi \right|}}{\delta }} \right)\frac{{\eta _n^2}}{{\left| \xi \right|}}$

(25)

将微势能 $\omega$ 代入式(4)可得应变能密度：

${W^g} = \frac{1}{2}\int_{{H_x}} {\frac{{\rm{1}}}{2}c\left( {0,\delta } \right)\left( {{\rm{1 - }}\frac{{\left| \xi \right|}}{\delta }} \right)\frac{{\eta _n^2}}{{\left| \xi \right|}}{\rm{d}}V'}$

(26)

作用域为 $\delta$ 的二维问题，式(26)可表示为：

${W^g} = \frac{{c\left( {0,\delta } \right)}}{4}\int_0^{2\pi } {\int_0^\delta {\left( {{\rm{1 - }}\frac{{\left| \xi \right|}}{\delta }} \right)\frac{{\eta _n^2}}{{\left| \xi \right|}}} } \left| \xi \right|{\rm{d}}\left| \xi \right|{\rm{d}}\theta$

(27)

将式(14)代入式(27)并假定 ${\varepsilon _{{\rm{11}}}}{\rm{ = }}{\varepsilon _{{\rm{22}}}}{\rm{ = }}\varepsilon$ 得：

${W^g} = \frac{{\pi c\left( {0,\delta } \right){\delta ^3}{\varepsilon ^2}}}{{{\rm{24}}}}$

(28)

由应变能密度等价，式(28)与式(17)相等得：

$c\left( {0,\delta } \right) = \frac{{{\rm{24}}E}}{{\pi {\delta ^3}\left( {1 - \nu } \right)}}$

(29)

将式(29)代入式(22)：

$c\left( \xi \right){\rm{ = }}\frac{{{\rm{24}}E}}{{\pi {\delta ^3}\left( {1 - \nu } \right)}}\left( {{\rm{1 - }}\frac{{\left| \xi \right|}}{\delta }} \right)$

(30)

从式(30)可看出，采用核函数修正项可以反映非局部作用程度，模量 $c\left( \xi \right)$ 越大，程度越强。

2.2 最小二乘方法优化

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

${W^{PD}} = \frac{1}{2}\int_{{H_x}} {\frac{{\rm{1}}}{2}\frac{{c\left( \xi \right)\eta _n^2}}{{\left| \xi \right|}}{\rm{d}}V'}$

(31)

式中，模量 $c\left( \xi \right)$ 为式(30)、式(31)对应变 ${\varepsilon _{kl}}$ 求导可得：

$\frac{{\partial {W^{PD}}}}{{\partial {\varepsilon _{kl}}}} = \frac{1}{2}\int_{{H_x}} {\frac{{c\left( \xi \right){\eta _n}}}{{\left| \xi \right|}}\frac{{\partial {\eta _n}}}{{\partial {\varepsilon _{kl}}}}{\rm{d}}V'}$

(32)

变形量 ${\eta _n}$ 对应变 ${\varepsilon _{kl}}$ 的导数可由变形前相对位置及应变 ${\varepsilon _{kl}}$ 表示：

$\frac{{\partial {\eta _n}}}{{\partial {\varepsilon _{kl}}}} = \frac{{{\xi _k}{\varepsilon _{kl}}{\xi _l}}}{{\left| \xi \right|}}$

(33)

代入式(32)并再次对应变求导得：

$\frac{{\partial {W^{PD}}}}{{\partial {\varepsilon _{ij}}\partial {\varepsilon _{kl}}}} = \int_{{H_x}} {\frac{{c\left( \xi \right)}}{2}\frac{{{\xi _i}{\xi _j}{\xi _k}{\xi _l}}}{{{{\left| \xi \right|}^3}}}{\rm{d}}V'}$

(34)

在经典连续介质力学中应变能密度与弹性四阶刚度张量关系为：

$\frac{{\partial {W^{CM}}}}{{\partial {\varepsilon _{ij}}\partial {\varepsilon _{kl}}}} = {\boldsymbol{C}}_{ijkl}^{CM}$

(35)

通过式(34)与式(35)相等可得：

${\boldsymbol{C}}_{ijkl}^{CM} = \int_{{H_x}} {\frac{{c\left( \xi \right)}}{2}\frac{{{\xi _i}{\xi _j}{\xi _k}{\xi _l}}}{{{{\left| \xi \right|}^3}}}{\rm{d}}V'}$

(36)

也可表示为张量形式：

${{\boldsymbol{C}}^{CM}} = \int_{{H_x}} {\frac{{c\left( \xi \right)}}{2}\frac{{\xi \otimes \xi \otimes \xi \otimes \xi }}{{{{\left| \xi \right|}^3}}}{\rm{d}}V'}$

(37)

对于二维问题式(36)中 $i,j,k,l$ 取值的范围为1到2，其离散形式为：

${\boldsymbol{C}}_{ijkl}^{CM} = \sum\limits_{N = 1}^M {\frac{{{c^N}\left( \xi \right)}}{2}\frac{{\xi _i^N\xi _j^N\xi _k^N\xi _l^N}}{{{{\left| {{\xi ^N}} \right|}^3}}}} \Delta {V_N}$

(38)

式中： $M$ 为点 $x$ 邻域 ${H_x}$ 内相互作用点总数； $N$ 为第 $N$ 作用点，例如 $C_{{\rm{1212}}}^{CM}$ 可以表示为：

${\boldsymbol{C}}_{{\rm{1212}}}^{CM} = \sum\limits_{N = 1}^M {\frac{{{c^N}\left( \xi \right)}}{2}\frac{{\xi _{\rm{1}}^N\xi _{\rm{2}}^N\xi _{\rm{1}}^N\xi _{\rm{2}}^N}}{{{{\left| {{\xi ^N}} \right|}^3}}}} \Delta {V_N}$

(39)

利用刚度矩阵对称性， $i,j,k,l$ 可采用Voigt标记法，将 $2 \times 2 \times 2 \times 2$ 四阶张量表示为 ${\rm{3}} \times {\rm{3}}$ 的二阶张量其中 ${\zeta _{\rm{1}}}{\rm{ = }}{\xi _{\rm{1}}}{\xi _{\rm{1}}}$ 、 ${\zeta _{\rm{2}}}{\rm{ = }}{\xi _{\rm{2}}}{\xi _{\rm{2}}}$ 及 ${\zeta _{\rm{3}}}{\rm{ = }}{\xi _{\rm{1}}}{\xi _{\rm{2}}}$ ，因此：

${\boldsymbol{C}}_{\alpha \beta }^{CM} = \sum\limits_{N = 1}^M {\frac{{{c^N}\left( \xi \right)}}{2}\frac{{\zeta _\alpha ^N\zeta _\beta ^N}}{{{{\left| {{\xi ^N}} \right|}^3}}}} \Delta {V_N}$

(40)

引入非局部作用核函数，将式(22)代入式(40)：

${\boldsymbol{C}}_{\alpha \beta }^{CM} = \sum\limits_{N = 1}^M {\frac{{{c^N}\left( {0,\delta } \right){g^N}\left( {\xi ,\delta } \right)}}{2}\frac{{\zeta _\alpha ^N\zeta _\beta ^N}}{{{{\left| {{\xi ^N}} \right|}^3}}}} \Delta {V_N}$

(41)

对于每个作用点对应核函数 ${g^N}\left( {\xi ,\delta } \right)$ 为反映非局部作用程度的固定值，因此将式(41)分离为：

${\boldsymbol{C}}_{\alpha \beta }^{CM} = \left[ {\sum\limits_{N = 1}^M {\frac{{{g^N}\left( {\xi ,\delta } \right)}}{2}\frac{{\zeta _\alpha ^N\zeta _\beta ^N}}{{{{\left| {{\xi ^N}} \right|}^3}}}} \Delta {V_N}} \right]{c^N}\left( {0,\delta } \right)$

(42)

刚度项 $C_{\alpha \beta }^{CM}$ 可由两向量相乘得到，例如 $C_{{\rm{12}}}^{CM}$ 可表示为：

$C_{{\rm{12}}}^{CM} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{{g^1}}}{2}\dfrac{{\zeta _{\rm{1}}^{\rm{1}}\zeta _{\rm{2}}^{\rm{1}}}}{{{{\left| {{\xi ^{\rm{1}}}} \right|}^3}}}\Delta {V_{\rm{1}}}}&{\dfrac{{{g^2}}}{2}\dfrac{{\zeta _{\rm{1}}^2\zeta _{\rm{2}}^{\rm{2}}}}{{{{\left| {{\xi ^2}} \right|}^3}}}\Delta {V_2}}& \ldots &{\dfrac{{{g^M}}}{2}\dfrac{{\zeta _{\rm{1}}^M\zeta _{\rm{2}}^{\rm{M}}}}{{{{\left| {{\xi ^M}} \right|}^3}}}\Delta {V_M}} \end{array}} \right] \left[ {\begin{array}{*{20}{c}} {{c^{\rm{1}}}\left( {0,\delta } \right)} \\ {{c^2}\left( {0,\delta } \right)} \\ \vdots \\ {{c^M}\left( {0,\delta } \right)} \end{array}} \right]$

(43)

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

$\left[ {\begin{array}{*{20}{c}} {{\sigma _{{\rm{11}}}}} \\ {{\sigma _{{\rm{22}}}}} \\ {{\sigma _{{\rm{12}}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{C_{11}}}&{{C_{12}}}&{{C_{13}}} \\ {{C_{21}}}&{{C_{22}}}&{{C_{23}}} \\ {{C_{31}}}&{{C_{32}}}&{{C_{33}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\varepsilon _{11}}} \\ {{\varepsilon _{22}}} \\ {{\varepsilon _{12}}} \end{array}} \right]$

(44)

中间即为刚度张量 $C_{\alpha \beta }^{CM}$ ，对于各向同性的弹性材料平面应力问题刚度张量为：

$C_{\alpha \beta }^{CM} = \frac{E}{{1 - {\nu ^2}}}\left[ {\begin{array}{*{20}{c}} 1&\nu &0 \\ \nu &1&0 \\ 0&0&{\left( {1 - \nu } \right)/2} \end{array}} \right]$

(45)

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

$\left[ {\begin{array}{*{20}{c}} {C_{{\rm{11}}}^{CM}} \\ {C_{{\rm{12}}}^{CM}} \\ {C_{{\rm{13}}}^{CM}} \\ \vdots \\ {C_{{\rm{33}}}^{CM}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\dfrac{{{g^1}}}{2}\dfrac{{\zeta _{\rm{1}}^{\rm{1}}\zeta _1^{\rm{1}}}}{{{{\left| {{\xi ^{\rm{1}}}} \right|}^3}}}\Delta {V_{\rm{1}}}}&{\dfrac{{{g^2}}}{2}\dfrac{{\zeta _{\rm{1}}^2\zeta _1^{\rm{2}}}}{{{{\left| {{\xi ^2}} \right|}^3}}}\Delta {V_{\rm{2}}}}&{\dfrac{{{g^{\rm{3}}}}}{2}\dfrac{{\zeta _{\rm{1}}^{\rm{3}}\zeta _1^{\rm{3}}}}{{{{\left| {{\xi ^{\rm{3}}}} \right|}^3}}}\Delta {V_{\rm{3}}}}& \cdots &{\dfrac{{{g^M}}}{2}\dfrac{{\zeta _{\rm{1}}^M\zeta _1^M}}{{{{\left| {{\xi ^M}} \right|}^3}}}\Delta {V_M}} \\ {\dfrac{{{g^1}}}{2}\dfrac{{\zeta _{\rm{1}}^{\rm{1}}\zeta _{\rm{2}}^{\rm{1}}}}{{{{\left| {{\xi ^{\rm{1}}}} \right|}^3}}}\Delta {V_{\rm{1}}}}&{\dfrac{{{g^2}}}{2}\dfrac{{\zeta _{\rm{1}}^2\zeta _{\rm{2}}^{\rm{2}}}}{{{{\left| {{\xi ^2}} \right|}^3}}}\Delta {V_{\rm{2}}}}&{\dfrac{{{g^{\rm{3}}}}}{2}\dfrac{{\zeta _{\rm{1}}^{\rm{3}}\zeta _{\rm{2}}^{\rm{3}}}}{{{{\left| {{\xi ^{\rm{3}}}} \right|}^3}}}\Delta {V_3}}& \cdots &{\dfrac{{{g^M}}}{2}\dfrac{{\zeta _{\rm{1}}^M\zeta _{\rm{2}}^M}}{{{{\left| {{\xi ^M}} \right|}^3}}}\Delta {V_M}} \\ {\dfrac{{{g^1}}}{2}\dfrac{{\zeta _{\rm{1}}^{\rm{1}}\zeta _3^{\rm{1}}}}{{{{\left| {{\xi ^{\rm{1}}}} \right|}^3}}}\Delta {V_1}}&{\dfrac{{{g^2}}}{2}\dfrac{{\zeta _{\rm{1}}^2\zeta _3^{\rm{2}}}}{{{{\left| {{\xi ^2}} \right|}^3}}}\Delta {V_{\rm{2}}}}&{\dfrac{{{g^{\rm{3}}}}}{2}\dfrac{{\zeta _{\rm{1}}^{\rm{3}}\zeta _3^{\rm{3}}}}{{{{\left| {{\xi ^{\rm{3}}}} \right|}^3}}}\Delta {V_3}}& \cdots &{\dfrac{{{g^M}}}{2}\dfrac{{\zeta _{\rm{1}}^M\zeta _3^M}}{{{{\left| {{\xi ^M}} \right|}^3}}}\Delta {V_M}} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ {\dfrac{{{g^1}}}{2}\dfrac{{\zeta _3^{\rm{1}}\zeta _3^{\rm{1}}}}{{{{\left| {{\xi ^{\rm{1}}}} \right|}^3}}}\Delta {V_1}}&{\dfrac{{{g^2}}}{2}\dfrac{{\zeta _3^2\zeta _3^{\rm{2}}}}{{{{\left| {{\xi ^2}} \right|}^3}}}\Delta {V_{\rm{2}}}}&{\dfrac{{{g^{\rm{3}}}}}{2}\dfrac{{\zeta _{\rm{3}}^{\rm{3}}\zeta _3^{\rm{3}}}}{{{{\left| {{\xi ^{\rm{3}}}} \right|}^3}}}\Delta {V_3}}& \cdots &{\dfrac{{{g^M}}}{2}\dfrac{{\zeta _3^M\zeta _3^M}}{{{{\left| {{\xi ^M}} \right|}^3}}}\Delta {V_M}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{c^{\rm{1}}}\left( {0,\delta } \right)} \\ {{c^2}\left( {0,\delta } \right)} \\ {{c^3}\left( {0,\delta } \right)} \\ \vdots \\ {{c^M}\left( {0,\delta } \right)} \end{array}} \right]$

(46)

式(46)可以简写为矩阵形式：

${{\boldsymbol{C}}^{CM}} = {\boldsymbol{Kc}}$

(47)

式中： ${{\boldsymbol{C}}^{CM}}$ 为所计算结构的刚度张量； ${\boldsymbol{K}}$ 为研究点 $x$ 与邻域 ${H_x}$ 相互作用点的常数项； ${\boldsymbol{c}}$ 为模量项。在模拟时 ${{\boldsymbol{C}}^{CM}}$ 、 ${\boldsymbol{K}}$ 为已知项， ${\boldsymbol{c}}$ 为未知项。而 ${{\boldsymbol{C}}^{CM}}$ 为 $n$ 行 $M$ 列矩阵，其中 $n$ 为刚度张量项数， $M$ 为 $x$ 邻域 ${H_x}$ 相互作用点个数。我们期望 $n$ 与 $M$ 相等，式(46)为适定方程组，模量项 ${\boldsymbol{c}}$ 具有唯一解。但对于二维或三维问题，作用点个数通常大于刚度张量项数，模量项 ${\boldsymbol{c}}$ 有无数组解，针对该方程组可以寻求最小范数及最小二乘误差意义上的最优解：

$\min {\rm{ }}\left\| {\boldsymbol{c}} \right\|\qquad\qquad$

(48)

$\min {\rm{ }}{\| {{{\boldsymbol{C}}^{CM}}{\rm{ - }}{\boldsymbol{Kc}}} \|^{\rm{2}}}$

(49)

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

${{\boldsymbol{A}}^ + } = \mathop {\lim }\limits_{\alpha \to 0} {( {{{\boldsymbol{A}}^{\rm{T}}}A + \alpha I} )^{ - 1}}{{\boldsymbol{A}}^{\rm{T}}}$

(50)

式中， ${\boldsymbol{A}}$ 为线性方程 $Ax = y$ 中非方矩阵，计算时一般采用：

${{\boldsymbol{A}}^ + } = {\boldsymbol{V}}{{\boldsymbol{D}}^ + }{{\boldsymbol{U}}^{\rm{T}}}\quad\qquad\qquad$

(51)

矩阵 ${\boldsymbol{U}}$ 、 ${\boldsymbol{D}}$ 与 ${\boldsymbol{ V}}$ 是矩阵 ${\boldsymbol{A}}$ 奇异值分解后得到的矩阵，对角矩阵 ${\boldsymbol{D}}$ 的伪逆 ${{\boldsymbol{D}}^ + }$ 是其非零元素取倒数之后再转置得到的。

除了寻求误差最小意义上的最优解外，还应该考虑其他定解条件，储存在键中的能量非负，所以模量项 ${\boldsymbol{c}}$ 中的每一项应：

${\boldsymbol{c}} = {\boldsymbol{c}}:{c_i} \geqslant 0\qquad\qquad\quad$

(52)

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

$\begin{split} & \min {\rm{ }}\frac{1}{2}{x^{\rm{T}}}Hx + {c^{\rm{T}}}x, \\& {\rm s.t.}\;\;\;\;\;\;Ax \geqslant b \end{split}$

(53)

矩阵H对称正定，A行满秩，拉格朗日函数为：

$L\left( {x,\lambda } \right) = \frac{1}{2}{x^{\rm{T}}}Hx + {c^{\rm{T}}}x - {\lambda ^{\rm{T}}}\left( {Ax - b} \right)$

(54)

可令：

${\nabla _x}L\left( {x,\lambda } \right) = 0,\;\;{\nabla _{\lambda} }L\left( {x,\lambda } \right) = 0$

(55)

将上述方程组写成分块矩阵形式：

$Hx - {A^{\rm{T}}}\lambda = - c ,\;\; - Ax = - b$

(56)

此方程组的系数矩阵为拉格朗日矩阵：

$\left[ {\begin{array}{*{20}{c}} H&{ - {A^{\rm{T}}}} \\ { - A}&0 \end{array}} \right]$

(57)

为了满足式(52)模量项非负的附加条件及利用QR求解最优解 ${\boldsymbol{c}}$ ，可构造 ${\boldsymbol{r}}$ 使其满足：

${\boldsymbol{Kr}} = {\boldsymbol{0}}\qquad\qquad$

(58)

将式(58)代入式(47)可得：

${{\boldsymbol{C}}^{CM}} = {\boldsymbol{K}}\left( {{\boldsymbol{c}}{\rm{ + }}{\boldsymbol{r}}} \right)$

(59)

此时可将问题转换为求解 ${\boldsymbol{c}}{\rm{ + }}{\boldsymbol{r}}$ 的最小范数最小二乘解，并且附加 ${\boldsymbol{c}}{\rm{ + }}{\boldsymbol{r}} \geqslant {\rm{0}}$ 的定解条件。该过程可分为两步执行，首先基于式(47)求解 ${\boldsymbol{c}}$ 的最小范数最小二乘解，随后引入满足式(58)的 ${\boldsymbol{r}}$ 求解 ${\boldsymbol{c}}{\rm{ + }}{\boldsymbol{r}}$ ，式(58)可保证：

$\min {\rm{ }}{\left\| {{{\boldsymbol{C}}^{CM}}{\rm{ - }}{\boldsymbol{K}}\left( {{\boldsymbol{c}}{\rm{ + }}{\boldsymbol{r}}} \right)} \right\|^{\rm{2}}}$

(60)

而 ${\boldsymbol{c}}{\rm{ + }}{\boldsymbol{r}}$ 的最小范数 $\min \left\| {{\boldsymbol{c}}{\rm{ + }}{\boldsymbol{r}}} \right\|$ 可等价于求解：

$\min {\rm{ }}\frac{{\rm{1}}}{{\rm{2}}}{\left( {{\boldsymbol{c}}{\rm{ + }}{\boldsymbol{r}}} \right)^{\rm{T}}}\left( {{\boldsymbol{c}}{\rm{ + }}{\boldsymbol{r}}} \right)\quad$

(61)

式(61)展开化简可得：

$\min {\rm{ }}\frac{{\rm{1}}}{{\rm{2}}}{{\boldsymbol{r}}^{\rm{T}}}{\boldsymbol{Ir}}{\rm{ + }}{{\boldsymbol{c}}^{\rm{T}}}{\boldsymbol{r}}\quad\quad$

(62)

${\boldsymbol{I}}$ 为 $M \times M$ 的单位矩阵，并 ${{\boldsymbol{c}}^{\rm T}}{\boldsymbol{c}}$ 是已知项对最小值无影响，至此可转化为二次规划问题求解：

$\begin{split} & \min {\rm{ }}\frac{{\rm{1}}}{{\rm{2}}}{{\boldsymbol{r}}^{\rm{T}}}{\boldsymbol{Ir}}{\rm{ + }}{{\boldsymbol{c}}^{\rm{T}}}{\boldsymbol{r}} ,\\& {\rm s.t.}\;\;\;{\boldsymbol{c}}{\rm{ + }}{\boldsymbol{r}} \geqslant {\boldsymbol{0}}, \\& \;\;\;\;\;\quad{\boldsymbol{Kr}} = {\boldsymbol{0}} \end{split} \quad\quad$

(63)

QR求解后的 ${\boldsymbol{c}}{\rm{ + }}{\boldsymbol{r}}$ 为研究点 $x$ 的正模量项，保证储存键中非负的应变能，此外若第一步求解 ${\boldsymbol{c}}$ 为非负则无需第二步计算。

3 数值验证

3.1 单轴荷载

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

图 3 单轴加载模型　/mm

Figure 3. Uniaxial loading model

下载: 全尺寸图片幻灯片

试件尺寸 $1\;{\rm{m}} \times 0.5\;{\rm{m}}$ ，采用平面应力方式计算，试件弹性模量取 $E = {\rm{20}}\;{\rm{GPa}}$ ，泊松比取 $\nu = 1/3$ ，密度取 $\rho = 2800\;{\rm{kg/}}{{\rm{m}}^{\rm{3}}}$ ，拉伸荷载取 ${p_0} = 20\;{\rm{MPa}}$ ，以中心为坐标原点，试件在荷载作用下 $x$ 及 $y$ 方向位移的理论解为：

${u_x}\left( {x,y} \right) = \frac{{{p_0}}}{E}x\;\;\;\;$

(64)

${u_y}\left( {x,y} \right) = - \nu \frac{{{p_0}}}{E}y$

(65)

将模型离散为 ${\rm{200}} \times {\rm{100}}$ 节点，节点间距为 $\Delta x{\rm{ = }}0.005\;{\rm{m}}$ ，半径取 $\delta {\rm{ = 3}}{\rm{.015}}\Delta x$ ，不采用任何矫正措施的结果与理论解的误差百分比如图4所示。

图 4 无矫正计算结果相对误差

Figure 4. Relative error of uncorrected calculation results

下载: 全尺寸图片幻灯片

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

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

图 5 应变能密度计算结果

Figure 5. Calculation results by strain energy density

下载: 全尺寸图片幻灯片

图 6 应变能密度相对误差

Figure 6. Relative error by strain energy density

下载: 全尺寸图片幻灯片

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

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

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

图 7 LSM单轴加载结果

Figure 7. Uniaxial calculation results by LSM

下载: 全尺寸图片幻灯片

图 8 LSM单轴相对误差

Figure 8. Uniaxial relative error by LSM

下载: 全尺寸图片幻灯片

表 1 最大相对误差对比表

Table 1. Strain softening stress-strain curve of rock

相对误差	加载方向/(%)	非加载方向/(%)
无矫正	21.33	32.05
应变能密度矫正	1.72	12.73
本文方法	1.69	4.35

下载: 导出CSV

| 显示表格

3.2 双轴荷载

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

图 9 双轴加载模型　/mm

Figure 9. Biaxial loading model

下载: 全尺寸图片幻灯片

两侧拉伸荷载取 ${p_x}= 20\;{\rm{ MPa}}$ ，上下拉伸荷载取 ${p_y}= 15\;{\rm{ MPa}}$ ，位移理论解可用叠加法得：

${u_x}\left( {x,y} \right) = \frac{{{p_x}}}{E}x - \nu \frac{{{p_y}}}{E}x$

(66)

${u_y}\left( {x,y} \right) = \frac{{{p_y}}}{E}y - \nu \frac{{{p_x}}}{E}y$

(67)

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

图 10 双轴应变能密度相对误差

Figure 10. Biaxial relative error by strain energy density

下载: 全尺寸图片幻灯片

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

图 11 LSM双轴加载结果

Figure 11. Biaxial calculation results by LSM

下载: 全尺寸图片幻灯片

图 12 LSM双轴相对误差

Figure 12. Biaxial relative error by LSM

下载: 全尺寸图片幻灯片

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

3.3 裂纹扩展

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

图 13 裂纹扩展模型^[29]　/mm

Figure 13. Crack growth model^[29]

下载: 全尺寸图片幻灯片

参数与文献[29]相同取弹性模量 $E = 72\;{\rm{GPa}}$ ，泊松比 $\nu = {\rm{0}}{\rm{.22}}$ ，密度 $\rho = 2440\;{\rm{kg/}}{{\rm{m}}^{\rm{3}}}$ ，两侧拉伸荷载 ${p_0} = 14\;{\rm{MPa}}$ ，断裂能 $G= 135\;{\rm{ J/}}{{\rm{m}}^{\rm{2}}}$ 。

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

图 14 不同时步损伤分布

Figure 14. Damage distribution at different time steps

下载: 全尺寸图片幻灯片

图 15 裂纹扩展速度对比

Figure 15. Crack growth speed comparison

下载: 全尺寸图片幻灯片

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

4 结论

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

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

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

图 1 物质点相互作用

Figure 1. Interaction of material points

下载: 全尺寸图片幻灯片

图 2 小变形假设

Figure 2. Small deformation assumption

下载: 全尺寸图片幻灯片

图 3 单轴加载模型　/mm

Figure 3. Uniaxial loading model

下载: 全尺寸图片幻灯片

图 4 无矫正计算结果相对误差

Figure 4. Relative error of uncorrected calculation results

下载: 全尺寸图片幻灯片

图 5 应变能密度计算结果

Figure 5. Calculation results by strain energy density

下载: 全尺寸图片幻灯片

图 6 应变能密度相对误差

Figure 6. Relative error by strain energy density

下载: 全尺寸图片幻灯片

图 7 LSM单轴加载结果

Figure 7. Uniaxial calculation results by LSM

下载: 全尺寸图片幻灯片

图 8 LSM单轴相对误差

Figure 8. Uniaxial relative error by LSM

下载: 全尺寸图片幻灯片

图 9 双轴加载模型　/mm

Figure 9. Biaxial loading model

下载: 全尺寸图片幻灯片

图 10 双轴应变能密度相对误差

Figure 10. Biaxial relative error by strain energy density

下载: 全尺寸图片幻灯片

图 11 LSM双轴加载结果

Figure 11. Biaxial calculation results by LSM

下载: 全尺寸图片幻灯片

图 12 LSM双轴相对误差

Figure 12. Biaxial relative error by LSM

下载: 全尺寸图片幻灯片

图 13 裂纹扩展模型^[29]　/mm

Figure 13. Crack growth model^[29]

下载: 全尺寸图片幻灯片

图 14 不同时步损伤分布

Figure 14. Damage distribution at different time steps

下载: 全尺寸图片幻灯片

图 15 裂纹扩展速度对比

Figure 15. Crack growth speed comparison

下载: 全尺寸图片幻灯片

表 1 最大相对误差对比表

Table 1 Strain softening stress-strain curve of rock

相对误差加载方向/(%) 非加载方向/(%)

无矫正 21.33 32.05
应变能密度矫正 1.72 12.73
本文方法 1.69 4.35

下载: 导出CSV

参考文献(29)

[1]	Haeri H, Shahriar K, Marji M F, et al. Experimental and numerical study of crack propagation and coalescence in pre-cracked rock-like disks [J]. International Journal of Rock Mechanics and Mining Sciences, 2014, 67: 20 − 28. doi: 10.1016/j.ijrmms.2014.01.008
[2]	Ma G, Li T, Wang Y, et al. Numerical simulations of nuclide migration in highly fractured rock masses by the unified network method [J]. Computers and Geotechnics, 2019, 111: 261 − 276. doi: 10.1016/j.compgeo.2019.03.024
[3]	赵延林, 万文, 王卫军. 类岩石材料有序多裂纹体单轴压缩破断试验与翼形断裂数值模拟[J]. 岩土工程学报, 2013, 35(11): 2097 − 2109. Zhao Yanlin, Wan Wen, Wang Weijun. Fracture experiments on ordered multi-crack body in rock-like materials under uniaxial compression and numerical simulation of wing cracks [J]. Chinese Journal of Geotechnical Engineering, 2013, 35(11): 2097 − 2109. (in Chinese)
[4]	郑颖人, 王乐, 孔亮. 钢材破坏条件与极限分析法在钢结构中应用探索[J]. 工程力学, 2018, 35(1): 55 − 65. doi: 10.6052/j.issn.1000-4750.2017.05.ST14 Zhen Yingren, Wang Le, Kong Liang. Steel damage condition and application of ultimate analysis method [J]. Engineering Mechanics, 2018, 35(1): 55 − 65. (in Chinese) doi: 10.6052/j.issn.1000-4750.2017.05.ST14
[5]	Kanno H, Moriguchi S. Development of a computational design optimization method for rockfall protection embankments [J]. Engineering Geology, 2021, 284: 105920-1 − 105920-13. doi: 10.1016/J.ENGGEO.2020.105920
[6]	Gao M T, Song Z Q, Duan H Q, et al. Mechanical properties and control rockburst mechanism of coal and rock mass with bursting liability in deep mining [J]. Shock and Vibration, 2020, 9: 1 − 15.
[7]	Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing [J]. International Journal for Numerical Methods in Engineering, 2015, 45(5): 601 − 620.
[8]	Cao W, Wang C. New primal-dual weak Galerkin finite element methods for convection-diffusion problems [J]. Applied Numerical Mathematics, 2021, 162: 171 − 191. doi: 10.1016/j.apnum.2020.12.012
[9]	Belhocine A, Abdullah O I. A thermomechanical model for the analysis of disc brake using the finite element method in frictional contact [J]. Journal of Thermal Stresses, 2020, 43(3): 305 − 320. doi: 10.1080/01495739.2019.1683482
[10]	王璇, 徐明. 胶结型含可燃冰砂土剪切特性的离散元模拟[J]. 工程力学, 2021, 38(2): 44 − 51. doi: 10.6052/j.issn.1000-4750.2020.03.0174 Wang Xuan, Xu Ming. Discrete element simulation of the shear behavior of cemented methane sands [J]. Engineering Mechanics, 2021, 38(2): 44 − 51. (in Chinese) doi: 10.6052/j.issn.1000-4750.2020.03.0174
[11]	Thoeni K, Giacomini A, Lambert C, et al. A 3D discrete element modelling approach for rockfall analysis with drapery systems [J]. International Journal of Rock Mechanics and Mining Sciences, 2014, 68: 107 − 119. doi: 10.1016/j.ijrmms.2014.02.008
[12]	Miotti D, Zamolo R. A fully meshless approach to the numerical simulation of heat conduction problems over arbitrary 3D geometries [J]. Energies, 2021, 14(5): 1351 − 1371. doi: 10.3390/en14051351
[13]	Davydov O, Safarpoor M. A meshless finite difference method for elliptic interface problems based on pivoted QR decomposition [J]. Applied Numerical Mathematics, 2021, 161: 489 − 509. doi: 10.1016/j.apnum.2020.11.018
[14]	Singh I, Mishra B, Bhattacharya S, et al. The numerical simulation of fatigue crack growth using extended finite element method [J]. International Journal of Fatigue, 2012, 36(1): 109 − 119. doi: 10.1016/j.ijfatigue.2011.08.010
[15]	He M, Yang Q, Li N, et al. An extended finite element method for heat transfer with phase change in frozen soil [J]. Soil Mechanics and Foundation Engineering, 2021, 57(6): 497 − 505. doi: 10.1007/s11204-021-09698-z
[16]	Dmitriev A I, Psakhie S G. Molecular-dynamics study of the initial stage of nanoscale deformation localization in the surface layers of a loaded solid [J]. Technical Physics Letters, 2004, 30(7): 578 − 579. doi: 10.1134/1.1783407
[17]	Silling S A. Reformulation of elasticity theory for discontinuities and long-range forces [J]. Journal of Mechanics Physics of Solids, 2000, 48(1): 175 − 209. doi: 10.1016/S0022-5096(99)00029-0
[18]	Silling S A, Epton M, Weckner O, et al. Peridynamic states and constitutive modeling [J]. Journal of Elasticity, 2007, 88(2): 151 − 184. doi: 10.1007/s10659-007-9125-1
[19]	黄丹, 章青, 乔丕忠. 近场动力学方法及其应用[J]. 力学进展, 2010, 40(4): 448 − 459. Huang Dan, Zhang Qing, Qiao Pizhong. A review on peridynamics method and its applications [J]. Advances in Mechanics, 2010, 40(4): 448 − 459. (in Chinese)
[20]	秦洪远, 黄丹, 刘一鸣, 等. 基于改进型近场动力学方法的多裂纹扩展分析[J]. 工程力学, 2017, 34(12): 31 − 38. doi: 10.6052/j.issn.1000-4750.2016.08.0634 Qin Hongyuan, Huang Dan, Liu Yiming, et al. An extended peridynamic approach for analysis of crack growth [J]. Engineering Mechanics, 2017, 34(12): 31 − 38. (in Chinese) doi: 10.6052/j.issn.1000-4750.2016.08.0634
[21]	Zhou X P, Gu X B, Wang Y T. Numerical simulations of propagation, bifurcation and coalescence of cracks in rocks [J]. International Journal of Rock Mechanics and Mining Sciences, 2015, 80: 241 − 254. doi: 10.1016/j.ijrmms.2015.09.006
[22]	王超, 熊伟鹏, 叶礼裕, 等. 冰桨接触过程中的遮蔽效应分析[J]. 华中科技大学学报, 2018, 46(6): 105 − 110. Wang Chao, Xiong Weipeng, Ye Liyu, et al. Analysis of shadowing effect during propeller-ice contact process [J]. Journal of Huazhong University of Science and Technology, 2018, 46(6): 105 − 110. (in Chinese)
[23]	Zhu Q Z, Ni T. Peridynamic formulations enriched with bond rotation effects [J]. International Journal of Engineering Science, 2017, 121: 118 − 129. doi: 10.1016/j.ijengsci.2017.09.004
[24]	王超聪, 刘齐文, 刘立胜, 等. 热防护材料烧蚀温度场的近场动力学模拟[J]. 科学技术与工程, 2017, 17(26): 172 − 176. doi: 10.3969/j.issn.1671-1815.2017.26.027 Wang Chaocong, Liu Qiwen, Liu Lisheng, et al. Numerical simulation of ablation temperature for thermal protective composites based on peridynamics [J]. Science Technology and Engineering, 2017, 17(26): 172 − 176. (in Chinese) doi: 10.3969/j.issn.1671-1815.2017.26.027
[25]	黄小华, 李双, 金艳丽, 等. 冲击荷载作用下泊松比对脆性材料破坏影响的近场动力学分析[J]. 振动与冲击, 2020, 39(20): 204 − 215. Huang Xiaohua, Li Shuang, Jin Yanli, et al. Effect of Poisson’s ratio on the fracture of brittle materials under impact loading via peridynamics [J]. Journal of Vibration and Shock, 2020, 39(20): 204 − 215. (in Chinese)
[26]	Silling S A, Askari E. A meshfree method based on the peridynamic model of solid mechanics [J]. Computers and Structures, 2005, 83(17): 1526 − 1535.
[27]	Prakash N, Seidel G D. Computational electromechanical peridynamics modeling of strain and damage sensing in nanocomposite bonded explosive materials (NCBX) [J]. Engineering Fracture Mechanics, 2017, 177: 180 − 202. doi: 10.1016/j.engfracmech.2017.04.003
[28]	Huang D, Lu G, Wang C, et al. An extended peridynamic approach for deformation and fracture analysis [J]. Engineering Fracture Mechanics, 2015, 141: 196 − 211. doi: 10.1016/j.engfracmech.2015.04.036
[29]	Ha Y D, Bobaru F. Studies of dynamic crack propagation and crack branching with peridynamics [J]. International Journal of Fracture, 2010, 162(1): 229 − 244.

施引文献(4)

期刊类型引用(4)

1.	杜帅奎，尚辰洋，石玥，张成龙，蒋昱，李立州. 双材料复杂界面性能的近场动力学研究. 航空精密制造技术. 2024(04): 13-16 . 百度学术
2.	韩康，张媛，王超，叶礼裕，郭春雨. 提升PD效率的粒子对方法及信息传递接口并行方法. 哈尔滨工程大学学报. 2024(08): 1451-1459 . 百度学术
3.	姚学昊，黄丹. 流固耦合问题的PD-SPH建模与分析. 工程力学. 2022(10): 17-25 . 本站查看
4.	李志远，黄丹，闫康昊. 基于近场动力学微分算子的变截面梁动力特性分析方法. 工程力学. 2022(12): 23-30 . 本站查看

其他类型引用(0)

资源附件(0)

图(15) / 表(1)

计量

文章访问数: 495
HTML全文浏览量: 188
PDF下载量: 122
被引次数: 4

1 理论与假定
1.1 近场动力学理论
1.2 微弹性模型
1.3 小变形假定
2 非局部最小二乘优化
2.1 非局部优化
2.2 最小二乘方法优化
3 数值验证
3.1 单轴荷载
3.2 双轴荷载
3.3 裂纹扩展
4 结论

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

通讯作者: 李树忱(1973−)，男，黑龙江人，教授，博士，博导，主要从事深部岩体破坏机理研究(E-mail: shuchenli@sdu.edu.cn)

计量

出版历程

PERIDYNAMIC METHOD BASED ON NONLOCAL LSM OPTIMIZATION AND DEFORMATION SIMULATION OF BRITTLE MATERIALS

1 理论与假定

1.1 近场动力学理论

1.2 微弹性模型

1.3 小变形假定

2 非局部最小二乘优化

2.1 非局部优化

2.2 最小二乘方法优化

3 数值验证

3.1 单轴荷载

3.2 双轴荷载

3.3 裂纹扩展

4 结论

期刊类型引用(4)

其他类型引用(0)

计量

出版历程

目录

1 理论与假定

1.1 近场动力学理论

1.2 微弹性模型

1.3 小变形假定

2 非局部最小二乘优化

2.1 非局部优化

2.2 最小二乘方法优化

3 数值验证

3.1 单轴荷载

3.2 双轴荷载

3.3 裂纹扩展

4 结论

通讯作者:
李树忱(1973−)，男，黑龙江人，教授，博士，博导，主要从事深部岩体破坏机理研究(E-mail: shuchenli@sdu.edu.cn)