动静区域耦合方法对血泵模拟结果的影响研究

林子祺， 于 勇， 张 宇

(1.北京理工大学 宇航学院流体力学实验室， 北京 100081；2.清华大学 医学院航天泰心科技有限公司人工心脏联合研究中心， 北京 100091)

引言

心室辅助装置(Ventricular Assist Device,VAD)以医工结合的方式帮助病人维持血液循环。对于VAD此类复杂泵体，计算流体力学(Computational Fluid Dynamics,CFD)能减少实验开销，实现对设备的评估和优化；但同时需要考虑CFD准确度的问题[1-3]。为方便进行CFD准确度的研究，美国食品与药品监管局(Food and Drug Administration, FDA)提出了“关键路径计划”(Critical Path Initiative, CPI)，发布了一个基准模型及围绕该模型获得的多个工况下的实验结果，以供与CFD仿真结果进行比对。2012年SANDY等[4]给出了该模型的详细几何参数。相较现有产品设计(如Impella、新一代HeartMate等)，标准血泵对红细胞的破坏性更强。但标准模型有详细的流动和溶血指数实验数据，是评价CFD仿真VAD准确性的合适选择[5]。

2017年，MALINAUSKAS等[6]总结了来自全球24个CFD仿真小组针对标准模型使用不同网格、湍流模型计算得出的结果，并与实验数据进行对比。所有小组对叶轮附近的流速仿真较为准确，但对泵出口扩散段的流速仿真结果较差。对于溶血估计的计算结果所有小组之间的差别很大，同时与实验差别也较大。

对常规工业泵的CFD模拟，所选湍流模型和定常/非定常方法的不同会影响CFD的结果[7-8]。在对标准血泵的模拟中，MALINAUSKAS[6]总结的24个小组中有14组明确使用定常仿真方法。近几年针对VAD的CFD仿真也大多使用定常方法，但是在出口扩散段流速与溶血率绝对值的估计准确度有待提升。使用非定常仿真一般能获得比定常仿真更准确的结果，但是计算量大，是否值得在工程上推广应用是一个值得讨论的问题。

本研究针对FDA标准血泵模型，对比了定常SRF和非定常滑移网格的仿真方法，评估了两种方法对速度、溶血率仿真准确性的影响，来讨论非定常计算的必要性。

1 材料与方法

1.1 FDA标模设计与工况

研究所用三维模型引自网站https://nciphub.org/wiki/FDA_CFD。转子直径52 mm，泵体内径60 mm，转子背面与泵体、转子叶刀顶部与泵体的间隔皆为1 mm。标准泵模型及对应的实验装置如图1所示。具体几何尺寸可下载网站中的Blood_Pump.zip文件。

图1 FDA发布的CPI实验装置

FDA为此标准血泵做了6个工况的实验[7]。本研究选用了其中实验数据完整的4个工况，工况参数见表1。研究所提及的“实验数据”皆指由FDA公布的流场或溶血数据。

表1 FDA血泵模拟工况参数

实验使用了粒子图像测速法(Particle Image Velocimetry, PIV)技术测量血泵内特定剖面的流场。实验用血为猪血液，密度1056 kg/m3，动力黏度0.0034 kg/(m·s)。实验时以雷诺数Re和流动系数φ为基准对实验时使用的流体流量与转速进行适当调整[9]。所用雷诺数的定义是

(1)

式中，ρ—— 血液密度

ω—— 转子转速

D—— 转子直径

μ—— 血液动力黏度

流动系数的定义是：

(2)

式中，Q—— 入口处血液质量流量

实验中定义坐标原点为转子无叶片一面的中心、z轴正方向与流体流入方向相反、x轴正方向与出口流体流动方向相同。同时定义了若干条采用线段(在使用PIV拍摄的剖面内)，采集线段上的流速，采样线段的命名沿用文献[7]，CS2，CS3，CS5分别表示不同PIV的拍摄平面，在本研究中B表示沿象限平分线的采样线段，L1～L4表示扩散段不同的采样线段。坐标轴、拍摄剖面示意及采样线段所在位置如图2所示。采样线段端点坐标见表2。

表2 仿真流速采集线段定义

图2 采样线段所在位置

1.2 流动基本方程

血液视为单相牛顿流体。数值求解RANS和URANS方程组。其中动量方程为：

(3)

