翟 越，李 楠，趙均海，王思維

（長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院，陜西 西安，710054）

巖石類材料在沖擊載荷作用下的強(qiáng)度、變形與破壞規(guī)律的研究對(duì)于國(guó)防工程、隧道工程與建筑抗震等方面具有重要的現(xiàn)實(shí)意義，但是由于沖擊載荷施加和測(cè)量的難度大，加之非均質(zhì)脆性材料的復(fù)雜性，使得這方面的研究進(jìn)展緩慢．各種沖擊荷載試驗(yàn)研究方法中，分離式霍普金森壓桿試驗(yàn)技術(shù)(Split Hopkinson Pressure Bar, 簡(jiǎn)稱SHPB)已成為研究中高應(yīng)變率（101-104s-1）下，材料動(dòng)態(tài)力學(xué)特性的最有效方法之一[1-4]．但是SHPB試驗(yàn)的測(cè)量精度與數(shù)據(jù)處理方面還存在準(zhǔn)確度較低和可靠性較差的問(wèn)題，因此隨著計(jì)算機(jī)技術(shù)的迅猛發(fā)展，對(duì)于SHPB試驗(yàn)的數(shù)據(jù)處理及數(shù)值模擬的研究成為熱點(diǎn)，并取得了一些成果[5-9]．這類問(wèn)題的研究中，應(yīng)力波在介質(zhì)中的傳播是一個(gè)核心問(wèn)題[10-11]．

本文構(gòu)建了應(yīng)力波在考慮損傷弱化及應(yīng)變率強(qiáng)化效應(yīng)的非線彈性巖石試件和高強(qiáng)度鋼制壓桿組成的桿系中傳播的波動(dòng)方程，并根據(jù)試件兩端與入射桿及透射桿接觸面的應(yīng)力波反射與透射關(guān)系，以及壓桿端部的自由邊界，給出初始條件和邊界條件，利用混合遺傳算法和有限差分法相結(jié)合，編制了優(yōu)化數(shù)值計(jì)算程序，反演分析出試件的彈性模量、分布參數(shù)等，進(jìn)而得到試件的應(yīng)力、應(yīng)變數(shù)值．通過(guò)計(jì)算應(yīng)力-應(yīng)變曲線和試驗(yàn)曲線對(duì)比，驗(yàn)證了非線彈性應(yīng)力波動(dòng)方程的可靠性，以及數(shù)值分析方法的適用性．

1 應(yīng)力波在細(xì)長(zhǎng)桿件中的波動(dòng)方程

1.1 應(yīng)力波在長(zhǎng)細(xì)桿件中的傳播理論

將SHPB試驗(yàn)的數(shù)值分析問(wèn)題，轉(zhuǎn)化為在物質(zhì)坐標(biāo)中研究應(yīng)力波在一組不同材料等直截面的細(xì)長(zhǎng)桿系中的縱向傳播問(wèn)題．取變形前(時(shí)間 0t= )質(zhì)點(diǎn)的空間位置作為物質(zhì)坐標(biāo)，桿系中心軸線為X軸，如圖1所示．設(shè)定桿在變形前的初始截面積0A(m2)、初始密度0ρ(kg/m3)和其它材料性能參數(shù)都與坐標(biāo)無(wú)關(guān)，截面形狀一般也無(wú)限制．

應(yīng)力波在各種介質(zhì)內(nèi)傳播的問(wèn)題中，波動(dòng)方程是最基本、最重要的理論基礎(chǔ)．波動(dòng)方程的建立可以通過(guò)聯(lián)立運(yùn)動(dòng)方程、動(dòng)力方程和介質(zhì)材料的物理方程，即本構(gòu)方程[12]，組成的方程組求解三個(gè)未知函數(shù)應(yīng)力),(tXσ，應(yīng)變),(tXε和質(zhì)子速度),(tXv，從而得到應(yīng)力波在該材料中傳播的波動(dòng)方程．下面分別探討應(yīng)力波在線彈性鋼制壓桿以及考慮損傷的非線彈性巖石試件中傳播的波動(dòng)方程．規(guī)定應(yīng)力和應(yīng)變均以拉為正，而質(zhì)點(diǎn)速度以沿X軸正方向?yàn)檎?，反之為?fù)．

