Fluid-Structure Interaction: Numerical Analysis of Biomagnetic Flow Inhibition on a Plaque in a Stenosed Bifurcation Artery

1. Introduction

The use of biomagnetic fluid flow as a treatment for atherosclerosis and other cardiovascular disorders has garnered interest in recent years [1,2]. Atherosclerosis is the accumulation of fatty acids, insulation, and other surpluses on the walls of veins, which can decrease blood flow and is a major cause of death in industrialized nations [1,2]. While the early stages of atherosclerosis may not have visible impacts, advanced stages can hinder blood flow and cause ischemia [3,4,5]. Atherosclerosis is closely associated with the carotid arteries, which are important for transporting oxygenated blood from the brain to the heart [6,7,8,9,10,11]. Computational fluid dynamics and numerical research of fluid-structure interaction have been used to study the relationship between hemodynamic effects and atherosclerotic infection [12,13,14].

The presence of a carotid sinus is associated with a higher risk of stroke due to easier recirculation and lower wall shear stress [15]. The size of the carotid sinus is a risk factor for atherosclerosis, with larger sinuses associated with higher Oscillatory Shear Index (OSI) [16,17]. High OSI is commonly found near bifurcations and can lead to atherosclerosis [18]. The use of a magnetic field can impact rheological models of blood control and flow, and has been shown to be useful in cases of hypertension [31]. A finite difference approach has been used to analyze non-isothermal flow in stenosed arteries exposed to a magnetic field [31]. Turbulent magnetohydrodynamic behavior of blood as it flows through narrowing arteries has also been investigated [32,33,34,35].

In this paper, we investigate the behavior of a Newtonian, viscous, steady, incompressible, laminar biomagnetic fluid passing through a two-dimensional stenotic bifurcated artery, incorporating an elastic wall model into the set of equations [36,37]. The problem is solved numerically using the ALE technique and strong coupling approaches, which are validated against the FSI-benchmark problem [38,40].

The upcoming content will follow this layout: Section 5 will focus on mathematical modeling and problem formulation, while Section 6 will provide a brief overview of the results and discussion. The conclusion will be presented in Section 7.

2. Problem Formulation

A typical Fluid-Structure Interaction (FSI) problem involves modeling the fluid and solid domains, imposing interface conditions, and enforcing boundary constraints. This study focuses on the interaction of an incompressible Newtonian fluid with a vertical elastic flap and elastic boundary.

2.1. Fluid Model

The characteristics of a Newtonian fluid are represented by the velocity ${\mathbf{v}}^f$ and pressure$p^f$, where the fluid’s density is denoted by ${\rho }^f$, and its kinematic viscosity is symbolized by ${\nu }^f$. The balance equation is then expressed as follows [40]:

$$\begin{array}{cc}\hfill \mathrm{div}{\mathsf{œ}}^f& ={\varrho }^f\frac{\partial {\mathbf{v}}^f}{\partial t}+{\varrho }^f(\nabla {\mathbf{v}}^f)({\mathbf{v}}^f-\frac{\partial u^f}{\partial t})\mathrm{in}{\mathrm{\Omega }}_t^f,\hfill \\ \hfill 0& =div{\mathbf{v}}^f\mathrm{in}{\mathrm{\Omega }}_t^f\hfill \end{array}$$

(1)

To solve the balance equation, it is important to establish the fundamental relationship for the stress tensor. In this case, a Newtonian fluid with constant density is being utilized:

$$\begin{array}{c}\hfill {\mathsf{œ}}^f=-p^f\mathbf{I}+\mu (\nabla {\mathbf{v}}^f+{(\nabla {\mathbf{v}}^f)}^T),\end{array}$$

(2)

the constant dynamic viscosity of the Newtonian fluid is represented by $\mu$, and the identity matrix is denoted by $\mathbf{I}$. In addition, the Lagrangian multiplier related to the incompressibility constraint in (1) is identified as $p^f$, and the stress tensor is denoted by ${\mathsf{œ}}^f$. The term $-p^f\mathbf{I}$ represents the inviscid reactive component of the Cauchy stress tensor.

2.2. Structure Model

