基于AUTODYN的火工品爆炸沖擊響應(yīng)仿真分析

薛 鋒，張 剛，王 飛，韓 銘，楊澤雨

0 引 言

航天器在飛行過程中，火工品裝置爆炸會在鄰近結(jié)構(gòu)上產(chǎn)生劇烈的瞬態(tài)機械響應(yīng)[1～4]。而火工品爆炸會在結(jié)構(gòu)上產(chǎn)生劇烈的瞬態(tài)沖擊加速度響應(yīng)。這種強烈的高頻沖擊所產(chǎn)生的高沖擊、寬頻帶的爆炸沖擊環(huán)境極易造成航天器結(jié)構(gòu)和儀器設(shè)備的損傷及故障，從而導(dǎo)致設(shè)備無法正常工作，甚至造成發(fā)射失敗。因此，爆炸沖擊環(huán)境對航天器系統(tǒng)的儀器設(shè)備是不可忽略的力學(xué)環(huán)境，必須進行相應(yīng)的地面試驗驗證，考核儀器設(shè)備的工作穩(wěn)定性。但由于航天器結(jié)構(gòu)尺寸復(fù)雜且比較大，每個部段都進行真實試驗不僅費時、費力，而且成本極大。目前通常采用火工式爆炸沖擊模擬試驗來考核產(chǎn)品，從而降低試驗成本及周期。

由于火工品爆炸沖擊是一種復(fù)雜的非線性瞬態(tài)響應(yīng)過程，伴隨著高溫、高壓的物理過程，同時還涉及到炸藥與空氣、炸藥與結(jié)構(gòu)及結(jié)構(gòu)與空氣的動態(tài)相互耦合作用[4]，因此相關(guān)的研究、文獻資料少。目前研究主要通過試驗提出一些經(jīng)驗規(guī)律及通過數(shù)值計算軟件進行模擬研究，但對于爆炸載荷作用下結(jié)構(gòu)的高頻響應(yīng)計算比較困難。

AUTODYN軟件是一個顯式非線性動力分析軟件，可以解決固體、流體、和液體的動態(tài)特性及他們之間的相互作用的高度非線性動力學(xué)問題[6]。在航天、航空及軍工行業(yè)中具有重要的應(yīng)用價值。爆炸沖擊是一個高度非線性的動態(tài)問題，屬于大變形、強非線性問題，同時還涉及到流體與固體的瞬態(tài)耦合作用。在火工品爆炸沖擊中，流體材料（空氣、炸藥）采用歐拉網(wǎng)格，其它結(jié)構(gòu)（響應(yīng)板）采用拉格朗日網(wǎng)格，然后通過通用接觸設(shè)置它們之間流固耦合相互作用。本文以AUTODYN軟件為計算平臺，對火藥激勵響應(yīng)板爆炸沖擊模擬試驗裝置的沖擊加載過程及響應(yīng)進行了仿真分析，并得到了響應(yīng)板上近場、中場、遠場的變化規(guī)律，進行了響應(yīng)的試驗驗證。通過仿真計算分析及相關(guān)試驗研究，最終實現(xiàn)了高量級火工品爆炸沖擊環(huán)境。為火工品激勵響應(yīng)板爆炸沖擊模擬試驗提供了一定的指導(dǎo)。

1 仿真計算模型

1.1 火藥激勵響應(yīng)板模型

仿真計算模型以火藥激勵響應(yīng)板爆炸沖擊模擬試驗裝置試驗為基礎(chǔ)，其中響應(yīng)板為1.2 m×0.6 m×0.01 m的鋼板，材料選用08F優(yōu)質(zhì)碳素結(jié)構(gòu)鋼。根據(jù)試驗研究和分析，火藥選用 PBX9010，質(zhì)量為10.2 g，結(jié)合試驗經(jīng)計算采用直徑為18.18 mm，厚度為5.5 mm的圓柱體火藥模型。試驗裝置如圖1所示。

圖1 火藥激勵響應(yīng)板爆炸沖擊模擬試驗裝置Fig.1 The Pyrotechnically Excitation Device Simulator

