On the Stability of MHD Boundary Layer Flow over a Stretching/Shrinking Wedge

The steady two dimensional magnetohydrodynamic (MHD) boundary layer flow and heat transfer over a stretching/shrinking permeable wedge is numerically investigated. The partial differential equations governing the flow and heat transfer are transformed into a system of ordinary differential equations using a similarity transformation. These equations are then solved numerically using the boundary value problem solver, bvp4c in Matlab software. It is found that dual solutions exist for a certain range of the shrinking strength. A stability analysis is performed to identify which solution is stable and physically reliable.

The steady two dimensional magnetohydrodynamic (MHD) boundary layer flow and heat transfer over a stretching/shrinking permeable wedge is numerically investigated. The partial differential equations governing the flow and heat transfer are transformed into a system of ordinary differential equations using a similarity transformation. These equations are then solved numerically using the boundary value problem solver, bvp4c in Matlab software. It is found that dual solutions exist for a certain range of the shrinking strength. A stability analysis is performed to identify which solution is stable and physically reliable.
One of the family of boundary layer similarity solutions was discovered by Falkner and Skan 1 , for the study of flow over a static wedge. The numerical solutions were then calculated by Hartree 2 . With a similarity transformation, the boundary layer equation reduces to an ordinary differential equation, which is well known as the Falkner-Skan equation. This problem reduces to the classical Blasius flow if the angle of the wedge is set to zero, while it reduces to the stagnation-point flow (Hiemenz flow) if the angle is 180°. The Falkner-Skan flow past a stretching wedge was considered by Postelnicu and Pop 3 and Su et al. 4 , while very recently Nadeem et al. 5 studied the Falkner-Skan problem for a static and moving wedge in the presence of induced magnetic field.
Magnetohydrodynamic (MHD) flow is the study of magnetic properties on electrically conducting fluids. Given its significance in engineering industrial applications, this has led to numerous studies of MHD flow in electrically conducting fluids, such as Ishak 6 , Jafar et al. 7 , Mat Yasin et al. 8 , Hayat et al. 9,10 , Lu et al. 11 , etc. Dual solutions were found for the unsteady flow over a shrinking sheet by Soid et al. 12 . Meanwhile, Sharma et al. 13 in their study of two dimensional MHD stagnation point flow over a stretching/shrinking sheet found dual solutions exist for a certain range of the shrinking parameter. The stability analysis revealed that the first solution is stable while the second solution is not.
With that, the present study examines the steady boundary layer flow and heat transfer over a stretching/ shrinking wedge immersed in a viscous and incompressible electrically conducting fluid, with uniform surface temperature. The analysis is conducted numerically by using bvp4c function and once the dual solutions are obtained, the stability analysis is performed to determine which solution is stable and physically reliable. Figure 1 depicts the physical model and coordinate system for the stretching/shrinking wedge, where x and y are respectively the Cartesian coordinates measured along the surface and normal to it, while u and v are the velocity components along the Cartesian coordinates x and y, respectively. It is assumed that the velocity of the stretching/ shrinking wedge is represented by

Mathematical Formulation
corresponds to stretching and < U 0 w corresponds to shrinking. It is also assumed that the velocity of the external flow outside the boundary layer is represented by where m and U ∞ are positive constants. A variable magnetic field of strength represented by B(x) is applied in the positive direction of y-axis. The induced magnetic field is assumed small and neglected. In Fig. 1, β is the Hartree pressure gradient parameter which corresponds to β = Ω/π for a total angle Ω of the wedge. We note that ≤ ≤ m 0 1 with m = 0 for the boundary-layer flow over a stationary flat plate (Blasius problem) and m = 1 for the flow near the stagnation point on an infinite wall.
The governing equations based on the above-mentioned assumptions are (Ishak et al. 14

)
Scientific REPORtS | (2018) 8:13622 | DOI:10.1038/s41598-018-31777-9 while the boundary conditions are where ν is the kinematic viscosity, T is the fluid temperature, T w is the uniform surface temperature, σ is the electrical conductivity, α is the thermal diffusivity of the fluid, ρ is the fluid density and v w (x) is the mass flux velocity with < v x ( ) 0 w for suction and > v x ( ) 0 w for injection. To obtain similarity solutions of equations (1) to (3) subject to boundary conditions (4), we introduce the following similarity variables (Ishak et al. 15,16 ): where ψ is the stream function, which is defined as u = ∂ψ/∂y and v = −∂ψ/∂x which satisfies the continuity equation (1). Thus, resulting to: where prime denotes differentiation with respect to η. At the boundary where η = 0, the transpiration rate is given by where s = f(0), a constant parameter with s > 0 for suction and s < 0 for injection. To obtain similarity solution, all parameters must be constant. For this purpose, we take = 17 ).
Substituting equations (5) and (6) into equations (2) and (3), the following ordinary (similarity) differential equations are obtained: which is subject to the boundary conditions: where λ is the stretching/shrinking parameter, β is the pressure gradient parameter, M is the magnetic parameter (Hartmann number), Pr is the Prandtl number and s is the suction/injection parameter, which are defined as: The physical quantities of interest in this study are the skin friction coefficient C f and the local Nusselt number Nu x which can be expressed as: where τ w is the wall shear stress along the stretching/shrinking surface and q w is the surface heat flux, which are defined as: w y w y 0 0 Substituting (6) into (13), one obtains where Re x = u e (x)x/ν is the local Reynolds number.

