火箭發(fā)動(dòng)機(jī)渦輪葉片動(dòng)態(tài)響應(yīng)計(jì)算及方法

任 眾，許開富，朱東華

(西安航天動(dòng)力研究所，陜西 西安 710100)

0 引言

渦輪葉片作為火箭發(fā)動(dòng)機(jī)的關(guān)鍵部件，要求其在高溫、高負(fù)荷、高轉(zhuǎn)速的工作環(huán)境下安全可靠地工作。局部進(jìn)氣的沖擊式渦輪為液體火箭發(fā)動(dòng)機(jī)中常見的一種渦輪形式，工作過程中葉片受到周期性高速氣流激振力的沖擊作用，其振動(dòng)及疲勞問題十分突出和復(fù)雜[1-2]。

現(xiàn)有的液體火箭發(fā)動(dòng)機(jī)渦輪由于工作時(shí)間短，設(shè)計(jì)過程中主要考慮氣動(dòng)性能及靜強(qiáng)度[3-5]，很少考慮由氣體激振引起的葉片振動(dòng)疲勞問題。但隨著可重復(fù)使用發(fā)動(dòng)機(jī)的研制，要求各零部件的服役時(shí)間成倍增加，因此，必須對(duì)渦輪葉片的高周疲勞問題開展深入的研究。而該問題的核心就是分析葉片在氣流激振力作用下的動(dòng)態(tài)響應(yīng)，獲得其振動(dòng)應(yīng)力[6-7]。

目前，葉片的動(dòng)態(tài)響應(yīng)計(jì)算分析主要有兩種方法。一種是理論模型的方法，將葉片簡(jiǎn)化成懸臂梁，根據(jù)不同的氣流激振力的形式獲得激振力的簡(jiǎn)化模型，再建立梁的彎曲振動(dòng)微分方程，求解方程[8-11]；另一種是三維仿真的方法，首先通過三維非定常流場(chǎng)仿真獲得激振力的時(shí)域分布，再進(jìn)行葉片在激振力作用下的振動(dòng)響應(yīng)仿真計(jì)算[12-13]。其中，第一種方法已在風(fēng)機(jī)葉片及航空發(fā)動(dòng)機(jī)壓氣機(jī)葉片中得到了成熟的應(yīng)用[14-16]，其激振力形式為連續(xù)周期性變化的，但局部進(jìn)氣的沖擊式渦輪的氣流激振力形式有著明顯的不同，其激振力形式為周期性脈沖變化，同時(shí)氣流的沖擊馬赫數(shù)也要高得多，由該激振力產(chǎn)生的動(dòng)態(tài)響應(yīng)規(guī)律會(huì)有所不同。目前，對(duì)于該種形式的激振力作用下的渦輪葉片動(dòng)態(tài)響應(yīng)的研究文獻(xiàn)較少；而第二種方法理論上可以得到準(zhǔn)確的葉片振動(dòng)響應(yīng)情況，但涉及到多學(xué)科的仿真分析，計(jì)算工作量大，占用資源多，計(jì)算時(shí)間長(zhǎng)，不利于工程設(shè)計(jì)階段對(duì)疲勞問題的快速評(píng)估。

本文提出了一套針對(duì)局部進(jìn)氣沖擊式渦輪葉片動(dòng)態(tài)響應(yīng)的理論計(jì)算模型，并以某液體火箭發(fā)動(dòng)機(jī)渦輪為研究對(duì)象，采用該方法進(jìn)行葉片在流體激振力下的動(dòng)態(tài)響應(yīng)計(jì)算，并與三維仿真結(jié)果對(duì)比，驗(yàn)證結(jié)果的準(zhǔn)確性，為發(fā)動(dòng)機(jī)渦輪高周疲勞設(shè)計(jì)提供依據(jù)。

1 計(jì)算方法

1.1 理論模型

1.1.1 受迫振動(dòng)微分方程

