王朝暉，葉小明

（北京航天動力研究所，北京100076）

0 引言

眾所周知，渦輪內(nèi)部流動的本質(zhì)是非定常的，包含了諸如尾跡、激波、泄露渦、二次流等多種復雜的非定?，F(xiàn)象。而當前的渦輪設計體系基本上都是基于定常假設，無法捕捉到內(nèi)部復雜的非定常流動現(xiàn)象，也就無法更清晰地揭示渦輪內(nèi)部的流動機理。隨著航空航天以及能源工業(yè)的快速發(fā)展，對渦輪性能的要求也越來越高，某些情況下，先前的定常假設體系已經(jīng)略顯欠缺，有必要開展非定常流動現(xiàn)象的研究，期望在深入理解非定常流動機理的基礎上進一步提高渦輪的性能。

早在20世紀80年代，國外就開始了葉輪機械非定常現(xiàn)象的研究，到現(xiàn)在已經(jīng)有許多學者對動靜葉排干涉、葉片非定常力等非定常問題做了大量有價值的研究工作。Dring等人在1982年的研究表明，當靜子和轉(zhuǎn)子之間間距為15%平均軸向弦長時，轉(zhuǎn)子前緣的非定常壓力脈動可高達相對動壓的80%左右[1]，這意味著它在葉片表面會產(chǎn)生很大的非定常載荷，這對葉輪機的性能、葉片強度等都會產(chǎn)生非常大的影響。1987年Giles用Lax-Wendroff時間推進格式求解了非定常歐拉方程組，分析了尾跡/轉(zhuǎn)子之間的干涉效應，定性地預測了氣動噪聲和葉片非定常壓力的產(chǎn)生[2]。1993年Aronoe用完全多重網(wǎng)格技術計算了轉(zhuǎn)子葉排間的流動[3]。此外，國內(nèi)方面也有學者對葉輪機械的非定?，F(xiàn)象進行了一系列的研究[4-6]。

液體火箭發(fā)動機渦輪的特點是小尺寸、小流量、高轉(zhuǎn)速、高壓力（可達幾十兆帕）、小展弦比和無冷卻，與航空發(fā)動機用的大展弦比、大流量的渦輪有著明顯的差別。過去眾多學者的研究大都集中在航空發(fā)動機渦輪或燃氣輪機，不論是設計、優(yōu)化還是內(nèi)部流動現(xiàn)象，對液體火箭發(fā)動機渦輪的研究并不多見，尤其對液體火箭發(fā)動機渦輪的非定常流動現(xiàn)象的研究更加少見。本文基于液體火箭發(fā)動機渦輪這個新的平臺，對其非定常流動現(xiàn)象進行深入的研究與分析。

1 計算網(wǎng)格

為了適應葉片通道空間拓撲結(jié)構(gòu)并保證較高的網(wǎng)格質(zhì)量，對計算域采用多塊結(jié)構(gòu)化網(wǎng)格技術。建立渦輪模型時考慮了輪緣輪轂處的加工倒角，考慮了密封，網(wǎng)格最小正交角度21°。表1給出了計算網(wǎng)格分布情況，圖1給出了子午流面網(wǎng)格劃分和密封位置網(wǎng)格的局部放大圖。圖2給出了計算結(jié)果得到壁面y+分布云圖，可見壁面y+最大值小于10，與所選用的Spalart-Allmaras湍流模型相一致。

表1 計算網(wǎng)格分布Tab.1 Distribution of computational grid

圖1 子午流面網(wǎng)格劃分圖Fig.1 Grid generation on meridional stream surface

圖2 壁面y+分布云圖Fig.2 Distribution of y+at wall

2 求解方法

