王圣業(yè) 符 翔 楊小亮 鄭皓榜 鄧小剛?

1 國防科技大學(xué)空天科學(xué)學(xué)院,長沙410073

2 軍事科學(xué)院,北京100190

1 引 言

湍流是流體在空間和時間上都具有急劇不規(guī)則和高度隨機性脈動的一種流動狀態(tài),在自然界和工程技術(shù)各領(lǐng)域中極為常見(鄭曉靜和王國華2020).湍流問題作為流體力學(xué)中有名的難題,經(jīng)過包括許多偉大科學(xué)家在內(nèi)的學(xué)者長達100 多年的努力,已形成獨特的研究體系(是勛剛1992).然而到今天,距離徹底弄清和完全精準地預(yù)測湍流仍有很長的路.流體力學(xué)中,對湍流問題的研究和預(yù)測包含三大手段: 理論、實驗和數(shù)值計算.隨著20 世紀70 年代計算機技術(shù)的飛速發(fā)展,數(shù)值計算(計算流體力學(xué),CFD) 在研究湍流運動規(guī)律,解決實際湍流問題中發(fā)揮著越來越重要的作用(Slotnick et al.2014).

湍流的數(shù)值模擬方法一般包括三大類: 雷諾平均模擬(RANS) 方法、直接數(shù)值模擬(DNS) 方法和大渦模擬(LES) 方法.DNS 方法不添加模型而試圖直接求解所有尺度的湍流.其網(wǎng)格分辨率及最大允許時間步長必須足夠小,從而捕捉最小湍流尺度(Kolmogorov 尺度).而RANS 方法則僅關(guān)心最大尺度湍流,通過建立模型去?；型牧髯V的特征.然而由于湍流的復(fù)雜性,很難尋找到一種尺度去近似模擬所有的流動情況(Wilcox 2006).LES 試圖僅對最小尺度湍流進行?；?來消除DNS 求解宏觀湍流問題時不必要的花費,因此相比RANS 方法可以建立更加簡潔和普適的模型.理論上,為了精準可靠地預(yù)測湍流,最好是直接將所有湍流波動求解直到最小尺度(即DNS).然而,由于湍流最小和最大尺度的比與湍流雷諾數(shù)有關(guān)(λ/L=),使得所研究湍流問題的雷諾數(shù)會受到計算機硬件水平的限制.因此DNS 目前僅用于低雷諾數(shù)的湍流機理研究,在湍流結(jié)構(gòu)生成和演化等方面取得了巨大的成就(傅德薰等2010).

在工程技術(shù)領(lǐng)域,RANS 方法扮演了中流砥柱的作用.RANS 方法的核心是對雷諾應(yīng)力進行封閉建模,其思想可追溯到一百多年前Boussinesq (1877) 提出著名的渦黏性假設(shè).其后,諸多流體力學(xué)先驅(qū)發(fā)展了一系列半經(jīng)驗理論(陳懋章2002).但由于只考慮了一階湍流統(tǒng)計量,這些渦黏模型也被稱為一階矩模型.過去的三十多年間,渦黏模型在預(yù)測穩(wěn)態(tài)流動問題方面獲得了很大的成功.典型代表為在航空工程領(lǐng)域得到廣泛應(yīng)用的Sparlart-Allmarsa(SA)一方程模型(Sparlart 1994,2007) 和剪切應(yīng)力輸運(SST) 兩方程模型(Menter 1994,2003).然而當流場中出現(xiàn)明顯的流動分離時,傳統(tǒng)渦黏模型往往會表現(xiàn)出不足,進而無法準確預(yù)測分離流動下飛行器的氣動特性.比如,在分離渦處會高估渦黏性,無法分辨不同尺度的渦結(jié)構(gòu)(周鑄等2017); 在部件結(jié)合處無法估計湍流正應(yīng)力各向異性,形成虛假分離(Rumsey et al.2016) 等.物理上更為完善的封閉建模方法是直接對湍流應(yīng)力項建立微分輸運方程.該方法最早可追溯到1940 年由我國著名力學(xué)大師周培源所建立的雷諾應(yīng)力輸運微分方程組(周培源1940).由于雷諾應(yīng)力屬于二階湍流統(tǒng)計量,而在建立輸運方程組的過程中會出現(xiàn)三階(更高階) 統(tǒng)計量.如果在二階范圍內(nèi)對更高階統(tǒng)計量建模,則稱為二階矩模型;如果對三階統(tǒng)計量繼續(xù)建立微分輸運方程,則為三階矩模型;依此類推.目前,學(xué)術(shù)界普遍將二階矩及其以上的湍流封閉模型稱為高階矩模型(Hanjali′c & Launder 2011).

在RANS 方法框架下,雷諾應(yīng)力(輸運)模型(RSM)主要考慮以二階矩封閉為主.Rotta (1951)進一步發(fā)展了周培源的工作,提出了完整的RSM 模型.其后,Donaldson 和Rosenabum (1968)提出了模型不變性概念,Lumley (1978) 發(fā)展了封閉性近似等,但都處于初級研究階段.到了20 世紀90年代,計算機水平的飛速發(fā)展和一階矩模型的不足使得人們又重新開始重視二階矩模型.Launder(1989) 對過去二階矩模型的發(fā)展進行了總結(jié),并認為在隨后幾年中將一些模型納入商業(yè)軟件是相當可行的.Speziale 等(1994)認為RSM 方法給出了用于復(fù)雜湍流流動計算的LES 方法和RANS 方法之間的橋梁.他進一步討論了雷諾應(yīng)力模型的可實現(xiàn)性問題,并提出了簡化的設(shè)計準則.其后,諸多學(xué)者基于上述準則,對雷諾應(yīng)力輸運方程中的未知項(如湍流耗散項、壓力應(yīng)變關(guān)聯(lián)項和擴散項) 的建模進行了研究.而在這其中,最為核心并且困難的是壓力應(yīng)變關(guān)聯(lián)項的建模.一方面由于壓力應(yīng)變關(guān)聯(lián)項的量級和生成項相當,其在很多實際流動問題中十分重要,并且也是區(qū)別于渦黏模型的關(guān)鍵(Wilcox 2006); 另一方面,由于壓力應(yīng)變關(guān)聯(lián)項很難通過實驗測量得到,如何構(gòu)造一個合理的模型需要高超的創(chuàng)造力.近年來,包括美國航空航天局(NASA)、德國宇航院(DLR)等機構(gòu)紛紛提出重點發(fā)展二階矩雷諾應(yīng)力模型的計劃.例如NASA 蘭利研究中心在2012 年指出了幾個關(guān)鍵研究熱點(Malik 2012): “改進RANS 結(jié)果計算機制,使二階矩封閉方法可實際應(yīng)用于剛性方程問題; 在雷諾應(yīng)力方程方面,改進長度尺度方程、壓力應(yīng)變模型、壁面處理以及壓縮性影響、動態(tài)激波相互作用效應(yīng).”而在2014 年又進一步指出(Slotnick et al.2014): “對RSM 模型研究工作的投資包括模型的不同變種在CFD 軟件中的應(yīng)用、驗證和確認.” DLR 更是走在前列,發(fā)展出的JHh-v2 雷諾應(yīng)力模型 (Jakirli′c & Hanjali′c (2002) 提出的 JHh 模型改進版) 基于先進的近壁耗散修正,對近壁湍流尤其是壓力梯度誘導(dǎo)分離等效應(yīng)具有很好的模擬效果;通過歐洲FLOMANIA 計劃發(fā)展的SSG/LRR-ω模型(Eisfeld&Brodersen 2005)借鑒了Menter 的混合函數(shù)思想,具有良好的魯棒性,在一些航空工程問題中取得了較好的模擬效果(Togiti et al.2014,Rumsey 2015,Lee-Rausch et al.2016,董義道等2016,Wang et al.2018).我國學(xué)者也在雷諾應(yīng)力模型發(fā)展方面提供了諸多開創(chuàng)性工作,例如符松教授提出的Fu-Launder-Tselepidakis 模型等(Fu 1987,1998,符松1994); 聶勝陽等將 SSG/LRR-ω模型與γ-Reθ轉(zhuǎn)捩模型結(jié)合 (Nie 2017),擴展了轉(zhuǎn)捩模型范圍; 作者等提出了基于新型渦黏尺度的雷諾應(yīng)力模型(EV-RSM),在魯棒性方面獲得了明顯提升,以此成功實現(xiàn)了高階矩模型的高精度離散及穩(wěn)定求解(Wang et al.2020).

