章錦威, 張小亮, 孫曉玲, 戚淑妮

(中國(guó)航空工業(yè)空氣動(dòng)力研究院, 沈陽 110034)

0 引 言

阻力高精度計(jì)算是民用飛機(jī)設(shè)計(jì)中計(jì)算流體力學(xué)仿真的主要研究方面，也是飛行器外形氣動(dòng)特性分析中非常關(guān)注的問題。為了系統(tǒng)評(píng)估CFD計(jì)算阻力的發(fā)展水平，國(guó)際上先后組織了許多專題研討會(huì)，比如美國(guó)AIAA應(yīng)用空氣動(dòng)力學(xué)會(huì)委員1998年成立了CFD阻力預(yù)測(cè)工作小組DPW(Drag Prediction Workshop)，已于2001年開始召開了六次阻力預(yù)測(cè)會(huì)議[1-6]。NASA還專門發(fā)布了NPARC Alliance CFD驗(yàn)證與確認(rèn)檔案官方網(wǎng)站[7]。近年來，國(guó)內(nèi)的科研單位[8-9]也陸續(xù)地開展了有關(guān)國(guó)際標(biāo)模驗(yàn)證和確認(rèn)的研究工作，初步設(shè)計(jì)了驗(yàn)證和確認(rèn)共享數(shù)據(jù)庫，并利用這些數(shù)據(jù)資源進(jìn)行了航空CFD的驗(yàn)證與確認(rèn)工作。鄧小剛[10]等對(duì)國(guó)內(nèi)外CFD驗(yàn)證與確認(rèn)的發(fā)展情況，并對(duì)CFD軟件的可信度分析工作提出了若干建議。陳樹生[11]等針對(duì)大型CFD軟件的驗(yàn)證和確認(rèn)，專門開發(fā)了具有標(biāo)準(zhǔn)共享的專業(yè)數(shù)據(jù)庫，適用于CFD軟件開發(fā)、測(cè)試、驗(yàn)證和評(píng)估等領(lǐng)域的自動(dòng)化驗(yàn)證確認(rèn)云平臺(tái)。

總體而言，國(guó)內(nèi)相關(guān)活動(dòng)缺乏完整性，且試驗(yàn)?zāi)Ｐ团c試驗(yàn)數(shù)據(jù)基本采用國(guó)外數(shù)據(jù)，CFD分析不夠深入。運(yùn)輸機(jī)標(biāo)模CHN-T1是由中國(guó)空氣動(dòng)力研究與發(fā)展中心自主設(shè)計(jì)的一個(gè)跨聲速運(yùn)輸機(jī)的典型算例，在FL-13風(fēng)洞、DNW-LLF風(fēng)洞進(jìn)行了低速測(cè)力試驗(yàn)，在FL-26風(fēng)洞、ETW風(fēng)洞進(jìn)行了高速測(cè)力試驗(yàn)，具有比較豐富可靠的試驗(yàn)數(shù)據(jù)[12]。

本文采用自研的非結(jié)構(gòu)混合網(wǎng)格解算器CARIA-OVERSET對(duì)運(yùn)輸機(jī)標(biāo)模CHN-T1進(jìn)行了驗(yàn)證確認(rèn)計(jì)算研究。根據(jù)第一屆航空CFD可信度研討會(huì)-CHN-T1標(biāo)模CFD可信度研討會(huì)活動(dòng)要求，對(duì)運(yùn)輸機(jī)標(biāo)模CHN-T1開展了網(wǎng)格收斂性研究和計(jì)算精度研究。主要研究了網(wǎng)格形式、網(wǎng)格尺度、模擬條件等對(duì)運(yùn)輸機(jī)標(biāo)模CHN-T1的總體氣動(dòng)特性、局部流場(chǎng)特征、壓力分布以及計(jì)算收斂性等的影響。

1 數(shù)值方法

