核电超长末叶气动稳定性的流固耦合数值研究

庾明达, 徐自力, 范志飞, 曹守洪, 范小平, 方宇

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

汽轮机中的颤振是一种持续、不衰减的自激振动现象,是由汽轮机流道内的气流力与叶片发生相互作用引起,所以研究叶片的气动稳定性需考虑流固耦合的作用。目前,研究颤振稳定性应用最多的流固耦合方法为能量法,仲继泽等通过能量法预测了相位角、流量等参数对压气机叶片的颤振稳定性的影响[1-2]。能量法思路清晰、简便,但由于假设过多,与实际情况不很相符[3]。另有很多学者采用流固耦合方法来模拟叶片的实际振动响应特征,Vahdati等用强耦合CFD方法研究了宽弦风扇的失速颤振现象,但此方法求解难度大,占用资源高[4];Dowell等基于带副翼操纵面的机翼3DOF模型,运用ROM降阶技术研究了临界音速条件下马赫数对机翼颤振速度和稳定性的影响[5],这种气动力降阶模型目前多局限于机翼等具有简单几何外形结构的流固耦合计算;文献[6-7]利用模态叠加法将叶片响应形式分解为单阶模态的线性振动,研究了跨音速压气机的颤振边界。还有不少学者使用二维模型模拟压气机的流场特征[8-10],但二维模型难以完全模拟高负载叶轮机械的流动状况。

目前,许多研究文献虽利用上述简化耦合模型研究了流场流动和结构振动特征,然而大多数文献集中于研究机翼或压气机叶片的颤振稳定性[11-14],对核电汽轮机末级叶片,尤其是本文所研究的超长核电末级叶片的气动特征鲜有涉及。

因此,本文采用三维流固耦合时域响应算法,研究了新开发的核电汽轮机超长末级叶片在设计背压3 600 Pa、高背压7 200 Pa和10 800 Pa工况下的气动稳定性问题。计算结果表明,在所计算的背压工况内叶片维持气动稳定,未发现颤振现象,但稳定性逐渐降低。

1 计算原理和方法

1.1 流体控制方程

为完全模拟核电叶片不同叶高的流动特征,采用三维Navier-Stokes(N-S)方程作为控制方程,该方程组表示流体运动所应满足的质量守恒,动量守恒和能量守恒三大物理学定律,可写为

(1)

式中:Q为瞬态矢量;F、G、H分别为3个方向上的通量;S为源性通量。

虽然瞬时N-S方程可以用来描述湍流,但其非线性特征不易求解,工程上常引入雷诺应力的方法将N-S方程中瞬态的变化量以时均化的形式表现,引入湍流控制方程把湍流的脉动值与时均值联系起来,方程变为雷诺平均纳维斯托克斯方程(RANS)。本文采用的k-ε湍流控制方程为

(2)

(3)

1.2 固体求解算法

固体域的动力学方程可写为

(4)

式中:M、C、K分别为质量、阻尼和刚度矩阵;q(t)为位移;p(t)为作用在流固耦合交界面上的气动压力。采用直接积分法求解叶片的振动方程得到叶片的有限元节点振动位移。

1.3 流固耦合的方法及计算流程

流固耦合计算的实质是流体和固体的交界面上发生相互作用,进而产生耦合并相互影响。由于流场网格与固体网格并不一一对应,所以在流固耦合交界面上,需运用样条方法,将流固交界面节点位移做转换

Δxa=GΔxs

(5)

式中:xa为流场位移;xs为固体位移;G为样条矩阵。反之,流场气动力需要施加到固体交界面上,根据做功等效原理有

qs=GTqa

(6)

式中:qa、qs为施加在流固交界面上的气动力。在流固耦合计算过程中,流场完成内部迭代计算后,按照式(6),与固体交界面交换压力等变量,之后固体产生位移,此时流场网格变形,以便进行下一次的迭代计算。在流固耦合交界面处,流场网格节点直接进行变形,其他较远处的网格变形遵循以下位移扩散方程

·(Γdispδ)=0

(7)

(8)

图1 时域响应流固耦合算法的流程图

2 末级汽轮机有限元网格模型

本文研究的对象为新开发的核电汽轮机,末级叶片长度超过2 m,静叶数为66,动叶数为76,蒸汽从末级叶片的静叶入口流入,随后经动叶出口流出,几何示意图如图2所示。

