跨聲速自然層流標(biāo)模CRM-NLF 風(fēng)洞試驗快速轉(zhuǎn)捩預(yù)測

劉銀輝，吳東潤，杜 璽，林大楷

(1.中國商用飛機(jī)有限責(zé)任公司 北京民用飛機(jī)技術(shù)研究中心，北京 102211；2.民用飛機(jī)設(shè)計數(shù)字仿真技術(shù)北京市重點實驗室，北京 102211)

0 引言

層流設(shè)計能給飛機(jī)帶來十分顯著的減阻收益[1]。研究表明，如果將層流設(shè)計成功應(yīng)用于機(jī)翼、發(fā)動機(jī)短艙、水平尾翼和垂直尾翼等部件，具有減少15%整機(jī)阻力的潛力，其中機(jī)翼采用層流設(shè)計帶來的收益最顯著，具有減少10%整機(jī)阻力的潛力。根據(jù)Arcara等[2]的估算，如果民機(jī)能在主翼上翼面、垂尾和平尾上實現(xiàn)50%弦長的層流覆蓋，在發(fā)動機(jī)短艙上實現(xiàn)40%弦長的層流覆蓋，巡航升阻比可提升14.7%，從而可以降低對應(yīng)方案起飛總重量的9.9%、操作空重的5.7%。現(xiàn)役大型民用飛機(jī)中，波音787 的短艙[3]以及波音737 Max[4]的翼梢小翼采用了層流設(shè)計，如圖1 所示。這標(biāo)志著先進(jìn)的飛機(jī)制造技術(shù)已經(jīng)可以一定程度上滿足層流設(shè)計的要求。開展層流設(shè)計研究，逐步推進(jìn)飛機(jī)各部件的層流應(yīng)用，是保證大型民用飛機(jī)具有氣動性能優(yōu)勢的重要手段之一。

圖1 層流飛機(jī)[4]Fig.1 Laminar aircraft[4]

邊界層轉(zhuǎn)捩預(yù)測是CFD 最具挑戰(zhàn)的難題之一。高效可靠的轉(zhuǎn)捩預(yù)測分析是飛機(jī)層流設(shè)計的關(guān)鍵。工程上常用的轉(zhuǎn)捩預(yù)測方法分為兩類：基于RANS 的轉(zhuǎn)捩模型方法和eN方法。對于第一類方法，隨著對轉(zhuǎn)捩物理機(jī)理認(rèn)知的逐漸深入和試驗數(shù)據(jù)的不斷豐富，轉(zhuǎn)捩模型已取得豐碩的研究成果。2021 年第一屆AIAA RANS 轉(zhuǎn)捩模型和預(yù)測會議中，針對標(biāo)模ERCOFTAC T3 平板、NLF(1)-0416 翼型、6∶1 橢球體和CRM-NLF 機(jī)翼，十余家單位采用不同的轉(zhuǎn)捩方法進(jìn)行了轉(zhuǎn)捩預(yù)測，但預(yù)測的轉(zhuǎn)捩位置偏差較大[5]。目前，文獻(xiàn)公開發(fā)表的基于RANS 方程的轉(zhuǎn)捩模型已多達(dá)幾十種。其中比較著名的有：基于間歇因子的轉(zhuǎn)捩模型方法，如 γ-Reθ[6-8]等；基于物理機(jī)理的層流動能轉(zhuǎn)捩模型，如kT-kL-ω[9-10]、k-ω-γ[11-13]等。總之，基于間歇因子的轉(zhuǎn)捩模型高度依賴于試驗數(shù)據(jù)的經(jīng)驗擬合，對轉(zhuǎn)捩前區(qū)域的物理機(jī)制考慮不足，適用范圍有限；基于層流動能的轉(zhuǎn)捩模型能較為充分地考慮轉(zhuǎn)捩過程的物理機(jī)理，然而，現(xiàn)有基于層流動能的轉(zhuǎn)捩模型難以解析地給出三維邊界層流動的轉(zhuǎn)捩起始位置和轉(zhuǎn)捩過程[14]。這些轉(zhuǎn)捩模型有待進(jìn)一步改進(jìn)，從而準(zhǔn)確預(yù)測工程上三維復(fù)雜外形的邊界層轉(zhuǎn)捩。