本文的流場(chǎng)計(jì)算主要采用中國(guó)航空工業(yè)空氣 動(dòng)力研究院自主研發(fā)的航空數(shù)值計(jì)算平臺(tái)CARIA_OVERSET，該平臺(tái)采用非結(jié)構(gòu)網(wǎng)格技術(shù)和基于對(duì)偶網(wǎng)格的有限體積方法，以EULER方程或NS方程加湍流模型作為主控方程，用于求解定常和非定常流動(dòng)的流場(chǎng)。無黏通量差分格式包括通量矢量分裂和通量差分分裂格式，其離散格式主要有完全迎風(fēng)格式、Frommer格式、中心格式和二階迎風(fēng)偏置格式。在時(shí)間推進(jìn)方面，采用隱式近似因子分解迭代方法，其中定常流動(dòng)采用一階格式，非定常流動(dòng)則采用隱式雙時(shí)間步長(zhǎng)方法，并通過隱式殘值光順和多重網(wǎng)格等技術(shù)加速收斂，縮短計(jì)算時(shí)間。

此外，CARIA_OVERSET包含豐富的湍流模型：零方程的Baldwin-Lomax模型，一方程的Baldwin-Barth、Spalart-Allmaras模型，兩方程的Wilcoxk-ω、SST、k-ε等十幾種湍流模型，以及豐富的可用于各種航空專業(yè)領(lǐng)域的計(jì)算邊界條件。另外，CARIA_OVERSET中包含了eN方法、γ-θ轉(zhuǎn)捩模型方法和固定轉(zhuǎn)捩預(yù)測(cè)方法等，能夠?qū)崿F(xiàn)對(duì)不同轉(zhuǎn)捩機(jī)理作用下的轉(zhuǎn)捩問題的準(zhǔn)確預(yù)測(cè)。

在本文計(jì)算中，對(duì)流項(xiàng)采用為二階迎風(fēng)Roe格式進(jìn)行離散，二階重構(gòu)采用的梯度方法為節(jié)點(diǎn)型Gauss方法。黏性項(xiàng)采用中心差分離散。時(shí)間推進(jìn)方法為隱式的LU-SGS方法，采用為二重“V”循環(huán)多重網(wǎng)格加速收斂，全局CFL數(shù)設(shè)置為1。計(jì)算假定流場(chǎng)為完全湍流流場(chǎng)，湍流模型主要采用一方程Spalart-Allmaras模型、兩方程的Wilcoxk-ω和SST模型。計(jì)算主要采用百核分區(qū)并行計(jì)算，迭代步數(shù)10 000步～20 000步之間。

γ-θ模型求解兩個(gè)變量的標(biāo)準(zhǔn)輸運(yùn)方程：

2 計(jì)算外形和網(wǎng)格

運(yùn)輸機(jī)CHN-T1標(biāo)模采用類似波音737、空客320、C919等飛機(jī)的布局形式，即窄體機(jī)身、下單翼形式的超臨界機(jī)翼和平尾、單立尾、翼吊式發(fā)動(dòng)機(jī)布局。CHN-T1標(biāo)模選取面積95.346 m2、展弦比9.3、梢根比0.298、1/4弦線后掠角25°的機(jī)翼設(shè)計(jì)參數(shù)、巡航馬赫數(shù)為0.78，巡航升力系數(shù)為0.5。該模型幾何尺寸合理，有利于開展風(fēng)洞試驗(yàn)?zāi)Ｐ驮O(shè)計(jì)加工和CFD計(jì)算的網(wǎng)格制作，可為各種CFD軟件代碼提供驗(yàn)證計(jì)算的標(biāo)準(zhǔn)算例, 其詳細(xì)參數(shù)可見文獻(xiàn)[15]。本文主要針對(duì)運(yùn)輸機(jī)標(biāo)模CHN-T1進(jìn)行了驗(yàn)證確認(rèn)計(jì)算研究，該模型的幾何特征參數(shù)見表1。

本文的網(wǎng)格收斂性研究，主要采用中國(guó)空氣動(dòng)力研究與發(fā)展中心提供的官網(wǎng)網(wǎng)格和自研網(wǎng)格。官網(wǎng)網(wǎng)格包括粗中細(xì)不同尺度的非結(jié)構(gòu)混合網(wǎng)格[16]和結(jié)構(gòu)多塊網(wǎng)格[17]，其中粗、細(xì)結(jié)構(gòu)多塊網(wǎng)格是以中等網(wǎng)格為基準(zhǔn)自主生成，網(wǎng)格生成參考組委會(huì)提供的網(wǎng)格生成標(biāo)準(zhǔn)和規(guī)范。網(wǎng)格具體信息參見表2，不同類型、中等尺度表面網(wǎng)格參見圖1。

