軸對稱短艙自然層流優化及轉捩敏感性分析

曹凡，張美芳，胡驍，唐智禮

南京航空航天大學 航空學院，非定常空氣動力學與流動控制工業和信息化部重點實驗室，南京 210016

發動機是飛機的核心部件之一，為了保護發動機，降低阻力和噪聲，短艙應運而生。短艙主要由進氣道、整流罩和內、外涵道噴管組成，是推進系統的重要組成部分[1]。隨著渦扇發動機的發展，涵道比不斷增大，短艙直徑也越來越大，勢必帶來阻力的上升，影響飛行效率。面積越大，短艙的設計難度就越高，大涵道比發動機短艙是綜合了氣動設計、結構設計、材料和工藝等多種學科的復雜系統。為了減小短艙面積增大帶來的摩擦阻力，通過氣動設計的手段對大涵道比發動機短艙進行外形層流化設計是近年來興起的有效手段[2]。

相同來流條件下，層流的摩阻遠小于湍流摩阻，如果將層流技術應用于機翼、尾翼和短艙，能夠減少整機阻力15%以上，因此，層流減阻技術在新型民機的氣動設計中應用前景廣闊。自然層流技術是一種層流主動控制方法，該技術通過對短艙型面的設計，使得短艙外表面保持較大的順壓梯度區域，從而抑制流向T-S(Tollmien-Schlichting)波的增長，推遲轉捩的發生保持大范圍層流區域[3]。目前，機翼和機身的優化減阻技術比較成熟，但對短艙的研究較少，在現有的經驗和技術基礎上對短艙進行優化，實現有效的減阻仍然是主要攻關問題。

近年來，國內學者對NLF短艙氣動外形設計也開展了一系列研究，但相關的風洞和飛行試驗卻很少。西北工業大學的何小龍等[12]研究了帶動力的自然層流短艙氣動優化設計；孟曉軒等[13]基于雙eN方法研究了CRM(Common Research Model)短艙流動轉捩的影響因素，認為馬赫數和攻角對轉捩位置影響較大,而雷諾數和湍流度相對較小。2019年,復旦大學Wang等[14]應用差分進化算法優化NLF短艙，將層流范圍提升了16.64%，獲得良好減阻效果;杜璽等[15]針對設計的自然層流短艙在德國ETW(European Transonic Wind)跨聲速風洞進行了對比驗證，結果表明低雷諾下轉捩預測結果吻合較好，而高雷諾數下轉捩提前發生與CFD(Computational Fluid Dynamics)結果差別明顯。目前來看，針對NLF短艙的氣動外形設計已經取得較好的進展，但實際設計過程中，識別對層流影響顯著的變量更加重要，因此需進行參數靈敏度分析。高雷諾數下，保持大范圍的層流流動極其困難，在短艙層流化設計中，外形幾何控制參數對轉捩位置的影響和氣動參數與層流范圍的關系，尚有待研究。

以CRM短艙為研究對象，在高雷諾數跨聲速狀態滿足幾何約束條件下利用發展的進化/確定性混合算法進行氣動外形優化，分別以總阻力最小化和層流最大化設計軸對稱NLF短艙，驗證層流短艙減阻的可能性。為了分析幾何控制參數、流量系數和馬赫數對層流范圍的影響，基于層流最大化的短艙構型開展轉捩敏感性研究，分析層流與阻力系數的關系，以期為后續自然層流短艙的設計和魯棒優化提供依據。

1 層流轉捩預測

(1)

Langtry和Menter[17]采用的關聯函數為

(2)

(3)

(4)

式中：U為當地速度；θ為動量厚度；s為流線弧長；Ui為速度張量分量；xi為笛卡爾坐標張量分量；μ為層流黏性系數；μt為湍流黏性系數；S為流線弧長；k為湍動能；Tu為湍流度；λθ為壓力梯度；ρ為密度；Pγ、Eγ為生成項、耗散項；σr=1。

(5)

(6)

(7)

式中：Fwake為確保Fθ t在翼型下游尾流區為0的參數；y為靠近壁面最近距離；δ為邊界層厚度；cθ t為源項系數；σθ t為擴散系數，取2.0；ce2為耗散項系數;t為時間。

