尹宇輝，李浩然，張宇飛，陳海昕

(清華大學(xué) 航天航空學(xué)院，北京 100084)

0 引 言

了解湍流結(jié)構(gòu)的形成和演化機(jī)理從而準(zhǔn)確預(yù)測(cè)湍流在現(xiàn)代工程問(wèn)題中具有重要意義。目前受限于計(jì)算機(jī)技術(shù)的發(fā)展，直接數(shù)值模擬（Direct Numerical Simulation, DNS）或大渦模擬（Large Eddy Simulation,LES）等高精度數(shù)值方法仍然難以用于解決實(shí)際工程問(wèn)題，雷諾平均N-S方程（RANS）仍然是氣動(dòng)設(shè)計(jì)中最常用的數(shù)值模擬方法。然而，許多復(fù)雜的流動(dòng)現(xiàn)象例如旋流、強(qiáng)壓力梯度流動(dòng)、流線大曲率問(wèn)題等，使用RANS方法無(wú)法給出符合工程設(shè)計(jì)精度要求的結(jié)果[1]。以飛機(jī)飛行中的流動(dòng)分離問(wèn)題為例，大迎角下失速、激波/邊界層干擾引起的分離以及機(jī)翼表面積冰后升阻力性能的惡化都對(duì)飛行安全有著決定性的影響。然而對(duì)于RANS方法，研究人員比較了幾種常用的湍流模型，發(fā)現(xiàn)它們都不能準(zhǔn)確預(yù)測(cè)翼型上表面由于強(qiáng)逆壓梯度而導(dǎo)致的流動(dòng)分離[2]。對(duì)于RANS方法的計(jì)算誤差，研究者一致認(rèn)為主要來(lái)自所預(yù)測(cè)的雷諾應(yīng)力的差異，即對(duì)于雷諾應(yīng)力的建模誤差。要改善預(yù)測(cè)效果，就要開發(fā)更準(zhǔn)確的雷諾應(yīng)力預(yù)測(cè)模型，即找到準(zhǔn)確的雷諾應(yīng)力與平均流動(dòng)之間的映射關(guān)系[3-4]。

最早出現(xiàn)并且至今依然經(jīng)常使用的雷諾應(yīng)力封閉方法是基于線性渦黏性假設(shè)的湍流模型（Linear Eddy Viscosity Model, LEVM），它假設(shè)雷諾應(yīng)力和平均流應(yīng)變率成正比，并基于經(jīng)驗(yàn)關(guān)聯(lián)式或求解控制方程來(lái)確定兩者的比例系數(shù)，即渦黏性?；贚EVM的湍流模型具有較好的計(jì)算魯棒性，計(jì)算復(fù)雜問(wèn)題時(shí)間成本較低，對(duì)于大部分邊界層附著流動(dòng)模擬效果較好，但是由于忽略了流線曲率的影響以及雷諾應(yīng)力的實(shí)際各向異性特性，LEVM無(wú)法預(yù)測(cè)復(fù)雜的流動(dòng)，例如方管流動(dòng)中的二次流[5]和周期山流動(dòng)中的分離區(qū)計(jì)算[6]。為了克服LEVM的局限性，研究者提出了非線性渦黏性模式和直接對(duì)雷諾應(yīng)力進(jìn)行建模的雷諾應(yīng)力模式。這些改進(jìn)的模型在某些復(fù)雜流動(dòng)中可以提高預(yù)測(cè)精度。然而，非線性渦黏模型中額外添加的項(xiàng)和新增系數(shù)可能隨流動(dòng)問(wèn)題變化且十分敏感，而雷諾應(yīng)力模式則提高了計(jì)算成本，并且計(jì)算的魯棒性較差[7]。綜上所述，目前仍未出現(xiàn)一種既能給出不同類型復(fù)雜流動(dòng)的精確結(jié)果，又能保持魯棒性和節(jié)省成本的湍流計(jì)算模型，而此問(wèn)題對(duì)于傳統(tǒng)湍流建模框架難以徹底解決。

幾十年的湍流研究中產(chǎn)生了大量的高置信度湍流數(shù)據(jù)，包括DNS和LES的數(shù)值計(jì)算結(jié)果以及實(shí)驗(yàn)觀測(cè)結(jié)果。與此同時(shí)，隨著高效統(tǒng)計(jì)推斷算法和機(jī)器學(xué)習(xí)技術(shù)的發(fā)展，應(yīng)用這些數(shù)據(jù)輔助湍流建模成為可能，這種框架被稱為數(shù)據(jù)驅(qū)動(dòng)建模，數(shù)據(jù)驅(qū)動(dòng)方法目前的發(fā)展思路可以歸納為以下三類[8-9]：1）基于統(tǒng)計(jì)方法對(duì)RANS方法中的湍流模型不確定度進(jìn)行評(píng)估。這種方法旨在對(duì)雷諾應(yīng)力張量中的不確定性[10-11]進(jìn)行量化，并通過(guò)統(tǒng)計(jì)推斷方法識(shí)別流場(chǎng)中不確定度高的區(qū)域[12]，從而有助于約束不確定性并明確RANS預(yù)測(cè)差異的具體來(lái)源。2）基于湍流數(shù)據(jù)推斷湍流模型方程中待標(biāo)定系數(shù)和生成、耗散項(xiàng)相對(duì)大小。這種方法選擇對(duì)于湍流計(jì)算有關(guān)鍵影響的關(guān)注量（Quantities of Interests, QoIs）作為推斷變量，可以是模型系數(shù)[13]或湍流模式控制方程[14-15]中的項(xiàng)（其初始值通常稱為先驗(yàn)結(jié)果）。然后將這些選定的變量視為全流動(dòng)分布的隨機(jī)場(chǎng)，并根據(jù)給定的高置信度數(shù)據(jù)（如表面壓力分布[15]或觀測(cè)位置處的速度值[16]）推斷出最大化所選目標(biāo)。如果選取的關(guān)注量合理，則使用修正后的關(guān)注量計(jì)算的流動(dòng)結(jié)果（稱為后驗(yàn)結(jié)果）將與目標(biāo)吻合程度較高。3）用機(jī)器學(xué)習(xí)對(duì)湍流模型的預(yù)測(cè)結(jié)果直接進(jìn)行修正或替代湍流模型進(jìn)行雷諾應(yīng)力的預(yù)測(cè)。這種方法的目的是建立平均流動(dòng)特征和雷諾應(yīng)力特征之間的映射。通過(guò)與大數(shù)據(jù)相結(jié)合，機(jī)器學(xué)習(xí)方法可以在不需要人為給出大量先驗(yàn)知識(shí)的前提下，提取多層次特征并給出回歸或分類結(jié)果，這種強(qiáng)大的特征提取和回歸能力使其適合于構(gòu)建湍流模式的映射。在上述三種方法中，前兩種方法不能產(chǎn)生預(yù)測(cè)模型，即無(wú)法用于改善訓(xùn)練和推斷過(guò)程中沒(méi)有出現(xiàn)過(guò)的算例數(shù)據(jù)集，而基于大數(shù)據(jù)訓(xùn)練的機(jī)器學(xué)習(xí)模型則可以被推廣到具有相同類型的流動(dòng)現(xiàn)象和流動(dòng)結(jié)構(gòu)的問(wèn)題之中。機(jī)器學(xué)習(xí)模型的這種泛化能力使其更具工程應(yīng)用價(jià)值，有助于降低進(jìn)行實(shí)驗(yàn)或DNS計(jì)算的成本。因此，本研究采取第三種研究思路，即使用機(jī)器學(xué)習(xí)的方法來(lái)改進(jìn)RANS模型的預(yù)測(cè)結(jié)果。

