馮浩洋, 蘇新兵, 馬斌麟, 王旭

(空軍工程大學(xué) 航空航天工程學(xué)院, 陜西 西安 710038)

?

多面體網(wǎng)格在靜氣動(dòng)彈性計(jì)算中的應(yīng)用

馮浩洋, 蘇新兵, 馬斌麟, 王旭

(空軍工程大學(xué) 航空航天工程學(xué)院, 陜西 西安 710038)

介紹了多面體網(wǎng)格技術(shù)在靜氣動(dòng)彈性問(wèn)題中的應(yīng)用。基于M6算例和徑向基函數(shù)插值法,探究了多面體網(wǎng)格的CFD計(jì)算效率和變形能力,使用AGARD445.6機(jī)翼進(jìn)行了靜氣動(dòng)彈性數(shù)值計(jì)算。結(jié)果表明,多面體網(wǎng)格具有計(jì)算效率高、變形能力強(qiáng)、計(jì)算結(jié)果準(zhǔn)確等優(yōu)點(diǎn),可適用于靜氣動(dòng)彈性問(wèn)題的數(shù)值計(jì)算。

多面體網(wǎng)格; 靜氣動(dòng)彈性計(jì)算; 網(wǎng)格效率; 網(wǎng)格變形

0　引言

隨著計(jì)算機(jī)性能的提高和相關(guān)學(xué)科的發(fā)展,基于計(jì)算流體力學(xué)(CFD)和計(jì)算結(jié)構(gòu)力學(xué)(CSD)的流固耦合仿真技術(shù)逐漸應(yīng)用于氣動(dòng)彈性數(shù)值模擬。在流固耦合仿真中,CFD計(jì)算占用了90%以上的時(shí)間,而且伴隨CSD的求解,CFD網(wǎng)格會(huì)產(chǎn)生變形。因此,CFD網(wǎng)格的計(jì)算效率和變形能力是仿真快速、準(zhǔn)確的關(guān)鍵。當(dāng)前,CFD網(wǎng)格主要有四面體、六面體(切割體)和多面體三種類(lèi)型。四面體網(wǎng)格對(duì)模型的自適應(yīng)性好,變形能力強(qiáng),但是計(jì)算量大,容易產(chǎn)生偽耗散。切割體網(wǎng)格生成質(zhì)量高,計(jì)算量小,但變形能力弱[1]。多面體網(wǎng)格作為近年來(lái)新興的一種網(wǎng)格技術(shù),具有生成方便快捷,對(duì)復(fù)雜外形具有良好的適應(yīng)能力等特點(diǎn),得益于其獨(dú)特的網(wǎng)格架構(gòu)和有更多的相鄰單元,數(shù)值模擬可以快速地收斂,梯度的計(jì)算和當(dāng)?shù)氐牧鲃?dòng)狀況預(yù)測(cè)更準(zhǔn)確。由于多面體對(duì)幾何的變形沒(méi)有四面體敏感,因此多面體網(wǎng)格可以更好地保證變形后的網(wǎng)格質(zhì)量[2-3]。

本文對(duì)四面體、切割體、多面體網(wǎng)格的計(jì)算效率進(jìn)行了比較,檢驗(yàn)了多面體網(wǎng)格的變形能力。在文獻(xiàn)[4]方法的基礎(chǔ)上,基于多面體網(wǎng)格和徑向基函數(shù)的網(wǎng)格變形方法,通過(guò)AGARD算例,對(duì)機(jī)翼的靜氣動(dòng)彈性現(xiàn)象進(jìn)行了仿真,驗(yàn)證了多面體網(wǎng)格的有效性。

1　網(wǎng)格計(jì)算效率

為了探究多面體網(wǎng)格的計(jì)算效率,本文分別采用四面體、切割體、多面體網(wǎng)格對(duì)M6機(jī)翼算例進(jìn)行了計(jì)算和對(duì)比。在文獻(xiàn)[2]中,已經(jīng)對(duì)網(wǎng)格總數(shù)不一致的情況進(jìn)行了討論,并得出了多面體網(wǎng)格收斂速度快、計(jì)算精度高的結(jié)論。本文進(jìn)一步對(duì)基于相同網(wǎng)格單元數(shù)的情況進(jìn)行了探討。在調(diào)整面網(wǎng)格尺寸并進(jìn)行網(wǎng)格獨(dú)立性驗(yàn)證后,可知三種網(wǎng)格數(shù)量在1 900 000左右可以滿(mǎn)足算例要求。剖分后的網(wǎng)格單元數(shù)量如下:四面體網(wǎng)格單元數(shù)目為1 984 521,切割體網(wǎng)格單元數(shù)目為1 948 756,多面體網(wǎng)格單元數(shù)目為1 956 546。剖分后的翼根前緣網(wǎng)格如圖1所示。

