FLASH INFORMATIQUE FI

High Performance Computing


Numerical simulation of hemolysis due to blood centrifugation




Mark L. SAWLEY

Alain SCHENKEL

Michel DEVILLE

Pascal HAGMANN

Jean-Denis ROCHAT


The separation of blood into its constituent components is an important procedure for a number of medical applications and various techniques and devices have been developed for this purpose. The yield efficiency of the separation process is of particular importance and depends on the requirement to minimize blood damage caused by the process. A numerical study of blood flow in a centrifugal separator has been undertaken to address this aspect. In this paper, the physical and numerical modelling of the blood flow is presented, with particular emphasis on the use of high performance computing to determine the complex flow fields.

Introduction

Different devices involving blood processing have been developed for various medical applications, such as rotary blood pumps, hemodialysis cannulas and ventricular assist devices. A common and much needed medical procedure is the separation of blood into its constituents for purposes of storage or transfusion. Since the blood constituents - red and white blood cells as well as platelets suspended in plasma - have different densities, a common separation technique is to accelerate their sedimentation by centrifugation and then collect the different constituents separately. One aspect that must be kept under control is the level of blood cell damage resulting from the separation process. The membrane of RBC (red blood cells), for instance, can break when subjected to conditions of excessive stress. This phenomenon, termed hemolysis, has various negative consequences: it reduces the hematocrit level, the proportion of the blood volume occupied by red cells, and thus the yield of a separator device; it also pollutes the blood plasma, as the hemoglobin contained in the RBC is released after breakage of the cell membrane. Strict regulations exist for the allowed level of hemoglobin in plasma samples destined for transfusion.
Over the past decade, increasing use has been made of CFD (Computational Fluid Dynamics) to understand and optimize blood processing devices. CFD involves the numerical resolution of the partial differential equations governing the, possibly unsteady, flow behaviour. It can be employed to reduce the time and cost of research and development involved with the manufacture of prototypes. CFD has become increasingly applied in recent years with the advent of advanced flow simulation software and relatively inexpensive commodity compute clusters. After validation of a CFD implementation, new device designs or modifications can quickly be analyzed for improvement. This is of particular interest when strong industrial constraints arise, such as the need for disposable devices that must be produced at low cost, as is common in medical applications for health safety reasons.
At the beginning of 2008, a joint project was initiated between Biofluid Systems, a Swiss company based in Nyon specialized in the conception of medical devices, and the Laboratory of Computational Engineering at EPFL. This project is aimed at improving the performance of a blood separator device, in particular, through minimization of the blood damage. The device is a centrifuge comprised of two parts: an upper part whose section is shown in Figure 1 and the centrifugation cylinder itself (lower part, not shown). The device is a few centimetres in diameter and rotates at a speed of several thousand revolutions per minute. Blood enters through the top and flows through the central canal to the bottom of the centrifugation cylinder where sedimentation occurs. During the sedimentation process, the concentrated RBC phase (shown in red) moves up along the internal face of the outer cylinder and reaches the upper part where collection of this phase takes place.
The collection mechanism is based on the pump principle. As the external sides of the collection chamber are part of the centrifuge, the fluid along the walls rotates at high velocity. When the fluid makes contact with the collectors (two stationary facing disks inside the collection chamber, shown in grey in Figure 1), its strong deceleration yields an increase in pressure, causing the fluid to exit by the centre (left on Figure 1).

JPEG - 13.3 kb
fig. 1
Top and collection chamber of blood separator device

This straightforward design allows, in particular, for low production cost. The collection mechanism, however, imposes harsh conditions on the RBC phase. In this project, the goal was to use CFD simulations to achieve a better physical understanding of the flow fields inside the collection chamber and the resulting hemolysis of the RBC phase. It should be noted that the flow in the collection chamber is biphasic involving a free surface interface between the immiscible RBC and air phases. While several CFD analyses of blood processing devices have been performed (see e.g. [1] and [2]), we are not aware of any study involving a free surface. The subsequent estimation of hemolysis is a challenging problem that requires the modelling of the interaction between the RBC and the plasma. As hemolysis is a result of the influence of excessive stress on the RBC, it is important to determine accurately the flow fields in the collector.

Physical and numerical modelling

