1. 清华大学 能源与动力工程系, 北京 100084;
2. 清华大学 第一附属医院, 北京 100016;
3. 清华大学 临床医学院, 心血管组织工程与生物材料实验室, 北京 100084;
4. 广州国家实验室, 广州 510005
收稿日期：2022-09-09
基金项目：广州实验室应急攻关项目(EKPG21-15)
作者简介：简萌(1999—), 女, 博士研究生
通讯作者：罗先武, 教授, E-mail: luoxw@tsinghua.edu.cn
摘要：为分析体外膜肺氧合技术(extracorporeal membrane oxygenation, ECMO)氧合器内膜丝排列方式对血液流动和气体输运的影响, 该文选取某商用氧合器的轴截面为简化模型, 采用基于沉浸边界(immersed boundary, IB)法的自编软件进行了稳态血液流动与多气体组分耦合输运的数值模拟。通过与商业计算流体动力学(computational fluid dynamics, CFD)软件Fluent计算结果、已有文献结果进行比较, 验证了所采用的数值方法和自编软件具有适用性。计算结果表明：在一定的孔隙率下改变膜丝排列方式和膜丝倾斜角将影响轴向和径向渗透率, 会明显改变血液流动和气体输运的性能; 具有径向交错和大倾斜角的膜丝排列能避免出口段大尺度旋涡, 降低血液损伤风险。进一步分析可知, ECMO气体输运性能主要受膜丝排列方式的影响, 且增大倾斜角有利于提升交错阵列的气体输运效率。该文研究可为实际工程中广泛存在的三维孔隙尺度多相流动数值模拟提供一定参考, 为ECMO结构优化、临床运行中参数调整提供科学依据。
关键词：体外膜肺氧合技术(ECMO)膜丝阵列孔隙尺度沉浸边界法多组分气体耦合求解器
Numerical simulation and analysis of multiphase flow through fiber array structure in extracorporeal membrane oxygenation
JIAN Meng1, ZHANG Mingkui2,3, HUANG Jianbing4, LUO Xianwu1,3
1. Department of Energy and Power Engineering, Tsinghua University, Beijing 100084, China;
2. First Hospital of Tsinghua University, Beijing 100016, China;
3. Cardiovascular Tissue Engineering and Biomaterials Laboratory, School of Clinical Medicine, Tsinghua University, Beijing 100084, China;
4. Guangzhou Laboratory, Guangzhou 510005, China
Abstract: [Objective] Extracorporeal membrane oxygenation (ECMO) is an effective life support treatment for severe cardiopulmonary failure and is widely used as supportive therapy for COVID-19. The design and assessment of ECMO oxygenators, which consist of thousands of 3D hollow fiber bundles, are essential to expanding their clinical applications. This study aims to investigate the effects of 3D fiber membrane array structure on the hemodynamic loss and gas transfer efficiency in ECMO. An axial slice of a commercial oxygenator is selected as a simplified model. [Methods] The immersed boundary (IB) method code was developed to simulate the two-dimensional steady laminar flow, and a segregated solver was implemented during the coupled multi-component gas transfer in ECMO. A grid independence test was carried out to ensure that the computational results were not influenced by the grid size. Results obtained from the IB method, commercial computational fluid dynamics (CFD) software Fluent with a body-fitted mesh, and the reference showed a good agreement, validating the accuracy of the IB method and the gas transfer solver. Seven array arrangement schemes with constant porosity were simulated at Re=5, and the results of permeance, wall shear stress, vortex distribution, and entropy generation rate were compared. [Results] Numerical results showed that when porosity was constant, different fiber array arrangements and angles between the odd and even row fibers could significantly change the flow state and gas transfer performance by affecting the relative value of axial and radial permeability. Inline arrangements and small angles between the odd and even row fiber arrangements deteriorated the uniformity of the flow state, consequently enlarging flow separation zones and causing peak wall shear stress. The array staggering from the axial direction and large angles between the odd and even rows could be used to avoid the large-scale vortex at the outlet of the ECMO and reduce the risk of blood damage. For all fiber array configurations, the head loss values predicted by entropy generation theory were smaller than the calculated results. As for gas transfer, in regions near the oxygenator inlet, outlet, and fluid retention zones, gas was mainly transferred through diffusion, controlled by the concentration gradient. In the middle stream regions, convection dominated the gas transfer. Further analysis showed that the gas transfer performance was mainly affected by the array arrangement. Compared with inline arrays, staggered arrangements increased the gas transfer rate on the upstream side and gaps of adjacent fibers. For the staggered arrangement where the flow passage had no evident periodic contraction and expansion, the flow retention area was small. Thus, the effective oxygenation area was the largest, and the gas was transported mainly by convection, which resulted in a high gas transfer rate. In addition, the gas transfer efficiency improved by increasing the angle between the odd and even rows. [Conclusions] The simulated pressure drop results of the simplified model are within the clinical operating range, but its gas transfer rate is lower than that of commercial oxygenators. Our results can serve as a reference for further 3D pore-scale numerical simulation and as a scientific basis for the structural optimization and clinical applications of ECMO.
Key words: extracorporeal membrane oxygenation (ECMO)fiber array arrangementpore-scaleimmersed boundary methodmulti-component gas transport coupling solver
体外膜肺氧合技术(extracorporeal membrane oxygenation, ECMO)是一种在体外通过对静脉血中氧气(O2)和二氧化碳(CO2)进行更新，从而实现短期替代人体肺功能的呼吸支持技术，在新型冠状病毒感染等导致的急性呼吸窘迫综合征治疗，肺移植、心脏瓣膜等心血管外科手术中发挥着重要作用[1]。ECMO氧合器需要满足血液灌注量小、压降低、气体交换效率高、生物相容性好等指标要求。因此，如何在降低血液损伤的前提下提升气体交换效率、延长呼吸支持周期是ECMO氧合器设计和优化的主要目标。目前通用的氧合器采用中空纤维膜结构，膜丝涂层改性、组件的血流动力学优化是提高氧合器生物相容性、优化血流模式、降低传质阻力的重要手段。Tabesh等[2]提出采用膜丝装填密度、O2输运量和压降损失率等指标评价氧合器性能，并通过实验对比了6款商用ECMO设备性能。
与临床实验相比，计算流体动力学(computational fluid dynamics, CFD) 可快速准确获得不同膜丝结构的流场参数，是研究ECMO内流动的有力工具[3]，精细化的流动模拟与流场诊断技术可揭示ECMO微观流动和宏观外特性之间的关系，为ECMO样机设计优化提供科学依据，提高改型优化的准确性[4]。由于商用ECMO设备膜丝外径为300~380 μm，膜表面积达1.2~2.5 m2，且排列紧密[5]，传统的计算方法难以满足贴体边界层网格数量与质量的要求。目前ECMO氧合器流动模拟主要采用多孔介质模型[6]和微尺度周期阵列模型[7]等简化方法，但各向同性多孔介质模型中血液均匀流动的假设与膜丝周围血流运行工况不符[8]，需要进行修正才能将计算误差控制在允许范围内[3]。
为表示不同膜丝排列下多孔介质的流动阻力，需要对表征单元体积(representative element volume)进行微尺度数值模拟。王东璞等[9]通过对有限尺寸固体颗粒阵列孔隙尺度直接进行数值模拟，研究了多孔介质结构与流动的相互作用对传热特性、流体质点输运的影响。Low等[10]通过微尺度二维数值模拟，获得了线性和交错周期阵列下血液沿轴向、横向流动时渗透率和气体扩散系数与孔隙率的关系，为预测实际情况下ECMO的压降和氧合效率提供了参考。文[11]利用沉浸边界(immersed boundary，IB) 法对二维随机和周期排列的圆形、方形和椭圆阵列绕流进行稳态模拟，得出了渗透率随孔隙率、Reynolds数、截面形状和长宽比的变化关系。
ECMO氧合器气体输运包括静脉血的O2富集和CO2清除，膜丝空腔内的气体在浓度差作用下在膜侧区域发生Knudsen扩散；在膜外壁血液区以溶解组分、血红蛋白结合组分、反应物组分[12]完成对流扩散输运[13]，其中血红蛋白与O2和CO2的亲和曲线具有Haldane-Bohr耦合效应[14]。该过程涉及多区域耦合及大量生化反应。Zierenberg等[15]通过血液侧O2和CO2耦合输运的数值模拟获得了不同二维阵列的瞬态气体输运特性，发现血液损伤风险与气体输运效率存在互相制约的关系。Taskin等[16]提出了一种多区域O2输运模型，同时考虑血液侧、膜侧传质阻力，能减小传质效率的过度预测。Harasek等[17]构建了耦合血液侧、纤维膜侧和纤维膜内腔的全区域CO2输运模型，CO2浓度的数值计算结果与实验相符。
综上所述，目前微尺度模型研究的膜丝阵列多为排列简单的二维圆柱阵列，其渗透率和气体输运的结果能否应用于实际复杂模型还需进一步验证；O2和CO2多区域耦合的输运研究较少，尚不能揭示血液流动与气体输运的联系。
考虑到ECMO氧合器具有复杂边界且流动和传质受边界层影响显著的特征，本文基于IB法高效构建三维周期性膜丝壁面边界和二维计算网格，以减少采用传统贴体非结构网格在网格生成中的时间成本，实现流动计算中近壁面网格与离散误差的精准控制，通过模拟某商用ECMO氧合器轴截面内二维流动，验证IB方法在模拟复杂二维流场的准确性；采用自编多区域、多组分气体输运分离式求解软件，分析膜丝排列及相邻层交叉角对流动、能量损失和气体输运的影响，并与实际ECMO运行参数对比，为ECMO氧合器设计和评价提供科学参考。
1 数值方法1.1 基于IB法的流动模拟本文将ECMO氧合器膜丝阵列的低Reynolds数流动视为二维稳态不可压缩层流，通过SIMPLE算法求解耦合的速度和压力场。血液流动的质量与动量守恒方程可表示如下：
| $\begin{gathered}\nabla \cdot \boldsymbol{u}=0, \end{gathered}$ | (1) |
| $\nabla \cdot(\boldsymbol{u} \boldsymbol{u})=-\frac{\nabla p_{\mathrm{b}}}{\rho}+\frac{\mu}{\rho} \nabla^{2} \boldsymbol{u} .$ | (2) |
ECMO氧合器内部流动满足Darcy定律，通常采用无量纲渗透率系数kp来描述多孔介质对流动的阻碍作用，具体可表示为
| $k_{\mathrm{p}}=\frac{\mu U f}{\nabla p_{\mathrm{b}}} \cdot \frac{A_{0}^{2}}{\left(4 {\rm{ \mathsf{ π} }} e_{0} b_{0}\right)^{2}} .$ | (3) |
计算中，通过设定流体域边界质量不平衡残差、流体网格速度和压力的累积残差阈值来控制迭代次数，保证计算结果收敛。
为求解上述基本控制方程，本文利用FOAM-extend 4.0平台自编了基于清晰界面IB方法的流动模拟软件。其中，边界通过三角形贴体网格离散，背景网格采用非结构化同位Cartesian网格。以图 1中二维单膜丝计算域为例，计算时背景网格被划分为流体网格、IB网格和固体网格，图中分别用圆形、三角形和菱形表示，其中流体网格和IB网格根据区域可进一步分成膜侧和血液侧网格，分别用蓝色和红色进行标记。固体网格位于膜丝内壁内侧，膜侧网格位于膜丝内外壁面之间。膜侧网格中与固体网格和血液侧网格存在交界面的为膜侧IB网格，其余为膜侧流体网格。血液侧网格位于膜丝外壁的外部，其中与膜侧IB网格存在交界面的为血液侧IB网格，其余网格为血液侧流体网格。图中，IB网格的边界面法线与沉浸边界的交点为IB点，记为ibp，血液侧ibp和膜侧ibp通常不重合，分别用红色和蓝色星形表示。为了计算边界物理量的梯度，每个膜侧ibp在面法向上设置1个采样点，对应图中红色矩形处。仅对流体网格进行离散化计算，固体网格的速度和压力可以进行预设。IB网格的速度和压力通过边界点和影响域内流体点经过加权的最小二乘法进行插值预估，以满足无滑移、无梯度等边界条件。流动计算时，膜侧气体扩散速率可忽略不计[16]，将膜侧流体网格速度预设为0。通过修正IB与流体网格交界面的流量通量实现流场迭代计算。
 |
