基于深度神經(jīng)網(wǎng)絡(luò)的橫流轉(zhuǎn)捩預(yù)測1)

胡震宇 王子路 陳堅(jiān)強(qiáng) 袁先旭 向星皓

(空氣動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室,四川綿陽 621000)

引言

轉(zhuǎn)捩是一個(gè)連續(xù)演化的流動過程,是一種從簡單分層穩(wěn)定狀態(tài)向復(fù)雜多變湍流狀態(tài)的過渡[1].轉(zhuǎn)捩研究不僅具有重要的理論意義,還具有重要的工程應(yīng)用價(jià)值.轉(zhuǎn)捩過程受多種模態(tài)影響[2],由于其復(fù)雜性,將長期作為現(xiàn)代流體力學(xué)重點(diǎn)關(guān)注的前沿領(lǐng)域而存在[3].

關(guān)于流體力學(xué)的研究,飛行試驗(yàn)、地面實(shí)驗(yàn)具有直觀、真實(shí)的特點(diǎn),能夠?yàn)榉€(wěn)定性理論、數(shù)值計(jì)算提供較為準(zhǔn)確的驗(yàn)證數(shù)據(jù)[4-6].直接數(shù)值模擬(direct numerical simulation,DNS)方法[7-10]直接求解N-S 方程,能夠準(zhǔn)確捕捉流場結(jié)構(gòu),但對計(jì)算資源要求較高;大渦模擬(large eddy simulation,LES)方法[11]對大尺度渦進(jìn)行解析、對小尺度渦進(jìn)行模化,計(jì)算結(jié)果較為準(zhǔn)確,但依然要求較高的網(wǎng)格分辨率;雷諾平均N-S (Reynolds averaged Navier-Stokes,RANS)方法[12-14]中,利用湍流模型封閉雷諾應(yīng)力項(xiàng),降低了對網(wǎng)格分辨率的要求,應(yīng)用較為廣泛.

基于RANS 方法的k-ε,k-ω,SSTk-ω等傳統(tǒng)湍流模型,能夠一定程度上滿足工程需要,但無法捕捉邊界層內(nèi)轉(zhuǎn)捩的過程.對于轉(zhuǎn)捩過程的描述,通常用到間歇因子γ[15]

式中,I(x,y,z,t)表征某一空間點(diǎn)在某一時(shí)刻的流動狀態(tài),處于湍流狀態(tài)時(shí)I(x,y,z,t)=1,處于層流狀態(tài)時(shí)I(x,y,z,t)=0.

間歇因子γ表征某一空間點(diǎn)處于湍流狀態(tài)的時(shí)間占總時(shí)間T的比例.通過求解間歇因子γ的輸運(yùn)方程,可以獲得間歇因子γ在流場中的分布情況.轉(zhuǎn)捩模型通過對間歇因子γ或其他參數(shù)進(jìn)行建模,并與湍流模型耦合求解,實(shí)現(xiàn)對轉(zhuǎn)捩過程的預(yù)測.

Steelant 等[15]針對bypass 轉(zhuǎn)捩,構(gòu)建了間歇因子γ的輸運(yùn)方程,計(jì)算結(jié)果能夠與實(shí)驗(yàn)數(shù)據(jù)較好地吻合.目前應(yīng)用較廣的轉(zhuǎn)捩模型,也是通過構(gòu)建間歇因子γ的輸運(yùn)方程來實(shí)現(xiàn)轉(zhuǎn)捩預(yù)測的.例如Langtry等[16]建立的γ-Reθ轉(zhuǎn)捩模型,王亮等[17]建立的k-ω-γ轉(zhuǎn)捩模型.

神經(jīng)網(wǎng)絡(luò)通過模仿生物的大腦神經(jīng)元結(jié)構(gòu)和功能,使計(jì)算機(jī)對數(shù)據(jù)信息能夠進(jìn)行智能處理,構(gòu)建復(fù)雜、非線性映射關(guān)系[18].相比于求解輸運(yùn)方程,深度神經(jīng)網(wǎng)絡(luò)(deep neural networks,DNN)所需的計(jì)算資源更少,顯著提高計(jì)算效率,在湍流模型化工作中具有良好的前景[19-20].

Templeton 等[21]提出了張量基神經(jīng)網(wǎng)絡(luò) (tensor based neural network,TBNN) 架構(gòu),并用于預(yù)測低速的管流、周期山流動.任海杰等[22]基于文獻(xiàn)[21]搭建的 TBNN 內(nèi)核構(gòu)建了神經(jīng)網(wǎng)絡(luò)模型,并針對Boussinesq 假設(shè)難以準(zhǔn)確捕捉雷諾應(yīng)力各向異性張量的問題,利用Pope[23]給出的各向異性張量b與基張量T(n)之間的本構(gòu)關(guān)系,構(gòu)建了張量不變量λ1～λ5與標(biāo)量函數(shù)g(n)間的映射關(guān)系,對高超聲速平板進(jìn)行預(yù)測,得到的各向異性張量分布與DNS 計(jì)算結(jié)果吻合較好.Wei 等[24]提出了基于深度強(qiáng)化學(xué)習(xí)的流體力學(xué)微分方程統(tǒng)一求解框架,求解了流體力學(xué)Burgers 方程、穩(wěn)態(tài)N-S 方程.Sekar 等[25]基于深度學(xué)習(xí)方法,針對具有不同外形的機(jī)翼,構(gòu)建了繞流流場預(yù)測模型.Zhu 等[26]通過徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)以及深度神經(jīng)網(wǎng)絡(luò),對渦黏系數(shù)進(jìn)行了預(yù)測,并用于封閉湍流模型.張珍等[27]提出了一種基于組合神經(jīng)網(wǎng)絡(luò)預(yù)測渦黏系數(shù)和雷諾應(yīng)力各向異性張量的修正方法,對雷諾應(yīng)力的線性部分進(jìn)行了求解.Beetham等[28]運(yùn)用傳統(tǒng) RANS 方法求解RANS 模型線性部分,運(yùn)用神經(jīng)網(wǎng)絡(luò)預(yù)測RANS 模型非線性部分,但缺乏足夠精度.

