轉(zhuǎn)捩對(duì)某民機(jī)超臨界機(jī)翼影響的數(shù)值模擬研究

李鵬飛, 冷智輝, 魏 闖, 任園軍, 梁 斌

(1.中國航空工業(yè)洪都航空工業(yè)集團(tuán) 飛機(jī)設(shè)計(jì)研究所, 江西 南昌 330024； 2.中國航空工業(yè)空氣動(dòng)力研究院 高速高雷諾數(shù)重點(diǎn)實(shí)驗(yàn)室, 遼寧 沈陽 110034)

0 引 言

超臨界機(jī)翼因具有優(yōu)越的氣動(dòng)性能得到了很大的發(fā)展，被廣泛應(yīng)用到民機(jī)機(jī)翼設(shè)計(jì)中，為民機(jī)的發(fā)展起到了巨大的推動(dòng)作用。與常規(guī)翼型相比，超臨界翼型[1]的上翼面平坦、下翼面后部內(nèi)凹，在跨聲速條件下，超臨界機(jī)翼表面不可避免的會(huì)出現(xiàn)邊界層發(fā)展與轉(zhuǎn)捩、激波/邊界層干擾、激波誘導(dǎo)分離等粘性支配的復(fù)雜流動(dòng)現(xiàn)象，決定了超臨界機(jī)翼氣動(dòng)性能對(duì)粘性的影響較為敏感。

轉(zhuǎn)捩表示邊界層由層流流動(dòng)向湍流流動(dòng)的轉(zhuǎn)變，是流體力學(xué)研究領(lǐng)域中一個(gè)十分重要而又非常復(fù)雜的難題?？缏曀贄l件下，轉(zhuǎn)捩對(duì)超臨界機(jī)翼的邊界層特性、激波位置、流動(dòng)分離位置等有十分重要的影響[2-4]，盡管19世紀(jì)80年代Reynolds就發(fā)現(xiàn)了轉(zhuǎn)捩現(xiàn)象，但至今理論上還無法準(zhǔn)確地確定轉(zhuǎn)捩點(diǎn)的位置，通過半個(gè)多世紀(jì)的理論與實(shí)驗(yàn)研究，對(duì)轉(zhuǎn)捩機(jī)理的認(rèn)識(shí)不斷深入，發(fā)展了許多邊界層轉(zhuǎn)捩判斷方法，同時(shí)隨著計(jì)算流體力學(xué)技術(shù)和計(jì)算機(jī)的發(fā)展以及工程應(yīng)用的需求，用CFD方法進(jìn)行轉(zhuǎn)捩數(shù)值模擬在飛行器氣動(dòng)設(shè)計(jì)、流動(dòng)細(xì)節(jié)及其機(jī)理分析中的權(quán)重日益提高，成為解決工程實(shí)際問題的重要手段[5-8]。預(yù)測附面層轉(zhuǎn)捩的方法很多，直 接 數(shù) 值 模 擬(DNS)[9]和大渦模擬(LES)[10]耗費(fèi)計(jì)算資源較大，目前主要用于簡單流場的模擬，工程上常用的是基于線性穩(wěn)定性理論的eN方法[11-12]和轉(zhuǎn)捩模型方法[13-14]。鑒于eN方法預(yù)測轉(zhuǎn)捩的成熟性，國外很多研究機(jī)構(gòu)將eN方法應(yīng)用其開發(fā)的CFD軟件中，如法國宇航院(ONERA)的ELSA[15-16]軟件，德國宇航院(DLR)的結(jié)構(gòu)網(wǎng)格求解器FlOWer[17-18]和非結(jié)構(gòu)混合網(wǎng)格的求解器TAU[19]中。國內(nèi)西北工業(yè)大學(xué)張坤和宋文萍等采用eN方法，通過精確求解線性穩(wěn)定性方程來判斷二維翼型和無限展長后掠機(jī)翼[20-21]邊界層轉(zhuǎn)捩，中國航空工業(yè)空氣動(dòng)力研究院董軍、高德峰、任園軍等應(yīng)用eN方法進(jìn)行了ONERA M6機(jī)翼的工程轉(zhuǎn)捩計(jì)算的研究[22]，國內(nèi)公開文獻(xiàn)中對(duì)eN方法的研究主要集中在二維翼型和三維機(jī)翼標(biāo)模上，而把eN方法預(yù)測轉(zhuǎn)捩用于飛機(jī)型號(hào)風(fēng)洞試驗(yàn)現(xiàn)象研究的文獻(xiàn)還不多見。

