三種典型翼型邊界層穩(wěn)定性對比分析

王 菲，額日其太，王 強，蘇沛然

(北京航空航天大學(xué)能源與動力工程學(xué)院，北京 100191)

0 引言

摩擦阻力可占亞音速運輸機總阻力的50%，而層流的阻力比相同雷諾數(shù)湍流的阻力要小90%［1］，因此減小表面摩擦阻力可以有效地降低飛行器的總阻力［2］。為了計算機翼的摩擦阻力從而進一步更合理地預(yù)測機翼飛行特性，轉(zhuǎn)捩位置的預(yù)測成為了重要前提工作。對于高空飛行外界擾動較小的特點，線性穩(wěn)定性方法的eN方法仍然是目前工程應(yīng)用中最常用的轉(zhuǎn)捩預(yù)測方法。L.J.Runyan和 D.George-Falvy通過修正 Mack的穩(wěn)定性方法計算了翼型在不同工況下的轉(zhuǎn)捩位置并與實驗對照，得到相應(yīng)的最大增長率的值［3］，為之后的相關(guān)研究提供了豐富的數(shù)據(jù)。朱萬琳等人在二維不可壓流動的穩(wěn)定性數(shù)值分析基礎(chǔ)上進行了轉(zhuǎn)捩的判定［4］。為了便于工程上翼型的選取及優(yōu)化，針對多個典型翼型邊界層穩(wěn)定性的對比分析具有重要的參考意義。首先本文采用萘升華實驗方法驗證對稱翼型SD8020數(shù)值模擬的準確性，并擬定穩(wěn)定性計算的幅值放大因子的值，然后針對常規(guī)翼型(NACA0012)，自然層流翼型(NACA64-204)和超臨界翼型(RAE2822)這三種典型翼型采用線性理論方法計算在不同弦長雷諾數(shù)情況下轉(zhuǎn)捩位置的變化，研究流向壓力梯度對轉(zhuǎn)捩點雷諾數(shù)、轉(zhuǎn)捩點弦長百分數(shù)和其相應(yīng)的擾動頻率的影響，并分析影響這三種翼型穩(wěn)定性增長的原因。

1 數(shù)值模擬方法

1.1 邊界層平均流求解

采用無粘流和邊界層方法求解基本流，利用F-S法變換基本流方程［5］，取，這里 ue表示自由流速度。相對于流函數(shù)ψ(x，y)的無量綱流函數(shù)為:

由此可得連續(xù)方程和動量方程為:

邊界條件為 η =0，f=f'=0;η =ηe，f'=1。其中導(dǎo)數(shù)表示對η求導(dǎo)，m表示無量綱梯度:

1.2 線性穩(wěn)定性方程

擾動的發(fā)展模式可以分成時間模式和空間模式，時間模式的擾動在空間某點僅隨時間規(guī)律變化，空間模式的擾動傳播則沿空間某一方向變化，而后者是可以通過實驗結(jié)果驗證的［6］，所以本文采用更接近實驗條件的空間模式方法進行擾動增長計算。將流函數(shù)定義為:

這里φ表示流函數(shù)的型函數(shù)。α和β分別表示流向和展向的擾動波數(shù)，ω表示擾動的圓頻率。對應(yīng)于空間模式 α ，β 為復(fù)數(shù)，α = αr+iαi，β = βr+iβi;ω 為實數(shù)。這里定義無量綱擾動相速度c:

利用線性穩(wěn)定性方程(O-S方程)計算擾動波在流場中的傳播，由于在二維基本流條件下三維擾動的失穩(wěn)點大于二維擾動的失穩(wěn)點，因此引入的二維擾動也是合理的［7］，即假定β=0。這里二維流無量綱化的O-S方程可以表示為一個四階常微分方程［1］:

邊界條件為 y=0，φ =0，φ'=0;y=δ，φ =0，φ'=0。R 表示弦長雷諾數(shù)，這里c表示弦長。

轉(zhuǎn)捩預(yù)測采用工程上最為常用的eN方法預(yù)測翼型上表面的轉(zhuǎn)捩位置。擾動增長率N為:

其中－αi為給定頻率下的空間放大率。

2 實驗驗證

本實驗在北京航空航天大學(xué)流體力學(xué)研究所低速風(fēng)洞中完成。實驗利用萘的易升華特性，將萘溶于丙酮溶液并噴涂于機翼模型表面。由于轉(zhuǎn)捩點前后氣流對壁面的剪切力變化很大，使萘在湍流區(qū)內(nèi)升華速度加快，因此模型表面湍流區(qū)的白色萘涂層先消失，從而在層流和湍流之間形成了明顯的邊界。選擇實驗狀態(tài)為Rec=7.8 ×105，Ma=0.087，攻角為 －2°以保持較長的層流區(qū)，該翼型轉(zhuǎn)捩位置在78%弦長左右(圖1)。

圖1 萘升華實驗Fig.1 Naphthalene sublimation experiment

