閆重陽(yáng)，張宇飛，陳海昕

清華大學(xué) 航天航空學(xué)院，北京 100086

精確模擬湍流流動(dòng)是許多工程應(yīng)用和學(xué)術(shù)研究問(wèn)題的關(guān)鍵。由于較低的計(jì)算成本和較高的穩(wěn)健性，在工業(yè)界的CFD仿真領(lǐng)域，雷諾平均(Reynolds Averaged Navier-Stokes, RANS)模型將在可見的未來(lái)繼續(xù)占有重要地位[1]。然而，傳統(tǒng)的RANS模型是針對(duì)簡(jiǎn)單、經(jīng)典的流動(dòng)而校準(zhǔn)的，在許多實(shí)際的場(chǎng)景下缺少足夠的仿真能力，例如翼型分離流、強(qiáng)逆壓梯度流，或者具有高流線曲率的流動(dòng)[1]。為了彌補(bǔ)傳統(tǒng)RANS模型在某些流動(dòng)場(chǎng)景下的不足，人們?cè)岢鲞^(guò)多種改進(jìn)。例如，非線性渦黏模型放棄了雷諾應(yīng)力張量與應(yīng)變率張量成線性關(guān)系的假設(shè)；雷諾應(yīng)力模型則放棄渦黏假設(shè)，直接對(duì)雷諾應(yīng)力建立輸運(yùn)方程。然而，這些更為復(fù)雜的RANS模型在計(jì)算準(zhǔn)確度提高的同時(shí)，產(chǎn)生的代價(jià)是計(jì)算穩(wěn)定性下降、計(jì)算成本升高。另外，如何校準(zhǔn)這些復(fù)雜RANS模型所帶來(lái)的諸多新的參數(shù)也是對(duì)湍流建模者的一個(gè)挑戰(zhàn)。

作為一種有別于傳統(tǒng)湍流建模思路的方法，數(shù)據(jù)驅(qū)動(dòng)的湍流建模近年來(lái)正在逐漸興起。這一變化得益于數(shù)十年來(lái)計(jì)算資源的大幅增長(zhǎng)，用于湍流模型校準(zhǔn)的實(shí)驗(yàn)數(shù)據(jù)與DNS(Direct Numerical Simulation)數(shù)據(jù)的增加，以及采樣、統(tǒng)計(jì)推斷、機(jī)器學(xué)習(xí)方法在湍流建模中的應(yīng)用。從數(shù)據(jù)驅(qū)動(dòng)的湍流建模角度來(lái)看，湍流模型方程的模化過(guò)程，包括雷諾應(yīng)力的計(jì)算、模型方程源項(xiàng)的建立、方程系數(shù)的選取等等，都會(huì)引入不確定性，這種不確定性通過(guò)RANS方程的前向傳播，最終會(huì)產(chǎn)生CFD預(yù)測(cè)結(jié)果的不確定性。數(shù)據(jù)驅(qū)動(dòng)則意味著使用統(tǒng)計(jì)方法來(lái)量化或減少湍流模型的不確定度[2]。圍繞這個(gè)問(wèn)題，根據(jù)不同的應(yīng)用目的，可以將主要的研究方向分為3類。

第1類研究是RANS模型不確定度量化(UQ，又稱前向方法，無(wú)數(shù)據(jù)方法)。這類研究旨在確定模型不確定性對(duì)QoI(Quantity of Interest)分布的影響。當(dāng)給定湍流模型中不確定性參數(shù)的變化范圍或先驗(yàn)概率分布時(shí)，通過(guò)不確定度量化，可以得到QoI的變化范圍或概率分布。此類研究通常使用蒙特卡羅方法。例如，Platteeuw等針對(duì)零壓力梯度平板流動(dòng)，使用概率分配法(PCM)研究k-ε模型系數(shù)的不確定性影響[3]；Dunn等針對(duì)后臺(tái)階流動(dòng)，使用拉丁超立方采樣方法研究了k-ε模型系數(shù)對(duì)再附點(diǎn)、摩擦系數(shù)等的影響[4]。

