基于數據驅動的流場控制方程的稀疏識別1)

江 昊 王伯福 盧志明

(上海大學,力學與工程科學學院,上海市應用數學和力學研究所,上海 200072)

引言

非線性動力系統廣泛存在于自然界及工業界中,對非線性動力系統的建模是亟待解決的問題.非線性動力系統(如電力網絡[1]、金融學[2]、神經科學[3]等)的控制方程大多是非線性偏微分方程,從傳統的守恒定律和物理原理等第一性原理出發去推導一些復雜系統的控制方程是困難的.現代流體力學雖然已經具備了較完備的理論方程和模型,但仍然有很多復雜問題的理論研究尚不完備,如帶化學反應流、多相流、非牛頓流、稀薄流等,這些系統的動力學方程或模型的理論推導難以實現[4].近些年來,隨著計算機硬件以及計算科學的發展,基于數據驅動方法對非線性復雜系統建模得到了研究者的青睞[5-13].其中稀疏識別方法[14]可以簡化復雜的數學模型,得到項數較少的簡單模型來表征系統的動力學特性,并可以綜合考慮過擬合和模型精確性之間的矛盾.稀疏識別方法已經被廣泛用于復雜非線性系統控制方程的識別[14-21].1996 年,Tibshirani[14]首次提出最小絕對收縮和選擇算子(least absolute shrinkage and selection operator,LASSO)方法,并且把LASSO 方法應用在Cox模型的稀疏識別中[15].2016 年,Brunton 等[18]通過稀疏識別在含有噪聲的非線性動力系統的數據中識別出了常微分方程.2017 年,Rudy 等[19]提出了偏微分方程函數識別(partial differential equations functional identification of nonlinear dynamics,PDE-FIND)方法可以通過測量或者數值模擬的時空數據來稀疏識別偏微分方程,正確地識別出了圓柱繞流流場的渦量輸運方程.胡軍等[20]將貝葉斯稀疏識別方法識別結果與PDE-FIND 方法識別結果相比較,表明了貝葉斯稀疏識別方法對偏微分方程具有非常強的稀疏恢復能力,但相比PDE-FIND 方法對噪聲更加敏感.

近年來,數據驅動方法在流場中的應用頗受研究者[4,22-27]的重視.2019 年,Chang 和Zhang[23]使用LASSO 方法從數據出發,識別出了單相地下水流動方程和污染物輸運方程.2020 年,Raissi 等[24]使用深度學習,對流場可視化的圖片進行學習,并反演了流場的壓力和速度分布場,可以幫助得到實驗難以測量的數據.Zhang 和Ma[25]采用PDE-FIND 方法,對基于直接模擬蒙特卡羅(direct simulation Monte Carlo,DSMC) 方法模擬產生的數據進行稀疏識別,反演了流體的對流、擴散以及渦量輸運方程,從數據驅動角度清晰地展示了玻爾茲曼方程和Navier-Stokes 方程之間的相關性.Schmelzer 等[26]采用SINDy 框架在數值模擬的數據中是識別出了代數雷諾應力模型.張亦知等[27]構建了基于數據驅動的湍流模型修正框架并在槽道流中進行了驗證.

本文采用基于數據驅動的稀疏識別方法(PDEFIND 方法和LASSO 方法)對圓柱繞流、頂蓋驅動方腔流、RB 對流和三維槽道湍流的控制方程進行稀疏識別,研究了不同流場,不同控制方程以及不同區域數據對于流場控制方程稀疏識別的影響.

1 流動控制方程的識別方法

非線性系統的控制方程大多是非線性偏微分方程(組),例如Burgers 方程(1),科爾特弗-德弗里斯(korteweg-de Vries,kdV)方程(2)

這些非線性偏微分方程可以寫成通用形式

其中,u表示時空變量,ut表示變量u的時間導數,下標x的個數表示變量u的空間n階導數,μ 表示系統的參數,N(·)表示時空變量u(x,t)及其時空各階導數的非線性函數.在非線性偏微分方程的稀疏識別中,將包含變量u及其導數相關的項組成矩陣U,將其他可能出現的項組成矩陣Q.將U,Q合并成矩陣Θ

對于離散的數據點,U和Q分別為(N×Mu) 和的二維矩陣,其中N=n×m為矩陣中每項對應的時空點數(n為空間網格總數,m為時間步數),Mu和Mq分別為U和Q的項數.將非線性偏微分方程表示為以下線性方程

