大斜度四面體網格在高升力構型中的應用

肖中云，劉 剛，周 鑄，江 雄

（中國空氣動力研究與發展中心，四川 綿陽621000）

0 引 言

大型客機和軍用運輸機設計關鍵技術之一是起飛和著陸狀態下的機翼高升力構型氣動設計，要求在著陸狀態下獲得高升力，起飛狀態獲得高升力和低阻力。在飛機設計階段，計算流體動力學（CFD）分析起到關鍵的技術驗證作用。高升力構型下機翼前緣縫翼向下向前伸展，后緣富勒襟翼后退展開，增升部件打開留下的空隙和支撐機構完全暴露在氣流當中，出現流動分離、邊界層摻混等復雜流動現象，給CFD數值模擬帶來很大挑戰。

高升力構型的CFD數值模擬經過多年的發展，在多段翼型、機翼增升裝置和全機模擬上都得到了應用，文獻［1］指出現有CFD手段尚很難準確模擬高升力構型的最大升力系數及失速攻角，特別是對于存在明顯分離區的復雜流動。準確預測分離流動的開始和發展，以及雷諾數效應依然是CFD的難點之一。在網格方法上，結構網格的模擬歷史較長，文獻［2］采用了多塊對接網格模擬了梯形機翼高升力構型，文獻［3］采用結構重疊網格模擬了波音777－200飛機著陸構型，不足之處是在處理復雜結構的時候網格生成難度極大，網格生成時間占計算周期的70%以上。非結構網格具有對復雜外形適應性強的特點，可以實現網格的自動生成，文獻［4］采用非結構網格方法模擬了高升力構型復雜繞流，模擬包括了襟翼托架、凹槽、縫隙等細小結構，發揮了網格在處理復雜外形方面的優勢。當然非結構網格也存在缺點，如存儲量較大等，新的計算方法在不斷的發展和完善當中［5－7］。

本文工作在“航空CFD可信度開放式專題研究活動”支持下得到開展，計算外形為NASA高升力梯形機翼構型，網格劃分采用全四面體粘性非結構網格，為了模擬粘性效應，物面附近采用層推進方法生成大伸展比邊界層網格。流場計算采用了基于消息傳遞（MPI）的并行分區計算方法，通過數值計算分析了高升力構型的升阻力系數和力矩系數的變化規律，并通過試驗數據進行驗證，考察了計算的網格收斂性質，最后指出數值模擬中存在的不足和改進方向。

1 計算方法

流動控制方程為雷諾平均N－S方程，湍流模型采用Spalart－Allmaras一方程模型，湍流模型和主控方程之間為松耦合。在計算方式上采用了基于消息傳遞（MPI）的分區并行計算方法，控制方程的時間積分采用LU－SGS隱式方法，單元面通量計算采用Roe通量差分分裂格式，二階精度通過原始變量的線性重構獲得，并采用Barth－Jespersen限制器使格式保持單調性。

1.1 單元線性重構

有限體積離散采用基于格心的離散格式［8］，這樣守恒率方程表述為格心流場變量的時間變化與邊界面通量之間的關系。為了計算邊界通量，首先需要得到邊界面的流場變量。假設流場變量在各網格單元內是分片線性分布，則面心的流場值可以表示為格心值和其梯度的表達式。

為此，梯度計算的格式選擇尤為重要。對于格心離散格式來說，目前發展的梯度計算方法主要分為最小二乘法和基于節點平均方法兩類［9］。最小二乘法在控制單元周圍選取一定數量相鄰單元作為梯度計算的模板，算法采用最小二乘法。模板選擇可以是面相鄰單元、點相鄰單元和介于兩者之間的選取模式，對于四面體網格來說，面相鄰單元只有4個，點相鄰單元有幾十個之多，不同選取模式之間計算量差別很大。本文采用基于節點平均的梯度計算方法，模板單元如圖1所示，首先將格心流場變量平均到格點，再用格點平均得到面心，格心流場變量的梯度由在單元體上運用格林－高斯公式得到。

這樣在計算格心梯度的時候，其計算模板間接包括了與其共節點的所有相鄰單元。

