体外膜肺氧合技术(extracorporeal membrane oxygenation, ECMO)是一种在体外通过对静脉血中氧气(O₂)和二氧化碳(CO₂)进行更新，从而实现短期替代人体肺功能的呼吸支持技术，在新型冠状病毒感染等导致的急性呼吸窘迫综合征治疗，肺移植、心脏瓣膜等心血管外科手术中发挥着重要作用^[1]。ECMO氧合器需要满足血液灌注量小、压降低、气体交换效率高、生物相容性好等指标要求。因此，如何在降低血液损伤的前提下提升气体交换效率、延长呼吸支持周期是ECMO氧合器设计和优化的主要目标。目前通用的氧合器采用中空纤维膜结构，膜丝涂层改性、组件的血流动力学优化是提高氧合器生物相容性、优化血流模式、降低传质阻力的重要手段。Tabesh等^[2]提出采用膜丝装填密度、O₂输运量和压降损失率等指标评价氧合器性能，并通过实验对比了6款商用ECMO设备性能。

与临床实验相比，计算流体动力学(computational fluid dynamics, CFD) 可快速准确获得不同膜丝结构的流场参数，是研究ECMO内流动的有力工具^[3]，精细化的流动模拟与流场诊断技术可揭示ECMO微观流动和宏观外特性之间的关系，为ECMO样机设计优化提供科学依据，提高改型优化的准确性^[4]。由于商用ECMO设备膜丝外径为300~380 μm，膜表面积达1.2~2.5 m²，且排列紧密^[5]，传统的计算方法难以满足贴体边界层网格数量与质量的要求。目前ECMO氧合器流动模拟主要采用多孔介质模型^[6]和微尺度周期阵列模型^[7]等简化方法，但各向同性多孔介质模型中血液均匀流动的假设与膜丝周围血流运行工况不符^[8]，需要进行修正才能将计算误差控制在允许范围内^[3]。

为表示不同膜丝排列下多孔介质的流动阻力，需要对表征单元体积(representative element volume)进行微尺度数值模拟。王东璞等^[9]通过对有限尺寸固体颗粒阵列孔隙尺度直接进行数值模拟，研究了多孔介质结构与流动的相互作用对传热特性、流体质点输运的影响。Low等^[10]通过微尺度二维数值模拟，获得了线性和交错周期阵列下血液沿轴向、横向流动时渗透率和气体扩散系数与孔隙率的关系，为预测实际情况下ECMO的压降和氧合效率提供了参考。文[11]利用沉浸边界(immersed boundary，IB) 法对二维随机和周期排列的圆形、方形和椭圆阵列绕流进行稳态模拟，得出了渗透率随孔隙率、Reynolds数、截面形状和长宽比的变化关系。

ECMO氧合器气体输运包括静脉血的O₂富集和CO₂清除，膜丝空腔内的气体在浓度差作用下在膜侧区域发生Knudsen扩散；在膜外壁血液区以溶解组分、血红蛋白结合组分、反应物组分^[12]完成对流扩散输运^[13]，其中血红蛋白与O₂和CO₂的亲和曲线具有Haldane-Bohr耦合效应^[14]。该过程涉及多区域耦合及大量生化反应。Zierenberg等^[15]通过血液侧O₂和CO₂耦合输运的数值模拟获得了不同二维阵列的瞬态气体输运特性，发现血液损伤风险与气体输运效率存在互相制约的关系。Taskin等^[16]提出了一种多区域O₂输运模型，同时考虑血液侧、膜侧传质阻力，能减小传质效率的过度预测。Harasek等^[17]构建了耦合血液侧、纤维膜侧和纤维膜内腔的全区域CO₂输运模型，CO₂浓度的数值计算结果与实验相符。

综上所述，目前微尺度模型研究的膜丝阵列多为排列简单的二维圆柱阵列，其渗透率和气体输运的结果能否应用于实际复杂模型还需进一步验证；O₂和CO₂多区域耦合的输运研究较少，尚不能揭示血液流动与气体输运的联系。

考虑到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)

其中: $\nabla$为Nabla算子, $\boldsymbol{u}$为流速, $p_{\mathrm{b}}$为血液压力, $\rho$为流体密度, $\mu$为流体动力黏度。