圖1 質(zhì)點(diǎn)空間位置的物質(zhì)坐標(biāo)Fig.1 Space coordinates of material point

1.2 應(yīng)力波在線彈性細(xì)長(zhǎng)桿件中的波動(dòng)方程

運(yùn)動(dòng)學(xué)方程，即連續(xù)方程或質(zhì)量守恒方程，由位移u的單值連續(xù)條件就可得到聯(lián)系應(yīng)變?chǔ)藕蛈的相容性方程，如式(1)所示．

動(dòng)力學(xué)方程，即動(dòng)量守恒方程，利用連續(xù)性條件可以得到：

材料本構(gòu)方程，即物理方程，可以表示為應(yīng)力與應(yīng)變的連續(xù)函數(shù)．SHPB裝置的撞擊桿、入射桿和透射桿都是由高強(qiáng)度不銹鋼制成，其材料在一般加載條件下可以看成理想線彈性體，因此其符合胡克定律，即：

式中，0E為材料的初始彈性模量，單位MPa．聯(lián)立求解式(1)、式(2)和本構(gòu)方程(3)，則應(yīng)力波在線彈性材料中的波動(dòng)方程可表示為：

1.3 應(yīng)力波在考慮損傷的動(dòng)態(tài)非線彈性桿件中傳播的波動(dòng)方程

1.3.1 考慮損傷的動(dòng)態(tài)非線彈性本構(gòu)模型

當(dāng)加載應(yīng)變率較大時(shí)（100s-1以上），巖石類材料的應(yīng)力-應(yīng)變曲線呈現(xiàn)出明顯的非線性，此時(shí)用線彈性本構(gòu)模型來(lái)描述顯然不太適合．國(guó)內(nèi)外學(xué)者提出多種的非線性本構(gòu)方程，其中，Saenz[13]建立的非線性本構(gòu)模型，在鋼筋混凝土的有限元分析中有廣泛的應(yīng)用[14]，其公式如下所示：

式中：B、C分別為材料的特征參數(shù)；fε為峰值應(yīng)力對(duì)應(yīng)的應(yīng)變．

為了考慮損傷對(duì)巖石類材料力學(xué)性能的影響，根據(jù)應(yīng)變等效原理，受損材料的本構(gòu)方程可由損傷后的有效應(yīng)力來(lái)取代無(wú)損材料本構(gòu)方程中的名義應(yīng)力，從而構(gòu)建出[15]，因此本文引入考慮閾值的統(tǒng)計(jì)損傷演化因子w，如下式．

式中，m和a為損傷因子的分布參數(shù)；sσ為損傷應(yīng)力閾值，單位Pa．當(dāng)應(yīng)力小于它時(shí)，不考慮損傷．

同時(shí)，為了考慮沖擊荷載作用下的應(yīng)變率相關(guān)性，引入應(yīng)變率強(qiáng)化因子dR，如式(7)所示：

式中，λ是材料參數(shù)，反映巖石類材料受壓時(shí)的率效應(yīng)；ε˙表示在動(dòng)荷載作用下反映材料動(dòng)態(tài)響應(yīng)的應(yīng)變率，單位為s-1，在SHPB試驗(yàn)中一般取平均應(yīng)變率值．

將式(6)和式(7)分別代入式(5)，則構(gòu)建出考慮損傷的動(dòng)態(tài)本構(gòu)方程，如式(8)所示：

式中，參數(shù)意義同前所述，其中sσ、fε和˙ε分別可由試驗(yàn)得到，其他特征參數(shù)可以在試驗(yàn)研究基礎(chǔ)上，通過(guò)優(yōu)化反演分析確定出來(lái)．

1.3.2 應(yīng)力波在非線性粘彈性材料中的波動(dòng)方程

