Analysis of Blood Turbulent Flow in Carotid Artery Including the Effects of Mural Thrombosis Using Finite Element Modeling

Arterial thrombosis is an extremely significant health problem. As a result of numerous factors involved in such problem, describing the role of hemodynamics in thrombogenesis has been asserted to be one of the most demanding and complicated challenges in biomechanics. An axisymmetric model considering fluid-structure interactions (FSI) was introduced and numerically solved for an artery with a thrombosis to perform flow and stress-strain analysis and investigate the probability of thrombus ruptures leading to embolization. Three models with different thrombus heights were considered and the Navier-Stokes equations were solved for the blood flow as the fluid domain. Results indicated that there are recirculation regions after thrombus bulk, which are susceptible to rethrombosis and stenosis. It was also shown that when the thrombus height increases, the shear stress magnitude on FSI boundary increases and the area near the thrombus peak is too susceptible to rupture. Besides, stress-strain distribution analysis demonstrated that by increasing the thrombus height, the region with high shear stress on the wall declines while the shear stress magnitude of the region under the peak increases up to 4 times. When the thrombus height is low enough (34% of artery diameter), its deformation is larger at the peak and a large area of its downstream side. However, by increasing the thrombus height, there are two sites of large deformation in thrombus at the peak and a small area at the leading edge (in compression site) of thrombus. These regions are vulnerable because of rupture probability.


INTRODUCTION
Arterial thrombosis can block blood flow to heart and brain tissues and cause a heart attack or stroke. Clinical studies confirm the link between thrombosis and atherosclerosis [1] . Thrombus superimposed on ruptured atherosclerotic plaque is commonly found in autopsy studies of heart disease [2,3,4] . Thrombosis is also associated with carotid artery plaque rupture in stroke and transient ischemic attack [5,6,7] . Since the relationship between thrombosis and fluid mechanics is complex and numerous factors are involved in this issue, arterial thrombosis is the subject of intense study and speculation in biomechanics. The process of thrombosis may be affected by a series of rheological and fluid dynamic parameters, including high rates of shear, areas of flow stagnation or recirculation and turbulence. In fact, Thrombosis is fundamentally linked to hemodynamics because blood transports cells and proteins to the thrombus and applies stress that may disrupt the thrombus. The convection of flow aggregates can be locally enhanced or diminished in arterial stenosis through localized regions of high shear, flow separation, recirculation, and reattachment [8] . The role of hemodynamics in thrombogenesis has been investigated on natural and artificial surfaces [9,10,11,12] . In two other studies by the same author, a mechanistic approach based on simple classical fluid dynamic concepts has been adopted to quantify the hemodynamic forces acting on mural thrombi under both steady and pulsatile flow conditions [13,14] . Arterial thrombosis is also investigated via numerical studies. However, crucial limiting assumptions such as Reynolds number much less than that accrue in cardiovascular system diseases as well as deviates from realistic boundary condition because of considering a very small numerical solution region are made in these studies [10] . Besides, several investigations have considered shear stress on the artery wall but there is no study that considers shear stress in flow domain where the most platelet activating process takes place.
Accurate modeling of arterial thrombosis has not been proposed yet and due to the importance of such modeling to get further insight into clinical aspects of thrombosis such modeling is extremely demanded. In this paper, for the first time, an axisymmetric model considering fluid-structure interactions (FSI) is introduced and numerically solved for common carotid artery with a thrombosis to perform flow and stress/strain analysis and investigate the probability of thrombus rupture which leads to embolization.

