拉格朗日方法对血液损伤预测的影响

拉格朗日方法用于心血管装置中的血液损伤评估——数值实现如何影响结果

摘要

本文评估了使用拉格朗日方法结合欧拉计算流体动力学预测心血管装置中血液损伤时，各种数值实现假设对结果的影响。所测试的实现假设包括不同的示踪粒子分布模式、随机游走模型以及基于路径线的简化轨迹计算。所评估的后处理实现选项包括单次通过和重复通过的应力累积以及时间平均。本研究表明，这些实现假设会显著影响最终的应力累积，即血液损伤模型的预测结果。在使用拉格朗日模型时应谨慎考虑。最终，应根据具体案例的物理特性选择适当的假设，并采用类似本文所示的敏感性分析。

关键词

血液损伤模型；计算流体动力学；拉格朗日方法；心血管装置；数值建模

1. 引言

已知心血管装置会引起需要抗凝治疗的血栓栓塞并发症。实验研究将这种血液损伤与升高的流体剪切应力水平以及血液成分对其的暴露时间相关联 [1,2]。大多数提出的血液损伤或溶血模型均遵循布莱克希尔等人 [3] 提出的幂律公式形式：

$$
D = C \cdot \tau^\alpha \cdot t^\beta
$$

其中，$D$ 是损伤指数，$\tau$ 是剪切应力，$t$ 是暴露时间，$C$、$\alpha$ 和 $\beta$ 是模型系数 [3–7]。

我们的研究小组采用了线性溶血模型（$C=\alpha=\beta=1$；[5]）来研究血小板活化，该模型已用于多种心血管装置中血流的数值模拟，例如机械心脏瓣膜(MHV) [8–10]和机械循环支持(MCS)设备[11,12]。为此，我们采用了多相离散相模型(DPM)，并计算了大量血小板的拉格朗日轨迹。尽管这些在装置中的仿真为评估其血栓形成潜力提供了宝贵的量化信息，但此类模型中数值实现方法的影响从未得到系统性的研究。本文旨在评估拉格朗日损伤模型在与传统欧拉CFD模型耦合时，其数值实现方法如何影响结果。这些影响通过使用我们的线性应力累积模型进行敏感性研究来评估。我们希望提醒注意，拉格朗日损伤模型预测准确血液损伤的能力在很大程度上取决于其复杂程度，且准确的预测能力可能会变得难以捉摸，需要通过一系列复杂的体外实验进行验证。然而，拉格朗日损伤模型的一个非常重要的用途在于比较性研究（例如，在相同数值方法下比较某种装置两种相似设计的血栓形成性能）。在这类研究中，比较可以消除先验实现假设的影响，同时仍保持模型的预测能力。

2. 方法

我们对多种数值实现方法进行了比较，并研究了它们对拉格朗日细胞损伤模型结果的影响。这些比较共使用了九种专用的仿真模型。这些模型采用两种流动状态：（i）通过冠状动脉支架的层流；（ii）通过机械心脏瓣膜的湍流。这两种装置代表了可能影响结果的典型实施假设情况。对于这两种情况，均采用离散相模型（DPM）结合粒子注入来计算拉格朗日轨迹。在装置入口处释放代表自然浮力血小板的球形颗粒，其直径为3 微米，密度与周围流体相同。在拉格朗日参考系中，通过计算作用在颗粒上的惯性和阻力之间的力平衡来确定血小板的轨迹。记录这些血小板颗粒在其相应流动轨迹上的应力历史，并据此计算每条血小板轨迹的线性应力累积（SA）。本文所使用的应力累积（SA）定义为应力对时间的积分，可通过瞬时应力幅值与暴露时间乘积的求和进行数值近似，即 [5]：

$$
SA = \int_0^T \tau(t) dt \approx \sum_{i=1}^{n} \tau_i \Delta t
$$

其中 $T$ 为轨迹持续时间，$t$ 为时间，$\Delta t$ 为时间步长。显然，也可采用其他非线性模型以类似方式测试不同实施方案的影响。各应力张量分量的贡献被转化为标量应力，这一方法最初由布卢德祖韦特提出[13] ，随后被阿佩尔等人[14]及我组[8]所应用。该标量应力基于单轴剪切流中单元变形的功与三维变形的功之间的比较：

