航空发动机宽弦空心风扇叶片以其良好的气动效率、轻量化设计和低噪的特点，显著提升了航空动力系统的经济性和环保性，对航空工业的发展具有重要意义。但宽弦空心风扇叶片在高速旋转过程中，其独特的薄壁特征和轻量化结构容易产生叶片振动问题，且该类型叶片在遭遇外物冲击时，面临着比传统实心叶片更高的结构损伤风险，直接威胁到飞机的安全性与可靠性。在各种外物高能冲击载荷中，鸟类撞击航空发动机风扇叶片的频率是最高的。当航空发动机风扇叶片受到鸟撞冲击后，轻微时可能造成航班颠簸，严重时会造成发动机空中停车，甚至可能会使风扇叶片掉块甩出进而危及飞行安全。我国的生态环境复杂且多样，不同类型的鸟类活动范围十分广阔，在军用航空和民用航空领域均面临着鸟撞航空发动机风扇叶片的危险。风扇叶片在运行过程中不仅需要避免发生共振，同时也需要具备良好的抗鸟撞性能，因此研究宽弦空心风扇叶片的抗鸟撞冲击响应具有重要意义。

在当前鸟撞风扇叶片研究中，国内外学者一般使用简化人工鸟体代替真实鸟体冲击风扇叶片，如球体、椭球体、圆柱体以及两端半球圆柱体等。Liu等使用圆柱鸟体对风扇叶片进行撞击实验仿真，研究了不同质量鸟体对叶片宏观损伤模式的影响；Zhou和Wang等采用两端半球圆柱鸟体研究了叶片受到鸟撞载荷后损伤行为和气动性能的影响。但由于过于简化的鸟体模型无法真实有效反映航发风扇叶片所受到的鸟撞载荷，因此近年来多数学者开始进行真实鸟体的研究。Li等构建了帆背潜鸭的几何简化鸟体模型，通过数值模拟方法研究了不同鸟撞方向对叶片损伤的影响；Wu等通过CT扫描获得野鸭的SPH粒子模型，开展了对不同鸟撞姿态、位置对风扇叶片的损伤研究，且通过Wilbeck实验证明了真实鸟体与常规人工鸟体在叶片撞击响应方面同样具有合理性。但CT扫描等方法对实验设备和场地有一定的要求，因此需要一种兼顾效率和保真的真实鸟体建模方法。值得注意的是，当前的风扇叶片逐渐向高效率轻量化发展，因此大多数都采用了空心设计。在鸟撞叶片研究中大多数学者仅分析了实心风扇叶片抗冲击特性，同时在鸟体撞击叶片时也常常忽略了气动力的存在。

如图1所示为本研究的研究思路整体框架，通过流固耦合方法下的预应力分析确定空心风扇叶片的最大转速，之后在流固耦合基础上对叶片进行振动分析，获得其合理转速范围。根据本研究所提出的近真实鸟体高效建模方法获得经过实验结果验证的有效近真实鸟体模型，在此基础上进行鸟撞空心风扇叶片仿真，评估其抗鸟撞冲击响应特性。

图1 本研究的研究思路整体框架

针对当前研究难以高效获取用于鸟撞仿真的近真实鸟体模型这一难题，本研究结合鸟类生物特征姿态图片提出了近真实鸟体高效建模方法。通过鸟撞平板实验结果验证了本研究所提出的近真实鸟体构建方法的有效性，相比于传统简化鸟体可以更加真实反映出鸟撞载荷。采用流固耦合方法，分析了宽弦空心风扇叶片的振动特性，并得到了风扇叶片的合理转速工作区间范围。最后，针对1.8kg真实灰鹤鸟体在相对速度100m/s以头部正面撞击转速为3344r/min的风扇叶片这一工况开展了叶片抗近真实鸟体撞击动态响应研究，以期为空心风扇叶片抗鸟撞设计提供参考。

1.1 鸟体有限元模型构建

