基于動網(wǎng)格的渦輪流場仿真方法及特性分析

徐學(xué)文，曲 凱，白 玉，2

（1.海軍航空大學(xué)，山東煙臺 264001；2. 61035部隊(duì)，北京 102205）

渦輪是直接帶動渦輪噴氣發(fā)動機(jī)壓氣機(jī)旋轉(zhuǎn)對空氣做功的重要部件。近年來，隨著無人機(jī)、導(dǎo)彈等裝備快速發(fā)展，渦輪噴氣發(fā)動機(jī)的應(yīng)用越來越廣泛，因此，對渦輪特性的研究也越來越多[1]。當(dāng)前，對渦輪流場特性分析主要基于實(shí)驗(yàn)的方法和氣體動力學(xué)分析[2-4]，但是，實(shí)驗(yàn)成本高昂，相似條件較難滿足，實(shí)驗(yàn)結(jié)果可參考性不大[5-6]，而氣體動力學(xué)分析也只能做簡單的量化計(jì)算和定性分析[7-10]。隨著計(jì)算機(jī)仿真技術(shù)發(fā)展，通過流場仿真計(jì)算來研究渦輪流場特性的越來越多，由于建模和計(jì)算的復(fù)雜性，以及三維仿真成本高昂[11-12]，故大多采用一維、二維流場穩(wěn)態(tài)仿真計(jì)算[1，13-14]，從而提高渦輪流場特性分析效率。本文基于動網(wǎng)格的二維流場仿真技術(shù)，探討渦輪旋轉(zhuǎn)條件下的流場仿真方法，計(jì)算渦輪流場參數(shù)變化，開展渦輪流場特性分析。

1 仿真計(jì)算模型

1.1 流場模型

本文選擇單級軸向渦輪作為研究對象，它由一級導(dǎo)向器（靜子）和一級渦輪（轉(zhuǎn)子）組成，從外界進(jìn)來的高壓氣體首先經(jīng)過導(dǎo)向器高速膨脹，在出口處獲得足夠速度后，直接沖擊在渦輪葉片上，推動渦輪旋轉(zhuǎn)做功，最后從渦輪出口高速噴出。

氣體在導(dǎo)向器和渦輪葉柵通道流動過程中，速度、壓力與溫度沿軸向不斷發(fā)生變化，而沿徑向受力變化較小[15]，因此，為簡化計(jì)算，建立二維流場仿真模型，假設(shè)流場參數(shù)沿葉片徑向均勻分布，選取導(dǎo)向器、渦輪各個(gè)葉片周邊氣體作為研究對象，并按照葉片型線分割流體，作為流場周期邊界，并適當(dāng)延長氣體進(jìn)口/出口邊界長度，如圖1所示。

圖1 渦輪二維流場仿真模型Fig.1 Two-dimensional flow field simulation model of turbine

1.2 流場計(jì)算方程

由于渦輪葉片是運(yùn)動的，因此，這里選用守恒型動網(wǎng)格流場計(jì)算方程[16]，對于邊界移動的任意控制體積V上的標(biāo)量φ（質(zhì)量ρ、速度u、能量E）的守恒型方程為：

式（1）中：V(t)為空間中大小和形狀都隨時(shí)間變化的控制體積；ρ為流體密度；?V(t)為控制體積的運(yùn)動邊界；ug為運(yùn)動網(wǎng)格的運(yùn)動速度；u為流體速度矢量；Γ為耗散系數(shù)；Sφ是標(biāo)量φ的源項(xiàng)。

湍流模型采用計(jì)算精度比較高的且應(yīng)用比較廣泛的k-ε二方程模型：

k控制方程：

1.3 邊界條件

本文選取渦輪在ω=1 500 r/min 轉(zhuǎn)速下的仿真計(jì)算，仿真模型邊界條件，如圖2所示。

圖2 渦輪二維流場計(jì)算網(wǎng)格及邊界條件Fig.2 Calculation grid and boundary conditions of two-dimensional flow field of turbine