Results and Discussions
In this study, the nonlinear ordinary differential equations (8) and (9) subject to the boundary conditions (10) were solved numerically using the bvp4c function available in Matlab software (see Shampine et al. 18 ). The relative error tolerance was set to10 −7 with the configured convergence criterion to achieve an accuracy of six decimal places. Dual solutions exist in the current study. These solutions are obtained by using different initial guesses for the values of f ″(0) and θ′(0) where all of the velocity and temperature profiles satisfy the far field boundary conditions (equation 10) asymptotically. For this purpose, we need to set different values of the boundary layer thickness, η ∞ . The range of η ∞ for the first solution is between 10 and 20 while for the second solution is between 30 and 40. The initial guesses and the boundary layer thickness vary for different values of parameters. In the numerical computations, we have to make sure that all profiles satisfy the infinity boundary conditions asymptotically, otherwise the numerical results are not valid. Besides, we compare our results with previously published one by other researchers. Table 1 shows the comparison of numerical results of the skin friction coefficient, f ″(0) with those reported by White 19 , Ishak et al. 14 and Postelnicu and Pop 3 for the case when λ = 0 (fixed surface), M = 0 (magnetic field is absent) and s = 0 (impermeable surface), which shows a favorable agreement. This gives confidence to the other numerical results that reported in the present study. For the sake of brevity, the value of Prandtl numer, Pr was set as 1 (such as ionized gases) and magnetic parameter, M was set as 0.2.  which equations (8-10) have no solutions for λ λ < c . The values of λ c are given in Figs 2 and 3. It is observed that dual solutions are found for the shrinking case, while for the stretching case, the solution is unique. Figure 2 also shows that the skin friction coefficient f ″(0) increases when the suction strength is increased. Meanwhile, Fig. 3 reveals that the local Nusselt number, −θ′(0) which represents the heat transfer rate at the surface, decreases as the effect of suction increases.     solution is lower than that of the second solution and the infinity boundary conditions (f → 1 and θ → 0 as η → ∞) are satisfied asymptotically which support the validity of the numerical results obtained. increases as the angle of the wedge increases. Besides that, the results also indicate that the skin friction coefficient f ″(0) increases as β increases for the first solutions, whereas it decreases for the second solution. Similar behaviour is observed for the absolute local Nusselt number, |−θ′(0)| as presented in Fig. 7.
The graphical results for the velocity profiles f ′(η) are presented in Fig. 8, while Fig. 9 graphically illustrates the temperature profiles θ(η) for different values of λ for β = 0.33, s = 1, M = 0.2 and Pr = 1. Both profiles reveal that the boundary layer thickness for the second solution is higher than that of the first solution and both satisfy the infinity boundary conditions asymptotically, which support the validity of the numerical results obtained.
Stability Analysis. Since the solutions are not unique for a certain range of the shrinking strength, a stability analysis is performed to determine which one of the solutions is stable and thus physically reliable. There have been several studies on stability analysis such as Merkin 20 , Weidman et al. 21 , Roşca and Pop 22 , Merril et al. 23 and Awaludin et al. 24 . An analysis of the temporal stability is performed to identify which solution is stable when time passes. The new governing equations for the unsteady state flow are given as where t denotes the time while equation (1) holds. The new similarity variables are introduced as where τ denotes dimensionless time.
Substituting equation (17) into equations (15) and (16) yields which is subject to the boundary conditions: 0 as (20) To test the stability of the steady flow solution f(η) = f 0 (η) and θ(η) = θ 0 (η) of equations (8) to (10) we write (see Merril et al. 23 ): where γ is an unknown eigenvalue; F(η) and G(η) are small relative to f 0 (η) and θ 0 (η). We take perturbation in exponential form since it will increase or decrease more rapidly compared to the power functions. This form of perturbation has also been considered by Merkin 20 , Weidman et al. 21 , Roşca and Pop 22 , Merril et al. 23 and Awaludin et al. 24 , among others. Equations (18) to (20) present an infinite set of eigenvalues, γ γ < < ... 1 2 . If the smallest eigenvalue is positive, there is an initial decay which is the flow is stable. Besides, there is an initial growth of disturbance if the smallest eigenvalue is negative which mean the flow is unstable. Substituting equation (21) into equations (18) and (19), the following linear eigenvalue problems are obtained: along with the boundary conditions: The stability of the steady state flow are determined by the smallest eigenvalue γ 1 . Since equations (22) and (23) satisfy homogeneous boundary and far-field conditions (24), without loss of generality we set F″(0) = 1 to determine the eigenvalues γ. There are more that one value of γ in equations (22) and (23) that give F″(0) = 1. A search for the lowest eigenvalues γ 1 satisfying equations (22) to (24) was carried out and the results are plotted in Fig. 10.
Based on equation (21), the unsteady solution f(η, τ) converges to the steady state solution F(η) as time passes (τ → ∞) if γ 1 is positive. The sign of γ 1 determines either there is an initial decay or initial growth of disturbance. Figure 10 shows that γ 1 is positive for the first solution, and negative for the second solution. Thus, the first solution is stable and physically reliable, whereas the second solution is unstable.

Conclusions
This study examined the steady two dimensional MHD flow and heat transfer of an incompressible and electrically conducting fluid over a permeable stretching/shrinking wedge. Numerical results, which were obtained using the bvp4c function in Matlab software, proved the existence of dual solutions for the shrinking case but unique solution for the stretching case. For the shrinking case, the existence of the solution depends on the shrinking strength λ and the angle of the wedge Ω. The range of λ (for which the solution exists) increases as the angle of the wedge Ω increases. Discussion on the skin friction coefficient and the local Nusselt number were carried out for the effect of all parameters involved. The boundary layer thickness for the second solution was found to be higher than that of the first solution, which gave a picture of its stability. To confirm this observation, a stability analysis was conducted and it was revealed that the first solution is stable and physically reliable while the second solution is not.