表1 幾何特征參數(shù)Table 1 Geometric characteristic parameters

表2 網(wǎng)格信息說明Table 2 Grid information description

圖1 不同類型、中等尺度計(jì)算網(wǎng)格示意圖Fig.1 Schematic diagram of different grid types

3 計(jì)算狀態(tài)介紹

計(jì)算狀態(tài)主要包括三組研究案例，分別為網(wǎng)格收斂性研究、支撐機(jī)構(gòu)及靜氣彈變形影響研究和雷諾數(shù)及轉(zhuǎn)捩影響研究，具體計(jì)算狀態(tài)參見表3。

表3 計(jì)算狀態(tài)說明Table 3 Calculation status specification

4 數(shù)值仿真結(jié)果分析

4.1 計(jì)算收斂歷程曲線

圖2是解的收斂歷史曲線。本文的氣動(dòng)力收斂原則，保證最后1000步升力系數(shù)波動(dòng)幅值小于0.001，阻力系數(shù)波動(dòng)幅值小于0.0001，俯仰力矩系數(shù)波動(dòng)幅值小于0.001。由圖可以看出，在保證氣動(dòng)力收斂準(zhǔn)則的前提下，所有狀態(tài)的密度殘差都下降了5個(gè)以上量級(jí)，網(wǎng)格尺度越小，收斂速度越慢。

(a) 非結(jié)構(gòu)混合網(wǎng)格

(b) 結(jié)構(gòu)網(wǎng)格

4.2 網(wǎng)格尺度影響分析

為了進(jìn)行網(wǎng)格收斂性的分析，至少要進(jìn)行三個(gè)不同網(wǎng)格尺度的計(jì)算，獲得S1(細(xì))、S2(中)和S3(粗)三個(gè)解，它們之間的差及其比率R為:

R=ε21/ε32(3)

網(wǎng)格收斂性條件為：0

為了展示網(wǎng)格密度對(duì)誤差估計(jì)計(jì)算結(jié)果的影響，本文對(duì)不同網(wǎng)格形式的粗中細(xì)各三套網(wǎng)格進(jìn)行收斂性分析，分析的數(shù)據(jù)分別是阻力系數(shù)CD、壓差阻力CD_PR、摩擦阻力CD_SF、俯仰力矩系數(shù)Cm和修正值。其中，網(wǎng)格誤差δ*(廣義Richardson外推法，忽略高階項(xiàng)的影響) 和修正值SC相關(guān)的計(jì)算公式為:

SC=S1-δ*(6)

式(4)、(5)中，r為不同網(wǎng)格尺度間的比值，本文取值為1.44。

圖3給出的是不同形式網(wǎng)格計(jì)算所得氣動(dòng)力系數(shù)隨網(wǎng)格數(shù)量變化對(duì)比曲線，其中，N為網(wǎng)格節(jié)點(diǎn)數(shù)。由圖可以看出，非結(jié)構(gòu)混合網(wǎng)格和結(jié)構(gòu)網(wǎng)格均有一定的網(wǎng)格收斂趨勢(shì)，且非結(jié)構(gòu)混合網(wǎng)格的收斂趨勢(shì)更明顯。

(a) 迎角和俯仰力矩

(b) 總阻力、壓差阻力及摩擦阻力

表4是在Ma=0.78、CL=0.500±0.001時(shí)，不同形式、不同尺度網(wǎng)格計(jì)算所得的氣動(dòng)力、力矩以及部分試驗(yàn)數(shù)據(jù)?？梢钥闯?，在固定升力條件下，計(jì)算迎角普遍小于試驗(yàn)迎角；網(wǎng)格尺度對(duì)摩擦阻力和壓差阻力均有一定影響，摩擦阻力差量為2 counts(非結(jié)構(gòu)混合網(wǎng)格)和0.7 counts(結(jié)構(gòu)網(wǎng)格)，壓差阻力差量為31 counts(非結(jié)構(gòu)混合網(wǎng)格)和6 counts(結(jié)構(gòu)網(wǎng)格)；在固定升力條件下，不同的網(wǎng)格尺度對(duì)力矩系數(shù)影響較小，計(jì)算得到的抬頭力矩普遍大于試驗(yàn)結(jié)果。