ECMO氧合器内部流动满足Darcy定律，通常采用无量纲渗透率系数k_p来描述多孔介质对流动的阻碍作用，具体可表示为

$ 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)

其中: $U$为表观流速, $f$为沿流动方向长度, $A_{0}$为湿周, $e_{0} 、b_{0}$分别为氧合器膜丝椭圆截面的长、短轴。

计算中，通过设定流体域边界质量不平衡残差、流体网格速度和压力的累积残差阈值来控制迭代次数，保证计算结果收敛。

为求解上述基本控制方程，本文利用FOAM-extend 4.0平台自编了基于清晰界面IB方法的流动模拟软件。其中，边界通过三角形贴体网格离散，背景网格采用非结构化同位Cartesian网格。以图 1中二维单膜丝计算域为例，计算时背景网格被划分为流体网格、IB网格和固体网格，图中分别用圆形、三角形和菱形表示，其中流体网格和IB网格根据区域可进一步分成膜侧和血液侧网格，分别用蓝色和红色进行标记。固体网格位于膜丝内壁内侧，膜侧网格位于膜丝内外壁面之间。膜侧网格中与固体网格和血液侧网格存在交界面的为膜侧IB网格，其余为膜侧流体网格。血液侧网格位于膜丝外壁的外部，其中与膜侧IB网格存在交界面的为血液侧IB网格，其余网格为血液侧流体网格。图中，IB网格的边界面法线与沉浸边界的交点为IB点，记为ibp，血液侧ibp和膜侧ibp通常不重合，分别用红色和蓝色星形表示。为了计算边界物理量的梯度，每个膜侧ibp在面法向上设置1个采样点，对应图中红色矩形处。仅对流体网格进行离散化计算，固体网格的速度和压力可以进行预设。IB网格的速度和压力通过边界点和影响域内流体点经过加权的最小二乘法进行插值预估，以满足无滑移、无梯度等边界条件。流动计算时，膜侧气体扩散速率可忽略不计^[16]，将膜侧流体网格速度预设为0。通过修正IB与流体网格交界面的流量通量实现流场迭代计算。

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)

其中: $i$为$\mathrm{O}_{2}$或$\mathrm{CO}_{2} ; \alpha_{i}$为气体溶解度, $\mathrm{mmHg}^{-1} ; p_{i}$为血液中气体组分的分压, $\mathrm{mmHg}$; $D_{\text {eff, } i}$为血液中气体有效扩散系数, $37^{\circ} \mathrm{C}, \mathrm{m}^{2} / \mathrm{s}$; $D_{\mathrm{rbc}}$为血浆中红细胞扩散系数, $\mathrm{m}^{2} / \mathrm{s} ; H_{\mathrm{b}, i}$表示红细胞压积为0.45时血红蛋白的携气能力, $\mathrm{mL} /$ $\mathrm{mL} ; S_{i}$为血液中气体饱和度，$\% ; t$为时间，$\mathrm{s}$; $D_{\mathrm{p}, i}$为膜侧的组分扩散系数, $\mathrm{m}^{2} / \mathrm{s}$。

本文使用分离式求解方法实现膜丝与血液交界面耦合。计算血液侧时，交界面为固定值边界，使用上一次迭代后膜丝外壁的气体分压值，此时第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)

其中: $a_{1} \sim a_{5}$为最小二乘法插值系数, $X_{j} 、Y_{j}$分别表示第$j$个IB网格中心到IB的距离对应$x 、y$方向的分量, $\omega_{j}$为加权系数, $\phi_{\mathrm{ibp}}$为边界上物理量的值。

计算膜丝侧时，血液侧气体分压数值不变，交界面为耦合通量边界，保证两个区域与交界面的气体扩散通量相等。耦合通量边界条件中ϕ_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)

其中: $c_{0} \sim c_{5}$为最小二乘法插值系数, $\boldsymbol{n}$为边界单位外法向向量, $\phi_{j, \mathrm{f}}$为第$j$个IB网格对应采样点物理量的值, $d$为采样点到沉浸边界的距离。当某侧区域迭代至收玫标准后再计算另一侧区域。

血红蛋白的携O₂、携CO₂能力存在耦合效应，2种气体仅在血液侧耦合输运，在膜侧以相互独立的气体组分存在。暂不考虑pH值、温度、化学物质以及各组分化学反应对输运的影响。O₂和CO₂的饱和度根据文[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)

