波轉(zhuǎn)子非定常泄漏流動(dòng)機(jī)理

劉琛源， 劉火星

北京航空航天大學(xué) 能源與動(dòng)力工程學(xué)院， 北京 100083

波轉(zhuǎn)子非定常泄漏流動(dòng)機(jī)理

劉琛源， 劉火星*

北京航空航天大學(xué) 能源與動(dòng)力工程學(xué)院， 北京 100083

針對(duì)制約波轉(zhuǎn)子性能的泄漏問題，通過提取波轉(zhuǎn)子中與泄漏相關(guān)的主要流動(dòng)現(xiàn)象并建立波轉(zhuǎn)子單通道泄漏模型，對(duì)非定常泄漏流動(dòng)機(jī)制進(jìn)行了詳細(xì)的數(shù)值研究，對(duì)本文給出的非定常泄漏損失預(yù)測(cè)模型進(jìn)行了數(shù)值驗(yàn)證。結(jié)果表明：轉(zhuǎn)子通道中存在不同程度的壓力波動(dòng)，波動(dòng)幅值與間隙寬度有關(guān)；連續(xù)反射膨脹波、周期性出現(xiàn)的弓形激波及其反射激波是壓力波動(dòng)的根本原因；間隙內(nèi)部泄漏過程存在3個(gè)主要的流動(dòng)階段；泄漏過程中通道激波傳播速度不變、波后時(shí)均壓力不變；在一定間隙寬度范圍內(nèi)，激波馬赫數(shù)、激波靜增壓比與無量綱間隙寬度均呈線性關(guān)系，當(dāng)間隙寬度從0增大到0.08時(shí)，激波馬赫數(shù)衰減7.3%，激波靜增壓比衰減10.1%；泄漏流動(dòng)通過泄漏產(chǎn)生的主膨脹波對(duì)激波傳播過程施加影響，通道激波衰減本質(zhì)上是理想激波與主膨脹波疊加效應(yīng)的結(jié)果；泄漏損失預(yù)測(cè)模型與數(shù)值結(jié)果吻合良好。

燃?xì)廨啓C(jī)； 波轉(zhuǎn)子； 單通道模型； 非定常泄漏流動(dòng)； 激波

20世紀(jì)初人們便已認(rèn)識(shí)到非定常流動(dòng)巨大的應(yīng)用前景，但是非定常流動(dòng)自身的復(fù)雜性在很長(zhǎng)一段時(shí)期內(nèi)阻礙了波轉(zhuǎn)子技術(shù)的發(fā)展和應(yīng)用[1]。近年來，隨著氣動(dòng)設(shè)計(jì)的不斷完善以及材料性能的限制，傳統(tǒng)燃?xì)廨啓C(jī)逐漸達(dá)到性能極限[2]。波轉(zhuǎn)子通過運(yùn)動(dòng)激波和膨脹波直接在不同狀態(tài)的氣體間交換能量，因此具備巨大的增壓潛力以及相對(duì)較高的增壓效率[3]；另外，由于能量交換過程不需要壓氣機(jī)和渦輪葉片，結(jié)構(gòu)相對(duì)簡(jiǎn)單，因此能夠有效降低發(fā)動(dòng)機(jī)結(jié)構(gòu)重量及制造成本[4]。

基于波轉(zhuǎn)子的工作特點(diǎn)，相比傳統(tǒng)燃?xì)廨啓C(jī)它還具有以下幾個(gè)優(yōu)勢(shì)：① 轉(zhuǎn)子轉(zhuǎn)速相對(duì)較低，因此材料應(yīng)力??；② 響應(yīng)迅速，響應(yīng)時(shí)間與壓力波的特征時(shí)間尺度為相同量級(jí)；③ 轉(zhuǎn)子通道內(nèi)工質(zhì)流速低，僅為典型的葉輪機(jī)械工質(zhì)流速的1/3左右[5]，因此通道壁面更耐腐蝕；④ 空氣和燃?xì)饨惶嫱ㄟ^轉(zhuǎn)子通道，具有自冷卻效應(yīng)，可以承受更高的燃?xì)鉁囟取?/p>

波轉(zhuǎn)子技術(shù)具有大幅提高現(xiàn)有發(fā)動(dòng)機(jī)性能的潛力，但是仍然存在一些技術(shù)障礙亟需克服，主要體現(xiàn)在3個(gè)方面：1)氣動(dòng)設(shè)計(jì)，即如何根據(jù)給定的應(yīng)用需求確定最佳設(shè)計(jì)參數(shù)、同時(shí)解決非設(shè)計(jì)點(diǎn)性能問題[4]；2)結(jié)構(gòu)設(shè)計(jì)，主要解決轉(zhuǎn)靜部件間的熱膨脹及密封問題[4]；3)損失機(jī)理分析及損失評(píng)估[6]，波轉(zhuǎn)子內(nèi)部損失主要包括：① 非定常波系損失，包括激波固有的熵增損失以及膨脹波后流動(dòng)不均勻帶來的摻混損失；② 轉(zhuǎn)子通道非瞬時(shí)開閉損失，通道非瞬時(shí)開閉使觸發(fā)的波系偏離理想波系，進(jìn)一步形成包含激波、旋渦在內(nèi)的復(fù)雜流動(dòng)結(jié)構(gòu)，流動(dòng)的紊亂會(huì)對(duì)波轉(zhuǎn)子性能產(chǎn)生不利影響；③ 非軸向進(jìn)氣損失，與通道非瞬時(shí)開閉具有類似的損失機(jī)制；④ 摩擦損失，轉(zhuǎn)子通道一般具有較大的相對(duì)表面積，壁面摩擦不可忽略；⑤ 泄漏損失，波系在端壁反射引起通道、端口、轉(zhuǎn)子容腔之間的壓力差，進(jìn)而引起泄漏。

波轉(zhuǎn)子氣動(dòng)設(shè)計(jì)方法研究較為完善，但是針對(duì)流動(dòng)損失機(jī)理的研究相對(duì)較少，這在一定程度上限制了波轉(zhuǎn)子技術(shù)的發(fā)展和應(yīng)用。Wilson[6]針對(duì)商用波轉(zhuǎn)子Comprex研究了各類損失對(duì)效率的影響，Wilson指出，泄漏損失在所有損失中占主導(dǎo)地位，特別對(duì)于小型波轉(zhuǎn)子，泄漏對(duì)性能的影響更為突出。Wilson[6]在三端口和四端口波轉(zhuǎn)子實(shí)驗(yàn)中發(fā)現(xiàn)泄漏對(duì)波轉(zhuǎn)子性能下降影響很大。