| 图 1 二维单膜丝计算区域及网格标记 |
| 图选项 |
1.2 Dirichlet-Neumann分离式气体输运求解方法当膜丝孔隙不发生润湿时，血液侧液相传质阻力和非湿润膜侧传质阻力远大于膜内侧气相传质阻力[18]，计算时可忽略膜内侧传质阻力，则血液侧和膜侧气体输运的控制方程可分别表示为：
| $\nabla \cdot\left[\boldsymbol{u}\left(\alpha_{i} p_{i}+H_{\mathrm{b}, i} S_{i}\right)\right]=\nabla \cdot\left(D_{\mathrm{eff}, i} \nabla\left(\alpha_{i} p_{i}\right)\right)+ \\\;\;\;\;\;\;\;\;\nabla \cdot\left(D_{\mathrm{rbc}} \nabla\left(H_{\mathrm{b}, i} S_{i}\right)\right), $ | (4) |
| $\frac{\partial p_{i}}{\partial t}=D_{\mathrm{p}, i} \nabla^{2} p_{i} .$ | (5) |
本文使用分离式求解方法实现膜丝与血液交界面耦合。计算血液侧时，交界面为固定值边界，使用上一次迭代后膜丝外壁的气体分压值，此时第j个IB网格的物理量?j插值预估满足如下公式：
| $\begin{gathered}\phi_{j}=\left(\phi_{\mathrm{ibp}}+a_{1} X_{j}+a_{2} Y_{j}+a_{3} X_{j} Y_{j}+\right. \\\left.a_{4} X_{j}^{2}+a_{5} Y_{j}^{2}\right) \omega_{j} .\end{gathered}$ | (6) |
计算膜丝侧时，血液侧气体分压数值不变，交界面为耦合通量边界，保证两个区域与交界面的气体扩散通量相等。耦合通量边界条件中?j插值预估满足如下公式：
| $\begin{gathered}\phi_j=\left(c_0+c_1 X_j+c_2 Y_j+c_3 X_j Y_j+\right. \\\left.c_4 X_j^2+c_5 Y_j^2\right) \omega_j, \end{gathered}$ | (7) |
| $\left.D_{\mathrm{p}, i} \frac{\partial \phi}{\partial \boldsymbol{n}}\right|_{\mathrm{ibp}}=D_{\mathrm{eff}, i} \frac{\phi_{j, \mathrm{f}}-\phi_{\mathrm{ibp}}}{d} .$ | (8) |
血红蛋白的携O2、携CO2能力存在耦合效应，2种气体仅在血液侧耦合输运，在膜侧以相互独立的气体组分存在。暂不考虑pH值、温度、化学物质以及各组分化学反应对输运的影响。O2和CO2的饱和度根据文[19]的公式来描述，具体可表示为：
| $S_i=\frac{\left(p_i\right)^{n_i}}{\left(p_{50, i}\right)^{n_i}+\left(p_i\right)^{n_i}}=\frac{\left(p_i / p_{50, i}\right)^{n_i}}{1+\left(p_i / p_{50, i}\right)^{n_i}}$ | (9) |
| $p_{50, i}=\beta_{1, i}\left(\frac{1+\beta_{2, i} p_k}{1+\beta_{3, { }_i} p_k}\right), \quad i \neq k \text {. }$ | (10) |
计算血液侧两种气体组分耦合输运时，采用分离式求解方法，即先计算某组分的对流扩散方程，收敛到一定残差后再对另一组分进行更新。采用局部壁面的Sherwood数评价对流输运与扩散输运的相对强度，i气体组分Sherwood数可表示为
| $S h_i=\frac{\left.r \alpha_i\left(\boldsymbol{n} \cdot \nabla p_i\right)\right|_{\text {wall }}}{C_{\text {wall, } i}-C_{\text {in }, i}} .$ | (11) |
2 计算模型2.1 计算域本文研究对象为某商用ECMO氧合器模型，其氧合区为同心圆筒状结构，膜丝轴切面为椭圆环，不同轴面上的孔隙率保持不变，但排列不同。为比较不同轴截面及缠绕角度下的流动，假设氧合区的周向流动可以忽略，只选取氧合区的一个轴截面进行二维定常流动模拟。为避免在ECMO的进出口出现回流，选定计算域时在膜丝阵列的进、出口分别延长了30、30.3 mm。图 2所示为计算域，其中：沿x方向为径向，沿y方向为轴向，沿z方向为垂直轴截面方向，血液从最左侧断面流入，从最右侧断面流出。为便于分析，选取Ⅰ—Ⅳ共4个代表性区域，其中：Ⅰ位于膜丝区入口近轴线的内侧圆角处，Ⅱ位于膜丝区中间段的局部区域，Ⅲ位于膜丝区出口远离轴线的外侧圆角处，Ⅳ位于膜丝区出口延长段区域。此外，还选取了包含膜丝区整个径向长度的区域A、B和C，其中：A、B、C分别位于膜丝区上壁面、中段和下壁面处。
 |
