微分雷諾應力模型在激波分離流中的應用

聶勝陽，高正紅，黃江濤

（西北工業大學 翼型葉柵空氣動力國防科技重點實驗室，陜西 西安710072）

0 引 言

基于Boussinesq渦粘性假設的湍流模型在很多工程流動類型上有很好的預測能力，但是對一些特殊的流動，它會給出與實驗較大的差別。平均應變率的突變導致這些湍流模型失效并不奇怪，雷諾應力的相應變化是與自由來流進程和時間尺度不相關的變化率衡量，所以導致Boussinesq渦粘性假設失效。同樣，當流動中因快速膨脹，劇烈變形，或者大的流線彎曲產生的額外應變率，都會產生不平衡的法向雷諾應力，渦粘性假設就失效了。

為了修正這些問題，在渦粘性模型的基礎上加入曲率修正［1］，應用于存在大的流線彎曲的情形如旋渦流動。另外一種更近似描述雷諾應力張量而不引入附加方程的方法是假定Boussinesq近似是雷諾應力級數展開項的主導項，這樣流動中引入各向異性的雷諾應力分量，這就是代數應力模型（ASM）。這種模型在流體各向異性特性顯著的流動中，預測較好，但它還是基于兩方程的渦粘性假設的模型，和它們的缺點一樣［2］。

本文將完全雷諾應力方程的模式化過程進行了詳細闡述，應用一系列模式將其封閉，得到微分雷諾應力模型（DRSM）；將其與廣泛使用的幾種湍流模型應用于較大攻角的跨聲速ONERA M6機翼的計算，表明其在激波分離流動中具有更好的預測能力；對各種模型的不足之處進行了討論，為模型的改進工作提供參考。

1 不可壓雷諾應力方程

雷諾應力的差分方程完整地描述了雷諾應力張量，τij≡－，其不可壓的形式為：

其中：

稱為擴散項；由式（1）可以看出，第一，方程中已經包含了τij的對流和擴散，所以二階矩模型包含流動發展歷程的影響。耗散項和湍流輸運項中的時間尺度是和平均流動時間尺度不相關的，所以二階矩模型的上游歷程的影響要比兩方程湍流模型中人為添加的上游歷程的影響更真實；第二，式（1）包含了對流項、產生項和體積力項，這些項對流線彎曲，流體旋轉等起主要作用。所以二階矩模型如果模化的好，也能更真實地包含這些產生不平衡雷諾應力的作用；第三，式（1）并沒有反應出法向應力相等，即使平均應變率為0，也就是說，它的值由初始條件和流動進程決定。因此，二階矩模型能夠處理平均應變率有突變的情況。

2 模式假設

為了封閉方程（1），必須對耗散項張量，再分布項張量，擴散項張量進行模化。耗散項：

因為耗散作用發生在最小尺度的渦，因此耗散項的建模就相對簡單，大家都采用Kolmogorov（1941）提出的當地各向同性假設，即：

這里

標量ε就是湍動能方程中的耗散率。而對耗散率的建模可參考k－ε兩方程模型中的耗散率方程的處理方法。

由于實際中的耗散率是各向異性的，尤其是靠近壁面處，所以為引入這一特性，后人進一步修正耗散項為：

其中bij為無量綱化的雷諾應力各向異性張量，即

fs是經驗性的添加的衰減函數，為

ReT為定義的湍流雷諾數，為

擴散項：

擴散項中包含壓力脈動項和三階速度脈動項，表示脈動壓力和雷諾應力對能量的輸運作用。現有的實驗數據無法找到有用的信息去模型化壓力脈動關聯項。最簡單而又常見的模化擴散項的方法是假定一個梯度輸運進程，并且引入分子粘性作用。這里使用Launder，Reece和 Rodi（1975）［2］提出來的一種近似方法：

這里他們建議Cs≈0.11。這也是渦粘性假設的湍流模型處理擴散項常用的方法。再分布項：

壓力張量再分布項也叫壓力張量關聯項，是最難模式化的部分。第一，它和產生項是同階的；第二，它包含大量不可測的關聯項，模式假設缺乏足夠的理論依據及實驗數據。再分布項有兩種主流的模型化方法，分別是Launder，Reece和 Rodi（1975）等提出的LRR壓應力模型以及Speziale，Sarkar和Gatski（1991）等提出的SSG壓應力模型，具體參考文獻［3］。本文選取較簡單的SSG模型來模化：

這里bij定義如式（5）。

Speziale，Sarkar和Gatski給出的模式參數為：

Willin ＆Johansson［4］給出了一組較簡單的模式參數為：

3 可實現的DRSM

3.1 引入壓縮性修正

