基于基因表達式編程的數據驅動湍流建模

趙耀民 *, 徐曉偉

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

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

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

引言

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

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

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

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

圖1 機器學習湍流建模示意圖Fig.1 Schematic for turbulence modelling with machine learning methods

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

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

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

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

1 基因表達式編程

1.1 基本方法

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

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

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

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

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

1.2 GEP 與湍流建模

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

1.2.1 顯式代數應力模型

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

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

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

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

式(2)中給出的顯式代數應力模型是Boussinesq線性渦黏假設的一種修正,也是許多數據驅動湍流建模框架的基礎[21,27-28].在GEP 湍流建模中,以標量不變量Ik、常數項、以及數學運算符號(如+,×)將作為輸入符號,通過符號回歸的方式訓練模型系數 ζk的代數表達式.

1.2.2 湍流傳熱模型

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

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

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

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

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

1.3 模型的先驗和后驗測試

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

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

1.4 湍流建模中的損失函數

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

1.4.1 凍結訓練

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

圖4 凍結訓練流程示意圖Fig.4 Schematic for ‘frozen’ training method

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

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

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

1.4.2 雙向耦合方法

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

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

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

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

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

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

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

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

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

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

2 應用算例

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

2.1 尾流混合問題的代數應力模型

2.1.1 統計穩態流動的RANS 模型

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

圖6 航空發動機內流葉柵尾流混合算例示意圖Fig.6 Simulation setup for cases

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

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

即在N個網格點上,GEP 模型得到的雷諾應力各向異性部分與高精度訓練數據的相對誤差.訓練得到的模型方程形式如下

對這一模型進行先驗測試,Weatheritt 等[47]發現相較于傳統Boussinesq 假設,該模型可以明顯改善對于尾流區域雷諾應力的預測.但是,由于存在大量高階非線性項,該模型在后驗測試中很難得到收斂的流場.

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

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

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

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

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

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

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

2.1.2 非穩態流動的unsteady RANS 建模策略

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

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

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

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

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

2.2 湍流傳熱模型

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

2.2.1 自然對流的湍流標量系數建模

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

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

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

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

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

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

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

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

2.2.2 橫向射流的湍流張量擴散系數建模

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

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

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

可以發現張量系數相關度的排序為

在此基礎上,GEP 模型訓練中符號回歸的方程設置為

此外,損失函數為

其中,V為三維流場中數據集所有離散點數.經過機器學習訓練之后,得出如下模型

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

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

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

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

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

2.3 GEP 方法在湍流建模中的其他應用

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

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

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

3 結論與展望

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

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