董義道，王東方，王光學(xué),2，鄧小剛

(1. 國防科技大學(xué) 航天科學(xué)與工程學(xué)院， 湖南 長沙　410073；2. 中山大學(xué) 物理學(xué)院， 廣東 廣州　510006)

?

雷諾應(yīng)力模型的初步應(yīng)用*

董義道1，王東方1，王光學(xué)1,2，鄧小剛1

(1. 國防科技大學(xué) 航天科學(xué)與工程學(xué)院， 湖南 長沙410073；2. 中山大學(xué) 物理學(xué)院， 廣東 廣州510006)

針對SSG/LRR-ω雷諾應(yīng)力模型，選取NASA湍流資源網(wǎng)站上的四個典型算例，即湍流平板邊界層流動、帶凸起管道流動、翼型尾跡區(qū)流動和NACA0012不同攻角繞流，開展初步的驗證與確認(rèn)工作，將部分結(jié)果和CFL3D進行對比。對于NACA0012翼型繞流，對比雷諾應(yīng)力模型和SA模型的升力系數(shù)，結(jié)果表明：在失速攻角附近，雷諾應(yīng)力模型明顯優(yōu)于SA模型。在此基礎(chǔ)上，將該模型應(yīng)用于DLR-F6翼身組合體的數(shù)值模擬，計算得到的機翼表面典型站位壓力分布和實驗值吻合良好，同時該模型捕捉到翼身交匯位置的小范圍分離。

雷諾應(yīng)力模型；驗證與確認(rèn)；復(fù)雜外形應(yīng)用

大渦模擬和直接數(shù)值模擬能夠精細捕捉流動特征，但是由于其計算量大，開展的相關(guān)研究還局限于學(xué)術(shù)領(lǐng)域。目前廣泛使用的求解器大多基于雷諾平均NS(Reynolds Averaged Navier Stokes, RANS)方程，即將湍流的貢獻抽象為雷諾應(yīng)力張量。為了求解雷諾應(yīng)力張量，需要引入湍流模型。絕大多數(shù)湍流模型基于線性Boussinesq假設(shè)，應(yīng)用比較廣泛的線性模型包括一方程SA模型[1]和兩方程SST模型[2]。

線性渦黏模型在很多簡單的流動問題中能夠得到可靠的模擬結(jié)果，但是無法捕捉復(fù)雜的流場特征，包括大范圍流動分離，流線彎曲效應(yīng)以及流場的各向異性。為此，發(fā)展了能夠更好反映流動物理特征的雷諾應(yīng)力模型。雷諾應(yīng)力模型又稱為二階矩封閉模型，不同于傳統(tǒng)的基于Boussinesq假設(shè)的線性模型，該模型直接求解雷諾應(yīng)力輸運方程，理論上能夠改善對分離流動、曲面流動以及各向異性效應(yīng)明顯的流動問題的模擬精度。但是，基本的雷諾應(yīng)力輸運方程無法封閉。Chou[3]和Rotta[4]最早研究雷諾應(yīng)力輸運方程的封閉性問題，之后，許多學(xué)者開展了方程封閉性研究，并形成了比較系統(tǒng)的設(shè)計準(zhǔn)則。Donaldson等[5]提出了“模型不變性”(invariant modeling)的概念，即封閉性近似應(yīng)當(dāng)嚴(yán)格滿足坐標(biāo)不變性。Lumley[6]嘗試發(fā)展了一套系統(tǒng)的步驟進行封閉性近似以保證流動的可實現(xiàn)性，即所有物理上為正的湍流相關(guān)變量在計算過程中保證正值。Speziale等[7]進一步討論了雷諾應(yīng)力模型的可實現(xiàn)性問題，并提出了簡化的設(shè)計準(zhǔn)則?；谏鲜鰷?zhǔn)則，不同學(xué)者對雷諾應(yīng)力輸運方程中的未知項，如湍流耗散項、壓力應(yīng)變關(guān)聯(lián)項的建模進行了研究，其中壓力應(yīng)變關(guān)聯(lián)項的建模是雷諾應(yīng)力模型研究的核心和關(guān)鍵。原因之一在于壓力應(yīng)變關(guān)聯(lián)項的量級和生成項相當(dāng)，其在很多工程應(yīng)用方面的流動問題中十分重要；此外，由于壓力應(yīng)變關(guān)聯(lián)項很難通過實驗測量得到，如何建立一個合理的模型需要極大的創(chuàng)造力。目前應(yīng)用比較廣泛的壓力應(yīng)變關(guān)聯(lián)項模型包括LRR模型和SSG[8]模型。此外，在對雷諾應(yīng)力輸運方程未知項建模過程中，不可避免地會涉及湍流尺度，因此為了最終實現(xiàn)模型封閉，還需要額外引入湍流尺度方程，如湍動能耗散率ε方程或比耗散率ω方程。基于不同的關(guān)聯(lián)項模型和尺度方程，發(fā)展了一系列雷諾應(yīng)力模型。包括Launder等最早發(fā)展的基于ε方程的LRR模型，在此基礎(chǔ)上，Wilcox使用比耗散率ω方程代替ε方程，發(fā)展了Stress-ω模型[9]，Eisfeld等基于Menter的思想，將LRR模型和SSG模型進行了混合，同時對尺度方程也進行了混合，發(fā)展了SSG/LRR-ω模型[10]。

