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

中山大學(xué)航空航天學(xué)院，廣東 廣州 510006

格子玻爾茲曼方法（LBM）［1］于20 世紀80 年代后期發(fā)展起來，為一種介觀尺度下的流體模擬方法，與傳統(tǒng)的通過求解納維斯托克斯（N-S）方程來對流體進行模擬的數(shù)值方法相比較，LBM 物理概念更清晰，它易于并行，且擅于處理復(fù)雜邊界，受到了很多科研工作者的關(guān)注，在流體力學(xué)，化學(xué)反應(yīng)，量子力學(xué)，電磁學(xué)等各個領(lǐng)域廣泛應(yīng)用。然而，LBM 存在著對于計算資源的要求較高的缺點。具體來說，LBM 作為一個全顯式算法，為了保證計算的穩(wěn)定性以及精度，其計算的時間步長很短，而較短的時間步長導(dǎo)致計算迭代次數(shù)增多，使得對計算資源的消耗較高。本文從這一點出發(fā)，構(gòu)造模型使得單次的模型計算代替多步的LBM迭代，以加快LBM的計算速度。

近些年來有很多嘗試提升流體模擬計算效率的研究工作。有些研究提供了代替原來的求解泊松方程的方法，如Molemaker等［2］提出的迭代正交投影框架；Lentine 等［3］提出的粗網(wǎng)格投影法等。它們計算效率高，但它們的精度較低且僅對低分辨率問題有較好的效果。對于復(fù)雜的數(shù)值模擬問題，數(shù)據(jù)驅(qū)動模型是比較合適的優(yōu)化方案，即通過給求解器供給足夠的相關(guān)數(shù)據(jù)，使得求解器學(xué)習(xí)到該數(shù)據(jù)當中的統(tǒng)計信息。早期的數(shù)據(jù)驅(qū)動優(yōu)化主要是基于模型降維（ROM，reduced order methods）的思想，如Treuille等［4］將N-S方程投影到低維的子空間上，以簡化模型并提高計算效率，使求解器的計算速度能達到即時模擬的速度。早期ROM方法結(jié)構(gòu)較為簡單，而這也限制了其性能。

隨著計算機硬件的發(fā)展以及機器學(xué)習(xí)相關(guān)理論和算法的成熟，越來越多的人將機器學(xué)習(xí)運用到流體數(shù)值模擬的加速當中。這類方法利用傳統(tǒng)求解器計算得到的數(shù)據(jù)，訓(xùn)練得到一個機器學(xué)習(xí)模型，以代替原本耗時較長的計算步。如Ladick等［5］利用隨機森林擬合拉格朗日描述下的N-S 方程，可以很大程度上簡化原來的求解模型；另外他們還提出了一種針對平滑粒子法（SPH）的隨機森林模型。上述機器學(xué)習(xí)方法較為傳統(tǒng)，近年來發(fā)展起來的深層學(xué)習(xí)算法被證實能夠更好地解決復(fù)雜問題［6］。其中，卷積神經(jīng)網(wǎng)絡(luò)（CNN，convolutional neural network）可以很好地對流體模擬這樣的對空間信息十分依賴的問題進行優(yōu)化。

研究者們還將原本復(fù)雜的迭代計算或求解過程轉(zhuǎn)變?yōu)閳D像生成問題，如Yang 等［7］提出了基于CNN 的壓力求解器模型，用以代替原本求解N-S方程過程中迭代求解壓強的步驟；以此為基礎(chǔ)，Tompson等［8］提出了一種非監(jiān)督訓(xùn)練壓力求解模型的方法，這種方法只需要較少的訓(xùn)練樣本便可以得到求解模型。它們對于N-S 方程的求解有較好的提速效果且精度較高，但并不適合全顯式的LBM算法。

在對于LBM 算法的加速中，Guo 等［9］訓(xùn)練CNN 模型直接由流場邊界條件預(yù)測流場的穩(wěn)定狀態(tài)，這種方法加速效果明顯，精度較高，但這對于非穩(wěn)態(tài)問題以及需要關(guān)注流場發(fā)展過程的穩(wěn)態(tài)問題并不適用。Hennigh［10］用一步的CNN 模型運算來代替LBM 的多步計算，這項研究關(guān)注于壓縮內(nèi)存的使用以實現(xiàn)更大幅度的加速，但模型存在一定的誤差，在反復(fù)迭代中誤差累積更為明顯，生成的流場精度不高。