第二類是eN方法。eN方法一般是對數(shù)值計算（如直接數(shù)值模擬、帶轉(zhuǎn)捩的RANS 等）得到的層流基本流進(jìn)行穩(wěn)定性分析獲得N值包絡(luò)線。基于線性穩(wěn)定性理論的半經(jīng)驗eN方法仍然是飛機(jī)設(shè)計計算中得到最廣泛認(rèn)可的轉(zhuǎn)捩預(yù)測方法[15]，它通過擾動幅值放大因子N判斷轉(zhuǎn)捩的發(fā)生。相比二維和軸對稱邊界層，針對復(fù)雜三維外形開展的研究相對較少。基于eN方法進(jìn)行三維邊界層轉(zhuǎn)捩分析的難點在于展向波參數(shù)和積分路徑的選取等，比較常用的方法包含時間模式包絡(luò)法、鞍點法、固定展向波數(shù)法和雙eN法等。時間模式包絡(luò)法的N值采用時間模式沿群速度方向積分獲得，如Malik 開發(fā)的轉(zhuǎn)捩軟件COSAL[16]。鞍點法最早由Cebeci 和Stewartson[17]提出，擾動沿群速度方向傳播，通過固定頻率對不同展向波數(shù)下的最大增長率進(jìn)行積分計算N值，然后再計算頻率的N值包絡(luò)，相對保守。天津大學(xué)團(tuán)隊在超聲速/高超聲速邊界層中開展了簡化鞍點法、基于內(nèi)波理論提出射線跟蹤法、一般三維邊界層廣義增長率守恒關(guān)系等[18-22]方面的工作，并開發(fā)了對應(yīng)的轉(zhuǎn)捩預(yù)測軟件。固定展向波數(shù)法，適用于無限展長后掠平板和后掠翼型等，擾動展向波數(shù)不變，得出不同展向波數(shù)和頻率下的N值包絡(luò)線。Malik 等[23]發(fā)現(xiàn)，假定無黏勢流方向的展向波數(shù)為實數(shù)，擾動增長率沿?zé)o黏勢流方向進(jìn)行N值計算，計算量小且精度足夠。基于Malik 等[23]的發(fā)現(xiàn)，可以進(jìn)一步分別對Tollmien-Schlichting（T-S）模態(tài)和橫流（CF）模態(tài)進(jìn)行N值計算，這種方法被稱作雙eN方法。T-S 模態(tài)無黏勢流方向和群速度方向夾角一般不超過5°，T-S 模態(tài)采用固定波角法（無黏勢流方向的展向波數(shù)為0）或無黏勢流方向的展向波數(shù)取較小的實數(shù)進(jìn)行N值計算，而CF 模態(tài)采用固定展向波數(shù)法進(jìn)行N值計算，如NASA 的eMalik[24]和LASTRAC[25]，以及西北工業(yè)大學(xué)的韓忠華等[26]發(fā)展的翼身組合體轉(zhuǎn)捩自動判斷方法等。為了提高計算效率，后續(xù)文獻(xiàn)基于eN方法提出了簡化的eN數(shù)據(jù)庫轉(zhuǎn)捩預(yù)測方法[27-31]。本文采用eN方法進(jìn)行跨聲速轉(zhuǎn)捩預(yù)測分析，簡化的eN數(shù)據(jù)庫方法不再贅述。

eN方法需CFD 求解器提供層流基本流計算的邊界層信息，一般地，基于RANS 的CFD 工具普遍存在計算周期過長、一些情況下解的精度或可靠性不足、CFD 用戶技能水平差異大[32]等問題。因此，中國商飛北研中心開發(fā)了面向工程應(yīng)用的快速CFD 分析工具CFAST，基于試驗轉(zhuǎn)捩位置設(shè)置CFD 轉(zhuǎn)捩點，轉(zhuǎn)捩點前的流場為層流基本流，單工況的計算時間僅為200 s，大幅縮短計算周期。本文研究對象是NASA的跨聲速高雷諾數(shù)自然層流標(biāo)模CRM-NLF[33]。首先，基本流場由CFAST 得出；然后，采用基于線性穩(wěn)定性理論的eN方法進(jìn)行轉(zhuǎn)捩預(yù)測分析，為層流機(jī)翼設(shè)計和優(yōu)化提供理論依據(jù)和參考。下面對研究模型、快速CFD 分析工具CFAST 和eN方法進(jìn)行簡介，并基于標(biāo)模NLF(2)-0415 對計算結(jié)果進(jìn)行驗證。

