楊金廣， 王春雪， 王大磊， 邵伏永， 楊晨， 吳虎

1.大連理工大學(xué) 能源與動(dòng)力學(xué)院， 大連 116024 2.中國(guó)航天科工集團(tuán) 北京動(dòng)力機(jī)械研究所， 北京 100074 3.西北工業(yè)大學(xué) 動(dòng)力與能源學(xué)院， 西安 710072

基于時(shí)間推進(jìn)的通流計(jì)算方法：現(xiàn)狀及展望

楊金廣1,*， 王春雪2， 王大磊2， 邵伏永2， 楊晨3， 吳虎3

1.大連理工大學(xué) 能源與動(dòng)力學(xué)院， 大連 116024 2.中國(guó)航天科工集團(tuán) 北京動(dòng)力機(jī)械研究所， 北京 100074 3.西北工業(yè)大學(xué) 動(dòng)力與能源學(xué)院， 西安 710072

為探究基于時(shí)間推進(jìn)的通流方法在葉輪機(jī)械設(shè)計(jì)分析中的應(yīng)用潛力，報(bào)告了國(guó)內(nèi)外相關(guān)研究現(xiàn)狀，詳細(xì)闡述了該方法在應(yīng)用過(guò)程中的若干關(guān)鍵問(wèn)題，包括：葉片堵塞的幾何描述、葉片力的模擬、葉片前后緣間斷問(wèn)題的處理和激波的捕獲等，并對(duì)不同模型進(jìn)行了比較?？偨Y(jié)可知，相比傳統(tǒng)通流方法，該方法具有捕捉激波，適應(yīng)堵塞，流場(chǎng)模擬更為精細(xì)以及計(jì)算擴(kuò)展性良好等諸多優(yōu)點(diǎn)，在先進(jìn)航空發(fā)動(dòng)機(jī)或燃?xì)廨啓C(jī)的部件設(shè)計(jì)和分析、整機(jī)穩(wěn)態(tài)、過(guò)渡態(tài)性能模擬中，都具有相當(dāng)?shù)膬?yōu)勢(shì)與應(yīng)用潛力；同時(shí)盡管該方法的研究取得了很大進(jìn)展，但在上述若干關(guān)鍵問(wèn)題上仍有提高空間。通過(guò)進(jìn)一步的發(fā)展與完善，該方法有望成為葉輪機(jī)械設(shè)計(jì)和分析的標(biāo)準(zhǔn)工具。

葉輪機(jī)械； 通流； 時(shí)間推進(jìn)； 葉片堵塞； 葉片力； 葉片前后緣間斷； 激波

吳仲華[1]于1952年提出三元流動(dòng)理論，將葉輪機(jī)械內(nèi)部全三維非定常黏性流動(dòng)近似為兩類(lèi)流面內(nèi)的準(zhǔn)三維定常流動(dòng)，即S1流面(葉片截面)和S2流面(子午面)，并在兩類(lèi)流面間迭代求解。三元流動(dòng)理論提出后被廣泛地應(yīng)用于葉輪機(jī)械的設(shè)計(jì)與性能預(yù)測(cè)中，極大地促進(jìn)了相關(guān)行業(yè)和學(xué)科的發(fā)展。而其中又以在S2流面上無(wú)葉和有葉區(qū)域內(nèi)都進(jìn)行計(jì)算的通流計(jì)算方法應(yīng)用最為廣泛，在風(fēng)扇、壓氣機(jī)以及渦輪等葉輪機(jī)械的設(shè)計(jì)體系中占據(jù)著重要地位。盡管全三維定常和非定常黏性計(jì)算迅速發(fā)展并且開(kāi)始應(yīng)用于葉輪機(jī)的設(shè)計(jì)中，發(fā)展先進(jìn)的通流計(jì)算方法仍然非常有必要。實(shí)際上，以20世紀(jì)90年代以來(lái)美、英等國(guó)的發(fā)動(dòng)機(jī)公司的風(fēng)扇、壓氣機(jī)設(shè)計(jì)體系為例，設(shè)計(jì)過(guò)程大致可以分為5個(gè)密切相關(guān)的步驟：初步設(shè)計(jì)、通流計(jì)算、葉片造型(二元?dú)鈩?dòng)優(yōu)化造型)、葉片造型(三元?dú)鈩?dòng)優(yōu)化造型)和放大尺寸的實(shí)驗(yàn)研究。這5個(gè)步驟環(huán)環(huán)相扣，每個(gè)階段采用不同層次的數(shù)學(xué)、物理模型和經(jīng)驗(yàn)數(shù)據(jù)，相互補(bǔ)充，相互交叉檢驗(yàn)，最終將設(shè)計(jì)風(fēng)險(xiǎn)降到最小[2]。相較于三維的黏性計(jì)算，通流計(jì)算所需要的時(shí)間要降低幾個(gè)數(shù)量級(jí)，尤其在多級(jí)計(jì)算中這一優(yōu)勢(shì)更為明顯。因此設(shè)計(jì)中可以將通流計(jì)算和三維計(jì)算作為互補(bǔ)的工具，后者可以在前者的計(jì)算結(jié)果基礎(chǔ)上進(jìn)一步驗(yàn)證；同時(shí)通流方法與三維計(jì)算相耦合，提供更加靈活的精度和計(jì)算周期選擇。此外通流方法在航空發(fā)動(dòng)機(jī)/燃?xì)廨啓C(jī)的整機(jī)模擬中具有很大的潛力，2015年P(guān)etrovic和Wiedermann[3]采用基于流函數(shù)的通流方法，結(jié)合簡(jiǎn)化的燃燒室模型，對(duì)整個(gè)燃?xì)廨啓C(jī)在S2流面上成功地進(jìn)行了性能模擬。

流行的通流計(jì)算方法有2種，Novak[4]提出的求解完全徑向平衡方程的流線曲率法，以及Marsh[5]提出的引入流函數(shù)的矩陣通流方法。這些方法具有快速穩(wěn)定的優(yōu)點(diǎn)，并且發(fā)展時(shí)間長(zhǎng)，理論較為完備，在工程實(shí)際中得到了廣泛的應(yīng)用。然而這2種方法均存在一個(gè)顯著的缺點(diǎn)，即由于激波的存在以及堵塞工況的出現(xiàn)，在跨、超聲速流動(dòng)中的應(yīng)用受到限制。因此在跨、超聲速風(fēng)扇、壓氣機(jī)和渦輪等葉輪機(jī)械設(shè)計(jì)中，需要一種實(shí)用的能夠處理高速流動(dòng)的通流方法，基于時(shí)間推進(jìn)求解Euler或Navier-Stokes方程的通流方法是就是一種非常合適的選擇。相比于流函數(shù)法或者流線曲率法，基于時(shí)間推進(jìn)的通流方法——計(jì)算流體力學(xué)(Computational Fluid Dynamics，CFD)通流方法具有以下優(yōu)點(diǎn)：① 質(zhì)量流量由流場(chǎng)計(jì)算本身得出，因而可以預(yù)測(cè)葉片排的堵塞工況；② 可以自動(dòng)捕獲特定形式的激波，從而進(jìn)行跨、超聲速流場(chǎng)計(jì)算；③ 可以應(yīng)用CFD領(lǐng)域發(fā)展成果，如有限體積方法、有限元方法以及高精度數(shù)值格式等，同時(shí)時(shí)間推進(jìn)的引入使得CFD通流方法可以自然地應(yīng)用于過(guò)渡態(tài)的計(jì)算。

本文對(duì)國(guó)內(nèi)外基于時(shí)間推進(jìn)的通流計(jì)算方法進(jìn)行總結(jié)，對(duì)該方法在具體實(shí)施中的關(guān)鍵問(wèn)題處理方法進(jìn)行討論，最后對(duì)其未來(lái)的發(fā)展方向進(jìn)行展望。

1 國(guó)內(nèi)外研究進(jìn)展

