基于基因表達式編程的數(shù)據(jù)驅(qū)動湍流建模

趙耀民 *, 徐曉偉

* (北京大學(xué)應(yīng)用物理與技術(shù)研究中心,高能量密度物理數(shù)值模擬教育部重點實驗室,北京 100871)

? (北京大學(xué)工學(xué)院力學(xué)與工程科學(xué)系,湍流與復(fù)雜系統(tǒng)國家重點實驗室,北京 100871)

** (澳大利亞墨爾本大學(xué)機械工程系,澳大利亞墨爾本 3010)

引言

湍流廣泛存在于自然界和工程應(yīng)用中,準確認識和預(yù)測湍流現(xiàn)象十分重要.湍流的研究有近百年的歷史,人們由理論、實驗、和數(shù)值模擬方向出發(fā)取得了一系列的研究成果[1].然而,在以航空航天為代表的工程應(yīng)用中,流動幾何邊界往往較為復(fù)雜,且流動雷諾數(shù)較高,研究的難度很大,相應(yīng)流動的基本規(guī)律和預(yù)測模型仍有待進一步深入研究.因此,湍流的機理、建模及控制仍然是航空航天等領(lǐng)域的重要瓶頸之一,針對湍流的研究有著重要的基礎(chǔ)和應(yīng)用價值.

計算流體動力學(xué)(computational fluid dynamics,CFD)是湍流的重要研究手段,常見的數(shù)值模擬方法可大致分為直接數(shù)值模擬(direct numerical simulation,DNS)、大渦模擬(large eddy simulation,LES)和雷諾平均模擬(Reynolds-averaged Navier-Stokes,RANS)等.其中,直接數(shù)值模擬和大渦模擬均需要較高的網(wǎng)格分辨精度.對于航空航天等工程應(yīng)用中的高雷諾數(shù)流動,直接數(shù)值模擬和大渦模擬所需要的網(wǎng)格量巨大,遠超當前超級計算機的算力極限[2].相對而言,RANS 雖然精度不及前兩種方法,但計算效率很高,是各種工程應(yīng)用中的主要研究手段.

雷諾平均模擬的結(jié)果往往受限于湍流模型的預(yù)測精度.以雷諾應(yīng)力輸運模型[3]為代表的高階矩封閉模型雖然精度較高,但計算量較大且收斂性較差,其應(yīng)用范圍較為有限.實際工程中的常用湍流模型包括 SA 模型[4]、k-ε模型[5]、k-ω模型[6]、SST 模型[7]等.這些模型均基于Boussinesq 渦黏假設(shè),即假設(shè)雷諾應(yīng)力的各向異性部分與應(yīng)變率張量存在線性關(guān)系.然而,當流動較為復(fù)雜,特別是出現(xiàn)如邊界層分離、強烈混合等現(xiàn)象時,雷諾應(yīng)力與應(yīng)變之間往往不再滿足簡單的線性關(guān)系,因此基于渦黏假設(shè)的模型將會出現(xiàn)較大誤差[8-9].

近年來,隨著相關(guān)算法的快速發(fā)展和可用數(shù)據(jù)的迅速增長,機器學(xué)習(xí)方法開始在流體力學(xué)領(lǐng)域得到廣泛關(guān)注[10],其中一個重要的研究方向是利用機器學(xué)習(xí)方法開發(fā)湍流模型[11-12].如圖1 所示,可以將數(shù)據(jù)驅(qū)動的湍流建模大致分為3 個步驟,即數(shù)據(jù)預(yù)處理、模型訓(xùn)練以及模型驗證和預(yù)測.在數(shù)值模擬或者實驗獲得的高可信度(high-fidelity,Hi-Fi)流場基礎(chǔ)上,首先通過數(shù)據(jù)預(yù)處理提取輸入特征變量(X表示)及輸出變量(Y表示);再由輸入及輸出變量出發(fā),采用不同的機器學(xué)習(xí)算法及策略訓(xùn)練模型;最后一般要對訓(xùn)練得到的模型(Y=f(X))進行一定的測試,驗證其預(yù)測能力.需要注意的是,最終模型的預(yù)測能力依賴于訓(xùn)練數(shù)據(jù)、特征變量的選擇、算法及策略等多種因素.

圖1 機器學(xué)習(xí)湍流建模示意圖Fig.1 Schematic for turbulence modelling with machine learning methods

根據(jù)輸出變量的不同,現(xiàn)有的數(shù)據(jù)驅(qū)動湍流建模研究大致可分為兩個類型.其中一部分研究關(guān)注于發(fā)展修正模型,即采用不同的機器學(xué)習(xí)方法修正現(xiàn)有的模型及其參數(shù),以改善原有模型的預(yù)測精度.Edeling 等[13]基于貝葉斯方法估計壓力梯度邊界層中模型誤差并嘗試修正模型參數(shù).Parish 和Duraisamy[14]提出了流場反演與機器學(xué)習(xí)(field inversion and machine learning,FIML)方法.此框架首先通過統(tǒng)計方法反演得到修正參數(shù)場作為訓(xùn)練目標,例如湍動能方程中湍流生成項的修正系數(shù)的空間分布;進一步地,采用機器學(xué)習(xí)方法重構(gòu)這一修正場,建立輸入特征量與修正場之間的映射關(guān)系.FIML 框架的優(yōu)勢在于較為靈活,其中流場反演部分可以使用貝葉斯方法[15],也可基于伴隨方法求解[16];另外,機器學(xué)習(xí)部分也可以由高斯過程、神經(jīng)網(wǎng)絡(luò)等不同方法實現(xiàn)[17-18].值得注意的是,這種框架中反演得到的修正場并不總是“可訓(xùn)練的”[19],而訓(xùn)練過程的缺陷會導(dǎo)致反演場的誤差,進而導(dǎo)致最終預(yù)測精度和泛化能力的不足.一種可能的改進方向是將反演和學(xué)習(xí)兩個過程耦合,但是相關(guān)的研究還較為初步[19].