In this study, blood is represented as a single-phase incompressible Newtonian fluid. Its viscosity is considered independent of the flow, varying only, as the density, with the hematocrit level of the sample being simulated. For typical operating rotational speed, the flow in the collection chamber reaches a turbulence regime, with a Reynolds number of the order of 105. The turbulence is modelled using the RNG k-e model, which takes into account the effect of swirl on turbulence and has better convergence properties for swirling flows than alternative models.
An important aspect of the flow in the collection chamber is the presence of an air-blood interface whose shape and dynamics are an integral part of the solution. The fluid occupying the full flow domain is modelled using the volume-of-fluid scheme, a computationally efficient interface-capturing method for free-surface flows.
Another aspect of the physical modelling concerns the estimation of hemolysis. In the present study, we use a mean field approximation, where the local increase in free hemoglobin is represented as a function of macroscopic quantities computed in the CFD simulations. The quantities that are naturally considered are the viscous shear stress and the Reynolds stress. Another important quantity is the exposure time, as hemolysis depends on the duration of the mechanical load. Several CFD hemolysis studies have been reported in the literature (see e.g. [1] and [2] and the references therein). Since, to the best of our knowledge, all the previously studied devices function on very different operating principles than the device considered here, it proved necessary to design our own specific model. A four parameter model was derived to determine the rate of production of hemolysis along a trajectory, involving local shear stress as well as load history. The particular dependence on shear stress follows from previous studies [3].
A commercial software package that solves the governing Reynolds-averaged Navier-Stokes equations using a finite volume discretization was employed in the present study. Particular care was taken in initializing the flow, as it involves high rotational speed as well as large velocity gradients. Most of the simulations were performed as fully unsteady problems to ensure reliability in the stability of the observed solutions.

Computational issues

To provide an indication of potential improvements in the device design, the CFD simulations and analysis must be sufficiently accurate to capture small effects. Computationally, this requires high resolution computations as well as the use of detailed - and therefore expensive - flow models. Two types of simulations were undertaken: high throughput parametric studies involving different device geometries and high resolution in-depth computations for a few selected cases. It was essential to use high performance computing for both of these types of simulations.
High performance computing nowadays generally involves the application of appropriate techniques to exploit available parallel computer systems. Two different techniques have been employed in the present study: computing in parallel and parallel computing. All simulations were performed on the Pleiades cluster at the EPFL [4], which consists of commodity processors interconnected by a Gigabit Ethernet switch. In its present configuration, the Pleiades cluster comprises a total of 960 compute cores with 120 Intel Xeon mono-processor nodes (2.8 GHz, 4 GB dual-access DDR memory) and 210 Intel Woodcrest bi-processor bi-core nodes (2.67 GHz, 8 GB dual-access DDR memory).

Computing in parallel:
For the parametric studies, high throughput simulations were performed based on axisymmetric swirl space, namely, full dependence of the velocity on the radial and axial directions while no dependence on the azimuthal direction. This involves computing a two-dimensional (2D) flow augmented by one scalar field for the azimuthal velocity. Using 2D meshes of about 110,000 cells, typical cases required 230 MB of RAM and satisfactory convergence to a stationary solution was achieved in about 30 hours of CPU time. The parametric studies required running hundreds of cases, to test various geometries and variable values. Each case was computed on a single processor core, using an automated procedure designed for the simultaneous submission of multiple cases covering a given range of parameters.
Parallel computing:
High resolution in-depth computations have also been performed of full 3D flows. Such simulations were aimed at exploring the eventual existence of azimuthal flow structures. Since such simulations are computationally very demanding, only a few selected cases were computed. The largest case that was simulated involved a computational mesh of about 4 million cells. A section of a mesh containing a total of 1 million cells is shown in Figure 2.
JPEG - 14.7 kb
fig. 2
Section of a 3D mesh used for high resolution computations

To undertake the parallel computations, domain decomposition was used to distribute the work load involved in the simulation on multiple processor cores. A case with 4 million mesh cells occupies 1.3 GB of disk space and requires up to 9 GB of RAM. It therefore cannot be run on a single node without a significant loss in performance due to the necessity of swapping. Simulations of large cases were performed on up to 8 quad-core nodes, that is, 32 cores.
Under certain circumstances, unsteady 3D solutions have been observed, which required a larger computational effort to simulate. As an indication, a case with 4 million elements running on 16 processors required a wall time of 55 hours to compute 5 milliseconds of real-time operation.

Simulation results

Figure 3 shows the result of an axisymmetric swirl simulation with the entire (2D) computational domain displayed. The inlet is at the lower right hand corner of the figure, while the outlet is on the left close to the axis of the centrifuge. At a sufficient distance from the collectors, the flow is seen to maintain a high speed, essentially solid-body rotation. As the collectors are stationary, the fluid is abruptly slowed down at their tips, and moves along their internal faces towards the axis of the centrifuge. A transition between a swirl-dominated flow and a radial-dominated flow occurs in the neighbourhood of the collectors’ tips. Local increase in hemolysis takes place in this region due to the presence of high shear.
To determine the hemolysis that is produced as the blood flows through the device, it is necessary to compute the trajectories of individual volume elements in order to track the flow conditions encountered by red blood cells along their path.
Displayed in Figure 3 is also the azimuthal projection of a selected single trajectory. In addition to the high-speed rotation in the swirl direction, a vortex-like motion is observed in the perpendicular plane. This indicates that recirculation takes place in the collection chamber. A similar vortex type structure also arises in the upper part of the chamber, and hence the majority of the RBC phase is comprised of recirculation zones. Since hemolysis depends on the exposure time of the RBC to high stress regions, the residence time of individual cells within such vortices is an important quantity.

