CN101908088A - 一种基于时域双向迭代的叶轮机叶片颤振应力预测方法 - Google Patents

Abstract

本发明公开了一种基于时域双向迭代的叶轮机叶片颤振应力预测方法，其特征在于，将叶片及其周围流场视为一个三维流固耦合系统，设计了一套时域的双向迭代方法，通过交替求解叶片变形和非定常流场得到叶片的颤振应力。该方法在计算机中设定以下模块：结构计算模块、流体计算模块、数据转换模块、颤振应力输出模块、初值计算模块和双向迭代模块。通过非线性迭代获得叶片静变形和稳态流场作为初值；交替调用结构计算模块和流体计算模块在时间上推进整个系统；通过数据转换模块传递流固边界信息；输出时间历程上的颤振应力。本发明实现了叶片与流场的一体化计算，考虑了耦合系统的非线性，能观察整个颤振发展进程，预测叶片的颤振应力。

Description

一种基于时域双向迭代的叶轮机叶片颤振应力预测方法

(一)技术领域

本发明涉及一种计算机辅助分析工具领域的方法，具体是涉及一种基于时域双向迭代的叶轮机叶片颤振应力预测方法。属于叶轮机械模拟技术领域。

(二)背景技术

目前，在叶轮机械中，存在着振动叶片与周围流体之间复杂的相互作用，弹性叶片通过其变形影响流场，流场的扰动反之改变叶片的形状。在这类流固耦合现象中，以颤振对结构可靠性的影响最大。颤振属于一种自激振动，叶片表面的非定常气动力来源于叶片本身的运动，非定常气动力和叶片的非定常运动互为因果、循环递增，使得短时间内叶片振幅急剧增大，甚至造成叶片断裂和机器损坏。由于进行颤振试验难度较高、成本巨大且具有相当的危险性，所以随着计算机及相关学科的发展，数值模拟方法因其成本低、效率高，成为颤振分析的重要工具。

在叶轮机械的颤振数值模拟中，通常使用能量法和特征值法，近年来也少量开始使用时间推进法。现有的颤振数值计算方法至少存在以下四个缺点：

首先，传统的能量法和特征值法都假定叶片以给定形式作运动，将流场和叶片运动互相孤立，忽略了弹性叶片本身变形对流场的影响，同时也极少考虑离心力对叶片变形和流场的影响，无法实现气动-结构一体化计算，存在较大的计算误差。

其次，由于流体和结构都存在强烈的非线性，这两种领域的非线性交替叠加，形成了更加复杂的耦合非线性。现有的颤振预测方法都引入了大量线性化假设，甚至只采用结构向流体的单向传递，弱化甚至消除了耦合非线性，无法准确地描述非线性现象，降低了模拟精度。

再次，现有的方法不能展示整个颤振发展的时间历程，也无法得到颤振应力，而颤振应力恰恰是颤振现象中非常重要的指标。

最后，以往的基于时间推进的叶轮机颤振数值模拟方法都采用固定时间步长，无法在计算中对步长进行有效的调整，造成计算速度和效率较低，甚至会导致收敛困难；而且大多采用商用的结构有限元计算软件，在计算中需要反复进行人工干预，不能在软件内部针对流固耦合计算的特性对整个计算流程和算法进行优化，增加了操作复杂度，降低了计算效率，使得现有的基于时间推进的预测方法计算费用较高，经济型较差。

(三)发明内容

1、目的：本发明的目的是克服常规的叶轮机颤振数值仿真方法的不足，提供一种基于时域双向迭代的叶轮机叶片颤振应力预测方法，使之能够准确描述流体与叶片之间的相互作用和耦合非线性现象，并对叶片的颤振应力进行预测。

2、技术方案：本发明一种基于时域双向迭代的叶轮机叶片颤振应力预测方法，它是在计算机上采用数值模拟的方法对叶轮机叶片的颤振应力进行预测。该方法具体步骤如下：

步骤一：在计算机中设定以下六个模块：结构计算模块、流体计算模块、数据转换模块、颤振应力输出模块、初值计算模块和双向迭代模块。

步骤二：调用初值计算模块通过非线性迭代获得叶片的初始静变形和稳态定常流场作为以后双向迭代计算的初值，用以加快整个系统的收敛。

步骤三：调用双向迭代模块，在时间上推进由叶片和周围流场组成的流固耦合系统。

步骤四：通过颤振应力输出模块，将结构的瞬态位移响应转化为瞬态颤振应力，并输出为文件，供后处理软件读取显示。

其中，步骤一所述的六个模块的具体结构如下：

模块一：结构计算模块：

结构计算模块的作用是求解三维结构动力方程，获得叶片的瞬态变形。该模块包含如下步骤：

步骤1.1：将叶片表面压力和结构有限元节点单元信息放到输入文件里，并在输入文件中对叶片的工况进行设置。

步骤1.2：调用Fortran语言编写的结构有限元计算程序求解三维结构动力方程

$M\stackrel{\·\·}{a}+C\stackrel{\·}{a}+(K-K_c+S)a=f+f_a+f_c---\left(1\right)$

