ABSTRACT

Regurgitation in the heart diastolic phase represents a critical flow condition associated with many heart valve design considerations. The finite volume method, the Low-Reynolds-Number k-ω turbulent model and sliding mesh model are employed to solve and compare the complex flow field and the torque in each case. The end results expected from a cardiovascular CFD analysis are not limited only to the flowfield investigations. More importantly, it needs an evalua- tion criterion to judge if the design is acceptable as considered from a broader blood cell damage or activation perspec- tive. In this study, blood cell damage index developed based on stress-time empirical rule and Lagrangian particle tracking is introduced to assess the viscous and turbulence-induced stresses effect to the blood cells.

Keywords:

Heart Valve; CFD; Blood Cell Damage Index; Stress; Lagrangian Particle Tracking

1. Introduction

Follow an ideal prosthetic heart valve is expected to satisfy the following requirements including: 1) low or subclinical thrombogenicity and blood cell damage; 2) low obstruction or impedance to the circulating blood flow; 3) minimal valve sound generation; 4) lifetime du- rability; 5) biocompatibility. It can be seen that except for the last requirement that is associated with the mate- rials, all the other demands are related to the fluid dy- namics of the transvalvular blood flow. The red blood cell damage and platelet activation, which in turn may lead to hemolysis and thrombosis, are highly related to the wall shear and turbulence Reynolds stresses [1,2]. Sutera and Mehrjardi [3] found that hemolysis occurs when Reynolds shear stress reaches 150 N/m² for an ex- posure time of 10² s, and loss of cell bioconcavity shows up at around 50 N/m². Platelet would also be activated or damaged by excessive shear stress. In-vitro experiments have shown that platelet lysis and the induced blood ag- gregation were caused by shear stress of the order of 10 - 50 N/m² [4,5]. Moreover, the endothelial cells on the vessel wall would be eroded by wall shear stress greater than 40 N/m² [6,7]. It is obvious that, in the design of prosthetic valves, minimizing valve-induced wall shear and Reynolds stresses should be kept as the prioritized hemodynamic design objective.

Blood damage is quantified by the hemolysis index (HI), which is calculated from photometric plasma he- moglobin measurements. Blood damage by fluid forces primarily depends on magnitude and duration of loading. Especially, shear forces within the flow field have dam- aging effects on blood corpuscles [8]. Prior experimental approaches for the quantification of flow-induced blood damage cover a wide range of shear levels and exposure times. The dependency from the variables exposure time and shear stress is characterized by [9] via a 2D regres- sion analysis of the data of Giersieppen et al. [10] A sig- nificant increase in blood damage can be observed for shear stresses of t ≥ 425 Pa and exposure times of t_exp ≥ 620 ms. Maximum hemolysis within the investigated range is HI = 3.5%.

The development of computational models of the fluid mechanics in heart valves is motivated by medical and economic concerns as well as efficiency. From 2D to 3D problem solution is challenging for the current generation of CFD (Computational Fluid Dynamic) software and computer hardware. However, the difficulty stems from the fact that the fluid mechanics of heart valves includes several phenomena that are difficult problems in their own right, such as fluid-structure interaction in which the motion of the valve leaflets drive most of the fluid mo- tion and the wide disparity in length scales from the large-scale motions with length scales in the valve orifice (∼ 25 mm ) to the flow in the hinge regions of mechanical valves (∼100 µm) [11]. 2D models of the flow through heart valves were developed in the 1970s using the vor- ticity-stream function formulation of the Navier-Stokes equation, which is convenient and efficient in two di- mensions but unwieldy in three dimensions. In these studies, the geometry of the native and prosthetic heart valves was greatly simplified because the numerical methods were limited to simple Cartesian grids [12,13]. To predict two-dimensional fluid-structure interactions, Peskin [14] developed an immersed boundary method to couple the blood flow to the structure of the heart and its valves. More recently, de Hart et al. [15] developed a 2D fictitious domain finite element code to predict the un- steady flow in the aortic valve.

