姚哲芳，任輝啟，沈兆武

（1.中國科學(xué)技術(shù)大學(xué)近代力學(xué)系，安徽 合肥 230027；2.西南科技大學(xué)土木工程與建筑學(xué)院，四川 綿陽 621010；3.總參工程兵科研三所，河南 洛陽 471023）

爆炸波模擬裝置是模擬爆炸沖擊波載荷環(huán)境的重要設(shè)備，對于防護(hù)工程與武器裝備毀傷評估方面的研究，具有十分重要的意義。爆室為爆炸波模擬裝置的重要部件，直接遭受炸藥爆炸，受力環(huán)境十分惡劣。由于工程的實(shí)際需要，爆室設(shè)計為一端封閉一端開口的厚壁圓筒。為使厚壁圓筒能夠在彈性范圍內(nèi)安全、重復(fù)地工作，需要深入研究它在柱狀裝藥內(nèi)爆炸條件下的強(qiáng)度問題。

對于上述問題，可參考火炮身管的設(shè)計經(jīng)驗以及爆炸容器的設(shè)計方法，而火炮發(fā)射藥為火藥，火藥與高能炸藥有較大區(qū)別。爆炸容器的研究始于1945年，中國于1963年建成了第一個爆炸容器，有關(guān)人員做了大量的研究工作［1－4］。盡管如此，由于爆炸容器設(shè)計涉及多個學(xué)科，動力學(xué)響應(yīng)過程非常復(fù)雜，迄今為止還沒有完整的爆炸容器設(shè)計標(biāo)準(zhǔn)。多數(shù)爆炸容器為兩端封閉的柱形殼體，研究主要集中在兩端封閉的柱形殼體（薄壁）在集團(tuán)裝藥內(nèi)爆炸載荷作用下的強(qiáng)度問題。而對一端封閉一端開口厚壁圓筒在柱狀裝藥爆炸條件下強(qiáng)度問題的研究未見報道，應(yīng)用于柱形殼體的一些研究方法（如等效靜載法等），不再適用于上述問題。

由于理論分析存在較大困難，實(shí)驗研究又受到場地條件、測試手段及經(jīng)費(fèi)等方面的限制，而計算機(jī)技術(shù)及數(shù)值算法有了很大發(fā)展，數(shù)值模擬成了上述問題研究的一個重要手段。本文中，主要采用數(shù)值模擬的方法，分析厚壁圓筒爆室內(nèi)柱狀裝藥爆炸非定常流場的演化過程以及筒體的動力響應(yīng)規(guī)律，同時與實(shí)驗結(jié)果及理論解進(jìn)行對比，得出一些有價值的結(jié)論，為爆室的設(shè)計提供參考。

圖1 厚壁圓筒爆室及裝藥剖面圖Fig.1 Cutaway view of the blast chamber and the charge

1 有限元計算模型

厚壁圓筒爆室及裝藥結(jié)構(gòu)如圖1所示。圓筒內(nèi)半徑RI＝170mm，外半徑RO＝340mm，筒體長度L1＝2 490mm，膛深L2＝2 198mm。TNT柱狀裝藥直徑68mm，長2 000mm，居中置放，左端位于筒底球心處，起爆點(diǎn)設(shè)置在爆室開口一端。

1.1 有限元模型及算法

由于爆室及裝藥結(jié)構(gòu)是對稱的，為了減少計算量，只取四分之一實(shí)體建模，在相關(guān)邊界上施加對稱邊界條件。

