徐佳敏，宋文濱

（上海交通大學 航空航天學院，上海 200240）

0 引 言

飛機增升裝置的設(shè)計直接關(guān)系到飛機的起飛和著陸場長要求，以及第二階段爬升率等性能參數(shù)，同時越來越高的適航噪聲標準逐漸成為飛機氣動設(shè)計中的重要約束因素。這使得在增升裝置設(shè)計中越來越需要考慮氣動效率指標和噪聲特性的多目標優(yōu)化，縫翼空腔內(nèi)的復(fù)雜流動對于以上兩點的影響都至關(guān)重要，是主要的研究問題之一［1－3］。深入理解縫翼流場的流動機理對實現(xiàn)高氣動效率和低噪聲水平的增升裝置設(shè)計尤為重要，通常需要將數(shù)值模擬和實驗研究緊密結(jié)合［4］。由于前緣縫翼幾何外形導(dǎo)致其流場結(jié)構(gòu)非常復(fù)雜，涉及很多復(fù)雜的物理流動現(xiàn)象，如大范圍的低速區(qū)域以及主翼前緣可能的高速區(qū)域，大壓力梯度，強三維展向效應(yīng)，邊界層和尾跡的混合等。這些流動對于準確的數(shù)值模擬而言都存在挑戰(zhàn)，本文的焦點在于前緣縫翼的非定常流動問題，為了避免襟翼流動對縫翼流場的影響，本文將研究對象限于有限展長的兩段翼型。

對機體氣動噪聲問題中聲源信息的數(shù)值分析依賴于對壓力脈動的準確模擬，可用的數(shù)值方法主要包括在定常CFD（Computational Fluid Dynamics）結(jié)果上增加聲源信息，例如隨機噪聲產(chǎn)生輻射（Stochastic Noise Generation and Radiation，SNGR）方法［5－6］，或者使用非定常方法。常用的非定常方法包括非定常雷諾平 均 （Unsteady Reynolds－Average Navier－Stokes，URANS），大渦模擬（Large Eddy Simulation，LES）方法和脫體渦模擬（Detached Eddy Simulation，DES）方法。非定常雷諾平均仍然是當前工業(yè)應(yīng)用較多的非定常數(shù)值模擬方法。基于非定常雷諾平均模型的多段翼型數(shù)值模擬已經(jīng)有了很多相關(guān)研究，Jenkins等［7］和 Khorrami等［8］對30P30N 多段翼構(gòu)型基于非定常雷諾平均方法在不同攻角下進行了計算，與PIV（Particle Image Velocimetry）實驗結(jié)果進行了對比，并且討論了縫翼噪聲特性隨攻角的變化；Khorrami等在用非定常雷諾平均方法捕捉到了有限厚度縫翼尾緣的高頻渦脫落［9］；K?nig等使用了大渦模擬方法對縫翼流場進行了非定常計算并且將聲源提取，再使用聲擾動方程（Acoustic Perturbation Equation，APE）方法計算了縫翼噪聲的傳播特性［10］；Ma和Zhang也使用了大渦模擬方法結(jié)合聲學擾動方程方法對聲襯在縫翼上的應(yīng)用進行了研究［11］。

