On New Numerical Techniques for the MHD Flow Past a Shrinking Sheet with Heat and Mass Transfer in the Presence of a Chemical Reaction

We use recent innovative solution techniques to investigate the problem of MHD viscous flow due to a shrinking sheet with a chemical reaction. A comparison is made of the convergence rates, ease of use, and expensiveness the number of iterations required to give convergent results of three seminumerical techniques in solving systems of nonlinear boundary value problems. The results were validated using a multistep, multimethod approach comprising the use of the shooting method, the Matlab bvp4c numerical routine, and with results in the literature.


Introduction
Boundary layer flow over a stretching surface occurs in several engineering processes such as hot rolling, wire drawing, and glass-fibre production.Materials that are manufactured by extrusion processes and heat-treated substances proceeding between a feed roll and a windup roll can be classified as a continuously stretching surface 1-3 .A shrinking film is useful in the packaging of bulk products since it can be unwrapped easily with adequate heat 4-7 .Shrinking problems can also be applied to the study of capillary effects in small pores and the hydraulic properties of agricultural clay soils 8 .Studies of flow due to a shrinking sheet with heat transfer and/or mass transfer have been considered by, among others, 7, 9 .
In recent years, several analytical or semianalytical methods have been proposed and used to find solutions to most nonlinear equations.These methods include the Adomian decomposition method ADM 10, 11 , differential transform method DTM 12 ,

Mathematical Formulation
We investigate the effect of chemical reaction, heat and mass transfer on nonlinear MHD boundary layer past a porous shrinking sheet with suction.The governing boundary layer equations of momentum, energy, and mass diffusion in terms of the velocity components u, v, and w are see Muhaimin et where α is the thermal conductivity of the fluid, B 0 is the magnetic field, κ is the thermal viscosity, K is the permeability of the porous medium, k 1 is the rate of chemical reaction, ν μ/ρ is the kinetic viscosity, μ is the dynamic viscosity, and σ is the electrical conductivity.
The applicable boundary conditions are where a > 0 is the shrinking constant and W is the suction velocity.The cases m 1 and m 2 correspond to shrinking sheets in the x-and y-directions, respectively.
Using the similarity transformations see Sajid and Hayat 32 : 2.1 are transformed to the system of nonlinear equations where Pr ν/κ is the Prandtl number, Sc ν/D is the Schmidt number, λ κ/Ka is the porosity, and γ is the chemical reaction parameter.We remark that 2.4 can be solved independently of equations of 2.5 -2.6 for f, but the solutions for θ and φ directly depend on the solution for f.To demonstrate how robust the proposed methods of solution are, the system of 2.4 -2.5 is solved simultaneously in the next section.Solving the equations simultaneously is also important when conducting the parametric study because some of the governing parameters such as Pr and m affect all three unknown variables.

Solution Methods
We solve 2.4 -2.6 using three recent innovative semi-numerical methods.Validation of the results is done by further solving the equations numerically using a shooting method and the Matlab bvp4c solver.For the last two methods we used a tolerance of 10 −6 .
We begin by transforming the domain 0, ∞ to −1, 1 , using the domain truncation method, the domain 0, ∞ is first approximated by the computational domain 0, L , where L is a fixed length that is taken to be larger than the thickness of the boundary layer.The domain 0, L is then transformed to −1, 1 using the algebraic mapping
The starting point is to assume that the independent variables f η , θ η , and φ η may be expanded as where the coefficient parameters

3.4
Starting from the initial approximations 3.5 which are chosen to satisfy the boundary conditions 2.7 , the subsequent solutions for F m , Θ m , and Φ m , m ≥ 1, are obtained by successively solving the linearized form of 3.3 which are

3.6
subject to the boundary conditions Once each solution for F i , Θ i , and Φ i i ≥ 1 has been found from iteratively solving 3.6 -3.7 for each i, the approximate solutions for f η , θ η , and φ η are obtained as In coming up with 3.8 , we have assumed that 3.9 Equations 3.6 -3.7 are integrated using the Chebyshev spectral collocation method Canuto et al. 33 and Trefethen 34 .The unknown functions are defined by the Chebyshev interpolating polynomials with the Gauss-Lobatto points defined as where N is the number of collocation points used.The unknown functions F i , Θ i , and Φ i are approximated at the collocation points by

