Numerical Investigation of MHD Pulsatile Flow of Micropolar Fluid in a Channel with Symmetrically Constricted Walls

: This article presented an analysis of the pulsatile ﬂow of non-Newtonian micropolar (MP) ﬂuid under Lorentz force’s effect in a channel with symmetrical constrictions on the walls. The governing equations were ﬁrst converted into the vorticity–stream function form, and a ﬁnite difference-based solver was used to solve it numerically on a Cartesian grid. The impacts of different ﬂow controlling parameters, including the Hartman number, Strouhal number, Reynolds number, and MP parameter on the ﬂow proﬁles, were studied. The wall shear stress (WSS), axial, and micro-rotation velocity proﬁles were depicted visually. The streamlines and vorticity patterns of the ﬂow were also sketched. It is evident from the numerical results that the ﬂow separation region near constriction as well as ﬂattening of the axial velocity component is effectively controlled by the Hartmann number. At the maximum ﬂow rate, the WSS attained its peak. The WSS increased in both the Hartmann number and Reynolds number, whereas it declined with the higher values of the MP parameter. The micro-rotation velocity increased in the Reynolds number, and it declined with increment in the MP parameter. The study found applications in understating the blood ﬂow, modeled as non-Newtonian micropolar ﬂuid, in stenotic arteries especially. The outcomes could be used in designing the biomedical devices and techniques for cardiovascular treatments.


Introduction
MP fluids are non-Newtonian fluids consisting of the dilute suspension with an individual motion of thin, rigid cylindrical macromolecules. Incompressible MP fluids have significance in the study of various phenomena such as blood rheology in medical sciences and melted plastic mechanics in industries. MP fluid theory explains the microrotation effects. Eringen [1] first described micro-inertia effects. Several numerical studies have been conducted by researchers to study the behavior of internal and external MP fluid flows. Agarwal et al. [2] examined MP fluid flow on a porous stationary surface with heat transfer. The 2D stagnation point flow of MP fluid for the steady case over a stretching sheet was examined by Nazar et al. [3]. Lok et al. [4] researched the steady mixed convection boundary layer flow of MP fluid on a double-infinite, vertical flat plate near the stagnation point. The flow behavior and heat transfer effects of mixed convection in MP fluid flow over a vertical flat plate with conduction were analyzed by Chang et al. [5]. Magyari et al. [6] examined the flow of quiescent MP fluid over a doubly infinite plate accelerated from rest to a constant velocity. The impacts of radiation and viscous dissipation on MP fluid stagnation-point flow to a nonlinearly stretching surface with suction and injection were reported by Babu et al. [7]. The flow of MP fluid over a porous stretch surface with heat transfer was analyzed by Turkyilmazoglu [8]. Waqas et al. [9] provided a mixed convection flow of MP liquid in the occurrence of the magnetic field on a nonlinear stretched surface. Ramadevi et al. [10] carried out an analysis of the nonlinear MHD radiative flow of MP fluid on a stretching surface.
MHD MP fluid over an oscillating, infinite vertical plate embedded in a porous medium was analyzed by Sheik et al. [11]. Hussanan et al. [12] analytically examined MP fluid flow over a vertical plate with Newtonian heating in the presence of the magnetic field and the absence of thermal radiation. Kumar et al. [13] examined the heat transfer mechanism with variable heat sink/source, a nonlinear approximation of Rosseland and Biot number over a stretched field. Shamshuddin et al. [14] used the finite-element approach for solving MHD, incompressible, dissipative, and chemically reacting MP fluid flow with heat transfer as well as mass transfer on an inclined heat source/sink plate. Nadeem et al. [15] examined the flow of MP fluid over the Riga plate with exponential surface temperature and heating effects.
Si et al. [16] investigated the behavior of MP fluid flow in a porous channel with mutable walls. Lu et al. [17] considered the 2D creeping flow of MP fluid in a thin permeable channel with a variable absorption rate. Fakour et al. [18] studied heat and mass transfer of MP fluid flow inside a channel with permeable walls. Tutty [19] investigated the nonuniform channel which is used as a simple model of a constricted arterial vessel. There are several studies regarding fluid motion with pulsation in a constricted channel. Peristalsis is a mechanism in which progressive transverse waves produced by flexible channel/tube boundary walls transport the fluid. Peristaltic pumping is also very effective in the design of several biomedical devices for maintaining blood supply during critical operations. Mekheimer et al. [20] studied the effect of an induced magnetic field on the peristaltic transport in a symmetric channel of an incompressible MP conductive fluid. Hayat and Ali [21] examined the peristaltic wave motion for the endoscope impact via the distance among two concentric tubes, finding the inner tube to be rigid when moving outwards to allow the MP fluid to flow.
Under certain physical situations, the behavior of the pulsatile flow of Newtonian and non-Newtonian fluids has been examined, usually with assumptions of long wavelength and low Reynolds number to simplification. A numerical study of the MHD pulsatile flow of Newtonian fluid was carried out by Bandyopadhyay and Layek [22] in a singleconstricted channel. Khair et al. [23] described the transition from laminar to the turbulent regime in a constricted channel for pulsatile flow. The steady and pulsatile flow of MHD Casson fluid in a constricted channel was studied by Ali et al. [24].
The present work's objective is to investigate the magnetohydrodynamic (MHD) pulsatile flow of non-Newtonian MP fluid in a channel having symmetrical constrictions on both the walls under the influence of the Lorentz force. The numerical method to solve the governing equations is based on the finite difference method on a Cartesian grid instead of the cylindrical one. The impacts of various parameters on the axial velocity, shear stress, and micro-rotation velocity are discussed. The streamlines and vorticity distributions of the pulsatile MP fluid flow are also shown. The flow separation region generated due to the constriction bumps is also discussed. The flow parameters under consideration for the study include the Hartmann number (M), Strouhal number (St), Reynolds number (Re), and MP parameter. The study finds applications in understating the blood flow, modeled as non-Newtonian micropolar fluid, in stenotic arteries especially. The outcomes can be used in designing the biomedical devices and techniques for cardiovascular treatments, e.g., evaluating the thrombogenic potential of implantable cardiac devices [25]. The rest of the article is structured as follows. Section 2 explains the mathematical formulation of the problem and method. Section 3 presents the results and discussion. Section 4 displays the conclusions.