$$
\tau = \sqrt{\frac{1}{2} \left( (\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i})^2 \right)}
$$

对于湍流情况，应力张量包含粘性应力和雷诺应力，而雷诺应力基于布辛涅斯克涡流粘度假设进行近似：

$$
\bar{\tau} {ij} = \mu \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i} \right) - \rho \overline{u’_i u’_j} = \mu \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i} \right) - \frac{2}{3} k \delta {ij} + 2 \mu_t S_{ij}
$$

其中，$\mu$ 是粘度，$\frac{\partial \bar{u}_i}{\partial x_j}$ 是平均速度梯度，$\rho$ 是密度，$\overline{u’_i u’_j}$ 是雷诺应力，$k$ 是湍流动能，且 $\mu_t$ 是湍流粘度。

数千个血小板在多条轨迹上达到的SA值的统计分布被归约为一个概率密度函数（ PDF）[9]。PDF方法描述了单个血小板轨迹达到特定SA值的可能性的统计分布，代表了装置的“血栓形成足迹”，并便于在不同装置或其设计迭代之间进行比较。尽管这种表示装置“血栓形成足迹”的方法最初是为不同器械设计之间的比较而提出的，

这里以一种 somewhat 非传统的方式使用，用于比较不同实现假设的影响。这种用法使我们能够专注于实现假设如何影响结果。释放的颗粒数量约为 11,500，该数值基于细化研究选定，以确保溶血指数的概率密度函数对此数量不敏感。

尽管这两类模型均通过求解瞬态流动方程进行计算，但为了简化并提高模型之间的可比性，边界条件设定为恒定值。在所有模型中，血液均被视为不可压缩牛顿流体，密度为 1060 m³/kg，粘度为0.035 帕·秒。层流情况模拟的是典型冠状动脉尺寸（直径为3.3毫米），入口速度为0.3 米每秒，对应的雷诺数为150[15,16]。支架长度为8 毫米，厚度为100 μm [17]，计算域长度为100 毫米，采用220万个多面体单元进行网格化——该网格划分基于针对回流区的敏感性研究（图1）。模拟持续时间为1.5 秒，依据是被捕获在支架梁之间的颗粒数量，时间步长采用固定的3 毫秒。所选择的颗粒数量能够以数百个在整个模拟持续时间内被捕获的颗粒来准确捕捉被捕获颗粒的现象。机械瓣膜中的湍流采用剪切应力输运（SST）k−ω模型并结合低雷诺数修正进行模拟。该模型因其能够准确描述各类非定常雷诺平均纳维‐斯托克斯（URANS）中的现象，适用于此类瓣膜内的流动[8]。直管和瓣膜的直径均为19 毫米。网格划分和时间离散化以及边界条件均基于阿莱穆和布鲁斯坦的工作[8] ，同时采用具有峰值收缩期幅度的恒定入口速度。计算域长度为100 毫米，且使用95万个多面体和15万个六面体单元进行网格化，第一层边界层的厚度根据湍流强度 [8]选择。仿真持续时间为3秒，固定时间步长为1毫秒[8]。所有模型均在ANSYS Fluent 16.1（ANSYS公司，宾夕法尼亚州卡农斯堡）中求解。

3. 直管的影响

在某些情况下，不仅需要将装置与其他设计进行比较，还需要与对照情况进行比较。一个典型的例子是用于治疗狭窄动脉的支架，可以将其与理想化血管进行对比。这种理想化血管可轻松建模为直管，显然应比带有支架的类似血管具有更低的血栓形成性。然而，这一假设无法通过简单的SA计算及两种情况下的SA分布比较直接得出。冠状动脉血流通过支架与直管对照情况的SA分布概率密度函数（PDF）比较（图2）说明了这一点。可以明显注意到，两种情况的SA分布主峰非常相似，但矛盾的是，在“尾部”较高的SA范围内，直管中出现更高SA值的概率反而更高。这种情况在SA值高于3.5帕·秒的范围内尤为显著（图2插图），而该范围有时被视为血小板活化的阈值[4,5,18]。

