Fluid-Structure Interaction in Coronary Stents: A Discrete Multiphysics Approach

: Stenting is a common method for treating atherosclerosis. A metal or polymer stent is deployed to open the stenosed artery or vein. After the stent is deployed, the blood ﬂow dynamics inﬂuence the mechanics by compressing and expanding the structure. If the stent does not respond properly to the resulting stress, vascular wall injury or re-stenosis can occur. In this work, a Discrete Multiphysics modelling approach is used to study the mechanical deformation of the coronary stent and its relationship with the blood ﬂow dynamics. The major parameters responsible for deforming the stent are sorted in terms of dimensionless numbers and a relationship between the elastic forces in the stent and pressure forces in the ﬂuid is established. The blood ﬂow and the stiffness of the stent material contribute signiﬁcantly to the stent deformation and affect its rate of deformation. The stress distribution in the stent is not uniform with the higher stresses occurring at the nodes of the structure. From the relationship (correlation) between the elastic force and the pressure force, depending on the type of material used for the stent, the model can be used to predict whether the stent is at risk of fracture or not after deployment.


Introduction
Atherosclerosis is a condition where arteries become clogged with fatty substances called plaque. The plaque is deposited on the arterial wall, and this leads to the narrowing of the artery and subsequently obstruction of the blood flow known as stenosis. This obstruction hinders the smooth transportation of blood through these arteries and consequently poses a serious health problem. When atherosclerosis affects an artery that transports oxygenated blood to the heart, it is called coronary artery disease. Coronary artery disease is the most common heart disease that becomes the leading cause of death globally. Worldwide, it is associated with 17.8 million death annually [1]. The healthcare service for coronary artery disease poses a serious economic burden even on the developed countries, costing about 200 dollars annually in the United States.
The obstruction of the flow line (stenosis) also alters the blood flow regime and causes a deviation from laminar to turbulence, or even transitional flow [2,3], a situation that signifies severely disturbed flow. Studies were carried out on the types of plaque and its morphologies as well as the flow type and its consequence [4][5][6][7] that occurred in human arteries. Although the disease is deadly, it is preventable. Therefore, it is paramount to manage or prevent it in order to restore normal blood flow in the affected artery.
One of the ways of managing coronary artery disease is restoring normal blood flow or revascularization in a patient with a severe condition using the percutaneous coronary intervention (with stent) [8]. Coronary stents are tubular scaffolds that are deployed to recover the shrinking size of a diseased (narrowed) arterial segment [9] and stenting is a primary treatment of a stenosed artery that hinders smooth blood flow [10].
The stents used in clinical practice come in differential geometry and design which implies varying stress distribution within the local hemodynamic environment as well as on the plaque and artery [11,12]. The stent structure also induces different levels of Wall Shear Stress (WSS) on the wall of the artery [9,13]. Many cases of stent failure due to unbearable stress were reported and therefore, a careful study on how these stresses are distributed is needed. In fact, stent fracture or failure often occurs after stent implantation, and it can be avoidable if the mechanical property and the performance of the material are predicted. An ideal stent should provide good arterial support after expansion by having high radial strength. It should also cause minimal injury to the artery when expanded and should have high flexibility for easy maneuvering during insertion [14,15].
Studies on the different stent designs and how they affect their mechanical performance were reported [12,13]. Stent deformation and fracture after implantation were also investigated [16][17][18][19][20]. Moreover, numerical modelling and simulations were also used in studying coronary stent and stent implantation. For instance, Di Venuta et al., 2017 carried out a numerical simulation on a failed coronary stent implant on the degree of residual stenosis and discovered that the wall shear stress increases monotonically, but not linearly with the degree of residual stenosis [8]. Simulation of hemodynamics in a stented coronary artery and for in-stent restenosis was performed by [21,22].
With a few exceptions [23][24][25], the mechanical properties of the stent and the blood fluid dynamics around the stent were studied separately. In this work, we propose a single Fluid-Structure Interaction (FSI) model that calculates the stress on the stent produced by the pulsatile flow around the stent; both the stent mechanics and the blood hydrodynamics are calculated at the same time. The model is based on the Discrete Multiphysics (DMP) framework [26], which has been used in a variety of FSI problems in biological systems such as the intestine [27], aortic valve [28,29], the lungs [28], deep venous valves [30,31]. In this study, therefore, we use the DMP framework to develop an FSI 3D coronary stent model coupled with the blood hydrodynamics and analyse the mechanical deformations produced by the flow hydrodynamics.

