A numerical method for MHD Stokes model with applications in blood flow

Chitranjan Pandey; B. V. Rathish Kumar

doi:10.1515/cmb-2023-0105

Open Access Published by De Gruyter Open Access November 15, 2023

A numerical method for MHD Stokes model with applications in blood flow

Chitranjan Pandey and B. V. Rathish Kumar

From the journal Computational and Mathematical Biophysics

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

Abstract

The magnetohydrodynamic (MHD) Stokes equations have several applications in the field of biofluid dynamics. In the present study, we propose the staggered finite volume method (S-FVM) for MHD Stokes equations and establish its equivalence to a nonconforming finite element approximation. We also theoretically establish the convergence of the proposed S-FVM. The error estimation is carried out in an unstructured grid framework which is known for its flexibility and robustness in dealing with complex domains. The apriori estimate shows that the L 2 -norm of the error for the pressure and velocity components is of order h , the spacial grid size. After validating the numerical performance of the scheme against benchmark test cases, we do numerical simulations for the blood flow through an injured arteriole and analyze the influence of the magnetic force on hemodynamics in the arteriole under an injured condition.

Keywords: Stokes equations; FVM; Galerkin expansion; error estimation; hemodynamics

MSC 2010: 74S05; 76M12; 46N60

1 Introduction

Blood flows in the human circulatory system is dependent on the heart’s pumping motion, which creates a pressure differential throughout the system. Blood is electrically conducting and has magnetohydrodynamic (MHD) flow properties because it is a suspension of red blood cells containing hemoglobin, which includes iron oxide [20]. Pirmohammadi and Ghassemi [18] studied natural-convection flow on a staggered grid system in the presence of a magnetic field using a hybrid scheme for the convection terms, which was a mix of the central difference scheme and the upwind finite volume scheme. Varshney et al. [21] have used a transverse magnetic field to study the flow of blood with several parameters like blood velocity, flow rate, wall shear stress, flow resistance, and flow acceleration in an artery with multiple stenosis. With the help of several simulations of the blood flow through an injured arteriole of realistic dimension, Ireland and Gent [11] studied the interaction of red blood cells with the surrounding plasma for different severity levels of the injury. On the other hand, MHD Stoke’s equations, Navier-Stokes equations deals with the interaction of magnetic fields with electrically conducting fluids, hence it has many biological and medical applications. So, it has been a topic of interest for many authors, and several studies have been done based on finite element method and finite volume method (FVM) [2,6,15]. Recently, Merve and Tezer-Sezgin [14] studied the lid-driven cavity flow and the backward-facing step flow under different magnetic intensities in the horizontal and vertical directions.

Furthermore, FVM is an efficient method used for the discretization of conservation laws of fluid flow. Consequently, FVM is widely used in solving fluid flow models numerically over regular domains as well as over complex domains [5,8]. Also, unstructured grid, due to its flexibility for fitting an arbitrary geometry, is used efficiently for the discretization of such equations over complex domains [13,22,23]. Staggering of velocity and pressure prognostic variables while solving the fluid flow models is well known to yield desired stability analogous to the up-winding methods, wherein the discrete models resulting from the numerical approximation will be the stable numerical counterparts of their continuous versions. Such stable schemes in an unstructured framework are known to efficiently solve fluid flow in complex geometries, as found in several biological flows like peristaltic flow in the esophagus or intestines, coronary arterial blood flow, or flow in aortic bifurcation [17,19].

A first-order accurate scheme, based on FVM, for steady-state Stokes equation was proposed by Eymard et al. [8], where the pressure and velocity coupling is done with such a staggered way that velocity components are located at cell center, while pressure field is expressed as a Galerkin expansion whose values are located at the vertices of the cell, which was analyzed by Alami-Idrissi and Atounti [1] on a general triangular meshes in point-value sense. Further, it was improved by Blanc et al. [3] on an equilateral mesh in average-value sense. Later on, Nicaise and Djadel [16] validated it is numerical convergence by mapping the scheme to a nonconforming finite element problem on a Ciarlet-regular mesh [4], satisfying the orthogonal property. Derivation of FVM schemes in an unstructured framework with the staggering of field variables to achieve stable numerical methods though complex are known for their flexibility and stability. Largely, the available literature on numerical methods based on FVM deal with the efficient implementations and conservation properties of the scheme. However, the question of convergence analysis and the theoretical establishment of these widely used schemes remain widely open.

Here, we theoretically analyze the proposed FVM scheme for a multidimensional MHD Stokes model on an unstructured grid made up of triangles/tetrahedrons and derive the corresponding aprior estimates in the functional analytic framework to theoretically prove the convergence of the scheme and also establish the nonconforming finite element approximation of the same scheme. While there are few studies on blood flow in arteriole [11], very little has been reported on injured arterioles. In particular, no investigation on magnetic force effect on the flow dynamics in such a small arteriole is available. Such a study on the MHD effect on blood flow, especially in an injured arteriole, is useful for providing insight into the effective drug delivery. Hence, here, we carry out numerical simulations for the blood flow through an injured arteriole under the influence of the magnetic forcing. The complete study is organized in the following way: in Section 2, we describe the governing equation; in Section 3, we define a regular finite volume mesh and propose a finite volume discretization of the governing equations (1)–(3); in Section 4, we analyze the convergence of the proposed staggered FVM scheme and discuss the nonconforming finite element approximation of the scheme; in Section 5, we present numerical validation of the scheme and later analyze the blood flow in an injured arteriole; and subsequently, we conclude in Section 6.