本文的理論模型以懸臂梁受迫振動(dòng)理論為基礎(chǔ)[17-19]。由于渦輪葉片主要是受氣動(dòng)彎矩的作用產(chǎn)生彎曲變形，同時(shí)液體火箭發(fā)動(dòng)機(jī)渦輪葉片大多數(shù)為較短的直葉片，剪切變形以及截面繞中性軸轉(zhuǎn)動(dòng)慣量對(duì)葉片振動(dòng)的影響較大。因此，將葉片簡(jiǎn)化成等截面的鐵木辛柯懸臂梁，并考慮離心力及阻尼的影響，推導(dǎo)得到梁的受迫振動(dòng)微分方程為

(1)

初始邊界條件

位移邊界條件

y(0,t)=0,y′(0,t)=0,y″(l,t)=0,y?(l,t)=0

其中

f(x,t)=f1(x)·f2(t)

式中：y為梁的位移；x為梁的軸向坐標(biāo)；t為時(shí)間；E為彈性模量；G為剪切模量；ρ為密度；I為截面主慣性矩；J為單位長(zhǎng)度梁對(duì)截面中性軸的轉(zhuǎn)動(dòng)慣量；k為截面形狀因子，對(duì)于矩形截面，k=5/6；A為截面面積；l為梁的長(zhǎng)度；c為黏性外阻尼系數(shù)；Fa為葉片離心力大小，由于實(shí)際工作過程中葉片的離心力沿葉高是變化的，為了使方程有理論解，假設(shè)葉片各截面的離心力均為平均離心力；ω為轉(zhuǎn)速；At為葉片截面面積；Rt為葉頂半徑；Rh為葉根半徑；f(x,t)為葉片單位高度所受的激勵(lì)力大??；f1(x)為氣流激振力沿葉高的分布函數(shù);f2(t)為單位氣流激振力隨時(shí)間的分布函數(shù)。

函數(shù)f1(x)通過圖1中的葉片通道的控制體動(dòng)量守恒可得[20]

ρg2(x)Ca2(x)Cu2(x)]

(2)

(3)

p=ρgRTg

(4)

式中：下標(biāo)1為葉片入口截面參數(shù)；下標(biāo)2為葉片出口截面參數(shù)；Zb為葉片數(shù)；ρg(x)為氣流密度的徑向分布；Ca為氣流軸向速度的徑向分布；Cu為氣流軸向速度的徑向分布；p為氣體壓力；Tg為氣體溫度；R為氣體常數(shù)。

圖1 葉片通道控制體示意圖Fig.1 Schematic diagram of blade channel control body

函數(shù)f2(t)的分布如圖2所示，其中局部進(jìn)氣1個(gè)周期內(nèi)流體激振力的作用時(shí)間

Δt=εT

(5)

式中：ε為渦輪局部進(jìn)氣度；T為1個(gè)旋轉(zhuǎn)周期的時(shí)間。

圖2 激振力的時(shí)間分布函數(shù)Fig.2 Time distribution function of excitation force

1.1.2 方程求解方法

采用模態(tài)疊加法求解式(1)，求解的形式可表示為[17]

(6)

式中：Yi(x)為梁的第i個(gè)主振型函數(shù)；qi(t)為第i個(gè)廣義坐標(biāo)權(quán)函數(shù)。在實(shí)際求解時(shí)，取前3到5階振型疊加就可得到足夠精確的結(jié)果。其中，各階主振型函數(shù)即為梁的自由振動(dòng)的解，令式(1)右側(cè)的激振力為零，再結(jié)合懸臂梁的位移邊界條件(3)，推導(dǎo)得到梁彎曲振動(dòng)的特征

(7)

式中λ1、λ2分別為特征方程A1X4+B1X2+C1=0的兩個(gè)非負(fù)解。其中

C1=-ω_i2ρA

求解上述方程的第i個(gè)解ω_i，即為梁的第i階固有頻率，對(duì)應(yīng)的振型函數(shù)

Yi(x)=[cosh(λ1x)-cos(λ2x)]+

(8)