1.2 網(wǎng)格密度選取

AUTODYN軟件提供了拉格朗日（Lagrange）、歐拉（Euler）、任意拉格朗日歐拉（Arbitrary Lagrangian-Eulerian，ALE）、光滑粒子流體動力（Smoothed Part icle Hydrodynamics，SPH）和多物質(zhì)流固耦合方法等多個求解器[7]。爆炸沖擊是一個高度非線性的動態(tài)問題，屬于大變形、強非線性問題，同時還涉及到空氣爆炸載荷、響應(yīng)板與工裝結(jié)構(gòu)等的相互耦合作用。由于多物質(zhì)流固耦合方法中，一些子域用歐拉法，另一些子域上用拉格朗日法，子域分離或者重疊[8]。在火工品爆炸沖擊中，流體材料（空氣、炸藥）采用歐拉和ALE算法，對其他的結(jié)構(gòu)（響應(yīng)板）采用拉格朗日算法，然后通過流固耦合方式來處理相互作用。將流固分開建模，并通過流固耦合的方式來處理相互作用，能方便的建立爆炸模型。這里我們選用多物質(zhì)流固耦合方法。

火藥激勵響應(yīng)板模型如圖 2所示，響應(yīng)板采用Lagrange方法，實際劃分的網(wǎng)格尺寸大致為5 mm×5 mm×2 mm，單元總數(shù)為36 000個，同時在模型中建立了 7個高斯觀測點進行數(shù)據(jù)的記錄，響應(yīng)板外圍與炸藥接觸處空氣采用Euler算法，Euler網(wǎng)格尺寸應(yīng)小于Lagrange網(wǎng)格，為了獲得準確的耦合相互作用，開展了多種工況的仿真計算分析：Euler和Lagrange網(wǎng)格尺寸比為1：1、1：2、1：4、1：5及1：7五種工況，經(jīng)計算分析發(fā)現(xiàn)，當Euler和Lagrange網(wǎng)格尺寸比例為1：1至1：4和1：5至1：7之間時，誤差較大，而在1：4和1：5之間時誤差較小，故Lagrange網(wǎng)格尺寸應(yīng)選為4～5倍的Euler網(wǎng)格尺寸。計算過程中發(fā)現(xiàn)Euler網(wǎng)格過大，甚至Euler網(wǎng)格尺寸大于Lagrange網(wǎng)格尺寸時，Euler網(wǎng)格會穿過Lagrange網(wǎng)格發(fā)生耦合泄露，經(jīng)過網(wǎng)格細化之后這種現(xiàn)象得到了很好的抑制。由于空氣域尺寸大細網(wǎng)格會帶來模型網(wǎng)格數(shù)量大幅度增加和計算時間步長變短，對于較大的模型，計算所需的時間耗費巨大，所以只對響應(yīng)板與炸藥接觸處的空氣網(wǎng)格加密，加密計算模型如圖3所示。

圖2 火藥激勵響應(yīng)板模型Fig.2 The Pyrotechnically Excitation Plate Models

圖3 火藥激勵響應(yīng)板流固耦合模型（1/4模型）Fig.3 The Fluid-solid Coupling Models of Pyrotechnically Excitation Plate

1.3 三維流固耦合模型

本文以火藥激勵響應(yīng)板爆炸沖擊模擬試驗裝置作為計算模型，仿真計算用AUTODYN顯式非線性動力分析軟件。根據(jù)模型的對稱性，計算模型只取1/4模型，在該模型中，響應(yīng)板結(jié)構(gòu)采用體單元而不是殼單元。

仿真計算用炸藥為 PBX9010，質(zhì)量為 10.2 g，PBX9010的配方為90%（質(zhì)量分數(shù)，下同）RDX（黑索金）和 10%Kel-F3700（三氟氯乙烯/偏氟乙烯共聚物）。炸藥、空氣采用3D多物質(zhì)Euler算法，空氣外圍施加壓力流出邊界條件來模擬無限空氣域。響應(yīng)板采用Lagrange求解器，結(jié)合火藥激勵響應(yīng)板爆炸沖擊模擬試驗裝置實際試驗情況定義響應(yīng)板沿寬度方向固支。模型采用歐拉-拉格朗日耦合作用進行計算，這里選用完全耦合進行計算。

