舒博文，杜一鳴，高正紅，夏露，陳樹生,4,*

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

2. 中國空氣動力研究與發(fā)展中心 空天技術研究所，綿陽 621000

3. 沈陽航空航天大學 航空宇航學院，沈陽 110000

4. 航空工業(yè)第一飛機設計研究院，西安 710089

飛行器設計中，邊界層分離對機翼失速、跨聲速抖振等特性都有著重要影響，同時，增升裝置、翼梢小翼等設計中都涉及到分離流動。隨著飛行器性能的不斷提升，飛行包線更大，所面臨的流動問題更為復雜，對于非設計點下的氣動性能更加關注。由此可見，準確預測分離剪切流動對飛行器設計十分關鍵。NASA CFD 2030遠景研究[1]中列舉了一些模擬成熟度不高的關鍵流動現(xiàn)象，“大迎角分離”、“激波/邊界層干擾” 以及“角區(qū)/接合區(qū)二次流動”等典型工程問題都在其中，這3類流動物理機理各不相同。低速大攻角分離是主要由逆壓梯度導致的分離，這類流動的模擬對增升裝置以及飛行器非設計點狀態(tài)評估十分重要。激波/邊界層分離則是來自于激波/邊界層干擾，跨聲速飛行時，機翼翼梢處極有可能出現(xiàn)明顯的激波誘導分離。準確的預測激波/邊界層擾，對飛行器跨聲速嗡鳴、抖振的預測有關鍵作用。二次流動則來源于雷諾正應力的各向異性，在航空工程中，這類流動的代表是角區(qū)/接合區(qū)流動(Corner/Junction flow)，它廣泛存在于許多外流(例如翼-身、短艙-掛架接合區(qū)等)和內流(例如發(fā)動機進氣道等)問題中[2]。如何更準確地預測邊界層分離仍是目前CFD應用面臨的主要挑戰(zhàn)之一。

基于雷諾平均納維-斯托克斯方程(Reynolds Averaged Navier-Stokes equations，RANS)的湍流模型經歷了長足的發(fā)展，以渦黏性假設為基礎的渦黏模型(Eddy Viscosity Model, EVM)已在中等壓力梯度下的附著流動中取得了巨大的成功，但其對分離和剪切流動的模擬能力仍不足。一般來說，EVM只有對二維剪切主導，正應力不重要的簡單運動才能給出可靠的結果，對于分離等復雜流動則缺陷明顯[3]。隨著計算機的快速發(fā)展，直接數(shù)值模擬(Direct Numerical Simulation, DNS)、大渦模擬(Large Eddy Simulation, LES)成為了研究湍流問題的重要工具。然而，受制于龐大的計算量，這些方法主要用于基礎研究[4]。隨著計算能力的提升，LES方法已經在很多工程問題中得到應用，但距離完全取代RANS方法還存在一定距離[1]。近年來快速發(fā)展的混合類方法(Hybrid LES/RANS)仍要依賴RANS方法對近壁面流動進行模擬。因此，基于RANS的湍流模型因其高效、魯棒性強的特點，仍將是長時間內工業(yè)界應用的主流[1,3-5]。

