陳辛陽，聶滋森，蔣子超，楊耿超，姚清河

中山大學航空航天學院，廣東廣州 510006

格子玻爾茲曼算法（LBM，Lattice Boltzmann method） 起源自格子氣自動機（LGA，lattice gas automaton），是一種新興的計算流體力學方法，相較于傳統計算流體力學（CFD， computational fluid dynamic） 方法，LBM 有易于處理邊界條件，適合復雜幾何形狀，程序結構簡單及可并行性高等優(yōu)點［1］。然而，LBM 存在計算資源消耗大以及對設備配置要求高等缺陷。具體來講，造成LBM 消耗大量計算資源的原因主要有兩點。第一，每一迭代內的計算資源消耗大，針對這一方面的相關研究主要是利用GPU 并行的方法降低計算資源消耗，目前已經形成了較為成熟的相關技術［2］；第二，由于LBM 是一個顯式算法，為了保障計算的穩(wěn)定性和精度，計算步長常常設置得較短，而較短的步長將導致迭代次數增多，大幅增加計算資源的消耗［3］。本文將從第二點出發(fā)，以單次模型計算代替較多計算步長小的普通LBM 迭代，從而減少計算資源的消耗。

在降低CFD 模型計算復雜度方面研究較多的是降維法（ROM，reduced order methods），比如：Xie 等的一種基于數據驅動的降維法（DDFROM），在構建非封閉ROM 后使用數據驅動模型矯正ROM，提高了模型的精度［4］；De Witt 和Lessig 等提出了改進ROM 法，提高計算效率并在小自由度下完成了可靠的流場可視化模擬［5］。ROM 方法在一些實際運用中也體現了其良好的性能與可觀的加速效果，Silva 等在空氣動力學計算中嘗試引入ROM 下的CFD 算法，取得了一定的成果［6］。然而，在早期的相關研究中大部分ROM 的結構較為簡單，影響了其性能。

隨著機器學習類方法的飛速發(fā)展，相關算法在各領域內應用廣泛，在CFD 領域中也是如此。利用ROM 類算法后，相關模型的復雜度得到了提升，模型性能也得到了提升。在熱對流問題的研究中，Yasin 等提出利用支持向量機（SVM， support vector machine） 對溫度分布和流場進行預測［7］。 而在針對粒子法的數據驅動研究中，Ladicky 等引入回歸森林（regression forests） 方法，近似模擬粒子行為模式，循環(huán)使用粒子預測模型完成對流場長時間運行的模擬［8］。

上述機器學習方法較為傳統，近年來的研究證實其效果不如包括深度學習算法在內的新型機器學習方法［9］。由Guo等提出的網絡模型完成從邊界條件直接到穩(wěn)態(tài)結果的預測過程［10］。這一直接預測穩(wěn)態(tài)結果的思路精度較高，但是對于非穩(wěn)態(tài)問題以及需要關注流場發(fā)展過程的穩(wěn)態(tài)問題并不適用。由Tompson 等運用卷積神經網絡（CNN，convolution neuron network） 取代壓強求解器，通過取代迭代過程降低了計算資源消耗［11］。然而，在顯式LBM 算法中并沒有一個迭代求解單個方程的子過程，也就導致這種思路不適合運用到LBM中。在對LBM 壓縮提速方面，Hennigh 提出了Lat-Net，利用深度學習模型對LBM 過程進行信息提取完成了一個神經網絡組成的流場模擬算法，在網絡中特別關注了對邊界條件信息的保持［12］。但是，由于網絡預測產生的誤差較大，在反復迭代中誤差累積較為明顯，導致結果有一定的偏差。同時對于邊界條件不鮮明的問題，該模型的效果有待驗證。

