A Fluid–Structure Interaction Analysis to Investigate the Influence of Magnetic Fields on Plaque Growth in Stenotic Bifurcated Arteries

Abstract

A finite element method is employed to examine the impact of a magnetic field on the development of plaque in an artery with stenotic bifurcation. Consistent with existing literature, blood flow is characterized as a Newtonian fluid that is stable, incompressible, biomagnetic, and laminar. Additionally, it is assumed that the arterial wall is linearly elastic throughout. The hemodynamic flow within a bifurcated artery, influenced by an asymmetric magnetic field, is described using the arbitrary Lagrangian–Eulerian (ALE) method. This technique incorporates the fluid–structure interaction coupling. The nonlinear system of partial differential equations is discretized using a stable P2P1 finite element pair. To solve the resulting nonlinear algebraic equation system, the Newton-Raphson method is employed. Magnetic fields are numerically modeled, and the resulting displacement, velocity magnitude, pressure, and wall shear stresses are analyzed across a range of Reynolds numbers (Re = 500, 1000, 1500, and 2000). The numerical analysis reveals that the presence of a magnetic field significantly impacts both the displacement magnitude and the flow velocity. In fact, introducing a magnetic field leads to reduced flow separation, an expanded recirculation area near the stenosis, as well as an increase in wall shear stress.

1. Introduction

Recently, there has been a notable surge in interest in looking to employ biomagnetic fluid flow as a therapy for illnesses such as atherosclerosis and other cardiovascular disorders. Researchers have been motivated to explore the topic of fluid flow analysis in numerous fields due to its significant importance. Atherosclerosis, once considered exclusively a cardiovascular disease, is now recognized as the buildup of fatty deposits and other substances in the walls of veins, which might potentially impede blood flow. Cardiovascular illnesses, such as stroke and heart disease, continue to be the primary causes of death in developed countries [1]. In this framework, particular attention is deserved to the rigidity of blood vessels [2,3] as well as the junction of the carotid veins situated on the outer edge of the blood vessels, which is the site of the most significant active infection at the moment; as a result, atherosclerosis is an important condition that is closely associated with this location [4]. It is widely recognized that two carotid veins, situated in a substantial region of the neck, convey oxygenated blood from the brain to the heart. In the literature, oxygenated blood flow is modeled as a Newtonian fluid that is stable, incompressible, biomagnetic, and laminar [5]. In [6], an effective use of the numerical research of the fluid–structure interaction (FSI) is made to enhance one’s understanding of the existing relationship between hemodynamic effects and atherosclerotic infection. This is accomplished due to the utilization of computational fluid dynamics, which not only replicates the flow of blood but also analyzes a study of the tension in the membrane.

The strain associated with the intense activity is the focus of the research presented in [7]. Using a three-dimensional FSI model of human atherosclerotic plaque to investigate the effect of plaque structure and material plaque characteristics on stress distribution, native improvements in stress/strain of plaque rupture, plaque area, and plaque rupture hazard estimation are reported in [8].

Many studies have been conducted to investigate the presence of bypass [9] or even more to analyze the association between the evolution of cardiovascular disorders and the dispersion of wall shear stress into the fluid domain [10,11,12]. The presence of the carotid sinus increases the risk of having a stroke [13] because it makes recirculation easier and lowers the wall shear stress. A larger carotid sinus is associated with increases in the OSI (oscillatory shear index) area, which is a risk factor for atherosclerosis [14,15]. A high OSI is frequently assumed near bifurcations [16]. When compared to the low-flared carotid sinus, the high-flared carotid sinus is more likely to be affected by atherosclerosis [17]. When considering stenosis, it is usual practice to make use of a magnetic field, see [18,19], since magnetic fields can impact rheological models of blood control and flow, and blood flow can be useful in some cases of hypertension.

The finite difference method is employed to examine the non-isothermal flow in stenosed arteries under the influence of a magnetic field. The findings provide strong evidence that the presence of particles in arteries can cause changes in both the temperature of the blood and the pattern of blood flow. Additionally, the magnetic field amplifies the flow characteristics of multiple narrowings in the arteries.

In general, blood behaves as non-Newtonian since the viscosity varies with the shear rate; nevertheless, the blood behaves as Newtonian when the shear rate is above 100 s⁻¹, and “the instantaneous shear rate over a cardiac cycle varies from zero to approximately 100 s⁻¹ in several large arterial vessels” [20]. In [20,21], a wider analysis to compare the various blood models is presented; in [20], the authors conclude that the Newtonian model of blood viscosity is a good approximation in regions of mid-range to high shear. In [22], the authors considered the problem without taking into account the elastic wall behavior. In the present paper, a Newtonian, viscous, steady, incompressible, laminar, biomagnetic fluid passing through a two-dimensional stenosis bifurcated artery is numerically investigated. As a novelty, we incorporated the elastic wall model into the set of equations derived and verified in [23] and simulated the corresponding system of equations for the elastic material under the effects of the considered fluid and magnetic effects. Moreover, the solution to the problem is investigated numerically by adopting the ALE technique along with strong coupling approaches, which are validated against the FSI benchmark problem, as detailed in [24,25].

