Three-dimensional computational model of a blood oxygenator reconstructed from micro-CT scans

Cardiopulmonary bypass procedures are one of the most common operations and blood oxygenators are the centre piece for the heart-lung machines. Blood oxygenators have been tested as entire devices but intricate details on the ﬂow ﬁeld inside the oxygenators remain unknown. In this study, a novel method is presented to analyse the ﬂow ﬁeld inside oxygenators based on micro Computed Tomography ( μ CT) scans. Two Hollow Fibre Membrane (HFM) oxygenator prototypes were scanned and three-dimensional full scale models that capture the device-speciﬁc ﬁbre distributions are set up for computational ﬂuid dynamics analysis. The blood ﬂow through the oxygenator is modelled as a non-Newtonian ﬂuid. The results were compared against the ﬂow solution through an ideal ﬁbre distribution and show the importance of a uniform distribution of ﬁbres and that the oxygenators analysed are not susceptible to ﬂow directionality as mass ﬂow versus area remain the same. However the pressure drop across the oxygenator is dependent on ﬂow rate and direction. By comparing residence time of blood against the time frame to fully saturate blood with oxygen we highlight the potential of this method as design optimisation tool. In conclusion, image-based reconstruction is found to be a feasible route to assess oxygenator perfor- mance through ﬂow modelling. It offers the possibility to review a product as manufactured rather than as designed, which is a valuable insight as a precursor to the approval processes. Finally, the ﬂow analysis presented may be extended, at computational cost, to include species transport in further studies.


Introduction
According to the National Health Service (NHS [19] ) and Health and Safety Executive (HSE [10] ), up to 25,0 0 0 people die of Chronic Obstructive Pulmonary Disease (COPD) every year and the disorder affects over a million individuals in Great Britain. In addition up to 400 people develop Acute Respiratory Distress Syndrome (ARDS) with a mortality rate of over 50%. Improving treatment techniques have a direct social and individual impact in saving lives, increasing life expectancy and reducing cost to the public health services. Two different methods are clinically approved and are used to transfer oxygen into the blood: • Mechanical ventilation, mainly controlled by pressure and/or volume (Chatburn et al. [4] ), transfers oxygen into the lungs of the patient for the gas exchange. Often the ventilation pressure and oxygen concentration are set very high to overcome * Corresponding author.
the impaired lung function. A potentially resulting barotrauma, volutrauma and oxygen toxicity may prevent or slow down the lung recovery. • Extra Corporeal Membrane Oxygenation (ECMO) is bypassing the cardiopulmonary cycle by oxygenating the blood outside the patient by pumping blood through a bed of micro-porous Hollow Fibres Membranes (HFM). For details on the development of ECMO, the reader is referred to Haworth [8] .
For the last decade ECMO devices based on HFM prevailed and has brought a new branch of studies, where the oxygenator itself is investigated. The oxygenator is one of the centre pieces of the ECMO circuit and the development of hollow fibres with membranes or micro porous walls is advancing rapidly with novel materials and coatings. The effort to optimize the HFM-assembly concentrates to reduce the pressure drop across the device, to minimize haemolysis, decrease the priming volume of blood, and increase biocompatibility and lifetime of HFM oxygenators.
Computational fluid dynamic (CFD) models are used to study blood flow path, heat exchange, pressure drop, stress analysis, mass transfer of oxygen and carbon dioxide, different convectiondiffusion models, blood stagnation, thrombogenicity, etc. These models use either a homogeneous porous media to describe the fibre bundle in bulk or a heterogeneous approach to model single fibre arrangements. The rheology is usually implemented as Newtonian flow field or more realistically described as shear thinning fluid. Non-Newtonian shear thinning numerical models for blood are well known and described by Johnston et al. [11] or Marcinkowska-Gapi ńska et al. [15] .
For example Gartner et al. [6] and Pelosi et al. [20] use a porous media approach to model thrombogenic depositions and study the thrombogenic potential of oxygenators. Zhang et al. [25] predict blood-gas exchange and pressure drop with a Newtonian flow field through a porous media. A comparison of different porous media models is presented by Khanafer et al. [13] with a set porosity ( = 0 . 75 ). Hormes et al. [9] present a novel diffusivity model to predict the O 2 and CO 2 mass transfer and partial pressure, which is validated by comparing the numerical solution with a purpose built oxygenator. The oxygenator and the CFD model use a uniform distribution of fibres.
Commonly the same partial pressure P O 2 is applied at the fibre surface for convection and diffusivity. Taskin et al. [23] presented a novel model to describe gas exchange as profile on the fibre surface. The study is comparing a single and multi fibre approach with uniform fibre distribution.
A study by Nagase et al. [18] presented mass transfer correlations based on the theory for heat transfer across a tube bank and concluded that the mass transfer performance of membrane oxygenators is attributable to the hollow fibre arrangement. The layouts considered are uniformly arranged staggered or squared (inline) fibres.
An oxygenator cross-section modelled with individual fibres in two dimensions is compared to the porous media approach by Mazaheri and Ahmadi [16] . To avoid computational costly simulations of individual fibres, novel porous models are developed to integrate the fibre distribution for a better local process mapping (Low et al. [14] , Zhang et al. [25] ) and an attempt to "adjust" the porosity based on fibre orientation is described by Bhavsar et al. [3] .
But all off these models have been developed with the assumption of a uniform distribution of fibres either staggered, squared or crossed (3D) or simplified with a porous media approach. Although a good agreement between numerical predictions and experimental results has been reported in the literature by Consolo et al. [5] and Pelosi et al. [20] , the porous media approach is inherently unable to capture and characterize the intricate details on the flow field within the fibre bundles.
Imaging modalities can be used as a diagnostic tool to investigate the micro-structural components in an oxygenator and, when combined with CFD, can even provide information on the localised functional behaviour (flow/oxygenation) of the device. In this paper we demonstrate a method for building a full-scale, image-based, three-dimensional computational model of a blood oxygenator. All individual fibres are reconstructed from the images of a μCT scan and CFD is performed using a Non-Newtonian fluid model to investigate the local flow field inside the device.

