鄭天韻 王圣業(yè)? 王光學(xué) 鄧小剛

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

2) (中山大學(xué)物理學(xué)院, 廣州 510275)

1 引 言

邊界層轉(zhuǎn)捩普遍存在于流體機(jī)械內(nèi)部或其表面, 是經(jīng)典物理學(xué)中亟待解決的挑戰(zhàn)性問題. 在飛行器的運行過程中, 流體受到包括來流湍流度、物面粗糙度、馬赫數(shù)等諸多因素的影響, 易由分層穩(wěn)定流動向混沌湍流轉(zhuǎn)變. 這一轉(zhuǎn)捩過程同時伴隨著壁面摩阻和傳熱特性的急劇變化[1]. 因此, 對轉(zhuǎn)捩流動進(jìn)行快速準(zhǔn)確的計算流體力學(xué)(computational fluid dynamics, CFD)模擬對飛行器性能和設(shè)計效率的提升有重要意義.

轉(zhuǎn)捩預(yù)測可以采用求解全尺度湍流脈動的直接數(shù)值模擬 (direct numerical simulation, DNS)或僅求解大尺度脈動的大渦模擬(large-eddy simulation, LES)方法, 但其計算量隨雷諾數(shù)Re呈指數(shù)型增長(相當(dāng)于Re9/4), 當(dāng)前計算機(jī)技術(shù)的發(fā)展仍舊難以滿足其工程應(yīng)用的計算需求[2]. 結(jié)合湍流模型的雷諾平均數(shù)值模擬(Reynolds averaged Navier-Stokes, RANS)方法憑借其易用性及高效性, 在工程實踐中仍占有舉足輕重的地位[3,4].

轉(zhuǎn)捩流動中, 流場在同一空間位置會間歇性地呈現(xiàn)層流或湍流狀態(tài), 稱為間歇現(xiàn)象. 若采用函數(shù)I(x,y,z,t)描述這一現(xiàn)象, 并定義層流時函數(shù)值為0, 湍流時為1, 那么間歇因子g即為該函數(shù)的時間平均值[5],

考慮間歇性的轉(zhuǎn)捩模式通過引入“間歇因子”定量描述湍流生成, 是目前工程轉(zhuǎn)捩預(yù)測中最流行的一類方法. 模式通過顯式方程或求解額外的輸運方程得到所需的間歇因子. 早期, Dhawan等[6]采用經(jīng)驗公式描述間歇因子的流向分布, 但得到的間歇因子僅取決于來流及物面條件, 不考慮流場結(jié)構(gòu), 僅適用于簡單流動. Libby[7]首先引入帶有間歇因子的輸運方程計算湍流, 其后許多研究工作聚焦于改善輸運方程對多種流動例如自由剪切流[8]、邊界層[9]轉(zhuǎn)捩的預(yù)測效果. 然而, 這些模式往往通過非局部變量判斷轉(zhuǎn)捩起始, 難以在并行計算和非結(jié)構(gòu)網(wǎng)格中得到應(yīng)用. Menter等[10?12]采用應(yīng)變率為底的雷諾數(shù)代替動量厚度雷諾數(shù)Req, 發(fā)展了完全基于局部變量的四方程SST-g-Req轉(zhuǎn)捩模型. 其中無量綱的g,Req輸運方程的守恒形式為

其中,Pg,Eg分別為g方程的生成項和耗散項,Pθt為方程的源項. 上述輸運方程僅用于獲取間歇因子, 還需耦合SST湍流模型以模擬轉(zhuǎn)捩過程. 國內(nèi)外許多研究[13?15]表明該模型對轉(zhuǎn)捩流動有較好的預(yù)測能力. Bas等[16,17]提出了一種基于當(dāng)?shù)仄骄康拇鷶?shù)轉(zhuǎn)捩模型(或稱為B-C模型), 提高了計算效率, 但不適用于湍流度在(0.5, 2.0)區(qū)間內(nèi)的情況[17]. 另外, 傳統(tǒng)轉(zhuǎn)捩模型受經(jīng)驗公式束縛, 對計算平臺敏感, 往往不具備很好的通用性.

要解決上述問題, 關(guān)鍵在于找到流場當(dāng)?shù)仄骄颗c間歇因子之間的映射關(guān)系, 而機(jī)器學(xué)習(xí)方法則提供了一種新的途徑. He等[18]于2015年提出殘差網(wǎng)絡(luò) (deep residual network, ResNet), 這種新的結(jié)構(gòu)在前饋神經(jīng)網(wǎng)絡(luò)的基礎(chǔ)上引入跨層連接, 解決了深層網(wǎng)絡(luò)的訓(xùn)練問題, 受到廣泛應(yīng)用.ResNet憑借其復(fù)雜的網(wǎng)絡(luò)結(jié)構(gòu), 具備從大量數(shù)據(jù)中描述復(fù)雜非線性映射的能力且易于訓(xùn)練, 已經(jīng)在各領(lǐng)域取得了舉世矚目的成就.

