基于雙eN方法的翼身組合體流動轉(zhuǎn)捩自動判斷

朱震，宋文萍，韓忠華

西北工業(yè)大學(xué) 航空學(xué)院 翼型葉柵空氣動力學(xué)國家級重點(diǎn)實(shí)驗(yàn)室，西安 710072

出于環(huán)保性和經(jīng)濟(jì)性的要求，通過減阻設(shè)計(jì)提升民機(jī)性能是設(shè)計(jì)師們長期追求的目標(biāo)。美國國家航空航天局(NASA)的ERA(Environmentally Responsible Aviation)項(xiàng)目[1]和歐洲航空研究咨詢委員會(Advisory Council for Aeronautical Research in Europe，ACARE)發(fā)布的《歐洲航空2050展望》[2]均指出，進(jìn)一步減少目前客機(jī)的阻力，減少二氧化碳排放和燃油的消耗是航空技術(shù)領(lǐng)域關(guān)注的熱點(diǎn)研究方向，對于提高新型客機(jī)性能具有重要的意義。目前，如果想要在減阻設(shè)計(jì)方面取得較大突破，自然層流機(jī)翼技術(shù)是一個(gè)值得期待的研究方向[3]。經(jīng)過多年發(fā)展，抑制Tollmien-Schlichting(T-S)波不穩(wěn)定擾動的自然層流翼型技術(shù)已經(jīng)較為成熟，但民機(jī)廣泛采用的后掠機(jī)翼設(shè)計(jì)給實(shí)現(xiàn)大范圍自然層流減阻設(shè)計(jì)帶來了挑戰(zhàn)。民機(jī)一般采用大展弦比后掠機(jī)翼設(shè)計(jì)推遲激波產(chǎn)生，減小激波阻力，增大巡航速度以提高巡航效率因子。但后掠機(jī)翼的邊界層由于受到橫向壓力梯度的影響，速度型存在橫流(Cross-Flow，CF)分量，橫流速度型具有拐點(diǎn)，很容易出現(xiàn)橫流不穩(wěn)定性(CFI)[4]導(dǎo)致的轉(zhuǎn)捩。而機(jī)身的三維位移效應(yīng)會進(jìn)一步增強(qiáng)機(jī)翼表面橫流流動，可能導(dǎo)致原本設(shè)計(jì)的具有較大范圍自然層流的后掠機(jī)翼因機(jī)身影響提前發(fā)生轉(zhuǎn)捩。因此面向翼身組合體外形，發(fā)展可靠實(shí)用的轉(zhuǎn)捩自動判斷方法，是當(dāng)前發(fā)展自然層流減阻設(shè)計(jì)技術(shù)的急需；通過翼身組合體流動轉(zhuǎn)捩判斷研究，進(jìn)一步分析機(jī)身對機(jī)翼轉(zhuǎn)捩的影響，對于自然層流后掠機(jī)翼的設(shè)計(jì)具有指導(dǎo)性意義。

目前國外研究機(jī)構(gòu)已經(jīng)圍繞翼型、機(jī)翼邊界層流動轉(zhuǎn)捩判斷開展了大量的研究工作[5-20]。而針對翼身組合體復(fù)雜三維外形，國外開展的邊界層轉(zhuǎn)捩判斷方法研究相對較少，且主要使用簡化的eN數(shù)據(jù)庫轉(zhuǎn)捩判斷方法。如法國宇航院(ONERA)在其雷諾平均Navier-Stokes(RANS)方程求解器elsA中加入了Arnal[9]和Casalis[12]等發(fā)展的eN數(shù)據(jù)庫方法，分別用來計(jì)算T-S波和CF波誘導(dǎo)的轉(zhuǎn)捩，并開展了帶增升裝置的翼身組合體外形的轉(zhuǎn)捩判斷研究[21]。德國宇航院(DLR)在其結(jié)構(gòu)化網(wǎng)格RANS求解器FLOWer[17]和非結(jié)構(gòu)/混合網(wǎng)格RANS求解器TAU[22]中加入了eN數(shù)據(jù)庫方法，開展了三維復(fù)雜外形的轉(zhuǎn)捩判斷研究。美國斯坦福大學(xué)的Lee和Jameson[23]開展了基于耦合eN數(shù)據(jù)庫方法的RANS求解器的層流機(jī)翼設(shè)計(jì)工作。eN數(shù)據(jù)庫方法是完全eN方法的近似和簡化。eN數(shù)據(jù)庫方法并不求解線性穩(wěn)定性方程，而只是根據(jù)邊界層的形狀因子、速度型拐點(diǎn)信息等部分特征參數(shù)，在數(shù)據(jù)庫中查詢事先求解的平板邊界層不穩(wěn)定擾動累積放大因子的包絡(luò)線。而完全eN方法則是通過求解線性穩(wěn)定性方程，得到真實(shí)機(jī)翼邊界層的不同頻率或波數(shù)的不穩(wěn)定擾動波的累積放大因子發(fā)展規(guī)律。最近，以γ-Reθ t轉(zhuǎn)捩模型為基礎(chǔ)發(fā)展橫流轉(zhuǎn)捩判斷方法的研究也逐漸開展起來。2016年，Grabe等[24]使用引入C1橫流判據(jù)和He橫流判據(jù)的γ-Reθ t轉(zhuǎn)捩模型開展了DLR-F4翼身組合體的機(jī)翼邊界層轉(zhuǎn)捩判斷研究。

目前，國內(nèi)也鮮有使用完全eN方法對翼身組合體復(fù)雜三維外形開展轉(zhuǎn)捩判斷的研究工作，國內(nèi)已開展的轉(zhuǎn)捩判斷研究工作主要針對翼型、機(jī)翼繞流進(jìn)行。在使用穩(wěn)定性理論和eN方法進(jìn)行的后掠機(jī)翼邊界層穩(wěn)定性及轉(zhuǎn)捩判斷研究方面，張坤和宋文萍[25]通過線性穩(wěn)定性理論(Linear Stability Theory，LST)，左歲寒等[26]通過線性拋物化穩(wěn)定性方程(Linear Parabolized Stability Equations，LPSE)開展了機(jī)翼邊界層轉(zhuǎn)捩判斷研究。孫朋朋[27]、靖振榮[28]和黃章峰[29]等通過LST開展了一系列的后掠機(jī)翼邊界層穩(wěn)定性及轉(zhuǎn)捩研究，以及使用擴(kuò)展的O-S方程(Extended Orr-Sommerfeld equation，EOS)和LPSE進(jìn)行了機(jī)翼邊界層的橫流穩(wěn)定性分析和轉(zhuǎn)捩判斷研究。韓忠華等[30]開展了基于線性穩(wěn)定性理論的雙eN方法耦合RANS求解器進(jìn)行自然層流后掠機(jī)翼優(yōu)化設(shè)計(jì)的工作。同樣，國內(nèi)近期也開展了以γ-Reθ t轉(zhuǎn)捩模型為基礎(chǔ)，發(fā)展橫流轉(zhuǎn)捩判斷方法的研究，如徐家寬[31-32]、史亞云[33]、戚瓊[34]、鞠勝軍[35]等使用添加了橫流判據(jù)的γ-Reθ t轉(zhuǎn)捩模型開展了三維流動中的橫流轉(zhuǎn)捩判斷研究。

