機(jī)器學(xué)習(xí)在湍流模型構(gòu)建中的應(yīng)用進(jìn)展

張偉偉,朱林陽,劉溢浪,寇家慶

(西北工業(yè)大學(xué) 航空學(xué)院,西安 710072)

0 引 言

湍流問題涉及到工程中的諸多領(lǐng)域，其重要性不言而喻。人們對于湍流的探索和研究也從未間斷?？焖倩兯枷隱1-2]和各種標(biāo)度率的提出以及擬序結(jié)構(gòu)[3]的演示促進(jìn)了人們對湍流的認(rèn)識。隨著計算機(jī)的迅速發(fā)展以及粒子圖像測速儀(PIV)的應(yīng)用，人們對湍流的研究更加深入細(xì)致。我國的湍流研究在前輩學(xué)者的引領(lǐng)下起步很早[4]，近些年也取得了可觀的研究成果。佘振蘇教授等[5]提出了基于結(jié)構(gòu)系綜理論(SED,Structural Ensemble Dynamics)的新模型(SED-SL)，陳十一院士團(tuán)隊(duì)[6]提出了帶約束的湍流大渦模擬方法，何國威院士等[7]綜述了在歐拉與拉格朗日兩種參考系下的各種時空關(guān)聯(lián)模型，等等。經(jīng)過前人的不斷努力，人們對湍流的認(rèn)識不斷加深，但湍流的本質(zhì)問題[8]以及在工程中的高效應(yīng)用仍未得到很好的解決，湍流仍然是航空航天航海等領(lǐng)域工程成敗的關(guān)鍵瓶頸之一。湍流的本質(zhì)、湍流演變的物理規(guī)律，以及如何在工程中更好地應(yīng)用湍流，一直是研究者們不斷探索的研究方向。

湍流研究的數(shù)值計算方法可根據(jù)網(wǎng)格分辨尺度分為直接數(shù)值模擬(DNS)、大渦模擬(LES)和雷諾平均應(yīng)力模型(RANS)。隨著計算機(jī)性能的提升和并行算法的開發(fā)，DNS和LES越來越多地在學(xué)術(shù)界和工業(yè)界得到應(yīng)用。但是，對于復(fù)雜的幾何外形和高雷諾數(shù)流動而言，急劇增長的網(wǎng)格量很大程度上限制了這兩種方法的應(yīng)用范圍，對于飛行器在飛行包線中的計算恐怕要到2030年才能實(shí)現(xiàn)[9]。相比之下，RANS模型雖然精度不及前兩種方法，但其易用性和高效性使之在工程實(shí)踐中得到廣泛應(yīng)用[10]。特別是在航空航天領(lǐng)域，流動的雷諾數(shù)普遍較高，湍流計算主要依賴于現(xiàn)有的RANS模型。

RANS模型大多基于渦黏假設(shè)，在附著流中容易取得較好的計算結(jié)果。對于分離流等復(fù)雜流動，雷諾應(yīng)力與應(yīng)變之間不再是簡單的線性相關(guān)，其各向異性使得湍流的模型化變得明顯困難[11]，因此，RANS模型的計算結(jié)果往往偏差較大而不可靠。另一方面，模型中的經(jīng)驗(yàn)性參數(shù)往往是依據(jù)某些特定的流動確定的，一定程度上也增加了模型的不確定度，繼而影響模型的適用性。目前常用的RANS模型還是20世紀(jì)提出的，主要有Baldwin-Lomax(BL)模型[12]、k-ε模型[13]、k-ω模型[14]和Spalart-Allmaras(SA)模型[15]等。另外，可視具體工程問題采用SST模型[16]、k-ω2模型[17]、k-τ模型[18]以及一些修正的k-ε模型[19]等。與一階矩封閉模式不同，雷諾應(yīng)力輸運(yùn)模型(RSTM)是基于雷諾應(yīng)力輸運(yùn)方程的二階矩封閉模式[20]。該模型雖然精度較高，但計算效率和魯棒性有待進(jìn)一步提高，尚未在航空工程中廣泛采用，許多CFD研究者也在致力于這一問題的改善。近年來，佘振蘇等人[21-22]基于李群理論提出了解析的應(yīng)力長函數(shù)表達(dá)式以描述壁湍流中的多層結(jié)構(gòu)。從歷史過程可以看出，RANS模型的發(fā)展已經(jīng)到了平臺期，難以取得新的突破來克服其固有的缺陷，湍流問題的求解亟待新思路、新方法的涌現(xiàn)。