| 图 2 计算域(单位：mm) |
| 图选项 |
2.2 ECMO膜丝阵列图 3所示为ECMO内2根交错排布膜丝的空间关系，其中一根膜丝的中心线与氧合器轴线平行，膜丝以蓝色线表示；另一根膜丝的中心线与氧合器轴线不平行，膜丝以红色线表示，而红色椭圆环表示该膜丝在轴截面上的形状。两膜丝中心线位于距氧合器轴线为R的圆柱平面上，膜丝中心线与氧合器轴线之间的倾斜角为θ(以逆时针方向为正)。
 |
| 图 3 2根交错排布的膜丝 |
| 图选项 |
图 4所示为ECMO膜丝阵列的7种排布方式，图中未表示膜丝内侧的结构。图 4a为线性阵列，所有膜丝的倾斜角均相等，θ0=52°。膜丝阵列的列间距lx=0.36 mm，行间距ly=0.788 1 mm。图 4b为交错阵列X1。将线性阵列中偶数列的膜丝倾斜角变为θ0，并改变两相邻交错膜丝的轴向间距。图中，交错膜丝中心距离a=0.281 mm。图 4c为交错阵列X2。在交错阵列X1的基础上，改变两相邻交错膜丝的轴向间距。图中，交错膜丝中心距b=0.188 mm。图 4d为交错阵列X3。在X1的基础上，改变两相邻交错膜丝的轴向间距和径向间距。图中，交错膜丝中心距c=0.25 mm。图 4e为交错阵列X4，在X1的基础上，改变膜丝倾斜角，θ1=15°。此时，为保证膜丝阵列孔隙率不变，将行间距改为l1=2.028 mm。图 4f为交错阵列X5，在X1的基础上，改变膜丝倾斜角，θ2=30°。此时，为保证膜丝阵列孔隙率不变，将行间距改为l2=1.050 mm。图 4g为交错阵列X6，在X1的基础上，改变膜丝倾斜角，θ3=60°。此时，为保证膜丝阵列孔隙率不变，将行间距改为l3=0.606 mm。
 |