區(qū)別于渦黏性湍流模型，雷諾應力模型(Reynolds Stress Model, RSM)作為“最完整的”RANS湍流模型[1]，直接求解雷諾應力的輸運方程，方程中自然的包含了流線彎曲、雷諾應力各向異性等非線性現(xiàn)象影響，理論上能夠更全面地反映復雜流動特性，對相關的流動進行模擬計算更具有優(yōu)勢。自雷諾應力輸運方程于1945年首次被周培源[6]提出,相繼出現(xiàn)了著名的LRR(Launder, Reece, Rodi)模型[7]、SSG(Speziale, Sarkar, Gatski)模型[8]。由于RSM要求解雷諾應力及湍流長度尺度共7個方程，很長時間內沒有受到重視。隨著計算能力的提升，RSM又進入了人們的視野。德國宇航院的Eisfeld等[9]類比Menter[10]提出的剪應力輸運(Shear Stress Transport, SST)模型的思想，將LRR模型和SSG模型相結合，提出了SSG-LRR-ω模型，并在一些問題中取得了較好的效果[11-14]。筆者所在團隊[15]討論了微分雷諾應力模型在激波分離流中的應用；董義道等[16]開展了SSG-LRR-ω模型的初步應用研究；王圣業(yè)等[17]將DES中的RANS部分選擇為SSG-LRR-ω模型，對基于雷諾應力模型的DES方法進行了研究。此外，美國國家航空航天局(National Aeronautics and Space Administration,NASA)在CFD2030遠景研究[1]中也將雷諾應力模型作為主要發(fā)展的RANS方法之一。為了解決ω在壁面缺少自然邊界條件的問題，Togiti和Eisfeld[18]提出SSG/LRR-g雷諾應力模型。雖然SSG/LRR-g在經典流動中取得了很好的效果，但對于該模型在典型航空分離流動中的具體表現(xiàn)以及與EVM的對比并沒有系統(tǒng)研究。因此，基于SSG/LRR-g模型，以NACA4412翼型大攻角分離、M6機翼大攻角激波/邊界層分離以及F6翼身組合體角區(qū)分離為例，通過對典型狀態(tài)計算分析，探討雷諾應力模型對逆壓梯度、激波誘導分離以及雷諾應力各項異性誘導分離流動的預測能力，以確認該模型處理常見航空分離流動的能力，為進一步改進模型奠定基礎，并通過與試驗以及渦黏模型的對比，顯示了RSM相對于常用線性渦黏模型處理航空分離流動問題的優(yōu)勢。

1 雷諾應力模型及數(shù)值方法

對Navier-Stokes方程進行Favre平均后會得到雷諾應力項，湍流脈動量對平均流場的作用即通過雷諾應力來體現(xiàn)。為了封閉方程需要對雷諾應力建模處理，即湍流模型。目前常用的渦黏模型基于Boussinesq假設，雷諾應力和當?shù)仄骄鶓兟食示€性變化且各向同性，沒有著重考慮流動歷史效應。雷諾應力模型則直接對雷諾應力的輸運方程進行建模，模型的湍流生成項是準確的，沒有進行模式化；雷諾應力的再分布效應、耗散效應以及擴散效應均進行建模，因此模型本身反映了流動的歷史效應、流線彎曲效應、應變率突變、雷諾應力各向異性等[19]。所有雷諾應力的值均依賴于流動初始條件和發(fā)展過程，是一種更完整的RANS方法。下面對本文采用的SSG/LRR-g模型進行介紹。

1.1 SSG/LRR-g雷諾應力模型

采用Favre平均，由動量方程推導可得到雷諾應力的輸運方程：

(1)

(2)

與湍動能k的輸運方程中[19]生成項只考慮應變率張量Sij不同，式(2)中包含了應變率張量Sij以及渦量Ωij：

(3)

耗散發(fā)生在最小的尺度，絕大部分模型都使用Kolmogorov的局部各向同性湍流假設模化耗散項：

(4)

(5)

湍動能耗散率ε的輸運方程十分復雜，難以實現(xiàn)逐項的?；?，對于擴散項的?；捎锰荻葦U散的思想，認為擴散導致的雷諾應力輸運和雷諾應力梯度成比例。根據模化形式不同，可分為簡單梯度擴散(Simple Gradient Diffusion Hypothesis, SGDH)和廣義梯度擴散(General Gradient Diffusion Hypothesis, GGDH)，采用廣義梯度擴散：

(6)

圖1 混合模型示意圖

(7)

表1 自由參數(shù)取值

(8)

(9)

為了封閉模型，還需要求解湍流長度尺度，目前廣泛使用的是湍流比耗散率ω的輸運方程，其邊界條件如式(10)所示。由于ω缺乏自然邊界條件，因此，采用Lakshmpitathy提出的ω的變形式[21]來處理,變形后的長度尺度為g，如式(11)所示:

(10)

(11)

式中:Δy表示物面法向首層網格距離;ν表示黏性系數(shù);β為SST模型中的系數(shù);D表示耗散系數(shù)。

1.2 數(shù)值方法