3）葉片壁面條件

假設(shè)氣體與渦輪葉柵通道壁面之間均為“無滑移”條件，壁面粗糙度為常數(shù)0.5。

4）渦輪—導(dǎo)向器界面條件

由于渦輪在發(fā)動機(jī)工作過程中以ω速度運(yùn)轉(zhuǎn)，對應(yīng)的流場仿真區(qū)域——渦輪葉片及其流場區(qū)域以圓周切向速度Vt上下做周期運(yùn)動，由此在導(dǎo)向器出口與渦輪進(jìn)口之間的燃?xì)饬魍ㄍǖ罊M截面時(shí)刻發(fā)生變化。為保證流場仿真過程中燃?xì)赓|(zhì)量守恒，在導(dǎo)向器出口與渦輪進(jìn)口之間采用渦輪—導(dǎo)向器間滑動邊界條件[16]，即在渦輪流場區(qū)域與導(dǎo)向器流場區(qū)域的邊界重合區(qū)域定義為內(nèi)部區(qū)域，邊界上其他區(qū)域（非重合區(qū)域）為周期區(qū)域，如圖2所示。渦輪與導(dǎo)向器間滑動

渦輪中氣體流動守恒就是從渦輪與導(dǎo)向器間重疊區(qū)域界面計(jì)算得到的。

5）周期邊界條件

不管是導(dǎo)向器還是渦輪，其上葉片均為均勻分布，因此，這里選取各1個(gè)葉片及其周邊流場區(qū)域作為仿真模型，葉片周邊流場兩側(cè)邊界設(shè)定為周期邊界。

2 動網(wǎng)格流場仿真方法

本文采用有限體積法[17]的無結(jié)構(gòu)化網(wǎng)格離散導(dǎo)向器—渦輪葉片周邊流場區(qū)域，由于本文渦輪與導(dǎo)向器間滑動速度與氣體進(jìn)口速度均比較低，為提高計(jì)算精度，僅在流場壁面、周期邊界、滑動界面采用加密網(wǎng)格，最終劃得網(wǎng)格單元數(shù)為7 917，最小單元面積為6.554×10-4m2，最大單元面積為2.709×10-2m2，如圖2所示。

在瞬態(tài)仿真計(jì)算中，流場控制方程的時(shí)間導(dǎo)數(shù)項(xiàng)用一階向后差分公式得到：

式（11）中，δVj是控制體上面j在時(shí)間步Δt內(nèi)掃出的體積。

3 仿真結(jié)果分析

本文仿真過程分2步：第1步，設(shè)定滑動邊界條件的滑動速度為0，應(yīng)用二階迎風(fēng)格式耦合算法求解流場控制方程，得到流場穩(wěn)態(tài)仿真結(jié)果，作為第2步瞬態(tài)流場仿真的初始化值；第2步，啟動導(dǎo)向器—渦輪間滑動邊界條件，采用二階隱式耦合算法計(jì)算瞬態(tài)流場控制方程。渦輪轉(zhuǎn)速不同，時(shí)間步長也不一樣，這里根據(jù)渦輪葉片通過1 個(gè)導(dǎo)向器葉柵通道跨距Δy的時(shí)間T來確定時(shí)間步長：

因此，渦輪葉片需要用18個(gè)時(shí)間步通過一個(gè)導(dǎo)向器葉柵通道跨距，這里選取流場仿真時(shí)間步長Δt為0.000 1 s，仿真時(shí)間為0.008 s，80 個(gè)時(shí)間步長，選取0.002 5～0.006 1 s 間計(jì)算結(jié)果進(jìn)行分析，研究流場變量沿軸向位置變化趨勢，每個(gè)流場變量在橫截面上取平均值。

0.002 5 s 時(shí)，導(dǎo)向器流場區(qū)域與渦輪流場區(qū)域在滑動界面處重合，如圖3 a）所示。