定常計算采用混合平面法，周向平均混合。非定常計算采用滑移面方法，交界面線性插值。本算例渦輪葉片數(shù)是23，50，31，48，為了減少非定常計算量，采用1989年Rai提出的區(qū)域縮放法[7]對轉(zhuǎn)靜子葉片數(shù)進行適當?shù)恼{(diào)整，調(diào)整后為25，50，25，50，約化為1:2:1:2，調(diào)整葉片數(shù)的同時對葉型進行相應比例的縮放，以保證節(jié)弦比和堵塞度不變。

物理時間步長的選擇起著非常重要的作用，在求解非定常流動情況時，時間步長太小，會使得計算量變大，推進時間變長，物理時間步長過大，又會導致不能捕捉到所關心變量的準確的變化情況，因此合適的物理時間步長是非定常計算中一個非常重要的參數(shù)。在葉輪機的非定常計算中，一般每個動葉柵距至少取10個物理時間步長以上。把動葉轉(zhuǎn)過一個靜葉柵距的時間定義為一個計算周期T，每個計算周期T選取20個計算點，計算25個周期，共500個物理時間步。

本算例為亞聲速流動情況，進口給定氣流總溫、總壓，軸向進氣，出口給定靜壓，壁面條件絕熱、無滑移。工質(zhì)選用氣態(tài)H2，cp=15 850 J/(kg·K)，γ=1.384，給出粘性系數(shù)隨溫度的變化曲線。湍流模型選擇魯棒性較好且收斂速度較快的Spalart-Allmaras模型。非定常計算采用Jameson提出的隱式雙重時間推進法[8]，并結(jié)合多重網(wǎng)格法、當?shù)貢r間步長和隱式殘差光順等加速收斂技術。空間離散采用二階中心差分，并添加二階和四階人工粘性項。虛擬時間推進采用四階龍格-庫塔方法，每個虛擬時間步內(nèi)均迭代計算50步。在進行非定常計算之前先進行定常求解，然后把定常計算的結(jié)果作為非定常計算的初始解。

為了驗證網(wǎng)格無關性，在進行非定常計算之前先對不同數(shù)量的網(wǎng)格進行穩(wěn)態(tài)計算，對比網(wǎng)格數(shù)量和計算效率的變化趨勢。發(fā)現(xiàn)當網(wǎng)格數(shù)達到182萬時，效率值基本無變化，所以接下來的計算都采用這一數(shù)量的網(wǎng)格進行。

3 結(jié)果分析

3.1 定常與非定?？傮w性能參數(shù)比較

為了了解模化（即葉片數(shù)和相應葉型的調(diào)整）對計算結(jié)果的影響，把?；昂蟮娜~型都進行了定常計算。原始葉型的定常計算結(jié)果、?；~型的定常和非定常時均結(jié)果的對比見表2。

表2 定常和非定常結(jié)果比較Tab.2 Comparison of steady and unsteady results

比較?；昂蟮亩ǔ＝Y(jié)果可以看出，在壓比相同的情況下，?；笮试黾恿?.13%，流量減小了0.18%，功率減小了0.4%，所有變化都小于0.5%。總體而言，所作的葉片?；瘜u輪的數(shù)值模擬結(jié)果的影響很小，基本可以忽略，可以認為所作的模化是合理的。

下面著重對比分析模化后定常與非定常計算結(jié)果的差別。表2顯示，?；蠖ǔＳ嬎愕馁|(zhì)量流量結(jié)果比非定常稍大，但都高于設計值。非定常時均效率要比定常結(jié)果高出0.4個百分點，這是因為定常計算在動靜交界面采用了周向平均模型，使本來不均勻的葉片排出口流場人為地混合均勻，增加了摻混損失。此外定常計算的流量高于非定常的結(jié)果，這主要是因為相同的出口背壓下，非定常出口尾跡的低能流體流到后面的通道，使之產(chǎn)生了一定程度的堵塞，流量降低，使得前排葉柵進口氣流攻角增大，進一步降低了流量，而定常計算情況下相應的通流能力較強。上述對比說明定常計算的結(jié)果只能是一定程度的近似，并不能真正反映渦輪的性能，想要得到渦輪的真實流動情況，還需要對其進行非定常計算。渦輪試驗是在西安交大葉輪機械研究所微燃機渦輪試驗臺上進行的，試驗渦輪為真實產(chǎn)品即?；皽u輪，由于試驗臺很難達到設計的高入口壓力，因此渦輪試驗是依據(jù)相似理論采用的模化試驗。圖3給出了渦輪性能參數(shù)試驗結(jié)果與數(shù)值結(jié)果的對比。