實(shí)際上，針對復(fù)雜靜/動力學(xué)系統(tǒng)模型化主要有兩個途徑。第一個是基于理論的模型架構(gòu)，即根據(jù)物理問題的控制方程，建立理想的系統(tǒng)描述。這種模型通常需要研究者對物理過程有深刻的理解，并將其轉(zhuǎn)化為數(shù)學(xué)模型，現(xiàn)行的湍流模型大都是采用這個途徑。第二個則是數(shù)據(jù)驅(qū)動方法，即根據(jù)系統(tǒng)仿真或試驗(yàn)中的樣本數(shù)據(jù)，直接構(gòu)造黑箱或灰箱模型。近年來，隨著計算機(jī)性能和精細(xì)化流動測試手段的發(fā)展，研究者逐步能夠獲得高精度、高時空分辨率的流場信息，或者直接通過開放的平臺(如http://turbulence.pha.jhu.edu/)獲取。如何高效地利用這些大數(shù)據(jù)，從中提取出關(guān)鍵信息，并指導(dǎo)流體力學(xué)的發(fā)展，已經(jīng)成為研究者關(guān)注的焦點(diǎn)。作為處理和分析數(shù)據(jù)的主要手段，數(shù)據(jù)挖掘、統(tǒng)計學(xué)習(xí)和機(jī)器學(xué)習(xí)等技術(shù)，則為開展此類研究提供了重要基礎(chǔ)。機(jī)器學(xué)習(xí)通過一些算法從數(shù)據(jù)中建立模型，使之具備一定的判斷和預(yù)測能力。常用的算法有徑向基神經(jīng)網(wǎng)絡(luò)(RBFNN)、隨機(jī)森林(RF)、支持向量機(jī)(SVM)和神經(jīng)網(wǎng)絡(luò)(NNs)等。這些算法已廣泛應(yīng)用于語音和圖像識別[23-24]、信號處理[25]以及降階處理[26]等領(lǐng)域。

流體力學(xué)研究者也已將這些方法成功用于各自領(lǐng)域的相關(guān)研究工作中，例如偏微分方程求解以及非線性動力學(xué)特征的模擬等。在湍流領(lǐng)域，國內(nèi)外的專家學(xué)者也加快了機(jī)器學(xué)習(xí)方法在湍流建模中的應(yīng)用步伐。早期的研究工作主要是探討與湍流相關(guān)的變量之間的關(guān)系，采用的方法也較為簡單[27-28]。實(shí)際上，正式采用數(shù)據(jù)驅(qū)動方法來改善或替代RANS模型的研究主要是近幾年開展的[29]。逐漸地，研究者開始側(cè)重于對RANS模型的改善以及更復(fù)雜機(jī)器學(xué)習(xí)方法的應(yīng)用。Tracey等人[30]針對二維及三維流動，構(gòu)建了替代SA模型中源項(xiàng)的神經(jīng)網(wǎng)絡(luò)模型。除了以RANS模型控制方程中的某些項(xiàng)作建模對象，研究者還針對修正系數(shù)分布或附加源項(xiàng)建立數(shù)據(jù)驅(qū)動模型并以此來改善原RANS模型的計算精度[31-33]。Wang等人[34]針對RANS模型和高分辨率數(shù)據(jù)之間雷諾應(yīng)力偏差構(gòu)建模型，提高了RANS模型的準(zhǔn)確性。上述工作主要是利用機(jī)器學(xué)習(xí)方法來完善RANS模型，以獲得更高的計算精度。Ling和Templeton等人[35]通過深度神經(jīng)網(wǎng)絡(luò)方法計算了雷諾應(yīng)力的各向異性分量，并顯示了二次流中的旋渦結(jié)構(gòu)和波形壁中的分離現(xiàn)象。除了應(yīng)用于RANS模型，數(shù)據(jù)驅(qū)動方法也已用于模型化LES中的湍流相關(guān)變量。Gamahara等人[36]針對亞格子應(yīng)力張量的分量分別建立了人工神經(jīng)網(wǎng)絡(luò)模型。Maulik等人[37]采用人工網(wǎng)絡(luò)模型預(yù)測了時空變化的湍流源項(xiàng)。Wang等人[38]用不同機(jī)器學(xué)習(xí)方法和特征研究了亞格子應(yīng)力封閉。這些研究工作僅以高分辨率的數(shù)據(jù)作為驅(qū)動，一定程度上降低了模型封閉或湍流相關(guān)變量模型化的難度，證實(shí)了純數(shù)據(jù)驅(qū)動的黑箱模型在湍流研究應(yīng)用中的可行性。另一方面，一些研究者還將機(jī)器學(xué)習(xí)用于描述和量化傳統(tǒng)模型計算結(jié)果的不確定度，對未來的建模工作具有很好的指導(dǎo)作用。目前主要的研究方向和建模流程可大致歸結(jié)為圖1所示。

圖1 機(jī)器學(xué)習(xí)應(yīng)用于湍流研究的主要研究方向及流程Fig.1 Main research direction and process of machine learning applied to turbulence study

機(jī)器學(xué)習(xí)在湍流研究中的快速發(fā)展在一定程度上也得益于諸多廣泛使用的開源平臺。目前，TensorFlow、Keras、Theano以及Matlab等平臺為開展學(xué)習(xí)工作提供了有力支持。這些平臺大都內(nèi)嵌了多種學(xué)習(xí)框架供用戶選用，如RBFNN、RF以及更復(fù)雜的深度神經(jīng)網(wǎng)絡(luò)(DNN)和卷積神經(jīng)網(wǎng)絡(luò)(CNN)等。同時，模型的超參數(shù)優(yōu)化可采用隨機(jī)梯度下降(SGD)、動量隨機(jī)梯度下降(SGDM),以及Adam算法等實(shí)現(xiàn)。此外，基于CPU或GPU的并行算法縮短了湍流大數(shù)據(jù)的模型優(yōu)化過程。這些開發(fā)平臺的使用大大提高了研究者的工作效率，便于迅速開展并推進(jìn)研究工作，而且，代碼的通用性也有助于研究者之間的相互交流和探討。

機(jī)器學(xué)習(xí)與湍流建模相結(jié)合的研究工作是流體力學(xué)領(lǐng)域新興的研究方向?，F(xiàn)有的研究成果有力地驗(yàn)證了其可行性，預(yù)示了機(jī)器學(xué)習(xí)在未來湍流模型應(yīng)用中的積極前景[39]。與此同時，研究者也面臨著諸多問題與挑戰(zhàn)。一方面，現(xiàn)有的研究主要是基于DNS和LES計算的高分辨率數(shù)據(jù)，因此，選取的算例多是針對簡單幾何體的中低雷諾數(shù)流動問題，如平板和管流等。然而，工程中的高雷諾數(shù)湍流，如機(jī)翼和翼型繞流，則難以獲取高分辨率數(shù)據(jù)，如何將機(jī)器學(xué)習(xí)應(yīng)用于這一領(lǐng)域的湍流計算值得進(jìn)一步探索。另一方面，模型的泛化能力和魯棒性以及穩(wěn)定性需要進(jìn)一步提高?，F(xiàn)有研究中，當(dāng)預(yù)測數(shù)據(jù)與訓(xùn)練數(shù)據(jù)差別較大時或者幾何外形發(fā)生變化時，模型的表現(xiàn)會不同程度地下降。這在一定程度上是數(shù)據(jù)驅(qū)動方法本身固有的缺陷，另外，模型構(gòu)建過程中的一些方面，如數(shù)據(jù)處理與特征構(gòu)建、模型對象的選取以及約束性條件的施加也會影響模型的最終性能。模型與求解器的耦合需要額外地考慮模型的穩(wěn)定性。一種耦合形式我們稱之為“松耦合”，其過程是先通過RANS基模型與N-S求解器求解出初始流場和雷諾應(yīng)力場，然后，所構(gòu)建的機(jī)器學(xué)習(xí)湍流模型根據(jù)已獲得的初始流場構(gòu)建模型輸入并展開預(yù)測。通過模型的預(yù)測值更新初始的雷諾應(yīng)力，將改善的雷諾應(yīng)力直接傳遞給求解器。另一種則與傳統(tǒng)RANS模型相同，所構(gòu)建的模型在從初始流場開始的每一迭代步都與求解器之間互相反饋，直至N-S方程求解器獲得收斂解，我們稱之為“緊耦合”，如圖2所示。其中，穩(wěn)定性問題在后者中表現(xiàn)得更明顯。

(a)松耦合

(b)緊耦合圖2 機(jī)器學(xué)習(xí)湍流模型與CFD求解器的耦合Fig.2 Coupling of machine learning turbulence model with CFD solver

1 機(jī)器學(xué)習(xí)建模過程

基于機(jī)器學(xué)習(xí)方法的模型構(gòu)建過程主要包含數(shù)據(jù)處理、特征選取以及模型框架的確定及參數(shù)優(yōu)化等幾個方面。針對每一方面，研究者又可以采取不同的方法，例如模型框架可選擇網(wǎng)絡(luò)模型或樹模型，激活函數(shù)可選擇tanh或ReLu等等。這些方面在不同程度上影響模型的性能。Ling等人[40]比較了不同的分類器在預(yù)測湍流模型不確定度中的應(yīng)用，Zhang等人[41]比較了不同輸入對結(jié)果的影響,等等。模型框架和參數(shù)優(yōu)化方法的種類繁多，在分類和回歸以及各自對應(yīng)的特定問題中通用性較差，難以界定某種模型或方法的具體優(yōu)劣性。但是無論采用什么方法，特征選擇和數(shù)據(jù)處理都是研究者普遍面臨的問題。

1.1 特征選擇

特征選擇主要是為了選出相關(guān)的和有用的特征，降低特征的冗余度，另外還有減少特征集，提高模型性能和數(shù)據(jù)理解等的作用。在實(shí)際應(yīng)用中，過多的特征不僅會減緩學(xué)習(xí)過程甚至?xí)蓴_學(xué)習(xí)算法。需要注意的是，特征選擇對不同機(jī)器學(xué)習(xí)算法的效果不同，有時反而可能會降低模型性能[40]。目前，常用的特征選擇算法有濾過法，打包法和嵌入法[42]，但在機(jī)器學(xué)習(xí)湍流研究中，這些算法的應(yīng)用還較少。湍流的物理機(jī)理復(fù)雜，特征選擇的難度較大，甚至被假定為最大的挑戰(zhàn)[30]。研究者主要是根據(jù)物理知識或已有的控制方程和表達(dá)式進(jìn)行經(jīng)驗(yàn)性的選取和比較。由于不同研究的建模對象不同，研究者根據(jù)自身經(jīng)驗(yàn)進(jìn)行的特征選取往往通用性差且存在一定的偶然性。Ling等人[43]提出具有不變量屬性的輸入特征，并與不具備該屬性的特征進(jìn)行了對比。Wang等人[34]則進(jìn)一步給出了特征的相對重要性。另外，由于不同的輸出與流場變量的相關(guān)性不同，模型的輸入特征也可予以區(qū)別選取[38]。特征選取過程中亦應(yīng)考慮特征的數(shù)據(jù)分布范圍，盡量保證每一特征的數(shù)據(jù)分布較為緊湊，避免輸入特征和輸出中出現(xiàn)極端異值。除此之外，為了避免不同特征間的量級差異，研究者可以對特征進(jìn)行線性歸一化或標(biāo)準(zhǔn)化等處理。

