洪俊武，王運(yùn)濤，李偉,2, *，楊小川

1. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000 2. 國防科技大學(xué) 空天科學(xué)學(xué)院，長沙 410073

目前，盡管計算流體力學(xué)(Computational Fluid Dynamics,CFD)方法和計算機(jī)技術(shù)都獲得了長足進(jìn)步與發(fā)展，但基于雷諾平均Navier-Stokes(Reynolds Averaged Navior-Stokes,RANS)方程模擬復(fù)雜高升力構(gòu)型的可信度依然有待提高，采用主流CFD方法獲得高升力構(gòu)型的最大升力系數(shù)及失速迎角，并準(zhǔn)確預(yù)測和模擬分離流動的產(chǎn)生和發(fā)展，依然是CFD面臨的巨大挑戰(zhàn)[1-2]。為了對高升力構(gòu)型的流動機(jī)理進(jìn)行研究并深入展開CFD軟件對此類構(gòu)型的驗(yàn)證與確認(rèn)工作，許多CFD可信度專題會議也以此作為研究主題[3-7]，其中，最具代表性的會議就是AIAA高升力構(gòu)型性能預(yù)測會議。

2017年6月，第3屆AIAA高升力構(gòu)型(HiLiftPW-3)預(yù)測會議在美國丹佛召開。HiLiftPW-3會議選擇的2組標(biāo)模分別是NASA(National Aeronautics and Space Administration)提供的HL-CRM(HiLift-Common Research Model)標(biāo)模和JAXA(Japan Aerospace eXploration Agency)提供的JSM(Jaxa hilift configuration Standard Model)標(biāo)模。會議組委會確定的研究工況分別是：針對HL-CRM標(biāo)模，開展網(wǎng)格收斂性研究(Case1)；針對JSM標(biāo)模，開展氣動力特性對比研究(Case2)；針對DSMA661二維翼型，開展湍流模型的標(biāo)定(Case3)。

Case1的主要內(nèi)容是對HL-CRM標(biāo)模進(jìn)行典型狀態(tài)的網(wǎng)格收斂性研究，組委會沒有提供計算狀態(tài)的風(fēng)洞試驗(yàn)結(jié)果，目的是對各CFD軟件開展驗(yàn)證工作。Case2的主要內(nèi)容是對JSM標(biāo)模有/無掛架短艙的2種構(gòu)型進(jìn)行模擬，并將計算結(jié)果與風(fēng)洞試驗(yàn)數(shù)據(jù)進(jìn)行對比，目的是考察掛架短艙對氣動特性的影響，確認(rèn)CFD結(jié)果的模擬精度。Case3的主要內(nèi)容是研究DSMA661翼型壓力分布、速度型分布和尾跡區(qū)流動，主要目的是驗(yàn)證所采用湍流模型的有效性。本文重點(diǎn)討論Case1和Case2的數(shù)值模擬結(jié)果。

采用TRIP3.0軟件，王運(yùn)濤等[8-10]研究了網(wǎng)格規(guī)模、湍流模型、轉(zhuǎn)捩模型等多種因素對TrapWing全展長襟翼高升力構(gòu)型數(shù)值模擬結(jié)果的影響，洪俊武等[11]采用拼接結(jié)構(gòu)網(wǎng)格技術(shù)模擬了TrapWing半展長襟翼高升力構(gòu)型的氣動特性和壓力分布特性。采用五階空間離散精度的WCNS(Weighted Compact Nonlinear Scheme)格式，王運(yùn)濤[12-13]和李松[14]等，開展了典型運(yùn)輸機(jī)高升力構(gòu)型數(shù)值模擬技術(shù)研究。此外，國內(nèi)一些CFD工作者[15-20]采用結(jié)構(gòu)網(wǎng)格技術(shù)、非結(jié)構(gòu)網(wǎng)格技術(shù)對高升力構(gòu)型開展了一系列數(shù)值模擬技術(shù)研究工作。