非定常雷諾平均方法雖然能夠用較高的效率得到大致準確的流場，但是由于求解依靠的是基于經(jīng)驗推導(dǎo)的湍流模型，因此其時均特性對小尺度的脈動會產(chǎn)生抑制作用，對流場中細節(jié)的小尺度流動特征均無法很好的捕捉；大渦模擬方法雖然能夠捕捉到流場中的小尺度的細節(jié)流動，但是對于現(xiàn)在的計算機水平而言，高雷諾數(shù)下伴隨薄邊界層所需的巨大網(wǎng)格量對計算資源的要求使得這種方法仍然無法進入實際工程應(yīng)用?；旌戏椒ǖ陌l(fā)展則能夠較好的解決這一矛盾，Imamura等［12］沿弦長方向?qū)φ麄€流域進行了雷諾平均區(qū)域以及大渦模擬區(qū)域的劃分，縫翼空腔以及主翼的前緣均屬于大渦模擬區(qū)域，因此能夠較好的抓住縫翼空腔及其與主翼間隙內(nèi)細節(jié)流動特征，通過頻譜分析得到的結(jié)果噪聲結(jié)果也與實驗結(jié)果比較吻合；法國宇航局（ONERA）的Deck發(fā)展了分區(qū)脫體渦模擬（Zonal Detached Eddy Simulation，Zonal－DES）方法用于復(fù)雜外形的非定常計算［3，13］，即自定義脫體渦／雷諾平均計算域的方法來計算三段翼型構(gòu)型的流場，對縫翼空腔氣動噪聲預(yù)測方面也有很好的表現(xiàn)。對于多段翼型的數(shù)值模擬國內(nèi)也開展了較多的工作，但較少涉及面向氣動噪聲問題的研究。宋科使用脫體渦模擬方法對帶襟翼的多段翼型進行了計算研究工作，但沒有涉及縫翼空腔內(nèi)的復(fù)雜流動的細致分析［14］；焦予秦同樣使用脫體渦模擬方法進行了30P30N多段翼型的非定常計算，但計算局限在二維構(gòu)型，無法得到三維效應(yīng)及影響［15］；王運濤利用基于雷諾平均方法的高精度數(shù)值格式對30P30N多段翼型進行了計算，得到結(jié)果與實驗符合良好，并且對轉(zhuǎn)捩位置進行了研究［16－17］；黃華使用Fluent中脫體渦模擬方法進行了基于前緣射流方法的縫翼噪聲控制技術(shù)研究［18］。

脫體渦模擬方法在是一種“弱耦合”的混合方法，因為模化的湍流能量與解析的湍流能量之間沒有轉(zhuǎn)換機制［19］。該方法主要存在的問題是由于?；瘧?yīng)力缺失（Modeled Stress Depletion，MSD）導(dǎo)致的網(wǎng)格誘導(dǎo)分離（Grid Induced Separation，GIS）。Spalart在原脫體渦模擬方法中加入了一個阻尼函數(shù)，形成了延遲脫體渦模擬方法（Delayed Detached Eddy Simulation，DDES）較好的解決了網(wǎng)格誘導(dǎo)分離這一非物理現(xiàn)象［20］；Deck［3］的分區(qū)脫體渦模擬方法通過自定義的湍流作用區(qū)域也能很好的克服這一問題。

本文借鑒Deck的分區(qū)脫體渦技術(shù)的思路，發(fā)展了分區(qū)LES／DES混合方法，并對典型構(gòu)型兩段翼型縫翼空腔內(nèi)的三維復(fù)雜細節(jié)流動及其噪聲特性進行了計算分析和討論。本文結(jié)構(gòu)安排為：第二部分給出了基本流動方程和主要的數(shù)值方法；第三部分介紹本文使用的混合方法，計算條件以及計算結(jié)果和討論；第四部分給出主要結(jié)論。

1 控制方程和數(shù)值方法

1.1 控制方程

在一般曲線坐標系下，無量綱化以后的三維可壓縮守恒形式Navier－Stokes方程：

式中，ξ、η、ζ為一般曲線坐標系下三個方向坐標，Q為守恒變量，E、，F(xiàn)、G 為無粘通量，Ev、Fv、Gv為粘性通量，Re為雷諾數(shù)。方程通過特征弦長L∞，自由來流密度ρ∞和速度u∞進行無量綱化處理。本文使用基于多塊結(jié)構(gòu)化網(wǎng)格的高精度有限差分方法來求解方程（1）。

1.2 空間離散

流場數(shù)值計算的空間離散方法主要包括有限差分和有限體積法，在有限差分格式中，與中心差分格式如Jameson中心差分格式［21］相比，迎風格式具有更好的數(shù)值穩(wěn)定性和穩(wěn)健性。本文使用基于通量差分分裂方法（Finite Difference Splitting，F(xiàn)DS）的 Roe格式，可表述為：

式中，F(xiàn)L和FR分別為控制體表面沿其法向左右兩側(cè)的無粘通量矢量，A為采用Roe平均變量計算得到的通量雅可比矩陣，QL和QR分別為控制體表面左右兩側(cè)的守恒變量。

空間離散的精度則通過高精度的插值及重構(gòu)來得到，本文采用的是五階WENO格式來重構(gòu)守恒變量QL和QR，其中：