近年來(lái)基于機(jī)器學(xué)習(xí)的湍流建模研究得到了廣泛的關(guān)注。其中對(duì)于訓(xùn)練目標(biāo)量的選擇反映了研究者對(duì)湍流預(yù)測(cè)誤差的主要來(lái)源的認(rèn)識(shí)以及應(yīng)當(dāng)如何處理的看法的不同，據(jù)此可以對(duì)相關(guān)研究進(jìn)行分類。

第一種觀點(diǎn)是以某一原始湍流模型為基準(zhǔn)，通過(guò)修改模型方程中的某些項(xiàng)的相對(duì)大小來(lái)實(shí)現(xiàn)改進(jìn)。Tracey等[17]利用機(jī)器學(xué)習(xí)對(duì)S-A模型中的源項(xiàng)進(jìn)行重構(gòu)，驗(yàn)證了該方法的可行性。Zhu等[18]直接以S-A模型計(jì)算的渦黏性為目標(biāo)，建立人工神經(jīng)網(wǎng)絡(luò)。用網(wǎng)絡(luò)代替原模型，提高了計(jì)算效率，降低了對(duì)網(wǎng)格密度的依賴性。以上為重現(xiàn)基準(zhǔn)模型結(jié)果的工作，而如果目標(biāo)是改進(jìn)原始湍流模型在復(fù)雜問(wèn)題中的效果，則應(yīng)首先使用貝葉斯反演方法來(lái)推斷修改后的源項(xiàng)，然后使用機(jī)器學(xué)習(xí)方法進(jìn)行訓(xùn)練[15,19-20]。在這些研究中，研究者通過(guò)在湍流模型方程的生成項(xiàng)或耗散項(xiàng)上乘上一個(gè)全場(chǎng)變化的乘子實(shí)現(xiàn)對(duì)模型方程和平均方程的修正。由于這種方法通過(guò)輸運(yùn)方程來(lái)影響雷諾應(yīng)力，因此更容易達(dá)到修正雷諾應(yīng)力場(chǎng)的光滑性要求。然而，基準(zhǔn)湍流模型使用的渦黏性假設(shè)限制了雷諾應(yīng)力的各向異性，無(wú)論如何改變幅值，其雷諾應(yīng)力主方向永遠(yuǎn)和應(yīng)變率張量一致，因此目前的理論最優(yōu)結(jié)果與真實(shí)雷諾應(yīng)力仍然存在一定的差異。

第二種觀點(diǎn)認(rèn)為應(yīng)選擇整個(gè)雷諾應(yīng)力作為訓(xùn)練目標(biāo)。在Ling的工作[21]中，雷諾應(yīng)力被表示為由平均流的應(yīng)變率張量和旋轉(zhuǎn)率張量構(gòu)成的10個(gè)整基的線性組合。神經(jīng)網(wǎng)絡(luò)的輸出變量是組合的各項(xiàng)系數(shù)，然后與相應(yīng)的張量基相結(jié)合，組合出雷諾應(yīng)力各向異性張量。這種類似非線性渦黏模型的想法啟發(fā)了后來(lái)的一些研究[22]。除線性組合外，另一種對(duì)雷諾應(yīng)力的表示方法是特征值分解方法。由于雷諾應(yīng)力為對(duì)稱二階張量，特征值分解將給出一組實(shí)特征值和相互正交的特征向量，其中特征向量可以進(jìn)一步用歐拉角[23]或四元數(shù)[24]來(lái)描述。對(duì)真實(shí)雷諾應(yīng)力和RANS預(yù)測(cè)的雷諾應(yīng)力進(jìn)行分解，得到各特征的差量場(chǎng)，并以此作為機(jī)器學(xué)習(xí)的目標(biāo)。

綜合對(duì)比兩種方法，以某個(gè)基準(zhǔn)湍流模型結(jié)果為基礎(chǔ)進(jìn)行修正的方法待修正變量數(shù)量少，預(yù)測(cè)難度較低，但無(wú)法預(yù)測(cè)完整的雷諾應(yīng)力特征；直接對(duì)雷諾應(yīng)力進(jìn)行表示和學(xué)習(xí)理論上能夠最大程度再現(xiàn)真實(shí)結(jié)果，但待預(yù)測(cè)量變多，同時(shí)由于預(yù)測(cè)量直接對(duì)平均方程產(chǎn)生影響，也帶來(lái)了光滑性的問(wèn)題。

