任 眾，許開富，朱東華

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

0 引言

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

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

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

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

1 計算方法

1.1 理論模型

1.1.1 受迫振動微分方程

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

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

函數(shù)f1(x)通過圖1中的葉片通道的控制體動量守恒可得[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個周期內(nèi)流體激振力的作用時間

Δt=εT

(5)

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

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

1.1.2 方程求解方法

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

(6)

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

(7)

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

C1=-ω_i2ρA

求解上述方程的第i個解ω_i，即為梁的第i階固有頻率，對應(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)的受迫振動方程，而右側(cè)可看做振幅為D2的周期激振力。將該激勵力看做一系列的微脈沖組合，利用杜哈梅積分，得到梁對該激振力的響應(yīng)

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

(10)

其中

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

(11)

式中yA為A點距截面彎曲主軸的距離。

1.2 三維仿真的方法

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

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

2 算例分析

2.1 研究對象

本文以某發(fā)動機(jī)渦輪為研究對象，計算葉片在氣體激振力下的動態(tài)響應(yīng)。該渦輪為單級局部進(jìn)氣沖擊式渦輪，工作時高溫高壓燃?xì)馔ㄟ^單個噴嘴加速，噴出的高速氣流沖擊渦輪葉片做功。其中，葉片的三維模型見圖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ù)截面慣性矩相等，將渦輪葉片簡化為矩形截面的懸臂梁，如圖4所示。簡化梁的長L、寬W、高H分別為29.5、9.45、5.44 mm。此外，由于液體火箭發(fā)動機(jī)渦輪葉片短，剛性大，一般將葉片和盤考慮為一個葉盤系統(tǒng)，因此本文單獨將葉片簡化為梁結(jié)構(gòu)，未考慮輪盤的阻抗作用，計算結(jié)果與實際存在一定的偏差。

圖4 葉片模型簡化Fig.4 Blade model simplification

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

表2 葉片前三階頻率對比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個旋轉(zhuǎn)周期內(nèi)的作用時間為0.3 ms；為了對兩種方法進(jìn)行精確比較，理論模型與三維仿真的阻尼系數(shù)均取0.2。對前三階的模態(tài)振型疊加，提取了梁的根部A點的應(yīng)力及頂部B點的位移結(jié)果，A點和B點的位置示如圖7所示。梁在3個周期內(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，同時由于阻尼的作用，響應(yīng)幅值逐漸衰減，經(jīng)歷7個振蕩后振幅基本衰減到零，每個旋轉(zhuǎn)周期內(nèi)該響應(yīng)重復(fù)1次。以該發(fā)動機(jī)為例，在單次工作時間范圍內(nèi)大約經(jīng)歷40萬次上述循環(huán)過程；由圖9可知，葉片根部的平均應(yīng)力約為295 MPa，在氣流沖擊后的第一個響應(yīng)動應(yīng)力幅值最大，約為155 MPa，然后由于阻尼的作用動應(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 三維仿真求解

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

圖10 流場計算模型及邊界Fig.10 CFD model and boundary

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

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

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

圖12 葉片氣動彎矩分布Fig.12 Distribution of aerodynamic bending moment of blade

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

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

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

3 結(jié)論

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

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

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

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