2 Governing equation

2.1 The continuous MHD Stokes equations

Let Ω be an open-bounded polygonal subset of R d , d = 2 or 3. The MHD Stokes equations with a magnetic reaction term is given as follows:

(1) − μ Δ u i ( x ) + ∂ p ( x ) ∂ x i + b i u i ( x ) = f i ( x ) , i = 1 , … , d

(2) ∇ ⋅ u = 0 , ∀ x ∈ Ω ,

with the Dirichlet boundary condition,

(3) u ( x ) = 0 , x ∈ ∂ Ω ,

where u = ( u 1 , … , u d ) is the fluid-flow velocity, p is the scalar pressure field, f i s are the source terms, and μ > 0 and b i ≥ 0 be the constants standing for the dynamic viscosity and magnetic force coefficient, respectively. The direction of the magnetic force, horizontal or vertical, depends on the values of b i s [14].

Assumption 2.1

Ω be an open-bounded polygonal subset of R d , d = 2 , or 3.
μ > 0 and b i ≥ 0 .
f i ∈ L 2 ( Ω ) ∀ i = 1 , … , d .

Under Assumption 2.1, the weak formulation of the system of equations (1)–(3),

a ( u , v ) + b ( p , v ) = ( f , v ) ∀ v ∈ V , b ( u , q ) = 0 ∀ q ∈ Q ,

where

V = ( H 0 1 ( Ω ) ) d , and Q = q ∈ L 2 ( Ω ) ∫ Ω q d x = 0

has unique solution [9], ( u , p ) ∈ ( V , Q ) , with

a ( u , v ) = μ ∑ i ∫ Ω ∇ u i . ∇ v i d x + ∑ i b i ∫ Ω u i v i d x ,

b ( u , q ) = − ∫ Ω q div ( u ) d x , and ( f , v ) = ∑ i ∫ Ω f i v i d x .

3 Finite volume discretization of governing equation

3.0.1 Admissible mesh

Let Ω be an open-bounded polygonal subset of R d , ( d = 2 or 3). An admissible finite volume mesh of Ω , denoted by T , is given by a family of control volumes (CVs), which are open polygonal convex subsets of Ω ¯ contained in hyperplanes of R d , denoted by ℰ , with strictly positive ( d − 1 ) -dimensional measure and a family of points of Ω denoted by P satisfying the following properties:

The closure of the union of all CVs is Ω ¯ .
For any K ∈ T , there exists a subset ℰ K of ℰ such that ∂ K = K ¯ \ K = ∪ σ ∈ ℰ K σ ¯ , and let ℰ = ∪ K ∈ T ℰ K .
For any ( K , L ) ∈ T 2 , with K ≠ L , either the d -dimensional Lebesgue measure of K ¯ ∩ L ¯ is 0 or K ¯ ∩ L ¯ = σ ¯ for some σ ∈ ℰ .
The family P = ( x K ) K ∈ T is such that x K ∈ K ¯ , and if σ = K ∣ L , it is assumed that x k ≠ x L and that the straight line D K , L going through x K and x L is orthogonal to K ∣ L .
For any σ ∈ ℰ such that σ ∈ ∂ Ω , let K be the CV such that σ ∈ ℰ K , and if x K ∉ σ , let D K , σ be the straight line going through x K and orthogonal to σ . Then, the condition D K , σ ∩ σ ≠ ∅ is assumed and let y σ = D K , σ ∩ σ .

Furthermore, the following notations are used:

size ( T ) = sup { diam ( K ) , K ∈ T } .
m ( K ) is the d -dimensional Lebesgue of K for any K ∈ T .
m ( σ ) is the ( d − 1 ) -dimensional Lebesgue of σ for any σ ∈ ℰ .
ℰ int = { σ ∈ ℰ , σ ⊄ ∂ Ω } and ℰ ext = { σ ∈ ℰ , σ ⊂ ∂ Ω } .
If σ ∈ ℰ int , σ = K ∣ L , then d σ = d K ∣ L = d ( x K , x L ) , and if σ ∈ ℰ K ∩ ℰ ext , then d σ = d K , σ = d ( x K , y σ ) .
For any σ ∈ ℰ , the transmissibility through σ is defined by τ σ = m ( σ ) d σ if, d σ ≠ 0 and τ σ = 0 if d σ = 0

Lemma 3.1

(Galerkin interpolation of the pressure field) Let Ω be the open bounded polygonal subset of R d , ( d = 2 or 3). and T be the admissible finite volume mesh (discussed in Section 3.0.1) made up of triangles/tetrahedrons which satisfy the required properties for finite element method with the regularity in Ciarlet’s sense, Nicaise and Djadel [16]. Let s T denote the set of vertices of T ( s K be that of a triangle/tetrahedron K ∈ T ) and ϕ s be the linear shape functions at the vertex s ∈ s T . If p K = ∑ s ∈ S K p s ϕ s be the local interpolation of pressure field over the triangle/tetrahedron K ∈ T , then the global interpolation error ‖ p − p T ‖ L 2 and ∣ p − p T ∣ H 1 will be given as follows [12]:

‖ p − p T ‖ L 2 ( Ω ) = ∑ K ∈ T ‖ p − p K ‖ L 2 ( K ) 2 1 ⁄ 2 ≤ C 1 h ∣ p ( x ) ∣ H 2 ( Ω ) ,

∣ p − p T ∣ H 1 ( Ω ) ≤ ∑ K ∈ T ∣ p − p K ∣ H 1 ( K ) 2 1 ⁄ 2 ≤ C 2 h ∣ p ( x ) ∣ H 2 ( Ω ) ,

for some C 1 , C 2 > 0 , where

∣ p ( x ) ∣ H r = ∑ ∣ α ∣ = r ∫ Ω ∣ D r p ( x ) ∣ 2 d x 1 ⁄ 2 .

The global pressure field is given as follows:

p T = ∑ S T p s ϕ s .

3.1 Discrete formulation

Let T = { K , σ } be the admissible finite volume mesh (Section 3.0.1) consisting of triangles satisfying the properties required for the finite element method. Then, integration of equation (1) over a CV K ∈ T gives

(4) − μ ∫ ∂ K ∇ u i . n K d γ ( x ) + ∫ K ∂ p ∂ x i d x + ∫ K b i u i d x = ∫ K f i d x .

A possible finite volume scheme using Galerkin expansion for the pressure field is defined by the following equations:

(5) − μ ∑ σ ∈ ℰ K F K , σ i + ∑ s ∈ S K p s ∫ K ∂ ϕ s ∂ x i d x + b i m ( K ) u K i = m ( K ) f K i ,

where

F K , σ i = T σ ( u L i − u K i ) , if σ ∈ ε int , σ = K ∣ L . − T σ u K i , if σ ∈ ε ext ∩ ε K ∀ K ∈ T .

(6) ∑ i = 1 d ∑ K ∈ T u K i ∫ K ∂ ϕ s ∂ x i d x = 0 ∀ s ∈ S T .

(7) ∫ Ω ∑ s ∈ S T p s ϕ s d x = 0 ,

where

f K i = 1 m ( K ) ∫ K f ( x ) d x ∀ K ∈ T , i = 1 , … , d .

The unknowns are u K i , ∀ K ∈ T , i = 1 , … , d , and p s , s ∈ S T .

4 A finite-volume error estimate

Definition 4.1

Let Ω be an open-bounded polygonal subset of R d , ( d = 2 or 3) and T be an admissible mesh (Section 3.0.1). Define X ( T ) to be the set of piecewise constant functions from Ω to R , i.e.,

u ∈ X ( T ) ⇒ u ( x ) = u K ∀ x ∈ K . ∀ K ∈ T

Definition 4.2

(Mesh-dependent norm) Let Ω be an open-bounded polygonal subset of R d , ( d = 2 , or 3) and T be an admissible mesh (Section 3.0.1). For u ∈ X ( T ) , L 2 -norm and discrete H 0 1 -norm are defined as follows:

‖ u ‖ L 2 ( Ω ) = ∫ Ω ∣ u ∣ 2 d x 1 ⁄ 2 = ∑ K ∈ T ∫ K ∣ u ∣ 2 d x 1 ⁄ 2 = ∑ K ∈ T m ( K ) ∣ u K ∣ 2 1 ⁄ 2

and,

‖ u ‖ 1 , T = ∑ σ ∈ ℰ τ σ ( D σ u ) 2 1 2 ,

where

D σ u = ∣ u K − u L ∣ if σ ∈ ℰ int σ = K ∣ L , D σ u = ∣ u K ∣ if σ ∈ ℰ ext ∩ ℰ K ,

and u K denotes the value taken by u on the CV K .

Lemma 4.1

(Discrete Poincaré inequality) Let Ω be an open-bounded polygonal subset of R d ( d = 2 or 3) and T be an admissible mesh (Section 3.0.1). Then, any u ∈ X ( T ) satisfies:

‖ u ‖ L 2 ( Ω ) ≤ diam ( Ω ) ‖ u ‖ 1 , T .

Theorem 4.2

Let Ω be an open-bounded polygonal subset of R d , d = 2 or 3, and T be the admissible mesh (Section 3.0.1) satisfying the assumptions in Lemma 3.1. Let u i ∈ C 2 ( Ω ¯ , R ) , i = 1 , … , d , and p ∈ C 1 ( Ω ¯ , R ) , and f i ∈ C 0 ( Ω ¯ , R ) is defined by f i = − μ Δ u i + p x i + b i u i , ∀ i = 1 , … , d . Then, equations (5)–(7) have unique solution u T i = u K i , K ∈ T , i = 1 , … , d , and p s , s ∈ S T such that u T i satisfies the following estimate:

∑ i = 1 d ‖ u T i ‖ L 2 2 1 ⁄ 2 = ‖ u T ‖ L 2 ≤ B

