基于VLBM的髂外动脉血流动力学分析

2022-03-03戎晨彬严微微俞慧丹许友生

浙江科技学院学报 2022年1期

戎晨彬,陈柔,严微微,俞慧丹,许友生

(1.浙江科技学院机械与能源工程学院,杭州 310023;2.中国计量大学计量测量工程学院,杭州 310018;3.印第安纳大学-普渡大学印第安纳波利斯分校机械工程系,印第安纳波利斯 46202)

下肢动脉粥样硬化狭窄是一种常见的动脉狭窄疾病,其主要特征是由下肢动脉内壁上形成沉积[1](斑块)而导致的动脉弹性下降,管径狭窄。临床表现早期有间歇性跛行、下肢发凉麻木和静息性下肢疼痛等,晚期则伴随着动脉狭窄可能会出现发绀、溃疡和坏疽等[2-3]。该病的发病率逐年增高,且随年龄的增长而增长,据调查[4-5],45岁以上人群的发病率约为5%,而70岁以上人群的发病率则大于15%。

髂外动脉是髂总动脉的终支之一,向下肢供应血液。准确的动脉血压测量对判断动脉粥样硬化狭窄程度十分重要,如脉搏波传导速度,该指标是测量动脉硬化的可靠方法,但测量数值与血压(舒张压和收缩压)的相关性较强[6-7]。有创动脉血压(invasive blood pressure,IBP)是监测探头经动脉穿刺置管后直接测量动脉血压,是测量血压的金标准,测量值准确可靠;但并发症较多,容易造成血肿、血栓、感染、血管损伤等,并且对穿刺技术要求极高[8]。计算机断层扫描血管造影(computed tomography angiography,CTA)是一种临床常用的无创动脉诊断方法,但它无法提供相关的流体力学信息。近年来,随计算机硬件水平的快速发展,CTA与计算流体力学(computational fluid dynamics,CFD)相结合来诊断动脉疾病的技术也愈发成熟,尤其在冠状动脉领域[9-10]已取得一定的成效,而在髂动脉领域鲜有相关研究。作为基于介观模拟尺度的CFD方法,格子玻尔兹曼方法(lattice Boltzmann method,LBM)具有一些独特的优势,例如易处理复杂的几何形状和边界[11-12],非常适合计算多相流等复杂流体[13-14],以及易实现大规模的图形处理器(graphics processing unit,GPU)并行计算[15]。体积格子玻尔兹曼方法(volumetric lattice Boltzmann method,VLBM)是Yu等[16]基于LBM提出的,相比流体粒子位于格点的传统LBM,VLBM假设粒子均匀分布在格子单元内,随碰撞将迁移至相邻格子单元,并将格子单元按流体、固体和边界进行分类。该方法可以更好地处理任意弯曲的复杂边界,并提高了GPU并行计算能力。

在上述研究的基础上,本研究提出一种基于GPU并行架构开发的VLBM,可在短时间内获得血流动力学分析结果,从而辅助医生临床诊断。为验证此方法的可行性,我们对狭窄髂外动脉的血流动力学进行了研究,首先介绍VLBM及其在模拟中的应用;然后对比模拟结果和IBP结果,验证VLBM在狭窄髂外动脉血流动力学分析中的可行性,并通过将支架植入狭窄处进行修复,获得恢复动脉的血流动力学信息;最后总结分析狭窄对血流的影响。本研究中狭窄髂外动脉的CTA影像、IBP测量结果等数据均由杭州市第一人民医院于2017年提供。

1 原理与方法

1.1 体积格子玻尔兹曼方法

体积格子玻尔兹曼方法(VLBM)假设粒子均匀分布在格子中,并将格子单元分成流体、固体和边界3种类型,如图1所示。引入参数P(x,t)=Vs(x,t)/V定义每个格子单元中固体体积的百分比,其中,P(x,t)=0的格子单元为流体单元,P(x,t)=1的为固态单元,以及0

图1 VLBM中格子单元的分类

我们以Qian等[17]于1992年提出的DdQm(d维空间,m个离散速度)系列LBM基础模型中的D2Q9模型为例,来讨论传统LBM和VLBM的差异。图2(a)为传统LBM模型,粒子在格点处碰撞后,沿其速度方向往相邻格点迁移,fi(x,t)为粒子密度分布函数。而在VLBM中,假设粒子均匀分布在格子单元内,随碰撞迁移至相邻格子单元中,如图2(b)所示。与传统LBM的fi(x,t)类似,VLBM引入ni(x,t)表示离散速度为ei的粒子在时间t占据格子单元的分布函数,i为离散速度数量,fi(x,t)与ni(x,t)的关系为

图2 传统LBM和VLBM粒子分布示意图

ni(x,t)=fi(x,t)[1-P(x,t)]V。

(1)

式(1)中:V取单位体积(V=1)。与LBM相似,VLBM的碰撞方程如下:

(2)

式(2)中:n′i(x,t)为碰撞后的偏分布函数,其中i=0,1,2,…,b;τ为无量纲弛豫时间;ni,eq(x,t)为平衡态分布函数,