有限元模型共有三部分，流體包括炸藥和空氣，結(jié)構(gòu)為厚壁圓筒爆室。流體部分使用多物質(zhì)ALE算法，建模時流體域取為足夠大的規(guī)則六面體，全部以空氣材料劃分網(wǎng)格，網(wǎng)格尺寸為5mm×5mm×5mm，相關(guān)邊界添加無反射邊界條件，炸藥材料利用在K文件中添加關(guān)鍵字＊INITIAL＿VOLUME＿FRACTION＿GEOMETRY方式在求解過程開始前自動填充到網(wǎng)格中。結(jié)構(gòu)部分采用Lagrange算法，由于筒體球底結(jié)構(gòu)的特殊性，劃分網(wǎng)格時經(jīng)過處理使筒體單元均為計算精度較高的八節(jié)點(diǎn)六面體單元，厚壁圓筒爆室有限元模型如圖2所示。流固的耦合方式對模擬結(jié)果的準(zhǔn)確性有重要影響，本文中采用基于罰函數(shù)的耦合算法完成流體與結(jié)構(gòu)間的相互作用。有關(guān)算法可參見文獻(xiàn)［5］。

1.2 測點(diǎn)布置

如圖3所示，在靠近圓筒內(nèi)壁處布置爆炸壓力測點(diǎn)A～Y，以考察筒體內(nèi)壁受到的爆炸載荷。為了全面考察筒體的動力響應(yīng)，筒體上沿徑向布置6排測點(diǎn)，建立原點(diǎn)位于筒底球面頂點(diǎn)的坐標(biāo)系，6排測點(diǎn)的r坐標(biāo)分別為0、RI／3、2RI／3、RI、（RI＋RO）／2和RO。筒體中部（與爆炸壓力測點(diǎn)O 的z 坐標(biāo)相等）三個測點(diǎn)編號為O1、O2和O3。

圖2 厚壁圓筒爆室有限元模型Fig.2 FEM of the thick－walled cylinder

圖3 測點(diǎn)布置Fig.3 Layout of the calculated points

1.3 材料模型及參數(shù)

TNT炸藥密度1.63g／cm3，CJ爆速6 970m／s，爆壓21GPa，采用JWL狀態(tài)方程［5］，壓力p 為比體積V與比內(nèi)能E的函數(shù)

式中：A、B、R1、R2、ω為實(shí)驗確定的常數(shù)。計算中，取A＝371GPa，B＝3.23GPa，R1＝4.15，R2＝0.95，ω＝0.3，V0＝1.0，E0＝7.0GPa。

空氣密度1.29mg／cm3，采用簡化的線性多項式狀態(tài)方程［5］

式中：γ為空氣的比熱比，E為比內(nèi)能，計算中取初始比內(nèi)能0.25MPa。

厚壁圓筒材料為炮鋼，密度7.85g／cm3，彈性模量200GPa，泊松比0.33。采用簡化的Johnson－Cook材料模型及Gruneisen狀態(tài)方程。Johnson－Cook材料模型［5］定義流動應(yīng)力為

2 爆炸流場分析

計算中，TNT藥量為11.833kg，圖4給出了不同時刻爆炸流場的壓力云圖?？梢郧逦乜闯觯河捎诒óa(chǎn)物的側(cè)向飛散，先導(dǎo)沖擊波拖著斜激波向筒底運(yùn)動，當(dāng)t＝50μs時，穩(wěn)定的爆轟波結(jié)構(gòu)已經(jīng)形成。爆轟波沿藥柱向筒底傳播，從藥柱邊緣發(fā)出的斜激波（激波角約35°）打在筒壁而發(fā)生規(guī)則反射，反射激波與筒壁間形成高壓區(qū)。同時，高壓的爆轟產(chǎn)物向筒外方向運(yùn)動。當(dāng)t＝210μs時，反射激波已在軸線處相遇形成了新的高壓區(qū)，該高壓區(qū)隨著透射激波的運(yùn)動而逐漸擴(kuò)大。當(dāng)t≈300μs時，柱狀炸藥爆轟完畢，先導(dǎo)沖擊波透射入筒底的空氣中，由于沒有了爆炸反應(yīng)區(qū)提供的能量，其后的壓力迅速降低。接近t＝330μs時，先導(dǎo)沖擊波碰到筒底發(fā)生反射，因為筒底是球面，反射激波匯聚形成t＝430μs時左端的高壓區(qū)。但該高壓氣團(tuán)壓力低于筒壁面反射激波在筒軸線匯聚形成的高壓氣團(tuán)，二者的運(yùn)動速度相反。后者向側(cè)向擴(kuò)散的同時向筒底運(yùn)動，形成了強(qiáng)大的高壓氣體“射流”，最終作用于筒底部，見t＝430～500μs時的壓力云圖，這一點(diǎn)可以從筒底處所受的壓力時程曲線（見圖5）看出。圖5中，第一個峰值為先導(dǎo)沖擊波作用于筒底產(chǎn)生的最大壓力，第二個峰值為反射激波在軸線處匯聚后作用于筒底產(chǎn)生的最大壓力。當(dāng)t＝500μs后筒底的波系結(jié)構(gòu)已變得十分復(fù)雜，難于分辨。