依據(jù)網(wǎng)格收斂性條件要求，非結(jié)構(gòu)混合網(wǎng)格計(jì)算所得氣動(dòng)力系數(shù)均為單調(diào)收斂，可以采用Richardson外推法初步估算，而結(jié)構(gòu)網(wǎng)格計(jì)算所得俯仰力矩系數(shù)為震蕩收斂，只能得到其誤差范圍，而其它氣動(dòng)力系數(shù)均發(fā)散，誤差和不確定度均不能進(jìn)行估算。非結(jié)構(gòu)混合網(wǎng)格采用Richardson外推法初步估算所得各項(xiàng)氣動(dòng)力參數(shù)與結(jié)構(gòu)網(wǎng)格線性插值計(jì)算結(jié)果接近，與文獻(xiàn)[17]提供的104億網(wǎng)格相比，除迎角存在一定差異外，其余結(jié)果接近。

表4 網(wǎng)格收斂性分析Table 4 Grid convergence analysis

圖4～圖5給出的是機(jī)翼和平尾不同展向站位及對(duì)應(yīng)壓力分布情況對(duì)比。由圖可以看出，非結(jié)構(gòu)混合網(wǎng)格計(jì)算所得機(jī)翼內(nèi)側(cè)站位的Cp分布差異較小，差異由內(nèi)向外逐漸增大，網(wǎng)格尺度越大，激波位置和強(qiáng)度的捕捉越精確，平尾內(nèi)側(cè)站位的吸力峰值差異明顯。結(jié)構(gòu)多塊網(wǎng)格計(jì)算所得機(jī)翼和平尾不同展向站位的Cp分布差異較小。以上結(jié)果與前文所得氣動(dòng)力系數(shù)分析結(jié)論一致。

圖4 機(jī)翼及平尾展向站位示意圖Fig.4 Schematic diagram of wing and tail deployment position

圖6給出的是固定升力條件下，非結(jié)構(gòu)混合網(wǎng)格采用不同網(wǎng)格尺度計(jì)算所得翼身結(jié)合處分離流場(chǎng)結(jié)果。由圖可以看出，在上述來流條件下，翼身結(jié)合處產(chǎn)生了很小的分離區(qū)，除粗網(wǎng)格模擬結(jié)果存在差異外，中細(xì)不同網(wǎng)格尺度在翼身結(jié)合處局部分離流的模擬上結(jié)果接近。

(a) 機(jī)翼

(b) 平尾

圖6 翼身結(jié)合處上表面流線(不同網(wǎng)格尺度)Fig.6 Upper surface streamline near wing body junction(different grid scales)

4.3 湍流模型影響分析

表5是在Ma=0.78、CL=0.500±0.001時(shí)，采用一方程的Spalart-Allmaras模型，兩方程的Wilcoxk-ω、SST共計(jì)三種湍流網(wǎng)格計(jì)算所得的氣動(dòng)力、力矩以及部分試驗(yàn)數(shù)據(jù)?？梢钥吹剑谕惶追墙Y(jié)構(gòu)混合網(wǎng)格，采用不同湍流模型計(jì)算所得結(jié)果存在一定差異，SA和SST兩種湍流模型的計(jì)算結(jié)果較為接近，SST計(jì)算所得阻力和迎角與試驗(yàn)結(jié)果更加接近。

表5 湍流模型對(duì)氣動(dòng)力的影響分析Table 5 Influence of turbulence models on aerodynamic forces

圖7給出的是固定升力條件下，非結(jié)構(gòu)混合網(wǎng)格不同湍流模型計(jì)算所得翼身結(jié)合處分離流場(chǎng)結(jié)果。由圖可以看出，在上述來流條件下，翼身結(jié)合處均產(chǎn)生了較小的分離區(qū)，不同湍流模型計(jì)算所得機(jī)翼弦向、展向分離邊界位置以及機(jī)身的分離泡中心、鞍點(diǎn)位置等均非常接近。

