趙雅甜，邵志遠(yuǎn)，閻超，向星皓

1.中南大學(xué) 交通運(yùn)輸工程學(xué)院，長沙 410083

2.北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京 100191

3.中國空氣動力研究與發(fā)展中心，綿陽 621000

流動分離與翼型和機(jī)翼氣動特性息息相關(guān)，飛行包線也時刻受到流動分離的影響?，F(xiàn)代飛行器愈發(fā)苛刻和復(fù)雜的飛行條件使得設(shè)計過程中對分離流動的準(zhǔn)確預(yù)測愈發(fā)關(guān)鍵。然而，對于較大逆壓梯度下的分離流動，精確的CFD（Com‐putational Fluid Dynamics）模擬仍十分困難［1-2］。包含展向流動的三維邊界層分離和高速下激波誘導(dǎo)邊界層分離等典型復(fù)雜流動現(xiàn)象對預(yù)測模型的精度提出了更高要求［3-4］。

湍流模型作為RANS（Reynalds Averaged Navier-Stokes Equations）方法的核心之一，對流場計算結(jié)果影響顯著。第2 次阻力預(yù)測會議（The Second DPW）的分析結(jié)果顯示，湍流模型對CFD 計算結(jié)果的影響最大，占到約15%［5］。湍流場中各參數(shù)脈動幾乎均與幾何邊界相關(guān)，因此，當(dāng)前尚不存在某種湍流模型可適用于全部流場。被廣泛使用的渦黏性模型由于Boussinesq假設(shè)的“先天缺陷”，導(dǎo)致對分離流動預(yù)測失準(zhǔn)［6］。不同于渦黏性模型，雷諾應(yīng)力模型從雷諾應(yīng)力輸運(yùn)方程出發(fā)，對雷諾應(yīng)力直接求解，考慮了雷諾應(yīng)力的對流和擴(kuò)散和流動歷史效應(yīng)［7］。因此該模型被認(rèn)為是最自然、最符合邏輯的模型［8］。NASA（National Aeronautics and Space Admin‐istration）在CFD 2030 遠(yuǎn)景中也將雷諾應(yīng)力模型作為主要發(fā)展的RANS 方法之一［9］。基于該類模型，Rumsey 等［10］通過對亞聲速駝峰算例和ONERA M6 機(jī)翼算例進(jìn)行計算后，發(fā)現(xiàn)雷諾應(yīng)力模型存在再附點(diǎn)附近流線后彎，但預(yù)測精度優(yōu)于渦黏性模型。對于小攻角來流下的三維機(jī)翼激波位置預(yù)測較為準(zhǔn)確。王圣業(yè)等［11］利用基于雷諾應(yīng)力模型的分離渦模擬方法對 24.6°迎角鈍前緣三角翼進(jìn)行數(shù)值仿真，發(fā)現(xiàn)模擬結(jié)果明顯好于傳統(tǒng)基于線性渦黏模型的分離渦模擬方法。舒博文等［12］發(fā)現(xiàn)雷諾應(yīng)力模型可在翼身接合區(qū)角區(qū)流動和三維強(qiáng)激波誘導(dǎo)分離等問題中得到正確的流動特征。

然而，雷諾應(yīng)力模型當(dāng)前也存在諸多問題。熊莉芳等［13］研究發(fā)現(xiàn)，雷諾應(yīng)力模型對邊界條件十分敏感，對計算網(wǎng)格要求苛刻，數(shù)值穩(wěn)定性較差。為解決數(shù)值剛性和邊界條件奇性問題，Men‐ter 和Egorov［14］提出Stress-BSL 模型，兼顧了計算效率，在常規(guī)流動預(yù)測中表現(xiàn)良好，但缺乏在強(qiáng)逆壓梯度下分離流動的具體研究及其與渦黏性模型的對比。本文基于Stress-BSL 模型，分析雷諾應(yīng)力模型在二維駝峰（Hump）逆壓梯度分離、跨聲速凸塊（Bump）激波-邊界層干擾分離、M6 跨聲速機(jī)翼三維邊界層分離中的預(yù)測表現(xiàn)?；谀Ｐ退从车奈锢憩F(xiàn)象和數(shù)學(xué)構(gòu)造，對預(yù)測分離流動時出現(xiàn)誤差的原因及機(jī)制進(jìn)行分析探究。通過與渦黏性模型進(jìn)行對比，預(yù)測其性能優(yōu)勢，為未來湍流模型的改進(jìn)和優(yōu)化提供參考。

