CHN-T1標(biāo)模的數(shù)值計算及氣動特性研究

李浩然, 李亞坤, 張宇飛, 陳海昕

(清華大學(xué) 航天航空學(xué)院, 北京 100084)

0 引 言

計算流體力學(xué)(Computational Fluid Dynamics，CFD)作為一項重要的流體分析手段，在近30年間得到了飛速的發(fā)展，尤其是在諸如飛行器氣動設(shè)計等外流問題中，CFD已成為初步設(shè)計、詳細(xì)設(shè)計和計算校核的主要方法，與風(fēng)洞試驗一起成為飛行器設(shè)計和研發(fā)的重要研究手段[1]。目前CFD技術(shù)在求解可壓縮流動的飛機(jī)流場問題中，各類算法、求解器、前處理與后處理軟件系統(tǒng)都已趨于成熟；在精細(xì)化湍流模擬[2]、氣動聲學(xué)[3]等問題中也有了長足的進(jìn)步。

在進(jìn)行飛行器氣動設(shè)計時，CFD 的可信度研究十分必要，也是當(dāng)前的研究熱點(diǎn)，其基本內(nèi)容和方法就是CFD的驗證和確認(rèn)(verification & validation)[4]。1998年，AIAA發(fā)布了世界上第一個關(guān)于CFD驗證與確認(rèn)的指南[5]。國際上組織了很多CFD的驗證和確認(rèn)專題研討會議，比如歐洲的計算空氣動力學(xué)研究項目ECARP(European Computational Aerodynamics Research Project)，該項目的目標(biāo)是評估各種湍流模型，并對軟件系統(tǒng)進(jìn)行統(tǒng)一的研究[6]。2001年6月至2016年7月期間，AIAA阻力計算的工作小組組織了六次阻力預(yù)測會議(Drag Prediction Workshop)[7]。其中DPW2會議選擇德國宇航局設(shè)計的跨聲速民用飛機(jī)模型DLR-F6作為實驗和CFD的基準(zhǔn)模型，其巡航馬赫數(shù)為0.75，巡航升力系數(shù)為0.5，2008年法國的S2MA風(fēng)洞對DLR-F6進(jìn)行了風(fēng)洞試驗； DPW4、DPW5、DPW6會議選擇CRM(Common Research Model)作為基準(zhǔn)模型，由于具有公認(rèn)的試驗數(shù)據(jù)，并且經(jīng)過了廣泛的驗證，CRM模型的計算已成為檢驗CFD程序?qū)Φ湫涂缏曀俟r預(yù)測能力的重要算例。此外，第一屆AIAA計算流體力學(xué)高升力裝置預(yù)測會議發(fā)布了標(biāo)準(zhǔn)的HILIFTPW-1高升力預(yù)測模型[8]；中國航空研究院在2012年設(shè)計了空氣動力學(xué)驗證模型CAE-AVM風(fēng)洞模型，該模型在荷蘭的DNW-HST跨聲速風(fēng)洞進(jìn)行了試驗測試，一方面驗證模型的氣動性能，另一方面也對CFD數(shù)值計算方法進(jìn)行校核[9]。

為了評估國內(nèi)CFD技術(shù)狀態(tài)，為飛行器研制提供技術(shù)參考，中國空氣動力研究與發(fā)展中心(CARDC)開展首屆CFD可信度專題研討會(AeCW-1)，并設(shè)計了單通道客機(jī)CHN-T1標(biāo)模[10]，用于風(fēng)洞試驗和CFD可信度研究。

本文工作結(jié)合AeCW-1研討會的計算任務(wù)，研究了SST與SA兩種湍流模式在標(biāo)模計算中的網(wǎng)格收斂性與抖振特性的計算能力；分別從標(biāo)模的抖振特性、靜氣彈效應(yīng)、雷諾數(shù)影響方面評估了CFD計算的可信度以及CHN-T1標(biāo)模的設(shè)計氣動特性。

1 計算狀態(tài)與計算方法

1.1 CHN-T1標(biāo)模簡介

AeCW-1的計算構(gòu)型有三種，其中Config.1為 “機(jī)身-機(jī)翼-平尾-垂尾”(BWHV)構(gòu)型,是基礎(chǔ)構(gòu)型，具體的參數(shù)見表1，其中Sref為全模的參考面積、Cref為參考弦長、xref/yref/zref為力矩參考點(diǎn)； Config.2在Config.1基礎(chǔ)上加上風(fēng)洞試驗的支撐；Config.3在Config.2基礎(chǔ)上考慮了機(jī)翼的靜氣彈變形。CHN-T1標(biāo)模的實驗馬赫數(shù)為0.78，雷諾數(shù)為3.3×106，升力系數(shù)為0.5。

