直升机在机动状态下的桨叶弹击数值模拟

林长亮，朱跃法,胡文刚

(航空工业哈尔滨飞机工业集团有限责任公司，哈尔滨 150066)

0 引 言

为了提高军用直升机战场生存力，国外对武装直升机旋翼桨叶的抗弹击提出了明确要求。例如，要求黑鹰直升机桨叶在被速度770 m/s、口径7.62 mm穿甲弹距离100 m击中时，不影响作战任务的概率不低于99%，在被速度500 m/s、口径12.7 mm穿甲弹距离800 m击中时，能安全飞行30 min的概率不低于97%，被一发口径23 mm炮弹击中不得损毁[1]；我国《军用直升机生存力要求GJB 3696-99》，明确规定了直升机旋翼桨叶必须具有一定的抗弹击能力[2]。

国内外对直升机旋翼桨叶抗弹击性能进行相关研究，A.R.Kiwan等[3]建立了冲击波对旋翼的毁伤模型，给出了旋翼毁伤概率曲线；B.Rasuo等[4]研究了7.9 mm弹丸对直升机主翼和尾翼的损伤效应，评估了损伤概率；M.Giglio等[5]采用7.62 mm弹丸，对直升机的尾桨轴进行了弹击试验，并使用扭转疲劳试验仪测试了弹击后尾桨轴的疲劳寿命，然后采用ABAQUS进行了裂纹扩展模拟对比了试验结果。上述研究结果表明：在裂纹扩展模拟中，精确描述弹道冲击现象，给出可靠的数值模型能够估算出较为精确的结果。B.Rasuo[6]使用7.9 mm孔径弹丸对静止状态的重型运输直升机尾桨叶根部位置进行了弹击，然后对损伤桨叶和非损伤桨叶进行了模态和疲劳试验对比。测试结果表明：损伤桨叶在最严重载荷工况下仍然能够飞行65小时，其振动测试结果与未损伤桨叶基本相当，充分的证明了复合材料结构优秀的抗弹击性能和高生存力。W.K.Chiu等[7]评估了一种直升机复合材料主旋翼叶片的战斗损伤修复能力，采用有限元分析方法对原始模型、受损模型和修复模型进行了比较，对所提出的修复概念的有效性进行了评估。

刘会平[8]建立了直升机超轻型陶瓷/复合材料装甲的有限元模型，分析了由不同材料和铺层组成的复合装甲的防弹性能，并对其进行了优化设计。研究结果表明：璃纤维和芳纶交替铺设具有较好的防弹性能。刘刚等[9]建立了旋翼等效剖面翼型模型，计算了碎片和冲击波对剖面翼型的联合毁伤，得到了冲击波和碎片联合作用下旋翼毁伤的临界判据。孙中涛等[10]建立了桨叶典型结构段的三维模型，在弹孔处采用节点移除法模拟了桨叶弹击损伤状况。郭超[11]研究了碎片式高射炮对武装直升机的毁伤效能，通过试验验证了碎片计算模型的结果。研究结果表明：炮弹从侧面对直升机造成毁伤的概率高于底部，高炮杀伤弹和AHEAD的毁伤效能影响因素不同。宫小泽等[12]以某武装直升机为研究对象，分析了AHEAD弹对目标的毁伤效能及毁伤目标所需的平均命中弹数。研究结果表明：毁伤目标所需平均命中弹数受弹丸速度、脱靶距离、弹丸俯仰角等因素影响。

综上所述，目前对于直升机弹击的研究大都使用概率计算公式或数值方法对直升机机体、直升机元件或旋翼典型剖面段进行损伤概率评估或数值模拟，尚未研究直升机在机动状态下，全尺寸桨叶结构受到载荷作用时，不同速度弹头对桨叶造成冲击损伤的数值模拟方法。本文采用柔性多体动力学、瞬态响应、自适应接触算法的联合仿真对某型直升机桨叶的弹击进行研究。首先介绍桨叶机动状态下的载荷计算方法，计算某型直升机旋翼桨叶的载荷与离心力；然后建立桨叶的全尺寸有限元模型，将桨叶载荷施加到有限元模型上，进行瞬态分析；阐述自适应接触算法原理，采用文献[13]中的试验值对该方法的正确性进行验证；最后，分析不同速度的圆柱弹头对桨叶造成的损伤。

