随着我国进入中度老龄化社会^［1］，动脉粥样硬化等中老年疾病占比也逐渐增加。动脉粥样硬化是一种由于胆固醇、脂肪和其他物质堆积导致动脉内壁狭窄的慢性疾病。堆积而成的斑块随时间不断增大，最终导致血管狭窄、闭塞甚至血栓^［2］。腘动脉作为人体下肢动脉系统的一条重要动脉，连接着大腿动脉和小腿动脉。腘动脉狭窄会限制下肢的血液供应，引起间歇性跛行、坏疽、溃疡等并发症^［3］，严重者面临截肢风险。研究表明^［4–5］，血管血流速度分布、壁面剪切力等血流动力学参数与动脉粥样硬化的形成及发展密切相关。

临床常用的计算机断层动脉造影（CTA）、核磁共振动脉造影（MRA）等虽然能显示病变部位和狭窄程度，但并不能提供较为全面的腘动脉血流动力学信息。计算流体力学的发展为血流流动力学带来了极大的便利，因其能提供精确的数值解，揭示血流的详细特征，并且有着无创的特点，成为研究和分析血流的重要工具之一^［6–7］。Gökgöl等^［8］研究了不同年龄段人群腘动脉血流动力学参数的变化。Desyatova等^［9］研究了腿部弯曲对腘动脉再狭窄的影响。以上研究把血液假设为牛顿流体，并将血管壁面设置为刚性壁面进行血流动力学研究。但在真实人体中，血管作为有弹性的壁面对血液也有着相互作用^{［10–11］}，血液由于其特殊的组成实际是一种非牛顿流体^{［12–13］}。本文对一例腘动脉狭窄病人进行三维血管重建，分别使用牛顿流体模型、Carreau-Yasuds非牛顿流体模型进行流固耦合(FSI)数值模拟，通过对比牛顿流体特性血液来探讨血液非牛顿流体特性对腘动脉狭窄血流动力学的影响。

1 材料和方法 1.1 几何模型

选用一例50岁腘动脉轻度狭窄病人的CT图像，采用MIMICS医学处理软件将二维切片图像重建为STL三维模型。将模型导入ANSYS-ICEM软件进行固体和流体网格划分，固体域采用四面体网格进行划分，流体域采用四面体^［14］与六面体混合网格划分，近壁处为7层六面体贴壁网格，网格模型如图1所示。最终通过ANSYS workbench软件中Fluent与瞬态结构模块(Transient structural)耦合进行计算。其中牛顿流体模型与非牛顿流体模型中均采用相同的计算网格、边界条件进行设置。

1.2 控制方程

假设血液为不可压缩流体，流动状态为层流，且忽略重力影响。不可压缩流体的Navier-Stokes方程为

$ \qquad \nabla {\text{·}} {{\boldsymbol{u}}} = 0 $

(1)

$ \qquad \rho \left( {\dfrac{{\partial {{\boldsymbol{u}}}}}{{\partial t}} + \left( {{{\boldsymbol{u}}} {\text{·}} \nabla } \right){{\boldsymbol{u}}}} \right) = - \nabla p - \nabla {\text{·}} {{\boldsymbol{\tau}}} $

(2)

式中：${\boldsymbol{u}}$为速度矢量；$\rho $为血液密度；$p$为压力；${\boldsymbol{\tau}}$为应力张量；t为时间。

固体(血管壁)控制方程为

$ \qquad \nabla {\text{·}}{\boldsymbol{\sigma}}_{\text{s}}={\rho }_{\text{s}}{\text{·}}{\boldsymbol{\alpha}}_{\text{s}} $

(3)

式中：${{\boldsymbol{\sigma}}_{{\mathrm{s}}}}$为血管壁应力张量；${\rho _{{\mathrm{s}}}}$为血管壁密度；${{\boldsymbol{\alpha}}_{{\mathrm{s}}}}$为血管壁加速度。

在流固耦合计算中，流体域和固体域通过交界面来传递速度和位移，因此交界面的控制方程为

$ \qquad {{{{\boldsymbol{d}}}}_{{\mathrm{s}}}} = {{{{\boldsymbol{d}}}}_{{\mathrm{f}}}} $

