颅内横窦静脉血流动力学计算机仿真

左 涛，彭红梅，张东威，刘宝治，赵洪明

（1.内蒙古民族大学 数理学院，内蒙古 通辽 028043；2.内蒙古民族大学 附属医院，内蒙古 通辽 028007）

近年来，人口老龄化水平加剧，根据《中国心血管健康与疾病报告2020概要》［1］，中国居民患脑血管疾病人群的数目和死亡率在增加，由此可见，脑血管疾病对人类健康危害越来越大。各发达国家从上世纪末开始极力重视从血流动力学的角度研究脑血管疾病的成因。目前我国也在大力支持和发展基础医学研究，利用计算机模拟血流动力学状况代替血流在体实验也取得了很大的成效，但是多数研究为脑动脉血管血流状况，对于脑静脉血管血流状况的研究较少，因此有必要对其进行研究。

大脑静脉包括：上矢状窦、下矢状窦、直窦、窦汇、横窦、乙状窦、海绵窦、岩上窦及岩下窦。下矢状窦向下连接直窦，直窦与上矢状窦向下汇入窦汇，海绵窦连接岩下窦和岩上窦（图1）。笔者所研究的横窦是颅内最大静脉窦之一，也是颅内静脉窦构造连接的重要枢纽。

静脉横窦又名侧窦或水平窦，位于枕骨内面的横窦沟内。横窦（图1）起于窦汇，延续为乙状窦，接受岩上窦和上矢状窦回流的血液［2］。因其部位存在弧状和乙状2种形状弯曲，其血流环境复杂，临床发现，在此处易发生颅内静脉窦血栓、静脉窦狭窄等脑血管疾病［3］。因此对静脉窦从血流动力学角度观察横窦处的血流状况，以分析出现血管堵塞、狭窄等疾病的原因，从而更好地为临床诊治和医务人员提供有效信息。

图1 脑静脉Fig.1 Cerebral vein

应用计算流体力学的方法，采用内蒙古民族大学附属医院神经内科提供的二维CT数据，利用医学建模软件MIMICS21.0进行血管重建，将得到横窦三维几何模型导入Geomagic Studio中进行切割分离和光滑处理，利用ICEM CFD进行网格划分，将模型结果导入FLUENT2020 R1中进行数值模拟，观察模拟结果显示在不同血液入口速度下血液流线分布、血液压力分布、血液壁面切应力分布［3-4］，通过模拟结果显示的血流动力学特性可以分析患者横窦的血流动力学原因。同时对横窦静脉血流状况进行计算机仿真，可以数值模拟分析出在颅内静脉窦处发生血栓、狭窄的血流动力学原因，为颅内脑血管疾病（如脑卒中）的诊断和治疗提供帮助。

1 材料和方法

1.1 三维模型几何构建

数据获取的图像采集方法为计算机断层扫描技术（CTA），图像数据参数为：平面分辨率为512×512，像素大小为0.488 3 mm，层间距为0.900 mm，共256张切片。医学影像三维重建技术是指将二维影像数据转换成三维可视化图像的技术，方便研究人员和医务人员对图像进行分析［4］。将CT数据（DICOM格式）导入Mimics后，首先阈值提取轮廓：选取合适的阈值范围使得轮廓最清晰，阈值范围为267~2 378，形成蒙面后进行3D计算得到骨骼、血管及组织模型（图2（a）），通过分离、分割去除骨骼和组织，从模型中提取出横窦静脉（图2（b）），对模型进行边缘切割、去除冗余结构，在3-matic中对其进行后处理，获得较为光滑的3D血管模型（图2（c））。

图2 血管模型Fig.2 Vascular models

1.2 计算模型与方法

1.2.1 计算模型网格划分

网格划分是计算流体力学的基础，网格的质量直接影响后面在fluent软件的数值模拟。选用在ICEM CFD软件中设定边界层网格，生成面网格。图（3）为横窦血管的网格划分模型图，在模型中共生成了553 975个网格和205 830个节点。

1.2.2 控制方程与边界条件

将横窦静脉中的血液假设为不可压缩牛顿流体，并且血液满足绝热、均匀等条件，血液流动满足的（Navier-Stokes）偏微分方程如下［5］：

图3 网格划分模型Fig.3 Grid division model

其中，t为时间，U为流体的速度矢量，P为流场压强，血液密度ρ=1 050 kg·m-3，血液黏性系数μ=0.003 5 Pa·s。