1 湍流模型與算法介紹

雷諾應(yīng)力模型直接求解雷諾應(yīng)力輸運(yùn)方程，本文選取 Stress-BSL 模型作為雷諾應(yīng)力模型的代表，將長度尺度ε與長度尺度ω進(jìn)行混合，該模型的構(gòu)造為

雷諾應(yīng)力輸運(yùn)方程為

式中：ρ為密度；為平均速度；Reij為雷諾應(yīng)力；μ為分子黏度性；μt為渦黏性系數(shù)。雷諾應(yīng)力再分配項作為輸運(yùn)方程的重中之重，該項的構(gòu)造為

式中：k為湍動能；δij為克羅內(nèi)克符號；C1、α、β、γ為模式常數(shù)；Sij為平均應(yīng)變率；Wij為角速度；aij為各向異性張量。Sij、Wij、aij3個參數(shù)的計算公式為

耗散項采用各向同性假設(shè)進(jìn)行模化，可表示為

ω輸運(yùn)方程為

本文空間離散方法針對有限體積法求解采用二階中心差分格式，針對黏性通量插值算法采用二階迎風(fēng) MUSCL（Monotone Upstreamcentred Schemes for Conservation Laws）重構(gòu)格式，無黏通量采用 HLLC（Harten-Lax-van Leer-Contact）格式。時間離散格式采用隱式LU-SGS（Lower-Upper Symmetric Gauss-Seidel）方法［15］。

2 數(shù)值計算與結(jié)果分析

2.1 二維駝峰（Hump）逆壓梯度分離

二維駝峰（Hump）模型［16-17］的前半部分上凸，流道收縮，流體在此處加速，壓強(qiáng)減小，后半部分突然下凹，流道擴(kuò)張，出現(xiàn)較強(qiáng)的逆壓梯度，邊界層發(fā)生分離。由于駝峰下游的壓力梯度接近于0，因此邊界層可以再附于下壁面，分離區(qū)位于駝峰后半部分下凹位置。

圖1 給出了Hump 模型的計算網(wǎng)格，對邊界層和駝峰后半部分進(jìn)行了網(wǎng)格加密，網(wǎng)格分辨率為1 633×433（流向×法向）。來流馬赫數(shù)為0.1，來流總溫與參考溫度的比值為1.007，來流總壓與參考壓強(qiáng)的比值為1.002；出口靜壓與參考壓強(qiáng)的比值為0.999 62。

圖1 Hump 模型計算網(wǎng)格Fig.1 Calculation grids of hump

圖2 給出了2 種模型計算的駝峰后緣流線和無量綱流向速度（U/U∞）分布，其中，U為當(dāng)前位置流速；U∞為自由來流速度；X、Y方向位置采用弦長（c）無量綱化?？梢园l(fā)現(xiàn)分離區(qū)上緣流動存在明顯的剪切，Stress-BSL 模型預(yù)測結(jié)果誤差更小，具體的流動分離和再附位置如表1。但是在再附點(diǎn)附近可以明顯地發(fā)現(xiàn)，Stress-BSL 模型預(yù)測的流線出現(xiàn)了不符合物理規(guī)律的向前彎曲，其他雷諾應(yīng)力模型也被發(fā)現(xiàn)存在類似的問題［4］。這可能是由于雷諾應(yīng)力模型在此處預(yù)測的湍流長度尺度增長率過大，導(dǎo)致剪切區(qū)域本應(yīng)向下游發(fā)展的流線受到分離區(qū)剪切影響突然增大，出現(xiàn)流線彎曲。

圖2 Hump 模型后緣無量綱流向速度與流線圖Fig.2 Nondimensional velocity and streamline behind hump

2.2 跨聲速凸塊（Bump）激波誘導(dǎo)邊界層分離

跨聲速凸塊（Bump）流動模型［18］由一個環(huán)形圓弧凸起構(gòu)成，氣體流經(jīng)凸起速度增加產(chǎn)生激波，在凸起的后半部分發(fā)生激波誘導(dǎo)分離，隨后邊界層再附并在凸起后緣產(chǎn)生分離區(qū)。圖3 給出了軸對稱跨聲速凸塊的計算網(wǎng)格，對邊界層、激波和環(huán)形凸塊后緣部分進(jìn)行了網(wǎng)格加密，網(wǎng)格分辨率為721×321（流向×法向）。為模擬軸對稱條件，需將二維平面網(wǎng)格繞X軸旋轉(zhuǎn)1°，并將橫向邊界面設(shè)置為周期邊界。來流馬赫數(shù)為0.875，來流單位雷諾數(shù)為2.76×106m?1。

