靖振榮，孫朋朋，黃章峰*，2

(1.天津大學 力學系，天津 300072;2.中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室，綿陽 621000)

機翼作為航空飛行器上產生升力的重要部件之一，其設計的水平很大程度上決定了飛機性能的優(yōu)劣.在飛機飛行阻力中，摩擦阻力是一個重要的組成部分.對于亞聲速飛機，表面摩擦阻力占總阻力的50%左右［1］.使飛機在飛行時機翼表面盡可能長地保持層流狀態(tài)是降低摩擦阻力的最有效的辦法之一.因此，深入研究轉捩機理、預測轉捩位置、防止或推遲轉捩對機翼的優(yōu)化設計和氣動減阻有著十分重要的意義.

飛機一般采用后掠機翼，其目的是提高飛行的臨界馬赫數(shù)、克服因接近音速飛行而急劇增大的空氣阻力、提高飛行穩(wěn)定性.在后掠機翼頭部邊界層內，由于壓力梯度和后掠角的共同作用，使得與勢流方向相垂直的方向有速度分量，稱為橫流［2］.橫流失穩(wěn)是引起后掠機翼邊界層流動轉捩的主要因素［3］.攻角是指機翼翼弦與自由來流方向的夾角，是機翼產生升力的主要因素之一，一般不為零.因此，有攻角的后掠翼成為人們的研究熱點.

在國外，Boltz等［4］在低湍流風洞中對后掠機翼邊界層進行了研究，發(fā)現(xiàn)隨著攻角的減小，同一位置出現(xiàn)轉捩所需的雷諾數(shù)逐漸增大，轉捩提前發(fā)生在機翼的迎風面.Haynes［5］采用非線性拋物化方程(NPSE)研究了橫流渦的非線性特性和擾動幅值飽和特性，發(fā)現(xiàn)在小攻角下，橫流擾動波的流向增長范圍隨著攻角的增加幾乎呈線性增長.Dagenhart和 Saric［6］實 驗 發(fā) 現(xiàn) － 4°攻 角 的NLF(2)-0415后掠機翼會產生較合適的壓力梯度，對Tollmien-Schlichting(T-S)波有穩(wěn)定作用，對橫流渦的增長卻有促進作用，采用eN方法預測的轉捩位置的 N 值范圍為6.4～6.8.Reibert和Saric［7］對后掠機翼邊界層橫流失穩(wěn)及轉捩的實驗研究進行了綜述.Bushnell等［8］總結了線性理論預測轉捩位置的成果，認為在粗糙度和背景湍流度較小時，eN方法可以應用到馬赫數(shù)、攻角、壓力梯度、后掠變化等情況，而且當來流湍流度小于0.05%時轉捩對應的 N 值為 9 ～11.Arnal等［9］利用eN方法研究了有無抽吸時的后掠翼的轉捩，發(fā)現(xiàn)當橫流轉捩占主導時，N值為6～10，并且抽吸可以減小N值.

在國內，孫朋朋和黃章峰［10］研究了后掠角對橫流不穩(wěn)定性的影響，發(fā)現(xiàn)在機翼頭部靠近前緣位置，后掠角為45°時橫流強度最強，駐波(ω=0°)在后掠角為30°～45°時不穩(wěn)定區(qū)域的范圍最大.李峰等［11］采用數(shù)值模擬方法研究了低速大攻角下機翼的繞流特性，發(fā)現(xiàn)在某些馬赫數(shù)和攻角下NACA0012翼型的湍流解具有周期性.吳鋆、等［12-13］采用實驗方法研究了攻角對翼型繞流結構的影響，發(fā)現(xiàn)分離點和分離剪切層形成旋渦的位置隨迎角的增大而向上游移動，同時翼型上表面流動分離后形成的回流區(qū)尺寸隨著翼型迎角的增加而增大.袁湘江等［14］采用直接數(shù)值模擬(DNS)和線性穩(wěn)定性(LST)分析方法，對高超聲速小攻角鈍錐邊界層的失穩(wěn)機制和轉捩特點進行了研究，發(fā)現(xiàn)攻角為2°的高超聲速鈍錐邊界層存在多支不穩(wěn)定模態(tài).對于機翼邊界層穩(wěn)定性問題，王斌等［15］采用風洞試驗研究了攻角對后掠機翼轉捩位置的影響，發(fā)現(xiàn)機翼上表面轉捩位置隨著攻角的增大而前移，當達到某一攻角時，轉捩快速接近前緣.