其中: $n_{i}$为经验系数; $p_{k}$为$k$气体组分分压, $\mathrm{mmHg} ; k$为$\mathrm{O}_{2}$或$\mathrm{CO}_{2}$成分; $p_{50, }$为饱和度为$50 \%$时对应的气体分压, $\mathrm{mmHg} ; \beta_{1, i}$为经验系数, $\mathrm{mmHg} ; \beta_{2, i} 、\beta_{3, i}$均为经验系数, $\mathrm{mmHg}^{-1}$。

计算血液侧两种气体组分耦合输运时，采用分离式求解方法，即先计算某组分的对流扩散方程，收敛到一定残差后再对另一组分进行更新。采用局部壁面的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)

其中：$r$为膜丝外径, $C_{\text {wall, } i} 、C_{\mathrm{in}, i}$分别为膜丝外壁和人口处的气体组分浓度。

2 计算模型 2.1 计算域

本文研究对象为某商用ECMO氧合器模型，其氧合区为同心圆筒状结构，膜丝轴切面为椭圆环，不同轴面上的孔隙率保持不变，但排列不同。为比较不同轴截面及缠绕角度下的流动，假设氧合区的周向流动可以忽略，只选取氧合区的一个轴截面进行二维定常流动模拟。为避免在ECMO的进出口出现回流，选定计算域时在膜丝阵列的进、出口分别延长了30、30.3 mm。图 2所示为计算域，其中：沿x方向为径向，沿y方向为轴向，沿z方向为垂直轴截面方向，血液从最左侧断面流入，从最右侧断面流出。为便于分析，选取Ⅰ—Ⅳ共4个代表性区域，其中：Ⅰ位于膜丝区入口近轴线的内侧圆角处，Ⅱ位于膜丝区中间段的局部区域，Ⅲ位于膜丝区出口远离轴线的外侧圆角处，Ⅳ位于膜丝区出口延长段区域。此外，还选取了包含膜丝区整个径向长度的区域A、B和C，其中：A、B、C分别位于膜丝区上壁面、中段和下壁面处。

2.2 ECMO膜丝阵列

图 3所示为ECMO内2根交错排布膜丝的空间关系，其中一根膜丝的中心线与氧合器轴线平行，膜丝以蓝色线表示；另一根膜丝的中心线与氧合器轴线不平行，膜丝以红色线表示，而红色椭圆环表示该膜丝在轴截面上的形状。两膜丝中心线位于距氧合器轴线为R的圆柱平面上，膜丝中心线与氧合器轴线之间的倾斜角为θ(以逆时针方向为正)。

图 4所示为ECMO膜丝阵列的7种排布方式，图中未表示膜丝内侧的结构。图 4a为线性阵列，所有膜丝的倾斜角均相等，θ₀=52°。膜丝阵列的列间距l_x=0.36 mm，行间距l_y=0.788 1 mm。图 4b为交错阵列X1。将线性阵列中偶数列的膜丝倾斜角变为θ₀，并改变两相邻交错膜丝的轴向间距。图中，交错膜丝中心距离a=0.281 mm。图 4c为交错阵列X2。在交错阵列X1的基础上，改变两相邻交错膜丝的轴向间距。图中，交错膜丝中心距b=0.188 mm。图 4d为交错阵列X3。在X1的基础上，改变两相邻交错膜丝的轴向间距和径向间距。图中，交错膜丝中心距c=0.25 mm。图 4e为交错阵列X4，在X1的基础上，改变膜丝倾斜角，θ₁=15°。此时，为保证膜丝阵列孔隙率不变，将行间距改为l₁=2.028 mm。图 4f为交错阵列X5，在X1的基础上，改变膜丝倾斜角，θ₂=30°。此时，为保证膜丝阵列孔隙率不变，将行间距改为l₂=1.050 mm。图 4g为交错阵列X6，在X1的基础上，改变膜丝倾斜角，θ₃=60°。此时，为保证膜丝阵列孔隙率不变，将行间距改为l₃=0.606 mm。

图 4中，阵列X1、X2和X3为孔隙率一定时不同膜丝间隙下的交错阵列；X1、X4、X5和X6为孔隙率一定时不同膜丝倾斜角下的交错阵列。

