改进键型近场动力学方法下的多裂纹板破坏分析

2020-05-25徐业鹏

科学技术与工程 2020年10期

孙杰, 徐业鹏

(河海大学力学与材料学院，南京 211100)

裂缝的存在通常会削弱工程结构的承载能力，甚至导致结构失效，因此含缺陷材料和结构的断裂问题及其数值模拟一直都是计算力学的研究重点。近年来，近场动力学方法以其在模拟裂纹扩展以及分析结构断裂方面的独特优势日益受到国内外学者的青睐。

近场动力学(peridynamics，PD)是基于非局部作用思想建立模型并通过求解空间积分方程来描述物质力学行为的一种方法。该方法由Silling[1]于2000年提出，随后，为证明求解不存在裂尖应力奇异性，Silling等[2]又对PD基于键的理论进行更深一步的推导，成功应用于研究微弹性材料[3-4]和混凝土材料[5]中的复杂裂纹问题。常规微弹性脆性模型(prototype microelastic brittle,PMB)在应用时存在限制，为此Gerstle等[6]在开展多晶断裂和PD断裂模拟的同时，引入了弯矩密度概念来解除泊松比的限制，引入了双参数微极模型。Huang等[7]、顾鑫等[8]为提高计算精度，又引入核函数修正项。黄丹等[9-10]进一步分析改进了核函数修正项对PD的作用，并开发出相应的静力动力的算法。钱剑等[11]在此基础上进行了动载作用下的复合型裂纹扩展模拟。

与此同时，多裂纹扩展问题中裂纹间的相互影响和连接一直都是较为复杂的难题，众多学者针对多裂纹问题进行过大量研究。Silling[12]用冲击物撞击了带有平行裂纹的圆盘，利用近场动力学成功模拟此算例，在试验中观察到了裂纹扩展的角度。Zhou等[13]研究了裂纹阵列对多裂纹扩展和聚结过程的影响。Zeng等[14]用扩展有限元方法，观察脆性岩石类材料的裂纹对其弹性性能和强度的影响情况。Vazic等[15]研究了小裂纹对宏观裂纹动力学扩展的影响，考虑了不同数量、不同位置和密度的小裂纹的各种组合，其大小取决于小裂纹的位置、密度和数量。张振南等[16]基于VMIB方法，通过对两条不同排列的平行裂纹开裂过程的模拟，探究多裂纹之间的相互影响。基于改进后的近场动力学方法，对含双裂纹的脆性板进行单轴拉伸破坏模拟，并将所得结果与已有文献结果进行对比，验证本文提出的本构模型和算法的可靠性。随后通过模拟不同初始形态的多裂纹板的裂纹扩展路径和分析承载能力，探究裂纹初始状态对结构破坏型式、扩展路径和承载能力的影响规律。

1 近场动力学基本理论

1.1 键型PD基本思想

近场动力学的基本思想如下：在假设物质点系统做刚体运动，保持物质系统构型不变并且物质点对之间的作用与时间无关的前提下，考虑某一个空间R内的所有物质点，领域内任意物质点x和一定范围大小δ(近场范围)内的其他物质点x′∈H:{|x′-x|≤δ,H∈R}存在关于位移u=u(x,t)的相互作用力f(图1)，那么对于该物质点就有如下关系：

f=f[x,x′,u(x,t),u(x′,t),t]

(1)

(2)

式中：H为近场范围；ρ为材料密度；b为外力密度。

图1 物质点对之间的相互作用

对于均质材料，点对力函数可以继续简化为与参考构型中两物质点的相对位置ξ=x′-x和相对位移η=u′-u相关的函数：

f=f(x′-x,u′-u)=f(ξ,η)

(3)

均匀材料的本构力函数的一般形式为

f(η,ξ)=F(ξ)(η+ξ)

(4)

根据Silling[1]以及黄丹等[17],本构力(点对力)可表示为

(5)

式(5)中：μ判断物质点对破坏情况。

(6)

式(6)中：s0为物质点对的临界伸长率，可取s0=[4πG0/(9Eδ)]1/2。

物质点对势能与传统应变能等效[1]求解微模量系数c，结果如下：

平面应力问题：

(7)

平面应变问题：

(8)

式中：E为弹性模量；ν为泊松比。

常规PMB模型在平面应力问题求解时受到泊松比限制，为突破该限制，Gerstle等[6]引入物质点对的转动自由度，提出双参数微极模型：

f(η,ξ)=D(ξ)η

(9)

令ξ=|ξ|，其中，

(10)

通过物质点对势能与传统应变能等效[10]，求出双参数微极模型的参数：

(11)

为了提高计算精度，加入考虑了长程力的尺寸效应，加入核函数修正项g=[1-(ξ/δ)2]2。将原来的常规微弹性脆性模型系数修正为c(ξ)=c(0,δ)g(ξ)。

核函数修正项必须满足长程力的变化规律：

(12)

引入核函数修正项g=[1-(ξ/δ)2]2后得到修正的微模量系数:

(13)

1.2 求解体系

在任意时刻，x由近场范围内物质点对作用力产生的内部合力可以表示为

∀x∈R,t≥0

(14)

PD是将固体离散为一系列具有信息的物质点，设物质点间距为|Δx|，则将式(2)改写为

(15)

通常在求解时加入人工阻尼来保证获得静力解，基本方程如下：