Prototypes and methods
The common method of extra-corporeal oxygenating blood is achieved by pumping ambient air or through the core of HFM's whereas the blood is flowing on the outside of the fibres. Molecular diffusion increases the oxygen level in the red blood cells and removes carbon dioxide from the venous blood.

Prototypes
Two manually assembled oxygenators are used for experimental evaluation. The blood samples, outdated packed cells, are sourced from the Welsh Blood Service [24] which is accredited by the Medicines and Healthcare products Regulatory Agency [7] . The experiments are conducted by Haemair Ltd. and follow their standard of operation procedure SOP2007v04. One of the prototypes is depicted in Fig. 1 . The core of the oxygenator is built with folded mats of hollow fibres produced by Membrana GmbH [17] . The fibres are state of the art ('OXYPLUS ® ') and have an outer diameter of approximately 380 μm and a nominal wall thickness of 90 μm. The fibre wall is highly porous ( ≥ 55%) to allow molecular diffusion but blocking any fluid. For more details see Membrana GmbH [17] . Resin seals both ends of the oxygenator around the fibres, forcing the blood to flow between inlet and outlet ( Fig. 1 ). Ambient air is pumped through the hollow core of the fibres. The nominal outer dimensions of the enclosure is 10 × 20 × 100 mm. Both prototypes are built to support flow along the fibres and to prevent crossflow between the fibres at the same time, to minimize clotting and haemolysis. Additionally a virtual idealised geometry is used to set up the CFD model. Size, volume and fibre count are comparable to the real world devices. The fibres in the virtual device are uniformly arranged to exploit symmetry and simplify meshing. The virtual model is used for comparison and to test the numerical solution. For all devices the key characteristics are listed in Table 1 .

Reconstructed geometry
μCT scans were performed on both prototypes using a Nikon XT H 225 with a voxel size of 22 μm. The source voltage and current were set to 55 kV and 174 μA, respectively. The exposure time for each radiograph was 2 s, with 720 radiographs being collected over 360 °.
A complete set of reconstructed two-dimensional images was computed with ∼1500 × 10 0 0 pixels and a displacement in zdirection of ∼30 μm. The images captured from the μCT scan show the fibres and their arrangement in great detail (see Fig. 2 ). The warp thread in OXYPLUS ® mats is a polyethylene terephthalate (PET) multifilament yarn (33 dtex 1 ) with 24 filaments. The threads, with an estimated diameter of ∼70 μm ± 10, are woven into the fibre bundle every 10 mm to keep the fibres at a minimum distance and hence allow blood flow in-between. Due to the small size of the threads, we assumed that any pressure drop or mixing affects can be neglected for this study. To verify that the fibres remain in a straight line (in z -direction, see Fig. 4 ) throughout the device a visual inspection was conducted on a 3D compilation [22] of the μCT images (see Fig. 2 ) as well as overlaying 2 images from the opposite ends. A cross-section image for each oxygenator ( Fig. 3 a and  b) was used to recreate a full size real world three-dimensional geometry to run CFD-simulations. Both prototypes were manually Static priming volume (mL) Blood volume oxygenated (mL) 17.95 12.10 13.53 assembled by folding HFM mats into the plastic enclosure. The device SRD3083 has added spacers (see Fig. 4 ) between the folds to support a uniform distribution of fibres throughout the device. The spacers are 100 μm thick and are completely embedded in the resin sealed ends.

