鸟撞飞机风挡数值模拟研究

李 旦，赵廷渝，王永虎

(中国民用航空飞行学院, 四川 广汉 618307)

鸟撞飞机风挡数值模拟研究

李 旦，赵廷渝，王永虎

(中国民用航空飞行学院, 四川 广汉 618307)

目前，鸟撞是威胁航空安全的重要因素之一，飞机风挡抗鸟撞是飞机安全飞行的重要保证。首先介绍了鸟撞飞机风挡的研究现状，利用LS-DYNA3D对鸟撞飞机风挡进行数值模拟，通过建立飞机全尺寸圆弧风挡模型及鸟体简化模型，计算得到风挡结构的变形、位移、有效应力、可能发生破坏的位置、鸟体水平与垂直方向速度、加速度等数据结果。仿真结果表明，鸟撞飞机风挡是发生在毫秒量级的非线性冲击动力学行为，整个撞击过程约5.6 ms，在T=1.8 ms时刻，风挡承受的有效应力最大，为8.304×107Pa，鸟体垂直方向加速度达到1.5228×104m/s2。同时，通过选取风挡三个网格单元，得到位移及有效应力变化历程，综合考虑位移和应变结果可知风挡正中心为受到鸟撞后最危险的部位，利用数值模拟方法进行鸟撞风挡分析，可减小时间成本，提高分析问题工作效率，为飞机风挡鸟撞适航验证提供更有效的方法。

圆弧风挡；冲击动力学；适航验证

引 言

航空运输业快速发展使得飞机数量剧增，在航线特别繁忙的区域，飞行活动愈趋频繁。由于飞机飞行速度快、鸟体体积小，机组要想在安全距离内规避与飞鸟发生撞击将会变得非常困难，使得飞鸟撞击飞机事件时有发生，轻者使飞机结构严重损坏，严重时导致机毁人亡。鸟撞飞机的严重程度主要取决于鸟体质量、撞击部位以及鸟与飞机相对撞速等因素。统计国内外鸟撞事件可知，飞机起降时发生鸟撞事件的概率最大，发生在风挡与发动机部位的撞击事件最多[1]，同时，机头、机翼前缘与起落架等部位发生的撞击事件较多，作为飞机比较薄弱和关键的部位，发生鸟撞时后果十分严重。风挡作为飞机上极其重要的部件，直接关系到飞行员的生命安全与仪表的正常使用。对于低空、高速飞行的军用飞机，保证风挡抗鸟撞性能合格相当重要，对飞行安全而言具有特殊的意义。

鸟撞风挡过程相当复杂，具有柔性撞击、大变形与高应变率等特点，属于典型的接触-碰撞弹塑性大变形动力学计算分析问题[2]，体现出几何非线性、物理非线性与接触状态非线性。由于目前的理论分析很难直接解决工程问题，因此，其分析结果不足以作为参考依据，这也使得试验研究和数值模拟成为解决这类问题的主要手段[3]。在试验方面，由于鸟撞飞机实验成本昂贵，试验周期较长，受外在条件限制，并不具有通用性。近年来，国内外学者将数值模拟作为研究鸟撞风挡问题的主要方法，例如，龚尧南[4]提出了关于鸟撞问题解藕的简化办法，避免了鸟撞问题流固耦合解法；Yao X H等[5]利用基于ALE的流固藕合方法分析了圆弧风挡受鸟撞击问题。本文利用LS-DYNA3D数值求解工具分析圆弧风挡受鸟撞击问题，提出了采用流固耦合算法分析鸟撞问题的新方法，通过仿真分析得到了一些有意义的结论。

1 数值建模与仿真

本文建立民用飞机全尺寸圆弧风挡有限元模型及鸟体的简化模型，通过关键参数选取进行数值模拟，数值求解最终得到的仿真结果与原有的实验数据结果[6]基本吻合。

1.1 建立模型

首先在Solidworks2014中建立全尺寸风挡模型以及鸟体模型。由于主要研究风挡被鸟体撞击后的动态性能和失效形式，可对风挡模型进行简化，以节约数值求解时间，避免结构复杂导致计算出错概率增多。鸟撞风挡有限元模型如图1所示，风挡模型采用壳单元，这样可以在数值模拟碰撞过程中更容易观察到风挡结构的变形、位移以及可能发生破坏位置等动态性能；鸟体采用两端带半球的圆柱体，模型采用体单元，然后导入LS-DYNA3D中进行网格划分。

