李金瑛, 戴玉婷,2, 楊 超

(1. 北京航空航天大學航空科學與工程學院, 北京 100083; 2. 天目山實驗室, 浙江杭州 310023)

引 言

動態(tài)失速是機翼做大幅度俯仰運動時產生的非線性、 非定常氣動現(xiàn)象。低Reynolds數(shù)下, 飛行器大幅俯仰運動引起分離, 誘導轉捩, 并可能進一步誘發(fā)失速顫振等氣動彈性失穩(wěn)現(xiàn)象[1-2]。國內外對動態(tài)失速開展了大量實驗與仿真研究[3-5]。Spentzos等[6]采用非定常Reynolds平均N-S(unsteady Reynolds-averaged Navier-Stokes, URANS)方程及k-ωSST二方程湍流模型對三維機翼進行了數(shù)值仿真計算。Wang等[7]嘗試了將URANS與分離渦(detached eddy simulation, DES)方法結合用于動態(tài)失速仿真。Kim等[8]采用了大渦模擬(large-eddy simulation, LES)方法對NACA0012翼型進行了動態(tài)失速仿真。

傳統(tǒng)湍流模型包括Spalart-Allmaras一方程湍流模型和k-ε、k-ω、k-ωSST二方程湍流模型等。其中, Menter[9]提出的k-ωSST模型具有良好的綜合性能, 但無法對轉捩過程進行準確計算?，F(xiàn)有轉捩過程計算方法主要有3種: 基于線性穩(wěn)定性理論的eN方法[10-12], 采用低Reynolds數(shù)湍流模型的仿真方法[13], 及采用引入轉捩間歇因子(γ)的湍流模型。其中, eN方法僅適用于自然轉捩的計算; 低Reynolds數(shù)湍流模型忽略了轉捩區(qū)物理特性, 易造成誤差[14]。采用轉捩間歇因子控制轉捩發(fā)展是目前主流的計算方式。國內外發(fā)展了多種γ的計算方法及湍流模型[15-18]: Dhawan等[19]提出了γ的代數(shù)計算方法。Suzen等[20]結合前人工作[21-22]提出了耦合RANS的γ輸運方程。Menter等[14]提出了γ-Reθ四方程湍流模型, 并將其簡化為SST-γ三方程湍流模型[23], 該模型在RANS方程中對轉捩過程模擬精度較高, 但增加了一個輸運微分方程, 計算效率低于k-ωSST二方程模型。目前, 須發(fā)展同時具備高精度與高效率的轉捩湍流模型。

利用數(shù)據(jù)驅動建模可以提高湍流模擬的效率。當前數(shù)據(jù)驅動建模方法包括直接辨識氣動力系數(shù)的黑箱模型、 耦合求解湍流方程的灰箱模型及辨識優(yōu)化湍流方程參數(shù)的白箱模型。在白箱模型方面, Singh等[24-25]使用機器學習修正了Spalart-Allmaras湍流模型參數(shù); Yang等[26]采用隨機森林與人工神經(jīng)網(wǎng)絡對k-ω-γ-Ar四方程湍流模型參數(shù)進行了修正; Zafar等[27]基于遞推神經(jīng)網(wǎng)絡與卷積神經(jīng)網(wǎng)絡建立了具有高泛化能力的轉捩湍流模型。此外, 國內外學者通過將神經(jīng)網(wǎng)絡與CFD湍流方程耦合求解[28-30]建立了多種灰箱模型。

本文基于SST-γ三方程湍流模型的流場計算結果, 訓練深度神經(jīng)網(wǎng)絡模型預測轉捩間歇因子, 建立流場參數(shù)與γ的關系, 據(jù)此修正k-ωSST二方程湍流模型, 建立將深度神經(jīng)網(wǎng)絡嵌入CFD工具的SST-γ-machine learning(下文簡稱為SST-γ-ML)耦合二方程湍流模型。利用該耦合二方程湍流模型對NACA0012翼型進行動態(tài)失速流場仿真。

1 耦合二方程湍流模型

1.1 耦合湍流模型建模

