唐明智 ，金東海 ，2，郭 昕 ，桂幸民 ，2

（1.北京航空航天大學(xué)能源與動力工程學(xué)院，北京100191；2.先進航空發(fā)動機協(xié)同創(chuàng)新中心，北京100191；3.中國航發(fā)四川燃氣渦輪研究院，成都610500）

0 引言

現(xiàn)階段3D數(shù)值模擬在風(fēng)扇/壓氣機設(shè)計中得到了廣泛使用。然而在設(shè)計的初級階段，通流模型由于一方面可以快速提供風(fēng)扇/壓氣機的性能及內(nèi)部流動特性，另一方面也易于將經(jīng)驗參數(shù)納入到數(shù)值模擬中[1-2]，因此仍具有重要作用[3-4]。目前應(yīng)用最廣泛的通流方法是流線曲率法，而基于周向平均Euler[5]或Navier-Stokes（N-S）方程[6-8]的通流模型也同樣被廣泛研究，并應(yīng)用于壓氣機的特性預(yù)測和流場模擬。這些通流模型適用于跨聲壓氣機的流場模擬和性能預(yù)測[9]。

在對N-S方程進行周向平均的過程中，不可避免地產(chǎn)生一些附加項，其中與表面力相關(guān)的項被模化為無黏葉片力項和葉片黏性項，而對于由于流場的周向不均勻性及控制方程的非線性所引發(fā)的高階項，即周向脈動應(yīng)力項，早期研究認為其對流場的影響可以忽略，即假設(shè)葉片通道進口前流動參數(shù)沿周向均勻分布，這主要是由于在早期研究中，相對于控制方程中的其他項，反映周向不均勻性的項的量級很??；另一個原因則是設(shè)計者缺乏模型來預(yù)測這些項的值[10]。然而，隨著當(dāng)前壓氣機負荷的提高，以及彎掠葉片技術(shù)的廣泛應(yīng)用，周向不均勻性的影響正被逐步認識到，有助于更為精確地應(yīng)用彎掠葉片技術(shù)。事實上，這些項的作用可能會比黏性項的作用更大，并且隨著負荷的增大，其作用也會增強[11]。此外，周向不均勻性會影響流動參數(shù)的軸向和展向的分布[9，12-13]，并能反映角區(qū)失速和徑向摻混等現(xiàn)象。對于掠葉片，葉片通道進口的流動平衡也會被周向不均勻性所影響，并會誘導(dǎo)流動參數(shù)的重新分配，而引發(fā)周向不均勻性的1個主要來源就是無黏葉片力[14]。

為了?；芟虿痪鶆蛐?，通常的方法是計算周向脈動應(yīng)力項，其中1種方法是開展S1流面上的計算，然后提取出周向脈動應(yīng)力項。吳仲華采用中心流線法來簡化對S1流面的計算[15]，并在之后得到了一定的發(fā)展[16-17]。由于周向脈動應(yīng)力項具有諧波特征[18]，因此非線性諧波平衡法也被用來?；芟虿痪鶆蛐訹19]，不過這種方法并未獲得廣泛應(yīng)用。

文獻[20]通過應(yīng)力輸運模型對周向脈動應(yīng)力項進行建模，分析了具有不同掠角的掠葉柵在0°迎角下流場中周向不均勻性的影響，結(jié)果表明周向不均勻性會改變掠葉片進口流動參數(shù)，并引發(fā)進口新的徑向平衡。本文在此基礎(chǔ)上應(yīng)用該應(yīng)力輸運模型分析迎角改變時掠葉柵進口流場中周向不均勻性的作用。

1 控制方程

1.1 周向平均通流模型

在相對柱坐標系中對N-S方程進行周向平均[10]，并將周向脈動應(yīng)力項整合為周向脈動源項P，可得到通流模型的主控方程

方程中各項定義為

式中：U為守恒量；F和G為對流（無黏）通量；Fv和Gv為擴散（黏性）通量；S為N-S方程組在相對柱坐標下導(dǎo)出的源項；FB為無黏葉片力；FF為黏性葉片力。

定義由于葉片切向厚度產(chǎn)生的堵塞系數(shù)b

式中：Δφ=φs-φp，為葉片通道周向?qū)挾龋籒為葉片數(shù)。