当前，简化的规则体结构如球体、圆柱体等无法有效表征不同种类鸟体真实的生物结构特征差别，而这些结构差异又会对仿真结果有着较大的影响。如图2展示了本研究所提出的近真实形状鸟体几何模型构建流程，本研究使用包含鸟体特定姿态和基本视图的摄影图片，采用图像长度测量工具对摄影图片中关键部位的特征参数进行测量。通过获得鸟体各个关键部位内特征参数和不同关键部位内特征参数之间的比例关系，并基于B样条曲线绘制各关键部位简化后的外轮廓线，采用旋转、扫掠等的方法以及等比例缩放方法绘制出真实鸟体结构实体模型。将鸟体分为喙、头部、颈部、躯干、翅膀和腿部六大关键部位，并根据对应的外形轮廓特征简化为包括圆锥、圆柱等在内的旋转几何体结构，并规定了每个关键部位中的特征参数。在此基础上，从而建立了只需要提供包含鸟体特定姿态和基本视图的摄影图片就可以高效绘制出任意种类、任意质量的真实鸟体实体结构模型方法。

最后本研究选择1.8kg的真实灰鹤鸟体进行后续的鸟撞叶片仿真，将几何模型导入ANSYS Workbench进行有限元网格划分，基于Lagrange单元质心位置生成对应的SPH粒子鸟体模型如图2所示。SPH粒子间距设置为4mm（略小于风扇叶片网格尺寸），既可以保留鸟体的生物特征，又确保计算时长在可接受范围内，此时SPH粒子数量为28200。

图2 近真实形状鸟体几何模型和有限元模型构建与流程图

1.2 叶片有限元模型构建

如图3所示为航发风扇叶片结构域和流体域网格模型构建流程。根据叶片实物模型构建出几何模型，并得出单个流体域和结构域几何实体。其中结构域叶片实体为空心叶片，内部包含4条加强筋，为方便划分网格这里对加强筋和空心叶片连接处设置为绑定接触。如图4所示为有限元模型的网格无关性验证，其中结构有限元主要模拟叶片在4000r/min离心力作用下的叶片最大应力，图中可以看到当网格尺寸为5mm时可以满足计算需求，最终选定5mm作为结构域网格尺寸，划分后的叶片单元数为15254。

针对流体域构建单个风扇叶片的周期性流道，为简化计算设置叶尖间隙为0，即不考虑叶尖泄漏流的影响。由于流场几何形状包含复杂曲面，选用多面体网格对流场进行划分，且尽可能使网格方向与气流方向保持一致。设定边界层以捕捉近壁面流动，同时满足SST湍流模型计算精度对网格尺寸的要求。流场网格无关性验证如图4，此时风扇转速设置为4000r/min，为了提高计算效率并减小计算离散误差，如图可以看到设置流场网格尺寸为3mm时可以保持计算结果稳定，此时划分后的流体域单元数为1061663。

图3 空心风扇叶片结构与流体域有限元模型构建

图4 有限元模型的网格无关性验证

1.3 材料模型及参数

鸟体软组织在超过100m/s的速度下撞击到硬物会表现出类流体特性，这主要是由于鸟体材料在撞击过程中所产生的应力远远大于材料强度，因此一般被视为柔性体。本研究使用LS-DYNA软件中的*MAT Null材料模型，结合多项式状态方程进而表征高速下鸟体的变形行为，状态方程如式1所示，鸟体状态方程参数如表1所示。

式中C₀~C₆为状态方程系数，E为单位体积内能，μ为相对体积。

表1 鸟体状态方程参数

风扇叶片转子为钛合金材料(TC4)，本研究选择Johnson-Cook弹塑性材料模型描述其动态本构关系，其强度方如式2所示，对应的材料参数如表2所示。

表2 钛合金TC4本构模型参数

