Non-isothermal Bingham fluid flow between two horizontal parallel plates with Ion-slip and Hall currents

This study presents the numerical solution of velocity and temperature fields based on mass conservation, momentum and energy balances for the time-dependent Couette-Poiseuille flow of Bingham materials through channels. The channel flow of Bingham fluid concerns the flow of cement paste in the building industry and the mudflow in the drilling industry. The specific aim is to introduce the magnetohydrodynamic (MHD) phenomena specified by both Ion-slip and Hall currents into the non-isothermal channel flow in a theoretical approach. The Bingham constitutive equation is formulated by the generalized Newtonian fluid technique and solved by employing the explicit Finite Difference Method (FDM) using the MATLAB R2015a and Compaq Visual FORTRAN 6.6a both. For the exactness of numerical performance estimations, the criteria for stabilization and the convergence factor are analyzed. The velocity and temperature profiles are discussed individually at the moving and stationary walls of the channel. It is observed that magnetohydrodynamic phenomena accelerate the flow, and the temperature distributions reach the steady-state situation earlier than velocity distributions. Furthermore, the dominance of MHD parameters on the velocity distributions, shear stress, temperature distributions, and Nusselt number are discussed.


Introduction
Bingham material is a viscoplastic fluid retaining a yield strength. Toothpaste is a well-known example that won't be fed until a certain pressure is employed in the tube. Enormous samples of Bingham fluid is observed in many geological and industrial materials such as handling of slurries, concrete paste, mudflow, lava, ceramic paste, painting oil, chocolate, etc. On the other hand, the MHD (magnetohydrodynamic) flow nature has applications in multiple flow-governing devices like-electrostatic precipitation, MHD power generators, accelerators, aerodynamic heating, etc., in the polymer connected technology and the petroleum industry.
The Bingham fluid was first introduced by Bingham [1], which shows a non-Newtonian behavior with a biviscous rheological model presented in [2]. Several theoretical studies were found in the literature discussing the Bingham fluid flow in different process conditions that address the channel flow. In this regard, the axisymmetric flow feature of Bingham fluid with the consideration of suddenly imposing an unchangeable pressure gradient was analyzed here [3]. On the basis of Couette-Poiseuille flow under slip conditions, a study on Bingham fluid was carried out by Chen and Zhu [4]. It was found that with the coexistence of double shear flow and plug flow, the increment in the slip parameter reveals an opposite nature for the shear rate and plug flow criteria. Analysis of the simultaneous approach of normal-mode of classical type and the approach of a non-modal kind on Bingham-Poiseuille flow was presented in [5]. Analytically, Eight different forms of velocity profiles were gained [6] by observing axially-orienting Couette-Poiseuille flow of Bingham fluid. By developing a single equation in pressure gradient explicitly for the flow of Bingham plastic through an annulus, an analytical solution was discussed, see [7]. A step was taken to inquire about a non-modal approach to Bingham plastic flow [8]. It was reported that the disturbance of non-axisymmetric phenomena has the critical energy Reynolds number in the lowest specification. The flow of dissipative fluid with the properties under temperature dependency was investigated [9]. As expected, this study observed that the thermal buoyancy parameter causes a hike in velocity as well as temperature. In [10], the start-up flow of Bingham plastic was presented by performing verification with small yield stress values. The analysis of the channel flow of Bingham rheological fluid by introducing a relaxation parameter with the aid of Lattice Boltzmann Method has performed by Ortega [11]. A semi-analytical solution based study of Bingham rheology with linearly varying yield stress along with plastic viscosity was presented in [12]. The use of a concentric annulus was conducted, see [13], to investigate velocity profiles of the analytical type of the pressure-driven flow of Bingham fluid. Some studies are recently performed where Bingham fluid was considered in different cases of the flow systems [14][15][16].
In order to manipulate the fluid flow in some real-life problems, the MHD effect is required to introduce. Relating to the MHD effect on Bingham fluid flow, the impact by raising the Hall and Ion-slip currents were theoretically discussed in several recent studies based on different process conditions of channel flow [17][18][19]. The Hall impact on MHD Couette flow of Bingham fluid along with the cases of suction, injection, Joule heating, viscous dissipation was analyzed in [20]. It was concluded that the velocity components fall when the yielding stress develops quantitatively, and viscosity influences the steady-state significantly. The MHD dissipative flow of Bingham fluid by concerning steady feature over a rotating porous disk with Ion-slip and Hall currents in a rotating system was considered in [21], along with that, the shear stress coefficient and the heat transfer coefficient were calculated. The Numerical simulation for Bingham fluids was pointed out in [22] while the Taylor Couette flows through two concentric cylinders are considered. Squeezed vortices were noticed by the influence of yield stress in the direction of the inner cylinder. The numerical determination of MHD flow through parallel plates for Bingham material with Hall current, viscous dissipation, and without suction was performed [23]. MHD Ciliary propulsion flow of microscopic organisms with Hall and Ion-slip effects was studied in [24]. It was mentioned that the development of the magnetic parameter tends down the bolus shaped by the flow.
The suction, injection, Joule heating, Hall current, Ionslip current, surface heat flux are the possible effects that could be considered in Bingham material flow in channels with different wall boundary conditions as well as different types of plates like the stationary plate, moving plate, porous plate, Riga plate, rotating, and oscillatory plate, etc. [25][26][27]. All the presented studies concerned one or more of the above effects with one of the mentioned channel conditions for Bingham material flow. However, the MHD effect on the Bingham fluid in channels is yet to understand in different process conditions.
The present research focuses on the MHD effect personalized by Hall and Ion-slip currents. Thus, considering all the above base studies, our specific purpose is to investigate the flow characteristics numerically for time-dependent MHD Bingham material flow between two horizontal parallel plates where the upper plate is assumed to have a constant velocity. The generalized Ohm's law [28] appears in the momentum balance, including both Hall and Ionslip currents. Also, the Couette flow is considered when a pressure gradient drives the flow. Furthermore, Energy dissipation and uniform magnetic field are involved. The explicit FDM (Finite Difference Method) framework was utilized to solve the non-isothermal problem equations. The obtained results are plotted and shown graphically. Section 2 describes the physical model and its mathematical illustration. Section 3 describes the representation of shear stress and Nusselt number based on the present problem. The numerical technique is concisely presented in Sect. 4. Section 5 presents the outcome of the stability analysis for this problem. The results are discussed in Sect. 6

