張 歡，劉天雄，李長(zhǎng)江，向樹紅，張慶明

（1.北京空間飛行器總體設(shè)計(jì)部，北京 100094；2.北京衛(wèi)星環(huán)境工程研究所，北京 100094；3.北京理工大學(xué) 爆炸科學(xué)與技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100081）

0 引言

火工裝置是由作動(dòng)器、裝藥和功能機(jī)構(gòu)組成，利用裝藥燃燒或爆炸產(chǎn)生的載荷通過功能機(jī)構(gòu)來完成特定功能的裝置[1]。航天器火工裝置在完成連接與釋放、切割與破碎、閥門打開與關(guān)閉等動(dòng)作時(shí)會(huì)產(chǎn)生高量級(jí)、寬頻帶、短時(shí)間的復(fù)雜爆炸沖擊環(huán)境，從而對(duì)航天器電子儀器、脆性材料、輕薄結(jié)構(gòu) 等產(chǎn)生破壞作用。典型火工沖擊故障模式包括晶體、陶瓷、環(huán)氧材料、玻璃封裝材料、焊點(diǎn)、電纜引線頭等的破裂、密封失效、污染粒子擴(kuò)散、繼電器和開關(guān)的抖動(dòng)及切換、微電子芯片結(jié)構(gòu)變形等[2-3]。

火工分離螺母用于衛(wèi)星與運(yùn)載火箭上面級(jí)的連接。本文使用LS-DYNA 程序?qū)鸸し蛛x螺母分離過程開展數(shù)值仿真分析，重點(diǎn)關(guān)注分離螺母解鎖時(shí)序的正確性、應(yīng)力波在結(jié)構(gòu)材料中的傳播過程以及關(guān)鍵位置的結(jié)構(gòu)加速度響應(yīng)和載荷輸出，以便對(duì)該螺母的試驗(yàn)結(jié)果做出預(yù)示。該仿真可作為提高火工分離螺母可靠性和開展優(yōu)化設(shè)計(jì)等工作的參考，其仿真結(jié)果可為航天器抗火工沖擊防護(hù)設(shè)計(jì)提供更加真實(shí)的載荷輸入，對(duì)縮短設(shè)計(jì)周期、提高分離方案有效性有重要意義。

1 數(shù)值仿真模擬算法

1.1 有限元控制方程

在動(dòng)力學(xué)有限元分析中，系統(tǒng)的控制方程可描述為[2]

式中：M為總質(zhì)量矩陣；C為阻尼系數(shù)矩陣；K為總剛度矩陣；f為節(jié)點(diǎn)載荷向量。

在顯式動(dòng)力學(xué)理論中，通常采用直接積分法中的中心差分格式對(duì)運(yùn)動(dòng)控制方程進(jìn)行積分。在中心差分格式中的加速度和速度可用位移表示為

式中：ut是t時(shí)刻的位移；ut-Δt是t-Δt時(shí)刻的位移；ut+Δt是t+Δt時(shí)刻的位移。將式(2)和式(3)代入式(1)可得用于求解各個(gè)離散時(shí)間點(diǎn)的位移的遞推公式：

在此遞推算法中存在一個(gè)起步問題，當(dāng)t=0時(shí)，為了計(jì)算uΔt，除了初始條件u0已知外，還需要知道u-Δt。由式(2)、式(3)可得

式中，(0)u˙ 可以從給定的初始條件中得到，而 (0)u˙ 則可以利用t=0 時(shí)的式(1)得到。

應(yīng)該注意到中心差分格式是條件穩(wěn)定算法，即增量步長(zhǎng)Δt必須小于某個(gè)臨界值Δtcr，否則算法將是不穩(wěn)定的。中心差分格式的解的穩(wěn)定性條件是[3]

式中Tn為有限元系統(tǒng)的最小固有振動(dòng)周期。

1.2 算法選擇

