Flow Characteristics in Elastic Arteries Using a Fluid-Structure Interaction Model

In this study the interaction of blood flow with arterial wall has been investigated using FSI (Fluid-Structure Interaction) modeling. Computer simulation of pulsatile blood flow was carried out on the basis of the time dependent axisymmetric Navier-Stokes equations for an incompressible Newtonian fluid flow. An elastic incompressible material with large deformation was considered for the arterial wall and momentum and continuity equations of elastodynamics have been solved. The specified boundary conditions for the Navier-Stokes equations were the pulsatile pressure waveforms of the brachial artery at inflow and outflow to the given pulse wave form of a cardiac cycle. Fluid and solid equations were solved with the ALE-based loose coupling method for FSI problems. Resultant flow, wall displacement, wall shear stress and wall circumferential strain waves, and their phase differences were determined. Stiffening of the arterial wall resulted in a significant decrease in the mean values of flow and wall shear stress and altered waveforms. A tenfold increase in wall stiffness caused 33% drop in flow and negative values of shear stress in 21% of the pressure pulse. For elastic moduli corresponding to wall displacements less than 1% the blood flow and wall shear stress were not sensitive to wall stiffness. Stress phase angle was altered by stiffening of the arterial wall. It was concluded that FSI modeling with pressure boundary conditions provides a proper evaluation of hemodynamic parameters that may determine endothelial injury.


INTRODUCTION
It has been well established that many cardiovascular diseases are closely associated with hemodynamic parameters and vessel wall mechanical characteristics. The arterial endothelial lining is exposed to both wall shear stress caused by pulsatile blood flow and circumferential stress caused by pulsating arterial pressure. These two stresses and their interaction are believed to play an important role in determining remodeling of the blood vessel and development of arterial diseases [1] . The human arterial system is formed by a network of distensible vessels that can be regarded as hollow tubes of variable diameters. Arteries deform under pressure fluctuation generated by the heart. Transmission of pressure waves through the blood is the result of the energy exchange between the blood and the vessel wall; therefore modeling of the time dependent arterial deformation or pulse propagation, is a fluid-solid interaction problem [2] .
Studies have modeled such phenomena. Witzing [3] , as a pioneer, solved the equations of motion for an inviscid fluid-filled elastic tube. Womersley [4] derived the coupled harmonic solution for a viscous fluid and wall motion in a deformable vessel, initiating research in this field. He solved the linearized Navier-Stokes equations for a thin walled isotropic elastic cylinder containing viscous Newtonian fluid.
Cox presented a complete review of this subject in 1969 [5] . This dealt with the extension of Womersley's theory to thick walled tubes, viscoelastic material behavior, anisotropic material behavior, initially stressed tubes and incorporation of non-linear terms in the Navier-Stokes equations. All computations of wall motion were geometrically and physically linear.
Development of numerical methods in moving boundaries was the first step in the application of fluidstructure interaction methods in cardiovascular systems. Reuderink [6] , Pietrabissa and Inzoli [7] used a decoupled method, in which for a given flow, the pressure wave was computed using linear wave propagation theory. Then, the equilibrium equations for the wall were solved with the previously computed pressure as natural boundary condition. The mesh for the fluid motion computation was deformed accordingly. Reuderink et al. deformed the outermost elements of their 3D carotid artery mesh [6] , whereas in their later publication on a 2D tube model, all elements were deformed [8] . In both cases (2D and 3D) the velocity of the nodal points with respect to stationary frame was neglected. Van de Vosse et al. computed the fluid motion using penalty function finite element method [9] .
Because a suitable wave propagation model has not been available due to material and/or geometrical complexities, decoupled methods have been used for simple models. Hilbert [10] , Steinman and Ethier [11] and Perktold and Rappithsch [12] solved the equations of motion with a weakly coupled method. At the entry of their models the flow was prescribed, at the ends the experimentally recorded pressures were used. In the computation for a periodic solution several flow cycles were applied. Steinman and Ethier [11] used a method similar to the one presented by Hilbert [10] . They solved the fluid motion equation with a penalty function method and computed the outlet pressure with a linear wave theory, rather than prescribing experimental data.
Perktold and Rappithsch [12] reported a 25% decrease in wall shear stress in an elastic model compared to a solid model. Rutten [13] used a loose coupling method to simulate pulsatile blood flow in 2D and 3D isotropic and anisotropic models. He used the strain energy density function for simulating nonlinear behavior of vessel wall. Qiu and Tarbell [1] used the free surface function of FIDAP®7.62 to simulate the moving walls of coronary artery based on experimental measurements in dogs. De Hart et al. [14] used a fictitious domain technique and loose coupling method to simulate 2D aortic valve motion in systolic phase. Their model was simplified with assumptions such as a solid vessel wall and neglected inertia of valve. They also used implicit coupling methods to simulate 3D aortic valve motion to the increase accuracy and stability of the algorithm [15] . Cebral et al. [16] simulated pulsatile flow of a real carotid artery obtained from MR images using an independent ring model for the elastic wall of the artery. They concluded that in rigid models, the regions with low shear stress (<10 dyn/cm 2 ) were smaller compare to those in elastic models. Lee and Xu [17] coupled ABAQUS and CFX4.2 to simulate blood flow in elastic arteries.
Most of the researchers have treated their models with pulsatile flow as a boundary condition. Since the pressure gradient is the cause of flow, application of pressure based boundary conditions might be more realistic in FSI modeling of arteries. In present study models of human brachial artery have been developed to evaluate effects of wall stiffness on mechanical parameters such as flow, wall shear stress and wall circumferential stress. Since critical values of those stresses contribute to endothelial damage the results may be applied in study of arterial pathology.
In addition to values of shear and circumferential stresses a quantitative parameter has been introduced to study effects of phase shift between wall shear stress caused by blood flow and wall circumferential stress caused by blood pressure. Brachial artery has been chosen due to the availability of pressure pulse data together with mechanical and geometrical parameters. Physiological pressure waves at the inlet and outlet of the model were used as boundary conditions. Results of the FSI algorithm were validated with published experimental data. The capability of the model to obtain hemodynamic parameters was discussed. The effects of stiffness of the arterial wall on flow rate and wall shear stress were studied and their application in endothelial pathology were discussed. Since aging is correlated with stiffening of the arterial wall the results are applicable in change of hemodynamic parameters such as pulsatile flow, pulsatile wall shear stress and circumferential stresses and their phase shift with age.

