Unsteady MHD Mixed Convection Slip Flow of Casson Fluid over Nonlinearly Stretching Sheet Embedded in a Porous Medium with Chemical Reaction, Thermal Radiation, Heat Generation/Absorption and Convective Boundary Conditions

Numerical results are presented for the effect of first order chemical reaction and thermal radiation on mixed convection flow of Casson fluid in the presence of magnetic field. The flow is generated due to unsteady nonlinearly stretching sheet placed inside a porous medium. Convective conditions on wall temperature and wall concentration are also employed in the investigation. The governing partial differential equations are converted to ordinary differential equations using suitable transformations and then solved numerically via Keller-box method. It is noticed that fluid velocity rises with increase in radiation parameter in the case of assisting flow and is opposite in the case of opposing fluid while radiation parameter has no effect on fluid velocity in the forced convection. It is also seen that fluid velocity and concentration enhances in the case of generative chemical reaction whereas both profiles reduces in the case of destructive chemical reaction. Further, increase in local unsteadiness parameter reduces fluid velocity, temperature and concentration. Over all the effects of physical parameters on fluid velocity, temperature and concentration distribution as well as on the wall shear stress, heat and mass transfer rates are discussed in detail.


Introduction
Boundary layer flow over a stretching surface has several engineering and industrial applications, for example, in cooling bath, Aerodynamic extrusion, plastic sheet, metallurgical process, glass blowing, manufacturing of rubber and plastic sheets, crystal growing, the sheets are continuously stretched. During manufacturing, sheets are stretched continuously in order to achieve the desired thickness. It shows that final product depends upon the stretching and On the other hand, non-Newtonian fluids has gained attraction due to its wide range application in various industries, such as design of solid matrix heat, nuclear waste disposal, chemical catalytic reactors, geothermal energy production, ground water hydrology, transpiration cooling, petroleum reservoirs, etc. These fluids are more complicated as compared to Newtonian fluids due to nonlinear relation between stress and strain rate. Several models have been proposed for the study of non-Newtonian fluids, but yet not a single model is developed that exhibits all properties of non-Newtonian fluids. In literature, the simplest model is the Maxwell model. Among numerous non-Newtonian fluids, there is another fluid known as Casson fluid. Casson fluid is a shear thinning fluid which is assumed to have an infinite viscosity at zero rate of shear, a yield stress below which no flow occurs and a zero viscosity at an infinite rate of shear. The Casson model also called rheological model was originally developed by Casson [35] for viscous suspension of cylindrical particles [36]. Common examples of Casson fluid are honey, jelly, soup, tomato sauce, concentrated fruit juices, etc. Also, it is the most appropriate rheological model for blood and chocolate. Furthermore, Casson fluid possesses yield stress and has great importance in polymer processing industries and biomechanics. The influence of thermal radiation on unsteady flow of Casson fluid caused by stretching sheet subjected to suction/blowing was presented by Mukhopadhyay [37]. Nadeem et al. [38] also analyzed the three dimensional hydromagnetic flow of Casson fluid in a porous medium. Numerical solutions of electrically conducting slip flow of Casson nanofluid generated during stretching sheet under the influence of convective boundary condition using similarity transformations were presented by Ibrahim and Makinde [39]. Motivated by its applications and interesting features, Benazir et al. [40] studied unsteady Casson flow past a vertical cone and flat sheet in the presence of magnetic field. Very recently, Oyelakin et al. [41] investigated unsteady electrically conducting flow of Casson nanofluid in the presence of slip and convective boundary conditions. Motivated by the above study on steady and unsteady stretching sheet and its widespread engineering and industry application, the objective of present investigation is to study the effects of chemical reaction and thermal radiation on unsteady electrically conducting flow of Casson fluid towards nonlinearly stretching sheet saturated in a porous medium in the presence of slip and convective boundary conditions. The highly nonlinear partial differential equations are transformed into system of ordinary differential equations via suitable transformations and then numerically solved using Keller-box method [42]. The numerical and graphical results are obtained by developed an algorithm in MATLAB software. The developed algorithm for present problem is validated through comparison with previously published results and excellent accuracy is achieved. The physical behavior of pertinent parameters on fluid velocity, temperature and concentration are highlighted and discussed.

