張亦知,程誠,范釔彤,李高華,李偉鵬
上海交通大學 航空航天學院,上海 200240
湍流邊界層廣泛存在于航空、航天、建筑、生物等領域,其主要的計算方法包括直接數(shù)值模擬(DNS)、大渦模擬(LES)和雷諾平均Navier-Stokes(RANS)模擬。直接數(shù)值模擬求解非定常流動中的湍流擾動,大渦模擬湍流小尺度結構而分辨湍流大尺度結構,具有較好的預測精度,但DNS和LES所需的計算網(wǎng)格規(guī)模與雷諾數(shù)的指數(shù)冪成正比,在高雷諾數(shù)下的計算成本巨大,使得利用DNS或LES開展高雷諾數(shù)流動仿真存在困難[1]。RANS通過引入湍流模型方程模化所有的時空多尺度湍流脈動,通過顯式表達式[2]或偏微分方程[3]求得雷諾應力,其計算量與雷諾數(shù)呈弱相關關系,計算效率高,在工程中得到廣泛應用。
然而,RANS的預測精度與DNS或LES存在差異,且在實際應用中湍流模型的選擇依賴于設計者的經(jīng)驗,通過調整或優(yōu)化湍流模型中的個別系數(shù)難以提高其預測精度,且不具有普適性。RANS計算最大的不確定性來源于湍流模型方程的固有形式[4-5],近年來,利用數(shù)據(jù)挖掘和機器學習方法修正RANS湍流模型顯示出巨大的潛力和前景[6-7]。
Duraisamy等[8-9]利用流場反演獲得湍流模型修正項,進而使用神經(jīng)網(wǎng)絡建立流場特征與修正項之間的映射關系,獲得嵌入式的湍流模型,經(jīng)重構后湍流模型計算結果更接近實驗數(shù)據(jù)。Ling和Templeton[10]利用DNS和LES的計算結果,發(fā)展了機器學習修正的RANS湍流模型,比較了不同機器學習方法,如支持向量機、決策樹以及隨機森林對流場不確定性的預測效果,結果表明這些方法在與訓練樣本相異的流場中均能準確地捕捉流動特征。Xiao等[11]使用貝葉斯方法量化RANS計算中模型形式的不確定性,并加入經(jīng)驗知識約束,結果表明該方法利用稀疏的先驗數(shù)據(jù)能較好地捕獲基準模型的誤差,獲得更準確的后驗平均速度分布。Wang等[12]的研究結果表明,雷諾應力的預測能力是影響RANS計算精度的重要因素,由此建立了基于雷諾應力張量的隨機森林機器學習重構方法,并針對湍流方管及典型分離流動進行了模型重構。Wu等[13]采用相似的方法構造了隨機森林模型,給出了RANS計算置信度的先驗估計。Zhang等[14]在翼型的近壁區(qū)、尾跡區(qū)等不同區(qū)域分別構造了徑向基神經(jīng)網(wǎng)絡,?;烁呃字Z數(shù)湍流邊界層中的渦黏項,并將其實時地引入流場求解過程,提高了計算準確度,且在翼型或流動條件變化時具有一定通用性。在湍流計算領域,基于大數(shù)據(jù)的機器學習方法具有廣闊的應用前景,文獻[15]對數(shù)據(jù)驅動式湍流建模理論和方法進行了綜述分析。
數(shù)據(jù)驅動式建模方法依賴于樣本數(shù)據(jù)的個數(shù)和質量,小樣本或樣本冗余誤差會導致模型的失真或失效。為克服數(shù)據(jù)驅動式建模對數(shù)據(jù)樣本個數(shù)和質量的依賴關系,需要引入基于物理知識約束的先驗條件,在建模過程中約束模型的邊界和搜索空間,有利于避免非物理現(xiàn)象的產(chǎn)生,提高模型的預測精度。近年來,基于物理知識的數(shù)據(jù)驅動式建模成為研究熱點,Raissi等[16]以自由振動的圓柱為研究對象,將有限的速度散點數(shù)據(jù)作為訓練樣本,利用神經(jīng)網(wǎng)絡重構的過程中,引入流體控制方程作為先驗的物理知識約束,對流動參數(shù)進行了預測,獲得準確的升/阻力預測結果。Raissi等[17]構造了基于物理信息的神經(jīng)網(wǎng)絡,用于求解非線性偏微分方程,利用少量的訓練數(shù)據(jù),成功捕捉了系統(tǒng)的非線性特征。針對湍流模擬的RANS方程,尚未發(fā)現(xiàn)基于物理知識約束的數(shù)據(jù)驅動式建模方法。
準確預測湍流摩擦阻力具有重要的工程價值,本文作者在利用DNS數(shù)據(jù)研究湍流摩擦阻力的分解中,發(fā)現(xiàn)湍流摩擦阻力的分解項在湍流邊界層的內層具有統(tǒng)一的比擬關系[18],與雷諾數(shù)大小無關,與可壓縮性無關[19],并在槽道和平板湍流邊界層中得到了驗證[18-20],結果表明該比擬關系是一種物理的、普適的規(guī)律,是一種先驗的物理知識。本文基于該物理知識開展了針對S-A湍流模型的數(shù)據(jù)驅動式建模及修正,在建模過程中引入物理知識約束,對比了有/無物理知識約束的修正對槽道湍流摩擦阻力預測精度的影響。
為封閉雷諾平均的Navier-Stokes方程,Spalart和Allmaras[21]構造了一方程湍流模型方程,通過求解標量運輸方程來計算湍流渦黏系數(shù)并獲得雷諾應力,具體形式為
(1)
(2)
(3)
(4)
固壁與流體相對運動而產(chǎn)生摩擦,摩擦阻力可以表征為固壁對流體做功的一種形式,對于槽道湍流而言,摩擦阻力做功,一部分能量通過黏性直接耗散轉換成熱能,一部分能量則通過湍動能生成項維持湍流態(tài)。基于這一思想,Renard和Deck[23]推導了湍流摩擦阻力系數(shù)Cf的分解公式:
(5)
(6)
式中:a1=2.644;b1=1.895;c1=-0.805;a2=1.777;b2=1.118;c2=1.642。式(6)中僅有一個自變量y+,圖1表明擬合結果與DNS數(shù)據(jù)的符合良好。本文將式(6)作為一種先驗物理知識約束引入到數(shù)據(jù)驅動建模中。
圖1 不同雷諾數(shù)下直接黏性耗散項的法向分布Fig.1 Direct viscous dissipation curves along wall-normal direction under different Reynolds numbers
保留S-A湍流模型方程的基本形式,對生成項P添加修正因子β,修正后的湍流模型方程為
(7)
修正因子β對流場的影響是全局性的,其取值隨著法向高度的變化而變化。本文通過添加修正因子β來校正S-A湍流模型,其本質是在S-A湍流模型中添加了源項的校正量δ,即δ=(β-1)P,其大小與生成項相關。初始時β=1,校正量δ均為0,式(7)即退化為原始的S-A方程。通過引入修正因子β,擬使得S-A湍流模型預測更準確,即校正后的S-A湍流模型計算結果更接近于真實解,真實解可取高精度的DNS、LES結果,也可利用試驗測量的數(shù)據(jù)進行標定,本文中的真實解為DNS計算結果[24]。
為獲得準確可靠的修正因子,需要設定合理的目標函數(shù),結合高效的優(yōu)化方法,快速獲得修正因子分布。傳統(tǒng)方法對模型中系數(shù)進行調整優(yōu)化,例如對馮卡門常數(shù)(κ)或湍流普朗特數(shù)(σ)的校正[24-25],但其難以克服重要流場物理特性(如雷諾數(shù))的影響。由此,考慮到湍流摩擦阻力的預測與平均流向速度剖面緊密相關,且引入式(6)作為物理知識約束,設定目標函數(shù)為
(8)
式中:(Uj-Uj,DNS)2表征對平均流向速度剖面的逼近程度;(Ωj-Ω(y+))2表征對近壁黏性耗散項分布的逼近程度;Nc為網(wǎng)格單元數(shù)。值得指出的是,修正因子β與直接黏性耗散和平均流向速度之間不存在直接的量化關系,僅存在隱式關聯(lián)與映射。
為對比分析有/無物理知識約束對湍流摩擦阻力的預測誤差,本文也采用了基于平均流向速度的目標函數(shù):
(9)
由于設計變量個數(shù)與網(wǎng)格數(shù)量相當,使用有限差分法求解梯度計算量巨大。伴隨方法(Adjoint Method)是一種高效的梯度求解方法,其計算規(guī)模與設計變量的數(shù)目基本無關,可大幅減小梯度的求解時間。伴隨方法以偏微分方程系統(tǒng)控制理論為基礎,Jameson[26]首次將伴隨方法應用于氣動設計中,把氣動設計轉換為一個滿足特定約束的最優(yōu)控制問題,梯度求解的計算量約為2倍的流場計算量。伴隨方法包括連續(xù)伴隨和離散伴隨,本文采用離散伴隨方法[27],伴隨方程和梯度的表達式為
(10)
(11)
圖2 湍流模型修正優(yōu)化流程Fig.2 Optimization process of turbulence model correction
以二維方程[29]來驗證伴隨方法所求得的梯度與有限差分法的一致性,圖3對比了有限差分法與伴隨方法所得的梯度分布,其偏差小于0.57%,屬于合理范圍內。
圖3 伴隨方法中的梯度驗證Fig.3 Gradient verification in adjoint method
圖4 直接黏性耗散項的法向分布Fig.4 Distribution of direct viscous dissipation in wall-normal direction
圖5 直接黏性耗散項相對誤差的法向分布Fig.5 Distribution of relative error of direct viscous dissipation in wall-normal direction
圖4和圖5中的比較結果表明,原始S-A湍流模型對近壁區(qū)直接黏性耗散項的預測存在明顯偏差;如果通過數(shù)據(jù)驅動模型僅修正流向速度剖面,對直接黏性耗散項的預測有所改善,但誤差仍然較大;在數(shù)據(jù)驅動建模的目標函數(shù)中,進一步引入關于直接黏性耗散的先驗物理知識作為約束條件,預測精度顯著改善,計算結果與真實結果的相對誤差小于3%。
圖6 平均流向速度的法向分布Fig.6 Distribution of averaged streamwise velocity in wall-normal direction
圖7 平均流向速度相對誤差的法向分布Fig.7 Distribution of relative error of averaged streamwise velocity in wall-normal direction
圖8比較了不同雷諾數(shù)下修正因子β的變化情況,注意原始S-A湍流模型中修正因子恒為1。當采用Ju為目標函數(shù)時,僅根據(jù)流向平均速度進行模型修正,修正因子β的法向變化趨勢基本與圖7中流向平均速度誤差的分布一致,說明修正因子β具有較好的活性,可根據(jù)當?shù)氐钠骄俣日`差對S-A模型中的湍流黏性生成項進行相應調整。當采用Jud為目標函數(shù)時,引入了直接黏性耗散項作為流向平均速度修正的補充,以Reτ=180的槽道為例,對于復合目標函數(shù)修正因子β做出相應的調整,在y+=4.65處達到峰值,說明目標函數(shù)Jud起到了復合校正效果。同時,圖8中β系數(shù)隨流場位置及當?shù)亓鲌鲎兞康淖兓l(fā)生演變,也說明了本文建立的伴隨優(yōu)化方法計算可以有效地建立對方程預測精度的調整和修正。需要注意的是,修正因子β的取值范圍在-10~20之間,說明原始S-A湍流模型中的生成項存在低估或者逆?zhèn)鬟f現(xiàn)象。修正因子的本質作用是在S-A湍流模型中引入δ=(β-1)P作為源項的修正項,該修正項的取值范圍根據(jù)預測結果與真實結果的偏差而定。
圖8 修正因子β的法向分布Fig.8 Distribution of correction factor β in wall-normal direction
湍流摩擦阻力是一個較小的壁面變量,目前針對大型客機的總阻力預測,利用RANS計算要確保阻力誤差在1 count以內仍是巨大的挑戰(zhàn),其中最為困難的是湍流摩擦阻力的準確預測,其約占總阻力的50%以上。本文針對槽道湍流這一簡單基礎構型進行了數(shù)據(jù)驅動的湍流模型修正,在低雷諾數(shù)下其預測精度顯著提高,而在較高雷諾數(shù)條件下,表1中的數(shù)據(jù)顯示無修正的湍流模型對摩擦阻力的預測誤差已經(jīng)小于0.48%,通過模型修正,這個誤差得到了進一步的減小。該結果表明,基于數(shù)據(jù)驅動式的湍流模型修正具有較好的應用前景。
表1 壁面摩擦阻力系數(shù)的相對誤差Table 1 Relative errors of skin friction coefficients
數(shù)據(jù)驅動式建模的本質是數(shù)據(jù)回歸和優(yōu)化問題引入物理知識的先驗約束,一定程度上可降低對數(shù)據(jù)樣本個數(shù)和質量的依賴,提高模型的預測精度。本文提出了一種基于物理知識的數(shù)據(jù)驅動式湍流模型修正方法,并以槽道湍流為研究對象驗證了該方法的優(yōu)越性,獲得的結論包括:
1) 基于湍流摩擦阻力分解的先驗物理知識,建立了一種基于物理知識約束的湍流模型修正方法,將物理知識顯性地引入目標函數(shù)的設定中,結合伴隨優(yōu)化方法,高效率地獲得修正因子的分布,以達到提高預測精度的目的。
2) 以槽道湍流為例,對比了有/無物理知識約束下的湍流模型預測精度,結果表明引入物理知識約束后,可提高對湍流摩擦阻力的預測精度。
本文研究存在如下不足,希望在后期研究中克服:
1) 引入的基于湍流摩擦阻力分解的先驗物理知識,適用于無逆壓梯度的湍流邊界層,而對于有逆壓梯度導致的流動分離問題,該物理知識的有效性需要進一步的驗證,或引入更為普適的物理知識。
2) 僅開展了基于DNS數(shù)據(jù)的先驗驗證,而沒有利用機器學習方法,如神經(jīng)網(wǎng)絡、決策樹等,建立物理變量與目標函數(shù)之間的映射關系,以開展無數(shù)據(jù)條件下的后驗驗證。
3) 建立目標函數(shù)時,用到了DNS平均速度剖面數(shù)據(jù),沒有完全避免對DNS數(shù)據(jù)的依賴,在后期工作中有待進一步完善。