了解人工肝中空纤维膜组件传质行为对其高效解毒具有重要指导意义。建立了中空纤维膜组件内流动传质的三维数值模拟模型，研究了小分子胆红素与大分子牛血清白蛋白（BSA）在组件内的传质行为，讨论了管程流量、壳程流量、中空纤维管长度对胆红素与BSA传质的影响。研究结果表明小分子胆红素的传质方式是对流传质和扩散传质二者兼具，而大分子BSA的传质方式以对流传质为主；胆红素和BSA的清除速率随着管程流量增大而提升，也随着中空纤维管长度增加而提升；随着壳程流量增大，胆红素的清除速率先快速提升再缓慢上升达到渐近值，而BSA的清除速率随壳程流量增大而逐渐减小。研究成果可以为人工肝中空纤维膜组件的结构优化设计和临床操作参数提供依据。
引用本文： 王子恒, 许少锋, 余一帆, 陆俊杰, 张学昌. 人工肝中空纤维膜组件内胆红素与牛血清白蛋白传质行为研究. 生物医学工程学杂志, 2024, 41(4): 742-750. doi: 10.7507/1001-5515.202311011 复制
0 引言
我国病毒性肝炎发病率很高，仅乙型肝炎一项，病毒感染者就有近1亿人，其中约7 000万为慢性乙型肝炎患者，约占全世界慢性乙型肝炎患者人数的1/5[1-3]。乙型肝炎病毒（hepatitis B virus，HBV）感染也是一个严重的世界性公共卫生问题，据世界卫生组织报道，全球70多亿人口中，超过20亿有乙型肝炎病毒感染证据，每年有110万人死于HBV和丙型肝炎病毒感染[4]。我国每年约有38.3万人死于肝癌，占全球肝癌死亡病例数的51%，全球每30 s就有1人死于病毒性肝炎相关疾病[5-6]。人工肝通过一个体外的机械、物理、生物装置，驱动患者血液流经体外人工肝器械，担负起暂时辅助严重病变肝脏的功能，清除各种有害物质，代偿肝脏的代谢功能，直至人体自身肝脏功能恢复或进行肝脏移植[7]。因此，进行人工肝技术的研究和开发具有重要的社会价值和科学意义。
目前关于人工肝中空纤维膜组件内的物质转运效率还不十分清楚，科研工作者对中空纤维膜组件的传质行为开展了大量的实验和数值模拟研究。Ilaria等[8]建立了对流强化型的中空纤维膜传质数学模型，得到了葡萄糖和氧气在中空纤维组件内的传输行为。Curcio等[9]结合实验和理论研究了聚醚砜膜和聚砜膜两种中空纤维膜组件的对流传质和扩散传质机制。Raff等[10]分析了人工肝中白蛋白修饰膜的传质特点，发现对流传质和扩散传质兼而有之；之后他们课题组的Donato等[11]又基于有限元方法建立了中空纤维透析膜的对流传质模型，结果显示增大中空纤维膜的压力系数和填充密度可以增强纤维膜的对流传质效果。Nedredal等[12]利用实验的方法研究了胆红素和免疫球蛋白在人工肝中空纤维组件内的传输行为，发现只有在强对流作用下才具有较强的毒素清除效果。Sorrell等[13]研究了醋氨酚肝毒素和氧气在中空纤维组件的内腔、纤维膜膜孔、外腔的传输行为，但该模型假设中空纤维管内物质传输为对流和扩散传质，膜厚区域和外腔区域只有扩散传质，而实际上膜厚区域和外腔区域也会发生对流传质。Lorenzin等[14]研究了中等分子截留量中空纤维膜对流过滤传质行为，发现在外腔不补充额外营养液情况下仍具有很高的对流传质效率，对流传质是物质传输的主要方式；但Macias等[15]通过半经验的方法研究却发现中等分子毒素主要是依靠扩散传质的方式去除的，得到相反的结果。可见，中空纤维膜组件内的传质行为还存在争议。
在中空纤维膜组件整体数值模拟方面，全组件计算流体动力学（computational fluid dynamics，CFD）模拟在中空纤维膜微/超滤等水处理方面有很多丰富的研究成果[16]，但由于人工肝中空纤维膜组件尺寸更小，流动问题更加复杂，故人工肝全组件CFD研究较少。Ding等[17-18]近十多年发展了一系列多孔介质模型，对人工肾中空纤维膜全组件流动传质进行CFD模拟。Sangeetha等[19]采用有限软件COMSOL对透析型中空纤维管的内腔和外腔流动进行了有限元分析，对比分析了直管和波浪形弯管的流动行为，发现弯管清除毒素的效果更好。Menshutina等[20]基于CFD的方法研究中空纤维膜组件内的流动行为，得到了沿纤维管方向的压力分布和速度分布，给出了营养液的最佳流量值。
由上述研究可以看到，目前并不十分清楚人工肝中空纤维膜组件的流动传质效率，基本上是结合实验和经验性或半经验性的数学模型，获得半经验性的结论。由于人工肝中流动问题的复杂性，关于中空纤维膜全组件传质行为的研究还有待进一步深入。本文通过建立人工肝中空纤维全组件的流动传质数值模型，研究小分子胆红素与大分子牛血清白蛋白（bovine serum albumin，BSA）在中空纤维膜组件内的流动传质规律，以期获得人工肝中的准确传质行为。
1 数值方法
1.1 CFD几何模型
人工肝中空纤维膜组件的结构示意图如图1左图所示，胆红素或BSA溶液在管程内流动，营养液（本文为纯水）在壳程内流动，两者流动方向相反，纤维膜区域为管程和壳程信息交换的界面，胆红素或者BSA等溶质分子从管程穿过纤维膜进入壳程被清除。根据中空纤维管的排列方式，可将其简化为一个代表性的单元，见横截面图中的虚线圆单元，该代表性单元包括中空纤维管内腔的管程区域、纤维膜膜厚区域以及中空纤维管外腔的壳程区域，三个区域的半径分别为中空纤维管的内侧半径
、外侧半径
和壳程半径
，虚线以外的壳程流动对虚线圆单元内的流动无影响。
图1
人工肝中空纤维膜组件结构示意图（左）与流动传质的CFD几何模型（右）
Figure1.
Schematic of hollow fiber membrane module (Left) and CFD model (Right) for artificial liver
根据代表性单元，可以构建出如图1右图所示的CFD几何模型。管程入口和壳程入口均采用流量入口边界，管程出口和壳程出口均采用压力出口边界，膜厚区域两侧也即中空纤维管内壁和外壁为内部边界（interior），壳程外表面为对称边界（symmetry），两端面膜厚区域（图中灰色区域）为壁面边界，并采用多孔介质模型对多孔纤维膜厚区域进行模拟，通过设置黏性阻力系数、惯性阻力系数以及纤维膜的孔隙率等参数来准确模拟纤维膜厚区域的传质行为。
1.2 控制方程
管程和壳程的连续性方程如下所示[21]：
 |
