CN105808829B

CN105808829B - 一种基于cpu+gpu异构并行计算的透平机械叶片固有频率特性分析方法

Info

Publication number: CN105808829B
Application number: CN201610117959.2A
Authority: CN
Inventors: 谢永慧; 刘天源; 张荻
Original assignee: Xian Jiaotong University
Current assignee: Xian Jiaotong University
Priority date: 2016-03-02
Filing date: 2016-03-02
Publication date: 2018-10-30
Anticipated expiration: 2036-03-02
Also published as: CN105808829A

Abstract

本发明提供了一种基于CPU+GPU异构并行计算的透平机械叶片固有频率特性分析方法，先按照待分析的透平机械叶片的三维模型和材料参数建立有限元模型，对叶片进行预应力分析；其次进行叶片网格数据预处理，在CPU和GPU上同时计算单元的刚度矩阵和质量矩阵，并组装为叶片的总刚度和质量矩阵；再对叶片和轮缘的约束条件进行设置，包括边界的刚性和弹性位移约束，叶根‑轮缘或连接件的接触耦合，修正总刚度矩阵；然后用CPU+GPU异构并行算法提取总刚度和质量矩阵的广义特征值和特征向量；之后将特征值和特征向量转化为叶片的频率以及振型并输出；最后判断固有振型的振动类型，据此绘制叶片的频率曲线分布，振动安全图或坎贝尔图。

Description

一种基于CPU+GPU异构并行计算的透平机械叶片固有频率特性分析方法

技术领域

本发明属于与工程设计与计算领域，具体涉及一种基于CPU+GPU异构并行计算的透平机械叶片固有频率特性分析方法。

背景技术

近年来，由于GPU强大的浮点计算能力并且以超过CPU摩尔定律的速度高速发展，使得基于GPU的并行算法迅速成为高性能计算的领域的研究热点之一。工程上常用的传统有限元并行计算方法，在技术层面上主要采用分布式计算、并行机或多线程等并行处理技术，或是在单个计算节点上利用GPU对线性方程组等通用问题的求解进行局部加速，对于工程领域的很多具体问题，加速效果非常有限。

目前尚未有专门针对透平机械叶片频率特性分析的软件，传统的通用有限元软件对于透平机械叶片这种复杂部件分析来说操作太过繁琐，而且随着工程设计对计算精度要求的提高，以及整圈失谐叶片研究的深入化，不得不使用更多的节点和单元(千万级别)去模拟实际叶片，通用有限元软件的计算速度已经不适于处理如此庞大的计算量。

发明内容

本发明的目的在于提供一种基于CPU+GPU异构并行计算的透平机械叶片固有频率特性分析方法，该方法专门针对透平机械叶片固有频率分析设计，且支持CPU+GPU异构并行的工程计算，能够方便工程设计人员操作，而且可以大幅度提升计算速度。

为了实现上述目的，本发明采用以下技术方案。

一种基于CPU+GPU异构并行计算的透平机械叶片固有频率特性分析方法，包括以下步骤：

1)叶片预应力分析：按照待分析的透平机械叶片的三维模型和材料参数，建立透平机械叶片的有限元模型，对透平机械叶片进行预应力分析，提取透平机械叶片结构的整体预应力场，同时，若透平机械叶片上带有连接件，则计算连接件接触部分的接触刚度，并存储获得的预应力场数据以及接触刚度的计算结果；

2)叶片数据预处理：根据步骤1)建立的有限元模型在CPU和GPU上并行计算，得到透平机械叶片中各单元的初始刚度矩阵{K₀}和质量矩阵{M}，利用透平机械叶片的转速数据以及步骤1)获得的预应力场数据，对透平机械叶片有限元模型的各个单元进行修正，获得考虑旋转软化和预应力效应的各单元的实际刚度修正矩阵{K}，并由CPU将各单元的实际刚度修正矩阵{K}组装为总体刚度矩阵[K]，将各单元的质量矩阵{M}组装为总体质量矩阵[M]；