1 研究模型和數(shù)值方法

1.1 標(biāo)模CRM-NLF 簡介

自然層流標(biāo)模CRM-NLF 是美國NASA 近年來采用跨聲速自然層流超臨界機(jī)翼設(shè)計的翼身組合體標(biāo)模，進(jìn)行了變迎角、變雷諾數(shù)的高雷諾數(shù)風(fēng)洞試驗驗證，得到超臨界自然層流機(jī)翼的邊界層轉(zhuǎn)捩特性，并公開了其幾何模型、CFD 計算結(jié)果、試驗數(shù)據(jù)等[33-34]。標(biāo)模CRM-NLF 轉(zhuǎn)捩由橫流駐波模態(tài)和T-S 模態(tài)主導(dǎo)，其設(shè)計目標(biāo)是在飛機(jī)的巡航馬赫數(shù)、雷諾數(shù)和后掠角范圍內(nèi)盡可能擴(kuò)大層流區(qū)。圖2 展示了風(fēng)洞試驗?zāi)Ｐ汀D3 展示了機(jī)翼沿展向的測試站位。圖4 展示了機(jī)翼典型壓力分布設(shè)計示意圖，壓力分布在前緣具有快速加速區(qū)，抑制橫流模態(tài)的增長，之后保持機(jī)翼中段為弱順壓梯度區(qū)，抑制T-S 模態(tài)增長。前緣橫流模態(tài)和機(jī)翼中段T-S 模態(tài)擾動抑制是層流設(shè)計的關(guān)鍵，標(biāo)模CRM-NLF 試驗結(jié)果表明以層流區(qū)弦向長度為特征長度的轉(zhuǎn)捩雷諾數(shù)最高可達(dá)1 000 萬。

圖2 CRM-NLF 風(fēng)洞試驗?zāi)Ｐ蚚34]Fig.2 Test model of CRM-NLF[34]

圖3 CRM-NLF 翼展位置示意圖[34]Fig.3 Spanwise stations of the CRM-NLF model[34]

圖4 機(jī)翼的典型壓力分布[34]Fig.4 Representative pressure distribution of the CRM-NLF model[34]

1.2 定常流場分析工具CFAST 簡介

CFAST 是中國商飛北研中心開發(fā)的面向工程應(yīng)用的CFD 快速分析工具。采用嵌套網(wǎng)格，機(jī)身構(gòu)造背景網(wǎng)格、機(jī)翼單獨劃分C 型網(wǎng)格，邊界層外流場采用歐拉控制方程求解。邊界層單獨劃分求解域，采用曲面非正交坐標(biāo)系，通過Keller-Box 轉(zhuǎn)化邊界層方程組為6 個常微分方程，利用預(yù)測-修正格式推進(jìn)求解。邊界層計算流程如下：在翼根和機(jī)翼駐線交點，首先按無限后掠翼假設(shè)給出初始值，之后沿翼根推進(jìn)求解駐線至翼梢，最后自駐點起始，沿流線完成邊界層全場計算。邊界層流場向無黏外流流場輸出溢出速度，使無黏外流在機(jī)翼表面的流向偏離物面，用于描述邊界層對流動的堵塞作用。CFAST 自駐點開始層流基本流計算，根據(jù)指定的強(qiáng)制轉(zhuǎn)捩位置轉(zhuǎn)捩，之后為湍流計算，湍流平均流動應(yīng)用零方程Cebeci-Smith 模型模擬。強(qiáng)制轉(zhuǎn)捩位置略大于試驗轉(zhuǎn)捩位置，保證試驗轉(zhuǎn)捩位置前的流場為層流。CFAST 在個人計算機(jī)上單一工況的計算時間僅約200 s，該計算速度主要是通過小量級網(wǎng)格的外流無黏與邊界層修正迭代完成的。由于邊界層計算域推進(jìn)求解一輪比無黏外流迭代求解一步的時間長數(shù)倍，因此在多步無黏外流迭代后才執(zhí)行一次邊界層計算域的計算，即邊界層修正。典型的翼身嵌套網(wǎng)格總數(shù)為100 萬量級，典型的計算收斂步數(shù)為60 步。圖5 給出了殘差收斂曲線，在網(wǎng)格殘差突變位置完成邊界層修正，殘差突變是邊界層輸出到外流的溢出速度更新導(dǎo)致的，該突變引起的殘差上升很快，在2～3 個迭代步后即可消除。