目前在世界范圍內(nèi)，針對(duì)泄漏流動(dòng)主要通過泄漏模型計(jì)入其對(duì)性能的影響，包括廣泛應(yīng)用于波轉(zhuǎn)子準(zhǔn)一維設(shè)計(jì)程序的單容腔模型[7]以及Akbari和Nalim[8]提出的分段泄漏模型。上述2種模型均未考慮具體流動(dòng)過程，特別是泄漏流動(dòng)和主流的相互作用。由于運(yùn)動(dòng)激波在波轉(zhuǎn)子的工作機(jī)制中起著重要作用，因此了解泄漏對(duì)激波傳播特性以及波后流場(chǎng)的影響非常重要。本文采用數(shù)值研究手段詳細(xì)分析了波轉(zhuǎn)子的泄漏流動(dòng)及損失機(jī)制，并提出了泄漏損失預(yù)測(cè)模型。

1 波轉(zhuǎn)子工作過程及泄漏流動(dòng)特征

圖1給出了一種典型的用于航空燃?xì)廨啓C(jī)循環(huán)的四端口通流外燃波轉(zhuǎn)子結(jié)構(gòu)[2]及實(shí)驗(yàn)裝置[9]，它由旋轉(zhuǎn)鼓筒及鼓筒兩側(cè)的靜止端板構(gòu)成，轉(zhuǎn)子通道位于旋轉(zhuǎn)鼓筒上。鼓筒高速旋轉(zhuǎn)，轉(zhuǎn)子通道與各個(gè)端口周期性地連通、閉合，從而在通道內(nèi)部觸發(fā)運(yùn)動(dòng)的激波和膨脹波。壓氣機(jī)出口空氣通過低壓空氣(LPA)端口進(jìn)入轉(zhuǎn)子通道，燃燒室出口燃?xì)鈩t通過高壓燃?xì)?HPG)端口進(jìn)入轉(zhuǎn)子通道，燃?xì)馀c空氣接觸后觸發(fā)運(yùn)動(dòng)激波，空氣經(jīng)激波增壓后通過高壓空氣(HPA)端口進(jìn)入燃燒室，最后燃?xì)馀蛎浐蠼?jīng)低壓燃?xì)?LPG)端口進(jìn)入渦輪，從而完成整個(gè)熱力循環(huán)過程。

通過波系的壓縮和膨脹作用，流入流出4個(gè)端口的氣體快速完成能量交換，最終使得來自壓氣機(jī)的空氣在更高的壓力和溫度下進(jìn)入燃燒室燃燒，從而提高整個(gè)燃?xì)廨啓C(jī)的熱力循環(huán)效率，同時(shí)保證膨脹后進(jìn)入渦輪的燃?xì)鈮毫Ω哂趬簹鈾C(jī)出口的空氣壓力，從而提高渦輪輸出功率。

波轉(zhuǎn)子發(fā)生泄漏的根本原因是轉(zhuǎn)子通道依次通過各個(gè)端口時(shí)，在每個(gè)通道末端觸發(fā)的激波或膨脹波之間存在一定時(shí)滯，因此在通道末端同一時(shí)刻不同通道內(nèi)的氣體壓力可能不同，高壓通道氣體通過轉(zhuǎn)靜間隙進(jìn)入低壓通道引起周向泄漏，通道與容腔之間的壓力差則引起徑向泄漏。

圖2給出了波轉(zhuǎn)子高壓循環(huán)的波系結(jié)構(gòu)。如圖2所示，當(dāng)入射激波(Incident shock wave)到達(dá)通道右端時(shí)，立即產(chǎn)生一道向左端傳播的反射激波(Reflected shock wave)，反射激波波后氣體此時(shí)位于高壓區(qū)域2中，經(jīng)過兩道激波的增壓作用，2區(qū)氣體壓力明顯高于1區(qū)以及轉(zhuǎn)子容腔壓力，2區(qū)通道內(nèi)的高壓氣體經(jīng)轉(zhuǎn)靜間隙分別泄漏到容腔以及1區(qū)通道內(nèi)部。

圖1 典型的四端口通流波轉(zhuǎn)子結(jié)構(gòu)及實(shí)驗(yàn)裝置[2,9]Fig.1 Typical structure of four-port through-flow wave rotor and experimental apparatus[2,9]

圖2 波轉(zhuǎn)子高壓循環(huán)波系結(jié)構(gòu) Fig.2 Wave diagram of high-pressure part of wave rotor

2 泄漏分析數(shù)值模型

數(shù)值模型通過構(gòu)造黎曼間斷形成通道內(nèi)的運(yùn)動(dòng)激波[10]。由于波轉(zhuǎn)子內(nèi)部存在通道非瞬時(shí)開閉、非軸向進(jìn)氣等多種損失，為了深入分析泄漏損失機(jī)理，通過建立二維單通道泄漏模型，隔離其他損失因素對(duì)泄漏流動(dòng)的影響，圖3給出了該模型的幾何示意圖。

激波強(qiáng)度、波前氣體狀態(tài)、氣體物性以及通道截面尺寸等均對(duì)激波傳播特性存在影響。Glass 和Patterson[11]通過實(shí)驗(yàn)研究發(fā)現(xiàn)，激波靜壓比小于3時(shí)，在一定距離內(nèi)，激波傳播速度不變；靜壓比大于3時(shí)，最大波速低于理論值，并且在傳播過程中波速逐漸減小。另外，隨著通道截面尺寸減小，氣體黏性效應(yīng)、稀薄效應(yīng)和壁面?zhèn)鳠嵝?yīng)越來越顯著，從而導(dǎo)致激波在傳播過程中不斷衰減[12]。

(1)

(2)

(3)

式中：pi為壓力；ai為聲速；γi為比熱比；Ma為入射激波馬赫數(shù)；i=1,2,4,5。

