楊 光，劉 波，李正舉，張 程

（中國(guó)空間技術(shù)研究院 通信衛(wèi)星事業(yè)部，北京 100094）

0 引言

航天器所經(jīng)受的沖擊環(huán)境主要來(lái)源于各種火工裝置的起爆，例如星箭分離、組合體分離、伸展部件的展開(kāi)等[1]。其特點(diǎn)為持續(xù)時(shí)間短、頻率范圍廣、加速度幅值高[2]。一般來(lái)說(shuō)，火工品起爆對(duì)航天器主結(jié)構(gòu)影響較小，但其引發(fā)的沖擊環(huán)境對(duì)石英晶體、陶瓷器件、精密電子部件和繼電器等元器件的影響很大，可能會(huì)造成這些元器件失效或損壞，甚至導(dǎo)致整個(gè)航天任務(wù)的失敗[3]。因此有必要進(jìn)行航天器沖擊響應(yīng)預(yù)示和試驗(yàn)驗(yàn)證工作。

沖擊響應(yīng)預(yù)示是航天工程中沖擊環(huán)境分析與試驗(yàn)條件制定的基礎(chǔ)，其難點(diǎn)之一在于沖擊源力函數(shù)的選取。沖擊載荷的高量級(jí)以及沖擊波傳遞的復(fù)雜性，使得沖擊源附近的力載荷難以準(zhǔn)確測(cè)量[4]。因此需要建立與沖擊源量級(jí)以及時(shí)間特性一致的力函數(shù)對(duì)沖擊源進(jìn)行模擬。本文以結(jié)構(gòu)星試驗(yàn)測(cè)試數(shù)據(jù)反推火工沖擊源力函數(shù)，擬合方法可為后續(xù)相似平臺(tái)衛(wèi)星沖擊響應(yīng)預(yù)示的輸入載荷設(shè)定提供參考。

1 沖擊預(yù)示有限元法

1.1 沖擊預(yù)示有限元法的選擇

沖擊預(yù)示有限元法將實(shí)際受力結(jié)構(gòu)和火工品簡(jiǎn)化為合理的有限單元?jiǎng)恿W(xué)模型，采用數(shù)值分析手段，通過(guò)瞬態(tài)響應(yīng)分析得到結(jié)構(gòu)的加速度響應(yīng)，基于此進(jìn)一步計(jì)算沖擊響應(yīng)譜，使其能夠合理預(yù)示航天器結(jié)構(gòu)復(fù)雜的沖擊環(huán)境。

有限元法可分為隱式和顯式2種，二者的區(qū)別為：在時(shí)間步長(zhǎng)方面，隱式有限元法無(wú)條件地穩(wěn)定，與積分時(shí)間步長(zhǎng)無(wú)關(guān)，而顯式有限元法則存在穩(wěn)定性條件[5-6]，因此隱式有限元法的時(shí)間步長(zhǎng)可在ms級(jí)，而顯式有限元法的時(shí)間步長(zhǎng)需要達(dá)到μs級(jí)。但這并不意味著顯式有限元法的計(jì)算效率更低。在顯式有限元法中，系統(tǒng)的質(zhì)量矩陣不包含耦合項(xiàng)，從而能夠提高計(jì)算效率；且顯式有限元法能夠更好地處理強(qiáng)非線性問(wèn)題。

本文將選擇顯式有限元法進(jìn)行計(jì)算，使用的軟件為MSC.Dytran。

1.2 顯式有限元法控制方程

動(dòng)力學(xué)系統(tǒng)運(yùn)動(dòng)微分方程可描述為

進(jìn)一步改寫為

其中：M為結(jié)構(gòu)質(zhì)量矩陣；C為結(jié)構(gòu)阻尼矩陣；K為結(jié)構(gòu)剛度矩陣；Fn,ext為外加載荷矢量；Fn,int為內(nèi)力矢量。

由式(3)可以看出，計(jì)算加速度需要進(jìn)行矩陣求逆，而在顯式有限元法的計(jì)算中，系數(shù)矩陣僅為結(jié)構(gòu)質(zhì)量矩陣，不包含耦合項(xiàng)，因此無(wú)須進(jìn)行矩陣分解。

采用中心差分法進(jìn)行時(shí)間迭代：

中心差分法需要穩(wěn)定性條件[7]

其中Tn為系統(tǒng)最小固有周期，即時(shí)間步長(zhǎng)Δt必須小于某一數(shù)值，算法才能穩(wěn)定。

2 沖擊預(yù)示有限元簡(jiǎn)化模型

