基于雷諾應力模型的高精度分離渦模擬方法?

王圣業 王光學 董義道 鄧小剛

1)(國防科學技術大學航天科學與工程學院,長沙 410073)2)(中山大學物理學院,廣州 510275)

基于雷諾應力模型的高精度分離渦模擬方法?

王圣業1)王光學2)董義道1)鄧小剛1)?

1)(國防科學技術大學航天科學與工程學院,長沙 410073)2)(中山大學物理學院,廣州 510275)

湍流流動,雷諾應力模型,分離渦模擬,加權緊致非線性格式

1 引 言

準確預測翼型和機翼在接近甚至超過失速迎角時的氣動特性,已成為現代飛行器設計的重要方面.然而,對于高雷諾數大迎角條件下的分離流動,精確的計算流體力學(computational fluid dynamics,CFD)模擬仍十分困難[1].

任何流動模擬的結果都依賴于所代表的物理模型的準度以及求解相關方程的數值方法的精度[2].為了精確預測流動物理,主要是湍流,最好是直接求解所有尺度的湍流脈動甚至是最小尺度,即直接數值模擬;或求解大尺度湍流而模型化最小尺度,即大渦數值模擬(large eddy simulation,LES).然而,在真實飛行雷諾數下,即使附著邊界層內的最大尺度脈動也會變得很小,這將使計算花費巨增而難以承擔.因此,基于雷諾平均納維-斯托克斯方程(Reynolds average Navier-Stokes,RANS)的湍流模型方法仍然是工業應用CFD的中堅力量[3].

混合雷諾平均/大渦模擬(hybrid RANS/LES)方法結合了大渦模擬方法在分離流動區域的高分辨率與雷諾平均方法在附著邊界層中的高效率,在大迎角分離模擬中受到了廣泛的關注[4?6].而在hybrid RANS/LES方法中,Spalart等[7]于1997年提出的分離渦模擬(detached eddy simulation,DES)方法由于其簡單、高效等特點得到了廣泛的發展和應用[8].然而,最初版本的DES方法基于一方程Spalart-Allmaras(SA)模型,在大規模分離流動支配的算例中表現良好,但對于由逆壓梯度造成的初始尾緣分離等問題卻很難準確估計.因此,為了更準確地模擬翼型和機翼的失速分離問題,尤其在失速迎角附近,DES的RANS部分需要采用更先進的湍流模型[9].2001年,Strelets[10]提出了基于二方程剪切壓力傳輸(shear stress transport,SST)模型的DES方法,對翼型最大升力迎角附近時的分離問題有一定的提高.然而,由于線性渦粘模型自身的限制,其效果仍不理想.2006年,Greschner等[11,12]發展了基于EASM LLk-ε模型(一類非線性渦粘模型)的DES方法,并在串聯圓柱-翼型算例中與傳統k-ε-DES方法進行了對比.結果表明,EASM-DES方法能夠更精細地分辨脫落渦結構,并且預測的氣動力更接近試驗值.但由于EASM模型仍基于Boussinesq假設框架,其對雷諾應力的求解仍然不直接,因此模擬效果對算例的敏感性較大.

雷諾應力模型(Reynolds stress model,RSM)是更高階矩封閉的湍流模型,采用六個方程求解雷諾應力項和一個方程求解尺度約束變量(例如湍動能耗散率ε或比耗散率ω),能夠直接估計雷諾應力各向異性、流線彎曲效應等[13].2011年,Probst等[9]基于ε型RSM模型,發展了εh-DES類方法,并在HGR-01翼型上取得了很好的結果.但由于ε型RSM模型魯棒性較差,其適應性將不如ω型RSM模型[14].

近年來,德國宇航局的Eisfeld等[13,15,16]發展了一種魯棒的RSM模型,Speziale-Sarkar-Gatski/Launder-Reece-Rodi(SSG/LRR)-ω模型,并在一系列復雜的航空分離流動中證明了其優勢.美國航空航天局的Rumsey等[2]在湍流模型資源(turbulence modeling resource,TMR)網站上將SSG/LRR-ω模型設為最推薦的RSM,并開展了廣泛的驗證與確認工作.國內,董義道等[17]基于TMR網站,成功開展了SSG/LRR-ω模型的初步應用研究,但在應用過程中發現基于RANS方法的SSG/LRR-ω模型對翼型大迎角狀態的模擬仍不理想.本文基于上述工作,在SSG/LRR-ω模型上發展了一類DES方法,并在典型航空算例(17?,45?及60?迎角的NACA0012翼型;12?迎角的NACA4412翼型;24.6?迎角的鈍前緣三角翼)中進行了驗證.為了進行對比研究,本文還同時采用了傳統基于線性渦粘模型的分離渦模擬方法,SST-DES方法.

對于DES方法,在RANS部分采用更先進的湍流模型,而在分離區LES部分則對數值精度的要求更高.相比傳統二階方法,高階方法在采用較少計算花費獲得更低誤差方面具有更大的潛力[18?20].由于在計算效率和計算精度上的優勢,高階方法受到了廣泛關注,并被推薦用于LES或hybrid RANS/LES方法中[21].本文中對RANS方程和模型方程的離散均采用了高階精度的加權緊致非線性格式(weighted compact nonlinear scheme,WCNS)[22,23],并且在計算網格導數時采用了滿足自由流保持的對稱網格導數算法(symmetrical conservative metric method,SCMM)[24].