目前人工智能在計(jì)算流體力學(xué)領(lǐng)域的研究,主要集中在重構(gòu)或修正湍流渦黏性和雷諾應(yīng)力[29].隨著學(xué)科間的交叉與發(fā)展,將神經(jīng)網(wǎng)絡(luò)應(yīng)用于邊界層轉(zhuǎn)捩問題求解,成為一個(gè)新的熱門研究方向.鄭天韻等[29]基于深度殘差網(wǎng)絡(luò)(deep residual network,ResNet),構(gòu)建了間歇因子γ與流場平均量間的映射關(guān)系,并與spallart allmaras (SA)湍流模型耦合,發(fā)展了一種類代數(shù)轉(zhuǎn)捩模型,性能接近SST-γ-Reθ模型,但收斂到同一精度節(jié)省了超過30% 的計(jì)算成本.Foroozan 等[30]運(yùn)用Johns Hopkins 湍流數(shù)據(jù)庫(Johns Hopkins turbulence databases,JHTDB)構(gòu)建了訓(xùn)練數(shù)據(jù)集,采用無監(jiān)督機(jī)器學(xué)習(xí)的方法,實(shí)現(xiàn)了對邊界層bypass 轉(zhuǎn)捩的識別.

本文DNN 基于Pytorch 平臺.數(shù)值計(jì)算基于Chant 2.0 平臺[31-32],該平臺部署了橫流拓展γ-Reθ轉(zhuǎn)捩模型[33].橫流拓展γ-Reθ轉(zhuǎn)捩模型在原始γ-Reθ轉(zhuǎn)捩模型[16]上進(jìn)行改進(jìn),在NLF(2)-0415 后掠翼、6 : 1 標(biāo)準(zhǔn)橢球體、鐮刀翼和DLR-F4 翼身組合體等典型橫流轉(zhuǎn)捩算例上均有著良好的計(jì)算效果,能做到與實(shí)驗(yàn)數(shù)據(jù)極大程度吻合[34-35].后文直接使用γ-Reθ轉(zhuǎn)捩模型,代指部署于Chant 2.0 平臺的橫流拓展γ-Reθ轉(zhuǎn)捩模型.

本文運(yùn)用γ-Reθ轉(zhuǎn)捩模型對橫流轉(zhuǎn)捩模態(tài)下間歇因子γ的數(shù)值與分布進(jìn)行計(jì)算,并對DNN 進(jìn)行監(jiān)督學(xué)習(xí),構(gòu)建從層流流場物理量到轉(zhuǎn)捩流場間歇因子γ間的映射關(guān)系,從而得到一種新的數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型.僅需將數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型與SSTk-ω湍流模型耦合,即可對橫流不穩(wěn)定性主導(dǎo)的邊界層轉(zhuǎn)捩進(jìn)行預(yù)測,簡化計(jì)算過程.

本文采用數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型與湍流模型耦合求解的手段,研究了人工智能轉(zhuǎn)捩計(jì)算方法對橫流轉(zhuǎn)捩的預(yù)測能力,以期為今后的人工智能算法與CFD耦合工作提供方法參考.

1 模型構(gòu)建

1.1 DNN 監(jiān)督學(xué)習(xí)框架

本文參照γ-Reθ轉(zhuǎn)捩模型求解間歇因子γ,修正SSTk-ω模型湍動能k輸運(yùn)方程的生成項(xiàng)與耗散項(xiàng),進(jìn)而實(shí)現(xiàn)轉(zhuǎn)捩預(yù)測的思路,利用DNN 構(gòu)建來流條件、壁面信息和流場結(jié)構(gòu)等物理信息與間歇因子γ間的映射關(guān)系,用新的數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型替代γ-Reθ轉(zhuǎn)捩模型,實(shí)現(xiàn)人工智能建模框架下的三維邊界層橫流轉(zhuǎn)捩的高效準(zhǔn)確預(yù)測.

目前已有許多基于層流解的邊界層轉(zhuǎn)捩研究工作[36-39].以本文所用的γ-Reθ轉(zhuǎn)捩模型構(gòu)造為例,在發(fā)生轉(zhuǎn)捩前流場處于全層流狀態(tài),通過迭代層流流場的物理量從而開啟橫流源項(xiàng)DSCF實(shí)現(xiàn)轉(zhuǎn)捩預(yù)測,即層流流場中存在能夠用于判斷橫流轉(zhuǎn)捩起始位置的點(diǎn).因此本文選擇使用層流流場物理量映射轉(zhuǎn)捩流場間歇因子,訓(xùn)練數(shù)據(jù)的具體獲取方式如下.首先對層流狀態(tài)進(jìn)行計(jì)算,快速得到收斂穩(wěn)定的層流流場,并提取訓(xùn)練所需的物理量.再通過SSTk-ω湍流模型與γ-Reθ轉(zhuǎn)捩模型計(jì)算轉(zhuǎn)捩流場,提取各個(gè)流場點(diǎn)對應(yīng)的γ值.將γ與上一步提取的層流流場物理量進(jìn)行匹配,用于監(jiān)督學(xué)習(xí).DNN 通過梯度下降法,計(jì)算γ預(yù)測值與監(jiān)督學(xué)習(xí)的γ真實(shí)值(即γ-Reθ轉(zhuǎn)捩模型計(jì)算得到的γ值)之間的損失函數(shù),并對DNN各節(jié)點(diǎn)的超參數(shù)進(jìn)行修正,使損失函數(shù)降低至設(shè)定值以下,從而獲得能夠滿足轉(zhuǎn)捩計(jì)算精度要求的映射關(guān)系,該映射關(guān)系即為本研究的數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型.DNN 監(jiān)督學(xué)習(xí)功能實(shí)現(xiàn)的具體過程見圖1.

