Unsteady Viscous Incompressible Bingham Fluid Flow through a Parallel Plate

Numerical investigation for unsteady, viscous, incompressible Bingham fluid flow through parallel plates is studied. The upper plate drifts with a constant uniform velocity and the lower plate is stationary. Both plates are studied at different fixed temperatures. To obtain the dimensionless equations, the governing equations for this study have been transformed by usual transformations. The obtained dimensionless equations are solved numerically using the explicit finite difference method (FDM). The studio developer Fortran (SDF) 6.6a and MATLAB R2015a are both used for numerical simulations. The stability criteria have been established and the system is converged for Prandtl number Pr≥0.08 with ΔY=0.05 and Δτ=0.0001 as constants. As a key outcome, the steady-state solutions have been occurred for the dimensionless time τ = 4.00. The influence of parameters on the flow phenomena and on shear stress, including Nusselt number, are explained graphically. Finally, qualitative and quantitative comparison are shown.


Introduction
A Bingham fluid is a non-Newtonian viscoplastic fluid that possesses a yield strength that must be outstripped before the fluid will flow. In many geological and industry materials, Bingham fluids are used as a general mathematical basis of mud flow in drilling engineering, including in the handling of slurries, granular suspensions, etc. Bingham fluid is named after Eugene C. Bingham, who declared its mathematical explanation.
In this regard, an exploration of the laws for plastic flow has been studied by Bingham [1]. Darby and Melson [2] formulated an empirical expression to guess the friction loss factor for the drift of a Bingham fluid. Bird et al. [3] analyzed the rheology and flow phenomena of viscoplastic materials. Convective heat transfer for Bingham plastic inside a circular pipe and the numerical approach for hydro-dynamically emerging flow and the simultaneously emerging flow were studied by Min et al. [4]. Liu and Mei [5] considered the slow spread of a Bingham fluid sheet on an inclined plane. For Bingham fluids, the Couette-Poiseuille flow between two porous plates, taking slip conditions into consideration, was investigated by Chen and Zhu [6]. Sreekala and Kesavareddy [7] mentioned the Hall impacts on magneto-hydro dynamics (MHD) Bingham plastic flow over a porous medium, including uniform suction and injection. For Bingham fluids, the MHD flow for an unsteady case considering Hall currents was described by Parvin et al. [8]. Rees and Bassom [9] considered Bingham fluids over a porous medium following a rapid modification of surface heat flux. Mollah et al. [10] numerically studied the ion-slip and Hall impacts for MHD Bingham fluid flowing through parallel plates, including the considerations of unsteady cases and suction. Gupta et al. [11] considered the MHD mixed convective stagnation point flow and the heat transfer of an incompressible nano fluid over an inclined stretching sheet. Mollah et al. [12] investigated Bingham fluid flow through an oscillatory porous plate with an ion-slip and hall current.
Along with these studies, our aim is to numerically investigate unsteady Bingham fluid flow between parallel plates. The current study is concerned with the generalized Ohm's law considering the ion-slip and Hall current is absent. Additionally, the Couette flow is considered where the viscous dissipation, pressure gradient, and rheology of Bingham fluids are involved. The explicit FDM has been used to demonstrate the dimensionless non-linear, coupled, partial differential equations. The achieved results are shown graphically.

Mathematical Formulation
The physical configuration of the considered problem that corresponds with the boundary conditions is shown in Figure 1. A laminar, incompressible, non-Newtonian Bingham fluid is considered to be flowing among two infinite plates situated at the y = ±h planes and lengthening from x = 0 to ∞. Mollah et al. [10] numerically studied the ion-slip and Hall impacts for MHD Bingham fluid flowing through parallel plates, including the considerations of unsteady cases and suction .Gupta et al. [11] considered the MHD mixed convective stagnation point flow and the heat transfer of an incompressible nano fluid over an inclined stretching sheet. Mollah et al. [12] investigated Bingham fluid flow through an oscillatory porous plate with an ion-slip and hall current.
Along with these studies, our aim is to numerically investigate unsteady Bingham fluid flow between parallel plates. The current study is concerned with the generalized Ohm's law considering the ion-slip and Hall current is absent. Additionally, the Couette flow is considered where the viscous dissipation, pressure gradient, and rheology of Bingham fluids are involved. The explicit FDM has been used to demonstrate the dimensionless non-linear, coupled, partial differential equations. The achieved results are shown graphically.

