末梢血管血流FSI

分类: 分析 | 综合版 2026-04-06

理论与物理的世界

末梢血管血流FSI的理论基础

概述

🧑🎓

老师，为什么在末梢血管的血流模拟中需要耦合流体和结构呢？

🎓

末梢动脉直径在数毫米以下，血管壁薄，因此血压引起的变形相对较大。壁面膨胀会增加横截面积从而改变流速，流速改变又会改变壁面剪切应力和压力。如果忽略这种双向作用，脉波传播速度和壁面应力的预测精度会大幅下降。

🧑🎓

具体有哪些医学用途呢？

🎓

动脉硬化的进展预测、支架置入后再狭窄风险评估、血管搭桥手术的设计支持等。壁面剪切应力（WSS）低的区域已知是斑块堆积风险高的区域，因此通过FSI获取准确的WSS分布在临床上很重要。

控制方程

🧑🎓

流体侧和结构侧分别使用什么方程？

🎓

流体侧是不可压缩Navier-Stokes方程。血液常被当作非牛顿流体处理。

$$ \rho_f \left(\frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} - \mathbf{v}_m) \cdot \nabla \mathbf{v}\right) = -\nabla p + \nabla \cdot \boldsymbol{\tau} $$

$$ \nabla \cdot \mathbf{v} = 0 $$

这里 $\mathbf{v}_m$ 是ALE（任意拉格朗日-欧拉）网格速度。血液粘度常用Carreau-Yasuo模型描述。

$$ \mu(\dot{\gamma}) = \mu_\infty + (\mu_0 - \mu_\infty)[1 + (\lambda\dot{\gamma})^2]^{(n-1)/2} $$

🧑🎓

结构侧如何描述？

🎓

血管壁被建模为超弹性体。代表性的有Mooney-Rivlin模型和各向异性的Holzapfel-Gasser-Ogden模型。

$$ \Psi = C_{10}(\bar{I}_1 - 3) + \frac{k_1}{2k_2}\sum_{i=4,6}\left\{\exp[k_2(\bar{I}_i - 1)^2] - 1\right\} $$

用两族纤维表示胶原纤维的取向，再现动脉壁的各向异性。

耦合界面条件

🧑🎓

流体和结构在哪里连接？

🎓

在血管内腔面满足三个条件。

1. 运动学条件: $\mathbf{v}_f = \dot{\mathbf{d}}_s$（界面处流体速度 = 结构位移速度）

2. 力学条件: $\boldsymbol{\sigma}_f \cdot \mathbf{n} = \boldsymbol{\sigma}_s \cdot \mathbf{n}$（牵引力的连续性）

3. 几何条件: 流体域形状跟随结构变形

严格满足这三个条件的是强耦合（strong coupling），血管FSI中强耦合是必须的。弱耦合会因附加质量效应而不稳定。

🧑🎓

附加质量效应是什么？

🎓

当流体与结构的密度比 $\rho_f/\rho_s$ 接近1时（血液和血管壁符合此情况），分离型的弱耦合会注入人为能量导致发散。这就是附加质量不稳定性。通过Robin-Neumann分割法或Quasi-Newton法等方法来避免此问题的研究正在进行。

Coffee Break 闲谈

血液是“非牛顿流体”——关于Carreau流体的讨论

在末梢血管血流分析中，初学者容易陷入“能否将血液当作牛顿流体（粘度恒定）处理”的问题。在大血管中大致可以，但在接近细小毛细血管的末梢血管中，当剪切速率降低时，血液粘度急剧上升的Carreau流体特性不可忽视。具体来说，当剪切速率低于1s⁻¹时，粘度会从约0.004Pa·s增加到0.04Pa·s，增长近10倍。忽略这种非线性效应会导致壁面剪切应力的分布出现很大偏差。

数值解法与实现

ALE公式化

🧑🎓

血管变形的话网格会动吧？具体怎么处理呢？

🎓

ALE法中引入网格速度 $\mathbf{v}_m$，修正流体的对流项。界面处网格跟随结构，内部通过拉普拉斯平滑或弹性体类比来维持网格质量。

$$ \nabla \cdot (k \nabla \mathbf{d}_m) = 0 $$

$k$ 是扩散网格位移的系数，诀窍是在界面附近设小（高刚度），远处设大。