In the mid-to-late 1990s, the computer resources avail- able to heart valve researchers increased to the point at which three-dimensional simulations of MHVs could be performed. King et al. [16] extended their two-dimen- sional simulations to three dimensions by using one quar- ter of the valve geometry. These simulations used a si- nusoidally varying inlet velocity to mimic the accelera- tion phase of systole with a peak Re = 3000. Peskin & McQueen extended their immersed boundary method to 3D and successfully computed the flow in the heart valves [17]; however, difficulties still exist with extend- ing their method to predict the turbulent flow that often occurs in diseased conditions or MHVs.

Because the MHV is moveable, interactions between blood flow and valve can involve a rich range of fluid- mechanical phenomena. The flow will affect movement of the valve and valve movements in turn influence the flowfield. Hence, simultaneous fluid-structure interac- tions (FSI) should be considered when studying the he- modynamics and biomechanics of MHV.

In this study, a fully-coupled fluid-structure interaction software system for a pulsatile flow across a moving tilt- ing-disc valve is developed and applied. A single-degree- of-freedom rotational occluder model is integrated simul- taneously with the CFD time-stepping. The blood dam- age index is introduced by the Lagrangian particle track- ing to assess the turbulence-induced stresses to the blood cells.

2. Numerical Model

The generic valve flow model considered consists of a flat-plate occluder and a straight duct of infinite extent which represents the vascular bed. The occluder has a chord length C, which nearly equals but differs with the duct height H by a clearance of 0.0015 C . This generic occluder has a thickness of 0.05 C and is pivoted at a point located at the quarter chord measured from the oc- cluder leading-edge. Physically, this generic model si- mulates the tilting-disc valve flows. It can be extended to study the bi-leaflet valve flow easily by a straightfor- ward addition of another pivoted occluder and duct into the computational domain.

The computational domain and coordinate system is illustrated in Figure 1. The coordinate origin is located on the lower wall with its abscissa passing through the pivot center. The opening angle is defined positive in the clockwise sense with zero-angle specified in the vertical position when occluder aligns with the abscissa. For simulating the valve flows, four blocks with a total di- mension of 20 duct heights are used presently. Except for Block III, all other blocks are filled by structured grids. In the vicinity of the occluder and wall surfaces there exist layers of structured grids that are employed to re- solve the boundary layer flows. This boundary layer grid system, either moving dynamically with the occluder or staying stationarily as attached to the duct wall, is fixed in shape. Adjacent to these structured grid systems, however, are unstructured grids which fill up the rest of the space in Block III and are updated in each time- marching step as a result of valve motion. This adaptive multi-block girding strategy reduces the labor of grid update required in the dynamic grid simulation whereas solves the problem of grid skewing and improper clus- tering or coarsening of the unstructured grids associated with large solid surface movement or topology change. Figure 2 depicts the hybrid structured and unstructured grid system at a typical valve opening angle. It is ob- served that appropriate grid clustering assures the flow details be properly resolved.

Specification of numerical boundary conditions is il- lustrated in Figure 3. No-slip condition is enforced at the solid boundary including occluder surface and vessel walls.

3. Mathematical Modeling

3.1. Dynamic Mesh Conservation Equations

The Navier-Stokes equations in integral form for an arbi- trary control volume Ω with differential surface area, whose boundary is moving can be written as:

(1)

where ρ is the fluid density, u is the flow velocity vector, u_g is the grid velocity of the moving mesh, Γ is the diffu- sion coefficient, S_f is the source term of f.

In order to satisfy the grid conservation law [18] the volume time derivative of the control volume is com- puted from

(2)

3.2. Structure Dynamic Equation

The valve rotation can be described by the relationship,

(3)