圖3 Bump 模型計算網(wǎng)格Fig.3 Calculation grids of bump

圖4 Bump 模型后緣無量綱流向速度與流線圖Fig.4 Nondimensional velocity and streamline behind bump

2 種模型預(yù)測所得流場流線和無量綱流向速度（U/U∞）分布見圖 4，可以明顯地發(fā)現(xiàn) SST 模型預(yù)測的流動再附點(diǎn)位置過于靠后，分離區(qū)過長，具體的流動分離點(diǎn)和再附點(diǎn)位置見表2。2 種模型的計算結(jié)果均出現(xiàn)了流動分離點(diǎn)提前，但Stress-BSL 模型的誤差更小。對于再附點(diǎn)，SST 模型預(yù)測的結(jié)果靠后，而Stress-BSL 模型的結(jié)果則較實(shí)驗值靠前。由于凸塊模型高度較小，再附點(diǎn)附近流線曲率很小，因此前述的原因?qū)е铝骶€向前彎曲，易使流動提前再附。整體來看 Stress-BSL 模型預(yù)測效果更好。

表2 Bump 模型分離、再附點(diǎn)位置Table 2 Separation and reattachment points of bump

表3 網(wǎng)格無關(guān)性驗證Table 3 Grid independence verification

2.3 跨聲速 M6 機(jī)翼三維邊界層分離

ONERA M6 機(jī)翼是跨聲速三維激波誘導(dǎo)邊界層分離的經(jīng)典算例［19］。在跨聲速來流條件下，機(jī)翼的上表面會出現(xiàn)激波，并且由于三維邊界層存在展向流動，在靠近翼尖的區(qū)域激波增強(qiáng)并出現(xiàn)激波誘導(dǎo)分離。研究表明，2 種湍流模型對于小攻角下的流動預(yù)測均與實(shí)驗吻合良好，因此本文采用大攻角6.06°進(jìn)行強(qiáng)激波誘導(dǎo)分離的預(yù)測能力對比，來流馬赫數(shù)為0.84，來流單位雷諾數(shù)為1.46×107m?1。

圖5 給出了跨聲速三維 ONERA M6 機(jī)翼模型的計算網(wǎng)格，采用 O-H 型拓?fù)?，由于沒有側(cè)滑角，故采用半模進(jìn)行計算，并對邊界層、機(jī)翼前緣部分和機(jī)翼翼尖段進(jìn)行了網(wǎng)格加密，總共設(shè)計了粗、中、細(xì)、極細(xì)4 種網(wǎng)格，半模網(wǎng)格量分別約為133 萬、265 萬、530 萬、795 萬，網(wǎng)格細(xì)節(jié)見表 3。

圖5 M6 計算網(wǎng)格Fig.5 Calculation grids of M6

圖6 為采用 Stress-BSL 湍流模型計算的4 套 ONERA M6 機(jī)翼網(wǎng)格在6.06°攻角下20%展長位置處的壓力系數(shù)（Cp）分布計算結(jié)果與實(shí)驗值的對比結(jié)果。發(fā)現(xiàn)細(xì)網(wǎng)格和極細(xì)網(wǎng)格對激波位置的捕捉基本一致，且與實(shí)驗值誤差較小，為節(jié)省計算資源，后續(xù)均采用細(xì)網(wǎng)格進(jìn)行數(shù)值模擬。

圖6 網(wǎng)格無關(guān)性驗證Fig.6 Grid independence verification