Mathematical Formulation
Consider hydromagnetic mixed convection flow of Casson fluid past an unsteady nonlinearly stretching sheet through porous medium under the influence of chemical reaction, thermal radiation, convective boundary conditions and partial slip. The x-axis is taken along the direction of stretching sheet and y-axis is perpendicular to the surface (see Fig 1). Initially at t 0 both the sheet and fluid are at rest and at the same temperature T 1 and concentration C 1 . Suddenly, the sheet is stretched with the nonlinear velocity of the form u w (x,t) = cx n /(1 − γt), where c, γ ! 0 are constants, t is time and n(>0) represents the nonlinearly stretching sheet parameter such that, n = 1 corresponds to linear stretching sheet and n 6 ¼ 1 represents nonlinear stretching sheet.
Moreover, a magnetic field Bðx; tÞ ¼ B 0 x ðnÀ 1Þ=2 ð1 À gtÞ À 1 2 = (being the function of x and t) is applied perpendicular to the stretching sheet with constant B 0 . Furthermore, it is also assumed that sheet wall is heated by temperature T f (x,t) and concentration C s (x,t). The temperature and nanoparticles concentration at free stream is T 1 and C 1 , respectively.
The rheological equation of an isotropic and incompressible flow of Casson nanofluid can be written as (see Mukhopadhyay [37] and Oyelakin et al. [41]) Here, π = e ij e ij and e ij is the (i,j)-th component of the deformation rate, π is the product of the component of deformation rate with itself, π 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. The governing equations for mixed convection flow of Casson fluid along with continuity equation are given as @u @x þ @u @y ¼ 0; In the above expressions u and υ denote the velocity components in x and y direction respectively, ν is kinematic viscosity, σ is the electrical conductivity, β is the Casson parameter, ρ is the fluid density, φ is the porosity, k 1 = k 0 (1 − γt)/x (n−1) is the variable permeability of porous medium, g is the gravitational force due to acceleration, β T is the volumetric coefficient of thermal expansion, β C the coefficient of concentration expansion, k is the thermal conductivity of the Casson fluid, q r is the radiative heat flux, Qðx; tÞ ð1À gtÞ is heat generation/absorption coefficient, D is the mass diffusivity and k c is the rate of chemical reaction.
The corresponding boundary conditions are written as follows: Here N 1 ðx; tÞ ¼ N 0 x À ðnÀ 1Þ= 2 ð1 À gtÞ The radiative heat flux q r described according to Rosseland approximation is give as where σ Ã is the Stefan-Boltzmann constant and k 1 Ã is the mean absorption coefficient. T 4 can be expressed as linear function of temperature. By expanding T 4 in a Taylor series about T 1 and neglecting higher terms, we can write Incorporating Eq (9) and Eq (10) in Eq (4), we obtain @T @t þ u @T @x þ u @T @y Now introduce the stream function ψ defined in its usual notation in terms of velocity, a variable η and the following transformations; The system of Eqs (3-7) and Eq (11) take the following form In the above expressions, A, M, K, l ¼ AE Gr x Re 2 x , N, δ, Pr, R d , Ec, ε, Bi 1 , Bi 2 , Sc and R are the local unsteadiness parameter, magnetic parameter, porosity parameter,thermal buoyancy parameter (λ > 0 corresponds to assisting flow, λ = 0 indicates no convection and λ < 0 denotes the opposing flow)(Gr x and Re x being Grashof number and local Reynold number respectively), buoyancy ratio parameter, slip parameter, Prandtl number, radiation parameter, Eckert number, heat generation/absorption parameter, Biot numbers, Schmidt number and chemical reaction parameter, and are defined as ; The wall skin friction, wall heat flux and wall mass flux, respectively, are defined by The dimensionless skin friction coefficient Cf aðT f À T 1 Þ and local Sherwood number Sh x ¼ xq s DðC s À C 1 Þ on the surface along x-direction, local Nusselt number Nu x and Sherwood number Sh x are given by

