叶片有限元分析中弹塑性过渡区应力奇异产生原因及解决方法

仲继泽,徐自力,方宇,范小平,赵世全



叶片有限元分析中弹塑性过渡区应力奇异产生原因及解决方法

仲继泽1,2,徐自力1,2,方宇3,范小平3,赵世全3

叶片在高温、高压蒸汽的推动下驱动汽轮机转子转动,将蒸汽的热能转化为机械能。为了提高单机的出力和机组的效率,长叶片和超长叶片不断应用于汽轮机末级,例如:国内研发的全转速钢制1 092 mm叶片已应用在1 000 MW汽轮机中[1],并开发出了全转速钢制1 200 mm叶片;日本三菱重工开发了转速为3 000 r/min的1 524 mm钢制长叶片[2]。叶片长度增加,加上转速高,叶片将承受更大的离心力,叶片根部会产生大应力,在应力集中部位甚至产生了塑性变形[3-4],由此降低了叶片疲劳寿命,且容易发生事故。因此,准确计算叶片,特别是叶根的应力水平,是叶片设计中至关重要的环节。

文献[5]研究了叶根接触应力计算的弹性有限元模型,即先用稀疏网格对整个叶片进行计算,然后对叶根接触面局部进行第2次加密计算。文献[6]研究了叶根和轮缘接触边界应力的典型特征。文献[7]采用弹性有限元方法对某汽轮机末级长叶片进行了优化设计,使背弧侧叶根圆角处最大应力由1 241 MPa下降到1 172 MPa,但该叶根圆实际上角最大应力位置已经发生屈服,进入塑性状态。文献[8]采用弹性和弹塑性有限元方法对某叶片进行了应力分析,得出弹性计算值和弹塑性计算值差别较大。文献[9]发现采用弹塑性有限元法进行寿命评估更加符合实际。综上所述,当叶片应力较大且局部进入塑性时,要准确计算叶片的应力必须采用弹塑性计算方法。采用弹塑性有限元法计算多个汽轮机末级长叶片应力时,在叶片根部弹塑性过渡区的应力计算值有时会大于塑性区的应力计算值,即应力奇异现象。这种应力奇异现象对叶片的安全评价会造成困扰,影响叶片疲劳、寿命评估的准确性。

针对叶片有限元分析中在弹塑性过渡区产生应力奇异的问题,本文研究了有限元法单元节点应力的计算过程,发现有限元法通常采用高斯积分点应力外推插值得到单元节点应力,当单元一部分位于弹性区、另一部分位于塑性区时,这种外插算法就会导致节点应力计算值大于结构的实际值,甚至超出弹塑性材料的屈服极限,使得节点应力计算值出现奇异,同时给出了消除叶片根部弹塑性过渡区应力奇异现象的方法。

1 弹塑性过渡区应力奇异现象

某叶片根部的三维实体及有限元网格见图1。假定材料为理想弹塑性的,弹性模量为2.1×105MPa,泊松比为0.3,屈服极限为800 MPa,密度为7 850 kg/m3,材料的应力应变曲线见图2。在叶根齿承载面(见图1b中的加黑部位)上,施加分布力使叶片在离心力的作用下满足力和力矩平衡条件。由于该叶片离心力近4.9×106N,在离心力作用下叶片背弧侧、内弧侧叶根圆角处部分区域的应力达到材料屈服极限,进入塑性区。

(a)三维实体

(b)有限元网格图1 叶根三维实体及有限元网格

图2 叶片材料应力-应变曲线

采用三维弹塑性有限元法计算了叶片的应力,叶片背弧侧叶根圆角处应力沿轴向的分布见图3。从图3中可以看出,在弹性和塑性过渡的2个区域出现了2个类似“猫耳朵”的突起,其所表示的应力明显大于理想弹塑性材料的屈服极限800 MPa,显然是不合理的,即出现了应力奇异现象。

图3 叶根圆角处应力沿轴向的分布

2 弹塑性过渡区应力奇异产生的原因

采用有限元法对叶片结构应力场进行分析时,通常先根据结构静力平衡的有限元方程,计算得到叶片上所有单元节点的位移,然后基于位移场计算得到各单元高斯积分点上的应力,最后将高斯积分点的应力进行外推插值得到单元节点的应力[10]。本文将以六面体单元为例说明产生应力奇异的原因。

