CN102799730B - 一种燃气轮机风扇叶片反扭过程的预估方法 - Google Patents

Abstract

本发明提供一种燃气轮机风扇叶片反扭过程的预估方法，目的是为了解决目前燃气轮机风扇叶片反扭设计方法不完善的问题，该方法对热态叶型进行ANSYS离心力静态变形分析，预测冷态叶型，再对其进行ANSYS离心力静态变形分析，过渡虚热态叶型，修正固体和流体网格，气动力变形的流固耦合分析等步骤最终得到冷态叶型。本发明为实现各个气动工况下预测的热态叶型与叶片实际工作态叶型的一致性，对精确预估燃气轮机风扇叶片反扭过程的方法进行研究是十分重要的，尤其对于发展大涵道比的跨音速风扇有着重要的工程应用价值。

Description

一种燃气轮机风扇叶片反扭过程的预估方法

技术领域

本发明涉及叶片叶型设计的技术领域，具体涉及一种燃气轮机风扇叶片反扭过程的预估方法。

背景技术

燃气轮机，特别是涡扇发动机风扇级的性能对发动机总性能有着至关重要的影响，而叶型则是决定风扇性能的最重要的参数之一，所以保证风扇叶片在工作点的形状和设计点的相同或接近是实现其设计性能指标的关键环节。燃气轮机风扇在工作过程中，叶片在气动力、离心力和热应力的共同作用下，会产生伸长、展开等复杂的变形，叶片会从加工的形状(冷态叶型)变为工作转速下的叶型(热态叶型)，而热态叶型又受转速和负载的制约，所以不同气动工况下热态叶型略有不同。在航空发动机的叶片设计中，一般根据气动性能、强度、寿命等要求设计出在设计工况下的叶片最佳气动几何形状(热叶片)，再利用结构分析对设计叶型进行反推以求得冷态加工叶型(热到冷)，即叶片反扭设计。另一方面，许多涉及到叶片变形方面的研究，如气动弹性的强迫响应、颤振等问题，往往需要叶片在非设计气动工况下的精确叶型。以冷态加工叶型为分析起点，得到工作热态下叶片形状的预估称为“冷到热”的过程。

大量研究表明，燃气轮机叶片形状的细微差别会影响叶尖间隙、气流角等，从而导致风扇叶片通道内的流动和设计工况有一定差别。特别对跨音速风扇，叶片变形可能会引起激波位置的移动，从而偏离其最佳设计点，进而影响燃气轮机的性能、效率，甚至带来气动弹性稳定性等问题。而现代发动机高推重比的需求，要求风扇叶片高负载，低重量。新型轻质复合材料、空心结构等措施的采用使得叶片相对刚性下降，在大的气动载荷下，较之常规风扇更易产生较大变形，气动力的非线性效应更明显。而在超、跨音速流动工况下叶片变形对激波位置有重要影响。叶片在气动力作用下的变形作用体现的更突出。而传统叶片反扭设计，只是根据结构强度需求，叶片安装角偏转一个角度。这种方法对新结构、新材料的设计很难得到优化的叶片制造叶型，从而导致叶片在设计转速下工作时与设计叶型不一致。

无论是“冷到热”还是“热到冷”的叶型计算，都要求反推的冷态叶型具有较高精度，才能实现各个气动工况下预测的热态叶型与叶片实际工作态叶型的一致性。因此对精确预估燃气轮机风扇叶片反扭过程的方法进行研究是十分重要的，尤其对于发展大涵道比的跨音速风扇有着重要的工程应用价值。

发明内容

发明目的是为了解决目前燃气轮机风扇叶片反扭设计方法不完善的问题，进而提供了一种燃气轮机风扇叶片反扭过程的预估方法。

本发明为解决上述技术问题采取的技术方案是：一种燃气轮机风扇叶片反扭过程的预估方法，该方法由以下步骤实现的：

步骤一：以设计时叶片几何数据为原始热态叶型，叶片材料为设计时选取的材料，获取弹性常数，泊松比性能参数，以设计点的气动工况为工作点；

步骤二：根据原始设计工况下热态叶片几何型值，构建风扇叶片的几何模型；

步骤三：将风扇叶片进行网格划分，并结合设计时给定的材料性能参数，构建风扇叶片的有限元模型；

步骤四：结合风扇叶片的设计点运行工况，即转速的设计要求，叶片位移约束边界条件，进行叶片在离心力作用下的静态变形分析，获得离心负载下叶片的静变形量；