圖1　翼根前緣網(wǎng)格示意圖Fig.1　The mesh of wing root leading edge

取N-S方程作為流場(chǎng)數(shù)值求解方程。空間離散采用中心差分格式,湍流模型為S-A模型,采用耦合隱式方法進(jìn)行求解。為了加速收斂,在不影響計(jì)算精度的前提下,設(shè)置庫(kù)朗數(shù)為50,松弛因子為0.5。計(jì)算狀態(tài)為:Ma=0.839 5;Re=11.72×106;α=3.06°。經(jīng)計(jì)算,得到如圖2所示的三種網(wǎng)格的機(jī)翼沿不同展向的壓力分布數(shù)值模擬計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)[5]的對(duì)比。

圖2　翼展不同站位的壓力系數(shù)分布Fig.2　Pressure coefficient distributions at different wingspan positions

由圖2可知,三種網(wǎng)格在計(jì)算準(zhǔn)確性上幾乎沒(méi)有差別,切割體網(wǎng)格在y/b=0.65處對(duì)前緣激波的模擬較另兩種網(wǎng)格差,這可能是由于其貼體性較差,在曲率較大的地方難以與幾何外形精確吻合導(dǎo)致的。

三種網(wǎng)格下的計(jì)算殘差曲線如圖3所示。由圖3可以看到,三者的收斂速度大致相同,多面體網(wǎng)格和四面體網(wǎng)格在150步左右收斂,切割體網(wǎng)格在180步左右收斂。四面體網(wǎng)格收斂速度快是因?yàn)槠涔?jié)點(diǎn)數(shù)最少,易于收斂,但其收斂精度最低。多面體網(wǎng)格在與四面體網(wǎng)格收斂速度相當(dāng)?shù)耐瑫r(shí),收斂精度遠(yuǎn)比其高,與切割體網(wǎng)格大致相同,但多面體網(wǎng)格收斂過(guò)程更加平穩(wěn)。

圖3　殘差收斂曲線Fig.3　The residual convergence curves

三種網(wǎng)格的收斂精度和計(jì)算時(shí)間如表1所示。由于切割體網(wǎng)格在計(jì)算過(guò)程中耗散最小,且M6機(jī)翼外形簡(jiǎn)單,因此其所用時(shí)間最短,而多面體網(wǎng)格計(jì)算時(shí)間小于四面體網(wǎng)格。

表1　收斂精度和計(jì)算時(shí)間Table1　Convergenceaccuracyandcomputationtime

由計(jì)算結(jié)果可知,多面體網(wǎng)格在流體計(jì)算過(guò)程中綜合了四面體網(wǎng)格和切割體網(wǎng)格的優(yōu)點(diǎn)。當(dāng)網(wǎng)格數(shù)相同時(shí),多面體網(wǎng)格在收斂速度快的同時(shí)還有較高的精度和平穩(wěn)度。因此,相較于切割體網(wǎng)格和四面體網(wǎng)格,采用多面體網(wǎng)格進(jìn)行靜氣動(dòng)彈性計(jì)算可以提高CFD的計(jì)算效率和準(zhǔn)確性,從而縮短整體計(jì)算時(shí)間。

2　網(wǎng)格變形能力

為了檢驗(yàn)多面體網(wǎng)格的變形能力[6],本文基于徑向基函數(shù)變形方法[7-10],對(duì)M6機(jī)翼的網(wǎng)格進(jìn)行了旋轉(zhuǎn)變形,并對(duì)變形前后機(jī)翼的氣動(dòng)參數(shù)進(jìn)行了計(jì)算和對(duì)比。

采用上文所剖分的多面體網(wǎng)格,以翼根前緣頂點(diǎn)為圓心,z軸為轉(zhuǎn)軸,對(duì)機(jī)翼進(jìn)行45°旋轉(zhuǎn)。為便于觀察,取流場(chǎng)的對(duì)稱(chēng)面網(wǎng)格旋轉(zhuǎn)變形,如圖4所示。對(duì)于航空外流問(wèn)題,計(jì)算結(jié)果的準(zhǔn)確性與邊界層的質(zhì)量息息相關(guān),因此本文著重對(duì)邊界層網(wǎng)格的變形質(zhì)量進(jìn)行分析。為了更好地顯示變形細(xì)節(jié),對(duì)變形前后翼根的前后緣(黑框內(nèi)部分)進(jìn)行局部放大,放大后的網(wǎng)格如圖5所示。

