牛笑天，李杰，周智鵬，楊釗，昌陌塵

西北工業(yè)大學(xué) 航空學(xué)院，西安 710072

隨著石油資源日益枯竭，環(huán)境保護(hù)和節(jié)能減排越來越被人們關(guān)注，國際航空運(yùn)輸協(xié)會提出了航空工業(yè)減少排放物和降低噪聲的新要求。對于民航客機(jī)，降低發(fā)動機(jī)功率是減少排放的重要途徑，發(fā)動機(jī)功率主要用來克服飛機(jī)飛行的阻力。在飛機(jī)的零升阻力中，型阻和摩擦阻力一般各占50%[1-2]。當(dāng)飛機(jī)的氣動布局確定后，壓差阻力可降空間已經(jīng)不多，然而摩擦阻力還有很大的下降空間[3-4]。由于層流附面層所引起的摩擦阻力比湍流附面層要小很多，減小摩阻非常有效的方法是：盡可能擴(kuò)大飛行器表面層流區(qū)范圍，延遲轉(zhuǎn)捩現(xiàn)象的發(fā)生；通過適當(dāng)外形設(shè)計或使用一些流動干預(yù)，形成有利于流動穩(wěn)定、抑制或推遲轉(zhuǎn)捩的壓力分布，即層流減阻技術(shù)[5-6]。

層流減阻技術(shù)是在目前諸多設(shè)計約束下民機(jī)減阻設(shè)計的重要可行技術(shù)之一[7]，因此成為飛行器設(shè)計者重點(diǎn)研究的方向。早期航空制造加工能力有限，機(jī)翼表面光潔度低，飛機(jī)層流特性不明顯。隨著航空工業(yè)設(shè)計技術(shù)和制造工藝的進(jìn)步，層流流動設(shè)計逐步成為可能。對于民機(jī)而言，自然層流技術(shù)(Natural Laminar Flow, NLF)可以在低速狀態(tài)下(馬赫數(shù)Ma=0.2附近)使翼面上具有60%左右的當(dāng)?shù)叵议L可以維持層流狀態(tài)，在高速狀態(tài)下(Ma=0.754附近)使翼面上具有40%左右的弦長可以維持層流狀態(tài)，降低摩擦阻力。全機(jī)采用NLF技術(shù)的話，可以降低摩擦阻力30%左右，減少整機(jī)的阻力達(dá)15%以上[8-10]。氣動收益非常明顯，更加經(jīng)濟(jì)。民機(jī)油耗每降低一個百分點(diǎn)，都可以為航空公司減少大量運(yùn)營成本。而且油耗降低帶來的好處還體現(xiàn)在環(huán)保方面，極大地減少了碳排放量，有利于減輕對大氣的污染。

美國國家航空航天局(National Aeronautics and Space Administration,NASA)、德國航空航天中心和法過宇航局等國外機(jī)構(gòu)早在20世紀(jì)80年代已經(jīng)開始對NLF技術(shù)進(jìn)行研究，并取得了一系列成果，應(yīng)用在Boeing757和A320上面[11-13]。日本“本田飛機(jī)”(Honda Jet)采用了自然層流機(jī)翼和自然層流機(jī)身頭部設(shè)計[14-15],在2003年的首飛中達(dá)到了預(yù)期的目標(biāo)和要求。

國內(nèi)目前對層流預(yù)測方法方面已經(jīng)開展了大量的研究工作。北京大學(xué)通過風(fēng)洞試驗(yàn)，對某錐形體進(jìn)行研究，結(jié)合理論分析，證明多孔滲透表面會有效抑制高超聲速邊界層氣動加熱并推遲表面轉(zhuǎn)捩位置[16-17]。西北工業(yè)大學(xué)開發(fā)出基于線性穩(wěn)定性理論的eN轉(zhuǎn)捩預(yù)測方法(N為擾動積累放大因子)，并用該方法進(jìn)行了無限展長機(jī)翼的研究和混合層流翼身組合體的研究[18-19]。中國空氣動力研究與發(fā)展中心基于γ-Reθ轉(zhuǎn)捩模型開展了一些標(biāo)定和應(yīng)用方面的研究[20]。北京航空航天大學(xué)分別將Walters和Menter的模型應(yīng)用到高超聲速流動的轉(zhuǎn)捩預(yù)測[21]。但現(xiàn)階段國內(nèi)在層流減阻技術(shù)的風(fēng)洞、飛行試驗(yàn)研究和工程應(yīng)用方面還相對滯后，利用現(xiàn)有飛機(jī)或針對性設(shè)計的飛機(jī)進(jìn)行層流技術(shù)飛行驗(yàn)證具有迫切的現(xiàn)實(shí)需求。層流技術(shù)核心為控制邊界層轉(zhuǎn)捩。邊界層轉(zhuǎn)捩對邊界層的發(fā)展、摩擦阻力、流動分離位置等有重大影響，同時，雷諾數(shù)對轉(zhuǎn)捩又有重要影響，對層流區(qū)域的預(yù)測以及阻力的精確評估，必然要考慮湍流模型以及轉(zhuǎn)捩的影響[22-23]?？缏曀亠w行狀態(tài)下各種因素均會對層流特性產(chǎn)生一定的影響，為了保證驗(yàn)證機(jī)達(dá)到既定的驗(yàn)證目標(biāo)，需要深入研究分析各種飛行參數(shù)對于層流流動特性的具體影響。然而目前針對跨聲速狀態(tài)下的層流特性敏感性研究相對較少。