為適配非穩(wěn)態(tài)問題并提高精度，本文構(gòu)造了一個代替原來多個LBM 時間步的機器學(xué)習(xí)模型。在U-Net［11］基礎(chǔ)上，引入殘差神經(jīng)網(wǎng)絡(luò)（Res-Net，residual network）［12］以優(yōu)化訓(xùn)練的效果，并在損失函數(shù)當中引入了流場的物理信息以進一步提升模型精度。

1 LBM數(shù)值模型

在LBM 數(shù)值模擬當中，流場被劃分成很多等大的格子，利用每個格點上粒子分布函數(shù)f來描述該處流場的狀態(tài)，在格點上用ft表示運動方向為ei的粒子數(shù)。本文所使用的針對二維問題的D2Q9 模型將粒子的運動方向離散為包括零向量的9 個方向，如圖1所示。

圖1 二維計算模型D2Q9示例Fig.1 Calculation model in 2D(D2Q9)

利用粒子分布函數(shù)可以得到該位置的其他物理信息，如宏觀密度ρ和速度u為

其中r為格點位置，t為時間。LBM 離散形式的控制方程為

其中δt為時間步長，τC為松弛因子，Ω為由τC計算得到的碰撞因子。控制方程當中包含了對粒子遷移和碰撞兩方面的描述，這兩部分在數(shù)值程序當中分開計算，兩者的表達式為

其中式（3）表述粒子的碰撞過程，式（4）表示粒子的遷移過程。

2 基于U-Net的LBM方法

2.1 網(wǎng)絡(luò)結(jié)構(gòu)

為了對LBM 進行加速，本文所采用的加速方案為壓縮其計算迭代步，即通過給設(shè)計好的CNN提供大量的數(shù)據(jù)，訓(xùn)練出模型來代替原本多個時間步的計算。其時間步推進可歸納為

那么，可設(shè)計該加速方案的回歸模型為

其中n為時間步長的跨越尺度，b為邊界條件信息。

由這個回歸模型，可以設(shè)計一個輸入為t時刻的分布函數(shù)，輸出為t+ Δt時刻分布函數(shù)的深層回歸網(wǎng)絡(luò)。對于本文所使用的D2Q9 模型，該網(wǎng)絡(luò)的輸入為10 個通道的一個“圖像”，包括了9 個方向的粒子分布函數(shù)以及一個邊界通道，輸出為9個通道的“圖像”，包括9 個方向的粒子分布函數(shù)。對于這類圖像生成問題，往往采用的是編碼-解碼格式：先編碼，即通過卷積（Conv）、池化（Pooling）進行深層次的特征提?。辉俳獯a，即通過反卷積（TransConv）或上采樣（Upsampling）利用特征計算得到目標輸出圖像。但在實踐當中，該問題極強的非線性使得這個模型的回歸十分困難。為了提高模型精度，本文做了一些嘗試，成功地將模型精度控制在較高的水平。本文構(gòu)造的卷積神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)如圖2所示。

圖2 基于U-Net對LBM加速的CNN網(wǎng)絡(luò)結(jié)構(gòu)Fig.2 CNN model to accelerate LBM based on U-Net

在實現(xiàn)這個編碼-解碼的結(jié)構(gòu)時，卷積操作使用3 × 3 的卷積核，反卷積操作使用2 × 2 的卷積核，訓(xùn)練當中使用的優(yōu)化算法為Adam［13］，卷積層激活函數(shù)的選取都為ReLu［14］，在訓(xùn)練過程中設(shè)置了學(xué)習(xí)率的衰減，使得模型在訓(xùn)練的后期收斂更穩(wěn)定。網(wǎng)絡(luò)結(jié)構(gòu)參考了Olaf 等［11］提出的U-Net 網(wǎng)絡(luò)結(jié)構(gòu)。其關(guān)鍵在于在解碼的過程中，將同等級的編碼信息補充到解碼信息當中去。該網(wǎng)絡(luò)結(jié)構(gòu)起初用于醫(yī)學(xué)圖像分割，起到了非常好的效果。之后，很多人將此網(wǎng)絡(luò)結(jié)構(gòu)用于其他的圖像生成問題當中去，很大程度地提高了模型精度。隨著神經(jīng)網(wǎng)絡(luò)層數(shù)的加深，一些關(guān)鍵的特征信息可能會丟失，而U-Net的連接層可以將這部分特征信息補充回來，很好地提升了圖像回歸問題的精度。