图1 鸟撞风挡有限元模型

1.2 选取关键参数

风挡采用塑性随动强化材料，材料具体参数见表1。鸟体材料一般采用塑性动力学材料、状态方程材料、流体材料三类，根据CCAR91.117规定，在距机场中心7.5 km范围内、离地高度700 m以下，飞机最大速度不能超过103 m/s，因此，本次仿真中撞速设为83.3 m/s(300 km/h)，鸟体采用塑性动力学材料[6]。有关研究表明，在低速飞行条件下，鸟体组织的力学性能不均匀，当速度增加时，可以忽略不均匀性，将鸟体看作是介于流体和固体之间的一种状态，同时，当鸟撞速度大于250 km/h时，鸟体基本上被撕成碎片，表现出流体特性。常用鸟体密度一般在900～950 kg/m3之间，本次仿真鸟体材料密度设为900 kg/m3，泊松比为0.3，国军标规定，鸟撞试验中鸟体长径比应为2:1，根据鸟体质量、密度和形状，可得鸟体总长度为250 mm、直径125mm。接触方式选用*CONTACT_ NODES_ TO_ SURFACE，时间步长为0.9，终止时间为10 ms。

表1 风挡材料参数

1.3 设置边界条件

民用飞机风挡边界条件介于完全固支与完全铰支之间，白金泽等[7]研究了风挡边界条件对鸟撞风挡动响应分析结果的影响，得出风挡四周完全固支与完全铰支条件下计算所获的位移、应变曲线基本相同，并与试验结果较为相符。本文对边界条件进行了简化，采用风挡下边缘完全固支，后弧框完全铰支的边界条件，同时，在仿真后期考虑了一种后弧框为自由边界的情况。

2 数值模拟结果与分析

首先对鸟撞风挡应力图进行分析，同时，对鸟撞风挡过程中鸟体的水平速度、垂直速度、水平加速度以及垂直加速度随时间变化曲线进行分析研究，最后通过选取风挡上三个网格单元，重点分析了网格单元承受的有效应力随时间变化曲线。

2.1 应力图分析

图2所示为鸟体以速度83.3 m/s撞击飞机风挡时的应力结果图。在T=0.0018 s时刻(图2(a))，由于鸟体接触风挡瞬间速度最大，作用于风挡的冲击能量最大，使得飞机风挡在极短的时间内承受的冲击载荷达到最大值[8]，从应力图上可以发现此时刻鸟撞风挡位置集中，此处风挡承受的冲击载荷达到应力峰值，其应力峰值为8.304×107Pa。大约在T=0.0042 s时刻(图2(b))，由于风挡后弧框铰支的缘故，在边缘处出现最大应力，达到3.872×107Pa。随着碰撞过程持续进行，鸟体的冲击能量一部分被风挡所吸收，同时，由于鸟体与风挡之间存在摩擦力作用，在摩擦过程中必然会消耗鸟体一定的冲击能量，综合上述两个因素可知鸟体速度将迅速下降。在T=0.0056 s(图2(c))时刻之后，鸟体与风挡的碰撞作用结束，在风挡发生碰撞的两边位置，出现了应力最大值，但此时应力峰值相比刚接触碰撞瞬间，数值减小了两倍多，在T=0.0085 s时刻之后，在风挡后弧处出现应力最大值。

图2 鸟撞风挡等效应力图

2.2 鸟体与风挡动态特性分析

主要选取鸟撞风挡过程中的水平和垂直方向速度、加速度数据结果进行动态特性分析，对数据进行滤波处理后绘制历程图，如图3和图4所示，根据历程图可知鸟撞风挡过程中速度与加速度剧烈波动过程主要发生在前4 ms，随着碰撞时间逐渐增长，鸟体速度趋于定值。

图3 鸟体水平速度变化历程

图4 鸟体垂直速度变化历程

观察图3可知，在T=0.2 ms时刻，鸟体与风挡发生碰撞，导致水平与垂直速度快速减小，到T=0.004 s时刻，鸟体水平速度从83.3 m/s减小到50 m/s，由于风挡采用塑性随动强化材料，受到鸟体冲击后会吸收大部分冲击能量，同时，由于鸟体与风挡之间存在摩擦力作用，也会消耗一部分鸟体能量，两者碰撞结束时水平速度减小为45 m/s。由图4可知，由于鸟撞风挡发生在毫秒级时间段内，碰撞瞬间风挡承受的冲击载荷相当大[9]，使得鸟体垂直方向速度迅速增大到46.072 m/s，极可能对风挡产生破坏。鸟体撞击风挡时，其水平与垂直加速度变化历程如图5与图6所示。

