航天器火工冲击源力函数模拟分析

杨 光，刘 波，李正举，张 程

（中国空间技术研究院 通信卫星事业部，北京 100094）

0 引言

航天器所经受的冲击环境主要来源于各种火工装置的起爆，例如星箭分离、组合体分离、伸展部件的展开等[1]。其特点为持续时间短、频率范围广、加速度幅值高[2]。一般来说，火工品起爆对航天器主结构影响较小，但其引发的冲击环境对石英晶体、陶瓷器件、精密电子部件和继电器等元器件的影响很大，可能会造成这些元器件失效或损坏，甚至导致整个航天任务的失败[3]。因此有必要进行航天器冲击响应预示和试验验证工作。

冲击响应预示是航天工程中冲击环境分析与试验条件制定的基础，其难点之一在于冲击源力函数的选取。冲击载荷的高量级以及冲击波传递的复杂性，使得冲击源附近的力载荷难以准确测量[4]。因此需要建立与冲击源量级以及时间特性一致的力函数对冲击源进行模拟。本文以结构星试验测试数据反推火工冲击源力函数，拟合方法可为后续相似平台卫星冲击响应预示的输入载荷设定提供参考。

1 冲击预示有限元法

1.1 冲击预示有限元法的选择

冲击预示有限元法将实际受力结构和火工品简化为合理的有限单元动力学模型，采用数值分析手段，通过瞬态响应分析得到结构的加速度响应，基于此进一步计算冲击响应谱，使其能够合理预示航天器结构复杂的冲击环境。

有限元法可分为隐式和显式2种，二者的区别为：在时间步长方面，隐式有限元法无条件地稳定，与积分时间步长无关，而显式有限元法则存在稳定性条件[5-6]，因此隐式有限元法的时间步长可在ms级，而显式有限元法的时间步长需要达到μs级。但这并不意味着显式有限元法的计算效率更低。在显式有限元法中，系统的质量矩阵不包含耦合项，从而能够提高计算效率；且显式有限元法能够更好地处理强非线性问题。

本文将选择显式有限元法进行计算，使用的软件为MSC.Dytran。

1.2 显式有限元法控制方程

动力学系统运动微分方程可描述为

进一步改写为

其中：M为结构质量矩阵；C为结构阻尼矩阵；K为结构刚度矩阵；Fn,ext为外加载荷矢量；Fn,int为内力矢量。

由式(3)可以看出，计算加速度需要进行矩阵求逆，而在显式有限元法的计算中，系数矩阵仅为结构质量矩阵，不包含耦合项，因此无须进行矩阵分解。

采用中心差分法进行时间迭代：

中心差分法需要稳定性条件[7]

其中Tn为系统最小固有周期，即时间步长Δt必须小于某一数值，算法才能稳定。

2 冲击预示有限元简化模型

2.1 卫星太阳电池阵解锁试验

本文以某平台工程星太阳电池阵解锁过程为研究对象，使用试验实测数据反推火工冲击源模拟函数，其中试验实测数据为地面整星太阳电池阵展开时的测试数据。卫星两侧的太阳电池阵分别通过8个火工品装置固定在卫星舱板上，展开时分5次起爆，由分布在火工品压紧点附近的冲击加速度传感器测量并记录火工品起爆时产生的加速度。8个火工品压紧点分布如图1所示，其中的黑色矩形块即为压紧点。展开过程中太阳电池阵板面垂直于地面放置，并通过吊索悬挂模拟失重状态，如图2所示。

图1 卫星舱板上火工品压紧点布局Fig. 1 The layout of initiating explosive devices on satellite board

图2 太阳电池阵解锁试验状态Fig. 2 Schematic diagram of solar array deployment test

2.2 有限元模型

卫星的舱板多采用蜂窝夹层结构板，由面板、蜂窝芯子和板芯胶组成。为避免拉弯耦合效应以及固化后翘曲变形，上、下面板一般选取相同的材料和厚度[8]。根据卫星通信舱南北板材料属性，建立舱板模型，其中上、下面板为铝合金材料，厚度0.3 mm，密度 2 700 kg/m3，弹性模量 70 GPa，泊松比0.3，屈服强度275 MPa；蜂窝芯子等效为二维正交各向异性材料，厚度 25 mm，密度 104.8 kg/m3，对其材料参数只定义垂直板面方向的剪切刚度为112 MPa。这种建模方式可以真实模拟舱板受力情况，使弯曲刚度由面板提供，剪切刚度由芯子提供。在起爆点处施加脉冲载荷，计算距其100 mm处的冲击响应。

2.3 有限元参数选取