為顯示U-Net對模型精度的提升作用，本文訓(xùn)練了另一組沒有采用U-Net 的編碼-解碼結(jié)構(gòu)以進行比較，發(fā)現(xiàn)精度的提升是巨大的：在訓(xùn)練300輪之后，U-Net 模型的平均百分比誤差可以達到0.138 9%，而普通的編碼-解碼結(jié)構(gòu)只能達到1.023 7%；且在實際的使用當中，普通的編碼-解碼結(jié)構(gòu)生成的模型并不能維持住邊界信息，且其生成的流場失真很嚴重。而僅僅應(yīng)用了U-Net的網(wǎng)絡(luò)結(jié)構(gòu)也存在著一定的問題，單純的深層網(wǎng)絡(luò)很容易產(chǎn)生梯度彌散或梯度爆炸的問題。加了正則化層［15］的深層網(wǎng)絡(luò)可能會解決這種問題，但它不能解決深層的網(wǎng)絡(luò)經(jīng)常會出現(xiàn)的模型退化問題。對此，本文參考了He 等［12］提出的ResNet 對卷積層進行了優(yōu)化，見圖2。

殘差神經(jīng)網(wǎng)絡(luò)被運用在了圖像識別問題上，很大程度上提高了模型的精度，其關(guān)鍵在于Shortcut 連接層。理論上，多層的神經(jīng)網(wǎng)絡(luò)可以擬合任意函數(shù)，但過于深層的網(wǎng)絡(luò)很容易導(dǎo)致模型退化而使得模型無法收斂。通過Shortcut 連接層，可以在層數(shù)“多余”的時候?qū)⒍嘤嗟膶佣探拥?，以此避免模型的退化問題。在訓(xùn)練的過程中，本文將殘差塊應(yīng)用到流場的編碼與解碼過程中，發(fā)現(xiàn)殘差神經(jīng)網(wǎng)絡(luò)對于訓(xùn)練效果的提升是很大的，提高了模型精度以及收斂速度。

2.2 損失函數(shù)

損失函數(shù)是模型訓(xùn)練過程中用來評價模型優(yōu)劣以幫助模型訓(xùn)練的工具。對于回歸問題，往往使用的是L2的損失函數(shù)（MSE，mean squared error）

其中i表示通道數(shù)，對應(yīng)D2Q9中9個方向，下標為“pred”的值為模型預(yù)測結(jié)果，下標為“true”的值為真值。

L2損失函數(shù)在回歸問題當中效果是比較不錯的，但它并不包含一些非常關(guān)鍵的物理信息。針對LBM，本文將一些關(guān)鍵的物理信息包含到了損失函數(shù)當中。針對不可壓粘性流體密度不變的性質(zhì)，可以寫lossρ為

另外，針對不可壓粘性流體流場速度散度為0的性質(zhì)，可以寫

其中

所以，本文的損失函數(shù)為

其中μ1，μ2，μ3為系數(shù)。

對于相同的訓(xùn)練集，利用不同的損失函數(shù)進行訓(xùn)練。如圖3所示，最終得到的兩模型中包含了物理信息的LBM損失函數(shù)（μ1= 1，μ2=μ3= 0.05）精度略高于MSE 損失函數(shù)。由于該模型在進行數(shù)值模擬時需要進行反復(fù)地迭代，這種微小的精度提升對模型的穩(wěn)定性會有一定的提升。

圖3 模型訓(xùn)練過程中損失函數(shù)值Fig.3 Loss function of the model during training

3 算例實驗

3.1 模型訓(xùn)練

本文設(shè)計用以證明模型有效性的算例為層流繞柱群的流場模擬。這種算例有比較復(fù)雜的漩渦結(jié)構(gòu)，能夠挑戰(zhàn)模型學(xué)習(xí)復(fù)雜流場信息的能力，也能較為直觀地看出模型的誤差。算例的左側(cè)為0.02 m/s的流入邊界，右側(cè)為流出邊界，上下設(shè)置為了周期性邊界，障礙物為無滑移碰撞邊界。訓(xùn)練集為20 組層流流過三個并排的幾何形狀不同（圓、橢圓和正方形）、大小及位置隨機的障礙物（如圖4所示），這樣的設(shè)置可以使模型對于障礙物的形狀以及位置更敏感，提高模型的泛化能力。其余的參數(shù)設(shè)置如表1。