式中
、
、
分别为x、y、z方向上的速度分量，
指跨膜流体产生的对流源相，
指膜的孔隙率。
管程和壳程内流动为定常流，其动量守恒方程如下[22]：
 |
 |
 |
式中
指模拟多孔膜引入的阻力源相。
对于管程以及壳程通道内的溶质浓度可以由以下方程控制[22]：
 |
式中
是溶质浓度，
是溶质扩散系数，
是膜表面所产生的溶质源相，方程连接了管程和壳程溶质的体积流量，至关重要。对于溶质源相由以下关系来确定：
 |
式中
指溶质跨膜通量，
指纤维膜厚度。
K-K方程是由Kedem和Katchalsky开发出来用于描述多孔膜内的传质行为[23]。传统的K-K方程忽略了管程和壳程两侧的质量传输阻力，因此用总传质系数取代溶质渗透率来修正K-K方程，其表达式为：
 |
 |
式中
指溶液穿过膜的对流速度；
指多孔膜的水力渗透率；
指穿过多孔膜的压力差；
指膜两侧的浓度差；σ指溶质的反射系数；R指气体常数；T指绝对温度；
指穿过多孔膜的溶质通量；
指多孔膜的孔内平均溶质浓度；
指溶质扩散穿过多孔膜的总传质系数。溶质的反射系数与溶质穿过膜总传质系数和膜孔直径以及溶质分子大小有关，Anderson等[24-25]给出了孔径与颗粒尺寸之间的关系：
 |
式中
为多孔膜的膜孔半径；
指溶质分子的斯托克斯半径；对于小颗粒，当膜孔半径约20 nm或更大时，σ值非常小可忽略[22]。因此，式(7)可以进一步简化为：
 |
对膜孔进行模拟时采用的是管道式作为孔隙的模型，目前多数文献采用圆形直管道孔来模拟膜孔，但实际很多膜孔可能是弯曲的。因此，为了考虑中空纤维膜孔的曲率，本文数值模型中耦合了弯曲孔隙扩散模型来修正溶质的跨膜输运，该模型中水力渗透系数
和扩散渗透率
的计算公式如下所示[26]：
 |
 |
 |
 |
 |
 |
式中
指溶质在多孔膜的有效扩散率，它明显小于溶质的体积扩散率
。摩擦系数
用来描述溶质分子与多孔膜的膜孔壁面之间的摩擦；
指扩散作用下孔隙入口的空间位阻因子；τ指膜孔的弯曲因子。
对于不同的物质扩散系数是不同的，扩散系数对模拟结果至关重要，模拟所用生物材料的扩散系数表达式为[27]：
 |
此处MW是相应溶质的分子量。根据文献[9]和[28]，本文胆红素的分子量为584.66 Da，斯托克斯半径
为0.5 nm，BSA分子量为67 000 Da，斯托克斯半径
为3.61 nm。
多孔介质模型是通过在多孔区域以动量源项的形式引入额外的流动阻力描述的。通过在多孔介质区域的动量方程中添加一个动量源相，该源相由两部分组成，即黏性损失项（Darcy）和惯性损失项，如以下方程所示：
 |
