葉輪機械全環非定常大規模并行模擬程序設計

張 健唐 靜邱 名鄧有奇龔小權

(中國空氣動力研究與發展中心 計算空氣動力研究所,綿陽 621000)

0 引 言

葉輪機械內部流動非常復雜,其本質為三維黏性非定常流動。計算流體力學(Computational Fluid Dynamics,CFD)技術飛速發展,其具備計算耗時短,成本相對較低,并且可以多維度地模擬出葉輪機內部復雜的非定常流動細節的優點,因而越來越多地用于葉輪機械的氣動設計和分析當中。

葉輪機內部非定常流動的來源有多種,動靜葉片排的相對轉動、激波邊界層干擾、尾跡傳播、二次流等都是造成流動不穩定的原因。通過數值模擬對上述非定常現象進行研究的關鍵在于如何建立模擬轉靜葉片相對運動過程的模型。目前工業上一般采用混合平面(Mixing Plane)模型[1],該模型將轉/靜界面上下游相同展向高度處的通量通過周向平均后進行交換,將非定常計算簡化為對一個葉片流道的定常計算,雖然大大減小了計算量,然而卻無法捕捉到轉靜子之間相互干擾等非定常現象。采用葉片約化技術[2],將葉片在周向方向上按照一定比例進行幾何縮放,能夠使轉靜葉片數具有較大公約數進而約化為少量幾個葉片通道進行非定常計算,但是由于改變了實際的幾何尺寸,這種方法會存在較大誤差。He和Ning提出的非線性諧波法(Non-linear Harmonic Method,NLH)[3]可以看作是一種定常/非定常混合方法,其基本假設是流場的主要擾動是由于葉片通過頻率(Blade Passing Frequence,BPF)引起的,從而將流場變量分解成為時間平均值和多個不同頻率諧波的擾動組成,BPF的整數倍數分別代表了不同諧波。為了得到最終的流場解,除了要求解時間平均方程外,還需要求解各個不同頻率的諧波方程。NLH模型優點是在計算消耗可承受狀態下得到相對精確的非定常流場,但缺點是需要預先給定頻率,對于未知頻率流動和多頻率的組合模擬存在困難,特別是旋轉失速等非定常流動,NLH不可能解決。而且,隨著頻率數目增加,NLH的求解量也會劇增。在未來的航空發動機設計流程中,全環非定常數值模擬多級壓氣機/渦輪內部的復雜三維流動必然會成為一種趨勢和常規手段。NASA支持開發的流場求解器TURBO[4]就是一個精確模擬三維非定常發動機流場的軟件。Chen等利用TURBO全環數值模擬了壓氣機失速開始狀態的流場,并研究了相關的失速控制技術[5]。北航邵飛等利用商用軟件CFX全環計算了某高壓渦輪非定常流動,捕捉到了流場中存在的低頻擾動[6]。

由于計算量巨大,展開葉輪機械全環非定常研究具有較大挑戰性。隨著高性能計算機的飛速發展,計算機硬件逐漸不再是限制CFD計算的主要因素。我國的神威太湖之光[7]、天河二號[8]等超級計算機已經達到了每秒億億次量級的浮點數計算能力。數千萬甚至上億網格量的計算已不再是難題。美國NASA關于2030年CFD發展展望的報告[9]中,多次強調了高性能計算(HPC)在未來葉輪設計的重要作用。基于此,本文提出了一個大規模并行模擬程序設計策略,并針對葉輪機械的三維非定常流動進行了全環數值模擬,從而對更加深入研究發動機內部復雜流動機理提供保障和基礎。

本文依托氣動中心千萬億次高性能計算集群,在自主研發的三維混合網格流場求解器MFlow[10]的基礎上,開發構建了葉輪機定常模塊、非定常流動數值模擬模塊。滑移網格及通量守恒插值模塊的研發實現了航發葉輪機的全環非定常模擬。研究先通過一個典型單級壓氣機進行了程序的驗證與確認,然后利用1.5級渦輪算例進行了程序并行效率的測試,最后對一個3.5級壓氣機超大規模算例進行了嘗試,分析了動靜葉片相互干擾過程,評估了程序對于大規模復雜工程應用的適用能力。

1 計算方法模型