1.2 輸入輸出數(shù)據(jù)處理

訓(xùn)練數(shù)據(jù)可以從湍流數(shù)值計算結(jié)果中直接獲得，也可以根據(jù)具體的建模需要間接求得。作為模型的定義域，訓(xùn)練數(shù)據(jù)所構(gòu)成的數(shù)據(jù)空間從根本上決定了模型的性能，其重要性不可忽視。如果訓(xùn)練數(shù)據(jù)量不大，建模者可以將全部數(shù)據(jù)用來訓(xùn)練模型。但是對于網(wǎng)格量較大的DNS、LES以及高雷諾數(shù)算例，可利用的訓(xùn)練數(shù)據(jù)往往過多而冗余。如果不加以篩選，過多的樣本會導(dǎo)致計算代價大而且訓(xùn)練過程周期長。為此，研究者大多在時間上采用快照，在空間上采用取截面的方法。這種方法雖然簡單高效，但是沒有擺脫網(wǎng)格劃分的影響，也不能直接體現(xiàn)出對流場關(guān)鍵區(qū)域的側(cè)重。另一種方法是根據(jù)訓(xùn)練數(shù)據(jù)間的相似性或歐式距離來剔除冗余的樣本。這種方法彌補(bǔ)了前者的缺點(diǎn)，但是篩選過程會隨著訓(xùn)練數(shù)據(jù)的增多而變得效率低下。研究者可依據(jù)具體問題結(jié)合著使用這兩種方法。對于模型的輸出，當(dāng)預(yù)測算例中出現(xiàn)超出模型泛化能力的流動特性時，模型的某些輸出值可能不滿足物理規(guī)律的約束或者屬于明顯錯誤的異值。這些值降低了模型的精度和魯棒性，因此，可以適當(dāng)?shù)剡M(jìn)行后處理。例如，為了確定數(shù)值穩(wěn)定性，Maulik等人[37]對模型輸出作如下截斷，