圖5 標(biāo)模CRM-NLF 殘差收斂歷程圖Fig.5 Residual convergence history diagram of the standard model CRM-NLF

在邊界層修正步，邊界層自駐點出發(fā)，沿流向求解邊界層方程并輸出邊界層速度、溫度、密度法向分布，為eN分析提供了所需的三維可壓縮邊界層信息。其中典型的密度和三維速度分布如圖6 所示。

圖6 三維可壓縮邊界層信息Fig.6 Boundary layer characteristics of three-dimensional compressible flow

1.3 轉(zhuǎn)捩預(yù)測的eN方法簡介

基于線性穩(wěn)定性理論的eN方法仍然是目前飛機(jī)設(shè)計計算中得到最廣泛認(rèn)可的轉(zhuǎn)捩預(yù)測方法。eN方法的基本思想是，邊界層內(nèi)存在的各種頻率的小幅值擾動波向下游傳播，擾動從開始放大處起，沿下游方向?qū)ζ渚€性增長率做積分，取其包絡(luò)線，從而得到N值，再通過試驗等手段對N值進(jìn)行標(biāo)定，確定轉(zhuǎn)捩判據(jù)。在所有頻率的擾動中，擾動的N值包絡(luò)最先達(dá)到轉(zhuǎn)捩判據(jù)的位置被認(rèn)為是轉(zhuǎn)捩位置。不穩(wěn)定模態(tài)擾動幅值滿足如下關(guān)系式：

式中：A為 擾動幅值；A0為擾動開始放大處的幅值；N與幅值增長倍數(shù)有關(guān)，一般認(rèn)為N是擾動頻率和流向位置的函數(shù)；αi為擾動增長率。

Malik 等[23]發(fā)現(xiàn)，假定無黏勢流方向的展向波數(shù)為實數(shù)，擾動增長率沿?zé)o黏勢流方向進(jìn)行N值計算，計算量小且精度足夠。本文eN積分路徑選取與無黏勢流方向保持一致。通常情況下，對于當(dāng)前跨聲速大型民用飛機(jī)典型后掠機(jī)翼表面三維流動問題，轉(zhuǎn)捩由T-S 模態(tài)和橫流駐波模態(tài)共同主導(dǎo)。本文基于線性穩(wěn)定性理論的雙eN方法對N進(jìn)行轉(zhuǎn)捩預(yù)測分析，雙eN方法的描述詳見文獻(xiàn)[26]。對于T-S 模態(tài)，采用固定波角法，沿著無黏勢流方向?qū)υ鲩L率進(jìn)行積分獲得不同頻率下的N值包絡(luò)線；對于橫流駐波模態(tài)，頻率取0，采用固定展向波數(shù)法，獲得不同展向波數(shù)下的N值包絡(luò)線。

1.4 算例驗證

選取標(biāo)模NLF(2)-0415，圖7 對比展示了Re=2.37×106、迎角為-4°時計算和試驗結(jié)果[35]，由圖可知，本文CFAST 計算結(jié)果與試驗結(jié)果基本吻合，說明CFAST 能準(zhǔn)確得出基本流結(jié)果。

圖7 以當(dāng)?shù)赝饩壦俣葻o量綱化的流向速度分布Fig.7 Wall-normal profiles of the streamwise velocity

2 計算參數(shù)與網(wǎng)格設(shè)置