圖4　機(jī)翼旋轉(zhuǎn)示意圖Fig.4　The wing rotation diagram

圖5　翼根切面前后緣網(wǎng)格變形比較Fig.5　Comparison of mesh deformation at wing root section leading and trailing edge

由圖可見(jiàn),在機(jī)翼旋轉(zhuǎn)變形45°的情況下,其前后緣網(wǎng)格都發(fā)生了不同程度的偏斜。后緣偏斜程度較大,這是因?yàn)楹缶壋叽巛^小,其邊界層網(wǎng)格的外邊界與后緣壁面的距離相對(duì)較遠(yuǎn)所致?；趶较蚧瘮?shù)的變形方法對(duì)遠(yuǎn)離壁面網(wǎng)格節(jié)點(diǎn)的控制力較弱,在大變形情況下,網(wǎng)格經(jīng)過(guò)多次迭代變形會(huì)產(chǎn)生較大的誤差,導(dǎo)致后緣處網(wǎng)格偏斜更為嚴(yán)重。

對(duì)于航空外流問(wèn)題,近壁面第一層網(wǎng)格高度應(yīng)當(dāng)滿(mǎn)足湍流模型的計(jì)算要求,本文選取的是S-A模型。為了提高升力、阻力計(jì)算結(jié)果的準(zhǔn)確程度,需要求解粘性子層,這就要保證壁面Y+值小于1[11]。取機(jī)翼y/b=0.5處的截面,在變形前后對(duì)其壁面Y+值進(jìn)行計(jì)算,結(jié)果如圖6所示。由圖6可以看到,該截面的Y+值幾乎沒(méi)有發(fā)生變化,只有后緣處變化幅度相對(duì)較大,這也與網(wǎng)格變形情況相符。圖中變形前后的Y+值均在(0,1)區(qū)間內(nèi),滿(mǎn)足湍流模型的計(jì)算要求。

圖6　y/b=0.5處Y+值變化Fig.6　The change of Y+value on y/b=0.5

依照上節(jié)邊界條件對(duì)變形后的流場(chǎng)網(wǎng)格進(jìn)行計(jì)算,得到如圖7所示的變形前后機(jī)翼沿不同展向的壓力分布。由圖7可知,盡管前緣和后緣處的網(wǎng)格發(fā)生了偏斜,但是計(jì)算結(jié)果仍然是準(zhǔn)確的,可以認(rèn)為局部網(wǎng)格的偏斜對(duì)計(jì)算沒(méi)有造成負(fù)面影響。

圖7　變形前后翼展不同站位的壓力系數(shù)分布Fig.7　Pressure coefficient distributions at different wingspan positions before and after deformation

網(wǎng)格變形前后的收斂曲線如圖8所示。由圖8可知,曲線都在260步左右達(dá)到收斂;但是變形后網(wǎng)格計(jì)算結(jié)果的收斂曲線不如變形前平穩(wěn),而且收斂精度比變形前略差。這是因?yàn)樽冃魏蟮木W(wǎng)格整體質(zhì)量有所下降,致使計(jì)算過(guò)程中出現(xiàn)偽耗散,使得收斂精度降低。盡管如此,由于計(jì)算結(jié)果的準(zhǔn)確性并沒(méi)有受到影響,所以仍然說(shuō)明多面體網(wǎng)格在該變形方法產(chǎn)生大變形后可以滿(mǎn)足計(jì)算要求。

圖8　變形前后殘差收斂曲線Fig.8　The residual convergence curves before and after deformation

3　流體、結(jié)構(gòu)計(jì)算和數(shù)據(jù)傳遞方法

靜氣動(dòng)彈性變形是由飛行器外部載荷與機(jī)體彈性結(jié)構(gòu)相互耦合作用產(chǎn)生的,需要耦合氣動(dòng)力方程和結(jié)構(gòu)靜平衡方程進(jìn)行求解。流體求解模型與上文相同,由于計(jì)算的是靜氣動(dòng)彈性問(wèn)題,所以不需要在時(shí)域上進(jìn)行求解。為了提高效率,對(duì)流體進(jìn)行穩(wěn)態(tài)求解。結(jié)構(gòu)求解采用柔度矩陣法,流體和結(jié)構(gòu)之間用松耦合的方式進(jìn)行數(shù)據(jù)交換,采用無(wú)限平板插值方法(IPS)完成兩者之間的數(shù)據(jù)傳遞[4,12]。

