基于數(shù)據(jù)驅(qū)動的流場控制方程的稀疏識別1)

江 昊 王伯福 盧志明

(上海大學(xué),力學(xué)與工程科學(xué)學(xué)院,上海市應(yīng)用數(shù)學(xué)和力學(xué)研究所,上海 200072)

引言

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

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

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

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

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

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

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

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

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

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

如果Θ內(nèi)所有的項都屬于該流場系統(tǒng)的控制方程,此時可以使用最小二乘法確定系數(shù)矩陣

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

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

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

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

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

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

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

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

2 識別結(jié)果

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

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

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

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

2.2 頂蓋驅(qū)動方腔流控制方程的識別

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

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

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

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

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

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

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

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

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

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

RB 對流流場通過Zhang[30]所發(fā)展的基于有限差分的直接數(shù)值模擬得到.不失一般性,模擬方形區(qū)域的RB 對流,Pr數(shù)固定為1,計算了4 個不同的Ra數(shù)的工況.Ra數(shù)分別取106,107,108,109,對應(yīng)的計算網(wǎng)格分別是128×128,512×512,720×720,1024×1024.當(dāng)Ra數(shù)從106增大到109時,RB 對流系統(tǒng)會從層流狀態(tài)轉(zhuǎn)變?yōu)閺姺蔷€性的湍流流動狀態(tài),整個系統(tǒng)的速度及溫度脈動劇烈,會出現(xiàn)局部梯度較大的數(shù)據(jù),對于系統(tǒng)控制方程的識別可能造成一定的困難.圖5 是t=80 時,不同Ra數(shù)的RB 對流系統(tǒng)的溫度云圖和速度矢量圖.

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

在整個計算區(qū)域選取數(shù)據(jù),采樣方法、導(dǎo)數(shù)計算方法以及Θ庫的構(gòu)建與2.1 節(jié)相同.

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

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

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