for some B ∈ R + , depending on Ω , μ , and f i , i = 1 , … , d .

Moreover, the error for the velocity component e T i ∈ X ( T ) such that e T i ( x ) = e K i = u ( x K ) − u K i ∀ a.e. x ∈ K ∈ T , i = 1 , … , d , satisfies the following estimate:

∑ i = 1 d ‖ e T i ‖ L 2 2 1 ⁄ 2 = ‖ e T ‖ L 2 ≤ C . s i z e ( T )

for some C ∈ R + , depending on Ω , u i , p , μ , and b i , where u T = ( u T i ) i = 1 d , and e T = ( e T i ) i = 1 d .

Proof

Let f K i = 0 and u K i = 0 ∀ K ∈ T , i = 1 , … , d . Now, multiplying equation (5) by u K i and summing over all i = 1 , … , d and K ∈ T , we have the following:

μ ‖ u T ‖ 1 , T 2 + ∑ i b i ‖ u T i ‖ L 2 2 = 0 .

As μ > 0 and b i ≥ 0 , above relation gives u K i = 0 , i = 1 , … , d . Now, if we write the right hand side (RHS) of equation (6) by g s ∀ s ∈ S T , then

∑ s g s = 0

gives the sufficient condition for the existence of solution for equations (5) and (6), and including mean pressure zero (7) will provide the unique solution for equations (5)–(7) [8, pp. 219–222].

Clearly, multiplying equation (5) by u i and summing over all i = 1 , … , d and K ∈ T and using equation (6), we have:

μ ‖ u T ‖ 1 , T 2 + ∑ i b i ‖ u T i ‖ L 2 2 = ∑ i ∑ K m ( K ) f K i u K i ,

which gives the estimate for u T on the application of discrete Holder’s inequality on the RHS of above equation and then using Lemma 4.1.

Now, adopting the notations,

F ¯ K , σ i = − ∫ σ ∇ u i ( x ) . n K , σ d γ ( x ) ,

ρ K i = u i ( x K ) − 1 m ( K ) ∫ K u i ( x ) d x ,

F K , σ * i = − T σ ( u i ( x L ) − u i ( x K ) ) , if σ ∈ ε int , σ = K ∣ L , T σ u i ( x K ) , if σ ∈ ε ext ∩ ε K ,

∀ K ∈ T , i = 1 , … , d .

Hence,

F ¯ K , σ i − F K , σ i = F ¯ K , σ i − F K , σ * i + F K , σ * i − F K , σ i = m ( σ ) R K , σ i + T σ ( e K i − e L i ) , if σ ∈ ε int , σ = K ∣ L , m ( σ ) R K , σ i + T σ e K i , if σ ∈ ε ext ∩ ε K ,

The numerical flux satisfies the following consistency estimate.

Lemma 4.3

(Consistency of the flux) Let T be the admissible mesh (Section 3.0.1) and u i ∈ C 2 ( Ω ¯ , R ) , i = 1 , … , d , then there exists a constant, say, c ∈ R + such that:

( F ¯ K , σ i − F K , σ * i ) 1 m ( σ ) = R K , σ i , w i t h ∣ R K , σ i ∣ ≤ c size ( T ) .

Moreover,

∣ ρ K i ∣ ≤ c 1 size ( T )

for some c 1 ∈ R + .

Now, equation (4) can be rewritten as follows:

(8) μ ∑ σ ∈ ℰ K F ¯ K , σ i + ∫ K ∂ p ( x ) ∂ x i d x + b i m ( K ) u i ( x K ) = m ( K ) f K i + b i m ( K ) ρ K i .

Subtracting equation (5) from equation (8) and using the consistency result from Lemma 4.3, we have

(9) μ ∑ σ ∈ ℰ K G K , σ i + ∫ K ∂ p ( x ) ∂ x i d x − ∑ s ∈ S K p s ∫ K ∂ ϕ s ∂ x i + b i m ( K ) e K i = b i m ( K ) ρ K i − μ ∑ σ ∈ ℰ K m ( σ ) R K , σ i ,

where

G K , σ i = − T σ ( e L i − e K i ) , if σ ∈ ε int , σ = K ∣ L , T σ e K i , if σ ∈ ε ext ∩ ε K ,

∀ K ∈ T , i = 1 , … , d .

Now, multiplying equation (9) by e K i and summing over all K ∈ T , then over all i = 1 , … , d gives-

(10) μ ∑ i = 1 d ∑ K ∈ T ∑ σ ∈ ℰ K G K , σ i e K i + ∑ i = 1 d ∑ K ∈ T ∫ K ∂ p ( x ) ∂ x i d x − ∑ s ∈ S K p s ∫ K ∂ ϕ s ∂ x i e K i + ∑ i = 1 d ∑ K ∈ T m ( K ) b i ( e K i ) 2 = ∑ i = 1 d ∑ K ∈ T b i m ( K ) ρ K i e K i − μ ∑ i = 1 d ∑ K ∈ T ∑ σ ∈ ℰ K m ( σ ) R K , σ i e K i .

Now, we solve equation (10) in several parts.