🧑🎓

大变形的话网格不会溃缩吗？

🎓

说得对。狭窄率高的病例（70%以上闭塞）可能出现ALE网格失效的情况。这时再网格化（remesh）或不使用网格的浸入边界法（Immersed Boundary法）成为备选方案。

浸入边界法

🧑🎓

浸入边界法是什么原理？

🎓

这是Peskin于1972年为心脏瓣膜模拟开发的方法。将结构嵌入固定的流体网格上，通过δ函数交换力和速度。

$$ \mathbf{f}(\mathbf{x}, t) = \int \mathbf{F}(s,t) \delta(\mathbf{x} - \mathbf{X}(s,t)) ds $$

$$ \frac{\partial \mathbf{X}}{\partial t}(s,t) = \int \mathbf{v}(\mathbf{x},t) \delta(\mathbf{x} - \mathbf{X}(s,t)) d\mathbf{x} $$

无需重新生成网格，因此对大变形有优势，但缺点是界面附近的精度依赖于δ函数的宽度。

1D-3D耦合法

🧑🎓

把末梢血管全部用3D计算太重了吧？

🎓

说得对。目标区域（例如颈动脉分叉部）用3D FSI详细求解，上游和下游用1D模型模拟全身循环。1D模型的控制方程是横截面积平均的守恒律。

$$ \frac{\partial A}{\partial t} + \frac{\partial (AU)}{\partial x} = 0 $$

$$ \frac{\partial U}{\partial t} + U\frac{\partial U}{\partial x} + \frac{1}{\rho}\frac{\partial p}{\partial x} = -\frac{8\pi\mu}{\rho A}U $$

压力-面积关系通过管定律闭合。

$$ p - p_{ext} = \beta(\sqrt{A} - \sqrt{A_0}), \quad \beta = \frac{\sqrt{\pi}Eh}{(1-\nu^2)A_0} $$

🧑🎓

3D和1D的接口怎么处理？

🎓

将3D侧的断面平均流量和断面平均压力作为1D的边界条件传递。在出口边界应用基于特征线法的RCR Windkessel模型是标准做法。Ansys Fluent和SimVascular支持这种耦合。

时间积分与收敛

🧑🎓

强耦合的迭代如何收敛？

🎓

在每个时间步内进行Dirichlet-Neumann迭代。基本做法是交替给结构施加Dirichlet条件（位移指定）、给流体施加Neumann条件（牵引力指定），但会因附加质量效应而发散。

IQN-ILS（基于逆最小二乘的界面拟牛顿）法对加速FSI收敛非常有效。它基于过去的界面残差和更新量构建雅可比矩阵的近似逆矩阵。preCICE库实现了这个方法。

方法	特点	收敛性
定点迭代	简单但收敛慢	在血管FSI中易发散
Aitken松弛	动态松弛加速	中等
IQN-ILS	拟牛顿法	5~10次迭代收敛
Anderson加速	混合法	与IQN-ILS相当

Coffee Break 闲谈

脉动流的边界条件——如何给定心拍波形

在末梢血管FSI分析中，意外令人困扰的是“入口边界条件的给定方法”。要再现与心拍同步的脉动流，使用多普勒超声实测的流速波形是理想选择，但每个患者的数据不同。实际工作中广泛使用“利用Womersley解（分析解）近似周期性速度剖面”的方法。Womersley数（α）超过10的大血管中剖面趋于平坦，3以下的细血管中接近抛物线——理解这个差异后再进行实现，可以减少边界条件的错误。

末梢血管血流FSI末梢血管血流FSI实践指南

🎓

典型的流程是这样的。

1. 获取医学影像: 从CT血管造影（CTA）或MRA获取血管形状

2. 分割: 使用ITK-SNAP、3D Slicer、Mimics等提取血管腔和壁

3. 形状修复与CAD化: 使用Geomagic Wrap等对表面网格进行平滑、缺损修复

4. 生成体网格: 血管壁用六面体单元（3层以上），管腔用棱柱边界层+四面体

5. 设置材料与边界条件: 超弹性材料、脉动入口速度、Windkessel出口条件

6. 执行FSI计算: 强耦合下计算心拍2~3个周期（为去除初始瞬态）

7. 后处理: 评估WSS、OSI（振荡剪切指数）、壁位移