ＳＳＴ湍流模型在強逆壓流模擬中的應(yīng)用

摘 要：本文對SST兩方程湍流模型提出了一些改進，并將之與Navier-Stokes方程的顯式有限元格式相結(jié)合，對高雷諾數(shù)可壓流進行數(shù)值模擬。流動方程和湍流方程以一種松散的關(guān)聯(lián)方式求解。時間方向上，流動方程用顯式多步Runge-Kutta格式推進，而湍流方程用多步點隱式格式推進以解除源項引起的時間步長限制。數(shù)值實踐中我們考核了跨音速激波、大攻角等引起的強逆壓流動的經(jīng)典算例，并與實驗結(jié)果進行了比較，結(jié)果是令人滿意的。

關(guān)鍵詞：湍流模型；有限元；N-S方程；強逆壓流動

中圖分類號：O355 文獻標(biāo)志碼：A

Applications of SST Turbulence Model in Prediction of Flows

with Strong Adverse Pressure Gradient

TANG Shimin

(Changzhou Radio and TV University， Changzhou Jinagsu 213001，China)

Abstract：In this paper， we introduce some improvements in the well-known SST turbulence model and combine it with an explicit finite element scheme for Navier-Stokes equations. The flow equations and turbulence equations are solved in loosely coupled manner. The flow equations are advanced using a multi-stage Runge-Kutta time stepping scheme， while the turbulence equations are advanced in a multi-stage point-implicit scheme to alleviate the time-step restriction due to the source term. We test the approach for strong adverse pressure gradient flows induced by transonic shock and large angle of attack. The comparison with experiment data shows that our numerical results are satisfactory.

Key words：turbulence model; explicit finite element; Navier-Stokes equations; adverse pressure flows

空氣動力學(xué)中Navier-Stokes方程的一個主要應(yīng)用就是模擬湍流。對于大規(guī)模的分離流、非定常流以及包含多長度尺度的流動，定義代數(shù)長度尺度是困難的，B-L、J-K等代數(shù)湍流模型的應(yīng)用變得非常復(fù)雜和不確定。顯然，數(shù)值方法的改進必須跟更一般的湍流模型的發(fā)展和應(yīng)用同步進行。

常用的兩方程湍流模型是k-ω和模型k-ε。k-ε模型不能準(zhǔn)確地推測逆壓流動，在粘性底層中也存在數(shù)值困難。Wilcox［1］提出的k-ω模型克服了k-ε模型的許多弱點，特別是在近壁面區(qū)域中，其性能比k-ε模型有較大的提高；然而，k-ω模型對自由來流的ω值非常敏感。Menter［2］通過混合函數(shù)將這兩種模型結(jié)合起來，既保留了k-ω模型在近壁面區(qū)域、k-ε模型在自由剪切層中各自的優(yōu)良特性，又克服了k-ω模型對來流ω值的敏感性；在此基礎(chǔ)上，通過對渦粘系數(shù)定義的限制性修正，考慮了逆壓邊界層中湍流剪切應(yīng)力的傳輸效應(yīng)，進一步設(shè)計出SST（Shear Stress Transport）模型。與k-ω和k-ε模型相比，SST模型的性能有了較大提高，尤其對強逆壓流中分離開始和分離程度的推測。

強逆壓邊界層流動是湍流模擬的難點之一。本文對SST模型提出了一些改進，并將之與N-S方程的顯式有限元解法相結(jié)合。湍流方程用多步點隱式格式進行時間推進以解除源項引起的時間步長限制。數(shù)值實驗中，我們考核了跨音速激波、大攻角等引起的強逆壓流動的經(jīng)典算例。

1 SST湍流模型及改進

1.1 SST模型

Menter首先把k-ε模型轉(zhuǎn)化為k-ω模型的形式。轉(zhuǎn)化形式和原始k-ω模型的差別在于ω方程右端多了一交叉擴散項，且方程的系數(shù)不一樣。轉(zhuǎn)化過程中還略去了一個對結(jié)果無影響小擴散項。然后應(yīng)用開關(guān)函數(shù)F1把兩個模型結(jié)合起來：

DρkDt=xj（μ+σkμt）kxj+τuixj-β*ρωk（1-1）

DρωDt=xj（μ+σωμt）ωxj+γvtτijuixj-βρω2+2ρ(1-F1)σω21ωkxjωxj

（1-2）

其中湍流剪切應(yīng)力

τij=μtuixj+ujxi-23ukxkδij-23ρkδij，渦粘系數(shù)vt=kω。

