王智环 贾雷明 何增 田宙

1)(清华大学工程物理系,北京 100084)

2)(西北核技术研究所,西安 710024)

基于线性硬化塑性本构模型,建立了冲击载荷作用下弹塑性球面应力波场的理论求解方法.首先,分析了冲击载荷卸载速率对球面应力波传播的影响,得到了3 种不同类型的应力波传播图像.在此基础上,建立了弹性阶段、塑性加载阶段以及卸载阶段球面波动方程的理论求解方法,给出了质点位移、质点速度、应力和应变等物理量的计算方案.与已有理论方法相比,该方法考虑了不同载荷卸载速率条件下应力波的不同传播情况,并且给出了卸载阶段应力波参量计算方法,具有更广的适用范围.利用上述方法计算了恒定冲击载荷和不同指数衰减冲击载荷作用下弹塑性球面应力波场参量,在弹性阶段和塑性加载阶段,理论计算得到的物理量与已有理论方法和数值模拟结果基本吻合,在卸载阶段,已有理论方法不再适用,而本文理论计算得到的物理量与数值模拟结果基本吻合,验证了该理论方法的正确性.

1 引言

无限大固体介质中球面应力波的传播问题一直被关注和研究.基于弹性[1-9]和黏弹性[10-15]本构关系,大量理论研究已经开展,分析了球面应力波传播特征和规律.当载荷较高时,介质发生塑性变形,需要利用弹塑性本构模型对介质进行描述,因此,开展弹塑性介质中球面应力波传播理论研究具有重要的意义.

针对弹塑性球面应力波,多位学者已经开展了部分理论研究.Yang[16]利用特征线理论研究了线性硬化塑性材料在线性衰减的冲击载荷作用下球面应力波的传播情况,分析了塑性剪切模量的影响,但对卸载区的求解借助了数值方法.Morland[17]通过假设跨越塑性波波阵面的应力间断的具体表达形式,限定了外加载荷的形式,在此条件下给出了弹塑性应力波参量的解析表达式.Milne 等[18,19]研究了恒定冲击载荷、连续卸载冲击载荷和渐增冲击载荷作用下弹塑性球面应力波的传播情况.Rapoport 等[20]对理想弹塑性介质中球面弹塑性波的传播进行了研究,针对球形空腔壁面施加恒定速度边界条件和恒定压力边界条件两种情况,给出了弹塑性波参量的解析表达式.Santos等[21]考虑了黏塑性Perzyna 型材料中空腔在恒定载荷作用下的动态响应,给出了介质中应力的计算公式并与有限元模拟结果进行了对比.肖建华等[22]对点爆炸震源产生的地震子波进行了求解.Chen等[23]研究了爆炸载荷加载下的弹塑性球面应力波问题,通过求解波动方程,给出了弹性区和塑性区应力波参量的解析解,得到了弹塑性边界上的应力时间历程.

目前针对弹塑性介质中球面应力波的理论研究,都没有考虑载荷卸载速率不同情况下应力波的传播情况.事实上,外加载荷卸载速率会影响弹塑性球面应力波传播情况,载荷卸载速率不同情况下应力波传播过程也存在区别.利用应力波理论对不同载荷卸载速率情况下的弹塑性球面应力波传播过程进行分析,可以加深对应力波传播衰减机理的认识,便于进一步分析弹塑性球面应力波的传播规律.此外,现有理论研究缺少对介质进入塑性后卸载区域的计算方法,也需要进一步补充研究.

本文基于线性硬化塑性本构模型,结合球面波波动方程,分析了不同卸载速率的冲击载荷加载下弹塑性球面应力波的传播过程,并通过解析求解波动方程,给出了应力波参量的精确解.通过与现有理论方法、数值模拟结果进行对比,验证了本文理论模型和方法的正确性,并对理论方法的适用性和计算过程中引入假设造成的误差进行了分析.

2 弹塑性球面应力波场的理论计算方法

2.1 物理模型