圖6為圓筒內(nèi)壁測點(diǎn)Y、T、O、J和E的壓力時程曲線。從圖上可以看出，靠近筒口的部位（測點(diǎn)Y）受到的爆炸載荷較小。由于穩(wěn)定的爆轟波傳播，斜激波打在筒內(nèi)壁后發(fā)生反射，反射點(diǎn)勻速地從筒口掃到筒底，筒中部（測點(diǎn)T、O和J）的爆炸載荷峰值壓力幾乎一致，峰值壓力在306MPa左右。由于裝藥爆轟完畢后斜激波的強(qiáng)度降低，所以靠近筒底的部位（測點(diǎn)E）的峰值壓力小于筒中部的峰值壓力。從測點(diǎn)O的壓力時程曲線可以看出，在1000μs內(nèi)，筒內(nèi)爆炸波在該測點(diǎn)位置發(fā)生了4次反射，分別對應(yīng)圖中的4個波峰。

圖4 爆炸流場壓力云圖Fig.4 Pressure contour of blast flow field

圖5 筒底（測點(diǎn)A）壓力時程曲線Fig.5 Pressure－time curve at the bottom of the cylinder

圖6 圓筒內(nèi)壁壓力時程曲線Fig.6 Pressure－time curves on inner wall of the cylinder

3 筒體動力響應(yīng)分析

圖7 厚壁圓筒的von－Mises等效應(yīng)力云圖Fig.7 Von－Mises stress contour of the thick－walled cylinder

圖8 圓筒內(nèi)壁面上爆炸壓力峰值及沖量的分布Fig.8 The distribution of pressure peak and impulse on the inner wall of the cylinder

厚壁圓筒爆室采用Johnson－Cook材料模型和斷裂準(zhǔn)則，在整個計算過程中，所有單元均沒有失效，并且沒有產(chǎn)生塑性應(yīng)變，說明在該藥量下厚壁圓筒爆室始終在彈性范圍內(nèi)工作。圖7給出了不同時刻的厚壁圓筒爆室的von－Mises等效應(yīng)力云圖，所有單元的von－Mises等效應(yīng)力最大值均小于或等于706MPa，而實(shí)際厚壁圓筒爆室所采用的材料許可應(yīng)力在1000MPa以上，由第四強(qiáng)度理論知，厚壁圓筒爆室是安全的。

圖9 測點(diǎn)O1的應(yīng)力時程曲線Fig.9 Stress histories at point O1

圖8給出了筒內(nèi)壁爆炸壓力峰值及沖量沿軸線的分布情況。從圖上可以看出，自筒口向筒底一定距離范圍內(nèi)峰值壓力逐漸增加，穩(wěn)定爆轟形成后，筒內(nèi)壁壓力峰值相等，為306MPa。爆轟完畢后，筒壁壓力峰值逐漸減小，由于壁面反射激波在軸線匯聚后形成高壓，筒底球頂處受到的壓力峰值最大，為411MPa。