图4 六面体单元及其高斯积分点示意

(1)

式中:Nij为形函数

(2)

不失一般性,先设定一个简单的应力场。该应力场包含塑性区和弹性区,并假设:弹性区与塑性区的交界面为平面;弹性区应力梯度大小(k=32 MPa)保持不变;塑性区为理想塑性,应力梯度为0,应力σs=800 MPa。那么,弹性区内任意一点P的应力

(3)

为方便分析,以弹性区应力梯度的方向由单元坐标系原点o指向单元节点1为例,研究了单元在应力场中的3种情况:单元整体位于弹性区;单元被弹塑性交界面分割成两部分,一部分位于弹性区,另一部分位于塑性区;单元整体位于塑性区。

2.1 单元整体位于弹性区

单元整体位于弹性区时,根据单元与弹塑性交界面的距离可以分成多种情况。为方便分析,以节点1位于交界面上为例,比较节点的应力有限元值和实际应力值。

在应力场中,根据单元节点到弹塑性交界面的距离、采用式(3)可计算得到单元节点的实际应力,见表1。根据高斯积分点到弹塑性交界面的距离,先采用式(3)计算出高斯积分点的实际应力,然后将高斯积分点的应力代入到式(1)中,可计算得到单元节点应力的有限元计算值,见表1。

表1 位于弹性区时单元节点应力

从表1可以看出,单元各个节点应力的有限元计算值与单元节点的实际应力值相等,在此情况下单元节点应力没有产生奇异现象。

2.2 单元跨过弹塑性交界面

单元一部分位于弹性区、另一部分位于塑性区,根据两部分高斯积分点数及单元节点数又可以分成很多种情况。为了简洁说明单元节点应力奇异的问题,结合式(3)、式(1)计算对以下4种情况进行分析。

(1)弹塑性交界面过高斯积分点Ⅰ。此时,只有单元节点1在塑性区,节点应力计算结果见表2。

表2 交界面过高斯积分点I时单元节点应力

从表2可以看出,采用有限元插值公式计算出的单元节点1的应力超出了实际应力40.7 MPa,显然是不合理的,即出现了应力奇异情况,这是由有限元节点应力的插值方式造成的。

(2)弹塑性交界面过高斯积分点Ⅱ、Ⅳ、Ⅴ。此时,单元节点1、2、4、5和高斯积分点Ⅰ在塑性区,单元高斯积分点Ⅱ、Ⅳ、Ⅴ在弹塑性交界面上,其余节点和高斯点在弹性区。节点应力计算结果见表3。

表3 交界面过高斯积分点Ⅱ、Ⅳ、Ⅴ时单元节点应力

从表3可以看出,节点1应力的有限元计算值比实际应力低16.7 MPa,单元节点2、4、5应力的有限元计算值超出实际应力38.8 MPa。这些计算结果是不合理的,即出现了应力奇异,这也是由有限元节点应力的插值方式造成的。

(3)弹塑性交界面过高斯积分点Ⅲ、Ⅵ、Ⅷ。此时,单元节点1、2、4、5和高斯积分点Ⅰ、Ⅱ、Ⅳ、Ⅴ在塑性区,高斯积分点Ⅲ、Ⅵ、Ⅷ在弹塑性交界面上,其余高斯积分点和节点在弹性区。节点应力计算结果见表4。

表4 交界面过高斯积分点Ⅲ、Ⅵ、Ⅷ时单元节点应力

从表4可以看出:单元节点2、4、5应力的有限元计算值比实际应力小6.8 MPa;单元节点3、6、8的应力有限元计算值超出实际应力38.8 MPa,应力奇异最为严重;单元节点7应力的有限元计算值比实际应力小16.8 MPa;单元节点1的应力没有出现奇异。

(4)弹塑性交界面过高斯积分点Ⅶ。此时,只有单元节点7位于弹性区,所有高斯积分点和其余节点在塑性区。节点应力计算结果见表5。

表5 交界面过高斯积分点Ⅶ时单元节点应力

从表5可以看出:单元节点7应力的有限元计算值超出实际应力40.6 MPa,产生应力奇异;其他节点应力没有出现奇异。

