CN106570204A - 一种基于cpu+gpu异构并行计算的透平机械叶片静强度特性分析方法 - Google Patents

Abstract

本发明公开了一种基于CPU+GPU异构并行计算的透平机械叶片静强度特性分析方法，先建立有限元模型，计算模型的总刚度矩阵，再计算透平机械叶片的离心载荷向量和气动载荷向量，然后对节点进行位移约束和耦合，修正总刚度矩阵；再用CPU+GPU并行求解总体刚度矩阵和载荷向量形成的方程组，获得节点位移向量，然后计算主应变和VonMises等效应力，绘制分布云图，最后进行安全考核。本发明针对透平机械叶片静强度分析设计，方便工程设计人员操作，同时采用的CPU+GPU并行算法可以有效地提高有限元方法的计算速度，为透平机械叶片设计提供准确、快速的叶片静强度特性分析结果，大幅度缩短透平机械叶片的设计周期。

Description

一种基于CPU+GPU异构并行计算的透平机械叶片静强度特性 分析方法

技术领域

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

背景技术

叶片在透平机械中承担着把热能转化为机械能的重要任务，是透平机械中最重要的零部件之一。叶片在运行中不仅承受着离心力、稳定气流力及交变应力的共同作用。随着透平功率的增大，末级叶片的离心载荷以及气流力载荷逐渐增大，叶片的工作环境趋于复杂恶劣，叶片损伤也是透平机械发生故障的主要原因，为了防止叶片在运行过程中损伤甚至破坏，在叶片设计前必需进行静强度校核。

叶片的静强度校核的主要任务是分析叶片在离心力、稳定气流力下的变形和应力分布。但由于透平机械叶片的模型复杂，边界条件以及载荷施加繁琐，而且基于CPU的有限元计算在处理叶片这种网格数量较大的模型(网格数量通常超过50万)往往计算速度缓慢，明显地加大了叶片的设计周期；此外，传统的有限元软件没有一个行之有效的针对叶片的高精度且操作简便的算法和流程，对透平设计人员提出了较高的技术要求，严重影响了透平设计计算的效率。

发明内容

本发明的目的在于提供一种基于CPU+GPU异构并行计算的透平机械叶片静强度特性分析方法，该方法针对透平机械叶片静强度分析设计，方便工程设计人员操作，同时采用的CPU+GPU并行算法可以有效地提高有限元方法的计算速度，能够大幅度缩短透平机械叶片设计周期。

为达到上述目的，本发明采用的技术方案为：

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

1)有限元模型建立：对于给定的透平机械叶片的三维模型和材料参数，在有限元网格剖分软件中建立透平机械叶片的有限元模型，将有限元模型的节点坐标与单元-节点对应关系存储在文件中，将叶根或轮缘要设置约束条件的节点以及要耦合的节点分别划分为集合，将要考核变形和应力部位的节点或单元输出到文件中，将透平机械叶片稳定工作时的气动数据记录在文件中；将透平机械叶片的预应力场数据记录在文件中，将透平机械叶片的材料参数以及工况参数输出到文件中；

2)计算总刚度矩阵[K]：读取步骤1)中透平机械叶片的有限元模型的节点和单元数据，先对单元进行染色，然后计算总刚度矩阵[K]的行索引和列索引向量以及在压缩存储中的索引值，再根据染色结果在CPU和GPU上同时计算所有单元刚度矩阵并组装为总刚度矩阵[K]；

3)计算总载荷向量{F}：根据步骤1)中存储的透平机械叶片的材料参数和工况参数，分别计算透平机械叶片因机械旋转引起的离心载荷向量以及因气流冲击导致的气动载荷向量，并合成为总载荷向量{F}；

4)边界条件处理：根据步骤1)中存储的叶根或轮缘要设置约束条件的节点以及要耦合的节点的集合，对该集合中的节点进行位移约束和耦合，并修正步骤2)中获得的总刚度矩阵和总载荷向量；

5)方程组求解，用CPU+GPU并行求解步骤4)修正后的总刚度矩阵和总载荷向量形成的方程组，获得节点位移向量；

6)应力结果处理：先将步骤5)获得的节点位移向量转化为单元高斯积分点上的应变和应力，并外插获得单元节点上的应变和应力，计算主应变和VonMises等效应力，再将单元节点的位移、应变和应力输出，绘制分布云图；

7)安全考核：统计步骤1)中考核变形和应力部位的节点的位移、主应变和Von-Mises等效应力，将其中的最大位移，最大主应变，最大等效应力值与叶片材料安全承载极限做对比，并输出考核结果，即完成了透平机械叶片静强度特性分析。

所述的步骤2)具体包括以下步骤：

2-1)根据步骤1)中文件存储的单元-节点对应关系，建立单元之间的连通矩阵[L]，其中L_ij表示第i个单元的第j个相邻单元编号；所述的相邻单元为有公共节点的单元；

2-2)依次遍历所有单元，根据连通矩阵[L]对单元进行染色分类，当两个单元相邻时，染上不同的颜色，获得颜色的个数以及每种颜色的单元编号；

2-3)根据步骤1)中文件存储的单元-节点对应关系，建立节点之间的关系矩阵[R]，其中R_ij表示与第i个节点相关的第j个节点的编号，并将每个节点的相关节点的编号按照从小到大排序，统计关系矩阵[R]中每个节点对应的相关节点的个数，组成相关和向量{M}，其中M_i表示第i个节点对应的相关节点总数；所述的相关节点为位于同一个单元的节点；

2-4)计算总刚度矩阵[K]的行索引向量和列索引向量以及在压缩存储中的索引值，对于第i个节点相关的第j个节点的刚度矩阵元素，其在总刚度矩阵的行索引IndexCol和列索引IndexRow分别为：

在压缩存储中的索引值IndexVector为：

其中U＝M₁+……+M_i-1；

2-5)选取一个颜色的单元集合，将该集合中单元分别分配给CPU和GPU，在CPU和GPU上分别计算各自分配到的单元的刚度矩阵；CPU上使用OpenMP进行多核并行计算，GPU上编写CUDA Kernel函数并行计算，计算方法如下：