本文利用RANS方程耦合eN-Database轉(zhuǎn)捩預(yù)測方法，針對(duì)某民機(jī)超臨界機(jī)翼，研究轉(zhuǎn)捩對(duì)其激波位置、分離位置等流動(dòng)特性的影響，分析其在自由轉(zhuǎn)捩試驗(yàn)中升力線小迎角下出現(xiàn)拐折現(xiàn)象的原因及其流動(dòng)機(jī)理，供超臨界機(jī)翼設(shè)計(jì)和氣動(dòng)特性分析參考使用。

1 數(shù)值計(jì)算方法

1.1 eN-Database轉(zhuǎn)捩預(yù)測方法

本文數(shù)值模擬以可壓縮雷諾平均Navier-Stokes方程為主控方程，湍流模型采用在航空計(jì)算中有較好精度的的S-A一方程模型，空間離散采用Jameson中心差分加人工粘性，時(shí)間推進(jìn)采用TVD型的顯示三步Runge-Kutta推進(jìn)，通過多重網(wǎng)格技術(shù)、隱式殘差光順等方法加速收斂，并將法宇航Perraul、Arnal和Casalis等提出的流向和橫向不穩(wěn)定性eN-Database判斷方法耦合到求解器中來實(shí)現(xiàn)對(duì)轉(zhuǎn)捩點(diǎn)位置的判斷。求解過程為，首先進(jìn)行全湍條件下的流場求解，然后將壓力分布和網(wǎng)格信息輸入到轉(zhuǎn)捩判斷程序中，獲得轉(zhuǎn)捩位置，最后將轉(zhuǎn)捩位置回傳到流場求解器中進(jìn)行求解，即可得到自由轉(zhuǎn)捩計(jì)算結(jié)果，也可人為給定轉(zhuǎn)捩位置回傳到流場求解器中實(shí)現(xiàn)固定轉(zhuǎn)捩流場計(jì)算，該方法詳細(xì)介紹見文獻(xiàn)[22]。

1.2 計(jì)算網(wǎng)格

模型為某民機(jī)翼身組合體巡航構(gòu)型，機(jī)翼為超臨界機(jī)翼。網(wǎng)格為結(jié)構(gòu)網(wǎng)格形式，H-O型網(wǎng)格拓?fù)浣Y(jié)構(gòu)，邊界層第一層網(wǎng)格高度為平均氣動(dòng)弦長的2×10-6倍，保證物面yplus<1，向外遞增的倍數(shù)不大于1.2，空間網(wǎng)格過度均勻，盡可能保證相鄰網(wǎng)格之間的比值小于1.4，總網(wǎng)格量約為1070萬節(jié)點(diǎn)，飛機(jī)表面網(wǎng)格見圖1。

圖1 計(jì)算網(wǎng)格Fig.1 Calculational grid

2 計(jì)算結(jié)果及分析

在中航工業(yè)氣動(dòng)院FL-3風(fēng)洞進(jìn)行某民機(jī)構(gòu)型的高速測力風(fēng)洞試驗(yàn)中，發(fā)現(xiàn)在自由轉(zhuǎn)捩試驗(yàn)狀態(tài)下，在0°迎角附近，升力曲線發(fā)生了明顯的拐折現(xiàn)象，而在固定轉(zhuǎn)捩試驗(yàn)時(shí)，升力曲線拐折現(xiàn)象消失，但同迎角時(shí)升力減小，為研究產(chǎn)生此異常現(xiàn)象的原因，應(yīng)用eN方法對(duì)風(fēng)洞試驗(yàn)?zāi)Ｐ瓦M(jìn)行數(shù)值模擬計(jì)算。

2.1 計(jì)算和試驗(yàn)結(jié)果對(duì)比分析

