湍流數據同化技術及應用

何創新，鄧志文，劉應征,*

1. 上海交通大學 機械與動力工程學院 動力機械與工程教育部重點實驗室，上海 200240 2.上海交通大學 燃氣輪機研究院，上海 200240

湍流普遍存在于航空航天、能源動力和海洋工程等領域，通常伴隨著傳熱、傳質、振動、噪聲等復雜的物理現象，因而相關研究對湍流流場的高精度實驗測量與數值模擬提出了很高的要求。高頻響的熱線風速儀(Hot Wire anemometer, HW)、高空間分辨率的粒子圖像測速技術(Particle Image Velocimetry, PIV)以及最新發展起來的層析PIV技術(Tomography Particle Image Velocimetry, Tomo-PIV)[1-3]為湍流流場分別提供了一維、二維和三維的測量手段。高精度大渦模擬(Large-Eddy Simulation, LES)和分離渦模擬(Detached-Eddy Simulation, DES)[4-6]方法為湍流多物理場的獲取提供了全三維高時空分辨率的解決方案。然而，單純實驗測量和數值模擬仍然無法滿足日益深入的研究需求。實驗測量方面：① 測量的空間維度有限，即使先進的三維PIV測量技術，其測量空間和流速范圍也受到硬件多方面的制約；② 測量區域受到幾何結構的遮擋與阻礙，復雜內流場及全場測量難以實現；③ 實驗測量誤差問題，尤其是近壁區的PIV測量誤差很大；④ 多物理場同步測量難度大。數值模擬方面：① 高精度LES和DES模擬技術是建立在高分辨率網格的基礎上，尤其對于含有邊界層轉捩的外流場模擬，計算量非常大；② 雷諾時均(Reynolds-Averaged Navier-Stokes, RANS)方法基于半理論半經驗的模型方程，對復雜流場模擬的準確度比較差；③ 數值計算需要給定準確的邊界條件，而這在很多工程流場中難以確定。

數據同化(Data Assimilation, DA)[7]是一種融合實驗測量與模型預測的數學方法，最早在氣象預報中得到發展與應用[8]，現階段還被廣泛應用于海洋、水文、地理等領域[9-11]。數據同化通過引入觀測站獲取的實時測量數據作為約束，植入流動與熱動力學方程中，從而預測未來時刻的天氣狀況。實驗觀測(Observation)、預測模型(Predictive model)和同化算法是數據同化的三要素。實驗觀測通過同化算法對預測模型進行修正，使之獲得更準確的預測結果。近年來，國外研究學者也開始將數據同化技術應用到湍流動力學研究領域[12-13]。數據同化方法可以很好地發揮實驗測量與數值模擬的各自優勢，在未來湍流研究中將具有很大的應用潛力。

在湍流數據同化研究中，實驗觀測可以是熱線風速儀、PIV等常規測試方法獲取的一維、二維乃至三維流場矢量分布，也可以是壓力傳感器、壓敏漆(Pressure Sensitive Paint, PSP)等方法獲取的壁面壓力分布；預測模型一般包括Navier-Stokes (N-S)方程、連續性方程、湍流封閉方程等；同化算法種類則比較多，常用的算法包括四維變分(4DVar)[14]、集合卡爾曼濾波(EnKF)[15]、伴隨(Adjoint formulation)[16]等。從廣義上講，最近興起的基于物理信息的機器學習(Physics informed machine learning)[17]方法也可歸納到同化算法當中，而不同的是其對歷史測量數據庫的依賴：通過對大樣本歷史數據庫的學習從而構建從實驗觀測到待測物理場的虛擬映射關系。從流場形式上來講，數據同化可分為穩態數據同化和非穩態數據同化[7]。顧名思義，穩態數據同化針對的流場一般為湍流雷諾時均場。非穩態數據同化與穩態相對，如4DVar方法。順序數據同化作為非穩態數據同化的一種，是指在預測模型的時間正向積分中，不斷引入動態觀測數據對模型進行約束，從而修正模型的運行軌跡。動態數據的約束可以在每個時間步進行，也可以在模型運行過程中的若干時刻進行，而后者對于觀測數據之間時刻的流場信息，依賴于模型本身的準確性以及前一觀測時刻同化結果帶來的時間相關性約束。在穩態數據同化中，預測模型一般包含RANS湍流模型方程[18-19]，通過數據同化算法對模型方程進行修正。而這種RANS模型修正方法在順序數據同化中效果不佳，因為湍流渦黏模型的修正不足以使流場得到快速響應，而一般在N-S方程中直接添加雷諾應力項以加快流場的改變。

鑒于數據同化的原理及優勢，單純的實驗測量或數值模擬方法的缺陷有望得到彌補，從而發展實驗測量與數值模擬的深度融合技術。本文結合國內外數據同化技術在湍流領域的研究成果，圍繞其基本數學原理、實現方法、發展現狀及應用進行綜述。