本文作者作為HiLiftPW-3的參研團(tuán)隊之一(參會序號為015)，利用自主開發(fā)的CFD軟件平臺TRIP3.0，根據(jù)HiLiftPW-3兩種高升力構(gòu)型標(biāo)模的外形和流場特征，采用拼接結(jié)構(gòu)網(wǎng)格技術(shù)和有限體積方法，通過求解三維曲線坐標(biāo)系下的RANS方程，對會議要求的Case1和Case2工況進(jìn)行了數(shù)值模擬。研究結(jié)果表明，本文方法對典型運(yùn)輸機(jī)三段翼布局的低速問題模擬具有良好的適用性，在失速迎角之前，數(shù)值模擬得到的氣動力系數(shù)和壓力分布均與試驗(yàn)結(jié)果高度吻合，同時也取得了合理的網(wǎng)格收斂性結(jié)果。

1 計算方法

本項研究采用有限體積法和結(jié)構(gòu)網(wǎng)格技術(shù)求解曲線坐標(biāo)系下的RANS方程。離散方程組的求解采用LU-SGS[21]方法；無黏項離散采用MUSCL格式[22]，并利用Roe格式進(jìn)行通量分裂；黏性項采用中心差分格式進(jìn)行計算。計算中，湍流模型采用目前應(yīng)用較為廣泛的SA(Spalart-Almaras)一方程模型[23]，該模型計算量相對較小，而且具有良好的魯棒性。

采用了多重網(wǎng)格技術(shù)、低速預(yù)處理技術(shù)、并行技術(shù)加速收斂。并行計算采用域分解法進(jìn)行數(shù)據(jù)劃分，并使用基于消息傳遞的并行編程模型來設(shè)計開發(fā)并行計算模塊，迭代過程中采用MPI消息傳遞庫，通過顯式地發(fā)送和接收消息來完成各處理器之間的數(shù)據(jù)交換，具有較好的通用性和可移植性，可以獲得較高的并行效率。實(shí)際計算中針對不同算例，在國產(chǎn)飛騰CPU上采用了16～640個CPU進(jìn)行并行運(yùn)算。

采用了拼接網(wǎng)格技術(shù)(Surface to Surface)來處理前緣縫翼和后緣襟翼偏轉(zhuǎn)帶來的剪刀縫，針對此類拓?fù)浣Y(jié)構(gòu)，拼接網(wǎng)格技術(shù)可以有效提高網(wǎng)格質(zhì)量并大大降低網(wǎng)格規(guī)模，同時能在拼接面上保證通量守恒，本文的拼接網(wǎng)格方法在此前的大量應(yīng)用中也已經(jīng)得到了充分驗(yàn)證。

2 HL-CRM標(biāo)模網(wǎng)格收斂性

HL-CRM是NASA提供的高升力構(gòu)型標(biāo)模，該構(gòu)型是翼身組合體布局，機(jī)翼采用三段構(gòu)型，多段翼之間無固定連接裝置。HL-CRM構(gòu)型的前緣縫翼與后緣襟翼的偏角分別為30°和37°，前緣縫翼的展向長度基本相當(dāng)于機(jī)翼展長，從翼根延伸到翼尖；后緣襟翼的展向長度大約是機(jī)翼展長的2/3，從翼根延伸到展向外側(cè)，該構(gòu)型為典型的著陸構(gòu)型。

2.1 計算網(wǎng)格

本文在研究中首先采用商業(yè)軟件生成了一套基準(zhǔn)的多塊拼接網(wǎng)格,并在這套基準(zhǔn)的中等規(guī)模(Medium)網(wǎng)格基礎(chǔ)上，對網(wǎng)格進(jìn)行重構(gòu)，獲得粗化的Tiny、Coarse網(wǎng)格和細(xì)化的Fine網(wǎng)格。這4套不同密度網(wǎng)格的單元規(guī)模分別為432萬，1 458萬、3 455萬和11 661萬，各套網(wǎng)格的具體參數(shù)見表1，中等網(wǎng)格表面網(wǎng)格及截面網(wǎng)格見圖1。