使用結構網格求解器，采用有限體積法求解三維可壓縮RANS方程組。其中，RANS方程無黏通量項采用通量差分分裂法(Flux Difference Splitting, FDS)，重構格式為MUSCL(Monotone Upstream Centered Scheme for Conservation Laws)格式；黏性通量采用二階中心差分格式離散。湍流模型方程的對流項統(tǒng)一采用一階差分格式。由于本文處理分離問題，因此使用全量Navier-Stokes方程組。時間推進采用隱式近似因子分解(Approximate Factorization method, AF)法。對于二維問題，將網格沿展向外推一層并設置對稱邊界條件，以滿足求解三維RANS方程組需求。

2 典型分離問題研究

2.1 NACA4412翼型后緣分離

NACA4412翼型在13.87°攻角時(接近失速攻角)的流動常被作為簡單高升力流動的驗證算例[12]。Coles和Wadcock[22]完成的試驗表明在翼型上表面后緣，存在一定的分離區(qū)。計算網格為NASA湍流模型網站[23]提供的C型拓撲密網格，網格分布為897×257，如圖2所示，計算狀態(tài)為：弦長c=1、Ma=0.09、基于弦長的雷諾數(shù)Rec=1.52×106、α=13.87°。

圖2 NACA4412計算網格

圖3展示了2種不同湍流模型的計算所得翼型表面壓力分布系數(shù)與試驗所得結果的對比。其中，橫坐標x/c表示用弦長c歸一化的翼型橫坐標?？梢钥吹?，采用2類模型的計算結果很相似，與試驗的差距主要在后緣分離泡處。Eisfeld等[12]研究表明SST模型預測結果通常分離點提前，再附滯后，導致分離泡長度過長。導致2個模型分離預測明顯區(qū)別的原因在于對于雷諾應力的預測不同，隨后討論。

圖3 壓力分布對比

圖4 測量站位示意圖

圖5 不同站位速度型及雷諾應力分布對比

可以看到，雖然2種模型均沒有準確地預測雷諾應力，但對于速度型等宏觀物理量都能取得較為滿意的結果。這是由于雷諾應力在黏性底層(y+≈5)附近遠小于分子黏性，使得邊界層平均速度型的漸進特性并不依賴于雷諾應力等湍流變量的漸進特性。因此，雖然雷諾應力與試驗有明顯區(qū)別，但模型仍然可以對平均流動速度型和宏觀力學特性進行相對正確的模擬[10]。

對于此類主要由逆壓梯度引起的分離，其準確預測的核心是雷諾應力的預測。RSM直接求解雷諾應力的輸運方程，相比于EVM用渦黏性代替雷諾應力，RSM能夠得到更準確的雷諾應力，因此對分離位置及分離區(qū)長度會更為準確。但是，由于RSM中關鍵項：壓力應變關聯(lián)項以及耗散項仍需要建模處理，這嚴重影響了RSM中雷諾應力的求解精確度。由于RSM得到的雷諾應力更為準確，其得到的速度型更加反映出了逆壓梯度的影響。但2個模型都低估了雷諾應力，這可能是RANS方法無法捕捉自由剪切層中的非定常湍流運動，這同時也影響著分離泡的預測精度。

2.2 M6機翼跨聲速激波誘導分離

M6機翼繞流是測試湍流模型處理跨聲速流動問題能力的經典算例，由Schmitt[24]在ONERA S2MA風洞完成試驗。M6機翼前緣后掠角30°，半展長1.196 3 m，具體的幾何參數(shù)如圖6[24]所示。激波給邊界層施加了一個強逆壓梯度，使邊界層變厚，產生分離。一方面，邊界層承受著激波產生的強逆壓梯度影響，另一方面，激波必須穿過無黏和有黏的多層流動結構。若流動不是層流的，則湍流被強化，黏性耗散被放大[25]。

圖6 M6機翼幾何參數(shù)[24]

大部分學者[9,11,18]研究了3.06°攻角和4.08°攻角激波誘導分離不強的流動，渦黏模型均能取得較好的結果。強激波誘導分離問題代表了跨聲速機翼表面的一類重要流動現(xiàn)象，但鮮有渦黏模型能夠準確預測[15]。因此，選擇試驗編號為Test 2565的6.06°攻角工況進行研究，討論雷諾應力模型對強激波誘導分離問題的模擬能力。設計了極粗、粗、中等、密網格4套網格，分別對應網格量為：74萬、134萬、330萬和510萬，4套網格均保證y+≈1，如圖7所示。計算馬赫數(shù)Ma=0.84，基于平均氣動弦長的雷諾數(shù)Remac=11.7×106，攻角6.06°。不同網格量預測的y/b=0.44(b為M6機翼展長)站位壓力分布如圖8所示，中等和密網格捕捉的激波位置和強度一致，且和試驗符合較好，之后的計算均采用密網格。