在传统简化鸟体中使用圆柱鸟体一般较多，该类型鸟体与本研究真实鸟体的几何外形存在着巨大的差异，高速撞击风扇叶片时造成的应力分布和损伤面积也会不同，因此需要对比分析鸟体冲击过程的动态响应。为了验证高速鸟撞过程中的柔性体以类流体状态进行变形的过程，一般常使用鸟撞平板实验对有限元方法进行验证。

为确保本研究所提出的真实鸟体的有效性，本研究通过Wilbeck鸟撞平板试验数据进行验证，Wilbeck通过空气压缩罐发射真实鸟体炮弹撞击到圆盘平板，并记录撞击中心处的压力变化。本研究针对鸟撞平板过程进行数值模拟验证，建立厚度为6cm，直径为60cm的金属平板模型，其中平板的材料更改为钛合金，以使得圆盘和风扇叶片的材料保持一致。真实鸟体通过几何模型缩放控制使其重量为1kg，同时圆柱鸟体长径比设置为2，两种鸟体重量相等。设置初始时两种鸟体距离平板均为1mm，速度为116m/s，撞击位置为平板中心。

图5展示了两种鸟体在初始撞击速度116m/s下的撞击平板中心压力-时间曲线和Wilbeck实验结果的对比。图中可以看到圆柱体鸟体模型对平板中心造成的的第一次峰值压力远远大于实验数据，而真实鸟体模型比圆柱体鸟体模型要更接近Wilbeck实验结果，尤其是在初始撞击阶段。这主要是由于本研究所构建的真实鸟体模型结构特征更加符合鸟体生物结构，而圆柱鸟体主要是对躯干部分的简化，这导致鸟体载荷过大。同时图5展示了两种鸟体在平板等效应力最大时的平板上应力分布情况，在应力峰值方面，圆柱鸟体对平板所造成的等效应力峰值大于真实鸟体，分别为66MPa和45MPa，前者比后者大1.5倍左右，这是由于在撞击时真实鸟体的颈部起到了缓冲作用。与此同时，在圆盘的应力分布方面，圆柱体人工鸟体对平板造成的应力变化更集中并产生明显的应力波，而真实鸟体则更加分散未显示明显的应力波，这主要是由于撞击时圆柱体鸟体与圆盘的接触面远大于真实鸟体。

图5 平板撞击中心压力曲线

为对比两种鸟体的类流变特性，图6给出了平板受到撞击后两种鸟体的动态铺展过程。图中可以看到两种鸟体都与Heimbs所总结的四个撞击变形阶段相对应，首先是初始撞击阶段鸟体最初接触到平板并产生应力波，该阶段的应力峰值被称为Hugoniot应力，随着应力波向鸟体后端传递，在鸟体侧面逐渐形成逆压梯度，进而加速SPH粒子向外加速，并释放卸载波进入冲击衰竭阶段。经过卸载波的传递后，鸟体变形逐渐以稳定的速度和压力进行，因此进入恒定流动阶段，而随着应力波不断被卸载波削弱，鸟体变形逐渐结束，到达撞击结束阶段。

在相同质量下，真实鸟体的长度要大于圆柱鸟体，因此撞击圆盘的持续时间也更长，进而造成真实鸟体的SPH模型飞溅程度大于圆柱体鸟体。真实鸟体与传统简化鸟体在冲击方面造成的结果截然不同，本研究所提出的真实鸟体可以更加准确表征鸟撞载荷，并通过Wilbeck实验数据验证了鸟体的有效性。

图6 鸟体撞击平板过程：

(a)理论变形过程;(b)圆柱鸟体变形过程;(c)近真实鸟体变形过程

3.1 流固耦合作用下叶片预应力分析

风扇叶片在正常运转中会受到以自身离心力与气动力为主的多种载荷以及结构约束的作用，经长时间运作后，叶片的力学性能和表现会降低，严重时甚至会导致叶片断裂。因此，叶片应当在能承受自身离心力载荷与气动力载荷的转速下运转，故本研究首先通过流固耦合强度分析确定叶片能承受的最大转速，并在此基础上进行叶片的振动特性分析。