步骤五：使用原始热态叶型减去离心静变形得到固体域中的预测冷态叶型①：首先提取叶片各节点的几何坐标即Xs0，Ys0，Zs0，然后根据离心静变形计算结果，在X、Y、Z轴方向上提取叶片各节点处变形值即UX1，UY1，UZ1，在步骤三中的风扇叶片几何外形的基础上实行叶片初步反扭计算，在X、Y、Z轴方向上使各节点坐标取值为Xs0-UX1，Ys0-UY1，Zs0-UZ1，从而，确定风扇叶片的预测冷态几何形式；

步骤六：将燃气轮机风扇叶片的几何模型导入计算流体力学分析软件的前处理器AUTOGRID，建立流体域，完成流体网格的划分；

步骤七：使用原始热态叶型减去离心静变形得到流体域中的预测冷态叶型①：首先提取叶片各节点的几何坐标即Xf0，Yf0，Zf0，然后使用所述有限元模型得到节点上的静变形量UX1，UY1，UZ1插值到叶片的流体网格节点上以获得UXf1，UYf1，UZf1；然后，使用动网格技术，将基于原始叶型下的流体域网格节点以叶片固面上节点变形量为边界条件，设流体域其他边界变形为0，进行调整，并最大程度保证流体网格质量不变；

步骤八：使用获得的冷态叶型，加载设计转速，材料性能参数，及叶片支撑边界条件构建新的风扇叶片有限元模型；进行叶片在离心力作用下的静态变形分析，获得离心负载下叶片的静变形量；在X、Y、Z轴方向上提取叶片各节点处变形值即UX2，UY2，UZ2；完成自由振动模态分析，求解动力学方程，其动力学方程如下：

$M\stackrel{\·\·}{U}+C\stackrel{\·}{U}+\mathrm{KU}=0---\left(1\right)$

其中M为质量矩阵，C为阻尼矩阵，K为刚度矩阵，为加速度，为速度，U为位移矢量；提取风扇叶片多阶固有频率，以及各阶频率下有限元各节点的振动模态值；

步骤九：在流体分析模型中，使用获得的冷态叶型加上离心静变形得到一个过渡虚热态叶型：首先提取叶片各节点的几何坐标即Xf，Yf，Zf，然后使用所述叶片有限元模型节点上的静变形量UX2，UY2，UZ2插值到叶片的流体网格节点上以获得UXf2，UYf2，UZf2；然后，使用动网格技术，以叶片固面上节点为移动边界，流体域其他边界变形为0，并最大程度保证流体网格质量不变，调整基于冷态叶型下的流体域网格节点；

步骤十：结合风扇叶片的运行工况进行叶片流固耦合的瞬态响应分析，获得气动力作用下的叶片变形量；建立基于RANS方程的流体动力学非定常流动计算模型，如下式：

其中Ω、是某一物理时刻的控制体的体积、表面面积；是守恒型流体状态向量，为对流通量，为粘性通量；是源项矢量，代表了外力和坐标系运动对守恒律的贡献；湍流流动的模拟采用了Sparlart-Allmas一方程湍流模型；求解RANS方程，从而获取流动速度及作用在风扇叶片上的气动力；

将有限元模型上节点的模态值插值到叶片的流体网格点上；组集动力方程所需的风扇叶片固有模态集；使用获得的固有模态集解耦并求解动力学方程，得到广义坐标下N个线性无关的常微分方程组如下，其中N为叶片模态数，

$M_g\left\{\stackrel{\·\·}{\mathrm{\ζ}}\right\}+C_g\left\{\stackrel{\·}{\mathrm{\ζ}}\right\}+K_g\left\{\mathrm{\ζ}\right\}=\left\{Q_g\right\}---\left(12\right)$

其中{ζ}为广义位移向量，是广义速度向量，是广义加速度向量，M_g、C_g和K_g分别是归一化后的广义质量矩阵、广义阻尼矩阵和广义刚度矩阵，{Q_g}为广义气动力，采用多步龙格库塔积分法分别求解各阶模态方程(12)，得到各个模态的广义位移，并线性叠加得到叶片的总物理位移；

叶片在非定常气动力作用下的变形使用流固耦合方法计算：在每个物理时间步，耦合求解流体运动方程和叶片变形的控制方程，并互相作为边界条件进行物理量交换；首先流体域求解了时间精确的非定常雷诺平均N-S方程，以得到每一步由于叶片变形而引起的流场变化，计算气动力作用下叶片表面受到的广义力，并求解广义坐标的结构动力学方程，得到广义位移{ζ}，再利用模态叠加法得到叶片表面每一点的位移{x}；采用投影-形函数插值法处理流体域和固体域的信息交换，然后以叶片外其他边界静止或滑动为约束条件，采用动网格技术，移动流场中所有网格点，以满足给定的叶片变形条件；以变化后网格为初值，迭代求解N-S方程，得到下一时刻的流场；进行以上循环达到收敛性目标为止；通过模态叠加的结构非线性瞬态响应分析，得到稳态响应的值，即为非定常气动力引起的变形，此刻叶型为虚热态叶型；