图2 末级叶片通道几何示意图

本文采用三维计算模型和循环对称边界的单通道流道作为计算对象,用以计算末级流场,这样可有效减少计算规模,提高计算效率。为了选取合适的网格,本文采用由疏到密具有不同单元数量、不同网格尺寸的4套流场网格,对设计背压工况的定常流动进行了计算,得到了流体施加在动叶片中部截面背弧面的压力分布,如图3所示。图3中,X为叶片中部展弦方向任意一点的位置坐标,C为叶片中部的弦长。

图3 设计背压工况下的动叶中部背弧面流体压力分布

由图3可知,随着网格单元数的不断增多,计算得到的背弧面静压力分布趋于定值,在流场网格单元数为707 100时,继续增加网格单元数,背弧面的静压力计算结果也基本保持不变。因此,对该核电末级叶片,单元数选取707 100。

图4给出了核电末级叶片的网格,其中叶片固体网格含有13 041个单元。动静叶流体计算域的网格均为结构化网格,对叶片周围的流动区域设置了边界层,内部计算域网格环绕边界层分布直至充满整个计算域。动、静叶网格的平均第一层边界层网格高度约为3.1×10-5m,对第一层网格高度做无量纲化处理,可得流场网格的y+≈4<5。

(a)末级流道整体 (b)叶片固体

(c)静叶叶顶流道 (d)动叶叶顶流道图4 核电末级叶片的计算网格

3 计算结果与讨论

3.1 流场特征分析

采用单通道模型,在流场计算域中,动静叶流道两侧均施加循环对称边界,将静叶与动叶流场域统一视为内部区域,以此模拟动叶与静叶之间的相对旋转运动,叶片表面、叶顶、叶毂为无滑移壁面,静叶进口设置蒸汽的流量和温度,在动叶出口设置压力值。对设计背压3 600 Pa和7 200、10 800 Pa两个高背压工况下的定常流动进行了计算,得到了3 600、7 200 Pa工况下末级流道截面流线分布,如图5所示,可知在设计背压时,流体从动叶流出后,流动保持平稳,而在提高背压后,流动未发生分离或是方向改变。而当背压提高至10 800 Pa时,动叶尾部流道内发生了旋涡现象,此时流线分布及压力分布如图6所示。由图6可知,流体是顺着压力梯度向下游流动,其流动趋势与压力梯度保持一致,然而在该背压工况下动叶流道中部的压力梯度发生了明显的波动,出现几个局部低压凸口,流体顺着压力梯度梯度流经这些区域时,由于动能已经减小,不足以克服压差导致流向改变,少量流体在动叶流道近出口处形成了旋涡。

(a)背压为3 600 Pa (b)背压为7 200 Pa图5 子午面流线分布

(a)流线分布 (b)压力分布图6 背压为10 800 Pa时子午面流线分布及压力分布

3.2 叶片气动稳定性分析

(a)3 600 Pa响应曲线 (b)3 600 Pa频谱图

(c)7 200 Pa响应曲线 (d)7 200 Pa频谱图

(e)10 800 Pa响应曲线 (f)10 800 Pa频谱图图7 各背压工况下的尾缘响应曲线及其频谱图

在不同背压工况下动叶顶尾缘点上的时间历程响应曲线及其频谱图如图7所示。由图7可知,3个背压工况下的响应曲线全部现为衰减趋势,响应曲线中的主频成分与动叶的第1阶固有动频ω1的比值接近1,即叶片是以动叶的第1阶固有动频为主进行自由衰减响应,对叶片进行的模态分析结果显示,叶片第2,3,4阶模态与第1阶模态的固有频率之比分别为1.61,2.57,3.39。由图7d可知,频谱图中含有叶片2阶固有频率附近的振动成分,表明在背压为7 200 Pa时,叶片表现为1阶、2阶模态的叠加振动。随着背压的逐步增大,叶片的幅值均分别呈现出逐渐减小的趋势,即叶片在这3个背压工况下都将保持气动稳定,不会发生颤振。

表1、表2给出了在3种工况下,叶片振动至最大位移和最小位移时流场的轴向力和周向力。由表1、表2可知,在3个背压条件下,当叶片从最大位移处振动至最小位移处时,流场的轴向力和周向力均发生了变化,即固体的振动影响了流体气动力的分布。

