基于多孔介质模型的膜式氧合器内部流场分析

2020-04-01叶非华廖虎易国斌

化工进展 2020年3期

叶非华，廖虎，易国斌

（1广东工业大学轻工化工学院，广东广州510006；2广东顺德工业设计研究院（广东顺德创新设计研究院），广东佛山528300；3武汉科技大学机械自动化学院，湖北武汉430081）

膜式氧合器是心肺手术中心肺机器的重要组成部分，可用于体外生命支持[1-2]。利用计算流体力学（CFD）对膜式氧合器中的血流进行数值模拟，可以预测其各项性能（例如压力分布、灌注动力学和气体传输性质等）[3-8]。由于氧合器纤维膜的不透明性，难以通过实验直接观察血液灌注动力学，相比于传统实验，CFD对氧合器性能的预测更为经济有效，是研究和优化氧合器设计的重要方法之一[9-10]。本文以一款市售国产膜式氧合器为研究对象，采用多孔介质模型，通过数值模拟氧合器内部流场特性，详细分析了流体的速度分布、压力分布、湍流强度分布等，并结合快速溶血数值预估模型计算得到了氧合器的标准溶血指数NIH。研究结果有利于产品开发人员了解膜式氧合器内部流体运动特性及其对产品性能的影响，为后续优化设计氧合器，进一步改善其性能提供一定的理论指导。

1 膜式氧合器工作原理

图1 所示为国产某型号膜式氧合器的工作原理。该膜式氧合器主要包含变温室、氧合室、中间管道，以及血液、温水和氧气的进出口等。两个腔室内部均由中空纤维膜填充，其中变温室用于对血液的温度进行调节，氧合室主要起排出血液二氧化碳，增加氧气含量的作用。膜式氧合器工作时，血液从变温室下底的进血口进入氧合器，由下而上从中空纤维管外通过变温室，与中空纤维膜管内从上往下流的温水在纤维表面进行热交换。变温后的血液通过中间管道进入氧合室，再由上往下从中空纤维膜管外通过氧合室，与中空纤维膜管内从上往下流的气体通过纤维壁上的膜孔进行气体交换。最后变温、氧合后的血液通过出血口离开氧合器。膜式氧合器避免了血液和水、氧气直接接触，生物相容性更好。

图1 膜式氧合器工作原理

2 数值模型

2.1 多孔介质模型

多孔介质模型是用于计算膜式氧合器内部流场的主要CFD 模型之一，即假设氧合器中的纤维束为多孔介质，通过渗透率和孔隙率估算纤维束对流体产生的动量损失，这种方法已被证明能够以有效的方式对氧合器性能进行模拟[11-15]。本研究采用Navier-Stokes方程和连续性方程描述笛卡尔坐标系中不可压缩黏性流体的稳态运动，如式(1)。

式中，ρ为密度；ui为流体速度；P为压力；τij为应力张量；gi为重力；S为源项。

多孔介质模型是在Navier-Stokes方程上叠加了一个动量源项S，来描述多孔介质区域中流体的动量损失，i方向上的源项为式(3)。

其中等式右边第一项代表黏性损失，第二项代表惯性损失。在流量Q＜6.00L/min时，惯性损失可以忽略不计[11,15-16]。因此多孔介质区域的动量损失仅取决于利用达西定律计算得到的黏性阻力，如式(4)。

式中，1/α为黏性阻力系数；A为横截面面积；ΔP为流体穿过纤维束的压降；μ为流体黏度；Q为流量；L为纤维束长度。

2.2 三维快速溶血数值预估模型

传统的溶血数值预估方法是基于拉格朗日粒子轨迹追踪的方法来进行溶血计算[17]，需要在设备入口处注入粒子并监测其在流场中的瞬时切应力和停留时间，累积计算沿着运动轨迹的所有时间步中的溶血值，最后求得总溶血系数[18]。由于多孔介质模型无法反应单个纤维周围的局部速度分布和切应力分布，使用拉格朗日法计算的溶血值会远大于实验值[19]。新型的溶血数值预估方法为三维快速溶血数值预估模型，这种模型考虑了流场整体平均效应，是对整个稳定流场的速度域与切应力域进行线性场积分计算[20-21]。

2.2.1 流体剪切应力计算模型

在湍流状态下，血液中的剪切应力包括由流体黏性引起的分子切应力和湍流引起的雷诺切应力。其张量表达式为式(5)，其中雷诺切应力张量为式(6)。

式中，μτ为湍流黏度；k为湍流动能；δij为Kronecker 数。利用米泽斯屈服准则（mises yield criterion）将剪切应力的张量形式简化为剪切应力标量形式，如式(7)。