圖1 HL-CRM標(biāo)橫網(wǎng)格Fig.1 Grid for HL-CRM model

表1 HL-CRM標(biāo)模網(wǎng)格參數(shù)Table 1 Grid parameters of HL-CRM model

2.2 網(wǎng)格收斂性

本節(jié)采用前述的4套計算網(wǎng)格開展了網(wǎng)格收斂性研究，來流條件為：馬赫數(shù)Ma=0.2，迎角α=8°、16°，基于平均氣動力弦長(MAC=7 005.32 mm)的雷諾數(shù)為Re=3.26×106。

圖2給出了極粗、粗、中、細(xì)4套不同網(wǎng)格密度下升力系數(shù)CL和平均殘差RESave的迭代收斂歷程，其中橫坐標(biāo)為迭代步數(shù)(Iteration)，所有網(wǎng)格均采用三重網(wǎng)格進(jìn)行計算?？梢钥吹?，4套網(wǎng)格的計算結(jié)果均獲得了良好的收斂性，極粗網(wǎng)格(Tiny)收斂速度最快，細(xì)網(wǎng)格(Fine)收斂最慢；極粗網(wǎng)格約3 000步殘差降到最低，粗網(wǎng)格需要5 000 步，中等網(wǎng)格需要接近10 000步，而細(xì)網(wǎng)格需要20 000步左右。由于網(wǎng)格越密，近壁面區(qū)的網(wǎng)格尺度也越小，導(dǎo)致附面層網(wǎng)格單元的當(dāng)?shù)貢r間步長也越小，這就使得流場收斂速度大大降低。

圖2 不同網(wǎng)格密度下的升力系數(shù)和平均殘差收斂歷程Fig.2 Convergence history of lift coefficient and average residual under different grid densities

圖3給出了采用不同密度網(wǎng)格得到的HL-CRM氣動特性隨網(wǎng)格單元數(shù)量的變化曲線，圖中橫坐標(biāo)中N是網(wǎng)格單元數(shù)目?？梢钥吹剑煌窍律ο禂?shù)(CL)、阻力系數(shù)(CD)和俯仰力矩系數(shù)(Cm)隨網(wǎng)格密度的增加都是單調(diào)變化的，計算結(jié)果具有良好的網(wǎng)格收斂性。另外，從圖3中可以看到，α=16°時粗網(wǎng)格的升力系數(shù)和力矩系數(shù)都和密網(wǎng)格結(jié)果相差較大，而α=8°時粗、密網(wǎng)格的升力系數(shù)和力矩系數(shù)差量則小得多。這是由于小迎角下，機(jī)翼周圍流場變化相對較小，此時的粗網(wǎng)格基本能捕捉到該狀態(tài)下空間流場的主要特征，所以粗網(wǎng)格的計算結(jié)果和密網(wǎng)格結(jié)果差距不大；而大迎角下，機(jī)翼周圍流場變化較為劇烈，一些區(qū)域甚至出現(xiàn)了小的分離氣泡，此時的粗網(wǎng)格由于尺度過大已經(jīng)無法再現(xiàn)真實(shí)流場的主要特征，因此和密網(wǎng)格結(jié)果差距明顯。

圖3 Case1氣動力系數(shù)的網(wǎng)格收斂性曲線Fig.3 Grid convergence curves of aerodynamic coefficients for Case1

3 JSM標(biāo)模有/無掛架短艙對比

