直升机旋翼桨叶鸟撞动态响应计算

林长亮，王益锋，王浩文，陈仁良，尚晓冬

(1.南京航空航天大学 直升机旋翼动力学国家级重点实验室，南京 210016;2.中国商用飞机有限责任公司 上海飞机设计研究院，上海 200232;3.清华大学 航天航空学院，北京 100084;4.哈尔滨飞机工业集团有限责任公司 飞机设计研究所，哈尔滨 150066)

鸟撞轻则使机体损伤，重则会造成灾难性后果，直接威胁人员的生命安全。Dolbeer等［1］在第8届鸟撞会议上指出“直升机鸟撞事故主要发生在距地面高度600 m以下，风挡和旋翼是鸟撞事故发生的主要部位”。低空飞行是直升机的显著使用特点，因此直升机发生鸟撞的可能性很大，并且随着飞行速度的提高，鸟撞事故的危害性也在逐步加大。旋翼是直升机的升力面、推力装置和操纵面，是区别于固定翼飞机的主要特征。在直升机飞行过程中旋翼一直处于高速旋转的状态，一旦发生鸟类撞击，就会对飞行安全会造成严重的威胁，是抗鸟撞研究中需要着重解决的关键技术。

对于鸟撞问题的研究，在飞机鸟撞方面，国内外学者作了大量的研究工作。Barber等［2-3］采用10%孔隙率的明胶代替真鸟进行试验研究，发现鸟撞过程可被描述成一个非恒定的流体动力学过程。Meguid等［4］通过计算研究了鸟体形状对计算结果的影响，发现鸟体和目标的最初接触面积对接触力的峰值有巨大影响。刘军等［5-6］进行了鸟撞平板试验研究，使用神经网络方法对试验中的鸟体参数进行反演计算。陈伟和关玉璞等［7-9］在发动机叶片鸟撞的载荷模型、瞬态响应计算、试验方法等方面开展了大量工作。Zhu等［10］对鸟撞飞机全尺寸风挡进行了试验研究和数值模拟。万小朋、李京菁等［11-12］对飞机机翼前缘的抗鸟撞性能进行了研究。

在直升机鸟撞研究方面，仅有美国西科斯基公司在S-92直升机抗鸟撞设计中，对垂尾前缘、尾桨鸟撞进行了动力学仿真，并开展了鸟撞试验验证研究［13-14］。目前，我国对直升机旋翼鸟撞问题进行的研究很少。王益锋等［15］采用Hertz接触理论处理局部弹性变形，运用Hamilton原理建立了直升机旋翼桨叶的弹性碰撞的动力学模型。温海涛、林长亮等［16-17］使用有限元软件对直升机主桨叶鸟撞的过程进行了数值模拟。

直升机旋翼在旋转过程中，桨叶载荷工况非常复杂，除桨叶自身重力外，不仅受到气动升力作用，还承受离心载荷。同时，由于采用桨毂结构形式的不同，桨叶的挥舞、摆振、扭转运动还存在不同的非线性耦合。采用有限元软件很难模拟出实际的旋翼工作特点。利用旋翼综合气弹分析程序LORA［19-20］，考虑旋翼桨毂结构特点和旋翼旋转过程中桨叶的挥舞、摆振、扭转运动的非线性耦合，求解出直升机旋翼桨叶在飞行过程中的稳态响应，以此作为鸟体撞击的初始状态，然后采用非线性流-固耦合算法，考虑载荷与响应之间耦合对桨叶响应的影响，利用Newmark数值积分的方法求解桨叶的动态响应，并分析了鸟体速度、质量、撞击位置、桨叶根部铰约束、离心力等系统参数对桨叶响应的影响，从而为直升机旋翼的抗鸟撞设计提供一定的理论依据。

1 桨叶鸟撞动力学建模

在桨叶鸟撞过程中，桨叶是可变形柔性靶体，其变形会影响撞击载荷的大小与分布。为准确预估桨叶鸟撞击响应，因此采用流-固耦合算法［8］计算桨叶鸟撞的动态响应。