表1 CHN-T1 BWHV構(gòu)型參數(shù)Table 1 The parameter of BWHV configuration of CHN-T1

1.2 計算狀態(tài)

本文涵蓋了部分AeCW-1研討會的計算任務(wù)，包括三個計算狀態(tài)，其中Case 1為網(wǎng)格收斂性研究，采用Config.1構(gòu)型，同時對比SST與SA湍流模式計算結(jié)果；Case 2b與Case2c為抖振特性研究，分別采用Config.2與Config.3構(gòu)型，并對比SST與SA的計算結(jié)果；Case 3為雷諾數(shù)的影響研究，采用SST湍流模式計算Config.3構(gòu)型。具體的計算狀態(tài)參見表2及文獻(xiàn)[11]。

表2 計算狀態(tài)Table 2 Computational State

1.3 計算方法

本文計算工作基于開源的CFL3D Version6.7結(jié)構(gòu)網(wǎng)格求解器[12]，該求解器由NASA蘭利研究中心開發(fā)于20世紀(jì)90年代，內(nèi)含多種湍流模式及計算格式，使用MPI進(jìn)行大規(guī)模并行計算。

對標(biāo)模的計算基于定常RANS的方法，時間離散采用隱式LU-ADI格式依當(dāng)?shù)貢r間步長進(jìn)行推進(jìn)；空間離散的重構(gòu)步采用三階迎風(fēng)MUSCL格式，其中限制器類型選擇iflim=4[12]，演化步采用無熵修正的Roe的FDS(Flux Differential Splitting，通量差分裂)格式[13]；采用三重網(wǎng)格W循環(huán)加速技術(shù)；采用自由來流邊界進(jìn)行全湍流計算，湍流模式選取航空領(lǐng)域廣泛應(yīng)用的SST[14]與SA[15]模式。

2 網(wǎng)格收斂性研究

網(wǎng)格收斂性是評價計算精度及網(wǎng)格質(zhì)量高低的一種策略，本文采用組委會提供的粗、中、細(xì)網(wǎng)格[11]進(jìn)行網(wǎng)格收斂性研究，三套網(wǎng)格及部件節(jié)點(diǎn)分布見圖1。機(jī)翼上的網(wǎng)格量對于捕捉激波的影響較大，三套網(wǎng)格在機(jī)翼弦向上的網(wǎng)格點(diǎn)數(shù)分別為61、81、121；y+的值對于摩擦阻力的計算尤為關(guān)鍵，在不含壁函數(shù)的SST與SA湍流模型中，由y+控制的第一層網(wǎng)格高度應(yīng)在粘性底層以內(nèi)，以便更好地描述邊界層速度型，三套網(wǎng)格y+的值分別約為1.0、0.6、0.4。隨著網(wǎng)格的加密，網(wǎng)格間距造成的數(shù)值耗散減小，對于機(jī)頭、機(jī)尾等部件的阻力預(yù)測更加準(zhǔn)確。

(a) Coarse grid

(b) Medium grid

(c) Fine grid

圖2展示了采用SST湍流模型計算得到的密度殘差收斂歷史，其中粗網(wǎng)格的收斂速度較慢，這是由于在該計算中使用了較小的CFL數(shù)以保證計算穩(wěn)定。

網(wǎng)格收斂性期望得到一條迎角AOA或力系數(shù)隨橫坐標(biāo)N-2/3/10-5的線性走勢，該式中N為網(wǎng)格量。圖3(a)展示了總阻力系數(shù)的網(wǎng)格收斂性，發(fā)現(xiàn)隨著網(wǎng)格的加密，SST與SA計算的阻力值都呈下降趨勢，在從粗網(wǎng)格到中網(wǎng)格時兩種湍流模式的阻力變化量基本一致，而從中網(wǎng)格到細(xì)網(wǎng)格時SA的變化量大于SST。相比之下SST計算的阻力隨網(wǎng)格量的變化更接近線性關(guān)系，并且利用外插得到的無限密網(wǎng)格上的收斂阻力值與組委會提供的104億網(wǎng)格參考值[11]較為接近。圖3(b)展示了壓差阻力CDp系數(shù)的網(wǎng)格收斂性，發(fā)現(xiàn)SST與SA都具有較好的收斂特性，并且兩者計算的壓差阻力較為接近。圖3(c)展示了摩擦阻力CDf的網(wǎng)格收斂性，在粗網(wǎng)格與中網(wǎng)格上SA的預(yù)測結(jié)果比SST大了約6 counts(1 count=0.0001)，而在細(xì)網(wǎng)格上兩種湍流模式的預(yù)測值接近。圖3(d)展示了力矩特性的網(wǎng)格收斂性，力矩特性并沒有體現(xiàn)出單調(diào)的收斂趨勢。總的來說對于組委會提供的網(wǎng)格而言，SST的網(wǎng)格收斂性表現(xiàn)優(yōu)于SA。