本節(jié)基于k-ωSST與數(shù)據(jù)驅動的轉捩間歇因子模型建立SST-γ-ML耦合二方程湍流模型, 在不增加輸運微分方程個數(shù)的前提下, 使耦合二方程湍流模型具備模擬轉捩過程的能力。

SST-γ-ML耦合二方程湍流模型的控制方程為

其中,S為Reynolds應力,Ω為渦量幅值,μ為動力黏度,μt為湍流黏度。其他常數(shù)與函數(shù)為

其中,ReT為湍流Reynolds數(shù)

(1)

控制方程組由k與ω的輸運微分方程及γ的神經(jīng)網(wǎng)絡代數(shù)方程組成。此方程組中, 前兩個輸運方程即構成由γ修正的k-ωSST二方程湍流模型。γ在SST-γ三方程模型中由輸運微分方程求出, 說明其與流場中U、 壓力p、k、ω等流場參數(shù)具有復雜的非線性關系, 因此, 采用神經(jīng)網(wǎng)絡代替輸運微分方程預測γ, 可提高計算效率, 該神經(jīng)網(wǎng)絡具體結構在1.2節(jié)詳細介紹。

SST-γ-ML耦合二方程湍流模型的訓練與嵌入策略如圖1所示。紅色虛線框內表示了訓練過程, 選取對轉捩過程呈現(xiàn)相對準確的SST-γ三方程模型計算的流場參數(shù)和轉捩間歇因子作為神經(jīng)網(wǎng)絡模型的輸入輸出數(shù)據(jù)集。

圖1 模型的訓練與使用Fig. 1 Training and use of model

神經(jīng)網(wǎng)絡與CFD平臺的耦合方式有松耦合與緊耦合。松耦合在平臺間傳遞流場數(shù)據(jù), 緊耦合將神經(jīng)網(wǎng)絡模型直接嵌入CFD平臺。對非定常流動, 如動態(tài)失速過程, 每個時間步間流場變化大, 須在每個時間步對CFD仿真計算與神經(jīng)網(wǎng)絡模型進行交互, 松耦合并不適用。本文采用緊耦合方式建立數(shù)據(jù)驅動的轉捩模型, 如圖1綠色線框所示。

基于SST-γ-ML耦合二方程湍流模型建模主要分5個步驟, 包括訓練神經(jīng)網(wǎng)絡、 凍結導出、 嵌入耦合模型、 模型編譯與CFD仿真計算。為將訓練完成的神經(jīng)網(wǎng)絡完整嵌入OpenFOAM開源CFD軟件中, 使用protocol buffer(pb)格式文件導出神經(jīng)網(wǎng)絡模型。

在CFD計算中對每個流體網(wǎng)格單元均獨立調用神經(jīng)網(wǎng)絡模型, 預測每個流體網(wǎng)格單元的轉捩間歇因子。雖然該建模方式計算效率低于直接預測全流場轉捩間歇因子, 但大幅簡化了神經(jīng)網(wǎng)絡結構, 訓練集數(shù)據(jù)量降低2～3個數(shù)量級。此外, 它對流場不同網(wǎng)格單元具有泛化性, 能夠直觀反映出流場參數(shù)與不同位置γ之間的關系。

1.2 神經(jīng)網(wǎng)絡模型和輸入?yún)?shù)設計

本文基于反向傳播的全連接結構人工神經(jīng)網(wǎng)絡建立流場參數(shù)與轉捩間歇因子之間的預測模型。反向傳播(back-propagation, BP)神經(jīng)網(wǎng)絡的參數(shù)及樣本數(shù)據(jù)預處理方法具體如下:

采用Adam算法作為誤差反向傳播的優(yōu)化算法, 均方誤差(mean squared error, MSE)作為誤差計算算法; 神經(jīng)網(wǎng)絡隱藏層層數(shù)為4層, 每層分別含有[40, 20, 10, 10]個神經(jīng)元, 批處理數(shù)為256, 學習率為0.001, 衰減率為0.001; 在輸入層到第1層隱藏層間加入批歸一化層, 從第4層隱藏層到輸出層選擇Sigmoid作為激活函數(shù), 以控制轉捩間歇因子輸出范圍為0～1。隱藏層間激活函數(shù)選取Tanh以擴大近壁面數(shù)據(jù)梯度; 訓練前對樣本數(shù)據(jù)進行打亂預處理, 并采用多個時間步流場參數(shù)作為訓練集數(shù)據(jù)。

通過對流場參數(shù)與轉捩間歇因子關系的物理分析及優(yōu)化設計, 選取神經(jīng)網(wǎng)絡模型的流場輸入?yún)?shù)組合為

(2)

其中,Ux為平行來流方向速度,Uy為垂直來流方向速度,dω為網(wǎng)格單元與壁面的距離。

6個輸入?yún)?shù)中,p,μt/ρ,Ux,Uy直接從流場中提取。k/ω與湍流Reynolds數(shù)ReT呈正相關(式(1)),k表征湍流運動分量具有的動能

流動從層流轉捩到湍流的過程中, 湍流動能逐漸增大, 與壁面間的剪切應力逐漸減小, 對應k增大,ω減小, 即輸入?yún)?shù)k/ω增大。

其中, dUy/dy為垂直來流方向速度導數(shù)。比例系數(shù)λθL表征流場壓力梯度, 壁面距離平方項與比例系數(shù)呈線性關系, 說明轉捩程度對壁面距離非常敏感。在接近壁面的位置, 極小距離下會發(fā)生顯著的轉捩發(fā)展, 故采用平方項輸入?yún)?shù)。

2 數(shù)據(jù)驅動的轉捩模型驗證

以標準T3A平板為研究對象, 在Re=6.12×105下進行穩(wěn)態(tài)轉捩仿真計算, 驗證數(shù)據(jù)驅動的耦合二方程湍流模型的有效性。

T3A平板長度為1.7 m, 流場高度1 m, 平板前緣半徑為0.75 mm, 流場主要區(qū)域為平板前區(qū)域(0～0.04 m)、 前緣區(qū)域與平板區(qū)域(0.04～1.70 m), 平板流場網(wǎng)格如圖2所示。在平板前緣附近進行網(wǎng)格加密, 網(wǎng)格量為26 820, 網(wǎng)格平均Y+為0.53。平板表面(wall)采用無滑移邊界, 流場頂部及平板前區(qū)域的流場底部(above)采用滑移邊界, 遠場設置自由來流速度U∞、 湍流度Tu∞與黏度比Rμ=μt/μ, 出口壓力為0。

圖2 T3A平板計算網(wǎng)格及邊界條件Fig. 2 Mesh and boundary conditions of T3A flat plate

采用如表1所示入口邊界條件, 分別采用SST-γ-ML耦合二方程湍流模型,k-ωSST二方程湍流模型與SST-γ三方程湍流模型計算T3A標準平板表面湍流度與摩阻系數(shù), 結果如圖3所示。

表1 T3A平板入口條件

(a) Turbulence intensity

結果表明, 耦合二方程湍流模型對T3A平板定常流動條件下的轉捩位置及對應摩阻系數(shù)預測準確, 與實驗值[31]及三方程湍流模型結果吻合較好。耦合二方程湍流模型摩阻系數(shù)預測結果相對三方程湍流模型結果誤差為1%, 驗證了數(shù)據(jù)驅動的SST-γ-ML耦合二方程湍流模型在轉捩預測中的有效性。

3 數(shù)值算例與結果討論

3.1 數(shù)值模型

以NACA0012二維翼型為研究對象, 應用第2節(jié)驗證的SST-γ-ML耦合二方程湍流模型在Re=1.35×105下分別進行靜態(tài)小迎角及動態(tài)失速過程的氣動力計算。

