李治龍, 周恒安, 翁俊輝, 歐洺余, 朱宏娜*

(1.西南交通大學(xué)物理科學(xué)與技術(shù)學(xué)院, 成都 610031; 2.西南交通大學(xué)信息科學(xué)與技術(shù)學(xué)院, 成都 610031)

非定常周期性流體流動(dòng)的研究是流體力學(xué)中熱門研究領(lǐng)域,其物理模型通常由一個(gè)或多個(gè)偏微分方程描述。常用的求解方法是通過離散偏微分控制方程建立數(shù)值求解方法,形成一個(gè)可求解的數(shù)學(xué)表達(dá)式[1-2]。基于納維-斯托克斯(Navier-Stokes, N-S)的計(jì)算流體力學(xué)(computational fluid dynamics, CFD)數(shù)值模擬是研究非定常周期性流體流動(dòng)的常用方法之一,但CFD對(duì)網(wǎng)格要求高,求解N-S方程的過程復(fù)雜,以致占用較多計(jì)算時(shí)間。特別在網(wǎng)格數(shù)量大且?guī)缀瓮庑螐?fù)雜等情況下,CFD方法會(huì)占用大量運(yùn)行時(shí)間[3]。

近年來(lái),深度學(xué)習(xí)技術(shù)發(fā)展迅速,廣泛運(yùn)用于計(jì)算機(jī)視覺、自然語(yǔ)言處理、語(yǔ)音識(shí)別等領(lǐng)域[4-5]。非定常周期性流體流動(dòng)問題具有非線性強(qiáng)且數(shù)據(jù)量大等特點(diǎn)使得深度學(xué)習(xí)在此領(lǐng)域具有重大應(yīng)用前景[6]?！叭f(wàn)能近似定理”[7-8]表明:由于神經(jīng)網(wǎng)絡(luò)在理論上可以近似任意函數(shù),所以深度學(xué)習(xí)在理論上也具備預(yù)測(cè)流體流動(dòng)的可行性。

基于深度學(xué)習(xí)在流體流動(dòng)預(yù)測(cè)領(lǐng)域中的應(yīng)用,也開展了相關(guān)研究。Guo等[9]以玻爾茲曼算法為訓(xùn)練數(shù)據(jù)集,對(duì)卷積神經(jīng)網(wǎng)絡(luò)(convolutional neural network, CNN)進(jìn)行訓(xùn)練,成功模擬了不同形狀結(jié)構(gòu)物的穩(wěn)態(tài)流場(chǎng)結(jié)果。Parish等[10]通過反演模型推斷方法得到平均流物理量和模型的映射關(guān)系,使用機(jī)器學(xué)習(xí)方法,通過仿真和實(shí)驗(yàn)數(shù)據(jù)對(duì)模型進(jìn)行修正,得到與數(shù)值模型相近的結(jié)果。Han等[11]通過深度神經(jīng)網(wǎng)絡(luò)從高維非定常流場(chǎng)數(shù)據(jù)中提取時(shí)空特征,能利用歷史數(shù)據(jù)預(yù)測(cè)未來(lái)時(shí)刻的數(shù)據(jù)。王臣等[12]運(yùn)用CNN對(duì)流體結(jié)構(gòu)相互作用問題進(jìn)行求解,求解了不同流速、不同結(jié)構(gòu)物形態(tài)及位置等條件下的流體結(jié)構(gòu)相互作用問題,實(shí)現(xiàn)結(jié)果的快速預(yù)測(cè)?；菪挠甑萚13]將條件生成對(duì)抗網(wǎng)絡(luò)與卷積神經(jīng)網(wǎng)絡(luò)結(jié)合,預(yù)測(cè)了不同時(shí)刻的流場(chǎng)變量,實(shí)現(xiàn)了結(jié)果的快速預(yù)測(cè)。Li等[14]在長(zhǎng)短期記憶網(wǎng)絡(luò)(long short-term memory, LSTM)中引入數(shù)值方法,將流場(chǎng)的時(shí)空數(shù)據(jù)輸入網(wǎng)絡(luò),構(gòu)成閉環(huán)序列結(jié)構(gòu),降低了測(cè)量和模型參數(shù)在長(zhǎng)期預(yù)測(cè)中的誤差累積,提高了模型的泛化性。

