On the Successive Linearisation Approach to the Flow of Reactive Third-Grade Liquid in a Channel with Isothermal Walls

1 School of Mathematical Sciences, University of KwaZulu-Natal, Private Bag X01, Scottsville, Pietermaritzburg 3209, South Africa 2 Institute for Advanced Research in Mathematical Modelling and Computations, Cape Peninsula University of Technology, P.O. Box 1906, Bellville 7535, South Africa 3 Department of Mathematics & Applied Mathematics, University of Venda, Private Bag X5050, Thohoyandou 0950, South Africa


Introduction
The rheological properties of many fluids used in industrial and engineering processes do exhibit non-Newtonian behaviour [1,2].Meanwhile, the study of heat transfer plays an important role during the handling and processing of non-Newtonian fluids [3][4][5].A complete thermodynamics analysis of the constitutive function for fluid of the differential type with the third-grade fluid being a special case has been performed by Fosdick and Rajagopal [6].Similar studies with respect to non-Newtonian fluid are also reported by Makinde [7][8][9].Moreover, the thermal boundary layer equations for non-Newtonian third-grade fluid constitute a nonlinear problem, and their solutions in space provide an insight into an inherently complex physical process in the system.In most cases, the nonlinear nature of the model equations precludes its exact solution.In recent time, several approximation techniques have been developed to tackle this problem [10][11][12][13], for example, the Adomian decomposition method, the variation iteration method, the improved finite differences method, the spectral method, and so forth.The ideas of developing new hybrids of numerical-analytical scheme to tackle nonlinear differential equations have experienced a revival.One such trend is the spectral homotopy analysis method that has recently been reported in [14,15] which is a hybrid between the standard homotopy analysis method [16] and the Chebyshev spectral collocation method [17][18][19].Other novel strategies involve using the Padé technique to improve the radius of convergence of the analytical methods of solution.Recent studies that make use of the Padé technique include the Hermite-Padé [7], the Homotopy-Padé [16,20,21], and the Hankel-Padé approaches.
The purpose of the present work is to present a new method, called the successive linearisation method (SLM), of solving nonlinear boundary value problems.We demonstrate the applicability of SLM in tackling nonlinear differentiation equations modeling the flow of a reactive third-grade liquid between two parallel isothermal plates.The mathematical formulation of the problem is established in Section 2. In Section 3 we introduce and apply some rudiments of SLM technique.Both numerical and graphical results are presented and discussed quantitatively with respect to various parameters embedded in the system in Section 4. A limited parametric study comparing numerical results generated using MAT-LAB's bvp4c boundary value solver is compared with the SLM results, and good agreement is observed.Using the SLM approach, multiple solutions which were theoretically proved to exist for such problems in [7] are also generated.