NACA0012二維翼型弦長c=0.15 m, 翼型流場網(wǎng)格如圖4所示。在翼型動態(tài)失速氣動力計算中, 旋轉域與固定域采用滑移網(wǎng)格以耦合不連續(xù)邊界, 適用于旋轉幾何體數(shù)值仿真。將全流場劃分為內部旋轉區(qū)域(紅圈內)與外部固定區(qū)域(紅圈外)。內部區(qū)域為596×100的O形結構網(wǎng)格, 最大Y+為0.93; 外流場網(wǎng)格量為9 779的混合網(wǎng)格。為確保尾跡區(qū)流動的充分發(fā)展, 壓力出口邊界設置在后緣下游20倍弦長位置。

圖4 NACA0012計算網(wǎng)格及邊界條件Fig. 4 Mesh and boundary conditions of NACA0012

遠場設置U∞=14 m/s,Tu∞=0.08%, 采用PIMPILE瞬態(tài)求解器進行壓力速度解耦計算, 計算時間步長固定為1×10-5s。對SST-γ三方程湍流模型, 控制方程中對流項U、k與ω的散度計算均采用線性迎風離散格式,γ的散度計算采用1階迎風離散格式。對SST-γ-ML耦合二方程湍流模型,γ無須進行控制方程計算, 其他參數(shù)離散格式與三方程湍流模型相同。

3.2 穩(wěn)態(tài)小迎角

本節(jié)應用SST-γ-ML耦合二方程湍流模型對NACA0012二維翼型進行小迎角穩(wěn)態(tài)氣動力計算。

應用SST-γ三方程湍流模型與SST-γ-ML耦合二方程模型對2°～8°迎角流場進行計算, 耦合二方程模型中神經(jīng)網(wǎng)絡訓練集迎角及數(shù)值仿真計算測試集迎角如表2所示。

表2 神經(jīng)網(wǎng)絡訓練集與測試集迎角設置

耦合二方程湍流模型與三方程湍流模型在各迎角下升力系數(shù)計算結果如圖5所示。

圖5 不同迎角下升力計算結果Fig. 5 Lift coefficient at different angles of attack

耦合二方程湍流模型對未經(jīng)訓練的α=5°, 6°, 7°內插迎角升力系數(shù)計算準確: 與三方程模型計算結果對比, 平均相對誤差為2.60%; 與實驗結果對比[32], 平均相對誤差為8.81%。未經(jīng)訓練的α=2°外插迎角結果, 對比實驗結果相對誤差為19.63%, 對比三方程結果相對誤差為23.81%。

僅采用4°, 8°迎角數(shù)據(jù)訓練的耦合二方程湍流模型在6°迎角下計算翼型表面壓力系數(shù)曲線如圖6所示。對未經(jīng)訓練的6°迎角情況, 耦合二方程湍流模型預測表面壓力系數(shù)結果與三方程模型準確對應。對4°與8°迎角算例, 耦合二方程湍流模型同樣可以準確預測壓力系數(shù)。表明該數(shù)據(jù)驅動的轉捩模型具有較好的迎角泛化能力。

圖6 6°迎角下壓力系數(shù)預測效果Fig. 6 Prediction of pressure coefficient at α=6°

3.3 動態(tài)失速

本節(jié)應用SST-γ-ML耦合二方程湍流模型對NACA0012翼型進行動態(tài)失速過程計算。設定翼型做給定正弦規(guī)律的俯仰運動

α(t)=10+15sin(18.67t)

其中,α(t)為t時刻翼型迎角。運動初始迎角為10°, 俯仰運動迎角范圍為α∈[-5°, 25°], 俯仰運動周期T=0.337 s。

計算總時長為10T, 取第6～10周期進行相平均, 得到SST-γ三方程湍流模型與SST-γ-ML耦合二方程湍流模型在一個穩(wěn)定周期內計算的升力系數(shù)。實驗[32]、 三方程湍流模型、 LES模型[8]及耦合二方程湍流模型所得升力系數(shù)滯回曲線如圖7所示。

圖7 升力系數(shù)隨迎角變化曲線Fig. 7 Curve of lift coefficient with angle of attack