圖7、表4 給出了2 種湍流模型的計算結(jié)果對比。由圖7（a）、圖7（b）可以發(fā)現(xiàn)，在機(jī)翼的背風(fēng)面存在明顯的“λ”型激波。對比2 種模型的預(yù)測結(jié)果，Stress-BSL 模型預(yù)測的激波結(jié)構(gòu)更為合理。對于激波誘導(dǎo)邊界層分離而言，波后高壓前傳導(dǎo)致波前邊界層增厚，流動受阻進(jìn)一步導(dǎo)致了激波的增強(qiáng)，因此流動分離的位置與激波位置幾乎重合。隨著逐漸靠近翼尖，激波增強(qiáng)，Stress-BSL 模型預(yù)測的結(jié)果始終與實(shí)驗值符合較好，而 SST 模型預(yù)測的激波位置開始前移，分離區(qū)過大（見圖7（c）），與實(shí)驗值的誤差增大。過大的分離區(qū)導(dǎo)致流動再附受阻，在機(jī)翼尾緣附近壓力系數(shù)明顯偏大（見圖8）。需要注意的是，圖8（a）、圖8（b）中Stress-BSL 模型預(yù)測的激波位置反而較SST 模型誤差偏大。這是由于風(fēng)洞實(shí)驗中翼根平面并不是對稱面，而CFD 計算中將翼根平面設(shè)置為對稱面，因此產(chǎn)生了邊界條件干擾，使得激波預(yù)測位置較實(shí)驗值后移。但由于SST 模型預(yù)測的激波位置靠前，因此導(dǎo)致了圖8（a）、圖8（b）的SST 模型預(yù)測激波位置更準(zhǔn)的巧合性結(jié)果。這種邊界條件干擾會隨著逐漸遠(yuǎn)離翼根平面而消失。Rumsey［4］研究也提到了這一點(diǎn)。

表4 升力系數(shù)與阻力系數(shù)對比Table 4 Comparison of lift and drag coefficient

圖7 不同湍流模型計算結(jié)果對比Fig.7 Comparison of simulation results with different turbulence models

圖8 不同站位壓力系數(shù)Fig.8 Pressure coefficient at different stations

另外，由于Stress-BSL 模型預(yù)測的流動分離更加準(zhǔn)確，分離區(qū)更小，背風(fēng)面后緣的壓力系數(shù)較SST 模型的預(yù)測結(jié)果更小，使得機(jī)翼上下表面間的壓強(qiáng)差更大，且壓差阻力更大，因此相較于SST 模型的預(yù)測結(jié)果，升力系數(shù)和阻力系數(shù)均更大（見表4）。

3 湍流模型性能機(jī)制

基于計算結(jié)果，本節(jié)著手從湍流模型的構(gòu)造機(jī)制出發(fā)，探明預(yù)測分離點(diǎn)提前、再附點(diǎn)滯后、分離區(qū)過大的具體原因。雷諾應(yīng)力的物理含義為流體微元表面上脈動動量輸運(yùn)的平均值，可以認(rèn)為對于湍流邊界層而言，更小的雷諾應(yīng)力意味著邊界層抵抗逆壓梯度的能力越差，越容易出現(xiàn)流動分離。因此，湍流模型對雷諾應(yīng)力的預(yù)測能力直接決定了對分離流動預(yù)測的準(zhǔn)確性。

3.1 Stress-BSL 模型

圖9 給出了Hump 模型、Bump 模型特點(diǎn)站位的Stress-BSL 模型雷諾應(yīng)力計算值與風(fēng)洞實(shí)驗值對比?？梢园l(fā)現(xiàn)，Stress-BSL 模型的預(yù)測值較實(shí)驗值偏低，導(dǎo)致了預(yù)測分離提前和再附滯后。雖然通過直接描述雷諾應(yīng)力的輸運(yùn)過程避免了引入假設(shè)的誤差，但雷諾應(yīng)力輸運(yùn)方程中難以實(shí)驗觀測和?；捻?，給建模帶來了未知性和誤差，使得 Stress-BSL 模型雖然具有優(yōu)勢，但仍存在較大誤差。分析雷諾應(yīng)力輸運(yùn)方程中的各項，可以得到 Stress-BSL 模型預(yù)測不準(zhǔn)的原因。

圖9 Stress-BSL 模型雷諾應(yīng)力計算結(jié)果與實(shí)驗值對比Fig.9 Comparison of stress-BSL and experimental Reynolds stress

1）生成項

由于雷諾應(yīng)力模型直接求解雷諾應(yīng)力輸運(yùn)方程，因此可以對生成項精確求解。但精確并不是絕對準(zhǔn)確，在求解雷諾應(yīng)力時由于模型建模誤差導(dǎo)致的計算誤差會發(fā)生累積。因而，所謂的精確是指在生成項求解過程中沒有引入新的變量和假設(shè)，不會產(chǎn)生新的誤差。

2）擴(kuò)散項