1　控制方程與數(shù)值方法

1.1控制方程

控制方程包括雷諾平均NS方程[9]和雷諾應(yīng)力模型方程。

(1)

耗散項一般建模為：

其中，ε表示湍動能耗散率，需要通過額外的湍流尺度方程求解得到。輸運項包括分子黏性輸運和湍流輸運，黏性輸運張量可以直接計算得到。

(2)

式(2)的變量為比耗散率ω，各向同性耗散率通過計算得到。

式(2)中的其他系數(shù)φ=αω,βω,σω,σd通過混合函數(shù)計算得到。

φ=F1φ(ω)+(1-F1)φ(ò)

(3)

混合函數(shù)定義為：

F1=tanh(ζ4)

(4)

其中：

對于壓力應(yīng)變關(guān)聯(lián)項中的系數(shù)φ′，采用同樣的混合策略，即

φ′=F1φ(LRR)+(1-F1)φ(SSG)

(5)

1.2數(shù)值方法

采用二階MUSCL格式，平均運動方程對流通量采用Roe通量差分分裂。黏性項離散采用二階中心格式。時間推進采用LU-SGS隱式時間推進，具體公式參考文獻[11]。

關(guān)于邊界條件處理，壁面邊界采用無滑移條件，即速度分量為0；雷諾應(yīng)力分量在壁面位置為0，比耗散率ω在壁面位置滿足：

(6)

其中，β1=0.075，Δd1表示第一層網(wǎng)格距壁面距離。遠場邊界的處理基于Riemann不變量，具體方法參考文獻[11]。

2　算例驗證與確認(rèn)

為了驗證SSG/LRR-ω模型程序?qū)崿F(xiàn)的正確性，選擇NASA湍流資源網(wǎng)站[12]上的四個標(biāo)準(zhǔn)算例，包括零壓力梯度(Zero Pressure Gradient, ZPG)湍流平板邊界層流動，帶凸起管道流動，翼型尾跡區(qū)流動和NACA0012翼型繞流。數(shù)值計算均采用全湍流模擬，計算結(jié)果和同等網(wǎng)格上相應(yīng)的CFL3D結(jié)果進行對比。

2.1零壓力梯度湍流平板邊界層流動

2.1.1算例說明

該算例主要研究零壓力梯度平板湍流邊界層的發(fā)展，來流馬赫數(shù)M=0.2，來流溫度Tref=300 K，參考長度L=1 m，基于參考長度的雷諾數(shù)ReL=5×106。不加特殊說明，下標(biāo)ref表示來流參考值。