图5 鸟体水平加速度变化历程

图6 鸟体垂直加速度变化历程

由图5可知，在鸟撞风挡瞬间，鸟体水平方向加速度迅速增大，T=0.001 s时刻，鸟体水平反向加速度达到最大值1.3243×104m/s2。在鸟体与风挡刚接触碰撞的瞬间，为风挡最易发生破坏作用的时间点[10]。塑性随动强化材料本身具有恢复塑性变形的能力[11]，在风挡承受变形极限后储备了大量的弹性势能，释放弹性势能的过程中造出鸟体与风挡二次冲击，但此时刻由于风挡吸能与摩擦作用，鸟体冲击能量迅速减小，因此，发生二次冲击时水平加速度峰值减小为6.7×103m/s2，垂直方向加速度峰值减小到7.8×103m/s2。在T=0.0056 s时刻之后，由于速度变化趋于定值，使得加速度数值变化为零。

文中风挡的动态特性分析，主要涉及风挡网格单元位移、有效应力变化历程[12]。如图7所示，分别选取风挡与鸟体最开始接触单元、中心单元、风挡与鸟体最后接触三个单元，将求解计算所得单元位移、有效应力数据进行滤波处理，得到单元位移、有效应力变化历程图，如图8与图9所示。

图7 风挡模型网格单元

图8 风挡网格单元位移变化历程

图9 风挡网格单元有效应力变化历程

由图8可知，鸟体在风挡上滑移时依次接触S581、S721、S881三个网格单元，由于碰撞过程中鸟体速度持续减小，造成三条曲线产生波峰时出现延迟。鸟体首先接触碰撞S581号网格并产生曲线波峰，再滑移到S721号网格，导致该网格位移变形量达到17.57 mm，鸟体滑移到S881网格最终离开风挡。鸟撞风挡过程也是塑性随动强化材料持续储能的过程，一旦风挡弹性势能大于鸟体本身的冲击能量，风挡必然会通过恢复塑性变形释放能量，造出网格产生正向位移变形[13]，最大变形量达到5.5374 mm。观察风挡模型可知，风挡中心距离风挡边界位置最远，随着作用于风挡模型的冲击应力波面积逐渐扩大，当冲击应力达到风挡中心位置时产生的剪切作用最强，造成S721号网格位移变形量最大。对比三条曲线位移波峰可知，风挡最大位移变形达到10-2量级[14]，与风挡厚度为同一量级。

由图9可知，由于风挡S581网格是最开始接触点，在两者接触瞬间鸟体的冲击能量最大，造成该风挡网格位置处承受的应力最大，最大值为7.5273×107Pa。综合风挡网格位移变形和有效应力动态特性可知，风挡正中心是鸟撞后极可能发生破坏的部位，因此，在鸟撞验证试验中应当着重考虑[15]。

3 结 论

本文基于飞机风挡实体模型，运用LS-DYNA建立能够有效模拟鸟撞风挡有限元模型，通过分析鸟撞风挡冲击动态响应，最终得到以下几点结论：

(1) 从等效应力图可以看出，在T=0.0018 s时刻，飞机风挡所承受的冲击载荷达到最大值。

(2) 根据网格位移变化历程图可得，接近中心位置附近位移变形量最大，达到17.57 mm。

(3) 综合考虑风挡网格位移和应力应变可知，风挡中心为鸟撞后最危险的位置，同时，风挡后圆弧处出现较大应力，造成后弧框附近也可能发生破坏，这在鸟撞风挡试验中应当重点考虑。

近年来，鸟撞风挡问题成为航空领域研究的新热点，由于传统有限元方法对鸟撞问题进行分析与研究时，鸟体网格单元经过碰撞后开始扭曲，模态失真导致精度下降[16]，当时间步长缩短时，很可能导致求解过程终止，且数值仿真中存在很多人为假设，不能真实地反映鸟撞问题的本质，因此，采取数值模拟与缩比模型试验相结合的方法，由仿真结果为试验提供支持，一定程度上减少不必要的破坏性试验，提高审定效率，将是未来研究的发展方向。

[1] 宋春艳,朱广荣.大型民机复杂结构抗鸟撞动力学分析与适航验证技术研究[C]//中国振动工程学会.第十届全国振动理论及应用学术会议论文集,2011:1101-1111.

