王智環(huán) 賈雷明 何增 田宙

1)(清華大學(xué)工程物理系,北京 100084)

2)(西北核技術(shù)研究所,西安 710024)

基于線性硬化塑性本構(gòu)模型,建立了沖擊載荷作用下彈塑性球面應(yīng)力波場的理論求解方法.首先,分析了沖擊載荷卸載速率對球面應(yīng)力波傳播的影響,得到了3 種不同類型的應(yīng)力波傳播圖像.在此基礎(chǔ)上,建立了彈性階段、塑性加載階段以及卸載階段球面波動方程的理論求解方法,給出了質(zhì)點位移、質(zhì)點速度、應(yīng)力和應(yīng)變等物理量的計算方案.與已有理論方法相比,該方法考慮了不同載荷卸載速率條件下應(yīng)力波的不同傳播情況,并且給出了卸載階段應(yīng)力波參量計算方法,具有更廣的適用范圍.利用上述方法計算了恒定沖擊載荷和不同指數(shù)衰減沖擊載荷作用下彈塑性球面應(yīng)力波場參量,在彈性階段和塑性加載階段,理論計算得到的物理量與已有理論方法和數(shù)值模擬結(jié)果基本吻合,在卸載階段,已有理論方法不再適用,而本文理論計算得到的物理量與數(shù)值模擬結(jié)果基本吻合,驗證了該理論方法的正確性.

1 引言

無限大固體介質(zhì)中球面應(yīng)力波的傳播問題一直被關(guān)注和研究.基于彈性[1-9]和黏彈性[10-15]本構(gòu)關(guān)系,大量理論研究已經(jīng)開展,分析了球面應(yīng)力波傳播特征和規(guī)律.當載荷較高時,介質(zhì)發(fā)生塑性變形,需要利用彈塑性本構(gòu)模型對介質(zhì)進行描述,因此,開展彈塑性介質(zhì)中球面應(yīng)力波傳播理論研究具有重要的意義.

針對彈塑性球面應(yīng)力波,多位學(xué)者已經(jīng)開展了部分理論研究.Yang[16]利用特征線理論研究了線性硬化塑性材料在線性衰減的沖擊載荷作用下球面應(yīng)力波的傳播情況,分析了塑性剪切模量的影響,但對卸載區(qū)的求解借助了數(shù)值方法.Morland[17]通過假設(shè)跨越塑性波波陣面的應(yīng)力間斷的具體表達形式,限定了外加載荷的形式,在此條件下給出了彈塑性應(yīng)力波參量的解析表達式.Milne 等[18,19]研究了恒定沖擊載荷、連續(xù)卸載沖擊載荷和漸增沖擊載荷作用下彈塑性球面應(yīng)力波的傳播情況.Rapoport 等[20]對理想彈塑性介質(zhì)中球面彈塑性波的傳播進行了研究,針對球形空腔壁面施加恒定速度邊界條件和恒定壓力邊界條件兩種情況,給出了彈塑性波參量的解析表達式.Santos等[21]考慮了黏塑性Perzyna 型材料中空腔在恒定載荷作用下的動態(tài)響應(yīng),給出了介質(zhì)中應(yīng)力的計算公式并與有限元模擬結(jié)果進行了對比.肖建華等[22]對點爆炸震源產(chǎn)生的地震子波進行了求解.Chen等[23]研究了爆炸載荷加載下的彈塑性球面應(yīng)力波問題,通過求解波動方程,給出了彈性區(qū)和塑性區(qū)應(yīng)力波參量的解析解,得到了彈塑性邊界上的應(yīng)力時間歷程.

目前針對彈塑性介質(zhì)中球面應(yīng)力波的理論研究,都沒有考慮載荷卸載速率不同情況下應(yīng)力波的傳播情況.事實上,外加載荷卸載速率會影響彈塑性球面應(yīng)力波傳播情況,載荷卸載速率不同情況下應(yīng)力波傳播過程也存在區(qū)別.利用應(yīng)力波理論對不同載荷卸載速率情況下的彈塑性球面應(yīng)力波傳播過程進行分析,可以加深對應(yīng)力波傳播衰減機理的認識,便于進一步分析彈塑性球面應(yīng)力波的傳播規(guī)律.此外,現(xiàn)有理論研究缺少對介質(zhì)進入塑性后卸載區(qū)域的計算方法,也需要進一步補充研究.

本文基于線性硬化塑性本構(gòu)模型,結(jié)合球面波波動方程,分析了不同卸載速率的沖擊載荷加載下彈塑性球面應(yīng)力波的傳播過程,并通過解析求解波動方程,給出了應(yīng)力波參量的精確解.通過與現(xiàn)有理論方法、數(shù)值模擬結(jié)果進行對比,驗證了本文理論模型和方法的正確性,并對理論方法的適用性和計算過程中引入假設(shè)造成的誤差進行了分析.

