基于计算流体力学的载瘤血管模型流体流动特性及数值仿真分析

Fluid Flow Characteristics and Numerical Simulation Analysis of Aneurysm-Bearing Vessel Models Based on Computational Fluid Dynamics

摘要: 颅内动脉瘤(IA)是脑动脉局部血管壁因结构异常或力学环境失衡，发生局限性异常扩张所形成的病理性病变，多数患者发病初期无明显症状，但破裂后易引发蛛网膜下腔出血等致命性并发症，致死率与致残率较高。深入探究IA的血流动力学机制，对破裂风险评估及干预策略制定具有重要意义。本研究采用计算流体力学(CFD)技术，基于载瘤血管几何模型，通过有限体积法(FVM)离散求解Navier-Stokes方程，在生理脉动血流条件下，对动脉瘤腔内的流体流动特性开展数值模拟研究，通过分析血流再循环、流速场分布规律、壁面压力及壁面切应力(WSS)等关键血流动力学参数，探究瘤腔内血流紊乱现象的诱发机制及其关键影响因素。数值仿真结果表明，动脉瘤瘤腔内存在显著血流回流区，流动分离效应与复杂涡流结构的形成进一步加剧局部血流不稳定性；低WSS区域与血流停滞区的空间分布与动脉瘤易破裂部位呈现高度相关性。

Abstract: Intracranial aneurysm (IA) is a pathological lesion characterized by localized abnormal dilation of the intracranial arterial wall, resulting from structural abnormalities or mechanical homeostasis imbalance. Most patients exhibit no obvious symptoms in the early stage; however, rupture of IA tends to induce life-threatening complications such as subarachnoid hemorrhage, leading to high mortality and morbidity rates. In-depth investigation into the hemodynamic mechanisms of IA holds significant implications for rupture risk assessment and the formulation of intervention strategies. This study employs computational fluid dynamics (CFD) technology and is based on the geometric model of tumor-bearing vessels. By discretizing and solving the Navier-Stokes equations using the finite volume method (FVM), numerical simulations of the hemodynamic characteristics within the aneurysm sac were conducted under physiological pulsatile flow conditions. Key hemodynamic parameters including blood flow recirculation, velocity field distribution patterns, wall pressure, and wall shear stress (WSS) were analyzed to explore the inducing mechanisms and critical influencing factors of flow disturbance in the aneurysm sac. The numerical simulation results demonstrate that there are significant blood flow recirculation zones within the aneurysm sac. The flow separation effect and the formation of complex vortex structures further exacerbate local flow instability. The spatial distribution of low WSS regions and blood flow stasis zones shows a high correlation with the prone-to-rupture sites of the aneurysm.

1. 引言

颅内动脉瘤(IA)的好发部位主要集中在血管分叉区域(终末动脉瘤)与动脉侧壁位置(侧壁动脉瘤)。对患者血管系统进行血流动力学分析，对了解和诊断其疾病具有重要价值[1]-[3]，如血管狭窄[4]、各类动脉瘤[5] (如梭形、囊状、假性及夹层动脉瘤)、血管瘤等病变。计算流体力学(CFD)模拟技术目前已在血流动力学分析领域得到广泛推广与应用，能够对流动结构、再循环流、壁面压力及壁面切应力(WSS)等关键参数展开相关研究，其计算域模型主要分为两类互补类型：一类是基准简化模型[6]，另一类是基于患者特异性血管系统构建的真实模型[7]。鉴于血流动力学参数可以为临床治疗策略评估提供重要参考，CFD已成为动脉瘤血流分析的主要技术方法，为相关研究提供可靠的数据支持。例如，Cebral等人通过数字减影血管造影(DSA)重建真实血管通道以显示颅内囊状动脉瘤内的血流结构，在15例脑动脉瘤患者的研究中，CFD模拟结果与DSA检查结果的一致性达78% [8]；Xiang等学者研究了血流动力学参数在动脉瘤破裂风险评估中的应用价值，结果表明基于真实血管数据的CFD模拟可有效预测动脉瘤的破裂风险[9]。Bauer等人针对梭形动脉瘤，采用脉动流与过渡流两种血流方案分析壁面切应力(WSS)分布特征，研究结合CFD模拟与磁共振测速(MRV)、激光多普勒测速(LDV)等实验方法，证实CFD与LDV检测结果高度吻合，同时发现涡结构的形成对WSS的时空分布具有显著影响，WSS最大值出现于涡初始生成区域，且动脉瘤未扩张区域及流动过渡区(层流–湍流转换区)表现为层流流动[10]。