在葉片區(qū)b＜1，在非葉片區(qū)b=1。周向脈動源項P中的高階項即為周向脈動應(yīng)力項，這些項是由于葉輪機的周向不均勻性及N-S方程的非線性所導(dǎo)致的，其與周向脈動源項均能夠反映周向不均勻性的大小。

與文獻[20]類似，為分析葉片通道進口徑向平衡的改變，對方程（1）中的徑向動量方程進行推導(dǎo)及簡化，可以得到無黏形式的完全徑向平衡方程

式中：等式左邊為徑向壓力梯度（RGP）；等式右邊分別為周向速度引發(fā)的離心加速度項（CENT_W）、由于子午流線的曲率而引發(fā)的離心加速度的徑向分量（CENT_M）、子午速度變化所引發(fā)的加速度的徑向分量（AC_M）、周向脈動源項的徑向分量 P'r，P'r與 Pr的區(qū)別在于二者的比值為密度，而下面將仍以Pr表示周向脈動源項的徑向分量，并不再重復(fù)說明；FBr為無黏葉片力的徑向分量，由于主要分析周向不均勻性對葉片通道進口流動平衡的影響，該處FBr=0，因此后面將不給出FBr的分布情況。

1.2 周向不均勻性的應(yīng)力輸運模型

類似雷諾平均后所產(chǎn)生的雷諾應(yīng)力項，式（1）中的周向脈動應(yīng)力項代表著周向不均勻性所帶來的影響，本節(jié)將簡要介紹周向脈動應(yīng)力項的應(yīng)力輸運模型，詳細建模過程可參考文獻[20]。類比湍流模型中的應(yīng)力輸運模型，可以推導(dǎo)出對于周向脈動應(yīng)力項的應(yīng)力輸運方程，以張量形式表示的無黏形式的應(yīng)力輸運方程為

式中：方程左側(cè)為對流項（CON）；右側(cè)分別為生成項（PRO）、輸運項（TRA）、速度 -壓力關(guān)聯(lián)項（PRE，簡稱壓力項）以及速度-源項（V-S）關(guān)聯(lián)項，其中生成項表征應(yīng)力輸運方程與平均后的軸對稱流場之間的關(guān)聯(lián)，速度-壓力關(guān)聯(lián)項表征由于壓力的脈動性而做的功，輸運項、壓力項及速度-源項關(guān)聯(lián)項均無法直接獲得，需要進行建模，而在這幾項中，壓力項為關(guān)鍵項，其余幾項則可忽略，因此下面簡要介紹壓力項的建模方法，詳細過程可參考文獻[20]。

假設(shè)速度和壓力可沿周向進行傅里葉級數(shù)展開，并保留到1階，則壓力項可進一步推導(dǎo)為

考慮到密度加權(quán)平均的定義，φ0和φ1之間滿足關(guān)系

從式中可見，當(dāng)確定3個速度和壓力的周向偏導(dǎo)數(shù)后，可解出壓力項。進一步引入S2流面無黏形式的周向動量方程以及絕對坐標系下流動無旋假設(shè)，并結(jié)合連續(xù)方程和無黏形式下的能量方程，可推導(dǎo)出壓力和3個速度分量的周向偏導(dǎo)數(shù)，其形式為

由此可實現(xiàn)對壓力項的求解，從而實現(xiàn)應(yīng)力輸運方程的封閉。另外從式（7）、（8）中可見，引發(fā)速度和壓力的周向不均勻性的1個主要因素是環(huán)量沿軸向和徑向的偏導(dǎo)數(shù)，即負荷的分配情況，而文獻[21]中也基于S2流面上周向動量方程分析了周向不均勻性產(chǎn)生的原因，認為對于定常問題周向不均勻性來源于黏性力的存在和氣流的環(huán)量沿流線的導(dǎo)數(shù)不為零，而在未發(fā)生大范圍分離時，相比黏性力，環(huán)量改變所反映的無黏葉片力是最主要的源項，因此環(huán)量改變對周向不均勻性的產(chǎn)生應(yīng)起主要作用。

環(huán)量的分布特征在設(shè)計問題中是給定量，在分析問題中則主要是源于無黏葉片力的作用，因此該模型同時適用于設(shè)計問題和分析問題。

2 應(yīng)力輸運模型的驗證及周向不均勻性的影響分析

2.1 研究對象及計算方法

