专利名称：一种基于时域双向迭代的叶轮机叶片颤振应力预测方法
技术领域：
本发明涉及一种计算机辅助分析工具领域的方法，具体是涉及一种基于时域双向 迭代的叶轮机叶片颤振应力预测方法。属于叶轮机械模拟技术领域。
背景技术：
目前，在叶轮机械中，存在着振动叶片与周围流体之间复杂的相互作用，弹性叶片 通过其变形影响流场，流场的扰动反之改变叶片的形状。在这类流固耦合现象中，以颤振对 结构可靠性的影响最大。颤振属于一种自激振动，叶片表面的非定常气动力来源于叶片本 身的运动，非定常气动力和叶片的非定常运动互为因果、循环递增，使得短时间内叶片振幅 急剧增大，甚至造成叶片断裂和机器损坏。由于进行颤振试验难度较高、成本巨大且具有相 当的危险性，所以随着计算机及相关学科的发展，数值模拟方法因其成本低、效率高，成为 颤振分析的重要工具。在叶轮机械的颤振数值模拟中，通常使用能量法和特征值法，近年来也少量开始 使用时间推进法。现有的颤振数值计算方法至少存在以下四个缺点首先，传统的能量法和特征值法都假定叶片以给定形式作运动，将流场和叶片运 动互相孤立，忽略了弹性叶片本身变形对流场的影响，同时也极少考虑离心力对叶片变形 和流场的影响，无法实现气动_结构一体化计算，存在较大的计算误差。其次，由于流体和结构都存在强烈的非线性，这两种领域的非线性交替叠加，形成 了更加复杂的耦合非线性。现有的颤振预测方法都引入了大量线性化假设，甚至只采用结 构向流体的单向传递，弱化甚至消除了耦合非线性，无法准确地描述非线性现象，降低了模 拟精度。再次，现有的方法不能展示整个颤振发展的时间历程，也无法得到颤振应力，而颤 振应力恰恰是颤振现象中非常重要的指标。最后，以往的基于时间推进的叶轮机颤振数值模拟方法都采用固定时间步长，无 法在计算中对步长进行有效的调整，造成计算速度和效率较低，甚至会导致收敛困难；而且 大多采用商用的结构有限元计算软件，在计算中需要反复进行人工干预，不能在软件内部 针对流固耦合计算的特性对整个计算流程和算法进行优化，增加了操作复杂度，降低了计 算效率，使得现有的基于时间推进的预测方法计算费用较高，经济型较差。
发明内容
1、目的本发明的目的是克服常规的叶轮机颤振数值仿真方法的不足，提供一种 基于时域双向迭代的叶轮机叶片颤振应力预测方法，使之能够准确描述流体与叶片之间的 相互作用和耦合非线性现象，并对叶片的颤振应力进行预测。2、技术方案本发明一种基于时域双向迭代的叶轮机叶片颤振应力预测方法，它 是在计算机上采用数值模拟的方法对叶轮机叶片的颤振应力进行预测。该方法具体步骤如 下
步骤一在计算机中设定以下六个模块结构计算模块、流体计算模块、数据转换 模块、颤振应力输出模块、初值计算模块和双向迭代模块。步骤二 调用初值计算模块通过非线性迭代获得叶片的初始静变形和稳态定常流 场作为以后双向迭代计算的初值，用以加快整个系统的收敛。步骤三调用双向迭代模块，在时间上推进由叶片和周围流场组成的流固耦合系 统。步骤四通过颤振应力输出模块，将结构的瞬态位移响应转化为瞬态颤振应力，并 输出为文件，供后处理软件读取显示。其中，步骤一所述的六个模块的具体结构如下模块一结构计算模块结构计算模块的作用是求解三维结构动力方程，获得叶片的瞬态变形。该模块包 含如下步骤步骤1. 1 将叶片表面压力和结构有限元节点单元信息放到输入文件里，并在输 入文件中对叶片的工况进行设置。步骤1. 2 调用Fortran语言编写的结构有限元计算程序求解三维结构动力方程
Md + Cd + (K-Kc + S)a = f + fa + fc(1)其中M，C和K分别是总体质量矩阵、阻尼矩阵和刚度矩阵。K。是旋转软化矩阵。S 是应力刚化矩阵(或称几何刚化矩阵)。f, fa和f。是节点外力列阵、气动力列阵和离心力 列阵。a是待求的有限元节点位移列向量。步骤1. 3 计算结束后，输出叶片表面有限元节点的位移。模块二 流体计算模块在流体计算模块中，提供了各种流体计算软件的接口，可以后台调用流体计算软 件获得每个时间步的三维非定常流场。该模块中根据用户需要调用各种流体动力学计算软 件，如Fluent，CFX等，在计算结束后输出叶片表面压力。模块三数据转换模块数据转换模块包括三个子模块，一是节点配对子模块，作用是找到耦合界面上每 个流体节点所对应的结构单元；二是载荷传递子模块，该子模块通过三维形函数插值，将叶 片表面的定常或非定常压力转换为叶片结构有限元模型上的节点力；三是变形传递子模 块，该子模块将叶片的变形转换为叶片周围流体域的边界运动。下面分别对这三个子模块 进行介绍节点配对子模块包含下列步骤步骤3. 1. 1 对耦合界面上每个流体节点f，遍历所有耦合界面上的结构节点，找 到距离f最近的结构节点S。步骤3. 1. 2 找到所有包含结构节点s的结构单元e2,…，em。步骤3. 1.3 对每个结构单元力(1 = 1，2，…，m)，用拟牛顿法求出f在该单元内
的局部坐标(ri; Si,、)，并求出f与该单元中心的相对距离式=拟+sf+tf。步骤3. 1. 4 找出与f相对距离最小的单元et，即dt = π η^^，.··，《)，该单元 即为流体节点f所对应的结构单元。
变形传递子模块包含下列步骤步骤3. 2. 1 设流体节点f在其所对应结构单元e上的投影点为f’，f’在整体 坐标系下的坐标为(xf，，yf，，zf，)，忽略流体节点与其投影点之间的距离，近似有xf, ^xf, yf, ^ yf，Zf, ^ zf，投影点f’在结构单元e中的局部坐标(r，s, t)通过求解下列方程组得 到 其中η为结构单元e中节点的数目，Xi，yi，Zi(i = 1,2,…，8)为结构单元e中8 个节点在整体坐标系下的坐标。投影点处的形函数NiG = 1,2,…，8)为N1 = (1-s) (1-t) (l-r)/8, N2 = (1+s) (l_t) (1-r) /8N3 = (1+s) (1+t) (1-r)/8, N4 = (1-s) (1+t) (1-r) /8
(3)N5 = (1-s) (1-t) (l+r)/8, N6 = (1+s) (1-t) (1+r) /8N7 = (1+s) (1+t) (1+r)/8, N8 = (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’的位移为uf,，Ui为结构有限元节点i的位移，通过下 列公式求解流体节点f的位移Uf 载荷传递子模块包含下列步骤步骤3. 3. 1 设流体表面网格c的中心点为f。，计算流体压力沿网格c的积分 Ffc = IlPdS^Pfc-Ac，其中Ω为网格c的表面，ρ为耦合界面上的流体压力，&为中心点f。 上的压力，Α。为网格c的表面积。步骤3. 3. 2 设流体表面网格中心点f。对应的结构单元为e，通过计算尺-况/^ (i =1,2,…，8)可以得到结构单元e上8个节点的节点力，其中Ni为&在结构单元e上投 影点处的形函数。步骤3. 3. 3 将结构单元e上8个节点的节点力集合入结构有限元节点载荷列阵。模块四颤振应力输出模块颤振应力输出模块通过叶片的瞬态位移响应计算出叶片的瞬态颤振应力，并输出 各种供颤振分析使用的后处理文件。该模块包括如下步骤
步骤4. 1 通过叶片的瞬态位移求出叶片的瞬态颤振应力，即O=De (5)其中ε为应变列阵，σ为应力列阵，D为弹性矩阵。步骤4. 2 将叶片随时间变化的位移、颤振应力和压力输出为可被Origin软件读 取的格式。步骤4. 3 将转子随时间变化的质量流量和压比输出为可被Origin软件读取的格式。步骤4. 4 将叶片瞬间的位移云图和颤振应力云图输出为可被Tecplot软件读取 的格式。步骤4. 5 如果整个计算结束，输出该级转子叶片在压气机特性图上的颤振边界。模块五初值计算模块初值计算模块通过反复的非线性迭代获得叶片的初始静变形和稳态定常流场，作 为以后双向迭代计算的初值，用以加快整个系统的收敛。该模块包含如下步骤步骤5. 1 首先通过三维模型设计软件UG建立叶片的三维实体模型；然后将实体 模型导入网格划分软件HyperMesh中，得到叶片的结构有限元网格；将实体模型导入网格 划分软件TurboGrid中，得到叶轮机内流场的流体计算网格。步骤5. 2 调用结构计算模块得到叶片在离心力作用下的初始位移其中Ktl 是结构的小位移刚度矩阵，Fc为叶片在转动中的离心力列阵。步骤5. 3 调用数据转换模块将叶片变形转化为流体域的边界运动。步骤5. 4 通过流体计算模块调用流体动力学计算软件，进行网格更新，然后计算 新的定常流场。步骤5. 5 调用数据转换模块将叶片表面的定常压力转换为结构有限元模型上的 节点力。步骤5.6 调用结构计算模块获得叶片的非线性变形，并对结构的刚度矩阵进行 更新。计算叶片非线性变形的具体步骤如下步骤5. 6. 1 根据步骤5. 2求得的初始位移％，求得应力刚化矩阵Ks、大变形刚度 矩阵Kl和反力Fm。步骤5. 6. 2 根据(Kq-Kc+Kq +KL) (Ia1 = Fc-Fnr 求出增量位移(! 。步骤5. 6. 3 计算总位移 B1 = a0+dalo步骤5. 6.4:如果求得的an满足收敛条件I |dan| |2< ε ||an||2，其中ε取0.001， 则非线性迭代收敛，计算结束，指向步骤5. 7 ；否则用ai代替％，前往步骤5. 3。步骤5. 7 求得新的结构刚度矩阵K = Kci-KfKfKp常规的有限元计算软件无法 同时考虑几何非线性和应力刚化的影响，本发明通过步骤5. 7即时对结构刚度矩阵进行修 正，从而克服了这一缺陷。此时的叶片静变形和定常流场就作为以后双向迭代计算的初始 值。模块六双向迭代模块双向迭代模块反复交替调用结构求解器和流体求解器，同时采用数据转换模块在 结构求解器和流体求解器之间传递边界数据，通过流体与结构双向的交互作用，在时间上 推进整个由叶片及其周围流场组成的流固耦合系统。双向迭代模块包含下列步骤
步骤6. 1 通过数据转换系统将叶片的变形转换为流体计算域的边界运动。步骤6. 2 通过流体计算模块调用流体动力学计算软件，进行网格更新，然后计算 新的非定常流场。步骤6. 3 通过数据转换模块将叶片表面的非定常压力转换为结构有限元模型上 的节点力。步骤6. 4 通过结构计算模块求得新的叶片变形。步骤6. 5 对双向迭代的计算时间步进行调整。本发明提供了一种新的流固耦合 自适应时间步技术，通过当前的收敛速度控制下一步计算的时间步长，随时对后续计算的 时间步长进行调整，从而有效地缩短收敛时间，提高计算效率。该流固耦合自适应时间步技 术包括下列步骤步骤6. 5. 1 计算在[(奶)2+(u,-u/ )2+…+Uru1/ )2]/("/+"/+···+"/)，其中
Ui为当前时间步计算的节点i在某个自由度的位移，Ui’为上个时间步计算的节点i在某个 自由度的位移，η为有限元节点的数目。步骤6. 5. 2 取各个自由度的d的平均值为Cltl，用以表征当前时间步的收敛速度。步骤6.5.3 用e代表收敛速度的阈值，用t代表当前的时间步长，如果Cltl <0.5e， 则新的时间步长tnew = 2t ；如果dQ>2e，则新的时间步长tnew = 0. 5t ；如果0. 5e<d0<2e, 则时间步长不变，tnew = t。步骤6. 5. 4 在获得流固耦合计算的物理时间步长后，对结构动力方程瞬态求 解的时间步长和流体非定常求解的时间步长作相应调整。步骤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 调用结构计算模块得到叶片在离心力作用下的初始位移列阵a^AT尺， 其中Ktl是结构的小位移刚度矩阵，Fc为叶片在转动中的离心力列阵。步骤5. 3 调用数据转换模块将叶片变形转化为流体域的边界运动。步骤5. 4 通过流体计算模块调用流体动力学计算软件，进行网格更新，然后计算 新的定常流场。步骤5. 5 调用数据转换模块将叶片表面的定常压力转换为结构有限元模型上的 节点力。步骤5.6 调用结构计算模块获得叶片的非线性变形，并对结构的刚度矩阵进行更新。如果非线性迭代收敛，计算结束，指向步骤5. 7;否则用新的叶片位移列阵取代旧的 叶片位移列阵，前往步骤5. 3。步骤5. 7 求得新的结构刚度矩阵K = Ktl-KfKdKp常规的有限元计算软件无法 同时考虑几何非线性和应力刚化的影响，本发明通过步骤5. 7即时对结构刚度矩阵进行修 正，从而克服了这一缺陷。此时的叶片静变形和定常流场就作为以后双向迭代计算的初始值。模块六双向迭代模块双向迭代模块反复交替调用结构求解器和流体求解器，同时采用数据转换模块在 结构求解器和流体求解器之间传递边界数据，通过流体与结构双向的交互作用，在时间上 推进整个由叶片及其周围流场组成的流固耦合系统。双向迭代模块包含下列步骤步骤6. 1 通过数据转换系统将叶片的变形转换为流体计算域的边界运动。步骤6. 2 通过流体计算模块调用流体动力学计算软件，进行网格更新，然后计算 新的非定常流场。步骤6. 3 通过数据转换模块将叶片表面的非定常压力转换为结构有限元模型上 的节点力。步骤6. 4 通过结构计算模块求得新的叶片变形。步骤6. 5 对双向迭代的计算时间步进行调整。本发明提供了一种新的流固耦合 自适应时间步技术，通过当前的收敛速度控制下一步计算的时间步长，随时对后续计算的 时间步长进行调整，从而有效地缩短收敛时间，提高计算效率。该流固耦合自适应时间步技 术包括下列步骤步骤6. 5. 1 计算oK("厂V+(U,-U/ )、..+ (" -"/ )2]/("/+"/+·..+"/)，其中
Ui为当前时间步计算的节点i在某个自由度的位移，Ui’为上个时间步计算的节点i在某个 自由度的位移，η为有限元节点的数目。步骤6. 5. 2 取各个自由度的d的平均值为Cltl，用以表征当前时间步的收敛速度。步骤6. 5. 3 用e代表收敛速度的阈值，用t代表当前的时间步长，如果Cltl <0. 5e, 则新的时间步长tnew = 2t ；如果dQ>2e，则新的时间步长tnew = 0. 5t ；如果0. 5e<d0<2e, 则时间步长不变，tnew = t。步骤6. 5. 4 在获得流固耦合计算的物理时间步长tnew后，对结构动力方程瞬态求 解的时间步长和流体非定常求解的时间步长作相应调整。步骤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输入数据转换模块，将叶片变形转化为流体域的 边界运动，生成叶片表面流体网格点的位移，写入文件blad印oint. txt ；步骤1.4:通过流体计算模块调用流体动力学计算软件Fluent，读入文件 bladepoint. txt，进行动网格计算，获得更新后的流场网格，并计算新的定常流场，将叶片 表面流体网格点的压力写入文件blad印res. txt ；步骤1. 5 将文件blacbpres. txt输入数据转换模块，将叶片表面的定常压力转换 为结构有限元模型上的节点力，写入文件bladeforce. txt ；步骤1.6 通过结构计算模块调用由标准Fortran语言编写的结构有限元分析程 序，读入节点力文件bladeforce. txt，进行非线性迭代获得叶片的变形和新的结构刚度矩 阵，将叶片表面有限元节点的变形写入文件bladenode. txt，将更新后的结构刚度矩阵写入 t^f牛 restartek. dat ；步骤1. 7 如果非线性迭代收敛，计算结束，所得叶片静变形和定常流场就作为以 后双向迭代计算的初始值；否则如果非线性迭代不收敛，指向步骤1. 3。步骤2:双向迭代调用双向迭代模块，在时间上推进整个流固耦合系统。与常规的单向迭代不同，本 发明通过流体与结构之间的双向迭代实现了气动-结构一体化计算，常规方法与本发明所 采用方法的区别见图5。双向迭代的具体步骤如下步骤2. 1 将文件bladenode. txt输入数据转换模块，采用三维形函数插值方 法，将叶片的变形转换为流体计算域的边界运动，将叶片表面流体网格点的位移写入文件 bladepoint. txt 步骤2. 2:通过流体计算模块调用流体动力学计算软件Fluent，读入文件 bladepoint. txt，进行动网格计算，获得更新后的流场网格，并计算新的三维非定常流场， 将叶片表面流体网格点的压力写入文件blacbpres. txt ；步骤2. 3 将文件blacbpres. 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。通过应力 云图可以获得颤振时叶片应力集中的危险位置，为后续分析提供依据。
权利要求
一种基于时域双向迭代的叶轮机叶片颤振应力预测方法，其特征在于，该方法具体步骤如下步骤一在计算机中设定以下六个模块结构计算模块、流体计算模块、数据转换模块、颤振应力输出模块、初值计算模块和双向迭代模块；步骤二调用初值计算模块通过非线性迭代获得叶片的初始静变形和稳态定常流场作为以后双向迭代计算的初值，用以加快整个系统的收敛；步骤三调用双向迭代模块，在时间上推进由叶片和周围流场组成的流固耦合系统；步骤四通过颤振应力输出模块，将结构的瞬态位移响应转化为瞬态颤振应力，并输出为文件，供后处理软件读取显示。
2.根据权利要求1所述的一种基于时域双向迭代的叶轮机叶片颤振应力预测方法，其 特征在于步骤一中所述的六个模块的具体结构如下模块一结构计算模块结构计算模块的作用是求解三维结构动力方程，获得叶片的瞬态变形；该模块包含如 下步骤步骤1. 1 将叶片表面压力和结构有限元节点单元信息放到输入文件里，并在输入文 件中对叶片的工况进行设置；步骤1. 2 调用Fortran语言编写的结构有限元计算程序求解三维结构动力方程 Md + Ca + (K-Kc + S)a = f + fa + fc⑴其中M，C和K分别是总体质量矩阵、阻尼矩阵和刚度矩阵；K。是旋转软化矩阵；S是应 力刚化矩阵即几何刚化矩阵；f，fa和f。是节点外力列阵、气动力列阵和离心力列阵，a是待 求的有限元节点位移列向量；步骤1. 3 计算结束后，输出叶片表面有限元节点的位移； 模块二 流体计算模块在流体计算模块中，提供了各种流体计算软件的接口，可以后台调用流体计算软件 获得每个时间步的三维非定常流场；该模块中根据用户需要调用流体动力学计算软件 Fluent, CFX，在计算结束后输出叶片表面压力； 模块三数据转换模块数据转换模块包括三个子模块，一是节点配对子模块，作用是找到耦合界面上每个流 体节点所对应的结构单元；二是载荷传递子模块，该子模块通过三维形函数插值，将叶片表 面的定常或非定常压力转换为叶片结构有限元模型上的节点力；三是变形传递子模块，该 子模块将叶片的变形转换为叶片周围流体域的边界运动； 节点配对子模块包含下列步骤步骤3. 1. 1 对耦合界面上每个流体节点f，遍历所有耦合界面上的结构节点，找到距 离f最近的结构节点S;步骤3. 1. 2 找到所有包含结构节点s的结构单元e2,…，effl ；步骤3. 1.3 对每个结构单元&(1 = 1,2,…，m)，用拟牛顿法求出f在该单元内的局部坐标(ri; Si,、)，并求出f与该单元中心的相对距离式=仏2矛；步骤3. 1.4:找出与f相对距离最小的单元et，即dt = min ((IljCl2,…，dm)，该单元即为流体节点f所对应的结构单元； 变形传递子模块包含下列步骤步骤3. 2. 1 设流体节点f在其所对应结构单元e上的投影点为f’，f'在整体坐标系 下的坐标为(Xf,，yf,，zf,)，忽略流体节点与其投影点之间的距离，近似有xf, xf，yf, ^yf, Zf, ^ zf，投影点f’在结构单元e中的局部坐标(r，s, t)通过求解下列方程组得到 ‘ 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’的位移为uf，，ui为结构有限元节点i的位移，通过下列公 式求解流体节点f的位移％ /=1载荷传递子模块包含下列步骤步骤3. 3. 1 设流体表面网格c的中心点为f。，计算流体压力沿网格c的积分Ffc = I1PdS^pfc-Ac，其中Ω为网格c的表面，ρ为耦合界面上的流体压力，&为中心点f。上的压力，Α。为网格c的表面积；步骤3. 3. 2 设流体表面网格中心点f。对应的结构单元为e，通过计算/(i = 1， 2，…，8)可以得到结构单元e上8个节点的节点力，其中Ni为f。在结构单元e上投影点 处的形函数；步骤3. 3. 3 将结构单元e上8个节点的节点力集合入结构有限元节点载荷列阵； 模块四颤振应力输出模块颤振应力输出模块通过叶片的瞬态位移响应计算出叶片的瞬态颤振应力，并输出各种 供颤振分析使用的后处理文件；该模块包括如下步骤步骤4. 1 通过叶片的瞬态位移求出叶片的瞬态颤振应力，即3δ =De (5)其中ε为应变列阵，σ为应力列阵，D为弹性矩阵；步骤4. 2 将叶片随时间变化的位移、颤振应力和压力输出为可被Origin软件读取的 格式；步骤4. 3 将转子随时间变化的质量流量和压比输出为可被Origin软件读取的格式； 步骤4. 4:将叶片瞬间的位移云图和颤振应力云图输出为可被Tecplot软件读取的格式；步骤4. 5 如果整个计算结束，输出该级转子叶片在压气机特性图上的颤振边界； 模块五初值计算模块初值计算模块通过反复的非线性迭代获得叶片的初始静变形和稳态定常流场，作为以 后双向迭代计算的初值，用以加快整个系统的收敛；该模块包含如下步骤步骤5. 1 首先通过三维模型设计软件UG建立叶片的三维实体模型；然后将实体模型 导入网格划分软件HyperMesh中，得到叶片的结构有限元网格；将实体模型导入网格划分 软件TurboGrid中，得到叶轮机内流场的流体计算网格；步骤5. 2 调用结构计算模块得到叶片在离心力作用下的初始位移S^ATFc,其中Ktl是 结构的小位移刚度矩阵，fc为叶片在转动中的离心力列阵；步骤5. 3 调用数据转换模块将叶片变形转化为流体域的边界运动； 步骤5. 4 通过流体计算模块调用流体动力学计算软件，进行网格更新，然后计算新的 定常流场；步骤5. 5 调用数据转换模块将叶片表面的定常压力转换为结构有限元模型上的节点力；步骤5. 6 调用结构计算模块获得叶片的非线性变形，并对结构的刚度矩阵进行更新； 计算叶片非线性变形的具体步骤如下步骤5. 6. 1 根据步骤5. 2求得的初始位移％，求得应力刚化矩阵Ks、大变形刚度矩阵 Kl和反力Fm;步骤 5. 6. 2 根据(KQ-KC+K。+IQda1 = Fc-Fnr 求出增量位移(Ia1 ； 步骤5. 6. 3 计算总位移a: = a0+dai ；步骤5. 6. 4:如果求得的an满足收敛条件I |dan| |2 < ￡||an||2，其中ε取0.001，则 非线性迭代收敛，计算结束，指向步骤5. 7 ；否则用ai代替％，前往步骤5. 3 ；步骤5. 7 求得新的结构刚度矩阵K = K0-KC+KS+KL ；此时的叶片静变形和定常流场就作 为以后双向迭代计算的初始值； 模块六双向迭代模块双向迭代模块反复交替调用结构求解器和流体求解器，同时采用数据转换模块在结构 求解器和流体求解器之间传递边界数据，通过流体与结构双向的交互作用，在时间上推进 整个由叶片及其周围流场组成的流固耦合系统；双向迭代模块包含下列步骤 步骤6. 1 通过数据转换系统将叶片的变形转换为流体计算域的边界运动； 步骤6. 2 通过流体计算模块调用流体动力学计算软件，进行网格更新，然后计算新的 非定常流场；步骤6. 3 通过数据转换模块将叶片表面的非定常压力转换为结构有限元模型上的节点力步骤6. 4 通过结构计算模块求得新的叶片变形；步骤6.5 对双向迭代的计算时间步进行调整；该流固耦合自适应时间步技术包括下 列步骤步骤 6. 5. 1 计算厂的，)2+(uru/ Y+…+(M1TU1/ )2]/(—/+.“+"力，其中 Ui 为当前时间步计算的节点i在某个自由度的位移，u/为上个时间步计算的节点i在某个自由 度的位移，η为有限元节点的数目；步骤6. 5. 2 取各个自由度的d的平均值为Cltl，用以表征当前时间步的收敛速度； 步骤6. 5. 3 用e代表收敛速度的阈值，用t代表当前的时间步长，如果Cltl < 0. 5e，则 新的时间步长tnew = 2t ；如果dQ > 2e，则新的时间步长tnew = 0. 5t ；如果0. 5e < d0 < 2e, 则时间步长不变，tnew = t ；步骤6. 5. 4 在获得流固耦合计算的物理时间步长后，对结构动力方程瞬态求解的 时间步长和流体非定常求解的时间步长作相应调整；步骤6.6 通过计算所得的叶片瞬态位移响应，观察颤振的发展进程；如果响应曲线衰 减，则不会发生颤振，计算结束；如果响应曲线增长，发生颤振，指向步骤6. 1 ；如果发生颤 振且计算达到预先指定时间点，计算结束，指向步骤6. 7输出颤振应力；步骤6. 7 通过颤振应力输出模块，将结构的瞬态位移响应转化为瞬态颤振应力，并输 出为文件，供后处理软件读取显示。
全文摘要
本发明公开了一种基于时域双向迭代的叶轮机叶片颤振应力预测方法，其特征在于，将叶片及其周围流场视为一个三维流固耦合系统，设计了一套时域的双向迭代方法，通过交替求解叶片变形和非定常流场得到叶片的颤振应力。该方法在计算机中设定以下模块结构计算模块、流体计算模块、数据转换模块、颤振应力输出模块、初值计算模块和双向迭代模块。通过非线性迭代获得叶片静变形和稳态流场作为初值；交替调用结构计算模块和流体计算模块在时间上推进整个系统；通过数据转换模块传递流固边界信息；输出时间历程上的颤振应力。本发明实现了叶片与流场的一体化计算，考虑了耦合系统的非线性，能观察整个颤振发展进程，预测叶片的颤振应力。
文档编号G06F17/50GK101908088SQ20101023721
公开日2010年12月8日 申请日期2010年7月22日 优先权日2010年7月22日
发明者徐可宁, 王延荣 申请人:北京航空航天大学