2 彈塑性球面應(yīng)力波場的理論計算方法

2.1 物理模型

設(shè)無限大的固體介質(zhì)中有一個半徑為r0的球形空腔,沖擊載荷均勻作用于空腔內(nèi)壁,t時刻空腔內(nèi)壁所受壓力為p(t).其中p(t)在t=0 時刻從0 突躍至峰值P0,隨后在t＞ 0 區(qū)間內(nèi)逐漸減小.

在球?qū)ΨQ應(yīng)力應(yīng)變狀態(tài)下,固體材料容變律可以寫為

其中σr和σθ分別為徑向應(yīng)力和環(huán)向應(yīng)力,εr和εθ分別為徑向應(yīng)變和環(huán)向應(yīng)變,λ和μ為拉梅常數(shù).畸變律采用線性硬化塑性模型:

此處忽略了反向塑性加載.式中θ0為到達初始屈服強度時的偏應(yīng)變;Y=2μθ0為屈服強度;θm為最大偏應(yīng)變;Gp為塑性剪切模量.畸變律曲線如圖1所示.

圖1 畸變律曲線Fig.1.Curve of the distortion law.

在小變形條件下,球面波傳播的波動方程為

其中,r和t分別為空間徑向坐標和時間坐標,u為介質(zhì)位移,ρ為介質(zhì)密度.幾何方程為

模型的初始條件和邊界條件分別為

2.2 載荷卸載速率對應(yīng)力波傳播的影響

對上述模型進行解析求解的難點在于,介質(zhì)本構(gòu)方程中畸變律在不同階段分別滿足方程(2)—(4),只有確定了介質(zhì)處于哪一階段,才能代入對應(yīng)的畸變律方程進行求解.

根據(jù)介質(zhì)本構(gòu)關(guān)系,彈性球面應(yīng)力波和塑性球面應(yīng)力波的傳播速度均為常數(shù).由于Gp＜μ,所以塑性波波速小于彈性波波速,在彈塑性球面波傳播過程中,強間斷彈性波波陣面和強間斷塑性波波陣面在介質(zhì)中依次向前傳播.盡管載荷是連續(xù)卸載的形式,但外加載荷對應(yīng)于邊界處的徑向應(yīng)力,而介質(zhì)的加卸載條件考察的是偏應(yīng)力,因此介質(zhì)中塑性波波陣面后仍可能有一段塑性加載過程,此加載過程的長短與外加載荷卸載速率相關(guān).

假設(shè)介質(zhì)中任意位置處完成塑性加載過程進入卸載后,不再進入塑性加載階段,則介質(zhì)中的彈塑性球面應(yīng)力波傳播過程可以分為3 種情況,在r-t平面上表示如圖2 所示.

圖2 應(yīng)力波傳播示意圖 (a)載荷卸載較快;(b)載荷卸載較慢;(c)載荷卸載速率介于上述二者之間Fig.2.Schematic diagrams of the stress wave propagation (a)Rapid unloading;(b)slow unloading;(c)middle unloading.

當外加載荷卸載較快時,介質(zhì)在塑性波波陣面后直接進入卸載過程,沒有塑性加載階段,在r-t平面上的應(yīng)力波傳播過程如圖2(a)所示.圖中的粗實線表示強間斷,細實線表示弱間斷.OE為強間斷彈性波波陣面,其斜率為彈性波波速的倒數(shù),在OE前方為未擾動區(qū).OP為強間斷塑性波波陣面,其斜率為塑性波波速的倒數(shù),跨越OP后立即進入卸載階段.當強間斷塑性波波陣面?zhèn)鞑ブ罰點時,偏應(yīng)力衰減到屈服強度,塑性波完全衰減,形成以彈性波波速向前后分別傳播的兩個弱間斷PC和PD,PD在邊界發(fā)生反射,成為弱間斷DF并向前傳播.根據(jù)上述分析,區(qū)域EOr為未擾動區(qū),區(qū)域EOPP1為彈性區(qū),在塑性波波陣面OP上發(fā)生塑性加載,區(qū)域P1POt為卸載區(qū).