2.1 衛(wèi)星太陽(yáng)電池陣解鎖試驗(yàn)

本文以某平臺(tái)工程星太陽(yáng)電池陣解鎖過(guò)程為研究對(duì)象，使用試驗(yàn)實(shí)測(cè)數(shù)據(jù)反推火工沖擊源模擬函數(shù)，其中試驗(yàn)實(shí)測(cè)數(shù)據(jù)為地面整星太陽(yáng)電池陣展開(kāi)時(shí)的測(cè)試數(shù)據(jù)。衛(wèi)星兩側(cè)的太陽(yáng)電池陣分別通過(guò)8個(gè)火工品裝置固定在衛(wèi)星艙板上，展開(kāi)時(shí)分5次起爆，由分布在火工品壓緊點(diǎn)附近的沖擊加速度傳感器測(cè)量并記錄火工品起爆時(shí)產(chǎn)生的加速度。8個(gè)火工品壓緊點(diǎn)分布如圖1所示，其中的黑色矩形塊即為壓緊點(diǎn)。展開(kāi)過(guò)程中太陽(yáng)電池陣板面垂直于地面放置，并通過(guò)吊索懸掛模擬失重狀態(tài)，如圖2所示。

圖1 衛(wèi)星艙板上火工品壓緊點(diǎn)布局Fig. 1 The layout of initiating explosive devices on satellite board

圖2 太陽(yáng)電池陣解鎖試驗(yàn)狀態(tài)Fig. 2 Schematic diagram of solar array deployment test

2.2 有限元模型

衛(wèi)星的艙板多采用蜂窩夾層結(jié)構(gòu)板，由面板、蜂窩芯子和板芯膠組成。為避免拉彎耦合效應(yīng)以及固化后翹曲變形，上、下面板一般選取相同的材料和厚度[8]。根據(jù)衛(wèi)星通信艙南北板材料屬性，建立艙板模型，其中上、下面板為鋁合金材料，厚度0.3 mm，密度 2 700 kg/m3，彈性模量 70 GPa，泊松比0.3，屈服強(qiáng)度275 MPa；蜂窩芯子等效為二維正交各向異性材料，厚度 25 mm，密度 104.8 kg/m3，對(duì)其材料參數(shù)只定義垂直板面方向的剪切剛度為112 MPa。這種建模方式可以真實(shí)模擬艙板受力情況，使彎曲剛度由面板提供，剪切剛度由芯子提供。在起爆點(diǎn)處施加脈沖載荷，計(jì)算距其100 mm處的沖擊響應(yīng)。

2.3 有限元參數(shù)選取

為對(duì)復(fù)雜的沖擊響應(yīng)進(jìn)行預(yù)示分析，需要合理選取有限元模型參數(shù)，以保證計(jì)算的精度與效率。本文計(jì)算模型尺寸為2.2 m×3 m，有限元計(jì)算網(wǎng)格采取四邊形單元（Quad4）。密集的網(wǎng)格可以提高計(jì)算的精度以及收斂一致性，但過(guò)于密集會(huì)導(dǎo)致計(jì)算時(shí)間的延長(zhǎng)，因此需選取合適的網(wǎng)格尺寸。分別選取邊長(zhǎng) 20、15、10、5 mm 的網(wǎng)格進(jìn)行計(jì)算，艙板的位移云圖如圖3所示。可以看出，隨著網(wǎng)格尺寸的減小，云圖的邊界逐漸光滑，收斂度提高。不同網(wǎng)格尺寸對(duì)應(yīng)的沖擊響應(yīng)譜如圖4所示，當(dāng)網(wǎng)格尺寸小于10 mm后，沖擊響應(yīng)譜曲線接近重合，可以保證計(jì)算的精度。

圖3 有限元計(jì)算網(wǎng)格尺寸對(duì)艙板位移云圖的影響Fig. 3 Influence of mesh size of FEM on the displacement of the cabin board

圖4 不同網(wǎng)格尺寸下的沖擊響應(yīng)譜及局部放大圖Fig. 4 SRS for different mesh dimensions (a part is enlarged)

