刘建明，蒋向华，王 东，马永峰

（1.中航工业沈阳发动机设计研究所，沈阳 110015；2.北京航空航天大学能源与动力工程学院北京 100191）

0 引言

鸟类与飞机等飞行器之间的碰撞称为鸟撞。文献[1]显示，发动机是飞机较易受鸟撞击的部位。随着航空发动机技术的发展，在未来的新发动机研制中，由于叶尖切向速度增加，叶片厚度日趋变薄，叶片采用新结构（整体叶盘、整体叶环等）和先进材料（金属基、树脂基复合材料等），使发动机在推力增加、质量减轻的同时，叶片的抗鸟撞能力也可能受到影响，从而影响发动机的可靠性和安全性[2]。基于中国航空发动机研制的实际情况，前期采用大量的数值模拟，最后进行少量的试验验证是比较合理的发展模式[3]。目前，用于鸟撞数值模拟的大型通用有限元程序有Pam-Crash[1,4]、Abaqus[5]和 MSC.Dytran[6-7]等，还有专门针对鸟撞发展起来的程序，如 NOSAPM[8]、PW/WHAM[9]、MAGNA[10]等。目前，在采用流固耦合算法研究叶片鸟撞的问题上，普遍采用壳元划分叶片，而壳元在模拟非薄壁结构问题上存在困难。如果采用实体元划分叶片，在失效单元的处理及耦合面的推进上，流固耦合程序还缺乏完善的理论。

本文针对目前通用流固耦合算法在模拟实体元结构破坏上存在的不足，基于鸟撞铝板的试验结果，采用MSC.Patran软件建立了鸟撞实体元平板的有限元模型，同时研究和验证了实体元平板鸟撞数值模拟方法，进一步模拟了考虑破坏的实体元空心叶片的鸟撞过程。

1 流固耦合方法研究及试验验证

为了进行实体元空心叶片的鸟撞数值模拟，本文采取以实体元自由表面的节点作为壳元节点的方法，建立了覆盖在实体元表面的1层4节点壳元。壳元起辅助作用，用来形成封闭的耦合面，并向实体元传递力，而自身对原结构的影响较小。

为了验证方法的可行性，针对文献[11]所做的鸟撞铝板试验，建立了鸟撞铝板模型。其中铝板尺寸为410mm×500mm×10mm，4边固支。铝板用8节点6面体实体元划分，沿厚度方向等厚划分2层单元；在厚度方向的4个面上建立1圈哑元；以铝板实体元的上下自由表面的节点为公共节点，分别建立1层4节点四边形壳元；哑元和壳元形成包围实体元的、封闭的耦合面。所建立的鸟撞铝板模型如图1所示，图中欧拉域2由包围铝板的欧拉网格组成，欧拉域1由欧拉域2之外不发生鸟撞的欧拉网格组成。

图1 鸟撞铝板模型

铝板用分段线性塑性材料模拟，材料参数见表1，表中：ρ为密度；E为弹性模量；σy为静态屈服强度；μ为泊松比。取壳元与实体元的材料参数一致进行初始计算。鸟体采用线性流体模型，质量M=1.8 kg，密度ρ=928.15 kg/m3，体积模量 K=2200 MPa，用长径比为2∶1的二端半球中间圆柱实体来模拟，以30°角斜撞向铝板中央。试验测得了铝板中心点的法向位移随时间的变化情况，以及撞击开始1.24 ms后沿铝板宽度方向的测点的位移。

表1 铝板材料参数

1.1 壳元厚度h的确定

考虑壳元厚度h对精度的影响，对此只改变h进行计算。分别计算了 h=0.01、0.1、1、0mm（实体元表面未铺壳元）4种情况下铝板中心点的法向位移随时间的变化情况，并与试验值进行对比，如图2所示。从图中可见，当h=0.01mm，即实体元与h之比为102数量级时，计算值与试验值的误差较小。因此，选取h=0.01mm。

图2 不同壳元厚度时板中心点的位移-时间曲线

1.2 壳元弹性模量E s的影响