有研究已經(jīng)將神經(jīng)網(wǎng)絡(luò)這類機(jī)器學(xué)習(xí)技術(shù)應(yīng)用于CFD領(lǐng)域. 許多研究者將機(jī)器學(xué)習(xí)方法應(yīng)用于重構(gòu)或修正渦黏系數(shù)或雷諾應(yīng)力. Ling等[19,20]構(gòu)建了張量基神經(jīng)網(wǎng)絡(luò)模型(TBNN) 以預(yù)測雷諾應(yīng)力各向異性張量, 改善了對管流角渦和波形壁流動分離的預(yù)測, 但其與求解器間的迭代收斂性還有待驗證. Zhu等[21]直接構(gòu)建了湍流渦黏性的純數(shù)據(jù)驅(qū)動“黑箱”代數(shù)模型, 但其對標(biāo)SA湍流模型,不具備對轉(zhuǎn)捩的預(yù)測能力. Duraisamy等[22,23]基于Ge等[24]的k-w-g轉(zhuǎn)捩模型, 利用神經(jīng)網(wǎng)絡(luò)(NN)和高斯過程(GP)重構(gòu)了g輸運方程中的源項, 改善了模型對于T3系列平板算例的預(yù)測精度. 該研究證明了機(jī)器學(xué)習(xí)方法在轉(zhuǎn)捩模型構(gòu)建中的可行性, 但沒有進(jìn)一步探究模型對于新幾何外形和來流條件的泛化能力.

與以往研究不同, 本文基于SST-g-Req模型的計算數(shù)據(jù), 利用深度殘差網(wǎng)絡(luò)ResNet重構(gòu)了當(dāng)?shù)仄骄颗c間歇因子g之間的純數(shù)據(jù)驅(qū)動“黑箱”映射模型, 并與SA模型相耦合, 發(fā)展了一種數(shù)據(jù)驅(qū)動的轉(zhuǎn)捩模型. 模型需要求解的微分方程數(shù)量僅為一, 相比于四方程的SST-g-Req模型具有更高的計算效率. 采用滿足伽利略不變性的當(dāng)?shù)亓鲌銎骄卣髁孔鳛檩斎胩卣? 便于與現(xiàn)代CFD求解器耦合,同時保證良好的泛化性能.

在相同網(wǎng)格條件下, 高精度格式相比于傳統(tǒng)二階格式能夠獲得更準(zhǔn)的氣動力結(jié)果[25,26]. 本文中將基于高階精度的加權(quán)緊致非線性格式(weighted compact nonlinear scheme, WCNS)[27,28]對 RANS方程和模型方程進(jìn)行數(shù)值離散.

2 模型構(gòu)造方法

2.1 整體思路

本文結(jié)合深度殘差網(wǎng)絡(luò)技術(shù)與SA全湍流模型[29], 參考基于局部變量的間歇轉(zhuǎn)捩模型構(gòu)造思路, 發(fā)展了一種數(shù)據(jù)驅(qū)動的轉(zhuǎn)捩預(yù)測方法. 模型的構(gòu)建過程可分為兩個部分: 數(shù)據(jù)學(xué)習(xí)和耦合求解.數(shù)據(jù)學(xué)習(xí)部分主要包括訓(xùn)練數(shù)據(jù)選擇, 神經(jīng)網(wǎng)絡(luò)模型框架和參數(shù)優(yōu)化. 在耦合求解部分, 首先將流場中間歇因子g的初值設(shè)為1(即采用全湍流模型),利用二階格式快速得到初始收斂場; 接著將流場的平均特征量 (q1,q2, …,qn) 代入神經(jīng)網(wǎng)絡(luò)模型得到g的預(yù)測結(jié)果; 最后利用g修正SA模型, 采用WCNS高精度格式計算收斂, 使其具有預(yù)期的轉(zhuǎn)捩預(yù)測效果. 具體過程如圖1所示.

圖 1 整體框架Fig. 1. Overall framework.

2.2 訓(xùn)練數(shù)據(jù)集的構(gòu)建

