劉 鈺，韓 峰，董 楠，陸希成，雷 鳴

（西北核技術(shù)研究所，陜西 西安710024）

爆炸容器用于封閉高能炸藥的爆炸效應(yīng)，限制爆炸沖擊波和爆炸產(chǎn)物的作用范圍，對(duì)人員和設(shè)備實(shí)現(xiàn)有效防護(hù)，便于實(shí)驗(yàn)觀測(cè)和測(cè)試的安全，并保護(hù)環(huán)境。爆炸容器的優(yōu)異性能使其廣泛應(yīng)用于爆炸加工、科學(xué)研究、特種物質(zhì)運(yùn)輸、新材料研制等多個(gè)領(lǐng)域［1］。由于爆炸容器需要承受一定當(dāng)量的爆炸沖擊波和封閉爆炸產(chǎn)物，所以應(yīng)關(guān)注爆炸容器的安全性。

當(dāng)前爆炸容器的設(shè)計(jì)主要采用工程設(shè)計(jì)、數(shù)值模擬、實(shí)驗(yàn)驗(yàn)證的思路。一般的容器設(shè)計(jì)方法，是在給定的爆炸載荷條件下，計(jì)算容器需要承載的最大載荷并增加一定的安全裕量，進(jìn)而確定容器的壁厚，最后通過(guò)實(shí)驗(yàn)評(píng)估容器的安全性。如常用的動(dòng)力系數(shù)法［2－3］，或基于塑性變形的設(shè)計(jì)方法［4］。但爆炸容器實(shí)驗(yàn)裝置是一個(gè)復(fù)雜系統(tǒng)，在研制和實(shí)驗(yàn)過(guò)程中不可避免地會(huì)受到各種隨機(jī)因素的影響，如爆炸載荷和屈服強(qiáng)度等［5－6］。上述確定性設(shè)計(jì)方法中，沒(méi)有考慮影響容器安全的不確定性因素，會(huì)使容器安全裕量設(shè)計(jì)得過(guò)大且實(shí)驗(yàn)成本較高。而基于隨機(jī)性分析的爆炸容器安全評(píng)判方法，可以分析并量化隨機(jī)因素對(duì)爆炸容器安全的影響，進(jìn)而估計(jì)爆炸容器的安全概率。

本文中，基于爆炸容器安全的概率評(píng)判方法，研究爆炸容器的安全概率區(qū)間估計(jì)問(wèn)題，并以柱形容器為例，給出計(jì)算容器安全概率區(qū)間估計(jì)的一般方法。

1 數(shù)值模擬結(jié)果的統(tǒng)計(jì)分析

1.1 實(shí)驗(yàn)對(duì)象及數(shù)值計(jì)算

在本文的數(shù)值計(jì)算中，爆炸裝置爆炸后的狀態(tài)及產(chǎn)生的力學(xué)效應(yīng)和爆炸產(chǎn)生的力學(xué)效應(yīng)與周?chē)橘|(zhì)的相互作用，被作為一個(gè)整體使用動(dòng)力學(xué)分析軟件進(jìn)行研究，爆炸及其力學(xué)效應(yīng)使用Euler方法描述，容器及其他結(jié)構(gòu)的變形采用Lagrangian方法描述，相互之間的作用采用流固耦合技術(shù)處理［7］。由于單層柱形容器結(jié)構(gòu)相對(duì)簡(jiǎn)單，有利于研究在不同當(dāng)量TNT 炸藥爆炸下變形和破壞的概率，選擇某單層柱形容器作為研究對(duì)象，其剖面結(jié)構(gòu)示意圖如圖1所示。

圖1 柱形容器結(jié)構(gòu)示意圖Fig.1Sketch of cylindrical vessel structure

炸藥的半徑和長(zhǎng)度取決于當(dāng)量，容器內(nèi)半徑為6cm、厚度1cm、半長(zhǎng)30cm。在8種不同TNT 炸藥當(dāng)量下，對(duì)質(zhì)量比內(nèi)能及切線模量進(jìn)行了拉丁超立方抽樣，獲得了100組計(jì)算參數(shù)，并計(jì)算得到了數(shù)值模擬結(jié)果，即變形量xr的隨機(jī)樣本，部分?jǐn)?shù)據(jù)的直方圖如圖2所示。

圖2 部分爆炸當(dāng)量下容器徑向變形量直方圖Fig.2 Histogram of radial distortion of the vessel with different TNT equivalentce

1.2 數(shù)值計(jì)算結(jié)果的概率分布檢驗(yàn)