LES 被普遍用于非穩(wěn)態(tài)流動的精細模擬,并且隨著計算機水平的發(fā)展而逐漸進入工程領(lǐng)域.諸多學(xué)者及機構(gòu)提出LES 方法作為下一代CFD 軟件核心的預(yù)測或計劃(Slotnick et al.2014).LES方法的核心是對亞濾波(或亞格子) 尺度應(yīng)力進行封閉建模,其思想首先由氣象學(xué)家Smagorinsky(1963) 提出.由于僅對最小尺度湍流進行?；?LES 框架下的湍流模型相比RANS 框架下更加普適.再加上LES 巨大的計算資源要求,使得模型被構(gòu)造得更為簡潔.目前,在實際流動問題中廣泛被應(yīng)用的亞濾波模型均是一階矩模型,通過Boussinesq 的渦黏性假設(shè)與濾波NS 方程建立聯(lián)系.然而,在近年來的一些精細流動研究中,很多學(xué)者指出了傳統(tǒng)渦黏性模型存在的不足(Georgiadis et al.2010).與 RANS 方法類似,LES 方法中也存在高階矩湍流模型.早在 1972 年便由氣象學(xué)家Deardoff (1973,1974) 提出并在氣象流動中進行了嘗試,但他的工作大大地走在了時代的前面.近年來計算機硬件水平及湍流結(jié)構(gòu)識別技術(shù)的進一步發(fā)展,為二階矩亞濾波模型的研究提供可能(Stoellinger et al.2015).

湍流模擬方法正經(jīng)歷從RANS 到LES 的關(guān)鍵性轉(zhuǎn)變,然而該過程也并非能一蹴而就(Slotnick et al.2014).工程問題中的特征雷諾數(shù)普遍較高,即使附著邊界層內(nèi)的最大尺度也會變得很小,采用LES 模型對網(wǎng)格尺度的要求也并不比DNS 減弱太多.如果在邊界層附近采用各向異性的模型,如RANS 模型,而在遠離壁面的區(qū)域采用LES 模型,則能大大降低對計算機硬件的要求.這便是混合RANS/LES 方法.Sparlart (1997) 基于SA 湍流模型提出了分離渦模擬方法(DES),是在實際工程中準確預(yù)測非定常湍流的里程碑工作.這之后,各種形式的混合RANS/LES 方法如雨后春筍般迅速發(fā)展(Frohlich & von Terzi 2008),同時在一定程度上被應(yīng)用于真實工程分離流動模擬當中,取得了令人鼓舞的成果.然而,受限于線性渦黏假設(shè),基于傳統(tǒng)RANS 方法的混合形式出現(xiàn)了諸多不足(Greschner et al.2008).正如NASA 報告(Slotnicket al.2014) 指出的: 首先,模擬任何從邊界層內(nèi)部發(fā)展出來的分離流動,都需要RANS 方法的不斷改進; 其次,為提高方法的魯棒性,在邊界層內(nèi),RANS 和 LES 需要自動進行光滑過渡.因此,改進混合 RANS/LES 方法,尤其是 “灰區(qū)” 問題,關(guān)鍵仍是RANS 模型.而高階矩模型在對過渡區(qū)弱分離流的捕捉方面具有比渦黏模型更大的優(yōu)勢.如Chaouat (2011) 將亞濾波尺度應(yīng)力模型與部分積分輸運模型(PITM) 結(jié)合; Probst (2011)將DES 思路引入εh-RSM 模型,在翼型尾跡小分離模擬中展示了相對渦黏性模型DES 的優(yōu)勢; 作者等(Wang et al.2017) 將IDDES 引入SSG/LRR-ω模型,在翼型和三角翼失速模擬中取得了較好的效果.

圖1

綜上,高階矩模型是湍流模式理論發(fā)展中舉足輕重的部分.尤其是進入新世紀后,隨著計算機硬件水平和數(shù)值求解技術(shù)的巨大提升,高階矩模型又迎來了發(fā)展的小高潮.多個國外科研機構(gòu)提出過復(fù)興高階矩模型的計劃,尤其是美國NASA 在CFD 的2030 愿景中(Slotnick et al.2014) 將高階矩模型作為RANS 部分研究和資助的主要方向,并將2020 年設(shè)置為關(guān)鍵里程碑節(jié)點(如圖1).基于此,本文對高階矩模型的發(fā)展情況進行了總結(jié),包括高階矩模型各項的建模過程、尺度提供方程的演化歷程、數(shù)值求解技術(shù)的研究需求、其相對于傳統(tǒng)渦黏模型的優(yōu)勢以及部分CFD 軟件的集成情況等.拋磚引玉,以求在促進我國的湍流模式理論研究中做出一點貢獻.

2 湍流量輸運方程

2.1 二階矩模型一般框架

如引言所述,周培源先生首次建立了一般湍流的雷諾應(yīng)力所滿足的微分輸運方程組.但在該方程組中又出現(xiàn)了未知的三階統(tǒng)計量(周培源1940).為了使方程封閉,周培源基于自由湍流條件引入了一些假設(shè),得到了二階矩模型的雛形(Chou 1945).其后,Rotta (1951)發(fā)展了周培源的思路,提出了更一般的二階矩模型.至今,RANS 方法下的二階矩模型基本都是基于他們的工作發(fā)展而來.本文設(shè)雷諾應(yīng)力項(“?”代表Favre 平均),六分量方程具有如下通用形式(Cecora et al.2015)

其中右端分別為生成項、壓力與應(yīng)變關(guān)聯(lián)項(再分配項)、耗散項、擴散項及質(zhì)量項.其中,生成項可精確得到

壓力與應(yīng)變關(guān)聯(lián)項、耗散項和擴散項在二階封閉范圍內(nèi),均是雷諾應(yīng)力、湍動能、耗散率、平均密度、平均速度等的函數(shù).而質(zhì)量項通常被忽略.

與RANS 方法中的二階矩模型類似,LES 方法中的二階矩模型(或稱為亞格子應(yīng)力模型)早在1973 年代就被Deardorff 提出,并被用于氣象流動.本文設(shè)亞格子應(yīng)力為則一般框架表示為(Chaouat et al.2017)

等式右端仍然包含生成項、壓力與應(yīng)變關(guān)聯(lián)項(再分配項)、耗散項、擴散項等.其中生成項也是通過亞格子應(yīng)力精確獲得.而其余右端項需要建模處理.

由于 LES 方法下的二階矩模型與 RANS 方法下有著相同的結(jié)構(gòu),如果引入過渡函數(shù)F使RANS 方法作用于邊界層附近,而在遠離邊界層區(qū)采用LES 方法,則可構(gòu)造混合RANS/LES 框架下的二階矩模型

作者曾基于Sanchez-Rocha 和Menon (2009) 在可壓縮流動中得到的混合RANS/LES 框架,通過選取SSG/LRR-ω模型作為RANS 部分,Deardorff 模型作為LES 部分,Menter 的F2函數(shù)作為過渡函數(shù)對塔形凸起物繞流進行了模擬.圖2 表明混合RANS/LES 框架下的二階矩模型具備準確分辨渦流動的潛力.

2.2 壓力與應(yīng)變關(guān)聯(lián)項建模

上述二階矩模型中除了生成項可以在二階湍流量范圍內(nèi)精確表示,其他右端項均需建立模型處理.而其中最關(guān)鍵的一項是壓力與應(yīng)變關(guān)聯(lián)項的建模.首先,該項與生成項階數(shù)相同,在大多數(shù)復(fù)雜湍流流動中扮演著關(guān)鍵角色,也是雷諾應(yīng)力模型區(qū)別于渦黏模型的核心.其次,由于該項涉及壓力脈動,而試驗中很難對其觀測,因此建立一個合理的封閉需要很大程度的創(chuàng)造性(Wilcox 2006).

圖2

Rotta (1951) 通過對壓力脈動泊松方程的推導(dǎo),將壓力與應(yīng)變關(guān)聯(lián)項分解成慢項和快項

慢項主要決定湍流是否會趨于各向同性,通過定義一個無量綱來度量雷諾應(yīng)力各向異性

并認為慢項直接與該度量成正比,即Aij=?C1εaij.Lumley 和Khajeh-Nouri (1974) 進一步提出了非線性的慢項表達式,可視為在線性項后添加了aij的張量積項