Applying Hölder’s inequality in the last term of equation (10), we obtain
∑ K ∈ T ∑ σ ∈ ℰ K m ( σ ) R K , σ i e K i ≤ ∑ σ ∈ ℰ K m ( σ ) R σ i ∣ D σ e K i ∣ ≤ ∑ σ ∈ ℰ K m ( σ ) d σ ∣ D σ e i ∣ 2 1 ⁄ 2 ∑ σ ∈ ℰ K m ( σ ) d σ ∣ R σ i ∣ 2 1 ⁄ 2 .
Now, using Lemma 4.3 and the fact that ∑ σ ∈ ℰ K m ( σ ) d σ = d . m ( Ω ) , d > 0 , the above inequality becomes-
∑ K ∈ T ∑ σ ∈ ℰ K m ( σ ) R K , σ i e K i ≤ c 3 h m ( Ω ) 1 ⁄ 2 ‖ e T i ‖ 1 , T ∀ i = 1 , … , d
for some c 3 ∈ R + .
Using Lemma 4.3 and Young’s inequality, we obtain
∑ K ∈ T b i m ( k ) ρ K i e K i ≤ ∑ K ∈ T c 1 b i h m ( K ) ∣ e K i ∣ ≤ c 1 2 h 2 b i m ( Ω ) 2 + b i ‖ e T i ‖ L 2 2 2 .
∫ K ∂ p ( x ) ∂ x i d x − ∑ s ∈ S K p s ∫ K ∂ ϕ s ∂ x i d x e K i ≤ ∫ K ∂ ( p − p K ) ∂ x i 2 d x 1 ⁄ 2 ∫ K ∣ e K i ∣ 2 d x 1 ⁄ 2 .

Hence, using Lemma 3.1, we have c 4 > 0 , depending on pressure field p so that

∑ i = 1 d ∑ K ∈ T ∫ K ∂ p ( x ) ∂ x i d x − ∑ s ∈ S K p s ∫ K ∂ ϕ s ∂ x i d x e K i ≤ c 4 h ‖ e T ‖ L 2 .

Now, using (A), (B), and (C) in equation (10), we obtain

(11) μ ‖ e T ‖ 1 , T 2 + ∑ i d b i ‖ e T i ‖ L 2 2 ≤ ∑ i d c 2 2 h 2 b i m ( Ω ) 2 + b i ‖ e T i ‖ L 2 2 2 + μ c 3 h m ( Ω ) 1 ⁄ 2 ‖ e T ‖ 1 , T + c 4 h ‖ e T ‖ L 2 .

Using Young’s inequality in the last term of RHS, we obtain

μ c 3 h m ( Ω ) 1 ⁄ 2 ‖ e T ‖ 1 , T ≤ μ c 3 2 h 2 m ( Ω ) 2 + μ ‖ e T ‖ 1 , T 2 2 .

So, equation (11) becomes (as b , μ ≥ 0 )

(12) ‖ e T ‖ 1 , T 2 ≤ ∑ i d c 2 2 h 2 b i m ( Ω ) μ + c 3 2 h 2 m ( Ω ) + 2 c 4 h ‖ e T ‖ L 2 μ .

Now, applying Lemma 4.1 and then Young’s inequality on the last term of equation (12), we have

2 c 4 h μ ‖ e T ‖ L 2 ≤ 2 c 4 h μ ( diam ( Ω ) ‖ e T ‖ 1 , T ) ≤ c 4 2 h 2 diam ( Ω ) 2 μ 2 ε + ε ‖ e T ‖ 1 , T 2

choosing ε = 1 ⁄ 2 and rewriting C 2 2 = max 2 ∑ i d c 2 2 b i m ( Ω ) μ , 2 c 3 2 m ( Ω ) , 4 c 4 2 h 2 diam ( Ω ) 2 μ 2 so that C 2 depends on μ , b i , u i , Ω , and p . Now, equation (12) can be written as follows:

(13) ‖ e T ‖ 1 , T 2 ≤ C 2 2 h 2 .

Applying the discrete Poincaré inequality 4.1 in equation (13), one can have the required results.□

4.1 Nonconforming finite element method approximation of the scheme

Nicaise and Djadel [16] validated the numerical convergence of FVM scheme for Stokes equations by mapping the scheme to a nonconforming finite element problem. Now, we consider a nonconforming approximation V h of the space V defined in Section 2. We approximate the problem (4) in the spaces V h and Q h defined by the following:

V h = { u h ∈ ( L 2 ( Ω ) ) d : u h ∣ K ∈ ( P 0 ( K ) ) d , ∀ K ∈ T }

and

Q h = { q h ∈ Q ∩ C ( Ω ¯ ) : q h ∣ K ∈ ( P 1 ( K ) ) , ∀ K ∈ T }

with the norm defined as follows:

‖ u h ‖ T = ∑ σ ∈ ℰ τ σ ( D σ u h ) 2 1 2 ,

where

D σ u h = ∣ u K − u L ∣ if σ ∈ ℰ int , σ = K ∣ L , D σ u h = ∣ u K ∣ if σ ∈ ℰ ext ∩ ℰ K .