这个特定案例说明了在直管流动仿真这种简单情况下，此类偏差是如何产生的。显然，这是一个特殊情况，所展示的偏差不一定能在其他更复杂的心血管装置中重现。在具有理想层流的直管中，血小板沿其对应的层作直线流动，并受到恒定的应力作用。那些靠近壁面流动的血小板速度最慢，同时受到最高的壁面剪切应力，因此它们的溶血指数不成比例地升高远高于在较低剪切应力层中流动更快的其他血小板。这种偏差更加明显，因为在近壁区域存在更多的血小板，这是由于周长更长以及在近壁区域采用了边界层网格划分（见第3.2 小节）。当引入支架等装置时，这种现象会被消除，因为支架梁会引发紊乱的流动模式，使得血小板不再被捕获在同一层中。另一方面，在支架模型中，流经支架梁之间回流区的颗粒也会被捕获，并暴露于更高应力下且持续时间更长。然而，混合效应会将这些颗粒中的大部分释放到其他层中，从而降低其总溶血指数。该数值实验得出的两个实际结论是： （i）溶血指数可用于比较不同设计，但不一定适用于与健康对照组进行比较；（ii）出于比较目的的SA计算应仅包含装置本身的流动区域，而不包括进出口管道，因为后者可能会掩盖差异。其他影响此现象并可用于更好地分析结果的实施方案将在接下来的小节中介绍。

3.1 单次与重复通过

血小板沿不同持续时间或暴露时间以及长度的多种轨迹流经装置。SA计算可假设为单次通过，即每个血小板仅通过装置一次，或假设其反复流经装置，但在每次连续通过时从随机位置进入装置。在使用概率密度函数对装置中的SA进行比较时，每种方法都有其优缺点，因此应根据比较的目的来实施相应的方法。

在单次通过法中[9]每个粒子具有不同的暴露时间（$T$ 公式(2)），因此SA值能够更好地区分模式，并有助于识别区域升高的溶血指数。因此，该选项通常用于装置特定几何特征的局部效应研究的设计优化中。此方法的主要缺点是，不同暴露时间的颗粒之间的比较不仅对边界层中的颗粒（如直管情况所示）存在偏差，而且对暴露于较低应力的慢速颗粒也存在偏差。一种可能的解决方法是在仍仅比较单次通过的情况下，按时间对溶血指数进行平均（如第3.3小节所述）。

重复通道法 [19]在本文中的计算方法是通过对单次通过结果进行后处理，将流出的粒子随机重新播种来实现的。必须保持播种模式均匀，否则位置将无法完全随机化。例如，如果边界层中的播种模式更密集，则该区域重新播种的颗粒概率会更高。溶血指数 SA 的计算基于相同的总持续时间（T）：

$$
SA = \sum_{i=1}^{N} \left( \int_0^{T_i} \tau(t) dt \right)
$$

另一种选择是对所有轨迹使用相同的高重复次数，并假设随机化最终会使所有重复通过的总持续时间趋于相等。如果采用此假设，则必须检查在求解结束时域内是否存在滞留的血小板，因为这些血小板会在下一次通道开始时从其最终位置被“传送”回装置入口。这两种重复通过的方式都能更好地模拟体外再循环实验，其中重复通过是流经装置流动的固有特征。因此，重复通过更适用于仿真验证或装置性能评估。该方法的主要缺点在于，由于大量重复带来的随机化和累积平均过程

在重复通过的情况下，概率密度函数（PDF）分布趋近于特征性的高斯‘钟形曲线’分布。这使得难以区分装置特定设计特征可能产生的不同模式。溶血指数（SA）的数值显然高于单次通过的溶血指数，但它更能代表流环中全部血小板的血栓形成潜力。

概率密度函数(PDF)的比较 来自颗粒单次通过直管和冠状动脉支架的流经过程)

0 展示了通过直管和冠状动脉支架的层流之间的此类比较，其中溶血指数 SA 使用重复通道法计算。总持续时间设置为模拟持续时间 1.5 秒，并计算了 100,000 条可能的轨迹。与其他重复通过法溶血指数 SA [19],类似，概率密度函数 PDFs 的形状几乎呈高斯分布。尽管两种情况下的峰值看起来相似，但重复通过法 SA 确实表明支架带来了额外的溶血指数。在这种情况下，支架的重复通过法 SA 尾部持续达到 163 Pa·s，而在直管中仅达到 55 Pa·s。

