Variational multiscale stabilized FEM for cardiovascular flows in complex arterial vessels under magnetic forces

Dipak Kumar Sahoo; Anil Rathi; B. V. Rathish Kumar

doi:10.1515/cmb-2023-0118

Article Open Access

Variational multiscale stabilized FEM for cardiovascular flows in complex arterial vessels under magnetic forces

Dipak Kumar Sahoo , Anil Rathi and B. V. Rathish Kumar

Published/Copyright: March 20, 2024

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Computational and Mathematical Biophysics Volume 12 Issue 1

Abstract

In this study, we present a variational multiscale stabilized finite element method for steady-state incompressible fluid flow under magnetic forces. In particular, an algebraic approach of approximating the subscales has been considered, and then, the stabilization parameters are derived using Fourier analysis. The proposed scheme is used to trace the blood flow dynamics in complex arterial vessels under multiple pathological conditions. We examine the pressure and stress distribution in addition to the flow pattern to assess the criticality of the diseased condition.

Keywords: stabilized finite element method; magnetohydrodynamic Navier-Stokes; diseased artery; wall shear stress

MSC 2010: 65M60; 76M10; 92B05; 74S05; 74S10; 76M12; 46N60; 65L60

1 Introduction

The intricate arterial system in the human body uses the heart to pump blood, which is then circulated throughout the entire body. Problems may occasionally arise as a result of high blood pressure, low blood pressure, an increase in cholesterol, etc. Hence, blood flow is impacted by an artery obstruction, which has an impact on the person’s entire body because of the complexity of the arterial network as seen in [18].

Figure 1

Vessels with aneurysm.

In the beginning, we recall a few arterial pathological conditions that are commonly noted in clinical contexts. Aneurysms manifest when the wall of an artery bulges due to the weakening of the vessel wall (Figure 1) [1]. The most serious threat of an aneurysm is that it can burst and be life-threatening. A large aneurysm can affect blood circulation and lead to downstream blood clot manifestation. And sometimes portions of an artery shrink to a thinner shape, affecting the flow, which we call stenosis (Figure 2) [17]. And stenosed cardiovascular vessels lead to health complications. A lot of work has been carried out to understand the altered flow conditions in diseased arteries. Ali et al. [2] have studied the unsteady pulsatile flow of blood through a tapered stenotic artery, where the rheology of blood is described using the Sisko model. The governing equations and the required boundary conditions are integrated across the whole arterial segment using the finite difference approach. In [20], viscous flow in stenotic elastic tubes with significant wall deformation and collapse has been investigated using generalized finite differences (GFDs) and a three-dimensional computer model with fluid-wall interactions [3,11,12]. The demonstration of an in vitro experimental setup modeling blood flow and arterial collapse in carotid arteries with stenosis has been carried out. In [21], under the impact of periodic body acceleration, the pulsatile blood flow via a porous-saturated stenotic artery has been investigated. The constitutive equation of the Cross model and a modified version of Darcy’s law that applies to the Cross model are both used to characterize the blood to study the flow over a W-shaped stenosis. The pathogenic condition in which blood flows via an artery blocked by blood clots and fatty cholesterol is modeled, and finite difference method (FDM) is used to obtain numerical solutions. The pulsatile flow of Casson nanofluid through a tapered, stenosed artery that is inclined in the presence of a magnetic field is explored by Ponalagusamy and Priyadharshini in [19], while iron oxide nanoparticles are permitted to flow alongside it. Using FDM, the momentum equations for temperature and concentration are coupled and solved to obtain the solutions for velocity, temperature, concentration, wall shear stress (WSS), and resistance to blood flow. Variations with respect to the Grashof number, Reynolds number, and Prandtl number are also investigated. A non-invasive, high-resolution, and quick method for the prediction of aneurysm-prone artery segments has been developed by Dolgov et al. [13] with a computer-based methodology using micro-computed tomography scanning and finite element analysis. In [16], Kumar et al. have investigated the influence of abdominal aortic aneurysms’ shape on local and global hemodynamic parameters and rupture predictions using finite element software packages. In [6], Cherkaoui et al. have studied the effect of magnetic field on blood flow, as it is a very significant topic for bio-medical engineers because of the wide use of ferrofluids in industry and medicine. In the field of magnetic separation tools, anti-cancer drug carriers, micro-valve applications, and chemotherapy ferrofluids are used. Besides, the magnetic field affects the blood flow characteristics through the arteries. The magnetic field affecting blood flow through arteries has been discussed.