1.2 湍流與轉(zhuǎn)捩模型

本文以RANS 方程對變雷諾數(shù)、變后掠角的NFL(2)-0415 后掠翼以及6 : 1 標(biāo)準(zhǔn)橢球體進(jìn)行求解,采用兩方程的剪切應(yīng)力輸運(yùn)SSTk-ω湍流模型,并與γ-Reθ轉(zhuǎn)捩模型耦合.γ-Reθ轉(zhuǎn)捩模型由間歇因子γ輸運(yùn)方程、轉(zhuǎn)捩動量厚度雷諾數(shù)輸運(yùn)方程構(gòu)成,具體構(gòu)造如下

橫流轉(zhuǎn)捩源項(xiàng)的常系數(shù)CCF參照文獻(xiàn)[35]的研究結(jié)果取0.2,其余常系數(shù)、參數(shù)與文獻(xiàn)[34]保持一致.

γ-Reθ模型與SSTk-ω湍流模型的耦合,通過有效間歇因子γeff對湍動能k方程生成項(xiàng)和耗散項(xiàng)修正實(shí)現(xiàn).SSTk-ω湍流模型的湍動能k方程被修改為

SSTk-ω湍流模型中的融合函數(shù)F1也需要修改

式中,Pk,Dk和Florig均為SSTk-ω湍流模型的原始形式[40].

1.3 模型計(jì)算能力驗(yàn)證

本文利用γ-Reθ轉(zhuǎn)捩模型對間歇因子γ進(jìn)行計(jì)算,并將計(jì)算結(jié)果作為γ真實(shí)值進(jìn)行監(jiān)督學(xué)習(xí),參與DNN 訓(xùn)練與測試.得到的數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型與SSTk-ω湍流模型耦合,用于預(yù)測變后掠角的NFL(2)-0415 后掠翼、6 : 1 標(biāo)準(zhǔn)橢球體的三維邊界層橫流轉(zhuǎn)捩.首先需要對γ-Reθ模型性能進(jìn)行校驗(yàn),以證明其計(jì)算結(jié)果的可靠性.

NFL(2)-0415 后掠翼在?4°攻角(attract of angle,AoA)、45°后掠角(sweep angle,Λ)下[41],上翼面的轉(zhuǎn)捩由橫流不穩(wěn)定性主導(dǎo),是低速橫流轉(zhuǎn)捩典型算例[42].選取該算例,運(yùn)用γ-Reθ轉(zhuǎn)捩模型對變雷諾數(shù)的NFL(2)-0415 后掠翼上翼面轉(zhuǎn)捩位置進(jìn)行預(yù)測,并將計(jì)算所得的γ分布與層流流場物理量進(jìn)行匹配,參與DNN 的訓(xùn)練、測試.具體計(jì)算工況見表1.

表1 NFL(2)-0415 后掠翼計(jì)算條件Table 1 Conditions of NFL(2)-0415 airfoil

如圖2 所示,C 型結(jié)構(gòu)化網(wǎng)格流向i共有400 網(wǎng)格點(diǎn),其中上、下翼面、上、下尾跡各分布100 網(wǎng)格點(diǎn),網(wǎng)格法向j共有151 網(wǎng)格點(diǎn),第1 層網(wǎng)格法向間距y+<1,網(wǎng)格展向k為二層網(wǎng)格,網(wǎng)格數(shù)量為120 800.計(jì)算采用周期邊界條件,以實(shí)現(xiàn)無線展長后掠翼的數(shù)值模擬.網(wǎng)格在邊界層內(nèi)及邊界層附近保證了良好的正交性,粗黑曲線與壁面間的網(wǎng)格點(diǎn)數(shù)據(jù)用于神經(jīng)網(wǎng)絡(luò)的訓(xùn)練.

圖2 NLF(2)-0415 后掠翼計(jì)算網(wǎng)格Fig.2 Computing mesh for NLF(2)-0415 airfoil

本文采用層流區(qū)域與轉(zhuǎn)捩區(qū)域Cf曲線延長線交點(diǎn)的無量綱化坐標(biāo)x/c作為轉(zhuǎn)捩位置xtr,如圖3所示.圖4 給出了本研究γ-Reθ轉(zhuǎn)捩模型的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果[41]、langtry-menter-CFH 模型[42]計(jì)算結(jié)果對比.

圖3 轉(zhuǎn)捩位置判定Fig.3 Determination of transition location

圖4 轉(zhuǎn)捩位置對比Fig.4 Comparison of transition location