为对复杂的冲击响应进行预示分析，需要合理选取有限元模型参数，以保证计算的精度与效率。本文计算模型尺寸为2.2 m×3 m，有限元计算网格采取四边形单元（Quad4）。密集的网格可以提高计算的精度以及收敛一致性，但过于密集会导致计算时间的延长，因此需选取合适的网格尺寸。分别选取边长 20、15、10、5 mm 的网格进行计算，舱板的位移云图如图3所示。可以看出，随着网格尺寸的减小，云图的边界逐渐光滑，收敛度提高。不同网格尺寸对应的冲击响应谱如图4所示，当网格尺寸小于10 mm后，冲击响应谱曲线接近重合，可以保证计算的精度。

图3 有限元计算网格尺寸对舱板位移云图的影响Fig. 3 Influence of mesh size of FEM on the displacement of the cabin board

图4 不同网格尺寸下的冲击响应谱及局部放大图Fig. 4 SRS for different mesh dimensions (a part is enlarged)

如前所述，显式有限元法需要满足稳定性条件，即时间步长必须小于临界值。本文算例采用模型的一阶固有频率为138.44 Hz，为满足稳定性条件，时间步长须小于2.3 ms。在满足稳定性条件的情况下，分别取 100、10、5、2 μs的时间步长进行计算分析，得到加速度响应的时域图（如图5所示）。可以看出，4条曲线的整体走势相同；局部放大显示：时间步长为100 μs时，由于时间步长较大，致使采样数据缺失，曲线平滑度降低；当时间步长小于10 μs后，3条曲线完全重合，此时时间步长的减小并不会提高计算的精确度。因此，综合考虑计算的精度与效率，选择10 μs的时间步长。

图5 时间步长对计算结果的影响及局部放大图Fig. 5 Calculation result under different time steps (a part is enlarged)

3 冲击源力函数简化模型对比分析

冲击响应与输入载荷的波形、幅值以及持续时长和固有周期的比值有关[2]，而固有周期一定，因此可分别从波形、幅值和持续时长3个方面对不同输入载荷进行对比分析，并采用改进的递归数字滤波法[9]进行冲击响应谱计算。

3.1 波形对比

大部分冲击载荷可由三角形脉冲、半正弦波脉冲和矩形脉冲3种基本脉冲（如图6所示）组合而成[10]。对这3种脉冲载荷进行对比分析，在冲击源处施加幅值 2000 N、持续时长 0.5 ms的 3种不同波形的冲击载荷，得到距冲击源100 mm位置处的冲击响应谱（如图7所示）。

图6 基本脉冲波形Fig. 6 Three kinds of simple pulse functions

图7的计算结果显示：在不同频段，不同波形的冲击响应谱曲线具有不同特点：在低频段，矩形脉冲的冲击响应谱曲线较高；在中频段（1000～4000 Hz），半正弦波脉冲的冲击响应谱曲线的增长减慢，出现极值；在高频段，3种波形对应的冲击响应谱曲线均出现最大值，且最大值对应的频率相近，其中三角形脉冲与矩形脉冲的最大值相近，大于半正弦波脉冲的。

图7 不同波形脉冲函数对应的冲击响应谱Fig. 7 SRS of three kinds of simple pulse input

3.2 幅值对比

在有限元分析模型中冲击源处施加持续时长0.5 ms，幅值分别为 1000、2000、3000、4000、5000 N的不同波形的冲击载荷，得到距冲击源100 mm位置处的冲击响应谱（如图8所示）。各冲击响应谱的最大值统计数据见表1。

图8 不同幅值脉冲函数对应的冲击响应谱Fig. 8 SRS of different pulse functions of various amplitudes

表1 不同幅值脉冲函数下的冲击响应谱最大值Table 1 Maximum of SRS for pulse functions of different amplitudes

根据图8与表1中的计算结果：对于波形相同、幅值不同的冲击载荷，其冲击响应谱曲线的形状相同，冲击响应谱最大值对应的频率相同，且冲击响应谱曲线最大值与输入载荷幅值呈正相关；三角形与矩形脉冲最大值对应的频率相同，且它们的最大值相近，大于半正弦波脉冲对应冲击响应谱的最大值。由此可见，冲击载荷的幅值不影响冲击响应谱形状，只在效果上造成冲击响应谱曲线的上下平移。

3.3 持续时长对比

在有限元分析模型中冲击源处施加幅值2000 N，持续时长分别为 0.25、0.50、0.75、1.00、1.25 ms的不同波形冲击载荷，得到距冲击源100 mm位置处的冲击响应谱（如图9所示）。通过观察冲击响应谱曲线，总结出不同持续时长脉冲下的冲击响应谱特性如表2所示。

图9 不同持续时长脉冲对应的冲击响应谱Fig. 9 SRS of pulse functions of different durations

表2 不同持续时长脉冲对冲击响应谱特性的影响Table 2 Effect of pulse durations on the SRS

根据上述分析，持续时长主要对半正弦波脉冲冲击响应谱曲线的最大值、极值以及它们对应的频率，而对三角形脉冲的影响较小，在一定范围内对矩形脉冲没有影响。

3.4 分析结果总结

总结3.1～3.3节的对比分析结果，冲击载荷的不同变量对冲击响应谱的影响如表3所示。