式中，ε的引入是為了避免分母為0的情況發(fā)生，ε應(yīng)當取一個足夠小的數(shù)。如果ε的取值越大，則離散格式的數(shù)值耗散越小［22］，由于本文計算工況為低速無激波狀態(tài)，因此ε可以取稍大值ε＝1.3。

1.3 時間推進

本文在時間推進方法上采用基于LU－SGS隱式格式的向后的三點二階時間推進方法，公式如下。

式中，R為方程右端項，Δt為時間步長。隱式LUSGS格式具體方法如下。

1.4 湍流模擬

基于SA 模型（Spalart－Allmaras Model）的脫體渦數(shù)值模擬方法是將SA方程中任意網(wǎng)格節(jié)點到最近的物面距離d用～d來表示：

式中，Δ＝max（Δx，Δy，Δz），為網(wǎng)格的單元最大尺度；CDES為常數(shù)取0.65。

Smagorinsky亞格子應(yīng)力模型是大渦模擬方法中相對易于實現(xiàn)以及計算高效的亞格子模型，解得湍流渦粘系數(shù)：

2 基于分區(qū)混合方法的數(shù)值模擬

2.1 改進的混合方法策略

Deck提出的分區(qū)脫體渦方法采用分區(qū)指定湍流模型的處理方法，通過在縫翼空腔內(nèi)使用脫體渦方法進行模擬［3］，而在其他區(qū)域均使用雷諾平均方法進行處理（圖1a），旨在保證能夠在所重點關(guān)注的縫翼空腔區(qū)域內(nèi)盡可能保證捕捉到非定常流動的細節(jié)，同時在流場其余部分使用雷諾平均方法來最大程度上降低計算成本。從圖2中的兩段翼型熵增云圖可以看到，除了縫翼空腔內(nèi)復(fù)雜流動外，來流流過縫翼尾緣與縫翼下表面的加速氣流混合產(chǎn)生的尾跡會流過整個主翼弦長，然后隨主翼尾跡一起開始耗散，縫翼尾緣產(chǎn)生的尾跡也是縫翼噪聲的一個重要組成部分，也是襟縫翼流動相互影響的原因之一。因此在計算能力相對提高的今天，在縫翼空腔內(nèi)使用大渦模擬方法進行模擬能夠幫助更好的來捕捉縫翼空腔內(nèi)的細節(jié)流動，在其余區(qū)域均選擇脫體渦方法進行計算（圖1b）相比之下能夠提高縫翼流場一些次要流動細節(jié)的捕捉能力，從而提高噪聲預(yù)測的準確性。在計算特性方面，本文方法對計算網(wǎng)格的要求要高于分區(qū)脫體渦方法，尤其在縫翼空腔內(nèi)，大渦模擬方法對于網(wǎng)格的要求相對于脫體渦方法要更高；在非核心區(qū)域的脫體渦網(wǎng)格也需要在有明顯非定常流動的區(qū)域進行加密，因此也會導(dǎo)致在計算量上相對分區(qū)脫體渦方法要更大，但是相比傳統(tǒng)的全場大渦模擬方法計算量依舊要小很多。

圖1 區(qū)域混合方法策略Fig.1 Zonal LES／DES method strategy

圖2 兩段翼型熵增云圖Fig.2 Entropycontour of two－element airfoil

由于產(chǎn)生氣動噪聲的壓力脈動的量級遠遠小于流場平均量的量級，對數(shù)值方法計算精度提出很高的要求，因此本文選取的是基于有限差分方法的高精度格式求解，有限體積方法因為其穩(wěn)定性和魯棒性而非常適合于復(fù)雜外形的計算，但是在高精度計算方面，計算更高效是有限差分方法的優(yōu)勢，特別是在采用高密度網(wǎng)格，針對流動或噪聲機理問題開展研究的問題中。

2.2 計算條件與計算網(wǎng)格

本文的研究重點在于前緣縫翼，選擇的構(gòu)型為兩段翼型［10］，直接選取兩段翼型的好處在于可以隔離后緣縫翼對流場特性的影響，突出前緣縫翼的噪聲機理。無量綱縫翼弦長為0.15，尾緣厚度為0.0007。計算工況為：自由來流馬赫數(shù)0.16，攻角13°，基于干凈機翼氣動弦長（0.4m）以及自由來流速度的雷諾數(shù)為1.4×106。

