李 松，王光學，王運濤，張玉倫，鄧小剛

（1.空氣動力學國家重點實驗室，四川 綿陽 621000；2.中國空氣動力研究與發(fā)展中心，四川 綿陽 621000；3.國防科學技術大學，湖南 長沙 410073）

WCNS格式在梯形翼高升力構型模擬中的應用研究

李 松1，2，王光學1，2，王運濤1，2，張玉倫1，2，鄧小剛2，3

（1.空氣動力學國家重點實驗室，四川 綿陽 621000；2.中國空氣動力研究與發(fā)展中心，四川 綿陽 621000；3.國防科學技術大學，湖南 長沙 410073）

采用5階精度的有限差分格式 WCNS-E-5，對梯形翼高升力構型進行了數(shù)值模擬。研究了網(wǎng)格收斂性，WCNS格式對網(wǎng)格的依賴性更小。得到了與實驗結果相一致的氣動特性和壓力分布。通過與一些著名CFD軟件的計算結果比較表明，相比二階格式，WCNS格式模擬高升力構型有明顯優(yōu)勢，特別在失速迎角附近對流動結構的刻畫更準確。同時，驗證了程序對復雜外形的模擬能力和魯棒性。

高升力構型；高階精度；數(shù)值模擬；失速迎角；WCNS

0 引 言

由于增升裝置對飛機的起飛著陸性能有重要影響，增升系統(tǒng)及其空氣動力特性的研究一直是航空界的前沿課題。一般地，對于雙發(fā)的大型商用飛機，起飛時的升阻比每提高1%，有效載荷能增加2800磅；著陸時的最大升力系數(shù)每增加1.5%，有效載荷能增加6600磅［1］。隨著計算機軟、硬件技術的迅猛發(fā)展和計算空氣動力學研究的不斷深入，CFD技術已經(jīng)成為飛行器設計與評估的重要手段?；诶字Z平均Navier-Stokes方程的數(shù)值模擬軟件幾乎可以模擬所有復雜外形飛行器的繞流流場［2］，但巡航狀態(tài)數(shù)值模擬技術相對成熟，而高升力構型的數(shù)值模擬方法仍有待評估，面臨諸多困難。通常，高升力構型是由縫翼、主翼和襟翼組成的多段翼系統(tǒng)，其搭接和縫道等形成的復雜性已經(jīng)使得計算網(wǎng)格的生成變得很困難。而要模擬其包含有邊界層轉捩、流動分離與再附、尾跡與邊界層相互作用等［3］復雜流動現(xiàn)象的繞流（如圖1所示），對于CFD而言更是巨大的挑戰(zhàn)。

為了正確評估并推進現(xiàn)有CFD模擬高升力系統(tǒng)的技術水平，國際國內先后組織了許多專題活動，如歐洲的EUROLIFT 項目［4］，2010年6月AIAA 的第一屆高升力預測研討會（1st AIAA CFD High Lift Prediction Workshop，HiLiftPW-1）［5］，國內在2009年組織了首屆航空CFD可信度開放式專題研究活動（Aeronautics High Credible CFD Simulations Open Workshop，HiCFD）。梯形翼高升力構型是Hi Lift-PW-1和 HiCFD所選的標準模型。

高升力構型的數(shù)值模擬也是國內外學者關注的重點。Rumsey等［6］對近十五年來計算高升力系統(tǒng)的CFD方法進行了評估。朱自強等［7］對高升力系統(tǒng)數(shù)值模擬取得的成果、當前的水平以及存在的問題進行了闡述。Rumsey等［8］對 Hi Lift PW-1中提交的39組結果進行了比較和統(tǒng)計，評估了高升力構型數(shù)值模擬的當前水平，指出了需要進一步努力的方向。從已查閱的文獻來看［4-10］，當前高升力構型的數(shù)值模擬以二階精度為主，高階精度的模擬結果尚未看到。

圖1 三段翼型上可能呈現(xiàn)的流動現(xiàn)象Fig.1 Sketch of flow phenomena on a three-element airfoil

