基于格子Boltzmann方法和大渦模擬的頸動脈分叉狹窄流動并行計算

張毅卓，葛 森，王良軍，謝 江，曹 潔，張 武，*

（1.上海大學上海市應用數學與力學研究所，上海200444；2.上海大學計算機工程與科學學院，上海200444；3.上海交通大學附屬仁濟醫院老年科，上海310101）

0 引言

頸動脈狹窄是頸動脈疾病的重要病因之一，常見的病變包括動脈粥樣硬化、頸動脈斑塊等，斑塊的破裂導致缺血性腦卒中的產生［1-2］。隨著病變區域的發展，頸動脈血管變得狹窄，從而阻礙血液流動，使得頸動脈血管內部的血流動力學發生復雜變化。

血管的壁面剪切應力（Wall Shear Stress，WSS）被認為是病變發生和發展的一個重要因素，并被廣泛研究［3］。血流動力學的異常變化將會導致斑塊和狹窄部位區域產生異常的WSS。較低的WSS是頸動脈粥樣硬化病變早期發展的主要因素之一，剪切力的作用導致血管壁受到損傷，血小板的活性被激活，在斑塊壁面出現局部聚集，引起頸動脈的進一步狹窄［4］。嚴重狹窄部位的斑塊肩部會產生高WSS，使斑塊可能出現裂縫、破裂和脫落，造成嚴重后果。

在傳統的有限元和有限體積計算方法中，為了根據特定患者的臨床數據計算WSS，首先基于超聲獲得患者的血流速度，將利用核磁共振成像（Magnetic Resonance Imaging，MRI）獲得的患者頸動脈血管幾何形狀轉換為標準有限元/有限體積求解器所需的三維網格［5］，并導入求解器進行求解，計算速度場并在血管壁面進行梯度計算。這種方法存在的問題是模型網格生成的預處理和求解器計算時間過長，限制了其在實際患者診斷中的適用性。

頸動脈處的血液流動是復雜的牛頓粘性流體運動，一般可以通過N-S（Navier-Stokes）方程來進行描述。格子Boltzmann 方法（Lattice Boltzmann Method，LBM）是一種針對流動問題求解的有效方法［6］。基于Boltzmann 方程，它的物理基礎主要包括兩部分：一是在介觀層面考慮粒子運動和對粒子之間的遷移、碰撞進行近似與簡化；二是通過時間、空間和粒子遷移方向的離散，使用離散速度模型來模擬流體運動。LBM 通過Boltzmann方程描述粒子分布函數隨時間的演化，通過粒子分布函數的求和可以得到流體的速度、密度和壓力等宏觀量。通過對時間、空間和粒子遷移方向的離散化，LBM具有天然的可并行性、容易編程和高精度等優點，能夠在并行計算機上進行有效計算，非常適合處理復雜幾何形狀和復雜流動的大規模計算問題［7］。

LBM在頸動脈數值模擬領域已有初步的研究。利用多松弛時間格子Boltzmann（Multi-Relaxation-Time Lattice Boltzmann，MRT-LB）方法模擬了三維狹窄和正常頸動脈分叉的脈動流模式［8］。對血液的牛頓模型和非牛頓模型應用LBM，并在穩態流動和脈動流動的背景下進行了研究。分析了WSS 對動脈粥樣硬化的影響，證明LBM 適用于研究動脈粥樣硬化發展引起變化的流體［9］。給出了LB 模型中三維WSS計算的公式，采用改進的算法，計算了人體腹主動脈壁剪切應力等［10］。

上述研究都是在較低雷諾數（Re=500～600）下進行的。考慮到血液是一種典型的脈動流動，動脈處的狹窄大大增加了產生復雜流動的可能性，包括考慮對應更高Re的流動。復雜流動會對狹窄部位下游區域的壓降和WSS 產生顯著影響，兩者都會對頸動脈壁面產生實質性的負面影響。從精確計算的角度來看，采用經典LBM 模擬這種流動面臨較大的挑戰［11］。

