王洪平, 高 琪,*, 魏潤杰, 王晉軍

(1. 北京航空航天大學(xué)流體力學(xué)研究所 流體力學(xué)教育部重點實驗室, 北京 100191; 2. 北京立方天地科技有限責任公司, 北京 100083)

基于層析PIV的湍流邊界層展向渦研究

王洪平1, 高 琪1,*, 魏潤杰2, 王晉軍1

(1. 北京航空航天大學(xué)流體力學(xué)研究所 流體力學(xué)教育部重點實驗室, 北京 100191; 2. 北京立方天地科技有限責任公司, 北京 100083)

層析粒子圖像測速;湍流邊界層;旋渦識別;展向渦

0 引 言

湍流中的擬序結(jié)構(gòu)被認為是湍流產(chǎn)生和發(fā)展的關(guān)鍵,因此針對擬序結(jié)構(gòu)的研究一直是湍流研究的一個重要方向。近期的湍流研究主要關(guān)注擬序結(jié)構(gòu)和湍流統(tǒng)計量之間的關(guān)系。Smits[1]在關(guān)于高雷諾數(shù)壁

湍流的綜述性文章中將擬序結(jié)構(gòu)主要分為4類,包括近壁區(qū)條帶,馬蹄渦,大尺度結(jié)構(gòu)(Large-Scale Motions,LSMs)和超大尺度結(jié)構(gòu)(Very-Large-Scale Motions,VLSMs)。這里馬蹄渦是對稱以及不對稱馬蹄渦、Ω形渦和拱形渦的統(tǒng)稱,這些渦一般由1個或2個流向渦腿和1個展向渦頭組成,展向渦頭的旋轉(zhuǎn)方向和湍流邊界層的平均剪切方向一致[2]。Adrian[3]用平面粒子圖像測速(Particle Image Velocimetry, PIV)實驗證明了在湍流邊界層外區(qū)存在馬蹄渦,同時這些馬蹄渦的渦腿附著在壁面上并沿流向排列形成了流向渦包結(jié)構(gòu)。目前,這種附著馬蹄渦模型已被廣泛接受[4]。但是,隨著數(shù)值計算以及實驗手段的不斷改進和完善,有學(xué)者提出了不同的意見。Yang[5]通過馬蹄渦發(fā)生器產(chǎn)生馬蹄渦并觀察馬蹄渦向下游發(fā)展的情況,發(fā)現(xiàn)由渦產(chǎn)生的噴射(Ejection)現(xiàn)象能夠?qū)е露务R蹄渦的產(chǎn)生,并且能夠誘導(dǎo)產(chǎn)生反向展向渦。Pirozzoli[6]通過條件平均進一步闡明了湍流中渦管和剪切層的相互作用,發(fā)現(xiàn)渦管基本分布在剪切層上并推斷前者可能是Kelvin-Helmholtz不穩(wěn)定性的產(chǎn)物。Lozano-Durán[7]在追蹤擬序結(jié)構(gòu)發(fā)展變化的過程中發(fā)現(xiàn)并不是所有的壁面附著結(jié)構(gòu)都產(chǎn)生于壁面。Eitel-Amor[8]重做了Zhou[9]的數(shù)值實驗,指出以前的研究[9-10]可能過高地估計了馬蹄渦存在的生命周期以及再生的重要性和魯棒性,在湍流耗散環(huán)境中馬蹄渦會很快衰減,不足以維持湍流大尺度運動。

通過上面的討論可知,馬蹄渦雖然能夠成功地描述湍流邊界層中的一些現(xiàn)象，但仍有學(xué)者對其在充分發(fā)展湍流中的重要性提出了疑問。在本文中,定義正展向渦和湍流邊界層中馬蹄渦頭的方向一致,用“p”表示(Prograde Spanwise Vortices),負展向渦用“r”表示(Retrograde Spanwise Vortices),與湍流邊界層平均剪切方向相反。眾所周知,在湍流邊界層中存在和馬蹄渦渦頭相反的負展向渦,而且總是和正展向渦成對出現(xiàn)[2,6,11]。隨著雷諾數(shù)的升高,負展向渦的作用將變得更加明顯[2]。目前,對負展向渦的產(chǎn)生機理主要有3種解釋:正展向渦誘導(dǎo)局部剪切失穩(wěn)[5],馬蹄渦的相互合并[12]以及Ω形渦的一部分[11]。然而,由于旋渦結(jié)構(gòu)復(fù)雜的三維特性,目前這些解釋仍缺少三維流場的進一步實驗驗證。