本文采用RANS方法結(jié)合基于當(dāng)?shù)刈兞康摩?Reθ轉(zhuǎn)捩預(yù)測模型，針對某特殊布局形式的層流機(jī)翼驗(yàn)證機(jī)展開跨聲速層流氣動特性和參數(shù)敏感性分析。通過三維機(jī)翼構(gòu)型DLR F-5對基于RANS的轉(zhuǎn)捩預(yù)測方法進(jìn)行算例驗(yàn)證分析。通過層流驗(yàn)證機(jī)中央驗(yàn)證段數(shù)值模擬結(jié)果與對應(yīng)試驗(yàn)數(shù)據(jù)的對比分析，來驗(yàn)證所采用轉(zhuǎn)捩預(yù)測方法的計算能力和精度。文中的數(shù)值模擬研究重點(diǎn)關(guān)注全機(jī)巡航狀態(tài)附近的氣動特性和中央驗(yàn)證段在不同飛行狀態(tài)下的表面轉(zhuǎn)捩位置及層流區(qū)長度。通過計算，分析層流流動對全機(jī)升阻及力矩特性的影響，并總結(jié)馬赫數(shù)、雷諾數(shù)、自由來流湍流度和迎角等關(guān)鍵流動參數(shù)對于翼段表面轉(zhuǎn)捩位置的影響規(guī)律。

1 計算方法

計算采用NASA的CFL3D求解器，該求解器采用格心格式的有限體積法對控制方程進(jìn)行離散。對于空間離散，無黏項(xiàng)采用Roe平均通量差分分裂格式(Flux Differences Splitting,FDS)，單元界面上差值模板為三階MUSCL(Monotonic Upstream-Centered Scheme for Conversation Laws)格式，采用min-mod限制器防止在間斷處出現(xiàn)數(shù)值震蕩。黏性項(xiàng)采用中心差分格式。時間推進(jìn)方式采用近似因子分解(AF)隱式時間推進(jìn)算法。

計算采用的湍流模型為基于k-ωSST(Shear Stress Transfer)模型的四方程γ-Reθ轉(zhuǎn)捩模型，該模型利用參數(shù)k-ω對標(biāo)準(zhǔn)兩方程k-ωSST模型中湍動能方程的生成項(xiàng)與破壞項(xiàng)進(jìn)行修正來模擬轉(zhuǎn)捩過程。標(biāo)準(zhǔn)形式的兩方程k-ωSST模型[24]為

(1)

(2)

(3)

式中：ρ為密度；k為湍動能；ω為比湍流耗散率；uj為各方向速度；xj為空間坐標(biāo)；τij為切向應(yīng)力；Sij為應(yīng)變率張量；μ為動力黏度；μt、νt為湍流黏度；Ω為渦量；β*、a1、σω2為模型常數(shù)，分別取值為0.09、0.31、0.86；F1、F2為模型中的混合函數(shù)；γ、β、σk、σω為模型參量，通過混合函數(shù)F1計算。

以兩方程k-ωSST模型為基礎(chǔ)，四方程轉(zhuǎn)捩模型添加間歇因子輸運(yùn)方程和動量厚度雷諾數(shù)輸運(yùn)方程，利用作用于湍動能方程的生成項(xiàng)與破壞項(xiàng)上的間歇因子γ實(shí)現(xiàn)模型結(jié)合[25]。

湍動能方程修改后的形式為

min(max(γeff,0.1),1.0)ρβ*ωk+

(4)

式中:γeff通過間歇因子γ計算，間歇因子輸運(yùn)方程為

(5)

式中:Pγ為生成項(xiàng)；Eγ為破壞項(xiàng)。二者計算公式為

Pγ=FlengthCa1ρS(γFonset)0.5(1-Ce1γ)

(6)

Eγ=Ca2ρΩFturb(Ce2γ-1)

(7)

其中：S為應(yīng)變率；Flength為控制轉(zhuǎn)捩長度的經(jīng)驗(yàn)函數(shù)；Fonset為控制轉(zhuǎn)捩起始位置函數(shù)；Ca1、Ca2、Ce1、Ce2為間歇方程常數(shù)。

動量厚度雷諾數(shù)輸運(yùn)方程

(8)

(9)

2 轉(zhuǎn)捩模型算例驗(yàn)證

2.1 計算模型

利用DLR F-5機(jī)翼,分別對標(biāo)準(zhǔn)兩方程γ-ReθSST模型和四方程γ-Reθ轉(zhuǎn)捩模型對于三維機(jī)翼跨聲速流動數(shù)值模擬和轉(zhuǎn)捩預(yù)測的能力進(jìn)行驗(yàn)證。如圖1所示，試驗(yàn)?zāi)Ｐ褪且粋€安裝在風(fēng)洞側(cè)壁的機(jī)翼模型。白色虛線為壓力分布提取截面，Y/Span為機(jī)翼橫截面所在展向位置百分比。翼根處Y/Span=0%，翼梢處Y/Span=100%。Sobieczky[26]于1994年開展了該試驗(yàn)的相關(guān)工作。機(jī)翼的根部平滑過渡至壁面，避免機(jī)翼根部出現(xiàn)馬蹄渦。試驗(yàn)測量手段包括在機(jī)翼不同的展向位置安裝的固定龍頭和使用升華技術(shù)的表面剪應(yīng)力流動顯示技術(shù)。