圖2 粗、中、細(xì)三套網(wǎng)格的殘差收斂歷史Fig.2 History of residual in coarse, medium and fine grid

圖3 氣動力系數(shù)的網(wǎng)格收斂性Fig.3 Grid convergence with aerodynamic coefficients

隨著網(wǎng)格的加密，數(shù)值模擬對于精細(xì)流動結(jié)構(gòu)的捕捉能力逐漸增強(qiáng)。圖4、圖5分別展示了SST與SA湍流模式計算η=40%、80%兩個展向位置的機(jī)翼表面壓力系數(shù)Cp分布，三套網(wǎng)格計算的壓力分布形態(tài)基本一致；隨著網(wǎng)格的加密，激波位置不變但激波處的壓力梯度越來越大；對η=80%位置吸力平臺的高度計算：用SST計算從粗到細(xì)網(wǎng)格時，平臺高度略微抬高，而采用SA模式計算的高度基本一致；采用SA模式計算的細(xì)網(wǎng)格下η=80%位置激波后的壓力恢復(fù)區(qū)出現(xiàn)明顯抖動。

對于翼身結(jié)合部位流動現(xiàn)象的模擬是標(biāo)模計算的一個重點(diǎn)，網(wǎng)格量的大小會影響該處分離流的計算[16]。圖6展示了翼身結(jié)合部位置的流線圖，發(fā)現(xiàn)兩種湍流模式在粗網(wǎng)格及細(xì)網(wǎng)格上都捕捉到了角隅渦，并且隨著網(wǎng)格的加密，分離區(qū)的大小幾乎不發(fā)生變化。

綜合考量氣動力系數(shù)以及流動結(jié)構(gòu)的捕捉情況，本文在后續(xù)的抖振特性、靜氣彈效應(yīng)、雷諾數(shù)影響研究中選用中網(wǎng)格進(jìn)行計算。

圖4 CHN-T1 Case1: SST計算的截面壓力分布的網(wǎng)格收斂性Fig.4 CHN-T1 Case1 grid convergence with sectional pressure distribution using SST

圖5 CHN-T1 Case1: SA計算的截面壓力分布的網(wǎng)格收斂性Fig.5 CHN-T1 Case1 grid convergence with sectional pressure distribution using SA

圖6 CHN-T1 Case1: 翼身結(jié)合部的角隅渦Fig.6 CHN-T1 Case1: horseshoe vortex in wing-body configuration

3 抖振特性研究

超臨界機(jī)翼在跨聲速飛行時可能會產(chǎn)生由激波/附面層相互干擾而引起抖振現(xiàn)象，這將嚴(yán)重影響飛行器的飛行品質(zhì)與結(jié)構(gòu)壽命[17]。發(fā)生抖振時機(jī)翼上表面的激波以及其后的分離泡位置變動，屬于典型的非定常現(xiàn)象，但是利用定常的RANS計算結(jié)合一定判據(jù)可以較好的確定飛機(jī)的抖振邊界[18]。

本文計算了在馬赫數(shù)固定的情況下，隨著迎角的不斷增大CHN-T1標(biāo)模逐漸進(jìn)入抖振的過程。圖7展示了利用SST與SA湍流模式計算得到的升力線，并與風(fēng)洞實驗值[19]對比，發(fā)現(xiàn)升力線都在迎角4°后偏離線性走勢。圖8展示了極曲線，其計算值與實驗值吻合。圖9反映了力矩系數(shù)隨升力系數(shù)的變化，在升力系數(shù)小于0.4時SA的計算結(jié)果與實驗值吻合較好，而升力系數(shù)大于0.4時SST的計算結(jié)果與實驗值吻合的更好。適航規(guī)章要求民用運(yùn)輸類飛機(jī)在正常使用狀態(tài)下不得超過抖振發(fā)生邊界[20]，當(dāng)Cm～CL曲線發(fā)生彎折即可判定發(fā)生抖振。后續(xù)對于抖振狀態(tài)的流場分析選擇SST計算的Case2b。