課題組此前基于周期山算例進(jìn)行了機(jī)器學(xué)習(xí)輔助湍流建模的應(yīng)用，改進(jìn)后的框架在預(yù)測(cè)結(jié)果光滑性和精度上有所提升[25]。本文將基于此前研究，研究將此框架應(yīng)用于積冰翼型的復(fù)雜繞流問(wèn)題這一類工程實(shí)際問(wèn)題里典型且關(guān)鍵的分離流動(dòng)中的可能性。研究主要關(guān)注在建立機(jī)器學(xué)習(xí)預(yù)測(cè)框架過(guò)程中的輸入量/輸出量特征選擇問(wèn)題，外流繞流問(wèn)題中數(shù)據(jù)分布不均勻帶來(lái)的模型預(yù)測(cè)困難。文中提出了輸入輸出特征的選擇準(zhǔn)則，同時(shí)針對(duì)翼型流場(chǎng)進(jìn)行了局部區(qū)域建模，驗(yàn)證了局部?jī)鼋Y(jié)的可行性，并將方法應(yīng)用于典型積冰翼型的流動(dòng)計(jì)算，提高了原始湍流模型的預(yù)測(cè)精度。

1 研究方法

1.1 機(jī)器學(xué)習(xí)輔助湍流建?？蚣?/h3>

首先對(duì)本文中所使用的機(jī)器學(xué)習(xí)輔助湍流建模的框架進(jìn)行闡述。三維可壓縮流動(dòng)的RANS方程組如式(1)所示：

其中，ρ、u、f、e、p、λ和μ分別為雷諾平均后的密度、速度、體積力、比總能、壓力、體積黏性和分子黏性。T代表總應(yīng)力張量，包括壓力、分子黏性應(yīng)力和由于雷諾平均而產(chǎn)生的雷諾應(yīng)力τ。在傳統(tǒng)湍流建模思路中，雷諾應(yīng)力是由湍流模型方程確定的。當(dāng)RANS方程和湍流模型方程均收斂時(shí)，求解過(guò)程結(jié)束。而機(jī)器學(xué)習(xí)在湍流建模中的作用是能夠在給定一系列輸入特征的情況下給出雷諾應(yīng)力的直接或間接預(yù)測(cè)結(jié)果，即提供一個(gè)有預(yù)測(cè)能力的回歸函數(shù)。要建立此回歸函數(shù)，需要基于大量的已存在的高置信度數(shù)據(jù)，這個(gè)過(guò)程稱為機(jī)器學(xué)習(xí)的訓(xùn)練過(guò)程。在訓(xùn)練完成之后，該函數(shù)則可以用于預(yù)測(cè)新的流動(dòng)過(guò)程，稱為測(cè)試過(guò)程，體現(xiàn)了機(jī)器學(xué)習(xí)模型在未知問(wèn)題上的泛化能力。預(yù)測(cè)目標(biāo)的選擇、輸入輸出特征和機(jī)器學(xué)習(xí)算法對(duì)預(yù)測(cè)性能都有決定性的影響。

預(yù)測(cè)目標(biāo)的選擇隨著學(xué)習(xí)真值類型的變化而變化。本文關(guān)注的應(yīng)用場(chǎng)景為高雷諾數(shù)下飛行器實(shí)際飛行過(guò)程中的典型分離流動(dòng)問(wèn)題，而此方面高置信度DNS和LES模擬結(jié)果較為缺乏，本文研究中所選擇的真值為課題組此前發(fā)展的SPF修正的三方程k-v2-ω（記為SPFk-v2-ω）湍流模式的計(jì)算結(jié)果[26]。該模型應(yīng)用于高雷諾數(shù)下典型結(jié)冰翼型的流動(dòng)預(yù)測(cè)結(jié)果和實(shí)驗(yàn)測(cè)量得到的力系數(shù)曲線、流場(chǎng)脈動(dòng)量和平均量測(cè)量結(jié)果符合較好[26]。該三方程模式依然基于線性渦黏性模式框架，因此以該模式迭代收斂結(jié)果的渦黏性場(chǎng)作為真值，記為。在得到真值之后，將預(yù)測(cè)目標(biāo)取為渦黏性真值場(chǎng)和基準(zhǔn)湍流模式計(jì)算收斂結(jié)果（記為）的差量場(chǎng)，為保證預(yù)測(cè)的渦黏性場(chǎng)恒為正值，差量場(chǎng)的形式采用先取對(duì)數(shù)再做差的形式，如式（2）所示：

選取差量場(chǎng)作為目標(biāo)能夠更好地利用基準(zhǔn)湍流模式對(duì)流動(dòng)的先驗(yàn)結(jié)果，提供更合理的初始猜測(cè)，使得回歸函數(shù)只關(guān)注差異顯著的區(qū)域，從而降低訓(xùn)練難度。

為使得機(jī)器學(xué)習(xí)輔助湍流建?？蚣芨贤牧饔?jì)算的合理性和可解釋性，對(duì)于機(jī)器學(xué)習(xí)中的預(yù)測(cè)目標(biāo)有以下兩點(diǎn)要求：

1）預(yù)測(cè)目標(biāo)應(yīng)和其對(duì)應(yīng)物理量具有同樣自由度；

2）預(yù)測(cè)目標(biāo)對(duì)于坐標(biāo)旋轉(zhuǎn)和平移應(yīng)具有不變性。

可以看出，本文選取的Δ ln(vt)自由度為1，和渦黏性的自由度一致；另由于預(yù)測(cè)目標(biāo)為標(biāo)量場(chǎng)，因此自動(dòng)滿足對(duì)于坐標(biāo)旋轉(zhuǎn)和平移的不變性。因此以上兩點(diǎn)要求得到滿足。

選取渦黏性的差量場(chǎng)作為目標(biāo)，其輸入量特征應(yīng)來(lái)自于基準(zhǔn)模式計(jì)算的結(jié)果，并且在代入平均方程計(jì)算時(shí)需要在迭代過(guò)程中保持凍結(jié)，直至平均方程殘差重新收斂，即意味著在第二輪重新迭代開始之前只調(diào)用一次機(jī)器學(xué)習(xí)模型進(jìn)行預(yù)測(cè)。

機(jī)器學(xué)習(xí)輔助的湍流建?？蚣苷砣缦?，并如圖1所示。

圖1 機(jī)器學(xué)習(xí)輔助湍流建?？蚣苁疽鈭DFig. 1 Schematic of the machine learning assisted turbulence modeling framework