3)叶片整体结构约束条件处理：首先按照透平机械叶片的结构特点对该叶片的轮缘或叶根部分进行边界位移约束，然后进行叶根和轮缘的接触边界耦合，同时，对于带有连接件的透平机械叶片，按照步骤1)获得的连接件接触部分的接触刚度进行连接件间的接触耦合，对步骤2)获得的总体刚度矩阵[K]和总体质量矩阵[M]进行修正；

4)广义特征值和特征向量的提取：采用CPU+GPU异构并行计算，对步骤3)修正后的透平机械叶片的总体刚度矩阵和总体质量矩阵进行广义特征值λ和特征向量{φ}的提取；

5)叶片固有频率/固有振型输出：将步骤4)获得的广义特征值λ转化为透平机械叶片的固有频率f，将获得的特征向量{φ}转化为透平机械叶片的固有振型，绘制透平机械叶片的固有振型云图；

6)绘制叶片振动安全图或坎贝尔图：根据步骤5)中固有振型云图的结果，判断透平机械叶片固有振型的振动类型，对于自由透平机械叶片绘制其频率分布曲线，对于带有连接件的透平机械叶片绘制其振动安全图，对于多转速下的振型结果，按照不同转速的固有频率绘制透平机械叶片振动的坎贝尔图，完成对透平机械叶片固有频率的特性分析。

所述步骤1)中计算连接件接触部分的接触刚度的具体步骤为：

A)提取连接件接触面上所有接触节点对的法向接触力和法向相对位移；

B)对于节点编号为i、j的接触节点对，连接件接触面上该接触节点对的法向接触刚度K_n通过公式得到，其中p_n为连接件接触面上该接触节点对的法向接触力，u_n为连接件接触面上该接触节点对的法向相对位移；

C)根据Hertz理论，计算接触面局部坐标系下的等效刚度矩阵{K′_ij}，

其中v为连接件材料的泊松比；

D)通过公式{K_ij}＝R^T{K′_ij}R，获得变换为整体主坐标下的接触刚度矩阵{K_ij}，其中R为主坐标位移与接触面局部坐标位移的变换矩阵；

E)计算连接件接触面上所有接触节点对的接触刚度矩阵{K_ij}，即得到连接件接触部分的接触刚度。

所述步骤2)中各单元的实际刚度修正矩阵{K}＝{K₀}+{K_σ}-{K_c}，其中{K_σ}为根据预应力场数据获得的刚度矩阵，{K_c}为透平机械叶片的转速矩阵。

所述步骤3)中的边界位移约束是对总体刚度矩阵[K]和总体质量矩阵[M]进行修正，得到修正后的总体刚度矩阵和修正后的总体质量矩阵其中 [A]为由原位移向量{δ}变换为新位移向量的变换系数矩阵。

所述步骤3)中按照连接件接触部分的接触刚度进行连接件间的接触耦合，是利用整体主坐标下的接触刚度矩阵{K_ij}对修正后的总体刚度矩阵进行再次修正，得到二次修正后的总体刚度矩阵

其中{K_i}和{K_j}为主对角线上的第i个节点和第j个节点位置的矩阵。

所述步骤4)中，当透平机械叶片带有连接件时，广义特征值λ和特征向量{φ}通过公式计算得到，当透平机械叶片不带有连接件时，广义特征值λ和特征向量{φ}通过公式计算得到。

所述步骤5)中通过公式将广义特征值λ转化为透平机械叶片的固有频率f。

所述步骤5)中将特征向量{φ}转化为透平机械叶片的固有振型，得到向量[X]、[Y]、[Z]、[S]，其中[X]为透平机械叶片结构的X方向固有振型，[Y]为透平机械叶片结构的Y方向固有振型，[Z]为透平机械叶片结构的Z方向固有振型，[S]为透平机械叶片结构的总固有振型。

与现有技术相比，本发明具有以下有益效果：