步骤十一：获取的虚热态叶型，判断其是否与原始热态叶型一致；即虚热态叶型与原始热态叶型之间的叶型差是否达到收敛要求，如果满足，则完成风扇叶片从热态叶型到冷态加工叶型的反扭设计；

步骤十二：如果不满足叶型差收敛条件，则需要在步骤五预测的冷态加工叶型基础上实施叶型修正：首先提取步骤十一所述的虚热态叶片各节点的几何坐标即X，Y，Z，然后提取步骤六中原始叶型各流体节点的几何坐标即Xf0，Yf0，Zf0，获得叶片各节点叶型差即UX3、UY3、UZ3；在步骤八的有限元模型中提取叶片各节点的几何坐标即Xs，Ys，Zs，将叶片流体网格各节点上的叶型差即UX3、UY3、UZ3插值到叶片有限元模型各节点上以获得UXs3，UYs3，UZs3；在步骤八中的风扇叶片几何外形的基础上实行叶型修正，在X、Y、Z轴方向上使各节点坐标取值为Xs+UXs3，Ys+UYs3，Zs+UZs3，从而，确定新的风扇叶片几何形式，即预测冷态叶型②；

步骤十三：在风扇叶片预测冷态叶型②的基础上，重复步骤八至步骤十二的过程，对新的风扇叶片几何外形实行叶型修正，直至满足最大位移值小于1×10^-2mm的叶型差收敛要求，完成风扇叶片从热态到冷态叶型的反扭设计，获取适合风扇叶片设计工况的冷态加工叶型。

进一步的，步骤二中使用三维实体构型软件UG构建风扇叶片的几何模型。

进一步的，步骤三中所述的将风扇叶片进行网格划分，并结合设计时给定的材料性能参数具体为：将风扇叶片的几何模型导入有限元软件ANSYS中，选取单元类型SOLID45对其进行网格划分，并结合设计时给定的材料性能参数，构建燃气轮机风扇叶片的有限元模型。

进一步的，步骤四中所述的进行叶片在离心力作用下的静态变形分析，获得离心负载下叶片的静变形量具体为：使用有限元软件ANSYS进行叶片在离心力作用下的静态变形分析，获得离心负载下叶片的静变形量。

进一步的，步骤一中的叶片材料为钛合金，密度为4440kg/m³，弹性常数为EX＝1.12E11Pa，泊松比为0.3。

进一步的，步骤十中所述的“流体域求解了时间精确的非定常雷诺平均N-S方程”可以使用通用的流体商业软件；所述的“使用获得的固有模态集解耦并求解动力学方程”为现有技术；所述的“采用投影-形函数插值法处理流体域和固体域的信息交换”可以采用现有插值技术；所述的“采用动网格技术，移动流场中所有网格点，以满足给定的叶片变形条件”可以采用现有动网格技术。

本发明与现有技术相比的优点在于：

(1)、本发明提出的叶片反扭过程的预估方法中同时考虑了离心力和非定常气动力对叶片变形的影响；通过耦合求解流体动力学方程和结构动力学方程，在叶片变形分析中考虑了气动力随叶型变化的非线性作用，而现有技术由于气动负载和叶片结构变形之间无耦合作用，不能包含气动力对叶片变形的非线性效应。

(2)、风扇叶片在超、跨音速流动工况下运行时，气动力的非线性效应更明显，本发明提出的叶片反扭过程的预估方法可以有效提高此类工况下叶片反扭过程的计算精度；对于新结构(如空心结构)、新材料(如新型轻质复合材料)的叶片设计，此方法与现有技术相比易于得到优化的叶片制造叶型。

(3)、本发明提出的叶片反扭过程的预估方法是一种普遍方法，不只局限于风扇叶片，同样适用于燃气轮机的压气机叶片，以及涡轮叶片的反扭过程计算。

附图说明

图1是本发明的一种燃气轮机风扇叶片反扭过程的预估方法流程图；

图2是本发明的投影-形函数插值流程图；

图3是本发明的流固耦合迭代流程图。

具体实施方式

下面结合附图和具体实施方式进一步说明本发明。

如图1～3所示，本实施方式的一种燃气轮机风扇叶片反扭过程的预估方法是由以下步骤实现的：

步骤一：以设计时叶片几何数据为原始热态叶型，叶片材料为设计时选取的材料，获取弹性常数，泊松比性能参数。如钛合金：弹性常数为EX＝1.12E11Pa，泊松比0.3，密度4440kg/m3。以设计点的气动工况为工作点；

步骤二：根据原始设计工况下热态叶片几何型值，使用三维实体构型软件UG构建风扇叶片的几何模型；