Numerical Scheme
The highly nonlinear governing Eqs (13-15) along with the associated boundary conditions (16) and (17) are solved numerically via Keller-box method. The details of this method can be found in the book of Cebeci and Bradshaw [42]. Like several other finite difference methods, the Keller-box method is very effective as it allow us to get the numerical solutions to system of nonlinear differential equations. Among many other numerical techniques, the finite difference method is more flexible for the reason that initial approximations control the convergence rate.
The four main steps are involved to obtain the numerical solutions via Keller-box method, and are as follow: 1. Reduce Eqs (13-15) into a system of first order equations; 2. Write the difference equations using central differences; 3. Linearize the algebraic equations via Newton's method, and write them in matrix-vector form; and 4. Finally, solve the obtained linear system by tridiagonal elimination technique.
It is worth mentioning that uniform step size of Δη = 0.01 is found to be satisfactory in obtaining the numerical solutions for each profile with an error tolerance 10 −5 .

Results and Discussion
In the present study, unsteady mixed convection flow of Casson fluid due to nonlinearly stretching sheet through porous medium under the influence of magnetic field and chemical reaction is explored. Moreover, the influence of partial slip and convective boundary conditions is also considered. Numerical computations are carried out for various values of local unsteadiness parameter A, Casson fluid parameter β, nonlinear stretching sheet parameter n, magnetic parameter M, porosity parameter K, thermal buoyancy parameter γ, concentration buoyancy ratio parameter N, Prandtl number Pr, radiation parameter R d , Schmidt number Sc, slip parameter δ and Biot numbers Bi 1 , Bi 2 .The present results are compared and validated in Tables (1-3) with some existing results in the available literature for limiting cases.
Tables 1 and 2 present the comparison of skin friction coefficient for different values ofβ, M and A, respectively, with the results of Nadeem et al. [37], Oyelakin et al. [43], Chamkha et al. [23], Sharidan et al. [21] and Mukhopadhyay et al. [36]. The results showed an excellent agreement. Table 3 describes the comparison of heat transfer rate for increasing values of Pr with the results of Grubka and Bobba [4], Aurangzaib et al. [30] and Mabood et al. [32], and revealed in a good agreement. . It is found that fluid velocity reduces with increase in A for both fluids. It is also observed that larger amplitude of A thinning the boundary layer. Fig 3 presents that fluid velocity decelerates as β increases for all cases of λ. It is noticed that velocity boundary layer thickness reduces faster in the case of opposing flow (λ < 0) when β approaches towards infinity. Physically it makes sense because plasticity of fluid is higher as β goes higher and fluid experiences a resistance. Consequently, the velocity boundary layer becomes thinner for larger values of β. Fig 4 demonstrates the variation of n on velocity profile for both A = 0 and A 6 ¼ 0. In both cases, fluid velocity falls with increase in n. However, the velocity boundary layer thickness rapidly reduces with n when A 6 ¼ 0. Fig 5 clears that increasing values of M decrease the fluid velocity in all the three cases of γ. It is obvious, because the strength of magnetic field presents a resistance to the flow and results a decline in velocity boundary layer thickness. Moreover, a rapid decrease is seen in the fluid velocity with increase in M when λ < 0. Fig 6 reveals the effect of K on velocity profile for n = 1 and n > 1. It is seen that fluid velocity is lower with increase in K. Physically, the holes inside porous medium become wider as K increases, in this case fluid experiences a drag force which act opposite to the direction of flow and results a reduction in velocity boundary layer thickness. The influence of λ on velocity profile for K = 0 and K 6 ¼ 0 is displayed in Fig 7. It is observed that fluid velocity is higher in the case of assisting flow (λ > 0) and lower in the case of opposing flow (λ < 0). This is physically realistic because in assisting flow buoyancy force dominates the viscous force which results an increment in the velocity, and in opposing flow the buoyancy force act opposite to the flow direction due to   which the thickness of velocity boundary layer declines. On the basis of this phenomenon, it may be concluded that assisting flow thickens the velocity boundary layer whereas the velocity boundary layer becomes thinner in opposing flow. A similar kind of trend is observed for the effect of N on velocity profile in the absence and presence of magnetic field (see Fig 8). The influence of δ on velocity profile for all cases of λ is described in Fig 9. It is observed that for all the three case of λ, velocity is a decreasing function of δ. However, in the case of assisting flow the fluid velocity initially decrease in the vicinity of stretching sheet and then increases near the free stream as δ increased. This phenomenon can be explained as that increasing δ permitting more fluid to slip over the sheet due to which the flow near the sheet reduced and the slip effect towards the free stream is less pronounced. Fig 10 illustrates the variation of fluid velocity for various values of R d when λ < 0, λ = 0 and λ > 0. It is interesting to notice that increasing values of R d decline fluid velocity when λ < 0, has no effect when λ = 0 and enhances when λ > 0. Physically, it can be explained as that in assisting flow the convective heat transfer effect   Next, Fig 11 shows the influence of the chemical reaction parameter R on velocity profile for various values λ. It is observed that fluid velocity slightly enhances in the case of generative chemical reaction (R < 0) and reverse effect in the case of destructive chemical reaction (R > 0). Hence velocity boundary layer slightly expands for generative chemical reaction and becomes thinner in case of destructive chemical reaction.
Figs 12-24 examined the variation of dimensionless temperature profile (θ(η)) for various values of A, β, n, M, K, γ, N, δ, Pr, R d , Ec, ε and Bi 1 . Fig 12 demonstrates the effect of A on temperature profile for M = 0 and M 6 ¼ 0. It can be easily noticed that temperature decrease as A increased. Also, thermal boundary thickness in the case of A 6 ¼ 0 is lower than A = 0. Clearly, Fig 13 reveals that increase in β lead to higher temperature across boundary region for both n = 1 and n 6 ¼ 1. This phenomenon indicates that decrease in yield stress thickens thermal boundary layer. It is also observed that the influence of β is more pronounced for unsteady linear stretching sheet. The effect of n on temperature profile for M = 0 and M 6 ¼ 0 is depicted in Fig 14. In both case the temperature is found to be reduced with increase in n. The thermal boundary layer becomes thinner when the sheet is stretched in a nonlinear way. Fig 15 elucidates the variation of temperature profile for various values of M when β = 0.6 (Casson fluid) and β ! 1 (Newtonian fluid). It is found that increasing values of M slightly enhance the temperature for both fluids. Since electromagnetic forces dominate viscous forces for higher values of M, it results in thickening thermal boundary layer. The influence of K on temperature profile for n = 1 and n 6 ¼ 1 is plotted in Fig 16. The temperature is found to be increasing function of K. The influence of K is more prominent for steady flow and thermal boundary layer thicker as compared to unsteady flow. Fig 17 displays influenced of buoyancy parameter λ on the temperature. The thermal boundary layer becomes smaller due to increase in buoyancy parameter. It is caused since the buoyancy forces enhance temperature gradient and consequently the thermal boundary layer thickness reduces. A similar pattern is observed for the effect of N on temperature profile, i.e. thickness of thermal boundary layer reduces as N increases (see Fig 18). It is also noticed from this figure that temperature is more influenced by Nin unsteady linear stretching sheet. The variation of temperature profile for various values of δ when A = 0 and A 6 ¼ 0 is displayed in Fig 19. It is interesting to note that temperature is higher for larger values of δ. It is obvious because increasing δ allow more fluid to slip past a sheet, this results a rise in thermal boundary layer thickness. It is worth mentioning here that the effect is more pronounced in steady flow. Fig 20 portrays the effect of Pr (i.e. Pr = 0.67, 0.71, 7, 11.4 for Argon at 25°C, condensed air, water and water at 4°C) on temperature profile for both A = 0 and A 6 ¼ 0. Prandtl number is the ratio of momentum diffusivity to thermal diffusivity. In both cases temperature reduces with increase in Pr. It is an agreement with the fact that higher values of Pr correspond to weaker thermal diffusivity which therefore thinning the thermal boundary layer. In addition, Pr is effectively used to control the both boundary thicknesses in heat transfer process. It is also noted from this figure that Argon at 25°C has thicker thermal boundary layer whereas water at 4°C has thinner boundary layer. The influence of R d on temperature profile for n = 1 and n 6 ¼ 1 is exhibited in Fig 21. It is observed that larger values of R d rise the temperature significantly. Physically, it is true due to the fact that higher values of R d lead to large radiation effect and thereby enhancing the thickness of thermal boundary layer. In the work [43] similar temperature pattern were obtained. Fig 22 portrays the variation of dimensionless temperature profile for various values of Ec when n = 1 and n 6 ¼ 1. It can be easily noticed from this figure that temperature is higher for larger values of Ec. It is well known that viscous dissipation generates heat due to frictional heating between the fluid particles and this extra heat lead to temperature and associated boundary layer thickness enhancement. The effect of ε on dimensional temperature for A = 0 and A 6 ¼ 0 is shown in Fig 23. It is worth mentioning here that ε > 0 represent heat generation and ε < 0 corresponds to heat absorption. Clearly, the existence of heat generation in the boundary layer leads to enhance the energy level. Consequently, thickens thermal boundary in the vicinity of sheet. On the other hand, opposite trend is observed in the presence of heat absorption, i.e. thermal boundary layer becomes thinner. It is also noticeable from this figure that the effect is more pronounced for steady flows. Fig 24  examined the effect of Bi 1 on temperature profile for both A = 0 and A 6 ¼ 0. It is noteworthy that Bi 1 << 1 shows that temperature field is uniform inside the surface whereas Bi 1 >> 1 indicates the non-uniformity of temperature field inside the body surface. Biot number Bi 1 is the ratio of internal thermal resistance of the sheet to the boundary layer thermal resistance of the hot fluid at the bottom of the surface. In both cases, temperature is an increasing function of Bi 1 . The reason behind this is that strength of Bi 1 increases the convection at the surface and thereby enhances the surface temperature. The reduction in concentration boundary layer thickness is found for non-linearly stretching sheet. The influence of M on concentration profile for A = 0 and A 6 ¼ 0 is depicted in Fig 28. It is noticeable that concentration distribution related boundary layer thickness enhance with increase in M. The reason already explained in Fig 5, i.e. strength of magnetic field enhances the Lorentz force which therefore reduces the velocity and results an increase in concentration boundary layer thickness. It is also observed that the effect of M is large when A = 0 (steady flow). A similar behavior of concentration profile is noted for increasing values of K when n = 1 and n 6 ¼ 1 (see Fig 29). It is also found from this figure that the effect is more pronounced in linear stretching sheet. From Fig 30, it is noticed that concentration distribution is a decreasing function of λ for both K = 0 and K 6 ¼ 0. However, the effect is very small due to the reason that the buoyancy term is only coupled with momentum equation. A same kind of effect has been observed for increasing values of N on concentration profile (see Fig 31). The influence of δ on concentration distribution for A = 0 and A 6 ¼ 0 is displayed in Fig 32. This figure reveals that in both cases concentration and associated boundary layer thickness rise as δ increases. Since a bulk of fluid is slipping over the sheet which reduces the flow field and thereby thickening the concentration boundary layer. It is also noticed from this figure that concentration is more influenced by δ in the case of steady flow. Fig 33 clears that concentration slightly reduces with increase in R d . The reason behind this phenomenon is that stronger R d lead to enhance the fluid velocity as well as temperature, and therefore concentration boundary layer becomes thinner. Fig 34 discusses   Unsteady MHD Mixed Convection Slip Flow of Casson Fluid over Nonlinearly Stretching when A = 0 and A 6 ¼ 0. In both cases, the concentration distribution reduces with increase in Sc. Physically, larger values of Sc correspond to low mass diffusivity which leads to a thinning of the concentration boundary layer. The influence of Sc on fluid velocity and temperature is very small because the physical parameter arises only in concentration equation and is not displayed here for the sake of brevity. Clearly, Fig 35 describes that concentration rises in generative chemical reaction (R < 0) and opposite to this is observed in the case of destructive chemical reaction (R > 0). The physical reason behind this behavior is that strength of chemical reaction alters the diffusion rates. It is also noticeable that the effect is major for steady flow. The variation of concentration profile for different values of Bi 2 when n = 1 and n 6 ¼ 1 is portrayed in Fig 36. It can be easily seen that increasing values of Bi 2 enhances the concentration distribution for both n = 1 and n 6 ¼ 1. Since the concentration distribution is driven by temperature field. Hence, higher values of Bi 2 would promote a deeper penetration of the concentration.  stress enhances with increase in λ and δ, while the wall shear stress initially increases with the increase in n for opposing flow (λ < 0) for larger slip rate and it reduces for assisting flow (λ > 0). Fig 39 presents the variation of local Nusselt number for various values of A, β and M. It is noticeable that increasing values of A enhances the heat transfer rate whereas it slightly reduces with increase in β and M. The influence of for different values of R d , δ and Pr on the local Nusselt number is displayed in Fig 40. It is observed that heat transfer rate rises with increase in R d and Pr whereas increasing values of δ decrease the heat transfer rate. Finally, Fig 41 presents the variation of the local Sherwood number for different values of Sc, A and R. It can be easily seen that Sherwood number, i.e., the mass transfer rate is higher for the increasing values of Sc, A and R.