本文使用多塊結(jié)構(gòu)化網(wǎng)格和RANS求解器，耦合邊界層方程求解和基于線性穩(wěn)定性理論的完全eN方法，發(fā)展了面向翼身組合體復(fù)雜三維外形流動的機(jī)翼轉(zhuǎn)捩自動判斷方法。采用雙eN方法，針對邊界層方程求解得到的邊界層速度型、溫度型等信息，通過求解線性穩(wěn)定性方程，分別得到機(jī)翼邊界層T-S波和CF波的擾動累積放大因子曲線，通過閥值判斷轉(zhuǎn)捩位置。該方法不僅可以判斷轉(zhuǎn)捩位置，還可以得到轉(zhuǎn)捩誘發(fā)機(jī)制，能夠?yàn)樽匀粚恿髟O(shè)計(jì)提供指導(dǎo)。針對翼身組合體外形繞流，為了便于實(shí)現(xiàn)轉(zhuǎn)捩自動判斷，本文采用一種通用的結(jié)構(gòu)化網(wǎng)格拓?fù)浣Y(jié)構(gòu)處理方法，對所有翼面類部件適用。使用上述方法，針對中短程民機(jī)，開展翼身組合體繞流轉(zhuǎn)捩判斷研究，并通過與單獨(dú)自然層流機(jī)翼的對比，進(jìn)一步研究了機(jī)身對后掠機(jī)翼邊界層橫流速度型、橫流不穩(wěn)定擾動放大率和轉(zhuǎn)捩位置的影響規(guī)律，可為考慮機(jī)身影響的自然層流后掠機(jī)翼減阻設(shè)計(jì)提供支持。

1 數(shù)值方法

1.1 完全雙eN轉(zhuǎn)捩判斷方法

對機(jī)翼邊界層的轉(zhuǎn)捩位置進(jìn)行準(zhǔn)確判斷是自然層流機(jī)翼減阻設(shè)計(jì)的關(guān)鍵技術(shù)之一。由層流到湍流的轉(zhuǎn)捩誘發(fā)機(jī)理復(fù)雜，影響因素眾多，轉(zhuǎn)捩判斷問題目前仍是流體力學(xué)中尚未完全解決的問題之一。基于線性穩(wěn)定性理論的eN方法是目前在工程實(shí)踐中得到最廣泛認(rèn)可的轉(zhuǎn)捩判斷方法[36]，因此本文基于線性穩(wěn)定性分析的完全eN方法判斷轉(zhuǎn)捩。針對自然層流翼身組合體外形，采用空間放大理論分析小擾動，進(jìn)行機(jī)翼邊界層轉(zhuǎn)捩自動判斷。eN方法由Smith[37]和Ingen[38]等首次提出，其思想是對邊界層流動中的小擾動進(jìn)行線性穩(wěn)定性分析，如果擾動逐漸衰減，則是穩(wěn)定的，如果擾動被放大，則是不穩(wěn)定的。對于不穩(wěn)定的擾動，從其開始放大處起，沿下游方向計(jì)算其累積的線性放大倍數(shù)，當(dāng)累積放大倍數(shù)到達(dá)擾動開始放大處振幅的eN倍時(shí)，認(rèn)為轉(zhuǎn)捩發(fā)生。基于平行流假設(shè)的三維小擾動形式為

(1)

(2)

圖1 空間放大理論中擾動波數(shù)矢量與擾動增長矢量 Fig.1 Disturbance wave number vector and disturbance amplification vector in spatial amplification theory

基于線性穩(wěn)定性理論的完全eN方法，對機(jī)翼邊界層信息求解線性穩(wěn)定性方程，進(jìn)行擾動累積放大因子的計(jì)算，而不對數(shù)據(jù)庫進(jìn)行查詢，因此不存在近似和簡化。完全eN方法進(jìn)行邊界層轉(zhuǎn)捩位置的判斷分兩步進(jìn)行：① 計(jì)算中性曲線，即擾動的穩(wěn)定區(qū)域和不穩(wěn)定區(qū)域的分界線。通過中性曲線就可以確定不穩(wěn)定擾動的頻率范圍、波數(shù)范圍和不同擾動各自開始被放大的位置；② 按照適當(dāng)?shù)姆e分策略，對這些頻率的擾動放大率從中性曲線下半支出發(fā)，沿著擾動傳播路徑進(jìn)行積分，得到擾動累積放大因子N。這兩步都需要求解三維可壓縮線性穩(wěn)定性方程，其數(shù)值求解方法為：① 采用連續(xù)法獲得具有較高精確度的特征變量初值；② 將高階的線性穩(wěn)定性控制方程組轉(zhuǎn)化為一階的控制方程組；③ 采用中心差分方法離散控制方程組；④ 采用牛頓法線化離散后的控制方程組；⑤ 采用塊矩陣消去法求解方程組。

1.1.1 固定波角的方法

對于T-S波不穩(wěn)定(TSI)擾動，原有的包絡(luò)線方法[39]始終尋找擾動放大率最大的波角方向進(jìn)行積分，如果存在較強(qiáng)的CF波，可能會錯(cuò)誤地把CF波擾動積分到T-S波擾動累積放大因子NTS中。本文對包絡(luò)線方法進(jìn)行了改進(jìn)，采用固定波角的方法計(jì)算T-S波擾動累積放大因子，即

(3)