因此攻角對橫流的產生和演化都有很大的影響，而且在后掠機翼頭部的迎風面與背風面邊界層內的壓力梯度有所不同，攻角對迎風面和背風面的橫流穩(wěn)定性影響也將不同.本文以適航的后掠機翼為研究對象，在給定后掠角及雷諾數(shù)的情況下，結合數(shù)值模擬和線性穩(wěn)定性理論，研究不同攻角對后掠機翼邊界層在機翼迎風面和背風面的流動穩(wěn)定性的影響，并采用轉捩預測的eN方法進行轉捩預測，為機翼的優(yōu)化減阻設計提供借鑒和參考.

1 數(shù)值計算方法

采用有限體積法，求解三維可壓縮無量綱Navier-Stokes方程:

式中:t為時間;x和y為流向和法向坐標;U為守恒型通量，由密度、3個方向的速度(流向、法向和展向)和溫度組成;E和F為非線性項通量，其中包含壓力項;Ev和 Fv為黏性項通量［16］.由于機翼在展向無限長，因此機翼和流動在不同的展向位置均相同，其在展向的導數(shù)恒為零.

無量綱黏性系數(shù)μ由Sutherland關系式確定:

式中:T 為溫度;C=110.4/T∞，T∞為來流溫度.

數(shù)值計算中，時間導數(shù)采用LU-SGS離散;對流項采用Steger-Warming格式進行通量分裂，采用3階WENO格式進行重構.網(wǎng)格外邊界采用遠場邊界條件;壁面采用無滑移絕熱邊界條件;展向采用周期邊界條件.

得到基本流場后，將Navier-Stokes方程寫成擾動形式，并進行線化就得到了經典的Orr-Sommerfeld方程.采用2點4階的Malik差分法，將Orr-Sommerfeld方程離散后得到特征關系式:

式中:α和β分別為擾動在流向和展向波數(shù);ω為擾動頻率;f為與基本流有關的量.對于空間模式;β、ω為實數(shù)，α為復數(shù)，其虛部 －αi為擾動的增長率.給定β、ω及f，采用Mueller法，就可以求出特征值α和特征函數(shù)φ.通過變動β、ω及f使αi=0，即擾動的增長率為零，從而得到中性曲線.對結果進行穩(wěn)定性分析，進而采用轉捩預測eN方法進行預測［17］.

2 計算模型和條件

采用NACA0012模型，弦長1 m，后掠角取25°.計算模型坐標示意圖如圖1所示.其中x、y、z為直角坐標，η、ξ為貼體坐標;Λ 為后掠角;為貼體坐標下的速度分量為勢流方向速度;為橫流速度.選取高空10 km處氣體參數(shù)作為參考量，溫度為223.3K，來流馬赫數(shù)為 0.7，單位雷諾數(shù)為 5.945 ×106/m.攻角 α 選取 0°、1°、2°、3°、4°、5°.若無注明，文中的量均為無量綱量，其中速度用來流速度U∞無量綱化，長度采用1mm無量綱化，相應的無量綱雷諾數(shù)為5945.

圖1 計算模型坐標示意圖Fig.1 Sketch map of the model coordinate

計算域及網(wǎng)格如圖2所示.其中，x方向為機翼弦長方向，y方向為機翼法向方向.機翼流向網(wǎng)格最小間距小于0.1mm，法向上一個邊界層內布有200個點以上，法向網(wǎng)格最小間距小于0.01mm.