2 湍流模型

2.1 SSG/LRR-ω模型

對Navier-Stokes方程進行Favre平均后出現雷諾應力項.傳統線性渦粘模型,如SA模型和SST模型,對雷諾應力項的求解均基于Boussinesq假設,而RSM則直接對其求解.

其中t,xk(k=1,2,3)分別為時間和空間坐標分量;,分別表示Favre平均的密度、速度分量和雷諾應力張量,而全雷諾應力不加特殊說明,下面出現的變量均為Favre平均后的量.(1)式右端第一項為生成項,可準確求解:

右端第三項為耗散項,

其中ε表示各向同性耗散率,需要通過額外的湍流尺度方程求解得到,δij為Kronecker符號.右端第四項為輸運項,

其中各向異性張量

表1 SSG/LRR-ω模型封閉系數中SSG和LRR部分的貢獻Table 1.Values of closure coefficients for the SSG and the LRR contributions to the SSG/LRR-ω redistribution term.

其他各項定義如下:

為了封閉上述模型,需要額外引入尺度方程,SSG/LRR-ω模型借鑒了Menter的思路,對ε方程和ω方程進行了混合[15]:

方程(9)的變量為比耗散率ω.各向同性耗散率通過下式計算:

混合函數定義如下:

其中,dw為距壁面法向距離,其他系數具體見表2.

表2 ω方程系數中ε和ω部分的貢獻Table 2.Values of coefficients of ω-equation corresponding to the ε and ω parts.

2.2 SSG/LRR-DES類方法

DES方法由Spalart等[7]提出,將LES方法和SA模型結合,并通過比較當地網格尺度與壁面距離實現兩種方法的自動切換.其后,Strelet[10]借鑒Spalart的思想,通過比較當地網格尺度與湍流長度尺度,將DES方法引入SST模型.本文采用的SSG/LRR-ω模型與SST模型均為基于比耗散率ω的湍流模型,因此在構造SSG/LRR-DES方法時將主要借鑒SST-DES的構造思路.

首先定義網格長度尺度

和湍流長度尺度

然后得到SSG/LRR-DES的限制器

對雷諾應力方程耗散項中的各向同性耗散率(10)式進行修正:

這里CDES為DES常數,通過Menter的混合函數(12)式得到

原始DES方法的限制器在復雜網格上處理RANS和LES的轉換過程中過于生硬,會造成模化應力損耗(modeled stress depletion,MSD)等問題[8].2006年,Spalart等[26]借鑒Menter的SST模型構造思想,采用“延遲LES函數”改善了原始版本的MSD問題,發展出了延遲分離渦模擬(delayed DES,DDES)方法.2008年,Shur等[27]將DDES方法和壁面模型大渦模擬(wall-modeled LES,WMLES)方法結合,發展了強化延遲分離渦模擬(improved DDES,IDDES)方法,并在非定常分離流動中得到了廣泛的應用.本文的計算中均采用IDDES形式,其構造過程如下.

首先采用WMLES中網格尺度的加權思想,定義新的IDDES限制器:

其中,fd為原DDES中的轉換函數,fd=1?tanh(8rd)3,fB為原WMLES中的轉換函數,fB=min(2exp(?9α2),1),這里α=0.25?dw/hmax;fe為壁面模擬的控制函數,其作用是保證在壁面附近網格分辨率滿足LES要求時,忽略過渡區雷諾應力的影響.

其中fe1定義為

fe2定義為

(22)式中rd,rdt和rdl均與文獻[27]中相同.

另外,在計算當地網格尺度時,Shur等[27]還引進了壁面距離的影響,新的Lg為

其中hmax=max(?x,?y,?z),hwn為壁面垂直方向的網格步長.IDDES中的系數為Cw=0.15,ct=1.87和cl=5.0.

3 數值方法

3.1 高精度數值方法

WCNS格式由Deng等[22]在2000年提出.其后,不同學者[20,28,29]發展了多種形式的WCNS格式,并在廣泛的流動問題中證明了該格式的優勢.本文采用的是WCNS系列格式中一種典型的五階顯式離散格式WCNS-E-5[23].WCNS-E-5格式由于其低耗散、高魯棒和優秀的自由流與渦保持特性,被廣泛應用于各種實際流動問題的高精度數值模擬中.其中岡敦殿等[30]將WCNS-E-5格式應用于平板圓臺突起物繞流的LES中,并和采用平面激光散射技術的試驗結果進行了對比,證明了格式精細模擬湍流流動的可行性.

另外本文在計算網格導數及雅克比時,采用了滿足幾何守恒律的對稱守恒網格導數算法[24],有利于提高高精度有限差分方法的魯棒性并減小數值誤差.本文中時間推進均采用子迭代步基于LU-SGS(lower-upper symmetric-Gauss-Seidal)方法的雙時間步法.

3.2 對雷諾應力方程的高精度離散

