微分雷諾應(yīng)力模型在激波分離流中的應(yīng)用

聶勝陽(yáng)，高正紅，黃江濤

（西北工業(yè)大學(xué) 翼型葉柵空氣動(dòng)力國(guó)防科技重點(diǎn)實(shí)驗(yàn)室，陜西 西安710072）

0 引 言

基于Boussinesq渦粘性假設(shè)的湍流模型在很多工程流動(dòng)類(lèi)型上有很好的預(yù)測(cè)能力，但是對(duì)一些特殊的流動(dòng)，它會(huì)給出與實(shí)驗(yàn)較大的差別。平均應(yīng)變率的突變導(dǎo)致這些湍流模型失效并不奇怪，雷諾應(yīng)力的相應(yīng)變化是與自由來(lái)流進(jìn)程和時(shí)間尺度不相關(guān)的變化率衡量，所以導(dǎo)致Boussinesq渦粘性假設(shè)失效。同樣，當(dāng)流動(dòng)中因快速膨脹，劇烈變形，或者大的流線(xiàn)彎曲產(chǎn)生的額外應(yīng)變率，都會(huì)產(chǎn)生不平衡的法向雷諾應(yīng)力，渦粘性假設(shè)就失效了。

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

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

1 不可壓雷諾應(yīng)力方程

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

其中：

稱(chēng)為擴(kuò)散項(xiàng)；由式（1）可以看出，第一，方程中已經(jīng)包含了τij的對(duì)流和擴(kuò)散，所以二階矩模型包含流動(dòng)發(fā)展歷程的影響。耗散項(xiàng)和湍流輸運(yùn)項(xiàng)中的時(shí)間尺度是和平均流動(dòng)時(shí)間尺度不相關(guān)的，所以二階矩模型的上游歷程的影響要比兩方程湍流模型中人為添加的上游歷程的影響更真實(shí)；第二，式（1）包含了對(duì)流項(xiàng)、產(chǎn)生項(xiàng)和體積力項(xiàng)，這些項(xiàng)對(duì)流線(xiàn)彎曲，流體旋轉(zhuǎn)等起主要作用。所以二階矩模型如果?；暮?，也能更真實(shí)地包含這些產(chǎn)生不平衡雷諾應(yīng)力的作用；第三，式（1）并沒(méi)有反應(yīng)出法向應(yīng)力相等，即使平均應(yīng)變率為0，也就是說(shuō)，它的值由初始條件和流動(dòng)進(jìn)程決定。因此，二階矩模型能夠處理平均應(yīng)變率有突變的情況。

2 模式假設(shè)

為了封閉方程（1），必須對(duì)耗散項(xiàng)張量，再分布項(xiàng)張量，擴(kuò)散項(xiàng)張量進(jìn)行?；?。耗散項(xiàng)：

因?yàn)楹纳⒆饔冒l(fā)生在最小尺度的渦，因此耗散項(xiàng)的建模就相對(duì)簡(jiǎn)單，大家都采用Kolmogorov（1941）提出的當(dāng)?shù)馗飨蛲约僭O(shè)，即：

這里

標(biāo)量ε就是湍動(dòng)能方程中的耗散率。而對(duì)耗散率的建?？蓞⒖糼－ε兩方程模型中的耗散率方程的處理方法。

由于實(shí)際中的耗散率是各向異性的，尤其是靠近壁面處，所以為引入這一特性，后人進(jìn)一步修正耗散項(xiàng)為：

其中bij為無(wú)量綱化的雷諾應(yīng)力各向異性張量，即

fs是經(jīng)驗(yàn)性的添加的衰減函數(shù)，為

ReT為定義的湍流雷諾數(shù)，為

擴(kuò)散項(xiàng)：