沿著流場周期邊界線（導(dǎo)向器葉片和渦輪葉片軸向長度段）葉柵流場通道橫截面逐漸減小，壓力變化呈現(xiàn)先降低后升高，在導(dǎo)向器流場出口處達(dá)到最大值；在渦輪葉柵通道流動中，壓力持續(xù)降低，低過出口氣體壓力P2，然后受P2的影響，增大到與外界壓力平衡，如圖3 b）所示。

氣體在導(dǎo)向器葉柵通道流動過程中，速度馬赫數(shù)持續(xù)增大，在滑動界面處，導(dǎo)向器流場出口速度等于渦輪流場進(jìn)口速度，如圖3 c）所示，在渦輪葉柵通道流動過程中，先呈現(xiàn)繼續(xù)升高，而后由于受出口氣體壓力P2（背壓）影響，氣體速度呈現(xiàn)降低的變化趨勢。

圖3 0.002 5 s時(shí)渦輪流場位置及參數(shù)變化曲線Fig.3 Turbine flow field position and parameter variation curve at 0.002 5 s

0.004 0 s時(shí)，渦輪向下移動了0.161 6 m，如圖4 a）所示。

氣體在導(dǎo)向器葉柵通道流動過程中壓力持續(xù)下降，在滑動界面處，達(dá)到最低值，由于受渦輪葉片阻力影響（導(dǎo)向器與渦輪流場出現(xiàn)非重合區(qū)域），氣體壓力在渦輪進(jìn)口處出現(xiàn)躍升，如圖4 b）所示，這是由于計(jì)算過程中保持氣體參量守恒[18]引起的，而不是由于激波引起的。

氣體壓力在渦輪葉柵通道中持續(xù)降低，受P2的影響出現(xiàn)壓力波峰振蕩，并且導(dǎo)向器與渦輪流場出現(xiàn)非重合區(qū)域越大，振蕩波峰越向上游擴(kuò)展。氣體速度在滑動界面處也出現(xiàn)不連續(xù)現(xiàn)象，并且在渦輪葉柵通道中也出現(xiàn)速度振蕩，如圖4 c）所示。這些流場參數(shù)在葉柵通道的變化趨勢與一維等熵流動理論分析結(jié)果[19]是一致的。

圖4 0.004 s時(shí)轉(zhuǎn)子流場位置及參數(shù)變化曲線Fig.4 Turbine flow field position and parameter variation curve at 0.004 s

0.004 2 s 時(shí)，渦輪葉片到達(dá)最下端，0.004 3 s 時(shí)，渦輪葉片離開導(dǎo)向器葉片流場區(qū)域，下一個(gè)渦輪葉片進(jìn)入導(dǎo)向器葉片流場區(qū)域。

4 結(jié)論

本文基于動網(wǎng)格二維流場仿真技術(shù)，采用流體動力學(xué)流場控制方程，應(yīng)用有限元體積法仿真計(jì)算了航空發(fā)動機(jī)單級渦輪某一轉(zhuǎn)速下氣體流場分布，并對渦輪流場特性進(jìn)行了分析，仿真過程及結(jié)果表明：

1）為提高流場仿真計(jì)算的收斂性和結(jié)果的精度，基于動網(wǎng)格的渦輪流場仿真分2步進(jìn)行：首先，計(jì)算穩(wěn)態(tài)條件下的流場仿真結(jié)果，將其作為下一步的初始化值；第2步，啟動某一轉(zhuǎn)速下渦輪轉(zhuǎn)動引起的瞬態(tài)流場仿真。

2）在渦輪轉(zhuǎn)動過程中，轉(zhuǎn)子—靜子間滑動界面時(shí)刻發(fā)生變化，引起界面處流場仿真參數(shù)發(fā)生階躍，這是由仿真計(jì)算而非激波引起的。

3）氣體在轉(zhuǎn)子—靜子葉柵通道流動過程中，流場參數(shù)變化主要受葉柵通道橫截面積和進(jìn)/出邊界條件影響，仿真結(jié)果與一維等熵流動理論分析結(jié)果一致。