A320撞击刚性靶体的数值模拟及冲击载荷工程模型验证

温丽晶， 张春明， 郭超， 高宇， 段璞， 欧卓成， 段卓平

(1.环境保护部核与辐射安全中心， 北京 100082 2.北京理工大学 爆炸科学与技术国家重点实验室， 北京 100081)

为防止大型商用飞机恶意撞击核电站的恐怖袭击事件造成安全壳损坏和放射性物质泄漏，给环境和人员带来危害，美国核管理委员会2009年修订了联邦法规10 CFR 50[1]，要求新建核电站申请时必须评估大型商用飞机撞击对核电厂的影响，中国国家核安全局2016年修订HAF102法规时在此方面也提出了要求[2]. 美国核管理委员会还在2011年发布的飞机撞击超设计基准事件评估相关的导则RG 1.127[3]中认可了美国电力研究院提出的新建核电厂设计的飞机撞击评估方法(NEI 07-13)[4]，该方法提出靶体的整体结构破坏可以使用载荷时间历史分析法来完成，不需要建立复杂的飞机模型，计算简便，受到研究者广泛使用，应用该方法首先要确定飞机的冲击载荷时间曲线.

1968年，Riera[5]假设靶体为刚性体首次提出了软冲击模型. 之后，多位学者对该模型进行了验证和修正[6-7]，在Riera模型冲击载荷计算式第二项中引入了修正系数α. 1993年，Sugano等[8]用F-4幻影战斗机完成了世界上第一次全尺度飞机碰撞试验，认为修正的Riera模型中系数α取0.9时与实验值最为接近；2017年，Wen等[9]参考国产C919设计飞机模型完成了撞击试验，提出了冲击载荷工程计算模型，确定了适用于试验的两种飞机模型的修正系数α值分别为1.0和0.8. 上述试验研究大都局限于一定撞击速度，对不同撞击速度下系数α值如何选取未见报道，同时，受到实验条件的限制，目前未开展过大型商用飞机撞击试验，系数α取F-4撞击试验确定的值[4]，由于大型商用飞机结构与F-4差异较大，其合理性值得商榷.

文中建立了空客A320撞击刚性靶体的有限元模型，计算飞机在不同撞击速度下的冲击载荷时间曲线，验证该模型对计算大型商用飞机恶意撞击的冲击载荷的适用性，并确定模型中系数α和γ与撞击速度的对应关系.

1 飞机及靶体模型建立

1.1 有限元模型

飞机结构复杂，包含的零部件繁多，建立完整的飞机模型非常困难，需要对飞机结构进行简化. 简化原则是保证飞机模型沿中轴线的总的外形和结构、线密度、质量分布以及强度与真实飞机的空间分布一致，重点考虑飞机的主梁、主框架、机翼横梁等主要承力框架结构和蒙皮结构. 各种复杂形状的梁简化为几何形状相对简单的工字梁、U型梁，不同厚度的飞机蒙皮取统一的厚度，忽略机舱内饰、电子设备、座椅、燃油、货物和附属设施等部分的几何特性，简化为质量点的形式分布在对应区域.

采用三维建模软件CATIA，依据空客A320的SRM手册[10]和AMM手册[11]对其主要结构进行建模. 使用HyperMesh几何型面网格自动划分模块对飞机几何模型进行网格划分，网格的尺寸为50 mm，网格的类型为混合型. 飞机的蒙皮和机身框等主要的支撑结构均为薄壁结构，采用壳单元(SHELL 163)来离散，飞机的长桁采用梁单元(BEAM188)进行离散，燃油和内饰等结构采用质量点(MASS166)的形式分布在相应区域. A320飞机的有限元模型如图1所示. A320飞机有限元模型的单元总数为120.36万，其中壳单元数量为69.53万，梁单元数量为6.35万，质量点数量为44.49万.

图1 A320飞机的有限元模型Fig.1 A320’s finite element model

