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

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

0 引言

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

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

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

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

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

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

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

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

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

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

1 格子Boltzmann方法

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

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

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

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

三維模型下所設(shè)的微觀速度為：

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

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

1）碰撞。

2）遷移。

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

相應的權(quán)重系數(shù)為：

2 壁面剪切力的數(shù)值解

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

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

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

其中，α，β，γ代表三維空間內(nèi)不同的維度。

3 LBM-LES并行算法

3.1 大渦模擬

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

LES 的基本思想是采用濾波器控制方程進行濾波，從而把復雜的流動現(xiàn)象分解為大尺度結(jié)構(gòu)和小尺度結(jié)構(gòu)。在數(shù)值模擬時對大尺度的流動現(xiàn)象進行直接模擬，針對小尺度脈動采用SGS模型進行模擬，模型的網(wǎng)格尺度可以設(shè)為大尺度，在格子內(nèi)采用SGS模型來進行更小尺度的模擬。在采用適當?shù)膩喐褡幽Ｐ偷那闆r下，這樣能夠模擬真實物理情況下的瞬時流場。

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

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

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

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

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

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

3.2 并行算法實現(xiàn)

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

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

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

步驟1 設(shè)定好頸動脈計算區(qū)域和最小網(wǎng)格尺寸，讀入面網(wǎng)格文件；

步驟2 判斷血管邊界，遍歷標記計算網(wǎng)格；

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

步驟4 遍歷網(wǎng)格點，進行計算數(shù)據(jù)初始化；

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

步驟6 模型遷移過程；

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

步驟8 若滿足設(shè)置的輸出迭代步數(shù)條件，輸出中間計算結(jié)果；

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

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

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

4 數(shù)值實驗與分析

4.1 實驗環(huán)境

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

4.2 計算模型

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

由于計算的真實血管形狀非常復雜，要精確描述血管壁面區(qū)域的不規(guī)則形狀，需要在適當?shù)奈恢貌捎秒A梯逼近或者插值處理，以保證滿足物理邊界上的條件。因此，邊界條件采用YMS 曲線插值處理［23］，在每次碰撞遷移的過程中，需要在每個與物理邊界相交的網(wǎng)格線上進行插值處理，血管壁面均采用無滑移邊界條件。

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

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

表1 脈動血流參數(shù)Tab.1 Pulsating blood flow parameter

4.3 數(shù)值實驗

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

4.3.1 流速矢量圖

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

4.3.2 剪切力矢量圖

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

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

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

4.3.3 并行時間和加速比

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

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

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

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

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

5 結(jié)語

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

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

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