为了验证壳元材料的Es对计算精度的影响，只改变壳元的Es，分别计算Es=108000、36000MPa时铝板中心点的法向位移随时间的变化情况，并与试验值进行对比，如图3所示。

计算表明，3种壳元弹性模量所对应的铝板中心点的法向位移的差别很小。

根据以上计算，当壳元厚度取为0.01 mm，材料参数与实体元一致时，撞击开始1.24ms后沿铝板宽度方向的测点的位移的计算值与试验值的对比，如图4所示。从图中可见，二者趋势比较一致。

1.3 壳元最大塑性应变εs的影响

鉴于文献[11]所做的鸟撞试验，铝板没有发生破坏，为了考察实体元（铺壳）结构发生破坏时、壳元最大塑性应变εs对计算结果的影响，建立了鸟体撞击钛合金实体元（铺壳）平板有限元模型。平板尺寸为100 mm×120 mm×1.2 mm，4周固支，材料参数见表2，εp为实体元最大塑性应变。鸟体半径r=15 mm，长度L=50mm，以600m/s的速度正撞击平板中央。

表2 钛合金材料参数

当 εs=0.05、0.01、0.1时，对实体元（铺壳）平板的失效情况进行计算，并与εs=0.05、壳元划分的平板的基准计算结果进行对比。某时刻2种划分时平板的失效如图5所示。由图可以清楚地看到，当εs=εp时，实体元(铺壳)和壳元的失效情况比较一致。因此，计算时取 εs=εp较为合理。

计算结果表明，对于实体元划分且厚度不大的叶片，可以采用在实体元的自由表面铺1层壳元来进行叶片鸟撞数值模拟，且实体元与壳元的厚度之比为102数量级，2种单元的材料参数应一致。壳元的作用是传递鸟体与实体元之间的相互作用力。如果实体元表面没有铺壳元，目前的Dytran还不具备模拟实体元失效的功能。若实体元在鸟撞过程中失效，则该单元被删除。

图5 某时刻实体元(铺壳)和壳元平板的失效

2 叶片鸟撞有限元建模

采用Patran建立了8节点6面体空心叶片鸟撞模型。为节省计算时间，只建立了发动机中与鸟体有初始碰撞行为的第1级转子叶片模型；鸟体采用欧拉单元模拟。为了进行考虑叶片破坏的鸟撞数值模拟，参考第1节中所用的方法，以叶片的叶盆和叶背的实体元自由表面上的节点为公共节点，建立了覆盖在实体元表面的1层4节点壳元，h=0.01mm（实体元与h比为102数量级）；叶片的叶尖、叶根以及2叶片之间建立了1层哑元，如图6所示。鸟撞有限元模型如图7所示。在图7中，欧拉域2由自适应网格组成，由包围相邻叶片区域的欧拉网格以及包围叶片的欧拉网格组成，如图8所示。欧拉域1包含了叶片转动范围之外的区域。欧拉域1和欧拉域2通过哑元加以分隔。

设鸟体质量 M=0.35 kg，密度 ρ=930 kg/m3，体积模量K=2200 MPa。文献[12]的研究表明：圆柱体鸟体和两端半球中间圆柱鸟体对计算结果的影响较小，参照航空发动机适航规定，用长径比为2∶1的圆柱体代替鸟体，其速度V=102.8m/s。叶片受离心载荷作用，转速为12000 r/min，初始状态为稳定旋转，由Nastran的预应力分析计算获得。

图8 自适应网格

叶片材料模型见参考文献[3]，实体元和壳元材料一致，材料参数见表3。D、P为表征动态硬化的材料系数，D=100、P=10，EH为硬化模量。

表3 叶片材料参数

3 不考虑叶片失效的鸟撞数值模拟

鸟体的材料特性是鸟撞仿真分析的重点和难点，真实的鸟体有骨有血有肉，鸟体的本构方程很难描述[12]。本文分析了实体元表面壳元厚度，鸟体的密度、体积模量，以及叶片材料的弹性模量、屈服应力、硬化模量的变化对叶片鸟撞响应的影响。叶片和鸟体材料计算以第2节中的数据为基准，只改变单个参数，其他不变。

3.1 壳元厚度及单元材料的验证

图9 鸟撞瞬间