2.2.2 溶血数值预估模型

Giersiepen 等[17]通过大量溶血实验，归纳了溶血值D与剪切应力τ以及暴露时间t的幂函数关系式，如式(8)。

式中，D为溶血值，是衡量溶血量的参数；ΔHb代表被损伤的血红蛋白浓度，g/L；Hb代表总血红蛋白浓度（Hb=140g/L）；τ为剪切应力；t为暴露时间；C、a、b是通过试验数据回归得到的经验常数。Giersiepen 等通过试验得到的一组常数值：C=3.62×10-7、a=2.416、b=0.785，被众多文献所引用，过往研究也证明了这组常数应用于溶血指数计算时的准确性。因此，在本研究中也采用这组常数数值进行溶血指数计算。

Garon 等[20]基于双曲输运方程开发了一种三维快速溶血数值模拟的方法来预测溶血。双曲输运方程如式(9)。

式中，V为速度矢量；DI为线性溶血指数；源项σ为单位时间溶血破坏率，具体计算如式(10)、式(11)。

在计算获得一个稳定流场后，膜式氧合器内部线性平均溶血指数为式(12)。

式中，Q为体积流量。通过指数换算得到溶血值，如式(13)。

将溶血值转换为最终的标准溶血参数值（NIH），如式(14)。

3 实验方案

通过实验测量膜式氧合器不同流量下变温室和氧合室的压降值，计算纤维束的黏性阻力。开展氧合器内部流场数值计算，验证模拟结果的准确性。实验流程见图2，本实验使用市售国产某型号成人膜式氧合器，以水作为流体介质，通过离心泵将水从进血口导入氧合器，测量离心泵在不同转速时氧合器进血口，出血口的流量值及压强值，其中流量值与压强值通过相应的传感器（速度传感器为SONOFLOW CO.56，压力传感器为NORa 有创压力传感器）测量得到。之后将氧合器的变温室和氧合室分离，并分别测量这两部分在不同流量下的压降值。

图2 实验流程

利用matlab软件将实验数据Q-ΔP线性插值拟合，如图3 所示，得到变温室和氧合室区域的ΔP/Q值，将其分别代入式(4)即可求得两区域各自的黏性阻力：氧合室黏性阻力1/α=9.210×108m-2；变温室黏性阻力1/α=1.069×109m-2。

图3 氧合室和变温室压降-流量拟合曲线

4 仿真计算设置

使用商业CAD 软件SolidWorks 2018 绘制膜式氧合器的流道三维模型。计算模型由两个计算域组成：环形纤维束区域（多孔介质域）和外部流体域（纤维束与壳体之间的间隙、入口管道、出口管道以及中间管道）。为让流动充分发展和避免回流，加长了出口管道，如图4所示，图5给出了氧合器各部分尺寸。

图4 膜式氧合器几何模型示意图

图5 氧合器结构尺寸（单位：mm）

由于流体域结构复杂，多孔介质域与流体域皆采用非结构网格划分，两计算域交界面设置为Interface。采用商业软件ICEM17.0 划分流场网格。为了验证网格无关性以及渗透率的正确性，分别以三种不同网格数（190万、296万、420万）计算入口流量为1.65～3.90L/min，流体介质为水时的进出口压降值，在低雷诺数（Re＜6000）下采用RNGk-ε湍流模型模拟，高雷诺数（Re＞6000）下使用标准k-ε湍流模型模拟，并与实验结果进行比较。

由图6可见，网格数逐渐增大对计算结果影响较小，296 万网格与420 万网格的计算结果基本一致，且与实验测试结果之间的误差较小，这证明了计算的准确性。综合考虑计算耗时和求解结果精度，在本论文中确定以296万网格数进行研究。

图6 网格无关性验证

图7 管道轴面速度云图

利用商业软件FLUENT17.0 进行流场模拟计算。入口采用速度入口边界条件（对应入口流量为1.65～6.00L/min），出口为静压0Pa 的压力出口，操作压力为默认的大气压101325Pa。湍流模型选择RNGk-ε模型，采用SIMPLE 算法进行压力-速度耦合，离散格式采用二阶迎风模式。流体介质血液被视为不可压缩的牛顿流体，具有恒定密度1055kg/m3和黏度0.0035kg/(m·s)，所有壁面均为无滑移固壁边界，将中空纤维膜作为各向同性多孔介质处理。当所有速度和湍流残差小于10-5时，计算视为收敛。

5 计算结果及分析

5.1 流场分析

