沖擊載荷下線性硬化材料中球面應力波場的理論計算方法研究

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

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

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

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

1 引言

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

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

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

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

2 彈塑性球面應力波場的理論計算方法

2.1 物理模型

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

在球對稱應力應變狀態(tài)下,固體材料容變律可以寫為

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

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

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

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

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

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

2.2 載荷卸載速率對應力波傳播的影響

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

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

假設介質中任意位置處完成塑性加載過程進入卸載后,不再進入塑性加載階段,則介質中的彈塑性球面應力波傳播過程可以分為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為強間斷彈性波波陣面,其斜率為彈性波波速的倒數(shù),在OE前方為未擾動區(qū).OP為強間斷塑性波波陣面,其斜率為塑性波波速的倒數(shù),跨越OP后立即進入卸載階段.當強間斷塑性波波陣面?zhèn)鞑ブ罰點時,偏應力衰減到屈服強度,塑性波完全衰減,形成以彈性波波速向前后分別傳播的兩個弱間斷PC和PD,PD在邊界發(fā)生反射,成為弱間斷DF并向前傳播.根據上述分析,區(qū)域EOr為未擾動區(qū),區(qū)域EOPP1為彈性區(qū),在塑性波波陣面OP上發(fā)生塑性加載,區(qū)域P1POt為卸載區(qū).

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

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

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

2.3 應力波參量求解

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

其中c為彈性波波速:

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

則波動方程(10)可化為

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

結合邊界條件可以得到f的表達式,從而計算得到介質位移及其他應力波參量.

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

其中cp為塑性波波速:

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

上式通解形式為

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

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

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

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

結合邊界條件易得g的表達式.

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

3 結果與討論

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

3.1 恒定載荷作用下的彈塑性應力波

Rapoport 等[20]給出了恒定沖擊載荷作用下理想彈塑性材料中應力波參量的解析計算方法.在本文模型中,取塑性剪切模量Gp=0,外加載荷p(t)=P0為恒定值,也可以得到此情況下應力波參量的解.同時,利用動力學軟件ANSYS AUTODYN ALE 求解器對相同的模型進行數(shù)值模擬.相關參數(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兩種網格方案進行數(shù)值模擬,圖3 給出了兩種網格方案條件下計算得到空腔壁面處的位移時間歷程和徑向應力的空間分布結果,兩種方案條件下的計算結果基本吻合,驗證了網格收斂性,后續(xù)的模擬均采用網格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]中理論方法和數(shù)值模擬方法得到的不同位置處徑向應力時間歷程和不同時刻徑向應力的空間分布.由于載荷沒有卸載,此情況下應力波傳播過程為圖2(b)所示的情況.Rapoport 等[20]沒有考慮弱間斷塑性應力波的傳播和后續(xù)卸載過程,因此其方法在圖4(b),4(c),4(e)和4(f)中不再適用.本文中的方法考慮了應力波傳播的整個過程,但在求解波動方程過程中引入了部分假設,因此得到的結果與Rapoport 等[20]的理論方法以及數(shù)值模擬得到的結果有些偏差.但本文的方法適用范圍更廣,對于Rapoport 等[20]的理論方法不適用的區(qū)域,本文方法也可以求解.在r和t較小的區(qū)域,計算結果不受本文假設近似的影響,理論方法得到的結果與文獻[20]中的理論方法結果完全相同.在其他區(qū)域,和數(shù)值模擬結果能夠基本吻合,驗證了理論方法的合理性.

圖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 指數(shù)衰減載荷作用下的彈塑性應力波

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

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

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

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

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

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

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

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

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

從表1 可以看出,大多數(shù)計算結果的相對偏差均小于10%.圖4 和圖6 中存在部分結果相對偏差較大,這兩個算例均為載荷卸載較慢的情況,對應的應力波傳播情況為圖2(b),由于直線化了弱間斷塑性波波陣面并忽略了卸載區(qū)弱間斷傳播的特征線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平面上,載荷卸載速率會影響塑性加載區(qū)的大小.當載荷卸載較快時,r-t平面上不存在塑性加載區(qū).隨著載荷卸載速率的降低,r-t平面上開始出現(xiàn)塑性加載區(qū),但塑性加載區(qū)較小時,加卸載邊界會與強間斷塑性波波陣面相交.載荷卸載速率進一步降低時,r-t平面上塑性加載區(qū)面積擴大,加卸載邊界不再與強間斷塑性波波陣面相交,因此在強間斷塑性波完全衰減后,還會以弱間斷的形式繼續(xù)傳播.

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