1 桨叶机动状态的数值模拟方法

1.1 结构介绍

桨叶剖面结构和展向铺层布置如图1所示[14]，桨叶展长5.275 m，主要结构包括:Kevlar预浸带大梁、前后Z型梁、上下翼面蒙皮、泡沫和不锈钢包铁。材料属性如表1所示，其中E11为横向弹性模量，E22为纵向弹性模量。

(a) 桨叶结构展向铺层示意图

(b) 桨叶结构剖面示意图图1 旋翼桨叶结构Fig.1 Rotor blade structure 表1 材料属性Table 1 Material property

材料名E11/MPaE22/MPa泊松比破坏值/MPa不锈钢210 000-0.331 080玻璃布18 00018 0000.25 180碳布23 00023 0000.25 345预浸带34 22034 2200.31 800

1.2 剖面特性计算

剖面特性是旋翼动力学分析的基础，目前对于直升机桨叶载荷以及动力学特性的计算，都是通过求解桨叶每个剖面位置的挥舞/摆振刚度、惯性矩、剪心等剖面特性数据，构建复合材料结构桨叶的等效模型。

本文使用Fortran语言程序，基于Euler-Bernoulli梁理论[15]和Vlasov薄壁梁理论[16]，对桨叶不同站位处的剖面进行有限元网格划分，如图2～图3所示，再使用VABS软件，对桨叶整体结构的剖面特性插值计算[17]。

图2 桨叶根部的有限元网格Fig.2 Finite element mesh of blade root

图3 桨叶段的有限元网格Fig.3 Finite element mesh of blade segment

1.3 桨叶载荷计算

直升机旋翼桨叶载荷包括静载荷和动载荷，计算状态包括超转、超扭、水平飞行、水平弯转、螺旋弯转、自由下滑等[18]。

本文选取了直升机前飞状态，旋翼转速360 rpm。使用CAMRIDII[19-20]软件，采用柔性多体动力学理论，建立了某型直升机桨叶的弹性结构模型，计算了桨叶结构的静载荷、动载荷，以及高速旋转产生的离心力，载荷计算流程如图4所示。

图4 桨叶载荷计算流程Fig.4 Blade load calculation process

动载荷计算采用包含机身结构的全机配平模型，取6阶模态与谐波，自由尾迹选取Johnson模型；静载荷采用孤立的旋翼桨叶模型，取6阶模态与0阶谐波，自由尾迹选取Johnson模型[21]。

桨叶弹性结构各处的载荷计算公式如式(1)～式(3)所示，离心力引起的应变为

(1)

挥舞力矩引起的应变为

(2)

摆振力矩引起的应变为

(3)

式中：ES,EIB,EIT为拉伸刚度，挥舞刚度，摆振刚度；FCN为剖面处离心力；MB为剖面处的挥舞力矩；MT为剖面处的摆振力矩；YN,ZN为剖面张力中心位置坐标；X,Z为剖面计算点位置坐标。

将上面计算的三部分应变叠加，并乘以各铺层材料的弹性模量得到各剖面站位的载荷。计算结果如图5～图6所示，可以看出：在前飞状态下，离心力对桨叶结构产生的影响较大，挥舞和摆振产生的弯矩较小。

图5 桨叶离心力Fig.5 Centrifugal force of blade

图6 桨叶挥舞/摆振力矩Fig.6 Blade swing/shimmy moment

2 旋翼桨叶弹击有限元建模

使用MSC.PATRAN建立了桨叶的全尺寸有限元模型，如图7所示。上下蒙皮和前后Z型梁结构用二维SHELL元模拟，Kevlar预浸带大梁用六面体体元模拟，其中二维壳单元数量为3万个，六面体单元数量为1.2万个。由于壳单元节点具有六个自由度，体单元节点有3个自由度，为了让壳单元的扭转载荷能够传递到体单元上，使用MPC对单元对应节点一一连接，以保证载荷传递的准确性。弹体直径8 mm，高6 mm。