由圖4 可知,γ-Reθ轉(zhuǎn)捩模型計(jì)算轉(zhuǎn)捩位置xtr在各工況下均小幅度先于Langtry-Menter-CFH 模型.γ-Reθ轉(zhuǎn)捩模型在低雷諾數(shù)下,計(jì)算轉(zhuǎn)捩位置xtr相較于實(shí)驗(yàn)結(jié)果略微提前,隨著雷諾數(shù)增大,計(jì)算轉(zhuǎn)捩位置xtr逐漸與實(shí)驗(yàn)結(jié)果貼近,最后幾乎重合.Langtry-Menter-CFH 模型在低雷諾數(shù)下計(jì)算轉(zhuǎn)捩位置xtr同樣存在略微提前的問題,隨著雷諾數(shù)增大,計(jì)算轉(zhuǎn)捩位置xtr出現(xiàn)滯后,雷諾數(shù)繼續(xù)增大,計(jì)算轉(zhuǎn)捩位置xtr逐漸與實(shí)驗(yàn)結(jié)果貼近并重合.

轉(zhuǎn)捩模型計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果間的差異,與模型本身構(gòu)建、模型中部分參數(shù)的標(biāo)定等諸多因素有關(guān),且不同研究中轉(zhuǎn)捩點(diǎn)坐標(biāo)xtr的判定方式同樣存在區(qū)別.盡管存在以上影響,γ-Reθ轉(zhuǎn)捩模型計(jì)算結(jié)果依然能與實(shí)驗(yàn)數(shù)據(jù)吻合較好,與Langtry-Menter-CFH 模型具有相仿的計(jì)算性能,能夠滿足預(yù)測精度和DNN 訓(xùn)練的需求.

2 深度神經(jīng)網(wǎng)絡(luò)

2.1 深度神經(jīng)網(wǎng)絡(luò)超參數(shù)

DNN 的輸入層為選定的層流流場物理量,并用γ-Reθ模型計(jì)算得到的間歇因子γ進(jìn)行監(jiān)督學(xué)習(xí).輸入層單個(gè)樣本的特征量為V[0]=(v1,v2,v3,···,vn),經(jīng)過隱藏層時(shí),在神經(jīng)元節(jié)點(diǎn)經(jīng)歷線性變換與激活函數(shù)的非線性變化,并將輸出值作為下一隱藏層的輸入,直至輸出最終結(jié)果,如圖5 所示.其中上角標(biāo)[j]表示第j個(gè)隱藏層,下角標(biāo)i表示該隱藏層的第i個(gè)神經(jīng)節(jié)點(diǎn).

圖5 深度神經(jīng)網(wǎng)絡(luò)正向傳播Fig.5 Positive propagation of DNN

訓(xùn)練樣本集中的單個(gè)樣本首先經(jīng)過如下線性變換與非線性變換

經(jīng)測試,在本研究中選取sigmoid 函數(shù)作為激活函數(shù)學(xué)習(xí)效果比ReLU 函數(shù)、tanh 函數(shù)更佳.非線性變換過程中的激活函數(shù)σ選取sigmoid 函數(shù),函數(shù)曲線如圖6 所示,定義如下

圖6 Sigmoid 函數(shù)Fig.6 Activation function of sigmoid

對于共有m個(gè)樣本的數(shù)據(jù)集,損失函數(shù)定義如下

通過梯度下降法,對神經(jīng)網(wǎng)絡(luò)的超參數(shù)進(jìn)行修正,使得損失函數(shù)L降低至設(shè)定值以下,從而實(shí)現(xiàn)神經(jīng)網(wǎng)絡(luò)的自主學(xué)習(xí)過程.

權(quán)重W修正如下

閾值b修正如下

其中函數(shù)f與Adam 優(yōu)化器有關(guān),在數(shù)據(jù)梯度較小的場景以較大的學(xué)習(xí)率進(jìn)行更新,在數(shù)據(jù)梯度較大的場景以較小的學(xué)習(xí)率進(jìn)行更新,比標(biāo)準(zhǔn)的SGD算法更有效地收斂.設(shè)置默認(rèn)學(xué)習(xí)率α=0.001,ε=1 × 10?10.

根據(jù)多次訓(xùn)練結(jié)果分析,當(dāng)隱藏層層數(shù)等于輸入層變量個(gè)數(shù),每層神經(jīng)元節(jié)點(diǎn)數(shù)等于輸入層變量個(gè)數(shù)3 倍時(shí),能獲得較好的訓(xùn)練效果.本文共選取7 個(gè)輸入變量,故DNN 隱藏層取7 層,每層均設(shè)置21 個(gè)神經(jīng)元節(jié)點(diǎn),輸出層為單一變量γpredict.

2.2 深度神經(jīng)網(wǎng)絡(luò)訓(xùn)練

本文DNN 訓(xùn)練基于Pytorch 平臺,構(gòu)建了層流流場物理量與間歇因子γ間的映射關(guān)系,目的在于實(shí)現(xiàn)數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型對三維流場橫流轉(zhuǎn)捩的精確預(yù)測.輸入層7 個(gè)變量V=R,Ωstreamwise),各變量定義如下.

用于尋找間歇因子在流向上增長位置的流向雷諾數(shù)

如圖7 所示,以NLF(2)-0415 后掠翼為例,特征距離s是翼型前緣點(diǎn)A到指定點(diǎn)B,壁面曲線的長度積分,點(diǎn)B處壁面法線L上所有點(diǎn)特征距離s與點(diǎn)B相同.將壁面曲線AB按照網(wǎng)格劃分進(jìn)行離散,點(diǎn)A作為第1 個(gè)離散點(diǎn),點(diǎn)B作為第n個(gè)離散點(diǎn),特征距離s如下

圖7 特征距離s 計(jì)算Fig.7 Calculation method of characteristic length s

用于尋找間歇因子在壁面法向上增長位置的無量綱壁面法向距離y+[34](通過程序Y+wall distance est imation 計(jì)算)

用于尋找間歇因子在壁面法向上增長位置的無量綱化的合速度