σθ t的大小決定了臨界動量厚度雷諾數對流動歷史和當地流動參數的依賴程度，σθ t越大,轉捩對流動歷史越不敏感，更多地取決于湍流度和壓力梯度等當地流動參數。擴散系數和源項系數的變化雖然會帶來轉捩位置的改變，但與氣動參數和幾何參數的變化相比，其影響程度相對較小[18-19]。

混合函數Fθ t在邊界層外為0，在邊界層內為1，起到關閉Pθ t的作用。為了耦合k-ωSST湍流方程，引入有效間歇因子γeff與湍動能方程進行交互，表達式為

γeff=max(γ,γsep)

(8)

Fθ t

(9)

式中：γsep用于捕獲并模擬分離誘導轉捩；Freattach為一種開關函數，表示當黏度比大到足以引起層流再附著時，禁止其修改；Rev為渦量雷諾數；Reθc為臨界動量厚度雷諾數。

(10)

(11)

圖1 NLF(1)-0416摩擦阻力系數曲線

表1 NLF(1)-0416翼型轉捩結果對比[20-21]

DLR-F5機翼前緣后掠角20°，采用超臨界對稱翼型。計算工況:Ma=0.82，α=2°，Re=1.5×106，參考弦長0.15 m，湍流強度0.5%，湍流黏性比10，物面表面網格y+<1，第1層網格高度L=1.5×10-6m，網格量約為247萬，如圖2所示。

圖2 DLR-F5計算網格

圖3 DLR-F5試驗結果與數值模擬結果對比

2 高雷諾數下軸對稱自然層流短艙優化設計

大型民用運輸機短艙的流動雷諾數一般高于2×107，為了實現高雷諾數下的層流流動，需要保持更加強烈的順壓梯度，這遠比在低雷諾數下的優化更加困難。通過精確的參數化方法、高精度流場求解、可靠的轉捩預測結合高效的優化算法進行高雷諾數下的自然層流短艙優化設計。

2.1 軸對稱短艙幾何設計

為了設計滿足幾何約束的短艙，同時對比優化前后的層流范圍，首先對CRM標模短艙[22]進行自然通氣狀態下仿真。設計工況：巡航高度11 000 m，基于短艙參考弦長雷諾數3.5×107，自然來流馬赫數Ma=0.85，攻角α=0°，湍流強度Tu=0.05%，湍流黏性比為10，流場計算采用有限體積方法，空間離散為二階精度，通量計算采用 Roe 格式，時間推進為 LU-SGS (Lower Upper Symmetric Gauss Seidel) 隱式方法。

計算網格如圖4所示，劃分六面體結構網格共850萬，其中物面流向和周向各201個網格節點，第1層網格高度4×10-6m。

圖4 CRM短艙物面網格

計算結果如圖5所示。CRM短艙在自然通氣狀態下，其層流范圍很小，且流動在靠近唇口的位置便開始轉捩發生流動分離。

圖5 CRM短艙計算結果

經測量CRM短艙，其最大弦長6.06 m，半徑1.75 m，為了滿足幾何約束(見表2)，對CRM短艙剖面進行CST(Class Shape Transformation)參數化[23]。

表2 短艙設計約束條件

CST方法由一個類別函數和一個形狀函數構成，具有精度高和設計參數少等特點，其設計參數具有明確的物理意義，能夠直接反映翼型的變化過程。其數學表達式為

(12)

式中：c為翼型弦長；Rle為翼型前緣半徑；Δzte/c為后緣點縱坐標；β為后緣傾角；bi為形狀函數控制參數(i=1,2,…,n-1)。采用基于Bernstein多項式的CST方法時，為了提高參數化精度,且防止過高的階數使得參數化趨于病態，取n=6。此時，參數化短艙外剖面共需要8個設計變量(Rle、Δzte、β、b1、b2、b3、b4、b5)。

對于流量系數的約束，采用了閆海津和杜璽[24]提出的在短艙出口處添加堵錐的方式代替動力短艙實現流量系數的連續調節，初始流量系數(Mass Flow Ratio,MFR)MFR=0.7，最終形成的滿足約束的軸對稱短艙半模如圖6所示。

圖6 帶堵錐的軸對稱短艙

2.2 自然層流短艙優化設計方法