對于快項,在 1970 年代前學(xué)者普遍認為其遠小于慢項,可以忽略 (Donaldson & Rosenbaum 1968,Daly & Harlow 1970).但其中也有不同的聲音,如Crow (1968) 和Reynolds (1970) 提供的湍流例子表明快速壓力應(yīng)變的影響遠遠大于慢壓力應(yīng)變.Launder,Reece 和Rodi (1975) 完全基于運動學(xué)的考慮,設(shè)計了一種巧妙的封閉逼近,即著名的LRR 模型(慢項仍基于Rotta 的線性模型)

其中

該封閉最顯著的特征之一是只存在一個封閉系數(shù),即0.40.6.同時,LRR 還提出了一種簡化版本

為便于區(qū)分,式 (8) 通常稱為 LRR-QI 模型,式 (10) 稱為 LRR-IP 模型.在計算壁面附近流動時,LRR 模型無法考慮壓力– 反射效應(yīng),需要進行修正.包括 Gibson 和 Launder (1978) 以及Craft 和Launder (1992) 的修正形式.

另一方面,Lumley (1978) 從可實現(xiàn)性的角度出發(fā),提出了Lumley 壓力與應(yīng)變關(guān)聯(lián)模型.其后,Speziale,Sarkar 和Gatski (1991) 提出了類似但更簡單的非線性模型,即SSG 模型

從形式上看,SSG 模型包括LRR 模型中的相關(guān)項,即,C2,為零.尤其是,SSG 模型可以不需要對壓力– 反射效應(yīng)進行修正,因此在壁面附近比LRR 表現(xiàn)更好.同一時期,還有Shih-Mansour-Chen 模型(Shih et al.1987),Fu-Launder-Tselepidakis 模型(Fu et al.1987)等,重點也均是對壓力與應(yīng)變關(guān)聯(lián)項的非線性封閉.

Wilcox (2006) 通過研究指出,LRR 模型的壓力– 反射問題與耗散尺度ε有關(guān).通過將ε尺度替換為ω尺度,可避免在壁面處添加有關(guān)修正,即 Stress-ω模型.由于基于ω的 LRR 模型與基于ε的 SSG 模型均無需壓力 – 反射修正,因此數(shù)學(xué)形式十分相似.Eisfeld 等 (2005) 將 Menter 在BSL(或 SST) 模型中的過渡函數(shù)引入到了壓力與應(yīng)變關(guān)聯(lián)項的建模中,得到了 SSG/LRR 混合模型.該模型與式(11) 相同,只是系數(shù)Ci和通過過渡函數(shù)F1獲得

圖3 展示了SSG/LRR-ω模型的構(gòu)造思路,其兼具SSG 模型和LRR 模型的優(yōu)勢,在航空復(fù)雜流動等實際應(yīng)用中獲得了良好的評價(Togiti et al.2014,Rumsey 2015,Lee-Rausch et al.2016,董義道等2016,Wang et al.2018).

圖3

LES 方法中的二階矩模型,關(guān)鍵也是對于壓力與應(yīng)變關(guān)聯(lián)項的封閉建模.Deardoff(1973,1974)提出的亞格子應(yīng)力模型與式(5) 類似,該項也被分為慢項和快項兩部分

第一項系數(shù)cm與雷諾應(yīng)力模型中Rotta 系數(shù)C1作用相同.第二項是針對均勻應(yīng)變湍流在初始各向同性狀態(tài)時獲得的.相對于RANS 框架下的壓力與應(yīng)變關(guān)聯(lián)項,LES 框架下的封閉要簡陋很多.這也是由于計算機硬件以及試驗觀測的限制造成的.當前,隨著DNS 的發(fā)展給壓力與應(yīng)變關(guān)聯(lián)項的研究和標定提供了新的途徑.如Jakirli′c 和Hanjali′c (2002)提出了改進近壁湍動能及應(yīng)力耗散模擬的JHh-v2 雷諾應(yīng)力模型,并利用完全發(fā)展的槽道流動(Rem= 5600) 進行了標定.Stoellinger 等(2015) 基于Hanjali′c等的工作得到了一種亞格子應(yīng)力模型,并在低雷諾數(shù)槽道流動中與DNS 結(jié)果進行了對比分析.

2.3 擴散項建模

式(1) 中擴散項包含兩部分: 湍流擴散和分子擴散.分子擴散與NS 方程一致,通常利用牛頓黏性應(yīng)力公式得到(陳懋章2002).而湍流擴散含有壓力脈動相關(guān)量和三階速度脈動關(guān)聯(lián)量,需要進行建模.在湍流擴散建模中最常用的方法是假設(shè)一個梯度輸運過程.Daly 和Harlow (1970) 認為,雷諾應(yīng)力的湍流擴散速度應(yīng)與雷諾應(yīng)力的梯度成正比,于是有

其中k2/ε為量綱決定系數(shù),由 Mellor 和 Herring (1973) 采用.這種形式在數(shù)學(xué)上很簡單,同時數(shù)值穩(wěn)定也很好,尤其在復(fù)雜外形應(yīng)用中獲得了成功.但是,該形式不符合旋轉(zhuǎn)不變性.Donaldson(1972) 提出了對稱形式,其后Launder 等(1975) 對此進行了完善,提出了更一般的形式

2.4 耗散項建模

因為耗散發(fā)生在最小的湍流尺度,大多數(shù)二階矩模型使用了Kolmogorov(1941)的局部各向同性假設(shè),即

其中ε被稱為各向同性耗散率,與標準兩方程模型(k-ε或k-ω模型)類似,通過再建立尺度方程得到(參見本文第3 節(jié)).基于該形式,可自然利用k=?τij/2 關(guān)系得到標準的湍動能輸運模型.

湍流耗散的準確建模對于整個雷諾應(yīng)力模型至關(guān)重要.在至今很長的一段時期,關(guān)于雷諾應(yīng)力模型的改進主要集中在湍流耗散上.一方面由于耗散在現(xiàn)實中是各向異性的,特別是在壁面附近.Hanjali′c 和 Launder (1976) 提出了低雷諾數(shù)耗散形式

其中bij為無量綱雷諾應(yīng)力各向異性張量,即

同時,fs是一個低雷諾數(shù)阻尼函數(shù) (或稱為過渡函數(shù)),根據(jù)經(jīng)驗選擇它隨湍流雷諾數(shù)變化fs=(1+ReT/10)?1,ReT≡/(εˉν).其后,Kebede 等(1985)、Gilbert 和Kleiser (1991)、Launder 和Tselepidakis (1993) 等均對該形式進行了標定和發(fā)展.尤其是,Jakirli′c 和 Hanjali′c (2002) 在其 1993年版本(Hanjali′c & Jakirli′c 1993)的基礎(chǔ)上提出了近壁耗散修正雷諾應(yīng)力模型(JHh-v2 模型),將各向異性耗散率分解成均勻部分和非均勻部分

其中前者為均勻部分,仍由式(17) 得到; 后者是對應(yīng)應(yīng)力分量的黏性擴散,根據(jù)應(yīng)力輸運方程的解計算得出.該模型在翼型小分離流動中取得了良好的計算效果(Probst et al.2011).

在LES 框架下,耗散項與濾波/網(wǎng)格尺度相關(guān).在Deardorf f提出的亞格子應(yīng)力模型中使用了各向同性耗散模型

該式也被用于k方程亞格子模型(Sanchez-Rocha & Menon 2009).Stoellinger 等(2015) 基于Hanjali′c等工作得到的亞格子應(yīng)力模型,同樣采用了各向異性的耗散模型

該模型在壁面處耗散不為0,在低雷諾數(shù)槽道流動中獲得了良好的效果,如圖4.

引言中提到,受限于RANS 方法的局限,高階矩RANS 模型在模擬非定常流動中也需要針對耗散項進行修正.得益于雷諾應(yīng)力模型和亞格子應(yīng)力模型在各向同性耗散項上的相似,發(fā)展高階矩模型下的混合RANS/LES 方法可借鑒渦黏模型的做法.如Probst(2011)將DES 思路引入εh-RSM模型,Wang 等 (2017) 將 IDDES 方法引入 SSG/LRR-ω模型

圖4

圖5

2.5 更高階矩模型一般框架

