SHPB 试验中应力波在损伤非线弹性材料中的传播

翟 越，李 楠，赵均海，王思维

（长安大学地质工程与测绘学院，陕西 西安，710054）

岩石类材料在冲击载荷作用下的强度、变形与破坏规律的研究对于国防工程、隧道工程与建筑抗震等方面具有重要的现实意义，但是由于冲击载荷施加和测量的难度大，加之非均质脆性材料的复杂性，使得这方面的研究进展缓慢．各种冲击荷载试验研究方法中，分离式霍普金森压杆试验技术(Split Hopkinson Pressure Bar, 简称SHPB)已成为研究中高应变率（101-104s-1）下，材料动态力学特性的最有效方法之一[1-4]．但是SHPB试验的测量精度与数据处理方面还存在准确度较低和可靠性较差的问题，因此随着计算机技术的迅猛发展，对于SHPB试验的数据处理及数值模拟的研究成为热点，并取得了一些成果[5-9]．这类问题的研究中，应力波在介质中的传播是一个核心问题[10-11]．

本文构建了应力波在考虑损伤弱化及应变率强化效应的非线弹性岩石试件和高强度钢制压杆组成的杆系中传播的波动方程，并根据试件两端与入射杆及透射杆接触面的应力波反射与透射关系，以及压杆端部的自由边界，给出初始条件和边界条件，利用混合遗传算法和有限差分法相结合，编制了优化数值计算程序，反演分析出试件的弹性模量、分布参数等，进而得到试件的应力、应变数值．通过计算应力-应变曲线和试验曲线对比，验证了非线弹性应力波动方程的可靠性，以及数值分析方法的适用性．

1 应力波在细长杆件中的波动方程

1.1 应力波在长细杆件中的传播理论

将SHPB试验的数值分析问题，转化为在物质坐标中研究应力波在一组不同材料等直截面的细长杆系中的纵向传播问题．取变形前(时间 0t= )质点的空间位置作为物质坐标，杆系中心轴线为X轴，如图1所示．设定杆在变形前的初始截面积0A(m2)、初始密度0ρ(kg/m3)和其它材料性能参数都与坐标无关，截面形状一般也无限制．