驗(yàn)證算例的計算條件：迎角α=2°，雷諾數(shù)Re=1.5×106，馬赫數(shù)Ma=0.82，自由來流湍流度FSTI=0.5%，湍流黏度與層流黏性之比μt/μ=10。

DLR F-5的幾何模型是一個帶有20°后掠角的后掠翼，平均氣動弦長為0.15 m，機(jī)翼剖面為對稱翼型，且在Ma=0.82時為超臨界狀態(tài)。三維機(jī)翼模型和其表面計算網(wǎng)格則如圖1中所示，第1層網(wǎng)格高度L=1.5×10-7m，滿足y+=1，附面層層數(shù)為41，總網(wǎng)格量約為225萬。

圖1 DLR F-5計算網(wǎng)格及壓力分布提取截面

2.2 計算結(jié)果與討論

采用自由轉(zhuǎn)捩預(yù)測方法、全湍流方法計算所得沿機(jī)翼展向各個剖面的壓力分布和試驗(yàn)結(jié)果的對比如圖2所示。其中縱坐標(biāo)Cp為壓力系數(shù)，橫坐標(biāo)x/c為單位化的機(jī)翼截面弦向距離，各截面位置已經(jīng)在圖1中給出。Exp-Upper表示機(jī)翼剖面上表面，Exp-Lowwer表示機(jī)翼剖面下表面，SST為全湍計算結(jié)果，γ-Reθ為轉(zhuǎn)捩模型計算結(jié)果。在弦向截面為0.11%、20.47%、64.58%處，γ-Reθ轉(zhuǎn)捩模型對激波的捕捉能力略差，預(yù)測到的激波強(qiáng)度偏強(qiáng)，激波位置，壓力分布計算結(jié)果與試驗(yàn)存在一定差異，但前緣10%至激波前的順壓梯度段計算結(jié)果整體與試驗(yàn)值吻合良好。

圖2 計算所得到的機(jī)翼展向各剖面壓力分布形態(tài)與試驗(yàn)值的對比

由轉(zhuǎn)捩模型預(yù)測得到的表面摩擦系數(shù)Cf等值線和表面流線如圖3(a)所示，全湍流計算結(jié)果如圖3(b)所示?；诹鲃涌梢暬蛪毫y量技術(shù)，構(gòu)建出了機(jī)翼周圍流場的簡圖，如圖3(c)所示。測量結(jié)果表明，邊界層層流區(qū)從前緣維持到60%弦長處，由于60%弦長附近激波產(chǎn)生的強(qiáng)逆壓梯度影響，層流邊界層轉(zhuǎn)捩并分離。通過分離區(qū)后，流動再附著，但由于流動穩(wěn)定性已經(jīng)被破壞，邊界層變?yōu)橥牧鬟吔鐚?。從表面摩擦系?shù)云圖中可以清晰地看到層流分離和湍流再附的位置，與試驗(yàn)所得到的簡圖上的位置相吻合，均是從Y/Span=20%處開始直至翼尖。

圖3 F5機(jī)翼表面摩擦系數(shù)云圖和表面流線

綜上所述，雖然激波預(yù)測結(jié)果與試驗(yàn)結(jié)果一定差異，但本文任務(wù)為研究飛機(jī)的層流特性。在跨聲速狀態(tài)下，所采用的轉(zhuǎn)捩預(yù)測方法能夠給出合理的轉(zhuǎn)捩發(fā)生位置，足以為研究任務(wù)提供支持。所采用的轉(zhuǎn)捩預(yù)測方法可以較好地模擬表面邊界層的發(fā)展，表面速度型和摩擦阻力系數(shù)的分布情況與試驗(yàn)吻合良好，為之后的層流機(jī)翼跨聲速層流特性驗(yàn)證及轉(zhuǎn)捩因素敏感性分析提供了很可靠的數(shù)據(jù)保證。

3 某層流機(jī)翼驗(yàn)證機(jī)跨聲速層流特性驗(yàn)證

利用某傳統(tǒng)翼型來驗(yàn)證高速條件下層流驗(yàn)證機(jī)中央翼段數(shù)值模擬的準(zhǔn)確性。試驗(yàn)采用的風(fēng)洞為風(fēng)雷FL-2號風(fēng)洞。試驗(yàn)采用1∶9.8縮比模型，馬赫數(shù)Ma=0.7，雷諾數(shù)Re=8×106，側(cè)滑角β=0°，自由來流湍流度FSTI=0.6%，湍流黏性與層流黏性之比μt/μ=10。圖4所示為風(fēng)洞試驗(yàn)圖。數(shù)值模擬計算條件與試驗(yàn)條件相同，數(shù)模及網(wǎng)格拓?fù)淙鐖D5所示，第1層網(wǎng)格高度L=7.5×10-7m，滿足y+=1，附面層層數(shù)為41，總網(wǎng)格量約為516萬。

圖4 層流驗(yàn)證機(jī)高速PSP試驗(yàn)

圖5 計算模型、拓?fù)浜途W(wǎng)格