本文中采用深度學習模型對原有LBM 中的計算過程進行替換，以一個模型計算代替原有算法的多個迭代過程，達到減少迭代次數的目的，本文中將該模型稱為壓縮LBM （ C-LBM， compressed LBM） 模型。本文所采用的網絡應用領域為圖像生成領域，在相關領域內已有一些成熟的相關研究［13-15］。而，本文采用了一個較為新穎的網絡結構。具體來講，C-LBM 模型是以卷積長短期記憶網絡（ConvLSTM， convolutional LSTM network）［16］組建的編碼網絡和以殘差神經網絡（Resnet， residual network）［17］組建的解碼網絡為主體的深度神經網絡模型，通過輸入帶有時間和空間雙重維度的流場信息取得對跨越多個迭代后流場情況的預測結果。該模型的優(yōu)勢在于具備適應非穩(wěn)態(tài)問題和穩(wěn)態(tài)問題過程模擬的能力，同時具有較好的精度。

1 Lattice Bhatnagar-Gross-Krook（LBGK）模型

LBM 中主要研究的流場物理量是粒子速度分布函數，通過刻畫粒子速度分布函數隨時空的變化，完成對各種宏觀物理量變化情況的求解。粒子速度分布函數在統計力學中是一個十分重要的基本物理量，常用符號f(x，ε，t)表示。該物理量的定義為：在不考慮粒子的轉動而只考慮其平動的前提下，在某一時刻t，在空間中以x為中心的微元dx內，速度處于ε與ε+ dε之間的分子數為fdxdε。在LBM 中粒子分布函數被離散至各個離散速度方向并記為fi，本文使用的針對二維的D2Q9模型中共有9個離散速度方向，各個方向的速度記為ci，fi與ci一一對應［18］。由此，相應的宏觀密度和速度為

LBGK 模型是基于Bhatnagar-Gross-Krook（BGK） 碰撞模型算子建立的一種簡單LBM 模型，其離散形式控制方程為

該方程即為LBGK 方程， 其中τc為松弛因子。LBGK 方程包括對粒子遷移和碰撞兩方面的描述，在LBM中將這兩個部分分開計算，其表達式為

其中公式（3）描述粒子碰撞過程，而公式（4）表示遷移過程［19］。

通過以上關于LBGK 模型的描述不難發(fā)現，在給定空間大?。▁，y） 下，整個計算過程中的核心物理量，即離散粒子速度分布函數，實際上是一個大小不變的三維數組，其大小為（x，y，9），其中最后一維對應離散速度數，而二維卷積運算所需要的輸入規(guī)格剛好與這一數組相對應，即前兩維表示空間上的長和寬，最后一維為通道數剛好對應離散速度數。由此，本文的網絡模型中將原有離散粒子速度分布函數信息通過卷積的形式提取并完成相應計算，而卷積運算速度較快且對原有信息的保全度高，完成運算資源消耗少且能保證一定精度的算法。

2 基于深度學習的壓縮LBM模型

2.1 C-LBM模型網絡結構

本文中的網絡模型主要由ConvLSTM 和Resnet組成。C-LBM 模型的編碼部分由ConvLSTM 網絡組成，ConvLSTM 對于空間和時間的雙重處理能力可以將其視為最適合本文建立預測模型的工具［17］。解碼部分由Resnet 組成，由于其擁有抗退化的能力，可以降低本模型的訓練難度，使得模型的擬合效果更好［18］。圖1 為本文的網絡結構圖，其中輸入記為F(T)和F(T+100)，輸出記為F(T+1100)，該粒子分布函數F(T)（與原LBM 計算完成T次迭代的結果對應） 為

圖1 C-LBM模型網絡結構示意圖Fig. 1 C-LBM model network diagram

該網絡的輸入為相隔100個迭代、兩個時刻的分布函數。對于二維問題，兩個輸入均為三維矩陣，第一維對應空間長度，第二維對應空間寬度，第三維對應神經網絡中通道數。整個預測網絡的結構大致包括輸入預處理（標準化輸入），通過堆累ConvLSTM 層的網絡提取特征及擬合，Resnet 解碼，卷積層輸出以及輸出結果還原（逆標準化輸出） 五個過程，即從已有分布函數結果抽取時空信息，得到對未來分布函數的預測結果。預處理內容為標準化過程，由于在訓練中輸出也進行了標準化處理，因此在使用時需對輸出結果進行還原處理。具體來講，這一處理為神經網絡數據集標準化中的Z-Score標準化