對應(yīng)上述實驗條件，采用eN方法對SD8020翼型進行轉(zhuǎn)捩預(yù)測。圖2為翼型上表面擾動增長曲線，計算20條不同頻率的擾動增長曲線，選取常用的擾動放大因子N=9。當(dāng)某一頻率增長率最先達到N=9時認為計算得到轉(zhuǎn)捩點，用虛線與頻率增長率曲線的交點對應(yīng)的弦長位置來表示。計算得到轉(zhuǎn)捩位置在79%左右，與實驗結(jié)果吻合，因此以下方法均采用e9方法進行轉(zhuǎn)捩預(yù)測。

圖2 翼型SD8020在弦長雷諾數(shù)Rec=7.8×105下不同頻率N因子曲線及轉(zhuǎn)捩位置Fig.2 N factor curves and transition locations of SD8020 airfoils for Rec=7.8 ×105

3 計算結(jié)果及分析

利用上述eN方法，選用三種典型翼型NACA0012，NACA64-204，RAE2822(圖3)，計算來流馬赫數(shù) Ma=0.3，攻角為0的條件下這三種翼型壓力分布及轉(zhuǎn)捩位置。圖4為三種翼型上表面壓力系數(shù)Cp分布曲線，其中NACA64-204和RAE2822分別對應(yīng)自然層流翼型和超臨界翼型，這兩種翼型壓力分布特點是都具有較長的順壓梯度區(qū)，在前緣附近劇烈的順壓梯度之后，有較長一段弱順壓區(qū)，最后達到最低壓力點。NACA64-204和RAE2822最低壓力點分別在弦長的40%和55%位置處。

圖5(a)(b)(c)為三種翼型上表面在相同弦長雷諾數(shù)(Rec=9×106)下擾動增長曲線，采用e9方法進行轉(zhuǎn)捩預(yù)測。三個翼型中NACA0012的轉(zhuǎn)捩點雷諾數(shù)最小，也即在相同弦長雷諾數(shù)下最早發(fā)生轉(zhuǎn)捩。由于在二維流動中T-S波對轉(zhuǎn)捩的影響起主導(dǎo)作用，而順壓梯度又對T-S波有穩(wěn)定作用。對應(yīng)于圖4壓力系數(shù)分布，NACA0012僅有10%弦長的順壓區(qū)，遠小于后兩者，因此其擾動也越不穩(wěn)定。三種翼型的高頻擾動增長起始位置總是比低頻更靠近前緣，而受順壓梯度影響，高頻增長速度也更容易減慢甚至衰減。相對的低頻擾動增長點起始位置更靠后，而其受順壓梯度影響的區(qū)域也減少。若低頻擾動增長率在進入逆壓區(qū)時并未開始衰減，則逆壓將會導(dǎo)致擾動迅速增長。這種現(xiàn)象與中性曲線的分布規(guī)律是一致的。正是由于逆壓梯度使速度剖面產(chǎn)生拐點，有拐點的速度剖面的穩(wěn)定性界限比無拐點的剖面低得多(拐點準則)［7］。

圖6(a)、圖5(b)和圖6(b)表示 NACA64-204翼型分別在弦長雷諾數(shù)為6×106、9×106和12×106不同頻率擾動增長率。由這三幅圖可以看出隨弦長雷諾數(shù)增長，在層流區(qū)內(nèi)由于順壓區(qū)與逆壓區(qū)的長度比變大，其轉(zhuǎn)捩點雷諾數(shù)也變大。當(dāng)弦長雷諾數(shù)增大到12×106時，在順壓區(qū)擾動增長率已經(jīng)達到了N=9，因此可以判定在最小壓力點之前轉(zhuǎn)捩已經(jīng)發(fā)生。NACA64-204翼型擾動增長的特點是擾動通過翼型前緣后部大范圍弱順壓梯度后，逆壓梯度對其影響非常明顯，即使很小的逆壓梯度也會使低頻擾動迅速增長。

弦長雷諾數(shù)的增長使高頻擾動增長率變大，因此使最不穩(wěn)定頻率(對應(yīng)轉(zhuǎn)捩位置的擾動頻率)也相應(yīng)提高(圖7)，而不同翼型相同弦長雷諾數(shù)下對應(yīng)最不穩(wěn)定頻率也存在很大差別。

圖5 三種翼型在弦長雷諾數(shù)Rec=9×106下不同頻率N因子曲線及轉(zhuǎn)捩位置Fig.5 N factor curves and transition locations of three airfoils for Rec=9×106

如圖8所示，同一翼型，隨弦長雷諾數(shù)增長，中性曲線位置向雷諾數(shù)更大頻率更高的方向偏移，中性曲線的不穩(wěn)定區(qū)域也變大。相比高弦長雷諾數(shù)，低弦長雷諾數(shù)下相同頻率的擾動更早進入中性曲線的不穩(wěn)定區(qū)域，即擾動增長的更早，且沿頻率增長方向擾動不穩(wěn)定區(qū)范圍更小。因此隨弦長雷諾數(shù)的增長，同一翼型最不穩(wěn)定擾動波頻率也變大。