与传统的计算方法相比，采用本发明的分析方法，操作人员无需关心步骤2)中的预应力场数据和连接件接触效应如何作用于整体叶片，本发明提供了准确且快速的方法，操作者只需直接导入步骤1)存储的预应力场数据和接触刚度的计算结果即可；步骤2)在CPU+GPU上同时进行单元矩阵并行计算和组装，可以有效提高计算速度；步骤3)中的约束条件处理可以避免操作人员进行复杂的节点选择，只需在初始的叶片有限元模型中选择好节点集合，计算过程中可以自行处理；步骤4)采用异构并行的方式提取广义特征值，合理地安排CPU和GPU的计算任务，极大地减少了计算时间；在步骤5)和步骤6)中，本发明提供了全面的计算数据输出方式，包括了透平叶片的固有振型云图、频率分布曲线、振动安全图、坎贝尔图，设计人员可以很方便地查看叶片的固有频率特性分析结果。

综上所述，本发明旨在向透平叶片设计方面的技术人员提供以下的便利：简化操作，设计人员不用进行专门的软件使用培训，就可以进行对透平叶片的固有频率特性计算和分析；利用CPU+GPU强大的并行计算能力，提高传统工程领域有限元方法的计算速度，为透平叶片的设计计算提供准确、快速的叶片频率特性分析结果，极大地缩短设计计算周期，使工程设计人员可以将更多精力集中在如何获得叶片设计的优化方案上，具有重要的工程应用价值。

附图说明

图1为透平叶片频率特性分析流程图；

图2为CPU-GPU计算安排以及数据交换过程示意图；

图3为实例一中的自由叶片的模型示意图；其中(a)为整圈叶片-轮缘几何模型，(b)为叶根-轮缘接触面示意图；(c)为叶片局部有限元网格模型。

图4为实例一中的自由叶片振型计算结果；其中(a)为1阶切向振动，(b)为1阶轴向振动，(c)为1阶扭转振动，(d)为2阶切向振动，(e)为2阶弯曲1阶扭转振动，(f)为3阶切向振动；

图5为实例一中的自由叶片动频值分布图；

图6为实例二中的带有连接件叶片的模型示意图；其中(a)为整圈叶片-轮缘几何模型，(b)为叶片连接件(阻尼块和围带)示意图；(c)为叶片局部有限元网格模型。

图7为实例二中的带有连接件叶片整圈振型计算结果；其中(a)为1阶节圆(1阶)，(b)为2阶节圆(16阶)，(c)为3阶节圆(29阶)，(d)为1节径1阶(2阶)，(e)为2节径1阶(4阶)，(f)为3节径1阶(6阶)，(g)为4节径1阶(8阶)，(h)为5节径1阶(10阶)，(i)为6节径1阶(12阶)，(j)为7节径1阶(14阶)，(k)为8节径1阶(17阶)，(l)为1节径2阶(21阶)，(m)为2节径2阶(38阶)，(n)为3节径2阶(46阶)，(o)为4节径2阶(54阶)，(p)为5节径2阶(67阶)，(q)为6节径2阶(69阶)，(r)为7节径2阶(78阶)，(s)为8节径2阶(90阶)，(t)为1节径3阶(30阶)，(u)为2节径3阶(34阶)，(v)为3节径3阶(42阶)，(w)为4节径3阶(50阶)，(x)为5节径3阶(61阶)，(y)为6节径3阶(71阶)，(z)为7节径3阶(84阶)，(a’)为8节径3阶(96阶)。

图8为实例二中的带有连接件叶片的振动安全图。

具体实施方式

下面结合附图对本发明做进一步详细说明。

本发明提供的基于CPU+GPU异构并行计算的透平机械叶片固有频率特性分析方法，公开了透平机械叶片的固有频率特性分析完整的计算流程以及具体实现方式，该方法首先需要按照待分析的透平机械叶片的三维模型和材料参数建立其有限元模型，对叶片进行预应力分析；其次进行叶片网格数据预处理，在CPU和GPU上同时计算单元的刚度矩阵和质量矩阵，并组装为叶片结构的总刚度和质量矩阵；再次对叶片和轮缘的约束条件进行设置，包括边界的刚性和弹性位移约束，叶根-轮缘或连接件的接触耦合，修正总刚度矩阵；然后采用CPU+GPU异构并行算法提取总刚度和质量矩阵的广义特征值和特征向量，提升计算效率；之后将特征值和特征向量转化为叶片的频率以及振型并输出；最后判断固有振型的振动类型，据此绘制叶片的频率曲线分布，振动安全图或坎贝尔图。最后通过两个实际叶片的算例验证了该方法的实用性。

