基于非結(jié)構(gòu)網(wǎng)格求解器的DES和DDES數(shù)值模擬研究

韓 濤，唐海龍

(1.海軍駐沈陽地區(qū)航空軍事代表室，沈陽 110034;2.中國航空工業(yè)空氣動力研究院 計算流體力學(xué)中心，沈陽 110034)

當(dāng)前工業(yè)界和研究中所使用的主要CFD手段仍然是求解RANS方程并附加湍流模型，這種根據(jù)統(tǒng)計學(xué)原理來重點研究流動平均量的研究方法的精準(zhǔn)度主要取決于湍流模型能否合理地描述被模擬的流動特征。目前湍流模型對于包含物面邊界的附著流或小分離流動大多效果較好、魯棒性強，而且模型大多通過固壁的物面率(Law of the wall)來對其中的經(jīng)驗性參數(shù)進(jìn)行校準(zhǔn),但對于自由剪切流動的模擬效果相對較差，至于對中等/大分離或強漩渦的非定常流動來說，RANS基本不能準(zhǔn)確模擬。

大渦模擬(LES)被認(rèn)為是唯一有希望解決大分離流動的數(shù)值模擬的技術(shù)，然而，盡管LES方法起源較早，但它在網(wǎng)格解析度和時間分辨率等方面的要求所帶來了巨大的計算量，使其在目前實際工業(yè)問題中無法應(yīng)用。

最初，有的學(xué)者為了降低LES的計算量，嘗試在粗網(wǎng)格下應(yīng)用LES(Very large LES-VLES)方法[2]，沒有獲得成功。但其以最小代價捕捉流場中大尺度的非定常現(xiàn)象和相干結(jié)構(gòu)的思想是非常誘人的，這也成為許多混合LES/RANS方法發(fā)展的目的。Hybrid LES/RANS方法就是為了搭建在工程應(yīng)用環(huán)境下RANS方法和LES方法所存鴻溝之間的橋梁，克服RANS模型對復(fù)雜分離流動模擬的能力不足和LES巨大計算量的問題，分別將各自的優(yōu)勢進(jìn)行巧妙的組合應(yīng)用于解決中等/大分離條件下的流動模擬問題。

基于非結(jié)構(gòu)網(wǎng)格的UNSMB計算平臺是具備大規(guī)模并行計算能力，包含多種二階精度的空間和時間格式及大量湍流模型的內(nèi)部代碼。當(dāng)前在該平臺的解算器框架下開展了以下工作：通過對湍流模擬方程的求解形式進(jìn)行改造，使其在同一的求解框架下既具備RANS分支，也包含LES分支，其湍流模擬方式的切換是通過對計算域中湍流尺度的判斷來完成的。本次混合LES/RANS方法的研究主要集中在基于DES類的SA-DES和SA-DDES的模塊開發(fā)，使其具備復(fù)雜大分離流動的模擬能力。下面將圍繞這兩個混合模型的模塊開發(fā)和相應(yīng)的算例驗證進(jìn)行介紹。

1 混合LES/RANS模擬方法介紹

1.1 基于Navier-Stokes方法的解算器

UMSMB求解平臺是基于非結(jié)構(gòu)對偶網(wǎng)格、采用格點有限體積法，求解可壓縮Euler/NS方程的CFD求解器。其解算的方程組可以簡單寫為

(1)

式(1)中，C和D分別代表對流和粘性通量項，而S表示源項；式(1)中的解矢量q為守恒變量，q=[ρ,ρu,ρv,ρw,ρE,ρφ]T，其中，ρ為密度，u、v、w為3個方向上的速度分量，E為總能量，φ代表所求解的湍流變量。

UNSMB求解器可應(yīng)用于非搭接的任意網(wǎng)格單元形式，為基于邊的迭代求解方式，通過代碼的前處理程序?qū)⒔Y(jié)構(gòu)、非結(jié)構(gòu)和雜交網(wǎng)格中基于單元點的信息處理成UNSMB能辨識的基于邊的信息系統(tǒng)，從而構(gòu)造了以原始網(wǎng)格格點為中心的離散控制體。

