短艙邊界層的穩(wěn)定性分析及轉(zhuǎn)捩預(yù)測

牛萌浩，文 豪，韓宇峰，蘇彩虹，*，孟曉軒

(1.天津大學(xué)機(jī)械工程學(xué)院，天津 300072；2.西北工業(yè)大學(xué)航空學(xué)院，西安 710072)

0 引言

準(zhǔn)確預(yù)測邊界層的轉(zhuǎn)捩位置是飛行器設(shè)計(jì)中需要考慮的重要問題。由于湍流邊界層表面摩阻顯著高于層流邊界層[1]，準(zhǔn)確預(yù)測邊界層的轉(zhuǎn)捩位置才能準(zhǔn)確計(jì)算飛行器的表面阻力。一般來說，民用渦扇發(fā)動(dòng)機(jī)短艙阻力占飛機(jī)總阻力的5%左右[2]，隨著發(fā)動(dòng)機(jī)涵道比不斷增大，短艙直徑也隨之增大，短艙表面阻力在全機(jī)阻力中所占的比重也明顯增加。所以在短艙設(shè)計(jì)時(shí)，希望盡可能優(yōu)化其外罩型面，使其表面邊界層最大程度的層流化，以提升經(jīng)濟(jì)性[3-4]。因此，準(zhǔn)確預(yù)測短艙表面邊界層的轉(zhuǎn)捩位置對(duì)于短艙的優(yōu)化設(shè)計(jì)具有重要意義。

目前，可用于工程上的轉(zhuǎn)捩預(yù)測方法可主要分為以下三類：（1）基于邊界層當(dāng)?shù)貐?shù)關(guān)系式的預(yù)測方法，例如廣泛采用的 Reθ/Me方法[5]；（2）以RANS為基礎(chǔ)的轉(zhuǎn)捩模型，如Langtry和Menter等[6-8]提出的γ-Reθt轉(zhuǎn)捩模型；（3）以線性穩(wěn)定性理論（linear stability theory,LST）為基礎(chǔ)的eN方法[9-11]。前兩種方法的最大優(yōu)點(diǎn)是，可以非常方便地植入到CFD求解器中，計(jì)算效率高，因而應(yīng)用廣泛。然而，其主要的不足是缺乏對(duì)流動(dòng)物理機(jī)制的考慮，需要通過大量實(shí)驗(yàn)數(shù)據(jù)校正模型中所包含的眾多參數(shù)[12]。第三種方法基于線性穩(wěn)定性理論，通過計(jì)算邊界層中不穩(wěn)定波在線性放大過程中幅值增長的倍數(shù)來預(yù)測轉(zhuǎn)捩。由于考慮了轉(zhuǎn)捩過程的基本機(jī)理，即邊界層的不穩(wěn)定性，eN方法被普遍認(rèn)為是更準(zhǔn)確、更可靠的轉(zhuǎn)捩預(yù)測方法[12-17]。但是，由于該方法并不基于局部變量，計(jì)算量比前兩種方法更大，而且對(duì)于復(fù)雜外形，在與多塊網(wǎng)格并行計(jì)算的CFD求解器結(jié)合時(shí)存在一定的困難。

eN方法在與復(fù)雜CFD流場求解器結(jié)合時(shí)主要存在兩個(gè)方面的困難[18]。一是，需要CFD求解器提供垂直于物面方向的物理量以確定線性穩(wěn)定性控制方程（即Orr-Sommerfeld方程）的系數(shù)。而一般情況下，求解流場時(shí)很難保證法向網(wǎng)格垂直于物面。二是，在應(yīng)用eN方法時(shí)，需要跟隨擾動(dòng)波的傳播方向積分增長率。同樣，網(wǎng)格線不可能剛好與擾動(dòng)波傳播方向重合。因此，其應(yīng)用要比僅依賴于當(dāng)?shù)刈兞康霓D(zhuǎn)捩模式方法復(fù)雜得多。但鑒于轉(zhuǎn)捩預(yù)測在短艙型面設(shè)計(jì)中的重要價(jià)值，為了獲取準(zhǔn)確、可靠的轉(zhuǎn)捩預(yù)測結(jié)果，有些情況下在計(jì)算資源上的投資是值得的。