從直方圖觀察，可以假設(shè)xr服從正態(tài)分布、對(duì)數(shù)正態(tài)分布、Beta分布以及Weibull分布，然后對(duì)給定的顯著性水平α，通過(guò)柯?tīng)柲缏宸驒z驗(yàn)法［8］（KS檢驗(yàn)）對(duì)徑向變形量數(shù)據(jù)進(jìn)行分布檢驗(yàn)。原假設(shè)為

式中：Xi表示樣本。

檢驗(yàn)方法是使用KS統(tǒng)計(jì)量Dn與檢驗(yàn)臨界值Dn，α進(jìn)行比較，若Dn＞Dn，α，則拒絕H0，否則接受H0。其中檢驗(yàn)統(tǒng)計(jì)量Dn是樣本順序統(tǒng)計(jì)量X（i）頻數(shù)分布與假設(shè)分布偏差的最大值，再通過(guò)計(jì)算得到臨界值Dn，α。在顯著性水平為α＝0.01下，對(duì)不同TNT 當(dāng)量WTNT徑向變形量數(shù)據(jù)進(jìn)行了分布檢驗(yàn)，如表1所示。

表1 置信度1－α ＝0.99和樣本量n ＝100時(shí)KS檢驗(yàn)結(jié)果Table 1 The results of KS test with the confidence 1－α ＝0.99 and sample size n ＝100

從表1可見(jiàn)，除了Weibull分布檢驗(yàn)P 值明顯較小外，其他3種分布都可以描述各當(dāng)量下徑向變形量數(shù)據(jù)的不確定性。但是，正態(tài)分布的檢驗(yàn)P 值是幾種分布中最大的，因此，徑向變形量數(shù)據(jù)服從正態(tài)分布的可能性最大。在此基礎(chǔ)上，可以進(jìn)一步通過(guò)Lilliefors檢驗(yàn)法［8］對(duì)徑向變形量數(shù)據(jù)進(jìn)行正態(tài)性檢驗(yàn)。原假設(shè)為

表2 置信度1－α ＝0.99和樣本量n ＝100時(shí)Lilliefors檢驗(yàn)結(jié)果Table 2 The results of Lilliefors test with the confidence 1－α ＝0.99 and sample size n ＝100

這樣，經(jīng)過(guò)KS檢驗(yàn)以及Lilliefors檢驗(yàn)，認(rèn)為徑向變形量數(shù)據(jù)是服從正態(tài)分布的。對(duì)于未知總體的概率分布估計(jì)問(wèn)題，還可以使用最大熵方法進(jìn)行估計(jì)［9］。

1.3 正態(tài)分布參數(shù)區(qū)間估計(jì)

設(shè)爆炸容器變形量xr服從N（xr，μ，σ）分布，其概率密度函數(shù)為

式中：μ 和σ 為分布參數(shù)，且σ ＞0［10］。

在考慮置信度要求下，一般使用參數(shù)的區(qū)間估計(jì)方法。對(duì)于正態(tài)分布，該總體的均值μ 和方差σ2都未知［11］，則在置信水平為1－α 下，μ 的區(qū)間估計(jì)為

σ2的區(qū)間估計(jì)為

1.4 安全概率P（D）的區(qū)間估計(jì)

對(duì)任意給定的變形量D，安全概率P（D）服從正態(tài)分布，分布參數(shù)為μ 和σ，可以通過(guò)參數(shù)μ 和σ的區(qū)間估計(jì)來(lái)分析P D（）的置信區(qū)間。

分析正態(tài)分布函數(shù)隨參數(shù)σ（或μ）的變化情況，如圖3（a）所示，如果隨機(jī)變量x ＞μ0，則概率是參數(shù)σ 的減函數(shù)；如果x ≤μ0，則概率）是參數(shù)σ 的增函數(shù)。同理可分析，圖3（b）顯示，概率P （x ，μ ，σ0）是參數(shù)μ 的減函數(shù)。

圖3 正態(tài)分布隨單參數(shù)變化規(guī)律示意圖Fig.3 Plot of normal CDF with the change of single parameter

在雙參數(shù)坐標(biāo)空間內(nèi)，對(duì)給定隨機(jī)變量x0，當(dāng)μ ＜x0時(shí)隨參數(shù)σ 增大而減小，隨參數(shù)μ 增大也減??；當(dāng)μ≥x0時(shí)，P （x0，μ，σ）隨參數(shù)σ增大而增大，隨參數(shù)μ增大仍減小，這與前面的分析結(jié)果是一致的。從而可以看出P （x ，μ，σ）在μ 和σ 的取值空間上關(guān)于x 是單調(diào)遞增的，如圖4所示。