其中,Θ是人為給定的一個過完備的候選項的庫,即Θ庫中包含所有可能屬于該非線性偏微分方程的項.ξ是稀疏向量,表示的是非線性偏微分方程系數矩陣,該稀疏向量ξ的非稀疏項等于該非線性偏微分方程的系數,則偏微分方程可以寫成線性加權和的形式

流體力學中非線性偏微分方程的稀疏識別就是通過流場的部分空間位置的時間序列數據,來找到稀疏系數矩陣ξ使得方程(6)滿足該流場動力學行為.時間序列數據可以通過實驗或者數值模擬來獲得,變量的時空導數通過多項式插值法或有限差分法獲得.對于不含有噪聲的數據通常采用有限差分法計算獲得,但對于有噪聲的數據,多項式插值法計算的導數比有限差分法計算的導數更精確[19].

如果Θ內所有的項都屬于該流場系統的控制方程,此時可以使用最小二乘法確定系數矩陣

然而,在過完備的庫Θ中通過最小二乘法求解系數矩陣ξ,會引起過擬合問題,解出的系數矩陣ξ中所有元素都是非零,這樣所識別出的非線性偏微分方程很復雜,比流場精確控制方程多出很多項,不符合流體動力學規律.

對于偏微分方程的稀疏識別,一種比較直觀的方法就是將所有項可能的組合都進行計算,比較所有項組合的誤差,得到最優的項組合,數學表示為加入L0正則化最小二乘法

但是求解該問題會遇到NP-hard 問題,即把所有的組合都進行計算和比較,這對于計算資源的消耗是巨大的.Tibshirani[14]提出的LASSO 算法是引入L1正則化的最小二乘法

該方法是對于L0正則化最小二乘法的一個凸優化問題,避免了L0正則化的NP-hard 問題,圖1 是LASSO方法的流程圖.但LASSO 方法對列相關性很高的數據矩陣的稀疏識別比較困難,例如對于KdV 方程的識別[19].

圖1 LASSO 方法流程圖Fig.1 Flow chart of LASSO scheme

Rudy 等[19]提出了一個新的稀疏識別方法:PDEFIND 方法.該算法是在加入L2正則化的最小二乘法(1) 中加入閾值篩選,使用Ridge 算法[28]求解L2正則化的最小二乘法,將小于閾值的系數所對應的項篩去,對剩余的項重復遞歸Ridge 算法求解和閾值篩選,直到所有項的系數都高于閾值,此時返回稀疏系數矩陣ξ,該過程被稱為STRidge 算法

顯然,閾值的選取對于STRidge 算法很重要,不同的閾值會導致識別的偏微分方程的稀疏程度不同,因此為了找到最佳的閾值,Rudy 等[19]提出了一個訓練閾值的算法(TrainSTRidge 算法).TrainSTRidge 算法和STRidge 算法相互耦合一起形成了PDE-FIND方法.圖2 是PDE-FIND 方法的流程圖.

圖2 PDE-FIND 方法流程圖Fig.2 Flow chart of PDE-FIND scheme

2 識別結果

2.1 圓柱繞流渦量輸運方程的識別

Rudy 等[19]使用PDE-FIND 方法對圓柱繞流的渦量輸運方程進行了稀疏識別,在未添加噪聲的流場數據中識別出的渦量輸運方程的平均系數誤差為1%±0.2%.本文同樣也對圓柱繞流問題的渦量輸運方程進行了識別.

采用LBM 方法對雷諾數為300 的圓柱繞流進行了直接數值模擬,計算網格為1000×250,圓柱圓心位置為(201,133),半徑為26,來流最大速度為u=0.1.圖3 給出了t=104時流場的渦量云圖.如圖3 所示將圓柱繞流后面的區域分為x∈ (250,500),x∈ (500,750),每個區域隨機選擇8000 個空間點并追蹤記錄60 個時間步長的速度數據和壓力數據,這樣每個區域就共獲得了48 萬個數據樣本,這對于稀疏識別所需的樣本是充足的,且采樣的方法不會影響識別的結果[19].采用5 次切比雪夫多項式插值[29]計算導數,導數階數最高保留到二階,變量次數最高保留到二階,非線性項最高保留到四階,將這些數據用來構建稀疏識別所需Θ庫(4).表1 給出了PDE-FIND 方法和LASSO方法對于x∈(250,500) 圓柱繞流數據的稀疏識別結果.

圖3 t=104 時圓柱繞流的渦量云圖Fig.3 Vorticity contour of flow past a cylinder at t=104

表1 Re=300,x ∈(250,500)圓柱繞流控制方程的識別Table 1 Identification of vorticity transport equation for flow past a cylinder with Re=300 in x ∈(250,500)