為了求解各廣義坐標(biāo)下的權(quán)函數(shù)，將式(6)代入式(1)，并利用振型函數(shù)的正交性，推導(dǎo)得

(9)

其中

式(9)左側(cè)為標(biāo)準(zhǔn)的阻尼系統(tǒng)的受迫振動(dòng)方程，而右側(cè)可看做振幅為D2的周期激振力。將該激勵(lì)力看做一系列的微脈沖組合，利用杜哈梅積分，得到梁對(duì)該激振力的響應(yīng)

qi(t)=eat[X1cos(bt)+X2sin(bt)]+

(10)

其中

系數(shù)X1、X2根據(jù)初始條件(2)確定。獲得了梁的位移響應(yīng)分布y(x,t)，可利用材料力學(xué)中的關(guān)系獲得某點(diǎn)A的應(yīng)力

(11)

式中yA為A點(diǎn)距截面彎曲主軸的距離。

1.2 三維仿真的方法

三維仿真的計(jì)算流程見圖3，由圖可知，整個(gè)計(jì)算流程包括CFD非定常計(jì)算、載荷傳遞、動(dòng)態(tài)響應(yīng)計(jì)算分析3部分。首先對(duì)渦輪整周模型進(jìn)行三維瞬態(tài)流場(chǎng)仿真，獲得某個(gè)葉片表面的壓力分布隨時(shí)間的變化，即葉片的氣體激振力；然后通過插值運(yùn)算，將流固交界面上的流場(chǎng)壓力映射到葉片表面的有限元網(wǎng)格節(jié)點(diǎn)上；再進(jìn)行葉片的動(dòng)態(tài)響應(yīng)計(jì)算計(jì)算，獲得葉片在激振力作用下的響應(yīng)情況；在此基礎(chǔ)上，結(jié)合材料的疲勞性能數(shù)據(jù)，可以做葉片的高周疲勞分析。整個(gè)計(jì)算過程在ANSYS Workbench平臺(tái)下完成。

圖3 三維仿真計(jì)算方法Fig.3 Calculation method of 3D simulation

2 算例分析

2.1 研究對(duì)象

本文以某發(fā)動(dòng)機(jī)渦輪為研究對(duì)象，計(jì)算葉片在氣體激振力下的動(dòng)態(tài)響應(yīng)。該渦輪為單級(jí)局部進(jìn)氣沖擊式渦輪，工作時(shí)高溫高壓燃?xì)馔ㄟ^單個(gè)噴嘴加速，噴出的高速氣流沖擊渦輪葉片做功。其中，葉片的三維模型見圖4(a)，由圖可知，葉片為直葉片，主要幾何參數(shù)如表1所示。

表1 渦輪葉片主要幾何參數(shù)Tab.1 Main geometric parameters of turbine blades

渦輪葉片材料為高溫合金，密度為8 240 kg/m3，工作溫度下的彈性模量為146 GPa，泊松比為0.33，屈服極限為860 MPa。

2.2 理論模型求解

根據(jù)截面慣性矩相等，將渦輪葉片簡(jiǎn)化為矩形截面的懸臂梁，如圖4所示。簡(jiǎn)化梁的長(zhǎng)L、寬W、高H分別為29.5、9.45、5.44 mm。此外，由于液體火箭發(fā)動(dòng)機(jī)渦輪葉片短，剛性大，一般將葉片和盤考慮為一個(gè)葉盤系統(tǒng)，因此本文單獨(dú)將葉片簡(jiǎn)化為梁結(jié)構(gòu)，未考慮輪盤的阻抗作用，計(jì)算結(jié)果與實(shí)際存在一定的偏差。

圖4 葉片模型簡(jiǎn)化Fig.4 Blade model simplification