圖3 性能曲線Fig.3 Performance curves

通過改變渦輪出口壓力，計算了設計工況點之外的10個工況點。對圖3給出的效率特性曲線對比分析，發(fā)現(xiàn)試驗結(jié)果顯示最高效率為75.88%，出現(xiàn)在設計點（即落壓比1.6） 處。對比定常效率和試驗效率，定常效率比試驗效率要高出約2個百分點，這是由于數(shù)值計算簡化了密封模型，認為工質(zhì)為理想氣體，因此得出的結(jié)果比試驗值大。不論是試驗結(jié)果還是數(shù)值結(jié)果，在額定轉(zhuǎn)速下渦輪效率隨落壓比均變化緩慢，且均維持在較高水平，說明該渦輪實際運行性能較好，在設計工況下具有良好的氣動性能表現(xiàn)。另外，從圖3還可以看出，在設計點工況，定常和非定常結(jié)果差別很?。辉谄x設計點工況，定常和非定常之間的差別開始變得明顯。

分析圖3給出的質(zhì)量流量隨落壓比變化曲線對比情況，發(fā)現(xiàn)數(shù)值計算結(jié)果和原設計值均比試驗結(jié)果大，這是因為渦輪試驗是一個?；囼?，試驗渦輪本身的流量及功率就非常小，軸系摩擦等阻力因素占的比重相對較大，在性能參數(shù)結(jié)果的后處理時相對保守而且存在一定的試驗誤差所致。定常和非定常結(jié)果非常接近，尤其在設計點處兩者差別非常小，在偏離設計點時開始稍有差別。從曲線的趨勢來看，不論是數(shù)值結(jié)果還是試驗結(jié)果，質(zhì)量流量都是隨著落壓比變大而增加且增加趨勢逐漸變緩并趨于穩(wěn)定值。這不難從理論上解釋，渦輪幾何形狀不變，進口總溫總壓不變，隨著出口壓力的減小，一級靜葉流速增加，流量變大，在還未達到臨界狀態(tài)前，流量隨著出口壓力的降低而增加，當達到或超過臨界狀態(tài)后，不論如何降低出口壓力，渦輪的流量都維持不變。

3.2 葉片表面非定常力分析

選取計算物理時間步的最后400步，得到監(jiān)測點的壓力隨時間變化情況。一級靜葉后緣吸力面一點（相對軸向弦長85%處，50%span） 壓力隨計算時間步的變化曲線如圖4所示，快速傅里葉變換得到頻譜圖5。

圖4 一級靜葉吸力面某點壓力變化曲線Fig.4 Static pressure at a point on suction surface of first-stage static stator

從圖5可以看出最大振幅對應頻率為37 916.67 Hz，幅值4.5 kPa，除以當?shù)貕毫ζ骄?，波動百分比?.57%。這說明一級靜葉的流動還是比較穩(wěn)定的，沒有出現(xiàn)明顯的非定常波動，靜葉尾緣部分的微小波動主要是因為受到了交界面下游動葉前緣附近的勢流場影響。從圖5中可以看到有4個峰值頻率，其中29 166.67 Hz和58 333.33 Hz的波動分別是由靜葉排和動葉排的周期性排列引起的，其余的峰值頻率明顯是由動靜葉排間相互干涉、葉排尾跡的非定常變化以及諸如二次渦流等非定常因素引起的。這些頻譜信息只有在非定常計算中才可能得到，在定常計算中是不能獲得這些頻譜信息的。