本文選擇間歇因子g作為殘差網(wǎng)絡(luò)的預(yù)測對象. 相較于直接預(yù)測雷諾應(yīng)力或渦黏性, 這樣做的好處在于避免神經(jīng)網(wǎng)絡(luò)這類模糊預(yù)測方法產(chǎn)生的結(jié)果直接作為最終結(jié)果, 降低極端異值和非物理解的影響, 結(jié)果更具可信度. 然而g并非流場中實際存在的物理量, 無法從DNS, LES等高精度數(shù)據(jù)中直接獲得其真實分布. 即便通過貝葉斯反演等方法利用高精度數(shù)據(jù)能夠反算出一組間歇因子g, 這種結(jié)果導(dǎo)向型的方法不能得到唯一確定的解[30]. 消耗大量計算代價得到的g場并不一定反映了真實的流動狀態(tài). Zhu等[21]基于SA模型的計算數(shù)據(jù)訓(xùn)練 NN 模型, 取得了不錯的效果. 鑒于此, 可采用更易獲得的g-Req轉(zhuǎn)捩模型得到的間歇因子g分布作為真值, 同時訓(xùn)練數(shù)據(jù)集中還包含了二階SA湍流模型的結(jié)果作為輸入特征.

本文主要針對自由流湍流強(qiáng)度遠(yuǎn)小于1.0%的自然轉(zhuǎn)捩. 傳統(tǒng)轉(zhuǎn)捩模型中, 函數(shù)或函數(shù)中的參數(shù)總是通過一系列零壓力梯度平板算例標(biāo)定得到. 在本工作中, 零壓力梯度自然轉(zhuǎn)捩平板(S&K和T3A-)扮演智能算法訓(xùn)練數(shù)據(jù)來源的角色.

2.3 湍流模型

一方程SA模型控制方程為

其中ν為分子運動黏性系數(shù),表示與湍流渦黏系數(shù)νt相關(guān)的修正渦黏系數(shù), 具體關(guān)系由下式定義:

生成項Pν與破壞項Dν定義為

其中,d為最小壁面距離.為渦量W的函數(shù):

上述方程中, 常數(shù)取值如下:

參考相關(guān)性間歇轉(zhuǎn)捩模型的思想, 本文將間歇因子作用于湍流模型以抑制轉(zhuǎn)捩前和轉(zhuǎn)捩過程中的湍流生成, 使得原全湍流結(jié)果向轉(zhuǎn)捩流場修正.考慮間歇因子描述流場某點處間歇狀態(tài)的客觀性定義, 通過神經(jīng)網(wǎng)絡(luò)描述其潛在函數(shù)關(guān)系并與SA湍流模型結(jié)合具備合理性. 對于SA湍流模型的輸運方程, 間歇因子用于抑制渦黏的生成項與破壞項:

其中,Pn和Dn分別是原SA模型的生成項和破壞項,b為破壞項修正系數(shù).

SA 模型不計算當(dāng)?shù)赝膭幽? 進(jìn)而無法考慮湍流度的影響. Medida等[31]將g-Req模型耦合至SA 模型時, 假定湍流度在整個流場中不改變. 本工作同樣采取這種方式, 并將自由流湍流度作為ResNet的輸入特征之一.

圖 2 轉(zhuǎn)捩平板計算網(wǎng)格Fig. 2. Computing mesh for transition plate.

圖 3 模型參數(shù)對壁面摩阻的影響 (a) nt∞; (b) bFig. 3. The influence of model parameters on wall friction:(a) nt∞; (b) b.

在上述修改的基礎(chǔ)上, 為保證間歇因子與SA模型的兼容性, 需要對破壞項系數(shù)b和SA模型流場變量遠(yuǎn)場值nt∞做相應(yīng)調(diào)整. 對于新的幾何外形和來流條件, 訓(xùn)練完成的ResNet模型能夠快速給出所需的g值, 而不需要再借助其他轉(zhuǎn)捩模型.

用于數(shù)值模擬的轉(zhuǎn)捩平板網(wǎng)格見圖2, 大小為324 × 108 (流向 × 法向), 平板上分布 291 網(wǎng)格單元. 如圖 3(a)所示, 在自然轉(zhuǎn)捩情況下, 過大的nt∞會導(dǎo)致摩阻曲線提早過渡到湍流狀態(tài), 因此遠(yuǎn)場值最終設(shè)置為 1.0 × 10–8. 如圖 3(b)所示, 原全湍流結(jié)果(b= 1)在湍流區(qū)域較實驗值整體偏低,破壞項修正系數(shù)b的作用就在于控制湍流區(qū)域的摩阻系數(shù)漲幅, 本文取b= 0.3.

3 間歇因子模型

上述過程通過引入間歇因子賦予SA模型預(yù)測轉(zhuǎn)捩的能力. 其中的關(guān)鍵在于, 對于一個全新的幾何外形和自由來流條件, 能否快速通過某種算法或模型給出足夠準(zhǔn)確的間歇因子預(yù)測場. 本節(jié)利用自由來流條件和局部無量綱特征量等信息, 試圖建立一種具備理想泛化能力的間歇因子模型.

