基于雙eN方法的短艙層流轉(zhuǎn)捩影響因素

孟曉軒，白俊強，張美紅，王美黎，何小龍，汪輝

1.西北工業(yè)大學 航空學院，西安 710072 2.西北工業(yè)大學 無人系統(tǒng)技術(shù)研究院，西安 710072 3.中國商用飛機有限責任公司 上海飛機設(shè)計研究院，上海 201210

現(xiàn)代民用飛機，尤其是民航客機，越來越追求經(jīng)濟性和環(huán)保性[1]，為提高大型民機的空氣動力學效率，大幅降低民機的油耗和污染排放，層流控制技術(shù)是一種極具減阻潛力和應(yīng)用前景的方法[2-4]。大涵道比渦扇發(fā)動機技術(shù)的快速發(fā)展帶來了發(fā)動機幾何尺寸的迅速增加，發(fā)動機阻力占全機阻力的比重也迅速增加，而短艙表面的摩擦阻力是短艙阻力的重要來源之一，因此對短艙外表面進行自然層流設(shè)計，能有效減小摩擦阻力，提升全機的經(jīng)濟性。此外發(fā)動機短艙對升力、力矩等氣動性能要求較低，約束較少；且不受附著線轉(zhuǎn)捩問題的影響，因此短艙的層流減阻設(shè)計具有較高的可行性[3]。經(jīng)自然層流(Natural Laminar Flow, NLF)設(shè)計的短艙稱之為自然層流短艙。

國外學者針對NLF短艙已開展過許多研究工作。美國蘭利研究中心[5]從20世紀80年代中期開始對發(fā)動機短艙的層流控制技術(shù)進行了可行性研究，經(jīng)分析認為層流控制技術(shù)有可能使短艙摩擦阻力減少約5%。Younghans和Lahti[6]在1984年利用實驗手段對單獨短艙和翼吊短艙進行了層流設(shè)計，在其外表面上均實現(xiàn)了大面積層流。德國宇航院(DLR)、羅爾斯-羅伊斯和MTU(Motorenund Turbinen-Union Friedrichshafen)[7]合作進行了自然層流短艙的飛行試驗，短艙表面的層流區(qū)達到了60%，且研究還表明噪聲和振動對層流控制的影響很小。Lin等[8]利用γ-Reθt轉(zhuǎn)捩模型預(yù)測了NLF短艙表面的轉(zhuǎn)捩現(xiàn)象，同試驗結(jié)果相比該方法還需進一步完善。目前，中國對自然層流的研究多集中在機翼上[9]，針對NLF短艙的相關(guān)研究較少，何小龍等[10]基于EFFD(Extended Free Form Deformation)參數(shù)化方法進行了NLF短艙和帶動力NLF短艙的優(yōu)化設(shè)計，取得了較好的結(jié)果。

發(fā)展高效可靠的轉(zhuǎn)捩預(yù)測方法有助于研究短艙的流動特點并進行層流減阻設(shè)計。常見的工程轉(zhuǎn)捩預(yù)測方法有基于湍流模式理論的γ-Reθt轉(zhuǎn)捩預(yù)測模式和基于線性穩(wěn)定性分析的eN轉(zhuǎn)捩預(yù)測方法。γ-Reθt轉(zhuǎn)捩預(yù)測模式是由Menter和Langtry[11-12]等提出的經(jīng)驗預(yù)測模型。經(jīng)典的γ-Reθt模型僅能考慮到沿流向的間歇因子變化規(guī)律，對流向T-S波不穩(wěn)定性引起的轉(zhuǎn)捩進行預(yù)測。為能對CF波引發(fā)的不穩(wěn)定性轉(zhuǎn)捩進行準確預(yù)測，Owen和Randall[13]于1952年率先提出Recf橫流轉(zhuǎn)捩判據(jù)，應(yīng)用起來比較復(fù)雜。Arnal等[14]在1984年提出了C1準則判據(jù)，應(yīng)用比較廣泛。之后德國宇航院的Grabe等[15-16]通過當?shù)鼗蠼鈾M流雷諾數(shù)，耦合C1準則進行橫流轉(zhuǎn)捩研究。中國學者[17-19]也先后展開了對橫流轉(zhuǎn)捩預(yù)測方法的研究。