4　算例驗(yàn)證

采用標(biāo)準(zhǔn)的AGARD445.6機(jī)翼對(duì)本文的氣動(dòng)彈性計(jì)算方法進(jìn)行驗(yàn)證。AGARD機(jī)翼為NACA65A004翼型,展弦比為1.65,梢根比為0.66,展長(zhǎng)為0.76 m,根弦長(zhǎng)為0.56 m,弦線后掠角為45°[13]。剖分的流體網(wǎng)格和結(jié)構(gòu)網(wǎng)格如圖9所示,流體網(wǎng)格數(shù)為2 401 689,結(jié)構(gòu)網(wǎng)格數(shù)為9 593。

圖9　AGARD機(jī)翼流體網(wǎng)格和結(jié)構(gòu)網(wǎng)格Fig.9　The fluid mesh and structure mesh of AGARD wing

根據(jù)文獻(xiàn)[12],設(shè)置結(jié)構(gòu)材料屬性為:E1=0.89 GPa,E2=1.54 GPa,ν=0.31,G=2.6 GPa,ρ=381.98 kg/m3。其中,E1為材料x(chóng)和z方向的彈性模量;E2為y方向的彈性模量(x沿弦向,y沿展向);ν為泊松比;G為每個(gè)方向的剪切模量;ρ為機(jī)翼模型的密度。

計(jì)算狀態(tài)選擇Ma=0.8,α=1°。AGARD 445.6機(jī)翼靜氣動(dòng)彈性迭代平衡后,相對(duì)于初始位置,翼梢前緣和后緣的變形量如表2所示。由表2可以看出,本文方法與文獻(xiàn)[14] 的計(jì)算結(jié)果基本一致。變形量有誤差是由于插值方法的限制和網(wǎng)格數(shù)量的不同,數(shù)據(jù)在傳遞過(guò)程中會(huì)有一定的偏差所致。由于結(jié)果相差都在3%以?xún)?nèi),所以這樣的誤差是可以接受的。

表2　AGARD 445.6機(jī)翼靜氣動(dòng)彈性結(jié)果比較Table2　ComparisonofAGARD445.6wingstaticaero-elasticcomputationalresults

5　結(jié)束語(yǔ)

基于多面體網(wǎng)格,本文對(duì)靜氣動(dòng)彈性計(jì)算的效率和網(wǎng)格變形問(wèn)題進(jìn)行了如下三方面的探究:

(1)在流場(chǎng)數(shù)值模擬的基礎(chǔ)上,對(duì)切割體、四面體和多面體網(wǎng)格的計(jì)算效率進(jìn)行了對(duì)比。M6機(jī)翼計(jì)算結(jié)果表明,在基于相同單元數(shù)的情況下,多面體網(wǎng)格相比另外兩種網(wǎng)格可以極大地縮短流體計(jì)算的時(shí)間。由于多面體網(wǎng)格能以更少的網(wǎng)格單元數(shù)量劃分流場(chǎng),因此,在實(shí)際計(jì)算過(guò)程中多面體網(wǎng)格擁有比本文所給數(shù)據(jù)更大的效率優(yōu)勢(shì)。

(2)基于徑向基函數(shù)的網(wǎng)格變形方法,本文對(duì)多面體網(wǎng)格的大變形能力進(jìn)行了檢驗(yàn)。結(jié)果表明,變形后整體網(wǎng)格質(zhì)量的下降導(dǎo)致了收斂精度的降低,但是邊界層網(wǎng)格的質(zhì)量變化不大,可以保證計(jì)算結(jié)果的準(zhǔn)確性。

(3)基于多面體網(wǎng)格對(duì)AGARD算例進(jìn)行了驗(yàn)證。計(jì)算結(jié)果表明,多面體網(wǎng)格可以對(duì)機(jī)翼的靜氣動(dòng)彈性問(wèn)題進(jìn)行較為準(zhǔn)確的仿真。

