Second Law Analysis of Unsteady MHD Viscous Flow over a Horizontal Stretching Sheet Heated Non-Uniformly in the Presence of Ohmic Heating: Utilization of Gear-Generalized Differential Quadrature Method

In this article, the entropy generation characteristics of a laminar unsteady MHD boundary layer flow are analysed numerically for an incompressible, electrically conducting and dissipative fluid. The Ohmic heating and energy dissipation effects are added to the energy equation. The modelled dimensional transport equations are altered into dimensionless self-similar partial differential equations (PDEs) through suitable transformations. The reduced momentum and energy equations are then worked out numerically by employing a new hybrid method called the Gear-Generalized Differential Quadrature Method (GGDQM). The obtained numerical results are incorporated in the calculation of the Bejan number and dimensionless entropy generation. Quantities of physical interest, like velocity, temperature, shear stress and heat transfer rate, are illustrated graphically as well as in tabular form. Impacts of involved parameters are examined and discussed thoroughly in this investigation. Exact and GGDQM solutions are compared for special cases of initial unsteady flow and final steady state flow. Furthermore, a good harmony is observed between the results of GGDQM and those given previously by the Spectral Relaxation Method (SRM), Spectral Quasilinearization Method (SQLM) and Spectral Perturbation Method (SPM).


Introduction
Physically, entropy is an assessment of molecular chaos or its randomness. As a thermally dynamic system becomes more disordered, the locations of the molecules become more and more uncertain and therefore their positions become less predictable and the entropy increases. The boost in the disorder of a thermodynamic system is termed entropy generation/production. Entropy generation determines the level of irreversibilities that accumulate during a process. With the increasing rate of entropy creation in any heated system, the quality of energy is reduced, that is, it destroys useful points [49][50][51][52]. Moreover, the present numerical results are portrayed and discussed thoroughly via various graphical and tabular illustrations, in order to examine the influences of several emerging key parameters, such as Prandtl number Pr, magnetic parameter M, Eckert number Ec, as well as the reduced dimensionless time ξ and the temperature difference parameter Ω.
The innovation of the current study lies essentially in the use of a new hybrid numerical method called Gear-Generalized Differential Quadrature Method (GGDQM) for studying thermodynamically the present unsteady boundary layer flow problem in the presence of viscous dissipation and Ohmic heating. The robustness and efficiency of the numerical results given by GGDQM are also compared analytically and numerically by considering the existing published results and introducing the notion of CPU time.

Flow and Heat Transfer Analysis
As schematically portrayed in Figure 1, we consider a two-dimensional unsteady laminar forced convective flow of an incompressible, viscous, and electrically conducting fluid driven by a linearly stretching horizontal surface, in the presence of a vertical applied magnetic field of constant strength B o , in such a way that the induced magnetic field is neglected under the assumption of a small magnetic Reynolds number. Initially (i.e., t = 0), the studied fluid and the sheet surface y = 0 are stationary and have a constant temperature T ∞ . After this time (i.e., t > 0), the sheet surface y = 0 is stretched linearly along the positive x-direction with a velocity U w = U o x and heated non-uniformly by an imposed nonlinear thermal boundary condition of the form T w = T ∞ + T o x 2 , where U o and T o are two dimensional constants that characterize the present unsteady boundary layer flow. The innovation of the current study lies essentially in the use of a new hybrid numerical method called Gear-Generalized Differential Quadrature Method (GGDQM) for studying thermodynamically the present unsteady boundary layer flow problem in the presence of viscous dissipation and Ohmic heating. The robustness and efficiency of the numerical results given by GGDQM are also compared analytically and numerically by considering the existing published results and introducing the notion of CPU time.

Flow and Heat Transfer Analysis
As schematically portrayed in Figure 1, we consider a two-dimensional unsteady laminar forced convective flow of an incompressible, viscous, and electrically conducting fluid driven by a linearly stretching horizontal surface, in the presence of a vertical applied magnetic field of constant strength  In the presence of an external magnetic field, viscous dissipation and Joule heating effects, the continuity, momentum, and energy equations governing the present unsteady boundary layer flow are written as follows In the presence of an external magnetic field, viscous dissipation and Joule heating effects, the continuity, momentum, and energy equations governing the present unsteady boundary layer flow are written as follows ∂u ∂x where t is the dimensional time, (u, v) are the tangential and normal fluid velocity components, respectively, µ is the dynamic viscosity of the incompressible fluid, ρ is the density, σ is the electrical conductivity of the conducting fluid, B o is the uniform magnetic field applied vertically along the y-direction, T is the temperature of the fluid throughout the boundary layer, k is the thermal conductivity of the working fluid and ρC p is the fluid heat capacitance. By introducing the following dimensionless quantities Equations (1)- (5) reduce to ∂θ ∂ξ Here, the subscripts η, ηη and ηηη used above for f and θ designate the first, second and third partial derivatives, respectively, with respect to the variable η.
The continuity equation described above by Equation (1) is satisfied identically by introducing the stream function ψ(t, x, y), such that Note that in the case of a non-uniform wall temperature heating condition with linear spatial variation (i.e., T w = T ∞ + T o x), the dimensionless energy equation takes the following form ∂θ ∂ξ In addition, the non-dimensional physical parameters ξ, Pr, Ec and M appearing above in Equations (7) and (8) are defined mathematically as Here, τ is the dimensionless time and ν is the kinematic viscosity of fluid, where τ = U o t and ν = µ/ρ. For the present boundary layer flow problem, the skin friction coefficient C f x and the local Nusselt number Nu x are defined as In dimensionless form, these engineering quantities reduce to

