平行双裂纹扩展的单键双参数近场动力学研究

黄小华， 李双， 金艳丽， 何小桥， 王家明

(1.广西大学 土木建筑工程学院，广西 南宁 530004；2. 工程防灾与结构安全教育部重点实验室，广西 南宁 530004；3. 香港城市大学 建筑学及土木工程学系，香港 999077)

材料在生产和使用过程中会产生一系列细小裂纹而带缝工作，这些裂缝极有可能导致材料的破坏，甚至结构的失效[1]。因此在材料中预制裂纹，研究裂纹扩展规律是解决含缺陷材料破坏问题的一种较为常见的手段。有限元等基于传统连续介质力学的数值方法在结构变形分析中取得了巨大的成功，但是它在不连续处没有定义，存在奇异性，在解决材料的破坏问题上举步维艰。新近兴起的近场动力学理论[2]采用位移的空间积分形式完全重构了传统连续介质力学理论，这种积分型方程在不连续处仍有定义，且形式不变，实现了对连续介质和非连续介质力学行为的统一描述，能模拟材料中裂纹自发地萌生和扩展，是一种特别适合模拟材料破坏的力学理论和方法。Silling[2]提出了键基近场动力学理论，它将物质点间的相互作用视为它们中心受力的键，这种中心力模型只有一个参数，即键刚度(或微观模量)，称为传统键基近场动力学模型(或单参数键基近场动力学模型，bond-based peridynamics，BPD)。

均匀各向同性线弹性材料的力学属性须由2个独立的弹性参数才能完全表述，采用单参数键基PD模型描述其力学行为时，2个独立弹性参数将被缩减为1个参数，必然导致某一个弹性参数将以恒定的形式固化到中心力模型中。例如泊松比参数被固化成常数，在平面应力状态下，泊松比只能是1/3，平面应变状态和三维应力状态下，泊松比只能是1/4[3]。毫无疑问，这种泊松比恒定的要求大大限制了传统键基近场动力学理论的适用范围。为此，Silling[4]又提出了态基近场动力学理论，它是通过考虑两物质点邻域内所有键的集合作用来增强材料的本构描述。与键基近场动力学模型中2物质点间的作用仅由它俩自身决定不同，态基近场动力学模型中2物质点间的作用须考虑2物质点邻域内所有物质点[5]，是一种典型的多体相互作用。因此，态基近场动力学方法的计算复杂性以及对资源的耗费都较键基近场动力学有了更高的要求。于是，为了保留传统键基近场动力学模型的简单性以及数值实施的稳定性，部分学者仍立足于键基近场动力学模型的两体间相互作用，通过引入切向键[6-8]、引入键的转角[9]、力量补偿方案[10]或以梁换键[11-12]等方式扩展键基近场动力学理论，突破了泊松比的限制。其中，通过引入切向键或键的转角等方式增加2物质点间的切向作用较为便捷，对传统键基近场动力学模型改动很小。Prakash等[9]通过引入法向弹簧，提出一种二维双参数键基近场动力学模型。Zhou等[6-7]基于超弹性理论提出一种三维两参数共轭键基近场动力学模型。Zhu[8]等通过增加键的转角，提出双参数键基近场动力学模型。当泊松比取上述恒定值时(如在平面应力状态下，泊松比取1/3)，两参数都能退化为单参数。但一些模型[6-7,9]退化后的参数与传统单参数键基近场动力学的参数并不相同。

本文在BPD模型的基础上，通过键内直接添加切向刚度系数，提出一种单键双参数键基近场动力学模型(TBPD)，推导模型的双参数与材料弹性参数之间的转换关系，建立TBPD模型单键断裂条件。然后探讨不同泊松比下BPD模型和TBPD模型的适用性，最后，采用TBPD模型，研究裂纹间距和随机分布的细小裂纹对平行双裂纹薄板破坏的影响。

1 单键双参数键基近场动力学方法

1.1 理论模型