高升力構型繞流是典型的粘性主導的問題，這就要求數(shù)值方法能夠準確地模擬物理粘性。二階精度格式的數(shù)值粘性與物理粘性同是二階項，數(shù)值粘性很容易污染物理粘性，進而造成較大的模擬誤差。高階精度格式耗散小，更適合模擬高升力構型流動問題。然而，現(xiàn)有的高階精度格式多限于網(wǎng)格簡單的問題，應用于復雜的工程問題仍有較大困難。原因在于，有限體積方法在復雜網(wǎng)格上容易實現(xiàn)，但是高階精度格式的構造困難而且計算量太大。有限差分方法可以逐維操作，高階精度格式容易構造且計算量小，但不容易滿足守恒律，對網(wǎng)格質量要求較高，在復雜網(wǎng)格上實施比較困難?？梢哉f，高階精度格式應用于復雜外形的模擬是CFD研究中一項具有挑戰(zhàn)性的工作。

WCNS格式［11-15］是鄧小剛提出的 一 類 非 線 性高階精度有限差分格式，具有高分辨率、低耗散、低色散和一致高階精度等特點。自20世紀90年代以來，該團隊就一直致力于高階精度格式相關問題的研究，在格式、算法和應用等方面取得了突破性進展。文獻［16］針對高階有限差分算法在滿足幾何守恒律方面的困難，發(fā)展了滿足幾何守恒律的高階精度的坐標變換導數(shù)算法，大大提高了 WCNS格式對復雜網(wǎng)格的適應性和魯棒性。

本文對流項采用5階精度的 WCNS-E-5格式離散，粘性項采用6階精度的中心格式離散，對梯形翼高升力構型進行了數(shù)值模擬。研究了網(wǎng)格收斂性；得到了與實驗結果一致的氣動特性和壓力分布；研究了初始條件對失速迎角的影響。通過與一些著名CFD軟件的計算結果比較表明，相比二階格式，WCNS格式模擬翼梢渦的能力更好，對流動結構的刻畫更準確。同時，也驗證了程序對復雜外形的模擬能力和魯棒性。

1 計算方法

控制方程為雷諾平均N-S方程+SST兩方程湍流模型，對流項的離散采用WCNS-E-5格式，粘性項的離散采用6階精度中心格式，離散方程采用LUSGS求解。計算中采用了并行技術加速。

1.1 控制方程

曲線坐標系下，雷諾平均N-S方程為：

1.2 對流項離散格式

對流項離散采用原始變量型的WCNS-E-5格式，記網(wǎng)格間距為h，以ξ方向為例：

1）內點格式

2）邊界和靠近邊界格式

以格點0表示左邊界，以格點N 表示右邊界，則左邊界和靠近左邊界格式如下：

在N 點附近的右邊界格式與上式類似。

由于WCNS格式在文獻［11-14］有詳細描述，相關的插值公式及粘性項離散格式在這里不再贅述。坐標變換導數(shù)的計算使用了滿足幾何守恒律的高階精度的算法，文獻［16］中有詳細描述。

2 計算網(wǎng)格及狀態(tài)

2.1 計算網(wǎng)格

高升力構型是安裝在機身上的大弦長、半展、三段構型，機翼沒有扭轉、沒有上反角，采用大弦長和相對較小展弦比。風洞實驗設計了兩種襟翼構型，一種為全展長襟翼，單縫襟翼從翼梢一直延伸到翼根，并融于機身；另一種為半展長襟翼，襟翼的展長大約占模型展長的1／2，并置于機翼的中間位置。以上兩種襟翼構型具有相同的前緣縫翼，縫翼從翼梢一直延伸到翼根，并融于機身。本文中計算的是全展長襟翼構型。

CFD計算有兩組稍微不同的數(shù)?？梢允褂?，一組是梯形翼的原始CAD數(shù)據(jù)（as-designed），一組是2002年質量保證（QA）重新檢測的數(shù)據(jù)（as-built）［17］。本文使用的數(shù)模采用2002年重新檢測的數(shù)據(jù)。