中國對流動穩(wěn)定性的研究以天津大學周恒院士團隊為代表[20]，做了大量的基于穩(wěn)定性和直接數(shù)值模擬(DNS)的研究。Lee和Wu[21-22]在平板轉(zhuǎn)捩實驗中揭示了非線性邊界層不穩(wěn)定性和湍流發(fā)生的一種特征流動結(jié)構(gòu)，提出了一種由類孤立波主導的轉(zhuǎn)捩物理過程。李新亮等[23]做了大量的基于DNS的轉(zhuǎn)捩過程研究，然而其主要研究焦點是高超聲速下的轉(zhuǎn)捩過程機理，即面向航天工程需求。面向航空問題，基于求解線性穩(wěn)定性方程而發(fā)展的完全eN方法，由Smith[24]和Ingen[25]等首次提出，是目前航空界最被接受和廣泛認可的方法。eN方法的思想是對邊界層流動中的小擾動進行線性穩(wěn)定性分析，如果擾動逐漸衰減，則是穩(wěn)定的，如果擾動被放大，則是不穩(wěn)定的。對于不穩(wěn)定的擾動，從其開始放大處起，沿下游方向計算其累積線性放大倍數(shù)，當累積放大倍數(shù)到達擾動開始放大處振幅的eN倍時，認為轉(zhuǎn)捩發(fā)生。在轉(zhuǎn)捩點上所對應(yīng)不穩(wěn)定波的幅值放大因子的包絡(luò)線值為轉(zhuǎn)捩關(guān)鍵N因子(Ncrit)。eN方法對N的積分策略有很多種，如固定β法、包絡(luò)線法、雙eN方法等。雙eN轉(zhuǎn)捩預(yù)測方法針對T-S波和CF波選擇不同的擾動放大率積分策略，可以同時計算T-S波和橫流不穩(wěn)定性擾動誘導轉(zhuǎn)捩。在飛行試驗數(shù)據(jù)分析中發(fā)現(xiàn)，NTS-NCF雙因子策略的一大好處是，在沿無黏流線附近方向T-S波的激發(fā)過程很相近，幅值增長率都很相似，因此在馬赫數(shù)不大于1的情況下，沿著邊界層無黏流線方向?qū)因子的計算可以反應(yīng)T-S波增長的過程。而對低湍流背景的后掠機翼上的橫流轉(zhuǎn)捩，橫流駐波的初始振幅比橫流行波的振幅要大得多，因此雙eN策略對航空外流的分析誤差不大。朱震等[26]利用雙eN方法開展了翼身組合體復(fù)雜三維外形的轉(zhuǎn)捩預(yù)測研究，得到了較好的預(yù)測結(jié)果。另外，應(yīng)用雙eN策略計算飛行試驗結(jié)果時所取得的關(guān)鍵值(NTS,NCF)的分布都求解的是不可壓穩(wěn)定性方程。為和已有的工程經(jīng)驗對比，本文也只求解不可壓穩(wěn)定性方程。

工程上除直接求解線性穩(wěn)定性方程的eN方法外，還發(fā)展出來有eN-數(shù)據(jù)庫[27-28]方法。eN-數(shù)據(jù)庫方法在計算時不再求解穩(wěn)定性方程，而是直接根據(jù)邊界層形狀因子、速度型等信息在數(shù)據(jù)庫中查找平板邊界層流動的不穩(wěn)定擾動累積放大因子的包絡(luò)線，從而確定N值。eN-數(shù)據(jù)庫方法由于不需要求解穩(wěn)定性方程，具有計算速度快的特點，但只適用于機翼這類構(gòu)型的繞流流動，對于非機翼類構(gòu)型，如橢球、短艙以及機身等構(gòu)型，eN-數(shù)據(jù)庫方法并不適用。基于非線性穩(wěn)定性分析的方法，如將拋物化穩(wěn)定性方程(PSE)和CFD耦合求解的方法[29]，近年來已有工程化的程序出現(xiàn)，能夠比較方便地考慮邊界層內(nèi)非平行流和非線性的影響，較好地描述邊界層內(nèi)的流動結(jié)構(gòu)和擾動形態(tài)，但其控制方程為拋物型偏微分方程，不是當?shù)靥卣髦祮栴}，需要給定入口邊界條件沿流動方向向下游推進求解，使用不便，且存在穩(wěn)定性方程求解計算量大的缺點。