圖7 M6機翼表面網格

圖8 不同網格預測y/b=0.44站位壓力分布

流動沿展向發(fā)展，激波逐漸由兩道匯聚成一道很強的激波，并引發(fā)激波誘導分離，在翼梢處達到最強，這一點可以直觀地從圖9得到。

圖9 SST模型壓力分布沿展向變化

從圖10的壓力分布對比中可以看到，SST模型在y/b=0.44站位處得到的激波位置很準確，但RSM模型預測的激波位置則靠后；從y/b=0.65站位起，SST模型對激波后的分離區(qū)預測結果出現(xiàn)明顯偏差，進而影響機翼后緣及下表面壓力分布。RSM模型預測的激波位置隨仍滯后于試驗，但對于壓力分布的形態(tài)、激波后分離區(qū)的預測，均比SST有巨大提升，能夠正確地預測流場。

圖10 不同站位壓力分布對比

不同模型得到的物面極限流線如圖11所示，可以清楚地看到M6機翼上表面的 型激波。RSM模型預測的分離區(qū)明顯小于SST模型。

圖11 極限流線對比

為了進一步確認大攻角時渦黏模型對激波誘導分離預測的失效原因，研究了隨網格加密不同模型的預測能力變化，如圖12所示。其中，SST模型在y/b=0.44站位處隨著網格加密，壓力分布和試驗結果的誤差減小，但在y/b=0.90站位處，壓力分布隨網格加密逐漸偏離試驗值，這與Rumsey和Vatsa[26]的結論一致。RSM則隨著網格加密,預測的激波更加陡峭，和試驗結果符合較好,表現(xiàn)出了應有的網格收斂趨勢。

圖12 典型站位網格收斂性對比

SST模型中關鍵的一項是引入的“應力限制器”:

(12)

由于已經處于SST預測的分離區(qū)內，因此SST模型計算得到的雷諾應力分量都很大?？梢钥吹秸瓜蚶字Z應力增長明顯；除流向-法向剪應力外，其余2個剪應力分量也增長明顯。三者基本具有相同量級，均不能被忽略。這也證明了在三維激波/邊界層分離流動中，對各向雷諾應力都需要比較準確(量級和相對大小至少要正確)的評估才能得到合理的分離和壓力分布計算結果。RSM并沒有“主雷諾應力”的概念，雷諾應力求解完全依賴于輸運方程，因此得到的雷諾應力能夠更為準確地反應流動特征，從而得到較為合理的結果。因此，準確預測激波/邊界層分離問題的核心首先在于湍流模型，而非數(shù)值格式和網格密度。

2.3 F6翼身接合處角區(qū)分離二次流動

翼身連接處流動十分復雜，除了有繞機翼前緣的馬蹄渦外，還存在雷諾應力梯度導致的二次角渦[27]，如圖14所示?，F(xiàn)有湍流模型對于此類物體連接處分離渦的預測常被認為是不可信的[28]。近年來，更精確的雷諾應力模型是預測分離接縫流動的最低要求正逐漸成為共識[29]，雷諾應力控制了近角區(qū)處應力誘導渦的發(fā)展，雷諾正應力間的區(qū)別是此類二次角渦的驅動力。壁面附近流動會在湍流非均勻性和各向異性作用下產生2個對轉的應力誘導渦(Stress-induced vortices)，在與馬蹄渦和角渦的相互作用下，與角渦轉向相同的應力誘導渦將與角渦融合，最終形成二次流渦系結構[30]。

圖14 平板-機翼接合構型流動現(xiàn)象[27]

