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

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

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

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

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

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

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

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

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

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

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

圖4 殼體上不同位置的彎曲應變曲線和不同時刻的殼體上彎曲波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)象，具備“駐波”的彎曲振動與膜振動疊加導致球殼上應變呈周期性空間分布。圖5 為L=10，62.5 mm 時球殼上應變峰值的空間分布情況，可見兩種情況下離擾動源正對位置較近時空間分布周期約為41 mm，離擾動源正對位置較遠時空間分布周期約為82 mm。

2 理論分析

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

圖5 2 種擾動源半徑下殼體應變峰值的分布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ù)學模型：

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

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

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

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

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

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

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

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

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

切力Q 與切應變γ 存在關系：

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

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

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

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

式（17）的通解形式為：

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

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

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

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

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

3 理論計算結果

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

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

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

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

另外，波的波速和波長還存在關系：

式中：f 為波的頻率。

殼體膜振動主頻的表達式為：

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

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

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

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

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

離擾動源正對位置較遠時其周期為:

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

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

4 結 論

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

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

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

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