1.1 流－固耦合鸟撞载荷模型

在鸟撞击载荷作用期间，在ti时刻，假设鸟体以初始速度，入射角θi(θi=90°时为正撞击)撞击桨叶，ti+1=ti+Δt时刻，桨叶发生变形，入射角变为 θi+1，剖面入射角示意图如图1所示。

图1 剖面撞击示意图Fig.1 Profile sketch of bird impact on the blade

入射角θi+1可以表示为:

以Δt时间内消耗的鸟体质量Δmi+1为研究对象，应用动量定理可得:

式中:ρ和Ai+1分别为鸟的密度和横截面积。撞击载荷的作用时间由撞击时鸟体在撞击过程中消耗的长度Si决定:

当Si等于鸟体的长度L时，撞击过程结束。

图2为面积随鸟体消耗长度的示意图，鸟体横截面积呈三角形形式变化，在S/L=0.2处面积为最大，a、b分别为椭圆形撞击区域的短轴与长轴，当正撞击时，a=b=1.5R，R为鸟体的初始半径。

图2 面积A随长度的变化Fig.2 Area“A”variation with length

1.2 载荷空间分布

载荷类型采用分布载荷。当斜撞击时，撞击区域为一个椭圆，方程为(Ysinθ)2/R2+Z2/R2=1。式中Y、Z分别沿椭圆的长、短轴方向，碰撞载荷在撞击区域上的平均分布力为:

鸟撞击载荷初始冲击压力虽然很高，但持续时间极短，而恒定流动压力虽然不高，但持续时间较长，鸟撞载荷的冲量主要是在此阶段传递给桨叶的［18］。因此，一般在建立鸟撞击载荷模型时根据恒定流动压力的分布情况来选择载荷的空间分布形式。本文中载荷空间分布形式为压力系数沿长轴线性分布，沿短轴均匀分布，载荷空间分布如图3所示。具体分布形式可以表示为:

沿椭圆短轴a，

沿椭圆长轴b，

图3 载荷空间分布示意图Fig.3 Load spatial distribution map

1.3 桨叶鸟撞动力学方程

根据Hamilton原理，引入Chopra有限元法对桨叶单元进行离散，组集动能、应变能、空气动力及撞击载荷产生的虚功，得到基于广义力形式的桨叶非线性动力学隐式微分方程:

式中:T、E、A、P表示动能、应变能、气动力以及撞击载荷引起的广义力，在计算过程中，由于撞击过程时间持续较短，假设撞击过程中入流和桨叶迎角没有发生变化。

1.4 桨叶动态响应求解步骤

本文采用有限转动梁理论处理桨叶的变形［19-20］，在翼型气动力计算时采用Leishman-Beddoes二维非定常和动态失速模型计算出桨叶剖面升力、阻力和俯仰力矩，旋翼诱导速度计算采用Glauert入流改进模型，采用Chopra的15自由度梁元对桨叶进行离散，根据Hamilton原理得出非线性隐式旋翼动力学方程，利用隐式Newmark数值算法对动力学方程进行了求解得出桨叶的稳态解，在每个时间步NTIME上，先预估得到叶片的响应(位移、速度和加速度)值，接着利用ASSMK模块求解桨叶的响应，当桨叶单元为NEIMP撞击单元时，进入ADDIT模块，计算撞击载荷值和作用时间，并将载荷和力矩施加到撞击单元上，这样通过ADDIT模块完成在桨叶动力学模型上施加鸟体分布载荷的过程，然后利用NMARK1模块计算出桨叶响应值，并且在每个时间步上进行迭代，直到响应收敛，得出这个时间步上的桨叶响应值，再将计算结果返回程序进行下一个时间步的计算。如此反复，直到撞击仿真时间结束为止。流程图如图4所示。

图4 桨叶响应计算流程图Fig.4 The calculation flowchart of Blade response

2 分析方法验证