1.1 控制方程

流場求解采用三維非定常雷諾平均Navier-Stokes方程(URANS)在旋轉坐標系下的形式:

式中:Ω為控制體單元；S代表單元面。U為守恒變量,F I和F V分別是無黏通量和黏性通量,S T為科里奧利力和離心力引起的源項,具體形式如下:

式中:ρ為密度,u為相對速度矢量,E為總能,p為壓力,q為熱通量,τ為剪切應力張量,ω為旋轉角速度,r為徑向矢量。

1.2 數值方法

對于方程式(1),時間推進采用雙重時間步法,外迭代采用歐拉后插顯式推進,每一個時間步子迭代利用LU-SGS(Lower-Upper Symmetric Gauss-Seidel)方法[11]隱式迭代。空間離散基于控制體的有限體積法,采用Roe迎風格式[12]計算無黏通量,利用格林-高斯公式計算控制體單元的梯度,選擇Venkataknshnan限制器[13]將控制體的原始變量通過線性重構插值到控制體單元面左右,達到更高的空間精度。采用中心格式計算黏性通量,湍流模型選取Spalart-Allmaras(SA)一方程模型[14]。

1.3 邊界條件

對于進口邊界,給定總溫T0,總壓p0以及來流速度方向。向外發出的黎曼不變量R－由第一個內部邊界點外插計算得到,利用R－計算出絕對速度大小V,再根據速度方向將V分解為各個速度分量。

式中:v d代表邊界內部第一個單元速度矢量,c d為當地聲速,γ為比熱比,n為法向量,C p為定壓比熱。密度ρ和壓力p利用T0、p0、V之間的代數關系以及等熵關系計算得到。

對于內流出口,給定靜壓,其他變量外插得到。同時引入徑向平衡方程式(5),通過給定展向高度r處的靜壓,其他位置的靜壓根據密度ρ和周向速度vθ積分得到。

對于黏性物面,為絕熱壁,采用無滑移邊界條件。對于轉靜交界面,在每一個物理時間步,上下游的守恒變量隨著網格滑移動態完成交換。動靜葉排之間利用滑移網格進行耦合,通過將對向網格的守恒變量插值到本方虛擬網格上完成邊界條件的更新,如圖1所示。

圖1 滑移網格技術Fig.1 Sliding interface technique

被插值單元的值通過包裹單元的節點值進行反距離權插值平均方法得到。如式(6)所示,ωi,j為距離權重,d i,j為兩點距離。

2 并行程序設計

2.1 并行程序框架

圖2給出了整個程序的流程圖,在每一個非定常時間步完成方程求解、網格旋轉、建立映射關系的迭代循環,在每一個子迭代需要完成并行網格分區之間以及轉/靜界面之間的數據交換。

圖2 并行程序框架Fig.2 Parallel program frame

并行區域的劃分是在前處理過程中完成的,分區完成后,程序將記錄下各個分區的幾何信息,同時記錄下并行左右分區的映射關系。轉/靜界面的映射關系是在每一個非定常步求解開始之前完成。數據的并行傳值利用消息傳遞接口(Message Passing Interface,MPI)實現。利用MPI_Comm_Split方法將不同葉片排分區網格分到不同的子通信域,子通信域內完成常規并行邊界的通量傳遞,子通信域之間完成滑移邊界的變量插值與傳遞。

2.2 并行網格分區

并行網格分區利用圖切分軟件包METIS完成,METIS是由美國密西根大學G.Karypis和V.Kumar編寫的用于圖的分區和稀疏矩陣排序的串行包,提供k路多層劃分算法對混合網格進行分區[15]。給定包含n個節點V和m個邊E的圖G＝(V,E),將V劃分到k個子集中V＝[V1,V2,...,V k],使得|V i|≈n/k,且相連頂點屬于不同子集的邊的數量要最小化。

圖3給出了k路多層劃分過程,整個劃分過程分成coarsening,initial partitioning和uncoarsening三個部分。coarsening將圖的大小逐漸縮小,G1-＞G2-＞G3-＞G4。在G4階段執行k路劃分,然后在uncoarsening階段將圖中的原始節點映射到G4劃分的子集中。由于本文采用格心格式的有限體積法,圖定義中的節點和邊分別指網格中的網格單元和網格面。采用k路多層圖劃分,保證了各個子區域的總網格數比較均衡的同時盡量減少了分區的邊切割,達到了負載均衡目的同時,最大限度減少了分區間的數據通信量。