Closed Form Solutions
Generally, for a non-conservative fluidic system, the viscous dissipation and Joule heating terms cannot be neglected in the energy equation (i.e., Ec = 0). Therefore, the closed-form analytical solutions for f (ξ, η) and θ(ξ, η) cannot be obtained for the present dynamical system. On the contrary, the exact solutions of the dimensionless stream function f (ξ, η) for the limiting cases of initial unsteady flow (i.e., ξ = 0) and final steady state flow (i.e., ξ = 1) can be found easily from Equation (7). Hence, by making use of the boundary conditions related to f (ξ, η), the special solutions for f (0, η) and f (1, η) are expressed formally as follows By virtue of these solution expressions and the boundary conditions of f (ξ, η) and θ(ξ, η), we obtain from Equation (8) the following results The validity of the numerical results presented in this paper is confirmed in the next section by comparing our findings with those developed analytically for the physical quantities f ηη (0, 0), f ηη (1, 0), θ ηη (0, 0) and θ ηη (1, 0), as shown in Equations (20) and (21).

Second Law Analysis
As is well known, the volumetric entropy production rate of an electrically conducting fluid flowing in the presence of an externally applied magnetic field is defined thermodynamically by .

S gen
Using the boundary layer approximations, Equation (22) reduce to As shown in Equation (24), the expression of . S gen depicts three bases of entropy production. These thermodynamic sources are the heat transfer, the viscous friction and the Ohmic heating, respectively.
The dimensionless form of Equation (24) is called entropy generation number Ns. This thermodynamic quantity is given by where .
S gen o represents a characteristic entropy generation of the studied system.
By employing the similarity transformations shown in Equation (6), we get Here, Ω denotes the temperature difference parameter, where Another interesting thermodynamic quantity called the Bejan number Be can be computed from the different entropic terms shown in Equation (26) as follows where After introducing the expressions of N h , N f and N m into Equation (29), we get From the definition of the Bejan number Be, it is obvious that Be is always comprised between 0 and 1. The zero value of Be implies that the combined contribution of fluid friction and magnetic field completely overrides the heat transfer effect, while the unit value of Be indicates that the heat transfer mechanism is the only cause of entropy creation. On the contrary, it is found that the heat transfer and the combined effects of magnetic field and viscous dissipation can make an equal entropic contribution in cases where Be = 0.5.
From the implementation point of view, the unsteady boundary layer flow problem under consideration can be further simplified by considering the following changes It is worth noting that the dimensionless space variable χ is introduced above instead of η for reducing the physical space domain from [0, ∞] to [0, 1], in which η ∞ represents the optimum value of the boundary layer thickness.

Solution Methodology
Due to the unsteadiness of the studied boundary layer flow problem and its nonlinear dynamical behaviour, the governing partial differential equations (PDEs) along with their associated boundary conditions (i.e., Equations (1)-(5)) are subjected to several necessary simplifications and suitable similarity transformations before being solved numerically by means of a powerful numerical tool. For this purpose, the resulting set of coupled nonlinear differential equations and boundary conditions (i.e., Equations (33)-(36)) is handled numerically using the Gear-Generalized Differential Quadrature Method (GGDQM), in order to reach a precision to the tenth decimal place as the standard of convergence (see Table 1).