在過(guò)去的幾十年里，基于時(shí)間推進(jìn)的通流計(jì)算方法有了一定程度的發(fā)展。對(duì)于國(guó)外而言，Spurr[6]最早于1980年在通流計(jì)算中采用時(shí)間推進(jìn)的方法對(duì)周向平均的Euler方程進(jìn)行了求解。英國(guó)劍橋大學(xué)Whittle實(shí)驗(yàn)室的Dawes[7]對(duì)三維Navier-Stokes方程求解器進(jìn)行了修改，使其能夠進(jìn)行Euler通流模型的求解，并以多級(jí)壓氣機(jī)為例，使用CFD通流方法提供多級(jí)環(huán)境，僅對(duì)單排葉片進(jìn)行三維計(jì)算，對(duì)CFD通流方法與三維黏性數(shù)值計(jì)算的結(jié)合進(jìn)行了探索性研究。比利時(shí)布魯塞爾自由大學(xué)的Yao和Hirsch[8]基于三維求解器修改實(shí)現(xiàn)了CFD通流方法計(jì)算，并細(xì)致分析了葉片區(qū)弦向氣流角和堵塞分布的影響。美國(guó)雪城大學(xué)的Damle等[9-10]發(fā)展了求解Euler方程的完全通流求解器，并實(shí)現(xiàn)了對(duì)跨聲速和超聲速葉片的設(shè)計(jì)，展示了其在跨、超聲速設(shè)計(jì)問(wèn)題中的潛力。瑞典查爾姆斯理工大學(xué)的Baralon等[11-13]對(duì)CFD通流方法進(jìn)行了深入的研究，提出了一種基于偽時(shí)間步長(zhǎng)的葉片力計(jì)算模型，并針對(duì)葉片前緣的間斷問(wèn)題提出了2種處理方法，使得CFD通流方法的求解穩(wěn)定性和實(shí)用性得到了提升；Baralon等同時(shí)對(duì)Euler通流模型的激波捕獲機(jī)理進(jìn)行了研究，提出了一種計(jì)算流面法向堵塞因子的新方法，并在Euler通流模型中引入了葉型損失、落后角和端壁損失，此外還引入了徑向摻混模型。Sturmayr和Hirsch[14]對(duì)通道激波的捕獲機(jī)理進(jìn)行了研究，并對(duì)設(shè)計(jì)問(wèn)題和分析問(wèn)題進(jìn)行了對(duì)比。進(jìn)入21世紀(jì)，CFD通流方法研究進(jìn)入高潮，比利時(shí)列日大學(xué)的Simon等[15-17]發(fā)展了求解周向平均Navier-Stokes方程的通流求解器，采用2個(gè)亞聲速算例進(jìn)行了驗(yàn)證，并與Euler通流模型進(jìn)行了對(duì)比，此外還對(duì)葉片力項(xiàng)、周向應(yīng)力、確定應(yīng)力等在周向平均通流模型中的相對(duì)重要性進(jìn)行了詳細(xì)研究。意大利米蘭理工大學(xué)的Persico和布雷西亞大學(xué)的Rebay[18]共同提出了一種基于罰函數(shù)的葉片力模型。都靈理工的Taddei和Larocca[19]對(duì)CFD通流方法的設(shè)計(jì)問(wèn)題進(jìn)行了持續(xù)的研究，并且提出了2種處理葉片前后緣間斷問(wèn)題的新方法[20-21]。佛羅倫薩大學(xué)的Pacciani等[22]提出了一種自適應(yīng)流面修改方法，以解決前后緣間斷問(wèn)題。值得注意的是俄羅斯中央航空發(fā)動(dòng)機(jī)研究院(Central Institute of Aviation Motors，CIAM)的Nigmatullin和Ivanov[23]基于高階精度的Godunov格式，采用貼體坐標(biāo)系，完全獨(dú)立地建立了一整套CFD通流計(jì)算理論模型，包括葉片力的模擬和葉片前后緣間斷問(wèn)題處理模型等，并且已經(jīng)成功應(yīng)用于航空發(fā)動(dòng)機(jī)整機(jī)流場(chǎng)計(jì)算。

相比于國(guó)外CFD通流方法的長(zhǎng)期發(fā)展，國(guó)內(nèi)相關(guān)研究起步較晚。袁寧等[24]對(duì)流函數(shù)法、流線曲率法以及求解Euler方程的CFD通流方法3種S2流面計(jì)算程序進(jìn)行了比較。季路成等[25]在國(guó)內(nèi)較早地將時(shí)間推進(jìn)思想應(yīng)用于通流計(jì)算，并設(shè)計(jì)了一軸向流動(dòng)跨聲速的激波轉(zhuǎn)子進(jìn)行驗(yàn)證。施發(fā)樹(shù)和劉興洲[26-27]與CIAM聯(lián)合開(kāi)發(fā)了渦輪發(fā)動(dòng)機(jī)整機(jī)通流計(jì)算程序，并對(duì)一臺(tái)小型雙涵道渦扇發(fā)動(dòng)機(jī)進(jìn)行了整機(jī)S2計(jì)算。清華的于龍江[28]、北航的曹志鵬[29]和金東海[30]等使用CFD通流方法進(jìn)行了整機(jī)流場(chǎng)仿真計(jì)算。李德英等[31]基于CFD通流方法開(kāi)展了多級(jí)渦輪的氣動(dòng)分析，并得到了與三維計(jì)算較為一致的結(jié)果。萬(wàn)科等[32-33]基于周向平均的Navier-Stokes方程通流計(jì)算方法考察了周向平均方法在風(fēng)扇/增壓級(jí)分析中的準(zhǔn)確性，并對(duì)周向不均勻項(xiàng)對(duì)于周向平均模型的影響進(jìn)行了研究。

盡管CFD通流方法已有幾十年的發(fā)展時(shí)間，理論上相比傳統(tǒng)的通流方法也更具優(yōu)勢(shì)，但要在工程實(shí)踐中成為一種成熟的設(shè)計(jì)手段仍然面臨諸多困難。本文將就各國(guó)研究者在CFD通流方法研究方面所做的工作進(jìn)行介紹，同時(shí)針對(duì)一些關(guān)鍵問(wèn)題展開(kāi)討論，分析影響計(jì)算準(zhǔn)確性的關(guān)鍵問(wèn)題，提出可能的解決方案，并展望未來(lái)的發(fā)展趨勢(shì)。

2 CFD通流方法基本計(jì)算模型

通流計(jì)算的基本思想在于認(rèn)為葉輪機(jī)中的流動(dòng)在無(wú)葉區(qū)滿足軸對(duì)稱假設(shè)，在葉片通道內(nèi)則沿S2流面流動(dòng)。因此，CFD通流計(jì)算模型一般通過(guò)對(duì)三維的Euler/Navier-Stokes方程在相應(yīng)方向進(jìn)行平均得到，從而實(shí)現(xiàn)從全三維的流動(dòng)問(wèn)題(圖1(a))向平均葉片S2流面求解(圖1(b))的降維簡(jiǎn)化[18]。如果此時(shí)進(jìn)一步在流場(chǎng)中施加適當(dāng)形式的體積力來(lái)模擬葉片的作用，迫使當(dāng)?shù)亓鲃?dòng)速度方向滿足與平均葉片流面相切，則問(wèn)題可以進(jìn)一步簡(jiǎn)化到子午平面內(nèi)的流場(chǎng)計(jì)算問(wèn)題(圖1(c))。

圖1 從全三維模型向通流模型的簡(jiǎn)化[18]Fig.1 Reduction from fully 3D to throughflow models[18]

不同的葉片力模型會(huì)決定計(jì)算模型是否需要對(duì)葉片區(qū)周向動(dòng)量方程進(jìn)行求解。因此，根據(jù)物理計(jì)算域的簡(jiǎn)化程度以及是否在有葉區(qū)求解周向動(dòng)量方程，目前的CFD通流方法可以分為3類(lèi)：

第1類(lèi)，在葉片區(qū)不需求解周向動(dòng)量方程(即下文式(1)中第4個(gè)方程)，周向動(dòng)量方程用來(lái)推導(dǎo)體積力形式的周向葉片力。該類(lèi)計(jì)算模型相當(dāng)于在葉片區(qū)求解了二維的流場(chǎng)問(wèn)題。該方法可以同時(shí)進(jìn)行分析問(wèn)題(正問(wèn)題)和設(shè)計(jì)問(wèn)題(反問(wèn)題)的求解。求解正問(wèn)題時(shí)，在葉片區(qū)每一時(shí)間步中通過(guò)子午速度計(jì)算滿足流面相切條件的周向速度，并對(duì)葉片力進(jìn)行更新，計(jì)算得到新的子午速度，直到得到流場(chǎng)的收斂解。在反問(wèn)題中，則首先規(guī)定載荷或者環(huán)量分布，通過(guò)流面相切條件計(jì)算得到新的平均葉片流面幾何；作為一個(gè)初值問(wèn)題，流面幾何的時(shí)間推進(jìn)計(jì)算需要設(shè)定邊界條件，通常通過(guò)確定前緣傾斜形狀或者積疊規(guī)律實(shí)現(xiàn)；根據(jù)新的流面幾何用來(lái)更新葉片力，直到得到收斂解。Dawes[7]、Damle[9]和Taddei[20]等均采用了這種計(jì)算模型，所不同的是Dawes的通流求解器由三維Navier-Stokes方程求解器修改而來(lái)，通過(guò)在周向僅布置一個(gè)網(wǎng)格單元結(jié)合周期邊界條件實(shí)現(xiàn)CFD通流方法求解；同時(shí)Dawes只進(jìn)行了通流正問(wèn)題的求解，主要用于為單葉片排三維計(jì)算提供多級(jí)環(huán)境，而Damle和Taddei均同時(shí)進(jìn)行了正、反問(wèn)題的求解。

