王 菲，額日其太，王 強(qiáng)，蘇沛然

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

0 引言

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

1 數(shù)值模擬方法

1.1 邊界層平均流求解

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

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

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

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

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

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

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

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

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

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

2 實(shí)驗(yàn)驗(yàn)證

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

圖1 萘升華實(shí)驗(yàn)Fig.1 Naphthalene sublimation experiment

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

圖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 計(jì)算結(jié)果及分析

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

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

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

弦長雷諾數(shù)的增長使高頻擾動(dòng)增長率變大，因此使最不穩(wěn)定頻率(對(duì)應(yīng)轉(zhuǎn)捩位置的擾動(dòng)頻率)也相應(yīng)提高(圖7)，而不同翼型相同弦長雷諾數(shù)下對(duì)應(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ù)下相同頻率的擾動(dòng)更早進(jìn)入中性曲線的不穩(wěn)定區(qū)域，即擾動(dòng)增長的更早，且沿頻率增長方向擾動(dòng)不穩(wěn)定區(qū)范圍更小。因此隨弦長雷諾數(shù)的增長，同一翼型最不穩(wěn)定擾動(dòng)波頻率也變大。

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

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

4 結(jié)論

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

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

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

［2］額日其太，沈遐齡，劉火星，等.飛機(jī)層流控制技術(shù)及其在大型運(yùn)輸機(jī)上的應(yīng)用［A］.大型飛機(jī)關(guān)鍵技術(shù)高層論壇暨中國航空學(xué)會(huì)2007年學(xué)術(shù)年會(huì)［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)捩計(jì)算及應(yīng)用［D］.西北工業(yè)大學(xué)，2002.

［5］陳懋章編著.粘性流體動(dòng)力學(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.