水轮机叶片毕业设计

一、工程背景及水轮机叶片简介

图 1、为某型水轮机叶片的CAD模型。在发电工作工程中水流由进水口流向出水口，叶片承受水流的冲刷从而开始运动，这种运动通过传动轴传递到发电机，从而带动发电机工作发电。但是水轮机在工作仅仅一年多时间以后，就有数片叶片发生了疲劳断裂事故，使得水轮机不能正常工作发电，造成了一定的经济损失，同时也说明水轮机叶片在结构的设计方面确实存在不完善之处。然而，由于水轮机在水下进行工作，很难通过测量得方法获得叶片上应力和位移的分布情况，也就无法知道叶片为何会断裂，无法有效的改善叶片的几何结构。在这种情况下，长江水利委员会陆水枢纽局的委托我们对LS591水轮机叶片的进行Ansys有限元模拟计算，获得叶片的应力场和位移场的分布，从而为叶片断裂事故分析提供技术支持，并对叶片结构的改进提供具体方案。

传动轴

进水口出水口

图1、CAD模型

二、ANSYS简介及解题步骤

1、ANSYS简介

对于大多数工程技术问题，由于物体的几何结构比较复杂或则问题的某些特征是非线性的，我们很难求得其解析解。这类问题的解决通常具有两种途径：一是引入简化假设，但这种方法只是在有限的情况下是可行的。也正是因为这样，有限元数值模拟的技术产生了。有限元方法通过计算机程序在工程中得到了广泛的应用。到80年代初期，国际上较大型的面向工程的有限元通用软件达到了几百种，其中著名的有：ANSYS,NASTRAN,ASKA, ADINA,SAP等。其中，以ANSYS为代表的工程数值模拟软件，即有限元分析软件，不断的吸取计算方法和计算机技术的最新进展，将有限元分析、计算机图形学和优化技术相结合，已成为解决现代工程问题必不可少的有力工具。尤其是在某些环境中，样机试验是不方便的或者不可能的，而利用ANSYS软件，对这个问题有了很好的解决。本文中水轮机叶片是在水下的环境进行工作，测量很难进行，利用有限元软件ANSYS这个问题得到了很好的解决。

2、ANSYS分析步骤

ANSYS分析可以分为三个步骤：

a、创建有限元模型

（1）创建或读入几何模型

根据实体模型按照给定的尺寸建立模型或者直接导

入已经生成的几何模型，并对其进行一定程度的修

复、简化等。

（2）定义单元类型，设定实常数、定义材料的属性定义单元类型：对于任何分析，必须在单元类型

库中选择一个或几个合适的单元类型，单元的类型决定了附加的自由度（位移、转角、温度等）。许多单

元还要设置一些单元的选项，诸如单元特性和假设

等。

设定实常数：有些单元的几何特性，不能仅用

其节点的位置充分表示出来时，就需要提供一些实

常数来补充几何信息。

定义材料属性：材料属性是与几何模型无关的本

构属性，例如杨氏模量、密度等。虽然材料属性并

不与单元类型联系在一起，但由于计算单元刚度矩

阵时需要材料属性，所以在此我们要对材料的属性

进行相关的定义。

（3）划分网格（节点及单元）

在做好上述的所有工作后，接下来就是对实体模型进

行网格划分，此步尤为关键，因为网格划分的好坏将

直接影响到计算结果的精确度与收敛性。根据模型的

拓扑结构决定采用映射网格还是自由网格，之后对网

格的尺寸进行设定，对关心的部位或者危险部位进行

必要的网格细化。总之,网格要足够细，才能保证结

果的精确性。

b、施加载荷并求解

（1）施加载荷及载荷选项、设定约束条件

施加约束根据具体情况对有限元模型进行约束设定

施加载荷包括集中载荷、面载荷、体载荷、惯性载荷

等。

（2）求解