传统键基近场动力学模型认为宏观连续体在其所在空间域Ω内由大量物质点组成，物质点以初始构型的位置X作为标记，携有体积VX和质量密度ρ。任一物质点X仅在其有限的“域”内通过键与其他物质点X′存在相互作用力f。域是以该物质点为中心，以δ为半径的范围，记域内所有其他物质点为Hδ={X′∈Ω:|X′-X|≤δ}，据牛顿第二定律，物质点X在t时刻的运动方程为：

(1)

传统键基近场动力学模型将物质点间的相互作用视为它们中心受力的键，键仅有一个刚度系数，表示法向作用，故这里称之为单键单参数键基近场动力学模型。

在传统键基近场动力学模型的基础上，通过键内直接添加切向刚度系数，使其具有法向和切向双刚度，并称之为单键双参数键基近场动力学模型，如图1所示。

图1 TBPD模型中键的运动示意Fig.1 Motion of a bond in TBPD model

(2)

此时，X和X′间键的微观弹性应变密度(单位体积的应变能密度)为:

(3)

将每个键的应变能都等分给键两端的物质点，则物质点X的宏观弹性应变能密度为：

(4)

对点力f是键的微观弹性应变能密度ω对该键相对位移矢量η的导数：

(5)

以示区别，传统BPD模型的对点力[14]为：

(6)

(7)

则式(6)可改写为：

(8)

(9)

式中δij为Kronecker符号。于是，Green应变张量E可表示为：

(10)

小变形时有|∂ui/∂Xj|≪1，略去位移梯度的乘积项，Green应变E变为Cauchy应变ε，则式(9)变为：

Fij=δij+εij+ωij

(11)

(12)

将方程(3)代入方程(4)，得物质点X的宏观弹性应变能密度：

(13)

将方程(12)代入方程(13)化简得：

(14)

式中h为材料厚度。在经典连续介质力学中，平面应力状态的应变能密度为：

(15)

(16)

式(16)中的法向刚度参数cn与传统键基近场动力学的微模量[14]完全相同。若方程(14)和方程(15)中采用平面应变状态下或者三维状态下的应变能密度，则可得到对应应力状态下TBPD模型的微模量与材料弹性参数间的转换关系，如表1所示。

表1 TBPD模型的微模量cn和ctTable 1 Micro-modulus cn and ct of TBPD model

不难验证，当切向刚度系数ct=0时，TBPD模型退化后参数cn与BPD模型参数完全相同，且有方程(5)右边表达式将退化成方程(8)右边表达式，则此时，TBPD模型将近似为BPD模型。微模量参数cn和ct作为键基PD的刚度系数，必须满足cn≥0且ct≥0。因此TBPD模型的泊松比适用范围为：平面应变和三维情形时-1<ν≤1/4，平面应力时-1<ν≤1/3。

值得指出的是，表1转换关系与Zhu等[8]相同，但两者引入切向刚度的方法和推导过程并不同。TBPD模型通过键内直接增加切向刚度系数，在小变形条件下，以均匀应变工况，得到模型的双参数与材料弹性参数的转换公式。Zhu等[8]未就双参数键基PD模型开展破坏分析，本文接下来将讨论TBPD模型的断裂准则及其破坏应用。

1.2 单键断裂条件

这里采用能量准则作为判断单键断裂的依据，即键中包含的微观应变能密度ω(η,ξ)超过某个极限值ωc(η,ξ)(简记为ωc)时，认为该键断裂，将不再传递力的作用。

考虑含裂纹CD的无限大板，板厚h，裂纹长度为2a。当裂纹CD穿过某一物质点的邻域时，所有被裂纹穿过的键都将断裂，存储在这些键中的应变能也被耗散。如图2所示，键AB将发生断裂，存储在该键中的应变能ωcVAVB将被耗散，其中VA和VB分别为物质点A和物质点B携带的体积，当VA和VB趋于无限小时，所有因裂纹CD而耗散的键的总应变能用积分表示为：

图2 断裂面未完全穿过物质点邻域的微势能积分域Fig.2 Integration domain of the micropotentials not completely crossing a fracture surface

