基于網(wǎng)格變形的葉片多目標氣動優(yōu)化設計*

劉浩　張雷　李霄琳

(吉林大學 機械科學與工程學院, 吉林 長春 130025)

基于網(wǎng)格變形的葉片多目標氣動優(yōu)化設計*

劉浩張雷?李霄琳

(吉林大學 機械科學與工程學院, 吉林 長春 130025)

摘要：壓氣機葉片的傳統(tǒng)優(yōu)化設計存在設計變量多和優(yōu)化周期長等不足.為此,文中利用非均勻有理B樣條基函數(shù)建立了轉子流場網(wǎng)格自由變形參數(shù)化方法,并結合改進的拉丁超立方試驗設計、Kriging響應面模型及NSGA-Ⅱ多目標遺傳算法構建了轉子葉片的氣動優(yōu)化設計體系.計算結果表明:在98%的堵塞質量流量工況下,優(yōu)化后的葉片總壓比提高了0.33%,等熵效率提高了0.83%;優(yōu)化后壓氣機轉子形狀為前傾型葉片,提升了轉子性能,降低了激波損失;與傳統(tǒng)優(yōu)化設計方法相比,文中優(yōu)化體系的設計變量明顯減少,縮短了優(yōu)化設計周期.

關鍵詞：壓氣機;優(yōu)化設計;網(wǎng)格變形;響應面;多目標遺傳算法

隨著航空發(fā)動機性能要求的不斷提升,在現(xiàn)有發(fā)動機型號的基礎上進行優(yōu)化改進和縮放是業(yè)內普遍采用的方法,優(yōu)化壓氣機葉片結構是提升發(fā)動機性能的有效途徑.計算機技術及計算流體動力學(CFD)技術的發(fā)展,使發(fā)動機葉片的數(shù)值優(yōu)化成為設計過程中必不可少的環(huán)節(jié).國內外學者對此進行了大量研究,Lian等[1]運用遺傳算法(GA)、梯度基算法和響應面方法對Rotor67的壓比、重量進行多目標優(yōu)化.Okui等[2]運用人工神經網(wǎng)格遺傳算法對壓氣機葉片進行優(yōu)化以提高轉子效率.李彬等[3]建立了基于自適應差分遺傳算法的長葉片透平多學科多目標優(yōu)化設計系統(tǒng).

上述葉片優(yōu)化設計方法是將葉片截面輪廓和積疊線相關的幾何變量作參數(shù)化,優(yōu)化迭代的每一步都會重復幾何形狀生成和網(wǎng)格劃分,降低了優(yōu)化的效率和網(wǎng)格的質量.特別是在復雜幾何特征位置處和邊界層的網(wǎng)格,容易出現(xiàn)網(wǎng)格劃分失敗.此外,傳統(tǒng)優(yōu)化方法對幾何模型的參數(shù)化要求高,幾何模型易更新失敗.

為解決上述問題,將網(wǎng)格變形[4]引入壓氣機葉片結構的優(yōu)化設計中,通過更新網(wǎng)格模型的網(wǎng)格節(jié)點位置產生新的形狀,不需要對幾何模型進行參數(shù)化及重復劃分網(wǎng)格,提高了計算效率.此外,變形后的網(wǎng)格與原網(wǎng)格具有一致的網(wǎng)格拓撲結構,保證了網(wǎng)格質量,消除了網(wǎng)格差異對計算流體動力學計算結果分析產生的干擾.

目前,運用網(wǎng)格變形技術進行結構優(yōu)化設計已成為國內外研究的熱點.王元元等[5]將自由變形(FFD)技術應用于運輸機的氣動優(yōu)化中,數(shù)值模擬和風洞試驗結果表明,優(yōu)化結果達到了改善后體流動、減小后體阻力的目的.Yusuke等[6]采用FFD技術對多體船進行多目標優(yōu)化設計,并對優(yōu)化后的結果進行了模型試驗驗證.Biancolini等[7]運用徑向基函數(shù)網(wǎng)格變形技術對船的帆型進行優(yōu)化設計,體現(xiàn)了網(wǎng)格變形技術的優(yōu)越性.

文中利用網(wǎng)格變形技術,結合改進拉丁超立方體試驗設計、Kriging響應面模型及多目標遺傳算法對NASA Rotor37軸流壓氣機葉片葉型進行多目標氣動優(yōu)化,以網(wǎng)格變形參數(shù)為設計變量,等熵效率、總壓比為目標函數(shù)進行優(yōu)化設計.