第2類(lèi)，在葉片區(qū)求解周向動(dòng)量方程，通過(guò)添加適當(dāng)形式的體積力，模擬無(wú)黏葉片力的強(qiáng)迫作用，使得在最終收斂的流場(chǎng)計(jì)算結(jié)果中，流體流動(dòng)方向滿足流面相切條件。該種計(jì)算模型在子午平面內(nèi)進(jìn)行求解，不同于第1類(lèi)模型。在正問(wèn)題求解中，該類(lèi)模型的周向速度完全通過(guò)求解得出，可以認(rèn)為具有準(zhǔn)三維的特征；同時(shí)葉片力模型的數(shù)學(xué)性質(zhì)更有利于計(jì)算迭代的穩(wěn)定和收斂，因此得到了廣泛使用。Sturmayr[14]、Baralon[11]、Simon[15]及Persico[18]等均采用了此種CFD通流方法基本計(jì)算模型。其中Sturmayr和Hirsch并沒(méi)有求解真實(shí)的葉型幾何，而是通過(guò)冪函數(shù)的形式模擬了中弧線以及堵塞因子的分布，使其符合真實(shí)葉型的主要特征，以此探索激波捕捉機(jī)理以及不同的參數(shù)分布對(duì)正問(wèn)題求解的影響。Baralon[11]和Simon[34]等均嘗試了采用法向堵塞因子，希望借此提高流場(chǎng)的計(jì)算精度。由于葉片力模型的不同，該類(lèi)方法在反問(wèn)題中的使用與第1類(lèi)模型有所不同。Sturmayr和Hirsch[14]對(duì)正問(wèn)題和反問(wèn)題進(jìn)行了比較，在反問(wèn)題中首先指定葉片尾緣的周向速度，然后根據(jù)來(lái)流周向速度，在葉片內(nèi)從前緣到尾緣按冪函數(shù)形式分布作為目標(biāo)周向速度，當(dāng)計(jì)算收斂時(shí)，計(jì)算得到的周向速度與目標(biāo)周向速度趨同，此時(shí)的流面即設(shè)計(jì)目標(biāo)。Simon[34]采用了相同的方法進(jìn)行反問(wèn)題的求解。

前兩類(lèi)CFD通流方法均在子午平面內(nèi)進(jìn)行流場(chǎng)計(jì)算。以周向平均后的軸對(duì)稱Euler方程為例，一般的控制方程表達(dá)形式為

(1)

式中：

其中：U為守恒量；F和G為對(duì)流(無(wú)黏)通量；S為源項(xiàng)；ρ為密度；p為壓力；z、r和φ分別為軸向、徑向和周向坐標(biāo)；u、v和w分別為軸向、徑向和周向速度分量；e為總能；ε(T)為溫度T時(shí)的單位質(zhì)量?jī)?nèi)能；b(r,z)為堵塞因子，模擬葉片厚度帶來(lái)的堵塞效應(yīng)；FB為無(wú)黏葉片力項(xiàng)；fBz、fBr和fBφ分別為FB在3個(gè)方向的分量；FL為黏性力項(xiàng)；fLz、fLr和fLφ分別為FL在3個(gè)方向的分量；在葉片區(qū)以外FB為0；t為時(shí)間；ω為角速度。

第3類(lèi)，直接在三維平均葉片流面上求解。目前的公開(kāi)文獻(xiàn)中，僅有Nigmatullin和Ivanov[23]通過(guò)對(duì)相對(duì)圓柱坐標(biāo)系下的三維控制方程進(jìn)行空間坐標(biāo)平均，建立了基于貼體坐標(biāo)系(ξ,η,ζ)的CFD通流方法基本計(jì)算模型。該計(jì)算模型在葉片區(qū)和無(wú)葉區(qū)采用了不同的控制方程，在無(wú)葉區(qū)滿足軸對(duì)稱條件，在葉片區(qū)對(duì)每個(gè)流面網(wǎng)格單元求解流面法向平均的控制方程。相比于前兩類(lèi)計(jì)算模型在子午平面進(jìn)行通量的求解，該計(jì)算模型通量的求解無(wú)疑帶有更多的三維特征，同時(shí)該模型的葉片力模擬也更為方便，通過(guò)理論推導(dǎo)即可得出。該葉片力模型能夠與流場(chǎng)的迭代同步進(jìn)行。該模型對(duì)于葉片前后緣間斷問(wèn)題采用了任意間斷分解方法進(jìn)行處理，較為完整地考慮了流動(dòng)中可能出現(xiàn)的各種情況，同時(shí)魯棒性較好，是非常具有應(yīng)用潛力的CFD通流方法。

3 CFD通流方法中的若干關(guān)鍵問(wèn)題

CFD通流方法在應(yīng)用中面臨諸多問(wèn)題，如葉片堵塞的描述、葉片力的模擬、葉片前后緣間斷問(wèn)題的處理以及激波的捕獲等，只有處理好這些問(wèn)題，CFD通流方法才可能作為一種快速且準(zhǔn)確的設(shè)計(jì)分析手段應(yīng)用于葉輪機(jī)械領(lǐng)域。下文將對(duì)這些問(wèn)題以及目前常見(jiàn)的處理方法進(jìn)行詳述。

3.1 葉片堵塞的幾何描述

通過(guò)對(duì)三維控制方程進(jìn)行平均，CFD通流方法實(shí)現(xiàn)了流場(chǎng)計(jì)算的降維，但與此同時(shí)，在葉片區(qū)需要使用堵塞因子對(duì)真實(shí)葉片厚度造成的流道堵塞效應(yīng)進(jìn)行模擬。堵塞因子的定義一般分為2種，一種為周向堵塞因子，另一種為法向堵塞因子。周向堵塞因子通過(guò)對(duì)三維Navier-Stokes方程進(jìn)行周向平均后直接導(dǎo)出，定義為

(2)

式中：φp-φs為葉片在周向以弧度計(jì)量的厚度；Nb為葉片數(shù)；π為圓周率。顯然在無(wú)葉區(qū)堵塞因子為1。Sturmayr和Hirsch[14]采用了周向堵塞因子進(jìn)行計(jì)算，但沒(méi)有使用真實(shí)的葉型幾何，而是近似計(jì)算葉型法向厚度最大點(diǎn)對(duì)應(yīng)的周向厚度，即

dmax≈dnmax/cosβnmax

(3)

式中：dnmax為葉型法向最大厚度；βnmax為法向厚度最大位置點(diǎn)對(duì)應(yīng)的葉片角。以dmax為基準(zhǔn)，沿流向按冪函數(shù)分布確定葉片周向厚度，并在葉片前后緣附近對(duì)堵塞因子分布進(jìn)行光滑。通過(guò)試取不同的冪指數(shù)和調(diào)整dmax的弦向位置實(shí)現(xiàn)不同的堵塞因子分布?；鶞?zhǔn)算例冪指數(shù)為2，dmax位于0.5cc處。他們研究了流場(chǎng)計(jì)算結(jié)果對(duì)于堵塞因子分布的敏感度，如表1[14],其中:c為弦長(zhǎng);Case 1 為基準(zhǔn)算例。由表1的計(jì)算結(jié)果發(fā)現(xiàn)，不同堵塞因子分布對(duì)堵塞流量以及對(duì)應(yīng)效率的計(jì)算結(jié)果有明顯影響。

為了更好地估計(jì)實(shí)際葉片厚度對(duì)流動(dòng)的影響，Baralon等[11]認(rèn)為采用法向堵塞因子比周向堵塞因子更為合適，并給出了相應(yīng)的法向堵塞因子計(jì)算方法，表達(dá)式為

(4)

表1 堵塞流量對(duì)堵塞因子分布的敏感度[14]

Table 1 Sensitivity of blocking mass flow to blockage

factor distribution[14]

CaseModificationMassflow/(kg·s-1)PressureratioEfficiency/%1dmaxat0.5c30.81.53782dmaxat0.75c32.71.52863dmaxat0.25c27.21.53664Nosmoothing28.91.53725dmax-20%32.21.54826dmax+20%29.31.51757dmaxexponent1.531.81.53818dmaxexponent329.11.5275

