爆炸下球殼變形空間周期分布的理論計算方法*

劉文祥，張德志，鐘方平，程 帥，張慶明

（1. 西北核技術(shù)研究院強動載與效應(yīng)實驗室，陜西 西安 710024；2. 北京理工大學(xué)爆炸科學(xué)與技術(shù)國家重點實驗室，北京 100081）

1976 年Buzukov[1]實驗發(fā)現(xiàn)了爆炸容器彈性變形范圍內(nèi)的應(yīng)變增長現(xiàn)象，在這種現(xiàn)象中容器殼體的最大變形沒有出現(xiàn)在第一個振動周期內(nèi)，而是出現(xiàn)在后續(xù)周期內(nèi)。應(yīng)變增長現(xiàn)象對容器安全來說是不利的，因為傳統(tǒng)的容器設(shè)計方法一般以第一個周期內(nèi)的變形峰值為依據(jù)。近幾年，應(yīng)變增長現(xiàn)象的研究有不少新發(fā)現(xiàn)。Dong 等[2]、Li 等[3]發(fā)現(xiàn)殼體不穩(wěn)定性導(dǎo)致彎曲振動進而形成應(yīng)變增長現(xiàn)象，并研究了相關(guān)機理；劉文祥等[4-5]實驗觀察到殼體塑性變形下的應(yīng)變增長現(xiàn)象，獲得了應(yīng)變增長系數(shù)超過6 的數(shù)據(jù)，甚至推測存在系數(shù)超過10 的情況。

振動疊加是形成應(yīng)變增長現(xiàn)象的重要原因[2,6-7]。Duffey 等[7]基于振動力學(xué)中“拍”形成機理，結(jié)合殼體振動頻率分析結(jié)果以及出現(xiàn)應(yīng)變增長現(xiàn)象的應(yīng)變曲線的特征，推測頻率相近的振動疊加形成應(yīng)變增長現(xiàn)象。Liu 等[8]提出了對應(yīng)變增長現(xiàn)象解剖式分析的方法，針對帶擾動源球殼上的應(yīng)變增長現(xiàn)象，分離出彎曲應(yīng)變（彎曲波引起的殼體彎曲變形導(dǎo)致的應(yīng)變）和膜應(yīng)變（殼體拉伸變形導(dǎo)致的應(yīng)變），直接證明了頻率相近的膜振動和彎曲振動疊加是導(dǎo)致應(yīng)變增長現(xiàn)象的主要原因，還發(fā)現(xiàn)了球殼變形呈空間周期分布。

解剖式分析是研究振動疊加形成應(yīng)變增長現(xiàn)象的新方法，早期研究雖提出了殼體彎曲波、變形空間周期性分布的相關(guān)機理，但缺乏理論的支撐。本文中參考Timoshenko 梁的彎曲理論[9-10]，理論分析球殼上彎曲波的波速以及殼體變形空間分布的周期，并與數(shù)值仿真結(jié)果進行對比，以驗證理論計算的可靠性。

1 球殼上的彎曲波和變形空間周期分布

球形爆炸容器往往存在人孔門、觀察窗等功能部件，在爆炸作用下這些部件與球殼的運動存在差異，由擾動源處產(chǎn)生彎曲波并在球殼上傳播，這些部件被稱為擾動源。完全約束擾動源是指絕對靜止?fàn)顟B(tài)的擾動源，其導(dǎo)致的彎曲波比其他約束條件擾動源的更強，導(dǎo)致的應(yīng)變增長現(xiàn)象也更嚴(yán)重。Liu 等[8]利用數(shù)值仿真分析了內(nèi)徑523 mm、厚度3 mm 的完全約束擾動源下球殼的響應(yīng)情況。殼體模型見圖1，其中擾動源半徑為L，應(yīng)變輸出位置與旋轉(zhuǎn)對稱軸的距離為r′，相應(yīng)的夾角為α，應(yīng)變輸出位置與擾動源正對位置的弧線長度為x。

圖2 為L=10 mm 時球殼上擾動源正對位置外壁和內(nèi)壁的應(yīng)變曲線，可見該位置出現(xiàn)了明顯的應(yīng)變增長現(xiàn)象。Liu 等[8]的研究表明，由于擾動源處產(chǎn)生的彎曲波在擾動源正對位置處匯聚，導(dǎo)致此位置出現(xiàn)了球殼上最嚴(yán)重的應(yīng)變增長現(xiàn)象。

Liu 等[8]根據(jù)薄殼變形機理（見圖3），提出了由殼體外壁應(yīng)變和內(nèi)壁應(yīng)變計算得到彎曲應(yīng)變和膜應(yīng)變的方法，其中彎曲應(yīng)變?yōu)閺澢ㄒ饸んw彎曲振動導(dǎo)致的應(yīng)變，膜應(yīng)變是殼體膜振動引起的殼體切向拉壓變形導(dǎo)致的應(yīng)變。