由于目前国内外对于直升机旋翼鸟撞问题的研究比较少，还没有公开的详细试验数据作为验证依据。并且在理论计算方面，直升机旋翼鸟撞问题属于软体瞬态冲击旋转梁范畴，而现有的国内外文献主要集中于刚性质量与梁的碰撞研究，对软体与梁的碰撞问题的研究非常少。因此，为了能够验证本文计算方法的可靠性，主要从以下两个方面进行了验证。

2.1 鸟体载荷模型验证

为验证本文采用鸟体载荷模型的正确性，针对文献［21］中平板叶片鸟撞试验进行了验证。文献［21］给出了质量为0.161 5 kg的鸟体垂直撞击平板叶片根部的高速摄影结果及叶尖Z向位移测量结果。在计算中，鸟体直径为50 mm，长度为80 mm，密度为1 028 kg/m3，速度为156.1 m/s垂直撞击平板叶片，撞击中心距叶片根部为60 mm。叶片材料参数见表1所示。叶片采用8节点壳单元进行划分，叶片有限元网格划分如图5所示。

表1 叶片材料参数Tab.1 Blade material parameters

图5 叶片有限元模型Fig.5 The Finite element model of blade

图6 叶尖位移时间历程曲线Fig.6 Tip displacement time history curve

图6给出了本文与试验的叶尖位移对比结果。可以看出，计算结果中含有高阶谐波成分，没有试验值曲线光滑，主要是由于本文没有对计算结果进行滤波处理。试验中，叶尖的最大位移为70.5 mm，本文计算结果中的最大位移为74.7 mm，误差为5.9%，总体上讲，计算结果与试验是比较吻合的，从而验证了本文所采用鸟撞载荷模型的有效性。

2.2 旋翼动力学模型验证

为验证旋翼动力学模型在计算瞬态响应时的正确性，本文选取 Keller实验作为数值算例［22］。Keller所采用的试件为与美国H-46运输直升机桨叶1/8Froude数相似的模型桨叶，桨叶长 1.006 m，重0.998 kg。实验装置如图7所示。整个实验过程为:桨叶初始处于静止状态，根部固支，挥舞限动角为0°，初始挥舞角为2°至9.7°，在自重作用下桨叶有初始位移。突然释放根部约束，桨叶下坠，当挥舞角到达0°时，根部与挥舞限动块发生碰撞。

图8～9分别给出了初始挥舞角为4°时，桨尖位移和距桨叶根部距离20%处桨叶上应变的变化曲线。与实验相比，本文桨尖最大位移误差为5.1%，20%处最大应变误差为6.56%。通过与试验值的比较验证了本文建模方法的正确性。

图7 Keller桨叶扬起下坠实验装置Fig.7 Configuration of blade droop stop impact test

图8 桨尖位移随时间的变化图Fig.8 Blade tip displacement changing with time

3 桨叶鸟撞动态响应计算

旋翼翼型为NACA0012，直径为3 m，弦长0.3 m，桨叶挥舞刚度3.78×105N·m2，摆振刚度9.45×106N·m2，扭转刚度 3.78 ×106N·m2，拉伸刚度 1.26 ×109N·m2，桨叶线密度48.6 kg/m。在鸟撞击瞬态响应研究中，把模拟鸟看做长径比为2的圆柱体［7，10］，鸟体直径0.05 m，质量 0.17 kg，鸟体以 120 m/s的相对速度与桨叶发生撞击，鸟体入射角为60°，撞击位置取桨叶弹性轴上径向位置x/L=0.4处，分析时间为0.04 s。