在對二階矩模型開展研究和應(yīng)用的同時,人們也一直在探索通過更高階的展開,建立更為準確的湍流方程.然而受限于計算花費和數(shù)值剛性等問題,二階矩模型尚未能在工程中被廣泛采用,而更高階矩模型目前更僅限于學(xué)術(shù)研究.圖5 是Jeyapaul 等(2014)統(tǒng)計的更高階矩模型所需的動量輸運方程數(shù).可以看到,隨著維度的增加,微分方程數(shù)的增長是巨大的.

本文中對更高階矩模型的討論僅限于三階和四階,由于對更高階項的建模還很不成熟,下面僅給出一般框架.對照式(1),三階矩模型的一般框架為

右端項依次為生成項、壓力應(yīng)變關(guān)聯(lián)項、耗散項和擴散項,其中后三項需要在不高于三階湍流量的范圍內(nèi)進行建模處理.

可以證明,該假設(shè)對高斯速度分布是精確的.而參考二階矩模型,Kawamura 等(1995) 提出湍流擴散速度應(yīng)與三階矩應(yīng)力的梯度成正比,并依據(jù)DNS 數(shù)據(jù)給出了比例系數(shù).除此之外,Grossman 和Narayan (1993),Cheng 等(2005),Gryanik 等(2005) 均提出了相關(guān)模型.

其次是耗散項.高階矩湍流模型中的大多數(shù)耗散模型均基于各向同性假設(shè),同時對其進行近壁修正.高階耗散模型的一個特點是大于二階矩的耗散壁值為零.這形成了一個簡單的耗散邊界條件.然而,近壁漸近線為0 需要通過阻尼或其他方式建模(Jeyapaul 等2014).一個比較準確的耗散模型由Dekeyser 和Launder (1985) 提出,其采用了廣義梯度擴散假設(shè).然而由于三階耗散并非恒為負,該假設(shè)在簡單槽道流動中存在問題.后針對該問題發(fā)展了改進版.Kawamura 等(1995) 考慮了近壁情況,提出了更為復(fù)雜的耗散模型.

最后是壓力應(yīng)變關(guān)聯(lián)項.在三階矩模型下,它仍然是最為困難的部分.Dekeyser 和 Launder(1985) 曾參照二階矩的思路提出過相應(yīng)的模型,其同樣包括快項和慢項.其后,Kurbatskii 和Poroseva (1997),Kawamura 等(1995) 等基于該思路進行了發(fā)展.

四階矩模型一般框架為

其中右端的壓力應(yīng)變關(guān)聯(lián)項、耗散項和擴散項需要在不高于四階湍流量的范圍內(nèi)進行建模處理.尤其是在四階矩輸運方程中會出現(xiàn)五階相關(guān)項,可以假設(shè)其是二階動量矩和三階動量矩的函數(shù)(Jovanovi′c et al.1993).該假設(shè)認為概率分布函數(shù)為Gram-Charlier 級數(shù),即截斷近似關(guān)于高斯分布.表1 匯總了RANS 框架下二至四階封閉模型的各項研究情況.隨著階數(shù)的提高,封閉建模的難度往往也更大.

3 尺度方程

經(jīng)過上節(jié)涉及的建模過程,RANS 框架下的雷諾應(yīng)力方程中往往仍包含一個未知量,即各向同性耗散率ε.在物理上實質(zhì)是需要一個標量的湍流長度尺度或時間尺度來封閉整個模型.而在LES 框架下,這一尺度是通過濾波或網(wǎng)格提供的.因此,本小節(jié)的討論均是在 RANS 框架下開展的.考慮到耗散率的計算方式以及相應(yīng)的附加輸運方程所求解變量存在不同,將這一類用來封閉模型系統(tǒng)的方程統(tǒng)稱為湍流尺度方程.

在建立耗散率輸運方程時需要考慮所模擬流動的湍流雷諾數(shù),即Ret=k2/νε.由于湍流耗散發(fā)生在最小尺度上,對于高雷諾數(shù)流動通常采用Kolmogorov 的當?shù)馗飨蛲约僭O(shè)(Kolmogorov 1941) 來模擬耗散率張量,即εij=2/3εδij.Hanjali′c 和 Launder (1976) 基于ε的精確輸運方程 (Daly1970) 建立了標準耗散率模型方程,即

表1 高階矩RANS 模型各建模項研究匯總

其中方程右側(cè)依次為擴散項、生成項和破壞項.該式與k-ε模型中的耗散率方程相比存在兩點不同,即: (1) 生成項中的P直接由雷諾應(yīng)力計算得到; (2) 擴散項包含一個各向異性擴散系數(shù)(Pope 2000).該模型的不足之處在于其生成項建模方式,即當湍流處于不平衡或者平均流動的時間尺度量級小于等于湍流的特征時間尺度時,小尺度湍流可能沒有足夠時間調(diào)整以適應(yīng)大尺度上的變化(Al-Sharif 2011).標準耗散率模型方程在 LRR 模型 (Launder et al.1975) 和 SSG 模型(Speziale et al.1991) 中得到了應(yīng)用.

與對k-ε模型的改進工作相似地,研究者們也對雷諾應(yīng)力模型中的耗散率模型方程提出改進建議,如Launder (1996) 建議對方程系數(shù)進行修改以提高總體模擬效果,Hanjali′c 等(1997) 建議在方程中增加新的源項以提高對再附分離流的預(yù)測能力,等等.這其中針對低雷諾數(shù)效應(yīng)所開展的改進工作具有代表性.比如,流動在接近壁面時,湍流雷諾數(shù)逐漸降低至0,耗散率張量的各向異性也越來越明顯,因此有必要對湍流模型進行相應(yīng)的修正以適應(yīng)流動特征的變化.考慮低雷諾數(shù)效應(yīng)時,由耗散率的精確輸運方程(Daly 1970) 出發(fā)對模型進行修正,其具體形式如下

圖6

式中右側(cè)的前三項分別對應(yīng)于耗散率的生成、破壞和擴散,模擬方式與高雷諾數(shù)模型相同; 在高雷諾數(shù)模型中被忽略的黏性擴散保留其精確形式;為了加入第四項的影響,可以將參數(shù)Cε1和Cε2相應(yīng)地設(shè)置為Ret的函數(shù); 第五項的一種常見建模方式為

作為代表,在 Hajali′c-Jakirli′c低雷諾數(shù)模型 (Hanjali′c 和 Launder 1976) 和 Shima 低雷諾數(shù)模型 (Shima 1998) 中采用了針對低雷諾數(shù)效應(yīng)改進的ε輸運方程來封閉整個模型系統(tǒng).

在對k-ε模型的研究(Speziale et al.1992) 中發(fā)現(xiàn),采用湍動能耗散率作為特征尺度存在兩個問題:一是對于耗散率缺乏自然邊界條件;二是在壁面附近需要處理高階關(guān)聯(lián).這兩個問題對于模型的穩(wěn)定求解和壁面附近的模擬精度產(chǎn)生了不利的影響.針對這兩個問題,研究者們采用了不同的尺度來構(gòu)造尺度方程,比如均勻耗散率Chien 1982,Kalitzin et al.1996)、比耗散率ω(Menter 1994,Wilcox 1988)、湍流耗散時間τ(Kalitzin et al.1996,Kok & Spekreijse 2000) 等.當然,并不是所有在模型中出現(xiàn)的尺度都會被應(yīng)用到雷諾應(yīng)力模型中.圖6 給出了應(yīng)用于雷諾應(yīng)力模型的尺度變量發(fā)展脈絡(luò).

為了解決湍動能耗散率缺乏自然邊界條件的問題,部分研究者(Craft & Launder 1996,Craft1998) 采用了均勻耗散率來作為湍流特征尺度,其形式為

在兩方程k ?ω模型大獲成功的基礎(chǔ)上,比耗散率ω也被引入到雷諾應(yīng)力模型的封閉處理中,其形式為

典型代表是 Wilcox 的 Stress-ω模型 (Wilcox 2006) 和 SSG/LRR-ω模型 (Eisfeld & Brodersen 2005).后者參考了SST 模型的設(shè)計思路(Menter,1994): 在近壁區(qū)采用ω方程,在遠離壁面的區(qū)域采用ε方程,兩者通過混合函數(shù)連接(參考圖3),有效地解決了方程對來流湍流度敏感的問題以及方程在壁面附近數(shù)值剛性較大的問題.