其中μ為數據集均值，σ為數據集標準差，x為處理前數據而z為處理后數據。同時，ConvLSTM 的狀態(tài)（stateful） 被設置為真（true）， 批次大?。╞atch size） 被設置為1。對于本文中的stateful 為true的循環(huán)神經網絡，其網絡訓練的batch size與模型使用的batch size 大小是相同且固定的。為了構建單次只利用一批數據且內容為兩個輸入到一個輸出的模型，在訓練時必須將batch size設置為1。

2.2 C-LBM模型損失函數與優(yōu)化方法

損失函數（loss function） 在深度學習中有著極其重要的意義，本文針對LBM 的物理背景結合神經網絡訓練中常用的損失函數，設計了一個針對C-LBM 模型的損失函數，命名為CLBMLOSS，其表達式為

其中表示預處理還原后預測結果，fi表示還原后對應真實值，表示還原前預測結果，表示還原前對應真實值。本文中，系數α1=90，系數α2=10，目的在于使各單項誤差數量級接近，保證其均能發(fā)揮作用。

CLBMLOSS的第一部分為在神經網絡訓練中常用的二范數誤差MSE，即預測值與真實值相減結果平方的平均值。第二部分被命名為LBMLOSS，其內容為利用真實值與預測值經過碰撞過程的計算結果做差。神經網絡中輸入和輸出的分布函數均是計算到每個迭代的末尾時的結果，即完成了碰撞和遷移兩個過程的計算后得到的結果。由此，可以將這兩個輸出結果代入再完成一次碰撞過程，這樣得到兩個碰撞一次后的分布函數結果，如果真實值和預測值完全相同則該誤差為零。通過增加該部分誤差，帶入了LBM 物理背景的信息，包括流場速度輸入和雷諾數等。與此同時，這一項誤差的引入使得預測值盡可能與LBM 相適應。具體來講LLBMLOSS計算公式為

其中為由真實值fi計算得來的平衡態(tài)分布函數，而為由預測值計算得來的平衡態(tài)分布函數。式（7） 第三項誤差內容為將預測值的9 個通道求和，將其與真實值9 個通道求和結果相減，9 個通道和的物理意義是在各個格點處的密度，該誤差的意義在于引入一個宏觀物理量對模型進行調控。

本文中訓練采用反向傳播算法（back propagation），該方法由Hinton 等提出。簡單的講，其內容包括激勵傳播和權重更新兩個部分，先正向傳播得到輸出再反向傳播得到響應誤差，之后利用誤差獲得梯度并對權重進行更新，達到訓練網絡的目的［20］。本文中利用的優(yōu)化算法為自適應動量的隨機優(yōu)化方法（adaptive moment estimation）［21］。

3 算例實驗

本文中采用方腔環(huán)流模型完成對模型的訓練與測試，該模型被廣泛運用于計算流體力學算法的驗證中［22］。本文選擇方腔環(huán)流模型作為實驗模型的原因主要有以下幾點：首先，方腔環(huán)流問題的初始條件和邊界條件極為簡單，程序設計較為簡便；其次，方腔環(huán)流模型下的流場特點鮮明，便于結果之間進行對比；最后，方腔環(huán)流模型的流場情況隨參數變化較為明顯，特別是在不同雷諾數下流場變化較為明顯，有利于測試模型的泛化能力。值得注意的是，雖然方腔環(huán)流問題是一個穩(wěn)態(tài)問題，但是在本文的研究中著重關注整個流場演化過程，而并非僅僅關注計算的最終穩(wěn)態(tài)結果，由此證明該模型對流場發(fā)展過程的還原能力以及在非穩(wěn)態(tài)問題中的潛力。