Figure 2

Vessels with stenosis.

Figure 3

Horizontal and vertical velocity when compared to Ghia et al. [17]: (a) horizontal velocity comparison and (b) vertical velocity comparison.

However, in the context of arterial vessels, multiple pathologies under external forces like those from magnetic fields have not been reported. More importantly, to analyze blood flow in vessels with multiple complexities, there is a need for a robust numerical method. We are going to address both of these needs. This study becomes vital since people are constantly exposed to a variety of magnetic fields, either mildly or severely. Academicians, scientists, and engineers exposed to magnetic fields include those working in labs, those undergoing computed tomography scans, magnetic resonance imaging, or X-rays in hospitals, and those operating heavy machinery. By studying these multiply-diseased arteries under magnetic forces, we are able to understand the influence of magnetic forces on the hemodynamics under such diseased conditions and thereby find suitable remedial measures. Also, we investigate the changes in the velocity field, re-circulation zones, pressure field, and WSS in a multiply-diseased arterial vessel with the increase in Hartman number, which is going to help during drug delivery. Furthermore, we have developed a robust numerical scheme that works well for flow simulations with high Reynolds numbers and vessels with complex geometries. It is well known that in convection-dominated flows, for which layers appear where velocity and pressure exhibit a rapid variation, the Galerkin approach leads to numerical oscillations in these layer regions, which pollute the entire solution domain. And the use of inappropriate combinations of interpolation functions to represent velocity and pressure fields can yield unstable schemes. So the classical Galerkin approach fails. Hence, we use a more stable scheme, namely, the subgrid-scale (SGS) stabilized finite element method. Generally, two approaches to SGS stabilized formulation, namely, the algebraic approach (ASGS), and the orthogonal projection approach (OSGS), have been studied. In [15], Hughes has introduced the concept of stabilized multi-scale subgrid method for the Helmholtz equation. Furthermore, in [8], Codina presents the SGS formulation of the ADR equation with constant coefficients. In [7], Chowdhury and Kumar have extended the SGS formulation for the advection diffusion reaction (ADR) equation with variable coefficients and performed a priori and a posteriori error analysis. In [9], Codina has studied a stabilized finite element method for solving systems of convection-diffusion-reaction equations again using the algebraic SGS formulation. However, the authors of [3–5,11,12] have worked with the OSGS method for elliptic partial differential equations. Here, we have derived a multi-scale stabilized finite element method that leads to highly stable solutions. In particular, the algebraic subgrid-scale method, i.e., the subgrid scale method with an algebraic approximation to the subscales, is considered, and the stabilization parameters are derived using a Fourier analysis, which is another key aspect of this study.

1.1 Organization of sections

In this article, in Section 2, governing equations are presented that model the blood flow in the arteries under the effect of magnetic forces. In Section 3, we derive a multi-scale variational stabilized finite element scheme for the governing equations, where the stabilization parameters are derived using Fourier analysis. This method is stable, and it gives more realistic data. In Section 4, we show numerical simulations and results where flows over arteries with multiple pathological conditions are studied. First, a flow simulation is carried out for a vessel with multiple aneurysms. Then, we discuss the flow when there are multiple stenoses in the artery. Then, flow over a combination of aneurysms and stenosis is demonstrated. In the end, the flow over a child’s aortic artery with multiple branches is presented.

2 Governing equation

We consider the steady flow in vessels under the effect of a magnetic field. The flow is assumed to be incompressible and Newtonian. Magnetic forces in the horizontal direction are modeled by the Navier-Stokes equation, with magnetic terms appearing in the Y-momentum equation.

Let Ω be an open, bounded, polygonal domain in R 2 with piece-wise smooth boundary ∂ Ω . The following three non-dimensional equations govern the flow.

X-momentum equation:

(1) u 1 ∂ u 1 ∂ x + u 2 ∂ u 1 ∂ y − 1 Re Δ u 1 + ∂ p ∂ x = f 1 , in Ω .

Y-momentum equation:

(2) u 1 ∂ u 2 ∂ x + u 2 ∂ u 2 ∂ y − 1 Re Δ u 2 + ∂ p ∂ y + Ha 2 Re u 2 = f 2 , in Ω .