為了準(zhǔn)確反映激波在波轉(zhuǎn)子通道中的傳播特性，需要選擇合適的氣動(dòng)參數(shù)和模型尺寸。四端口通流波轉(zhuǎn)子增壓比通常在2.5～3.5之間，此時(shí)通道進(jìn)口觸發(fā)的激波靜壓比小于2。本文給定泄漏模型低壓段初始?jí)毫?01 325 Pa、高壓段初始?jí)毫?05 300 Pa，初始溫度均為298.15 K，初始時(shí)刻氣體均靜止，工質(zhì)為采用理想氣體模型的空氣。設(shè)初始時(shí)刻低壓段為1區(qū)、高壓段為4區(qū)，入射激波形成后波后為2區(qū)，反射激波波后為5區(qū)，根據(jù)式(1)～式(3)，泄漏流量為零時(shí)的理論參數(shù)如下：入射激波壓比p2/p1=1.93；反射激波壓比p5/p2=1.82；總增壓比p5/p1=3.51。數(shù)值模型幾何參數(shù)選擇參考表1[13]，具體尺寸在圖3中給出。

圖3 波轉(zhuǎn)子二維單通道泄漏模型 Fig.3 Two-dimensional model for single passage leakage of wave rotor

表1 典型波轉(zhuǎn)子構(gòu)型的主要幾何參數(shù)[13]Table 1 Main geometric parameters of typical configurations of wave rotor[13]

ConfigurationNASA(3?Port)NASA&AllisonKentfieldGEComprex(ABB)TokyoUniversityLengthL/mm230，46015228030093．2186Meanradiusr/mm15081．5102584860CellwidthW/mm6．35，12．78．7516．81098，16CellheightH/mm10．22255．923．11016Clearanceδ/mm0．127?0．3810．130．180．640．150．5?1．0

3 數(shù)值模型求解方法及驗(yàn)證

3.1 數(shù)值模型

數(shù)值模型各算例泄漏間隙寬度最小值為0.1 mm，間隙內(nèi)部流動(dòng)Kn數(shù)約為5.6×10-4，因此所有算例均可采用基于連續(xù)介質(zhì)假設(shè)的Navier-Stokes方程進(jìn)行數(shù)值求解[14]。

選擇FLUENT 12.1密度基隱式求解器，流動(dòng)控制方程對(duì)流項(xiàng)離散采用AUSM+(Advection Upstream Splitting Method)格式[15]，對(duì)流通量計(jì)算進(jìn)行三階MUSCL(Monotonic Upwind Scheme for Conservation Laws)插值，黏性項(xiàng)離散采用中心差分格式，時(shí)間項(xiàng)離散采用二階隱式格式，湍流模擬選擇Spalart-Allmaras一方程模型。如圖3所示，Symmetry采用對(duì)稱面邊界條件[16]，其他模型邊界均選擇無滑移絕熱壁面。所有算例物理時(shí)間步最大值為10-8s，當(dāng)激波傳播到通道末端時(shí)調(diào)整為5×10-10s；虛擬時(shí)間步設(shè)置保證方程殘差不再減小；數(shù)值模擬總時(shí)間為0.42 ms。計(jì)算域離散采用四邊形結(jié)構(gòu)化網(wǎng)格，距壁面第一層網(wǎng)格高度y+最大值為0.001 mm，從而保證所有算例y+值小于1；平行壁面方向網(wǎng)格最大間距為0.1 mm；在間隙內(nèi)部對(duì)網(wǎng)格進(jìn)行局部加密，垂直壁面方向最大網(wǎng)格間距為0.01 mm；網(wǎng)格總數(shù)為22萬，經(jīng)網(wǎng)格無關(guān)性驗(yàn)證，證明上述網(wǎng)格可以很好地捕捉激波。工質(zhì)采用理想氣體模型，空氣動(dòng)力黏性系數(shù)根據(jù)薩瑟蘭定律確定，定壓比熱計(jì)算采用分段多項(xiàng)式法。

3.2 數(shù)值方法驗(yàn)證

針對(duì)本文采用的數(shù)值方法，需要重點(diǎn)考察以下2點(diǎn)：① 存在黏性效應(yīng)的情況下求解器對(duì)一維運(yùn)動(dòng)激波的捕捉能力；② 數(shù)值格式對(duì)非定常流場(chǎng)中復(fù)雜二維波系及其相互作用的模擬能力。

為了考察數(shù)值格式對(duì)一維運(yùn)動(dòng)激波的模擬精度，選擇日本東北大學(xué)無膜激波管實(shí)驗(yàn)對(duì)數(shù)值方法進(jìn)行校驗(yàn)。圖4給出了實(shí)驗(yàn)裝置[17]和相應(yīng)的數(shù)值模型，其中被驅(qū)動(dòng)段初始?jí)毫?0 000 Pa，入射激波馬赫數(shù)為1.2，計(jì)算域初始溫度為297 K，入射激波進(jìn)入微尺度通道后，5個(gè)等間距布置的壓力測(cè)點(diǎn)(P1～P5)依次記錄激波通過后流場(chǎng)的壓力變化，圖5給出了實(shí)驗(yàn)及數(shù)值結(jié)果，p為壓力，采用被驅(qū)動(dòng)段初始?jí)毫M(jìn)行無量綱化，t為時(shí)間。對(duì)于激波強(qiáng)度、波速以及激波的沿程衰減，數(shù)值和實(shí)驗(yàn)結(jié)果吻合良好。

為了驗(yàn)證數(shù)值格式對(duì)二維非定常流動(dòng)的模擬能力，考察入射激波在泄漏間隙內(nèi)的繞射過程。激波繞射是一種典型的自相似現(xiàn)象[18]，因此可以將激波在間隙內(nèi)的繞射流場(chǎng)及無量綱流動(dòng)參數(shù)與激波繞射實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，校驗(yàn)本文數(shù)值方法。數(shù)值校驗(yàn)共考察了入射激波馬赫數(shù)Ma=1.5和Ma=2.0兩個(gè)算例。

圖4 無膜激波管實(shí)驗(yàn)裝置簡(jiǎn)圖及數(shù)值模型 Fig.4 Experimental setup of diaphragmless shock tube and numerical model

圖5 激波管各測(cè)點(diǎn)壓力隨時(shí)間變化 Fig.5 Variation of pressure with time at different test locations in shock tube

(1) 流動(dòng)結(jié)構(gòu)驗(yàn)證