大部分關(guān)于短艙轉(zhuǎn)捩的研究都是從自然層流化設(shè)計(jì)的角度開展，例如文獻(xiàn)[19-22]，直接關(guān)注短艙轉(zhuǎn)捩預(yù)測的工作非常少，有關(guān)層流短艙的實(shí)驗(yàn)結(jié)果也很少。Lin等[23]采用γ-Reθ模型對(duì)短艙表面邊界層的轉(zhuǎn)捩位置進(jìn)行預(yù)測（馬赫數(shù)范圍為0.8～0.88，攻角從-2°～0°），發(fā)現(xiàn)短艙周向大部分區(qū)域與實(shí)驗(yàn)結(jié)果的比較是合理的。杜璽等[2]在德國ETW風(fēng)洞開展了層流短艙實(shí)驗(yàn)，采用TSP方法測量轉(zhuǎn)捩。實(shí)驗(yàn)觀察到短艙表面存在明顯的楔形層流區(qū)，猜測是風(fēng)洞中灰塵雜質(zhì)帶來的影響。盡管無法明確給出轉(zhuǎn)捩線沿短艙周向的分布，但總體來看，馬赫數(shù)從0.8增加到0.87，短艙層流區(qū)面積明顯增加；攻角從0°增加到4°，短艙下部（迎風(fēng)區(qū)附近）層流區(qū)增加。孟曉軒等[24]分析了短艙轉(zhuǎn)捩的影響因素，但缺乏對(duì)短艙邊界層穩(wěn)定性的細(xì)致分析。

本文基于線性穩(wěn)定性理論細(xì)致分析短艙邊界層的穩(wěn)定性特征，給出短艙邊界層中起主導(dǎo)作用的失穩(wěn)機(jī)制，并采用eN方法給出轉(zhuǎn)捩線的分布，分析攻角、巡航馬赫數(shù)對(duì)轉(zhuǎn)捩的影響規(guī)律。目的是為短艙優(yōu)化設(shè)計(jì)及開展實(shí)驗(yàn)測試提供重要參考。

1 計(jì)算模型和流動(dòng)條件

短艙表面流動(dòng)狀態(tài)主要與外罩型線有關(guān)，因此，可將短艙的進(jìn)排氣系統(tǒng)簡化，且忽略安裝構(gòu)件的影響[25]，得到的通氣短艙如圖1所示。

圖1 短艙示意圖（U∞位于x-z平面內(nèi)）Fig.1 Schematic diagram of the nacelle (U∞liesin the x-z plane)

短艙的軸向長度為5.5 m，坐標(biāo)原點(diǎn)位于短艙進(jìn)氣道的中心，x軸沿短艙軸線方向。標(biāo)記短艙頂部為θ =0°（背風(fēng)面），底部為θ = 180°（迎風(fēng)面）。考慮巡航狀態(tài)下的流動(dòng)條件，飛行高度為11887 m，來流溫度為216.65 K，壁面溫度條件為絕熱壁，飛行速度為250.8 m/s，以自由來流物理量定義的單位雷諾數(shù)為5.58×106/m，側(cè)滑角為0°，來流攻角為5°，馬赫數(shù)為0.85。在本文第4節(jié)中，我們還分析了攻角分別為1°和3°，馬赫數(shù)分別為0.8和0.89的情況。

2 基于線性穩(wěn)定性理論的e N方法

高空中的自然轉(zhuǎn)捩是由邊界層中不穩(wěn)定波的持續(xù)增長和放大導(dǎo)致的。以線性穩(wěn)定性理論為基礎(chǔ)的eN方法本質(zhì)上是基于一種模態(tài)競爭的思想，即所有可能的小擾動(dòng)各自獨(dú)立增長，最先達(dá)到預(yù)設(shè)轉(zhuǎn)捩閾值的擾動(dòng)就是觸發(fā)轉(zhuǎn)捩的關(guān)鍵擾動(dòng)，而其對(duì)應(yīng)的位置就是預(yù)測的轉(zhuǎn)捩位置。

根據(jù)線性穩(wěn)定性理論，小擾動(dòng)可以寫成如下的行進(jìn)波形式：