Gear-Generalized Differential Quadrature Method (GGDQM)
For realizing a fine spatial discretization for the variable χ, it is more useful to use GDQM with the modified Gauss-Lobatto grid points χ i , which are given by Here, N is the total number of Gauss-Lobatto collocation points, where 1 ≤ i ≤ N. Accordingly, the functions F (ξ, x) and Θ (ξ, x) defined above in Equation (32) can be approximated at a collocation point χ = χ i by In addition, the weighting coefficients d ij (n) appearing in Equation (42) are given by Shu [45] as follows Here, n represents the order of differentiation with respect to the variable χ.
After substituting the discretized form of F(ξ, χ) and Θ(ξ, χ) with their partial derivatives into Equations (33)-(36), we get the following semi-discrete system where 0 < ξ < 1. The linear and nonlinear parts L F ξ , L Θ ξ , NL F ξ and NL Θ ξ arising from Equations (33) and (34) are given by It is worth pointing out that the solutions of the initial unsteady flow (i.e., ξ = 0) and final steady state flow (i.e., ξ = 1) can be found numerically by solving successively the following nonlinear algebraic systems where ξ = 1. In the above limiting cases, the linear and nonlinear parts L F 0 , L F 1 , L Θ 0 , L Θ 1 , NL F 0 , NL F 1 and NL Θ 0 , NL Θ 1 shown above are expressed by By utilizing the Newton-Raphson Method (NRM), the nonlinear algebraic systems (S 0 ) and (S 1 ) can be handled and then solved accurately, in order to find the numerical estimate values of the solutions {(F i (0), Θ i (0)) and (F i (1), Θ i (1)) / 1 ≤ i ≤ N}. In this unsteady boundary layer flow problem, the solutions {(F i (0), Θ i (0))/ 1 ≤ i ≤ N} of the initial unsteady flow (i.e., ξ = 0) corresponding to the algebraic system (S 0 ) are taken as the initial conditions for the problem under consideration. Therefore, for generating numerically the general solutions {(F i (ξ), Θ i (ξ)) / 1 ≤ i ≤ N and 0 < ξ < 1} with η ∞ = 15, the non-autonomous differential system S ξ along with the initial conditions {(F i (0), Θ i (0))/ 1 ≤ i ≤ N} is integrated temporarily using an auto-adaptative implicit algorithm based on the Gear Method (GM). Furthermore, to achieve an absolute accuracy of the order of 10 −10 , it is found that the dimensionless time-step size ∆ξ and the number of collocation points N must be selected as ∆ξ = 10 −5 and N = 150 in all subsequent analyses.
Under the above convergence criterion, the dimensionless physical quantities of interest Re 1/2 x C f x , Re −1/2 x Nu x , Ns and Be can then be computed numerically from the solutions {(F i (ξ), Θ i (ξ)) / 1 ≤ i ≤ N and 0 ≤ ξ ≤ 1} as follows Also, the discretized forms of the non-dimensional velocity f η (ξ, η) and temperature θ(ξ, η) can be deduced as follows where

Validation of the Numerical Results
To authenticate the exactness of our solution methodology developed in this investigation by the Gear-Generalized Differential Quadrature Method (GGDQM), several numerical simulations are carried out using the Matlab software, in order to compare the numerical findings with our analytical solutions found for f ηη (0, 0), f ηη (1, 0), θ ηη (0, 0) and θ ηη (1, 0) (i.e., Equations (20) and (21)) and those obtained numerically by Motsa et al. [53] for f ηη (ξ, 0) via the Spectral Relaxation Method (SRM) and the Spectral Quasilinearization Method (SQLM), as shown in Tables 2-5. Also, GGDQM is further tested by computing the values of θ ηη (ξ, 0) from Equation (12) by considering the absence of viscous dissipation and Joule heating effects (i.e., Ec = 0). These values are evaluated for different values of Pr and compared in Table 1 with those obtained by Agbaje and Motsa [54] with the help of the Spectral Perturbation Method (SPM) and the Spectral Relaxation Method (SRM). As expected, it is found from Tables 1-5 that there is a good agreement between the compared results. Hence, the accuracy of our GGDQM numerical code is strengthened by validating our findings against the analytical and numerical results of some limiting cases. Furthermore, the numerical results listed in Table 5 for CPU time prove that GGDQM is a fast implementation method compared with SRM and SQLM, where the time required for GGDQM to generate the results shown in Table 5 is less than 7 s.