步骤三：将风扇叶片的几何模型导入有限元软件ANSYS中，选取单元类型SOLID45对其进行网格划分，并结合设计时给定的材料性能参数，构建风扇叶片的有限元模型；

步骤四：结合风扇叶片的设计点运行工况，即转速的设计要求，叶片位移约束边界条件，使用有限元软件ANSYS进行叶片在离心力作用下的静态变形分析，获得离心负载下叶片的静变形量；

步骤五：在固体分析模型中，使用原始热态叶型减去离心静变形得到固体域中的预测冷态叶型①：首先在ANSYS软件中提取叶片各节点的几何坐标即Xs0，Ys0，Zs0，然后根据离心静变形计算结果，在X、Y、Z轴方向上提取叶片各节点处变形值即UX1，UY1，UZ1，在步骤三中的风扇叶片几何外形的基础上实行叶片初步反扭计算，在X、Y、Z轴方向上使各节点坐标取值为Xs0-UX1，Ys0-UY1，Zs0-UZ1，从而，确定风扇叶片的预测冷态几何形式；

步骤六：将燃气轮机风扇叶片的几何模型导入计算流体力学分析软件的前处理器AUTOGRID，建立流体域，完成流体网格的划分；

步骤七：使用原始热态叶型减去离心静变形得到流体域中的预测冷态叶型①：首先在流体分析模型中提取叶片各节点的几何坐标即Xf0，Yf0，Zf0，然后使用投影-形函数插值程序将叶片有限元模型节点上的静变形量UX1，UY1，UZ1，插值到叶片的流体网格节点上，获得UXf1，UYf1，UZf1。然后，使用动网格技术，将基于原始叶型下的流体域网格节点以叶片固面上节点变形量为边界条件，设流体域其他边界变形为0，进行调整，并最大程度保证流体网格质量不变；

步骤八：在有限元分析软件中，使用获得的冷态叶型，加载设计转速，材料性能参数，及叶片支撑边界条件构建新的风扇叶片有限元模型。使用ANSYS软件进行叶片在离心力作用下的静态变形分析，获得离心负载下叶片的静变形量。在X、Y、Z轴方向上提取叶片各节点处变形值即UX2，UY2，UZ2。利用ANSYS软件完成自由振动模态分析，求解动力学方程，其动力学方程如下：

$M\stackrel{\·\·}{U}+C\stackrel{\·}{U}+\mathrm{KU}=0---\left(1\right)$

步骤九：在流体分析模型中，使用获得的冷态叶型加上离心静变形得到一个过渡虚热态叶型：首先在流体分析模型中提取叶片各节点的几何坐标即Xf，Yf，Zf，然后使用投影-形函数插值程序将叶片有限元模型节点上的静变形量UX2，UY2，UZ2，插值到叶片的流体网格节点上，获得UXf2，UYf2，UZf2。然后，使用动网格技术，以叶片固面上节点为移动边界，流体域其他边界变形为0，并最大程度保证流体网格质量不变，调整基于冷态叶型下的流体域网格节点；

步骤十：结合风扇叶片的运行工况，即风扇级的进口气流总压，总温，出口静压，固壁面气流速度等设计要求，进行叶片流固耦合的瞬态响应分析，获得气动力作用下的叶片变形量。建立基于RANS方程的流体动力学非定常流动计算模型，如公式(2)：

求解RANS方程，从而获取流动速度及作用在风扇叶片上的气动力。

推导公式(2)的过程如下：在大多数流固耦合问题中，由于流体域的计算网格随时间运动和变形，所以必须求解任意拉格朗日-欧拉系(Arbitrary Lagrangian Eulerian,ALE)下的流体运动方程。由质量连续方程、动量守恒方程和能量守恒方程组成的非定常Navier-Stokes方程，写为矢量积分形式，如公式(2)。

其中Ω、是某一物理时刻的控制体的体积、表面面积。是守恒型流体状态向量，为对流通量；是源项矢量，代表了外力和坐标系运动对守恒律的贡献。

$\begin{array}{ccc}\stackrel{\→}{Q}=\left[\begin{array}{c}\mathrm{\ρ}\\ \mathrm{\ρu}\\ \mathrm{\ρv}\\ \mathrm{\ρw}\\ \mathrm{\ρE}\end{array}\right]& {\stackrel{\→}{F}}_c=\left[\begin{array}{c}\mathrm{\ρV}\\ \mathrm{\ρuV}+n_{x}p\\ \mathrm{\ρvV}+n_{y}p\\ \mathrm{\ρwV}+n_{z}p\\ \mathrm{\ρHV}\end{array}\right]+V_t\left[\begin{array}{c}\mathrm{\ρ}\\ \mathrm{\ρu}\\ \mathrm{\ρv}\\ \mathrm{\ρw}\\ \mathrm{\ρE}\end{array}\right]& \stackrel{\→}{S}=\left[\begin{array}{c}0\\ {\mathrm{\ρf}}_{e,x}\\ {\mathrm{\ρf}}_{e,y}\\ {\mathrm{\ρf}}_{e,z}\\ \mathrm{\ρ}{\stackrel{\→}{f}}_e\·\stackrel{\→}{v}+{\stackrel{\·}{q}}_h\end{array}\right]\end{array}---\left(3\right)$