在风扇叶片表面气动力计算中，流场的基本控制方程采用定常可压缩雷诺平均(RANS)N-S方程。为了提高对逆压梯度边界层流动以及分离状态的计算精度，本研究采用了k-ωSST(shear stress transport)两方程模型作为湍流RANS方程组封闭模型。

航发风扇在工作时，通过单个叶片的空气流量可达14kg/s，此时气动力对叶片的影响已经不可忽略，因此需要分析叶片表面气动力的分布以及对叶片的影响。经Fluent计算，风扇叶片在4000r/min转速时所受气动压强分布如图7所示。当空气流经叶片表面由于压力面的叶根部分几何形状下凹，同时该叶片的扭转角较大，叶片前缘正对来流方向经流空气量较大,受到的气动压强较大,达到了1.10×10⁵Pa。吸力面叶片中部位置由于几何形状上凸，导致该处空气流速快，因此吸力面叶中部位所受的气动压强较小，仅为9.90×10⁴Pa左右。

图7 风扇叶片表面气动压强分布

在ANSYS Workbench中通过数值插值法将得到的气动力施加在风扇叶片的结构网格上，而后通过流固耦合算法进行离心力分析，即可得到风扇叶片在离心力-气动力耦合下的力学响应。经计算，空心叶片蒙皮和加强筋在来流空气流量为14kg/s，转速为4000r/min时的应力分布如图8所示。

由图8可知，当风扇叶片转速为4000r/min时，在离心力与气动力的综合作用下，叶片的最大应力集中在吸力面中下部位处，达到了837MPa，但却不超过钛合金的屈服应力950MPa；此外，风扇叶片的最大应力分布在靠近叶根的蒙皮处，且无论是吸力面还是压力面，最大应力均分布在靠近叶根处，因此在设计时可适当增大叶根部位蒙皮和加强筋的厚度，以增加叶片承受载荷的能力。

由于在4000r/min时，风扇叶片结构的最大应力已达到了叶片材料屈服应力的88%，因此本研究主要研究风扇叶片在转速为0~4000r/min时的振动特性，进而确定合适的转速工作区间。

图8 空心风扇叶片流固耦合下的应力分布

3.2 叶片振动特性分析

通过节点特征方程解得特征根和其对应的向量矩阵，可获得叶片的固有模态。结合激振力计算式和叶片固有频率可绘制Campbell图用于叶片频率裕度的分析，从而寻找叶片合理的工作转速区间，本研究对转速为0r/min至4000r/min的叶片进行模态分析。

选取K=1，2，3，4，5代表低频激振力得出航发风扇叶片流固耦合下的Campbell图如图9所示，图中可以看到在2000r/min转速后，激振力与前三阶的固有频率中共有5个交点，与第四阶固有频率没有测试。在2000r/min以下主要是风扇叶片启动和停车阶段的运行转速，一般不会停留过长，因此本研究主要考虑2000r/min至4000r/min之间的振动特性。当风扇叶片在这5个共振点处进行长时间工作时，会发生共振现象，产生危险。为了让叶片运转时避开共振点并寻找合理的工作区间，需要进一步计算叶片的共振频率裕度。

图9 风扇叶片流固耦合下的Campbell图

风扇叶片在0~4000r/min转速下一至四阶的频率裕度图线如图10所示，图中可以看到在转速为1800~2900r/min时，一阶和二阶频率裕度小于10%;转速超过2900r/min后，一阶频率裕度大于10%。而在2900~3800r/min之间，二阶频率裕度也大于10%,并且三阶裕度不小于5%。对于风扇叶片，模态阶次越低，振动能量越大，产生的振动最危险，因此工程上最关注前三阶阶模态的频率裕度，一般只考虑前三阶的裕度，且要求第一、二阶频率裕度大于10%。由于高阶振动产生的能量较低，所以不考虑第四阶以上振动的频率裕度。因此，该风扇叶片可以迅速加速至2900r/min，快速通过低阶共振点，然后在转速为2900~3800r/min的区间内工作，这样较为安全和稳定。