其中M，C和K分别是总体质量矩阵、阻尼矩阵和刚度矩阵。K_c是旋转软化矩阵。S是应力刚化矩阵(或称几何刚化矩阵)。f，f_a和f_c是节点外力列阵、气动力列阵和离心力列阵。a是待求的有限元节点位移列向量。

步骤1.3：计算结束后，输出叶片表面有限元节点的位移。

模块二：流体计算模块：

在流体计算模块中，提供了各种流体计算软件的接口，可以后台调用流体计算软件获得每个时间步的三维非定常流场。该模块中根据用户需要调用各种流体动力学计算软件，如Fluent，CFX等，在计算结束后输出叶片表面压力。

模块三：数据转换模块：

数据转换模块包括三个子模块，一是节点配对子模块，作用是找到耦合界面上每个流体节点所对应的结构单元；二是载荷传递子模块，该子模块通过三维形函数插值，将叶片表面的定常或非定常压力转换为叶片结构有限元模型上的节点力；三是变形传递子模块，该子模块将叶片的变形转换为叶片周围流体域的边界运动。下面分别对这三个子模块进行介绍：

节点配对子模块包含下列步骤：

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

步骤3.1.2：找到所有包含结构节点s的结构单元e₁，e₂，…，e_m。

步骤3.1.3：对每个结构单元e_i(i＝1，2，…，m)，用拟牛顿法求出f在该单元内的局部坐标(r_i，s_i，t_i)，并求出f与该单元中心的相对距离

步骤3.1.4：找出与f相对距离最小的单元e_t，即d_t＝min(d₁，d₂，…，d_m)，该单元即为流体节点f所对应的结构单元。

变形传递子模块包含下列步骤：

步骤3.2.1：设流体节点f在其所对应结构单元e上的投影点为f’，f’在整体坐标系下的坐标为(x_f’，y_f’，z_f’)，忽略流体节点与其投影点之间的距离，近似有x_f′≈x_f，y_f′≈y_f，z_f′≈z_f，投影点f’在结构单元e中的局部坐标(r，s，t)通过求解下列方程组得到

