基于流固耦合的燃機壓氣機葉片結構動力學分析

史久志，覃海起，尹洪，鄧小文, 呂玉坤

(1.廣東電網(wǎng)有限責任公司電力科學研究院，廣東 廣州 510080; 2.華北電力大學 能源與動力工程學院，河北 保定 071000)

基于流固耦合的燃機壓氣機葉片結構動力學分析

史久志1,2，覃海起1，尹洪1，鄧小文1, 呂玉坤2

(1.廣東電網(wǎng)有限責任公司電力科學研究院，廣東 廣州 510080; 2.華北電力大學 能源與動力工程學院，河北 保定 071000)

摘要：采用流固耦合的方法對某壓氣機轉(zhuǎn)子葉片進行氣動彈性數(shù)值模擬計算，將定常計算結果與實驗結果進行對比，驗證計算方法的正確性。在此基礎上，采用單向流固耦合計算對葉片進行結構動力學分析，把非定常計算得出的氣動力載荷轉(zhuǎn)化為結構動力學分析中的載荷壓力場形式，分析葉片在氣動力作用下的瞬態(tài)響應，并進行模態(tài)分析。所得結果表明：葉片在氣動力作用下承受較大的交變應力幅值，容易發(fā)生高周疲勞；產(chǎn)生了較大的動應力，最大應力值位于葉根前緣。通過分析葉片的各階模態(tài)振型和危險點的位置，可以為葉片的疲勞設計分析提供依據(jù)，加強壓氣機轉(zhuǎn)子葉片運行的安全可靠性。

關鍵詞：流固耦合；氣動力；瞬態(tài)響應；模態(tài)分析；高周疲勞

葉片高周疲勞引起的斷裂故障是嚴重影響壓氣機設計質(zhì)量和研制定型周期的關鍵問題之一[1-3]，而作用在葉片上的氣流擾動力是激發(fā)振動的最主要因素。因此，對壓氣機葉片的非定常氣動力分析是很有必要的。

葉片在氣流激振力作用下的動力響應問題比較復雜，它涉及到空氣動力學和結構動力學內(nèi)容。目前分析方法主要有：弱耦合法，包括時域的瞬態(tài)響應法和頻域的穩(wěn)態(tài)響應法[4]；強耦合法，包括順序耦合法和任意拉格朗日坐標系法[5-6]。

弓三偉等[7]借助CFX和ANSYS分析軟件，采用單向順序耦合法對某汽輪機末級大展弦比葉片進行了流固耦合計算，結果表明在非設計狀態(tài)下流場對結構的變形具有更強的敏感性。徐可寧等[8]采用雙向流固耦合法對某壓氣機葉片進行氣動彈性數(shù)值計算，結果表明交替迭代算法是可行的，能得到葉片的瞬態(tài)響應，從而判斷葉片是否發(fā)生顫振。張大義等[9]研究了某葉片的瞬態(tài)動力響應，指出葉片最大變形位于葉尖前緣，最大應力點位于葉根前緣。

本文以某壓氣機葉片為研究對象，首先采用CFX軟件對流體域進行計算流體力學分析，提取葉片表面氣體激振力結果，然后采用有限元分析軟件ANSYS將該壓力載荷加載到固體域表面，對其進行結構動力學分析，獲得葉片的模態(tài)振型，并進行強度計算。最后對該轉(zhuǎn)子進行雙向流固耦合計算，研究葉片的振動特性，為葉片的設計和應用提供分析依據(jù)。

1基本控制方程及數(shù)值算法

本文對NASA Rotor67轉(zhuǎn)子的單通道進行單向流固耦合和雙向流固耦合計算，并對兩種計算結果進行分析。

1.1流體動力學模型

