CN105787199A - 一种叶轮机转子叶片从热态到冷态的叶型转换方法 - Google Patents

Abstract

本发明涉及一种叶轮机转子的叶型转换方法，特别涉及一种叶轮机转子叶片从热态到冷态的叶型转换方法。所述的方法在已知叶轮机转子叶片热态叶型的前提下，基于有限元逆向迭代的叶型转换方法，获取叶片的冷态叶型。所述的方法提高了计算精度，并且节省了人力，提高了效率。

Description

一种叶轮机转子叶片从热态到冷态的叶型转换方法

技术领域

本发明涉及一种叶轮机转子的叶型转换方法，特别涉及一种叶轮机转子叶片从热态到冷态的叶型转换方法。

背景技术

当前叶轮机转子叶片设计中，经常需要在热态和冷态两种叶型之间进行转换。所谓冷态叶型，即工程图上给出的机械加工叶型，提供给叶片制造加工方和结构强度振动分析人员；而热态叶型指的是气动设计人员在压气机设计阶段给出的工作状态叶型，即高速旋转的柔性叶片在离心载荷、表面气动载荷和温度载荷的多重作用下形成的变形后的叶型。由于叶轮机转子叶片的冷态和热态叶型之间差别较大，存在弹性恢复角、径向拉伸、轴向偏移和周向挠度，因此作为叶轮机完整设计链中必不可少的一环，热态和冷态叶型之间的转换精度尤其需要得到设计者的重视。当前主流方法是采用近似线性化方法处理叶型从热态到冷态的转换问题，这种方法将叶片冷热态之间的变形视为完全线性过程，通过对热态叶型实施强度分析获得近似的变形，再将近似得到的变形从热态叶型中减去来获得冷态叶型。现有的方法至少存在以下三个缺点：

首先，现有方法无法考虑叶片的非线性作用，这对于过去刚性较强的压气机叶片是没有问题的，但对于新型的高负荷高性能叶片而言，其气动载荷更高而刚性更差，不能忽视由几何非线性导致的叶片大变形。因此，压气机尤其是低压压气机叶片的变形应视为非线性，不适宜进行简单的线性化处理。

其次，现有方法通过试探获得了“伪”冷态叶片，使用“伪”冷态叶片变形而获得的最终热态叶片方案与初始的气动方案存在一定差别，导致最终方案无法完全忠实于设计者的原始设计思路。

最后，现有方法的操作流程繁琐，需要多次在气动设计和结构设计间反复修正，必须进行大量的人工干预，实施效率较低。

发明内容

本发明解决的技术问题为：克服常规的从热态到冷态的叶型转换方法的不足，提供一种基于有限元逆向迭代的叶型转换方法，以准确获得叶片的冷态叶型，使工厂根据本方法提供的冷态叶型加工得到的叶片在投入实际使用时，其工作状态的热态叶型能够与气动设计人员给定的原始气动热态叶型一致，从而保证实际叶片能够忠实于设计者的意图。

本发明的技术方案为：所述的方法在已知叶轮机转子叶片热态叶型的前提下，基于有限元逆向迭代的叶型转换方法，获取叶片的冷态叶型。

作为本技术方案的一种改进，所述的方法包括如下步骤：

步骤一，获取叶片热态叶型坐标，并计算叶片的容差ε1、ε2和ε3；

步骤二，由流体计算软件对热态叶型进行三维计算流体动力学分析，获得叶片表面静压；

步骤三，将叶片热态叶型的流体网格和叶片表面静压转化为叶片热态叶型的结构网格上的气动力；

步骤四，由已知的热态叶型坐标和气动力，将叶片热态叶型结构网格转化为叶片冷态叶型的结构网格；

步骤五：将叶片由热态叶型坐标到模拟热态叶型坐标的有限元网格的变形，通过三维贴体插值得到模拟热态叶型的流体网格；

步骤六：由流体计算软件对所获得的模拟热态叶型的流体网格进行三维计算流体动力学分析；