將式(1)所表示的基本控制方程對任意控制體(Ω)及其的邊界(S)進(jìn)行積分，得到基本的計算格式為

(2)

時間相關(guān)的推進(jìn)格式是通過雙時間步長方法[5]來完成，包含一個全局的物理時間步長(Δt)和一個應(yīng)用于三階龍格庫塔格式內(nèi)迭代的虛擬當(dāng)?shù)貢r間步長(Δτ)。將式(2)所表示的流動控制方程改寫為時間導(dǎo)數(shù)項和和殘值項R(q)，它包括對流、擴散通量項和源項，即R(q)={C(q)+D(q)-Sq}。通過引入虛擬時間(τ)，非定常的N-S方程改寫為

(3)

在當(dāng)前的計算中，物理時間導(dǎo)數(shù)的離散通過三點后差的格式得到，其二階精度的計算格式為

(4)

在每個時間步Δt中，都是通過三階龍格庫塔推進(jìn)算法來迭代求解，但物理擴散項僅在第一步進(jìn)行計算。虛擬時間步長的選擇需要受到物理真實步長的限制，為

(5)

此處，Δτ0表示為滿足穩(wěn)定性問題所確定的當(dāng)?shù)貢r間步長；CFLv是粘性CFL數(shù)，為UNSMB中用戶自己輸入的參數(shù)(CFLVIS)。此外，還通過隱式殘值光順和多重網(wǎng)格來加速收斂，并通過MPI通信技術(shù)來實現(xiàn)大型并行計算。

1.2 Spalart-Allmaras RANS模型

(6)

此處不考慮轉(zhuǎn)捩相關(guān)項，至于詳細(xì)的模型常數(shù)和經(jīng)驗函數(shù)就不在此贅述了。

1.3 SA-DES模型

大渦模擬在模擬湍流的細(xì)節(jié)和流動分離等方面都有很好的精準(zhǔn)度，但其對高雷諾數(shù)下物面流動的計算量是當(dāng)前計算機水平所不能企及的，故其距工程應(yīng)用還有相當(dāng)大的距離。RANS方法對湍流的模擬是通過湍流模型來完成的，不能解析任何湍流細(xì)節(jié)，反映的是平均流動的變化情況，但RANS方法的高效性和魯棒性使其在工業(yè)界得到廣泛的應(yīng)用。并且兩種模型方法的求解方程相似度很高，以雷諾平均速度〈ui〉的輸送方程為例，即

(7)

(8)

通過對比可以明顯看出公式(7)和公式(8)非常相似，可以在同一解算器中應(yīng)用相同的程序框架來實現(xiàn)。再對LES方法中應(yīng)用渦粘性假設(shè)的SGS湍流模型進(jìn)行分析會發(fā)現(xiàn)大多數(shù)的SGS模型的推導(dǎo)都來源于RANS模型，為RANS模型的對應(yīng)物。因此，LES方法和RANS方法不僅主控方程結(jié)構(gòu)相似，而且湍流模型也相似。于是有人提出了分別應(yīng)用RANS和LES優(yōu)勢的混合LES/RANS模型方法，即在近物面應(yīng)用RANS模型，遠(yuǎn)離物面的分離區(qū)域應(yīng)用LES模型進(jìn)行解析湍流。

圖1 DES方法的示意草圖

Spalart等人通過對SA模型進(jìn)行如下的修改完成了最早DES格式的設(shè)計(如圖1所示)，也即被廣泛稱之為DES97格式[7]。

Δ=max{Δx;Δy;Δz}

(9)

通過各項同性湍流算例的校準(zhǔn)得到模型參數(shù)Cdes=0.65，有文獻(xiàn)研究表明數(shù)值結(jié)果對該參數(shù)不是很敏感。

(2)對模型方程中長度尺度進(jìn)行了如下修改

(10)

圖2 DES方法中RANS和LES模式切換圖示