通常更大范圍的層流要求短艙表面維持較大范圍的順壓梯度，這會導致短艙后半段的壓力快速恢復產生強激波，增大壓阻和波阻從而對總阻力的減小產生負面效應。本文著眼于自然層流短艙的優化設計，期望在短艙的表面實現更大范圍的層流，但對于具體的可優化程度沒有經驗參考，伴隨摩擦阻力減小對激波阻力與總阻力的影響也未知。因此，為了更好地分析自然層流短艙的局部敏感性，分別以轉捩位置最大化和總阻力最小化為優化目標進行優化，分析不同參數影響下的層流范圍，驗證層流減阻的可能性。希望獲得一款能夠保持較大層流范圍的短艙進行敏感性分析，有利于研究層流范圍變化規律。

2.2.1 進化/確定性混合優化算法

解決氣動優化問題可以用3種不同的方法：① 無導數搜索算法(Derivative-Free Search Algorithns,DFSAs)；② 基于梯度的算法(Gradient-Based Algorithms,GBAs)；③ 進化算法(Evolntionan Algorithms,EAs)。每一種算法都有各自的優缺點，可以結合不同的算法[25]，把其優點同時發揮出來，比如:EAs的魯棒性和全局搜索,以及GBAs的快速收斂能力。通常使用基于梯度的搜索方法，計算目標函數的導數以確定最佳的下降方向，梯度法已與遺傳算法(Genetic Algorithm,GA)結合用來設計翼型和機翼[26-27]。

以GA為主算法，用擬牛頓方法(BFGS算法進行)輔助加速收斂，發展了一種新的多級混合優化算法。文獻[28]對該混合算法的原理進行了詳細介紹。進化和確定性算法在這種混合優化方法中交替進行，并成功應用于高超聲速飛行器氣動外形優化，提升了優化效率。GA/BFGS算法流程如下：

Step1初始化種群，并計算適應值。定義初始種群數量Np進化代數Ke，局部搜索選擇頻率為f。

Step2執行GA算法，產生較優種群。產生一個隨機整數Ns(Ns≤Np)，并從當前種群中選中并標記Ns個個體。

Step3對標記的Ns個個體執行確定性優化算法BFGS，BFGS將這Ns個個體每個迭代Kb步，產生Ns個新個體,并將其加入種群替代原有個體。

Step4最優選擇(即從GA種群中選擇前Ns個最優個體)，新的搜索空間簡化為收斂的較小空間。隨著搜索空間的縮小，種群規模和搜索頻率也可以減少。

Step5繼續對簡化的搜索空間和局部搜索頻率進行優化，直至收斂。

多級方法在這個交替的過程中耦合。為了保持基因的多樣性，可以在進化過程開始時使用大種群規模，以便選擇更好的基因并傳給下一代。該算法交替過程如圖7所示。

圖7 GA/BFGS混合優化算法示意圖

為驗證算法的有效性，對(GA/BFGS算法)與純梯度(BFGS)算法[29]、差分進化算法(Differential Evdution,DE)[30]、GA算法[31]在Ackley多峰函數上進行了對比驗證，如圖8所示。種群個體均為100、代數5 000、交叉概率0.9、變異概率0.1，差分算法中縮放因子取0.5。

Ackley多峰函數為

(13)

Ackley函數取值范圍為x∈[-32, 32]，維度d=30，具有很多局部極小值，雖然梯度算法進化效率很高，但容易陷入局部最優解，而GA、DE這種進化算法能夠利用交叉、變異算子跳出局部最優解，進化效率相對較慢。本文所使用的混合算法在進化前期效率不如梯度算法，如圖8所示，但避免了陷入局部最優解，很好地平衡了進化效率和尋優能力。

圖8 Ackley 函數收斂曲線

2.2.2 軸對稱短艙自然層流設計

自然層流短艙設計過程中，目標是在激波前保持順壓梯度，抑制T-S波轉捩，從而獲得大范圍層流區域。采用GA/BFGS混合優化算法結合六階CST(形狀類別變化函數)參數化方法，以總阻力最小化和轉捩位置最大化為目標分別進行單點優化。