另外一個大類的工作可歸類為替代模型[12]的開發(fā).研究者們不滿足于單純修正現(xiàn)有的模型參數(shù),而是探索構(gòu)建滿足一定物理約束的RANS 模型以替代原有模型.針對RANS 和高精度流場的雷諾應(yīng)力偏差,Wang 等[20]采用隨機森林方法構(gòu)建模型,在周期山測試算例中可以明顯改善雷諾應(yīng)力場的預(yù)測精度,但是模型在后驗測試中對于平均流場的預(yù)測精度還需要進一步驗證.Ling 等[21]采用深度神經(jīng)網(wǎng)絡(luò)方法訓(xùn)練滿足伽利略不變性的雷諾應(yīng)力模型.訓(xùn)練基于速度梯度基張量和不變量[22],得到的模型在后驗測試中可以改善對于方腔流的預(yù)測精度.Yin等[23]同樣基于神經(jīng)網(wǎng)絡(luò)方法,深入研究輸入特征的選擇標準以優(yōu)化模型訓(xùn)練框架.優(yōu)化后的方法應(yīng)用于周期山算例,顯著改善了流場的預(yù)測精度.張偉偉等[12,24-25]在機器學(xué)習(xí)湍流建模領(lǐng)域也開展了一系列工作,例如利用神經(jīng)網(wǎng)絡(luò)方法直接構(gòu)建基于平均流場的渦黏場預(yù)測模型[26],針對近壁區(qū)、尾流區(qū)以及遠場分別建模,成功地預(yù)測了機翼的升阻力系數(shù)、阻力分布等.Xie 等[27-29]在一系列工作中采用不同方法針對大渦模擬的亞格子應(yīng)力進行建模,在后驗測試中可以更為準確地預(yù)測湍流能譜等.

在眾多數(shù)據(jù)驅(qū)動湍流建模研究中,澳大利亞墨爾本大學(xué)Weatheritt 和Sandberg[30]采用基因表達式編程方法(gene-expression programming,GEP)[30],在最近的一些研究取得了一系列進展.GEP 方法[31]是遺傳算法的一種,可以較為方便地實現(xiàn)符號回歸,近年來已在非穩(wěn)態(tài)RANS[32]、湍流傳熱[33]、大渦模擬亞格子應(yīng)力[34]以及燃燒[35]等多個領(lǐng)域的建模中得到了初步應(yīng)用.相較于神經(jīng)網(wǎng)絡(luò)等方法的“黑箱”模型,GEP 方法可以顯式給出模型方程,因此可以很方便地應(yīng)用于工程軟件中[36-37].

本文旨在對基于GEP 方法的湍流顯式模型相關(guān)工作進行簡要總結(jié).文章的結(jié)構(gòu)安排如下:第1 部分介紹基于GEP 的模型訓(xùn)練方法及其框架;第2 部分,介紹方法的幾種典型應(yīng)用;最后第3 部分是結(jié)論與展望.

1 基因表達式編程

1.1 基本方法

GEP 方法[28]是遺傳算法的一種,基于達爾文自然選擇理論,通過種群的演化實現(xiàn)優(yōu)化過程.這里僅對GEP 方法進行簡要介紹,方法的具體細節(jié)可參考文獻[30].GEP 方法中的種群往往包含大量個體,而每個個體均代表一個候選模型.如圖2 所示,模型的基因型(genotype)是一個字符串,又被稱作染色體.該字符串由不同符號組成,包括運算符(如 +,×等)、常數(shù)項及變量符號(x,y).在將字符串表達為方程形式的過程中,首先將不同符號按照樹結(jié)構(gòu)(即表達樹,expression tree)編碼,將運算符放置在非終端節(jié)點,而常數(shù)項和變量符號放置在終端節(jié)點.進一步地,基于此樹結(jié)構(gòu)將其解碼為表現(xiàn)型(phenotype),該表達型就是一個候選模型方程.這一基于樹結(jié)構(gòu)的表達過程保證公式符合語法規(guī)則,可以被應(yīng)用于實際運算.

圖2 GEP 方法的基因型及其表達過程Fig.2 Genotype and expression process in GEP

如圖3 所示,在訓(xùn)練的開始,首先隨機生成一個包含大量個體的種群,種群中個體的優(yōu)劣可以由使用者定義的損失函數(shù)得到.類似于自然界中的自然選擇過程,種群中損失函數(shù)較低的個體的“健壯度”較好,因而有更大的概率進入下一代.進一步地,通過種群繁殖和基因變異等操作,可以得到更新一代的種群.正是通過這樣的迭代過程,整個種群向著更適應(yīng)自然選擇標準的方向演化,演化最終得到的最優(yōu)個體即為訓(xùn)練得到的模型.

圖3 GEP 方法流程圖Fig.3 Schematic for GEP method

相較于其他機器學(xué)習(xí)方法,GEP 方法的一個顯著特點是可以通過符號回歸得到顯式的模型方程,因而具有模型可解釋、易應(yīng)用的特點.另外,訓(xùn)練過程的尋優(yōu)過程主要由基因的遺傳和變異實現(xiàn).相較于梯度下降等方法,這種訓(xùn)練過程更接近全局搜索,不易收斂于局部最優(yōu)解;但同時收斂過程具有一定的隨機性,收斂速度無法得到保證,特別是可能很難得到模型中常系數(shù)的最優(yōu)解[31].

1.2 GEP 與湍流建模

GEP 方法近年來被應(yīng)用于不同領(lǐng)域的湍流建模問題.其中,既包括RANS 計算中的雷諾應(yīng)力封閉模型,也包括對湍流傳熱問題的建模.本節(jié)分別介紹這兩種典型的建模過程.

1.2.1 顯式代數(shù)應(yīng)力模型

