基于離散伴隨的流場反演在湍流模擬中的應用

閆重陽，張宇飛，陳海昕

清華大學 航天航空學院，北京 100086

精確模擬湍流流動是許多工程應用和學術研究問題的關鍵。由于較低的計算成本和較高的穩健性，在工業界的CFD仿真領域，雷諾平均(Reynolds Averaged Navier-Stokes, RANS)模型將在可見的未來繼續占有重要地位[1]。然而，傳統的RANS模型是針對簡單、經典的流動而校準的，在許多實際的場景下缺少足夠的仿真能力，例如翼型分離流、強逆壓梯度流，或者具有高流線曲率的流動[1]。為了彌補傳統RANS模型在某些流動場景下的不足，人們曾提出過多種改進。例如，非線性渦黏模型放棄了雷諾應力張量與應變率張量成線性關系的假設；雷諾應力模型則放棄渦黏假設，直接對雷諾應力建立輸運方程。然而，這些更為復雜的RANS模型在計算準確度提高的同時，產生的代價是計算穩定性下降、計算成本升高。另外，如何校準這些復雜RANS模型所帶來的諸多新的參數也是對湍流建模者的一個挑戰。

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

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

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

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

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

(1)

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

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

1 流場反演框架

1.1 流場反演問題構建

Spalart和Allmaras構造的經典SA一方程湍流模型中的標量輸運方程為[18]

(2)

(3)

P為生成項，D為破壞項，其具體形式為

(4)

(5)

(6)

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

在本文中，無需求出β(x)完整的后驗概率分布函數，只需對β(x)進行最大后驗推斷即可。該優化問題可以用梯度優化算法結合離散伴隨進行求解。把所有空間點上的β值組裝為一個向量β，β是一個隨機向量。同理，所有空間點上的觀測值組成了觀測值向量d。假設β的先驗p(β)、觀測值d都服從高斯分布，則求解β的最大后驗值變為如下的優化問題：

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

(7)

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

(8)

(9)

式中：1為所有分量都為1的向量。此時的目標函數可以寫為

(10)

σobs和σβ的相對大小決定了觀測值和先驗對β后驗影響的權重。至此，對于SA模型的流場反演問題歸結為以β為設計變量，以式(10)為目標函數的優化問題。

1.2 離散伴隨方法

伴隨方法是一種可以高效求解目標函數對設計變量梯度的方法，常結合梯度類優化算法用于解決優化問題。從Pironneau[20]、Jameson[21]開始，伴隨方法廣泛被用于飛機的氣動優化設計。伴隨方法的優勢在于其計算量與設計變量的個數無關，因此特別適合于設計變量是全場分布的參數場的流場反演問題。

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

氣動目標函數J、RANS方程的殘差R可以寫為空間網格坐標x、空間流場變量w的函數：J(w,x)，R(w,x)。離散伴隨的求解目的是在所有網格點上RANS方程的殘差R(w,x)=0的約束下，得到目標函數J(w,x)對設計變量x的全導數：

(11)

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

(12)

故

(13)

引入伴隨變量φ，滿足：

(14)

則

(15)

1.3 約束增廣伴隨方法

一些流場優化問題具有特定的氣動約束，例如翼型/機翼定升力系數優化，內流道問題的定質量流量優化等等。對于每個氣動約束，往往有一個增廣設計變量(如攻角、體積力等)，通過調整該增廣設計變量來滿足該氣動約束。

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

(16)

式中：

式(16)展開為

(17)

步驟2逐行求出中間矩陣M。對每一個約束RC,i：

(18)

步驟3求增廣流場變量的伴隨φC:

(19)

步驟4求原始流場變量的伴隨φ：

(20)

(21)

相比普通伴隨方法，約束增廣伴隨方法只需針對每個氣動約束再求解一次伴隨方程，所增加的計算量與普通伴隨相當。由此可以大大提高氣動約束下的伴隨優化效率。

2 結果驗證

2.1 結冰翼型繞流算例

2.1.1 算例介紹

結冰翼型的失速預測是飛機安全性十分關注的一個問題。在大攻角下，結冰翼型在位于前緣的冰型上會發生流動分離，并在后方形成較大的分離泡。現有的SA模型對于這樣的分離流動不能準確預測。本文選取加蓋了角狀冰型944的GLC305翼型作為研究對象，對其在特定攻角下的流場進行了流場反演。該冰型的實驗數據由NASA格倫研究中心獲得，并被作為本文流場反演的觀測數據[28]。

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

圖1 GLC305_944結冰翼型的計算網格Fig.1 Computation grid for iced airfoil GLC305_944

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

圖2 原始SA模型的流場圖Fig.2 Flow field of original SA model

圖3 計算得到的Cp分布與實驗結果的比較Fig.3 Comparation of Cp distribution by computation and experiment

2.1.2 結果分析

在最簡單的情形下，結冰翼型流場反演所用到的實驗數據僅為CL，修正量β的先驗分布是均值為1、方差全場相等的高斯分布。此時的目標函數可寫為

(22)