JSM是JAXA提供的高升力構(gòu)型標(biāo)模，該構(gòu)型也是翼身組合體布局，機(jī)翼采用三段翼構(gòu)型，前緣縫翼通過8個小連接塊與主翼固定，后緣襟翼通過3個大的滑軌與主翼固定。另外，該構(gòu)型分為無發(fā)動機(jī)掛架短艙的Case2a和有發(fā)動機(jī)掛架短艙的Case2c，本次會議主要研究該構(gòu)型在有/無短艙情況下氣動力特性與試驗(yàn)的對比以及有/無短艙情況下二者的氣動特性差異。該構(gòu)型的前緣縫翼與后緣襟翼的偏角分別是為25°和30°，前緣縫翼的展向長度基本相當(dāng)于機(jī)翼展長，從翼根延伸到翼尖；后緣襟翼的展向長度大約是機(jī)翼展長的3/4，從翼根延伸到展向外側(cè)，該構(gòu)型也是典型的著陸構(gòu)型。

3.1 計算網(wǎng)格

Case2計算網(wǎng)格的拓?fù)浣Y(jié)構(gòu)與Case1類似，Case2a(無短艙)和Case2c(有短艙)2套網(wǎng)格的具體參數(shù)見表2。圖4給出了Case2a構(gòu)型在風(fēng)洞中的安裝照片和計算網(wǎng)格，圖5給出了Case2c構(gòu)型在風(fēng)洞中的安裝照片和計算網(wǎng)格。由于Case2包含滑軌和多個連接部件，構(gòu)型復(fù)雜程度比Case1高得多，因此網(wǎng)格量也大得多。

表2 JSM標(biāo)模網(wǎng)格參數(shù)Table 2 Grid parameters of JSM model

圖4 Case2a的風(fēng)洞安裝照片和計算網(wǎng)格Fig.4 Installation picture in wind tunnel and computational grid for Case2a

圖5 Case2c的風(fēng)洞安裝照片和計算網(wǎng)格Fig.5 Installation picture in wind tunnel and computational grid for Case2c

3.2 Case2a無掛架短艙構(gòu)型計算

Case2的來流條件為：Ma=0.172，α=4.36°,10.47°,14.54°,18.58°,20.59°,21.57°，基于平均氣動力弦長(MAC=529.2 mm)的雷諾數(shù)為Re=1.93×106。圖6給出了JSM無掛架短艙(Nacelle/Pylon off)構(gòu)型氣動特性隨迎角變化曲線。

從升力曲線來看，計算和試驗(yàn)的失速迎角十分吻合，都在20.59°附近；在19.59°之前，二者幾乎重合；失速迎角附近，計算得到的升力系數(shù)略高。作者認(rèn)為，這個差異是試驗(yàn)?zāi)Ｐ偷膹椥宰冃卧斐傻?。由于試?yàn)?zāi)Ｐ驮诖笥菚r翼尖上翹明顯，造成升力系數(shù)偏低；而計算模型采用的是剛性外形，這使得計算值相對試驗(yàn)結(jié)果偏高。

從阻力曲線來看，計算結(jié)果比試驗(yàn)值明顯偏高，一方面本文采用了SA模型進(jìn)行全湍流模擬，而SA的數(shù)值模擬結(jié)果通常會導(dǎo)致阻力略微偏大；另一方面這也是本屆高升力會議出現(xiàn)的一個普遍現(xiàn)象，圖7是本屆高升力會議組織方發(fā)布的Case2兩組升阻比曲線統(tǒng)計結(jié)果，可以看到，除個別明顯異常的結(jié)果外，30余家科研機(jī)構(gòu)提供的40余組升阻比曲線全部位于試驗(yàn)曲線的右側(cè)，阻力普遍偏大，所有數(shù)值結(jié)果均存在不同程度的平移。與這一現(xiàn)象不同的是，第1屆高升力構(gòu)型會議中主流軟件的阻力計算結(jié)果均與試驗(yàn)值高度吻合[7]；第2屆會議中，絕大多數(shù)軟件獲得的阻力計算結(jié)果也明顯超過試驗(yàn)值，而且迎角越大差異越明顯[24]。雖然高升力構(gòu)型的研究重點(diǎn)并不在于阻力的精確模擬，但是從數(shù)據(jù)對比來看，這個問題也是CFD工作者需要密切關(guān)注的重點(diǎn)問題之一。