對(duì)于Ma=1.5算例，圖6(a)給出的泄漏間隙激波繞射密度等值線圖和圖6(b)[19]給出的實(shí)驗(yàn)結(jié)果均清晰地顯示出轉(zhuǎn)角處局部膨脹波和剪切層上緣激波構(gòu)成的的λ波結(jié)構(gòu)；對(duì)比圖6(c)和圖6(d)[19]可以看到，數(shù)值模擬較好地表現(xiàn)了繞射激波波后接觸面環(huán)繞旋渦并與剪切層前緣下表面融合的現(xiàn)象。

隨著入射激波馬赫數(shù)提高到2.0，繞射激波波后流動(dòng)結(jié)構(gòu)發(fā)生局部變化，轉(zhuǎn)角處λ波的膨脹波分支演變?yōu)榕c普朗特邁耶流類似的中心膨脹波，λ波的激波分支沿剪切層上表面向后移動(dòng)到剪切層尾部，同時(shí)強(qiáng)度大大增強(qiáng)，形成一道S形激波；圖7(a)的數(shù)值結(jié)果和圖7(d)[20]的實(shí)驗(yàn)結(jié)果均可觀察到典型的S形激波。

Skews[18]指出，對(duì)于Ma=2.0的入射激波，剪切層形狀基本保持為直線，僅在其尾部向旋渦一側(cè)彎曲，對(duì)于繞射激波波后間斷面，與低馬赫數(shù)入射激波情況相同的是間斷面均在靠近旋渦時(shí)環(huán)繞旋渦，但在接近剪切層時(shí)，間斷面突然轉(zhuǎn)向并以直線形式連接轉(zhuǎn)角頂點(diǎn)。對(duì)比圖6(c)和圖7(c)給出的流場(chǎng)熵增分布，可以觀察到入射激波馬赫數(shù)變化時(shí)剪切層和間斷面形狀的改變；另外，Ma=2.0時(shí)，從90° 轉(zhuǎn)角頂點(diǎn)處沿徑向形成一道新的間斷面，即中心膨脹波波尾，該間斷面與S形激波上緣相交。

圖6 Ma=1.5入射激波直角繞射數(shù)值和實(shí)驗(yàn)結(jié)果對(duì)比 Fig.6 Comparison between numerical and experimental results for incident shock of Mach number 1.5 diffracting at corner of ninety degrees

圖7 Ma=2.0入射激波直角繞射數(shù)值和實(shí)驗(yàn)結(jié)果對(duì)比 Fig.7 Comparison between numerical and experimental results for incident shock of Mach number 2.0 diffracting at corner of ninety degrees

(2) 流場(chǎng)參數(shù)驗(yàn)證

圖8為激波初始折轉(zhuǎn)角和剪切層角度示意圖。圖中：EW為膨脹波波頭；IS為入射激波；SL為剪切層；V為旋渦；W為壁面?？疾烊肷浼げǔ跏颊坜D(zhuǎn)角m0、剪切層角度e0、壁面激波馬赫數(shù)Maw、旋渦移動(dòng)速度u隨入射激波馬赫數(shù)變化，評(píng)估數(shù)值方法對(duì)激波膨脹波相互作用、氣流分離與剪切、旋渦運(yùn)動(dòng)等非定常流動(dòng)現(xiàn)象的模擬精度。圖9給出了Ma=1.5、Ma=2.0算例數(shù)值結(jié)果和相應(yīng)的實(shí)驗(yàn)結(jié)果，除Ma=2.0算例剪切層角度預(yù)測(cè)偏差略大外，其他參數(shù)均吻合較好。

通過詳細(xì)的數(shù)值校驗(yàn)，發(fā)現(xiàn)本文采用的數(shù)值方法可以準(zhǔn)確模擬激波傳播、低雷諾數(shù)下的激波衰減以及激波、膨脹波相互作用、旋渦運(yùn)動(dòng)、滑移間斷等非定常流動(dòng)現(xiàn)象。

圖8 激波初始折轉(zhuǎn)角和剪切層角度示意圖 Fig.8 Schematic diagram of shock diffraction angle and shear layer angle

圖9 典型流場(chǎng)參數(shù)隨入射激波馬赫數(shù)變化數(shù)值和實(shí)驗(yàn)結(jié)果對(duì)比 Fig.9 Comparison between numerical and experimental results for change of typical flow parameters with incident shock Mach number

4 計(jì)算結(jié)果分析

4.1 非定常泄漏流動(dòng)特征

為了深入分析非定常泄漏流動(dòng)的發(fā)生機(jī)制及其對(duì)通道內(nèi)激波傳播過程的影響，先考慮不含通道激波的泄漏過程。

(1) 不含通道激波的情況

選擇泄漏間隙0.8 mm算例進(jìn)行分析。圖10 給出了泄漏模型初始流場(chǎng)，主通道內(nèi)氣體狀態(tài)與零間隙寬度下反射激波波后狀態(tài)相同。在大于零時(shí)刻瞬間，一道激波在間隙入口形成并通過泄漏間隙進(jìn)入副通道，膨脹波形成后在主通道(以下簡(jiǎn)稱通道)中傳播。

圖11給出了不同時(shí)刻通道內(nèi)壓力分布，壓力數(shù)據(jù)采用通道初始?jí)毫M(jìn)行無量綱化。觀察圖11 發(fā)現(xiàn)，在壁面和通道中心截面(Symmetry)之間存在連續(xù)反射的二維膨脹波和周期性出現(xiàn)的二維激波及其反射激波。

圖10 不含通道激波時(shí)泄漏模型初始流場(chǎng) Fig.10 Initial flow field of leakage model without channel shock wave

觀察圖11(a)，泄漏流動(dòng)開始瞬間，在間隙入口處形成了一道扇形膨脹波，由于在膨脹波傳播方向流道幾何非對(duì)稱，膨脹波寬度沿波面出現(xiàn)不均勻性，引起波后氣流折轉(zhuǎn)，因此膨脹波后又形成一道弓形激波；膨脹波到達(dá)通道中心截面時(shí)，與另一側(cè)間隙處形成的膨脹波相交并相互穿過(可以理解為波在虛擬固壁上的反射)，在0.014 ms時(shí)可以清晰地觀察到對(duì)稱面附近的反射膨脹波，緊接著，激波到達(dá)對(duì)稱面并發(fā)生反射，在0.024 ms時(shí)可以清晰地觀察到反射激波。

隨著第一道反射膨脹波向壁面方向傳播，第一道弓形激波傳播過程中形成的波后高壓區(qū)域逐漸被第一道膨脹波及其反射波限制在靠近壁面的三角形區(qū)域內(nèi)，稱為下側(cè)高壓區(qū)，如圖11(c)所示。