Materials and Methods
A two-dimensional pulsatile flow of MP fluid was analyzed by a uniform magnetic field applied perpendicular to the flow direction, as shown in Figure 1. The geometry under consideration was a constricted channel. The center of the constriction was placed at x = 0 with a total width of constriction as 2x 0 , as depicted in Equation (1). The constrictions on the walls of the channel are formulated and implemented as: where y 1 and y 2 define the lower and upper walls with constriction heights h 1 and h 2 , respectively.

Materials and Methods
A two-dimensional pulsatile flow of MP fluid was analyzed by a uniform magneti field applied perpendicular to the flow direction, as shown in Figure 1. The geometr under consideration was a constricted channel. The center of the constriction was place at = 0 with a total width of constriction as 2 , as depicted in Equation (1). The con strictions on the walls of the channel are formulated and implemented as: where and define the lower and upper walls with constriction heights ℎ and ℎ respectively. The momentum equations of the unsteady flow are given by: The continuity equation is given by: Here the velocity components along the -and -axis are and , respectively. , and represent pressure, density, and kinematic viscosity, respectively. represent the micro-rotation velocity, represents vortex viscosity, ≡ , , the curren density, ≡ (0, , 0) the magnetic field with uniform strength , the electric con ductivity, and the dynamic viscosity. = ( + ) 2 ⁄ represents the spin gradien viscosity, where defines the micro-inertia density. If ≡ , , indicates th electric field directed along the normal to the flow plane, then ≡ (0, 0, ). In addition using Ohm's law: = 0, = 0, = ( + ) (6 Maxwell's equation ∇ × = 0 for stationary flow implies that = , where is constant, assumed to be zero for simplicity. Then, = . Therefore, applying × = − , Equation (2) becomes: The momentum equations of the unsteady flow are given by: The continuity equation is given by: Here the velocity components along thex-andŷ-axis areû andv, respectively.p, ρ, and ν represent pressure, density, and kinematic viscosity, respectively.N represents the micro-rotation velocity, k represents vortex viscosity, J ≡ J x , J y , J z the current density, B ≡ (0, B 0 , 0) the magnetic field with uniform strength B 0 , σ the electric conductivity, and µ the dynamic viscosity. γ = j(µ + k)/2 represents the spin gradient viscosity, where j defines the micro-inertia density. If E ≡ E x , E y , E z indicates the electric field directed along the normal to the flow plane, then E ≡ (0, 0, E z ). In addition, using Ohm's law: Maxwell's equation ∇ × E = 0 for stationary flow implies that E z = a, where a is a constant, assumed to be zero for simplicity. Then, J z = σûB 0 . Therefore, applying J × B = −σûB 2 0 , Equation (2) becomes: We define the dimensionless quantities: Here, T is the period of flow pulsation, L the maximum channel width, Re the Reynolds number, M the Hartmann Number, the micro-rotation velocity, and K the MP parameter.