如前所述，顯式有限元法需要滿足穩(wěn)定性條件，即時(shí)間步長(zhǎng)必須小于臨界值。本文算例采用模型的一階固有頻率為138.44 Hz，為滿足穩(wěn)定性條件，時(shí)間步長(zhǎng)須小于2.3 ms。在滿足穩(wěn)定性條件的情況下，分別取 100、10、5、2 μs的時(shí)間步長(zhǎng)進(jìn)行計(jì)算分析，得到加速度響應(yīng)的時(shí)域圖（如圖5所示）。可以看出，4條曲線的整體走勢(shì)相同；局部放大顯示：時(shí)間步長(zhǎng)為100 μs時(shí)，由于時(shí)間步長(zhǎng)較大，致使采樣數(shù)據(jù)缺失，曲線平滑度降低；當(dāng)時(shí)間步長(zhǎng)小于10 μs后，3條曲線完全重合，此時(shí)時(shí)間步長(zhǎng)的減小并不會(huì)提高計(jì)算的精確度。因此，綜合考慮計(jì)算的精度與效率，選擇10 μs的時(shí)間步長(zhǎng)。

圖5 時(shí)間步長(zhǎng)對(duì)計(jì)算結(jié)果的影響及局部放大圖Fig. 5 Calculation result under different time steps (a part is enlarged)

3 沖擊源力函數(shù)簡(jiǎn)化模型對(duì)比分析

沖擊響應(yīng)與輸入載荷的波形、幅值以及持續(xù)時(shí)長(zhǎng)和固有周期的比值有關(guān)[2]，而固有周期一定，因此可分別從波形、幅值和持續(xù)時(shí)長(zhǎng)3個(gè)方面對(duì)不同輸入載荷進(jìn)行對(duì)比分析，并采用改進(jìn)的遞歸數(shù)字濾波法[9]進(jìn)行沖擊響應(yīng)譜計(jì)算。

3.1 波形對(duì)比

大部分沖擊載荷可由三角形脈沖、半正弦波脈沖和矩形脈沖3種基本脈沖（如圖6所示）組合而成[10]。對(duì)這3種脈沖載荷進(jìn)行對(duì)比分析，在沖擊源處施加幅值 2000 N、持續(xù)時(shí)長(zhǎng) 0.5 ms的 3種不同波形的沖擊載荷，得到距沖擊源100 mm位置處的沖擊響應(yīng)譜（如圖7所示）。

圖6 基本脈沖波形Fig. 6 Three kinds of simple pulse functions

圖7的計(jì)算結(jié)果顯示：在不同頻段，不同波形的沖擊響應(yīng)譜曲線具有不同特點(diǎn)：在低頻段，矩形脈沖的沖擊響應(yīng)譜曲線較高；在中頻段（1000～4000 Hz），半正弦波脈沖的沖擊響應(yīng)譜曲線的增長(zhǎng)減慢，出現(xiàn)極值；在高頻段，3種波形對(duì)應(yīng)的沖擊響應(yīng)譜曲線均出現(xiàn)最大值，且最大值對(duì)應(yīng)的頻率相近，其中三角形脈沖與矩形脈沖的最大值相近，大于半正弦波脈沖的。

圖7 不同波形脈沖函數(shù)對(duì)應(yīng)的沖擊響應(yīng)譜Fig. 7 SRS of three kinds of simple pulse input

3.2 幅值對(duì)比

在有限元分析模型中沖擊源處施加持續(xù)時(shí)長(zhǎng)0.5 ms，幅值分別為 1000、2000、3000、4000、5000 N的不同波形的沖擊載荷，得到距沖擊源100 mm位置處的沖擊響應(yīng)譜（如圖8所示）。各沖擊響應(yīng)譜的最大值統(tǒng)計(jì)數(shù)據(jù)見(jiàn)表1。

圖8 不同幅值脈沖函數(shù)對(duì)應(yīng)的沖擊響應(yīng)譜Fig. 8 SRS of different pulse functions of various amplitudes

表1 不同幅值脈沖函數(shù)下的沖擊響應(yīng)譜最大值Table 1 Maximum of SRS for pulse functions of different amplitudes

根據(jù)圖8與表1中的計(jì)算結(jié)果：對(duì)于波形相同、幅值不同的沖擊載荷，其沖擊響應(yīng)譜曲線的形狀相同，沖擊響應(yīng)譜最大值對(duì)應(yīng)的頻率相同，且沖擊響應(yīng)譜曲線最大值與輸入載荷幅值呈正相關(guān)；三角形與矩形脈沖最大值對(duì)應(yīng)的頻率相同，且它們的最大值相近，大于半正弦波脈沖對(duì)應(yīng)沖擊響應(yīng)譜的最大值。由此可見(jiàn)，沖擊載荷的幅值不影響沖擊響應(yīng)譜形狀，只在效果上造成沖擊響應(yīng)譜曲線的上下平移。

3.3 持續(xù)時(shí)長(zhǎng)對(duì)比