建立應(yīng)力波在非線彈性材料中的波動(dòng)方程時(shí)，采用上述式(8)作為巖石類材料的本構(gòu)方程，將其與運(yùn)動(dòng)學(xué)方程(1)和動(dòng)力學(xué)方程(2)聯(lián)立求解，可得考慮損傷的非線彈性應(yīng)力波波動(dòng)方程微分形式，如式(9)所示：

2 巖石類材料動(dòng)態(tài)特性的數(shù)值分析

2.1 數(shù)值計(jì)算基本假定

為了便于數(shù)值計(jì)算，在可靠性可以接受的情況下，作如下假定：首先，桿件在變形時(shí)橫截面始終保持為平面，沿橫截面只有均布的軸向應(yīng)力，因此各材料參數(shù)都只是位置X和時(shí)間t的函數(shù)，整個(gè)問(wèn)題簡(jiǎn)化為一維問(wèn)題，而忽略桿中質(zhì)點(diǎn)橫向運(yùn)動(dòng)的慣性作用，即不計(jì)橫向收縮或膨脹對(duì)動(dòng)能的影響．根據(jù)桿中彈性波的精確解可知，只要波長(zhǎng)比桿的橫向尺寸大很多時(shí)，這一近似假定所引起的誤差是允許忽略的[2]；其次，由于應(yīng)力波波速很高，在通過(guò)微元體的瞬間，微元體還來(lái)不及和鄰近的其它微元體進(jìn)行熱量交換，因此可近似地認(rèn)為沖擊加載過(guò)程是絕熱的，無(wú)需列出能量守恒方程，就可得到關(guān)于變量σ、ε和v的封閉的控制方程組．?dāng)?shù)值分析及試驗(yàn)所用的SHPB裝置示意圖如圖2所示．

圖2 SHPB結(jié)構(gòu)圖Fig.2 Sketch map of SHPB apparatus

2.2 壓桿和試件中應(yīng)力波波動(dòng)方程的有限差分解

有限差分是一種常用的數(shù)值解法，它是用差商代替偏導(dǎo)數(shù)，得到在微分方程的差分形式，通過(guò)解差分方程得到微分方程解的近似值．應(yīng)力波在鋼制壓桿中的波動(dòng)方程式(4)，可用有限差分法的中心差分格式表述如下：

試件中的應(yīng)力波波動(dòng)方程式（9）的有限差分解可按上述中心差分方法得到．

2.3 初始條件和邊界條件

2.3.1 初始條件

假設(shè)在加載過(guò)程中，子彈和入射桿始終保持接觸狀態(tài)，此時(shí)在接觸面上的軸向方向，兩桿的質(zhì)點(diǎn)位移相同，即初始條件為：時(shí)間點(diǎn)為1時(shí)（j=1），入射桿初始端與子彈接觸面上質(zhì)點(diǎn)的位移 Wz(i,1)等于接觸面上子彈質(zhì)點(diǎn)速度 Vi( 1)與時(shí)間的乘積：

由此，將試驗(yàn)測(cè)得的子彈撞擊波形整形器改進(jìn)后的入射桿時(shí)產(chǎn)生的半正弦加載應(yīng)力波（如圖3）通過(guò)初始條件直接加載在入射桿的初始端．

2.3.2 邊界條件

要數(shù)值模擬應(yīng)力波在整個(gè)SHPB壓桿和試件中的傳播，必然要確定壓桿兩端的邊界條件，以及試件和壓桿的接觸面上的邊界條件．這里將試件在動(dòng)態(tài)沖擊荷載作用下的慣性和端頭摩擦效應(yīng)忽略不計(jì)．

(1) 桿件端部邊界條件

當(dāng)端部固定時(shí)，由于入射波 R (C0,t)與反射波F(C0,t)的疊加，其邊界條件為端頭處的位移 u0始終為零，而應(yīng)力σ0為入射時(shí)的2倍．其中C0為應(yīng)力波波速．