按照HiLiftPW-1給出的網(wǎng)格指導［18］，本文使用Ansys公 司的 ICEM 軟 件生 成了 相同 拓撲 的 粗（Coarse）、中 （Medium）、細 （Fine）三 套對 接網(wǎng)格（point-to-point），均為369塊。網(wǎng)格使用H 型拓撲，邊界層網(wǎng)格采用O型拓撲。圖2給出了中等網(wǎng)格的網(wǎng)格拓撲、表面網(wǎng)格分布以及細節(jié)處理。三套網(wǎng)格的詳細統(tǒng)計信息如表1所示。

圖2 高升力構型中等網(wǎng)格Fig.2 Medium grids for high-lift configuration

表1 網(wǎng)格統(tǒng)計表Table 1 Grid statistics

2.2 計算狀態(tài)

梯形翼模型分別于98／99年和02／03年做了兩組實驗［17］，如圖3所 示。前一組實驗 稱為梯形翼實 驗（Trapezoidal Wing Experiment），涵蓋多個雷諾數(shù)下的一系列的構型。后一組實驗稱為3D高升力流動現(xiàn)象實驗（3-D High-Lift Flow Physics Experiment），更關注特定雷諾數(shù)（4.3×106）下的詳細的流動現(xiàn)象。

圖3 梯形翼高升力構型模型Fig.3 Fullspan model of high-lift trapwing configuration

兩次實驗為CFD研究高升力構型的流動現(xiàn)象提供了豐富的數(shù)據(jù)，在本文的計算中來流參數(shù)如下：

3 計算結果

按照HiLift-1的要求，本文中初始流場均使用自由來流，并采用全湍流進行模擬。本文計算的結果與T513次實驗以及多個著名CFD軟件計算的結果進行了比較（相關數(shù)據(jù)均來自 Hi Lift PW-1官方網(wǎng)站［18］），需要指出的是CFX使用轉捩模型進行模擬。

圖4給出了中等網(wǎng)格下典型迎角的收斂歷程曲線。四個計算狀態(tài)的殘差都迅速地下降了四個量級。升力系數(shù)基本穩(wěn)定，計算已經(jīng)收斂。這表明計算中幾何守恒得到了滿足。同時，也表明本文的計算方法對復雜網(wǎng)格具有良好的適應性。

圖4 高升力構型的收斂歷程曲線Fig.4 Convergence histories of high-lift configuration

3.1 網(wǎng)格收斂性

為了排除網(wǎng)格對計算結果的影響，對13°和28°兩個迎角進行了網(wǎng)格收斂性分析。圖5給出了網(wǎng)格收斂性的結果?？梢钥闯觯?3°迎角下本文計算的升力隨網(wǎng)格加密是單調收斂的，而28°迎角下本文計算的升力隨網(wǎng)格加密有一個波動，是振蕩收斂的。

圖5 高升力構型的網(wǎng)格收斂性曲線Fig.5 Grid-convergence of high-lift configuration

為了研究網(wǎng)格對計算結果的可信度的影響，表2給出了13°迎角下的Richardson外插值、基于外插值的相對誤差和網(wǎng)格收斂指標（Grid Convergence Index，GCI）。Richardson外插方法是Richardson于1910年首次使用［19］，基于已有網(wǎng)格的計算結果（要單調收斂）可以外推得到空間步長h趨于0時的結果。網(wǎng)格收斂指標由 Roache提出［20］，它是一種網(wǎng)格細化效果的統(tǒng)一度量方法，度量參數(shù)中引入了網(wǎng)格細化率r和收斂階p。相關計算公式可見文獻［21］。