Stress-BSL 模型中擴(kuò)散項將雷諾應(yīng)力的擴(kuò)散速率與其在流場中的梯度表示為正比。De‐muren 和Sarkar［20］通過對平面通道湍流研究后發(fā)現(xiàn)：擴(kuò)散項在對數(shù)層對雷諾正應(yīng)力的計算結(jié)果影響不大，但可以對邊界層外層湍流轉(zhuǎn)化為各向同性產(chǎn)生影響。這一項并不是雷諾應(yīng)力模型預(yù)測帶有強(qiáng)烈剪切的分離流動準(zhǔn)確與否的決定性因素，絕大多數(shù)文獻(xiàn)并未對其影響進(jìn)行專門研究。

3）耗散項

湍流耗散的情況十分復(fù)雜，包括由大渦和小渦拉伸、壓強(qiáng)作用、湍流輸運(yùn)、分子擴(kuò)散造成的耗散作用，以及湍動能耗散自身的耗散作用。由于耗散作用基本發(fā)生在湍流的最小尺度中，此時的耗散作用除近壁區(qū)等特殊區(qū)域外，可以認(rèn)為是各向同性的。Stress-BSL 模型中擴(kuò)散項采用了各向同性的耗散模型［21］。由圖9 也可以發(fā)現(xiàn)，在近壁區(qū)雷諾應(yīng)力的計算值和實(shí)驗誤差很小，因此該假設(shè)基本合理。但需要注意的是，耗散項中的模式封閉參數(shù)事實(shí)上需要根據(jù)雷諾數(shù)不同進(jìn)行經(jīng)驗修正，在不同的流動區(qū)域內(nèi)取值不同，且耗散項與湍流長度尺度息息相關(guān)，對雷諾應(yīng)力模型的預(yù)測能力影響較大。

4）再分配項

雷諾應(yīng)力再分配項是脈動壓強(qiáng)和應(yīng)變率張量相關(guān)的平均值［22-24］。由于湍動能輸運(yùn)方程中并不存在再分配項，說明再分配項對于湍動能的產(chǎn)生沒有貢獻(xiàn)，而只是在湍流脈動速度的各個分量之間進(jìn)行調(diào)節(jié)。不同于擴(kuò)散作用依賴于明確的梯度，再分配項是一種能量轉(zhuǎn)移，通過將量值較大的速度分量轉(zhuǎn)移給量值較小的速度分量，以達(dá)到各方向分量的平衡。為明確其作用，給出二維湍流場的雷諾應(yīng)力輸運(yùn)方程

可以發(fā)現(xiàn)，Y方向的雷諾正應(yīng)力生成項為0。假設(shè)在均勻剪切場中初始時只有正應(yīng)力，沒有切應(yīng)力。由于流場的剪切作用，雷諾切應(yīng)力開始增長，導(dǎo)致X方向正應(yīng)力的生成項大于0，X方向正應(yīng)力開始增長。但Y方向的雷諾正應(yīng)力生成項為0，這將使得Y方向正應(yīng)力被耗散進(jìn)而導(dǎo)致切應(yīng)力生成項消失。再分配項則可以將X方向分量的一部分能量轉(zhuǎn)移給Y方向，使其保持一定強(qiáng)度，從而保持了切應(yīng)力的生成和強(qiáng)度，形成了一個維持湍流場的閉環(huán)機(jī)制。由于再分配項的重要作用，同時與生成項有相同的量級，且缺乏模化依據(jù)，因此是預(yù)測誤差的主要來源。Stress-BSL 模型將其表示為快速和慢速項之和?？焖夙椫邪骄俣忍荻龋@說明雷諾應(yīng)力再分配不是局部過程，不能僅取決于當(dāng)?shù)匾稽c(diǎn)的參數(shù)。但模型構(gòu)造卻基于同一物理位置的湍流波動量的相關(guān)性完成模型閉合，這會導(dǎo)致一部分湍流信息的丟失，使得快速項的建模失準(zhǔn)，在速度梯度越大的區(qū)域，誤差也會越大。圖9 的結(jié)果證明，在邊界層內(nèi)剪切最強(qiáng)烈（速度梯度大）的區(qū)域，雷諾應(yīng)力計算值與實(shí)驗值誤差最大。其中，圖中縱坐標(biāo)Y?Y0、Z?Z0分別為Hump 模型和Bump 模型的當(dāng)?shù)乇诿娣ㄏ蚓嚯x，并使用弦長進(jìn)行無量綱化。由于邊界層外側(cè)的速度梯度小于內(nèi)側(cè)，因此邊界層內(nèi)側(cè)的雷諾應(yīng)力預(yù)測誤差遠(yuǎn)大于外側(cè)。另外，當(dāng)前快速項的模化是基于快速畸變近似，非線性模型可認(rèn)為是利用張量函數(shù)性質(zhì)在原線性模型中添加了各項異性張量的二次項［25］，而不是在物理現(xiàn)象層面的對流動信息的補(bǔ)充。