較差的數值穩定性是限制RSM使用的障礙,尤其是結合高精度數值方法時.在迭代的過程中,雷諾應力項很可能不滿足Schumann[31]提出的現實性限制.Chassaing等[32]提出一種魯棒的隱式格式,依賴于使用當地雙時間步技術和顯式增強現實性限制.Yossef[33]發展了一種無條件正收斂隱式格式和一種現實性限制,保證了任意時間步內雷諾正應力項為正值.本文中,在雷諾應力方程的迭代過程中添加如下限制條件:

除了時間迭代方面,也應同時關注雷諾應力方程的空間離散.對于加權型有限差分格式,非線性權是影響數值穩定性的主要方面.需要指出的是,在對雷諾應力方程采用WCNS類格式離散時,加權公式中的小量εIS應取為10?18,而非文獻[23]中對Navier-Stokes方程離散時的10?6.

4 計算結果

4.1 NACA0012翼型大迎角氣動特性模擬

NACA0012翼型大迎角分離算例是推動DES發展的經典算例.2001年,Strelets等[7]在提出SST-DES類方法時,就在該算例上與原始的SADES類方法進行了對比,但發現并無本質的差異.2016年,Yang和Zha[34]結合高精度WENO(weighted essentially nonoscillatory)格式,采用改進的SA-IDDES方法對該算例進行了模擬,但發現在中等迎角下SA-IDDES方法還沒有SA-URANS方法準確.本文為了更好地對SSG/LRR-IDDES方法進行對比研究,計算條件和網格均參考Yang和Zha的工作[34],并同時采用了SST-IDDES方法.

圖1 NACA0012翼型三維粗糙網格Fig.1.Three-dimensional coarse grid of the NACA0012 airfoil.

NACA0012翼型基于弦長c的雷諾數為1.3×106,基于自由流速度的馬赫數為0.5.在0?—90?迎角范圍內,選取17?(大范圍邊界層分離),45?和60?(大規模脫體分離)三個典型狀態進行模擬.由于該算例在高雷諾數(Re>1.0×105)下受雷諾數影響很小,因此升力和阻力系數的參考值選取Re=2.0×106的試驗值[34,35].本文計算采用粗和密兩套O型網格,網格單元數分別為192×102×30和288×102×30.第一層網格小于1;遠場長度約為100c;展向長度為1.0c.圖1為NACA0012翼型三維網格示意圖.時間推進采用雙時間步方法,周期為T=c/U∞,時間步長0.01T.時間平均統計從50T時開始,持續50T.

圖2展示了3個典型迎角下統計時間的平均升、阻力系數.在17?迎角時,翼型處于失速狀態,所有方法得到的阻力系數相近且與試驗吻合較好;但對于升力系數,SSG/LRR-IDDES在粗網格上略低于試驗值,而在密網格上與試驗值吻合很好;而其他方法均明顯高于試驗值.Yang和Zha[34]在17?迎角的計算中,得出了WENO+SA-IDDES結果明顯差于WENO+SA-URANS的現象,但并未給出解釋.而本文采用的線性渦粘模型、SST模型并未出現該現象.在45?和60?迎角時,翼型處于過失速狀態,所有IDDES方法的結果均明顯優于URANS方法.同時對比兩類URANS方法,SSG/LRR-ω模型略優于SST模型,而該現象與相關湍流模型在分離渦處對雷諾應力預測的能力有關.

圖3展示了17?迎角時升、阻力系數在100T內的變化過程,其中AoA表示迎角.可以看到在25T后,SST-URANS和SSG/LRR-URANS得到的升、阻力系數出現類似簡諧振蕩的發展過程.而SSTIDDES和SSG/LRR-IDDES方法得到的升、阻力系數均為無周期振蕩,表明其均能夠描述分離湍流的隨機性.對比試驗值,SST-IDDES得到的升力系數明顯偏高,并且隨網格加密無明顯改善,反映了基于線性渦粘模型的DES方法在翼型最大升力迎角附近模擬的局限.

經典的渦拉伸原理[36]表明,三維性是湍流最本質的特性之一.圖4給出的粗網格上100T時刻Q判據為0的三維等值面圖中,SST-URANS預測的渦脫落過程明顯為二維過程;SSG/LRRURANS能夠預測出上翼面附近較小尺度的渦,但整個渦脫落過程仍未表現出明顯的三維性;而SST-IDDES和SSG/LRR-IDDES均得到了高度混亂的三維渦結構,但SSG/LRR-IDDES對前緣渦拉伸-彎曲-破裂的過程則模擬的更為精細.

圖5對比了不同方法在粗網格上100T時刻0.5展長處的展向速度,其中SST-URANS得到的展向流動十分微弱,也印證了其為二維渦脫落過程;SSG/LRR-URANS得到的展向流動稍強于SSTURANS,但也未表現出明顯的三維性;而SSTIDDES和SSG/LRR-IDDES均預測出了明顯的展向流動.

圖2 (網刊彩色)0?—90?范圍內3個典型迎角下的升阻力結果對比Fig.2.(color online)Comparisons of lift and drag coefficients at 3 typical attack angles of among 0?–90?.

圖3 (網刊彩色)17?迎角時升、阻力系數變化過程Fig.3.(color online)Lift and drag coe ff cient history at attack angle of 17?.