除了模型系數(shù)以外，不確定性還可能來(lái)源于模型方程形式(例如方程源項(xiàng))或者?；?xiàng)本身(例如渦黏性ν，雷諾應(yīng)力τ)。Emory等針對(duì)平面槽道流、方腔二次流和跨聲速凸塊流動(dòng)，進(jìn)行了基于雷諾應(yīng)力不確定性的UQ分析：通過(guò)將雷諾應(yīng)力推向滿足可實(shí)現(xiàn)性約束的極限情況，來(lái)向模型注入雷諾應(yīng)力不確定性，從而得到QoI的分布范圍[5]。

第2類研究被稱為后向方法，即根據(jù)可觀測(cè)數(shù)據(jù)推斷模型不確定性參數(shù)的后驗(yàn)分布。假設(shè)模型不確定性參數(shù)為θ，θ的先驗(yàn)分布已知，此類研究的目的是在給定實(shí)驗(yàn)或DNS觀測(cè)數(shù)據(jù)G的情況下，計(jì)算出θ的后驗(yàn)分布p(θ|G)。這本質(zhì)上是一個(gè)貝葉斯推斷問(wèn)題：

p(θ|G)∝p(G│θ)p(θ)

(1)

第3類研究則引入機(jī)器學(xué)習(xí)方法以進(jìn)行模型改進(jìn)[13-14]。作為一種研究思路，Parish和Duraisamy提出了FIML(Field Inversion and Machine Learning)用于數(shù)據(jù)驅(qū)動(dòng)建模的范式[15]：首先通過(guò)流場(chǎng)反演得到空間分布的修正場(chǎng)，然后利用機(jī)器學(xué)習(xí)建立修正值與當(dāng)?shù)亓鲃?dòng)特征的映射，最后把該機(jī)器學(xué)習(xí)模型加入到RANS模型中，在每一步迭代中動(dòng)態(tài)給出修正量的值，這樣就得到了增強(qiáng)的RANS模型。Parish和Duraisamy[15]將FIML成功地應(yīng)用于一個(gè)一維傳熱問(wèn)題和k-ω模型計(jì)算的一維槽道流。Singh等也利用神經(jīng)網(wǎng)絡(luò)建立了當(dāng)?shù)亓鲃?dòng)特征到β(x)的映射，從而得到了能夠更準(zhǔn)確預(yù)測(cè)翼型分離流動(dòng)的增強(qiáng)版SA模型[16]。Duraisamy等則分別針對(duì)SA模型繞凸壁邊界層流動(dòng)、k-ω-γ三方程模型平板邊界層旁路轉(zhuǎn)捩問(wèn)題應(yīng)用了FIML[17]。

本文使用基于離散伴隨的流場(chǎng)反演，對(duì)2個(gè)典型的湍流分離問(wèn)題進(jìn)行了研究，即結(jié)冰翼型繞流和周期山分離流。利用有限的實(shí)驗(yàn)數(shù)據(jù)，對(duì)SA模型輸運(yùn)方程生成項(xiàng)的乘子進(jìn)行貝葉斯推斷，以達(dá)到對(duì)流場(chǎng)及湍流模型的修正。在本文中，針對(duì)周期山的流場(chǎng)反演可被歸結(jié)為一個(gè)約束優(yōu)化問(wèn)題，而常規(guī)的伴隨方法在解決約束優(yōu)化問(wèn)題上較為低效。本文發(fā)展了約束增廣的伴隨方法，從而使其處理周期山流場(chǎng)反演問(wèn)題的效率大大提高。計(jì)算結(jié)果表明，基于伴隨優(yōu)化的流場(chǎng)反演能以相當(dāng)高的精度擬合指定的實(shí)驗(yàn)數(shù)據(jù)。進(jìn)一步的分析展示了流場(chǎng)反演作為一種數(shù)據(jù)驅(qū)動(dòng)方法的意義：一方面通過(guò)對(duì)湍流模型不確定參數(shù)的貝葉斯推斷，把有限的實(shí)驗(yàn)數(shù)據(jù)中的信息泛化到整個(gè)流場(chǎng)，從而使修正后的整個(gè)流場(chǎng)在一定程度上更接近真實(shí)流場(chǎng)；另一方面給出了具有明確物理意義的修正量分布，可以直接作為下一步的機(jī)器學(xué)習(xí)增強(qiáng)模型的學(xué)習(xí)目標(biāo)，也可為人工湍流建模提供重要參考。

1 流場(chǎng)反演框架

1.1 流場(chǎng)反演問(wèn)題構(gòu)建