式中：f為頻率；x0為中性曲線計(jì)算出的小擾動開始放大的x位置。本文采用固定波角方法，選擇擾動波波角(波數(shù)矢量與x軸的夾角)ψ= 0° 方向，即邊界層外邊界速度方向，計(jì)算不同頻率f擾動的放大率并進(jìn)行積分，因此能夠保證始終對T-S波擾動進(jìn)行積分，不會混入CF波擾動。

1.1.2 同時(shí)固定βr和f的方法

對于CF波不穩(wěn)定(CFI)擾動誘導(dǎo)的轉(zhuǎn)捩，本文采用同時(shí)固定βr和f的方法計(jì)算CF波擾動累積放大因子，即

(4)

式中：NCF為CF波擾動累積放大因子；βr為橫流展向波數(shù)；f=0 Hz時(shí)為橫流駐波。針對本文中短程民機(jī)巡航狀態(tài)，選擇對不同展向波數(shù)βr的橫流駐波進(jìn)行擾動放大率的積分，計(jì)算橫流擾動累積放大因子NCF。對于CF波不穩(wěn)定擾動，基于平行流假設(shè)的線性穩(wěn)定性理論沒有考慮擾動幅值本身的變化存在一定偏差。根據(jù)黃章峰等[30]的研究，線性穩(wěn)定性方程計(jì)算得到的擾動放大增長因子與考慮擾動幅值本身變化的拋物化穩(wěn)定性方程相比，形態(tài)趨勢一致但相對偏小。因此本文采用線性穩(wěn)定性方程仍可較為有效地預(yù)測CF波不穩(wěn)定誘導(dǎo)的轉(zhuǎn)捩，但其對應(yīng)的轉(zhuǎn)捩閥值相對較小。

1.1.3 雙eN轉(zhuǎn)捩判斷方法

雙eN方法是一種轉(zhuǎn)捩判斷方法的構(gòu)架，可以針對T-S波和CF波選擇不同的擾動放大率積分策略構(gòu)建不同的轉(zhuǎn)捩判斷方法。對于T-S波的擾動放大增長因子NTS，采用1.1.1節(jié)描述的固定波角的方法計(jì)算；對于CF波的擾動放大增長因子NCF，采用1.1.2節(jié)描述的同時(shí)固定βr和f的方法計(jì)算。根據(jù)事先給定的轉(zhuǎn)捩閥值[(NTS)tr,(NCF)tr]，確定機(jī)翼邊界層的轉(zhuǎn)捩位置。如果NTS先到達(dá)其對應(yīng)的閥值(NTS)tr，就認(rèn)為T-S波誘導(dǎo)轉(zhuǎn)捩發(fā)生。如果NCF先到達(dá)其對應(yīng)的閥值(NCF)tr，則認(rèn)為邊界層流動轉(zhuǎn)捩由CF波不穩(wěn)定性觸發(fā)。

1.2 RANS求解器與轉(zhuǎn)捩判斷方法耦合流程

采用耦合完全雙eN方法的三維RANS求解器PMNS3D[39]開展面向翼身組合體外形的機(jī)翼邊界層轉(zhuǎn)捩自動判斷方法研究。求解器主要由3個(gè)模塊構(gòu)成：① 三維RANS求解器模塊。本模塊采用多塊結(jié)構(gòu)化網(wǎng)格、有限體積法、LU-SGS(Lower-Upper Symmetric Gauss-Seidel)時(shí)間推進(jìn)、多重網(wǎng)格加速收斂技術(shù)。湍流流動模擬采用Spalart-Allmaras(S-A)湍流模型。此模塊對翼身組合體黏性繞流進(jìn)行數(shù)值模擬，并從計(jì)算結(jié)果中提取機(jī)翼邊界層外邊界流向速度ue、展向速度we、溫度Te等流動參數(shù)，作為三維層流邊界層方程求解的邊界條件。② 三維層流邊界層方程求解模塊。因直接從RANS方程解提取的層流邊界層信息難以滿足線性穩(wěn)定性分析的精度要求，此模塊通過非正交貼體坐標(biāo)系對機(jī)翼邊界層方程進(jìn)行求解，可在對網(wǎng)格量沒有過高要求的條件下，得到滿足三維邊界層穩(wěn)定性分析所需的高精度速度型u、溫度型T以及它們在物面法向的一階導(dǎo)數(shù)u′、T′和二階導(dǎo)數(shù)u″、T″。③ 基于線性穩(wěn)定性理論的雙eN轉(zhuǎn)捩判斷模塊。通過對邊界層內(nèi)部的速度型、溫度型進(jìn)行線性穩(wěn)定性分析，使用完全雙eN方法得到邊界層轉(zhuǎn)捩的位置xtr，并反饋給三維RANS求解模塊。重復(fù)迭代至轉(zhuǎn)捩位置收斂，依據(jù)收斂的轉(zhuǎn)捩信息繼續(xù)RANS方程的迭代求解，直至最終達(dá)到流場收斂標(biāo)準(zhǔn)。求解器中不同模塊耦合流程如圖2所示。

1.3 翼身組合體結(jié)構(gòu)化網(wǎng)格拓?fù)?/h3>

本文使用多塊結(jié)構(gòu)化網(wǎng)格，開展面向翼身組合體的雙eN轉(zhuǎn)捩自動判斷方法研究。對于三維復(fù)雜外形，如果使用非結(jié)構(gòu)網(wǎng)格，雖沒有嚴(yán)格的網(wǎng)格拓?fù)浣Y(jié)構(gòu)要求，但在計(jì)算擾動累積放大因子時(shí)，因其相鄰網(wǎng)格單元的無序性，需要采取專門的積分路徑確定方法[36]。本文針對復(fù)雜外形流動，建立了一種便于實(shí)現(xiàn)轉(zhuǎn)捩自動判斷的，且通用的翼身組合體結(jié)構(gòu)化網(wǎng)格拓?fù)浣Y(jié)構(gòu)處理方法，對所有翼面類部件均適用，在計(jì)算擾動累積放大因子時(shí)，積分路徑可以直接使用機(jī)翼表面網(wǎng)格線，相對非結(jié)構(gòu)網(wǎng)格較為簡單便利。

具體方法如下：包裹整個(gè)機(jī)翼(不包含翼梢端面和機(jī)翼后緣厚度端面)邊界層生成一塊完整的C型單塊網(wǎng)格，網(wǎng)格塊示意圖如圖3中線框所示。