1）在訓(xùn)練集和測(cè)試集的算例上生成計(jì)算網(wǎng)格，并且計(jì)算基準(zhǔn)湍流模式的結(jié)果；

2）從訓(xùn)練集的基準(zhǔn)模型計(jì)算結(jié)果平均值x中構(gòu)造出平均流輸入特征q(x)|train并提取出基準(zhǔn)渦黏性場(chǎng)；

3）獲得待學(xué)習(xí)目標(biāo)差量場(chǎng)Δln(vt)|train；

4）根據(jù)訓(xùn)練集數(shù)據(jù)使用隨機(jī)森林算法構(gòu)建映射函數(shù)f:qΔln(vt)；

5）將訓(xùn)練后的隨機(jī)森林模型應(yīng)用于測(cè)試集，基于測(cè)試集平均流輸入特征q(x)|test給出Δln(vt)|test，和測(cè)試集基準(zhǔn)模型給出的結(jié)合得到經(jīng)過(guò)機(jī)器學(xué)習(xí)方法修正后的渦黏性場(chǎng)；

6）將預(yù)測(cè)出的渦黏性場(chǎng)代入CFD求解器中并凍結(jié)，重新迭代RANS方程直至平均流重新收斂，得到修正后的平均場(chǎng)結(jié)果x′|test，求解結(jié)束。

1.2 基于張量分析和物理特征分析的輸入特征選擇

根據(jù)1.1節(jié)的描述，本文中所使用的機(jī)器學(xué)習(xí)框架輸入特征應(yīng)當(dāng)由基準(zhǔn)湍流模型計(jì)算得到的平均流特征而構(gòu)造出。而這些輸入特征還應(yīng)當(dāng)滿足以下三條基本原則：

1）完整性，即輸入特征集應(yīng)包括與雷諾應(yīng)力分布相關(guān)的所有可能信息；

2）緊湊性，即應(yīng)識(shí)別并刪除輸入特征集中的無(wú)效或冗余信息；

3）可實(shí)現(xiàn)性，即要求只要流動(dòng)結(jié)構(gòu)保持不變，輸入特征集的變量在不同情況下（例如不同的參考坐標(biāo)系或流向）應(yīng)保持一致。

滿足前兩個(gè)原則有助于提高機(jī)器學(xué)習(xí)模型的性能，使其更具外推性和準(zhǔn)確性；第三個(gè)原則保證了數(shù)據(jù)驅(qū)動(dòng)模型能夠具有和傳統(tǒng)湍流模型一樣的物理合理性。關(guān)于輸入量特征的選取已有許多研究[27-28]。本文提出了兩種分別基于張量分析和物理特征識(shí)別的選擇準(zhǔn)則，以往研究中出現(xiàn)的輸入量均可以被分類入這兩種角度，同時(shí)還可以在準(zhǔn)則的指導(dǎo)下構(gòu)造出新的有效特征并添加到輸入集中。

1.2.1 基于雷諾應(yīng)力表示的特征量選擇

首先考慮通過(guò)尋找雷諾應(yīng)力所依賴的張量或矢量，并構(gòu)造一個(gè)完整的張量基來(lái)表示雷諾應(yīng)力。非線性渦黏性假設(shè)就是一種典型的基于雷諾應(yīng)力表示方法的湍流模式，其將雷諾應(yīng)力假定為平均流應(yīng)變率S和旋轉(zhuǎn)率Ω的多項(xiàng)式組合。Pope[29]利用張量函數(shù)表示理論，推導(dǎo)出包含10個(gè)分量的整基，從而得到多項(xiàng)式的緊致形式。完整性基礎(chǔ)保證了由S和Ω構(gòu)成的每一個(gè)對(duì)稱二階偏張量都可以表示為這10個(gè)張量的線性組合。受此啟發(fā)，Wu等[24]在此基礎(chǔ)上進(jìn)一步增加壓力梯度Δp和湍動(dòng)能梯度Δk，以補(bǔ)充強(qiáng)壓力梯度和湍流強(qiáng)度非均勻分布的影響。為了方便表示，將這兩個(gè)梯度轉(zhuǎn)化為相應(yīng)的反對(duì)稱張量Ap和Ak。同樣根據(jù)張量表示定理得到由此這四個(gè)張量組成的47個(gè)分量的整基。由于整基的分量都是張量，不能直接用作輸入特征，進(jìn)一步選擇張量的第一不變量，即每個(gè)張量的跡作為輸入特征。針對(duì)本文中的翼型繞流問(wèn)題，首先從以上特征中去掉二維下為零以及絕對(duì)值相同的重復(fù)特征（此時(shí)剩余17個(gè)特征），另外由于結(jié)冰翼型網(wǎng)格較為復(fù)雜，難以實(shí)現(xiàn)較高的光滑性，因此進(jìn)一步去掉剩余特征中對(duì)于網(wǎng)格光滑性十分敏感的部分高階特征，最終剩余6個(gè)特征，如表1所示。表中符號(hào)說(shuō)明如下：，其中，表 示向量范數(shù)，表示Frobenius矩陣范數(shù)，變量ε是湍流耗散率，nsym和nanti表示在構(gòu)造整基時(shí)中使用的對(duì)稱/反對(duì)稱張量的數(shù)量。

表1 基于雷諾應(yīng)力表示選取的特征量Table 1 Input features selected based on the Reynolds stress representation

1.2.2 基于流動(dòng)辨識(shí)的特征量選擇