V_t为控制体表面网格运动速度，

$V_t=\stackrel{\→}{g}\·\stackrel{\→}{n}=n_x\frac{\∂x}{\∂t}+n_y\frac{\∂y}{\∂t}+n_z\frac{\∂z}{\∂t}---\left(4\right)$

V为控制体表面的逆变速度，

$V\≡\stackrel{\→}{v}\·\stackrel{\→}{n}=n_{x}u+n_{y}v+n_{z}w---\left(5\right)$

为粘性通量：

${\stackrel{\→}{F}}_v=\left[\begin{array}{c}0\\ n_x{\mathrm{\τ}}_{\mathrm{xx}}+n_y{\mathrm{\τ}}_{\mathrm{xy}}+n_z{\mathrm{\τ}}_{\mathrm{xz}}\\ n_x{\mathrm{\τ}}_{\mathrm{yx}}+n_y{\mathrm{\τ}}_{\mathrm{yy}}+n_z{\mathrm{\τ}}_{\mathrm{yz}}\\ n_x{\mathrm{\τ}}_{\mathrm{zx}}+n_y{\mathrm{\τ}}_{\mathrm{zy}}+n_z{\mathrm{\τ}}_{\mathrm{zz}}\\ n_x{\mathrm{\θ}}_x+n_y{\mathrm{\θ}}_y+n_z{\mathrm{\θ}}_z\end{array}\right]---\left(6\right)$

其中，

${\mathrm{\θ}}_x={\mathrm{u\τ}}_{\mathrm{xx}}+{\mathrm{v\τ}}_{\mathrm{xy}}+{\mathrm{w\τ}}_{\mathrm{xz}}+k\frac{\∂T}{\∂x}$

${\mathrm{\θ}}_y={\mathrm{u\τ}}_{\mathrm{yx}}+{\mathrm{v\τ}}_{\mathrm{yy}}+{\mathrm{w\τ}}_{\mathrm{yz}}+k\frac{\∂T}{\∂y}---\left(7\right)$

${\mathrm{\θ}}_z={\mathrm{u\τ}}_{\mathrm{zx}}+{\mathrm{v\τ}}_{\mathrm{zy}}+{\mathrm{w\τ}}_{\mathrm{zz}}+k\frac{\∂T}{\∂z}$

对于牛顿流体，切应力项和速度梯度相关，应用Stokes假设，各切应力项可写为：

${\mathrm{\τ}}_{\mathrm{xx}}=2{\mathrm{\μu}}_x-\frac{2}{3}\mathrm{\μ}(u_x+v_y+w_z)$

${\mathrm{\τ}}_{\mathrm{yy}}={2\mathrm{\μv}}_y-\frac{2}{3}\mathrm{\μ}(u_x+v_y+w_z)$

${\mathrm{\τ}}_{\mathrm{zz}}=2{\mathrm{\μw}}_z-\frac{2}{3}\mathrm{\μ}(u_x{+v}_y+w_z)---\left(8\right)$

τ_xy＝τ_yx＝μ(u_y+v_x)

τ_xz＝τ_zx＝μ(u_z+w_x)

τ_yz＝τ_zy＝μ(v_z+w_y)

这里k为导热系数。对于工程湍流问题，根据Boussinesq提出的涡粘性假设，以涡粘性系数考虑湍流效应。动力粘性系数μ和导热系数k可以写成下列形式：

μ＝μ_l+μ_t

$k=\mathrm{Cp}\·(\frac{{\mathrm{\μ}}_l}{{\mathrm{Pr}}_l}+\frac{{\mathrm{\μ}}_t}{{\mathrm{Pr}}_t})---\left(9\right)$

这里下标l代表层流流动，t代表湍流流动。湍流粘性系数μ_t由涡粘性变量确定：

${\mathrm{\μ}}_t=f_{v1}\mathrm{\ρ}\tilde{v}---\left(10\right)$

建立Spalart-Allmaras模式输运方程，涡粘性变量表述如下：