Discrete Multiphysics
Discrete Multiphysics framework combines together particle-based techniques such as Smooth Particle Hydrodynamics (SPH) [32,33], Discrete Element Method (DEM) [34,35], Lattice Spring Model (LSM) [36,37], PeriDynamics (PD) [38], and even Artificial Neural Networks (ANN) [39,40]. In this case, the model couples SPH and LSM. The computational domain is divided into the liquid domain and the solid domain. The liquid domain represents the blood, and it is modelled with SPH particles; the solid domain represents the stent and the arterial walls and it is modelled with LSM particles. Details on SPH theory can be found in [41], and of LSM in [42,43]. Here a brief introduction of the equations used in SPH and LSM is provided.

Smooth Particle Hydrodynamics (SPH)
This section provides a basic introduction to SPH; additional details can be found in [41,44]. The general idea of SPH is to approximate a partial differential equation over a group of movable computational particles that are not connected over a grid or a mesh [37]. Newton's second law is integrated to give an approximate motion of the particles characterized by their own properties such as mass, velocity, pressure, and density expressed by the fundamental identity: where f (r) is a generic function defined over the volume, r is the position where the property is measured, and δ(r) is the delta function which is approximated by a smoothing (Kernel) function W over a characteristic with h (smoothing length). In this study, we use the so-called Lucy Kernel [41]. This approximation gives rise to which can be discretised over a series of particles of mass m = ρ(r)dr obtaining where m i and ρ i are the mass and density of the ith particles, and i ranges over all particles within the smoothing Kernel. Using this approximation, the Navier-Stoke equation can be discretised over a series of particles to obtain: where v is the particle velocity, t the time, m is the mass, ρ the density, and P the pressure associated with particles i and j. The term f i is the volumetric body force acting on the fluid and Π i,j introduces the viscous force as defined by [45]. An equation of state is required to relate pressure and density. In this paper, Tait's equation of state is used: where c 0 and ρ 0 are a reference sound speed and density. To ensure weak compressibility, c 0 is chosen to be at least 10 times larger than the highest fluid velocity.

Lattice Spring Model (LSM)
Elastic objects can be simulated using lattice spring models. As already discussed in [29], the main element of this model is composed of a mass point and linear spring which exerts forces at the nodes connected by a linear spring and placed on a lattice. Any material point of the body can be referred to by its position vector r = (x, y, z) [46]; when the body undergoes deformation its position changes and the displacement is related to the applied force as: where F is the force, l 0 is the initial distance between two particles, l is the instantaneous distance, and k the spring constant (or Hookean constant). According to [42], in a regular cubic lattice structure, the spring constant is related to the bulk modulus of the material by and where K is the bulk modulus, E the young modulus, l 0 is the initial particle distance and k the spring constant. From Equations (7) and (8) the spring constant is then related to the Young modulus of the material by ChemEngineering 2021, 5, 60 4 of 11