本发明提供的基于CPU+GPU异构并行计算的透平机械叶片固有频率特性分析方法，适用于透平机械叶片固有频率的分析计算。本发明专门为透平机械叶片的固有频率特性分析所设计，提高工程设计人员操作的便捷性；通过与CPU+GPU异构并行计算平台结合，大幅度提升透平机械叶片频率特性分析的计算速度。

本发明提供的基于CPU+GPU异构并行计算的透平机械叶片固有频率特性分析方法，其步骤如下：

1.叶片预应力分析：按照待分析的透平机械叶片的三维模型和材料参数建立其有限元模型，对叶片进行预应力分析，提取叶片结构的整体预应力场，对于带有连接件的叶片，还需要考虑连接件的接触效应，并计算连接件接触部分的接触刚度，存储预应力计算结果以及接触刚度的结果；

2.叶片网格数据预处理：根据有限元模型在CPU和GPU上并行计算得到透平机械叶片的单元刚度矩阵和单元质量矩阵，利用叶片的转速数据以及第1步的预应力场，对叶片有限元模型的各个单元进行修正，获得考虑旋转软化和预应力效应的单元刚度矩阵修正，并由CPU将之组装为总体刚度矩阵和质量矩阵；

3.叶片整体结构约束条件处理：修正第二步获得的总体刚度矩阵:首先按照叶片的结构特点对叶片的轮缘或叶根部分进行边界位移约束，然后处理叶根和轮缘的接触边界耦合问题，对于带有连接件的叶片还需按照第1步获得的连接件的接触刚度处理连接件间的接触耦合问题；

4.广义特征值和特征向量提取：采用CPU+GPU异构并行计算对叶片结构的总体刚度矩阵和质量矩阵进行广义特征值和特征向量的提取(该步骤的计算时间一般可以占到总计算时间的80％以上，利用CPU+GPU异构并行计算可以大幅度提升计算效率，有效地减少计算时间)；

5.叶片固有频率/振型输出：将第4步获得的广义特征值和特征向量转化为叶片的固有频率和固有振型，(获得的固有频率和振型输出为Tecplot等软件可读取的数据格式)，绘制叶片的固有振型云图；

6.绘制叶片振动安全图或坎贝尔图：根据第5步中振型云图的结果，判断固有振型的振动类型，对于自由叶片绘制频率分布曲线，对于带有连接件的叶片绘制叶片的振动安全图，对于多转速下的振型结果，按照不同转速的固有频率可以绘制叶片振动的坎贝尔图。

按照上述技术方案，本发明提供了一种完善的透平机械叶片固有频率特性分析的计算流程和操作实现方式，其过程如图1所示，并在计算过程中采用了CPU+GPU异构并行优化。下面对其中各步骤的实现方式进行具体说明。

1.给定透平机械叶片的三维模型和材料参数，在网格剖分软件ANSA中建立叶片有限元模型，叶片的预应力分析可以在专用的有限元软件(如ANSYS/ABAQUS)中完成，将结果中的预应力场数据以二进制文件格式导出并存储，提取连接件接触面上所有接触节点对的法向接触力和法向相对位移后，做如下处理：

A.对于其中某节点对(节点编号为i，j)，以接触面法向接触力p_n与接触面的法向相对位移u_n的比值作为接触面的法向接触刚度K_n

B.根据Hertz摩擦理论，计算接触面局部坐标系下等效刚度矩阵{K′_ij}为

其中，ν为接触面的接触摩擦系数。

C.获得变换为整体主坐标下的刚度矩阵{K_ij}，式中R为主坐标位移与接触面局部坐标位移的变换矩阵。

{K_ij}＝R^T{K′_ij}R

D.将获得接触面上所有节点对的接触刚度矩阵{K_ij}存储在文件中。