盡管采用ω方程的SSTk-ω模型和 SSG/LRR-ω模型已經(jīng)取得了較為廣泛的應(yīng)用,但ω方程本身在理論上仍然存在一些問題.Togiti 和Eisfield (2015) 就指出: 因為ω在壁面附近存在奇異行為且缺乏自然邊界條件,基于ω方程的雷諾應(yīng)力模型的精度在某種程度上取決于近壁網(wǎng)格分辨率.為了解決這兩個問題,Togiti 和Eisfield 將兩方程模型框架中的尺度g(Speziale et al.1992) 引入SSG/LRR 模型,其具體形式為

數(shù)值模擬結(jié)果表明,采用g方程的雷諾應(yīng)力模型可以得到與ω方程相當?shù)慕Y(jié)果,并且其對于近壁網(wǎng)格分辨率的敏感性降低.

本文作者(Wang et al.2020)基于Menter 提出的k-kL模型思路,發(fā)展了渦黏尺度νt并用于構(gòu)建新型雷諾應(yīng)力模型.其具體形式為

該尺度不僅有自然的壁面邊界條件(壁面處為0),同時具備利用Boussinesq 假設(shè)和Bradshaw 近似獨立提供湍流應(yīng)力的能力.以此能夠?qū)崿F(xiàn)傳統(tǒng)渦黏模型和雷諾應(yīng)力模型的結(jié)合,在計算穩(wěn)定性和效率方面具備一定的優(yōu)勢.

利用雷諾應(yīng)力模型模擬轉(zhuǎn)捩問題,其實質(zhì)也可歸為湍流尺度問題.目前針對渦黏模型的轉(zhuǎn)捩模型以Langtry 和Menter (2009) 提出的γ-Reθ模型為代表,通過在已有的湍流模型中附加間歇因子輸運方程,達到控制湍流生成和耗散的目的.Nie 等(2017)成功將γ-Reθ模型擴展至SSG/LRR-ω模型,給雷諾應(yīng)力模型求解轉(zhuǎn)捩問題提供了新思路.

4 高階矩模型數(shù)值求解

高階矩湍流模型相比線性渦黏模型在數(shù)值求解時也更為困難.當采用線性渦流黏性模型時,線性雷諾應(yīng)力張量與(平均流) 黏性應(yīng)力張量成比例,有助于NS 方程解算器的穩(wěn)定; 相反地,當使用高階矩模型時,非物理的雷諾應(yīng)力張量很容易引起失穩(wěn)效應(yīng).因此,數(shù)值求解高階矩模型方程時要保證獲得的雷諾應(yīng)力張量合理(即可實現(xiàn)性),這對相關(guān)物理限制、數(shù)值離散方法、邊界條件和網(wǎng)格因素等提出了更高的挑戰(zhàn).圖7 給出了影響高階矩模型數(shù)值求解的主要因素,下面分別進行討論.

圖7

4.1 可實現(xiàn)性限制

Schumann (1977) 首先提出了可實現(xiàn)性這一概念.他要求雷諾應(yīng)力模型得到的能量分量應(yīng)為正,且非對角分量滿足Schwarz 不等式.即

其后,Lumley (1978,1983) 指出在雷諾應(yīng)力模型設(shè)計中應(yīng)系統(tǒng)地使用可實現(xiàn)性限制.然而,隨著實際應(yīng)用逐漸增多,人們發(fā)現(xiàn)即使在模型的設(shè)計過程中滿足可實現(xiàn)性,在數(shù)值求解時也可能產(chǎn)生違反物理的雷諾應(yīng)力結(jié)果.

為解決該問題,Chassaing (2003) 提出了一個顯式的強制可實現(xiàn)方法: 在每個計算網(wǎng)格點上檢查可實現(xiàn)性約束,并且在違反的情況下,雷諾應(yīng)力模型的局部解被強制設(shè)為零.而后,Chaouat(2006) 改進了這一措施,對正應(yīng)力分量強制實施同性湍流,切應(yīng)力分量仍強制為零.如下

該限制措施由于簡單高效在復(fù)雜問題中獲得了廣泛的應(yīng)用.如作者(Wang et al.2018) 基于該可實現(xiàn)性限制,采用三階精度WCNS 格式對SSG/LRR-ω模型進行了離散求解,在跨聲速渦輪葉片、翼身組合構(gòu)型等真實外形求解中仍能保證良好的數(shù)值穩(wěn)定性.

Mor-Yossef(2014)將最初針對兩方程模型開發(fā)的無條件正收斂隱式格式(UPC)拓展至雷諾應(yīng)力模型.UPC 方法具有任意時間步正應(yīng)力保正性的特點,但其他可實現(xiàn)性約束仍有可能會違反.因此,針對式33(b) 定義了和滿足的臨近點平均值并施加了如下限制條件

4.2 數(shù)值離散方法

對雷諾應(yīng)力模型的數(shù)值離散包括空間方面和時間方面.空間離散方面最受關(guān)注的是對流項的處理.目前,基于近似黎曼解的數(shù)值方法應(yīng)用最為廣泛,但在NS 方程中表現(xiàn)較好的低耗散通量格式可能會促進數(shù)值不穩(wěn)定(Brun et al.2000,Gerolymos & Vallet 2007).Naser 等(2014) 指出,數(shù)值不穩(wěn)定性應(yīng)歸因于對接觸間斷的錯誤處理.對于雷諾應(yīng)力模型的對流通量由更具耗散性的質(zhì)量通量來近似更為合適.而NS 方程仍可采用低耗散通量.另一方面,為提升黎曼求解的精度,插值過程需要更多的模板點.而在不光滑區(qū)需要在插值模板中增加限制器(包括加權(quán)等方式).由于雷諾應(yīng)力各分量在無量綱后相比平均流動量具有量級小、變化大的特點,在插值過程中要格外注意經(jīng)驗參數(shù)的選擇(Wang et al.2018).

時間求解方面,傳統(tǒng)的顯式格式雖然魯棒,但在真實流動問題中很難突破時間步長的限制.Chassaing (2003) 提出了一種高效、魯棒的隱式格式.該格式依賴于雙時間步方法,用于定常流動時雙時間步是人工時間增量(LDT).同時結(jié)合可實現(xiàn)性限制,實現(xiàn)了二階精度雷諾應(yīng)力模型的穩(wěn)定求解.在隱式時間過程中,如何處理源項對收斂性和穩(wěn)定性有較大的影響.正如Chaouat(2011)提出的建議: 將時間推進格式分裂,其中控制方程的對流項、擴散項在時間上顯式處理,而雷諾應(yīng)力方程的源項隱式處理.同時,他又針對性的設(shè)計了隱式格式來保持湍流耗散和法向雷諾應(yīng)力的正性.Mor-Yossef (2014) 發(fā)展的UPC 格式同樣是基于對源項的隱式處理,并在推廣到雷諾應(yīng)力模型時獲得了很好的收斂性表現(xiàn).

4.3 邊界條件

邊界條件是影響模型穩(wěn)定性的重要因素,尤其是壁面邊界條件.對于雷諾應(yīng)力各分量,在壁面處應(yīng)保證為0.但尺度方程由于變量形式不同,其邊界條件處理方式和數(shù)值特性有明顯差異.如第3 節(jié)所述,早期的RSM 通常基于耗散率尺度ε,但由于數(shù)值剛性相關(guān)問題,基于ε的RSM 很難用于復(fù)雜外形模擬.目前在真實工程模擬中,尤其航空領(lǐng)域,最常用的尺度提供變量是ω.特別是Eisfeld等(2005)引用了Menter 的混合思想提出的SSG/LRR-ω雷諾應(yīng)力模型,在許多實際的航空流動中,如三維跨音速流動中的激波– 邊界層相互作用、高升力下的自由渦和拐角流動中都得到了證明.然而,耗散率ω在壁面邊界為無窮大值.Wilcox 及Menter 等采用的近似條件ωwall=60μ/(βy21)(非自然邊界條件),在一定程度上會影響數(shù)值穩(wěn)定性和靈敏度.這點在高精度離散時尤為明顯.針對該問題,Bassi 等(2011) 提出了ω對數(shù)形式轉(zhuǎn)換,有效避免了非物理狀態(tài); Hartmann 等(2011) 提出了新的ω邊界條件,并在k方程中添加人工黏性.他們均實現(xiàn)了k-ω模型的高精度離散求解.而針對雷諾應(yīng)力模型,作者(Wang et al.2020) 發(fā)展的νt尺度具備自然的邊界條件(壁面處為0),基于經(jīng)典的LU-SGS 方法和可實現(xiàn)性限制成功地實現(xiàn)了雷諾應(yīng)力模型5-7 階的高精度離散和穩(wěn)定求解,其中源項均采用了Chaouat 建議的隱式處理.圖8 展示了不同離散精度下翼型尾跡模擬結(jié)果對比.對RANS 方程和RSM 方程采用一致7 階精度離散求解,在該粗糙網(wǎng)格上就能與實驗吻合較好.