圖3 METIS的k路多層圖劃分過程Fig.3 Multilevel k-way partitioning of METIS

2.3 分區間數據交換

本文并行算法基于分布式大規模并行計算機系統,采用所有CPU都參與迭代計算的對等并行模式,MPI數據傳遞采用非阻塞通信,盡量加大計算和通信的重疊,每個分區網格只保存必需的局部網格數據和并行分區邊界上網格單元及頂點的對應關系。多級葉輪網格內部有兩種并行邊界,一種是同一葉片排內部,在采用METIS并行網格分區時形成的交界面,另一種是不同葉片排之間的滑移網格交界面。對于前者,并行邊界面左右網格的對應關系在前處理分區時就已經確定,所以只需要計算一次并記錄下來,以空間換時間,提高并行效率。對于后者,每一非定常步轉子葉片網格將進行旋轉,滑移邊界左右的對接信息發生更新,需要重新計算。插值網格之間的映射關系通過ADT(Alternating Digital Tree)方法[16]建立,該方法將空間單元進行有序排序后,構建為二叉樹數據,以便待插值點快速搜索查找[17],減少計算時間。

程序中,對于每一個分區所在的進程中,將該進程的并行邊界左右單元(物理單元以及虛擬單元)編號按照向相鄰分區發送和從相鄰分區接收的順序連續存儲在數組內存中(圖4),交換數據時只需要將相應數據按照單元編號存入或取出,并調用MPI數據傳遞接口完成。

圖4 交換數據存儲結構Fig.4 Transfer data's storage structure

3 單級壓氣機算例驗證

NASA Stage 35作為典型跨聲速壓氣機常用于壓氣機內流數值模擬驗證。表1列出了Stage 35的具體的設計參數。Stage 35在設計轉速17188.7 rpm、流量20.2 kg/s下產生1.8的總壓比,Reid和Moore早期對其做過試驗,可以得到詳細的幾何參數、工作條件和試驗數據[18]。本文對Stage 35算例分別利用混合平面法進行了單通道定常計算[19]、利用非線性諧波方法進行了“準”非定常計算、利用滑移面法進行了360°全環非定常計算,其中非線性諧波法是通過商用軟件Numeca完成的,因為只有一組轉靜干擾,所以周期性擾動個數設為1,每個擾動采用的諧波數為3個。最后將計算結果同試驗數據以及NASA TURBO軟件計算的360°全環非定常數據進行了對比分析。

表1 Stage 35設計參數Table 1 Design parameters of Stage 35

本次計算的網格如圖5所示,首先完成單通道的網格生成,第一層網格間距3×10-6m,y＋取1,同時考慮了葉尖間隙的影響,間隙展向分布了17個網格點。然后旋轉復制得到Stage 35全環網格。總網格量大概在3300萬左右。

本次的計算在中國空氣動力研究與發展中心的千萬億次集群上完成。并行計算采用512個核,集群芯片的處理器為國產的飛騰FT-1500A處理器,主頻2 GHz。每個時間步葉片旋轉1°,總共完整計算3圈。圖6給出了計算穩定后某一時刻50%葉高全環剖面的壓力分布云圖。

為了獲取壓氣機的整體性能并且評估數值模擬的精度,圖7和圖8分別給出了壓氣機在100%轉速下不同方法和模型計算得到的總壓比和絕熱效率隨流量變化的特性曲線,同時給出試驗數據進行對比。從圖看出,數值模擬得到的壓氣機整體性能和試驗數據吻合良好。相比較而言,本文的全環非定常計算的總壓比相對于文獻[20]中TURBO軟件全環的計算結果與試驗數據更加接近,尤其是在接近阻塞(試驗中的4004工況)時。本文計算得到最大流量與Reid和Moore試驗數據非常接近,表2給出了該狀態下的總壓比、總溫比和絕熱效率與試驗及TURBO計算結果的對比。可以看到MFlow的計算結果與試驗的誤差基本都在Reid和Moore提出的試驗不確定度范圍之內。文獻中沒有給出TURBO計算的絕熱效率特性曲線,相比于混合平面法和非線性諧波法,MFlow全環非定常計算的絕熱效率與試驗更加接近。在小流量工況下計算得到的壓比相對混合面法和NLH低,其原因可能上下游網格信息傳遞時由于插值精度不夠,導致數值耗散大[21]。從整體性能來看,本文程序的計算結果具有較高精度,滿足工程應用需求。