擴(kuò)散項(xiàng)中包含壓力脈動(dòng)項(xiàng)和三階速度脈動(dòng)項(xiàng)，表示脈動(dòng)壓力和雷諾應(yīng)力對(duì)能量的輸運(yùn)作用?，F(xiàn)有的實(shí)驗(yàn)數(shù)據(jù)無(wú)法找到有用的信息去模型化壓力脈動(dòng)關(guān)聯(lián)項(xiàng)。最簡(jiǎn)單而又常見(jiàn)的?；瘮U(kuò)散項(xiàng)的方法是假定一個(gè)梯度輸運(yùn)進(jìn)程，并且引入分子粘性作用。這里使用Launder，Reece和 Rodi（1975）［2］提出來(lái)的一種近似方法：

這里他們建議Cs≈0.11。這也是渦粘性假設(shè)的湍流模型處理擴(kuò)散項(xiàng)常用的方法。再分布項(xiàng)：

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

這里bij定義如式（5）。

Speziale，Sarkar和Gatski給出的模式參數(shù)為：

Willin ＆Johansson［4］給出了一組較簡(jiǎn)單的模式參數(shù)為：

3 可實(shí)現(xiàn)的DRSM

3.1 引入壓縮性修正

可壓流動(dòng)的特征是密度的顯著變化。在馬赫數(shù)小于5（另一說(shuō)為8）時(shí)［2］，壁面流動(dòng)受壓縮性影響很小，因?yàn)橥ㄟ^(guò)激波，在小距離內(nèi)壓力梯度變化不劇烈。Tw／Te在邊界層內(nèi)壓縮性對(duì)渦的作用也不明顯?；谶@樣的事實(shí)，Morkovin（1962）認(rèn)為湍流的密度脈動(dòng)量相對(duì)與自由來(lái)流是個(gè)微量，因此在湍流模式中使用不可壓的形式［2］。但 Morkovin假設(shè)在即使不是高超音速情形下也有其局限性，因?yàn)槊芏让}動(dòng)量相對(duì)自由來(lái)流的密度并不是小量。所以壓縮性效應(yīng)明顯的流動(dòng)情形，應(yīng)將密度脈動(dòng)量考慮進(jìn)去，一種有效的數(shù)學(xué)處理方法就是采用Favre（1965）的質(zhì)量平均速度的概念，定義為：

某時(shí)刻瞬時(shí)速度為：

這樣在湍流模型中也就反映了密度的脈動(dòng)。這里直接給出Favre平均的雷諾應(yīng)力方程，具體推導(dǎo)過(guò)程見(jiàn)參考文獻(xiàn)［2］：

其中：

產(chǎn)生項(xiàng)：Pij＝－

再分布項(xiàng)：Πij＝

耗散項(xiàng)：εij＝

擴(kuò)散項(xiàng)：Dij＝－

壓力作用項(xiàng)：PWw＝

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

由于壓力作用項(xiàng)是脈動(dòng)量和平均量的耦合，還沒(méi)有合適的模型去表征它的特點(diǎn)，并且單個(gè)脈動(dòng)量在時(shí)間平均中是可以忽略掉的，本文去掉壓力作用項(xiàng)。

3.2 耗散項(xiàng)εij的改進(jìn)

前面介紹二階矩的模型化過(guò)程，但沒(méi)有給出耗散率ε的方程。這是因?yàn)?，在壁面上，耗散率為ε?，所以在壁面上，式（1.7）、（1.8）中均出現(xiàn)奇點(diǎn)。在計(jì)算中，ε會(huì)在壁面附近劇烈變化，因此在壁面附近需將網(wǎng)格劃分的足夠的小。因此本文采用Wilcox（1988）等［5］提出的比耗散率ω的輸運(yùn)方程來(lái)?；纳㈨?xiàng)。輸運(yùn)方程為：

模式參數(shù)為：

耗散率補(bǔ)充方程為：

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

4 模型驗(yàn)證

選用ONERA M6機(jī)翼，計(jì)算狀態(tài)為Ma＝0.84，α＝6.06，Re＝11.7×106。這是一個(gè)跨音速大迎角狀態(tài)，實(shí)驗(yàn)數(shù)據(jù)來(lái)自參考文獻(xiàn)［6］。