為研究2類湍流模型對此類流動預測結果的影響，本節(jié)選取DPW-II[31]未經整流設計的DLR-F6機翼-機身-掛架-短艙(Wing-Body-Pylon-Nacell, WBPN)構型,其前緣后掠角27.1°，試驗在ONERA S2MA風洞中完成。風洞油流試驗如圖15 所示，機翼上表面翼根后緣處存在一個清晰的角區(qū)分離。此外，外翼段后緣有一道清晰的分離線,短艙-掛架上的流動向后掠方向偏轉(圖15(a))，而在機翼下表面掛架兩側也會出現(xiàn)旋渦流動和明顯的流線偏轉(圖15(b))。計算狀態(tài)為Ma=0.75，Remac=3×106,CL=0.5的巡航狀態(tài)。網格采用組委會提供的中等網格，如圖16 所示，網格量830萬。

圖15 DLR-F6 WBPN構型巡航狀態(tài)風洞油流流動顯示[31]

圖16 F6-WBPN中等網格

RSM和SST模型得到的機翼上下表面及短艙流動形態(tài)對比如圖17所示。能夠看到RSM模型預測的翼根角區(qū)分離形態(tài)和試驗符合較好,而SST得到的角區(qū)分離則明顯比試驗大，這可從圖18中更直觀地看到。RSM對機身分離的抑制作用強，得到的短艙后方機翼上表面的流動偏轉也有所減弱，符合試驗觀測。表2給出了2個模型和試驗氣動力系數(shù)的對比，由于RSM預測的角區(qū)分離較低，因此相同升力系數(shù)對應的攻角比SST模型小，另外RSM對短艙掛架上方機翼流動偏轉也有一定影響，綜合作用下得到的阻力系數(shù)和試驗符合更好。

表2 氣動力系數(shù)對比

圖17 機翼上下表面及短艙流動形態(tài)對比

圖18 翼根后緣分離形態(tài)對比

截取了翼根處y/b=0.150和y/b=0.239截面的壓力分布，如圖19所示。y/b=0.239截面處已離開分離區(qū)，因此RSM和SST模型預測的結果基本一致，和試驗符合的較好。但在y/b=0.150截面處，2個模型的計算結果從x/c=0.4開始出現(xiàn)明顯區(qū)別。SST模型由于得到的分離區(qū)過大，后緣壓力損失較大，同時截面中部的壓力分布也和試驗有一定差距。RSM在截面中段的壓力分布和試驗符合良好，后緣壓力分布較SST模型有明顯改進，這也說明后緣分離是流動流向發(fā)展的結果，特別是與x/c=0.4后的流動特性有關。

圖19 翼根區(qū)域展向截面壓力分布對比

進一步對翼身接合區(qū)的流動發(fā)展特性進行分析，圖20展示了不同模型得到的翼身接合區(qū)的流向渦量分量演化。可以看到，對翼根分離有顯著減弱作用的RSM模型在角區(qū)混合邊界層處模擬出了較大的流向渦量。這個較大的渦量使得角區(qū)流動從勢流區(qū)獲得了更多的能量，從而減小分離。此時雖然渦黏系數(shù)較小，但在非線性雷諾應力本構關系的作用下，雷諾應力各向異性能夠被準確模擬。法向-展向雷諾剪應力將顯著增強，進而促進應力誘導渦的強度，增強應力誘導渦的強度，強化了角渦的能量交換并抑制了分離。

圖20 翼身接合區(qū)的流向渦量分量演化對比

渦黏模型的基礎是Boussinesq假設，這意味著雷諾應力張量和平均應變率張量是線性關系，雷諾應力各向同性。在平衡態(tài)附著流動中EVM表現(xiàn)良好，而在此類二次流動中，EVM無法預測雷諾應力各向異性以及流動的歷史效應，導致了很大的誤差。RSM求解雷諾應力輸運方程，自然包含了雷諾應力的對流和擴散，進而考慮了流動的歷史效應。此外，對流項和生成項的精確求解則反映了流線彎曲、旋轉等現(xiàn)象。即使在平均應變率張量為零的流動中，RSM也沒有雷諾正應力相等的判斷，模型中的雷諾應力的值完全依賴于初始條件和流動發(fā)展過程。RSM的上述優(yōu)點使其特別適合于此類二次流動。