For any v h = ( v h 1 , … , v h d ) , w h = ( w h 1 , … , w h d ) ∈ V h , and q h ∈ Q h , we define

a h ( v h , w h ) ≔ − μ ∑ K ∈ T ∑ σ ∈ ℰ K F K , σ ( v h ) . w h + ∑ i = 1 d b i v h i w h i

b h ( v h , q h ) ≔ ∫ Ω v h . ∇ q h d x

( f , w h ) ≔ ∫ Ω f . w h d x ,

where

F K , σ ( v h ) = T σ ( v L − v K ) , if σ ∈ ε int , σ = K ∣ L , − T σ v K , if σ ∈ ε ext ∩ ε K , ∀ K ∈ T

Lemma 4.4

{ u K = ( u K 1 , … , u K d ) } K ∈ T and { p s } s ∈ S is a solution of the system of equations (5)–(7) if and only if u h = ∑ K ∈ T u K χ K ∈ V h and p h = ∑ s ∈ S p s ϕ s is a solution of

(14) a h ( u h , w h ) + b h ( w h , p h ) = ( f , w h ) , ∀ w h ∈ V h b h ( u h , q h ) = 0 , ∀ q h ∈ Q h .

5 Numerical simulation

5.1 Code validation

We consider the lid-driven cavity flow under the influence of the vertically applied magnetic forcing of different magnitudes. We validate the proposed scheme (system of equations (5)–(7)) or, equivalently, the scheme mentioned in Lemma 4.4 [16] by comparing the numerical results with that of Merve and Tezer-Sezgin [14]. All numerical simulations are done on FreeFem++ [10].

We set a square domain of unit sides and examine Stokes flow under the influence of the magnetic field applied in the vertical direction. The domain is triangulated with a grid resolution of 100 × 100 . We set the boundary conditions as follows: the upper wall is kept at a constant horizontal velocity of 1 unit and all the other three walls are at no-slip condition. b 1 = 0 , f 1 = 0 , f 2 = 0 , and b 2 = H a 2 as taken by Merve and Tezer-Sezgin [14]. We did three simulations for different Hartmann numbers ( H a ) , namely, Ha = 1 , Ha = 30 , and Ha = 80 . We compare the numerically simulated horizontal velocity ( u 1 ) at the vertical line segment x = 0.5 with that of Merve and Tezer-Sezgin [14]. Figure 1 shows that the numerical results are in good agreement with the results obtained by Nicaise and Djadel [14].

Figure 1

Comparison of the horizontal velocity at x = 0.5 .

5.2 Hemodynamic simulations

We take a realistic geometry to mimic an arteriole as it was taken by Ireland and Gent [11] to analyze the hemodynamic process around the injury. The arteriole’s diameter is fixed at 0.05 mm. The arteriole’s entire length is 1.00 mm. The arteriole damage was depicted by a rectangular laceration in the middle of the length of the arteriole. The laceration was estimated to be 0.10 mm in length and 0.02 mm in depth, so it represents an injury inflicted by a razor of the thickness range of 0.1 mm. In addition, the no leakage condition will represent a sharp stenosis case. We use a grid resolution of 80 × 10 for all the hemodynamic simulations presented in Section 5.2. We set inlet velocity and the Reynolds number according to Doutel et al. [7]. We did several numerical simulations to observe the influence of magnetic force near the injury. We set the inlet as follows: u has a parabolic profile with maximum velocity 1 mm/s and v = 0 . While at the outlet, ∂ u ∂ n = 0 , v = 0 , and p = 0 . The Reynolds number (Re) is chosen to be 0.7 [7]. We observe the effect of the magnitude of magnetic field applied in the positive x -direction in terms of the following variables: the horizontal velocity, the pressure field, and the wall-shear-stress (WSS) on the top and bottom walls of the arteriole.

Figures 2 and 6 present the u -velocity contour and pressure field, respectively, in the presence of no magnetic force. The horizontal flow velocity is found to be relatively higher around the injury to maintain the flow. Furthermore, in the presence of very small magnetic force in the horizontal direction, Ha = 10 , i.e., b 1 = 1 0 2 and b 2 = 0 , the flow slows down, and the pressure increases a little, which can be seen in Figures 3 and 7. A slight decay in the magnitude of horizontal velocity is observed in the presence of relatively high-intense magnetic force and can be seen in Figures 2–5. As the horizontal magnetic force generates a drag force against the flow, the pressure inside the arteriole as well as the WSS on the walls of the arteriole increase, as shown in Figures 6–9. Figures 10 and 11 illustrate the WSS plots for the top and bottom walls, respectively, for various magnetic force intensities in the horizontal direction. The WSS, on the top and the bottom wall, is found to be relatively higher around the injury region even in the absence of the magnetic force due to the blockage. Furthermore, in the presence of magnetic force, a considerable increment in WSS is observed at the walls of the arteriole. The stronger the magnetic force, the stronger the generated drag force opposite to the flow. Hence, WSS is increased on the entire top and bottom wall. A smooth variation in WSS is observed on the bottom wall around the injury region, while the variation is sharp on the top wall due to the sharp edges of the damage.