在有限元分析模型中沖擊源處施加幅值2000 N，持續(xù)時(shí)長(zhǎng)分別為 0.25、0.50、0.75、1.00、1.25 ms的不同波形沖擊載荷，得到距沖擊源100 mm位置處的沖擊響應(yīng)譜（如圖9所示）。通過(guò)觀察沖擊響應(yīng)譜曲線，總結(jié)出不同持續(xù)時(shí)長(zhǎng)脈沖下的沖擊響應(yīng)譜特性如表2所示。

圖9 不同持續(xù)時(shí)長(zhǎng)脈沖對(duì)應(yīng)的沖擊響應(yīng)譜Fig. 9 SRS of pulse functions of different durations

表2 不同持續(xù)時(shí)長(zhǎng)脈沖對(duì)沖擊響應(yīng)譜特性的影響Table 2 Effect of pulse durations on the SRS

根據(jù)上述分析，持續(xù)時(shí)長(zhǎng)主要對(duì)半正弦波脈沖沖擊響應(yīng)譜曲線的最大值、極值以及它們對(duì)應(yīng)的頻率，而對(duì)三角形脈沖的影響較小，在一定范圍內(nèi)對(duì)矩形脈沖沒(méi)有影響。

3.4 分析結(jié)果總結(jié)

總結(jié)3.1～3.3節(jié)的對(duì)比分析結(jié)果，沖擊載荷的不同變量對(duì)沖擊響應(yīng)譜的影響如表3所示。

表3 不同變量對(duì)沖擊響應(yīng)譜的影響Table 3 Influence of pulse function factors on the SRS

4 組合脈沖函數(shù)分析與試驗(yàn)數(shù)據(jù)擬合

4.1 不同組合脈沖函數(shù)計(jì)算對(duì)比分析

為分析組合函數(shù)中各基本脈沖函數(shù)對(duì)計(jì)算結(jié)果的影響，采用4組組合脈沖函數(shù)（見(jiàn)表4）進(jìn)行計(jì)算對(duì)比分析。

表4 組合脈沖函數(shù)Table 4 Combined pulse function

組合脈沖函數(shù)以半正弦波脈沖為基礎(chǔ)，疊加其他脈沖函數(shù)，對(duì)比不同波形脈沖疊加以及不同參數(shù)脈沖疊加后的沖擊響應(yīng)，得到?jīng)_擊響應(yīng)譜（如圖10所示）。

圖10 組合脈沖函數(shù)沖擊響應(yīng)譜對(duì)比分析Fig. 10 Comparison of SRS for different pulse function combinations

圖10(a)為將單一的半正弦波脈沖及其分別與三角形脈沖、矩形脈沖疊加的組合脈沖進(jìn)行對(duì)比。結(jié)果顯示：與單一半正弦波脈沖相比，當(dāng)在半正弦波脈沖上疊加三角形脈沖或矩形脈沖后，沖擊響應(yīng)譜曲線仍然存在極值，但曲線最大值增大；三角形脈沖與矩形脈沖對(duì)曲線最大值的影響相近，但是矩形脈沖會(huì)增大曲線在低頻段的值。

圖10(b)為將不同參數(shù)的半正弦波脈沖與三角形脈沖疊加后的組合函數(shù)進(jìn)行對(duì)比。結(jié)果顯示：三角形脈沖幅值的增大會(huì)使沖擊響應(yīng)譜曲線的最大值增大，半正弦波脈沖的持續(xù)時(shí)長(zhǎng)會(huì)影響曲線極值以及極值點(diǎn)的位置。

綜合兩方面的對(duì)比，脈沖參數(shù)對(duì)組合脈沖函數(shù)沖擊響應(yīng)譜的影響與其對(duì)單一脈沖函數(shù)沖擊響應(yīng)譜的影響規(guī)律相同。

4.2 試驗(yàn)數(shù)據(jù)擬合

采用3種基本脈沖函數(shù)的組合模擬沖擊源特性開(kāi)展仿真計(jì)算，并通過(guò)試驗(yàn)數(shù)據(jù)進(jìn)行反推。在進(jìn)行沖擊源函數(shù)模擬時(shí)，首先根據(jù)試驗(yàn)數(shù)據(jù)的走勢(shì)，確定組合脈沖函數(shù)的類型，再根據(jù)試驗(yàn)數(shù)據(jù)曲線的最大值、極值等特性，調(diào)整組合脈沖函數(shù)中各基本脈沖函數(shù)的參數(shù)，從而對(duì)試驗(yàn)曲線進(jìn)行擬合。