Continuity equation:

(3) ∇ ⋅ u = 0 , in Ω .

Boundary condition:

(4) u = 0 , on ∂ Ω ,

where, u = [ u 1 , u 2 ] represents the velocity of blood particles inside the vessels, p denotes the pressure inside the vessels, f 1 and f 2 represent the source terms that can represent external body forces, Re denotes the Reynolds number, and Ha represents the effect of the magnetic field. Our goal is to find u and p when the speed of flow and intensity of the magnetic field vary.

3 Weak formulation and stabilization scheme derivation

3.1 Weak formulation

We wish to find u ∈ V = [ H 0 1 ( Ω ) ] 2 and p ∈ Q = L 2 0 ( Ω ) ∕ R such that the following equations are satisfied for all v = [ v 1 , v 2 ] ∈ [ H 0 1 ( Ω ) ] 2 and q ∈ L 2 0 ( Ω ) ∕ R such that

(5) u 1 ∂ u 1 ∂ x , v 1 + u 2 ∂ u 1 ∂ y , v 1 − 1 Re ( ∇ u 1 , ∇ v 1 ) − p , ∂ v 1 ∂ x = ( f 1 , v 1 ) ∀ v 1 ∈ H 0 1 ( Ω ) ,

(6) u 1 ∂ u 1 ∂ x , v 2 + u 2 ∂ u 1 ∂ y , v 2 − 1 Re ( ∇ u 2 , ∇ v 2 ) − p , ∂ v 2 ∂ y + Ha 2 Re ( u 2 , v 2 ) = ( f 2 , v 2 ) ∀ v 2 ∈ H 0 1 ( Ω ) ,

(7) ( ∇ ⋅ u , q ) = 0 ∀ q ∈ L 2 0 ( Ω ) ∕ R .

3.2 Discrete formulation

We use finite element space discretization to create the discrete variational formulation. The domain Ω is discretized using non-overlapping triangles K i . And the parameter h denotes: h = max K i ∈ T h ( h K i ) , where h K i is the diameter of subdomain K i of Ω k . Now, we construct the appropriate finite dimensional subspaces V h ⊂ V and Q h ⊂ Q . The standard Galerkin finite element formulation is given as: find u h = ( u h 1 , u h 2 ) ∈ V h and p h ∈ Q such that

(8) u 1 h ∂ u 1 h ∂ x , v 1 h + u 2 h ∂ u 1 h ∂ y , v 1 h − 1 Re ( ∇ u 1 h , ∇ v 1 h ) − p h , ∂ v 1 h ∂ x = ( f 1 , v 1 h ) ∀ v 1 h ∈ H 0 1 ( Ω ) ,

(9) u 1 h ∂ u 1 h ∂ x , v 2 h + u 2 h ∂ u 1 h ∂ y , v 2 h − 1 Re ( ∇ u 2 h , ∇ v 2 h ) − p , ∂ v 2 h ∂ y + Ha 2 Re ( u 2 h , v 2 h ) = ( f 2 , v 2 h ) ∀ v 2 h ∈ H 0 1 ( Ω ) ,

(10) ( ∇ ⋅ u h , q h ) = 0 ∀ q h ∈ L 2 0 ( Ω ) ∕ R .

3.3 Stabilization scheme

The Galerkin formulation alone might not produce a stable result since the flow would be modeled at a high Reynolds number. Here, the velocity and pressure components are divided into coarse and fine scales using the algebraic variational subgrid-scale approach. The impacts of the fine scale are then stated in terms of the coarse scale, giving us the whole equation, which proves to be a fairly stable system. In this case, Fourier analysis is used to obtain the stabilization parameters. First, we decompose the space V = V h ⊕ V ˜ , and similarly, Q = Q h ⊕ Q ˜ , where V h and Q h are called the course scale velocity and pressure components, respectively, and V ˜ and Q ˜ are called the fine-scale velocity and pressure components, respectively. Here, V h is the finite-dimensional subspace of V . Now, we decompose the elements of the space as v = v h ⊕ v ˜ and q = q h ⊕ q ˜ .

We will demonstrate the subgrid-scale stabilized variational finite element formulation for the convection-diffusion reaction equation:

(11) − k Δ u + a ⋅ ∇ u + s u = f , in Ω ,

(12) u = 0 , on ∂ Ω ,