1 湍流測量與數值模擬的挑戰

湍流問題是困擾人類的世紀難題，其時間與空間上的多尺度特性決定了它與其他力學問題有著本質的區別。常規的湍流研究方法包含實驗測量和數值模擬兩大類，以獲取流場速度矢量和其他衍生物理信息的時空分布為主要目標。

1.1 實驗測量技術

湍流測量已由單點或者陣列式的傳感器測量發展到了多維度的激光無侵入式測量技術。雖然高頻響低噪聲的點式測量技術在現階段仍發揮著重要作用，但是高時空分辨率的多維和多場測量技術具有更重要的實際意義。PIV技術通過在流場中播撒示蹤粒子，使用相機捕捉激光平面上粒子位移圖像從而計算得到高分辨率二維速度矢量分布，已經在眾多的工程領域得到應用[20-22]，而后發展起來的二維三分量(2D3C)的Stereo-PIV和三維三分量(3D3C)的Tomo-PIV也越來越受到關注[23-24]。然而，這些測量技術也面臨著諸多挑戰，包括測量誤差的影響、測量范圍的限制以及多場同步測量的困難。以Tomo-PIV為例，由于相機景深的限制，測量區域在深度方向上尺寸不能太大，而激光光強限制了其能夠測量的最大流速，相機的數量限制也導致了三維重構中會產生很多幽靈粒子圖像[25]，從而引起較大的測量誤差。

流場中的標量(溫度、濃度、壓力等)分布也是湍流的重要特性，常用的測量手段有溫敏漆(Temperature Sensitive Paint, TSP)[26]、熱敏液晶[27]等壁面溫度場測量，以及激光誘導熒光(Laser Induced Fluorescence, LIF)[28]等流域中溫度或濃度場測量。LIF雖然能夠實現三維標量測量[29]，但同樣受到硬件條件的限制。目前除了壁面的PSP測量技術[30]之外，對于流域內部的高分辨率壓力場還沒有直接測量的方法，而是間接地通過PIV測得的流場和控制方程逆向求解[31]，受到PIV測量區域大小和誤差的嚴重影響[32]。在眾多研究當中，同步獲取全三維的多物理場分布是分析摻混、傳熱、振動、噪聲及化學反應的重要前提，比如發動機燃燒室中燃料的摻混和燃燒，這對測量技術帶來了重大挑戰。目前的成熟方法能夠做到二維兩物理場的同步測量[33-35]，而對于三維或者更多物理場的同步測量還難以實現。

1.2 數值模擬技術

實驗能夠測量的湍流信息是有限的，而更全面的流場信息獲取一般通過數值模擬的方式。由于湍流的多尺度特性，其數值模擬相對于其他物理過程的模擬難度更大。直接數值模擬(Direction Numerical Simulation, DNS)必須解析所有的湍流尺度，而無法應用于工程計算，常用的數值模擬方式為雷諾時均模擬RANS、大渦模擬LES和兩者混合的分離渦模擬DES。

RANS計算量小而被廣泛應用于工程計算，通過使用半理論半經驗的封閉方程(RANS湍流模型)，模化所有湍流尺度，僅求解雷諾時均場，因此，模擬結果的正確性很大程度上取決于RANS模型的表現。RANS模型的不確定度主要來自4個層面[36]：① 雷諾時均過程中湍流信息的丟失，導致湍流對時均場的影響只能通過時均場的相關函數來表示；② 雷諾應力的模化形式，比如線性的Boussinesq渦黏假設將雷諾應力張量模化成渦黏標量場；③ RANS模型方程的選取，也會對湍流模擬造成一系列適用性問題，不同的RANS模型最佳適應場合不一樣，并且也無法準確復現復雜的湍流時均場；④ RANS模型常數都是根據簡單典型的流場進行標定的，無法滿足普適性要求。因此，除了用于模型標定的湍流附面邊界層等流動之外，RANS方法對于復雜湍流場(甚至很多簡單流場)的模擬，都會產生很大的偏差[37]。

RANS的誤差來源于對全部湍流尺度的數學模化，而通過計算網格和N-S方程直接解析部分湍流尺度是提高數值模擬預測準確性的有效手段，因而發展了LES和DES方法。LES直接解析湍流中大尺度旋渦結構，而使用亞格子應力模型模化小尺度旋渦對大尺度結構的影響。在遠壁區網格尺度滿足LES要求[38]的前提下，因為絕大部分湍動能已經通過解析方法直接求得，小尺度結構的影響并不是特別明顯，因此流場分布對亞格子模型并不敏感。而近壁邊界層中湍流尺度非常小，這就大大提高了LES對網格的要求。在邊界層網格分辨率不夠時，LES結果會出現明顯偏差，并影響遠壁區的流場分布[6]。針對此問題，混合DES方法應運而生[5-6]，通過使用RANS方法準確模化附面邊界層流動，而使用LES解析分離和遠壁區的大尺度湍流結構。盡管如此，DES方法仍然存在以下問題而制約其工程應用，包括RANS/LES耦合位置產生灰色區域、耦合之后原來的RANS和LES模型性能都發生了很大的變化、RANS往LES模式轉換延遲、進口邊界難以給定等問題，同時邊界層轉捩也是一個棘手的問題，通過補丁來修正模型方程，而不具備流場模擬的普適性[39-41]。