圖8

4.4 網(wǎng)格因素

如同數(shù)值方法,雷諾應(yīng)力模型對網(wǎng)格方面的要求也更為嚴苛.雖然目前尚未發(fā)現(xiàn)雷諾應(yīng)力模型存在針對網(wǎng)格的特殊要求,但下面幾方面是要格外注意的.首先,基于結(jié)構(gòu)網(wǎng)格的有限差分格式是最早試圖求解雷諾應(yīng)力方程的方法(Wilcox 2006).但近年來,隨著離散格式精度的不斷提高,網(wǎng)格導(dǎo)數(shù)的幾何守恒性對離散準確性,尤其是穩(wěn)定性的影響也更加顯著.針對NS 方程,Visbal 和Gaitonde (2002) 發(fā)現(xiàn)通過采用守恒形式的網(wǎng)格導(dǎo)數(shù)可大幅減小該誤差.而多數(shù)CFD 軟件中,雷諾應(yīng)力模型往往采用與NS 方程相似的求解過程.因此減小雷諾應(yīng)力模型在網(wǎng)格轉(zhuǎn)換過程中的誤差也是十分必要的.本文作者(Wang et al.2018) 將Deng 等(2011,2013) 提出的對稱守恒網(wǎng)格導(dǎo)數(shù)算法(SCMM) 擴展至雷諾應(yīng)力模型,并獲得了良好的效果.

非結(jié)構(gòu)網(wǎng)格相比結(jié)構(gòu)網(wǎng)格更適用于工程問題,但由于高階矩模型方程數(shù)較多,無論基于有限體積方法還是有限元方法都需要巨大的存儲和計算花費.Mor-Yossef (2016) 成功地將其基于結(jié)構(gòu)網(wǎng)格上(有限差分法) 提出UPC 格式擴展至非結(jié)構(gòu)網(wǎng)格(有限體積法),但并未像基于有限差分法時開展三階精度的嘗試.主要原因是有限體積方法需要更復(fù)雜的模板點.采用DG 為代表的有限元方法可以克服模板向外擴展的困難,但又會引入新的問題,即壁面距離的準確計算.一種方式是通過求解程函方程(Eikonal-equation) 得到任意點到最近壁面邊界的距離.Liu 等(2010) 提出了應(yīng)用與內(nèi)部流動問題的壁面距離算法.Schoenawa 和Hartmann (2014) 將該算法擴展至外流問題,并成功用于SST 模型的高精度求解.但在面對復(fù)雜外形時,還需添加人工黏性來保證穩(wěn)定.而在非結(jié)構(gòu)網(wǎng)格上,高于三階的雷諾應(yīng)力模型離散還少有研究.

圖9

5 高階矩模型對比線性渦黏模型

Boussinesq 假設(shè)是將湍流脈動運動類比分子熱運動產(chǎn)生黏性的過程.對于剪切應(yīng)變占主導(dǎo)的流動具有良好的適用性,這也是工程中最關(guān)心的邊界層等為代表的流動類型.但隨著CFD 技術(shù)求解問題的日漸復(fù)雜,基于 Boussinesq 假設(shè)的線性渦黏模型在眾多流動問題中暴露出了不足.Wilcox(2006)總結(jié)了一些線性渦黏模型失敗的并且值得注意的應(yīng)用類型,包括: 平均應(yīng)變率突然變化的流動、曲壁面上的流動、旋轉(zhuǎn)流動(或渦流動)、三維流動等.

5.1 平均應(yīng)變率突然變化的流動

考慮Tucher 和Reynolds (1968) 實施的實驗.該實驗中平均應(yīng)變率維持在一個有限距離中,然后超過臨界點后突然去除.湍流由均勻拉伸而變成各向異性,并在停止拉伸點(臨界點) 的下游逐漸趨向各向同性.Wilcox 和Rubesin (1980) 采用他們提出的k-ω2模型對該流動進行了模擬,結(jié)果如圖9 所示.當平均應(yīng)變率突然在x ≈2.3 m 去除后,線性渦黏模型預(yù)測了一個各向同性的瞬時,即所有法向的雷諾應(yīng)力變得平等.進而,湍流以有限速率接近各向同性.

更值得注意的是圖9 中K值的變化,定義如下

如果x ≈2.3 m 的下游沒有消除平均應(yīng)變率,線性渦黏模型預(yù)測的K值應(yīng)為連續(xù)的.但前面提到,由于Boussinesq 假設(shè)是對平均運動黏性的比擬,因此當平均應(yīng)變率去除后,K值會產(chǎn)生突變.同時由于Boussinesq 假設(shè)對正應(yīng)力的估計是各向同性的,因此K值在突變后將為0.后面,Wilcox(2006) 利用其發(fā)展的ω基雷諾應(yīng)力模型對該問題又進行了模擬,所得結(jié)果很好地克服了線性渦黏模型的缺陷.

5.2 曲壁面上的流動

這里考慮一個經(jīng)典的彎曲管道流動,由Smits 等(1979) 進行的實驗.實驗使用了一個高度為0.127 m 的恒截面積方形風管,中部具有30?彎曲.CFD 模擬的主要焦點是評估凸壁曲率的湍流模型(圖10a 中為下壁),被NASA 湍流模型網(wǎng)站(Rumsey 2020) 選用并收錄了多個模型的計算結(jié)果.

本文引用該網(wǎng)站上SA,SA-RC 和 SSG/LRR-ω模型的結(jié)果進行對比分析.觀察壁面摩阻,SA模型在凸壁面附近的預(yù)測明顯高于實驗值.真實的物理過程中,凸壁面附近的湍流水平往往低于直壁面.而SA 模型缺失了捕捉該現(xiàn)象的能力.同時對于凸壁面附近的雷諾應(yīng)力,SA 模型也給出了較大的偏差.對于基于k的兩方程模型,這種情況也依然存在.通過引入判斷函數(shù)在曲面邊界附近修正渦黏模型方程(稱為旋轉(zhuǎn)/曲線修正,如SA-RC (Shur et al.2000)) 是一種常用的工程方法,然而它們是不能用于任意流動的特殊修正方法(Wilcox 2006).馮·卡門曾提出關(guān)于曲壁面上流動的經(jīng)典穩(wěn)定性論點,認為k方程應(yīng)該是的方程更合適.而一些 RC 修正也正是基于概論點構(gòu)建的,如Wilcox 和Chambers (1977) 修正等.從這個角度講,雷諾應(yīng)力模型自然具備預(yù)測流線曲率影響的能力.圖10 中SSG/LRR-ω模型的結(jié)果也印證了這點.對于凹壁面(圖10a 中為上壁),由于存在Gortler 渦導(dǎo)致所有二維模擬都較差.但雷諾應(yīng)力模型往往仍好于線性渦黏模型.

5.3 旋轉(zhuǎn)/渦流動

考慮一個航空工業(yè)中最典型的翼梢渦流動,實驗在 NASA 蘭利研究中心的低速風洞中開展(Chow 1993,1997).實驗?zāi)Ｐ褪且粋€半展長機翼(翼型為NACA0012),翼根固定于風洞壁面,安裝角 10?(即攻角 10?).如圖11(a)所示,RANS 模擬的主要難點在于翼梢渦的保持,即過大的物理黏性會導(dǎo)致翼梢渦在向下游的過程中被逐漸耗散.圖11(b)中給出了DLR 采用TAU 軟件開展的CFD 模擬 (Cecora et al.2015),RANS 模型包括原始版 SA,SST,SSG/LRR-ω以及 JHh-v2 模型,后兩者是雷諾應(yīng)力模型.橫坐標表示以翼梢后緣為原點的x方向坐標,縱坐標表示x對應(yīng)截面上渦核中心流向速度.可以看到,渦黏模型由于過高估計渦核處的渦黏系數(shù),導(dǎo)致翼梢渦迅速衰減.而兩個雷諾應(yīng)力模型無此現(xiàn)象.

