血栓栓塞的综合性数值模型：通过耦合计算流体动力学（CFD）与周围动力学（Peridynamics）框架研究流体与血栓之间的相互作用

《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》：A Comprehensive Numerical Model of Thrombus Embolization: Fluid–Thrombus Interactions Through a Coupled Computational Fluid Dynamics–Peridynamics Framework

　　阿比谢克·卡尔马克卡尔（Abhishek Karmakar）、格雷格·W·伯格林（Greg W. Burgreen）、奥利维尔·德贾尔丹（Olivier Desjardins）、詹姆斯·F·安塔基（James F. Antaki）
康奈尔大学梅宁生物医学工程学院，伊萨卡，纽约州14850，美国

**摘要**
血栓栓塞性疾病约占全球死亡人数的三分之一，其特点是血栓的栓塞——即血栓在流体驱动力作用下发生凝聚/粘附性断裂的过程。对这一多物理现象进行全面的数值建模一直具有挑战性，且可靠的计算框架也非常稀缺。本研究通过提出一种将拉格朗日非平凡状态量子动力学（Lagrangian non-ordinary state-based peridynamics）与欧拉计算流体动力学（Eulerian computational fluid dynamics）相结合的框架，显著推动了血栓栓塞数值建模的进步。该框架能够：（i）纳入复杂的宏观尺度本构定律来描述血栓；（ii）捕捉血栓与接触表面之间的剥离现象；（iii）将血栓的结构动力学与在任意非结构化网格上离散化的流体流动结合。后一项能力是通过开发一种用于多边形网格与量子动力学粒子云之间快速双向数据传输的稳健插值方案实现的。这也有助于扩展量子动力学与流体耦合数值模型的研究范围，因为现有的模型大多局限于基于键合的量子动力学和笛卡尔网格。该框架通过五个具有计算独特性的基准测试进行了严格验证。验证后的框架被用于再现托宾（Tobin）等人关于在不同流速下聚碳酸酯圆柱管内预形成血栓栓塞的实验结果。为了证明该数值框架能够支持空间异质材料属性（这是其他方法所不具备的能力），我们还报告了一个合成测试案例。

**引言**
血栓栓塞——即血栓在血流作用下的碎裂——是多种危及生命的心血管疾病（如肺栓塞[1]、心肌梗死[2]和缺血性中风[3]）的先兆。这些疾病占全球死亡人数的约三分之一[4]。尽管统计数据如此严峻，但对于理解和建模这一关键现象的关注相对较少。这一研究空白源于将流体的非线性力学、软组织断裂和动态血栓生长耦合起来的复杂性。文献中少数研究血栓变形的相关内容如下：
- 福格尔森（Fogelson）和盖伊（Guy）[5]提出了一种基于笛卡尔网格的二维浸没边界方法，该方法能够捕捉血小板与血管壁之间的微观尺度粘附和凝聚相互作用。他们通过弹簧来模拟血小板间的相互作用，但这不足以反映血栓的复杂本构响应[6]。
- 徐（Xu）等人[7]提出了一个用于血栓形成和栓塞的二维相场模型，其中血栓被建模为Kelvin–Voigt材料。
- 托宾（Tobin）等人基于多相框架提出了一个三维血栓栓塞数值模型，将血栓建模为Phien–Thanner粘弹性材料[8]。这两种基于连续介质的方法都受到将血栓视为欧拉场这一假设的限制。
- 使用欧拉框架求解血栓结构动力学方程存在问题：（i）质量守恒受到破坏；（ii）无法结合血栓固有的（通常是点状的）异质材料属性；（iii）张量量（如变形梯度和应力）的传输误差分析十分复杂。
- 为克服这些缺点，需要高度精细的网格以最小化数值扩散，但这严重限制了基于欧拉的数值框架在求解大变形结构动力学问题时的适用性。因此，需要采用拉格朗日公式来模拟血栓的结构响应，这样就可以直接赋予空间异质材料属性，自然实现质量守恒，并完全消除与结构方程求解相关的数值扩散。

本文提出的数值框架是一个全面的血栓栓塞模型，具备四个关键能力：（i）能够容纳任意复杂的宏观尺度超弹性本构模型以捕捉血栓的非线性力学响应；（ii）表示血栓的任意空间材料异质性；（iii）显式建模血栓与壁面的剥离现象；（iv）对流体网格离散化方式具有通用性。我们通过使用量子动力学（peridynamics, PD）求解固体力学控制方程，并将其与在非结构化多边形网格上使用有限体积方法求解的流体流动方程相结合来实现这些能力。接下来我们将论证选择PD作为捕捉血栓结构响应的建模框架的合理性。