將式（1）代入線性化的N-S方程中，減去平均流滿足的方程，得到所謂的Orr-Sommerfeld方程（O-S方程）。對(duì)于二維情況，β =0。已知物理量沿法向的分布，給定擾動(dòng)頻率ω，求解O-S方程的特征值問題可以得到不穩(wěn)定波的波數(shù)αr和放大率-αi（αi＜0表示擾動(dòng)波是增長的）以及特征函數(shù)。從擾動(dòng)開始放大的位置即中性點(diǎn)開始，沿?cái)_動(dòng)傳播方向（二維情況中擾動(dòng)向下游傳播），對(duì)不穩(wěn)定波的放大率積分，可得擾動(dòng)波的幅值為：

式中，x0為中性點(diǎn)的位置，A0為擾動(dòng)在中性點(diǎn)處的初始幅值。

對(duì)于一般的三維邊界層，擾動(dòng)往往是三維的，采用O-S方程求解β（空間模式中β為復(fù)數(shù)）時(shí)，需要補(bǔ)充額外的條件。由于α、β滿足一定的相容關(guān)系[26]，可取βi=0。Cebeci和Stewartson[27]通過研究一個(gè)波包的演化，采用最速下降法提出一個(gè)補(bǔ)充條件：

eN方法需要沿著擾動(dòng)傳播方向?qū)υ鲩L率積分，計(jì)算幅值的放大倍數(shù)。擾動(dòng)沿群速度Vg的方向傳播：

這需要計(jì)算時(shí)對(duì)每個(gè)頻率在每個(gè)積分站位求解群速度方向，再求解O-S方程得到增長率，計(jì)算量非常大。一般認(rèn)為，對(duì)于以流向流動(dòng)為主的情況，擾動(dòng)群速度方向與邊界層外緣的勢流方向很接近[28-30]。為了節(jié)省計(jì)算量，我們沿著無黏流線計(jì)算不穩(wěn)定波的幅值增長。這樣，N值由下式給出：

s為無黏流線上的積分站位。假設(shè)轉(zhuǎn)捩判據(jù)為N = Ncri，不同頻率的擾動(dòng)達(dá)到轉(zhuǎn)捩閾值Ncri的位置不同，其中最靠上游的就是預(yù)測的轉(zhuǎn)捩位置。Ncri的具體取值需要依賴實(shí)驗(yàn)數(shù)據(jù)和經(jīng)驗(yàn)確定。

穩(wěn)定性方程的求解采用天津大學(xué)自主開發(fā)的SAYR程序，該程序采用四點(diǎn)四階Malik格式離散O-S方程，采用Muller法求邊界層不穩(wěn)定模態(tài)的特征值。目前，該程序已成功應(yīng)用于從亞聲速到高超聲速多種構(gòu)型邊界層的穩(wěn)定性分析中，如后掠翼型[31-32]、超聲速平板[33]、高超聲速圓錐[34-36]以及組合體外形[37]等。

3 短艙表面邊界層的轉(zhuǎn)捩預(yù)測

3.1 計(jì)算基本流場

在短艙四周和內(nèi)部構(gòu)建結(jié)構(gòu)網(wǎng)格，計(jì)算域入口和出口距短艙約為26～27個(gè)短艙長度。短艙表面的邊界層內(nèi)布置70個(gè)左右的網(wǎng)格點(diǎn)，總網(wǎng)格量約為5300萬。借助CFD開源程序，利用加速收斂的多重網(wǎng)格技術(shù)，采用二階精度的有限差分法求解RANS方程，得到短艙的基本流場。求解RANS方程時(shí)，采用Spalart-Allmaras（S-A）湍流模型，取固定轉(zhuǎn)捩位置為x =3.7 m，此位置之前的流場為層流解，之后采用S-A湍流模式求解。

邊界層中T-S波的穩(wěn)定性受壓力梯度的影響。通常情況下，順壓梯度使邊界層更穩(wěn)定；反之，逆壓梯度使邊界層更不穩(wěn)定。圖2給出了4個(gè)周向位置處壓力系數(shù)沿x的變化。可以看出，在背風(fēng)面θ = 0°處，只有很短的順壓區(qū)，逆壓區(qū)長度較長，在x =1.6 m附近出現(xiàn)激波。在迎風(fēng)面θ =180°處，順壓區(qū)相對(duì)較長，在x = 2.4 m處出現(xiàn)激波。側(cè)面的情況介于兩者之間。由于逆壓梯度的存在，短艙表面會(huì)出現(xiàn)邊界層分離。在分離區(qū)附近，邊界層局部平行性假設(shè)不再成立，線性穩(wěn)定性理論不再適用。因此，我們關(guān)注的主要是短艙前半部分，即邊界層分離前的流場。