3.2 SST 模型

相較于雷諾應(yīng)力模型，渦黏性模型對雷諾應(yīng)力的求解是基于Boussinesq 假設(shè)，通過渦黏性系數(shù)（μt）完成的，因此SST 與Stress-BSL 模型對分離流動預(yù)測失準(zhǔn)的原因大相徑庭。SST 湍流模型的具體構(gòu)造細(xì)節(jié)可參見文獻(xiàn)［26］。

1）流動分離

起初的兩方程模型（k-ε模型、k-ω模型）對于逆壓梯度的預(yù)測能力很大程度上局限于邊界層對數(shù)區(qū)域［27-29］，預(yù)測效果并不好。Menter 和Egorov［14］發(fā)現(xiàn)邊界層尾跡區(qū)域的渦黏性水平最終決定了渦黏性模型對逆壓梯度流動的預(yù)測能力，因此在構(gòu)造SST 模型時引入了 Bradshaw 假設(shè)，在渦黏性系數(shù)表達(dá)式中強(qiáng)制保證湍流剪切應(yīng)力與邊界層尾跡區(qū)的湍動能成正比，從而獲得了很好的附著流和弱逆壓梯度流動的預(yù)測結(jié)果。其中渦黏性系數(shù)表達(dá)式為

對于存在很大速度梯度剪切很強(qiáng)的邊界層內(nèi)，渦黏性系數(shù)的表達(dá)式中分母則由平均應(yīng)變率（S）代替，使得湍動能生成項保持與耗散項相等，這也對應(yīng)了Bradshaw 假設(shè)中包含的湍流能量平衡。但對比雷諾應(yīng)力計算結(jié)果與實(shí)驗值（見圖10），顯然模型對于雷諾應(yīng)力的限制過度了，導(dǎo)致雷諾應(yīng)力預(yù)測值較實(shí)際情況偏小。這是由于Bradshaw 假設(shè)是通過Klebanoff 零壓力梯度平板邊界層試驗數(shù)據(jù)總結(jié)得出的［30］，其基礎(chǔ)為無壓力梯度的附著流動，當(dāng)逆壓梯度非常大時，流動的平均應(yīng)變率很高，使得湍動能生成項明顯大于耗散項，但 SST 模型中的控制器限制其生成項仍等于耗散項，使得湍動能和渦黏性系數(shù)小于真實(shí)值，造成了預(yù)測的雷諾應(yīng)力偏小。

圖10 SST 模型雷諾應(yīng)力計算結(jié)果與實(shí)驗值對比Fig.10 Comparison of SST and experimental Reyn‐olds stress

還需要注意的是，基于二維平板邊界層數(shù)據(jù)提出的Bradshaw 假設(shè)在推廣到三維邊界層時引入了“主雷諾應(yīng)力”的概念，將三維邊界層內(nèi)雷諾應(yīng)力張量其他各向分量的變化描述為對主雷諾應(yīng)力的變化的貢獻(xiàn)。但在強(qiáng)逆壓梯度下，雷諾切應(yīng)力各方向分量（UV、VW、UW Stress）間的差值很小，雷諾正應(yīng)力（UU、VV、WW Stress）與切應(yīng)力達(dá)到同一量級（見圖11）。顯然主雷諾應(yīng)力的的概念不再成立。

圖11 雷諾應(yīng)力對比Fig.11 Reynolds stress comparison

考慮到BSL 模型與SST 模型的構(gòu)造基本相同，唯一區(qū)別在于BSL 模型未對湍動能生成項進(jìn)行限制，因此可以通過BSL 模型，對基于Brad‐shaw 假設(shè)的控制器所帶來的影響進(jìn)行對比分析。通過BSL 模型與SST 模型的計算結(jié)果對比（見圖8、圖12），發(fā)現(xiàn) BSL 模型的計算結(jié)果不同程度上抑制了流動分離的提前，證實(shí)了上述分析的SST 模型預(yù)測失準(zhǔn)的原因。