表1 訓(xùn)練集LBM仿真參數(shù)設(shè)置Table 1 Parameter setting for LBM simulation training set

圖4 部分訓(xùn)練集Fig.4 Part of the training data set

設(shè)置的跳步數(shù)量為n= 200，即模型的回歸模型為

對于本文所涉及的U-Net 模型，輸入為t時刻的粒子分布函數(shù)和流場邊界信息，輸出為t+200Δt時刻的粒子分布函數(shù)。其中，對于本文后面所涉及的算例，流場中僅包含了碰撞邊界，那么輸入方式為二進制數(shù)組，在有障礙物的地方該數(shù)組的值為1，其他地方為0. 圖5為一個二維的正方形障礙物邊界例。

圖5 幾何邊界輸入Fig.5 Geometry boundary input

作為數(shù)據(jù)的預(yù)處理，將輸入到訓(xùn)練當中的所有數(shù)據(jù)都進行了放大處理，因為流體系統(tǒng)是一個非常精細的系統(tǒng)，其臨近點之間的差異可能并不是很大，如果直接將沒有經(jīng)過處理的數(shù)據(jù)輸入到訓(xùn)練當中去，有很大的可能會造成模型最終收斂得到整個流場為一個平均值。

由于計算的前1 000 個時間步流場的變化往往比較小，可能會影響訓(xùn)練的效果，所以訓(xùn)練集采用的為第1 000 個時間步之后的計算結(jié)果，那么在使用模型的時候也是先利用LBM 計算1 000 步之后再調(diào)用訓(xùn)練好的模型。訓(xùn)練時所使用的批大小為50，訓(xùn)練300輪得到最終模型。

3.2 誤差分析

首先，對本文所提出的損失函數(shù)lossLBM做誤差分析。包含一定的物理信息勢必會使得模型更貼近真實的物理狀況，在實際的訓(xùn)練當中也可以看出精度的提升，如表2所示。

表2 訓(xùn)練至300輪時模型精度比較Table 2 Accuracy of the model after 300 epochs of training%

若系數(shù)適當，該損失函數(shù)的效果略微優(yōu)于MSE，通過對系數(shù)的進一步優(yōu)化，可達到更好的結(jié)果。而系數(shù)設(shè)置不合理可能會使得模型的訓(xùn)練更為困難，這是因為本文所引入的兩個物理信息在數(shù)量級上略大于MSE，比較大的系數(shù)會使得模型更難以收斂。

下面在實際的算例當中調(diào)用模型，來看該模型的性質(zhì)。對于訓(xùn)練集內(nèi)部有的算例（算例1）來說，U-Net 模型的表現(xiàn)非常不錯，計算結(jié)果如圖6所示。直觀來看，U-Net 模型的計算結(jié)果幾乎與LBM 的計算結(jié)果完全一致。為更清晰地看出模型的誤差發(fā)展，利用下式計算模型的誤差，即其中i表示通道數(shù)，對應(yīng)D2Q9中9個方向，j和k分別表示點x方向位置以及y方向位置，nx和ny分別表示x方向格點數(shù)及y方向格點數(shù)，nlattice為總格點數(shù)。

圖6 算例1流場速度圖（第一行為LBM計算結(jié)果，第二行為U-Net模型的計算結(jié)果）Fig.6 Velocity field of example 1(the first line is generated by LBM,the second line is generated by U-Net model)

由圖7所示，隨著對于模型的迭代，結(jié)果的誤差增高，但其精度仍控制在較高水平。經(jīng)過30 次的迭代后，計算得到的粒子分布函數(shù)的相對誤差在1%以下，絕對誤差在0.000 7 左右。在第30 個迭代步（對應(yīng)的LBM 模型的計算步為第7 000 步），取流場中心一列做速度分布圖，見圖8。從圖中可以看出，CNN 模型可以很好地學(xué)習(xí)到流場的流動規(guī)律，計算所得流場的速度水平分量和豎直分量都能較好地與LBM 數(shù)值模擬的結(jié)果相吻合，U-Net模型對于訓(xùn)練集當中所有的算例精度很高。

圖7 算例1誤差發(fā)展圖Fig.7 Development of error of example 1

圖8 算例1迭代30次后流場中心線上的速度分布Fig.8 Distribution of flow field centerline velocity of example 1 after 30 times of iteration