大渦模擬（Large Eddy Simulation，LES）方法已被證明對復雜流動具有較強的數值模擬能力［12］。在LES 中，用非穩態的Navier-Stokes 方程來直接模擬大尺度渦，而較小尺度的渦對大渦的影響則通過近似模型來模擬，采用的模擬方法稱為亞網格尺度（SubGrid Scale，SGS）模型［11，13］。

LES 有許多優點，在SGS 模型的幫助下，空間和時間分辨率可以進一步增加。LES 提供了關于流動的大渦的精確信息，能夠提供更好的物理洞察力，并有潛力成為更準確的預測工具［14］。通過LES 計算主動脈內的擾動流場和壁面剪切應力，使用磁共振成像（MRI）的臨床數據作為輸入的模型，可以很容易地找到主動脈壁上較高WSS 的位置［15］。Vergara 等［16］通過研究腹主動脈瘤向湍流的過渡，把臨床和計算結果加以對比，表明LES模型具有較好的適用性。因此，在頸動脈狹窄流動數值模擬中考慮大渦模擬，將LBM 與LES 相結合，充分利用LBM 的并行特性和LES 對復雜流動的精確描述，以減少計算時間，提高計算精度，有利于更好地理解頸動脈狹窄后部的血液流動，可能會對臨床上提供更準確的診斷指導。

LBM的計算時間通常隨著網格尺度的減小和雷諾數的增加而增加。但隨著多核處理器的快速發展，在大規模并行架構上采用OpenMP 并行編程模型能夠顯著減少程序求解的運行時間，提高模擬速度［17-18］。

本文的主要工作是提出了LBM 結合LES 的并行算法——LBM-LES頸動脈模擬算法來對頸動脈狹窄部位的脈動流動進行數值模擬；并針對計算時間較長的挑戰，采取OpenMP模型對該算法進行并行加速［19］。血管幾何模型由MRI臨床數據建模得到，頸動脈脈動血流不同時刻的速度則由B 超檢查得到，根據實測參數和幾何模型，通過LBM-LES 并行算法模擬出不同脈動時刻的三維頸動脈的血流速度和WSS。最后分析了實驗結果，評估了不同核數下的計算時間和加速比。

1 格子Boltzmann方法

LBM 源于格子氣自動機，通過Boltzmann 方程，對時間和空間進行離散化，用離散的模型從介觀上模擬復雜的流動現象。LBM 的離散模型通常定義為DmQn，下標表示m 維空間和格子上n個離散速度方向。

LBM 在相空間將粒子的速度u 簡化為有限維的速度空間{e0，e1，…，en-1}，引入了n 個分布函數fi(x，t)(i=0，1，…，n-1)來描述離散時間步長t 下沿ei方向流動的粒子。經典的LBGK（Lattice Bhatnagar-Gross-Krook）演化方程［6，20］具有以下形式：

其中：Rg為流體常數，T 為流體溫度。無量綱松弛時間τ 與聲速cs和流體的運動粘度ν 有關，可以用表示。格子單位Δx、時間步長Δt和流體密度ρ 的單位分別為m、s 和kg/m3。

以D3Q19格子模型為例來說明LBM。如圖1所示，圖中0，1，…，18 對應表示e0，e1，…，e18，其中e0位置是流體粒子的初始位置，e1，e2，…，e18指定粒子在碰撞后流向的方向。

三維模型下所設的微觀速度為：

圖1 D3Q19離散速度模型Fig.1 D3Q19 discrete velocity model

根據LBGK方程，LBM的演化可分為以下兩個步驟：

1）碰撞。

2）遷移。

在格子單元中，時間步長和空間步長通常相同，因此c=Δx/Δt=1。格子聲速取值，當速度較低時，局部平衡態分布函數可以采用Maxwell 分布的二階Taylor 展式，模型的分布函數如下：

相應的權重系數為：

2 壁面剪切力的數值解