| 图 4 7种不同的膜丝阵列 |
| 图选项 |
图 4中，阵列X1、X2和X3为孔隙率一定时不同膜丝间隙下的交错阵列；X1、X4、X5和X6为孔隙率一定时不同膜丝倾斜角下的交错阵列。
膜丝外径r=0.3 mm，内径r0=0.24 mm，各阵列孔隙率ε=0.683。在膜丝内的O2、CO2经由膜丝壁面与血液进行质量交换。
2.3 计算条件与网格血液的流变呈现黏弹性和剪切稀化特征。Taskin等[16]通过计算氧合器二维阵列O2输运结果对比了改进Casson模型[19]和考虑血液黏弹性的Oldroyd模型[20]的预测效果，发现Re=5时2种模型的壁面O2分压的预测相对误差在10%以内，血液黏弹性对气体输运的影响可以忽略不计。Karimi等[21]通过数值计算对比了9种血液流变模型对主动脉内血流动力学的预测结果，其中Carreau模型和改进Casson模型在较低和较高的切应力水平下与血液的实际流变特点相符。本文采用交错阵列X5在Re=5的条件下比较Carreau模型和改进Casson模型对流动和输运的预测效果，2种模型预测的渗透率相对误差为0.135%。图 5展示了2种流变模型O2分压沿轴线分布的预测结果，其中L表示轴向坐标，膜丝区入口位于L=-0.031 m处。2种模型计算的O2分压沿轴线的累积相对误差为7.7%，本文数值模拟采用改进Casson模型，黏度μ与剪切率φ的关系可表示为
 |
| 图 5 X5阵列2种流变模型下O2分压沿轴向分布对比 |
| 图选项 |
| $\mu=\left(\left(\frac{\tau_{0}}{\varphi+1}\right)^{1 / 2}+\mu_{c}^{1 / 2}\right)^{2}$ | (12) |
血液流动计算时, 人口边界条件为恒定速度, 满足
计算域的背景网格为正交网格，膜丝采用三角面网格表示，依据不同缠绕方式所得阵列的膜丝数量为4 509~8 172根。为确保网格分辨率，对膜丝浸没区域2个坐标方向进行等比例局部加密，并进行网格无关性验证。表 1显示了在不同网格下，Re=5、θ=45°、ε=0.735时无量纲渗透率系数kp、壁面平均切应力τave的计算结果。当网格节点数大于1.1×107时，采用IB方法预测的kp、τave均变化较小，因此选取1.1×107网格进行数值模拟。
表 1 网格无关性验证
| 网格节点数/107 | kp/10-2 | τave/Pa |
| 0.900 | 1.23 | 7.83 |
| 1.100 | 1.18 | 8.29 |
| 1.300 | 1.13 | 8.52 |
表选项
3 结果分析3.1 模型验证为了验证IB方法在计算复杂多孔介质微间隙流动的准确性，本文对自编软件和商业CFD软件Fluent的模拟结果进行了对比。需要说明的是，Fluent模拟的计算域、边界条件与IB方法相同，网格数为0.689×107，计算得到kp=1.13×10-2，2种方法预测kp的相对误差为3.5%，说明本文数值模拟方法具有较好的适应性。
图 6为Re分别为5、20、40时，采用IB方法、商业软件Fluent以及文[16]计算的O2在血液和单膜丝区域内输运的结果，图中数据表示膜丝外壁O2分压随周向角γ分布，膜丝前缘点为γ=0°，后缘点为γ=180°。通过比较可知，当Re=5时，O2分压处于710~750 mmHg范围内，且在0°~180°范围内单调增加，3组数据具有较好的一致性。当Re=20时，3组数据均显示在140°附近出现最大的O2分压。此时发生分离流动，并使角度在进一步增大时的对流传质减弱，O2分压下降。当Re=40时，IB的计算结果与Fluent的模拟结果均表明绕膜丝壁面的O2流动在60°附近发生分离，然后在130°附近再次发生分离，在第2次分离点处O2分压出现极大值，之后随着角度增大O2分压下降。图中结果验证了本文多区域O2输运计算方法的准确性。
 |
| 图 6 单膜丝外壁面的O2分压分布 |
| 图选项 |
为进一步验证多组分气体输运效果，模拟了O2和CO2在血液侧的耦合输运。图 7为Re=10、ε=0.3时二维交错膜丝阵列中2种气体组分沿膜丝壁面的Sherwood数随周向角分布情况，计算结果与文[10]数据吻合，且均表明在0°~180°范围内出现2个明显的对流传质峰值。由此可知，本文所采用的数值方法在模拟膜丝阵列多相流动方面具有良好的适应性。
 |
| 图 7 二维交错膜丝阵列壁面气体Sherwood数(Sh)分布 |
| 图选项 |
3.2 不同膜丝阵列的血流特性图 8展示了采用不同阵列时沿氧合器轴向的血液压力分布。图中，曲线斜率绝对值越大，所对应的渗透率越低。结果表明：线性阵列的进出口压降高达7×105 Pa，远高于其余阵列；对比交错阵列X1和X2可知，随着轴向间隙增加，渗透率增大；对比交错阵列X1和X3可知，相邻膜丝间设置径向交错间隙较轴向交错间隙更利于血液通过。对比膜丝倾斜角不同时的交错阵列X1、X4、X5和X6可知，膜丝倾斜角为15°时(即X4)压降较大，而其他阵列的压降均较小。因此，当膜丝阵列的截面孔隙率和列间距保持一定时，随着倾斜角增大，膜丝阵列行间距减小，渗透率呈现先增大后减小的趋势。
 |