Results and Discussion
In this paper, the behaviors of velocity and temperature fields, skin friction coefficient, Nusselt number, entropy generation and Bejan number toward the involved pertinent parameters are examined numerically using the Gear-Generalized Differential Quadrature Method (GGDQM). In addition, the present numerical outputs are validated and discussed clearly and wittily via several graphical and tabular illustrations as shown in Figures 2-17 and Tables 1-6. Table 6 shows the effects of all physical parameters over the skin friction coefficient Re 1/2 x C f x and Nusselt number Re −1/2 x Nu x for the unsteady flow. It is observed that both Re 1/2 x C f x and Re −1/2 x Nu x decreases as time parameter ξ increases. Skin friction coefficient Re 1/2 x C f x increases; however, the Nusselt number Re −1/2 x Nu x decreases when increasing the Hartman number. Nusselt x Nu x increases for increasing values of Prandtl number, whereas the opposite behavior is seen for rising values of Eckert number. Figure 2 is sketched to see the variation of time parameter ξ. As ξ increases the velocity f η (ξ, η) and momentum boundary layer thickness decrease. The impact of magnetic interaction parameter M on the speed shape is observed in Figure 3. Boundary layer thickness and velocity decline as M increases. Physically, the magnetic field is correlated with the electrically conducting fluid and creates a Lorentz force that is opposite to the direction of fluid flow; consequently, fluid velocity decreases. Figure 4 shows that temperature θ(ξ, η) lessens by increasing ξ. Figure 5 exhibits the impact of Prandtl number on temperature θ(ξ, η). It is clear that temperature dropped with increasing Prandtl number. When increasing the Prandtl number (decreasing the thermal conductivity of fluid), the heat flow rate from the stretching boundary trims down and consequently the thermal boundary layer descends. Growing M decelerates the fluid flow; therefore, the friction between the fluid layers increases and generates frictional heating that raises the temperature (see Figure 6). The effect of Eckert number (Ec) on the temperature is shown in Figure 7. As expected, temperature mounts with mounting values of Ec. By increasing Ec, the resistance sandwiched between the fluid adjoining layers increases and leads to a change of the kinetic energy into thermal energy. Furthermore, as Ec increases, the width of the thermal boundary layer also increases. The influence of ξ on Ns is exposed in Figure 8. At the boundary and near to it entropy increases with ξ. However, after reaching a maximum value a decreasing trend is observed. Figure 9 displays the effects of Pr on entropy production number Ns. It is found that the entropy creation number is a rising function of Pr (due to high-temperature gradients) in the boundary layer flows. Figure 10 specifies the deviation of entropy creation number with dimensionless temperature Ω. It is observed that entropy decreases with increasing values of Ω. Hence, one can attain the main goal, that is, entropy creation minimization by reducing the working temperature disparity (T w − T ∞ ). Figure 11 specifies that entropy near the surface increases with M but after a certain η entropy decreases Ns. Figure 12 represents the effect of Ec on entropy production number Ns and increasing behavior is observed. Figure 13 demonstrates the impacts of ξ on Be. The figure shows that Be decreases with an increase in ξ. The impact of Prandtl on Bejan number Be is illustrated in Figure 14. With a large Prandtl number Pr Bejan number decreases. Figures 15 and 16, shows that Be declines when raising the dimensionless temperature difference Ω and Hartmann number M, respectively. Figure 17 portrays the Be for different values of Ec. Note that for Ec = 0, entropy creation is only due to heat transport. As the Eckert number is inversely proportional to the temperature difference, for high temperature dissimilarity between the surface and the ambient fluid, the viscous dissipation parameter becomes zero (Ec = 0.0). Therefore, the heat transfer irreversibility in the entire flow region is completely dominant, i.e., (Be = 1.0). Table 1. Comparison of our numerical results with those obtained by Agbaje and Motsa [54] for θ η (ξ, 0) at different values of Pr, when ξ = 0.5, Ec = 0 and M = 1.

Closing Remarks
The impacts of energy dissipation and magnetic field on heat transfer and entropy generation are analyzed by utilizing a new hybrid numerical technique called gear-generalized differential quadrature method (GGDQM). The flow driven by a stretching boundary is assumed to be unsteady. Following are the key findings: Fluid decelerates with the enhancement of reduced dimensional time and magnetic parameter. Increase in reduced dimensional time and Prandtl number reduced the fluid temperature and reverse behavior was observed with rising values of Eckert number and magnetic parameter.

Closing Remarks
The impacts of energy dissipation and magnetic field on heat transfer and entropy generation are analyzed by utilizing a new hybrid numerical technique called gear-generalized differential quadrature method (GGDQM). The flow driven by a stretching boundary is assumed to be unsteady. Following are the key findings: Fluid decelerates with the enhancement of reduced dimensional time and magnetic parameter. Increase in reduced dimensional time and Prandtl number reduced the fluid temperature and reverse behavior was observed with rising values of Eckert number and magnetic parameter.
Entropy generation number Ns, rises with enhancing values of reduced dimensionless time at the boundary and its vicinity. The effects become opposite after the certain distance from the stretching boundary.
Entropy generation number Ns, enhances with magnetic parameter, Prandtl and Eckert number. Reduction in Ns is observed with rising values of Ω.
The effects of emerging parameters on Ns are significantly prominent at the surface of stretching boundary. Decrement in Bejan number Be is observed at the surface of stretching sheet with enhancing values of Eckert number, magnetic parameter, Prandtl number and temperature difference parameter. Funding: This research received no external funding. The APC was given by Ton Duc Thang University, Ho Chi Minh City, Vietnam. However, no grant number is available from source.