短艙計算工況：巡航高度11 000 m，基于短艙最大長度雷諾數Re=2.5×107，流量系數MFR=0.7，湍流度Tu=0.05%，黏性比為10。劃分結構網格如圖9所示，第1層網格高度5×10-3m，網格前緣用C型網格保證正交性，并進行加密處理，軸向分布200個網格點，網格數量約11萬。

圖9 軸對稱短艙剖面網格

多層混合交替初始階段種群規模Np=40，局部搜索選擇頻率為f=2，進化代數100代，保證優化過程中的基因多樣性，隨著搜索空間的縮小，種群規模也可以減少，這里減少到20個，局部搜索頻率減少到1，繼續進化10代，確保最后幾代個體能夠收斂到全局最優解。

計算結果如表3所示，Baseline為短艙設計工況，minCD是以總阻力系數最小為優化目標的工況，maxXtr是以轉捩位置最大為優化目標的工況。ΔCD為短艙內外剖面的總阻力系數相對于基準外形的變化量，由于此次設計不改變內型線，其值變化可近似代表外剖面的阻力變化；Xtr為短艙轉捩位置；Cf-up為短艙外剖面的摩擦阻力系數，參考面積為383.7 m2；Dmax為短艙外罩最大直徑；Rle為短艙前緣半徑。

表3 軸對稱短艙優化前后計算結果對比

對比分析發現在總阻力系數最小時，轉捩位置推遲了3.6%，摩擦阻力也略有降低，但短艙外罩最大直徑大大降低，短艙剖面面積減小約28.5%，前緣半徑變小，如圖10所示。當層流面積最大化時，短艙外剖面的摩擦阻力減小24%，表明延遲轉捩，增大層流范圍極大地降低摩擦阻力，但也使得總阻力明顯增加。這表明合理的平衡層流面積最大化和總阻力最小化2個目標，才能實現真正的短艙減阻目的。

圖10 優化前后短艙外形對比

歸一化后Cp、Cf曲線如圖11、圖12所示，原始外形由于唇口吸力峰的存在過早轉捩，層流范圍最大化的外形在前半段維持了良好的順壓梯度，可達到總長的27%左右，但是在短艙尾部出現由激波導致的強逆壓梯度，壓力快速恢復，出現相對強烈的激波，這是導致阻力略有增加的原因。阻力最小化外形后半段外剖面迅速向內剖面靠近，壓縮了短艙面積。

圖11 優化前后壓力系數曲線對比

圖12 優化前后摩阻系數曲線對比

3 軸對稱自然層流短艙局部敏感性分析

基于單點優化設計，獲得了2種不同的短艙構型，研究表明如果以總阻力為優化目標最后的優化結果的確會使總阻力減小，但是層流轉捩位置推遲較小，層流范圍增大不明顯。單點優化的設計結果在偏離設計點時，往往不能保證飛行性能。同時，為了提取層流短艙的設計規律，分析了幾何參數和飛行參數的敏感性。為了更好地分析單變量影響下軸對稱自然層流短艙的轉捩變化規律，選取層流范圍最大的短艙構型為初始構型，較大范圍的層流更容易對比分析不同參數帶來的影響。

3.1 局部單變量幾何敏感性

由分析可知，短艙剖面采用六階CST參數化方法進行參數化，考慮幾何約束去掉后緣厚度Δzte后還有7個參數，分別為前緣半徑Rle、后緣傾角β、短艙剖面形狀函數的控制參數(b1、b2、b3、b4、b5),控制范圍如圖13所示。在保持飛行狀態不變的情況下，對CST參數化設計變量初始值的(-15%, 15%)區間如表4所示，取6組數據進行單一變量的幾何敏感性分析，并給出Cp以及歸一化后轉捩位置,如表5、圖14所示。

表4 設計變量取值范圍

通過對比表5的數據，按對轉捩位置影響大小進行排序，轉捩位置隨b1、b2、β的變化最為明顯，其次為b5、Rle，影響最小的參數為b3、b4。其中b1對層流范圍的影響最為顯著，當其降低10%時，層流范圍降至10%以下。由圖13可知,b1主要控制短艙截面的頭部區域，在短艙的設計加工中需要重點考慮。b2、β變化下，層流范圍在10%、20%內波動較大，是次要考慮的設計點。而b3、b4的變化對轉捩位置的影響不敏感，層流范圍都可以保持在20%以上。