为研究膜式氧合器内部流体运动特性，利用速度云图、矢量图、流线图进行同步分析。图7显示了入口流量为4.50L/min 时的速度分布云图，其中A、B、C三个截面分别为膜式氧合器入口管道、中间管道和出口管道轴中间面。图8为全局速度流线图。

图8 全局速度流线图

从图7 和图8 可见，由于入口和流动腔室设计，流体在进入变温室的纤维束之前，轴向射流受到纤维束区域的流动阻力，先分流到流体域中，之后流体在膜前后压差的作用下突破流动阻力径向穿过纤维束区域，向上汇集到中间管道，由中间管道进入氧合室，最后在出口处重新汇集流出氧合器。提取不同流量下中间管道和出口管道的轴切面平均速度变化，如图9和图10所示，在中间管道起点处由于流体的汇集，流速激增。随着流体的继续流动，平均流速逐渐降低，直至到达管道中部流速才保持稳定，但在管道终点处由于流体的分流，使得切面平均流速有所增加。出口管道流速分布与中间管道类似。

图9 中间管道轴切面平均速度变化

图10 出口管道轴切面平均速度变化

为对比不同流量下流场分布情况，进一步研究了膜式氧合器径向中间截面的速度云图和出口管道区域的速度矢量图。如图11 所示，流量在1.65～6.00L/min 时，随着入口流量的增加，膜式氧合器内部速度梯度分布形式基本保持不变。从图12 可见，在出口管道上壁面出现较大涡旋，其涡流区约占流道的1/6～1/3，过大的涡旋产生了血液再循环流动现象，会导致高剪切应力的产生，同时增加了血红细胞在氧合器中的暴露时间，加剧了血红细胞破坏的可能性。随着流量的增加，涡旋区域面积逐渐增大，且涡旋中心逐渐向出口方向移动，回流区域面积增大，使得倒流量逐渐增加。

图11 氧合器径向中间截面的速度云图

图12 出口管道速度矢量图

图13 氧合器压力分布

图13显示了入口流速为0.955m/s（流量4.50L/min）时全局和径向中间截面的总压云图。由于入口管道、中间管道和出口管道的非轴对称设计，膜式氧合器内部压力分布不同于Gage 等[11]和Pelosi 等[16]预测的同心均匀模式，而是呈倾斜状态逐渐减小。出入口总压降ΔP 约为36677Pa，由于纤维束的高黏性阻力，大部分压力损失位于纤维束内，其中53.3%位于氧合室，42.6%位于变温室。入口管道和中间管道对压降的影响可忽略不计，出口管道提供的压力损失约占4.1%，且由于血液进入出口管道时流通面积有较大变化，在出口管道内产生了负压区，负压区与涡旋区一致。

图14 壁面剪切应力分布云图

5.2 血液损伤分析

本文将三维快速溶血数值预估模型通过用户自定义函数（UDF）编译到FLUENT 17.0中，在获得稳定的流场后计算膜式氧合器的标准溶血指标，同时将壁面剪切应力、湍流强度和红细胞停留时间作为辅助参数，用来分析氧合器中血液的损伤程度。

如图14 所示，流量为4.50L/min 时的最大壁面剪切应力约为118Pa，位于氧合室入口流体域外壁面。高壁面剪切应力区域存在于氧合室和变温室的入口及出口部分，与高流速区对应，而其他流域的剪切应力较低，不太可能产生过量流动引起血液损伤。图15 显示流体域中湍流强度较大，主要出现在分流或回流处，而多孔介质域中也存在小区域的湍流强度较大区，位于与流体域的相连之处。氧合器纤维束内的血流速度较小，施加在血液上的剪切应力低，对血细胞破坏小。而氧合器流体域，尤其是氧合室和变温室的入口及出口部分，皆出现较高的速度、壁面剪切应力和湍流强度分布，这些区域是血细胞容易遭到破坏的高发区域。

图15 湍流强度分布云图

图16 流量1.65～6.00L/min时红细胞停留时间

图17 较大流量时红细胞停留时间分布

计算收敛后，将直径为7μm的球形颗粒代表红细胞，利用离散相模型（DPM）计算粒子从氧合器入口到出口所经历的时间。图16 显示随着流量的增加，红细胞的停留时间及其标准差逐渐减小。图17 显示在较大流量时，红细胞在氧合器中的停留时间分布趋于平稳。当流量为3.50L/min 时，红细胞平均停留时间依然达到9.13s，停留时间过长。这是由于该膜式氧合器采用的是分离式单通道结构设计，血液预存量大，通过氧合器流程长。此外流场中形成的涡旋、回流等不规则流动现象也是造成红细胞停留时间变长的原因。红细胞在氧合器停留时间长，血液损伤的风险高，不利于氧合器的长期使用。