何創(chuàng)新，鄧志文，劉應征,*

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

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

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

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

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

1 湍流測量與數(shù)值模擬的挑戰(zhàn)

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

1.1 實驗測量技術

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

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

1.2 數(shù)值模擬技術

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

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

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

1.3 前景與挑戰(zhàn)

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

2 基于雷諾時均的湍流穩(wěn)態(tài)數(shù)據(jù)同化

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

2.1 模型常數(shù)修正

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

(1)

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

(2)

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

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

(3)

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

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

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

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

xf=φ(x0,v)

(4)

yexp=Hxf+w

(5)

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

(6)

(7)

(8)

K=PHT(HPHT+R)-1

(9)

(10)

(11)

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

(12)

(13)

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

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

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

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

2.2 模型形式誤差修正

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

(14)

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

(15)

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

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

(16)

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

(17)

(18)

(19)

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

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

(20)

(21)

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

(22)

可求得伴隨方程組為

(23)

(24)

(25)

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

圖4 三維方柱計算域和實驗觀測數(shù)據(jù)的植入[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渦黏假設的框架之內(nèi)，盡管使用空間分布的β修正系數(shù)可以顯著減小模型方程的形式誤差，而其適用性仍受到標量渦黏νt的制約。為解決此問題，需拋開渦黏假設，發(fā)展對雷諾應力項修正的數(shù)據(jù)同化方法。Foures[56]和Symon[57-58]等使用連續(xù)伴隨直接對雷諾應力項F進行同化，

(26)

(27)

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

(28)

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

(29)

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

(30)

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

圖6 圓管的各向異性數(shù)據(jù)同化結(jié)果Fig.6 Results of anisotropic data assimilation of a circular jet

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

3 湍流的非穩(wěn)態(tài)數(shù)據(jù)同化

3.1 4DVar變分數(shù)據(jù)同化

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

(31)

(32)

Yt=H(Xt(x))+t

(33)

J(η,qt)=J0+Jt

(34)

(35)

(36)

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

(37)

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

(38)

(39)

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

(40)

由于伴隨變量滿足

(41)

λtf=0

(42)

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

(43)

(44)

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

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

3.2 順序數(shù)據(jù)同化

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

圖7 順序數(shù)據(jù)同化原理圖Fig.7 Schematic diagram of sequential data assimilation

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

(45)

(46)

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

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

(47)

(48)

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

(49)

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

(50)

(51)

(52)

V(T)=0

(53)

其中：ρ為流體密度。

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

(54)

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

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

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

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

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

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

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

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

(55)

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

圖10 N-S方程驅(qū)動的神經(jīng)網(wǎng)絡模型[84]Fig.10 N-S equation informed neural networks[84]

4 湍流數(shù)據(jù)同化的應用

4.1 物理約束的輔助測量

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

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

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

4.2 數(shù)據(jù)驅(qū)動的數(shù)值模擬

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

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

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

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

4.3 基于大數(shù)據(jù)的流場預測

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

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

5 結(jié) 語

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

穩(wěn)態(tài)數(shù)據(jù)同化一般結(jié)合RANS模型方程，從重新標定模型常數(shù)、修正渦黏模型方程形式誤差、修正雷諾應力項等方面著手。基于EnKF的模型常數(shù)標定方法是一種顯著提升模型預測精度的方便高效的穩(wěn)態(tài)數(shù)據(jù)同化方法，而其方程的形式?jīng)Q定了其有限的適用性。基于伴隨的數(shù)據(jù)同化能夠?qū)δＰ头匠痰男问竭M行修正，確定最優(yōu)的修正系數(shù)空間分布，但是其仍處于渦黏假設的框架之內(nèi)，對于某些復雜流場并不適用。而對雷諾應力項進行修正的數(shù)據(jù)同化適用性更廣，理論上沒有特定流場的局限。

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

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