Impact of Lorentz force on the pulsatile flow of a non-Newtonian Casson fluid in a constricted channel using Darcy’s law: a numerical study

The present paper examines the flow behavior and separation region of a non-Newtonian electrically conducting Casson fluid through a two-dimensional porous channel by using Darcy’s law for the steady and pulsatile flows. The vorticity-stream function approach is employed for the numerical solution of the flow equations. The effects of various emerging parameters on wall shear stress and stream-wise velocity are displayed through graphs and discussed in detail. It is noticed the increasing values of the magnetic field parameter (Hartman number) cause vanishing of the flow separation region and flattening of the stream-wise velocity component. The study also reveals that the non-Newtonian character of Casson fluid bears the potential of controlling the flow separation region in both steady and pulsating flow conditions.


Mathematical formulation
We consider a non-Newtonian electrically conducting Casson fluid through a two-dimensional porous channel. The channel has constrictions on the upper and lower walls, which are placed L units apart from each other. The flow is subjected to a uniform magnetic field B perpendicular to the channel walls and an electric field J normal to the plane of flow. We consider a Cartesian coordinate system ∼ x, ∼ y such that the direction of the flow is along the ∼ x-axis and the direction of β is along the ∼ y -axis. The constrictions are spanned from x = −x 0 to x = x 0 with its center at x = 0, as shown in Fig. 1. It is assumed that the magnetic Reynolds number for the flow is very small, i.e., Re m ≪ 1.
The governing equations of the non-Newtonian fluids are highly non-linear and more complicated as compared to those of Newtonian fluids. Due to complexity, no single constitutive equation exhibiting all properties of such fluids is available. The Casson fluid model is often considered to be better than the general viscoplastic model in fitting the rheological data that is closely related to the blood flow. The rheological equation of state for an incompressible flow of Casson fluid model is as follows 15,33,34 : www.nature.com/scientificreports/ where π is the product of the component of deformation rate with itself, π = e ij e ij with e ij being the (i, j)th component of the deformation rate, π c is a critical value of this product based on the non-Newtonian model, µ B is the plastic dynamic viscosity of the non-Newtonian fluid, and p y is the yield stress of the fluid. When π < π c then Eq. (1)   x -and ∼ y -axes, respectively, ∼ p is the pressure, ρ is the density, ν is the kinematic viscosity, β is the Casson fluid parameter, J ≡ J x , J y , J z is current velocity, B ≡ (0, B o , 0) is the magnetic field, B o is the strength of the uniform magnetic field, σ is electric conductivity, µ m is the magnetic permeability of the medium. Assume that E ≡ E x , E y , E z = (0, 0, E z ) denotes the electric field. From Ohm's law In the case of the steady flow, Maxwell's equation ∇ × E = 0 implies that E z is constant. For the present study, we assume that E z is zero. Then, Eq. (4) gives, To obtain the dimensionless form of the system of Eqs. (7), (4), and (5), the following quantities are introduced: Here T is the period of flow pulsation, Re is the Reynolds number of the flow, St is the Strouhal number, M is the Hartman number of the flow, U is the characteristic flow velocity.
Using the quantities from Eq. (8) in Eqs. (7), (4), and (5) ∂u ∂x + ∂v ∂y = 0  Transformation of coordinates. The transformation for the coordinates is defined as: Due to the transformation (24), η = 0 and η = 1 correspond to the lower and upper walls, respectively. Eqs. (22) and (23) are reduced to the following forms: The boundary conditions for ψ and ω at the walls become: Here ǫ is the pulsating amplitude. For the steady flow, ǫ = 0 , and for the pulsatile flow ǫ = 1.