**经典有限元方法（FEM）**的控制方程以偏微分方程的形式表述，在位移不连续处需要特殊的数值处理。节点分割（nodal splitting [9]）和粘聚区方法（Cohesive Zone Method, CZM [10]）扩展了FEM框架以处理此类不连续性。但在血栓栓塞中，裂纹的起始位置和传播路径完全受流体动力载荷的影响，在模拟开始前无法预先知晓——这使得基于节点分割的方法无法实现裂纹的网格独立表示[11]，且CZM的计算开销非常高，尤其是在血栓栓塞导致的广泛碎裂情况下更为显著。扩展的有限元方法（XFEM [12]）消除了对裂纹符合网格的要求，但需要显式跟踪裂纹前沿，通常通过水平集方法（level-set methods [13], [14]）来实现。水平集的数量与需要解决的裂纹数量线性相关，这使得其在涉及血栓断裂的复杂三维碎裂场景中的应用变得困难。
- 相场方法（Phase Field Methods, PFMs [15], [16]）通过将断裂建模为损伤积累过程来规避裂纹跟踪问题，因此裂纹自然表现为扩散的损伤场。然而，PFMs需要在潜在断裂区域使用高度精细的网格来解决正则化长度尺度问题，这需要预先知道断裂可能发生的位置[17]。更广泛地说，所有基于网格的断裂方法都需要在裂纹尖端周围进行网格细化，这一限制在考虑血栓生长时更加明显，因为不断变化的血栓几何形状要求在每个时间步长都重新生成网格，从而使这类方法在计算上变得不可行。这些考虑促使我们采用无网格框架来模拟血栓栓塞的结构物理。

在无网格方法中，物质点方法（Material Point Method, MPM [18]）和平滑粒子流体动力学（Smoothed Particle Hydrodynamics, SPH [19]）是模拟大变形和断裂问题的突出方法。在MPM中，断裂要么通过允许裂纹表面切割的节点处存在多个速度场来模拟（类似于FEM中的节点分割方法），要么通过不显式表示裂纹表面的连续方法来模拟[20]。前者需要跟踪不断变换的裂纹表面，对于复杂的裂纹模式来说实现起来特别复杂。后者在裂纹方向与背景网格单元方向一致时才能获得准确的解[21]。SPH存在拉伸不稳定性问题[22]以及在施加关键边界条件时遇到困难[23]，从而限制了其在固体断裂力学中的准确性。因此，我们选择使用PD。

PD是对经典连续介质力学（Continuum Mechanics, CCM）的非局部重构形式，它以积分方程而非偏微分方程的形式表述[24]，这使得控制方程在位移不连续处能够得到良好的定义。断裂和损伤通过粒子间的键合失效自然地得到模拟，更重要的是，它不需要裂纹前沿跟踪、显式的裂纹分支和合并准则或网格拓扑更新。此外，PD通过参数（δ，也称为视界）的影响半径具有内在的非局部性，从而隐式地规范了断裂过程区域，避免了经典损伤模型的网格依赖性损伤定位[25]。因此，它非常适合捕捉在动态流体动力载荷作用下的血栓断裂，因为断裂位置事先是未知的。自2000年引入以来，PD已被用于模拟多种材料的断裂[26], [27], [28]，包括生物材料[29], [30]。PD有多种变体：（i）基于键合的（Bond-based, BB-PD）；（ii）基于平凡状态的（Ordinary State-based, OSB-PD）；（iii）基于非平凡状态的（Non-Ordinary State-based, NOSB-PD）[24], [31], [32]。在这些公式中，BB-PD将三维泊松比限制为ν=0.25，这与血栓的近乎不可压缩性（ν>0.4）不相容。OSB-PD部分放宽了这一限制，但将力状态限制在键合方向上，无法直接纳入任意的CCM本构定律——这对于需要模拟血栓复杂材料响应的数值框架来说是一个不可妥协的要求。NOSB-PD通过对应原理[31], [33]克服了这两个限制，提供了将任意CCM本构定律转换为PD框架的系统方法，无需重新表述PD的控制方程。因此，我们选择NOSB-PD作为求解血栓结构动力学的数值框架。在我们之前的工作中[30]，我们开发了第一个专门针对血栓本构建模设计的NOSB-PD框架，并通过标准基准测试进行了严格验证，随后根据已发表的实验数据进行了校准。本研究直接采用了该经过验证的框架。