計算中得到了筒體各測點(diǎn)的應(yīng)力歷史（主要考察等效應(yīng)力），測點(diǎn)的應(yīng)力變化范圍應(yīng)該可以反映整個筒體在動力響應(yīng)過程中的應(yīng)力歷史。圖9為筒體中部測點(diǎn)O1單元的應(yīng)力歷史，可以看出，由于筒壁體內(nèi)應(yīng)力波的傳播很復(fù)雜，應(yīng)力出現(xiàn)增長現(xiàn)象，即第一個應(yīng)力峰值不是應(yīng)力的最大值。

為了直觀比較筒體各部位的等效應(yīng)力峰值，圖10給出了筒體上6排測點(diǎn)等效應(yīng)力峰值的分布情況?？梢钥闯?，危險點(diǎn)在r＝0，z＝0，即筒底球面頂點(diǎn)，等效應(yīng)力峰值達(dá)到706MPa。筒體材料許可應(yīng)力取1000MPa，由第四強(qiáng)度理論得到爆室總體動態(tài)安全系數(shù)為1000／706≈1.4。筒底封堵部分三排測點(diǎn)（r＝0，RI／3，2RI／3）的等效應(yīng)力水平相對整個筒體較高，均在550～706MPa范圍內(nèi)，說明筒底封堵部分是最危險的部位，實(shí)際使用時應(yīng)適當(dāng)保護(hù)，加入緩沖材料可以降低作用在筒底的峰值超壓，降低該處的應(yīng)力。另外，應(yīng)考慮盡可能將柱狀裝藥外移，增加裝藥與筒底的距離。

圖10 筒體等效應(yīng)力峰值的分布Fig.10 Von－Mises stress distribution of the cylinder

圖11 測點(diǎn)O1、O2和O3環(huán)向應(yīng)變時程曲線Fig.11 Hoop strain－time curves at points O1，O2and O3

圖12 筒壁環(huán)向應(yīng)變最大值的分布Fig.12 The distribution of hoop strain maximum

圖13 測點(diǎn)O1環(huán)向應(yīng)變的頻譜Fig.13 Frequency spectrum of hoop strain at point O1

圖11為筒中部測點(diǎn)O1、O2和 O3的環(huán)向應(yīng)變波形，圖12給出了筒壁環(huán)向應(yīng)變最大值沿圓筒軸線的分布情況。對圖11中的三個應(yīng)變波形進(jìn)行快速傅里葉變換分析（圖13為內(nèi)壁處測點(diǎn)O1應(yīng)變波形的頻譜），發(fā)現(xiàn)三者的頻譜基本一致，徑向運(yùn)動的能量均主要集中在f＝3.516kHz附近的一個很小的區(qū)域內(nèi)。

4 計算結(jié)果的驗證

4.1 柱狀裝藥的爆速

圖4中的爆轟波結(jié)構(gòu)類似于有限直徑柱狀裝藥的定常二維爆轟波結(jié)構(gòu)，此時的斜激波運(yùn)動速度與柱狀裝藥的爆速相同。由激波關(guān)系式可知，爆速越高，斜激波波后的壓力越高，打在筒內(nèi)壁上形成的反射激波波后的壓力也越高，即筒內(nèi)壁受到的壓力越高。另外，爆速實(shí)際上反映了炸藥能量釋放的快慢，對流場的非定常演化過程有重要影響。因此，為了準(zhǔn)確地計算爆炸載荷與筒體的相互作用，必須準(zhǔn)確地模擬柱狀裝藥的爆速。柱狀裝藥爆速D與直徑d的經(jīng)驗關(guān)系式［7］為

式中：Di為理論爆速，A和df為擬合系數(shù)。由上式計算得到柱狀裝藥爆速的經(jīng)驗值6942.4m／s。筒內(nèi)壁壓力測點(diǎn)等間距布設(shè)，由壓力峰值的到達(dá)時間差可得到數(shù)值計算的柱狀裝藥爆速6945.3m／s，與經(jīng)驗值吻合較好，說明數(shù)值計算結(jié)果是可信的。

4.2 計算結(jié)果與實(shí)測結(jié)果的比較