Model and Geometry
A three-dimensional stent model including blood flow hydrodynamics and stent mechanics is developed. The model simulates the blood dynamics in a 3D channel similar to a coronary artery with a 1.5 × 10 −3 m internal radius, including a PS-shape stent of 4 struts in the x-direction and 4 struts in the z-direction (circumference). The stent has a thickness of 100 µm, and 7.5 mm length, the size that is within the range of the stent used in clinical practice [47]. We choose a PS-Shaped because it performs better compare to most commercially shaped stents as reported by [25].
The geometry was created using CAD. From the geometry, we used MATLAB script to generate the coordinate of the computational particles as the points are created with MATLAB (details can be found in [29]). The script also generates a LAMMPS data file for the simulations and the simulations were run with LAMMPS, an open-source software [48]. The three-dimensional model consists of 1,862,804 particles: 1,609,452 particles for the fluid, 46,336 particles for the stent and 207,016 particles for the arterial wall. The fluid has a density (ρ) of 1056 kg m −3 and viscosity (µ) 0.0035 Pa·s. Figure 1 shows the section geometry of the stent within and outside the arterial wall. Local acceleration term g 0 was included to force the fluid to flow at a particular velocity. The inclusion of the local velocity is due to the unsteady or pulsatile flow existing in the cardiovascular system [49]. Womersley parameter α, which is the ratio of unsteady force to viscous force, was used in the model to induce the velocity profile of the flow. The particle spacing is 3.33 × 10 −5 m. The optimal spacing value is obtained after several simulations with different particle spacings to make sure the results are independent of the particle resolution. Stress and deformation and other postprocessing calculations were done with the visualization software OVITO [50]. The arterial wall is assumed to be rigid with a no-slip condition [25] whereas the stent is elastic. The no-slip and no penetration boundary condition is imposed at the interface of the solid-liquid interaction as discussed in our previous work [51]. Figure 1a,b shows a section of the solid geometry which includes the stent and the wall and the complete geometry of the stent respectively.
A three-dimensional stent model including blood flow hydrodynamics and stent mechanics is developed. The model simulates the blood dynamics in a 3D channel similar to a coronary artery with a 1.5 × 10 m internal radius, including a PS-shape stent of 4 struts in the x-direction and 4 struts in the z-direction (circumference). The stent has a thickness of 100 μm, and 7.5 mm length, the size that is within the range of the stent used in clinical practice [47]. We choose a PS-Shaped because it performs better compare to most commercially shaped stents as reported by [25].
The geometry was created using CAD. From the geometry, we used MATLAB script to generate the coordinate of the computational particles as the points are created with MATLAB (details can be found in [29]). The script also generates a LAMMPS data file for the simulations and the simulations were run with LAMMPS, an open-source software [48]. The three-dimensional model consists of 1,862,804 particles: 1,609,452 particles for the fluid, 46,336 particles for the stent and 207,016 particles for the arterial wall. The fluid has a density ( ) of 1056 kg m −3 and viscosity ( ) 0.0035 Pa•s. Figure 1 shows the section geometry of the stent within and outside the arterial wall. Local acceleration term was included to force the fluid to flow at a particular velocity. The inclusion of the local velocity is due to the unsteady or pulsatile flow existing in the cardiovascular system [49]. Womersley parameter α, which is the ratio of unsteady force to viscous force, was used in the model to induce the velocity profile of the flow. The particle spacing is 3.33 × 10 m. The optimal spacing value is obtained after several simulations with different particle spacings to make sure the results are independent of the particle resolution. Stress and deformation and other postprocessing calculations were done with the visualization software OVITO [50]. The arterial wall is assumed to be rigid with a no-slip condition [25] whereas the stent is elastic. The no-slip and no penetration boundary condition is imposed at the interface of the solid-liquid interaction as discussed in our previous work [51]. Figure 1a,b shows a section of the solid geometry which includes the stent and the wall and the complete geometry of the stent respectively.
The flow is being driven by a sinusoidal (pulsatile) acceleration ( ) of the flow in the axial direction to simulate a heartbeat of 60/min with a period ( ) of 1 s; the same approach was used in a previous publication [29], where the reader can find more details.
Eighteen (18) sets of simulations were run with a combination of 1 − 3 and 1 − 6 (see Section 4 for reasons). The University of Birmingham Bluebear (super-computer) was used for the computation (simulation) where ninety (90) processors were assigned for 72 h to run each simulation. For each simulation, a dump file (result file) of 9-10 GB of memory is obtained.  The flow is being driven by a sinusoidal (pulsatile) acceleration (G) of the flow in the axial direction to simulate a heartbeat of 60/min with a period (T) of 1 s; the same approach was used in a previous publication [29], where the reader can find more details.
Eighteen (18) sets of simulations were run with a combination of v1 − v3 and k1 − k6 (see Section 4 for reasons). The University of Birmingham Bluebear (super-computer) was used for the computation (simulation) where ninety (90) processors were assigned for 72 h to run each simulation. For each simulation, a dump file (result file) of 9-10 GB of memory is obtained.
To better understand the hydrodynamics of the fluid and to be able to access the deformation of the stent, the discussion is carried out in terms of dimensionless numbers. According to the Buckingham π theorem, a physically meaningful equation involving n physical variables can be rewritten in terms of a set of p = n − k dimensionless parameters Π 1 , Π 2 , . . . , Π p , where k is the number of physical dimensions involved. In this case, the physiochemical properties of blood are constant, and the geometry is fixed. Therefore, we want to express the resulting stress on the stent as a function f of the type In this case, the dynamic pressure can be written in terms of fluid average velocity v and density r, Moreover, k [kg s −2 ] and E are related in Equation (9). Therefore, assuming that the lattice spacing is fixed, we can replace Equation (10) with Since we have 5 variables and 3 units, we can rewrite Equation (12) based on two dimensionless numbers The first dimensionless number can be defined as [elastic forces that contrast deformation (in the solid)] [pressure forces that tend to deform the stent (from the liquid)] Knowing the typical ranges of E and d for the stent, and ρ and v for the blood, we can calculate the typical range of Π 1 . The second parameter Π 2 can be defined as We can have different types of Π 2 according to the type of stress we use in Equation (15). The stress tensor has 6 independent components that can be composed in different ways to provide different types of information. One possibility is to use the Frobenius norm.
In this case, we have a Π 2 based on the Frobenius norm that expresses, in dimensionless form, the total stress in the stent. Another possibility is the von Mises stress which provides another Π 2 number defined as Physically, Π F 2 and Π V 2 are dimensionless stresses and can have both a local form Π 2 (x, y, z) (when we calculate them at each x, y, z position), and a global form <Π 2 > (when we average them over the whole stent). Table 1 shows the parameters used in the simulation.

