郭 鹏， 刘志远， 张桂昌， Reza Hedayati， 张俊红

(1. 天津理工大学 天津市先进机电系统设计与智能控制重点实验室，天津 300384 2.天津大学 内燃机燃烧学国家重点实验室，天津 300350；3. 中国民航大学 航空工程学院，天津 300300；4. 代尔夫特理工大学 航空航天工程学院，荷兰 代尔夫特 2600AA)

鸟撞对飞机运行安全有重要影响。根据美国联邦航空管理局(federal aviation administration, FAA)统计[1]，从1988年到2018年间，野生动物与飞机相撞已造成超过263架飞机损毁，287人丧生。图1(a)为从1990-2017年间野生动物撞击飞机的事故统计，2017年较前年增加1 069起，其中的鸟撞所占的比例达到95%。图1(b)为飞机各部件发生撞鸟的比例[1]，其中发动机比例最高达到37.21%。鸟撞时，风扇最先受到撞击，其抗鸟撞能力成为影响飞行安全的重要因素。研究鸟撞风扇过程对于发展风扇安全强度设计理论具有重要的科学意义，对于风扇叶片适航认证具有主要的工程意义。

图1 野生动物撞击飞机次数统计及各部件发生撞鸟比例

国内外在鸟撞模拟中常使用圆柱体、半球端圆柱体、椭圆体和球体来替代真实的鸟体，该方法可以反映鸟躯干的主要质量和形状对碰撞过程的影响[2]，但忽略了其他部位对碰撞过程的影响。近年来，国内外开始使用更加复杂的鸟模型进行研究。Hedayati等[3]建立了真实绿头鸭的SPH模型，研究分析鸟的几何形状和撞击方向对鸟撞平板过程的影响，该模型相比于传统的替代鸟模型更接近于试验结果。此后Hedayati等[4]使用多个简单几何体对白枕鹊鸭的形状进行简化，建立了一种在几何形状类似于白枕鹊鸭的新的替代鸟模型，并且发现半球端圆柱体模型更接近真实鸟模型以尾部方向撞击平板的计算结果，椭圆体模型可以作为腹部撞击方向的替代模型。McCallum等[5-6]利用SPH方法，基于生物识别和已发布的CT扫描数据建立了一种多材料鸟模型，该模型的数值模拟结果与经典流体动力学理论非常吻合。Zhang等[7-8]建立带有更多几何特征的帆背潜鸭模型，该模型冲击刚性靶板的计算结果与试验吻合较好。并使用帆背潜鸭模型撞击静止风扇，其研究表明鸟的几何以及撞击方向对鸟撞风扇有着重要且不可忽略的影响。寇建峰等[9]针对雷达罩以及机翼前缘开展不同姿态鸟体的鸟撞研究，发现鸟体姿态对结构撞击响应有很大的影响，并且提出结构的抗鸟撞安全评估必须考虑不同鸟体姿态的影响。杨杰[10]模拟鸟撞整级旋转风扇的不同位置，结果表明撞击位置越高造成的损伤越大，损伤叶片损失的能量越大。

鸟撞叶片过程是一个高度非线性的瞬态冲击过程，上述研究表明撞击姿态以及撞击位置对叶片损伤有着不可忽略的影响，同时叶片的扭转结构使得鸟撞风扇叶片过程更加复杂。现有研究中，对真实鸟与真实航空发动机风扇叶片碰撞过程的研究较少，且简化鸟体撞击叶片过程中仅考虑简化鸟体的吸入速度，未考虑风扇旋转速度造成鸟体-叶片相对运动速度的影响。本文基于真实绿头鸭SHP模型，模拟鸟撞击固定转速下的某航空涡扇发动机整级风扇过程，分析不同撞击位置以及不同撞击姿态对风扇的动力学响应及损伤变形情况的影响。

1 数值模型

1.1 网格模型建立

FAA统计数据表明[1]，在1990—2017年绿头鸭撞击飞机的总次数达到1 003次，为与民航发动机发生碰撞次数最频繁的鸟类之一。因而本文将绿头鸭作为鸟撞模拟的原型，鸭重820 g，体长40 cm，翼展45 cm。在医学成像中心对绿头鸭进行CT扫描，如图2(a)，获取1 566张医学数字成像(DICOM)图像，如图2(b)。利用DICOM图像及SPH方法对绿头鸭进行建模，建模过程中考虑绿头鸭的内部空腔(空腔部位没有SPH单元)，最终模型由41 685个SPH单元组成，每个单元的质量为0.019 1 g，如图2(c)。

图2 绿头鸭建模过程