2.叶片网格信息预处理：首先按照1中获得的叶片有限元模型计算有限元模型中各个单元的初始刚度矩阵{K₀}和质量矩阵{M}，根据预应力场数据获得刚度矩阵{K_σ}，考虑叶片的转速矩阵{K_c}，各单元的实际刚度修正矩阵{K}表示为

{K}＝{K₀}+{K_σ}-{K_c}

获取单元矩阵可以同时在CPU和GPU上组织算法，CPU上调用线性代数库blas加快速度，GPU可以每次同时计算多个单元刚度矩阵在多个位置的元素；考虑到CPU和GPU负载均衡问题；在16核32线程节点机上CPU可同时进行30个单元矩阵计算，经过测试，此时单块GPU(Nvidia Tesla k20)计算12个单元可以达到延迟最小。

最后将{K}和{M}按照有限元模型中单元和节点的对应信息分别组装为总体刚度矩阵[K]和总体质量矩阵[M]，由于该过程极易导致线程冲突，且一般计算时间很短，故只需用CPU进行串行计算即可。

3.叶片结构约束条件处理：首先是对叶根承载面或轮缘的进/出气边节点施加位移弹性约束或刚性约束，一般位移约束问题为：

其中，δ_i为第i个节点的被约束位移自由度，C_i为第i个节点的约束系数，C₀为被约束自由度之间的耦合常数；该约束问题可以通过一个线性变换统一转化为下式的约束问题，假设前n个自由度被约束；选取一个主自由度δ₁，约束方程化为

构造新的节点位移向量为

新位移向量与原位移向量{δ}变换关系为

简记为对第2步得到的[K]和[M]做变换，获得修正后的刚度矩阵和质量矩阵：

这样处理的好处是叶片频率分析问题中常见边界约束统一为一个形式，方便了程序的编写，而且[A]矩阵接近于一个单元矩阵，采用CSR格式存储后存储量极小，在GPU上借助相应的线性代数库cublas可以快速运算。

在叶根-轮缘接触面耦合问题中，可以事先将叶根以及轮缘部分有限元模型网格剖分成完全相同的形式，接触面耦合问题简化为两组几乎重合的节点组的所有自由度绑定问题，然后将两部分节点在一定容差范围内合并，这样在第1步的有限元模型中即可完成该耦合问题。采用这样的处理方式既可以获得很高的计算精度，并节省编程的工作量。

对于连接件的接触面耦合问题，利用第1步获得的第i个节点和第j个节点的接触刚度矩阵{K_ij}二次修正获得的总刚度矩阵其中{K_i}和{K_j}为主对角线上的第i个节点和第j个节点位置的矩阵。二次修正后的刚度矩阵为

由于一般连接件上的接触节点并不多，因此这步的计算时间很短，可以直接在CPU上组织计算。

4.根据第3步获得的刚度矩阵和质量矩阵得到广义特征值问题

叶片的固有频率问题需要提取前s(大于100)个特征值λ和特征向量{φ}，可以采用重启动的Block Lanczos算法，每次Lanczos算法循环内包括的运算按照耗时程度排列依次为：

1)一次线性方程组(解为向量组)求解

2)两次Lanczos向量完全正交化

3)两次稀疏矩阵-向量组乘积

4)三次密集矩阵-向量组乘积

5)两次三元向量组运算

6)一次小规模矩阵的Cholesky分解

分析其主要运算量，第1步和第2步一般占据了总运算的80％以上时间。其中，第2-5步均为大规模矩阵-向量组运算，适合在GPU上计算，最终的求解特征向量为一次矩阵-向量组运算，也应在GPU上计算；其次，我们可以通过保留中间矩阵，将稀疏矩阵-向量组相乘运算3)改为一个小规模密集矩阵-向量组乘积运算4)提高计算效率；此外，而6)中由于计算矩阵规模非常小，在CPU上获得更高的效率。如上述并行方法优化传统的Lanczos方法可以获得明显的速度提升。

特别地，对第1)步，如果叶片是刚度较大的短叶片，求解方程组主要由GPU承担(迭代法)求解效率更高；对于刚度较小(柔度较大)的中长型叶片，利用CPU求解主要部分(消元法)求解更快；因此第1)步的安排需要按具体问题的特点分配计算任务。