F1在緊靠壁面處為1，激活模型k-ω；離開壁面向0趨近，激活轉(zhuǎn)化的模型k-ε。這種混合模型的性能在邊界層內(nèi)50%左右的區(qū)域接近模型，而流場其它部分接近于k-ω模型。用1、2分別表示k-ω模型和轉(zhuǎn)化k-ε模型中的常數(shù)，則混合模型中的常數(shù)為：

=F11+(1-F1)2（1-3）

對原k-ω始模型

σk1=0.5，σω1=0.5，β1=0.0750，β*=0.09，κ=0.41，γ1=β1/β*-σω1κ2/β*

對轉(zhuǎn)化k-ε模型

σk2=1.0，σω2=0.856，β2=0.0828，β*=0.09，κ=0.41，γ2=β2/β*-σω2κ2/β*

混合函數(shù)F1的定義為

F1=tanh(arg41)(1-4)

arg1=minmaxk0.09ωy，500vy2ω，4ρσω2kCDkωy2(1-5)

CDkω=max2ρσw21ωkxjωxj，10-20

y表示當(dāng)前點到物面的最短距離。

渦粘模型和全雷諾應(yīng)力模型的一個主要差別在于后者包含了傳輸項，從而考慮了重要的湍流應(yīng)力傳輸效應(yīng)。為考慮這種效應(yīng)，Menter在湍流動能的生成超過耗散的地方，對渦粘值進行了限制，從而構(gòu)建了SST模型。重新定義的渦粘系數(shù)公式為

vt=a1kmax(a1ω，|Ω|F2)（1-7）

|Ω|為渦張量的標(biāo)量度量，常數(shù)a1=0.31。F2在邊界層中等于1，在自由剪切層中等于0。自由剪切層中Bradshaw假設(shè)未必成立，應(yīng)恢復(fù)使用原始渦粘公式。在逆壓邊界層中，生成項大于耗散項，即Ω＞a1ω。開關(guān)函數(shù)F2的定義為

F2=tanharg22(1-8)

arg2=max2k0.09ωy，500vy2ω（1-9）

另外，SST模型中重新定義σk1=0.85。

1.2 邊界條件

無窮遠處，自由來流值取為

ω∞=(1~10)U∞L，vt∞=10-(2-5)v∞，k∞=vt∞ω∞

其中L為計算域長度。

在無滑移壁面上，除ω外所有湍流量均取零。壁面上的ω值按如下簡單方式確定

ω=1060vβ1(Δy)2（1-10）

其中Δy為壁面到最近一點的距離；只要第一點處y+小于2.5，該方式是合適的［3］。

入流邊界上湍流量取確定值，出流邊界上則假定是零梯度的。

1.3 改進

（1-6）中的CDkω代表交叉擴散項的正值部分。數(shù)值實踐表明以（1-6）式表達CDkω會使（1-5）最后一項的數(shù)值特性較差，導(dǎo)致F1在邊界層的外邊緣以外升高到1。應(yīng)用下式取代（1-6）可獲得較好的效果

CDkω=max(2ρσω2)1ωkxj

ωxj，10-8CDmax)（1-11）

對繞流問題，由于后緣駐點，為使湍流方程時間推進收斂，ω應(yīng)修正為［4］