To obtain the flow field, the valve-closing process can be divided into incremental time steps. In each step, the blood pressure and the friction force computed from the previous time step is integrated over the leaflet surface to compute the moment due to the pressure and shear stress. The computational process shows in Figure 4. For fluid solver, a commercial CFD package, Fluent, is employed. The above described method was programmed as a user- defined function in Fluent. And Fluent-subroutine “Dy- namic Grid” was used for computational domain re- meshing.

The motion of the solid-body can be specified as a UDF code. The valve rotation Equation (3) can be de- scribed by the relationship,

(4)

Equation (4) is a linear function so we increment the “adaptive time step method” in the UDF code as show in Figure 4. Using the adaptive time step method lets the determined pitch moment the same as in fluid and struc- ture solver.

3.3. Numerical Method

A pressure based finite volume method was used to solve the Equation (1). The contributions from the transient, total flux and source term in Equation (1) are assembled as flows. For each control volume the convection contri- butions from the adjacent volume as given by Equation (2) are summed and with the diffusion contributions and source contributions are substituted into Equation (1) to yield a system of equation of the from

(5)

which may be solved for the dependent variable f using any suitable solver.

The pressure correction equation is derived from the continuity equation where the velocity terms are substi- tuted with the velocity correction terms [19], given by

(6)

where the “start” velocities u and v are the guessed ve- locities at the previous iteration. Integrating the continu- ity equation over an arbitrary control volume and apply- ing Green’s first theorem gives

(7)

On application of the geometric conservation law [18]

(8)

Using the corrected pressures by Equation (6) the ve- locity corrections are calculated and hence the corrected velocity is calculated and applied to the equation for convective mass flux and rearranging gives the system of equations which for the control volume surrounding node P may be expressed as

(9)

Equation (9) may be solved using any suitable solver with the resulting solution used to update the variables in Equation (1).

3.4. Turbulent Model

Numerical and experimental results show that the time averaged Reynolds number in arteries is not high (300 - 500); however, when the waveform is up to peak value, Re_max is near 2000 - 2300. Thus, the flow may become turbulent and hence a modified version of the Low- Reynolds-Number (LRN) k-ω model of Wilcox (1998) was selected [20].

4. Boundary Conditions

For the present pulsatile flow simulation, a prescribed inflow volume flowrate history is specified at the inflow plane. This aortic volume flow rate waveform, as shown in Figure 5, is adopted from an in-vitro experiment car- ried out by Simon et al. [21]. The present incompressible flow simulation requires a given velocity distribution be specified at the entrance plane. This velocity profile is unknown in advance although the volume flow rate is given. It is assumed in the present simulations that fully- developed pulsatile inflow is established in block I. For all the subsequent unsteady simulations, the inflow ve- locity distribution is calculated in a sub-iteration manner over block I before marching to the next time step. At each instant, a fully-developed parabolic laminar pipe flow velocity distribution corresponding to the specified volume flow rate is taken as the initial guess. Periodic boundary condition is enforced over the inflow and out- flow planes of block I during the sub-iterated pseudo time-marching until certain convergence criterion is sat- isfied.

Fully-coupled fluid/structure interaction is carried out by integrating together the fluid and structure equations. Since the volume flowrate is given, the occluder motion will be determined as a result interacted with the passing pulsatile flow. How would be the maximum opening angle reached for the occluder in heart systole is deter- mined by the specified inflow volume flowrate condition. In order to make the interaction simulation tractable, the impact dynamic modeling is not attempted presently. The closure criterion, therefore, is set by neglecting the asso- ciated occluder rebounding that might occur in the real valve flow.

5. Results and Discussion

5.1. Regurgitant Flow in the Diastolic Phase

Regurgitation in the diastolic phase represents a critical flow condition associated with many valve design con- siderations. Hemolysis and cavitation inception are two serious problems found with the regurgitant flow for mechanical heart valves. The high-speed jet emanating from the tiny clearance between valve edge and housing rim constitutes a high-shear and low-pressure flow char- acteristic which is prone to cavitation inception. The burst of the cavitation bubbles is very detrimental to blood cells and the valve structure as well. Experimental investigation of this regurgitant flow is difficult to per- form and CFD analysis becomes the alternative re- searchers are interested in. It was argued that the dimen- sion disparity, namely, the clearance size is of orders of magnitude smaller than the occluder diameter, allows