1數(shù)值模型

文中基于Fluent軟件平臺,以Rotor37為研究對象,對雷諾平均N-S控制方程進行三維定常計算.采用Spalart-Allmaras湍流模型,對流項采用二階精度的迎風格式進行離散.為保證網(wǎng)格在葉片前、后緣附近有較好的正交性,采用了分區(qū)、貼體的網(wǎng)格拓撲結構.流道葉身網(wǎng)格為77×57×209的結構型網(wǎng)格,葉頂間隙采用25×33×209的蝶形網(wǎng)格,網(wǎng)格總數(shù)約為155萬,并且保證了計算結果的網(wǎng)格無關性.根據(jù)試驗條件設定邊界條件,入口邊界給定平均分布的總壓、總溫;出口邊界給定輪轂半徑處的靜壓,其余出口壓力按徑向平衡方程計算得到;壁面設定為絕熱、固壁、無滑移邊界條件.

2網(wǎng)格變形及參數(shù)化

文中在圓柱坐標系下構造了FFD控制體,控制體的形狀與轉子流場外形相似,控制體是連續(xù)體,控制頂點可以任意分布.

Rotor37流場網(wǎng)格變形的基本步驟如下:

(1)合并3個參數(shù)控制體(流場入口至葉片前緣區(qū)域、葉片前緣至后緣區(qū)域、葉片后緣至出口區(qū)域)之間的交界面成為一個控制體.控制體的參數(shù)空間為Q(u,v,w)、局部坐標系O-UVW為圓柱坐標系,X0(O)是局部坐標系的原點,U、V、W是軸矢量,則任意一點X在局部坐標系中的坐標表示為

X(u,v,w)=X0+uU+vV+wW

(1)

參數(shù)空間坐標為

(2)

(3)

(4)

且滿足0

(2)添加控制頂點Pi,j,k,并將待變形網(wǎng)格映射到參數(shù)空間,計算每個網(wǎng)格節(jié)點X在參數(shù)空間中的參數(shù)坐標,則任意網(wǎng)格節(jié)點在參數(shù)空間中的坐標為

(5)

(0≤u,v,w≤1)

(6)

式中, Wi,j,k為Pi,j,k的權因子,Bi,p(u)、Bj,q(v)、Bk,o(w)分別為定義在節(jié)點向量上的p階、q階、o階非均勻有理B樣條基函數(shù),u={u0,u1,…,up,…,ul,…, ul+p},v={v0,v1,…,vq,…,vj,…,vj+q},w={w0,w1,…,wo,…,wk,…,wk+o}.

(4)保持每個網(wǎng)格節(jié)點在參數(shù)空間中相應的參數(shù)坐標不變,由于控制體的形狀發(fā)生改變,故每個網(wǎng)格節(jié)點的空間位置也會發(fā)生改變,從而實現(xiàn)網(wǎng)格的變形.設網(wǎng)格節(jié)點X移動到X′,則X′的坐標可表示為

(7)

圖1給出了Rotor37流場網(wǎng)格和參數(shù)化控制體構造過程.圖中的4個控制點接近葉片葉頂部位,改變其位置葉片隨之發(fā)生明顯的變形.相關研究表明,Rotor37的傾斜和掠型葉片有助于提高轉子的氣動性能[8-10].葉片隨4個控制點沿不同方向移動而產生的傾斜和掠型效果如圖2所示.表1列出了控制點經過許用最大位移后,得到變形后網(wǎng)格的質量數(shù)據(jù),對比發(fā)現(xiàn)變形前后的網(wǎng)格質量參數(shù)變化很小,變形后的網(wǎng)格與原網(wǎng)格保持同一網(wǎng)格質量精度.

圖1　網(wǎng)格變形參數(shù)控制體Fig.1　Parametric deformable control volume of mesh

圖2　變形后的傾斜和掠型葉片F(xiàn)ig.2　Leaned and swept blades after deformation

表1　原網(wǎng)格和變形后網(wǎng)格的質量對比Table 1　Comparison of quality between the original mesh and the deformed mesh

3優(yōu)化方法與設計