對抖振特性的研究從分離區(qū)及激波強(qiáng)度兩個方面開展。首先探索分離區(qū)的變化，選取迎角為3°、3.75°、4.25°觀察機(jī)翼表面流線分布，如圖10所示，在迎角3°時機(jī)翼表面幾乎不出現(xiàn)分離；當(dāng)迎角增大到3.75°時能明顯觀察到激波/邊界層干擾從而產(chǎn)生的分離，其中分離線的位置緊挨激波位置，再附線的位置約為50%弦長；當(dāng)迎角到4.25°時機(jī)翼表面出現(xiàn)了大面積分離，機(jī)翼中部的流線分布較為紊亂，此時已進(jìn)入了危險的抖振區(qū)域。以圖11中4個典型的截面為例分析進(jìn)入抖振狀態(tài)的壓力分布特點(diǎn)，如圖12所示，發(fā)現(xiàn)隨著迎角增大，各截面的前緣吸力峰值不斷提高，進(jìn)而導(dǎo)致了激波的增強(qiáng)；外翼段的后加載顯著下降，使得進(jìn)入抖振狀態(tài)的低頭力矩減小，Cm～CL曲線發(fā)生彎折，操縱特性變差，在展向η=40%、70%位置處的后緣的壓力恢復(fù)劇烈，使得邊界層更容易出現(xiàn)分離。

圖7 CHN-T1 Case2b/ Case2c抖振研究:升力線Fig.7 CHN-T1 Case2b/ Case2c buffet study: lift curve

圖8 CHN-T1 Case2b/ Case2c抖振研究:極曲線Fig.8 CHN-T1 Case2b/ Case2c buffet study: drag polar

圖9 CHN-T1 Case2b/ Case2c抖振研究:力矩特性線Fig.9 CHN-T1 Case2b/ Case2c buffet study: pitching moment curve

探索抖振狀態(tài)下激波隨迎角的變化，圖13展示了迎角為0°、2°、3°、4.25°時表面壓力分布及提取的空間激波，其中激波判據(jù)為:

式中:U為速度矢量，a為聲速，p為壓力梯度矢量，為壓力梯度矢量的二范數(shù)，該式表示了無量綱速度在壓力梯度方向上的投影，當(dāng)S≥1.0則認(rèn)為出現(xiàn)激波[21]。 從S=1.0的等值面圖中可以發(fā)現(xiàn)在0°迎角時機(jī)翼上表面無激波出現(xiàn)，而迎角2°時機(jī)翼中段出現(xiàn)激波，隨著迎角不斷增加，激波范圍擴(kuò)展至內(nèi)外翼段。

(a) AOA=3.0°

(b) AOA=3.75°

(c) AOA=4.25°

圖11 CHN-T1 Config.2機(jī)翼展向典型截面位置示意圖Fig.11 CHN-T1 Config.2 typical sections on the wing platform

4 靜氣彈效應(yīng)研究

AeCW-1組委會提供了經(jīng)過靜氣動彈性變形的標(biāo)模構(gòu)型Config.3，圖14展示了經(jīng)過靜氣彈變形后的機(jī)翼形變及扭轉(zhuǎn)角沿展向分布的情況，迎角自-2°變化至4.25°時，變形出現(xiàn)在Kink至翼尖區(qū)域，變形后外翼段扭轉(zhuǎn)角減小。

圖12 CHN-T1 Case2b:典型截面壓力分布隨迎角的變化Fig.12 CHN-T1 Case2b: the change of pressure distribution on typical sections with AOA

圖13 CHN-T1 Case2b:表面壓力分布及激波隨迎角變化Fig.13 CHN-T1 Case2b : the change of pressure distribution and shock with AOA

圖14 CHN-T1 Case3: 靜氣彈變形后的機(jī)翼幾何變化Fig.14 CHN-T1 Case3: the geometry of wing after the elastic deformation

圖15為2°、3°、3.75°、4.25°迎角下靜氣彈效應(yīng)對機(jī)翼展向升力系數(shù)分布的影響，在定迎角計算時，靜氣彈變形使得整機(jī)升力系數(shù)有所降低，同時由于Kink以外位置的扭轉(zhuǎn)角減小，在2°、3°、3.75°迎角下變形后的機(jī)翼外翼段升力系數(shù)明顯降低，而進(jìn)入抖振狀態(tài)的迎角4.25°下升力系數(shù)的變化并不明顯。

截取機(jī)翼展向85%位置(η=85%)研究靜氣彈效應(yīng)對壓力分布的影響，如圖16所示，在2°、3°、3.75°

圖15 CHN-T1 Case2b/Case2c:靜氣彈效應(yīng)對機(jī)翼展向升力系數(shù)分布的影響Fig.15 CHN-T1 Case2b/Case2c: the elastic effect for distribution of lift coefficient in span direction