圖2給出了計(jì)算和風(fēng)洞試驗(yàn)結(jié)果升力系數(shù)對(duì)比曲線，風(fēng)洞試驗(yàn)是在中航工業(yè)氣動(dòng)院FL-3風(fēng)洞完成，固定轉(zhuǎn)捩計(jì)算和試驗(yàn)給定機(jī)翼上翼面轉(zhuǎn)捩位置為7%當(dāng)?shù)叵议L處，下翼面自由轉(zhuǎn)捩，計(jì)算和試驗(yàn)狀態(tài)：Ma=0.8，基于機(jī)翼平均氣動(dòng)弦長Re=1.75×106。從中可以看出，全湍計(jì)算升力系數(shù)最小，同一升力下阻力最大，固定轉(zhuǎn)捩計(jì)算升力線相比全湍計(jì)算結(jié)果基本成向上平移趨勢，自由轉(zhuǎn)捩和固定轉(zhuǎn)捩下機(jī)翼氣動(dòng)性能差別很大，與固定轉(zhuǎn)捩相比，自由轉(zhuǎn)捩試驗(yàn)升力系數(shù)較大，阻力系數(shù)較小，計(jì)算結(jié)果與此趨勢一致，只是計(jì)算與試驗(yàn)在量值上略有區(qū)別，說明本文計(jì)算方法是可靠的。

從風(fēng)洞試驗(yàn)結(jié)果可以看出，自由轉(zhuǎn)捩時(shí)升力曲線約在迎角-1°到0°之間發(fā)生拐折，自由轉(zhuǎn)捩計(jì)算預(yù)測到了這種現(xiàn)象，約在迎角0°到1°之間升力曲線發(fā)生拐折，固定轉(zhuǎn)捩試驗(yàn)和計(jì)算升力曲線在迎角0°附近拐折現(xiàn)象消失。

(a) α～CL

(b) CD～CL

2.2 自由轉(zhuǎn)捩升力線拐折流動(dòng)機(jī)理分析

圖3給出了迎角-2°～2°固定轉(zhuǎn)捩和自由轉(zhuǎn)捩計(jì)算機(jī)翼上翼面層流區(qū)范圍。從可以看出，迎角-2°～0°自由轉(zhuǎn)捩計(jì)算層流區(qū)范圍明顯大于固定轉(zhuǎn)捩，亦即轉(zhuǎn)捩位置靠后，結(jié)合圖2可以看出，此時(shí)自由轉(zhuǎn)捩和固定轉(zhuǎn)捩計(jì)算升力系數(shù)差值較大；迎角0°和1°之間，自由轉(zhuǎn)捩計(jì)算轉(zhuǎn)捩位置發(fā)生突變前移，迎角1°和2°時(shí)自由轉(zhuǎn)捩計(jì)算層流區(qū)接近固定轉(zhuǎn)捩結(jié)果，此時(shí)兩者計(jì)算升力系數(shù)差值較小。綜合以上分析，機(jī)翼上翼面層流區(qū)越大，升力系數(shù)越大，迎角0°和1°之間，機(jī)翼上翼面轉(zhuǎn)捩位置發(fā)生突變，升力曲線因而發(fā)生拐折。