4.4 支撐機(jī)構(gòu)及靜氣彈變形影響研究

圖8(a～c)分別給出了干凈模型、帶支撐機(jī)構(gòu)模型和靜氣彈變形機(jī)翼帶支撐機(jī)構(gòu)模型計(jì)算所得全機(jī)的升力、俯仰力矩以及總阻力系數(shù)對(duì)比曲線。由圖可以看出，本文方法計(jì)算所得全機(jī)升阻力系數(shù)與試驗(yàn)數(shù)據(jù)吻合度較好，俯仰力矩曲線規(guī)律與試驗(yàn)數(shù)據(jù)較為相似。在本文計(jì)算條件下，不同構(gòu)型計(jì)算所得全機(jī)升阻力系數(shù)及俯仰力矩系數(shù)有不同程度變化：

(1) 帶支撐機(jī)構(gòu)模型相比干凈模型，升力系數(shù)最大增加0.0172，阻力系數(shù)最多減小9 counts(壓差阻力-6.5 counts和摩擦阻力-2.4 counts)，最大升阻比增加1，低頭力矩系數(shù)最大增加0.0487；

(2) 靜氣彈變形機(jī)翼帶支撐機(jī)構(gòu)模型相比干凈模型，升力系數(shù)最多減小0.0092，阻力系數(shù)最多減小22.5 counts(壓差阻力-23.5 counts和摩擦阻力+1 counts)，最大升阻比增加0.1， 低頭力矩系數(shù)最多減小0.008 53。

圖9～圖10給出了當(dāng)AOA=3.75°時(shí)，干凈模型、帶支撐機(jī)構(gòu)模型和靜氣彈變形機(jī)翼帶支撐機(jī)構(gòu)模型計(jì)算所得機(jī)翼和平尾不同展向站位Cp分布曲線。由圖可以看出，支撐機(jī)構(gòu)主要影響平尾(與支撐方式有關(guān))，支撐機(jī)構(gòu)導(dǎo)致截面Cp環(huán)量減小，且越靠?jī)?nèi)側(cè)影響越大；靜氣彈變形主要影響機(jī)翼，靜氣彈變形機(jī)翼上表面吸力減小，激波位置前移，且越靠外側(cè)影響越大。

圖8 不同構(gòu)型計(jì)算所得全機(jī)氣動(dòng)力系數(shù)曲線Fig.8 Aerodynamic coefficient curves calculated by different configurations

圖9 AOA=3.75°，機(jī)翼不同展向站位Cp分布變化對(duì)比Fig.9 Comparison of Cp distribution at different wing deployment stations and AOA=3.75°

圖11～圖12給出了不同迎角條件下，機(jī)翼上表面流線分布及展向中間站位Cp分布曲線。當(dāng)迎角AOA從3°增加到3.75°時(shí)，激波位置后移，翼中激波位置附近開始出現(xiàn)分離；當(dāng)迎角AOA從3.75°增加到4.25°時(shí)，激波位置前移，機(jī)翼后緣分離位置迅速前

移。激波位置的變化引起局部流動(dòng)分離特征的變化，這是誘發(fā)跨聲速機(jī)翼抖振問題的主要因素。

圖11 不同迎角下，不同構(gòu)型計(jì)算所得機(jī)翼表面流線對(duì)比Fig.11 Comparison of wing surface streamlines calculated from different configurations at different angles of attack

圖12 不同迎角下，η=0.50機(jī)翼展向站位的Cp分布對(duì)比Fig.12 Comparison of Cp distribution at η=0.50 wing deployment station and different angles of attack

4.5 雷諾數(shù)及轉(zhuǎn)捩影響研究