步骤七：如果步骤六和步骤二两次流体动力学分析获得的压气机特性间的差异不满足容差，则用模拟热态叶型坐标代替热态叶型坐标，回到步骤三，否则计算结束，此时得到的叶片冷态叶型的结构网格即为所需要求的叶轮机转子叶片的冷态叶型的坐标。

作为本技术方案的一种改进，通过三维贴体插值，将叶片热态叶型的流体网格和叶片表面静压转化为叶片热态叶型的结构网格上的气动力。

作为本技术方案的一种改进，通过迭代法将叶片热态叶型结构网格转化为叶片冷态叶型的结构网格。

本发明的有益效果为：(1)常规方法都采用直接近似法处理叶型从热态到冷态的转换问题，对热态叶型实施强度分析获得变形，再将其从热态叶型中减去来获得近似的冷态叶型，这大幅降低了计算精度。本发明通过反复逆向迭代的方法随时修改结构刚度矩阵和气动载荷，克服了这一缺陷。

(2)常规方法采用了完全线性化假设，这对于弯扭掠复合的先进压气机叶片是不适合的，针对此类型叶片计算出的弹性恢复角的误差极大，这将导致所生产出的叶片在工作状态的热态叶型严重偏离原始气动设计。本发明在计算中引入了几何非线性效应，采用了非线性方法描述叶型从热态到冷态的转换过程，避免了大量的线性化假设对计算结果精度的影响。

(3)本发明采用了三维贴体形函数插值方法来传递气动载荷和结构变形，相对于常规的近似距离权函数插值法，可以提高插值精度和插值效率。

(4)常规方法的操作流程繁琐，需要多次在气动设计和结构设计间反复修正，必须进行大量的人工干预，实施效率较低，而本发明采用Fortran程序对整个转换过程进行了自动化操作，大大提高了计算效率，更加简便易行。

附图说明

图1是本发明的实现流程图；

图2是本发明的三维贴体插值流程图；

图3是本发明的三维贴体插值示意图；

图4是本发明的热态到冷态叶型迭代流程图；

图5是本发明实施例中的流体表面网格；

图6是本发明实施例中的热态到冷态迭代历程图；

图7是本发明实施例中的叶片变形曲线。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明是叶轮机转子叶片从热态到冷态的叶型转换方法，方法的实现流程如图1所示。所述的方法在已知叶轮机转子叶片热态叶型的前提下，基于有限元逆向迭代的叶型转换方法，获取叶片的冷态叶型。

步骤一：给定叶片热态叶型(其坐标为X_H)、计算容差ε1、ε2和ε3。

步骤二：由流体计算软件对热态叶型X_H进行三维计算流体动力学分析，获得叶片表面静压。

步骤三：将叶片热态叶型的流体网格和叶片表面静压通过三维贴体插值转化为叶片热态叶型的结构网格上的气动力F_p，实施流程如图2所示，其具体步骤如下：

步骤3.1：遍历所有耦合界面上的流体节点，找到距离结构节点s最近的流体节点f。

步骤3.2：找到所有包含流体节点f的流体表面网格单元c₁,c₂,…,c_m，并将这些流体表面网格单元沿垂直网格面方向进行拉伸，构造成为拥有8个节点的六面体临时单元h₁,h₂,…,h_m。

步骤3.3：对每个临时单元h_i，求出s相对于单元h_i的局部坐标(r_i,s_i,t_i)，并求出s与单元h_i中心的相对距离

步骤3.4：找到与s相对距离最小的单元h_t，即d_t＝min(d₁,d₂,…,d_m)，单元h_t即为结构节点s所对应的临时单元。