從表2中可以看出，本文計算的結果是單調收斂的，實際收斂精度階高于其它計算結果，基于Richardson外插值的相對誤差遠遠小于其他計算結果。對于網(wǎng)格收斂指標 GCI，細網(wǎng)格指標GCI12均小于中網(wǎng)格指標GCI23，本文計算的結果小于其它計算結果的對應指標。這表明本文使用的網(wǎng)格對計算結果影響較小，計算結果可信度較高。

表2 13°迎角的Richardson外插和網(wǎng)格收斂指標Table 2 Richardson extrapolation and grid convergence index atα=13°

3.2 力特性比較

圖6是中等網(wǎng)格計算的氣動力系數(shù)。由圖中可見，本文的計算結果與實驗一致性很好。相對于其它典型軟件的計算結果，更接近實驗值，在非線性段特別是失速迎角附近，WCNS格式優(yōu)勢更明顯。

相對于實驗數(shù)據(jù)，本文計算的升力系數(shù)和阻力系數(shù)整體偏小，這與計算中沒有使用轉捩模型有關。根據(jù)文獻［22］的結果，使用轉捩模型后，邊界層厚度減小，升力系數(shù)和阻力系數(shù)都會增加。圖6中易見，使用轉捩模型的CFX的計算結果，在線性段更接近實驗值。這也表明了開展高升力構型模擬的轉捩研究的潛在必要性。

在飛行器／機翼的俯仰振蕩、動失速等過程中，其流場特征不僅與特定狀態(tài)有關還與到達該狀態(tài)的動態(tài)過程（如角速率等）有關，這被稱為遲滯現(xiàn)象（hysteresis）。靜失速過程中特定狀態(tài)的流場的特征會依賴于其歷史狀態(tài)［23］，迎角增加過程和迎角減小過程得到的同一迎角的（定常）流動的特征不相同，這也被稱為遲滯現(xiàn)象。高升力構型的實驗中在靜失速迎角附近出現(xiàn)了后一種遲滯現(xiàn)象。數(shù)值計算中發(fā)現(xiàn)，在失速迎角附近，不同的初場或算法可能導致不同的流場演化過程，進而導致不同流態(tài)的結果，文獻［22］也提到了這種現(xiàn)象。圖6中同時給出了 WCNS格式使用自由來流和小迎角解作為初場得到的結果。從圖中可見，自由來流初場計算的失速迎角為35°，小迎角解初場計算的失速迎角為37°，分別與T513實驗中迎角減小（r106）、迎角增大（r105）過程所測得的失速迎角是一致的。這表明本文的結果對失速迎角的預測更為準確。

圖6 高升力構型的氣動系數(shù)曲線Fig.6 Aerodynamics characteristic of high-lift configuration

3.3 表面壓力分布及空間流場模擬

圖7給出了13°、21°、28°和34°四個迎角表面壓力分布云圖以及表面流線。由圖中可見，各迎角下，主翼上弦向流動是最主要的，只在翼梢附近的小范圍內出現(xiàn)展向流動；隨著迎角的增加，主翼上的吸力峰顯著提高，襟翼上表面的分離區(qū)則逐漸減小以至完全消失，從而升力增大。由圖中還可以看到34°迎角時主機翼和襟翼上都能夠維持附著流動的狀態(tài)，沒有大范圍的分離。

圖7 高升力構型的表面壓力分布云圖及表面流線Fig.7 Surface pressure distribution and streamlines of high-lift configuration

圖8給出了34°迎角計算的翼梢渦的形態(tài)。圖9給出了28°和34°迎角三個典型站位上的壓力分布。從壓力分布對比來看，在機翼中段站位上，各計算結果的壓力分布與實驗值都符合得較好。在靠近翼梢的位置，二階格式的計算結果嚴重偏離實驗值，而高階的WCNS格式仍能給出與實驗值基本相符的預測，包括是接近失速的34°迎角。

圖8 34°迎角計算的翼梢渦Fig.8 Wing-tip vortex under 34°AOA

圖9 典型站位表面壓力系數(shù)的比較Fig.9 Cpdistribution for typical wing stations