當外加載荷卸載較慢時,介質(zhì)在塑性波波陣面后繼續(xù)發(fā)生塑性加載行為,且塑性加載區(qū)較大,加卸載邊界不會與強間斷塑性波波陣面相交,在r-t平面上的應(yīng)力波傳播過程如圖2(b)所示.初始的應(yīng)力波傳播情況與圖2(a)相同,但跨越OP后介質(zhì)繼續(xù)發(fā)生塑性加載.當強間斷塑性波波陣面?zhèn)鞑ブ罰點時,偏應(yīng)力衰減到屈服強度,塑性波繼續(xù)以弱間斷的形式向前傳播.因為幾何擴散效應(yīng)會使?jié)M足屈服條件的塑性應(yīng)力波隨傳播而衰減,從而轉(zhuǎn)化為彈性波,所以弱間斷塑性波的傳播速度小于強間斷塑性波[24],且其傳播速度不再是常數(shù).弱間斷塑性波波陣面沿PD傳播至D點后完全衰減為彈性波.邊界在C點處開始進入卸載,并向前傳播弱間斷CG,CD為加卸載邊界.

當外加載荷卸載速率介于上述兩種情況之間時,介質(zhì)在塑性波波陣面后繼續(xù)發(fā)生塑性加載行為,但塑性加載區(qū)較小,加卸載邊界與強間斷塑性波波陣面相交,在r-t平面上的應(yīng)力波傳播過程如圖2(c)所示.此時加卸載邊界CD與塑性波波陣面相交于點D.跨越OD段后,介質(zhì)繼續(xù)進行塑性加載,跨越DP段后介質(zhì)直接進入卸載階段.

由此確定了不同情況下r-t平面上彈性區(qū)、塑性加載區(qū)和卸載區(qū)的分布情況,就可以在各區(qū)域?qū)?yīng)的本構(gòu)方程代入波動方程進行解析求解,得到應(yīng)力波傳播過程中參量的解.

2.3 應(yīng)力波參量求解

在不同區(qū)域求解對應(yīng)的應(yīng)力波波動方程,可以得到應(yīng)力波參量的解.彈性區(qū)波動方程為

其中c為彈性波波速:

引入位移勢函數(shù)?,滿足:

則波動方程(10)可化為

由于彈性前驅(qū)波中只存在右行波,上式的通解形式為

結(jié)合邊界條件可以得到f的表達式,從而計算得到介質(zhì)位移及其他應(yīng)力波參量.

塑性加載區(qū)波動方程為

其中cp為塑性波波速:

利用位移勢函數(shù)可得:

上式通解形式為

結(jié)合邊界條件可以得到f1和f2的解.

卸載區(qū)波動方程為

將位移u進行分解[16]:

其中U滿足彈性區(qū)波動方程,可以利用求解彈性區(qū)方程的方法進行求解,但此時需同時考慮左行波和右行波.g滿足方程:

結(jié)合邊界條件易得g的表達式.

利用上述方法,在應(yīng)力波傳播的r-t圖中不同區(qū)域分別求解對應(yīng)的方程,區(qū)域交界處采用位移連續(xù)條件,就可以得到全場位移的解.將位移對時間求導(dǎo)即為質(zhì)點速度.利用幾何關(guān)系(6)和(7)以及本構(gòu)關(guān)系(1)—(4),即可得到應(yīng)變和應(yīng)力的解.需要注意的是,為便于求解,本文引入了部分假設(shè),主要包括將圖2(b)中弱間斷塑性波波陣面PD和加卸載邊界CD近似為直線,忽略特征線DH和HK,以及將圖2(c)中加卸載邊界CD近似為直線,假設(shè)CD與左行特征線平行等.

3 結(jié)果與討論

本節(jié)利用彈塑性球面應(yīng)力波傳播過程分析方法和應(yīng)力波參量解析計算方法,計算彈塑性應(yīng)力波參量在傳播過程中的變化情況,并與現(xiàn)有理論方法以及動力學(xué)軟件AUTODYN 計算結(jié)果進行對比,驗證本文方法的合理性.最后對本文方法的適用性及計算過程中引入假設(shè)造成的誤差進行了分析.

3.1 恒定載荷作用下的彈塑性應(yīng)力波

Rapoport 等[20]給出了恒定沖擊載荷作用下理想彈塑性材料中應(yīng)力波參量的解析計算方法.在本文模型中,取塑性剪切模量Gp=0,外加載荷p(t)=P0為恒定值,也可以得到此情況下應(yīng)力波參量的解.同時,利用動力學(xué)軟件ANSYS AUTODYN ALE 求解器對相同的模型進行數(shù)值模擬.相關(guān)參數(shù)取值為:密度ρ=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兩種網(wǎng)格方案進行數(shù)值模擬,圖3 給出了兩種網(wǎng)格方案條件下計算得到空腔壁面處的位移時間歷程和徑向應(yīng)力的空間分布結(jié)果,兩種方案條件下的計算結(jié)果基本吻合,驗證了網(wǎng)格收斂性,后續(xù)的模擬均采用網(wǎng)格dr=2.8×10—4m 進行計算.