采用CFAST 完成無黏外流和邊界層耦合計算，機(jī)翼下表面全湍流，上表面基于CRM-NLF 的低湍流風(fēng)洞試驗結(jié)果[33]，選取15 個不同展向位置，結(jié)合試驗轉(zhuǎn)捩點給定強(qiáng)制轉(zhuǎn)捩位置。計算設(shè)置的強(qiáng)制轉(zhuǎn)捩點較試驗轉(zhuǎn)捩點向后延長了2%弦長，便于eN分析計算試驗轉(zhuǎn)捩點的N值。不同展向網(wǎng)格對應(yīng)站位η的轉(zhuǎn)捩點由給定的強(qiáng)制轉(zhuǎn)捩位置插值擬合得出，其中η表示展向站位距機(jī)身軸線距離與翼梢距機(jī)身軸線距離之比。具體計算工況參見表1。

表1 計算工況Table 1 Calculated parameters

本文無黏外流網(wǎng)格設(shè)置如圖8（a）所示，網(wǎng)格總數(shù)110 萬。機(jī)翼邊界層第一層網(wǎng)格高度1.1×10-7m，法向增長率1.26，邊界層內(nèi)法向網(wǎng)格點數(shù)沿弦向逐漸增長，前緣處為21，尾緣處增至46，以滿足邊界層厚度沿弦向的增長需求，機(jī)身不計算邊界層。Venkatachari等[36]基于RANS 計算網(wǎng)格設(shè)置如圖8（b）所示，網(wǎng)格總數(shù)1.2 億，邊界層第一層網(wǎng)格高度2.54×10-6m，法向增長率1.02，法向網(wǎng)格數(shù)200，邊界層內(nèi)網(wǎng)格數(shù)約40，機(jī)身計算邊界層。

圖8 網(wǎng)格分布Fig.8 Grid distribution

3 結(jié)果與分析

3.1 基本流場

圖9 展示了CFAST 帶轉(zhuǎn)捩計算獲得的俯視表面壓力系數(shù)Cp云圖，并與Paredes 等[37]的計算結(jié)果進(jìn)行了對比，CFAST 半模計算給出了左側(cè)機(jī)身，Paredes等[37]給出了右側(cè)機(jī)身，兩者的結(jié)果具有相同的變化特征，翼根先向外伸出第一道激波，機(jī)翼拐折后，外翼段存在第二道激波。固定迎角變雷諾數(shù)（Case 1 和Case 2）時，激波的形態(tài)保持一致，激波線基本不變；固定雷諾數(shù)變迎角（Case 1 和Case 3）時，隨迎角的增大，內(nèi)翼段展向向外延伸出的激波更長，激波有較明顯的變化。Case 2 和Case 3 的結(jié)果圖類似。

圖9 帶轉(zhuǎn)捩計算的表面壓力系數(shù)云圖（Case 1）Fig.9 Contours of surface pressure coefficient in the imposed transition front solutions for Case 1

圖10 進(jìn)一步展示了不同翼展位置的表面壓力系數(shù)Cp的弦向分布，并與NASA 的試驗結(jié)果[33]和Paredes 等[37]的RANS 結(jié)果進(jìn)行了對比。本文與Paredes等[37]基于RANS 的轉(zhuǎn)捩計算結(jié)果基本重合，機(jī)翼激波位置較試驗略有后移，機(jī)翼吸力面激波前壓力整體略低于試驗，其余與試驗吻合較好。圖11 展示了表面摩擦系數(shù)分布，對于摩擦系數(shù)突增的轉(zhuǎn)捩點和摩擦系數(shù)陡降的激波位置，二者計算誤差在5%弦長以內(nèi)，說明CFAST 基本流計算結(jié)果準(zhǔn)確可靠。需要指出，Paredes 等[37]的網(wǎng)格總數(shù)1.2 億，本文網(wǎng)格總數(shù)110 萬，單工況計算時間僅為200 s，CFAST 能準(zhǔn)確捕捉基本流場的同時計算效率明顯提高。

圖10 不同翼展位置的表面壓力系數(shù)（Case 1）Fig.10 Experimental and predicted pressure coefficients at different sections of the wing for Case 1