图7 旋翼桨叶剖面Fig.7 Rotor blade profile

3 桨叶预应力加载

将计算的载荷施加到桨叶有限元模型上，使用RBE3完成载荷加载，加载位置选取为桨叶轴线与剖面的交点，如图8所示。

图8 桨叶加载Fig.8 Blade loading

在桨叶模型上加载主要是需要得到桨叶在载荷作用下的位移响应状态，故使用NASTRAN的SOL106经典非线性静力求解器完成计算，通过迭代计算得到目标加载载荷时，桨叶结构的位移变形。然后生成后缀为SOL的预应力文件，将此文件读入DYTRAN中，完成桨叶预应力加载模拟。

4 桨叶弹击损伤数值模拟

4.1 材料失效判据

由于旋翼桨叶结构的有限元模型由二维单元和三维单元构成，根据DYTRAN的材料库特点[22-23]，使用两种复合材料本构和一种金属失效本构作为失效判据。

针对上下蒙皮结构和前后Z型梁等二维复合材料单元，采用Hashin失效准则的本构[12]。Kevlar预浸大梁使用Dytran中的DMATOR正交各向异性三维本构。不锈钢前缘包铁使用DMATEP各向同性弹塑性本构[24]。

Hashin失效准则[23]表述如下:

纤维拉伸失效

(4)

纤维压缩失效

(5)

基体拉伸失效

(6)

基体压缩失效

(7)

式中:XT和YT为纤维轴向和横向的拉伸强度；XC和YC为纤维轴向和横向的压缩强度；SL为纤维轴向剪切强度。

正交各向异性三维本构表述如下:

σ=Cε

(8)

式中:σ为应力矩阵;C为刚度矩阵;ε为应变矩阵。

C=TTCLT,其中T为材料坐标系(a,b,c)与基本坐标系之间的变换矩阵，CL为定义在材料坐标系中的本构矩阵，如式(9)所示:

(9)

4.2 自适应接触算法

自适应接触法的基本原理是:在每一时间步，首先检查各从节点是否穿透主面，如果发生穿透，接触面之间会产生相互作用力，使它们相互分开，当组成接触面的单元因承受过大力而发生破坏时，这些单元将被从接触面中删除，从而改变了接触面的构型[23]。

本文将子弹与局部冲击部位作为接触区定义接触，子弹使用刚体元模拟，桨叶根部约束XYZ方向。

旋翼桨叶在360 rpm转速下，弹击位置线速度为180 m/s，弹头速度为300～600 m/s。假设旋翼处于悬停状态，弹头相对速度为480～780 m/s。

4.3 算例验证

为了验证接触算法的可靠性，本文引用文献[13]的试验结果，模拟了平头弹冲击5 mm厚(0°/90°)5正交Kevlar层合板，弹头尺寸参阅文献[13]，层合板为直径160 mm的圆板，模型边界条件示意图如图9所示。经计算，计算值与文献中的试验值误差在10%以内，验证了算法的精度，结果对比如图10所示。

图9 模型边界条件示意图Fig.9 Schematic diagram of model boundary conditions

图10 计算值与试验值对比Fig.10 The calculated value compared with the experimental value

计算结果的试验值与计算值对比如表2所示。

4.4 桨叶弹击数值模拟流程

本文介绍的弹击数值模拟流程如图11所示。在准备好桨叶需要的输入数据后，使用自编程序按不同剖面站位划分剖面的有限元网格，再使用VABS软件完成剖面特性计算；将剖面特性数据输入到CAMRIDII中建立弹性桨叶结构的载荷分析模型，完成桨叶载荷计算；在此基础上，使用NASTRAN软件对桨叶的全尺寸有限元模型进行瞬态分析，得到桨叶施加载荷的预应力状态；使用DYTRAN软件，基于自适应接触法，模拟桨叶弹击模拟损伤。

图11 桨叶弹击流程Fig.11 Blade bullet penetration flow

5 桨叶弹击计算结果

5.1 桨叶的结构损伤部位