设无限大的固体介质中有一个半径为r0的球形空腔,冲击载荷均匀作用于空腔内壁,t时刻空腔内壁所受压力为p(t).其中p(t)在t=0 时刻从0 突跃至峰值P0,随后在t＞ 0 区间内逐渐减小.

在球对称应力应变状态下,固体材料容变律可以写为

其中σr和σθ分别为径向应力和环向应力,εr和εθ分别为径向应变和环向应变,λ和µ为拉梅常数.畸变律采用线性硬化塑性模型:

此处忽略了反向塑性加载.式中θ0为到达初始屈服强度时的偏应变;Y=2µθ0为屈服强度;θm为最大偏应变;Gp为塑性剪切模量.畸变律曲线如图1所示.

图1 畸变律曲线Fig.1.Curve of the distortion law.

在小变形条件下,球面波传播的波动方程为

其中,r和t分别为空间径向坐标和时间坐标,u为介质位移,ρ为介质密度.几何方程为

模型的初始条件和边界条件分别为

2.2 载荷卸载速率对应力波传播的影响

对上述模型进行解析求解的难点在于,介质本构方程中畸变律在不同阶段分别满足方程(2)—(4),只有确定了介质处于哪一阶段,才能代入对应的畸变律方程进行求解.

根据介质本构关系,弹性球面应力波和塑性球面应力波的传播速度均为常数.由于Gp＜µ,所以塑性波波速小于弹性波波速,在弹塑性球面波传播过程中,强间断弹性波波阵面和强间断塑性波波阵面在介质中依次向前传播.尽管载荷是连续卸载的形式,但外加载荷对应于边界处的径向应力,而介质的加卸载条件考察的是偏应力,因此介质中塑性波波阵面后仍可能有一段塑性加载过程,此加载过程的长短与外加载荷卸载速率相关.

假设介质中任意位置处完成塑性加载过程进入卸载后,不再进入塑性加载阶段,则介质中的弹塑性球面应力波传播过程可以分为3 种情况,在r-t平面上表示如图2 所示.

图2 应力波传播示意图 (a)载荷卸载较快;(b)载荷卸载较慢;(c)载荷卸载速率介于上述二者之间Fig.2.Schematic diagrams of the stress wave propagation (a)Rapid unloading;(b)slow unloading;(c)middle unloading.

当外加载荷卸载较快时,介质在塑性波波阵面后直接进入卸载过程,没有塑性加载阶段,在r-t平面上的应力波传播过程如图2(a)所示.图中的粗实线表示强间断,细实线表示弱间断.OE为强间断弹性波波阵面,其斜率为弹性波波速的倒数,在OE前方为未扰动区.OP为强间断塑性波波阵面,其斜率为塑性波波速的倒数,跨越OP后立即进入卸载阶段.当强间断塑性波波阵面传播至P点时,偏应力衰减到屈服强度,塑性波完全衰减,形成以弹性波波速向前后分别传播的两个弱间断PC和PD,PD在边界发生反射,成为弱间断DF并向前传播.根据上述分析,区域EOr为未扰动区,区域EOPP1为弹性区,在塑性波波阵面OP上发生塑性加载,区域P1POt为卸载区.

当外加载荷卸载较慢时,介质在塑性波波阵面后继续发生塑性加载行为,且塑性加载区较大,加卸载边界不会与强间断塑性波波阵面相交,在r-t平面上的应力波传播过程如图2(b)所示.初始的应力波传播情况与图2(a)相同,但跨越OP后介质继续发生塑性加载.当强间断塑性波波阵面传播至P点时,偏应力衰减到屈服强度,塑性波继续以弱间断的形式向前传播.因为几何扩散效应会使满足屈服条件的塑性应力波随传播而衰减,从而转化为弹性波,所以弱间断塑性波的传播速度小于强间断塑性波[24],且其传播速度不再是常数.弱间断塑性波波阵面沿PD传播至D点后完全衰减为弹性波.边界在C点处开始进入卸载,并向前传播弱间断CG,CD为加卸载边界.