圖2 壓力系數(shù)沿x的分布Fig.2 Variation of the pressure coefficient along x

為了驗(yàn)證流場計(jì)算的可靠性，特別是從不同流向位置開始采用S-A模型計(jì)算對(duì)短艙上游流場的影響，我們將固定轉(zhuǎn)捩位置選在短艙末端，即全場用層流求解，這時(shí)計(jì)算全部流場需要較長的時(shí)間才能收斂。選取激波前的位置，比較層流解和固定轉(zhuǎn)捩位置的層流解加S-A湍流模型所得速度剖面，結(jié)果如圖3所示。縱坐標(biāo)η表示垂直于壁面的方向。速度用飛行速度進(jìn)行無量綱化?？梢钥闯?，在激波出現(xiàn)前，邊界層剖面是一致的。這表明，在短艙前部，即邊界層分離前，固定轉(zhuǎn)捩位置的湍流模式可以給出定常收斂的基本流場，可作為線性穩(wěn)定性分析的基礎(chǔ)。之后的其他工況均采用該方法進(jìn)行基本流場計(jì)算。

為了驗(yàn)證邊界層內(nèi)的網(wǎng)格量是足夠的，我們將邊界層內(nèi)的網(wǎng)格從70個(gè)增加到130個(gè)。圖4給出了網(wǎng)格加密前后θ =0°和θ =180°子午面上不同流向位置速度剖面的對(duì)比?？梢钥闯觯Y(jié)果吻合很好。這說明，邊界層內(nèi)的網(wǎng)格量是足夠的。在之后其他工況的計(jì)算中，也采用相同的網(wǎng)格布置。

3.2 獲取線性穩(wěn)定性分析所需的邊界層剖面

采用LST進(jìn)行分析時(shí)，需要提供沿壁面法向方向的物理量及導(dǎo)數(shù)，以得到O-S方程的系數(shù)矩陣。但是，通常CFD的計(jì)算網(wǎng)格并不垂直于壁面，為此，需在垂直于壁面方向上重構(gòu)網(wǎng)格，采用插值的方法得到新網(wǎng)格下的流場信息。

以圖5（a）所示的壁面某一位置P(i,k)為例，說明從該點(diǎn)出發(fā)提取流場信息的步驟，具體方法可參考文獻(xiàn)[38]。i、k分別表示貼體坐標(biāo)系下的流向和展向指標(biāo)，網(wǎng)格線單位矢量分別記為ai、ak。

圖3 不同計(jì)算方法所得流向速度剖面的比較Fig.3 Comparison of streamwise velocity profiles obtained by different methods

（1）求出P點(diǎn)沿壁面法向的單位矢量 aj=ak×ai；

（2）沿著aj方向，以指數(shù)型分布布置新網(wǎng)格點(diǎn)，并確保邊界層內(nèi)有足夠的點(diǎn)；

（3）采用四點(diǎn)三階Lagrange插值方法，得到新網(wǎng)格點(diǎn)上的各物理量；

（4）確定邊界層外緣處的速度方向，即無黏流線的方向。從P點(diǎn)開始，沿流線方向在壁面網(wǎng)格上插值出下一個(gè)積分站位。重復(fù)上述步驟，依次向下游推進(jìn)，最終可得到由無黏流線和壁面法向確定的曲面內(nèi)的流場，即圖5（b）所示新網(wǎng)格上的流場。

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