1.3 前景與挑戰

單純的實驗測量和數值模擬存在各自的問題與缺陷，這2種方法的融合技術備受關注。數據同化使用到的實驗觀測數據和數值預測模型多種多樣，比如在三維流場測量中，通過質量守恒約束來修正速度場的散度[42]也可以看成是數據同化的簡單應用。而在工程應用領域，數據同化技術具有廣闊的應用前景，同時也具有很多技術挑戰。例如，通過二維PIV測量旋轉機械若干截面上的局部流場分布而重建全三維流場和壓力場，對其性能分析及優化設計具有非常重要的意義；通過快響應PSP測得飛行器壁面非穩態壓力分布并融合數值模型，重構空間三維瞬態流場和壓力場，從而獲取聲源信息分析遠場噪聲傳播特性；在高超聲速邊界層模擬中，融合敏感位置壁面壓力傳感器數據和DES模型，使用觀測數據驅動RANS與LES的動態耦合以及實時調整模型參數，顯著提升DES模型對邊界層轉捩的預測精度等等。而針對不同的測量方法，使用怎樣的算法來融合實驗測量與數值模型使實測數據對數值模型產生最有效的控制與約束并降低計算量，以及如何考慮實測數據與數值模型各自的誤差，則是數據同化方法的主要挑戰。

2 基于雷諾時均的湍流穩態數據同化

鑒于RANS方法的工程應用優勢，針對RANS模型進行數據同化研究具有廣闊的應用前景。結合1.2節陳述的RANS模型4個層面的不確定度，數據同化工作主要從以下幾個方面展開。

2.1 模型常數修正

在RANS模型提出初期，一般通過幾種典型的流場來標定模型中的常數，而對于復雜的流場，通過修正模型常數可以獲得模型性能明顯的提升。Li等[43]提出了一種數據驅動的D-DARK(Data-Driven Adaptive RANSk-ω)方法來獲取最優的模型常數。定義RANS模型的常數矩陣為θ，實驗觀測為φo，相對應的模型預測結果為φp，D-DARK的目的為將以下目標函數最小化：

(1)

式中：N為數據點個數；Wi為權函數矩陣。計算J對θ的導數Γ便可迭代求得最優系數矩陣。由于實際中無法獲得φp(θ)函數的顯示表達式，導數的獲取通過單獨擾動每個常數來獲得。對每個模型常數相對于默認參數施加微小擾動，可得導數為

(2)

其中：lθ為模型常數的個數。模型常數便可通過式(3)來更新:

θ(n+1)=θ(n)-ξnΓ(θ)

(3)

其中：ξn為自適應迭代步長。通過多次迭代，模型常數θ便可收斂到J(θ)的極小值點θ*。值得注意的是，此方法并不能使J(θ)收斂到最小值，因此最終預測結果可能與觀測值還存在較大偏差。

D-DARK方法實現起來不太方便，而另一種常用的方法是基于貝葉斯推斷(Bayesian inference)或集合卡爾曼濾波(Ensemble Karman Filter, EnKF)的最優常數估計[13,44-46]。事實上，EnKF可認為是貝葉斯推斷的一個特例，其基本的實現流程如圖1[46]所示。首先對需要優化的常數圍繞某參考值進行擾動并采樣，再使用每組樣本常數進行RANS計算獲得預測空間的樣本分布并與實驗觀測進行對比。EnKF的目的就是確定預測數據與實驗觀測吻合時對應的模型常數。下面以EnKF為例介紹其數據同化原理。

圖1 基于貝葉斯推斷模型常數優化流程圖[46]Fig.1 Flow chart of model constant optimization based on Bayesian inference[46]

在EnKF數據同化中，使用xf表示模型預測的狀態變量，可以理解為更新以后的流場信息；x0為初始狀態變量，即初始流場信息；φ為流場時間演化的非線性函數；v為系統誤差；w為實驗觀測誤差；yexp為實驗觀測；H為狀態變量空間向觀測空間的映射函數。存在以下關系:

xf=φ(x0,v)

(4)

yexp=Hxf+w

(5)

EnKF的算法流程分為預測、分析和確認3個步驟。確認步驟有時可以省略，但對于N-S方程這樣的強非線性預測模型而言，一般建議保留。在預測步驟中，對擾動之后的所有樣本常數有