(17)

式中：ΩA和ΩB分别表示位于裂纹CD两侧的物质点A和物质点B的取值范围，它们须满足一定的断键条件，即线段AB的距离rAB≤δ，且线段AB和线段CD有交点。

ΩB1+ΩB2=ΩB

(18)

则对于所有距裂纹CD垂直距离为z的物质点，因裂纹CD而耗散的键的总应变能为：

(19)

式中：dz为物质点A携带体积在z轴的长度，S′A为裂纹CD在直线AA2所在水平面上的投影区。

记Gc为产生每单位裂纹表面所需的能量，则：

(20)

与传统BPD模型相类似，将ωc假定为ξ的线性函数，即:

ωc(η,ξ)=ωc0ξ

(21)

将式(21)代入式(20)，得:

(22)

故

(23)

由式(3)、式(21)和式(23)知，TBPD模型的单键断裂条件为：

(24)

在采用TBPD方法计算过程中，一旦键满足式(24)，则该键断裂且永久失效。不难验证，当切向刚度系数ct=0时，TBPD模型的单键断裂条件可退化后为BPD模型的单键断裂条件，即：

(25)

式中c和sc分别为BPD模型中键的微观模量以及键的临界伸长率。

在TBPD模型中，键断裂与否采用状态函数μ(ξ,t)来记录。物质点的损伤通过物质点邻域内键的断裂程度，用函数φ(X,t)来计算，μ(ξ,t)和φ(X,t)具体的表达式可参考文献[15]。

2 模型验证

算例1尺寸为1 m×0.5 m的长方形薄板，板厚0.01 m，其弹性模量200 GPa，密度7 850 kg/m3，沿两短边作用p=200 MPa的均匀法向拉力，采用BPD和TBPD模型分别计算，物质点间距Δ=0.01 m，邻域半径δ=3.015Δ，拉伸荷载通过转变成体力荷载p/Δ，施加在两短边最外一层的真实物质点上[14]，如图3所示。

图3 单轴拉伸荷载作用下的薄板近场动力学离散示意Fig.3 Geometry of a plate under uniaxial tension and its PD discretization

1)当泊松比ν=1/3时，如图4所示，采用BPD模型或TBPD模型计算的位移结果都与解析解很好地吻合。采用TBPD方法计算的位移ux的最大相对误差为3.25%，分布在薄板4个角点附近，这是由近场动力学的边界效应导致的。

图4 2种模型计算的位移ux的相对误差Fig.4 Relative error on ux of two models

2)当泊松比ν∈[0,1/3]时，采用反算泊松比方法[8-9]和计算物质点的位移ux的相对误差来验证BPD和TBPD模型的合理性，结果列于表2和表3。

表2 BPD模型与TBPD模型的泊松比ν反算结果 Table 2 Poisson′s ratio inverse calculation result of BPD and TBPD model

表3 BPD模型与TBPD模型计算的ux最大相对误差Table 3 Maximum relative error on ux calculated by BPD model and TBPD model %

由表2可知，BPD模型无论材料泊松比怎么取值，反算泊松比总是0.332 8，当ν≠1/3时相对误差大，这证实了它仅适合用来描述前面所说的特定泊松比材料。而采用TBPD模型时，反算泊松比值与真实泊松比非常接近，相对误差为0.16%～0.70%。故它较适合用来描述泊松比在一定范围内的材料力学行为。此外，采用TBPD模型计算，在ν∈[0,1/3]时，位移ux的误差全部小于5.0%，且随着泊松比的减小，误差逐渐减小；而采用BPD模型计算，在ν∈[0,0.20]时，位移ux的误差全部大于5.0%，且泊松比越小，相对误差越大。

图5 剪切边界条件几何模型示意Fig.5 Geometric model of shear boundary conditions

如图6所示，采用TBPD模型计算的位移结果与解析解高度吻合，位移ux和uy的最大相对误差都为0.202%。

