CN110162821B - 一种计算鸟撞高速旋转发动机叶片的方法 - Google Patents

Abstract

本发明涉及一种计算鸟撞发动机叶片响应的数值仿真方法，基于光滑粒子法(Smoothed Particle Hydrodynamics Method,SPH)和有限元法(Finite Element Method,FEM)，建立了考虑鸟体姿态角的撞击环境因素影响的数值仿真模型，在此模型下进行动力学响应的仿真。此外，本发明使得本领域内技术人员在进行鸟撞动力学响应研究的过程中能够更加简单的修改不同参数进行研究，为研究航空发动机结构的设计和优化提供参考依据。

Description

一种计算鸟撞高速旋转发动机叶片的方法

技术领域

本发明属于冲击动力学建模领域，具体涉及一种计算鸟撞发动机叶片响应的仿真方法。

背景技术

SPH方法是一种无网格方法，这种方法是模拟流体流动的一种拉格朗日型粒子方法，不像有限元方法，需要对计算域进行精确的网格划分。SPH法利用一系列带有质量、密度、速度和加速度等物理量的粒子，来求解具有各种边界条件的积分方程或偏微分方程，得到精确解。

鸟撞问题一直严重威胁航空飞行器的安全，为了降低财产损失，提高飞行器的安全性，利用数值计算和实验，能够飞机结构的抗撞性。利用数值分析方法建立研究对象的模型，通过划分网格单元，施加边界条件和载荷，然后利用求解器计算，分析响应结果。

航空发动机作为飞机的心脏，直接影响着飞机的性能和安全性，进行鸟撞击发动机叶片的响应分析，降低其引起的损坏程度，保证航空安全，有十分重要的意义。

Frederik在《Numerical and experimental investigation of the shock andsteady state pressures in the bird material during bird strike》一文中分别进行了三组不同的对比实验，使用了圆形平板，楔形和分束器，他发现明胶模型和真实鸟体，对平板和楔形物的撞击，冲击力和动量几乎相同，同时，他介绍了另一种测量动量转换的方法，进一步说明了明胶替代鸟体模型的可行性和准确性，但是他未考虑靶体的运动情况。文颖娟在《叶片鸟撞击响应分析模化技术与验证研究》一文中开展了鸟碰撞叶片平板的瞬态响应实验，测量了叶尖的位移时间曲线，并将近似弹塑性材料模型和流体动力学材料模型进行了对比，分析了每种材料参数变化对叶片响应的敏感性，并采用有限元方法比较了不同的鸟体接触模式，评价了不同的参数对叶片响应的结果，但是在她的研究中并没有涉及到考虑曲率的叶片。

发明内容

本发明的目的在于提供了一种计算鸟撞发动机叶片响应的仿真方法，以冲击动力学为理论基础，目的在于分析考虑鸟体姿态角的鸟撞发动机叶片数值仿真模型的动力学响应。

实现本发明目的的技术解决方案为：一种计算鸟撞发动机叶片响应的仿真方法，方法步骤如下：

步骤1、设定鸟体的材料参数和速度参数，发动机叶片的速度参数、材料参数和状态方程参数，建立鸟体数值仿真模型，转入步骤2；

步骤2、基于FEM对鸟体数值仿真模型进行离散，利用SPH和FEM转换的方法，将鸟体数值仿真模型的FEM单元转换成SPH粒子，转入步骤3；

步骤3、计算SPH粒子和FEM单元之间的接触力，其中SPH粒子采用跳蛙格式求解纳维斯托克方程，有限元采用中心差分法求解显式动力学方程，转入步骤4；

步骤4、根据求解后的显式动力学方程，得到鸟体数值仿真模型的能量变化，转入步骤5；

步骤5、输出鸟体数值仿真模型能量随时间变化图。

本发明与现有技术相比，其显著优点：

(1)、本发明考虑了鸟体模型姿态角的影响，使计算更加贴近实际情况。

(2)、在国际先进计算方法SPH和FEM法的基础上建立了数值计算模型，与国际接轨。

(3)、本发明计算出鸟体数值模型动能的变化，为工程应用提供了一定的指导作用。