(1)

Ling等人[35]對模型預(yù)測的各向異性值進(jìn)行如下約束，

(2)

1.3 模型框架

機(jī)器學(xué)習(xí)方法的多樣以及相對應(yīng)的超參數(shù)極大地豐富了模型框架的選擇性。主要的模型框架及其在湍流中的不同應(yīng)用如表1所示。這些模型框架都能夠用于構(gòu)建復(fù)雜的非線性映射關(guān)系，研究者可根據(jù)具體問題具體選擇。模型的超參數(shù)可通過不斷調(diào)整和對比結(jié)果來確定，也可借鑒于某些超參數(shù)優(yōu)化算法，如Spearmint、Hyperopt等等。在超參數(shù)的優(yōu)化過程中，除了考慮模型在訓(xùn)練集上的精度外，亦應(yīng)考慮其在驗(yàn)證集或后驗(yàn)測試中的表現(xiàn)。

表1 主要的模型框架及應(yīng)用Table 1 Main model framework and appliance

2 機(jī)器學(xué)習(xí)在湍流模型中的應(yīng)用

經(jīng)過近幾年的努力，機(jī)器學(xué)習(xí)方法在湍流研究中的應(yīng)用得以迅速發(fā)展。在應(yīng)用范圍上，目前已開展了平板、管流及翼型等幾何體的剪切湍流以及分離流、二次流等復(fù)雜流動現(xiàn)象的建模工作。在應(yīng)用程度上，機(jī)器學(xué)習(xí)不僅用來構(gòu)建流場變量之間的復(fù)雜映射關(guān)系，還直接應(yīng)用于雷諾應(yīng)力的計算。特別地，根據(jù)模型作用可將現(xiàn)有的研究工作大致分為三類：模型的不確定度分析，改善湍流模型，以及替代湍流模型。

2.1 湍流模型的不確定度分析