Vorticity-Stream Function Formulation
The dimensionless stream function (ψ) and vorticity function (ω) for the flow under consideration are as follows: Some manipulations with Equations (9) and (10) produce: Using the quantities in Equation (13), we obtained the following vorticity transport equation as: Again, using the quantities in Equation (13), Equation (11) becomes: Here, u, v, and N are primitive variables, and ω and ψ are non-primitive variables.

Boundary Conditions
The steady case of the fluid flow from Equation (7) is considered to obtain the boundary conditions for the current problem: where J × B = −σ(E z +ûB 0 )B 0 . By substituting in Equation (18) and rearranging: Using the dimensionless variables from Equation (8) and some manipulations results in the following: Here C = (1 + K). Approximating the term on the right-hand side of Equation (20): Solving Equation (20) gives: The inlet velocity profile for M = 0 is: where u(y) represents the steady velocity profile given by Equations (22) and (23). A sinusoidal time-dependent flow is considered for pulsatile flow: Further, u = 0 and v = 0 (i.e., no-slip conditions) are considered on the walls. The proper boundary conditions for N on both the walls are: where 0 ≤ s ≤ 1. s = 0, s = 1/2, and s = 1 are for the flow with high concentration, weak concentration, and turbulence, respectively. N = 0 is considered for the inlet boundary condition of the micro-rotation velocity function. The outlet boundary conditions are set considering the flow fully developed.

Coordinates Transformation
Consider the following relation for transforming the coordinates: For computation purposes, we mapped the constriction to a straight channel which resulted in mapping the domain [y 1 , y 2 ] to [0, 1]. Equations (15)-(17) on applying Equation (26) result as follows: where: The velocity components u and v becomes: The boundary conditions at the walls, in the (ξ, η) coordinate system for ψ, ω, and N are: The value of determines the nature of the flow, where 0 and 1 represent the steady and pulsatile flows, respectively.

Numerical Method
The finite difference method was employed to acquire the numerical solution of Equations (27)-(29) over a uniform structured Cartesian grid ξ i , η j . The solution at time level l + 1 = l + ∆t, for l = 0, 1, 2, · · · , was computed using the known solution at time level l. To obtain the solution at the time level l + 1, firstly, the space derivatives of Equation (29) were discretized using the central difference, and the resulting linear system was solved for ψ = ψ(ξ, η) by the tri-diagonal matrix algorithm (TDMA) method. Then, Equations (27) and (28) were solved for the vorticity function ω = ω(ξ, η) and micro-rotation function N = N(ξ, η) by the alternating direction implicit (ADI) method. The over-relaxation parameter used for computations was λ = 1.4. The execution time of the calculations could be reduced by parallel implementation of the computer program, Ali and Syed [26]. However, developing a parallel solution on any shared, distributed, or hybrid memory programming paradigms is not a trivial task.
Equation (29) is discretized for the solution at advanced time level l + 1, and for l = 0, 1, 2, · · · , is given by: For the sake of simplicity, the l + 1 superscript from ψ is removed. Rearranging, Equation (33) results as: where A(j), B(j), C(j), and S(j) are given as: The solution of Equation (27) was computed at l + 1 2 time level by incorporating the solution of level l in the ADI method's first half. The explicit and implicit schemes at time levels l and l + 1 2 in ξ-direction and η-direction, respectively, were used while discretizing the derivatives of ω.
Equation (36) can be rearranged as: where A 1 (j), B 1 (j), C 1 (j), and S 1 (j) are given as: The ω at both the walls is given as: In the second step of the ADI method, using the solution computed at l + 1/2 level, the solution was obtained at the l + 1 time level. The explicit and implicit schemes at time levels l + 1/2 and l + 1 in the η-direction and ξ-direction, respectively, were used while discretizing the derivatives of ω.
Equation (38) can be written as: where A 2 (j), B 2 (j), C 2 (j), and S 2 (j) are given as: In a similar way, using the ADI method, the solution of Equation (28) was computed.