膜丝外径r=0.3 mm，内径r₀=0.24 mm，各阵列孔隙率ε=0.683。在膜丝内的O₂、CO₂经由膜丝壁面与血液进行质量交换。

2.3 计算条件与网格

血液的流变呈现黏弹性和剪切稀化特征。Taskin等^[16]通过计算氧合器二维阵列O₂输运结果对比了改进Casson模型^[19]和考虑血液黏弹性的Oldroyd模型^[20]的预测效果，发现Re=5时2种模型的壁面O₂分压的预测相对误差在10%以内，血液黏弹性对气体输运的影响可以忽略不计。Karimi等^[21]通过数值计算对比了9种血液流变模型对主动脉内血流动力学的预测结果，其中Carreau模型和改进Casson模型在较低和较高的切应力水平下与血液的实际流变特点相符。本文采用交错阵列X5在Re=5的条件下比较Carreau模型和改进Casson模型对流动和输运的预测效果，2种模型预测的渗透率相对误差为0.135%。图 5展示了2种流变模型O₂分压沿轴线分布的预测结果，其中L表示轴向坐标，膜丝区入口位于L=-0.031 m处。2种模型计算的O₂分压沿轴线的累积相对误差为7.7%，本文数值模拟采用改进Casson模型，黏度μ与剪切率φ的关系可表示为

$ \mu=\left(\left(\frac{\tau_{0}}{\varphi+1}\right)^{1 / 2}+\mu_{c}^{1 / 2}\right)^{2} $

(12)

其中: 对于人体血液而言, Casson屈服应力$\tau_{0}=$ $0.005 \mathrm{~Pa}$; Casson黏度$\mu_{\mathrm{c}}=0.0035 \mathrm{~Pa} \cdot \mathrm{s}$, 密度$\rho=1056 \mathrm{~kg} / \mathrm{m}^{3}$。

血液流动计算时, 人口边界条件为恒定速度, 满足$R e=5$; 出口给定恒定总压; 膜丝外壁、氧合器壁面均满足无滑移条件。气体输运计算时, 忽略膜丝内腔的气体输运阻力, 血液人口$\mathrm{O}_{2}$分压为$56 \mathrm{mmHg}(1 \mathrm{mmHg}=133.3224 \mathrm{~Pa}), \mathrm{CO}_{2}$分压为$45 \mathrm{mmHg}$, 膜侧内壁$\mathrm{O}_{2}$分压为$760 \mathrm{mmHg}, \mathrm{CO}_{2}$分压为0, 其余经验系数依据文[10]取值。

计算域的背景网格为正交网格，膜丝采用三角面网格表示，依据不同缠绕方式所得阵列的膜丝数量为4 509~8 172根。为确保网格分辨率，对膜丝浸没区域2个坐标方向进行等比例局部加密，并进行网格无关性验证。表 1显示了在不同网格下，Re=5、θ=45°、ε=0.735时无量纲渗透率系数k_p、壁面平均切应力τ_ave的计算结果。当网格节点数大于1.1×10⁷时，采用IB方法预测的k_p、τ_ave均变化较小，因此选取1.1×10⁷网格进行数值模拟。

3 结果分析 3.1 模型验证

为了验证IB方法在计算复杂多孔介质微间隙流动的准确性，本文对自编软件和商业CFD软件Fluent的模拟结果进行了对比。需要说明的是，Fluent模拟的计算域、边界条件与IB方法相同，网格数为0.689×10⁷，计算得到k_p=1.13×10^-2，2种方法预测k_p的相对误差为3.5%，说明本文数值模拟方法具有较好的适应性。

图 6为Re分别为5、20、40时，采用IB方法、商业软件Fluent以及文[16]计算的O₂在血液和单膜丝区域内输运的结果，图中数据表示膜丝外壁O₂分压随周向角γ分布，膜丝前缘点为γ=0°，后缘点为γ=180°。通过比较可知，当Re=5时，O₂分压处于710~750 mmHg范围内，且在0°~180°范围内单调增加，3组数据具有较好的一致性。当Re=20时，3组数据均显示在140°附近出现最大的O₂分压。此时发生分离流动，并使角度在进一步增大时的对流传质减弱，O₂分压下降。当Re=40时，IB的计算结果与Fluent的模拟结果均表明绕膜丝壁面的O₂流动在60°附近发生分离，然后在130°附近再次发生分离，在第2次分离点处O₂分压出现极大值，之后随着角度增大O₂分压下降。图中结果验证了本文多区域O₂输运计算方法的准确性。

