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

1.大連理工大學 能源與動力學院， 大連 116024 2.中國航天科工集團 北京動力機械研究所， 北京 100074 3.西北工業(yè)大學 動力與能源學院， 西安 710072

基于時間推進的通流計算方法：現狀及展望

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

1.大連理工大學 能源與動力學院， 大連 116024 2.中國航天科工集團 北京動力機械研究所， 北京 100074 3.西北工業(yè)大學 動力與能源學院， 西安 710072

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

葉輪機械； 通流； 時間推進； 葉片堵塞； 葉片力； 葉片前后緣間斷； 激波

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

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

本文對國內外基于時間推進的通流計算方法進行總結，對該方法在具體實施中的關鍵問題處理方法進行討論，最后對其未來的發(fā)展方向進行展望。

1 國內外研究進展

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

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

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

2 CFD通流方法基本計算模型

通流計算的基本思想在于認為葉輪機中的流動在無葉區(qū)滿足軸對稱假設，在葉片通道內則沿S2流面流動。因此，CFD通流計算模型一般通過對三維的Euler/Navier-Stokes方程在相應方向進行平均得到，從而實現從全三維的流動問題(圖1(a))向平均葉片S2流面求解(圖1(b))的降維簡化[18]。如果此時進一步在流場中施加適當形式的體積力來模擬葉片的作用，迫使當地流動速度方向滿足與平均葉片流面相切，則問題可以進一步簡化到子午平面內的流場計算問題(圖1(c))。

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

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

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

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

前兩類CFD通流方法均在子午平面內進行流場計算。以周向平均后的軸對稱Euler方程為例，一般的控制方程表達形式為

(1)

式中：

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

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

3 CFD通流方法中的若干關鍵問題

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

3.1 葉片堵塞的幾何描述

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

(2)

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

dmax≈dnmax/cosβnmax

(3)

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

為了更好地估計實際葉片厚度對流動的影響，Baralon等[11]認為采用法向堵塞因子比周向堵塞因子更為合適，并給出了相應的法向堵塞因子計算方法，表達式為

(4)

表1 堵塞流量對堵塞因子分布的敏感度[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為相應的葉片厚度。該表達式僅限于計算點A和點B之間的堵塞因子，如圖2[11]所示，從上方葉片前緣點向平均流面(在S1面上為平均流線)作垂線，垂線與平均流面交點即為A點；同理，點B由下方葉片尾緣點向平均流面作垂線得到；點P為點A和點B之間的平均流線上任一點。對于前緣到點A和尾緣到點B的部分，堵塞因子定義為從bnA和bnB線性變化到1?？梢园l(fā)現，該方法定義的法向堵塞因子與周向堵塞因子差別較為明顯，對于某一展向葉片截面，計算得到的流道面積大小和喉道軸向位置均不相同，這些差異會對流場計算結果造成較為明顯的影響，尤其在激波存在的跨、超聲速流動中。該種堵塞因子定義方法存在實際操作復雜、前后緣部分分布過于簡單、精確度較差的缺點。

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

bn=bcosβ

(5)

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

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

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

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

3.2 葉片力

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

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

(6)

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

單位法向量

則FB可表示為

FB=B(r,z)·ng

(7)

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

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

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

(8)

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

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

(9)

式中：λ為典型響應時間平方的倒數，s-2。ρfBφ的響應時間應該和流場計算中的當地網格時間步長處于同一量級，因此λ可以表達為

式中：Δx為典型網格尺寸；Δt為當地網格時間步長。該方程可以理解為一種壓力梯度，一旦流面相切條件不符合就會產生，該方法相當于對周向葉片力在時間上的一種離散。Simon和Léonard[15]采用了類似的葉片力模擬方法，只不過采用葉片力的模作為求解變量。萬科等[33]指出該種葉片力模型在實際使用過程中，λ的取值對迭代收斂的影響較大，取值不當會影響收斂速度，嚴重時甚至會導致迭代不收斂。

(10)

式中：C為松弛因子。

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

B(r,z)=-KW·ng

(11)

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

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

(12)

式中：J為坐標系變換的雅克比行列式；ζz、ζr和ζφ為與ζ=const網格面對應的軸向、徑向和周向的度量系數；S2、S3和S4分別為軸向、徑向和周向動量方程中不考慮無黏葉片力項的殘差。該方法簡單通用，不需要引入偽時間項，亦避免了因此帶來的常系數取值問題，同時具有較好的收斂性。

對于黏性力，一般可以寫成如式(13)所示的表達式。

(13)

式中：黏性力方向與當地速度方向相反。L(r,z)為一正值，代表黏性損失的大小。對于黏性損失的計算，Bosman和Marsh[38]提出通過給定流場的熵分布進行計算，該損失模型簡單且使用廣泛。由此L(r,z)可以和熵分布聯系起來，即

(14)

式中：T為當地靜溫；s為熵。

3.3 葉片前、后緣間斷問題

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

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

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

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

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

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

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

(15)

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

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

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

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

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

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

3.4 激波的捕獲

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

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

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

3.5 周向應力

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

4 結論及展望

CFD通流計算方法可以完成傳統流線曲率通流方法的設計和分析功能；同時還具有更多優(yōu)點：精確捕捉激波，良好的堵塞適應性，能夠提供更為精細的流場結構，更加正確的幾何描述，也可以更方便地擴展到帶有雙涵道、引氣、冷卻以及葉尖間隙等復雜結構的情況。

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

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

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

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

[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] 蔣浩興. 國外發(fā)展風扇/壓氣機設計體系的一些經驗和啟示[J]. 航空發(fā)動機, 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ā)動機壓氣機三種S2流面計算程序的比較[J]. 推進技術, 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] 季路成, 孟慶國, 周盛. 葉輪機通流計算的時間推進方法[J]. 航空動力學報, 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ā)樹, 劉興洲. 多部件模型在全尺寸小型雙函道渦扇發(fā)動機氣流數值模擬中的應用[J]. 推進技術, 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ā)樹. 一體化彈用小渦扇發(fā)動機系統的氣動熱力數值模擬[D]. 南京: 南京航空航天大學, 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ā)動機整機準三維流場仿真[J]. 航空動力學報, 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ā)動機二維數值仿真[J]. 航空動力學報, 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ā)動機考慮級間引氣的二維數值模擬[J]. 航空動力學報, 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] 李德英， 宋彥萍， 陳浮, 等. 任意曲線坐標系Euler方程S2流面的計算方法[J]. 西安交通大學學報, 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] 萬科, 朱芳, 金東海, 等. 周向平均方法在某風扇/增壓級分析中的應用[J]. 航空學報, 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.

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

*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； 網絡出版時間： 2017-03-09 09:05

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

中央高?；究蒲袠I(yè)務費專項資金 (DUT15RC(3)035)

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

楊金廣， 王春雪， 王大磊， 等． 基于時間推進的通流計算方法： 現狀及展望 [J]． 航空學報, 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