牛彦泽，徐业鹏，黄 丹
(河海大学力学与材料学院，江苏，南京，211100)
摘 要:在双参数微极近场动力学弹脆性模型基础上，引入反映长程力尺寸效应的核函数修正项以提高定量计算精度和收敛稳定性。通过定量变形计算和双轴动载作用下含中心裂纹准脆性板的裂纹扩展过程模拟，验证了本文模型与算法。应用本文模型能准确反映动载作用下的裂纹起裂、扩展、分叉与二次分叉等现象。进一步分析了双轴动载作用下含不同倾角、不同间距平行双裂纹脆性板的破坏模式、起裂时间与承载能力，以及含不同倾角多裂纹脆性板的破坏机制。结果表明:初始裂纹倾角与间距对脆性板承载能力有显著影响，在平行双裂纹倾角一定时，间距越大脆性板承载能力越强；双轴动载作用下多裂纹扩展时，裂纹之间首先相互连接且无分叉，在初始裂纹连接贯通后向边界扩展时出现分叉及二次分叉。
关键词:断裂力学；动力破坏；近场动力学；裂纹扩展；核函数
动载作用下材料和结构的断裂问题及其数值模拟一直是计算力学以及相关工程领域关注的重点和难点。工程实际和实验中常见的裂纹在动载作用下的起裂、扩展、分叉与二次分叉等现象，在数值模拟中却常常难以准确进行描述。尤其是工程中常见的脆性裂纹动态扩展，因其断裂的突然性和危害性，一直是关注的重点。为了分析动载作用下的脆性破坏规律，国内外学者开展了大量试验和数值研究。典型如 Sawamoto等[1]采用离散元法分析了冲击荷载作用下混凝土结构的损伤破坏；刘凯欣等[2]基于离散元法研究了冲击荷载作用下正交各向异性脆性板的破坏及应力波传播；刘海峰和宁建国[3]采用动态本构模型来描述混凝土材料的冲击特性；闫秋实[4]研究了装配式钢筋混凝土梁的抗冲击力学性能；易伟建和刘彪[5]利用偏微分方程求解在冲切作用下圆柱节点承载力的解析解以及方柱节点承载力的级数解。
近年来兴起的近场动力学(Peridynamics，PD)方法基于非局部作用思想建模并通过求解空间积分方程来描述材料的力学行为，以其在解决不连续力学问题方面的独特优势而成为国内外相关学者关注的重点[6－7]。Zhou 等[8]和沈峰等[9]采用 PD 方法分别模拟了含单裂纹、多裂纹岩石板和混凝土构件的单轴破坏过程；王富伟和黄再兴[10]和谷新保和周小平[11]基于PD理论对层合板的冲击损伤破坏和对含圆孔板的破坏问题进行了数值模拟。Wildman等[12]和Zaccariotto等[13]将PD方法与其他方法进行耦合模拟裂纹扩展路径。然而，早期广泛使用的常规微弹脆性 (Prototype Microelastic Brittle，PMB)单参数近场动力学本构模型使用范围受限[14]，其描述的材料泊松比为定值[6－7]，且定量计算精度不高、收敛稳定性差。为此，一方面，Gerstle等[15]构建了双参数PD微极模型(Micropolar Peridynamics)以部分突破常规模型的泊松比限制；另一方面，Huang等[16]、Kilic[17]则在常规近场动力学本构模型中引入反映长程力尺寸效应的核函数修正项，来提高 PD方法的定量计算精度和收敛稳定性。
在已有研究基础上，本文基于PD微极模型引入反映长程力尺寸效应的核函数修正项，构建定量计算精度更高且收敛性更稳定的双参数PD微极模型。通过对方板变形定量计算及双轴动载作用下含中心裂纹脆性板的裂纹扩展过程模拟，验证了本文的模型和算法。进一步分析了双轴动载作用下含不同倾角、不同间距平行双裂纹脆性板的破坏模式、破坏荷载及其受裂纹倾角与裂纹间距的影响规律，以及含不同倾角多裂纹脆性板的破坏机制。
如图1所示，设在某时刻t，某空间域R内任一物质点x与它周围一定范围尺寸δ内的任意物质点x′之间存在相互作用力f:
式中:
δ称为近场范围尺寸；u为位移。
图1 物质点间的相互作用力
Fig.1 Interaction of material points
则对任意物质点，有:
式中:ρ为密度；u˙为加速度；b为体力项。
设近场范围内物质点间的相互作用力与时间t无关，对于无记忆均匀材料，作用力f可简化为[6]:
式中:ξ=x′-x为相对位置矢量；η=u′-u为相对位移矢量。
根据保守势场的性质，存在一个标量函数ω（η,ξ）使得:
则近场域H内的物质点对势能:
基于此，早期广泛使用的微弹脆性近场动力学模型[14](PMB)将物质点对等效成类似于中心弹簧的铰接杆单元，其本构力函数一般形式为:
式中:c为物质点对的拉伸微模量系数(表征其拉伸刚度)；s为物质点对的伸长率:
μ表征物质点对的破坏情况:
s0表示临界伸长率[14]，根据断裂能等效可表示为:
式中:G0为断裂能密度；E为弹性模量。
将局部损伤定义为近场范围内物质点对的断裂数量与其总数之比:
式中，0≤φ≤1表征局部损伤程度。
常规微弹脆性(PMB)模型由于仅考虑物质点对的轴向变形，应用范围受到限制，所能描述的材料泊松比为定值[6,14]。若增加考虑物质点对的弯曲刚度，类似于将原来只考虑轴向拉伸变形的杆单元变换为可考虑轴向和弯曲变形的梁单元，可得到双参数PD微极模型[15]，将本构力函数修正为:
式中，D(ξ)为微弹性矩阵:
式中，d为物质点对的弯曲微模量系数。由此则近场域H内的物质点对势能:
通过能量等效[6,15]，可求出双参数PD微极模型中的微模量系数c和d(以平面应力问题为例):
式中:ν为泊松比。当ν取1/3时，d=0，模型即退化为PMB模型。
在已有研究基础上，本文为了提高双参数 PD微极模型的定量计算精度和收敛稳定性，引入可反映长程力尺寸效应的核函数修正项g(ξ,)δ，则微弹性矩阵D(ξ)变为:
核函数选取满足长程力随物质点对变化规律[16]:
式中，Δ(ξ)为狄拉克函数。
借鉴文献[16]中对常规PMB模型的修正思路，根据式(16)取g(ξ,δ)=[ 1- (ξ/δ)2]2，同样通过物质点对势能与传统理论的应变能密度等效，可得到修正后的双参数(微模量系数):
对于本文研究的动载问题，采用与文献[18]相同的求解思路，将式(2)离散为:
式中:n为时间步长序号；Vj为j处物质点体积，采用中心正交的均匀空间离散，三维物质点体积为Vj=|Δx|3，二维物质点体积为Vj=| Δx|2，|Δx|为离散间距。
在本文模型中，式(18)可表示为:
离散格式采用中心差分格式[19]:
则迭代公式为:
通过初始时刻物质点的位置、边界条件与荷载信息，可以求解任意时间步物质点的位置、速度、位移和合力等。本文基于 Fortran语言编程实现近场动力学迭代计算，后处理采用Paraview实现。
本节通过对经典方板变形问题进行定量计算，验证核函数修正后模型的定量计算精度和收敛稳定性。取如图2所示方板，边长为l=0.4 m ，离散间距 Δx=4 mm ，单向受均匀分布的拉伸荷载，大小为0.25 MPa。材料参数为:弹性模量E=30 GPa，泊松比ν=0.2，密度ρ=2400 kg/m3。为便于定量比较分析，对方板内部点(距左上角点竖直、水平距离均为0.08 m的点1)进行采样观测，改变近场范围尺寸，比较有、无核函数修正项时方板的水平和竖向位移计算结果。
图2 方板模型
Fig.2 Square plate model
图3 所示为有、无核函数作用下观测点水平和竖向位移定量计算结果随近场范围尺寸δ/Δx的变化曲线。由图可知，在无核函数修正时，定量计算结果会随着近场范围尺寸的微小变化而产生突变，且波动很大，收敛稳定性很差。在有核函数作用时，定量计算结果会随着近场范围尺寸的变化连续变化，且趋于稳定，不会产生突变，表明引入核函数修正项的模型计算精度与收敛稳定性更好。
图3 定量计算结果比较
Fig.3 Comparison of quantitative calculation results
文献[20]曾采用基于有限元法的 RFPA (Rock Failure Process Analysis)软件对图4所示的含中心单裂纹脆性板在双轴拉伸动载作用下的破坏过程进行模拟。RFPA软件运用物理统计的方法描述材料的非均匀性，对于某一特定的结构面的非均匀处理，可以采用不同分布参数来定义其性质。它由应力分析和破坏分析两部分组成，应力分析采用线弹性有限元法，破坏判据依据修正的库仑准则。为验证本文的模型和算法，采用与文献[20]相同的模型和工况:方板边长l=254 mm，中心裂纹倾角45 °(本文倾角均指与水平向右方向的夹角)，初始裂纹长度28.3 mm。材料弹性模量E=47.5 GPa ，密度ρ=2500 kg/m3，泊松比ν=0.25，抗拉强度ft=85 MPa。在近场动力学模拟中，方板均匀离散为 64516个物质点，离散间距 Δx=1 mm ，Δt=1.2× 10-8s ，s0取 1.8× 10-3(临界伸长率s0=ft/E)，方板四边施加沿板边均匀分布的动载，采用分级加载方法，加载速率σx=7.583 MPa/s，σy=6.32 MPa/s。
图4 受双轴动载作用含中心单裂纹方板
Fig.4 Square plate with single central crack under biaxial dynamic load
图5 所示为本文模拟得到的裂纹扩展全过程。由图可知，在0.26 ms时裂纹开始起裂，起裂方向沿45 °；随着动载的持续增大，裂纹继续扩展；在0.34 ms时裂纹上下部首次出现分叉；在0.37 ms时，裂纹出现二次分叉；在0.38 ms时，裂纹扩展至边缘，方板完全断裂且裂纹的上下部出现明显的分叉，裂纹基本沿45 °扩展分叉。
图6为本文模拟所得裂纹扩展路径与文献[20]结果的比较。由图可见，二者总体吻合较好，验证了本文模型和算法在模拟动载作用下脆性裂纹扩展问题的可行性。此外，由于采用空间积分方程求解来模拟裂纹扩展，在近场动力学模拟过程中不存在常规连续介质方法模拟时的网格二次剖分或裂尖应力奇异性等问题，随着荷载的施加，初始裂纹自然起裂、扩展和分叉。
图5 中心单裂纹板裂纹扩展
Fig.5 Crack propagation in center-notched plate
图6 中心单裂纹板的破坏路径
Fig.6 Crack paths in a plate with single central crack
双裂纹情形在工程中较为常见，且由于裂纹间的相互影响，破坏模式各异，因而含双裂纹板的破坏过程及其影响因素也颇受关注。图7所示为受双轴动载拉伸作用的含平行双裂纹板，尺寸和材料参数与2.2节一致。两条裂纹长度均为30 mm，加载速率为σx=σy=6 MPa/s 。分别取平行裂纹倾角α为 0°、15°、30°、45°、60°、75°、90°以及两裂纹中点距离d与方板边长l之比d/l为0.2、0.4、0.6的情况进行分析。
图7 受双轴动载作用的含平行双裂纹方板
Fig.7 Square plate with two parallel cracks under biaxial dynamic load
图8 含不同倾角平行双裂纹方板的裂纹扩展
Fig.8 Crack propagation of square plates with two different parallel cracks
图8 所示为d/l=0.4时，含不同初始倾角的平行双裂纹方板在同样双轴动载作用下的裂纹扩展路径与破坏模式。由图可见，动力破坏模式主要有三种:当平行双裂纹初始倾角小于15°时，平行裂纹快速连接汇聚，在扩展过程中不会出现分叉现象；倾角在15°～75°之间，平行双裂纹在扩展和交汇连接过程中会出现明显的分叉、汇集与二次分叉等现象，平行裂纹之间的相互影响非常明显；裂纹倾角大于75°时，平行裂纹之间相互影响较小，裂纹各自扩展并分叉，但初始裂纹之间不交汇。
图9所示为含不同初始间距和倾角的平行双裂纹板受双轴动载作用的起裂时间比较。由于逐渐增大的动载作用下初始裂纹起裂后即迅速失稳扩展，起裂时间可以一定程度上表征结构的承载能力。由图中结果可知，平行双裂纹间距一定时，随着倾角增加，起裂时间先减小后增大，表明脆性板承载能力先减弱后增强，三种情况下的规律均一致。初始倾角 30°时起裂时间最短，即板的承载能力最弱。若平行双裂纹倾角不变时，随着双裂纹间距增加，起裂时间逐渐增大，裂纹之间的相互影响逐渐减小，脆性板承载能力逐渐增强。
图9 含不同间距和倾角平行双裂纹方板起裂时间
Fig.9 Cracking time of square plates with two parallel cracks(with different distances and angles)
考虑到多裂纹之间的相互影响，构建如图10所示的模型，尺寸和材料参数与前文一致，裂纹中心在方板中线上，从左至右四条裂纹倾角分别为 60°、15°、135°和 75°。
图10 受双轴动载作用的多裂纹方板
Fig.10 Square plate with multi-cracks under biaxial dynamic load
图11 所示为含不同倾角多裂纹方板的裂纹扩展与破坏过程。第二条(15°)裂纹尖端在0.288 ms时开始起裂，并向两端扩展；在0.326 ms时，可见另外三条裂纹尖端也各自开始扩展；在0.355 ms时，四条初始裂纹交汇成一条裂纹，并沿着位于边缘的第一条和第四条裂纹的扩展路径继续向板边界扩展，在此扩展过程中出现裂纹分叉和二次分叉；在0.416 ms时，裂纹由多条路径扩展至边缘。
为了进一步分析多裂纹脆性板破坏过程中最先起裂的裂纹与裂纹排布方式的关系，改变上述算例中最先起裂的倾角15°裂纹的位置进行计算。图12所示为改变四条裂纹排布方式的多裂纹方板的破坏结果。结合图11、图12可见，无论怎样改变四条裂纹的排布方式，倾角15°的裂纹尖端均最先开始起裂。由此可见，含不同倾角多裂纹板在动载作用下，最先起裂的裂纹应受其初始倾角影响最大，受裂纹排布顺序影响较小。在多条裂纹之间互相连接贯通过程中，无分叉现象，且有的初始裂纹裂尖出现损伤但不扩展(如图11(c))。在初始裂纹汇聚后向边界扩展时出现分叉及二次分叉等现象。
图11 多裂纹板裂纹扩展过程
Fig.11 Crack propagation in plate with multiple cracks
图12 不同裂纹排布方式的多裂纹板破坏结果
Fig.12 Cack paths of multi-cracked plate with different crack configurations
(1)本文在双参数 PD微极模型基础上，引入反映长程力尺寸效应的核函数修正项以提高定量计算精度和收敛稳定性。通过定量变形计算及双轴动载作用下含中心裂纹准脆性板中的破坏过程模拟，验证了本文的模型和算法。应用本文模型能准确反映双轴动载作用下的脆性裂纹起裂、扩展、分叉与二次分叉等现象。
(2)通过对双轴动载作用下含不同倾角、不同间距平行双裂纹准脆性板的破坏模式、起裂时间与承载能力进行模拟分析，结果表明:双向等加载速率动载作用下，平行双裂纹间距一定时，随着裂纹倾角增加，起裂时间先减小后增大，脆性板承载能力先减弱后增强，在倾角为30°时板承载能力最弱。初始裂纹倾角相同时，板的承载能力随着双裂纹间距增大而增强。
(3)通过模拟双轴等加载速率动载作用下含不同倾角多裂纹脆性板的破坏过程，表明最先起裂的裂纹受其初始倾角影响较大，受裂纹排布顺序影响较小。多裂纹之间先相互连通交汇，再向边缘扩展。在多条裂纹之间互相连接贯通过程中，无分叉现象，且部分初始裂纹裂尖出现损伤但不扩展。在初始裂纹汇聚后向边界扩展时出现分叉及二次分叉等现象。在双轴动载作用下，多裂纹的随机分布情况、倾角等因素对板破坏机制的影响，以及影响多裂纹交汇的关键因素等，有待后续深入研究。
参考文献:
[1]Sawamoto Y, Tsubota H, Kasai Y. Analytical studies on local damage to reinforced concrete structures under impact loading by discrete element method [J]. Nuclear Engineering & Design, 1998, 179(179): 157－177.
[2]刘凯欣, 郑文刚, 高凌天. 脆性材料动态破坏过程的数值模拟[J]. 计算力学学报, 2003, 20(2): 127－132.Liu Kaixin, Zheng Wengang, Gao Lingtian. Numerical simulation for the dynamic failure process in brittle materials [J]. Chinese Journal of Computational Mechanics, 2003, 20(2): 127－132. (in Chinese)
[3]刘海峰, 宁建国. 冲击荷载作用下混凝土动态本构模型的研究[J]. 工程力学, 2008, 25(12): 135－140.Liu Haifeng, Ning Jianguo. Dynamic constitutive model of concrete subjected to impact loading [J]. Engineering Mechanics, 2008, 25(12): 135－140. (in Chinese)
[4]闫秋实, 邵慧芳, 李亮. 冲击荷载作用下装配式钢筋混凝土梁力学性能研究[J]. 工程力学, 2017, 34(4):196－205.Yan Qiushi, Shao Huifang, Li Liang. Study on the behavior of precast reinforced concrete beams under impact loading [J]. Engineering Mechanics, 2017, 34(4):196－205. (in Chinese)
[5]易伟建, 刘彪. 混凝土板柱节点冲切承载力的极限分析[J]. 工程力学, 2017, 34(8): 125－132.Yi Weijian, Liu Biao. Ultimate strength analysis of punching shear capacity of slab-column connections [J].Engineering Mechanics, 2017, 34(8): 125－132. (in Chinese)
[6]Silling S A. Reformulation of elasticity theory for discontinuities and long-range forces [J]. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175－209.
[7]黄丹, 章青, 乔丕忠, 等. 近场动力学方法及其应用[J].力学进展, 2010, 40(4): 448－459.Huang Dan, Zhang Qing, Qiao Pizhong, et al. A review on peridynamics method and its application [J]. Advancein Mechanics, 2010, 40(4): 448－459. (in Chinese)
[8]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 &Mining Sciences, 2015, 80: 241－254.
[9]沈峰, 章青, 黄丹, 等. 基于近场动力学理论的混凝土轴拉破坏过程模拟[J]. 计算力学学报, 2013, 30(增刊1): 79－83.Shen Feng, Zhang Qing, Huang Dan, et al. Damage and failure process of concrete structure under uniaxialtension based on peridynamics modeling [J].Chinese Journal of Computational Mechanics, 2013,30(Suppl1): 79－83. (in Chinese)
[10]王富伟, 黄再兴. 复合材料层合板冲击损伤近场动力学模型与分析[J]. 计算力学学报, 2014, 31(6): 709－713.Wang Fuwei, Huang Zaixing. Peridynamic model and analysis for impact damage of composite laminate [J].Chinese Journal of Computational Mechanics, 2014,31(6): 709－713. (in Chinese)
[11]谷新保, 周小平. 含圆孔拉伸板的近场动力学数值模拟[J]. 固体力学学报, 2015, 36(5): 376－383.Gu Xinbao, Zhou Xiaoping. The numerical simulation of tensile plate with circular hole using peridynamic theory[J]. Chinese Journal of Solid Mechanics, 2015, 36(5):376－383. (in Chinese)
[12]Wildman R A, O’Grady J T, Gazonas G A. A hybrid multiscale finite element/peridynamics method [J].International Journal of Fracture, 2017, 207(2): 1－13.
[13]Zaccariotto M, Tomasi D, Galvanetto U. An enhanced coupling of PD grids to FE meshes [J]. Mechanics Research Communications, 2017, 84: 125―135.
[14]Silling S A, Askari E. A meshfree method based on the peridynamic model of solid mechanics [J]. Computers &Structures, 2005, 83(17/18): 1526－1535.
[15]Gerstle W, Sau N, Sakhavand N. On peridynamic computational simulation of concrete structures [J].American Concrete Institute, 2009, 265(11): 245－264.
[16]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.
[17]Kilic B. Peridynamic theory for progressive failure prediction in homogeneous and heterogeneous materials[D]. Tucson: The University of Arizona, 2008.
[18]郁杨天, 章青, 顾鑫. 含单边缺口混凝土梁冲击破坏的近场动力学模拟[J]. 工程力学, 2016, 33(12): 80－85.Yu Yangtian, Zhang Qing, Gu Xin. Impact failure simulation of a single-edged notched concrete beam based on peridynamics [J]. Engineering Mechanics,2016, 33(12): 80－85. (in Chinese)
[19]Vogler T J, Lee M Y, Grady D E. Static and dynamic compaction of ceramic powders [J]. International Journal of Solids & Structures, 2007, 44(2): 636－658.
[20]Yang Y F, Tang C A, Xia K W. Study on crack curving and branching mechanism in quasi-brittle materials under dynamic biaxial loading [J]. International Journal of Fracture, 2012, 177(1): 53－72.
PERIDYNAMIC MODELLING AND ANALYSIS FOR CRACK PROPAGATION IN BRITTLE MATERIALS SUBJECTED TO BIAXIAL DYNAMIC LOAD
NIU Yan-ze , XU Ye-peng , HUANG Dan
(College of Mechanics and Materials, Hohai University, Nanjing, Jiangsu, 211100, China)
Abstract:A kernel function reflecting the internal length effect of long-range forces has been introduced into a regular micropolar prototype elastic brittle peridynamic model, in order to improve the calculation accuracy and numerical stability of convergence. The proposed model and algorithms were validated through deformation calculation and simulation of crack propagation in a center-notched quasi-brittle plate subjected to a biaxial dynamic load. The crack initiation, propagation and branching can be described naturally by using the proposed approach. The failure mechanism as well as the influence of the direction and distance between two parallel cracks on the failure mode，cracking time, and the extreme failure load of a brittle plate subjected to biaxial dynamic loads have been analyzed further. Numerical results show that the directions of pre-existed cracks as well as the distance between cracks take an important effect on the failure mechanism of the plate. With given crack orientation, the plate shows higher loading capacity when the distance between two parallel pre-existed cracks increases. In the multiple-crack propagation process under biaxial dynamic loads, the pre-existed cracks connect to each other firstly without bifurcation, and then the converged cracks will bifurcate and re-branch whenpropagate to the boundary.
Key words:fracture mechanics; dynamic failure; peridynamics; crack propagation; kernel function
徐业鹏(1982―)，男，湖北人，讲师，博士，主要从事断裂力学研究(E-mail: xuyepeng285@sina.com).
牛彦泽(1993―)，男，黑龙江人，硕士生，主要从事近场动力学研究(E-mail: hhuyanze_niu@foxmail.com);
作者简介:
通讯作者:黄 丹(1979―)，男，湖北人，教授，博士，主要从事计算力学与工程仿真研究(E-mail: danhuang@hhu.edu.cn).
基金项目:国家自然科学基金项目(51679077, 11132003)；中央高校基本科研业务费项目(2015B18314，2017B13014)；江苏省自然科学基金项目(BK20130822)；NSFC-广东联合基金(第二期)超级计算科学应用研究专项项目(U1501501)
收稿日期:2017-06-09；修改日期:2017-11-16
文章编号:1000-4750(2018)10-0249-08
doi:10.6052/j.issn.1000-4750.2017.06.0450
文献标志码:A
中图分类号:O346.1+1