圖6為迎角α=2°時風(fēng)洞試驗(yàn)結(jié)果與計算結(jié)果對比。從圖6(a)的試驗(yàn)熱敏成像圖可以看出，中央翼段的轉(zhuǎn)捩區(qū)呈現(xiàn)出鋸齒狀特征，與圖6(b)計算結(jié)果所呈現(xiàn)的光滑轉(zhuǎn)捩區(qū)存在一定差異，該鋸齒狀形態(tài)是由于用于試驗(yàn)的模型翼段表面存在一定程度的粗糙度，導(dǎo)致局部轉(zhuǎn)捩提前，因而試驗(yàn)轉(zhuǎn)捩位置取層流區(qū)所能達(dá)到最大弦向距離。

圖6 α=2°時試驗(yàn)與計算結(jié)果對比

表1為不同迎角下計算與試驗(yàn)轉(zhuǎn)捩位置的對比。從對比數(shù)據(jù)可以看出，在高速條件下，計算所得的轉(zhuǎn)捩位置與試驗(yàn)結(jié)果存在一定差異，而隨著迎角的增加，計算誤差逐漸減小，而且計算與試驗(yàn)的轉(zhuǎn)捩位置隨迎角的變化趨勢一致，α=2°時誤差僅有2%。考慮到后續(xù)對機(jī)翼進(jìn)行轉(zhuǎn)捩敏感性分析時，重點(diǎn)關(guān)注在不同條件下飛機(jī)的轉(zhuǎn)捩變化趨勢，對計算結(jié)果的絕對數(shù)值不做太高要求，計算結(jié)果定性分析準(zhǔn)確，能準(zhǔn)確反映轉(zhuǎn)捩現(xiàn)象的變化趨勢即可。

表1 不同迎角的計算與試驗(yàn)轉(zhuǎn)捩位置對比

圖7所示為α=2°時上翼面的試驗(yàn)PSP (Pressure Sensitive Paint)壓力系數(shù)云圖，根據(jù)PSP云圖沿展向做空間平均得出不同迎角下試驗(yàn)的Cp曲線。圖8所示為α=2°時Cp試驗(yàn)曲線與計算曲線對比，從圖中可以看出，在前緣附近Cp計算值比試驗(yàn)值略低，而從弦向位置約16%c處開始，Cp計算值逐漸高于試驗(yàn)值，分析試驗(yàn)與計算數(shù)據(jù)差異的來源為試驗(yàn)?zāi)Ｐ图庸r的翼型厚度存在的一定差異。

圖7 α=2°時試驗(yàn)PSP壓力云圖

圖8 α=2°時壓力系數(shù)試驗(yàn)值與計算值對比

從Cp分布趨勢來看，計算曲線與試驗(yàn)曲線幾乎保持一致，包括壓力峰值與翼面上壓力拐點(diǎn)的弦向站位均吻合良好。故采用的計算方法對于層流機(jī)翼壓力分布預(yù)測較為準(zhǔn)確。

本文主要研究層流驗(yàn)證機(jī)的中段翼，如圖5所示，該段機(jī)翼后掠角β=0°，且機(jī)翼兩側(cè)有機(jī)身，有利于抑制翼尖流動的三維效應(yīng)，故整個驗(yàn)證機(jī)的中段翼流動接近于二維翼型模型，這也正是本驗(yàn)證機(jī)的主要設(shè)計思路。如圖9所示，在本文雷諾數(shù)Re=107附近和β=0°情況下的轉(zhuǎn)捩模式主要為T-S波失穩(wěn)(Tollmien-Schlichting Instability)主導(dǎo)轉(zhuǎn)捩，未出現(xiàn)橫流不穩(wěn)定性轉(zhuǎn)捩。試驗(yàn)風(fēng)洞為低湍流度風(fēng)洞，故不考慮橫流行波失穩(wěn)。對于橫流駐波主導(dǎo)的失穩(wěn)，其擾動源主要是壁面粗糙度。本次層流驗(yàn)證機(jī)風(fēng)洞層流試驗(yàn)?zāi)Ｐ捅绕胀y力測壓的試驗(yàn)?zāi)Ｐ途哂懈叩墓鉂嵍?，故具有比較高的橫流駐波失穩(wěn)臨界雷諾數(shù)，不容易發(fā)生該類型的轉(zhuǎn)捩。因此，分析不考慮橫流不穩(wěn)定性轉(zhuǎn)捩。綜上所述，采用第1節(jié)提到的計算方法足夠反映研究的問題，后續(xù)將據(jù)此方法展開詳細(xì)計算并對計算結(jié)果進(jìn)行系統(tǒng)分析。

圖9 失穩(wěn)類型與前緣后掠角、雷諾數(shù)之間的關(guān)系

4 計算模型與網(wǎng)格

本文研究對象為某層流機(jī)翼驗(yàn)證機(jī)，計算模型的三視圖如圖10所示。模型采用雙發(fā)、雙機(jī)身、雙斜置垂尾的布局形式，兩側(cè)的半機(jī)身通過中央層流機(jī)翼驗(yàn)證段相連接，驗(yàn)證翼段設(shè)置8°的前緣后掠角。該機(jī)雙機(jī)身的設(shè)計，可以保證層流驗(yàn)證機(jī)在飛行時飛機(jī)中段翼盡可能減少橫流轉(zhuǎn)捩，排除無關(guān)變量的干擾。同時，中央層流機(jī)翼驗(yàn)證段可以進(jìn)行更換，對于層流翼型在飛行雷諾數(shù)下的研究，經(jīng)濟(jì)高效。計算模型主要幾何參數(shù)為：參考面積為5 m2，參考展長為 5.9 m，機(jī)身長度為 4.4 m，平均氣動弦長為 1.5 m，質(zhì)心距離機(jī)頭的距離為0.5 m。