LS-DYNA 以Lagrange 算法為主，兼有Euler和ALE（Arbitrary Lagrange Eulerian）算法[4-5]。作為L(zhǎng)S-DYNA的主要算法，Lagrange算法用Lagrange坐標(biāo)系來描述物體變形，點(diǎn)的坐標(biāo)固結(jié)在變形體的內(nèi)部，并跟蹤質(zhì)點(diǎn)的運(yùn)動(dòng)，它的質(zhì)量自動(dòng)守恒，能夠清晰地顯示求解區(qū)域內(nèi)部多種物質(zhì)的界面和自由界面，明確地定義及直觀地處理邊界條件；但對(duì)于大變形情況，網(wǎng)格可能發(fā)生嚴(yán)重畸變而導(dǎo)致計(jì)算終止。Euler 算法用Euler 坐標(biāo)系描述物體運(yùn)動(dòng)，坐標(biāo)系固定在空間里，當(dāng)物體運(yùn)動(dòng)變形時(shí)，坐標(biāo)系不變，適用于描述材料有大扭曲變形的問題，然而它不顯式地描述接觸面和材料邊界，因此對(duì)各類固體邊界和接觸面的定義很不方便，也不便于描述復(fù)雜的材料本構(gòu)關(guān)系。ALE 算法是為解決流體問題而在吸收Lagrange 算法和Euler 算法的基礎(chǔ)上發(fā)展起來的一種混合算法，它可以克服單元嚴(yán)重畸變引起的數(shù)值計(jì)算困難，很好地處理整個(gè)物體發(fā)生空間位移及本身發(fā)生大變形的問題，并實(shí)現(xiàn)流-固耦合的動(dòng)態(tài)分析。由于火工分離螺母包含火工單元，若采用Lagrange 網(wǎng)格，炸藥及空氣在數(shù)值模擬的過程中會(huì)產(chǎn)生大的變形，畸變的網(wǎng)格可能會(huì)減慢運(yùn)算速度甚至?xí)?dǎo)致計(jì)算終止[4-7]，所以本文對(duì)空氣及炸藥采用ALE 網(wǎng)格，螺栓等其他結(jié)構(gòu)采用Lagrange 網(wǎng)格，并采用*CONSTRAINED_LAGRANGE_IN_SOLID 工具完成ALE 與Lagrange 網(wǎng)格的耦合。

1.3 材料本構(gòu)及狀態(tài)方程

火工分離螺母采用帶失效模式的塑性隨動(dòng)模型[4]，彈塑性材料模型包括隨動(dòng)硬化、各向同性硬化及其組合硬化3 種模式，硬化模式由硬化參數(shù)β（0≤β≤1）控制，當(dāng)β=0 時(shí)材料為隨動(dòng)硬化，β=1時(shí)則為各向同性硬化，0＜β＜1 時(shí)為混合硬化?？刂脐P(guān)鍵字為*MAT_PLASTIC_KINEMATIC，輸入數(shù)據(jù)主要包括：密度ρ，彈性模量E，泊松比υ，初始屈服極限σ0，切線模量ET。材料參數(shù)見表1。

表1 火工分離螺母的拉格朗日域材料參數(shù)[4]Table 1 Material parameters of Lagrange region for pyroshock separation nut

1.3.1 線性狀態(tài)方程

空氣采用空白材料模型，用線性狀態(tài)方程描述，材料參數(shù)見表2[4]。

線性多項(xiàng)式狀態(tài)方程為：Pa=C0+C1μ+C2μ2+C3μ3+(C4+C5μ+C6μ2)E，其中：Pa為空氣壓力；E為單位體積內(nèi)能；μ=ρ/ρ0- 1，ρ為當(dāng)前氣體密度，ρ0為初始空氣密度；C0～C6是狀態(tài)方程參數(shù)，對(duì)于空氣，C0=C1=C2=C3=C6=0，C4=C5=γ-1，γ是比熱系數(shù)。

表2 空氣材料模型的狀態(tài)方程參數(shù)Table 2 Parameters of material and state equation of air

1.3.2 爆轟產(chǎn)物狀態(tài)方程

火工品采用高能炸藥燃燒材料模型，由JWL狀態(tài)方程描述。

爆轟產(chǎn)物狀態(tài)方程是炸藥爆轟C-J 狀態(tài)之后，爆轟產(chǎn)物系統(tǒng)中各物理量（壓力、體積、溫度等）之間的關(guān)系式，它體現(xiàn)了炸藥的做功能力。LS-DYNA 中凝聚炸藥及爆轟產(chǎn)物材料模型采用了JWL 狀態(tài)方程，控制關(guān)鍵字為*EOS_JWL。JWL狀態(tài)方程的形式為