雖然傳統(tǒng)的深度學(xué)習(xí)方法也能對(duì)流體流動(dòng)進(jìn)行預(yù)測(cè),但由于其“獨(dú)立同分布”的模型假設(shè)和只能處理歐幾里得數(shù)據(jù)的缺陷導(dǎo)致預(yù)測(cè)精度不高。而圖卷積循環(huán)網(wǎng)絡(luò)是一種基于圖結(jié)構(gòu)的深度學(xué)習(xí)方法,通過學(xué)習(xí)圖結(jié)構(gòu)數(shù)據(jù)特征,能有效處理流體流動(dòng)數(shù)據(jù)。流體流動(dòng)數(shù)據(jù)作為典型的時(shí)空數(shù)據(jù),具有時(shí)空依賴性。對(duì)流體流動(dòng)數(shù)據(jù)進(jìn)行建模,需要考慮未來(lái)時(shí)刻的流體流動(dòng)數(shù)據(jù)對(duì)歷史的流體流動(dòng)數(shù)據(jù)時(shí)空依賴。若在建立模型時(shí)不考慮歷史的流體流動(dòng)數(shù)據(jù),會(huì)丟失流體流動(dòng)的時(shí)空連續(xù)性關(guān)系,導(dǎo)致預(yù)測(cè)結(jié)果出現(xiàn)較大偏差。因此,本文采用圖卷積循環(huán)網(wǎng)絡(luò)對(duì)流體流動(dòng)進(jìn)行預(yù)測(cè),建模歷史數(shù)據(jù)以及未來(lái)數(shù)據(jù)間的依賴關(guān)系,避免求解復(fù)雜物理方程,實(shí)現(xiàn)流場(chǎng)結(jié)果的快速預(yù)測(cè)。相比于復(fù)雜物理方程求解,能達(dá)到較高的精度并明顯縮短計(jì)算耗時(shí),為快速實(shí)時(shí)預(yù)測(cè)提供支撐。

1 圖卷積循環(huán)網(wǎng)絡(luò)

1.1 流體流動(dòng)預(yù)測(cè)任務(wù)

{y1,y2,…,yh}=model(X1,X2,…,Xt)

(1)

1.2 圖卷積神經(jīng)網(wǎng)絡(luò)

圖卷積神經(jīng)網(wǎng)絡(luò)(graph convolution neural network, GCN)是一種多層神經(jīng)網(wǎng)絡(luò)[15]。GCN可直接在圖上作用,屬于譜域圖卷積,其核心思想是中心節(jié)點(diǎn)對(duì)其鄰居節(jié)點(diǎn)進(jìn)行信息聚合,并將鄰居節(jié)點(diǎn)的特征信息進(jìn)行編碼[16]。譜域圖卷積通過傅里葉變換將時(shí)域信號(hào)轉(zhuǎn)換到頻域,過濾后將結(jié)果通過傅里葉反變換變換到時(shí)域。在譜域圖卷積中,圖由對(duì)應(yīng)的拉普拉斯矩陣L表示為

L=D-A

(2)

式(2)中:D為圖中節(jié)點(diǎn)的度矩陣;A為圖中節(jié)點(diǎn)的鄰接矩陣。通常由正則化的拉普拉斯矩陣Lsym表示為

(3)

式(3)中:I為單位矩陣。將L進(jìn)行特征值分解,得

L=UΛUT

(4)

式(4)中:Λ為對(duì)角矩陣L;對(duì)角線的值為L(zhǎng)的特征值;U為L(zhǎng)的特征向量矩陣。