国内外鸟撞数值模拟中替代鸟体模型目前普遍采用半球端圆柱体模型，但是半球端圆柱体模型只是对鸟体躯干的简化，而没有将鸟体其他部位考虑进来，不是对鸟体真实外形的建模。为了验证真实鸟模型相比于传统替代鸟体模型具有更准确的数值模拟结果，本文选取传统替代鸟体模型(半球端圆柱体模型)作为真实绿头鸭模型的对比。半球端圆柱体的长径比选取为2[11]，根据鸟的质量0.8 kg以及鸟的密度938 kg/m3可以得到半球端圆柱体的直径D=0.086 9 m和长度L=0.173 8 m。利用SPH法对半球端圆柱体离散化，粒子距离为3.2 mm，模型由32 038个SPH粒子组成每个粒子质量为0.024 97 g,如图3。

图3 传统替代鸟体SPH模型

本文研究的发动机风扇由24块窄弦叶片组成，叶身高度为603.2 mm,初始扭转角为61.3°。用3 mm六面体单元对整级风扇进行有限元网格划分，获得1 520 208网格单元，如图4所示。

图4 涡扇发动机风扇有限元网格模型

1.2 材料模型

鸟材料采用Null材料模型结合Gruneisen状态方程来描述；在空材料参数中提供模型的本构关系计算黏性应力，使用状态方程来计算压力，其中材料密度设为为938 kg/m3。Gruneisen状态方程按下式定义压缩材料的压力：

(1)

对于膨胀材料，Gruneisen状态方程按下式定义压力：

p=ρ0C2μ+(γ0+aμ)E

(2)

在上述Gruneisen状态方程中C是vs(vp)曲线的截距；S1，S2，和S3是vs(vp)曲线的斜率系数[12]；γ0是Gruneisen常数，a是的一阶体积修正量。

(3)

式中：C=1 480，S1=1.92，S2=0，S3=0。

风扇叶片采用的材料是Ti-6Al-4V钛合金材料，采用塑性动力学模型，材料参数见表1。

表1 叶片材料参数

1.3 边界条件

沿叶根向上选取1/6、2/6、3/6、4/6、5/6叶高五个撞击高度，研究撞击位置对鸟撞风扇叶片过程的影响，如图5所示。

图5 鸟体撞击位置

分别选取15种鸟体吸入姿态研究鸟体撞击姿态对撞击过程的影响。如图6所示，自然飞行头向风扇为0°姿态，鸟体绕Y轴依次旋转45°得到Y-45°、Y-90°、Y-135°、Y-180°、Y-225°、Y-270°、Y-315°姿态；鸟体绕Z轴依次旋转45°得到Z-45°、Z-90°、Z-135°、Z-180°、Z-225°、Z-270°、Z-315°姿态。

图6 撞击姿态示意图

模拟中，鸟的飞行速度为116 m/s，风扇转速为3 800 r/min。以往鸟撞风扇研究大多未考虑风扇的转动，鸟沿风扇轴线撞向叶片，鸟撞过程中只与一或两个叶片接触，发生正向碰撞。实际上，鸟撞向旋转的风扇，会与多个叶片发生接触，并且风扇叶片和鸟的相互作用包括撞击和对鸟体的切割两部分。为更加真实的模拟鸟撞风扇的情况，鸟的速度设为该位置处相对于旋转风扇的相对速度[13]，如图7。相对速度的大小计算公式为：

图7 鸟相对叶片速度

(4)

v叶片=rω

(5)

式中：ω为风扇的转速；r为撞击位置旋转半径；v相对为风扇叶片撞击位置处的切向速度。

2 模型验证

本文通过与Wilbeck试验数据[14]进行对比来进行模型准确性验证。Wilbeck采用压缩空气炮装置将质量为1 kg的鸟体分别以116 m/s、225 m/s和253 m/s的速度射出，与刚性目标板撞击，在刚性平板中心位置安装压力传感器以获取撞击中心位置的压力时间曲线，实验简图，见图8。

1.气罐及空气炮； 2.鸟弹及弹底板； 3.炮管； 4.压力传感器； 5.刚性目标板； 6.弹底板止动装置。

进行半球端圆柱体鸟模型与绿头鸭真实鸟模型撞击平板的模拟，分别获取116 m/s、225 m/s、253 m/s的不同速度撞向固定平板时平板中心的压力曲线，并与对应速度下Wilbeck的试验数据进行比较。

本文建立的平板模型为直径0.6 m厚度0.06 m的圆盘，采用六面体实体网格对圆盘进行离散化，获得18 000网格单元。该数值模拟中，绿头鸭以尾部垂直撞向平板。通过固连接触将壳单元节点约束、限定在圆盘外表面节点上，以获取撞击中心的压力曲线。平板模型和壳单元的布置位置及其相对于鸟模型的位置，如图8。

图9 鸟体撞击平板