| 图 8 不同膜丝阵列方式下沿轴向的血液压力分布 |
| 图选项 |
图 9展示了线性阵列与交错阵列X1、X2和X3的流速分布及速度矢量图。对于每种排列方式，从左至右分别对应图 2Ⅰ—Ⅳ局部区域的流场。由图中结果可知，在膜丝区进口区域Ⅰ、出口区域Ⅲ和Ⅳ的流动非均匀性较强，而中间区域Ⅱ则可视为顺轴向的周期性绕流。
 |
| 图 9 4种膜丝阵列方式在不同局部区域的流速分布及速度矢量 |
| 图选项 |
对于线性阵列，在入口区域Ⅰ处的轴向流动阻力大，血液接触膜丝后主要体现为径向流动，并伴随绕椭圆阵列的流动；在区域Ⅱ处的间隙流速高，膜丝前缘与尾缘附近存在明显分离流动；而在出口区域Ⅲ和Ⅳ处，流线在圆角外侧汇聚，在圆角区下方存在较大的低速区。
对于交错阵列X1、X2和X3，流动特征如下：
1) 在区域Ⅰ处，阵列X2的流动冲击最大，阵列X3的流动最均匀，但阵列X3的轴向渗透率显著高于径向。
2) 在区域Ⅱ处，血液顺轴向绕流膜丝时，在膜丝下游存在尾流。其中，阵列X2的尾流区最大，而阵列X3的尾流区最小。
3) 在区域Ⅲ和Ⅳ处，由于流道曲率和过流截面减小的双重效应，存在高速射流，并在出口圆角下游发生流动分离，各膜丝层径向出流与高速射流交汇时会产生较大的剪切掺混，径向出流被卷吸，剪切层随射流向下游发展，射流迅速卷吸下壁面附近流体，使下壁面附近出现回流。其中阵列X2在出口圆角处流动分离点更靠前，下游产生的涡结构尺度更大。此时，旋涡区存在高剪切力，容易造成溶血和血小板激活。阵列X3由于射流冲击角度小，下游旋涡基本消失。
为获取流场中涡结构的位置和尺度，图 10展示了具有不同倾斜角的4种交错阵列X1、X4、X5和X6的局部流速矢量和涡量分布。可以看出，当孔隙率一定时随着膜丝倾斜角增大，膜丝间的轴向距离减小，膜丝轴截面的长宽比减小，膜丝装填量增加，流动趋于均匀。
 |
| 图 10 4种交错阵列的局部流速矢量及涡量分布(Re=5) |
| 图选项 |
对于倾斜角15°的阵列X4，在区域Ⅰ处主要为径向流动，受流道曲率影响，流动汇聚在内圆角处。在区域Ⅰ和Ⅲ处，由于轴截面上膜丝的迎流尺度大，形成钝体绕流，背流面流动分离区大，背流侧相邻膜丝之间出现涡量高值区，形成旋涡，轴向绕流受阻，该阵列流动损失大，渗透率远小于其他阵列形式。
对于阵列X5，由于膜丝的倾斜角增至30°，轴截面的长短轴比减小，此时流动阻力减小，血液在区域Ⅰ处基本沿径向通过，在内侧圆角流动汇聚得到缓解，膜丝背流侧的旋涡基本消失。在区域Ⅲ处血液绕流阻力减小，流道扩散形成的滞流区范围较大。
对于倾斜角为52°的阵列X1、倾斜角为60°的阵列X6，在区域Ⅰ处为斜向入流，沿轴向流动逐渐加强，过流通路中的流量分配均匀，滞流区小。
在区域Ⅳ处：阵列X4与X1的出口下游出现旋涡，且X4的涡量高于X1，涡尺度更大，流动均匀性差；而阵列X5和X6则无旋涡区，流动较均匀。
ECMO氧合器的能量损失反映了血液流动受到的机械剪切力水平，与血细胞的运动和损伤评估密切相关，可作为调节ECMO循环系统中血液泵功率的依据。本文将血液流动的能量损失与熵产联系起来，考虑壁面熵产的修正，并对比了压降损失与熵产理论损失。表 2展示了Re=5时不同膜丝阵列的氧合区进出口压降Δpb、无量纲渗透率系数kp、壁面平均切应力τave、壁面最大切应力τmax、总熵产Sg[22]、壁面熵产Sgwall[23]、氧合器进出口总压降计算的扬程损失hp、熵产计算的扬程损失hSg。Sgwall、Sg、hp和hSg可分别表示为：
| $S g_{\text {wall }}=\int_\mathit{\Omega} \frac{\boldsymbol{\tau} \cdot \boldsymbol{u}}{T} \mathrm{~d}\mathit \Omega, $ | (13) |
| $\begin{gathered}S g=\frac{\mu}{T}\left[2\left(\frac{\partial u_x}{\partial x}\right)^2+2\left(\frac{\partial u_y}{\partial y}\right)^2+\left(\frac{\partial u_x}{\partial x}+\frac{\partial u_y}{\partial y}\right)^2\right]+ \\S g_{\text {wall }}, \end{gathered}$ | (14) |
| $h_{\mathrm{p}}=\frac{\int_{\text {in }}\left(p_{\mathrm{b}}+0.5 \rho \boldsymbol{u}^2\right) \mathrm{d} q-\int_{\text {out }}\left(p_{\mathrm{b}}+0.5 \rho \boldsymbol{u}^2\right) \mathrm{d} q}{\rho q g}, $ | (15) |
| $h_{S g}=\frac{T S g}{\rho q g} .$ | (16) |
由表 2中结果可知：
表 2 不同膜丝排布形式的流动参数对比(Re=5)
| 排列 | Δpb/mmHg | kp/10-3 | τave/Pa | τmax/Pa | Sg/(10-3 W·K-1) | Sgwall/(10-6 W·K-1) | hp/m | hSg/m |
| 线性阵列 | 553.00 | 1.93 | 21.28 | 127.81 | 79.37 | 10.15 | 7.14 | 4.86 |
| 交错阵列X1 | 70.45 | 16.59 | 3.69 | 12.03 | 7.70 | 0.84 | 0.99 | 0.72 |
| 交错阵列X2 | 104.00 | 9.88 | 4.20 | 15.32 | 8.54 | 1.58 | 2.02 | 1.65 |
| 交错阵列X3 | 35.38 | 22.09 | 2.83 | 12.39 | 5.85 | 0.22 | 0.49 | 0.38 |
| 交错阵列X4 | 218.30 | 0.53 | 18.22 | 297.05 | 34.81 | 1.38 | 2.89 | 2.13 |
| 交错阵列X5 | 115.24 | 5.97 | 10.39 | 42.07 | 15.61 | 4.24 | 1.67 | 1.18 |
| 交错阵列X6 | 76.85 | 11.84 | 8.29 | 58.10 | 7.13 | 5.34 | 0.98 | 0.72 |
表选项
1) 对于不同膜丝阵列，壁面熵产率在总熵产率中占比均较小。渗透率越高，扬程损失越低。熵产理论预测的扬程损失小于实际扬程损失hp。
2) 与具有相同孔隙率的交错阵列相比，线性阵列由于径向间隙较小，所以流动的压差大，表明流动阻力大；最大切应力τmax大，说明间隙处血液损伤风险也大。
3) 对比线性阵列、交错阵列X1和X2的结果可以看出，随着相邻膜丝最小间隙的增大，渗透率增大，平均切应力、总熵产率和壁面熵产率均有所降低，血液损伤风险降低。
4) 对比交错阵列X1、X2和X3的结果可知，X1的轴向、径向间隙最接近，各向均匀性较好，虽然X1的平均切应力高于X3，但其最大切应力最小，说明各向均匀排布有利于避免单方向装填密度过大导致的局部切应力高值。交错阵列X3具有最小的压降和平均切应力，流动性能较优。
5) 对比交错阵列X4、X5和X6以及X1可知，随着倾斜角增加，平均切应力先减小后增大，在倾斜角为52°时达到最小。随着倾斜角增大，膜丝轴截面椭圆长宽比减小，膜丝之间的流道较通畅，流动分离较小，渗透率提高；进一步增大倾斜角，膜丝壁面对流动的约束增加，渗透率降低，总熵产增加，损失主要集中于旋涡和混合流动较强的区域。
3.3 膜丝阵列的气体输运特征根据单膜丝的数值模拟结果可知，膜丝外壁气体分压反映了对流输运和扩散输运的相对强度。图 11展示了Re=5时交错阵列X5在图 2中A、B、C区域的O2和CO2分压云图。总体而言，数值模拟结果揭示了血液以较大流速进入氧合器后，沿径向、轴向随流氧合并排出CO2，使O2分压沿程增大，CO2分压沿程降低。
 |