$Figure 2 Contour representation of u u -velocity for Ha = 0 {\rm{Ha}}=0 .$

Figure 2

Contour representation of u -velocity for Ha = 0 .

$Figure 3 Contour representation of u u -velocity for Ha = 10 {\rm{Ha}}=10 .$

Figure 3

Contour representation of u -velocity for Ha = 10 .

$Figure 4 Contour representation of u u -velocity for Ha = 50 {\rm{Ha}}=50 .$

Figure 4

Contour representation of u -velocity for Ha = 50 .

$Figure 5 Contour representation of u u -velocity for Ha = 100 {\rm{Ha}}=100 .$

Figure 5

Contour representation of u -velocity for Ha = 100 .

$Figure 6 Contour representation of pressure ( × 1 0 − 3 \times 1{0}^{-3} ) field for Ha = 0 {\rm{Ha}}=0 .$

Figure 6

Contour representation of pressure ( × 1 0 − 3 ) field for Ha = 0 .

$Figure 7 Contour representation of pressure ( × 1 0 − 3 \times 1{0}^{-3} ) field for Ha = 10 {\rm{Ha}}=10 .$

Figure 7

Contour representation of pressure ( × 1 0 − 3 ) field for Ha = 10 .

$Figure 8 Contour representation of pressure ( × 1 0 − 3 \times 1{0}^{-3} ) field for Ha = 50 {\rm{Ha}}=50 .$

Figure 8

Contour representation of pressure ( × 1 0 − 3 ) field for Ha = 50 .

$Figure 9 Contour representation of pressure ( × 1 0 − 3 \times 1{0}^{-3} ) field for Ha = 100 {\rm{Ha}}=100 .$

Figure 9

Contour representation of pressure ( × 1 0 − 3 ) field for Ha = 100 .

$Figure 10 WSS on the top wall for different Ha values. Horizontal axis represents the length of the arteriole and the vertical axis, WSS × 1 0 − 2 . {\rm{WSS}}\times 1{0}^{-2}.$

Figure 10

WSS on the top wall for different Ha values. Horizontal axis represents the length of the arteriole and the vertical axis, WSS × 1 0 − 2 .

$Figure 11 WSS on the bottom wall for different Ha values. Horizontal axis represents the length of the arteriole and the vertical axis, WSS × 1 0 − 2 . {\rm{WSS}}\times 1{0}^{-2}.$

Figure 11

WSS on the bottom wall for different Ha values. Horizontal axis represents the length of the arteriole and the vertical axis, WSS × 1 0 − 2 .

6 Conclusion

A finite volume scheme for MHD Stokes equations on an unstructured regular mesh of triangles/tetrahedron has been proposed, and furthermore, its theoretical convergence is studied. The equivalence of a nonconforming finite-element approximation with the proposed finite volume scheme has been established. The numerical validation of the scheme is done with the lid-driven cavity flow for the different intensities of the magnetic forces. Furthermore, hemodynamic simulations in an injured arteriole have been carried out under the influence of different magnetic intensities, and the blood flow parameters have been studied. The flow slows down in the presence of magnetic force but both the blood pressure and the wall-shear-stress tend to increase. Although the blood flow deceleration is helpful in the drug intake at the injured site of the arteriole, the increase in WSS is not favorable for the health of the blood cell.

Funding information: This research received no specific grant from any funding agency and commercial or nonprofit sectors.
Conflict of interest statement: Authors state no conflict of interest.
Ethical approval: This research did not required ethical approval.
Data availability statement: Data sharing is not applicable.

References

[1] Alami-Idrissi, A., & Atounti, M. (2002). An error estimate for finite volume methods for the Stokes equations. Journal of Inequalities in Pure and Applied Mathematics, 3(1), 17.Search in Google Scholar

[2] Anwar, O., & Ghosh, B. (2011). Applied magnetofluid dynamics: Modelling and computation. Germany: Lap Lambert Academic Publ. Search in Google Scholar

[3] Blanc, P., Eymard, R., & Herbin, R. (2004). An error estimate for finite volume methods for the Stokes equations on equilateral triangular meshes. Numerical Methods for Partial Differential Equations: An International Journal, 20(6), 907–918. 10.1002/num.20020Search in Google Scholar

[4] Ciarlet, P. G. (1978). Interpolation error estimates for the reduced Hsieh-Clough-Tocher triangle. Mathematics of Computation, 32(142), 335–344. 10.1090/S0025-5718-1978-0482249-1Search in Google Scholar

[5] Dalal, A., Eswaran, V., & Biswas, G. (2008). A finite-volume method for Navier-Stokes equations on unstructured meshes. Numerical Heat Transfer, Part B: Fundamentals, 54(3), 238–259. 10.1080/10407790802182653Search in Google Scholar

[6] Das, C., Wang, G., & Payne, F. (2013). Some practical applications of magnetohydrodynamic pumping. Sensors and Actuators A: Physical, 201, 43–48. 10.1016/j.sna.2013.06.023Search in Google Scholar

[7] Doutel, E., Galindo-Rosales, F. J., & Campo-Deannno, L. (2021). Hemodynamics challenges for the navigation of medical microbots for the treatment of CVDs. Materials, 14(23), 7402. 10.3390/ma14237402Search in Google Scholar PubMed PubMed Central

