基于cfd的水轮机叶片动应力分析

基于cfd的水轮机叶片动应力分析

在预测设计（cd）和计算流态动力学（cfd）等技术的支持下，该工艺的功率性能和抗空蚀性指标显著提高，混合流量模型的最大效率为95%。而机组的动力学品质随着机组容量增加、转轮尺寸加大及比转速的提高越来越突现,运行中的振动及稳定性问题越来越多,使能量性能和抗空蚀性能相当优秀的水轮机在部分负荷工况下运行出力摆动加大、机组振动加剧、叶片试运行阶段产生裂纹。裂纹主要发生在叶片出水边靠近上冠和下环的过渡处。大型水电机组专家调研组通过调查分析国内水电站转轮裂纹得出:疲劳强度设计不足是转轮叶片裂纹的主要原因。因此,研究混流式水轮机转轮叶片的疲劳设计校核,应成为整个水轮机设计的一个组成部分,对解决转轮叶片裂纹问题具有重要意义。鉴此,本文通过统计分析方法,对水轮机叶片某一危险点在运行过程中产生的疲劳应力进行分析,建立了水轮机叶片应力的概率密度函数,并进行假设检验。1动应力在cfd过程中的应用水力因素引起的叶片动应力计算是进行疲劳强度校核的前提,利用CFD计算各种工况下叶片上的动压力成为叶片动应力计算的基础。利用CFD软件计算水轮机叶片的动载面力,采用FEA法计算叶片动应力。1.1材料性能参数设计及后处理模块前后处理过程包括叶片建模、参数设定、网格划分、载荷施加和求解及后处理模块。采用叶片材料性能参数:弹性模量E=2.1×1011Pa,泊松比μ=0.3,密度=7850kg/m3;单元类型为三维实体单元Solid187;网格划分采用自由网格。1.2载荷加载问题由于叶片表面为复杂的空间曲面,使用有限元计算强度时,叶片表面的载荷数据仅能通过CFD获得,但因计算要求不同,FEM和CFD分析所用的网格不同,故CFD计算所用节点不能作为FEM分析的网格节点,必须进行数据转换。这就需要实现CFD与FEM软件的数据接口,直接根据CFD软件的输出数据文件产生FEM软件的载荷输入数据文件,采用VisualC++程序编写接口程序。需解决以下问题。①叶片上施加载荷面的选取。采用Solid187单元如图1所示。由Solid187单元的属性可知,单元面和节点相互对应排列,如面I由节点J、I、K组成,节点顺序不能改变。在对叶片进行有限元网格剖分时,单元节点的次序按顺序排列。因此,选择单元时要求该单元有3个点在叶片表面上,利用ANSYS输出文件寻找对应的坐标值,查找与该点水力学数据距离最近点的坐标,再利用水力学数据将该3个点的压力找出。判断3个点组成面的序号,再用参数化设计语言(APDL)语句进行叶片面上的载荷加载。②坐标原点的重合。由于水力学数据和有限元数据的坐标原点可能不同,因此需进行坐标点的重合。选择叶片上一个特征点,比较其在水力学数据和有限元数据中的坐标位置,从而确定两个坐标系间的差别。1.3局部危险部位图2为叶片在某一工况的VonMises等效应力分布图,其最大应力分布在叶片出水边靠近上冠和下环处。其局部危险部位MX为最大等效应力分布点所在位置节点2。通过APDL参数化程序,求出叶片在某一时间段内10000个动载面力作用下的应力,再求出危险部位MX在该时间段内的等效应力,如图3所示。2叶片荷载谱的准备2.1删除最大变程阈值5%10%的点在实际应力载荷历程中,不能产生疲劳损伤的小幅载荷称无效载荷,如非峰谷值的点及幅值小于最大变程阈值5%～10%的点(通常取最大变程阈值的10%)。对非峰谷值点,若满足下式,则删除该点。(xi-xi-1)(xi-xi+1)≤0(1)式中,xi为i点的应力值;xi-1、xi+1分别为i点相邻的两个点的应力值。2.2叶片应力时间历程的雨流计数结果雨流法是考虑材料应力应变行为提出的一种计数法,认为塑性是疲劳损伤的必要条件,由应力应变的迟滞回线表现之。同时,假定一个大的变程所引起的损伤不受一个小的迟滞回线而截断的影响,可逐次将构成较小迟滞回线的较小循环从整个应变时间历程中提取重新组合。通常认为两者对材料的损伤等效。根据|xi-xi+1|≤|xi+1-xi+2|判断xi、xi+1是否满足要求,若满足,将其值取出计算应力幅值:Sa=|xi-xi+1|/2、均值Sm=|xi+xi+1|/2,同时统计频数,最后形成均值—幅值矩阵;否则读取下个数据。将叶片危险部位的应力时间历程进行雨流计数,可得到应力均值、幅值和频数信息。将应力幅值数据分为10组、均值数据分为8组,可得到循环次数矩阵N。图4为叶片危险部位MX节点2的雨流计数结果3307个全循环中均值、幅值的直方图。由图4的偏态外形,可初步认为载荷幅值服从半正态分布、载荷均值服从正态分布,且两者之间相互独立。2.3计算模型的求解使用χ2检验观察数据是否遵循某种母体分布。根据分组数据计算:D=n∑i=1∑i=1n(vi-Ei)2Ei=(vi−Ei)2Ei=j∑i=1∑i=1j(vi-npi)2npi(2)(vi−npi)2npi(2)式中,Vi为样本值落在第i区间的实际频数;Ei为理论频数;pi为母体落入第i区间的理论概率;n为总样本数;npi为随机变量落在该区间的理论频数。利用子样均值ˉXX¯¯¯和标准差S分别代替母体的均值μ和方差σ,即可确定理论频率函数:Ei=n∫ti+1titi+1tif(x,a1,a2,…,am)dx(3)式中,f(x,a1,a2,…,am)为母体的概率密度函数;a1,a2,…,am为母体的分布参数;ti、ti+1分别为第i组母体的上限和下限。将各组的理论频数Ei和实测频数vi分别代入式(2),可得D值。将D值与χ2α2α分布数值表查对应的χ2α2α(j-1-r)值比较(r为需要估计母体的参数个数),判断计算结果是否落在对应的接受区间内,计算结果大于χ2α2α值,则拒绝假设;否则接受假设。当将f(x,a1,a2,…,am)看为关于母体分布参数a1,a2,…,am的随机变量函数时,则计算D就可归纳为:minD=n∑i=1(vi-Ei)2Ei=n∑i=1(vi-n∫ti+1tif(x,a1,a2,⋯,am)dx)2n∫ti+1tif(x,a1,a2,⋯,am)dx(4)minD=∑i=1n(vi−Ei)2Ei=∑i=1n(vi−n∫ti+1tif(x,a1,a2,⋯,am)dx)2n∫ti+1tif(x,a1,a2,⋯,am)dx(4)这样当D取最小值时,确定的母体参数为参数的最佳估计量。将图4直方图均值数据分为12组,组间距22.096;幅值数据分为10组,组间距2.14。将3307个数据分配到各组的个数即为实测频数。假设均值数据服从正态分布,幅值数据服从半正态分布,其概率密度为:f(sm)=1√2πσmexp[-(sm-μm)22σ2m](5)f(sm)=12π√σmexp[−(sm−μm)22σ2m](5)f(sa)=√2√πσaexp[-(sa-μa)22σ2a]f(sa)=2√π√σaexp[−(sa−μa)22σ2a]其中μa<sa<∞(6)均值数据各组的理论频数为:Ei=n∫ti+1titi+1tif(sm)dx=3307∫ti+1ti1√2πσmexp∫ti+1ti12π−−√σmexp[-(sm-μm)22σ2m]dx(7)[−(sm−μm)22σ2m]dx(7)由于第1、2、11、12组实测频数组中数据小于5,因此需将其合并,将理论频数和实测频数代入式(2),求得D=12.497时,μm=86.076、σm=2.70。根据χ2α2α分布数值表查出χ20.0520.05(10-1-2)=14.067>D,可认为均值数据服从正态分布,接受假设。其概率密度函数为:f(sm)=12.70√2πexp[-(sm-86.076)22×2.72](8)f(sm)=12.702π√exp[−(sm−86.076)22×2.72](8)采用传统方法计算的母体参数估计值:μm=E(sm)=86.094MPaσm=√Var(sm)=2.64ΜΡaσm=Var(sm)−−−−−−−√=2.64MPaD=13.311认为3307个幅值数据服从半正态分布,采用最优化方法计算:μa=0,σa=6.254,D=11.486。采用传统方法计算的参数为:μa=0,σa=6.276,D=12.574。实测频数中,由于第10项数据小于5,需将其合并到第9项,根据分布数值表查出χ20.0520.05(9-1-2)=12.592>D,因此可认为幅值数据服从半正态分布,接受假设。采用最优化方法计算的D值比传统方法小,因此得到的母体参数为参数的最佳估计量。该叶片所受疲劳应力载荷的二维联合概率密度函数为:f(sa,sm)=1σaσmexp{-[(sm-μm)22σ2m+(sa-μa)22σ2a]}=f(sa,sm)=1σaσmexp{−[(sm−μm)22σ2m+(sa−μa)22σ2a]}=12.7×6.254exp{-[(sm-86.076)22×2.72+s2a2×6.2542]}通过叶片疲劳应力的概率密度函数及叶片疲劳强度的概率密度函数的干涉模型,就可实现水轮机叶片的疲劳可靠性校核。3疲劳应力的计数结