對韓步璋等［18］實驗工況的基本流進行數(shù)值模擬，得到的壓力系數(shù)與實驗結果的對比如圖3所示.可以看出，計算結果與實驗結果符合得很好，說明本研究計算的結果可靠.

圖2 計算網(wǎng)格示意圖Fig.2 Sketch map of the computational domain

圖3 壓力系數(shù)沿弦長方向分布Fig.3 Distribution of pressure coefficient along chord

圖4給出了在橫流坐標系下，21.6萬、40萬、60萬3種網(wǎng)格模型在不同位置處的橫流速度剖面，其中，為邊界層外沿處的勢流速度，方向與邊界層外沿勢流方向垂直.可以看出，3種網(wǎng)格模型計算得到的結果基本一致，說明模型網(wǎng)格已經夠密.

圖4 橫流速度剖面Fig.4 Profile of the crossflow velo city

3 數(shù)值結果分析

3.1 基本流分析

圖5給出了不同攻角情況下機翼的壓力系數(shù)曲線.從曲線中可以看出，隨著攻角的增大，機翼兩側面的壓力系數(shù)之差逐漸增大，說明產生的升力逐漸增加.當攻角α為4°和 5°時，在 x/c=0.1～0.3位置，機翼背風面壓力系數(shù)出現(xiàn)突降現(xiàn)象.這是由于當攻角α≥4°時，在機翼的背風面一定范圍內會有激波出現(xiàn)［18］.

圖5 不同攻角下的壓力系數(shù)Fig.5 Pressure coefficients with different attack angles

選取數(shù)值模擬計算的流場截面，建立新的坐標系:ξ的大小定義為在翼型剖面上距離翼型前緣的壁面弧長距離，方向沿翼型剖面的弧線方向;η的大小定義為該點到壁面的距離，方向與壁面法向方向相同.重新劃分網(wǎng)格并加密處理，得到可用于穩(wěn)定性分析的基本流場.

圖6(a)和圖6(b)分別給出了不同攻角工況下，在ξ=100位置處迎風面和背風面的流向速度剖面.由圖可知，隨著攻角的增大，邊界層附近的流向速度在迎風面逐漸減小，在背風面逐漸增大.在η＞1處流向速度都接近于常值.在迎風面邊界層中流向速度小于背風面邊界層中的流向速度.

圖6 在ξ=100處的流向速度 剖面Fig.6 Profile of flow velocity atξ=100

圖7 在ξ=100處橫流速度剖面Fig.7 Profile of crossflow velocity at ξ=100

圖7(a)和圖7(b)對比了迎風面和背風面在橫流坐標系下，不同攻角工況時ξ=100位置處的橫流速度剖面.可以看出，橫流速度在其拐點處達到峰值，且隨著攻角的增大，峰值的大小在迎風面逐漸增大，在背風面逐漸減小.攻角不為零時，迎風面的橫流速度大于背風面的橫流速度.在迎風面，橫流速度的峰值均在η=0.2mm位置左右，而在背風面，橫流速度的峰值所在位置隨著攻角的增加而增加.

3.2 中性曲線

圖8(a)和圖8(b)給出了背風面和迎風面在不同攻角情況下，展向波數(shù)β=3的擾動波在(ξ，ω)的中性曲線(ω為擾動頻率).可以看出，在迎風面，隨著攻角的增大，展向波數(shù)β=3的擾動波的中性曲線的流向范圍逐漸增大.頻率較大的擾動波對攻角的變化不太敏感，而頻率較小的擾動波受到攻角的影響很大，而且對于攻角小于4°的小攻角，中性曲線在流向的最大位置隨著攻角的增加而線性增加，與 Haynes［5］的結果一致.在背風面，中性曲線的流向范圍隨攻角的增大而減小，而且中性曲線在流向的最大位置隨著攻角的增加而線性減小.隨著攻角的增加，在背風面出現(xiàn)了高頻的不穩(wěn)定波.迎風面和背風面在(ξ，ω)平面上的中性曲線與ξ坐標軸都有交點，而且不同的攻角時交點的位置不同，這說明不穩(wěn)定區(qū)域存在頻率ω=0的駐波，而且攻角對駐波的穩(wěn)定性影響很大.