Mathematical Formulation
The physical configuration of the considered problem that corresponds with the boundary conditions is shown in Figure 1. A laminar, incompressible, non-Newtonian Bingham fluid is considered to be flowing among two infinite plates situated at the = ℎ planes and lengthening from = 0 ∞. The upper plate has a uniform velocity, and the lower plate is stationary. Two constant temperatures, and are considered respectively for both upper and lower plates, where > . A constant pressure gradient is driven on the fluid along -direction and the body force are neglected.
Within the above considerations, the model is governed by the following equations: Continuity equation: Momentum equation: Energy equation: where the apparent viscosity of the Bingham fluids: Where Κ represents the plastic viscosity of the Bingham fluid and represents the yield stress. The corresponding initial and boundary conditions can be given as follows: The upper plate has a uniform velocity, U 0 and the lower plate is stationary. Two constant temperatures, T 2 and T 1 are considered respectively for both upper and lower plates, where T 2 > T 1 . A constant pressure gradient dP dx is driven on the fluid along X-direction and the body force are neglected. Within the above considerations, the model is governed by the following equations: Momentum equation: Energy equation: where the apparent viscosity of the Bingham fluids: where K represents the plastic viscosity of the Bingham fluid and τ 0 represents the yield stress.
The corresponding initial and boundary conditions can be given as follows: For the numerical solution, the above Equations (1)-(5) are required to transfer dimensionless form. The dimensionless quantities (6) that have been used which are considered as follows: The non-dimensional variables (6) are used in Equations (1)-(5) to obtain the dimensionless governing Equations (7)- (11): where The non-dimensional parameters are as follows: Eckert number, The subjected dimensionless conditions can be noted as follows:

Shear Stress and Nusselt Number
The shear stresses at both the stationary and moving plates are studied from the velocity profile. The local shear stress in X-direction for stationary plate is τ L1 = µ ∂U ∂Y Y=−1 and for the moving plate . Moreover, the rate of heat transfer, or the Nusselt number, for the stationary and moving plates are studied from the temperature fields. The Nusselt number in the X-direction for the stationary plate is

Numerical Technique
A set of finite difference approaches are required to solve dimensionless coupled non-linear partial differential Equations (7)- (10). The dimensionless coupled non-linear partial differential Equations (7)-(10) have been solved by using explicit FDM, subject to the boundary conditions (11). There are number of experiments to solve this kind of unsteady time-dependent problem using FDM where the pressure gradient is constant. The stability and convergence criteria of the explicit FDM are established for finding the restriction of the values of various parameters to get converged solutions (Mollah et al. [10,12] and Akter et al. [13]).
Here, the region inside the boundary layer is distributed into grid spaces of lines parallel to the X and Y axes; where the X axis is towards the plate and Y axis is normal to the X axis. It is determined that the height of the plate is X max = (40), i.e., X extends from 0 to 40 and Y max = (2) which corresponds to Y → ∞ , i.e., Y limits from 0 to 2 (length). There are m = 40 and n = 40 grid spacing lines in the X and Y directions, respectively, as shown in Figure 2. It is shown that ∆X and ∆Y are constant mesh or grid space lines towards X and Y directions, respectively, and is noted as follows, ∆X = 1.0 (0 ≤ x ≤ 40) and, with the smaller time-step of ∆τ = 0.0001. Here, the region inside the boundary layer is distributed into grid spaces of lines parallel to the and axes; where the axis is towards the plate and axis is normal to the axis. It is determined that the height of the plate is = 40 , i.e., extends from 0 to 40 and = 2 which corresponds to → ∞, i.e., limits from 0 to 2 (length). There are = 40 = 40 grid spacing lines in the and directions, respectively, as shown in Figure 2. It is shown that ∆ and ∆ are constant mesh or grid space lines towards and directions, respectively, and is noted as follows, ∆ = 1.0 0 ≤ ≤ 40 and, with the smaller time-step of ∆ = 0.0001. denote the values of and at the last step of time, respectively. By applying the explicit FDM, the set of finite difference equations for the model are expressed as follows: Further, the finite difference subjected conditions may be described as follows: U and θ denote the values of U and θ at the last step of time, respectively. By applying the explicit FDM, the set of finite difference equations for the model are expressed as follows: Inventions 2019, 4, 51 5 of 11 Further, the finite difference subjected conditions may be described as follows:

Stability and Convergence Analysis
The stability and convergence criteria of the explicit FDM are established for finding the restriction of the values of various parameters to get converged solutions. Excluding the stability criteria and convergence analysis of the acquired finite difference method, the analysis will remain incomplete unless an explicit process is being used. For the constant grid spaces, the stability analyses of the design have been performed as follows.
Equations (12) and (15) will be disregarded, because ∆τ does not exist. At a time arbitrarily called t = 0, the Fourier expansion for U and θ at are all e iαX e iβY apart from a constant, where i = √ −1. The expressions for U and θ at time t = τ can be defined as follows: After the time-step, these terms have taken the form as follows: Substituting Equations (15) and (17) into Equations (13) and (14), concerning the coefficients U and V as constants over any one time-step, the obtained equations can be described as follows: The simplification of Equations (19) and (20) can be denoted in matrix form, and are expressed as follows: where, . The Eigen values of the above matrix T have been computed to acquire the stability condition and finally the stability conditions of the problem can be expressed in the following Equation (22): Using ∆Y = 0.05, ∆τ = 0.0001 and the initial conditions, the above equation gives P r ≥ 0.08.

Results and Discussion
In order to assess the physical characteristics of this developed mathematical model, the steady-state numerical solution have been simulated for the dimensionless primary velocity U and temperature fields θ within the boundary layer. The influences of Bingham number τ D , Reynolds number R e , Prandtl number P r and Eckert number E c on velocity and temperature fields, shear stress, and Nusselt number, where R e = 2.0, P r = 0.08, E c = 0.10 and τ D = 0.001 and following the procedures described by Akter et al. [13]. The total results were discussed as the following five parts: 1.
Examine the mesh sensitivity 2.
Examine the sensitivity of MATLAB and FORTRAN coding output.

3.
Finding the steady state solution 4.
Effect of parameters 5.
Comparison with the published results

Mesh Sensitivity:
To obtain the appropriate grid space for m and n the computations have been carried out for three different grids, (m = 20, 20); (m = 30, 30) and (m = 40, n = 40), which are shown in Figure 3. It is seen from this figure that the graph for (m = 40, n = 40) is more closed than others. Also it has been verified to enlarge the quoted point of the Figure 3a

Sensitivity of MATLAB R2015a and SDF 6.6a
The accuracy is determined for the simulated results by MATLAB R2015a and SDF 6.6a tools and explicit finite difference is used as a solution technique. Using MATLAB and Fortran code, the Figure 4a  The accuracy is determined for the simulated results by MATLAB R2015a and SDF 6.6a tools and explicit finite difference is used as a solution technique. Using MATLAB and Fortran code, the Figure 4a The accuracy is determined for the simulated results by MATLAB R2015a and SDF 6.6a tools and explicit finite difference is used as a solution technique. Using MATLAB and Fortran code, the Figure 4a,b are shown for velocity field and, Figure 4c,d are shown for the shear stress respectively. The impact of velocity field and shear stress at moving plate have been illustrated for several values of Reynolds number . The identical results are obtained for both the tools (see Figure 4a-d). Thus, the sensitivity of coding achieved accuracy.

Steady-State Solution:
To acquire the steady-state solution of this developed mathematical model, the computations for velocity and temperature have been continued for different dimensionless times, such as = 0.20, 0.70, 1.20, 3.00, 4.00, 5.00, 6.00, etc. It is shown in Figure 5a,b that the computations for have shown little change after = 3.00 and also shown negligible change after = 4.00. Thus, the solutions of all variables for = 4.00 are taken essentially as the steady-state solutions. Figure 5a,b show that the both velocity filed and temperature field attain steady-state monotonically. It also should be noted that the temperature profiles are reached steady state position than the velocity profiles.