短艙邊界層中存在的不穩(wěn)定模態(tài)有T-S波和橫流渦，需要分析在轉(zhuǎn)捩中起主導(dǎo)作用的模態(tài)。由于橫流失穩(wěn)主要發(fā)生在短艙側(cè)面，我們?cè)诙膛搨?cè)面選取三條典型的流線，分別是從唇口處θ =50°、θ = 90°和θ =120°出發(fā)的流線，分別記為Line 1、Line 2和Line 3，如圖6所示。在每條流線上選取一個(gè)離唇部不遠(yuǎn)的 位 置：Q1(0.501 m,-1.303 m,1.180 m)、Q2(0.506 m,-1.751 m,0.112 m)和Q3(0.506 m,-1.564 m,-0.793 m)。將Q1、Q2、Q3處的速度分別投影到無黏流線方向和垂直于流線方向上，得到勢流方向的速度分量us和橫流速度分量uc。圖7（a、b）分別給出了us和uc沿法向的分布?？梢园l(fā)現(xiàn)，橫流速度分量很小，不到邊界層外緣速度的3%。

圖4 加密前后流向速度剖面的比較（加密的工況只顯示了部分網(wǎng)格點(diǎn)）Fig.4 Comparison of streamwise velocity profiles before and after the mesh refinement (only a part of the grid nodes is shown for the denser grid case)

圖5 用于線性穩(wěn)定性分析的新網(wǎng)格Fig.5 New gridsused for the linear stability analysis

圖6 短艙側(cè)面流線Fig.6 Streamlineson the side of the nacelle

圖7 速度剖面Fig.7 Velocity profiles

采用線性穩(wěn)定性理論計(jì)算三條流線上的T-S波和定常橫流渦（stationary vortex）的增長率。圖8（a～c）分別給出了三條流線上不同頻率的T-S波和定常橫流渦（頻率為0）的增長率沿流向的變化?？梢园l(fā)現(xiàn)，T-S波開始失穩(wěn)的位置更靠前，且增長率顯著大于定常橫流渦。利用式（5），從中性點(diǎn)開始沿著無黏流線積分增長率得到N值曲線，結(jié)果如圖9所示。可見，T-S波的N值顯著大于定常橫流渦的N值。因此，對(duì)于巡航狀態(tài)下的小攻角情況，T-S波比橫流的不穩(wěn)定性更強(qiáng)，是轉(zhuǎn)捩的主導(dǎo)模態(tài)。在之后的轉(zhuǎn)捩預(yù)測中，我們只考慮T-S波失穩(wěn)引起的轉(zhuǎn)捩。

圖8 定常橫流渦與T-S波增長率的比較Fig.8 Comparison of growth ratesof the stationary crossflow vortex and the T-Swave

圖9 定常橫流渦和T-S波積分所得N值的比較Fig.9 Comparison of N factorsof the stationary crossflowvortex and the T-Swave

3.4 結(jié)合判據(jù)預(yù)測轉(zhuǎn)捩位置

前面提到，eN方法中轉(zhuǎn)捩判據(jù)的確定需要依賴實(shí)驗(yàn)和經(jīng)驗(yàn)。但對(duì)于寬體客機(jī)的發(fā)動(dòng)機(jī)短艙，雷諾數(shù)高達(dá)3×107左右，可獲得的實(shí)驗(yàn)數(shù)據(jù)很少。并且，由于構(gòu)型和流動(dòng)條件與航空翼型顯著不同，也無法直接借用翼型常用的轉(zhuǎn)捩判據(jù)。例如，Mack[39]曾提出一個(gè)與來流湍流度有關(guān)的經(jīng)驗(yàn)的轉(zhuǎn)捩判據(jù)，即：

式中，Tu表示來流湍流度。Lin等曾將其用于自然層流化短艙的轉(zhuǎn)捩預(yù)測中[23]。然而，由于目前并不清楚實(shí)驗(yàn)中湍流度的具體情況，無法給出可靠的Ncri。因此，本文暫時(shí)不討論具體的轉(zhuǎn)捩位置，僅給出不同判據(jù)下轉(zhuǎn)捩線的分布。待獲得進(jìn)一步的實(shí)驗(yàn)數(shù)據(jù)后，可標(biāo)定Ncri值，從而確定最終的轉(zhuǎn)捩位置。

以從θ =90°處出發(fā)的流線為例，不同頻率T-S波的N值曲線在圖9（b）中給出。暫時(shí)將轉(zhuǎn)捩閾值Ncri定為5～10的整數(shù)，可以得到不同頻率T-S波達(dá)到轉(zhuǎn)捩閾值對(duì)應(yīng)的位置，如圖10所示。對(duì)于某一Ncri，最靠上游的位置就是預(yù)測的轉(zhuǎn)捩位置，其對(duì)應(yīng)的頻率就是該判據(jù)下觸發(fā)轉(zhuǎn)捩的關(guān)鍵頻率?？梢姡S著Ncri增大，轉(zhuǎn)捩位置靠后，關(guān)鍵頻率降低。Ncri從5增加到10時(shí)，轉(zhuǎn)捩位置從0.57 m推遲至0.99 m，關(guān)鍵頻率從4.5 kHz降到2.9 kHz。