圖4 正態(tài)分布隨雙參數(shù)變化規(guī)律示意圖Fig.4Surface plot of normal CDF with the change of both parameters

由于x 在參數(shù)μ的兩側(cè)時(shí)，P （x）與參數(shù)的單調(diào)性關(guān)系是不同的，參數(shù)μ又在區(qū)間上變化，根據(jù)x 和μ 的大小關(guān)系，可以分段給出P（x）的上、下限函數(shù)

則由式（6）～（8）可計(jì)算得到，在置信度1－α 下，安全概率P D（）的置信區(qū)間為［PLD（），PUD（）］。

P（x）的函數(shù)圖像如圖5所示。圖5（a）為不同參數(shù)組合下P D（）函數(shù)的圖像，其中aL、aU分別表示分布參數(shù)的上、下界。圖5（b）是圖5（a）中曲線的包絡(luò)，圖中標(biāo)出了在D 取值不同時(shí)，不同參數(shù)組合下P D（）函數(shù)的上、下限。同樣，也可以通過(guò)優(yōu)化方法計(jì)算P D（）的區(qū)間估計(jì)。

圖5 在置信度為1－α ＝0.95下，P（D）的上、下限函數(shù)圖Fig.5 The upper and lower limits curve of the P（D）with the confidence intervals 1－α ＝0.99

2 柱形容器安全概率的區(qū)間估計(jì)

根據(jù)以上統(tǒng)計(jì)分析方法，可以對(duì)柱形容器在不同TNT 爆炸當(dāng)量下，100組容器徑向變形量數(shù)值模擬數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析。

對(duì)于給定變形量要求D ＝15%，在置信度1－α ＝0.95 下，考察徑向變形量xr＜D 的概率的置信區(qū)間為，結(jié)果如表3所示。

表3 置信度1－α ＝0.95時(shí)安全概率P（15%）的置信區(qū)間Table 3 The confidence intervals of the security probability P（15%）corresponds to1－α ＝0.95

顯然，對(duì)于給定變形量D＝15%，在置信度1－α＝0.95下，隨爆炸TNT 當(dāng)量變化，柱形容器徑向變形量xr＜15%的概率具有明顯變化。那么可以反過(guò)來(lái)考察，如果要求在置信度1－α＝0.99下，對(duì)不同當(dāng)量TNT 爆炸，柱形容器徑向變形量xr不超過(guò)Di（i＝1，2，…，8）的概率下限PL（Di）＝0.999，則可以計(jì)算出不同當(dāng)量下的Di，如表4所示。

表4 安全概率PL（Di）＝0.999對(duì)應(yīng)的給定變形量Table 4 The radial distortion threshold corresponding to the security probability PL（Di）＝0.999

此時(shí)各當(dāng)量下，徑向變形量以最小有99.9%的概率小于Di，且這個(gè)結(jié)論的置信度為99%，這樣的結(jié)論已經(jīng)是相當(dāng)可靠的了，那么此時(shí)計(jì)算得到給定變形量的下限可以作為衡量實(shí)驗(yàn)安全的判據(jù)。

3 結(jié) 論

綜上所述，在爆炸容器設(shè)計(jì)完成后，該容器安全概率的區(qū)間估計(jì)與爆炸載荷具有一定的相關(guān)性：

（1）TNT 當(dāng)量與Di是近似線性的，說(shuō)明炸藥當(dāng)量越大，徑向變形量上限也越大，容器安全性越低，如圖6所示。

（2）隨著TNT 當(dāng)量的增大，容器的安全概率將逐漸降低，對(duì)應(yīng)的（在一定置信度下）概率上下限也逐漸降低，且分散性（不確定性）增大，如圖7所示。

圖6 TNT 當(dāng)量與徑向變形量的關(guān)系圖Fig.6 The relation between the radial distortion and TNT equivalence

圖7 TNT 當(dāng)量與安全概率的關(guān)系圖Fig.7 The relation between the security probability and TNT equivalence

本文在爆炸容器數(shù)值計(jì)算中考慮了主要計(jì)算參數(shù)的隨機(jī)性，并通過(guò)Monte－Carlo方法獲得了計(jì)算結(jié)果的隨機(jī)樣本，取代了以往單一的、確定性數(shù)值計(jì)算結(jié)果，能夠更客觀地反映真實(shí)的物理過(guò)程。應(yīng)用統(tǒng)計(jì)分析方法對(duì)柱形容器爆炸后徑向變形量的不確定性進(jìn)行了定量描述，給出了在一定置信度下容器安全概率的區(qū)間估計(jì)。