對于這塊網(wǎng)格，執(zhí)行圖2中右側(cè)虛線框內(nèi)的轉(zhuǎn)捩自動判斷流程。因?yàn)闄C(jī)翼表面流向網(wǎng)格線完整連續(xù)，這樣在計(jì)算擾動累積放大因子時(shí)，積分路徑可以直接使用機(jī)翼表面網(wǎng)格線。如果需要對平尾、垂尾等翼面的邊界層進(jìn)行轉(zhuǎn)捩自動判斷，上述方法同樣適用。而機(jī)身部分則不存在網(wǎng)格拓?fù)浣Y(jié)構(gòu)的特殊要求，可按照常規(guī)方法自由生成網(wǎng)格。

圖2 RANS求解器耦合轉(zhuǎn)捩判斷模型的流程Fig.2 Coupling structure of RANS solver and transition prediction model

圖3 DLR-F4翼身組合體外形機(jī)翼邊界層轉(zhuǎn)捩 判斷結(jié)構(gòu)網(wǎng)格拓?fù)銯ig.3 Topology of structured grid for boundary layer transition prediction of wing of DLR-F4 wing-body configuration

2 方法驗(yàn)證

采用本文方法，對DLR-F4翼身組合體外形機(jī)翼邊界層進(jìn)行自由轉(zhuǎn)捩數(shù)值模擬，以驗(yàn)證本文方法的正確性。計(jì)算狀態(tài)參考DLR在歐洲跨聲速風(fēng)洞(European Transonic Wind tunnel，ETW)使用溫敏漆(Temperature Sensitive Paint，TSP)技術(shù)進(jìn)行的翼身組合體機(jī)翼邊界層轉(zhuǎn)捩測量試驗(yàn)[40]：馬赫數(shù)Ma= 0.785，雷諾數(shù)Re= 6.0×106，迎角α=-0.87°。計(jì)算網(wǎng)格為多塊結(jié)構(gòu)化網(wǎng)格，y+約為0.8，整個(gè)翼身組合體網(wǎng)格量為420萬，其中包裹機(jī)翼的網(wǎng)格塊機(jī)翼流向分布185個(gè)網(wǎng)格點(diǎn)，展向分布49個(gè)網(wǎng)格點(diǎn)，法向分布105個(gè)網(wǎng)格點(diǎn)。采用S-A湍流模型模擬湍流流動。因風(fēng)洞試驗(yàn)湍流度高于真實(shí)飛行時(shí)高空安靜大氣湍流度，放大因子轉(zhuǎn)捩閥值[(NTS)tr, (NCF)tr]根據(jù)經(jīng)驗(yàn)取[10.5, 7.5]。

使用本文方法預(yù)測的機(jī)翼邊界層轉(zhuǎn)捩位置與試驗(yàn)對比，如圖4所示，其中云圖為TSP試驗(yàn)結(jié)果[40]，淺色部分為層流區(qū)域，深色部分為湍流區(qū)域；黑色實(shí)線為本文方法計(jì)算得到的機(jī)翼上表面轉(zhuǎn)捩線，和試驗(yàn)轉(zhuǎn)捩位置基本一致。試驗(yàn)結(jié)果TSP照片顯示機(jī)翼中外段轉(zhuǎn)捩線較為光滑，而機(jī)翼內(nèi)段的轉(zhuǎn)捩線呈鋸齒狀。Fey等[40]分析認(rèn)為這可能是因?yàn)闄C(jī)翼內(nèi)段為橫流駐波不穩(wěn)定性誘發(fā)的轉(zhuǎn)捩。此外，Grabe等[24]使用引入C1橫流判據(jù)和He橫流判據(jù)的γ-Reθ t轉(zhuǎn)捩模型對此狀態(tài)的機(jī)翼邊界層轉(zhuǎn)捩判斷研究也表明，機(jī)翼內(nèi)段為橫流不穩(wěn)定性主導(dǎo)的轉(zhuǎn)捩。而本文計(jì)算結(jié)果同樣表明，翼根區(qū)域?yàn)闄M流駐波不穩(wěn)定性擾動所導(dǎo)致的轉(zhuǎn)捩：圖 5為展向25%處站位A和47%處站位B的機(jī)翼上表面T-S波(虛線)和CF波(實(shí)線)不穩(wěn)定性擾動放大積分因子N發(fā)展規(guī)律。由圖 5(a)可見，展向站位A轉(zhuǎn)捩位置(x/c)tr在流向x/c=0.34位置處，c為弦長，CF波擾動累積放大因子先于T-S波達(dá)到其對應(yīng)閥值，即CF波不穩(wěn)定性在此站位誘發(fā)轉(zhuǎn)捩。而圖 5(b)的展向站位B處為T-S波不穩(wěn)定性誘發(fā)轉(zhuǎn)捩，轉(zhuǎn)捩位置(x/c)tr在流向x/c=0.45位置處。

圖4 本文方法對DLR-F4機(jī)翼邊界層轉(zhuǎn)捩判斷結(jié)果 與TSP試驗(yàn)照片[40]對比 Fig.4 Comparison of transition prediction result of DLR-F4 wing boundary layer by proposed method and TSP test image[40]

圖5 機(jī)翼上表面T-S波和CF波不穩(wěn)定性擾動 放大積分因子Fig.5 Disturbance amplification integral factors of TSI and CFI on upper surface of wing

此算例表明，本文發(fā)展的完全雙eN方法不僅能夠?qū)C(jī)翼邊界層轉(zhuǎn)捩位置進(jìn)行較為準(zhǔn)確的預(yù)測，同時(shí)能夠判斷出機(jī)翼內(nèi)段轉(zhuǎn)捩是由橫流不穩(wěn)定性主導(dǎo)的機(jī)理，與TSP試驗(yàn)結(jié)果分析[40]和使用引入橫流判據(jù)的γ-Reθ t轉(zhuǎn)捩模型[24]的研究結(jié)論一致。驗(yàn)證了本文方法的正確性，說明本文發(fā)展的方法能夠考慮機(jī)身的影響，實(shí)現(xiàn)對機(jī)翼邊界層的轉(zhuǎn)捩自動判斷，能夠支持下文開展的中短程民機(jī)翼身組合體流動轉(zhuǎn)捩判斷研究。

3 中短程民機(jī)翼身組合體流動轉(zhuǎn)捩判斷

3.1 翼身組合體的幾何外形