與真實(shí)的物理解或高分辨率結(jié)構(gòu)相比，模型的計算結(jié)果往往存在一定的不確定度，其主要來源于NS方程的系綜平均計算、模型封閉函數(shù)形式、雷諾應(yīng)力表達(dá)式和模型經(jīng)驗(yàn)參數(shù)等。在具體應(yīng)用中，由于模型之間的封閉形式和經(jīng)驗(yàn)參數(shù)等不同，對不同流動的不確定度也難以把握。對于模型參數(shù)和模型形式導(dǎo)致的不確定度，可以采用靈敏度分析、概率分析以及貝葉斯方法等加以量化[47-53]。近年來，機(jī)器學(xué)習(xí)方法開始逐漸被應(yīng)用于模型的不確定度分析。構(gòu)建RANS模型參數(shù)與參數(shù)對應(yīng)的偏差之間的數(shù)據(jù)驅(qū)動模型，通過計算模型的最小值來確定RANS模型的最優(yōu)參數(shù)有利于提高精度[54]。模型研究者可以通過構(gòu)建分類器來預(yù)測流場中RANS模型的不確定區(qū)域，進(jìn)一步地還可以針對模型所確定的不確定區(qū)域采用更好的計算方法[40,55-56]。Singh等人[57]結(jié)合流場反演和貝葉斯方法來量化RANS模型的不確定度。更多關(guān)于RANS模型或LES的不確定度分析可參考文獻(xiàn)綜述[58]。

2.2 改善湍流模型

RANS模型由于存在結(jié)構(gòu)和參數(shù)上的不確定度，計算精度及適用范圍會受到不同程度的限制。對于某些存在激波-邊界層干擾的流動或分離流，RANS模型的計算結(jié)果往往存在一定的偏差甚至不準(zhǔn)確，因此，研究者希望根據(jù)高分辨率數(shù)據(jù)來減小RANS模型計算的偏差，或者使之能夠用于分離流的計算。這種采用機(jī)器學(xué)習(xí)方法來改善RANS模型的研究思路大致可分為兩類：一種是通過改變模型的控制方程形式，如乘以修正系數(shù)或給方程增加源項(xiàng)；另一種是在RANS模型基礎(chǔ)上構(gòu)造偏差函數(shù)，然后將RANS模型和偏差函數(shù)的計算結(jié)果疊加，作為最終的雷諾應(yīng)力值。

前者主要是針對基于渦黏假設(shè)的一方程或兩方程RANS模型。研究者首先根據(jù)實(shí)驗(yàn)結(jié)果或高分辨率數(shù)據(jù)反演計算出引入的修正系數(shù)或增加的源項(xiàng)，然后以此作為輸出，構(gòu)建出數(shù)據(jù)驅(qū)動的模型并將其與求解器耦合，流程圖如圖3所示。譬如，Singh等人[44]通過修正SA模型中的生成項(xiàng)來改變控制方程形式，使得修正后的模型能夠用于分離流的計算并取得與實(shí)驗(yàn)更相符的計算結(jié)果。SA模型的原方程為：

(3)

引入空間變量β(x)作為生成項(xiàng)的修正系數(shù)，方程變?yōu)椋?/p>

(4)

圖3 數(shù)據(jù)驅(qū)動的RANS求解器的構(gòu)建過程[31]Fig.3 Process for building a data-driven RANS solver[31]

通過實(shí)驗(yàn)結(jié)果對RANS模型的結(jié)果進(jìn)行優(yōu)化，利用伴隨方法對場反演后得到β(x)，然后基于機(jī)器學(xué)習(xí)方法建立流場變量與β(x)之間的映射。所構(gòu)建的模型具有很好的遷移性,并保證了原模型的收斂性，在不同程度上改善了原模型的計算結(jié)果。其中，14°攻角下S809翼型算例的壓力對比如圖4所示。

相似的研究工作還分別圍繞k方程的生成項(xiàng)和轉(zhuǎn)捩模型的源項(xiàng)或匯項(xiàng)等展開[31]。

第二種研究思路不改變控制方程形式，而是對模型計算后的結(jié)果加以修正。

τij=τij,RANS+Δτij,Model

(5)

由于基于渦黏假設(shè)的RANS模型不能準(zhǔn)確地計算各向異性的雷諾應(yīng)力，因此對二次流和分離流等流動現(xiàn)象的計算偏差較大。Xiao等人[34,59-60]針對RANS模型計算結(jié)果和高分辨率數(shù)據(jù)之間的差異進(jìn)行建模以修正原模型的計算結(jié)果(圖5)。在采用數(shù)據(jù)驅(qū)動方法的基礎(chǔ)上，研究者強(qiáng)調(diào)了雷諾應(yīng)力的物理約束，提出了“基于物理的機(jī)器學(xué)習(xí)”概念。當(dāng)訓(xùn)練模型與預(yù)測模型相同時，模型能夠明顯改善RANS模型的計算結(jié)果，如圖6所示。但是，即便是基于相似的流動現(xiàn)象(分離流)，對于幾何模型發(fā)生變化的區(qū)域(周期山的收縮段)，修正模型甚至惡化了原RANS模型的計算結(jié)果，如圖7所示。

(a)反演求得的β場

(b)神經(jīng)網(wǎng)絡(luò)預(yù)測的β場

(c)不同計算結(jié)果的對比

圖4 S809 翼型在α=14°和Re=2×106時的模型預(yù)測結(jié)果對比[44]
Fig.4 Comparison of inverse and NN-augmented predictions (using data-set P)for S809 airfoil atα=14°andRe=2×106[44]