式中：P為爆轟產(chǎn)物的壓力；V是相對(duì)比容，V=v/v0，量綱為1，v= 1/ρ是比容，v0是初始比容；e是比內(nèi)能，可由炸藥的絕對(duì)內(nèi)能es和初始比內(nèi)能e0得出，es=mCVT，e0=es/v0=ρ0mCVT，e=es/v=e0/V；A、B、R1、R2、ω這5 個(gè)常數(shù)需通過試驗(yàn)數(shù)據(jù)擬合，當(dāng)炸藥密度大于1 g/cm3時(shí)，可通過γ擬合法，而不需做試驗(yàn)。

本文采用的材料參數(shù)見表3。

表3 火工品材料C-J 參數(shù)及JWL 狀態(tài)方程參數(shù)[8]Table 3 Parameters of material C-J and JWL state equation of detonater

1.4 預(yù)緊力的施加

根據(jù)文獻(xiàn)[9]，擰緊力矩T與軸向預(yù)緊力p的換算公式分別為

式中：K是擰緊力矩系數(shù)，一般取0.2；d是螺栓直徑；F0是軸向力；S是螺栓截面積。本模型中，T= 70 N·m，d= 0.012 m。經(jīng)計(jì)算p= 258 020 759 Pa。

使用如圖1的預(yù)緊力載荷曲線，施加位置為螺栓中部截面。該曲線使預(yù)緊力得到緩慢施加而不產(chǎn)生明顯的應(yīng)力波效應(yīng)，脈寬的選取為計(jì)算開始至預(yù)緊力開始釋放的時(shí)間。

圖1 預(yù)緊力曲線Fig.1 Initial stress curve

1.5 接觸類型選擇

使用*CONTACT_SURFACE_TO_SURFACE 工具定義零件間的接觸，同時(shí)使用*CONTACT_ INTERIOR 關(guān)鍵字定義Lagrange 域網(wǎng)格的自體接觸，這樣定義的接觸可以降低因連續(xù)的兩層單元出現(xiàn)負(fù)體積而使計(jì)算停止的概率，使動(dòng)力松弛部分預(yù)緊力的計(jì)算更加穩(wěn)定。

2 有限元模型建立

仿真分兩步進(jìn)行：先建立分離螺母模型，實(shí)現(xiàn)螺母的正確解鎖；然后建立分離螺母與星箭接口局部結(jié)構(gòu)組裝模型，并輸出星箭界面載荷，為衛(wèi)星的防護(hù)設(shè)計(jì)提供載荷輸入。

2.1 火工分離螺母有限元模型

采用TRUEGRID 和LS-PREPOST 軟件作為前處理軟件，Lagrange網(wǎng)格均采用SOLID實(shí)體單元，ALE 網(wǎng)格采用SOLID_ALE 單點(diǎn)多物質(zhì)單元。單元尺寸最小0.02 cm，最大0.1 cm，單元總數(shù)454 619。劃分好的有限元模型剖面圖見圖2。

圖2 火工分離螺母有限元模型剖面圖Fig.2 Sections of pyroshock separation nut FEA mode

2.2 火工分離螺母-航天器結(jié)構(gòu)一體化建模

火工分離螺母和航天器結(jié)構(gòu)的組合體模型見圖3，其中衛(wèi)星接頭采用SOLID 單元，運(yùn)載火箭上面級(jí)為實(shí)體SOLID 單元加SHELL 殼單元，SHELL 厚度為2 mm，連接界面采用共節(jié)點(diǎn)。計(jì)算時(shí)提取典型位置A1、A2、A3 的加速度響應(yīng)數(shù)據(jù)。

圖3 火工分離螺母-航天器結(jié)構(gòu)一體化模型Fig.3 Integration FEA model of pyroshock separation nut and spacecraft structure

3 計(jì)算結(jié)果與分析

3.1 動(dòng)態(tài)解鎖過程

火工分離螺母的解鎖過程見圖4，解鎖過程按正確的工作時(shí)序完成，在做功元件活塞的推動(dòng)作用下鎖緊瓣發(fā)生較大變形。

圖4 火工分離螺母的動(dòng)態(tài)解鎖過程Fig.4 Dynamic separation process of pyroshock separation nut

3.2 火工品爆炸過程