图6 TBPD方法计算的位移ux与解析解间的相对误差Fig.6 Relative error on ux for TBPD method

算例3边长为50 mm的正方形薄板，中心有一长度2a=10 mm的水平裂纹，弹性模量192 GPa，密度8 000 kg/m3，泊松比ν=1/3，临界断裂能83 kJ/m2，在薄板上、下边界作用大小20 m/s的均匀法向速度荷载。分别采用BPD和TBPD模型进行计算，物质点间距Δ=0.25 mm，邻域半径δ=3.015Δ，计算时步取1.3367E-8 s，上、下两边的速度边界条件是施加在薄板上下边以外宽度为δ的3层虚拟物质点上[14]。

采用BPD和TBPD模型计算时，薄板都是在裂纹左右两端同时启裂，然后各自向左右两侧沿水平方向扩展，最终贯通薄板。薄板在启裂和贯通时刻稍有差异，但薄板破坏耗时基本相同，相关数据列于表4中。如图7所示，取t=6.68 μs时刻的位移进行对比，两模型计算的位移uy差异甚微。如图8所示，取t=24.09 μs时刻的薄板破坏形态进行对比，两模型计算的结果几乎相同。

图7 泊松比ν=1/3时2种模型在t=6.68 μs时刻的uy云图Fig.7 Comparison of the uyat ν=1/3, t=6.68 μs of two models

图8 2种模型在t=24.09 μs时刻的破坏形式对比Fig.8 Comparison of the failure modes at t=24.32 μs of two models

表4 薄板启裂和贯通时刻Table 4 Initiating and through time of plate μs

算例4边长200 mm的正方形薄板，其弹性模量为30 GPa，密度2 265 kg/m3，Ⅰ型能量释放率GΙc=110 J/m2，Ⅱ型能量释放率GIIc=10GIc，泊松比0.2。在试件中部的左右两侧各有一条长25 mm，宽5 mm的预制裂纹[16]， 如图9所示。首先在薄板左侧上半部分和右侧下部分都水平施加合力为5 kN的均布力p，然后在薄板上下边界施加v=10 mm/s的速度荷载。采用TBPD模型计算，物质点间距Δ=4 mm，邻域半径δ=3.015Δ，计算步长1.0×10-7s，左右两侧的均布力p通过转变成体力荷载p/Δ后，施加在薄板该部位的最外一层，宽度为Δ的真实物质点上，速度荷载v作用在薄板以外宽度为δ的3层虚拟物质点上。

图9 混合破坏模式断裂试样的几何形状Fig.9 The geometry of the mixed mode fracture specimen

图10(a)给出了采用TBPD模型计算得到的薄板的最终破坏形态，它与OŽBOLT等采用基于微平面模型的有限元代码模拟的结果[16]和试验结果[17]基本相同，如图10(b)、(c)所示。

图10 混合破坏模式断裂试样的模拟结果Fig.10 Simulation results of the mixed mode fracture specimen

上述算例分别从薄板连续变形的定量计算和薄板破坏形态模拟2个方面都验证了TBPD模型的合理性和有效性。

3 平行双裂纹薄板的受拉破坏

如图11所示，一边长为50 mm的正方形薄板，中部预制有长2a=10 mm的2条水平裂纹，记为主裂纹①和主裂纹②，两裂纹纵向间距为d。薄板的弹性模量为203 GPa，泊松比0.3，密度为7 850 kg/m3，在薄板上、下边界作用大小20 m/s的均匀法向速度荷载。所用参数与赵金海等[18]相同，采用TBPD模型进行计算，物质点间距Δ=0.5 mm，邻域半径δ=3.015Δ，计算步长取1.336 7×10-8s。按以下2种情形分别计算和讨论：

图11 平行双裂纹薄板的裂纹布置Fig.11 Micro-crack distribution of plates with parallel double cracks

1)平行双裂纹间距d取3 mm、5mm和7 mm时，对薄板破坏形态的影响；