式中：S為柵距；β為平均流面角；β1和β2為葉片角；t1和t2為相應(yīng)的葉片厚度。該表達(dá)式僅限于計(jì)算點(diǎn)A和點(diǎn)B之間的堵塞因子，如圖2[11]所示，從上方葉片前緣點(diǎn)向平均流面(在S1面上為平均流線)作垂線，垂線與平均流面交點(diǎn)即為A點(diǎn)；同理，點(diǎn)B由下方葉片尾緣點(diǎn)向平均流面作垂線得到；點(diǎn)P為點(diǎn)A和點(diǎn)B之間的平均流線上任一點(diǎn)。對(duì)于前緣到點(diǎn)A和尾緣到點(diǎn)B的部分，堵塞因子定義為從bnA和bnB線性變化到1?？梢园l(fā)現(xiàn)，該方法定義的法向堵塞因子與周向堵塞因子差別較為明顯，對(duì)于某一展向葉片截面，計(jì)算得到的流道面積大小和喉道軸向位置均不相同，這些差異會(huì)對(duì)流場(chǎng)計(jì)算結(jié)果造成較為明顯的影響，尤其在激波存在的跨、超聲速流動(dòng)中。該種堵塞因子定義方法存在實(shí)際操作復(fù)雜、前后緣部分分布過(guò)于簡(jiǎn)單、精確度較差的缺點(diǎn)。

Simon[34]借鑒Baralon的思路，提出了一種簡(jiǎn)化的法向堵塞因子計(jì)算模型

bn=bcosβ

(5)

式中：b為周向堵塞因子，cosβ為對(duì)應(yīng)葉片角余弦值。Simon使用Rotor67轉(zhuǎn)子作為算例，采用2種堵塞因子分別進(jìn)行了計(jì)算，并對(duì)得到的結(jié)果進(jìn)行了比較，發(fā)現(xiàn)采用法向堵塞因子的堵點(diǎn)流量計(jì)算結(jié)果與實(shí)驗(yàn)值符合得更好，并且激波更為靠前。

圖2 法向堵塞因子[11]Fig.2 Normal blockage factor[11]

Nigmatullin和Ivanov[23]的計(jì)算模型將周向堵塞因子與貼體坐標(biāo)系和圓柱坐標(biāo)系之間的變換進(jìn)行了關(guān)聯(lián)，定義b=1/ζφ，應(yīng)用于度量系數(shù)的計(jì)算當(dāng)中，從而模擬葉片堵塞效應(yīng)?？紤]到ζ=const對(duì)應(yīng)于平均S2流面，因此該定義可以作為法向堵塞因子的一種標(biāo)準(zhǔn)定義。

從流動(dòng)的物理特征上理解，基于周向平均的堵塞因子定義不能真切反應(yīng)葉片真實(shí)存在造成的物理堵塞效應(yīng)，而在這方面法向堵塞因子更具優(yōu)勢(shì)，但其計(jì)算方法仍有待改進(jìn)。另外堵塞因子的定義要與控制方程相匹配，如果采用周向平均堵塞因子，則可以采用周向堵塞因子；如果沿平均流面法向進(jìn)行平均，則采用法向堵塞因子更好，這也是很多文獻(xiàn)中忽略的問(wèn)題。總之，目前對(duì)于堵塞因子的研究仍然不夠系統(tǒng)和全面，堵塞因子的定義以及堵塞效應(yīng)的數(shù)學(xué)描述需要從理論上更進(jìn)一步完善。

3.2 葉片力

通流計(jì)算中，除了需要考慮實(shí)際葉片厚度帶來(lái)的影響，另一個(gè)主要問(wèn)題便是如何在葉片區(qū)模擬葉片對(duì)流體的作用力。Marble[35]最早提出采用分布的體積力模擬葉輪機(jī)流場(chǎng)，這種理念也被拓展應(yīng)用于通流計(jì)算中，用分布的體積力模擬葉片力[36-37]。通常將葉片力分為無(wú)黏葉片力和黏性力兩部分進(jìn)行模擬，從而達(dá)到簡(jiǎn)化的目的。對(duì)于無(wú)黏葉片力，其作用方向通常定義為葉片平均流面法向，即

FBW=uρfBz+vρfBr+wρfBφ=0

(6)

式中：FB為無(wú)黏葉片力；W為相對(duì)速度矢量；定義平均流面α=φ-g(r,z)，g為r和z的函數(shù)，流面單位法向量

單位法向量

則FB可表示為

FB=B(r,z)·ng

(7)

葉片力的軸向和徑向分量為

由此FB的方向已經(jīng)確定，B(r,z)為葉片力。

一般存在2種思路去模擬無(wú)黏葉片力，一種通過(guò)周向動(dòng)量方程直接計(jì)算得到無(wú)黏葉片力周向分量ρfBφ，如式(8)所示，通過(guò)投影即可得到B(r,z)及各分量的大小，Dawes[7]和Damle[9]等均采用了此種方法。

(8)

式中：m為子午方向；qm為質(zhì)量流量；Vφ為周向速度。

另一種思路為將葉片力看做時(shí)間相關(guān)的未知量，通過(guò)求解時(shí)間相關(guān)的微分方程對(duì)無(wú)黏葉片力進(jìn)行計(jì)算，相比于求解式(6)的剛性約束，該方法數(shù)學(xué)性質(zhì)上更易于穩(wěn)定求解。該方法基本思想在于周向葉片力的表達(dá)必須滿足流面相切條件，基于此Baralon等[11]提出的一種采用偽時(shí)間步長(zhǎng)的葉片力計(jì)算模型，時(shí)間相關(guān)的微分方程定義為

(9)

式中：λ為典型響應(yīng)時(shí)間平方的倒數(shù)，s-2。ρfBφ的響應(yīng)時(shí)間應(yīng)該和流場(chǎng)計(jì)算中的當(dāng)?shù)鼐W(wǎng)格時(shí)間步長(zhǎng)處于同一量級(jí)，因此λ可以表達(dá)為

式中：Δx為典型網(wǎng)格尺寸；Δt為當(dāng)?shù)鼐W(wǎng)格時(shí)間步長(zhǎng)。該方程可以理解為一種壓力梯度，一旦流面相切條件不符合就會(huì)產(chǎn)生，該方法相當(dāng)于對(duì)周向葉片力在時(shí)間上的一種離散。Simon和Léonard[15]采用了類(lèi)似的葉片力模擬方法，只不過(guò)采用葉片力的模作為求解變量。萬(wàn)科等[33]指出該種葉片力模型在實(shí)際使用過(guò)程中，λ的取值對(duì)迭代收斂的影響較大，取值不當(dāng)會(huì)影響收斂速度，嚴(yán)重時(shí)甚至?xí)?dǎo)致迭代不收斂。

(10)

式中：C為松弛因子。

Persico和Rebay[18]同樣基于流面相切條件，提出了一種基于罰函數(shù)的葉片力模型為

B(r,z)=-KW·ng

(11)

式中：K為懲罰因子。該模型可以理解為將流面相切條件施加于有彈性的葉片上，葉片的剛度即懲罰因子K。在實(shí)踐中發(fā)現(xiàn)當(dāng)K≥104時(shí)，最終的計(jì)算結(jié)果中流面相切條件的偏差不大于0.1°。同時(shí)采用上述模型的好處在于不需要求解額外的微分方程。此外可以在計(jì)算進(jìn)程中調(diào)節(jié)剛度K，在開(kāi)始迭代階段K取小值，從而容許流面相切條件有較大的偏差，防止求解的發(fā)散，隨著計(jì)算的進(jìn)行，剛度K不斷增大，從而收斂到滿意的解。該模型的根本目的在于將式(6)的剛性約束轉(zhuǎn)化為柔性約束，利于求解。

不同于上述2種思路，Nigmatullin和Ivanov[23]采用了另外一種顯式的葉片力計(jì)算方法。他們首先計(jì)算出軸向、徑向和周向動(dòng)量方程中不考慮無(wú)黏葉片力項(xiàng)的殘差，對(duì)貼體坐標(biāo)系下離散后的動(dòng)量方程進(jìn)行理論推導(dǎo)，并利用流面相切條件，可將無(wú)黏葉片力通過(guò)計(jì)算出的殘差項(xiàng)進(jìn)行表達(dá)，這樣隨著流場(chǎng)計(jì)算的收斂，葉片力的分布也一同收斂。葉片力具體表達(dá)式為

(12)

式中：J為坐標(biāo)系變換的雅克比行列式；ζz、ζr和ζφ為與ζ=const網(wǎng)格面對(duì)應(yīng)的軸向、徑向和周向的度量系數(shù)；S2、S3和S4分別為軸向、徑向和周向動(dòng)量方程中不考慮無(wú)黏葉片力項(xiàng)的殘差。該方法簡(jiǎn)單通用，不需要引入偽時(shí)間項(xiàng)，亦避免了因此帶來(lái)的常系數(shù)取值問(wèn)題，同時(shí)具有較好的收斂性。