由于輸入特征完全基于基準(zhǔn)計(jì)算平均場(chǎng)進(jìn)行構(gòu)建，當(dāng)平均場(chǎng)真值和基準(zhǔn)值差距較大時(shí)難以給出合理表示結(jié)果，為進(jìn)一步補(bǔ)充輸入特征引入第二種思路，即通過(guò)尋找或構(gòu)造標(biāo)志性變量用以定位和定量描述基準(zhǔn)場(chǎng)和真值場(chǎng)之間較大的差異。與1.2.1節(jié)的表示方法相比，尋找識(shí)別差異的標(biāo)記并將其作為輸入更直接并且和預(yù)測(cè)目標(biāo)相關(guān)性更高。除了標(biāo)記功能外，所選特征還應(yīng)具有明確的物理意義，以保證在具有相似流結(jié)構(gòu)的流中具有相似的性能。一方面可以通過(guò)構(gòu)造能夠識(shí)別關(guān)鍵流動(dòng)特征的準(zhǔn)則量來(lái)實(shí)現(xiàn)，因?yàn)榛鶞?zhǔn)值和真值的差異往往出現(xiàn)在復(fù)雜的流動(dòng)結(jié)構(gòu)中。以所積冰型為明冰的翼型為例，冰后產(chǎn)生強(qiáng)烈分離形成流動(dòng)分離泡并在機(jī)翼上再附，基準(zhǔn)湍流模型預(yù)測(cè)的誤差主要來(lái)自分離泡內(nèi)剪切層、強(qiáng)逆壓力梯度和旋流的預(yù)測(cè)，通過(guò)識(shí)別以上流動(dòng)結(jié)構(gòu)就可以找到待修正區(qū)域，在傳統(tǒng)建模中也采用了這種思路。另一方面，可以構(gòu)造出湍流相關(guān)物理量的相對(duì)大小比率來(lái)反映出基準(zhǔn)模型對(duì)于關(guān)鍵流動(dòng)現(xiàn)象的辨識(shí)程度，將其加入輸入量就可以更好定位基準(zhǔn)湍流模式的誤差以及給出修正方向。綜合兩方面考慮，以及和上面相同的考慮特征光滑性對(duì)于網(wǎng)格光滑性的敏感程度之后，整理出如下所示的6個(gè)標(biāo)量特征。

1）流動(dòng)剪切層和旋渦流動(dòng)的識(shí)別函數(shù)[14]：f1=d2Ω/(vt+v)，其中d為到壁面的距離，vt為基準(zhǔn)湍流模式預(yù)測(cè)的渦黏性，Ω=；

2）邊界層識(shí)別函數(shù)（DDES使用的保護(hù)函數(shù)）[30]：f2=1-tanh(8rd)3，其中；

3）基準(zhǔn)模式計(jì)算湍動(dòng)能和平均流動(dòng)能之比：f3=

4）基于壁面距離的湍流雷諾數(shù)[31]：

5）湍流時(shí)間尺度和平均流時(shí)間尺度之比：f5=S k/ε；

6）基準(zhǔn)模式計(jì)算渦黏性和分子黏性之比：f6=vt/v 。

綜上所述，從平均流變量中構(gòu)造出了12個(gè)標(biāo)準(zhǔn)化的不變量作為機(jī)器學(xué)習(xí)框架的輸入特征。

1.3 針對(duì)翼型繞流問(wèn)題的局部建模和代入計(jì)算

此前大部分機(jī)器學(xué)習(xí)輔助湍流建?？蚣艿难芯慷家詢?nèi)流流動(dòng)作為研究對(duì)象，內(nèi)流和外流在關(guān)鍵流動(dòng)結(jié)構(gòu)上可能比較接近（例如周期山流動(dòng)分離和結(jié)冰翼型后流動(dòng)分離），但是在進(jìn)行外流流動(dòng)的模擬時(shí)為了保證遠(yuǎn)場(chǎng)邊界條件對(duì)繞流物體附近流場(chǎng)的影響足夠小，遠(yuǎn)場(chǎng)的特征長(zhǎng)度需要是繞流物體特征長(zhǎng)度的50～100倍。遠(yuǎn)離壁面的流動(dòng)數(shù)據(jù)即使網(wǎng)格變疏，但仍然占據(jù)流場(chǎng)數(shù)據(jù)的一大部分，并且此部分流動(dòng)數(shù)據(jù)由于沒(méi)有受到過(guò)多擾動(dòng)，大部分輸入特征保持一致沒(méi)有波動(dòng)，預(yù)測(cè)的雷諾應(yīng)力和真值區(qū)別也較小，因此對(duì)于機(jī)器學(xué)習(xí)模型的訓(xùn)練帶來(lái)了不利影響。為解決這一問(wèn)題，Zhu等[18]在全場(chǎng)重構(gòu)S-A模型時(shí)，使用了分區(qū)建模的方法，根據(jù)到壁面距離對(duì)整個(gè)流場(chǎng)進(jìn)行劃分，不同區(qū)域分別訓(xùn)練。在此問(wèn)題中，本文進(jìn)一步考慮只針對(duì)機(jī)翼附近的區(qū)域進(jìn)行訓(xùn)練，而其他區(qū)域依然使用基準(zhǔn)湍流模式進(jìn)行計(jì)算，區(qū)分于原來(lái)的分區(qū)域建模方法，本文的思路可以稱為局部區(qū)域建模。和分區(qū)域建模對(duì)比，局部區(qū)域建模首先大大降低了機(jī)器學(xué)習(xí)模型訓(xùn)練的成本，因?yàn)檫h(yuǎn)場(chǎng)并不需要進(jìn)行訓(xùn)練；另外局部區(qū)域建模能夠合理限制機(jī)器學(xué)習(xí)算法的影響區(qū)域，機(jī)器學(xué)習(xí)模型影響的范圍基本被約束在壁面附近，對(duì)于復(fù)雜組合的流動(dòng)問(wèn)題，這種思路能夠保證機(jī)器學(xué)習(xí)的誤差不會(huì)影響到原本區(qū)別很小的區(qū)域造成額外的誤差。

接下來(lái)從模型訓(xùn)練和代入計(jì)算角度詳細(xì)說(shuō)明局部區(qū)域建模的思路。進(jìn)行局部區(qū)域建模從機(jī)器學(xué)習(xí)模型訓(xùn)練角度容易處理，本文根據(jù)預(yù)測(cè)目標(biāo)量（差量場(chǎng)）在空間中的顯著區(qū)域大小，確定了拆分的壁面距離范圍，并按照所生成翼型結(jié)構(gòu)網(wǎng)格的ijk指標(biāo)范圍進(jìn)行了拆分，在機(jī)器學(xué)習(xí)模型的訓(xùn)練和預(yù)測(cè)中，均只使用該區(qū)域內(nèi)的流動(dòng)數(shù)據(jù)。