3.1 模型訓練

表1為本文生成訓練集和測試集時LBM的相關基本物理參數的設定。

表1 生成流場的物理參數Table 1 Physical parameters of the flow field

表1 中L0、t0、v0分別為單位長度、單位時間和單位速度。本算例中的訓練集一共包含345組，測試集一共包括5組，測試集不參與訓練過程，每組訓練集和測試集中所包含的樣本批數相同，均為5批，一組樣本內打開stateful 關聯。由表達式表示整個訓練集，一組訓練集或測試集的輸入集合GIN和輸出集合GOUT為

其中j表示組數，j=0，1，2，…，350，Tj表示第j組內第一批輸入的第一個時間步對應的原LBM 計算中所需的迭代次數。

為驗證ConvLSTM 與Resnet 網絡結構以及損失函數CLBMLOSS 對模型性能的影響，本文訓練了3個使用不同網絡結構和訓練損失函數的模型，分別 以“MSE+ConvLSTM+Resnet”，“CLBMLOSS+ConvLSTM”和“CLBMLOSS+ConvLSTM+Resnet”進行標記，在3. 2 的誤差分析中也將沿用這三個模型進行相關測試。圖2中展示了這三個模型的訓練過程中對應損失函數值的下降過程，橫坐標為訓練輪數。

圖2 訓練中損失函數值下降情況Fig. 2 Loss function value decline during training

圖2中，使用ConvLSTM+Resnet 結構后損失函數最終收斂數值明顯較低，初步證明了這種結構的模型具有更好的學習能力和抗退化能力。相較于MSE 訓練的模型，CLBMLOSS 訓練的模型損失函數最終收斂數值略低，表明引入CLBMLOSS 對模型的學習能力沒有產生嚴重影響。同時CLBMLOSS 包括MSE 部分，同一結果下前者必然大于后者，而圖2中在模型收斂后使用CLBMLOSS 訓練的模型損失函數值更低。由此可知，若使用相同評價標準，利用CLBMLOSS訓練的模型精度更高。

3.2 誤差分析

誤差分析部分首先驗證模型在Re=2 000 的測試集上的精度。本文模型的重要創(chuàng)新點的第一部分是結合ConvLSTM 與Resnet 的網絡結構，第二部分是針對LBM 物理背景設計的損失函數CLBMLOSS，以下也將檢驗這些創(chuàng)新點對模型精度帶來的增益效果。

Re=2 000 的實驗集共6 組，每組5 批樣本，總共30 批樣本，將這些樣本依次編號為1 到30。由此1 到5 為一組，6 到10 為一組，依此類推，各組之間的狀態(tài)（state） 不能傳遞，即每隔5 批就將state 清零。以下利用均方差總體評價模型輸出精度情況，令函數F (x，y，i) 表示位置為(x，y) 處，全局均方差計算公式為

圖3為各個不同網絡結構和優(yōu)化方法下的全局均方差值，橫軸為各批次預測輸出對應普通LBM結果的迭代次數，即第一組樣本為1 800、1 900、2 000、2 100和2 200次迭代結果，依此類推。

圖3 Re=2 000全局均方差誤差Fig. 3 Global mean square error when Re=2 000

圖3中，添加了Resnet 解碼器的網絡在誤差上明顯低于未添加解碼器的網絡，表明了Resnet 解碼器對模型精度提升的貢獻。同時，利用CLBMLOSS對模型預測精度產生了一定幅度的積極影響，這也證明了使用CLBMLOSS 在提升預測準確性上的意義。從圖3中不難看出，每組的第一批樣本的預測精度弱于同組內其他批次的，這是由于第一批預測模型中還沒有狀態(tài)量（state） 的數據，導致預測的根據不足，從第二批預測開始模型同時使用本批數據和state 量，精度較高，這也提示在使用模型時最好用第二批及以后的數據作為最終輸出結果。