圖5 頻譜圖Fig.5 Distribution of frequencies

圖6 給出了一級動葉吸力面前緣某點（相對軸向弦長15%位置，50%span） 的壓力隨計算時間步的變化情況。可以看到，該點的壓力呈現(xiàn)出明顯的周期性波動，并且呈現(xiàn)出類似諧波函數(shù)的變化特點，說明非定常計算可以反映出渦輪內(nèi)部真實的流動情況。壓力最大值與最小值相差約90 kPa，波動范圍達1.2%，表現(xiàn)出較強的非定常現(xiàn)象，說明葉片會受到非定常激振力。

圖6 一級動葉吸力面某點壓力變化曲線Fig.6 Static pressure at a point on suction surface of first-stage moving rotor

對曲線進行快速傅里葉變換，得到該點靜壓的頻譜圖7，峰值頻率為29 166.67 Hz，這與葉排通過頻率BPF是相一致的。

圖7 頻譜圖Fig.7 Distribution of frequencies

在一級靜葉和一級動葉排相互干擾中，葉排通過頻率BPF與動靜葉片數(shù)及動葉轉(zhuǎn)速有關。葉排通過頻率BPF是一片動葉通過一個靜葉通道所需要時間的倒數(shù)：

式中：Nstator為靜葉排的葉片數(shù)；Ω為轉(zhuǎn)速。按照原始的一級靜葉數(shù)目23來說，計算可得葉片通過頻率BPF=26 833.33 Hz，這與頻譜圖上的峰值頻率并不吻合，兩者之間相差約2 300 Hz。這是因為在計算過程中人為更改了一級靜葉葉片數(shù)目，由23改為了25，同時對一級靜葉進行了縮放，由此造成了流動參數(shù)頻譜的移動。從上述結(jié)果可以說明，區(qū)域縮放法雖然可以在很大程度上減少計算量，但是由于人為地改變了葉片數(shù)目，將會導致流動參數(shù)頻譜的微小移動，因此在利用區(qū)域縮放法研究葉片非定常氣動力、葉片疲勞和葉片顫振等方面的問題時，需要注意縮放前后相應流動參數(shù)對應頻譜圖的微小偏移。

葉片力定義為葉片吸力面和壓力面所受到壓力引起的合力，葉片力除以平均值得到葉片力系數(shù)。圖8描繪了一級動葉葉片力圍繞其平均值的波動情況，由于動靜葉排的相對周向運動，前排尾跡周期性地掃過動葉，尾跡與邊界層的相互作用周期性地改變著動葉表面的壓力分布情況，進而使得葉片受力也跟著發(fā)生周期性波動。

圖8 一級動葉葉片力系數(shù)隨時間變化曲線Fig.8 Variation of blade force coefficient of first-stage moving rotor with time

從圖8可看出葉片力呈現(xiàn)類似簡諧波的、顯著的、有規(guī)律的周期性波動，葉片力圍繞平均值的波動幅值達15.78%，波動的頻率為29 166.67 Hz，和動葉片通過靜葉柵距的頻率相一致。

動葉葉片所受到的沿轉(zhuǎn)軸方向的力矩直接決定著渦輪的做功能力。葉片力矩定義為葉片力周向分量對轉(zhuǎn)軸的力矩，葉片力矩除以平均值得到葉片力矩系數(shù)。通過非定常計算考察動葉在旋轉(zhuǎn)過程中所受到的沿轉(zhuǎn)軸方向的力矩變化情況對理解和評估渦輪轉(zhuǎn)子做功性能和穩(wěn)定性等將是非常有幫助的。圖9給出了一級動葉葉片力矩系數(shù)隨時間變化曲線。

圖9 一級動葉葉片力矩系數(shù)隨時間變化曲線Fig.9 Variation of blade torque coefficients of first-stage moving rotor with time

