張亦知，程誠(chéng)，范釔彤，李高華，李偉鵬

上海交通大學(xué) 航空航天學(xué)院，上海 200240

湍流邊界層廣泛存在于航空、航天、建筑、生物等領(lǐng)域，其主要的計(jì)算方法包括直接數(shù)值模擬(DNS)、大渦模擬(LES)和雷諾平均Navier-Stokes(RANS)模擬。直接數(shù)值模擬求解非定常流動(dòng)中的湍流擾動(dòng)，大渦模擬湍流小尺度結(jié)構(gòu)而分辨湍流大尺度結(jié)構(gòu)，具有較好的預(yù)測(cè)精度，但DNS和LES所需的計(jì)算網(wǎng)格規(guī)模與雷諾數(shù)的指數(shù)冪成正比，在高雷諾數(shù)下的計(jì)算成本巨大，使得利用DNS或LES開展高雷諾數(shù)流動(dòng)仿真存在困難[1]。RANS通過(guò)引入湍流模型方程模化所有的時(shí)空多尺度湍流脈動(dòng)，通過(guò)顯式表達(dá)式[2]或偏微分方程[3]求得雷諾應(yīng)力，其計(jì)算量與雷諾數(shù)呈弱相關(guān)關(guān)系，計(jì)算效率高，在工程中得到廣泛應(yīng)用。

然而，RANS的預(yù)測(cè)精度與DNS或LES存在差異，且在實(shí)際應(yīng)用中湍流模型的選擇依賴于設(shè)計(jì)者的經(jīng)驗(yàn),通過(guò)調(diào)整或優(yōu)化湍流模型中的個(gè)別系數(shù)難以提高其預(yù)測(cè)精度，且不具有普適性。RANS計(jì)算最大的不確定性來(lái)源于湍流模型方程的固有形式[4-5]，近年來(lái)，利用數(shù)據(jù)挖掘和機(jī)器學(xué)習(xí)方法修正RANS湍流模型顯示出巨大的潛力和前景[6-7]。

Duraisamy等[8-9]利用流場(chǎng)反演獲得湍流模型修正項(xiàng)，進(jìn)而使用神經(jīng)網(wǎng)絡(luò)建立流場(chǎng)特征與修正項(xiàng)之間的映射關(guān)系，獲得嵌入式的湍流模型，經(jīng)重構(gòu)后湍流模型計(jì)算結(jié)果更接近實(shí)驗(yàn)數(shù)據(jù)。Ling和Templeton[10]利用DNS和LES的計(jì)算結(jié)果，發(fā)展了機(jī)器學(xué)習(xí)修正的RANS湍流模型，比較了不同機(jī)器學(xué)習(xí)方法，如支持向量機(jī)、決策樹以及隨機(jī)森林對(duì)流場(chǎng)不確定性的預(yù)測(cè)效果，結(jié)果表明這些方法在與訓(xùn)練樣本相異的流場(chǎng)中均能準(zhǔn)確地捕捉流動(dòng)特征。Xiao等[11]使用貝葉斯方法量化RANS計(jì)算中模型形式的不確定性，并加入經(jīng)驗(yàn)知識(shí)約束，結(jié)果表明該方法利用稀疏的先驗(yàn)數(shù)據(jù)能較好地捕獲基準(zhǔn)模型的誤差，獲得更準(zhǔn)確的后驗(yàn)平均速度分布。Wang等[12]的研究結(jié)果表明，雷諾應(yīng)力的預(yù)測(cè)能力是影響RANS計(jì)算精度的重要因素，由此建立了基于雷諾應(yīng)力張量的隨機(jī)森林機(jī)器學(xué)習(xí)重構(gòu)方法，并針對(duì)湍流方管及典型分離流動(dòng)進(jìn)行了模型重構(gòu)。Wu等[13]采用相似的方法構(gòu)造了隨機(jī)森林模型，給出了RANS計(jì)算置信度的先驗(yàn)估計(jì)。Zhang等[14]在翼型的近壁區(qū)、尾跡區(qū)等不同區(qū)域分別構(gòu)造了徑向基神經(jīng)網(wǎng)絡(luò)，模化了高雷諾數(shù)湍流邊界層中的渦黏項(xiàng)，并將其實(shí)時(shí)地引入流場(chǎng)求解過(guò)程，提高了計(jì)算準(zhǔn)確度，且在翼型或流動(dòng)條件變化時(shí)具有一定通用性。在湍流計(jì)算領(lǐng)域，基于大數(shù)據(jù)的機(jī)器學(xué)習(xí)方法具有廣闊的應(yīng)用前景，文獻(xiàn)[15]對(duì)數(shù)據(jù)驅(qū)動(dòng)式湍流建模理論和方法進(jìn)行了綜述分析。