A.进行网格质量检查，从步骤1)中存储的文件中分别读取CPU和GPU分配到的单元内的节点坐标，计算所述单元在所有高斯节点上的Jaccobi值，如果小于等于0，则退出，返回步骤1)重新进行有限元模型建立，否则进行下一步；

B.遍历CPU和GPU分配到的单元中的所有节点，计算几何子矩阵[B_L]_i，设N_i为第i个节点的插值函数，用表示对第k个局部坐标的偏微分，则第i个节点对应的几何子矩阵[B_L]_i为：

其中

C.计算第i个节点相关的第j个节点的线性刚度子矩阵[K_L]_ij和预应力刚度子矩阵[K_σ]_ij，

[K_L]_ij＝[B_L1i ^TEB_L1j+GB_L2i ^TB_L2j]，

其中，[I]为3阶单位矩阵，σ＝[σ_x,σ_y,σ_z,τ_xy,τ_xz,τ_yz]为预应力张量，E为叶片材料参数中的弹性模量，G为叶片材料参数中的剪切模量；

D.计算刚度子矩阵[K]_ij＝[K_L]_ij+[K_σ]_ij，并将所有的高斯节点位置的刚度子矩阵[K]_ij相加，获得最终的刚度子矩阵；

E.在GPU和CPU上按照索引向量IndexCol，IndexRow，IndexVector步骤将D中计算的刚度子矩阵加到总刚度矩阵的对应位置，利用MPI_WAIT()进行线程协调，直到GPU和CPU都完成分配的计算任务；

2-6)若遍历完所有颜色集合，则进行步骤2-7)，否则继续进行步骤2-5)；

2-7)将CPU计算获得的刚度矩阵和GPU获得的刚度矩阵相加为总刚度矩阵[K]。

所述的步骤3)具体包括以下步骤：

3-1)选择一种颜色的单元集合，将该集合中单元分别分配给CPU和GPU，由GPU和CPU分别计算各自分配到的单元的离心载荷向量其中ρ为叶片材料参数中的材料密度，ω为透平机械叶片工况参数中的转速，{r}为单元内高斯节点到透平转轴的垂直矢量，[N]为单元型函数；

3-2)计算CPU和GPU上分配到的单元的所有节点在总载荷向量的索引IndexForce，对于第i个节点：

将所述单元的离心载荷向量按照IndexForce的索引加入总载荷向量；

3-3)若遍历完所有颜色集合，则进行步骤3-4)，否则继续进行步骤3-2)；

3-4)选取步骤1)存储的透平机械叶片稳定工作时的气动数据中的一个单元表面S^e，用CPU计算单元气动载荷向量其中p为单元表面受到的气流冲击压强，{n}为单元表面的单位法矢量，[N]为单元型函数；

3-5)将所述单元的气动载荷向量按照IndexForce的索引加入总载荷向量；

3-6)遍历完S集合，获得总载荷向量{F}；其中S集合为透平机械叶片稳定工作时的气动数据中的所有单元表面。

所述的步骤4)具体包括以下步骤：

4-1)提取步骤1)存储的叶根或轮缘要设置约束条件的节点的集合的数据，对节点被约束自由度对应刚度矩阵位置进行修正，对被约束自由度的对角线元素赋值1，并将其余行和列元素全部置0，并将总载荷向量{F}的对应位置取为约束值，获得修正后的总刚度矩阵[K]^m和修正后的总载荷向量

4-2)提取步骤1)存储的叶根或轮缘要耦合的节点的集合的数据，对于第i个节点和第j个节点的耦合节点对，其法向和两个切向耦合刚度分别为K_n，K_t1，K_t2；根据节点的坐标计算主坐标位移与耦合局部坐标位移的变换矩阵T，然后计算主坐标系下的耦合刚度矩阵其中

4-3)利用二次修正总刚度矩阵[K]^m，修正后的总刚度矩阵为：

其中[K]_i和[K]_j为矩阵[K]^m主对角线上的第i个节点和第j个节点对应位置的子矩阵。

所述的步骤5)具体为：修正后的总刚度矩阵节点位移向量{δ}和修正后的总载荷向量构成的方程组为：

对于展弦比小于10的叶片，具体算法如下：

5-1)构造预处理矩阵[P]^-1：当展弦比小于3时，[P]^-1取修正后的总刚度矩阵的主对角线block矩阵，block的尺寸取5-10之间；当展弦比大于3小于10时，取修正后的总刚度矩阵的不完全cholesky分解；

5-2)在CPU上计算残差向量组{δ₀}为与节点位移向量{δ}维数相同的随机向量，构造预处理向量{z₀}＝[P]^-1{r₀}，并对辅助向量组{q₀}赋值，{q₀}＝{r₀}；

5-3)进行迭代循环，迭代数k＝0,1,2……

5-3-1)用CPU计算向量组内积{z_k}^T{r_k}，同时GPU上计算获得

5-3-2)用CPU计算{δ_k+1}＝{δ_k}+α_k{q_k}，同时GPU计算

5-3-3)判断||{r_k+1}||＜10^-7，满足则进入步骤d)，否则继续进行步骤5-3-4)；

5-3-4)求解预处理方程{z_k+1}＝[P]^-1{r_k+1}，当展弦比小于3时，在GPU上进行求解，当展弦比大于3时，在CPU和GPU上共同求解；

5-3-5)用CPU计算并获得{q_k+1}＝{z_k+1}+β_k+1{q_k}，返回步骤5-3-1)；

5-4)将获得的向量{δ_k+1}作为原方程组的解向量{δ}；

对于展弦比大于10的叶片，采用CPU+GPU异构式的稀疏矩阵的高斯消元法求解。

所述的步骤6)具体包括以下步骤：

6-1)根据步骤5)中获得的节点位移向量{δ}，选择一种颜色的单元集合，同时由GPU和CPU计算该集合内单元的高斯节点上的应变分量{E}^G和应力分量{S}^G，{E}^G＝[B]^e{δ}^e，{S}^G＝[D][B]^e{δ}^e，