當(dāng)端部自由時(shí)，同樣由于入射波 R (C0,t)與反射波F(C0,t)的疊加，邊界條件為自由端的應(yīng)力為零，而位移為入射時(shí)為2倍．

(2) 兩種介質(zhì)接觸面上的邊界條件

由于壓桿和巖石試件的阻抗(00Cρ)存在較大差異，根據(jù)一維應(yīng)力波在不同介質(zhì)的傳播理論[11]，SHPB試驗(yàn)中的應(yīng)力波會(huì)在兩種介質(zhì)的接觸面上發(fā)生反射與透射．試驗(yàn)中測(cè)得的入射波和反射波是在入射桿上測(cè)得的，而透射波是在透射桿上測(cè)得的，也就是第一個(gè)界面的透射波通過(guò)試件傳播到第二個(gè)界面產(chǎn)生的透射波，因此試驗(yàn)測(cè)量的三個(gè)波并不是在同一界面上發(fā)生的．在相同壓桿和入射波R(C0,t)的條件下，反射波 F (C0,t)和透射波 T (C0,t)的大小是由試件的阻抗所決定的，計(jì)算公式如下所示：

式中，F(xiàn)和T分別為應(yīng)力波由壓桿向巖石類材料傳播時(shí)的反射系數(shù)和透射系數(shù)，可由下式計(jì)算：

式中，T′為應(yīng)力波由巖石試件向透射桿透射的系數(shù)；L為應(yīng)力波由入射壓桿向試件透射后，再由試件向透射桿透射的兩次透射系數(shù)，可表示為：

2.4 數(shù)值分析程序及計(jì)算結(jié)果分析

2.4.1 花崗巖SHPB試驗(yàn)及數(shù)值分析程序

為了給數(shù)值分析提供試驗(yàn)基礎(chǔ)，對(duì)花崗巖試件進(jìn)行了SHPB試驗(yàn)，所用裝置如圖2所示．SHPB系統(tǒng)的壓桿和子彈的材料均為高強(qiáng)度不銹鋼，直徑皆為50 mm，入射桿、透射桿和吸收桿的長(zhǎng)度皆為2 000 mm，子彈長(zhǎng)度為300 mm，并采用了直徑為18 mm、厚度為1 mm的圓形銅片作為波形整形器，將應(yīng)力波改進(jìn)為半正弦波[16]，如圖3所示．試驗(yàn)中控制氣壓為0.8 MPa，子彈沖擊速度為12.84 m/s，測(cè)得平均應(yīng)變率為50 s-1．

圖3 加脈沖整形器后的波形圖Fig.3 Pulse curve with pulse shaper

基于上述有限差分法求解波動(dòng)方程、初始條件和邊界條件所組成的方程組，應(yīng)用VC++編寫(xiě)了巖石類材料的動(dòng)態(tài)分析程序．在數(shù)值計(jì)算時(shí)，應(yīng)得到非線性本構(gòu)方程中的特征參數(shù)．為此編制了自適應(yīng)混合遺傳算法[17]，并嵌入到有限差分程序中，實(shí)現(xiàn)每一次數(shù)值分析應(yīng)力波傳播的同時(shí)進(jìn)行優(yōu)化反演分析，從而確定波動(dòng)方程中的待定參數(shù)．?dāng)?shù)值分析中的物理參數(shù)選取如表1所示．

表1 數(shù)值模擬程序的計(jì)算參數(shù)選取Tab.1 Count parameter of numerical simulation program

2.4.2 數(shù)值分析結(jié)果

基于試驗(yàn)數(shù)據(jù)，通過(guò)有限差分?jǐn)?shù)值計(jì)算和優(yōu)化反演的混合程序，確定出花崗巖的動(dòng)態(tài)本構(gòu)模型的特征參數(shù)如表2所示．

表2 花崗巖非線彈性動(dòng)態(tài)本構(gòu)模型反演分析結(jié)果Tab.2 Inversion Analysis results of nonlinear elastic constitutive model