撞击目标靶体尺寸为60 m×24.6 m，厚度为4.5 m. 目标靶体采用实体单元(Solid 164)进行离散，单元尺寸为300 mm，单元总数为24.6万. 飞机碰撞过程数值模拟采用典型的有限元显式算法，该数值模拟采用的计算方法、材料动力学模型等已通过飞机模型撞击混凝土靶体的数值模拟计算验证[12].

1.2 材料模型和参数

由于文中重点考虑飞机对刚性靶体的冲击载荷特性，不关注飞机撞击破坏过程的细节,考虑到整个计算规模，文中将撞击目标靶体设置为刚体，对于飞机使用的金属材料，忽略其温度软化效应，选用简化的Johnson_Cook模型，即Simplified_Johnson_Cook模型. 该模型形式如下：

(1)

表1 材料模型参数[13]Tab.1 Parameters in the material model[13]

2 数值计算结果分析

模拟了A320飞机分别以150，200，250和300 m/s 4种不同撞击速度V0垂直与刚性靶体的碰撞过程. 图2为A320飞机以200 m/s速度垂直撞击刚性靶体的典型图像.计算结果显示，飞机碰撞过程是一个渐进式压溃的过程，机体结构发生高度非线性的变形、断裂和飞散等撞毁现象. 图3显示了中央翼盒的碰撞破坏过程. 取不同时刻靶体的撞击面上所有单元力的总和可以得到撞击过程中靶体受到的冲击载荷曲线. 4种撞击速度下计算得到的冲击载荷

图2 飞机以200 m/s速度撞击刚性靶体的典型图像Fig.2 Typical images of an A320 airplane impacting upon a rigid target under the initial impact velocity of 200 m/s

曲线如图4所示，可以看出，飞机撞击速度V0越大，冲击载荷的峰值越大，峰值出现的时间也越早.

图3 中央翼盒的碰撞过程Fig.3 Impact process of the central wing box

图4 不同撞击速度下的冲击载荷数值计算结果对比Fig.4 Numerical results for the impact load time curve under different initial impact velocities

3 冲击载荷工程模型验证

在前期的研究中，通过分析飞机模型撞击试验结果，并与Riera理论模型计算的冲击载荷进行对比，提出了冲击载荷的工程计算方法，可以在不知道飞机压损载荷分布的情况下计算飞机的冲击载荷，计算表达式为[9]

γ=Pc(t)/F(t)=IPc(t)/IF(t),

(2)

式中:F(t)为冲击载荷;α为考虑了飞机压碎部分结构破坏特性和压碎质量在靶体表面飞散分布而引入的修正系数;γ为冲击载荷载荷系数，为压损载荷Pc(t)占冲击载荷F(t)的比例;μ为沿飞机轴线分布的线质量密度;x(t)为自飞机头部算起的飞机破坏长度;V为飞机撞击破坏过程的碰撞速度;Pc(t)为压损载荷;IPc(t)为压损载荷的冲量;IF(t)为冲击载荷的冲量.

3.1 冲击载荷工程模型相关系数α和γ的确定

从A320飞机撞击刚性靶体的数值模拟计算结果中，可以提取出飞机撞击过程中未变形部分的加速度变化曲线. 根据Riera理论模型的假设[5]，当已知飞机的质量、线密度及初始速度，由飞机未变形部分的加速度，通过迭代方法即可求出在不同时刻飞机的压损载荷Pc(t)[9]. 图5为典型撞击速度下计算得到的压损载荷Pc(t)与数值模拟计算得到的冲击载荷Fa(t)的对比，可以看出压损载荷的峰值出现在发动机与靶体撞击的阶段，这是由于飞机机体结构的强度较低，而发动机的强度较大，与实际情况符合. 对压损载荷时间曲线Pc(t)和冲击载荷曲线Fa(t)进行积分，可以得到压损载荷的冲量时间曲线IPc(t)及冲击载荷的冲量时间曲线IFa(t)，进而可得到撞击过程压损载荷冲量IPc(t)与冲击载荷冲量IFa(t)之比γ(t)，如图6所示. 可以看出，在整个撞击过程中γ(t)值除撞击的初始时刻有一个短时间的尖峰外，其它时间都处于一个稳定值0.17，因此，200 m/s撞击速度下，冲击载荷工程模型中γ可以取常数0.17.