優(yōu)化步驟如下:①對所要優(yōu)化的問題進行數(shù)值描述,劃分網(wǎng)格;②參數(shù)化模型,確立設計變量;③通過試驗設計生成樣本點;④對樣本點進行數(shù)值模擬,獲得響應值數(shù)據(jù)庫;⑤采用Kriging響應面方法構造近似模型;⑥確立目標函數(shù)和約束條件,并采用多目標遺傳算法進行全局尋優(yōu);⑦對最優(yōu)設計點進行數(shù)值流場驗證,判斷數(shù)值模擬結果與近似模型的響應值偏差是否在誤差范圍內;⑧將不滿足條件的數(shù)值模擬解加入數(shù)據(jù)庫,返回步驟⑤直至誤差在允許范圍內,得到優(yōu)化結果,最終輸出網(wǎng)格文件和幾何文件.

3.1　響應面法

選用非線性程度較高和對局部響應突變問題具有較好擬合效果的Kriging模型[11]作為響應面近似模型.該模型已在航空航天、船舶、汽車等工程的優(yōu)化設計領域得到廣泛的應用.

根據(jù)Kriging模型原理,可將未知函數(shù)關系表示為

y(x)=f(x)+Z(x)

(8)

式中,未知函數(shù)f(x)表示全局近似,Z(x)表示服從正態(tài)分布N(0,σ2)的隨機誤差模型,用于局部近似.

隨機函數(shù)Z(x)的協(xié)方差矩陣為

cov[Z(x(i)),Z(x(j))]=σ2R[R(x(i),x(j))]

(9)

式中:i,j∈{1,2,…,ns};R(x(i),x(j))為樣本空間中任意兩個樣本點x(i)和x(j)的空間相關函數(shù);R為ns×ns相關系數(shù)矩陣,其元素由相關函數(shù)R(x(i),x(j))計算得到,R(x(i),x(j))對模擬精度起關鍵作用.文中采用高斯函數(shù)求解R(x(i),x(j)),即

(10)

(11)

(12)

(13)

3.2　多目標優(yōu)化算法

文中采用改進的非支配排序方法NSGA-Ⅱ[12]進行多目標優(yōu)化.多目標優(yōu)化問題的最優(yōu)解不再是給定約束條件下使所有目標函數(shù)最大或最小,而是Pareto最優(yōu)解集.Pareto解也叫非支配解,當

(14)