2 材料計算模型參數(shù)

2.1 PBX9010炸藥的狀態(tài)方程

爆轟產(chǎn)物中的壓力用JWL狀態(tài)方程[8]來描述：

式中 E為單位質(zhì)量內(nèi)能；υ為比容；A，B，R1，R2，ω為常數(shù)。其中，式（1）右端第1項在高壓段起主要作用，第2項在中壓段起主要作用，第3項代表低壓段[6,9]。在爆轟產(chǎn)物膨脹的后期，方程式前 2項的作用可以忽略。PBX9010炸藥爆轟產(chǎn)物JWL狀態(tài)方程參數(shù)，如表1所示。

表1 PBX9010爆轟產(chǎn)物JWL狀態(tài)方程參數(shù)Tab.1 The Parameters of JWL Equation of State Detonation Product for PBX9010

2.2 響應(yīng)板材料模型

響應(yīng)板的材料為08F優(yōu)質(zhì)碳素結(jié)構(gòu)鋼。在本文中，08F優(yōu)質(zhì)碳素結(jié)構(gòu)鋼的強度模型采用 Johnson-Cook。該強度模型可以考慮高溫高壓下材料應(yīng)變硬化、應(yīng)變率硬化、溫度軟化的影響，適用于描述金屬材料從低應(yīng)變率到高應(yīng)變率下的動態(tài)行為[8]。

Johnson-Cook強度模型表達式為

式中 C為屈服強度；D，n為應(yīng)變的強化系數(shù)；E1為應(yīng)變率的敏感系數(shù)；npε為等效塑性應(yīng)變；*pε為塑性應(yīng)變率；m為溫度軟化參數(shù)。

2.3 空氣材料模型

空氣采用理想氣體模型，模型如下：

式中 p為氣體壓力；pA為外部環(huán)境壓力；ρ為初始空氣密度；R為氣體常數(shù)；θ為初始氣體溫度；θZ為絕對零度。計算中，初始空氣域初始溫度為300 K，密度為1.177 kg/m3,氣體常數(shù)為 287 J/（kg?K），外界空氣壓力pA為 101 340 Pa。

3 仿真計算結(jié)果分析

根據(jù)前期大量經(jīng)驗積累，使用直徑為36.36 mm，厚度為5.5 mm的圓柱體火藥激勵響應(yīng)板中心位置，其中炸藥為圓柱形，高5.5 mm，半徑為18.18 mm，爆心（Detonation）坐標為（0，0，-5.5）位于響應(yīng)板中心處，對響應(yīng)板的瞬態(tài)沖擊響應(yīng)過程進行了數(shù)值模擬。計算時，在響應(yīng)板的不同位置設(shè)置了監(jiān)測點，研究響應(yīng)板的近區(qū)、中區(qū)、遠區(qū)等不同位置的爆炸沖擊響應(yīng)及其傳遞規(guī)律，獲取沖擊環(huán)境變化規(guī)律，從而指導(dǎo)火工品爆炸沖擊試驗，優(yōu)化裝置的設(shè)計，縮短調(diào)試時間，降低試驗成本，實現(xiàn)高量級、高頻率火工品爆炸沖擊環(huán)境的可控性及工程實用性。響應(yīng)板 7個監(jiān)測點的具體位置見表2。

表2 測點位置區(qū)域劃分及特點Tab.2 Region Partition and Characteristics of Measuring Points