想要完美地解釋負展向渦的產(chǎn)生機理以及和正展向渦的關(guān)系,就必須得到時間解析的瞬時三維三分量(3D3C)速度場。隨著科學(xué)技術(shù)的發(fā)展和實驗手段的不斷完善,層析PIV技術(shù)[13-14]的出現(xiàn)使這種想法成為可能。從本質(zhì)上講,層析PIV和平面PIV的原理相同,它通過多相機多視角成像,然后利用倍增代數(shù)重構(gòu)技術(shù)(Multiplicative Algebraic Reconstruction Technique,MART)得到流場示蹤粒子的三維空間分布,進而互相關(guān)分析得到3D3C速度場。經(jīng)過近10年的發(fā)展,層析PIV已經(jīng)被用于測量高低雷諾數(shù)湍流邊界層[15-16]、零質(zhì)量射流[17]、激波邊界層相互作用[18]、三角翼前緣渦[19]以及渦環(huán)撞壁實驗[20]等,得到了復(fù)雜流場的三維結(jié)構(gòu)。

本文運用層析PIV技術(shù)測量了Reτ≈1700的湍流邊界層,并且從統(tǒng)計角度對正負展向渦進行了分析,希望對湍流邊界層中的擬序結(jié)構(gòu)有更深入的理解。本文的組織結(jié)構(gòu)如下:第1節(jié)和第2節(jié)分別介紹了層析PIV的具體實驗方法以及旋渦辨識的方法,這2部分是數(shù)據(jù)分析的基礎(chǔ);第3節(jié)統(tǒng)計了正負展向渦沿法向的變化規(guī)律;第4節(jié)進一步討論了正負展向渦的空間位置對流場速度的影響;第5節(jié)討論了幾個典型的三維流場結(jié)構(gòu);最后給出總結(jié)和結(jié)論。

1 實驗方法

本實驗在北京航空航天大學(xué)大水洞完成,整個測量工況如圖 1所示,其中x,y,z分別代表流向、法向和展向。為了生成滿足要求的湍流邊界層,在大水洞中豎直安裝了一塊大小為7m×1m的有機玻璃板并且前緣為橢圓倒角,直徑為3mm的拌線布置在距平板前緣0.1m的位置。水洞底壁的寬度為1m,玻璃板被安裝在距側(cè)壁0.6m的位置,實驗時水位高度也為0.6m。測量體處在展向中心,距玻璃板前緣6.7m,下邊界距平板壁面的高度為5.5mm。流場中投放平均粒徑為20μm的空心玻璃微珠作為示蹤粒子,其密度為1.05×10-3g·mm-3。4臺分辨率為2456pixel×2048pixel的12位CCD相機呈“×”形布置在水洞一側(cè),相機之間的夾角約為60°,各臺相機配置尼康f=45mm定焦體視鏡頭。為了減少由于視角引起的光學(xué)畸變,鏡頭的光圈f#被調(diào)到8和11之間,并且相機、鏡頭和物面滿足Scheimpflug條件。光源采用500mJ雙脈沖激光(鐳寶光電Vlite-500 PIV專用雙脈沖激光器),波長532nm,經(jīng)過擴束光路被擴展成厚度為16mm的體光源。相機和激光之間通過同步器(北京立方天地SM-MicroPulse725同步器)進行控制,實驗時激光頻率為4Hz,2次曝光的時間間隔Δt=1.2ms。每臺相機采樣了300幅圖像,采樣時長約為38s,最終得到了150個速度場。

圖1 層析PIV實驗工況示意圖