在非定常渦流動中,基于高階矩模型的混合RANS/LES 方法依然具備潛力.圖12 給出了SSTURANS,SST-IDDES,SSG/LRR-URANS 和SSG/LRR-IDDES 方法對24.6?迎角鈍前緣三角翼的模擬情況(Wang et al.2017).實驗在NASA Langley 中心的風洞實施(Chu & Luckring 1996),來流馬赫數(shù)為0.85,基于三角翼氣動弦長的雷諾數(shù)達6 百萬.實驗迎角增加至24.6?時,主渦由于激波誘導(dǎo)作用而發(fā)生破裂,破裂位置位于翼根弦長60%～80%間.SST-URANS 得到的展向壓力分布,始終存在兩個峰值(主渦和二次渦),說明其推遲了渦破裂的發(fā)生.SSG/LRR-URANS 在80%占位處峰值消失,說明捕捉到了渦破裂.但對破裂后非定常流動的模擬較差.而添加IDDES 修正可彌補這點,即SSG/LRR-IDDES 模型.而對比傳統(tǒng)的SST-IDDES,SSG/LRR-IDDES 的結(jié)果更接近實驗值.

圖10

圖11

5.4 三維流動

仍然考慮一個航空工業(yè)中的典型流動,翼身結(jié)合處流動(典型的三維T 型區(qū),如圖13a).該問題來源于第五屆AIAA 阻力預(yù)測會議(DPW-5),計算模型為波音公司設(shè)計的通用(寬體客機)研究模型CRM (Vassberg et al.2008).該飛機模型的“設(shè)計巡航速度” 為Ma= 0.85,在美國NASA 的國家跨聲速風洞以及埃姆斯研究中心風洞中均進行了詳細的風洞實驗.實驗中在小迎角下(4?附近),翼身結(jié)合處并未發(fā)現(xiàn)分離渦.而在DPW-5 會議提交的CFD 模擬中,渦黏模型普遍給出了分離渦,并且導(dǎo)致了“升力突降” 的現(xiàn)象(Levy et al.2013).隨后,該問題被廣泛關(guān)注,包括NASA 專門組織了拐角流動系列實驗(Rumsey 2016)等等.圖13 為德國宇航院Togiti 等(2014)對該問題開展的研究,采用SA 和SST 等渦黏模型在翼身連接處出現(xiàn)分離渦.主要原因是Boussinesq 假設(shè)無法考慮雷諾正應(yīng)力的各向異性,即在三維流動中存在潛在的缺陷.而雷諾應(yīng)力模型自然無此問題.參考雷諾應(yīng)力模型,Spalart 等(2000)提出了對Boussinesq 假設(shè)的修正,即二次本構(gòu)關(guān)系(QCR),很好地解決了該問題.但如同RC 修正,這是一種不通用的方式.

本節(jié)通過數(shù)個典型流動案例說明了高階矩湍流模型相比于傳統(tǒng)渦黏模型存在的幾方面優(yōu)勢.然而,優(yōu)勢的獲得往往伴隨著一定的代價.最主要的是計算花費和魯棒性要求的提升.對于前者,多位學(xué)者(Mor-Yossef 2014,C′ecora et al.2015,Wang et al.2020)針對不同的算例開展過評估.雖然在相同網(wǎng)格上,高階矩模型需要高出渦黏模模型30%～60%的計算花費.但考慮到相同網(wǎng)格時高階矩模型獲得更準確結(jié)果的潛力,尤其針對渦流動.因此在達到同等準度條件下,其對網(wǎng)格量的需求往往不會高于渦黏模型.而對于魯棒性,則仍需要在模式理論、數(shù)值算法等多個方面繼續(xù)開展研究.這也是目前影響高階矩模型廣泛應(yīng)用最大的挑戰(zhàn).

圖12

6 CFD 軟件中的高階矩湍流模型

CFD 軟件作為湍流模型應(yīng)用于實際流動的載體,對于一些模型的發(fā)展起到了巨大的推動作用.例如由Menter(1994)提出的SST 湍流模型,以及其衍生的γ-Reθ轉(zhuǎn)捩模型(Menter et al.2006)、尺度自適應(yīng) (SAS) 模型(Menteret al.2003) 等,借助ANSYS 公司旗下的CFX,Fluent 等 CFD 軟件在世界范圍內(nèi)得到了廣泛的應(yīng)用.而上述模型之所以被納入該軟件,除了Menter 本人作為ANSYS公司技術(shù)顧問外,模型本身也必須具備魯棒性強、適應(yīng)范圍廣等優(yōu)秀的性能.在CFD 軟件發(fā)展的早期,集成的湍流模型基本上為低階矩模型.直到1990 年代后,隨著計算機硬件水平的提升和二階矩模式理論在一定程度上的成熟,其中一些表現(xiàn)較好的模型才開始納入商業(yè)CFD 軟件(Launder 1989).

圖13

在ANSYS Fluent 中(以2019 R3 版本為例),二階矩湍流模型被分為兩大類: 基于ε的模型和基于ω的模型.前者默認采用Fu 等(1987)在LRR-IP 模型基礎(chǔ)上改進的線性壓力–應(yīng)變模型(也稱 Fu-Launder-Leschziner 模型),還包括二次壓力 – 應(yīng)變模型 (即 SSG 模型).后者則包括 Stress-ω模型和Stress-BSL 模型.這兩種模型均基于LRR-QI 壓力– 應(yīng)變模型,區(qū)別在于尺度方程ω分別采用的Wilcox (2006) 版本和Menter (1994) 版本.在軟件中對模型的可選項還包括邊界條件、壁面反射效應(yīng)和近壁處理等內(nèi)容.同時針對特殊問題,可添加SAS 修正(分離渦問題)、低雷諾數(shù)數(shù)修正或剪切流修正等.

ANSYS CFX 與Fluent 相似,雷諾應(yīng)力模型也被分為了兩大類,即基于ε的模型和基于ω的模型.前者包括了LRR-IP 模型、LRR-QI 模型和SSG 模型,后者則包括了Stress-ω模型和Stress-BSL模型.值得注意的是,在該軟件中對所有模型均采用了先進壁面函數(shù).

NUMECA FINE/Marine 中也集成了雷諾應(yīng)力模型,即Rij-ω模型.該模型是在 Shima (1998)的低雷諾數(shù)模型基礎(chǔ)上將ε方程替換為ω方程(Deng & Visonneau,1997) 并采用了新的近壁模型而得到的.而壓力與應(yīng)變關(guān)聯(lián)項是在LRR-QI 模型基礎(chǔ)上的改進版.

除了商業(yè)CFD 軟件,專用CFD 軟件也為一些二階矩模型在該領(lǐng)域的發(fā)展和推廣起到了積極作用.例如由美國NASA 主導(dǎo)開發(fā)的CFL3D,FUN3D 軟件、德國DLR 開發(fā)的TAU 軟件均集成了SSG/LRR-ω模型.通過系統(tǒng)的驗證、確認及應(yīng)用,使該模型成為針對航空航天氣動問題中高階矩模型的標桿(Eisfeld et al.2016).表2 給出了部分CFD 軟件中高階矩模型集成情況.

表2 部分CFD 軟件中高階矩模型情況

7 高階矩湍流模型挑戰(zhàn)與展望

高階矩模型是湍流模式理論發(fā)展中舉足輕重的部分,對其物理機制和數(shù)值特性的研究不僅是認識湍流的基礎(chǔ)科學(xué)問題,更是解決復(fù)雜湍流的工程需求.從上述第2～6 節(jié)的進展情況看,高階矩模型未來面臨的最大挑戰(zhàn)是如何確保在不低于渦黏模型準確度的同時,實現(xiàn)魯棒、高效的求解,進而得到更廣泛的應(yīng)用.這需要在湍流建模理論、尺度提供方程、數(shù)值離散方法、CFD 軟件集成及驗證與確認等方面持續(xù)開展工作,如圖14.