圖11 表面摩擦系數(shù)分布（Case 3）Fig.11 Friction coefficient distributions at different sections of the wing for Case 3

3.2 轉(zhuǎn)捩預(yù)測分析

標(biāo)模CRM-NLF 轉(zhuǎn)捩由橫流駐波模態(tài)和T-S 模態(tài)主導(dǎo)。圖12 展示了標(biāo)模CRM-NLF 機(jī)翼上表面T-S波（紅實線）和橫流駐波（藍(lán)實線）不穩(wěn)定性擾動放大因子N增長曲線，圖中黑實線表示風(fēng)洞試驗的轉(zhuǎn)捩位置，黑點畫線表示CFAST 計算獲得的機(jī)翼上表面壓力系數(shù)分布。由圖可知，在展向37%處站位RowC（機(jī)翼拐折）以內(nèi)，轉(zhuǎn)捩由T-S 模態(tài)主導(dǎo)，站位I在接近激波附近的區(qū)域T-S 模態(tài)快速增長，可能是由激波誘導(dǎo)層流分離導(dǎo)致的。需要指出，T-S 模態(tài)和橫流駐波模態(tài)的N均不超過7，說明標(biāo)模CRM-NLF 機(jī)翼的壓力分布通過機(jī)翼前緣抑制橫流駐波模態(tài)和機(jī)翼中段抑制T-S 模態(tài)可以兼顧控制T-S 模態(tài)和橫流駐波模態(tài)的增長幅度，使機(jī)翼表面維持較長的層流區(qū)。

圖12 不穩(wěn)定模態(tài)沿流向的變化（Case 1）Fig.12 Streamwise variations of the N-factors of T-S and stationary crossflow modes for Case 1

eN方法是一種半經(jīng)驗的轉(zhuǎn)捩預(yù)測方法，轉(zhuǎn)捩N值可以根據(jù)試驗轉(zhuǎn)捩位置進(jìn)行標(biāo)定，也可以根據(jù)事先給定的轉(zhuǎn)捩閾值 [(NTS)tr,(NCF)tr]，計算機(jī)翼的轉(zhuǎn)捩位置，然后與試驗結(jié)果進(jìn)行比較。轉(zhuǎn)捩由T-S 模態(tài)和橫流駐波模態(tài)中先達(dá)到轉(zhuǎn)捩閾值的不穩(wěn)定模態(tài)主導(dǎo)。風(fēng)洞試驗來流湍流度為Tu=0.24%，根據(jù)Mack[38]的經(jīng)驗關(guān)系式N=-8.43-2.4ln(Tu)，轉(zhuǎn)捩閾值為[(NTS)tr,(NCF)tr]=[6.0,6.0]。風(fēng)洞試驗研究了巡航馬赫數(shù)下迎角和雷諾數(shù)對轉(zhuǎn)捩的影響，本文在風(fēng)洞試驗的基礎(chǔ)上還分析了馬赫數(shù)對轉(zhuǎn)捩的影響。

3.2.1 迎角對轉(zhuǎn)捩的影響分析