數(shù)據(jù)驅(qū)動(dòng)式建模方法依賴于樣本數(shù)據(jù)的個(gè)數(shù)和質(zhì)量，小樣本或樣本冗余誤差會(huì)導(dǎo)致模型的失真或失效。為克服數(shù)據(jù)驅(qū)動(dòng)式建模對(duì)數(shù)據(jù)樣本個(gè)數(shù)和質(zhì)量的依賴關(guān)系，需要引入基于物理知識(shí)約束的先驗(yàn)條件，在建模過(guò)程中約束模型的邊界和搜索空間，有利于避免非物理現(xiàn)象的產(chǎn)生，提高模型的預(yù)測(cè)精度。近年來(lái)，基于物理知識(shí)的數(shù)據(jù)驅(qū)動(dòng)式建模成為研究熱點(diǎn)，Raissi等[16]以自由振動(dòng)的圓柱為研究對(duì)象，將有限的速度散點(diǎn)數(shù)據(jù)作為訓(xùn)練樣本，利用神經(jīng)網(wǎng)絡(luò)重構(gòu)的過(guò)程中，引入流體控制方程作為先驗(yàn)的物理知識(shí)約束，對(duì)流動(dòng)參數(shù)進(jìn)行了預(yù)測(cè)，獲得準(zhǔn)確的升/阻力預(yù)測(cè)結(jié)果。Raissi等[17]構(gòu)造了基于物理信息的神經(jīng)網(wǎng)絡(luò)，用于求解非線性偏微分方程，利用少量的訓(xùn)練數(shù)據(jù)，成功捕捉了系統(tǒng)的非線性特征。針對(duì)湍流模擬的RANS方程，尚未發(fā)現(xiàn)基于物理知識(shí)約束的數(shù)據(jù)驅(qū)動(dòng)式建模方法。

準(zhǔn)確預(yù)測(cè)湍流摩擦阻力具有重要的工程價(jià)值，本文作者在利用DNS數(shù)據(jù)研究湍流摩擦阻力的分解中，發(fā)現(xiàn)湍流摩擦阻力的分解項(xiàng)在湍流邊界層的內(nèi)層具有統(tǒng)一的比擬關(guān)系[18]，與雷諾數(shù)大小無(wú)關(guān)，與可壓縮性無(wú)關(guān)[19]，并在槽道和平板湍流邊界層中得到了驗(yàn)證[18-20]，結(jié)果表明該比擬關(guān)系是一種物理的、普適的規(guī)律，是一種先驗(yàn)的物理知識(shí)。本文基于該物理知識(shí)開展了針對(duì)S-A湍流模型的數(shù)據(jù)驅(qū)動(dòng)式建模及修正，在建模過(guò)程中引入物理知識(shí)約束，對(duì)比了有/無(wú)物理知識(shí)約束的修正對(duì)槽道湍流摩擦阻力預(yù)測(cè)精度的影響。

1 基于物理知識(shí)約束的數(shù)據(jù)驅(qū)動(dòng)式建模方法

1.1 經(jīng)典S-A湍流模型

為封閉雷諾平均的Navier-Stokes方程，Spalart和Allmaras[21]構(gòu)造了一方程湍流模型方程，通過(guò)求解標(biāo)量運(yùn)輸方程來(lái)計(jì)算湍流渦黏系數(shù)并獲得雷諾應(yīng)力，具體形式為

(1)

(2)

(3)

(4)

1.2 源于湍流摩擦阻力分解的物理知識(shí)