對于訓(xùn)練集以外，但是相似度與訓(xùn)練集當中算例較高的算例，以三個障礙物均為正方形為例（算例2）。該算例是訓(xùn)練集當中所沒有的，計算結(jié)果見圖9。對于該算例來說，雖然可以直觀地看出U-Net模型生成流場的誤差，但其生成的流場還是比較真實的。

圖9 算例2流場速度圖（第一行為LBM計算結(jié)果，第二行為U-Net模型的計算結(jié)果）Fig.9 Velocity field of example 2(the first line is generated by LBM,the second line is generated by U-Net model)

該算例的誤差發(fā)展如圖10 所示，對于更加泛化的算例，模型的計算結(jié)果誤差相較于訓(xùn)練集中的算例（算例1）稍高一些，經(jīng)過30次的迭代，其相對誤差積累到了1.5%左右。在第30 個迭代步（LBM 的計算步為第7 000 步），取流場中心一列做速度分布圖（見圖11）。從圖11可以看出，模型計算得到的水平速度分量精度較高，與LBM 計算結(jié)果吻合地較好；豎直的速度分量與LBM 計算結(jié)果之間有一定的誤差，但二者趨勢一致，且誤差在數(shù)值上較小，對整體精度的影響有限。綜合以上，該模型精度較高，且有一定的泛化能力。

圖10 算例2誤差發(fā)展圖Fig.10 Development of error of example 2

圖11 算例2迭代30次后流場中心線上的速度分布Fig.11 Distribution of flow field centerline velocity of example 2 after 30 times of iteration

3.3 加速效果分析

本文所提出神經(jīng)網(wǎng)絡(luò)回歸模型的目的為提升計算效率，將它利用在一些不需要很高精度且對計算速度有一定要求的算例當中。因此，我們需要對U-Net模型的計算速度進行一個評估。

對LBM 方法運算表現(xiàn)的評判標準為每秒百萬格子更新數(shù)，其表達式為

其中nl為格子數(shù)，T為計算耗時。

本文的測試平臺參數(shù)為：CPU：Intel（R）Xeon （R） W-3175X CPU @ 3.10GHz； GPU：RTX2080Ti 。實現(xiàn)算法的編程語言為Python 3.7，神經(jīng)網(wǎng)絡(luò)回歸模型的訓(xùn)練和預(yù)測是通過深度學(xué)習(xí)開源庫Keras 所實現(xiàn)的。LBM 的浮點數(shù)計算為雙精度計算，而由于雙精度計算的CNN 模型幾乎不會有任何的精度提升，且會大幅降低計算速度，本文CNN模型采用的是單精度浮點數(shù)。

在調(diào)用模型時進行計時，測定神經(jīng)網(wǎng)絡(luò)模型的加速效果。以格子數(shù)為1282的二維算例為測試算例。調(diào)用模型，并將其與LBM計算200個時間步進行比較；與一些其他算例下利用GPU 并行計算對LBM 進行加速的計算效率進行比較，結(jié)果見表3?？芍海?）根據(jù)這個標準計算， U-Net 模型的計算速度相對于GPU 的并行計算還是要快一些。但需要注意到的是，表3中兩個GPU并行的測試算例與本文所選取的算例不一致。在進一步的測試中發(fā)現(xiàn)，U-Net模型對于更大尺度上的算例有著更好的加速效果，測試結(jié)果如圖12 所示。（2）該模型在更大尺度的流場上加速效果會更好，MLUPs 最高大致在1 600左右。

圖12 計算效率與模型尺度的關(guān)系Fig.12 Relationship between calculation efficiency and scale of the model

表3 計算效率Table 3 Calculation efficiency

綜上，該模型能顯著加快LBM 的計算，尤其在大尺度的數(shù)值模擬當中，得到的加速效果比較好。

4 結(jié) 論

本文提出了一種基于U-Net 的對LBM 進行加速的卷積網(wǎng)絡(luò)神經(jīng)結(jié)構(gòu)，引入殘差神經(jīng)網(wǎng)絡(luò)預(yù)防了深層神經(jīng)網(wǎng)絡(luò)的退化，并設(shè)計了一個包含物理信息的損失函數(shù)，在一定程度上提高了模型的精度。為驗證該模型的有效性，本文進行了數(shù)值實驗。其結(jié)果表明，本文所提出的模型能夠在保證精度較高的同時，實現(xiàn)很大程度的加速。另外，該模型是具有一定的泛化能力的。