圖5 湍流預(yù)測建模的基于物理信息的機(jī)器學(xué)習(xí)(PIML)框架[59]Fig.5 Proposed physics-informed machine learning (PIML)framework for predictive turbulence modeling[59]

改善原湍流模型以提高其精度的建模思想首先肯定了RANS模型的借鑒意義，本質(zhì)上是給模型的控制方程增加了一個源項(xiàng)或修正項(xiàng)?；诖硕⒌哪Ｐ臀覀兎Q之為修正模型。這種修正模型確實(shí)在一定程度上提高了原模型的精度，但同時也增加了額外的計算量，降低了整體的計算效率。另一方面，修正模型是針對某個特定的RANS模型所構(gòu)建的，然而，不同的RANS模型存在不同的不確定度。因此，所構(gòu)建的模型一般只適于某個特定的RANS模型，這會在很大程度上限制修正模型的使用范圍。

圖6 湍流基模型,DNS和PIML預(yù)測的湍動能的輪廓對比圖(用周期山作訓(xùn)練數(shù)據(jù))[34]Fig.6 Profiles of turbulence kinetic energy corresponding to baseline,PIML-predicted and DNS results (flows over periodic hills as training database)[34]

圖7 湍流基模型,DNS和PIML預(yù)測的湍動能的輪廓對比圖(用波形管道和彎曲后向臺階作訓(xùn)練數(shù)據(jù))[34]Fig.7 Profiles of turbulence kinetic energy corresponding to baseline,PIML-predicted and DNS results (flows over wavy channels and curved backward facing step as training database)[34]

2.3 替代湍流模型

與構(gòu)建修正模型的思想不同，替代湍流模型的出發(fā)點(diǎn)是希望直接從數(shù)據(jù)中找出流場變量與湍流之間映射關(guān)系?；诖硕⒌哪Ｐ臀覀兎Q之為替代模型。Ling等人[35]基于Pope推導(dǎo)的基張量和不變量構(gòu)建了雷諾應(yīng)力各向異性的張量基神經(jīng)網(wǎng)絡(luò)模型(TBNN)，如圖8所示。其中，雷諾應(yīng)力各向異性張量b可由基張量T線性表達(dá)為[61]:

圖8 張量基神經(jīng)網(wǎng)絡(luò)架構(gòu)示意圖[35]Fig.8 Schematic of tensor basis neural network architectures[35]

(6)

基張量與不變量的表達(dá)式分別如下，

T(1)=S

T(2)=SR-RS

T(5)=RS2-S2R

T(7)=RSR2-R2SR

T(8)=SRS2-S2RS

T(10)=RS2R2-R2S2R

(7)

λ1=Tr(S2),λ2=Tr(R2),λ3=Tr(S3),

λ4=Tr(R2S),λ5=Tr(R2S2)

(8)

其中，I表示單位陣，Tr表示跡，S表示平均應(yīng)變率張量，R表示平均旋轉(zhuǎn)張量。

該模型能夠預(yù)測出管流中的角渦和波形壁面的流動分離(圖9)，這一點(diǎn)是渦黏模型無法實(shí)現(xiàn)的。然而，TBNN與求解器之間進(jìn)行迭代收斂的可行性還有待進(jìn)一步驗(yàn)證。

圖9 波形壁測試算例中的流向速度云圖(分離區(qū)用灰線描繪)[35]Fig.9 Contours of streamwise velocity normalized by bulk velocity in the wavy wall test case,zoomed into the near-wall region[35],(separated regions outlined in grey)

Zhu等人[46]探討了直接構(gòu)建純數(shù)據(jù)驅(qū)動的湍流黑箱代數(shù)模型,并實(shí)現(xiàn)了與求解器之間的耦合計算(圖10)。研究結(jié)果表明，基于NACA0012翼型的三個亞聲速狀態(tài)算例，所構(gòu)建的模型可以實(shí)現(xiàn)與SA模型相當(dāng)?shù)木群透叩挠嬎阈?，并對計算狀態(tài)(圖11)和幾何外形(圖12)具備一定的泛化能力，驗(yàn)證了替代模型的可行性。與修正模型不同，替代模型的性能對最終摩擦阻力的計算結(jié)果起決定作用，而且對求解器的收斂性也有重要影響。換言之，模型構(gòu)建過程對模型輸出的精度和光滑性以及模型的穩(wěn)定性都有很高要求。對于回歸問題，損失函數(shù)決定了量級大的值比量級小的值容易取得更高的擬合精度，然而，邊界層內(nèi)(尤其是高雷諾數(shù)時)的劇烈變化容易導(dǎo)致湍流變量的明顯的尺度效應(yīng)。因此，替代模型很難在整個的訓(xùn)練數(shù)據(jù)空間中都保證較高的精度。模型性能欠佳的“邊緣區(qū)域”容易出現(xiàn)異值，繼而影響模型輸出的光滑性。此外，將所構(gòu)建的模型嵌入求解器后，模型預(yù)測在迭代過程中會存在一定程度的不穩(wěn)定，殘差在下降過程中出現(xiàn)震蕩，如圖13所示。