在雷諾平均模擬中,對Navier-Stokes 方程進行雷諾平均,動量方程中將產(chǎn)生一個非封閉項,即雷諾應(yīng)力.這里ui是湍流脈動速度分量.Boussinesq(1877) 假定雷諾剪切應(yīng)力與平均速度梯度間為線性相關(guān),該線性假定可推廣到各項異性雷諾應(yīng)力張量ai j與平均應(yīng)變率張量Sij之間,即

其中,k=為湍動能,δij為 Kronecker 符號,μt為湍流動力學(xué)黏度系數(shù),νt=μt/ρ 為渦黏系數(shù).

Boussinesq 渦黏假設(shè)在大多數(shù)情況下并不適用.Pope[22]基于量綱分析和坐標變換不變量,進一步添加了非線性項,將其擴展為

式中,Ui,j為平均速度梯度,湍流時間尺度τ=1/ω=0.09k/ε ,ε 為湍動能耗散率,ω 為比耗散率.在三維流場中,存在10 個線性無關(guān)的張量基與5 個獨立的標量不變量,具體如下

式(2)中給出的顯式代數(shù)應(yīng)力模型是Boussinesq線性渦黏假設(shè)的一種修正,也是許多數(shù)據(jù)驅(qū)動湍流建?？蚣艿幕A(chǔ)[21,27-28].在GEP 湍流建模中,以標量不變量Ik、常數(shù)項、以及數(shù)學(xué)運算符號(如+,×)將作為輸入符號,通過符號回歸的方式訓(xùn)練模型系數(shù) ζk的代數(shù)表達式.

1.2.2 湍流傳熱模型

與雷諾應(yīng)力相似,對能量方程進行雷諾平均模擬時也將產(chǎn)生另一個非封閉項,即湍流熱通量傳統(tǒng)湍流傳熱的建模一般根據(jù)標準梯度擴散假設(shè)(standard gradient diffusion hypothesis,SGDH)

式中 Θ,j為平均溫度梯度,αt=νt/Prt為湍流擴散系數(shù),Prt為湍流普朗特數(shù).實際工程計算中往往將湍流普朗特數(shù)取常數(shù).例如,空氣中Prt取值為0.90 左右.然而這種常數(shù)湍流普朗特數(shù)的假設(shè)會引入較大的誤差,因此在機器學(xué)習(xí)湍流建模中,人們往往將Prt(或者 αt)作為訓(xùn)練目標.通過建立Prt隨平均溫度梯度等變量的分布函數(shù),可以提升湍流傳熱的預(yù)測精度.

另外,SGDH 中的擴散系數(shù)為各向同性,且假設(shè)湍流熱通量與平均溫度梯度線性相關(guān),這也是模型誤差的另外一個可能來源.針對這一問題,一種更為通用的模型可以將 αt替換為湍流擴散系數(shù)張量Di j[38]

為進一步提高湍流傳熱模型的計算精度,在機器學(xué)習(xí)建模中往往采用更為泛化的湍流擴散系數(shù)張量模型.本文采用 Weatheritt 等[33]對Younis 等[40]提出的顯式湍流擴散系數(shù)張量模型的簡潔表述

與雷諾應(yīng)力相似,在湍流傳熱建模中,機器學(xué)習(xí)的目標是建立模型方程(7)中g(shù)m的代數(shù)表達式.其中,輸入代數(shù)符號為標量不變量Ik和Jγ、常數(shù)項以及數(shù)學(xué)運算符號(如+,×).

1.3 模型的先驗和后驗測試

機器學(xué)習(xí)方法訓(xùn)練得到的模型在應(yīng)用前必須進行充分的測試.在湍流建模領(lǐng)域,一般將測試分為先驗(a priori)和后驗(a posteriori).所謂先驗測試,即將訓(xùn)練或者測試的高精度流場數(shù)據(jù)作為輸入特征變量代入模型方程,進而將模型預(yù)測的輸出變量與高精度數(shù)據(jù)作對比,以檢測模型對于雷諾應(yīng)力(或是湍流傳熱)的預(yù)測精度.與之相對應(yīng),后驗測試則需要將訓(xùn)練得到的模型應(yīng)用于實際RANS 計算,因此模型的輸入特征變量由RANS 計算的流場提供,而模型的預(yù)測精度則根據(jù)RANS 計算的結(jié)果進行判定.需要指出的是,一些數(shù)據(jù)驅(qū)動的湍流模型盡管可以給出令人滿意的先驗測試結(jié)果,但是在實際計算中可能導(dǎo)致較大誤差[41].

在機器學(xué)習(xí)建模中,先驗測試評價模型本身對于高精度訓(xùn)練數(shù)據(jù)的擬合和預(yù)測能力,而后驗測試則關(guān)注模型能否在實際計算中更好地預(yù)測流場信息.對于湍流建模而言,盡管先驗測試十分重要,但是后驗測試的CFD 計算中的預(yù)測效果往往才是最終的目標.

1.4 湍流建模中的損失函數(shù)

機器學(xué)習(xí)算法中損失函數(shù)的定義往往對于訓(xùn)練結(jié)果有決定性作用.對于湍流建模問題而言,也同樣需要強調(diào)損失函數(shù)的重要性.湍流建模問題中一個顯著的特點是訓(xùn)練數(shù)據(jù)與最終模型的應(yīng)用環(huán)境可能有很大區(qū)別.一方面,研究者一般認為直接數(shù)值模擬或者大渦模擬的結(jié)果為“真實值”,作為訓(xùn)練和測試數(shù)據(jù);另一方面,模型的最終需要應(yīng)用于雷諾平均模擬,而雷諾平均模擬中的湍動能、湍流耗散時間尺度等往往與訓(xùn)練數(shù)據(jù)存在較大區(qū)別.這會給損失函數(shù)的定義帶來一定的困難.本節(jié)介紹兩種不同的損失函數(shù)定義方法.