圖4 (網刊彩色)17?迎角下100T時刻Q判據為0的等值面圖,顏色由馬赫數標識Fig.4.(color online)Iso-surface of the Q-criterion=0 at 100T and attack angle of 17?,colored by Mach number.

圖5 (網刊彩色)17?迎角下100T時刻,0.5展長處展向速度云圖Fig.5.(color online)Spanwise velocity contours of 50%span at 100T and attack angle of 17?.

圖6展示了45?迎角時升、阻力系數在100T內的變化過程.可以看到在25T后,SST-URANS得到的升、阻力系數出現類似簡諧振蕩的發展過程,而SSG/LRR-URANS并未得到有規律的簡諧振蕩過程,這與RSM能夠估計雷諾應力的各向異性有關.而SST-IDDES和SSG/LRR-IDDES方法得到的升、阻力系數均在試驗值附近無規律振蕩.

圖6 (網刊彩色)45?迎角時升、阻力系數變化過程Fig.6.(color online)Lift and drag coe ff cient history at attack angle of 45?.

圖7 (網刊彩色)45?迎角下100T時刻Q判據為0的等值面圖,顏色由馬赫數標識Fig.7.(color online)Iso-surface of the Q-criterion=0 at 100T and attack angle of 45?,colored by Mach number.

來流流過機翼后形成脫體的分離湍流,并將產生明顯的展向流動.圖7給出了粗網格上45?迎角,100T時刻,Q判據為0的三維等值面圖.SSTURANS僅預測出了大尺度的展向渦,其渦脫落過程仍為二維過程;SSG/LRR-URANS能夠預測出少量的流向和法向渦,但仍未表現出明顯的三維過程;而SST-IDDES和SSG/LRR-IDDES得到了高度混亂的大小渦結構,預測出了明顯的三維效應. 圖8對比了不同方法在粗網格上100T時刻0.5展長處的展向速度,更加清晰地表明各方法對于分離湍流的預測能力. 這也是圖3中SST-IDDES和SSG/LRR-IDDES的結果優于SSG/LRR-URANS、更優于SST-URANS的原因.

圖8 (網刊彩色)45?迎角下100T時刻,0.5展長處展向速度云圖Fig.8.(color online)Spanwise velocity contours of 50%span at 100T and attack angle of 45?.

圖9 (網刊彩色)60?迎角時升、阻力系數收斂過程Fig.9.(color online)Lift and drag coe ff cient history at attack angle of 60?.

圖10 (網刊彩色)60?迎角下100T時刻Q判據為0的等值面圖,顏色由馬赫數標識Fig.10.(color online)Iso-surface of the Q-criterion=0 at 100T and attack angle of 60?,colored by Mach number.

圖11 (網刊彩色)60?迎角下100T時刻,0.5展長處展向速度云圖Fig.11.(color online)Spanwise velocity contours of 50%span at 100T and attack angle of 60?.

圖9展示了60?迎角時升、阻力系數的變化過程.與45?迎角時的結果類似,在25T后,SSTURANS得到的升、阻力系數出現明顯的周期性;而SSG/LRR-URANS，SST-IDDES和SSG/LRRIDDES呈現出無規律振蕩,但SSG/LRR-URANS得到的升、阻力系數平均值明顯高于試驗值.圖10給出了60?迎角100T時刻下,Q判據為0的三維等值面圖.SST-URANS預測的渦脫落過程仍為二維過程;SSG/LRR-URANS得到的結果也仍未表現出明顯的三維過程;而SST-IDDES和SSG/LRRIDDES得到了高度混亂的大小渦結構,預測出了明顯的三維效應.圖11給出了60?迎角100T時刻0.5展長處的展向速度云圖,結果也與45?迎角時類似.

4.2 NACA4412翼型尾緣分離模擬

NACA4412翼型是經典的低速湍流驗證算例[37],有多種條件下的詳細試驗數據進行對比.本文選擇的是以Wadcock[38]實施的最大升力構型(12?迎角)試驗為參考的算例.試驗條件為:翼型弦長0.9 m,雷諾數1.64×106,自由流速度29.1 m/s以及馬赫數0.085.該條件下翼型尾緣開始出現分離,對任何hybrid RANS/LES方法都是巨大的挑戰[39].

本文計算的條件設置與試驗相同.采用粗和密兩套O型網格,展向拉伸0.1 m,網格單元數分別為169×123×30和249×123×30.第一層網格均小于1.圖12為NACA4412翼型三維網格示意圖.時間推進采用雙時間步法,周期為T=c/U∞,時間步長0.01T.時間平均統計從40T時開始,持續40T.

表3列出了平均升力系數和分離位置的計算值與試驗值.可以看到,SST-URANS、SSG/LRRURANS和SST-IDDES得到的結果均與試驗值偏差較大,并且隨著網格加密未有改善. 而SSG/LRR-IDDES的結果與試驗吻合較好.

圖12 NACA4412翼型三維粗網格Fig.12.Three-dimensional coarse grid of the NACA4412 airfoil.

表3 NACA4412翼型平均升力系數和分離位置的計算值與試驗值對比Table 3.Comparation of computional results of lift coefficient and location of separation for NACA4412 airfoil case.