The potential methods for numerically solving the coupled fluid–structure interaction (FSI) problem encompass solutions involving both the fluid and structural components through either partitioned or monolithic methodologies. Combining these elements often entails adopting a mixed representation, commonly referred to as the arbitrary Lagrangian–Eulerian (ALE) description, introducing added nonlinearities to the resultant equations. Constructing tools to effectively model diverse FSI scenarios remains a formidable task.

Approaches to solving fluid–structure interaction problems can be categorized into weakly (separately) and strongly coupled strategies, contingent on their manner of data interchange. In the separated coupling approach, the interconnected problem is divided into fluid and structural segments, each addressed independently. The concept of partitioned analysis for coupled systems was introduced by Park and Felippa [25].

In contrast, strongly coupled strategies entail the simultaneous monolithic solution of both FSI components. This involves formulating a singular set of equations upon discretizing the governing fluid and structural equations while incorporating interface boundary conditions. Consequently, the comprehensive FSI challenge is resolved using a unified, monolithic ALE framework.

Clearly, each approach carries its own merits and demerits, primarily related to flexibility, stability, and programming robustness. Prioritizing stability and robustness, we favor the monolithic ALE-FSI modeling approach [26] over partitioned methodologies, which offer a certain degree of flexibility [27].

A universal FSI problem comprises fluid description, solid description, interface conditions, and constraints for remaining boundaries. In this article, the incompressible Newtonian fluid flow interacts with a vertical elastic flap and elastic boundary are taken. Moreover, by employing the ALE technique, our study goes beyond a simplified one-way interaction model and achieves a more realistic representation of the complexities inherent in fluid–structure interaction. This choice aligns with our goal of providing a comprehensive and accurate analysis of the problem, enabling us to uncover insights that might otherwise be overlooked.

In summary, the ALE-FSI technique [25,26] was chosen for its capacity to capture large deformations, dynamic mesh movement, and accurate interface tracking, all of which contribute to a deeper understanding of fluid–structure interaction. Its advantages lie in its stability, accuracy, and real-world applicability, making it a pivotal tool in addressing the complexities of our study.

In the present paper, a study on the impact of magnetic fields on plaque growth in stenotic bifurcated arteries is presented. A novel approach integrating fluid–structure interaction analysis with computational modeling techniques is used. By considering the influence of magnetic fields on blood flow characteristics and wall tension, the research provides valuable insights into the biomechanical factors affecting plaque development.

2. FSI Problem Formulation

2.1. Fluid Model

In the present paper, the velocity $v^f$ and pressure p^f refer to fluid, as fluid is supposed to be Newtonian with constant density ${\rho }^f$ and dynamic viscosity $\mu$. Then the balance equation takes form as follows [27]:

$${\sigma }^f={\rho }^f\frac{𝜕v^f}{𝜕t}+{\rho }^f\nabla v^f\left(v^f-\frac{𝜕u}{𝜕t}\right), \mathrm{in} {\mathsf{\Omega }}_t^f,$$

(1)

$$v^f=0, \mathrm{in} {{\mathsf{\Omega }}_t^f}^{\prime }.$$

(2)

In Equation (1), the Lagrangian multiplier corresponding to the incompressibility constraint in (1) is denoted by ρ^f and ${\sigma }^f$ is a Cauchy stress tensor. To gain a deeper comprehension of the relationship, the existing denotes the stress tensor. For the solution of the balance equation, the basic relation for the stress tensor has to be introduced. A Newtonian fluid with constant density is being used here, yielding the following:

$${\sigma }^f={-p}^{f}I+\mu \left(\nabla v^f+{\left(\nabla v^f\right)}^T\right),$$

(3)

where μ is the dynamic viscosity and I is the identity matrix. The negative term −p^f I is the inviscid reactive part of the Cauchy stress tensor.

The magnetic field B is typically expressed in vector form as B = Bx i ˆ + By j ˆ + Bz k ˆ, where Bx, By, and Bz represent the components of the magnetic field in the x, y, and z directions, respectively. The unit vectors i ˆ, j ˆ, and k ˆ correspond to the directions of the x, y, and z axes, respectively.

Under the magnetohydrodynamic (MHD) approximation, the electric field can be considered insignificant compared to the field created by motion [28]. Therefore, we obtain the expression $\gamma =B^2\mathbf{u}$.

2.2. Structure Model

Let us assume that a solid object is elastic and has a bearing variable density. Let u^s be the displacement, and let $v^s$ be its velocity. Then, the balance equation is given by the following:

$${\rho }^s\frac{𝜕v^s}{𝜕t}+{\rho }^s\left(\nabla v^s\right)v^s-{\sigma }^s=0, \mathrm{in} {\mathsf{\Omega }}_t^s.$$

(4)