層析PIV實驗前需要標定物理空間和像素空間的關(guān)系。本文采用多項式標定(x,z方向采用三次多項式,y向采用一次多項式),標定板在測量體厚度范圍內(nèi)以2mm等間距移動并成像。初始得到的映射函數(shù)需要進行自修正[21],將映射函數(shù)的誤差控制在0.1pixel以內(nèi)，具體處理過程可以參考文獻[22]。在層析PIV中,圖像前處理也是一個非常重要的環(huán)節(jié),可以很大程度上提升重構(gòu)的精度和效率。本文對原始圖像進行了時序最小背景剔除和空間滑動背景剔除,最終的圖像進行了3×3的高斯平滑,高斯核σ=0.5。經(jīng)過圖像前處理后,粒子濃度約為0.06ppp(particle per pixel)。結(jié)合映射函數(shù),空間粒子灰度場通過MART算法重構(gòu)得到。本文進行10次MART迭代,重構(gòu)體的大小為80mm×16mm×45mm,圖像放大率為0.055mm/pixel,因此重構(gòu)體素(voxel)的個數(shù)為1454voxel×290voxel×818voxel。

重構(gòu)后的粒子灰度場直接通過互相關(guān)分析確定粒子的位移。本文采用多層迭代體變形算法[23],最終窗口大小為48voxel×48voxel×48voxel,重疊率為75%。為了提高計算效率,第1層互相關(guān)采用快速傅里葉變換算法(Fast Fourier Transform, FFT),其余互相關(guān)在得到速度參考矢量后進行窗口預(yù)偏置,然后采用小位移直接互相關(guān)算法。在體變形時只采用線性插值,這已經(jīng)能滿足計算精度的要求。另外,每一層互相關(guān)得到的速度場需要進行歸一化壞點剔除[24]和高斯濾波,才能成為下一層互相關(guān)迭代的參考速度矢量。最終得到的速度場在剔除壞點的基礎(chǔ)上進行了高斯平滑,平滑窗的大小為3×3×3,平滑強度σ=1。

流場的品質(zhì)和參數(shù)通過激光多普勒測速(LaserDopplerVelocimetry,LDV)進行標定,表 1給出了流場的一些基本信息。實驗時水溫為18℃,對應(yīng)的運動粘性系數(shù)ν為1.055mm2/s。當自由來流速度U∞約為0.413m/s時,流場的湍流度為0.97%。通過Musker方法[25]計算得到壁面摩阻速度uτ為15.57mm/s,邊界層厚度δ為120mm,對應(yīng)的基于壁面摩阻速度和動量厚度的雷諾數(shù)分別為1769和5176。根據(jù)公式用壁面粘性尺度(WallUnits,WU)對物理量進行無量綱化。

(1)

因此,粘性尺度的測量體大小為1176WU×235WU×661WU。速度場的法向測量范圍為100

表1 湍流邊界層的基本參數(shù)

2 旋渦及渦心辨識方法

2.1 三維流場的旋渦辨識

在湍流邊界層中,存在渦量的地方可能存在旋渦,也可能是剪切層,因此,必須從湍流背景流場中提取真正的旋渦。根據(jù)Robinson[26]關(guān)于旋渦的定義,旋渦通常會以一定的速度對流,需要考慮伽利略不變性(Galilean Invariant),但在湍流邊界層中旋渦的對流速度往往與所處位置的流動結(jié)構(gòu)有關(guān)。目前的旋渦識別方法大多通過分析局部速度梯度張量來實現(xiàn),例如Q準則、Δ準則、λ2準則以及λci準則,這些判據(jù)之間存在一定的聯(lián)系,某些情況下是相互等價的[27]。本文采用基于當?shù)孛}動速度u的λci準則,即旋渦強度判據(jù),該判據(jù)用局部速度梯度的復(fù)特征值的虛部來表示旋轉(zhuǎn)的強度。由于分析的是脈動速度,湍流邊界層當?shù)貢r均速度U(僅考慮流向速度)已經(jīng)從流場中剔除?？紤]到PIV數(shù)據(jù)存在噪音以及空間分辨率的問題,速度梯度的計算往往存在誤差。同時,旋渦強度在湍流邊界層中并不是均勻分布,而是與法向高度y+有關(guān)。因此,閾值的設(shè)置不僅要考慮誤差,還要考慮y+的影響。根據(jù)其他學(xué)者的研究成果[2,11],本文閾值的設(shè)定如公式(2)所示。

(2)

(3)

2.2 三維流場中展向渦位置

本文對展向渦的定位采用了類似文獻[28]的方法,即找出展向渦所對應(yīng)的連通區(qū)域,以該連通區(qū)域的質(zhì)心作為三維流場中展向渦的位置。具體計算過程為:

3 展向渦沿法向分布規(guī)律

3.1 展向渦強度

(a) 代表平均值

(b) 代表脈動值