| 图 11 交错阵列X5在局部区域的气体分压(Re=5) |
| 图选项 |
在A区域，前排膜丝的流动分离使绕后排膜丝流动的驻点发生偏移。在上壁面附近的膜丝尾缘滞流区对流作用较弱，气体以扩散为主，出现O2分压高值区和CO2低值区，该区域膜丝表面浓度边界层迅速增厚，由于气体组分在膜侧顺浓度梯度扩散，气体输运量低。在B区域，截面上的流动呈现周期性。而对于接近出口的区域C，在圆角附近对流较强，壁面附近流速缓慢，气体以扩散输运为主。
选取线性阵列、交错阵列X1和X3，分析相同孔隙率下阵列轴向、径向排布均匀性对气体输运特征的影响。图 12为3种膜丝阵列方式在区域Ⅰ、Ⅱ和Ⅲ的O2分压云图。由对比可知：
 |
| 图 12 3种膜丝阵列在不同区域的O2分压(Re=5) |
| 图选项 |
1) 线性阵列在区域Ⅰ的进口圆角处汇聚流动明显，而膜丝背流面的滞流区很大，导致低O2浓度血液仅与少量膜丝表面接触并快速通过。在区域Ⅱ，膜丝前缘和尾缘较大区域以扩散输运为主；在区域Ⅲ，主流部分O2分压低于300 mmHg，以纯扩散输运模式为主。此时有效对流氧合面积为膜丝上下缘的少量区域，沿径向膜丝前后缘浓度边界层增厚，氧合器底部流速低，以扩散输运为主。
2) 交错阵列X1在区域Ⅰ的入口圆角处流动汇聚不明显，而在区域Ⅱ的气体以对流输运为主，在膜丝滞流区有一定的扩散输运，且下层膜丝尾缘扩散区较大；区域Ⅲ主要表现为对流扩散输运模式。相邻交错膜丝间的流动分离区发生扩散输运，分离点前后的区域发生对流扩散，此时的有效氧合区大于线性阵列，主流部分O2分压约为700 mmHg。
3) 交错阵列X3的流动无明显的周期性扩散和收缩，相邻膜丝之间的流动较顺畅，有效氧合面积最大。在区域Ⅲ处，主要为对流扩散输运模式。血流在出口圆角处无明显流动汇聚，膜丝表面以对流输运为主，扩散输运的流动分离区较X1更小，主流部分O2分压超过700 mmHg。
图 13中展示了Ⅱ区域线性阵列和交错阵列X1中上下层膜丝表面气体Sherwood数沿周向的分布。图中，顺流向取膜丝的上缘点为γ=0°，下缘点为γ=180°，以Sherwood数反映气体的对流传质与扩散传质之比。可以看出，沿膜丝壁面O2和CO2的输运具有相似规律，但CO2对流输运更强，这是因为CO2在血液中的溶解度更高且扩散系数更小。线性阵列仅在γ=100°附近产生较大对流输运，绕流前缘点处几乎完全为扩散输运；交错阵列中上、下层膜丝存在不同的传质特点，上侧膜丝(见图 13b中a1)在前缘点附近、最小径向间隙处气体对流输运均较强，在迎流面流动扩散区、膜丝尾缘附近对流作用较弱。对于下层膜丝(见图 13b中a2)，由于其迎流面与上层膜丝构成间隙流动，对流输运增加，而背流面处流动分离，传质性能较弱。
 |