为进一步验证多组分气体输运效果，模拟了O₂和CO₂在血液侧的耦合输运。图 7为Re=10、ε=0.3时二维交错膜丝阵列中2种气体组分沿膜丝壁面的Sherwood数随周向角分布情况，计算结果与文[10]数据吻合，且均表明在0°~180°范围内出现2个明显的对流传质峰值。由此可知，本文所采用的数值方法在模拟膜丝阵列多相流动方面具有良好的适应性。

3.2 不同膜丝阵列的血流特性

图 8展示了采用不同阵列时沿氧合器轴向的血液压力分布。图中，曲线斜率绝对值越大，所对应的渗透率越低。结果表明：线性阵列的进出口压降高达7×10⁵ Pa，远高于其余阵列；对比交错阵列X1和X2可知，随着轴向间隙增加，渗透率增大；对比交错阵列X1和X3可知，相邻膜丝间设置径向交错间隙较轴向交错间隙更利于血液通过。对比膜丝倾斜角不同时的交错阵列X1、X4、X5和X6可知，膜丝倾斜角为15°时(即X4)压降较大，而其他阵列的压降均较小。因此，当膜丝阵列的截面孔隙率和列间距保持一定时，随着倾斜角增大，膜丝阵列行间距减小，渗透率呈现先增大后减小的趋势。

图 9展示了线性阵列与交错阵列X1、X2和X3的流速分布及速度矢量图。对于每种排列方式，从左至右分别对应图 2Ⅰ—Ⅳ局部区域的流场。由图中结果可知，在膜丝区进口区域Ⅰ、出口区域Ⅲ和Ⅳ的流动非均匀性较强，而中间区域Ⅱ则可视为顺轴向的周期性绕流。

对于线性阵列，在入口区域Ⅰ处的轴向流动阻力大，血液接触膜丝后主要体现为径向流动，并伴随绕椭圆阵列的流动；在区域Ⅱ处的间隙流速高，膜丝前缘与尾缘附近存在明显分离流动；而在出口区域Ⅲ和Ⅳ处，流线在圆角外侧汇聚，在圆角区下方存在较大的低速区。

对于交错阵列X1、X2和X3，流动特征如下：

1) 在区域Ⅰ处，阵列X2的流动冲击最大，阵列X3的流动最均匀，但阵列X3的轴向渗透率显著高于径向。

2) 在区域Ⅱ处，血液顺轴向绕流膜丝时，在膜丝下游存在尾流。其中，阵列X2的尾流区最大，而阵列X3的尾流区最小。

3) 在区域Ⅲ和Ⅳ处，由于流道曲率和过流截面减小的双重效应，存在高速射流，并在出口圆角下游发生流动分离，各膜丝层径向出流与高速射流交汇时会产生较大的剪切掺混，径向出流被卷吸，剪切层随射流向下游发展，射流迅速卷吸下壁面附近流体，使下壁面附近出现回流。其中阵列X2在出口圆角处流动分离点更靠前，下游产生的涡结构尺度更大。此时，旋涡区存在高剪切力，容易造成溶血和血小板激活。阵列X3由于射流冲击角度小，下游旋涡基本消失。

为获取流场中涡结构的位置和尺度，图 10展示了具有不同倾斜角的4种交错阵列X1、X4、X5和X6的局部流速矢量和涡量分布。可以看出，当孔隙率一定时随着膜丝倾斜角增大，膜丝间的轴向距离减小，膜丝轴截面的长宽比减小，膜丝装填量增加，流动趋于均匀。

对于倾斜角15°的阵列X4，在区域Ⅰ处主要为径向流动，受流道曲率影响，流动汇聚在内圆角处。在区域Ⅰ和Ⅲ处，由于轴截面上膜丝的迎流尺度大，形成钝体绕流，背流面流动分离区大，背流侧相邻膜丝之间出现涡量高值区，形成旋涡，轴向绕流受阻，该阵列流动损失大，渗透率远小于其他阵列形式。