3.2 展向渦對流速度

(a)

(b)

圖3 正負展向渦的流向?qū)α魉俣?a)和法向?qū)α魉俣?b)沿高度的變化。(a)中的黑色虛線代表速度剖面,(b)中的黑色虛線代表V+=0。

Fig.3 (a) The mean streamwise velocity of prograde and retrograde spanwise vortices as a function ofy+,theblackdottedlinerepresentsthemeanvelocityprofile. (b)Themeanwall-normalvelocityofprogradeandretrogradespanwisevorticesasafunctionofy+,theblackdottedlinerepresentsV+=0.

(a) 正展向渦 (b) 負展向渦 (c) 全流場

圖4 流向速度和法向速度的聯(lián)合概率密度分布

Fig.4 The joint probability density function of streamwise velocity and wall-normal velocity. (a) prograde spanwise vortices, (b) retrograde spanwise vortices, (c) the whole flow field.

4 正負展向渦的空間關(guān)系

通過第3節(jié)的分析可知:正負展向渦在法向的區(qū)別主要集中在強度、空間流場結(jié)構(gòu)以及對流速度上,本節(jié)主要討論正負展向渦在x-y平面內(nèi)的空間關(guān)系。以正展向渦的渦心為中心,在半徑為200WU的區(qū)域內(nèi)統(tǒng)計周圍正負展向渦出現(xiàn)的概率,在展向z允許有±2Δz+(±19.4WU)的偏差。在統(tǒng)計過程中,同一個展向渦如果同時被不同的中心渦選中,則選擇其中距離最近的一個。圖5為統(tǒng)計得到的結(jié)果,圖5(a)表示中心為正展向渦時周圍正展向渦的分布情況，而圖5(b)表示中心是正展向渦,負展向渦在其周圍的分布情況，橫縱坐標分別表示離中心渦的流向和法向距離,云圖代表周圍渦出現(xiàn)的概率,顏色越紅說明在該位置出現(xiàn)的概率越高。在圖5中,正展向渦基本上沿流向水平排列,間距約為60WU，由于樣本數(shù)量有限導(dǎo)致統(tǒng)計的分辨率較低。對于圖5(b)中的負展向渦,基本上分布在正展向渦的流向和法向。Natrajan[11]在文章中展示了法向分布的這種情況,但并沒有展示流向分布這種情況,在實際流場中,后者往往也是存在的，其出現(xiàn)的概率幾乎和法向出現(xiàn)的概率相等。另外,仔細觀察負展向渦法向出現(xiàn)的位置可以發(fā)現(xiàn)：在正展向渦底部時其偏向上流,而在正展向渦上部時偏向下游,這與Natrajan[11]的結(jié)果一致。

進一步分析正負展向渦的空間位置對流場速度的影響。如圖6(a)所示,以正展向渦為中心,負展向渦的位置可以用他們之間的連線與流向的夾角θ表示(半徑變化較小,一般在100個粘性尺度以內(nèi),如圖5(b)所示)。θ的取值范圍為0°到360°,180°和0°表示負展向渦在正展向渦的上下游,而90°和270°表示在上層和下層。連線中心的脈動速度可以由插值得出,用該速度可以比較容易和準確地判斷出流場速度與展向渦空間位置的關(guān)系。圖6(b)給出了連線中心的脈動速度(u+,v+)隨角度θ的變化規(guī)律,圓點代表流向速度u+,三角代表法向速度v+。隨著角度的變化,流向速度呈現(xiàn)出正弦曲線的變化規(guī)律,而法向速度表現(xiàn)出與之相對的余弦變化趨勢。當負展向渦分布在正展向渦的流向,即θ為0°或180°時,流向速度的平均值幾乎為0,而法向速度的平均值卻最大,約為1.2;當負展向渦分布在正展向渦的法向時,即θ為90°或270°時,流向速度表現(xiàn)出高低速區(qū)域,而法向速度卻維持在0左右,如圖6(a)中的紅藍箭頭所示。該規(guī)律完全背離了圖4所描述的速度之間的依賴關(guān)系,高低速區(qū)的產(chǎn)生無法簡單地用掃掠和噴射來解釋。需要注意的是,這種和渦結(jié)構(gòu)有很強關(guān)系的高低速區(qū)域的尺度都比較小,與圖5中展向渦之間的距離相當。一般來說,噴射產(chǎn)生于近壁區(qū)并逐漸向上運動,相反,掃掠產(chǎn)生于較高位置朝壁面運動[7]。這種大范圍相對運動容易產(chǎn)生大尺度結(jié)構(gòu),而中小尺度的流動結(jié)構(gòu)更多是由于不穩(wěn)定性和旋渦運動造成的。