2)平行双裂纹间距d=3 mm时，考虑3种随机分布的细小裂纹对薄板破坏形态的影响。

3.1 裂纹间距的影响

d取3 mm、5mm和7 mm时，裂纹扩展路径如图12(b)、(c)、(d)所示。薄板破坏过程都是在裂纹左右两端同时启裂，水平延伸一定距离后，裂纹①斜向上方、裂纹②斜向下方继续扩展，两裂纹纵向间距越来越大，直至贯通薄板。裂纹间距对薄板破坏模式影响不大，但对于裂纹水平延伸距离、扩展角度、裂纹贯通后纵向间距、起裂和贯通时间都有一定影响，如表5所示。随着裂纹间距的增大，裂纹起裂时刻提前，水平延伸长度逐渐减小，但裂纹偏移角度、贯通间距和贯通时刻却随之增大。

图12 TBPD模型计算的不同间距d下平行双裂纹扩展Fig.12 Crack propagation paths of parallel double cracks calculated by TBPD model with longitudinal spacing d

表5 不同纵向间距d的平行双裂纹扩展情况Table 5 Propagation of parallel double cracks with different longitudinal spacing d

值得注意的是，图12(b)破坏模式与赵金海等基于传统BPD模型所得结果[18]有很大差异，文献[18]给出的破坏模式如图13所示，其双裂纹水平延伸一定距离后，都斜向薄板水平中心交汇成一条裂纹，然后水平延伸，直至贯通薄板。鉴于传统BPD模型只能模拟ν=1/3的材料的力学行为[3]，所以将它用于模拟其他泊松比材料时，裂纹扩展路径将可能出现失真或失效的情形。

图13 d=3 mm时BPD模型计算泊松比0.3的薄板破坏Fig.13 Crack propagation path of plate with Poisson′s ratio 0.3 calculated by BPD model at d=3 mm

3.2 3种随机细小裂纹的影响

在平行双裂纹薄板内，裂纹间距d=3 mm时考虑3种随机细小裂纹的分布，细小裂纹长度取0～2 mm，细小裂纹布置方式和破坏形态见图14。

如图14(a)～(d)所示，当细小裂纹主要分布在某一主裂纹附近时，该主裂纹将沿着垂直于加载方向扩展直至贯通，而另一主裂纹几乎没有扩展。如图14(e)、(f)所示，当细小裂纹分布在主裂纹的某一端附近，则主裂纹沿着分布有细小裂纹的这一端在荷载作用下扩展贯通，而该主裂纹的另一端只发生较小变化。对比图12(b)可以看出，随机细小裂纹的分布在一定程度上会抑制某一主裂纹在单侧或两侧的扩展，导致薄板的破坏形态发生显著的改变。

图14 含随机细小裂纹的平行双裂纹板及其扩展路径Fig.14 Parallel double cracks plate with small cracks and its′ propagation paths

4 结论

1)在不计刚体转动和小变形条件下，以均匀应变工况，一次得到了模型的双参数与材料弹性参数的转换公式。且切向刚度系数的引入，使双参数键的性能较单参数键更加完备，解决问题的能力更强。BPD模型只能模拟特定泊松比材料的力学行为，否则材料的反算泊松比、位移结果存在较大误差，破坏模式也可能失真；TBPD模型突破了特定泊松比的限制，能够模拟泊松比在一定范围内的材料力学行为，退化后的模型参数与传统BPD的参数完全相同。

2) 以能量准则作为判断单键断裂的依据，详细给出了TBPD模型单键断裂条件，在特定泊松比情形下，它可退化为BPD模型的单键断裂条件。

3)在一定条件下，平行双裂纹的纵向间距对薄板的破坏过程产生一定的影响，随着双裂纹纵向间距的增大，裂纹起裂时刻提前，水平延伸长度逐渐减小，但裂纹偏移角度、贯通间距和贯通时刻却随之增大。随机细小裂纹的分布在一定程度上会抑制某一主裂纹在单侧或两侧的扩展，导致薄板的破坏形态发生显著的改变。