可壓流動的特征是密度的顯著變化。在馬赫數小于5（另一說為8）時［2］，壁面流動受壓縮性影響很小，因為通過激波，在小距離內壓力梯度變化不劇烈。Tw／Te在邊界層內壓縮性對渦的作用也不明顯。基于這樣的事實，Morkovin（1962）認為湍流的密度脈動量相對與自由來流是個微量，因此在湍流模式中使用不可壓的形式［2］。但 Morkovin假設在即使不是高超音速情形下也有其局限性，因為密度脈動量相對自由來流的密度并不是小量。所以壓縮性效應明顯的流動情形，應將密度脈動量考慮進去，一種有效的數學處理方法就是采用Favre（1965）的質量平均速度的概念，定義為：

某時刻瞬時速度為：

這樣在湍流模型中也就反映了密度的脈動。這里直接給出Favre平均的雷諾應力方程，具體推導過程見參考文獻［2］：

其中：

產生項：Pij＝－

再分布項：Πij＝

耗散項：εij＝

擴散項：Dij＝－

壓力作用項：PWw＝

對比可壓和不可壓的完全雷諾應力方程，可以發現，方程多了壓力作用項，這個由壓縮性引起的項在不可壓流動中是零。其他項的物理意義與不可壓方程相同，但在形式上由于引入質量平均概念而有所不同，對再分布項，耗散項和擴散項的建模思路和不可壓的二階矩模型一樣，自由來流的速度項均需采用式（12）中定義的質量平均速度。

由于壓力作用項是脈動量和平均量的耦合，還沒有合適的模型去表征它的特點，并且單個脈動量在時間平均中是可以忽略掉的，本文去掉壓力作用項。

3.2 耗散項εij的改進

前面介紹二階矩的模型化過程，但沒有給出耗散率ε的方程。這是因為，在壁面上，耗散率為ε＝0，所以在壁面上，式（1.7）、（1.8）中均出現奇點。在計算中，ε會在壁面附近劇烈變化，因此在壁面附近需將網格劃分的足夠的小。因此本文采用Wilcox（1988）等［5］提出的比耗散率ω的輸運方程來模化耗散項。輸運方程為：

模式參數為：

耗散率補充方程為：

由前述的式子就將SSG－ε模型發展成了SSG－ω模型。

4 模型驗證

選用ONERA M6機翼，計算狀態為Ma＝0.84，α＝6.06，Re＝11.7×106。這是一個跨音速大迎角狀態，實驗數據來自參考文獻［6］。

本文計算采用的離散格式為中心格式添加適當的人工粘性，對時間推進采用顯式三步Runge－Kutta推進。采用多重網格和隱式殘差光順等措施加速收斂。計算網格為非結構網格，粘性層用25層三棱柱，第一層的高度為1e－6，外推比例為1.3。為了對不同湍流模型計算效果進行比較，計算過程采用同一套分布合理的網格，使用相同的時間推進手段和CFL數，自由來流的湍流度和湍流粘性比相同。

圖1和圖2給出了表面極限流線計算結果，可以看出M6機翼在該狀態時，出現 形激波，在翼展60%處，激波匯合，形成更強的激波，波后由于強的逆壓梯度，出現大范圍的分離渦。比較發現，SA模型和SST模型的分離區域最大，在機翼中段激波后就出現了輕度分離。匯合后的激波靠近前緣，尤其是SA模型。SST模型預測出來的激波不平直。而基于Wilcox提出的標準k－ω兩方程，以及基于標準k－ω兩方程的SSG－ω模型和Willin ＆Johansson提出的改進SSG－ω模型預測出相似的極限流線分布和激波形狀和位置，W＆J改進的SSG－ω模型預測出的分離區域，流線更紊亂，標準k－ω模型預測出的激波位置更靠近下游，激波后的流動更光順。

從壓力分布與實驗數據的對比上，80%展向位置之前，SST模型給出最接近實驗的激波位置，但波后的分離較大，整個后半區域都是低壓區；在80%展向外，SA和SST模型對大的分離流動和激波位置預測失真，而基于SSG壓力關聯項假設的DRSM預測的激波位置最接近實驗，波后的分離流動雖然不是很貼合，但也可以接受。而 Wilcox的標準k－ω模型給出的激波位置靠后。

圖1 計算網格、表面壓力系數分布和表面極限流線圖（SA、SST、k－ω、SSG＋ω和W＆J SSG＋ω）Fig.1 Surface and symmetry grid ＆surface pressure coefficient distributions ＆surface skin friction lines（SA、SST、k－ω、SSG＋ω ＆ W＆J SSG＋ω）

圖2 等壓線分布圖Fig.2 Iso－surface pressure coefficient lines distributions

圖3 各站位壓力分布對比Fig.3 Pressure coefficient distributions at different standings