附图说明

图1为本发明一种计算鸟撞发动机叶片响应的仿真方法的流程图。

图2为本发明建立的鸟体数值仿真模型图。

图3为实施例中计算出鸟体的动能随时间变化曲线图。

图4实施例中不同姿态角鸟体的动能随时间变化曲线图。

具体实施方式

结合图1、图2，一种计算鸟撞发动机叶片响应的仿真方法，包括以下步骤：

步骤1、设定鸟体的材料参数和速度参数，发动机叶片的速度参数、材料参数和状态方程参数，建立鸟体数值仿真模型，转入步骤2；

坐标系中的x₀，y₀，z₀分别是鸟体数值仿真模型的位置坐标，鸟体数值仿真模型的密度ρ，弹性模量为E，泊松比为γ，ω_e为叶片的旋转角速度，C为材料声速，S₁、S₂、S₃为体积应力应变曲线的倾斜系数，γ₀为Mie-Gruneisen系数， a为体积修正因子,E为弹性模量。

步骤2、基于FEM对鸟体数值仿真模型进行离散，即将鸟体数值仿真模型的结构离散化，画出有限元网格图，对节点、FEM单元分别编号，定出节点的x、 y坐标即FEM单元i、j、m三节点对应的整体号，由节点坐标计算常数b_i、c_i及 FEM单元的面积，再计算单元刚度矩阵[k]^e；再利用SPH和FEM转换的方法，将鸟体数值仿真模型的FEM单元转换成SPH粒子，转入步骤3；

有限元法的形函数为：

其中，Q_i、Q_j、Q_j分别对应不同的三个节点，I为单位矩阵。

体积力为：

其中，{G}为单位面积上的体积力，tdxdy为微分体积。

节点力为：

F^e＝[k]^eδ^e (15)

其中，[k]^e为单元刚度矩阵，δ^e为节点位移列阵。

FEM单元转换成SPH粒子，利用以下三种方法实现：

1、FEM单元转换成Ⅰ型SPH粒子；

2、FEM单元转换成Ⅱ型SPH粒子；

3、FEM单元转换成Ⅰ型SPH粒子,再转换成Ⅱ型SPH粒子；

步骤3、计算SPH粒子和FEM单元之间的接触力，其中SPH粒子采用跳蛙格式求解纳维斯托克方程，有限元采用中心差分法求解显式动力学方程，转入步骤4；

对于两个物体因接触变形，它们的运动方程以及边界条件可表示成：

运动方程：

▽σ+ρb＝ρa在(Ω＝Ω_A∪Ω_B域内) (16)

位移边界条件：

位移边界条件：

(在Γ_u边界上) (18)

其中，a表示加速度，n为方向矢量，

表示法向矢量,u为位移，Γ为计算域。两个物体之间发生接触的条件可表示成：

其中，Γ_C＝Γ_A∩Γ_B，Γ_t∪Γ_t0∪Γ_u，Γ_t∩Γ_t0＝0，Γ_t0∩Γ_u，Γ_t∩Γ_u＝0。

将失效单元节点物理量赋给新的SPH粒子，具体包括：

位移：

x_p＝x_n (20)

质量：

速度：

v_p＝v_n (22)

应力：

声速：

c_p＝c_n (24)

其中，下角标p和n表示SPH粒子和有限元单元的节点，x表示位移，m为质量，N_e为与节点相关的有限单元总数，N_n为一个单元的节点个数，ρ_ei为单元材料的密度，V_ei为单元的体积，σ_ei为单元i的应力张量，w_gj和σ_gj表示单元i 的内部高斯积分点g处的加权系数和应力张量，N_g为一个单元高斯积分点数。

新的SPH粒子的光滑长度、应变率张量以及扭转率张量可以通过以下公式计算：

其中，ρ₀是初始的密度，r₀为初始单元尺寸，α为定义的光滑长度与单元尺寸的比值，d为单元维数。

由公式(10)，利用SPH的邻近搜索方法，计算进入粒子支持域范围内未失效单元节点产生的接触力。

