基于abaqus的SPH算法鸟撞分析研究

龙思海，滕春明，张海东

(航空工业洪都，江西 南昌，330024)

0 引言

随着飞机低空高速飞行任务的增加，鸟撞飞机的概率也随之增加。鸟撞飞机带来的危害已越来越为人们所认识。由于工程试验法试验周期长、次数多，导致费用很高，例如在某型飞机风挡的研制过程中，全尺寸鸟撞试验就进行了多次，制造了大量的试验件，试验费用极高。因此，鸟撞动响应仿真分析成为重点的研究方向，型号设计前期准确的鸟撞动响应仿真分析能够节省大量的试验费用。

张志林[1]等着重考虑了应变率对透明件材料性能的影响和几何非线性对刚度矩阵的影响，结果显示：应变率对位移、应变影响较大，考虑应变率相关性分析所得结果比不考虑应变率相关性分析结果更接近试验结果；几何非线性分析所获得的风挡最大法向位移比线性分析得到的值大，几何非线性对飞机风挡鸟撞动响应分析结果的影响不可忽略，并以此提出了鸟撞击载荷柔性靶理论。

王昌银[2]等基于ANASYS/LS-DYNA[3-8]采用Lagrange网格和SPH模拟鸟体进行风挡抗鸟撞性能对比分析，鸟体采用可压碎泡沫本构模型。计算表明：SPH分析得出的结果与传统的Lagrange网格方法得出的结果以及试验实测结果基本一致。

本文基于abaqus[9-11]的SPH算法进行鸟撞动响应分析，分别采用壳单元和体单元模拟风挡，对比分析试验实测结果与计算结果之间的差异，表明abaqus的SPH算法在鸟撞分析中都是可行的，根据迭代计算得到的鸟体及玻璃的本构可用于后续型号研制分析。

1 风挡结构简介

某型飞机风挡为圆弧整体风挡，如图1所示。风挡玻璃采用航空YB-DM-3有机玻璃，骨架材料为铝合金，骨架通过螺栓与机身连接。

2 SPH方法

SPH(Smoothed Particle Hydrodynamics)方法的全称为“光滑粒子流体动力学”[12]，是Lucy和Gingold等人于1977年提出的一种无网格化拉格朗日算法[13]。SPH方法的基本思想是把材料看成一组流体粒子的集合，每个粒子具有自己的质量、速度和能量特征，并描述为一个和流体性质相关的插值点，用规则的内插函数计算所有粒子的场函数，近似描述整个问题的场分布。

图1 风挡结构示意图

考虑到鸟撞风挡动态响应的变形特点，本文使用SPH方法和有限元方法，鸟体的模拟采用SPH方法，用流动的粒子描述鸟体的大变形、破碎及飞散，其它区域使用有限元 Lagrange方法。SPH粒子与Lagrange单元的接触采用基于罚函数的接触方法。

3 基于abaqus的SPH算法仿真分析

本文采用abaqus软件建立风挡抗鸟撞有限元模型，通过修改关键字定义玻璃的破坏模式，然后通过abaqus计算求解。

3.1 鸟体模型

对于鸟体模型的形状，一般有圆柱形鸟体和两端半圆球的圆柱形鸟体。将鸟体简化为长径比为2：1的圆柱体。本文中鸟的密度为970kg/m3，鸟的重量为1.8kg，根据计算得到鸟体模型长度为212mm、直径为106mm。鸟体采用SPH粒子模拟。鸟体模型如图2所示，鸟体材料参数见表1。

图2 鸟体模型

表1 鸟体材料参数

3.2 风挡模型

风挡玻璃材料为3号定向有机玻璃，风挡玻璃单元沿厚度方向不少于3层，以模拟风挡玻璃内外表面受撞击过程中的拉压应力。鸟撞过程是一个高速的变化过程，一般在几毫秒内完成，在瞬态撞击过程中，在撞击区附近产生剧烈的应力、应变的变化，而弧框和边框连接的区域还无法响应其剧烈的应力变化，因此，约束风挡玻璃周边三个方向平动自由度作为边界条件，以模拟支持。玻璃采用弹塑性材料模拟，风挡玻璃材料参数见表2，风挡的壳单元 （S4R）和体单元（C3D8R）模型如图3所示。

3.3 分析结果

表2 风挡玻璃材料参数

图3 风挡有限元模型

有限元模型中的接触算法采用通用接触算法，鸟撞试验结果表明，该风挡可承受1.8kg鸟体正面撞击不发生破坏的临界速度为450km/h。因此，在仿真分析时，计算分析了鸟体以450km/h正面撞击的情况，取风挡对称线中点为被撞击中心点。分析结果中取被撞击中心点附近的单元应力或者节点位移进行对比分析，计算得出的被撞击点处应力、位移、应变-时间曲线以及动能-时间曲线见图4。

图4 两种风挡单元模型下应力、位移、应变、动能-时间曲线

3.4 分析与试验结果对比

图4（a）给出了两种风挡模型分别采用壳单元和体单元时的被撞击点处应力数值计算结果，比较图中3条曲线，可以看出，两种模型的计算结果与试验结果趋势一致，体单元模型得到的应力峰值比壳单元应力更大，体单元模型应力计算结果与试验测量结果吻合的更好，两种模型的应力峰值比试验结果出现的更早。

图4（b）给出了被撞击点处位移数值计算结果，并结合试验结果进行对比。可以看出三条数值曲线的变化趋势比较一致，壳单元模型计算结果更接近试验结果。图4（c）给出了被撞击点处应变数值计算结果，可以看出两种计算结果趋势一致，壳单元应变计算结果更大，体单元计算结果小。

两种模型中鸟体的速度、质量和密度一致，但由于风挡由不同种类单元建模，导致在撞击过程中鸟体损失的能量不同，如图4（d）和4（e）所示的撞击过程中动能和塑性耗散随时间的变化曲线，撞击瞬间两种模型的动能均为14091J，撞击过程中壳单元和体单元吸收的能量分别为：2191J和2681J。可见壳单元模拟的风挡模型刚度更弱，耗散了更多的能量。

两种风挡模型计算得到的应力、位移结果与试验结果对比见表3，由表3可以看出，两种模型仿真分析结果与试验测试结果吻合较好，从应力和位移峰值结果来看，壳单元模型的计算结果与试验测量结果更加接近，表明在鸟撞分析中，使用壳单元、体单元模拟风挡进行鸟撞分析都是可行的。

表3 计算结果对比

4 结论

通过鸟撞性能仿真对比分析表明：

1）基于abaqus的SPH算法进行鸟撞动响应分析是可行的；

2）仿真分析中采用壳单元和体单元模拟风挡得到的分析结果与试验结果趋势是一致的；

3）壳单元模拟风挡得到的计算结果更加接近试验测量结果且单元规模较小，在分析过程中可以有效减少求解时间；

4）通过试验结果迭代分析得到的玻璃以及鸟体的本构模型可以用于后续型号设计。