clearance flow be treated as a local phenomenon. The interaction with the much larger transvalvular main stream is assumed weak so almost all the investigations, to the author’s knowledge, were conducted based on a local viewpoint only. The boundary conditions specified on the two sides across the occluder are quiescent flows with given steady or transient pressure differential. The occluder was assumed either stationary or moving with a prescribed closing velocity. This regurgitant valve flow is herein revisited using the fully coupled FSI code. The appropriateness of taking regurgitant flow as of local and quasi-steady nature is re-evaluated.

There are 4 instants that are taken over the diastolic phase of the volume flow rate trace, representing differ- ent stages of the regurgitant flow. Two observation planes are selected to explore the regurgitant flow details. The first cross-sectional plane A-A is taken one channel height upstream of the occluder. The second plane, B-B, however, is defined in the clearance gap close to the jet exit. Steady-state simulation is carried out using the last instant (t = 800 ms) volume flowrate as the inflow condi- tion. The axial velocity distributions obtained for all in- stants, either simulated using unsteady or steady flow assumptions, show almost identical profiles at the B-B cross sectional plane, as illustrated in Figure 6. This im- plies that in the clearance gap, indeed, flow properties are local and quasi-steady. However, away from the jet exit, for example, at the downstream A-A plane, unsteady and quasi-steady simulations yield completely different re- sults.

Shown in Figure 7 is the axial velocity distribution. Similarity is found among the solution group simulated using the unsteady, coupled fluid/structure solver. Sharp contrast, however, exists between results computed under quasi-steady and unsteady assumptions. Intuitively it would be reasonable to expect that the unsteady result obtained at the late valve closure instant (t = 800 ms)

would agree more with the quasi-steady result because the transient induced by valve closing should largely diminish or die out then. This physical intuition is proved to be wrong in Figure 8 using velocity vector plots. It is seen that the vortices obtained by the coupled fluid/ structure simulations encompass completely different fun- damental structures, suggesting that the flow characteris- tics before valve closure is essential.

Note that the quasi-steady simulation has a uniform, quiescent background flow to start with. In consequence, the clearance jets develop into a pair of symmetric, flat- tened vortices with high energy concentrated in the near wall regions. The wall shear stress which is important for predicting endothelial cell lesion is hence overestimated in the quasi-steady simulation.

Comparing to the quasi-steady result, the unsteady simulations show a more round and dominant vortex structure swirling around the occluder. All unsteady simulations, pertaining no matter to early or late instant in the diastolic phase, possess similar vortex structures. These vortices, in fact, are initiated long before the com- pletion of valve closure, attributed to the occluder closure motion. Right before the valve closure, a concentrated vortex core is produced in the upstream vicinity of the occluder trailing-edge. Lai et al. [14] also found a similar low-pressure vorticity zone occurring at similar place during valve closure, with vorticity strength dependent on the closure speed. When the occluder comes to a com- plete stop at the closure position, the clearance jet flows

develop and interact with this low-pressure vortex struc- ture. The upper jet energy was absorbed by this vortex, resulting in an even larger vortex resident next to the regurgitant flow side. The persistence of this lingering vortex diffuses the velocity shear originally contained in the wall jet regions. This explains why the unsteady simulation results have a mild velocity gradient distribu- tion which persists without decay over the entire diastolic phase.

This closure vortex and wall jet interaction phenome- non is newly found in the present fluid/structure interac- tion simulation. The previous concept that clearance flow characteristics can be taken as local and quasi-steady should be re-evaluated.