Mathematical formulation of the problem
The transient MHD Bingham fluid is considered to be streaming between two infinitely elongated horizontal parallel plates planned at y = ±h planes also ranging from x = 0 to ∞ and from z = 0 to ∞ . The upper plate is set to be occupied a fixed velocity U 0 and lower one is planned at motionless condition, while two individual unchanged temperatures T 2 and T 1 respectively are specified for both the upper and lower plates under the conditioning of T 2 > T 1 . Such Couette-Poiseuille flow was studied in Ref. [20,21]. Along the X-axis, a fixed pressure gradient dp dx is employed on the fluid, also an external magnetic field B 0 of uniform nature is acted vertical to the X-axis. The generation of secondary velocity at the Z-axis direction is expectedly happening due to the Hall and Ion-slip currents. The velocity vector for the governing fluid is given as, = u + v + w . Figure 1 presents the schematic flow configuration together with the boundary conditions of the problem.
The generalized Ohm's law describes the current density , with the electric field , whose representation follows: where is the fluid's electric conductivity, e is the Hall parameter and i is the Ion-slip parameter. The simultaneous expressions of the Hall and Ion-slip currents were studied from Ref. [21]. The simplified form of the above equation is: Thus the momentum equation, for this study, takes the following vector form: Here, for Bingham fluid, the viscosity is not a constant, as it is a non-Newtonian fluid. The initial approach for this fluid was described in [1]. Many of the other scholars are investigated later. The apparent viscosity is given by: = K + 0 |̇| , where K is the plastic viscosity of Bingham materials that gives resistance to the flow, and 0 is the yield stress below which Bingham materials behave like a solid, the rheological equation for the yield stress was described in [2]. This is represented by, (see [20] and [23]): With the considerations described above, the study, for the unsteady two-dimensional (2D) non-newtonian Couette-Poiseuille flow of Bingham fluid considering Hall and Ion-slip currents, is governed by the mass conservation, momentum and energy balances [20,23]. Finally, the governing equations take the following form: Continuity equation: Momentum equations: Energy equation: where x, y are the Cartesian coordinates; u, v are the x, y components of flow velocity respectively; p is pressure; is the density of the fluid; is thermal conductivity and c p is specific heat at the constant pressure and the viscosity The schematic flow characteristics are subjected to the following initial and boundary estimations: The retrieved dimensionless coupled nonlinear PDEs are given below: and the corresponding dimensionless boundary conditions are shown below: The parameters used in non-dimensional operations of Eqs. (10-13) are given below:

Shear stress and Nusselt number
Some impressive physical quantities which are relevant to this study are presented here. The shear stress along with the Nusselt number, are investigated as a function of different parameters where both are calculated at the plates, upper and lower (moving and stationary). The mathematical explanations of shear stress and Nusselt number are taken from Ref. [19 and 25]. The shear stress along X-direc- The Nusselt number along X-direction for stationary

Numerical technique
Under the boundary settings of Eqs. (14 and 15), the explicit FDM was chosen to solve the obtained dimensionless nonlinear coupled PDEs in Eqs. (9)(10)(11)(12)(13). For which, the boundary layer regime was essential to divide into tiny mesh spaces. For this division, several lines, parallel to both the axes, were taken. A particular test (mesh sensitivity) was performed to choose these number of lines, see Sect. 6.1. In this study, the direction of X-axis is maintained along the plate while Y-axis is maintained perpendicular to the plate-direction, shown in Fig. 1. Along the X-axis, parallel to Y-axis, the appropriate number of lines were chosen m = 40 . Along the Y-axis, parallel to X-axis, the appropriate number of lines were chosen n = 40 . Thus it was arbitrarily decided that the size of each of the plates was, X max = 40 i.e. X varied from 0 to 40 and the distance between the plates was, Y max = 2 i.e. Y varying values from 0 to 2. Figure 2 represents the mesh spacing for the current study.
Th e u n i fo r m ce l l e d g e s a re ΔX a n d ΔY along X and Y paths respectively which are t a k e n a s ΔX = with the smaller (in magnitude) time-step, Δ = 0.0001.

Stability and convergence analysis
The criteria of stability and convergence were analyzed along with the FDM due to demonstrate the expected order of convergence and the stable numeric computations. In the convergence test, Fourier exponential expansions for U , W and are used. With the use of created uniform mesh sizes, the stability and convergence criteria are simplified as follows: Therefore, the converged solutions will be found after the adoption of the above two inequalities being perfect. Using ΔY = 0.05, Δ = 0.0001 and the initial conditions ( U = 0 and V = 0 ), the above inequalities give the restrictions for the dimensionless parameters. Finally, the constraints for the dimensionless parameters are obtained as: P r ≥ 0.08, i ≥ 2, H a ≤ 20 and R e ≥ 0.0166 with the constant value e = 0.10 and E c = 0.10.

Results and discussion
Due to exploring the inner physical situation of the flow model, the numerical values in steady-state situations are reckoned for the non-dimensional primary velocity U , secondary velocity W , and temperature within each plate's boundary layer region. First, establishing a compatible mesh zone is ascertained inside the boundary layer regime through the sensitivity test of meshes for the numerical computations, see Sect. 6.1. Secondly, a code validation test is performed by two different tools, MATLAB R2015a and COMPAQ VISUAL FORTRAN 6.6a (CVF 6.6a), to verify the accuracy of the coding, see Sect. 6.2. Thirdly, the sensitivity test of time is ascertained to assure the steadystate solutions for all profiles, see Sect. 6.3. Fourthly, to test the qualitative track of the present study, a comparison between the current result and the published result is described, see Sect. 6.4. Finally, the effect of magnetohydrodynamic parameters, Ion-slip parameter i and Hartmann number H a on velocity profiles-U and W , temperature field also on shear stress L and Nusselt number Nu L are discussed at both the stationary plate (lower plate) and moving plate (upper plate), see Sect. 6.5.