從圖9可見，力矩的變化趨勢和葉片力的情況類似，波動幅值稍大一些，達到了17.31%。不論是力還是力矩的波動，都對渦輪整體性能有一定的影響，過于強烈的葉片力波動將會對葉片的強度和壽命造成影響，嚴重時可能造成葉片的損壞，不穩(wěn)定的力矩會影響渦輪功率的穩(wěn)定輸出。

二級動葉吸力面某一點（相對軸向弦長15%處，50%span）壓力變化情況如圖10所示。頻譜圖如圖11所示，頻譜圖上峰值頻率29 166.67 Hz的波動源自二級靜葉的影響。從波動曲線可以看到該點的最大靜壓差為88 kPa，波動幅值達1.44%。圖12和圖13給出了二級動葉葉片力和力矩系數(shù)的變化情況，可以看出，與一級動葉類似，力和力矩系數(shù)均呈現(xiàn)明顯的規(guī)律性的周期性波動。二級動葉葉片力的大小變化幅度5.9%，二級動葉沿轉(zhuǎn)軸方向的力矩波動幅值8.23%，波動頻率為29 166.67 Hz，與動葉通過靜葉柵距頻率相一致。與一級動葉相比，二級動葉葉片力和力矩的波動幅值略小。

圖10 二級動葉吸力面某點壓力變化曲線Fig.10 Static pressure at a point on suction surface of second-stage moving rotor

圖11 頻譜圖Fig.11 Distribution of frequencies

圖12 二級動葉葉片力系數(shù)隨時間變化曲線Fig.12 Variation of blade force coefficient of second-stage moving rotor with time

圖13 二級動葉葉片力矩系數(shù)隨時間變化曲線Fig.13 Variation of blade torque coefficients of second-stage moving rotor with time

通過分析可以發(fā)現(xiàn)本渦輪算例的葉片力和力矩的波動較為明顯。不穩(wěn)定現(xiàn)象產(chǎn)生的根源主要是渦輪內(nèi)部存在著周向不均勻性以及由于端壁二次流、通道渦等因素。如何減小這種非定常因素是降低渦輪非定常葉片力和力矩的關鍵，由于本渦輪葉片葉高相對較小，屬于小展弦比渦輪，因此端壁二次流在主流中所占的比重相對較大，可以考慮采用端壁筋方法抑制端壁二次流的產(chǎn)生和發(fā)展，即在葉柵通道的端壁面上沿流向設置一條高度和邊界層厚度相當?shù)募毿〉耐蛊穑宰柚贡诿孢吔鐚恿黧w從壓力面到吸力面的橫向流動，達到減少端壁二次流的目的，從而減弱渦輪內(nèi)部的非定常因素。另外，還可以充分論證各葉排間的軸向間距，篩選出最優(yōu)的葉排軸向間距，也可以起到減弱渦輪非定常效應的作用。

葉片表面壓力分布是衡量葉片氣動性能的重要參數(shù)，它不僅決定了葉片作功能力，而且在一定程度上會影響葉片的強度和壽命等。一級動葉處于兩排靜子中間，受到尾跡與勢流的雙重干擾，非定常性較強。圖14給出了各排葉片表面壓力變化情況的對比。通過對比可見，一級動葉的壓力波動幅度最顯著，與上文結(jié)論一致。從圖14（b）一級動葉的表面壓力曲線可以看出，靜子尾跡被轉(zhuǎn)子葉片切割時，在一級動葉葉片前緣有較大幅度的壓力擾動。這表明尾跡明顯改變了來流氣流角，導致一級動葉葉片前緣承受著較大幅度的非定常負荷。尾跡在向下游運動的過程中，對動葉吸力面的影響比對壓力面的影響更為強烈，而且一級動葉吸力面的波動范圍和波動幅度都比壓力面更為強烈。雖然壓力面在尾跡掃過時同樣引起表面負荷波動，但相對于吸力面而言，其脈動幅度較小，主要是因為壓力面前緣附近的逆壓梯度小，大部分區(qū)域為順壓梯度，流動保持加速流動狀態(tài)。這種情況下，邊界層流動對外界擾動并不敏感，尾跡的影響相對較弱。