采用1個前掠20°葉柵作為研究對象，該葉柵基元的厚度分布采用文獻 [22]中的設(shè)計進口馬赫數(shù)為0.62的可控擴散葉型，中弧線采用任意造型的方式獲得。葉柵的基本參數(shù)見表1。由于缺乏試驗數(shù)據(jù)作為驗證，本文采用3D數(shù)值模擬的結(jié)果作為基準，并應(yīng)用了NUMECA FINE/Turbo3D數(shù)值模擬軟件對各葉柵算例進行數(shù)值模擬研究。通過對某單級跨聲速風(fēng)扇、某風(fēng)扇/增壓級以及某單級離心壓氣機的數(shù)值計算結(jié)果與試驗結(jié)果的對比校核，NUMECA FINE/Turbo的計算精度能夠滿足壓氣機設(shè)計過程中的流場驗算和分析[23]。邊界條件的設(shè)置為：在計算域的進口處給定總溫、總壓以及氣流方向，通過改變進口氣流角來獲得不同的名義迎角，這里的名義迎角是指遠前方均勻來流氣流角與葉片幾何進口角的夾角，后面簡稱為迎角，給定的名義迎角為0°、±2°和±4°；采用當(dāng)?shù)赜堑母拍睿慈~片通道進口處氣流角與葉片幾何進口角的夾角。在出口給定流量以控制進口馬赫數(shù)，對于各個迎角，進口馬赫數(shù)均為0.62；環(huán)壁和葉片表面采用絕熱無滑移邊界。湍流模型采用Spalart-Allmaras（S-A）方程模型，網(wǎng)格劃分保證固壁面第1層網(wǎng)格單元y≤10。3D數(shù)值模擬計算網(wǎng)格采用AUTOGRID模塊自動生成，如圖1所示。對于不同迎角算例，均采用同一套網(wǎng)格，總網(wǎng)格數(shù)約為162萬。

表1 葉柵基本參數(shù)

圖1 3D數(shù)值模擬網(wǎng)格

另外，為分析周向不均勻性的影響，開展了以下2種類型的周向平均通流計算：

（1）CAM：未考慮周向不均勻性的周向平均通流計算；

（2）CAM+MODEL：加入了以應(yīng)力輸運模型計算的周向脈動應(yīng)力項的周向平均通流計算。

周向平均通流模型的計算網(wǎng)格如圖2所示，總網(wǎng)格數(shù)為3185。通流模型采用基于時間推進的有限體積方法來求解控制方程，時間離散格式為顯式格式為Jameson4步Runge-Kutta法。對無黏通量的空間離散采用Edwards[24]提出的低耗散通量分裂格式（LDFSS）。對于邊界條件的設(shè)置，進出口邊界條件與3D數(shù)值模擬相同，環(huán)壁采用對稱邊界條件。

圖2 通流計算網(wǎng)格

本文主要關(guān)注的是周向不均勻性對葉片通道進口流場的影響，此處黏性作用不明顯，雖然無黏葉片力在此處為0，但壓力勢已經(jīng)會以周向不均勻性的形式造成一定影響，下面將分別從葉片通道進口氣流角和徑向平衡的角度進行分析。

2.2 周向不均勻性對葉片通道進口流動參數(shù)的影響

為考察迎角改變時周向不均勻性對葉片通道進口流場的影響，首先分析在引入周向不均勻性后葉片通道進口氣流角的變化，本節(jié)將靠近葉片前緣處的氣流角與計算域進口氣流角做差以便于分析這一變化，結(jié)果如圖3所示。從圖3（a）中可見，對于20%展高處，隨著迎角的增大，3D和CAM預(yù)測的差值有一定的增長趨勢，而CAM+MODEL則呈現(xiàn)出下降趨勢，這主要是由于隨著迎角增大，Pu的作用增強（圖4（a）），導(dǎo)致周向速度變化更大，但總體來說CAM+MODEL所預(yù)測的差值CAM要更接近3D的，在名義迎角從-4°增大到4°的過程中，CAM+MODEL與CAM所預(yù)測的氣流角差值間的偏差量從約0.1°增加到了約0.7°，預(yù)測精度提高了25%以上；對于50%和80%展高處，在加入周向不均勻性的影響后，對氣流角的預(yù)測精度同樣提高，并且隨著名義迎角的增大，3D和CAM+MODEL預(yù)測出的αLE-αinlet的值均呈現(xiàn)出下降趨勢；與20%展高處的結(jié)果類似，加入周向不均勻性模型后通流模型對氣流角改變量的預(yù)測精度也同樣有所提高，提高50%以上。