其中E为叶片材料参数中的弹性模量，G为剪切模量，μ为泊松比，[B]^e为单元的几何矩阵，{δ}^e为单元的位移向量；

6-2)计算该单元内所有高斯节点的插值函数方阵[NG]，对每一个应变分量{E}^G和应力分量{S}^G分别外插，获得单元节点上的应变分量{E}^e和应力分量{S}^e，{E}^e＝[NG]^-1{E}^G，{S}^e＝[NG]^-1{S}^G；

6-3)依照步骤6-1)和6-2)对所有单元内的所有节点进行应力分量和应变分量的计算，对重复计算的节点的应力分量和应变分量分别取算术平均值，获得所有节点的实际应力分量{S}和实际应变分量{E}；

6-4)根据每个单元节点或高斯节点的实际应力分量和实际应变分量计算三个主应力S₁，S₂，S₃和三个主应变E₁，E₂，E₃；对每个节点计算其VonMises等效应力，VonMises等效应力的计算方法为

6-5)根据节点位移向量{δ}、实际应力分量{S}、实际应变分量{E}和VonMises等效应力，绘制叶片的位移、应变、应力分布云图以及变形图。

相对于现有技术，本发明的有益效果为：

本发明提供了一种支持CPU+GPU异构并行计算的透平机械叶片静强度特性有限元分析方法，该方法首先按照待分析的透平机械叶片建立有限元模型，将叶根或轮缘需要设置约束条件以及耦合的节点划分为集合，将需要考核变形和应力部位的节点或单元输出到文件中；计算模型的总刚度矩阵，通过对单元的染色，提高在CPU和GPU上计算单元刚度矩阵和组装总刚度矩阵的效率；施加载荷，计算透平机械叶片因机械旋转引起的离心载荷向量以及气流冲击导致的气动载荷向量，并合成为总载荷向量；边界条件处理，根据叶根或轮缘需要设置约束条件以及耦合的节点集合，对节点进行位移约束和耦合，并修正总刚度矩阵；然后用CPU+GPU并行求解总体刚度矩阵和载荷向量形成的方程组，获得节点位移向量，对于结构柔度不同(展弦比的不同)的叶片采用不同的算法求解；应力结果处理，将节点位移向量转化为单元高斯积分点上的应变和应力，并外插获得单元节点上的应变/应力，计算主应变和VonMises等效应力，最后将单元节点的位移、应变、应力输出，绘制分布云图；安全考核，统计叶片重点考核位置的位移，主应变，Von-Mises等效应力，将获得的最大位移、最大主应变、最大等效应力值与材料安全承载极限做对比，并输出考核结果。本发明给出了透平机械叶片的静强度分析完整的计算流程以及具体实现方式，专门针对透平机械叶片静强度分析设计，方便工程设计人员操作，同时采用的CPU+GPU并行算法可以有效地提高有限元方法的计算速度，为透平机械叶片设计提供准确、快速的叶片静强度特性分析结果，大幅度缩短透平机械叶片的设计周期。

附图说明

图1是本发明的透平机械叶片静强度有限元分析方法的流程图；

图2是总刚度矩阵的计算流程图；

图3是实例一中的悬臂梁模型示意图；其中(a)是示意图，(b)是网格剖分结果图；

图4是实例一中的悬臂梁模型静强度分析结果分布云图；其中(a)为悬臂梁静强度分析中X向位移变形分布(左侧为TBFEA结果，右侧为ANSYS计算结果)，(b)为悬臂梁静强度分析中Y向位移变形分布(左侧为TBFEA结果，右侧为ANSYS计算结果)，(c)为悬臂梁静强度分析中VonMises等效应力分布(左侧为TBFEA结果，右侧为ANSYS计算结果)；

图5是实例一中的悬臂梁模型静强度分析结果对比图；其中(a)为悬臂梁的部分节点VonMises等效应力分布曲线，(b)为悬臂梁模型采用不同并行方法时总刚度矩阵计算时间对比，(c)为悬臂梁模型TBFEA和ANSYS总计算时间对比；

图6是实例二中的长叶片模型示意图；其中(a)为长叶片模型整体网格(左侧为压力面网格，右侧为吸力面网格)，(b)为长叶片模型叶根局部网格，(c)为长叶片模型拉筋局部网格，(d)为长叶片模型围带局部网格；

图7是实例二中的长叶片结算结果分布云图；其中(a)为长叶片模型VonMises等效应力分布(TBFEA计算结果)，(b)为长叶片模型VonMises等效应力分布(ANSYS计算结果)，(c)为长叶片模型叶根平台圆角VonMises等效应力分布(TBFEA计算结果)，(d)为长叶片模型叶根平台圆角VonMises等效应力分布(ANSYS计算结果)，(e)为长叶片模型拉筋圆角VonMises等效应力分布(TBFEA计算结果)，(f)为长叶片模型拉筋圆角VonMises等效应力分布(ANSYS计算结果)，(g)为长叶片模型围带圆角VonMises等效应力分布(TBFEA计算结果)，(h)为长叶片模型围带圆角VonMises等效应力分布(ANSYS计算结果)。

具体实施方式

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

本发明提供了一种完善的透平机械叶片静强度特性分析的计算流程和操作实现方式，其过程如图1所示，并在计算过程中采用了CPU+GPU异构并行优化。下面对其中各步骤的实现方式进行说明。

1)有限元模型建立：

对于给定透平机械叶片的三维模型和材料参数，在有限元网格剖分软件中建立叶片有限元模型，将有限元模型的节点编号以及节点坐标压缩，即使节点和单元的编号从1开始，每次增加1，直到达到节点和单元个数，节点编号以及节点坐标保存在二进制文件Node中，单元编号以及单元包含的节点编号保存在二进制文件Element中，将叶根或轮缘需要设置约束条件的部分节点划分为集合B₁,B₂,…,B_Nb；将叶根或轮缘需要耦合的节点划分为集合C₁,C₂,…,C_Nc；将需要考核变形和应力部位的节点划分为集合D₁,D₂,…,D_Nd；将透平机械叶片稳定工作时的气动数据记录在文件Flow中；将叶片的预应力场数据记录在文件Stress中；将叶片的材料数据及工况数据记录在文件Blade中。