Spalart和Allmaras構(gòu)造的經(jīng)典SA一方程湍流模型中的標(biāo)量輸運(yùn)方程為[18]

(2)

(3)

P為生成項(xiàng)，D為破壞項(xiàng)，其具體形式為

(4)

(5)

(6)

式中：β為空間坐標(biāo)x的函數(shù)。從數(shù)據(jù)驅(qū)動(dòng)湍流建模的角度而言，β(x)是一個(gè)空間上的隨機(jī)場(chǎng)，遵循一定的概率分布。在獲得實(shí)驗(yàn)數(shù)據(jù)或DNS數(shù)據(jù)之前，可以為β人為指定一個(gè)先驗(yàn)分布p(β)，它代表了研究者對(duì)湍流建模的先驗(yàn)知識(shí)[19]。在得到實(shí)驗(yàn)數(shù)據(jù)或DNS數(shù)據(jù)d之后，可以根據(jù)貝葉斯定理推斷得到β的后驗(yàn)分布。其中d可以是積分量，例如升力系數(shù)；也可以是空間或物體表面某些位置的流場(chǎng)數(shù)據(jù)，例如表面摩擦系數(shù)Cf分布等。

在本文中，無(wú)需求出β(x)完整的后驗(yàn)概率分布函數(shù)，只需對(duì)β(x)進(jìn)行最大后驗(yàn)推斷即可。該優(yōu)化問(wèn)題可以用梯度優(yōu)化算法結(jié)合離散伴隨進(jìn)行求解。把所有空間點(diǎn)上的β值組裝為一個(gè)向量β，β是一個(gè)隨機(jī)向量。同理，所有空間點(diǎn)上的觀測(cè)值組成了觀測(cè)值向量d。假設(shè)β的先驗(yàn)p(β)、觀測(cè)值d都服從高斯分布，則求解β的最大后驗(yàn)值變?yōu)槿缦碌膬?yōu)化問(wèn)題：

(β-βprior)TCβ-1(β-βprior)]

(7)

式中：h(β)為模型在特定β取值下所預(yù)測(cè)的QoI；βprior為β先驗(yàn)分布的均值；Cm、Cβ為觀測(cè)值d以及修正量β先驗(yàn)分布的協(xié)方差矩陣。式(7)右邊實(shí)際上為p(β|d)的對(duì)數(shù)，其中第1項(xiàng)代表了模型預(yù)測(cè)的QoI與實(shí)驗(yàn)值之間的偏差，第2項(xiàng)代表了β與其先驗(yàn)均值之間的偏差。

(8)

(9)

式中：1為所有分量都為1的向量。此時(shí)的目標(biāo)函數(shù)可以寫為

(10)

σobs和σβ的相對(duì)大小決定了觀測(cè)值和先驗(yàn)對(duì)β后驗(yàn)影響的權(quán)重。至此，對(duì)于SA模型的流場(chǎng)反演問(wèn)題歸結(jié)為以β為設(shè)計(jì)變量，以式(10)為目標(biāo)函數(shù)的優(yōu)化問(wèn)題。

1.2 離散伴隨方法

伴隨方法是一種可以高效求解目標(biāo)函數(shù)對(duì)設(shè)計(jì)變量梯度的方法，常結(jié)合梯度類優(yōu)化算法用于解決優(yōu)化問(wèn)題。從Pironneau[20]、Jameson[21]開始，伴隨方法廣泛被用于飛機(jī)的氣動(dòng)優(yōu)化設(shè)計(jì)。伴隨方法的優(yōu)勢(shì)在于其計(jì)算量與設(shè)計(jì)變量的個(gè)數(shù)無(wú)關(guān)，因此特別適合于設(shè)計(jì)變量是全場(chǎng)分布的參數(shù)場(chǎng)的流場(chǎng)反演問(wèn)題。

早期的氣動(dòng)優(yōu)化設(shè)計(jì)所使用的伴隨方法被稱為連續(xù)伴隨。在該方法中，需要首先推導(dǎo)伴隨方程，并施加適當(dāng)?shù)倪吔鐥l件；然后將伴隨方程離散化并迭代求解[22]。然而，連續(xù)伴隨在方程推導(dǎo)(特別是黏性項(xiàng)、邊界條件的處理)、代碼實(shí)現(xiàn)上較為復(fù)雜，限制了它在工程上的應(yīng)用[23]。隨后，Elliott和Peraire[24], Nielsen和Anderson[25]發(fā)展了離散伴隨方法，首先將Navier-Stokes方程和目標(biāo)函數(shù)離散化，然后在此基礎(chǔ)上直接構(gòu)造離散形式的伴隨方程。相比連續(xù)伴隨，離散伴隨的計(jì)算速度更快，梯度求解的精度更高，并且大大減少了代碼開發(fā)的時(shí)間[22]。