實(shí)驗布置如圖14所示，厚壁圓筒能夠在支撐墻套管內(nèi)滑動，封閉端頂有反后座裝置防止后座行程過大，裝藥形式、尺寸及藥量與數(shù)值模擬的情況一致。采用中北大學(xué)放入式電子測壓器DT－FDCYQ對筒內(nèi)底部靠近內(nèi)壁面處的膛壓進(jìn)行測試，沿筒體軸線距離筒口不同位置上安放美國PCB公司137A系列傳感器測試爆炸沖擊波自由場壓力。同時，在圓筒外壁面上距離開口端200mm處以及距離封閉端500mm處各對稱粘貼一對應(yīng)變片，測試環(huán)向應(yīng)變。

圖14 實(shí)驗布置圖Fig.14 Experimental layout

共進(jìn)行了10炮次，發(fā)現(xiàn)筒體內(nèi)外壁面無裂紋等損傷現(xiàn)象，內(nèi)外徑無變化，測得的數(shù)據(jù)重復(fù)性很好，說明在此工況下該厚壁圓筒爆室能夠在彈性范圍內(nèi)可安全、重復(fù)地工作。筒內(nèi)底部靠近內(nèi)壁面處典型的膛壓波形如圖15（a）所示，圖15（b）～（c）分別給出了沿筒體軸線距離筒口2.0和5.0m處的壓力時程曲線。圖16（a）～（b）分別給出了圓筒外壁面上距離開口端200mm處及距離封閉端500mm處環(huán)向應(yīng)變，其中實(shí)測應(yīng)變曲線經(jīng)過了20kHz低通濾波處理。由圖15～16可見，實(shí)測結(jié)果與數(shù)值計算結(jié)果吻合較好，證明本文中的有限元計算模型與方法是合理的。

圖15 爆炸壓力計算結(jié)果與實(shí)測結(jié)果的比較Fig.15 The comparison of blast pressure between simulation and experiment

圖16 環(huán)向應(yīng)變計算結(jié)果與實(shí)測結(jié)果的比較Fig.16 The comparison of hoop strain between simulation and experiment

4.3 軸對稱平面應(yīng)變假設(shè)條件下動力響應(yīng)的理論解

在柱狀裝藥內(nèi)爆炸作用下，厚壁圓筒在彈性范圍內(nèi)工作，實(shí)際上是一個瞬態(tài)振動問題。圓筒長度比徑向尺寸大得多，假設(shè)位于圓筒軸線上的柱狀裝藥瞬間爆轟完畢，則筒壁所受的爆炸載荷沿長度方向不變化，問題簡化為軸對稱平面應(yīng)變問題。

由彈性力學(xué)基本方程，易得厚壁圓筒徑向振動微分方程為

式中：ρ、E和μ分別為材料的密度、彈性模量和泊松比。

問題的邊界條件與初始條件分別為

文獻(xiàn)［8］中給出了位移解u（r，t），再利用彈性方程求得徑向應(yīng)力σr（r，t）、環(huán)向應(yīng)力σθ（r，t）和縱向應(yīng)力σz（r，t），進(jìn)而得到von－Mises等效應(yīng)力σe（r，t）。

從圖10可以看出，厚壁圓筒的管體部分（z／L2＞0）最大等效應(yīng)力均發(fā)生在內(nèi)壁（r＝RI），中部一段（0.3＜z／L2＜0.7）內(nèi)壁的等效應(yīng)力峰值基本相等，為570MPa左右。因此，將數(shù)值計算得到的圓筒中部測點(diǎn)O的爆炸壓力曲線代入文獻(xiàn)［8］中的動力響應(yīng)理論公式，借助Maple數(shù)學(xué)計算軟件，最終計算得到等效應(yīng)力時程曲線，最大等效應(yīng)力亦發(fā)生在內(nèi)壁，為688MPa?？梢姡谳S對稱平面應(yīng)變假設(shè)前提下，厚壁圓筒動力響應(yīng)的理論計算值偏大，可以作為保守估算。