Results and discussion
The numerical solutions of Eqs. (25) and (26) subjected to the boundary conditions (Eq. 29) are computed using the finite difference method. The numerical method follows a standard approach, as used by Bandyopadhyay and Layek 32 . The computational domain is taken as The grid points ξ i , η j , for i = 1, 2, . . . , n and j = 1, 2, . . . , m , are obtained using the step sizes of �ξ and �η in ξ-direction and η-direction, respectively. A fixed step-size t is used to advance the solution in time. The solution at time level l is known whereas the solution at each time level l + 1 , for l = 0, 1, 2, . . . , is computed. For a time level l + 1 , firstly, Eq. (26) is solved for ψ = ψ(ξ , η) using the central differences for the discretization of the space derivatives and the well-known Tri-Diagonal Matrix Algorithm (TDMA) for the linear system. Secondly, Eq. (25) is solved for ω = ω(ξ , η) using the well-known Alternating Direction Implicit (ADI) method. For this, the time derivative is discretized using the forward/backward difference, and the space derivatives are discretized using the central differences. The resulting linear systems at each of the two half-steps of the ADI method are solved using TDMA. For ensuring the stability of the numerical scheme, the values of the step sizes are taken as follows: t = 0.00025 for the steady case, and t = 0.00005 for the pulsatile flow with �ξ = 0.05 , �η = 0.02 and Re = 700 . Although the results for the present work are computed in serial on a high-performance multicore machine, the solutions can be obtained through parallel computing to save the execution time 35 .
Distinct fluid parameters that describe the flow characteristics, e.g., Hartman number M , Strouhal number St , porosity parameter D a and Casson fluid parameter β are considered for obtaining the flow profiles. In the present study, we used the computational domain as, for −10 ≤ ξ ≤ 10 and 0 ≤ η ≤ 1 is discretized by a grid of 400 × 50 elements by setting step sizes of �ξ = 0.05 in ξ-direction and �η = 0.02 in η-direction. Further, we consider the flow Reynolds number Re = 700 , and the length of the constriction 2x 0 = 4. A symmetric constriction is considered on each of the upper and lower walls at the same location with heights h1 = h2 = 0.35 ; hence, the channel width remains 30% of the channel width at the throat of the constriction. The time step is considered as t = 0.00025 and 0.00005 for steady and pulsatile flows, respectively. As the pulsatile motion is modeled by adding the sinusoidal time-dependent function sin(2πt) in the inflow boundary condition, the effects of flow parameters are shown for time levels t = 0.0, 0.25, 0.50, 0.75 . These time levels at every cycle correspond to the start of the pulsating motion, the maximum flow rate, the minimum flow rate, and the instantaneous zero net flow, respectively.
The present numerical scheme is validated by comparison of the pulsatile flow results with the relevant ones obtained by Bandyopadhyay and Layek 32 , which discussed a Newtonian fluid without porosity effect. Figure 2 St   However, at the upstream of the constriction, the shear stress is linearly distributed on the wall. During the first quarter period of pulsatile time 0 ≤ t ≤ 0.25 , the flow is accelerating. Consequently, the peak value of the shear stress rises at the throat of the constriction and reaches its maximum at t = 0.25 . The effects of M during this phase of time are significant. As M is increased, the shear stress on the wall increases. Further, during this phase of time, the separation region reduces with increasing M . Thus, M has a substantial impact on the shear stress for the pulsatile flow as well. In the next half period of pulsatile time 0.25 < t ≤ 0.75 , the flow begins to decelerate. Consequently, the shear stress declines, whereas the flow separation region grows. The distribution of shear stress oscillates as well as alternates its sign at t = 0.75 . This behavior is not observed during the rest of the cycle. The results at the lower wall are symmetric about the zero wall stress line to those at the upper wall.
The effects of Strouhal number on the wall shear stress at the upper wall are shown in Fig. 5 for St = 0.02, 0.05, 0.08 with M = 5 , β = 5 , and D a = 0.002 . The distribution of shear stress increases with the increasing St at the start of the pulsatile motion t = 0, but the effects are not significant as compared to that of M where considerable variation is observed. The peak values of the shear stress coincide with the variation of St when the flow is maximum ( t = 0.25 ). In the decelerating phase particularly 0.25 < t < 0.5 , the two peak values of the wall shear stress with the opposite sign are observed one is on the middle of the constriction and other is near the downstream of the constriction, further, the peak values of the wall shear stress decreasing with increasing St , when the flow rate is minimum ( t = 0.5 ), but the situation is reversed and ultimately changes its sign at zero net flow ( t = 0.75 ). However, the region of separation decays with increments in St and grows with time. The results at the lower wall are symmetric about the zero wall stress line to those at the upper wall.
The influence of the Casson fluid parameter (β) on the wall shear stresses are shown in Fig. 6 for β = 10, 5, 1, 0.5 with M = 5 , St = 0.02 , and D a = 0.002 . The wall shear stress increases with β during a complete cycle. However, the wall shear stress distribution behaves differently not only in the different phases of a cycle but also for different values of β . The maximum peak values of the wall shear stress appear in the vicinity of the constriction (−x 0 < x < x 0 ) except for the time phase when the net flow rate is zero at t = 0.75 . At the start of the motion (t = 0) , an oscillating behavior of the wall shear stress is observed for large values of β equal to 10 and 5 in the upstream of the constriction (x ≤ −x 0 ) and tend to be straighter in the downstream of the constriction (x ≥ x 0 ) . In the acceleration phase (0 < t < 0.25) , the peak value of the wall shear stress increases,   Fig. 9, the streamlines run smoothly over the constriction at the start of the motion (t = 0) . However, with the start of flow acceleration, a small vortex (or bubble) appears near the walls in the lee of the constriction due to the boundary layer effects. The vortices get enlarged with time. In Fig. 10, the vortices on the walls are not the same in size and shape; a kind of asymmetric flow is found. The rotational behavior of streamlines vanishes as the value of M increases at time t = 0 . Furthermore, the rotation of flow in the lee of constriction grows in size with time and decay with increasing M . Thus, the effect of boundary on the fluid can be controlled by applying the external uniform magnetic field. Moreover, the vorticity contours are found to be anti-symmetric with respect to the central line during the whole of a period of the pulsatile flow.  Fig. 11, the streamline show smoothness over the constriction at the start of the motion (t = 0) . However, as the flow starts to accelerate, a small vortex (or bubble) appears near the walls in the lee of the constriction, and the size of vortices increases with time; a kind of similar behavior of streamlines are found as that is noted for different value of M .  • The peak value of the shear stresses on walls increases with time and reaches its maximum value when the incoming flow is maximum. At the start of motion, the boundary layer effects reduce with increasing M . Thus, the magnetic field parameter M significantly affects the shear stress on the walls and results in smoother flow even in the throat of constriction. However, the effects of the boundary layer continue to grow during the pulsatile flow. The boundary layer effects are higher at the minimum flow rate.

Data availability
The pictorial and graphical data used to support the findings of this study are included in the manuscript.