基于U-Net的格子玻爾茲曼方法*

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

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

格子玻爾茲曼方法（LBM）［1］于20 世紀(jì)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ù)當(dāng)中的統(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ù)值模擬的加速當(dāng)中。這類方法利用傳統(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ù)當(dāng)中引入了流場的物理信息以進一步提升模型精度。

1 LBM數(shù)值模型

在LBM 數(shù)值模擬當(dāng)中，流場被劃分成很多等大的格子，利用每個格點上粒子分布函數(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計算得到的碰撞因子?？刂品匠坍?dāng)中包含了對粒子遷移和碰撞兩方面的描述，這兩部分在數(shù)值程序當(dāng)中分開計算，兩者的表達式為

其中式（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）進行深層次的特征提??；再解碼，即通過反卷積（TransConv）或上采樣（Upsampling）利用特征計算得到目標(biāo)輸出圖像。但在實踐當(dāng)中，該問題極強的非線性使得這個模型的回歸十分困難。為了提高模型精度，本文做了一些嘗試，成功地將模型精度控制在較高的水平。本文構(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)練當(dāng)中使用的優(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)鍵在于在解碼的過程中，將同等級的編碼信息補充到解碼信息當(dāng)中去。該網(wǎng)絡(luò)結(jié)構(gòu)起初用于醫(yī)學(xué)圖像分割，起到了非常好的效果。之后，很多人將此網(wǎng)絡(luò)結(jié)構(gòu)用于其他的圖像生成問題當(dāng)中去，很大程度地提高了模型精度。隨著神經(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%；且在實際的使用當(dāng)中，普通的編碼-解碼結(jié)構(gòu)生成的模型并不能維持住邊界信息，且其生成的流場失真很嚴(yán)重。而僅僅應(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個方向，下標(biāo)為“pred”的值為模型預(yù)測結(jié)果，下標(biāo)為“true”的值為真值。

L2損失函數(shù)在回歸問題當(dāng)中效果是比較不錯的，但它并不包含一些非常關(guān)鍵的物理信息。針對LBM，本文將一些關(guān)鍵的物理信息包含到了損失函數(shù)當(dāng)中。針對不可壓粘性流體密度不變的性質(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)練當(dāng)中的所有數(shù)據(jù)都進行了放大處理，因為流體系統(tǒng)是一個非常精細的系統(tǒng)，其臨近點之間的差異可能并不是很大，如果直接將沒有經(jīng)過處理的數(shù)據(jù)輸入到訓(xùn)練當(dāng)中去，有很大的可能會造成模型最終收斂得到整個流場為一個平均值。

由于計算的前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)練當(dāng)中也可以看出精度的提升，如表2所示。

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

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

下面在實際的算例當(dāng)中調(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)練集當(dāng)中所有的算例精度很高。

圖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)練集當(dāng)中算例較高的算例，以三個障礙物均為正方形為例（算例2）。該算例是訓(xùn)練集當(dāng)中所沒有的，計算結(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ò)回歸模型的目的為提升計算效率，將它利用在一些不需要很高精度且對計算速度有一定要求的算例當(dāng)中。因此，我們需要對U-Net模型的計算速度進行一個評估。

對LBM 方法運算表現(xiàn)的評判標(biāo)準(zhǔ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ù)這個標(biāo)準(zhǔn)計算， 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ù)值模擬當(dāng)中，得到的加速效果比較好。

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)很大程度的加速。另外，該模型是具有一定的泛化能力的。