本文計算采用的湍流模型為剪切應力輸運(SSTκ-ω)模型[10]。SSTκ-ω模型是由Menter提出和發(fā)展的一種二方程湍流模型，其在近壁面邊界層采用κ-ω模型，邊界層邊緣和自由剪切層采用κ-ε模型，通過系數(shù)疊加統(tǒng)一成κ-ω形式?？缏曀賶簹鈾C流場存在附面層分離、激波、二次流等復雜流動，近壁面壓力分布的模擬結果直接影響到流固耦合計算的準確性，而該模型能夠充分發(fā)揮κ-ω方程對于邊界層流動計算和κ-ε方程對于自由剪切層計算的優(yōu)勢，進而保證流固耦合計算的準確性。

SSTκ-ω湍流模型方程為：

(2)

式中:κ為湍動能；ω為比耗散率；Γκ、Γω分別為κ和ω的有效擴散率；Gκ為平均速度梯度產(chǎn)生的湍動能；Gω為ω產(chǎn)生的湍流耗散率；Fκ、Fω分別為κ和ω的耗散項；Dω為橫向耗散導數(shù)項；Sκ、Sω為源項;xi、xj都是變量；ρ為密度；νi為速度。

1.2結構動力學模型

多自由度系統(tǒng)動力學方程為

(3)

式中：m為質(zhì)量矩陣；C為結構阻尼矩陣；K為剛度矩陣；x為位移矢量；F為力載荷量，如氣動載荷、離心載荷等。

結構阻尼矩陣采用Rayleigh阻尼，表示為質(zhì)量矩陣和剛度矩陣的線性組合，即

(4)

α和β是兩個系數(shù)，可由兩階較重要的振型決定，取兩階振型阻尼比分別為ζ1和ζ2，動頻分別為ω1和ω2,即

(5)

(6)

得出α和β后，可以確定任意第i階振型的阻尼比

(7)

1.3流固耦合計算

1.3.1單向流固耦合計算

單向流固耦合計算不考慮變形葉片對流場的影響，流場計算得到的氣動壓力載荷在其內(nèi)部進行迭代，當計算達到收斂容差后，把氣動力加載到葉片上，得到了對應時刻下的葉片形變結果，再進行下一個時間步求解。本文研究葉片在氣流激振力作用下的強迫響應，具體方法：將流場分析得到的葉片表面氣動壓力場分布結果加以處理，作為結構瞬態(tài)分析的力邊界條件加載到有限元模型上，以此來求解葉片的瞬態(tài)響應。計算流程如圖1所示。

圖1　單向流固耦合計算流程

1.3.2雙向流固耦合計算

雙向流固耦合計算同時考慮了流場作用對固體結構的變形影響以及結構變形后對流場的反作用，結構變形改變了流體域的邊界及耦合面上流體附面層的厚度和分離等。首先進行流體域求解，得到流場參數(shù)的分布結果，通過提取流體域耦合面上的壓力載荷數(shù)據(jù)并將其加載至固體域的耦合面上，得到了該壓力載荷下的固體形變結果，將固體結構變形后的節(jié)點位移反饋至流體域，修正流體域耦合面附近的網(wǎng)格后再次進行流場計算，如此反復求解。當流體域與固體域計算都達到收斂容差限后，再進行下一時間步的求解。計算流程如圖2所示。

圖2　雙向流固耦合流程

2流固耦合數(shù)值模擬

2.1數(shù)值模擬

基于NASA Rotor67轉(zhuǎn)子具有旋轉(zhuǎn)周期對稱性特性，選取轉(zhuǎn)子全周區(qū)域的1/22(轉(zhuǎn)子葉片有22個葉片)為研究對象，采用專業(yè)網(wǎng)格劃分軟件對流體域模型進行網(wǎng)格劃分，網(wǎng)格節(jié)點數(shù)為63萬，經(jīng)網(wǎng)格無關性驗證，該網(wǎng)格數(shù)量很好地兼顧了計算精度與計算時間的影響。葉片材料參數(shù)及壓氣機設計參數(shù)見表1。

表1材料參數(shù)和設計參數(shù)

名稱數(shù)值楊氏模量E/Pa1.2×1011泊松比0.34密度ρ/(kg·m-3)4850轉(zhuǎn)速n/(r·min-1)16043壓比1.63