GCN利用定義在頻域中對(duì)角化的線性算子,來(lái)代替經(jīng)典卷積算子實(shí)現(xiàn)卷積操作[17]。GCN通過傅里葉變換,將節(jié)點(diǎn)特征映射到頻域,經(jīng)過濾操作后,再通過傅里葉逆變換獲得時(shí)域節(jié)點(diǎn)特征,從而得到譜域圖卷積[18]。將卷積核gθ參數(shù)化為切比雪夫多項(xiàng)式TK的K-1階展開,表示為

(5)

式(5)中:θ∈RK為切比雪夫系數(shù)的向量。使用卷積核gθ對(duì)圖G進(jìn)行卷積(*),即

gθ*x=UgθUTx

(6)

譜域圖卷積采用的層間傳播通式為

Hl+1=f(Hl,A)

(7)

式(7)中:l為圖卷積的層數(shù);Hl為第l層圖卷積的激活值,初始的輸入表示為H0=x∈RN×n,N為圖中的節(jié)點(diǎn)數(shù),n為節(jié)點(diǎn)的特征數(shù),基本實(shí)現(xiàn)為

(8)

1.3 門控循環(huán)單元

門控循環(huán)單元(gate recurrent unit, GRU)是LSTM的變體[19]。GRU旨在解決循環(huán)神經(jīng)網(wǎng)絡(luò)(recurrent neural network, RNN)中梯度消失的問題,并能提高訓(xùn)練效率[19]。

GRU使用門控機(jī)制控制輸入和記憶信息,并在當(dāng)前時(shí)間步輸出預(yù)測(cè)信息。GRU有兩個(gè)門,分別為一個(gè)重置門和一個(gè)更新門,其結(jié)構(gòu)如圖1所示。重置門定義了前面單元結(jié)構(gòu)傳入的記憶信息與本單元輸入信息相結(jié)合的方式,更新門定義了前面單元結(jié)構(gòu)傳入本單元信息保留的量。

圖1 GRU結(jié)構(gòu)圖Fig.1 The structure of GRU

網(wǎng)絡(luò)的前向傳播公式為

(9)

GRU相互連接,形成循環(huán)結(jié)構(gòu),當(dāng)前單元結(jié)構(gòu)的輸入是由之前單元結(jié)構(gòu)的最終記憶信息和當(dāng)前單元結(jié)構(gòu)的輸入信息的融合結(jié)果。將具有時(shí)間關(guān)聯(lián)性的數(shù)據(jù)按時(shí)間順序輸入到GRU串聯(lián)組成的網(wǎng)絡(luò)中,經(jīng)過不斷學(xué)習(xí),就能捕捉到輸入時(shí)間序列數(shù)據(jù)中的時(shí)間相關(guān)性。

1.4 圖卷積循環(huán)網(wǎng)絡(luò)

考慮到流體流動(dòng)具有時(shí)空特征,本文中設(shè)計(jì)基于圖卷積循環(huán)網(wǎng)絡(luò)流體流動(dòng)預(yù)測(cè)模型來(lái)完成非定常周期性流體流動(dòng)的預(yù)測(cè)任務(wù),模型框架如圖2所示。模型的輸入為T個(gè)時(shí)間步的歷史數(shù)據(jù),通過GCN捕捉空間相關(guān)性,得到流體節(jié)點(diǎn)網(wǎng)絡(luò)的固定拓?fù)浣Y(jié)構(gòu)和學(xué)習(xí)序列之間的自適應(yīng)關(guān)聯(lián)關(guān)系,然后將具有空間特征的時(shí)間序列輸入到GRU中,經(jīng)GRU之后信息傳遞獲取序列數(shù)據(jù)的變化動(dòng)態(tài),使用ReLU激活函數(shù),最終得到未來(lái)時(shí)刻的預(yù)測(cè)結(jié)果。其計(jì)算過程為

圖2 模型框架Fig.2 Model architecture diagram

(10)