圖10 不同頻率T-S波對(duì)應(yīng)的轉(zhuǎn)捩位置Fig.10 Transition location for T-S waveswith different frequencies

至此，對(duì)于某一轉(zhuǎn)捩閾值Ncri，我們得到了一個(gè)轉(zhuǎn)捩點(diǎn)。由于無黏流線并不完全在θ =90°子午面內(nèi)，因此得到的轉(zhuǎn)捩點(diǎn)也會(huì)偏離該子午面。對(duì)短艙的其他周向位置都進(jìn)行類似的分析，可以得到轉(zhuǎn)捩線沿周向角度的分布，結(jié)果如圖11（a）所示。圖11（b）給出了轉(zhuǎn)捩線在短艙上的分布?？梢?，在θ =0°～90°范圍內(nèi)，轉(zhuǎn)捩位置變化不大，而在θ =90°～180°，即靠近迎風(fēng)面處，轉(zhuǎn)捩位置隨周向角度增加明顯后移。迎風(fēng)面處轉(zhuǎn)捩位置最靠后，而背風(fēng)面轉(zhuǎn)捩靠前。這與通常認(rèn)識(shí)的小攻角對(duì)轉(zhuǎn)捩的影響結(jié)果一致[23]。

圖11 短艙轉(zhuǎn)捩預(yù)測結(jié)果Fig.11 Transition prediction result for the nacelle

從圖11中注意到，當(dāng)Ncri比較大時(shí)，轉(zhuǎn)捩位置在迎風(fēng)面附近存在空檔。為了具體考察迎風(fēng)面處的流動(dòng)情況，圖12（a）給出了θ =180°子午面上的流線分布?？梢姡吔鐚釉趚 =2.4 m處開始分離，壁面附近出現(xiàn)分離泡。圖12（b）給出了對(duì)應(yīng)的壁面摩擦系數(shù)沿x方向的分布。Cf曲線與x軸的第一個(gè)交點(diǎn)代表分離開始的位置，第二個(gè)交點(diǎn)代表流動(dòng)再附的位置，圖13給出了不同頻率T-S波的N值曲線。前面曾提到，在邊界層分離附近，線性穩(wěn)定性理論不再適用，因此，N值曲線在分離位置附近停止。從圖中可以看出，轉(zhuǎn)捩閾值Ncri大于5時(shí)，T-S波無法達(dá)到給定的轉(zhuǎn)捩閾值。因此，圖11（a、b）所示的轉(zhuǎn)捩線在迎風(fēng)面存在空檔。

圖12 迎風(fēng)面處的流線和壁面摩擦系數(shù)沿x的分布Fig.12 Streamlines and skin friction distributions along x in the windward side

圖13 迎風(fēng)面不同頻率T-S波的N值曲線Fig. 13 N factors of T-Swaves with different frequencies in the windward side

對(duì)于層流分離泡，通常認(rèn)為Cf到達(dá)最小值的位置可作為轉(zhuǎn)捩位置[40]，即圖12（b）所示的x =2.8 m處。在高雷諾數(shù)下，流動(dòng)分離會(huì)促使轉(zhuǎn)捩快速發(fā)生，對(duì)于當(dāng)前工況，也可近似地把分離發(fā)生的位置當(dāng)作轉(zhuǎn)捩位置。在本文之后的分析中，轉(zhuǎn)捩線存在空檔意味著eN方法無法給出轉(zhuǎn)捩位置，實(shí)際情況中的轉(zhuǎn)捩由邊界層分離引起。

4 關(guān)鍵參數(shù)對(duì)短艙轉(zhuǎn)捩的影響

4.1 攻角

對(duì)于巡航馬赫數(shù)為0.85的情況，我們還計(jì)算了攻角分別為1°、3°的情況，并與前面分析的5°攻角的結(jié)果進(jìn)行比較。