文獻(xiàn)[23]指出對(duì)于一個(gè)光滑平板的邊界層，由層流轉(zhuǎn)變?yōu)橥牧鞯霓D(zhuǎn)捩Re數(shù)一般在2×106～5×106之間，具體數(shù)值取決于噪聲水平和湍流度。這樣，當(dāng)試驗(yàn)?zāi)Ｐ偷娘L(fēng)洞Re數(shù)處于這個(gè)量級(jí)以及翼面具有平坦的壓力分布時(shí)，翼面上會(huì)存在明顯的層流邊界層。圖4分別給出了機(jī)翼三個(gè)不同站位剖面0°和1°迎角下固定轉(zhuǎn)捩和自由轉(zhuǎn)捩計(jì)算摩阻系數(shù)和壓力系數(shù)對(duì)比曲線。從中可以看出，迎角0°時(shí)機(jī)翼上翼面激波上游壓力分布較為平坦，自由轉(zhuǎn)捩條件下為層流，內(nèi)側(cè)機(jī)翼因上翼面后部出現(xiàn)激波誘導(dǎo)轉(zhuǎn)捩為湍流，其轉(zhuǎn)捩位置比固定轉(zhuǎn)捩位置明顯靠后，而外側(cè)機(jī)翼上翼面并未出現(xiàn)激波，轉(zhuǎn)捩位置更靠后，轉(zhuǎn)捩位置越靠前，激波位置就越靠前，超聲速區(qū)面積隨之減小，升力因而減小，這也是迎角0°時(shí)，固定轉(zhuǎn)捩下升力系數(shù)較自由轉(zhuǎn)捩小很多的原因。迎角1°時(shí)，機(jī)翼上翼面前緣出現(xiàn)明顯超聲速區(qū)，壓力分布成尖峰型，前緣附近具有很強(qiáng)的逆壓梯度，在低Re數(shù)下層流邊界層就能很快轉(zhuǎn)變?yōu)橥牧鳎滢D(zhuǎn)捩點(diǎn)接近前緣，略靠后于固定轉(zhuǎn)捩位置，激波位置同樣略靠后，但差別很小，各個(gè)站位壓力分布同樣差別很小，因而升力系數(shù)差別很小，這樣就造成自由轉(zhuǎn)捩條件下，升力線在0°發(fā)生拐折。

圖3 固定轉(zhuǎn)捩和自由轉(zhuǎn)捩計(jì)算機(jī)翼上翼面層流區(qū)范圍Fig.3 Comparison of laminar region and turbulentregion between fixed and free transition

圖4 α=0°和1°固定轉(zhuǎn)捩和自由轉(zhuǎn)捩摩阻和壓力分布對(duì)比曲線Fig.4 Comparison of skin friction and pressure distribution between fixed and free transition

2.3 轉(zhuǎn)捩對(duì)超臨界機(jī)翼邊界層流動(dòng)影響分析

圖5分別給出了迎角-2°～1°固定轉(zhuǎn)捩和自由轉(zhuǎn)捩計(jì)算機(jī)翼上翼面極限流線對(duì)比結(jié)果。固定轉(zhuǎn)捩時(shí)產(chǎn)生湍流邊界層/激波相互作用，迎角-2°時(shí)固定轉(zhuǎn)捩內(nèi)翼段后緣有小分離，隨迎角增加分離區(qū)大小逐漸增加，并向外翼段擴(kuò)展。自由轉(zhuǎn)捩機(jī)翼上翼面流動(dòng)變化較為復(fù)雜，迎角-2°～0°時(shí)，激波前為層流區(qū)域，產(chǎn)生轉(zhuǎn)捩邊界層/激波相互作用，在0°和1°之間，轉(zhuǎn)捩位置突變前移至靠近前緣位置，此時(shí)產(chǎn)生湍流邊界層/激波相互作用。迎角-2°時(shí)，內(nèi)翼段后緣有小分離，激波強(qiáng)度較固定轉(zhuǎn)捩強(qiáng)而且轉(zhuǎn)捩邊界層抵抗激波逆壓梯度的能力較弱，從而誘導(dǎo)產(chǎn)生激波后分離氣泡，隨迎角增加，激波誘導(dǎo)分離氣泡和后緣分離弦向范圍逐漸增加，至0°時(shí)，兩者匯合產(chǎn)生激波后緣完全分離。迎角1°時(shí)，轉(zhuǎn)捩位置突變前移至機(jī)翼前緣附近，湍流邊界層承受激波逆壓梯度能力較強(qiáng)，激波后分離氣泡消失，僅存在后緣分離，此時(shí)自由和固定轉(zhuǎn)捩機(jī)翼上翼面流動(dòng)基本相同。

圖5 固定和自由轉(zhuǎn)捩計(jì)算機(jī)翼上翼面極限流線Fig.5 Limited streamline distribution alongthe upper surface of the wing

3 結(jié) 論