圖11 不同時(shí)刻通道內(nèi)壓力分布 Fig.11 Pressure distribution in channel at different times

第一道反射膨脹波到達(dá)壁面后，再次發(fā)生反射，形成第二道反射膨脹波，第一道反射膨脹波的二次反射觸發(fā)了間隙上游壁面附近第二道弓形激波的形成，圖11(b) 0.032 ms、0.036 ms時(shí)的壓力分布云圖分別給出了第二道反射膨脹波和第二道弓形激波的波形。

隨著第二道反射膨脹波向?qū)ΨQ面方向傳播，第一道反射激波形成的波后高壓區(qū)域逐漸被第一、第二反射膨脹波共同限制在對(duì)稱面附近，稱為上側(cè)高壓區(qū)，如圖11(c)所示。

通過上述作用機(jī)制，泄漏發(fā)生后的第一道膨脹波(稱為主膨脹波)以近一維的形式在通道內(nèi)傳播，泄漏初始時(shí)刻主膨脹波在對(duì)稱面的反射波則在非定常泄漏過程中不斷發(fā)生反射，并周期性地觸發(fā)泄漏間隙上游靠近壁面處弓形激波的形成，該弓形激波同樣在對(duì)稱面發(fā)生反射；弓形激波及其反射激波形成的高壓區(qū)域分別被連續(xù)反射膨脹波限制在通道壁面及對(duì)稱面附近。圖11(c) 0.15 ms 時(shí)的壓力云圖給出了膨脹波的反射路徑。由于膨脹波、激波在通道壁面間的反射作用，通道內(nèi)氣體流動(dòng)方向不再與壁面平行，而是出現(xiàn)周期性折轉(zhuǎn)，如圖12所示。

圖13給出了通道中心截面上不同監(jiān)測(cè)點(diǎn)壓力隨時(shí)間的變化曲線，其中P3、P4、P5距通道末端分別為3 cm、2 cm、1 cm?？梢钥吹剑髋蛎洸ㄍㄟ^后各監(jiān)測(cè)點(diǎn)壓力均發(fā)生大幅下降，隨后在激波、局部高壓區(qū)與連續(xù)反射膨脹波的交替作用下壓力出現(xiàn)周期性震蕩，隨著非定常泄漏過程的發(fā)展，各監(jiān)測(cè)點(diǎn)壓力震蕩幅值逐漸減小。

圖14進(jìn)一步給出了各個(gè)監(jiān)測(cè)點(diǎn)壓力波動(dòng)頻率隨波動(dòng)周期的變化。由于局部高壓區(qū)范圍在傳播過程中不斷縮小，因此各個(gè)監(jiān)測(cè)點(diǎn)壓力波動(dòng)頻率隨著泄漏過程的進(jìn)行逐漸降低。

(2) 包含通道激波的情況

對(duì)于存在通道激波的泄漏過程，激波在到達(dá)通道末端時(shí)與間隙泄漏流動(dòng)存在非定常相互作用，間隙內(nèi)部流動(dòng)大致可分為3個(gè)階段：① 通道激波的繞射、反射及傳播；② 間隙出口非定常膨脹波的逆向傳播及與間隙激波的非定常相互作用；③ 準(zhǔn)穩(wěn)態(tài)流動(dòng)階段。間隙寬度在一定范圍內(nèi)變化時(shí)，前兩個(gè)階段流動(dòng)特征基本不變，僅第3階段有所區(qū)別。

圖12 0.126 ms通道內(nèi)流線分布 Fig.12 Distribution of streamlines in channel at 0.126 ms

圖13 通道中心截面不同監(jiān)測(cè)點(diǎn)壓力變化曲線 Fig.13 Pressure curves of different monitor points at centerline of channel

圖14 通道中心截面不同監(jiān)測(cè)點(diǎn)壓力波動(dòng)頻率隨波動(dòng)周期變化 Fig.14 Variation of pressure oscillation frequency with time of different monitor points at centerline of channel

選擇泄漏間隙0.2 mm、0.4 mm、0.8 mm 3個(gè)算例進(jìn)行分析，設(shè)通道激波形成瞬時(shí)為零時(shí)刻，并以間隙寬度為零時(shí)通道反射激波波后壓力為參考值對(duì)流場(chǎng)壓力數(shù)據(jù)進(jìn)行無量綱化。

圖15給出了通道激波在間隙入口繞射、反射的流動(dòng)發(fā)展過程。激波繞射使通道內(nèi)氣體流動(dòng)方向發(fā)生折轉(zhuǎn)，氣體開始泄漏到間隙內(nèi)部；當(dāng)入射激波到達(dá)通道末端發(fā)生反射時(shí)，入射激波消失，繞射波部分形成一道獨(dú)立激波在間隙內(nèi)部傳播；反射激波形成后波面向間隙一側(cè)彎曲，待其傳播到間隙左側(cè)壁面附近時(shí)，與旋渦相互作用同時(shí)波面發(fā)生分裂，其中上半支越過間隙繼續(xù)在通道內(nèi)部傳播，下半支則形成一道獨(dú)立激波在間隙內(nèi)兩側(cè)壁面間來回反射同時(shí)向間隙出口傳播。

圖15 入射激波在通道末端間隙繞射、反射靜壓云圖和流線分布 Fig.15 Static pressure and streamline distribution of incident shock wave diffraction, reflection in leakage gap of channel end

綜合以上分析，間隙內(nèi)共形成了兩道激波，正是這兩道波觸發(fā)了初始泄漏流動(dòng)，其中第一道為入射激波的繞射波分支，第二道為反射激波與旋渦相互作用引起波面分裂后的下半分支；第一道激波強(qiáng)度較弱，第二道激波則較強(qiáng)，并通過在壁面間的反射增加了第一道波的強(qiáng)度。

圖16給出了上述兩道激波在間隙內(nèi)的傳播過程?？梢钥吹?，隨著第二道激波在壁面間的不斷反射，第一道激波強(qiáng)度不斷增強(qiáng)，同時(shí)間隙左側(cè)前緣分離泡脫離壁面并向下游移動(dòng)，新的分離泡又在前緣處重新生成。

當(dāng)?shù)谝坏兰暗诙兰げㄏ嗬^到達(dá)間隙出口后，分別在出口處反射膨脹波并向上游傳播，圖17 給出了非定常膨脹波在間隙內(nèi)的逆向傳播過程。

