Dynamic Simulation and Hemolysis Evaluation of the Regurgitant Flow over a Tilting-Disc Mechanical Heart Valve in Pulsatile Flow

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 evaluation 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.


Introduction
Follow the 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 durability; 5) biocompatibility. It can be seen that except for the last requirement that is associated with the materials, all the other demands are related to the fluid dynamics 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 2 for an exposure time of 10 2 s, and loss of cell bioconcavity shows up at around 50 N/m 2 . Platelet would also be activated or damaged by excessive shear stress. In-vitro experiments have shown that platelet lysis and the induced blood aggregation were caused by shear stress of the order of 10 -50 N/m 2 [4,5]. Moreover, the endothelial cells on the vessel wall would be eroded by wall shear stress greater than 40 N/m 2 [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 hemoglobin measurements. Blood damage by fluid forces primarily depends on magnitude and duration of loading. Especially, shear forces within the flow field have damaging 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 regression analysis of the data of Giersieppen et al. [10] A significant increase in blood damage can be observed for shear stresses of τ ≥ 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 motion 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 vorticity-stream function formulation of the Navier-Stokes equation, which is convenient and efficient in two dimensions 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 unsteady flow in the aortic valve.
In the mid-to-late 1990s, the computer resources available 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-dimensional simulations to three dimensions by using one quarter of the valve geometry. These simulations used a sinusoidally varying inlet velocity to mimic the acceleration 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 extending 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 fluidmechanical phenomena. The flow will affect movement of the valve and valve movements in turn influence the flowfield. Hence, simultaneous fluid-structure interactions (FSI) should be considered when studying the hemodynamics and biomechanics of MHV.
In this study, a fully-coupled fluid-structure interaction software system for a pulsatile flow across a moving tilting-disc valve is developed and applied. A single-degreeof-freedom rotational occluder model is integrated simultaneously with the CFD time-stepping. The blood damage index is introduced by the Lagrangian particle tracking to assess the turbulence-induced stresses to the blood cells.

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 occluder leading-edge. Physically, this generic model simulates the tilting-disc valve flows. It can be extended to study the bi-leaflet valve flow easily by a straightforward 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 dimension 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 resolve 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 timemarching 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 clustering 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 observed that appropriate grid clustering assures the flow details be properly resolved.
Specification of numerical boundary conditions is illustrated in Figure 3. No-slip condition is enforced at the solid boundary including occluder surface and vessel walls.

Dynamic Mesh Conservation Equations
The Navier-Stokes equations in integral form for an arbitrary control volume Ω with differential surface area , whose boundary is moving can be written as:  where ρ is the fluid density, u is the flow velocity vector, u g is the grid velocity of the moving mesh, Γ is the diffusion coefficient, S φ is the source term of φ.
In order to satisfy the grid conservation law [18] the volume time derivative of the control volume is computed from

Structure Dynamic Equation
The valve rotation can be described by the relationship, 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 userdefined function in Fluent. And Fluent-subroutine "Dynamic Grid" was used for computational domain remeshing.
The motion of the solid-body can be specified as a UDF code. The valve rotation Equation (3) 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 structure solver.

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 contributions 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  which may be solved for the dependent variable φ using any suitable solver. The pressure correction equation is derived from the continuity equation where the velocity terms are substituted with the velocity correction terms [19], given by where the "start" velocities u and v are the guessed velocities at the previous iteration. Integrating the continuity equation over an arbitrary control volume and applying Green's first theorem gives On application of the geometric conservation law [18] ( ) Using the corrected pressures by Equation (6)

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].

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 carried 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 velocity 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 outflow planes of block I during the sub-iterated pseudo time-marching until certain convergence criterion is satisfied.
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 determined 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 associated occluder rebounding that might occur in the real valve flow.

Regurgitant Flow in the Diastolic Phase
Regurgitation in the diastolic phase represents a critical flow condition associated with many valve design considerations. 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 characteristic 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 perform and CFD analysis becomes the alternative researchers are interested in. It was argued that the dimension 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 different 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 condition. The axial velocity distributions obtained for all instants, 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 implies 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 results.
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 fundamental structures, suggesting that the flow characteristics 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, flattened 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 completion 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 complete stop at the closure position, the clearance jet flows develop and interact with this low-pressure vortex structure. 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 phenomenon is newly found in the present fluid/structure interaction 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 regurgitant flow analysis. Otherwise, the structure of the valve leakage flow will be misinterpreted and the clinical implication could be wrongly overestimated.

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][23][24][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 experiments 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]): 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 Cells are taken as passive markers carried by the stream. Lagrangian tracking is employed to define the instantaneous cell position and hence the stress experienced at the instant of observation. In tracking the particles, 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 residence 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, respectively. There are totally 40 particles that are released ciated with each particle path line. Large HI values are found in general for particles drifted through the highshear regions.
Since regurgit from the inlet plane and the accumulated HI values asso-ant flow is of particular significance to bl

Conclusions
ucture interaction solver is deflowrate history obtained fr ood trauma, particle tracers are also released at the regurgitant 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 influenced 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.
A fully-coupled fluid/str 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 om 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 different transient dynamics is involved in the valve opening and closure phases. Flowfield details are disclosed using turbulent intensity, velocity, streamline, and vorticity snapshots corresponding to various instants of a pulsatile cycle. Complex vortex structures are found in the valve opening and closure phases. It is found in the present fully-coupled fluid/structure simulation that the vor- tices generated in the valve closure and regurgitant period encompass important mechanism which was overlooked 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 observed 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 simulation can yield correct description of a pulsatile transvalvular flow, in particular for analyzing the regurgitant flow in the diastolic phase.