本文采用RANS方程耦合eN-Database轉(zhuǎn)捩預(yù)測方法，研究了轉(zhuǎn)捩對(duì)某民機(jī)超臨界機(jī)翼氣動(dòng)特性的影響。與固定轉(zhuǎn)捩相比，自由轉(zhuǎn)捩狀態(tài)下計(jì)算得到的升力系數(shù)更大，相同升力系數(shù)下計(jì)算得到的阻力系數(shù)較小，與試驗(yàn)結(jié)果趨勢一致，進(jìn)一步驗(yàn)證了本文計(jì)算方法的可靠性。對(duì)于文中的民機(jī)構(gòu)型，可以得出以下結(jié)論：相同迎角下，轉(zhuǎn)捩位置越靠前，激波越靠前，超聲速區(qū)面積減小，升力隨之減??；在迎角0°～1°之間，機(jī)翼上翼面自由轉(zhuǎn)捩位置發(fā)生急劇移動(dòng)，是自由轉(zhuǎn)捩升力曲線發(fā)生拐折的原因。

參考文獻(xiàn)：

[1]Whitcom B, Richard T.Review of NASA supercritical airfoils[R].ICAS Paper No.74-10, 1974.

[2]沈克揚(yáng).關(guān)于超臨界機(jī)翼Re數(shù)效應(yīng)的機(jī)理研究[J].民用飛機(jī)設(shè)計(jì)研究, 1998, 3: 7-16.

[3]Chen Y C, Si J T, Han X L, et al.Investigation of transition effect on the pressure distribution of super-critical wing[J].Acta Aerodynamica Sinica, 2003, 21(4): 470-475.(in Chinese)陳迎春, 司江濤, 韓先鋰, 等.轉(zhuǎn)捩對(duì)超臨界機(jī)翼壓力分布的影響分析[J].空氣動(dòng)力學(xué)學(xué)報(bào), 2003, 21(4): 470-475.

[4]Xu X, Liu D W, Chen D H, et al.Numerical investigation on shock-induced separation structure of supercritical airfoil[J].Advanced Materials Research, 2013, 756-759: 4502-4505.

[5]Wang F, Eriqitai, Wang Q, et al.Boundary layers stability analysis of three typical airfoils[J].Acta Aerodynamica Sinica, 2011, 29(4): 481-485.(in Chinese)王菲, 額日其太, 王強(qiáng), 等.三種典型翼型邊界層穩(wěn)定性對(duì)比分析[J].空氣動(dòng)力學(xué)學(xué)報(bào), 2011, 29(4): 481-485.

[6]Deng L, Qiao Z D, Yang X D, et al.Computational aerodynamic analysis of thick flatback airfoils coupling RANS and transition prediction code[J].Acta Aerodynamica Sinica, 2011, 29(5): 613-618.(in Chinese)鄧?yán)? 喬志德, 楊旭東, 等.基于N-S方程和自由轉(zhuǎn)捩預(yù)測耦合求解的鈍后緣翼型氣動(dòng)性能計(jì)算[J].空氣動(dòng)力學(xué)學(xué)報(bào), 2011, 29(5): 613-618.

[7]Huang Z F, Lu X Z, Yu G T.Crossflow instability analysis and transition prediction of airfoil boundary layer[J].Acta Aerodynamica Sinica, 2014, 32(1): 14-20.(in Chinese)黃章峰, 逯學(xué)志, 于高通.機(jī)翼邊界層的橫流穩(wěn)定性分析和轉(zhuǎn)捩預(yù)測[J].空氣動(dòng)力學(xué)學(xué)報(bào), 2014, 32(1): 14-20.

[8]Feng W L, Chen B, Lyu L Y.Research Reynolds number effect correction for air supercritical laminar complex wing[J].Acts Aerodynamica Sinica, 2015, 33(4): 470-474.(in Chinese)馮文梁, 陳斌, 呂凌英.超臨界層流復(fù)合翼飛機(jī)Re數(shù)效應(yīng)修正方法研究[J].空氣動(dòng)力學(xué)學(xué)報(bào), 2015, 33(4): 470-474.

[9]Muppidi S, Mahesh K.DNS of transition in supersonic boundary layers.AIAA-2011-3564[R].Reston: 2011.