表3 不同变量对冲击响应谱的影响Table 3 Influence of pulse function factors on the SRS

4 组合脉冲函数分析与试验数据拟合

4.1 不同组合脉冲函数计算对比分析

为分析组合函数中各基本脉冲函数对计算结果的影响，采用4组组合脉冲函数（见表4）进行计算对比分析。

表4 组合脉冲函数Table 4 Combined pulse function

组合脉冲函数以半正弦波脉冲为基础，叠加其他脉冲函数，对比不同波形脉冲叠加以及不同参数脉冲叠加后的冲击响应，得到冲击响应谱（如图10所示）。

图10 组合脉冲函数冲击响应谱对比分析Fig. 10 Comparison of SRS for different pulse function combinations

图10(a)为将单一的半正弦波脉冲及其分别与三角形脉冲、矩形脉冲叠加的组合脉冲进行对比。结果显示：与单一半正弦波脉冲相比，当在半正弦波脉冲上叠加三角形脉冲或矩形脉冲后，冲击响应谱曲线仍然存在极值，但曲线最大值增大；三角形脉冲与矩形脉冲对曲线最大值的影响相近，但是矩形脉冲会增大曲线在低频段的值。

图10(b)为将不同参数的半正弦波脉冲与三角形脉冲叠加后的组合函数进行对比。结果显示：三角形脉冲幅值的增大会使冲击响应谱曲线的最大值增大，半正弦波脉冲的持续时长会影响曲线极值以及极值点的位置。

综合两方面的对比，脉冲参数对组合脉冲函数冲击响应谱的影响与其对单一脉冲函数冲击响应谱的影响规律相同。

4.2 试验数据拟合

采用3种基本脉冲函数的组合模拟冲击源特性开展仿真计算，并通过试验数据进行反推。在进行冲击源函数模拟时，首先根据试验数据的走势，确定组合脉冲函数的类型，再根据试验数据曲线的最大值、极值等特性，调整组合脉冲函数中各基本脉冲函数的参数，从而对试验曲线进行拟合。

分析多个起爆点附近100 mm处的冲击响应谱，如图11中的3条测点曲线，可以看到：试验数据冲击响应谱曲线初始值较低，并在低频阶段以较大斜率上升；在2000 Hz左右出现极值，随后曲线斜率增大，在4000～8000 Hz范围内达到最大值。

图11 模拟曲线与试验数据对比Fig. 11 Comparison of simulated curve and test data

根据试验数据曲线特性，若选取单一的三角形脉冲或矩形脉冲，则无法拟合出极值这一特性；若选取单一的半正弦波脉冲，则在满足极值大小的情况下曲线最大值无法达到试验值。因此需要选取组合脉冲函数对冲击源载荷进行准确模拟。根据上述不同脉冲函数对冲击响应谱的影响规律分析，首先使用半正弦波脉冲拟合试验曲线在中频段的极值大小与极值点对应的频率，由试验数据曲线可得，冲击响应谱曲线的极值出现在1000～2000 Hz范围内，极值大小约1000g；经过仿真计算，幅值为4000 N、持续时长为0.6 ms的半正弦波脉冲产生的冲击响应谱曲线极值点对应的频率为1500 Hz，极值大小为1200g，与试验曲线极值范围接近。但此时冲击响应谱的最大值与试验数据相差1500g左右，需要叠加三角形脉冲或矩形脉冲来增大曲线的最大值。考虑到矩形脉冲会增大低频段的冲击响应谱，而试验曲线在低频段的冲击响应谱较小，因此只叠加三角形脉冲。经过仿真计算，叠加幅值为1000 N、持续时长为0.6 ms的三角形脉冲可以使冲击响应谱曲线最大值增大到4000g，得到与试验结果相近的冲击响应谱（见图11）。

在时域方面，计算得到冲击源处冲击加速度为16 230g，距冲击源 100 mm处冲击响应最大值为1204g。由试验实测数据可得，多个距冲击源100 mm处的冲击响应最大值分布在900g～1600g，模拟值与试验数据相符。

由上述分析可知，组合函数可以近似地对冲击源载荷进行模拟，得到与试验实测数据相近的冲击响应。

5 结束语

本文应用显式有限元法对不同脉冲载荷在航天器上产生的冲击响应进行分析和仿真计算。通过合理建立有限元模型，使用简单脉冲载荷的叠加来模拟冲击响应，可避免对复杂火工品装置的建模。利用试验数据反推冲击源载荷，得到与试验实测数据相近的冲击响应预示结果，验证了本文方法正确可行，可为相同平台卫星冲击响应预示分析的输入载荷设定提供依据，方便后续冲击响应分析工作。

冲击环境具有高复杂性与不确定性，有限元法在响应预示方面仍有很大的局限性，因此在不断探索与完善有限元法的同时，应将有限元法与其他方法（经验法、模拟试验等）相结合，更加可靠地预示冲击响应。