Abstract

Translate

Atherosclerosis is an accumulation of plaque, which can result in changes in blood flow in the vicinity, leading to severe heart attack. This paper presents a phenomenological fluid–structure interaction study of plaque rupture in stenosed bifurcated elastic arteries. We use the coupled monolithic Arbitrary Lagrange Euler (ALE) formulation for fluids and solids. We consider the Navier–Stokes equation to govern the non-Newton blood flow and linear elastic model for walls. We treat the interface as a continuum. We utilize the stable $P_{2} P_{1}$ finite element pair for velocity and pressure discretization in space. The nonlinear discretized algebraic system is tackled using the Newton method, with the Jacobian matrices approximated via a divided differences approach. The resulting linear systems are addressed using the direct solver MUltifrontal Massively Parallel Sparse direct Solver (MUMPS). We then determine the wall shear stress (WSS) for both minimum and maximum times, accounting for elastic walls. The study’s findings enhance our understanding of the mechanisms behind plaque rupture and aid in developing better diagnostic and therapeutic strategies.

Full text

Turn on search term navigation

Translate

1. Introduction

Atherosclerosis is a disease that is influenced by various risk factors, including hypercholesterolemia, modified lipoproteins, high blood pressure, diabetes, infections, and smoking. Researchers have been studying atherosclerosis for over 100 years, with a focus on the role of cholesterol in its development, as first reported by Anitschkow and Chalatov. Atherosclerotic plaques can be classified as either ‘vulnerable’ or ‘stable’, with vulnerable plaques having a significant lipid core and a thin bottom cap that separates tissue factor thrombogenic macrophages from the blood. Owasit et al. [1] found that the stenosis geometry of certain cardiovascular patients cannot be described by a vertically symmetric function across the stenosis. In their study, Zuhaila et al. [2] analyzed the dynamic response of blood flow heat transmission in stenotic conditions through a bifurcated artery with a stiff wall (2D) and Newtonian, incompressible, laminar, and constant blood flow.

Li et al. [3] examined the impact of different computational models on hydrodynamic factors when an incompressible fluid interacts with two symmetric elastic or poroelastic structures. They conducted numerical experiments on blood flow using the Carreau–Yasuda model to simulate viscosity and investigate the influence of non-Newtonian blood rheology and poroelasticity on a benchmark vessel while also presenting a two-dimensional simulation of blood flow in an axisymmetric stenosis artery that considered both non-Newtonian fluid properties and fluid–structure interactions. The results showed that blood flow exhibits non-Newtonian behavior in small vessels and in settings with complicated geometry.

Golam et al. [4] numerically investigated the effects of non-Newtonian modeling on unsteady periodic flows in a 2D pipe with two idealized stenoses of 75 percent and 50 percent degrees, respectively, to handle complex geometries. The study examined various non-Newtonian blood constitutive equations, including the (i) Carreau, (ii) Cross, (iii) Modified Casson, and (iv) Quemada models, and compared the results with the Newtonian viscosity model. Wang [5] proposed a mathematical model to convert heat to blood flow in a narrow tube, while Liu et al. [6] examined the modeling of blood flow in pulmonary arterial fluid–structure interactions with a single, multidisciplinary, and continuous variant formula. A general framework for patient-specific simulations for Blood Flow FSI was developed, which included medical image mesh production, the low-order variational multiscale model (VMS) finite element, compressible and incompressible material formulation, and boundary conditions such as downstream coupling of closed-loop lumped-parameter network (LPN) circulation and time integration system. The five crucial parameters influencing blood concentrations are plasma viscosity, hematocrit, red blood cell adhesion, the accumulation of red blood cells, and temperature.

Ali et al. and Raza et al. (2024) investigated magnetohydrodynamic (MHD) hybrid nanofluid flow in porous microchannels using fractal-fractional derivatives, analyzing velocity, heat transfer, and pressure gradients while comparing simple and hybrid suspensions with numerical validation [7,8]. Also, fluid–structure interaction analysis of biomagnetic blood flow and plaque growth in stenotic bifurcated arteries was reported in [9,10].

Carlo Palombo et al. [11] examined the mechanism for arterial stiffness development and progression, summarized the evidence, caution, and clinical applications of stiffness as a substitute marker for the risk of cardiovascular disease, artery strength, atherosclerotic burden, and cardiovascular disruption. The phenomena associated with the transit of hemodynamic bioheat have been found to have a growing impact on the development of atherogenic processes. It is essential to understand the distortion of temperature distribution as a function of vessel diameter in order to create adequate bioheat transport models. In physiological settings, when the diameter of the blood vessel is large, the thermal distribution is disrupted [12]. It is worth noting that localized areas of cooling can be detected in heated hyperthermic tissues when blood flow is present according to a prior study [13,14]. The typical human blood temperature is approximately 37 °C, and irreversible harmful effects can occur in blood proteins at such high temperatures, leading to cell death [15].

This paper is structured as follows: Section 2 details the formulation of the equations using the Arbitrary Lagrangian–Eulerian (ALE) method. Section 3 focuses on the discretization of the flow problem in two spatial dimensions, employing the $P_{2} P_{1}$ finite element pair within a standard finite element method (FEM) framework. Section 4 addresses the resulting discrete problems in a fully coupled, monolithic manner for displacement, velocity, and pressure $(u, v, p)$ utilizing outer Newton iterations and an inner direct solver. In Section 5, the problem configuration is described. Finally, numerical results and discussions for plaque rupture bifurcated stenosis are presented.