總的來(lái)說(shuō),圖卷積循環(huán)網(wǎng)絡(luò)模型既具有GCN處理圖數(shù)據(jù)的能力,又具有GRU捕捉時(shí)間序列的能力,能處理流體流動(dòng)數(shù)據(jù)中復(fù)雜的時(shí)空依賴關(guān)系。一方面,GCN被用來(lái)捕獲流體數(shù)據(jù)的拓?fù)浣Y(jié)構(gòu)中節(jié)點(diǎn)之間的空間依賴關(guān)系。另一方面,GRU被用來(lái)捕獲圖中邊信息的動(dòng)態(tài)變化以得到節(jié)點(diǎn)自身的時(shí)間依賴性。通過以上模塊的有效結(jié)合,充分建立圖神卷積神經(jīng)網(wǎng)絡(luò)對(duì)流體流動(dòng)時(shí)空依賴性的模型,完成對(duì)流場(chǎng)參數(shù)的預(yù)測(cè)。

2 實(shí)驗(yàn)及結(jié)果

2.1 數(shù)據(jù)準(zhǔn)備

圓柱繞流是CFD中常見的非定常周期性流體流動(dòng)案例,因此本文中采用COMSOL多物理場(chǎng)仿真軟件對(duì)流域長(zhǎng)度L=0.5 m、流域?qū)挾萕=0.2 m的圓柱繞流場(chǎng)進(jìn)行了數(shù)值模擬,并通過改變圓柱的半徑、位置以及入口速度得到500個(gè)圓柱繞流的流場(chǎng)變量數(shù)據(jù)集。計(jì)算網(wǎng)格如圖3所示,每個(gè)網(wǎng)格平均有1 852個(gè)節(jié)點(diǎn),2 972條邊。

圖3 計(jì)算網(wǎng)格Fig.3 Sketch of computation grid

模型輸出的時(shí)間為3 s,間隔步長(zhǎng)為0.01 s,共300個(gè)時(shí)間步,流場(chǎng)變量為3個(gè),分別為流體流動(dòng)的流速、渦旋以及壓力,按照8∶2的比例將數(shù)據(jù)集劃分成訓(xùn)練集和測(cè)試集。

2.2 實(shí)驗(yàn)設(shè)置

在實(shí)驗(yàn)中,采用Python語(yǔ)言編寫代碼,使用Pytorch深度學(xué)習(xí)框架實(shí)現(xiàn),用前200個(gè)時(shí)間步長(zhǎng)的流場(chǎng)變量數(shù)據(jù)分別預(yù)測(cè)了未來(lái)1、50、100個(gè)時(shí)間步長(zhǎng)的流場(chǎng)變量。

本文中主要涉及的超參數(shù)包括批數(shù)據(jù)大小、衰減率、學(xué)習(xí)率和迭代次數(shù)。網(wǎng)絡(luò)設(shè)置的批數(shù)據(jù)大小為16,衰減率為0.001,學(xué)習(xí)率為0.01,迭代次數(shù)為2 000輪,圖卷積為1層,隱藏維度為32,輸出維度16,GRU為2層,隱藏層單元數(shù)為64,選擇MSE loss作為損失函數(shù),Adam為優(yōu)化器,ReLU函數(shù)為激活函數(shù)。所有實(shí)驗(yàn)均在具有16 GB RAM的單個(gè)NVIDIA 2080Ti GPU上進(jìn)行。

2.3 評(píng)價(jià)指標(biāo)

為了驗(yàn)證圖卷積循環(huán)網(wǎng)絡(luò)的預(yù)測(cè)性能,選取平均絕對(duì)百分比誤差(mean absolute percentage error, MAPE)和平均絕對(duì)誤差(mean absolute error, MAE)作為評(píng)價(jià)本文模型預(yù)測(cè)精度的指標(biāo)。MAPE通過評(píng)價(jià)預(yù)測(cè)值和真實(shí)值的偏離程度,從而評(píng)價(jià)模型的精度。MAE通過平均所有預(yù)測(cè)值與真實(shí)值的偏差的絕對(duì)值,避免誤差相互抵消,從而評(píng)價(jià)模型的精度。MAPE和MAE的計(jì)算公式為