当外加载荷卸载速率介于上述两种情况之间时,介质在塑性波波阵面后继续发生塑性加载行为,但塑性加载区较小,加卸载边界与强间断塑性波波阵面相交,在r-t平面上的应力波传播过程如图2(c)所示.此时加卸载边界CD与塑性波波阵面相交于点D.跨越OD段后,介质继续进行塑性加载,跨越DP段后介质直接进入卸载阶段.

由此确定了不同情况下r-t平面上弹性区、塑性加载区和卸载区的分布情况,就可以在各区域将对应的本构方程代入波动方程进行解析求解,得到应力波传播过程中参量的解.

2.3 应力波参量求解

在不同区域求解对应的应力波波动方程,可以得到应力波参量的解.弹性区波动方程为

其中c为弹性波波速:

引入位移势函数ϕ,满足:

则波动方程(10)可化为

由于弹性前驱波中只存在右行波,上式的通解形式为

结合边界条件可以得到f的表达式,从而计算得到介质位移及其他应力波参量.

塑性加载区波动方程为

其中cp为塑性波波速:

利用位移势函数可得:

上式通解形式为

结合边界条件可以得到f1和f2的解.

卸载区波动方程为

将位移u进行分解[16]:

其中U满足弹性区波动方程,可以利用求解弹性区方程的方法进行求解,但此时需同时考虑左行波和右行波.g满足方程:

结合边界条件易得g的表达式.

利用上述方法,在应力波传播的r-t图中不同区域分别求解对应的方程,区域交界处采用位移连续条件,就可以得到全场位移的解.将位移对时间求导即为质点速度.利用几何关系(6)和(7)以及本构关系(1)—(4),即可得到应变和应力的解.需要注意的是,为便于求解,本文引入了部分假设,主要包括将图2(b)中弱间断塑性波波阵面PD和加卸载边界CD近似为直线,忽略特征线DH和HK,以及将图2(c)中加卸载边界CD近似为直线,假设CD与左行特征线平行等.

3 结果与讨论

本节利用弹塑性球面应力波传播过程分析方法和应力波参量解析计算方法,计算弹塑性应力波参量在传播过程中的变化情况,并与现有理论方法以及动力学软件AUTODYN 计算结果进行对比,验证本文方法的合理性.最后对本文方法的适用性及计算过程中引入假设造成的误差进行了分析.

3.1 恒定载荷作用下的弹塑性应力波

Rapoport 等[20]给出了恒定冲击载荷作用下理想弹塑性材料中应力波参量的解析计算方法.在本文模型中,取塑性剪切模量Gp=0,外加载荷p(t)=P0为恒定值,也可以得到此情况下应力波参量的解.同时,利用动力学软件ANSYS AUTODYN ALE 求解器对相同的模型进行数值模拟.相关参数取值为:密度ρ=2630 kg·m—3,弹性模量E=50 GPa,泊松比ν=0.16,屈服强度Y=40 MPa,空腔半径r0=0.2 m,载荷峰值P0=0.2 GPa.

采用dr=5.6×10—4m 和dr=2.8×10—4m两种网格方案进行数值模拟,图3 给出了两种网格方案条件下计算得到空腔壁面处的位移时间历程和径向应力的空间分布结果,两种方案条件下的计算结果基本吻合,验证了网格收敛性,后续的模拟均采用网格dr=2.8×10—4m 进行计算.

图3 网格无关性 (a)空腔壁面位移时间历程;(b)径向应力空间分布Fig.3.Gird independence:(a)Time history of the cavity surface displacement;(b)radial stress distribution.