(a)

(b)

Fig.5 The distribution of prograde (a) and retrograde (b) spanwise vortices around a prograde vortex

(a)

圖6 (a)高低速區(qū)域與正負展向渦位置關(guān)系示意圖,(b)速度與正負展向渦空間位置的關(guān)系

Fig.6 (a) The schematic diagram about the relationship between the low (high) momentum region and the spatial arrangement of spanwise vortices, (b) The corresponding variations inu+(v+)andθ

5 三維流場

6 結(jié) 論

本文通過層析PIV測量了Reτ≈1700的湍流邊界層,并統(tǒng)計分析了100

(1) 對于旋渦強度,正負展向渦表現(xiàn)出不同的性質(zhì)。隨著法向高度的增加正展向渦逐漸變?nèi)醵撜瓜驕u基本保持不變。

(2) 流向速度和法向速度呈現(xiàn)出很強的相關(guān)性,導(dǎo)致流場中的大尺度高低速區(qū)域成為雷諾應(yīng)力-uv的主要來源。

(3) 正展向渦在流向呈現(xiàn)出等間距分布,而負展向渦大多分布在正展向渦的流向和法向，但并沒有統(tǒng)計這種分布在流向持續(xù)的長度。

(4) 在小尺度范圍內(nèi)(約為100WU),高低速區(qū)域并不受噴射和掃掠機制的控制,而是更多和渦結(jié)構(gòu)的空間位置以及局部不穩(wěn)定性有關(guān)。

[1] Smits A J, Mckeon B J, Marusic I. High-Reynolds number wall turbulence[J]. Annual Review of Fluid Mechanics, 2011, 43(1): 353-375.

[2] Wu Y, Christensen K T. Population trends of spanwise vortices in wall turbulence[J]. J Fluid Mech, 2006, 568: 55-76.

[3] Adrian R J, Meinhart C D, Tomkins C D. Vortex organization in the outer region of the turbulent boundary layer[J]. J Fluid Mech, 2000, 422: 1-54.

[4] Adrian R J. Hairpin vortex organization in wall turbulence[J]. Physics of Fluids, 2007, 19(4).

[5] Yang W, Meng H, Sheng J. Dynamics of hairpin vortices generated by a mixing tab in a channel flow[J]. Exp Fluids, 2001, 30(6): 705-722.

[6]Pirozzoli S. Flow organization near shear layers in turbulent wall-bounded flows[J]. Journal of Turbulence, 2011: N41.

[7] Lozano-Durán A, Jiménez J. Time-resolved evolution of coherent structures in turbulent channels: characterization of eddies and cascades[J]. J Fluid Mech, 2014, 759: 432-471.

[8] Eitel-Amor G, ?rlü R, Schlatter P, et al. Hairpin vortices in turbulent boundary layers[J]. Physics of Fluids, 2015, 27(2): 025108.

[9] Zhou J, Adrian R J, Balachandar S, et al. Mechanisms for generating coherent packets of hairpin vortices in channel flow[J]. J Fluid Mech, 1999, 387: 353-396.

[10] Kim K, Sung H J, Adrian R J. Effects of background noise on generating coherent packets of hairpin vortices[J]. Physics of Fluids, 2008, 20(10): 105107.

[11] Natrajan V K, Wu Y, Christensen K T. Spatial signatures of retrograde spanwise vortices in wall turbulence[J]. J Fluid Mech, 2007, 574: 155-167.

[12] Adrian R J, Balachandar S, Lin Z C. Spanwise growth of vortex structure in wall turbulence[J]. KSME International Journal, 2001, 15(12): 1741-1749.

[13] Gao Q, Wang H P, Shen G X. Review on development of volumetric particle image velocimetry[J]. Chinese Science Bulletin, 2013, 58(36): 4541-4556.

[14] Elsinga G E, Scarano F, Wieneke B, et al. Tomographic particle image velocimetry[J]. Exp Fluids, 2006, 41(6): 933-947.

