HiLiftPW-3高升力構(gòu)型數(shù)值模擬

洪俊武，王運濤，李偉,2, *，楊小川

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

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

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

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

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

1 計算方法

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

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

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

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

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

2.1 計算網(wǎng)格

本文在研究中首先采用商業(yè)軟件生成了一套基準的多塊拼接網(wǎng)格,并在這套基準的中等規(guī)模(Medium)網(wǎng)格基礎(chǔ)上，對網(wǎng)格進行重構(gòu)，獲得粗化的Tiny、Coarse網(wǎng)格和細化的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標橫網(wǎng)格Fig.1 Grid for HL-CRM model

表1 HL-CRM標模網(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給出了極粗、粗、中、細4套不同網(wǎng)格密度下升力系數(shù)CL和平均殘差RESave的迭代收斂歷程，其中橫坐標為迭代步數(shù)(Iteration)，所有網(wǎng)格均采用三重網(wǎng)格進行計算?？梢钥吹?，4套網(wǎng)格的計算結(jié)果均獲得了良好的收斂性，極粗網(wǎng)格(Tiny)收斂速度最快，細網(wǎng)格(Fine)收斂最慢；極粗網(wǎng)格約3 000步殘差降到最低，粗網(wǎng)格需要5 000 步，中等網(wǎng)格需要接近10 000步，而細網(wǎng)格需要20 000步左右。由于網(wǎng)格越密，近壁面區(qū)的網(wǎng)格尺度也越小，導(dǎo)致附面層網(wǎ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ù)量的變化曲線，圖中橫坐標中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ù)差量則小得多。這是由于小迎角下，機翼周圍流場變化相對較小，此時的粗網(wǎng)格基本能捕捉到該狀態(tài)下空間流場的主要特征，所以粗網(wǎng)格的計算結(jié)果和密網(wǎng)格結(jié)果差距不大；而大迎角下，機翼周圍流場變化較為劇烈，一些區(qū)域甚至出現(xiàn)了小的分離氣泡，此時的粗網(wǎng)格由于尺度過大已經(jīng)無法再現(xiàn)真實流場的主要特征，因此和密網(wǎng)格結(jié)果差距明顯。

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

3 JSM標模有/無掛架短艙對比

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

3.1 計算網(wǎng)格

Case2計算網(wǎng)格的拓撲結(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標模網(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)型氣動特性隨迎角變化曲線。

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

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

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

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

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

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

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

圖8 Case2a所有參會結(jié)果與試驗的比較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分布計算值與試驗結(jié)果對比，可以看到，該迎角下位于翼根和機翼中部3個不同站位的試驗和計算值高度吻合，幾乎重合；而接近翼尖的89%站位試驗值與計算值差異明顯。產(chǎn)生這一現(xiàn)象的原因前面已經(jīng)闡述了，這是由于風(fēng)洞試驗使用的是彈性模型，數(shù)值計算模擬的是剛性模型，試驗的彈性模型在大迎角下翼尖變形量較大，因此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)型氣動特性隨迎角變化曲線。從升力曲線來看，計算和試驗的失速迎角較為吻合，試驗結(jié)果在20.09°附近，計算結(jié)果略為提前，失速迎角為19.59°。在19.59°之前，二者幾乎重合。從阻力曲線來看，計算結(jié)果也比試驗值偏高，與Case2a規(guī)律一致；從力矩曲線來看，在失速迎角之前，計算結(jié)果也與試驗值吻合良好。圖12為組委會給出的Case2c所有參會計算結(jié)果與試驗結(jié)果的氣動特性對比，圖中綠色曲線為本文結(jié)果?？梢钥吹剑疚慕Y(jié)果具有良好的可信度。

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

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

圖12 Case2c所有參會結(jié)果與試驗結(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)型氣動力差異的試驗和計算結(jié)果對比。為方便比較，圖中給出了組委會提供的所有參會結(jié)果，綠色曲線為本文結(jié)果。

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

圖15 Case2兩種構(gòu)型氣動力差量所有參會結(jié)果與試驗結(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ù)，通過求解曲線坐標系下的RANS方程，數(shù)值模擬了第3屆高升力構(gòu)型會議提供的2組標模。通過理論分析和與試驗結(jié)果的對比，結(jié)論如下：

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

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

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