圖12 BSL 與SST 模型計算結(jié)果對比Fig.12 Comparison of BSL and SST models

究其原因，是因為渦黏性模型所基于的Boussinesq 假設(shè)包含的線性關(guān)系、各向同性等“先天缺陷”，導(dǎo)致了渦黏性模型在預(yù)測存在突然的平均應(yīng)變率變化、額外應(yīng)變率和明顯的各向異性等特征的流動時出現(xiàn)失準(zhǔn)。

2）流動再附

如上所述，SST 模型對邊界層雷諾應(yīng)力的低估導(dǎo)致了流動分離的提前。相應(yīng)的，湍流模型對分離區(qū)上緣與主流之間的雷諾應(yīng)力的預(yù)測則影響了流動的再附位置的計算結(jié)果［31］。由于雷諾應(yīng)力表征脈動動量的輸運(yùn)，因此在分離區(qū)上緣預(yù)測的雷諾應(yīng)力越大，動量輸運(yùn)越多，湍流能量越高，流動再附位置越靠前。相反地，預(yù)測的雷諾應(yīng)力越小，則流動再附位置的計算結(jié)果越靠后?？紤]到SST 模型中雷諾切應(yīng)力主要由渦黏性系數(shù)（μt）和平均應(yīng)變率（S）決定。因此希望通過對比 Stress-BSL 模型與 SST 模型的計算結(jié)果找出主要影響因素（見圖13、圖14）。

圖13 渦黏性系數(shù)對比Fig.13 Comparison of eddy viscosity coefficient

圖14 平均應(yīng)變率對比Fig.14 Comparison of mean strain rate

發(fā)現(xiàn)相較于Stress-BSL 模型，SST 模型預(yù)測的μt偏小，S則偏大，但圖9、圖10 的結(jié)果表明SST 模型預(yù)測的雷諾應(yīng)力更小。據(jù)此推測，相較于S，μt的計算結(jié)果對于雷諾應(yīng)力的預(yù)測影響更大。考察計算μt的SST 模型輸運(yùn)方程，由于雷諾應(yīng)力被低估，導(dǎo)致湍動能輸運(yùn)方程的生成項被低估，擴(kuò)散項同樣偏小，上述共同作用使得湍動能k偏小，但表達(dá)式為隱式，影響較為復(fù)雜。對于比耗散率ω，利用量綱分析，可以得到湍流特征長度尺度δl的表達(dá)式

發(fā)現(xiàn)δl與ω之間為倒數(shù)關(guān)系。由于雷諾應(yīng)力表征湍流內(nèi)的動量輸運(yùn)，因此雷諾應(yīng)力的預(yù)測與湍流能量緊密相關(guān)，而湍流能量又與湍流長度尺度直接相關(guān)。因此雷諾應(yīng)力被低估，說明湍流能量和長度尺度被低估，則ω被高估。考慮ω輸運(yùn)方程，則耗散項同樣會被高估，這就要求生成項偏高更多，由于已經(jīng)分析過湍動能生成項被低估，因此可以認(rèn)為幾乎完全由位于分母位置的μt導(dǎo)致了ω生成項被高估。這說明，雖然k、ω輸運(yùn)方程相互耦合，且均與μt相關(guān)，但顯然μt、ω的相互影響更為直接。據(jù)此推測，與ω相關(guān)的模型長度尺度的建模失準(zhǔn)、忽略耗散過程的各向異性特征是造成μt預(yù)測結(jié)果偏小的主要原因。

3）關(guān)鍵參數(shù)再標(biāo)定

基于Bradshaw 假設(shè)的限制器強(qiáng)制湍動能生成項與耗散項相等，使得SST 模型低估了雷諾應(yīng)力，但該假設(shè)豐富了模型的物理信息，因此并不能直接摒棄該假設(shè)。Zhao 等［32］研究后發(fā)現(xiàn)，系數(shù)a1對于分離激波具有很高的敏感度，并且a1是激波誘導(dǎo)分離的分離區(qū)內(nèi)壓力系數(shù)不確定性的主要來源，因此可以通過增大a1放寬對雷諾應(yīng)力的限制。經(jīng)過試算，發(fā)現(xiàn)a1=0.35 可以取得相對滿意的結(jié)果。