为了验证第1节中所用方法的准确性，计算了壳元厚度h=0.01、0.1 mm的2种情况下叶片应力峰值最大的单元应力（简称叶片应力）随时间的变化，并与未铺壳实体元叶片（h=0）的结果对比。某时刻的鸟撞瞬间如图9所示，不同h的叶片应力曲线如图10所示。从图9中可见，由于鸟撞产生的应力值超过了叶片的屈服应力，叶片产生局部塑性变形，巨大的相互作用力使得鸟体发生流变。从图10中可见，当h=0.01mm时，结果相对h=0的误差较小，为1.3%，从而验证了实体元与h之比为102数量级、材料一致的准确性。以下取h=0.01mm进行计算。文献[13]将鸟撞分为4个阶段：即初始撞击、压力衰减、恒定流动、流动终止。叶片应力响应趋势与鸟体正撞击刚性靶的应力趋势比较一致（图10）。

图10 不同壳元厚度的叶片应力-时间曲线

3.2 鸟体密度ρ的影响

选取鸟体密度 ρ=1300、930、500 kg/m3进行计算。3种ρ鸟体撞击下叶片应力曲线如图11所示。从图中可见，ρ对叶片响应的影响较大。叶片应力峰值随ρ的增大而增大，因为ρ越大，与叶片相撞的鸟体的质量越大，使得撞击能量越大，应力峰值也越大。

图11 不同鸟体密度的叶片应力-时间曲线

3.3 鸟体体积模量的影响

选取鸟体体积模量K=2200、6000、10000 MPa进行计算。3种K鸟体撞击下叶片应力曲线如图12所示。从图中可见，鸟体K对叶片响应的影响较小。

图12 不同鸟体体积模量的叶片应力-时间曲线

3.4 叶片弹性模量的影响

选取弹性模量 E=200000、112500、56250 MPa进行计算。3种E的叶片应力曲线如图13所示。从图中可见，E对叶片初始撞击的应力峰值影响较小，但对恒定流动的应力峰值的影响较大，E越大，恒定流动的应力峰值也越大。因此，在其他条件相同的情况下，E大的叶片能吸收更多的鸟体能量。

图13 不同叶片弹性模量的叶片应力-时间曲线

3.5 叶片屈服应力σs的影响

选取屈服应力σs=400、900、1300 MPa进行计算。3种σs的叶片应力曲线如图14所示。从图中可见，σs对叶片应力响应的影响较大，σs越大，叶片初始撞击和恒定流动的应力峰值也越大。

图14 不同叶片屈服应力的叶片应力-时间曲线

3.6 叶片硬化模量E H的影响

选取硬化模量 EH=14286、7000、3000 MPa进行计算。3种EH的叶片应力曲线如图15所示。从图中可见，叶片EH对叶片应力响应的影响较大，EH越大，初始撞击的应力峰值也越大，而恒定流动的应力峰值越小。

4 考虑叶片失效的鸟撞数值模拟

由于叶片的转速很高，撞击过程时间较短，因此，鸟体撞击叶片产生的应力和应变很大，从而易造成叶片的损伤和破坏。为了更加准确真实地模拟鸟体对叶片的撞击过程，给叶片添加最大塑性应变破坏准则进行计算，最大塑性应变εb=0.05。叶片破坏时的鸟撞瞬间如图16所示。从图中可见，叶片前缘由于受到撞击而造成破坏。

5 结论

（1）基于鸟撞铝板的试验结果，验证了在实体元表面铺一层壳元进行流固耦合计算的可行性。计算结果表明，当实体元与壳元的厚度之比为102数量级，且2种单元的材料参数一致时，铺1层壳元对计算结果的影响较小。

（2）鸟体密度对叶片应力响应的影响较大，叶片应力峰值随鸟体密度的增大而增大；鸟体体积模量对叶片应力响应的影响较小。

（3）叶片弹性模量对叶片初始撞击的应力峰值影响较小，但对恒定流动的应力峰值的影响较大，弹性模量越大，恒定流动的应力峰值也越大。叶片屈服应力、硬化模量对叶片应力响应的影响较大，屈服应力越大，叶片初始撞击和恒定流动的应力峰值也越大；硬化模量越大，初始撞击的应力峰值越大，而恒定流动的应力峰值越小。