(6)

(7)

(8)

K=PHT(HPHT+R)-1

(9)

(10)

(11)

式中：R為擾動樣本的協方差矩陣，其定義可見文獻為分析步驟之后的狀態變量樣本；H為狀態空間向觀測空間的映射矩陣。確認步驟與式(10)和式(11)類似，只是對卡爾曼增益增加了一個松弛因子γ以增強計算的穩定性，即

(12)

(13)

其中：j為確認步驟中的迭代次數。經過EnKF同化之后的模型常數樣本會分布在一個非常窄的范圍內，進行樣本平均之后即可得到最優化的模型系數矩陣。

通過不同RANS模型和EnKF數據同化常數優化獲得的圓管射流結果如圖2和圖3所示[47]。射流的雷諾數為Re=6 000, PIV實驗在水槽中進行，入口條件為充分長的圓管流動。圖2展示了通過不同的RANS模型得到的截然不同的射流時均場，圖中U為射流時均流向速度，U0為噴口流速，D為噴嘴直徑，r/D為徑向位置。k-ω模型預測的射流時均場與PIV測量結果的偏差最明顯，但是其他的RANS模型結果偏差也比較大。圖3顯示，通過EnKF數據同化之后，所有RANS的預測結果均得到了很大的改善，尤其是k-ω模型。這說明，對于特定流場下標定的RANS模型常數，并不一定適用于非標定的流場情形，但是對于多種流場，可以通過重新標定模型系數來降低模型的不確定度。基于EnKF的數據同化方法提供了一種方便可行的常數標定方法，當然其他的一些方法[43, 48]也能達到類似的目的。

圖2 不同RANS模型預測的射流時均場與PIV測量結果的對比[47]Fig.2 Mean flows of circular jet predicted using different RANS models and comparison with PIV results [47]

圖3 基于不同RANS模型的常數數據同化結果及與PIV結果的對比[47]Fig.3 Data assimilation results based on constant optimization of different RANS models and comparison with PIV results[47]

2.2 模型形式誤差修正

對于特定的RANS模型，其誤差來源除了模型中的調節常數以外，另一個重要方面是模型固有的方程形式。對于湍動能k的計算，其確切的控制方程為

(14)

式中：U為速度矢量場；T′、P和ε分別為方程的傳輸項、生成項及耗散項，且T′和ε均與湍流的脈動統計特性有關，使得方程無法封閉。在k-ε方程中，使用梯度-擴散假設來對T′進行模化

(15)

其中：νt為湍流渦黏；σk為需要標定的模型常數。而對ε模化的整個方程則被完全經驗式地認為與k形式相似。這種模化的表達形式給RANS模型帶來了很大的方程形式誤差[49]。在這方面工作中，Oliver等[50-51]基于貝葉斯估計在Boussinesq渦黏假設中增加了一個修正項，用來彌補RANS模型預測的誤差，但這一修正項的選取及對觀測數據的依賴關系，是該方法的關鍵。

另一種常用的方法是針對渦黏模型，在湍流生成項中增加空間(用x表示空間位置)分布的修正系數β(x)，以修正模型的形式誤差，而最優的β(x)分布可通過貝葉斯反演算法來獲得[52]。以Spalart-Allmaras (SA)模型[53]為例, 有

(16)

對于最優修正系數β的獲取，更常用的一種方法是采用伴隨方程的反演算法[54]。定義目標函數

(17)

(18)

(19)

通過式(18)可計算得到全導數，從而使用梯度下降法逐步逼近獲得最優的β分布。Singh等[54]使用此方法對光滑壁面流動分離做了數據同化研究，選用壁面壓力系數為實驗觀測進行數據同化，在反演修正系數β之后，結果相對于默認的SA模型有了非常明顯的提升。值得注意的是，此伴隨方程完全通過離散的矩陣計算得到，需要很大的計算量且不方便實現，因而到目前為止，這種方法僅在二維流場計算中得到了應用。

為解決離散伴隨計算量大、實施困難等問題，作者團隊[16]提出了連續伴隨數據同化(ABDA)模型。連續伴隨方法早期被應用于拓撲優化研究，相比離散伴隨方法表現出更高的工程應用價值[55]。為獲取最優的修正系數β分布，定義拉格朗日函數

(20)

(21)

選擇合適的伴隨變量，使得式(21)的后3項為0，即

(22)

可求得伴隨方程組為

(23)

(24)

(25)

對應的變量含義及邊界條件見文獻[16]。ABDA模型建立了方程驅動的連續伴隨系統，聯立主方程組(N-S方程、連續性方程和SA模型方程)和伴隨方程組(式(23)～式(25))迭代求解即可獲得最優β分布。相對于離散伴隨方法，連續伴隨方法避免了大型矩陣計算，更能適用于復雜三維流場的數據同化，如圖4和圖5所示的三維方柱分離與再附流動[16]，其中使用了熱線測得的若干位置流向速度Uexp為觀測量，圖中U0為來流速度。