(11)

(12)

式中:N為評(píng)價(jià)模型的樣本數(shù)量;yi為第i個(gè)樣本的真實(shí)值;ypred為第i個(gè)樣本的預(yù)測(cè)值。

2.4 結(jié)果討論

將COMSOL軟件仿真的數(shù)據(jù)作為真實(shí)樣本,將圖卷積循環(huán)網(wǎng)絡(luò)預(yù)測(cè)的1、50、100個(gè)時(shí)間步長(zhǎng)的流體流動(dòng)流速、渦旋以及壓力等流場(chǎng)變量數(shù)據(jù)進(jìn)行可視化對(duì)比展示,驗(yàn)證預(yù)測(cè)模型的準(zhǔn)確性,對(duì)比結(jié)果如圖4～圖6所示。

圖4 速度預(yù)測(cè)結(jié)果對(duì)比Fig.4 Comparison of speed prediction results

圖5 渦旋預(yù)測(cè)結(jié)果對(duì)比Fig.5 Comparison of vortex prediction results

圖6 壓力預(yù)測(cè)結(jié)果對(duì)比Fig.6 Comparison of pressure prediction results

圖4～圖6中,將真實(shí)樣本與圖卷積循環(huán)網(wǎng)絡(luò)對(duì)流場(chǎng)變量預(yù)測(cè)結(jié)果進(jìn)行比較,可見圖卷積循環(huán)網(wǎng)絡(luò)預(yù)測(cè)效果較為準(zhǔn)確,隨著預(yù)測(cè)時(shí)間的增加,部分小尺度結(jié)構(gòu)的偏差逐漸增大,但總體獲得了較好的預(yù)測(cè)結(jié)果。從圖4～圖6中可以看出,網(wǎng)絡(luò)捕捉到三個(gè)流體特征的絕大部分規(guī)律,較為準(zhǔn)確地預(yù)測(cè)出未來(lái)時(shí)刻三個(gè)流場(chǎng)變量特征。從預(yù)測(cè)的未來(lái)不同時(shí)刻的結(jié)果來(lái)看,預(yù)測(cè)1個(gè)時(shí)間步τ的預(yù)測(cè)結(jié)果最好,50個(gè)時(shí)間步的預(yù)測(cè)效果較為準(zhǔn)確,100個(gè)時(shí)間步的預(yù)測(cè)結(jié)果可以觀察到等值線的抖動(dòng),這是由于預(yù)測(cè)結(jié)果與真實(shí)樣本之間的預(yù)測(cè)誤差逐步累積,從而與真實(shí)樣本出現(xiàn)一定的偏差,導(dǎo)致預(yù)測(cè)結(jié)果出現(xiàn)誤差,并隨著預(yù)測(cè)時(shí)間的增加,流場(chǎng)變量的預(yù)測(cè)結(jié)果出現(xiàn)變形。從圖4～圖6中還可以看出,圖卷積循環(huán)網(wǎng)絡(luò)能夠較好地預(yù)測(cè)流體流動(dòng)的不同流場(chǎng)變量,捕捉流體流動(dòng)的特征,經(jīng)可視化后,較好地還原了流線及渦的形態(tài),符合流體運(yùn)動(dòng)物理規(guī)律。以上實(shí)驗(yàn)表明,圖卷積循環(huán)網(wǎng)絡(luò)具有模擬復(fù)雜流體流動(dòng)的潛力。

表1給出了本模型在測(cè)試集上預(yù)測(cè)1、50、100個(gè)時(shí)間步的流場(chǎng)變量預(yù)測(cè)的精度。

表1 不同時(shí)間步長(zhǎng)的預(yù)測(cè)精度