血液在血管中流動時，會對血管壁表面施加作用力，壁面剪切應力表示血液對血管壁面單位面積上所施加的剪切力。對于頸動脈粥樣硬化和斑塊部位，血管被設定為剛性壁面，血流可以被精確地建模為均勻介質、不可壓縮的牛頓流體，對于直徑較大的血管，如頸動脈，這通常被認為是可接受的［9］。

當物體在血管壁面流動時，壁面速度變為零，無滑移邊界條件成立。對于不可壓縮的牛頓流體，還需要血管壁面的三維空間信息和血管內血流速度，以及相應的壁面法向量［3］。WSS可以通過剪切力張量σ和邊界法向單位矢量n確定：

D3Q19模型的壁面剪切應力計算公式［10］如下：

其中，α，β，γ代表三維空間內不同的維度。

3 LBM-LES并行算法

3.1 大渦模擬

大渦模擬的精度很高，能夠模擬瞬時流場運動。考慮到脈動血流的影響，同時為了精確描述流體運動，本文采用大渦模擬結合LBGK模型來對頸動脈狹窄流動進行數值模擬。

LES 的基本思想是采用濾波器控制方程進行濾波，從而把復雜的流動現象分解為大尺度結構和小尺度結構。在數值模擬時對大尺度的流動現象進行直接模擬，針對小尺度脈動采用SGS模型進行模擬，模型的網格尺度可以設為大尺度，在格子內采用SGS模型來進行更小尺度的模擬。在采用適當的亞格子模型的情況下，這樣能夠模擬真實物理情況下的瞬時流場。

在LBGK 模型中加入LES 的SGS 模型，采用經典的Smagorinsky 模型［11-12，21］。在LBM-LES 方法中，Smagorinsky 模型的引入是通過改變原有LBM 碰撞模型中的無量綱松弛時間τ來實現的。

引入LES后的有效粘性可以表示為：

其中：v0是層流粘性，vt是渦粘性。由Smagorinsky 模型給出渦粘性即為SGS模型。

在Smagorinsky模型中，渦粘性vt由應變率張量

和濾波尺度Δ確定，如式（13）所示：

其中：Cs是Smagorinsky 參數，|S|為應變率張量Sαβ的模。|S|在LBM模型中可以通過下式獲得：

3.2 并行算法實現

LBM-LES 并行算法具有LBM 的可并行化優點，可以采用OpenMP 模型進行加速。本文采用的OpenMP 是用于共享內存機器的一套并行編程庫。程序在運行時采用fork/join 的并行模式。首先主程序執行串行代碼，在執行到包含#parallel 編譯指令的代碼段時，主進程派生出子線程分配給多處理器的不同計算核來處理計算任務，實現線程級別的并行計算，子線程執行完畢后銷毀，繼續執行主進程的任務。

為了更好地描述LBM-LES 并行算法，給定程序輸入為：網格模型和計算參數；程序的輸出為：網格格點數據文件，包括格點坐標、流速分量、剪切力分量和壓力信息。

本文所采取的LBM-LES并行算法的具體計算步驟如下：

步驟1 設定好頸動脈計算區域和最小網格尺寸，讀入面網格文件；

步驟2 判斷血管邊界，遍歷標記計算網格；

步驟3 讀入數值計算參數，包括速度向量、雷諾數、特征長度和迭代步數等；

步驟4 遍歷網格點，進行計算數據初始化；

步驟5 在LBGK模型碰撞過程中計算Smagorinsky參數，引入LES；

步驟6 模型遷移過程；

步驟7 宏觀量計算過程，得到宏觀速度、壓力、剪切力、以及對應速度三維分量、剪切力三維分量；

步驟8 若滿足設置的輸出迭代步數條件，輸出中間計算結果；

步驟9 判斷是否收斂，若不滿足收斂條件，重復步驟5～7，否則結束，輸出程序計算信息，終止計算。