ω=maxω，(3/(8|S|2)-1/2(1-12)

|S|為變形率張量的標(biāo)量度量。該修正對最后的收斂結(jié)果沒有影響。

（6）中的|Ω|可用|S|代替

|S|=2SijSij，Sij=12uixj+uixj（1-13）

Boussinesq近似是基于變形率而非渦量，這個替代應(yīng)是合理的；對簡單邊界層計算，上式與（1-7）幾無差別。

由于Bradshaw假設(shè)在非常靠近壁面的粗糙層中不成立，設(shè)計一函數(shù)F3防止SST限制在該區(qū)域激活［5］，該函數(shù)可設(shè)計成

F3=1-tahn150vωy24（1-14）

相應(yīng)地，（1-7）修正為

vt=a1kmax(a1ω，|Ω|F2F3)（1-15）

2.數(shù)值方法

2.1 控制方程

流動控制方程為Favre平均的N-S方程，雷諾應(yīng)力和熱通量項通過Boussinesq假設(shè)模擬。守恒形式下的流動方程可寫成

wt+fcx+gcay=fx+gyy（2-1）

其中為流動解的守恒變量，fc、gc為對流通量，fv、gv為粘性通量。

湍流方程（1-1）、（1-2）也可相似地寫成

wt+fcx+gcy=fvx+gyy+h（2-2）

這里w=［ρk，ρω］T，而fc、gc為相應(yīng)的對流項，fv、gv為擴散項，h為源項。

求解過程包括對流動方程和湍流方程在三角形網(wǎng)格上的空間離散，然后將離散方程在時間方向上推進至定常狀態(tài)。

2.2 空間離散

流動方程應(yīng)用Galerkin有限元離散，流動變量記于三角形節(jié)點上。對流通量在節(jié)點上計算，并假定在三角形單元上是線性變化的。對于粘性通量，由于流動變量本身在網(wǎng)格單元上線性變化，因此粘性應(yīng)力中的速度梯度在單元形心上計算。此外還需附加人工粘性項以保證對流項的穩(wěn)定，該項由拉普拉斯算子和雙調(diào)和算子混合構(gòu)建而成。邊界處的人工粘性項需進行了適當(dāng)修正［6］，同時高雷諾數(shù)繞流計算中還考慮了網(wǎng)格的壓縮性影響［7］；這樣使人工粘性項與物理粘性項相比盡可能小，又保證邊界層外捕捉激波的需要。除激波附近外，空間近似在全流場是二階精度的。

對湍流方程，擴散項的離散相似地采用Galerkin有限元離散，變量在三角形單元上線性變化。源項中的速度梯度也通過線性單元假設(shè)計算。然而，湍流方程中的對流通量項是通過一階逆風(fēng)格式近似的。盡管僅有一階精度，這種方式有助于保證時間推進過程中湍流變量的穩(wěn)定性和正值性。與N-S方程不同，湍流方程中的對流項并非主導(dǎo)項；所以湍流解對該項的近似精度并不敏感。

2.3 時間推進

離散的流動方程以顯式五步R-K格式進行時間推進，并應(yīng)用當(dāng)?shù)貢r間步長、隱式殘值光順加速迭代收斂。原則上，湍流方程也可以同樣的格式進行時間推進。然而，源項的存在會強加更嚴格的時間步長限制。如果流動方程和湍流方程以完全關(guān)聯(lián)的顯式進行時間推進，必須取流動方程和湍流方程的最小時間步長，在源項占優(yōu)的區(qū)域會引起很慢的收斂。另一方面，如果以不關(guān)聯(lián)的顯式方式推進，湍流方程可能會大大滯后于流動方程，從而不能推進到定常狀態(tài)。為了使湍流方程和流動方程以同樣的速度推進，必須對源項進行隱式處理。這里并非簡單地對源項進行隱式處理，湍流方程以一種點隱式的方式進行時間方向上的推進。首先，離散的湍流方程可以寫成如下形式

ΔwiΔt=r(wi)+h(wi)(2-3）

這里r(wi)代表離散的對流項和擴散項，取決于點及其周圍節(jié)點上w的值；h(wi)代表離散的源項，取決于點上的值。然后，將上述方程進行關(guān)于點的線性化處理，從而得到關(guān)于Δwi的方程

Δwi=1Δt-rw-hw-1(r(wi)+h(wi))(2-4)

顯式R-K格式在湍流方程的時間推進中以多步點隱式格式代替，第q步可寫成

w(q)=w(0)+1αqΔt-rw-hw-1（2-5）

(r(q-1)(wi)+hq-1(wi))

其中αq為R-K系數(shù)，Δt為當(dāng)?shù)貢r間步長。如此隱式處理使得由于源項引起的時間步長限制被解除。

在源項的線化處理中，下列近似給出了很好的數(shù)值性質(zhì)：

(ρk)(Pk-Dk)≈2Dkρk

（2-6）

(ρω)(Pω+Cω-Dω)≈|Cω|+2Dwρω(2-8)

這里P、D、C分別是產(chǎn)生項、擴散項、附加的交叉擴散項。上述表達式改變符號后將移到格式的左端，從而增加對角占優(yōu)性。

兩方程湍流模型的數(shù)值實踐表明，在ω值較低的區(qū)域，剪切變形率小的擾動會導(dǎo)致遠場和邊界層邊緣附近渦粘系數(shù)巨變。為理解這種效應(yīng)，從不可壓流的k-ω模型導(dǎo)出下列渦粘傳輸方程：

DvtDt=(1-γ)vtωuiaxj+ujxiuixj-(β*-β)k+……（2-8）

如果ω趨近于零且vt為有限值，變形率小擾動區(qū)域的渦粘產(chǎn)生項將趨近無窮。簡單處理該問題的方式是將湍流動能產(chǎn)生項限制為

Pk=min(Pk，20Dk)(2-9)

該限制器對收斂結(jié)果不產(chǎn)生影響，只是在變形率出現(xiàn)波動的區(qū)域消除了渦粘系數(shù)的振蕩。

3 數(shù)值實驗

3.1 RAE2822翼型例10