其中，B_i集合为约束节点编号以及节点在X/Y/Z三个方向上的自由度约束值；C_i集合为耦合节点对的节点编号以及法向和两个切向的耦合刚度大小；D_i集合为考核位置的单元编号或节点编号；Flow包括了需要施加气流载荷的单元表面节点编号以及其受到气流冲击压强大小；Stress为叶片有限元模型节点编号与对应的初始应力张量分量；Blade文件中包括了叶片材料参数：密度，弹性模量，剪切模量，泊松比，最大允许位移，安全承载极限应力、安全承载极限应变，工作转速，展弦比。

2)计算总刚度矩阵[K]：

读取叶片有限元模型的节点和单元数据，先对单元进行染色，然后对各个单元进行网格质量检查，根据染色结果在CPU和GPU上同时计算所有单元刚度矩阵并组装为总刚度矩阵；通过单元染色的方法将具有共同节点的单元区分开，减少并行计算的线程冲突，可以提高有限元方法的总刚度矩阵计算速度，流程可参见图2，具体的计算过程如下：

2-1)根据Element中单元和节点的对应关系，建立单元之间的连通矩阵[L]，其中L_ij表示第i个单元的第j个相邻单元(有公共节点的单元即为相邻单元)编号；

2-2)依次遍历所有单元，根据连通矩阵[L]对单元进行染色分类，当两个单元相邻时，染上不同的颜色，获得颜色的个数Ncolor以及每种颜色的单元编号；

2-3)根据Element中单元节点的对应关系，建立节点之间的关系矩阵[R]，其中R_ij表示与第i个节点相关的第j个节点的编号(位于同一个单元的节点为相关节点)，并将每个节点的相关节点的编号按照从小到大排序，统计关系矩阵[R]中每个节点对应的相关节点个数，组成相关和向量{M}，M_i表示第i个节点对应的相关节点总数；

2-4)计算总刚度矩阵[K]的行索引向量和列索引向量以及在压缩存储的索引值，对于第i个节点相关的第j个节点的刚度矩阵元素，其在总刚度矩阵的行索引IndexCol和列索引IndexRow分别为

在压缩存储的索引值IndexVector为

其中U＝M₁+……+M_i-1

2-5)选取一个颜色的单元集合，将该集合中单元分别分配给CPU和GPU，CPU和GPU计算单元个数的比例需要根据CPU和GPU计算能力进行测试然后确定，原则是保证各自计算负载均衡，以使时间延迟最小(例如测试GPU：1块Tesla k20c，CPU：双处理器Intel XeonE2650，GPU分配的单元个数为CPU的2-3倍)；然后在CPU和GPU上分别计算各自分配到的单元的刚度矩阵；CPU上使用OpenMP进行多核并行计算，GPU上编写CUDA Kernel函数并行计算，计算方法如下：

A.进行网格质量检查，从Node文件中分别读取CPU和GPU分配到的单元内所有节点坐标，计算所述单元在所有高斯节点上的Jaccobi值，如果小于等于0，则退出，返回步骤1)重新进行有限元模型建立，否则进行下一步；

B.遍历CPU和GPU分配到的单元中所有节点，计算几何子矩阵[B_L]_i，设N_i为第i个节点的插值函数，用表示对第k个局部坐标的偏微分，则第i个节点对应的几何子矩阵[B_L]_i为

其中

C.计算第i个节点相关的第j个节点的线性刚度子矩阵[K_L]_ij和预应力子刚度矩阵[K_σ]_ij

其中，[I]为3阶单位矩阵，σ＝[σ_x,σ_y,σ_z,τ_xy,τ_xz,τ_yz]为预应力张量，E为叶片材料参数中的弹性模量，G为叶片材料参数中的剪切模量；

D.计算刚度子矩阵[K]_ij＝[K_L]_ij+[K_σ]_ij，并将所有的高斯节点位置的矩阵[K]_ij相加，获得最终的刚度子矩阵；

E.在GPU和CPU上分别将单元刚度矩阵组装为总刚度矩阵，按照第4)步获得的索引向量IndexCol，IndexRow，IndexVector，然后将步骤D中计算的子刚度矩阵加到总刚度矩阵的对应位置，利用MPI_WAIT()进行线程协调，直到GPU和CPU都完成分配的计算任务；

2-6)如果遍历完所有颜色集合，则进行步骤2-7)，否则继续进行步骤2-5)；

2-7)将CPU计算获得的刚度矩阵和GPU获得的刚度矩阵相加为总刚度矩阵[K]。

3)计算总载荷向量{F}：

利用在步骤1)中的叶片材料参数和工况参数：材料密度ρ，旋转角速度ω，叶片表面单元受到的气流冲击压强p，分别计算透平机械叶片因机械旋转引起的离心载荷向量以及气流冲击导致的气动载荷向量，并合成为总载荷向量{F}。设{r}为单元内高斯节点到透平转轴的垂直矢量，{n}为单元表面的单位法矢量，[N]为单元的型函数。计算流程如下：

3-1)选择一种颜色的单元集合，将该集合中单元分别分配给CPU和GPU，由GPU和CPU分别计算各自分配到的单元的离心载荷向量

3-2)计算CPU和GPU上分配到的单元的所有节点在总载荷向量的索引IndexForce，对于第i个节点：

将单元中的离心载荷向量按照IndexForce的索引加入总载荷向量；

3-3)若遍历完所有颜色集合，则进行步骤3-4)，否则继续进行步骤3-2)；

3-4)选取步骤1)存储的透平机械叶片稳定工作时的气动数据中的一个单元表面S^e，用CPU计算单元气动载荷向量

3-5)将单元的气动载荷向量按照IndexForce的索引加入总载荷向量；

3-6)遍历完S集合，获得总载荷向量{F}；其中S集合为透平机械叶片稳定工作时的气动数据中的所有单元表面。

4)边界条件处理：