圖10 計算模型三視圖

流場網(wǎng)格拓?fù)浜万?yàn)證機(jī)表面網(wǎng)格情況分別如圖11、圖12所示。圖12網(wǎng)格量為310萬，對圖12的網(wǎng)格進(jìn)行整體的網(wǎng)格量放縮，最終得到6套網(wǎng)格。網(wǎng)格量從疏到密分別為310萬、640萬、1 300萬、2 700萬、5 600萬、12 000萬，對以上6套網(wǎng)格進(jìn)行編號，依次命名為1號、2號、3號、4號、5號、6號。采用表2邊界條件對6套網(wǎng)格共12個計算狀態(tài)進(jìn)行CFD(Computational Fluid Dynamics)計算。

表2 網(wǎng)格無關(guān)性驗(yàn)證邊界條件

圖11 流場網(wǎng)格拓?fù)?/p>

圖12 計算網(wǎng)格

圖13為網(wǎng)格無關(guān)性驗(yàn)證氣動力系數(shù)曲線。圖13(a)為α=2°時的升力系數(shù)曲線，升力系數(shù)隨著網(wǎng)格量的增大而減小，圖中基于SST計算模型的升力系數(shù)最大差量為8×10-4，平均相對誤差為0.10%。基于γ-Reθ計算模型的升力系數(shù)最大差量為8.1×10-4，平均相對誤差為0.11%。圖13(b) 為α=2°時的阻力系數(shù)曲線，阻力系數(shù)隨著網(wǎng)格量的增大而增大，圖中基于SST計算模型的阻力系數(shù)最大差量為9.7×10-5，平均相對誤差為0.05%?；讦?Reθ計算模型的阻力系數(shù)最大差量為6.5×10-5，平均相對誤差為0.03%。圖13(c) 為2°迎角時的力矩系數(shù)曲線，力矩系數(shù)隨著網(wǎng)格量的增大而增大，圖中基于SST計算模型的力矩系數(shù)最大差量為9.7×10-5，平均相對誤差為0.06%。基于γ-Reθ計算模型的力矩系數(shù)最大差量為6.5×10-5，平均相對誤差為0.07%。對比情況見表3。

表3 氣動力系數(shù)隨網(wǎng)格量變化對比

圖13 網(wǎng)格無關(guān)性驗(yàn)證氣動力系數(shù)曲線(迎角為2°)

綜上所述，網(wǎng)格量變化引起的相對誤差量級在10-3～10-4，網(wǎng)格量變化對于計算結(jié)果影響較小。選擇3號網(wǎng)格作為驗(yàn)證機(jī)計算網(wǎng)格。3號網(wǎng)格第1層網(wǎng)格高度L=1.5×10-6m。附面層層數(shù)為41，附面層沿物面法向向外的增長比為1.14，總網(wǎng)格量1 300萬左右，第1層網(wǎng)格間距滿足y+=1。

5 計算結(jié)果與討論

5.1 自由轉(zhuǎn)捩與全湍流計算結(jié)果對比

計算狀態(tài)見表4中所列共計14個計算狀態(tài)。分別采用自由轉(zhuǎn)捩模式和全湍流計算模式對飛行驗(yàn)證機(jī)巡航馬赫數(shù)下的迎角序列進(jìn)行計算，得到全機(jī)的升力系數(shù)曲線、極曲線和力矩系數(shù)曲線如圖14所示。

表4 全機(jī)氣動特性計算狀態(tài)和方法

從圖14(a)～圖14(c)可以看出，采用自由轉(zhuǎn)捩模式或全湍流模式對于全機(jī)的升力和力矩特性影響不大，對阻力特性會產(chǎn)生明顯影響，但影響的幅度不大，因?yàn)樽杂赊D(zhuǎn)捩的層流區(qū)主要存在于中央翼段的部分區(qū)域，影響范圍相對有限。阻力差異主要集中在高速狀態(tài)的巡航升力系數(shù)附近(見圖14(c))。進(jìn)一步對阻力系數(shù)差異進(jìn)行分析，以自由轉(zhuǎn)捩方法計算結(jié)果的中段翼壓力分布和摩阻云圖為研究對象，以全湍流方法計算結(jié)果的摩阻云圖為對照，Ma=0.7時，如圖14(d)所示，機(jī)翼表面未出現(xiàn)激波，轉(zhuǎn)捩發(fā)生于弱逆壓梯度起始點(diǎn)，轉(zhuǎn)捩方式為T-S波突然失穩(wěn)形成的自然轉(zhuǎn)捩。升力系數(shù)較小時如CL=0.12，順壓梯度區(qū)較大，占中段翼弦長約50%左右，自由轉(zhuǎn)捩其層流區(qū)較大，自由轉(zhuǎn)捩摩擦阻力較全湍流模式要小，差異明顯。而升力系數(shù)較大時如CL=0.60，機(jī)翼上表面順壓梯度區(qū)較小，占中段翼弦長約10%左右，圖14(d)機(jī)翼上表面逆壓梯度區(qū)較大，轉(zhuǎn)捩位置較為靠近前緣，層流區(qū)很小，自由轉(zhuǎn)捩計算結(jié)果的摩擦阻力接近全湍流計算結(jié)果，差異很小。

圖14 全機(jī)氣動力系數(shù)曲線

5.2 自由來流湍流度的影響