用于在流場結(jié)構(gòu)層面確定間歇因子增長位置、無量綱速度對應(yīng)的應(yīng)變率張量、旋轉(zhuǎn)率張量的模S,R[34]

在γ-Reθ轉(zhuǎn)捩模型中用于構(gòu)造無量綱橫流強(qiáng)度Hcrossflow、無量綱速度對應(yīng)的流向渦強(qiáng)度Ωstreamwise[34]

本文D N N 所用的訓(xùn)練數(shù)據(jù)為4 5°后掠角NLF(2)-0415 后掠翼層流狀態(tài)的計(jì)算數(shù)據(jù),以γ-Reθ轉(zhuǎn)捩模型計(jì)算得到的間歇因子γ進(jìn)行監(jiān)督學(xué)習(xí).前文已經(jīng)給出45°后掠角下5 種不同雷諾數(shù)工況,并對γ-Reθ轉(zhuǎn)捩模型計(jì)算的準(zhǔn)確性進(jìn)行了校驗(yàn).對于CFD 計(jì)算結(jié)果而言,影響占比較大的是邊界層內(nèi)以及邊界層附近的間歇因子γ分布情況,因此每種雷諾數(shù)下,取上翼面近壁面處i·j=101 × 81=8181 個(gè)網(wǎng)格點(diǎn)的數(shù)據(jù)用作訓(xùn)練以及測試,如圖8 所示,即圖中白色邊線范圍內(nèi)網(wǎng)格點(diǎn)的計(jì)算數(shù)據(jù),由速度分布可知,該區(qū)域包含了上翼面全部邊界層.

圖8 訓(xùn)練數(shù)據(jù)所處區(qū)域Fig.8 The zone of data for DNN

2.3 深度神經(jīng)網(wǎng)絡(luò)擬合校驗(yàn)

隨著隱藏層層數(shù)以及隱藏層節(jié)點(diǎn)數(shù)的增加,神經(jīng)網(wǎng)絡(luò)的學(xué)習(xí)能力增強(qiáng),容易過度挖掘訓(xùn)練樣本中特有的映射關(guān)系,并將其視為普適性規(guī)律,亦或是在已經(jīng)挖掘出恰當(dāng)映射關(guān)系后,進(jìn)一步學(xué)習(xí)并擬合噪音樣本,從而導(dǎo)致訓(xùn)練得到的DNN 泛化能力差,出現(xiàn)過擬合現(xiàn)象.現(xiàn)采用5 折交叉驗(yàn)證的方法,以確保本研究訓(xùn)練的DNN 沒有發(fā)生過擬合,5 折交叉驗(yàn)證方法如圖9 所示.

圖9 5 折交叉驗(yàn)證Fig.9 5-fold cross validation

5 種工況下共提取了40905 個(gè)數(shù)據(jù)樣本,隨機(jī)均勻劃分為5 組,每組包含8181 個(gè)數(shù)據(jù)樣本.依次將每組數(shù)據(jù)樣本作為測試集,余下的4 組數(shù)據(jù)樣本作為訓(xùn)練集,訓(xùn)練得到DNN1～ DNN5,并計(jì)算訓(xùn)練集、測試集的損失函數(shù)以及測試集損失函數(shù)的均值.5 折交叉驗(yàn)證的結(jié)果見表2.

表2 交叉驗(yàn)證結(jié)果Table 2 Result of cross validation

由于測試集中的數(shù)據(jù)樣本沒有參與DNN 訓(xùn)練,且DNN 必然存在對特有映射關(guān)系的構(gòu)建以及對訓(xùn)練集噪音樣本的學(xué)習(xí),交叉驗(yàn)證中測試集損失函數(shù)高于訓(xùn)練集是完全合理的.測試集損失函數(shù)、損失函數(shù)均值均維持在較低的水平,且與訓(xùn)練集損失函數(shù)間的差異遠(yuǎn)未到跨越一個(gè)或多個(gè)量級的程度,說明DNN 并沒有發(fā)生過擬合,能夠較為精準(zhǔn)地預(yù)測間歇因子γ的分布.

3 算例分析

3.1 變雷諾數(shù)NLF(2)-0415 后掠翼測試

將NLF(2)-0415 后掠翼Re=2.37 × 106以外的其他4 種工況共32724 個(gè)數(shù)據(jù)樣本作為訓(xùn)練集,對DNN 進(jìn)行訓(xùn)練.損失函數(shù)設(shè)置為7.50 × 10?3,完成30000步訓(xùn)練后,深度神經(jīng)網(wǎng)絡(luò)模型實(shí)際損失函數(shù)為2.0 ×10?4,將其作為本研究算例分析所用的數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型.以數(shù)據(jù)樣本是否參與DNN 訓(xùn)練作為內(nèi)插與外推的判別標(biāo)準(zhǔn),分別用Re=2.73 × 106,Re=2.37 ×106工況下的上翼面近壁區(qū)域數(shù)據(jù)樣本對數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型進(jìn)行內(nèi)插與外推測試.圖10 給出了兩種工況下數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型對間歇因子預(yù)測值與γ-Reθ轉(zhuǎn)捩模型計(jì)算值的對比.

圖10 間歇因子預(yù)測值與計(jì)算值對比Fig.10 Comparison between predicted value and calculated value of intermittency