繞翼型的跨音速流動，為評價湍流模型估算邊界層沿曲面固壁在包括由激波引起的逆壓梯度下發(fā)展的能力，提供了很好的檢驗。由于流動在激波的下游出現(xiàn)了分離，該例是對湍流模型一個嚴厲的考核。來流條件為：M∞=0.75、Rec=6.2×106、α=3.19°，計算中考慮到洞壁干擾，攻角修正為2.8°。轉(zhuǎn)捩點由湍流模型自動確定。

圖1給出的表面壓強系數(shù)和摩擦力系數(shù)分布與實驗結(jié)果［8］ 的比較表明：SST很好地計算了激波位置，而模型沒有計算出流動分離，J-K模型計算的激波位置不夠準(zhǔn)確。所有模型都沒有計算出前緣上表面的壓力峰值，這也許是實驗中用于轉(zhuǎn)捩的糙粒造成的。

圖2給出的三個位置的速度型表明：由于考慮了湍流動能的傳輸，SST模型比k-ω模型更好地反應(yīng)了逆壓作用，從而較準(zhǔn)確地計算了流動的分離現(xiàn)象。坐標(biāo)原點位于弦線中點。

3.2 NACA4412翼型繞流

該算例的來流參數(shù)為：α=13.87°、Rec=1.5×106，來流馬赫數(shù)M∞根據(jù)實驗［9］條件推算為0.078。該流動已非常接近于失速狀態(tài)。

圖3給出了計算得到翼型后緣流場狀況：弦向速度分量等值線和流線圖。SST模型計算得到的渦尺度和位置與實驗結(jié)果接近。

圖4給出的多個站位的速度型與實驗結(jié)果也很接近。

圖1 RAE2822翼型表面壓強系數(shù)（左）、摩擦力系數(shù)分布（右）

（點：實驗，實線：SST，虛線：k-ω，點劃線：J-K）

Fig.1 Surface pressure and shear stress distribution for RAE2822 airfoil

dot:experiment， solid line:SST， dashed line:， dot line:J-K

圖2 RAE2822翼型上表面速度型（點：實驗，實線：SST，虛線： k-ω）

Fig.2 Velocity profile on the upper surface of RAE2822 airfoil

dot:experiment， solid line:SST， dashed line:k-ω

圖3 NACA4412翼型局部弦向速度分量等值線和流線圖

Fig.3 Partial display of chord-wise mean velocity contour and stream lines for NACA4412 airfoil

圖4 NACA4412翼型上表面速度型，（點：實驗，實線：SST）

Fig.4 Velocity profile on the upper surface of NACA4412 airfoil (dot:experiment， solid line:SST)

4.結(jié)論

分析以及數(shù)值結(jié)果與實驗結(jié)果的比較表明，與代數(shù)模型或一般兩方程模型相比，SST湍流模型能較好地模擬強逆壓流動。本文所采取的流動方程和湍流方程松散關(guān)聯(lián)的求解方式也具有較好的收斂性。

參考文獻

［1］ Wilcox D C. Reassessment of the Scale-Determining Equation for Advanced Turbulence Models ［J］. AIAA J， Vol. 26， No. 11， 1988:pp1299-1310.

［2］Menter F R. Zonal Two EquationTurbulence Models for Aerodynamic Flows ［C］. AIAA Paper 1993(6):2009

using the Wilcox， EASM and SST Models ［C］. 41st Aerospace Sciences Meeting and Exhibit Reno， Nevada， 6-9 Jan 2003.

［3］Goodheart K， Dykas S and Schnerr G H. Numerical Insights into the Solution of Transonic Flow Test Cases

［4］Michelassi V. Shock-boundary layer interaction and transition modeling in turbo-machinery flows ［J］. J. of power and energy， Vol. 211， 1997:pp 225-234.

［5］Hellsten A. Some Improvement in Menter’sSST Turbulence Model ［C］. AIAA Paper 98-2554，June 1998.

［6］周春華.非結(jié)構(gòu)網(wǎng)格上Euler方程有限體積解法的改進 ［J］.南京航空航天大學(xué)學(xué)報，Vol. 23， No. 3， 2003:pp313-317.

［7］D.J. Mavriplis and A. Jameson. Multigrid solution of the Navier-Stokes equations on triangular meshes［J］. AIAA J， Vol.28， No. 8， 1990:pp1415-1425.

［8］P.H. Cook， M.A. McDonald and M.C.P. Firmin. Airfoil RAE2822 pressure distributions and boundary layer and wake measurement， AGARD Advisory Report No. 138， May 1979.

［9］Donald Coles and Alan J Wadcock， Flying-Hot-Wire Study of Flow Past an NACA 4412 Airfoil at Maximum Lift ［J］. AIAA J， Vol. 17， No. 4， 1979:pp321-329.