氣動(dòng)目標(biāo)函數(shù)J、RANS方程的殘差R可以寫為空間網(wǎng)格坐標(biāo)x、空間流場(chǎng)變量w的函數(shù)：J(w,x)，R(w,x)。離散伴隨的求解目的是在所有網(wǎng)格點(diǎn)上RANS方程的殘差R(w,x)=0的約束下，得到目標(biāo)函數(shù)J(w,x)對(duì)設(shè)計(jì)變量x的全導(dǎo)數(shù)：

(11)

對(duì)R(w,x)取微分，得：

(12)

故

(13)

引入伴隨變量φ，滿足：

(14)

則

(15)

1.3 約束增廣伴隨方法

一些流場(chǎng)優(yōu)化問(wèn)題具有特定的氣動(dòng)約束，例如翼型/機(jī)翼定升力系數(shù)優(yōu)化，內(nèi)流道問(wèn)題的定質(zhì)量流量?jī)?yōu)化等等。對(duì)于每個(gè)氣動(dòng)約束，往往有一個(gè)增廣設(shè)計(jì)變量(如攻角、體積力等)，通過(guò)調(diào)整該增廣設(shè)計(jì)變量來(lái)滿足該氣動(dòng)約束。

作為一種解決方法，可以用梯度優(yōu)化算法(如SQP、SNOPT)處理約束優(yōu)化問(wèn)題。對(duì)于每一個(gè)氣動(dòng)約束，可以把它直接作為優(yōu)化算法的約束，同時(shí)增加一個(gè)對(duì)應(yīng)的增廣設(shè)計(jì)變量。伴隨算法則需要分別給出目標(biāo)函數(shù)和約束對(duì)原始及增廣設(shè)計(jì)變量的導(dǎo)數(shù)。然而，由于增廣設(shè)計(jì)變量的存在，設(shè)計(jì)空間過(guò)大，其中的大部分區(qū)域都是不滿足氣動(dòng)約束的，因此優(yōu)化效率較低。如果離散伴隨能直接提供目標(biāo)函數(shù)在滿足氣動(dòng)約束下對(duì)設(shè)計(jì)變量的全導(dǎo)數(shù)，那么就只需處理原始設(shè)計(jì)變量，可以大大提高約束優(yōu)化的效率。

(16)

式中：

式(16)展開為

(17)

步驟2逐行求出中間矩陣M。對(duì)每一個(gè)約束RC,i：

(18)

步驟3求增廣流場(chǎng)變量的伴隨φC:

(19)

步驟4求原始流場(chǎng)變量的伴隨φ：

(20)

(21)

相比普通伴隨方法，約束增廣伴隨方法只需針對(duì)每個(gè)氣動(dòng)約束再求解一次伴隨方程，所增加的計(jì)算量與普通伴隨相當(dāng)。由此可以大大提高氣動(dòng)約束下的伴隨優(yōu)化效率。

2 結(jié)果驗(yàn)證

2.1 結(jié)冰翼型繞流算例

2.1.1 算例介紹

結(jié)冰翼型的失速預(yù)測(cè)是飛機(jī)安全性十分關(guān)注的一個(gè)問(wèn)題。在大攻角下，結(jié)冰翼型在位于前緣的冰型上會(huì)發(fā)生流動(dòng)分離，并在后方形成較大的分離泡?，F(xiàn)有的SA模型對(duì)于這樣的分離流動(dòng)不能準(zhǔn)確預(yù)測(cè)。本文選取加蓋了角狀冰型944的GLC305翼型作為研究對(duì)象，對(duì)其在特定攻角下的流場(chǎng)進(jìn)行了流場(chǎng)反演。該冰型的實(shí)驗(yàn)數(shù)據(jù)由NASA格倫研究中心獲得，并被作為本文流場(chǎng)反演的觀測(cè)數(shù)據(jù)[28]。