圖14 一個計算周期T內(nèi)葉片表面壓力變化情況Fig.14 Pressure distribution on blade surface in one calculation period T

渦輪的另一個非常重要的非定常特性就是葉片進口氣流角會發(fā)生周期性的變化，這主要是由于前排葉柵通道出來的不均勻流周期性地經(jīng)過下游葉片入口所引起的。圖15給出了中葉展截面動葉進口氣流角的變化情況。

圖15 不同時刻中葉展截面動葉進口氣流角沿柵距的分布Fig.15 Distribution of inflow angle at mid-span section of rotor at different time

從圖15可以看到，一級動葉進口氣流角在一個周期內(nèi)變化幅度大約為5°左右，二級動葉進口氣流角波動幅度達10°左右。由此可見非定常因素在本算例中是非常明顯的。

4 結(jié)論

通過對液體火箭發(fā)動機兩級軸流氫燃料渦輪進行非定常數(shù)值模擬，著重分析了渦輪動靜葉排干涉、非定常葉片力以及渦輪性能曲線等問題，得到以下結(jié)論：

1）區(qū)域縮放法可以有效地計算渦輪內(nèi)部非定常流動現(xiàn)象，通過快速傅里葉變換可以獲得到流場中壓力波動的頻譜信息。計算結(jié)果表明：本算例中一級動葉葉片力和力矩存在較大的非定常波動，波動幅值達15.78%和17.31%。

2）定常計算忽略了渦輪內(nèi)部流動周向不均勻性等非定常因素，與真實流動情況有一定差別，不能捕捉到葉排間干涉和葉片非定常力等現(xiàn)象，無法更加準確全面地反映渦輪的性能。某些情況下，只有定常計算還不夠，還需要對其進行非定常計算作為補充和校核。

3） 在非定常條件下，渦輪靜子尾跡和位勢作用將直接影響并改變下游動葉葉片表面負荷分布，特別是動葉前緣和吸力面擴壓段等位置受到的影響更為強烈，從而對渦輪性能產(chǎn)生較為顯著的影響。非定常因素導致了動葉進口氣流角發(fā)生較大的波動，一級和二級動葉進口氣流角波動范圍分別可達5°和10°左右。

4） 就總體性能參數(shù)而言，在設計點工況，定常和非定常結(jié)果相差不大；在偏離設計點工況，差別開始變得稍加明顯，而非定常結(jié)果更加符合實際情況。

[1]DRING R P,JOSLYN H D,HARDIN L W,et al.Turbine rotor-stator interaction,ASME 82-GT-3[R].USA:ASME,1982.

[2]GILES M B.Calculation of unsteady wake/rotor interactions,AIAA-87-0006[R].USA:AIAA,1987.

[3]ARONOE A.Viscous analysis of three-dimensional rotor flows using a multigrid method,NASA-TM-106266[R].USA:NASA,1993.

[4]鄒正平,葉建,張永新,等.非定常流動對葉片表面負荷分布影響的數(shù)值模擬研究[J].燃氣渦輪試驗與研究,2006,19(1):21-26.

[5]祈明旭,豐鎮(zhèn)平.動靜干涉效應對軸流透平級氣動性能的影響[J].工程熱物理學報,2003,24(1):39-42.

[6]吳先鴻,陳懋章.非定常動靜葉干涉作用的二維數(shù)值模擬[J].航空動力學報,1998,13(2):123-128.

[7]RAI M M.Three-dimensional Navier-Stokes simulations of turbine rotor-stator interaction:Methodology[J].Journal of Propulsion and Power,1989,5(3):305-311.

[8]JAMESON A.Time dependent calculations using multi-grid with applications to unsteady flows past airfoils and wings,AIAA 91-1596[R].USA:AIAA,1991.