1.4.1 凍結(jié)訓(xùn)練

Akolekar 等[37]基于GEP 方法提出了凍結(jié)訓(xùn)練方法.圖4 以二方程雷諾平均模型為例介紹凍結(jié)訓(xùn)練的流程.由直接數(shù)值模擬等高可信度流場,可以得到平均速度、平均速度梯度張量、雷諾應(yīng)力、湍動能等數(shù)據(jù),并以此作為“真實值”.在此基礎(chǔ)上,將“凍結(jié)”的平均速度和湍動能代入模型輸運方程,如k-ω模型中的ω方程,得到適用于雷諾平均模擬的湍流耗散時間尺度.采用這種湍流特征時間尺度對速度梯度張量進行無量綱化,可以得到1.2.1 小節(jié)中的張量基函數(shù)和標量不變量作為模型的輸入變量.進一步地,使用GEP 方法訓(xùn)練基于無量綱速度梯度張量基函數(shù)和標量不變量的顯式代數(shù)應(yīng)力模型.

圖4 凍結(jié)訓(xùn)練流程示意圖Fig.4 Schematic for ‘frozen’ training method

在這一訓(xùn)練過程中,通過求解k-ω模型中的ω方程得到湍流耗散時間尺度,用于進一步的速度梯度張量無量綱化.在求解過程中,平均速度場和湍動能來自于“凍結(jié)的”高精度訓(xùn)練數(shù)據(jù),在迭代過程中不發(fā)生變化.這里僅對ω進行迭代更新,直至收斂.因此該方法被稱為凍結(jié)訓(xùn)練.

凍結(jié)訓(xùn)練方法中采用求解模型方程得到的湍流耗散時間尺度進行無量綱化,而不是直接使用高精度訓(xùn)練數(shù)據(jù)中的特征時間尺度.這是因為直接數(shù)值模擬等方法得到的湍流耗散時間尺度與實際雷諾平均模擬運算中的結(jié)果往往存在較大差異.這種差異可以認為是高精度訓(xùn)練數(shù)據(jù)與雷諾平均模擬之間存在一定的相容性問題.換而言之,完全依據(jù)直接數(shù)值模擬等高精度流場訓(xùn)練得到的雷諾應(yīng)力模型可能并不適用于雷諾平均模擬,這種情況下盡管模型可能在先驗測試中準確度較高,但在后驗測試中仍然存在較大誤差.因此,必須在模型訓(xùn)練過程中考慮這種差異,通過凍結(jié)方法增強模型在實際雷諾平均模擬中的適用性,提升其后驗精度.

凍結(jié)訓(xùn)練以修正模型預(yù)測的雷諾應(yīng)力為目標,對應(yīng)的損失函數(shù)是通過對比高精度訓(xùn)練數(shù)據(jù)中的雷諾應(yīng)力和模型給出的雷諾應(yīng)力來得到.通過求解模型輸運方程,凍結(jié)訓(xùn)練引入適用于雷諾平均模擬的特征時間尺度進行無量綱化,從而可以在一定程度上緩解高精度訓(xùn)練數(shù)據(jù)與RANS 模擬的相容性問題.然而,最近一些學(xué)者研究發(fā)現(xiàn)[42-43],即便將直接數(shù)值模擬得到的雷諾應(yīng)力代入RANS 計算中,得到的平均流場也可能會產(chǎn)生較大的誤差.Wu 等[41]還發(fā)現(xiàn),這種預(yù)測誤差在雷諾數(shù)增高時會進一步增大.這種現(xiàn)象導(dǎo)致數(shù)據(jù)驅(qū)動的RANS 模型很難在后驗驗證中準確預(yù)測平均流動.

1.4.2 雙向耦合方法

針對湍流建模中常見的高精度訓(xùn)練數(shù)據(jù)和RANS計算的相容性問題,Zhao 等[44]首次提出了機器學(xué)習(xí)和CFD 雙向耦合的模型訓(xùn)練方法.雙向耦合方法的最大特點是訓(xùn)練過程中模型的優(yōu)劣通過實際RANS 計算的結(jié)果來評估,而不是僅僅基于高精度訓(xùn)練數(shù)據(jù)判定對于雷諾應(yīng)力的預(yù)測效果.如圖5 所示,對于每一代的候選模型,將一系列模型方程讀入到RANS 求解器中,并基于該模型實現(xiàn)RANS 計算;進一步將計算的結(jié)果返回GEP 訓(xùn)練過程,以評估模型方程的損失函數(shù).由于在訓(xùn)練過程中,每一代模型的評估均需要GEP 算法和RANS 計算共同作用,因此稱之為機器學(xué)習(xí)和CFD 雙向耦合的訓(xùn)練方法.雙向耦合訓(xùn)練方法的步驟可以總結(jié)為以下幾個步驟:

圖5 雙向耦合方法訓(xùn)練流程示意圖Fig.5 Schematic for CFD-driven machine learning method

(1)初始化種群,得到一系列模型方程;

(2)種群中的一系列候選模型均由CFD 運算進行評估,這一過程可以并行操作,具體步驟如下:

①模型方程寫入文件,再由RANS 求解器讀入此文件;

② RANS 求解器基于GEP 算法提供的候選模型運行CFD 模擬;

③RANS 計算的結(jié)果分兩類.如果RANS 收斂,模型的損失函數(shù)可由收斂后流場計算得到;如果RANS 不收斂,模型損失函數(shù)設(shè)為無窮大;

④ 將CFD 計算得到的模型損失函數(shù)返回GEP 算法,以評估候選模型;

⑤ GEP 算法迭代,模型種群繁殖、變異,進入下一代,直到得到期望的模型.