圖16 CHN-T1 Case2b/Case2c: η=85%位置截面壓力分布Fig.16 CHN-T1 Case2b/Case2c: the sectional pressure distribution on η=85%

以及迎角4.25°下，靜氣彈變形將導(dǎo)致前緣吸力峰、吸力平臺下降，從一定程度上削弱了激波。

5 雷諾數(shù)影響研究

雷諾數(shù)不同，通常會對轉(zhuǎn)捩點(diǎn)的位置、邊界層厚度以及激波位置等產(chǎn)生影響，從而導(dǎo)致飛機(jī)氣動特性的變化，文獻(xiàn)[22]中指出運(yùn)輸機(jī)模型的翼型厚度和彎度較大，其所受雷諾數(shù)的影響較大。本文以AeCW-1組委會要求的Case 3為例，研究雷諾數(shù)對標(biāo)模氣動特性的影響。

表3展示了Case 3的結(jié)果，雷諾數(shù)分別取為3.3×106和15.0×106，采用SST湍流模式進(jìn)行定升力CL=0.5的計算，得到的總阻力與實驗的誤差在6 counts之內(nèi)；統(tǒng)計了計算得到的壓差阻力與摩擦阻力，發(fā)現(xiàn)當(dāng)雷諾數(shù)增加4.5倍后總阻力下降了47.8 counts，其中壓差阻力下降了17.1 counts，摩擦阻力下降了29.7 counts。

表3 Case 3結(jié)果Table 3 Computational Results for Case 3

分別對壓差阻力和摩擦阻力進(jìn)行分析。圖17展示了主翼展向η=20%、40%、70%、85%截面的壓力分布，發(fā)現(xiàn)在大雷諾數(shù)下后加載Cp包圍的面積增大，同時機(jī)翼上表面超聲速區(qū)的吸力平臺有所下降，進(jìn)而導(dǎo)致了激波強(qiáng)度的減弱、壓差阻力減小。

圖17 CHN-T1 Case3:不同雷諾數(shù)下典型截面的壓力分布Fig.17 CHN-T1 Case3: typical sectional pressure distribution with different Reynolds number

圖18展示了兩雷諾數(shù)下的壁面摩擦系數(shù)(Cf)分布圖，不同雷諾數(shù)下Cf分布的差別主要體現(xiàn)在機(jī)身、主翼、平尾、垂尾的前緣。大雷諾數(shù)下前緣的摩擦系數(shù)明顯下降。

圖18 CHN-T1 Case3:不同雷諾數(shù)下壁面摩擦系數(shù)分布Fig.18 CHN-T1 Case3: the contour of skin friction coefficient with different Reynolds number

6 結(jié) 論

本文采用結(jié)構(gòu)網(wǎng)格求解器CFL3D計算了CHN-T1標(biāo)模的典型工況，研究了SST與SA湍流模式在網(wǎng)格收斂性以及抖振特性上的表現(xiàn)，分析了靜氣彈效應(yīng)及雷諾數(shù)對于標(biāo)模氣動特性的影響，主要結(jié)論如下：

(1) 在網(wǎng)格收斂性方面，采用CFL3D的SST模式表現(xiàn)優(yōu)于SA模式；隨著網(wǎng)格加密激波位置不變，激波強(qiáng)度逐漸增大；采用粗網(wǎng)格與細(xì)網(wǎng)格同樣能夠捕捉到翼身結(jié)合部的角隅渦。

(2) 對抖振特性的研究，采用SST模式計算的Case2b與實驗值吻合最好，CHN-T1標(biāo)模的抖振源自機(jī)翼中段出現(xiàn)的激波/邊界層干擾產(chǎn)生的分離，當(dāng)迎角不斷增大時，各截面壓力分布的前緣吸力峰值不斷提高，導(dǎo)致了激波的增強(qiáng)，外翼段的各截面壓力分布后加載下降，使得進(jìn)入抖振狀態(tài)的低頭力矩減小，Cm～CL曲線發(fā)生彎折。

(3) 靜氣彈效應(yīng)主要影響了標(biāo)模Kink位置以外的升力系數(shù)沿展向分布，靜氣彈變形后的機(jī)翼外翼段升力系數(shù)明顯降低，同時該處激波被削弱致使阻力降低。

(4) 雷諾數(shù)增大時，標(biāo)模壓差阻力的減小主要來源于壓力系數(shù)后加載區(qū)包圍的面積增大以及超聲速吸力平臺區(qū)高度的降低，對于摩擦阻力的影響主要表現(xiàn)在標(biāo)模各個部件前緣的摩擦系數(shù)下降。