綜合來看，RSM能夠在典型航空分離流動中取得比渦黏模型更好的結果，具有一定的優(yōu)勢，但仍有改進的空間。首先是耗散項ε的建模：耗散項在近壁面是高度各向異性的，而目前包括本文采用的雷諾應力模型在內的大部分模型均采用各向同性假設來建模，影響到雷諾應力沿法向分布，這也是NACA4412分離流動中雷諾正應力預測與試驗偏離較大的原因之一。因此，應關注耗散率的各向異性建模。其次，RANS框架下封閉模型需要求解長度尺度方程，目前常用的是比耗散率ω的輸運方程，輸運方程基于量綱分析，并不像ε方程精確推導得到。采用的RSM模型也是基于ω長度尺度方程。長度尺度方程對邊界層對數(shù)區(qū)計算結果影響很大，也直接決定著分離的預測精度。進一步研究應關注長度尺度方程的建模，例如k-kl湍流模型中出現(xiàn)了von-Karmen長度尺度，可以自動調整湍流模擬尺度滿足復雜流動的模擬需求，且自動滿足壁面對數(shù)率，在分離流動中體現(xiàn)了很好的結果。將RSM中的長度尺度方程替換為更為準確且便于實現(xiàn)的長度尺度方程也是今后的發(fā)展方向。

對于EVM，從與RSM的對比中可以看出其主要缺點在于缺乏反應雷諾應力各向異性的源項。因此，對于EVM的改進應著眼于添加反應雷諾應力各向異性的修正。其中，繼續(xù)完善改進二次本構關系以及通過DNS、LES數(shù)據，基于機器學習[32-33]構造雷諾應力修正項是值得發(fā)展的方向。

3 結 論

采用SSG/LRR-g雷諾應力模型，選取NACA4412 低速大攻角流動、M6跨聲速分離和F6翼身接合處角區(qū)二次流動進行數(shù)值模擬，并和常用的渦黏模型進行對比，討論了RSM對典型航空外流分離流動的預測能力，得到了以下結論：

1) 對于NACA4412后緣分離代表的逆壓梯度主導分離流動，SST模型預測的分離泡長度更長，由于RSM模型預測的雷諾應力更大，得到了更小的分離泡。RSM能夠更為準確地反映逆壓梯度的影響，得到的速度型分布和雷諾應力分布和試驗符合的更好。

2) 對于M6機翼跨聲速強激波誘導分離流動，RSM在激波較弱的區(qū)域得到的激波位置略滯后于試驗，但在激波誘導分離很強的區(qū)域得到的激波位置和分離區(qū)均和試驗結果符合較好，更為合理的雷諾應力分布得到了更準確的激波后分離區(qū)流動，從而更準確地模擬整個流動。在此類流動中，“主雷諾應力分量”的概念不再存在，因此SST模型過度限制了雷諾應力，錯誤地預測了大范圍的分離區(qū)，進而導致預測激波位置十分靠前。隨著網格加密，RSM呈現(xiàn)出了網格收斂的結果，而SST模型的結果則愈發(fā)偏離試驗值。

3) 對于F6翼身接合處角區(qū)分離流動，RSM得到的分離范圍和試驗結果符合較好，SST模型過度估計了分離區(qū)大小。RSM能夠更加準確地預測雷諾應力各向異性，而這正是渦黏模型的缺點之一，因此雷諾應力模型對于此類角區(qū)分離流動比渦黏模型有明顯的優(yōu)勢。

4) SSG/LRR-g模型中耗散項的模型化仍采用各向同性假設，降低了壁面附近雷諾應力預測精度，今后應開展各向異性耗散率建模的研究。

5) 長度尺度方程是分離流動預測的關鍵因素，目前SSG/LRR-g模型仍基于傳統(tǒng)兩方程模型的長度尺度輸運方程。今后應將新的長度尺度方程引入RSM模型建模當中，k-kl模型中的von-Karmen長度尺度是很好的發(fā)展方向。

6) 計算效率和魯棒性仍是制約RSM發(fā)展的關鍵因素。在開展RSM各項建模的研究基礎上，也應該大力開展相關計算方法研究，提升計算魯棒性與效率。

7) 針對EVM的改進，除了改進長度尺度方程外，重點在于添加反應雷諾應力各向異性的修正?；诙伪緲嬯P系以及機器學習方法構建修正項是值得關注的方向。