圖13 和圖14 分別展示了迎角為1.5 和2.0 時擾動放大因子N的云圖，并對比展示了Paredes 等[37]依據(jù)三維PSE 得出的計算結(jié)果，圖中黑實線表示風(fēng)洞試驗測量得到的轉(zhuǎn)捩線，黑虛線表示數(shù)值計算得出的激波線，實心圓表示轉(zhuǎn)捩閾值[ (NTS)tr,(NCF)tr]為 [6.0,6.0]得出的轉(zhuǎn)捩位置，與Lynde 等[34]和Paredes 等[37]進(jìn)行轉(zhuǎn)捩預(yù)測分析得出的轉(zhuǎn)捩閾值一致。由圖可知，本文計算得出的轉(zhuǎn)捩位置和試驗轉(zhuǎn)捩位置基本一致。Lynde 等[34]和Paredes 等[37]分析風(fēng)洞試驗結(jié)果TSP照片認(rèn)為機(jī)翼內(nèi)段可能為T-S 不穩(wěn)定模態(tài)誘發(fā)轉(zhuǎn)捩。本文計算結(jié)果同樣表明，機(jī)翼內(nèi)段的轉(zhuǎn)捩由T-S 不穩(wěn)定模態(tài)主導(dǎo)。在機(jī)翼外段，橫流駐波模態(tài)的N值明顯較大，而此時轉(zhuǎn)捩線與激波線基本保持一致，轉(zhuǎn)捩可能是由激波導(dǎo)致的。對于T-S 模態(tài)，本文計算的N值與Paredes 等[37]計算結(jié)果基本保持一致，吻合很好；對于橫流駐波模態(tài)，本文計算的N值在機(jī)翼中段整體略大于Paredes 等[37]的結(jié)果，而在機(jī)翼外段整體略小于Paredes 等[37]的結(jié)果，演化趨勢一致，這可能是由PSE 考慮非平行性和機(jī)翼表面曲率導(dǎo)致的。隨著迎角的增大，在翼展的大部分區(qū)域，橫流駐波模態(tài)擾動放大因子N明顯減小，轉(zhuǎn)捩主要由T-S 模態(tài)主導(dǎo)。總之，本文基于CFD 快速分析工具CFAST進(jìn)行線性穩(wěn)定性分析能較為準(zhǔn)確的進(jìn)行轉(zhuǎn)捩預(yù)測分析。

圖13 橫流駐波和T-S 波N 值的等值線云圖（Case 1）Fig.13 N-factor contours of stationary CF waves and T-S waves for Case 1

圖14 橫流駐波和T-S 波N 值的等值線云圖 (Case 3)Fig.14 N-factor contours of stationary CF waves and T-S waves for Case 3

3.2.2 雷諾數(shù)對轉(zhuǎn)捩的影響分析

圖15 和圖16 分別展示了雷諾數(shù)為12.5×106、15.0×106、17.5×106和20.0×106時橫流駐波模態(tài)和T-S 模態(tài)擾動放大因子N的等值線云圖，其他參數(shù)與表1 中的Case 1 一致。圖中黑實線表示試驗轉(zhuǎn)捩位置，黑虛線表示激波線，實心圓表示轉(zhuǎn)捩閾值[(NTS)tr,(NCF)tr]為 [6.0,6.0]得出的轉(zhuǎn)捩位置。當(dāng)雷諾數(shù)不超過17.5×106時，本文計算得出的轉(zhuǎn)捩位置和試驗轉(zhuǎn)捩位置基本一致，雷諾數(shù)為20.0×106時，本文計算得出的轉(zhuǎn)捩位置和試驗轉(zhuǎn)捩位置偏差增大。隨雷諾數(shù)從12.5×106增大到15.0×106，轉(zhuǎn)捩位置變化不明顯，機(jī)翼內(nèi)段的轉(zhuǎn)捩由T-S 不穩(wěn)定模態(tài)主導(dǎo)，機(jī)翼外段的轉(zhuǎn)捩可能是由激波導(dǎo)致的；而隨雷諾數(shù)繼續(xù)增大到17.5×106和20.0×106，轉(zhuǎn)捩位置明顯前移，橫流駐波模態(tài)N值明顯增大，而轉(zhuǎn)捩位置處的T-S 模態(tài)N值有所減小，機(jī)翼外段的轉(zhuǎn)捩可能由橫流駐波模態(tài)主導(dǎo)。造成風(fēng)洞試驗轉(zhuǎn)捩位置前移的原因可能是：一方面隨雷諾數(shù)增大，流動容易不穩(wěn)定，轉(zhuǎn)捩位置會前移；另一方面，邊界層厚度隨雷諾數(shù)的增大而變薄，風(fēng)洞品質(zhì)和機(jī)翼表面瑕疵等因素對轉(zhuǎn)捩的影響會明顯增強(qiáng)，當(dāng)雷諾數(shù)增大到20.0×106時，可能在前緣直接誘發(fā)轉(zhuǎn)捩。

圖15 橫流駐波模態(tài)N 值云圖Fig.15 N-factor contours of cross-flow modes at different ReMAC

圖16 T-S 模態(tài)N 值云圖Fig.16 N-factor contours of T-S modes at different ReMAC