Effects of Parameters:
To attain the clear conception of physical properties of the present study, the impact of parameters, namely Reynolds number and Bingham number on velocity and temperature, shear stress and Nusselt number at Prandtl number = 0.08 and Eckert number = 0.10 are illustrated in Figures 6 and 7. For brevity, the impact of Prandtl number and Eckert number are not shown.

Steady-State Solution:
To acquire the steady-state solution of this developed mathematical model, the computations for velocity and temperature have been continued for different dimensionless times, such as τ = 0.20, 0.70, 1.20, 3.00, 4.00, 5.00, 6.00, etc. It is shown in Figure 5a,b that the computations for U and θ have shown little change after τ = 3.00 and also shown negligible change after τ = 4.00. Thus, the solutions of all variables for τ = 4.00 are taken essentially as the steady-state solutions. Figure 5a,b show that the both velocity filed and temperature field attain steady-state monotonically. It also should be noted that the temperature profiles are reached steady state position than the velocity profiles.

Steady-State Solution:
To acquire the steady-state solution of this developed mathematical model, the computations for velocity and temperature have been continued for different dimensionless times, such as = 0.20, 0.70, 1.20, 3.00, 4.00, 5.00, 6.00, etc. It is shown in Figure 5a,b that the computations for have shown little change after = 3.00 and also shown negligible change after = 4.00. Thus, the solutions of all variables for = 4.00 are taken essentially as the steady-state solutions. Figure 5a,b show that the both velocity filed and temperature field attain steady-state monotonically. It also should be noted that the temperature profiles are reached steady state position than the velocity profiles.

Effects of Parameters:
To attain the clear conception of physical properties of the present study, the impact of parameters, namely Reynolds number and Bingham number on velocity and temperature, shear stress and Nusselt number at Prandtl number = 0.08 and Eckert number = 0.10

Effects of Parameters:
To attain the clear conception of physical properties of the present study, the impact of parameters, namely Reynolds number R e and Bingham number τ D on velocity and temperature, shear stress and Nusselt number at Prandtl number (P r = 0.08) and Eckert number (E c = 0.10) are illustrated in Figures 6 and 7. For brevity, the impact of Prandtl number P r and Eckert number E c are not shown.
both the shear stress and the Nusselt number at moving plate reduce with the increment of Reynolds number . The effects of Bingham number on velocity field, temperature field, shear stress and Nusselt number at moving plate are shown in Figure 7a  The influence of Reynolds number R e on velocity field, temperature field, shear stress and Nusselt number at moving plate are shown in Figure 6a-d. It is seen from Figure 6a,b, the velocity and temperature distributions decrease with the increase of Reynolds number R e . Also it is seen from Figure 6b that the temperature distribution has minor decreasing effect. It has been checked to observe the effects by enlargement of the marked point on the curves of Figure 6b. The enlarge sub-figure of Figure 6b shows clearly the minor decreasing effects. It is observed from Figure 6c

Comparison
Finally, qualitative and quantitative comparison of the study subject to the published results of Parvinet al. [8] are discussed in Figures 8 and 9.

Comparison
Finally, qualitative and quantitative comparison of the study subject to the published results of Parvin et al. [8] are discussed in Figures 8 and 9.
Thus, in Figures 8 and 9, the obtained results are quantitatively not identical but quite similar qualitatively with the published results of Parvin et al. [8]. The problem has been solved numerically by explicit FDM. Since the present estimation has been verified by SDF 6.6a and MATLAB R2015a both software. Hence, our model is less flawed concerning published results. 0.08 and = 0.10 at time = 4.00 (steady state).

Comparison
Finally, qualitative and quantitative comparison of the study subject to the published results of Parvinet al. [8] are discussed in Figures 8 and 9.   Thus, in Figures 8 and 9, the obtained results are quantitatively not identical but quite similar qualitatively with the published results of Parvin et al. [8]. The problem has been solved numerically by explicit FDM. Since the present estimation has been verified by SDF 6.6a and MATLAB R2015a both software. Hence, our model is less flawed concerning published results.

Conclusions
In this study, the unsteady viscous incompressible Bingham fluid flow through a nonconducting parallel plate has been investigated numerically by explicit finite difference method. The