根據LBM-LES 模型和OpenMP 編程模型給出整體的并行算法，如圖2 所示。首先，程序讀入光固化立體造型術（STereoLithography，STL）網格文件，在網格的初始化過程中，由于每個格子的不相關性，每個格子類型的判斷以及格點上微觀量的初始化都可以進行并行操作，加快處理速度。在迭代過程中，網格的碰撞和遷移過程都可以采用合適的OpenMP指令來實現。同時，宏觀量的計算只和格子自身存儲數據有關，因此也可以采用并行算法來進行加速。直到最后滿足程序終止條件，輸出迭代結果。

圖2 本文并行算法流程Fig.2 Flowchart of the proposed parallel algorithm

4 數值實驗與分析

4.1 實驗環境

實驗在上海大學高性能計算集群“自強4000”上進行并行性能測試，所用共享內存節點具體配置信息如下：4 顆Intel Sandy Bridge 架構的E5-4650 CPU（2.7 GHz/8-core），512 GB RDIMM DDR3 1 600 MHz 共享內存，計算隊列共計32 個計算核數，編譯命令采用“-O3”選項進行優化。

4.2 計算模型

程序讀入的網格文件為真實頸動脈血管模型，通過從合作醫院臨床采集得到的MRI 圖像數據，應用醫學后處理軟件MimiCS 軟件進行三維重構，采用閾值分割的方法分離血管周圍組織，從而得到頸動脈血管真實幾何模型［22］，并導出為STL格式的網格文件。頸動脈三維幾何模型如圖3 所示，血流方向從下至上，血流從頸總動脈（Common Carotid Artery，CCA）入口流入，從頸外動脈（Externl Carotid Artery，ECA）、頸內動脈（Internal Carotid Artery，ICA）流出。計算區域的高度為CCA入口和ECA、ICA出口之間70.0 mm的區域。

由于計算的真實血管形狀非常復雜，要精確描述血管壁面區域的不規則形狀，需要在適當的位置采用階梯逼近或者插值處理，以保證滿足物理邊界上的條件。因此，邊界條件采用YMS 曲線插值處理［23］，在每次碰撞遷移的過程中，需要在每個與物理邊界相交的網格線上進行插值處理，血管壁面均采用無滑移邊界條件。

為了模擬脈動流動，在一個脈動周期內，在不同的時間點來改變CCA 入口的血流速度，分別模擬最小舒張速度（Velocity Minimum，VM）、最 大 收 縮 期 速 度（Velocity Maximum，VMAX）、雙側切速度（Velocity in Dicrotic Notch，VDN）和舒張流動后期速度（Velocity in Diastolic Flow pHase，VDFH）。在具體計算過程中，分別在參數初始化（迭代步數為0），迭代步數為800、1 200、1 600 的時刻改變頸動脈入口的流速，來模擬脈動流動的情況。表1 給出了三維算例不同時間點的具體計算數據。

圖3 三維血管模型Fig.3 Three-dimensional carotid model

表1 脈動血流參數Tab.1 Pulsating blood flow parameter

4.3 數值實驗

通過對不同雷諾數下的頸動脈狹窄流動計算，圖4、5 分別給出了雷諾數為500（LBGK模型）、2 000（LBM-LES模型）情況下在不同時刻（迭代步數800、1 200、1 600、2 000）的流速云圖和剪切力矢量圖。

4.3.1 流速矢量圖

如圖4 所示，圖（a）～（d）分別代表四個不同時刻的頸動脈分叉部位X軸方向的截面圖。對于存在斑塊區域的狹窄頸動脈血管模型，血流在流過彎曲部位導致血流形態發生變化，在ICA 處具有較大的血流速度，同時斑塊后段產生一個小的血流較低的區域，這種流動情況會導致斑塊的繼續生長，病變區域繼續向下游擴展，隨著斑塊進一步進行擴大，頸動脈的狹窄變得更加嚴重。當Re=2 000 時，結合LES，在不同時刻的截面中，血液流動變得復雜，在收縮期間，斑塊附近的流動區域具有較大的流速，在狹窄段的下游區域，血流流過狹窄區域之后流速上升。

4.3.2 剪切力矢量圖

圖4 流速云圖Fig.4 Velocity cloud chart