對(duì)于黏性力，一般可以寫(xiě)成如式(13)所示的表達(dá)式。

(13)

式中：黏性力方向與當(dāng)?shù)厮俣确较蛳喾?。L(r,z)為一正值，代表黏性損失的大小。對(duì)于黏性損失的計(jì)算，Bosman和Marsh[38]提出通過(guò)給定流場(chǎng)的熵分布進(jìn)行計(jì)算，該損失模型簡(jiǎn)單且使用廣泛。由此L(r,z)可以和熵分布聯(lián)系起來(lái)，即

(14)

式中：T為當(dāng)?shù)仂o溫；s為熵。

3.3 葉片前、后緣間斷問(wèn)題

葉片前緣沖角以及尾緣落后角在CFD通流計(jì)算中是需要被重點(diǎn)關(guān)注的問(wèn)題。由于采用了軸對(duì)稱假設(shè)，CFD通流計(jì)算方法無(wú)法在無(wú)葉區(qū)感受到葉片存在而提前產(chǎn)生旋流速度；而一旦進(jìn)入葉片區(qū)，流動(dòng)則必須滿足流面相切條件，周向速度會(huì)被平均葉片流面角度所約束，因此在葉片前緣網(wǎng)格界面會(huì)產(chǎn)生速度的間斷。如果采用葉片力源項(xiàng)模擬前緣的速度轉(zhuǎn)折，就會(huì)產(chǎn)生體積力的突變，即便不斷加密網(wǎng)格也無(wú)濟(jì)于事。在實(shí)踐中發(fā)現(xiàn)，只要前緣來(lái)流沖角大于2°，就會(huì)產(chǎn)生比較明顯的數(shù)值損失[11]，而這種損失所對(duì)應(yīng)的熵增是非物理的。對(duì)于葉片尾緣，同樣會(huì)產(chǎn)生類(lèi)似的間斷問(wèn)題。

上述問(wèn)題在很長(zhǎng)一段時(shí)間內(nèi)沒(méi)有得到足夠的重視，隨著CFD通流方法的發(fā)展，Baralon等[11]較早關(guān)注和深入研究了該問(wèn)題，并給出了2種處理辦法。第1種通過(guò)在流場(chǎng)計(jì)算過(guò)程中人為修改葉片前緣部分流面，使平均葉片流面角度與來(lái)流角度相吻合，沖角變?yōu)榱悖@種方法相當(dāng)于引入了長(zhǎng)度尺度，從而消除了間斷問(wèn)題。對(duì)于尾緣可以采用類(lèi)似的方法，根據(jù)輸入的落后角大小對(duì)尾緣部分流面進(jìn)行修改。這種方法在過(guò)去的十幾年得到了較為廣泛的應(yīng)用，包括Simon[34]和Persico[18]等學(xué)者均采用了此種前后緣處理方法。圖3給出了修改葉片前緣流面和不進(jìn)行處理計(jì)算所產(chǎn)生的熵增對(duì)比[34]，其中l(wèi).e.(leading edge)指葉片前緣。

盡管簡(jiǎn)單明了且應(yīng)用廣泛，此種方法仍存在明顯的缺點(diǎn)：首先，在沖角較大時(shí)，對(duì)平均葉片流面的修改程度較大，對(duì)于流場(chǎng)計(jì)算結(jié)果的影響會(huì)比較大；其次，在超聲速流場(chǎng)計(jì)算時(shí)對(duì)激波的捕獲會(huì)產(chǎn)生明顯影響。由于前緣流面的修改會(huì)導(dǎo)致葉片流道面積的變化，對(duì)激波是否產(chǎn)生以及激波位置均有決定性影響，進(jìn)一步對(duì)流場(chǎng)的速度分布以及損失預(yù)估的精確度產(chǎn)生重要影響。為此，Baralon等[11]同時(shí)給出了第2種方法：動(dòng)量間斷處理。這種方法通過(guò)在前后緣兩側(cè)建立質(zhì)量守恒、能量守恒關(guān)系，同時(shí)在動(dòng)量方程中引入另外的目標(biāo)間斷力項(xiàng)。所謂“間斷力”用來(lái)控制前緣下游氣流角的間斷以及所產(chǎn)生的熵增。該方法主要分為2個(gè)求解步驟：第1步求解目標(biāo)間斷力，通過(guò)引入流面相切、前緣展向連續(xù)以及給定前緣總壓損失3個(gè)條件，與上述守恒關(guān)系聯(lián)立求解得到；第2步即在間斷力已知的情況下求解兩側(cè)流動(dòng)參數(shù)，并對(duì)前緣網(wǎng)格界面通量進(jìn)行修改。該方法通過(guò)間斷關(guān)系式對(duì)前、后緣問(wèn)題進(jìn)行了解決，但在實(shí)際應(yīng)用中仍然存在問(wèn)題，即流動(dòng)參數(shù)的迭代求解在數(shù)學(xué)上存在困難，尤其在馬赫數(shù)接近1時(shí)，因此無(wú)法應(yīng)用于跨聲速流動(dòng)的情況，如跨聲速轉(zhuǎn)子中。此外該方法還需要人為輸入額外的參數(shù)，即總壓損失。

基于對(duì)前后緣流面進(jìn)行修改的思路，Taddei和Larocca[20]提出了采用反問(wèn)題修正前緣流面的方法。相比于人為地直接對(duì)平均葉片流面進(jìn)行修改，該方法通過(guò)給定來(lái)流到前緣下游部分的旋流速度分布，求解出新的前緣部分流面。由于直接對(duì)速度分布進(jìn)行操作，該方法對(duì)于前緣流場(chǎng)的控制更好，彌補(bǔ)了Baralon第1種方法的一些缺點(diǎn)，但仍然帶有一定的任意性，并且在來(lái)流沖角較大且葉片扭曲程度較強(qiáng)時(shí)容易計(jì)算發(fā)散。為此Taddei和Larocca[21]提出了一種基于激盤(pán)AD

圖3 葉片前緣的修改[34]Fig.3 Adaptation of leading edge [34]

(Actuator Disk)模型的前后緣間斷問(wèn)題處理方法。通過(guò)在葉片前緣或尾緣采用如圖4[21]所示的激盤(pán)模型，使得前緣或尾緣網(wǎng)格線i+1/2兩側(cè)保證質(zhì)量守恒、能量守恒和熵守恒，同時(shí)集成了沖角或者落后角關(guān)系。其中：x為空間軸;t為時(shí)間軸;Vx為x方向速度；a為當(dāng)?shù)芈曀?i為網(wǎng)格標(biāo)號(hào)。通過(guò)對(duì)激盤(pán)模型的求解，得到1、2兩側(cè)的質(zhì)量、能量以及間斷的動(dòng)量通量，用于上下游流場(chǎng)的求解。由于激盤(pán)模型的建立包含了非堵塞工況的假設(shè)，因此該前后緣間斷處理不能適用于堵塞工況。同時(shí)由于沒(méi)有包含激波關(guān)系式，該模型同樣不能對(duì)激波進(jìn)行描述。

與Taddei和Larocca[20]的做法類(lèi)似，Pacciani等[22]對(duì)葉片前后緣間斷問(wèn)題進(jìn)行了研究。他們基于對(duì)平均葉片流面進(jìn)行修改的思想，提出一種時(shí)間相關(guān)的適應(yīng)性公式對(duì)流面進(jìn)行更新，如式(15)所示。

(15)

圖4 前后緣的激盤(pán)模型[21]Fig.4 Model of AD at leading or trailing edge[21]

Nigmatullin和Ivanov[23]同樣較早關(guān)注這個(gè)問(wèn)題，并對(duì)其進(jìn)行了系統(tǒng)且深入的研究，提出了葉片前后緣任意間斷分解的處理模型，較好地解決了該問(wèn)題。如圖5和圖6所示[23]，由于在無(wú)葉區(qū)(軸對(duì)稱區(qū))和葉片區(qū)采用了不同的控制方程，因此前后緣網(wǎng)格界面的法線方向在軸對(duì)稱側(cè)和葉片側(cè)并不相同。前后緣網(wǎng)格線在S1流面上可分解為CA、CB。圖5和圖6中：x為CA面的法向量;ξ為CB面法向量；N為平均流面的法向量，與ξ相互垂直；M為AB的法向量；l為AB的切向量。在前緣，2處和1處分別表示軸對(duì)稱側(cè)和葉片側(cè)，在尾緣則相反。2處參數(shù)與L處參數(shù)，R處參數(shù)與1處參數(shù)均通過(guò)精確的波特征關(guān)系相聯(lián)系；在前緣，L處與R處通過(guò)質(zhì)量守恒、能量守恒以及熵守恒關(guān)系關(guān)聯(lián)；在尾緣處則需要多補(bǔ)充一個(gè)方程，一般為AB向的動(dòng)量方程。