圖3給出了計算中采用的計算域和網(wǎng)格分布，其中圖3（a）顯示整個流域的計算網(wǎng)格；圖3（b）給出了兩段翼型近場網(wǎng)格拓撲和分布；由于在縫翼空腔內(nèi)使用的是大渦模擬方法，因此需要盡可能保證此處的網(wǎng)格尺度在各個方向上的一致性（如圖3c所示）。有厚度的后緣會產(chǎn)生周期性的渦脫落，是縫翼高頻噪聲的主要來源，因此縫翼后緣網(wǎng)格需要一定程度上的加密；縫翼的前緣尖端處幾何變化很大，從流動角度考慮此處為縫翼下緣加速氣流與空腔內(nèi)低速回流形成的非定常剪切層區(qū)域，需要足夠好的幾何過渡來保證網(wǎng)格變換矩陣的連續(xù)性，以此保證此處剪切層區(qū)域計算的精確性。展向網(wǎng)格尺度的選取應(yīng)參考縫翼空腔內(nèi)的平面網(wǎng)格尺度，盡可能保證縫翼空腔內(nèi)的網(wǎng)格三個方向尺度一致，考慮到流動的展向效應(yīng)稍次要與平面內(nèi)的渦結(jié)構(gòu)捕捉，展向網(wǎng)格尺度選取，共包含21個節(jié)點，模型展向無量綱長度0.02，約為縫翼弦長的10%。流域整體網(wǎng)格遠場無量綱尺度半徑為20，截面網(wǎng)格的網(wǎng)格點為33萬，全場網(wǎng)格量約為735萬。

圖3 計算網(wǎng)格Fig.3 Computational grid

時間推進采用向后三點二階精度方法，取無量綱時間步長為0.0005s，從而保證在縫翼空腔內(nèi)除邊界層外其余CFL數(shù)均保持在1左右。翼型表面采用無滑移固壁邊界條件，展向兩端為周期性邊界條件，遠場采用基于一維黎曼不變量的無反射遠場邊界條件。

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

本節(jié)將計算結(jié)果和可用的風洞試驗數(shù)據(jù)進行了對比。圖4給出了本文計算得到的50%展長位置非定常時間平均壓力系數(shù)與實驗結(jié)果對比［10］，計算結(jié)果與實驗吻合良好。圖5為縫翼空腔處非定常時間平均馬赫數(shù)云圖，截面A至D的速度剖面（線段的法向速度）與PIV實驗的對比如圖6所示，兩者吻合良好，對本文非定常計算進行了驗證。在截面B與C速度剖面的近壁附近因為原試驗中PIV粒子數(shù)量較少，所以測量可能存在一定誤差［10］。圖7兩段翼型瞬態(tài)渦量分布云圖，可以清晰的看到計算捕捉到了縫翼后緣尾跡渦量從比較集中的狀態(tài)發(fā)展到離散狀態(tài)的過程。圖8為縫翼空腔內(nèi)瞬態(tài)馬赫數(shù)分布云圖，表明縫翼空腔是一個有很強非定常效應(yīng)的低速區(qū)域，同時可以看到縫翼內(nèi)側(cè)表面處二次分離的痕跡，二次分離的現(xiàn)象往往與計算模型的展向長度密切相關(guān)。圖9為縫翼尖端的壓力云圖，使用高精度的計算格式（圖9b）相比低階格式（圖9a）可以更好的捕獲到此處的壓力脈動脫落的過程。

圖4 非定常時均壓力系數(shù)分布Fig.4 Time averaged pressure coefficient distribution

圖5 不同速度剖面截面位置分布Fig.5 Velocity profile at different section

圖6 A～D截面位置速度剖面Fig.6 Velocity profile at section A～D

圖7 渦量分布云圖Fig.7 Contour of vorticity magnitude

圖8 縫翼空腔瞬態(tài)馬赫數(shù)云圖（t＝0.05s）Fig.8 Mach number contour in slat cove region

圖9 縫翼下緣渦脫落Fig.9 Pressure shedding on the slat cusp