圖1 格心計算格式的模板單元Fig.1 Stencils of cell centered scheme

公式（2）適用于任意形狀單元，對于單形網格單元（如四面體、棱柱、六面體等），邊界面流場變量可以表示成格心和格點值的簡單表達式［10］，其格心梯度在實質上等同于式（2）。式（1）中α代表限制器函數，目的是為了在單元重構中不出現新的極值點，滿足計算格式的保單調性要求。本文采用Barth－Jespersen［11］提出的非結構網格的斜率限制器，該方法限制邊界面上的重構值不超出該單元所有緊相鄰單元（包括單元本身）的最小值和最大值范圍。具體表達形式如下：

1.2 關于大伸展比拉伸網格

為了模擬粘性效應，網格生成通常沿物面或尾跡層生成大伸展比的拉伸網格，如圖2所示，為此網格單元在流向和法向上的尺度可能相差幾個量級。舉例來說，如果網格單元的高寬比h／L＝1∶100，此時對應的最小內角＝tg－1（1／100）＝0.57°，也就是說在粘性四面體網格中存在大量的這種尖劈形網格單元，需要在計算方法上慎重處理。

圖2 用于粘性計算大伸展比拉伸網格Fig.2 High aspect ratio meshes for viscous effect

式（2）中要用到面心的流場變量值φij，最簡單的方法是取相鄰單元格心值的平均，該方法的優勢在于構造簡單，不需要新增加數據結構，易于在復雜網格單元中實現。對高斜度的非結構網格來說，相鄰格心的連線不通過面心，所以這種簡單平均的方法存在誤差，并且斜度越大誤差越大，以圖2所示的四面體拉伸網格可以看到，面心a與格心連線ij相差很遠，所以在各向異性的非結構四面體網格中不能采用這種近似算法，而應代之以上節介紹的最小二乘法或基于節點平均的方法。

在應用限制函數的時候，最傳統的做法是用兩側單元的流場值對邊界面重構值進行限制，避免重構產生新的極值點。同樣因為格心連線與面心距離很遠的緣故，這種方法不能適用于粘性四面體網格，Barth－Jespersen斜率限制器提出用單元所有緊相鄰單元（包括單元本身和面相鄰單元）的最小值和最大值作為限制范圍，放寬了限制條件，避免限制過大帶來較大的數值耗散。

1.3 網格并行分區

計算外形為NASA高升力梯形機翼全展長襟翼模型［12］，該機翼沒有扭轉、沒有上反角，采用了大弦長和相對較小展弦比（AR＝4.56），目的是獲得較高的試驗雷諾數。該構型為典型的著陸構型，前緣縫翼與后緣襟翼的偏角分別為30°和25°。模型平均氣動弦長c＝1.0058m，半展長2.1615m，半模參考面積2.046m2，襟翼弦長和縫翼弦長分別為0.3m和0.127m。參考的風洞試驗是1998年在NASA Ames 12英尺增壓風洞（PWT）中完成的。

計算網格為全四面體非結構網格，網格在物面附近采用“層推進”方法生成的大伸展比拉伸網格，控制參數包括第一層網格物面距離、生長比、最大推進層數等，網格的法向尺度異于其他兩個方向，又稱為各向異性網格。第一層網格點的物面高度為2.5×10－5m，滿足格心的無量綱高度y＋≈1。粘性層網格的生長比1.25，計算域的遠場邊界取為100倍平均氣動弦長。

為了評價計算結果的網格收斂性，在同一外形上生成了粗網格、中等密度網格和密網格。表1給出了三套網格的對比，三套網格的第一層網格高度、生長比、遠場邊界等參數保持不變，主要差別體現在物面的網格密度不同，如圖3所示。圖4給出的是中等密度網格在機翼截面上的網格分布，從中可以看到沿物面法向生成的是大伸展比拉伸網格，目的是滿足粘性附面層計算的要求；此外，網格在縫翼和主翼之間分布合理、光滑過渡，基本上滿足高升力構型復雜空間流場對網格的要求。并行計算采用了Metis分區軟件［13］對網格進行分區，目標是保證各分區間負載平衡和通信量最小，以中等密度網格下32個分區為例，每個分區的網格單元數383016個，相差只有1個單元。信息交換量可以用每個分區的分區邊界面個數表示，其平均值為11696個，最大相差56%。圖5給出的是物面和對稱面上的網格分區情況，分區對接面傳遞信息包括對接單元和對應頂點的流場信息，在對接面上的流通量保持守恒。