首先采用流-固耦合和流-固非耦合两种算法分别计算了桨叶的动态响应。图10给出了桨叶鸟撞载荷曲线。当考虑载荷与响应之间的耦合效应时，碰撞载荷峰值为1.91×105N，载荷作用时间为 1.06 ms，当不考虑耦合响应时，碰撞载荷峰值为2.0×105N，载荷作用时间为1.0 ms;图11给出了桨叶碰撞位置位移随时间的变化曲线，可以看出，考虑耦合效应的局部弹性变形要稍小于不考虑耦合效应的局部变形，随着时间的推移，两种算法几乎同时到达向下位移峰值。通过以上比较可以看出，采用流固耦合算法时，虽然载荷作用时间延长，但是由于桨叶变形起到了卸载作用，因此使得撞击载荷、桨叶的整体变形要小于不考虑耦合效应时的结果;同时可以看出，算法的选择对桨叶到达最大变形的时间以及桨叶的振动形式影响不大。

图9 x/L=0.2处应变随时间的变化图Fig.9 The strain changes with time of the position x/L=0.2

图10 碰撞载荷曲线Fig.10 The impact load curve

图11 撞击位置位移曲线Fig.11 The displacement curve of impact position

4 系统参数对桨叶响应的影响

4.1 速度影响分析

鸟体的质量和速度，决定了鸟撞载荷冲量的大小，为研究速度对桨叶响应的影响，保持其他参数不变，单独改变鸟体的初始速度，对桨叶的响应进行了分析.图12给出了不同撞击速度的载荷曲线，图13为不同速度下的碰撞位置位移曲线。可以看出随着鸟体速度的增加，碰撞载荷变大，载荷作用时间减小，桨叶撞击位置局部弹性变形和整体弹性变形增大。

4.2 质量影响分析

分别选取质量0.15 kg，0.17 kg，0.19 kg，采用流 -固耦合算法，对桨叶鸟撞响应进行了分析，计算结果如图14、15所示。

可以看出鸟体质量对桨叶的鸟撞响应影响较明显，响应的差异基本与质量成比例关系，单独增加鸟体质量时，碰撞载荷增加，撞击位置位移增大。

4.3 撞击位置影响分析

为比较桨叶对不同加载位置的响应，分别计算了三组算例。撞击位置分别取在x/L=0.1，x/L=0.4处，x/L=0.7处。计算结果如图16、17所示。

图12 不同速度下的碰撞载荷曲线Fig.12 The impact load curve with different speed

图13 不同速度下的碰撞位置位移曲线Fig.13 The displacement curve of impact location under different speed

图14 不同质量的碰撞载荷曲线Fig.14 The impact load curve with different quality

图15 不同质量的碰撞位置位移曲线Fig.15 The displacement curve of impact location under different quality

图16 不同撞击位置的碰撞位置位移曲线Fig.16 The displacement curve of impact location with different impact position

图17 不同撞击位置的桨尖位移曲线Fig.17 The displacement curve of blade tip with different impact position

图18 不同约束的碰撞位置位移曲线Fig.18 The displacement curve of impact location with different constraint condition

图19 不同约束的桨根挥舞弯矩曲线Fig.19 The flap moment curve of blade root with different constraint condition

图20 不同转速的碰撞载荷曲线Fig.20 The impact load curve with different rotation speed

从图16、17碰撞位置位移、桨尖位移时间历程中可以看出，在撞击载荷作用期间，桨叶局部弹性变形相差不大，但是随着时间的推移，当撞击点离桨根越远时，碰撞位置、桨尖位移越大。对于x/L=0.1处的撞击，碰撞位置处的位移振动后迅速返回平衡位置，但桨尖位移存在明显的滞后，这是由于从撞击力在撞击处产生加速度到桨尖呈现位移，中间经过了桨叶的弹性变形及加速度两次对时间的积分过程，因此桨叶尖位移出现了滞后。

4.4 根部约束影响分析

为了研究根部铰约束对桨叶响应的影响，保持其他参数不变，单独改变根部约束形式，采用铰接式同时根部施加角弹簧，角弹簧刚度为6.0×102N/rad。从图18碰撞位置位移时间历程中可以看出，由于根部铰约束刚度的降低，桨叶刚性运动影响增强，整体弹性位移增大。从图19桨根挥舞弯矩随时间的变化曲线可以看出，虽然桨叶整体位移增大，但是由于桨叶刚性运动的影响，桨根挥舞弯矩减小。