选择求解器 ANSYS提供了两个直接求解器：波前求解器、稀疏矩阵求解器，同时还提供了三个迭代求解器：PCG、JCG、ICCG。因此，在前根据具体情况选择合适的求解器，这样直接影响求解的速度和结果的精确度。

进行求解

c、后处理

（1）查看结果

静态分析的结果写入结果文件，结果由以下数据构成：

基本数据——节点位移（UX,UY,UZ,ROTX,ROTY,ROTZ）

导出数据——节点单元应力、单元应变、单元集中力、节点反力等。

（2）分析结果

可直接LIST结果数据，也可通过等值线、矢量图等形式

对结果进行观察分析。

三、叶片几何模型的建立

叶片的几何模型根据相关的图纸建立，CAD模型如图1、图2所示。

图2、几何模型

我们把整个叶片离散成474个小六面体分别建立，其目的有三：

●建立模型的需要

从图2、知道叶片的几何结构不是很规则，其上下两个表面都是形式比较复杂的超曲面。基于这种情况，很难建立一个和实际叶片一模一样、丝毫不差的模型，只要把叶片离散成474个足够小的小六面体逼近实际模型，这个问题就得以解决。为了保证叶

片几何模型上下两个表面光滑，我们在建模时采用了以下的方案：

A

B

C

D

E

F

G

H

图3、建模方案

其中A,B,C,D,E,F,G,H分别为相邻小六面体的顶点，首先通过ABCD,EFGH建立两条三次样条插值曲线，然后建立曲线ED和曲线ABCD,EFGH相切，这样就使得ABCDEFGH成为一条光滑的曲线。使用这种方法，也就保证了叶片的上下两个表面都是光滑的曲面。

由于叶片的上下两个面是通过叶片上的一些离散点三次样条插值得到，所以几何模型和实际模型的逼近效果较好。

●划分高精度六面体映射网格的需要

映射网格比自由网格具有更高的计算精度。然而划分六面体映射网格对模型的拓扑结构有严格的限制：只有形状较规则的六面体和三棱柱才能划分映射网格。我们把叶片分成474个小六面体，这些小体形状不很奇异、比较规则，正好可以满足划分映射网格的要求。

●在叶片上加载的需要

接下来就要在叶片的上下两个表面上加载边界压强，但是由于叶片表面上承受的压强是随着曲面变化的，从而叶片上每一点的压强都不一样，这就为我们加载添加了困难。由于上下两个表面都被分成474个小面，我们就可以按照等效的原则把载荷平均加载到各个小面上。

四、网格的划分

1、单元模型（三维20节点单元）介绍

采用高精度的solid95单元对叶片进行离散。Solid95如图4、所示。

图4、三维20节点Solid95单元

它是三维8节点Solid单元Solid45的高阶形式，它能够容忍不规则

的形状而保持足够的精度。Solid95单元具有协调的形函数并且能够很好的模拟曲线边界，对于叶片的上下两个曲面的几何模型来说，这种单元非常合适。该单元有20个节点，每个节点有三个自由度：x, y, z 方向的位移。

一个20节点的等参单元由图5所示。在母单元中建立ξηζ坐标系，起原点在母单元的形心处，也可以将ξηζ理解为实际单元的局部坐标系。

ζ

ξ

η

(a)

(b) 图5、20节点等参单元

坐标变换式和位移模式可统一写成如下的形式：

∑==n i i i x N x 1

∑==n i i i y N y 1 ∑==n i i i z N z 1 （1）

∑==n i i i u N u 1 ∑==n i i i v N v 1 ∑==n

i i i w N w 1

（2）

式中 n —单元的节点数。

当n=8时，指的是8节点等参单元，首先写出它的形函数