采用本文發(fā)展的基于雙eN的翼身組合體流動轉(zhuǎn)捩自動判斷方法，進(jìn)行中短程民機(jī)翼身組合體流動轉(zhuǎn)捩自動判斷研究，并與單獨(dú)自然層流(NLF)機(jī)翼對比，進(jìn)一步研究了機(jī)身對機(jī)翼邊界層轉(zhuǎn)捩的影響規(guī)律，可為中短程民機(jī)自然層流后掠機(jī)翼設(shè)計(jì)提供指導(dǎo)。本文研究的翼身組合體由自然層流后掠機(jī)翼和典型民機(jī)機(jī)身外形組成，翼身組合體外形見圖6。自然層流后掠機(jī)翼展弦比為10.5，梢根比為0.3，前緣掠角為19°。在機(jī)翼翼根翼梢均配置自然層流超臨界翼型NPU-LSC-72613[41]，翼梢無扭轉(zhuǎn)。機(jī)翼上反角為0°，安裝角為2°。

圖6 自然層流翼身組合體幾何外形Fig.6 Geometry shape of NLF wing-body configuration

自然層流超臨界翼型NPU-LSC-72613的設(shè)計(jì)思想[42]是在超臨界條件下通過弱激波抬高翼型后部的負(fù)壓，使得前緣到激波位置維持順壓梯度，抑制T-S波不穩(wěn)定性，以保證在設(shè)計(jì)條件下獲得較大范圍的層流。圖7為NPU-LSC-72613翼型在Ma=0.72、Re=2.0×107、升力系數(shù)CL=0.6狀態(tài)下的表面壓力系數(shù)Cp分布及上下表面轉(zhuǎn)捩位置，從圖中可以看到翼型上下表面都存在有約50%～60%范圍的層流區(qū)域。

圖7 NPU-LSC-72613自由轉(zhuǎn)捩壓力系數(shù)分布(上表面轉(zhuǎn)捩位置(x/c)tr=0.577, 下表面轉(zhuǎn)捩位置(x/c)tr=0.507)Fig.7 Pressure coefficient distribution of NPU-LSC-72613 with free transition (transition location on upper surface: (x/c)tr=0.577; transition location on lower surface: (x/c)tr=0.507)

3.2 計(jì)算結(jié)果與分析

對翼身組合體繞流進(jìn)行了轉(zhuǎn)捩自動判斷研究，并通過與單獨(dú)機(jī)翼進(jìn)行對比，研究分析機(jī)身的存在對機(jī)翼邊界層橫流強(qiáng)度、CF波擾動累積放大因子和轉(zhuǎn)捩位置的影響規(guī)律。計(jì)算狀態(tài)為：Ma=0.78、Re=2.5×107，迎角為0°。計(jì)算網(wǎng)格為多塊結(jié)構(gòu)化網(wǎng)格，y+約為0.8，整個(gè)翼身組合體網(wǎng)格量為420萬，在包裹機(jī)翼的單塊網(wǎng)格上，機(jī)翼流向分布185個(gè)網(wǎng)格點(diǎn)，展向分布49個(gè)網(wǎng)格點(diǎn)，法向分布105個(gè)網(wǎng)格點(diǎn)。單獨(dú)機(jī)翼網(wǎng)格采取和翼身組合體網(wǎng)格中的機(jī)翼表面完全相同的網(wǎng)格分布。采用S-A湍流模型模擬湍流流動。放大因子轉(zhuǎn)捩閥值[(NTS)tr, (NCF)tr]取為[10.5, 7.5]。

圖8為有/無機(jī)身時(shí)，機(jī)翼上下表面壓力分布云圖和轉(zhuǎn)捩線。可見，單獨(dú)機(jī)翼采取自然層流翼型有效抑制T-S波，較小的前緣后掠角抑制CF波，上表面維持了較大范圍(60%)的層流，下表面維持了一定范圍(30%)的層流。而對于翼身組合體外形，由于機(jī)身的存在，機(jī)翼翼根下表面轉(zhuǎn)捩位置大幅提前，雙eN方法結(jié)果表明此處為CF波不穩(wěn)定誘導(dǎo)轉(zhuǎn)捩。

圖8 有/無機(jī)身時(shí)機(jī)翼表面壓力分布及轉(zhuǎn)捩線Fig.8 Pressure coefficient contours and transition lines of wing surface with/without fuselage

為了詳細(xì)分析機(jī)身對機(jī)翼邊界層轉(zhuǎn)捩和CF波不穩(wěn)定性擾動的影響，圖9對翼身組合體和單獨(dú)機(jī)翼展向14%處站位A和36%處站位B的CF波擾動累積放大因子進(jìn)行對比。從圖中可看出，由于機(jī)身的存在，翼身組合體相對于單獨(dú)機(jī)翼，橫流不穩(wěn)定性有所增強(qiáng)。對比圖9(a)和圖9(b)可以看出上表面兩個(gè)站位處CF波擾動累積放大因子增強(qiáng)的幅度沿展向減小。對比圖9(c)和圖9(d)可見下表面規(guī)律同樣如此。

值得注意的是，圖9(c)中展向站位A下表面，因有機(jī)身導(dǎo)致CF波擾動累積放大因子曲線在前緣處快速增長，包絡(luò)線形成的峰值從6增長到8.4，超過閥值(NCF)tr=7.5誘發(fā)轉(zhuǎn)捩， 轉(zhuǎn)捩位置從(x/c)tr=0.55大幅提前到(x/c)tr=0.03的前緣處。而隨著與機(jī)身距離的增大，圖9(d)中展向站位36%處CF波不穩(wěn)定擾動放大率的增強(qiáng)幅度變小，CF波擾動累積放大因子曲線未在前緣處到達(dá)閥值，轉(zhuǎn)捩位置僅從(x/c)tr=0.34提前到(x/c)tr=0.31。

圖9 有/無機(jī)身時(shí)機(jī)翼上表面和下表面CF波擾動累積放大因子及(NCF)peak定義Fig.9 N factors of CFI on upper surface and lower surface of wing with/without fuselage and definition of (NCF)peak

通過對各站位尤其是靠近翼根的展向站位A處下表面的對比可以看出，對于自然層流后掠機(jī)翼的設(shè)計(jì)，CF波擾動累積放大因子曲線在前緣處的峰值尤為關(guān)鍵，因?yàn)榇酥狄坏┏^(NCF)tr，橫流不穩(wěn)定性誘發(fā)的轉(zhuǎn)捩就會在前緣區(qū)域發(fā)生，必然無法實(shí)現(xiàn)較大的層流范圍。為了達(dá)到自然層流減阻的目的，這種情況應(yīng)當(dāng)盡量避免。