現(xiàn)已發(fā)展的諸多工程轉(zhuǎn)捩預(yù)測方法，無論是經(jīng)典的γ-Reθt轉(zhuǎn)捩預(yù)測模式或者考慮橫流不穩(wěn)定性轉(zhuǎn)捩而改進的γ-Reθt轉(zhuǎn)捩預(yù)測模式，還是eN-數(shù)據(jù)庫方法或者基于求解線性穩(wěn)定性方程的eN方法，多集中在對翼型、機翼或翼身組合體等構(gòu)型的轉(zhuǎn)捩預(yù)測上[30-31]，側(cè)重于機翼繞流的轉(zhuǎn)捩預(yù)測研究。針對短艙這類具有顯著三維效應(yīng)的類圓柱形構(gòu)型，在求解過程中，由于短艙不能同機翼一樣通過求解邊界層方程的方法獲得層流速度型，只能直接從CFD流場中提取邊界層速度型進行穩(wěn)定性分析，因此計算量較大。國內(nèi)外也鮮有相關(guān)文獻研究層流轉(zhuǎn)捩預(yù)測方法在短艙上的應(yīng)用，在這類構(gòu)型上有關(guān)層流轉(zhuǎn)捩的試驗數(shù)據(jù)儲備也不是很多。

由于NLF短艙在大型民機減阻設(shè)計上具有較高的可行性和可觀的發(fā)展前景，針對短艙表面的層流優(yōu)化設(shè)計就日趨重要，開展短艙層流轉(zhuǎn)捩的影響因素研究對NLF短艙的優(yōu)化設(shè)計具有較強的指導意義。本文基于線性穩(wěn)定性分析理論，從CFD流場中直接提取邊界層速度型信息進行穩(wěn)定性分析，將NTS-NCF雙因子轉(zhuǎn)捩判斷模塊同高精度雷諾平均Navier-Stokes(RANS)求解器[32]耦合，建立了一套可同時計算T-S波和橫流不穩(wěn)定性擾動誘導轉(zhuǎn)捩的流動轉(zhuǎn)捩自動判斷方法，對影響短艙層流轉(zhuǎn)捩的因素進行研究。

1 雙eN轉(zhuǎn)捩預(yù)測方法

本文所用雙eN轉(zhuǎn)捩預(yù)測方法的具體求解過程如圖1所示。首先設(shè)定初始轉(zhuǎn)捩位置，使其盡量靠近模型后緣從而得到足夠大的層流區(qū)域，利用RANS求解器計算流場，使用SA湍流模型模擬湍流流動。在CFD計算結(jié)果中提取邊界層信息，如速度、速度梯度、溫度、溫度梯度等，傳遞進入線性穩(wěn)定性分析求解器，分別計算幅值放大因子NTS和NCF的包絡(luò)線，并利用轉(zhuǎn)捩判據(jù)判斷轉(zhuǎn)捩位置，將計算得到的轉(zhuǎn)捩位置反饋到CFD求解器中，進行再一次的循環(huán)計算，直到循環(huán)得到的轉(zhuǎn)捩位置不變?yōu)橹埂?/p>

圖1 RANS求解器耦合eN穩(wěn)定性分析方法求解流程Fig.1 Solution for RANS solver coupling eN stability analysis method

1.1 邊界層信息提取

在邊界層內(nèi)分布超過80個網(wǎng)格點，模型前緣附近物面周向均布32個點。設(shè)定初始轉(zhuǎn)捩位置，將位于轉(zhuǎn)捩位置上游的模型表面和距表面0.01 m厚的區(qū)域都定義為層流，關(guān)掉湍流模型的產(chǎn)生項，進行CFD計算可得到層流區(qū)。沿著32個點依次向上游和下游尋找無黏流線：向上游尋找，直到駐點；向下游尋找，直到設(shè)定的轉(zhuǎn)捩位置或者流動分離點。將無黏流線投影到物面上，每隔一段距離記錄一個坐標，作為穩(wěn)定性分析的坐標點存儲在程序中。

沿指定點的法向方向，在每個點上提取邊界層信息，利用可壓縮伯努利方程計算邊界層頂部的速度Ue，取95%Ue處的壁面法向位置為邊界層頂部，記錄當?shù)氐乃俣仁噶浚瑢⒃擖c位置信息和法向方向與壁面交匯的點對應(yīng)。該方法對于壓力梯度不大或者負壓梯度區(qū)域的流動，近似效果較好。而對靠近前緣這些順壓梯度很大的位置，需要進行適當?shù)男拚涸搮^(qū)域的邊界層內(nèi)是有黏有旋流動，渦量和平均應(yīng)變率都很大，比值接近1，而在邊界層外是無黏無旋流動，渦量很小，平均應(yīng)變率較大，因此可以利用渦量和平均應(yīng)變率的比值是否為1來修正邊界層外緣的位置。