Mathematical representation:
The fluid model: To model the arterial system with the FSI model, the mechanical properties of blood and arterial wall should be studied. Blood is a complex suspension of different corpuscles. In large arteries, these corpuscles have dimensions much smaller than the diameter of the vessel and blood can be considered a homogeneous and incompressible fluid. Furthermore, it can be assumed to behave as a Newtonian fluid at high shear stress [1,16,18,19,20,21] .
Continuity and momentum equations of Newtonian fluid with negligible body forces are described as Eq. (1): is the deformation rate tensor. Character "I" stands for the unit matrix.
The wall model: The arterial wall as the solid part of the FSI model is comprised of three layers: intima, media and adventitia, each of them being constituted of different materials with different mechanical features such as collagen fibers, elastin fibers, smooth muscle cells, ground substances, and water [22,23] . Elastin in the media bears most of the pressure load at low strains. The collagen fiber network limits the radial deformability at higher blood pressures, and causes the steep rise in wall stiffness at higher strains, resulting in a material with nonlinear elasticity.
In the physiological range of blood pressure, the collagen fibers are stretched and the stress-strain relationship of the arterial wall may be considered as linear during the pressure pulse. Large deformation theory was used to model arterial wall deformation.
The continuity and momentum equations of elastodynamics for the wall are written as shown in Eq. The gradient in Eq. (2) is defined based on moving coordinates and is different from the gradient in Eq. (1). If solid is assumed as elastic and isotropic material obeying hook's law, Cauchy stress tensor can be used as Eq. (3): where C 4 is a fourth order elasticity tensor, : is the product sign and ) ) ( is the strain tensor for small deformations. Due to large deformations, the above strain is not applicable for the description of deformations that occur in pressurized arteries. The Green-Lagrange strain is objective and will be used instead. To use this, corresponding stress measure that relates stresses to the reference situation must be adopted. The resulting invariant stress tensor S is called the 2 nd Piola-Kirchhoff stress tensor [13] . Since S is invariant, it is a proper description for large deformation computations. To determine the stress in the material, Cauchy stress has to be computed afterwards as Eq. (4): is the deformation gradient tensor with respect to a reference configuration of the solid and J is the determinant of F. This tensor describes both rotation and displacement. Therefore, the Green-Lagrange strain tensor E is defined as Eq. (5): This tensor is reduced to the linear strain tensor for small deformations. Continuity equation for large deformation is defined as Eq. (6): Coupling conditions: In order to obtain a complete FSI model of coupled system, two coupling conditions on the fluid-solid interface are required [18,19,[24][25][26][27][28][29][30][31][32][33][34] . The first condition describes the fluid-solid interface as a Dirichlet boundary for the fluid i.e. the preset velocity values of the fluid must be equal to the structural nodal velocities describing by Eq. (7a): The second condition indicates that the fluid-solid interface is treated as a Neumann boundary for a solid, i.e. the fluid boundary tractions t f are surface loads for the structure and is given by Eq. (7b): where characters 0 and t show reference and changed configurations of fluid-solid interface, respectively. The negative sign refers to opposite directions of the normal vectors on the interface of fluid and structure with respect to a reference configuration. The formulation implies that the fluid boundary moves with the structure. The coupling conditions ensure conservation of mass, momentum and mechanical energy on the interface.