圖13 (網刊彩色)NACA4412翼型80T時刻,Q判據為0等值面圖,顏色由馬赫數標識Fig.13.(color online)Iso-surface of the Q-criterion=0 at 80T for NACA4412 airfoil case,colored by Mach number.

圖13給出了粗網格上80T時刻Q判據為0的等值面圖.在翼型尾緣處,流動開始分離并產生湍流尾跡渦.SSG/LRR-IDDES成功模擬了該過程,并且分離位置與試驗吻合.圖14給出了該網格上得到的平均壓力系數分布和試驗值的對比.SSG/LRR-IDDES在前緣得到的吸力值明顯低于另三種方法,并且與試驗值更接近,這也是其得到更準確升力系數的原因.在翼型尾緣處,由于SSG/LRR-IDDES成功模擬了尾跡渦結構,其得到的壓力系數也更接近試驗值.圖15展示了翼型尾緣(x/c=0.952)和尾跡(x/c=1.282)兩個典型站位處的流向速度分布對比.在x/c=0.952處,SSG/LRR-IDDES成功捕捉到了逆向速度,計算結果與試驗吻合.表明SSG/LRR-IDDES能在邊界層附近更好地感受逆壓梯度的影響,對于弱非定常流動較傳統SST-IDDES方法有一定提高.而僅采用URANS方式會引入過多湍流黏性,使弱非定常運動的發展受到抑制.在x/c=1.282處,SSG/LRR-IDDES得到的流向速度略大于試驗值,且與其他三種方法無明顯區別.

圖14 (網刊彩色)NACA4412翼型表面壓力系數分布Fig.14.(color online)Distribution of pressure coefficient on the surface of NACA4412 airfoil.

圖15 (網刊彩色)NACA4412翼型流向速度分布比較Fig.15.(color online)Comparison of streamwise velocity pro files for NACA4412 airfoil.

4.3 鈍前緣三角翼渦破裂模擬

現代戰斗機和導彈多采用三角翼布局,以獲得良好的飛行品質和機動性能.然而在大迎角下,三角翼會產生渦破裂現象對氣動特性造成影響.本文采用的計算模型為65?后掠三角翼,是NASA Langley中心為研究該外形的雷諾數、馬赫數影響而完成的試驗模型[40].該模型分為四部分:前緣、平板部分、后緣及整流罩.前緣部分提供4種可替換外形,本文選用中等半徑鈍前緣外形.來流條件為:Ma=0.85,Re=6×106.迎角選擇24.6?,為渦破裂現象發生的臨界迎角.

本文采用的計算網格自主生成,網格單元數約6百萬,網格結構為多塊對接網格.圖16給出了三角翼的表面網格.網格拓撲采用C-H型,以避免翼尖奇性軸的產生.計算區域的遠場邊界取為50倍根弦長.壁面的第一排網格達到了10?6弦長,網格在背風區和剪切層附近均進行了適當的加密,以保證分離區、附面層內和剪切層的數值模擬精度.時間推進采用雙時間步法,周期為T=cref/U∞,其中cref為參考氣動弦長,時間步長取0.01T.時間平均統計從10T時開始,持續20T.

圖16 鈍前緣三角翼計算網格Fig.16. Computational mesh for the blunt-edge deltawing.

圖17展示了鈍前緣三角翼不同站位處壓力分布與試驗[41]的對比,其中η=2z/bl,bl為當地展長,cr為翼根弦長.回顧鈍前緣三角翼分離的基本特性.相比傳統尖前緣,鈍前緣三角翼在前緣處有轉捩過程[42],但本文計算的馬赫數較高,迎角也較大,因此轉捩過程很短,分離渦為全湍流狀態.主渦旋轉產生展向流動,在旋渦下部、三角翼的上翼面加速,出現負壓力峰值,該點沿展向到翼邊緣是逆壓梯度,將誘導邊界層分離,產生二次分離和二次旋渦,在主渦負壓力峰值外側又出現二次負壓力峰值.觀察x/cr=0.4站位的壓力分布,SST-IDDES,SSG/LRR-URANS和SSG/LRR-IDDES均成功捕捉到了二次渦結構,而其中SSG/LRR-IDDES與試驗值吻合最好.

當迎角大到一定程度時,主渦開始破裂.本文選擇的24.6?迎角,為渦破裂現象發生的臨界迎角,而x/cr=0.6是該迎角下渦破裂發生前的臨界站位.因此該站位的模擬對CFD方法是個挑戰.觀察本文計算值與試驗值的對比,四種方法均產生了不同程度的偏差.其中SSG/LRR-IDDES在吸力峰處與試驗吻合很好,但在翼根處吻合較差.x/cr=0.6站位的翼根處臨近整流罩頭部(存在速度駐點),壓力系數由負值迅速變為接近1從而形成較大的壓力梯度,因此給數值模擬帶來很大困難.

圖17 (網刊彩色)鈍前緣三角翼不同站位處壓力分布對比Fig.17.(color online)Comparisons of surface pressure with experiments for blunt-edge deltawing at typical stations.

圖18 (網刊彩色)三角翼20T時刻Q判據等值面圖,顏色由壓力系數標識Fig.18.(color online)Iso-surface of the Q-criterion at 20T for deltawing case,colored by pressure coefficient.

