老师，为什么在末梢血管血流模拟中需要将流体与结构进行耦合？

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

具体有哪些医学应用？

主要用于动脉硬化进展预测、支架植入后再狭窄风险评估、血管搭桥手术设计支持等。已知壁面剪切应力（WSS）较低的区域斑块沉积风险较高，因此通过流固耦合（FSI）精确计算WSS分布具有重要的临床意义。

末梢血管血流FSI

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

流体侧采用不可压缩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} $$

Q: 结构侧如何描述？

血管壁建模为超弹性体，常用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\} $$ 该模型通过两组纤维方向表征胶原纤维排列，再现动脉壁的各向异性特性。

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

理论与物理的世界

理论与物理

概述

🧑‍🎓

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

🎓

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

🧑‍🎓

具体有哪些医学用途呢？

🎓

动脉硬化的进展预测、支架置入后再狭窄风险评估、血管搭桥手术的设计支持等。壁面剪切应力（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倍。忽略这种非线性效应会导致壁面剪切应力的分布出现很大偏差。

各项的物理意义

结构-热耦合项：温度变化引起的热膨胀诱发结构变形，变形又影响温度场。$\sigma = D(\varepsilon - \alpha \Delta T)$。【日常例子】夏天铁轨伸长导致间隙变窄——温度上升→热膨胀→产生应力的典型例子。电子基板在焊接后翘曲也是不同材料热膨胀率差异导致的。发动机的气缸体因高温部和低温部的温差产生热应力，最坏情况下会导致裂纹。
流体-结构耦合（FSI）项：流体压力·剪切力使结构变形，结构变形改变流体域的双向相互作用。【日常例子】强风下吊桥缆索振动（涡激振动）——风力使结构摇晃，摇晃的结构改变风流，进而放大振动。心脏血流与血管壁的弹性变形、飞机机翼的颤振（气动弹性不稳定性）也是典型的FSI问题。有时单向耦合即可，但变形大时双向耦合是必须的。
电磁-热耦合项：焦耳发热 $Q = J^2/\sigma$ 引起温度上升，温度变化改变电阻的反馈循环。【日常例子】电暖炉的镍铬丝通电后发热（焦耳热）变红——温度上升电阻改变，电流分布也变化。IH电磁炉的涡流发热、输电线路温度上升导致的垂度增加也是这种耦合的例子。
数据传递项：通过插值解决不同物理场间网格不一致的问题。【日常例子】天气预报中结合“气温数据”和“风数据”计算体感温度时，如果各自的观测地点不同就需要插值——CAE的耦合分析中，结构网格和CFD网格通常也不一致，因此界面处的数据传递（插值）精度直接关系到结果的可信度。

假设条件与适用范围

弱耦合假设（单向耦合）：一方物理场影响另一方但反向可忽略时有效
需要强耦合的情况：FSI中的大变形、电磁-热耦合中温度依赖性强的场合
时间尺度分离：各物理场特征时间差异大时，可通过子循环提高效率
界面条件一致性：需确认耦合界面处能量·动量守恒在数值上得到满足
不适用的场合：三个以上物理场同时强耦合时，有时需要整体式方法

量纲分析与单位制

变量	SI单位	注意事项·换算备忘
热膨胀系数 $\alpha$	1/K	钢：约12×10⁻⁶，铝：约23×10⁻⁶
耦合界面力	N/m²（压力）或N（集中力）	确认流体侧与结构侧的力平衡
数据传递误差	无量纲（%）	插值精度依赖于网格密度比。5%以下是目标

数值解法与实现

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以下的细血管中接近抛物线——理解这个差异后再进行实现，可以减少边界条件的错误。

整体式方法

将所有物理场作为一个联立方程组同时求解。对强耦合稳定，但实现复杂，内存消耗大。

分区法（分离迭代法）

各物理场独立求解，在界面交换数据。易于实现，可利用现有求解器。适用于弱耦合。

界面数据传递

最近邻法（最简单但精度低）、投影法（守恒性好）、RBF插值（对网格不一致鲁棒性强）。需要在守恒性和精度间取得平衡。

子迭代

在每个耦合步内进行充分迭代，确保界面条件的一致性。残差基准基于各物理场的典型值进行缩放。

Aitken松弛

自动调整耦合迭代的松弛系数。防止过度松弛导致发散，是加速收敛的自适应方法。

稳定性条件

注意附加质量效应（流体-结构耦合中结构密度≈流体密度时）。不稳定时可应用Robin型界面条件或IQN-ILS法。

Aitken松弛的比喻

Aitken松弛类似于“平衡跷跷板”。一方推得太用力，另一方就会弹起，反作用力又导致推得更用力——为了抑制这种振荡，自动调整推力大小的就是Aitken松弛。当耦合迭代振荡不收敛时，根据上次的修正量自动调整下次修正量的自适应方法。

实践指南

🎓

典型的流程是这样的。

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（振荡剪切指数）、壁位移