$\frac{\∂\tilde{v}}{\∂t}+\frac{\∂\left(\tilde{v}{v}_j\right)}{\∂t}=C_{b1}(1-f_{t2})\tilde{S}\tilde{v}+\frac{1}{\mathrm{\σ}}\left\{\frac{\∂}{{\∂x}_j}\right[(v_L+\tilde{v})\frac{\∂\tilde{v}}{{\∂x}_j}]+C_{b2}\frac{\∂\tilde{v}}{{\∂x}_j}\frac{\∂\tilde{v}}{\∂x_j}\}-[C_{\mathrm{\ω}1}{f}_{\mathrm{\ω}}-\frac{C_{b1}}{{\mathrm{\κ}}^2}{f}_{t2}]{\left(\frac{\tilde{v}}{d}\right)}^2+f_{t1}{\left|\right|\mathrm{\Δ}\stackrel{\→}{v}\left|\right|}_2^2---\left(11\right)$

根据以上Boussinesq涡粘性假设及相应的一方程输运方程，RANS方程就可以完成封闭求解，从而获取流动速度及作用在风扇叶片上的气动力；

使用投影-形函数插值程序将有限元模型上节点的模态值插值到叶片的流体网格点上。组集动力方程所需的风扇叶片固有模态集。使用获得的固有模态集解耦并求解动力学方程，得到广义坐标下N个线性无关的常微分方程组如下，其中N为叶片模态数：

$M_g\left\{\stackrel{\·\·}{\mathrm{\ζ}}\right\}+C_g\{\stackrel{\·}{\mathrm{\ζ}}\}+K_g\left\{\mathrm{\ζ}\right\}=\left\{Q_g\right\}---\left(12\right)$

采用多步龙格库塔积分法分别求解各阶模态方程(12)，得到各个模态的广义位移，并线性叠加得到叶片的总物理位移。公式(12)推导过程如下：

基于结构动力学建立的叶片动力学方程为：

$M\left\{\stackrel{\·\·}{x}\right\}+C\left\{\stackrel{\·}{x}\right\}+K\left\{x\right\}=Q---\left(13\right)$

其中{x}是振动系统的物理位移向量，是速度向量，是加速度向量。M为系统的质量矩阵，C为阻尼矩阵，K为刚度矩阵，Q为叶片受到的外力。在燃气轮机中流体的密度远小于叶片的密度，叶片自身的机械阻尼相对于流体的气动阻尼来说要小的多，一般可以认为叶片的固有频率和振型不受气动力的影响。定义系统的广义坐标与物理坐标之间的关系式(14)，带入方程(13)，得到解耦为N个线性无关的常微分方程组(12)。其中，Φ为质量归一化的模态，{Q_g}为广义气动力，如公式(15)所示。

{ζ}=Φ{x}  (14)

叶片在非定常气动力作用下的变形使用流固耦合方法计算：在每个物理时间步，耦合求解流体运动方程和叶片变形的控制方程，并互相作为边界条件进行物理量交换。首先流体域求解了时间精确的非定常雷诺平均N-S方程，以得到每一步由于叶片变形而引起的流场变化，计算气动力作用下叶片表面受到的广义力，并求解广义坐标的结构动力学方程，得到广义位移{ζ}，再利用模态叠加法得到叶片表面每一点的位移{x}。采用投影-形函数插值法处理流体域和固体域的信息交换，然后以叶片外其他边界静止或滑动为约束条件，采用动网格技术，移动流场中所有网格点，以满足给定的叶片变形条件；以变化后网格为初值，迭代求解N-S方程，得到下一时刻的流场；进行以上循环达到收敛性目标为止。通过模态叠加的结构非线性瞬态响应分析，得到稳态响应的值，即为非定常气动力引起的变形。此刻叶型为虚热态叶型。

步骤十一：根据流固耦合迭代计算获取的虚热态叶型，判断其是否与原始热态叶型一致。即虚热态叶型与原始热态叶型之间的叶型差是否达到收敛要求。如果满足，则完成风扇叶片从热态叶型到冷态加工叶型的反扭设计；

步骤十二：如果不满足叶型差收敛条件，则需要在步骤五预测的冷态加工叶型基础上实施叶型修正：首先提取流固耦合计算终止后叶片各节点的几何坐标即X，Y，Z，然后提取步骤六中原始叶型各流体节点的几何坐标即Xf0，Yf0，Zf0，获得叶片各节点叶型差即UX3、UY3、UZ3。在步骤八的有限元模型中提取叶片各节点的几何坐标即Xs，Ys，Zs，通过投影-形函数插值程序将叶片流体网格各节点上的叶型差即UX3、UY3、UZ3，插值到叶片有限元模型各节点上，获得UXs3，UYs3，UZs3。在步骤八中的风扇叶片几何外形的基础上实行叶型修正，在X、Y、Z轴方向上使各节点坐标取值为Xs+UXs3，Ys+UYs3，Zs+UZs3，从而，确定新的风扇叶片几何形式，即预测冷态叶型②；