[10]Pan H L, Ma H D, Shen Q.LES application to unsteady flat plate shock wave/turbulent boundary layer interaction[J].Acts Aeronautica et Astronautica Sinica, 2011, 32: 242-248.(in Chinese)潘宏祿, 馬漢東, 沈清.基于LES方法的平板非定常激波/湍流邊界層干擾研究[J].航空學(xué)報(bào), 2011, 32(2): 242-248.(in Chinese)[11]Smith A M O, Gamberoni N.Transition pressure gradient and stability theory[R].Douglas AirCraft Co., CA Rept.ES 26388, Long Beach, CA, 1956.

[12]VanIngen J L.A suggested semi-empirical method for the calculation of the Boundary-Layer transition region[R].Univ.of Delft.Rept.VTH-74, Dept.of Aerospace Engineering, 1956.

[13]Wang G, Liu Y, Wang G Q, et al.Transitional flow simulation based onγ-Reθttransition model[J].Acta Aeronautica et Astronautica Sinica, 2014, 01: 70-79.(in Chinese)王剛, 劉毅, 王光秋, 等.采用γ-Reθt模型的轉(zhuǎn)捩流動(dòng)計(jì)算分析[J].航空學(xué)報(bào), 2014, 01: 70-79.

[14]Mend D H, Zhang Y L, Wang G X, et al.Application ofγ-Reθttransition model to two-dimensional low speed flows[J].Acts Aeronautica et Astronautica Sinica, 2011, 32(5): 792-801.(in Chinese)孟德虹, 張玉倫, 王光學(xué), 等.γ-Reθt轉(zhuǎn)捩模型在二維低速問題中的應(yīng)用[J].航空學(xué)報(bào), 2011, 32(5): 792-801.

[15]Jean-Pierre Archambaud, Raffaele Donelli.Automatic transition prediction using simplified methods[C]//47th AIAA Aerospace Sciences Meeting Including The New Horizons Forum and Aerospace Exposition, 2009.

[16]Cambier L, Heib S, Plot S.ONERA else CFD input from research and feed-back from industry[R].28th International Congress of The aeronautical Sciences, 2012.

[17]Andreas Krumbein.Automatic transition prediction and application to 3D wing configuration [J].Journal of Aircraft, 2007, 44(1): 119-133.

[18]Andreas Krumbein.Automatic transition prediction and application to hight-lift multi-element configurations[J].Journal of Aircraft, 2004, 4(5): 1362-1366.

[19]Normann Krimmelbein, Andreas Krumbein.Automatic transition prediction for three-dimensional configurations with focus on industrial application[C]//40th Fluid Dynamics Conference and Exhibit, 2010.

[20]Zhang K, Song W P.Application of the full eNtransition method to the infinite swept-wing’s transition prediction[J].Journal of Northwestern Polytechnical University, 2011, 01: 142-147.(in Chinese)張坤, 宋文萍.eN方法在無限展長后掠翼邊界層轉(zhuǎn)捩判斷中的初步應(yīng)用[J].西北工業(yè)大學(xué)學(xué)報(bào), 2011, 01: 142-147.

[21]Song Wenping, Zhu Zhen, Zhang Kun, et al.Simulations of the viscous flow around swept wings and optimization design using the RANS solver with automatic transition prediction[J].Aeronautical Science & Technology, 2015, 26(11): 23-29.(in Chinese)宋文萍, 朱震, 張坤, 等.耦合轉(zhuǎn)捩自動(dòng)判斷的機(jī)翼粘性繞流計(jì)算與優(yōu)化設(shè)計(jì)[J].航空科學(xué)技術(shù), 2015, 26(11): 23-29.

[22]Dong J, Gao D F, Ren Y J, et al.Coupling and appli-cation of eN-database method to transition prediction in a 3D unstructured solver[J].Journal of Shenyang Aerospace University, 2015, 32(2): 11-17.(in Chinese)董軍, 高德峰, 任園軍, 等.eN-Database轉(zhuǎn)捩預(yù)測方法在三維非結(jié)構(gòu)求解器中的耦合與應(yīng)用[J].沈陽航空航天大學(xué)學(xué)報(bào), 2015, 32(2): 11-17.

[23]Elsenaar A, Binion TW, Stamewsky E.Reynolds number effects in transonic flows[R].AGARD AG-303, December 1988.