All methods of extracorporeal respiratory support, including any form of extracorporeal membrane oxygenation (ECMO), will require vascular access to drain deoxygenated venous blood and subsequently reinfuse oxygenated blood back to the patient. Moreover, an ambulatory, long-term paracorporeal artificial lung system for end-stage chronic lung diseases such as chronic obstructive pulmonary disease (COPD) would require a dependable, less traumatic vascular access cannula. But traditionally, venous access requires multiple sites of cannulation with a long and complicated circuit. Our patented Wang-Zwische double-lumen cannula (DLC) was developed for long-term respiratory support and allows both venovenous ECMO and an ambulatory artificial lung to be applied by a single-site, percutaneous jugular vein cannulation (Figure 1).1 Previous pediatric DLCs were prone to kinking, and design characteristics resulted in significant recirculation which led to deterioration in gas exchange performance, but the Wang-Zwische DLC has been characterized by high-efficiency gas exchange with minimal recirculation.1 Although the DLC allows for shortening of the extracorporeal circuit, it is also the portion of the circuit with the highest resistance, shear stress, and consequently, potential for blood trauma. To investigate the potential severity of blood trauma contributed by the cannula, blood flow characteristics in the Wang-Zwische DLC have been studied.
Bench testing and in vivo large animal tests are often used to evaluate blood trauma grossly through measurements of hemoglobin molecules from ruptured red blood cells. However, computational fluid dynamics (CFD) can provide more detail regarding shear stress distribution by identifying areas with the highest values, thereby suggesting regions to address in future product improvement. Other researchers have applied CFD toward calculating steady and unsteady shear-induced hemolysis through a finite element technique and further used CFD in predicting hemolysis associated with various cardiovascular devices by incorporating factors contributing to blood cell damage, such as shear stress and exposure time.2–7 However, because of the nature of the governing equations, only experimentally validated CFD can concretely predict the reality of the blood flow pattern. This work aims to analyze the flow field inside the 27-Fr DLC to evaluate its potential for causing blood trauma and to provide insight for future device improvement. To this purpose, three-dimensional (3D) numerical simulations were performed, with corresponding pressure and mass flow rate calculations measured on a steady-flow bench setup.
Materials and Methods
The 27-Fr Wang-Zwische DLC (marketed as the Avalon Elite, Avalon Laboratories LLC, Rancho Dominguez, CA) was chosen for the flow dynamic study. A CentriMag pump (Levitronix LLC, Waltham, MA) was connected to the DLC by 3/8 × 3/32 Tygon tubing. The DLC was placed in an open reservoir. The circuit and reservoir were primed with a 50% aqueous glycerol solution and maintained at 42°C by a VWR heater (PolyScience, Niles, IL) to mimic human blood viscosity.8 The CentriMag pump was set to a flow rate of 2 L/min, which is the flow used in our developmental ambulatory artificial lung experiments, and is directly applicable to most patients with end-stage chronic lung diseases. Furthermore, an elevated temperature was tested as this scenario may be used in some patients in extreme situations. A Transonic flow sensor (HT110 flowmeter, Transonic System Inc., Ithaca, NY) acquired flow rate measurements while two TruWave pressure transducers (Edwards Lifesience, Irvine, CA) were applied to monitor inlet and outlet pressure of the DLC. An Ultra-Miniature Millar Mikro-Tip catheter (Millar Instruments Inc., Houston, TX) was placed from the end of the DLC into the opening of the drainage and infusion lumen and progressively advanced to both the inlet and outlet of the DLC to trace the changes in pressure along the length of each lumen. A PCU-2000 pressure control unit (Millar Instruments Inc.) was used to transfer the pressure signal from the Millar pressure catheter to the data acquisition system at a rate of 1000 samples/s through a NI BNC-2110 module linked to an NI DAQCard-AI-16XE-50 data acquisition card (National Instruments, Austin, TX). All recording and data compilation were performed in MATLAB (The MathWorks Inc., Natick, MA).
The numerical analysis of the flow through the catheter was carried out separately for each lumen. Numerical computations were performed using the steady three-dimensional Reynolds-averaged Navier-Stokes (RANS) equations and the low-Reynolds k-ω turbulence model.9 Discretization of the governing equations was based on a cell-centered finite volume method (FLUENT v.6.2.16, 2005, Fluent Inc., Lebanon, NH); while a second-order accuracy in space for the diffusive and convective terms was used. In particular, the convective fluxes were computed adopting an upwind interpolation, while the diffusive fluxes were discretized using a centered-differencing scheme. Gradients of the solution reconstruction and the velocity components for the evaluation of the viscous fluxes have been computed using the Node-based Green-Gauss method. The discrete system of RANS equations was solved using a steady segregated implicit solver (SIMPLE).
Grid System and Boundary Conditions
The computational domain coincides with the real geometry of the DLC. In addition, two intake plenums were added to simulate the intake through the side holes as well as through the inlet section of the drainage lumen. Furthermore, to simulate the flow expansion at the exit of the infusion lumen, an exhaust plenum was added around the exit section.
Two types of body-fitted unstructured grids were used for DLC discretization: a nonconformal hexahedral grid for the drainage lumen (HEXPRESS v.2.8-4, Numeca International, Brussels, Belgium) and a conformal hexahedral grid for the infusion lumen (GAMBIT v.2.2.30, 2005, Fluent Inc.). In both meshes, the near wall grid was fine enough to resolve the laminar sublayer (y+ = 1/2). To obtain a mesh-independent solution, three different grid sizes were used for each lumen; norms were defined as maximum absolute differences in wall shear stress (WSS) and velocity magnitude at multiple sites of interest and compared for respective grid refinements. A 6% difference was calculated between the coarsest grid and the first grid refinement, and less than a 1% difference was determined between the first refined grid and the third grid with the highest number of elements. The grid with the highest refinement was used for all calculations and figures and consisted of 4.2 million elements for the drainage lumen and 1.05 million elements for the infusion lumen.
Blood was assumed to be a homogeneous, incompressible Newtonian fluid with a density of 1050 kg/m3 and a dynamic viscosity of 0.0035 kg/m·s. Computations on both lumens were carried out by considering a fixed flow rate equal to the one used during the test rig experiments (2 L/min) and by imposing very low values of turbulence intensity. In particular, as the flow field is laminar in the region surrounding the inlet sections of the two intake plenums of the drainage lumen, turbulence values were set to zero.
The open-air reservoir static pressure was imposed at the exit of the infusion lumen as well as at the entrance of the drainage lumen. By means of an iterative procedure, a pressure that resulted in a flow rate equal to the one measured in the experimental rig was selected at the outlet section of the drainage lumen. Furthermore, the walls were imposed with a no-slip condition. Table 1 shows flow and turbulent values imposed at the boundaries of each lumen. These values were chosen by considering the flow regimen at the inlet boundaries of each lumen and turbulent length scale relationships.
Numerical and Experimental Comparison
The comparison between experimental and numerical static pressure along the axis of each lumen was carried out to assess the ability of the numerical model to simulate the 3D flow field inside the DLC. Experimental measurements and simulated results are compared in Figure 2, A and B for the drainage and the infusion lumen, respectively. The data characterize pressure trends along the length of the cannula and segment the cannula into various regions, each revealing a unique pressure gradient.
Figure 2 also shows how well the numerical model reproduces the experimental static pressure behavior along each lumen, with the largest differences being localized in the radial holes of the drainage lumen (zones 1 and 5 in Figure 2C) and near the exit section of the infusion lumen (zone 4 in Figure 2D). In the former case, the Boussinesq eddy-viscosity approximation may fail due to 3D effects within the DLC. As a consequence, the strong curvature of the streamlines could not correctly account for single-point closures. Furthermore, the accuracy of the measurements may fail due to the interaction of the radial jets produced by the side holes of the DLC and the pressure probe. At the outlet section of the infusion lumen, difficulties were encountered in controlling the probe position due to flow separation from the surface of the infusion lumen's elbow-shaped exit.
Numerical results for each lumen were focused on the analysis of the 3D flow field. Distributions of WSSs and fluid dynamic variables were presented for a qualitative analysis, and peak values of computed WSS were compared with the hemolysis threshold reported by Grigioni et al.3 for quantitative purposes. Furthermore, the recirculation flows predicted by the numerical model within each lumen were investigated. Key parameters and simulated results are summarized in Table 1.
Figure 3 shows the 3D flow structures inside the drainage lumen in terms of path lines (note that the mean flow is directed along the negative x-axis of the lumen). Figure 3, A and B shows the path lines isolated to those that specifically originate from corresponding side holes, highlighting that flow through each side hole forms a circular jet that expands within the lumen and flows toward the exit. Specifically, Figure 3A shows the first (distal) group of side holes positioned in the inferior vena cava (IVC), while Figure 3B shows the second (proximal) group in the superior vena cava (SVC) (see Figure 2C, zones 1 and 5, respectively). Figure 3C illustrates the distribution of path lines near the drainage lumen inlet (corresponding to the first group of side holes); the upstream path lines (originating from the opening at the base of the cannula positioned in the IVC) flow around the inlet jets produced by the first set of lateral holes. Indeed, the path lines from both sets of lateral side holes and the opening at the base of the cannula become well mixed further upstream near the lumen exit as is shown in Figure 3D. In this latter region, the interaction between flow and lumen walls creates a wide recirculation region due to the quick increase of the lumen cross-section (see Figure 2C, zones 7 and 8).
Figure 4A shows the shear stress levels on the drainage lumen wall. A closeup of regions A and B is shown in Figures 4, B and C, which illustrates the WSS contours of the first and the second group of drainage side holes, respectively. The maximum value of shear stress (900 Pa) is localized at the sharp circular edge of the proximal set of lateral holes (section B) positioned further upstream from the drainage lumen inlet (section A). Furthermore, the analysis showed that the WSS values quickly reduced to lower levels, reaching maximum only near the sharp edges.
The simulations demonstrate that 70% of the blood drained from the IVC (28% through drainage end opening and 42% through the first group of eight lateral holes) and 30% flows through the second group of three side holes positioned in the SVC. Thus, the intake jets of the second group of side holes have higher kinetic energy than those of the first group due to a higher amount of flow per individual side hole.
Figure 5 depicts the drainage flow field around two inlet sections belonging to the first and the second group of side holes. Specifically, the plots refer to the velocity and pressure contours on a vertical plane containing the axis of the lumen as well as the diameter of the two inlet side holes. Observing Figure 5, A and B, one can recognize the dependency of the kinetic energy of the jet expanding within the drainage lumen from the intensity of the axial pressure gradients. The analysis shows that the radial pressure gradients near the wall form an accelerating flow path responsible for an increase of shear stress at the sharp corners of the holes. Furthermore, Figure 5B shows that low negative pressure exists where the strongest curvature of the streamlines takes place, that is, around the jet (−800 Pa) and near the sharp corners (−900 Pa). The value at the corner is not shown as it is localized in a very small area. Figure 5, C and D highlights that the same effects are much more evident for the second group of side holes.
Figure 6 describes the flow field within the region of the transition tube (i.e., diffusion duct in Figure 2C, zone 7) of the drainage lumen. Figure 6A shows the pressure recovery along the wall of this portion of the lumen. Upstream and downstream of this transition, static pressures are a constant −4842 Pa, but because of the device geometry, localized pressure gradients are observed. Figure 6, B-E illustrates isolines of the x-velocity component on four cross-sections along the diffusion region and the circular exit duct. The plots highlight the presence of an extended and complex recirculation region (see also Figure 3D). The separation starts near the beginning of the diffusion region due to the adverse pressure gradient and the wall friction that, at the corner of the cross-section A-A (Figure 6B), induces velocity to fall well below the superficial velocity flow. The computations show that the separation increases progressively within the diffuser (Figure 6, B and C) as the adverse pressure gradient becomes stronger. The size of the separation bubble reaches its maximum near the beginning of the circular exit duct (Figure 6D). However, downstream of this section, the recirculation reduces, and after the reattachment point (Figure 6E), pressure becomes uniform across the tube (Figure 6A).
Figure 7 shows the 3D flow structure of path lines inside the infusion lumen (note that the mean flow is directed along the positive x-axis). Figure 7, A and B illustrates the path line distribution at the S-shaped curve of the infusion lumen (zone 2 in Figure 2D). In this region, interaction between the inflow and the infusion lumen wall creates secondary flows. Proceeding from the end of the S-shaped curve toward the lumen outlet, the intensity of the secondary flow decreases as the streamlines align parallel to the infusion lumen axis. Figure 7, C and D depicts the flow field at the elbow-shaped exit (zone 4 in Figure 2D). The interaction between the flow and the wall at the elbow creates two contra-rotating vortices (Figure 7C) and a separation zone (Figure 7D), which usually lends to areas of stagnancy.
Figure 8A characterizes the WSS of the infusion lumen. Along the length of the infusion lumen, shear stress levels remain lower than the reported hemolysis threshold of 400 Pa,3 reaching peak values in the proximity of the lumen's exit (region D in the current Figure). Figure 8, B and C illustrate the WSS contours surrounding the lumen exit and S-shaped section, with the maximum shear stress value (233 Pa) identified near the lumen exit.
Figure 9 illustrates the characteristics of the flow field surrounding the S-shaped region and lumen exit of the infusion portion of the DLC. Figure 9, A-D refers to the elbow-shaped exit, whereas Figure 9, E-H refers to the S-shaped curve, highlighting the flow behavior occurring in each section. One can observe that just before the flow expansion at the lumen's exit, the flow path includes a wide separation zone at the base of the elbow-shaped exit (Figure 9A). Conversely, the flow through the S-shaped curve results in an accelerating attached flow (Figure 9E). A comparison of z-velocity components on the cross-sections AA and BB of Figure 9B and corresponding Figure 9F shows more intense secondary flows in the elbow-shaped exit than those in the region of the S-shaped curve. Furthermore, we notice that secondary flows are responsible for carrying high-energy fluid from the middle of the lumen to the walls, thereby increasing the shear stress at the wall.
Figure 9, C and D illustrates the pressure field at the elbow-shaped exit. Figure 9C shows the adverse pressure gradient that causes the flow separation of the boundary layer at the elbow, and Figure 9D indicates that the pressure increase extends for a wide region and low negative pressures exist in correspondence with the axis of the two contra-rotating vortices. Figure 9, G and H illustrates the pressure field at the S-shaped curve. Figure 9G shows the favorable pressure gradient that accelerates the flow, and Figure 9H characterizes a reduction in pressure over the cross-section BB corresponding to secondary flows but still remains positive and almost uniform as compared with the pressure distribution presented in Figure 9D.
We chose the 27-Fr DLC in this flow study because it allows for percutaneous cannulation and easily accommodates flows up to 2 L/min, which provides adequate gas exchange for most end-stage COPD patients. The experimental and numerical fluid dynamic analysis of the DLC was performed at a fluid flow rate of 2 L/min. Typical transitional flow Reynolds numbers (Re = 2000-2500) were computed within the lumens, and the computations were carried out by using the low Reynolds k-ω turbulence model. For comparison, the laminar solver was used (data not shown), but as expected, it overestimated the recirculation regions predicted within the drainage and the infusion lumen.10,11 We remark that the presence of these flow structures shows the need of local experimental data to assess the accuracy of the k-ω model in predicting the flow field in the complex geometric regions of the cannula. Although a more accurate approach might be suggested by using large eddy simulation (LES), detached eddy simulation (DES), or Reynolds stress model (RSM), the single-point closure model used in our study was selected due to the high robustness and lower computational cost for this complex computational domain. The numerical model allowed us to predict the 3D flow structures inside each lumen and identified circular jets expanding through each side hole as well as the presence of a complex 3D recirculation region inside the drainage lumen. In addition, a recirculation region and secondary flow was computed at the S-shaped curve and at the elbow-shaped exit of the infusion lumen.
Previous simulations performed on the drainage lumen by imposing velocity conditions at the lateral holes suggested that the maximum values of shear stress occurred at the inlet of the side holes.12 As a consequence, two plenums were added to the computational domain surrounding both sets of side holes to properly simulate flow in these regions. Likewise, an exhaust plenum was added surrounding the exit section of the infusion lumen to simulate the flow expansion that would occur experimentally.
Previous authors have performed studies using mathematical models to characterize the damage incurred to red blood cells after the exposure to various flow conditions.3,6,13,14 Particularly, Grigioni et al. proposed a discussion of a hemolysis threshold by applying stress analysis concepts to data originally quantified by Sallam and Hwang. In addition, Goubergrits and Affeld studied a hemolysis model and the relations between blood damage, shear loading, and exposure time, thus demonstrating that it is necessary to take into account other parameters such as the influence of turbulence, near-wall treatment, and the influence of increasing or decreasing shear stress. Each of these models can be used for comparative analysis, but they should be experimentally validated to optimize the design of artificial organs. For these reasons, this study can only provide qualitative assessments on the risk of hemolysis in the Wang-Zwische DLC.
The presence of side holes along the drainage lumen improves the fluid flow inside the cannula; in fact, ∼72% of the total mass flow rate occurs through the holes. Inside the drainage and infusion lumens, different peak values of WSS were computed. Specifically, the sharp edge of the second group of lateral holes of the drainage lumen produced a peak WSS value around 900 Pa. Lower peak values (233 Pa) were computed for the infusion lumen in proximity of the elbow-shaped exit duct probably due to the influence of the secondary flows. In comparison of the maximum WSS in both lumens of the Wang-Zwische DLC, it appears that the drainage lumen is more likely to contribute to blood hemolysis. In particular, along the infusion lumen, the shear stress levels were lower than the hemolysis threshold calculated by Grigioni et al. (400–600 Pa),3 but the drainage lumen WSS exceeded this range. Nevertheless, low values of exposure time within the high WSS regions of the drainage lumen are expected on the basis of the computed flow field.
A recirculation region was predicted in proximity of the exit sections of both the infusion and drainage lumen due to adverse pressure gradients caused by the geometry of the catheter, which could be a concern for thrombosis risk. This behavior is due to the significantly wider drainage lumen than that of the infusion lumen. From our CFD study, the Wang-Zwische DLC could be modified to provide less blood trauma and reduced thrombogenicity for long-term percutaneous paracorporeal artificial lung applications. The enlargement of the drainage holes in the SVC will increase the correspondent blood flow, thereby balancing blood drainage between the IVC and SVC. This in turn will decrease the whole drainage lumen resistance and subsequently improve overall DLC performance. Furthermore, rounding the edges of the drainage holes in the IVC and SVC will also decrease WSS. Finally, reducing the curvature of drainage/infusion lumen junctions and the angle of the elbow-shaped exit will provide streamline blood flow within the DLC, resulting in minimal recirculation and potential for thrombosis formation in those areas, as well as lower WSS and related blood trauma.
We acknowledge a few limitations in this study. First, we assumed blood to be a simple viscous fluid (Newtonian) and omitted transitional and turbulence considerations on blood. Furthermore, applying simulations with a closure turbulent model does not take into account transitional flow due to the low predictive capabilities of RANS-based models; particularly for cases of strong nonequilibrium effects, such as unsteady separation, where the simulation of recirculation flow regions may not be accurate. These limitations impact all areas of fluid dynamics, including the pressure distribution, WSSs, and the length of the separation zone predicted within the DLC cannula. However, despite these limitations, our implementation of a low-Reynolds turbulence model agreed well with our experimental data, and therefore we are confident that our conclusions would not change qualitatively on the basis of more accurate turbulence treatments or by considering the influence of a non-Newtonian blood.
In this 3D CFD study, we observed higher peak WSS values for the drainage lumen, which may potentially cause blood trauma. Furthermore, recirculation regions were predicted in the proximity of the exit sections of both the infusion and drainage lumens, which may contribute to thrombosis formation. This study provides insight for future DLC modifications in minimizing cannula-induced blood trauma and thrombogenicity in long-term applications.
1. Wang D, Zhou X, Liu X, et al
: Wang-Zwische double lumen cannula—Toward a percutaneous and ambulatory paracorporeal artificial lung. ASAIO J
54: 606–611, 2008.
2. De Wachter D, Verdonck P: Numerical calculation of hemolysis levels in peripheral hemodialysis catheters. Artif Organs
26: 576–582, 2002.
3. Grigioni M, Daniele C, D'Avenio G, Barbaro V: A discussion on the threshold limit for hemolysis related to Reynolds shear stress. J Biomech
23: 1107–1112, 1999.
4. Jegger D, Sundaram S, Shah K, et al
: Using computational fluid dynamics to evaluate a novel venous cannula (Smart Cannula) for the use in cardiopulmonary bypass operating procedures. Perfusion
22: 257–265, 2007.
5. Giersiepen M, Wurzinger LJ, Opitz R, Reul H: Estimation of shear stress related blood damage in heart-valve prostheses: In vitro comparison of 25 aortic valves. Int J Artif Organs
13: 300–306, 1990.
6. Gu L, Smith W: Evaluation of compuational models for hemolysis estimation. ASAIO J
51: 202–207, 2005.
7. Leverett L, Hellums JD, Alfrey CP, Lynch EC: Red blood cell damage by shear stress. Biophys J
12: 257–273, 1972.
8. Wang S, Rider AR, Kunselman AR, et al
: Effects of the pulsatile flow settings on pulsatile waveforms and hemodynamic energy in a PediVAS centrifugal pump. ASAIO J
55: 271–276, 2009.
9. Wilcox D: Turbulence Modeling for CFD
, 3rd ed. 2002: DCW Industries.
10. Ahmed S, Giddens D: Velocity measurements in steady flow through axisymmetric stenoses at moderate Reynolds numbers. J Biomech
16: 509–516, 1983.
11. Ghalichi F, Deng X, De Champlain A, et al
: Low Reynolds number turbulence modeling of blood flow in arterial stenoses. Biorheology
35: 281–294, 1998.
12. Colacino F, De Bartolo C, Fragomeni G, et al
: Numerical and experimental flow analysis of the Wang-Zwische double lumen cannula (Abstract). ASAIO J
55: 174, 2009.
13. Goubergrits L, Affeld K: Numerical estimation of blood damage in artificial organs. Artif Organs
28: 499–507, 2004.
14. Sallam A, Hwang N: Human red blood cell hemolysis in a turbulent shear flow: Contribution of Reynolds shear stresses. Biorheology
21: 783–797, 1984.