非定常膨脹波使間隙內(nèi)氣體壓力不斷降低，間隙內(nèi)部最終會(huì)形成如圖18所示的準(zhǔn)穩(wěn)態(tài)流動(dòng)結(jié)構(gòu)，圖中δ為間隙寬度。當(dāng)間隙寬度δ變化時(shí)，在壁面前緣處形成大小不同的翼型分離泡，分離泡上表面與間隙右側(cè)壁面構(gòu)成收斂擴(kuò)張通道。分離泡前緣附近靠近壁面的氣流折轉(zhuǎn)最為劇烈，因此該處氣流在由亞聲加速到超聲過程中會(huì)在分離泡前緣產(chǎn)生一束較強(qiáng)的定常膨脹波，該膨脹波貫穿泄漏間隙并使進(jìn)入間隙的氣流在收斂擴(kuò)張通道中加速到超聲，在圖19中可以觀察到大范圍的超聲流動(dòng)區(qū)。由于在分離泡尾緣流道相對(duì)收縮，因此在尾緣處產(chǎn)生一道定常斜激波，該斜激波在間隙右側(cè)壁面反射后引起邊界層氣流分離，在圖18(a)、圖18(b)中可以看到右側(cè)壁面超聲流動(dòng)區(qū)下游形成的小分離泡。

圖16 激波在泄漏間隙內(nèi)的反射及傳播過程 Fig.16 Progress of shock wave reflection and propagation in leakage gap

圖17 非定常膨脹波在間隙內(nèi)部的逆向傳播過程 Fig.17 Backward propagation of nonsteady expansion waves in leakage gap

圖18 泄漏間隙準(zhǔn)穩(wěn)態(tài)流動(dòng)階段 Fig.18 Quasi-steady flow period in leakage gap

圖20給出了同一時(shí)刻間隙寬度δ從0～0.8 mm 變化時(shí)通道靜壓及流線分布。當(dāng)間隙寬度為零時(shí)，反射激波波后壓力分布較為均勻，與理論計(jì)算結(jié)果一致。隨著間隙寬度增加，波后出現(xiàn)了連續(xù)反射膨脹波、周期性出現(xiàn)的弓形激波及其反射波在通道內(nèi)傳播引起的壓力波動(dòng)以及由此形成的局部高壓區(qū)和低壓區(qū)，且波動(dòng)幅值隨著間隙寬度增加而增大。圖21給出了壓力波動(dòng)幅值隨間隙寬度的變化關(guān)系，間隙寬度為0.2 mm時(shí)，波動(dòng)幅值為2%，間隙寬度為0.4 mm時(shí)，波動(dòng)幅值為3.2%，間隙寬度為0.8 mm時(shí)，波動(dòng)幅值為4.4%。

圖19 0.8 mm間隙0.256 ms時(shí)刻流場(chǎng)馬赫數(shù)分布 Fig.19 Mach number distribution at 0.256 ms in 0.8 mm leakage gap

圖20 0.336 ms間隙寬度變化時(shí)通道靜壓及流線分布 Fig.20 Static pressure and streamline distribution in channel with varying gap widths at 0.336 ms

進(jìn)一步分析圖20發(fā)現(xiàn)，隨著間隙寬度增加，泄漏流量增大，激波波后平均壓力下降，從而導(dǎo)致激波強(qiáng)度減弱，激波傳播速度也隨之降低。泄漏還導(dǎo)致了激波波后附面層分離形式的改變，除了由于逆壓梯度在激波分岔腳形成的主分離泡，在膨脹波壁面反射引起的逆壓梯度作用下又形成了新的次分離泡。

觀察圖21、圖22，對(duì)于激波波后流場(chǎng)，隨著泄漏過程的進(jìn)行，同一位置不同時(shí)刻壓力均值不變，不同位置時(shí)均壓力相同。流場(chǎng)平均壓力保持不變，這是非定常泄漏過程的典型流動(dòng)特征。

圖21 不同泄漏間隙寬度P4測(cè)點(diǎn)壓力隨時(shí)間變化 Fig.21 Pressure variation with time at P4 test point with varying leakage gap widths

圖22 間隙0.8 mm算例通道中心截面不同監(jiān)測(cè)點(diǎn)壓力變化曲線 Fig.22 Pressure curves of different monitor points at centerline of channel for case with 0.8 mm gap width

對(duì)比分析不含通道激波的非定常泄漏過程發(fā)現(xiàn)，無論是否存在通道激波，通道內(nèi)均存在周期性的壓力波動(dòng)，這說明兩者具有類似的壓力波動(dòng)機(jī)制，不同的是，通道內(nèi)的主膨脹波被一道反射激波所取代。

4.2 泄漏對(duì)激波波速和激波強(qiáng)度的影響

激波波速和激波強(qiáng)度是波轉(zhuǎn)子重要設(shè)計(jì)參數(shù)，考察泄漏對(duì)它們的影響具有重要意義。

圖23給出了間隙寬度0 mm、0.8 mm 2個(gè)算例通道中心截面不同時(shí)刻的靜壓分布，圖中x為沿通道方向坐標(biāo)?？梢钥闯?，無論是否存在泄漏，激波運(yùn)行軌跡均為直線，這說明在泄漏初始時(shí)刻激波波速已經(jīng)下降到特定值，在激波傳播過程中波速并不隨著泄漏的發(fā)生而改變。因此對(duì)于給定的間隙寬度，激波馬赫數(shù)不變。

圖24給出了泄漏間隙進(jìn)出口條件相同時(shí)激波絕對(duì)馬赫數(shù)隨無量綱間隙寬度δ/H的變化，兩者近似為線性關(guān)系，δ/H=0.08時(shí)，激波絕對(duì)馬赫數(shù)相比零間隙工況衰減了7.3%。雖然激波波后壓力出現(xiàn)波動(dòng)，但壓力均值基本不變，選擇激波通過測(cè)點(diǎn)后3個(gè)完整的波動(dòng)周期對(duì)壓力進(jìn)行時(shí)均運(yùn)算，得到激波靜增壓比隨無量綱間隙寬度的變化關(guān)系，如圖25所示，兩者呈線性關(guān)系，δ/H=0.08 時(shí)，激波靜增壓比相比零間隙工況衰減了10.1%。