對(duì)于局部區(qū)域建模，將預(yù)測(cè)模型和CFD平均流計(jì)算過(guò)程耦合起來(lái)是一個(gè)重要問(wèn)題。本文采用的基本思路是依然在全場(chǎng)計(jì)算湍流模型方程，當(dāng)每一迭代步需要將渦黏性代入平均流方程中時(shí)，在指定需要凍結(jié)的區(qū)域的渦黏性場(chǎng)替換為凍結(jié)的值，其余部分代入湍流模式計(jì)算的結(jié)果。這種方法能夠自由指定需要替換的區(qū)域而不必修改湍流模型方程的求解區(qū)域和相應(yīng)的邊界條件，容易在程序中實(shí)現(xiàn)；并且經(jīng)過(guò)計(jì)算驗(yàn)證，在收斂的情況下平均場(chǎng)能夠得到光滑的平均速度場(chǎng)。但是這種局部?jī)鼋Y(jié)的方法對(duì)CFD計(jì)算的收斂帶來(lái)了挑戰(zhàn)。由于凍結(jié)的渦黏性和基準(zhǔn)湍流模型計(jì)算的渦黏性存在較大差異，導(dǎo)致在分區(qū)的邊界上存在間斷，驗(yàn)證計(jì)算中表明如果直接進(jìn)行凍結(jié)會(huì)導(dǎo)致湍流模型方程在迭代過(guò)程中很快發(fā)散，并且無(wú)法通過(guò)調(diào)整CFL數(shù)或增加迭代步數(shù)來(lái)改善。

本文為解決局部區(qū)域凍結(jié)CFD計(jì)算容易發(fā)散的問(wèn)題，采用了混合的思路：在局部?jī)鼋Y(jié)區(qū)域?qū)⒆x入的機(jī)器學(xué)習(xí)模型預(yù)測(cè)的渦黏性和基準(zhǔn)湍流模式計(jì)算的渦黏性按照一定比例進(jìn)行混合，混合系數(shù)記為λhybrid，如式（3）所示，當(dāng)λhybrid為0時(shí)表示完全為凍結(jié)量，為1時(shí)則表示未代入全部為基準(zhǔn)計(jì)算結(jié)果。

在實(shí)際計(jì)算中，首先使用基準(zhǔn)湍流模式計(jì)算出待修正構(gòu)型的收斂結(jié)果，再采用如上的混合凍結(jié)方法計(jì)算若干步，之后再完全替換為凍結(jié)渦黏性，計(jì)算驗(yàn)證結(jié)果表明了這種處理能夠有效地解決凍結(jié)區(qū)域內(nèi)外初始差距過(guò)大帶來(lái)的湍流模型方程發(fā)散的問(wèn)題。

2 數(shù)值方法驗(yàn)證

本部分內(nèi)容將針對(duì)上面研究方法部分提到的全場(chǎng)和局部?jī)鼋Y(jié)渦黏性的計(jì)算方法進(jìn)行驗(yàn)證，在此過(guò)程中暫不引入機(jī)器學(xué)習(xí)模型，直接使用真值進(jìn)行計(jì)算，即獲得當(dāng)前方法能夠得到的理論上限，考察方法的合理性。所選取的算例采用了帶944冰型的GLC305翼型繞流問(wèn)題[32]（記為GLC305/944），模擬和實(shí)驗(yàn)中的來(lái)流馬赫數(shù)為0.12，基于來(lái)流速度和翼型弦長(zhǎng)的雷諾數(shù)為3.5×106，在本部分的驗(yàn)證工作中，選取迎角6°的工況作為研究對(duì)象，圖2給出了所生成的結(jié)構(gòu)網(wǎng)格，采用C型網(wǎng)格拓?fù)洌硇捅砻婢W(wǎng)格點(diǎn)282個(gè)（尾跡區(qū)網(wǎng)格點(diǎn)88個(gè)），垂直壁面的法向網(wǎng)格點(diǎn)97個(gè)，總計(jì)網(wǎng)格量為44 329；所選擇的需要進(jìn)行局部?jī)鼋Y(jié)的區(qū)域如紅色框所示，凍結(jié)范圍的確定目前根據(jù)對(duì)流動(dòng)狀態(tài)的觀察結(jié)果確定，選取包括速度梯度、分離泡范圍、渦黏性場(chǎng)變化等平均場(chǎng)特征明顯、梯度較大的區(qū)域，該區(qū)域內(nèi)總計(jì)網(wǎng)格量為20 952。本文使用適用于可壓縮流動(dòng)求解的開源程序CFL3D v6.7[33]求解流動(dòng)，并且在原程序的基礎(chǔ)之上添加了凍結(jié)全場(chǎng)/局部渦黏性以及進(jìn)行渦黏性混合的功能。

圖2 GLC305/944翼型計(jì)算網(wǎng)格和局部?jī)鼋Y(jié)區(qū)域示意Fig. 2 Schematic of the computation grid of GLC305/944 airfoil and the partial frozen area

首先使用SST模式和SPFk-v2-ω對(duì)該算例進(jìn)行模擬，前者的結(jié)果將作為后續(xù)機(jī)器學(xué)習(xí)訓(xùn)練和預(yù)測(cè)中的基準(zhǔn)結(jié)果，給出渦黏性的初始猜測(cè)以及組合出輸入特征。后者的渦黏性結(jié)果將作為機(jī)器學(xué)習(xí)訓(xùn)練的標(biāo)簽值，并在預(yù)測(cè)集上作為真值和預(yù)測(cè)結(jié)果對(duì)比。圖3（a、b）分別給出了SST和SPFk-v2-ω的平均速度計(jì)算結(jié)果。從圖中可以看出，SST模式計(jì)算得到的分離泡分布在整個(gè)翼型上表面，而SPFk-v2-ω的計(jì)算結(jié)果和文獻(xiàn)中的實(shí)驗(yàn)結(jié)果相近，分離泡在43%左右弦長(zhǎng)處再附。兩種模式計(jì)算結(jié)果的區(qū)別主要在于對(duì)冰角后的自由剪切層流動(dòng)模擬的差異。文獻(xiàn)[26]指出，SST模式由于基于湍流平衡假設(shè)，在自由剪切層中湍動(dòng)能生成項(xiàng)始終和破壞項(xiàng)賦值相近，從而無(wú)法模擬出自由剪切層發(fā)展過(guò)程中高低速流動(dòng)的充分混合，產(chǎn)生了更大的分離區(qū)，而SPFk-v2-ω模式通過(guò)設(shè)置自由剪切層標(biāo)識(shí)函數(shù)，對(duì)原始三方程源項(xiàng)進(jìn)行修改并設(shè)置閾值函數(shù)，最終得到了和實(shí)驗(yàn)更為符合的湍動(dòng)能分布和生成耗散特性。