圖14（a、b）分別給出了攻角為1°和3°轉(zhuǎn)捩線在短艙上的分布情況。與圖11（b）比較可見，隨著攻角減小，短艙背風(fēng)面和迎風(fēng)面轉(zhuǎn)捩位置的差別減小，除靠近迎風(fēng)面的區(qū)域外，轉(zhuǎn)捩線整體更靠近下游。此外，比較3°攻角和5°攻角的情況，可以發(fā)現(xiàn)，3°攻角下，轉(zhuǎn)捩線更加分散。這意味著，轉(zhuǎn)捩判據(jù)不確定性引起的轉(zhuǎn)捩位置的誤差更大。

圖15（a、b）分別給出了Ncri= 5和Ncri= 8時(shí)不同攻角下轉(zhuǎn)捩位置隨周向角度的變化?？梢姡?dāng)Ncri= 5時(shí)，在靠近迎風(fēng)面的θ = 150°～180°區(qū)域內(nèi)，攻角越小，轉(zhuǎn)捩位置越靠前；而在此之外的背風(fēng)面和側(cè)面區(qū)域，攻角減小使得轉(zhuǎn)捩位置更加靠后?？偟膩砜?，攻角的存在使得迎風(fēng)面附近的區(qū)域更加穩(wěn)定，而背風(fēng)面及側(cè)面的區(qū)域更加不穩(wěn)定。

圖14 不同攻角的轉(zhuǎn)捩位置Fig.14 Transition locationsunder different anglesof attack

圖15 不同攻角轉(zhuǎn)捩位置的比較Fig.15 Comparison of transition locationsunder different anglesof attack

當(dāng)Ncri= 8時(shí)，對(duì)比3°攻角和5°攻角的結(jié)果可以發(fā)現(xiàn)，攻角減小，背風(fēng)面及側(cè)面大部分區(qū)域轉(zhuǎn)捩靠后。當(dāng)然，還有一種可能是，需要根據(jù)實(shí)際情況調(diào)整轉(zhuǎn)捩判據(jù)Ncri的取值。但Ncri取值不同并不影響定性的規(guī)律。對(duì)于1°攻角的工況，在周向θ =90°～180°范圍內(nèi)，在Ncri≥6條件下eN方法無法給出轉(zhuǎn)捩發(fā)生的位置，轉(zhuǎn)捩位置由邊界層分離決定。圖16（a、b）分別給出了攻角為1°和3°情況下，迎風(fēng)面θ =180°子午面上的流線分布，可以看到下游均存在流動(dòng)分離。

圖16 迎風(fēng)面處的流線分布Fig.16 Streamlinesin the windward side

為了具體說明不同攻角下T-S波的失穩(wěn)情況，圖17（a～c）分別給出了周向θ =45°處不同頻率T-S波增長率的等值線圖。可以發(fā)現(xiàn)，隨著攻角減小，失穩(wěn)區(qū)范圍和增長率都明顯減小。因此，減小攻角抑制了不穩(wěn)定波的增長，流場變得穩(wěn)定，轉(zhuǎn)捩推遲發(fā)生。前面曾提到，文獻(xiàn)[2]通過實(shí)驗(yàn)觀察到，攻角從0°增加到4°，短艙下部層流區(qū)明顯增加，與我們的結(jié)果是一致的。

4.2 馬赫數(shù)

在來流攻角為5°時(shí)，我們還分析了馬赫數(shù)分別為0.80和0.89的工況，并與之前分析的馬赫數(shù)為0.85的工況進(jìn)行比較，以考察巡航馬赫數(shù)對(duì)轉(zhuǎn)捩的影響。

圖18（a、b）分別給出了馬赫數(shù)為0.80、0.89的轉(zhuǎn)捩預(yù)測結(jié)果。對(duì)比圖11（b）可以發(fā)現(xiàn)，馬赫數(shù)減小，轉(zhuǎn)捩線整體更靠上游。并且，對(duì)于馬赫數(shù)為0.8的情況，不同判據(jù)下的轉(zhuǎn)捩線分布更為集中。圖19（a、b）分別給出了Ncri= 5和Ncri= 8時(shí)不同馬赫數(shù)下轉(zhuǎn)捩位置隨周向角度的變化?？梢?，馬赫數(shù)增加使得所有子午面上的轉(zhuǎn)捩位置均向下游移動(dòng)。