圖4 三維方柱計算域和實驗觀測數據的植入[16]Fig.4 Computational domain of 3D cube flow and implantation of experimental data[16]

圖5 三維方柱流動反演得到的流場和β分布[16]Fig.5 Inversed flow and β distribution of 3D cube using data assimilation[16]

2.3 雷諾應力項修正

上述的模型修正方法均建立在Boussinesq渦黏假設的框架之內，盡管使用空間分布的β修正系數可以顯著減小模型方程的形式誤差，而其適用性仍受到標量渦黏νt的制約。為解決此問題，需拋開渦黏假設，發展對雷諾應力項修正的數據同化方法。Foures[56]和Symon[57-58]等使用連續伴隨直接對雷諾應力項F進行同化，

(26)

(27)

對雷諾應力的數據同化方法雖然可以不受渦黏假設的約束而獲得非常好的結果，但是最明顯的缺點是計算不穩定，以致對高雷諾數的流場同化幾乎無法實現。作者團隊對此方法進行了改進，提出了基于各項異性的湍流數據同化方法。

(28)

式中：νeff為流體分子黏度ν與湍流渦黏中各項同性部分νt0之和；與式(26)不同的是，這里F僅表示各向異性渦黏部分νta的貢獻：

(29)

其中： “.*”為向量對應元素相乘符號。F通過使用連續伴隨方程進行同化，而νeff則直接通過常規的RANS渦黏模型進行模化。F的分布會根據RANS模型的選取而有所變化，但是對于任何流場，均能同化獲得最優F分布來復現真實湍流中雷諾應力的貢獻，即

(30)

此方法相對于式(26)最明顯的優勢為，在稍微增加常規RANS方程計算量的前提下，利用各向同性渦黏明顯提高了計算收斂的穩定性，而方法的適用性與式(26)一致。

圖6 圓管的各向異性數據同化結果Fig.6 Results of anisotropic data assimilation of a circular jet

與DA-SA時均流場進行對比，其速度等值線幾乎重疊。而如圖6(b)～圖6(c)所示，該數據同化方案中，雷諾應力項的還原非常準確，這相比于常規的SA模擬和各項同性ABDA模型(式(23)～式(25))有非常明顯的改善。值得注意的是，式(16)所示模型方程形式的修正方案雖然能夠顯著提升流場的還原準確度，但是并不是建立在對雷諾應力項的準確還原之上，因此即使流場有顯著改善，但是雷諾應力項可能與真實值(LES)間的偏差依然很大。

3 湍流的非穩態數據同化

3.1 4DVar變分數據同化

4DVar是一種基于伴隨的非穩態數據同化方法，它的目的是根據一套初始狀態值和實驗觀測來改善模型的預測軌跡，并通過當前和過去的觀測數據準確預測未來狀態。其基本的方程為[59]

(31)

(32)

Yt=H(Xt(x))+t

(33)

J(η,qt)=J0+Jt

(34)

(35)

(36)

其中：B為背景誤差協方差矩陣；R為觀測誤差協方差矩陣。最小化目標函數J(η,qt)通過伴隨系統計算J沿δn=(δqt,δη)方向的導數來實現，而最小化的評估區間為[t0,tf]時間段。對控制方程進行時間積分并引入伴隨變量λt：

(37)

其中：?XM為對應于非線性算子M的切線線性算子，引入伴隨算子(?XM)*和(?XH)*滿足：

(38)

(39)

并對式(37)進行部分積分可得

(40)

由于伴隨變量滿足

(41)

λtf=0

(42)

將伴隨方程式(41)代入式(40)并使用dXt0=δη及dXt=(?Xt/?n)δn可得目標函數的導數為

(43)

(44)

式中：Q為預測模型誤差的協方差矩陣，一般認為誤差為中心高斯隨機分布。

通過式(43)～式(44)求得的目標函數導數可以迭代確定最優的η和qt分布。具體的求解可采用差分-離散或者離散-差分的方法。前者為連續伴隨方法，需要推導完整的伴隨方程組和邊界條件，而且使用到一些假設條件；而后者避免了方程的推導，直接通過數據驅動來進行系統求解，精度更高，但是如前所述，其計算量更大。由于伴隨變量的λt的初始值未知，而只知道其終值λtf=0，這就意味著求解式(40)所示的伴隨方程需要進行時間逆向積分，而對于主控制方程式(31)進行時間正向積分。由此，求解整個方程系統需要保存所有計算時間步下的結果以完成正向和逆向時間積分，并以此為一個迭代步數，整個系統經過若干個迭代步直至方程組完全收斂。4DVar變分同化是使用最廣泛的數據同化方法之一，國內在數值天氣預報、地理、核電等領域都有廣泛的研究及應用[60-63]，但尚未涉及到湍流流場的深入研究，國外近幾年開始了對湍流流場的4DVar數據同化研究[64-67]。Mons等[64]對非穩態來流情況下圓柱繞流進行了4DVar同化，取圓柱壁面的壓力系數為實驗觀測，使用4DVar數據同化獲取最優的初始場及非穩態入流條件參數。結果顯示，經過同化之后的流場和壓力分布準確反映了周期性非穩態來流的影響。