where k > 0 is the diffusion coefficient, s ≥ 0 is the reaction coefficient, a ∈ R 2 is the advection coefficient, and f is a given function.

3.3.1 Variational formulation

Find u ∈ V = H 0 1 ( Ω ) such that B ( u , v ) = ( f , v ) , where

(13) B ( u , v ) = k ( ∇ u , ∇ v ) + ( a ⋅ ∇ u , v ) + s ( u , v ) .

Now, we split the unknown u into components u h and u ˜ with the former being called course-scale and latter fine-scale component. Then, equation (11) can be decomposed into two variational subequations, We seek u h ∈ V h and u ˜ ∈ V ˜ such that

(14) B ( u h , v h ) + B ( u ˜ , v h ) = ( f , v h ) , ∀ v h ∈ V h ,

(15) B ( u h , v ˜ ) + B ( u ˜ , v ˜ ) = ( f , v ˜ ) , ∀ v ˜ ∈ V ˜ .

These equations can further be written as:

(16) B ( u h , v h ) + ⟨ u ˜ , ℒ * v h ⟩ = ⟨ f , v h ⟩ , ∀ v h ∈ V h ,

(17) ⟨ ℒ u h , v ˜ ⟩ + ⟨ ℒ u ˜ , v ˜ ⟩ = ⟨ f , v ˜ ⟩ , ∀ v ˜ ∈ V ˜ .

From equation (17), we can deduce that ℒ u ˜ = ℛ , i.e.,

(18) − k Δ u ˜ + a · ∇ u ˜ + s u ˜ = ℛ ,

where ℛ = f − ℒ u h .

Then, using the Fourier analysis, we approximate u ˜ by u ˜ ≈ τ ℛ , where τ is a parameter to be determined.

Fourier transform of a function g is defined on each triangle K ⊂ of Ω by:

(19) g ˆ ≔ ∫ K e − i k . x h g ( x ) d Ω x ,

where i = − 1 , h is diameter of element K , and k = ( k 1 , k 2 ) is the dimensionless wave number.

For functions with high wave numbers, the first-order derivative can be approximated by:

(20) ∂ ˆ g ∂ x j ( k ) ≈ i k j h ( k ) .

Now, taking the Fourier transform of (13), we will have

(21) u ˜ ˆ ≈ T ( k ) R ˆ ( k ) , T ( k ) ≔ k ∣ k 2 ∣ h 2 + i a . k h + s − 1 .

Using Plancherel’s formula and mean value theorem there exists a wave number k 0 for which we can deduce as in [10]:

(22) ‖ u ˜ ‖ ≈ ∣ T ( k 0 ) ∣ ‖ R ‖ k .

Now, identifying T with ∣ T ( k 0 ) ∣ , we can find [10]:

(23) T = c 1 k h 2 + s 2 + c 2 ∣ a ∣ h 2 − 1 ∕ 2 ,

where c 1 and c 2 are constants.

Now, substituting the value of u ˜ in equation (16), we will have our stabilization scheme, which is wholly in terms of course scale and is computable:

(24) B ( u h , v h ) + ⟨ T ℛ , ℒ * v h ⟩ = ⟨ f , v h ⟩ , ∀ v h ∈ V h .

3.4 Stabilization scheme for our equation

For our problem, a similar thing can be followed after linearizing the equations:

(25) B ( U , V ) = u 1 ∂ u 1 ∂ x , v 1 + u 2 ∂ u 1 ∂ y , v 1 − 1 Re ( ∇ u 1 , ∇ v 1 ) − p , ∂ v 1 ∂ x + u 1 ∂ u 1 ∂ x , v 2 + u 2 ∂ u 1 ∂ y , v 2 − 1 Re ( ∇ u 2 , ∇ v 2 ) − p , ∂ v 2 ∂ y + Ha 2 Re ( u 2 , v 2 ) + ( ∇ ⋅ u , q ) ,

where U = [ u 1 , u 2 , p ] and V h = [ v 1 , v 2 , q ] and

(26) F ( V ) = ( f 1 , v 1 ) + ( f 2 , v 2 ) .

The residues are found to be as follows:

(27) ℛ h = f 1 + 1 Re Δ u 1 h − u h ⋅ ∇ u 1 h − ∂ p h ∂ x f 2 + 1 Re Δ u 2 h − u h ⋅ ∇ u 2 h − ∂ p h ∂ y − Ha 2 Re u 2 h − ∇ ⋅ u h .