GLC305_944結(jié)冰翼型的計(jì)算網(wǎng)格如圖1所示。流向和法向分別有457和97個(gè)網(wǎng)格點(diǎn)，同時(shí)在冰型附近的加速區(qū)也做了一定網(wǎng)格加密。為保證壁面y+<1，壁面第1層網(wǎng)格的高度小于1×106c(c為弦長(zhǎng))。計(jì)算的來(lái)流條件為：來(lái)流馬赫數(shù)Ma∞=0.12,雷諾數(shù)Re=3.5×106，攻角α=6°。

圖1 GLC305_944結(jié)冰翼型的計(jì)算網(wǎng)格Fig.1 Computation grid for iced airfoil GLC305_944

原始SA模型在該工況下計(jì)算出的升力系數(shù)CL=0.49, 而實(shí)驗(yàn)CL=0.67。圖2展示了原始SA模型的流場(chǎng)分布圖，其中速度u是在ADFlow求解器內(nèi)無(wú)量綱化后的值。從中可以看出預(yù)測(cè)結(jié)果的翼型上表面全部處于分離狀態(tài)，這可以解釋計(jì)算得到的升力系數(shù)偏低的原因。圖3展示了計(jì)算得到的升力系數(shù)Cp分布與實(shí)驗(yàn)結(jié)果的比較。由圖可知，原始SA模型計(jì)算得到的吸力平臺(tái)遠(yuǎn)低于實(shí)驗(yàn)觀測(cè)值。

圖2 原始SA模型的流場(chǎng)圖Fig.2 Flow field of original SA model

圖3 計(jì)算得到的Cp分布與實(shí)驗(yàn)結(jié)果的比較Fig.3 Comparation of Cp distribution by computation and experiment

2.1.2 結(jié)果分析

在最簡(jiǎn)單的情形下，結(jié)冰翼型流場(chǎng)反演所用到的實(shí)驗(yàn)數(shù)據(jù)僅為CL，修正量β的先驗(yàn)分布是均值為1、方差全場(chǎng)相等的高斯分布。此時(shí)的目標(biāo)函數(shù)可寫為

(22)

圖4展示了流場(chǎng)反演的目標(biāo)函數(shù)、翼型升力系數(shù)隨迭代次數(shù)的變化。在僅僅14步以內(nèi)，SA模型經(jīng)過(guò)流場(chǎng)反演就能以～10-8的精度擬合實(shí)驗(yàn)的CL，這說(shuō)明SA模型在本文的修正框架下具有準(zhǔn)確預(yù)測(cè)結(jié)冰翼型分離流動(dòng)的潛力，同時(shí)顯示了所使用的離散伴隨方法的高效性。

圖4 流場(chǎng)反演的目標(biāo)函數(shù)、翼型升力系數(shù)隨 迭代次數(shù)的變化Fig.4 Variation of object function and lift coefficient of airfoil with iterations during field inversion

圖5(a)展示了流場(chǎng)反演得到的最終流場(chǎng)。由圖可知，流場(chǎng)反演之后所計(jì)算出的分離泡長(zhǎng)度由整個(gè)翼型上表面縮短為約43%弦長(zhǎng)。圖5(b)展示了流場(chǎng)反演所得到的修正量β的分布。從圖中可以看出，修正區(qū)域主要位于翼型上表面的分離剪切層中，且分為2個(gè)較為明顯的部分。第1部分是在分離泡前部、剪切層起始的位置，該區(qū)域內(nèi)的生成項(xiàng)乘子被減小；第2部分位于靠后的位置，此區(qū)域內(nèi)的生成項(xiàng)乘子被增大。此修正結(jié)果表明，在結(jié)冰翼型分離流動(dòng)預(yù)測(cè)中，誤差最主要來(lái)源于分離剪切層。在剪切層的起始部分，SA模型傾向于高估渦黏性；而在剪切層的后半部分，SA模型則傾向于低估渦黏性。這與李浩然等[29]的結(jié)論一致：減小剪切層前部的渦黏性有助于模擬分離泡內(nèi)從層流向湍流轉(zhuǎn)捩的過(guò)程，而增大剪切層后部的渦黏性有助于模擬該處湍流的非均衡特性，增加剪切層內(nèi)的湍流摻混，從而促進(jìn)分離泡再附。這表明，流場(chǎng)反演能夠以數(shù)據(jù)驅(qū)動(dòng)的方式，僅利用有限的觀測(cè)信息，正確推斷出整個(gè)流場(chǎng)范圍內(nèi)對(duì)湍流模型不確定量所應(yīng)做出的修正。