圖4展示了流場反演的目標函數、翼型升力系數隨迭代次數的變化。在僅僅14步以內，SA模型經過流場反演就能以～10-8的精度擬合實驗的CL，這說明SA模型在本文的修正框架下具有準確預測結冰翼型分離流動的潛力，同時顯示了所使用的離散伴隨方法的高效性。

圖4 流場反演的目標函數、翼型升力系數隨 迭代次數的變化Fig.4 Variation of object function and lift coefficient of airfoil with iterations during field inversion

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

圖5 以CL作為目標進行流場反演得到的 最終流場及修正量分布Fig.5 Final flow field and distribution of corrections of field inversion with CL as object

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

圖6 以CL作為目標進行流場反演得到的翼型表面 壓力分布與實驗結果的比較Fig.6 Comparation of surface pressure distributions of airfoil between field inversion with CL as object and experiment results

(23)

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

圖7 以CL、Cp作為目標進行流場反演得到的 壓力分布與實驗結果的比較Fig.7 Comparation of pressure distributions between field inversion with CL and Cp as objects and experiment results

圖8 以CL、Cp作為目標進行流場反演得到的 最終流場及修正量分布Fig.8 Final flow field and distribution of corrections of field inversion with CL and Cp as objects

2.2 周期山分離算例

2.2.1 算例介紹

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

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

周期山的上下壁面使用等溫固壁邊界條件，進出口使用周期邊界條件。為了保證雷諾數一定，需要保證山頂處的質量流量一定。為此，給全場施加一個均勻的體積力場，體積力沿x軸正方向。在CFD計算過程中，由程序根據當前的質量流量自動調整體積力的大小，使其達到目標流量。

圖10 原始SA模型的計算結果與DNS結果的對比Fig.10 Comparation between computation results of original SA model and DNS

(24)

本文的周期山算例具有質量流量約束，因此需要使用1.3節中的約束增廣伴隨方法。為了驗證約束增廣伴隨方法的高效性，本文同時測試了使用和不使用該方法進行流場反演的2種方法，并進行結果對比。方法1使用該方法，由CFD程序進行定質量流量計算，由離散伴隨給出定質量流量下目標函數對β的梯度，由優化算法更新β。方法2使用普通伴隨，把體積力作為一個設計變量，CFD程序在給定體積力下進行計算，由離散伴隨給出目標函數對β和體積力的梯度，把質量流量作為優化問題的約束，由優化算法更新β和體積力。

2.2.2 結果分析

圖11展示了使用2種方法流場反演的收斂結果。由圖可知，使用約束增廣伴隨方法能在1/4的迭代步數內達到與使用普通伴隨同樣的收斂效果。這說明對于帶物理約束的優化問題，使用約束增廣伴隨方法的確能大大提高優化效率。從質量流量的變化曲線來看，約束增廣伴隨方法的質量流量始終等于目標流量，而普通伴隨的質量流量在很長一段時間內在目標流量上下波動，這說明普通伴隨的很多迭代都位于不滿足約束的設計空間內，因此大大降低了優化效率。

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

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

圖12 原始SA模型、DNS以及2種流場反演方法得到的周期山表面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種方法得到的流場與修正量云圖。由圖可知，與原始SA模型相比，流場反演得到的分離泡再附位置均有所提前，這與DNS的結果一致。然而，流場反演得到的分離泡形狀與DNS結果相比有一定的差別。另一方面，雖然2種方法的修正區域大致相同，但具體修正量也存在一定的區別。

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

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

3 結 論

傳統的SA模型在計算結冰翼型分離流、周期山流動等具有大分離的流動時往往不能較為準確地進行預測。本文基于數據驅動的思路，假設SA模型的誤差來源于其渦黏性輸運方程的生成項所乘的系數，使用基于離散伴隨的流場反演框架，根據有限的實驗數據對該系數在全場的分布進行修正。對結冰翼型分離流、周期山流動2個流動問題進行了研究，得到了以下結論：

1) 無論指定擬合的觀測數據是升力系數，壓力分布，還是壁面摩阻系數分布，流場反演均能以相當高的精度修正模型使其預測值與指定的實驗值吻合。

2) 流場反演能夠將有限的觀測數據信息泛化到整個流場，從而得到對整個流場的修正。更多的觀測數據有助于流場反演得到更接近真實值的流場。

3) 作為一種數據驅動的方法，流場反演能夠僅使用極其有限的實驗數據，正確推斷出對SA模型所要進行修正的區域及修正量。這為進一步的湍流模型改進提供了參考。

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

5) 周期山算例展現了流場反演的多解性。這表明，流場反演本質上得到的是概率最大的可能結果，這一結果總體上更接近真實流場，但其具體取值則取決于給定的先驗及觀測數據，同時也受流場反演所用的優化算法特性的影響。

本文研究結果表明，作為一種數據驅動的湍流建模方法，基于離散伴隨的流場反演能夠從有限的實驗數據中推斷得到整體流場以及修正量分布。這可以進一步為基于機器學習的模型增強提供學習目標，也能對傳統的湍流建模提供重要參考。