步骤十三：在风扇叶片预测冷态叶型②的基础上，重复步骤八至步骤十二的过程，对新的风扇叶片几何外形实行叶型修正，直至满足最大位移值小于1×10^-2mm的叶型差收敛要求，完成风扇叶片从热态到冷态叶型的反扭设计，获取适合风扇叶片设计工况的冷态加工叶型。

本发明未详细阐述的技术内容属于本领域技术人员的公知技术。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (6)

1.一种燃气轮机风扇叶片反扭过程的预估方法，其特征在于：该方法由以下步骤实现的：

步骤一：以设计时叶片几何数据为原始热态叶型，叶片材料为设计时选取的材料，获取弹性常数，泊松比性能参数，以设计点的气动工况为工作点；

步骤二：根据原始设计工况下热态叶片几何型值，构建风扇叶片的几何模型；

步骤三：将风扇叶片进行网格划分，并结合设计时给定的材料性能参数，构建风扇叶片的原始有限元模型；

步骤四：结合风扇叶片的设计点运行工况，即转速的设计要求，叶片位移约束边界条件，进行叶片在离心力作用下的静态变形分析，获得离心负载下叶片的静变形量；

步骤五：使用原始热态叶型减去离心静变形得到固体域中的预测冷态叶型①：首先提取叶片各节点的几何坐标即Xs0，Ys0，Zs0，然后根据离心静变形计算结果，在X、Y、Z轴方向上提取叶片各节点处变形值即UX1，UY1，UZ1，在步骤三中的风扇叶片的有限元模型的基础上实行叶片初步反扭计算，在X、Y、Z轴方向上使各节点坐标取值为Xs0-UX1，Ys0-UY1，Zs0-UZ1，从而，确定风扇叶片的预测冷态几何形式；

步骤六：将燃气轮机风扇叶片的几何模型导入计算流体力学分析软件的前处理器AUTOGRID，建立流体域，完成流体网格的划分；

步骤七：使用原始热态叶型减去离心静变形得到流体域中的预测冷态叶型①：首先提取叶片各节点的几何坐标即Xf0，Yf0，Zf0，然后使用所述有限元模型得到节点上的静变形量UX1，UY1，UZ1插值到叶片的流体网格节点上以获得UXf1，UYf1，UZf1；然后，使用动网格技术，将基于原始叶型下的流体域网格节点以叶片固面上节点变形量为边界条件，设流体域其他边界变形为0，进行调整，并最大程度保证流体网格质量不变；

步骤八：使用获得的冷态叶型，加载设计转速，材料性能参数，及叶片支撑边界条件构建新的风扇叶片有限元模型；进行叶片在离心力作用下的静态变形分析，获得离心负载下叶片的静变形量；在X、Y、Z轴方向上提取叶片各节点处静变形值即UX2，UY2，UZ2；完成自由振动模态分析，求解动力学方程，其动力学方程如下：

$M\stackrel{\·\·}{U}+C\stackrel{\·}{U}+\mathrm{KU}=0---\left(1\right)$

其中M为质量矩阵，C为阻尼矩阵，K为刚度矩阵，为加速度，为速度，U为位移矢量；提取风扇叶片多阶固有频率，以及各阶频率下有限元各节点的振动模态值；

步骤九：在流体分析模型中，使用获得的冷态叶型加上步骤八中获得的离心静变形得到一个过渡虚热态叶型：首先提取叶片各节点的几何坐标即Xf，Yf，Zf，然后使用步骤八中所述叶片有限元模型节点上的静变形量UX2，UY2，UZ2插值到叶片的流体网格节点上以获得UXf2，UYf2，UZf2；然后，使用动网格技术，以叶片固面上节点为移动边界，流体域其他边界变形为0，并最大程度保证流体网格质量不变，调整基于冷态叶型下的流体域网格节点；

步骤十：结合风扇叶片的运行工况进行叶片流固耦合的瞬态响应分析，获得气动力作用下的叶片变形量；建立基于RANS方程的流体动力学非定常流动计算模型，如下式：

其中Ω、是某一物理时刻的控制体的体积、表面面积；是守恒型流体状态向量，为对流通量，为粘性通量；是源项矢量，代表了外力和坐标系运动对守恒律的贡献；湍流流动的模拟采用了Sparlart-Allmas一方程湍流模型；求解RANS方程，从而获取流动速度及作用在风扇叶片上的气动力；

将步骤八中获得的有限元模型上节点的模态值插值到叶片的流体网格点上；组集动力方程所需的风扇叶片固有模态集；使用获得的固有模态集解耦并求解动力学方程，得到广义坐标下N个线性无关的常微分方程组如下，其中N为叶片模态数，

