A Spectral Relaxation Approach for Unsteady Boundary-Layer Flow and Heat Transfer of a Nanofluid over a Permeable Stretching / Shrinking Sheet

This paper introduces two novel numerical algorithms for the efficient solution of coupled systems of nonlinear boundary value problems. The methods are benchmarked against existing methods by finding dual solutions of the highly nonlinear system of equations that model the flow of a viscoelastic liquid of Oldroyd-B type in a channel of infinite extent. The methods discussed here are the spectral relaxation method and spectral quasi-linearisation method. To verify the accuracy and efficiency of the proposed methods a comparative evaluation of the performance of the methods against established numerical techniques is given.


Introduction
Exact solutions to a wide class of problems in engineering and science are generally available only for a limited range of problems.For this reason the quest for new techniques and the improvement of existing techniques for finding solutions of nonlinear equations are an ongoing challenge in engineering and science.In addition to the classical numerical methods, such as those based on finite differences, finite elements, and finite volume techniques, there is currently a wide variety of methods for nonlinear equations, such as, among others, linearisation methods [1][2][3][4] and the transform methods of Fokas [5][6][7].Advances in decomposition and variational methods in recent years have added to the repertoire of available techniques for finding solutions of nonlinear BVPs.In recent years these methods have been further augmented by He's [8,9] and Liao's [10,11] homotopy based methods and their various variants such as the spectralhomotopy analysis of Motsa et al. [12,13] and Liao's [14] generalized boundary element method.
This paper introduces two novel techniques based on a combination of linearisation techniques and spectral methods and that allow for simple and straightforward integration of systems of ordinary differential equations on finite and infinite domains.We present an overview of these techniques and provide a comparative evaluation of the two methods against results in the literature.
A recent study by Bachok et al. [15] investigated the twocomponent convection in a viscoelastic liquid of Oldroyd-B type occupying a horizontal channel of infinite extent and depth .The flow is governed by the coupled system of equations where (), (), and () represent, respectively, the nondimensional stream function, temperature, and nanoparticle volume fraction, Pr is the Prandtl number, Le is the Lewis number, Nb is the Brownian motion number, and  and Nt are constant dimensionless parameters.
In the present study we revisit the solution of the system of nonlinear equations ( 1)-( 5) using the spectral relaxation method (Motsa [16]) and the spectral quasi-linearization method (SQLM).The objective of this study is to give a comparative analysis of the performance of the two techniques in finding solutions of coupled highly nonlinear problems in fluid mechanics.We determine, inter alia, the accuracy of each method and demonstrate how dual solutions can be obtained using the SQLM.In addition, we present a systematic approach of obtaining critical parameter values and multiple solutions of the governing equations.