下面，通過速度場標量圖4比對真實結果，并展示來自真實值與預測值之間的殘差圖，同時附加對居中行和居中列結果的比較。以下誤差分析均在實驗集5 上（對應普通LBM 算法完成2 200 次迭代結果） 進行，輸出結果對應普通LBGK 算法完成2 000 次迭代。其中，圖4 （a1） 與4 （a2） 為輸入數據所求得的速率標量場；圖4 （b） 為輸出結果對應的常規(guī)LBM 速率標量場即真實值；圖4（c1）、4 （d1） 和4 （e1） 為各模型預測的速率標量場；圖4 （c2）、4 （d2） 和4 （e2） 為各模型預測值與真實值之間的絕對誤差。

圖4 Re=2 000時全流場的可視化結果與誤差分析Fig. 4 Visualization results and error analysis of the flow field when Re=2 000

從圖4 中不難看出，同時利用了ConvLSTM+Resnet結構和CLBMLOSS優(yōu)化函數的模型有更高的精度。 首先， 從對比中可以看出， 具有Conv-LSTM+Resnet 結構的模型最大誤差和整體誤差都明顯小于沒有這一結構的模型，這體現了Resnet 提高了模型復雜度和學習能力；其次，具有CLBMLOSS 優(yōu)化函數的模型零誤差及極低誤差區(qū)域面積更大，這是由于該損失函數的引入增加了物理信息，引導模型加強識別流場中的低速部分，同時提高各個通道之間的關聯性，從而有效地降低了誤差，尤其是低速區(qū)域的誤差。為了進一步驗證，以下對比在居中行與居中列處，各個模型預測的速率結果和真實速率結果，同時加入絕對誤差圖像，即真實值與預測值的絕對差值。圖5（a1）與5（a2）展示了結果曲線，圖5（b1）與5（b2）則是誤差曲線圖。

圖5 Re=2 000時居中行與列可視化結果與誤差分析Fig. 5 Middle row and column visualization results and error analysis when Re=2 000

圖5進一步證明了從二范數誤差圖中所得出結論的正確性。加入Resnet 降低了最大誤差，而CLBMLOSS 的引入對低速區(qū)域誤差降低有較好作用，在高速區(qū)有無CLBMLOSS 對精度的影響并不明顯，不過在高速區(qū)所有模型的誤差相較于真實值都很小，影響并不明顯?？偠灾?，該實驗體現了ConvLSTM 層和Resnet 層組合網絡以及CLBMLOSS損失函數的強大功效。

為了進一步驗證本文中各模型在泛化問題中的表現，本文選用了Re=4 000 的非訓練雷諾數測試集，對該測試集下的精度進行了校驗。同樣的，實驗集共6組，每組5批樣本，總共30批樣本，將這些樣本依次編號為1 到30。圖6 是不同網絡結構和優(yōu)化方法下的全局均方差值，橫軸為各批次預測輸出對應普通LBM 結果的迭代次數，即第一組樣本對應1 900、2 000、2 100、2 200 和2 300 次迭代結果，依此類推。

從圖6 中可以看出，CLBMLOSS 的引入在提升模型泛化能力中的效果較為明顯，引入CLBMLOSS損失函數的兩個模型預測精度更高，這是由于加入CLBMLOSS 后在優(yōu)化過程中引入了雷諾數信息，有效增強了模型物理泛用性。在這一部分實驗中，加入Resnet 的網絡只在約半數的批次中有更高精度。在這一點上，我們認為可能網絡結構復雜化增大了過擬合的程度，不過基于上一個實驗中Resnet的增強效果，我們依然認可引入Resnet的提升意義。

圖6 Re=4 000全局均方差誤差Fig. 6 Global mean square error when Re=4 000