3.2.3 馬赫數(shù)對轉(zhuǎn)捩的影響分析

在風(fēng)洞試驗基礎(chǔ)上，本節(jié)分析了馬赫數(shù)對轉(zhuǎn)捩的影響，其他參數(shù)與表1 中的Case 1 一致。圖17 展示了不同馬赫數(shù)時表面壓力系數(shù)Cp分布。由圖可知，隨馬赫數(shù)從0.856 降為0.820 和0.800，Cp分布變化明顯，機(jī)翼中段由弱順壓梯度區(qū)變?yōu)槟鎵禾荻葏^(qū)，激波消失。隨馬赫數(shù)從0.856 增加到0.870，激波位置后移，激波強(qiáng)度有所增強(qiáng)。圖18 和圖19 分別展示了橫流駐波模態(tài)和T-S 模態(tài)擾動放大因子N值云圖。與圖13（Ma=0.856）相比，隨馬赫數(shù)增大，橫流駐波模態(tài)N值逐漸增大，而T-S 模態(tài)N值明顯減小。圖中黑點劃線表示 [(NTS)tr,(NCF)tr]=[6.0,6.0]作為轉(zhuǎn)捩閾值得出的轉(zhuǎn)捩位置。隨馬赫數(shù)增大，內(nèi)翼段轉(zhuǎn)捩位置明顯后移，而外翼段層流區(qū)長度超過50%弦長，變化不明顯。

圖17 不同翼展位置的表面壓力系數(shù)Fig.17 Pressure coefficients at different sections of the wing

圖18 橫流駐波模態(tài)N 值云圖Fig.18 N-factor contours of stationary cross-flow modes

圖19 T-S 模態(tài)N 值云圖Fig.19 N-factor contours of T-S mode

4 結(jié)論

標(biāo)模CRM-NLF 的后掠機(jī)翼邊界層是典型的三維邊界層，其轉(zhuǎn)捩由橫流駐波模態(tài)和T-S 模態(tài)主導(dǎo)。本文基于跨聲速自然層流標(biāo)模CRM-NLF 的風(fēng)洞試驗數(shù)據(jù)，結(jié)合內(nèi)部快速CFD 分析工具CFAST 和eN轉(zhuǎn)捩計算工具進(jìn)行了轉(zhuǎn)捩預(yù)測分析。總結(jié)如下：

1）通過轉(zhuǎn)捩N值標(biāo)定，研究發(fā)現(xiàn)，在內(nèi)翼段，轉(zhuǎn)捩由T-S 模態(tài)主導(dǎo)；在外翼段，轉(zhuǎn)捩由橫流駐波模態(tài)和激波共同主導(dǎo)。隨著迎角和雷諾數(shù)的增大，轉(zhuǎn)捩位置會有所前移。采用雙eN方法，通過湍流度結(jié)合Mack公式確定T-S 模態(tài)和橫流駐波模態(tài)的轉(zhuǎn)捩N值為6，預(yù)測的轉(zhuǎn)捩位置與風(fēng)洞試驗結(jié)果吻合很好，能較為準(zhǔn)確的進(jìn)行轉(zhuǎn)捩預(yù)測分析；

2）在風(fēng)洞試驗的基礎(chǔ)上本文分析了馬赫數(shù)對轉(zhuǎn)捩的影響。馬赫數(shù)在0.800～0.870 之間時，隨馬赫數(shù)增大，內(nèi)翼段的轉(zhuǎn)捩位置明顯后移，而外翼段層流區(qū)長度超過50%弦長，變化不明顯；

本文旨在開發(fā)一套適用于復(fù)雜三維外形的快速轉(zhuǎn)捩預(yù)測工具，經(jīng)CRM-NLF 標(biāo)模驗證CFAST 在個人計算機(jī)上工況的基準(zhǔn)分析時間約200s，且轉(zhuǎn)捩預(yù)測精度足夠，初步達(dá)成研究目標(biāo)，在后續(xù)層流設(shè)計研究中，有助于加速后掠翼層流優(yōu)化設(shè)計的迭代效率。后續(xù)也將進(jìn)一步研究來流湍流度和表面制造因素對轉(zhuǎn)捩的影響。