4.5 旋翼转速影响分析

当桨叶旋转时，离心力对桨叶的响应有很大的影响。保持其他参数不变，单独改变旋转角速度，对桨叶鸟撞响应进行了计算。计算结果如图20、21所示。

从图20碰撞载荷变化的时间历程可以看出，当旋转角速度增加时，碰撞载荷增大，这是由于旋转角速度的增加导致鸟体相对桨叶的撞击速度增大，但是迎角相对减小。

图21 不同转速的碰撞位置位移曲线Fig.21 The displacement curve of impact location with different rotation speed

从图21碰撞位置位移时间历程可以看出，当旋转角速度增加时，由于离心力的作用，桨叶的局部刚度和整体刚度得到增强，桨叶的局部弹性变形变小，碰撞位置位移减小，当载荷作用结束后，碰撞位置位移明显减小。

5 结论

(1)以旋翼综合气弹分析程序LORA为基础，采用流固算法，建立了直升机桨叶鸟撞动力学方程，并利用数值积分求解出了鸟撞桨叶的动态响应，为研究直升机旋翼鸟撞问题提供了一种参考方法。

(2)采用流固耦合算法时，作用时间延长，但碰撞载荷以及桨叶局部弹性变形相对要小一些，同时发现算法的选择对桨叶到达最大变形的时间以及振动形式影响不大。

(3)研究了鸟体速度、质量、撞击位置、桨叶根部约束和离心力等参数对桨叶动态响应的影响，从而为直升机桨叶抗鸟撞设计提供一些理论依据。

［1］Dolbeer R A，Wright S E，Cleary E C.Bird strikes to civil helicopters in the United States［C］.1990-2005.the 8th meeting ofBirdStrikeCommitteeUSA/Canada，2006:45-50.

［2］ Barber J P，Taylor H R，Wilbeck J S.Characterization of bird impact on a rigid plate:partⅠ［R］.AFFDL-TR-75-5，AD A021142，1975:1-44.

［3］ Barber J P，Wilbeck J S，Taylor H R.Bird impact forces and pressures on rigid and compliant targets［R］.AFFDL-TR-77-60，AD A 061313，1978:1-78.

［4］Meguid S，Mao R，Ng T.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.

［5］刘 军，李玉龙.鸟体本构模型参数反演Ⅰ:鸟撞平板试验研究［J］.航空学报，2011，32(5):802-811.

LIU Jun， LIYu-long. Parameters inversion on bird constitutive model PartⅠ:study on experiment of bird striking on plate［J］.Acta Aeronautica et Astronautica Sinica，2011，32(5):802-811.

［6］刘 军，李玉龙，石霄鹏，等.鸟体本构模型参数反演Ⅱ:模型参数反演研究［J］.航空学报，2011，32(5):812-821.

LIU Jun，LI Yu-long，SHI Xiao-peng，et al.Parameters inversion on bird constitutive model PartⅡ:study on model parameters inversion［J］.Acta Aeronautica et Astronautica Sinica，2011，32(5):812-821.

［7］陈 伟，关玉璞，高德平.发动机叶片鸟撞击瞬态响应的数值模拟［J］，航空学报，2003，24(6):531-533.

CHEN Wei， GUANYu-pu， GAODe-ping. Numerical simulation of the transient response of blade due to bird impact［J］.Acta Aeronautica et Astronautica Sinica，2003，24(6):531-533.

［8］陈 伟，高德平.载荷与响应耦合下叶片鸟撞击响应分析［J］，航空动力学报，1998，13(1):93-96.

CHEN Wei，GAO De-ping.Transient responses of blades to bird impact coupled with deformations of blade［J］.Journal of Aerospace Power，1998，13(1):93-96.

［9］关玉璞，陈 伟.粒子分离器涡流叶片鸟撞击损伤试验［J］，航空动力学报，2007，22(12):2094-2095.