圖8 中性穩(wěn)定曲線(β=3)Fig.8 Neutral stability curves(β =3)

圖9(a)和圖9(b)給出了迎風面和背風面在不同攻角下駐波(ω=0)在(ξ，β)平面上的中性曲線.可以看出，展向波數(shù)β較大的駐波的流向增長區(qū)域隨攻角的增大變化較小，而展向波數(shù)較小的波(1＜β＜4)的流向增長區(qū)域隨攻角的增大變化較大.隨著攻角的增大，迎風面駐波的流向增長區(qū)域逐漸增大，背風面駐波的流向增長區(qū)域逐漸減小，且在迎風面的增長范圍大于背風面.說明相對于背風面，迎風面更不穩(wěn)定.類似(ξ，ω)平面上的中性曲線，小攻角時(ξ，β)平面上迎風面/背風面的中性曲線在流向的最大位置隨著攻角的增加而線性增加/減小.

圖9 中性穩(wěn)定曲線(ω=0)Fig.9 Neutral stability curves(ω =0)

3.3 N值曲線及轉捩預測

圖10給出了迎風面和背風面在不同攻角工況下采用轉捩預測eN方法計算的N值曲線.可以看出，迎風面的N值相比背風面增長較快，且在迎風面的N值要遠大于背風面的N值.對比可知，當機翼攻角不為零時，擾動波在迎風面的N值最大，說明擾動波在迎風面?zhèn)鞑r被放大的倍數(shù)最大.如圖10所示，當機翼攻角為3°時，擾動波的N值可達到 10，即擾動波可被放大 e10倍(約2.2萬倍)，表明迎風面更容易發(fā)生轉捩.

攻角增大時，在迎風面，來流速度方向與迎風面壁面切線方向的夾角增大，氣流在壁面附近的偏轉角也增大，因此氣流沿著流向的速度減小，更多的流體流向展向方向，導致迎風面的橫流速度增強，而在背風面反之.從圖6可以看出，隨著攻角的增加，迎風面的流向速度逐漸減小，而背風面的流向速度逐漸增大，從圖7可以看出隨著攻角的增加，迎風面的橫流速度明顯增大，而背風面的橫流速度明顯減小.橫流不穩(wěn)定性是由于拐點產生的，其不穩(wěn)定性的強度與橫流速度的最大值及其位置息息相關.橫流速度的最大值越大，橫流不穩(wěn)定性越強.圖8和圖9中迎風面的中性曲線范圍所包含的面積要大于背風面，而且圖10中迎風面的N值要明顯大于背風面的N值，說明相對于背風面，迎風面更不穩(wěn)定.

圖10 eN方法的N值曲線Fig.10 N factor in eN method

圖11為在不同位置處，擾動波N值隨攻角的變化曲線，其中攻角為正值表示為迎風面，攻角為負值表示為背風面.可以看出，在同一位置處迎風面的N值隨著攻角的增大而增大，背風面的N值隨攻角的增大而減小.這說明在攻角不為零的情況下，擾動波的增長在迎風面得到了加強，在背風面受到抑制.攻角的大小與擾動波受到增強和抑制的程度有關，攻角越大，擾動波受到抑制和得到加強的程度越大.

圖11 N值隨攻角的變化曲線Fig.11 Curves of N factor changing with attack angle

研究表明，當擾動波的無量綱幅值達到來流的20%就會發(fā)生轉捩［19］.因此當飛行器在空中飛行時，機翼邊界層流動中擾動的初始幅值A0的大小對流動能否轉捩有著極其重要的影響.圖12給出了不同攻角下ξ位置處發(fā)生轉捩時的初始位置處(ξ=25)的擾動幅值A0.可以看出，隨著攻角的增大，迎風面引起轉捩所需的初始擾動幅值逐漸減小，說明轉捩越容易發(fā)生;而在背風面，隨攻角的增大引起轉捩所需的初始擾動幅值逐漸增大，說明轉捩越難以發(fā)生.