式中，U—— 速度张量

τ—— 瞬时剪切应力

f—— 动力黏度附加力

湍流模型采用RNG k-epsilon 模型：

(4)

其中,C1=1.42,C2=1.68。

湍流黏性系数

(5)

其中，Cμ=0.0845。

在计算中，Y+为10左右，壁面采用增强壁面函数(Enhanced Wall Treatment)

(6)

仿真中模型入口直径为12 mm，入口边界条件见表3。出口使用自由流边界(outflow boundary)。

表3 模拟速度入口边界参数

1.3 流场仿真可信度分析

本研究使用了3个参数衡量仿真所得到流场与实验结果的偏差：确定系数D、极值误差α与极值点误差α0，其中确定系数越接近1表示仿真与实验的整体偏差越小，极值与极值点误差用以定量描述仿真与实验曲线在极点处的偏离程度。

(7)

(8)

(9)

式中，vi—— PIV方法所测一个空间点在xOy平面内的速度

vmax—— 使用PIV方法所测的特定线段上在xOy平面内的速度最大值

x0—— 取得vmax时特定线段上的位置参数

xmax,xmin—— 分别为采样线段坐标的最大值与最小值

1.4 溶血预测模型与方法

CFD用于评估溶血的绝对水平非常具有挑战性[10]。GIERSIPEN等[11]提出了一种基于应力的幂律模型，该方法将血浆游离血红蛋白的产生与剪切应力大小和暴露时间联系起来：

(10)

H是溶血率，即血浆中的自由血红蛋白和总的血红蛋白的百分比。Hct是血细胞比容，fHb是自由血浆中的血红蛋白(mg/dl plasma)，Hb是总的血红蛋白浓度(mg/dl blood)，C，a和b是经验系数。C的单位是Pa-1s-1，a和b是无量纲的，系数通常由在均匀剪切条件下在Couette型剪切装置中的溶血测量确定[12]。有许多实验团队以此方法针对不同的动物血液整定配套的实验参数C，a，b[13]。对于猪血，使用SONG等[14]整定的参数C=1.8×10-8，a=0.765，b=1.991。

借助离散相模型(Discrete Phase Model，DPM)对溶血率进行估计。对于计算已收敛的流场，在入口出撒入无质量的颗粒并进行追踪，统计其所经历的流场单元的剪切应力，给出溶血的预测值，即统计第p个粒子在第i个时间段内的溶血率：

(11)

随后再计算单个粒子在第i个时间段之前的累计溶血率：

Dp,i=Dp,i-1+(1-Dp,i-1)dp,i

(12)

最后，根据所有粒子的溶血率进行统计平均

(13)

其中，N为布置的DPM粒子个数，P为粒子编号。参考Richard的处理方式[7]，相关溶血系数将获得的绝对溶血率与工况C5的绝对溶血率相除，以呈现同一方法对不同工况溶血率变化趋势的估计能力。

1.5 网格及无关性验证

定常与非定常仿真使用相同的网格进行。网格采用组合的方式，由转子周围流体域及壳体壁面附近流体域组成，不同流体域之间使用interface进行流场信息交换，如图3所示。图3中外部流场网格以半透明灰色表示，近转子流场网格以显示面网格来表示，紫色部分为近转子网格的内部固壁边界(转子)。

图3 网格组合方式

对于定常仿真，对转子附近的流体域定义旋转坐标系；对非定常仿真，转子附近流体域有旋转运动，其他流体域静止，动、静区域通过interface进行耦合，在保证库朗数小于1的前提下选择合适的时间步长：为保证进行C6工况(流速最大)仿真时库朗数小于1，选取的时间步长为5e-5 s，在进行C1/C2/C5的仿真时沿用此时间步长。

进行网格无关性验证的不同网格属性及进行部分工况的定常计算结果见表4。

使用表4中的网格对C5工况也进行了仿真，获取表3中各个线段的流况数据并与实验数据进行对比，得到确定系数相差小于0.001。三种网格对工况C1/C2进出口压头的计算结果偏差在2%以内。特别地，根据FDA实验所获得的C1工况压头为172±12 mmHg，C2工况压头为353±18 mmHg，网格M2与M3获得的压头均在此范围内且十分接近。综合考虑计算资源的消耗与计算资源的提升，使用M2网格进行后续的计算。

表4 多种参与测试的网格属性