Conclusions
The objective of present study is to analyze unsteady hydromagnetic mixed convection flow of Casson fluid towards nonlinearly stretching sheet saturated in a porous medium under the influence of chemical and convective boundary conditions. The influence of viscous dissipation, joule heating, heat generation/absorption and partial slip is also considered in the present problem. The resulting governing equations are solved numerically using Keller-box method. The numerical as well as graphical results are obtained by developing an algorithm in MATLAB software. Physically, the effect of local unsteadiness parameter A, Casson parameter β, nonlinear stretching parameter n, magnetic parameter M, porosity parameter K, thermal buoyancy parameter γ, buoyancy ratio parameter N, Prandtl number Pr, radiation parameter R d , Eckert number Ec, heat generation/absorption parameter ε, Schmidt number Sc, chemical reaction parameter R, slip parameter δ and Biot numbers Bi 1 , Bi 2 on fluid velocity, temperature, concentration as well as wall shear stress, heat and mass transfer rates are investigated. The main findings of this study can be summarized as follows: 1. The magnitude of wall shear stress, heat and mass transfer rates are found to be enhanced with the increase in A.

2.
For increasing values of β, the velocity boundary layer becomes thinner whereas both thermal and concentration boundary layer become thicker with increase in β.
3. The fluid velocity reduces in the case of opposing flow (λ < 0) while enhances in the case of assisting flow (λ > 0) and opposite to this is observed for temperature and concentration distributions.
4. The fluid velocity is observed as increasing function of R d when λ > 0 whereas decreasing function of R d when λ < 0 and has no effect on velocity when λ = 0.
5. The velocity as well as concentration boundary layer thicknesses are noticed to be decreased with increase in R.
6. The small variations of temperature and concentration distributions are observed for the effect of M and K in the case of unsteady flow.
7. The variation of temperature for increasing values of ε is more pronounced for the steady flow.