[8] Eymard, R., Gallouët, T., & Herbin, R. (2000). Finite volume methods. Handbook of Numerical Analysis, 7, 713–1018. 10.1016/S1570-8659(00)07005-8Search in Google Scholar

[9] Girault, V., & Raviart, P.-A. (1979). Finite element approximation of the Navier-Stokes equations (Vol. 749). Berlin: Springer. 10.1007/BFb0063447Search in Google Scholar

[10] Hecht, , F. (2012). New development in FreeFem++. Journal of Numerical Mathematics, 20(3–4), 251–266. https://freefem.org/. 10.1515/jnum-2012-0013Search in Google Scholar

[11] Ireland, S. L., & Gent, S. P. (2016). Computational fluid dynamics simulations of blood flow through an injured arteriole. Journal of Medical Devices, 10(3), 030939. 10.1115/1.4033868Search in Google Scholar

[12] Johnson, C. (2012). Numerical solution of partial differential equations by the finite element method. North Chelmsford: Courier Corporation. Search in Google Scholar

[13] Kang, S., & Kim, Y. (2002). Pressure-based unstructured-grid finite-volume method for simulating laminar reacting flows. Numerical Heat Transfer: Part B: Fundamentals, 41(1), 53–72. 10.1080/104077902753385009Search in Google Scholar

[14] Merve, G., & Tezer-Sezgin, M. (2015). MHD Stokes flow in lid-driven cavity and backward-facing step channel. European Journal of Computational Mechanics, 24(6), 279–301. 10.1080/17797179.2016.1181031Search in Google Scholar

[15] Naceur, S., Kadid, F. Z., & Abdessemed, R. (2014). A solution of two-dimensional magnetohydrodynamic flow using the finite volume method. SJEE, 11(2), 201–211. 10.2298/SJEE130520017NSearch in Google Scholar

[16] Nicaise, S., & Djadel, K. (2005). Convergence analysis of a finite volume method for the Stokes system using non-conforming arguments. IMA Journal of Numerical Analysis, 25(3), 523–548. 10.1093/imanum/dri007Search in Google Scholar

[17] Pandey, R., Kumar, M., & Srivastav, V. K. (2020). Numerical computation of blood hemodynamic through constricted human left coronary artery: Pulsatile simulations. Computer Methods and Programs in Biomedicine, 197, 105661. 10.1016/j.cmpb.2020.105661Search in Google Scholar PubMed

[18] Pirmohammadi, M., & Ghassemi, M. (2009). Effect of magnetic field on convection heat transfer inside a tilted square enclosure. International Communications in Heat and Mass Transfer, 36(7), 776–780. 10.1016/j.icheatmasstransfer.2009.03.023Search in Google Scholar

[19] Politis, A., Stavropoulos, G., Christolis, M., Panagopoulos, P., Vlachos, N., & Markatos, N. (2008). Numerical modelling of simulated blood flow in idealized composite arterial coronary grafts: Transient flow. Journal of Biomechanics, 41(1), 25–39. 10.1016/j.jbiomech.2007.08.007Search in Google Scholar PubMed

[20] Tzirtzilakis, E. (2005). A mathematical model for blood flow in magnetic field. Physics of Fluids, 17(7), 077103. 10.1063/1.1978807Search in Google Scholar

[21] Varshney, G., Katiyar, V., & Kumar, S. (2010). Effect of magnetic field on the blood flow in artery having multiple stenosis: A numerical study. International Journal of Engineering, Science and Technology, 2(2), 967–82. 10.4314/ijest.v2i2.59142Search in Google Scholar

[22] Vidović, D., Segal, A., & Wesseling, P. (2004). A superlinearly convergent finite volume method for the incompressible Navier-Stokes equations on staggered unstructured grids. Journal of Computational Physics, 198(1), 159–177. 10.1016/j.jcp.2004.01.005Search in Google Scholar

[23] Xie, B., Ii, S., Ikebata, A., & Xiao, F. (2014). A multi-moment finite volume method for incompressible navier-stokes equations on unstructured grids: Volume-average/point-value formulation. Journal of Computational Physics, 277, 138–162. 10.1016/j.jcp.2014.08.011Search in Google Scholar

Received: 2023-07-31

Revised: 2023-10-04

Accepted: 2023-10-19

Published Online: 2023-11-15

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

A numerical method for MHD Stokes model with applications in blood flow

Abstract

1 Introduction

2 Governing equation

2.1 The continuous MHD Stokes equations

Assumption 2.1

3 Finite volume discretization of governing equation

3.0.1 Admissible mesh

Lemma 3.1

3.1 Discrete formulation

4 A finite-volume error estimate

Definition 4.1

Definition 4.2

Lemma 4.1

Theorem 4.2

Proof

Lemma 4.3

4.1 Nonconforming finite element method approximation of the scheme

Lemma 4.4

5 Numerical simulation

5.1 Code validation

5.2 Hemodynamic simulations

6 Conclusion

References

Journal and Issue

Articles in the same Issue