從力矩曲線來看，計算結(jié)果也與試驗(yàn)值吻合良好。產(chǎn)生這一現(xiàn)象的主要原因是該構(gòu)型沒有平尾，所以全機(jī)力矩的主要來源變成了機(jī)翼，這也意味著只要計算得到的機(jī)翼壓力分布與試驗(yàn)吻合的話，升力和力矩系數(shù)都應(yīng)該是吻合的。下文的壓力分布對比也證明了這一點(diǎn)。

圖6 Case2a氣動力系數(shù)計算與試驗(yàn)結(jié)果的比較Fig.6 Comparison of aerodynamic coefficients of calculation and test results for Case2a

圖8為組委會給出的Case2a所有參會結(jié)果與試驗(yàn)結(jié)果的氣動特性對比，圖中黑色為試驗(yàn)結(jié)果，紅色為參會的全部數(shù)值模擬結(jié)果，綠色曲線為本文結(jié)果?？梢钥吹?，本文結(jié)果具有良好的可信度。

圖9給出了α=10.47°時，典型機(jī)翼站位的壓力系數(shù)(Cp)分布計算值與試驗(yàn)結(jié)果對比，4個展向站位分別是η=16%(A-A)，41%(D-D)，56%(E-E)和89%(H-H)?？梢钥吹?，該迎角下，不同站位前緣縫翼、主翼和后緣襟翼上壓力分布的試驗(yàn)和計算值高度吻合，除個別試驗(yàn)跳點(diǎn)以外，幾乎重合在一起。

圖7 Case2升力隨阻力變化曲線統(tǒng)計結(jié)果與試驗(yàn)結(jié)果的比較Fig.7 Comparison of lift-drag curves between statistical and test results for Case2

圖8 Case2a所有參會結(jié)果與試驗(yàn)的比較Fig.8 Comparison of results for Case2a including all simulations and test

圖9 Case2a α=10.47°典型站位Cp分布對比Fig.9 Comparison of distribution of Cp at typical station α=10.47° for Case2a

圖10給出了α=18.58°典型站位的Cp分布計算值與試驗(yàn)結(jié)果對比，可以看到，該迎角下位于翼根和機(jī)翼中部3個不同站位的試驗(yàn)和計算值高度吻合，幾乎重合；而接近翼尖的89%站位試驗(yàn)值與計算值差異明顯。產(chǎn)生這一現(xiàn)象的原因前面已經(jīng)闡述了，這是由于風(fēng)洞試驗(yàn)使用的是彈性模型，數(shù)值計算模擬的是剛性模型，試驗(yàn)的彈性模型在大迎角下翼尖變形量較大，因此Cp差異較為明顯；而小迎角下的所有站位和大迎角下的非翼尖站位變形量都較小，Cp分布都吻合良好。

圖10 Case2a α=18.58°典型站位Cp分布對比Fig.10 Comparison of distribution of Cp at typical station α=18.58° for Case2a

3.3 Case2c有掛架短艙構(gòu)型計算

圖11給出了Case2c有掛架短艙(Nacelle/Pylon on)構(gòu)型氣動特性隨迎角變化曲線。從升力曲線來看，計算和試驗(yàn)的失速迎角較為吻合，試驗(yàn)結(jié)果在20.09°附近，計算結(jié)果略為提前，失速迎角為19.59°。在19.59°之前，二者幾乎重合。從阻力曲線來看，計算結(jié)果也比試驗(yàn)值偏高，與Case2a規(guī)律一致；從力矩曲線來看，在失速迎角之前，計算結(jié)果也與試驗(yàn)值吻合良好。圖12為組委會給出的Case2c所有參會計算結(jié)果與試驗(yàn)結(jié)果的氣動特性對比，圖中綠色曲線為本文結(jié)果?？梢钥吹剑疚慕Y(jié)果具有良好的可信度。