圖5和圖6顯示了應(yīng)力波在衛(wèi)星接頭和運(yùn)載火箭上面級(jí)的傳播過程，這一過程在100 μs 左右開始，“波”的形狀已不太明顯，這主要是由應(yīng)力波 過界面后的反射、折射和衍射引起的。因?yàn)樾l(wèi)星接頭的剛度比運(yùn)載上面級(jí)的剛度大，所以應(yīng)力波在衛(wèi)星接頭中的傳播更快，應(yīng)力波過界面后迅速衰減。

圖5 衛(wèi)星接頭上的應(yīng)力波傳播Fig.5 Stress wave propagation in satellite tie-in material

圖6 運(yùn)載上面級(jí)上的應(yīng)力波傳播Fig.6 Stress wave propagation in upper stage of launch vehicle

提取結(jié)構(gòu)靠近分離螺母近場(chǎng)典型位置的加速度響應(yīng)曲線（見圖7）可看出：在1 ms 內(nèi)按照爆炸、預(yù)緊力釋放、結(jié)構(gòu)撞擊的時(shí)序出現(xiàn)了3 個(gè)峰值，其中爆炸過程峰值明顯，并且很快衰減；預(yù)緊力釋放過程的峰值最大，能量釋放最多；碰撞過程的峰值不明顯，能量釋放較少，與產(chǎn)品手冊(cè)中提供的數(shù)據(jù)一致。

圖7 近場(chǎng)典型位置的加速度響應(yīng)Fig.7 Acceleration response at typical position in near field

3.3 組合體計(jì)算結(jié)果

典型位置A2 的y向加速度響應(yīng)曲線（見圖8），與試驗(yàn)1 ms 數(shù)據(jù)對(duì)比，加速度響應(yīng)與試驗(yàn)數(shù)據(jù)量級(jí)相同，且按照解鎖過程爆炸、預(yù)緊力釋放及結(jié)構(gòu)撞擊的時(shí)序出現(xiàn)了3 個(gè)峰值，與試驗(yàn)情況相符。

圖8 典型位置A2 的y 向加速度響應(yīng)Fig.8 y direction acceleration response at typical position A2

提取了星箭界面衛(wèi)星接頭一側(cè)的載荷曲線（見圖9～圖11），可以看出：x向和z向的載荷曲線相似，均為沿時(shí)間軸上下振蕩；y向載荷有3 個(gè)較大的離散峰值，曲線在時(shí)間軸一側(cè)，其值均為正。3個(gè)曲線可作為星箭接口結(jié)構(gòu)火工沖擊響應(yīng)預(yù)示的輸入條件。

圖9 星箭界面x 向界面力Fig.9 x orientation load of satellite and launch vehicle interface

圖10 星箭界面y 向界面力Fig.10 y orientation load of satellite and launch vehicle interface

圖11 星箭界面z 向界面力Fig.11 z orientation load of satellite and launch vehicle interface

3.4 結(jié)構(gòu)加速度響應(yīng)預(yù)示

以星箭界面提取的力曲線作為載荷輸入，通過MPC 將星箭界面節(jié)點(diǎn)與一個(gè)節(jié)點(diǎn)連接，在該點(diǎn)施加3 向載荷，以NASTRAN SOL112 SOLVER 為求解器，對(duì)圖3中典型位置A1、A2、A3 的y向結(jié)構(gòu)加速度響應(yīng)進(jìn)行預(yù)示，并與試驗(yàn)結(jié)果進(jìn)行對(duì)比（如圖12～圖14所示），可以看出，典型位置的結(jié)構(gòu)加速度沖擊響應(yīng)與試驗(yàn)值基本吻合。預(yù)示結(jié)果滿足制定組件沖擊試驗(yàn)條件的要求，火工分離螺母動(dòng)作時(shí)產(chǎn)生的瞬態(tài)沖擊載荷可以作為火工沖擊防護(hù)設(shè)計(jì)的載荷輸入。

圖12 A1 位置y 向的加速度響應(yīng)與試驗(yàn)值對(duì)比Fig.12 Comparison between calculation and test data of acceleration response in y direction at position A1

圖13 A2 位置y 向的加速度響應(yīng)與試驗(yàn)值對(duì)比Fig.13 Comparison between calculation and test data of acceleration response in y direction at position A2

圖14 A3 位置y 向的加速度響應(yīng)與試驗(yàn)值對(duì)比Fig.14 Comparison between calculation and test data of acceleration response in y direction at position A3

4 結(jié)論