沿無黏流線記錄的坐標沿法線向物面投影，確定分離點、壓力最小值點：分離點用物面摩擦系數(shù)確定，表面極限流線方向和無黏流線方向的內(nèi)積為負的第1個點即為分離點；而壓力最小值點為沿無黏流線方向的壓力最小值點。將分離點、壓力最小值點依照計算的是T-S波主導的轉(zhuǎn)捩還是CF波主導的轉(zhuǎn)捩，作為計算相應(yīng)N值時積分的終止位置。具體計算時，若僅對CF波轉(zhuǎn)捩進行分析，可僅對幅值增長率積分至壓力最小值點。

在CFD流場中，根據(jù)無黏流線上記錄的坐標點和沿法向投影到壁面上的點，在兩點確定的區(qū)間上對流場信息進行插值，提取邊界層信息，包括速度、溫度、密度、壓強等。將速度矢量沿無黏流線方向投影得到流線方向的流向速度U和橫向方向的橫流速度W，進而得到速度梯度的二階空間導數(shù)，利用當?shù)氐臏囟取⒚芏鹊刃畔ⅲ梢杂嬎惝數(shù)乩字Z數(shù)。這些信息將會在求解線性穩(wěn)定性方程中使用。

1.2 線性穩(wěn)定性分析

穩(wěn)定性方程的分析方法有時間模式(Temporal Theory)和空間模式(Spatial Theory)兩種，分別從不同角度分析問題，利用群速度的關(guān)系可以導出相互變換關(guān)系。本文選擇用時間模式分析穩(wěn)定性方程，具有編程簡單、易于實現(xiàn)的優(yōu)點。

對于給定的擾動，指定其擾動頻率，隨著邊界層增厚，擾動幅度在中性點x0的下游開始增加，之后在下游某一位置xl下降，如圖2所示，圖中x為沿擾動波傳播方向的坐標，f為擾動波頻率，A為擾動波的幅值，N為擾動波的幅值放大因子。因此，從x0到xl處沿擾動波傳播方向的積分即為中性點總的幅值增長率，定義為

(1)

式中：αi為空間擾動波的幅值放大率。對于一系列含有不同頻率的擾動，可以得到總的幅值增長率的曲線系，其包絡(luò)線定義為

N=max[ln(A/A0)]f

(2)

式中：N為幅值放大因子，是一個比值，包含的是最容易被放大的擾動波的集合。

圖2 中性曲線(上)、幅值增長曲線與包絡(luò)線(下)Fig.2 Neutral curve (up), amplitude growth curves, and envelope lines (low)

首先，沿著每一條無黏流線方向，初步計算沿該方向的一條中性曲線，得到擾動的穩(wěn)定區(qū)域和不穩(wěn)定區(qū)域的分界線。通過中性曲線確定不穩(wěn)定擾動的頻率范圍、波數(shù)范圍和不同擾動各自開始被放大的位置，從而確定沿無黏流線方向所需考慮的點，即僅從中性曲線所包含的不穩(wěn)定區(qū)范圍內(nèi)對幅值放大率進行積分。在每個站位進行穩(wěn)定性分析時，對確定的波數(shù)范圍分成40個流向波數(shù)α′和橫向波數(shù)β′的組合，對每個點求解線性穩(wěn)定性方程，對這些頻率的擾動放大率從中性曲線下半支出發(fā)，沿著擾動傳播路徑(近似在勢流方向)進行積分，積分至下游的分離點或者壓力最小值點或者指定點，得到擾動累積放大因子N。

短艙上常見的轉(zhuǎn)捩樣式是自然轉(zhuǎn)捩，包括自由來流中的擾動穿過邊界層在邊界層里觸發(fā)的T-S波引起的T-S轉(zhuǎn)捩和因為橫流速度存在引起的橫流轉(zhuǎn)捩。在亞聲速時，T-S波的傳播方向和x方向夾角幾乎為0°，而橫流駐波的傳播方向和x方向夾角為85°～90°。因此本文采用固定波角法分別計算NTS和NCF，對于NTS，固定波角ψ=0°，即邊界層外邊界速度方向；對于NCF，采用同時固定其頻率和波角的方法計算CF波擾動累積放大因子，固定其頻率f=0 Hz，固定波角ψ=90°。