3.2. 示踪粒子分布模式

如前所述，在近壁区域释放过多颗粒是不可取的，因为在直管中，这些颗粒对溶血指数 SA 的贡献会远超其他区域。当使用某些商业 CFD 求解器的默认选项时，这一问题尤为突出，因为在大多数情况下，颗粒会在网格节点上释放，导致壁面边界层中的播种密度更高。如果采用重复通过法（见第 3.1 节），这种播种模式可能会对任何类型的流动仿真产生偏差。以下通过三个相同直管仿真但不同示踪粒子分布模式的情况进一步说明了这种偏差。

本文介绍的全部三个模型均采用专门的示踪粒子分布模式，这些示踪粒子并非从网格节点释放。同心圆形种子模式（在图4中以蓝色标记）

是本文中使用的模式。该模式基于颗粒在径向和周向的均匀间距。笛卡尔播种模式（图 4中的红色）与同心模式具有相同的均匀间距，但定义在垂直方向和水平方向上。边界层种子模式（图4中的绿色）与同心圆形种子模式类似，但靠近壁面的两个环被五个边界层取代，这些边界层由1.2的增长因子定义。这三种情况的颗粒数量相近，约为11,500个。

图4的上图展示了这三种播种模式下溶血指数（SA）概率密度函数（PDF）的比较。如预期所示，边界层中溶血指数大于3.5 Pa·s的概率最高。两种均匀分布的示踪粒子分布模式的主要峰值相似，而显然在边界层播种模式下，血小板位于核心流中的概率较低。笛卡尔播种模式下溶血指数大于3.5 Pa·s的概率最低。这可能是因为靠近壁面的第一层血小板并非以等距环状排列，而是根据其笛卡尔坐标具有可变的距离。

3.3. 基于时间平均的轨迹比较

另一种在考虑不同轨迹时间长度的情况下比较其溶血指数的方法是对公式(2)中的积分进行时间平均：

$$
\overline{SA} = \frac{1}{T} \int_0^T \tau(t) dt \approx \frac{1}{T} \sum_{i=1}^{n} \tau_i \Delta t
$$

尽管有人认为该方法在概念上与重复通过法相似，但此方法可用于区分轨迹，并且无需使用概率密度函数。一个典型的例子是将溶血指数与Hellums等人[2]的血小板激活水平阈值进行比较（图5）。由于这些激活水平是在不同情况下测得的暴露时间和剪切应力水平，这种比较方法同时考虑了每颗血小板的轨迹持续时间和 $\overline{SA}$。最终得到的 $\overline{SA}$ 点云分布（直管 – 绿色，支架 – 红色）可以显示每颗血小板在各自轨迹结束时达到的SA值距离激活水平阈值有多远。对于几乎任何暴露时间，支架模型中的最大 $\overline{SA}$均大于直管模型。这些情况明显使血小板更接近Hellums等人的[2]激活水平阈值。与传统的PDF图相比，此类比较图形化地更显著地表示了“血栓形成潜力”，至少在此支架和直管模型的特殊情况下是如此。此外，与无法与任何激活阈值进行比较的SA不同，将 $\overline{SA}$和时间分离使得剪切应力和暴露时间的轨迹可以进行比较。这种方法的主要缺点是忽略了时间变化，但需要强调的是，这些时间变化在常规积分中也被忽略，只有在比较完整的应力‐时间加载波形时才能被考虑进去。

进一步对SA值进行时间平均，可以使其无量纲化。隆格斯特和克莱因斯特罗尔 [20] 提出了以下他们称之为血小板刺激历史（PSH）的无量纲SA：

$$
PSH = \frac{1}{\tau_w T_c} \int_0^T \tau(t) dt
$$

