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

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

0 引 言

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

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

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

1 數(shù)值方法

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

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

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

γ-θ模型求解兩個變量的標(biāo)準輸運方程：

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

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

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

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

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

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

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

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

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

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

4.1 計算收斂歷程曲線

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

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

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

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

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

R=ε21/ε32(3)

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

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

SC=S1-δ*(6)

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

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

(a) 迎角和俯仰力矩

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

表4是在Ma=0.78、CL=0.500±0.001時，不同形式、不同尺度網(wǎng)格計算所得的氣動力、力矩以及部分試驗數(shù)據(jù)?？梢钥闯觯诠潭ㄉl件下，計算迎角普遍小于試驗迎角；網(wǎng)格尺度對摩擦阻力和壓差阻力均有一定影響，摩擦阻力差量為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)格尺度對力矩系數(shù)影響較小，計算得到的抬頭力矩普遍大于試驗結(jié)果。

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

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

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

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

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

(a) 機翼

(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時，采用一方程的Spalart-Allmaras模型，兩方程的Wilcoxk-ω、SST共計三種湍流網(wǎng)格計算所得的氣動力、力矩以及部分試驗數(shù)據(jù)?？梢钥吹?，基于同一套非結(jié)構(gòu)混合網(wǎng)格，采用不同湍流模型計算所得結(jié)果存在一定差異，SA和SST兩種湍流模型的計算結(jié)果較為接近，SST計算所得阻力和迎角與試驗結(jié)果更加接近。

表5 湍流模型對氣動力的影響分析Table 5 Influence of turbulence models on aerodynamic forces

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

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

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

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

(2) 靜氣彈變形機翼帶支撐機構(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°時，干凈模型、帶支撐機構(gòu)模型和靜氣彈變形機翼帶支撐機構(gòu)模型計算所得機翼和平尾不同展向站位Cp分布曲線。由圖可以看出，支撐機構(gòu)主要影響平尾(與支撐方式有關(guān))，支撐機構(gòu)導(dǎo)致截面Cp環(huán)量減小，且越靠內(nèi)側(cè)影響越大；靜氣彈變形主要影響機翼，靜氣彈變形機翼上表面吸力減小，激波位置前移，且越靠外側(cè)影響越大。

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

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

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

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

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

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

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

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

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

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

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

圖14 Re=3.30×106，不同模擬條件計算所得機翼不同站位Cf分布對比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，不同模擬條件計算所得機翼不同站位Cf分布對比Fig.15 Re=1.50×107，Comparison of Cf distribution of different wing deployment stations calculated from different simulation conditions

5 結(jié) 論

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

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

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

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

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

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