得到兩種模式的計(jì)算結(jié)果之后，如1.3節(jié)所述進(jìn)行兩種凍結(jié)計(jì)算。首先將SPFk-v2-ω模式計(jì)算得到的整個(gè)計(jì)算域的渦黏性提取出來(lái)，并且代入SST模式的收斂計(jì)算結(jié)果進(jìn)行續(xù)算，并直至重新收斂；第二種凍結(jié)計(jì)算是將如圖2所示紅色區(qū)域內(nèi)的渦黏性進(jìn)行凍結(jié)，并且采用在前55 000步迭代中進(jìn)行混合（取λhybrid= 0.5），再全部?jī)鼋Y(jié)為讀入的渦黏性進(jìn)行重新計(jì)算直至流動(dòng)收斂，兩種迭代得到的平均速度云圖（以迎角6°下的結(jié)果為例）結(jié)果如圖3（c、d）所示，對(duì)應(yīng)的升力系數(shù)的收斂曲線如圖4所示。從收斂情況可以看出，局部?jī)鼋Y(jié)的開始階段力系數(shù)的震蕩較明顯，但是混合方法保證了計(jì)算能夠繼續(xù)，經(jīng)過(guò)一定步數(shù)迭代后替換為全凍結(jié)就可以保證流動(dòng)收斂到合理的結(jié)果。從圖3（b～d）對(duì)比可以看出，無(wú)論使用全場(chǎng)凍結(jié)還是局部?jī)鼋Y(jié)，最終得到的平均流場(chǎng)和SPFk-v2-ω的計(jì)算結(jié)果一致性很高，證明了當(dāng)前方法的有效性和機(jī)器學(xué)習(xí)模型能夠達(dá)到的理論上限足以實(shí)現(xiàn)對(duì)平均流動(dòng)的正確修正。

此外，本文還針對(duì)所提出的渦黏性混合方法進(jìn)行了進(jìn)一步計(jì)算驗(yàn)證，以說(shuō)明不同的混合系數(shù)對(duì)于計(jì)算的影響。同樣，迎角6°算例，取混合系數(shù)Ratio分別為0 （相當(dāng)于不進(jìn)行混合）、0.25、0.50、0.75進(jìn)行4次計(jì)算。所有計(jì)算均在使用基準(zhǔn)SST湍流模型計(jì)算40 000步結(jié)果收斂之后，在其基礎(chǔ)上進(jìn)行續(xù)算得到。4個(gè)計(jì)算結(jié)果的升力系數(shù)結(jié)果如圖5所示。當(dāng)混合系數(shù)取0時(shí)，可以看到流動(dòng)在達(dá)到收斂結(jié)果之前就計(jì)算發(fā)散；當(dāng)采用了混合系數(shù)之后，在進(jìn)行混合計(jì)算的40 000步至80 000步之間，計(jì)算均可收斂。而當(dāng)80 000步以后，此時(shí)所有計(jì)算都結(jié)束混合并凍結(jié)為雷諾應(yīng)力全量，混合系數(shù)為0.25和0.5的結(jié)果能夠最終順利收斂至同一結(jié)果，對(duì)應(yīng)的升力系數(shù)為真值的升力系數(shù)，而混合系數(shù)0.75則在此過(guò)程中出現(xiàn)了計(jì)算發(fā)散，說(shuō)明在混合計(jì)算階段沒(méi)有起到完全消除計(jì)算不穩(wěn)定性的作用。以上計(jì)算結(jié)果說(shuō)明，混合系數(shù)分別影響“原始計(jì)算-部分凍結(jié)”和“部分凍結(jié)-全部?jī)鼋Y(jié)”的兩個(gè)過(guò)程，混合系數(shù)越小，第一個(gè)過(guò)程的震蕩和發(fā)散的可能越大，反之則第二個(gè)過(guò)程震蕩和發(fā)散的可能性更大。只要確保流動(dòng)收斂之后，不同混合系數(shù)可以達(dá)到同樣的計(jì)算結(jié)果。

圖3 使用不同渦黏性計(jì)算方法得到的平均場(chǎng)結(jié)果Fig. 3 Mean flow results obtained with different viscosity computation methods

圖4 兩種凍結(jié)計(jì)算方式的升力系數(shù)收斂曲線Fig. 4 Lift coefficient convergence histories of two different frozen methods

圖5 不同混合系數(shù)對(duì)于計(jì)算結(jié)果的影響Fig. 5 Effect of different hybrid ratios on the computation results

3 模型訓(xùn)練與預(yù)測(cè)結(jié)果

本文通過(guò)改變來(lái)流迎角構(gòu)造出訓(xùn)練集和測(cè)試集，分別計(jì)算迎角為4°、5°、6°和7°四個(gè)工況，選取4°和6°的結(jié)果為訓(xùn)練集，5°和7°的結(jié)果為測(cè)試集。由于流動(dòng)分離現(xiàn)象的變化隨迎角增加具有一致性，因此迎角5°和7°相當(dāng)于機(jī)器學(xué)習(xí)預(yù)測(cè)的內(nèi)插和外推情況，可以更好地考察機(jī)器學(xué)習(xí)模型的泛化能力。

本文使用的機(jī)器學(xué)習(xí)算法為隨機(jī)森林方法，基于Python平臺(tái)的scikit-learn包進(jìn)行程序?qū)崿F(xiàn)[34]。隨機(jī)森林方法基于一組經(jīng)過(guò)訓(xùn)練的決策樹，能夠在進(jìn)行回歸或分類的同時(shí)給出每個(gè)輸入特征在所有決策樹中的重要性評(píng)價(jià)。模型超參數(shù)包括決策樹的最大特征數(shù)和決策樹數(shù)量，前者取為5，從而滿足1 + log2n的關(guān)系式，其中n= 12為輸入特征數(shù)量；后者取為300，和文獻(xiàn)[25]保持一致。