The above equation can be written as a Lagrangian description with respect to a certain initial state ${\mathsf{\Omega }}^s$ [27], thus yielding the following:

$${\rho }^s\frac{{𝜕}^2{u}^s}{𝜕t^2}-\left(J{\sigma }^s{F}^{-T}\right)=0, \mathrm{in} {\mathsf{\Omega }}^s.$$

(5)

In Equation (5), $F^{-T}$ = (F − 1)^T, where F − 1 is the inverse deformation gradient, and J = det F. The basic equations for the stress tensor for compressible and incompressible structures can be dealt with in the same manner. The undeformed structure density is ${\rho }^s$, while the Poisson ratio ${\nu }^s$, the Lame coefficient $\mathsf{\lambda },$ and Young’s modulus E define the elasticity of structure as follows:

$${\nu }^s=\frac{{\lambda }^s}{2\left({\mu }^s+{\lambda }^s\right)}; E=\frac{{\mu }^s\left(2{\mu }^2+3{\lambda }^s\right)}{{\mu }^s+{\lambda }^s}, {\mu }^s=\frac{E}{2\left({\nu }^s+1\right)}; {\lambda }^s=\frac{E{\nu }^s}{\left({\nu }^s+1\right)\left(1-2{\nu }^s\right)}.$$

(6)

For an incompressible structure, the Poisson ratio is 1/2, while for a compressible structure, it is less than 1/2. The basic relation between stress and strain depends upon the 2nd Piola–Kirchhoff stress tensor S and the Green Lagrangian strain tensor E as a function of the Green Lagrangian strain tensor. The 2nd Piola–Kirchhoff stress tensor holz may be derived from the Cauchy stress tensor ${\sigma }^s$ as follows:

$$S^s=J{F}^{-1}{\sigma }^s{F}^{-T},$$

(7)

and E can be rewritten as follows:

$$E=\frac{1}{2}{F}^{T}F-I.$$

(8)

In the present paper, the Cauchy stress tensor for the St. Venant–Kirchhoff material model is defined as follows:

$${\sigma }^s=\frac{1}{J}F{S}^s{F}^T,{ S}^s={\lambda }^{s}tr\left(E\right)I+2{\mu }^{s}E,$$

(9)

where J represents the determinant of the deformation gradient tensor F, described as follows:

$$F=I+\nabla u^s.$$

(10)

2.3. ALE Form of Fluid–Structure Interaction Equations

A pseudo-solid mapping method is used to construct this fluid–structure interaction problem [27]. The complete dimensionless equation system with the mentioned material relation is given by the following:

$$\frac{𝜕u}{𝜕t}=\left\{\begin{array}{c}∆u\in {\mathsf{\Omega }}^f\\ v\in {\mathsf{\Omega }}^s\end{array}\right.$$

(11)

$$\frac{𝜕v}{𝜕t}=\left\{\begin{array}{c}-\left(v-\frac{𝜕u}{𝜕t}\right)F^{-1}\nabla v+\nabla ·\left(J{p}^f{F}^{-T}+J\mu \nabla \left(v{F}^{-1}{F}^{-T}\right)\right)\in {\mathsf{\Omega }}^f\\ \frac{1}{\beta }\nabla ·\left(-J{p}^s{F}^{-T}\right)\in {\mathsf{\Omega }}^s\end{array}\right.$$

(12)

$$0=\left\{\begin{array}{c}\nabla ·\left(Jv{F}^{-T}\right)\in {\mathsf{\Omega }}^f\\ J-1\in {\mathsf{\Omega }}^s\end{array}\right.$$

(13)

where β is simply the solid-to-fluid density ratio, i.e., β = ϱ^s/ϱ^f. The following Neumann and Dirichlet constraints are used to preserve no-slip conditions along the fluid structure interface:

$$v^f=v^s;{\sigma }^{f}n={\sigma }^{s}n .$$

(14)

3. FEM Discretization

A common choice of finite element space is the Taylor–Hood pair P2P1. The Taylor–Hood element is an optimal convergent, and this element satisfies the Ladyzhenskaya–Babuška–Brezzi (LBB) stability condition [29]. The Taylor–Hood P2P1 element is chosen for its accuracy, stability, conservation properties, and compatibility with solid mechanics elements, making it well-suited for studying the complexities of fluid–structure interaction in the specific FSI problem at hand. Its advantages align with the complexities of the problem, allowing for a more robust and insightful analysis of the fluid–structure system. Going from discontinuous space to continuous space, the dimension of pressure space is reduced, which saves the degrees of freedom (see Figure 1 for the degrees of freedom).

The basis function ϕ_i of the P1 element is given by the following [29]:

$${\varphi }_i\left(x\right)=a_i+b_{i}x+c_{i}y ,$$

(15)

and the basis vector set is <1, x, y>. For the P2 element, we have the following:

$${\varphi }_i\left(x\right)=a_i+b_{i}x+c_{i}y+d_i{x}^2+e_i{y}^2+f_{i}xy ,$$

(16)

and the basis vector set is <1, x, y, x², xy, y²>.

Weak Formulation

The standard Galerkin finite element technique is used for discretization in space on the time interval [0, T]. Equations (11)–(13) are multiplied by the test functions ζ, ξ, and γ and integrated over the space and the time interval [0, T], obtaining the following:

$${\int }_0^T{\int }_{\mathsf{\Omega }}\frac{𝜕u}{𝜕t}·\zeta dVdt={\int }_0^T{\int }_{{\mathsf{\Omega }}^s}v·\zeta dVdt-{\int }_0^T{\int }_{{\mathsf{\Omega }}^f}\nabla v·\nabla \zeta dVdt$$

(17)

$${\int }_0^T{\int }_{{\mathsf{\Omega }}^f}J\frac{𝜕v}{𝜕t}·\xi dVdt+{\int }_0^T{\int }_{{\mathsf{\Omega }}^s}\beta J\frac{𝜕v}{𝜕t}·\xi dVdt={\int }_0^T{\int }_{\mathsf{\Omega }}Jp{F}^{-T}·\nabla \xi dVdt$$

$$-{\int }_0^T{\int }_{{\mathsf{\Omega }}^f}J\nabla v{F}^{-1}\left(v-\frac{𝜕u}{𝜕t}\right)·\xi dVdt-{\int }_0^T{\int }_{{\mathsf{\Omega }}^f}J\mu \nabla v{F}^{-1}{F}^{-T}·\nabla \xi dVdt$$

(18)

and, finally:

$$0={\int }_0^T{\int }_{{\mathsf{\Omega }}^s}\left(J-1\right)\gamma dVdt+{\int }_0^T{\int }_{{\mathsf{\Omega }}^f}\nabla ·\left(Jv{F}^{-T}\right)\gamma dVdt.$$

(19)

In Equations (17)–(19), the integration by parts together with the boundary conditions have been applied. Indeed, after FEM space discretization, we arrive at a nonlinear algebraic system of equations:

$$\left(\begin{array}{ccc}{S}_{uu}& S_{uv}& 0\\ S_{vu}& S_{vv}& kB\\ c_u{B}_s^T& c_v{B}_f^T& 0\end{array}\right)\left(\begin{array}{c}{u}_h\\ v_h\\ p_h\end{array}\right)=\left(\begin{array}{c}rhsu\\ rhsv\\ rhsp\end{array}\right),$$

(20)

where S represents the diffusive, reactive, and convective operators from the governing equations, and B and $B^T$ are the discrete gradient and the discrete divergence operators, respectively.

4. Numerical Solution

The system of nonlinear algebraic equations of saddle point type in (20) is solved by using the Newton iteration, which can exhibit quadratic convergence for a sufficient close solution to find a root of the residual R(X) = 0, by using the already known function value and its first derivative. The formula for Newton iteration with damping reads as follows:

$${\mathit{X}}^{n+1}={\mathit{X}}^n+{\omega }^n{\left[\frac{𝜕\mathit{R}\left({\mathit{X}}^n\right)}{𝜕\mathit{X}}\right]}^{-1}\mathit{R}\left({\mathit{X}}^n\right),$$

(21)

where X$\left({\mathit{u}}_{\mathit{h}},{\mathit{v}}_{\mathit{h}},{\mathit{p}}_{\mathit{h}}\right).$ The Jacobian matrix $\frac{𝜕\mathit{R}\left({\mathit{X}}^{\mathit{n}}\right)}{𝜕\mathit{X}}$ is computed by taking the finite differences and the residual matrix as follows:

$${\left[\frac{𝜕\mathit{R}\left({\mathit{X}}^n\right)}{𝜕\mathit{X}}\right]}_{ij}\approx \frac{{\left[\mathit{R}\right]}_i\left({\mathit{X}}^n+{\alpha }_j{e}_j\right)-{\left[\mathit{R}\right]}_i\left({\mathit{X}}^n-{\alpha }_j{e}_j\right)}{2{\alpha }_j} ,$$

(22)

where the coefficients α_j > 0 are increments at each iteration step n, and e_j are the unit basis vectors in Rn of Equation (21) (as detailed in [29]). The parameter ωⁿ ∈ (−1, 0) is taken in such a way that the error measure decreases:

$$\mathit{R}\left({\mathit{X}}^{n+1}\right){\mathit{X}}^{n+1}\le \mathit{R}\left({\mathit{X}}^n\right){\mathit{X}}^n .$$

(23)

This damping greatly improves the robustness of the Newton iteration in the case when the current approximation ${\mathit{X}}^{\mathit{n}}$ is not close enough to the final solution, as particularly explained in [25]. In this considered two-dimensional problem, a direct solver for sparse systems, like MUMPS, is used. A description of the steps of the Newton iteration algorithm [30] is reported in

Table 1. Algorithm: Newton iteration and line search.

Step	Description
1	Nonlinear tolerance parameter input
2	Initialized n = 0. Take $X^n$ as a starting guess
3	Build the residual vector $R\left(X^n\right)=A{X}^n-b$
4	Compute the Jacobian matrix $\mathbf{J}\left(X^n\right)$ = $\mathsf{\partial }\mathbf{R}/\mathsf{\partial }\mathbf{X}\left(X^n\right)$
5	Go to the linear system for correction of $\mathsf{\delta }\mathbf{X} \mathrm{and} \mathrm{solve} \mathbf{J}\left(X^n\right)$$\mathsf{\delta }{X}^n$ = $\mathrm{R}\left(X^n\right)$
6	Find an optimal step length ${\omega }^n$ ∈ (−1, 0]
7	Update the solution $X^{n+1}$ = $X^n$ + ${\omega }^n$ $\mathsf{\delta }{X}^n$

.

Specifically, we specify the stopping criteria that determine the iteration process’s convergence. This includes detailing factors such as the maximum allowable number of iterations and the tolerance level for the residual norm, which is 10⁻⁶. Where the Newton iteration algorithm may fail to converge or encounter difficulties during the solution process. This may involve techniques like adaptive step sizes, modifying the initial guess, or refining the numerical discretization. By sharing these strategies, we hope to demonstrate the robustness of our solution approach and its potential applicability to a range of fluid–structure interaction problems.

5. Mathematical Model

To model the blood flow [5,21], let us refer to a biomagnetic fluid that moves through a stenosis in a bifurcated artery; it has a two-dimensional structure, is viscous, remains stable, and cannot be compressed. In addition, let us assume that the magnetic field has a minor impact on the viscosity of the fluid, even if this assumption deserves further investigation [31], that the walls of the geometry under investigation are insulated, and that the electric field is relatively weak. It is generally accepted that the walls of the arteries possess a degree of linear flexibility [32]. The problem geometry and a coarse mesh are shown in Figure 2.

A magnetic field of constant strength is applied along the Y axis. Both the inflow velocity profile and the outflow pressure are assumed to follow a parabolic shape. The outflow pressure is expected to remain unchanged at zero. The governing equations for flow in the aforementioned framework, according to the arbitrary Lagrangian–Eulerian (ALE) formulation, are as follows:

• Conservation of mass

$$\frac{𝜕u^{*}}{𝜕x^{*}}+\frac{𝜕v^{*}}{𝜕y^{*}}=0;$$

(24)

• Conservation of momentum

$$\frac{-𝜕p^{*}}{𝜕x^{*}}-{\sigma }^{*}{u}^{*}{B}^2+\mu {\nabla }^2{u}^{*}={\rho }_f\left(u^{*}-w^{*}\right)\frac{𝜕u^{*}}{𝜕x^{*}}+{\rho }_f\left(v^{*}-w^{*}\right)\frac{𝜕u^{*}}{𝜕y^{*}};$$

(25)

$$\frac{-𝜕p^{*}}{𝜕y^{*}}+\mu {\nabla }^2{v}^{*}={\rho }_f\left(u^{*}-w^{*}\right)\frac{𝜕v^{*}}{𝜕x^{*}}+{\rho }_f\left(v^{*}-w^{*}\right)\frac{𝜕v^{*}}{𝜕y^{*}};$$

(26)

In Equations (24)–(26), $w^{*}$ represents the velocity in the mesh coordinate system, while $u^{*}$ and $v^{*}$ are the velocity components in the various dimensions. The symbol ${\sigma }^{*}$ represents the electrical conductivity of the biomagnetic fluid. The amount ${\sigma }^{*}{u}^{*}{B}^2$ in Equation (25), generated by the fluid’s electric conductivity, represents the Lorentz force per unit volume.

The governing equation for dimensional solid displacement is as follows:

• Conservation of momentum

$$\nabla {\sigma }^{*}=-F_v^{*};$$

(27)

We define the strain tensor by assuming that the wall is linearly elastic, i.e.,

$${\sigma }^{*}=\frac{1}{J}FS{F}^T.$$

(28)

The relationship between the second Piola–Kirkhoff stress tensor S and the strain ϵ is described through the equation:

$$ϵ=\frac{1}{2}\left[{\left(\nabla u_s\right)}^T+\nabla u_s+{\left(\nabla u_s\right)}^T·\nabla u_s\right].$$

(29)

The first step towards the numerical solution of the system is to convert it into dimensionless form by introducing the following dimensionless variables and parameters:

$$\left(x,y\right)=\left(\frac{x^{*}}{h},\frac{y^{*}}{h}\right);\left(u,v,w\right)=\left(\frac{u^{*}}{u_0},\frac{v^{*}}{u_0},\frac{w^{*}}{u_0}\right);p=\frac{p^{*}}{{\rho }_f{u}_0^2}; \mathrm{R}\mathrm{e}=\frac{h{\rho }_f{u}_0}{\mu };\mathrm{H}\mathrm{a}=\frac{{\sigma }^{f*}{h}^2{B}^2}{\mu };$$

(30)

where h is the shortest distance between two walls of the artery and $u_0$ is the maximum velocity of the blood at the inlet, Re is the Reynolds number, and Ha is the Hartmann number.

On account of the dimensionless quantities (30), Equations (24)–(28) can be rewritten as follows:

$$\frac{𝜕u}{𝜕x}+\frac{𝜕v}{𝜕y}=0;$$

(31)

$$\frac{-𝜕p}{𝜕x}-\frac{{Ha}^2}{Re}u+\frac{1}{Re}{\nabla }^{2}u=\left(u-w\right)\frac{𝜕u}{𝜕x}+\left(v-w\right)\frac{𝜕u}{𝜕y} ;$$

(32)

$$\frac{-𝜕p}{𝜕y}+\frac{1}{Re}{\nabla }^{2}v=\left(u-w\right)\frac{𝜕v}{𝜕x}+\left(v-w\right)\frac{𝜕v}{𝜕y};$$

(33)

$$\nabla \sigma =-F_v.$$

(34)

Problem Description

The mesh reported in Figure 2 represents the typical two-dimensional representation of the model. The computational domain can be thought of as a symmetric bifurcation or stenosis. An isotropic linear elastic material with a known Young’s modulus and Poisson ratio is used to construct the walls. The Lame coefficient, λ, and the shear modulus, μ, are therefore used to define them. They take on the following responsibilities as a group:

$$\nu =\frac{\lambda }{2\left(\mu +\lambda \right)}; E=\frac{\mu \left(2\mu +3\lambda \right)}{\mu +\lambda };\mu =\frac{E}{2\left(\nu +1\right)} ; \lambda =\frac{E\nu }{\left(\nu +1\right)\left(1-2\nu \right)}.$$

(35)

Materials presenting ν < 0.5 are compressible, while materials having ν = 0.5 are incompressible [33]. The parameters Young’s modulus E = 5 × 10⁵ and Poisson ratio ν = 0.49 are utilized in the present analysis.

Figure 2 shows a computational domain that has stenosis. The dimensions of this domain are L = 6.34, the diameter of the parent artery is h = 1, and the diameter of the stenosis site is reduced by 50%. The diameter of a daughter arterial is h₁ = 0.37, and w = 0.08 is the elastic wall width of this artery. We determine that the angle of the bifurcation artery is 37°. At the point where t = 0 and t = 1, the beginning of line A can be found at the coordinates (1.8, t). If line C is the center line along which pressure measurements are taken, then the distance d between line A and the stenosis is calculated to be d = 0.8. The coarse mesh that is displayed in Figure 2 contains a total of 1478 domain elements in addition to 338 boundary elements. In

Table 2. Mesh refinement levels, depending on the number of components and wall shear loads.

Mesh Level	Mesh Elements	Wall Shear Stress	Abs Error
0	1478	0.072341	-
1	1922	0.072244	0.000097
2	2432	0.072347	0.000103
3	7528	0.072557	0.000210
4	22153	0.072601	0.000040
5	24322	0.072600	0.000001

, we can see the absolute error of the WSS as a function of the mesh refinement level and the number of elements. We can also see the WSS on the upper elastic wall in this table. In the computation for Level 1, we use a total of 1922 finite elements.

The inlet velocity profile is assumed to have a parabolic shape, namely:

$$u_f\left(x,y\right)=2y\left(1-y\right).$$

(36)

The pressure determines the parameters for the circumstances of the outflow borders. At the point where the pressure is discharged, it is equal to zero.

The comparison for the velocity magnitude between FSI and CFD cases is shown in

Table 3. Comparison of FSI and CFD cases for the velocity magnitude for different values of Ha and fixed values of Re.

	Re = 500	Re = 500	Re = 1000	Re = 1000
Ha	FSI	CFD	FSI	CFD
0	0.5595	0.5846	0.5774	0.5959
8	0.5443	0.5884	0.5624	0.5914
10	0.5358	0.5783	0.5504	0.5865
12	0.5169	0.5709	0.5424	0.5825

. This numerical comparison is taken for the varying values of Ha for fixed Re. The table shows that the numerical results are in good agreement with the existing experimental and benchmark results [34,35,36].

6. Results and Discussion

In order to solve the formulated equation system, the ALE approach is used, paired with appropriate boundary conditions. This produces physical insight into the problem. In order to convert the governing nonlinear PDEs into a dimensionless form, well-known parameters such as the Reynolds number and the Hartmann number are utilized, in the range 500 to 2000 for the Reynolds number and 0 to 12 for the Hartmann number, respectively, for the fluid flow simulation within the simulation domain. In detail, this section is organized as follows. First, the velocity profiles in the stenosis are discussed as a function of Re and Ha. Then, attention is paid to the total displacement of the upper wall of the bifurcated artery due to the intensity of the magnetic field. By simplifying artery walls as elastic, the analysis of blood flow can become more manageable. This reduction facilitates comprehension of fundamental flow phenomena such as velocity profiles, shear loads, and potential locations for flow separation. Elastic models provide precise computations of stress and strain distributions inside the artery wall, which are crucial for comprehending mechanical stability and identifying possible locations of failure.

Section 6.2 highlights the velocity magnitude profile both before and after the stenosis at points A and B as an additional measure for the analysis of the primary impact that the stenosis has on the velocity of the fluid flow. In Section 6.3 and Section 6.4, attention is paid to the pressure field and to the wall shear stress, respectively.

In Figure 3, we show how a symmetric, stenosed bifurcated artery reacts to an asymmetric magnetic field with linearly elastic walls. This reaction is assessed by the behavior of the dimensionless velocity distribution of the artery’s blood flow.

In an entirely hydrodynamic situation, with Ha = 0, the maximum velocity is reached. Almost immediately after the stenosis, the flow separates, which is a strong indication that the wall shear stress has diminished. A little depression emerges in the region where the two flows are separated. This can lead to the development of atherosclerosis (plaque-induced artery wall degeneration). Because of the walls’ elasticity, it was discovered that there was very little deformation in the structure. When an external magnetic field is introduced, the flow separation area becomes less extensive and the velocity of the artery slows down, but this latter change is largely inconsequential. This results in the hollow becoming less noticeable. As the strength of the magnetic field increases, the zone of separation shrinks to a much more manageable size. The condition known as thrombosis can be remedied with this procedure (a medical term for a blood clot). When MHD flow is present, one might see extreme wall deformation, which ultimately leads to the creation of the WSS. Due to the predominance of viscous forces, significant wall deformation occurs at a low Reynolds number, such as Re = 300, for the magnetic parametric values Ha examined.

Figure 4 illustrates the velocity curve throughout the artery at a Reynolds number of 1000. A cavity that might be considered stable was produced in the case where hydrodynamic forces were absent (Ha = 0). Under the influence of a magnetic field, on the other hand, the cavity will contract, resulting in a smaller overall size. Since the value of the magnetic parameter was increased from 8 to 10 and then to 12, the cavity was no longer there. A rise in the magnetic field parameter Ha creates a significant degree of wall deformation, which in turn results in a reduction in the amount of flow separation.

Figure 5 depicts the patterns that emerge when velocity is plotted against Re = 1500. Shortly after the stenosis had developed at position 11, this massive recirculation is easily observed in the pure hydrodynamic scenario, where Ha = 0. When Reynolds is equal to 1500, the size of the cavity steadily shrinks when the value of Ha is increased to greater levels. This indicates that the Reynolds number, which equals 1500, is directly proportional to the amount of recirculation, whereas the Hartmann number, which equals n, is inversely related to the amount of recirculation.

Figure 6 highlights the increase in arterial velocity that occurs after the value of Re = 2000 has been applied. When Ha = 0, a scenario characterized as being fully hydrodynamic, a vast, extended hollow is produced. The cavity, on the other hand, contracted when a magnetic field was injected. Even though the value of the magnetic parameter is increased to Ha = 8, 10, and 12, the cavity continues to exist because of the huge value of Re = 2000 in comparison to the cases, which have Re values of 500 and 1000, respectively. When the parameter Ha, which represents the magnetic field, approaches the values 8, 10, and 12, the cavity area decreases, and the wall deformation is reduced to a minimum.

6.1. Displacement

Figure 7 illustrates the total displacement of the upper wall (the same will be for the lower wall due to the symmetric nature of the geometry) of the bifurcated artery under the effect of the Hartmann number, i.e., Ha = 0, 8, 10, and 12. A decreasing trend is evident in the total displacement with an increase in Reynolds number, Re. The maximum displacement occurs immediately before the stenosis, and it begins to decrease after the stenosis.

6.2. Velocity Profile at Positions A and B

In this section, we take a closer look at the behavior of the velocity profile at positions A and B. When compared side by side, the relative speeds of point A and point B reveal a stunning discrepancy. When the Reynolds number is held constant, the relationship between velocity and the inverse of the magnetic field holds, whereas the relationship between the magnetic field and the velocity profile holds the opposite.

The velocity magnitude profiles at locations A and B are reported in Figure 8. Moreover,

Table 4. Maximum velocity magnitude at locations A and B and gain in velocity at B.

		Re = 500			Re = 1000
Ha	A	B	Gain %	A	B	Gain %
0	0.49	0.69	41	0.50	0.70	40
8	0.42	0.62	24	0.45	0.61	36
10	0.37	0.56	23	0.42	0.56	33
12	0.34	0.40	18	0.37	0.51	38
		Re = 1500			Re = 2000
Ha	A	B	Gain %	A	B	Gain %
0	0.50	0.71	42	0.51	0.72	39
8	0.46	0.65	41	0.47	0.66	40
10	0.45	0.61	36	0.46	0.64	39
12	0.42	0.57	36	0.43	0.60	40

shows the maximum velocity before stenosis at location A and after stenosis at location B. The percent increment is also given for different values of the Reynolds number.

6.3. Pressure

In Figure 9, pressure attitude is analyzed versus the magnetic field parameter, i.e., Ha = 0, 8, 10, and 12, and values are plotted for different Reynolds numbers, i.e., Re = 500, 1000, 1500, and 2000. The pure hydrodynamic case, i.e., Ha = 0, yields the minimum pressure. The magnetic field parameter Ha displays a direct relation with the pressure and an inverse relation with the Reynolds number Re. Furthermore, Re = 500 represents the maximum pressure for a higher value of the magnetic field parameter, i.e., Ha = 8, 10, and 12. Theorem-type environments (including propositions, lemmas, corollaries, etc.) can be formatted as follows:

6.4. Wall Shear Stress

The total value of the wall shear stress along the lower wall of stenosed bifurcation is calculated and listed in

Table 5. Wall shear stress along the lower wall as a function of Re and Ha.

Ha	Re = 500	Re = 1000	Re = 1500	Re = 2000
0	0.05476	0.03018	0.02161	0.01716
8	0.06673	0.03522	0.02457	0.01918
10	0.07223	0.03807	0.02623	0.02030
12	0.07776	0.04128	0.02824	0.02165

under the influence of the magnetic parameter Ha = 0, 8, 10, and 12, and for the different values of the Reynolds numbers, i.e., Re = 500, 1000, 1500, and 2000, respectively.

Furthermore, results depict that wall shear stress is directly influenced by the magnetic parameter, Ha, and has an inverse relation with the Reynolds numbers, Re. An increasing behavior can be observed for the small values of the Reynolds number, Re, and at the higher value of the Hartmann number, Ha. An artery with a small Reynolds number flow, i.e., Re = 500, may experience reasonable stress in the presence of a higher value of the magnetic field parameter, i.e., Ha = 12, as shown in Figure 10.

Wall shear stress is closely associated with the growth of atherosclerosis. The wall shear stress affects the endothelial cell, which causes a change in the gene pattern [37]. Plaque formation occurs due to wall shear stress or endothelial cells. Low wall shear stress is a relevant factor for the formation of the plaque. Usually, atherosclerosis appears in the bifurcated area. Here, it is evident that the values of the wall shear stress remain the same for the upper wall and lower wall, due to the symmetric nature of the considered geometry.

There are some findings from the numerical analysis, as follows:

• The presence of a magnetic field, described through the parameter Ha, in increasing order, reduces the size of the cavity in the same order adjacent to stenosis.
• The presence of a magnetic field decreases the gain in velocity after stenosis as compared to a purely hydrodynamic case, i.e., Ha = 0. For instance, in the case of Re = 500, when Ha = 0, velocity gain is 41%, and when Ha = 12, velocity gain is only 18%.
• The pressure becomes higher with the increase in Hartmann number values, i.e., Ha = 0, 8, 10, 12.
• A flow with a small Reynolds number, i.e., Re = 500, and a higher value of the magnetic field parameter, i.e., Ha = 12, gives significant rise to wall shear stress.
• The wall displacement is large when Re = 500 and Ha = 12, and the maximum wall displacement is observed just before the stenosis.

7. Conclusions

In the present paper, we investigate the impact of magnetic fields on plaque growth in stenotic bifurcated arteries. The novelty is given by the approach, i.e., an integration of fluid–structure interaction analysis with computational modeling techniques. By considering the influence of magnetic fields on blood flow characteristics and wall tension, we delve into the biomechanical factors affecting plaque development.

The blood is modeled as a biomagnetic fluid that is two-dimensional, stable, incompressible, and laminar as it travels through an artery that has symmetric stenosis and elastic wall surfaces. An arbitrary Lagrangian–Eulerian method is employed in conjunction with a bi-quadratic linear stable finite element pair to solve the governing nonlinear partial differential equations of fluid and solid (elasticity). The Newton-Raphson method is utilized to find a solution to the nonlinear algebraic problem obtained.

The analysis of stress distribution, plaque rupture risk, and wall shear stress highlights the importance of considering multiple factors in assessing cardiovascular health outcomes.

In order to analyze the main features of the obtained solutions, charts and tables are presented and discussed. Overall, the comparison with the existing literature confirms the novelty of the results, validating the choice of numerical methods and convergence strategies. The findings not only advance the understanding of FSI problems with complex geometries [38], but also contribute to the development of efficient and reliable numerical techniques in this critical field of engineering research.

In conclusion, this numerical analysis suggests that the magnetic field has crucial managerial behavior patterns toward blood flow. As a result, it is possible that the magnetic field could be implemented in a wide variety of treatments for conditions such as plaque, atheroma, hypertension, blood pressure, and other similar conditions. Further investigations will focus on modeling the blood flow as non-Newtonian, investigating time periodicity, and considering non-symmetric domains. Regarding the inlet velocity profile, the effects of different inlet conditions may be further investigated [39,40]. Moreover, a comparison with experimental results will further validate this analysis.