圖3 掠葉柵不同展高處αLE-αinlet隨名義迎角的變化

由于葉片通道進口處無黏葉片力為0，因此上述氣流角的改變是由于周向不均勻性所誘發(fā)的，周向脈動源項的軸向和周向分量Px、Pu會分別影響葉片通道進口軸向和周向速度，從而造成氣流角的改變，并且隨著迎角的增大和負荷的提高，周向不均勻性的影響也在增大。

由于Pu會對環(huán)量造成影響，因此對其預(yù)測是否準確很重要。葉柵不同展高截面葉片上游距前緣約4%弦長處無量綱Pu隨迎角的變化情況如圖4所示，其中Pu是由計算域進口對應(yīng)流動參數(shù)進行的無量綱化。對于3個不同展高處，Pu的分布呈現(xiàn)出相近的特征，都是隨著迎角的增大Pu的絕對值增加，說明其作用將增強。此外，對于不同迎角，應(yīng)力輸運模型計算得到的Pu與3D結(jié)果的偏差在不同截面都很小，都在20%以內(nèi)，說明對于非設(shè)計工況，應(yīng)力輸運模型能夠很好地預(yù)測周向脈動源項。

在加入應(yīng)力輸運模型后，通流模型的計算時間比未加入應(yīng)力輸運模型時的約增加40%，這是由于在主控方程中加入6個微分方程，但計算時間仍遠遠少于3D數(shù)值模擬。各算法的計算時間見表2。

圖4 周向分量Pu隨名義迎角的變化掠葉柵不同展高處周向脈動源項

2.3 周向不均勻性對葉片通道進口流動徑向平衡的影響

上述分析表明，在加入周向不均勻性后，葉片通道進口氣流角發(fā)生了改變，即當(dāng)?shù)赜前l(fā)生了改變。對于掠葉柵，各展向截面當(dāng)?shù)赜堑淖兓瘯淖內(nèi)~片負荷沿展向的分布情況，進而改變進口徑向流動平衡。為進一步考察名義迎角和當(dāng)?shù)赜歉淖儗较蚱胶獾挠绊懀o出了不同名義迎角下葉柵上游距前緣約4%弦長處徑向平衡方程（3）中各項沿展向的分布情況（如圖5所示），并以流道進口的流動參數(shù)對這些項進行了無量綱化。圖中的命名規(guī)則以“-2_3D”為例，其中-2表示名義迎角為-2°，3D指3維計算的結(jié)果。

對比圖 5（a）～（e）中各項，從量上來看，與 RGP處于同一量級的項為CENT_M和Pr，而CENT_W和AC_M則要比其余3項小1～2個量級，表明周向不均勻性明顯影響葉片通道進口徑向平衡及徑向壓力梯度。

表2 葉柵基本參數(shù)

圖5 掠葉柵進口徑向平衡參數(shù)沿展向的分布隨名義迎角的變化

對于RGP隨名義迎角的變化情況（圖5（a）），對于 3D 結(jié)果，迎角從 -2°～0°增大時，RGP 有所減小，但并不明顯；而從0°～2°增大時，RGP的值則明顯減小。這說明名義迎角的增大會使徑向壓力梯度得到加強，而CAM的結(jié)果中RGP隨著名義迎角的增大則呈現(xiàn)出相反的趨勢，對徑向壓力梯度變化的描述不足。在加入周向不均勻性的影響之后，CAM+MODEL計算得到的各名義迎角下RGP的分布都更接近3D結(jié)果，隨迎角的變化趨勢也要更接近3D結(jié)果。結(jié)合圖5（b）～（e）可見，在加入周向不均勻性的影響后，RGP的變化主要來源于CENT_M和Pr，并且Pr對于RGP的變化起主要作用。