圖19 (網刊彩色)三角翼20T時刻流線圖,顏色由馬赫數標識Fig.19.(color online)Streamlines at 20T for deltawing case,colored by Mach number.

在x/cr=0.8處,渦破裂已發生,上翼面吸力峰消失.SST-IDDES,SSG/LRR-URANS和SSG/LRR-IDDES均成功地預測了該現象,而SST-URANS仍得到了主渦及二次渦結構.結合其他站位上SST-URANS得到的壓力分布,可以看出在24.6?迎角下,其并未預測到渦破裂現象,即推遲了渦破裂的發生.而該現象的產生明顯由SSTURANS在分離區過大的湍流黏性導致.由于來流為跨聲速,前緣渦引起的上翼面流動加速,形成了局部超音速區域并產生了激波.對于跨聲速流動,Rumsey[14]在ONERO M6機翼和NASA CRM翼身組合構型中均表明,SSG/LRR-ω模型能更好地預測激波附近雷諾應力的劇烈變化,對激波誘導分離的模擬能力明顯優于SST等渦粘模型.本文發展的SSG/LRR-IDDES方法繼承了RSM的優勢,在非定常區(x/cr=0.8和x/cr=0.95)得到了優于SST-IDDES方法的結果.

圖18和圖19分別展示了20T時刻三角翼Q判據等值面圖和流線圖.SST-URANS得到的仍為定常狀態的前緣分離渦,在24.6?迎角時未發生渦破裂.SSG/LRR-URANS雖然預測到了渦破裂,但渦破裂后的流動未表現出明顯的非定常特性.表明對渦破裂后形成的非定常運動,需采用諸如DES,LES等的尺度模擬方法.SST-IDDES和SSG/LRR-IDDES均較好地模擬了渦破裂后的非定常流動狀態,但就精細度而言,SSG/LRRIDDES要更為優秀.

5 結 論

本文借鑒SST-IDDES方法的構造方式,在SSG/LRR-ωRSM上發展了SSG/LRR-IDDES分離渦模擬方法.通過結合高精度WCNS格式在NACA系列翼型及鈍前緣三角翼算例上進行了驗證,并和傳統SST-IDDES方法以及SSG/LRRURANS和SST-URANS方法進行了對比.

為說明SSG/LRR-IDDES方法在翼型大迎角氣動力模擬方面的提升,選取NACA0012翼型17?,45?和60?三個典型狀態進行模擬.在17?迎角時,翼型邊界層分離擴大到前緣附近,升力系數接近最小.此時,采用SST-URANS和SSG/LRR-URANS均會過高估計升力,而采用傳統基于線性渦粘模型的SST-IDDES方法也并未明顯改善.這說明對于傳統基于線性渦粘模型的IDDES方法,由于邊界層附近RANS部分本身對逆壓梯度引起的分離的模擬能力較差,使得整個方法較URANS并無明顯提高.而基于RSM,可以更準確地估計逆壓梯度區的雷諾應力變化,使得結合IDDES時得到了很好的結果.到45?和60?迎角時,翼型產生大范圍脫體渦,升力系數回升.此狀態是傳統IDDES方法模擬的優勢方面,而SSG/LRR-IDDES也繼承了該特性,并較SSG/LRR-URANS有較大提高.

為說明SSG/LRR-IDDES方法在翼型最大升力迎角附近的模擬能力,選取NACA4412翼型12?狀態進行模擬.此時,翼型尾跡由于逆壓梯度影響開始分離并產生湍流尾跡渦.采用SST-URANS和SSG/LRR-URANS均會過高估計升力,并延遲了尾跡分離;而采用SST-IDDES也未明顯改善.但是采用SSG/LRR-IDDES方法給出了與試驗吻合很好的氣動力特性和準確的分離過程,表明該方法對模擬翼型失速迎角附近特性的能力提高.

為進一步說明SSG/LRR-IDDES方法在三維機翼失速迎角附近的模擬能力,選取鈍前緣三角翼24.6?狀態進行模擬.該迎角下,三角翼主渦由于激波誘導的影響,在x/cr=0.6站位后發生非定常渦破裂現象.采用SST-URANS方法無法成功預測渦破裂現象;而采用SSG/LRR-URANS方法雖然預測到了渦破裂,但渦破裂后的流動未表現出明顯的非定常特性.采用SST-IDDES方法和SSG/LRR-IDDES方法均較好地模擬了渦破裂后的非定常流動狀態,但在表面壓力分布、流場精細度等方面,SSG/LRR-IDDES方法更加優秀.

另外,本文采用的WCNS-E-5格式能夠在粗糙網格上獲得較高的保真度,體現了高精度數值方法在模擬分離流動中的效率優勢.但在加密網格計算時,也出現了網格收斂性差的現象,如NACA0012翼型60?迎角時,SSG/LRR-IDDES在粗網格下得到的平均氣動力系數略優于密網格.該現象的產生與IDDES方法基于網格長度尺度來調整湍流黏性有關,也是DES類方法需要克服的問題之一[8].下一步,將繼續結合高精度格式對SSG/LRR-IDDES方法在更廣泛的流動中進行驗證并提高網格收斂性.