经计算，弹击后主要失效部位为上下蒙皮与大梁的过渡区域和大梁、前缘蒙皮，如图12所示。桨叶蒙皮的主要失效形式为层间，前缘弹击点处有一定程度的侵彻损伤。

图12 弹击主要失效部位剖面图Fig.12 Section diagram of the main failure site of missile strike

弹头在击中桨叶前缘后，大梁和前缘蒙皮首先吸收了弹头的初始动能，当弹头侵彻大梁到一定深度时，对于上下蒙皮与大梁连接的过渡区域，大梁瞬间产生的挤压载荷对薄弱的蒙皮过渡区造成大面积损伤。这是由于大梁刚度较大，而上下蒙皮刚度较弱，弹头在击中高刚度大梁产生的冲击载荷对于刚度较弱的过渡区域会产生明显的应力集中载荷，从而造成损伤。因此在桨叶抗弹击设计中，需要考虑到结构之间的刚度匹配问题，抗弹击部位的刚度与周边结构要基本一致或差距较小，才能够保证吸收的冲击能量在结构中传递时不会出现明显的应力集中现象，避免周边结构出现损伤。

5.2 弹头剩余速度

弹头冲击桨叶后的剩余速度如图13所示。

图13 弹头剩余速度Fig.13 Warhead residual velocity

从图13可以看出：当子弹出现反弹时，速度为负；当子弹发生侵彻现象时，速度为正；当子弹动能被桨叶完全吸收时，速度为0。从计算结果可以看出，弹头在击中桨叶后基本出现反弹，剩余速度与初始速度单调线性递增，与文献[13]的计算结果趋势基本一致，表明本文计算模型的结果合理。

5.3 桨叶根部的应变

弹头初始动能与桨叶结构的吸能对比如图14所示，不同冲击速度下的旋翼桨叶根部应变值如图15所示。由于弹头与桨叶的距离是固定的，所以不同速度的弹头撞击到桨叶的时间不同。在480～780 m/s的冲击速度下，桨叶根部应变峰值响应时间基本在2.8 ms左右。在弹头动能基本被桨叶完全吸收的情况下，应变峰值没有随弹击速度线性增大。

图14 弹头初始动能与桨叶结构的吸能对比Fig.14 The comparison between the initial kinetic energy of the warhead and the energy absorption of the blade structure

图15 不同冲击速度下的旋翼桨叶根部应变值Fig.15 The strain value of rotor blade root at different impact velocities

5.4 桨叶结构损伤程度对比

为了对比桨叶的损伤程度，对三维体元，采用单元体积作为对比参数；对于二维壳元，采用单元面积乘层合板厚度作为对比参数，不同冲击速度下的结构损伤程度如图16所示。

图16 不同冲击速度下的结构损伤程度Fig.16 Structural damage degree at different impact velocities

从图16可以看出：在大多数情况下，桨叶大梁的损伤度高，表明大梁是抗弹击的主要吸能结构，但桨叶结构的损伤程度与弹击速度是非线性的关系；当蒙皮的损伤程度高于大梁时，桨叶根部会产生较大的应变值，说明蒙皮损伤对于桨叶结构安全的危害性要高于大梁的损伤。

6 结 论

(1) 在高速冲击中，桨叶结构的抗弹性能、根部的应变峰值、峰值响应时间以及结构的损伤程度与弹击速度是非线性关系。

(2) 除了弹击点，高冲击应力会对冲击点周边结构造成损伤，冲击应力波在桨叶结构上的传递存在迟滞效应，大梁是抗弹击的主要吸能结构。

(3) 在桨叶结构的抗弹击设计中，需要保证各结构部件的刚度分布合理，以确保在弹击中桨叶结构不会出现明显的应力集中现象。

(4) 本文在桨叶弹击的数值模拟中，预先计算了桨叶在机动状态下的拉力、挥舞力矩、摆振力矩，并将其施加到全尺寸有限元模型上。通过这种对桨叶施加预应力的方法，能够更准确的模拟桨叶在实际工作状态中受到弹击时的应力分布状态和破坏形式。