文獻［8］指出，靠近主翼翼梢區(qū)域的流場很難預測，而準確、成功地捕捉翼梢渦是預測這些位置的表面壓力分布的非常重要的因素。對于翼梢附近的流動，二階精度格式難以準確捕捉翼梢附近流場的急劇變化。而高階精度格式的精度更高、數(shù)值粘性更小，可以顯著地減小二階精度格式的不足，對翼梢附近流動的模擬更準確，能夠準確地捕捉翼梢渦，對這些位置的表面壓力分布預測得較好。

4 結 論

本文采用高階精度WCNS-E-5格式，對梯形翼高升力構型進行了數(shù)值模擬，研究了網(wǎng)格收斂性，并對計算結果進行了較為細致的分析。相比二階精度格式，WCNS格式有以下特點：

（1）得到了網(wǎng)格收斂解的結果；對網(wǎng)格的依賴性更??；

（2）氣動特性與實驗值具有很好的一致性；模擬大迎角流動更有優(yōu)勢，對失速迎角的預測更準確；

（3）對翼梢附近的流動結構模擬得更準確，表面壓力分布與實驗值符合得更好；

（4）高階精度的WCNS格式更適合模擬高升力構型繞流問題；

（5）程序對復雜網(wǎng)格的適應能力比較強，收斂性比較好，比較魯棒。

［1］ MEREDITH P T.Viscous phenomena affecting high-lift systems and suggestions for future CFD development，high-lift systems aerodynamics［R］.AGARD CP 315，1993，191-198

［2］ TINOCO E N，BOGUE D R.Progress toward CFD for Full flight envelope［J］.Aeronautical Journal，2005，109（1100）：451-460.

［3］ SMITH A M O.High lift aerodynamics［J］.Journal of Aircraft，1975，12（6）：501-530.

［4］ HANSEN H ，THIEDE P，MOENS F，et al.Overview about the Europe high lift research program EUROLIFT［R］.AIAA 2004-0767.

［5］ SLOTNICK J P，HANNON J A，CHAFFIN M.Overview of the first AIAA CFD high lift prediction workshop （invited）［R］.AIAA 2011-862.

［6］ RUMSEY C L，YING S X.Prediction of high lift：review of present CFD capability［J］.Progress in Aerospace Science，2002，38（2）：145-180.

［7］ ZHU Z Q，CHEN Y C，WU Z C，et al.Numerical simulation of high lift system configration［J］.ACTA Aeronautica et Astronautica Sinica，2005，26（3）：257-262.（in Chinese）朱自強，陳迎春，吳宗成，等.高升力系統(tǒng)外形的數(shù)值模擬計算［J］.航空學報，2005，26（3）：257-262.

［8］ RUMSEY C L，LONG M，STUEVER R A，et al.Summary of the first AIAA CFD high lift prediction workshop （invited）［R］.AIAA 2011-939.

［9］ HONG J W，WANG Y T，PANG Y F，et al.Numerical research of high-lift configurations by structured mesh method ［J］.ACTA Aerodynamica Sinica，2013，31（1）：75-81.（in Chinese）洪俊武，王運濤，龐宇飛，等.結構網(wǎng)格方法對高升力構型的應用研究［J］.空氣動力學學報，2013，31（1）：75-81.

［10］WANG Y T，HONG J W，MENG D H.The influence of turbulent models to trapwing simulation［J］.ACTA Aerodynamica Sinica，2013，31（1）：52-55.（in Chinese）王運濤，洪俊武，孟德虹.湍流模型對梯形翼高升力構型的影響［J］.空氣動力學學報，2013，31（1）：52-55.

［11］DENG X G，MAEKAWA H，SHEN C.A class of high order dissipative compact schemes［R］.AIAA 96-1972.

［12］DENG X G，MAO M L.Weighted compact high-order nonlinear schemes for the Euler equations［R］.AIAA 97-1941.

［13］DENG X G，ZHANG H X.Developing high-order weighted compact nonlinear schemes［J］.Journal of Computational Physics，2000，165：24-44.

