Cardiovascular diseases, mainly due to atherosclerosis, are causes with the highest percentage leading to death worldwide.
Curiosity of finding cures for these diseases has enthused mathematicians to study them from their mathematical viewpoint. Consequently, many computational techniques and models have been developed. These tools aim at describing the blood flow and are effictive in studying the response of the arterial wall under certain conditions, the characteristics of the blood components
in addition to those of the heart [3, 16, 17, 21, 22, 24].
Indeed, mathematical models with numerical analysis and simulations play an important role in providing knowledge and insights that are
The rheological behavior of blood is captured by deriving constitutive models that constitute a constructive tool in the diagnosis of the pathologies, investigating appropriate remedies and proposing preventive therapies [1, 14].
Recently, the lumen-wall modeling has been adopted using fluid-structure interaction (FSI) model. In the FSI, the behavior of the blood and the arterial wall are taken into consideration, so that one is capable of representing them by their appropriate dynamics and models. A most commonly used method when dealing with FSI systems is the Arbitrary Lagrangian-Eulerain (ALE) method  that is effective when combining the fluid formulation in the Eulerian description and the structure formulation in the Lagrangian description.
An overview of FSI in biomedical applications has been considered in the book . In particular, modeling of cardiovascular diseases has been highlighted in [9, 18, 25].
In lumen-wall modeling some difficulties encounter
due to the complexity of the arterial wall formed of several layers, each with its own unique mechanics and thickness.
Assuming that the arterial wall is negligibly thin, or the ratio of the arterial wall thickness to the aretry raduis is small, reduced shell or membrane models have been employed [4, 13].
Introduction of computational model with FSI in order to investigate the wall shear stresses, blood flow field and recirculation zones in stenosed arteries have been studied in  where the blood is considered to be an incompressible Newtonian fluid, whereas in the case of a compressible non-Newtonian fluid these factors have been analyzed in .
Currently, numerous computational models are simple when describing the cardiovascular diseases and the related processes such as inflammations, coagulation, plaque growth, clot formation, etc.. Indeed, they are managed to capture only some essential features of the processes that take place in the cardiovascular system. Further, these models neglect some of the blood components with their characteristics as well as the arterial wall layers and their own mechanics.
Consequently, more suitable models are needed through which the physiological parameters associated to the blood, atherosclerosis and clots must be investigated clinically.
In addition, the arterial wall must be considered as a multi-component structure taking into account the effects of the plaque growth and the clot formation on their mechanism. Further, timescales of the biological phenomena, the pulse duration and the time of the plaque growth must be analyzed.
In the present work, the interaction between the blood modeled by a modified Carreau’s model and the hyperelastic incompressible arterial wall has been considered. In a first step, we introduce the FSI system which is composed of the incompressible Navier-Stokes equations representing the blood flow dynamics and the quasi-static equilibrium equations describing the elastic large deformation of the arterial wall. In addition, coupling conditions that ensure a global energy balance of the FSI system have been imposed on the common interface. Variational formulation has been presented and its discrete formulation has been derived. In fact, the Navier-Stokes equations have been semi-discretized in time, while, the nonlinear material equilibrium equations have been solved using the Newton-Raphson method. Numerical simulations have been performed using FreeFem++ . A deep analysis has been made for better understanding of the behavior of the blood flow, the blood viscosity, the maximum shear stress and the recirculation zones. Based on the numerical results, a location where the blood is thought to accumulate and solidify has been identified. Finally, the factors affecting this zone have been investigated, in particular, this zone is subjected to forces exerted by the artery wall and the blood flow. Consequently, numerical simulations for the deformation of this zone have been performed in order to understand its rupture and thus the release of a clot that will lead to the occlusion of small arterioles.
2. The Fluid-Structure Interaction Problem
The total domain representing the artery in the actual configuration at time is composed of two sub-domains and representing the lumen of the artery and the arterial wall, respectively.
We denote by the reference configuration of the structure of density . Its deformation is described by the displacement field that satisfies the quasi-static incompressible equilibrium equations. The evolution of the structure domain can also be given by the deformation map defined in terms of the displacement as . Its deformation gradient
which is a second order tensor is given by. Its associated Jacobian is .
On the other hand, we describe the blood flow dynamics by the incompressible Navier-Stokes equations on the sub-domain . We denote by
the velocity of the blood and its pressure, respectively. Further, the blood is assumed to be an incompressible fluid with a constant density .
The sub-domain of moving boundaries evolves from some reference configuration according to an ALE map given by
that is, .
The ALE map is considered to be an extension of the displacement of the interface , that is to say
The operator stands for an extension of the displacement of the boundary . Possible extensions can be found in [19, Section 5.3], [5, Chapter 2] (harmonic, biharmonic, wislow, etc.). In particular we consider the harmonic extension as we will see in Subsection 2.1. The deformation gradient associated to is defined by where the symbol indicates the gradient with respect to the variable .
Its Jacobian is .
Here and throughout the context, denotes the displacement of the domain which we set to be . Formulating the Navier-Stokes equations in the ALE frame results a new variable that describes the velocity of the domain . It is related to the displacement by the relation . It is worth to point out that . One must distinguish between the physical velocity of the particles and the velocity of the fluid domain .
In what follows, we refer to the elements in the reference configuration by ” ”. In fact the velocity and the pressure of the blood are given on the reference configuration by
The Cauchy stress tensor is expressed in terms of the strain tensor as
where represents the blood viscosity that will be detailed in the sequel. In the reference configuration , the stress tensor is given by
The arterial wall is assumed to be a hyperelastic material then it is characterized by the existence of an energy density function such that the first Piola-Kirchhoff stress tensor . Further, due to the incompressible behavior of the material its Piola-Kirchhoff stress tensor is modified to the form
The variable , called the hydrostatic pressure, plays the role of the Lagrange multiplier associated to the incompressibility condition .
On the fluid domain , a volumetric force is applied. Moreover, a velocity is enforced on the inlet of the artery . On the contrary, a free-exit condition given by is enforced on the outlet .
On the other hand, a volumetric force is applied on the structure domain which is assumed to be fixed on the boundary , that is to say, on .
On the interface , surface forces and are exerted from the fluid domain and the structure domain, respectively.
The FSI model describing the blood-wall interaction is obtained by the coupling
between the incompressible Navier-Stokes equations which are formulated in the ALE frame and the quasi-static incompressible elasticity equations formulated in the Lagrangian frame on the reference configuration . The FSI system is
From Expression (2) we get that the ALE map and the structure deformation coincide on the interface , that is to say,
Due to the incompressibility condition the -th component of is
where is the Kronecker symbol. The shear stress components are and , whereas and are the normal stress components.
In a two dimensional space the maximum shear stress- an effective parameter in studying the forces exerted on a fluid- is given by the expression 
The variational formulation associated to System (5) is
, , ,
The coupling conditions on the interface are given in the strong form as
The spaces and are respectively
2.1. The Discrete Variational Formulation of the FSI Problem
The variational formulation (7)-(10) of the FSI stands for the incompressible homogeneous Navier-Stokes equations coupled with the quasi-static incompressible equilibrium equations. We assume that no external forces are exerted on neither the fluid domain nor the structure domain, i.e, and . Consider a time step and finite element partitions and for the fluid and the solid sub-domains respectively, of a maximum diameter denoted by . Our aim is to approximate the solution at time , for , in the finite element spaces. The approximation of the solution at time is denoted by .
Semi-Discretization in Time of the Fluid Sub-Problem
In order to guarantee the existence and uniqueness of the solution of the discrete fluid sub-problem when performing the numerical simulations, we use the penalty method . This method consists of replacing the natural weak formulation by a regular one by adding a term multiplied by a sufficiently small parameter . Indeed, writing the modified formulation in a matrix form results a positive definite matrix, which assures the existence and the uniqueness of the solution of the discrete sub-problem. The weak formulation associated to the Navier-Stokes equations obtained upon adding a negligible parameter is then semi-discretized in time, that is, the convective term and the viscosity are considered at the instant , whereas other terms are considered at time . The discrete formulation reads
where the non-linear convective term can be approximated by [12, Section 9.5, p. 267]
Notice that, since the strain rate tensor is symmetric, then we have which gives
Newton-Raphson Method for the Structure Sub-Problem
Regarding the structure sub-problem, at the time iteration we will solve the non-linear problem (9) using Newton-Raphson method. The variational formulation corresponding to the structure sub-problem at the iteration is
The method depends on linearizing the structure sub-problem (13) with respect to the unknowns and .
We start initialization by considering a suitable choice of the initial values (). In particular, we link the iterations of the Newton-Raphson method with the time iteration by considering and . Then, we solve iteratively the obtained system corresponding to the Newton-Raphson method until its solution converges to a solution of the non-linear System (9). To ensure the existence of the solution of the structure problem (13) we use the penalty method by modifying System (13) through adding
the penalized term with . For simplicity of notation, in what follows we omit the subscript of the deformation , that is, we write .
We proceed to derive the formulation of the structure sub-problem corresponding to the Newton-Raphson method. Let us define the following space
Given , a tolerance and
we construct iteratively the two sequences and by solving for the following system:
Set . Repeat: for , while , find (,) in satisfying
for all .
Set and .
When the condition is fulfilled then convergence of the Newton-Raphson method is achieved. Thus, the solution of the structure sub-problem (13) is given by , for the last value of for which the Newton-Raphson method converges.
Space Discretization of the FSI Problem
Space discretization of the variational formulation is carried out using the finite element method (FEM) . We consider the two finite element spaces associated to the fluid weak formulation
and those associated to the structure weak formulation
where , , and are finite dimensional subspaces.
The functional spaces associated to the velocity and displacement fields are considered to be , whereas those associated to the pressures (fluid and hydrostatic) are considered to be .
In what follows, all terms are discretized in space as mentioned above, so that the approximation of solution in finite element spaces is verifying (11) and (14).
Finally, the discrete variational formulation reads:
Given and a tolerance , find such that
where the non-linear convective term
is approximated by the Expression (12).
The coupling conditions on are
Fix . Set and . Repeat: for , find (,) in satisfying