圖5 以CL作為目標(biāo)進(jìn)行流場(chǎng)反演得到的 最終流場(chǎng)及修正量分布Fig.5 Final flow field and distribution of corrections of field inversion with CL as object

圖6展示了流場(chǎng)反演得到的翼型表面壓力分布與實(shí)驗(yàn)結(jié)果的比較。由圖可知，雖然流場(chǎng)反演得到的升力系數(shù)與實(shí)驗(yàn)結(jié)果高度一致，但壓力分布形態(tài)仍存在一定差別。從Cp的分布形態(tài)上可以推測(cè)得出，流場(chǎng)反演得到的分離泡長(zhǎng)度更短，且分離泡再附位置的形狀更陡峭。出現(xiàn)這種差異的原因是，流場(chǎng)反演本質(zhì)上是基于數(shù)據(jù)驅(qū)動(dòng)，把當(dāng)前的先驗(yàn)和觀測(cè)數(shù)據(jù)信息向整個(gè)流場(chǎng)進(jìn)行泛化，所得到的流場(chǎng)僅是基于當(dāng)前信息下概率最大的流場(chǎng)。相比于未修正的計(jì)算結(jié)果，這個(gè)流場(chǎng)與真實(shí)流場(chǎng)更為接近，但仍存在一定的差別。

圖6 以CL作為目標(biāo)進(jìn)行流場(chǎng)反演得到的翼型表面 壓力分布與實(shí)驗(yàn)結(jié)果的比較Fig.6 Comparation of surface pressure distributions of airfoil between field inversion with CL as object and experiment results

(23)

式中：λ為一個(gè)超參數(shù)，由經(jīng)驗(yàn)給定。用于反演的Cp數(shù)據(jù)點(diǎn)全部位于壓力分布光順的部分，而排除了翼型前緣不光滑的結(jié)冰區(qū)域。圖7、圖8展示了該情況下流場(chǎng)反演的結(jié)果。從圖中可以看出，除了翼型上表面前緣的小部分區(qū)域外，反演得到的流場(chǎng)在翼型大部分位置的Cp都與實(shí)驗(yàn)符合較好。與之對(duì)應(yīng)的是分離泡的長(zhǎng)度也有所加長(zhǎng)，且再附位置的分離泡形狀更加平緩。這表明，在為流場(chǎng)反演增加了新的可觀測(cè)數(shù)據(jù)信息之后，流場(chǎng)反演的結(jié)果能夠更加接近真實(shí)流場(chǎng)。

圖7 以CL、Cp作為目標(biāo)進(jìn)行流場(chǎng)反演得到的 壓力分布與實(shí)驗(yàn)結(jié)果的比較Fig.7 Comparation of pressure distributions between field inversion with CL and Cp as objects and experiment results

圖8 以CL、Cp作為目標(biāo)進(jìn)行流場(chǎng)反演得到的 最終流場(chǎng)及修正量分布Fig.8 Final flow field and distribution of corrections of field inversion with CL and Cp as objects

2.2 周期山分離算例

2.2.1 算例介紹

周期山算例是另一個(gè)具有大分離流動(dòng)特征的算例，由于其有大量實(shí)驗(yàn)和DNS數(shù)據(jù)，故被廣泛用于測(cè)試CFD模型[6]。本文所計(jì)算的周期山模型幾何及網(wǎng)格如圖9所示，流向和法向的網(wǎng)格數(shù)分別為81和77。在本文所研究的工況下，基于山頂高度H和山頂截面上的平均速度Ub的雷諾數(shù)Reb=5 600。

圖9 周期山模型幾何及網(wǎng)格Fig.9 Geometry model and computation grid of periodic hill

周期山的上下壁面使用等溫固壁邊界條件，進(jìn)出口使用周期邊界條件。為了保證雷諾數(shù)一定，需要保證山頂處的質(zhì)量流量一定。為此，給全場(chǎng)施加一個(gè)均勻的體積力場(chǎng)，體積力沿x軸正方向。在CFD計(jì)算過(guò)程中，由程序根據(jù)當(dāng)前的質(zhì)量流量自動(dòng)調(diào)整體積力的大小，使其達(dá)到目標(biāo)流量。