本小節對血管WSS 進行分析。如圖5 所示，根據壁面剪切力的矢量云圖可以很明顯看出：在收縮期之前，斑塊處WSS較小（圖（a）），隨著收縮期，血流速度升高，斑塊肩部產生較大的WSS（圖（b）），在狹窄段的下游區域兩側，具有高WSS；隨著收縮期結束，斑塊處WSS 開始減小，在剛開始進入舒張流動后期，斑塊處WSS 較高。而在較高雷諾數的情況下，同樣參數條件的WSS 比低雷諾數情況下要高，在斑塊肩部產生高WSS，在斑塊尾部產生一個WSS 尖點，高WSS 會導致血管壁面的損傷，進而促進斑塊的擴張。

圖5 剪切力矢量圖Fig.5 WSS vector diagram

4.3.3 并行時間和加速比

本文計算程序在高性能集群共享內存節點上進行并行性能的測試。在提交任務隊列時設置運行參數，分別采用4 個、8個、16個、24個和32個核來迭代計算2 000步，計算網格量固定，數值實驗算例的網格總數固定為4 000 萬，所得運行時間數據為三次計算結果的平均值，包括網格生成、初始化、碰撞遷移和宏觀量計算過程。計算時間并未考慮計算結果的輸出時間。圖6給出了不同計算核數下程序的運行時間，圖7給出了不同計算核數的程序加速比。

圖6 LBM和LBM-LES模型在不同核數下的運行時間對比Fig.6 Running time comparison of LBM and LBM-LES models with different cores

圖7 LBM和LBM-LES模型在不同核數下的加速比對比Fig.7 Speedup comparison of LBM and LBM-LES models with different cores

從圖6 可以看出，LBM-LES 的計算時間要略高于原有LBM 模型，這是由于額外計算松弛因子所增加的計算開銷。隨著計算核心數量的增加，程序的運行時間逐漸減少。從圖7 可以看出，加速比的增長率在計算核數4～8，8～16，16～32 的階段呈現出先上升后下降的趨勢，隨著核數的增加，性能增益相對降低，兩個不同模型之間的加速比相對較符合。這很可能是由于集群共享內存節點一塊E5-4650 CPU具有8個核心，不同CPU 之間的計算核之間數據通信會對程序的加速效果產生影響，數據的傳遞和通信是耗時的，結果低于理想加速比。同時，隨著計算網格量的增加，OpenMP 并行編程模型的并行性能逐步降低，采用MPI進行跨節點的并行加速，將會進一步提高并行性能。

并行計算實驗結果顯示了在高性能集群上采用并行加速算法能夠有效縮短LBM-LES 算法的計算時間，如何充分發揮大型超級計算機的并行計算能力，進行跨節點MPI 編程模式的大規模問題的LBM-LES 并行計算，深入開展頸動脈狹窄血液流動精確計算，是今后進一步的研究內容。

5 結語

本文介紹了LBM 結合大渦模擬模型在模擬動脈狹窄血液流動研究中的應用，驗證了大渦模擬模型適合于頸動脈復雜脈動流動的模擬和分析，然后對模擬的程序進行了優化與并行。

針對臨床數據，采用三維模型重構方法構建頸動脈模型，基于LBM-LES 模型，計算了得到的三維真實血管，在一個脈動周期內，給出了特征時間點的速度云圖、WSS矢量云圖的瞬時分布，模擬的流動在不同位置顯示出許多特征，這些信息對于理解復雜的動脈血流和診斷動脈疾病方面可能是非常有幫助的。利用OpenMP 編程模型對程序進行并行加速，明顯地縮短了程序計算時間，并實現了較好的加速比。

通過快速計算實際脈動流動的復雜狀況，為個體頸動脈內進行體外非侵入式仿真模擬血流動力學分析提供了一種可行的方法，通過并行加速能夠在較短時間為醫生診斷頸動脈斑塊疾病提供一定的理論依據和輔助作用，為頸動脈斑塊體外輔助診治提供幫助。