如圖2所示：(1)近壁面的長度尺度的增長不能超出RANS的計算區(qū)域，在貼近物面處dCdesΔ，模型將轉(zhuǎn)化為SGS模式。湍流的渦粘性系數(shù)υt的過渡是連續(xù)和光滑的，這是因為兩種模式的變化僅僅是流動附加方程中源項的平滑變化。事實上，RANS模型項的應(yīng)用并沒有阻礙流動變化為非定常流動，即使是近壁面DES流場解在統(tǒng)計上是穩(wěn)態(tài)也基本上會轉(zhuǎn)化為非定常，這是由外場流動的非定常脈動所導(dǎo)致。

Spalart等人創(chuàng)造“脫體渦模擬”(Detached Eddy Simulation-DES)術(shù)語來命名這種方法，其主要目的是將流場的求解方法進(jìn)行混合，在外部流動區(qū)域應(yīng)用LES，應(yīng)用LES模式來解析遠(yuǎn)離所有邊界的脫體渦，與此同時在近壁面流動區(qū)域應(yīng)用RANS模型。應(yīng)用RANS模型的目的是為了得到近壁面流動在統(tǒng)計意義上的合理描述，其計算網(wǎng)格僅需在物面法向加密，而在切向可應(yīng)用相對較粗的網(wǎng)格。此外，這兩種方法的切換是通過一個自動的判據(jù)來完成，從而避免需要用戶事先指定。

公式(9)的提出是基于在LES方法中最粗網(wǎng)格步長決定了可解析渦結(jié)構(gòu)最小尺度。這個問題在LES方法中比混合LES/RANS方法更突出，它影響LES/RANS的交界面。文獻(xiàn)[8]中對公式(9)式進(jìn)行了詳細(xì)的討論，圖1～2所示網(wǎng)格類型是一種針對DES方法來說比較安全的網(wǎng)格形式。

1.4 SA-DDES模型

自1997年Spalart等人提出DES方法以來，混合LES/RANS方法就一直被當(dāng)作工程CFD計算中解決復(fù)雜漩渦運動和分離流動問題的最合適的手段，引起了工業(yè)界和學(xué)術(shù)界的廣泛研究，其不完善的地方還有很多，遠(yuǎn)未成熟。

早期DES方法在薄的邊界層可以很好地模擬真實的流動現(xiàn)象，但在模擬厚邊界層和淺分離區(qū)域時可能會出現(xiàn)一些錯誤的行為，主要出現(xiàn)在當(dāng)平行于物面的網(wǎng)格間距Δ‖小于邊界層厚度δ時(如圖3所示)，這種現(xiàn)象既可能是由網(wǎng)格加密引起也可能是邊界層本身就比較厚引起的。當(dāng)DES網(wǎng)格加密使得長度尺度進(jìn)入LES分支時(渦粘性低于RANS的水平)，由速度脈動分辨得到的解析雷諾應(yīng)力不能替代模化的雷諾應(yīng)力。雷諾應(yīng)力的虧空[8]降低了表面的摩擦力，可能導(dǎo)致由網(wǎng)格所誘導(dǎo)的流動提前分離。

為了避免在附著流中出現(xiàn)非物理現(xiàn)象，消除DES方法對網(wǎng)格的限制，在邊界層內(nèi)保存RANS特性，物面的網(wǎng)格間距Δ‖與邊界層厚度δ無關(guān)，減小灰色區(qū)域。Spalart等人[9]通過加入一個函數(shù)fd來重新定義耗散長度尺度

圖3 模糊的DES網(wǎng)格類型

(11)

此處

fd=1-tanh[(8rd)3]

(12)

及

(13)

2 算例測試和驗證

前面介紹了我院開發(fā)SA-DES和SA-DDES程序模塊的相關(guān)情況，下面將通過一系列的測試算例和工程算例來檢測和驗證該模塊的數(shù)值模擬能力。

2.1 NACA0012翼型的算例

翼型的動態(tài)失速是一個伴隨流動分離的典型非定常不穩(wěn)定過程，失速渦在前緣附近已經(jīng)形成，并沿著翼型表面移動和增長，最終在翼型后緣附近產(chǎn)生分離，脫離翼型表面，這種分離流動的動態(tài)變化影響了機翼升力的產(chǎn)生。