對于CENT_W（圖5（b）），3D、CAM+MODEL和CAM之間的偏差反映出周向分速度的變化，即在靠近葉片通道進口處環(huán)量已經(jīng)由于周向不均勻性的作用而發(fā)生改變。CAM+MODEL相比CAM偏離3D更多，說明加入周向不均勻性模型后通流計算得到的vu所受到的影響要大于3D數(shù)值模擬，但從圖中橫坐標的數(shù)值可見，CAM+MODEL與CAM和3D之間的實際區(qū)別并不大。而從趨勢上來看，CAM與3D之間確實存在偏差，而且在加入周向不均勻性模型后，CAM+MODEL的結(jié)果向著3D的方向移動，表現(xiàn)出了一定的一致性，說明在葉片通道進口vu確實會由于周向不均勻性的作用而發(fā)生改變；此外，當(dāng)迎角從-2°～2°增大時，CAM與3D之間的偏差增大，反映出了周向不均勻性的作用在增強，與此同時，CAM與CAM+MODEL之間的偏差也隨著迎角的變化呈現(xiàn)出了相似變化，說明加入周向不均勻性模型后通流模型能夠與3D數(shù)值模擬一樣反映出負荷變化對葉片通道進口流場的影響。

CENT_M（圖 5（c））在名義迎角從 -2°～2°增大后，其值呈現(xiàn)出減小趨勢，說明迎角增大后子午流線曲率增大。而CAM對CENT_M的預(yù)測要高于3D的結(jié)果，趨勢也相反；在引入周向不均勻性的影響后，CAM+MODEL所預(yù)測的CENT_M的精度有所提高，可達60%以上，趨勢也與3D一致，說明加入周向不均勻性影響后通流模型能夠更好地刻畫子午流線曲率。

AC_M（圖5（d））在名義迎角發(fā)生改變時變化并不是很明顯，其值略有增加。在引入周向不均勻性的影響后，CAM+MODEL的分布相比CAM更接近3D。

從圖5（e）中可見，隨著名義迎角的增大，周向脈動源項的徑向分量Pr的絕對值增大，反映出負荷沿展向有增大的趨勢。此外，應(yīng)力輸運模型計算得到的Pr的分布與3D結(jié)果基本一致，最大偏差在50%以內(nèi)。

綜上所述，CAM+MODEL能夠得到與3D相接近的分布，說明應(yīng)力輸運模型在迎角改變的情況下有能力反映出與3D計算相近的周向不均勻性的影響，提高通流計算的預(yù)測精度。同時，加入周向不均勻性后，通流模型能夠更好地反映出掠葉柵在不同迎角下流場中周向不均勻性的影響。

3 結(jié)論

研究表明，周向不均勻性會對葉輪機流場造成影響，因此有必要在通流模型中計入周向不均勻性的影響，以使其更好地計入真實3D流動效應(yīng)，提高設(shè)計、分析通流模型的預(yù)測精度。為?；⒎治鲇歉淖儠r掠葉柵進口流動周向不均勻性的影響，本文將1種?；芟蚱骄髂Ｐ椭兄芟蛎}動應(yīng)力項的應(yīng)力輸運模型結(jié)合到周向平均通流模型中，針對掠葉柵在不同迎角下的結(jié)果對比分析了該模型的可行性及周向不均勻性的作用，得出了以下結(jié)論：

（1）隨著迎角的增大，周向脈動源項的絕對值及其對葉片通道進口流場的影響都隨之提高，在本文算例中受周向不均勻性的影響葉片通道進口氣流角的改變量可達0.7°。

（2）在加入周向不均勻性模型后，在迎角改變的條件下通流計算預(yù)測的進口氣流角均更接近3D結(jié)果，對葉片通道進口氣流角的改變量的預(yù)測精度普遍提升30%以上，對葉片通道進口流動徑向平衡的描述精度提高60%以上。引入周向不均勻性有助于提升通流分析軟件對采用掠設(shè)計的壓氣機流場的計算精度。

（3）對于處于不同迎角下的葉柵流場，該應(yīng)力輸運模型都能較好的預(yù)測出周向脈動源項，在前緣前模型預(yù)測結(jié)果與3D結(jié)果的相對偏差普遍在20%之內(nèi)，最大不超過50%。

（4）由于應(yīng)力輸運模型在主控方程中增加了6個微分方程，因此對通流模型的計算時間影響較大，加入該模型后計算時間在原通流模型的基礎(chǔ)上增加了約40%，但仍遠遠少于3D數(shù)值模擬所用時間。該模型將周向不均勻性與環(huán)量分布關(guān)聯(lián)，可以為設(shè)計或分析問題提供周向脈動應(yīng)力項的預(yù)測。