应力波在各种介质内传播的问题中，波动方程是最基本、最重要的理论基础．波动方程的建立可以通过联立运动方程、动力方程和介质材料的物理方程，即本构方程[12]，组成的方程组求解三个未知函数应力),(tXσ，应变),(tXε和质子速度),(tXv，从而得到应力波在该材料中传播的波动方程．下面分别探讨应力波在线弹性钢制压杆以及考虑损伤的非线弹性岩石试件中传播的波动方程．规定应力和应变均以拉为正，而质点速度以沿X轴正方向为正，反之为负．

图1 质点空间位置的物质坐标Fig.1 Space coordinates of material point

1.2 应力波在线弹性细长杆件中的波动方程

运动学方程，即连续方程或质量守恒方程，由位移u的单值连续条件就可得到联系应变ε和u的相容性方程，如式(1)所示．

动力学方程，即动量守恒方程，利用连续性条件可以得到：

材料本构方程，即物理方程，可以表示为应力与应变的连续函数．SHPB装置的撞击杆、入射杆和透射杆都是由高强度不锈钢制成，其材料在一般加载条件下可以看成理想线弹性体，因此其符合胡克定律，即：

式中，0E为材料的初始弹性模量，单位MPa．联立求解式(1)、式(2)和本构方程(3)，则应力波在线弹性材料中的波动方程可表示为：

1.3 应力波在考虑损伤的动态非线弹性杆件中传播的波动方程

1.3.1 考虑损伤的动态非线弹性本构模型

当加载应变率较大时（100s-1以上），岩石类材料的应力-应变曲线呈现出明显的非线性，此时用线弹性本构模型来描述显然不太适合．国内外学者提出多种的非线性本构方程，其中，Saenz[13]建立的非线性本构模型，在钢筋混凝土的有限元分析中有广泛的应用[14]，其公式如下所示：

式中：B、C分别为材料的特征参数；fε为峰值应力对应的应变．

为了考虑损伤对岩石类材料力学性能的影响，根据应变等效原理，受损材料的本构方程可由损伤后的有效应力来取代无损材料本构方程中的名义应力，从而构建出[15]，因此本文引入考虑阈值的统计损伤演化因子w，如下式．

式中，m和a为损伤因子的分布参数；sσ为损伤应力阈值，单位Pa．当应力小于它时，不考虑损伤．

同时，为了考虑冲击荷载作用下的应变率相关性，引入应变率强化因子dR，如式(7)所示：

式中，λ是材料参数，反映岩石类材料受压时的率效应；ε˙表示在动荷载作用下反映材料动态响应的应变率，单位为s-1，在SHPB试验中一般取平均应变率值．

将式(6)和式(7)分别代入式(5)，则构建出考虑损伤的动态本构方程，如式(8)所示：

式中，参数意义同前所述，其中sσ、fε和˙ε分别可由试验得到，其他特征参数可以在试验研究基础上，通过优化反演分析确定出来．

1.3.2 应力波在非线性粘弹性材料中的波动方程

建立应力波在非线弹性材料中的波动方程时，采用上述式(8)作为岩石类材料的本构方程，将其与运动学方程(1)和动力学方程(2)联立求解，可得考虑损伤的非线弹性应力波波动方程微分形式，如式(9)所示：

2 岩石类材料动态特性的数值分析

2.1 数值计算基本假定

为了便于数值计算，在可靠性可以接受的情况下，作如下假定：首先，杆件在变形时横截面始终保持为平面，沿横截面只有均布的轴向应力，因此各材料参数都只是位置X和时间t的函数，整个问题简化为一维问题，而忽略杆中质点横向运动的惯性作用，即不计横向收缩或膨胀对动能的影响．根据杆中弹性波的精确解可知，只要波长比杆的横向尺寸大很多时，这一近似假定所引起的误差是允许忽略的[2]；其次，由于应力波波速很高，在通过微元体的瞬间，微元体还来不及和邻近的其它微元体进行热量交换，因此可近似地认为冲击加载过程是绝热的，无需列出能量守恒方程，就可得到关于变量σ、ε和v的封闭的控制方程组．数值分析及试验所用的SHPB装置示意图如图2所示．

图2 SHPB结构图Fig.2 Sketch map of SHPB apparatus

2.2 压杆和试件中应力波波动方程的有限差分解

有限差分是一种常用的数值解法，它是用差商代替偏导数，得到在微分方程的差分形式，通过解差分方程得到微分方程解的近似值．应力波在钢制压杆中的波动方程式(4)，可用有限差分法的中心差分格式表述如下：

试件中的应力波波动方程式（9）的有限差分解可按上述中心差分方法得到．

2.3 初始条件和边界条件

2.3.1 初始条件

假设在加载过程中，子弹和入射杆始终保持接触状态，此时在接触面上的轴向方向，两杆的质点位移相同，即初始条件为：时间点为1时（j=1），入射杆初始端与子弹接触面上质点的位移 Wz(i,1)等于接触面上子弹质点速度 Vi( 1)与时间的乘积：

由此，将试验测得的子弹撞击波形整形器改进后的入射杆时产生的半正弦加载应力波（如图3）通过初始条件直接加载在入射杆的初始端．

2.3.2 边界条件

要数值模拟应力波在整个SHPB压杆和试件中的传播，必然要确定压杆两端的边界条件，以及试件和压杆的接触面上的边界条件．这里将试件在动态冲击荷载作用下的惯性和端头摩擦效应忽略不计．

(1) 杆件端部边界条件

当端部固定时，由于入射波 R (C0,t)与反射波F(C0,t)的叠加，其边界条件为端头处的位移 u0始终为零，而应力σ0为入射时的2倍．其中C0为应力波波速．

当端部自由时，同样由于入射波 R (C0,t)与反射波F(C0,t)的叠加，边界条件为自由端的应力为零，而位移为入射时为2倍．

(2) 两种介质接触面上的边界条件

由于压杆和岩石试件的阻抗(00Cρ)存在较大差异，根据一维应力波在不同介质的传播理论[11]，SHPB试验中的应力波会在两种介质的接触面上发生反射与透射．试验中测得的入射波和反射波是在入射杆上测得的，而透射波是在透射杆上测得的，也就是第一个界面的透射波通过试件传播到第二个界面产生的透射波，因此试验测量的三个波并不是在同一界面上发生的．在相同压杆和入射波R(C0,t)的条件下，反射波 F (C0,t)和透射波 T (C0,t)的大小是由试件的阻抗所决定的，计算公式如下所示：

式中，F和T分别为应力波由压杆向岩石类材料传播时的反射系数和透射系数，可由下式计算：

式中，T′为应力波由岩石试件向透射杆透射的系数；L为应力波由入射压杆向试件透射后，再由试件向透射杆透射的两次透射系数，可表示为：

2.4 数值分析程序及计算结果分析

2.4.1 花岗岩SHPB试验及数值分析程序

为了给数值分析提供试验基础，对花岗岩试件进行了SHPB试验，所用装置如图2所示．SHPB系统的压杆和子弹的材料均为高强度不锈钢，直径皆为50 mm，入射杆、透射杆和吸收杆的长度皆为2 000 mm，子弹长度为300 mm，并采用了直径为18 mm、厚度为1 mm的圆形铜片作为波形整形器，将应力波改进为半正弦波[16]，如图3所示．试验中控制气压为0.8 MPa，子弹冲击速度为12.84 m/s，测得平均应变率为50 s-1．

图3 加脉冲整形器后的波形图Fig.3 Pulse curve with pulse shaper

基于上述有限差分法求解波动方程、初始条件和边界条件所组成的方程组，应用VC++编写了岩石类材料的动态分析程序．在数值计算时，应得到非线性本构方程中的特征参数．为此编制了自适应混合遗传算法[17]，并嵌入到有限差分程序中，实现每一次数值分析应力波传播的同时进行优化反演分析，从而确定波动方程中的待定参数．数值分析中的物理参数选取如表1所示．

表1 数值模拟程序的计算参数选取Tab.1 Count parameter of numerical simulation program

2.4.2 数值分析结果

基于试验数据，通过有限差分数值计算和优化反演的混合程序，确定出花岗岩的动态本构模型的特征参数如表2所示．

表2 花岗岩非线弹性动态本构模型反演分析结果Tab.2 Inversion Analysis results of nonlinear elastic constitutive model

为了和试验结果比较，取数值分析中压杆中心处的应力应变值，作为花岗岩试件的平均应力和应变．试验应力-应变曲线和计算应力应曲线的对比如图4所示．图中，两曲线有较好的一致性，表明所采用考虑损伤弱化和应变率强化的非线弹性本构方程和在此基础上推导出来的非线弹性应力波波动方程都是适用和可靠的，同时也验证了这种数值分析方法是合理的．

图4 花岗岩的试验曲线和再生曲线（应变率50s-1）Fig.4 Test curve and reproduction curve of granite

这种方法相对于SHPB传统数据处理方法的两个最突出的优点，一是采用的数据可以采自压杆的任意位置；二是不再要求入射杆和透射杆要满足足够的长度．该方法适用于各种岩石类材料和低阻抗材料的动态特性的研究，为SHPB装置的试验结果分析与数据处理提出一种有效的方法．

3 结 论

(1) 基于考虑损伤弱化和应变率强化的非线弹性动态本构方程，及相应的运动方程和动力方程，建立了应力波在岩石类材料杆件中传播的波动方程；

(2) 在上述波动方程的基础上，利用有限差分法编制了可以分析了冲击加载波由入射杆传播到试件再到透射杆的整个加载过程的数值计算程序；

(3) 将数值程序与混合遗传算法相结合，辨识出方程的特征参数，将理论计算结果与试验结果比较得到了可靠性验证．

References

[1] DAI feng, HUANG Sheng, XIA Kaiwen, et al Some fundamental issues in dynamic compression and tension tests of rock using split Hopkinson pressure bar[J]. Rock Mech Rock Eng, 2010, 43: 657-666.

[2] EZIO Cadoni. Dynamic characterization of orthogenesis rock subjected to intermediate and high strain rates in tension[J]. Rock Mech Eng, 2010, 43: 667-676.

[3] YAN fei, FENG Xiating, CHEN Rong, et al. Dynamic tensile failure of the rock interface between tuff and basalt[J]. Rock Mech Rock Eng, 2012, 45: 341-348.

[4] 翟毅, 许金余, 王鹏辉. 纤维混凝土动态压缩力学性能的SHPB试验研究[J]. 西安建筑科技大学学报: 自然科学版, 2009, 41(1): 141-148.ZHAI Yi, XU Jinyu, WANG Penghui. Dynamic compressive testing and mechanical behavior of fiber reinforced concrete using a split Hopkinson Pressure Bar[J]. J. Xi’an Univ. of Arch. &Tech: Natural Science Edition, 2009, 41(1): 141-148.

[5] 江见鲸, 贺小岗. 工程结构计算机仿真分析[M]. 北京:清华大学出版社, 1996.JANG Jianjing, HE Xiaogang. Computer simulation analysis of engineering structure[M]. Bejing: Tsinghua University press, 1996.

[6] 唐志平, 王礼立. SHPB试验的电脑化数据处理系统[J].爆炸与冲击, 1986, 6(4): 320-327.TANG Zhiping, WANG Lili. Computerized data processing system of SHPB experiment[J]. Explosion and Shock Waves, 1986, 6(4): 320-327.

[7] LIFSHITZ J M, LEBER H. Data processing in the split Hopkinson pressure bar tests[J]. International Journal of Impact Engineering, 1994, 15(6): 723-733.

[8] LI X B, LO T S, Zhao J, et al. Oscillation Elimination in the Hopkinson bar apparatus and resultant complete dynamic stress-strain curve for rocks[J]. International Journal of Rock Mechanics & Mining Sciences, 2000, 37:1055-1060.

[9] HAO Y, HAO H. Numerical investigation of the dynamic compressive behaviour of rock materials at high strain rate[J]. Rock Mech Rock Eng, 2013, 46: 373 -388.

[10] FOLLANSBEE P S, FRANTZ C. Wave propagation in the split Hopkinson pressure bar[J]. Journal of Engineering Materials Technology, 1983, 105: 61-66.

[11] FAN L F, REN F, MA G W. Experimental study on viscoelastic behavior of sedimentary[J]. Rock Mech Rock Eng, 2012, 45: 433-438.

[12] 王礼立. 应力波基础[M]. 北京: 国防工业出版社, 2005.WANG Lili. Foundations of Stress Waves[M]. Beijing:National Defense Press, 2005.

[13] SAENZ L P. Discussion of Equation for the Stress-strain Curve of Concrete by Desay and Krishnan[J]. Journal of ACI, 1964, 61(9): 1229-1235.

[14] 吕西林, 金国芳, 吴晓涵. 钢筋混凝土结构非线性有限元理论与应用[M]. 上海: 同济大学出版社, 1996.LÜ Xilin, JIN Guofang, WU Xiaohan. The structure of the nonlinear finite element theory and application of reinforced concrete[M]. Shanghai: Tongji University press, 1996.

[15] 商怀帅, 杨鲁生. 基于损伤理论的混凝土双轴压本构模型[J]. 中南大学学报: 自然科学版, 2013, 44(1): 340-344.SHANG Huaishuai, YANG Lusheng. Constitutive model of damage of concrete under biaxial compression[J].Journal of Central South University: Science and Technology, 2013, 44(1): 340-344.

[16] 翟越, 马国伟, 赵均海, 等. 花岗岩和混凝土在冲击荷载下的动态性能比较研究[J]. 岩石力学与工程学报,2007, 26(4): 762-768.ZHAI Yue, MA Guowei, ZHAO Junhai, et al. Dynamic capability of granite and concrete under impact compressive loading[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(4): 762-768.

[17] 翟越, 赵均海. 基于自适应混合遗传算法的岩石类材料动态参数反演分析[J]. 地球科学与环境学报, 2008,30(3): 286-291.ZHAI Yue, ZHAO Junhai. Inverse analysis based on adaptive hybrid genetic algorithms for dynamic characteristic parameters of rock materials[J]. Journal of Earth Sciences and Environment, 2008, 30(3): 286-291.