ＳＳＴ湍流模型在強逆壓流模擬中的應用

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

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

中圖分類號：O355 文獻標志碼：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

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

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

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

1 SST湍流模型及改進

1.1 SST模型

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

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）

其中湍流剪切應力

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

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

=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/β*

對轉化k-ε模型

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

混合函數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表示當前點到物面的最短距離。

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

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

|Ω|為渦張量的標量度量，常數a1=0.31。F2在邊界層中等于1，在自由剪切層中等于0。自由剪切層中Bradshaw假設未必成立，應恢復使用原始渦粘公式。在逆壓邊界層中，生成項大于耗散項，即Ω＞a1ω。開關函數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ω代表交叉擴散項的正值部分。數值實踐表明以（1-6）式表達CDkω會使（1-5）最后一項的數值特性較差，導致F1在邊界層的外邊緣以外升高到1。應用下式取代（1-6）可獲得較好的效果

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

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

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

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

|S|為變形率張量的標量度量。該修正對最后的收斂結果沒有影響。

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

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

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

由于Bradshaw假設在非?？拷诿娴拇植趯又胁怀闪ⅲO計一函數F3防止SST限制在該區(qū)域激活［5］，該函數可設計成

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

相應地，（1-7）修正為

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

2.數值方法

2.1 控制方程

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

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為相應的對流項，fv、gv為擴散項，h為源項。

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

2.2 空間離散

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

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

2.3 時間推進

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

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

這里r(wi)代表離散的對流項和擴散項，取決于點及其周圍節(jié)點上w的值；h(wi)代表離散的源項，取決于點上的值。然后，將上述方程進行關于點的線性化處理，從而得到關于Δ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系數，Δt為當地時間步長。如此隱式處理使得由于源項引起的時間步長限制被解除。

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

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

（2-6）

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

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

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

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

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

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

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

3 數值實驗

3.1 RAE2822翼型例10

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

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

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

3.2 NACA4412翼型繞流

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

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

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

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

（點：實驗，實線：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.結論

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

參考文獻

［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］周春華.非結構網格上Euler方程有限體積解法的改進 ［J］.南京航空航天大學學報，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.