Mesh sensitivity test
To find out the precise grid space lines (i.e. values of m and n ), several trial values of mesh spaces are experimented, out of these trial values three separate grid spacing as m = 40, n = 40; m = 100, n = 100 and m= 150, n= 150 are shown in Fig. 3 for the primary velocity field. The obtained curves are smooth and have a minor change. It is also checked for a quoted point of Fig. 3 by enlarging that part. The sub-figure (enlargement) of Fig. 3 shows that the doted, dash-dot and dashes lines are closed to each other that clarify the mesh independency of the model. Thus, m= 40 and n= 40 can be chosen as the appropriate mesh space lines. Therefore, further computations is continued for the mesh space lines, m= 40 and n= 40.

Code validation test
The exactitude of numerical codes is simulated for propriety, by MATLAB R2015a and COMPAQ VISUAL FORTRAN (CVF) 6.6a tools, with explicit FDM technique. Using MAT-LAB code, the obtained primary velocity field and shear stress are shown in Fig. 4-a and c respectively for several values of Ion-slip parameter i . Again, by using FORTRAN code, the obtained results of the same profiles, described above, are shown in Fig. 4-b and d respectively for the same values of Ion-slip parameter i . The rigorously identical results are discerned for both the calculating tools MATLAB and CVF (see by relating Fig. 4(a) with 4(c), as well as Fig. 4(b) with 4(d)). Thus, the sensitivity of coding achieved accuracy. However, in all the following sections, it is decided to proceed with the results from MATLAB as a simulation tool.

Time sensitivity test (steady-state solutions)
To meet the steady-state solutions of the configured mathematical flow scheme developed here, the numeric computations for both the velocities (primary, secondary) and temperature profiles are sustained for multifarious sequential dimensionless times, such as = 0.20, 0.70, 1.20, 2.00, 2.50 and 3.00. Figure 5-a, b and c are the representative graphs for the computed results of dimensionless primary velocity U, dimensionless secondary velocity W and dimensionless temperature show a little change after τ = 2.00 and also show a negligible change between τ = 2.50 and τ = 3.00. Depending on this numerical experimentation, the time value τ = 3.00 is specified as the ultimate steady-state for all variables in this study. Figure 5-a, b and c show that both velocity distributions and temperature distributions attain steady-state monotonically. Note that the temperature distribution is reached steady-state position faster than the velocity distributions.

Validity test (comparison)
The quantitative and qualitative measurement of the upshots of current study is analyzed with the published results of [23] through Figs. 6 and 7. It is shown from the Figs. 6 and 7, the reported results are not in quantitative similarity rather they are in an excellent qualitative similarity with the published results of [23]. Since the present results is verified by both CVF 6.6a and MATLAB R2015a tools, thus it can be claimed that the reported model is less blemished with the concern of published results.

Effect of magnetohydrodynamic parameters
To quest and clarify the physical modulation of the flow structure, the influences of magnetohydrodynamic parameters-Ion-slip parameter i and Hartmann number H a on velocity fields and temperature fields, shear stress and Nusselt number at the moving as well as stationary plates are illustrated via graphs in Figs. 8,9,10,11,12,13,14,15,16,17. With a general use of Hall parameter e = 0.10 , Reynolds number R e = 3.0 , Bingham number D = 0.1 , Prandtl number P r = 0.08 , and Eckert number E c = 0.10 that satisfies the stability condition of the present model discussed in Sect. 5. The magnitude of the dimensionless parameters that have been used, in this study, are chosen as investigated in some recent studies by Jha et al. [29] and Ibrahim and Anbessa [30]. For ensuring a concise discussion, the impacts of rest of the parameters namely, Hall parameter e , Prandtl number P r , Eckert number E c , Reynolds number R e and Bingham number D are not discussed. To observe the Couette-Poiseuille flow phenomena more clearly and to observe the flow properties adjacent to the plates, the Figs. 8,9,10,11,12,13,14,15,16,17 are presented in a different way, where the flow is chosen to discuss at both the upper and lower plates, in X direction. Thus, for both the plates, i.e. the lower stationary plate and the upper moving plate, the varying feature of the velocity distributions and temperature distributions are considered in the X direction.