JPEG - 9 kb
fig. 3
Computed swirl velocity (background colour) in the collection chamber. The volume occupied by the air phase is shown in grey. The azimuthal projection of a recirculating trajectory is indicated in black

The 3D flow simulations have revealed the existence of various types of pathlines. In Figure 4(a) is presented a top view of a 3D flow trajectory starting near the inlet and integrated for 5 seconds. The trajectory keeps rotating within the recirculation region without entering the collector interspace, even though the average residence time in the collection chamber, determined by the ratio of volume to flux, is about 1 second.
Other trajectories, such as that shown in Figure 4(b), exit quickly after a few cycles within the outer part of the rotating chamber. Once between the collectors, the fluid quickly slows down and exits the device in a spiral-like motion. The trajectory shown in Figure 4(b) is coloured according to the level of hemolysis accumulated in the corresponding volume element. As expected, hemolysis abruptly increases near the collector tip where the fluid slows down.
Other trajectories determined from the 3D simulations display more complex flow behaviour. The example presented in Figure 4(c) shows a pathline that enters the collector interspace where it follows a chaotic path before returning to the collection chamber rather than exiting the device through the central outlet. This behaviour suggests that some trajectories undergo multiple accelerations and decelerations, thus resulting in an increase in hemolysis.

JPEG - 23 kb
fig. 4
Top view of 3D flow trajectories in the collection chamber: (a) recirculating trajectory, (b) directly exiting trajectory, and (c) complex trajectory. The device outlet is on the border of the central region in white. The collector disks are in dark grey, while the outer region in light grey is the volume where essentially solid-body rotation occurs

Conclusions

The numerical simulation results obtained in the present study suggest that when flowing through the device, individual RBC may encounter quite different flow conditions and exposure times. An accurate estimate of the global hemolysis thus requires an adequate sampling of the trajectory space. This in turn requires the capability to reproduce accurately not only local but also global flow properties. The associated issue of performing accurate Lagrangian tracking in unsteady complex flows is also important; this is a very challenging problem from a computational point of view.
While the ability to perform reliable CFD analysis as part of the design process may yield significant cost reduction for an industrial product development, the present study illustrates that it may also require significant computational resources. While the capability of high performance computers will continue to increase, gaining the competitive edge through computer-assisted engineering is likely to remain a competence-demanding task. This emphasizes the continued importance of appropriate partnerships between academia and industrial developers.

Acknowledgements

This work has been performed within the framework of a project funded by the Swiss Commission for Technology and Innovation (CTI).

References

[1] M. Behbahani, M. Behr, M. Hormes, U. Steinseifer, D. Arora, O. Coronado and M. Pasquali, A review of computational fluid dynamics analysis of blood pumps, European Journal of Applied Mathematics, 20, 363-397 (2009). [2] L. Gu and W.A. Smith, Evaluation of computational models for hemolysis estimation, ASAIO Journal. 51, 202-207 (2005). [3] S.A. Jones, A relationship between Reynolds stresses and viscous dissipation: implications to red cell damage, Annals of Biomedical Engineering, 23, 21-28 (1995). [4] Pleiades cluster; see pleiades.epfl.ch.



Cherchez ...

- dans tous les Flash informatique
(entre 1986 et 2001: seulement sur les titres et auteurs)
- par mot-clé

Avertissement

Cette page est un article d'une publication de l'EPFL.
Le contenu et certains liens ne sont peut-être plus d'actualité.

Responsabilité

Les articles n'engagent que leurs auteurs, sauf ceux qui concernent de façon évidente des prestations officielles (sous la responsabilité du DIT ou d'autres entités). Toute reproduction, même partielle, n'est autorisée qu'avec l'accord de la rédaction et des auteurs.


Archives sur clé USB

Le Flash informatique ne paraîtra plus. Le dernier numéro est daté de décembre 2013.

Taguage des articles

Depuis 2010, pour aider le lecteur, les articles sont taggués:
  •   tout public
    que vous soyiez utilisateur occasionnel du PC familial, ou bien simplement propriétaire d'un iPhone, lisez l'article marqué tout public, vous y apprendrez plein de choses qui vous permettront de mieux appréhender ces technologies qui envahissent votre quotidien
  •   public averti
    l'article parle de concepts techniques, mais à la portée de toute personne intéressée par les dessous des nouvelles technologies
  •   expert
    le sujet abordé n'intéresse que peu de lecteurs, mais ceux-là seront ravis d'approfondir un thème, d'en savoir plus sur un nouveau langage.