The Spectral Relaxation Method
In this section we describe the development of the spectral relaxation method (SRM) for the solution of the nonlinear system (1)-(3).The system of equations is decoupled and the resulting subsystems are solved in a sequential manner.The method appears to be particularly effective for nonlinear systems of differential equations in which some of the unknown functions have exponentially decaying profiles.The algorithm for the method when applied to the system (1)-(3) may be summarized as follows.
(1) Reduce the order of the momentum equation from three to one by using the transformation (2) Assume that  is known from some initial guess and arrange the transformed equations in a particular order, placing the equations with the least unknowns at the top of the equations list.This gives with ( 2) and (3) placed below (8) in their original form.
(4) In the equation for  (1st equation), the iteration scheme is developed by assuming that only linear terms in  are to be evaluated at the current iteration level (denoted by  + 1) and all other terms (linear and nonlinear) in , , and  are assumed to be known from the previous iteration (denoted by ).In addition nonlinear terms in  are also evaluated at the previous iteration.
(5) In developing the iteration scheme for the next equation , only linear terms in  are evaluated at the current iteration level ( + 1) with all other terms evaluated at the previous level, except  which is now known from the solution of the first equation.
(6) This process is repeated in the equations for  using the updated solutions for , .The same procedure is effected on the equation for , now using the updated solutions for , , and .
The strategy used for decoupling the system of equations is analogous to the Gauss-Seidel relaxation method which is normally used in solving linear algebraic system of equations.In the context of the algorithm described above, we obtain the following iteration scheme: subject to the boundary conditions Given a set of initial approximations  0 ,  0 ,  0 ,  0 , (9) can be solved in a sequential manner, one after the other, using an appropriate numerical method.It must be noted that (9) now forms a sequence of linear differential equations with variable coefficients which can easily be solved using spectral collocation methods.The spectral methods are preferred here because of their remarkably high accuracy and relative simplicity in discretizing and subsequent solution of variable-coefficient linear differential equations with smooth solutions over simple domains.Spectral methods, such as the Chebyshev spectral collocation methods, have been found to be very efficient in discretizing and solving other iteration schemes for solving nonlinear equations of boundary layer type.Examples of recently developed spectral based method for solving boundary layer type equations include the spectral homotopy analysis method [12,13] and the successive linearisation method [17,18].The SLM has been successfully applied to different nonlinear boundary value problems in fluid mechanics; see Awad et al. [19], Motsa and Sibanda Advances in Mathematical Physics 3 [20], and Sibanda et al. [17,18].A large volume of literature exists on the practical implementation of spectral collocation methods including the books [21,22] which the interested reader may find useful.
In applying the spectral collocation method, we find the unknown function at  + 1 collocation points by requiring that (9) be satisfied exactly at these points.A convenient set of collocation points is the Gauss-Lobatto points defined on [−1, 1] by We approximate the derivatives of the unknown functions, say   (), using the so-called differentiation matrix  of size ( + 1) × ( + 1) which is computed as at the collocation points as the matrix vector product where  is a finite value used to numerically approximate the conditions at infinity, D = 2/, and f = [( 0 ), ( 1 ), . . ., (  )]  is the vector function at the collocation points.The variable  is defined from the linear transformation  = ( + 1)/2 used to convert the truncated interval [0, ] into the interval [−1, 1] on which the spectral method can be implemented.The second order derivatives are obtained as powers of D; that is Applying the Chebyshev spectral collocation method on (9) we obtain the following sequence of matrix equations: where where I is an ( + 1) × ( + 1) diagonal matrix, diag [ ] denotes a diagonal matrix, and f, g, Θ, and Φ correspond to the approximation of the functions , , , and  when evaluated at the collocation points.The SRM iteration scheme ( 14)-( 15) is solved for  = 0, 1, 2, . .., starting from suitable initial guesses  0 ,  0 ,  0 , and  0 .We choose the following initial approximations which satisfy the boundary conditions of the governing equations:

The Spectral Quasi-Linearisation Method
In this section, we present an alternative method for solving the governing system of ( 1)-( 3).We employ the quasilinearization method (QLM) which is a generalization of Newton-Raphson method.This method was initially proposed by Bellman and Kalaba [1] for solving nonlinear boundary value problems.Here we extend the application of the QLM to the governing equations ( 1)-( 3) and discretize the QLM equations using the spectral collocation method as described in the previous section.In this work, the spectral method based quasi-linearization method is referred to as the spectral quasi-linearisation method (SQLM).
To develop the SQLM scheme it is convenient to write the governing system as a sum of its linear L and nonlinear components  as for  = 1, 2, 3, where (1), (2), and (3) are labeled using  = 1, 2, and 3, respectively.Thus, we have The QLM is equivalent to multivariable Taylor series expansion of the nonlinear functions   assuming that the difference in the values of each unknown function at the current iteration denoted by, say,  +1 and the previous iteration   is small.Thus, the QLM scheme corresponding to (1)-( 3) is given by

= Nb𝜃
subject to To derive the above QLM scheme, it was noted that the solution for  can be resolved independently of ( 2) and (3).Thus, in deriving the iteration schemes for  and , the solution for  at the current iteration ( + 1) was assumed to be known.
Applying the Chebyshev pseudo-spectral method on (22) we obtain the following SQLM scheme in matrix form: subject to the boundary conditions where where O denotes an ( + 1) × 1 matrix of zeros and [ ]  denotes diag( ); that is, the vector elements are placed on the main diagonal of a matrix whose entries everywhere else are zero.Starting from the initial approximations (19), the approximate SQLM solutions for (), (), and () are obtained by solving (24) subject to the boundary conditions (25).