固壁與流體相對(duì)運(yùn)動(dòng)而產(chǎn)生摩擦，摩擦阻力可以表征為固壁對(duì)流體做功的一種形式，對(duì)于槽道湍流而言，摩擦阻力做功，一部分能量通過(guò)黏性直接耗散轉(zhuǎn)換成熱能，一部分能量則通過(guò)湍動(dòng)能生成項(xiàng)維持湍流態(tài)。基于這一思想，Renard和Deck[23]推導(dǎo)了湍流摩擦阻力系數(shù)Cf的分解公式：

(5)

(6)

式中：a1=2.644；b1=1.895；c1=-0.805；a2=1.777；b2=1.118；c2=1.642。式(6)中僅有一個(gè)自變量y+，圖1表明擬合結(jié)果與DNS數(shù)據(jù)的符合良好。本文將式(6)作為一種先驗(yàn)物理知識(shí)約束引入到數(shù)據(jù)驅(qū)動(dòng)建模中。

圖1 不同雷諾數(shù)下直接黏性耗散項(xiàng)的法向分布Fig.1 Direct viscous dissipation curves along wall-normal direction under different Reynolds numbers

1.3 基于物理知識(shí)約束的湍流模型修正方法

保留S-A湍流模型方程的基本形式，對(duì)生成項(xiàng)P添加修正因子β，修正后的湍流模型方程為

(7)

修正因子β對(duì)流場(chǎng)的影響是全局性的，其取值隨著法向高度的變化而變化。本文通過(guò)添加修正因子β來(lái)校正S-A湍流模型，其本質(zhì)是在S-A湍流模型中添加了源項(xiàng)的校正量δ，即δ=(β-1)P，其大小與生成項(xiàng)相關(guān)。初始時(shí)β=1，校正量δ均為0，式(7)即退化為原始的S-A方程。通過(guò)引入修正因子β，擬使得S-A湍流模型預(yù)測(cè)更準(zhǔn)確，即校正后的S-A湍流模型計(jì)算結(jié)果更接近于真實(shí)解，真實(shí)解可取高精度的DNS、LES結(jié)果，也可利用試驗(yàn)測(cè)量的數(shù)據(jù)進(jìn)行標(biāo)定，本文中的真實(shí)解為DNS計(jì)算結(jié)果[24]。

為獲得準(zhǔn)確可靠的修正因子，需要設(shè)定合理的目標(biāo)函數(shù)，結(jié)合高效的優(yōu)化方法，快速獲得修正因子分布。傳統(tǒng)方法對(duì)模型中系數(shù)進(jìn)行調(diào)整優(yōu)化，例如對(duì)馮卡門常數(shù)(κ)或湍流普朗特?cái)?shù)(σ)的校正[24-25]，但其難以克服重要流場(chǎng)物理特性(如雷諾數(shù))的影響。由此，考慮到湍流摩擦阻力的預(yù)測(cè)與平均流向速度剖面緊密相關(guān)，且引入式(6)作為物理知識(shí)約束，設(shè)定目標(biāo)函數(shù)為

(8)

式中：(Uj-Uj,DNS)2表征對(duì)平均流向速度剖面的逼近程度；(Ωj-Ω(y+))2表征對(duì)近壁黏性耗散項(xiàng)分布的逼近程度；Nc為網(wǎng)格單元數(shù)。值得指出的是，修正因子β與直接黏性耗散和平均流向速度之間不存在直接的量化關(guān)系，僅存在隱式關(guān)聯(lián)與映射。

為對(duì)比分析有/無(wú)物理知識(shí)約束對(duì)湍流摩擦阻力的預(yù)測(cè)誤差，本文也采用了基于平均流向速度的目標(biāo)函數(shù)：

(9)