[2] 刘兵,王红红,邓勇军.基于ANSYS/LS-DYNA分析平头弹侵彻间隙式双层靶的时效模型[J].四川理工学院学报:自然科学版,2015,28(3):11-15.

[3] 王礼立,朱锡雄,施绍裘,等.鸟撞高速飞机风挡若干问题的冲击动力学研究[J].航空学报,1991,12(2):27-33.

[4] 龚尧南,许素强.飞机风挡透明件的鸟撞分析[J].航空学报,1991,12(2):73-78.

[5] YAO X H,ZHAO L M,LIU X M.Numerical simulation for arc windshields of aircrafts subjected to bird impact[J].上海交通大学学报,2004,38(zl):147-151.

[6] 姚小虎,韩强,张晓晴,等.飞机圆弧风挡抗鸟撞试验研究[J].振动与冲击,2005,25(5):417-422.

[7] 白金泽,孙秦.飞机风挡结构抗鸟撞一体化设计技术研究[J].力学与实践,2005,27(1):14-18.

[8] MEGUID S A,MAO R H,NG T Y.FE analysis of geometry effects of an artificial bird striking an aeroengine fan blade[J].International Journal of Impact Engineering,2008,35(6):487-498.

[9] 王文智,万小朋,郭葳.民机风挡结构抗鸟撞仿真分析与设计[J].西北工业大学学报,2009,27(4):481-485.

[10] 朱书华.鸟撞飞机风挡动响应分析与仿真试验平台研究[D].南京:南京航空航天大学,2009.

[11] 贾建东.飞机典型结构抗鸟撞设计与分析[D].南京:南京航空航天大学,2010.

[12] 张大海,费庆国,刘宏月.基于真实形状鸟体的撞击方向对鸟撞分析影响研究[J].振动与冲击,2015,34(22):103-108.

[13] 谢灿军,童明波,刘富,等.民用飞机平尾前缘鸟撞数值分析及试验验证[J].振动与冲击,2015,34(14):172-178.

[14] LIU L Z,MA Z Z,PANG S H,et al.Investigation of the bird shape effects on the simulation results of bird striking laminated aircraft windshields[J].Advanced Materials Research,2012,418-420(6):72-76.

[15] 李成,郭文辉,铁瑛.采用流固耦合方法的复合材料层合板鸟撞分析[J].郑州大学学报:工学版,2015,36(6):76-79.

[16] 刘永强,冯震宙,唐长红,等.不同铺层角度GLARE层板复合材料的抗鸟撞性能[J].机械工程材料,2016,40(5):105-110.

Numerical Simulation of Bird Impact on Aircraft Windshield

LIDan,ZHAOTingyu,WANGYonghu

(Civil Aviation Flight University of China, Guanghan 618307, China)

At present, the bird strikes is one of the important factors that threaten aviation safety, and aircraft windshield anti-bird strike is an important guarantee for aircraft safety flight. Firstly, the research status of the bird's windshield is introduced. The numerical simulation of the bird's windshield is carried out by using LS-DYNA3D. The deformation, displacement and effective stress of the windshield structure, the possible damage to the location, bird body level, vertical direction speed, acceleration and other data results are calculated by establishing the full-size circular windshield model and the bird body simplified model. The simulation results show that the wind impact of the aircraft is in the order of millisecond nonlinear dynamic behavior. The whole impact process is about 5.6 ms; whileT=1.8 ms, the effective stress on the windshield, which reaches 8.304×107Pa, is maximum , and the vertical acceleration of birds is 1.5228×104m/s2. At the same time, by selecting the three grid elements of the windshield, the displacement and the effective stress change process are obtained. The displacement and strain results show that the center of the windshield is the most dangerous part of the bird; and analyzing the windshield by the numerical simulation method can reduce the cost of time, improve the efficiency of the problem, and provide a more effective method for airworthiness.

bird impact; arc windshield; impact dynamics; airworthiness verification

2017-03-14

国家自然科学基金民航联合研究基金项目(U1333133);四川省教育厅科研项目(16ZB030)

李 旦(1992-),男,甘肃天水人,硕士生,主要从事民机水上迫降数值方法方面的研究,(E-mail)1171085449@qq.com; 赵廷渝(1965-),男,教授,硕士,主要从事民用航空发动机性能及控制技术方面的研究,(E-mail)2494238988@qq.com

1673-1549(2017)04-0036-05

10.11863/j.suse.2017.04.07

O347

A