同樣的，圖7為預測結果精度和真實值與預測值之間的殘差，同時附加對居中行和居中列結果的比較。以下誤差分析均在實驗集25 上進行（對應普通LBM 算法完成4 200個迭代）。同樣的，圖7（a1） 與7 （a2） 為輸入數據所求得的速率標量場，圖7 （b） 為輸出結果所對應的常規(guī)LBM 速率標量場即真實值，圖7 （c1）、7 （d1） 和7 （e1） 為各模型預測的速率標量場，圖7 （c2）、7 （d2） 和7（e2） 為各模型預測值與真實值之間絕對誤差。為了進一步驗證，同樣對比在居中行與居中列處各個模型的預測速率和真實速率，同時加入絕對誤差圖像，即真實值與預測值的絕對差值。圖8（a1） 與8 （a2） 展示了結果曲線，圖8 （b1） 與8（b2） 則是誤差曲線圖。

圖8 Re=4 000居中行與列可視化結果與誤差分析結果Fig. 8 Middle row and column visualization results and error analysis when Re=4 000

從圖7中可以明顯看出CLBMLOSS 在模型泛化能力方面的提升，總體誤差和最大誤差都明顯下降，同時在上個實驗中提到的對低速區(qū)精度的優(yōu)化效果也依然存在。

圖7 Re=4 000時全流場可視化結果與誤差分析結果Fig.7 Visualization results and error analysis of the flow field when Re=4 000

從圖7 和圖8 中的誤差數據可以看出，Conv-LSTM 和Resnet 組成的網絡結構以及CLBMLOSS 損失函數對網絡精度的提升有積極意義；同時，本模型具備一定程度的多雷諾數泛化能力，在不利用Re=4 000 訓練集訓練的前提下，仍然可以應用于相關算例。

3.3 計算壓縮效果分析

本文測試平臺為Intel （R） Core （TM） i7-9750H CPU@ 2. 60Hz 以及NVIDIA GeForce GTX 1660 Ti。這一平臺屬于典型的家用機配置，該測試環(huán)境有利于模擬本模型算法在受限制的計算機能下的計算能力，這也是本算法的設計目標應用領域，即針對計算能力有限的設備完成具備嚴謹物理背景的流體仿真。表2中展示了不同雷諾數下不同時段利用模型完成計算和原有算法計算的耗時對比，時段由完成的迭代數量表示。

表2 耗時對比Table 2 Time-consuming comparison

從表2中不難看出，計算耗時與計算區(qū)間和雷諾數并無關聯，原串行算法耗時約為模型計算的15倍。在評價LBM 的計算效率時引入一個專門量，即每秒百萬格子更新數（MLUPs， million lattice updates per second），其計算公式為

其中nlattice為格子單元的數量，compute time 為計算時間。本文中的算例下nlattice=25 600，表3 列出了C-LBM 模型下連續(xù)預測計算過程的各項耗時及由此計算得到的MLUPs。對于C-LBM 模型由于其計算一次相當于普通LBM 計算1 000 次，所以在計算其對應的MLUPs時還需要考慮到這一對應關系。

表3 計算效率實驗結果Table 3 Calculation efficiency experiment results

在家用級別平臺上，本文新模型的計算效率為串行LBM 程序的15 倍左右。由此可見，相較于串行LBM 程序，本文提出的新模型在減小計算資源消耗方面效果顯著，該算法的實驗結果達到預期。

4 結 論

為提高傳統LBM 的計算效率，本文實現了對LBM 原有計算過程的簡化壓縮，克服LBM 在消耗計算資源方面的劣勢?；谝陨夏繕?，本文利用ConvLSTM 層和Resnet 層等深度學習網絡結構設計了一個預測網絡模型并命名為C-LBM 模型，該模型通過輸入早期的粒子分布函數數據對跨越多個迭代后的粒子分布函數進行預測。完成模型訓練后，本文對該模型預測精度以及計算效率進行了分析與驗證，得出了該模型具備基本可靠性和一定泛化能力的結論，并達到了大幅壓縮計算資源消耗的目標。