時,解x0支配或優(yōu)于x1(x0

(15)

Pareto最優(yōu)前沿指所有Pareto最優(yōu)解對應的目標函數(shù)空間的區(qū)域PF,

(16)

NSGA-Ⅱ最突出的特征是采用了非支配排序和適應度共享,前者在優(yōu)化過程中收斂到Pareto最優(yōu)前沿,后者則保證了Pareto最優(yōu)解的多樣性.由于文中的目標函數(shù)之間存在矛盾,不可能使各個目標函數(shù)同時達到最優(yōu)解,因此優(yōu)化問題的最終解是從所有Pareto最優(yōu)解中選擇一個最優(yōu)折衷解.

3.3　設計變量與目標函數(shù)

如圖3所示,將控制網(wǎng)格變形的控制點作為設計變量,文中按層選取控制點,對應的徑向比例系數(shù)分別為0.00、0.25、0.60、1.00.將16個控制點作為設計變量,由于網(wǎng)格模型是周期對稱的,為保證周期網(wǎng)格交界面在變形后仍能保持一致,必須將控制點成對分組,將靠近葉片前緣的控制點劃為組a、c、e、g,葉片尾緣的控制點劃為組b、d、f、h.葉片根部控制點組g、h固定不動,控制點組a、b、c、d、e、f分別沿轉子葉片的軸向和切向變化以控制網(wǎng)格變形,因此共有12個有效設計變量.與考慮葉型截面及積疊線的傳統(tǒng)優(yōu)化設計方案(文獻[13]使用了23個設計變量,文獻[14]使用了18個設計變量)相比,文中方案的設計變量明顯減少,計算時間相應減少,從而縮短優(yōu)化設計周期.

圖3　控制點的分組Fig.3　Grouping of control points

葉片氣動性能的多目標優(yōu)化設計,以等熵效率η最大為目標函數(shù)、總壓比r1變化不超過0.5%為約束條件進行優(yōu)化.

目標函數(shù)為

(17)

式中，γ是設計變量,rref是參考總壓比.

4結果分析

4.1　數(shù)值計算的有效性

利用Denton[15]給出的試驗數(shù)據(jù)校驗數(shù)值計算的精度和可靠性.湍流模型為Spalart-Allmaras模型,并且保證數(shù)值結果的網(wǎng)格無關性,計算堵塞質量流量是20.87kg/s,在試驗堵塞質量流量(20.93±0.14)kg/s的偏差范圍內.計算和試驗得到的轉子特性曲線如圖4所示，其中r2為流量比.計算得到的總壓比在分布趨勢和數(shù)值上與試驗值都有很好的一致性,等熵效率在分布趨勢上與試驗值一致,但在數(shù)值上有差別,峰值等熵效率比試驗值低近2%,但這些差別與文獻[15]中給出的仿真計算結果基本一致,試驗失速點流量比約為0.919 0,數(shù)值計算失速點流量比為0.913 5,相差0.59%,因此文中的數(shù)值計算結果是可信的.相比其他湍流模型,Spalart-Allmaras(S-A)模型在保證計算精度的同時可獲得較高的計算效率,針對優(yōu)化問題使用Spalart-Allmaras模型是較好的選擇.

圖4　轉子特性的數(shù)值計算值和試驗值對比Fig.4　Comparison of rotor performance between numerical calculation and experimental values

4.2　優(yōu)化結果分析

在Rotor37的優(yōu)化設計中,Benini[13]的優(yōu)化使轉子在特定工況(總壓比不變)下的等熵效率提高了1.5%;王祥鋒等[14]的優(yōu)化使轉子等熵效率提高了0.7%,總壓比提高了0.66%;Samad等[16]的研究使轉子等熵效率提高了0.51%,總壓比提高了1.25%.

文中采用12個設計變量,試驗設計(DOE)方法采用改進的拉丁超立方抽樣法,即最優(yōu)空間填充設計,共生成281個樣本點.在98%的堵塞質量流量工況下,采用多目標遺傳算法得到的Pareto線如圖5所示.圖5表明,隨著等熵效率的提高,進出口總壓比逐漸降低.要確定目標函數(shù)最優(yōu)設計點時,需要先從Pareto分布中挑選折衷的優(yōu)化候選點,并對每個設計點分別進行CFD計算,再進行分析比較,進而確定最優(yōu)設計點.在98%的堵塞質量流量計算工況下,等熵效率從優(yōu)化前的0.863 57變?yōu)?.870 73,增加了0.83%,總壓比從優(yōu)化前的2.082變?yōu)?.089,增加了0.33%.

圖5　多目標遺傳算法的Pareto前沿分布Fig.5　Pareto front distribution of multi-objective genetic algorith

圖6對比了優(yōu)化前、后葉片的外形變化,優(yōu)化后的葉片形狀可以看作是前傾型葉片.在不同設計轉速(100%、80%、60%)下優(yōu)化前、后轉子的性能曲線如圖7所示,從圖可以看出,在所有工況下轉子的等熵效率都有實質性的提升.以100%設計轉速為例,在95%~99%的堵塞質量流量工況范圍內,轉子效率平均提高近0.7%,而總壓比變化不明顯.在最大等熵效率工況下出口轉子性能參數(shù)沿徑向的分布曲線如圖8所示.優(yōu)化后轉子葉片的等熵效率在20%葉高以上得到提升,其總壓比在30%~85%葉高之間有所提高.出口轉子性能的變化是由于葉片曲率的改變引起三維激波結構發(fā)生了變化,而流場的激波現(xiàn)象對轉子的局部性能有主要的影響.

圖6　優(yōu)化前、后葉片形狀對比Fig.6　Comparison of blade profile between before and after optimization

圖7　優(yōu)化前、后轉子的性能對比Fig.7　Comparison of performance between original and optimized rotor

圖8　最大效率工況下出口氣動參數(shù)展向分布Fig.8　Aerodynamic performance radial plots at outlet at the peak efficiency

圖9　最大效率工況下馬赫數(shù)等值線分布對比Fig.9　Comparison of Mach-number isoline distributions at peak efficiency