圖3 網(wǎng)格無關(guān)性 (a)空腔壁面位移時間歷程;(b)徑向應(yīng)力空間分布Fig.3.Gird independence:(a)Time history of the cavity surface displacement;(b)radial stress distribution.

圖4 給出了利用本文理論方法、文獻[20]中理論方法和數(shù)值模擬方法得到的不同位置處徑向應(yīng)力時間歷程和不同時刻徑向應(yīng)力的空間分布.由于載荷沒有卸載,此情況下應(yīng)力波傳播過程為圖2(b)所示的情況.Rapoport 等[20]沒有考慮弱間斷塑性應(yīng)力波的傳播和后續(xù)卸載過程,因此其方法在圖4(b),4(c),4(e)和4(f)中不再適用.本文中的方法考慮了應(yīng)力波傳播的整個過程,但在求解波動方程過程中引入了部分假設(shè),因此得到的結(jié)果與Rapoport 等[20]的理論方法以及數(shù)值模擬得到的結(jié)果有些偏差.但本文的方法適用范圍更廣,對于Rapoport 等[20]的理論方法不適用的區(qū)域,本文方法也可以求解.在r和t較小的區(qū)域,計算結(jié)果不受本文假設(shè)近似的影響,理論方法得到的結(jié)果與文獻[20]中的理論方法結(jié)果完全相同.在其他區(qū)域,和數(shù)值模擬結(jié)果能夠基本吻合,驗證了理論方法的合理性.