| 图 13 两种阵列相邻层膜丝表面气体Sherwood数沿周向分布 |
| 图选项 |
图 14为4种不同倾斜角膜丝阵列在Ⅰ、Ⅱ、Ⅲ区域的气体分压云图。由图可知，倾斜角对气体输运的影响存在双重效应。一方面，随着倾斜角减小，血液在进出口附近的Ⅰ、Ⅲ区域径向流动增强，产生大量旋涡与滞留区，气体在这些区域以扩散输运为主，不利于血液氧合；另一方面，在Ⅱ区域处，倾斜角减小使尾缘处的流动分离点向尾部移动，膜丝间隙氧合面积增大、对流输运增强，有利于氧合作用。
 |
| 图 14 不同倾斜角交错排列的局部O2分压云图(Re=5) |
| 图选项 |
表 3列出了不同膜丝阵列及常见商用氧合器的气体输运性能数据，其中Δp为氧合器进出口气体分压差的绝对值，i组分气体体积输运速率mi可表示为
表 3 二维排列及商用氧合器的气体输运性能
| Δpb/mmHg | ΔpO2/mmHg | ΔpCO2/mmHg | mO2/(mL·min-1) | mCO2/(mL·min-1) | ||
| 排布形式 | 线性阵列 | 553.2 | 107 | 4.8 | 134 | 86 |
| 交错阵列X1 | 70.5 | 623 | 33.9 | 222 | 135 | |
| 交错阵列X3 | 35.4 | 655 | 37.1 | 227 | 149 | |
| 交错阵列X4 | 218.3 | 523 | 26.5 | 206 | 105 | |
| 交错阵列X5 | 115.2 | 603 | 30.8 | 219 | 122 | |
| 交错阵列X6 | 76.9 | 621 | 32.2 | 222 | 128 | |
| 氧合器型号 | HMO2000[24] | 34 | — | — | 272 | 234 |
| Sorin Avant[25] | 95 | — | — | 265 | 200 | |
| Medos Hilite7000[26] | 90 | — | — | 300 | 250 | |
| Medtronic Affinity[27] | 45 | — | — | 265 | 225 | |
| Terumo Rx15[28] | 90 | — | — | 260 | 180 |
表选项
| $m_i=q\left[\alpha_i\left(p_{i, \text { out }}-p_{i, \text { in }}\right)+H_{b, i}\left(S_{i, \text { out }}-S_{i, \text { in }}\right)\right] .$ | (17) |
4 结论本文采用基于IB法的自编软件对氧合器内不同膜丝阵列的多相流动进行数值模拟，分析了不同排布方式、倾斜角的膜丝阵列对血液的渗透率、流动特征、能量损失以及气体输运特性的影响。主要结论如下：
1) 基于IB法计算的复杂多孔介质微间隙流动、多区域气体耦合输运结果与商业CFD软件Fluent模拟结果、相关文献中的结果均相符，验证了本文数值模拟方法在预测氧合器内二维稳态多相流动方面具有准确性。
2) 在孔隙率一定的条件下，氧合器的氧合区渗透率同时受膜丝排布方式和倾斜角的影响，而气体输运性能则主要受排布方式的影响。
3) 交错阵列中小倾斜角膜丝的径向渗透率远大于轴向渗透率，易在膜丝区出口形成大尺度涡，流动损失大；而膜丝倾斜角较大时交错阵列的轴向渗透率略高于径向渗透率，流道无周期性收缩扩张，血液流动滞流区较小、出口无旋涡，且气体对流输运增强。
4) 基于膜丝阵列二维轴截面的数值模拟可较好预测交错阵列中膜丝倾斜角较大时的多相流动特征，并准确预测ECMO的压降。但当交错阵列倾斜角较小时，由于血液流动存在较强的周向运动，采用二维模型模拟ECMO中的多相流动则有较大局限性。
参考文献
| [1] | MACLAREN G, BUTT W, BEST D, et al. Extracorporeal membrane oxygenation for refractory septic shock in children: One institution's experience[J]. Pediatric Critical Care Medicine, 2007, 8(5): 447-451. DOI:10.1097/01.PCC.0000282155.25974.8F |
| [2] | TABESH H, AMOABEDINY G, POORKHALIL A, et al. A theoretical model for evaluation of the design of a hollow-fiber membrane oxygenator[J]. Journal of Artificial Organs, 2012, 15(4): 347-356. DOI:10.1007/s10047-012-0655-3 |
| [3] | LOW K W Q, VAN L R, ROLLAND S A, et al. Pore-scale modeling of non-newtonian shear-thinning fluids in blood oxygenator design[J]. Journal of Biomechanical Engineering, 2016, 138(5): 051001. DOI:10.1115/1.4032801 |
| [4] | 罗先武, 叶维祥, 宋雪漪, 等. 支撑"双碳"目标的未来流体机械技术[J]. 清华大学学报(自然科学版), 2022, 62(4): 678-692. LUO X W, YE W X, SONG X Y, et al. Future fluid machinery supporting "double-carbon" targets[J]. Journal of Tsinghua University (Science and Technology), 2022, 62(4): 678-692. (in Chinese) |
| [5] | KIM J J, JANG T S, KWON Y D, et al. Structural study of microporous polypropylene hollow fiber membranes made by the melt-spinning and cold-stretching method[J]. Journal of Membrane Science, 1994, 93(3): 209-215. DOI:10.1016/0376-7388(94)00070-0 |
| [6] | ZHANG J T, NOLAN T D C, ZHANG T, et al. Characterization of membrane blood oxygenation devices using computational fluid dynamics[J]. Journal of Membrane Science, 2006, 288(1-2): 268-279. |
| [7] | HORMES M, BORCHARDT R, MAGER I, et al. A validated CFD model to predict O2 and CO2 transfer within hollow fiber membrane oxygenators[J]. The International Journal of Artificial Organs, 2011, 34(3): 317-325. DOI:10.5301/IJAO.2011.6494 |
| [8] | MAZAHERI A R, AHMADI G. Uniformity of the fluid flow velocities within hollow fiber membranes of blood oxygenation devices[J]. Artificial Organs, 2006, 30(1): 10-15. DOI:10.1111/j.1525-1594.2006.00150.x |
| [9] | 王东璞, 王子奇, 刘爽, 等. 复杂边界和极端条件对单相和多相湍流结构和输运的影响[J]. 清华大学学报(自然科学版), 2022, 62(4): 758-773. WANG D P, WANG Z Q, LIU S, et al. Effects of complicated boundaries and extreme conditions on the flow structure and transport of single-and multi-phase turbulence[J]. Journal of Tsinghua University (Science and Technology), 2022, 62(4): 758-773. DOI:10.16511/j.cnki.qhdxxb.2022.25.014 (in Chinese) |
| [10] | LOW K W Q, VAN LOON R, ROLLAND S A, et al. Formulation of generalized mass transfer correlations for blood oxygenator design[J]. Journal of Biomechanical Engineering, 2017, 139(3): 031007. DOI:10.1115/1.4035535 |
| [11] | MATSUMURA Y, JACKSON T L. Numerical simulation of fluid flow through random packs of cylinders using immersed boundary method[J]. Physics of Fluids, 2014, 26(4): 043602. DOI:10.1063/1.4870246 |
| [12] | SVITEK R G, FEDERSPIEL W J. A mathematical model to predict CO2 removal in hollow fiber membrane oxygenators[J]. Annals of Biomedical Engineering, 2008, 36(6): 992-1003. DOI:10.1007/s10439-008-9482-3 |
| [13] | LYSAGHT M, WEBSTER T J. Biomaterials for artificial organs[M]. Cambrige: Woodhead Publishing, 2011. |
| [14] | DASH R K, BASSINGTHWAIGHTE J B. Blood HbO2 and HbCO2 dissociation curves at varied O2, CO2, pH, 2, 3-DPG and temperature levels[J]. Annals of Biomedical Engineering, 2004, 32(12): 1676-1693. DOI:10.1007/s10439-004-7821-6 |
| [15] | ZIERENBERG J R, FUJIOKA H, HIRSCHL R B, et al. Oxygen and carbon dioxide transport in time-dependent blood flow past fiber rectangular arrays[J]. Physics of Fluids, 2009, 21(3): 033101. DOI:10.1063/1.3056413 |
| [16] | TASKIN M E, FRASER K H, ZHANG T, et al. Micro-scale modeling of flow and oxygen transfer in hollow-fiber membrane bundle[J]. Journal of Membrane Science, 2010, 362(1-2): 172-183. DOI:10.1016/j.memsci.2010.06.034 |
| [17] | HARASEK M, LUKITSCH B, ECKER P, et al. Fully resolved computational (CFD) and experimental analysis of pressure drop and blood gas transport in a hollow fibre membrane oxygenator module[J]. Chemical Engineering Transactions, 2019, 76: 193-198. |
| [18] | 许梦菲, 梁亚静, 臧慧, 等. 聚-4-甲基-1-戊烯中空纤维膜式人工肺组件的氧气传质性能的研究[J]. 南京大学学报(自然科学), 2021, 57(4): 641-647. XU M F, LIANG Y J, ZANG H, et al. Study on oxygen mass transfer performance of poly-4-methyl-1-pentene hollow fiber membrane artificial lung module[J]. Journal of Nanjing University (Science and Technology), 2021, 57(4): 641-647. (in Chinese) |
| [19] | HEWITT T J, HATTLER B G, FEDERSPIEL W J. A mathematical model of gas exchange in an intravenous membrane oxygenator[J]. Annals of Biomedical Engineering, 1998, 26(1): 166-178. DOI:10.1114/1.53 |
| [20] | OLDROYD J G. Non-newtonian effects in steady motion of some idealized elastico-viscous liquids[J]. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 1958, 245(1241): 278-297. |
| [21] | KARIMI S, DABAGH M, VASAVA P, et al. Effect of rheological models on the hemodynamics within human aorta: CFD study on CT image-based geometry[J]. Journal of Non-Newtonian Fluid Mechanics, 2014, 207: 42-52. DOI:10.1016/j.jnnfm.2014.03.007 |
| [22] | BEJAN A. Entropy generation minimization-the method of thermodynamic optimization of finite-size systems and finite-time processes[M]. New York: CRC Press, 1995. |
| [23] | 张翔, 王洋, 徐小敏, 等. 低比转数离心泵叶轮内能量转换特性[J]. 农业机械学报, 2011, 42(7): 75-81. ZHANG X, WANG Y, XU X M, et al. Energy conversion characteristic within impeller of low specific speed centrifugal pump[J]. Transactions of the Chinese Society for Agricultural Machinery, 2011, 42(7): 75-81. (in Chinese) |
| [24] | GETINGE INC. Extracorporeal life support[EB/OL]. [2022-05-08]. https://www.Getinge.com/int/product-catalog/pls-system/. |
| [25] | SORIN INC. EOS ECMO intended for long-duration procedures: hollow fiber oxygenator. [EB/OL][2022-10-03]. https://assets.livanova.com/media/com/documents/therapeutic%20areas/cardiopulmonary/brochures/ous_only_livanova_livan-ova_eos_ecmo_09295-85b.pdf?ext=.pdf. |
| [26] | XENIOS INS. Hilite product family[EB/OL]. [2022-05-08]. https://www.xenios-ag.com/medos/products/hilite-oxygenators/. |
| [27] | MEDTRONIC INC. Medtronic extracorporeal support and circulation technologies[EB/OL]. [2022-05-08]. https://global.medtronic.com/xg-en/healthcare-professionals/products/cardiovascular/extracorporeal-life-support/nautilus-smart-ecmo-module.html. |
| [28] | TERUMO INC. Use of terumo cardiovascular's CAPIOX. FX25oxygenator in extracorporeal membrane oxygenation therapy[EB/OL]. [2022-05-08]. https://www.terumocv.com/doc/CAPIOX-FX25-Oxygenator-in-ECMO-Therapy_April2020_FINAL.pdf. |