其中 $\tau_w$ 是血管的平均壁面剪切应力，$T_c$ 是输入波形周期。这两个参数代表的是整个流场，而非血小板本身及其轨迹，因此PSH相较于流体力学中的传统无量纲化，更接近于对溶血指数的归一化。一种可能的方法是为这种无量纲SA（NDSA）参数用于表示施加在血小板上的外部流动应力与血小板内部结构应力之间的比值：

$$
NDSA = \frac{1}{\sigma_c} \int_0^T \tau(t) dt
$$

显然，当前的拉格朗日方法无法直接计算 $\sigma_c$。未来可采用多尺度模拟来建模血小板的轨迹、其变形以及与周围血浆的流‐固相互作用[21] ，以用于计算 $\sigma_c$。作为一种近似方法，可以假设 $\sigma_c$ 与血小板的极限应力数量级相同。血小板的极限应力约为3000 帕，该值由未经处理的人类血小板的最大伸长长度[22]和线性应力应变曲线[23]的实验确定。

4. 结合湍流模型的拉格朗日方法

当需要对湍流进行数值 modeling 时，计算效率最高的方法是采用多个输运方程来近似雷诺应力的非定常雷诺平均纳维‐斯托克斯（URANS）。这种平均化可能会影响流动应力的计算方式（方程(4)）以及由此产生的拉格朗日轨迹。其他不基于平均和脉动分量来近似流动的湍流模型，对于拉格朗日损伤模型而言可能计算成本过高。直观上，雷诺应力应对溶血指数（SA）有非常显著的贡献。图6证实了在机械心脏瓣膜中湍流情况下的这一观点。应注意的是，使用雷诺应力预测血液损伤存在一定的争议，从认为应忽略雷诺应力的推测 [24]到明确表明雷诺应力与血细胞损伤相关[25,26]的实验测量结果不等。

在轨迹计算中考虑湍流脉动的一种可能方法是使用离散随机游走模型[27]来模拟与湍流涡旋的相互作用。该方法此前已被我们小组用于SA计算[8,10,28]。在此方法中，粒子与流体涡旋的相互作用可通过速度脉动进行估算，而速度脉动与局部湍流动能和涡旋寿命相关。涡旋寿命与拉格朗日积分时间尺度有关 。对于当前机械心脏瓣膜的情况，在入口处假设湍流强度为10%，入口处的比耗散率可近似为 ，时间尺度为 。图7比较了使用和不使用随机游走模型的SA概率密度函数。随机模型提高了SA水平，并降低了SA处于第一模式的概率，但两条PDF曲线较为接近。尽管随机模型为轨迹增加了一些“噪声”，但该过程也可能降低SA，因为它以类似于重复通道法中随机化的方式平滑了应力水平。当应用于MHV中的脉动流等情况时，随机模型可能发挥更大的作用。特别是在晚期收缩期，阀门关闭前减速阶段出现向湍流转变时，这种效应可能尤为显著。然而，这并非本文工作的重点，本文关注的是实现假设的一般性影响，而非MHV中脉动流等特定情况。

5. 耦合流体‐颗粒相互作用

根据所使用的数值求解器，有时无法采用双向颗粒‐血液相互作用耦合。代表流体颗粒轨迹的路径线可以在后处理中由每个时间步长的速度场计算得出，也可以在求解过程中通过单向流体到颗粒耦合进行计算。直观上，可能会认为这两种耦合方式对诸如代表血小板的微小中性浮力颗粒影响很小。本节通过比较双向DPM耦合与路径线（此处使用单向DPM耦合计算）来验证这一假设（图8）。这些情况还与最后时间步的瞬时迹线进行了比较。瞬时迹线（在某些商业CFD软件包中称为“pathlines”）是仅基于单个时间步速度场计算出的轨迹。显然，瞬时迹线的SA值与瞬态轨迹有很大不同，但在此展示是为了说明这种差异。完全耦合DPM的SA值最低，这可能是由于包含了将颗粒从高剪切应力区域移开的阻力所致。

专家评论

本文评估了不同实现假设对预测血液损伤的应力累积（SA）拉格朗日计算的影响。在本项特定的比较研究中，SA的计算基于线性SA叠加法则，但当采用其他任何幂律模型公式时（其中累积损伤由应力和时间的不同幂次乘积计算得出），预计也将得到类似的结果。实现假设包括不同的示踪粒子分布模式、随机游走模型，以及使用路径线的简化轨迹计算。评估的后处理实现选项包括单次通过和重复通过溶血指数，以及时间平均。