Results and Discussion
Three flow velocities of 0.4 ms −1 , 0.23 ms −1 and 0.16 ms −1 were chosen to represent the normal coronary artery and the baseline flow, respectively. The value of k was chosen from 0.5 to 5 to cover materials with the lowest to highest Young modulus. The blood flow velocity observed within the stent ranges from 0.23 ms −1 to 0.4 ms −1 , whereas the minimum flow velocity which may occur due to stenosis is 0.16 ms −1 and taken to be the baseline [52]. The velocity profile at different viewpoints is shown in Figure 2. Physically, Π F 2 and Π V 2 are dimensionless stresses and can have both a local form Π2 (x, y, z) (when we calculate them at each x, y, z position), and a global form <Π2> (when we average them over the whole stent). Table 1 shows the parameters used in the simulation.

Results and Discussion
Three flow velocities of 0.4 ms , 0.23 ms and 0.16 ms were chosen to represent the normal coronary artery and the baseline flow, respectively. The value of was chosen from 0.5 to 5 to cover materials with the lowest to highest Young modulus. The blood flow velocity observed within the stent ranges from 0.23 ms to 0.4 ms , whereas the minimum flow velocity which may occur due to stenosis is 0.16 ms and taken to be the baseline [52]. The velocity profile at different viewpoints is shown in Figure 2.  The dimensionless von Mises stress is shown in Figure 3 for two stents at different Π 1 . Different values of Π 1 means different k and v. It is shown that the stress is more severe at the nodes (joints) with higher stress at higher Π 1 as clearly shown in Figure 3b. This may lead to potential stent failure (rapture) at the joins or size change of the stent resulting from compression or expansion. The expansion characteristics of a stent are the main causes of vascular wall injuries [53]. Either of these conditions (failure or size change) will cause severe pain and damage to the patient and lead to restenosis and or stent redeployment.
The dimensionless von Mises stress is shown in Figure 3 for two stents at different 〈 〉. Different values of 〈 〉 means different and . It is shown that the stress is more severe at the nodes (joints) with higher stress at higher 〈 〉 as clearly shown in Figure  3b. This may lead to potential stent failure (rapture) at the joins or size change of the stent resulting from compression or expansion. The expansion characteristics of a stent are the main causes of vascular wall injuries [53]. Either of these conditions (failure or size change) will cause severe pain and damage to the patient and lead to restenosis and or stent redeployment. The result is first presented in Figure 4 and shows how the stress varies with and . This is then sorted in dimensionless form and presented in Figures 5 and 6, which shows the average stress <Π F 2>, and <Π V 2> versus <Π1> in dimensionless form. If we use dimensionless numbers the three curves of Figure 4 collapse in only one curve (Figures 5 and 6).  The result is first presented in Figure 4 and shows how the stress varies with k and v. This is then sorted in dimensionless form and presented in Figures 5 and 6, which shows the average stress <Π F 2 >, and <Π V 2 > versus <Π 1 > in dimensionless form. If we use dimensionless numbers the three curves of Figure 4 collapse in only one curve (Figures 5 and 6).
The dimensionless von Mises stress is shown in Figure 3 for two stents at different 〈 〉. Different values of 〈 〉 means different and . It is shown that the stress is more severe at the nodes (joints) with higher stress at higher 〈 〉 as clearly shown in Figure  3b. This may lead to potential stent failure (rapture) at the joins or size change of the stent resulting from compression or expansion. The expansion characteristics of a stent are the main causes of vascular wall injuries [53]. Either of these conditions (failure or size change) will cause severe pain and damage to the patient and lead to restenosis and or stent redeployment. The result is first presented in Figure 4 and shows how the stress varies with and . This is then sorted in dimensionless form and presented in Figures 5 and 6, which shows the average stress <Π F 2>, and <Π V 2> versus <Π1> in dimensionless form. If we use dimensionless numbers the three curves of Figure 4 collapse in only one curve ( Figures 5 and 6).   The plot confirms that the stress can be effectively sorted out with two dimensionless numbers based on the Buckingham π theorem. The three curves of Figure 4 can be fit by the same function as indicated in Equation (13). This function can be approximated by the following correlation (dotted line in) The same approach can be used for the von Mises stress which also has a correlation Numerically, we identified that the dimensionless numbers computed can be used as the fundamental group of the system in which the stress can be express in terms of Π1 and Π2.
Due to the pulsatile flow, the stent contracts and expands during the simulation. This causes the diameter of the stent to change with the flow. Figure 7 shows that the percentage change in the stent's diameter is fluctuating. This is because the arterial blood flow, which contributed to the stent deformation, is pulsatile in nature [54,55], therefore, it is expected to have a nonlinear change in the diameter. Note that the deformation is not only a function of the pressure forces from the liquid (blood) but also the elastic forces from the solid (stent). For that reason, several oscillation modes occur at the same time and the diameter change is not a simple repetition of the pulsatile flow. The plot confirms that the stress can be effectively sorted out with two dimensionless numbers based on the Buckingham π theorem. The three curves of Figure 4 can be fit by the same function as indicated in Equation (13). This function can be approximated by the following correlation (dotted line in) The same approach can be used for the von Mises stress which also has a correlation Numerically, we identified that the dimensionless numbers computed can be used as the fundamental group of the system in which the stress can be express in terms of Π 1 and Π 2 .
Due to the pulsatile flow, the stent contracts and expands during the simulation. This causes the diameter of the stent to change with the flow. Figure 7 shows that the percentage change in the stent's diameter is fluctuating. This is because the arterial blood flow, which contributed to the stent deformation, is pulsatile in nature [54,55], therefore, it is expected to have a nonlinear change in the diameter. Note that the deformation is not only a function of the pressure forces from the liquid (blood) but also the elastic forces from the solid (stent). For that reason, several oscillation modes occur at the same time and the diameter change is not a simple repetition of the pulsatile flow. Another likely incidence of vascular wall injury associated with stent expansion which can be quantified using the model is the so-called dogboning (DB) effect/ratio. This phenomenon occurs when the stent expands at the ends, resulting in increased stress and injuries at the arterial wall [53]. This occurs when the diameter expands at both ends of Another likely incidence of vascular wall injury associated with stent expansion which can be quantified using the model is the so-called dogboning (DB) effect/ratio. This phenomenon occurs when the stent expands at the ends, resulting in increased stress and injuries at the arterial wall [53]. This occurs when the diameter expands at both ends of the stent and contracts at the centre. The dogboning ratio is defined as, where D max,end is the maximum stent diameter at the end (distal and proximal), and D min,central is the minimum stent diameter at the centre. Our model is capable of analysing the dogboning ratio which could be used to assess and reduce the potential risk of vascular wall injury, and therefore, in this study, DB was calculated to be 4.4%. This is less than the 6.3% reported by [25] for PS-shaped stent.