圖10 有/無機(jī)身時(shí)機(jī)翼下表面展向站位14%處 邊界層內(nèi)橫流速度型Fig.10 CF velocity profiles inside boundary layers at spanwise 14% station on lower surface of wing with/without fuselage

為了更好地理解機(jī)身對機(jī)翼橫流不穩(wěn)定擾動的影響，針對機(jī)身影響較大的展向站位14%，對其橫流擾動快速增長的前緣區(qū)域，即圖9(c)中a、b、c、d (x/c=0.019 0, 0.048 0, 0.064 4, 0.093 7) 這4個(gè)流向站位，進(jìn)行了有/無機(jī)身時(shí)的機(jī)翼下表面邊界層內(nèi)橫流速度型對比，見圖10，橫坐標(biāo)W/Ue為無量綱展向速度，縱坐標(biāo)Y/Ye為無量綱法向位置。橫流速度型拐點(diǎn)速度越大，拐點(diǎn)與物面距離越小，則橫流強(qiáng)度越大，對應(yīng)的CF波不穩(wěn)定擾動放大率越大。從圖中可見，在前緣加速區(qū)，隨流向速度的增大，4個(gè)流向站位處的橫流強(qiáng)度沿流向遞減。同時(shí)，每個(gè)流向站位處的橫流強(qiáng)度都因?yàn)闄C(jī)身三維位移效應(yīng)，相對單獨(dú)機(jī)翼增強(qiáng)。這就是導(dǎo)致圖9(c)中展向站位14%處CF波擾動累積放大因子在前緣區(qū)域快速增長，且因有機(jī)身進(jìn)一步增大的原因。

為了進(jìn)一步研究機(jī)身對機(jī)翼的橫流不穩(wěn)定擾動的影響沿著機(jī)翼展向的變化規(guī)律，對有/無機(jī)身時(shí)機(jī)翼下表面CF波擾動累積放大因子的包絡(luò)線進(jìn)行對比，如圖11所示，其中黑色曲線為單獨(dú)機(jī)翼，紅色曲線為翼身組合體。從圖中可以看出，展向所有CF波擾動累積放大因子的包絡(luò)線形成的包絡(luò)面，因機(jī)身三維位移效應(yīng)對機(jī)翼的影響整個(gè)鼓起，在機(jī)翼前緣區(qū)域形成一道沿機(jī)翼展向走向的山峰狀包絡(luò)面，距機(jī)身越近，其峰值越高，在翼根區(qū)域超過(NCF)tr，導(dǎo)致了翼根前緣區(qū)域CF波不穩(wěn)定誘導(dǎo)轉(zhuǎn)捩的發(fā)生。

圖11 有/無機(jī)身時(shí)機(jī)翼下表面橫流擾動累積 放大因子包絡(luò)線Fig.11 NCF envelope curves on lower surface of wing with/without fuselage

圖12 有/無機(jī)身時(shí)機(jī)翼下表面(NCF)peak沿展向 變化曲線及差值Fig.12 (NCF)peak on lower surface of wing with/without fuselage and delta value at spanwise station location

為更清晰觀察機(jī)身對機(jī)翼橫流轉(zhuǎn)捩位置的影響，圖12給出了NCF前緣區(qū)域峰值(NCF)peak(如圖9(d)中所示)沿展向的變化，同時(shí)給出了翼身組合體與單獨(dú)機(jī)翼的(NCF)peak的差值Δ(NCF)peak沿展向的變化。圖12中(NCF)peak,wing為單獨(dú)機(jī)翼下表面CF波擾動累積放大因子峰值沿展向變化曲線，因?yàn)閱为?dú)機(jī)翼翼根連接對稱面邊界條件，所以沿展向的橫向流動有一個(gè)發(fā)展過程：起初會呈現(xiàn)增強(qiáng)趨勢，很快達(dá)到最大值，然后隨著和對稱面的距離增大，機(jī)翼自身三維位移效應(yīng)減弱，橫流不穩(wěn)定擾動強(qiáng)度減小，(NCF)peak逐漸減小。(NCF)peak,wingbody為翼身組合體的機(jī)翼下表面CF波擾動累積放大因子峰值沿展向變化曲線，它處處高于(NCF)peak,wing曲線，即機(jī)身的三維位移效應(yīng)增強(qiáng)了整個(gè)機(jī)翼的橫流強(qiáng)度，翼根處橫流不穩(wěn)定擾動最強(qiáng)，隨著與機(jī)身距離的增大，強(qiáng)度沿展向逐漸減弱。Δ(NCF)peak為前兩者的差值，代表機(jī)身三維位移效應(yīng)對機(jī)翼的干擾導(dǎo)致的(NCF)peak的增長量，因此沿展向單調(diào)遞減，此遞減趨勢呈現(xiàn)出兩個(gè)階段：機(jī)身影響在翼根區(qū)域較為明顯，展向站位14%的翼根處(NCF)peak相對單獨(dú)機(jī)翼增強(qiáng)40%，在展向站位30%處迅速減小到8%；機(jī)身影響在機(jī)翼中外段相對小很多，呈現(xiàn)出基本不變的趨勢。

圖12數(shù)據(jù)表明，單獨(dú)機(jī)翼展向所有(NCF)peak都未超過(NCF)tr，即7.5，故不會發(fā)生橫流不穩(wěn)定性誘導(dǎo)的前緣轉(zhuǎn)捩，但由于機(jī)身三維位移效應(yīng)的干擾，導(dǎo)致機(jī)翼下表面在展向站位35%內(nèi)超過7.5，導(dǎo)致原本設(shè)計(jì)的自然層流后掠機(jī)翼，其內(nèi)段部分轉(zhuǎn)捩位置提前到前緣區(qū)域，如展向站位14%處轉(zhuǎn)捩位置從(x/c)tr=0.55大幅提前到(x/c)tr=0.03的前緣處。

4 結(jié) 論

本文將求解線性穩(wěn)定性方程的完全雙eN轉(zhuǎn)捩自動判斷方法推廣到了三維復(fù)雜外形。將針對實(shí)際應(yīng)用設(shè)計(jì)的自然層流機(jī)翼配置到典型民機(jī)機(jī)身上，對此翼身組合體流動轉(zhuǎn)捩特征進(jìn)行了研究，分析了機(jī)身的存在對機(jī)翼的影響。