正如本研究明确展示的那样，实现假设会显著影响溶血指数（SA），即血液损伤模型的预测结果。部分结果在某种程度上出乎意料，例如显示出随机游走模型对SA的影响非常有限。显然，在使用针对特定案例或情景的拉格朗日模型时，应谨慎考虑。在更复杂的模型可能计算成本过高而对结果影响有限的情况下，可以采用一些简化假设。最终，应根据具体案例的物理特性和敏感性分析来确定适当的假设，类似于本文所呈现的案例研究中采用的方法。

本研究未包含一些额外的考虑因素。由于颗粒尺寸大于亚微米尺度[29], ，因此未考虑布朗运动，尽管仍可认为3 μm的颗粒尺寸接近该阈值，且在低雷诺数层流中布朗力可能变得显著。颗粒之间的碰撞也被忽略。其他对拉格朗日轨迹的影响可能源于血管运动。例如，被捕获在冠状动脉支架支架梁之间的颗粒可能在动脉收缩时被释放。若不采用更复杂的流固耦合（FSI）考虑或实验获得的数据，这种效应无法进行数值模拟，但它可能对概率密度函数（PDF）中溶血指数（SA）较高“尾部”范围产生显著影响，因为这些颗粒具有更长的暴露时间，并承受最高的应力水平。

五年展望

在接下来的五年中，随着高性能计算（HPC）资源的到来，使得在合理时间内求解更复杂模型成为可能，有若干关键方向值得探索。关于使用URANS以及如何正确考虑雷诺应力和随机运动的争议，可能会通过直接数值模拟（DNS）得以解决，其中纳维‐斯托克斯方程在无平均或湍流近似的情况下被求解。该方法要求网格分辨率达到微小的柯尔莫哥洛夫湍流能量级联尺度，并以相似的时间和空间分辨率计算拉格朗日轨迹。无网格模型，例如耗散粒子动力学（DPD），将整个连续体在分子尺度上离散为颗粒，并通过力定律定义颗粒间相互作用，也可能被采用。与DNS类似，这类模型无需对湍流进行特殊处理。此外，可以直接从颗粒计算拉格朗日轨迹，而无需引入额外相；如前所述，还可模拟血小板变形及其内部应力，用于计算无量纲SA（NDSA）。最后，目前可用于离散相模型的侵蚀模型可被修改以纳入血液损伤模型。在此类应用中，血液损伤将与CFD解耦合，而非在后处理中计算，从而使更复杂模型的应用成为可能。例如，若血液成分发生形变，可能会影响流动；或者颗粒可能分裂为多个更小的颗粒，正如溶血或血小板微颗粒脱落过程中所发生的那样。

关键问题

• 在将装置中单次通过的溶血指数与对照情况进行比较时，例如支架植入动脉与无支架动脉相比，由于靠近壁面的颗粒被捕获在速度更慢的边界层流动中并承受更高的剪切应力水平，导致对照情况下的溶血指数可能出现偏高。而在存在装置的情况下，紊乱的流动模式可防止这种偏差的发生。
• 对于直管中的层流情况，重复通过的溶血指数有助于消除溶血指数偏差，但对于优化装置设计以降低血栓形成性（即降低溶血指数）的效果有限。
• 播种模式可能加剧直管溶血指数的偏差，因为在高密度种子区域这种偏差会被放大。在缺乏流动混合的情况下，即使进行重复通过也不一定能消除该偏差。
• 对溶血指数进行时间平均化和无量纲化处理，可能有助于将具有不同暴露时间的路径线与完整的血小板活化水平轨迹范围进行比较。
• 引入随机游走模型以表征湍流效应对轨迹的影响，会使溶血指数水平相对于未考虑脉动的轨迹略有增加。
• 耦合离散相模型轨迹与路径线显示出相似的壁面剪切应力分布，但其幅度不同。
• 即使边界条件保持恒定，瞬时迹线也不能用于湍流条件下装置之间的比较。