本文提出的数值框架将NOSB-PD粒子云与在任意非结构化多边形网格上求解的流体流动相结合。这比现有的PD-流体耦合框架有了显著进步，因为现有的框架中的流体控制方程是在笛卡尔网格上离散化和求解的[34], [35]，从而将其适用性限制在简单几何形状上。这一限制并非PD-流体框架所独有，因为基于直接激励浸没边界方法的一般粒子-流体耦合策略[36], [37], [38], [39], [40]也主要在笛卡尔网格上得到验证。尝试将这些框架推广到非结构化网格（如通过修改标准δ函数扩散器和积累算子[41]或通过切割单元方法[42]）会引入显著的计算开销。前者需要求解大型矩阵系统，而后者涉及复杂的接口重建算法，对于动态演变的固体几何结构（如断裂的血栓）来说计算负担很重。此外，使用笛卡尔网格限制了经典耦合技术在大多数临床相关几何结构（如心室辅助装置（VADs）或患者特定血管中的应用，这些结构需要使用非结构化多边形网格进行准确表示。由于必须在每个时间步长执行流体-固体力交换，因此插值过程的计算效率至关重要。为了解决这些限制，我们提出了一种计算效率高、适用于将粒子云与任意流体网格耦合的无矩阵双向插值技术。虽然这里是为NOSB-PD粒子云展示的，但所提出的插值方案具有通用性，适用于任何需要与任意多边形网格上的有限体积求解器耦合的拉格朗日粒子系统。此外，在笛卡尔网格假设下，该方案可以恢复为经典耦合方法的特殊情况，证明了与其现有方法的一致性。进一步使用所提出的插值算子，明确的推导表明该耦合方案主要是耗散的，这对于涉及严重断裂的流体-结构相互作用（FSI）问题的数值稳定性至关重要。

**剩余文本的结构如下：**
第2.1节“基于非平凡状态的量子动力学（NOSB-PD）”和第2.2节“流体流动方程”介绍了PD的基本概念，以及实现流体-结构方程耦合所需的修改。第2.4节“插值”和第2.5节“数值算法”详细介绍了插值和数值算法过程。该框架通过一系列基准测试进行了严格验证，检验了其在以下方面的能力：（i）在CFD网格和粒子云之间插值时保持空间流动结构；（ii）将流体限制在PD粒子云内；（iii）强制粒子云速度并有效再现周围流体的响应；（iv）预测钝体周围的流体阻力；（v）处理浸没在矩形通道中的大变形物体的流动。经验证的框架随后被用于再现托宾等人的实验观察结果，即在不同流速下聚碳酸酯圆柱管内预形成血栓的栓塞现象。最后，讨论了该框架所依据的假设含义，并通过一个合成测试案例展示了该框架表示和保持空间异质材料属性的能力。我们通过一个三维半环形管道（即弯曲管道）的流动数值测试案例来验证我们的插值程序。该计算流体动力学（CFD）网格由六面体元素组成，包含了边界层和自适应元素。其几何形状和尺寸如图4(a)所示。为了在完整域中获得非平凡的不对称流动结构，在蓝色阴影部分（图4(a)）将流入边界条件设置为（0 0 0）m/s，而在红色阴影部分，则将流速设置为（讨论部分本文首次提出了一个数值框架，该框架将基于任意非结构化多面体网格的CFD分析与粒子动力学（PD）相结合，用于解决和预测在流体动力作用下的血栓破碎和分层现象。在本节中，我们讨论了构成该数值框架的各个组成部分的准确性。结论本文通过将非正常状态基础上的粒子动力学与流体流动方程相结合，提出了一个全面的血栓栓塞数值模型。该框架能够解析血栓栓塞过程中的两种竞争性物理现象：(a) 血栓内部的凝聚性断裂；(b) 血栓从固定表面的粘附性断裂。这种耦合方法不依赖于流体网格的拓扑结构，这得益于一种高效且稳健的插值算法。RediT作者贡献声明Abhishek Karmakar：撰写——初稿、软件开发、方法论设计、概念框架构建；Greg W. Burgreen：撰写——审阅与编辑、软件实现；Olivier Desjardins：撰写——审阅与编辑、方法论设计；James F. Antaki：撰写——审阅与编辑、项目指导及资金申请。利益冲突声明作者声明他们没有已知的财务利益冲突或可能影响本文研究工作的个人关系。致谢本项工作得到了美国国立卫生研究院（资助编号R01 HL089456）的支持。OD衷心感谢ONR通过高超音速空气热力学、高速推进和材料计划（N000142512177项目）提供的资金支持。JFA则感谢美国陆军医学研究采购活动（HT9425-25-1-0142项目）的资助。