航天器火工分離螺母的火工沖擊環境數值仿真研究

張 歡，劉天雄，李長江，向樹紅，張慶明

（1.北京空間飛行器總體設計部，北京 100094；2.北京衛星環境工程研究所，北京 100094；3.北京理工大學 爆炸科學與技術國家重點實驗室，北京 100081）

0 引言

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

火工分離螺母用于衛星與運載火箭上面級的連接。本文使用LS-DYNA 程序對火工分離螺母分離過程開展數值仿真分析，重點關注分離螺母解鎖時序的正確性、應力波在結構材料中的傳播過程以及關鍵位置的結構加速度響應和載荷輸出，以便對該螺母的試驗結果做出預示。該仿真可作為提高火工分離螺母可靠性和開展優化設計等工作的參考，其仿真結果可為航天器抗火工沖擊防護設計提供更加真實的載荷輸入，對縮短設計周期、提高分離方案有效性有重要意義。

1 數值仿真模擬算法

1.1 有限元控制方程

在動力學有限元分析中，系統的控制方程可描述為[2]

式中：M為總質量矩陣；C為阻尼系數矩陣；K為總剛度矩陣；f為節點載荷向量。

在顯式動力學理論中，通常采用直接積分法中的中心差分格式對運動控制方程進行積分。在中心差分格式中的加速度和速度可用位移表示為

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

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

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

應該注意到中心差分格式是條件穩定算法，即增量步長Δt必須小于某個臨界值Δtcr，否則算法將是不穩定的。中心差分格式的解的穩定性條件是[3]

式中Tn為有限元系統的最小固有振動周期。

1.2 算法選擇

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

1.3 材料本構及狀態方程

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

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

1.3.1 線性狀態方程

空氣采用空白材料模型，用線性狀態方程描述，材料參數見表2[4]。

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

表2 空氣材料模型的狀態方程參數Table 2 Parameters of material and state equation of air

1.3.2 爆轟產物狀態方程

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

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

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

本文采用的材料參數見表3。

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

1.4 預緊力的施加

根據文獻[9]，擰緊力矩T與軸向預緊力p的換算公式分別為

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

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

圖1 預緊力曲線Fig.1 Initial stress curve

1.5 接觸類型選擇

使用*CONTACT_SURFACE_TO_SURFACE 工具定義零件間的接觸，同時使用*CONTACT_ INTERIOR 關鍵字定義Lagrange 域網格的自體接觸，這樣定義的接觸可以降低因連續的兩層單元出現負體積而使計算停止的概率，使動力松弛部分預緊力的計算更加穩定。

2 有限元模型建立

仿真分兩步進行：先建立分離螺母模型，實現螺母的正確解鎖；然后建立分離螺母與星箭接口局部結構組裝模型，并輸出星箭界面載荷，為衛星的防護設計提供載荷輸入。

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

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

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

2.2 火工分離螺母-航天器結構一體化建模

火工分離螺母和航天器結構的組合體模型見圖3，其中衛星接頭采用SOLID 單元，運載火箭上面級為實體SOLID 單元加SHELL 殼單元，SHELL 厚度為2 mm，連接界面采用共節點。計算時提取典型位置A1、A2、A3 的加速度響應數據。

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

3 計算結果與分析

3.1 動態解鎖過程

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

圖4 火工分離螺母的動態解鎖過程Fig.4 Dynamic separation process of pyroshock separation nut

3.2 火工品爆炸過程

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

圖5 衛星接頭上的應力波傳播Fig.5 Stress wave propagation in satellite tie-in material

圖6 運載上面級上的應力波傳播Fig.6 Stress wave propagation in upper stage of launch vehicle

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

圖7 近場典型位置的加速度響應Fig.7 Acceleration response at typical position in near field

3.3 組合體計算結果

典型位置A2 的y向加速度響應曲線（見圖8），與試驗1 ms 數據對比，加速度響應與試驗數據量級相同，且按照解鎖過程爆炸、預緊力釋放及結構撞擊的時序出現了3 個峰值，與試驗情況相符。

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

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

圖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 結構加速度響應預示

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

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

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

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

4 結論

本文使用LS-DYNA對火工分離螺母的動態分離過程進行了數值仿真，根據仿真結果，得到以下主要結論：

1）LS-DYNA 軟件能夠很好地模擬火工分離螺母的動作，準確反映火工品爆炸、預緊力釋放及零件結構撞擊過程。應力波在剛度大的結構材料中傳播更快，過界面后迅速產生衰減。此傳播信息可以作為沖擊環境防護設計的參考。

2）關鍵點的加速度響應隨動作時序產生離散的峰值，火工品爆炸、結構預緊力釋放及零件結構撞擊時將分別產生加速度響應峰值。螺栓結構預緊力釋放過程引起的結構加速度響應峰值大于爆炸應力波傳播引起的結構加速度響應峰值。結構撞擊引起的加速度響應峰值與應力波引起的結構響應峰值相當。

3）火工分離螺母-航天器結構一體化建模所提取的界面力曲線可以作為火工沖擊防護設計的載荷輸入；本文使用的數值仿真方法可用于航天器結構響應預示。

（References）

[1]張歡, 劉天雄, 李長江, 等.航天器火工沖擊環境防護技術現狀與應用[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]周杰梁, 鄭百林.基于中心差分法的纖維結構碰撞動力學分析[J].力學季刊, 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]劉繼先, 黃光萍.沖擊響應譜試驗參數的設置[J].現代雷達, 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].長沙:國防科學技術大學, 2009:24-25

[5]陳敏.宇航線式火工分離裝置數值模擬及優化設計[D].北京:北京工業大學, 2007:5-8

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

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

[8]趙錚, 陶鋼, 杜長興.爆轟產物JWL 狀態方程應用研究[J].高壓物理學報, 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]李志廣.鈦合金螺紋連接結構預緊力、應力、可靠性分析[D].長沙:國防科學技術大學, 2007:7

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