叶轮机械计入叶间相位角的气弹稳定性流固耦合预测方法，采用能量法，通过不同叶间相位角的模态气动阻尼比预测叶轮机械的气弹稳定性。该方法对建立的单扇区有限元模型采用含预应力的模态分析，通过数据传递方法实现流固耦合交界面上的振动位移传递，并采用动网格技术得到考虑叶间相位角的网格文件，用于CFX中Junction Box模块的非定常计算流体力学分析，进而获得各叶间相位角的模态气动阻尼比，用于预测气弹稳定性。该基于流固耦合方法，充分利用了计算结构动力学和计算流体动力学的前沿技术，在每个子系统内保证了计算精度，而固定域和可动域的划分提高了计算效率。它在叶轮机械模拟技术领域具有良好的实用价值和广阔的应用前景。
1.叶轮机械计入叶间相位角的气弹稳定性流固耦合预测方法，其特征在于，该方法包括如下步骤：步骤一：建立单扇区有限元模型；
针对给定的实体模型在有限元软件MARC中建立叶轮机械的单扇区有限元模型；
步骤二：对建立的有限元模型进行含预应力的模态分析；
首先，将MARC中建立的有限元模型导入到有限元软件ANSYS中，定义相应的材料参数，并给定转速以及位移约束条件，通过静力分析得到离心力引起的预应力；
然后，采用振型归一化的模态分析，并计入离心力引起的预应力，获得叶片的固有振动特性，提取叶片振动固有频率以及模态振型，得到固体域叶片表面的节点振幅；
首先，通过给定的叶轮机械叶型数据，在计算流体力学仿真工具CFX的TurboGrid模块中建立叶栅的单通道流场模型，其中叶片表面附近为O型网格，其他区域为H型网格；
然后，将建立的单通道流场模型导入到ICEM CFD中，通过绕旋转轴旋转、复制(N-1)份网格得到全环叶栅模型，并输出流体域网格点和实体单元信息以及流体域各叶片表面网格点和表面单元信息；
步骤四：采用数据传递方法获得插值到流体域叶片表面的振动位移；
首先，基于恒定位置矢量差方法，采用有限单元法中局部坐标形式的形函数进行三维线性插值，将固体域中模态分析得到的第P阶模态对应的节点振动位移插值到流固耦合交界面上的表面单元网格点上，获得流体域中叶片表面网格点的振幅；
然后，通过流体域中各叶片扇区之间的角度关系，即可依次逐个得到全环各叶片表面网格点的振幅；
步骤五：采用动网格技术获得每一个叶片扇区的可动域网格点振幅；
根据流体域中网格点编号顺序，对每一个叶片沿O型域网格的法向逐层递推搜索，得到与该叶片表面网格点编号顺序相同的各层网格点，并按照各层网格点与叶片表面网格点之间法向距离的比值大小来线性插值，得到该叶片可动域的网格点振幅；
步骤六：获得非定常计算流体力学(CFD)分析所需的考虑叶间相位角的网格文件；
首先，将叶片一个振动周期分为K步，叶片以简谐规律振动，给定各叶片相对于第一个叶片的初始相位角，通过初始坐标以及各时间步的振动位移，可以得到一个周期内每一时间步上各叶片扇区可动域的坐标值；
然后，通过固定域坐标值和每一时间步上各叶片扇区可动域的坐标值，可以得到一个周期内每一时间步的流体网格点坐标组成的网格文件；
步骤七：采用CFX中Junction Box模块进行非定常分析；
首先，将步骤三建立的流场模型的网格点坐标值换成步骤六中获得的第一个时间步的网格点坐标值，对流场模型加载边界条件和初始条件；通道入口给定总温和总压，出口给定平均静压，轮缘、轮毂给定无滑移、光滑壁面边界条件，叶片表面给定动网格壁面边界条件；
以叶片表面给定无滑移、光滑壁面边界条件的定常解作为非定常分析的初始条件；
然后，指定叶片的振动模态以及振动幅值，其运动周期为固有频率的倒数，在每一个振动时间步上求解采用k-ε湍流模型封闭的Reynolds平均Navier-Stokes方程，并写出每一个求解时间步的瞬时结果文件；
步骤八：获得各叶间相位角对应的模态气动阻尼比，预测振荡叶栅的气弹稳定性；
首先，在不同叶间相位角情况下，根据一个稳定周期内各时间步的结果文件，基于等效粘性阻尼推导的模态气动阻尼比获得叶片一个振动周期的气动功和模态气动阻尼比；
其中气动功是流场各网格点气动力在一个稳定振动周期内所做的总功；
基于等效粘性阻尼推导的模态气动阻尼比公式为 其中W为非定常气aero
动力在一个周期内所做的气动功，q 为叶片的振动幅值；
然后，通过不同叶间相位角情况下的模态气动阻尼比，判断某阶模态下某一工况的气弹稳定性，在各叶间相位角下，只要有某一个叶间相位角对应的模态气动阻尼比为负，则该阶模态下，该工况气弹失稳。
2.根据权利要求1所述的叶轮机械计入叶间相位角的气弹稳定性流固耦合预测方法，其特征在于：步骤四中所述的数据传递方法的具体实施步骤如下：步骤1：获得与第一个叶片表面的流体网格点F距离最近的固体节点S；
对流体域中第一个叶片表面的任意网格点F，在固体域中找到与其距离最近的固体表面节点S，与节点S相关的单元至多有4个，分别为S1、S2、S3和S4，并从固体域中获得各单元对应的结构单元的表面编号；
步骤2：获得流体网格点F在不同固体单元中的局部坐标值；
流体网格点F的坐标为(x0，y0，z0)，固体单元八个节点的坐标分别为(xi，yi，zi)，(i＝
1，...，8)，每个节点的形函数在局部坐标系下可以表示为Ni(ε，η，ζ)，(i＝1，...，8)，则通过流体网格点坐标和固体单元节点坐标之间的关系：可以得到流体网格点F在与节点S相关的不同单元的局部坐标值；
分别对比不同单元对应的局部坐标值，找出各局部坐标值的平方和最小值对应的单元M以及该单元的表面编号，即为流体网格点F在固体域中的投影单元和单元面，将网格点F投影到单元M的表面得到F′，由于流体域网格和固体域网格对原始几何构型的不同精度导致网格点F和投影点F′之间存在坐标差，采用恒定位置矢量差法，使流体网格点F的函数值始终等于其投影点F′的函数值；
首先，获得流体网格点F在投影单元M上的局部坐标值，并进一步得到单元M的八个形函数Ni(ε，η，ζ)＝(1±ε)(1±η)(1±ζ)，(i＝1，…，8)；
然后，通过固体域中实体单元的八个节点振动位移(ui，vi，wi)(i＝1，...，8)可以插值获得流固耦合交界面上流体网格点F的位移(u0，v0，w0)，具体的位移插值为步骤5：获得其他叶片的表面网格点振动位移；
首先，获得与第一个叶片表面网格点编号相对应的其他各叶片的表面网格点编号；
在旋转坐标轴为z的情况下，各叶片逆时针编号，第i个叶片表面网格点坐标(xi，yi，zi)与第一个叶片表面网格点坐标(x1，y1，z1)之间满足如下关系：其中当第一个叶片在第一象限时， 当第一个叶片在其他象限时，按照对应的角度关系给出θ1；
然后，将第一个叶片表面网格点的振幅分别赋值给其他各叶片的对应网格点，获得全环所有叶片表面网格点的振幅。
3.根据权利要求1所述的叶轮机械计入叶间相位角的气弹稳定性流固耦合预测方法，其特征在于：步骤六中所述生成考虑叶间相位角的网格文件实施步骤如下：步骤(1)：将叶片的一个振动周期分为K步，其中一个周期对应的频率为叶片在某阶模态下的固有频率ω，给定第一个叶片正弦振动的初始相位角为零，则第i个叶片在第k个时间步的相位角为其中，l为节径数，N为全环叶片个数；
其中，Ui，Vi，Wi分别为第i个叶片扇区的某一个网格点在x，y，z三个方向的振幅；
步骤(2)：通过动网格技术搜索获得第i个叶片对应可动域的网格点编号，并从全域网格点坐标值可以得到该叶片扇区的可动域网格点初始坐标值(xi，yi，zi)，而该叶片扇区的每一个网格点在第k个时间步的x，y，z三个方向振动位移分别为uik，vik，wik，则该叶片扇区的每一个网格点在第k个时间步的坐标值为步骤(3)：通过固定域坐标值和每一时间步上各叶片扇区可动域的坐标值，可以得到一个周期内每一时间步的流体网格点坐标，并将每一个时间步的网格点坐标值输出为网格文件。
叶轮机械计入叶间相位角的气弹稳定性流固耦合预测方法
(一)技术领域
[0001] 本发明属于叶轮机械模拟技术领域，具体涉及基于能量法的一种叶轮机械计入叶间相位角的气弹稳定性流固耦合预测方法。(二)背景技术
[0002] 科学技术的进步使得叶轮机械的叶尖切线速度和压比等不断提高，导致风扇/压气机叶片气弹稳定性问题突显。随着叶轮机械向着大功率、高性能的发展，风扇/压气机叶片的工作条件越来越恶劣，叶尖切线速度越来越高，叶片的刚性相对越来越小，叶片发生气弹不稳定的可能性越来越大。国内外许多发动机均发生过风扇/压气机叶片的气弹稳定性故障，叶轮机械叶片气弹稳定性问题成为现代航空发动机设计所面临的主要难点之一。
[0003] 20世纪80年代以来，随着计算机的发展和计算技术的进步，出现了数值方法，而计算流体力学的发展更进一步推动了振荡叶栅气弹稳定性问题的研究。目前，研究最多的数值方法主要为能量法和特征值法。由于叶轮机械内部的流场呈现着复杂的三维非定常特性，而且分离流动、旋转失速、激波、激波与附面层的相互干扰、激波与激波间的相互干扰等非线性因素的影响使得气动载荷很难用准确的函数来表达，通过特征值法来定量地准确描述气弹稳定性几乎不可能实现，因此，采用能量法能够很好地保证流场计算精度。
[0004] 能量法忽略了流体和固体之间的相互作用而仅仅考虑固体振动对流场的影响，因此，分别对流体域和固体域单独求解，并通过发展的基于有限元形函数的概念实现流体域和固体域之间的数据传递，这也是气弹稳定性数值预测技术的关键。
[0005] 另一方面，由于叶片的振荡会引起叶栅之间非定常流场的相互作用，对叶片气弹稳定性有很大的影响。早期对叶轮机械的气动弹性问题大多假设叶栅中的所有叶片均以相同的振幅振荡，并且不计叶间相位角的影响，但是随着对气动弹性研究的深入，发现叶间相位角是影响叶片气动弹性的一个相当关键的参数，它对叶片气弹稳定性起着重要的作用，而且实际叶片并不可能满足上述的假设。
[0006] 因此，对于叶轮机械叶片气动弹性这样的流固耦合问题，有必要发展和完善适用于工程设计需求的气弹稳定性数值预测方法，并考虑叶间相位角等关键参数的影响。其中最大的困难是如何准确高效地实现流固耦合交界面上的数据传递，并在振荡叶栅中计入叶间相位角。目前还没有一项较为成熟的技术来预测计入叶间相位角影响的叶轮机械叶片的气弹稳定性。(三)发明内容
[0007] 本发明的目的是提供叶轮机械计入叶间相位角的气弹稳定性流固耦合预测方法，它解决了现有技术的不足。
[0008] 本发明选择精度和效率平衡的能量法作为模型，它满足叶片小幅振动的假设，并且忽略流体对振动系统的影响，而仅仅考虑固体域对流体域的影响，能够降低计算成本。
[0009] 本发明叶轮机械计入叶间相位角的气弹稳定性流固耦合预测方法，其具体步骤如下：
[0010] 步骤一：建立单扇区有限元模型；
[0011] 针对给定的实体模型在有限元软件MARC中建立叶轮机械的单扇区有限元模型；
[0012] 步骤二：对建立的有限元模型进行含预应力的模态分析；
[0013] 首先，将MARC中建立的有限元模型导入到有限元软件ANSYS中，定义相应的材料参数，并给定转速以及位移约束条件，通过静力分析得到离心力引起的预应力；
[0014] 然后，采用振型归一化的模态分析，并计入离心力引起的预应力，获得叶片的固有振动特性，提取叶片振动固有频率以及模态振型，得到固体域叶片表面的节点振幅；
[0015] 步骤三：建立计算流体力学模型；
[0016] 首先，通过给定的叶轮机械叶型数据，在计算流体力学仿真工具CFX的TurboGrid模块中建立叶栅的单通道流场模型，其中叶片表面附近为O型网格，其他区域为H型网格；
[0017] 然后，将建立的单扇区流场模型导入到ICEM CFD中，通过绕旋转轴旋转、复制(N-1)份网格得到全环叶栅模型，并输出流体域网格点和实体单元信息以及流体域各叶片表面网格点和表面单元信息；
[0018] 步骤四：采用数据传递方法获得插值到流体域叶片表面的振动位移；
[0019] 首先，基于恒定位置矢量差方法，采用有限单元法中局部坐标形式的形函数进行三维线性插值，将固体域中模态分析得到的第P阶模态对应的节点振动位移插值到流固耦合交界面上的表面单元网格点上，获得流体域中叶片表面网格点的振幅；
[0020] 然后，通过流体域中各叶片扇区之间的角度关系，即可依次逐个得到全环各叶片表面网格点的振幅；
[0021] 步骤五：采用动网格技术获得每一个叶片扇区的可动域网格点振幅；
[0022] 根据流体域中网格点编号顺序，对每一个叶片沿O型域网格的法向逐层递推搜索，得到与该叶片表面网格点编号顺序相同的各层网格点，并按照各层网格点与叶片表面网格点之间法向距离的比值大小来线性插值，得到该叶片可动域的网格点振幅；
[0023] 步骤六：获得非定常计算流体力学(CFD)分析所需的考虑叶间相位角的网格文件；
[0024] 首先，将叶片一个振动周期分为K步，叶片以简谐规律振动，给定各叶片相对于第一个叶片的初始相位角，通过初始坐标以及各时间步的振动位移，可以得到一个周期内每一时间步上各叶片扇区可动域的坐标值；
[0025] 然后，通过固定域坐标值和每一时间步上各叶片扇区可动域的坐标值，可以得到一个周期内每一时间步的流体网格点坐标组成的网格文件；
[0026] 步骤七：采用CFX中Junction Box模块进行非定常分析；
[0027] 流场分析通过计算流体力学仿真工具CFX实现；
[0028] 首先，将步骤三建立的流场模型的网格点坐标值换成步骤六中获得的第一个时间步的网格点坐标值，对流场模型加载边界条件和初始条件；通道入口给定总温和总压，出口给定平均静压，轮缘、轮毂给定无滑移、光滑壁面边界条件，叶片表面给定动网格壁面边界条件；以叶片表面给定无滑移、光滑壁面边界条件的定常解作为非定常分析的初始条件；
[0029] 然后，指定叶片的振动模态以及振动幅值，其运动周期为固有频率的倒数，在每一个振动时间步上求解采用k-ε湍流模型封闭的Reynolds平均Navier-Stokes方程，并写出每一个求解时间步的瞬时结果文件；
[0030] 步骤八：获得各叶间相位角对应的模态气动阻尼比，预测振荡叶栅的气弹稳定性；
[0031] 首先，在不同叶间相位角情况下，根据一个稳定周期内各时间步的结果文件，基于等效粘性阻尼推导的模态气动阻尼比获得叶片一个振动周期的气动功和模态气动阻尼比；
[0032] 其中气动功是流场各网格点气动力在一个稳定振动周期内所做的总功；
[0033] 基于等效粘性阻尼推导的模态气动阻尼比公式为 其中W为非定aero
常气动力在一个周期内所做的气动功，q 为叶片的振动幅值；
[0034] 然后，通过不同叶间相位角情况下的模态气动阻尼比，判断某阶模态下某一工况的气弹稳定性，在各叶间相位角下，只要有某一个叶间相位角对应的模态气动阻尼比为负，则该阶模态下，该工况气弹失稳。
[0035] 这种基于能量法的模拟方法适用于叶轮机械全环叶栅考虑叶间相位角影响的气弹稳定性分析，在保证了一定计算精度的基础上，能够很好地提高计算效率，节约计算资源和减少时间成本。这种方法中所用到的基于有限元形函数的数据传递方法准确并且效率高。而为了保证计算效率，将流体域分为固定域和每个叶片表面附近的可动域，并在每一个可动域采用动网格技术，摒弃CFX本身只能运动表面网格节点的约束，使得叶片能够在更大的振幅下振动，同时通过节径数给出全环叶栅各叶片的叶间相位角，并通过不同叶间相位角下的模态气动阻尼比判断叶轮机械叶片的气弹稳定性。
[0036] 下面分别介绍该数值预测方法中用到的数据传递方法和考虑叶间相位角的网格文件生成方法。
[0037] 数据传递方法的具体实施步骤如下：
[0038] 步骤一：获得与第一个叶片表面的流体网格点F距离最近的固体节点S；
[0039] 对流体域中第一个叶片表面的任意网格点F，在固体域中找到与其距离最近的固体表面节点S，与节点S相关的单元至多有4个，分别为S1、S2、S3和S4，并从固体域中获得各单元对应的结构单元的表面编号；
[0040] 步骤二：获得流体网格点F在不同固体单元中的局部坐标值；
[0041] 流体网格点F的坐标为(x0，y0，z0)，固体单元八个节点的坐标分别为(xi，yi，zi)，(i＝1，...，8)，每个节点的形函数在局部坐标系下可以表示为Ni(ε，η，ζ)，(i＝1，..，8)，则通过流体网格点坐标和固体单元节点坐标之间的关系：
[0042]
[0043] 可以得到流体网格点F在与节点S相关的不同单元的局部坐标值；
[0044] 其中局部坐标系下的形函数可以表示为
[0045] N(ε，η，ζ)＝(1±ε)(1±η)(1±ζ)
[0046] 步骤三：获得流体网格点F对应固体域的投影单元M；
[0047] 分别对比不同单元对应的局部坐标值，找出各局部坐标值的平方和最小值对应的单元M以及该单元的表面编号，即为流体网格点F在固体域中的投影单元和单元面，将网格点F投影到单元M的表面得到F′，由于流体域网格和固体域网格对原始几何构型的不同精度导致网格点F和投影点F′之间存在坐标差，采用恒定位置矢量差法，使流体网格点F的函数值始终等于其投影点F′的函数值；
[0048] 步骤四：插值得到流体网格点F的位移；
[0049] 首先，获得流体网格点F在投影单元M上的局部坐标值，并进一步得到单元M的八个形函数Ni(ε，η，ζ)＝(1±ε)(1±η)(1±ζ)，(i＝1，...，8)；
[0050] 然后，通过固体域中实体单元的八个节点振动位移(ui，vi，wi)(i＝1，...，8)可以插值获得流固耦合交界面上流体网格点F的位移(u0，v0，w0)，具体的位移插值为[0051]
[0052] 步骤五：获得其他叶片的表面网格点振动位移；
[0053] 首先，获得与第一个叶片表面网格点编号相对应的其他各叶片的表面网格点编号；
[0054] 在旋转坐标轴为z的情况下，各叶片逆时针编号，第i个叶片表面网格点坐标(xi，yi，zi)与第一个叶片表面网格点坐标(x1，y1，z1)之间满足如下关系：
[0055]
[0056] 其中当第一个叶片在第一象限时， 当第一个叶片在其他象限时，按照对应的角度关系给出θ1；
[0057] 其中， N为全环叶片个数；
[0058] 然后，将第一个叶片表面网格点的振幅分别赋值给其他各叶片的对应网格点，获得全环所有叶片表面网格点的振幅。
[0059] 考虑叶间相位角的网格文件生成方法具体实施步骤如下：
[0060] 步骤一：将叶片的一个振动周期分为K步，其中一个周期对应的频率为叶片在某阶模态下的固有频率ω，给定第一个叶片正弦振动的初始相位角为零，则第i个叶片在第k个时间步的相位角为
[0061]
[0062] 其中，l为节径数，N为全环叶片个数；
[0063] 则该叶片扇区的每一个网格点在该时间步的振动位移为
[0064]
[0065] 其中，Ui，Vi，Wi分别为第i个叶片扇区的某一个网格点在x，y，z三个方向的振幅；
[0066] 步骤二：通过动网格技术搜索获得第i个叶片对应可动域的网格点编号，并从全域网格点坐标值可以得到该叶片扇区的可动域网格点初始坐标值(xi，yi，zi)，而该叶片扇区的每一个网格点在第k个时间步的x，y，z三个方向振动位移分别为uik，vik，wik，则该叶片扇区的每一个网格点在第k个时间步的坐标值为
[0067]
[0068] 步骤三：通过固定域坐标值和每一时间步上各叶片扇区可动域的坐标值，可以得到一个周期内每一时间步的流体网格点坐标，并将每一个时间步的网格点坐标值输出为网格文件。
[0069] 本发明的优点在于：
[0070] (1)基于能量法的流固耦合预测方法充分借助当前计算结构动力学和计算流体动力学的前沿技术来实现各部分的单独计算，并通过设计的数据传递方法来进行流体域和固体域的数据传递，从而使得耦合计算大大简化，并在每个子系统内均能保持计算的高精度；
[0071] (2)基于恒定位置矢量差方法，采用有限元法中局部坐标形式的形函数进行三维线性插值，方法简单，虽然会在流体网格造成一定的横向剪切变形，但是计算效率较高；
[0072] (3)将全环叶栅模型分为叶片表面附近的可动域和其他区域的固定域能够很好地保证计算效率，而在通过节径数给出的叶间相位角情况下计算全范围内的模态气动阻尼比能够更准确地用于判断叶轮机械叶片的气弹稳定性。(四)附图说明
[0073] 图1是本发明计入叶间相位角的气弹稳定性流固耦合预测方法流程图；
[0074] 图2是现有MARC中的NASA67叶片有限元模型；
[0075] 图3是现有NASA67叶片第一阶弯曲模态；
[0076] 图4是现有NASA67叶片全环叶栅流场模型；
[0077] 图5是本发明数据传递方法流程图；
[0078] 图6是现有NASA67叶片单扇区叶栅的可动域和固定域；
[0079] 图7是本发明动网格技术示意图；
[0080] 图8是本发明模态气动阻尼比与叶间相位角的关系图；(五)具体实施方式
[0081] 下面将结合附图和实施例对本发明作进一步的详细说明。本发明是叶轮机械计入叶间相位角的气弹稳定性流固耦合预测方法，方法的流程如图1所示。
[0082] 步骤一：建立单扇区有限元模型；
[0083] 针对给定的实体模型在有限元软件MARC中建立叶轮机械的单扇区有限元模型；
[0084] 步骤二：对建立的有限元模型进行含预应力的模态分析；
[0085] 首先，将MARC中建立的有限元模型导入到有限元软件ANSYS中，定义相应的材料参数，并给定转速以及位移约束条件，通过静力分析得到离心力引起的预应力；
[0086] 然后，采用振型归一化的模态分析，并计入离心力引起的预应力，获得叶片的固有振动特性，提取叶片振动固有频率以及模态振型，得到固体域叶片表面的节点振幅；
[0087] 步骤三：建立计算流体力学模型；
[0088] 首先，通过给定的叶轮机械叶型数据，在计算流体力学仿真工具CFX的TurboGrid模块中建立叶栅的单通道流场模型，其中叶片表面附近为O型网格，其他区域为H型网格；
[0089] 然后，将建立的单扇区流场模型导入到ICEM CFD中，通过绕旋转轴旋转、复制(N-1)份网格得到全环叶栅模型，并输出流体域网格点和实体单元信息以及流体域各叶片表面网格点和表面单元信息；
[0090] 步骤四：采用数据传递方法获得插值到流体域叶片表面的振动位移；
[0091] 首先，基于恒定位置矢量差方法，采用有限单元法中局部坐标形式的形函数进行三维线性插值，将固体域中模态分析得到的第P阶模态对应的节点振动位移插值到流固耦合交界面上的表面单元网格点上，获得流体域中叶片表面网格点的振幅；
[0092] 然后，通过流体域中各叶片扇区之间的角度关系，即可依次逐个得到全环各叶片表面网格点的振幅；
[0093] 步骤五：采用动网格技术获得每一个叶片扇区的可动域网格点振幅；
[0094] 根据流体域中网格点编号顺序，对每一个叶片沿O型域网格的法向逐层递推搜索，得到与该叶片表面网格点编号顺序相同的各层网格点，并按照各层网格点与叶片表面网格点之间法向距离的比值大小来线性插值，得到该叶片可动域的网格点振幅；
[0095] 步骤六：获得非定常计算流体力学(CFD)分析所需的考虑叶间相位角的网格文件；
[0096] 首先，将叶片一个振动周期分为K步，叶片以简谐规律振动，给定各叶片相对于第一个叶片的初始相位角，通过初始坐标以及各时间步的振动位移，可以得到一个周期内每一时间步上各叶片扇区可动域的坐标值；
[0097] 然后，通过固定域坐标值和每一时间步上各叶片扇区可动域的坐标值，可以得到一个周期内每一时间步的流体网格点坐标组成的网格文件；
[0098] 步骤七：采用CFX中Junction Box模块进行非定常分析；
[0099] 流场分析通过计算流体力学仿真工具CFX实现；
[0100] 首先，将步骤三建立的流场模型的网格点坐标值换成步骤六中获得的第一个时间步的网格点坐标值，对流场模型加载边界条件和初始条件；通道入口给定总温和总压，出口给定平均静压，轮缘、轮毂给定无滑移、光滑壁面边界条件，叶片表面给定动网格壁面边界条件；以叶片表面给定无滑移、光滑壁面边界条件的定常解作为非定常分析的初始条件；
[0101] 然后，指定叶片的振动模态以及振动幅值，其运动周期为固有频率的倒数，在每一个振动时间步上求解采用k-ε湍流模型封闭的Reynolds平均Navier-Stokes方程，并写出每一个求解时间步的瞬时结果文件；
[0102] 步骤八：获得各叶间相位角对应的模态气动阻尼比，预测振荡叶栅的气弹稳定性；
[0103] 首先，在不同叶间相位角情况下，根据一个稳定周期内各时间步的结果文件，基于等效粘性阻尼推导的模态气动阻尼比获得叶片一个振动周期的气动功和模态气动阻尼比；
[0104] 其中气动功是流场各网格点气动力在一个稳定振动周期内所做的总功；
[0105] 基于等效粘性阻尼推导的模态气动阻尼比公式为 其中W为非定aero
常气动力在一个周期内所做的气动功，q 为叶片的振动幅值；
[0106] 然后，通过不同叶间相位角情况下的模态气动阻尼比，判断某阶模态下某一工况的气弹稳定性，在各叶间相位角下，只要有某一个叶间相位角对应的模态气动阻尼比为负，则该阶模态下，该工况气弹失稳。
[0107] 这种基于能量法的模拟方法适用于叶轮机械全环叶栅考虑叶间相位角影响的气弹稳定性分析，在保证了一定计算精度的基础上，能够很好地提高计算效率，节约计算资源和减少时间成本。这种方法中所用到的基于有限元形函数的数据传递方法准确且高效。而为了保证计算效率，将流体域分为固定域和叶片表面附近的可动域，并在每一个可动域保证叶片能够在更大的振幅下振动，同时通过节径数给出全环叶栅各叶片的叶间相位角，并通过不同叶间相位角下的模态气动阻尼比判断叶轮机械叶片的气弹稳定性。
[0108] 下面通过一个实例来说明此数值方法。选用NASA67转子叶片，其基本的设计气动参数见表1。基于给定的叶型数据在有限元软件MARC中建立单个叶片的有限元模型如图2，其中单元选取八节点六面体实体单元。
[0109] 表1NASA 67转子叶片基本设计参数
[0110]
[0111] 将MARC中建立的有限元模型导入到有限元软件ANSYS中进行含预应力的模态分析。其中材料选取钛合金，边界条件给定叶根固支。首先，通过静力分析得到离心力引起的预应力，然后，进行含预应力的模态分析，得到叶片的各阶动频和模态，如图3为叶片的第一阶弯曲模态。
[0112] 通过给定的叶型数据在CFX TurboGrid中建立流体域单通道模型，并通过旋转复制得到全环叶栅模型，如图4，其中，模型忽略了叶尖间隙的作用。
[0113] 采用恒定位置矢量差法，基于有限元局部坐标系下形函数的概念发展了数据传递方法，将单个叶片的第一阶弯曲模态插值到流体域的流固耦合交界面的叶片表面网格点上，其数据传递的流程如图5，具体实施步骤如下：
[0114] 步骤一：获得与第一个叶片表面的流体网格点F距离最近的固体节点S；
[0115] 对流体域中第一个叶片表面的任意网格点F，在固体域中找到与其距离最近的固体表面节点S，与节点S相关的单元至多有4个，分别为S1、S2、S3和S4，并从固体域中获得各单元对应的结构单元的表面编号；
[0116] 步骤二：获得流体网格点F在不同固体单元中的局部坐标值；
[0117] 流体网格点F的坐标为(x0，y0，z0)，固体单元八个节点的坐标分别为(xi，yi，zi)，(i＝1，...，8)，每个节点的形函数在局部坐标系下可以表示为Ni(ε，η，ζ)，(i＝1，..，8)，则通过流体网格点坐标和固体单元节点坐标之间的关系：
[0118]
[0119] 可以得到流体网格点F在与节点S相关的不同单元的局部坐标值；
[0120] 其中局部坐标系下的形函数可以表示为
[0121] N(ε，η，ζ)＝(1±ε)(1±η)(1±ζ)
[0122] 步骤三：获得流体网格点F对应固体域的投影单元M；
[0123] 分别对比不同单元对应的局部坐标值，找出各局部坐标值的平方和最小值对应的单元M以及该单元的表面编号，即为流体网格点F在固体域中的投影单元和单元面，将网格点F投影到单元M的表面得到F′，由于流体域网格和固体域网格对原始几何构型的不同精度导致网格点F和投影点F′之间存在坐标差，采用恒定位置矢量差法，使流体网格点F的函数值始终等于其投影点F′的函数值；
[0124] 步骤四：插值得到流体网格点F的位移；
[0125] 首先，获得流体网格点F在投影单元M上的局部坐标值，并进一步得到单元M的八个形函数Ni(ε，η，ζ)＝(1±ε)(1±η)(1±ζ)，(i＝1，...，8)；
[0126] 然后，通过固体域中实体单元的八个节点振动位移(ui，vi，wi)(i＝1，...，8)可以插值获得流固耦合交界面上流体网格点F的位移(u0，v0，w0)，具体的位移插值为[0127]
[0128] 步骤五：获得其他叶片的表面网格点振动位移；
[0129] 首先，获得与第一个叶片表面网格点编号相对应的其他各叶片的表面网格点编号；
[0130] 在旋转坐标轴为z的情况下，各叶片逆时针编号，第i个叶片表面网格点坐标(xi，yi，zi)与第一个叶片表面网格点坐标(x1，y1，z1)之间满足如下关系：
[0131]
[0132] 其中当第一个叶片在第一象限时， 当第一个叶片在其他象限时，按照对应的角度关系给出θ1；
[0133] 其中， N为全环叶片个数；
[0134] 然后，将第一个叶片表面网格点的振幅分别赋值给其他各叶片的对应网格点，获得全环所有叶片表面网格点的振幅。
[0135] 在振荡叶栅作用下的非定常流场分析时，仅将振荡叶片附近的流体域指定为可动域，用O型网格来实现，O型可动域的网格点坐标在每一求解时刻均随着叶片的振荡实时更新，远离叶片的其他区域均被定义为固定域，网格点坐标始终保持不变，用H型网格来实现。在振荡叶栅作用下的流场分析中，将每一个叶栅通道的流体域网格划分为叶片表面附近的O型网格和其他区域的H型网格如图6。
[0136] 从网格点的编号来看，每个叶片扇区的O型可动域的各层网格点编号规律相同，而且对应网格点位于同一直线上，如图7，当叶片振动时，流体域中固体边界附近的网格点坐标均会按法向距离的大小比例变化，且始终能够保证流体域表面一层的网格尺寸大小满足精度要求。当获得每个叶片附近可动域的振幅后，生成考虑叶间相位角的网格文件，生成方法的具体实施步骤如下：
[0137] 步骤一：将叶片的一个振动周期分为K步，其中一个周期对应的频率为叶片在某阶模态下的固有频率ω，给定第一个叶片正弦振动的初始相位角为零，则第i个叶片在第k个时间步的相位角为
[0138]
[0139] 其中，l为节径数，N为全环叶片个数；
[0140] 则该叶片扇区的每一个网格点在该时间步的振动位移为
[0141]
[0142] 其中，Ui，Vi，Wi分别为第i个叶片扇区的某一个网格点在x，y，z三个方向的振幅；
[0143] 步骤二：通过动网格技术搜索获得第i个叶片对应可动域的网格点编号，并从全域网格点坐标值可以得到该叶片扇区的可动域网格点初始坐标值(xi，yi，zi)，而该叶片扇区的每一个网格点在第k个时间步的x，y，z三个方向振动位移分别为uik，vik，wik，则该叶片扇区的每一个网格点在第k个时间步的坐标值为
[0144]
[0145] 步骤三：通过固定域坐标值和每一时间步上各叶片扇区可动域的坐标值，可以得到一个周期内每一时间步的流体网格点坐标，并将每一个时间步的网格点坐标值输出为网格文件。
[0146] 得到计入叶间相位角影响的各时间步网格文件后，在CFX的Junction Box动网格模块中进行非定常分析，其中以叶片不振荡、考虑叶间相位角的全环叶栅模型在设计状态的边界条件下计算得到的定常解作为初始条件，并设置设计状态的非定常边界条件，在每一个求解时间步读取网格文件，实现叶片的振荡作用。
[0147] CFX非定常分析后，读取一个稳定周期内各时间步的叶片表面气动力以及相应的振荡位移，进而求得一个周期内的非定常气动功，并得到给定叶间相位角情况下基于等效粘性阻尼推导的模态气动阻尼比，进而曲线拟合得到不同叶间相位角对应的模态气动阻尼比，如图8。图中模态气动阻尼比与叶间相位角近似呈简谐规律，且在接近0°时模态气动阻尼比达到最小值，在接近180°时模态气动阻尼比达到最大值，而且在整个叶间相位角范围内，模态气动阻尼比均为正，表示NASA67在叶片的一阶弯曲模态下不会出现气弹失稳的现象，这与实验值是吻合的。
[0148] 对于其他的叶轮机械叶片的气弹稳定性判断，可以通过以上的方法，在叶片的某一阶模态下，得到不同叶间相位角情况下的模态气动阻尼比，只有当整个叶间相位角范围内模态气动阻尼比均为正时，叶轮机械叶片才气弹稳定。