计算粒子i的速度时，所有在其两倍光滑长度范围内即在支持域内的粒子，如物体A内n₁、n₂等，以及物体B内的n₁、n₂和n₃粒子等，SPH方法动量方程式：

其中，m_Bj为在物体B中的粒子的质量，ρ_Bj和σ_Bj表示物体B中粒子密度和应力，N_B表示在物体B中处于粒子i支持域内的粒子总数。

将鸟体模型与有限元模型进行接触计算，即光滑粒子法与有限元法以及流体与固体的耦合计算。未被转换的有限元单元网格的节点作为背景粒子，用于被搜索，当SPH粒子进入作为背景粒子的有限元影响域内，产生接触力，有限元节点进入SPH粒子的支持域内也会产生接触力。

步骤4、根据求解后的显式动力学方程，得到鸟体数值仿真模型的能量变化，转入步骤5；

其中局部凹陷的薄膜变形能可表示为：

弯曲变形能为：

N_x'和N_z'为叶片中面内力。

为x'y'平面内的曲率，

为 y'z'平面内的曲率，M_x'为x'方向的弯矩，M_z'为z'方向的弯矩。

叶片塑性形变所消耗的能量为ΔW_shear，可表示成：

其中2s为剪切过渡区的宽度，ε_eff表示极限拉伸的应变，γ₀表示最大剪切应变，

σ_y为材料的动态屈服应力，H为叶片厚度，Δz'为在t时刻撕裂材料的长度，μ为剪切位移。

弯曲变形所消耗的能量ΔW_bend为：

其中，R是圆柱形外物半径，κ为叶片上弯曲部分的曲率半径，与弹道角θ有关。

在碰撞中的能量方程可表示为：

E_p＝E₀-U_P-∑W_i-∑W_p (33)

其中E₀为碰撞前的外物动能，E_p为碰撞后的外物剩余动能，U_P为外物变形消耗的能量，∑W_i为叶片在碰撞中所耗能量，包括变形能或生成裂纹的能量，∑W_p为碰撞产生的热以及冲击波传播消耗的能量和，忽略不计。

步骤5、输出鸟体数值仿真模型能量随时间变化图。

利用计算软件求解方程后，得到的动能变化数据，利用ORIGIN软件绘制动能变化曲线。

以上所述为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，都可以轻易的变化或者替换，这些都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

实施例1

下面结合附图以及具体实施例对本发明进行进一步介绍，显然所描述的实施例仅仅是本发明的某一具体实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例公开了一种计算鸟撞发动机叶片响应的仿真方法，具体方法如下：

步骤1、根据公式(1)建立鸟体数值仿真模型，鸟体速度116m/s，密度 937kg/m³，泊松比3.3，叶片的材料参数和状态方程参数选取如下表所示：

表1本实施例采用叶片材料模型参数

表2本实施例采用叶片的状态方程参数

步骤2、基于FEM对鸟体数值仿真模型进行离散，利用SPH和FEM转换的方法，将鸟体数值仿真模型的FEM单元转换成SPH粒子，转入步骤3；

步骤3、计算SPH粒子和FEM单元之间的接触力，其中SPH粒子采用跳蛙格式求解纳维斯托克方程，有限元采用中心差分法求解显式动力学方程，转入步骤4；

步骤4、根据求解后的显式动力学方程，得到鸟体数值仿真模型的能量变化，转入步骤5；

步骤5、通过ORIGIN输出鸟体数值仿真模型能量随时间变化图如图3，改变初始坐标(x₀，y₀，z₀)，输出鸟体数值仿真模型能量随时间变化图如图4所示，由图3、图4可知，随着时间的增加，鸟体动能因碰撞引起动能值发生改变，逐渐减少，这种现象是由碰撞导致速度减小造成的。改变姿态角，动能衰减趋势和范围不同。

Claims (1)

1.一种计算鸟撞发动机叶片响应的仿真方法，其特征在于，包括以下步骤：

步骤1、设定鸟体的材料参数和速度参数，发动机叶片的速度参数、材料参数和状态方程参数，建立鸟体数值仿真模型，转入步骤2；

步骤2、基于FEM对鸟体数值仿真模型进行离散，利用SPH和FEM转换的方法，将鸟体数值仿真模型的FEM单元转换成SPH粒子，转入步骤3；