間歇因子預(yù)測值與計(jì)算值間誤差越低,散點(diǎn)分布越接近線性函數(shù)y=x.以此為評判標(biāo)準(zhǔn),兩種工況下數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型對間歇因子的預(yù)測均取得了較好的效果.由于沒有參與DNN 訓(xùn)練,用作外推測試的工況Re=2.37 × 106間歇因子預(yù)測結(jié)果在線性函數(shù)y=x附近離散程度較高,這與第2.3 節(jié)中測試集損失函數(shù)高于訓(xùn)練集損失函數(shù)的原因一致.由圖10 可知,兩種工況下近壁區(qū)域間歇因子散點(diǎn),均較為均勻地分布在線性函數(shù)y=x兩側(cè),說明兩種工況下間歇因子預(yù)測值雖然存在不同程度的離散,但與計(jì)算值變化規(guī)律基本一致,具有良好的預(yù)測效果.

運(yùn)用數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型對Re=2.73 × 106,Re=2.37 × 106工況下全流場的間歇因子γ進(jìn)行預(yù)測,并與SSTk-ω湍流模型耦合,計(jì)算NLF(2)-0415 后掠翼上翼面橫流轉(zhuǎn)捩.定義Δγ為數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型間歇因子預(yù)測值減去γ-Reθ轉(zhuǎn)捩模型間歇因子計(jì)算值,用于表征數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型計(jì)算結(jié)果在空間分布上相較于γ-Reθ轉(zhuǎn)捩模型的偏差.圖11給出了上翼面近壁區(qū)域Δγ的空間分布.

圖11 近壁區(qū)域Δγ 分布Fig.11 Distribution of Δγ near the wall

觀察圖11,在轉(zhuǎn)捩區(qū)域(即圖11 中ZONEtr)貼近壁面處,少部分網(wǎng)格點(diǎn)Δγ大于0,說明數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型存在將近壁區(qū)域間歇因子低值點(diǎn)預(yù)測出高值的問題.其中工況Re=2.37 × 106的數(shù)據(jù)未參與DNN 訓(xùn)練,數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型預(yù)測出現(xiàn)誤差的網(wǎng)格點(diǎn)數(shù)量多于工況Re=2.73 × 106.

以上翼面的Cf曲線作為數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型預(yù)測精度的判斷標(biāo)準(zhǔn),兩種工況下耦合測試的結(jié)果如圖12 所示.

由圖12 可以看出,兩種工況下數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型計(jì)算得到的Cf上翼面分布曲線與γ-Reθ轉(zhuǎn)捩模型計(jì)算結(jié)果吻合較好,其中數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型計(jì)算的轉(zhuǎn)捩位置略微提前.

圖12 數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型與γ-Reθ 轉(zhuǎn)捩模型Cf 計(jì)算結(jié)果Fig.12 Cf of data driven transition model and γ-Reθ transition model

對數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型預(yù)測轉(zhuǎn)捩位置出現(xiàn)提前的問題進(jìn)行分析.由圖12 可知,工況Re=2.37 ×106下,外插測試的轉(zhuǎn)捩位置提前幅度偏大.前文已經(jīng)對間歇因子預(yù)測值的誤差進(jìn)行了分析,將近壁區(qū)域內(nèi)少量間歇因子低值點(diǎn)預(yù)測出較高的值,會導(dǎo)致渦黏性系數(shù)μt的高值區(qū)域范圍增大,使得μt在網(wǎng)格流向與壁面法向上的增長提前,從而導(dǎo)致預(yù)測的轉(zhuǎn)捩位置更為靠前.

以工況Re=2.37 × 106下的NLF(2)-0415 后掠翼轉(zhuǎn)捩預(yù)測過程為例,采用CPU 時(shí)間作為參照標(biāo)準(zhǔn),使用本研究的預(yù)測方法,層流流場計(jì)算2 萬步至收斂需要3467 s,數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型輸出間歇因子的時(shí)間忽略不計(jì)(不足10 s),耦合SST 模型計(jì)算5 萬步至收斂需要14796 s,共需要計(jì)算7 萬步,耗時(shí)約18000 s;使用γ-Reθ轉(zhuǎn)捩模型與SST 模型耦合計(jì)算5 萬步至收斂,耗時(shí)22660 s.對比可知,本研究提出的轉(zhuǎn)捩預(yù)測方法可以節(jié)約近20%的CPU 時(shí)間,一定程度提高了計(jì)算效率.

3.2 變后掠角NLF(2)-0415 后掠翼預(yù)測

NLF(2)-0415 后掠翼在后掠角足夠小時(shí),邊界轉(zhuǎn)捩主要為流向轉(zhuǎn)捩,轉(zhuǎn)捩位置靠近機(jī)翼后緣;隨著后掠角的增大,邊界層橫流轉(zhuǎn)捩位置逐漸前移,并在后掠角為55°左右時(shí)橫流轉(zhuǎn)捩位置最靠近機(jī)翼前緣,隨后逐漸后移[43].延用圖2 的計(jì)算網(wǎng)格,保持其他工況不變,運(yùn)用數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型預(yù)測后掠角為45°～65°,Re=2.73 × 106工況下NLF(2)-0415 后掠翼上翼面的橫流轉(zhuǎn)捩過程,上翼面流線如圖13 所示.

圖13 上翼面流線Fig.13 Streamline of upper wing

表3 給出了數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型與Zlow 平臺γ-Reθt-CF-模型對轉(zhuǎn)捩位置的計(jì)算結(jié)果[43]的對比.

表3 數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型與γ-Reθt-CF 轉(zhuǎn)捩模型對轉(zhuǎn)捩位置的計(jì)算結(jié)果對比Table 3 Comparison of xtr between data driven transition model and γ-Reθt-CF transition model

由表3 可知,數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型在不同后掠角工況下計(jì)算得到的轉(zhuǎn)捩位置與γ-Reθt-CF 轉(zhuǎn)捩模型基本一致,證明本研究的數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型對改變來流條件的橫流轉(zhuǎn)捩預(yù)測具有一定的泛化能力.