Heimbs研究结果表明，鸟撞刚性目标板撞击中心压力响应分为四个阶段。第一阶段为初始撞击阶段，初始撞击阶段的峰值压力称为雨贡纽压力，第二阶段为压力衰减阶段，第三阶段为恒定流动阶段，第四阶段为压力结束阶段。不同撞击速度下，碰撞中心压力曲线仿真与试验对比如图10所示，可以清晰的看出相比于半球端圆柱体模型，绿头鸭模型各阶段压力响应与Wilbeck实验结果吻合更好。尤其是在初始撞击阶段吻合更好，对比116 m/s、225 m/s以及253 m/s撞击速度下的撞击中心峰值压力，绿头鸭鸟体模型的峰值压力都要比半球端圆柱体模型更加接近Wilbeck实验的峰值压力。造成这个结果的主要原因是绿头鸭模型与实际鸟外形相似，并且同时考虑到绿头鸭的内部空腔，对空腔进行无粒子化处理，使得模型更加接近真实的鸟，而半球端圆柱体只是对鸟躯干的一个近似的模拟，并且半球端圆柱体模型内部并没有空腔，因而绿头鸭SPH模型结果相比于传统的半球端圆柱体模型更加接近Wlibeck的实验结果，真实绿头鸭SPH模型可更真实反应鸟撞动力学特性，可进行下一步研究。

图10 撞击中心处的压力曲线

3 结果与讨论

参考Zhang等分析鸟撞风扇过程中风扇损伤的判据，为了分析撞击位置以及撞击姿态对叶片不同位置损伤的影响，本文也将讨论分析叶根及前缘的等效应力、撞击力及前缘最大位移。叶根等效应力用于判断叶根处得损伤状况，前缘等效应力以及前缘最大位移用于判断前缘损伤状态，前缘是撞击过程中出现大变形的位置因而增加前缘最大位移来表示前缘的变形情况。

3.1 撞击位置对风扇叶片损伤的影响

图11为撞击过程中不同时刻鸟体状态，撞击过程中鸟体被叶片切割成鸟块，随后鸟块与叶片发生碰撞，鸟块在与叶片高速的碰撞过程中发生流变。

图11 鸟撞风扇过程中不同时刻鸟的状态

图12(a)为鸟与风扇的接触力时间曲线，每条接触力-时间曲线有多个波峰且每个波峰的高度不相同，接触力-时间曲线波峰数量与鸟被切割成鸟块的数量相同，且每条接触力-时间曲线每个波峰的大小与叶片撞击的鸟块大小相关。撞击位置为4/6、5/6叶高时整体接触力要小于撞击位置为1/6、2/6、3/6叶高的接触力，如表2所示，1/6、2/6、3/6、4/6、5/6位置接触力峰值为119 kN、155 kN、109 kN、67.2 kN、35.2 kN。从叶根位置至2/6叶高位置叶片的扭转幅度很小，接触力主要受撞击速度影响，撞击高度越大，撞击位置叶片线速度越大，相对撞击速度越大，撞击接触力越大。从2/6叶高位置至叶尖位置，叶片扭转幅度增大，接触力主要受叶片扭转角度影响，撞击高度越大,叶片扭转角越大，叶片对鸟体的切割作用更加明显，叶片法向的拍击减弱使得接触力逐渐减小。

表2 不同撞击位置的最大接触力、叶根等效应力和前缘等效应力

文献[15]计算了圆柱鸟模型以102 m/s速度撞击转速3 500 r/min风扇3/6叶高过程中叶片瞬态冲击响应特性，结果表明，叶片峰值应力1.396 GPa, 在本文研究中鸟体以116 m/s速度撞击转速3 800 r/min风扇3/6叶高过程中叶片峰值应力为1.53 GPa。由于叶片模型、鸟体模型及计算工况的差异，本研究叶片峰值应力计算结果与文献[15]计算结果略有差异但处于同一量级且数值接近。因此可以认为本研究中等效应力计算结果合理可信，可用于进一步讨论分析。

图12(b)为不同撞击位置下叶根应力-时间曲线，1/6、2/6、3/6、4/6、5/6叶高位置叶根最大等效应力为1.23 GPa、1.24 GPa、1.10 GPa、0.969 GPa、与0.479 GPa。图12(c)为不同撞击位置最大鸟块撞击叶片的前缘等效应力-时间曲线，1/6、2/6、3/6、4/6、5/6叶高位置前缘最大等效应力为1.18 GPa、1.62 GPa、1.53 GPa、1.03 GPa与0.479 GPa。不同撞击位置下叶根及叶片前缘处最大等效应力变化规律与接触力变化规律相似，其产生原因也与接触力变化原因相同。

图12 接触力、叶根应力及前缘应力

除5/6叶高位置外，其他撞击位置叶根及叶片前缘处最大等效应力都大于材料的屈服极限0.95 GPa，会在叶根及叶片前缘处产生塑性变形造成不可逆损伤，撞击位置为2/6叶高时，叶片叶根、前缘处损伤达最大。