$M_g\left\{\stackrel{\·\·}{\mathrm{\ζ}}\right\}+C_g\left\{\stackrel{\·}{\mathrm{\ζ}}\right\}+K_g\left\{\mathrm{\ζ}\right\}=\left\{Q_g\right\}---\left(12\right)$

其中{ζ}为广义位移向量，是广义速度向量，是广义加速度向量，M_g、C_g和K_g分别是归一化后的广义质量矩阵、广义阻尼矩阵和广义刚度矩阵，{Q_g}为广义气动力，采用多步龙格库塔积分法分别求解各阶模态方程(12)，得到各个模态的广义位移，并线性叠加得到叶片的总物理位移；

叶片在非定常气动力作用下的变形使用流固耦合方法计算：在每个物理时间步，耦合求解流体运动方程和叶片变形的控制方程，并互相作为边界条件进行物理量交换；首先流体域求解了时间精确的非定常雷诺平均N-S方程，以得到每一步由于叶片变形而引起的流场变化，计算气动力作用下叶片表面受到的广义力，并求解广义坐标的结构动力学方程，得到广义位移{ζ}，再利用模态叠加法得到叶片表面每一点的位移{x}；采用投影-形函数插值法处理流体域和固体域的信息交换，然后以叶片外其他边界静止或滑动为约束条件，采用动网格技术，移动流场中所有网格点，以满足给定的叶片变形条件；以变化后网格为初值，迭代求解N-S方程，得到下一时刻的流场；进行以上循环达到收敛性目标为止；通过模态叠加的结构非线性瞬态响应分析，得到稳态响应的值，即为非定常气动力引起的变形，此刻叶型为虚热态叶型；

步骤十一：获取的虚热态叶型，判断其是否与原始热态叶型一致；即虚热态叶型与原始热态叶型之间的叶型差是否达到收敛要求，如果满足，则完成风扇叶片从热态叶型到冷态加工叶型的反扭设计；

步骤十二：如果不满足叶型差收敛条件，则需要在步骤五预测的冷态加工叶型基础上实施叶型修正：首先提取步骤十一所述的虚热态叶片各节点的几何坐标即X，Y，Z，然后提取步骤六中原始叶型各流体节点的几何坐标即Xf0，Yf0，Zf0，获得叶片各节点叶型差即UX3、UY3、UZ3；在步骤八的有限元模型中提取叶片各节点的几何坐标即Xs，Ys，Zs，将叶片流体网格各节点上的叶型差即UX3、UY3、UZ3插值到叶片有限元模型各节点上以获得UXs3，UYs3，UZs3；在步骤八中的新的风扇叶片有限元模型的基础上实行叶型修正，在X、Y、Z轴方向上使各节点坐标取值为Xs+UXs3，Ys+UYs3，Zs+UZs3，从而，确定新的风扇叶片几何形式，即预测冷态叶型②；

步骤十三：在风扇叶片预测冷态叶型②的基础上，重复步骤八至步骤十二的过程，对新的风扇叶片几何外形实行叶型修正，直至满足最大位移值小于1×10^-2mm的叶型差收敛要求，完成风扇叶片从热态到冷态叶型的反扭设计，获取适合风扇叶片设计工况的冷态加工叶型。

2.根据权利要求1所述的一种燃气轮机风扇叶片反扭过程的预估方法，其特征在于：步骤二中使用三维实体构型软件UG构建风扇叶片的几何模型。

3.根据权利要求1所述的一种燃气轮机风扇叶片反扭过程的预估方法，其特征在于：步骤三中所述的将风扇叶片进行网格划分，并结合设计时给定的材料性能参数具体为：将风扇叶片的几何模型导入有限元软件ANSYS中，选取单元类型SOLID45对其进行网格划分，并结合设计时给定的材料性能参数，构建燃气轮机风扇叶片的有限元模型。

4.根据权利要求1所述的一种燃气轮机风扇叶片反扭过程的预估方法，其特征在于：步骤四中所述的进行叶片在离心力作用下的静态变形分析，获得离心负载下叶片的静变形量具体为：使用有限元软件ANSYS进行叶片在离心力作用下的静态变形分析，获得离心负载下叶片的静变形量。

5.根据权利要求1所述的一种燃气轮机风扇叶片反扭过程的预估方法，其特征在于：步骤一中的叶片材料为钛合金，密度为4440kg/m³，弹性常数为EX＝1.12E11Pa，泊松比为0.3。

6.根据权利要求1所述的一种燃气轮机风扇叶片反扭过程的预估方法，其特征在于：步骤十中所述的“流体域求解了时间精确的非定常雷诺平均N-S方程”可以使用通用的流体商业软件。