圖13 形狀控制參數控制范圍

表5 不同幾何參數變化下的轉捩位置

通過分析層流短艙的初始構型壓力分布如圖14 中所示，在短艙截面前半部分保持較長范圍的順壓梯度，緊接著出現由激波導致的強逆壓梯度，壓力迅速恢復。跨聲速狀態下，流動從翼型前緣開始，穿過一系列膨脹波加速到超聲速狀態，而后保持超聲速流動直到后緣附近。由于壓力需要恢復到后駐點處的壓力值，同時由超聲速狀態變為亞聲速狀態，因此在層流短艙截面后緣附近產生一道激波，通過逆壓梯度的大小可以評估激波的強度。

圖14 自然層流短艙幾何參數變化下的Cp曲線

隨著b1的增大，曲線吸力峰值減小，外形表面的壓力變化更為平緩，轉捩位置變化較小，但b1減小時吸力峰變強使得提前發生轉捩，層流范圍變小；b2對前緣吸力峰的影響也較大。同時由于轉捩位置變化較大，導致后緣處激波位置和強度變化也較大；b3、b4、b5、β在吸力峰后的壓力變化趨于平緩，對壓力分布的影響較小；隨著前緣半徑Rle增大前緣吸力峰變強，對激波強度的影響較小。

3.2 流量系數敏感性

MFR是指進入進氣道的實際空氣質量流量與同馬赫數下進入進氣道的理論最大質量流量之比[24]。本文優化后的軸對稱層流短艙MFR=0.699，在(0.67，0.73)范圍內通過移動堵錐實現流量系數的調節，轉捩位置變化如圖15所示。

圖15 轉捩位置隨流量系數變化關系

流量系數增加的過程中，轉捩位置的大小波動明顯，MFR=0.678時，轉捩位置前移到0.4 m以下。可以發現，巡航狀態下進氣流量對層流影響很大，此構型下盡量保持MFR>0.69才能獲得較大范圍層流。

3.3 馬赫數敏感性

短艙能否在不同馬赫數下產生較為強烈的順壓梯度以保持穩定的層流范圍以適應飛行狀況的變化是非常關鍵的。此次馬赫數的單變量影響分析圍繞巡航點依次在(0.8, 0.9)內取了10個工況，并在亞聲速區間取4個工況,分析該自然層流短艙的變化規律，結果如圖16～圖18所示。

圖16 轉捩位置隨馬赫數變化關系

圖17 Cp隨馬赫數變化曲線

圖18 Cf隨馬赫數變化曲線

壓力和阻力系數曲線都對馬赫數變化體現出了明顯的敏感性。需要特別指出的是當Ma=0.87 時，轉捩位置出現強烈的波動，經重復計算排除了計算意外。除此之外，馬赫數越大順壓梯度越強，吸力峰小幅降低，激波增強且后移。尤其是在飛行速度相差較大的情況，轉捩位置在亞音速區間比跨聲速區間更加敏感。跨聲速狀態下，隨著馬赫數增大，短艙整體上能夠保持層流范圍大于20%，層流范圍波動較小。

4 結 論

在幾何約束條件下,針對短艙剖面開展軸對稱自然層流短艙優化設計研究，并利用CST參數化后帶有物理意義的幾何參數，分析了幾何擾動、流量系數和馬赫數的轉捩敏感性。主要得到以下結論：

2) 基于層流最大化的軸對稱短艙構型，開展了局部幾何控制參數敏感性研究。外形控制參數b1、b2對壓力分布和層流范圍影響最大，而前緣半徑Rle、后緣傾角β對壓力分布影響很小，是層流范圍的次要影響因素。b3、b4、b5的敏感性相對較低，其主要影響短艙外表面湍流段的壓力分布。為保持大范圍層流，要重點關注短艙剖面頭部外形的設計。形成了基于CST控制參數的短艙剖面外形設計準則。

3) 轉捩位置會隨著流量系數和馬赫數的增大呈波動狀態，敏感性較高。在設計點(MFR=0.7，Ma=0.85)附近，層流范圍能夠保持在20%以上，降低了摩擦阻力，除了Ma=0.87以外，設計點附近均能保持較大范圍的層流，為后續魯棒性優化設計提供指導。