在最大等熵效率工況下葉片通道內馬赫數(shù)等值線分布如圖9所示,葉片型面的變化對激波空間形態(tài)的改變有明顯的作用,而各性能參數(shù)及流場結構的變化主要來源于前緣激波的結構變化.在50%葉高處激波并沒有明顯變化,在70%葉高處激波變?yōu)楹舐蛹げ?相應的激波損失減小,流道內流體馬赫數(shù)從原來的0.9增加到1.1;在90%葉高處激波是后掠激波,前緣激波無論在強度還是在波形上都發(fā)生了明顯變化,流道內流體馬赫數(shù)由1.0增加到1.2,由強正激波變?yōu)槿跣奔げ?阻力減少,因此提升了轉子效率.

圖10顯示了葉片優(yōu)化前、后的壁面極限流線.原葉片因激波和吸力面邊界層之間的相互作用而在葉片型面接近中間處產生一條分離線,在邊界層分離區(qū)域靠近葉根的部位分離線之后出現(xiàn)一條附著線,在附著線之后又出現(xiàn)了分離線,隨后在葉片上部區(qū)域靠近尾緣的附近出現(xiàn)附著線,說明原葉片下游流場復雜,出現(xiàn)二次分離,邊界層分離嚴重,分離泡多,損失大.優(yōu)化后葉片的壁面極限流線與原葉片相比,分離線在50%葉高以上開始明顯向吸力面下游傾斜,說明激波由正激波逐漸變?yōu)樾奔げ?分離區(qū)域變小,并且在葉片上部靠近尾緣的區(qū)域沒有出現(xiàn)附著線,說明分離區(qū)域內的流場損失比原葉片小,因此優(yōu)化后葉片性能得到了提高.

圖10　最大效率工況下壁面極限流線對比Fig.10　Comparison of the wall surface limiting streamline at the peak efficiency

5結論

文中針對NASARotor37葉片的優(yōu)化獲得以下結果:

(1)運用網(wǎng)格變形技術結合響應面模型及多目標優(yōu)化算法對壓氣機轉子葉片進行氣動優(yōu)化,葉片性能均得到一定的提升;

(2)基于自由變形控制體的網(wǎng)格變形方法的計算效率高,變形后的網(wǎng)格質量好,節(jié)省了重新劃分網(wǎng)格的時間,同時相比葉型截面及積疊線的優(yōu)化設計方法,輸入變量的數(shù)量明顯減少,縮短了優(yōu)化設計周期;

(3)在98%的堵塞質量流量工況下,優(yōu)化后轉子的等熵效率提高了0.83%,總壓比提高了0.33%,優(yōu)化后的葉型為前傾葉片,結果表明前傾型葉片能顯著改善葉片前緣的激波結構,削弱激波強度,降低激波損失,提升轉子性能.

參考文獻:

[1]LianY,LiouMS.Multi-objectiveoptimizationoftranso-niccompressorbladeusingevolutionaryalgorithm[J].JournalofPropulsionandPower,2005,21(6):979-987.

[2]OkuiH,VerstraeteT,vandenBraembusscheRA,etal.Three-dimensionaldesignandoptimizationofatransonicrotorinaxialflowcompressors[J].JournalofTurbomachinery,2013,135(3):031009.

[3]李彬,宋立明,李軍,等.長葉片透平級多學科多目標優(yōu)化設計 [J].西安交通大學學報,2014,48(1):1- 6.

LiBin,SongLi-ming,LiJun,etal.Multidisciplinaryandmulti-objectiveoptimizationdesignoflongbladeturbinestage[J].JournalofXi’anJiaotongUniversity,2014,48(1):1- 6.

[4]張湘玉,廖文和,劉浩.保細節(jié)的基于曲面控制的網(wǎng)格變形術 [J].華南理工大學學報:自然科學版,2010, 38(7):39- 44.

ZhangXiang-yu,LiaoWen-he,LiuHao.Detailpreservingmeshdeformationbasedonsurfacecontrol[J].JournalofSouthChinaUniversityofTechnology:NaturalScienceEdition,2010,38(7):39- 44.

[5]王元元,張彬乾,郭兆電.基于FFD技術的大型運輸機上翹后體氣動優(yōu)化設計 [J].航空學報,2013,34(8):1806-1814.

WangYuan-yuan,ZhangBin-qian,GuoZhao-dian.AerodynamicoptimizationdesignforlargeupsweptafterbodyoftransportaircraftbasedonFFDtechnology[J].ActaAeronauticaetAstronautiaSinica,2013,34(8):1806-1814.