图10 风扇叶片前四阶频率裕度图

3.3 叶片抗鸟撞冲击响应分析

取该风扇叶片在巡航阶段的转速3344r/min为工作转速，此时叶片的前三阶频率裕度均大于10%，设定鸟体初始速度为100m/s。根据马力等人建立的概率经验公式可以得到风扇叶片转子外围叶尖区域位置被撞击的概率，在距风扇叶片转子轴线径向距离为0.77R处的叶尖区域被鸟体撞击概率最大，造成的风险最高。

在进行鸟撞叶片显式动力学分析之前，首先对风扇叶片赋予预应力，按照3.1章节的方法获得风扇叶片的离心力和气动力作为预应力载荷。设定鸟体撞击位置为77%叶高处，针对叶片转子整体设置运行转速，同时设置了SPH鸟体粒子与风扇叶片有限元网格之间的接触方式为撞击接触，并采用了FEM-SPH耦合算法进而同时描述鸟体SPH粒子和叶片拉格朗日单元之间的载荷传递。如图11所示为鸟撞航空发动机风扇叶片有限元模型和边界条件设置。

图11 鸟撞空心风扇叶片有限元模型和边界条件

在航空发动机风扇叶片后方还存在其他转子部件后，因此在鸟撞冲击下需要主要关注受撞击叶片的轴向位移情况。在鸟体受到切割过程中，颈部最先接触到叶片，此时叶片的受力和变形较小，可以先不考虑。因此本研究选取主要切割鸟体躯干部分的5块叶片进行研究，获取叶片后缘叶尖处在鸟撞过程中的轴向位移变化如图12所示。图中可以看到5块叶片发生变形时均产生了振荡回弹现象。其中，1号叶片和5号叶片在鸟撞过程中叶片位移波动较小，这是由于1号叶片主要切割鸟体颈部，而5号叶片主要解除到散开的零碎鸟体，因此所受影响程度较小。2、3、4号叶片位移要明显大于其他叶片，其中3号叶片在切割鸟体过程中的轴向位移最大，在6.8ms时一次峰值达到了50.8mm。2号叶片的轴向位移仅次于3号叶片，但2号叶片最先切割鸟体躯干部分，因此其轴向位移一次峰值最先出现，在6.4ms时达到了45.7mm。

值得注意的是，鸟体在8.0ms时已经被切割完成，因此在8.0ms后叶片不再受到鸟体冲击载荷。而在图中可以看到受到鸟体冲击影响较大的2、3、4号叶片在8.0ms后其轴向位移开始剧烈振荡，3号叶片的轴向位移在经过减小后增大逐渐达到二次峰值，在12.1ms时达到了78.3mm，表明该叶片可能会与后面的转子叶片发生进一步碰撞。同时2号和4号叶片的轴向位移也经过了先减小后增大进而达到了二次峰值，同时二次峰值均大于一次峰值。形成这种现象主要是由于叶片经过鸟体撞击后开始变形同时叶片刚度降低，但叶片仍在旋转，因此离心力依旧存在，这使得叶片运行不稳定呈现出沿轴向往返移动的周期变形趋势。

图12 叶片后缘轴向位移随时间变化曲线

图13展示了空心风扇叶片受到鸟撞载荷后的等效应力演化过程，在5ms时，1号最先接触到鸟体，1号和5号叶片主要在前缘位置产生局部应力。在6ms时，1号和2号叶片应力水平下降，3号叶片开始出现大应力区域，最大应力达到2000MPa，随着冲击过程的持续进行，应力波也由局部扩散至整个叶片，从而导致了整个叶片应力的大幅增加。在7ms时，3号叶片的应力水平不再增加，但整体叶片开始变形且较为严重，在此期间4号叶片前缘开始出现局部应力并开始变形。在8ms时，4号叶片整体开始变形，此时鸟体已经切割完毕，4号叶片主要受到鸟体碎片的冲击，而不再切割鸟体，因此前缘部分未出现高应力区域，而是逐渐转移到叶片压力面，此时主要是鸟体碎片粒子的溅射作用所造成的。