由于設(shè)計(jì)變量個(gè)數(shù)與網(wǎng)格數(shù)量相當(dāng)，使用有限差分法求解梯度計(jì)算量巨大。伴隨方法(Adjoint Method)是一種高效的梯度求解方法，其計(jì)算規(guī)模與設(shè)計(jì)變量的數(shù)目基本無(wú)關(guān)，可大幅減小梯度的求解時(shí)間。伴隨方法以偏微分方程系統(tǒng)控制理論為基礎(chǔ)，Jameson[26]首次將伴隨方法應(yīng)用于氣動(dòng)設(shè)計(jì)中，把氣動(dòng)設(shè)計(jì)轉(zhuǎn)換為一個(gè)滿足特定約束的最優(yōu)控制問(wèn)題，梯度求解的計(jì)算量約為2倍的流場(chǎng)計(jì)算量。伴隨方法包括連續(xù)伴隨和離散伴隨，本文采用離散伴隨方法[27]，伴隨方程和梯度的表達(dá)式為

(10)

(11)

圖2 湍流模型修正優(yōu)化流程Fig.2 Optimization process of turbulence model correction

以二維方程[29]來(lái)驗(yàn)證伴隨方法所求得的梯度與有限差分法的一致性，圖3對(duì)比了有限差分法與伴隨方法所得的梯度分布，其偏差小于0.57%，屬于合理范圍內(nèi)。

圖3 伴隨方法中的梯度驗(yàn)證Fig.3 Gradient verification in adjoint method

2 計(jì)算結(jié)果與分析

圖4 直接黏性耗散項(xiàng)的法向分布Fig.4 Distribution of direct viscous dissipation in wall-normal direction

圖5 直接黏性耗散項(xiàng)相對(duì)誤差的法向分布Fig.5 Distribution of relative error of direct viscous dissipation in wall-normal direction

圖4和圖5中的比較結(jié)果表明，原始S-A湍流模型對(duì)近壁區(qū)直接黏性耗散項(xiàng)的預(yù)測(cè)存在明顯偏差；如果通過(guò)數(shù)據(jù)驅(qū)動(dòng)模型僅修正流向速度剖面，對(duì)直接黏性耗散項(xiàng)的預(yù)測(cè)有所改善，但誤差仍然較大；在數(shù)據(jù)驅(qū)動(dòng)建模的目標(biāo)函數(shù)中，進(jìn)一步引入關(guān)于直接黏性耗散的先驗(yàn)物理知識(shí)作為約束條件，預(yù)測(cè)精度顯著改善，計(jì)算結(jié)果與真實(shí)結(jié)果的相對(duì)誤差小于3%。

圖6 平均流向速度的法向分布Fig.6 Distribution of averaged streamwise velocity in wall-normal direction

圖7 平均流向速度相對(duì)誤差的法向分布Fig.7 Distribution of relative error of averaged streamwise velocity in wall-normal direction