Modelling technology and mesh dependency
We use the cell centre based finite volume formulation to solve the incompressible Navier-Stokes equation for momentum and mass conservation. The blood flow through the devices is modelled with ANSYS Inc. [1] -FLUENT. The solver was set for a pressure-based steady state solution with absolute velocity. The problem was discretised using a least-square cell based gradient formulation combined to a second order upwind discretisation scheme for momentum and a second order scheme for pressure. The SIMPLE algorithm was used for pressure coupling. A representative section of fibres from the device SRD3078 was used to build a smaller version of blood oxygenator with a manually refined mesh with approximately 2 million elements. A grid convergence study was conducted to find a suitable set of parameters to create an appropriate unstructured mesh for all devices. With a minimum size restriction for elements set to 1 × 10 −6 m and a biased sweep mesh for the main body, we have a confidence level of ±2% error in regards to pressure drop for the numerical solution. To limit the amount of elements further, a minimum gap size of ࣛ 10 μm is imposed between individual fibres or between fibres and enclosure. This is achieved by merging fibres or joining fibres with the enclosure wall if the gap is smaller (see Fig. 5 ). A total of 7 fibres are merged for the device SRD3078 and 8 fibres for the device SRD3083. For all meshes a minimum of five elements is enforced between the fibres. All three geometries have a final unstructured mesh with ∼20 million cells. The inlet and outlet ( φ 3 mm) tube length is set to 20 mm to support a fully developed fluid flow entering the oxygenator. The outlet is configured as pressure outlet and is set to 0 Pa and the fibres and the enclosure are defined as no-slip walls.

Simulation and visualisation of blood flow
A study comparing three rheological models (Casson, Ree-Eyring and Quemada [15] ) concluded that the Quemada model fitted most precisely with their experimental findings. To solve the blood flow as incompressible Non-Newtonian fluid, the model (Quemada [21] , Eq. (1) ) is implemented as a user defined function Fig. 3. (a,b) μCT images of a cross-section for two manually assembled blood oxygenator prototypes. (b) Spacers between the folds help to create a more uniform distribution of fibres. (c) Idealised CFD-geometry to provide a comparable device with uniformly distributed fibres. (D) "Spacers" inserted on both ends for the device SRD3083 to support a uniform distribution of fibres (cf. Fig. 3 a and b).  [15] and listed in Table 2 . Shear strain rate, denoted as ˙ γ is obtained from the flow field and the density is set to 1050 kg m −3 for average whole blood. With the above, viscosity can be described as The boundary condition for the blood flow is set as mass flow rate to 8 . 75 × 10 −5 kg s −1 ( 5 mL min −1 ) and 6 . 65 × 10 −3 kg s  same boundary conditions and parameters applied. The idealised geometry is symmetrical, hence there is no need to calculate the reverse flow for comparison. To visualise mass flow versus area ( Fig. 7 ) an ANSYS plug-in 'MassFlowAreaRatio' was developed to interrogate the results interactively. The plug-in is using the built-in API to Perl and the Power Syntax from CFX-Post (ANSYS CFX Reference Guide release 14.5) and allows to set arbitrary numbers of cross-sections in any direction with a pre-defined threshold for mass flow. Each step of the integration may be saved and stacks of images can be automatically saved to create an animation (see supplementary videos online-10.1016/j.medengphy.2017.06.035 ).

Oxygenator shunt fraction
Perfusion mismatch or impaired lung function is usually calculated using the Berggren Equation for shunt fraction ( Q s / Q t [2] ) to express the ratio of blood bypassing oxygenation ( Q s ) to the total cardiac output ( Q t ). In this study we use a similar approach to de- The blood being oxygenated is defined as V b * and is derived in Eq. (3) and rearranged in Eq. (4) . In essence it describes the volume of annuli with thickness d around the fibres. See Table 1 for parameters length ( l ), count of fibres ( N f ) and priming volume V b .
Assuming that only blood within a distance of d = 10 μm to the outer fibre wall is active in gas exchange (approximately the size of human capillaries), Eq. (2) is used to find the ratio of blood not being exposed to the hollow fibre membrane and hence bypassing the oxygenation. O sf is only valid for non-turbulent flow regimes with a perfect flow in longitudinal direction ( z ).