由文獻(xiàn)［8］中式（2.12）計算得到厚壁圓筒第一振型自由振動頻率f＝3.540kHz。而由前面對應(yīng)變波形的頻譜分析，筒中部徑向運(yùn)動的能量主要集中在f＝3.516kHz附近的一個很小區(qū)域內(nèi)。二者基本一致，說明在內(nèi)柱狀裝藥爆炸條件下，筒中部的徑向運(yùn)動因受端部效應(yīng)的影響小以第一振形為主。

5 結(jié) 論

ANSYS／LS－DYNA有限元程序能夠較好地模擬厚壁圓筒爆室在柱狀裝藥爆炸作用下的動力響應(yīng)過程，可以得到爆炸流場以及動力響應(yīng)較全面的描述，彌補(bǔ)了實(shí)驗的不足。本文中的研究方法可為工程提供參考。

在本文算例裝藥條件下，厚壁圓筒爆室在彈性范圍內(nèi)工作，由于先導(dǎo)沖擊波以及筒壁面反射激波在軸線匯聚形成高壓氣團(tuán)的作用，筒底球頂處的等效應(yīng)力最大，為706MPa，總體動態(tài)安全系數(shù)為1.4，強(qiáng)度設(shè)計是安全的。

爆速實(shí)際上反映了炸藥能量釋放的快慢，對流場的非定常演化過程以及爆炸載荷與筒體的相互作用有重要影響。本文中數(shù)值計算得出柱狀裝藥的爆速為6 945.3m／s，經(jīng)驗公式得出的為6 942.4m／s，二者吻合得較好。

筒體中部環(huán)向應(yīng)變波形的頻譜分析結(jié)果證明，此處的徑向運(yùn)動受端部效應(yīng)影響小以第一振形為主。

軸對稱平面應(yīng)變假設(shè)下厚壁圓筒動力響應(yīng)的理論解可以作為保守估算。

［1］Baker W E，Hu W C L，Jackson T R.Elastic response of thin spherical shells to axisymmetric blast loading［J］.Journal of Applied Mechanics，1966，33（4）：800－806.

［2］Duffey T A，Romero C.Strain growth in spherical explosive chambers subjected to internal blast loading［J］.International Journal of Impact Engineering，2003，28（9）：967－983.

［3］胡八一，柏勁松，張明，等.真實(shí)爆炸容器殼體動力響應(yīng)的強(qiáng)度分析［J］.應(yīng)用力學(xué)學(xué)報，2001，18（3）：136－139.HU Ba－yi，BAI Jin－song，ZHANG Ming，et al.The dynamic response analysis of a real explosion－container vessel［J］.Chinese Journal of Applied Mechanics，2001，18（3）：136－139.

［4］馬圓圓，鄭津洋，陳勇軍，等.橢圓封頭圓柱形爆炸容器動力響應(yīng)的數(shù)值模擬［J］.爆炸與沖擊，2009，29（3）：249－254.MA Yuan－yuan，ZHENG Jin－yang，CHEN Yong－jun，et al.Numerical simulation on dynamic response of the cylindrical explosion containment vessel with an elliptical cover［J］.Explosion and Shock Waves，2009，29（3）：249－254.

［5］Hallquist J O.LS－DYNA theory manual［Z］.Livermore：Livermore Software Technology Corporation，2006.

［6］Rusinek A，Rodríguez－Martínez J A，Arias A，et al.Influence of conical projectile diameter on perpendicular impact of thin steel plate［J］.Engineering Fracture Mechanics，2008，75（10）：2946－2967.

［7］孫承緯，衛(wèi)玉章，周之奎.應(yīng)用爆轟物理［M］.北京：國防工業(yè)出版社，2000.

［8］林祖森.受沖擊內(nèi)壓作用的厚壁圓筒的動力分析［J］.兵工學(xué)報，1986，（1）：57－64.LIN Zu－sen.Dynamic analysis of thick cylinders subjected to internal impulsive pressure［J］.Acta Armamentarii，1986（1）：57－64.