圖5 前緣任意間斷模型[23]Fig.5 Arbitrary discontinuity model at leading edge[23]

圖6 尾緣任意間斷模型[23]Fig.6 Arbitrary discontinuity model at trailing edge[23]

該前后緣間斷處理模型實(shí)際上相當(dāng)于對(duì)一般黎曼問(wèn)題中的中間波進(jìn)行修改，這里的中間波不再是接觸間斷，而是在L和R之間考慮氣流轉(zhuǎn)折和流動(dòng)圖景的一種“中間”波，在正常工況下具有熵守恒的性質(zhì)。同時(shí)Nigmatullin和Ivanov[23]還給出了前緣存在脫體激波、前緣堵塞、尾緣堵塞以及出現(xiàn)回流等各種特殊情況的處理方法，使得在任意流態(tài)下前后緣間斷問(wèn)題均能被較好地處理。

目前在CFD通流方法前后緣間斷問(wèn)題的處理中，Baralon等[11]的平均葉片流面修改方法雖然存在一些缺點(diǎn)，但由于較為簡(jiǎn)單和易于實(shí)現(xiàn)，得到了較為廣泛的應(yīng)用。Nigmatullin和Ivanov[23]提出的前后緣任意間斷分解模型是對(duì)前后緣間斷問(wèn)題在物理流動(dòng)圖景下的一種解析解，可以適用于各種工況的計(jì)算。Taddei和Larocca[21]的激盤(pán)處理可以看做是任意間斷分解方法中的正常工況處理。因此任意間斷分解法是一種較為完善的前后緣處理模型，但由于較為復(fù)雜，采用較少。

3.4 激波的捕獲

通常認(rèn)為CFD通流方法相比于流線曲率法和流函數(shù)法最大的優(yōu)勢(shì)在于捕獲激波的能力。然而，在正問(wèn)題和反問(wèn)題中，CFD通流方法表現(xiàn)出來(lái)的激波捕獲能力以及捕獲激波的性質(zhì)并不相同。Sturmayr和Hirsch[14]采用CFD通流方法對(duì)Rotor67分別進(jìn)行正問(wèn)題和反問(wèn)題求解，并對(duì)流場(chǎng)馬赫數(shù)進(jìn)行了對(duì)比，發(fā)現(xiàn)反設(shè)計(jì)得到的流場(chǎng)分布中沒(méi)有出現(xiàn)激波。Baralon等[11]對(duì)通流計(jì)算中的激波問(wèn)題進(jìn)行了詳細(xì)分析。結(jié)果發(fā)現(xiàn)，對(duì)于正問(wèn)題，須給定連續(xù)光滑的流動(dòng)角度分布，此時(shí)朗金-雨貢紐關(guān)系描述的是流向的正激波關(guān)系，因此當(dāng)流向馬赫數(shù)大于1時(shí)，CFD通流方法能夠捕捉到如圖7(a)所示的流向正激波。

對(duì)于設(shè)計(jì)問(wèn)題，需要給定連續(xù)光滑的周向速度分布，此時(shí)朗金-雨貢紐關(guān)系式描述的是軸向正激波，因此CFD通流方法只有在軸向馬赫數(shù)(考慮徑向速度分量時(shí)為子午馬赫數(shù))大于1時(shí)才能捕捉到如圖7(b)所示的軸對(duì)稱正激波。由于在傳統(tǒng)葉輪機(jī)中，子午馬赫數(shù)很難達(dá)到1以上，因此可以認(rèn)為在反問(wèn)題計(jì)算中，CFD通流方法不具備自動(dòng)捕獲激波的能力，因此需要添加相應(yīng)的激波損失模型。

圖7 S1流面的流向正激波，軸對(duì)稱正激波 Fig.7 Normal shock and axisymmetric shock on S1stream surface

3.5 周向應(yīng)力

在通流計(jì)算中，為了得到更加精準(zhǔn)的流場(chǎng)模擬結(jié)果，需要對(duì)流動(dòng)的三維特征以及非定常效應(yīng)進(jìn)行模擬。傳統(tǒng)方法均采用經(jīng)驗(yàn)關(guān)系式，適用的普遍性較差。為實(shí)現(xiàn)CFD通流方法的精細(xì)化預(yù)測(cè)，減小對(duì)經(jīng)驗(yàn)關(guān)系的依賴，必須進(jìn)一步對(duì)CFD通流方法進(jìn)行完善。Adamczyk[39]針對(duì)多級(jí)葉輪機(jī)非定常流動(dòng)提出了通道平均理論，通過(guò)對(duì)非定常Navier-Stokes方程進(jìn)行一系列平均，利用雷諾平均消除湍流脈動(dòng)，用時(shí)間平均消除非定常影響，用通道平均消除非周期性的影響，依次產(chǎn)生了雷諾應(yīng)力、確定應(yīng)力和通道平均應(yīng)力等附加項(xiàng)。在此基礎(chǔ)上進(jìn)一步進(jìn)行周向平均，便會(huì)產(chǎn)生周向應(yīng)力項(xiàng)，該項(xiàng)反映了流動(dòng)周向非均勻性的平均效應(yīng)。Baralon等[13]以Rotor67為例對(duì)周向應(yīng)力項(xiàng)進(jìn)行了研究，發(fā)現(xiàn)周向應(yīng)力項(xiàng)對(duì)流場(chǎng)的影響大于黏性應(yīng)力項(xiàng)，產(chǎn)生周向不均勻的主要原因包括間隙泄漏渦、通道激波、激波-邊界層干涉以及輪轂角區(qū)失速等物理現(xiàn)象。Baralon等還對(duì)不同的應(yīng)力項(xiàng)進(jìn)行了具體的對(duì)比，并在通流計(jì)算中引入周向應(yīng)力項(xiàng)，計(jì)算結(jié)果表明相比一般的CFD通流方法，考慮周向應(yīng)力的CFD通流方法與三維結(jié)果吻合得更好。Simon等[17]以單級(jí)壓氣機(jī)為算例，對(duì)不同通流模型進(jìn)行計(jì)算比較，包括只計(jì)入無(wú)黏葉片力、計(jì)入無(wú)黏葉片力和黏性應(yīng)力、計(jì)入無(wú)黏葉片力和周向應(yīng)力以及計(jì)入全部模型。研究結(jié)果表明，在CFD通流計(jì)算中無(wú)黏葉片力的作用占主要地位，周向應(yīng)力作用大于黏性應(yīng)力并且隨著負(fù)荷的提高而增大。Thomas和Léonard[40-42]首先對(duì)一低速壓氣機(jī)級(jí)中間展向位置處的S1面流動(dòng)進(jìn)行了研究，證實(shí)周向應(yīng)力項(xiàng)具有諧波特征。通過(guò)使用擾動(dòng)場(chǎng)的傅里葉級(jí)數(shù)對(duì)周向應(yīng)力進(jìn)行模擬，并以附加項(xiàng)的形式代入到了一維計(jì)算中，在一維流場(chǎng)中重現(xiàn)了二維流動(dòng)特征。Thomas和Léonard進(jìn)一步將周向應(yīng)力通過(guò)三維擾動(dòng)場(chǎng)的傅里葉級(jí)數(shù)進(jìn)行近似，并通過(guò)諧波對(duì)近似周向應(yīng)力進(jìn)行重構(gòu)，與CFD通流方法進(jìn)行耦合求解，成功地反映流動(dòng)的周向非均勻性，驗(yàn)證了諧波重構(gòu)的周向應(yīng)力在CFD通流方法中對(duì)流動(dòng)三維特征的再現(xiàn)能力。Thomas和Léonard還使用非線性諧波結(jié)合浸入邊界方法對(duì)周向應(yīng)力進(jìn)行了建模。浸入邊界方法用來(lái)構(gòu)建連續(xù)的力場(chǎng)以取代真實(shí)葉片，解決葉片表面無(wú)滑移邊界條件與傅里葉級(jí)數(shù)表達(dá)對(duì)葉片表面流動(dòng)連續(xù)性要求的沖突。該模型被應(yīng)用于無(wú)黏圓柱算例，并與解析解進(jìn)行了對(duì)比，結(jié)果吻合較好。該方法模擬周向不均勻項(xiàng)時(shí)，諧波法的級(jí)數(shù)越高，獲得的周向脈動(dòng)動(dòng)能與實(shí)際狀態(tài)越接近。

4 結(jié)論及展望