平板長2 m，前緣位于x=0位置，入口邊界設(shè)置在平板前緣上游x=-1/3 m位置。對于該算例的條件，最大邊界層厚度大約為0.03 m，計算域縱向高度取H=1 m是合理的。沿平板流向分布545個網(wǎng)格點，法向分布385個網(wǎng)格點，圖1為網(wǎng)格及邊界條件示意圖。

圖1　湍流平板網(wǎng)格及邊界條件示意圖Fig.1　Grid and boundary conditions for ZPG flat-plate

對于入口邊界(Inflow，x=-1/3 m)，總壓Pt=1.028 28Pref，總溫Tt=1.008Tref；出口邊界(Outflow，x=2 m)，背壓P=Pref，計算域上部邊界采用遠場黎曼條件(Riemann farfield)，平板壁面(x≥0,y=0)設(shè)置為絕熱壁(Adiabatic wall)，平板前緣上游(x<0,y=0)采用對稱面條件(Sym)。

2.1.2相關(guān)變量定義

相關(guān)變量的定義及其計算公式為：

(7)

其中：下標(biāo)∞表示來流，w表示壁面；ρ,u,θ,μ,τ,Reθ,Cf分別表示密度、流向速度、動量厚度、黏性系數(shù)、切應(yīng)力、基于動量厚度的雷諾數(shù)和摩阻系數(shù)；uτ,u+,y+分別表示摩擦速度、基于摩擦速度的無量綱速度和旋渦典型雷諾數(shù)。

測試主要關(guān)注以下兩條曲線：

1)壁面摩阻系數(shù)Cf隨雷諾數(shù)Reθ的變化；

2)u+隨y+的變化(Reθ=10 000)。

對于曲線2，基于Coles平均速度型的壁面律公式(Coles理論[13])為：

(8)

其中，κ為馮卡門常數(shù)，Π為型面參數(shù)。

2.1.3數(shù)值計算結(jié)果對比與分析

1)壁面摩阻對比。圖2給出了當(dāng)前程序計算的壁面摩阻系數(shù)Cf和CFL3D計算結(jié)果的對比，二者吻合得很好。

圖2　平板壁面摩阻對比Fig.2　Comparison of skin-friction for ZPG flat-plate

2)u+隨y+的變化曲線。圖3給出了線性-對數(shù)律層壁面律曲線，橫坐標(biāo)表示對y+取對數(shù)。從圖3中可以看出，當(dāng)前程序計算得到的分布和CFL3D吻合得很好，同時和Coles理論擬合的曲線基本一致。

圖3　平板對數(shù)律對比Fig.3　Comparison of logarithmic law for ZPG flat-plate

3)典型站位速度型對比。圖4給出了x=0.970 08和x=1.903 34兩個站位的速度型對比，其中橫坐標(biāo)u/Uinf表示無量綱流向速度。當(dāng)前程序計算結(jié)果和CFL3D吻合得很好。

圖4　平板典型站位速度型對比Fig.4　Comparison of velocity distribution in typical stations for ZPG flat-plate

2.2帶凸起管道流動

2.2.1算例說明

該算例下壁面為曲線邊界，因此流動存在壓力梯度。壁面曲線函數(shù)表達式為：

(9)

顯然，凸起位置區(qū)間0.3≤x≤1.2，壁面起始和終止位置分別為xwall_st=0和xwall_ed=1.5 。沿著壁面起始和終止位置分別向上游和下游延伸至xin=-25 m 和xout=25 m位置，延伸部分設(shè)置為對稱面邊界(Sym)，計算域高度取H=5 m，上邊界同樣設(shè)置為對稱面，入口及出口邊界條件設(shè)置和湍流平板邊界層算例相同。沿平板流向分布1409個網(wǎng)格點，法向分布641個網(wǎng)格點。圖5為網(wǎng)格分布及邊界條件示意圖。

圖5　凸起管道流動網(wǎng)格及邊界條件示意圖Fig.5　Grid and boundary conditions for bump flow

2.2.2數(shù)值計算結(jié)果對比與分析