SA和SST都過度地預測了激波后的分離流動，這可能是在這種出現強的激波附面層干擾和強逆壓梯度時，兩種模型都沒有很好地限制空間湍流粘性系數的非物理性增加，而基于k－ω模型的計算在激波位置的預測上顯得滯后，這與附面層內湍流粘性系數的過度增長有關，但在激波后，又很好地限制了這種增長。而DRSM模型直接對雷諾應力分量建模，在各向同性流體中就退化成了Wilcox標準k－ω模型。

可以總結出，DRSM模型較一方程和兩方程模型，如本文第一部分所述，對快速膨脹，流線彎曲等流體各向異性性質顯著的流動預測更加合理。Willin＆Johanson的改進參數在激波位置的預測上稍好，但改進余地有限。

在流線比較規則的區域，壓力梯度比較和緩，如機翼80%展向之前，SA模型和Menter的SST模型對激波位置的預測比本文使用的其他湍流模型更好，尤其是SST模型。單獨比較Menter的SST模型和Wilcox的標準k－ω模型，SST模型是在近壁處采用Wilcox的標準k－ω模型，在邊界層邊緣和自由剪切層采用k－ε模型（k－ω形式），其間通過一個混合函數過渡。這就是Menter的BSL模型改善對自由來流ω值的敏感性的思路。其應用表現和Wilcox的標準k－ω模型類似。故通過修改ω的輸運方程達到改進DRSM的目的是不可行的。數學上分析兩方程模型在有逆壓梯度流動中的表現很大程度上被限定在對數區。盡管模型在對數區的模擬能力很重要，尤其是有強的逆壓梯度，但尾流區的渦粘性系數卻決定了渦粘性模型對強逆壓梯度的預測能力。為此Menter參照Johanson ＆King模型，考慮了湍流剪切應力的輸運而修改了渦粘性的表達式，使模型對強逆壓梯度流動的預測更加合理（具體見參考文獻［7］），對太強的逆壓梯度分離的預測又顯得太敏感。而標準的kω模型不能準確預測壓力誘導的分離，使得激波位置靠后，而又限制了分離的發展。完全雷諾應力模型，也就是二階矩模型考慮了湍流剪切應力的輸運，理論上說能預測壓力誘導的分離流動，但計算結果表明DRSM較兩方程的標準k－ω模型只是在激波位置上較好，這讓人不由得想到在對二階矩模型壓力再分布項模化的時候使用的SSG模型。SSG壓應力模型較LRR壓應力模型簡單，它不需要LRR模型那樣，為獲得較好的對數層的解，而加壁面反射作用修正，這也許就是在強逆壓梯度區SSG模型未能較好地預測出激波位置的原因。因此，更為合理的微分雷諾應力模型應當綜合SSG模型的簡潔性和LRR模型的準確性，在邊界層內使用LRR－ω模型，在邊界層外使用SSG－ε模型，通過類似與Menter的BSL模型中的混合函數來過渡。

5 結 論

渦粘性假設和代數應力模型都是基于雷諾應力方程而簡化建立的模式假設。雷諾應力方程是對雷諾平均假設下最完整的描述湍流的控制方程，因此本文著眼于直接對雷諾應力方程的建模，也就是利用微分雷諾應力模型求解激波分離流這種復雜的流動；由于雷諾應力方程引入了22個未知量，需大量合理的模式假設將其封閉。本文總結前人工作的基礎上，應用成熟的一系列假設將雷諾應力方程封閉，并針對本文應用的流場特點，引入了壓縮性修正和改進的耗散率方程，構造出可實現的微分雷諾應力模型（DRSM）；將此模型和傳統的一方程模型（SA）和兩方程模型（SST、Wilcox標準k－ω）應用于較大攻角的跨音速ONERA M6機翼的計算，結果表明，SA和SST在流場未出現大的分離前能較好地預測激波位置但分離程度都過大，出現激波誘導的大分離后分離流態和激波位置都嚴重失真，而DRSM在激波誘導的分離區域表現較好。

［1］SPALART P R，SHUR M L.On the sensitization of turbulence models to rotation and curvature［J］.AerospaceScienceandTechnology，1997，1（5）：297－302.

［2］WILCOX D C.Turbulence modeling for CFD［R］.DCW Industries，Inc.1994.

［3］SPEZIALE C G，Sarkar S，Gatski T G.Modelling the pressure－strain correlation of turbulence［J］.Fluid Mech，1991，227：245－272.

［4］WALLIN S，JOHANSSON A.An explicit algebraic Reynolds stress model for incompressible and compressible flows［J］.FluidMech，2000，403：89－132.

［5］WILCOX D C.Reassessment of the scale determining equation for advanced turbulence models［J］.AIAA Journal，1988，26（11）：1299－1310.

［6］SCHMITT V，CHARPIN F.Pressure distributions on the ONERA－M6－wing at transonic Mach numbers［R］.AGARD－AR－138.May，1978.

［7］MENTER F R.Two－equation eddy－viscosity turbulence models for engineering applications［J］.AIAAJournal，1994，32（8）：1598－1605.