CFD通流計(jì)算方法可以完成傳統(tǒng)流線曲率通流方法的設(shè)計(jì)和分析功能；同時(shí)還具有更多優(yōu)點(diǎn)：精確捕捉激波，良好的堵塞適應(yīng)性，能夠提供更為精細(xì)的流場(chǎng)結(jié)構(gòu)，更加正確的幾何描述，也可以更方便地?cái)U(kuò)展到帶有雙涵道、引氣、冷卻以及葉尖間隙等復(fù)雜結(jié)構(gòu)的情況。

在正問(wèn)題中，CFD通流方法對(duì)于堵塞因子以及流動(dòng)角度的分布非常敏感；而對(duì)于反問(wèn)題則沒(méi)有此類(lèi)問(wèn)題；目前的CFD通流方法在傳統(tǒng)的葉輪機(jī)正問(wèn)題中可以捕捉流向正激波，在反問(wèn)題中只能捕獲軸向正激波。除此之外，無(wú)黏葉片力模型對(duì)于CFD通流方法求解得到的流場(chǎng)參數(shù)分布具有決定性影響，目前的葉片力模型仍然稍顯粗糙，計(jì)算結(jié)果與三維真實(shí)流場(chǎng)相差較大?；诟男屠杪鼏?wèn)題的任意間斷分解方法為解決前后緣間斷問(wèn)題提供了很好的基礎(chǔ)，但在理論上還需進(jìn)一步驗(yàn)證和完善。上述問(wèn)題是目前制約著CFD通流方法在工程實(shí)際中大規(guī)模應(yīng)用的主要問(wèn)題。

未來(lái)下列問(wèn)題的研究與解決將是建立高精度CFD通流方法的關(guān)鍵：首先，需要建立法向堵塞因子的標(biāo)準(zhǔn)計(jì)算方法，排除人為調(diào)整因素，得到準(zhǔn)確的流場(chǎng)幾何描述，以期更加正確地反應(yīng)真實(shí)幾何對(duì)流場(chǎng)的物理堵塞效應(yīng)。其次，通過(guò)添加邊界層厚度實(shí)現(xiàn)對(duì)堵塞更為精確的描述，同時(shí)需要研究邊界層厚度的添加對(duì)于葉片力計(jì)算的影響；其三，對(duì)跨、超聲速流動(dòng)中激波的三維圖景進(jìn)行建模描述，克服CFD通流方法正反問(wèn)題中分別只能捕獲流向、軸對(duì)稱正激波損失的缺點(diǎn)；其四，將基于改型黎曼問(wèn)題的前后緣任意間斷分解模型與流面修正相結(jié)合，實(shí)現(xiàn)通流計(jì)算中對(duì)葉片流動(dòng)分離工況的處理。此外，以通道平均理論為基礎(chǔ)的高階CFD通流方法也是值得關(guān)注的方向。

與傳統(tǒng)的流線曲率法類(lèi)似，影響CFD通流方法計(jì)算精度的最主要因素在于所采用的損失和落后角(滑移)模型。相關(guān)的研究較多，并且目前尚無(wú)統(tǒng)一的高精度模型存在，因此沒(méi)有在本文中進(jìn)行討論。需要指出的是，高精度的通流CFD計(jì)算是本文所討論的堵塞、激波、葉片力、前后緣間斷等諸因素與合適的損失和落后角模型全局耦合匹配的結(jié)果。

CFD通流方法相對(duì)于傳統(tǒng)流線曲率通流方法的優(yōu)點(diǎn)決定了其不僅在傳統(tǒng)的葉輪機(jī)部件設(shè)計(jì)和分析中能夠占有一席之地，而且在下一代先進(jìn)航空發(fā)動(dòng)機(jī)或燃?xì)廨啓C(jī)的部件設(shè)計(jì)和分析、整機(jī)穩(wěn)態(tài)數(shù)值模擬，甚至是過(guò)渡態(tài)性能模擬中，都具有相當(dāng)?shù)膬?yōu)勢(shì)與應(yīng)用潛力。目前中國(guó)在這方面取得了一些進(jìn)展，但主要還是對(duì)世界研究水平的跟隨，CFD通流方法的潛力還沒(méi)有被完全發(fā)揮出來(lái)。通過(guò)進(jìn)一步的發(fā)展與完善，該方法有希望成為葉輪機(jī)械以及燃?xì)廨啓C(jī)設(shè)計(jì)和分析的標(biāo)準(zhǔn)工具。

[1] WU Z H. A general theory of three-dimensional flow in subsonic and supersonic turbomachines of axial-, radial-, and mixed-flow types: NACA TN 2604[R]. Washington， D.C.: NASA, 1952.

[2] 蔣浩興. 國(guó)外發(fā)展風(fēng)扇/壓氣機(jī)設(shè)計(jì)體系的一些經(jīng)驗(yàn)和啟示[J]. 航空發(fā)動(dòng)機(jī), 2001(2): 45-51.

JIANG H X. The experiences and enlightenments of foreign fan/compressor design system development[J]. Aeroengine, 2001(2): 45-51 (in Chinese).

[3] PETROVIC M V, WIEDERMANN A. Fully coupled through-flow method for industrial gas turbine analysis: GT2015-42111[R]. New York: ASME, 2015.

[4] NOVAK R A. Streamline curvature computing procedures for fluid-flow problems[J]. Journal of Engineering for Power, 1967, 89(4): 478-490.

[5] MARSH H. A digital computer program for the through-flow fluid mechanics in an arbitrary turbomachine using a matrix method: RM 3509[R]. London: National Gas Turbine Establishment, 1968.

[6] SPURR A. The prediction of 3D transonic flow in turbomachinery using a combined throughflow and blade-to-blade time marching method[J]. International Journal of Heat and Fluid Flow, 1980, 2(4): 189-199.

[7] DAWES W N. Toward improved throughflow capability: The use of three-dimensional viscous flow solvers in a multistage environment[J]. Journal of Turbomachinery, 1992, 114(1): 8-17.

[8] YAO Z, HIRSCH C. Throughflow model using 3D Euler or Navier-Stokes solver[J]. VDI Berichte, 1995, 1185: 51.

[9] DAMLE S V, DANG T Q, REDDY D R. Throughflow method for turbomachines applicable for all flow regimes[J]. Journal of Turbomachinery, 1997, 119(2): 256-262.

[10] DAMLE S V. Throughflow method for turbomachines using Euler solvers: AIAA-1996-0010[R]. Reston, VA: AIAA, 1996.

[11] BARALON S, ERIKSSON L E, H?LL U. Validation of a throughflow time-marching finite-volume solver for transonic compressors: 98-GT-47[R]. New York: ASME, 1998.

[12] BARALON S, H?LL U, ERIKSSON L E. Viscous throughflow modelling of transonic compressors using a time-marching finite volume solver[C]//the 13th International Symposium on Air Breathing Engines. Reston: AIAA, 1997: 502-510.

[13] BARALON S, ERIKSSON L E, H?LL U. Evaluation of higher-order terms in the throughflow approximation using 3D Navier-Stokes computations of a transonic compressor rotor: 99-GT-74[R]. New York: ASME, 1999.

[14] STURMAYR A, HIRSCH C. Throughflow model for design and analysis integrated in a three-dimensional Navier-Stokes solver[J]. Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, 1999, 213(4): 263-273.

[15] SIMON J F, LéONARD O. A throughflow analysis tool based on the Navier-Stokes equations[C]//Proceedings of ETC 6th European Conference on Turbomachinery. Florence: European Turbomachinery Society, 2005: 7-11.

[16] SIMON J F, LéONARD O. Modeling of 3-D losses and deviations in a throughflow analysis tool[J]. Journal of Thermal Science, 2007, 16(3): 208-214.

[17] SIMON J F, THOMAS J P, LéONARD O. On the role of the deterministic and circumferential stresses in throughflow calculations[J]. Journal of Turbomachinery, 2009, 131(3): 031019.

[18] PERSICO G, REBAY S. A penalty formulation for the throughflow modeling of turbomachinery[J]. Computers & Fluids, 2012, 60(10): 86-98.

[19] TADDEI S R, LAROCCA F. Axisymmetric design of axial turbomachines: An inverse method introducing profile losses[J]. Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, 2008, 222(6): 613-621.

[20] TADDEI S R, LAROCCA F. CFD-based analysis of multistage throughflow surfaces with incidence[J]. Mechanics Research Communications, 2013, 47(47): 6-10.

[21] TADDEI S R, LAROCCA F. An actuator disk model of incidence and deviation for RANS-based throughflow analysis[J]. Journal of Turbomachinery, 2014, 136(2): 021001.

[22] PACCIANI R, RUBECHINI F, MARCONCINI M, et al. A CFD-based throughflow method with an explicit body force model and an adaptive formulation for the S2 streamsurface[J]. Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, 2016,230(1): 16-28.

[23] NIGMATULLIN R Z, IVANOV M J. The mathematical models of flow passage for gas turbine engines and their components: LS 198[R]. Neuillysurseine: AGARD, 1994.