3.1 模型選擇

在轉(zhuǎn)捩問題中, 流場某一空間位置的間歇狀態(tài)與當(dāng)?shù)氐钠骄兞块g存在復(fù)雜的映射關(guān)系, 且這種關(guān)系是非線性的, 難以建立解析的表達(dá)式. 為了滿足工程需求, 往往需要建立數(shù)學(xué)模型以近似描述這類映射.

傳統(tǒng)間歇轉(zhuǎn)捩模型的構(gòu)建往往基于將轉(zhuǎn)捩的起始與某個或某些特征量相關(guān)聯(lián), 最終也都證明模型具有良好的效果. 這些設(shè)計的關(guān)鍵特征與流場轉(zhuǎn)捩行為的關(guān)聯(lián)性直接影響到模型的準(zhǔn)確程度和適用范圍. Bas等[16]就在此類工作上付出過努力, 他們構(gòu)建的局部量Term1比較了當(dāng)?shù)貏恿亢穸萊eq是否達(dá)到轉(zhuǎn)捩的臨界動量厚度Reqc, 以判斷轉(zhuǎn)捩起始位置; Term2控制邊界層內(nèi)的間歇因子生成;最后, 通過指數(shù)函數(shù)生成間歇因子. 這種簡單的代數(shù)模型在部分情況下能夠達(dá)到接近g-Req模型的效果, 但對于湍流度在(0.5, 2.0)區(qū)間內(nèi)的轉(zhuǎn)捩問題并不適用. 不過該研究足以表明, 構(gòu)建合適的局部特征對轉(zhuǎn)捩預(yù)測起到關(guān)鍵作用.

由于轉(zhuǎn)捩問題的復(fù)雜性和混沌性, 人工構(gòu)建普適的相關(guān)特征量以構(gòu)建模型絕非易事. 幸運的是,當(dāng)前高性能計算的發(fā)展使得神經(jīng)網(wǎng)絡(luò)模型的訓(xùn)練更加高效, 能夠借助其強(qiáng)大的特征學(xué)習(xí)和表示能力構(gòu)建模型. 具體來說, 神經(jīng)網(wǎng)絡(luò)的每一層都對上一層傳遞來的信號進(jìn)行加工, 這些本與轉(zhuǎn)捩不直接相關(guān)的初始特征經(jīng)過多隱藏層的處理, 逐步轉(zhuǎn)換為與轉(zhuǎn)捩密切相關(guān)的特征量. 輸入層和所有隱藏層可看成“特征學(xué)習(xí)”的過程, 而最后一層的輸出層僅僅是以生成的特征量為輸入的簡單模型[32], 最終完成間歇因子的預(yù)測. 這種多層激活函數(shù)嵌套的結(jié)構(gòu)普遍存在于深層的網(wǎng)絡(luò)結(jié)構(gòu)之中. 這便是深度神經(jīng)網(wǎng)絡(luò)為什么適合于本研究的重要原因.

本工作希望能夠利用已有高精度計算數(shù)據(jù)中的信息, 通過深層次的網(wǎng)絡(luò)結(jié)構(gòu)和輔助優(yōu)化算法提取到與轉(zhuǎn)捩密切相關(guān)的流場特征, 從而構(gòu)建出當(dāng)?shù)靥卣髁颗c流場間歇性之間的類代數(shù)映射模型, 以獲得高精度和高效性并存的效果.

3.2 殘差神經(jīng)網(wǎng)絡(luò)

深度學(xué)習(xí)作為機(jī)器學(xué)習(xí)技術(shù)的典型代表, 能夠挖掘數(shù)據(jù)集中潛藏的深層次非線性映射并建立模型, 實現(xiàn)判斷和預(yù)測. 殘差神經(jīng)網(wǎng)絡(luò)在一般神經(jīng)網(wǎng)絡(luò)的基礎(chǔ)上引入跨層連接(skip connection), 使得誤差梯度得以跨層傳遞, 是一種典型的深度學(xué)習(xí)算法. 為保證理論的完整性, 引入跨層連接的優(yōu)勢將在附錄中展開說明. 基本網(wǎng)絡(luò)結(jié)構(gòu)可由一個輸入層、若干隱藏層、一個輸出層和跨層連接組成, 如圖4所示.