為和本文提出的局部區(qū)域建模凍結(jié)的計(jì)算方法進(jìn)行對(duì)比，建立兩個(gè)隨機(jī)森林模型，分別使用全場(chǎng)數(shù)據(jù)和局部數(shù)據(jù)進(jìn)行訓(xùn)練和預(yù)測(cè)。訓(xùn)練在工作站上進(jìn)行，CPU型號(hào)為Intel Xeon E5-2 670，主頻為2.30 GHz，內(nèi)存為64 GB。兩個(gè)模型訓(xùn)練過(guò)程的基本信息如表2所示。從對(duì)比可以看出，采用局部區(qū)域凍結(jié)能夠極大降低訓(xùn)練成本，對(duì)于大數(shù)據(jù)量的情況下優(yōu)勢(shì)十分顯著。

表2 兩個(gè)隨機(jī)森林模型訓(xùn)練過(guò)程基本信息Table 2 Basic information for the training process using two random forest models

機(jī)器學(xué)習(xí)的結(jié)果在作為訓(xùn)練集的迎角4°和6°上均和真值符合很好，為展現(xiàn)機(jī)器學(xué)習(xí)模型的泛化能力，下面集中于預(yù)測(cè)集上的結(jié)果，即迎角5°和7°。圖6以迎角5°為例給出渦黏性對(duì)數(shù)差量場(chǎng)的真值和預(yù)測(cè)值的云圖對(duì)比，其中圖6（a）為真值，圖6（b）為對(duì)全場(chǎng)渦黏性差量場(chǎng)進(jìn)行訓(xùn)練和預(yù)測(cè)的結(jié)果，圖6（c）則為針對(duì)局部渦黏性差量場(chǎng)進(jìn)行建模的結(jié)果，因此并沒(méi)有全場(chǎng)的結(jié)果。首先可以從渦黏性差量場(chǎng)的真值分析基準(zhǔn)SST模式和作為真值的SPFk-v2-ω模式計(jì)算結(jié)果在平均場(chǎng)上的差別：在流動(dòng)剪切層剛開始發(fā)展的位置降低了渦黏性，保證分離泡初始發(fā)展、高度等特征和真值更加符合，在翼型弦長(zhǎng)30%以后相較于基準(zhǔn)SST模式具有更高的渦黏性，保證流動(dòng)在正確的位置附著。對(duì)比基于全場(chǎng)和局部數(shù)據(jù)訓(xùn)練模型可以看出，局部訓(xùn)練數(shù)據(jù)將范圍約束在翼型附近之后，對(duì)于剪切層開始位置的擬合程度更高，而全場(chǎng)數(shù)據(jù)訓(xùn)練模型則沒(méi)有完全反映出冰后的渦黏性分布趨勢(shì)。

圖6 渦黏性對(duì)數(shù)差量場(chǎng)對(duì)比Fig. 6 Comparison of the eddy viscosity difference in the log scale

圖7 翼型上表面速度剖面對(duì)比Fig. 7 Comparison of the velocity profiles on the upper surface of the airfoil

用預(yù)測(cè)出的對(duì)數(shù)差量場(chǎng)結(jié)合SST模式的渦黏性場(chǎng)重構(gòu)出預(yù)測(cè)的渦黏性值，并且將渦黏性代入CFD求解器中進(jìn)行凍結(jié)，直至平均流重新收斂。其中基于局部數(shù)據(jù)訓(xùn)練的模型，其凍結(jié)計(jì)算按照數(shù)值方法部分的說(shuō)明進(jìn)行混合計(jì)算，混合步長(zhǎng)和系數(shù)均保持一致。圖7分別展示了迎角5°和7°兩個(gè)預(yù)測(cè)算例從流場(chǎng)中截取特定流向位置剖面的速度數(shù)據(jù)進(jìn)行對(duì)比，黑色實(shí)線、紅色實(shí)線和藍(lán)色虛線分別為原始SST模式、作為真值的SPFk-v2-ω模式以及使用局部機(jī)器學(xué)習(xí)模型預(yù)測(cè)的結(jié)果?？梢钥闯鰺o(wú)論是內(nèi)插還是外推結(jié)果，使用機(jī)器學(xué)習(xí)方法修正后的渦黏性代入計(jì)算的平均流效果能夠和真值結(jié)果符合，分離泡的大小和再附位置都和真值十分接近。

進(jìn)一步對(duì)比翼型表面的壓力系數(shù)結(jié)果，如圖8所示?？梢钥闯鰤毫ο禂?shù)也和真值符合程度很好。綜上所述，機(jī)器學(xué)習(xí)應(yīng)用于分離流動(dòng)的平均流預(yù)測(cè)具有較好的擬合能力和泛化能力，能夠切實(shí)改進(jìn)預(yù)測(cè)效果。

圖8 翼型壓力系數(shù)對(duì)比Fig. 8 Comparison of the pressure coefficients of the airfoil

4 結(jié) 論

本文討論了機(jī)器學(xué)習(xí)輔助湍流建?？蚣艿慕ⅲ◤妮斎肓考蠘?gòu)建、特征選擇和機(jī)器學(xué)習(xí)方法選擇。針對(duì)將機(jī)器學(xué)習(xí)模型應(yīng)用于外流流動(dòng)這一場(chǎng)景，采用了僅針對(duì)流場(chǎng)部分區(qū)域進(jìn)行建模的方式以及相應(yīng)提升計(jì)算收斂性的雷諾應(yīng)力混合計(jì)算方法。將這種方法和全場(chǎng)建模在時(shí)間成本和預(yù)測(cè)效果上進(jìn)行了對(duì)比，結(jié)果表明局部建模方法降低了機(jī)器學(xué)習(xí)框架應(yīng)用于復(fù)雜外流流動(dòng)問(wèn)題的時(shí)間成本，提高了訓(xùn)練效率和準(zhǔn)確性，僅針對(duì)關(guān)鍵區(qū)域進(jìn)行修正同樣能夠?qū)崿F(xiàn)對(duì)于平均流動(dòng)較好的修正效果。進(jìn)一步工作將探討以LES或DNS結(jié)果為真值情況下局部建模方法的效果。