(16)

采用中心差分格式

(17)

按时间顺序进行离散,代入(15)化简，得到：

(18)

由Fortran程序完成求解。

2 算例分析

2.1 模型验证

为验证本文本构模型和算法的正确性，首先对含双裂纹混凝土板进行了单轴拉伸破坏模拟。具体模型如图2所示，弹性模量E=30 GPa，泊松比ν=0.33，材料密度为ρ=2 400 kg/m3，离散为10 200个物质点，离散间距为0.000 5 m,近场尺寸δ=3Δx=0.001 5 m。在模拟过程中采用位移加载控制方式对其进行加载，每一步荷载增量为0.04 mm，迭代步长为Δt=1×10-7s。

为模拟不同裂纹初始布置下板的断裂破坏模式，采用如下三种方案(均保持两条裂纹平行)。

(1)方案I：2b/2a=1，即b=5 mm，α=45°。

(2)方案II：2b/2a=2，即b=10 mm，α=45°。

(3)方案III：2b/2a=1，即b=5 mm，α=60°。

模拟结果如图3(方案I～III的裂纹扩展过程)所示，裂纹沿着最易汇合的方向扩展，一起沿水平方向开裂，在两条裂纹之间形成破裂区，进而贯通，同时两条裂纹也向板的两侧扩展，最后形成一条贯通性裂纹。三个方案中PD所得裂纹扩展路径与VMIB[16]方法所得裂纹扩展路径基本一致。

图2 双裂纹板模型

图3 裂纹扩展路径比对

2.2 含三裂纹脆性板单轴拉伸破坏分析

主要研究在含三裂纹脆性板在受到单轴拉伸荷载作用时，裂纹的初始布置是否会改变裂纹扩展的路径以及对临界破坏荷载是否有影响。本模型具体材料参数和2.1节一致，q=1 000 kPa，迭代时间步长Δt=1×10-7s。模型如图4所示。

图4 三裂纹板模型

观察裂纹的扩展路径主要有两种形式，即以中心裂纹为主的开裂与以两边裂纹为主的开裂，最后都是沿着竖直方向扩展直到板完全破坏。以α=30°为例，当β=0°时，脆性板以两边裂纹为主开裂，直至破坏，当β逐渐增大后，裂纹扩展路径开始发生变化，由以周边裂纹为主的开裂变为以中心裂纹为主的开裂，如图5(a)所示。

观察有无周边裂纹时临界破坏荷载，如“β=0°且无周边裂纹”与“β=α=0°”对比可得，周边裂纹的出现明显削弱了板的承载能力，但同时对比“β=15°且无周边裂纹”与“α=0°，β=15°”，板承载能力又有稍许提升，可见裂纹的增加对板的承载能力存在复杂的影响。进一步观察图5(b)可知，临界破坏荷载随着中心裂纹角度与周边裂纹角度的增加而减小。当三条裂纹角度都在0°时，板的承载能力最强，即此时板处在最稳定的状态，而当有一条裂纹达到90°时，临界破坏荷载达到最小值。在周边裂纹角度值一定的情况下，板的承载能力随着中心裂纹角度的增加而减小。同样，在中心裂纹角度一定的情况下。临界破坏荷载值也随着α的增大而减小。即板整体承载能力随着中心裂纹角度与周边裂纹角度的增加而减小。

图5 PD模拟结果

2.3 含多裂纹脆性板单轴拉伸破坏分析

含初始裂纹板的破坏问题[18-20]在实际工程与航空航天领域都有重要的作用，在文献[21]中，由于航空事故对航空问题中的铝板进行了裂纹破坏问题分析，而该类问题在脆性板中也具有深远研究意义。采用的模型如图6所示，由于初始裂纹布置对板的承载能力有着不可忽略的影响，通过改变角度α与中心裂纹角度β，观察初始裂纹布置对临界破坏荷载的影响。材料数据如下：弹性模量E=0.69 GPa，泊松比ν=0.3，材料密度为ρ=2 400 kg/m3，取物质点间距为0.001 m，物质点总数为30 200，近场尺寸δ=3Δx=0.003 m。对其进行均布拉伸，q=1 000 kPa，迭代步长为Δt=1×10-7s。具体模型如图6所示。

图6 多裂纹板模型

由图7易知，在α与β同时达到90°时，临界破坏荷载达到最大值，几乎接近板的抗拉强度。而在α=β=0°时，临界破坏荷载最小，即板处在最脆弱的状态。

在角度β不变的情况下，临界破坏荷载随着周边裂纹角度α的增加而增加，同样，在α不变的情况下，临界破坏荷载也几乎随着β的增加而增加，仅在β=45°时，α由15°～30°时临界破坏荷载出现了较为明显的下降，整体规律为随着α与β的增大，临界破坏荷载越大，即板的承载能力越强，越稳定。反之当α与β的越小，临界破坏荷载越小，即板越易破坏。

此外，在不同阶段，临界破坏荷载随裂纹角度变化幅度也有差异。在中心裂纹角度β<45°的时候，板的承载能力受周边裂纹角度α影响较小，例如在β=15°时，图7中折线较为平稳，即临界破坏荷载随角度α变化波动不大，当β>45°的情况，整体折线相对较陡，即临界破坏荷载随角度α变化波动较为明显。