3.2 順序數據同化

4DVar數據同化實現了在整個時間段內的參數最優化分布，而最大的問題是需要進行時間正向和逆向積分，需事先獲取整個時間段內的觀測數據及保存每時間步的計算結果，存儲量大且計算復雜。順序數據同化(Sequential DA)實現了時間上的正向順序推進，使用實時動態數據修正模型預測軌跡，其計算過程與常規的數值模擬類似，計算量相對4DVar顯著降低。圖7所示為順序數據同化的原理示意圖，通過在整個時間段內的正向和逆向積分迭代，4DVar數據同化獲得了連續的時間序列并逼近觀測數據。為降低計算量和存儲量，順序數據同化通過時間正向積分，同樣使同化結果逼近實驗觀測，但是在有觀測數據約束的時刻，通過對模型預測軌跡進行修正，使得結果出現不連續現象。這種不連續現象隨著觀測數據時間間隔的減小而減弱，當觀測時間間隔等于計算時間步長時，同化結果理論上與4DVar一致。

圖7 順序數據同化原理圖Fig.7 Schematic diagram of sequential data assimilation

Meldi和Poux[68]使用卡爾曼濾波構建了順序數據同化方法，其在計算流體力學底層PISO (Pressure Implicit with Splitting of Operators) 算法中，將半離散形式的動量方程作為預測模型。

(45)

(46)

此方法的關鍵在于卡爾曼增益K的確定，如式(9)所示，K與實驗測量誤差、模型預測誤差的協方差矩陣直接相關，這些協方差矩陣一般可以通過解析或者統計的方法來確定。統計的方法需要計算大量樣本[69-70]，這給順序數據同化方法帶來了很大的困難，而解析的方法需要計算預測模型誤差傳遞特性，這對于N-S方程來說難以實現。

上述方法基于強烈假設來構建卡爾曼降階模型，從而確定卡爾曼增益，其適用性和同化的準確性還有待提升。作者團隊[32]基于4DVar方法進行改進，提出了基于連續伴隨的順序數據同化方法，避免了時間正向和逆向積分，較4DVar方法節省約90%的存儲量。其基本原理與式(28)類似，在N-S方程中加入同化控制變量F

(47)

(48)

式中：流體黏度ν也可以用等效黏度νeff替代，從而引入亞格子尺度湍流模型。湍流渦黏可以使用LES或DES模型來獲得。此問題的拉格朗日函數為

(49)

其中：與4DVar相同，[t0,T]為數據同化的時間窗口；為控制方程式(47)和式為在時間窗口[t0,T]中評估的目標函數，即

(50)

(51)

(52)

V(T)=0

(53)

其中：ρ為流體密度。

連續伴隨方法通過推導確切的偏微分方程組及邊界條件(見參考文獻[32])降低了數據計算量。為了避免時間的逆向積分，將時間窗口縮短為一個計算時間步長，即[t0,t1]，并使用歐拉格式對式(51)中的時間偏導進行差分，從而得到

(54)

式中：Δt為計算時間步長。變換之后的伴隨方程式(54)消除了時間非穩態項，避免了時間逆向積分，在存在觀測數據的時刻，可隨主控制方程組同步求解，從而實現流場的順序數據同化，其計算流程如圖8[32]所示。圖9[32]展示了射流流場的數據同化結果，觀測區域選為0.2

圖8 連續伴隨順序數據同化計算迭代流程[32]Fig.8 Iteration procedure of adjoint-based sequential DA[32]

圖9 射流流場的伴隨順序數據同化[32]Fig.9 Adjoint-based sequential DA of jet flow[32]

3.3 基于物理信息的機器學習

機器學習方法在湍流研究中快速發展，推動了湍流的預測往人工智能領域邁進。數據驅動的機器學習模型直接通過流場數據本身進行學習訓練，可以很好地構建相關物理量之間的映射關系，得到了廣泛的研究與應用。此類模型算法在湍流特征提取[71-72],流場預測與重構[73-76],邊界層分類[77],湍流模型優化[78],流場的模態識別[79]與時空分辨率增強[80-82]等諸多領域取得了很好的應用效果。但此類模型本身缺乏相應物理方程的約束，模型內部的權重與偏置系數的獲取需要進行大量的數據訓練擬合，同時模型的泛化能力有待提升。

