張 健唐 靜邱 名鄧有奇龔小權(quán)
(中國空氣動力研究與發(fā)展中心 計算空氣動力研究所,綿陽 621000)
葉輪機械內(nèi)部流動非常復(fù)雜,其本質(zhì)為三維黏性非定常流動。計算流體力學(xué)(Computational Fluid Dynamics,CFD)技術(shù)飛速發(fā)展,其具備計算耗時短,成本相對較低,并且可以多維度地模擬出葉輪機內(nèi)部復(fù)雜的非定常流動細節(jié)的優(yōu)點,因而越來越多地用于葉輪機械的氣動設(shè)計和分析當(dāng)中。
葉輪機內(nèi)部非定常流動的來源有多種,動靜葉片排的相對轉(zhuǎn)動、激波邊界層干擾、尾跡傳播、二次流等都是造成流動不穩(wěn)定的原因。通過數(shù)值模擬對上述非定?,F(xiàn)象進行研究的關(guān)鍵在于如何建立模擬轉(zhuǎn)靜葉片相對運動過程的模型。目前工業(yè)上一般采用混合平面(Mixing Plane)模型[1],該模型將轉(zhuǎn)/靜界面上下游相同展向高度處的通量通過周向平均后進行交換,將非定常計算簡化為對一個葉片流道的定常計算,雖然大大減小了計算量,然而卻無法捕捉到轉(zhuǎn)靜子之間相互干擾等非定常現(xiàn)象。采用葉片約化技術(shù)[2],將葉片在周向方向上按照一定比例進行幾何縮放,能夠使轉(zhuǎn)靜葉片數(shù)具有較大公約數(shù)進而約化為少量幾個葉片通道進行非定常計算,但是由于改變了實際的幾何尺寸,這種方法會存在較大誤差。He和Ning提出的非線性諧波法(Non-linear Harmonic Method,NLH)[3]可以看作是一種定常/非定?;旌戏椒?其基本假設(shè)是流場的主要擾動是由于葉片通過頻率(Blade Passing Frequence,BPF)引起的,從而將流場變量分解成為時間平均值和多個不同頻率諧波的擾動組成,BPF的整數(shù)倍數(shù)分別代表了不同諧波。為了得到最終的流場解,除了要求解時間平均方程外,還需要求解各個不同頻率的諧波方程。NLH模型優(yōu)點是在計算消耗可承受狀態(tài)下得到相對精確的非定常流場,但缺點是需要預(yù)先給定頻率,對于未知頻率流動和多頻率的組合模擬存在困難,特別是旋轉(zhuǎn)失速等非定常流動,NLH不可能解決。而且,隨著頻率數(shù)目增加,NLH的求解量也會劇增。在未來的航空發(fā)動機設(shè)計流程中,全環(huán)非定常數(shù)值模擬多級壓氣機/渦輪內(nèi)部的復(fù)雜三維流動必然會成為一種趨勢和常規(guī)手段。NASA支持開發(fā)的流場求解器TURBO[4]就是一個精確模擬三維非定常發(fā)動機流場的軟件。Chen等利用TURBO全環(huán)數(shù)值模擬了壓氣機失速開始狀態(tài)的流場,并研究了相關(guān)的失速控制技術(shù)[5]。北航邵飛等利用商用軟件CFX全環(huán)計算了某高壓渦輪非定常流動,捕捉到了流場中存在的低頻擾動[6]。
由于計算量巨大,展開葉輪機械全環(huán)非定常研究具有較大挑戰(zhàn)性。隨著高性能計算機的飛速發(fā)展,計算機硬件逐漸不再是限制CFD計算的主要因素。我國的神威太湖之光[7]、天河二號[8]等超級計算機已經(jīng)達到了每秒億億次量級的浮點數(shù)計算能力。數(shù)千萬甚至上億網(wǎng)格量的計算已不再是難題。美國NASA關(guān)于2030年CFD發(fā)展展望的報告[9]中,多次強調(diào)了高性能計算(HPC)在未來葉輪設(shè)計的重要作用。基于此,本文提出了一個大規(guī)模并行模擬程序設(shè)計策略,并針對葉輪機械的三維非定常流動進行了全環(huán)數(shù)值模擬,從而對更加深入研究發(fā)動機內(nèi)部復(fù)雜流動機理提供保障和基礎(chǔ)。
本文依托氣動中心千萬億次高性能計算集群,在自主研發(fā)的三維混合網(wǎng)格流場求解器MFlow[10]的基礎(chǔ)上,開發(fā)構(gòu)建了葉輪機定常模塊、非定常流動數(shù)值模擬模塊。滑移網(wǎng)格及通量守恒插值模塊的研發(fā)實現(xiàn)了航發(fā)葉輪機的全環(huán)非定常模擬。研究先通過一個典型單級壓氣機進行了程序的驗證與確認(rèn),然后利用1.5級渦輪算例進行了程序并行效率的測試,最后對一個3.5級壓氣機超大規(guī)模算例進行了嘗試,分析了動靜葉片相互干擾過程,評估了程序?qū)τ诖笠?guī)模復(fù)雜工程應(yīng)用的適用能力。
流場求解采用三維非定常雷諾平均Navier-Stokes方程(URANS)在旋轉(zhuǎn)坐標(biāo)系下的形式:
式中:Ω為控制體單元;S代表單元面。U為守恒變量,F I和F V分別是無黏通量和黏性通量,S T為科里奧利力和離心力引起的源項,具體形式如下:
式中:ρ為密度,u為相對速度矢量,E為總能,p為壓力,q為熱通量,τ為剪切應(yīng)力張量,ω為旋轉(zhuǎn)角速度,r為徑向矢量。
對于方程式(1),時間推進采用雙重時間步法,外迭代采用歐拉后插顯式推進,每一個時間步子迭代利用LU-SGS(Lower-Upper Symmetric Gauss-Seidel)方法[11]隱式迭代??臻g離散基于控制體的有限體積法,采用Roe迎風(fēng)格式[12]計算無黏通量,利用格林-高斯公式計算控制體單元的梯度,選擇Venkataknshnan限制器[13]將控制體的原始變量通過線性重構(gòu)插值到控制體單元面左右,達到更高的空間精度。采用中心格式計算黏性通量,湍流模型選取Spalart-Allmaras(SA)一方程模型[14]。
對于進口邊界,給定總溫T0,總壓p0以及來流速度方向。向外發(fā)出的黎曼不變量R-由第一個內(nèi)部邊界點外插計算得到,利用R-計算出絕對速度大小V,再根據(jù)速度方向?qū)分解為各個速度分量。
式中:v d代表邊界內(nèi)部第一個單元速度矢量,c d為當(dāng)?shù)芈曀?γ為比熱比,n為法向量,C p為定壓比熱。密度ρ和壓力p利用T0、p0、V之間的代數(shù)關(guān)系以及等熵關(guān)系計算得到。
對于內(nèi)流出口,給定靜壓,其他變量外插得到。同時引入徑向平衡方程式(5),通過給定展向高度r處的靜壓,其他位置的靜壓根據(jù)密度ρ和周向速度vθ積分得到。
對于黏性物面,為絕熱壁,采用無滑移邊界條件。對于轉(zhuǎn)靜交界面,在每一個物理時間步,上下游的守恒變量隨著網(wǎng)格滑移動態(tài)完成交換。動靜葉排之間利用滑移網(wǎng)格進行耦合,通過將對向網(wǎng)格的守恒變量插值到本方虛擬網(wǎng)格上完成邊界條件的更新,如圖1所示。
圖1 滑移網(wǎng)格技術(shù)Fig.1 Sliding interface technique
被插值單元的值通過包裹單元的節(jié)點值進行反距離權(quán)插值平均方法得到。如式(6)所示,ωi,j為距離權(quán)重,d i,j為兩點距離。
圖2給出了整個程序的流程圖,在每一個非定常時間步完成方程求解、網(wǎng)格旋轉(zhuǎn)、建立映射關(guān)系的迭代循環(huán),在每一個子迭代需要完成并行網(wǎng)格分區(qū)之間以及轉(zhuǎn)/靜界面之間的數(shù)據(jù)交換。
圖2 并行程序框架Fig.2 Parallel program frame
并行區(qū)域的劃分是在前處理過程中完成的,分區(qū)完成后,程序?qū)⒂涗浵赂鱾€分區(qū)的幾何信息,同時記錄下并行左右分區(qū)的映射關(guān)系。轉(zhuǎn)/靜界面的映射關(guān)系是在每一個非定常步求解開始之前完成。數(shù)據(jù)的并行傳值利用消息傳遞接口(Message Passing Interface,MPI)實現(xiàn)。利用MPI_Comm_Split方法將不同葉片排分區(qū)網(wǎng)格分到不同的子通信域,子通信域內(nèi)完成常規(guī)并行邊界的通量傳遞,子通信域之間完成滑移邊界的變量插值與傳遞。
并行網(wǎng)格分區(qū)利用圖切分軟件包METIS完成,METIS是由美國密西根大學(xué)G.Karypis和V.Kumar編寫的用于圖的分區(qū)和稀疏矩陣排序的串行包,提供k路多層劃分算法對混合網(wǎng)格進行分區(qū)[15]。給定包含n個節(jié)點V和m個邊E的圖G=(V,E),將V劃分到k個子集中V=[V1,V2,...,V k],使得|V i|≈n/k,且相連頂點屬于不同子集的邊的數(shù)量要最小化。
圖3給出了k路多層劃分過程,整個劃分過程分成coarsening,initial partitioning和uncoarsening三個部分。coarsening將圖的大小逐漸縮小,G1->G2->G3->G4。在G4階段執(zhí)行k路劃分,然后在uncoarsening階段將圖中的原始節(jié)點映射到G4劃分的子集中。由于本文采用格心格式的有限體積法,圖定義中的節(jié)點和邊分別指網(wǎng)格中的網(wǎng)格單元和網(wǎng)格面。采用k路多層圖劃分,保證了各個子區(qū)域的總網(wǎng)格數(shù)比較均衡的同時盡量減少了分區(qū)的邊切割,達到了負(fù)載均衡目的同時,最大限度減少了分區(qū)間的數(shù)據(jù)通信量。
圖3 METIS的k路多層圖劃分過程Fig.3 Multilevel k-way partitioning of METIS
本文并行算法基于分布式大規(guī)模并行計算機系統(tǒng),采用所有CPU都參與迭代計算的對等并行模式,MPI數(shù)據(jù)傳遞采用非阻塞通信,盡量加大計算和通信的重疊,每個分區(qū)網(wǎng)格只保存必需的局部網(wǎng)格數(shù)據(jù)和并行分區(qū)邊界上網(wǎng)格單元及頂點的對應(yīng)關(guān)系。多級葉輪網(wǎng)格內(nèi)部有兩種并行邊界,一種是同一葉片排內(nèi)部,在采用METIS并行網(wǎng)格分區(qū)時形成的交界面,另一種是不同葉片排之間的滑移網(wǎng)格交界面。對于前者,并行邊界面左右網(wǎng)格的對應(yīng)關(guān)系在前處理分區(qū)時就已經(jīng)確定,所以只需要計算一次并記錄下來,以空間換時間,提高并行效率。對于后者,每一非定常步轉(zhuǎn)子葉片網(wǎng)格將進行旋轉(zhuǎn),滑移邊界左右的對接信息發(fā)生更新,需要重新計算。插值網(wǎng)格之間的映射關(guān)系通過ADT(Alternating Digital Tree)方法[16]建立,該方法將空間單元進行有序排序后,構(gòu)建為二叉樹數(shù)據(jù),以便待插值點快速搜索查找[17],減少計算時間。
程序中,對于每一個分區(qū)所在的進程中,將該進程的并行邊界左右單元(物理單元以及虛擬單元)編號按照向相鄰分區(qū)發(fā)送和從相鄰分區(qū)接收的順序連續(xù)存儲在數(shù)組內(nèi)存中(圖4),交換數(shù)據(jù)時只需要將相應(yīng)數(shù)據(jù)按照單元編號存入或取出,并調(diào)用MPI數(shù)據(jù)傳遞接口完成。
圖4 交換數(shù)據(jù)存儲結(jié)構(gòu)Fig.4 Transfer data's storage structure
NASA Stage 35作為典型跨聲速壓氣機常用于壓氣機內(nèi)流數(shù)值模擬驗證。表1列出了Stage 35的具體的設(shè)計參數(shù)。Stage 35在設(shè)計轉(zhuǎn)速17188.7 rpm、流量20.2 kg/s下產(chǎn)生1.8的總壓比,Reid和Moore早期對其做過試驗,可以得到詳細的幾何參數(shù)、工作條件和試驗數(shù)據(jù)[18]。本文對Stage 35算例分別利用混合平面法進行了單通道定常計算[19]、利用非線性諧波方法進行了“準(zhǔn)”非定常計算、利用滑移面法進行了360°全環(huán)非定常計算,其中非線性諧波法是通過商用軟件Numeca完成的,因為只有一組轉(zhuǎn)靜干擾,所以周期性擾動個數(shù)設(shè)為1,每個擾動采用的諧波數(shù)為3個。最后將計算結(jié)果同試驗數(shù)據(jù)以及NASA TURBO軟件計算的360°全環(huán)非定常數(shù)據(jù)進行了對比分析。
表1 Stage 35設(shè)計參數(shù)Table 1 Design parameters of Stage 35
本次計算的網(wǎng)格如圖5所示,首先完成單通道的網(wǎng)格生成,第一層網(wǎng)格間距3×10-6m,y+取1,同時考慮了葉尖間隙的影響,間隙展向分布了17個網(wǎng)格點。然后旋轉(zhuǎn)復(fù)制得到Stage 35全環(huán)網(wǎng)格??偩W(wǎng)格量大概在3300萬左右。
本次的計算在中國空氣動力研究與發(fā)展中心的千萬億次集群上完成。并行計算采用512個核,集群芯片的處理器為國產(chǎn)的飛騰FT-1500A處理器,主頻2 GHz。每個時間步葉片旋轉(zhuǎn)1°,總共完整計算3圈。圖6給出了計算穩(wěn)定后某一時刻50%葉高全環(huán)剖面的壓力分布云圖。
為了獲取壓氣機的整體性能并且評估數(shù)值模擬的精度,圖7和圖8分別給出了壓氣機在100%轉(zhuǎn)速下不同方法和模型計算得到的總壓比和絕熱效率隨流量變化的特性曲線,同時給出試驗數(shù)據(jù)進行對比。從圖看出,數(shù)值模擬得到的壓氣機整體性能和試驗數(shù)據(jù)吻合良好。相比較而言,本文的全環(huán)非定常計算的總壓比相對于文獻[20]中TURBO軟件全環(huán)的計算結(jié)果與試驗數(shù)據(jù)更加接近,尤其是在接近阻塞(試驗中的4004工況)時。本文計算得到最大流量與Reid和Moore試驗數(shù)據(jù)非常接近,表2給出了該狀態(tài)下的總壓比、總溫比和絕熱效率與試驗及TURBO計算結(jié)果的對比??梢钥吹組Flow的計算結(jié)果與試驗的誤差基本都在Reid和Moore提出的試驗不確定度范圍之內(nèi)。文獻中沒有給出TURBO計算的絕熱效率特性曲線,相比于混合平面法和非線性諧波法,MFlow全環(huán)非定常計算的絕熱效率與試驗更加接近。在小流量工況下計算得到的壓比相對混合面法和NLH低,其原因可能上下游網(wǎng)格信息傳遞時由于插值精度不夠,導(dǎo)致數(shù)值耗散大[21]。從整體性能來看,本文程序的計算結(jié)果具有較高精度,滿足工程應(yīng)用需求。
圖5 Stage 35計算網(wǎng)格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接近阻塞狀態(tài)氣動特性對比Table 2 Aerodynamic performance of Stage 35 near choke
圖9和圖10分別給出了利用不同模型計算得到的50%葉高切面壓力分布以及熵值云圖。通過對比可以看到不同模型計算得到的波系結(jié)構(gòu)基本一致,但是混合平面法由于在轉(zhuǎn)/靜交界面進行了周向平均,所以壓力和熵值并不連續(xù),尤其是熵增有一個很明顯的階躍。NLH方法的結(jié)果與非定常模擬的結(jié)果更加接近,能夠捕捉到轉(zhuǎn)/靜干擾現(xiàn)象。但是仔細觀察轉(zhuǎn)子產(chǎn)生的熵值進入靜子區(qū)域,在幾條高亮的尾跡中間還有兩條強度較小的尾跡。這是因為NLH方法是通過有限的幾個頻率諧波方程的定常解合成一個近似的非定常解,不同頻率會帶來不同波長和幅值的正弦振蕩,只有在諧波數(shù)足夠多的情況下才能完全消除擾動[22]。相比之下全環(huán)模擬的結(jié)果與實際物理現(xiàn)象更加接近。
圖9 不同方法計算壓力云圖比較Fig.9 Static pressure comparison of different methods
圖10 不同方法計算熵值云圖比較Fig.10 Entropy comparison of different methods
為研究相鄰葉片排之間的周期性干擾過程,在計算過程中對靜子尾跡接近出口某處設(shè)置一個壓力探測點,壓力隨時間周期變化結(jié)果如圖11所示,可以看到監(jiān)測點的靜壓在一個周期后呈現(xiàn)出良好的周期性,說明非定常計算已經(jīng)達到穩(wěn)定。對其進行頻譜分析結(jié)果如圖12所示,可以看到計算很好的捕捉到了葉片通過頻率(BPF)及其諧波。
圖11 探測點壓力變化曲線Fig.11 Static pressure at the probe
圖12 探測點壓力頻譜分析結(jié)果Fig.12 Pressure spectrum
為測試程序的并行計算效率,同時進一步檢驗程序?qū)Σ煌N類葉輪機械模擬的適應(yīng)性,本節(jié)采用了一個稍大網(wǎng)格規(guī)模的1.5級渦輪算例進行了計算和分析。其基本幾何構(gòu)型如圖13所示,全環(huán)轉(zhuǎn)子葉片共41個,轉(zhuǎn)速3500 rpm,靜子葉片36個[23]。圖14給出了計算所采用的其中一個葉片通道網(wǎng)格,整個3排葉片全環(huán)網(wǎng)格量共計約4000萬。
圖13 Aachen 1.5級渦輪幾何參數(shù)Fig.13 Geometry parameters of Aachen 1.5 turbine stage
圖14 Aachen 1.5級渦輪計算網(wǎng)格Fig.14 Computational grid of Aachen 1.5 turbine stage
計算不考慮渦輪的傳熱過程,每一個非定常步轉(zhuǎn)子大約轉(zhuǎn)動1度,每一個非定常步下子迭代200次。對同一狀態(tài)分別利用180、360、720以及1440個CPU核心數(shù)分別計算,從而進行并行效率的測試。圖15給出了50%展向環(huán)面的熵值云圖計算結(jié)果,反映了非定常流動損失的輸運過程。第一排靜葉產(chǎn)生熵的尾跡,遇到轉(zhuǎn)子葉片前緣被切斷后轉(zhuǎn)向進入轉(zhuǎn)子葉片通道并向下游傳輸,再與轉(zhuǎn)子葉片尾跡產(chǎn)生的熵混合,一起進入下游的第二排靜葉通道混合。從熵增的流動過程和在轉(zhuǎn)/靜界面的連續(xù)性來看,計算結(jié)果合理。圖16是利用不同CPU核數(shù)進行了一定步數(shù)非定常計算后殘差隨時間收斂曲線對比,可以明顯看到隨核數(shù)成倍增加,收斂時間幾乎成倍減少。圖17進一步給出了加速比曲線,可以看到720核時并行效率有90%,1440核并行效率依然能夠接近80%。利用1440核完成本算例一個狀態(tài)的計算,轉(zhuǎn)子完整旋轉(zhuǎn)一周計算時間不到24小時,顯示出程序具有較高的并行效率和加速比。
圖15 50%葉高環(huán)面熵值云圖Fig.15 Entropy contours at 50%span annulus
圖16 不同核數(shù)計算殘差隨時間收斂曲線Fig.16 Residual convergence using different number of cores
圖17 并行效率測試結(jié)果Fig.17 Parallel efficiency
在完成程序的驗證與確認(rèn)及性能分析后,利用程序?qū)σ粋€超大規(guī)模算例進行了嘗試。該算例是某壓氣機的前3.5級,設(shè)計轉(zhuǎn)速15000 rpm。圖18給出了本次計算所用的網(wǎng)格,網(wǎng)格總量接近2億,共分4096個核進行計算。
圖18 3.5級壓氣機全環(huán)網(wǎng)格Fig.18 Full annulus grid of 3.5 stage compressor
表3給出了在接近設(shè)計點的一個狀態(tài)下,分別利用混合平面法和全環(huán)非定常方法計算得到的壓氣機氣動性能對比??梢钥吹?由于混合平面法方法模型本身的缺陷,進出口流量誤差隨著級數(shù)的增多累積明顯,達到2.5%,而全環(huán)非定常的差別只有0.7%,其中可能是插值的精度損失造成的。除流量外,定常模型計算得到的總壓比、絕熱效率與非定常比差別較明顯,說明當(dāng)壓氣機級數(shù)越多,定常模型計算結(jié)果的可信度就越低。
表3 3.5級壓氣機設(shè)計點計算結(jié)果Table3 Aerodynamic performance of 3.5 stage compressor at design point
圖19給出了計算得到的整機壓力分布云圖,圖20給出了三個不同時刻在前兩級葉片排在70%葉高切面的焓值局部放大云圖,可以幫助分析壓氣機內(nèi)部復(fù)雜的動靜干擾過程??梢钥吹?轉(zhuǎn)子葉柵外伸激波會射入上游靜子葉柵通道,外伸激波在壓力面干擾強、傳播速度快、并形成反射,從而使得靜葉壓力面的附面層快速增厚,靜葉表面同時也會受到上游轉(zhuǎn)子渦脫落尾跡的影響,使得靜葉表面壓力存在劇烈的非定常波動。受上下游靜子尾跡及壓力波傳播干擾,轉(zhuǎn)子通道內(nèi)的激波形態(tài)也會不斷發(fā)生變化。說明全環(huán)模擬能夠揭示更多葉輪機械內(nèi)部復(fù)雜的三維非定常流動細節(jié)。
圖19 3.5級壓氣機壓力云圖Fig.19 Pressure contours of 3.5 stage compressor
本文研究開發(fā)了基于滑移網(wǎng)格技術(shù)、METIS網(wǎng)格分區(qū)技術(shù)和并行邊界處理及虛擬單元技術(shù)等的多級葉輪機械全環(huán)非定常大規(guī)模并行模擬程序,以此為基礎(chǔ),開展了數(shù)值模擬驗證和確認(rèn)、并行效率測試、大規(guī)模多級非定常算例計算等工作。數(shù)值計算結(jié)果表明:
圖20 前兩級葉片70%葉高焓值云圖Fig.20 Enthalpy contours of 70%span of first two balde row
1)全環(huán)非定常數(shù)值模擬方法得到的壓氣機整體氣動性能與試驗數(shù)據(jù)吻合良好。該方法可以消除周期性模型的限制,能夠保證物理量在葉輪機械轉(zhuǎn)靜葉片交界面處的連續(xù),保證交界面處物理量的連續(xù),削弱非物理熵增,捕捉到更精細的非定常流場結(jié)構(gòu),從而為更加深入研究揭示發(fā)動機內(nèi)部復(fù)雜流動機理提供技術(shù)支撐。
2)本文所采用的METIS的網(wǎng)格分區(qū)方法以及MPI并行策略使程序具有良好的負(fù)載均衡和并行效率,同時具備針對不同種類葉輪機械模擬的適應(yīng)性。
3)對億級規(guī)模網(wǎng)格的超大規(guī)模算例計算表明,本文所建立的計算方法和平臺對于超大規(guī)模復(fù)雜工程應(yīng)用問題具備良好的可拓展性,能夠?qū)崿F(xiàn)葉輪機械內(nèi)部復(fù)雜的三維非定常流動細節(jié)模擬,可以滿足實際復(fù)雜工程問題的流動分析和葉輪機械的精細流動研究要求。
綜上,隨著未來高性能計算技術(shù)的不斷發(fā)展,本文所采用的全環(huán)數(shù)值模擬模型以及大規(guī)模并行程序策略將為揭示葉輪機械內(nèi)部復(fù)雜三維非定常流動機理細節(jié)發(fā)揮更重要作用。實踐證明湍流模型、插值方法等對程序模擬精度有較大影響,下一步工作中將繼續(xù)針對此對程序改進提升。同時,將繼續(xù)借助強大的計算平臺,結(jié)合全環(huán)計算方法的優(yōu)勢重點進行葉輪機械非定常流動機理、壓氣機流動失穩(wěn)等方面的研究。