探究血流参数对血管疾病血流动力学特性的影响，对临床研究与诊疗实践具有重要意义[11] [12]。Iimuro等人针对1例75岁女性门静脉系统动脉瘤患者在血液流速1000 mL/min的工况下进行仿真计算。结果显示，该患者动脉瘤腔内出现湍流状态的再循环流，且瘤壁区域的壁面切应力(WSS)显著升高；此类异常血流动力学参数可能促进血栓形成或进一步增加动脉瘤破裂风险[13]。

现有统计结果表明，全球50岁左右人群中未破裂颅内动脉瘤的患病率约为3%；而在中国35~75岁人群中，该疾病患病率达7% [14]，发病率介于2%~7%之间，破裂病例的死亡率约为60% [15]。尽管近年来不同地区的死亡率呈现下降趋势，但多数存活患者仍会遗留显著的神经功能障碍，这一现状给患者家庭及社会带来了沉重负担。Souza等人针对颅内动脉瘤模型开展了数值研究，重点分析血流特性、结构响应及血管壁力学行为等关键参数，这些参数为动脉瘤破裂风险的量化评估提供了关键参考。该研究采用四种不同雷诺数(Re = 1, 100, 500, 1000)对同一模型进行CFD模拟。结果显示，随着雷诺数增大，动脉瘤区域逐渐出现血流再循环现象；且壁面切应力(WSS)、血管壁位移及应变均在最高雷诺数(Re = 1000)条件下达到最大值[16]。Shishir等人针对弯曲动脉段伴囊状动脉瘤的血流动力学特性展开研究，构建了不同动脉瘤直径(5~8 mm)的简化几何模型并进行CFD模拟。模拟设定参数如下：雷诺数(Re)为110~850，血液密度1008 kg/m⁻³，体积流量3.6 mL/s，入口流速0.2365 m/s，动力粘度0.0035 Pa∙s。CFD模拟结果表明，动脉瘤区域出现的再循环流(涡流)特性受动脉瘤尺寸的显著影响[17]。

现有研究在几何建模与流动条件设定上各有侧重：或采用简化模型与定常流条件探讨普遍规律，或基于真实模型但聚焦于特定时相的稳态分析。本研究旨在整合上述方法的优势，基于患者特异性影像数据构建的个体化终末动脉瘤几何模型，施加符合生理条件的入口脉动流速波形进行瞬态仿真模拟。系统分析在一个完整心动周期内，血流再循环、壁面压力、壁面切应力(WSS)等血流动力学关键参数的动态演变过程及其空间分布特征，揭示流动结构与动脉瘤壁面力学环境之间的关联，为评估动脉瘤的稳定性和破裂风险提供更具时变分辨性与个体针对性的血流动力学依据。

2. 数值方法

2.1. 控制方程和边界条件设定

血液运动由瞬态条件下的Navier-Stokes方程控制，其数学表达式如下：

$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u_x\right)}{\partial x}+\frac{\partial \left(\rho u_y\right)}{\partial y}+\frac{\partial \left(\rho u_z\right)}{\partial z}=0$ (1)