2.3 单元整体位于塑性区

全部节点和高斯积分点都在塑性区。单元节点应力的有限元值和实际应力均等于屈服应力,单元节点应力没有产生奇异。

综上所述,单元整体处于弹性区或塑性区时,采用有限元法计算单元节点应力不会产生应力奇异现象。当单元跨过弹塑性交界面,即一部分处于弹性区域内、另一部分处于塑性区域内时,单元节点应力计算值会超出实际应力值,此时产生了应力奇异现象。分析认为,应力奇异是由有限元法计算单元节点应力时采用的外推插值算法造成的。

3 弹塑性过渡区应力奇异解决方法

在有限元方法中,不管是在弹性区和塑性区,还是在弹塑性过渡区,高斯积分点应力的有限元计算值都不会产生奇异。采用单元节点附近不同单元内的高斯积分点应力的加权平均值作为单元节点的应力,可以消除单元节点应力奇异现象。加权平均方法计算单元节点应力的计算式为

(4)

以六面体单元的高斯积分点I位于弹塑性交界面上为例,采用有限元节点应力公式计算得到节点1的应力为840.7 MPa,超出节点实际应力800 MPa,即产生了应力奇异。根据高斯积分点到弹塑性交界面的距离,采用式(3)计算出节点1附近高斯积分点的应力,然后将各个高斯积分点应力代入到式(4)中,计算出节点1的应力为800 MPa,与节点1的实际应力800 MPa一致,即节点1应力的计算值没有发生奇异。

以图1中所示的叶根为例对叶片进行应力分析,得到叶片根部圆角处应力的分布曲线,见图5。从图5可以看出,在叶片根部弹塑性过渡区,采用加权平均方法计算出的应力未出现奇异现象。因此,对叶片进行弹塑性有限元分析时节点应力可采用与节点相邻的高斯积分点应力的加权平均进行计算。

图5 采用加权平均方法得到的应力沿轴向的分布

4 结 论

本文基于理想弹塑性模型,对叶片有限元分析中根部弹塑性过渡区应力奇异的问题进行了研究。研究表明,当单元一部分位于弹性区、另一部分位于塑性区时,有限元节点应力的外插算法会造成节点应力计算值高于结构的实际应力,甚至超出理想弹塑性材料的屈服极限,使应力产生奇异。研究也发现,在叶片有限元分析中,采用节点相邻的高斯积分点应力的加权平均方法计算节点应力,能有效避免叶片弹塑性过渡区的应力奇异现象。该研究可为大型汽轮机、重型燃机及航空发动机叶片强度设计提供参考。

[1] 王为民, 潘家成. 东方-日立型超超临界1 000 MW汽轮机可靠性设计 [J]. 热力透平, 2006, 35(1): 8-12. WANG Weimi, PAN Jiacheng. Reliability of ultra-supercritical 1 000 MW steam turbine designed by DFSTW and Hitachi [J]. Thermal Turbine, 2006, 35(1): 8-12.

[2] HISASHI F, HIROHARU O, TOSHIHIRO M, et al. Development of 3,600-rpm 50-inch/3,000-rpm 60-inch ultra-long exhaust end blades [J]. Mitsubishi Heavy Industries Technical Review, 2009, 46(2): 18-25.

[3] RAO J S, SURESH S. Blade root shape optimization [J]. Proceedings of Future of Gas Turbine Technology, 2006, 13(2): 1-11.

[4] 徐自力, 李辛毅, 安宁, 等. 相关参数对汽轮机低压叶片疲劳寿命的影响 [J]. 动力工程, 2004, 24(1): 33-36. XU Zili, LI Xinyi, AN Ning, et al. Effect of relative factors on fatigue life for low pressure blades of steam turbines [J]. Power Engineering, 2004, 24(1): 33-36.

[5] SINCLAIR G B, CORMIER N G, GRIFFIN J H, et al. Contact stresses in dovetail attachments: finite element modeling [J]. ASME Journal of Engineering for Gas Turbines and Power, 2002, 124: 182-189.

[6] SINCLAIR G B, CORMIER N G. Contact stresses in dovetail attachments: alleviation via precision crowning [J]. ASME Journal of Engineering for Gas Turbines and Power, 2003, 125: 1033-1041.