表6給出了不同模擬條件計(jì)算所得各項(xiàng)氣動(dòng)參數(shù)及其差量值?？梢钥闯?，雷諾數(shù)從330萬增加到1500萬，全湍模擬計(jì)算所得全機(jī)總阻力減小48 counts(壓差阻力減小20 counts、摩擦阻力減小28 counts)迎角減小0.249°，抬頭力矩增加0.0149；自由轉(zhuǎn)捩模型計(jì)算所得全機(jī)總阻力減小16 counts(壓差阻力減小7 counts、摩擦阻力減小9 counts)，迎角減小0.053°，抬頭力矩增加0.0176。全湍計(jì)算結(jié)果與試驗(yàn)結(jié)果更接近，自由轉(zhuǎn)捩模型計(jì)算所得差量偏小。

圖13～圖15給出了全湍及自由轉(zhuǎn)捩模型計(jì)算所得表面Cf分布云圖、機(jī)翼不同展向站位Cp及Cf曲線?？梢钥闯?，雷諾數(shù)越大，轉(zhuǎn)捩位置越靠前，全湍與自由轉(zhuǎn)捩模型計(jì)算結(jié)果的差異越小，這與上文所得氣動(dòng)力變化規(guī)律一致。在固定升力條件下，機(jī)翼不同展向站位Cp曲線差異小。

表6 CL=0.500，不同模擬條件計(jì)算所得氣動(dòng)力差量對(duì)比Table 6 CL=0.500, aerodynamic coefficients calculated from different simulation conditions

圖13 不同模擬條件計(jì)算所得機(jī)翼表面Cf分布云圖對(duì)比Fig.13 Comparison of Cf distribution of wing surface calculated by different simulation conditions

圖14 Re=3.30×106，不同模擬條件計(jì)算所得機(jī)翼不同站位Cf分布對(duì)比Fig.14 Re=3.30×106，Comparison of Cf distribution of different wing deployment stations calculated from different simulation conditions

圖15 Re=1.50×107，不同模擬條件計(jì)算所得機(jī)翼不同站位Cf分布對(duì)比Fig.15 Re=1.50×107，Comparison of Cf distribution of different wing deployment stations calculated from different simulation conditions

5 結(jié) 論

按照第一屆航空CFD可信度研討會(huì)-CHN-T1標(biāo)模CFD可信度研討會(huì)活動(dòng)要求，本文基于自研的非結(jié)構(gòu)混合網(wǎng)格解算器CARIA-OVERSET對(duì)運(yùn)輸機(jī)標(biāo)模CHN-T1進(jìn)行了初步的驗(yàn)證確認(rèn)計(jì)算研究，得到了以下基本結(jié)論：

(1) 同等網(wǎng)格節(jié)點(diǎn)數(shù)量，非結(jié)構(gòu)混合網(wǎng)格和結(jié)構(gòu)多塊網(wǎng)格計(jì)算所得氣動(dòng)力各項(xiàng)參數(shù)的差異很小。

(2) 對(duì)于網(wǎng)格尺度引起的誤差，非結(jié)構(gòu)混合網(wǎng)格大于結(jié)構(gòu)多塊網(wǎng)格，但非結(jié)構(gòu)網(wǎng)格得到的修正值與結(jié)構(gòu)網(wǎng)格接近，不同網(wǎng)格尺度對(duì)翼身結(jié)合處的分離特性模擬差異較小。

(3) 基于同一套網(wǎng)格，本文所采用的三種湍流模型計(jì)算結(jié)果表明，在固定升力條件下，湍流模型引起的摩擦阻力差量大于壓差阻力，且對(duì)迎角和俯仰力矩的影響較大，翼身結(jié)合處的分離模擬差異較小。

(4) 帶支撐構(gòu)型的全機(jī)升力和低頭力矩增加、阻力減小，力矩拐點(diǎn)提前，機(jī)翼靜氣彈變形后升力、阻力和低頭力矩均減小，力矩拐點(diǎn)后移。全湍和自由轉(zhuǎn)捩 模型的計(jì)算結(jié)果存在一定差異，但差異隨著雷諾數(shù)的增大而變小。

本文計(jì)算分析所得大部分結(jié)論具有一定的參考價(jià)值，但也存在些許不足，特別是基于自由轉(zhuǎn)捩模型的模擬計(jì)算，相關(guān)結(jié)論還需要更多可靠的試驗(yàn)數(shù)據(jù)作進(jìn)一步的驗(yàn)證分析。