由于目前還沒(méi)有合適的方法對(duì)多面體網(wǎng)格的變形質(zhì)量進(jìn)行評(píng)價(jià),因此本文只對(duì)網(wǎng)格變形前后的流場(chǎng)計(jì)算結(jié)果進(jìn)行了比較,并由此判定網(wǎng)格變形質(zhì)量滿(mǎn)足計(jì)算要求,并沒(méi)有對(duì)網(wǎng)格質(zhì)量的變化進(jìn)行精確的量化。如何對(duì)網(wǎng)格變形前后的質(zhì)量進(jìn)行量化和比較,從而在切割體、四面體和多面體網(wǎng)格中找到變形能力最強(qiáng)的網(wǎng)格以及各自適宜的變形方式是今后需要做的工作。

[1]王明強(qiáng),朱永梅,劉文欣.有限元網(wǎng)格劃分方法應(yīng)用研究[J].機(jī)械設(shè)計(jì)與制造,2004,2(1):22-24.

[2]李海峰,吳冀川,劉建波,等.有限元網(wǎng)格質(zhì)量剖分與網(wǎng)格質(zhì)量判定指標(biāo)[J].中國(guó)機(jī)械工程,2012,23(3):368-377.

[3]許曉平,周洲.多面體網(wǎng)格在CFD中的應(yīng)用[J].飛行力學(xué),2009,27(6):87-89.

[4]陳大偉,楊國(guó)偉.靜氣動(dòng)彈性計(jì)算方法研究[J].力學(xué)學(xué)報(bào),2009,41(4):469-479.

[5]Schmitt V,Charpin F. Pressure distributions on the ONERA-M6-Wing at transonic Mach numbers,experimental data base for computer program assessment[R].France:Report of the Fluid Dynamics Panel Working Group 04,AGARD AR 138,1979.

[6]張偉偉,高傳強(qiáng),葉正寅.氣動(dòng)彈性計(jì)算中網(wǎng)格變形方法研究進(jìn)展[J].航空學(xué)報(bào),2014,35(2):303-319.

[7]Hardy R L.Theory and applications of the multiquadric-biharmonic method[J].Computers & Mathematics with Applications,1990,19(8-9):163-208.

[8]Faul A C,Goodsell G,Powell M J D.A Krylov subspace algorithm for multiquadric interpolation in many dimensions[J].IMA Journal of Numerical Analysis,2005,25(1):1-24.

[9]Gumerov N A,Duraiswami R.Fast radial basis function interpolation via preconditioned Krylov iteration[J].SIAM Journal on Scientific Computing,2010,29(5):1876-1899.

[10]Beatson R K,Powell M J D,Tan A M.Fast evaluation of polyharmonic splines in three dimensions[J].IMA Journal of Numerical Analysis,2006,27(3):427-450.

[11]胡坤,李振北.ANSYS ICEM CFD工程實(shí)例詳解[M].北京:人民郵電出版社,2014:210-215.

[12]盧學(xué)成,葉正寅,張陳安.基于ANSYS/CFX耦合的機(jī)翼顫振分析[J].計(jì)算機(jī)仿真,2010,27(9):88-91.

[13]Yates E C Jr.AGARD standard aeroelastic configurations for dynamic response Ⅰ:wing 445.6 [R].NASA TM-100492,1987.

[14]Cai J,Liu F.Static aero-elastic computation with a coupled CFD and CSD method[R].AIAA-2000-0717,2000.

(編輯：姚妙慧)

Application of polyhedron grid in static aero-elastic computation

FENG Hao-yang, SU Xin-bing, MA Bin-lin, WANG Xu

(Aeronautics and Astronautics Engineering College, AFEU, Xi’an 710038, China)

The static aero-elastic computation with polyhedron grid was introduced. The CFD computational efficiency and deformation ability of polyhedron grid were studied based on M6 numerical example and interpolation method of radial basis function. The static aero-elastic numerical value was calculated with the use of AGARD 445.6 wing. The results show that the polyhedron grid possesses the advantages of computational efficiency, deformation ability and accurate calculation results, so it is suitable to the static aero-elastic computation.

polyhedron grid; static aero-elastic computation; grid efficiency; grid deformation

2015-10-21;

2016-03-14; 網(wǎng)絡(luò)出版時(shí)間:2016-04-22 09:52

馮浩洋(1991-)，男，河北高邑人，碩士研究生，研究方向?yàn)榍奥右盹w機(jī)的靜氣動(dòng)彈性現(xiàn)象及規(guī)律。

V211.3

A

1002-0853(2016)04-0024-05