圖23 間隙寬度0 mm、0.8 mm通道中心截面不同時(shí)刻靜壓分布 Fig.23 Distribution of static pressure at centerline of channel for gap widths being 0 mm and 0.8 mm at different times

對(duì)于理想正激波，激波靜增壓比與激波相對(duì)馬赫數(shù)滿足式(2)。根據(jù)數(shù)值計(jì)算結(jié)果可以得到通道內(nèi)反射激波相對(duì)馬赫數(shù)，進(jìn)而由式(2)得到激波靜增壓比理論解；同樣，由圖21給出的反射激波波后壓力波動(dòng)曲線可以求出波后平均壓力，將其與波前壓力比值作為反射激波靜增壓比的數(shù)值解。圖26給出了不同間隙寬度下數(shù)值解相對(duì)理論解的偏差，可以看到，所有算例偏差均在0.01以內(nèi)，因此在分析性能參數(shù)時(shí)可以將存在泄漏時(shí)的通道反射激波視為一道理想正激波。

圖24 反射激波絕對(duì)馬赫數(shù)隨泄漏無量綱泄漏間隙寬度變化 Fig.24 Absolute Mach number of reflected shock wave with non-dimensional leakage gap width

圖25 反射激波靜增壓比隨泄漏無量綱泄漏間隙寬度變化 Fig.25 Change of static pressure ratio of reflected shock wave with non-dimensional leakage gap widths

圖26 反射激波靜增壓比計(jì)算結(jié)果相對(duì)理論偏差 Fig.26 Deviation of calculation results from theory for static pressure ratio of reflected shock wave

4.3 存在泄漏時(shí)的通道激波衰減機(jī)制

考察不同間隙寬度下通道內(nèi)不含運(yùn)動(dòng)激波時(shí)的非定常泄漏過程，將主膨脹波通過給定測(cè)點(diǎn)后的壓力做時(shí)均運(yùn)算，進(jìn)而得到時(shí)均壓力相對(duì)通道初始?jí)毫Φ乃p幅值；同樣，可以得到不同間隙寬度下反射激波波后時(shí)均壓力相對(duì)零間隙寬度時(shí)的衰減幅值。

圖27給出了上述2種情形的對(duì)比結(jié)果，其中p0為零間隙寬度時(shí)的通道初始?jí)毫?，p為存在泄漏時(shí)的通道壓力?？梢钥吹剑瑝毫λp幅值曲線幾乎重合，這說明存在泄漏時(shí)通道激波的衰減本質(zhì)上是泄漏產(chǎn)生的主膨脹波在通道內(nèi)的傳播引起的，存在泄漏時(shí)的通道激波可視為理想激波與泄漏引起的主膨脹波的線性疊加。

綜上所述，泄漏流動(dòng)通過泄漏引起的主膨脹波對(duì)激波傳播施加影響，泄漏間隙寬度一定時(shí)，主膨脹波強(qiáng)度不變，因此通道激波的衰減幅值在傳播過程中維持不變。

圖27 有/無反射激波兩種情形下不同間隙寬度通道內(nèi)壓力衰減幅值對(duì)比 Fig.27 Comparison of attenuation amplitude of channel pressure with different leakage gap widths with and without reflected shock

4.4 非定常泄漏損失預(yù)測(cè)模型

基于對(duì)通道激波衰減機(jī)制的認(rèn)識(shí)，通過對(duì)不含通道激波時(shí)非定常泄漏過程進(jìn)行?；⒎嵌ǔＰ孤p失預(yù)測(cè)模型。該模型主要基于以下幾個(gè)重要假設(shè)：① 間隙泄漏流動(dòng)等效為沿通道的均勻流動(dòng)；② 膨脹波系簡(jiǎn)化為一道膨脹波；③ 以膨脹波波頭、波尾平均速度代表簡(jiǎn)化后的膨脹波波速；④ 泄漏流動(dòng)為等熵流動(dòng)；⑤ 泄漏過程中膨脹波強(qiáng)度及波后氣體參數(shù)均保持不變。

在上述假設(shè)基礎(chǔ)上，最終得到無量綱泄漏間隙寬度G為

(4)

(5)

由于本文所有算例泄漏間隙均處于壅塞狀態(tài)，即間隙上下游壓力滿足式(6)，根據(jù)式(4)，激波衰減幅值ep僅是無量綱間隙寬度G與流量系數(shù)C0的函數(shù)。根據(jù)CFD計(jì)算結(jié)果確定C0、ep，再由泄漏損失預(yù)測(cè)模型可求得G值，從而驗(yàn)證預(yù)測(cè)模型的可靠性。表2給出了預(yù)測(cè)模型數(shù)值驗(yàn)證結(jié)果，Gmodel為預(yù)測(cè)模型G值?？梢钥闯?，兩者吻合良好。

(6)

表2泄漏流動(dòng)損失預(yù)測(cè)模型數(shù)值驗(yàn)證結(jié)果

Table2Numericalverificationresultsoflosspredictionmodelforleakageflow

GC0epGmodel0．010．80640．98730．00980．020．82320．97400．01980．040．82510．94850．03970．080．83280．89930．0792

5 結(jié) 論

1) 闡明了非定常泄漏過程轉(zhuǎn)子通道的壓力波動(dòng)機(jī)制，即連續(xù)反射膨脹波、周期性出現(xiàn)的弓形激波及其在通道中心截面的反射波引起了通道壓力的周期性波動(dòng)。

2) 闡明了間隙內(nèi)部泄漏流動(dòng)的3個(gè)發(fā)展階段及其流動(dòng)特征：① 通道激波的繞射、反射及傳播；② 間隙出口非定常膨脹波的逆向傳播及與間隙激波的非定常相互作用；③ 準(zhǔn)穩(wěn)態(tài)流動(dòng)階段。

3) 在泄漏過程中通道激波傳播速度不變、波后時(shí)均壓力不變。

4) 激波馬赫數(shù)、激波靜增壓比與無量綱間隙寬度均呈線性關(guān)系，當(dāng)間隙寬度從0增大到0.08時(shí)，激波馬赫數(shù)衰減7.3%，激波靜增壓比衰減10.1%。

5) 泄漏引起的激波衰減本質(zhì)上是主膨脹波與通道激波疊加效應(yīng)的結(jié)果。

6) 本文給出的泄漏損失預(yù)測(cè)模型可以準(zhǔn)確預(yù)測(cè)泄漏引起的激波衰減。