步骤3、计算SPH粒子和FEM单元之间的接触力，其中SPH粒子采用跳蛙格式求解纳维斯托克方程，有限元采用中心差分法求解显式动力学方程，转入步骤4；

步骤4、根据求解后的显式动力学方程，得到鸟体数值仿真模型的能量变化，转入步骤5；

步骤5、输出鸟体数值仿真模型能量随时间变化图；

所述步骤1中的鸟体数值仿真模型建立方程为：

其中，x₀，y₀，z₀分别是鸟体数值仿真模型的坐标；

鸟体的材料参数分别为鸟体数值仿真模型的密度ρ，弹性模量为E，泊松比为γ，碰撞点处的速度V_e表示为：

V_e＝V_br+V_bz＝ω_e×r_e+V_bz＝[r_eyω_e -r_exω_e v_bz] (2)

其中，V_br是碰撞点在叶片上的切向速度，ω_e为叶片的旋转角速度，r_e是碰撞点距离叶片的旋转轴z的旋转半径，V_bz是叶片相对z轴的平移速度，外物对于碰撞点的相对速度V_r表示成：

V_r＝V_i-V_e＝[V_rx V_ry V_rz]＝[v_ix-r_eyω_e v_iy-r_exω_e v_iz-v_bz] (3)

鸟体数值仿真模型材料本构方程式为：

其中，σ_ij为应力，P为压力，δ_ij为克罗内克函数，

为交错符号，C为体积应力应变曲线的间断系数，即材料声速取C＝1483m/s，S₁、S₂、S₃为体积应力应变曲线的倾斜系数，S₁＝1.75，S₂＝S₃＝0，γ₀为Mie-Gruneisen系数，γ₀＝0，a为体积修正因子，变量μ＝ρ/ρ₀-1，E为弹性模量，ρ₀为初始密度；

所述步骤2中使用有限元方法将鸟体数值仿真模型离散，即将鸟体数值仿真模型的结构离散化，画出有限元网格图，对节点、FEM单元分别编号，定出节点的x、y坐标即FEM单元i、j、m三节点对应的整体号，由节点坐标计算常数b_i、c_i及FEM单元的面积，再计算单元刚度矩阵[k]^e；再把FEM单元转换成SPH粒子，利用以下三种方法实现：

1、FEM单元转换成Ⅰ型SPH粒子；

2、FEM单元转换成Ⅱ型SPH粒子；

3、FEM单元转换成Ⅰ型SPH粒子，再转换成Ⅱ型SPH粒子；

不同的类型需要根据单元失效情况来判断，若失效单元的周围单元并没有全部失效，则在失效与并未失效的单元交接处，生成粒子，定义为Ⅰ型SPH粒子；若失效单元的周围单元也全部失效，则在此失效单元处，生成Ⅱ型SPH粒子；

所述步骤3中，粒子之间的接触力为：

其中，N为粒子i在SPH方法中的形函数，b_c为接触部分的体力；

形函数N为：

其中，j是粒子i支持域内的粒子，m为粒子质量，N_p是粒子i支持域内的总粒子数量，W_j(x_i)是粒子i的权函数；ρ_j是SPH粒子的密度；

至于b_c，利用接触势函数，将其定义为：

其中，φ为接触势函数，NCONT是对粒子i有影响的粒子总数量，包括不同物体的粒子，有限元单元网格节点也计入粒子，K和n是自定义的参数，K的取值与材料性质和接触速度有关，r_ij是粒子之间的距离，Δp_avg则是所有有影响的粒子的光滑长度的平均值；

带入公式(6)中，得到：

采用SPH粒子近似得到：

SPH法中的光滑核函数的梯度决定接触力的方向；

在SPH粒子上施加接触力，需要对SPH动量方程修正：

在有限元节点上施加接触力，需要把有限元动力学方程中的节点载荷替换为外力f_c(x_i)：

其中，M为鸟体数值仿真模型的质量，C为鸟体数值仿真模型的阻尼，Q为鸟体数值仿真模型的刚度矩阵。

Images