本節(jié)將采用SA-DES和SA-DDES來模擬翼型的動態(tài)失速運動，測試程序模塊的有效性，并與URANS的計算結(jié)果進(jìn)行對比分析。

2.1.1 計算模型及網(wǎng)格

本次測試計算應(yīng)用NACA12翼型作為計算模型，展向拉伸1倍弦長，弦長為0.1 m。使用POINTWISE生成O-型結(jié)構(gòu)網(wǎng)格，網(wǎng)格數(shù)量為225×191×31，展向均布31層，為了精確捕捉附面層內(nèi)的流動，網(wǎng)格第一層高度為1×10-6，y+均小于1。圖4為部分網(wǎng)格視圖。

圖4 NACA0012翼型的表面及空間網(wǎng)格

2.1.2 數(shù)值模擬的時間歷程分析

本算例將開展NACA0012翼型在45°失速攻角下的研究，并將計算結(jié)果與試驗數(shù)據(jù)[10,11]和Hong-silk IM[12]的DDES計算結(jié)果進(jìn)行對比分析。計算狀態(tài)為Ma=0.50，Alpha=45°，基于弦長的雷諾數(shù)Re=1.3×106。分別進(jìn)行了基于URANS、DES和DDES的3種非定常模擬方法研究，其物理時間步長均設(shè)置為0.02T0，為1.2×10-5，T0為來流的特征時間尺度。由于DDES模擬方法對初場敏感，因此，本次計算中均以自由來流作為初場值，均大致計算400T0時長。

圖5描述了UNSMB代碼DDES數(shù)值計算的殘值收斂過程。該圖表明計算的收斂性很好，呈現(xiàn)了明顯的周期性，其表示為上翼面漩渦的脫離周期，不過該 DDES的幅值隨時間變化，有一定的混沌現(xiàn)象。

圖5 DDES計算的密度收斂歷程

2.1.3 數(shù)值模擬的時均結(jié)果分析

圖6表明3種模擬結(jié)果從100T0開始對瞬態(tài)阻力系數(shù)進(jìn)行平均，直至400T0，其時均阻力系數(shù)最后均趨于常數(shù)，不過URANS的阻力波動幅度大，時均值變化較為緩慢，而DES和DDES的變化速度較快，大致在250T0時就基本不變。NACA0012翼型45°攻角下試驗[10,11]的時均阻力為1.075， UNSMB計算的為CDURANS=0.979，CDDES=1.110，CDDDES=1.078。通過與試驗結(jié)果的對比表明,DDES計算結(jié)果與試驗值的吻合度非常高，僅高0.28%，而DES結(jié)果大了3.3%，且URANS的計算結(jié)果與試驗值的差值也不是特別多，充分表明UNSMB代碼對大分離流動的模擬準(zhǔn)度較好，其中DES模擬結(jié)果的稍差可能是由于翼型表面的弦向網(wǎng)格分布過密所導(dǎo)致。

圖6 時均阻力系數(shù)比較

2.1.4 數(shù)值模擬的瞬態(tài)流場分析

下面將對3種模擬方法的瞬態(tài)流場進(jìn)行簡要地分析，以展向50%長度中截面作為研究對象，分別選取兩個典型時刻進(jìn)行比較。

圖7為3種模擬方法的瞬態(tài)流場圖，通過應(yīng)用湍流渦粘性系數(shù)對剖面進(jìn)行染色，反映了尾跡區(qū)的運動。通過3組圖像的對比可以看出,收斂后的URANS瞬態(tài)流動變化較小，更多地體現(xiàn)RANS的特性，而DES和DDES可以清晰地分辨出尾跡漩渦的運動過程，解析的漩渦尺度更小，流動是在翼型的后緣卷起再脫離上表面。此外，通過對等值面的量值比較可以發(fā)現(xiàn)URANS計算的湍流渦粘系數(shù)比DES和DDES大得多，而DDES的渦粘系數(shù)最小，這說明URANS的數(shù)值耗散太大，尾跡中流動的脈動量被快速抹平，但DES和DDES中較小的湍流渦粘性系數(shù)就可以較好地模擬尾跡區(qū)的流動特征。