圖12 初始擾動幅值A0隨攻角的變化曲線Fig.12 Curves of amplitude A0 changing with attack angle

在工程應用中，判斷轉捩發(fā)生的預設值需要由實驗測得并需要進一步假定.若采用eN方法計算的擾動的N值達到預設值時，則認為轉捩發(fā)生了.根據(jù)圖13給出了當計算的擾動波的N值達到預設值發(fā)生轉捩的位置.可以看出，隨著攻角的增大迎風面轉捩位置向機翼前緣移動，背風面則與之相反.

圖13 轉捩位置ξ隨攻角的變化曲線Fig.13 Curves of transition position ξ changing with attack angle

本文的來流參數(shù)選取為高空10 km處，來流馬赫數(shù)為0.7，預計來流中擾動速度約為0.1m/s(約為來流速度的0.05%)，當擾動波的無量綱幅值達到來流的20%就會發(fā)生轉捩［19］，因此預計轉捩發(fā)生的N值為6左右.從圖13可以看出，此時轉捩發(fā)生的位置在0.1～0.2個弦長之間.實際發(fā)生的轉捩位置與來流擾動的幅值息息相關.

4 結論

在采用數(shù)值模擬方法得到無限展長后掠機翼基本流場的基礎上，通過求解O-S方程，研究了攻角對后掠機翼邊界層穩(wěn)定性的影響并進行了轉捩預測，得到結論如下:

1)隨著攻角的增大，流向速度在迎風面逐漸減小，而橫流速度線性增大，在背風面反之.

2)攻角對中性駐波的展向波數(shù)影響較小，但是駐波的中性曲線在流向的最大位置對攻角很敏感，在迎風面/背風面隨著攻角的增加呈線性增加/減小.

3)擾動波的增長在迎風面得到加強，在背風面受到抑制.

4)迎風面轉捩最容易發(fā)生，隨著攻角的增大，轉捩位置向機翼前緣移動.

5)當擾動速度為來流速度的0.05%時，轉捩發(fā)生的 N值在6左右，轉捩發(fā)生的位置在0.1 ～0.2 個弦長之間.

References)

［1］ Joslin R D.Overview of laminar flow control［M］.Virginia:National Aeronautics and Space Administration，Langley Research Center，1998.

［2］周恒，趙耕夫.流動穩(wěn)定性［M］.北京:國防工業(yè)出版社，2004:157-158.Zhou H，Zhao G F.Hydrodynamic stability［M］.Beijing:National Defense Industry Press，2004:157-158(in Chinese).

［3］吳永健.橫流不穩(wěn)定性實驗研究［D］.南京:南京航空航天大學，2002.Wu Y J.Experimental study on crossflow instabilities in the boundary-layer of swept wing［D］.Nanjing:Nanjing University of Aeronautics and Astronautics，2002(in Chinese).

［4］ Boltz FW，Kenyon G C，Allen C Q.Effects of sweep angle on the boundary-layer stability characteristics of an untapered wing at low speeds，Technical Note:D-338［R］.Moffett Field:National Aeronautics and Space Administration，1960.

［5］ Haynes T S.Nonlinear stability and saturation of crossflow vortices in swept-wing boundary layers［D］.Tempe:Arizona State University，1996.

［6］ Dagenhart JR，Saric W S.Crossflow stability and transition experiments in swept-wing flow［M］.Virginia:National Aeronautics and Space Administration，Langley Research Center，1999:1，7.

［7］ Reibert M S，Saric W S.Review of swept-wing transition，AIAA-1997-1816［R］.Reston:AIAA，1997.

［8］ Bushnell D M，Malik M R，HarveyW D.Transition prediction in external flows via linear stability theory［C］∥Symposium Transsonicum Ⅲ.Berlin Heidelberg:Springer-Verlag，1988:225.

