多面體網格在靜氣動彈性計算中的應用

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

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

?

多面體網格在靜氣動彈性計算中的應用

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

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

介紹了多面體網格技術在靜氣動彈性問題中的應用。基于M6算例和徑向基函數插值法,探究了多面體網格的CFD計算效率和變形能力,使用AGARD445.6機翼進行了靜氣動彈性數值計算。結果表明,多面體網格具有計算效率高、變形能力強、計算結果準確等優點,可適用于靜氣動彈性問題的數值計算。

多面體網格; 靜氣動彈性計算; 網格效率; 網格變形

0　引言

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

本文對四面體、切割體、多面體網格的計算效率進行了比較,檢驗了多面體網格的變形能力。在文獻[4]方法的基礎上,基于多面體網格和徑向基函數的網格變形方法,通過AGARD算例,對機翼的靜氣動彈性現象進行了仿真,驗證了多面體網格的有效性。

1　網格計算效率

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

圖1　翼根前緣網格示意圖Fig.1　The mesh of wing root leading edge

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

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

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

三種網格下的計算殘差曲線如圖3所示。由圖3可以看到,三者的收斂速度大致相同,多面體網格和四面體網格在150步左右收斂,切割體網格在180步左右收斂。四面體網格收斂速度快是因為其節點數最少,易于收斂,但其收斂精度最低。多面體網格在與四面體網格收斂速度相當的同時,收斂精度遠比其高,與切割體網格大致相同,但多面體網格收斂過程更加平穩。

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

三種網格的收斂精度和計算時間如表1所示。由于切割體網格在計算過程中耗散最小,且M6機翼外形簡單,因此其所用時間最短,而多面體網格計算時間小于四面體網格。

表1　收斂精度和計算時間Table1　Convergenceaccuracyandcomputationtime

由計算結果可知,多面體網格在流體計算過程中綜合了四面體網格和切割體網格的優點。當網格數相同時,多面體網格在收斂速度快的同時還有較高的精度和平穩度。因此,相較于切割體網格和四面體網格,采用多面體網格進行靜氣動彈性計算可以提高CFD的計算效率和準確性,從而縮短整體計算時間。

2　網格變形能力

為了檢驗多面體網格的變形能力[6],本文基于徑向基函數變形方法[7-10],對M6機翼的網格進行了旋轉變形,并對變形前后機翼的氣動參數進行了計算和對比。

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

圖4　機翼旋轉示意圖Fig.4　The wing rotation diagram

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

由圖可見,在機翼旋轉變形45°的情況下,其前后緣網格都發生了不同程度的偏斜。后緣偏斜程度較大,這是因為后緣尺寸較小,其邊界層網格的外邊界與后緣壁面的距離相對較遠所致。基于徑向基函數的變形方法對遠離壁面網格節點的控制力較弱,在大變形情況下,網格經過多次迭代變形會產生較大的誤差,導致后緣處網格偏斜更為嚴重。

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

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

依照上節邊界條件對變形后的流場網格進行計算,得到如圖7所示的變形前后機翼沿不同展向的壓力分布。由圖7可知,盡管前緣和后緣處的網格發生了偏斜,但是計算結果仍然是準確的,可以認為局部網格的偏斜對計算沒有造成負面影響。

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

網格變形前后的收斂曲線如圖8所示。由圖8可知,曲線都在260步左右達到收斂;但是變形后網格計算結果的收斂曲線不如變形前平穩,而且收斂精度比變形前略差。這是因為變形后的網格整體質量有所下降,致使計算過程中出現偽耗散,使得收斂精度降低。盡管如此,由于計算結果的準確性并沒有受到影響,所以仍然說明多面體網格在該變形方法產生大變形后可以滿足計算要求。

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

3　流體、結構計算和數據傳遞方法

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

4　算例驗證

采用標準的AGARD445.6機翼對本文的氣動彈性計算方法進行驗證。AGARD機翼為NACA65A004翼型,展弦比為1.65,梢根比為0.66,展長為0.76 m,根弦長為0.56 m,弦線后掠角為45°[13]。剖分的流體網格和結構網格如圖9所示,流體網格數為2 401 689,結構網格數為9 593。

圖9　AGARD機翼流體網格和結構網格Fig.9　The fluid mesh and structure mesh of AGARD wing

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

計算狀態選擇Ma=0.8,α=1°。AGARD 445.6機翼靜氣動彈性迭代平衡后,相對于初始位置,翼梢前緣和后緣的變形量如表2所示。由表2可以看出,本文方法與文獻[14] 的計算結果基本一致。變形量有誤差是由于插值方法的限制和網格數量的不同,數據在傳遞過程中會有一定的偏差所致。由于結果相差都在3%以內,所以這樣的誤差是可以接受的。

表2　AGARD 445.6機翼靜氣動彈性結果比較Table2　ComparisonofAGARD445.6wingstaticaero-elasticcomputationalresults

5　結束語

基于多面體網格,本文對靜氣動彈性計算的效率和網格變形問題進行了如下三方面的探究:

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

(2)基于徑向基函數的網格變形方法,本文對多面體網格的大變形能力進行了檢驗。結果表明,變形后整體網格質量的下降導致了收斂精度的降低,但是邊界層網格的質量變化不大,可以保證計算結果的準確性。

(3)基于多面體網格對AGARD算例進行了驗證。計算結果表明,多面體網格可以對機翼的靜氣動彈性問題進行較為準確的仿真。

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

[1]王明強,朱永梅,劉文欣.有限元網格劃分方法應用研究[J].機械設計與制造,2004,2(1):22-24.

[2]李海峰,吳冀川,劉建波,等.有限元網格質量剖分與網格質量判定指標[J].中國機械工程,2012,23(3):368-377.

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

[4]陳大偉,楊國偉.靜氣動彈性計算方法研究[J].力學學報,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]張偉偉,高傳強,葉正寅.氣動彈性計算中網格變形方法研究進展[J].航空學報,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工程實例詳解[M].北京:人民郵電出版社,2014:210-215.

[12]盧學成,葉正寅,張陳安.基于ANSYS/CFX耦合的機翼顫振分析[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; 網絡出版時間:2016-04-22 09:52

馮浩洋(1991-)，男，河北高邑人，碩士研究生，研究方向為前掠翼飛機的靜氣動彈性現象及規律。

V211.3

A

1002-0853(2016)04-0024-05