圖8比較了不同雷諾數(shù)下修正因子β的變化情況，注意原始S-A湍流模型中修正因子恒為1。當(dāng)采用Ju為目標(biāo)函數(shù)時(shí)，僅根據(jù)流向平均速度進(jìn)行模型修正，修正因子β的法向變化趨勢(shì)基本與圖7中流向平均速度誤差的分布一致，說(shuō)明修正因子β具有較好的活性，可根據(jù)當(dāng)?shù)氐钠骄俣日`差對(duì)S-A模型中的湍流黏性生成項(xiàng)進(jìn)行相應(yīng)調(diào)整。當(dāng)采用Jud為目標(biāo)函數(shù)時(shí)，引入了直接黏性耗散項(xiàng)作為流向平均速度修正的補(bǔ)充，以Reτ=180的槽道為例，對(duì)于復(fù)合目標(biāo)函數(shù)修正因子β做出相應(yīng)的調(diào)整，在y+=4.65處達(dá)到峰值，說(shuō)明目標(biāo)函數(shù)Jud起到了復(fù)合校正效果。同時(shí)，圖8中β系數(shù)隨流場(chǎng)位置及當(dāng)?shù)亓鲌?chǎng)變量的變化而發(fā)生演變，也說(shuō)明了本文建立的伴隨優(yōu)化方法計(jì)算可以有效地建立對(duì)方程預(yù)測(cè)精度的調(diào)整和修正。需要注意的是，修正因子β的取值范圍在-10～20之間，說(shuō)明原始S-A湍流模型中的生成項(xiàng)存在低估或者逆?zhèn)鬟f現(xiàn)象。修正因子的本質(zhì)作用是在S-A湍流模型中引入δ=(β-1)P作為源項(xiàng)的修正項(xiàng)，該修正項(xiàng)的取值范圍根據(jù)預(yù)測(cè)結(jié)果與真實(shí)結(jié)果的偏差而定。

圖8 修正因子β的法向分布Fig.8 Distribution of correction factor β in wall-normal direction

湍流摩擦阻力是一個(gè)較小的壁面變量，目前針對(duì)大型客機(jī)的總阻力預(yù)測(cè)，利用RANS計(jì)算要確保阻力誤差在1 count以內(nèi)仍是巨大的挑戰(zhàn)，其中最為困難的是湍流摩擦阻力的準(zhǔn)確預(yù)測(cè)，其約占總阻力的50%以上。本文針對(duì)槽道湍流這一簡(jiǎn)單基礎(chǔ)構(gòu)型進(jìn)行了數(shù)據(jù)驅(qū)動(dòng)的湍流模型修正，在低雷諾數(shù)下其預(yù)測(cè)精度顯著提高，而在較高雷諾數(shù)條件下，表1中的數(shù)據(jù)顯示無(wú)修正的湍流模型對(duì)摩擦阻力的預(yù)測(cè)誤差已經(jīng)小于0.48%，通過(guò)模型修正，這個(gè)誤差得到了進(jìn)一步的減小。該結(jié)果表明，基于數(shù)據(jù)驅(qū)動(dòng)式的湍流模型修正具有較好的應(yīng)用前景。

表1 壁面摩擦阻力系數(shù)的相對(duì)誤差Table 1 Relative errors of skin friction coefficients

3 結(jié)論與展望

數(shù)據(jù)驅(qū)動(dòng)式建模的本質(zhì)是數(shù)據(jù)回歸和優(yōu)化問(wèn)題引入物理知識(shí)的先驗(yàn)約束，一定程度上可降低對(duì)數(shù)據(jù)樣本個(gè)數(shù)和質(zhì)量的依賴，提高模型的預(yù)測(cè)精度。本文提出了一種基于物理知識(shí)的數(shù)據(jù)驅(qū)動(dòng)式湍流模型修正方法，并以槽道湍流為研究對(duì)象驗(yàn)證了該方法的優(yōu)越性，獲得的結(jié)論包括：

1) 基于湍流摩擦阻力分解的先驗(yàn)物理知識(shí)，建立了一種基于物理知識(shí)約束的湍流模型修正方法，將物理知識(shí)顯性地引入目標(biāo)函數(shù)的設(shè)定中，結(jié)合伴隨優(yōu)化方法，高效率地獲得修正因子的分布，以達(dá)到提高預(yù)測(cè)精度的目的。

2) 以槽道湍流為例，對(duì)比了有/無(wú)物理知識(shí)約束下的湍流模型預(yù)測(cè)精度，結(jié)果表明引入物理知識(shí)約束后，可提高對(duì)湍流摩擦阻力的預(yù)測(cè)精度。

本文研究存在如下不足，希望在后期研究中克服：

1) 引入的基于湍流摩擦阻力分解的先驗(yàn)物理知識(shí)，適用于無(wú)逆壓梯度的湍流邊界層，而對(duì)于有逆壓梯度導(dǎo)致的流動(dòng)分離問(wèn)題，該物理知識(shí)的有效性需要進(jìn)一步的驗(yàn)證，或引入更為普適的物理知識(shí)。

2) 僅開展了基于DNS數(shù)據(jù)的先驗(yàn)驗(yàn)證，而沒(méi)有利用機(jī)器學(xué)習(xí)方法，如神經(jīng)網(wǎng)絡(luò)、決策樹等，建立物理變量與目標(biāo)函數(shù)之間的映射關(guān)系，以開展無(wú)數(shù)據(jù)條件下的后驗(yàn)驗(yàn)證。

3) 建立目標(biāo)函數(shù)時(shí)，用到了DNS平均速度剖面數(shù)據(jù)，沒(méi)有完全避免對(duì)DNS數(shù)據(jù)的依賴，在后期工作中有待進(jìn)一步完善。