［14］DENG X G.High-order accurate dissipative weighted compact nonlinear schemes［J］.Science in China （Series A），2002，45 （3）：356-370.

［15］DENG X G，MAO M M，TU G H，et al.Extending the fifthorder weighted compact nonlinear scheme to complex grids with characteristic-based interface conditions［J］.AIAA Journal，2010，48（12）：2840-2851.

［16］DENG X G，MAO M M，et al.Geometric conservation law and applications to high-order finite difference schemes with stationary grids［J］.Journal of Computational Physics，2011，230：1100-1115.

［17］HANNON J A，WASHBURN A E，JENKINS L N，et al.Trapezoidal wing experimental repeatability and velocity profiles in the 14-by-22-foot subsonic tunnel（invited）［R］.AIAA 2012-0706.

［18］http：／／hiliftpw.larc.nasa.gov／index-workshop1.html

［19］RICHARDSON L F.The approximate arithmetical solution by finite differences of physical problems involving differential equations with an application to the stresses in a masonry dam ［J］.Transactions of the royal society of London（Series A），1908，210：307-357.

［20］ROACHE P J.Perspective：a method for uniform reporting of grid refinement studies［J］.ASME Journal of Fluids Engineering，1994（116）：405-413.

［21］OBERKAMPF W L，ROY C J.Verification and validation in scientific computing［M］.Cambridge University Press，2010.

［22］SCLAFANI A J，SLOTNICK J P，VASSBERG J C.Extended OVERFLOW analysis of the NASA trap wing wind tunnel model［R］.AIAA 2012-2919.

［23］YANG Z F，IGARASHI H，MARTIN M，et al.An experimental investigation on aerodynamic hysteresis of a low-Reynolds number airfoil［R］.AIAA 2008-0315.

Numerical simulation of high lift trapezoidal wing configuration with WCNS-E-5 scheme

LI Song1，2，WANG Guangxue1，2，WANG Yuntao1，2，ZHANG Yulun1，2，DENG Xiaogang2，3

（1.State Key Laboratory of Aerodynamics，Mianyang Sichuan 621000，China；2.China Aerodynamics Research and Development Center，Mianyang Sichuan 621000，China；3.National University of Defense Technology，Changsha Hunan 410073，China）

High order accuracy numerical simulations for flows of high lift trapezoidal wing configuration were performed by solving Navier-Stokes equation together with SST two-equation turbulence model.In the study，WCNS-E-5 scheme was utilized in the discretization of convection terms，while 6thorder central difference scheme for viscous term.Grid convergence result was obtained，which was in good agreement with test data，and WCNS scheme was less dependent on mesh density than 2ndorder scheme.Compared with other computational results obtained by typical CFD solvers，both the force coefficients and the surface pressure distributions given in this paper shown better accordance with the test data，and more reasonable prediction of flow structure near stalling was obtained.As revealed from these results，WCNS scheme took higher accuracy order and was more suitable for simulation of flow around high lift configurations than 2ndorder scheme.The capability and the robustness of the program to simulate complex high lift configurations were confirmed.

high lift configuration；high order scheme；numerical simulation；stalling angle；WCNS

V211.3

Adoi：10.7638／kqdlxxb-2012.0188

0258-1825（2014）04-0439-07

2012-11-12；

2013-09-05

國家重點基礎研究發(fā)展計劃（973計劃，2014CB744803.2）

李 松（1982-），男，助理工程師。研究方向：計算空氣動力學.E-mail：lisonic＠foxmail.com

李松，王光學，王運濤，等.WCNS格式在梯形翼高升力構型模擬中的應用研究［J］.空氣動力學學報，2014，32（4）：439-445.

10.7638／kqdlxxb-2012.0188. LI S，WANG G X，WANG Y T，et al.Numerical simulation of high lift trapezoidal wing configuration with WCNS-E-5 scheme［J］.ACTA Aerodynamica Sinica，2014，32（4）：439-445.