后掠角從45°向65°增大,數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型計(jì)算的轉(zhuǎn)捩位置出現(xiàn)了先向機(jī)翼前緣移動,再向機(jī)翼后緣移動的現(xiàn)象,其中后掠角55°時(shí)轉(zhuǎn)捩位置最為靠前,與文獻(xiàn)[43]給出的轉(zhuǎn)捩位置變化規(guī)律完全一致.相較于后掠角45°,后掠角65°下轉(zhuǎn)捩位置略微靠前,盡管區(qū)別很小,但數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型依然準(zhǔn)確識別了這個(gè)細(xì)微的差異.

以上分析表明,僅用NLF(2)-0415 后掠翼45°后掠角計(jì)算數(shù)據(jù)訓(xùn)練得到的數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型能夠用于不同后掠角下的橫流轉(zhuǎn)捩預(yù)測,具備一定的泛化能力,下一步將針對不同外形的橫流轉(zhuǎn)捩進(jìn)行模型測試.

3.3 6 : 1 標(biāo)準(zhǔn)橢球體預(yù)測

6 : 1 標(biāo)準(zhǔn)橢球體是典型低速三維轉(zhuǎn)捩標(biāo)模.在迎角為15°,單位雷諾數(shù)為6.5 × 106/m,馬赫數(shù)為0.136 的工況下,橢球體表面邊界層轉(zhuǎn)捩由橫流轉(zhuǎn)捩主導(dǎo)[42],具體計(jì)算工況見表4.現(xiàn)將數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型用于預(yù)測6 : 1 標(biāo)準(zhǔn)橢球體表面的轉(zhuǎn)捩過程,以進(jìn)一步校驗(yàn)其泛化能力.

表4 橢球體計(jì)算工況Table 4 Conditions of ellipsoid

本文共對橢球體劃分3 套計(jì)算網(wǎng)格,如圖14所示,網(wǎng)格量分別為43 萬、180 萬和330 萬,在保證網(wǎng)格法向j上第1 層網(wǎng)格y+≈3 的前提下,通過改變流向i、法向j、周向k的網(wǎng)格數(shù)量、增長率來改變網(wǎng)格總量,分別命名為Mesh A,Mesh B,Mesh C.

圖14 橢球體計(jì)算網(wǎng)格Fig.14 Computing mesh for ellipsoid

圖15 給出了本文3 套網(wǎng)格下,橢球體上下表面交界線上,層流流場、湍流流場和轉(zhuǎn)捩流場的Cf分布.

由圖15 可知,對于層流狀態(tài)和湍流狀態(tài),3 套網(wǎng)格計(jì)算得到的Cf分布并無明顯差異,但求解γ-Reθ轉(zhuǎn)捩模型后,高分辨率Mesh B,Mesh C 計(jì)算得到的Cf分布基本重合,與低分辨率Mesh A 存在較大差異.以上分析表明,對于層流與湍流流場的計(jì)算,低分辨率的Mesh A 已經(jīng)滿足網(wǎng)格無關(guān)性,但運(yùn)用γ-Reθ轉(zhuǎn)捩模型求解轉(zhuǎn)捩流場,要求網(wǎng)格分辨率達(dá)到Mesh B 以上.由于數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型輸入數(shù)據(jù)為層流流場物理量,輸出量γ直接與湍流模型耦合,因此層流流場與湍流流場計(jì)算結(jié)果準(zhǔn)確與否,對數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型預(yù)測結(jié)果影響較大.

圖15 不同流動狀態(tài)的Cf 分布Fig.15 Cf distribution in different flow states

現(xiàn)采用低分辨率Mesh A 計(jì)算層流流場,并提取層流流場物理量輸入數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型,用于預(yù)測間歇因子γ.圖16 給出了y+≈30 處數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型預(yù)測的間歇因子分布,并與中等分辨率Mesh B 計(jì)算的轉(zhuǎn)捩流場間歇因子分布進(jìn)行對比.

觀察圖16,數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型對橢球體迎風(fēng)面橫流轉(zhuǎn)捩的起始位置預(yù)測略微提前,但仍具有較高預(yù)測精度.在2 次預(yù)測出橫流轉(zhuǎn)捩后,均發(fā)生了不同程度的再層流化.造成這個(gè)問題的可能原因如下.

圖16 y+≈30 處間歇因子分布Fig.16 Distribution of γ in y+≈30

DNN 僅依靠單一網(wǎng)格點(diǎn)的輸入物理量對該點(diǎn)間歇因子進(jìn)行孤立預(yù)測,在預(yù)測出大范圍轉(zhuǎn)捩后,無法進(jìn)一步考慮擾動在流向上的傳遞、疊加,進(jìn)而修正轉(zhuǎn)捩區(qū)域后側(cè)流場間歇因子的預(yù)測值.

基于上述推論,本研究對數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型補(bǔ)充了考慮擾動傳遞、疊加的算法限制,使其輸出的間歇因子預(yù)測值在自由來流方向上不會錯(cuò)誤地產(chǎn)生衰減,如圖17 所示.

圖17 算法限制原理Fig.17 The theory of algorithm