A viable leakage flow study cannot ignore the initial flow structure which is already in existence prior to the complete valve closure. The role of the trapped, lingering closure vortex should be included in the future regurgi- tant flow analysis. Otherwise, the structure of the valve leakage flow will be misinterpreted and the clinical im- plication could be wrongly overestimated.

5.2. Blood Cell Damage Evaluation

The results expected from a cardiovascular CFD analysis are not limited only to the flowfield investigations. More importantly, it needs an evaluation criterion to judge if the design is acceptable as considered from a broader blood cell damage or activation perspective. Researches [22-25] found that both shear stress and exposure time are important factors for accounting the breakdown of red blood cell membrane, which would cause protein hemoglobin to be released into the surrounding plasma. Although there are some disagreements existing in the calibration of the parameters involved in the empirical rule.

A general consensus is reached in controlled experi- ments that the ratio of damaged hemoglobin content to the total hemoglobin content can be expressed by a power law (proposed by Giersiepen et al. [9]):

(10)

In the present CFD simulation a cumulative blood cell damage assumption is adopted. The hemolysis index HI is therefore defined using the following integral estimate [10],

(11)

Cells are taken as passive markers carried by the stream. Lagrangian tracking is employed to define the instantaneous cell position and hence the stress experi- enced at the instant of observation. In tracking the parti- cles, residence time can also be obtained by adding up the time steps a particle has traveled since released from the inlet until discharged out of the exit. Longer resi- dence time usually correlates with higher potential of blood damage and thrombogenesis.

The inflow and exit planes for evaluating blood trauma are selected to be located at X = −10H and 10H, respec- tively. There are totally 40 particles that are released from the inlet plane and the accumulated HI values asso- ciated with each particle path line. Large HI values are found in general for particles drifted through the high- shear regions.

Since regurgitant flow is of particular significance to blood trauma, particle tracers are also released at the re- gurgitant phase. The results are shown in Figure 9. The regurgitated particles are squeezed through the upper and lower clearance gaps and the pathlines are greatly influ- enced by the vortex motion on the two sides of the leaflet. The vortex entrainment effect is best seen around t = 500 ms, and all particles leave the closure vortices eventually, as seen from the snapshot presented corresponding to t = 800 ms. During the particle traveling period the closure vortices do not change too much in their strength and position (see Figure 8), working as a cleaner to wash out the particles trapped in the leeward vicinity of the leaflet.

6. Conclusions

A fully-coupled fluid/structure interaction solver is de- veloped to simulate the transvalvuar flow passing over a tilting-disc prosthesis valve. Numerical simulation of a flat-plate occluder in an infinite straight duct is then conducted to explore the flowfield as interacted by an immersed moving occluder.

A pulsatile aortic volume flowrate history obtained from in-vitro test is specified as the inflow boundary condition. The resultant occluder motion is solved and the pitching moment trajectory shows a hysteresis loop when plotted against the opening angle, indicating dif- ferent transient dynamics is involved in the valve open- ing and closure phases. Flowfield details are disclosed using turbulent intensity, velocity, streamline, and vortic- ity snapshots corresponding to various instants of a pul- satile cycle. Complex vortex structures are found in the valve opening and closure phases. It is found in the pre- sent fully-coupled fluid/structure simulation that the vor-

tices generated in the valve closure and regurgitant pe- riod encompass important mechanism which was over- looked in the past. During valve closure, the swing of the leaflet in the decelerated background stream generates different vortices around the two sides of the occluder. A low-pressure vortex core is further induced on the back side around the trailing-edge when leaflet is fully closed. This closure vortex is trapped around the occluder while interacting with the regurgitant jet flow. Contrary to the unsteady results, symmetric flattened vortex pair is ob- served in the quasi-steady simulation of the regurgitant flow, and the shear stress is overestimated. In summary, it is concluded that only fully-coupled fluid/structure si- mulation can yield correct description of a pulsatile transvalvular flow, in particular for analyzing the regur- gitant flow in the diastolic phase.

References