圖6給出了凸起位置壓力系數(shù)Cp分布和摩阻系數(shù)Cf分布，從圖6中可以看出，當(dāng)前程序計算得到的壓力系數(shù)分布和CFL3D幾乎完全重合，計算得到的摩阻系數(shù)在峰值位置略小，但總體分布吻合得很好。此外，由于壁面前緣位置對應(yīng)對稱面邊界向絕熱壁邊界的過渡，因此解在前緣位置存在奇性，導(dǎo)致摩阻系數(shù)一定程度的振蕩，這種振蕩在CFL3D以及當(dāng)前程序計算結(jié)果中均可以觀察到。而在壁面后緣位置，存在同樣的奇性。

圖6　壁面壓力系數(shù)和摩阻分布對比Fig.6　Comparison of pressure and skin-friction for bump flow

2.3二維翼型尾跡區(qū)流動

2.3.1算例說明

該算例主要考察湍流模型對于翼型尾跡區(qū)流動的模擬能力，翼型選擇10%厚度的非對稱常規(guī)翼型(Model-A翼型)。Nakayama的實驗[14]對該翼型的流場進行了測量，需要指出的是，實驗中對翼型上下表面的轉(zhuǎn)捩位置進行了控制，而數(shù)值計算采用全湍流模擬。

來流馬赫數(shù)M=0.088 ，來流溫度T=300 K，翼型弦長c=1 m，基于弦長的雷諾數(shù)Rec=1 200 000。數(shù)值計算時遠場取20倍弦長，網(wǎng)格C型拓?fù)洌芟蚍植?121個網(wǎng)格點，法向分布193個網(wǎng)格點。圖7為壁面附近網(wǎng)格分布示意圖，壁面設(shè)置為絕熱壁，遠場采用黎曼遠場條件。主要關(guān)注尾跡區(qū)典型站位(x/c=1.01,1.05,1.20,1.40,1.80,2.19)的流向速度分布。

圖7　翼型尾跡區(qū)流動近壁面網(wǎng)格示意圖Fig.7　Near-wall grid of airfoil near-wake flow

圖8　翼型尾跡區(qū)速度型分布Fig.8　Velocity distribution of airfoil near-wake flow

2.3.2數(shù)值計算結(jié)果對比與分析

圖8給出了尾跡區(qū)典型站位的速度型分布，其中實線表示當(dāng)前程序計算結(jié)果，虛線表示CFL3D計算結(jié)果，實驗測量值用符號表示，橫坐標(biāo)u/Uinf表示無量綱流向速度。從圖8中可以明顯看出，當(dāng)前程序計算結(jié)果和CFL3D結(jié)果基本吻合，且在峰值區(qū)和實驗值吻合得更好。

2.4NACA0012翼型繞流

2.4.1算例說明

該算例主要考察NACA0012翼型在不同攻角下的表面壓力及升阻力特性，來流馬赫數(shù)M=0.15，參考長度c=1 m，基于參考長度的雷諾數(shù)Re=6×106，來流溫度300 K，來流攻角分別取0,10,15三個狀態(tài)。計算域取500倍弦長以消除遠場的影響，翼型流向分布897個網(wǎng)格點，法向分布257個網(wǎng)格點，翼型表面采用絕熱壁條件，遠場采用黎曼遠場條件。翼型表面壓力分布實驗值來源文獻[15]，不同攻角下的升阻力來源文獻[16]。圖9為近壁面網(wǎng)格分布示意圖。

圖9　NACA0012近壁面網(wǎng)格示意圖Fig.9　Near-wall grid of NACA0012 airfoil

2.4.2數(shù)值計算結(jié)果對比與分析

1)不同攻角下的表面壓力系數(shù)分布。從圖10中可以看出，當(dāng)前程序計算得到的壓力系數(shù)分布曲線和CFL3D的幾乎完全重合，并且和實驗值吻合得很好。

圖10　NACA0012翼型表面壓力系數(shù)分布Fig.10　Surface pressure coefficient of NACA0012