Cross flow ratio coefficient
The presented oxygenator prototypes have been built to support flow along the fibres to minimize clotting and haemolysis as well as to study flow rate and residence time. The goal is a perfect flow in longitudinal direction z (parallel to the fibres) in order to reduce priming volume and minimise residence time with complete oxygen saturation of the blood. Due to an asymmetrical arrangement of inlet and outlet, as well as a deviation from a uniform fibre distribution in the prototypes, we propose to use the coefficient O cf ( Eq. (5) ) to quantify the amount of cross-flow present.

Results
As shown in Figs. 2 and 3 , the images reconstructed from the μCT scan represent the oxygenator and the individual fibres in great detail. This allows to rebuild an accurate three-dimensional geometry with a CAD software which subsequently can be used to set up a CFD model. A visual inspection of two cross-section images from opposite ends and of slices of the 3D reconstruction shows very little deviation of the fibres in longitudinal direction ( z ). The spacers built into SRD3083 increase the uniformity of the fibre distribution significantly. A two-dimensional cross-section plane is defined at position l /2 and for each cell in the mesh, mass flow and area are extracted and plotted (see Fig. 6 ). A formal description is found in Eq. (6) . If Q bc is the prescribed mass flow rate through the device and A CS is the cross-sectional area of the oxygenator then an area A can be defined for every value of the flow fraction Q , where the flow fraction is defined as 0 ≤ Q ≤ 1.
The full data set for the cross-section is shown in Fig. 6 . Note that the graph is a reflection on how uniformly the flow is distributed across the cross-section of the oxygenator. This is illustrated by a snap-shot of mass flow to area ratio at 70% in Fig. 7 . The simulations show that large flow channels absorb most of the flow. For the device SRD3078 70% mass flow goes through only 10% area of the cross-section whereas the uniform arrangement of fibres in the idealised geometry yields a much more even mass flow versus area ratio. It can be further seen that the distribution of the fibres within the device has little impact on the flow directionality as mass flow versus area remain the same. Figs. 6 and 7 highlight the impact of a non-uniform distribution of fibres within the oxygenator.
The pressure drop however is flow rate and direction dependent. Eight experiments on three different builds of this prototype model are conducted. With an average flow rate of 380 mL min −1 a mean pressure drop of 3792 Pa is measured. The computational model for SRD3078 shows good agreement with 4053 Pa and 3853 Pa [R]everse. As these devices have all been manually assembled, they have potentially large flow channels as shown in Fig. 7 . SRD3083 (with spacers) has a calculated pressure drop of 6443 Pa and 6218 Pa [R] and the idealised geometry 10,186 Pa. The pressure drop for 5 mL min −1 varies between 40 Pa and 220 Pa, but there is no experimental data available to compare with.
Cross-flow coefficients ( O cf ) for all devices and flow rates are depicted in Fig. 8 . Results are extracted at 14 cross-sections throughout the main body of the oxygenators at locations corresponding to the mesh elements. Increased cross-flow activity is visible close to inlet and outlet and is dependent on flow rate and fibre distribution. Due to a grid refinement at the inlet/outlet Eq. (5) is bound to overestimate O cf and needs to be accounted for. The data suggest that decreasing the uniformity of fibres (SRD3078) leads to flow direction dependency whereas a uniform fibre distribution has the same O cf in both directions.
With Eq. (2) and setting d = 10 μm (representative of the size of capillaries in the human lung) we can calculate the shunt fraction O sf for each device (see Table 3 ). The physical layout of the devices support flow parallel to the fibres, but as shown in Fig. 8 crossflow is present. A parameter to adjust O sf should be imposed to account for the cross-flow; otherwise Eq. (2) is overestimating the dead space. We propose to use O cf as a relaxation factor for the shunt fraction O sf ( Eq. (7) ).  Table 3 Shunt fractions independent of ( O sf ) and dependent on ( O s f * ) q (flow rate).  Table 4 Streamline ( S l ) in seconds and standard residence ( τ ) time in seconds for flow rate q .
Thus we can account for a device dependent flow behaviour occurring from different fibre arrangements as well as for different flow rates. Table 3 lists the improved shunt fractions. The lower limit of O s f * is V b whereas the upper limit is defined by V b * . The blood volume around the fibres being oxygenated, V b * , is adjusted by the distance d (see Eq. (2) ) and needs to be chosen such that Residence time of blood within the oxygenator was calculated in two ways; (i) by extracting streamlines ( S l ) from the numerical solution and (ii) by the generic equation for residence time τ = V/q with V = V b and q = Q d , set to the corresponding flow rate.
The standard residence time is valid for steady laminar flow only. Reynolds number extracted for the numeric models are R e < 1 for 5 mL min −1 flow rate and R e < 80 for 380 mL min −1 hence flow can be assumed laminar. Streamlines were extracted in the main body of the oxygenator starting from 10 0 0 equally spaced seed points. For the results in Table 4 only the streamlines reaching the outlet ( n > 95 % for all data sets) are considered. Fig. 9 a and b depicts the streamline time distribution for two different flow rates. For better visualisation, only data points ≤ S l med (median) are plotted. The bin size for 5 mL min −1 flow rate is 1 s and for 380 mL min −1 0.01 s. Standard residence time in seconds (s) for 5 mL min −1 is calculated with τ = (V b × 60) / 5 and for 380 mL min −1 with τ = (V b × 60) / 380 . The results are listed in Table 4 .