由圖15 可以發(fā)現(xiàn)，對a1進(jìn)行調(diào)整可以在一定程度上取得更加接近實(shí)驗數(shù)據(jù)的結(jié)果，再標(biāo)定后的模型對于M6 機(jī)翼的背風(fēng)面激波預(yù)測更加合理，兩道激波合并為一道的位置向翼尖移動，且合成后的強(qiáng)激波的位置向下游移動。由壁面流線可以發(fā)現(xiàn)再標(biāo)定后模型預(yù)測的分離區(qū)更小，沒有出現(xiàn)分離起始位置過于靠近翼根位置的現(xiàn)象，整體對激波位置的捕捉和激波后區(qū)域的壓力系數(shù)的預(yù)測均更為準(zhǔn)確。

圖15 a1標(biāo)定前后壓力對比Fig.15 Comparison of pressure with different a1

4 結(jié)論

經(jīng)過試算，對比篩選后選用了SST 模型、Stress-BSL 模型作為渦黏性模型和雷諾應(yīng)力模型的代表對湍流分離流動進(jìn)行了計算。選取了NACA4412 翼型、二維駝峰、二維跨聲速凸塊、跨聲速三維ONERA M6 機(jī)翼，從數(shù)學(xué)構(gòu)造和物理現(xiàn)象等方面入手對這2 個湍流模型的分離流動的預(yù)測性能進(jìn)行了對比，給出了較為深入的機(jī)制分析，并基于分析所得原因?qū)δＰ瓦M(jìn)行了常數(shù)再標(biāo)定，提出了改進(jìn)方向，得到如下主要結(jié)論：

1）對于強(qiáng)逆壓梯度下的分離流動，Stress-BSL 模型和SST 模型均出現(xiàn)預(yù)測失準(zhǔn)，但Stress-BSL 模型的預(yù)測結(jié)果與實(shí)驗值誤差更小、性能更好。在激波誘導(dǎo)邊界層分離的復(fù)雜流動中，Stress-BSL 模型預(yù)測的激波特征和壓力系數(shù)均與實(shí)驗值吻合較好，而SST 模型則明顯預(yù)測分離區(qū)過大。在三維邊界層分離中，物理基礎(chǔ)更加可靠的Stress-BSL 模型的計算結(jié)果與實(shí)驗值誤差很小，而SST 模型計算的壓力系數(shù)與雷諾應(yīng)力則與實(shí)驗值偏差較大。

2）Stress-BSL 模型對分離流動預(yù)測的結(jié)果存在誤差，主要是由于雷諾應(yīng)力再分配項的?；粶?zhǔn)確，其基于快速畸變近似和張量函數(shù)性質(zhì)，缺乏與實(shí)際物理過程的聯(lián)系。雷諾應(yīng)力輸運(yùn)方程中耗散項同樣復(fù)雜，影響重要但程度不如再分配項。其他如擴(kuò)散項等對雷諾應(yīng)力計算的影響不大。另外，計算中發(fā)現(xiàn)，數(shù)值穩(wěn)定性仍是雷諾應(yīng)力模型亟待解決的問題，對條件要求較為苛刻，需根據(jù)其特點(diǎn)開發(fā)相應(yīng)的數(shù)值方法。

3）SST 模型對分離流動預(yù)測的結(jié)果取決于對雷諾應(yīng)力的?；７蛛x提前是因為基于Brad‐shaw 假設(shè)的控制器，使得本應(yīng)遠(yuǎn)大于耗散項的湍動能生成項被限制，導(dǎo)致雷諾應(yīng)力預(yù)測偏低。而耗散尺度建模失準(zhǔn)導(dǎo)致的主導(dǎo)雷諾應(yīng)力預(yù)測結(jié)果的渦黏性系數(shù)被低估則是再附滯后的主要原因。

4）對于流動分離提前的問題，本文通過對湍動能生成項限制放寬后，預(yù)測結(jié)果得到改善。另外，渦黏性模型對分離流動計算不準(zhǔn)的問題很大程度上可以通過對湍流尺度的建模進(jìn)行改善。對于湍流各向異性特性的補(bǔ)充，則應(yīng)聚焦于對雷諾應(yīng)力本構(gòu)關(guān)系添加非線性項進(jìn)行修改。

5）模型的長度尺度對預(yù)測準(zhǔn)確性起到關(guān)鍵作用，可以在消除邊界奇性的基礎(chǔ)上，補(bǔ)充對長度尺度的控制，以獲得更好的預(yù)測結(jié)果。