[1] WEBER H E. Shock wave engine design[M]. New York: John Wiley & Sons, Inc, 1995: 1-4.

[2] DEMPSEY E, MULLER N, AKBARI P, et al. Performance optimization of gas turbines utilizing four-port wave rotors: AIAA-2006-4152[R]. Reston: AIAA, 2006.

[3] AKBARI P, KHARAZI A A, MULLER N. Utilizing wave rotor technology to enhance the turbo compression in power and refrigeration cycles: IMECE-2003-44222[R]. New York: ASME, 2003.

[4] AKBARI P, NALIM R, MULLER N. A review of wave rotor technology and its applications: IMECE-2004-60082[R]. New York: ASME, 2004.

[5] KENTFIELD J A C. Nonsteady, one-dimensional, internal, compressible flows: theory and applications[M]. Oxford: Oxford University Press, 1993: 127-128.

[6] WILSON J. An experimental determination of losses in a three-port wave rotor[J]. Journal of Engineering for Gas Turbines and Power, 1998, 120(4): 833-842.

[7] PAXSON D E. Comparison between numerically modeled and experimentally measured wave-rotor loss mechanisms[J]. Journal of Propulsion and Power, 1995, 11(5): 908-914.

[8] AKBARI P, NALIM M R. Leakage assessment of pressure-exchange wave rotors[J]. Journal of Propulsion and Power, 2008, 24(4): 732-740.

[9] WILSON J, WELCH G E, PAXSON D E. Experimental results of performance tests on a four-port wave rotor: AIAA-2007-1250[R]. Reston: AIAA, 2007.

[10] ELEUTERIO F T. Riemann solvers and numerical methods for fluid dynamics[M]. Heidelberg: Springer, 2009: 115-119.

[11] GLASS I I, PATTERSON G N. A theoretical and experimental study of shock-tube flows[J]. Journal of the Aeronautical Sciences, 1955, 22(2): 73-100.

[12] ZHANG G, SETOGUCHI T, KIM H D. Numerical simulation of flow characteristics in micro shock tubes[J]. Journal of Thermal Science, 2015, 24(3): 246-253.

[13] OKAMOTO K, NAGASHIMA T. Visualization of wave rotor inner flow dynamics[J]. Journal of Propulsion and Power, 2007, 23(2): 292-300.

[14] ZHANG G, KIM H D. Numerical simulation of shock wave and contact surface propagation in micro shock tubes[J]. Journal of Mechanical Science and Technology, 2015, 29(4): 1689-1696.

[15] 閻超. 計(jì)算流體力學(xué)方法及應(yīng)用[M]. 北京: 北京航空航天大學(xué)出版社, 2006: 171.

YAN C. Methods and applications of computational fluid dynamics[M]. Beijing: Beihang University Press, 2006: 171 (in Chinese).

[16] SAITO T, TAKAYAMA K. Numerical simulations of nozzle starting process[J]. Shock Waves, 1999, 9(2): 73-79.

[17] DESHPANDE A, PURANIK B. Effect of viscosity and wall heat conduction on shock attenuation in narrow channels[J]. Shock Waves, 2015, 25(1): 1-11.

[18] SKEWS B W. The perturbed region behind a diffracting shock wave[J]. Journal of Fluid Mechanics, 1967, 29(4): 705-719.

[19] TAKAWAMA K, INOUE O. Shock wave diffraction over a 90 degree sharp corner-posters presented at 18thISSW[J]. Shock Waves, 1991, 1(4): 301-312.

[20] IGRA O, FALCOVITZ J, REICHENBACH H, et al. Experimental and numerical study of the interaction between a planar shock wave and a square cavity[J]. Journal of Fluid Mechanics, 1996, 313: 105-130.

(責(zé)任編輯： 張晗)

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

Mechanismofnonsteadyleakageflowofwaverotors

LIUChenyuan,LIUHuoxing*

SchoolofEnergyandPowerEngineering,BeihangUniversity,Beijing100083,China

Theproblemofleakageflowcanrestricttheperformanceofawaverotor.Thispaperextractsmainflowphenomenarelatedtoleakageinawaverotor,andestablishesatwo-dimensionalsinglepassageleakagemodel.Themechanismofnonsteadyleakageflowisanalyzednumerically,andthepredictionmodelforlossofunsteadyleakageflowisfurtherverified.Theresultsshowthatthereexistspressureoscillationofvaryingextentinrotorpassage,andtheoscillationamplitudeisrelatedtothegapwidth.Thecontinuousreflectedexpansionwave,periodicemergingbowshockanditsreflectedshockarereasonsforpressureoscillation.Thereexistthreemainflowstagesintheleakagegapduringleakageflowprocessandduringwhichthepropagationspeedofpassageshockandpost-shocktime-averagingpressurestayunchanged.Withinacertainrangeofgapwidth,shockMachnumberandshockstaticpressureratioarebothlinearlydependentonthenondimensionalleakagegapwidth.Whenthegapwidthisincreasedfrom0to0.08,theshockMachnumberdecreasesby7.3%,andtheshockstaticpressureratiodecreasesby10.1%.Thepropagationofchannelshockwaveisinfluencedbyleakageflowprocessviamainexpansionwave,andshockwaveattenuationisessentiallytheresultofsuperpositionofidealshockwaveandexpansionwave.TheresultsobtainedwithleakagelosspredictionmodelisfoundtobeconsistentwithCFDresults.

gasturbine;waverotor;singlepassagemodel;nonsteadyleakageflow;shockwave

2016-07-12；Revised2016-08-24；Accepted2016-10-08；Publishedonline2016-10-101621

.E-mailliuhuoxing@buaa.edu.cn

2016-07-12；退修日期2016-08-24；錄用日期2016-10-08； < class="emphasis_bold">網(wǎng)絡(luò)出版時(shí)間

時(shí)間：2016-10-101621

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

.E-mailliuhuoxing@buaa.edu.cn

劉琛源， 劉火星． 波轉(zhuǎn)子非定常泄漏流動(dòng)機(jī)理J． 航空學(xué)報(bào),2017,38(5):120606.LIUCY,LIUHX.Mechanismofnon-steadyleakageflowofwaverotorsJ.ActaAeronauticaetAstronauticaSinica,2017,38(5):120606.

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

10.7527/S1000-6893.2016.0269

V231.3

A

1000-6893(2017)05-120606-14