式中，
是第i个（x、y或z方向）动量方程中的源项，
是速度大小的绝对值；
是粘度。增加的源相有助于描写多孔介质单元动量梯度，产生一个与流体速度成正比的压力降。
对于简单、均匀的多孔介质，可将上述方程简化为：
 |
式中
为多孔介质的渗透率，
称为黏性阻力系数，
称为惯性阻力系数，在各向异性多孔介质中，若某个方向阻力系数远远大于其他方向，无需设置很大的值，只需将其阻力系数设置为超过其他方向的2～3个量级。对于通过多孔介质的层流，常数
可以忽略不计，此时压降正比于速度。忽略对流加速与扩散项，可以用Darcy定律描绘多孔介质的流动：
 |
管程内流动被认为是完全发展的层流，管程入口处的速度分布为：
 |
式中
是管程内流体的体积流量，
是中空纤维管的内侧半径，
是中空纤维膜组件内纤维管的数量。
通过计算管程进出口溶质的平均浓度确定中空纤维膜组件清除速率，其计算公式为[29]：
 |
式中
指管程入口流量，
指管程出口流量，
指管程入口溶质平均浓度，
指管程出口溶质平均浓度。
采用ANSYS-FLUENT求解器对上述控制方程进行求解，分析胆红素以及BSA的传质行为。参照文献[30]的中空纤维膜组件的结构和操作参数，本文模拟中取管程入口流量为300 mL/min，壳程入口流量为500 mL/min，中空纤维管数量为12 000，其他参数分别取：中空纤维管长度L = 270 mm，管程半径
，外侧半径
，壳程半径
，纤维膜的孔隙率
，膜孔直径
（对应膜孔半径
），膜孔弯曲因子
。考虑到人体内白蛋白的浓度以及患者体内胆红素的浓度，本文模拟取胆红素和BSA入口浓度为0.000 1 mol/L和0.000 597 mol/L。
2 结果与讨论
2.1 网格划分敏感性测试
采用ICEM对CFD几何模型进行网格划分，网格划分采用六面体网格，因纤维膜区域以及膜表面与管程和壳程交界处流动的复杂性，需要对纤维膜边界处以及纤维膜区域的网格进行加密（具体网格划分图参见附件1）。采用三种不同数量的网格对中空纤维膜组件进行数值模拟，选用管程出口处沿管程径向上的胆红素浓度分布进行比较，进行网格划分敏感性测试。图2给出了不同网格数量时胆红素沿管程径向的浓度分布图，可以看到网格数量612 864和1 677 900的浓度分布几乎一致。为排除网格数量因素对本文结果的影响，同时节约计算资源，本文选用612 864个网格对中空纤维膜组件进行求解。
图2
不同网格数量下胆红素在管程出口处浓度分布
Figure2.
The bilirubin concentration distributions at the tube-side outlet for different total elements
2.2 中空纤维膜组件内胆红素传质行为
采用本文的数值模型分析胆红素溶液在中空纤维膜组件内流动时的压力分布与浓度分布。图3为管程入口流量为300 mL/min、壳程入口流量为500 mL/min时中空纤维膜组件内中截面处的压力分布，管程与壳程内的流动互为逆流，可以看出管程或壳程的入口处压力最大，并沿轴向方向不断减小。因此在纤维膜两侧存在跨膜压差，跨膜压差会导致在纤维膜两侧发生对流传质。
图3
中空纤维膜组件内压力分布（中截面）
Figure3.
Pressure distribution in hollow fiber membrane module (middle section)
图4为中空纤维膜组件内胆红素浓度分布以及5个横断面处的浓度分布图。可以看到，胆红素溶液在管程入口处进入流体域，但在壳程以及纤维膜上均有出现，其浓度随着轴向距离的增加而不断减小，这说明胆红素发生了从管程到壳程的传递，且胆红素的清除主要发生在管程入口附近。一方面，由图3可知纤维膜两侧存在跨膜压差会导致胆红素发生对流传质，且管程入口附近压力大于壳程出口附近压力，故胆红素的对流传质主要发生在管程入口附近（此处跨膜压差最大）。另一方面，管程与壳程侧也存在胆红素浓度梯度，会发生扩散传质。由图4可以看到在壳程入口附近，膜厚区域和壳程也有胆红素出现，壳程入口附近压力大于管程出口附近压力，也即壳程入口附近不会发生胆红素从管程到壳程的对流传质，这说明胆红素也会以扩散传质方式从管程传递到壳程。
图4
中空纤维膜组件内胆红素浓度分布
Figure4.
The bilirubin concentration distribution in hollow fiber membrane module
图5和图6分别为不同管程流量（对应壳程流量为500 mL/min）、不同壳程流量（对应管程流量为300 mL/min）、不同中空纤维管长度时的胆红素清除速率和中截面处浓度分布。从图5a可以看到，随着管程流量增大，胆红素的清除速率逐渐提升，胆红素会通过对流和扩散两种传质方式从管程经纤维膜传递到壳程。由于流量与压力降成正比，因此管程流量越大，管程入口附近处跨膜压差就越大，对流传质就越强，胆红素的清除速率就越大。同时由于管程流量提高也会提高纤维膜两侧的浓度梯度，使得胆红素的扩散传质也得到了提升。由图6a可见胆红素在管程出口的浓度不断增加，这是因为管程流量增加使得胆红素单位时间内进入管程的总量增加了，部分胆红素还未以对流与扩散作用穿过纤维膜便流出管程。
图5
胆红素清除速率
a. 不同管程流量；b. 不同壳程流量；c. 不同中空纤维管长度
Figure5. The bilirubin clearance ratesa. varying tube-side flow rates; b. varying shell-side flow rates; c. hollow fiber lengths
图6
胆红素浓度分布图
a. 不同管程流量；b. 不同壳程流量；c. 不同中空纤维管长度
Figure6. The bilirubin concentration distributionsa. varying tube-side flow rates; b. varying shell-side flow rates; c. hollow fiber lengths
从图5b可以看到，随着壳程流量增加，胆红素清除率先快速增加，然后缓慢增加至一渐近值。这是因为随着壳程流量增大，壳程中胆红素被带走得比较快，进而提高了纤维膜两侧的浓度梯度，所以清除速率得到提高；另一方面随着壳程流量增大，壳程的压力增加使得管程与壳程的压力差减小，胆红素对流通量减少。随着壳程流量增大，胆红素的对流通量越来越少，而扩散通量也接近极限值。从图6b胆红素浓度分布也可以看到，随着壳程流量增加，管程出口胆红素浓度一开始减小幅度比较大，然后减小幅度变缓慢，胆红素清除速率达到渐近值。
从图5c和6c可以看到，随着中空纤维管长度增加，胆红素的清除速率增大。一方面由于中空纤维管长度的增加导致溶质与膜的接触面积增加，增强了扩散传质；另一方面由于纤维管长度越长，管程和壳程的压力降也越大，则跨膜压差也越大，增强了对流传质。
2.3 中空纤维膜组件内BSA传质行为
同样分析了管程流量、壳程流量以及中空纤维管长度等参数对大分子BSA在中空纤维膜组件内传质行为的影响，模拟参数设置同胆红素一样。从图7a可以到，随着管程流量增大，BSA的清除速率逐渐提升；从图8a浓度分布也可以看到，在管程入口附近纤维膜区域与壳程区域胆红素浓度随管程流量增大而增加。与图5a胆红素的清除速率对比，BSA的清除速率明显较小，这是由于BSA分子尺寸较大，根据式(12)～(16)可以得到BSA的摩擦系数
以及空间位阻因子
值较小，进而BSA的扩散传质较少，所以BSA清除速率的提升不如胆红素。这也说明了BSA的传质方式以对流传质为主，而胆红素的传质方式是对流传质和扩散传质二者兼具。
图7
BSA清除速率
a. 不同管程流量；b. 不同壳程流量；c. 不同中空纤维管长度
Figure7. The BSA clearance ratesa. varying tube-side flow rates; b. varying shell-side flow rates; c. hollow fiber lengths
图8
BSA浓度分布
a. 不同管程流量；b. 不同壳程流量；c. 不同中空纤维管长度
Figure8. The BSA concentration distributionsa. varying tube-side flow rates; b. varying shell-side flow rates; c. hollow fiber lengths
由图7b和图8b可以看到，随着壳程流量增大，在管程入口附近纤维膜区域与壳程区域BSA浓度也逐渐减小，管程出口处BSA浓度逐渐增大，BSA的清除速率逐渐降低。根据前述结果得知，大分子BSA主要以对流作用进行传质，而壳程流量增大是减小了管程向壳程的跨膜压差，因此BSA的清除速率随着壳程流量增大而逐渐降低。尽管纤维膜两侧也存在BSA浓度梯度，但BSA的清除速率还是逐渐降低，这进一步说明了BSA的传质行为是以对流作用为主。
从图7c和图8c可以看到，随着中空纤维管长度增加，在管程入口附近纤维膜区域与壳程区域BSA浓度也逐渐增加，管程出口处BSA浓度逐渐减小，BSA的清除速率逐渐提升。由于BSA以对流传质为主，而胆红素兼具对流和扩散传质，因此BSA的清除速率比胆红素的清除速率低，由图5c和图7c可印证这一点。
3 结论
本文建立了人工肝中空纤维膜组件内流动传质数值模型，研究了小分子胆红素和大分子BSA在组件内的传质行为。对流作用和扩散作用控制着胆红素的传质行为，而BSA则以对流传质为主。增大管程流量提高了溶质对流传质以及扩散传质作用，从而可以提高胆红素与BSA的清除速率。增加壳程流量提高了溶质的扩散传质作用但却降低了溶质的对流传质作用，因此随着壳程流量的增大，胆红素的清除速率先快速提升然后缓慢上升达到一渐近值，而BSA的清除速率则随壳程流量增大而逐渐减小。增加中空纤维管长度可以提高溶质的对流传质以及扩散传质强度，胆红素与BSA的清除速率均随着中空纤维管长度的增加而得到提升。本研究可以指导人工肝器械中空纤维膜组件的结构设计和工艺参数设计，有助于深入了解中空纤维膜组件内不同溶质的传质机制。本文并没有考虑纤维膜表面的孔径大于膜内的孔径以及浓差极化等情况的影响，未来可开展膜厚区域的梯度多孔介质传质模型研究以及考虑浓差极化的影响，这将有助于更加准确地获取中空纤维膜组件传质规律。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突。
作者贡献声明：王子恒负责论文的设计模拟与计算、数据分析，许少锋负责本研究的构思和设计以及完善了论文的修改。论文的初稿由王子恒和许少锋撰写，余一帆完成文章最终校对，陆俊杰和张学昌对文章数据处理与撰写提供了宝贵的建议。
本文附件见本刊网站的电子版本（biomedeng.cn）。
0 引言
我国病毒性肝炎发病率很高，仅乙型肝炎一项，病毒感染者就有近1亿人，其中约7 000万为慢性乙型肝炎患者，约占全世界慢性乙型肝炎患者人数的1/5[1-3]。乙型肝炎病毒（hepatitis B virus，HBV）感染也是一个严重的世界性公共卫生问题，据世界卫生组织报道，全球70多亿人口中，超过20亿有乙型肝炎病毒感染证据，每年有110万人死于HBV和丙型肝炎病毒感染[4]。我国每年约有38.3万人死于肝癌，占全球肝癌死亡病例数的51%，全球每30 s就有1人死于病毒性肝炎相关疾病[5-6]。人工肝通过一个体外的机械、物理、生物装置，驱动患者血液流经体外人工肝器械，担负起暂时辅助严重病变肝脏的功能，清除各种有害物质，代偿肝脏的代谢功能，直至人体自身肝脏功能恢复或进行肝脏移植[7]。因此，进行人工肝技术的研究和开发具有重要的社会价值和科学意义。
目前关于人工肝中空纤维膜组件内的物质转运效率还不十分清楚，科研工作者对中空纤维膜组件的传质行为开展了大量的实验和数值模拟研究。Ilaria等[8]建立了对流强化型的中空纤维膜传质数学模型，得到了葡萄糖和氧气在中空纤维组件内的传输行为。Curcio等[9]结合实验和理论研究了聚醚砜膜和聚砜膜两种中空纤维膜组件的对流传质和扩散传质机制。Raff等[10]分析了人工肝中白蛋白修饰膜的传质特点，发现对流传质和扩散传质兼而有之；之后他们课题组的Donato等[11]又基于有限元方法建立了中空纤维透析膜的对流传质模型，结果显示增大中空纤维膜的压力系数和填充密度可以增强纤维膜的对流传质效果。Nedredal等[12]利用实验的方法研究了胆红素和免疫球蛋白在人工肝中空纤维组件内的传输行为，发现只有在强对流作用下才具有较强的毒素清除效果。Sorrell等[13]研究了醋氨酚肝毒素和氧气在中空纤维组件的内腔、纤维膜膜孔、外腔的传输行为，但该模型假设中空纤维管内物质传输为对流和扩散传质，膜厚区域和外腔区域只有扩散传质，而实际上膜厚区域和外腔区域也会发生对流传质。Lorenzin等[14]研究了中等分子截留量中空纤维膜对流过滤传质行为，发现在外腔不补充额外营养液情况下仍具有很高的对流传质效率，对流传质是物质传输的主要方式；但Macias等[15]通过半经验的方法研究却发现中等分子毒素主要是依靠扩散传质的方式去除的，得到相反的结果。可见，中空纤维膜组件内的传质行为还存在争议。
在中空纤维膜组件整体数值模拟方面，全组件计算流体动力学（computational fluid dynamics，CFD）模拟在中空纤维膜微/超滤等水处理方面有很多丰富的研究成果[16]，但由于人工肝中空纤维膜组件尺寸更小，流动问题更加复杂，故人工肝全组件CFD研究较少。Ding等[17-18]近十多年发展了一系列多孔介质模型，对人工肾中空纤维膜全组件流动传质进行CFD模拟。Sangeetha等[19]采用有限软件COMSOL对透析型中空纤维管的内腔和外腔流动进行了有限元分析，对比分析了直管和波浪形弯管的流动行为，发现弯管清除毒素的效果更好。Menshutina等[20]基于CFD的方法研究中空纤维膜组件内的流动行为，得到了沿纤维管方向的压力分布和速度分布，给出了营养液的最佳流量值。
由上述研究可以看到，目前并不十分清楚人工肝中空纤维膜组件的流动传质效率，基本上是结合实验和经验性或半经验性的数学模型，获得半经验性的结论。由于人工肝中流动问题的复杂性，关于中空纤维膜全组件传质行为的研究还有待进一步深入。本文通过建立人工肝中空纤维全组件的流动传质数值模型，研究小分子胆红素与大分子牛血清白蛋白（bovine serum albumin，BSA）在中空纤维膜组件内的流动传质规律，以期获得人工肝中的准确传质行为。
1 数值方法
1.1 CFD几何模型
人工肝中空纤维膜组件的结构示意图如图1左图所示，胆红素或BSA溶液在管程内流动，营养液（本文为纯水）在壳程内流动，两者流动方向相反，纤维膜区域为管程和壳程信息交换的界面，胆红素或者BSA等溶质分子从管程穿过纤维膜进入壳程被清除。根据中空纤维管的排列方式，可将其简化为一个代表性的单元，见横截面图中的虚线圆单元，该代表性单元包括中空纤维管内腔的管程区域、纤维膜膜厚区域以及中空纤维管外腔的壳程区域，三个区域的半径分别为中空纤维管的内侧半径、外侧半径和壳程半径，虚线以外的壳程流动对虚线圆单元内的流动无影响。
图1
人工肝中空纤维膜组件结构示意图（左）与流动传质的CFD几何模型（右）
Figure1.
Schematic of hollow fiber membrane module (Left) and CFD model (Right) for artificial liver
根据代表性单元，可以构建出如图1右图所示的CFD几何模型。管程入口和壳程入口均采用流量入口边界，管程出口和壳程出口均采用压力出口边界，膜厚区域两侧也即中空纤维管内壁和外壁为内部边界（interior），壳程外表面为对称边界（symmetry），两端面膜厚区域（图中灰色区域）为壁面边界，并采用多孔介质模型对多孔纤维膜厚区域进行模拟，通过设置黏性阻力系数、惯性阻力系数以及纤维膜的孔隙率等参数来准确模拟纤维膜厚区域的传质行为。
1.2 控制方程
管程和壳程的连续性方程如下所示[21]：
|
式中、、分别为x、y、z方向上的速度分量，指跨膜流体产生的对流源相，指膜的孔隙率。
管程和壳程内流动为定常流，其动量守恒方程如下[22]：
|
|
|
式中指模拟多孔膜引入的阻力源相。
对于管程以及壳程通道内的溶质浓度可以由以下方程控制[22]：
|
式中是溶质浓度，是溶质扩散系数，是膜表面所产生的溶质源相，方程连接了管程和壳程溶质的体积流量，至关重要。对于溶质源相由以下关系来确定：
|
式中指溶质跨膜通量，指纤维膜厚度。
K-K方程是由Kedem和Katchalsky开发出来用于描述多孔膜内的传质行为[23]。传统的K-K方程忽略了管程和壳程两侧的质量传输阻力，因此用总传质系数取代溶质渗透率来修正K-K方程，其表达式为：
|
|
式中指溶液穿过膜的对流速度；指多孔膜的水力渗透率；指穿过多孔膜的压力差；指膜两侧的浓度差；σ指溶质的反射系数；R指气体常数；T指绝对温度；指穿过多孔膜的溶质通量；指多孔膜的孔内平均溶质浓度；指溶质扩散穿过多孔膜的总传质系数。溶质的反射系数与溶质穿过膜总传质系数和膜孔直径以及溶质分子大小有关，Anderson等[24-25]给出了孔径与颗粒尺寸之间的关系：
|
式中为多孔膜的膜孔半径；指溶质分子的斯托克斯半径；对于小颗粒，当膜孔半径约20 nm或更大时，σ值非常小可忽略[22]。因此，式(7)可以进一步简化为：
|
对膜孔进行模拟时采用的是管道式作为孔隙的模型，目前多数文献采用圆形直管道孔来模拟膜孔，但实际很多膜孔可能是弯曲的。因此，为了考虑中空纤维膜孔的曲率，本文数值模型中耦合了弯曲孔隙扩散模型来修正溶质的跨膜输运，该模型中水力渗透系数和扩散渗透率的计算公式如下所示[26]：
|
|
|
|
|
|
式中指溶质在多孔膜的有效扩散率，它明显小于溶质的体积扩散率。摩擦系数用来描述溶质分子与多孔膜的膜孔壁面之间的摩擦；指扩散作用下孔隙入口的空间位阻因子；τ指膜孔的弯曲因子。
对于不同的物质扩散系数是不同的，扩散系数对模拟结果至关重要，模拟所用生物材料的扩散系数表达式为[27]：
|
此处MW是相应溶质的分子量。根据文献[9]和[28]，本文胆红素的分子量为584.66 Da，斯托克斯半径为0.5 nm，BSA分子量为67 000 Da，斯托克斯半径为3.61 nm。
多孔介质模型是通过在多孔区域以动量源项的形式引入额外的流动阻力描述的。通过在多孔介质区域的动量方程中添加一个动量源相，该源相由两部分组成，即黏性损失项（Darcy）和惯性损失项，如以下方程所示：
|
式中，是第i个（x、y或z方向）动量方程中的源项，是速度大小的绝对值；是粘度。增加的源相有助于描写多孔介质单元动量梯度，产生一个与流体速度成正比的压力降。
对于简单、均匀的多孔介质，可将上述方程简化为：
|
式中为多孔介质的渗透率，称为黏性阻力系数，称为惯性阻力系数，在各向异性多孔介质中，若某个方向阻力系数远远大于其他方向，无需设置很大的值，只需将其阻力系数设置为超过其他方向的2～3个量级。对于通过多孔介质的层流，常数可以忽略不计，此时压降正比于速度。忽略对流加速与扩散项，可以用Darcy定律描绘多孔介质的流动：
|
管程内流动被认为是完全发展的层流，管程入口处的速度分布为：
|
式中是管程内流体的体积流量，是中空纤维管的内侧半径，是中空纤维膜组件内纤维管的数量。
通过计算管程进出口溶质的平均浓度确定中空纤维膜组件清除速率，其计算公式为[29]：
|
式中指管程入口流量，指管程出口流量，指管程入口溶质平均浓度，指管程出口溶质平均浓度。
采用ANSYS-FLUENT求解器对上述控制方程进行求解，分析胆红素以及BSA的传质行为。参照文献[30]的中空纤维膜组件的结构和操作参数，本文模拟中取管程入口流量为300 mL/min，壳程入口流量为500 mL/min，中空纤维管数量为12 000，其他参数分别取：中空纤维管长度L = 270 mm，管程半径，外侧半径，壳程半径，纤维膜的孔隙率，膜孔直径（对应膜孔半径），膜孔弯曲因子。考虑到人体内白蛋白的浓度以及患者体内胆红素的浓度，本文模拟取胆红素和BSA入口浓度为0.000 1 mol/L和0.000 597 mol/L。
2 结果与讨论
2.1 网格划分敏感性测试
采用ICEM对CFD几何模型进行网格划分，网格划分采用六面体网格，因纤维膜区域以及膜表面与管程和壳程交界处流动的复杂性，需要对纤维膜边界处以及纤维膜区域的网格进行加密（具体网格划分图参见附件1）。采用三种不同数量的网格对中空纤维膜组件进行数值模拟，选用管程出口处沿管程径向上的胆红素浓度分布进行比较，进行网格划分敏感性测试。图2给出了不同网格数量时胆红素沿管程径向的浓度分布图，可以看到网格数量612 864和1 677 900的浓度分布几乎一致。为排除网格数量因素对本文结果的影响，同时节约计算资源，本文选用612 864个网格对中空纤维膜组件进行求解。
图2
不同网格数量下胆红素在管程出口处浓度分布
Figure2.
The bilirubin concentration distributions at the tube-side outlet for different total elements
2.2 中空纤维膜组件内胆红素传质行为
采用本文的数值模型分析胆红素溶液在中空纤维膜组件内流动时的压力分布与浓度分布。图3为管程入口流量为300 mL/min、壳程入口流量为500 mL/min时中空纤维膜组件内中截面处的压力分布，管程与壳程内的流动互为逆流，可以看出管程或壳程的入口处压力最大，并沿轴向方向不断减小。因此在纤维膜两侧存在跨膜压差，跨膜压差会导致在纤维膜两侧发生对流传质。
图3
中空纤维膜组件内压力分布（中截面）
Figure3.
Pressure distribution in hollow fiber membrane module (middle section)
图4为中空纤维膜组件内胆红素浓度分布以及5个横断面处的浓度分布图。可以看到，胆红素溶液在管程入口处进入流体域，但在壳程以及纤维膜上均有出现，其浓度随着轴向距离的增加而不断减小，这说明胆红素发生了从管程到壳程的传递，且胆红素的清除主要发生在管程入口附近。一方面，由图3可知纤维膜两侧存在跨膜压差会导致胆红素发生对流传质，且管程入口附近压力大于壳程出口附近压力，故胆红素的对流传质主要发生在管程入口附近（此处跨膜压差最大）。另一方面，管程与壳程侧也存在胆红素浓度梯度，会发生扩散传质。由图4可以看到在壳程入口附近，膜厚区域和壳程也有胆红素出现，壳程入口附近压力大于管程出口附近压力，也即壳程入口附近不会发生胆红素从管程到壳程的对流传质，这说明胆红素也会以扩散传质方式从管程传递到壳程。
图4
中空纤维膜组件内胆红素浓度分布
Figure4.
The bilirubin concentration distribution in hollow fiber membrane module
图5和图6分别为不同管程流量（对应壳程流量为500 mL/min）、不同壳程流量（对应管程流量为300 mL/min）、不同中空纤维管长度时的胆红素清除速率和中截面处浓度分布。从图5a可以看到，随着管程流量增大，胆红素的清除速率逐渐提升，胆红素会通过对流和扩散两种传质方式从管程经纤维膜传递到壳程。由于流量与压力降成正比，因此管程流量越大，管程入口附近处跨膜压差就越大，对流传质就越强，胆红素的清除速率就越大。同时由于管程流量提高也会提高纤维膜两侧的浓度梯度，使得胆红素的扩散传质也得到了提升。由图6a可见胆红素在管程出口的浓度不断增加，这是因为管程流量增加使得胆红素单位时间内进入管程的总量增加了，部分胆红素还未以对流与扩散作用穿过纤维膜便流出管程。
图5
胆红素清除速率
a. 不同管程流量；b. 不同壳程流量；c. 不同中空纤维管长度
Figure5. The bilirubin clearance ratesa. varying tube-side flow rates; b. varying shell-side flow rates; c. hollow fiber lengths
图6
胆红素浓度分布图
a. 不同管程流量；b. 不同壳程流量；c. 不同中空纤维管长度
Figure6. The bilirubin concentration distributionsa. varying tube-side flow rates; b. varying shell-side flow rates; c. hollow fiber lengths
从图5b可以看到，随着壳程流量增加，胆红素清除率先快速增加，然后缓慢增加至一渐近值。这是因为随着壳程流量增大，壳程中胆红素被带走得比较快，进而提高了纤维膜两侧的浓度梯度，所以清除速率得到提高；另一方面随着壳程流量增大，壳程的压力增加使得管程与壳程的压力差减小，胆红素对流通量减少。随着壳程流量增大，胆红素的对流通量越来越少，而扩散通量也接近极限值。从图6b胆红素浓度分布也可以看到，随着壳程流量增加，管程出口胆红素浓度一开始减小幅度比较大，然后减小幅度变缓慢，胆红素清除速率达到渐近值。
从图5c和6c可以看到，随着中空纤维管长度增加，胆红素的清除速率增大。一方面由于中空纤维管长度的增加导致溶质与膜的接触面积增加，增强了扩散传质；另一方面由于纤维管长度越长，管程和壳程的压力降也越大，则跨膜压差也越大，增强了对流传质。
2.3 中空纤维膜组件内BSA传质行为
同样分析了管程流量、壳程流量以及中空纤维管长度等参数对大分子BSA在中空纤维膜组件内传质行为的影响，模拟参数设置同胆红素一样。从图7a可以到，随着管程流量增大，BSA的清除速率逐渐提升；从图8a浓度分布也可以看到，在管程入口附近纤维膜区域与壳程区域胆红素浓度随管程流量增大而增加。与图5a胆红素的清除速率对比，BSA的清除速率明显较小，这是由于BSA分子尺寸较大，根据式(12)～(16)可以得到BSA的摩擦系数以及空间位阻因子值较小，进而BSA的扩散传质较少，所以BSA清除速率的提升不如胆红素。这也说明了BSA的传质方式以对流传质为主，而胆红素的传质方式是对流传质和扩散传质二者兼具。
图7
BSA清除速率
a. 不同管程流量；b. 不同壳程流量；c. 不同中空纤维管长度
Figure7. The BSA clearance ratesa. varying tube-side flow rates; b. varying shell-side flow rates; c. hollow fiber lengths
图8
BSA浓度分布
a. 不同管程流量；b. 不同壳程流量；c. 不同中空纤维管长度
Figure8. The BSA concentration distributionsa. varying tube-side flow rates; b. varying shell-side flow rates; c. hollow fiber lengths
由图7b和图8b可以看到，随着壳程流量增大，在管程入口附近纤维膜区域与壳程区域BSA浓度也逐渐减小，管程出口处BSA浓度逐渐增大，BSA的清除速率逐渐降低。根据前述结果得知，大分子BSA主要以对流作用进行传质，而壳程流量增大是减小了管程向壳程的跨膜压差，因此BSA的清除速率随着壳程流量增大而逐渐降低。尽管纤维膜两侧也存在BSA浓度梯度，但BSA的清除速率还是逐渐降低，这进一步说明了BSA的传质行为是以对流作用为主。
从图7c和图8c可以看到，随着中空纤维管长度增加，在管程入口附近纤维膜区域与壳程区域BSA浓度也逐渐增加，管程出口处BSA浓度逐渐减小，BSA的清除速率逐渐提升。由于BSA以对流传质为主，而胆红素兼具对流和扩散传质，因此BSA的清除速率比胆红素的清除速率低，由图5c和图7c可印证这一点。
3 结论
本文建立了人工肝中空纤维膜组件内流动传质数值模型，研究了小分子胆红素和大分子BSA在组件内的传质行为。对流作用和扩散作用控制着胆红素的传质行为，而BSA则以对流传质为主。增大管程流量提高了溶质对流传质以及扩散传质作用，从而可以提高胆红素与BSA的清除速率。增加壳程流量提高了溶质的扩散传质作用但却降低了溶质的对流传质作用，因此随着壳程流量的增大，胆红素的清除速率先快速提升然后缓慢上升达到一渐近值，而BSA的清除速率则随壳程流量增大而逐渐减小。增加中空纤维管长度可以提高溶质的对流传质以及扩散传质强度，胆红素与BSA的清除速率均随着中空纤维管长度的增加而得到提升。本研究可以指导人工肝器械中空纤维膜组件的结构设计和工艺参数设计，有助于深入了解中空纤维膜组件内不同溶质的传质机制。本文并没有考虑纤维膜表面的孔径大于膜内的孔径以及浓差极化等情况的影响，未来可开展膜厚区域的梯度多孔介质传质模型研究以及考虑浓差极化的影响，这将有助于更加准确地获取中空纤维膜组件传质规律。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突。
作者贡献声明：王子恒负责论文的设计模拟与计算、数据分析，许少锋负责本研究的构思和设计以及完善了论文的修改。论文的初稿由王子恒和许少锋撰写，余一帆完成文章最终校对，陆俊杰和张学昌对文章数据处理与撰写提供了宝贵的建议。
本文附件见本刊网站的电子版本（biomedeng.cn）。