Mathematical Formulation
Figure 1 depicts the problem geometry.We consider the steady flow of an incompressible third-grade reactive fluid placed between two parallel isothermal plates.It is assumed that the flow is hydrodynamically and thermally fully developed under the action of a constant axial pressure gradient.
Following [1,4,5,[7][8][9], the dimensionless governing equations for the momentum and energy balance can be written as 2 where  is the dimensionless velocity component,  is the dimensionless normal coordinate,  is the dimensionless temperature, and , , , and  represent the Frank-Kamenetskii parameter, activation energy parameter, the dimensionless non-Newtonian parameter, and the viscous heating parameter, respectively.The additional Arrhenius kinetics term in energy balance equation ( 2) is due to [3].
The appropriate boundary conditions in dimensionless form are given as follows: the surface of the channel is fixed, impermeable, and maintained at a given temperature: and the symmetry condition along the centerline, that is, We have employed the following nondimensional quantities in (1)-( 4): where  is the absolute temperature,  is the fluid characteristic velocity,  0 is the plate temperature,  is the thermal conductivity of the material,  is the heat of reaction,  is the rate constant,  is the activation energy,  is the universal gas constant,  0 is the initial concentration of the reactant species,  is the channel half width,  3 is the material coefficient,  is the modified pressure, and  is the fluid dynamic Third-grade reactive fluid viscosity coefficient.In the following sections, ( 1)-( 4) are solved numerically using the new successive linearization technique.

Successive Linearisation Method (SLM)
In this section, we apply the proposed linearisation method of solution, hereinafter referred to as the successive linearisation method (SLM), to solve the governing equations ( 1) and ( 2).Before applying the SLM, we first note that the problem can be significantly simplified by finding the explicit analytical solution for the derivative /.First, ( 1) is rewritten as The above equation is then integrated on both sides, and the symmetry boundary condition /(0) = 0 is used to evaluate the resulting integrating constant to give We note that (7) is a cubic equation which can have either one or three real solutions.If only positive values of  are considered, (7) will have a unique real solution which can be computed using Maple and is given by where The analytical result is very important because, when evaluated at  = 1, it gives an explicit analytical expression for the skin friction coefficient (  ).A close inspection of (8) indicates that if   has a critical point, then (  ,   ) = (−3/2, −2/27).This critical point was reported as a bifurcation point in [7].It is worth noting that /, and hence the solution (), is only valid when  ≥ −2/27.
Since the momentum equation ( 1) is decoupled from the energy equation ( 2), we solve for the velocity () first then substitute the result in the energy equation to obtain ().To solve (), we write (1) as where () is a known explicit function of  given by with / given by (8).Equation (10) can easily be integrated using any numerical method.The SLM is based on the assumption that the unknown function () can be expanded as where   are unknown functions and   are successive approximations whose solutions are obtained recursively from solving the linear part of the equation that results from substituting (12) in the governing equations (2) using  0 () as an initial approximation.The linearisation technique is based on the assumption that   becomes increasingly smaller as  becomes larger, that is, lim Substituting ( 12) in (2) gives where () is a known function (from ( 8) and ( 9)) given by We choose  0 = 0 as an initial approximation which is chosen to satisfy the boundary conditions.The subsequent solutions for   ( ≥ 1) are obtained by successively solving the linearised form of ( 14) which are given as subject to the boundary conditions where Once each solution for   ( ≥ 1) has been found from successively solving (11) and for each , the approximate solutions for () are obtained as where  is the order of SLM approximation.We remark that the coefficient parameter  −1 and the right-hand side  −1 of ( 16) for  = 1, 2, 3, . . .are known (from previous iterations).Thus, (16) can easily be solved using analytical means (whenever possible) or any numerical methods such as finite differences, finite elements, Runge-Kuttabased shooting methods, or collocation methods.In this work, (10) and ( 16) are solved using the Chebyshev spectral collocation method.This method is based on approximating the unknown functions by the Chebyshev interpolating polynomials in such a way that they are collocated at the Gauss-Lobatto points defined as where  is the number of collocation points used (see, e.g., [17,19]).The derivatives are approximated at the collocation points by where D is the Chebyshev spectral differentiation matrix (see, e.g., [17,19]).Substituting ( 21) in ( 10) and ( 16) leads to matrix equations given by in which A is a ( + 1) × ( + 1) square matrix and T, W, F, and R are ( + 1) × 1 column vectors defined by W = [ ( 0 ) ,  ( 1 ) , . . .,  ( −1 ) ,  (  )]  , (27) In the above definitions, a −1 is a diagonal matrix of size (+1)×(+1) and the superscript  denotes transpose.After modifying the matrix system ( 22) and ( 23) to incorporate boundary conditions, the solutions for () and   () are obtained as (29)

Results
In this section, we present the results showing the velocity distribution and temperature distribution, for different values of the governing parameters.To check the accuracy of the proposed successive linearisation method (SLM), comparison is made with numerical solutions obtained using the MATLAB routine bvp4c, which is an adaptive Lobatto quadrature scheme.All the SLM results presented in this work were generated using  = 60 collocation points.
Figure 2 depicts the effect of non-Newtonian parameter () on both the velocity and temperature profiles.Generally, both fluid velocity and temperature profiles attained their maximum values along the channel centerline and minimum at the walls satisfying the boundary conditions.Moreover, a gradual decrease in the magnitude of fluid velocity and temperature profiles is noticed with an increase in the value of .This can be attributed to the fact that as  increases, the fluid viscosity increases leading to a decrease in the flow rate.In Figure 3, we observed that the fluid temperature generally increases with an increase in the value of the Frank-Kamenetskii parameter () due to the Arrhenius kinetics.The effect of viscous dissipation parameter () on the fluid temperature is displayed in  generation due to viscous heating increases as the parameter value of  increases leading to a general increase in the fluid temperature.Meanwhile, the possibility of a lower and upper solution branches is also highlighted in Figures 3 and 4.This can be attributed to the nonlinear nature of the Arrhenius kinetics in the governing thermal boundary layer equation (2).It is noteworthy that the fluid temperature decreases with an increase in the activation energy parameter () as illustrated in Figure 5.As  increases, the fluid becomes less volatile, and its activation energy decreases.A slice of the bifurcation diagram for  > 0 in the (, −  (1)) plane is shown in Figure 6.It represents the variation of wall heat flux with the Frank-Kamenetskii parameter ().In particular, for every 0 ≤  ≤ 0.1 there is a critical value   (a turning point) such that, for 0 ≤  <   there are two solutions.This result is in perfect agreement with the one reported in Makinde [7].
In Table 1, we show the computations illustrating the comparison between the SLM results at different orders and the bvp4c numerical results for wall heat flux Nu for various values of , , , and .It can be seen from the table that the SLM results are in very good agreement with the bvp4c numerical results.

Conclusion
In this work, we employed a very powerful new linearisation technique, known as the successive linearisation method  (SLM), to investigate the flow of reactive third-grade liquid in a channel with isothermal walls.The SLM results for the governing flow parameters were compared with results obtained using MATLAB's bvp4c function, and excellent agreement was observed.From the results obtained in the study, the following was observed.
(i) An increase in the non-Newtonian parameter () leads to a gradual decrease in the magnitude of fluid velocity and temperature profiles.
(ii) The fluid temperature generally increases with an increase in the value of the Frank-Kamenetskii parameter ().
(iii) The internal heat generation due to viscous heating increases as the parameter value  increases leading to a general increase in the fluid temperature.
It was also shown that the governing nonlinear equations admit multiple solutions.Using the SLM approach, lower and upper branch solutions were obtained and discussed.

Figure 1 :
Figure 1: Geometry of the problem.

Figure 4 .Figure 4 :
scheme.All the SLM results presented in this work were generated using  = 60 collocation points.Figure2depicts the effect of non-Newtonian parameter () on both the velocity and temperature profiles.Generally, both fluid velocity and temperature profiles attained their maximum values along the channel centerline and minimum at the walls satisfying the boundary conditions.Moreover, a gradual decrease in the magnitude of fluid velocity and temperature profiles is noticed with an increase in the value of .This can be attributed to the fact that as  increases, the fluid viscosity increases leading to a decrease in the flow rate.In Figure3, we observed that the fluid temperature generally increases with an increase in the value of the Frank-Kamenetskii parameter () due to the Arrhenius kinetics.The effect of viscous dissipation parameter () on the fluid temperature is displayed in Figure4.The internal heat

Table 1 :
Comparison between the SLM results at different orders and the bvp4c numerical results for wall heat flux Nu for various values of , , , and .