(4)

$ \qquad \boldsymbol{\sigma}_{\mathrm{s}} {\text{·}} \boldsymbol{n}_{\mathrm{s}}=\boldsymbol{\sigma}_{\mathrm{f}} {\text{·}} \boldsymbol{n}_{\mathrm{f}} $

(5)

$ \qquad {{\boldsymbol{u}}_{{\mathrm{s}}}} = {{\boldsymbol{u}}_{{\mathrm{f}}}} $

(6)

式中：d_s、d_f分别为固体域和流体域血管壁面位移；${\boldsymbol{n}}$_s、${\boldsymbol{n}}$_f分别为固体域和流体域血管壁边界法向；u_s、u_f分别固体域和流体域速度；$ {{\boldsymbol{\sigma}}}_{{\mathrm{s}}} $、$ {{\boldsymbol{\sigma}}}_{{\mathrm{f}}} $分别为固体域和流体域应力张量。

1.3 流体模型 1.3.1 牛顿流体模型

牛顿流体是黏度恒定，剪切应力与速度梯度成正比的一种流体。计算时设置血液为牛顿流体，其密度为1 060 kg·m⁻³，黏度为3.5 mPa·s。

1.3.2 非牛顿流体模型

非牛顿流体的黏度不符合牛顿流体定律，其黏度随流动条件的变化而变化。血液主要由红细胞、白细胞、血小板和血浆组成，在流动过程中这些成分相互作用，使血液表现出非牛顿性质。Carreau-Yasuds模型^［15］是一种常用的非牛顿血液流动模型，该模型中血液黏度服从本构方程^［16］，即

$ \qquad \mu = {\mu _\infty } + \left( {{\mu _0} - {\mu _\infty }} \right){\left[ {{1} + {{\left( {\lambda x} \right)}^a}} \right]^{\tfrac{{n - 1}}{a}}} $

(7)

式中：$\mu $为非牛顿流体黏度；${\mu _0}$为剪切率为0时的流体黏度，取${\mu _0} = {22}$ mPa·s；${\mu _\infty }$为剪切率趋于无穷大时的流体黏度，取${\mu _\infty } = {2}{.2}$ mPa·s；$\lambda $为时间常数，取$\lambda = {11}$；$a$为模型指数，取$a = {0}{.392}$，$\dfrac{{n - {1}}}{a}$为幂律指数，取0.644，n约为1.252；x为剪切率。

1.4 边界条件及材料属性

利用多普勒超声测得的腘动脉血流速度波形如图2所示。在1个标准心动周期0.7 s内，采用傅里叶函数拟合得到的入口速度波形如图3所示，并将其作为入口边界条件。压力出口作为出口边界条件，取出口压力p=0 Pa^［17］。计算时长取3个心动周期共2.1 s，并取最后一个心动周期作为研究对象。

设置腘动脉血管壁为各向同性无滑移的线弹性材料，密度为1 150 kg·m⁻³，杨氏模量为1.32 MPa，泊松比为0.45^［18］，出入口管壁设置为固定支撑，血管壁内壁面为流固耦合交界面。

1.5 计算设置及网格验证

为保证计算结果的准确性，进行网格尺度及时间步长无关性验证。分别采用0.6、0.5、0.4、0.3 mm四套不同尺度的网格进行计算，其他条件不变。结果发现：当最小网格为0.4 mm时，即流体域和固体域网格数量分别为335 605、333 484时，其速度和壁面切应力( WSS)不再随网格数量增加而发生变化；采用0.01、0.008、0.005 s三种时间步长进行时间步长无关性验证时发现，当时间步长取0.01 s时，速度和WSS不再随时间步长的减小而发生变化。

2 数值模拟结果 2.1 速度分布

分别选择收缩期加速点（0.04 s）、收缩期峰值处（0.08 s）、逆流最低值处（0.28 s）以及心动周期终点（0.70 s）四个代表性时间点进行分析。图4为各时刻不同模型速度分布。