圖17 不同攻角的增長率等值線圖Fig.17 Growth rate contoursunder different anglesof attack

從圖18（b）中可以看出，當(dāng)馬赫數(shù)為0.89時(shí)，eN方法無法給出迎風(fēng)面附近區(qū)域的轉(zhuǎn)捩位置。此時(shí)，轉(zhuǎn)捩位置與邊界層分離有關(guān)。圖20給出了馬赫數(shù)為0.89工況中θ =180°子午面上的流線分布，可以看出下游存在流動(dòng)分離。還有一種可能是，轉(zhuǎn)捩仍然是T-S波觸發(fā)的，而我們選用的轉(zhuǎn)捩判據(jù)閾值Ncri太大。這就需要根據(jù)實(shí)驗(yàn)測試的具體情況調(diào)整轉(zhuǎn)捩判據(jù)的閾值。但是，不管判據(jù)的選擇為何，馬赫數(shù)對(duì)轉(zhuǎn)捩位置影響的規(guī)律不會(huì)改變。

圖18 不同馬赫數(shù)的轉(zhuǎn)捩位置Fig.18 Transition locationsunder different Mach numbers

圖20 迎風(fēng)面處的流線Fig.20 Streamlines in the windward side

圖21 不同馬赫數(shù)的增長率等值線圖Fig.21 Growth rate contours under different Mach numbers

類似地，以周向θ =90°為例，圖21（a～c）分別給出了不同馬赫數(shù)下T-S波增長率的等值線圖。可見，隨著馬赫數(shù)的增大，T-S波開始失穩(wěn)的位置更靠近下游，失穩(wěn)區(qū)逐漸縮小，最大增長率也逐漸降低，流場更穩(wěn)定，因此導(dǎo)致轉(zhuǎn)捩推遲發(fā)生。前面曾提到，杜璽等的實(shí)驗(yàn)結(jié)果表明[2]，馬赫數(shù)從0.8增加到0.87，短艙表面層流區(qū)明顯增加，這與我們的結(jié)論是一致的。

5 結(jié)論

本文以民航發(fā)動(dòng)機(jī)翼吊短艙為研究對(duì)象，采用CFD求解器獲得基本流場，對(duì)短艙邊界層進(jìn)行了穩(wěn)定性分析，并采用基于線性穩(wěn)定性理論的eN方法進(jìn)行了轉(zhuǎn)捩預(yù)測；給出了不同轉(zhuǎn)捩判據(jù)下短艙轉(zhuǎn)捩線的分布，研究了關(guān)鍵參數(shù)如攻角、巡航馬赫數(shù)對(duì)轉(zhuǎn)捩位置的影響規(guī)律。所得結(jié)論如下：

1）巡航狀態(tài)下，短艙攻角很小，橫流速度分量很小，橫流不穩(wěn)定性對(duì)轉(zhuǎn)捩的影響不大，起主導(dǎo)作用的不穩(wěn)定波是T-S波；

2）轉(zhuǎn)捩起始位置對(duì)來流攻角和飛行馬赫數(shù)非常敏感。針對(duì)本文研究的攻角范圍，結(jié)果表明，攻角增加，迎風(fēng)面附近的轉(zhuǎn)捩顯著推遲，而背風(fēng)面及側(cè)面大部分區(qū)域的轉(zhuǎn)捩提前。在本文討論的巡航馬赫數(shù)范圍內(nèi)，隨著飛行馬赫數(shù)增大，短艙表面邊界層更加穩(wěn)定，轉(zhuǎn)捩位置整體后移。

致謝：感謝西北工業(yè)大學(xué)汪輝副教授、白俊強(qiáng)教授、631所顏洪教授、清華大學(xué)張宇飛副教授在流場計(jì)算和分析中給予的極大支持和幫助；感謝中國商飛上海飛機(jī)設(shè)計(jì)院張美紅研究員、劉凱禮研究員對(duì)本文工作給予的寶貴建議。此項(xiàng)工作是在國家超級(jí)計(jì)算天津中心的“天河一號(hào)”超級(jí)計(jì)算機(jī)上完成，感謝“天河一號(hào)”的大力支持。