近年來，基于物理信息的機器學習模型備受研究學者關注。Raissi等[83-84]提出基于物理信息的機器學習模型率先將傳統的N-S方程與人工神經網絡相結合并應用于渦激振動的圓柱繞流中，成功地從流場的某單一物理量信息(如速度場或濃度場)中提取出與之對應的其他物理量信息(如壓力場)。其設計的基本網絡結構如圖10[84]所示, 其中，左邊灰色區域部分為傳統的全連接神經網絡，網絡輸入為流場的時間與空間信息，輸出為相對應的物理量信息，如速度(u,v),壓力p和圓柱振動位移η。白色圓圈代表偏微分算子，用于求解計算相應物理量的時間與空間偏微分。右邊區域為N-S方程的物理損失項約束，用于約束模型訓練，使得模型預測結果滿足物理方程的規律。

圖10神經網絡中的u、v、p、η這些物理量信息可以通過求解損失函數Jloss(式(55))的極小值獲得。其中Jloss的前2項分別為實測的流場和圓柱位移與模型預測的流場和圓柱位移的誤差，Jloss的最后一項為物理誤差損失；當e1、e2、e3的值在模型訓練中逼近于0時，則意味著相關物理量u、v、p、η滿足了特定物理方程的約束。

|v(tn,xn,yn)-vn|2)+

(55)

現階段，基于物理信息的機器學習方法還多局限于層流的預測，但是其基本的數學原理和約束方程同樣適用于湍流問題。從基本數學思維上講，基于物理信息的機器學習方法同樣是為了最小化目標函數，且圖10中物理約束方程與連續伴隨方程(式(49))中的控制方程具有異曲同工之處。隨著技術的發展，基于物理信息的機器學習方法也將在湍流問題預測和重構中發揮重要作用。

圖10 N-S方程驅動的神經網絡模型[84]Fig.10 N-S equation informed neural networks[84]

4 湍流數據同化的應用

4.1 物理約束的輔助測量

數據同化技術作為一種實驗測量與數值模擬的耦合方法，在湍流研究領域發揮著越來越重要的作用。數據同化的主要應用之一可以概括為基于物理約束的輔助測量。在控制方程的物理約束下，數據同化輔助完成實驗測量，獲得更全面的流場數據。通過PIV數據來同化獲得全場的流場信息是比較常見的應用，比如Dovetta等[85]通過部分區域PIV測量的時均流場獲取翼型周圍全場的流場分布，其還通過變分數據同化方法并以PIV為觀測，獲取了管內的全場流動信息，并進一步準確計算了剪切應力分布[86]。作者團隊[87]展示了數據同化能夠修正PIV測量的誤差，獲得更加準確的近壁方柱繞流全局流動分布。如圖11[87](圖中d為方柱截面的尺寸)所示，由于分辨率的原因，方柱下方的近壁射流區域無法被PIV準確測得，因而產生了很大的測量誤差，而這被數據同化準確地還原了。除此以外，數值同化也可以用來預測流場中無法測量的衍生物理信息，比如溫度場、濃度場、壓力場等。

基于PIV流場測量的壓力場計算是數據同化的重要應用[32, 66]。常規的方法是通過N-S方程求得壓力梯度之后，再對壓力梯度積分求得壓力分布[31]，但一方面要求獲得三維的流場，另一方面受PIV測量誤差影響，壓力重構精度并不高。圖12[32]為射流的壓力場重構，以LES模擬的三維湍流流場為實驗觀測，數據同化獲得的流場在圖9中給出，可以看到非穩態壓力場也被準確還原。而常規的壓力梯度積分重構并沒有得到比較好的結果，原因是在積分中容易導致誤差積累，即使此案例中使用到的觀測為LES純凈流場，圖12(c)顯示數值誤差的積累也是一個非常嚴重的問題。

鑒于梯度積分法中需要對流場求時間導數，這就需要先獲得時間解析的非穩態流場。而當流場的采樣頻率較低時，計算得到的壓力脈動通常會偏大[88]。在數據同化中，不需要每一時間步都進行流場觀測數據的植入，這就使得實驗測量中流場的采樣時間間隔可以100倍(甚至更高)于求解時間導數所需要的流場采樣間隔。最近，作者團隊[89]結合線性隨機估計(Linear Stochastic Estimation, LSE)，實現了從1 Hz的PIV測量數據同化獲取1 kHz的非穩態壓力分布。其首先使用1 Hz的PIV測量和1 kHz的動態壓力傳感器同步測量，獲得了粗糙的1 kHz流場數據，再通過數據同化進行誤差過濾，準確重構了1 kHz的流場和非穩態壓力場。

4.2 數據驅動的數值模擬