图4 给出了利用本文理论方法、文献[20]中理论方法和数值模拟方法得到的不同位置处径向应力时间历程和不同时刻径向应力的空间分布.由于载荷没有卸载,此情况下应力波传播过程为图2(b)所示的情况.Rapoport 等[20]没有考虑弱间断塑性应力波的传播和后续卸载过程,因此其方法在图4(b),4(c),4(e)和4(f)中不再适用.本文中的方法考虑了应力波传播的整个过程,但在求解波动方程过程中引入了部分假设,因此得到的结果与Rapoport 等[20]的理论方法以及数值模拟得到的结果有些偏差.但本文的方法适用范围更广,对于Rapoport 等[20]的理论方法不适用的区域,本文方法也可以求解.在r和t较小的区域,计算结果不受本文假设近似的影响,理论方法得到的结果与文献[20]中的理论方法结果完全相同.在其他区域,和数值模拟结果能够基本吻合,验证了理论方法的合理性.

图4 径向应力时间历程 (a)r=0.5 m;(b)r=0.8 m;(c)r=1.4 m 和径向应力空间分布 (d)t=0.1 ms;(e)t=0.2 ms;(f)t=0.3 ms (恒定冲击载荷)Fig.4.Time history of the radial stress:(a)r=0.5 m;(b)r=0.8 m;(c)r=1.4 m and the radial stress distribution:(d)t=0.1 ms;(e)t=0.2 ms;(f)t=0.3 ms (constant impact loading).

3.2 指数衰减载荷作用下的弹塑性应力波

取塑性剪切模量Gp=12.5 MPa,外加载荷p(t)=P0exp(—t/T0)为指数衰减形式,其中P0为载荷峰值,T0为载荷衰减指数.T0的大小会影响载荷卸载速率,进而影响应力波传播过程.

给定T0,首先需要判断应力波传播过程属于图2 中的哪种情况.首先按照图2(b)所示情况进行计算并得到C点的时间坐标tC,若tC=0,则应力波传播中不包含塑性加载区,为图2(a)所示的情况.若tC＞ 0,再判断加卸载边界CD是否与强间断塑性波波阵面OP相交,若不相交,则为图2(b)所示的情况,否则为图2(c)所示的情况.

按照上述方法计算,载荷衰减指数T0分别取0.01 ms,0.2 ms 和0.03 ms 时,对应的应力波传播过程分别为图2(a)—2(c)所示的情况.图5—图7分别给出了衰减指数T0的3 个取值情况下,本文理论方法和数值模拟方法得到的不同位置处径向应力时间历程和不同时刻径向应力的空间分布.

图5 径向应力时间历程 (a)r=0.3 m;(b)r=0.5 m;(c)r=0.6 m 和径向应力空间分布 (d)t=0.05 ms;(e)t=0.2 ms;(f)t=0.6 ms (T0=0.01 ms)Fig.5.Time history of the radial stress:(a)r=0.3 m;(b)r=0.5 m;(c)r=0.6 m and the radial stress distribution:(d)t=0.05 ms;(e)t=0.2 ms;(f)t=0.6 ms (T0=0.01 ms).

图7 径向应力时间历程 (a)r=0.4 m;(b)r=0.6 m;(c)r=0.7 m 和径向应力空间分布 (d)t=0.1 ms;(e)t=0.2 ms;(f)t=0.6 ms (T0=0.03 ms)Fig.7.Time history of the radial stress:(a)r=0.4 m;(b)r=0.6 m;(c)r=0.7 m;and the radial stress distribution:(d)t=0.1 ms;(e)t=0.2 ms;(f)t=0.6 ms (T0=0.03 ms).

从上述结果可以看出,3 种情况下本文理论方法计算得到的结果与数值模拟结果能够基本吻合,验证了应力波传播过程分类和参数理论计算方法的合理性.在冲击载荷作用下,介质中首先产生弹塑性双波结构,随着应力波的传播,塑性波逐渐衰减,最终强间断完全消失,成为弹性波向前传播.由于几何扩散效应,在传播过程中弹性波峰值也逐渐衰减.由于发生了塑性变形,在应力波传播过后,塑性区仍会存在残余应力,残余应力不会衰减,而是永久留存于介质中.当载荷衰减指数增大时,介质中传播的应力波脉宽以及残余应力均会增大.

3.3 理论方法适用性及误差分析