圖5 Stage 35計算網格Fig.5 Computational grid of Stage 35

圖7 Stage 35總壓特性曲線Fig.7 Total pressure ratio of Stage 35

圖8 Stage 35絕熱效率特性曲線Fig.8 Adiabatic efficiency of Stage 35

表2 Stage 35接近阻塞狀態氣動特性對比Table 2 Aerodynamic performance of Stage 35 near choke

圖9和圖10分別給出了利用不同模型計算得到的50%葉高切面壓力分布以及熵值云圖。通過對比可以看到不同模型計算得到的波系結構基本一致,但是混合平面法由于在轉/靜交界面進行了周向平均,所以壓力和熵值并不連續,尤其是熵增有一個很明顯的階躍。NLH方法的結果與非定常模擬的結果更加接近,能夠捕捉到轉/靜干擾現象。但是仔細觀察轉子產生的熵值進入靜子區域,在幾條高亮的尾跡中間還有兩條強度較小的尾跡。這是因為NLH方法是通過有限的幾個頻率諧波方程的定常解合成一個近似的非定常解,不同頻率會帶來不同波長和幅值的正弦振蕩,只有在諧波數足夠多的情況下才能完全消除擾動[22]。相比之下全環模擬的結果與實際物理現象更加接近。

圖9 不同方法計算壓力云圖比較Fig.9 Static pressure comparison of different methods

圖10 不同方法計算熵值云圖比較Fig.10 Entropy comparison of different methods

為研究相鄰葉片排之間的周期性干擾過程,在計算過程中對靜子尾跡接近出口某處設置一個壓力探測點,壓力隨時間周期變化結果如圖11所示,可以看到監測點的靜壓在一個周期后呈現出良好的周期性,說明非定常計算已經達到穩定。對其進行頻譜分析結果如圖12所示,可以看到計算很好的捕捉到了葉片通過頻率(BPF)及其諧波。

圖11 探測點壓力變化曲線Fig.11 Static pressure at the probe

圖12 探測點壓力頻譜分析結果Fig.12 Pressure spectrum

4 并行效率測試

為測試程序的并行計算效率,同時進一步檢驗程序對不同種類葉輪機械模擬的適應性,本節采用了一個稍大網格規模的1.5級渦輪算例進行了計算和分析。其基本幾何構型如圖13所示,全環轉子葉片共41個,轉速3500 rpm,靜子葉片36個[23]。圖14給出了計算所采用的其中一個葉片通道網格,整個3排葉片全環網格量共計約4000萬。

圖13 Aachen 1.5級渦輪幾何參數Fig.13 Geometry parameters of Aachen 1.5 turbine stage

圖14 Aachen 1.5級渦輪計算網格Fig.14 Computational grid of Aachen 1.5 turbine stage

計算不考慮渦輪的傳熱過程,每一個非定常步轉子大約轉動1度,每一個非定常步下子迭代200次。對同一狀態分別利用180、360、720以及1440個CPU核心數分別計算,從而進行并行效率的測試。圖15給出了50%展向環面的熵值云圖計算結果,反映了非定常流動損失的輸運過程。第一排靜葉產生熵的尾跡,遇到轉子葉片前緣被切斷后轉向進入轉子葉片通道并向下游傳輸,再與轉子葉片尾跡產生的熵混合,一起進入下游的第二排靜葉通道混合。從熵增的流動過程和在轉/靜界面的連續性來看,計算結果合理。圖16是利用不同CPU核數進行了一定步數非定常計算后殘差隨時間收斂曲線對比,可以明顯看到隨核數成倍增加,收斂時間幾乎成倍減少。圖17進一步給出了加速比曲線,可以看到720核時并行效率有90%,1440核并行效率依然能夠接近80%。利用1440核完成本算例一個狀態的計算,轉子完整旋轉一周計算時間不到24小時,顯示出程序具有較高的并行效率和加速比。