常規的數值模擬結果與實驗結果一般存在比較大的誤差，結合數據同化技術，使用少量實驗結果作為觀測，可以有效提高模型預測的準確性。Meldi和Poux[68]在平板分離再附流場模擬中，使用基于卡爾曼濾波的數據同化方法在平板上方截面植入若干時刻的觀測數據。圖13[68](圖中Hp為平板高度)中由于觀測數據非常有限，數值模擬提升的幅度不大，這可以增加觀測數據植入的時間頻率來進一步改善模型預測。作者團隊[16]提出的ABDA模型根據熱線測得的若干位置速度分布降低湍流模型的不確定度。Singh等[54]使用類似的方法，通過數據同化反演獲得翼型不同工況下的模型修正系數之后，結合機器學習預測任意翼型工況下的修正系數分布，從而建立了離線的數據驅動數值模擬方法(如圖14[54]所示，圖中C為翼型弦長)。

圖13 平板繞流數據驅動模擬中不同下游位置的時均流向速度和橫向速度分布[68]Fig.13 Time-averaged distributions of streamwise velocity at different downstream locations in data-driving simulation of flow over bluff plate[68]

圖14 翼型表面壓力系數Cp的數據同化反演與機器學習預測[54]Fig.14 Data-assimilation inversing and machine-learning simulation of pressure coefficient distribution Cp on a aerofoil[54]

使用數據同化技術也可以為DNS或LES數值模擬提供進口邊界條件，從而顯著縮小模擬區域，簡化幾何模型。Gronskis等[65]對圓柱繞流做了DNS模擬，其計算域為不包含圓柱在內的尾跡區域，同時使用4DVar技術，結合PIV測量結果，對DNS的入口速度分布進行數據同化。從其結果看出，在獲取進口速度條件之后，DNS對下游尾跡區域流場的模擬與PIV測量結果非常吻合，同時相比PIV結果具有更高的分辨率和更低的噪聲水平。

4.3 基于大數據的流場預測

隨著數據同化技術的發展，流場數據的精度和廣度不斷增加，而歷史數據庫的積累，為湍流的智能化鋪平了道路。未來的世界是大數據的世界，數字孿生[90]也成為了湍流研究的發展方向，而更高效更準確地使用這些數據來預測未來的狀態，是這項研究中的關鍵問題。Renganathan等[91]嘗試了面向數字孿生范式的空氣動力學數據融合，基于貝葉斯推斷，結合數值模擬數據庫和實測數據，重建機身表面高分辨率壓力分布。其重建方案如圖15所示，具體的變量定義和算法策略見參考文獻[91]，這里不做詳細展開。類似的研究還包括眾多基于機器學習的流場預測[92-94]，而這些僅僅是邁出了大數據時代的一小步。純數據驅動的人工智能使湍流喪失了物理本質，而物理信息的機器學習和數據同化技術在這一點上本質相同。

圖15 基于傳感器數據和貝葉斯推斷的 飛行器氣動力學預測[91]Fig.15 Aerodynamics prediction of aircrafts based on probe data and Bayesian inference[91]

5 結 語

本文綜述了現階段湍流的主要數據同化技術及應用。數據同化作為實驗測量和數值計算的耦合方法，彌補了單純實驗測量和數值計算中的缺點和不足，提高了實驗測量的精度和廣度，改善了數值模擬的預測性能。

穩態數據同化一般結合RANS模型方程，從重新標定模型常數、修正渦黏模型方程形式誤差、修正雷諾應力項等方面著手。基于EnKF的模型常數標定方法是一種顯著提升模型預測精度的方便高效的穩態數據同化方法，而其方程的形式決定了其有限的適用性。基于伴隨的數據同化能夠對模型方程的形式進行修正，確定最優的修正系數空間分布，但是其仍處于渦黏假設的框架之內，對于某些復雜流場并不適用。而對雷諾應力項進行修正的數據同化適用性更廣，理論上沒有特定流場的局限。

非穩態的數據同化一般包括4DVar等時間連續的數據同化方式以及順序數據同化。4DVar通過時間正向和逆向積分迭代，在整個時間段內對預測模型進行修正，但其缺點是存儲量和計算量都非常大。順序數據同化不需要時間逆向積分，可以在若干時刻對實驗觀測進行間斷性地植入，不斷地修正模型運行軌跡，其計算量和存儲量與常規數值模擬相似。

數據同化在促進湍流研究發展方面具有非常重要的意義，成為融合實驗測量和數值計算的重要手段，開拓了物性約束的輔助測量及數據驅動的數值模擬兩條研究道路。隨著機器學習、人工智能的飛速發展，湍流研究也向智能化邁進。基于物理信息的機器學習雖然現階段還主要應用于層流預測，但是鑒于其基本的數學原理，也將在湍流領域獲得廣泛的應用前景。