(3)

式(3)中:ωi为权重系数;cs为格子声速;U为宏观速度;N(x,t)=∑ni(x,t)为当前格子单元密度。本研究选用i=0,1,2,…,18的D3Q19格子模型进行三维模拟,离散速度ei在i=0时为(0,0,0),i=1～6时为(±1,0,0)、(0,±1,0)、(0,0,±1),以及i=7～18时为(±1,±1,0)、(0,±1,±1)、(±1,0,±1)。相对应的权重系数在i=0时ωi=1/3,i=1～6时ωi=1/18,以及i=7～18时ωi=1/36。为将非滑移边界条件整合到格子模型中,流动部分分为两部分:1)从迎风方向的邻近格子单元中迁移出的粒子数[1-P(x,t)]n′i(x-eiδt,t);2)从逆风方向的格子单元反弹的粒子数P(x-eiδt,t)n′i*(x,t)。此处δt为时间步长,n′i*(x,t)为对应n′i(x,t)从非滑移壁面反弹后的粒子分布函数。因此,迁移方程为

n″i(x,t+δt)=[1-P(x,t)]n′i(x-eiδt,t)+P(x-eiδt,t)n′i*(x,t)。

(4)

宏观参数密度和速度分别从粒子分布函数的0、1阶中获得:

ρ(x,t)=∑ni(x,t)/[1-P(x,t)];

(5)

u(x,t)=∑eini(x,t)/∑ni(x,t)。

(6)

压强可从密度中获得:

(7)

式(7)中:p0为参考压强;ρ0为参考密度。

1.2 计算模型设置

通过对CTA影像进行处理,提取狭窄髂外动脉形态结构模型。CTA是一种医学检查方法,它将CT扫描与注射造影剂相结合,产生人体特定部位的多层等距二维数字图像,图像由大量不同灰度值的像素按矩阵排列方式构成。3D Slicer软件可以设定灰度阈值从而分割出二维图像中髂动脉部分,提取该部分的轮廓线,并通过几何单元拼接来拟合动脉壁,最终达到重构3D动脉模型的效果。重构支架植入前的3D髂外动脉狭窄模型如图3(a)所示,该髂外动脉为重度狭窄,狭窄率约为70%。狭窄率表示动脉的狭窄程度,常用的计算公式为:狭窄率=(狭窄远端正常直径-狭窄段的最窄直径)/狭窄远端正常直径。支架植入后如图3(b)所示,髂动脉模型通过3D Slicer软件拉伸狭窄处壁面后导出,将狭窄处恢复成健康动脉。

图3 3D髂外动脉结构模型

为简便起见,假设血液为不可压缩牛顿流体,密度为1 060 kg/m3,运动黏度为3.3×10-6m2/s,动脉壁为刚性壁,且壁面采用无滑移条件。将脉动流作为髂总动脉剖面的速度入口边界条件,在髂内外动脉末端设置压力出口边界条件,具体边界条件方案如下:

1)速度入口边界条件。多普勒超声(Doppler ultrasound,DUS)可测量一个心动周期内的血流速度曲线。在模拟中,需将速度曲线插值成1 000个时间点,如图4(a)中的蓝点所示。本研究选取了6个具有代表性的红点作为采样时间点,以表示血流在各个阶段(舒张期和收缩期)流动时的不同现象。在指定时刻对髂总动脉剖面设定速度边界时,将DUS同一时刻测得的速度作为髂总动脉中心处的最大速度,然后用抛物线分布到壁面上,以此贴近实际情况,如图4(b)所示。值得注意的是,速度值是脉动的,每个时间点的入口边界应按上述方法设置。为将速度边界条件引入VLBM,我们采用Guo等[18]于2002年提出的非平衡态外推格式,该方法具有运算简单、数值稳定等优点。

图4 髂总动脉入口血流速度

2)压力出口边界条件。压力出口使用与电路相似的三元windkessel模型[19],图5为该模型应用于动脉中时的示意图。血流速度(Q)由近端血管阻力(r)、远端血管顺应性(C)和每个出口下游远端血管阻力(R)的综合效应控制,因此,三元windkessel模型也被称为RCR模型。基于血流和电路之间的相似性,流速Q和平均压力p之间的关系满足以下常微分方程[20]:

图5 动脉中三元windkessel模型示意图

(8)

求解式(8)可得到压力的解析解

(9)

式(9)中:p0为出口的初始压力。R、C和r这3个参数的关系和取值范围可通过参考文献[21]来确定。同样,采用非平衡态外推格式将压力边界条件引入VLBM。

1.3 算法流程设计

使用C++编程语言开发基于GPU并行架构的VLBM算法,算法的设计流程具体可分为以下几个步骤:

1)读取髂外动脉模型,根据模型尺寸确定计算区域;

2)设置格子单元长度l,对计算区域进行格子划分,计算每个格子单元中固体体积的百分比P(x,t),并设置边界条件及给定所有格点初始宏观参量;