$\{\begin{array}{l}\frac{\partial \left(\rho u_x\right)}{\partial t}+\nabla \cdot \left(\rho u_x\overrightarrow{u}\right)=-\frac{\partial p}{\partial x}+\frac{\partial {\tau }_{xx}}{\partial x}+\frac{\partial {\tau }_{yx}}{\partial y}+\frac{\partial {\tau }_{zx}}{\partial z}+F_x\\ \frac{\partial \left(\rho u_y\right)}{\partial t}+\nabla \cdot \left(\rho u_y\overrightarrow{u}\right)=-\frac{\partial p}{\partial y}+\frac{\partial {\tau }_{xy}}{\partial x}+\frac{\partial {\tau }_{yy}}{\partial y}+\frac{\partial {\tau }_{zy}}{\partial z}+F_y\\ \frac{\partial \left(\rho u_z\right)}{\partial t}+\nabla \cdot \left(\rho u_z\overrightarrow{u}\right)=-\frac{\partial p}{\partial z}+\frac{\partial {\tau }_{xz}}{\partial x}+\frac{\partial {\tau }_{yz}}{\partial y}+\frac{\partial {\tau }_{zz}}{\partial z}+F_z\end{array}$ (2)

式中，x、y、z三个方向的速度分量分别以u_x、u_y、u_z表示，t代表时间变量，$\rho$ 为流体密度(单位：kg/m³)，p表示流体单元所受压强，$\tau$ 为粘性应力，${\tau }_{ij}$ 中，i代表正应力的作用方向，j为作用面对外法线方向，F为体积力，下标字母表示微元体体积力方向。

已有研究结果[18]证实，当血管内径大于0.5 mm时，采用牛顿流体替代非牛顿流体开展数值模拟，所得误差能够控制在2%范围之内。本文研究的颅内动脉属于大动脉，血管内径远大于0.5 mm。因此本研究将血液假设为不可压缩牛顿流体，其密度取值1060 kg/m³，动力粘度为0.0035 Pa∙s；动脉壁边界设定为刚性和非滑移边界条件，流入口采用随心动周期变化的血流速度边界条件(具体变化规律如图1所示)，出口边界条件设定为压力边界，对应的静压力取值为13,325 pa ≈ 100 mmgh。

Figure 1. Unsteady boundary conditions

图1. 非稳态边界条件

2.2. 计算求解

在非稳态数值迭代计算过程中，设定单个心动周期时长为0.8 s，将其划分为80个时间步数，时间步长$\Delta t=0.01\text{s}$ ，为保障计算收敛性，每个时间步长内执行50次内部迭代；为降低初始条件引发的计算误差，本研究共进行两个心动周期的数值模拟计算，最终选取第二个心动周期的计算结果作为后续的分析依据，基于该周期所提取的参数开展血流动力学特性分析。

2.3. 网格划分

本研究使用Fluent Meshing (ANSYS公司，美国)生成非结构化网格模型。鉴于壁面相关参数在血液流体动力学研究中的重要性，本研究针对网格模型采用混合划分方案：流体域中心区域选用六面体网格，近壁面边界层则划分为多面体网格。在进行网格划分前，先完成血管入口、出口及壁面的命名定义；为平衡网格数量与网格质量，设定多面体单元的尺寸范围为0.1~0.9 mm (最小尺寸0.1 mm，最大尺寸0.9 mm)；为保障数值计算精度，将边界层设置为5层，网格增长率取1.2。网格数为62,071，节点数为180,964。划分好的体网格如图2所示。

3. 结果与讨论

为揭示单个心动周期内不同时刻各物理量的分布规律及动态变化特征，本研究基于入口平均速度脉动曲线，选取四个典型时刻展开分析：0.88 s (心缩期上升阶段)、0.93 s (心缩期峰值时刻)、0.98 s (心缩期下降阶段)及1.5 s (舒张期后期)。其中，时间为0.93 s时，血管内血流速度达到心动周期的最高值。同时，为直观观测各物理量的动态演化过程，本研究在关键区域设置六个观测点(如图3所示)：观测点1设于入口侧血管，观测点2、点3位于瘤顶区，观测点4置于瘤体侧面，观测点5、点6分别设在瘤颈区，后续

Figure 2. Mesh model

图2. 网格模型

Figure 3. Distribution positions of each observation point on the surface of tumor-bearing vessels

图3. 各观测点在载瘤血管表面的分布位置

分析将基于这些观测点的监测数据展开分析。

3.1. 流体流动状态分析

动脉瘤形成后，会导致再循环流(涡流)形成，对WSS、壁压和出口流态产生很大影响。在图4中，心缩期上升阶段瘤腔内的流线稀少，血液流动较为紊乱，此时动脉瘤内的血液并未充满瘤腔，随着血液流速达到最大；当达到心缩期峰值时刻(0.93 s)时，血流充分填充动脉瘤腔，腔内流线呈螺旋状分布；从瘤体左侧观测，可发现顺时针漩涡结构，表明瘤腔内存在局部涡流现象。进入心缩期下降阶段后，入口血流速度逐渐降低，动脉瘤内最大流速同步下降，但受血液惯性效应影响，瘤腔内仍保持较高的流动速度，且在瘤腔血液入口侧形成高速流动带。该高速流体与瘤体中部的低速流体相互作用，导致腔内流场呈现复杂紊乱特征，这一流动状态可能会增加动脉瘤的破裂风险。至舒张期后期(1.5 s)，随着整体血流速度的降低，血液对瘤壁的冲击作用减弱，腔内流线数量也随之减少，但瘤腔内仍存在低流速涡流区域。此外，可以看出在远离动脉瘤，即血管入口和出口处的流动为层流。

(a) 心缩期上升阶段 (b) 心脏收缩期峰值点 (c) 心缩期下降阶段 (d) 舒张期后期

Figure 4. Velocity streamline diagram of tumor-bearing vessels over one cardiac cycle

图4. 载瘤血管在一个心动周期内的流速流线图

在图5(a)、图5(b)中，瘤腔中的两个平面分别表示了心动周期内峰值流速时刻垂直入口血流平面和瘤颈平切处平面的速度幅度等值线。在垂直入口血流平面上，最大速度出现在载瘤血管主干内，血流进入动脉瘤腔后流速明显减小，且靠近瘤壁的区域流速进一步降至0~0.2 m/s，符合壁面黏性作用下的流速边界特征，同时瘤腔内血流呈滞缓状态，为局部涡流的形成提供了流体力学条件，而瘤颈与瘤腔之间2~4倍的流速差异，体现了该区域血流动力学的显著异质性。在图5(c)中，速度分布矢量显示了流动与瘤壁的相互作用，导致在瘤顶方向上形成涡流的反向流动。在瘤颈平切截面处，血流呈现显著的区域性分布特征：一侧为高速血流区域，另一侧则为低速血流区域，这是由于在流动–壁面作用区附近存在回流流动，如图5(b)、图5(d)所示。随着血液流动速度的增加，瘤体入口侧形成高速血流集中区，且该区域的速度梯度显著；在此过程中，血流对动脉瘤壁的冲击区域也随之扩大。

图6为各观测点的速度曲线图。如图所示，0.93 s时，除了点2和点3的血流速度最大值由于惯性作用稍微延后，其余参考点处的血液流速均达到最大值。其中，入口侧血管的观测点1处的血流速度最大，

(a) (b)

(c) (d)

Figure 5. Contour plots of velocity distribution and vector diagrams of velocity distribution on two planes in the aneurysm lumen

图5. 瘤腔中两个平面的速度分布等值线及速度分布矢量图

Figure 6. Characteristic curves of velocity variation with time at each observation point

图6. 各观测点位置的速度随时间变化特征曲线

其数值为0.032 m/s。瘤内点2、点3和点6处由于血流冲击，血液流速分别为0.013 m/s、0.012 m/s和0.016 m/s，其余各点的血液流速均接近零，这表明瘤顶区在整个心动周期内血流处于缓慢流动状态；这一特征表明动脉瘤腔内的血液流动整体处于低速状态。

3.2. 壁面力学特性分析

图7为一个心动周期内载瘤血管的壁面压力分布图。在心缩期上升阶段，瘤腔表面压力分布相对均匀，未出现明显波动；随着血流速度的提升，血流冲击侧逐渐形成高压区域，此时高速血流持续冲击薄弱的动脉瘤壁面，促使壁面扩张，该区域的高流速特性也会增加动脉瘤的破裂风险。在心缩期峰值点时，血管壁面整体压力升至周期内最高水平(压力峰值达14,424.6 Pa)，瘤体被血流完全充盈，血管近心端与远端的压力差进一步扩大。进入心缩期下降阶段，压力呈缓慢回落趋势，瘤体压力虽有所降低但仍维持中高水平，局部已出现压力分布不均的现象，可能伴随血流再循环的初始变化。到舒张期后期，整体压力进一步下降，瘤体压力同步回落且分布更趋平缓，与舒张期血流动力减弱的生理状态相符。在整个心动周期内，血管远端压力始终相对偏低，呈现出血液推进过程中的压力梯度；同时，载流动脉弯曲处的外侧区域压力偏高，这一分布特征与血液的流动规律相一致。

牛顿内摩擦定律指出壁面切应力(WSS)与速度梯度呈正相关特性，结合3.1小节对血流状态的分析

(a) 心缩期上升阶段 (b) 心脏收缩期峰值点

(c) 心缩期下降阶段 (d) 舒张期后期

Figure 7. The wall pressure distribution cloud diagram of tumor-bearing vessels over one cardiac cycle

图7. 载瘤血管在一个心动周期内的壁面压力分布云图

可以发现，动脉瘤腔内血流速度在单个心动周期内未产生显著变化，速度梯度处于较低水平，因此瘤腔壁面的WSS分布整体呈低值特征。已有多项研究证实，壁面切应力(WSS)的数值大小与分布特征，是影响动脉瘤破裂风险的关键因素之一，因此壁面切应力分布规律的研究对临床风险评估具有重要意义。如图8所示，壁面切应力(WSS)在动脉瘤瘤顶区域的分布较为均匀，且整体WSS值偏低，在整个心动周期内未出现明显变化。

(a) 心缩期上升阶段 (b) 心脏收缩期峰值点 (c) 心缩期下降阶段 (d) 舒张期后期

Figure 8. WSS distribution cloud diagram of tumor-bearing vessels over one cardiac cycle

图8. 载瘤血管在一个心动周期内的WSS分布云图

如图9(a)和图9(b)所示，动脉瘤膨胀区域入口附近的壁面切应力(WSS)显著高于出口区域。已有研究结果显示，WSS的分布特征与变化规律对动脉瘤破裂位点的预测具有重要参考意义，且多项研究结果表明，低壁面切应力区域为动脉瘤破裂的高发区域[19]。此外，血管横截面的变化会对WSS产生显著影响，血管狭窄的存在同样会对壁面切应力(WSS)产生影响，且狭窄区域的WSS值会随管腔狭窄程度的增加而升高；如图9(c)所示，在血管狭窄区域，随着血管横截面积的减小，壁面切应力(WSS)值呈上升趋势。

Figure 9. WSS distribution cloud of tumor-bearing vessels

图9. 载瘤血管壁面切应力分布云图

为进一步探究载瘤血管壁不同部位的WSS大小及动态变化规律，基于六个观测点绘制的单个心动周期内各点位的WSS变化曲线图(见图10)。由图可知，瘤顶区及瘤体侧面的观测点(点2、点3、点4) WSS值处于0~3.5 Pa区间；瘤颈区观测点(点5、点6)的WSS值分别为5.3~16.2 Pa、4.2~15.0 Pa；而正常血管段的观测点1 (入口侧血管)在全周期内WSS值最大，范围为7.4~20.2 Pa。对比可见，正常血管壁的WSS值整体高于动脉瘤腔壁，且动脉瘤壁各区域的WSS分布存在明显差异。在心缩期峰值时刻(0.93 s)，瘤颈区域(点5、点6)的WSS值约为瘤顶区域(点2、点3)及瘤体侧面(点4)的3~5倍，即瘤颈处WSS显著高于后两者。

Figure 10. Characteristic curves of dynamic variation of WSS values at each observation point over time

图10. 各观测点WSS值随时间的动态变化特征曲线

瘤顶区长期处于低壁面切应力(WSS)微环境中，这一力学刺激会诱发血管内皮功能紊乱[20]；同时结合前文流场分析结果，瘤顶处血流在整个心动周期内均维持低速状态(流速趋近于0)，导致血小板等血液有形成分易黏附于内皮功能受损的瘤顶内壁，进而形成附壁斑块。该病理变化会进一步导致瘤顶区域血管壁的弹性模量下降，使其力学承载能力减弱。

4. 结论

本研究采用计算流体动力学(CFD)方法，针对载瘤血管数值模型，开展脉动血流工况下的流场特性数值仿真与分析。结果显示，在整个心动周期内，高速血流进入动脉瘤腔后即刻发生流动分离现象，形成显著的回流区域与涡旋结构；此类复杂流场结构不仅加剧了局部血流的动力学不稳定性，还可能通过持续的动力学作用影响血管壁的生物学响应。低壁面切应力(WSS)区域通常与血管内皮细胞功能异常相关，可能促进血管壁的进一步重塑与薄弱化；而高压区的周期性力学作用可能加剧瘤壁的结构疲劳与损伤累积。上述血流动力学参数的异常分布特征，共同揭示了动脉瘤在血流作用下潜在的力学环境，为理解其进展机制提供了重要的流体力学依据。

NOTES

^*通讯作者。