4-1)提取步骤1)中获得的集合B₁,B₂,…,B_Nb的数据，约束节点的编号以及在X/Y/Z三个方向上的自由度约束值，对节点被约束自由度对应刚度矩阵位置进行修正，仅对被约束自由度的对角线元素赋1，并将其余行和列元素全部置0，并将总载荷向量{F}的对应位置取为约束值，以此修正刚度矩阵[K]和{F}，获得修正后的总刚度矩阵[K]^m和修正后的总载荷向量

4-2)提取步骤1)中获得的集合C₁,C₂,…,C_Nc的数据，对于第i个节点和第j个节点的耦合节点对，其法向和两个切向耦合刚度分别为K_n，K_t1，K_t2；根据节点的坐标计算主坐标位移与耦合局部坐标位移的变换矩阵T，然后计算主坐标系下的耦合刚度矩阵其中

4-3)利用二次修正总刚度矩阵[K]^m，则修正后的总刚度矩阵为

其中[K]_i和[K]_j为总刚度矩阵[K]^m主对角线上的第i个节点和第j个节点对应位置的子矩阵。

5)方程组求解：

用CPU+GPU并行求解修正后的总刚度矩阵和总载荷向量形成的方程组，获得节点位移向量{δ}；对于结构柔度不同(展弦比不同)的叶片采用不同的算法求解。

本方法中求解的线性方程组中总刚度矩阵是稀疏的，半正定的且对称，则总刚度矩阵节点位移向量{δ}和总载荷向量构成的方程组形式如下：

对于展弦比小于10的叶片，具体的算法如下：

5-1)构造预处理矩阵[P]^-1：当展弦比小于3时，[P]^-1取修正后的总刚度矩阵的主对角线block矩阵，block的尺寸取5-10之间；当展弦比大于3小于10时，取修正后的总刚度矩阵的不完全cholesky分解；

5-2)在CPU上计算残差向量组{δ₀}为与节点位移向量{δ}维数相同的随机向量，构造预处理向量{z₀}＝[P]^-1{r₀}，并对辅助向量组{q₀}赋值，{q₀}＝{r₀}；

5-3)进行迭代循环，迭代数k＝0,1,2……

5-3-1)用CPU计算向量组内积{z_k}^T{r_k}，同时GPU上计算获得

5-3-2)用CPU计算{δ_k+1}＝{δ_k}+α_k{q_k}，同时GPU计算

5-3-3)判断||{r_k+1}||＜10^-7，满足则进入步骤d)，否则继续进行步骤5-3-4)；

5-3-4)求解预处理方程{z_k+1}＝[P]^-1{r_k+1}，当展弦比小于3时，在GPU上进行求解，当展弦比大于3时，在CPU和GPU上共同求解；

5-3-5)用CPU计算并获得{q_k+1}＝{z_k+1}+β_k+1{q_k}，返回步骤5-3-1)；

5-4)将获得的向量{δ_k+1}作为原方程组的解向量{δ}。

对于展弦比大于10的叶片，采用CPU+GPU异构式的稀疏矩阵的高斯消元法求解。

6)应力结果处理：

6-1)根据步骤5)中获得的节点位移向量{δ}，遍历步骤2)中的单元染色集合，对同一种颜色集合的单元，同时由GPU和CPU计算该集合内单元的高斯节点上的应变分量{E}^G和应力分量{S}^G，GPU和CPU各自分配的单元数量比例按照第2)步中的分配比例，最后将该单元编号与对应的高斯节点的应变分量和应力分量保存在文件Str_Element中。高斯节点的应变分量和应力分量的计算方式如下：

{E}^G＝[B]{δ}^e，{S}^G＝[D][B]{δ}^e

其中E为叶片材料参数中的弹性模量，G为剪切模量，μ为泊松比，[B]^e为单元的几何矩阵，{δ}^e为单元的位移向量；

6-2)计算该单元内所有高斯节点的插值函数方阵[NG]，对每一个应变分量{E}^G和应力分量{S}^G分别外插，获得单元节点上的应变分量{E}^e和应力分量{S}^e，并将其保存在文件Str_Node中。单元节点的应变和应力的外插方式如下：

{E}^e＝[NG]^-1{E}^G，{S}^e＝[NG]^-1{S}^G

6-3)在遍历单元的过程中若下一个单元含有之前单元的节点，则再对该节点的应力分量和应变分量进行计算，直至遍历完所有节点，累加所有该节点在不同单元中获得的应力分量和应变分量分别取算术平均值，获得所有节点的实际应力分量{S}和实际应变分量{E}。

6-4)根据每个单元节点或高斯点的实际应力分量和实际应变分量可以计算其三个主应力S₁，S₂，S₃和三个主应变E₁，E₂，E₃；对每个节点计算其VonMises等效应力，单元内高斯节点的等效应力结果保存在Str_Element文件，单元节点的等效应力结果保存在Str_Node文件。VonMises等效应力的计算方法如下：

6-5)将Node，Element，Str_Node中的数据合并输出在文件Plot中，根据节点位移向量{δ}、实际应力分量{S}、实际应变分量{E}和VonMises等效应力，可以在Tecplot绘制叶片的位移，应变，应力分布云图以及变形图。

7)安全考核：

统计步骤1)中考核变形和应力部位的节点的位移、主应变和Von-Mises等效应力，将其中的最大位移，最大主应变，最大等效应力值与叶片材料安全承载极限做对比，并输出考核结果，即完成了透平机械叶片静强度特性分析。

首先遍历集合D₁,D₂,…,D_Nd，在CPU和GPU上同时进行比较位移值max(δ)和等效应力值max(S_vonMises)和最大主应变max(E₁)；考核单元应力时，从Str_Element中提取高斯节点的计算结果比较，考核节点应力时，从Str_Node中提取单元节点的计算结果；比较后获得D_i中最大位移，等效应力值，最大主应变以及最大值对应的单元编号/节点编号/节点坐标，并输出在结果文件Result中。