[15] Elsinga G E, Adrian R J, Van Oudheusden B W, et al. Three-dimensional vortex organization in a high-Reynolds-number supersonic turbulent boundary layer[J]. J Fluid Mech, 2010, 644: 35-60.

[16] Gao Q. Evolution of eddies and packets in turbulent boundary layers[D]. Minnesota: The University of Minnesota, 2011.

[17] 高琪, 王洪平. 層析PIV技術(shù)及其合成射流測量[J]. 中國科學(xué):技術(shù)科學(xué), 2013, (7): 828-835.

Gao Q, Wang H P. The tomographic PIV and it's application on synthetic jet measurement[J]. Science China: Technological Sciences, 2013, (7): 828-835.

[18] Humble R A, Elsinga G E, Scarano F, et al. Three-dimensional instantaneous structure of a shock wave/turbulent boundary layer interaction[J]. J Fluid Mech, 2009, 622: 33-62.

[19] 王成躍, 高琪, 魏潤杰, 等. 三角翼前緣渦三維流動顯示和層析PIV測量[C]. 第八屆全國流體力學(xué)學(xué)術(shù)會議, 2014.

Wang C Y, Gao Q, Wei R J, et al. 3D flow visualization and tomographic particle image velocimetry for vortex breakdown of a delta wing[C]. The 8th National Congress of Fluid Mechanics, 2014.

[20] Violato D, Ianiro A, Cardone G, et al. Investigation on circular and chevron impinging jets by IR thermography and time-resolved tomographic PIV[C]. ASME-JSME-KSME 2011 Joint Fluids Engineering Conference, 2011.

[21] Wieneke B. Volume self-calibration for 3D particle image velocimetry[J]. Exp Fluids, 2008, 45(4): 549-556.

[22] Gao Q, Wang H P, Wang J J. A single camera volumetric particle image velocimetry and its application[J]. Science China: Technological Sciences, 2012, 55(9): 2501-2510.

[23] Scarano F. Iterative image deformation methods in PIV[J]. Meas Sci Technol, 2002, 13(1): R1-R19.

[24] Westerweel J, Scarano F. Universal outlier detection for PIV data[J]. Exp Fluids, 2005, 39(6): 1096-1100.

[25] Kendall A, Koochesfahani M. A method for estimating wall friction in turbulent wall-bounded flows[J]. Exp Fluids, 2008, 44(5): 773-780.

[26] Robinson S K. Coherent motions in the turbulent boundary layer[J]. Annual Review of Fluid Mechanics, 1991, 23(1): 601-639.

[27] Chakraborty P, Balachandar S, Adrian R J. On the relationships between local vortex identification schemes[J]. J Fluid Mech, 2005, 535: 189-214.

[28] Del Alamo J C, Jiménez J, Zandonade P, et al. Self-similar vortex clusters in the turbulent logarithmic region[J]. J Fluid Mech, 2006, 561: 329-358.

(編輯：張巧蕓)

Study of spanwise vortices in turbulent boundary layer flow based on tomographic PIV

Wang Hongping1, Gao Qi1,*, Wei Runjie2, Wang Jinjun1

(1. Fluid Mechanics Key Laboratory of Ministry of Education, Beijing University of Aeronautics and Astronautics, Beijing 100191, China; 2. MicroVec., Inc, Beijing 100083, China)

tomographic PIV;turbulent boundary layer;vortex identification;spanwise vortices

1672-9897(2016)02-0059-08

10.11729/syltlx20150086

2015-06-08;

2015-07-30

國家自然科學(xué)基金(11472030,11327202,11490552)

WangHP,GaoQ,WeiRJ,etal.StudyofspanwisevorticesinturbulentboundarylayerflowbasedontomographicPIV.JournalofExperimentsinFluidMechanics, 2016, 30(2): 59-66. 王洪平, 高 琪, 魏潤杰, 等. 基于層析PIV的湍流邊界層展向渦研究. 實驗流體力學(xué), 2016, 30(2): 59-66.

V211.7;O

A

王洪平(1987-),男，四川廣元人,博士。研究方向:實驗流體力學(xué)。通信地址: 北京航空航天大學(xué)航空科學(xué)與工程學(xué)院流體力學(xué)研究所(100191)。E-mail: h.p.wang@foxmail.com。

*通信作者E-mail:qigao@buaa.edu.cn