CHEMICAL ENGINEERING TRANSACTIONS VOL. 76, 2019 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Guest Editors: Petar S. Varbanov, Timothy G. Walmsley, Jiří J. Klemeš, Panos Seferlis Copyright © 2019, AIDIC Servizi S.r.l. ISBN 978-88-95608-73-0; ISSN 2283-9216 Fully Resolved Computational (CFD) and Experimental Analysis of Pressure Drop and Blood Gas Transport in a Hollow Fibre Membrane Oxygenator Module Michael Haraseka,*, Benjamin Lukitscha,c, Paul Eckera,b, Christoph Janeczekb,c, Martin Elenkovb,c, Thomas Keckb,c, Bahram Haddadia, Christian Jordana, Susanna Neudlc,d, Claus Krennc,d, Roman Ullrichc,d, Margit Gfoehlerb aInstitute of Chemical, Environmental and Bioscience Engineering, TU Wien, Getreidemarkt 9 1060 Vienna, Austria bInstitute of Engineering Design and Product Development, TU Wien, Getreidemarkt 9 1060 Vienna, Austria cCCORE Technology GmbH, Argentinierstraße 35 1040 Vienna, Austria dDepartment of Anaesthesia, Critical Care and Pain Medicine, MedUni Vienna, Währinger Gürtel 18 1090 Vienna, Austria michael.harasek@tuwien.ac.at Membrane oxygenators or artificial lungs have become a reliable and live saving clinical technique. To limit side effects on blood platelet parameters the high extrinsic surface introduced by the hollow fiber membrane packing has to be decreased. This adjustment must be accompanied with an improvement of gas exchange performance to guarantee sufficient blood gas transfer. Computational fluid dynamics (CFD) provides a spatial and temporal resolution of the membrane oxygenation process and enables systematic optimization of artificial lungs. In scope of this research a new CFD approach to examine the gas exchange performance of oxygenators was developed. Blood and sweep gas flow in the fiber packing as well as blood gas exchange through the membrane between blood and sweep fluid are fully resolved and simulated. The results were compared to in vitro experiments comprising determination of blood side pressure loss and CO2 exchange performance of a prototype membrane module. While flow simulations were conducted for the complete prototype module, gas exchange simulations were limited to a representative fiber packing segment. The presented simulation approach provides a basis for the design of future artificial lungs. 1. Introduction Artificial lungs or extracorporeal membrane oxygenators (ECMOs) are medical devices to supplement the respiratory function during cardiopulmonary bypass, or support patients suffering from respiratory failure (Federspiel and Henchir, 2008). In state-of-the-art oxygenators, blood gas exchange is enabled by hollow fiber membranes. In these devices blood flows on the shell side of the packing, while the fiber lumen is swept with O2. CO2 and O2 are exchanged through the membrane following the partial pressure gradient. To provide sufficient gas exchange a membrane area of up to 2.5 m2 is built into oxygenators. This large extrinsic surface area has serious effects on platelets and clotting factors (Jaffer et al., 2018). Therefore, the focus in optimization of membrane oxygenators lies in minimizing membrane area while maximizing the gas exchange rates. Using experimental methods for the necessary development iterations is usually very time consuming, expensive and limited to pointwise data. In contrast, computational fluid dynamics (CFD) can provide a rather fast and flexible platform for the investigation of different designs and operational parameters. CFD also delivers a spatial and temporal resolution of the process. Various CFD simulation approaches to describe the gas exchange process in oxygenators were proposed. Saving computational costs Zhang et al. (2007) modeled the fiber packing as homogenous porous medium using Darcy’s law to modify the momentum equation. Gas exchange was then implemented using a uniform source term in the packing. The source term is computed based on Sherwood correlations which in turn have to be determined experimentally for different packing densities and flow arrangements (Low et al., 2017). As the DOI: 10.3303/CET1976033 Paper Received: 15/03/2019; Revised: 14/04/2019; Accepted: 29/04/2019 Please cite this article as: Harasek M., Lukitsch B., Ecker P., Janeczek C., Elenkov M., Keck T., Haddadi B., Jordan C., Neudl S., Krenn C., Ullrich R., Gfoehler M., 2019, Fully Resolved Computational (CFD) and Experimental Analysis of Pressure Drop and Blood Gas Transport in a Hollow Fibre Membrane Oxygenator Module, Chemical Engineering Transactions, 76, 193-198 DOI:10.3303/CET1976033 193 blood flow through the fiber packing is not resolved, the main blood gas transport resistance in oxygenators – blood gas concentration polarization layer attached to the membrane surface (Federspiel and Henchir, 2008) – cannot be examined. Due to the high number of computational cells which are necessary to represent the geometry of a fiber packing, simulations that compute the flow and blood gas transport in the packing are computationally expensive. Additionally, the strong concentration gradients of the blood gases in the polarization layer make a high cell refinement aligned with the fiber surface mandatory (Dzhonova-Atanasova et al., 2018). Gas transport simulations are limited to microscopic geometries such as single fibers and periodic fiber arrangements (Low et al., 2017) or parts of packings (Hormes et al., 2011). To reduce computational effort and simulation complexity gas transport in the membrane and sweep gas flow are neglected. Instead a fixed uniform partial pressure is applied on the blood side fiber surfaces as boundary condition for the specie balances. Ambitions to simulate the gas transport in the lumen or the membrane are rare. To the authors best knowledge only one study was exploiting this research field. Taskin et al. (2010) additionally resolved the membrane wall and simulated the oxygen transport in a small packing of an experimental oxygenator. With this approach, they were capable to consider the membrane gas transport resistance and could show that the common boundary condition method leads to an overprediction of oxygen exchange rates by the factor of two due to the negligence of the membrane resistance. Moreover, gas transport in the sweep flow was not simulated. In this work a novel CFD approach is presented. In addition to blood side and membrane wall also the fiber lumen side is added to the simulation domain. Sweep gas and blood flow as well as gas transport in blood, membrane and sweep fluid are then simulated simultaneously. This allows to evaluate the CO2 enrichment of sweep flow and provides the possibility to assess different fiber geometries and membrane materials. In order to carry out these simulations the inhouse developed multi-region solver membrane Foam (Haddadi et al., 2018) was utilized. To limit computational costs, the simulated geometry was reduced to a representative packing segment of a prototype hollow fiber membrane module. To determine adequate velocity boundary conditions blood flow simulations of the whole prototype module were conducted in a preliminary study. Pressure drop in dependency from the blood flow rate as well as gas exchange performance of the membrane module were evaluated experimentally during in vitro tests using bovine blood. Results of both simulations are compared with the experimental data. 2. Materials and methods 2.1 In vitro tests CO2 removal performance tests of a prototype module were conducted in an in vitro recirculation loop. This loop consists of a centrifugal pump (BPX-80, Medtronic), an extracorporeal membrane oxygenator (ECMO Adult, Eurosets) and the prototype module (Figure 1a). Bovine blood was provided by the teaching and research farm of the University of Veterinary Medicine, Vienna. All tests were performed under the ethics proposal GZ: 66.009/0007/V/3b/2018. Blood (1,000 mL/min) was pumped through the shell side of the ECMO and the prototype module. A clamp-on ultrasonic flow probe (SONOFLOW CO.55/080) monitored the blood flow rate. Blood temperature was maintained at 37 °C using the ECMO heat exchanger. The ECMO fiber lumen were swept with a N2/CO2 saturation stream to enrich the blood with CO2, setting a pathological partial pressure level of 50 mmHg. CO2 was then removed from blood in the prototype module. Three blood samples were taken before and after the prototype module and analyzed with a blood gas analyzer (ABL500 FLEX, Radiometer Medical A/S). The prototype module was swept with pure O2 (1 L/min) for CO2 removal. To regulate the sweep gas flows of the ECMO and the prototype module, three mass flow controllers (GF40, Brooks) were facilitated. Outgoing sweep gas flow rates were recorded using a piston stroke volumetric measurement device (Defender 510, Bios DryCal). To determine the gas exchange performance, CO2 concentration in the outgoing sweep gas flow of the prototype module was measured (BINOS 100 M, Emerson). To gain the pressure characteristics of the modules all pressures were measured using miniaturized pressure transmitters (AMS 4711, Analog Microelectronics). Besides CO2 removal performance tests, blood side pressure loss over the prototype module was measured for different blood flows ranging from 1,200 to 2,900 mL/min. The membrane module (Figure 1b) contains a hollow fiber packing shaped in the form of a hollow cylinder. The packing has an outer diameter of 14 mm and an inner diameter of 6 mm. It has an active fiber length of 100 mm and is positioned coaxially in the module, shaping an outer ring gap of 3 mm. Blood inlet and outlet pipes are also positioned coaxially at the ends of the module. As flow guidance a plug placed at the middle of the module blocks the blood flow in the central channel, forcing a transverse flow through the packing. Approximately 410 fibers (Membrana Oxyplus® 90/200 PMP, 3M) were incorporated into the module. Tangential and radial spacing between the fibers measures 220 and 200 µm. 194 Figure 1: a) Scheme of in vitro recirculation loop with prototype module, oxygenator (Oxy) and ECMO pump. b) Blood guiding parts of the membrane module with cut through the membrane packing. 2.2 Computational Fluid Dynamics For all fluid dynamic simulations, the opensource toolbox OpenFOAM® 4.1 (ESI Group, Paris, France) was used. As a preliminary study, blood flow through the whole prototype module was simulated to gain pressure loss and velocity distribution. Therefore, the finite volume formulation of the steady incompressible Navier- Stokes equation was solved applying the PIMPLE algorithm. Due to the narrow space between fibers and a Reynolds number of about 1,600 in the inlet pipe the flow was assumed laminar and no turbulence model was applied. Second order discretization schemes (Van Leer) were utilized. A uniform inlet velocity at the beginning of the inlet pipe was set and varied to match inlet flow rates ranging from 800 to 3,000 mL/min. At all walls, including membrane surfaces a no-slip condition was applied. At the end of the outlet pipe a fixed uniform value of 0 Pa was set. For the remaining boundaries zero gradient conditions were imposed. While there are efforts to model the multiphase character of blood (Gonzalez & Rojas-Solorzano, 2016), a single-phase approach was chosen here to reduce computational costs. The influence of the shear thinning rheological behavior of blood on the flow field was examined by conducting simulations with two viscosity models, Newtonian and power law. For Newtonian viscosity, the kinematic viscosity of bovine blood was used. For the power law model coefficients for human blood were chosen. Flow and gas transport in a representative packing segment of the prototype module were simulated using membraneFoam. The multi-region solver was developed based on OpenFOAM®. The flows in the single regions which are separated by the membrane surface are computed separately but are coupled by transmembrane flux (Ji) which is implemented as volumetric source term at the membrane surface and calculated based on partial pressure difference (Δpi), membrane permeance (P) and membrane surface area (A) of the computational cell (Haddadi et al., 2018). Due to the generic architecture of membraneFoam the solver can be applied to any (environmental relevant) membrane processes including biogas upgrading, CO2 separation and pervaporation of ABE (acetone, butanol, ethanol) fermentation products. 𝐽𝑖 = 𝑃 ∙ 𝐴 ∙ Δ𝑝𝑖 (1) Blood region – gas transport simulation Solver, discretization schemes, turbulence model and boundary conditions for the continuous shell side blood region were set equally to the preliminary flow simulations described above. To model the CO2 transport, the different CO2-related species, dissolved CO2 and bicarbonate were summarized to one total CO2 specie. The CO2 partial pressure (𝑝𝐶𝑂2 ) was then calculated based on the concentration of this single specie (𝑐𝐶𝑂2,𝑡𝑜𝑡𝑎𝑙) (Svitek and Federspiel, 2008). 𝑐𝐶𝑂2,𝑡𝑜𝑡𝑎𝑙 = 𝑞 ∙ 𝑝𝐶𝑂2 𝑡 𝑤𝑖𝑡ℎ 𝑞 = 0.128 𝑎𝑛𝑑 𝑡 = 0.369 (2) Assuming a constant dissociation slope 𝜆 = 𝛿𝑐𝐶𝑂2 𝛿𝑝𝐶𝑂2⁄ at clinically relevant CO2 partial pressures of approximately 50 mmHg a diffusion coefficient of total CO2 (𝐷𝐶𝑂2,𝑡𝑜𝑡𝑎𝑙) can be expressed using the solubility of dissolved CO2 (𝛼𝐶𝑂2 ) and diffusion coefficients of dissolved CO2 (𝐷𝐶𝑂2 ) and bicarbonate (𝐷𝐻𝐶𝑂3− ): 𝐷𝐶𝑂2,𝑡𝑜𝑡𝑎𝑙 = 𝐷𝐻𝐶𝑂3− + (𝐷𝐶𝑂2 − 𝐷𝐻𝐶𝑂3− ) 𝛼𝐶𝑂2 𝜆 (3) Lumen and membrane region – gas transport simulation Compressible gas flow in the fiber lumen is simulated using the PIMPLE algorithm. To reduce computational cost the sweep fluid manifold and the sweep fluid collector at the end of the fibers were not included in the simulation domain. Each simulated fiber has a single inlet and outlet. Discretization and boundary conditions a) b) Plug as flow guidance 195 are set up equal to the blood region. Due to the small lumen diameter no turbulence model is applied. The CO2/O2 mixture is modeled as an ideal mixture. Gas viscosity was calculated based on Sutherland (Sutherland, 1893). The membrane region was set up similar to the lumen region. As the permeating gases are confined in the porous substructure of the membrane wall it was assumed that gas transport is solely of diffusive and not advective nature. Hence velocity in the membrane was set uniformly to zero. To account for the pore structure of the membrane the diffusion coefficient of the permeating species in the gas filled pores 𝐷𝑖,𝑝𝑜𝑟𝑒𝑠 was corrected according to Lu et al. (2008). In Eq(4) 𝜏 is the tortuosity and 𝜖 the porosity of the membrane wall. Diffusion through the membrane material was neglected because of the significantly lower diffusion coefficient. 𝐷𝑖 = 𝐷𝑖,𝑝𝑜𝑟𝑒𝑠 𝜖 𝜏 (4) 3. Results 3.1 In vitro results Blood side pressure loss in the prototype module is low and ranges from 40 mmHg at 850 mL/min to approximately 270 mmHg at 2,890 mL/min. The correlation between blood side pressure loss in the prototype module and blood flow rate is illustrated in Figure 2a. For a blood flow rate of 1,280 mL/min and an O2 sweep flow of 983 mL STP/min, a CO2 removal rate of 10.8 mL STP/min was recorded. This equals a specific transmembrane CO2 flux of 220 mL STP/min/m2. A hematocrit of approximately 30 % was determined for the bovine blood which matches well with literature data (Windberger et al., 2003) but is significantly lower than the average for human blood (45 %) (Galdi et al., 2008). Figure 2: a) Pressure drop in the prototype module in dependency of blood flow rate, determined by experiments (x) and CFD using different viscosity models: Newtonian (█) and power law (▲). b) Radial velocity within the packing determined with CFD at four different radial positions (ϕ) in dependency of the fibre packing length. 3.2 Computational fluid dynamics results Blood side pressure loss over the membrane module was computed based on flow simulations. Applying the power law with model coefficients for human blood leads to an underprediction of the pressure loss (approximately 12 mmHg at 800 mL/min, 20 mmHg at 1800 mL/min). Using Newtonian viscosity of bovine blood allows an accurate prediction of the pressure loss at flow rates lower than 1,800 mL/min (Figure 2a). A deviation of 11 % must be observed at a flow of 3000 mL/min. Based on the velocity distribution the radial velocity was computed. The average radial velocity amounts to 0.03 m/s. Radial velocities vary more before the plug (0 – 45 mm) than after the plug (55 – 100 mm) (Figure 2b). Therefore, the blood flow appears to be better distributed after the plug. This can also be seen in Figure 3 which shows the radial velocity plots at different longitudinal positions in the packing. Due to the more homogenous distribution after the plug (50 mm) the flow fields at 60 and 80 mm deviate less from each other than those at 20 and 40 mm. Figure 4 shows the results of the micro scale CO2 transport simulation. An ideal cross flow set-up for eight fibres in staggered and non-staggered arrangement was simulated. Since the macroscopic flow simulations show that the fibers are mostly positioned in the slipstream of each other this work is focusing on non-staggered arrangement only. Peak velocity (0.03 m/s) between fibres matches the radial velocity of the macroscopic flow simulation. The parabolic velocity profiles of the permeate gas flow show a maximum velocity of 2.5 m/s at permeate flow rates of 1 L/min. As reported in literature (Federspiel & Henchir, 2008), the concentration polarization layer is accountable for the main drop in driving force for transmembrane flux. 0 40 80 120 160 200 240 280 320 800 1800 2800 p re s s u re d ro p m e m b ra n e m o d u le [m m H g ] blood flow rate [mL/min] a) CFD Newtonian (Bovine) CFD Power Law (Human) Experimental Reynolds Number = 2300 0.0 0.1 0.2 0.3 0.4 0.5 0 10 20 30 40 50 60 70 80 90 100 ra d ia l v e lo c it y [ m /s ] fiber packing length [mm] ϕ = 45° ϕ = 135° ϕ = 225° ϕ = 315° p lu g p lu g a s f lo w g u id a n c e b) 196 Figure 3: Radial velocity profiles at different distances to the inlet pipe entry into the membrane packing. CO2 partial pressure in blood decreases from 50 mmHg outside the boundary layer to around 15 mmHg at the membrane surface. CO2 partial pressure drop over the selective membrane layer totals to roughly 12 mmHg. An additional partial pressure drop of 2 mmHg occurs in the membrane. No partial pressure gradient in radial direction can be observed in the lumen. Specific CO2 flux variates from approximately 1000 mL STP/min/m2 at the front of the first fibre to 600 - 700 mL STP/min/m2 at the fibre sides and declines to 200 mL STP/min/m2 at the fiber’s back side. The specific CO2 flux based on simulation results is 570 mL STP/min/m2. Figure 4: Velocity magnitude profiles (left) and CO2 partial pressure profiles (middle) in two quantitative scales applied uniformly to shell, membrane and lumen region. Specific CO2 flux on outer membrane surface (right). 4. Discussion A clear difference in blood flow pressure loss over the membrane can be observed for the two applied viscosity models. Using a power law viscosity model with coefficients for human blood gives a lower effective viscosity and causes an underprediction of the pressure loss while Newtonian viscosity of bovine blood predicts the pressure loss reasonably. The low average hematocrit of bovine blood (30 %) seems to allow the neglection of shear-thinning effects. The 11 % deviation from numerical to experimental determined pressure drop at blood flow rates of 2,960 mL/min could be explained by the inlet pipe Reynolds number of about 2,900. At this flow rates a turbulence model should be applied, however, flow is expected to turn laminar shortly after the pipe exit. Micro scale simulation of the CO2 transport overpredicts the specific flux by 350 mL STP/min/m2 (160 %). This overprediction was to be expected due to the idealized flow past a non-staggered packing. Furthermore, the CO2 solubility model used in the CFD simulations was developed for human blood and might not apply for bovine blood leading to additional deviations from experimental results. Comparing the effective viscosity of human (3.5 mPa s) and bovine blood (5.8 mPa s) indicates a presumably higher protein concentration in bovine blood plasma (Windberger et al., 2003). A higher viscosity increases the boundary layer built up at the membrane surface and could reduce gas exchange performance. This should be considered when determining the blood gas exchange performance of an oxygenator using non-human blood during in vitro trials. At higher protein concentrations CO2 species diffusion is reduced additionally (Geers and Gros, 2000). This further amplifies the built up of the concentration polarization layer at the membrane surface. While bovine blood viscosity was chosen for the simulations no data for CO2 diffusion coefficients in bovine blood plasma was available. 5. Conclusions A new CFD approach was developed which allows to determine the gas exchange performance of an oxygenator by fully resolving the blood gas transport through the membrane as well as flow and species Plug as flow guidance (see Figure 1b) r a d ia l v e lo c it y [ m /s ] velocity magnitude [m/s] Detail a a CO2 partial pressure [mmHg] b Detail b specific CO2 flux [ml STP/min/m2] 197 transport on blood and sweep gas side. In scope of this work the new CFD approach was applied to a packing segment of a prototype hollow fiber membrane module. The numerically determined specific transmembrane CO2 flux amounts to 570 mL STP/min/m2. The results were compared to in vitro trials in which the CO2 removal performance of the prototype module was determined for bovine blood. A specific CO2 flux of 220 mL STP/min/m2 was measured at blood flowrates of 1280 mL/min and an O2 sweep flow of 983 mL STP/min. The numerical overprediction (350 mL STP/min/m2) was to be expected due to the reduction of the module geometry to a packing segment with an idealized cross flow set-up. Additionally, the higher blood viscosity of bovine blood suggests that CO2 diffusion could be lower than in human blood. For further improvement of the simulation it is recommended to adjust solubility and diffusion models. By exploiting symmetries of the prototype module larger parts of the fiber packing will be simulated in future. The spatial and temporal resolution provided by the proposed CFD approach will give close insights in the membrane oxygenation process and will help to improve the design of new artificial lungs. Acknowledgments This research is supported by the Austrian Research Promotion Agency (FFG). Project: MILL – Minimal Invasive Liquid Lung. FFG project number: 857859. References Dzhonova-Atanasova D. B., Tsibranska I. H., Paniovska S. P., 2018, CFD Simulation of Cross-Flow Filtration, Chemical Engineering Transactions, 70, 2041–2046, DOI: 10.3303/CET1870341 Federspiel W. J., Henchir K. A., 2008, Lung, Artificial: Basic Principles And Current Applications, In Encyclopedia of Biomaterials and Biomedical Engineering, Second Edition (Vols. 1–0, pp. 1661–1672). CRC Press, DOI: 10.1081/E-EBBE2-120007349 Galdi G., Robertson A., Rannacher R., Turek S., 2008, Hemodynamical Flows: Modeling, Analysis and Simulation, Birkhäuser, Basel, CHE, DOI: 10.1007/978-3-7643-7806-6 Geers C., Gros G., 2000, Carbon dioxide transport and carbonic anhydrase in blood and muscle, Physiological Reviews, 80(2), 681–715, DOI: 10.1152/physrev.2000.80.2.681 Gonzalez A. B., Rojas-Solorzano L., 2016, Hemodynamic Design Optimization of a Ventricular Cannula Applying Computational Fluid Dynamics (CFD), Chemical Engineering Transactions, 49, 601–606, DOI: 10.3303/CET1649101 Haddadi B., Jordan C., Miltner M., Harasek M., 2018, Membrane modeling using CFD: Combined evaluation of mass transfer and geometrical influences in 1D and 3D, Journal of Membrane Science, 563, DOI: 10.1016/j.memsci.2018.05.040 Hormes M., Borchardt R., Mager I., Rode T. S., Behr M., Steinseifer U., 2011, A validated CFD model to predict O₂ and CO₂ transfer within hollow fiber membrane oxygenators, The International Journal of Artificial Organs, 34(3), 317–325, Jaffer I. H., Reding M. T., Key N. S., Weitz J. I., 2018, Hematologic Problems in the Surgical Patient, In Hematology (pp. 2304-2312.e4). Elsevier, DOI: 10.1016/B978-0-323-35762-3.00159-1 Low K. W. Q., Loon R. V., Rolland S. A., Sienz J., 2017, Formulation of Generalized Mass Transfer Correlations for Blood Oxygenator Design, Journal of Biomechanical Engineering, 139(3), 031007, DOI: 10.1115/1.4035535 Lu J.-G., Zheng Y.-F., Cheng M.-D., 2008, Wetting mechanism in mass transfer process of hydrophobic membrane gas absorption, Journal of Membrane Science, 308(1), 180–190, DOI: 10.1016/j.memsci.2007.09.051 Sutherland W., 1893, The viscosity of gases and molecular force, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 36(223), 507–531, DOI: 10.1080/14786449308620508 Svitek R. G., Federspiel W. J., 2008, A Mathematical Model to Predict CO2 Removal in Hollow Fiber Membrane Oxygenators, Annals of Biomedical Engineering, 36(6), 992–1003, DOI: 10.1007/s10439-008-9482-3 Taskin M. E., Fraser K. H., Zhang T., Griffith B. P., Wu Z. J., 2010, Micro-scale Modeling of Flow and Oxygen Transfer in Hollow Fiber Membrane Bundle, Journal of Membrane Science, 362(1–2), 172–183, DOI: 10.1016/j.memsci.2010.06.034 Windberger U., Bartholovitsch A., Plasenzotti R., Korak K. J., Heinze G., 2003, Whole Blood Viscosity, Plasma Viscosity and Erythrocyte Aggregation in Nine Mammalian Species: Reference Values and Comparison of Data, Experimental Physiology, 88(3), 431–440, DOI: 10.1113/eph8802496 Zhang J., Nolan T. D. C., Zhang T., Griffith B. P., Wu Z. J., 2007, Characterization of membrane blood oxygenation devices using computational fluid dynamics, Journal of Membrane Science, 288(1), 268–279, DOI: 10.1016/j.memsci.2006.11.041 198