2. Bifurcated Stenosis Fluid-Structure Interaction Problem

We formulate the bifurcated fluid–structure interaction problem of the fluid and solid domains and the conditions to the interfaces. We define $Ω_{t}^{f}$ as a fluid, i.e., blood, and $Ω_{t}^{s}$ for elastic blood walls at time $\in [0, T]$ . Let $Γ_{t}^{0} = {\bar{Ω}}_{t}^{f} \cap {\bar{Ω}}_{t}^{s}$ be the boundary part where blood flow interacts with the elastic walls, and let $Γ_{t}^{i}, i = 1, 2, 3$ be the other boundaries of the structure and the fluid. The mapping $Ψ^{s}$ is defined as

$\begin{matrix} Ψ^{s} & : Ω^{s} \times [0, T] \mapsto R^{3}, \end{matrix}$

and the displacement

u^{s}

and the velocity

v^{s}

are

$u^{s} (X, t) = Ψ^{s} (X, t) - X, v^{s} (X, t) = \frac{\partial Ψ^{s}}{\partial t} (X, t) .$

The blood flow velocity field $v^{f}$ on the Eulerian current configuration $Ω_{t}^{f}$ is

$v^{f} (x, t) : Ω_{t}^{f} \times [0, T] \mapsto R^{3},$

with auxiliary sufficient smooth mappings satisfying

(1) $\begin{matrix} ζ^{f} : Ω^{f} \times [0, T] \mapsto R^{3}, u^{f} (X, t) = ζ^{f} (X, t) - X, \end{matrix}$

(2) $\begin{matrix} ζ^{f} (X, t) = Ψ^{s} (X, t), \forall (X, t) \in Γ^{0} \times [0, T] . \end{matrix}$

Equation (2) is the artificial mesh deformation for the fluid. The balance equations of the fluid (blood) are

(3) $\begin{matrix} ϱ^{f} \frac{\partial v^{f}}{\partial t} + ϱ^{f} (\nabla v^{f}) (v^{f} - \frac{\partial u^{f}}{\partial t}) = div σ^{f} & Ω_{t}^{f}, \end{matrix}$

(4) $\begin{matrix} div v^{f} = 0 & Ω_{t}^{f}, \end{matrix}$

where

\frac{\partial u^{f}}{\partial t}

is the mesh moving velocity. Similarly, the balance equations of the solid (elastic) domain are

(5) $\begin{matrix} ϱ^{s} \frac{\partial v^{s}}{\partial t} + ϱ^{s} (\nabla v^{s}) v^{s} = div σ^{s} & Ω_{t}^{s}, \end{matrix}$

(6) $\begin{matrix} div v^{s} = 0 & Ω_{t}^{s} . \end{matrix}$

The interaction arises from the transfer of momentum across the shared part of the boundary $Γ_{t}^{0}$ , where the forces are in balance, and the no-slip boundary condition for the fluid is

(7) $\begin{matrix} σ^{f} n = σ^{s} n Γ_{t}^{0}, & v^{f} = v^{s} Γ_{t}^{0} . \end{matrix}$

A natural boundary condition on part $Γ_{t}^{1}$ is

$\begin{matrix} σ^{f} n = p_{B} n & Γ_{t}^{1}, \end{matrix}$

with

p_{B}

being an assigned value. A Dirichlet boundary condition on

Γ_{t}^{1}

can also be prescribed as

$\begin{matrix} v^{f} = v_{B} & Γ_{t}^{1}, \end{matrix}$

where

v_{B}

is given and the solid displacement at

Γ_{t}^{2}

defined as

$\begin{matrix} u^{s} = 0 & Γ_{t}^{2}, \end{matrix}$

and the stress-free boundary condition at

Γ_{t}^{3}

is defined as

$\begin{matrix} σ^{s} n = 0 & Γ_{t}^{3} . \end{matrix}$

Considering $Ω = Ω^{f} \cup Ω^{s}$ , where $Ω^{f}, Ω^{s}$ are fluid and solid domains at the initial state, i.e.,

$\begin{matrix} u & : Ω \times [0, T] \to R^{3}, & v & : Ω \times [0, T] \to R^{3}, \end{matrix}$

$v$ is the velocity, and $u$ is the displacement on the respective domains:

$v = \{\begin{matrix} v^{s} Ω^{s}, \\ v^{f} Ω^{f}, \end{matrix} u = \{\begin{matrix} u^{s} Ω^{s}, \\ u^{f} Ω^{f} . \end{matrix}$

By the Equations (2) and (7), $v$ and $u$ are continuous across the interface $Γ_{t}^{0}$ and global quantities on $Ω$ as the deformation gradient and its determinant, which are defined as follows:

(8) $\begin{matrix} F = I + Grad u, & J = det F . \end{matrix}$

Using the Equation (8), the Equations (5) and (6) become in the Lagrangian frame:

(9) $\begin{matrix} J ϱ^{s} \frac{d v}{d t} = Div P^{s} & Ω^{s}, \end{matrix}$

(10) $\begin{matrix} J = 1 & Ω^{s} . \end{matrix}$

The Equations (3) and (4) are in the arbitrary Lagrangian–Eulerian formulation in $Ω_{t}^{f}$ , now bringing the equations to $Ω^{f}$ by the mapping $ζ^{f}$ defined by (1):

$\begin{matrix} ϱ^{f} \frac{\partial v}{\partial t} + ϱ^{f} (Grad v) F^{- 1} (v - \frac{\partial u}{\partial t}) = J^{- 1} Div (J σ^{f} F^{- T}) & Ω^{f}, \\ Div (J v F^{- T}) = 0 & Ω^{f} . \end{matrix}$

The relation for the mapping $ζ^{f}$ for displacement $u^{f}$ , together with the Dirichlet boundary conditions required by (2), is defined as follows:

$\begin{matrix} \frac{\partial u}{\partial t} = Δ u Ω^{f}, & u = u^{s} Γ^{0}, & u = 0 Γ^{1} . \end{matrix}$

The full set of the equations is

(11) $\begin{matrix} \frac{\partial u}{\partial t} = & \{\begin{matrix} v & Ω^{s}, \\ Δ u & Ω^{f}, \end{matrix} \end{matrix}$

(12) $\begin{matrix} \frac{\partial v}{\partial t} = & \{\begin{matrix} \frac{1}{J ϱ^{s}} Div P^{s} & Ω^{s}, \\ - (Grad v) F^{- 1} (v - \frac{\partial u}{\partial t}) + \frac{1}{J ϱ^{f}} Div (J σ^{f} F^{- T}) & Ω^{f}, \end{matrix} \end{matrix}$

(13) $\begin{matrix} 0 = & \{\begin{matrix} J - 1 & in Ω^{s}, \\ Div (J v F^{- T}) & Ω^{f}, \end{matrix} \end{matrix}$

with the initial conditions (ICs) defined as

$\begin{matrix} u (0) = 0 Ω, & v (0) = v_{0} Ω, \end{matrix}$

and the boundary conditions (BCs) defined as

$u = 0, v = v_{B} Γ^{1}, u = 0 Γ^{2}, σ^{s} n = 0 Γ^{3} .$

The ALE method is highly effective in fluid–structure interaction (FSI) problems, as it allows the computational mesh to adapt dynamically to deformations in both the fluid and the structure. By continuously adjusting the mesh while maintaining the physical integrity of the interface, the ALE method ensures stable and accurate simulations. This approach strikes a balance between computational efficiency and precision, making it possible to model complex FSI scenarios that would be difficult to handle using purely Lagrangian or Eulerian methods [16,17,18].

Constitutive Equations

We consider the laminar flow of an incompressible viscous fluid in a bifurcated artery with a stenosis. The stress tensor for the fluid is defined as follows:

$\begin{matrix} σ^{f} = - p^{f} I + 2 μ^{f} D, \end{matrix}$

where

D = \frac{1}{2} (\nabla v^{f} + ({\nabla v^{f})}^{T}

is the symmetric part of the gradient velocity, and

p^{f}

is the Lagrange multiplier. To incorporate the non-Newtonian effect in our analysis, we use the Casson–Papanastasiou model, which includes an exponential term suggested by Papanastasiou [19]. This term eliminates the need for a stress threshold in this situation. The regularized Casson model is expressed as follows:

$μ^{f} = {[\sqrt{α_{T} μ_{p}} + \sqrt{\frac{τ_{y}}{\dot{γ}}} (1 - e^{- \sqrt{m_{p} \dot{γ}}})]}^{2},$

where

\dot{γ} = \sqrt{2 D : D}

. In addition, we consider the blood vessel walls to be elastic, using a mechanically homogeneous neo-Hookean isotropic hyperelastic material model. This model is suitable for both compressible materials (where

ν^{s} < 1 / 2

) and incompressible materials (where

ν^{s} = 1 / 2

). The Cauchy stress tensor

σ^{s}

is defined by one of the following constitutive laws: The neo-Hooke compressible material (

J = det F

) is defined as

$\begin{matrix} σ^{s} = λ^{s} (J - \frac{1}{J}) I + \frac{μ^{s}}{J} (F F^{T} - I) . \end{matrix}$

Similarly, an alternate compressible neo-Hooke model is possible (for instance, $σ^{s} = λ^{s} log (J) I + \frac{μ^{s}}{J} (F F^{T} - I)$ or $σ^{s} = λ^{s} (J - 1) I + \frac{μ^{s}}{J} (F F^{T} - I)$ ), which is reported in the literature as leading to a similar behavior for small volumetric changes.

The neo-Hooke incompressible model is

$\begin{matrix} σ^{s} = - p^{s} I + μ^{s} (F F^{T} - I), det F = 1, \end{matrix}$

and the first Piola-Kirchhoff stress tensor is

$\begin{matrix} P^{s} = - J p^{s} F^{- T}, \end{matrix}$

where

p^{s}

is the Lagrange multiplier for the incompressibility constraint (10). The material’s elasticity is characterized by two parameters: the Poisson ratio

ν^{s}

and the Young’s modulus E. These parameters allow us to determine the Lamé constants

λ^{s}

and the shear modulus

μ^{s}

using the following relations:

$ν^{s} = \frac{λ^{s}}{2 (λ^{s} + μ^{s})} E = \frac{μ^{s} (3 λ^{s} + 2 μ^{s})}{(λ^{s} + μ^{s})}$

$μ^{s} = \frac{E}{2 (1 + ν^{s})} λ^{s} = \frac{ν^{s} E}{(1 + ν^{s}) (1 - 2 ν^{s})} .$

The stenosis walls can show anisotropic behavior in the case of more complex blood flow. More realistic constitutive relations can be employed in this case (for details, see, for example, [20,21]).

3. Finite Element Method (FEM) Discretization

The Taylor–Hood stable finite element pair $P_{2} P_{1}$ is utilized for its optimal convergence properties; it is the most accurate and robust finite element pair for incompressible viscous flow and satisfies the Ladyzhenskaya-Babuška–Brezzi (LBB) stability condition [22]. The velocity and displacement are discretized by a quadratic element $P_{2}$ and the pressure by $P_{1}$ . Figure 1 illustrates the corresponding degrees of freedom.

The basis function and basis vector for $P_{2}$ and $P_{1}$ , respectively, (for details, see [23]) are defined as follows:

$ϕ_{i} (x) = a_{i} + b_{i} x + c_{i} y + d_{i} x^{2} + e_{i} y^{2} + f_{i} x y, 〈 1, x, y, x^{2}, x y, y^{2} 〉$

$ϕ_{i} (x) = a_{i} + b_{i} x + c_{i} y, 〈 1, x, y 〉 .$

Variational Formulation

The standard Galerkin finite element method is applied to the Equations (11)–(13) by multiplying the test functions $ζ, ξ, γ$ and integrating over the space and time interval I. Through integration by parts and using the boundary conditions, we reached the following:

$\begin{matrix} \int_{0}^{T} \int_{Ω} \frac{\partial u}{\partial t} \cdot ζ d V d t = & \int_{0}^{T} \int_{Ω^{s}} v \cdot ζ d V d t - \int_{0}^{T} \int_{Ω^{f}} Grad u \cdot Grad ζ d V d t, \end{matrix}$

$\begin{matrix} \int_{0}^{T} \int_{Ω^{f}} J \frac{\partial v}{\partial t} \cdot ξ d V d t & + \int_{0}^{T} \int_{Ω^{s}} β J \frac{\partial v}{\partial t} \cdot ξ d V d t \\ = & - \int_{0}^{T} \int_{Ω^{f}} J Grad v F^{- 1} (v - \frac{\partial u}{\partial t}) \cdot ξ d V d t \\ + \int_{0}^{T} \int_{Ω} J p F^{- T} \cdot Grad ξ d V d t \\ - \int_{0}^{T} \int_{Ω^{f}} J μ Grad v F^{- 1} F^{- T} \cdot Grad ξ d V d t, \end{matrix}$

$0 = \int_{0}^{T} \int_{Ω^{s}} (J - 1) γ d V d t + \int_{0}^{T} \int_{Ω^{f}} Div (J v F^{- T}) γ d V d t .$

where

β = \frac{ρ^{f}}{ρ^{s}}

. The resulting nonlinear algebraic system of equations reads as follows:

(14) $(\begin{matrix} A_{u u} & A_{u v} & 0 \\ A_{v u} & A_{v v} & k B \\ c_{u} B_{s}^{T} & c_{v} B_{f}^{T} & 0 \end{matrix}) (\begin{matrix} u_{h} \\ v_{h} \\ p_{h} \end{matrix}) = (\begin{matrix} {rhs}_{u} \\ {rhs}_{v} \\ {rhs}_{p} \end{matrix}),$

where

A

denotes the operator for diffusion, reaction, and convection from the governing equations, whereas

B

and

B^{T}

represent the discrete gradient and divergence operators, respectively.

4. Solution Algorithm

The Newton iteration method exhibiting the quadratic convergence is used for the nonlinear algebraic equations of the saddle point type (14). We find the root of the residual as

$R (X) = 0,$

The formula for the Newton iteration with damping is given as follows:

$X^{n + 1} = X^{n} + ω^{n} {[\frac{\partial R (X^{n})}{\partial X}]}^{- 1} R (X^{n}),$

where

X = (u_{h}, v_{h}, p_{h})

is the solution vector, n is the iteration number, and

ω^{n}

is the damping parameter. The Jacobian matrix

\frac{\partial R (X^{n})}{\partial X}

is computed by the finite differences as

(15) ${[\frac{\partial R (X^{n})}{\partial X}]}_{i j} \approx \frac{{[R]}_{i} (X^{n} + α_{j} e_{j}) - {[R]}_{i} (X^{n} - α_{j} e_{j})}{2 α_{j}},$

and in Equation (15), the coefficients

α_{j} > 0

represent the increments at each iteration step n taken as

1 E - 6

, while

e_{j}

denotes the unit basis vectors in

R^{n}

; for details, see [24]. For a more detailed explanation, please refer to [25,26]. Furthermore, the parameter

ω^{n} \in (- 1, 0)

is chosen such that the error measure decreases with each iteration as

$R (X^{n + 1}) \cdot X^{n + 1} \leq R (X^{n}) \cdot X^{n},$

and the damping parameter in the Newton iteration greatly enhances its robustness, particularly when the current approximation

X^{n}

is not close enough to the final solution. For more information, please refer to [25,27]. In this

2 D

problem, we employ a direct solver for sparse systems such as the MUltifrontal Massively Parallel sparse direct Solver (MUMPS) [28]. Algorithm 1 presents the Newton iteration algorithm used in this study. Algorithm 1 outlines the Newton iteration and line search method used in this study.

Algorithm 1 Newton iteration with line search

1: Tolerance parameter for nonlinear input 2: Start with

n = 0

and take

X^{n}

as the starting guess. 3: Residual vector construction

R (X^{n}) = A X^{n} - b

. 4: The Jacobian matrix computation

J (X^{n}) = \frac{\partial R}{\partial X} (X^{n})

. 5: Go for linear system correction of

δ X

$J (X^{n}) δ X^{n} = R (X^{n})$

6: Obtain the optimal step length

ω^{n} \in (- 1, 0]

. 7: Solution update:

X^{n + 1} = X^{n} + ω^{n} δ X^{n}

5. Problem Configuration

A parabolic flow profile is prescribed at the inlet:

$v_{i n}^{f} (x, y) = 8 y (1 - y)$

The outlet boundary conditions are prescribed on pressure. The pressure sets equal zero at the outlet. The boundary condition for the solid domain is a fixed point constraint at the start and end, meaning that those points do not move due to flow. Other parts of the domain are free to move. The outlet boundary conditions are prescribed on pressure, and the pressure is set to zero at the outlet.

Initially, the velocity and displacement are sets equal to zero. The parameters for the Casson–Papanastasiou models are $α_{T} = 1$ , $μ_{p} = 1 / R e = 0.005$ , $τ_{y} = 0.001$ , and $m_{p} = 100$ . Unless otherwise specified, we used Young’s modulus $E = 1 \times 10^{10}$ and Poisson’s ratio $ν = 0.49$ .

Computation domain: As shown in Figure 2, a prototype geometric model with stenosis and a bifurcation with a symmetrical arrangement is considered for the computational domain. The full length of the region is $L = 6.34$ . The diameter h is considered to be 1 and shrinks to $50 %$ at the region of stenosis, which implies that stenosis restricts the artery about $50 %$ of the time. The elastic wall’s width w is estimated to be $0.1$ . The diameter of the daughter artery is $0.37$ . The bifurcation artery has an inclination of $37^{°}$ . C is the central line along which pressure is recorded. A and B are selected to predict the velocity profile’s behavior before and after stenosis. The distance between lines a and b is calculated to be $d = 0.8$ . Table 1 shows the absolute error of wall shear stresses (WSS) as a function of the mesh refinement level and the number of elements. Additionally, it presents the WSS on the upper elastic wall. For Level 0 computations, we used a total of 884 finite elements.

6. Results and Discussions

In this section, the implications and significance of the findings are explained and interpreted in the context of the research hypotheses. The limitations of the study are discussed, and potential directions for future research are suggested. Moreover, the results are compared and contrasted with previous studies or theories where applicable.

6.1. Velocity Profile

Figure 3 displays the blood flow behavior at $R e = 200$ nonturbulent flow. As observed, the blood flow velocity decreased by $26 %$ , from about 2.3 units to 1.7 units, as it approached the stenotic region labeled ‘B’. The gradual increase in the arterial diameter by approximately $31.5 %$ at the beginning part of the stenotic region is the cause of this drop in velocity. Moreover, the lower wall of the artery expanded more than the upper wall due to its Young’s modulus being around $25 %$ less than the upper wall, which is more rigid because of the plaque’s inherent property. The decrease in blood velocity in the expanded lower wall caused a region of relatively high pressure normal to the wall surface, resulting in arterial wall expansion.

It is worth noting that backflow nucleation occurred at the bifurcation point as the blood divided into branches with smaller diameters, causing an increase in velocity compared to the main artery by around $30 %$ . These findings provide insight into the effects of stenosis on blood flow behavior and arterial wall expansion. The results also highlight the significance of considering the heterogeneous properties of arterial walls when analyzing blood flow behavior in stenotic arteries.

6.2. Blood Backflow

The velocity vector simulation results of blood flow at Reynolds number $R e = 200$ are presented in Figure 4. The left image shows the normalized velocity vectors of the blood flow as it moves across the artery, revealing significant backflow in the trailing edge of the plaque, even at a low Reynolds number. However, it is not evident from this image whether there is any stagnation or other irregularity in the blood flow patterns. A low-velocity values, high-pressure effects are showcased in Figure 5.

The right image in Figure 4 shows the actual velocity of the blood flow, which is indicated by the length of the arrows. It is apparent from this image that the relative velocity of blood in the trailing edge of the plaque is much lower, by over $90 %$ . This low flow rate can result in blood stagnation and the formation of micro-vortices, which can lead to activation and stress on blood platelets, further exacerbating the process of atherosclerosis. Such changes in the flow pattern of the blood flow can increase the likelihood of blood particle deposition at the stenosis region and subsequent thrombosis.

6.3. Higher Reynolds Number

Figure 6 shows the lower wall deformation for a Reynolds number of $R e = 2000$ , which is relatively small due to the high velocity and low normal pressure. As the Reynolds number increased from $R e = 200$ to 2000, the effects of backflow, stagnation, vortices, and low lateral velocities were further amplified, as expected. Figure 6 and Figure 7 illustrate these effects.

Moreover, as the velocity of the blood flow increased, the normal pressure on the arterial wall decreased. This resulted in a small deformation of the arterial wall, which is apparent when comparing Figure 3 with Figure 6 and Figure 7.

Figure 7 also shows that the stagnation region grew significantly with a Reynolds number of $R e = 2000$ .

6.4. Dynamic Energized Hotspots

At the backflow region labeled ‘H’ in Figure 8, the blood flow characteristics can create an environment conducive to the formation of highly Dynamic Energized Hot Spots (DEHSs).

DEHSs are micro- to mini-sized vortices with low lateral speed (stagnation) but high rotational inertia. The vortices increase both the normal and shear stress in this region, causing further lesions to the plaque. This, in turn, activates the platelets and recruits more blood cells to the plaque, causing it to grow larger. Due to the low lateral speed, the damage to the artery wall and plaque can be severe. If these DEHSs are frequently in contact with a vulnerable region of the plaque, it can cause it to rupture, leading to thrombosis.

The severity of the damage caused by DEHSs emphasizes the need to understand the dynamics of blood flow and the associated fluid mechanics in atherosclerotic arteries. Figure 8 provides a higher magnified view of the trailing end of the artery where there is significant backflow, which can serve as a breeding ground for the formation of DEHSs.

6.5. Shear Stress

The slight narrowing of the arterial wall in the stenotic region caused the blood velocity to increase as the blood flow moved across the region. Figure 3 shows that the average velocity changes from high (2.2) to low (1.6) and then high (2.5) as the blood flowed across the artery. It is worth noting that the high-velocity resulted in greater shear stress at the plaque region. Specifically, as the blood flowed at its highest velocity through the narrowest region or the cap of the plaque, the shear stress was the greatest, as depicted in Figure 9.

The high shear stress can cause the fracture of the plaque cap, leading to thrombus formation and even to myocardial infarction. Understanding the effects of shear stress on plaque vulnerability is essential to prevent and manage atherosclerosis-related complications. Figure 9 highlights the importance of assessing and monitoring shear stress levels in patients with stenotic plaques.

6.6. Wall Shear Stress

Observing wall shear stress (WSS) is essential to reduce the recirculation area in stenosed arteries, which in turn can lower the risk of atherosclerosis. Figure 10 illustrates the behavior of WSS in the stenosis for different Reynolds numbers ( $R e$ ). The figure shows that as $R e$ decreased, the magnitude of WSS increased. This indicates that a decrease in $R e$ can reduce the chance of atherosclerosis in the stenosed region by increasing the total WSS experienced by the stenosis. Furthermore, it has been observed that WSS can play a crucial role in the formation and development of atherosclerotic plaques. A low or disturbed WSS can cause the endothelium to become dysfunctional and lead to the initiation of atherosclerosis. On the other hand, a high and unidirectional WSS can maintain the endothelium in a healthy state and protect against atherosclerosis. Therefore, the measurement and analysis of the WSS can provide valuable insights into the progression and treatment of cardiovascular diseases.

6.7. Wall Displacement

The deformation of the arterial wall is an important factor to consider, as it can affect the overall health of the artery. In Figure 11, the behavior of wall displacement is shown. It can be observed that there is a comparatively small deformation at the stenosis position. This is due to the difference in elastic modulus between the stenotic and nonstenotic walls. As the value of Reynolds number $R e$ increased, the displacement or deformation of the wall decreased.

It is important to note that excessive deformation of the arterial wall can lead to the rupture of the plaque and cause further damage to the artery. Therefore, monitoring the wall displacement can provide insight into the overall health of the artery and can be used as an early warning sign for potential health problems. Further research is needed to determine the optimal range of wall displacement for healthy arterial function.

7. Summary and Conclusions

In summary, the behavior of blood flow in a stenosed artery can have significant impacts on the development and progression of atherosclerosis. High shear stress and wall displacement are key factors to observe when analyzing the behavior of blood flow. Shear stress, caused by the high velocity of blood flow in the narrowest region or the cap of the plaque, can result in the fracture of the cap of the plaque, leading to thrombosis. Wall displacement or deformation, caused by the different elastic modulus at the stenosis and the other wall, can reduce as the value of the Reynolds number increases. Wall shear stress is another important factor to observe, as it helps to reduce the recirculation area in a stenosed artery. The magnitude of the WSS increases with a decrease in Reynolds number, which in turn, reduces the chance of atherosclerosis.

The computer model and simulation results have demonstrated a notable alteration in the blood flow dynamics as blood moves through the stenotic region of the artery. The observed changes encompass the expansion of the arterial wall, modifications in the velocity profiles, and the presence of blood stagnation and backflow. Based on the obtained results, it has been suggested that highly dynamic energized hot spots (DEHSs) may form, and this can lead to the development of a larger plaque, which can further narrow the artery lumen. Moreover, the simulation results have also revealed that the high shear stress at the narrowest part of the artery, which is located at the cap of the plaque, may cause it to rupture.

In conclusion, a better understanding of the behavior of blood flow in a stenosed artery can aid in the prevention and treatment of atherosclerosis. Factors such as shear stress and wall displacement should be considered when analyzing the risk of atherosclerosis in a patient. Future research could focus on developing more accurate and comprehensive models of blood flow behavior in a stenosed artery, as well as investigating new methods for preventing and treating atherosclerosis.

Author Contributions

Conceptualization, M.R., K.I. and M.A.A.; methodology, M.R., K.I. and M.A.A.; validation, M.R., M.A.A., I.H. and M.G.; formal analysis, M.R., K.I., I.H. and M.A.A.; writing—original draft preparation, M.R., K.I. and M.A.A.; writing—review and editing, M.R. and M.A.A.; supervision, M.R. and M.G. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data will be made available by the authors on request.

Acknowledgments

The authors sincerely appreciate the anonymous referees for their valuable comments, which greatly improved the presentation of our results.

Conflicts of Interest

The authors state that they did not have any conflicts of interest related to this study.

Abbreviations and Nomenclature

The following abbreviations are used in this manuscript:

ALE	Arbitrary Lagrange Euler
FSI	Fluid Structure Interaction
WSS	Wall Shear Stress
SVK	Saint Venant–Kirchhof
MUMPS	MUltifrontal Massively Parallel Sparse direct Solver
DEHSs	Dynamic Energized Hot Spots
MHD	Magnetohydrodynamic
u	Displacement vector (elastic or fluid domains)
v	Velocity vector (fluid or solid)
p	Pressure in the fluid domain
$σ$	Stress tensor (Cauchy stress)
$μ$	Dynamic viscosity
$ρ$	Density of the material (fluid or solid)
$τ$	Shear stress
$γ$	Shear rate
J	Determinant of the deformation gradient
F	Deformation gradient tensor
$λ$	First Lamé parameter
$ν$	Poisson’s ratio
D	Rate of deformation tensor
P	First Piola–Kirchhoff stress tensor
A	Coefficient matrix for diffusion, reaction, and convection
B	Gradient operator in discretized equations
$B^{T}$	Transpose of the gradient operator
R	Residual vector
$α$	Regularization or damping parameter
$ω$	Step size or weighting factor in iterations
E	Young’s modulus (elastic modulus)
$μ^{s}$	Shear modulus (second Lamé parameter)

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. The position of the degrees of freedom for the [Forumla omitted. See PDF.] element.

Figure 2. A schematic diagram of the problem and coarse mesh.

Figure 3. Velocity contour maps with Reynolds number [Forumla omitted. See PDF.].

Figure 4. Normalized and actual velocity vector results showing blood backflow in the trailing end of the plaque.

Figure 5. Low-velocity–high-pressure effect causing the arterial wall to expand.

Figure 6. [Forumla omitted. See PDF.] lower wall deformation is very small due to high velocity and low normal pressure.

Figure 7. With [Forumla omitted. See PDF.], the stagnation region grew significantly.

Figure 8. Higher magnified region of the trailing end of the artery where there is significant backflow.

$View Image - Figure 9. High shear stress can result in the fracture of the cap of the plaque.$

Figure 9. High shear stress can result in the fracture of the cap of the plaque.

Figure 10. Total wall shear stress experienced by the stenosis from A to B.

Figure 11. Total displacement of the upper wall.

Table 1

For $R e = 200$ : Absolute error for the WSS against refinement level showing mesh independence.

Level	Number of Elements	WSS	Abs. Error
0	884	0.487301	-
1	1149	0.487350	0.000091
2	1459	0.487389	0.000105
3	4523	0.487470	0.000109
4	13,297	0.487500	0.000030
5	14,360	0.487500	0.000003

References

1. Owasit, P.; Sriyab, S. Mathematical modeling of non-Newtonian fluid in arterial blood flow through various stenoses. Adv. Differ. Equ.; 2021; 2021, 340. [DOI: https://dx.doi.org/10.1186/s13662-021-03492-9]

2. Jamali, M.S.A.; Ismail, Z. Simulation of Heat Transfer on Blood Flow through a Stenosed Bifurcated Artery. J. Adv. Res. Fluid Mech. Therm. Sci.; 2019; 60, pp. 310-323.

3. Li, T.; Wang, X.; Yotov, I. Non-Newtonian and poroelastic effects in simulations of arterial flows. arXiv; 2020; arXiv: 2010.14072

4. Mir, G.R.; Shupti, S.P.; Molla, M.M. Pulsatile Non-Newtonian Laminar Blood Flows through Arterial Double Stenoses. J. Fluids; 2014; 2014, 757902. [DOI: https://dx.doi.org/10.1155/2014/757902]

5. Wang, C.Y. Heat Transfer to Blood Flow in a Small Tube. J. Biomech. Eng.; 2008; 130, pp. 024501-1-024501-3. [DOI: https://dx.doi.org/10.1115/1.2898722]

6. Liu, J.; Yang, W.; Lan, I.S.; Marsden, A.L. Fluid-structure interaction modeling of blood flow in the pulmonary arteries using the unified continuum and variational multiscale formulation. Mech. Res. Commun.; 2020; 107, pp. 103556-103563. [DOI: https://dx.doi.org/10.1016/j.mechrescom.2020.103556] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32773906]

7. Ali, A.H.; Raza, A.; Batiha, B.; Abed, A.M.; Abduljabbar, Z.A. Exact and fractional solution of MHD generalized Couette hybrid nanofluid flow with Mittag-Leffler and power law kernel. Int. J. Thermofluids; 2024; 24, 100837. [DOI: https://dx.doi.org/10.1016/j.ijft.2024.100837]

8. Raza, A.; Stadoleanu, O.V.; Abed, A.M.; Ali, A.H.; Sallah, M. Heat transfer model analysis of fractional Jeffery-type hybrid nanofluid dripping through a poured microchannel. Int. J. Thermofluids; 2024; 22, 100656. [DOI: https://dx.doi.org/10.1016/j.ijft.2024.100656]

9. Iqbal, K.; Rossi di Schio, E.; Anwar, M.A.; Razzaq, M.; Shahzad, H.; Valdiserri, P.; Fabbri, G.; Biserni, C. A Fluid–Structure Interaction Analysis to Investigate the Influence of Magnetic Fields on Plaque Growth in Stenotic Bifurcated Arteries. Dynamics; 2024; 4, pp. 572-591. [DOI: https://dx.doi.org/10.3390/dynamics4030030]

10. Anwar, M.A.; Iqbal, K.; Razzaq, M. Analysis of biomagnetic blood flow in a stenosed bifurcation artery amidst elastic walls. Phys. Scr.; 2021; 96, 085202. [DOI: https://dx.doi.org/10.1088/1402-4896/abf67b]

11. Carlo, P.; Kozakova, M. Arterial stiffness, atherosclerosis and cardiovascular risk: Pathophysiologic mechanisms and emerging clinical indications. Vasc. Pharmacol.; 2016; 77, pp. 1-7.

12. Zhang, B.; Gu, J.; Qian, M.; Niu, L.; Zhou, H.; Ghista, D. Correlation between quantitative analysis of wall shear stress and intima-media thickness in atherosclerosis development in carotid arteries. Biomed. Eng. Online; 2017; 16, 137. [DOI: https://dx.doi.org/10.1186/s12938-017-0425-9] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29208019]

13. Chinyoka, T.; Makinde, O.D. Computational dynamics of arterial blood flow in the presence of magnetic field and thermal radiation therapy. Adv. Math. Phys.; 2014; 2014, 915640. [DOI: https://dx.doi.org/10.1155/2014/915640]

14. Siegel, R.; Perlmutter, M. Heat transfer for pulsating laminar duct flow. J. Heat Transf.; 1962; 84, pp. 111-122. [DOI: https://dx.doi.org/10.1115/1.3684307]

15. Chato, D. Heat transfer to blood vessels. J. Biomech. Eng.; 1980; 102, pp. 110-118. [DOI: https://dx.doi.org/10.1115/1.3138205]

16. Turek, S.; Hron, J.; Razzaq, M.; Wobker, H.; Schäfer, M. Fluid-Structure Interaction II-Modelling, Simulation. Optimization; 2010; 73, pp. 413-424. [DOI: https://dx.doi.org/10.1007/978-3-642-14206-2_8]

17. Razzaq, M. Numerical Techniques for Solving Fluid-Structure Interaction Problems with Applications to Bioengineering. Ph.D. Thesis; TU Dortmund: Dortmund, Germany, 2009.

18. Razzaq, M.; Hron, J.; Turek, S. Numerical Simulation of Laminar Incompressible Fluid-Structure Interaction for Elastic Material with Point Constraints. Advances in Mathematical Fluid Mechanics; Rannacher, R.; Sequeira, A. Springer: Berlin/Heidelberg, Germany, 2010; pp. 451-472. [DOI: https://dx.doi.org/10.1007/978-3-642-04068-9]

19. Papanastasiou, T.C. Flows of materials with yield. J. Rheol.; 1987; 31, pp. 385-404. [DOI: https://dx.doi.org/10.1122/1.549926]

20. Ciarlet, P.G. Mathematical Elasticity: Three-Dimensional Elasticity; Volume 20 of Studies in Mathematics and its Applications Society for Industrial and Applied Mathematics: Amsterdam, The Netherlands, 1994; Volume I.

21. Tallec, P.L.; Mani, S. Numerical analysis of a linearised fluid-structure interaction problem. Num. Math.; 2000; 87, pp. 317-354. [DOI: https://dx.doi.org/10.1007/s002110000183]

22. Ervin, V.J.; Jenkins, E.W. The LBB Condition for the Taylor-Hood P2-P1 and Scott-Vogelius P2-discP1 Element Pairs in 2-D; Technical Report TR2011 04 EJ Clemson University: Clemson, SC, USA, 2011.

23. Chaskalovic, J. Mathematical and Numerical Methods for Partial Differential Equations; Springer: Berlin/Heidelberg, Germany, 2014; pp. 63-109.

24. Turek, S.; Hron, J.; Razzaq, M.; Wobker, H.; Schäfer, M. Numerical benchmarking of fluid-structure interaction: A comparison of different discretization and solution approaches. Fluid-Structure Interaction II: Modelling, Simulation, Optimisation; Bungartz, H.-J.; Mehl, M.; Schäfer, M. Volume 73 of Lecture Notes in Computational Science and Engineering Springer: Berlin/Heidelberg, Germany, 2010; pp. 413-424.

25. Hron, J.; Turek, S. A Monolithic FEM/Multigrid Solver for ALE Formulation of Fluid-Structure Interaction with Application in Biomechanics; Lecture Notes in Computational Science and Engineering Springer: Berlin/Heidelberg, Germany, 2006; Volume 53, pp. 146-170.

26. Turek, S.; Rivkind, L.; Hron, J.; Glowinski, R. Numerical study of a modified time-steping θ-scheme for incompressible flow simulations. J. Sci. Comput.; 2006; 28, pp. 533-547. [DOI: https://dx.doi.org/10.1007/s10915-006-9083-y]

27. Turek, S. Efficient Solvers for Incompressible Flow Problems, An Algorithmic and Computational Approach; Lecture Notes in Computational Science and Engineering Springer: Berlin/Heidelberg, Germany, 1999; Volume 6.

28. Amestoy, P.; L’Excellent, J.-Y.; Mary, T.; Puglisi, C. MUMPS: MUltifrontal Massively Parallel Solver for the direct solution of sparse linear equations. Proceedings of the EMS/ECMI Lanczos Prize, ECM Conference; Sevilla, Spain, 15–19 July 2024.

© 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

A Phenomenological Fluid–Structure Interaction Study of Plaque Rupture in Stenosed Bifurcated Elastic Arteries

Content area

Abstract

Full text

2. Bifurcated Stenosis Fluid-Structure Interaction Problem

3. Finite Element Method (FEM) Discretization

4. Solution Algorithm

5. Problem Configuration

6.1. Velocity Profile

6.2. Blood Backflow

6.3. Higher Reynolds Number

6.4. Dynamic Energized Hotspots

6.5. Shear Stress

6.6. Wall Shear Stress

6.7. Wall Displacement

7. Summary and Conclusions