孫 巖, 黃 勇, 王運(yùn)濤,*, 孟德虹, 王 昊
(1. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 四川 綿陽(yáng) 621000; 2. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所, 四川 綿陽(yáng) 621000)
TRIP軟件的靜氣動(dòng)彈性計(jì)算模塊開發(fā)及精度驗(yàn)證
孫 巖1,2, 黃 勇2, 王運(yùn)濤2,*, 孟德虹2, 王 昊2
(1. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 四川 綿陽(yáng) 621000; 2. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所, 四川 綿陽(yáng) 621000)
在結(jié)構(gòu)網(wǎng)格流場(chǎng)求解軟件TRIP基礎(chǔ)上,發(fā)展了一套靜氣動(dòng)彈性計(jì)算模塊。首先簡(jiǎn)要介紹了靜氣動(dòng)彈性計(jì)算模塊總體架構(gòu)、構(gòu)成單元和耦合策略,然后詳細(xì)介紹了構(gòu)成單元中采用的一些數(shù)值計(jì)算方法。最后利用大展弦比機(jī)翼、DLR-F6翼身組合體模型和HIRENASD機(jī)翼模型對(duì)靜氣動(dòng)彈性計(jì)算模塊進(jìn)行了測(cè)試和精度驗(yàn)證。大展弦比機(jī)翼計(jì)算結(jié)果表明柔度矩陣、模態(tài)疊加和有限元數(shù)值求解在模態(tài)選取合適的情況下能夠得到一致的變形計(jì)算結(jié)果。DLR-F6和HIRENASD模型不同計(jì)算軟件的結(jié)果一致性表明本文靜氣動(dòng)彈性計(jì)算模塊采用了正確的算法流程,而計(jì)算與試驗(yàn)的結(jié)果一致性表明靜氣動(dòng)彈性計(jì)算模塊具有很好的預(yù)測(cè)精度。
TRIP;靜氣動(dòng)彈性;柔度矩陣;弱耦合;HIRENASD;DLR-F6
為了減輕結(jié)構(gòu)重量、提高飛行機(jī)動(dòng)性,飛行器設(shè)計(jì)中傾向采用柔性的結(jié)構(gòu)與材料,使得飛行器氣動(dòng)彈性效應(yīng)日益明顯,需要更加精確的氣動(dòng)彈性數(shù)值模擬工具用于飛行器的精細(xì)氣動(dòng)設(shè)計(jì)[1]。
此外,風(fēng)洞模型的完全剛性假設(shè)在新型先進(jìn)飛機(jī)設(shè)計(jì)中受到了挑戰(zhàn)。越來越多的研究表明機(jī)翼在風(fēng)洞氣動(dòng)載荷作用下的變形對(duì)模型氣動(dòng)特性存在明顯的影響,而這種影響在增壓試驗(yàn)中表現(xiàn)的更為顯著[2]。例如波音的超聲速運(yùn)輸機(jī)模型在美國(guó)蘭利的NTF風(fēng)洞開展高雷諾數(shù)風(fēng)洞試驗(yàn),翼梢扭轉(zhuǎn)角可以達(dá)到-2.2°,由此造成不同速壓下的巡航阻力系數(shù)差異最大接近50個(gè)阻力單位[3]。風(fēng)洞模型靜彈性變形對(duì)氣動(dòng)特性的影響已經(jīng)引起了全世界研究者的關(guān)注與重視[4],高效高精度的靜氣動(dòng)彈性計(jì)算程序逐漸成為研究這一問題的重要手段。
為了提高計(jì)算結(jié)果的可信度水平,國(guó)外廣泛開展了氣動(dòng)彈性計(jì)算軟件的驗(yàn)證與確認(rèn)研究,其中影響范圍較大的是AIAA發(fā)起的氣動(dòng)彈性預(yù)測(cè)研討會(huì)AePW(Aeroelastic Prediction Workshop)[5],得到了世界范圍內(nèi)大量學(xué)者的積極響應(yīng)。
近些年,隨著計(jì)算機(jī)硬件技術(shù)的快速發(fā)展,國(guó)內(nèi)高校與科研機(jī)構(gòu)基于各自求解器開發(fā)了多種氣動(dòng)彈性數(shù)值計(jì)算平臺(tái)[6-8],并廣泛應(yīng)用于解決工程實(shí)際問題,但靜氣動(dòng)彈性數(shù)值計(jì)算程序精度驗(yàn)證方面的系統(tǒng)研究還比較少。
本文在課題組開發(fā)的流場(chǎng)計(jì)算軟件TRIP上發(fā)展了一套高效的靜氣動(dòng)彈性計(jì)算模塊,下文將簡(jiǎn)要介紹計(jì)算模塊的架構(gòu)、組成、流程和采用的數(shù)值方法。然后通過2個(gè)算例對(duì)計(jì)算模塊的精度進(jìn)行驗(yàn)證。
圖1給出了TRIP軟件靜氣動(dòng)彈性計(jì)算模塊的總體架構(gòu),該模塊包含四個(gè)主要功能單元:氣動(dòng)數(shù)值計(jì)算、結(jié)構(gòu)靜變形計(jì)算、耦合界面數(shù)據(jù)插值以及動(dòng)態(tài)網(wǎng)格方法,四個(gè)主要功能單元采用弱耦合方式串聯(lián)在一起,如圖1所示。
圖1的左側(cè)給出了氣動(dòng)載荷和結(jié)構(gòu)位移的傳遞路線。首先CFD計(jì)算出氣動(dòng)載荷,利用界面數(shù)據(jù)插值將氣動(dòng)載荷傳遞給CSM,CSM利用氣動(dòng)輸入載荷計(jì)算得到結(jié)構(gòu)的位移,通過界面數(shù)據(jù)插值再將結(jié)構(gòu)位移傳遞給流場(chǎng),并利用動(dòng)態(tài)網(wǎng)格技術(shù)生成流場(chǎng)物面變形后的新網(wǎng)格,CFD再計(jì)算新網(wǎng)格下的流場(chǎng)。如此反復(fù)迭代,直至載荷和位移均達(dá)到收斂。圖1的右側(cè)給出了每個(gè)功能單元采用的主要計(jì)算方法,下節(jié)將對(duì)這些計(jì)算方法進(jìn)行更加詳細(xì)的介紹。
2.1流場(chǎng)計(jì)算方法
流場(chǎng)解算器采用課題組開發(fā)多年的亞跨超流場(chǎng)解算軟件TRIP。TRIP軟件通過有限體積法離散求解歐拉(Euler)和雷諾平均方程(RANS),方程對(duì)流項(xiàng)離散采用二階MUSCL型的Roe通量差分分裂格式,粘性項(xiàng)離散采用二階中心差分格式,湍流模型采用工程常用的單方程SA模型[9]與兩方程SST模型[10],計(jì)算加速方面采用多重網(wǎng)格和基于MPI信息傳遞的并行計(jì)算。課題組針對(duì)TRIP已經(jīng)開展了大量的算例驗(yàn)證與測(cè)試,軟件的精度和可信度得到了廣泛的認(rèn)可[11-12]。
2.2結(jié)構(gòu)靜變形計(jì)算方法
結(jié)構(gòu)靜變形計(jì)算采用柔度矩陣、模態(tài)疊加和有限元數(shù)值求解三種方法。
柔度矩陣方法中,結(jié)構(gòu)網(wǎng)格點(diǎn)的位移通過下式計(jì)算得到:
式中,us為結(jié)構(gòu)網(wǎng)格點(diǎn)的位移矩陣,C為結(jié)構(gòu)的柔度矩陣,fs為結(jié)構(gòu)網(wǎng)格點(diǎn)上的輸入載荷。本文中結(jié)構(gòu)柔度矩陣C通過有限元模型施加單位載荷的位移響應(yīng)組裝得到。
模態(tài)疊加法中,結(jié)構(gòu)網(wǎng)格點(diǎn)的位移通過下式計(jì)算得到[13]:
式中,us含義同式(1),Φs為結(jié)構(gòu)的廣義模態(tài)矩陣,通過結(jié)構(gòu)模態(tài)分析得到,q為系統(tǒng)的廣義位移矢量。對(duì)于結(jié)構(gòu)靜變形有:
式中,K為結(jié)構(gòu)的剛度矩陣,fs含義同式(1),式(2)代入式(3),可以得到:
式中:
式中,Ω為結(jié)構(gòu)模態(tài)剛度矩陣,是一個(gè)對(duì)角陣,對(duì)角線各元素為廣義模態(tài)圓頻率的平方。式(4)代入式(2)即可求解得到結(jié)構(gòu)網(wǎng)格點(diǎn)的位移:
有限元數(shù)值求解采用了自主編寫的一套結(jié)構(gòu)靜變形求解程序,在方法上借鑒了商業(yè)軟件Nastran的數(shù)據(jù)格式和眾多開源程序的設(shè)計(jì)架構(gòu)。
2.3耦合界面數(shù)據(jù)傳遞
氣動(dòng)/結(jié)構(gòu)之間的數(shù)據(jù)傳遞可以通過一個(gè)矩陣進(jìn)行描述,比如對(duì)于位移變量可以表示為:
式中,u為網(wǎng)格點(diǎn)位移矢量,下標(biāo)a表示氣動(dòng),下標(biāo)s表示結(jié)構(gòu),H為轉(zhuǎn)換矩陣。為了滿足傳遞過程中的能量守恒,有虛功原理可以得到:
式中,f為網(wǎng)格點(diǎn)上的載荷。
本文靜氣動(dòng)彈性計(jì)算模塊采用RBF、TPS和IPS三種方法構(gòu)建傳遞矩陣H。一般形式的徑向基函數(shù)插值可以表示為[14]:
式中,f(x)為被插函數(shù),x為插值點(diǎn)位置坐標(biāo),xi為插值基點(diǎn)的位置坐標(biāo),xj(j=1,…,m)為位置坐標(biāo)的分量,‖x-xi‖表示兩點(diǎn)之間的距離,N為插值基點(diǎn)的數(shù)量,m為描述問題的維數(shù),二維m為2,三維m為3,βi(i=0,…,m)、γi(i=1,…,N)為插值系數(shù),φ為徑向基函數(shù)。
利用RBF方法得到的傳遞矩陣H可以表示為:
其中:
式(11)中的M為氣動(dòng)網(wǎng)格點(diǎn)的數(shù)目,式(12)中右側(cè)的0為(m+1)×N的零矩陣,I為N×N的單位矩陣。
不同的基函數(shù)φ可以得到不同的傳遞矩陣H,TPS和IPS是RBF插值的一種特殊形式,對(duì)應(yīng)不同的基函數(shù)類型,基函數(shù)的選取可以參考文獻(xiàn)[14]。
2.4動(dòng)態(tài)網(wǎng)格生成
本文動(dòng)態(tài)網(wǎng)格生成采用網(wǎng)格變形的思路,可以避免引起新的數(shù)值離散誤差。目前計(jì)算模塊主要采用面向結(jié)構(gòu)網(wǎng)格的RBF+TFI方法,針對(duì)一些特殊問題,也采用RBF+Delaunay和RBF+HBGM的網(wǎng)格變形方法,有關(guān)這些方法可以參考文獻(xiàn)[15-16]。
本節(jié)將首先通過一個(gè)大展弦比機(jī)翼的變形結(jié)果分析三種不同結(jié)構(gòu)靜變形計(jì)算方法的差異,然后利用兩個(gè)風(fēng)洞模型的試驗(yàn)數(shù)據(jù)結(jié)果對(duì)靜氣動(dòng)彈性計(jì)算模塊精度進(jìn)行驗(yàn)證。下文如不作特殊說明,流場(chǎng)計(jì)算中湍流模擬采用兩方程SST模型,耦合數(shù)據(jù)傳遞采用TPS方法,動(dòng)態(tài)網(wǎng)格變形采用RBF+TFI方法。
3.1大展弦比機(jī)翼變形結(jié)果分析
本小節(jié)將分析三種不同結(jié)構(gòu)靜變形求解方法結(jié)果的差異,圖2給出了采用的大展弦比機(jī)翼的外形、網(wǎng)格和有限元模型。其中有限元模型采用不同厚度的殼單元構(gòu)建,因此耦合數(shù)據(jù)傳遞采用IPS方法。
圖3給出了三種不同方法計(jì)算得到的大展弦比機(jī)翼變形結(jié)果,圖中Moden表示采用結(jié)構(gòu)前n階模態(tài)位移的變形計(jì)算結(jié)果。
圖3可以看出:柔度矩陣方法與有限元分析結(jié)構(gòu)一致性很好,模態(tài)疊加方法的變形結(jié)果與選取的模態(tài)有關(guān)。但隨著模態(tài)數(shù)量的增加,模態(tài)方法的結(jié)果與前兩種方法的結(jié)果也很快趨向一致。另一方面,不同模態(tài)對(duì)彎曲變形和扭轉(zhuǎn)變形的貢獻(xiàn)不同,對(duì)于彎曲變形,前2階的結(jié)構(gòu)模態(tài)貢獻(xiàn)了主要的變形,而對(duì)于扭轉(zhuǎn)變形,后續(xù)模態(tài)仍然有著明顯的貢獻(xiàn)。因此,在選擇結(jié)構(gòu)模態(tài)計(jì)算結(jié)構(gòu)靜變形時(shí),應(yīng)當(dāng)首先開展結(jié)構(gòu)模態(tài)對(duì)靜變形計(jì)算結(jié)果的影響分析。如不作特殊說明,本文后續(xù)的算例中均采用柔度矩陣方法計(jì)算結(jié)構(gòu)靜變形。
3.2DLR-F6模型靜氣動(dòng)彈性模擬
DLR-F6是德國(guó)宇航院設(shè)計(jì)的一類現(xiàn)代民用運(yùn)輸機(jī)標(biāo)模,在多座風(fēng)洞開展了不同類型的風(fēng)洞試驗(yàn),具有豐富的試驗(yàn)數(shù)據(jù),有關(guān)DLR-F6模型的詳細(xì)可以參考文獻(xiàn)[17]。本文采用DLR-F6的翼身組合體構(gòu)型開展精度驗(yàn)證,圖4給出了該構(gòu)型的外形、網(wǎng)格拓?fù)浜陀邢拊P汀?/p>
流場(chǎng)計(jì)算Ma數(shù)為0.75,基于平均氣動(dòng)弦長(zhǎng)的雷諾數(shù)Re為3.0×106,模型迎角為-3°~1.5°,來流速壓q與模型材料彈性模量E的比值q/E為2.0×10-7。圖5給出了該計(jì)算參數(shù)條件下DLR-F6翼身組合體模型機(jī)翼變形展向分布結(jié)果。圖中Test表示在NASA的NTF風(fēng)洞利用光學(xué)方法測(cè)量的機(jī)翼變形,TAU為DLR的流場(chǎng)解算程序。
可以看出:不同升力系數(shù)下,試驗(yàn)Test、TAU和本文Present得到的彎曲變形結(jié)果在變形分布和量值上都比較一致,表明本文靜氣動(dòng)彈性計(jì)算模塊采用了正確的算法流程且具有較好的彎曲變形預(yù)測(cè)精度。由于DLR-F6是實(shí)體金屬模型,來流的速壓并不大,因此由氣動(dòng)載荷引起的機(jī)翼最大彎曲變形并不明顯,約為6.5 mm(半翼展為585.6 mm)。
3.3HIRENASD模型靜氣動(dòng)彈性模擬
HIRENASD是亞琛工大設(shè)計(jì)的一類典型運(yùn)輸機(jī)機(jī)翼外形,在歐洲跨聲速風(fēng)洞ETW開展了大量試驗(yàn),有關(guān)HIRENASD模型的詳細(xì)參數(shù)可以參考文獻(xiàn)[18]。為了消除壁面邊界層對(duì)機(jī)翼流場(chǎng)的影響,風(fēng)洞中給HIRENASD機(jī)翼模型安裝了一個(gè)假機(jī)身,本文采用包含假機(jī)身的外形開展靜氣動(dòng)彈性模擬。圖6給出了HIRENASD機(jī)翼模型的外形、網(wǎng)格拓?fù)浜陀邢拊P汀F渲辛鲌?chǎng)計(jì)算網(wǎng)格從氣動(dòng)彈性預(yù)測(cè)會(huì)議AePW的官網(wǎng)上下載得到,為ICEM軟件生成的多塊對(duì)接網(wǎng)格,如圖6(a);有限元的固支約束添加在機(jī)翼與天平連接的端面上,如圖6(b)。
流場(chǎng)計(jì)算的Ma數(shù)為0.80,基于平均氣動(dòng)弦長(zhǎng)的雷諾數(shù)Re為1.4×107,模型迎角α為3°,來流速壓q與模型材料彈性模量E的比值q/E為4.7×10-7。
圖7給出了該參數(shù)條件下的機(jī)翼彎曲變形展向分布結(jié)果,可以看出,試驗(yàn)Test、TAU和本文Present得到的機(jī)翼變形結(jié)果十分吻合。
圖8給出了HIRENASD機(jī)翼展向相對(duì)位置η等于0.95時(shí)的翼剖面壓力系數(shù)分布,可以看出,采用剛性外形,計(jì)算的壓力系數(shù)分布與試驗(yàn)結(jié)果之間存在顯著的差異,而采用靜氣動(dòng)彈性耦合計(jì)算得到的壓力系數(shù)分布與試驗(yàn)結(jié)果十分吻合,進(jìn)一步驗(yàn)證了本文靜氣動(dòng)彈性計(jì)算模塊的精度。
本文在TRIP軟件基礎(chǔ)發(fā)展了靜氣動(dòng)彈性計(jì)算模塊,并通過算例對(duì)計(jì)算結(jié)果的精度進(jìn)行了驗(yàn)證,說明該計(jì)算模塊具有較好的靜氣動(dòng)彈性預(yù)測(cè)精度。同時(shí),計(jì)算結(jié)果表明風(fēng)洞模型結(jié)構(gòu)靜變形對(duì)壓力系數(shù)分布有著明顯影響,應(yīng)當(dāng)?shù)玫皆囼?yàn)數(shù)據(jù)使用單位的重視。
致謝:本文得到了課題組張玉倫、洪俊武、王光學(xué)、張書俊、李偉和楊小川等人的幫助,在此表示感謝。
[1]Mouton S, Sant Y, Lyonnet M. Predication of the aerodynamic effect of model deformation during transonic wind tunnel tests[J]. International Journal of Engineering Systems Modelling and Simulation, 2013, 5(1-3): 44-56.
[2]孫巖, 張征宇, 鄧小剛, 等. 風(fēng)洞模型靜彈性變形對(duì)氣動(dòng)力影響研究[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2013, 31(3): 294-300.
[3]Owens L R, Wahls R A, Rivers S M. Off-design Reynolds number effects for a supersonic transport[J]. Journal of Aricraft, 2005, 42(6): 1427- 441.
[4]Levy D W, Laflin K R, Tinoco E N, et al. Summary of data from the fifth computational fluid dynamics drag prediction workshop[J]. Journal of Aircraft, 2014, 51(4): 1194-1213.
[5]Heeg J, Chwalowski P, Florance J P, et al. Overview of the aeroelastic predication workshop[R]. AIAA 2013-0783, 2013.
[6]楊超, 王立波, 謝長(zhǎng)川, 等. 大變形飛機(jī)配平與飛行載荷分析方法[J]. 中國(guó)科學(xué)(E輯: 技術(shù)科學(xué)), 2012, 42(10): 1137-1147.
[7]徐敏, 邱滋華, 裴曦. 基于CFD/CSD/CAA耦合的聲氣動(dòng)彈性仿真技術(shù)在壁板顫振中的應(yīng)用[J]. 強(qiáng)度與環(huán)境, 2013, 40(5): 16-24.
[8]黃煒, 陸志良, 唐迪, 等. 基于多點(diǎn)約束的大展弦比機(jī)翼靜氣動(dòng)彈性計(jì)算[J]. 北京航空航天大學(xué)學(xué)報(bào), 2014, 40(12): 1666-1671.
[9]Spalart P R, Allmaras S R. A one-equation turbulence model for aerodynamic flows[R]. AIAA-92-0439, 1992.
[10]Menter F R. Two-equation eddy-viscosity turbulence models for engineering application[J]. AIAA Journal, 1994, 32(8): 1598-1605.
[11]王運(yùn)濤, 王光學(xué), 張玉倫. DPWⅢ機(jī)翼和翼身組合體構(gòu)型數(shù)值模擬[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2011, 29(3): 264-269.
[12]王運(yùn)濤, 張書俊, 孟德虹. DPW4翼/身/平尾組合體的數(shù)值模擬[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2013, 31(6): 739-744.
[13]Ritter M. Static and forced motion aeroelastic simulations of the HIRENASD wind tunnel model[R]. AIAA 2012-1633, 2012.
[14]Rendall T C S, Allen C B. Unified fluid-structure interpolation and mesh motion using radial basis functions[J]. International Journal for Numerical Methods in Engineering, 2008, 74: 1519-1559.
[15]孫巖, 鄧小剛, 王運(yùn)濤, 等. RBF_TFI結(jié)構(gòu)動(dòng)網(wǎng)格技術(shù)在風(fēng)洞靜氣動(dòng)彈性修正中的應(yīng)用[J]. 工程力學(xué), 2014, 31(10): 228-233.
[16]孫巖, 鄧小剛, 王光學(xué), 等. 基于徑向基函數(shù)改進(jìn)的Delaunay圖映射動(dòng)網(wǎng)格方法[J]. 航空學(xué)報(bào), 2014, 35(3): 727-735.
[17]Burner A W, Goad W K, Massey E A, et al. Wind deformation measurements of the DLR-F6 transport configuration in the national transonic facility[R]. AIAA 2008-6921, 2008.
[18]Ritter M, Hassan D. Assesment of the ONERA/DLR numerical aeroelastics predication capabilities on the HIRENASD configuration[R]. IFASD-2011-109, 2011.
DevelopmentandprecisionvalidationofstaticaeroelasticcomputationalmoduleonflowsolverTRIP
(1.StateKeyLaboratoryofAerodynamics,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China; 2.ComputationalAerodynamicsInstituteofChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China)
A static aeroelastic computational module is developed on the flow solver TRIP. Firstly, the frame, main component units and coupling strategy of the static aeroelastic computational module are introduced briefly. Then numerical methods used in the component units are described in detail. Finally, the static aeroelastic computational module is tested and validated through three cases: high respect-ratio wing, DLR-F6 wing-body model and HIRENASD wing model. The computational results of high respect-ratio wing show that flexible matrix method, modal superposition method and finite element method can have same deformation results when the structure modes are selected properly. The result consistency of DLR-F6 or HIRENASD model among different solvers demonstrates that the static aeroelastic computational module has adopted a correct algorithm procedure. The results consistency of DLR-F6 or HIRENASD model between computations and tests shows that the present static aeroelastic computational module has a good prediction precision.
TRIP; static aeroelasticity; flexibility matrix; weak coupling; HIRENASD; DLR-F6
V211.3
A
10.7638/kqdlxxb-2015.0154
0258-1825(2017)05-0620-05
2015-08-26;
2015-10-14
國(guó)家重點(diǎn)研發(fā)計(jì)劃(2016YFB0200700)
孫巖(1986-),男,安徽鳳陽(yáng)人,博士,助理研究員,研究方向:計(jì)算氣動(dòng)彈性力學(xué). E-mail:supersunyan@163.com
王運(yùn)濤*(1967-),男,博士,研究員,博士生導(dǎo)師,主要研究方向:計(jì)算空氣動(dòng)力學(xué). E-mail: ytwang@skla.cardc.cn
孫巖, 黃勇, 王運(yùn)濤, 等. TRIP軟件的靜氣動(dòng)彈性計(jì)算模塊開發(fā)及精度驗(yàn)證[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2017, 35(5): 620-624.
10.7638/kqdlxxb-2015.0154 SUN Y, HUANG Y, WANG Y T, et al. Development and precision validation of static aeroelastic computational module on flow solver TRIP[J]. Acta Aerodynamica Sinica, 2017, 35(5): 620-624.
SUN Yan1,2, HUANG Yong2, WANG Yuntao2,*, MENG Dehong2, WANG Hao2