步骤3.5：假设结构节点s在流体表面网格单元上的投影点为s'，如图3所示。投影点s'在整体坐标系下的坐标为(x_s',y_s',z_s')，忽略流体节点与其投影点之间的距离，近似有x_s'≈x_s,y_s'≈y_s,z_s'≈z_s，基于有限元形函数理论，投影点s'相对于临时单元h_t的局部坐标(r,s,t)可以通过形函数插值得到

其中(x_i,y_i,z_i),1≤i≤8为临时单元h_t中8个节点在整体坐标系下的坐标，N_i为投影点处的形函数

这个方程为三元三次方程，直接求解比较困难，可以通过近似法求解，如拟牛顿法和梯度法等。

步骤3.6：假设结构节点s上的流体静压为p_s，投影点s'的位移为p_s'，近似有p_s≈p_s'，投影点s'的静压可以通过形函数插值得到，即式中p_i是临时单元h_t的8个节点i上的静压，由于构造临时单元h_t的8个节点中的4个临时节点上的形函数N＝0，可得由此就通过流体表面网格上各节点的静压插值获得了结构节点s的静压p_s，并将其转换为结构节点上的气动力F_p。

步骤四：由已知的热态叶型坐标X_H和气动力F_p，通过迭代将叶片热态叶型结构网格转化为叶片冷态叶型的结构网格，假设所计算得到的冷态叶型的坐标为X_C，实施流程如图4所示，具体步骤如下：

步骤4.1：设叶片的转速为Ω，则引入矩阵

步骤4.2：令i＝0，并假设当前冷态叶型坐标

步骤4.3：根据当前冷态叶型坐标求得结构的小位移刚度矩阵K₀、离心刚度矩阵K_C和所受的离心力F_C,分别为

$K_0=\underset{e}{\Σ}{K_0}^e=\underset{e}{\Σ}{\∫}_e{B_0}^T{\mathrm{DB}}_{0}dV$

$K_C=\underset{e}{\Σ}{K_C}^e=\underset{e}{\Σ}{\∫}_e{\mathrm{\ρN}}^T{\mathrm{\Ω}}^T\mathrm{\Ω}NdV$

$F_C=\underset{e}{\Σ}{F_C}^e={\∫}_e{\mathrm{\ρN}}^T{\mathrm{\Ω}}^T\mathrm{\Ω}{\left[\begin{array}{ccc}{x}_0& y_0& z_0\end{array}\right]}^{T}dV$

其中所有的积分都是在结构单元e上进行体积积分，式中B₀代表小位移刚度矩阵，D为弹性矩阵，ρ为密度，N为形函数矩阵，[x₀y₀z₀]^T为初始坐标。

步骤4.4：令n＝0，并假设初始位移

步骤4.5：根据a_n求得应力刚化矩阵K_σ、大变形刚度矩阵K_L和反力F^nr，分别为：

$K_{\mathrm{\σ}}=\underset{e}{\Σ}{K_{\mathrm{\σ}}}^e=\underset{e}{\Σ}{\∫}_e{G}^{T}MGdV$

$K_L=\underset{e}{\Σ}{K_L}^e=\underset{e}{\Σ}{\∫}_e({B_0}^T{\mathrm{DB}}_L+{B_L}^T{\mathrm{DB}}_0+{B_L}^T{\mathrm{DB}}_L)dV$

$F^{nr}=\∫{\stackrel{\‾}{B}}^T\mathrm{\σ}dV=\∫{\left(B+B_L\right)}^T\mathrm{\σ}dV$

其中：

$I_3=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$

$G=\left[\begin{array}{ccc}\frac{\∂N_1}{\∂x}{I}_3& \frac{\∂N_2}{\∂x}{I}_3& \mathrm{...}\\ \frac{\∂N_1}{\∂y}{I}_3& \frac{\∂N_2}{\∂y}{I}_3& \mathrm{...}\\ \frac{\∂N_1}{\∂z}{I}_3& \frac{\∂N_2}{\∂z}{I}_3& \mathrm{...}\end{array}\right]$

$M=\left[\begin{array}{ccc}{\mathrm{\σ}}_x{I}_3& {\mathrm{\τ}}_{xy}-I_3& {\mathrm{\τ}}_{xz}-I_3\\ {\mathrm{\τ}}_{yx}{I}_3& {\mathrm{\σ}}_y{I}_3& {\mathrm{\τ}}_{yz}-I_3\\ {\mathrm{\τ}}_{zx}{I}_3& {\mathrm{\τ}}_{zy}-I_3& {\mathrm{\σ}}_z{I}_3\end{array}\right]$

${\mathrm{\θ}}_x^T=\left[\begin{array}{ccc}\frac{\∂u}{\∂x}& \frac{\∂v}{\∂x}& \frac{\∂w}{\∂x}\end{array}\right]$

${\mathrm{\θ}}_y^T=\[\begin{array}{ccc}\frac{\∂u}{\∂y}& \frac{\∂v}{\∂y}& \frac{\∂w}{\∂y}\end{array}\]$

${\mathrm{\θ}}_z^T=\left[\begin{array}{ccc}\frac{\∂u}{\∂z}& \frac{\∂v}{\∂z}& \frac{\∂w}{\∂z}\end{array}\right]$

$A=\left[\begin{array}{ccc}{\mathrm{\θ}}_x^T& 0& 0\\ 0& {\mathrm{\θ}}_y^T& 0\\ 0& 0& {\mathrm{\θ}}_z^T\\ 0& {\mathrm{\θ}}_z^T& {\mathrm{\θ}}_y^T\\ {\mathrm{\θ}}_z^T& 0& {\mathrm{\θ}}_x^T\\ {\mathrm{\θ}}_y^T& {\mathrm{\θ}}_x^T& 0\end{array}\right]$

σ＝[σ_xσ_yσ_zτ_xyτ_yzτ_xz]^T

B_L＝AG

步骤4.6：根据(K₀-K_C+K_σ+K_L)da_n+1＝F_C+F_p-F^nr求出增量位移da_n+1。

步骤4.7：计算新位移a_n+1＝a_n+da_n+1。

步骤4.8：如果满足收敛条件||da_n+1||₂＜ε1||a_n+1||₂，则获得总增量刚度矩阵和位移a⁽ⁱ⁾＝a_n+1，否则令n＝n+1，前往步骤4.5。

步骤4.9：计算模拟热态叶型坐标和下一步新的冷态叶型坐标

步骤4.10：如果满足收敛条件则获得冷态叶型坐标和模拟热态叶型坐标否则令i＝i+1，前往步骤4.3。

步骤五：将叶片由热态叶型坐标X_H到模拟热态叶型坐标X'_H的有限元网格的变形通过三维贴体插值得到模拟热态叶型的流体网格，其方法与步骤三相似。

步骤六：由流体计算软件对所获得的模拟热态叶型X'_H的流体网格进行三维计算流体动力学分析。

步骤七：如果步骤六和步骤二两次流体动力学分析获得的压气机特性间的差异不满足容差ε3，则用X'_H代替X_H，回到步骤三，否则计算结束，此时得到的X_C即为所需要求的叶轮机转子叶片的冷态叶型的坐标。

现采用某航空发动机压气机转子叶片为实例，对本发明作进一步的详细说明。本例采用自编程序的形式实现发明内容，自编程序的开发平台是工程计算语言Fortran90，编译平台是IntelVisualFortran。例中流体动力学计算采用了商用软件NUMECA。

步骤1：由气动设计人员给出该叶片热态叶型的叶型坐标表。设定计算容差ε1＝1e-6、ε2＝1e-6和ε3＝1e-3。同时给定叶片材料参数和设计转速。采用自编程序的trans_NUMECA函数将叶型坐标表转化为geomTurbo几何文件。

步骤2：由商用软件NUMECA读入步骤1给出的geomTurbo几何文件，采用NUMECA的AutoGrid5模块进行网格划分得到流体网格，然后采用NUMECA的Fine模块进行参数设定和流体动力学计算，得到设计转速下的叶片表面静压分布和压气机特性(流量-总压比关系曲线、流量-效率关系曲线等)。由NUMECA的CFView模块给出叶片流体表面网格的网格坐标和表面静压，输出至文件surfpre.dat中。叶片流体表面网格如图5所示。

步骤3：调用自编程序的trans_fmesh_smesh函数，读入surfpre.dat文件，通过插值获得叶片的8节点六面体结构网格，存入数组coord(num_node,3)和elem(num_elem,8)中，其中num_node为结构网格的节点数，每个节点存入x、y和z三个方向的坐标值；num_elem为结构网格的单元数，每个单元存入8个节点的编号。然后采用三维贴体形函数插值将流体网格的表面静压插值到结构网格的节点上，形成结构节点力，存入列阵gfp中。

步骤4：采用自编程序的hot_cold函数，经过迭代获得了冷态叶型和模拟热态叶型，本例中经过了8次迭代达到收敛标准，收敛历程如图6所示。

步骤5：调用自编程序的trans_smesh_fmesh函数，采用三维贴体形函数插值，将模拟热态叶型的结构网格转化为气动造型输入的geomTurbo几何文件。

步骤6：由商用软件NUMECA读入步骤5给出的geomTurbo几何文件，采用NUMECA的AutoGrid5模块进行网格划分得到流体网格，然后采用NUMECA的Fine模块进行参数设定和流体动力学计算，得到设计转速下的叶片表面静压分布和压气机特性。由NUMECA的CFView模块给出叶片流体表面网格的网格坐标和表面静压，输出至文件surfpre_new.dat中。

步骤7：将步骤6计算获得的压气机特性与步骤2的计算结果进行比较。如果不满足容差，则用步骤4获得的模拟热态叶型代替给定热态叶型，返回步骤3继续进行计算，直至满足收敛条件为止。本例中经过2次迭代达到收敛标准。此时所输出的gongyi.dat文件即为所需要输出的冷态叶型的叶型坐标表，经过一定数据标准化处理后则成为提供给工厂加工使用的工艺文件。同时还输出了hot_cold.dat文件，在文件中给出了叶片从热态到冷态的几何变形，包括叶片的径向伸长量、周向挠度变化、轴向位置变化和弹性恢复角。该文件可以输出给Excel或Origin进行后处理，获得叶片的变形曲线，如图7所示。图中横坐标为变形量，纵坐标为叶高。

Claims (4)

1.一种叶轮机转子叶片从热态到冷态的叶型转换方法，其特征为：所述的方法在已知叶轮机转子叶片热态叶型的前提下，基于有限元逆向迭代的叶型转换方法，获取叶片的冷态叶型。

2.根据权利要求1所述的一种叶轮机转子叶片从热态到冷态的叶型转换方法，其特征为，所述的方法包括如下步骤：

步骤一，获取叶片热态叶型坐标，并计算叶片的容差ε1、ε2和ε3；

步骤二，由流体计算软件对热态叶型进行三维计算流体动力学分析，获得叶片表面静压；

步骤三，将叶片热态叶型的流体网格和叶片表面静压转化为叶片热态叶型的结构网格上的气动力；

步骤四，由已知的热态叶型坐标和气动力，将叶片热态叶型结构网格转化为叶片冷态叶型的结构网格；

步骤五：将叶片由热态叶型坐标到模拟热态叶型坐标的有限元网格的变形，通过三维贴体插值得到模拟热态叶型的流体网格；

步骤六：由流体计算软件对所获得的模拟热态叶型的流体网格进行三维计算流体动力学分析；

步骤七：如果步骤六和步骤二两次流体动力学分析获得的压气机特性间的差异不满足容差，则用模拟热态叶型坐标代替热态叶型坐标，回到步骤三，否则计算结束，此时得到的叶片冷态叶型的结构网格即为所需要求的叶轮机转子叶片的冷态叶型的坐标。

3.根据权利要求2所述的一种叶轮机转子叶片从热态到冷态的叶型转换方法，其特征为：通过三维贴体插值，将叶片热态叶型的流体网格和叶片表面静压转化为叶片热态叶型的结构网格上的气动力。

4.根据权利要求2所述的一种叶轮机转子叶片从热态到冷态的叶型转换方法，其特征为：通过迭代法将叶片热态叶型结构网格转化为叶片冷态叶型的结构网格。