One of the primary causes of device failure in membrane oxygenators is the development of thrombi, which is partially attributed to perfusion dynamics and biocompatibility complications.1–4 Over the past couple of decades, many groups have used computational fluid dynamics (CFDs) in evaluating various artificial lung designs.4–14 For membrane oxygenators, numerical analytic assessments of the fluid flow using CFD can be a useful tool in providing information, such as pressure distribution, perfusion dynamics, and gas transport properties. Computational modeling allows rapid design iterations and cost-effective insight into device performance compared with conventional manufacturing and in vitro testing. Furthermore, adjustments in simulated boundary and flow conditions allow optimization of parameters, such as gas exchange efficiency and minimization of recirculation zones. However, before the computed data can be reliably used, CFD results must be experimentally validated because simplifying assumptions must often be used, and the nature of the mathematical formulations may yield multiple solutions.
Experimental validation of CFD results for membrane oxygenators is usually based on correlation of numerical and experimental pressure distributions in these devices.11–14 For example, Funakubo et al.13 compared overall numerical and experimental pressure drop, whereas both Gage et al.12 and Zhang et al.11 sampled pressure distribution from multiple sites drilled along the exterior housing of devices being tested. However, neither of these studies acquired data validating the actual velocity field inside the fiber bundle.
Unfortunately, for membrane devices such as the artificial lung, direct observation of the perfusion dynamics is difficult to accomplish experimentally due to the opacity of the fiber bed, so limited emphasis is placed in confirmation of the velocity field. While only a few studies have experimentally analyzed the perfusion dynamics in membranous devices, most of these data either lack spatial resolution or are restricted to unidirectional flow fields.15–20
Because pressure is an integral property in incompressible Navier–Stokes equation solutions, nonunique velocity fields might exist for a given pressure distribution. Therefore, although it may be shown that calculated and experimental pressure measurements are in good agreement, assertions regarding the fluid flow behavior might be mathematically unjustified and possibly invalid. This article summarizes our efforts to identify a novel means for validating simulated perfusion data through fiber membrane devices such as artificial lungs and to be able to characterize multidimensional velocity fields using an inexpensive, noninvasive, and reliable technique.
This methods section has been divided into four subsections. First, we discuss the experimental setup, addressing the angiography equipment and the bench circuit. Second, techniques for analyzing the acquired angiographic images are described, including key assumptions, digital filters, a correlation-based tracking algorithm, and methods for removing user biases when conducting the various computational tasks. Next, an overview of the computational domain is presented; and finally, settings for the CFD simulations of blood flow in the artificial lung including governing equations, boundary conditions, and modeling assumptions are explored.
Biplane Angiography Setup
A conventional clinical x-ray system was used for acquiring experimental perfusion data due to its availability, high temporal resolution, and relatively inexpensive application. Two-dimensional (2D) biplane digital subtraction angiography (DSA) images were acquired by a Siemens Artis Zee floor-mount system (Siemens, Malvern, PA) at 7.5 frames per second for a contrast-enhanced flow through a commercial membrane oxygenator (Affinity NT Oxygenator, Medtronic, Minneapolis, MN). The oxygenator was oriented such that orthogonal x-ray acquisitions were bisecting the inlet flow channel (Figure 1). The fluid circuit consisted of a 2 L/min flow rate (CentriMag; Levitronix, Zurich, Switzerland) and a 37% aqueous-glycerin by weight solution at 21°C. One hundred fifty milliliters of an iodinated contrast agent iohexol (Omnipaque 300 mgI/ml; GE Healthcare Inc., Waukesha, WI) was diluted to a room temperature viscosity of 3.3 cP and infused through an inline port 1 m upstream of the oxygenator inlet at 33 ml/s by a Mark V Plus infusion pump (MEDRAD Inc., Warrendale, PA). The infusion rate was matched to the circuit flow rate to adequately represent flow through the fiber bundle.21 Biplane images were acquired for 20 s to allow full visualization of the dye entering and leaving the device.
Analysis of Angiographic Images
All image analysis and program coding were performed in MATLAB (v.7.10 R2010a, The MathWorks, Inc., Natick, MA). The contrast agent was assumed homogeneously dispersed, and the presence or absence of dye was represented in a binary fashion.22 Before analyzing measured flow field data, a median filter23 using the nearest 3 × 3 neighboring points was used to provide spatial smoothing to remove noise and preserve edges of the binary data leading to improved handling of the velocity estimates.23,24 In addition, the outer edge of the perfusion projection was assumed to be in the bisecting image plane.
Propagation of the contrast wave through the membrane was tracked using a maximum cross-correlation (MCC) method. Although the MCC method is a nontraditional application of x-ray image analysis, the method has been shown to be a reliable tool for assessing flow progression.25–27 In this approach, a grid is mapped over the image and a template of a given size is defined at each grid point. The computer searches the subsequent time step (viz., x-ray image) for a geometrical pattern with the highest correlation to the template of the previous time step. However, the size of the template and the search region can influence the final output, as false-positive correlations may exist. Therefore, to refine the accuracy of the method, various tests for convergence and filtering algorithms were applied.
To remove the subjectivity in defining template and search window sizes, iterations were run using rectangular template sizes ranging from 5 to 53 pixels per side (pps), and search window sizes ranging from 13 to 149 pps (digital resolution was 5.333 pixels/mm) centered around fixed grid points. For each template size, the average velocity in the fiber bundle calculated by each search window was averaged and is graphed in Figure 2, where the error bars indicate standard deviation between the averages computed by the different search windows. As might be expected from a pattern-based correlation method, larger templates will be more robust against false positives; as represented in Figure 2, increasing template size led to more consistent and constant correlations. The maximum template (53 pps) and search window (149 pps) sizes were used for the analysis in this study, and all search window and template pairs were verified to ensure that the full waveform displacement was captured in subsequent image projections.
Because it is possible that a template pattern may produce a high correlation with a random pattern in the following image, spatially incoherent vectors were filtered using various algorithms. Resultant vectors in the inner flow channel, for example, clearly should not correlate with positions upstream of the flow field and were therefore removed. Similarly, vectors with magnitudes that differed by >2 standard deviations from surrounding vectors were declared inaccurate. Filtering rules such as these were run tangentially with the cross-correlation code such that if an erroneous vector was detected, the next highest correlation that was not in violation of one of the filters was accepted. Finally, only vectors with a cross-correlation >75% were graphed. No interpolation scheme was used, and therefore, sparsities in the vector field will not indicate a region with no perfusion, but rather a nodal point for which no acceptable correlation was found.
Most researchers using CFD to model flow in membrane oxygenators treat blood as a Newtonian fluid undergoing incompressible laminar flow with use of the Ergun approximation for simplifying the computation of the porous fiber membrane.4–7,11–14 It is important to verify that such results are both accurate and useful for design iterations. Geometrical representations of the Affinity membrane oxygenator were drawn in a commercial CAD program (SolidWorks 2007; SolidWorks Corp., Concord, MA) and exported to a grid-generation package (ANSYS Meshing 12.1; ANSYS Inc., Canonsburg, PA). The corresponding problem was then run with commercial CFD software (ANSYS Fluent v12.1; ANSYS Inc.) for approximation of the continuity (mass conservation) and momentum (balance) equations of fluid motion in the absence of any turbulence model. The geometric model consisted of three fluid domains: the inner flow (inlet tube and inner flow channels), an annular fiber bundle, and the outer flow (an outer gap between the fiber bundle and housing and the outlet port); the heater was not included in the computational domain (Figure 3). Due to the symmetrical nature of the design, half of the oxygenator was modeled, bisecting the outlet port, to reduce computational expense.
Governing Equations, Boundary Conditions, and Assumptions
The fiber membrane was treated as a lumped porous medium incorporated into the Navier–Stokes equations as a single momentum sink according to the Ergun equation.11,28 Conservation of mass and momentum balance are expressed as
where u ≡ (u1,u2,u1)T is the velocity vector; D/Dt is the substantial derivative; p is pressure; μ is viscosity; FB is the body force; and S is the source term. For the ith direction, the source term is expressed by the Ergun approximation as
with coefficients α and C2 defined as
where ε is porosity and DO the outer fiber diameter. The superficial velocity, U, is defined as
where Q is the volumetric flow rate; and A is the cross-sectional area of the fiber bundle, independent of porosity. The Ergun approximation prescribes pressure loss through a fiber bundle through simultaneous viscous and kinetic energy losses as factors of individual fiber diameter and fiber bundle porosity. We remark that the expression relating to α, ε, and DO is not unique and is very sensitive to geometric details.
Boundary conditions and modeling parameters for the study are listed in Table 1. Symmetry (slip condition) was assigned to the problem domain boundaries formed from bisecting the oxygenator, and a no-slip condition was assigned to all solid surfaces. Gravitational effects were not included in the model (hence, FB ≡ 0), and a second-order discretization method was used with a convergence criterion of 10–4 for continuity and velocity components using a steady-state solution technique with ANSYS Fluent. Blood was modeled as a Newtonian fluid, with constant viscosity.
The grids used in our investigation were generated using unstructured 3D tetrahedral/hybrid grid cells (Figure 4). Convergence tests of the numerical solutions were performed with six grid sizes ranging from 20K to 1.72M cells, using the grid with the highest resolution as the theoretical exact solution, f. The error between f and grid functions from coarser grids, f*, was quantified using a Euclidian norm composed of the dimensionless pressure drop (p*) and maximum velocity (u*) across the device; this norm was defined as
where f* ≡ (p*,u*)T and n is number of grid cells. Pressure and velocity components were made dimensionless by diving by the dynamic pressure and inlet superficial velocity, Uinlet, respectively, as
The error norms were fitted using a least squares correlation, which characterized the convergence rate of the second-order numerical method as 1.69 (Figure 5). The grid with the highest resolution was selected for all data representation, where the pressure drop and velocity vectors for the finest mesh are estimated to be within an error band 3.3% and 1.6% of the continuum values, respectively.
Progression of the contrast-enhanced flow through the device at 2 L/min is shown in Figure 6. Once the DICOM (Digital Imaging and Communication in Medicine) images had been converted to binary images and spatially smoothed, the cumulative perfusion from previous images was subtracted from each new image, revealing only the net new flow for the given time step (i.e., the current location of the advancing contrast bolus). The flow travels up the inlet tube (frames 1–3), into the tapering channels (frame 5), through the fiber bundle (frames 5–23), into the outer gap (frame 19), and exits through the outlet port (bottom left of image, frames 21–23).
The design of the flow channels along with the packing density of the fibrous media contributes to the mostly radial perfusion dynamics of the fluid flow. For the 2 L/min flow rate, no preferential flow was observed due to the position of the outlet port. However, the progression of the contrast wave front is not fully uniform, as is characterized by the uneven distribution of the waveform in the vertical direction. This image series indicates that flow reaches the corners of the device less readily (see specifically, frames 17–23), including the corner closest to the outlet of the device. At 2 L/min flow, most of the flow traveled through the bundle in approximately 2.5 s (frames 5–19).
Under the given assumptions and modeling parameters, the numerical results correlated well with previously published data for flow rate through the same device at 5 L/min.11 Both pressure distribution and superficial velocity fields were comparable at a test run of 5 L/min (data not shown), and values calculated from the 2 L/min flow rate are proportionally the same, as expected. The computational and experimental (MCC) velocity vectors for corresponding planes are shown in Figure 7. Both images display the physical velocities, where the magnitude in the CFD image is characterized by a color scale and the MCC vector magnitudes according to its vector length. Five representative points were taken to compare the two measurements: two points sampled in the open channels of the device and three points from various sites within the fiber bundle. The vector magnitudes in the MCC figure are the averages of nine vector magnitudes (center plus eight nearest grid points) surrounding the sample location in the corresponding CFD image.
In portions of the device preceding the fiber bundle (inlet tube and inner flow channels), the velocity values calculated for the 2 L/min flow by CFD and MCC methods are consistent. However, in the fiber bundle, where the Ergun approximation was used to characterize the fiber membrane, CFD velocity measurements severely underestimated the experimental velocity.
In our study, the acquisition of high-speed densitometric images were analyzed with a MCC tracking algorithm and showed that fluid velocity through the membrane oxygenator is nonuniform in propagation, which was not indicated in the CFD model. In addition, the experimental data characterized flow speed through the fiber membrane as much higher than previously thought. The higher flow rate through the fiber bundle is significant because this suggests higher shear stresses imposed on the blood constituents, and it is also important in calculating mass transfer characteristics of the device.
When comparing blood flow velocities between the numerical approximations and experimental data, it is important to recognize that the velocities observed experimentally are the physical fluid velocities. Yet, most studies that characterize CFD flow properties in the fiber bundle report the superficial velocity (which does not account for fiber bundle porosity). For example, in a study by Funakubo et al.,13 velocity values were approximated to be 2–6 mm/s for a 3 L/min flow rate through a fiber bundle with 40% porosity. Likewise, for a 5 L/min flow through a 25% porous fiber bundle, Sato et al.14 determined fiber membrane velocities ranging from 7 to 8 mm/s. Similarly, Zhang et al.11 calculated velocities of 5–6 mm/s for a 5 L/min flow through a fiber bundle with 45% porosity. However, without distinguishing between superficial and physical velocities, there exists the potential for misinterpretation of the data. For instance, according to the work by Gartner et al.,4 which correlated thrombus deposition in a fiber membrane with low CFD velocity values, assertions regarding a minimum flow speed “threshold” are significantly different if determined based on the superficial versus the physical, velocity.
The experimental velocities were measured using a cross-correlation method based on binary patterns between sequential image frames. Unfortunately, the MCC technique in its simplest construction does possess a degree of subjectivity with regard to definitions of the template and search window sizes, and it is commonly validated with other forms of measurement such as thermal imaging. However, in these previous studies, supplemental experimental data were required because efforts to identify a nonsubjective means for establishing correlation cutoff values were not identified.26 Alternatively, we believe that through the convergence tests performed on our data set, user biases have likely been reduced.
The major requirement for application of the MCC method was to guarantee that the leading edge of the dye was in the projected image plane. Because densitometric images are integral projections, no information exists for through-plane velocities. However, based on the annular structure of the Affinity Oxygenator and the woven fiber bundle, the shortest flow path, when observing the device’s cross-section, would be direct lateral perfusion (i.e., radial perfusion). This was confirmed by comparing the orthogonal biplane projections, which showed nearly identical densitometric profiles. Because our assumption of visually 2D flow provided no details of through-plane flow, the images were converted to a binary interpretation of the contrasted flow to enhance the geometrical patterning of the MCC method.
Limitations also exist in the experimental data acquisition methods. Advancements in temporal resolution would improve flow field resolution allowing more accurate characterization of the perfusion behavior and minimize errors in tracking algorithms. Image distortions arising from a point-source x-ray beam may necessitate preprocessing modifications to the projection data set, although new flat panel detectors and internal calibrations can help minimize spatial errors. Furthermore, the oxygenator was placed in the center of the x-ray image field to reduce parallax between the biplane projections. Occasionally, alignment of biplane projections did require rotation of the image ± 1–2° about the z-axis. Finally, based on known dimensions of the oxygenator, distortions were measured to be <8%, and therefore, no adjustments were made to reshape the digital x-ray images.
Although many design strategies have been tested over the years, virtually all previous studies confer accurate modeling based solely on pressure distribution because the densely packed fiber bundle precludes conventional methods of experimental flow validation such as particle image velocimetry. Yet, for projection forms of the incompressible Navier–Stokes equations (some form of which is essentially always used in commercial CFD software, including Fluent), the velocity field can be calculated from the momentum equation in the absence of any pressure gradient terms, leaving the divergence-free condition for computing pressure. Furthermore, because pressure is not explicitly stated in the divergence-free equation, and only the gradient of pressure appears in the momentum equation, nonunique solutions may exist for given pressure distributions.29 Although this condition is usually present at higher Reynolds numbers, the various assumptions traditionally used in calculations of blood flow through membrane devices may result in similar errors. In fact, it can be shown that pressure fields computed in this way are least accurate at low Reynolds numbers (Re).29
Further research is required to understand which modeling parameters are specifically responsible for the diminished velocity predictions, but the most likely source is the assumption regarding treatment of the fiber bundle. Woven fiber bundles and higher membrane density result in more homogenous structures, but high blood-side pressures could lead to deformation of the fiber bundle, leading to flow channeling or other phenomena not predictable by CFD simulations. Although directly modeling individual fibers with present-day computing is impractical, the Ergun approximation neglects direct calculation of local perfusion dynamics around the fibers, which could lead to compounding errors in the velocity field.
Recently, Mazaheri and Ahmadi7 simulated fluid flow through a membrane oxygenator with 60% porosity at 4 L/min, for which they compared a porous media assumption versus modeling individual fibers. Although their results indicated nearly a six-fold increase in velocity when using the porous media approximation and was restricted to only a 2D application, they did note that the porous media model did not correctly characterize the nonuniform flow distribution.7 While validity of the porous media approximation is qualitative, at best, there currently exist no practical alternatives for computing flow through an entire 3D device. Although studies that analyze flow around individual fibers are useful for assessing optimal packing arrangements based on local mass transfer, fiber resistance, and shear stresses,30–32 these results do not translate into global perfusion dynamics such as geometrically induced eddies or regional mass flux.
Summary and Conclusion
Previously, direct observation of fluid flow properties was difficult to achieve experimentally due to the opacity of the fiber bed. Although computational simulations are useful for providing insight into blood flow characteristics in membrane oxygenators, assertions regarding internal flow behavior might be mathematically unjustified when based exclusively on correlation of numerical and experimental pressure distributions.
To validate the CFD velocity field, a novel approach using clinical x-ray DSA and a pattern-searching tracking algorithm provided quantitative measurements of flow perfusion inside a membrane artificial lung. Using conventional methods for modeling blood flow through the artificial lung, we were able to perform direct analysis of numerical and experimental velocity fields. Based on the results of this research, it was shown that when using the Ergun equation for approximating momentum loss in the fiber bundle, velocity predictions were significantly lower than corresponding experimental measurements. These data suggest the need for experimental validation of flow characteristics, which may play a significant role in mass transfer calculations in fiber bundles, or in design of the housing geometry.
The authors thank Scott Kratzwald, Brett Wand, and Sean Brown for their assistance with the x-ray acquisitions.
1. Eash HJ, Jones HM, Hattler BG, Federspiel WJ. Evaluation of plasma resistant hollow fiber membranes for artificial lungs. ASAIO J. 2004;50:491–497
2. Snyder TA, Eash HJ, Litwak KN, et al. Blood biocompatibility assessment of an intravenous gas exchange device. Artif Organs. 2006;30:657–664
3. Lehle K, Alois P, Otto G, et al. Efficiency in extracorporeal membrane oxygenation—Cellular deposits on polymethylpentene membranes increases resistance to blood flow and reduce gas exchange capacity. ASAIO J. 2008;54:612–617
4. Gartner MJ, Wilhelm CR, Gage KL, Fabrizio MC, Wagner WR. Modeling flow effects on thrombotic deposition in a membrane oxygenator. Artif Organs. 2000;24:29–36
5. Goodin MS, Thor EJ, Haworth WS. Use of computational fluid dynamics in the design of the Avecor Affinity oxygenator. Perfusion. 1994;9:217–222
6. Asakawa Y, Funakubo A, Fukunaga K, et al. Development of an implantable oxygenator with cross-flow pump. ASAIO J. 2006;52:291–295
7. Mazaheri AR, Ahmadi G. Uniformity of the fluid flow velocities within hollow fiber membranes of blood oxygenation devices. Artif Organs. 2006;30:10–15
8. Verdonck P. The role of computational fluid dynamics for artificial organ design. Artif Organs. 2002;26:569–570
9. Graefe R, Borchardt R, Arens J, Schlanstein P, Schmitz-Rode T, Steinseifer U. Improving oxygenator performance using computational simulation and flow field-based parameters. Artif Organs. 2010;34:930–936
10. Fill B, Gartner M, Johnson G, Horner M, Ma J. Computational fluid flow and mass transfer of a functionally integrated pediatric pump-oxygenator configuration. ASAIO J. 2008;54:214–219
11. Zhang J, Nolan T, Zhang T, Griffith B, Zhongjun J. Characterization of membrane blood oxygenation devices using computational fluid dynamics. J Membr Sci. 2007;288:268–279
12. Gage KL, Gartner MJ, Burgreen GW, Wagner WR. Predicting membrane oxygenator pressure drop using computational fluid dynamics. Artif Organs. 2002;26:600–607
13. Funakubo A, Taga I, McGillicuddy JW, Fukui Y, Hirschl RB, Bartlett RH. Flow vectorial analysis in an artificial implantable lung. ASAIO J. 2003;49:383–387
14. Sato H, Taga I, Kinoshita T, Funakubo A, Ichiba S, Shimizu N. In vitro
evaluation of a newly developed implantable artificial lung. Acta Med Okayama. 2006;60:113–119
15. Chen V, Li H, Fane A. Non-invasive observation of synthetic membrane processes—A review of methods. J Membr Sci. 2004;241:23–44
16. Heese F, Robson P, Hall L. Quantification of fluid flow through a clinical blood filter and kidney dialyzer using magnetic resonance imaging. IEEE Sens J. 2005;5:273–276
17. Hirano A, Yamamoto K, Matsuda M, et al. Flow uniformity in oxygenators with different outlet port design. ASAIO J. 2009;55:209–212
18. Lu J, Lu W. Blood flow velocity and ultra-filtration velocity measured by CT imaging system inside a densely bundled hollow fiber dialyzer. Int J Heat Mass Transfer. 2010;53:1844–1850
19. Poh CK, Hardy PA, Liao Z, Huang Z, Clark WR, Gao D. Effect of flow baffles on the dialysate flow distribution of hollow-fiber hemodialyzers: A nonintrusive experimental study using MRI. J Biomech Eng. 2003;125:481–489
20. Ronco C, Brendolan A, Crepaldi C, Rodighiero M, Scabardi M. Blood and dialysate flow distributions in hollow-fiber hemodialyzers analyzed by computerized helical scanning technique. J Am Soc Nephrol. 2002;13(suppl 1):S53–S61
21. Dorsaz PA, Doriot PA, Dorsaz L, Chatelain P, Rutishauser W. A new densitometric approach to the assessment of mean coronary flow. Invest Radiol. 1997;32:198–204
22. Prause GM, Onnasch DW. Binary reconstruction of the heart chambers from biplane angiographic image sequences. IEEE Trans Med Imaging. 1996;15:532–546
23. Pratt WB, Krishna P, Olsen LJ. Hsp90-binding immunophilins in plants: The protein movers. Trends Plant Sci. 2001;6:54–58
24. Condell J, Scotney B, Morrow P. The effect of presmoothing image sequences on the computation of optical flow. Lect Notes Comput Sci. 2006;4141:780–791
25. Dransfeld S, Larnicol G, Le Traon PY. The potential of the maximum cross correlation technique to estimate surface currents from thermal AVHRR global area coverage data. IEEE Geosci Remote Sens Lett. 2006;3:508–511
26. Bowen M, Emery W, Wilkin J, Tildesley P, Barton I, Knewtson R. Extracting multiyear surface currents from sequential thermal imagery using the maximum cross-correlation technique. J Atmos Oceanic Technol. 2002;19:1665–1676
27. Crocker R, Matthews D, Emery W, Baldwin D. Computing costal ocean surface currents from infrared and ocean color satellite imagery. IEEE Trans Geosci Remote Sens. 2007;45:435–447
28. Ergun S. Fluid flow through packed columns. Chem Eng Prog. 1952;48:89–94
30. Dierickx PW, De Wachter D, Verdonck PR. Blood flow around hollow fibers. Int J Artif Organs. 2000;23:610–617
31. Chan KY, Fujioka H, Bartlett RH, Hirschl RB, Grotberg JB. Pulsatile flow and mass transport over an array of cylinders: Gas transfer in a cardiac-driven artificial lung. J Biomech Eng. 2006;128:85–96
32. Zierenberg JR, Fujioka H, Cook KE, Grotberg JB. Pulsatile flow and oxygen transport past cylindrical fiber arrays for an artificial lung: Computational and experimental studies. J Biomech Eng. 2008;130:031019