圖15為壓力分布展向截面示意圖，在馬赫數(shù)Ma=0.7，雷諾數(shù)Re=1.1×107，迎角α=2°，自由來流湍流度序列為FSTI=0.2%,0.6%,1.0%情形下的自由轉(zhuǎn)捩計算結(jié)果見圖16～圖18。

圖15 中央翼段展向站位

從圖16的壓力分布計算結(jié)果可以看出，自由來流湍流度的變化對中段翼層流壓力分布特征基本無影響；相反，從圖17、圖18的計算結(jié)果中可以看出，上表面層流區(qū)長度受自由來流湍流度的影響較大。自由來流湍流度為0.2%時，轉(zhuǎn)捩位置在逆壓梯度起始點(diǎn)附近(距中段翼前緣50%c位置處)，自由來流湍流度為0.6%時，轉(zhuǎn)捩發(fā)生在弱順壓梯度區(qū)(距中段翼前緣約30%c位置處)，自由來流湍流度為1%時，轉(zhuǎn)捩發(fā)生在強(qiáng)順壓梯度截止點(diǎn)(距中段翼前緣約18%c位置處)，轉(zhuǎn)捩形式均為T-S波急劇失穩(wěn)的自然轉(zhuǎn)捩。故可知順壓梯度對于層流流動的維持并非充分條件，湍流度的增加直接導(dǎo)致展向各站位處轉(zhuǎn)捩位置的大幅提前，轉(zhuǎn)捩位置的變化特征表現(xiàn)出一定的規(guī)律性：自由來流中的湍流流動對機(jī)翼層流區(qū)的維持具有較強(qiáng)的破壞作用，湍流流動越強(qiáng)，越易導(dǎo)致T-S波失穩(wěn)。

圖16 不同自由來流湍流度下中央翼段展向各站位壓力分布系數(shù)曲線對比

圖17 不同自由來流湍流度下全機(jī)表面摩擦阻力系數(shù)分布

圖18 不同自由來流湍流度下中央翼段上表面展向各站位轉(zhuǎn)捩位置

5.3 雷諾數(shù)的影響

馬赫數(shù)Ma=0.70，自由來流湍流度FSTI=0.2%，迎角α=2°，雷諾數(shù)Re=5.50×106, 8.25×106, 1.10×107, 1.35×107, 1.65×107, 2.00×107情形下的自由轉(zhuǎn)捩計算結(jié)果如圖19～圖21所示。

對于一般的附著流動來說，雷諾數(shù)的大小影響模型表面上附面層的性質(zhì)，從而改變附面層的厚度、附面層轉(zhuǎn)捩位置、表面摩擦阻力,以及與氣體黏性有關(guān)的氣流分離情況。從圖19中可以看出，Re>8.25×106之后計算所得到的壓力分布基本都是一致的，Re=5.5×106下的壓力分布特征在60%弦長處表現(xiàn)為輕微的差異，在雷諾數(shù)的一定區(qū)間內(nèi)(Re=107左右)，雷諾數(shù)對于壓力分布的影響很小。從圖19、圖20的計算結(jié)果中可以看出，雷諾數(shù)的變化對上表面層流區(qū)長度的影響具有明顯的非線性特性。當(dāng)雷諾數(shù)介于5.50×106～1.35×107時，雷諾數(shù)的增加直接導(dǎo)致展向各站位處轉(zhuǎn)捩位置有規(guī)律的明顯前移，隨著雷諾數(shù)的增加轉(zhuǎn)捩位置前移的幅度也會先有所增大，轉(zhuǎn)捩發(fā)生在逆壓梯度起始點(diǎn)和弱順壓梯度區(qū)。這主要是因?yàn)殡S著雷諾數(shù)增大，機(jī)翼附面層變薄，流動穩(wěn)定性降低，機(jī)翼表面附近的流動遇到擾動更容易失穩(wěn)變成湍流，轉(zhuǎn)捩位置逐漸向強(qiáng)順壓梯度區(qū)靠近。轉(zhuǎn)捩方式為T-S波快速失穩(wěn)的自然轉(zhuǎn)捩。從圖21可以看到，在Re>1.35×107之后，變化相同程度的雷諾數(shù)，中段翼上表面轉(zhuǎn)捩位置前移程度很小，移動量為%c的個位數(shù)。這主要是此時轉(zhuǎn)捩位置已逐步接近中段翼上表面的強(qiáng)順壓梯度區(qū)，即使由于雷諾數(shù)增大附面層已經(jīng)變薄，流動穩(wěn)定性有降低趨勢,但在較強(qiáng)的順壓梯度影響下，附面層流動獲得了足夠維持層流穩(wěn)定的能量，熵增程度有限，主導(dǎo)轉(zhuǎn)捩的因素由雷諾數(shù)變?yōu)閺?qiáng)順壓梯度，不易發(fā)生轉(zhuǎn)捩，由雷諾數(shù)影響的轉(zhuǎn)捩位置提前程度逐漸達(dá)到極限，此時轉(zhuǎn)捩方式依舊為T-S波急劇失穩(wěn)導(dǎo)致的自然轉(zhuǎn)捩。雷諾數(shù)對轉(zhuǎn)捩位置的影響容易被其他敏感性因素所覆蓋，敏感性優(yōu)先級較低。

圖19 不同雷諾數(shù)下中央翼段展向各站位壓力分布系數(shù)曲線對比

圖20 不同雷諾數(shù)下全機(jī)表面摩擦阻力系數(shù)分布