本文的理论计算方法是在小变形假设下得到的,因此载荷峰值P0不能过高,以保证介质满足小变形条件.当载荷峰值P0很高时,介质发生大变形,控制方程形式将发生改变,本文方法不再适用.

在对应力波传播过程分析中,本文假设介质中任意位置处完成塑性加载过程进入卸载后,不再进入塑性加载阶段,因此要求冲击载荷P(t)在峰值后为连续卸载形式,不能出现二次加载.此外,对于介质本构模型,本文忽略了反向塑性加载阶段,实际问题中介质可能出现反向塑性加载,但只存在于边界r=r0附近,因此当出现反向塑性加载现象时,本文理论方法的计算结果在边界r=r0附近不再适用,但对其他区域仍然适用.

在满足上述物理模型假设条件下,本文结果与数值模拟结果仍存在一定偏差,这是由于求解过程中引入的假设造成的.本文引入假设导致应力波峰值后的计算结果与实际有一定的偏差,因此定义径向应力误差:

其中σr为本文计算得到的径向应力,σr,N为数值模拟得到的径向应力,对于时间历程取x=t,空间分布则取x=r,则e可以反映本文计算结果与数值模拟结果在[x1,x2]区间内的平均相对偏差大小.对图4—图7 中的所有结果,按 (22)式计算,得到的误差如表1 所示.

从表1 可以看出,大多数计算结果的相对偏差均小于10%.图4 和图6 中存在部分结果相对偏差较大,这两个算例均为载荷卸载较慢的情况,对应的应力波传播情况为图2(b),由于直线化了弱间断塑性波波阵面并忽略了卸载区弱间断传播的特征线DH和HK,导致应力波峰值后的计算结果与实际有所偏差.图7(e)结果相对偏差较大,这是由于在其对应的应力波传播情况中,假设了加卸载边界CD与左行特征线平行,导致边界r=r0附近的计算结果与实际有所偏差.r=r0附近的计算结果与实际有所偏差.

图6 径向应力时间历程 (a)r=0.5 m;(b)r=0.75 m;(c)r=1.25 m 和径向应力空间分布 (d)t=0.1 ms;(e)t=0.2 ms;(f)t=0.6 ms (T0=0.2 ms)Fig.6.Time history of the radial stress:(a)r=0.5 m;(b)r=0.75 m;(c)r=1.25 m and the radial stress distribution:(d)t=0.1 ms;(e)t=0.2 ms;(f)t=0.6 ms (T0=0.2 ms).

表1 图4—图7 中各曲线相对偏差Table 1.Relative errors of the curves in Figs.4-7.

4 结论

在球形空腔内部冲击载荷作用下,介质中产生弹塑性球面应力波并向前传播.本文从理论分析的角度研究了弹塑性球面应力波的传播过程,并建立了应力波参量的解析求解方法.

外加冲击载荷的卸载速率对弹塑性球面应力波传播过程有很大影响.在r-t平面上,载荷卸载速率会影响塑性加载区的大小.当载荷卸载较快时,r-t平面上不存在塑性加载区.随着载荷卸载速率的降低,r-t平面上开始出现塑性加载区,但塑性加载区较小时,加卸载边界会与强间断塑性波波阵面相交.载荷卸载速率进一步降低时,r-t平面上塑性加载区面积扩大,加卸载边界不再与强间断塑性波波阵面相交,因此在强间断塑性波完全衰减后,还会以弱间断的形式继续传播.

在上述应力波传播分析结果基础上,利用解析方法对应力波参量进行求解.在不同载荷加载条件下,利用本文方法得到应力波参量的解并与Rapoport等[20]理论方法、数值模拟得到的结果进行对比.结果显示,本文方法与现有理论方法、数值模拟方法得到的结果基本吻合,能够反映应力波传播过程中的主要变化规律,且绝大多数结果相对数值模拟结果偏差均小于10%.现有理论方法没有考虑载荷卸载速率对应力波传播的影响,也没有给出卸载区域的求解方法,本文理论方法考虑了上述问题,相对现有方法具有更广泛的适用范围.