從表1 可得出,PDE-FIND 方法和LASSO 方法都正確地識別出了渦量輸運方程,并且系數相同.兩種算法對于渦量輸運方程的稀疏識別的能力一致,且平均系數誤差為3.21%.渦量輸運方程是一個線性偏微分方程,PDE-FIND 方法和LASSO 方法對線性偏微分方程的稀疏識別能力一致.對于x∈(500,750) 區域的數據,兩種算法也都識別出了正確的控制方程,但所得平均系數誤差有所增加,為8.75%.這是因為越靠近圓柱區域的流場信息越豐富,因而識別出的控制方程系數越精確.這與文獻[20]的結論是吻合的,即選取具有豐富動力學行為的時空演化數據可以提高對動力學控制方程的識別的準確性.

2.2 頂蓋驅動方腔流控制方程的識別

本節對二維頂蓋驅動方腔流的NS 方程組進行稀疏識別.相比于上一節渦量輸運方程的識別,原始變量形式的NS 方程組是一組偏微分方程組

其中,方程(12)是連續性方程,方程(13)是動量方程.二維頂蓋驅動方腔流場數據使用基于有限差分的直接數值模擬得到.方腔寬高比為2,方腔頂部的壁面以大小為1 的速度勻速向左運動,流動雷諾數為100.模擬計算網格數為200×100,圖4 給出了t=40 時流場的壓力云圖和流線圖.

圖4 t=40 時頂蓋驅動方腔流的壓力云圖和速度矢量圖Fig.4 Pressure contour and velocity vector of lid-driven cavity flow at t=40

注意到連續性方程(12),不含有時間偏導項,將連續性方程可改寫為

將x看成時間t,就可以采用稀疏識別方法來識別連續性方程了.

在整個方腔區域選取數據,采樣方法、導數計算方法以及Θ庫的構建與2.1 節相同.表2 給出了PDE-FIND 方法和LASSO 方法對于Re=100 頂蓋驅動方腔流控制方程的稀疏識別結果.從表2 可得出,對于頂蓋驅動方腔流的數據,PDE-FIND 方法正確地識別出了連續性方程和NS 方程的各項,且平均系數誤差為2.77%,而LASSO 方法只能正確地識別出連續性方程,沒有識別出動量方程,例如在u方程中未選擇出uux,可能是因為LASSO 算法對含強非線性項的方程的識別會失效,具體的原因將在下一節進行討論.

表2 Re=100,頂蓋驅動方腔流控制方程的識別Table 2 Identification of vorticity transport equation for lid-driven cavity flow with Re=100

2.3 Rayleigh-B′enard 對流控制方程的識別

下面考慮識別Rayleigh-B′enard(RB)對流系統中的控制方程,其物理模型是一個下底板加熱,上底板冷卻的對流系統,系統中充滿流體介質,由于上下底板溫度不一樣,導致不同高度位置的流體之間存在密度差,在浮力的驅動下系統中流體產生運動.在OberbeckBoussinesq (OB) 近似下的無量綱控制方程組如下

其中,方程(15)是連續性方程,方程(16)是含有熱浮力項的NS 方程,方程(17)是熱輸運方程,Rayleigh 數(Ra)和Prandtl 數(Pr)是RB 系統的重要控制參數.

RB 對流流場通過Zhang[30]所發展的基于有限差分的直接數值模擬得到.不失一般性,模擬方形區域的RB 對流,Pr數固定為1,計算了4 個不同的Ra數的工況.Ra數分別取106,107,108,109,對應的計算網格分別是128×128,512×512,720×720,1024×1024.當Ra數從106增大到109時,RB 對流系統會從層流狀態轉變為強非線性的湍流流動狀態,整個系統的速度及溫度脈動劇烈,會出現局部梯度較大的數據,對于系統控制方程的識別可能造成一定的困難.圖5 是t=80 時,不同Ra數的RB 對流系統的溫度云圖和速度矢量圖.

圖5 Pr=1,不同Ra 數的RB 對流系統的溫度云圖和速度矢量圖Fig.5 Temperature contour and velocity vector of RB convection for different Rayleigh number with Pr=1

在整個計算區域選取數據,采樣方法、導數計算方法以及Θ庫的構建與2.1 節相同.