3)使用引入VLBM的D3Q19模型的控制方程,通过给定的初始宏观参量计算出所有格点各个速度方向的平衡态分布函数,作为模拟的初场;

4)通过粒子碰撞迁移规则求解控制方程;

5)使用非平衡态外推格式对相应边界处格点进行边界处理;

6)计算所有格点的宏观参量,如速度、密度、压力等;

7)判断宏观参量是否收敛;

8)如果宏观参量收敛,输出计算结果;否则重回第4步继续计算,直至计算结果收敛。

使用该算法模拟狭窄髂外动脉,将狭窄处上下端的血压设置为计算输出结果进行格子无关性检验,发现在达到格子单元长度l=2.99×10-4m后,计算输出结果的变化很小,因此可用该格子单元长度模拟狭窄髂外动脉,此时对应的格子数量为142×198×258。对于一个典型的髂外动脉模拟任务,在计算机CPU为Intel Xeon E5-2683 v3,GPU为NVIDIA Tesla P100的配置下,需要约20 min的时间可完成模拟任务。

2 结果与讨论

2.1 验证VLBM的可行性

验证VLBM用于狭窄髂外动脉模拟的可行性是第一步,也是至关重要的一步。我们将模拟获得的狭窄处上下端血压与IBP结果做对比以进行验证,IBP结果来自临床测量,如图6的红色虚线所示。在使用计算结果之前,需要确认时间收敛性,当宏观参数为每个相邻的心动周期收敛时,设置收敛条件。图6中的蓝色实线为使用自编的VLBM算法计算的髂外动脉狭窄处上下端血压曲线。从对比结果可以看出,2条代表模拟血压和IBP的曲线具有较高的吻合度。定量地检测血压的具体参数,即收缩压Ps、舒张压Pd和平均血压Pmean。收缩压是血压曲线中的最大值,舒张压是血压曲线中的最小值,平均血压定义为Pmean=(Ps+2Pd)/3[22]。表1为髂外动脉狭窄处上下端IBP与VLBM血压参数对比及相对误差,其中相对误差为0.455%～1.585%。因此,无论是定性还是定量的数据结果均表明将VLBM用于狭窄髂外动脉的血流动力学分析是可行的。

图6 髂外动脉狭窄处上下端IBP与VLBM血压对比

表1 髂外动脉狭窄处上下端IBP与VLBM血压参数对比及相对误差

2.2 支架植入前后血流动力学分析

分别选取支架植入前狭窄率为70%和支架植入后还原的髂外动脉流场,观察狭窄对流场的影响。图4(a)的6个抽样点中,t=0.20 s时的入口血流流速最大,狭窄处上下端的流场变化最明显,可以更好地显示流场的变化。因此以t=0.20 s为例,图7显示了支架植入前后髂外动脉的流线,将狭窄处放大后做支架植入前后的速度截面图如图8所示,从中可以看出,与支架植入后还原的髂外动脉相比,血流流经狭窄处后,流速明显增大并向动脉中心处汇聚。与髂外动脉狭窄处血流速度有明显的增加相比,髂内动脉血流速度的变化可以忽略不计。

图7 髂外动脉支架植入前后的流线图

图8 髂外动脉支架植入前后的血流速度截面图

壁面切应力(wall shear stress,WSS)是血流作用于动脉内壁产生的机械力之一[23],对动脉粥样硬化狭窄的形成起着重要作用。有研究表明[24-27],在早期低WSS血流分布区域更易于动脉粥样硬化斑块的形成,这可能导致动脉管腔逐渐变窄,斑块侵入腔内,WSS会先恢复正常,之后开始升高,WSS越高,对动脉壁施加的拉伸力越大,从而增大动脉破裂的风险。将图7中髂外动脉狭窄(还原)部分放大后如图9和图10所示。图9(a)显示了狭窄髂外动脉的WSS分布,其中时间点取自图4(a)中的6个采样时间点。狭窄部位WSS较大,高WSS的分布范围随血流速度的变化而变化。当t=0.20 s血流接近收缩期峰值时,高WSS的分布范围更广,随后逐渐减小,与一个心动周期内血流的变化趋势一致。支架植入后,每个时间点的WSS均减小。综上所述,狭窄将导致高WSS,增加动脉破裂的风险。

图9 髂外动脉支架植入前后的WSS分布

图10 髂外动脉支架植入前后的压力云图

图10为髂外动脉支架植入前后的压力云图。支架植入后血压随时间变化,还原处上下端压降约为0.047 6 mmHg,可以忽略不计。相反,支架植入前的狭窄处有一个明显的压降,此外,血流量越大,压降越大。经计算,支架植入前狭窄处的平均压降约为8.838 5 mmHg。从定量角度来分析,如图11所示,支架植入后,收缩压在狭窄(还原)处上端降低,同时在下端升高,以此减小压降,而支架植入前后舒张压变化不大。因此,狭窄将导致压降增大,尤其是在流速较高的收缩期。