设定计算模型为稳态流动，颅内静脉窦的血流入口速度设为0.12 m·s-1和0.22 m·s-1［6］。血管的入口直径d=8.56 mm，根据雷诺数公式：，在血液入口速度为0.12 m·s-1时，血液的雷诺数为313，在血液入口速度为0.22 m·s-1时，血液的雷诺数为565，雷诺数Re<2 300，判断血液流动为层流［4］。在计算过程中，设置出口压强为0 Pa，假设血管壁刚性无滑移，不考虑重力的影响。

1.2.3 计算数值模拟

使用计算流体力学软件Fluent采用有限体积法进行数值模拟，有限体积法是将模型区域划分为一系列控制区域，每个控制区域选一个节点作为代表，通过守恒型的控制方程，对控制体积做积分导出离散方程。

计算收敛条件时，将残差设置为10-5，最大迭代次数设置为200。分别设置血流入口速度为0.12 m·s-1和0.22 m·s-1，计算颅内横窦静脉的血液流场分布、壁面切应力及壁面压力分布等血流动力学特性。同时可以将Fluent软件数值模拟的结果代入CFD-post软件进行可视化处理，得到横窦血管的血液流场分布、壁面压力及壁面切应力等云图。

2 计算结果

2.1 血液流场分布

本研究选用健康成年人脑CT数据，得到的结果为正常血液流动状态分布，流场的空间分布见图4。图4（a）是入口速度为0.12 m·s-1血液在静脉横窦中血液流场分布图，除血管弯曲区域血液流动模式为直线流动，在横窦静脉乙状弯曲处（箭头1所指）出现明显涡旋流动现象，并且血液流速很低，血流方向错综复杂，在横窦静脉弧状弯曲处（箭头2所指）其流动模式为螺旋流动，血液流速缓慢且流动不规律，箭头（1、2）所指处均出现低流速区，血液流速在此区域存在最小值。弧状弯曲血液流入端（箭头3所指）血液流速存在最大值。图4（b）是入口速度为0.22 m·s-1血液在静脉横窦中血液流场分布图，在静脉横窦中血液流动模式大部分为层流流动，在横窦静脉乙状弯曲处（箭头1所指）血液流动不规律并伴有紊流现象，但血液流速明显加快，在横窦静脉弧状弯曲处（箭头2所指）血液流速增加且流动不规律，血管弯曲处（箭头1、箭头2处）低流速区域减小，但仍为低流速区域，在此区域存在最小速度，弧状弯曲血液流入端（箭头3所指）血液流速存在最大值。

图4 血液流场分布Fig.4 Blood flow field distribution

2.2 血管壁面压力分布

本例CT数据模拟显示的颅内横窦静脉血管壁面压力分布图见图5。图5（a）是入口速度为0.12 m·s-1时血管壁面压力分布图，血液从模型入口到出口壁面压力逐渐递减，在壁面压力云图中各个压力值范围分布规则，压力差（图中最大压力与最小压力的差值）范围为0~81.250 Pa，在乙状弯曲处（箭头1所指）压力值较小，属于低压力区，在弧状弯曲处（箭头2所指）虽然血管壁面压力值相对于乙状弯曲处较大，但依然属于低压力区。图5（b）是入口速度为0.22 m·s-1时的血管壁面压力分布图，在云图中各压力值分布范围规则，压力差范围为0~130.000 Pa，从模型中可以看出血管壁面压力从入口到出口逐级递减，在乙状弯曲处（箭头1所指）增大血液入口速度，血管壁面压力明显增大，但仍属于低压力区，在弧状弯曲处（箭头2所指）增大血流入口速度，血管壁面压力显著增大，属于高压力区。如图5（b）所示，由于血管弯曲处血液流动不规律并伴有紊流现象，血液呈三维螺旋流动，相对于附近直线流动血管壁面压力较大。

图5 血管壁面压力分布图Fig.5 Distribution diagram of vascular wall pressure

2.3 壁面切应力分布

血管壁面切应力又称内皮剪切应力，是腔内流动的黏性血流对其流经的血管壁表面产生的摩擦力［7］。本例CT数据模拟颅内横窦静脉血管壁面切应力分布见图6。图6（a）为入口速度为0.12 m·s-1时血管壁面切应力分布图，血管乙状弯曲（箭头1所指）壁面切应力很小，属于低切应力区域，弧状弯曲（箭头2所指）的内侧和外侧壁面切应力都很小，都属于低切应力区域。图6（b）为入口速度为0.22 m·s-1时血管壁面切应力分布图，横窦静脉乙状弯曲（箭头1所指处）内外侧壁面切应力明显增大，属于高壁面切应力区域。弧状弯曲（箭头2所指处）外侧壁面切应力变化明显，存在高切应力区域，内侧壁面切应力变化很小，属于低切应力区域，外侧壁面切应力明显大于内侧。