相比于凍結(jié)方法,雙向耦合方法具有以下特點.首先,雙向耦合方法的訓(xùn)練目標靈活,可以根據(jù)RANS 計算的結(jié)果設(shè)定不同的損失函數(shù)作為優(yōu)化目標;相對而言,凍結(jié)方法只能對比模型預(yù)測的雷諾應(yīng)力與高精度訓(xùn)練數(shù)據(jù).其次,雙向耦合方法對于訓(xùn)練數(shù)據(jù)的需求量較少,實驗測量得到的少許平均流場數(shù)據(jù)即可作為訓(xùn)練數(shù)據(jù);相對而言,凍結(jié)訓(xùn)練往往需要高精度數(shù)值模擬得到的全場雷諾應(yīng)力等高階湍流統(tǒng)計量.最后,由于雙向耦合方法得到的模型在訓(xùn)練過程中就已經(jīng)經(jīng)過RANS 計算的檢驗,因此模型的收斂性和后驗預(yù)測精度均有一定保證;相對而言,凍結(jié)方法得到的模型在應(yīng)用于真實RANS 計算時,對于平均流場的預(yù)測精度可能存在較大誤差.雙向耦合方法面臨的主要挑戰(zhàn)在于計算效率方面.與凍結(jié)訓(xùn)練相比,盡管雙向耦合方法在GEP 訓(xùn)練中所需的種群代數(shù)相差不大(一般約為數(shù)百代),但是由于其迭代過程中每一代候選模型的評估均需運行RANS 計算,因此訓(xùn)練所需的計算量往往遠大于凍結(jié)方法.

2 應(yīng)用算例

第1 部分介紹了基因表達式編程方法以及如何基于此方法訓(xùn)練湍流和傳熱模型.接下來,介紹GEP 方法在幾種湍流建模問題中的具體應(yīng)用.

2.1 尾流混合問題的代數(shù)應(yīng)力模型

2.1.1 統(tǒng)計穩(wěn)態(tài)流動的RANS 模型

在航空發(fā)動機內(nèi)流中的葉片尾沿下游,葉片壓力側(cè)和吸力側(cè)的來流在強剪切層的主導(dǎo)下產(chǎn)生強烈混合,這是典型的尾流混合問題.尾流混合往往導(dǎo)致顯著的流動損失,因此準確模擬尾流混合問題對于內(nèi)流預(yù)測和葉片設(shè)計十分重要,然而基于Boussinesq假設(shè)的常用工程湍流模型對于尾流損失的預(yù)測精度往往不甚理想[45].最近幾年來,GEP 方法在一系列的工作中[37,44,46-47]被應(yīng)用于尾流混合問題,針對如圖6所示的高壓渦輪和低壓渦輪平面葉柵算例,目標為發(fā)展具有較高預(yù)測精度和較強泛化能力的尾流混合模型.

圖6 航空發(fā)動機內(nèi)流葉柵尾流混合算例示意圖Fig.6 Simulation setup for cases

表1 給出了一系列渦輪葉柵算例的參數(shù).基于直接數(shù)值模擬或是大渦模擬提供的高精度流場數(shù)據(jù),這些算例被用于模型的訓(xùn)練和測試.以典型航空發(fā)動機工況下的高壓渦輪(表1 算例A) 為研究目標,Weatheritt 等[47]采用1.2.1 小節(jié)中介紹的凍結(jié)方法訓(xùn)練顯式代數(shù)應(yīng)力模型.訓(xùn)練中損失函數(shù)設(shè)置為

表1 渦輪葉柵尾流混合算例Table 1 Parameters of turbine wake mixing cases

即在N個網(wǎng)格點上,GEP 模型得到的雷諾應(yīng)力各向異性部分與高精度訓(xùn)練數(shù)據(jù)的相對誤差.訓(xùn)練得到的模型方程形式如下

對這一模型進行先驗測試,Weatheritt 等[47]發(fā)現(xiàn)相較于傳統(tǒng)Boussinesq 假設(shè),該模型可以明顯改善對于尾流區(qū)域雷諾應(yīng)力的預(yù)測.但是,由于存在大量高階非線性項,該模型在后驗測試中很難得到收斂的流場.

為進一步發(fā)展魯棒性更強的尾流混合模型,Zhao 等[44]提出了1.2.2 小節(jié)中介紹的雙向耦合訓(xùn)練方法,并將其應(yīng)用于高壓渦輪算例A.不同于式(10)中給出的凍結(jié)訓(xùn)練損失函數(shù),雙向耦合方法訓(xùn)練過程中損失函數(shù)的設(shè)置較為靈活,不僅可以對比雷諾應(yīng)力等高階統(tǒng)計量,更可以選擇實驗、直接數(shù)值模擬等高可信度數(shù)據(jù)中可提供的任意感興趣的平均流場信息.以渦輪尾流混合問題為例,這里選擇根據(jù)RANS 模擬結(jié)果中的平均尾流損失設(shè)置損失函數(shù).具體形式如下

這里 Ω (y) 為尾流損失曲線,,po和pt(y) 分別代表來流總壓、出口靜壓和當?shù)乜倝?Δx是GEP 模型預(yù)測的尾流損失曲線與高可信度流場數(shù)據(jù)的相對誤差,而最終損失函數(shù)JCFD則是兩個不同位置的相對誤差之和.值得注意的是,在這樣的雙向耦合訓(xùn)練過程中,僅從高可信度流場中提取少數(shù)平均尾流損失曲線作為訓(xùn)練數(shù)據(jù),而不需要凍結(jié)方法中必須的全場雷諾應(yīng)力信息,因此對于訓(xùn)練數(shù)據(jù)的要求大大降低.

采用雙向耦合方法訓(xùn)練并設(shè)置式(12)中的損失函數(shù),Zhao 等[44]得到的尾流混合模型如下

將此模型與凍結(jié)方法獲得的模型式(10) 作對比,發(fā)現(xiàn)雙向耦合方法得到的模型方程中高階非線性項大大減少.其原因在于這些高階項往往導(dǎo)致實際RANS 計算難以收斂,因此在雙向耦合訓(xùn)練過程中被自然選擇所淘汰.由此可見,雙向耦合方法獲得的模型一般具有更強的收斂性和魯棒性.