如果问题规模巨大，单个设备(CPU/GPU)计算时申请空间超出了内存/显存，计算效率会明显降低，此时可以优先将计算参与度小的大型密集型向量放入计算任务小的设备上，而不是直接利用硬盘交换数据；当然如果问题规模超过了所有内存与显存的总量，此时必须使用硬盘做数据交换。

由于硬盘与GPU无法直接进行数据交换，所有的刚度/质量矩阵以及中间变量的生成均需现在内存中存储，在需要GPU计算时再传递，传递后释放掉内存空间，总结广义特征值提取的并行算法计算任务分配如图2所示。

5.对第4步获得的s个特征值做如下处理转化为叶片结构的固有频率，并输出到文件中：

对第4步的s个特征向量做如下操作转化为叶片的振型向量并输出：

1)分别提取特征向量{φ}中每个节点对应3个自由度处的值，组成3个向量[X],[Y],[Z]，其中N为叶片有限元模型总节点数

2)计算总的振型向量[S]

3)将[X],[Y],[Z],[S]与有限元模型节点编号对应并输出，其中[X]即为叶片结构的X方向固有振型，[Y]即为叶片结构的Y方向固有振型，[Z]即为叶片结构的Z方向固有振型，[S]即为叶片结构的总固有振型。

该步骤过程的转化适用于在GPU上并行计算，对于节点数较小的模型加速不明显，但对于节点数很多的模型计算速度可以提升几十倍；处理完毕后将固有频率和振型输出到文件中。

6.绘制固有频率分布曲线或振动安全图以及坎贝尔图：对于自由叶片，不仅需要查看整圈振型图，还需查看其中单个叶片的振型云图，最后按照总振型云图判断单个叶片的阶次和弯曲扭转形式，绘制固有频率分布曲线；对于存在连接件的叶片，主要查看整圈的振型云图，按照Z方向振型判断整圈叶片振动的阶次和节径数，将固有频率整理为振动安全图，如果计算了叶片在多个转速下的固有振型，可以将不同转速下的固有频率结果绘制为坎贝尔图，为叶片振动安全性评估提供依据。

下面结合两个具体的实例对本发明提供的基于CPU+GPU异构并行计算的透平机械叶片固有频率特性分析方法进行具体说明。

实例一

实例一以某整圈自由叶片(不存在连接件)为例，三维模型为具有燕尾型叶根和轮缘，共43个叶片，图3中(a)给出了该叶片的整圈几何模型，(b)为燕尾型的叶根-轮缘接触面示意图，叶片及轮缘的材料参数见表1。网格剖分在ANSA软件中完成，图3中(c)局部的有限元网格模型，网格主体采用8节点六面体单元，在叶根过渡部分采用四面体以及退化单元进行网格划分，整体共具有760885个单元，728400个节点。

表1某自由叶片和轮缘的材料属性

对叶片施加工作转速为3000r/min的离心载荷，进行预应力分析；单元矩阵在CPU(Intel Xeon E5-2650)和GPU(Nvidia Tesla K20c)上计算，并在CPU上组装总矩阵；施加边界条件为轮缘进气端面的切向和轴向刚性位移约束，出气端面的切向位移刚性约束和轴向位移弹性约束，轮缘接触面和叶根接触面网格节点进行合并；提取前250阶特征值和特征向量，计算固有频率值并判断其振动类型如表2，绘制振型云图如图4，获得其中单个自由叶片的振动类型，最后绘制动频值分布图如图5。

该问题在32G内存，Intel Xeon E5-2650/Nvidia Tesla K20c单节点计算机上运行，ANSYS计算时间2763s，约为45分钟，VIB程序计算时间为954s，约16分钟；相对于ANSYS，VIB程序计算的加速比可以达到2.90倍，该方法加速效果良好。

表2整圈自由叶片固有动频(前250阶)

实例二