首先對(duì)梁進(jìn)行模態(tài)計(jì)算，計(jì)算得到的前三階固有頻率如表2所示，相應(yīng)的振型如圖5所示。由表2可知，理論模型計(jì)算的固有頻率比三維仿真值偏小，一方面是由于理論模型中假設(shè)葉片各截面的離心力均為平均離心力，從而導(dǎo)致旋轉(zhuǎn)剛化作用對(duì)固有頻率影響的差異；另一方面是由于理論模型將葉片截面簡(jiǎn)化為矩形截面，截面形狀的不同也會(huì)導(dǎo)致固有頻率的差異。

表2 葉片前三階頻率對(duì)比Tab.2 Comparison of the first three frequencies of blades

圖5 梁的前三階振型Fig.5 The first three modes of beam

然后采用模態(tài)疊加法求解梁在氣流激振力作用下的位移響應(yīng)及應(yīng)力響應(yīng)，其中，由理論方法得到的氣體激振力載荷沿葉高的分布規(guī)律如圖6所示，該激振力在1個(gè)旋轉(zhuǎn)周期內(nèi)的作用時(shí)間為0.3 ms；為了對(duì)兩種方法進(jìn)行精確比較，理論模型與三維仿真的阻尼系數(shù)均取0.2。對(duì)前三階的模態(tài)振型疊加，提取了梁的根部A點(diǎn)的應(yīng)力及頂部B點(diǎn)的位移結(jié)果，A點(diǎn)和B點(diǎn)的位置示如圖7所示。梁在3個(gè)周期內(nèi)的頂部位移響應(yīng)如圖8中的虛線所示，根部的應(yīng)力響應(yīng)如圖9中的虛線所示。

圖6 梁的激振力分布圖Fig.6 Distribution of excitation force on beam

由圖8可知，在一次氣流激振力的沖擊作用下，葉片的響應(yīng)頻率約為4 066 Hz，同時(shí)由于阻尼的作用，響應(yīng)幅值逐漸衰減，經(jīng)歷7個(gè)振蕩后振幅基本衰減到零，每個(gè)旋轉(zhuǎn)周期內(nèi)該響應(yīng)重復(fù)1次。以該發(fā)動(dòng)機(jī)為例，在單次工作時(shí)間范圍內(nèi)大約經(jīng)歷40萬次上述循環(huán)過程；由圖9可知，葉片根部的平均應(yīng)力約為295 MPa，在氣流沖擊后的第一個(gè)響應(yīng)動(dòng)應(yīng)力幅值最大，約為155 MPa，然后由于阻尼的作用動(dòng)應(yīng)力幅值逐漸衰減到零。

圖7 位置示意圖Fig.7 Location diagram

圖8 頂部位移響應(yīng)結(jié)果Fig.8 Response results of top displacement

圖9 根部應(yīng)力響應(yīng)結(jié)果Fig.9 Response results of root stress

2.3 三維仿真求解

渦輪的流場(chǎng)三維仿真模型及邊界條件見圖10。由圖可知，計(jì)算域由噴嘴通道與動(dòng)葉通道兩部分組成。整個(gè)計(jì)算域均采用結(jié)構(gòu)網(wǎng)格，單元總數(shù)475萬。噴嘴入口為總壓邊界，動(dòng)葉出口為靜壓邊界，級(jí)間采用瞬態(tài)轉(zhuǎn)靜交接面，其余均為無滑移絕熱壁面邊界。利用商用軟件CFX開展瞬態(tài)仿真計(jì)算，每個(gè)葉片通道設(shè)置5個(gè)時(shí)間步，一共計(jì)算10個(gè)旋轉(zhuǎn)周期，采用SST湍流模型及高階差分格式對(duì)方程進(jìn)行離散求解。

圖10 流場(chǎng)計(jì)算模型及邊界Fig.10 CFD model and boundary

計(jì)算得到的葉片通道內(nèi)的馬赫數(shù)分布如圖11所示，由圖11可知，一周大概8個(gè)葉片受到高速氣流的沖擊作用。單個(gè)渦輪葉片受到的氣動(dòng)彎矩分布如圖12所示，由圖12可知，單個(gè)葉片在一個(gè)旋轉(zhuǎn)周期內(nèi)受到氣流的作用時(shí)間比例約為0.155，與式(3)計(jì)算的作用時(shí)間基本吻合，載荷的形狀接近圖2的方波。