1.3 轉(zhuǎn)捩判斷

1.3.1 T-S轉(zhuǎn)捩判據(jù)

眾多試驗表明來流湍流度對轉(zhuǎn)捩有很大的影響，根據(jù)Dryden湍流度為0.1%、1%的試驗數(shù)據(jù)，Mark[33]總結(jié)出放大因子與湍流度相關(guān)聯(lián)的表達式，即Mark關(guān)系式：

N=-8.43-2.4ln(Tu)

(3)

如果0.1%

圖3為Mark關(guān)系式給出的轉(zhuǎn)捩雷諾數(shù)RexT與湍流度相關(guān)曲線，圖中，湍流度Tu采用對數(shù)坐標。可以發(fā)現(xiàn)，式(3)在特別低的湍流度下，不具有參考價值。而一般實際飛行條件下的湍流度都小于0.05%，此時，NTS.crit的值依賴于飛行試驗確定。

圖3 Mark關(guān)系式和不可壓平板轉(zhuǎn)捩預(yù)測試驗對比Fig.3 Comparison of Mark formula and incompressible plate transition prediction test

1.3.2 CF轉(zhuǎn)捩判據(jù)

對于CF波轉(zhuǎn)捩閾值NCF.crit，其關(guān)鍵值依賴于模型表面的加工狀況，即表面粗糙度和具體的飛行條件。在缺少相關(guān)試驗的條件下，參照風洞中NLF0415機翼橫流轉(zhuǎn)捩的標定研究，采用NCF.crit=6.4[35]這一固定值進行計算。

當不穩(wěn)定擾動在下游某點處累積的線性放大倍數(shù)達到擾動開始放大處振幅的eN,crit倍時，認為轉(zhuǎn)捩發(fā)生，該點即為對應(yīng)無黏流線的轉(zhuǎn)捩位置。將轉(zhuǎn)捩位置反饋給RANS求解器開始新一輪的循環(huán)，反復(fù)迭代直到轉(zhuǎn)捩位置收斂時終止，輸出轉(zhuǎn)捩位置。

2 算例驗證

2.1 數(shù)值方法準確性驗證

為驗證雙eN轉(zhuǎn)捩預(yù)測方法的可靠性和適應(yīng)性，下面針對6∶1標準橢球體(橢球體長矩L=2.4 m，短矩D=0.6 m)模型進行轉(zhuǎn)捩預(yù)測。該模型由德國宇航院進行了多個迎角和雷諾數(shù)的風洞試驗，試驗數(shù)據(jù)比較詳實和全面。CFD計算網(wǎng)格為800萬，附面層第1層網(wǎng)格為1.2×10-6m，法向增長率為1.1，物面網(wǎng)格的y+均小于1。迎風端對稱面的角度Φ=0°，背風端對稱面的角度Φ=180°。自由來流湍流度為0.1%，黏性比為10。根據(jù)經(jīng)驗關(guān)系式設(shè)定eN方法中的N值為：NTS=8，NCF=5.5。利用該方法對橢球的典型轉(zhuǎn)捩現(xiàn)象進行了數(shù)值模擬。圖4為標準橢球在Ma=0.136、α=15°、Tu=0.1%、Re=6.5×106狀態(tài)下試驗與eN方法計算所得摩阻系數(shù)Cf分布云圖的對比結(jié)果。

圖4 試驗和eN轉(zhuǎn)捩預(yù)測方法所得橢球轉(zhuǎn)捩位置對比Fig.4 Comparison of ellipsoid transition position obtained by test and eN transition prediction method

從計算結(jié)果可以看出，本文建立的雙eN穩(wěn)定性分析轉(zhuǎn)捩預(yù)測方法能夠較好地預(yù)測出橢球體邊界層內(nèi)出現(xiàn)的橫流不穩(wěn)定性轉(zhuǎn)捩現(xiàn)象，具有較高的預(yù)測精度，由此驗證了雙eN線性穩(wěn)定性分析轉(zhuǎn)捩預(yù)測方法對于類圓柱構(gòu)型表面的層流轉(zhuǎn)捩現(xiàn)象具有較強的模擬和預(yù)測能力。

2.2 網(wǎng)格收斂性驗證

為了對短艙層流轉(zhuǎn)捩的影響因素進行研究，本文選擇CRM模型的單獨短艙為研究對象，并對該短艙構(gòu)型進行網(wǎng)格收斂性驗證。本文選擇的CRM模型由NASA的亞聲速固定翼空氣動力技術(shù)研究小組和DPW的組織委員會合作設(shè)計開發(fā)，為CFD的驗證和確認工作提供基準外形[36]。