圖4為響應(yīng)板在t=1.998 ms時的火藥激勵位置附近的變形情況與應(yīng)力云圖，在爆炸沖擊作用下響應(yīng)板產(chǎn)生了明顯的塑性變形。在該時刻屬于爆轟波高壓階段，高頻沖擊較大，炸藥起爆后響應(yīng)板在 t=1.998 ms時，響應(yīng)板就已經(jīng)產(chǎn)生了明顯的變形。說明炸藥爆炸瞬間，爆炸物質(zhì)在有限的空間和極短的時間內(nèi)釋放出大量能量，產(chǎn)生的大量高溫、高壓產(chǎn)物迅速擠壓周圍介質(zhì)，使介質(zhì)的壓強、密度、溫度等物理量迅速上升從而形成爆炸沖擊波，爆炸沖擊波與響應(yīng)板發(fā)生相互作用，使響應(yīng)板發(fā)生塑性變形。圖5為在t=1.998 ms時的仿真模型變化情況。

圖4 響應(yīng)板應(yīng)力云圖Fig.4 The Stress Nephogram of Plate

圖5 t=3.563ms時的仿真模型變化情況Fig.5 The Changes of Simution Models When t=3.563ms

3.1 仿真計算結(jié)果數(shù)據(jù)處理分析

該模型共設(shè)置了7個Gauge點，每個監(jiān)測點距離火藥激勵位置不同，通過計算得到了近場、中場的數(shù)據(jù)。針對仿真計算數(shù)據(jù)，通過數(shù)據(jù)處理軟件Matlab編程進行了處理，得到了相應(yīng)的時域圖及沖擊響應(yīng)譜曲線圖。G1～G4的加速度時程曲線如圖 6所示。由圖 6可知，隨著測點離炸藥激勵位置距離的增大，響應(yīng)板加速度響應(yīng)幅值整體呈遞減趨勢。在炸藥爆炸瞬間，響應(yīng)板上瞬間激起一高頻、高量級沖擊信號，隨著時間的增加，該信號瞬間衰減。

圖6 仿真計算數(shù)據(jù)處理結(jié)果Fig.6 The Data Processing Results of Simulation

續(xù)圖6

3.2 仿真計算結(jié)果對比分析

針對仿真計算結(jié)果，采用Matlab編程分別對近區(qū)、中區(qū)的仿真計算結(jié)果進行了沖擊響應(yīng)譜曲線對比分析，結(jié)果如圖7～ 9所示。

圖7 測點G1、G2、G3、G4沖擊響應(yīng)譜曲線對比（近區(qū)與中區(qū)）Fig.7 The Curve Comparison of the Response Spectrum for the G1、G2、G3、G4 Points(near-field and middle field)

圖8 測點G5、G6沖擊響應(yīng)譜曲線對比（近區(qū)與中區(qū)）Fig.8 The Curve Comparison of the Response Spectrum for the G5、G6 Points (near-field and middle field)

圖9 測點G1、G5、G7沖擊響應(yīng)譜曲線對比（近區(qū)）Fig.9 The Curve Comparison of the Response Spectrum for the G1、G5、G7 Points（near-field）

由圖7～9可知，近區(qū)主要由應(yīng)力波傳播支配，其高頻響應(yīng)量級高，沖擊響應(yīng)譜拐點頻率在2000～3000 Hz之間，最大值約為25 000 g左右；隨著測點離炸藥激勵位置距離的增大，沖擊響應(yīng)譜整體呈衰減趨勢，且其拐點頻率趨勢變化很?。?000 Hz～3000 Hz之間），低頻部分由100 g左右減小至30 g，在進行試驗時可通過計算結(jié)果選擇合適的激勵點，避免低頻響應(yīng)對產(chǎn)品造成損壞。隨著距離的增加，低頻減小較明顯，沖擊響應(yīng)譜最大值呈減小趨勢。由圖9可知，近區(qū)的G1、G5、G7位置略有不同，但近區(qū)的高低頻量級、沖擊響應(yīng)譜趨勢，拐點頻率一致，從一個方面驗證了仿真計算結(jié)果的準確性。

可根據(jù)仿真計算結(jié)果。選擇合適的火藥激勵位置、火藥藥量等來實現(xiàn)不同的試驗要求，減少試驗調(diào)試時間，減低試驗成本。

4 仿真結(jié)果與試驗結(jié)果對比分析

4.1 試驗設(shè)計