图5 200 m/s撞击速度下飞机的压损载荷和冲击载荷对比Fig.5 Comparison of the crushing load of the airplane with the impact load under the initial impact velocity of 200 m/s

图6 200 m/s撞击速度下γ随时间变化曲线Fig.6 Parameter γ under the initial impact velocity of 200 m/s

对于飞机以200 m/s速度撞击的典型工况，采用修正的Riera模型计算了α在取不同值时的冲击载荷时间曲线和冲击载荷的冲量时间曲线，分别与数值模拟计算结果对比，如图7所示，可以看出：当α取1.05时，使用Riera模型计算的结果与数值模拟计算的结果吻合较好.

图7 α取不同值时的冲击载荷及冲击载荷冲量曲线与数值计算结果对比Fig.7 Comparison of the two numerical results of both the impact loads and the impulses

同理，根据上述方法可以确定飞机以150，250和300 m/s速度撞击时冲击载荷工程模型中系数α和γ的合理取值. 表2列出了4种撞击速度下冲击载荷工程模型中系数α和γ的取值.

表2 不同撞击速度下冲击载荷工程模型中α和γ的取值

Tab.2 Values ofαandγin the impact load engineering model under different initial impact velocities

序号撞击速度V0/(m·s-1)γα11500.1891.0022000.1701.0532500.1531.1043000.1421.20

图8为系数α和γ与飞机撞击速度V0的关系图，可知，当飞机撞击速度在150～300 m/s范围内，工程公式中的系数α、γ与飞机撞击速度V0的关系表示式为

(3)

γ=0.127+7.0×10-7(V0-449)2.

(4)

图8 α、γ与飞机撞击速度V0的关系Fig.8 Relationship of α and γ with V0

3.2 工程计算结果与数值计算结果对比

使用工程公式迭代计算飞机以不同速度垂直碰撞刚性靶体的冲击载荷及其冲量，分别与数值模拟计算结果进行对比，如图9所示，两种方法获得的冲击载荷时间曲线差异较小，曲线特征基本一致. 对比具体计算数据表明：冲击载荷峰值误差较小，撞击速度为200 m/s时，两个峰值的平均误差最大，为10.01%；总冲量误差很小，撞击速度为150 m/s时，总冲量误差最大，为1.69%. 因此，α、γ与飞机撞击速度的关系表示式是合理的，在不知道飞机压损载荷的情况下，可以使用工程公式近似计算飞机的冲击载荷时间曲线.

图9 200 m/s撞击速度下工程公式和数值模拟计算结果对比Fig.9 Comparison of the calculated result by the engineering model with that calculated by the numerical simulation under the impact velocity of 200 m/s

4 结 论

文中建立了详细的空客A320有限元模型，模拟了不同撞击速度下飞机与刚性靶体的碰撞过程，结果表明：

① 飞机碰撞过程是一个即撞即坏的渐进式压溃过程，飞机撞击速度越大，靶体上承受的冲击载荷(力)的峰值越大，峰值出现的时间也越早.

② 得到了Riera模型修正系数α及工程模型冲击载荷系数γ与飞机撞击速度V0之间的对应关系；不同撞击速度下数值模拟与工程模型计算的冲击载荷曲线特征和幅值吻合较好，验证了冲击载荷工程模型的合理性.

③ 利用冲击载荷工程模型可以在飞机压损载荷分布未知的条件下，计算得到任意速度下的飞机冲击载荷曲线，避免了利用Riera模型计算冲击载荷曲线时需提前确定飞机压损载荷分布的难题，为确定大型商用飞机撞击核电厂整体结构安全评估用冲击载荷时间曲线提供了更为便捷的计算方法.