GUAN Yu-pu，CHEN Wei. Experimental study on bird damage of vortex vane in the particle separator［J］.Journal of Aerospace Power，2007，22(12):2094-2095.

［10］ Zhu S H，Tong M B，Wang Y Q.Experiment and numerical simulation of a full scale aircraft windshield subjected to bird impact［R］.AIAA，2009:2009-2575.

［11］万小朋，赵美英.基LS-DYNA的飞机机翼前缘抗鸟撞分析［J］.西北工业大学学报，2007，25(2):285-289.

WAN Xiao-peng，ZHAO Mei-ying.A method for calculating anti-bird impact capability of an aircraft wing using LS-DYNA Software［J］. Journal of Northwestern Polytechnical University，2007，25(2):285-289.

［12］李京菁，赵美英.机翼前缘抗鸟撞设计仿真方法验证［J］.机械科学与技术，2011，30(10):1757-1760.

LI Jing-jing，ZHAO Mei-ying.Verification of the simulation methods to bird impact resistance design of wing leading edge［J］.Mechanical Science and Technology for Aerospace Engineering，2011，30(10):1757-1760.

［13］ Dobyns A，Federici F，Young R.Bird strike analysis and test of a spinning S-92 tail rotor［C］.American Helicopter Society 57th Annual forum.May 9-11，2001.

［14］ Dobyns A.Bird strike analysis of S-92 vertical tail cover using DYTRAN［R］.American Helicopter Society.October 7，1998.

［15］王益锋，王浩文，高 正.直升机旋翼桨叶的弹性碰撞动力学建模［J］.航空动力学报，2009，24(5):2046-2050.

WANG Yi-feng，WANG Hao-wen，GAO Zheng.Dynamic model of elastic impact of helicopter rotor blades［J］.Journal of Aerospace Power，2009，24(5):2046-2050.

［16］温海涛，关玉璞，高德平.直升机主桨叶的鸟撞有限元数值模拟［J］.航空动力学报，2009，24(5):1150-1158.

WEN Hai-tao，GUAN Yu-pu，GAO De-ping.Numerical simulations of bird impact on helicopter rotor blades［J］.Journal of Aerospace Power，2009，24(5):1150-1158.

［17］林长亮，王浩文，陈仁良.采用流固耦合方法的直升机桨叶鸟撞数值模拟［J］.科学技术与工程，2012，12(1):1-6.

LIN Chang-liang， WANG Hao-wen， CHEN Ren-liang.Numerical simulation of bird impact on the helicopter blade by fluid-structure coupling method［J］.Science Technology and Engineering，2012，12(1):1-6.

［18］尹 晶.鸟撞击的载荷因素对叶片响应的影响［J］.航空动力学报，1992，7(1):51-54.

YIN Jing.Effect of some load factors of bird impact on blade reponse［J］.Journal of Aerospace Power，1992，7(1):51-54.

［19］王益锋，王浩文，高 正.一种复合材料柔性梁的响应计算方法［J］.南京航空航天大学学报，2008，40(2):180-184.

WANG Yi-feng，WANG Hao-wen，GAO Zheng.Method for calculating response of composite flexible beams［J］.Journal of Nanjing University of Aeronautics and Astronautics，2008，40(2):180-184.

［20］王浩文，高 正，郑兆昌.前飞状态下直升机旋翼系统气弹响应及稳定性分析［J］.振动工程学报，1999，12(4):521-528.

WANG Hao-wen， GAO Zheng， ZHENG Zhao-chang.Aeroelastic response and stability of helicopter rotor blades in forward flight［J］.Journal of Vibration Engineering，1999，12(4):521-528.

［21］文颖娟.叶片鸟撞击响应分析模化技术与验证研究［D］.南京:南京航空航天大学，2006.

［22］Keller J A，Smith E C.Experiment/theoretical correlation of analysis for helicopter rotor blade/droop stop impacts［R］.AIAA-97-1094，1997.