設(shè)計了一套火藥激勵響應(yīng)板爆炸沖擊模擬試驗裝置（圖1）。該裝置是以柔性炸藥條激勵響應(yīng)板來模擬爆炸沖擊環(huán)境的試驗系統(tǒng)。被試產(chǎn)品安裝在響應(yīng)板的合適位置，在響應(yīng)板背面安置了柔性炸藥條，產(chǎn)品的輸入點上利用加速度傳感器監(jiān)測載荷條件。通過調(diào)節(jié)火藥激勵位置、形狀、波形發(fā)生器、響應(yīng)板邊界等來實現(xiàn)不同的沖擊響應(yīng)譜，滿足不同的試驗要求。

4.2 仿真計算與試驗結(jié)果對比分析

如上所述，對于高壓、高頻瞬態(tài)爆炸與結(jié)構(gòu)動態(tài)耦合作用的計算比較困難，對結(jié)構(gòu)爆炸沖擊響應(yīng)高精度計算就更加困難。因此，根據(jù)計算進行了相應(yīng)的試驗研究，應(yīng)用高頻響、高量級沖擊傳感器對響應(yīng)板上特征點響應(yīng)過程進行記錄，沖擊傳感器方向為軸向，并將試驗結(jié)果與數(shù)值計算結(jié)果進行了對比，結(jié)果如下圖10示。圖10為某產(chǎn)品的爆炸沖擊模擬試驗結(jié)果，該試驗的激勵點位置與仿真計算中的G7一致，試驗中火藥形狀為圓環(huán)形，導(dǎo)爆索使用80RDXXHV，藥量為10.2 g，仿真計算時的火藥等效為圓柱狀，藥量與試驗使用相同。仿真計算結(jié)果與試驗結(jié)果的沖擊響應(yīng)譜型整體趨勢基本一致，最高量級約為20 000 g，拐點頻率在1000 Hz～2000 Hz之間和4000 Hz左右，驗證了仿真計算結(jié)果的正確性。

通過仿真計算可以得到不同火藥形狀、藥量，激勵位置，不同區(qū)域的沖擊響應(yīng)譜及拐點頻率，在進行試驗時可根據(jù)不同的試驗要求，選擇合適的火藥激勵位置、火藥形狀及藥量、響應(yīng)板邊界等實現(xiàn)爆炸沖擊模擬試驗。

圖10 仿真計算與試驗結(jié)果響應(yīng)譜對比Fig.10 The Response Spectrum Comparison of Simulation and Test

5 結(jié) 論

本文利用AUTODYN接觸響應(yīng)板柔性導(dǎo)爆索爆炸沖擊響應(yīng)問題進行了數(shù)值模擬及試驗研究，得出了以下結(jié)論：

a）計算過程中發(fā)現(xiàn)歐拉網(wǎng)格過大，甚至歐拉網(wǎng)格尺寸大于拉格朗日網(wǎng)格尺寸時，歐拉網(wǎng)格會穿過朗格朗日網(wǎng)格發(fā)生耦合泄露，Euler網(wǎng)格尺寸與 Lagrange網(wǎng)格尺寸比例為1：4至1：5之間時，誤差最小，過大或者過小都會造成誤差過大，經(jīng)過網(wǎng)格細化之后這種現(xiàn)象得到了很好的抑制。

b）近區(qū)，隨著測點離炸藥激勵位置距離的增大，低頻減小較明顯，沖擊響應(yīng)譜最大值呈減小趨勢，高低頻量級、沖擊響應(yīng)譜趨勢，拐點頻率一致。

c）從近區(qū)到中區(qū)，隨著測點與激勵點距離的增大，沖擊譜量級明顯減小，拐點頻率呈減小趨勢，且低頻減小趨勢明顯。

d）通過仿真計算結(jié)果和試驗結(jié)果的對比，發(fā)現(xiàn)其沖擊譜量級及不同位置測點的沖擊譜衰減趨勢基本一致，仿真計算模擬可以得到多個測點的變化規(guī)律，可以對不同火藥量、火藥形狀、激勵位置、響應(yīng)板尺寸及邊界等進行數(shù)值模擬，選擇合適的激勵位置及測點位置，減少調(diào)試時間，降低試驗成本。