Effects on primary velocity
The impacts of Hartmann number H a and Ion-slip parameter i on the profiles primary velocity are shown in Figs. 8 and 9. Figure 8(a, b), reveals that, H a opposes the primary velocity subject to the growth of H a at both plates. It is happened quantitatively as expectations because the orthogonally applied magnetic field to the flow direction generates a force of resistance naming Lorentz force. It is found from Fig. 9(a, b) that i influences the primary velocity at both the plates. The inclusion of the Ion-slip parameter enhances the operative conductivity for which the primary velocity grows up. Viewing the primary velocities for both the plates, it is noted that the primary velocity characterized at moving plate is speedier than the stationary one.

Effects on secondary velocity
The varying estimations of secondary velocity with Hartmann number H a and Ion-slip parameter i are shown in Figs. 10 and 11. From Figs. 10(a, b) and 11(a, b), it is clear that H a and i both oppose the secondary velocity subject to their growth at both plates. By comparing the secondary velocity profiles among both plates, it is noted that the secondary velocity characterized at moving plate is speedier than the secondary velocity at stationary one.

Effects on temperature
The impacts of Hartmann number H a and Ion-slip parameter i on temperature are graphed in Figs. 12 and 13. Figure 12(a, b), focused that, the temperature profile falls when the values of H a develops at both the plates. It is marked from Fig. 13(a, b) that the temperature influences by i at stationary plate while i occurs diminution on temperature by the rising values of i at moving plate. By comparing temperatures among both plates, it is found that the moving plate's temperature is greater than the stationary plate's temperature. It is also verified by the close-up view of the quoted point of the curves of Figs. 12 and 13. The enlargement sub-figures of Figs. 12(a, b) and 13(a, b) show that the effects are clear.

Effects on shear stress
The impacts of Hartmann number H a and Ion-slip parameter i on shear stress are shown in Figs. 14 and 15. From Fig. 14(a, b), it is found that H a opposes the shear stress subject to the growth of H a at both the plates. Figure 15(a, b) marks that, i influences the shear stress estimations at stationary plate as well as moving plate. By comparing shear stresses among both plates, it is noted that the shear stress characterized at the plate of motion is larger than the shear stress at the plate of stationary position.

Effects on Nusselt number
The impacts of Hartmann number H a and Ion-slip parameter i on Nusselt number are shown in Figs. 16 and 17. Figure 16(a, b) marks that, the Nusselt number influences with the growth of H a at both the plates. By observing, Fig. 17(a, b), it is noted that, i opposes the Nusselt number (heat transfer rate) at both the plates. By comparing Nusselt numbers among both plates, it is checked that the Nusselt number characterized at the plate of motion (upper plate) is greater than the Nusselt number at the plate of stationary position. Also, it is verified by enlarging the quoted point of the curves of Fig. 17(a, b). The enlargement sub-figures of Fig. 17(a, b) show that the effects are clear.

Conclusion
In this work, the unsteady viscous incompressible flow through the inner zone between two horizontally placed parallel plates incorporating the magnetohydrodynamic impact by both the Hall and Ion-slip currents was investigated numerically by using explicit FDM technique. The flow is involving the industrial-viscoplastic/Bingham fluids between two infinitely extended horizontal plates. For the exactitude of the estimations of numerical performance, the criteria for stabilisation and convergence were analyzed. The obtained converging ranges for respective materials are P r ≥ 0.08,H a ≤ 20 , i ≥ 2, and R e ≥ 0.0166 with the constant value e = 0.10 and E c = 0.10 . The simulated outcomes are delineated graphically for several values of major parameters as Ion-slip parameter i and Hartmann number H a . For conciseness, the impacts of rest of the parameters, namely, the effects of Hall parameter e , Reynolds number R e , Prandtl number P r , Eckert number E c and Bingham number D are not discussed. Some significant outcomes of the study are given to the following section: In future work, the presented numerical model of the Bingham material could be used to predict the magnetohydrodynamic flow phenomena in different process conditions, for example, a rotating channel is of great interest to the authors.

Compliance with ethical standards
Conflict of interest The authors declare no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.