圖13給出了α=10.47°典型機(jī)翼站位的Cp分布計算值與試驗(yàn)結(jié)果對比?？梢钥吹?，該迎角下不同站位的試驗(yàn)和計算值高度吻合，幾乎重合在一起。圖14給出了α=18.58°典型機(jī)翼站位的Cp分布計算值與試驗(yàn)結(jié)果對比，可以看到，該迎角下位于翼根和機(jī)翼中部3個不同站位的試驗(yàn)和計算值也高度吻合，而接近翼尖的站位試驗(yàn)值與計算值差異明顯，產(chǎn)生這一現(xiàn)象的原因與Case2a完全相同，此處不再贅述。

圖11 Case2c氣動力系數(shù)計算與試驗(yàn)結(jié)果的比較Fig.11 Comparison of aerodynamic coefficients of calculation and test results for Case2c

圖12 Case2c所有參會結(jié)果與試驗(yàn)結(jié)果的比較Fig.12 Comparison of results Case2c including all simulations and test

圖13 Case2c α=10.47°典型站位Cp分布對比Fig.13 Comparison of distribution of Cp at typical station α=10.47° for Case2c

圖14 Case2c α=18.58°典型站位Cp分布對比Fig.14 Comparison of distribution of Cp at typical station α=18.58° for Case2c

3.4 有/無掛架短艙構(gòu)型的氣動力差量對比

圖15給出了在失速迎角之前，有/無掛架短艙2種構(gòu)型氣動力差異的試驗(yàn)和計算結(jié)果對比。為方便比較，圖中給出了組委會提供的所有參會結(jié)果，綠色曲線為本文結(jié)果。

由于2種構(gòu)型的氣動力差量是小量，所以圖中本文的升力差量曲線與試驗(yàn)值雖然不完全吻合，但是2條曲線的平均差量只有0.02左右，二者實(shí)際上是很接近的。可以看到，小迎角情況下，試驗(yàn)和計算得到的各氣動力系數(shù)差量都相當(dāng)吻合，大迎角時差距有所增大，但也在合理范圍之內(nèi)。這也進(jìn)一步說明了，本文的CFD方法對于此類高升力構(gòu)型的模擬是非常適用的，無論是氣動力系數(shù)還是不同構(gòu)型之間氣動力系數(shù)的差量，計算和試驗(yàn)結(jié)果都具有良好的一致性。

圖15 Case2兩種構(gòu)型氣動力差量所有參會結(jié)果與試驗(yàn)結(jié)果的比較Fig.15 Comparison of aerodynamic dispersion between two configurations for Case2 including all simulations and test results

4 結(jié) 論

本文采用拼接結(jié)構(gòu)網(wǎng)格技術(shù)，通過求解曲線坐標(biāo)系下的RANS方程，數(shù)值模擬了第3屆高升力構(gòu)型會議提供的2組標(biāo)模。通過理論分析和與試驗(yàn)結(jié)果的對比，結(jié)論如下：

1) 本次對比研究中，風(fēng)洞模型翼尖在大迎角下的變形是導(dǎo)致試驗(yàn)結(jié)果和計算結(jié)果差異的主要來源之一。

2) 對此類三段翼典型布局的高升力構(gòu)型，數(shù)值方法在失速迎角之前，可以獲得可信度相當(dāng)高的流場結(jié)果。對比結(jié)果表明，無論是壓力分布還是氣動力系數(shù)，試驗(yàn)和計算都取得了良好的一致性。

3) 本文的數(shù)值實(shí)踐表明，TRIP軟件平臺采用的拼接結(jié)構(gòu)網(wǎng)格技術(shù)和計算方法對此類外形和低速流場具有良好的適用性；本文的計算結(jié)果在與全球40余組數(shù)值結(jié)果的各項對比中均名列前茅，在各大科研機(jī)構(gòu)和CFD軟件的同臺競技中表現(xiàn)優(yōu)秀，計算結(jié)果與試驗(yàn)結(jié)果合理地吻合，可以為大型飛機(jī)增升裝置的氣動設(shè)計提供可靠的技術(shù)支持。