Results and Discussion
A grid of 400 × 50 was found to be suitable for the current work after a grid independence test was carried out for multiple grids with −10 ≤ ξ ≤ 10 and 0 ≤ η ≤ 1. For ξ and η directions, we considered the step length 0.05 and 0.02, respectively. The constriction length, i.e., x 0 , was 2. The height of the constriction on both walls was considered as 0.35. The time step, ∆t, was taken as 0.0001. We considered t = 0, 0.25, 0.50, 0.75 to show the influence of the flow controlling parameters in a pulsatile cycle. These four time levels were corresponding to the specific states of the flow pulsation: t = 0 corresponded to the start of pulsation motion, t = 0.25 corresponded to the maximum flow rate, t = 0.50 corresponded to the minimum flow rate, and t = 0.75 corresponded to the instantaneous zero flow rate. The magnitude of the WSS was the same for the upper and lower walls. Therefore, the WSS distribution was depicted only on the upper wall in the study. Figure 2 presented the axial velocity (u) profile and micro-rotation velocity (N) profile for the four pulsation cycles at different values of η and at the center of constriction (x = 0) with M = 5, St = 0.02, K = 0.6, and Re = 700. The phase-amplitude of u profile increased, whereas a decrease in the shifting phase was observed as the distance from the bottom wall increased. An opposite behavior for N was observed. For validity of the present scheme, Figure 3 compares the present study with Bandyopadhyay and Layek [22] for the WSS by varying magnetic field strength. The results were found to be promising in the comparison.  Figure 2 presented the axial velocity ( ) profile and micro-rotation velocity ( ) profile for the four pulsation cycles at different values of and at the center of constriction ( = 0) with = 5, = 0.02, = 0.6, and = 700. The phase-amplitude of profile increased, whereas a decrease in the shifting phase was observed as the distance from the bottom wall increased. An opposite behavior for was observed. For validity of the present scheme, Figure 3 compares the present study with Bandyopadhyay and Layek [22] for the WSS by varying magnetic field strength. The results were found to be promising in the comparison.  The WSS on the upper wall in a pulsatile cycle at different times for = 0.02, 0.04, 0.06, and 0.08 with = 5, = 0.6, and = 700 is shown in Figure 5. During 0 <   The WSS at the four time levels for K = 0.3, 0.6, 0.9, and 1.2 with St = 0.02, M = 5, and Re = 700 is shown in Figure 6. The WSS fell with the increasing values of K during a complete cycle. During 0 < t < 0.25, the flow separation region had an inverse relation with K. At t = 0, a decrease in the WSS was witnessed with increasing K. At t = 0.25, the WSS reached the maximum peak for all values of K. During 0.25 < t < 0.75, the flow started to decelerate, and a decrease in the WSS for all the K's was observed. The flow separation region slightly expanded with the increasing values of K. The WSS altered its sign at t = 0.75.
The WSS at the four time levels for Re = 500, 700, 900, and 1100 with St = 0.02, M = 5, and K = 0.6 is shown in Figure 7. The WSS has an inciting trend towards Re. At t = 0, an increase in the WSS was witnessed with increasing Re. At t = 0.25, the WSS reached the maximum peak for all values of Re. The flow started to decelerate during 0.25 < t < 0.75 and a decrease in the WSS for all the Re's was observed. The flow separation region expanded with the increasing values of Re. The WSS altered its sign at t = 0.75.   during a complete cycle. During 0 < < 0.25, the flow separation region had an inverse relation with . At = 0, a decrease in the WSS was witnessed with increasing . At = 0.25, the WSS reached the maximum peak for all values of . During 0.25 < < 0.75, the flow started to decelerate, and a decrease in the WSS for all the 's was observed. The flow separation region slightly expanded with the increasing values of . The WSS altered its sign at = 0.75.        Figure 9b,c, it can be seen that near the lower wall, the micro-rotation velocity boundary layer had a declining behavior towards Re and M. In contrast, it had inclining behavior towards Re and M near the upper wall.   Figure 10 presents the influence of M on the streamline at the four time levels. The flow separation region had a direct relation with M. At t = 0, the streamlines ran smoothly over the constriction. It can be seen that at t = 0.25 and t = 0.50, the flow separation region was decreased for larger values of M. At t = 0.75, vortices took up the largest portion of the channel. A symmetric behavior can be seen in the streamlines. Figure 11 presents the effects of varying the Reynolds number on the streamline at the four time levels for Re = 500, 700, 900, and 1100 with M = 5, K = 0.6, and St = 0.02. The streamlines near the constriction were smooth. It is noted as well that the disturbance in the flow tended to grow, leading towards the turbulence, as the value of Re was increased. The flow separation region was maximum at t = 0.25 for Re = 1100, as can be seen in Figure 11d.

Conclusions
The pulsatile flow, in a constricted channel, of non-Newtonian MP fluid under the impact of the applied magnetic field was examined numerically on a Cartesian grid. The effects of M, St, Re, and K on the WSS, u, and N profiles were studied. The key outcomes of the present study are listed as follows:

•
The direct relation between the WSS and M was observed, and the WSS attained its peak value at t = 0.25. The flow began to decelerate as the flow rate tended to be zero. The sign of the WSS changed when the net flow rate was zero;