分析多個(gè)起爆點(diǎn)附近100 mm處的沖擊響應(yīng)譜，如圖11中的3條測(cè)點(diǎn)曲線，可以看到：試驗(yàn)數(shù)據(jù)沖擊響應(yīng)譜曲線初始值較低，并在低頻階段以較大斜率上升；在2000 Hz左右出現(xiàn)極值，隨后曲線斜率增大，在4000～8000 Hz范圍內(nèi)達(dá)到最大值。

圖11 模擬曲線與試驗(yàn)數(shù)據(jù)對(duì)比Fig. 11 Comparison of simulated curve and test data

根據(jù)試驗(yàn)數(shù)據(jù)曲線特性，若選取單一的三角形脈沖或矩形脈沖，則無(wú)法擬合出極值這一特性；若選取單一的半正弦波脈沖，則在滿足極值大小的情況下曲線最大值無(wú)法達(dá)到試驗(yàn)值。因此需要選取組合脈沖函數(shù)對(duì)沖擊源載荷進(jìn)行準(zhǔn)確模擬。根據(jù)上述不同脈沖函數(shù)對(duì)沖擊響應(yīng)譜的影響規(guī)律分析，首先使用半正弦波脈沖擬合試驗(yàn)曲線在中頻段的極值大小與極值點(diǎn)對(duì)應(yīng)的頻率，由試驗(yàn)數(shù)據(jù)曲線可得，沖擊響應(yīng)譜曲線的極值出現(xiàn)在1000～2000 Hz范圍內(nèi)，極值大小約1000g；經(jīng)過(guò)仿真計(jì)算，幅值為4000 N、持續(xù)時(shí)長(zhǎng)為0.6 ms的半正弦波脈沖產(chǎn)生的沖擊響應(yīng)譜曲線極值點(diǎn)對(duì)應(yīng)的頻率為1500 Hz，極值大小為1200g，與試驗(yàn)曲線極值范圍接近。但此時(shí)沖擊響應(yīng)譜的最大值與試驗(yàn)數(shù)據(jù)相差1500g左右，需要疊加三角形脈沖或矩形脈沖來(lái)增大曲線的最大值?？紤]到矩形脈沖會(huì)增大低頻段的沖擊響應(yīng)譜，而試驗(yàn)曲線在低頻段的沖擊響應(yīng)譜較小，因此只疊加三角形脈沖。經(jīng)過(guò)仿真計(jì)算，疊加幅值為1000 N、持續(xù)時(shí)長(zhǎng)為0.6 ms的三角形脈沖可以使沖擊響應(yīng)譜曲線最大值增大到4000g，得到與試驗(yàn)結(jié)果相近的沖擊響應(yīng)譜（見(jiàn)圖11）。

在時(shí)域方面，計(jì)算得到?jīng)_擊源處沖擊加速度為16 230g，距沖擊源 100 mm處沖擊響應(yīng)最大值為1204g。由試驗(yàn)實(shí)測(cè)數(shù)據(jù)可得，多個(gè)距沖擊源100 mm處的沖擊響應(yīng)最大值分布在900g～1600g，模擬值與試驗(yàn)數(shù)據(jù)相符。

由上述分析可知，組合函數(shù)可以近似地對(duì)沖擊源載荷進(jìn)行模擬，得到與試驗(yàn)實(shí)測(cè)數(shù)據(jù)相近的沖擊響應(yīng)。

5 結(jié)束語(yǔ)

本文應(yīng)用顯式有限元法對(duì)不同脈沖載荷在航天器上產(chǎn)生的沖擊響應(yīng)進(jìn)行分析和仿真計(jì)算。通過(guò)合理建立有限元模型，使用簡(jiǎn)單脈沖載荷的疊加來(lái)模擬沖擊響應(yīng)，可避免對(duì)復(fù)雜火工品裝置的建模。利用試驗(yàn)數(shù)據(jù)反推沖擊源載荷，得到與試驗(yàn)實(shí)測(cè)數(shù)據(jù)相近的沖擊響應(yīng)預(yù)示結(jié)果，驗(yàn)證了本文方法正確可行，可為相同平臺(tái)衛(wèi)星沖擊響應(yīng)預(yù)示分析的輸入載荷設(shè)定提供依據(jù)，方便后續(xù)沖擊響應(yīng)分析工作。

沖擊環(huán)境具有高復(fù)雜性與不確定性，有限元法在響應(yīng)預(yù)示方面仍有很大的局限性，因此在不斷探索與完善有限元法的同時(shí)，應(yīng)將有限元法與其他方法（經(jīng)驗(yàn)法、模擬試驗(yàn)等）相結(jié)合，更加可靠地預(yù)示沖擊響應(yīng)。