汽轮机末级转子叶片流固耦合计算

THERMALTURBINE

Vol󰀁36No󰀁3

Sept2007汽轮机末级转子叶片流固耦合计算

弓三伟,任丽芸,刘火星,李󰀁琳,邹正平

1

1

1

1,2

(1.北京航空航天大学能源与动力工程学院,航空发动机气动热力重点实验室,北京100083;

2.航空发动机数值仿真研究中心,北京100083)

摘󰀁要:󰀁对于大展弦比的叶片,流场与结构的相互作用不可忽略。本文利用商业软件CFX和ANSYS对某汽轮机末级转子叶片进行了设计和非设计状态下的流固耦合计算,对两种状态下耦合前后的流场结果进行了对比分析,叶片结构计算采用了双向与单向两种耦合方式。结果表明,非设计状态下流场对结构变形具有更强的敏感性;在稳态载荷作用下,两种耦合方式下的叶片位移、应力值差别不大。关键词:󰀁流固耦合;汽轮机;数值模拟;动网格

中图分类号:TK262󰀁󰀁󰀁󰀁文献标识码:A󰀁󰀁󰀁󰀁文章编号:1672-59(2007)03-0153-05

Fluid󰀂SolidInteractionofRotorBladeinLastStageofSteamTurbine

GONGSan󰀂wei,RENLi󰀂yun,LIUHuo󰀂xing,LILin,ZOUZheng󰀂ping

SmiulationResearchCenteronAero󰀂Engines,Beijing100083,China)

1

1

1

1

1,2