輸入層由一組代表流場不同屬性的當(dāng)?shù)仄骄縬= (q1,q2, …,qn)組成. 這些輸入量經(jīng)歷帶權(quán)重的連接到達(dá)下一隱藏層, 新的輸入值在層中節(jié)點將與閾值比較, 然后通過非線性激活函數(shù)轉(zhuǎn)換. 特征信息每經(jīng)過一個隱藏層都會被變換到新的特征空間, 直到輸出結(jié)果. 以圖4所示的神經(jīng)網(wǎng)絡(luò)模型為例, 對于一個新的輸入xi, 第一隱藏層的輸出分量可表示為

其中, 上標(biāo)(m)表示第m個隱藏層,?為非線性激活函數(shù), 本文采用 ReLU 函數(shù);wij為上層與本層間的連接權(quán);cj為神經(jīng)元節(jié)點處的閾值. 同樣, 第二、三隱藏層的輸出分量可分別寫為

其中m1,m2分別表示第一、二隱藏層的神經(jīng)元節(jié)點數(shù). 可以看到, 第一隱藏層的輸出通過跨層連接與第三隱藏層未經(jīng)激活的輸出相加, 合并的數(shù)據(jù)流經(jīng)過激活函數(shù)后向下傳遞, 這種跨層連接實際上就是在層間增加一個恒等映射, 最終第四隱藏層的輸出分量可表示為

圖 4 殘差神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)示意圖Fig. 4. Structure of residual neural network.

得益于上述多層次、非線性的激活變換, 網(wǎng)絡(luò)得以具備描述輸入特征與輸出量間的深層次非線性映射的能力. 常用的非線性變換函數(shù)(或稱激活函數(shù))有Sigmoid, ReLU和TANH等, 具體形式如下.

Sigmoid函數(shù)定義為

ReLU函數(shù)定義為

Tanh函數(shù)定義為

激活函數(shù)能夠在網(wǎng)絡(luò)中引入非線性變換, 且連續(xù)可導(dǎo), 便于后續(xù)采用梯度下降等優(yōu)化算法訓(xùn)練神經(jīng)網(wǎng)絡(luò), 曲線如圖 5所示. 對于回歸預(yù)測, 輸出層往往不需添加非線性函數(shù), 直接按照一定權(quán)值將特征矢量線性組合并輸出結(jié)果.

損失函數(shù)是直接反映模型精度的重要指標(biāo), 本文采用損失函數(shù)設(shè)置為

對于給定的訓(xùn)練數(shù)據(jù)集, 可以通過梯度下降法對網(wǎng)絡(luò)節(jié)點的權(quán)重wij與閾值cj進(jìn)行優(yōu)化:

其中f由不同優(yōu)化算法決定, 這里采用Adam算法, 具體可參考文獻(xiàn)[33]. 一旦確定了網(wǎng)絡(luò)節(jié)點的最優(yōu)參數(shù), 模型的預(yù)測效果則完全取決于隱藏層、每層節(jié)點數(shù)、激活函數(shù)、權(quán)重與閾值等超參數(shù), 而與訓(xùn)練數(shù)據(jù)的規(guī)模無關(guān).

圖 5 激活函數(shù)曲線 (a) Sigmoid; (b) ReLU; (c) TanhFig. 5. Activation function curves: (a) Sigmoid; (b) ReLU;(c) Tanh.

本研究的神經(jīng)網(wǎng)絡(luò)部分基于TensoFlow框架.輸入層由自由來流條件、當(dāng)?shù)亓鲌銎骄考捌鋵?dǎo)數(shù)等特征構(gòu)成, 輸入特征量的詳細(xì)描述列于表1. 此外, 經(jīng)過反復(fù)的實驗權(quán)衡計算成本與損失函數(shù)之間的平衡, 采用的網(wǎng)絡(luò)結(jié)構(gòu)由淺到深依次由6個全連接層、兩個內(nèi)含8個隱藏層的殘差塊和8個全連接層組成, 每個隱藏層有24個神經(jīng)元節(jié)點, 激活函數(shù)為ReLU. 輸出層為無激活的單一變量g.

神經(jīng)元和隱藏層數(shù)的增加伴隨著過擬合的風(fēng)險, 研究通過5次交叉驗證確定超參數(shù), 確保模型沒有出現(xiàn)過擬合. 訓(xùn)練數(shù)據(jù)被隨機(jī)地均分為5組,其中一組作為驗證集, 剩余的分組作為訓(xùn)練集,5次驗證的平均絕對誤差列于表2. 輸入特征采用當(dāng)?shù)責(zé)o量綱量以保證泛化性. 為了降低各特征間的相關(guān)性, 降低輸入冗余, 在網(wǎng)絡(luò)訓(xùn)練前將輸入數(shù)據(jù)進(jìn)行了白化(whiting)處理, 使得其均值為0, 方差為1.