圖15 50%葉高環面熵值云圖Fig.15 Entropy contours at 50%span annulus

圖16 不同核數計算殘差隨時間收斂曲線Fig.16 Residual convergence using different number of cores

圖17 并行效率測試結果Fig.17 Parallel efficiency

5 超大規模算例

在完成程序的驗證與確認及性能分析后,利用程序對一個超大規模算例進行了嘗試。該算例是某壓氣機的前3.5級,設計轉速15000 rpm。圖18給出了本次計算所用的網格,網格總量接近2億,共分4096個核進行計算。

圖18 3.5級壓氣機全環網格Fig.18 Full annulus grid of 3.5 stage compressor

表3給出了在接近設計點的一個狀態下,分別利用混合平面法和全環非定常方法計算得到的壓氣機氣動性能對比。可以看到,由于混合平面法方法模型本身的缺陷,進出口流量誤差隨著級數的增多累積明顯,達到2.5%,而全環非定常的差別只有0.7%,其中可能是插值的精度損失造成的。除流量外,定常模型計算得到的總壓比、絕熱效率與非定常比差別較明顯,說明當壓氣機級數越多,定常模型計算結果的可信度就越低。

表3 3.5級壓氣機設計點計算結果Table3 Aerodynamic performance of 3.5 stage compressor at design point

圖19給出了計算得到的整機壓力分布云圖,圖20給出了三個不同時刻在前兩級葉片排在70%葉高切面的焓值局部放大云圖,可以幫助分析壓氣機內部復雜的動靜干擾過程。可以看到,轉子葉柵外伸激波會射入上游靜子葉柵通道,外伸激波在壓力面干擾強、傳播速度快、并形成反射,從而使得靜葉壓力面的附面層快速增厚,靜葉表面同時也會受到上游轉子渦脫落尾跡的影響,使得靜葉表面壓力存在劇烈的非定常波動。受上下游靜子尾跡及壓力波傳播干擾,轉子通道內的激波形態也會不斷發生變化。說明全環模擬能夠揭示更多葉輪機械內部復雜的三維非定常流動細節。

圖19 3.5級壓氣機壓力云圖Fig.19 Pressure contours of 3.5 stage compressor

6 總結與展望

本文研究開發了基于滑移網格技術、METIS網格分區技術和并行邊界處理及虛擬單元技術等的多級葉輪機械全環非定常大規模并行模擬程序,以此為基礎,開展了數值模擬驗證和確認、并行效率測試、大規模多級非定常算例計算等工作。數值計算結果表明:

圖20 前兩級葉片70%葉高焓值云圖Fig.20 Enthalpy contours of 70%span of first two balde row

1)全環非定常數值模擬方法得到的壓氣機整體氣動性能與試驗數據吻合良好。該方法可以消除周期性模型的限制,能夠保證物理量在葉輪機械轉靜葉片交界面處的連續,保證交界面處物理量的連續,削弱非物理熵增,捕捉到更精細的非定常流場結構,從而為更加深入研究揭示發動機內部復雜流動機理提供技術支撐。

2)本文所采用的METIS的網格分區方法以及MPI并行策略使程序具有良好的負載均衡和并行效率,同時具備針對不同種類葉輪機械模擬的適應性。

3)對億級規模網格的超大規模算例計算表明,本文所建立的計算方法和平臺對于超大規模復雜工程應用問題具備良好的可拓展性,能夠實現葉輪機械內部復雜的三維非定常流動細節模擬,可以滿足實際復雜工程問題的流動分析和葉輪機械的精細流動研究要求。

綜上,隨著未來高性能計算技術的不斷發展,本文所采用的全環數值模擬模型以及大規模并行程序策略將為揭示葉輪機械內部復雜三維非定常流動機理細節發揮更重要作用。實踐證明湍流模型、插值方法等對程序模擬精度有較大影響,下一步工作中將繼續針對此對程序改進提升。同時,將繼續借助強大的計算平臺,結合全環計算方法的優勢重點進行葉輪機械非定常流動機理、壓氣機流動失穩等方面的研究。