Conclusions
In this paper, a Discrete Multiphysics model is used to simulate a coronary stent accounting for both its hemodynamics and mechanical stress. The model is three-dimensional and includes both the fluid (blood) and the solid structures (arterial wall and stent) and it is used to study the link between the flow dynamics and the mechanical deformation of the stent. The mechanical stress is computed using dimensionless numbers and a relationship between elastic forces and pressure forces was established. The results show that the blood flow contributes significantly to the stent deformation and the stiffness of the stent material affected the rate of deformation. Nonuniform stress distributions are observed. In particular, high stresses are observed at the nodes of the stent.
Given a specific Π 1 and from the corresponding Π V 2 , the maximum stress can be obtained. However, a fracture is not directly accounted for. This will depend on the pressure (stress) and type of material used for the stent as reported by [19]. Different materials have different yield stress (above which the stent fracture occurs). Our model can be used for all types of material by comparing the maximal stress in the structure against the material yield stress. Using the correlation in Equation (21), maximum Π V 2 . can be found when a specific stent material's property (Young modulus) is known. That means, from the value of Π 1 , we will be able to predict the Von Mises stress from Π V 2 (which can be compared against the yield stress) as well as the flow velocity. For clinical purposes, one just needs to know the material and the blood flow velocity in the diseased artery to be able to predict whether the stent is at risk of fracture or not. Therefore, the model can be used to predict the deformation of the stent once in place and the conditions that can potentially cause its failure.