從CRM翼身組合體中分離出單獨短艙，通過平移旋轉(zhuǎn)變換使短艙軸線與x軸重合，以軸線前緣點為坐標原點。以CRM翼身組合體的參考面積、展長為短艙的參考面積、展長，以短艙的軸向長度為短艙參考弦長。短艙的參考面積S=191.85 m2，參考展長B=29.38 m，參考弦長c=5.966 m，以平均氣動弦1/4弦線位置(1.492 m, 0 m, 0 m)為短艙力矩參考點。

參考CRM翼身組合體設(shè)計工況，選擇短艙基準工況如下：高度H=11.277 6 km，自由來流馬赫數(shù)Ma=0.85，迎角α=2°，側(cè)滑角β=0°，基于短艙參考弦長的雷諾數(shù)Re=36.779×106，來流初始湍流度Tu=0.1%，黏性比Rt=10。強迫湍流度在上游遠場保持不變，在短艙前緣開始自由衰減。采用SA湍流模型模擬湍流流動，轉(zhuǎn)捩閾值根據(jù)經(jīng)驗取NTS.crit=8.1，NCF.crit=6.4。

選擇90萬、300萬、700萬、2 000萬和5 000萬5個不同密度的網(wǎng)格，計算短艙基準工況，進行網(wǎng)格收斂性驗證，各網(wǎng)格的基本信息如表1所示。所有網(wǎng)格壁面第1層網(wǎng)格高度均為1×10-6，滿足y+<1的條件。圖5給出了網(wǎng)格量為700萬的短艙表面網(wǎng)格。圖6為不同密度網(wǎng)格對短艙阻力系數(shù)CD的影響，圖7為不同密度網(wǎng)格計算的短艙頂部截面壓力系數(shù)Cp分布曲線。

從圖7所示不同密度網(wǎng)格下計算的壓力分布曲線來看，基準狀態(tài)下短艙外表面的主要流動特點是：前緣附近很快達到吸力峰值，下游完全處于逆壓梯度區(qū)，中間某處存在一個激波，在上表面激波位于45%短艙長度，在下表面激波位于15%短艙長度，左右兩側(cè)激波較弱，位于20%短艙長度。90萬網(wǎng)格(Tiny)沒有很好地模擬壓力系數(shù)沿流向的變化趨勢，也沒能準確捕捉到短艙上表面存在的激波。結(jié)合圖6中網(wǎng)格量對計算精度的影響，權(quán)衡計算精度和計算效率后，選擇700萬網(wǎng)格量(Medium)對短艙層流轉(zhuǎn)捩的影響因素進行分析。

表1 不同密度網(wǎng)格的基本信息Table 1 Basic information for different density meshes

圖5 短艙表面網(wǎng)格Fig.5 Surface mesh of nacelle

圖6 不同密度網(wǎng)格對短艙阻力系數(shù)的影響Fig.6 Effect of different density meshes on nacelle drag coefficient

圖7 不同密度網(wǎng)格所得短艙頂部截面壓力系數(shù)分布曲線Fig.7 Pressure coefficient distribution curves of nacelle top section obtained from different density meshes

3 短艙層流轉(zhuǎn)捩影響因素分析

3.1 變馬赫數(shù)分析

所研究的馬赫數(shù)變化范圍為0.6～0.87，覆蓋高亞聲速到亞聲速的工況。在馬赫數(shù)0.8～0.87之間研究了4個狀態(tài)點，詳細分析了巡航馬赫數(shù)附近短艙外表面轉(zhuǎn)捩位置對飛行速度的敏感性。

隨馬赫數(shù)變化的最顯著特征為:馬赫數(shù)升高，轉(zhuǎn)捩位置xtrans向下游移動，如圖8所示,圖中，V∞為來流速度。從壓力分布隨馬赫數(shù)變化的關(guān)系中可以看出，馬赫數(shù)增加引起短艙2個典型站位的壓力分布發(fā)生很大變化，且具有相同的發(fā)展趨勢，即隨著馬赫數(shù)增加，氣流壓縮性越來越強，負壓峰值難以維持很高的值，負壓峰值下降，下游的壓力梯度隨之減小，如圖9所示，激波被推遲但激波強度增強，圖中：壓力系數(shù)采用對數(shù)坐標。在圖10所示的壓力分布局部放大圖中，可以通過壓力分布隨位置變化的陡峭程度看出前緣附近壓力梯度的變化規(guī)律。由于T-S波對壓力梯度變化敏感，在飛行雷諾數(shù)下，越強的逆壓梯度引起T-S轉(zhuǎn)捩越早發(fā)生，因此得到轉(zhuǎn)捩位置隨馬赫數(shù)升高而向下游變化的規(guī)律。