（4）考虑叶片失效的计算进一步模拟了鸟撞对叶片的真实损伤过程。

[1]刘军,李玉龙,刘元镛.基于SPH方法的叶片鸟撞数值模拟研究[J].振动与冲击,2008,27（9）：90-93.

LIU Jun,LI Yulong,LIU Yuanyong.Numerical simulation study of bird-impact on a blade using SPH method[J].Journal of Vibration and Shock,2008,27（9）：90-93.(in Chinese)

[2]陈伟,高德平,尹晶.航空发动机叶片的鸟撞击损伤研究[J].燃气涡轮试验与研究,1998,11（4）：34-39.

CHEN Wei,GAO Deping,YIN Jin.Damage investigation of bird impact on areoengine blades[J].Gas Turbine Experiment and Research,1998,11（4）：34-39.(in Chinese)

[3]蒋向华,王延荣.采用流固耦合方法的整级叶片鸟撞击数值模拟[J].航空动力学报,2008,23（2）：299-304.

JIANG Xianghua,WANG Yanrong.Numerical simulation of bird impact on bladed rotor stage by fluid-solid coupling method[J].Journalof Aerospace Power,2008,23（2）：299-304.(in Chinese)

[4]刘军,李玉龙,徐绯.基于PAM-CRASH的鸟撞飞机风挡动响应分析[J].爆炸与冲击,2009,29（1）：80-84.

LIU Jun,LI Yulong,XU Fei.Dynamic response analysis of bird-impact aircraft windshields based on PAM-CRASH[J].Explosion and Shock Waves,2009,29（1）：80-84.(in Chinese)

[5]ZHU Shuhua,TONGMingbo,Wang Yuequan.Experiment and numerical simulation of a full-scale aircraft windshield subjected to bird impact[R].AIAA-2009-2575.

[6]刘建明,蒋向华.材料参数对叶片鸟撞动响应影响数值模拟[J].航空发动机,2010,36（5）：36-38.

LIU Jianming,JIANG Xianghua.Numerical simulation of blade material effect on dynamic response of bird impact on flat blade[J].Aeroengine,2010,36（5）：36-38.(in Chinese)

[7]刘建明,蒋向华,秦银雷等.采用流固耦合的实体元空心叶片鸟撞数值模拟 [J].航空动力学报,2010,25(10),2211-2216.

LIU Jianming,JIANG Xianghua,QIN Yinlei,et al.Numerical simulation of bird impact on solid-element hollow blades by using fluid-solid coupling method[J].Journal of Aerospace Power,2010,25(10),2211-2216.

[8]Storace A F,Nimmer R P,Ravenhall R.Analytical and experiment investigation of bird impact on fan and compressor blading[J].Journalof Aircraft,1984,21（7）:520-527.

[9]Teichman H C,Tadros R N.Analytical and experimental simulation of fan blade behaviour and damage under bird impact[J].Journal of Engineering for Gas Turbines and Power,1991,113（3）：582-594.

[10]McCarty R E,Hart J L.Validation of the MAGNA computer program for nonlinear finite element analysis of aircraft transparency bird impact[R].AD-A-140701,1983.

[11]白金泽.基于神经网络方法的鸟撞飞机风挡反问题研究[D].西安:西北工业大学,2003.

BAI Jinze.Inverse issue study of bird-impact to aircraft windshield based on neural network method[D].Xi'an:Northwestern Polytechnical University,2003.(in Chinese)

[12]朱书华,童明波.鸟体形状对飞机风挡鸟撞动响应的影响[J].南京航空航天大学学报,2008,40（4）:551-555.

ZHU Shuhua,TONG Mingbo.Bird shape sensitivity to dynamic response of bird strike on aircraft windshield[J].Journal of Nanjing University of Aeronautics and Astronautics,2008,40（4）：551-555.(in Chinese)

[13]Barber JP.Taylor H R.Wilbeck JS.Bird impact forces and pressures on rigid and compliant targets[R].AD-A-061313,1978.