由图4可以看出，无论是牛顿流体模型还是非牛顿流体模型，在收缩期加速点(0.04 s)时，血流速度不断加快，主腘动脉和左侧血管狭窄处均出现速度较大的流动区域，牛顿流体模型与非牛顿流体模型的速度场整体相似，但非牛顿流体模型的高速流动区域略小于牛顿流体区域，并存在较大面积的低速区域；收缩期峰值处(0.08 s)时，动脉血流速度达到峰值，受左侧腘动脉狭窄影响而产生的高速血流冲击下游血管壁面，容易造成动脉壁损伤；逆流最低值处(0.28 s)流速显著降低，壁面速度梯度更明显；舒张末期(0.70 s)时，由于较低的流速，非牛顿流体模型的结果表现出更高的黏度，牛顿流体模型与非牛顿流体模型的速度场差异更加明显，非牛顿流体模型的速度分布没有牛顿流体模型的均匀，并出现部分低流速区域。

2.2 WSS分布

WSS是衡量血流与动脉内皮细胞之间摩擦力的血流动力学参数，过高和过低的WSS都会对血管产生不利影响^［19］，例如高WSS会损害内皮细胞、低WSS更容易聚集脂质物质产生动脉粥样硬化等。图5为各时刻不同模型壁面剪切力云图。

如图5所示，在收缩期加速点(0.04 s)时，除小部分高WSS区域外，大部分表面均处于低WSS；随着血流速度不断增加达到峰值(0.08 s)，狭窄处的高WSS会促进巨噬细胞产生铁蛋白酶^［20］，使斑块纤维帽减弱，进而演变为血栓；在逆流最低值处(0.28 s)，血流相对较平缓，狭窄和弯曲处的WSS有所降低；舒张末期(0.70 s)时，牛顿流体模型与非牛顿流体模型的WSS分布差异最大，非牛顿流体模型中由于黏性力作用，对壁面产生较大剪切力，WSS高于牛顿流体模型的。

高入口速度时牛顿流体模型与非牛顿流体模型的WSS分布相似，低入口速度时，牛顿流体模型的WSS低于非牛顿流体模型的。以上发现也验证了Dubey等^［21］的研究结果，进一步确认了本文结果的可靠性。相比于牛顿流体模型，非牛顿流体模型结果在血管内显示更大的低WSS区域，并且其WSS梯度分布更加平滑，也更贴近实际血液特性。

2.3 流线分布

血管的弯曲狭窄常伴随着流动分离，使得层流的流动状态被破坏，因此选择在流速最高，即收缩期峰值处(0.08 s)分析血管流线分布。图6为收缩期峰值处(0.08 s)不同模型的流线分布。由图可知，牛顿流体模型的流线分布在狭窄及壁面处出现紊流，非牛顿流体模型的流线分布在狭窄处速度增大，但整体流动较平滑，也更均匀。

为了进一步探讨截面上牛顿流体与非牛顿流体模型的流动特点，在狭窄的上游和下游选取三个截面(A−A、B−B、C−C)，收缩期峰值处(0.08 s)各剖面速度分布如图7所示。相比于牛顿流体模型，非牛顿流体模型的流速较低，这与Frolov等^［22］对于牛顿流体模型比非牛顿流体模型的截面速度高出22.6% ～ 75.3%的研究结果一致。非牛顿流体模型中由于黏度变化，在一定程度上减少了血液的流动分离和扰动，因此非牛顿流体模型的流线分布和速度分布均更加稳定，更符合正常人体的血液流动。

3 结论

针对腘动脉狭窄，分别使用牛顿流体模型和Casson-Yasuds非牛顿流体模型进行流固耦合数值模拟分析。结果表明：①由于非牛顿流体模型中的黏性力作用，其速度分布出现低流速区域。②非牛顿流体模型的血管壁面产生较大剪切力，壁面剪切力数值高于牛顿流体模型。③牛顿流体模型的流线分布在狭窄处出现紊流，而非牛顿流体模型的流线分布较牛顿流体模型的更平滑、均匀。通过以上研究发现，非牛顿流体模型比牛顿流体模型更符合真实人体血流情况，可为临床诊断和治疗提供分析依据。因此，在未来下肢动脉血流研究中，血液的非牛顿特性应被考虑在内。