圖8 CRM短艙轉(zhuǎn)捩位置隨馬赫數(shù)變化Fig.8 Variation of CRM nacelle transition position with Mach number

圖9 短艙前緣壓力梯度和轉(zhuǎn)捩位置隨馬赫數(shù)變化Fig.9 Variation of pressure gradient and transition position of naceller leading edge with Mach number

圖10 不同馬赫數(shù)下CRM短艙各截面壓力系數(shù)分布Fig.10 Pressure coefficient distribution of each crosssection of CRM nacelle with different Mach numbers

3.2 變雷諾數(shù)分析

雷諾數(shù)覆蓋了(5～40)×106的變化范圍。由于國內(nèi)目前沒有可靠的運行高雷諾數(shù)條件下靜聲風洞的條件和經(jīng)驗，因此進行了中等雷諾數(shù)到高雷諾數(shù)的計算與分析，研究轉(zhuǎn)捩位置對雷諾數(shù)的敏感性，對標風洞試驗可能出現(xiàn)的所有工況。

隨著雷諾數(shù)增加，轉(zhuǎn)捩位置有明顯向前移動的趨勢，如圖11所示。雷諾數(shù)增加到10×106以上，移動位置變得不明顯。而圖12的壓力系數(shù)分布對比結(jié)果表明，在相同馬赫數(shù)、迎角條件下，短艙上的壓力系數(shù)分布隨雷諾數(shù)變化很小。因此，在相同的壓力梯度作用下，雷諾數(shù)增加，層流更難以保持，帶來轉(zhuǎn)捩位置向上游移動。

圖11 CRM短艙轉(zhuǎn)捩位置隨雷諾數(shù)變化Fig.11 Variation of transition position of CRM nacelle with Reynolds number

3.3 變湍流度分析

考慮的自由來流湍流度變化范圍為0.5%～0.02%，包括0.5%、0.2%、0.1%、0.05%和0.02%。在巡航飛行工況下，對應(yīng)飛行高度上的來流湍流度很低，一般認為在0.05%以下。在起飛和著陸過程中，在不同的風況和飛機飛行速度下，來流湍流度變化范圍為0.05%～0.1%。對常規(guī)跨聲速風洞，一般湍流度在0.2%～0.5%，而對靜聲風洞，湍流度可以下降到0.05%以下。中國目前缺乏靜聲風洞，因此將湍流度敏感性分析擴展到試驗和飛行中所能遇到的所有工況，覆蓋大湍流度的飛行條件和背景噪聲非常低的安靜飛行條件。

湍流度表征自由來流中擾動的強度，擾動傳遞到邊界層產(chǎn)生不穩(wěn)定波，在低湍流度下，擾動線性增長，當擾動波幅值增長到一定程度時，非線性影響增強，發(fā)生二次失穩(wěn)或三波共振，很快發(fā)生轉(zhuǎn)捩，稱為自然轉(zhuǎn)捩。在高湍流度下，較強的初始擾動將跳過線性增長階段直接發(fā)生旁路轉(zhuǎn)捩(Bypass)；在本文研究的0.02%～0.5%湍流度范圍內(nèi)，擾動較弱，其轉(zhuǎn)捩途徑一般為自然轉(zhuǎn)捩，此時雙eN方法能夠反映擾動的線性增長過程，較好地進行轉(zhuǎn)捩預(yù)測。

如圖13所示，隨著湍流度增加，擾動變強，轉(zhuǎn)捩更容易發(fā)生，轉(zhuǎn)捩位置向上游移動，從短艙長度的10%位置向前移動到5%位置，變化范圍很小。流動的壓力系數(shù)分布不隨湍流度的變化而變化，如圖14所示。

圖13 CRM短艙轉(zhuǎn)捩位置隨湍流度變化Fig.13 Variation of transition position of CRM nacelle with degrees of turbulence

圖14 不同湍流度下CRM短艙頂部截面壓力系數(shù)分布Fig.14 Pressure coefficient distribution of top section of CRM nacelle with different degrees of turbulence