Determining System Critical Values
The system (1)-(3) was reported in Bachok et al. [15] to have multiple solutions and bifurcation curves on the planes (Re In this section, we demonstrate how the spectral quasilinearisation method (SQLM) may be used to determine the critical values of  beyond which solutions do not exist.
To determine the critical values   , we regard  and   (0) to be additional unknown variables.We define   (0) =  and differentiate (1) with respect to  to obtain subject to where ℎ = /.Criticality of the curve  = () corresponds to / = 0; thus, ℎ  (0) = 0. Applying quasilinearisation on ( 1) and ( 28)-( 29) gives subject to Equations ( 30) and (31) form a system of four equations which are solved iteratively for the unknowns , ℎ, ,  starting from suitable initial approximations  0 , ℎ 0 ,  0 ,  0 .Applying the spectral collocation method on (30) and (31) gives subject to the boundary conditions Advances in Mathematical Physics where

Results and Discussion
The governing systems of (1) and (3) were solved using the spectral relaxation method (SRM) and spectral quasilinearisation method (SQLM).In this section, we present the results of the numerical computations for the velocity, temperature, and concentration profiles for the various input parameters.The convergence and stability of the iteration schemes can be assessed by considering the norm of the difference in the values of functions between two successive iterations.For a converging iteration scheme, the difference in the norm of the value of the functions computed at successive iterations is expected to decrease with an increase in the number of iterations.The accuracy of the present results was verified by comparing with the results of Bachok et.al. [15].
In the calculations presented here, the values of the governing physical parameters were chosen deliberately to match the results of [15] in order to enable effective comparison.The effect of the number of collocation points  was examined in order to select the smallest value of  which gives a consistent solution.This is achieved by repeatedly solving the governing equations using the proposed iteration schemes with different values of  until the consistent solution within a small error tolerance level  was reached.In this work, the error tolerance was set to be  = 10 −9 .Unless otherwise specified, the number of collocation points used in this study is  = 100.The , that is,  at infinity, was also chosen in such a way that further changes in its magnitude do not produce changes in the values of  and its derivatives within the tolerance level .A value of  = 20 was found to be appropriate in all cases discussed in this investigation.Figures 1 and 2 show the variation of the error   , the difference of the norm of () computed using the spectral relaxation method at different iterations for different system parameter values.In all three cases it can be seen that   decreases as the number of iterations increases.This shows that the SRM converges for the selected choice of parameters.The convergence of the SRM improves with an increase in the unsteadiness parameter  and mass suction parameter  but decreases with an increase in the stretching parameter .
Figure 3 gives the variation of the SQLM generated error   with the number of iterations.The convergence of the SQLM is seen to be rapid for the first few iterations then plateaus after about five iterations.This result presents the main difference between the SRM and SQLM.The SRM results in Figures 1 and 2 indicate strict convergence of the method with the error progressively reduced with an increase in the number of iterations.For the same number of collocations points, the SQLM will converge faster than the SRM but is less accurate.This observation suggests that the convergence rate of the SQLM approach is superior.However, it turns out that the SRM yields a better accuracy for the three ODEs system.
To obtain the multiple solutions we choose an initial guess of the form where  is an unknown constant which must be carefully selected.It is observed that when varying  between negative and positive values, the concavity of the velocity profile   () changes in the region near  = 0.By fixing the values of the governing physical parameters and varying the values of  in the SQLM implementation, the values of  that give the multiple solutions can be identified from the maximum of the residual of the momentum equation (1). Figure 4 shows the variation of the maximum residual of the momentum equation ( 1) against  at different SQLM iterations when  =  =  = 1.It can be seen that the graph develops two local minima.The multiple solutions were obtained by setting  to be equal to the value of  at which the local minima occurs.For example, when  =  =  = 1, two solutions were obtained when  = −0.5 and  = 1.It is worth mentioning here that the SRM approach did not yield the second solution.
Figure 5 shows the two solution profiles () and   (), respectively, for different values of  and fixed values of  and .The corresponding temperature and concentration profiles are depicted in Figure 6.
In Table 1 we give a comparison between the present SQLM results and the results of [15] for the critical values of  for different values of  and .We note that the present results are in good agreement with the results reported in Bachok et al. [15].This shows the accuracy of the SQLM when compared to the shooting method as used in [15].
Figure 7 shows the variation of the local skin friction and local Nusselt and Sherwood numbers with  for different values of  and  and corresponding to the other parameters values reported in Bachok et al. [15].The SQLM easily resolves the dual solutions when  >   .These results obtained using the SQLM are qualitatively and quantitatively the same as those of Bachok et al. [15].The results show the two solutions branches for each value of the unsteadiness parameter  and the mass flux parameter .

Conclusion
In this study we have used two methods, the spectral relaxation method and the spectral quasi-linearisation method to solve the highly nonlinear equations that describe the unsteady heat transfer in a nanofluid over a permeable stretching or shrinking surface.The results obtained by Bachok et al. [15] have been confirmed.In terms of the accuracy, convergence rates, and robustness of the techniques used, we have demonstrated that (i) The SQLM converges faster than SRM; (ii) The SRM is more accurate than the SQLM;

Figure 1 :Figure 2 :
Figure 1: Effect of varying the parameters  and  on the error   from the SRM.

Figure 3 :
Figure 3: Effect of varying parameters on the SQLM error   .
The parameters   , Nu  , and Sh  are, respectively, skin friction coefficient, local Nusselt number, and the local Sherwood number and are defined in terms of the governing parameters as Nu  , ), and (Re −1/2  Sh  , ).