圖11 葉片通道馬赫數(shù)分布Fig.11 Mach number distribution of blade passage

通過上文的非定常流場(chǎng)仿真，獲得了葉片的激振力，在此基礎(chǔ)上開展動(dòng)態(tài)響應(yīng)仿真。有限元計(jì)算模型及邊界如圖13所示。由圖13可知，模型由葉片及底座兩部分組成，均采用六面體結(jié)構(gòu)網(wǎng)格進(jìn)行劃分，網(wǎng)格節(jié)點(diǎn)1.25萬。首先對(duì)葉片進(jìn)行模態(tài)分析，通過對(duì)比各階頻率的大小來驗(yàn)證理論梁模型的準(zhǔn)確性。然后利用商用軟件ANSYS DYNA進(jìn)行葉片的動(dòng)態(tài)響應(yīng)仿真，獲得葉片的位移及動(dòng)應(yīng)力響應(yīng)分布。

圖12 葉片氣動(dòng)彎矩分布Fig.12 Distribution of aerodynamic bending moment of blade

圖13 有限元模型及邊界Fig.13 FEM model and boundary

模態(tài)計(jì)算得到的葉片前三階固有頻率如表2所示，由表2可知，三維仿真計(jì)算得到的固有頻率比理論模型偏大，且隨著階次的增大，偏差逐漸增大，但最大誤差不超過10%，說明采用理論梁模型對(duì)葉片進(jìn)行簡(jiǎn)化的準(zhǔn)確性。

計(jì)算得到的葉片根部應(yīng)力響應(yīng)如圖9中的實(shí)線所示，由圖9可知，三維仿真得到的平均應(yīng)力約為305 MPa，響應(yīng)頻率為4 000 Hz，響應(yīng)動(dòng)應(yīng)力的最大幅值為175 MPa。與理論模型方法相比，動(dòng)應(yīng)力平均值、頻率及最大幅值誤差分別為3.39%、1.62%、11.4%，從而驗(yàn)證了理論模型進(jìn)行葉片動(dòng)態(tài)響應(yīng)計(jì)算的準(zhǔn)確性。

3 結(jié)論

本文建立了一套針對(duì)局部進(jìn)氣沖擊式渦輪葉片的動(dòng)態(tài)響應(yīng)的理論計(jì)算模型，并以某液體火箭發(fā)動(dòng)機(jī)渦輪為研究對(duì)象，分別采用理論模型的方法與三維仿真的方法計(jì)算了葉片在氣流及激振力作用下的動(dòng)態(tài)響應(yīng)。主要結(jié)論如下：

1)根據(jù)液體火箭發(fā)動(dòng)機(jī)渦輪葉片的特點(diǎn)，將葉片簡(jiǎn)化成等截面的鐵木辛柯懸臂梁，同時(shí)考慮離心力及阻尼的影響，推導(dǎo)得到梁的受迫振動(dòng)微分方程，并采用模態(tài)疊加法對(duì)方程進(jìn)行求解。

2)采用流固耦合的方法對(duì)葉片動(dòng)態(tài)響應(yīng)進(jìn)行三維仿真，首先通過非定常全流場(chǎng)仿真獲得氣流激振力，在此基礎(chǔ)上開展葉片的三維有限元模態(tài)分析及動(dòng)態(tài)響應(yīng)分析。

3)與三維仿真計(jì)算方法對(duì)比，采用理論模型方法計(jì)算得到的葉片固有頻率誤差不超過10%，同時(shí)根部平均動(dòng)應(yīng)力低3.39%，響應(yīng)頻率高1.62%，響應(yīng)最大幅值低11.4%。驗(yàn)證了該理論模型的準(zhǔn)確性，可以采用該方法開展發(fā)動(dòng)機(jī)渦輪葉片高周疲勞設(shè)計(jì)。