对于阵列X5，由于膜丝的倾斜角增至30°，轴截面的长短轴比减小，此时流动阻力减小，血液在区域Ⅰ处基本沿径向通过，在内侧圆角流动汇聚得到缓解，膜丝背流侧的旋涡基本消失。在区域Ⅲ处血液绕流阻力减小，流道扩散形成的滞流区范围较大。

对于倾斜角为52°的阵列X1、倾斜角为60°的阵列X6，在区域Ⅰ处为斜向入流，沿轴向流动逐渐加强，过流通路中的流量分配均匀，滞流区小。

在区域Ⅳ处：阵列X4与X1的出口下游出现旋涡，且X4的涡量高于X1，涡尺度更大，流动均匀性差；而阵列X5和X6则无旋涡区，流动较均匀。

ECMO氧合器的能量损失反映了血液流动受到的机械剪切力水平，与血细胞的运动和损伤评估密切相关，可作为调节ECMO循环系统中血液泵功率的依据。本文将血液流动的能量损失与熵产联系起来，考虑壁面熵产的修正，并对比了压降损失与熵产理论损失。表 2展示了Re=5时不同膜丝阵列的氧合区进出口压降Δp_b、无量纲渗透率系数k_p、壁面平均切应力τ_ave、壁面最大切应力τ_max、总熵产Sg^[22]、壁面熵产Sg_wall^[23]、氧合器进出口总压降计算的扬程损失h_p、熵产计算的扬程损失h_Sg。Sg_wall、Sg、h_p和h_Sg可分别表示为：

$ 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)

其中: $\boldsymbol{u}_{x} 、\boldsymbol{u}_{y}$分别为$\boldsymbol{u}$在$x 、y$方向的速度, $\boldsymbol{\tau}$为壁面切应力, $T$为温度, $\mathit \Omega$为壁面表面积, in、out分别表示计算域的进口与出口, $g$为重力加速度, $q$为单位时间内血液的体积流量。

由表 2中结果可知：

1) 对于不同膜丝阵列，壁面熵产率在总熵产率中占比均较小。渗透率越高，扬程损失越低。熵产理论预测的扬程损失小于实际扬程损失h_p。

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区域的O₂和CO₂分压云图。总体而言，数值模拟结果揭示了血液以较大流速进入氧合器后，沿径向、轴向随流氧合并排出CO₂，使O₂分压沿程增大，CO₂分压沿程降低。

在A区域，前排膜丝的流动分离使绕后排膜丝流动的驻点发生偏移。在上壁面附近的膜丝尾缘滞流区对流作用较弱，气体以扩散为主，出现O₂分压高值区和CO₂低值区，该区域膜丝表面浓度边界层迅速增厚，由于气体组分在膜侧顺浓度梯度扩散，气体输运量低。在B区域，截面上的流动呈现周期性。而对于接近出口的区域C，在圆角附近对流较强，壁面附近流速缓慢，气体以扩散输运为主。

选取线性阵列、交错阵列X1和X3，分析相同孔隙率下阵列轴向、径向排布均匀性对气体输运特征的影响。图 12为3种膜丝阵列方式在区域Ⅰ、Ⅱ和Ⅲ的O₂分压云图。由对比可知：

1) 线性阵列在区域Ⅰ的进口圆角处汇聚流动明显，而膜丝背流面的滞流区很大，导致低O₂浓度血液仅与少量膜丝表面接触并快速通过。在区域Ⅱ，膜丝前缘和尾缘较大区域以扩散输运为主；在区域Ⅲ，主流部分O₂分压低于300 mmHg，以纯扩散输运模式为主。此时有效对流氧合面积为膜丝上下缘的少量区域，沿径向膜丝前后缘浓度边界层增厚，氧合器底部流速低，以扩散输运为主。

2) 交错阵列X1在区域Ⅰ的入口圆角处流动汇聚不明显，而在区域Ⅱ的气体以对流输运为主，在膜丝滞流区有一定的扩散输运，且下层膜丝尾缘扩散区较大；区域Ⅲ主要表现为对流扩散输运模式。相邻交错膜丝间的流动分离区发生扩散输运，分离点前后的区域发生对流扩散，此时的有效氧合区大于线性阵列，主流部分O₂分压约为700 mmHg。