考核位置一般取为叶根平台圆角，拉筋圆角，围带圆角位置，叶顶位置，将各个位置的最大主应变和第1)步中的安全承载极限应变做对比，将最大等效应力与第1)步中的安全承载极限应力做对比，如果在第1)步设置了允许位移，则将最大位移与允许位移对比，超过允许值则认为该叶片存在不安全因素，将各个位置的安全考核结论输出在结果文件Result中。

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

本说明书的算例的计算性能在工作站(1块Tesla k20c/双处理器Intel XeonE2650/128G物理内存)进行测试。基于CUDA架构按照本发明提供的方法编制程序(TBFEA程序)计算结果以及计算速度与商用软件ANSYS进行对照，ANSYS软件采用双CPU处理器16核并行模式。

实例一

采用悬臂梁模型进行测试，悬臂梁尺寸为X向50mm，Y向4mm，Z向3mm，悬臂梁的几何模型如图3中(a)，有限元网格模型如图3中(b)，材料参数为弹性模量210GPa，泊松比0.3，密度7850kg/m³。边界条件为YOZ平面上节点全部约束，载荷为离心力载荷，绕Z轴角速度为100rad/s。为了验证并行计算在不同网格密度下的加速效果，计算的网格数量按照表1取值。

表1测试计算悬臂梁网格数量

图4中对比了该程序计算悬臂梁模型时TBFEA与ANSYS结果的分布云图，为X/Y向的位移分布云图和VonMises等效应力分布云图；图5的(a)中给出了TBFEA与ANSYS中相同节点的VonMises等效应力值，可以看出本发明的计算结果与现有商用软件基本相同，达到了工程计算要求的精度。

图5的(b)中对比了不同网格数量下采用2CPU+GPU与一般采用2CPU的总刚度矩阵/载荷向量计算修正时间，在500k网格规模以上加速比基本达到了3-4倍左右，在4000k网格数量时，1CPU计算时间为152s，2CPU计算时间为81s，2CPU+GPU的计算时间为25s。图5的(c)中对比了采用ANSYS与TBFEA的整个静强度分析过程的总时间，在50万网格规模以上加速比基本达到了4倍左右，在4000k网格数量时，采用ANSYS计算时间为3161s，采用TBFEA的计算时间为759s，加速比高达4.16，本方法明显提高了计算效率，减少了总的计算时间。

实例二

对于某带有围带/拉筋结构的长叶片模型，主体采用八节点六面体单元，在叶根和拉筋过渡部分采用四面体以及退化单元进行网格划分，单元共具有206292个单元，216798个节点，图6中(a)给出了长叶片的整体有限元网格模型，图6的(b),(c),(d)给出了叶根平台/拉筋/围带局部的网格剖分状况。本实例不考虑轮缘部分，将枞树型叶根四个承载齿面全固定，对拉筋围带部分的接触面施加100N的正压力，模拟拉筋/围带接触面因接触产生的压力，最后施加工作转速为3000r/min的离心力载荷。

表2长叶片各部分最大Von-Mises应力状况汇总

本实例的考核部位为：叶根平台圆角，围带部分圆角，拉筋部分圆角；从计算获得的等效应力分布云图(图7)以及叶片各考核部分最大Von-Mises等效应力对比(表2)可以看出，TBFEA与ANSYS计算结果基本一致，达到了工程计算的精度。对比ANSYS与TBFEA的计算时间，TBFEA程序计算时间为246s，ANSYS计算时间为1459s，加速比可以达到6倍左右，针对真实叶片模型，本方法显著地缩短了静强度分析计算时间。

以上所述，仅是本发明的较佳实施例，并非对本发明作任何限制，凡是根据本发明技术实质对以上实施例所作的任何简单修改、变更以及等效结构变换，均仍属于本发明技术方案的保护范围内。

Claims (6)

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

1)有限元模型建立：对于给定的透平机械叶片的三维模型和材料参数，在有限元网格剖分软件中建立透平机械叶片的有限元模型，将有限元模型的节点坐标与单元-节点对应关系存储在文件中，将叶根或轮缘要设置约束条件的节点以及要耦合的节点分别划分为集合，将要考核变形和应力部位的节点或单元输出到文件中，将透平机械叶片稳定工作时的气动数据记录在文件中；将透平机械叶片的预应力场数据记录在文件中，将透平机械叶片的材料参数以及工况参数输出到文件中；

2)计算总刚度矩阵[K]：读取步骤1)中透平机械叶片的有限元模型的节点和单元数据，先对单元进行染色，然后计算总刚度矩阵[K]的行索引和列索引向量以及在压缩存储中的索引值，再根据染色结果在CPU和GPU上同时计算所有单元刚度矩阵并组装为总刚度矩阵[K]；

3)计算总载荷向量{F}：根据步骤1)中存储的透平机械叶片的材料参数和工况参数，分别计算透平机械叶片因机械旋转引起的离心载荷向量以及因气流冲击导致的气动载荷向量，并合成为总载荷向量{F}；

4)边界条件处理：根据步骤1)中存储的叶根或轮缘要设置约束条件的节点以及要耦合的节点的集合，对该集合中的节点进行位移约束和耦合，并修正步骤2)中获得的总刚度矩阵和总载荷向量；

5)方程组求解，用CPU+GPU并行求解步骤4)修正后的总刚度矩阵和总载荷向量形成的方程组，获得节点位移向量；

6)应力结果处理：先将步骤5)获得的节点位移向量转化为单元高斯积分点上的应变和应力，并外插获得单元节点上的应变和应力，计算主应变和VonMises等效应力，再将单元节点的位移、应变和应力输出，绘制分布云图；

7)安全考核：统计步骤1)中考核变形和应力部位的节点的位移、主应变和Von-Mises等效应力，将其中的最大位移，最大主应变，最大等效应力值与叶片材料安全承载极限做对比，并输出考核结果，即完成了透平机械叶片静强度特性分析。

2.根据权利要求1所述的基于CPU+GPU异构并行计算的透平机械叶片静强度特性分析方法，其特征在于，所述的步骤2)具体包括以下步骤：

2-1)根据步骤1)中文件存储的单元-节点对应关系，建立单元之间的连通矩阵[L]，其中L_ij表示第i个单元的第j个相邻单元编号；所述的相邻单元为有公共节点的单元；