表 1 作為神經(jīng)網(wǎng)絡(luò)輸入的流場局部平均特征量Table 1. The local average flow features used as the inputs of neural network.

表 2 5 次交叉驗證結(jié)果Table 2. Results of fivefold cross validation.

4 算例分析

計算基于自主研發(fā)的WCNS軟件平臺. WCNS格式于2000年由Deng等[28]提出, 其優(yōu)勢在廣泛的流動問題中得到證明. 本文采用的WCNS-E6E5格式具備良好的渦保持特性和很強(qiáng)的魯棒性, 在各類工程流動問題中得到廣泛應(yīng)用. 其中, 王圣業(yè)等[35]將WCNS-E6E5格式應(yīng)用于三角翼的分離渦模擬中, 并和傳統(tǒng)基于線性渦黏模型的分離渦模擬方法進(jìn)行了對比, 證明了格式在分離湍流模擬的精細(xì)度方面更加優(yōu)秀.

4.1 T3系列平板

T3系列低速平板實驗將用于測試上述數(shù)據(jù)驅(qū)動間歇因子模型對流場間歇性的預(yù)測能力并結(jié)合SA模型驗證最終的預(yù)測效果. 實驗包括S&K,T3A, T3A-和 T3B, 通常用于標(biāo)定轉(zhuǎn)捩模型中的經(jīng)驗系數(shù)和經(jīng)驗關(guān)系式. 本文主要關(guān)注于低湍流度自然轉(zhuǎn)捩實驗S&K 和T3A-的模擬. 采用的計算網(wǎng)格見圖 1, 網(wǎng)格為 324 × 108(流向 × 法向), 平板上分布 291 網(wǎng)格單元. 為適應(yīng)低速流動條件, 計算采用預(yù)處理技術(shù)[36]. 設(shè)置的進(jìn)口條件如表3.

表 3 平板算例入口條件Table 3. The entry condition of plate cases.

湍流邊界層具有更大的摩阻, 由此能夠通過摩阻曲線判斷邊界層的轉(zhuǎn)捩情況. 圖6給出了壁面摩阻的對比結(jié)果, 可以發(fā)現(xiàn)SA-ResNet模型和SST-g-Req模型能夠較好地預(yù)測自然轉(zhuǎn)捩, 而SA模型未能預(yù)測出轉(zhuǎn)捩, 且全湍流計算的壁面摩阻與實驗值相差很大. BC模型能夠預(yù)測到轉(zhuǎn)捩現(xiàn)象, 但過早地預(yù)測了T3A-平板算例的轉(zhuǎn)捩位置. 原因在于模型很大程度上依賴人為標(biāo)定的經(jīng)驗公式, 這在一定范圍內(nèi)增加了模型的不確定度, 影響了適用性.同為一方程的SA-ResNet模型在保證了求解效率的同時, 避免了上述問題.

值得注意的是, ResNet模型完全基于SST-g-Req的計算數(shù)據(jù)進(jìn)行學(xué)習(xí), 但預(yù)測結(jié)果并不完全相同. 這是由于兩者基于的全湍流模型存在原理上的差異, 再者SA-ResNet模型中的破壞項系數(shù)和流場變量初值對轉(zhuǎn)捩位置、轉(zhuǎn)捩區(qū)長度和湍流區(qū)的摩阻幅值存在些許影響, 這在第2節(jié)中進(jìn)行過說明.

圖 6 轉(zhuǎn)捩平板壁面摩阻曲線 (a) S&K; (b) T3AFig. 6. Wall friction curve of transition plate cases: (a) S&K; (b) T3A-.

圖 7 T3A-平板間歇因子和湍流黏性分布 (a) SST-g-Req 預(yù)測g 場; (b) SA-ResNet預(yù)測g 場; (c) SST-g-Req 和 SA-ResNet預(yù)測g 的差異; (d) SA-ResNet預(yù)測的湍流黏性Fig. 7. Intermittency and turbulent viscosity distribution of T3A- case: (a) g from SST-g-Req; (b) g from SA-ResNet; (c) discrepancy of g between SST-g-Req and SA-ResNet; (d) turbulent viscosity from SA-ResNet.