流場計算模型及葉片、輪轂表面網(wǎng)格分布如圖3(a)所示。固體域單元類型選取8節(jié)點Solid 185單元，采用自由剖分四面體網(wǎng)格，單元總數(shù)為6 740個，節(jié)點總數(shù)8 218個，固體域計算網(wǎng)格如圖3(b)所示。

圖3　流場及固體計算模型及網(wǎng)格

2.2計算方法的正確性驗證

對NASA Rotor 67轉(zhuǎn)子湍流SSTκ-ω模型做定常計算，計算參數(shù)設置見表2。

表2計算設置參數(shù)

名稱數(shù)值時間步長/s1.7×10-5總時間/s0.017進口總溫度T/K288進口總壓pi/kPa101.325出口平均靜壓po/kPa115.325

將定常計算結果與NASA報告[11]中的實驗結果進行對比分析，所得結果見表3。

表3定常計算結果與實驗結果對比

名稱計算值實驗值相對誤差/%質(zhì)量流量/(kg·s-1)33.6633.660壓比1.611.63-1.2絕熱效率0.880.8780.2

由表3可知，在近設計工況下，定常計算結果與實驗結果相對誤差較小，結果符合良好，證明采用的計算模型和方法是可行的。

3葉片結構動力學分析

3.1模態(tài)分析

葉片根部采用固支約束，分別計算葉片的靜頻和在預應力作用下的動頻，所得結果見表4。

表4葉片前6階靜頻和動頻計算結果

階數(shù)靜頻/Hz動頻/Hz振型1675.47764.33一階彎振21832.11930.7二階彎振32467.42496.9一階扭轉(zhuǎn)44075.04140.9復合振動55125.95195.9二階彎扭66045.46116.5復合振動

將離心載荷加載到葉片的每個節(jié)點單元上，同時將氣動力以預應力的形式加載到葉片表面每個網(wǎng)格節(jié)點上，得到葉片的固有頻率和相應的模態(tài)振型，分析各階振型下葉片的振動情況。圖4為葉片的前六階模態(tài)振型。

圖4　葉片的前6階振型

關注葉片前兩階的頻率，提取前兩階的動頻，求得Rayleigh阻尼系數(shù)分別為0.000 50、0.000 52， 依次將其帶入方程(7)得到各階的阻尼比。

由表4及圖4可知，在氣動力和離心力載荷作用下，葉片的固有頻率變大。葉片第一階和第二階的振動形式為彎曲振動，第三階振動形式為一階扭轉(zhuǎn)振動，第四、第五、第六階振動形式為彎曲和扭轉(zhuǎn)的復合振動。當葉片處于各階共振模態(tài)時，隨著共振階次的增加，振動形式越來越復雜。

3.2單向流固耦合計算及強度分析

如何將非定常計算得到的氣動壓力分布數(shù)據(jù)轉(zhuǎn)化成適合于結構動力學分析中的壓力載荷形式并引入有限元分析，是強度計算中的一個難題。在ANSYS瞬態(tài)分析中，壓力或者任何其他形式的表面力都不能直接加載到葉片單元節(jié)點上，因此必須把氣動計算結果作一定形式的轉(zhuǎn)化。本文首先在固體域表面創(chuàng)建一層面單元，將非定常計算得到的壓力載荷輸出到面單元的每個節(jié)點上，再以預應力的形式加載到葉片的每個單元上，分析葉片在某一時刻的強迫響應。

圖5(a)為某一時刻面單元在流場中的模型，利用該面單元提取激振力并將其加載到葉片表面單元上；圖5(b)為某一時刻葉片在氣動力作用下的壓力分布云圖。

圖5　面單元模型及葉片壓力分布云圖

在氣動力的作用下，葉片會產(chǎn)生位移響應。t時刻葉片的位移、應力分布如圖6所示。

圖6　葉片的位移分布云圖及應力分布云圖