圖1 球殼數(shù)值模型Fig. 1 The numerical model for the spherical shell

圖2 擾動源半徑為10 mm 時球殼上擾動源正對位置的應(yīng)變曲線Fig. 2 Total strain-time curves of spherical shellat the pole opposite to the site of perturbation when the disturbance source radius is 10 mm

圖3 薄殼變形機理Fig. 3 Deformation mechanism of thin shell

圖4 為分離出來的L=10 mm 時殼體上不同位置的彎曲應(yīng)變曲線，可以觀察到彎曲波由擾動源產(chǎn)生并在球殼上的傳播。其中彎曲波A 為頻率與球殼膜振動頻率相近的第一個彎曲波，彎曲波波群的頭部為波長最短的彎曲波。由于彎曲波頻率越高，波速越快，波長短的波將逐漸超過波長長的波。根據(jù)彎曲波到達(dá)的殼體位置以及相應(yīng)時刻，可計算出彎曲波A 以及最短彎曲波的傳播速度，如表1所示，其中彎曲波A 的傳播速度為420～450 m/s，最短彎曲波的傳播速度為2 900～3 200 m/s。

圖4 殼體上不同位置的彎曲應(yīng)變曲線和不同時刻的殼體上彎曲波Fig. 4 Bending strain-time curves at different positions on spherical shell and bend waves at different moments

表1 彎曲波的速度Table 1 Velocities of bending waves propagating along shell.

Liu 等[8]研究發(fā)現(xiàn)，彎曲波運動到擾動源正對位置后反射，而反射返回的彎曲波與入射的彎曲波由于頻率相近，兩者相互疊加將形成“駐波”現(xiàn)象，具備“駐波”的彎曲振動與膜振動疊加導(dǎo)致球殼上應(yīng)變呈周期性空間分布。圖5 為L=10，62.5 mm 時球殼上應(yīng)變峰值的空間分布情況，可見兩種情況下離擾動源正對位置較近時空間分布周期約為41 mm，離擾動源正對位置較遠(yuǎn)時空間分布周期約為82 mm。

2 理論分析

Timoshenko 在研究梁彎曲時考慮了剪切效應(yīng)，推導(dǎo)了梁彎曲波波長與波速之間的關(guān)系[9-10]。本文中將參考相關(guān)理論，推導(dǎo)球殼彎曲波波長和波速的關(guān)系。取圖1 中球殼的微元段dx 為對象，其彎曲變形如圖6 所示，Q 表示作用在任一截面的切力，M 表示相應(yīng)的彎矩，α 為球殼的微元段轉(zhuǎn)角。需要說明的是，球殼微元除了受到外力Q、M，還受到切向拉壓力以及垂直于紙面的拉壓力的作用，但這兩個外力對微元彎曲變形沒有貢獻，因此未顯示在圖中。

圖5 2 種擾動源半徑下殼體應(yīng)變峰值的分布Fig. 5 Total strain along spherical shell for two kinds of disturbance source radii

圖6 殼體彎曲變形示意圖Fig. 6 Schematic diagram of spherical shell bending deformation

進行了以下假設(shè)和限制，將問題簡化為可以進行較易求解的數(shù)學(xué)模型：

(1)彎曲波理論分析是基于平截面假定的，即假設(shè)與圖6 中虛線垂直的平截面在變形后仍為平面，并保持和變形后的虛線垂直。該假設(shè)可把平面彎曲問題簡化為一維問題。

(2)殼體發(fā)生較小的彎曲變形，即tan α 遠(yuǎn)小于1。該限制可以將后續(xù)分析的某些方程簡化為線性方程。

按照動量定理和動量矩定理，可得到兩個動力學(xué)公式：

式中：ρ0為變形前殼體的密度，A0為變形前殼體截面面積，I 為截面對中性軸的轉(zhuǎn)動慣量，V 為殼體的橫向移動速度，ω 為截面轉(zhuǎn)動的角速度。

在平截面假定下，可得到：式中：w 為殼體的橫向位移，K 為曲率。

當(dāng)殼體僅發(fā)生較小的彎曲變形時，即tan α 遠(yuǎn)小于1，可得到：

另外，彎矩M 按其定義是截面上法應(yīng)力對中性軸的合力矩，為：

式中：E 為材料彈性模量。

Timoshenko 認(rèn)為，當(dāng)彎曲波波長較短時，必須考慮剪切效應(yīng)。對彎曲變形考慮剪切應(yīng)力時，原來的平截面將因為剪切應(yīng)變而發(fā)生撓曲，即殼體的橫向位移w 實際上是由于彎矩M 作用下的轉(zhuǎn)角α 所對應(yīng)的wM和由于切力Q 作用下的剪切應(yīng)變γ 所對應(yīng)的wQ兩部分組成：

切力Q 與切應(yīng)變γ 存在關(guān)系：