首先利用數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型對流場中所有網(wǎng)格點(diǎn)的間歇因子γ進(jìn)行預(yù)測.在自由來流方向,即圖17中MN指向上,若某網(wǎng)格點(diǎn)的間歇因子γ預(yù)測值低于前側(cè)距離最近的網(wǎng)格點(diǎn),則使用前側(cè)網(wǎng)格點(diǎn)的γ預(yù)測值對其重新賦值,從而保證自由來流方向上間歇因子γ預(yù)測值始終遞增.考慮到橢球體最前端存在少量γ接近1 的高值點(diǎn),直接啟用此算法會導(dǎo)致大量網(wǎng)格點(diǎn)間歇因子被錯(cuò)誤地賦為高值,因此算法僅對x/a＞0.05 (橢球體長軸a=2.4 m,即x＞0.12 m)的網(wǎng)格點(diǎn)預(yù)測值進(jìn)行限制.圖18 給出了補(bǔ)充算法限制后y+≈30 處間歇因子預(yù)測值.

圖18 補(bǔ)充算法限制后間歇因子分布Fig.18 Distribution of γcorrect

補(bǔ)充算法限制后的間歇因子分布與圖16(b)基本一致,具有較高的預(yù)測精度.繼續(xù)采用低分辨率Mesh A,將數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型預(yù)測的間歇因子γ以及補(bǔ)充算法限制后的間歇因子γcorrect分別與SSTk-ω湍流模型耦合,計(jì)算橢球體表面轉(zhuǎn)捩過程,并與本研究使用的γ-Reθ轉(zhuǎn)捩模型、Chant 平臺橫流模型[35]、BCDF-CF 代碼[34]、OVERFLOW-CF 代碼[34]計(jì)算結(jié)果、實(shí)驗(yàn)結(jié)果[44]對比,如圖19 所示.

圖19 不同算法對Cf 的預(yù)測結(jié)果以及實(shí)驗(yàn)結(jié)果Fig.19 Cf of different methods and test

觀察圖19(a),對數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型補(bǔ)充限制修正前,橢球體表面Cf分布依然能夠反映出轉(zhuǎn)捩陣面的形狀,但轉(zhuǎn)捩與湍流區(qū)域Cf峰值較低,且第2 次預(yù)測出橫流轉(zhuǎn)捩后發(fā)生的再層流化較為明顯.觀察圖19(b),對數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型補(bǔ)充算法限制后,對轉(zhuǎn)捩陣面形狀的預(yù)測能力與本研究使用的γ-Reθ轉(zhuǎn)捩模型、Chant 平臺橫流轉(zhuǎn)捩模型[35]、Langtry 等[34]給出的兩套代碼相仿,但在轉(zhuǎn)捩位置預(yù)測上依然存在一定程度的提前;層流區(qū)域Cf谷值約為0.001,轉(zhuǎn)捩與湍流區(qū)域Cf峰值大于0.003,與本研究使用的γ-Reθ轉(zhuǎn)捩模型、Chant 平臺橫流轉(zhuǎn)捩模型[35]、Langtry等[34]使用的兩套代碼計(jì)算結(jié)果基本一致,Cf值均略低于實(shí)驗(yàn)結(jié)果[44].

以上結(jié)果表明,數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型在顯著降低計(jì)算成本的同時(shí),具有良好的橫流轉(zhuǎn)捩預(yù)測能力,實(shí)現(xiàn)了對橢球體表面橫流轉(zhuǎn)捩的預(yù)測.

4 結(jié)論

本文運(yùn)用數(shù)據(jù)驅(qū)動對三維橫流轉(zhuǎn)捩進(jìn)行了預(yù)測,結(jié)論如下.

(1) 選取層流流場轉(zhuǎn)捩相關(guān)物理量對DNN 進(jìn)行訓(xùn)練,獲得的數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型能夠預(yù)測變雷諾數(shù)的NLF(2)-0415 后掠翼上翼面橫流轉(zhuǎn)捩過程,預(yù)測轉(zhuǎn)捩位置有小幅度提前,具有較高預(yù)測精度.

(2) 由變雷諾數(shù)的NLF(2)-0415 后掠翼計(jì)算數(shù)據(jù)訓(xùn)練得到的數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型具有較好的泛化能力,可以預(yù)測變后掠角NLF(2)-0415 后掠翼上翼面的轉(zhuǎn)捩位置,以及橢球體的三維邊界層橫流轉(zhuǎn)捩.

(3) 數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型與湍流模型耦合的轉(zhuǎn)捩預(yù)測方法,需保證層流、湍流流場的計(jì)算精度,在減少整個(gè)計(jì)算過程所需求解輸運(yùn)方程數(shù)量的同時(shí),進(jìn)一步降低了對網(wǎng)格分辨率的要求,實(shí)現(xiàn)高效轉(zhuǎn)捩預(yù)測.使用同一套計(jì)算網(wǎng)格,相比于γ-Reθ轉(zhuǎn)捩模型,數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型能夠節(jié)省近20%的計(jì)算時(shí)間.

各種因素與邊界層轉(zhuǎn)捩間具有復(fù)雜的影響機(jī)理,如提高粗糙度會導(dǎo)致橫流轉(zhuǎn)捩起始位置前移,如流動分離也會引發(fā)邊界層轉(zhuǎn)捩.針對某個(gè)表面粗糙度設(shè)計(jì),僅具備較強(qiáng)橫流轉(zhuǎn)捩模態(tài)識別能力的數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型,難以完全滿足工程應(yīng)用或科研工作的需求.在下一步工作中,將針對變粗糙度、多轉(zhuǎn)捩模態(tài)下間歇因子的預(yù)測展開研究工作,嘗試在保證數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型橫流轉(zhuǎn)捩預(yù)測精度的前提下,實(shí)現(xiàn)對包括橫流轉(zhuǎn)捩在內(nèi)的多轉(zhuǎn)捩模態(tài)識別,并適用于不同加工精度的表面,提高數(shù)據(jù)驅(qū)動轉(zhuǎn)捩模型的應(yīng)用范圍.