2 结果

2.1 流场

先后使用SRF与滑移网格方法对表2中4个工况进行仿真，在表3所示部分采样位置下分析流场速度，按照1.3可信度分析方法获得的结果见表5。

表5 仿真的定量比较记录

图4表示了针对采样线段CS2B，四种工况下分别使用SRF与滑移网格方法获得的仿真数据与实验数据比较结果。可以看出，滑移网格方法对实验结果的拟合比SRF方法更加准确。

图4 采样线段CS2B四个工况下定常SRF与非定常滑移网格仿真与CPI实验流速对比

图5表示了工况C2下，对出口扩散段的4条采样线段，使用SRF与滑移网格方法进行的仿真数据与实验数据比较结果。相比于泵内的流动，出口扩散段的仿真距离实验结果仍有一定差距。但是对工况C2的CS5L3/CS5L4两条采样线段而言，滑移网格的拟合效果有明显改善，但在别的工况、别的出口扩散段采样线段中，滑移网格方法并没有获得更接近实验数据的结果。

图5 工况C2出口扩散段采样线段CPI实验和仿真结果对比

2.2 溶血系数

图6所示为仿真获得的相关溶血系数σ与实验数据C的对比。前人仿真数据引自MALINAUSKAS的整理工作[6]。可以看出在使用基于应力的幂律模型预测溶血时，无论使用定常SRF还是非定常滑移网格方法计算都不能获得较准确的效果。

图6 相对溶血系数对比图

3 讨论

研究使用的网格及计算方法进行过网格无关性验证，在泵内流场(采样线段CS2B/CS3B)的仿真上对比实验数据的确定系数在0.9以上，C6工况有略低的确定系数，而滑移网格方法能获得比SRF方法更高确定系数的结果。

在出口扩散段流场(采样线段CS5L1～4)信息的仿真上，离出口越远，仿真数据与实验数据的差距越大，在C6工况差距最大，且在此处非定常滑移网格方法也不能保证比定常SRF方法更接近实验数据。对于溶血的预测亦然。针对溶血预测的准确性，WU等人认为基于应力的幂律模型对溶血的预测本身并不准确[17]，如果采用基于WU等人提出的能量耗散幂律模型会准确很多[18]，这值得进一步关注。从目前的结果看来，非定常滑移网格方法在溶血方面没有显著的优势可能与此有关。

在数值仿真中，也在同样的部位出现了与实验值差异较大与溶血预测不准确的问题[10]。且对于偏离设计的C6工况，各项参数也会出现数据离散[5]。如何把出口扩散段的流场信息仿真准确及获得准确的溶血预测结果，是接下来的工作难点。对于CFD仿真的准确度问题，使用定常方法还是非定常方法只是一个讨论的侧面。HUO等人也指出湍流模型、网格类型与密度、时间步长等因素对FDA血泵流场的预测结果都有很大影响[16]。在改变这些因素之后，在难以准确预测的区域非定常方法是否会有比定常方法更显著的优势，有待进一步研究。

SRF方法不改变网格结构，通过施加旋转坐标系来进行定常的迭代计算。此方法在模拟轴流泵(如HearMate II)时[15]，由于出口法向量与转子转轴同向、流体总是向着转轴运动，可以认为是与实际情况吻合的。但对于本研究中的FDA标准血泵，出口法向量垂直于转子转轴，SRF方法模拟的准确性尚待讨论。而滑移网格方法通过改变转子角位置进行非定常仿真，可以保证与实验实际情况吻合，但会消耗更多的计算资源。本研究亦表明，滑移网格方法确实在一定程度上获得了比SRF方法更准确的仿真结果。

4 结论

非定常滑移网格方法同定常SRF方法相比，结果更贴合实验结果。对于泵内流速场的求解，非定常方法能有更高的确定系数，尤其对于峰值速度、对取得峰值所在极值点的仿真比定常方法明显更加准确；对于出口扩散段的求解，两种方法对于越靠外侧的特征线段求解与实验结果差异越大，其中非定常方法在外侧的求解能优于定常结果，但仍然不能从本质上改善计算准确性。

综上所述，非定常求解结果在一定程度上会比定常求解更准确，但考虑到工程应用时的计算效率问题，尤其是在溶血估计上的不良表现，可以考虑选择定常方法在损失少量精度的情况下较快地获得计算结果。