大斜度四面體網(wǎng)格在高升力構(gòu)型中的應(yīng)用

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

（中國空氣動(dòng)力研究與發(fā)展中心，四川 綿陽621000）

0 引 言

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

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

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

1 計(jì)算方法

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

1.1 單元線性重構(gòu)

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

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

這樣在計(jì)算格心梯度的時(shí)候，其計(jì)算模板間接包括了與其共節(jié)點(diǎn)的所有相鄰單元。

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

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

1.2 關(guān)于大伸展比拉伸網(wǎng)格

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

圖2 用于粘性計(jì)算大伸展比拉伸網(wǎng)格Fig.2 High aspect ratio meshes for viscous effect

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

在應(yīng)用限制函數(shù)的時(shí)候，最傳統(tǒng)的做法是用兩側(cè)單元的流場(chǎng)值對(duì)邊界面重構(gòu)值進(jìn)行限制，避免重構(gòu)產(chǎn)生新的極值點(diǎn)。同樣因?yàn)楦裥倪B線與面心距離很遠(yuǎn)的緣故，這種方法不能適用于粘性四面體網(wǎng)格，Barth－Jespersen斜率限制器提出用單元所有緊相鄰單元（包括單元本身和面相鄰單元）的最小值和最大值作為限制范圍，放寬了限制條件，避免限制過大帶來較大的數(shù)值耗散。

1.3 網(wǎng)格并行分區(qū)

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

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

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

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

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

圖4 機(jī)翼橫截面上的網(wǎng)格分布Fig.4 Grid distribution on wing’s cross section

圖5 并行計(jì)算的網(wǎng)格分區(qū)Fig.5 Grid partition for parallel computation

2 計(jì)算結(jié)果與分析

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

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

圖7 俯仰力矩系數(shù)與試驗(yàn)值比較Fig.7 Pitching moment compared with experimental data

圖8給出的是攻角α＝11.02°時(shí)的物面壓力云圖和流線分布，可以看到當(dāng)前狀態(tài)下襟翼后緣開始有分離出現(xiàn)，并且在機(jī)翼翼尖存在強(qiáng)烈的翼尖渦流動(dòng)，給翼尖區(qū)域的準(zhǔn)確數(shù)值模擬帶來一定的難度。

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

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

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

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

圖10 升力和阻力系數(shù)的網(wǎng)格收斂性Fig.10 Grid convergence of lift and drag coefficients

3 結(jié) 論

基于全四面體粘性非結(jié)構(gòu)網(wǎng)格，采用格心格式求解雷諾平均N－S方程，對(duì)梯形機(jī)翼高升力構(gòu)型流場(chǎng)進(jìn)行了數(shù)值模擬，計(jì)算結(jié)果通過試驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比驗(yàn)證，得到結(jié)論如下：

（1）高伸展比粘性四面體網(wǎng)格存在面心與相鄰單元格心不共線的特點(diǎn)，在面心流場(chǎng)變量重構(gòu)當(dāng)中，限制器函數(shù)不能僅用到相鄰單元的格心值進(jìn)行限制，否則會(huì)帶來限制過強(qiáng)的問題。

（2）數(shù)值計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)基本一致，驗(yàn)證了本文方法的可行性，同時(shí)數(shù)值計(jì)算還存在翼尖渦數(shù)值耗散過大的問題，影響了翼尖區(qū)域壓力分布的預(yù)測(cè)精度。

（3）非結(jié)構(gòu)網(wǎng)格面向工程應(yīng)用需要解決好空間網(wǎng)格不易控制的問題，提高分離渦模擬的精度。

［1］朱自強(qiáng)，陳迎春，等.高升力系統(tǒng)外形的數(shù)值模擬計(jì)算［J］.航空學(xué)報(bào)，2005，26（3）：257－262.

［2］王運(yùn)濤，洪俊武，等.采用TRIP軟件數(shù)值模擬梯形翼高升力構(gòu)型［C］.中國航空學(xué)會(huì)2009年學(xué)術(shù)年會(huì)，北京：中國航空學(xué)會(huì)，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.