[6]YusukeT,DanieleP,EmilioFC,etal.Singleandmultiobjectivedesignoptimizationofafastmultihullship:numericalandexperimentalresults[J].JournalofMarineScienceandTechnology,2011,16(4):412- 433.

[7]BiancoliniME,ViolaIM,RiotteM.SailstrimoptimisationusingCFDandRBFmeshmorphing[J].Computers&Fluids,2014,93(10):46-60.

[8]BeniniE,BiolloR.Aerodynamicsofsweptandleanedtransoniccompressor-rotors[J].AppliedEnergy,2007, 84(10):1012-1027.

[9]JangCM,SamadA,KimKY.Optimaldesignofswept,leanedandskewedbladesinatransonicaxialcompressor[C]∥ASMETurboExpo2006:PowerforLand,Sea,andAir.Barcelona：AmericanSocietyofMechanicalEngineers,2006:1279-1288.

[10]BiolloR,BeniniE.Recentadvancesintransonicaxialcompressoraerodynamics[J].ProgressinAerospaceSciences,2013,56:1-18.

[11]SimpsonTW,MaueryTM,KorteJJ,etal.Krigingmodelsforglobalapproximationinsimulation-basedmultidisciplinarydesignoptimization[J].AIAAJournal, 2001,39(12):2233-2241.

[12]DebK,PratapA,AgarwalS,etal.Afastandelitistmultiobjectivegeneticalgorithm:NSGA-II[J].IEEETransactionsonEvolutionaryComputation,2002,6(2):182-197.

[13]BeniniE.Three-dimensionalmulti-objectivedesignoptimizationofatransoniccompressorrotor[J].JournalofPropulsionandPower,2004,20(3):559-565.

[14]王祥鋒,韓萬金.跨音壓氣機轉子多目標氣動優(yōu)化設計 [J].吉林大學學報:工學版,2010,40(1):299-304.

WangXiang-feng,HanWan-jin.Multi-objectiveaerodynamicoptimaldesignforatransoniccompressorrotor[J].JournalofJilinUniversity:EngineeringandTechnologyEdition,2010,40(1):299-304.

[15]DentonJD.Lessonsfromrotor37 [J].JournalofThermalScience,1997,6(1):1-13.

[16]SamadA,KimKY.Shapeoptimizationofanaxialcompressorbladebymulti-objectivegeneticalgorithm[J].ProceedingsoftheInstitutionofMechanicalEngineers,PartA:JournalofPowerandEnergy,2008,222(6):599- 611.

Multi-Objective Aerodynamic Optimum Design of a Blade Based on Mesh Deformation

LiuHaoZhangLeiLiXiao-lin

(School of Mechanical Science and Engineering, Jilin University, Changchun 130025, Jilin, China)

Abstract:There exist the defects of multi design variables and long optimization cycle in the traditional optimum design of compressor blades. In order to solve these problems, this paper proposes a free-form mesh deformation parameterization method of fluid grids by using a non-uniform rational B-spline basis function, and constructs an design system of aerodynamic optimization of rotor blades by combining the advanced design of Latin hypercube sampling experiments, the Kriging response surface model and the NSGA-Ⅱ multi-objective genetic algorithm. Calculation results show that (1) the optimized blade has a total pressure ratio improvement by 0.33% and an isentropic efficiency improvement by 0.83% at a choke mass flow of 98%; (2) the optimized blade is a forward-leaned blade, which helps reduce the shock loss and improve the performance of the rotor; and (3) in comparison with the traditional optimization design method, the proposed optimization system reduces design variables and shortens optimized cycles.

Key words:compressor; optimum design; mesh deformation; response surface; multi-objective genetic algorithm

中圖分類號：V232.4

doi：10.3969/j.issn.1000-565X.2015.09.022

作者簡介：劉浩(1985-),男,博士生,主要從事葉片氣動優(yōu)化設計研究.E-mail: adj1212@163.com? 通信作者： 張雷(1970-),男,教授,博士生導師,主要從事先進加工制造技術研究.E-mail: zhanglei@jlu.edu.cn

*基金項目：國家“863計劃”項目(2012AA041304)；高等學校博士學科點專項科研基金資助項(20110061110022)

收稿日期：2014-10-28

文章編號：1000-565X(2015)09-0141-08

Foundation items： Supported by the National High-Tech R & D Program of China(2012AA041304) and the Specialized Research Fund for the Doctoral Program of Higher Education(20110061110022)