純數(shù)據(jù)驅(qū)動的建模工作亦在二維及三維亞格子應(yīng)力的模型封閉中開展。Sarghini等人[62]用神經(jīng)網(wǎng)絡(luò)方法提高了LES中亞格子應(yīng)力的計算效率。針對時空變化的湍流源項(xiàng)Π，Maulik等人[37]構(gòu)建了人工神經(jīng)網(wǎng)絡(luò)的亞格子模型，與之相關(guān)的渦量-流函數(shù)控制方程為:

圖10 構(gòu)建學(xué)習(xí)器和代理器的流程圖[46](q和μt分別表示輸入特征和渦黏)Fig.10 Flow chart for building the learning machine and surrogate machine[46],(q and μt mean the input features and eddy viscosity,respectively)

圖11 訓(xùn)練和預(yù)測算例中SA (方框)和RBFNN (三角形)計算的力系數(shù)結(jié)果對比[46]Fig.11 Comparison of the force coefficient calculated by SA (square)and RBFNN (delta)for both training and predicting cases[46]

(a)NACA0014 翼型

(b)RAE2822翼型

圖12 預(yù)測算例P1、P2和P5的摩擦阻力系數(shù)分布對比[46](圖中數(shù)據(jù)表示SA/RBFNN計算的摩擦阻力系數(shù)CD,f)
Fig.12 Predictions for NACA0014 and RAE2822 airfoil at P1,P2 and P5 cases[46].(the data inside areCD,fvalues calculated by SA/RBFNN model)

(9)

除了將模型與Smagorinsky和Leith模型對比外，Maulik等還進(jìn)一步研究了不同輸入特征和鄰居網(wǎng)格信息以及后驗(yàn)信息對模型性能和模型架構(gòu)選擇的影響，指出超參數(shù)優(yōu)化過程中需要后驗(yàn)信息的耦合。針對槽道湍流，Gamahara等人[36]采用人工神經(jīng)網(wǎng)絡(luò)對亞格子應(yīng)力的每一分量分別建立模型并對比了不同輸入特征的結(jié)果，缺點(diǎn)是模型容易受到濾過尺度的限制。Wang等人[38]用隨機(jī)森林和人工神經(jīng)網(wǎng)絡(luò)基于不同輸入特征研究了大渦模擬中亞格子應(yīng)力的數(shù)據(jù)驅(qū)動封閉，并對構(gòu)建的模型與Smagorinsky模型(SM)和動態(tài)Smagorinsky模型(DSM)以及其他方法計算的結(jié)果進(jìn)行了比較，如圖14所示。

(a)動能

(b)耗散率圖14 不同計算方法之間可分辨的動能和耗散率的時間演化圖[38]Fig.14 Temporal evolution of resolved kinetic energy and dissipation by different simulations[38]

2.4 模型的魯棒性與穩(wěn)定性

模型訓(xùn)練過程中，建模者都希望實(shí)現(xiàn)更高的模型精度，然而，單純地追求精度的提高往往是以損失魯棒性為代價的[63]。由于存在過擬合的風(fēng)險，訓(xùn)練階段較低的損失函數(shù)對應(yīng)的模型框架或模型參數(shù)并不意味著較高的預(yù)測能力，因此，可以在損失函數(shù)中引入L1或L2范數(shù)約束，或者采用Dropout等方法。此外，模型選擇時也可以耦合后驗(yàn)分析或者進(jìn)行交叉驗(yàn)證[35,37]。如果需要將所構(gòu)建的模型與求解器耦合，則在模型構(gòu)建過程中必須要考慮到模型的穩(wěn)定性。單純的基于數(shù)據(jù)的訓(xùn)練與預(yù)測完全忽略了模型偏差對流場的影響。事實(shí)上，對于顯式地改變雷諾應(yīng)力或亞格子應(yīng)力的模型，任何微小的模型擾動都有可能被累積放大為速度場中明顯偏差，而對于改變RANS模型控制方程中某一項(xiàng)表達(dá)式的修正模型，穩(wěn)定性問題并不明顯[45,64]。Wu等人[64-65]提出了衡量模型條件的方法，將雷諾應(yīng)力分解為線性和非線性部分并分別建模以增強(qiáng)模型的穩(wěn)定性。針對這一方面的研究工作還有待進(jìn)一步的探索和完善。

3 總結(jié)與展望