本文使用LS-DYNA對(duì)火工分離螺母的動(dòng)態(tài)分離過程進(jìn)行了數(shù)值仿真，根據(jù)仿真結(jié)果，得到以下主要結(jié)論：

1）LS-DYNA 軟件能夠很好地模擬火工分離螺母的動(dòng)作，準(zhǔn)確反映火工品爆炸、預(yù)緊力釋放及零件結(jié)構(gòu)撞擊過程。應(yīng)力波在剛度大的結(jié)構(gòu)材料中傳播更快，過界面后迅速產(chǎn)生衰減。此傳播信息可以作為沖擊環(huán)境防護(hù)設(shè)計(jì)的參考。

2）關(guān)鍵點(diǎn)的加速度響應(yīng)隨動(dòng)作時(shí)序產(chǎn)生離散的峰值，火工品爆炸、結(jié)構(gòu)預(yù)緊力釋放及零件結(jié)構(gòu)撞擊時(shí)將分別產(chǎn)生加速度響應(yīng)峰值。螺栓結(jié)構(gòu)預(yù)緊力釋放過程引起的結(jié)構(gòu)加速度響應(yīng)峰值大于爆炸應(yīng)力波傳播引起的結(jié)構(gòu)加速度響應(yīng)峰值。結(jié)構(gòu)撞擊引起的加速度響應(yīng)峰值與應(yīng)力波引起的結(jié)構(gòu)響應(yīng)峰值相當(dāng)。

3）火工分離螺母-航天器結(jié)構(gòu)一體化建模所提取的界面力曲線可以作為火工沖擊防護(hù)設(shè)計(jì)的載荷輸入；本文使用的數(shù)值仿真方法可用于航天器結(jié)構(gòu)響應(yīng)預(yù)示。

（References）

[1]張歡, 劉天雄, 李長(zhǎng)江, 等.航天器火工沖擊環(huán)境防護(hù)技術(shù)現(xiàn)狀與應(yīng)用[J].航天器工程, 2014, 23(2):104-113 Zhang Huan, Liu Tianxiong, Li Changjiang, et al.Status and application analysis of spacecraft pyroshock protecting techniques[J].Spacecraft Engineering, 2014, 23(2):104-113

[2]周杰梁, 鄭百林.基于中心差分法的纖維結(jié)構(gòu)碰撞動(dòng)力學(xué)分析[J].力學(xué)季刊, 2011, 32(3):466-472 Zhou Jieliang, Zheng Bailin.Impact dynamics analysis of fiber-composed structure based on centrai difference method[J].Chinese Quarterly of Mechanics, 2011, 32(3):466-472

[3]劉繼先, 黃光萍.沖擊響應(yīng)譜試驗(yàn)參數(shù)的設(shè)置[J].現(xiàn)代雷達(dá), 2010, 32(2):91-94 Liu Jixian, Huang Guangping.Test parameters setup of shock response spectrum[J].Modern Radar, 2010, 32(2):91-94

[4]夏曉宇.火工分離裝置工作過程性能分析[D].長(zhǎng)沙:國(guó)防科學(xué)技術(shù)大學(xué), 2009:24-25

[5]陳敏.宇航線式火工分離裝置數(shù)值模擬及優(yōu)化設(shè)計(jì)[D].北京:北京工業(yè)大學(xué), 2007:5-8

[6]白金澤.LS-DYNA3D 理論基礎(chǔ)與實(shí)例分析[M].北京:科學(xué)出版社, 2005:150-151

[7]趙海鷗.LS-DYNA 動(dòng)力分析指南[M].北京:兵器工業(yè)出版社, 2003:6-10

[8]趙錚, 陶鋼, 杜長(zhǎng)興.爆轟產(chǎn)物JWL 狀態(tài)方程應(yīng)用研究[J].高壓物理學(xué)報(bào), 2009, 23(4):277-282 Zhao Zheng, Tao Gang, Du Changxing.Application research on JWL equation of state of detonation products[J].Chinese Journal of High Pressure Physics, 2009, 23(4):277-282

[9]李志廣.鈦合金螺紋連接結(jié)構(gòu)預(yù)緊力、應(yīng)力、可靠性分析[D].長(zhǎng)沙:國(guó)防科學(xué)技術(shù)大學(xué), 2007:7

[10]LSTC.LS-DYNA keyword user’s manual[G].Version 971, 2010:428