將式(11)中的雙向耦合模型應(yīng)用于后驗測試,并將其與去除高階項后的簡化凍結(jié)模型的預(yù)測結(jié)果作對比.如圖7(a) 所示,相較于作為基準模型的k-ωSST 模型和凍結(jié)模型,雙向耦合方法在后驗測試中可以顯著提高對于式(12)中尾流損失 的預(yù)測精度.進一步對模型進行分析,Zhao 等[44]指出雙向耦合模型對于尾流損失預(yù)測的提升,其主要原因在于方程中占主導(dǎo)的項的系數(shù)幅值更大(2 .57 ＞ 1 .334).考慮到,這說明該模型在尾流區(qū)域?qū)е赂鼜姷臄U散,因而雙向耦合模型預(yù)測的尾流損失的峰值和尾流寬度均更接近于大渦模擬的結(jié)果.由此可見,GEP 方法訓(xùn)練得到的湍流模型是顯式方程,這樣可以分析、解釋模型預(yù)測精度提升的作用機制.

圖7 不同渦輪葉柵算例(見表1)中的尾流損失剖面Fig.7 Kinetic wake loss profiles from test cases in Table 1

為進一步驗證雙向耦合模型的泛化能力,還將其應(yīng)用于表1 所示的算例B,C 和D 中.這些算例與訓(xùn)練算例A 相比,其雷諾數(shù)、出流馬赫數(shù)、攻角以及幾何外形均有很大差別.如圖7(b)～圖7(d)所示,在這一系列后驗測試中,這種雙向耦合方法訓(xùn)練得到的顯式代數(shù)應(yīng)力模型可以較為準確地預(yù)測尾流混合導(dǎo)致的流動損失.因此該模型具有較強的泛化能力.

2.1.2 非穩(wěn)態(tài)流動的unsteady RANS 建模策略

上一節(jié)討論了針對統(tǒng)計穩(wěn)態(tài)流動的建模過程,可以看到GEP 方法能有效提升RANS 對于尾流混合的預(yù)測精度.但是當流動處于非穩(wěn)態(tài)時,往往應(yīng)采用非定常雷諾平均(unsteady RANS)模擬,而建模策略可能也需要相應(yīng)調(diào)整.

航空發(fā)動機內(nèi)流領(lǐng)域一種常見的非穩(wěn)態(tài)現(xiàn)象是兩排葉片之間存在的相對運動導(dǎo)致的流動周期性變化.Akolekar 等[48]針對來流為周期性尾流擾動的低壓渦輪葉柵展開研究,探討機器學(xué)習(xí)湍流建模的策略.圖8 給出了鎖相平均的葉柵湍動能云圖.將一個來流擾動周期分為20 個不同的相,由圖8 可發(fā)現(xiàn)周期性來流尾流對渦輪葉柵流場造成顯著影響,不同相的流動之間呈現(xiàn)顯著差別.

圖8 來流尾流擾動下低壓渦輪葉柵的鎖相平均湍動能云圖[48]Fig.8 Phase-lock averaged TKE contours for LPT flow disturbed by incoming wakes[48]

針對這種非穩(wěn)態(tài)流動,湍流建模策略可有幾種:(1)針對不同相的流場單獨建模;(2)根據(jù)時間平均后的流場建模;(3)根據(jù)一定的流動物理分組建模.顯然,針對不同相分別建模會給出適用于各相的一系列湍流模型,但是模型間的切換會給使用者的實際模擬帶來不便.Akolekar 等[46]進一步發(fā)現(xiàn),針對時間平均的流場建模,雖然可以一定程度改善預(yù)測精度,但是由于忽略了各相間的差異導(dǎo)致其精度仍然有限.相對而言,最優(yōu)的策略是根據(jù)流動物理將不同相分組,進而針對各組具有相似物理的流動分別建模.

在這種根據(jù)流動物理對不同相進行分組的過程中,GEP 方法可以起到關(guān)鍵作用.針對低頻來流尾流擾動下的低壓渦輪葉柵流動,Akolekar 等[48]首先對20 相分別建模.得到這些單相模型后,一方面可以分別測試各模型的預(yù)測精度,發(fā)現(xiàn)單相訓(xùn)練得到的模型在不同相的測試中表現(xiàn)并不理想.另一方面,也可以根據(jù)GEP 方法顯式給出的方程形式,對不同相的模型進行分析.在分析過程中,Akolekar 等[48]發(fā)現(xiàn)不同相的模型方程大致可分為兩組.對于其中一組(包含大約10 相),其模型方程中項的系數(shù)幅值遠大于另外一組(包含剩余10 相),這說明這些相中需要引入相對更強的湍流擴散.進一步分析流場特征,可以發(fā)現(xiàn)第一組的流動主要特征是來流尾流與葉片邊界層相互作用,進而導(dǎo)致葉柵尾流混合強度顯著增強;而另外一組流動中來流尾流與葉片則沒有明顯的相互作用,相應(yīng)需要的湍流擴散修正較弱.這一分組過程首先依據(jù)GEP 訓(xùn)練的單相模型方程得到,進一步得到了流場分析的驗證.在此分組的基礎(chǔ)上,可以針對兩組流動分別建模,得到適用于對應(yīng)組的兩個模型.在先驗測試中,這種基于流動物理的模型可顯著提升模型的預(yù)測精度.

2.2 湍流傳熱模型

近年來,GEP 方法也被應(yīng)用于溫度或標量場的建模之中.本文首先針對豎直平板間自然對流這一較為簡單的算例,介紹基于GEP 的標量系數(shù)湍流傳熱模型;然后針對更為復(fù)雜的三維橫向流中的射流算例,介紹張量系數(shù)湍流傳熱建模.