MATERIALS AND METHODS
Solid and Fluid Models:In order to investigate the effects of thrombosis on the artery hemodinamics and stress distribution of the artery wall, the model geometry is considered (Fig. 1). Due to the high shear rate, thrombosis in arteries is formed from platelets and does not contain red blood cells [15] . Therefore, thrombus is assumed to have the same properties as platelets do and to be linear elastic, isotropic, incompressible and homogenous material. Young's modulus of 3 kPa and Poison's ratio of 0.48 is considered for thrombus [16] . Equilibrium equations and boundary conditions for thrombus are: ij are displacement and stress tensors for solid and fluid, respectively. Comma denotes the derivative with respect to the corresponding coordinate. The flow is assumed to be laminar, Newtonian, viscous and incompressible. The incompressible Navier-Stokes equations in Arbitrary Lagrangian-Eulerian (ALE) formulation are used as governing equations [18] . This yields the following: where u is flow velocity, g u is mesh velocity, p is pressure, in p is inlet pressures, out p is outlet pressures and ρ is fluid density. We assumed three models with different thrombus heights (1.25, 1.75 and 1.8 mm respectively) and the fixed outlet pressure of 20 mmHg [17] . Fig1. Geometry and dimensions of artery and thrombus (d varies in three models 1.25, 1.75 and 1.8 mm respectively and the ureter diameter is assumed 7 mm).
An axisymmetric model with fluid-structure interactions (FSI) is introduced and solved using the ADINA code (ADINA TM , version 8.3, Automatic Dynamic Incremental Nonlinear Analysis, Watertown, MA). For thrombus, as shown in Fig. 2a, a rule-based mesh of 4470 axisymmetric rectangular solid elements of the plane strain type with 9 nodes in each element is used. It should be mentioned that are arranged so that the aspect ratios are less in elements which are under compression to prepare a more compression strength. This is done because in an ALE mesh, when deformation in compression direction is imposed to fluid elements, overlapping of solid and fluid elements will occur [18] . Besides, due to large deformation of the compression site of the thrombus, larger elements are considered in this site. The mesh is also arranged to avoid distortion in solid elements due to tension or compression forces imposed from fluid domain on solid elements. For fluid domain, as shown in Fig. 2b, a rule-based mesh of 83030 axisymmetric triangle solid elements of the planar fluid type with 3 nodes in each element is used. ALE formulation and leader and follower points are used to keep the mesh quality in fluid elements during the solution time. It should be mentioned that considering the selection of the leader and follower point's method and the deformation of the fluid d elements due to displacement of the solid body, a mesh with triangle elements will lead to a good convergence.
To impose the loads, we considered 30 steps of 1s and increasing the inlet velocity and outlet pressure with a ramp function. In order to reach a steady state solution with the established inlet velocity and, more importantly, the correct pressure drop along the artery, 10 time steps of .001s with the fixed ultimate load are considered.
Solution Methods Used in the Software:In order to solve the governing equations of the fluid, Upwind and finite element methods and Sparse solver (which uses the standard Galerkin procedure for incompressible fluid flow) are used. Both the fluid and solid elements satisfy the inf-sup (infimum supremum) or Babus ka-Brezzi condition for stability and optimality in nearly or exactly incompressible analysis [18] . Newtonian iteration method and FSI conditions are used in solving process. The finite element equations for both the structure and the fluid were solved by the Newton-Raphson iterative method. Convergence is considered achieved if: (9) where f denotes the variables being solved (flow velocity, pressure and wall displacement), i is iteration a. b. u are displacement and acceleration vectors, respectively and f is the external force vector. In the implicit Newmark algorithm, the following assumptions are made: α α (12) where α and δ are parameters to control the precision and stability of integration and are considered to have a value of α =0.25(0.5+ δ ) 2 and δ = 0.54 which are proved to be suitable for the problem of the present paper by many evaluation examples. The solution obtained was fully coupled in the sense that the kinematical conditions of displacement, velocity, and acceleration continuity across the no-slip fluid-solid (artery-blood) interface are satisfied at all the time during the analysis.

RESULTS AND DISCUSSION
This study is regarded as an initial step toward the modeling of the effects of thrombosis on arteries considering fluid-structure interactions (FSI) The main objective of this study is to investigate the effects of a thrombus on hemodinamics and stress distribution of the artery considering the model assumptions (FSI, load and boundary condition), input parameters (severity, nonsymmetrical and imposed pressure and velocity) and output variables (flow characteristics). In order to investigate the probability of thrombus rupture which leads to embolization, blood velocity, shear stress and pressure gradient distribution as well as stress and strain distribution in the thrombus are discussed in the following.