1) 利用DLR-F4翼身組合體模型對所發(fā)展的方法進(jìn)行了檢驗(yàn)。結(jié)果表明該方法判斷出的轉(zhuǎn)捩位置與試驗(yàn)結(jié)果吻合較好，驗(yàn)證了方法的正確性。

2) 本文方法捕捉到了DLR-F4翼根處由橫流駐波導(dǎo)致的轉(zhuǎn)捩現(xiàn)象。表明本文發(fā)展的針對三維復(fù)雜外形的轉(zhuǎn)捩自動判斷方法不僅能夠較為準(zhǔn)確地判斷轉(zhuǎn)捩位置，還能夠正確判斷轉(zhuǎn)捩誘發(fā)機(jī)制。

3) 對自然層流中短程民機(jī)外形的轉(zhuǎn)捩特征與橫流不穩(wěn)定性的定量分析表明，機(jī)身的三維位移效應(yīng)對流動產(chǎn)生了橫向的排擠，這增強(qiáng)了機(jī)翼邊界層的橫流強(qiáng)度，使得前緣區(qū)域的橫流不穩(wěn)定性變強(qiáng)，導(dǎo)致翼根提前發(fā)生轉(zhuǎn)捩。

[1] COLLIER F, THOMAS R, BURLEY C, et al. Environmentally responsible aviation—Real solutions for environmental challenges facing aviation[C]∥27th Congress of the International Council of the Aeronautical Sciences 2010. Stockholm: ICAS Secretariat, 2010: 300-315.

[2] DARECK M, EDELSTENN C, ENDER T, et al. Flightpath 2050 Europe’s vision for aviation[R]. Belgium: Advisory Council for Aeronautics Research in Europe, 2011.

[3] HEPPERLE M. MDO of forward swept wings[C]∥KATnet II Multi-Disciplinary Design and Configuration Optimization Workshop, 2008: 28-30.

[4] ARNAL D, COUSTOLS E, JUILLEN J C. Experimental and theoretical study of transition phenomena on an infinite swept wing[J]. La Recherche Aerospatiale (English Edition), 1984, 1984(4): 39-54.

[5] DRELA M. Newton solution of coupled viscous/inviscid multielement airfoil flows[C]∥21st Fluid Dynamics, Plasma Dynamics and Lasers Conference, 1990: 1470.

[6] CEBECI T, STEWARTSON K. On stability and transition in three-dimensional flows[J]. AIAA Journal, 1980, 18(4): 398-405.

[7] MALIK M R. COSAL: A black-box compressible stability analysis code for transition prediction in three-dimensional boundary layers: NASA-CR-165925[R]. Washington，D.C.: NASA, 1982.

[8] MACK L M. Stability of three-dimensional boundary layers on swept wings at transonic speeds[C]∥IUTAM. Symposium Transsonicum III. Berlin Heidelberg: Springer, 1989: 209-223.

[9] ARNAL D. Transition prediction in transonic flow[C]∥IUTAM. Symposium Transsonicum III. Berlin Heidelberg: Springer, 1989: 253-262.

[10] ARNAL D, CASALIS G, JUILLEN J C. Experimental and theoretical analysis of natural transition on “infinite” swept wing[C]∥IUTAM. Laminar-Turbulent Transition. Berlin Heidelberg: Springer, 1990: 311-325.

[11] RADESPIEL R, GRAAGE K, BRODERSEN O. Transition predictions using Reynolds-averaged Navier-Stokes and linear stability analysis methods: AIAA-1991-1641[R]. Reston, VA: AIAA, 1991.

[12] CASALIS G, ARNAL D. ELFIN II subtask 2.3: Database method—Development and validation of the simplified method for pure cross-flow instability at low speed: ELFIN II-145[R]. Toulouse: ONERA-CERT, 1996.

[13] STOCK H W, HAASE W. Feasibility study of eNtransition prediction in Navier-Stokes methods for airfoils[J]. AIAA Journal, 1999, 37(10): 1187-1196.

[14] NEBEL C, RADESPIEL R, WOLF T. Transition prediction for 3D flows using a Reynolds-averaged Navier-Stokes code andN-factor methods: AIAA-2003-3593[R]. Reston, VA: AIAA, 2003.

[15] KRUMBEIN A. Automatic transition prediction and application to high-lift multi-element configurations: AIAA-2004-2543[R]. Reston, VA: AIAA, 2004.

[16] STOCK H W. Infinite swept wing Navier-Stokes computations with eNtransition prediction[J]. AIAA Journal, 2005, 43(6): 1221-1229.

[17] KRUMBEIN A. Automatic transition prediction and application to 3D wing configurations: AIAA-2006-0914[R]. Reston, VA: AIAA, 2006.

[18] CLIQUET J, HOUDEVILLE R, ARNAL D. Application of laminar-turbulent transition criteria in Navier-Stokes computations: AIAA-2007-0515[R]. Reston, VA: AIAA, 2007.

[19] KRUMBEIN A. eNtransition prediction for 3D wing configurations using database methods and a local, linear stability code[J]. Aerospace Science and Technology, 2008, 12(8): 592-598.

[20] PERRAUD J, ARNAL D, CASALIS G, et al. Automatic transition predictions using simplified methods[J]. AIAA Journal, 2009, 47(11): 2676-2684.

[21] PERRAUD J, CLIQUET J, HOUDEVILLE R, et al. Transport aircraft three-dimensional high-lift wing numerical transition prediction[J]. Journal of Aircraft, 2008, 45(5): 1554-1563.

[22] KRIMMELBEIN N, RADESPIEL R, NEBEL C. Numerical aspects of transition prediction for three-dimensional configurations: AIAA-2005-4764[R]. Reston, VA: AIAA, 2005.

[23] LEE J D, JAMESON A. Natural-laminar-flow airfoil and wing design by adjoint method and automatic transition prediction: AIAA-2009-0897[R]. Reston, VA: AIAA, 2009.

[24] GRABE C, SHENGYANG N, KRUMBEIN A. Transition transport modeling for the prediction of crossflow transition: AIAA-2016-3572[R]. Reston, VA: AIAA, 2016.

[25] 張坤, 宋文萍. eN方法在無限展長后掠翼邊界層轉(zhuǎn)捩判斷中的初步應(yīng)用[J]. 西北工業(yè)大學(xué)學(xué)報(bào), 2011, 29(1): 142-147.