2.2.1 自然對流的湍流標量系數(shù)建模

自然對流在室內(nèi)通風(fēng)環(huán)境(如雙層玻璃窗)、電子設(shè)備冷卻等場景中十分常見.圖9 給出了以浮力驅(qū)動豎直充分發(fā)展槽道流中的自然對流計算域示意圖,其中兩個豎直板沿x和z方向采用周期性邊界,y方向左右兩側(cè)為不同溫度的等溫壁面.兩壁面中,左側(cè)高溫無量綱化后為 Θh=0.5 ,右側(cè)低溫為 Θc=-0.5,有限寬度H=2h.對于豎直槽道流中空氣流動,流場特性由Rayleigh 數(shù)Ra=gβΔΘH3/(υκ) 決定.其中ΔΘ=Θh-Θc是兩個豎直板間無量綱的溫度差,g為重力加速度,β 為體積膨脹系數(shù),ν 為運動黏度,κ 為擴散系數(shù),空氣普朗特數(shù)Pr=ν/κ=0.709 .

圖9 以豎直槽道流中的自然對流算例的計算域Fig.9 Computational domain of natural convection case in a differentially heated vertical channel

Ng 等[49]通過直接數(shù)值模擬(DNS)中得到了該算例 的平均流場、雷諾應(yīng)力和湍流熱通量.Xu 等[50]以此數(shù)據(jù)集中的算例為訓(xùn)練數(shù)據(jù),提取湍流渦黏系數(shù),并以垂直壁面方向的熱通量為學(xué)習(xí)目標.基于GEP方法,通過最小化損失函數(shù)

得到模型方程f1=0.969+2J0,嵌入到SGDH 中為

根據(jù) αt=vt/Prt,這里湍流普朗特系數(shù)等價于Prt=1/(0.969+2J0).與此對比,本文以常用的常湍流普朗特數(shù)為基準模型(baseline model),其中Prt=0.9 .

圖10 湍流普朗特數(shù)的先驗預(yù)測Fig.10 Prediction of turbulent Prandtl number

在后驗測試中,Xu 等[50]在給定DNS 的速度場相關(guān)數(shù)據(jù)的基礎(chǔ)上,求解溫度的對流-擴散輸運方程,來測試湍流傳熱模型的預(yù)測效果.該模型被應(yīng)用于3 個不同Ra算例,圖11 分別展示了后驗測試中基準模型和GEP 模型預(yù)測的平均溫度剖面和垂直壁面方向的熱通量.可以看到,基準模型在平均溫度預(yù)測上有較大的誤差,且該誤差隨著Ra增大而增大.與此對比,GEP 模型能夠有效的改善熱通量的預(yù)測,進而修正平均溫度剖面,并且計算結(jié)果與DNS 的真實值趨于一致.需要強調(diào)的是,上述先驗和后驗測試均采用了Ra=2.0×107的數(shù)據(jù)訓(xùn)練得到的GEP 模型,然后在不同下測試該模型.以上結(jié)果表明,GEP 模型不僅在訓(xùn)練算例中表現(xiàn)優(yōu)秀,而且能夠在其他數(shù)據(jù)集的交叉驗證中驗證其其有效性.因此,該模型具有較好的泛化能力.

圖11 平均溫度剖面和垂直壁面方向的熱通量的后驗預(yù)測Fig.11 Prediction of mean temperature profiles and wall-normal heat flux

2.2.2 橫向射流的湍流張量擴散系數(shù)建模

三維橫向流中的射流 (jet in cross-flow,JICF)在航空發(fā)動機等渦輪葉片中有廣泛應(yīng)用,準確預(yù)測其流場和溫度場是工業(yè)界亟待解決的難題.圖12 給出了典型JICF 的算例設(shè)置示意圖.其中,x為流向,y為垂直壁面方向,z為展向,射流孔直徑為d,橫向主流無量綱溫度 Θin,冷卻射流無量綱溫度 Θj.流場的主要特征參數(shù)為雷諾數(shù)Rej=Ujd/ν ,入射角 φ,吹風(fēng)比(blowing ratio)BR=Uj/Uin.算例中高溫橫向主流(freestream inflow)與從氣膜孔流出的冷卻氣體射流(coolant inflow)相互作用,形成摻混區(qū)域.由于摻混區(qū)域存在逆梯度輸運(counter-gradient transport)、擬序結(jié)構(gòu)形成與破碎等復(fù)雜三維非定常流動現(xiàn)象,基于RANS 對其湍流傳熱的預(yù)測難度非常大.

圖12 橫流中的射流示意圖Fig.12 The schematic description of jet in crossflow case

Weatheritt 等[33]基于Bodart 等[51]的JICF 大渦模擬數(shù)據(jù)進行GEP 湍流傳熱建模,其中展向?qū)挾葹?d,豎向高度為 8.62d,φ=30°,BR=1.0,Rej=5.4×103,這里稱該算例為Bodart 算例.由于流場的復(fù)雜性,標準梯度擴散假設(shè)難以準確預(yù)測傳熱過程,因此必須采用湍流擴散系數(shù)張量模型對三維流場中的熱通量進行數(shù)據(jù)驅(qū)動建模.

可以發(fā)現(xiàn)張量系數(shù)相關(guān)度的排序為

在此基礎(chǔ)上,GEP 模型訓(xùn)練中符號回歸的方程設(shè)置為

此外,損失函數(shù)為

其中,V為三維流場中數(shù)據(jù)集所有離散點數(shù).經(jīng)過機器學(xué)習(xí)訓(xùn)練之后,得出如下模型

該模型被應(yīng)用于Bodart 算例.圖13 給出了流向x/d=20處截面(見圖12)上垂直壁面方向的熱通量預(yù)測結(jié)果,并將GEP 模型、基準模型(SGDH) 與LES 數(shù)據(jù)作對比.在不同展向位置(z/d=0.50,0.67,1.10),與SGDH 相比,GEP 模型提高了對兩個熱通量分量vθ和wθ 的預(yù)測精度.此外,沿垂直壁面方向上,GEP模型預(yù)測的峰值點位置也與LES 結(jié)果較好吻合.