式中：G 為剪切彈性模量，η 為殼體截面的形狀系數(shù)。矩形面η 的計算公式為[11]：

泊松比υ=0.29，則η=0.85。

把式（3）和（7）～(11）代入式（1）和式(2)，則得到：

兩式消去含wM的項，得到四階偏微分方程：

式（17）的通解形式為：

式中：D 為彎曲波振幅， ψ =2π/λ 為波數(shù)， p=ψC 為圓頻率， λ 為波長，C 為波速。則得到：

式（19）存在兩個解，即：

當(dāng)λ 無限小時，式（20）中右邊根號內(nèi)取正號“+”時，得到C=C0；右邊根號內(nèi)取負(fù)號“-”時，得到C=CQ。顯然，C=CQ才是真解，因此式（19）的真解取：

對于矩形截面，旋轉(zhuǎn)半徑計算公式為[12]：

式（23）給出了球殼上彎曲波波速和波長的關(guān)系。

3 理論計算結(jié)果

由公式（23）得到球殼上彎曲波波速與彎曲波波長的關(guān)系曲線，見圖7。

圖7 球殼上彎曲波波速與波長的關(guān)系Fig. 7 Relation between length and velocity of bending wave

可見彎曲波的波長越短，波速越快；當(dāng)波長λ 無限短時，彎曲波波速趨于極限值。取彈性模量E=200 GPa，密度ρ=78 00 kg/m3，泊松比υ=0.29，其計算為：

可見，最短彎曲波的波速為聲速的0.574 倍，約為2 895 m/s，與表1 中由數(shù)值仿真得到的結(jié)果2 965～3 188 m/s 相近。

另外，波的波速和波長還存在關(guān)系：

式中：f 為波的頻率。

殼體膜振動主頻的表達(dá)式為：

彎曲波A 為膜振動頻率相近的第一個彎曲波，其頻率近似取為殼體膜振動主頻，取球殼半徑a=261.5 mm，即可計算得到彎曲波A 的頻率約為5 173 Hz。

由式（23）和式（25）均可繪制出彎曲波A 的波速和波長之間的關(guān)系曲線，如圖8 所示，兩曲線的交點可以得到彎曲波A 的波速。

圖8 彎曲波A 的波長與波速Fig. 8 Length and velocity of bending wave A

由圖8 可知，頻率與呼吸振動頻率相近的彎曲波A 的波速約為376 m/s。而表1 中數(shù)值仿真得到的波速為427～443 m/s，計算結(jié)果與數(shù)值仿真結(jié)果相差(11.9～15.1)%。

Liu 等[8]給出了殼體變形空間周期分布的計算公式，離擾動源正對位置較近時其周期為：

離擾動源正對位置較遠(yuǎn)時其周期為:

式中：C 為波速，取彎曲波A 的計算結(jié)果376 m/s；f 為相應(yīng)的頻率，取計算結(jié)果5 173 Hz。則式(27)～(28)分別計算得到在離擾動源正對位置較近時周期為36.3 mm，離擾動源正對位置較遠(yuǎn)時其周期為72.7 mm，與數(shù)值仿真結(jié)果41、82 mm 相差約11.5%。

可見，理論分析結(jié)果與數(shù)值仿真結(jié)果存在一定的差異，該差異主要來源于兩個因素：一個因素是對理論分析的簡化處理，理論分析的兩個主要假設(shè)或者限制條件在數(shù)值仿真中并不完全成立；另一個因素是對數(shù)值仿真數(shù)據(jù)的讀取方法，計算彎曲波波速要讀取彎曲波在球殼上運動的距離和時間，殼體變形空間分布周期要讀取球殼上每個周期的長度，這些數(shù)據(jù)的讀取方法本身必然帶來誤差。

4 結(jié) 論

參考Timoshenko 梁的彎曲理論，基于平截面假定和殼體發(fā)生較小的彎曲變形的假設(shè)，推導(dǎo)出球殼上彎曲波波速和波長的關(guān)系，計算得到了最短彎曲波和與膜振動頻率相近的彎曲波的波速，再結(jié)合早期研究提出的殼體變形分布周期與彎曲波波速的關(guān)系，計算得到了殼體變形空間分布的周期。主要結(jié)論有：

（1）理論計算結(jié)果與數(shù)值仿真結(jié)果基本吻合，其中彎曲波波速的計算結(jié)果與數(shù)值仿真結(jié)果相差在15%以內(nèi)，殼體變形空間分布周期的計算結(jié)果與數(shù)值仿真結(jié)果相差在12%以內(nèi)。

（2）彎曲波的波長越短，波速越快。當(dāng)波長無限短時，波速趨于一極限值，約為聲速的0.574 倍。

本文的理論研究為解剖式分析應(yīng)變增長現(xiàn)象的新方法提供了一定的理論依據(jù)，進一步驗證了該分析方法的合理性。