圖10 原始SA模型的計(jì)算結(jié)果與DNS結(jié)果的對(duì)比Fig.10 Comparation between computation results of original SA model and DNS

(24)

本文的周期山算例具有質(zhì)量流量約束，因此需要使用1.3節(jié)中的約束增廣伴隨方法。為了驗(yàn)證約束增廣伴隨方法的高效性，本文同時(shí)測(cè)試了使用和不使用該方法進(jìn)行流場(chǎng)反演的2種方法，并進(jìn)行結(jié)果對(duì)比。方法1使用該方法，由CFD程序進(jìn)行定質(zhì)量流量計(jì)算，由離散伴隨給出定質(zhì)量流量下目標(biāo)函數(shù)對(duì)β的梯度，由優(yōu)化算法更新β。方法2使用普通伴隨，把體積力作為一個(gè)設(shè)計(jì)變量，CFD程序在給定體積力下進(jìn)行計(jì)算，由離散伴隨給出目標(biāo)函數(shù)對(duì)β和體積力的梯度，把質(zhì)量流量作為優(yōu)化問(wèn)題的約束，由優(yōu)化算法更新β和體積力。

2.2.2 結(jié)果分析

圖11展示了使用2種方法流場(chǎng)反演的收斂結(jié)果。由圖可知，使用約束增廣伴隨方法能在1/4的迭代步數(shù)內(nèi)達(dá)到與使用普通伴隨同樣的收斂效果。這說(shuō)明對(duì)于帶物理約束的優(yōu)化問(wèn)題，使用約束增廣伴隨方法的確能大大提高優(yōu)化效率。從質(zhì)量流量的變化曲線來(lái)看，約束增廣伴隨方法的質(zhì)量流量始終等于目標(biāo)流量，而普通伴隨的質(zhì)量流量在很長(zhǎng)一段時(shí)間內(nèi)在目標(biāo)流量上下波動(dòng)，這說(shuō)明普通伴隨的很多迭代都位于不滿足約束的設(shè)計(jì)空間內(nèi)，因此大大降低了優(yōu)化效率。

圖11 使用/不使用約束增廣伴隨的周期 山流場(chǎng)反演收斂情況Fig.11 Convergence history of flow field inversion of periodic hill with/without constraint augmented adjoint method

圖12(a)比較了原始SA模型、DNS，以及2種流場(chǎng)反演方法得到的周期山表面Cf分布。由圖可知，原始SA模型與DNS的計(jì)算結(jié)果存在較大差別，而無(wú)論是否使用約束增廣伴隨方法，流場(chǎng)反演都能準(zhǔn)確擬合指定的觀測(cè)數(shù)據(jù)。圖12(b)展示了原始SA模型、DNS，以及2種流場(chǎng)反演方法得到的空間不同截面的速度分布，其中x和y都以山頂高度H進(jìn)行無(wú)量綱化；馬赫數(shù)的橫向尺度被放大6倍以便于觀察。由圖可知，流場(chǎng)反演得到的速度分布在流場(chǎng)底部與DNS符合較好，但在流場(chǎng)上部則與DNS有一定的差別，且2種伴隨方法的結(jié)果之間也有一定誤差，但相比SA模型有較大的改進(jìn)。

圖12 原始SA模型、DNS以及2種流場(chǎng)反演方法得到的周期山表面Cf分布和空間不同截面速度型分布Fig.12 Cf distributions on surface of periodic hill and velocity profile distributions at different sections computed with original SA model, DNS and two methods of field inversion

圖13展示了2種方法得到的流場(chǎng)與修正量云圖。由圖可知，與原始SA模型相比，流場(chǎng)反演得到的分離泡再附位置均有所提前，這與DNS的結(jié)果一致。然而，流場(chǎng)反演得到的分離泡形狀與DNS結(jié)果相比有一定的差別。另一方面，雖然2種方法的修正區(qū)域大致相同，但具體修正量也存在一定的區(qū)別。

圖13 2種方法得到的流場(chǎng)與修正量云圖Fig.13 Contours of flow field and corrections computed with two methods