[24] 袁寧, 張振家, 顧中華, 等. 渦噴發(fā)動(dòng)機(jī)壓氣機(jī)三種S2流面計(jì)算程序的比較[J]. 推進(jìn)技術(shù), 1998, 19(1): 51-57.

YUAN N, ZHANG Z J, GU Z H, et al. A comparison of three kinds of calculation program of S2 stream surface in the compressor of aero-engine[J]. Journal of Propulsion Technology, 1998, 19(1): 51-57 (in Chinese).

[25] 季路成, 孟慶國(guó), 周盛. 葉輪機(jī)通流計(jì)算的時(shí)間推進(jìn)方法[J]. 航空動(dòng)力學(xué)報(bào), 1999, 14(1): 23-26.

JI L C, MENG Q G, ZHOU S. Time-marching method for through-flow computation of turbomachinery[J]. Journal of Aerospace Power, 1999, 14(1): 23-26 (in Chinese).

[26] 施發(fā)樹(shù), 劉興洲. 多部件模型在全尺寸小型雙函道渦扇發(fā)動(dòng)機(jī)氣流數(shù)值模擬中的應(yīng)用[J]. 推進(jìn)技術(shù), 1998, 19(4): 22-26.

SHI F S, LIU X Z. Multicomponent models in application to numerical simulation of a small full-sized by-pass turbofan engine[J]. Journal of Propulsion Technology, 1998, 19(4): 22-26 (in Chinese).

[27] 施發(fā)樹(shù). 一體化彈用小渦扇發(fā)動(dòng)機(jī)系統(tǒng)的氣動(dòng)熱力數(shù)值模擬[D]. 南京: 南京航空航天大學(xué), 1999.

SHI F S. Aero-thermodynamic numerical simulation of integrated small turbofan engine system[D]. Nanjing: Nanjing University of Aeronautics and Aerospace, 1999 (in Chinese).

[28] 于龍江, 陳美寧, 樸英. 航空發(fā)動(dòng)機(jī)整機(jī)準(zhǔn)三維流場(chǎng)仿真[J]. 航空動(dòng)力學(xué)報(bào), 2008, 23(6): 1008-1013.

YU L J, CHEN M N, PIAO Y. Quasi-3D simulation of aero engine full flow field[J]. Journal of Aerospace Power, 2008, 23(6): 1008-1013 (in Chinese).

[29] 曹志鵬, 劉大響, 桂幸民, 等. 某小型渦噴發(fā)動(dòng)機(jī)二維數(shù)值仿真[J]. 航空動(dòng)力學(xué)報(bào), 2009, 24(2): 439-444.

CAO Z P, LIU D X, GUI X M, et al. Two dimensional numerical simulation of small turbojet engine[J]. Journal of Aerospace Power, 2009, 24(2): 439-444 (in Chinese).

[30] 金東海, 桂幸民. 某渦扇發(fā)動(dòng)機(jī)考慮級(jí)間引氣的二維數(shù)值模擬[J]. 航空動(dòng)力學(xué)報(bào), 2011, 26(6): 1346-1351.

JIN D H, GUI X M. Two dimensional numerical simulation of a turbofan engine with air bleeding in compressor[J]. Journal of Aerospace Power, 2011, 26(6): 1346-1351 (in Chinese).

[31] 李德英， 宋彥萍， 陳浮, 等. 任意曲線坐標(biāo)系Euler方程S2流面的計(jì)算方法[J]. 西安交通大學(xué)學(xué)報(bào), 2015, 49(7): 42-48.

LI D Y, SONG Y P, CHEN F, et al. Euler S2 stream surface calculation for arbitrary curvilinear coordinate system[J]. Journal of Xi’an Jiaotong University, 2015, 49(7): 42-48 (in Chinese).

[32] WAN K, JIN H, JIN D, et al. Influence of non-axisymmetric terms on circumferentially averaged method in fan/compressor[J]. Journal of Thermal Science, 2013, 22(1): 13-22.

[33] 萬(wàn)科, 朱芳, 金東海, 等. 周向平均方法在某風(fēng)扇/增壓級(jí)分析中的應(yīng)用[J]. 航空學(xué)報(bào), 2014, 35(1): 132-140.

WAN K, ZHU F, JIN D H, et al. Application of circumferentially averaged method in fan/booster[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(1): 132-140 (in Chinese).

[34] SIMON J F. Contribution to throughflow modelling for axial flow turbomachines[D]. Liege: University of Liege, 2007.

[35] MARBLE F. Three-dimensional flow in turbomachines[J]. High Speed Aerodynamics and Jet Propulsion, 1964, 10: 83-166.

[36] DENTON J D. Throughflow calculations for transonic axial flow turbines[J]. Journal of Engineering for Power, 1978, 100(2): 212-218.

[37] DANG T Q, WANG T. Design of multi-stage turbomachinery blading by the circulation method: actuator duct limit: 92-GT-286[R]. New York: ASME, 1992.

[38] BOSMAN C, MARSH H. An improved method for calculating the flow in turbo-machines, including a consistent loss model[J]. Journal of Mechanical Engineering Science, 1974, 16(1): 25-31.

[39] ADAMCZYK J J. Model equation for simulating flows in multistage turbomachinery[J]. Lecture Series—van Kareman Institute for Fluid Dynamics, 1996, 5: N1-N28.

[40] THOMAS J P, LéONARD O. Investigating circumferential non-uniformities in throughflow calculations using a harmonic reconstruction: GT2008-50328[R]. New York: ASME, 2008.

[41] THOMAS J P, LéONARD O. Towards a high order throughflow: Part Ⅰ—Investigating the effectiveness of a harmonic reconstruction for 3D flows: GT2010-22841[R]. New York: ASME, 2010.

[42] THOMAS J P, LéONARD O. Toward a high order throughflow—Investigation of the nonlinear harmonic method coupled with an immersed boundary method for the modeling of the circumferential stresses[J]. Journal of Turbomachinery, 2012, 134(1): 011017.

(責(zé)任編輯： 鮑亞平， 張晗)

*Corresponding author. E-mail: jinguang_yang@dlut.edu.cn

Time marching based throughflow method: Current status and future development

YANG Jinguang1,*, WANG Chunxue2, WANG Dalei2, SHAO Fuyong2, YANG Chen3, WU Hu3

1.SchoolofEnergyandPower,DalianUniversityofTechnology,Dalian116024,China2.BeijingPowerMachineryResearchInstitute,ChinaAerospaceScienceandIndustryCorporation,Beijing100074,China3.SchoolofPowerandEnergy,NorthwesternPolytechnicalUniversity,Xi’an710072,China

Aiming to have an exploration on application potential of the throughflow computation method for turbomachines based on the time marching, the research progress home and abroad is concluded, and several key issues of application of this method are generalized, including geometric description of blade blockage, simulation of blade force, discontinuity problems in blade leading and trailing edges, and shock capture. Different computation models are also compared. It is summarized that the throughflow computation method based on the time marching has many advantages compared to the traditional throughflow method, and thus has the potential to be a competitive tool in advanced gas turbine engine components design and analysis, as well as in total engine steady and transient state simulation. This method is expected to have much more improvements in the aspects mentioned above, although it has achieved great progress. It is believed that this method will be adopted as a new standard tool in turbomachinery design and analysis with further developments and improvements.

turbomachine； throughflow； time marching； blade blockage； blade force； discontinuity at blade leading edge and trailing edge； shock

2016-11-28； Revised： 2016-12-19； Accepted： 2017-01-11； Published online： 2017-03-09 09:05

URL： www.cnki.net/kcms/detail/11.1929.V.20170309.0905.004.html

the Fundamental Research Funds for Central Universities (DUT15RC(3)035)

V231.3

A

1000-6893(2017)09-520996-13

2016-11-28； 退修日期： 2016-12-19； 錄用日期： 2017-01-11； 網(wǎng)絡(luò)出版時(shí)間： 2017-03-09 09:05

www.cnki.net/kcms/detail/11.1929.V.20170309.0905.004.html

中央高?；究蒲袠I(yè)務(wù)費(fèi)專(zhuān)項(xiàng)資金 (DUT15RC(3)035)

*通訊作者.E-mail: jinguang_yang@dlut.edu.cn

楊金廣， 王春雪， 王大磊， 等． 基于時(shí)間推進(jìn)的通流計(jì)算方法： 現(xiàn)狀及展望 [J]． 航空學(xué)報(bào), 2017, 38(9): 520996. YANG J G, WANG C X, WANG D L, et al. Time marching based throughflow method: Current status and future development[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(9): 520996.

http://hkxb.buaa.edu.cn hkxb@buaa.edu.cn

10.7527/S1000-6893.2017.620996