圖4 徑向應(yīng)力時間歷程 (a)r=0.5 m;(b)r=0.8 m;(c)r=1.4 m 和徑向應(yīng)力空間分布 (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 指數(shù)衰減載荷作用下的彈塑性應(yīng)力波

取塑性剪切模量Gp=12.5 MPa,外加載荷p(t)=P0exp(—t/T0)為指數(shù)衰減形式,其中P0為載荷峰值,T0為載荷衰減指數(shù).T0的大小會影響載荷卸載速率,進而影響應(yīng)力波傳播過程.

給定T0,首先需要判斷應(yīng)力波傳播過程屬于圖2 中的哪種情況.首先按照圖2(b)所示情況進行計算并得到C點的時間坐標tC,若tC=0,則應(yīng)力波傳播中不包含塑性加載區(qū),為圖2(a)所示的情況.若tC＞ 0,再判斷加卸載邊界CD是否與強間斷塑性波波陣面OP相交,若不相交,則為圖2(b)所示的情況,否則為圖2(c)所示的情況.

按照上述方法計算,載荷衰減指數(shù)T0分別取0.01 ms,0.2 ms 和0.03 ms 時,對應(yīng)的應(yīng)力波傳播過程分別為圖2(a)—2(c)所示的情況.圖5—圖7分別給出了衰減指數(shù)T0的3 個取值情況下,本文理論方法和數(shù)值模擬方法得到的不同位置處徑向應(yīng)力時間歷程和不同時刻徑向應(yīng)力的空間分布.

圖5 徑向應(yīng)力時間歷程 (a)r=0.3 m;(b)r=0.5 m;(c)r=0.6 m 和徑向應(yīng)力空間分布 (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 徑向應(yīng)力時間歷程 (a)r=0.4 m;(b)r=0.6 m;(c)r=0.7 m 和徑向應(yīng)力空間分布 (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).

從上述結(jié)果可以看出,3 種情況下本文理論方法計算得到的結(jié)果與數(shù)值模擬結(jié)果能夠基本吻合,驗證了應(yīng)力波傳播過程分類和參數(shù)理論計算方法的合理性.在沖擊載荷作用下,介質(zhì)中首先產(chǎn)生彈塑性雙波結(jié)構(gòu),隨著應(yīng)力波的傳播,塑性波逐漸衰減,最終強間斷完全消失,成為彈性波向前傳播.由于幾何擴散效應(yīng),在傳播過程中彈性波峰值也逐漸衰減.由于發(fā)生了塑性變形,在應(yīng)力波傳播過后,塑性區(qū)仍會存在殘余應(yīng)力,殘余應(yīng)力不會衰減,而是永久留存于介質(zhì)中.當載荷衰減指數(shù)增大時,介質(zhì)中傳播的應(yīng)力波脈寬以及殘余應(yīng)力均會增大.

3.3 理論方法適用性及誤差分析

本文的理論計算方法是在小變形假設(shè)下得到的,因此載荷峰值P0不能過高,以保證介質(zhì)滿足小變形條件.當載荷峰值P0很高時,介質(zhì)發(fā)生大變形,控制方程形式將發(fā)生改變,本文方法不再適用.

在對應(yīng)力波傳播過程分析中,本文假設(shè)介質(zhì)中任意位置處完成塑性加載過程進入卸載后,不再進入塑性加載階段,因此要求沖擊載荷P(t)在峰值后為連續(xù)卸載形式,不能出現(xiàn)二次加載.此外,對于介質(zhì)本構(gòu)模型,本文忽略了反向塑性加載階段,實際問題中介質(zhì)可能出現(xiàn)反向塑性加載,但只存在于邊界r=r0附近,因此當出現(xiàn)反向塑性加載現(xiàn)象時,本文理論方法的計算結(jié)果在邊界r=r0附近不再適用,但對其他區(qū)域仍然適用.

在滿足上述物理模型假設(shè)條件下,本文結(jié)果與數(shù)值模擬結(jié)果仍存在一定偏差,這是由于求解過程中引入的假設(shè)造成的.本文引入假設(shè)導(dǎo)致應(yīng)力波峰值后的計算結(jié)果與實際有一定的偏差,因此定義徑向應(yīng)力誤差:

其中σr為本文計算得到的徑向應(yīng)力,σr,N為數(shù)值模擬得到的徑向應(yīng)力,對于時間歷程取x=t,空間分布則取x=r,則e可以反映本文計算結(jié)果與數(shù)值模擬結(jié)果在[x1,x2]區(qū)間內(nèi)的平均相對偏差大小.對圖4—圖7 中的所有結(jié)果,按 (22)式計算,得到的誤差如表1 所示.

從表1 可以看出,大多數(shù)計算結(jié)果的相對偏差均小于10%.圖4 和圖6 中存在部分結(jié)果相對偏差較大,這兩個算例均為載荷卸載較慢的情況,對應(yīng)的應(yīng)力波傳播情況為圖2(b),由于直線化了弱間斷塑性波波陣面并忽略了卸載區(qū)弱間斷傳播的特征線DH和HK,導(dǎo)致應(yīng)力波峰值后的計算結(jié)果與實際有所偏差.圖7(e)結(jié)果相對偏差較大,這是由于在其對應(yīng)的應(yīng)力波傳播情況中,假設(shè)了加卸載邊界CD與左行特征線平行,導(dǎo)致邊界r=r0附近的計算結(jié)果與實際有所偏差.r=r0附近的計算結(jié)果與實際有所偏差.

圖6 徑向應(yīng)力時間歷程 (a)r=0.5 m;(b)r=0.75 m;(c)r=1.25 m 和徑向應(yīng)力空間分布 (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 結(jié)論

在球形空腔內(nèi)部沖擊載荷作用下,介質(zhì)中產(chǎn)生彈塑性球面應(yīng)力波并向前傳播.本文從理論分析的角度研究了彈塑性球面應(yīng)力波的傳播過程,并建立了應(yīng)力波參量的解析求解方法.

外加沖擊載荷的卸載速率對彈塑性球面應(yīng)力波傳播過程有很大影響.在r-t平面上,載荷卸載速率會影響塑性加載區(qū)的大小.當載荷卸載較快時,r-t平面上不存在塑性加載區(qū).隨著載荷卸載速率的降低,r-t平面上開始出現(xiàn)塑性加載區(qū),但塑性加載區(qū)較小時,加卸載邊界會與強間斷塑性波波陣面相交.載荷卸載速率進一步降低時,r-t平面上塑性加載區(qū)面積擴大,加卸載邊界不再與強間斷塑性波波陣面相交,因此在強間斷塑性波完全衰減后,還會以弱間斷的形式繼續(xù)傳播.

在上述應(yīng)力波傳播分析結(jié)果基礎(chǔ)上,利用解析方法對應(yīng)力波參量進行求解.在不同載荷加載條件下,利用本文方法得到應(yīng)力波參量的解并與Rapoport等[20]理論方法、數(shù)值模擬得到的結(jié)果進行對比.結(jié)果顯示,本文方法與現(xiàn)有理論方法、數(shù)值模擬方法得到的結(jié)果基本吻合,能夠反映應(yīng)力波傳播過程中的主要變化規(guī)律,且絕大多數(shù)結(jié)果相對數(shù)值模擬結(jié)果偏差均小于10%.現(xiàn)有理論方法沒有考慮載荷卸載速率對應(yīng)力波傳播的影響,也沒有給出卸載區(qū)域的求解方法,本文理論方法考慮了上述問題,相對現(xiàn)有方法具有更廣泛的適用范圍.