3.11
where T k is the kth Chebyshev polynomial defined as The derivatives at the collocation points are represented as where a is the order of differentiation and D 2/L D with D being the Chebyshev spectral differentiation matrix.Substituting 3.13 in 3.6 -3.7 leads to the matrix equation

3.15
In the above definitions, a k,i−1 , and b k,i−1 , c k,i−1 k 1, . . ., 4 are diagonal matrices of size N 1 × N 1 .After modifying the matrix system 3.14 to incorporate the boundary conditions, the solution is obtained as 3.16

Spectral-Homotopy Analysis Method (SHAM)
The spectral-homotopy analysis method SHAM has been used by Motsa et al. 25,26 .It is also convenient to first ensure that the boundary conditions are made homogeneous by using the transformations where f 0 η , and θ 0 η , φ 0 η are chosen to satisfy the boundary conditions 2.7 of the governing equations 2.4 -2.6 .From 3.1 and the chain rule, we have that

3.18
Substituting 3.1 and 3.17 -3.18 in the governing equations and boundary conditions gives

3.19
where prime now denotes derivative with respect to ξ and The initial guesses used are

3.21
Solving the linear part of the equation system 3.19 , that is, will yield the initial SHAM approximate solution.Applying the Chebyshev pseudospectral method on equations 3.22 -3.23 yields the matrix form where

3.26
The superscript T denotes the transpose, diag is a diagonal matrix, and I is an identity matrix of size N 1 × N 1 .The boundary conditions 3.23 are implemented in matrix B and vector R of equation 3.24 .The values of Y 0 ξ 1 , Y 0 ξ 2 , . . ., Y 0 ξ N−1 are then determined from the following equation: 27 which provides us with the initial approximation for the solution of the governing equations 3.19 .With the initial approximate solution, we then find approximate solutions for the nonlinear equations 3.19 .We start by defining the following linear operators: where q ∈ 0, 1 is the embedding parameter and F ξ; q , Θ ξ; q , and Φ ξ; q are unknown functions.The zeroth-order deformation equations are given by

3.29
where is the nonzero convergence controlling auxiliary parameter and N F , N Θ , and N Φ are nonlinear operators given by

3.30
The m-th order deformation equations are given by

3.31
subject to the boundary conditions where

3.34
Applying the Chebyshev pseudospectral transformation to equations 3.31 -3.33 gives rise to the matrix equation subject to the boundary conditions where B and R are as defined in 3.25 and

3.37
Applying the boundary conditions 3.32 on the right-hand side of 3.35 yields the following recursive formula for higher-order approximations Y m ξ for m ≥ 1:

Improved Spectral-Homotopy Analysis Method (ISHAM)
Details of the improved spectral-homotopy analysis method ISHAM can be found in Makukula et al. 30 .The main objective is to improve the convergence rate of the spectralhomotopy analysis method by using an optimal initial approximation.Hence, instead of a random solution choice a systematic approach is employed to find the optimal initial approximation.This is achieved by first assuming that the solutions f η , θ η , and φ η can be expanded into

3.39
where F i , Θ i , and Φ i are unknown functions whose solutions are obtained using the SHAM approach at the ith iteration and F m , Θ m , and Φ m m ≥ 1 are known from previous iterations.We use the same initial guesses as with the SHAM solution in Sections 3.1 and 3.2.Substituting 3.39 into the governing equations gives

subject to the boundary conditions
and r 3,i−1 are as defined in equation 3.4 .Starting from the initial guesses 3.5 , the subsequent solutions F i , Θ i , and Φ i i ≥ 1 are obtained by recursively solving 3.40 using the SHAM approach.To find the SHAM solutions of 3.40 , we start by defining the following linear operators:

Mathematical Problems in Engineering 13
The zeroth-order deformation equations are given by

3.43
N F , N Θ , and N Φ are nonlinear operators given by

3.44
The mth order deformation equations are

subject to the boundary conditions
where

3.47
The initial approximations F i,0 , Θ i,0 , and Φ i,0 that are used in the higher-order equations 3.45 -3.47 are obtained by solving the linear part of 3.40 given by

3.49
In a similar manner, we apply the spectral methods to solve for the initial approximate solutions F i,0 , Θ i,0 , and Φ i,0 , and the higher-order deformation equations 3.45 -3.47 for higher order approximate solutions F i,m , Θ i,m , and Φ i,m for m ≥ 1.The solutions for F i , Θ i , and Φ i are then generated using the solutions for F i,m , Θ i,m , and Φ i,m as follows:

3.50
The i, m approximate solutions for f η , θ η , and φ η are then obtained by substituting F i , Θ i , and Φ i from 3.50 into 3.39 , where i is the ith iteration of the higher-order deformation equation and m is the mth iteration of the initial approximation.

Results and Discussion
Equations 2.4 -2.6 subject to boundary conditions 2.7 have been solved using three recent semi-numerical techniques as described above.To validate our results, we have compared the skin friction coefficient, the Nusselt number, and the Sherwood number with the theoretical results of Muhaimin et al. 9 .We have further compared our results with the full numerical solutions obtained using the shooting method and the Matlab bvp4c routine.The comparison is given in Tables 1-3.Tables 1-3 give values of the skin friction, heat transfer rate, and the mass transfer rate, respectively, for different porosity values.The convergence to the two numerical results of the SLM is achieved at the third order of approximation, at the sixth order for the SHAM, and at  1 shows an increase in the surface shear stress f 0 with an increase in the porosity parameter λ.The increase in the skin friction with the porosity may be accounted for by the fact that the velocity gradient increases with porosity Takhar et al. 35 .Tables 2 and 3 show an increase in the surface heat transfer rate −θ 0 and the mass transfer rate −θ 0 with the porosity parameter for large suction values s 3 , suggesting an increase in temperature and concentration gradients with increasing porosity.
Figure 1 serves two purposes: a to give sense of the accuracy of the improved spectral homotopy analysis ISHAM by means of a comparison between the numerical results and the second-order improved spectral-homotopy analysis results and b to demonstrate the effects of the suction parameter s and the Hartmann number M on the velocity profiles f η .
Firstly we observe an excellent agreement between the second-order ISHAM and the numerical bvp4c results for all parameter values used.Secondly we note that these results are qualitatively similar to those reported in Noor et al. 6 for the case of one-direction shrinking m 1 and show that increasing the suction parameter s and the Hartmann number M leads to an increase in the velocity.This in turn leads to a decrease in the boundary layer thickness as fluid is sucked out of the flow region.

Conclusions
We have successfully solved the nonlinear system of equations governing MHD boundary layer past a porous shrinking sheet with a chemical reaction and suction.We demonstrated three recent innovative methods, namely, the successive linearisation method SLM , the spectral-homotopy analysis method SHAM , and the improved spectral-homotopy analysis method ISHAM , and compared the performance of the three methods with regard to the speed of convergence of the solution the number of iterations required , computational efficiency, and the ease of application of the method.The results were compared with those obtained using the well-known shooting method and the Matlab bvp4c solver.We found that the ISHAM converged at second order.The magnitude of the parameter values used did not affect its performance under the same conditions with the SLM and SHAM.Nevertheless, the ISHAM does not come cheap in terms of the size of the code and computer time, taking about three times as long as the SLM to compute the same result and about double the time taken with the SHAM.The SLM converged at third order, is easy to implement, and has

Mathematical Problems in Engineering
shown a good level of stability when solving highly nonlinear problems.The SHAM gives good convergence under the same conditions but poor convergence with highly nonlinear problems.It is easy to implement but not as easy as with the SLM.Results from simulations revealed an excellent agreement between results from the shooting method and the bvp4c.Our findings indicate that the ISHAM is the best approach of the three methods in terms of the accuracy of the results and speed of convergence.Parametric studies for effects of different parameter values in the problems agreed with results present in the literature.

Figure 1 :
Figure 1: On the comparison between the 2nd-order ISHAM solution figures and the bvp4c numerical solution solid line for f η and θ η at different values of λ when M 1, m 1, Pr 3, λ 1 2, s 1, L 30, and N 150.

Table 1 :
Comparison of the approximate solutions of f 0 at different orders of the SLM, SHAM, and ISHAM against the numerical solutions at different values of λ when s 3, M 1, m 1, Sc 0.62, γ 3, Pr 1, λ 0, −1, L 30, and N 150.

Table 2 :
Comparison of the approximate solutions of −θ 0 at different orders of the SLM, SHAM, and ISHAM against the numerical solutions at different values of λ when s 3, M 1, m 1, Sc 0.62, γ 3,
second order for the ISHAM.Comparison with results reported in Muhaimin et al. 9 shows an excellent agreement.Table