表1 不同背压下叶片振动至最大、最小位移时的流场轴向力

表2 不同背压下叶片振动至最大、最小位移时的流场周向力

对数阻尼比的计算公式为

(9)

式中:m为周期间隔;A为振动幅值。计算时,由于该核电叶片的材料阻尼很小,式(7)的阻尼比主要为气动阻尼比。气动阻尼比与背压的关系如图8所示,可知随着背压的提高,3个工况下气动阻尼均为正,但提高背压后阻尼比随之减小,叶片的气动稳定性随之减弱。

图8 气动阻尼比与背压的关系

背压工况下振动历程的相平面图如图9所示。由图9可知:在背压为3 600、7 200 Pa时叶片振动的速度和位移均呈现出杂乱的变化趋势,相平面图呈现出多个频率和幅值的叠加运动,气动阻尼对结构有显著影响,相平面呈现不规则的收敛趋势,体现了强烈的非线性特征;背压升至10 800 Pa时,频率成分变得单一,此时的阻尼比相当小,有限时间段内叶片在相平面图上收敛趋势变得更加缓慢,随着时间的进一步推进,叶片最终收敛到静平衡位置。

(a)背压为3 600 Pa (b)背压为7 200 Pa

(c)背压为10 800 Pa图9 各背压下振动响应的相平面图

80%叶高截面流场马赫数云图和截面展弦向叶片压力分布如图10、11所示。由图10、11可知,叶栅流道内出现了激波,激波不断朝叶片前缘移动,在背压为3 600、7 200 Pa时,流场的背弧面均可观察到稳定的激波;到背压为10 800 Pa时,叶片背弧面上出现了低能流团,同时边界层出现了分离,并促使激波结构发生了变化,此时的激波分为两支,一支为前缘正激波,一支与流道平行使整个激波呈“人”字形分布。

(a)背压为3 600 Pa (b)背压为7 200 Pa

(c)背压为10 800 Pa图10 80%叶高截面流场马赫数云图

图11 80%叶高截面展弦向叶片压力分布

文献[15-16]指出,流道内激波的位置变化是导致叶片气动弹性稳定性差异的主要诱因,当激波从叶尖前缘脱体后,叶片从气流中吸收能量导致振动发散并迅速丧失其气动稳定性。本文研究中,背压为10 800 Pa工况时动叶流道内虽已出现靠近前缘的分岔激波,但并没有脱离动叶前缘叶栅流道,背压的提高使动叶栅内的激波向前移动,但激波没有发生脱体,这在一定程度上解释了3个工况下核电叶片保持气动稳定的原因;在叶片振动相对剧烈的叶高中上部流道内,流动保持稳定,没有紊流产生。当背压为10 800 Pa及叶片开始振动后的最大、最小位移时,提取叶片0.5、0.7叶高截面和靠近叶尖振动观察点的0.95叶高截面的流动分布情况,如图12所示。由图12可知,叶片最大、最小位移时叶片3个截面的动叶叶栅的气流分别沿着叶片内、背弧面向下游流动,未发生附面层分离,叶片保持气动稳定,同时流场在各个截面处的流动平稳。

(a)最大位移时 (b)最小位移时图12 背压为10 800 Pa时叶片不同振动时刻的瞬态流线分布

4 结 论

本文采用三维时域响应流固耦合算法,研究了新开发核电汽轮机末级叶片从设计背压为3 600 Pa逐渐提高到高背压为7 200、10 800 Pa条件下的气动稳定性问题。计算结果表明,随着背压的升高,在叶栅流道内的近出口处引起了流体压力的扰动,在压差的影响下逐渐形成了旋涡现象。在从设计背压提高至两个高背压工况下,核电汽轮机叶片的时域响应均以1阶动频为主分别进行自由衰减,表明该核电超长叶片在所计算的背压范围内将保持气动稳定。背压提高后,动叶栅和流道内的激波随着背压的增大不断向前缘移动,但并未发生脱体,叶片振动至不同时刻时在叶高各截面流道内流动没有发生分离,这解释了叶片保持气动稳定的原因。但是,气动阻尼比随背压的提高逐渐消减,气动稳定性随背压的提高而逐渐减弱。