实例二以某带有连接件(围带和阻尼块)的整圈叶片为例，三维模型采用菌型叶根，通过阻尼块以及围带发生接触作用将整圈叶片连接起来，共92个叶片，图6中(a)给出了该叶片的整圈几何模型，(b)为连接件(围带和阻尼块)的示意图，叶片及轮缘的材料参数见表3。网格剖分在ANSA软件中完成，图6中的(c)为叶片局部的有限元网格模型，主体采用8节点六面体单元，在叶根和围带的过渡部分采用四面体以及退化单元进行网格划分，整体模型的节点总数为1176877，单元总数为1337312。

表3某带有连接件叶片和轮缘的材料属性

对叶片施加工作转速为3000r/min的离心载荷，设置围带和阻尼块的接触摩擦系数0.2，进行预应力分析；单元矩阵在CPU(Intel Xeon E5-2650)和GPU(Nvidia TeslaK20c)上计算，并在CPU上组装总矩阵；施加边界条件为轮缘进气端面的切向和轴向刚性位移约束，出气端面的切向位移刚性约束和轴向位移弹性约束，轮缘接触面和叶根接触面网格节点进行合并，阻尼块和围带间的接触刚度按照接触面间的接触力和相对位移进行修正；提取前100阶特征值和特征向量，计算固有频率并判断振动阶次和节径数，汇总如表4，计算固有振型并输出振型云图如图7；最后绘制工作转速下3000r/min的振动安全图如图8。与实例一采用相同配置计算机，ANSYS计算时间为3234s，VIB程序计算时间为873.57s，并行程序计算加速比可以达到3.7左右，该方法加速效果明显。

表4某整圈带连接件叶片固有频率值

Claims

1.一种基于CPU+GPU异构并行计算的透平机械叶片固有频率特性分析方法，其特征在于，包括以下步骤：

2.根据权利要求1所述的基于CPU+GPU异构并行计算的透平机械叶片固有频率特性分析方法，其特征在于，所述步骤1)中计算连接件接触部分的接触刚度的具体步骤为：

其中v为连接件材料的泊松比；

3.根据权利要求1所述的基于CPU+GPU异构并行计算的透平机械叶片固有频率特性分析方法，其特征在于，所述步骤2)中各单元的实际刚度修正矩阵{K}＝{K₀}+{K_σ}-{K_c}，其中{K_σ}为根据预应力场数据获得的刚度矩阵，{K_c}为透平机械叶片的转速矩阵。

4.根据权利要求1所述的基于CPU+GPU异构并行计算的透平机械叶片固有频率特性分析方法，其特征在于，所述步骤3)中的边界位移约束是对总体刚度矩阵[K]和总体质量矩阵[M]进行修正，得到修正后的总体刚度矩阵和修正后的总体质量矩阵其中[A]为由原位移向量{δ}变换为新位移向量的变换系数矩阵。

5.根据权利要求4所述的基于CPU+GPU异构并行计算的透平机械叶片固有频率特性分析方法，其特征在于，所述步骤3)中按照连接件接触部分的接触刚度进行连接件间的接触耦合，是利用整体主坐标下的接触刚度矩阵{K_ij}对修正后的总体刚度矩阵进行再次修正，得到二次修正后的总体刚度矩阵

6.根据权利要求5所述的基于CPU+GPU异构并行计算的透平机械叶片固有频率特性分析方法，其特征在于，所述步骤4)中，当透平机械叶片带有连接件时，广义特征值λ和特征向量{φ}通过公式计算得到，当透平机械叶片不带有连接件时，广义特征值λ和特征向量{φ}通过公式计算得到。

7.根据权利要求1所述的基于CPU+GPU异构并行计算的透平机械叶片固有频率特性分析方法，其特征在于，所述步骤5)中通过公式将广义特征值λ转化为透平机械叶片的固有频率f。

8.根据权利要求1所述的基于CPU+GPU异构并行计算的透平机械叶片固有频率特性分析方法，其特征在于，所述步骤5)中将特征向量{φ}转化为透平机械叶片的固有振型，得到向量[X]、[Y]、[Z]、[S]，其中[X]为透平机械叶片结构的X方向固有振型，[Y]为透平机械叶片结构的Y方向固有振型，[Z]为透平机械叶片结构的Z方向固有振型，[S]为透平机械叶片结构的总固有振型。