2-2)依次遍历所有单元，根据连通矩阵[L]对单元进行染色分类，当两个单元相邻时，染上不同的颜色，获得颜色的个数以及每种颜色的单元编号；

2-3)根据步骤1)中文件存储的单元-节点对应关系，建立节点之间的关系矩阵[R]，其中R_ij表示与第i个节点相关的第j个节点的编号，并将每个节点的相关节点的编号按照从小到大排序，统计关系矩阵[R]中每个节点对应的相关节点的个数，组成相关和向量{M}，其中M_i表示第i个节点对应的相关节点总数；所述的相关节点为位于同一个单元的节点；

2-4)计算总刚度矩阵[K]的行索引向量和列索引向量以及在压缩存储中的索引值，对于第i个节点相关的第j个节点的刚度矩阵元素，其在总刚度矩阵的行索引IndexCol和列索引IndexRow分别为：

$IndexCol=\left[\begin{array}{ccc}3(i-1)+1& 3(i-1)+1& 3(i-1)+1\\ 3(i-1)+2& 3(i-1)+2& 3(i-1)+2\\ 3(i-1)+3& 3(i-1)+3& 3(i-1)+3\end{array}\right]$

$IndexRow=\left[\begin{array}{ccc}3(R_{ij}-1)+1& 3(R_{ij}-1)+2& 3(R_{ij}-1)+3\\ 3(R_{ij}-1)+1& 3(R_{ij}-1)+2& 3(R_{ij}-1)+3\\ 3(R_{ij}-1)+1& 3(R_{ij}-1)+2& 3(R_{ij}-1)+3\end{array}\right]$

在压缩存储中的索引值IndexVector为：

$IndexVector=\left[\begin{array}{ccc}3U+3(j-1)+1& 3U+3(j-1)+2& 3U+3(j-1)+3\\ 3U+M_i+3(j-1)+1& 3U+M_i+3(j-1)+2& 3U+M_i+3(j-1)+3\\ 3U+2M_i+3(j-1)+1& 3U+2M_i+3(j-1)+2& 3U+2M_i+3(j-1)+3\end{array}\right]$

其中U＝M₁+……+M_i-1；

2-5)选取一个颜色的单元集合，将该集合中单元分别分配给CPU和GPU，在CPU和GPU上分别计算各自分配到的单元的刚度矩阵；CPU上使用OpenMP进行多核并行计算，GPU上编写CUDA Kernel函数并行计算，计算方法如下：

A.进行网格质量检查，从步骤1)中存储的文件中分别读取CPU和GPU分配到的单元内的节点坐标，计算所述单元在所有高斯节点上的Jaccobi值，如果小于等于0，则退出，返回步骤1)重新进行有限元模型建立，否则进行下一步；

B.遍历CPU和GPU分配到的单元中的所有节点，计算几何子矩阵[B_L]_i，设N_i为第i个节点的插值函数，用表示对第k个局部坐标的偏微分，则第i个节点对应的几何子矩阵[B_L]_i为：

${\[B_L\]}_i={\left[\begin{array}{cccccc}{\∂}_1{N}_i& & & & {\∂}_3{N}_i& {\∂}_2{N}_i\\ & {\∂}_2{N}_i& & {\∂}_3{N}_i& & {\∂}_1{N}_i\\ & & {\∂}_3{N}_i& {\∂}_2{N}_i& {\∂}_1{N}_i& \end{array}\right]}^T={\[B_{L1i},B_{L2i}\]}^T$

其中

C.计算第i个节点相关的第j个节点的线性刚度子矩阵[K_L]_ij和预应力刚度子矩阵[K_σ]_ij，[K_L]_ij＝[B_L1i ^TEB_L1j+GB_L2i ^TB_L2j]，

其中，[I]为3阶单位矩阵，σ＝[σ_x,σ_y,σ_z,τ_xy,τ_xz,τ_yz]为预应力张量，E为叶片材料参数中的弹性模量，G为叶片材料参数中的剪切模量；

D.计算刚度子矩阵[K]_ij＝[K_L]_ij+[K_σ]_ij，并将所有的高斯节点位置的刚度子矩阵[K]_ij相加，获得最终的刚度子矩阵；

E.在GPU和CPU上按照索引向量IndexCol，IndexRow，IndexVector步骤将D中计算的刚度子矩阵加到总刚度矩阵的对应位置，利用MPI_WAIT()进行线程协调，直到GPU和CPU都完成分配的计算任务；

2-6)若遍历完所有颜色集合，则进行步骤2-7)，否则继续进行步骤2-5)；

2-7)将CPU计算获得的刚度矩阵和GPU获得的刚度矩阵相加为总刚度矩阵[K]。

3.根据权利要求2所述的基于CPU+GPU异构并行计算的透平机械叶片静强度特性分析方法，其特征在于，所述的步骤3)具体包括以下步骤：

3-1)选择一种颜色的单元集合，将该集合中单元分别分配给CPU和GPU，由GPU和CPU分别计算各自分配到的单元的离心载荷向量其中ρ为叶片材料参数中的材料密度，ω为透平机械叶片工况参数中的转速，{r}为单元内高斯节点到透平转轴的垂直矢量，[N]为单元型函数；

3-2)计算CPU和GPU上分配到的单元的所有节点在总载荷向量的索引IndexForce，对于第i个节点：

$IndexForce=\left[\begin{array}{c}3(i-1)+1\\ 3(i-1)+2\\ 3(i-1)+3\end{array}\right]$

将所述单元的离心载荷向量按照IndexForce的索引加入总载荷向量；

3-3)若遍历完所有颜色集合，则进行步骤3-4)，否则继续进行步骤3-2)；

3-4)选取步骤1)存储的透平机械叶片稳定工作时的气动数据中的一个单元表面S^e，用CPU计算单元气动载荷向量其中p为单元表面受到的气流冲击压强，{n}为单元表面的单位法矢量，[N]为单元型函数；

3-5)将所述单元的气动载荷向量按照IndexForce的索引加入总载荷向量；