[7] 李彬, 宋立明, 李军, 等. 长叶片透平级多学科多目标优化设计 [J]. 西安交通大学学报, 2014, 48(1): 1-5. LI Bin, SONG Liming, LI Jun, et al. Multidisciplinary and multi objective optimization design of long blade turbine stage [J]. Journal of Xi’an Jiaotong University, 2014, 48(1): 1-5.

[8] RAO J S, KISHORE C B, MAHADEVAPPA V. Weight optimization of turbine blades [C]∥Proceedings of the 12th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery. Honolulu, Hawaii, USA: Altair Engineering Inc., 2008: 1-19.

[10]ZIENKIEWICZ O C, TAYLOR R L, ZHU J Z. The finite element method: its basis and fundamentals [M]. 6th Edition. Oxford, UK: Butterworth-Heinemann, 2005: 462-467.

[11]COOK R D, MALKUS D S, PLESHA M E, et al. Concepts and applications of finite element analysis [M]. 4th ed. New York, USA: Wiley, 2001: 198-200.

(编辑 苗凌)

(1.西安交通大学航天航空学院,710049,西安;2.西安交通大学机械结构强度与振动国家重点实验室,710049,西安;3.东方汽轮机有限公司,618000,四川德阳)

采用弹塑性有限元法、借助大型商业有限元软件对汽轮机叶片进行应力分析时,弹塑性过渡区应力的计算值有时会高于塑性区应力的计算值,即会产生应力奇异现象。为分析产生这一现象的原因,以8节点六面体单元为例,研究了有限元法计算应力的过程,并在理想弹塑性的条件下,采用有限元法和解析法计算了弹塑性过渡区单元节点应力。研究发现,有限元法通常采用高斯积分点应力值外推插值法得到单元节点应力,当单元一部分位于弹性区、另一部分位于塑性区时,这种外插算法会导致节点应力计算值高于结构的实际应力,甚至超出理想弹塑性材料的屈服极限,从而造成应力奇异。研究表明,在叶片弹塑性的有限元分析中,采用相邻高斯积分点应力加权平均的方法计算单元节点应力,可有效避免弹塑性过渡区应力产生奇异的现象。

叶片;有限元;弹塑性过渡区;应力奇异

Origin and Elimination of Stress Singularity in Blade Elasto-Plastic Transition Region in Finite Element Analysis

ZHONG Jize1,2,XU Zili1,2,FANG Yu3,FAN Xiaoping3,ZHAO Shiquan3

(1. School of Aerospace, Xi’an Jiaotong University, Xi’an 710049, China; 2. State Key Laboratory for Strength and Vibration of Mechanical Structures, Xi’an Jiaotong University, Xi’an 710049, China; 3. Dongfang Turbine Co., Ltd., Deyang, Sichuan 618000, China)

When stresses of steam turbine blades are computed with the finite element method, sometimes the calculated stress in elasto-plastic transition region gets higher than the calculated stress in plastic region, namely there exits stress singularity. To explain this fact, the calculation process of nodal stress with finite element method is discussed in detail, where a kind of 8-node hexahedron element and an ideal elasto-plastic material model are chosen as the example. Nodal stresses of the elements in elasto-plastic transition region are comparatively calculated. It is found that in the finite element method the nodal stress is calculated via extrapolation of stresses at Gauss integration point in an element, and the calculated nodal stress maybe exceed the actual stress even the yield limit of ideal elasto-plastic material when one part of an element is located in elastic zone and the other part remains in plastic zone, thus stress singularity is brought out by the extrapolation algorithm. It is suggested that the nodal stresses are calculated in terms of weighted average stress at Gauss integration points of neighboring elements to effectively eliminate stress singularity.

blade; finite element method; elastic-plastic transition region; stress singularity

2014-12-17。 作者简介:仲继泽(1988—),男,博士生;徐自力(通信作者),男,教授,博士生导师。 基金项目:国家“973计划”资助项目(2011CB706505);国家自然科学基金资助项目(51275385)。

时间:2015-06-17

http:∥www.cnki.net/kcms/detail/61.1069.T.20150617.0902.010.html

10.7652/xjtuxb201509009

TK263.3;TB125

A

0253-987X(2015)09-0047-05