圖21 不同雷諾數(shù)下中央翼段上表面展向各站位轉(zhuǎn)捩位置

5.4 馬赫數(shù)的影響

雷諾數(shù)Re=1.10×107，迎角α=0°，自由來流湍度FSTI=0.2%，馬赫數(shù)Ma=0.61,0.64,0.67,0.70,0.73,0.76,0.79情形下的自由轉(zhuǎn)捩計算結(jié)果如圖22～圖24所示。

從圖22的計算結(jié)果中可以看出，馬赫數(shù)的變化對中央層流翼段壓力分布特征有非常明顯的影響，隨著馬赫數(shù)增大，除駐點(diǎn)和后緣外，上下翼面壓力分布負(fù)壓值整體上移。當(dāng)Ma=0.76,0.79時，翼面弱激波強(qiáng)度大幅增加，激波前后壓力分布形態(tài)有顯著的變化。

圖22 不同馬赫數(shù)下中央翼段展向各站位壓力分布系數(shù)曲線對比

從圖23、圖24的整個馬赫數(shù)序列的計算結(jié)果中可以發(fā)現(xiàn)，當(dāng)Ma=0.70時，翼面上表面層流區(qū)長度最長，轉(zhuǎn)捩位置最為靠后；當(dāng)馬赫數(shù)減小時，翼面轉(zhuǎn)捩位置也隨之提前，馬赫數(shù)由0.70降低至0.67時層流區(qū)長度有輕微的減短，但隨著馬赫數(shù)繼續(xù)降低為0.64時，翼面的轉(zhuǎn)捩位置提前到將近20%c左右，層流區(qū)長度只有原先的1/2；馬赫數(shù)繼續(xù)降低，轉(zhuǎn)捩位置表現(xiàn)為小范圍的波動。而當(dāng)馬赫數(shù)從0.70增加時，翼面層流區(qū)長度和轉(zhuǎn)捩位置的變化均相對較小，只是在小范圍內(nèi)波動。轉(zhuǎn)捩位置突變的馬赫數(shù)Ma=0.64，當(dāng)Ma>0.67以后，轉(zhuǎn)捩位置變化較小。巡航馬赫數(shù)Ma=0.70，故轉(zhuǎn)捩位置的不穩(wěn)定變化對實(shí)際飛行影響較小，巡航馬赫數(shù)的層流區(qū)在高速區(qū)間內(nèi)最長，證明本文層流機(jī)翼的設(shè)計較為優(yōu)秀。

圖23 不同馬赫數(shù)下全機(jī)表面摩擦阻力系數(shù)分布

如圖24所示，馬赫數(shù)對層流區(qū)的影響具有較強(qiáng)的非線性特性。隨著馬赫數(shù)增大，轉(zhuǎn)捩位置先逐漸向后緣靠近，當(dāng)轉(zhuǎn)捩位置到達(dá)40%c后，在40%c近波動。隨著馬赫數(shù)從0.6增大到0.7，逆壓梯度起始點(diǎn)略有后移，轉(zhuǎn)捩位置相應(yīng)向后移動,轉(zhuǎn)捩形式為T-S波急劇失穩(wěn)產(chǎn)生的自然轉(zhuǎn)捩。Ma=0.70 時，逆壓梯度起始點(diǎn)附近壓力分布形式接近于壓力平臺，轉(zhuǎn)捩發(fā)生于壓力分布向逆壓梯度發(fā)展的過渡階段，T-S波逐漸失穩(wěn)形成的自然轉(zhuǎn)捩。隨著馬赫數(shù)從0.70增大到0.79，層流驗(yàn)證機(jī)中段翼表面負(fù)壓峰值逐漸升高，激波開始形成，在30%c～70%c壓力分布形態(tài)變得不規(guī)則，順壓梯度和逆壓梯度交替出現(xiàn)，在馬赫數(shù)增大導(dǎo)致的當(dāng)?shù)乩字Z數(shù)變化和弱激波的相互干擾下，層流驗(yàn)證機(jī)中段翼上表面轉(zhuǎn)捩位置變化沒有明顯的規(guī)律，呈現(xiàn)高度非線性。由于壓力分布梯度不規(guī)則，出現(xiàn)多個局部小范圍逆壓梯度，導(dǎo)致整體轉(zhuǎn)捩位置要比Ma=0.70時要靠前。Ma=0.79時，轉(zhuǎn)捩均發(fā)生在弱順壓梯度區(qū)，而非激波之后的強(qiáng)逆壓梯度區(qū)，轉(zhuǎn)捩形式依然為自然轉(zhuǎn)捩。另一方面也反映出層流區(qū)對于逆壓梯度非常敏感，局部小范圍的逆壓梯度很可能導(dǎo)致T-S波失穩(wěn)而發(fā)生轉(zhuǎn)捩。

圖24 不同馬赫數(shù)下中央翼段上表面展向各站位轉(zhuǎn)捩 位置

5.5 來流迎角的影響

馬赫數(shù)Ma=0.70，雷諾數(shù)Re=1.1×107，自由來流湍流度FSTI=0.2%，來流迎角α=0°,2°,4°情形下,自由轉(zhuǎn)捩計算結(jié)果如圖25～圖27所示。