本文計(jì)算采用的離散格式為中心格式添加適當(dāng)?shù)娜斯ふ承裕瑢?duì)時(shí)間推進(jìn)采用顯式三步Runge－Kutta推進(jìn)。采用多重網(wǎng)格和隱式殘差光順等措施加速收斂。計(jì)算網(wǎng)格為非結(jié)構(gòu)網(wǎng)格，粘性層用25層三棱柱，第一層的高度為1e－6，外推比例為1.3。為了對(duì)不同湍流模型計(jì)算效果進(jìn)行比較，計(jì)算過(guò)程采用同一套分布合理的網(wǎng)格，使用相同的時(shí)間推進(jìn)手段和CFL數(shù)，自由來(lái)流的湍流度和湍流粘性比相同。

圖1和圖2給出了表面極限流線(xiàn)計(jì)算結(jié)果，可以看出M6機(jī)翼在該狀態(tài)時(shí)，出現(xiàn) 形激波，在翼展60%處，激波匯合，形成更強(qiáng)的激波，波后由于強(qiáng)的逆壓梯度，出現(xiàn)大范圍的分離渦。比較發(fā)現(xiàn)，SA模型和SST模型的分離區(qū)域最大，在機(jī)翼中段激波后就出現(xiàn)了輕度分離。匯合后的激波靠近前緣，尤其是SA模型。SST模型預(yù)測(cè)出來(lái)的激波不平直。而基于Wilcox提出的標(biāo)準(zhǔn)k－ω兩方程，以及基于標(biāo)準(zhǔn)k－ω兩方程的SSG－ω模型和Willin ＆Johansson提出的改進(jìn)SSG－ω模型預(yù)測(cè)出相似的極限流線(xiàn)分布和激波形狀和位置，W＆J改進(jìn)的SSG－ω模型預(yù)測(cè)出的分離區(qū)域，流線(xiàn)更紊亂，標(biāo)準(zhǔn)k－ω模型預(yù)測(cè)出的激波位置更靠近下游，激波后的流動(dòng)更光順。

從壓力分布與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比上，80%展向位置之前，SST模型給出最接近實(shí)驗(yàn)的激波位置，但波后的分離較大，整個(gè)后半?yún)^(qū)域都是低壓區(qū)；在80%展向外，SA和SST模型對(duì)大的分離流動(dòng)和激波位置預(yù)測(cè)失真，而基于SSG壓力關(guān)聯(lián)項(xiàng)假設(shè)的DRSM預(yù)測(cè)的激波位置最接近實(shí)驗(yàn)，波后的分離流動(dòng)雖然不是很貼合，但也可以接受。而 Wilcox的標(biāo)準(zhǔn)k－ω模型給出的激波位置靠后。

圖1 計(jì)算網(wǎng)格、表面壓力系數(shù)分布和表面極限流線(xiàn)圖（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 等壓線(xiàn)分布圖Fig.2 Iso－surface pressure coefficient lines distributions

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

SA和SST都過(guò)度地預(yù)測(cè)了激波后的分離流動(dòng)，這可能是在這種出現(xiàn)強(qiáng)的激波附面層干擾和強(qiáng)逆壓梯度時(shí)，兩種模型都沒(méi)有很好地限制空間湍流粘性系數(shù)的非物理性增加，而基于k－ω模型的計(jì)算在激波位置的預(yù)測(cè)上顯得滯后，這與附面層內(nèi)湍流粘性系數(shù)的過(guò)度增長(zhǎng)有關(guān)，但在激波后，又很好地限制了這種增長(zhǎng)。而DRSM模型直接對(duì)雷諾應(yīng)力分量建模，在各向同性流體中就退化成了Wilcox標(biāo)準(zhǔn)k－ω模型。

可以總結(jié)出，DRSM模型較一方程和兩方程模型，如本文第一部分所述，對(duì)快速膨脹，流線(xiàn)彎曲等流體各向異性性質(zhì)顯著的流動(dòng)預(yù)測(cè)更加合理。Willin＆Johanson的改進(jìn)參數(shù)在激波位置的預(yù)測(cè)上稍好，但改進(jìn)余地有限。

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

5 結(jié) 論

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

［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.