ZHANG K, SONG W P. Application of the full eNtransition method to the infinite swept-wing’s transition prediction[J]. Journal of Northwestern Polytechnical University, 2011, 29(1): 142-147 (in Chinese).

[26] 左歲寒, 楊永, 李棟. 基于線性拋物化穩(wěn)定性方程的后掠翼邊界層內(nèi)橫流穩(wěn)定性研究[J]. 計(jì)算物理, 2010, 27(5): 665-670.

ZUO S H, YANG Y, LI D. Investigation on cross-flow instabilities in swept-wing boundary layers with linear parabolized stability equations[J]. Chinese Journal of Computational Physics, 2010, 27(5): 665-670 (in Chinese).

[27] 孫朋朋, 黃章峰. 后掠角對后掠機(jī)翼邊界層穩(wěn)定性及轉(zhuǎn)捩的影響[J]. 北京航空航天大學(xué)學(xué)報(bào), 2015, 41(7): 1313-1321.

SUN P P, HUANG Z F. Effect of sweep angle on stability and transition in a swept-wing boundary layer[J]. Journal of Beijing University of Aeronautics and Astronautics, 2015, 41(7): 1313-1321 (in Chinese).

[28] 靖振榮, 孫朋朋, 黃章峰. 小攻角對后掠機(jī)翼邊界層穩(wěn)定性及轉(zhuǎn)捩的影響[J]. 北京航空航天大學(xué)學(xué)報(bào), 2015, 41(11): 2177-2183.

JING Z R, SUN P P, HUANG Z F. Effect of attack angle on stability and transition in a swept-wing boundary layer[J]. Journal of Beijing University of Aeronautics and Astronautics, 2015, 41(11): 2177-2183 (in Chinese).

[29] 黃章峰, 逯學(xué)志, 于高通. 機(jī)翼邊界層的橫流穩(wěn)定性分析和轉(zhuǎn)捩預(yù)測[J]. 空氣動力學(xué)學(xué)報(bào), 2014, 32(1): 14-20.

HUANG Z F, LU X Z, YU G T. Cross-flow instability analysis and transition prediction of airfoil boundary layer[J]. Acta Aerodynamica Sinica, 2014, 32(1): 14-20 (in Chinese).

[30] HAN Z H, CHEN J, ZHU Z, et al. Aerodynamic design of transonic natural-laminar-flow (NLF) wing via surrogate-based global optimization: AIAA-2016-2041[R]. Reston, VA: AIAA, 2016.

[31] 徐家寬, 白俊強(qiáng), 喬磊, 等. 后掠翼邊界層橫流不穩(wěn)定轉(zhuǎn)捩預(yù)測模型[J]. 航空動力學(xué)報(bào), 2015, 30(4): 927-935.

XU J K, BAI J Q, QIAO L, et al. Prediction model of cross-flow instability transition in swept wing boundary layers[J]. Journal of Aerospace Power, 2015, 30(4): 927-935 (in Chinese).

[32] 徐家寬, 白俊強(qiáng), 喬磊, 等. 橫流不穩(wěn)定性轉(zhuǎn)捩預(yù)測模型[J]. 航空學(xué)報(bào), 2015, 36(6): 1814-1822.

XU J K, BAI J Q, QIAO L, et al. Transition model for predicting crossflow instabilities[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(6): 1814-1822 (in Chinese).

[33] 史亞云, 白俊強(qiáng), 華俊, 等. 基于當(dāng)?shù)刈兞康臋M流轉(zhuǎn)捩預(yù)測模型的研究與改進(jìn)[J]. 航空學(xué)報(bào), 2016, 37(3): 780-789.

SHI Y Y, BAI J Q, HUA J, et al. Study and modification of cross-flow induced transition model based on local variables[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(3): 780-789 (in Chinese).

[34] 戚瓊, 韓慶. 基于Spalart-Allmaras-γ-Reθ t轉(zhuǎn)捩模型的橫流不穩(wěn)定性轉(zhuǎn)捩預(yù)測方法[J]. 氣體物理, 2016, 1(3): 19-24.

QI Q, HAN Q. Prediction method of cross-flow instabilities-induced transition based on Spalart-Allmaras-γ-Reθ ttransition model[J]. Physics of Gases, 2016, 1(3): 19-24 (in Chinese).

[35] 鞠勝軍, 閻超, 葉志飛.γ-Reθ t-CF轉(zhuǎn)捩模型在Spalart-Allmaras湍流模型中的推廣及驗(yàn)證[J]. 航空學(xué)報(bào), 2017, 38(4): 120383.

JU S J, YAN C, YE Z F. Genevalization and validation ofγ-Reθ t-CF transition modeling in combination with Spalart-Allmaras turbulence model[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(4): 120383 (in Chinese).

[36] KRUMBEIN A, KRIMMELBEIN N, GRABE C. Streamline-based transition prediction techniques in an unstructured computational fluid dynamics code[J]. AIAA Journal, 2017, 55(5): 1548-1564.

[37] SMITH A M O, GAMBERONI N. Transition, pressure gradient and stability theory[M]. Long Beach: Douglas Aircraft Company, 1956.

[38] INGEN J L V. A suggested semi-empirical method for the calculation of the boundary layer transition region: VTH-74[R]. Delft: Delft University of Technology, 1956.

[39] SONG W P, ZHU Z, YANG H, et al. Laminar flow wing’s optimization design by RANS solver with automatic transition prediction[C]∥28th Congress of the International Council of the Aeronautical Sciences 2012. Stockholm: ICAS Secretariat, 2012: 716-725.

[40] FEY U, EGAMI Y, ENGLER R. High Reynolds number transition detection by means of temperature sensitive paint: AIAA-2006-0514[R]. Reston, VA: AIAA, 2006.

[41] HAN Z H, DENG J, LIU J, et al. Design of laminar supercritical airfoils based on Navier-Stokes equations[C]∥28th Congress of the International Council of the Aeronautical Sciences 2012. Stockholm: ICAS Secretariat, 2012: 706-715.

[42] 喬志德. 自然層流超臨界翼型的設(shè)計(jì)研究[J]. 流體力學(xué)實(shí)驗(yàn)與測量, 1998, 12(4): 23-30.

QIAO Z D. Design of supercritical airfoils with natural laminar flow[J]. Experiments and Measurements in Fluid Mechanics, 1998, 12(4): 23-30 (in Chinese).