Blood Velocity, Shear Stress and Pressure Gradient
Distribution: Fig. 3 shows the velocity distribution of blood in artery with different thrombus heights. As it can be seen, the maximum magnitude of velocity occurs after the thrombus region and has a value of about 2 times of the inlet velocity. Results show jet flow along the axis of the artery and recirculation regions near the wall, after the thrombus. Fig. 4 shows recirculation regions and flow reattachment points. These regions are susceptible to rethrombosis and stenosis [19] . Figure 5 shows Turbulence-k distribution in artery in three models which tell us about the flow kinetic energy and as we see from the figure, by increasing the thrombus height, the high magnitude region is reduced near its peak. Fig. 6 and Fig. 7 show fluid wall shear stress as a function of longitudinal position in two cases of 35% and 51% area reduction in site of thrombus respectively. Results show that when the thrombus height is low (1.25mm) we have three peaks of shear stress on the artery wall. The first region is before the thrombus bulk that is a susceptible area for forming the plaque, the maximum is near the thrombus peak and the third region which is much lower than the others is directly after the thrombus bulk on the wall. By enlarging the thrombus in length and height, the high shear stress before thrombus bulk vanished. Furthermore, the shear stress maximum value near the thrombus peak increased so that the maximum value of shear stress in third model (h=1.8mm) was about 3.6 times the value of the first model (h=1.25mm). Pressure gradient variation as a function of longitudinal position for the case of thrombus 51% area reduction (third model) is shown in Fig. 8. It is realized that pressure gradient variation on FSI boundary, for all the cases is higher (about 5 to 8 times) than that on the artery symmetry line. It should be mentioned that by obtaining the shear rates according to shear stress, our results are in the range of shear rates introduced in a theoretical and experimental study of a thrombus problem with conditions close to that considered in our study and this may validate the correction of the results of our study [14] . Its also worthwhile to mention that, using the embolization probability map as a function of wall shear rate and thrombus height introduced in the mentioned theoretical and experimental reference, it can be concluded that in lower magnitudes of thrombus heights (h=1.25mm), regions near the top of the thrombus (with shear rate of the value more than 1200 s -1 ) are exposed to high probability of rupture and embolization, but the regions of the thrombus with height less than the thrombus half height and at the right hand side of thrombus bulk (with shear rate of the value more than 100 s -1 and less than 1200 s -1 ), have average embolization probability (about 25% to 75%). In addition, when the thrombus height increases (h=1.8mm) shear rate on the top of the thrombus have a value of more than 20000 s -1 and all regions of thrombus are exposed to high probability of embolization.

Stress and Strain Distribution in Thrombus:
Thrombus stress analysis is important for a better understanding of its rupture process. Fig.9 shows the shear stress distribution in thrombus bulk in three models with different thrombus heights. As we see from the figure, when the thrombus height is low, the shear stress has high magnitudes at the middle of the thrombus bulk which placed on artery wall and a region a little under the thrombus peak. By increasing the thrombus height, the region with high shear stress on the wall reduced while the shear stress magnitude of the other region increased. The beginning and ending parts of thrombus have lower beginning and ending parts of thrombus have lower magnitudes of shear stress than the peak side of thrombus bulk which is susceptible to be ruptured. Figure 10 shows Strain distribution in thrombus bulk in three models. When the thrombus height is low, its deformation is larger at the peak and a large area downstream side of it but by increasing the thrombus height, there are two sites of high deformation in thrombus that are at the peak and a small area at the ending (in site of compression) of thrombus. These regions are in danger because of rupture probability.

CONCLUSION
In this study, an axisymmetric model considering fluid-structure interactions (FSI) is introduced and numerically solved for an artery with a thrombosis. Effects of a thrombus on hemodinamics and stress distribution of the artery and the probability of thrombus rupture which leads to embolization are investigated. Our results indicate that thrombus is very sensitive to hemodynamic variations. This study is regarded as an initial step toward the modeling of the effects of thrombosis on arteries considering fluid-structure interactions (FSI) and may help future efforts to be directed in this regard. For a more accurate analyzing of hemodynamic conditions in an artery with thrombus and to get further insight into clinical aspects of thrombosis and embolization, 3-D modeling with proper model assumptions based on real patient data and accurate morphology and material properties of thrombus as well as considering elastic properties of artery wall and pulsatile flow condition are essential to reach more reliable results. Such an accurate simulation can be extremely useful to determine critical conditions and to identify thrombus prone to rupture.