And the adjoint ℒ * can be found to be as follows:

(28) ℒ * = − 1 Re Δ v 1 h − u h ⋅ ∇ v 1 h − ∂ q h ∂ x − 1 Re Δ v 2 h − u h ⋅ ∇ v 2 h − ∂ q h ∂ y − Ha 2 Re v 2 h − ∇ ⋅ v h .

Using Fourier analysis as in the above, the stabilization parameters can be found to be

(29) τ 1 = c 1 Re ⋅ h 2 2 + c 2 ∣ u h ∣ 2 h 2 − 1 ∕ 2 ,

(30) τ 2 = c 3 Re ⋅ h 2 + Ha 2 Re 2 + c 4 ∣ u h ∣ 2 h 2 − 1 ∕ 2 ,

(31) τ 3 = 1 Re .

Now, the final stabilization scheme can be given as:

(32) B ( U h , V h ) + ⟨ T ℛ , ℒ * V h ⟩ = ⟨ F , V h ⟩ , ∀ V h ∈ V h ,

where U h = [ u 1 h , u 2 h , p h ] and V h = [ v 1 h , v 2 h , q h ] ,

which in its expanded form is as follows:

u 1 h ∂ u 1 h ∂ x , v 1 h + u 2 h ∂ u 1 h ∂ y , v 1 h − 1 Re ( ∇ u 1 h , ∇ v 1 h ) − p , ∂ v 1 h ∂ x − p h , ∂ v 2 h ∂ y + u 1 h ∂ u 1 h ∂ x , v 2 h + u 2 h ∂ u 1 h ∂ y , v 2 h − 1 Re ( ∇ u 2 h , ∇ v 2 h ) + Ha 2 Re ( u 2 h , v 2 h ) + ∑ K T 1 ( − 1 Re Δ v 1 h − u h ⋅ ∇ v 1 h − ∂ q h ∂ x ) , f 1 + 1 Re Δ u 1 h − u h ⋅ ∇ u 1 h − ∂ p h ∂ x + ∑ K T 2 − 1 Re Δ v 2 h − u h ⋅ ∇ v 2 h − ∂ q h ∂ y , f 2 + 1 Re Δ u 2 h − u h ⋅ ∇ u 2 h − ∂ p h ∂ y − Ha 2 Re u 2 h + ∑ K ( T 3 ( − ∇ ⋅ v h , − ∇ ⋅ u h ) + ( ∇ ⋅ u h , q h ) = ( f 1 , v 1 h ) + ( f 2 , v 2 h ) .

4 Simulations and results

The aforementioned new scheme is validated by considering the flow over a lid-driven cavity, which is elaborately described in Section 4.1. Newton’s quasi-linearization approach is used to handle the nonlinear term. For our study, we have considered four different complex models that could occur in clinical contexts. At first, a model with two stenoses is considered. Second, a model with an aneurysm and stenosis is taken, and finally, as a third model, flow over a child’s aorta is studied. Then, in Section 4.2, we have performed grid validation tests for all our models and chosen the optimal grid size for our study. Then, separately for each complex model, we analyzed the flow dynamics, regions of maximum and minimum WSS, and pressure contours. We have considered continuous piecewise quadratic finite element space (P2) for approximating velocity and continuous piecewise linear finite element space (P1) for approximating pressure.

4.1 Code validation

To validate our code, we have considered Navier-Stokes flow in a lid-driven cavity, where the upper wall is moving in the right direction at a constant speed (1,0) and a no-slip boundary condition is applied to the other three walls. In Figure 3, we have plotted the horizontal velocity on the line X = 0.5 and the vertical velocity on the line Y = 0.5 at different Reynolds numbers. It shows the similarity with the data taken from Ghia et al. [14]. Hence, our scheme is validated.

Figure 4

Grid validation: (a) horizontal velocity and (b) vertical velocity.

4.2 Grid validation

To establish the optimal grid size, we first carry out the grid validation test. For model 1, we have tested the code using five different grid sizes, which are characterized by different degrees of freedom (DOF): 7,597, 9,637, 11,873, 14,961, and 17,769, respectively. Then, in Figure 4, the horizontal and vertical velocity plots are displayed on the lines X = 3.2 and Y = 0.4 , respectively. For horizontal velocity, we observe a minimal change in the graphs. For vertical velocity with a DOF of 7,597, there is a variation in the graph; however, for all other grid sizes, the change is minimal. So, the grid with a DOF of 11,873 is best chosen for the whole study.

Figure 5

Velocity streamlines: (a) Re = 50, without magnetic effect (Ha = 0); (b) Re = 50, under mild magnetic influence (Ha = 10); and (c) Re = 50, under high magnetic influence (Ha = 100).

Similarly, for models 2 and 3, grid validation tests were carried out. For our study in this article, grid sizes with a DOF of 13,069 and 14,333 are chosen, respectively, for model vessel 2 and model vessel 3.

4.3 Model vessel 1

First, a model vessel with two stenoses with varying degrees of stenosis is considered. A steady flow with Re = 50 is simulated. Such a study represents the flow in arteries with a radius of 1 mm. Flow pattern and WSS results are studied under varying magnetic force effects.

From Figure 5, we can observe that the flow velocity is maximum near the stenosed regions, which was expected because of the narrowing of the artery. Also, we observe that with the increase in magnetic field intensity, the flow has slowed down and the area of the recirculation zone has increased. This would be helpful for drug delivery while treating the stenosed region. From Figure 6, it is evident that with the increase in magnetic field intensity, the pressure has come down. Consequently, the intensity of the circulation will decrease, and thereby, the danger of blood element damage will be reduced. From WSS plots in Figure 7, we observe that near the stenosed regions, WSS is maximum. And the WSS near the first stenosed region is greater than that near the second one due to flow and geometric variation. It may be noted that while the first constriction is one-sided, the second one is due to the lumen reduction from both walls. Also, the maximum WSS observed near the first stenosed region is found to significantly raise as Ha increases from 10 to 100 due to flow intensity variation resulting from raising magnetic forces.

Figure 6

Pressure contours: (a) Re = 50, without magnetic effect (Ha = 0); (b) Re = 50, under mild magnetic influence (Ha = 10); and (c) Re = 50, under high magnetic influence (Ha = 100).

Figure 7

WSS plots: Re = 50.

Figure 8

Velocity plots: (a) Re = 100, without magnetic effect (Ha = 0); (b) Re = 100, under small magnetic effect (Ha = 10); and (c) Re = 100, under high magnetic effect (Ha = 100).

4.4 Model vessel 2

Now, we consider the model with one stenosis and an aneurysm, where a steady flow with Re = 50,100 is simulated. This would normally represent the flow in typical arteries with a radius of 2–4 mm. Flow velocity and WSS results are studied under the influence of magnetic fields.

From Figure 8, we observe that the flow velocity is minimum near the aneurysm region and maximum near the stenotic region. This was expected because of the areal increase in the aneurysm region and the narrowing near the stenotic region. From Figure 9, it is evident that with an increase in magnetic field intensity, the pressure has come down, which is evident from the intensity of the green color before the stenotic region. From the plot in Figure 10, we observe that near the head and toe of the aneurysm, WSS is relatively large, and near the stenotic region, WSS is observed to be maximum. Also, with an increase in magnetic field, the shear stress at both ends of the aneurysm and in the stenotic region is increasing. So there is an increased chance for the wall to rupture near the head and the toe region of the aneurysm.

Figure 9

Pressure plots: (a) Re = 100, without magnetic effect (Ha = 0); (b) Re = 100, under small magnetic effect (Ha = 10); and (c) Re = 100, under high magnetic effect (Ha = 100).

Figure 10

Vessels with aneurysm: (a) Re = 50, under the influence of magnetic effect and (b) Re = 100, under the influence of magnetic effect.

Figure 11

Velocity plots: (a) Re = 3,000 (Ha = 0), (b) Re = 3,000 (Ha = 10), and (c) Re = 3,000 (Ha = 100).

4.5 Model vessel 3

Finally, we have considered flow across the aorta in our third model (Figure 2). Here, the aortic artery, with several branches carrying oxygen-rich blood to different parts of the body, is considered. It travels across the chest and belly after beginning in the bottom-left corner of the heart. Blood vessels leave the aorta along the route and reach out to various organs and supporting tissues. Re = 3,000 is used in the simulation as this is the average value for flow along the aorta.

Figure 11 displays the flow velocity across the aorta. Maximum velocity is seen in the vicinity of the first bifurcation, and it progressively declines as the flow moves past. And when the magnetic effects grow, the velocity magnitudes decrease. From Figure 12, pressure is observed to be relatively large near all the bifurcation regions and decreases gradually. Such a pressure variation is realistic with the systemic flow circulation in the aortic artery with bifurcations.

Figure 12

Pressure plots with Re = 3,000 and under varying magnetic fields: (a) Re = 3,000 (Ha = 0), (b) Re = 3,000 (Ha = 10), and (c) Re = 3,000 (Ha = 100).

These studies clearly indicate the capability of the proposed stabilized finite element method to successfully simulate the flow in multiply-affected complex arterial geometries with bulges, tapering, constrictions, and bifurcations.

5 Conclusion

We have derived a new variational multi-scale subgrid stabilized finite element method for simulating convection-dominated flows in complex geometries. Flow in several complex arterial vessels has been successfully simulated with physiologically meaningful results. The scheme works well with both high and low blood velocities, i.e., for tracing the flow in the aorta, arteries, and arterioles with anatomically realistic geometrical variations. This study also provides input into the use of magnetic forces to facilitate effective remedial measures while treating multiple affected arterial vessels, like in drug delivery. Overall, with the increase in magnetic field intensity, the flow velocity and pressure come down. The hemodynamic parameters such as WSS provide insight into the potential rupture sites, especially in the case of vessels with aneurysms. Subtle features such as flow separation in the vicinity of stenotic regions are successfully captured by the new computational scheme.

Acknowledgment

Authors thank Chitranjan Panndey and Sanehlata, research scholars of the department of mathematics, IIT Kanpur, for their support during the simulations.

Funding information: This research received no specific grant from any funding agency, commercial, or nonprofit sectors.
Conflict of interest: The authors have no conflicts of interest to disclose.
Ethics statement: This research did not require ethical approval.
Data availability statement: The authors confirm that the data supporting the findings of this study are available within the article.

References

[1] Acosta, J. M., Cayron, A. F., Dupuy, N., Pelli, G., Foglia, B., Haemmerli, J., … Morel, S. (2021). Effect of aneurysm and patient characteristics on intracranial aneurysm wall thickness. Frontiers in Cardiovascular Medicine, 8, 775307.10.3389/fcvm.2021.775307Search in Google Scholar PubMed PubMed Central

[2] Ali, N., Zaman, A., & Sajid, M. (2014). Unsteady blood flow through a tapered stenotic artery using Sisko model. Computers & Fluids, 101, 42–49. 10.1016/j.compfluid.2014.05.030Search in Google Scholar

[3] Baiges, J., & Codina, R. (2017). Variational multiscale error estimators for solid mechanics adaptive simulations: An orthogonal subgrid scale approach. Computer Methods in Applied Mechanics and Engineering, 325, 37–55. 10.1016/j.cma.2017.07.008Search in Google Scholar

[4] Bayona, C., Baiges, J., & Codina, R. (2018). Variational multiscale approximation of the one-dimensional forced Burgers equation: The role of orthogonal subgrid scales in turbulence modeling. International Journal for Numerical Methods in Fluids, 86(5), 313–328. 10.1002/fld.4420Search in Google Scholar

[5] Castillo, E., & Codina, R. (2019). Dynamic term-by-term stabilized finite element formulation using orthogonal subgrid-scales for the incompressible Navier-Stokes problem. Computer Methods in Applied Mechanics and Engineering, 349, 701–721. 10.1016/j.cma.2019.02.041Search in Google Scholar

[6] Cherkaoui, I., Bettaibi, S., Barkaoui, A., & Kuznik, F. (2022). Magnetohydrodynamic blood flow study in stenotic coronary artery using lattice Boltzmann method. Computer Methods and Programs in Biomedicine, 221, 106850. 10.1016/j.cmpb.2022.106850Search in Google Scholar PubMed

[7] Chowdhury, M., & Kumar, B. R. (2020). On subgrid multiscale stabilized finite element method for advection-diffusion-reaction equation with variable coefficients. Applied Numerical Mathematics, 150, 576–586. 10.1016/j.apnum.2019.10.021Search in Google Scholar

[8] Codina, R. (1998). Comparison of some finite element methods for solving the diffusion-convection-reaction equation. Computer Methods in Applied Mechanics and Engineering, 156, 185–210. 10.1016/S0045-7825(97)00206-5Search in Google Scholar

[9] Codina, R. (2000). On stabilized finite element methods for linear systems of convection-diffusion-reaction equations. Computer Methods in Applied Mechanics and Engineering, 188(1–3), 61–82. 10.1016/S0045-7825(00)00177-8Search in Google Scholar

[10] Codina, R. (2002). Stabilized finite element approximation of transient incompressible flows using orthogonal subscales. Computer Methods in Applied Mechanics and Engineering, 191(39–40), 4295–4321. 10.1016/S0045-7825(02)00337-7Search in Google Scholar

[11] Codina, R. (2008). Analysis of a stabilized finite element approximation of the Oseen equations using orthogonal subscales. Applied Numerical Mathematics, 58(3), 264–283. 10.1016/j.apnum.2006.11.011Search in Google Scholar

[12] Codina, R., & Blasco, J. (2002). Analysis of a stabilized finite element approximation of the transient convection-diffusion-reaction equation using orthogonal subscales. Computing and Visualization in Science, 4, 167–174. 10.1007/s007910100068Search in Google Scholar

[13] Dolgov, V. Y., Klyshnikov, K. Y., Ovcharenko, E. A., Glushkova, T. V., Batranin, A. V., Agienko, A. S., & Kutikhin, A. G. (2019). Finite element analysis-based approach for prediction of aneurysm-prone arterial segments. Journal of Medical and Biological Engineering, 39, 102–108. 10.1007/s40846-018-0422-xSearch in Google Scholar

[14] Ghia, U. K. N. G., Ghia, K. N., & Shin, C. T. (1982). High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method. Journal of computational physics, 48(3), 387–411. 10.1016/0021-9991(82)90058-4Search in Google Scholar

[15] Hughes, T.J.R. (1995). Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods. Computer Methods in Applied Mechanics and Engineering, 127, 387–401. 10.1016/0045-7825(95)00844-9Search in Google Scholar

[16] Kumar, S., Kumar, B. V., & Rai, S. K. (2023). Influence of abdominal aortic aneurysm shape on hemodynamics in human aortofemoral arteries: A transient open-loop study. Physics of Fluids, 35(4), 041903. 10.1063/5.0139085Search in Google Scholar

[17] Li, X., Sun, B., Zhao, H., Ge, X., Liang, F., Li, X., & Liu, X. (2018). Retrospective study of hemodynamic changes before and after carotid stenosis formation by vessel surface repairing. Scientific Reports, 8(1), 5493. 10.1038/s41598-018-23842-0Search in Google Scholar PubMed PubMed Central

[18] Nguyen, P. H., Coquis-Knezek, S. F., Mohiuddin, M. W., Tuzun, E., & Quick, C. M. (2015). The complex distribution of arterial system mechanical properties, pulsatile hemodynamics, and vascular stresses emerges from three simple adaptive rules. American Journal of Physiology-Heart and Circulatory Physiology, 308(5), H407–H415. 10.1152/ajpheart.00537.2014Search in Google Scholar PubMed

[19] Ponalagusamy, R., & Priyadharshini, S. (2017). Numerical modelling on pulsatile flow of Casson nanofluid through an inclined artery with stenosis and tapering under the influence of magnetic field and periodic body acceleration. Korea-Australia Rheology Journal, 29, 303–316. 10.1007/s13367-017-0030-2Search in Google Scholar

[20] Tang, D., Yang, C., Kobayashi, S., & Ku, D. N. (2001). Generalized finite difference method for 3-D viscous flow in stenotic tubes with large wall deformation and collapse. Applied numerical mathematics, 38(1–2), 49–68. 10.1016/S0168-9274(00)00062-3Search in Google Scholar

[21] Zaman, A., Ali, N., & Sajid, M. (2017). Numerical simulation of pulsatile flow of blood in a porous-saturated overlapping stenosed artery. Mathematics and Computers in simulation, 134, 1–16. 10.1016/j.matcom.2016.09.008Search in Google Scholar

Received: 2023-08-02

Revised: 2023-10-31

Accepted: 2023-12-29

Published Online: 2024-03-20

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/cmb-2023-0118

Keywords for this article

stabilized finite element method; magnetohydrodynamic Navier-Stokes; diseased artery; wall shear stress

Creative Commons

BY 4.0