表1 稀網格、中等網格和密網格比較Table 1 Comparison of coarse，medium and fine grids

圖3 三種網格在翼尖附近的比較Fig.3 Surface grid comparison at wing tip

圖4 機翼橫截面上的網格分布Fig.4 Grid distribution on wing’s cross section

圖5 并行計算的網格分區Fig.5 Grid partition for parallel computation

2 計算結果與分析

計算狀態選擇與試驗條件相同，馬赫數M＝0.15，攻角α＝－3.4°～36°，基于平均氣動弦長的雷諾數Re＝1.5×107，流動假設為全湍流。圖6給出的是采用中等密度網格計算得到的極曲線與試驗值的比較，從圖中可以看到，在所有計算點上計算值與試驗值基本吻合一致，不足是在高升力狀態下升阻比的計算值略低于試驗值，在低升力狀態下也存在一定偏差。在α＝11.02°時還分別用稀網格和密網格進行了計算，可以看到三種網格得到的升力系數和試驗值都比較一致。圖7給出的是俯仰力矩系數隨升力系數變化曲線與試驗值的比較，可以看到本文得到的俯仰力矩系數變化規律和試驗值比較一致，但在量值上存在差別，在升力系數（1.5＜CL＜2.5）范圍內預測的力矩系數絕對值比試驗值偏小。從單個點上稀網格和密網格的計算結果來看，網格加密使計算值更接近于試驗。

圖6 極曲線與試驗值的比較Fig.6 Polar curve compared with experimental data

圖7 俯仰力矩系數與試驗值比較Fig.7 Pitching moment compared with experimental data

圖8給出的是攻角α＝11.02°時的物面壓力云圖和流線分布，可以看到當前狀態下襟翼后緣開始有分離出現，并且在機翼翼尖存在強烈的翼尖渦流動，給翼尖區域的準確數值模擬帶來一定的難度。

圖8 物面壓力云圖與流線分布（α＝11°）Fig.8 Wall pressure contour and streamlines（α＝11°）

圖9給出了機翼不同展向截面上的壓力系數分布與試驗值的比較，其中計算網格采用中等密度網格，橫坐標表示測壓點位置（為和試驗一致單位用英寸表示），縫翼和襟翼位置為展開狀態下的坐標值，展向站位表示為機翼截面橫向坐標與半展長的比值。從圖中可以看到，在前面四個站位上壓力系數的計算值和試驗值都吻合較好，準確反映了各氣動部件的負壓峰值及壓力突變位置，并且在一些復雜流動區域（如前緣縫翼的內凹區域、主翼后緣存放襟翼的區域）計算值和試驗也基本一致。在第4個截面位置（η＝0.95）上，后緣襟翼的負壓峰值明顯低于試驗值，主翼的后半段背風區的壓力也出現偏差，其主要原因是在靠近翼尖位置存在有強烈的翼尖渦流動，在翼尖集中渦的影響下，上表面負壓增加，襟翼的負壓峰值顯著上升，本文在數值模擬中存在對集中渦數值耗散過大的缺陷，即數值模擬得到的渦強度偏小，相當于弱化了集中渦的影響，導致上翼面的壓力系數預測偏低。

圖9 機翼Y向截面上的壓力系數分布比較（α＝11.02°）Fig.9 Comparison of Cpdistribution on cross sections（α＝11.02°）