圖13 流向位置x/d=20 的熱通量預(yù)測Fig.13 The prediction of heat flux vector at x/d=20

此外,該模型還被應(yīng)用于Sakai 等[52]的JICF 大渦模擬算例(Sakai 算例),以測試基于Bodart 算例的湍流張量擴散系數(shù)模型的泛化能力.與Bodart 算例相比,Sakai 算例的幾何尺寸和主要特征參數(shù)均不同,其中,展向?qū)挾葹?3d,豎向高度為 5d,φ =35°,Rej=1.5×104.圖14 比較了兩組不同吹風(fēng)比BR=0.5 和BR=1.0算例中壁面絕熱效率η的流向分布,其中,壁面絕熱效率定義為.當BR=1.0(與Bodart 算例一致),GEP 模型對橫向平均和中線處 η 符合高精度數(shù)據(jù)的總體趨勢,在摻混區(qū)域與參考數(shù)據(jù)非常接近.當BR=0.5,GEP 模型對橫向平均和中線處 η 的預(yù)測精度均優(yōu)于基準模型.

圖14 壁面絕熱效率的流向分布Fig.14 Wall adiabatic effectiveness streamwise distribution

綜上所述,無論是針對湍流標量系數(shù)(如湍流普朗特數(shù))還是對湍流張量擴散系數(shù)進行建模,GEP 方法均能提升對湍流傳熱的預(yù)測精度,且具有一定程度的泛化能力.

2.3 GEP 方法在湍流建模中的其他應(yīng)用

自Weatheritt 和Sandberg[30]首次將GEP 應(yīng)用于湍流建模以來,此方法已在多個相關(guān)領(lǐng)域得到了成功應(yīng)用.2.1 節(jié)和2.2 節(jié)重點介紹了針對代數(shù)應(yīng)力模型和湍流傳熱問題的典型建模過程.在此基礎(chǔ)上,在本節(jié)簡要介紹GEP 方法在湍流建模領(lǐng)域的一些拓展和其他應(yīng)用.

與1.2.1 小節(jié)中的RANS 建模過程非常類似,GEP 方法還可應(yīng)用于大渦模擬亞格子應(yīng)力建模.Schoepplein 等[35]首次將GEP 方法應(yīng)用于湍流預(yù)混火焰的亞格子應(yīng)力建模.與RANS 建模中以平均流場為訓(xùn)練數(shù)據(jù)不同,亞格子應(yīng)力的建?；跒V波后的直接數(shù)值模擬瞬時流場.該訓(xùn)練采用凍結(jié)方法,所得到的亞格子應(yīng)力代數(shù)表達式在先驗測試中具有較高的準確度.但是,該模型并未經(jīng)過后驗測試,且模型的泛化能力也需要進一步驗證.最近,Reissmann 等[34]嘗試將雙向耦合方法[44]應(yīng)用于泰勒-格林渦的亞格子應(yīng)力建模中.在后驗測試中,所得模型可較為準確預(yù)測湍動能等統(tǒng)計量演化過程.當然,這一模型對充分發(fā)展湍流的適用性還有待一步驗證.

在最近的工作中,GEP 方法還被嘗試應(yīng)用于邊界層轉(zhuǎn)捩問題的建模[53].基于層流動能(laminar kinetic energy,LKE)方法[54],Akolekar 等[53]針對LKE 方程中的層流渦黏系數(shù)等關(guān)鍵參數(shù)提出修正的代數(shù)模型.在后驗測試中,該模型實現(xiàn)了對于低壓渦輪葉柵中分離引起的邊界層轉(zhuǎn)捩的準確預(yù)測.雖然該模型的魯棒性等問題還有待進一步研究,但該結(jié)果顯示了GEP 方法在邊界層轉(zhuǎn)捩等復(fù)雜流動的建模中具有良好的應(yīng)用前景.

3 結(jié)論與展望

近年來,機器學(xué)習(xí)輔助的湍流建模研究發(fā)展迅速,然而在訓(xùn)練數(shù)據(jù)和訓(xùn)練方法、模型可解釋性和泛化能力等方面仍然存在較大挑戰(zhàn).本文討論的一系列工作基于GEP 方法在這一領(lǐng)域取得了一定進展,實現(xiàn)了多種類型的湍流建模.在訓(xùn)練數(shù)據(jù)方面,這些工作關(guān)注具有較強實際應(yīng)用背景的復(fù)雜流動,以發(fā)展實用工程湍流模型為目標.在訓(xùn)練方法方面,通過引入雙向耦合方法等模型訓(xùn)練框架,考慮RANS計算與高精度訓(xùn)練數(shù)據(jù)的差別,強調(diào)模型的實際預(yù)測能力.另外,由于GEP 方法可以提供顯式模型方程,因此模型可解釋性較強.

值得注意的是,目前大多數(shù)機器學(xué)習(xí)方法得到的湍流模型距離真實工程應(yīng)用還存在明顯的差距.接下來的工作重點之一是提高模型的泛化能力.一種具有較強應(yīng)用前景的數(shù)據(jù)驅(qū)動模型應(yīng)該可以較好地預(yù)測多種類型流動(如遠離壁面的自由剪切流動和近壁面湍流)及多種流動物理(如強混合、轉(zhuǎn)捩、流動分離等).此外,針對雙向耦合方法,還可以通過優(yōu)化GEP 算法的收斂速度或種群大小、改善CFD求解器等方式來提升訓(xùn)練的計算效率.以上述工作重點為目標,需要進一步完善現(xiàn)有的訓(xùn)練數(shù)據(jù)和訓(xùn)練方法.