圖7(a),(b)比較了SA-ResNet模型和SST-g-Req模型對于T3A-轉(zhuǎn)捩平板間歇因子預(yù)測結(jié)果.兩模型結(jié)果的差異見圖7(c). 可以看出, ResNet模型能夠模擬出平板表面的層流邊界層發(fā)展和轉(zhuǎn)捩的過程, 模擬得到的邊界層厚度和轉(zhuǎn)捩位置與SST-g-Req模型基本一致. 在全湍流區(qū)上游壁面附近ResNet模型的預(yù)測結(jié)果出現(xiàn)不光滑的情況.原因在于壁面附近網(wǎng)格較為密集, 同時層流與湍流的分界區(qū)域間歇因子急劇變化, 這種情況下, 神經(jīng)網(wǎng)絡(luò)這類模糊預(yù)測方法有時難以清晰地給出分界線. 但這并沒有影響結(jié)果的正確性, 間歇因子正確地控制了流場中湍流的生成, 如圖7(d). 這正是不直接以渦黏性或雷諾應(yīng)力作為神經(jīng)網(wǎng)絡(luò)預(yù)測對象的優(yōu)勢所在.

4.2 S809翼型

S809翼型的亞聲速定常繞流模擬將作為驗證模型泛化性能的典型算例. 該翼型為厚度21%c的層流翼型, 專門為橫軸風(fēng)力渦輪機(jī)設(shè)計. 翼型的低速試驗在戴爾福特科技大學(xué)(Delft University of Technology)的低湍流度風(fēng)洞中進(jìn)行[37]. 計算網(wǎng)格如圖 8 所示, 采用 C 型拓?fù)浣Y(jié)構(gòu), 共劃分約 6.6 萬網(wǎng)格單元, 壁面首層網(wǎng)格距離達(dá)到 1 × 10–6c, 遠(yuǎn)場邊界取 120c, 進(jìn)口馬赫數(shù)為 0.1, 雷諾數(shù)為 2.0 × 106,翼型前緣湍流度設(shè)為0.2%.

圖9對比了S809翼型氣動特性計算值和實驗值的結(jié)果. 由圖中可以看出, 結(jié)合了ResNet的SA模型與SST-g-Req模型在升阻力特性上都與實驗值更加相符. 不考慮轉(zhuǎn)捩的原SA模型預(yù)測的升力系數(shù)偏小, 且總是過多地預(yù)測了阻力. 結(jié)合ResNet的SA模型則很大程度上修正了這一情況.

圖 8 S809 翼型計算網(wǎng)格Fig. 8. Computing mesh for S809 airfoil.

圖 9 S809 翼型氣動特性曲線 (a) Cl; (b) CdFig. 9. Aerodynamic characteristics of the S809 airfoil: (a) Cl;(b) Cd.

由于訓(xùn)練集中不包含逆壓梯度和大分離的例子, 這里僅討論S809翼型迎風(fēng)面的轉(zhuǎn)捩位置情況.轉(zhuǎn)捩位置隨迎角變化曲線見圖10, 計算中的自然轉(zhuǎn)捩位置取摩阻最小值點, 分離轉(zhuǎn)捩的位置取摩阻由負(fù)值變正值時摩阻為零點. 由圖中可以看出, 在0°—3°范圍內(nèi)SA-ResNet模型的轉(zhuǎn)捩位置較實驗略微靠前, 其他迎角下計算值與實驗值符合較好.

圖 10 S809 翼型迎風(fēng)面轉(zhuǎn)捩位置隨迎角變化曲線Fig. 10. S809 airfoil transition position changes with the angle of attack.

圖11給出了不同迎角下翼型表面的摩阻系數(shù)分布, 由圖可以看出 SA 模型在 1°, –6°和 9°迎角下都未能預(yù)測出轉(zhuǎn)捩. 而SA-ResNet模型與SST-g-Req轉(zhuǎn)捩模型的結(jié)果十分貼近. 在 1°迎角下, 翼型上下翼面分別在0.550和0.526附近發(fā)生分離轉(zhuǎn)捩. 上翼面轉(zhuǎn)捩位置隨迎角增大不斷前移, 在9°迎角時轉(zhuǎn)捩位置到達(dá)0.01附近, 壁面摩阻沿流動方向不斷下降, 流動再度層流化. 下翼面轉(zhuǎn)捩位置隨迎角不斷向后緣移動.

上述結(jié)果驗證了SA-ResNet模型的轉(zhuǎn)捩預(yù)測能力. 在此基礎(chǔ)上, 模型的另一優(yōu)勢在于求解效率.相較于四方程SST-g-Req轉(zhuǎn)捩模型, SA-ResNet模型通過代數(shù)“黑箱”模型給出間歇因子分布, 不需引入額外的微分方程, 降低了計算耗時. 效率提升在大網(wǎng)格量的算例上更為明顯. 表4列出了S&K、T3A-轉(zhuǎn)捩平板和S809翼型3°迎角下的計算時間.對于S&K平板算例, SA-ResNet模型的計算耗時為SST-g-Req轉(zhuǎn)捩模型的85.6%, 然而在網(wǎng)格量更大的S809翼型上, 前者所需CPU時間僅為后者的67.2%. Wang等[38]結(jié)合了BC代數(shù)轉(zhuǎn)捩模型和WALE子網(wǎng)格模型對三維圓柱繞流問題進(jìn)行了模擬, 結(jié)果發(fā)現(xiàn)每次迭代所需CPU時間不到SST-g-Req轉(zhuǎn)捩模型的30%. 這表明在流動拓展到三維后, 代數(shù)轉(zhuǎn)捩模型帶來的效率優(yōu)勢十分可觀.