(1)應(yīng)力輸運方程中各源項(尤其壓力與應(yīng)變關(guān)聯(lián)項)的物理機制及建模理論仍然將是高階矩模型研究的核心.例如,1980 年代末開始興起的重整化速度群(Renormalization group)理論.Yakhot&Orszag(1986)首先將其用來分析湍流并成功地建立了基于該理論的k-ε模型(RNGk-ε模型).該理論基于統(tǒng)計學(xué)方法,通過將自然空間中湍流方程轉(zhuǎn)化為傅里葉空間(譜空間) 進行分析,得到獨立于經(jīng)驗標定的湍流方程參數(shù).其后,Rubinstein & Barton(1990) 將重整化速度群理論用于發(fā)展二階矩模型,取得了令人鼓舞的效果,表明了該理論在推導(dǎo)高階矩模型中具有很好的潛力.系統(tǒng)性逆尺度(Systematic upscaling) 是一種新的多尺度計算方法,用于精確推導(dǎo)在越來越大尺度上給定的物理系統(tǒng)的控制方程(Brandt 2005).該方法結(jié)合了過去出現(xiàn)的兩種多層計算模式的優(yōu)勢,即應(yīng)用數(shù)學(xué)中的多重網(wǎng)格和理論物理中的重整化速度群.美國NASA 在CFD 的2030 愿景報告(Slotnick et al.2014) 中將其列為有希望推動湍流模式理論變革的方法之一.機器學(xué)習(xí)(Machine learning) 是通過構(gòu)建算法從數(shù)據(jù)中建立模型.近年來,得益于人工智能發(fā)展的高潮,機器學(xué)習(xí)也被成功用于湍流建模,并迅速成為本領(lǐng)域的熱門(Lav et al.2019,張偉偉等2019).Tracey 等(2015) 通過神經(jīng)網(wǎng)絡(luò)模型構(gòu)建了SA 方程中新的源項,而模型訓(xùn)練所采用的數(shù)據(jù)可通過實驗結(jié)果或高分辨數(shù)據(jù)反演得到.相似的思路也可用于k方程模型(Zhang 2015) 甚至高階矩模型.雖然壓力與應(yīng)變關(guān)聯(lián)作用很難通過實驗直接測定,但基于DNS 數(shù)據(jù)對壓力與應(yīng)變關(guān)聯(lián)項進行機器學(xué)習(xí)建模仍有相當大的潛力.

圖14

(2) 尺度提供方程一直是RANS 理論中的關(guān)鍵.從圖6 的發(fā)展趨勢來看,不同的尺度方程通過互補、融合可以實現(xiàn)更優(yōu)秀的特性,如更符合物理的邊界條件、更小的數(shù)值剛性制約以及更好的分離渦求解能力等等.Wang 等(2020) 推導(dǎo)的νt尺度方程,通過引入Menter 的k-kL模型近壁耗散修正,提升了雷諾應(yīng)力模型的近壁漸進特性.Eisfeld 和Rumsey (2020) 將Yap (1987)在ε方程采用的長度尺度修正引入到了ω方程,使新的SSG/LRR-ω模型能更好地模擬后體分離問題.

(3) 數(shù)值離散方法是實現(xiàn)高階矩模型廣泛應(yīng)用的關(guān)鍵.從軟件集成的情況來看,二階精度下的雷諾應(yīng)力方程離散已比較成熟,無論針對結(jié)構(gòu)網(wǎng)格采用有限差分離散還是非結(jié)構(gòu)網(wǎng)格采用有限體積離散.但上升到高精度離散時,穩(wěn)定性和收斂性則遠不如低階矩模型.而高階矩模型的優(yōu)勢,又恰好是針對旋轉(zhuǎn)、分離這類復(fù)雜的流動,需要高精度格式的低耗散特性作支撐.因此,未來的關(guān)鍵和難點是如何實現(xiàn)高階矩模型的高階精度離散及穩(wěn)定求解.

(4) CFD 軟件能夠為湍流模型的應(yīng)用與推廣提供平臺.但從表2 總結(jié)的情況看,國外開發(fā)的商業(yè)軟件極少會采用我國原創(chuàng)的湍流模式.未來對于湍流模式的推廣,一種是基于開源CFD 軟件平臺(如Openfoam 等),但主要適合于學(xué)術(shù)研究和交流; 一種是走自主發(fā)展的道路.目前,諸多科研院所自主開發(fā)了CFD 軟件或代碼,同時集成了各自發(fā)展的湍流模式.但由于軟件規(guī)模與國外商業(yè)CFD 軟件差距懸殊,影響力仍十分有限.只有從國家層面加以推動整合,才能形成有競爭力的軟件產(chǎn)品,為我國原創(chuàng)的湍流模式提供良好的應(yīng)用平臺.

(5) 無論是開展新型高階矩模型研究,還是將已有模型納入CFD 軟件作為產(chǎn)品模塊,都缺少不了系統(tǒng)的驗證與確認.例如,SSG/LRR-ω模型通過參與美國NASA 組織的湍流模型資源(TMR)網(wǎng)站(Rumsey 2020)、第五次阻力預(yù)測會議(DPW5) (Togiti et al.2014)、第二次高升力預(yù)測會議(HiLiftPW-2)(Thompson & Hassan 2015) 等的驗證與確認,得到了廣泛的好評和迅速的推廣.我國在近些年也加大了這方面的投入,如空氣動力研究與發(fā)展中心組織的第一屆航空CFD 可信度研討會 (AeCW-1) (王運濤等2019),已經(jīng)有高階矩模型參與對比(聶勝陽等2019).且標準模型和實驗數(shù)據(jù)均為自主設(shè)計和測試.因此,未來為支撐原創(chuàng)湍流模式和自主CFD 軟件的發(fā)展,驗證與確認工作必須走在前面.

8 結(jié)論

本文從CFD 由RANS 到LES 的發(fā)展趨勢談起,對高階矩湍流模型進行了綜述.重點總結(jié)了高階矩模型各項的建模過程、尺度提供方程的演化歷程和數(shù)值求解技術(shù)的研究需求.同時,通過典型湍流問題展示了其相對于傳統(tǒng)渦黏模型的優(yōu)勢,并且給出了部分CFD 軟件對高階矩模型的集成情況.最后對高階矩湍流模型未來的研究與發(fā)展方向進行了展望.本文得到的主要結(jié)論如下:

(1) 進入新世紀后,隨著計算機硬件水平和數(shù)值求解技術(shù)的發(fā)展,高階矩模型又再次迎來了發(fā)展的高潮.多個國家的科研機構(gòu)均提出過復(fù)興高階矩模型的計劃.然而,在工業(yè)領(lǐng)域,高階矩模型仍未被廣泛推廣,因此需要進一步加強應(yīng)用性研究.

(2) 應(yīng)力輸運方程中各源項(尤其壓力與應(yīng)變關(guān)聯(lián)項) 的建模理論是高階矩模型研究的核心.近年來隨著DNS 的不斷突破,高分辨實驗觀測技術(shù)的發(fā)展給方程建模提供了更準確、豐富的信息.而諸如重整化速度群、系統(tǒng)性逆尺度、機器學(xué)習(xí)等理論的發(fā)展,有望為高階矩建模提供新思路.

(3) 尺度項/方程對高階矩模型的實際表現(xiàn)有著綜合性的影響.根據(jù)歷史上尺度方程的演化規(guī)律可以得出其研究重點包括更符合物理的邊界條件、更小的數(shù)值剛性制約、更好的分離渦求解能力以及對轉(zhuǎn)捩模擬的能力等.

(4)當前專門針對高階矩模型的數(shù)值方法研究還很少.從學(xué)術(shù)的角度,高階矩模型還很難做到像渦黏模型一樣,與NS 方程保持一致的任意階精度求解.這將導(dǎo)致在湍流量的計算上產(chǎn)生過多誤差,甚至有時比渦黏模型的效果差.從工程的角度,高階矩模型仍然缺乏魯棒性.這是制約高階矩模型廣泛應(yīng)用的關(guān)鍵問題.

(5) 高階矩模型相比渦黏模型的優(yōu)勢仍然不夠清晰,這需要持續(xù)的開展驗證與確認研究.而隨著湍流模擬方法經(jīng)歷由RANS 到LES 的轉(zhuǎn)變,LES 框架下Boussinesq 假設(shè)是否更普適性以及建立輸運方程模型是否有優(yōu)勢等問題也需要開展探索.

(6) 高階矩模型的推廣需要優(yōu)秀的CFD 軟件平臺,而國外開發(fā)的商業(yè)軟件極少會采用我國原創(chuàng)的湍流模式.對于我國而言,只有堅持走自主創(chuàng)新之路,形成有競爭力的軟件產(chǎn)品,才能真正為原創(chuàng)的湍流模式提供良好的應(yīng)用平臺.

致 謝國家自然科學(xué)基金(12002379)、湖南省自然科學(xué)基金(2020JJ5648)、國防科技大學(xué)科研計劃(ZK20-43) 和國家專項工程(GJXM92579) 資助項目.