3-6)遍历完S集合，获得总载荷向量{F}；其中S集合为透平机械叶片稳定工作时的气动数据中的所有单元表面。

4.根据权利要求3所述的基于CPU+GPU异构并行计算的透平机械叶片静强度特性分析方法，其特征在于，所述的步骤4)具体包括以下步骤：

4-1)提取步骤1)存储的叶根或轮缘要设置约束条件的节点的集合的数据，对节点被约束自由度对应刚度矩阵位置进行修正，对被约束自由度的对角线元素赋值1，并将其余行和列元素全部置0，并将总载荷向量{F}的对应位置取为约束值，获得修正后的总刚度矩阵[K]^m和修正后的总载荷向量

4-2)提取步骤1)存储的叶根或轮缘要耦合的节点的集合的数据，对于第i个节点和第j个节点的耦合节点对，其法向和两个切向耦合刚度分别为K_n，K_t1，K_t2；根据节点的坐标计算主坐标位移与耦合局部坐标位移的变换矩阵T，然后计算主坐标系下的耦合刚度矩阵其中

${\[K^{\′}\]}_{ij}=\left[\begin{array}{ccc}{K}_n& 0& 0\\ 0& K_{t1}& 0\\ 0& 0& K_{t2}\end{array}\right]$

4-3)利用二次修正总刚度矩阵[K]^m，修正后的总刚度矩阵为：

$\[\stackrel{\‾}{K}\]=\left[\begin{array}{ccccc}...& ...& ...& ...& ...\\ ...& {\[K\]}_i+{\[\stackrel{\‾}{K}\]}_{ij}& ...& -{\[K\]}_{ij}& ...\\ ...& ...& ...& ...& ...\\ ...& -{\[K\]}_{ij}& ...& {\[K\]}_j+{\[\stackrel{\‾}{K}\]}_{ij}& ...\\ ...& ...& ...& ...& ...\end{array}\right]$

其中[K]_i和[K]_j为矩阵[K]^m主对角线上的第i个节点和第j个节点对应位置的子矩阵。

5.根据权利要求4所述的基于CPU+GPU异构并行计算的透平机械叶片静强度特性分析方法，其特征在于，所述的步骤5)具体为：修正后的总刚度矩阵节点位移向量{δ}和修正后的总载荷向量构成的方程组为：

对于展弦比小于10的叶片，具体算法如下：

5-1)构造预处理矩阵[P]^-1：当展弦比小于3时，[P]^-1取修正后的总刚度矩阵的主对角线block矩阵，block的尺寸取5-10之间；当展弦比大于3小于10时，取修正后的总刚度矩阵的不完全cholesky分解；

5-2)在CPU上计算残差向量组{δ₀}为与节点位移向量{δ}维数相同的随机向量，构造预处理向量{z₀}＝[P]^-1{r₀}，并对辅助向量组{q₀}赋值，{q₀}＝{r₀}；

5-3)进行迭代循环，迭代数k＝0,1,2……

5-3-1)用CPU计算向量组内积{z_k}^T{r_k}，同时GPU上计算获得

5-3-2)用CPU计算{δ_k+1}＝{δ_k}+α_k{q_k}，同时GPU计算

5-3-3)判断||{r_k+1}||＜10^-7，满足则进入步骤d)，否则继续进行步骤5-3-4)；

5-3-4)求解预处理方程{z_k+1}＝[P]^-1{r_k+1}，当展弦比小于3时，在GPU上进行求解，当展弦比大于3时，在CPU和GPU上共同求解；

5-3-5)用CPU计算并获得{q_k+1}＝{z_k+1}+β_k+1{q_k}，返回步骤5-3-1)；

5-4)将获得的向量{δ_k+1}作为原方程组的解向量{δ}；

对于展弦比大于10的叶片，采用CPU+GPU异构式的稀疏矩阵的高斯消元法求解。

6.根据权利要求5所述的基于CPU+GPU异构并行计算的透平机械叶片静强度特性分析方法，其特征在于，所述的步骤6)具体包括以下步骤：

6-1)根据步骤5)中获得的节点位移向量{δ}，选择一种颜色的单元集合，同时由GPU和CPU计算该集合内单元的高斯节点上的应变分量{E}^G和应力分量{S}^G，{E}^G＝[B]^e{δ}^e，{S}^G＝[D][B]^e{δ}^e，

$\[D\]=\left[\begin{array}{cccccc}\frac{E(1-\mathrm{\μ})}{(1+\mathrm{\μ})(1-2\mathrm{\μ})}& E-2G& E-2G& & & \\ E-2G& \frac{E(1-\mathrm{\μ})}{(1+\mathrm{\μ})(1-2\mathrm{\μ})}& E-2G& & & \\ E-2G& E-2G& \frac{E(1-\mathrm{\μ})}{(1+\mathrm{\μ})(1-2\mathrm{\μ})}& & & \\ & & & G& & \\ & & & & G& \\ & & & & & G\end{array}\right]$

其中E为叶片材料参数中的弹性模量，G为剪切模量，μ为泊松比，[B]^e为单元的几何矩阵，{δ}^e为单元的位移向量；

6-2)计算该单元内所有高斯节点的插值函数方阵[NG]，对每一个应变分量{E}^G和应力分量{S}^G分别外插，获得单元节点上的应变分量{E}^e和应力分量{S}^e，{E}^e＝[NG]^-1{E}^G，{S}^e＝[NG]^-1{S}^G；

6-3)依照步骤6-1)和6-2)对所有单元内的所有节点进行应力分量和应变分量的计算，对重复计算的节点的应力分量和应变分量分别取算术平均值，获得所有节点的实际应力分量{S}和实际应变分量{E}；

6-4)根据每个单元节点或高斯节点的实际应力分量和实际应变分量计算三个主应力S₁，S₂，S₃和三个主应变E₁，E₂，E₃；对每个节点计算其VonMises等效应力，VonMises等效应力的计算方法为

6-5)根据节点位移向量{δ}、实际应力分量{S}、实际应变分量{E}和VonMises等效应力，绘制叶片的位移、应变、应力分布云图以及变形图。