從圖6可看到，t+nΔt為葉片位移和應力最大時對應的時刻。此時葉片最大位移為1.683 mm，位于葉尖前緣處，最大應力為354 MPa，位于葉根前緣處，如圖6(a)所示。將圖6(a)與圖4中的模態(tài)振型進行對比可以看出，葉片的振動形式是一階彎曲振動，最大位移發(fā)生在葉尖前緣處，沿徑向逐漸減小，轉(zhuǎn)子葉片在氣動力作用下承受較大的交變應力幅值如圖6(b)所示，容易發(fā)生高周疲勞；轉(zhuǎn)子葉片在氣流壓力波動下產(chǎn)生了較大的動應力，最大應力值位于葉根前緣。通過分析葉片危險點的位置，可以為疲勞設計提供依據(jù)。

3.3雙向流固耦合計算

通過雙向流固耦合計算，可獲取葉片在流場中的振動特性。在葉片表面近壁區(qū)選取6個監(jiān)測點來監(jiān)測葉片關鍵點的振動位移，節(jié)點分別為葉尖T1、T2、T3和葉片中部(M1、M2、M3點)共6個點，分布如圖7所示。由于壓力面和吸力面對應網(wǎng)格節(jié)點的位移基本相同，因此在弦向上只選取壓力面上的節(jié)點進行監(jiān)測。

圖7　監(jiān)測點的位置示意圖

圖8為雙向流固耦合計算在近設計工況下得到的葉片關鍵點位移響應曲線。圖9為葉片在近失速工況下關鍵點的位移響應曲線，時間步長設置為5.667×10-6s。

圖8　葉片在近設計工況下關鍵點的位移響應曲線

圖9　葉片在近失速工況下關鍵點的位移響應曲線

從圖9可以看出，葉片在阻尼作用下，各個監(jiān)測點的振動周期一致，振幅隨時間逐漸減小趨于某一定值，通過各個測點最后的位移可以描繪出葉片的形狀，得出葉片的振動頻率為356.8 Hz，遠小于葉片的固有頻率，表明該工況下葉片不會發(fā)生共振。

4結論

a) 對壓氣機轉(zhuǎn)子葉片進行分析，證明采用流固耦合計算方法對葉片進行氣動彈性數(shù)值模擬可行。

b) 分析了葉片在預應力作用下的各階振型，當葉片處于各階共振模態(tài)時，隨著共振階次的增加，振動形式越來越復雜。

c) 采用流固耦合計算，可以獲得不同時刻由氣動力激振引起的葉片動應力分布，進一步進行結構動力學分析，得到葉片的瞬態(tài)位移響應和應力響應，可以確定葉片的振動形式以及該工況下應力集中的位置，作為進一步分析的依據(jù)。通過雙向流固耦合計算，可以得到葉片位移響應曲線，在葉片阻尼的作用下，葉片的振動位移幅值隨時間逐漸減小。

參考文獻：

[1] SRINIVASAN A V. Flutter and Resonant Vibration Characteristic of Engine Blades[J]. Engineering for Gas Turbines and Power, 1997, 119(8): 742-775.

[2] KULKARNI A, TISSERANT D, HOSNY D. Turbine Wheel High Cycle Fatigue Reliability Prediction[C] //Proceedings of ASME Turbo Expo 2010. London: ASME, 2010:171-181.

[3] 宋兆泓，陳光，張景武，等.航空發(fā)動機典型故障分析[M]. 北京：北京航空航天大學出版社，1993.

[4] SHIPMAN W C, PALA ZOTTO A N. The Consideration of High Cycle Fatigue in a Turbine Blade [R]. Washington D.C.:AIAA, 2002.

[5] HIRT C W， AMSDEN A A， COOK J L. An Arbitrary Lagrangian-eulerian Computing Method for All Flow Speeds[J]. Journal of Computational Physsics， 1974(14): 227-253.