［9］ Arnal D，Gasparian G，Salinas H.Recent advances in theoretical methods for laminar-turbulent transition prediction［J］.AIAA，1998:98-0223.

［10］孫朋朋，黃章峰.后掠角對后掠機翼邊界層穩(wěn)定性及轉捩的影響［J］.北京航空航天大學學報，2015，41(7):1313-1321.Sun P P，Huang Z F.Effect of the sweep angle on the stability and transition in a swept-wing boundary layer［J］.Journal of Beijing University of Aeronautics and Astronautics，2015，41(7):1313-1321(in Chinese).

［11］李鋒，汪翼云，崔爾杰.翼型大攻角繞流的數(shù)值模擬［J］.航空學報，1992，13(1):17-22.Li F，Wang Y Y，Cui E J.The numerical simulation of compressible flow around an airfoil at high angle of attack［J］.Acta Aeronautica et Astronautica Sinica，1992，13(1):17-22(in Chinese).

［12］吳鋆，王晉軍，李天.NACA0012翼型低雷諾數(shù)繞流的實驗研究［J］.實驗流體力學，2013，27(6):32-38.Wu J，Wang J J，Li T.Experimental investigation on low Reynolds number behavior of NACA0012 airfoil［J］.Journal of Experiments in Fluid Mechanics，2013，27(6):32-38(in Chinese).

［13］吳鋆，李天，王晉軍.低Reynolds數(shù)NACA0012翼型繞流的流動特性分析［J］.實驗力學，2014，29(3):265-272.Wu J，Li T，Wang J J.Characteristic analysis of flow around NACA0012 airfoil in a low-Reynold-number media［J］.Journal of Experimental Mechanics，2014，29(3):265-272(in Chinese).

［14］袁湘江，李國良，劉智勇，等.小攻角高超聲速鈍錐邊界層失穩(wěn)特性［J］.航空動力學報，2011，26(12):2805-2811.Yuan X J，LiG L，Liu ZY，etal.Study of the instability characteristic in the boundary layer of a hypersonic blunt cone at low angle of attack［J］.Journal of Aerospace Power，2011，26(12):2805-2811(in Chinese).

［15］王斌，白存儒，楊廣郡，等.后掠機翼低速流動轉捩位置的升華法測量［J］.實驗力學，2009，24(3):197-201.Wang B，Bai C R，Yang G J，et al.Measurement of transition location change of swept wing in a low speed flow based on sublimation method［J］.Journal of Experimental Mechanics，2009，24(3):197-201(in Chinese).

［16］孫朋朋.馬赫數(shù)、攻角及后掠角對后掠機翼邊界層穩(wěn)定性及轉捩的影響［D］.天津:天津大學，2015.Sun P P.Effect of Mach number，attack angle and sweep angle on the stability and transition in a swept-wing boundary layer［D］.Tianjin:Tianjin University，2015(in Chinese).

［17］黃章峰，逯學志，于高通.機翼邊界層的橫流穩(wěn)定性分析和轉捩預測［J］.空氣動力學學報，2014，32(1):14-20.Huang Z F，Lu X Z，Yu G T.Cross-flow instability analysis and transition prediction of airfoil boundary layer［J］.ACTA Aerodynamic Sinica，2014，32(1):14-20(in Chinese).

［18］韓步璋，黃奕裔，張其威，等.NACA0012翼型跨音速測壓實驗研究［J］.南京航空航天大學學報，1987，19(2):92-102.Han B Z，Hang Y Y，Zhang QW，et al.An experiment of pressure measurement for NACA0012 airfoil in a transonic wind tunnel［J］.Journal of Nanjing Aeronautical Institute，1987，19(2):92-102(in Chinese).

［19］ Huang Z F，Cao W，Zhou H.The mechanism of breakdown in laminar-turbulent transition of a supersonic boundary layer on a flat plate-temporal mode［J］.Science in China Series G:Mechanics and Astronomy，2005，48(5):614-625.