[1]Slotnick J,Khodadoust A,Alonso J,Darmofal D,Gropp W,Lurie E,Mavriplis D 2014CFD Vision 2030 Study:A Path to Revolutionary Computational Aerosciences(Washington,DC:Langley Research Center,NASA)Tech.Rep.NASA/CR-2014-218178

[2]Eisfeld B,Rumsey C,Togiti V 2016AIAA J.54 1524

[3]Rumsey C 201452nd Aerospace Sciences MeetingNational Harbor,Maryland,January 13–17,2014 AIAA 2014-0201

[4]Tucker P 2006Int.J.Numer.Meth.Fluids51 261

[5]Richez F,Pape A,Costes M 2015AIAA J.53 3157

[6]Xu G,Jiang X,Liu G 2016Acta Mech.Sin.32 588

[7]Spalart P,Jou W H,Strelets M,Allmaras S 1997Comments on the Feasibility of LES for Wings,and on Hybrid RANS/LES Approach(Columbus:Greyden Press)

[8]Spalart P 2009Annu.Rev.Fluid Mech.41 181

[9]Probst A,Radespiel R,Knopp T 201120st AIAA Computational Fluid Dynamics ConferenceHonolulu,Hawaii,June 27–30,2011 AIAA 2011-3206

[10]Strelets M 200139th AIAA Aerospace Sciences Meeting and ExhibitReno,NV,8-11 January 2001,AIAA 2001-0879

[11]Greschner B,Thiele F,Gurr A,Casalino D,Jacob M200612th AIAA/CEAS Aeroacoustics ConferenceCambridge,Massachusetts,May 8–10,2006 AIAA 2006-2628

[12]Greschner B,Thiele F,Jacob M,Casalino D 2008Comput.Fluids37 402

[13]Cécora R D,Radespiel R,Eisfeld B,Probst A 2015AIAA J.53 739

[14]Rumsey C 2015 in Eisfeld B(ed.)Differential ReynoldsStress Modeling for Separating Flows in Industrial Aerodynamics(Springer Tracts Mechanical Engineering)p19

[15]Eisfeld B,Brodersen O 200523rd AIAA Applied Aerodynamics ConferenceToronto,Ontario Canada,June 6–9,2005 AIAA 2005-4727

[16]Togiti V,Eisfeld B,Brodersen O 2014J.Aircraft51 1331

[17]Dong Y D,Wang D F,Wang G X,Deng X G 2016J.National Univ.Defense Technol.38 46(in Chinese)[董義道,王東方,王光學,鄧小剛2016國防科技大學學報38 46]

[18]Shu C 2003Int.J.Comput.Fluid D17 107

[19]Wang Z,Fidkowski K,Abgrall R,Bassi F,Caraeni D,Cary A,Deconinck H,Hartmann R,Hillewaert K,Huynh H,Kroll N,May G,Persson P O,van Leer B,Visbal M 2013Int.J.Numer.Meth.Fluids72 811

[20]Wang S,Deng X,Wang G,Xu D,Wang D 2016Int.J.Comput.Fluid D30 469

[21]Georgiadis N,Rizzetta D,Fureby C 2010AIAA J.48 1772

[22]Deng X,Zhang H 2000J.Comput.Phys.165 22

[23]Deng X,Liu X,Mao M,Zhang H 200517th AIAA Computational Fluid Dynamics ConferenceToronto,Ontario Canada,June 6–9,2005 AIAA 2005-5246

[24]Deng X,Mao M,Tu G,Liu H,Zhang H 2011J.Comput.Phys.230 1100

[25]Bellot G,Corrsin S 1971J.Fluid Mech.48 273

[26]Spalart P,Deck S,Shur M,Squires K,Strelets M,Travin A 2006Theor.Comput.Fluid Dyn.20 181

[27]Shur M,Spalart P,Strelets M,Travin A 2008Int.J.Heat Fluid Fl.29 1638

[28]Nonomura T,Fujii K 2009J.Comput.Phys.228 3533

[29]Liu H,Ma Y,Yan Z,Mao M,Deng X 20148th International Conference on Computational Fluid DynamicsChengdu,China,July 14–18,2014 ICCFD8-2014-0082

[30]Gang D D,Yi S H,Zhao Y F 2015Acta Phys.Sin.64 054705(in Chinese)[岡敦殿,易仕和,趙云飛2015物理學報64 054705]

[31]Schumann U 1977Phys.Fluids20 721

[32]Chassaing J,Gerolymos G,Vallet I 2003AIAA J.41 763

[33]Yossef Y 2014J.Comput.Phys.276 635

[34]Yang Y,Zha G 201646th AIAA Fluid Dynamics ConferenceWashington,D.C.,USA,June 13–17,2016 AIAA 2016-3185

[35]Shur M,Spalart P,Strelets M,Travin A 1999Proceedings of the 4th International Symposium on Engineering Turbulence Modelling and MeasurementsCorsica,France,May 24–26,1999 p669

[36]Chen M Z 2002Fundamentals of Viscous Fliud Dynamics(Beijing:Higher Education Press)p239(in Chinese)[陳懋章 2002 粘性流體動力學基礎 (北京:高等教育出版社)第239頁]

[37]Chen Y,Guo L D,Peng Q,Chen Z Q,Liu W H 2015Acta Phys.Sin.64 134701(in Chinese)[陳勇,郭隆德,彭強,陳志強,劉衛紅2015物理學報64 134701]