綜合以上結(jié)果可知，流場(chǎng)反演具有多解性，即在準(zhǔn)確擬合指定的觀測(cè)數(shù)據(jù)的基礎(chǔ)上，可能存在不同的全場(chǎng)修正量分布。本文中的多解性是由伴隨優(yōu)化的局部最優(yōu)特性造成的。流場(chǎng)反演的多解性同時(shí)表明，流場(chǎng)反演本質(zhì)上是數(shù)據(jù)驅(qū)動(dòng)的對(duì)先驗(yàn)和觀測(cè)信息的泛化，其目的是得到當(dāng)前信息下概率最大的一種可能的流場(chǎng)，它可能與真實(shí)流場(chǎng)仍存在一定差別。由于本文所用到的實(shí)驗(yàn)數(shù)據(jù)僅為周期山下壁面的Cf分布，因此反演得到的下半部分流場(chǎng)相比上半部分更能擬合真實(shí)流場(chǎng)。如何讓流場(chǎng)的反演得到的結(jié)果從整體上更接近真實(shí)流場(chǎng)，而不僅僅是擬合指定的觀測(cè)數(shù)據(jù)，仍是一個(gè)有待研究的問(wèn)題。

3 結(jié) 論

傳統(tǒng)的SA模型在計(jì)算結(jié)冰翼型分離流、周期山流動(dòng)等具有大分離的流動(dòng)時(shí)往往不能較為準(zhǔn)確地進(jìn)行預(yù)測(cè)。本文基于數(shù)據(jù)驅(qū)動(dòng)的思路，假設(shè)SA模型的誤差來(lái)源于其渦黏性輸運(yùn)方程的生成項(xiàng)所乘的系數(shù)，使用基于離散伴隨的流場(chǎng)反演框架，根據(jù)有限的實(shí)驗(yàn)數(shù)據(jù)對(duì)該系數(shù)在全場(chǎng)的分布進(jìn)行修正。對(duì)結(jié)冰翼型分離流、周期山流動(dòng)2個(gè)流動(dòng)問(wèn)題進(jìn)行了研究，得到了以下結(jié)論：

1) 無(wú)論指定擬合的觀測(cè)數(shù)據(jù)是升力系數(shù)，壓力分布，還是壁面摩阻系數(shù)分布，流場(chǎng)反演均能以相當(dāng)高的精度修正模型使其預(yù)測(cè)值與指定的實(shí)驗(yàn)值吻合。

2) 流場(chǎng)反演能夠?qū)⒂邢薜挠^測(cè)數(shù)據(jù)信息泛化到整個(gè)流場(chǎng)，從而得到對(duì)整個(gè)流場(chǎng)的修正。更多的觀測(cè)數(shù)據(jù)有助于流場(chǎng)反演得到更接近真實(shí)值的流場(chǎng)。

3) 作為一種數(shù)據(jù)驅(qū)動(dòng)的方法，流場(chǎng)反演能夠僅使用極其有限的實(shí)驗(yàn)數(shù)據(jù)，正確推斷出對(duì)SA模型所要進(jìn)行修正的區(qū)域及修正量。這為進(jìn)一步的湍流模型改進(jìn)提供了參考。

4) 周期山算例表明，本文采用的約束增廣伴隨優(yōu)化的效率是普通伴隨的4倍以上。這為今后其他帶有物理約束的流場(chǎng)反演或氣動(dòng)優(yōu)化問(wèn)題提供了高效的解決思路。

5) 周期山算例展現(xiàn)了流場(chǎng)反演的多解性。這表明，流場(chǎng)反演本質(zhì)上得到的是概率最大的可能結(jié)果，這一結(jié)果總體上更接近真實(shí)流場(chǎng)，但其具體取值則取決于給定的先驗(yàn)及觀測(cè)數(shù)據(jù)，同時(shí)也受流場(chǎng)反演所用的優(yōu)化算法特性的影響。

本文研究結(jié)果表明，作為一種數(shù)據(jù)驅(qū)動(dòng)的湍流建模方法，基于離散伴隨的流場(chǎng)反演能夠從有限的實(shí)驗(yàn)數(shù)據(jù)中推斷得到整體流場(chǎng)以及修正量分布。這可以進(jìn)一步為基于機(jī)器學(xué)習(xí)的模型增強(qiáng)提供學(xué)習(xí)目標(biāo)，也能對(duì)傳統(tǒng)的湍流建模提供重要參考。