METHOD OF SOLUTION
Models of brachial artery were developed for analysis. Compared to other same level arteries wall stiffness of brachial artery is lower, hence effects of arterial stiffening on hemodynamic parameters can be studied with more details in a wide range of elastic moduli. Furthermore the geometrical conditions of brachial artery such as no branching between measurement sites make brachial artery a proper model for the analysis. The availability of pressure data on two close sites together with available data on mechanical and geometrical parameters were other reasons for choosing brachial artery.
Two different experimental pressure waveforms from previously published data were used as fluid boundary conditions [23] . These data were obtained from measurements along the human brachial artery at two points five centimeters apart. The time difference between the pressure waves was derived from Moens-Korteweg equation for wave velocity using the mechanical properties of the arterial wall (Fig. 1). Because of differing wave velocities for models with different elastic moduli, alterations in calculated pressure gradient wave were observed while the mean value for the pressure gradient pulse remained constant. Therefore stiffening of the arterial wall led to an altered pressure gradient wave. Mechanical properties and model dimensions used for the model are summarized in Table 1.   [19] and pressure gradient for differing wall stiffness (bottom) Axial displacement of the arterial wall was restricted at inlet and outlet boundaries. The tethering effect (axial pre-stretch) was not considered in the model. For the symmetrical model the adopted cylindrical coordinates also show the principal directions of stress tensor. Therefore considering boundary conditions the pre-stretch in axial direction does not affect the circumferential stress as a major model output.
Fluid and solid equations together with coupling conditions were solved by ALE-based loose coupling algorithm using FIDAP8.5 software in an axisymmetric model of the artery. Quadrilateral 4-node elements were used for fluid and solid mesh. The element numbers for the fluid and wall models were 2800 and 1120 respectively. Maximum aspect ratio and skewness were 3.9 and 0.01 respectively. That the results were independent of the computational mesh was verified. The governing finite element equations for both the solid and fluid were solved by segregated iterative algorithm. Each period was divided into 87 uniform time steps and convergence criteria were set to 0.001 for both velocity and displacement.
Using finite element method to reach a good convergence, the model was loaded to 80 mmHg at 1 s duration, then problem solved for five cardiac cycles with 0.863 s period corresponding to 70 bpm for a healthy human. Computational results showed that wall displacements reached the time-average absolute difference of 0.001 mm after the second period and velocities reached the time-average absolute difference of 0.001 m/s after the fourth period. These values were chosen as convergence criteria with an order of magnitude significantly smaller than the average outputs. The results were recorded from the fifth cycle at the mid cross section of the arterial model. The results included wall shear and circumferential stress waves, pulsatile flow, pressure pulse at different sites, and pulsatile radial displacement. Effects of arterial wall stiffness on the above parameters were studied. The values of the wall elastic moduli were chosen in the range from a healthy to an aged brachial artery.
The FSI algorithm was validated by comparing the results of the FSI model with a very high stiffness to a solid model without the FSI algorithm. Velocity convergence was reached after the eighth cycle in the model without the FSI algorithm. Similar results were obtained with negligible differences in velocity values at some stage of arterial pulse.
Results of FSI modeling were validated using FSI models of canine thoracic aorta with geometrical and material parameters similar to those of published experimental data [35] and obtaining similar results. The experimental data included the measurements of pressure waveform and radial displacements. The flow waveform obtained from experimental data with applications of long wavelength approximation theory.  [20,23,[36][37][38] . The FSI modeling results showed a good agreement with published data.

Effect of elastic modulus on hemodynamic parameters:
Flow rate: The effects of wall stiffness on flow rate are shown in Fig. 2. The stiffer the arterial wall, the faster pressure wave velocity. Since the pressure gradient wave is derived from the difference between pressure waves at proximal and distal sites and their time difference, it can be concluded that the pressure gradient waveform is highly affected by the wall stiffness ( Fig. 1(bottom)). The pressure gradient wave is dependent on the wave velocity which is in turn related to the wall stiffness. But this change is only for the shape of pressure gradient and not the average value (or the total value in an arterial pulse). This results in an altered flow waveform. Results show a significant drop in peak systolic flow by stiffening of the arterial wall and a further marked flow rate drop in the rest of the systolic phase. It can be seen that for a model with a highly stiffened wall (Y=4.5 MPa) flow becomes negative in late systole.
Effects of wall stiffening on mean blood flow rate is depicted in Fig. 3 and shows that there is a significant decrease in average flow rate by elevation of the wall stiffness: a two fold increase in wall stiffness results in a 17% decrease in average flow rate and a ten fold increase causes a 33% drop in average flow rate. Results show that an increase in stiffness to more than 4.5 MPa has no significant effect on flow due to the very low distensibility of the arterial wall. Published experimental data [19] point to a significant decrease in blood flow by aging which is correlated with wall stiffening.
The dependency of flow to the elastic modulus value is due to the fact that the stiffened arterial wall is less distensible with smaller radial displacements. Assuming the same geometry for depressurized normal and stiffened models, the elastic model expands less after loading compared to the stiff model resulting to limited flow for stiff models.
Since aging correlates with stiffening of the arterial wall, the results may be used to study the effects of aging. By stiffening of the arterial wall, the flow rate decreases markedly, and to maintain the flow level the pressure level needs to be evaluated in such a way that the artery expand to the previous level. This might be considered as one of mechanisms for the increase of mean arterial pressure by aging. The correlation of aging with arterial wall stiffening might lead to altered pressure gradient wave and backflow in late systole in aged arteries.   Figure 4 compares velocity profiles for models with two different wall elastic moduli at different stages of pressure pulse. As can be seen there are marked differences between velocity profiles in systolic phase (t/T=0.08, 0.14) while in diastolic phase (t/T=0.5) there are no significant differences. Backward flow near the arterial wall is observed for models with stiff wall in systolic phase (t/T=0.14) as already discussed. Backward velocity profiles result in negative wall shear stress in systolic phase for models with stiff wall.
Wall Shear Stress: Wall shear stress variation with stiffening of the arterial wall is shown in Fig. 5. The WSS wave shows a major fluctuation during systole. The systolic peaks are highly sensitive to the stiffness of the arterial wall which markedly decreases systolic peak values. A 10 fold increase of wall stiffness causes maximum WSS values to decrease 51%. This causes a negative minimum shear stress. Further increase of the elastic modulus has no significant effect on WSS values.  Since endothelial damage is correlated with WSS values, FSI modeling may provide useful information on WSS values and its relationship with arterial stiffening. The results may be applied in aging. Figure 5 indicates a significant decrease in the WSS minimum value with wall stiffening. It can be seen that for a model with a stiff wall the WSS remains negative in 21% of the pressure pulse duration, but for a distensible artery model the WSS remains positive throughout. While for distensible model with appropriate mechanical and geometrical parameters WSS remains in the positive domain, for a model with stiffened wall WSS fluctuates in negative and positive regions. Phase difference between biaxial loads on endothelial lining: The FSI method provides a suitable model to study the effect of the phase difference between wall shear stress and circumferential stress during the pressure pulse representing a biaxial pulsatile loading system, on the endothelial lining. In a study of endothelial pathology not only the critical values of shear and circumferential stresses are important but also the phase shift between them that might be a major determinant [1] .
The value of the circumferential stress is obtained by the quantification of the circumferential strain on the inner wall and the value of the endothelial cell elastic modulus. The CFD modeling of rigid vessel will not present the wall circumferential strain. However, FSI modeling provides a simultaneous solution for shear stress and circumferential strain and their phase shift.   To analyze the phase difference between these two waves a mathematical algorithm based on the discrete Fourier transform is done and the phase angle between two waves is calculated. Figure 6 shows the stress phase angle (SPA) between CS and WSS for different harmonic frequencies. The first harmonic frequency is 1.16 Hz corresponding to 70 bpm for an average healthy heart. The results are shown for models with different elastic moduli. Figure 7 shows the first harmonic values of SPA and the impedance phase angle between pressure and flow (IPA) for different wall elastic moduli. The results show a non-linear decrease in both SPA and IPA with stiffening of the arterial wall. Experimental results have shown that the more negative the SPA values the higher risk of arterial diseases [1]. As can be seen this trend occurs with stiffening of the arterial wall. Since mechanical and geometrical parameters affect the values of CS and WSS, FSI modeling may be used to study atherogenesis at different arterial sites.

CONCLUSION
FSI modeling is used to analyze the interaction between the elastic wall and blood flow in a model of brachial artery. The pressure based boundary conditions are utilized using physiological data and wave velocity. Evaluation of the model shows the capability to predict flow and wall mechanical parameters.
Results showed that FSI modeling with pressure boundary conditions provides a proper evaluation of hemodynamic shear stress and wall circumferential stress values which are major determinants of endothelial injury and atherogenesis. The availability of non-invasive pressure recording at different sites facilitates application of pressure boundary conditionbased FSI modeling. FSI modeling showed marked influence of wall stiffness on mechanical parameters based on the fact that wall stiffness effects markedly on the pressure gradient wave shape considering wave velocity. A tenfold increase in wall stiffness caused 33% drop in flow and negative values of shear stress in 21% of the pressure pulse resulting in fluctuation of shear stress in negative and positive regions, a high risk of endothelial.
For elastic moduli corresponding to wall displacements less than 1% the blood flow and wall shear stress were not sensitive to wall stiffness. The stress phase angle was also sensitive to wall stiffness. Due to the correlation of aging with arterial wall stiffening (arteriosclerosis) the results can be extended in study of aging effects. By aging arterial wall stiffens resulting in a decrease of blood flow rate. Due to metabolic requirements of organs a sustainable blood flow in the physiological range is needed. This is achieved by control mechanisms such as elevation of mean arterial pressure and arterial wall remodeling.