圖10為縫翼空腔Q標準等值面圖，從圖中可以看到在縫翼下緣處剪切層的Kelvin－Helmholtz不穩(wěn)定性產(chǎn)生的二維展向渦（圖10b）；展向渦系隨剪切層流動發(fā)展空間尺度會變大，單位能量尺度變小，但同時也會產(chǎn)生展向效應(yīng)；渦撞擊到縫翼下表面后部分進入回流區(qū)，到達縫翼尖端區(qū)域，與新生剪切層混合形成一個反饋機制；撞擊到縫翼下表面后另一部分的渦沿著縫翼下表面形成一個加速，同時演化為流向渦（圖10c），而后與縫翼后緣上表面的加速流混合形成尾跡。

圖10 縫翼空腔Q標準等值面圖（馬赫數(shù)云圖分布）Fig.10 ISO surface of Q－criterion（contoured with Mach number）

圖11中給出縫翼空腔周圍設(shè)置的壓力監(jiān)測點分布示意圖，這些壓力監(jiān)測點分別布置在在縫翼后緣、內(nèi)壁，以及縫翼空腔內(nèi)部等感興趣的位置。圖12為縫翼后緣附近的監(jiān)測點1～3，從頻譜上可以非常清晰的看到高頻壓力脈動的存在，這是由縫翼尾緣上下表面兩股加速氣流混合帶來的高頻渦脫落引起的；圖14為縫翼空腔右下方的壓力監(jiān)測點9～11的噪聲頻譜分析圖，此處的壓力脈動表現(xiàn)出寬頻特性，由于監(jiān)測點11離縫翼尖端距離較近，因此能夠看到由于縫翼尖端脫渦引起的高頻特性；監(jiān)測點4～8放置在了整個回流區(qū)周圍（圖13），監(jiān)測點7和8分布在縫翼尖端產(chǎn)生的強剪切層上，在剪切層中各個尺度的壓力脈動能量均比較大，因為流動本身的壓力脈動與聲場的壓力脈動均混合在了一起。監(jiān)測點7和監(jiān)測點11一樣距離縫翼尖端較近，但是監(jiān)測點7處于剪切層上，各個頻率段的能量都非常大，因此高頻段在監(jiān)測點7沒有顯得很突出，并且真實情況下縫翼下緣的厚度非常小，因此這種高頻效應(yīng)在真實飛行情況下并不明顯。

圖11 縫翼空腔周圍壓力監(jiān)測點分布Fig.11 Pressure observation at points near the slat cove

圖12 縫翼尾緣壓力監(jiān)測點1～3頻譜圖Fig.12 Pressure observation points 1～3 near the slat trailing edge

圖13 縫翼空腔壓力監(jiān)測點4～8頻譜圖Fig.13 Pressure observation at points 4～8in the slat cove

圖14 縫翼空腔周圍壓力監(jiān)測點9～11頻譜圖Fig.14 Pressure observation at points 9～11near the slat cove

3 結(jié) 論

本文提出一種高精度的分區(qū)LES／DES數(shù)值模擬方法，并針對有限展長的兩段翼型進行數(shù)值模擬和分析，在縫翼空腔內(nèi)的復(fù)雜分離流動區(qū)域采用大渦模擬方法進行求解，能夠明顯解決湍流模型對小尺度流動發(fā)展的壓制作用，進而捕捉到縫翼空腔內(nèi)流動小尺度流動，以及大能量尺度的渦發(fā)展演繹為小能量尺度等細節(jié)過程；而在非核心區(qū)域內(nèi)使用脫體渦數(shù)值模擬方法相比雷諾平均方法能夠更好的捕捉到流場中的縫翼尾跡流動，解得的流場更接近真實物理背景。通過計算結(jié)果與已發(fā)表實驗數(shù)據(jù)的比較，證實了方法的有效性，最后通過縫翼附近的壓力觀測點數(shù)據(jù)對縫翼壓力脈動特點進行了分析。

［1］K?NIG D，KOH S R，SCHRODER W，et al.Slat noise source identification［C］.15thAIAA／CEAS Aeroacoustics Conference，May 2009，Miami，AIAA 2009－3100.