表1統(tǒng)計了三種翼型在弦長雷諾數(shù)分別為6×106、9×106、12×106下的轉(zhuǎn)捩點雷諾數(shù)和對應(yīng)弦長百分數(shù)。從表中可以看出，工程上常采用雷諾平均方法計算全湍流流場［8］或固定轉(zhuǎn)捩點雷諾數(shù)來假設(shè)轉(zhuǎn)捩位置的方法是不準確的。同一翼型隨著弦長雷諾數(shù)的增長，轉(zhuǎn)捩點雷諾數(shù)也隨之變大，但并未按線性比例增長，其對應(yīng)的弦長百分數(shù)反而減小，也就是說轉(zhuǎn)捩點相對提前。由于轉(zhuǎn)捩位置的變化會對翼型性能計算產(chǎn)生很大的影響［9］，可見準確計算翼型的轉(zhuǎn)捩點并分別求解層流區(qū)和湍流區(qū)是十分必要的。

圖8 NACA64-204不同弦長雷諾數(shù)下的中性曲線Fig.8 The neutral curves for different Reynolds numbers for NACA64-204

4 結(jié)論

本文首先采用萘升華實驗方法驗證對稱翼型SD8020數(shù)值模擬的準確性，并擬定eN方法的幅值放大因子的值，然后采用該數(shù)值方法分別計算了在不同的弦長雷諾數(shù)下三種典型翼型的壓力分布，擾動增長率以及最不穩(wěn)定擾動波頻率，計算了翼型上表面轉(zhuǎn)捩位置，并分析了不同翼型壓力分布對擾動發(fā)展和轉(zhuǎn)捩位置的影響。結(jié)果表明:(1)在相應(yīng)實驗條件下采用幅值放大因子N=9的轉(zhuǎn)捩預(yù)測位置與實驗結(jié)果吻合。(2)在三種典型翼型轉(zhuǎn)捩預(yù)測結(jié)果中，相同弦長雷諾數(shù)下，NACA0012最先發(fā)生轉(zhuǎn)捩，而NACA64-204和RAE2822都保持著較長的層流區(qū)。(3)壓力梯度對擾動增長有很大影響，順壓梯度對二維翼型擾動有明顯的穩(wěn)定作用。(4)同一翼型隨弦長雷諾數(shù)增長，順壓區(qū)內(nèi)擾動增長更快，因此轉(zhuǎn)捩點相對提前;又因為順壓區(qū)占弦長百分比固定，因此層流區(qū)內(nèi)順壓區(qū)與逆壓區(qū)的長度比變大。又因為中性曲線變寬并向下游方向移動，最不穩(wěn)定頻率也隨之增大。由于轉(zhuǎn)捩點雷諾數(shù)并不是工程上常用的固定值，而轉(zhuǎn)捩位置的變化會對翼型性能參數(shù)計算產(chǎn)生很大的影響，因此準確計算翼型的轉(zhuǎn)捩點，可以得到更合理的翼型氣動性能參數(shù)(如摩擦阻力)，也為工程上選取和優(yōu)化翼型提供了基本的參考數(shù)據(jù)。

表1 三種翼型不同弦長雷諾數(shù)下轉(zhuǎn)捩點雷諾數(shù)及轉(zhuǎn)捩點弦長百分數(shù)Table 1 Reynolds numbers and the percentage of chord length at transition points for three a irfoils

［1］周恒，趙耕夫.流動穩(wěn)定性［M］.北京:國防工業(yè)出版社，2004:226.

［2］額日其太，沈遐齡，劉火星，等.飛機層流控制技術(shù)及其在大型運輸機上的應(yīng)用［A］.大型飛機關(guān)鍵技術(shù)高層論壇暨中國航空學(xué)會2007年學(xué)術(shù)年會［C］.中國廣東深圳:2007.

［3］RUNYAN L J，GEORGE-FALVY D.Amplification factors at transition on an unswept wing in free flight and on a swept wing in wind tunnel［A］.17th Aerospace sciences meeting［C］.New Orleans，Louisiana:AIAA，1979.

［4］朱萬琳.二維層流轉(zhuǎn)捩計算及應(yīng)用［D］.西北工業(yè)大學(xué)，2002.

［5］陳懋章編著.粘性流體動力學(xué)基礎(chǔ)［M］.北京:高等教育出版社，2002.

［6］CRIMINALE W O，JCAKSON T L，JOSLIN R D.Theory and computation of hydrodynamic stability［M］.Cambridge University Press，2003.

［7］SCHLICHTING H.Boundary-layer theory［M］.7th ed.MCGraw-Hill，1979.

［8］DRIVER J，ZINGG D W.Numerical aerodynamic optimization incorporating laminar-turbulent transition prediction［J］.AIAA Journal.2007，45(8):1810 －1818.

［9］KRUMBEIN A.Automatic transition prediction and application to 3d wing configurations［A］.44th AIAA Aerospace sciences meeting and exhibit［C］.Reno，Nevada:2006.