8/)1)(1)(1(000ζηξ+++=i N ( i =1,2,……8) （3） 其中，ξξξi =0 ，ηηηi =0，和ζζζi =0，而i ξ，i η，i ζ是节点i 的局部坐标，对于角节点它们分别为+1和-1。

观察形函数（3），其右端的每一项正好是距节点i 距离为2的三个平面方程的函数。将其他7个角节点代入结果等于零，将节点i 代入正好等于1，因此系数八分之一是按形函数要求而确定的。依照这个办法，能写出节点9–20的各个对应的形函数：

4/)1)(1)(1(002ζηξ++-=i N （i =9，10，11，12）

4/)1)(1)(1(002ξηη++-=i N （i =13，14，15，16）

4/)1)(1)(1(002ξηζ++-=i N （i =17，18，19，20） （4） 其中，ξξξi =0 ，ηηηi =0，和ζζζi =0，对于节点9到20，i ξ，i η，i ζ分别取0或+1和-1，例如09=ξ ，19-=η，和19-=ζ。

对于20节点等参单元，其边上节点形函数如式（4）所示，其角节点由如下的线性组合来表示 ∑=++++-=12

18888),,(i i i i i i i i N N N N ζηξ （5） 其中，i N 即式（3）表示的形函数。

如果增加一个约定：在形函数（4）和（5）中令某一个形函数或某几个形函数恒等于零，即表示20节点单元由相应的一个或几个边上的节点不存在。有了这个约定，则（4）和（5）就可以表示为8-20等参单元的形函数。这种单元由实用价值。

按几何关系和式（2），应变计算公式为：

[]{}[]????????????????????????==??????????????????????+????+????+???????=202120321:::/}{δδδδεB B B B B z u x w y w z v x v y u z v y v x u e （6）

其中

[]?????????

???????????=x i z i y i z i x i y

i z i y i x i i N N N N N N N N N B ,,,,,,,,,000000000 {}??????????=i i i i w v u δ （i =1，2，…20）（7） 记号x i N ,、y i N ,、和z i N ,分别表示i N 对x 、y 和z 的偏导数。根据

复合函数求导法则，他们与ξ,i N 、η,i N 、和ζ,i N 有如下的关系式

[]??????????=???????

?????????????=??????????z i y i x i z i y i x i i i i N N N J N N N z y x z y x z y x N N N ,,,,,,,,,,,,,,,,,,ζζζηηηξξξζηξ （8） 其中

∑==n i i i x N x 1,,ξξ……∑==n

i i N z 1,,ξξ （9）

[]??

????????=??????????-ζηξ,,,1,,,i i i z i y i x i N N N J N N N （10）

应力的计算公式为：

{}[]T xz yz xy z y x τττσσσσ=

=[][]{}{}e e DB DB DB B D δδ][2021?= (11)

而

[][]?

?

?

?

??

?

?

??

??????????=x i z i y i z i x

i y i z i y i x i z

i y i x i z i y i x i i N A N A N A N A N A N A N N A N A N A N N A N A N A N A B D ,2,2,2,2,2,2,,1,1,1,,1,1,1,3000 （i =1，2，…20） （12）

其中[D]为弹性矩阵，)1(1μμ-=A )1(2)21(2μμ--=A )]21)(1[()1(3μμμ-+-=E A E 为弹性模量，μ为泊松比。

单元刚度矩阵可以分成2020?个子矩阵，典型的子矩阵是： [][][][][][][]ζηξd d d J B D B dxdydz B D B k j T

i j T i ij ??????---==111111 （13） 而

[][][]????

??????++++++++=z

j y i y j z i z j x i x j z i x j x i z j z i y j y i y j x i x j y i x j y i y j x i z j z i y j y i x j x i j T

i N N A N N A N N A N N A N N N N A N N N N A N N A N N A N N A N N N N

A N N

B D B ,,2,,1,,2,,1,,,,2,,,,2,,1,,2,,1,,,,2,,)

()(???

?

?

++++)(,,,,2,,,,2,,1,,2,,1y j y i x j x i z j z i y j z i z j y i x j z i z j x i N N N N A N N N N A N N A N N A N N A (14)

以上为20节点单元有限元基础知识的一些简要介绍。

2、网格模型

我们考虑了两种有限元网格划分。其中，较粗的一种网格由3264个20节点的SOLID95单元组成，沿叶片厚度方向分为4层（如图5）；较细的一种网格由11592个错误！链接无效。组成，沿叶片厚度方向分为5层，为了更真实反映断裂部位的应力集中，在实际断裂部位，有限元网格进行了必要的加密（如图6）。

图6、有限元网格之一（叶片部分含3260个单元）

图7、有限元网格之二（叶片部分含11592个单元）

五、约束和外载

根据实际工况，我们把传动轴x,y,z三个方向的位移全都约束为0。叶片在水下工作，其上下两个面都受到压强的作用，但是我们无法知道叶片面上每个点的压强，只能够通过数值的方法计算出各个

离散点上压强的大小。由于我们在建立模型的时候整个叶片分成了

474个小六面体，因此上下两个面也就各被分为474个小面，我们认为这些小面足够的小，从叶片上取出一个典型的小面（如图8），面上的压强可以通过下面这个式子差值得到： 4

l k j i p p p p p +++=面 （15）

i j

l k

图8、叶片上的一个面

其中i p j p k p l p 为小面各个顶点处的压强（通过流体力学数值

方法求得）。用这种差值方法计算出各个小面上的压强，然后我们使用FORTRAN 编制了一个加载小程序产生各个小面上的加载命令行： SFA, 1805, 1, pres, 250380

SFA, 1793, 1, pres, 214049

SFA, 1766, 1, pres, 124498

………

其中SFA 是加面载荷的命令，1805表示要加载的面在ANSYS 中的编号，1,pres 表示加压强，250380表示压强的大小，单位为国际标准单位Pa 。然后，Ansys 读入这些加载命令行 ，外加载荷就加到叶片上了。约束和外载都加上去了后，如图9所示：

图9、加载后的模型

六、计算及结果分析

由于叶片采用的是一种各向同性的铝合金材料，根据金属材料手册，叶片材料的机械性能参数如表一所示：

表一：材料参数表

我们计算了三种不同的工况，即三种不同的压力边界条件（由

清华大学提供）。由于第一种工况是最危险的工况，而且另外两种工

况下的应力场和位移场和第一种工况除了在数量上有差别外，其分布

特点基本和第一种工况下一致。因此，本文把第一种工况作为代表工

况，分析叶片的刚度和强度。其他两种工况下的计算结果在结论中给

出，其具体分析和第一种工况一样，这里就不再赘述。

1、粗网格的计算结果

图10给出的是网格较粗时，在第一种工况（最危险工况）下

的有限元计算结果。从图中可见：

(1)最大等效应力水平为167MPa，叶片的最大等效应力发生在

实际断裂部位,这和实际工作中出现的情况吻合得很好。但

似乎在这么低的应力水平下不会发生在实际工作中如此快

的断裂破坏，所以，我们有必要对网格进行细化以求得更加

精确的解。

(2)叶片沿Z方向的最大位移分别为1.82mm，最大位移发生在

叶片的外缘底部。

(a)叶片上的等效应力分布

(b)叶片上Z向的位移分布

图10、有限元计算结果（对应于粗网格）

2、网格较细时的计算结果

由于传动轴和叶片连接的部位，在几何形式上变化比较急剧，这就相当于轴上的台阶，在这种地方特别容易产生应力集中。粗网格的计算结果中最大等效应力发生在这个部位就很好的说明了这一点。所以我们有必要使用第二种在此处加密了的网格进行计算（为了反映应力在叶片厚度方向的梯度变化，在叶片厚度方向也加了一层网格）。图11给出的是网格较细时，最危险工况下的有限元计算结果。从图可见：

（1）叶片的最大等效应力也发生在实际断裂部位，最大等效应力水平约为556MPa，已超过材料的屈服应力550MPa。我们发现加

密后的计算结果，最大等效应力比粗网格时增加了2倍左右，在这么大的应力水平下发生很快发生断裂破坏是合理的。同时

在叶片断裂部位附近，等效应力下降的梯度也明显变大，表明

该处应力集中较明显，易于发生疲劳破坏。同时这表明，计算

时对该部位进行网格加密细分是十分必要的。网格加密后的计

算结果更精确，更加符合实际情况。

（2）叶片外缘沿Z方向的最大位移分别为1.835mm。对比图10-图11可知,网格的细分加密对位移的计算结果影响不大（不加密时，Z方向的最大位移分别为1.82mm）。

（a）叶片上的等效应力分布（对应于细网格）

（b）叶片上Z向的位移分布

图11、有限元计算结果（对应于粗网格）

3、应力水平的降低

通过有限元计算我们知道叶片的最大等效应力基本上处于屈服应力状态，在这么高的应力水平循环作用下叶片容易发生疲劳断裂。下面我们就要想办法改进叶片的结构使得应力水平降低，避免断裂破坏事故的过早发生。材料力学和机械设计有关书籍指出：可以通过在

机械零件的台阶处加卸载槽或者过渡圆角来降低应力集中系数。由于

传动轴要和其他零件相互配合，所以在这里使用过渡圆角将会对传动轴的配合产生不利的影响，考虑到这个因素我们使用卸载槽来缓解应力集中。应力集中系数和卸载槽的半径有关,理论上半径大到一定的程度，应力集中系数就会不再下降了。而半径越大，对叶片的刚度削弱的就越厉害，叶片的位移也就越大，对水轮机的工作性能就会产生不好的影响。所以我们需要找到一个最佳半径，使得应力水平有较大的降低而不使刚度明显下降。我们作了r=5mm r=10mm r=15mm三种半径的卸载槽，分析其应力水平和位移场，从而得到我们需要的最佳卸载槽半径。

（1）加R=5mm卸载槽的计算结果

图12给出的是在叶片应力集中处加R=5mm卸载槽时的有限元计算结果。从图中可见：

●叶片的最大等效应力有一定的降低，叶片上的最大等效应力由

556MPa降到352MPa,约为屈服应力550MPa的3左右，此时应

力水平仍然很高容易发生断裂破坏。

●在叶片断裂部位附近，等效应力下降的梯度有所下降，但应力

集中依然比较显著，这表明加卸载槽对于降低应力集中有一定

的效果，同时说明R=5mm还不能达到我们的要求。

叶片外缘沿Z方向的最大位移分别由不加卸载槽时的1.835mm、变为 1.849mm。由此可知,加入卸载槽使得叶片的刚度有所下

降，但最大位移只增加了0.014mm。

（a）叶片上的等效应力分布(倒角r=5mm)

（b）叶片上Z向的位移分布

图12、有限元的计算结果（对应倒角r=5mm）

（2）加R=10mm卸载槽的计算结果

图13给出的是在叶片应力集中处加R=10mm卸载槽时的有限元计算结果。从图可见：

●叶片的最大等效应力有显著的降低，叶片上的最大等效应力由

556MPa降到261MPa,约为屈服应力550MPa的一半。R=10mm比

R=5mm的卸载槽对应力水平又降低了91Mpa,效果要更好。

●在叶片断裂部位附近，等效应力的变化梯度下降明显，此时应

力分布很均匀，应力集中得到了很好的缓解。表明加R=10mm 卸载槽对于降低应力集中可以产生很好的效果。

叶片外缘沿Z方向的最大位移由不加倒角时的1.835mm变为

1.879mm。位移的计算结果比不加倒角增加了0.044mm。叶片的

刚度没有明显的降低。

（a）叶片上的等效应力分布(倒角r=10mm)

（b）叶片上Z向的位移分布

图13、有限元的计算结果（对应倒角r=10mm）

（3）加R=15mm卸载槽的计算结果

图14给出的是在叶片应力集中处加R=15mm卸载槽时的有限元计算结果。从图中可见：