［2］AGARWAL A，MORRIS P J.Investigation of the physical mechanisms of tonal sound generation by slats［C］.8thAIAA／CEAS Aeroacoustics Conference，June 2002，Breckenridge，AIAA 2002－2575.

［3］DECK S.Zonal－detached－eddy simulation of the flow around a high－lift configuration［J］.AIAA Journal，2005，43（11）.

［4］張斌，詹浩，朱軍.飛機增升裝置的數(shù)值模擬研究［J］.航空工程進展，2011，2（1）.

［5］FRANCESCANTONIO P，F(xiàn)ERRANTE P，DECONINCK T，et al.Assessment of SNGR method for robust and efficient simulations of flow generated noise［C］.19thAIAA／CEAS Aeroacoustics Conference，May 2013，Berlin，AIAA 2013－2264.

［6］YAO H D，DAVIDSON L，ERIKSSON L E，GRUNDESTAM O，et al.Surface integral analogy approaches to computing noise generated by a 3dhigh－lift wing configuration［C］.50thAIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition，January 2012，Nashvill，AIAA 2012－0386.

［7］JENKINS L N.，KHORRAMI M R，CHOUDHARI M M.Characterization of unsteady flow structures near leading－edge slat：Part I.PIV measurements［C］.10thAIAA／CEAS Aeroacoustics Conference， May 2004， Manchester，AIAA 2004－2801.

［8］KHORRAMI M R，CHOUDHARI M M，JENKINS L N.Characterization of unsteady flow structures near leading－edge slat：Part II.2Dcomputations［C］.10thAIAA／CEAS Aeroacoustics Conference，May 2004，Manchester，AIAA 2004－2802.

［9］KHORRAMI M R，BERKMAN M E，CHOUDHARI M M.Unsteady flow computations of a slat with a blunt trailing edge［J］.AIAA Journal，2000，38（11）：2050－2058.

［10］K?NIG D，KOH S R，MEINKE M，et al.Two－step simulation of slat noise［J］.Computers ＆ Fluids，2010，（39）：512－524.

［11］MA Z K，ZHANG X.Numerical investigation of broadband slat noise attenuation with acoustic liner treatment［J］.AIAA Journal，2009，47（12）.

［12］IMAMURA T，ENOMOTO S，YOKOKAWA Y，et al.Threedimensional unsteady flow computations around a conventional slat of high－lift devices［J］.AIAA Journal，2008，46（5）：1045－53.

［13］DECK S.Recent improvements in the zonal detached eddy simu－lation（ZDES）formulation［J］.Theor.Comput.Fluid Dyn.，2012，26：523－550.

［14］宋科，喬志德.多段翼型大迎角分離流動的DELAYED RANS／LES混合算法［J］.航空計算技術(shù)，2009，39（3）：42－47.

［15］焦予秦，范娟莉，羅淞.多段翼型流動非定常性計算研究［J］.科學技術(shù)與工程，2011，11（13）：2994－2998.

［16］王運濤，王光學，張玉倫.30P－30N多段翼型復(fù)雜流場數(shù)值模擬技術(shù)研究［J］.空氣動力學學報，2010，28（1）：99－103.

［17］王運濤，孟德虹，鄧小剛.多段翼型高精度數(shù)值模擬技術(shù)研究［J］.空氣動力學學報，2013，31（1）：88－93.

［18］黃華.基于前緣射流的縫翼噪聲控制研究［D］.上海：上海交通大學.2013.

［19］AUPOIX A，ARNAL D，BEZARD H，et al.Transition and turbulence modeling［J］.Aerospace Lab.Journal，2011，（2）.

［20］SPALART P R，DECK S，SHUR，et al.A new version of detached eddy simulation resistant to ambiguous grid densities［J］.Theoretical and Computational Fluid Dynamics，2006，20（3）：181－195.

［21］JAMESON A，SCHMIDT W，TURKEL E.Numerical solutions of the Euler equations by finite volume methods using Runge－Kutta time－stepping schemes［C］.14thFluid and Plasma Dynamics Conference，1981，California，AIAA 1981－1259.

［22］SHEN Y Q，WANG B Y，ZHA G C.Implicit WENO scheme and high order viscous formulas for compressible flows［C］.25thAIAA Applied Aerodynamics Conference，June 2007，Miami，AIAA 2007－4431.