当叶片遭受鸟体冲击后，其前缘受撞击区域所出现凹陷、卷边变形等现象将会降低叶片的气动性能和稳定性，进而影响整个发动机性能，而且撕裂和碎片脱落还可能引发二次损伤，对整机的安全性构成威胁。文献研究表明当风扇叶片受损叶片为3片时，风扇最大推力将下降20%左右，图中可以看到本研究所研究鸟撞叶片过程中至少有4片叶片受损较为严重，因此风扇推力将受到大幅度下降。

图13 叶片在鸟撞不同时刻的等效应力演化过程

空心风扇叶片受到鸟撞载荷后的等效应变演化过程如图14所示，图中可以看到5ms时，1号叶片前缘只有小区域产生应变，这是由于该叶片主要切割鸟体颈部。在6ms时2号和3号叶片前缘局部出现卷曲，此时接触鸟体的叶片前缘区域等效应变达到最大值。在7ms时2号和3号叶片前缘整体开始卷曲变形，同时4号叶片前缘也开始变形。在8ms时5号叶片受到鸟体碎片粒子的溅射，从叶片压力面开始变形，同时可以看到3号和4号叶片沿周向继续变形并发生碰撞，在碰撞接触区域产生了高应变区域。

为降低风扇叶片所受鸟撞载荷，提升叶片存活度，可以考虑针对叶片几何外形进行优化，例如增加前缘厚度提高抗鸟撞能力或调整前缘角度降低鸟撞载荷。同时也可以考虑增加外部防护装置，如防护索等结构，在鸟体接触到叶片之前将其打散。

图14 叶片在鸟撞不同时刻的等效应变演化过程

图15展示了鸟体在撞击叶片是被切割的全过程，在0-5.0ms期间主要是颈部受到叶片切割，5.0ms时鸟体颈部基本被切割完毕，6.0ms时鸟体躯干被切割了11/2。其中躯干部分主要是被三个叶片完成了切割过程，同时对应图12中轴向位移最为明显的2号、3号和4号叶片。在8.0ms时已经基本完成了鸟体的切割过程，鸟体碎片粒子开始飞溅开进而均匀接触到周围叶片，同时在自身轴向速度的作用下逐渐离开风扇叶片流道。

图15 鸟体被叶片切割的全过程 

本研究通过鸟体姿态图像构建了包含生物特征形状的真实鸟体模型，采用流固耦合方法对某型航空发动机的宽弦空心风扇叶片的振动特性与鸟撞冲击响应进行了研究。主要研究结论如下。

1）钛合金宽弦空心风扇叶片在转速为2900r/min至3800r/min时，第一、二阶的频率裕度均在10%之上，第三阶的频率裕度不低于5%，具有良好的共振裕度；在该转速范围内可以避免共振引发的叶片结构损伤，确保发动机的安全性和稳定性。

2）在灰鹤真实形状鸟体以100m/s正面撞击转速为3344r/min的风扇叶片过程中，2、3、4号叶片位移明显大于其他叶片，其中3号叶片在切割鸟体过程中的轴向位移最大，达到了78.3mm，表明该叶片可能会与后面的转子叶片发生进一步碰撞。

3）叶片应力分布和塑性变形表明，在切割鸟体过程中，应力集中首先在叶片前缘位置处出现，随后扩散至整个叶片，导致叶片整体产生显著变形；叶片切割鸟体后，由于受到鸟体碎片的冲击，应力峰值从叶片前缘部分转移至叶片压力面区域。