3.2 撞击姿态对风扇叶片损伤的影响

图13为鸟以Y-315°姿态撞击叶片过程中不同时刻的应力分布，叶片最大应力主要出现在前后缘、前后叶根位置。相应区域在撞击过程中最可能发生塑性变形、撕裂等损伤[16]。

图13 不同时刻叶片应力分布图

图14(a)，为绕Y轴旋转得到的7个姿态的鸟撞击风扇过程中的接触力时间曲线。绕Y-0°姿态的接触力峰值为105 kN，在所有姿态中最小。Y-270°和Y-315°姿态的接触力峰值要选大于其他姿态，接触力峰值都是229 kN。图14(b)，为绕Z轴旋转得到的7个姿态的鸟撞击风扇过程中的接触力时间曲线，鸟体Z轴姿态对接触力峰值影响较小，Z-45°、90°、135°、180°、225°、270°、315°姿态的接触力峰值接近，分别为153 kN、138 kN、165 kN、145 kN、144 kN、158 kN、145 kN。对比15个姿态的接触力峰值，Y-270°和Y-315°姿态的接触力最大。

图14 不同姿态接触力

图15为各撞击姿态下，与最大鸟块相撞的叶片前叶根的等效应力，各撞击姿态下前叶根最大等效应力如表3所示，Y-135°、Y-270°、Y-315°、Z-135°和Z-315°撞击姿态下，撞击过程中前叶根最大等效应力分别为1.99 GPa、1.98 GPa、2.00 GPa、1.98 GPa和2.04 GPa，高于其他姿态。图16为Y轴8种姿态和Z轴8种姿态下，与最大鸟块相撞的叶片后叶根的等效应力时间曲线，各撞击姿态下后叶根在撞击过程中最大等效应力如表3所示，Z-135°姿态在撞击过程中后叶根受到的等效应力最大，为1.94 GPa。在撞击过程中前叶根受到的应力普遍大于后叶根，因而在撞击过程中前叶根受到的损伤大于后叶根。

表3 不同姿态最大接触力、前后叶根等效应力与前缘接触处最大位移

图15 不同姿态前叶根等效应力

图16 不同姿态后叶根等效应力

选用与最大鸟块相撞的叶片前缘接触处节点最大位移表征鸟撞风扇过程中与最大鸟块撞击的叶片前缘的变形情况。通过比较发现，Y-270°姿态叶片前缘的变形要明显大于其他所有位置，前缘节点位移最大值为52.3 mm。

在鸟撞风扇叶片过程中，鸟体姿态主要影响鸟体被切割次数、鸟体被切割后各鸟块的质量、鸟块与叶片的接触位置以及接触面积。这些因素进一步影响叶片的应力响应及撞击过程对叶片的损伤大小。Y-135°、Y-270°、Y-315°、Z-135°及Z-315°撞击姿态对前叶根损伤最大，Z-135°撞击姿态对后叶根损伤最大，Y-270°撞击姿态对前缘的损伤最大。0°撞击姿态下，接触力、前叶根应力、后叶根应力及前缘位移为最小，对叶片前后叶根及前缘的损伤最小。在撞击过程中前叶根受到的应力要普遍大于后叶根，因而在撞击过程中前叶根相对于后叶根受到的损伤要更大。

4 结 论

本文采用真实绿头鸭SPH模型，模拟鸟撞整级旋转风扇叶片过程，对鸟撞位置和撞击姿态对叶片损伤的影响进行研究，得到如下结论：

(1) 鸟撞击高速旋转风扇的过程中会与多个叶片发生相互作用，被多个叶片切割为不同大小的鸟块，鸟体被切割后与叶片发生碰撞发生流变。与最大鸟块发生碰撞的叶片应力最大。鸟撞过程中，叶片最大等效应力出现位置为叶片前缘、后缘以及叶根。

(2) 从叶根到2/6叶高区域，撞击的位置越高，撞击力、叶根处应力以及叶片前缘应力越大。从2/6叶高到叶尖区域，撞击的位置越低，叶片受到的撞击力、叶根处及前缘处的应力越大，叶根及前缘受到损伤越大。

(3)Y-135°、Y-270°、Y-315°、Z-135°及Z-315°撞击姿态对前叶根损伤最大，Z-135°撞击姿态对后叶根损伤最大，Y-270°撞击姿态对前缘的损伤最大。0°撞击姿态下，接触力、前叶根应力、后叶根应力及前缘位移为最小，对叶片前后叶根及前缘的损伤最小。

(4) 在撞击过程中前叶根应力要普遍大于后叶根，因而在撞击过程中前叶根相对于后叶根受到的损伤更大。