(1.NationalKeyLab.onAero󰀂Engines,SchoolofJetProp.,BeijingUniv.ofAero.&Astro.,Beijing100083,China;2.Numerical

Abstract:󰀁Theinteractionbetweenfluidfieldandstructureshouldnotbeignoredforthebladewithhighaspectratio.Inthispaper,fluid󰀂solidinteractionofrotorbladeinthelaststageofsteamturbineisstudiedondesignandoff󰀂designconditionsbymeansofthecommercialcodeCFXandANSYS,andtheflowfieldresultswithcoupledanduncoupledmethodsarecomparativelyanalyzedonthetwoconditions.CoupledandsinglecodeapproachesareusedinNumericalSimulationofbladestructureapplication.Theresultsindicatethat,theflowfieldonoff󰀂designconditionhasintensivesensibilitytothestructuredeformation,andthedisplacementandstressbetweencoupledandsinglecodeapproacheshavelittledifferencesunderthestableload.Keywords:󰀁

fluid󰀂solidinteraction;steamturbine;numericalsimulation;movingmesh

到全三维非定常粘性数值模拟阶段

[4]

,非线性有

0󰀁前言

在叶轮机内部,流体与结构叶片之间存在着相互作用

[1]

限元分析法已相当成熟。流固耦合理论也从早期

的线性化已经发展到了非线性阶段,对于强迫响应问题,SaymaA.I.等人

[5]

早已提出了非线性的

:一方面,叶片在气动力的作用下产流固耦合数值模型,即便理论和技术方面都有所突破,但要真正求解流固耦合问题还相当困难。通常,国内外大多数研究方法还都是通过规定叶片或叶栅的振动规律(实际上就是把叶片看成是振动着的绝对刚体,它不受流场的影响),在某种振动规律下来进行流场的非定常计算,或者通过单方面给定流体动力来分析固体的动力响应和稳定性问题,文献[6,7]基本上是采用这样的做法,文献[8]虽然同时考虑了流体与结构之间的相互作用,但是结构计算采用静强度方法,而流固耦合

生振动变形;另一方面,变形的叶片会改变流场的边界,使流动状况发生改变,进而改变叶片所受到的气动力。对于非定常流动强烈、流固耦合作用强的流体机械,将流体与结构简单解耦单独分析,可能导致分析结果与真实情况差别巨大,此时将流体与结构进行耦合计算是非常必要的。

流固耦合涉及流体和固体两个方面,在过去的十几年里,计算流体力学(CFD)和计算结构力学(CSM)都取得了长足的发展,CFD技术已发展

收稿日期:2007-04-23

作者简介:弓三伟(1983-),男,陕西武功人,北京航空航天大学能源与动力工程学院硕士生,从事叶轮机内部流固耦合研究。1汽轮机末级转子叶片流固耦合计算

[9]

实际是一种瞬态的相互作用,周新海等人对0󰀁0001s,计算结果表明这个时间步长完全满足瞬态动力分析的积分时间步长要求,当然时间步长可以取得更小一些,这样在每个物理时间步里流场计算会更容易收敛,但相应地也会加大耦合的计算步数。在每个物理时间步里,根据CFX和ANSYS的求解顺序又将双向耦合分为两种方式:同时求解和顺序求解,本文采用顺序求解方式,先求解CFX,再求解ANSYS。

1.2󰀁计算网格和边界条件

如图1所示为耦合计算的流体和叶片结构网格。通道内流体域采用H型网格,网格节点总数为39040,单元总数为33600,网格近壁区y+约为17󰀁4,满足所使用湍流模型的要求。叶片结构网格采用SOLID45单元,它是8节点的六面体单元,对于叶片结构的计算,8节点的单元已经足以满足计算的要求,叶片节点总数为4080,单元总数为26,叶片表面有2856个节点用来构成耦合交接面。

叶栅间的非定常粘性流动以及叶栅振动特性

[10]

进行了研究,袁新等人发展了一套高精度适合大尺度分离流诱发振动的流固耦合数值分析方法。

现阶段的耦合求解一般有直接耦合和间接耦合两种。本文采用商业软件CFX、ANSYS进行的流固耦合计算是间接耦合的做法,直接利用CFX进行全三维非定常粘性数值模拟,利用ANSYS进行结构瞬态动力分析,在每个物理时间步上进行耦合迭代,各自收敛后再瞬态向前推进,流场与结构的信息交换是通过ANSYS的双向耦合功能MFX(Multi󰀂fieldsolver󰀂multiplecodecoupling)来实现的,结构变形引起的流场网格位移由CFX内部的动网格技术来处理。

1󰀁耦合模型的建立

本文以某汽轮机末级为模型,一方面为了探索流固耦合的计算途径,另一方面是要考虑流场变化对结构响应的影响和结构变形对流场的影响,所以只对转子域的单通道模型进行了计算,没有考虑上游尾迹的作用。

1.1󰀁计算方法

流固耦合计算涉及流场的非定常计算和结构的瞬态动力分析。

流场非定常计算采用CFX商用软件,流动工质选取CFX软件内部的湿蒸汽模块,控制方程为非定常雷诺平均N󰀂S方程,空间采用二阶迎风格式,时间采用二阶欧拉后差格式(TheSecondOr󰀂derBackwardEulerTransientScheme),耦合计算使用动网格。非定常计算中采用了隐式双重时间步法

[11]

图1󰀁单通道耦合计算网格示意图

耦合计算前,先对整级进行定常计算,将转子域的气动参数(总温、总压和流动方向等)作为耦合计算的边界值。耦合计算过程中,流场计算是一个非定常过程,在开始计算前给定一个初场,初场是利用现有网格进行定常计算的结果。非定常计算的边界条件和定常初场计算一样:出口给定反压,进口给定总温、总压值,对整级定常计算的转子进口速度进行拟合,达到与多级环境中转子进口相一致的速度分布型。结构计算规定叶根截面的所有节点位移为零,在瞬态分析中考虑离心力的作用。

,湍流模型采用k󰀂Epsilon模型。

叶片的结构计算采用ANSYS软件中的瞬态

动力分析方法(Transient),瞬态分析的阻尼模型采用比例阻尼,为了加快耦合计算的收敛速度,适当加大了刚度阻尼的值。

流场的非定常计算和结构瞬态分析采用相同的物理时间步长,时间步长是根据叶片结构的固有频率来确定的,具体的,先对叶片进行模态分析,叶片一阶频率对应的最大周期约为0󰀁0077s,再将这个时间离散为很多步,最终取时间步长为2󰀁结果分析

耦合计算是在设计和非设计两个状态下分别第3期󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁

热力透平

155

进行的。

2.1󰀁设计状态下耦合计算结果分析2.1.1󰀁流场计算结果分析

由于耦合计算的模型是转子域的单个通道,没有上游尾迹的作用,这相当于没有考虑上游的非定常来源,所以最终将会收敛为一个稳定的流场。

耦合计算过程中,在靠近叶片近壁区的流场中取了6个监测点,分别位于压力面和吸力面的根、中、尖位置。图2给出了压力面三个监测点的网格位移曲线,计算初期短时间内,离心载荷从零加到最大,而且叶片的气动载荷并不是从零开始加载,而是以定常流场为初载荷的,所以计算初期叶片存在大幅值振动,这个会显著改变流场边界,使得叶片表面附近的网格也跟着振动起来,随着时间的推移,由于阻尼的作用,网格的振动幅值越来越小,在计算的终点时刻,网格位移趋于零,这表示网格不再运动。

󰀁󰀁图3给出了6个监测点的静压变化曲线,可以很明显地看到,在耦合计算的初期,流场的静压有一定量的波动,这是计算初期叶片大幅值振动引起的,另外,我们可以看出静压的波动量相对不是很大,这说明设计状态下叶片结构的变形对流场影响较小。随后在很长的时间内,流场的静压变化越来越小,最终趋于稳定。

图3󰀁叶片表面近壁区监测点的静压变化曲线

表1对耦合收敛流场与定常初场的计算结果进行了对比,从中可以看出,在保证膨胀比基本一致的情况下,耦合计算总温比基本不变,效率上升了1󰀁0795%,喉道附近通道面积减少了1󰀁0998%,喉道面积收缩使得出口流量减小,出口的相对马赫数也加速到更大,这些变化说明,设计状态下结构变形对流场参数有一定量的影响,但

图2󰀁叶片表面近壁区监测点的网格位移曲线图

不是很大。

效率0󰀁8847570󰀁43080󰀁0095511󰀁0795%

出口相对Ma

0󰀁9268060󰀁9352260󰀁0084200󰀁9085%

喉道附近通道面积(m2)

0󰀁04367870󰀁0431983-0󰀁0004804-1󰀁0998%

表1󰀁耦合收敛流场与定常初场计算结果对比

定常初场

耦合收敛流场绝对误差相对误差

总温比1󰀁0674331󰀁0674690󰀁0000360󰀁0034%

膨胀比

2󰀁77052󰀁7353-0󰀁000352-0󰀁0133%

流量(kgs-1)

91󰀁31691󰀁25760-0󰀁38556-0󰀁4207%

󰀁󰀁这里涉及到的参数

T1*

流量:ma,总温比:󰀁=*,膨胀比:󰀂=T

T2

P1h1-h2

=*,效率: isT

P2h1-h2

其中h表示静焓,*表示总参数,is表示等熵过程的焓值,下标1表示转子进口参数,下标2表示转子出口参数,与效率计算有关的水蒸汽热力性质参照文献[12]。

由于叶片结构变形从叶根到叶尖逐渐增大,叶根和叶中的相对马赫数分布差别较小,图4给出了定常计算和耦合计算在S1流面0󰀁9叶高位

*

*

置的相对马赫数等直线图,由于耦合计算叶片尾缘位置的网格变化最大,所以在很大程度上改变了尾迹的方向,使尾迹移向通道中部。本文的算例是汽轮机末级,如果是中间级,尾迹方向发生变

图4󰀁S1流面0󰀁9叶高位置相对马赫数分布

156汽轮机末级转子叶片流固耦合计算

化会改变下一级进气攻角,很可能会加强下游通道的二次流。2.1.2󰀁结构计算结果分析

图5和6分别给出了叶片上最大周向位移节点的位移时间曲线和最大应力节点的应力时间曲线。初始时刻离心载荷与定常初载使叶片出现大幅值的振动,叶片的周向位移和应力很快都达到了峰值,随着时间的推移振动幅值越来越小,最大应力值也趋于稳定。转子单通道流场的非定常计算没有考虑上游尾迹的非定常效应,随着流场的进一步收敛稳定,叶片最终处于一个稳定的流场中,受到的载荷也趋于稳定,所以结构位移和应力都趋于定值。

根尾缘处。在单独的定常气动载荷作用下,叶片的最大应力只有17󰀁319MPa,比离心载荷产生的应力要小得多;将耦合计算得到的气动载荷提取出来,再利用单向耦合的方法对叶片进行加载,比较发现,耦合后的气动载荷引起的最大应力为16󰀁912MPa,小于定常气动载荷引起的最大应力,说明耦合后的气动载荷减小了;从量级上来看,双向耦合和单向耦合(离心+定常气动载荷)的位移和应力都相当于单独的离心载荷与单独的气动载荷作用的叠加,由于气动载荷相对离心载荷较小,而且耦合前后变化小,所以双向耦合与单向耦合的位移、应力值差别很小。

表2󰀁各种载荷下计算结果的最大位移和应力对比

最大Von

Mises应力

离心载荷

单向耦合(定常气动载荷)单向耦合(耦合气动载荷)

单向耦合(离心+定常气动载荷)双向耦合(离心+耦合气动载荷)

430󰀁29MPa17󰀁319MPa16󰀁912MPa441󰀁47MPa441󰀁31MPa

最大位移0󰀁0052522m

0󰀁0006382m0󰀁0006173m0󰀁0058788m0󰀁0058285m

2󰀁2󰀁非设计状态下耦合计算结果分析

非设计状态的计算与设计状态采用相同的耦合模型和计算方法,只是流场的进出口边界条件不一样,具体是根据非设计状态下整级的定常结

图5󰀁最大周向位移节点的位移时间曲线

果来确定的。2.2.1󰀁流场计算结果分析

通过对非设计状态下耦合前后流场计算结果的分析,可以得出与设计状态下基本一致的结论,但相比之下,非设计状态也有区别于设计状态的一些特点。图7给出了流场6个监测点的静压变化曲线,可以看出,从计算初期到收敛稳定的静压波动量要比设计状态大一些。

图6󰀁最大应力节点的应力时间曲线

󰀁󰀁为进一步分析MFX双向耦合计算结果的叶片结构应力和位移情况,又通过单向耦合方式对叶片进行了静强度计算,以此作为对比。单向耦合采用与前面MFX完全相同的模型,流场和叶片结构网格不变,单向耦合的气动载荷由前面的定常计算以及耦合后的流场来提供,叶片结构采用ANSYS软件中的静力分析方法(Static)。另外,还对叶片单独进行了离心载荷的静应力计算。

表2给出了各种载荷下计算结果的最大位移和应力值,叶片的最大应力为441󰀁31MPa,位于叶图7󰀁叶片表面近壁区监测点的静压变化曲线

表3对耦合收敛流场与定常初场的计算结果进行了对比,从中可以看出,在保证膨胀比基本一致的情况下,流量、效率和出口相对马赫数在耦合

前后都有很大的差异,通道面积和设计状态一样仍第3期󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁

热力透平

157

然是减小的,但是流量却增大了。具体分析可以发现,在非设计状态下通道中有很强的旋涡流动,如图8所示,它占据了通道很大一部分空间,耦合计算考虑了结构变形对流场的影响;从图中可以看出叶片的结构变形使分离涡的涡核减小了,因此主流的流通面积相对扩张了,所以流量增大了,主流的

总温比1󰀁0132741󰀁0132830󰀁0000090󰀁0009%

膨胀比

1󰀁2006741󰀁2007870󰀁0001130󰀁0094%

流量(kgs-1)

88󰀁45368-190󰀁034561󰀁580881󰀁7872%

相对扩张也导致了出口相对马赫数的减小。另外,这里所考虑的效率是用静焓计算得到的内效率,流场静压、速度变化对静焓影响很大,效率在耦合前后的较大差异正是由于流场静压和速度的波动引起的。从以上分析可以说明,在非设计状态下,流场对结构变形具有更强的敏感性。

效率0󰀁5213250󰀁414382-0󰀁106943-20󰀁5137%

出口相对Ma

0󰀁5715030󰀁5378-0󰀁026125-4󰀁5713%

喉道附近通道面积(m2)

0󰀁04367870󰀁0432167-0󰀁0004620-1󰀁0577%

表3󰀁耦合收敛流场与定常初场计算结果对比

定常初场

耦合收敛流场绝对误差相对误差

3󰀁结论

利用商业软件CFX和ANSYS对转子域的单

定常计算󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁耦合计算

图8󰀁转子通道叶中处的分离涡

个通道进行了设计状态和非设计状态下的流固耦合计算,耦合从双向和单向两个角度进行,并对结果进行了对比分析。

(1)流场计算结果表明:设计状态下的结构变形对流场的影响不是很大,耦合前后流场总体参数存在一定量的差别,由于计算选取的叶片展弦比比较大,结构变形对流场的边界(特别是在叶尖位置)存在显著的影响,迫使尾迹的方向发生改变;非设计状态下,耦合前后流场总体参数存在较大差异,通道中的分离涡受到叶片结构变形的影响,整体上流场对结构变形具有更强的敏感性。

(2)结构计算结果表明:两种计算状态下,气动载荷相对离心载荷较小,而且耦合前后变化小,所以双向耦合与单向耦合的位移、应力值差别很小;在没有非定常来源的稳态气动载荷作用下,非设计状态下气动载荷小,由气动载荷引起的位移、应力相应地小于设计状态。

参考文献:

[1]刑景堂,周盛,崔尔杰.流固耦合力学概述[J].力学进展,

1997,27,(1):19~30.

[2]孟庆国,周盛.叶轮机械非定常流动研究进展3[J].力学进

展,1997,27,(2):232~244.

[3]MarshallJG,ImregunM.AReviewofAeroelasticityMethodswith

EmphasisonTurbomachineryApplications[J].JournalofFluidsandStructures,1996,Vo.l10,pp237~267.[4]刘宝杰,邹正平,严明,等.叶轮机计算流

2.2.2󰀁结构计算结果分析

非设计状态下结构分析仍然采用两种耦合方式。表4给出了各种载荷下计算结果的最大位移

和应力值,和设计状态一样,非设计状态下双向耦合和单向耦合(离心+定常气动载荷)的位移、应力值差别也很小;在非设计状态下利用单向耦合计算出气动载荷引起的最大应力只有4󰀁12MPa,这主要是因为选取的非设计点膨胀比只有1󰀁2,气体对叶片加功量小,所以气动载荷小,图9具体比较了两种状态下叶片受到的气动载荷。

表4󰀁各种载荷下计算结果的最大位移和应力对比

最大Von

Mises应力430󰀁29MPa4󰀁1211MPa3󰀁4298MPa434󰀁40MPa433󰀁02MPa

最大位移0󰀁0052522m

0󰀁0001434m0󰀁0001235m0󰀁0055348m0󰀁0053581m

离心载荷

单向耦合(定常气动载荷)单向耦合(耦合气动载荷)

单向耦合(离心+定常气动载荷)双向耦合(离心+耦合气动载荷)

设计状态󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁非设计状态图9󰀁叶片法向力周向分量的收敛过程

(下转第163页)

第3期󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁󰀁

热力透平

163

基于上述计算结果可以判断:转速低于2220r/min时,910叶片不会出现严格意义上的整圈节径振动;转速处于2220r/min~20r/min之间时,910叶片从󰀂围带自锁的整圈振动 向󰀂围带和拉筋均自锁的整圈振动 过渡;转速高于20r/min时,ILB叶片完全锁紧,达到󰀂围带和拉筋整圈自锁 。

需要补充说明的是,由于实际运转过程中拉筋的接触对围带的接触具有干涉作用,故而状态C、E和F仅能作为计算的参考。如果将从开始耦合直到完全锁紧的󰀂过渡区 扩大为:1105~20r/min,可能更加稳妥。在这个过渡区里,ILB叶片的整圈振动特性既不同于完全锁紧的整圈振动,也不同于自由叶片的振动。从计算角度,这个区域里的振动特性是很难精确预测的。如前文所述,在对910mmILB叶片进行动频

[6]

测试时,在转速2400r/min附近发现了不符合整圈节径振动特性的共振点,这一现象似乎可以应用本文提出的观点进行解释,因为2400r/min刚好落在2220~20r/min的过渡区范围之内。

综上所述,由于对ILB叶片的所有振动特性的计算都建立在整圈自锁的假定之上,故而应当特别留意该假定在实际工作状态中(尤其是在额定转速下)是否能够保证。在叶片的基本结构和尺寸已经确定之后,还需要对围带(凸台)的间隙

进行调整,并控制其公差,以此达到叶片锁紧力的要求。最后,在样机试验过程中,可以由频率测量的结果对叶片自锁情况作进一步的确认。

4󰀁结论

ILB叶片自身带有围带和拉筋有助于叶片形成整圈自锁的结构,提高刚度,降低应力。然而在设计过程中应当注意围带和拉筋的锁紧条件是否满足,如有必要则应对其进行合理调整,这是确保ILB叶片的振动特性符合设计要求的前提。

参考文献:

[1]苏振华,周宏利,周克澄.火力发电设备技术手册第二卷汽轮

机[M].北京:机械工业出版社,1998.

[2]谢永慧,孟庆集.汽轮机叶片三维有限元模型的建立及静态与

动态应力分析[J].热力透平,2003(1):36~42.

[3]秦飞,李英,张晓峰.基于有限元方法的汽轮机叶片动强度评

估[J].北京工业大学学报,2006,32(8):765~768.[4]LFWagner.ComputationoftheVibratoryCharacteristicsof

GroupedandInterlockedBlades[A].Orlando:TechnicalReportofWestinghouse,1996.

[5]褚玉喜,陆伟.整圈自锁(ILB)的设计及其动静频试验[J].

上海汽轮机,2002(3):15~20.

[6]肖俊峰,等.910mm叶片动频率试验初步分析结果[Z].上

海汽轮机有限公司内部资料,2007.

(上接第157页)

󰀁󰀁体动力学技术现状与发展趋势[J].航空学󰀁󰀁报,2002,23,(5):394~403.

[5]SaymaAI,VahdatMi&ImregunM.AnIntegratedNonlinearAp󰀂

proachforTurbomachineryForcedResponsePrediction.PartI:Formulation[J].JournalofFluidsandStructures,2000,Vo.l14,pp87~101.

[6]孟越,李琳,李其汉.尾流激振情况下叶片强迫响应瞬态分析

方法[J].北京航空航天大学学报,2006,32,(6):671~674.[7]赵书军.压气机三维非定常流计算与叶片强迫响应研究[D].

西安:西北工业大学,2003.

[8]曾强.压气机转子叶片流固耦合计算及软件集

成研究[D].南京:南京航空航天大学,2006.

[9]杨青真,施永强,肖军,周新海.气固耦合振动叶栅非定常流动

分析研究[J].,2005,23,(2):167~171.

[10]袁新,金琰,畅国勇.叶轮机械中流体激振问题的流固耦合研

究[A].第八届全国空气弹性学术交流会会议论文集[C],2003.

[11]JamesonA.TimeDependedCalculationsUsingMulti󰀂gridwith

ApplicationstoUnsteadyFlowsPastAirfoilsandWings[R].AIAA91󰀂1596,1991.

[12][德]W.瓦格纳,A.克鲁泽著,项红卫译.水和蒸汽的性质

[M].北京:科学出版社,2003,pp58~59.