表3 給出了Ra數為109的識別結果,PDE-FIND方法和LASSO 方法都正確識別出了Ra=109的RB對流系統中連續性方程和熱輸運方程,并且系數一致.而且PDE-FIND 方法仍然正確識別出了RB 對流的動量方程,對應Ra從106到109,識別的控制方程組平均系數誤差分別為0.36%,0.06%,0.40%,0.46%.而LASSO 方法對NS 方程仍沒有正確識別,在u方程識別中未選擇uxx.原因是由變量及變量各階導數組成的過于完備庫中,回歸系數趨于相近的項之間會產生分組效應(grouping effect)[31,32].而LASSO 方法對于具有分組效應的候選項的篩選可能會失敗,表現為在有分組效應的候選項中只選擇一項,而其他的候選項的系數會設為0.在Ra=106,Ra=107,Ra=108的RB 對流控制方程稀疏識別中,也存在同樣的現象.

表3 Ra=109,Pr=1RB 對流控制方程識別Table 3 Identification of governing equations for Rayleigh-B′enard convection with Ra=109,Pr=1

2.4 三維槽道湍流控制方程的識別

最后對三維槽道湍流進行了稀疏識別,數據通過離散統一的氣動格式(discrete unified gas-kinetic scheme,DUGKS)[33]模擬得到,計算域為4π×2×4π/3,槽道半高為1,計算網格為128×128×128,垂向采用非均勻網格,流向和展向采用均勻網格,壁面摩擦速度定義的雷諾數Reτ=180,此時的槽道流處于充分發展湍流狀態.圖6 給出了槽道展向中心截面的速度云圖,各種尺度的復雜結構清晰可見.在整個槽道中選取數據,采樣方法、導數計算方法以及Θ 庫的構建與2.1 節相同.

圖6 Reτ=180,三維槽道湍流展向中心截面速度云圖Fig.6 Velocity contour of 3D turbulent channel flow with Reτ=180 at z=0

表4 給出了三維槽道湍流控制方程的識別結果,PDE-FIND 方法對于正確地識別出了控制方程,并且平均系數誤差為0.37%,說明PDE-FIND 可以在在湍流狀態的流場數據中識別控制方程.而LASSO方法只識別出了連續性方程,這與前兩節的結果一致.說明在過完備候選庫中列相關性強弱對于PDEFIND 方法的稀疏識別沒有影響,PDE-FIND 方法不存在分組效應.本節還研究了不同y+的數據對識別的影響,取y+∈(0.21,19.3),y+∈(19.3,59.7),y+∈(59.7,182.9),都正確識別出了控制方程,系數平均誤差分別為0.42%,0.48%,0.39%.說明在三維槽道湍流控制方程的稀疏識別中,誤差與壁面距離幾乎無關.

表4 Reτ=180 三維槽道湍流控制方程識別結果Table 4 Identification of governing equations for 3D turbulent channel flow with Reτ=180

3 結論

本文采用PDE-FIND 方法及LASSO 方法對二維圓柱繞流、頂蓋驅動方腔流動、RB 對流和三維槽道湍流控制方程進行了稀疏識別,主要得到以下結論:

(1)使用PDE-FIND 方法識別出了流場中的控制方程,得到了正確的方程及精確的系數.

(2) LASSO 方法正確識別出了其中的線性偏微分控制方程如連續性方程,渦量輸運和溫度輸運方程,但對含非線性的偏微分方程識別效果差,這是由于LASSO 方法識別過程中侯選庫存在分組效應,降低了識別能力.而在PDE-FIND 方法中是不存在分組效應的,所以PDE-FIND 方法對強非線性方程識別有一定的優勢.

(3)豐富的流動結構對于數據驅動的流場控制方程的稀疏識別有重要作用,選擇具有豐富流動結構的區域的數據可以提升稀疏識別的有效性和精確性.

從2.3 節和2.4 節結果可知,PDE-FIND 方法對層流和湍流兩種流動狀態都可以精確的識別出控制方程.這對于復雜流動建模提供了一種新的研究方法,例如將PDE-FIND 方法用于湍流脈動量控制方程和湍流統計量的控制方程的系數識別.未來考慮將稀疏識別算法與k-ε 等湍流模式相結合,提高湍流模式在數值模擬中的精確性.另外,本文的流場數據都是通過精細的直接數值模擬得到的,PDE-FIND 都識別出了正確的方程形式,得到了準確的方程系數,但從實驗或者實測數據往往具有噪聲,對于含有噪聲的數據,PDE-FIND 方法識別能力往往變弱,如何提高PDE-FIND 方法在含有噪聲數據的識別能力是值得深入研究的課題.可行的方法之一是結合數據去噪技術和尋找最優的插值點數和多項式插值階數等,提高高階導數項的計算精度,從而提高稀疏識別能力.