將數(shù)據(jù)驅(qū)動方法應(yīng)用于湍流計算的研究工作在近幾年得到了飛速發(fā)展。一方面，機(jī)器學(xué)習(xí)新方法的涌現(xiàn)為研究者提供了更多的技術(shù)手段，有利于更復(fù)雜關(guān)系的實(shí)現(xiàn)。另一方面，所構(gòu)建的模型逐漸從湍流變量間的簡單函數(shù)映射關(guān)系演變?yōu)橹苯拥睦字Z應(yīng)力計算，在RANS和LES的湍流計算中起到越來越重要的作用?，F(xiàn)有的研究成果顯示了這一研究方向的廣闊前景和發(fā)展?jié)摿?，給研究者以開拓新局面的鼓舞和振奮。與此同時，研究過程中也同樣暴露出許多問題與挑戰(zhàn)。首先，由于模型是數(shù)據(jù)驅(qū)動的，因此，模型的性能在本質(zhì)上是由所選取的訓(xùn)練數(shù)據(jù)決定的，模型往往在那些與訓(xùn)練數(shù)據(jù)差異較大的預(yù)測數(shù)據(jù)中表現(xiàn)較差。這就要求研究者謹(jǐn)慎地選取模型構(gòu)建對象和方法以盡量提高模型的泛化能力。第二，受數(shù)據(jù)集的分布和損失函數(shù)等的影響，模型對整個數(shù)據(jù)集的性能也不盡相同。對于那些位于模型性能欠佳區(qū)域的數(shù)據(jù)，計算結(jié)果往往出現(xiàn)異值，繼而影響模型的光滑性和魯棒性。即便基于相同的訓(xùn)練數(shù)據(jù)，采用不同的機(jī)器學(xué)習(xí)方法所構(gòu)建的模型也會有不同的魯棒性。模型的異值會在一定程度上降低收斂性，如果出現(xiàn)在流場中的關(guān)鍵區(qū)域也會降低求解精度，所以應(yīng)尤其注意模型在流場關(guān)鍵區(qū)域的性能。第三，樣本特征的構(gòu)建和選取是影響模型性能的重要因素。湍流的復(fù)雜和無序造成了樣本構(gòu)建和選取的困難。研究者借助于特征選取算法的同時，也需要對湍流本身有深刻的物理見解，特別地，應(yīng)根據(jù)流場中的不同流動特性有針對性地構(gòu)建和選擇模型的輸入特征。目前，系統(tǒng)的針對湍流機(jī)器學(xué)習(xí)的特征選取的研究還相對較少。

在現(xiàn)有工作積累的經(jīng)驗(yàn)基礎(chǔ)上，機(jī)器學(xué)習(xí)在未來的湍流模型化中必將扮演著更加重要的角色。研究者在未來的研究工作中需要關(guān)注以下幾點(diǎn)：

1)湍流建模輸入特征的提取與優(yōu)化。單純地依靠海量數(shù)據(jù)作為直接輸入，通過簡單增加神經(jīng)網(wǎng)絡(luò)模型的深度和寬度來構(gòu)建復(fù)雜映射關(guān)系并不是一種好的建模策略。如何對輸入數(shù)據(jù)進(jìn)行特征提取，特別是結(jié)合前人在湍流理論中已經(jīng)取得的物理規(guī)律和經(jīng)驗(yàn)，如選擇不變量、利用各種標(biāo)度率、量綱分析方法等，精簡輸入信息、簡化網(wǎng)絡(luò)模型的維度和復(fù)雜性，從而實(shí)現(xiàn)模型的緊湊性和泛化能力的提升。當(dāng)然，也可以借鑒人工智能方法開展輸入信息的特征構(gòu)建與提取，如Featuretools等。

2)湍流模型與N-S方程耦合求解過程中的穩(wěn)定性與收斂性。隨著機(jī)器學(xué)習(xí)在湍流模型化中的作用越來越明顯，模型的輸出結(jié)果對速度場的影響也越大。一方面，機(jī)器學(xué)習(xí)方法所構(gòu)建的模型本身就有可能出現(xiàn)局部峰值和光滑性差的輸出值，降低計算精度。另一方面，模型訓(xùn)練大多是基于某些時刻或流場中某些位置的數(shù)據(jù)進(jìn)行的，模型難以包含流場在迭代過程中的全部信息。當(dāng)模型參與求解器的每一步迭代時，模型的擾動會導(dǎo)致N-S方程求解的流場擾動，繼而反饋為模型的輸入偏差。這種反饋會造成迭代過程的發(fā)散，即使對于單向耦合策略也有可能造成求解過程的失穩(wěn)。如果模型的穩(wěn)定性較差，特別是在非線性很強(qiáng)的高雷諾數(shù)及大分離等復(fù)雜流動中，模型的誤差隨著迭代過程的累積容易在速度場中被放大，導(dǎo)致N-S方程求解器的最終計算結(jié)果變差。特別是采用高維深度學(xué)習(xí)網(wǎng)絡(luò)，所帶來的局部極值對耦合過程的穩(wěn)定性和收斂性都帶來了不利影響。

3)N-S方程和湍流模型耦合求解過程中的多平衡態(tài)問題?；谑諗繝顟B(tài)的計算結(jié)果構(gòu)建的模型在耦合求解器后有時還會出現(xiàn)多態(tài)解的問題。由于N-S方程和湍流的非線性和復(fù)雜性，模型的計算結(jié)果也會受初始條件等的影響，將訓(xùn)練后的模型嵌入N-S求解器后，最終的計算結(jié)果有可能會收斂為其它平衡態(tài)的解。關(guān)于模型與求解器的耦合問題值得進(jìn)一步研究。