Discussion
The presented μCT scans provide a feasible way to reconstruct an HFM-oxygenator with computer-aided design software (CAD) based on individual hollow fibres. The geometry can subsequently be used to run CFD models to study relevant design parameters for blood oxygenators. As an example our model shows the insensitivity of flow direction in regards to mass flow versus area but highlights the importance of a uniform fibre distribution. In the case of SRD3078 70% of blood is flowing through the device with very little contact to fibres, suggesting limited oxygenation efficiency ( Fig. 7 ). By comparing pressure drops and residence times we can highlight the dependency on flow direction.
Especially residence time of blood is an important factor in the design of a blood oxygenator, to provide optimal mass transfer and blood oxygen saturation. We have shown that residence time based on velocity vectors (streamlines) in a full scale model reveal a distribution which cannot be represented by the standard equation for fluid resident time in a reservoir ( Table 4 and Fig. 9 ). According to Kang et al. [12] the time of oxygen saturation of red blood cells is between 60 ms (rest) and 120 ms (exercise). This time scale is at least a magnitude shorter than the minimum residence time of blood in the studied oxygenator. Hence our results suggest that even with a significantly reduced length of the main body of the oxygenator sufficient oxygenation and a potential decrease of the device's thrombogenicity could be achieved.
Furthermore we have presented a method to quantitatively assess cross-flow activity and a device dependent shunt fraction. Despite the prototypes physically supporting longitudinal flow, we could reveal cross-flow activity probably due to the positioning of the blood inlet and outlet tube and a non-uniform distribution of fibres. Oxygenator shunt fraction, O s f * , takes flow rate, direction and cross-flow activity into account and can be adjusted using a single parameter d to describe the blood volume being oxygenated. We have used a nominal value of d = 10 μm to reflect human capillaries, but this may be different for other hollow fibre membrane oxygenators and needs further investigation. Combined they show the potential as a qualitative design optimization parameter.
So far, to the best of our knowledge, the existing literature reports studies with a porous zone or a uniformly distributed arrangement of fibres to model the blood flow through the oxygenator. Using either method assumes a homogeneous uniform permeability or a pre-calculated permeability based on a uniform distribution of fibres. In this study, we have shown that such a distribution of fibres is unlikely and even a small variability in gap size between the fibres changes the flow behaviour significantly.

Conclusions
μCT scans provide a feasible tool to rebuild the geometry of a hollow fibre membrane oxygenator. We have shown how to quan-titatively and qualitatively assess prototype oxygenators and can conclude that a significant change in transport phenomena and mass exchange is to be expected if large gaps between fibres are present, i.e. uniform distribution of fibres are paramount for blood oxygenators.
To solve a full scale meshed oxygenator is computationally more expensive than a porous media approach, but yields an insight into the oxygenator otherwise not achievable. We were able to solve the flow problem on the model presented ( ∼20 million cells) within 24 h on a single workstation computer (Xeon E3-1240 v3, 4 cores at 3.4 GHz, 16 GB RAM).
We have shown the capability of using a full scale threedimensional numerical model to realistically simulate flow behaviour using accurate geometrical descriptions of the oxygenator rather than idealised geometries. As the boundaries of the fibres are accurately meshed, further analysis of mass transfer is possible by solving a convection-diffusion problem based on the flow solutions in this work. The reconstruction of individual fibres for example, could be used to apply the transport models presented by Taskin et al. [23] to calculate individual P O 2 profiles and study the mass exchange in detail. Such a model is therefore a feasible tool to study design optimisations of HFM oxygenators and may facilitate the development of a wearable ECMO device and/or an implantable artificial lung.

Ethical approval
Not required.