$\left\{\begin{array}{c}{x}_{f^{\′}}=\underset{i=1}{\overset{8}{\mathrm{\Σ}}}{N}_i{x}_i\\ y_{f^{\′}}=\underset{i=1}{\overset{8}{\mathrm{\Σ}}}{N}_i{y}_i\\ z_{f^{\′}}=\underset{i=1}{\overset{8}{\mathrm{\Σ}}}{N}_i{z}_i\end{array}\right.---\left(2\right)$

其中n为结构单元e中节点的数目，x_i，y_i，z_i(i＝1，2，…，8)为结构单元e中8个节点在整体坐标系下的坐标。投影点处的形函数N_i(i＝1，2，…，8)为

N₁＝(1-s)(1-t)(1-r)/8，N₂＝(1+s)(1-t)(1-r)/8

N₃＝(1+s)(1+t)(1-r)/8，N₄＝(1-s)(1+t)(1-r)/8

                                             (3)

N₅＝(1-s)(1-t)(1+r)/8，N₆＝(1+s)(1-t)(1+r)/8

N₇＝(1+s)(1+t)(1+r)/8，N₈＝(1-s)(1+t)(1+r)/8

步骤3.2.2：在求出投影点f’在结构单元e中的局部坐标(r，s，t)后，如果投影点在面r＝1上，忽略流体节点与其对应投影点之间的距离，强制令r＝1；如果投影点在面s＝1上，忽略流体节点与其对应投影点之间的距离，强制令s＝1；如果投影点在面t＝1上，忽略流体节点与其对应投影点之间的距离，强制令t＝1。

步骤3.2.3：假设投影点f’的位移为u_f’，u_i为结构有限元节点i的位移，通过下列公式求解流体节点f的位移u_f：

$u_f\≈u_{f^{\′}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{N}_i{u}_i---\left(4\right)$

载荷传递子模块包含下列步骤：

步骤3.3.1：设流体表面网格c的中心点为f_c，计算流体压力沿网格c的积分

其中Ω为网格c的表面，p为耦合界面上的流体压力，

为中心点f_c上的压力，A_c为网格c的表面积。

步骤3.3.2：设流体表面网格中心点f_c对应的结构单元为e，通过计算

(i＝1，2，…，8)可以得到结构单元e上8个节点的节点力，其中N_i为f_c在结构单元e上投影点处的形函数。

步骤3.3.3：将结构单元e上8个节点的节点力集合入结构有限元节点载荷列阵。

模块四：颤振应力输出模块：

颤振应力输出模块通过叶片的瞬态位移响应计算出叶片的瞬态颤振应力，并输出各种供颤振分析使用的后处理文件。该模块包括如下步骤：

步骤4.1：通过叶片的瞬态位移求出叶片的瞬态颤振应力，即

σ＝Dε    (5)

其中ε为应变列阵，σ为应力列阵，D为弹性矩阵。

步骤4.2：将叶片随时间变化的位移、颤振应力和压力输出为可被Origin软件读取的格式。

步骤4.3：将转子随时间变化的质量流量和压比输出为可被Origin软件读取的格式。

步骤4.4：将叶片瞬间的位移云图和颤振应力云图输出为可被Tecplot软件读取的格式。

步骤4.5：如果整个计算结束，输出该级转子叶片在压气机特性图上的颤振边界。

模块五：初值计算模块：

初值计算模块通过反复的非线性迭代获得叶片的初始静变形和稳态定常流场，作为以后双向迭代计算的初值，用以加快整个系统的收敛。该模块包含如下步骤：

步骤5.1：首先通过三维模型设计软件UG建立叶片的三维实体模型；然后将实体模型导入网格划分软件HyperMesh中，得到叶片的结构有限元网格；将实体模型导入网格划分软件TurboGrid中，得到叶轮机内流场的流体计算网格。

步骤5.2：调用结构计算模块得到叶片在离心力作用下的初始位移

其中K₀是结构的小位移刚度矩阵，F_c为叶片在转动中的离心力列阵。

步骤5.3：调用数据转换模块将叶片变形转化为流体域的边界运动。

步骤5.4：通过流体计算模块调用流体动力学计算软件，进行网格更新，然后计算新的定常流场。

步骤5.5：调用数据转换模块将叶片表面的定常压力转换为结构有限元模型上的节点力。

步骤5.6：调用结构计算模块获得叶片的非线性变形，并对结构的刚度矩阵进行更新。计算叶片非线性变形的具体步骤如下：

步骤5.6.1：根据步骤5.2求得的初始位移a₀，求得应力刚化矩阵K_s、大变形刚度矩阵K_L和反力F^nr。

步骤5.6.2：根据(K₀-K_C+K_σ+K_L)da₁＝F_C-F^nr求出增量位移da₁。

步骤5.6.3：计算总位移a₁＝a₀+da₁。

步骤5.6.4：如果求得的a_n满足收敛条件||da_n||₂＜ε||a_n||₂，其中ε取0.001，则非线性迭代收敛，计算结束，指向步骤5.7；否则用a₁代替a₀，前往步骤5.3。

步骤5.7：求得新的结构刚度矩阵K＝K₀-K_C+K_S+K_L。常规的有限元计算软件无法同时考虑几何非线性和应力刚化的影响，本发明通过步骤5.7即时对结构刚度矩阵进行修正，从而克服了这一缺陷。此时的叶片静变形和定常流场就作为以后双向迭代计算的初始值。

模块六：双向迭代模块：

双向迭代模块反复交替调用结构求解器和流体求解器，同时采用数据转换模块在结构求解器和流体求解器之间传递边界数据，通过流体与结构双向的交互作用，在时间上推进整个由叶片及其周围流场组成的流固耦合系统。双向迭代模块包含下列步骤：

步骤6.1：通过数据转换系统将叶片的变形转换为流体计算域的边界运动。

步骤6.2：通过流体计算模块调用流体动力学计算软件，进行网格更新，然后计算新的非定常流场。

步骤6.3：通过数据转换模块将叶片表面的非定常压力转换为结构有限元模型上的节点力。

步骤6.4：通过结构计算模块求得新的叶片变形。

步骤6.5：对双向迭代的计算时间步进行调整。本发明提供了一种新的流固耦合自适应时间步技术，通过当前的收敛速度控制下一步计算的时间步长，随时对后续计算的时间步长进行调整，从而有效地缩短收敛时间，提高计算效率。该流固耦合自适应时间步技术包括下列步骤：

步骤6.5.1：计算其中u_i为当前时间步计算的节点i在某个自由度的位移，u_i’为上个时间步计算的节点i在某个自由度的位移，n为有限元节点的数目。

步骤6.5.2：取各个自由度的d的平均值为d₀，用以表征当前时间步的收敛速度。

步骤6.5.3：用e代表收敛速度的阈值，用t代表当前的时间步长，如果d₀＜0.5e，则新的时间步长t_new＝2t；如果d₀＞2e，则新的时间步长t_new＝0.5t；如果0.5e＜d₀＜2e，则时间步长不变，t_new＝t。

步骤6.5.4：在获得流固耦合计算的物理时间步长t_new后，对结构动力方程瞬态求解的时间步长和流体非定常求解的时间步长作相应调整。

步骤6.6：通过计算所得的叶片瞬态位移响应，观察颤振的发展进程。如果响应曲线衰减，则不会发生颤振，计算结束；如果响应曲线增长，发生颤振，指向步骤6.1；如果发生颤振且计算达到预先指定时间点，计算结束，指向步骤6.7输出颤振应力。

步骤6.7：通过颤振应力输出模块，将结构的瞬态位移响应转化为瞬态颤振应力，并输出为文件，供后处理软件读取显示。

3、本发明的优点在于：

(1)在计算开始前，加入一个初值计算的过程。初值计算中引入了非设计状态时叶片离心力对结构刚度矩阵和流场的影响，通过反复的非线性迭代获得初值，从而加快了整个计算的收敛速度；同时，常规的有限元计算软件无法同时考虑几何非线性和应力刚化的影响，本发明通过非线性迭代随时对结构刚度矩阵进行修改，克服了这一缺陷。

(2)常规方法都假定叶片以指定形式作运动，将流场和叶片运动互相孤立，忽略了弹性叶片本身变形对流场的影响，本发明将流场和叶片运动结合起来，成为一个整体的流固耦合系统，采用一个双向迭代的流程，实现了气动-结构一体化计算，减小了计算误差。

(3)通过双向迭代可以准确地描述结构非线性、流体非线性和流固耦合非线性，避免了单向传递或大量的线性化假设对模拟结果精度的影响。

(4)本发明通过时间推进提供整个颤振发展的时间历程图，并得到各个时间点的颤振应力。一方面可以验证实际中颤振现象的发生发展过程，为颤振机理的阐释提供依据；另一方面可以为颤振试验提供指导和依据，预测达到预期颤振应力所需的时间，指导试验及时停车，防止由于叶片损毁甚至飞脱碰撞造成严重的试验事故。

(5)以往的基于时域的叶轮机颤振数值模拟方法都采用固定时间步长，不仅造成收敛困难，而且计算速度和效率较低，本发明中发展了一种流固耦合自适应时间步技术，通过当前的收敛速度控制下一步计算的时间步长，随时对后续计算的时间步长进行调整，从而有效地缩短收敛时间，提高计算效率。同时，采用全自动化的操作流程，通过简单的输入文件进行控制，无需在计算过程中进行大量人工干预，避免了繁琐的操作。另外，通过自编的结构有限元求解程序，在程序内部针对流固耦合计算的特性对整个计算流程和算法进行优化，大大提高了计算效率和经济性。

(四)附图说明

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

图2是本发明的系统模块图；

图3是本发明的初值计算模块的流程图；

图4是本发明所用某航空发动机压气机叶片及流场模型；

图5是本发明的双向迭代与现有方法的区别示意图；

图6是本发明的流固耦合自适应时间步技术的示意图；

图7是本发明叶根某点颤振应力的时间历程图；

图8是本发明预测压气机颤振特性图；

图9是本发明叶根某点静压的时间历程图；

图10是本发明叶片瞬间颤振应力分布图。

(五)具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。本发明是基于时域双向迭代的叶轮机叶片颤振应力预测方法，方法的实现流程如图1所示。

步骤一：在计算机中设定以下六个模块：结构计算模块、流体计算模块、数据转换模块、颤振应力输出模块、初值计算模块和双向迭代模块。

步骤二：调用初值计算模块通过非线性迭代获得叶片的初始静变形和稳态定常流场作为以后双向迭代计算的初值，用以加快整个系统的收敛。

步骤三：调用双向迭代模块，在时间上推进由叶片和周围流场组成的流固耦合系统。

步骤四：通过颤振应力输出模块，将结构的瞬态位移响应转化为瞬态颤振应力，并输出为文件，供后处理软件读取显示。

步骤一中提到的六个模块包括结构计算模块、流体计算模块、数据转换模块、颤振应力输出模块、初值计算模块和双向迭代模块，如图2所示。这六个模块并非互相孤立，而是互相联系，在使用其中一个模块时，可能会在该模块中调用其他模块。六个模块的结构如下：

模块一：结构计算模块：

结构计算模块的作用是求解三维结构动力方程，获得叶片的瞬态变形。该模块按照下列步骤进行运算求解：

步骤1.1：将叶片表面压力和结构有限元节点单元信息放到输入文件里，并在输入文件中对叶片的工况进行设置。

步骤1.2：调用Fortran语言编写的结构有限元计算程序求解三维结构动力方程获得叶片的瞬态变形。

步骤1.3：计算结束后，输出叶片表面有限元节点的位移。

模块二：流体计算模块：

在流体计算模块中，提供了各种流体计算软件的接口，可以后台调用流体计算软件获得每个时间步的三维非定常流场。该模块包含如下步骤：

步骤2.1：按照需要选择求解使用的流体动力学计算软件，如Fluent，CFX等。

步骤2.2：在后台调用所选择的流体动力学计算软件进行求解。

步骤2.3：计算结束后，输出叶片表面的压力。

模块三：数据转换模块：

数据转换模块包括三个子模块，一是节点配对子模块，作用是找到耦合界面上每个流体节点所对应的结构单元；二是载荷传递子模块，该子模块通过三维形函数插值，将叶片表面的定常或非定常压力转换为叶片结构有限元模型上的节点力；三是变形传递子模块，该子模块将叶片的变形转换为叶片周围流体域的边界运动。

模块四：颤振应力输出模块：

颤振应力输出模块通过叶片的瞬态位移响应计算出叶片的瞬态颤振应力，并输出各种供颤振分析使用的后处理文件。该模块包括如下步骤：

步骤4.1：通过叶片的瞬态位移求出叶片的瞬态颤振应力。

步骤4.2：将叶片随时间变化的位移、颤振应力和压力输出为可被Origin软件读取的格式。

步骤4.3：将转子随时间变化的质量流量和压比输出为可被Origin软件读取的格式。

步骤4.4：将叶片瞬间的位移云图和颤振应力云图输出为可被Tecplot软件读取的格式。

步骤4.5：如果整个计算结束，输出该级转子叶片在压气机特性图上的颤振边界。

模块五：初值计算模块：

初值计算模块通过反复的非线性迭代获得叶片的初始静变形和稳态定常流场，作为以后双向迭代计算的初值，用以加快整个系统的收敛。该模块包含如下步骤：

步骤5.1：首先通过三维模型设计软件UG建立叶片的三维实体模型；然后将实体模型导入网格划分软件HyperMesh中，得到叶片的结构有限元网格；将实体模型导入网格划分软件TurboGrid中，得到叶轮机内流场的流体计算网格。

步骤5.2：调用结构计算模块得到叶片在离心力作用下的初始位移列阵

其中K₀是结构的小位移刚度矩阵，F_C为叶片在转动中的离心力列阵。

步骤5.3：调用数据转换模块将叶片变形转化为流体域的边界运动。

步骤5.4：通过流体计算模块调用流体动力学计算软件，进行网格更新，然后计算新的定常流场。

步骤5.5：调用数据转换模块将叶片表面的定常压力转换为结构有限元模型上的节点力。

步骤5.6：调用结构计算模块获得叶片的非线性变形，并对结构的刚度矩阵进行更新。如果非线性迭代收敛，计算结束，指向步骤5.7；否则用新的叶片位移列阵取代旧的叶片位移列阵，前往步骤5.3。

步骤5.7：求得新的结构刚度矩阵K＝K₀-K_C+K_S+K_L。常规的有限元计算软件无法同时考虑几何非线性和应力刚化的影响，本发明通过步骤5.7即时对结构刚度矩阵进行修正，从而克服了这一缺陷。此时的叶片静变形和定常流场就作为以后双向迭代计算的初始值。

模块六：双向迭代模块：

双向迭代模块反复交替调用结构求解器和流体求解器，同时采用数据转换模块在结构求解器和流体求解器之间传递边界数据，通过流体与结构双向的交互作用，在时间上推进整个由叶片及其周围流场组成的流固耦合系统。双向迭代模块包含下列步骤：

步骤6.1：通过数据转换系统将叶片的变形转换为流体计算域的边界运动。

步骤6.2：通过流体计算模块调用流体动力学计算软件，进行网格更新，然后计算新的非定常流场。

步骤6.3：通过数据转换模块将叶片表面的非定常压力转换为结构有限元模型上的节点力。

步骤6.4：通过结构计算模块求得新的叶片变形。

步骤6.5：对双向迭代的计算时间步进行调整。本发明提供了一种新的流固耦合自适应时间步技术，通过当前的收敛速度控制下一步计算的时间步长，随时对后续计算的时间步长进行调整，从而有效地缩短收敛时间，提高计算效率。该流固耦合自适应时间步技术包括下列步骤：

步骤6.5.1：计算

其中u_i为当前时间步计算的节点i在某个自由度的位移，u_i’为上个时间步计算的节点i在某个自由度的位移，n为有限元节点的数目。

步骤6.5.2：取各个自由度的d的平均值为d₀，用以表征当前时间步的收敛速度。

步骤6.5.3：用e代表收敛速度的阈值，用t代表当前的时间步长，如果d₀＜0.5e，则新的时间步长t_new＝2t；如果d₀＞2e，则新的时间步长t_new＝0.5t；如果0.5e＜d₀＜2e，则时间步长不变，t_new＝t。

步骤6.5.4：在获得流固耦合计算的物理时间步长t_new后，对结构动力方程瞬态求解的时间步长和流体非定常求解的时间步长作相应调整。

步骤6.6：通过计算所得的叶片瞬态位移响应，观察颤振的发展进程。如果响应曲线衰减，则不会发生颤振，计算结束；如果响应曲线增长，发生颤振，指向步骤6.1；如果发生颤振且计算达到预先指定时间点，计算结束，指向步骤6.7输出颤振应力。

步骤6.7：通过颤振应力输出模块，将结构的瞬态位移响应转化为瞬态颤振应力，并输出为文件，供后处理软件读取显示。

现采用某航空发动机压气机叶片为实例，对本发明作进一步的详细说明。本例中使用的设备为一台Intel内核的微型计算机，主频为2.8GHz，内存为2G，操作系统为Microsoft Windows XP Professional，Service Pack 2。

步骤1：初值计算

调用初值计算模块以加快收敛，其流程如图3所示，具体步骤如下：

步骤1.1：首先通过三维模型设计软件UG建立叶片的三维实体模型；然后将实体模型导入网格划分软件HyperMesh中，得到叶片的结构有限元网格；将实体模型导入网格划分软件TurboGrid中，得到叶轮机内流场的流体计算网格；叶片的结构有限元网格和流体网格见图4，为清晰起见，流体只显示了轮毂和叶片表面的网格；

步骤1.2：通过结构计算模块调用Fortran语言编写的结构有限元分析程序，得到叶片在离心力作用下的变形，将叶片表面有限元节点的位移写入文件bladenode.txt；

步骤1.3：将文件bladenode.txt输入数据转换模块，将叶片变形转化为流体域的边界运动，生成叶片表面流体网格点的位移，写入文件bladepoint.txt；

步骤1.4：通过流体计算模块调用流体动力学计算软件Fluent，读入文件bladepoint.txt，进行动网格计算，获得更新后的流场网格，并计算新的定常流场，将叶片表面流体网格点的压力写入文件bladepres.txt；

步骤1.5：将文件bladepres.txt输入数据转换模块，将叶片表面的定常压力转换为结构有限元模型上的节点力，写入文件bladeforce.txt；

步骤1.6：通过结构计算模块调用由标准Fortran语言编写的结构有限元分析程序，读入节点力文件bladeforce.txt，进行非线性迭代获得叶片的变形和新的结构刚度矩阵，将叶片表面有限元节点的变形写入文件bladenode.txt，将更新后的结构刚度矩阵写入文件restartek.dat；

步骤1.7：如果非线性迭代收敛，计算结束，所得叶片静变形和定常流场就作为以后双向迭代计算的初始值；否则如果非线性迭代不收敛，指向步骤1.3。

步骤2：双向迭代

调用双向迭代模块，在时间上推进整个流固耦合系统。与常规的单向迭代不同，本发明通过流体与结构之间的双向迭代实现了气动-结构一体化计算，常规方法与本发明所采用方法的区别见图5。双向迭代的具体步骤如下：

步骤2.1：将文件bladenode.txt输入数据转换模块，采用三维形函数插值方法，将叶片的变形转换为流体计算域的边界运动，将叶片表面流体网格点的位移写入文件bladepoint.txt：

步骤2.2：通过流体计算模块调用流体动力学计算软件Fluent，读入文件bladepoint.txt，进行动网格计算，获得更新后的流场网格，并计算新的三维非定常流场，将叶片表面流体网格点的压力写入文件bladepres.txt；

步骤2.3：将文件bladepres.txt输入数据转换模块，采用三维形函数插值方法，将叶片表面的非定常压力转换为结构有限元模型上的节点力，写入文件bladeforce.txt；

步骤2.4：通过结构计算模块调用Fortran语言编写的结构有限元分析程序，读入节点力文件bladeforce.txt，通过隐式Newmark法求解三维结构动力方程，获得叶片的瞬态变形，并将叶片表面有限元节点的位移写入文件bladenode.txt；

步骤2.5：在双向迭代中，运用本发明提供的流固耦合自适应时间步技术，通过当前的收敛速度控制下一步计算的时间步长，随时对后续计算的时间步长进行调整，从而有效地缩短收敛时间，提高计算效率。本发明采用的流固耦合自适应时间步技术的示意图见图6。

步骤2.6：通过计算所得的叶片瞬态位移响应，观察颤振的发展进程。如果响应曲线衰减，则不会发生颤振，计算结束；如果响应曲线增长，则颤振发生，指向步骤2.1；如果颤振发生且计算达到预先指定的时间点，计算结束，指向步骤3输出颤振应力。以叶根某点为例，在图7中给出了颤振应力的时间历程图。由图中可以看出，在极短的时间内叶片振幅急剧增大，严重威胁到了发动机的运行安全，这符合以往观察到的颤振的特征。由于本发明可以展示整个颤振发展的时间历程，并得到颤振应力，这突破了现有的颤振数值模拟中的局限。其具体意义如下：通过本发明提供的模拟颤振发展历程图，可以观察整个颤振发生和发展的细节，从而通过数值模拟方法验证了实际中的颤振现象，为颤振机理的阐释提供了依据；通过颤振发展历程图和计算所得颤振应力，可以为颤振试验提供指导和依据，预测达到预期颤振应力所需的时间，指导试验及时停车，防止由于叶片损毁甚至飞脱碰撞造成严重的试验事故。

步骤3：结果输出

通过颤振应力输出模块，将结构的瞬态位移响应转化为瞬态颤振应力，同时输出一系列供颤振分析使用的后处理文件。所有输出的数据文件包括：

步骤3.1：输出该级转子叶片在压气机特性图上的颤振边界，如图8所示，同时在图中附上了试验所得的颤振边界以作比较。由图中可以看出，数值预测的结果与试验的结果比较吻合；

步骤3.2：输出叶片节点位移、叶片节点颤振应力、叶片表面静压、压气机质量流量、压气机压比的时间历程图；以叶根某点静压为例，其时间历程图见图9；

步骤3.3：输出每个时间步叶片表面的位移云图和颤振应力云图，并将所有时间步的云图连接起来形成动画；以某时刻叶片表面颤振应力为例，其云图见图10。通过应力云图可以获得颤振时叶片应力集中的危险位置，为后续分析提供依据。

Claims (2)

1.一种基于时域双向迭代的叶轮机叶片颤振应力预测方法，其特征在于，该方法具体步骤如下：

步骤一：在计算机中设定以下六个模块：结构计算模块、流体计算模块、数据转换模块、颤振应力输出模块、初值计算模块和双向迭代模块；

步骤二：调用初值计算模块通过非线性迭代获得叶片的初始静变形和稳态定常流场作为以后双向迭代计算的初值，用以加快整个系统的收敛；

步骤三：调用双向迭代模块，在时间上推进由叶片和周围流场组成的流固耦合系统；

步骤四：通过颤振应力输出模块，将结构的瞬态位移响应转化为瞬态颤振应力，并输出为文件，供后处理软件读取显示。

2.根据权利要求1所述的一种基于时域双向迭代的叶轮机叶片颤振应力预测方法，其特征在于：步骤一中所述的六个模块的具体结构如下：

模块一：结构计算模块：

结构计算模块的作用是求解三维结构动力方程，获得叶片的瞬态变形；该模块包含如下步骤：

步骤1.1：将叶片表面压力和结构有限元节点单元信息放到输入文件里，并在输入文件中对叶片的工况进行设置；

步骤1.2：调用Fortran语言编写的结构有限元计算程序求解三维结构动力方程

$M\stackrel{\·\·}{a}+C\stackrel{\·}{a}+(K-K_c+S)a=f+f_a+f_c---\left(1\right)$

其中M，C和K分别是总体质量矩阵、阻尼矩阵和刚度矩阵；K_c是旋转软化矩阵；S是应力刚化矩阵即几何刚化矩阵；f，f_a和f_c是节点外力列阵、气动力列阵和离心力列阵，a是待求的有限元节点位移列向量；

步骤1.3：计算结束后，输出叶片表面有限元节点的位移；

模块二：流体计算模块：

在流体计算模块中，提供了各种流体计算软件的接口，可以后台调用流体计算软件获得每个时间步的三维非定常流场；该模块中根据用户需要调用流体动力学计算软件Fluent，CFX，在计算结束后输出叶片表面压力；

模块三：数据转换模块：

数据转换模块包括三个子模块，一是节点配对子模块，作用是找到耦合界面上每个流体节点所对应的结构单元；二是载荷传递子模块，该子模块通过三维形函数插值，将叶片表面的定常或非定常压力转换为叶片结构有限元模型上的节点力；三是变形传递子模块，该子模块将叶片的变形转换为叶片周围流体域的边界运动；

节点配对子模块包含下列步骤：

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

步骤3.1.2：找到所有包含结构节点s的结构单元e₁，e₂，…，e_m；

步骤3.1.3：对每个结构单元e_i(i＝1，2，…，m)，用拟牛顿法求出f在该单元内的局部坐标(r_i，s_i，t_i)，并求出f与该单元中心的相对距离

步骤3.1.4：找出与f相对距离最小的单元e_t，即d_t＝min(d₁，d₂，…，d_m)，该单元即为流体节点f所对应的结构单元；

变形传递子模块包含下列步骤：

步骤3.2.1：设流体节点f在其所对应结构单元e上的投影点为f’，f’在整体坐标系下的坐标为(x_f’，y_f’，z_f’)，忽略流体节点与其投影点之间的距离，近似有x_f′≈x_f，y_f′≈y_f，z_f′≈z_f，投影点f’在结构单元e中的局部坐标(r，s，t)通过求解下列方程组得到

$\left\{\begin{array}{c}{x}_{f^{\′}}=\underset{i=1}{\overset{8}{\mathrm{\Σ}}}{N}_i{x}_i\\ y_{f^{\′}}=\underset{i=1}{\overset{8}{\mathrm{\Σ}}}{N}_i{y}_i\\ z_{f^{\′}}=\underset{i=1}{\overset{8}{\mathrm{\Σ}}}{N}_i{z}_i\end{array}\right.---\left(2\right)$

其中n为结构单元e中节点的数目，x_i，y_i，z_i(i＝1，2，…，8)为结构单元e中8个节点在整体坐标系下的坐标；投影点处的形函数N_i(i＝1，2，…，8)为

N₁＝(1-s)(1-t)(1-r)/8，N₂＝(1+s)(1-t)(1-r)/8

N₃＝(1+s)(1+t)(1-r)/8，N₄＝(1-s)(1+t)(1-r)/8

                                            (3)

N₅＝(1-s)(1-t)(1+r)/8，N₆＝(1+s)(1-t)(1+r)/8

N₇＝(1+s)(1+t)(1+r)/8，N₈＝(1-s)(1+t)(1+r)/8

步骤3.2.2：在求出投影点f’在结构单元e中的局部坐标(r，s，t)后，如果投影点在面r＝1上，忽略流体节点与其对应投影点之间的距离，强制令r＝1；如果投影点在面s＝1上，忽略流体节点与其对应投影点之间的距离，强制令s＝1；如果投影点在面t＝1上，忽略流体节点与其对应投影点之间的距离，强制令t＝1；

步骤3.2.3：假设投影点f’的位移为u_f’，ui为结构有限元节点i的位移，通过下列公式求解流体节点f的位移u_f：

$u_f\≈u_{f^{\′}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{N}_i{u}_i---\left(4\right)$

载荷传递子模块包含下列步骤：

步骤3.3.1：设流体表面网格c的中心点为f_c，计算流体压力沿网格c的积分

其中Ω为网格c的表面，p为耦合界面上的流体压力，

为中心点f_c上的压力，A_c为网格c的表面积；

步骤3.3.2：设流体表面网格中心点f_c对应的结构单元为e，通过计算(i＝1，2，…，8)可以得到结构单元e上8个节点的节点力，其中N_i为f_c在结构单元e上投影点处的形函数；

步骤3.3.3：将结构单元e上8个节点的节点力集合入结构有限元节点载荷列阵；

模块四：颤振应力输出模块：

颤振应力输出模块通过叶片的瞬态位移响应计算出叶片的瞬态颤振应力，并输出各种供颤振分析使用的后处理文件；该模块包括如下步骤：

步骤4.1：通过叶片的瞬态位移求出叶片的瞬态颤振应力，即

δ＝Dε    (5)

其中ε为应变列阵，σ为应力列阵，D为弹性矩阵；

步骤4.2：将叶片随时间变化的位移、颤振应力和压力输出为可被Origin软件读取的格式；

步骤4.3：将转子随时间变化的质量流量和压比输出为可被Origin软件读取的格式；

步骤4.4：将叶片瞬间的位移云图和颤振应力云图输出为可被Tecplot软件读取的格式；

步骤4.5：如果整个计算结束，输出该级转子叶片在压气机特性图上的颤振边界；

模块五：初值计算模块：

初值计算模块通过反复的非线性迭代获得叶片的初始静变形和稳态定常流场，作为以后双向迭代计算的初值，用以加快整个系统的收敛；该模块包含如下步骤：

步骤5.1：首先通过三维模型设计软件UG建立叶片的三维实体模型；然后将实体模型导入网格划分软件HyperMesh中，得到叶片的结构有限元网格；将实体模型导入网格划分软件TurboGrid中，得到叶轮机内流场的流体计算网格；

步骤5.2：调用结构计算模块得到叶片在离心力作用下的初始位移其中K₀是结构的小位移刚度矩阵，f_C为叶片在转动中的离心力列阵；

步骤5.3：调用数据转换模块将叶片变形转化为流体域的边界运动；

步骤5.4：通过流体计算模块调用流体动力学计算软件，进行网格更新，然后计算新的定常流场；

步骤5.5：调用数据转换模块将叶片表面的定常压力转换为结构有限元模型上的节点力；

步骤5.6：调用结构计算模块获得叶片的非线性变形，并对结构的刚度矩阵进行更新；计算叶片非线性变形的具体步骤如下：

步骤5.6.1：根据步骤5.2求得的初始位移a₀，求得应力刚化矩阵K_S、大变形刚度矩阵K_L和反力F^nr；

步骤5.6.2：根据(K₀-K_C+K_σ+K_L)da₁＝F_C-F^nr求出增量位移da₁；

步骤5.6.3：计算总位移a₁＝a₀+da₁；

步骤5.6.4：如果求得的a_n满足收敛条件||da_n||₂＜ε||a_n||₂，其中ε取0.001，则非线性迭代收敛，计算结束，指向步骤5.7；否则用a₁代替a₀，前往步骤5.3；

步骤5.7：求得新的结构刚度矩阵K＝K₀-K_C+K_S+K_L；此时的叶片静变形和定常流场就作为以后双向迭代计算的初始值；

模块六：双向迭代模块：

双向迭代模块反复交替调用结构求解器和流体求解器，同时采用数据转换模块在结构求解器和流体求解器之间传递边界数据，通过流体与结构双向的交互作用，在时间上推进整个由叶片及其周围流场组成的流固耦合系统；双向迭代模块包含下列步骤：

步骤6.1：通过数据转换系统将叶片的变形转换为流体计算域的边界运动；

步骤6.2：通过流体计算模块调用流体动力学计算软件，进行网格更新，然后计算新的非定常流场；

步骤6.3：通过数据转换模块将叶片表面的非定常压力转换为结构有限元模型上的节点力；

步骤6.4：通过结构计算模块求得新的叶片变形；

步骤6.5：对双向迭代的计算时间步进行调整；该流固耦合自适应时间步技术包括下列步骤：

步骤6.5.1：计算

其中u_i为当前时间步计算的节点i在某个自由度的位移，u_i’为上个时间步计算的节点i在某个自由度的位移，n为有限元节点的数目；

步骤6.5.2：取各个自由度的d的平均值为d₀，用以表征当前时间步的收敛速度；

步骤6.5.3：用e代表收敛速度的阈值，用t代表当前的时间步长，如果d₀＜0.5e，则新的时间步长t_new＝2t；如果d₀＞2e，则新的时间步长t_new＝0.5t；如果0.5e＜d₀＜2e，则时间步长不变，t_new＝t；

步骤6.5.4：在获得流固耦合计算的物理时间步长t_new后，对结构动力方程瞬态求解的时间步长和流体非定常求解的时间步长作相应调整；

步骤6.6：通过计算所得的叶片瞬态位移响应，观察颤振的发展进程；如果响应曲线衰减，则不会发生颤振，计算结束；如果响应曲线增长，发生颤振，指向步骤6.1；如果发生颤振且计算达到预先指定时间点，计算结束，指向步骤6.7输出颤振应力；

步骤6.7：通过颤振应力输出模块，将结构的瞬态位移响应转化为瞬态颤振应力，并输出为文件，供后处理软件读取显示。

Images