以攻角α＝11.02°為例，本文采用三種不同密度網格分別對流場進行了模擬，以探討數值計算的網格收斂性。其中稀網格、中等網格和密網格分別對應4.66M、12.26M和26.74M空間單元，圖10給出的是升力系數和阻力系數隨網格規模的變化曲線，橫坐標用網格指數（網格單元數的－2／3次方）表示，圖中直線表示試驗值。從圖上可以看到，在固定攻角下，升阻力系數的計算值比試驗值偏小，隨網格密度增加，差距縮小，計算值逐漸靠近試驗值。對于該圖形來說，越靠近橫軸左側網格密度越大，如果沿變化曲線畫一條趨勢線，則趨勢線與縱軸的交點代表網格無窮大時的計算值，稱此時對應的計算解為網格無關解。從圖10來看，升阻力系數并不完全隨網格指數線性變化，說明計算網格還沒有達到求網格無關解的要求。從圖9可以看到，在展向截面靠近翼尖位置（η＝0.95），計算的壓力系數較試驗值偏小，主要原因是翼尖位置的空間網格密度不夠，不能準確預測翼尖渦的影響。對于非結構網格來說，空間網格是自動填充而成的，其網格密度比較不易進行控制，這也是非結構網格面向工程應用需要解決的關鍵問題之一。

圖10 升力和阻力系數的網格收斂性Fig.10 Grid convergence of lift and drag coefficients

3 結 論

基于全四面體粘性非結構網格，采用格心格式求解雷諾平均N－S方程，對梯形機翼高升力構型流場進行了數值模擬，計算結果通過試驗數據進行了對比驗證，得到結論如下：

（1）高伸展比粘性四面體網格存在面心與相鄰單元格心不共線的特點，在面心流場變量重構當中，限制器函數不能僅用到相鄰單元的格心值進行限制，否則會帶來限制過強的問題。

（2）數值計算結果與試驗數據基本一致，驗證了本文方法的可行性，同時數值計算還存在翼尖渦數值耗散過大的問題，影響了翼尖區域壓力分布的預測精度。

（3）非結構網格面向工程應用需要解決好空間網格不易控制的問題，提高分離渦模擬的精度。

［1］朱自強，陳迎春，等.高升力系統外形的數值模擬計算［J］.航空學報，2005，26（3）：257－262.

［2］王運濤，洪俊武，等.采用TRIP軟件數值模擬梯形翼高升力構型［C］.中國航空學會2009年學術年會，北京：中國航空學會，2009.

［3］ROGERS S E，ROTH K.Advances in overset CFD processes applied to subsonic high－lift aircraft［R］.AIAA Paper 2000－4216，2000.

［4］CHAFFIN M S，PIRZADEH S.Unstructured Navier－Stokes high－lift computations on a trapezoidal wing［R］.AIAA Paper 2005－5084，2000.

［5］DISKIN B，THOMAS J L.Comparison of node－centered and cell－centered unstructured finite volume discretizations：inviscid fluxes［J］.AIAAJournal，2011，49（4）：836－854.

［6］BERGER M，AFTOSMIS M J.Analysis of slope limiters on irregular grids［R］.NAS Technical Report，NAS－05－007，2005.

［7］RUMSEY C，SUSAN X Y.Prediction of high lift：review of present CFD capability［J］.ProgressinAerospaceSciences，2002，38：145－180.

［8］BLAZEK J.Computational fluid dynamics：principles and applications［M］.Elsevier publisher，2001.

［9］MAVRIPLIS D J.Unstructured mesh discretizations and solvers for computational aerodynamics［R］.AIAA Paper 2007－3955，2007.

［10］FRINK N T.Recent progress toward a three－dimensional unstructured Navier－Stokes flow solver［R］.AIAA Paper 1994－0061，1994.

［11］BARTH T J，JESPERSEN D C.The design and application of upwind schemes on unstructured meshes［R］.AIAA Paper 1989－0366，1989.

［12］JOHNSON P L，JONES K M，MADSON M D.Experimental investigation of a simplified 3Dhigh lift configuration in support of CFD validation［R］.AIAA Paper 2000－4217，2000.

［13］KARYPIS G，KUMAR V METIS.A software pachage for partitioning unstructured graphs，partitioning meshes，and computing fill－reducing orderings of sparse matrices，version 4.0［M］.University of Minnesota，1998.