［1］ 朱文輝，薛鴻陸，韓鈞萬(wàn)，等.爆炸容器動(dòng)力學(xué)研究進(jìn)展評(píng)述［J］.力學(xué)進(jìn)展，1996，26（1）：68－78.Zhu Wen－h(huán)ui，Xue Hong－lu，Han Jun－wan，et al.The research advances in the dynamics of explosive chambers［J］.Advance in Mechanics，1996，26（1）：68－78.

［2］ 胡八一，柏勁松，劉大敏，等.爆炸容器的工程設(shè)計(jì)方法及應(yīng)用［J］.壓力容器，2000，17（2）：39－41.Hu Ba－yi，Bai Jing－song，Liu Da－min，et al.The engineering design method of explosion containment vessels and its application［J］.Pressure Vessel Technology，2000，17（2）：39－41.

［3］ 王定賢，胡昊，王萬(wàn)鵬，等.動(dòng)力系數(shù)法在爆炸容器設(shè)計(jì)中的應(yīng)用［J］.四川兵工學(xué)報(bào)，2011，32（4）：24－26.Wang Ding－xian，Hu Hao，Wang Wan－peng，et al.Application of the dynamical coefficient method in explosion vessel design［J］.Journal of Sichuan Ordnance，2011，32（4）：24－26.

［4］ 胡永樂(lè)，喻名德，崔云霄.爆炸容器塑性失穩(wěn)和變形吸能［J］.兵工學(xué)報(bào)，2009，30增刊2：28－32.Hu Yong－le，Yu Ming－de，Cui Yun－xiao.Plastic instabilities and deformation energy absorption of explosion containment vessels［J］.Acta Armamentarii，2009，30suppl 2：28－32.

［5］ 陳國(guó)華，戴樹(shù)和.概率安全評(píng)定方法的計(jì)算程序及其工程應(yīng)用［J］.南京航空航天大學(xué)學(xué)報(bào)，1999，31（3）：346－350.Chen Guo－h(huán)ua，Dai Shu－h(huán)e.Calculating procedure of probability safety assessment method and its engineering application［J］.Journal of Nanjing University of Aeronautics ＆Astronautics，1999，31（3）：346－350.

［6］ 趙建平.壓力容器概率安全評(píng)定技術(shù)進(jìn)展［J］.化工設(shè)備與管道，2001，38（2）：17－19.Zhao Jian－ping.Survey of probabilistic safety assessment technology for pressure vessels［J］.Process Equipment ＆Piping，2001，38（2）：17－19.

［7］ 浦錫鋒，王仲琦，白春華，等.用于爆炸流場(chǎng)與結(jié)構(gòu)間相互作用分析的Euler－Lagrange耦合模擬技術(shù)［J］.爆炸與沖擊，2011，31（1）：6－10.Pu Xi－feng，Wang Zhong－qi，Bai Chun－h(huán)ua，et al.An Euler－Lagrange coupling method for numerical simulation of explosion－structure interaction［J］.Explosion and Shock Waves，2011，31（1）：6－10.

［8］ 吳翊，李永樂(lè)，胡慶軍.應(yīng)用數(shù)理統(tǒng)計(jì)［M］.長(zhǎng)沙：國(guó)防科技大學(xué)出版社，2005：112－120.

［9］ 劉鈺，韓峰，王玉恒，等.一種基于密度核估計(jì)的最大熵方法［J］.工程數(shù)學(xué)學(xué)報(bào)，2011，28（3）：285－292.Liu Yu，Han Feng，Wang Yu－h(huán)eng，et al.A new maximum entropy method based on kernel density estimate［J］.Chinese Journal of Engineering Mathematics，2011，28（3）：285－292.

［10］ Casella G，Berger R L.Statistical inference［M］.Beijing：China Machine Press，2010：417－427.

［11］ 韓峰，王建國(guó)，喬登江.對(duì)數(shù)正態(tài)分布場(chǎng)合下產(chǎn)品加固性能的Bayes評(píng)估方法［J］.應(yīng)用概率統(tǒng)計(jì)，2009，25（4）：433－440.Han Feng，Wang Jian－guo，Qiao Deng－jiang.Bayesian assessment for product’s radiation hardening performance under lognormal distribution［J］.Chinese Journal of Applied Probability and Statistics，2009，25（4）：433－440.

［12］ 袁亞湘，孫文瑜.最優(yōu)化理論與方法［M］.北京：科學(xué)出版社，2001：373－391.