2)升力曲線和極曲線對比。圖11給出了升力曲線(CL-a)和極曲線(CL-Cd)，包括失速攻角附近的升力系數(shù)。當(dāng)前程序計算得到的升阻力特性和實驗值以及CFL3D吻合得很好。此外，圖11還給出了SA一方程模型的計算結(jié)果，從局部放大圖可以看出，在失速攻角附近，雷諾應(yīng)力模型計算結(jié)果明顯優(yōu)于SA模型。

圖11　升力曲線和極曲線對比Fig.11　Comparison of lift and polar curve

3　雷諾應(yīng)力模型在DLR-F6翼身組合體中的應(yīng)用

3.1計算構(gòu)型及網(wǎng)格說明

圖12　DLR-F6翼身組合體對稱面及表面網(wǎng)格示意圖Fig.12　Surface and symmetry grid of DLR-F6

本算例幾何外形為AIAA第二次阻力預(yù)測工作組(DPW Ⅱ)選擇的DLR-F6翼身組合體。該翼身組合體展長b=1.171 3 m，平均氣動弦長0.141 2 m，全模參考面積為0.145 4 m2，坐標(biāo)選擇X軸沿機身軸線向后，Y軸位于飛機縱向?qū)ΨQ平面內(nèi)，Z軸垂直于對稱面，滿足右手法則。

計算網(wǎng)格由ICEM軟件生成的多塊對接結(jié)構(gòu)網(wǎng)格，網(wǎng)格塊數(shù)為148，圖12為對稱面、機身和機翼表面網(wǎng)格示意圖。

3.2結(jié)果對比與分析

圖13給出了翼身交匯位置的流線，從圖13中可以看出，采用的雷諾應(yīng)力模型能夠捕捉到該位置的小范圍分離。

圖13　DLR-F6翼身交匯處流線Fig.13　Streamline in the wing-body intersection

圖14為機翼典型站位y/b分別為0.239,0.331,0.377,0.411的壓力分布Cp和實驗值Exp的對比曲線。從圖14中可以看出，數(shù)值計算結(jié)果和實驗值總體吻合得較好。

圖14　DLR-F6翼身組合體機翼典型站位壓力分布Fig.14　Pressure distribution in typical stations of DLR-F6

4　結(jié)論

1)主要針對SSG/LRR-ω雷諾應(yīng)力模型，選擇四個典型算例開展了初步驗證與確認(rèn)工作。為了驗證模型實現(xiàn)的正確性，將本文計算結(jié)果和同等網(wǎng)格上的CFL3D計算結(jié)果進行了對比，數(shù)值模擬結(jié)果表明了該模型實現(xiàn)正確。對于部分算例，數(shù)值計算結(jié)果和實驗進行了對比，證明了該模型能夠很好地捕捉流場特征。

2)對DLR-F6翼身組合體開展了數(shù)值模擬，典型站位壓力分布和實驗值吻合得很好，同時該模型能夠捕捉到翼身交匯位置的小范圍分離。

3)對于NACA0012翼型繞流，對比了雷諾應(yīng)力模型和SA模型的升力系數(shù)。在失速攻角附近，雷諾應(yīng)力模型計算的升力系數(shù)和實驗更加吻合。從而為后續(xù)將雷諾應(yīng)力模型應(yīng)用于大攻角大范圍分離流動，例如大攻角三角翼流動奠定基礎(chǔ)。

References)

[1]Spalart P R, Allmaras S R. A one-equation turbulence model for aerodynamic flows [C]//Proceedings of 30th Aerospace Sciences Meeting and Exhibit, 1994.

[2]Menter F R. Two-equation eddy-viscosity turbulence models for engineering applications [J]. AIAA Journal, 1994, 32(8): 1598-1605.

[3]Chou P Y. On the velocity correlations and the solutions of the equations of turbulent fluctuation [J]. Quarterly of Applied Mathematics, 1945, 3: 38.

[4]Rotta J C. Statistische theorie nichthomogener turbulenz[J]. Zeitschrift für Physik, 1951, 29: 547-572.