[38]Wadcock A 1987Investigation of Low Speed Turbulent Scparatcd Flow Around Airfoils(Washington,DC:Ames Research Center,NASA)Tech.Rep.NASA-CR-177450

[39]Roy R,Stoellinger M 201553rd AIAA Aerospace Sciences MeetingKissimmee,Florida,January 5–9,2015 AIAA 2015-1982

[40]Luckring M,Hummel D 2013Aerosp.Sci.Technol.24 77

[41]Chu J,Luckring M 1996Experimental Surface Pressure Data Obtained on 65 deg Delta Wing Across Reynolds Number and Mach Number Ranges.Vol.3:Medium-Radius Leading Edge(Washington,DC:Ames Research Center,NASA)NASA-TM-4645-Vol-3

[42]Luckring M 2013Aerosp.Sci.Technol.24 10

High-order detached-eddy simulation method based on a Reynolds-stress background model?

Wang Sheng-Ye1)Wang Guang-Xue2)Dong Yi-Dao1)Deng Xiao-Gang1)?

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

20 March 2017;revised manuscript

14 May 2017)

Referring to the construction of shear stress transport-improved delayed detached-eddy simulation(SST-IDDES)method,a variant of IDDES method based on the Speziale-Sarkar-Gatski/Launder–Reece–Rodi(SSG/LRR)-ωReynoldsstress model(RSM)as Reynolds-averaged Navier-Stokes(RANS)background model,is proposed.Through combining high-order weighted compact nonlinear scheme(WCNS),the SSG/LRR-IDDES method is applied to three aeronautic cases and compared with traditional methods:SST-unsteady Reynolds-averaged Navier-Stokes(URANS),SSG/LRRURANS,and SST-IDDES.To verify the SSG/LRR-IDDES method in simulating airfoil stalled flow,NACA0012 airfoil is adopted separately at attack angles of 17?,45?and 60?.At the attack angle of 17?,SST-URANS,SSG/LRR-URANS,and SST-IDDES methods each predict a higher lift coefficient than the experimental data,while the SSG/LRR-IDDES method obtains a better lift coefficient result and a higher fidelity vortical flow structure.It indicates that the RSM can improve the prediction of RANS-mode for pressure-induced separations on airfoil surfaces in detached-eddy simulation.At the attack angles of 45?and 60?,the SSG/LRR-IDDES method captures the massively separated flow with threedimensional vortical structures and obtains a good result,which is the same as that from the traditional SST-IDDES method.To indicate the improvement of the SSG/LRR-IDDES method in simulating airfoil trailing edge separation,NACA4412 airfoil is adopted.At the attack angle of 12?(maximum lift),the trailing edge separation is mainly induced by pressure gradient.The SSG/LRR-IDDES method can predict the separation process reasonably and obtains a good lift coefficient and location of separation compared with experimental results.However,none of other methods can predict trailing edge separation.It confirms that when RSM is adopted as RANS background model in detached-eddy simulation,the ability to predict pressure-induced separation on airfoil surface is improved.For further verifying the SSG/LRR-IDDES method for simulating three-dimensional separated flow,blunt-edge deltawing at the attack angle of 24.6?is adopted.At this attack angle,the primary vortex will break,which is difficult to predict by using the SSTURANS method.For the SSG/LRR-URANS method,it predicts the vortex breakdown successfully,but the breakdown process does not show any significant unsteady characteristic.The SST-IDDES and the SSG/LRR-IDDES methods both predict a significant unsteady vortex breakdown.But in terms of the accuracy of surface pressure and the fidelity of unsteady flow,the result obtained by the SSG/LRR-IDDES method is better than by the SST-IDDES method.

turbulence flows,Reynolds stress model,detached eddy simulation,weighted compact nonlinear scheme

(2017年3月20日收到;2017年5月14日收到修改稿)

基于Speziale-Sarkar-Gatski/Launder-Reece-Rodi(SSG/LRR)-ω雷諾應力模型發展了一類分離渦模擬方法,結合高精度加權緊致非線性格式在典型翼型及三角翼算例中進行了驗證,并和傳統基于線性渦粘模型的分離渦模擬方法進行了對比.結果表明:基于SSG/LRR-ω模型的分離渦模擬方法,提高了原雷諾應力模型對非定常分離湍流的模擬能力;同時相比于傳統基于線性渦粘模型的分離渦模擬方法,尤其是在翼型最大升力迎角和三角翼渦破裂迎角附近,該方法在平均氣動力預測的準確度、分離湍流模擬的精細度等方面更加優秀.

10.7498/aps.66.184701

?國防科學技術大學科研計劃(批準號:ZDYYJCYJ20140101)資助的課題.

?通信作者.E-mail:xgdeng2000@vip.sina.com

感謝中山大學國家超級計算廣州中心在計算資源方面提供的幫助.

PACS:47.27.–i,47.27.em,47.27.ep,47.11.BcDOI:10.7498/aps.66.184701

*Project supported by the Foundation of the National University of Defense Technology of China(Grant No.ZDYYJCYJ20140101).

?Corresponding author.E-mail:xgdeng2000@vip.sina.com