圖 11 S809 翼型不同迎角下的摩阻系數(shù)曲線Fig. 11. Friction coefficient of S809 airfoil at different angles of attack.

表 4 模型計算時間對比 (殘差收斂至 O(10–4))Table 4. Comparison of transition model’s computing time.

5 結(jié)論和展望

本文基于 WCNS-E6 E5格式, 結(jié)合深度殘差神經(jīng)網(wǎng)絡(luò)方法和SA湍流模型, 參考基于當(dāng)?shù)亓鲌鲎兞康南嚓P(guān)性間歇模型的構(gòu)造思路, 重構(gòu)了間歇因子映射模型, 發(fā)展了一種基于數(shù)據(jù)驅(qū)動的SAResNet高精度轉(zhuǎn)捩模型. 對S&K、T3A-轉(zhuǎn)捩平板試驗和S809翼型不同迎角下的亞聲速定常繞流進(jìn)行了模擬, 并與SST-g-Req轉(zhuǎn)捩模型計算值和實驗值進(jìn)行了對比.

轉(zhuǎn)捩平板算例表明, 深度殘差網(wǎng)絡(luò)方法能夠通過計算數(shù)據(jù)描述流場平均量與間歇因子間的復(fù)雜非線性映射; 得到的間歇因子能夠較好地與SA模型兼容, 準(zhǔn)確捕捉流場中的轉(zhuǎn)捩現(xiàn)象. 在有限訓(xùn)練數(shù)據(jù)的情況下, SA-ResNet對于未包含于訓(xùn)練集的S809翼型的不同迎角情況能夠得出較為準(zhǔn)確的預(yù)測. 這表明深度殘差網(wǎng)絡(luò)能夠做到的并不僅僅是對數(shù)據(jù)進(jìn)行簡單的插值, 其具備探索流場變量與轉(zhuǎn)捩過程間復(fù)雜關(guān)系的強(qiáng)大潛力.

在S809翼型亞聲速定常繞流的模擬中, SAResNet的性能接近SST-g-Req, 但收斂到同一精度時節(jié)省了超過30%的計算成本. 下一步工作在于向多種轉(zhuǎn)捩機(jī)理的流動和不同來源的計算數(shù)據(jù)拓展, 強(qiáng)化模型的工程應(yīng)用能力.

附錄

一般而言, 更多的神經(jīng)元單元和隱藏層能夠提升網(wǎng)絡(luò)性能同時加速訓(xùn)練收斂. 然而, 網(wǎng)絡(luò)增加到一定層數(shù)后, 在利用反向傳播算法訓(xùn)練網(wǎng)絡(luò)的過程中容易出現(xiàn)梯度消失/爆炸問題[39,40], 導(dǎo)致訓(xùn)練無法收斂. 例如對于一般的全連接神經(jīng)網(wǎng)絡(luò), 有

其中w為連接權(quán),c為偏差,z為網(wǎng)絡(luò)層的輸出,s為激活函數(shù)輸入值,?為激活函數(shù). 因此, 在誤差梯度反向傳播過程中, 對權(quán)重有

對閾值有

對于多層網(wǎng)絡(luò):

誤差梯度中存在激活函數(shù)導(dǎo)數(shù)和權(quán)重的連乘, 當(dāng)網(wǎng)絡(luò)層數(shù)增長到一定規(guī)模以后, 淺層隱藏層接收到的梯度由于一系列連乘而劇烈衰減/增長, 從而引發(fā)梯度消失或梯度爆炸問題. 這種問題的出現(xiàn)阻礙了訓(xùn)練的收斂, 甚至導(dǎo)致網(wǎng)絡(luò)退化. 這種退化是由過多的隱藏層造成而非過擬合.

由于跨層連接的存在, 被連接包圍的網(wǎng)絡(luò)層(殘差塊)僅學(xué)習(xí)連接前后的殘差映射(因此稱為殘差網(wǎng)絡(luò)), 誤差梯度能夠越過殘差塊反向傳播, 減少傳播路徑, 從而降低梯度消失或梯度爆炸的風(fēng)險.