3) 交错阵列X3的流动无明显的周期性扩散和收缩，相邻膜丝之间的流动较顺畅，有效氧合面积最大。在区域Ⅲ处，主要为对流扩散输运模式。血流在出口圆角处无明显流动汇聚，膜丝表面以对流输运为主，扩散输运的流动分离区较X1更小，主流部分O₂分压超过700 mmHg。

图 13中展示了Ⅱ区域线性阵列和交错阵列X1中上下层膜丝表面气体Sherwood数沿周向的分布。图中，顺流向取膜丝的上缘点为γ=0°，下缘点为γ=180°，以Sherwood数反映气体的对流传质与扩散传质之比。可以看出，沿膜丝壁面O₂和CO₂的输运具有相似规律，但CO₂对流输运更强，这是因为CO₂在血液中的溶解度更高且扩散系数更小。线性阵列仅在γ=100°附近产生较大对流输运，绕流前缘点处几乎完全为扩散输运；交错阵列中上、下层膜丝存在不同的传质特点，上侧膜丝(见图 13b中a1)在前缘点附近、最小径向间隙处气体对流输运均较强，在迎流面流动扩散区、膜丝尾缘附近对流作用较弱。对于下层膜丝(见图 13b中a2)，由于其迎流面与上层膜丝构成间隙流动，对流输运增加，而背流面处流动分离，传质性能较弱。

图 14为4种不同倾斜角膜丝阵列在Ⅰ、Ⅱ、Ⅲ区域的气体分压云图。由图可知，倾斜角对气体输运的影响存在双重效应。一方面，随着倾斜角减小，血液在进出口附近的Ⅰ、Ⅲ区域径向流动增强，产生大量旋涡与滞留区，气体在这些区域以扩散输运为主，不利于血液氧合；另一方面，在Ⅱ区域处，倾斜角减小使尾缘处的流动分离点向尾部移动，膜丝间隙氧合面积增大、对流输运增强，有利于氧合作用。

表 3列出了不同膜丝阵列及常见商用氧合器的气体输运性能数据，其中Δp为氧合器进出口气体分压差的绝对值，i组分气体体积输运速率m_i可表示为

$ 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)

从表 3可以看出，相同排列方式和孔隙率条件下，随着倾斜角增大，气体输运性能略有提高。径向间隙不均匀的交错阵列X3的气体输运性能最优，两种气体输运速率均高于X1，约为线性排列的1.7倍。因此，交错阵列的气体输运性能远优于线性阵列，且增大倾斜角也有助于提高气体输运量。与商用氧合器相比，膜丝阵列轴截面的压降数据处于商用氧合器运行范围内，但气体输运性能则低于商用氧合器。因为二维膜丝阵列未考虑垂直于轴面的周向流动，导致预测的血液氧合和CO₂清除性能均低于商用氧合器。

4 结论

本文采用基于IB法的自编软件对氧合器内不同膜丝阵列的多相流动进行数值模拟，分析了不同排布方式、倾斜角的膜丝阵列对血液的渗透率、流动特征、能量损失以及气体输运特性的影响。主要结论如下：

1) 基于IB法计算的复杂多孔介质微间隙流动、多区域气体耦合输运结果与商业CFD软件Fluent模拟结果、相关文献中的结果均相符，验证了本文数值模拟方法在预测氧合器内二维稳态多相流动方面具有准确性。

2) 在孔隙率一定的条件下，氧合器的氧合区渗透率同时受膜丝排布方式和倾斜角的影响，而气体输运性能则主要受排布方式的影响。

3) 交错阵列中小倾斜角膜丝的径向渗透率远大于轴向渗透率，易在膜丝区出口形成大尺度涡，流动损失大；而膜丝倾斜角较大时交错阵列的轴向渗透率略高于径向渗透率，流道无周期性收缩扩张，血液流动滞流区较小、出口无旋涡，且气体对流输运增强。

4) 基于膜丝阵列二维轴截面的数值模拟可较好预测交错阵列中膜丝倾斜角较大时的多相流动特征，并准确预测ECMO的压降。但当交错阵列倾斜角较小时，由于血液流动存在较强的周向运动，采用二维模型模拟ECMO中的多相流动则有较大局限性。