從表1中可見,對(duì)于MAPE,三個(gè)流場(chǎng)變量的預(yù)測(cè)誤差都隨著預(yù)測(cè)時(shí)間步的增加而增大,隨著預(yù)測(cè)時(shí)間步的增加,誤差逐步累積;從數(shù)值上看,預(yù)測(cè)1個(gè)時(shí)間步的MAPE最小,預(yù)測(cè)50個(gè)時(shí)間步的MAPE也較小,預(yù)測(cè)100個(gè)時(shí)間步的MAPE最大,三個(gè)預(yù)測(cè)時(shí)間步的MAPE都在15%以下,說(shuō)明預(yù)測(cè)的精度較高。對(duì)于MAE,三個(gè)流場(chǎng)變量的預(yù)測(cè)誤差同樣隨著預(yù)測(cè)時(shí)間步的增加而增大,但不同流場(chǎng)變量的MAE不同,這是由于不同流場(chǎng)變量的預(yù)測(cè)數(shù)值范圍不同,所以導(dǎo)致不同流場(chǎng)變量的MAE差別較大?？傮w上看,三個(gè)流場(chǎng)變量的預(yù)測(cè)誤差都較小,具有一定的模型精度。在定量對(duì)比中可知,本文模型預(yù)測(cè)結(jié)果誤差較小,說(shuō)明本文的預(yù)測(cè)方法效果良好,能夠較好的預(yù)測(cè)未來(lái)時(shí)刻的流場(chǎng)變量。

在計(jì)算耗時(shí)方面,本文中統(tǒng)計(jì)了CFD仿真與圖卷積循環(huán)網(wǎng)絡(luò)預(yù)測(cè)第1、50和100時(shí)間步的平均時(shí)間消耗,如表2所示。結(jié)果表明,網(wǎng)絡(luò)訓(xùn)練完成后,相比于CFD數(shù)值模擬,網(wǎng)絡(luò)預(yù)測(cè)流體流場(chǎng)變量的計(jì)算時(shí)間大大減少,縮短了5倍以上,大大降低了計(jì)算時(shí)間成本,縮短計(jì)算周期,可實(shí)時(shí)預(yù)測(cè)未來(lái)時(shí)刻的流體流場(chǎng)變量。

表2 仿真時(shí)間消耗

3 結(jié)論

基于圖卷積循環(huán)網(wǎng)絡(luò)對(duì)非定常周期性流體流動(dòng)狀態(tài)進(jìn)行預(yù)測(cè),通過神經(jīng)網(wǎng)絡(luò)的訓(xùn)練避免求解復(fù)雜物理方程,實(shí)現(xiàn)流體流場(chǎng)變量的快速預(yù)測(cè)。通過分析歷史數(shù)據(jù),該模型能夠?qū)崿F(xiàn)流體流動(dòng)狀態(tài)的實(shí)時(shí)預(yù)測(cè),預(yù)測(cè)結(jié)果表明,在較低預(yù)測(cè)誤差下預(yù)測(cè)時(shí)間步可達(dá)到100。針對(duì)圓柱繞流非定常周期性流動(dòng)問題,通過模擬不同流速條件、不同圓柱位置條件等情況數(shù)據(jù),提高圖卷積循環(huán)網(wǎng)絡(luò)的泛化性能,使其學(xué)習(xí)并擬合流體流動(dòng)規(guī)律,在面對(duì)新的輸入條件時(shí)也能預(yù)測(cè)出較為準(zhǔn)確的流場(chǎng)變量,避免求解復(fù)雜的物理方程。相比于CFD仿真模擬,本文預(yù)測(cè)模型計(jì)算耗時(shí)較短,具有明顯的優(yōu)勢(shì)。在數(shù)據(jù)集的選取方面,本文中選取的數(shù)據(jù)集具有一定局限性,只針對(duì)了圓柱繞流問題進(jìn)行了討論,沒有討論其他非定常流場(chǎng)的情況,未來(lái)將對(duì)其他非定常流場(chǎng)進(jìn)行研究,并將二維流場(chǎng)擴(kuò)展到三維流場(chǎng)。