[6] RAMMAKRISHNAN K, LAWLESS P B, FLEETER S. Finite Element Simulation of Turbo Machine Blade Row Viscous Interactions Vane Vibratory Stress Predictions[R]. Washington D.C.:AIAA，2004.

[7] 弓三偉，任麗蕓，劉火星，等.汽輪機末級轉(zhuǎn)子葉片流固耦合計算[J]. 熱力透平，2007,36(3)：154-157.

GONG Sanwei, REN Liyun, LIU Huoxing, et al. Fluid Solid Interaction of Rotor Blade in Last Stage of Steam Turbine[J]. Thermal Turbine, 2007, 36(3): 154-157.

[8] 徐可寧，王延榮.壓氣機轉(zhuǎn)子葉片的氣動彈性數(shù)值模擬[J]. 航空動力學報，2010,25(10)：2206-2210.

XU Kening, WANG Yanrong. Numerical Simulation of Aero-elastic Response in Compressor Rotor blades[J]. Journal of Aerospace Power， 2010,25(10):2006-2010.

[9] 張大義，馬艷紅，洪杰.氣流激勵下葉片動力響應分析方法[J]. 航空動力學報，2009,24(7)：1523-1527.

ZHANG Dayi， MA Yanhong， HONG Jie. Dynamical Response Analysis on Blades Arisen by Aerodynamic Forces Through Sequential Method [J]. Journal of Aerospace Power, 2009, 24(7):1523-1527.

[10] MENTER F R. Two Equation Eddy Viscosity Turbulence Models for Engineering Applications[J]. AIAA Journal, 1994, 32(8): 1598-1605.

[11] STRAZISAR A J, WOOD J R, HATHAWAY M D, et al. Laser Anemometer Measurements in a Transonic Axial-flow Fan Rotor [R]. Washington D.C.:NASA,1989.

Structural Dynamics Analysis on Gas Turbine Compressor Blade Based on Fluid-solid Coupling

SHI Jiuzhi1,2, QIN Haiqi1, YIN Hong1, DENG Xiaowen1, Lü Yukun2

(1. Electric Power Research Institute of Guangdong Power Grid Co., Ltd., Guangzhou, Guangdong 510080, China; 2. School of Energy and Power Engineering, North China Electric Power University, Baoding, Hebei 071000, China)

Abstract：Fluid-solid coupling method is used for numerical simulating calculation on aeroelasticity of rotor blades of some compressor, and result of steady-state computation and experimental result are compared to verify veracity of this computation method. By means of one-way fluid-solid coupling calculation, this paper carries out structural dynamics analysis on blades which is to transfer aerodynamic load obtained from unsteady computation to the form of load pressure field for structural dynamics analysis. It analyzes transient response of blades under aerodynamic force and conducts modal analysis as well. Results indicate that blades may bear larger amplitude of alternating stress under aerodynamic force and be easy to produce high-cycle fatigue. In addition, larger dynamic stress is produced and the maximum stress occurs at the leading edge of blade root. Analysis on modal shape of each order and locations of dangerous points may provide basis for design and analysis on fatigue of blades as well as strengthen operational security and reliability of blades.

Key words：fluid-solid coupling; aerodynamic force; transient response; modal analysis; high-cycle fatigue

收稿日期：2016-01-22

基金項目：中國博士后科學基金面上資助項目(2015M570696)；中國南方電網(wǎng)有限責任公司科技項目(GDKJ00000005)

doi:10.3969/j.issn.1007-290X.2016.05.001

中圖分類號：V232.4

文獻標志碼：A

文章編號：1007-290X(2016)05-0001-05

作者簡介：

史久志(1990)，男，河南商丘人。在讀碩士研究生，主要從事燃氣輪機等發(fā)電設備研究工作。

覃海起(1990)，男，廣西柳州人。在讀研碩士研究生，主要從事壓氣機穩(wěn)定性研究。

尹洪(1987)，男，重慶人。在站博士后，工學博士，主要從事燃氣輪機等發(fā)電設備研究工作。

(編輯王夏慧)