圖7 URANS(A、B)和DES(C、D)及DDES(E、F)的瞬態(tài)流場圖

2.2 鈍前緣三角翼的算例驗證

在2001年，北約RTO AVT-113工作團(tuán)隊提出了第二屆“國際漩渦流動試驗”(VFE-2)[13]的研究項目，進(jìn)行了大量的試驗和CFD計算研究，目的是為了拓展對三角翼流動拓?fù)浜臀锢憩F(xiàn)象的認(rèn)識及對計算代碼的驗證和確認(rèn)。下面我們將通過對中等前緣鈍度的VFE-2構(gòu)型進(jìn)行DES計算來驗證代碼的可靠性。

2.2.1 計算模型和網(wǎng)格

VFE-2模型為一個帶尾撐的薄三角翼，包含4個可換的前緣(1個尖前緣和3個鈍前緣)，如圖8所示。

圖8 計算模型的正視圖

根據(jù)國外對VFE-2的研究表明，1200萬單元的網(wǎng)格就足夠的。本次計算使用ICEM CFD生成結(jié)構(gòu)網(wǎng)格，第一層高度為1×10-6，附面層內(nèi)網(wǎng)格保持在30～50重左右；大部分y+都在1左右，僅很少部分的y+>1，且最大值為2.9。

圖9 計算網(wǎng)格全局視圖和對稱面部分視圖

2.2.2 DES模擬結(jié)果

本算例主要使用SA-DES來進(jìn)行驗證計算，其計算狀態(tài)為Ma=0.4，Alpha=13.3°，基于平均氣動弦長的雷諾數(shù)Rmac=3×106，其時間步長設(shè)置為Δt=2.5×10-5。該計算狀態(tài)為高亞音速小分離流動狀態(tài)，在距離頂點一段距離的位置流動產(chǎn)生分離，且二次渦的旋轉(zhuǎn)方向與主渦相同。模擬該類流動最大的難點在于前緣分離點的位置，以及主渦及二次渦的展向位置及大小。

圖10為對DES瞬態(tài)計算結(jié)果在某一時間歷程內(nèi)進(jìn)行時均后的結(jié)果分析，其流動結(jié)構(gòu)和定常計算結(jié)果類似，右側(cè)圖可以看出內(nèi)側(cè)二次渦較早的形成，但強度很低。

圖10 SA-DES計算得到上表面Cp分布和等截面的總壓恢復(fù)等值線圖

2.2.3 流動拓?fù)浣Y(jié)構(gòu)分析

圖11顯示的是亞聲速下三角翼的壓力云圖和表面摩擦力線圖，該圖可以增加我們對鈍前緣三角翼流動拓?fù)涞睦斫?。分離線是從前緣開始的，但并沒有延伸到機翼的頂點位置，一旦分離產(chǎn)生，自由剪切層就對流至機翼上表面，從而在近物面得到一個壓力較低和旋度較強的區(qū)域，導(dǎo)致了外側(cè)主渦的產(chǎn)生。

圖11 表面極限流線和壓力分布

3 結(jié)論

本文主要闡述了UNSMB代碼的混合LES/RANS模塊的開發(fā)研究，其湍流模型為SA-DES和SA-DDES，并分別介紹了混合模塊相關(guān)的理論和驗證算例分析。計算結(jié)果表明,UNSMB的DES和DDES模塊可以應(yīng)用于大分離流動及復(fù)雜漩渦運動的數(shù)值模擬計算，并具備良好的計算精準(zhǔn)度，可以解析和捕捉由流場內(nèi)部自身所誘導(dǎo)的強非定常運動和漩渦結(jié)構(gòu)。該混合模塊在CFD理論研究和工程計算中都具備廣闊的應(yīng)用前景，可以滿足該項目所要求的飛機結(jié)冰后氣動力性能評估的要求。