A three-dimensional computer model of the human heart for studying cardiac fluid dynamics
Authors: McQueen David M. , Peskin Charles S. 
Published: November 2005

In all areas of computational fluid dynamics (CFD), proper treatment of the boundary conditions is essential to computing fluid behavior correctly. In many engineering problems, CFD is simplified by a priori knowledge of the motion of the boundary. The well-known parabolic velocity profile in fully-developed flow of an incompressible Newtonian fluid in a pipe of circular cross-section is easily computed because the boundary (the pipe wall) is known to be in a fixed location. Even in more complex settings, such as flow around a ship's propeller, the motion of the boundary (the propeller) can be specified in advance.

By contrast, in most biological fluid dynamics problems the boundaries are not rigid and their motions are the result of forces imposed on them by the motion of the surrounding fluid. The motion of the fluid, of course, cannot be known without knowledge of the boundary motion. The motion of the boundary and the motion of the fluid form a coupled system; both motions must be computed simultaneously, which makes biological CFD difficult.

A particular problem of interest is the flow of blood in the chambers of the human heart. The heart is an organ whose muscular contractions pump blood around the body. Simplifying somewhat, the heart consists of two main pumping chambers that contract simultaneously. One chamber, the left ventricle, accepts oxygen-enriched blood from the lungs and pumps it to the body. The other chamber, the right ventricle, accepts oxygen-depleted blood from the body and pumps it to the lungs. The inlet and outlet of each ventricle are guarded by valves whose opening and closing guarantee one-directional flow around the circulatory system. There are a total of four valves. The valves generally consist of two or three leaflets - membranes made of very flexible but inextensible material. Familiar examples of materials with this property would be paper or fabric which can be easily bent or twisted but which are not easily stretched. One edge of each valve leaflet is securely attached to the wall of the heart, but the other edge is free of attachment and can move with the flow. Structures analogous to a valve leaflet are a shirt pocket, with one edge (three sides of a rectangular patch pocket) securely stitched to the shirt and one edge free of attachment, or a flag, one edge attached to the flag pole, the other edge free to wave in the wind. When flow is passing through the valve in the forward direction, the valve's leaflets are positioned out of the way, permitting flow. When flow attempts to pass in the reverse direction, the leaflets come together, their free edges pressing against the free edges of their neighbors to occlude the flow passage.

The motion of the leaflets is not caused by muscles in the valve. The outflow valves are entirely passive structures with no muscular tissue whatsoever. Even in the case of the inflow valves, whose free edges are connected to the heart muscle by a sparse network of tendons, the opening and closing motions result from an interaction with the surrounding fluid. The forward motion of the fluid pushes the leaflets aside out of the main flow stream, but the inextensibilty of the leaflet material prevents free motion of the fluid near the leaflet, affecting the entire flow field. Reverse motion of the fluid causes the leaflets to move back into the flow passage where contact between neighboring leaflets and the inextensibilty of the leaflet material halts the flow. The highly interactive nature of the fluid and leaflet motions makes this an especially interesting and challenging CFD problem.

Commercially available software packages intended for engineering CFD are not equipped to handle this type of dynamic interaction between boundary and fluid. We have developed a numerical method (the "Immersed Boundary Method") which simultaneously computes the motion of a fluid and the motion of an elastic boundary immersed in, and interacting with, that fluid. In the Immersed Boundary Method, the fluid is represented by Eulerian velocities and pressures that are stored on a regular three-dimensional computational lattice. The scale of the heart chambers is such that blood can be treated as a Newtonian fluid. Fluid dynamics is computed by numerical solution of the Navier-Stokes equations, including a body force. The boundary is represented by elastic structures that are free to move continuously in the space sampled by the computational lattice. The essence of the method is to replace the elastic boundary by the forces that result from its deformations. These forces are applied to the lattice in the neighborhood of the elastic boundary with the aid of a numerical approximation to the Dirac delta function. The fluid moves under the action of this body force. The numerical delta function is then used again, to interpolate the newly computed lattice velocities to the locations of the boundary, and then the boundary is moved at the interpolated velocity to a new location (the no-slip condition). The process of computing forces, then fluid motion and then new boundary location is repeated cyclically in a time-stepping procedure with a suitably chosen time step. The only requirements for the method are the physical properties of the fluid, the (possibly time-dependent) elastic properties of the boundary, and the initial geometry of the boundary. A complete description of the Immersed Boundary Method can be found in [1, 2].