由圖25～圖27的計算結(jié)果可知，迎角的變化對中央層流翼段壓力分布形態(tài)特征和上表面層流區(qū)長度的影響非常明顯，且具有很強(qiáng)的規(guī)律性。隨著迎角增大，壓力分布所圍成的面積逐漸增大，局部升力增大。然而層流驗(yàn)證機(jī)中段翼上表面從前緣駐點(diǎn)到60%c處的順壓梯度整體變小，上表面轉(zhuǎn)捩位置逐漸提前。隨著迎角從0°增大到4°，轉(zhuǎn)捩位置約從57%c處提前到30%，流動迎角的改變直接影響著壓力分布曲線前部的順壓梯度形態(tài)，從而進(jìn)一步影響到表面層流特性的維持。從圖25的15%、25%展向站位的壓力分布可以看出，隨著迎角增大，到4°迎角時，壓力分布出現(xiàn)了2個逆壓梯度區(qū)，而轉(zhuǎn)捩均發(fā)生于第1個逆壓梯度起始點(diǎn)附近。再次說明層流區(qū)對于逆壓梯度非常敏感，局部小范圍的逆壓梯度很可能導(dǎo)致T-S波失穩(wěn)而發(fā)生轉(zhuǎn)捩。降低迎角對于拓展層流區(qū)長度是有利的影響，同時會降低升力。迎角的增加會直接導(dǎo)致展向各站位處轉(zhuǎn)捩位置提前，升力增大。飛行器設(shè)計者需要做好兩者之間的權(quán)衡，在滿足升力需求的情況下，盡可能拓展機(jī)翼表面層流區(qū)的長度以求降低阻力。

圖25 不同來流迎角下中央翼段展向各站位壓力分布系數(shù)曲線對比

圖26 不同來流迎角下全機(jī)表面摩擦阻力系數(shù)分布

圖27 不同迎角下中央翼段上表面展向各站位轉(zhuǎn)捩位置

總之在一定范圍內(nèi)，跨聲速層流轉(zhuǎn)捩由多種因素共同制約且相互影響而造成。圖28為各敏感性因素影響轉(zhuǎn)捩的過程和相互關(guān)系。其中箭頭方向代表所影響的對象，線性代表影響權(quán)重。轉(zhuǎn)捩的本質(zhì)是流動熵增導(dǎo)致流動穩(wěn)定性被破壞。自由來流湍流度影響全局流動穩(wěn)定性，翼面壓力梯度影響局部流動穩(wěn)定性，兩者均直接影響流動穩(wěn)定性。雷諾數(shù)、馬赫數(shù)、迎角間接影響流動穩(wěn)定性，其中雷諾數(shù)容易被其他因素所覆蓋。

圖28 各敏感性因素影響轉(zhuǎn)捩的過程及相互關(guān)系

6 結(jié) 論

針對某特殊布局形式的層流機(jī)翼驗(yàn)證機(jī)開展了跨聲速層流特性和參數(shù)敏感性分析，重點(diǎn)關(guān)注全機(jī)巡航狀態(tài)附近的氣動特性和中央驗(yàn)證段在不同飛行狀態(tài)下的表面轉(zhuǎn)捩位置及層流區(qū)長度。結(jié)果表明：馬赫數(shù)、雷諾數(shù)、自由來流湍流度、來流迎角等均會對中央驗(yàn)證段表面轉(zhuǎn)捩位置產(chǎn)生明顯的影響。通過以上參數(shù)敏感性的分析，對于該層流機(jī)翼驗(yàn)證機(jī)，可以得出以下結(jié)論：

1) 自由來流湍流度、雷諾數(shù)、馬赫數(shù)、來流迎角對上表面層流區(qū)長度的影響都非常明顯，并且各自表現(xiàn)出的規(guī)律不同。

2) 自由來流湍流度、雷諾數(shù)、來流迎角對于轉(zhuǎn)捩位置的影響均具有明顯的規(guī)律性，轉(zhuǎn)捩位置隨自由來流湍流度、雷諾數(shù)、來流迎角對于轉(zhuǎn)捩位置的增加而提前。

3) 自由來流湍流度對于轉(zhuǎn)捩位置具有較強(qiáng)的影響。雷諾數(shù)改變轉(zhuǎn)捩位置受順壓梯度影響比較大。迎角對轉(zhuǎn)捩位置的影響主要通過改變壓力分布形態(tài)而實(shí)現(xiàn)。

4) 馬赫數(shù)對于層流驗(yàn)證機(jī)轉(zhuǎn)捩位置的影響呈現(xiàn)高度非線性。馬赫數(shù)增大使得機(jī)翼上表面出現(xiàn)激波，壓力分布形態(tài)不規(guī)則，產(chǎn)生的局部小范圍逆壓梯度會導(dǎo)致轉(zhuǎn)捩發(fā)生。

5) 總之在一定范圍內(nèi)，跨聲速層流轉(zhuǎn)捩由多種因素共同制約且相互影響而造成。計算模型的選擇將決定是否能夠捕捉到轉(zhuǎn)捩現(xiàn)象。轉(zhuǎn)捩的本質(zhì)是流動熵增導(dǎo)致流動穩(wěn)定性被破壞。自由來流湍流度影響全局流動穩(wěn)定性，翼面壓力梯度影響局部流動穩(wěn)定性，兩者均直接影響流動穩(wěn)定性，且影響程度都很大。雷諾數(shù)通過影響附面層性質(zhì)間接影響流動穩(wěn)定性，本身不直接影響流動穩(wěn)定性，且影響容易被其他因素所覆蓋。