图6 血管壁面切应力分布Fig.6 Distribution of vascular wall shear stress

3 讨论

通过对比图4（a）入口速度为0.12 m·s-1和图4（b）入口速度为0.22 m·s-1的血液流场数值模拟结果可得，弧状弯曲和乙状弯曲区域血液流动环境复杂，图4（a）中弧状弯曲和乙状弯曲位置（箭头1、箭头2所指）血液流速很低，均为低流速区域，图4（b）相对于图4（a）血液入口速度增大，血液在此位置（箭头1、箭头2所指）血液流速虽然有增大的现象，但是仍然属于低流速区域。流速持续很低，血液中血小板和脂类等成分附壁堆积，易造成血管狭窄，影响血液流动［8］。当增大血液入口速度时，在弧状弯曲和乙状弯曲位置（箭头1、箭头2所指）血液流速明显增大，但仍有低流速区域。

通过对比图5（a）入口速度为0.12 m·s-1和图5（b）入口速度为0.22 m·s-1的血管壁面压力数值模拟结果可得，从模型入口到出口血管壁面压力逐渐递减，在图5（a）中弧状弯曲和乙状弯曲处压力值较小，属于低压力区域，图5（b）相对于图5（a）血液入口速度增大，弧状弯曲和乙状弯曲虽然血管壁面压力明显增大，但乙状弯曲（箭头1所指）属于低压力区域，弧状弯曲（箭头2所指）属于高压力区域。高压力环境（图5（b）中箭头2所指）会导致血管结构改变，损伤壁面［8-9］。图5（a）血管壁面压力差范围为0~81.250 Pa，图5（b）血管壁面压力差范围为0~130.000 Pa，图5（b）相对于图5（a）血流入口速度增大，血管壁面压力差范围在变大。

通过对比图6（a）入口速度为0.12 m·s-1和图6（b）入口速度为0.22 m·s-1的血管壁面切应力数值模拟结果可得，血液入口速度为0.12 m·s-1和0.22 m·s-1的分布规律相似，从切应力空间分布来看，切应力分布与血液流速和血管管径有关。图6（a）弧状弯曲（箭头2所指）和乙状弯曲（箭头1所指）内外侧的壁面切应力值都很小，均为低切应力区域，图6（b）相对于图6（a）血液入口速度增大，乙状弯曲（箭头1所指）内侧切应力明显大于外侧切应力，而弧状弯曲（箭头2所指）外侧切应力明显大于内侧切应力，随着血液入口速度的增大，弧状弯曲和乙状弯曲区域产生高切应力区域，横窦静脉持续存在高壁面切应力区域会增大血液与血管的摩擦，损伤血管内膜。横窦静脉持续存在低壁面切应力区域，血管内脂类颗粒及其他物质容易滞留、堆积，造成血管局部狭窄，易发生病变。

4 结论

本研究通过医学影像软件MIMICS21.0基于CT数据成功建立了颅内横窦的三维血管模型，将得到的三维血管模型代入到计算流体力学软件FLUENT中进行血流动力学数值模拟，同时计算出在入口速度为0.12 m·s-1和0.22 m·s-1时的横窦静脉的血液流场、壁面压力、壁面切应力等血流动力学特征，数值模拟分析显示在静脉横窦弧状弯曲和乙状弯曲处血流速度缓慢，出现涡旋流动等血流现象，在此区域壁面切应力值较小，为低切应力区，此现象说明了在静脉横窦弯曲处容易发生淤滞，血管内流通区域减少，影响血液正常流动，增大血液入口速度后，弧状弯曲处出现高压力区域，高压力环境会导致血管结构改变，损伤血管壁面，在乙状弯曲内侧和弧状弯曲外侧出现高壁面切应力区域，高壁面切应力容易使血管受损。横窦由于其构造的独特性，在弯曲区域血流环境复杂，极易发生疾病，通过计算机仿真模拟真实地反映了横窦静脉血管内血液流动，可以帮助医务人员更清楚地认识横窦静脉的血流动力学特点，通过血流动力学参数对患病风险进行评估，为颅内横窦疾病的治疗提供理论支持。