eN轉(zhuǎn)捩預(yù)測方法計算無黏流線方向N增長的包絡(luò)線，在前緣附近是一條以沿流向位置為函數(shù)的斜率很大的直線，在不同湍流度時，利用Mark關(guān)系式建立轉(zhuǎn)捩判據(jù)，因而能夠分辨出不同湍流度的影響。由于在很高的雷諾數(shù)下，N增長的包絡(luò)線是一條斜率很大的直線，在很大的湍流度變化范圍內(nèi)，轉(zhuǎn)捩位置的變化范圍在5%短艙長度以內(nèi)，表明在高雷諾數(shù)工況下，該短艙上的轉(zhuǎn)捩位置對湍流度的變化很不敏感。

3.4 變迎角分析

對于迎角來說，若偏離設(shè)計點過大，流動帶來的分離問題導致氣動力變化引起問題的嚴重性遠遠大于轉(zhuǎn)捩位置移動帶來的阻力變化。因此，僅考慮0°～6°的小迎角范圍。

隨著迎角的變化，轉(zhuǎn)捩位置在短艙的迎風面和背風面上出現(xiàn)較小變化，如圖15所示。迎角的變化引起了壓力分布的明顯變化，如圖16所示：在上表面，隨著迎角的增加，吸力峰值增加，壓力梯度緩慢減弱，激波增強。計算的轉(zhuǎn)捩位置在上表面向下游有小幅的移動，移動的范圍在5%短艙長度以內(nèi)，如圖17所示。主要原因是壓力梯度的變化幅度不大，所以作用效果不顯著。在下表面，隨著迎角的增加，外表面的壓力梯度增加，吸力峰值下降，和上表面變化趨勢相反，轉(zhuǎn)捩位置向前緣唇口移動，但是在6°迎角，逆壓梯度增加，逆壓梯度變化量的絕對值很小，不足以迅速激發(fā)T-S波發(fā)展至不穩(wěn)定態(tài)，引起在下表面的轉(zhuǎn)捩位置在0°～4°工況時隨迎角增加向前緣移動，在6°工況時向下游移動，表現(xiàn)出復(fù)雜的變化關(guān)系。

圖15 CRM短艙轉(zhuǎn)捩位置隨迎角變化Fig.15 Variation of transition position of CRM nacelle with angle of attack

總體而言，在小迎角變化范圍內(nèi)，短艙表面上的壓力分布出現(xiàn)較大的變化，逆壓梯度出現(xiàn)變化，eN轉(zhuǎn)捩預(yù)測方法對逆壓梯度的小幅度變化反應(yīng)敏感，但是轉(zhuǎn)捩位置的變化范圍仍在5%短艙長度范圍內(nèi)。

圖16 不同迎角下CRM短艙各截面壓力分布Fig.16 Pressure distribution of each cross-section of CRM nacelle with different angle of attack

圖17 短艙前緣壓力梯度和轉(zhuǎn)捩位置隨迎角變化Fig.17 Variation of pressure gradient and transition position of naceller leading edge with angle of attack

4 結(jié) 論

自由來流馬赫數(shù)、雷諾數(shù)、湍流度以及迎角均對所研究短艙的轉(zhuǎn)捩位置有一定影響，具體變化情況如下：

1) 隨著馬赫數(shù)增加，短艙外表面吸力峰值下降，激波位置向后移動，逆壓梯度減小，轉(zhuǎn)捩位置后移。隨著雷諾數(shù)和自由來流湍流度增加，短艙各截面的壓力分布幾乎不變，轉(zhuǎn)捩位置更靠近唇口。

2) 隨著迎角增加，短艙外表面迎風區(qū)吸力峰值下降，壓力梯度增加，但壓強變化量變小，轉(zhuǎn)捩位置向上游移動(壓強變化量特別小的工況除外)。背風區(qū)吸力峰值增加，壓力梯度減小，轉(zhuǎn)捩位置向下游移動。

3) 在所研究的工況下，馬赫數(shù)和迎角帶來的壓力梯度變化明顯，引起轉(zhuǎn)捩位置發(fā)生變化；而雷諾數(shù)和湍流度在高雷諾數(shù)工況下影響相對較小，轉(zhuǎn)捩位置的移動不超過5%短艙長度。因此在設(shè)計階段，對保持層流設(shè)計來說，要盡量避免較大的逆壓梯度，保持順壓梯度。