Structure properties

Suppose the solid object has an elastic property and its density varies. Let ${\mathbf{u}}^{\mathbf{s}}$ represent its displacement and ${\mathbf{v}}^{\mathbf{s}}$ represent its velocity. The balance equation can then be expressed as follows:

$${\rho }^s\frac{\partial {\mathbf{v}}^s}{\partial t}+{\rho }^s(\nabla {\mathbf{v}}^s){\mathbf{v}}^s-\mathrm{div}{\sigma }^s=0+{\rho }^s\mathbf{g},\mathrm{in}{\mathrm{\Omega }}_t^s.$$

(3)

The above equation can be expressed in a Lagrangian framework with respect to a specific initial state ${\mathrm{\Omega }}^s$. More information can be found in [40]

$${\rho }^s\frac{{\partial }^2{\mathbf{u}}^{\mathbf{s}}}{\partial t^2}-\mathrm{div}\left(J{\sigma }^s{F}^{-T}\right)=0+{\rho }^s\mathbf{g},\mathrm{in}{\mathrm{\Omega }}^s.$$

(4)

Here, $F^{-T}={\left(F^{-1}\right)}^T$, $F^{-1}$ represents the inverse deformation gradient, and $J=detF$. The equations for the stress tensor can be applied to both compressible and incompressible structures using the same approach. The undeformed structure density is denoted by ${\rho }^s$, while the elasticity of the structure is defined by the Poisson ratio ${\nu }^s$ and Young’s modulus E.

$${\nu }^s=\frac{{\lambda }^s}{2\left({{\mu }^s+\lambda }^s\right)}E=\frac{{\mu }^s(2{\mu }^2+3{\lambda }^s)}{({\lambda }^s+{\mu }^s)}$$

(5)

$${\mu }^s=\frac{E}{2({\nu }^s+1)}{\lambda }^s=\frac{E{\nu }^s}{({\nu }^s+1)(1-2{\nu }^s)},$$

(6)

The Poisson ratio is $1/2$ for incompressible structures, while it is less than $1/2$ for compressible structures. The relationship between stress and strain is expressed in terms of the $2^{\mathrm{nd}}$ Piola-Kirchhoff stress tensor S and the Green-Lagrangian strain tensor E. The $2^{2\mathrm{nd}}$ Piola-Kirchhoff stress tensor can be derived from the Cauchy stress tensor ${\sigma }^s$, as

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

(7)

and E can be written as

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

(8)

In this present work, the Cauchy stress tensor ${\sigma }^s$ for the St.Venant-Kirchhoff material model is defined as follows

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

(9)

where, J represent the determinant of the deformation gradient tensor F, described by

$$F=I+\nabla {\mathbf{u}}^{\mathbf{s}}.$$

(10)

2.3. Fluid Structure Interaction Set of Equations

The technique of constructing the fluid structure interaction problem is through the application of a pseudo-solid mapping method. The dimensionless system, which incorporates the mentioned material relation, is presented in the following manner:

$$\begin{array}{cc}\hfill \frac{\partial \mathbf{u}}{\partial t}=& \left\{\begin{array}{cc}\mathrm{\Delta }\mathbf{u}\hfill & \mathrm{in}{\mathrm{\Omega }}^f,\hfill \\ \mathbf{v}\hfill & \mathrm{in}{\mathrm{\Omega }}^s,\hfill \end{array}\right.\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \frac{\partial \mathbf{v}}{\partial t}=& \left\{\begin{array}{c}-(\mathbf{v}-\frac{\partial \mathbf{u}}{\partial t}){\mathbf{F}}^{-1}\left(\mathrm{Grad}\mathbf{v}\right)\hfill \\ +\mathrm{Div}\left(-J{p}^f{\mathbf{F}}^{-T}+J\mu \mathrm{Grad}\mathbf{v}{\mathbf{F}}^{-1}{\mathbf{F}}^{-T}\right)\hfill & \mathrm{in}{\mathrm{\Omega }}^f,\hfill \\ \frac{1}{\beta }\mathrm{Div}\left(-J{p}^s{\mathbf{F}}^{-T}\right)\hfill & \mathrm{in}{\mathrm{\Omega }}^s,\hfill \end{array}\right.\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill 0=& \left\{\begin{array}{cc}\mathrm{Div}\left(J\mathbf{v}{\mathbf{F}}^{-T}\right)\hfill & \mathrm{in}{\mathrm{\Omega }}^f,\hfill \\ J-1\hfill & \mathrm{in}{\mathrm{\Omega }}^s,\hfill \end{array}\right.\hfill \end{array}$$

(13)

the solid to fluid density ratio, denoted by $\beta$, is defined as $\beta =\frac{{\varrho }^s}{{\varrho }^f}$, where ${\varrho }^s$ and ${\varrho }^f$ are the densities of the solid and fluid, respectively. To ensure that the no-slip conditions are maintained at the fluid structure interface, the following Neumann and Dirichlet constraints are applied:

$$\begin{array}{cc}\hfill {\mathsf{œ}}^f\mathbf{n}& ={\mathsf{œ}}^s\mathbf{n}\hfill \\ \hfill {\mathbf{v}}^f& ={\mathbf{v}}^s.\hfill \end{array}$$

(14)

3. The Finite Element Method (FEM) Discretization

The Taylor-Hood pair $P_2{P}_1$ is a widely used finite element space for fluid-structure interaction problems. This choice of element is optimal in terms of convergence, and it also satisfies the LBB stability condition [41]. The Ladyzhenskaya–Babuska–Brezzi (LBB) condition is a stability condition for mixed finite element methods that is necessary to ensure convergence of the numerical solution to the correct solution of the partial differential equation being solved. By transitioning from a a discontinuous to a continuous space, the dimension of the pressure space is reduced, which leads to a reduction in the number of degrees of freedom required for the simulation. The degrees of freedom for this element are illustrated in Figure 1.

The basis function ${\varphi }_i$ of $P_1$ element is given by:

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

(15)

and basis vector set is

$$\langle 1,x,y\rangle .$$

(16)

For $P_2$ element, we have

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

(17)

and basis vector set is

$$\langle 1,x,y,x^2,xy,y^2\rangle$$

(18)

The weak formulation of the problem is obtained using the standard Galerkin finite element method, which involves discretizing the space over the time interval $I=[0,T]$. To achieve this, the equations given in Equations 11–13 are multiplied by the test functions $\mathsf{ı},\mathsf{¸},\gamma$, and then integrated over the space and the time interval I. Integration by parts is then applied to some of the terms, and the boundary conditions are used to obtain a final expression.

$$\begin{array}{cc}\hfill {\int }_0^T{\int }_{\mathrm{\Omega }}\frac{\partial \mathbf{u}}{\partial t}·\mathsf{ı}dVdt=& {\int }_0^T{\int }_{{\mathrm{\Omega }}^s}\mathbf{v}·\mathsf{ı}dVdt-{\int }_0^T{\int }_{{\mathrm{\Omega }}^f}\mathrm{Grad}\mathbf{u}·\mathrm{Grad}\mathsf{ı}dVdt,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\int }_0^T{\int }_{{\mathrm{\Omega }}^f}J\frac{\partial \mathbf{v}}{\partial t}·\mathsf{¸}dVdt& +{\int }_0^T{\int }_{{\mathrm{\Omega }}^s}\beta J\frac{\partial \mathbf{v}}{\partial t}·\mathsf{¸}dVdt\hfill \\ \hfill =& -{\int }_0^T{\int }_{{\mathrm{\Omega }}^f}J\mathrm{Grad}\mathbf{v}{\mathbf{F}}^{-1}(\mathbf{v}-\frac{\partial \mathbf{u}}{\partial t})·\mathsf{¸}dVdt\hfill \\ \hfill & +{\int }_0^T{\int }_{\mathrm{\Omega }}Jp{\mathbf{F}}^{-T}·\mathrm{Grad}\mathsf{¸}dVdt\hfill \\ \hfill & -{\int }_0^T{\int }_{{\mathrm{\Omega }}^f}J\mu \mathrm{Grad}\mathbf{v}{\mathbf{F}}^{-1}{\mathbf{F}}^{-T}·\mathrm{Grad}\mathsf{¸}dVdt,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill 0=& {\int }_0^T{\int }_{{\mathrm{\Omega }}^s}(J-1)\gamma dVdt+{\int }_0^T{\int }_{{\mathrm{\Omega }}^f}\mathrm{Div}\left(J\mathbf{v}{\mathbf{F}}^{-T}\right)\gamma dVdt.\hfill \end{array}$$

(21)

After applying FEM space discretization, we obtained a nonlinear algebraic system of equations as follows:

$$\begin{array}{c}\hfill \left(\begin{array}{ccc}{\mathbf{S}}_{\mathbf{u}\mathbf{u}}& {\mathbf{S}}_{\mathbf{u}\mathbf{v}}& \mathbf{0}\\ {\mathbf{S}}_{\mathbf{v}\mathbf{u}}& {\mathbf{S}}_{\mathbf{v}\mathbf{v}}& k\mathbf{B}\\ c_{\mathbf{u}}{\mathbf{B}}_{\mathbf{s}}^{\mathbf{T}}& c_{\mathbf{v}}{\mathbf{B}}_f^T& \mathbf{0}\end{array}\right)\left(\begin{array}{c}{\mathbf{u}}_h\\ {\mathbf{v}}_h\\ p_h\end{array}\right)=\left(\begin{array}{c}\mathrm{rhs}\mathbf{u}\\ \mathrm{rhs}\mathbf{v}\\ \mathrm{rhs}p\end{array}\right),\end{array}$$

(22)

The nonlinear algebraic system of equations involves the diffusive, reactive, and convective operator denoted by $\mathbf{S}$ from the governing equations, as well as the discrete gradient and the discrete divergence operators represented by $\mathbf{B}$ and ${\mathbf{B}}^{\mathbf{T}}$, respectively.

4. Solution Algorithm

In order to obtain a solution for the system of nonlinear algebraic equations of saddle point type presented in (22), the Newton iteration method is employed. This method is known for its quadratic convergence property when the solution is sufficiently close. The ultimate goal is to find a root of the residual.

$$\mathbf{R}\left(\mathbf{X}\right)=\mathbf{0},$$

(23)

The Newton iteration method with damping is utilized, which involves using the function value and its first derivative that are already known. The corresponding formula for this iteration is as follows:

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

(24)

Here, $\mathbf{X}=({\mathbf{u}}_h,{\mathbf{v}}_h,p_h)$, and the Jacobian matrix $\frac{\partial \mathbf{R}\left({\mathbf{X}}^{\mathbf{n}}\right)}{\partial \mathbf{X}}$ is calculated through finite differences, using the residual matrix $\mathbf{R}\left(\mathbf{X}\right)$.

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

(25)

The coefficients ${\alpha }_j>0$ correspond to increments at each iteration step n, and ${\mathbf{e}}_j$ represents the unit basis vectors in ${\mathbf{R}}^n$ for Equation 24. For more detailed explanations, refer to [42]. The parameter ${\omega }^n\in (-1,0)$ is chosen such that the error measure decreases.

$$\mathbf{R}\left({\mathbf{X}}^{n+1}\right)·{\mathbf{X}}^{n+1}\le \mathbf{R}\left({\mathbf{X}}^n\right)·{\mathbf{X}}^n,$$

The Newton iteration algorithm in Algorithm 1 is used to solve the nonlinear algebraic system of equations of saddle point type in (22). To improve the robustness of the Newton iteration, damping is used when the current approximation ${\mathbf{X}}^n$ is not sufficiently close to the final solution. The damping parameter ${\omega }^n\in (-1,0)$ is chosen to decrease the error measure. Further details on this technique can be found in [38,44]. A direct solver for sparse systems, such as MUMPS [31], is employed in this $2D$ problem.

Algorithm 1 Newton iteration and line search

1: Input the nonlinear tolerance parameter 2: Initialize n to zero and take a starting guess ${\mathbf{X}}^n$ 3: Compute the residual vector $\mathbf{R}\left({\mathbf{X}}^n\right)=\mathbf{A}{\mathbf{X}}^n-\mathbf{b}$ 4: Compute the Jacobian matrix $\mathbf{J}\left({\mathbf{X}}^n\right)=\frac{\partial \mathbf{R}}{\partial \mathbf{X}}\left({\mathbf{X}}^n\right)$ 5: Solve the linear system for the correction of $\delta \mathbf{X}$: $$\mathbf{J}\left({\mathbf{X}}^n\right)\delta {\mathbf{X}}^n=\mathbf{R}\left({\mathbf{X}}^n\right)$$ 6: Find an optimal step length ${\omega }^n\in (-1,0]$ using line search 7: Update the solution as ${\mathbf{X}}^{n+1}={\mathbf{X}}^n+{\omega }^n\delta {\mathbf{X}}^n$

5. Mathematical Modeling

The focus of the study is on the flow of biomagnetic fluid through a stenosis in a bifurcated artery. The fluid has a two-dimensional structure, is viscous, stable, and non-compressible. The researchers assumed that the magnetic field has a negligible impact on fluid viscosity, and that the geometry being studied has insulating walls with a relatively weak electric field. It is generally accepted that the walls of arteries have some degree of linear flexibility. The problem geometry and a coarse mesh are presented in Figure 2.

To analyze the flow of biomagnetic fluid, a constant magnetic field is applied along the $\overline{Y}$ axis. The inflow velocity profile and the outflow pressure are assumed to follow a parabolic shape, with the outflow pressure expected to remain at zero. The governing equations for flow are formulated using the Arbitrary Lagrangian Eulerian (ALE) approach.

In summary, the study focuses on the flow of biomagnetic fluid through a stenosis in a bifurcated artery. The fluid is assumed to be two-dimensional, viscous, stable, and non-compressible. The researchers assume that the magnetic field has a negligible effect on fluid viscosity and that the walls of the artery being studied are insulating with a weak electric field. A constant magnetic field is applied along the $\overline{Y}$ axis, and the inflow velocity profile and outflow pressure are assumed to follow a parabolic shape. The governing equations are derived using the ALE formulation. The problem geometry and mesh are shown in Figure 2.

Conservation of mass

$$\frac{\partial u^{*}}{\partial x^{*}}=-\frac{\partial v^{*}}{\partial y^{*}},$$

(26)

conservation of momentum

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

(27)

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

(28)

The velocity components in different dimensions are represented by $u^{*}$ and $v^{*}$, while $w^{*}$ denotes velocity in the mesh coordinate system. The electrical conductivity of the biomagnetic fluid is denoted by ${\sigma }^{*}$. The Lorentz force per unit volume generated by the fluid’s electric conductivity is represented by ${\sigma }^{*}{u}^{*}{\left(B^{*}\right)}^2$ in equation 27.

The governing equations for dimensional solid displacement can be expressed as follows: conservation of Momentum:

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

(29)

Now that we have reached this point, we are able to define the strain tensor by assuming that the wall is linearly elastic:

$$\sigma =(1/J)FS{F}^T,$$

(30)

The value J of the determinant is employed in the equation $F=1+\nabla u_S$, where F is the value of the determinant. The equation $S=C:ϵ$ describes the relationship between the second Piola-Kirkhoff stress tensor, S, and the strain, $ϵ$ and

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

(31)

To simplify the numerical solution of the system, the first step is to convert it into a dimensionless form by introducing the following non-dimensional variables:

$$\frac{x^{*}}{h}=x,\frac{y^{*}}{h}=y,\frac{u^{*}}{u_0^{*}}=u,\frac{v^{*}}{u_0^{*}}=v,\frac{w^{*}}{u_0^{*}}=w,\frac{p^{*}}{{\rho }_f{u}_0^{*2}}=p,$$

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 and the non dimensional parameters:

$$Re=\frac{h{\rho }_f^{*}{u}_0^{*}}{{\mu }^{*}}\mathrm{Reynolds}\mathrm{number},$$

$$Ha=\frac{{\sigma }^{*}{h}^2{\left(B^{*}\right)}^2}{{\mu }^{*}}\mathrm{Hartmann}\mathrm{number}.$$

The system is transformed into a dimensionless form for simplicity by introducing the aforementioned non-dimensional parameters and variables, resulting in the following set of equations:

$$\nabla ·u=0,$$

(32)

$$-\frac{\partial p}{\partial x}-\frac{H{a}^2}{Re}u+\frac{1}{Re}\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)=(u-w)\frac{\partial u}{\partial x}+(v-w)\frac{\partial u}{\partial y},$$

(33)

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

(34)

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

(35)

5.1. Problem Configuration

The mesh serves as a common 2-dimensional representation of the model, which can be observed in Figure 2. The computational domain can be viewed as a symmetrical bifurcation and stenosis. The walls are constructed using an isotropic linear elastic material with a known Poisson ratio and Young’s modulus. The group of Lame coefficient, $\lambda$, and shear modulus, $\mu$, are used to define the walls and have the following responsibilities:

$$\begin{array}{cc}\hfill \nu & =\frac{\lambda }{2(\lambda +\mu )},E=\frac{\mu (3\lambda +2\mu )}{\lambda +\mu },\hfill \\ \hfill \mu & =\frac{E}{2(1+\nu )},\lambda =\frac{\nu E}{(1+\mu )(1-2\nu )},\hfill \end{array}$$

(36)

Materials with a Poisson ratio $\nu <0.5$ are considered compressible, while materials with $\nu =0.5$ are considered incompressible, as explained in [39]. In the present analysis, a Young’s modulus of $E=5\times {10}^5$ and a Poisson ratio of $\nu =0.49$ were utilized.

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

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

Figure 2 displays a computational domain with stenosis. The parent artery diameter is $h=1$, and the diameter of the stenosis site is reduced by $50\%$ while the daughter artery diameter is $h_1=0.37$. The elastic wall width of this artery is $w=0.08$, and the bifurcation artery angle is ${37}^{\circ }$. Line A starts at the coordinates $(1.8,t)$ for $t=0$ and $t=1$. The centerline for pressure measurements is line C, and the distance between line A and the stenosis is $d=0.8$.

The coarse mesh, shown in Figure 2, has 1478 domain elements and 338 boundary elements. Table 1 shows the absolute error of the WSS with respect to the mesh refinement level and the number of elements. The WSS on the upper elastic wall is also provided in the table. In Level 1 computation, 1922 finite elements are used.

Assuming a parabolic shape for the inlet velocity profile, the following equation is used:

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

(37)

The pressure conditions at the outflow boundaries are determined by the pressure being zero at the discharge point.

6. Results and Discussion

To obtain a solution for the formulated system of equations, the ALE approach is utilized along with appropriate boundary conditions, which helps to gain a better understanding of the problem. The governing nonlinear PDEs are transformed into a dimensionless form using well-known parameters like the Reynolds number and Hartmann number ($Ha$). The Reynolds number ranges from 500 to 2000, and the Hartmann number ranges from 0 to 12 for the fluid flow simulation within the domain. At low Reynolds numbers, such as $Re$ 300 for the magnetic parameter values examined by Ha, significant wall deformation occurs due to the dominance of viscous forces. The study utilizes various parameters like displacement, velocity, pressure, and Wall Shear Stress (WSS) to provide more clarity about the results. Moreover, plotting the velocity magnitude profile at points A and B before and after the stenosis is performed to analyze the primary impact of the stenosis on the fluid flow velocity.

The dimensionless velocity distribution of the blood flow in a symmetric stenosed bifurcated artery subjected to an asymmetric magnetic field with linearly elastic walls is analyzed in Figure 3. The effect of the magnetic field on the flow behavior is assessed. In the absence of a magnetic field, the flow reaches its maximum speed but separates almost immediately after the stenosis, indicating a decrease in wall shear stress. This separation can lead to the development of atherosclerosis. However, the elastic walls show minimal deformation. Upon introducing an external magnetic field, the separation area decreases, resulting in a decrease in velocity and a less noticeable hollow. With an increase in magnetic field strength, the separation zone reduces, leading to the management of thrombosis. In the presence of MHD flow, significant wall deformation may occur, leading to the creation of Wall Shear Stress (WSS).

The velocity profile along the artery at a Reynolds number of 1000 is depicted in Figure 4. When there were no hydrodynamic forces $(Ha=0)$, a stable cavity was formed. However, when a magnetic field was applied, the cavity contracted, resulting in a smaller size. As the magnetic parameter Ha increased from 8 to 10 and then to 12, the cavity disappeared altogether. Increasing the value of Ha led to a considerable amount of wall deformation, which reduced the amount of flow separation.

In Figure 5, one can observe the velocity patterns as they relate to a Reynolds number of 1500. A large recirculation can be seen soon after the stenosis forms at position 11 in the purely hydrodynamic scenario where $Ha=0$. With Reynolds at 1500, as the value of Ha is increased, the size of the cavity gradually reduces. This implies that the amount of recirculation is directly proportional to the Reynolds number of 1500, while it is inversely related to the Hartmann number of n.

Figure 6 demonstrates the change in arterial velocity as Re is set to 2000. In a fully hydrodynamic scenario where Ha equals zero, a large and extended cavity is created shortly after the stenosis. However, the cavity shrinks when a magnetic field is introduced. Even with increasing values of the magnetic parameter, such as Ha equal to 8, 10, and 12, the cavity persists due to the high value of Re equaling 2000 in comparison to Re values of 500 and 1000. As the magnetic parameter Ha approaches values of 8, 10, and 12, the cavity area decreases and wall deformation is minimized. This suggests that the size of the cavity area is directly proportional to the Reynolds number, which is 2000 in this case, whereas the Hartmann number is inversely proportional to the wall deformation and cavity size.

Figure 3. Velocity magnitude for $Re=500$ and $Ha=0,8,10,12$

Figure 3. Velocity magnitude for $Re=500$ and $Ha=0,8,10,12$

Figure 4. Velocity magnitude for $Re=1000$ and $Ha=0,8,10,12$

Figure 4. Velocity magnitude for $Re=1000$ and $Ha=0,8,10,12$

Figure 5. Velocity magnitude for $Re=1500$ and $Ha=0,8,10,12$

Figure 5. Velocity magnitude for $Re=1500$ and $Ha=0,8,10,12$

Figure 6. Velocity magnitude for $Re=2000$ and $Ha=0,8,10,12$

Figure 6. Velocity magnitude for $Re=2000$ and $Ha=0,8,10,12$

6.1. Displacement

The results shown in Figure 7 provide insight into the total displacement of the upper wall (and lower wall, due to the symmetric nature of the geometry) of a bifurcated artery when subjected to different Hartmann numbers, specifically $Ha=0,8,10,12$. It can be observed that as the Reynolds number, $Re$, increases, there is a decreasing trend in the total displacement. The maximum displacement is seen just before the stenosis, and it begins to decrease after the stenosis. This reduction in the total displacement is a result of the contraction of the recirculation cavity caused by the application of the magnetic field. As the Hartmann number is increased, the wall deformation is reduced, leading to a decrease in the size of the recirculation cavity and a subsequent reduction in the total displacement.

6.2. The Velocity Profile at Position A and B

This section focuses on the velocity profile behavior at positions A and B. Comparing the relative velocities of points A and B shows a significant discrepancy. When the Reynolds number is constant, the velocity profile is inversely proportional to the magnetic field parameter, while the magnetic field parameter is directly proportional to the velocity profile. Table 2 summarizes the maximum velocity at location A before stenosis and location B after stenosis, along with the percentage increase, for various Reynolds numbers. The results demonstrate that the maximum velocity at location A before stenosis increases with increasing Reynolds number, while the maximum velocity at location B after stenosis decreases with increasing Reynolds number. This trend indicates that high Reynolds numbers tend to induce a more uniform velocity profile, whereas low Reynolds numbers tend to create a larger discrepancy between the velocity profiles at positions A and B.

Figure 8. The velocity magnitude profile at location A (left) and B (right)

Figure 8. The velocity magnitude profile at location A (left) and B (right)

6.3. Pressure at the centre

The pressure behavior in the bifurcated artery is an important aspect of studying its fluid dynamics. In Figure 9, the pressure attitude is analyzed against the magnetic field parameter, $Ha=0,8,10,12$, and the values are plotted against different Reynolds numbers, i.e., $Re=500,1000,1500,2000.$

The pure hydrodynamic scenario with $Ha=0$ yields the minimum pressure in the artery. As the value of $Ha$ increases, the pressure in the bifurcation region increases as well. This indicates a direct relation between the magnetic field parameter and the pressure.

In addition, the pressure behavior is also related to the Reynolds number. As the Reynolds number increases, the pressure decreases for a fixed value of the magnetic field parameter. This shows an inverse relation between the Reynolds number and the pressure in the bifurcation region.

Moreover, it is observed that at a lower Reynolds number, i.e., $Re=500$, the maximum pressure occurs for a higher value of the magnetic field parameter, i.e., $Ha=8,10,12$. This result is particularly interesting as it suggests that the behavior of the pressure in the bifurcation region is affected by both the Reynolds number and the magnetic field parameter.

Figure 9. Pressure values at the centre C of parent artery

Figure 9. Pressure values at the centre C of parent artery

6.4. Wall Shear Stress

In this study, we have investigated the effects of the magnetic field and Reynolds number on the wall shear stress in a stenosed bifurcated artery. Our results, presented in Table 3 and Figure 10, show that the wall shear stress is directly influenced by the magnetic parameter $Ha$ and has an inverse relation with the Reynolds number $Re$.

Interestingly, we found that for small values of $Re$, an increasing behavior can be observed for higher values of $Ha$. Specifically, our results suggest that an artery with a small Reynolds number flow, i.e., $Re=500$, may experience significant stress in the presence of a higher value of the magnetic field parameter, i.e., $Ha=12$.

It is important to note that 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 [43]. The plaque formation occurs due to the 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.

Additionally, we observed that the values of wall shear stress remain the same for the upper wall and lower wall due to the symmetric nature of the considered geometry in Figure 10. Therefore, our results for wall shear stress are valid for both the upper and lower walls of the bifurcated artery. Overall, our study provides important insights into the effects of magnetic field and Reynolds number on wall shear stress and their potential implications for the development of atherosclerosis in stenosed bifurcated arteries.

Figure 10. The relationship between the Hartmann number and the total wall shear stress (WSS) along the upper wall

Figure 10. The relationship between the Hartmann number and the total wall shear stress (WSS) along the upper wall

7. Conclusions

This computational study investigated the flow of a biomagnetic fluid through a two-dimensional, stable, incompressible, and laminar arterial stenosis with elastic wall surfaces. To solve the governing nonlinear partial differential equations of fluid and solid (elasticity), an arbitrary Lagrangian Eulerian method was used in combination with a bi-quadratic linear stable finite element pair. The nonlinear algebraic problem was solved using the Newton-Raphson method. Charts and tables were included and analyzed to investigate the key features of the obtained solutions. The following are some concluding remarks based on the numerical analysis findings:

• The presence of a magnetic field parameter, $Ha$, in increasing order reduces the size of the cavity adjacent to the stenosis.
• A higher value of the magnetic field parameter decreases the velocity gain after the stenosis compared to the purely hydrodynamic case ($Ha=0$). For example, when $Re=500$, the velocity gain is 41
• The pressure increases with higher values of Hartmann number ($Ha=0,8,10,12$).
• A flow with a small Reynolds number ($Re=500$) and a higher magnetic field parameter value ($Ha=12$) results in a significant increase in wall shear stress.
• The wall displacement is larger when $Re=500$ and $Ha=12$, and the maximum wall displacement is observed just before the stenosis.

In conclusion, this numerical analysis suggests that the magnetic field has significant effects on blood flow behavior. This suggests that magnetic fields could potentially be implemented in a wide range of treatments for conditions such as plaque, atheroma, hypertension, blood pressure, and other similar conditions. However, to further validate these findings, more research based on patient-specific data is needed. Hence, further investigations are always necessary to ratify the conclusions drawn in this study.