為了和試驗(yàn)結(jié)果比較，取數(shù)值分析中壓桿中心處的應(yīng)力應(yīng)變值，作為花崗巖試件的平均應(yīng)力和應(yīng)變．試驗(yàn)應(yīng)力-應(yīng)變曲線和計(jì)算應(yīng)力應(yīng)曲線的對(duì)比如圖4所示．圖中，兩曲線有較好的一致性，表明所采用考慮損傷弱化和應(yīng)變率強(qiáng)化的非線彈性本構(gòu)方程和在此基礎(chǔ)上推導(dǎo)出來(lái)的非線彈性應(yīng)力波波動(dòng)方程都是適用和可靠的，同時(shí)也驗(yàn)證了這種數(shù)值分析方法是合理的．

圖4 花崗巖的試驗(yàn)曲線和再生曲線（應(yīng)變率50s-1）Fig.4 Test curve and reproduction curve of granite

這種方法相對(duì)于SHPB傳統(tǒng)數(shù)據(jù)處理方法的兩個(gè)最突出的優(yōu)點(diǎn)，一是采用的數(shù)據(jù)可以采自壓桿的任意位置；二是不再要求入射桿和透射桿要滿足足夠的長(zhǎng)度．該方法適用于各種巖石類材料和低阻抗材料的動(dòng)態(tài)特性的研究，為SHPB裝置的試驗(yàn)結(jié)果分析與數(shù)據(jù)處理提出一種有效的方法．

3 結(jié) 論

(1) 基于考慮損傷弱化和應(yīng)變率強(qiáng)化的非線彈性動(dòng)態(tài)本構(gòu)方程，及相應(yīng)的運(yùn)動(dòng)方程和動(dòng)力方程，建立了應(yīng)力波在巖石類材料桿件中傳播的波動(dòng)方程；

(2) 在上述波動(dòng)方程的基礎(chǔ)上，利用有限差分法編制了可以分析了沖擊加載波由入射桿傳播到試件再到透射桿的整個(gè)加載過(guò)程的數(shù)值計(jì)算程序；

(3) 將數(shù)值程序與混合遺傳算法相結(jié)合，辨識(shí)出方程的特征參數(shù)，將理論計(jì)算結(jié)果與試驗(yàn)結(jié)果比較得到了可靠性驗(yàn)證．

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] 翟毅, 許金余, 王鵬輝. 纖維混凝土動(dòng)態(tài)壓縮力學(xué)性能的SHPB試驗(yàn)研究[J]. 西安建筑科技大學(xué)學(xué)報(bào): 自然科學(xué)版, 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] 江見(jiàn)鯨, 賀小崗. 工程結(jié)構(gòu)計(jì)算機(jī)仿真分析[M]. 北京:清華大學(xué)出版社, 1996.JANG Jianjing, HE Xiaogang. Computer simulation analysis of engineering structure[M]. Bejing: Tsinghua University press, 1996.

[6] 唐志平, 王禮立. SHPB試驗(yàn)的電腦化數(shù)據(jù)處理系統(tǒng)[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] 王禮立. 應(yīng)力波基礎(chǔ)[M]. 北京: 國(guó)防工業(yè)出版社, 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] 呂西林, 金國(guó)芳, 吳曉涵. 鋼筋混凝土結(jié)構(gòu)非線性有限元理論與應(yīng)用[M]. 上海: 同濟(jì)大學(xué)出版社, 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] 商懷帥, 楊魯生. 基于損傷理論的混凝土雙軸壓本構(gòu)模型[J]. 中南大學(xué)學(xué)報(bào): 自然科學(xué)版, 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] 翟越, 馬國(guó)偉, 趙均海, 等. 花崗巖和混凝土在沖擊荷載下的動(dòng)態(tài)性能比較研究[J]. 巖石力學(xué)與工程學(xué)報(bào),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] 翟越, 趙均海. 基于自適應(yīng)混合遺傳算法的巖石類材料動(dòng)態(tài)參數(shù)反演分析[J]. 地球科學(xué)與環(huán)境學(xué)報(bào), 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.