升力系數(shù)的滯回曲線趨勢與實驗及三方程結果均符合較好, 與三方程模型計算結果相對誤差為11.95%。從曲線上看, 升力系數(shù)曲線在翼型上仰過程的線性段計算準確; 在20°～25°迎角失速區(qū)間, 耦合二方程模型與三方程模型幾乎同時達到第1次升力系數(shù)峰值, 耦合二方程模型計算的第2次升力系數(shù)峰值出現(xiàn)時間略晚于三方程模型; 在下俯過程中, 耦合二方程模型計算所得升力系數(shù)相對三方程結果波動小, 且取值偏低。

結合渦量云圖, 對比翼型在上仰過程中耦合二方程湍流模型與三方程湍流模型對動態(tài)失速過程中典型流動狀態(tài)的仿真結果, 如圖8所示。

圖8 翼型上仰過程典型流動狀態(tài)渦量場: (a) 三方程湍流模型; (b) 耦合二方程湍流模型Fig. 8 Vortex fields of typical flow conditions via upstroke: (a) 3 equation turbulence model; (b) coupled 2 equation turbulence model

在翼型上仰過程中,α∈[-5°, 14°]迎角下, 耦合二方程模型對附著在翼型上表面的湍流邊界層及后緣脫落渦仿真準確(圖8(a-1; b-1)); 隨迎角增大, 耦合二方程模型與三方程模型幾乎同時出現(xiàn)前緣渦并迅速增長(圖8(a-2, 3; b-2, 3))。前3個典型流動狀態(tài)上, 耦合二方程模型數(shù)值仿真結果與三方程模型均準確對應, 該現(xiàn)象解釋了升力系數(shù)曲線(圖7)在上仰過程的α∈[-5°, 20°]迎角區(qū)間耦合二方程與三方程模型曲線的準確對應, 耦合二方程模型曲線的第1個升力系數(shù)峰與三方程曲線同時出現(xiàn)。在α∈[20°, 25°]區(qū)間, 前緣渦持續(xù)生成并脫落, 耦合模型對此期間前緣渦生成與脫落規(guī)律、 流動狀態(tài)與渦量大小等實現(xiàn)了準確仿真, 但速率略慢于三方程仿真結果(圖8(a-4; b-4)), 該速率差解釋了升力系數(shù)曲線(圖7)的α∈[20°, 25°]上仰區(qū)間里, 第2個升力系數(shù)峰的出現(xiàn)時間略晚于三方程曲線的誤差原因?？傮w而言, 耦合二方程模型對翼型上仰過程仿真較為準確。

圖9對比了翼型下俯過程中兩個湍流模型對渦量場的仿真情況。

圖9 翼型下俯過程典型流動狀態(tài)渦量場: (a) 三方程湍流模型; (b) 耦合二方程湍流模型Fig. 9 Vortex fields of typical flow conditions via downstroke: (a) 3 equation turbulence model; (b) coupled 2 equation turbulence model

在翼型下俯過程中, 交替發(fā)生前緣渦與后緣渦的脫落(圖9(a-1, 2, 3; b-1, 2, 3)); 在小迎角下, 上翼面尾緣流動再次附著(圖9(a-4; b-4)); 此后, 邊界層均附著在翼型表面(圖9(a-5; b-5))。耦合二方程模型在下俯過程中對典型流動狀態(tài)的仿真結果與三方程均能準確對應。

4 結論

本文建立了由湍動能k、 單位湍動能耗散率ω的輸運方程及預測轉捩間歇因子γ的神經(jīng)網(wǎng)絡構成的耦合二方程湍流模型(SST-γ-ML模型), 用于翼型動態(tài)失速過程的高精度方程。

1) 利用基于數(shù)據(jù)驅動的轉捩間歇因子預測模型對T3A平板進行了流場仿真, 耦合二方程湍流模型對包含轉捩過程的平板表面預測摩阻系數(shù)相對于SST-γ三方程模型結果平均誤差為1%;

2) 對NACA0012翼型進行了低Reynolds數(shù)動態(tài)失速過程的流場仿真, 耦合二方程湍流模型預測的非定常氣動升力與SST-γ三方程模型結果相比, 誤差小于12%。