[5]Donaldson C, Rosenbaum H. Calculation of the turbulent shear flows through closure of the Reynolds equations by invariant modeling[R]. ARAP Report 127,Aeronautical Research Associates of Princeton, Princeton, NJ, 1968.

[6]Lumley J L. Computational modeling of turbulent flows[J]. Advances in Applied Mechanics, 1978, 18: 123-176.

[7]Speziale C G, Abid R, Durbin P A. On the realizability of Reynolds stress turbulence closures[J]. Journal of Scientific Computing, 1994, 9: 369-403.

[8]Speziale C G, Sarkar S, Gatski T B. Modelling the pressure-strain correlation of turbulence: an invariant dynamical systems approach [J]. Journal of Fluid Mechanics, 1991, 227: 245-272.

[9]Wilcox D C. Turbulence modeling for CFD [M]. La Canada,CA: DCW Industries,Inc, 1998.

[10]Cecora R D, Eisfeld B, Probst A, et al. Differential Reynolds stress modeling for aeronautics [J].AIAA Journal, 2015, 53(3): 739-755.

[11]張毅鋒. 高階精度格式(WCNS)加速收斂和復(fù)雜流動數(shù)值模擬的應(yīng)用研究 [D].綿陽:中國空氣動力研究與發(fā)展中心, 2008.

ZHANG Yifeng. Investigations of convergence acceleration and complex flow numerical simulation for high-order accurate scheme(WCNS)[D]. Mianyang:China Aerodynamics Research and Development Center, 2008.(in Chinese)

[12]Rumsey C L. NASA langley research center turbulence modeling resource[DB/OL]. [2015-08-28]. http://turbmodels.larc.nasa.gov/index.html.

[13]Coles D. The law of the wake in the turbulent boundary layer[J]. Journal of Fluid Mechanics, 1956, 1: 191-226.

[14]Nakayama A. Characteristics of the flow around conventional and supercritical airfoils[J]. Journal of Fluid Mechanics, 1985, 160: 155-179.

[15]Gregory N, O′Reily C L. Low-speed aerodynamic characteristics of NACA0012 aerofoil sections, including the effects of upper-surface roughness simulation hoar frost[R].NASA R&M 3716, 1970.

[16]Ladson C L. Effects of independent variation of Mach and Reynolds numbers on the low-speed aerodynamic characteristics of the NACA0012 airfoil section[R].NASA TM 4074, 1988.

Preliminary application of Reynolds stress model

DONG Yidao1, WANG Dongfang1, WANG Guangxue1,2, DENG Xiaogang1

(1. College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China；2. School of Physics, Sun Yat-sen University, Guangzhou 510006, China)

For the verification and validation of SSG/LRR-ωReynolds stress model, four typical two dimensional cases from NASA turbulence resources website were chosen, including zero pressure gradient flat plate, bump-in-channel, airfoil near-wake and flow over NACA0012 airfoil. A part of numerical results were in good agreement with that of CFL3D. For flow over NACA0012 airfoil, lift coefficients of Reynolds stress model and SA model were compared. It is obvious that near the stall angle of attack, Reynolds stress model has advantages over SA model. Based on these results, SSG/LRR-ωReynolds stress model was applied to the simulation of complex DLR-F6 wing-body configuration. Pressure coefficient in typical stations is comparable to that of experiment. Besides, small range of separation in the wing-body intersection is well captured.

Reynolds stress model; verification and validation; application in complex configuration

10.11887/j.cn.201604008http://journal.nudt.edu.cn

2015-09-18

國防科學(xué)技術(shù)大學(xué)科研計劃資助項目(ZDYYJCYJ20140101)

董義道(1991—)，男，江蘇淮安人，博士研究生，E-mail：tianyatingxiao@163.com；鄧小剛(通信作者)，男，教授，博士，博士生導(dǎo)師，E-mail：xgdeng2000@vip.sina.com

TN95

A

1001-2486(2016)04-046-08