Tau – Homotopy Analysis Method for Solving Micropolar Flow due to a Linearly Stretching of Porous Sheet

A modification of the homotopy analysis method (HAM) for solving a system of nonlinear boundary value problems (BVPs) in semi–infinite domain, micropolar flow due to a linearly stretching of porous sheet, is proposed. This method is based on operational matrix of exponential Chebyshev functions to construct the derivative and product of the unknown function in the matrix form. In addition, by using Tau method the problem converts to algebraic equations to obtain the solution iteratively. In whole we can say, the computer oriented of this method is the most important aspect of it. During comparison between our methods and those previously reported, significant consequences are demonstrated.


Introduction
The homotopy analysis method (HAM) [21] is a powerful tool and has been already used for several nonlinear problems [1,2,3,4,5,12,13,14,15,22,28,30,34,35].The replacement of a nonlinear equation by a system of ordinary differential equations is the basic idea of the homotopy analysis method which prepare us to solve this system with the symbolic computation software such as Maple, Mathematica and Matlab.The solution of this system of ODEs is used to form a convergent series which, as proved in [21,23], is the solution of the original nonlinear equation.In using the HAM, and in order to

A brief discussion of problem
The dynamics of micropolar fluids, originated from the theory of Eringen [8], refers to the flow of a class of non-Newtonian fluids, can display the effects of local rotary inertia and couple stresses arising from practical micro-rotation.Physically they display many industrially important liquids including randomly oriented particles suspended in a viscous medium.The classical theories of continuum mechanics are inappropriate to expound the microscopic aspect of such complex hydrodynamic behavior.Theory of Eringen [8] is able to explain the experimentally displayed phenomenon of drag reduction in skin friction near the surface of a rigid body in fluids containing small amount of additives, when compared with the skin friction in the same fluids without additives.The Newtonian fluid mechanics cannot explain this phenomenon.This theory is also capable of explaining the flow of colloidal fluids, liquid crystals, animal blood, etc.A commonly noticed situation in industrial problems is a plate moving through a station-ary fluid.Valeggaar [33], Gupta and Gupta [11] have analyzed the stretching problem with constant surface temperature.Since the physical properties of the ambient fluid effectively influence the boundary layer characteristics, the investigation of non-Newtonian fluid over a moving sheet has attracted considerable attention.Fox et al. [10] have studied the flow of a power law fluid on a moving surface.The similarity solution of steady flow of micropolar fluid past a stretching surface, have investigated by Heraska and Watson [16].Also Hassanien and Gorla [17] analyzed the fluid mechanics and heat transfer characteristics of a steady, incompressible, micropolar fluid flowing past a non-isothermal stretching sheet with suction and blowing.The same problem was considered by Hady [18] Using a method of successive approximation.During another study, Flow over a porous stretching sheet with strong suction or injection was examined by Kelson and Farrell [20], they obtained analytical results for the shear stress and the microrotation at the surface for the limiting cases of large suction and injection.According to the theory, an additional transport equation representing the principle of conservation of local angular momentum must be solved, as well as the usual transport equations for the conservation of mass and momentum.Lukaszewicz provided last theories and an extensive survey of the literature [24].These equations involve a micro-rotation (angular velocity) field, in addition to the familiar velocity field.The two fields are coupled through a viscosity parameter.The limiting behavior, as the coupling parameter vanishes, is explored when the two fields are uncoupled.The velocity field is, then, governed by the Newtonian flow equations.This condition was analyzed by El-Mistikawy [9].

Problem formulation
Using cartesian coordinates, we consider an incompressible, steady micropolar fluid flowing over a flat sheet located at y = 0 and with a fixed end at x = 0. We suppose that the sheet is linearly stretched away from a fixed point o-as the origin of coordinates systemin x-direction.To investigate the flow, we begin with the dimensionless boundary-layer equations for steady self-similar flow of mentioned micropolar fluid in an otherwise quiescent medium.Then there are [9] In the above, a prime denotes differentiation with respect to η, the similarity coordinate measuring distances normal to the sheet.F (η) is the streamfunction, H(η) is the microrotation whose direction of rotation is normal to the plane of the flow, ε implies microrotational parameter, as 0 ≤ ε < 1 is in general likely to be function of the material properties of the working fluid and the flow conditions.Also β refers to a constant characteristic of the fluid.And satisfying boundary conditions are ) where ϕ > 0 corresponds to suction, and ϕ < 0 refers to injection.The present paper deals with the boundary layer flow of a micropolar fluid due to a linearly stretching flat sheet.By considering governing system of conservation equations in their self-similar form, we investigate the variation of sheer stress on the sheet by implementing THAM.Also velocity field and micro-rotation field are obtained.Finally during comparison with numerical results [6] efficiency of present method is obtained and excellent agreement between the numerical scheme is employed.
3 Properties of exponential Chebyshev functions

Exponential Chebyshev functions
The Chebyshev polynomials, denoted by T i (z) are orthogonal with Chebyshev weight function: where δ i,j is the Kronecker function, c 0 = π and c i = π/2 for i > 0.
The Chebyshev polynomials are generated from the three-term recurrence relation: where T 0 (z) = 1 and T 1 (z) = z.Other relations which we use here are In order to use these polynomials on the interval x ∈ [0, ∞) we define the exponential Chebyshev functions by introducing the change of variable z = 1 − 2 exp(−x/l).Let the exponential Chebyshev functions T i (1 − 2 exp(−x/l)) be denoted by ET i (x).The exponential Chebyshev functions are orthogonal with respect to the weight function ω(x) = 1/l √ exp(x/l) − 1 in the interval [0, ∞) with the orthogonality property: Then ET i (x) can be obtained as follows: A function f (x) defined over the interval [0, ∞) can be expanded as the below formula: where the coefficients a i are given by: In practice, only the first m-terms exponential Chebyshev functions are considered.Then we have

The derivative operational matrix of shifted Chebyshev
In first, we described some properties of the exponential Chebyshev functions which are needed here by using Eqs.(3.6) and (3.7) as The derivatives operator of the vector Φ(x) can be expressed by where D can be generated by using Eqs.(3.9), (3.10) and (3.11) as that can be considered in the following form B is the banded non-singular matrix and its inverse can be simply found.Therefore, Φ′ (x) is given as By adding ET ′ 0 (x) to Φ′ (x) and deleting ET ′ m (x) from Φ′ (x), Φ ′ (x) is constructed.At last, the derivative operational matrix of the exponential Chebyshev functions is obtained by adding the vector 0 in the first raw of 1/l B −1 C and deleting the end raw of it as As it can be seen, the derivative operational matrix of the exponential Chebyshev functions is lower-triangular, so there exists more zero elements in it.Consequently, It can be computed fast.
Remark 3.1.The derivative of order p for the vector Φ(x) can be expressed by Subsequently, the derivative of order p for the function f (x) is obtained as

The product operational matrix of the exponential Chebyshev functions
The following property of the product of two exponential Chebyshev functions vectors will also be applied.
where Ṽ is an m×m product operational matrix for the vector . Using above equation and by the orthogonal property equation (3.8) the elements where g ijk is given by To achieve g ijk more simply, we described the product of two the exponential Chebyshev functions, ET i (x)ET j (x) as Eq.(3.7) and then we have: Thus, the product operational matrix for the vector V is given by

The methodology
Prior to an outline of the Tau-homotopy analysis method (THAM) let us present a brief description of the standard homotopy analysis method.Readers may reffer to [23] for more general description of the HAM.

The homotopy analysis method
We consider the following differential equation where F is a nonlinear operator, x denotes the independent variables, f (x) is an unknown function, respectively.While f 0 (x) denotes an initial guess of the f (x), ̸ = 0 an auxiliary parameter, H(x) ̸ = 0 an auxiliary function, L an auxiliary linear operator, q ∈ [0, 1] as an embedding parameter and ψ(x; q) is an unknown function, Liao [23] constructs the so-called zero-order deformation equation, It should be emphasized that we have great freedom to choose the initial guess, the auxiliary linear operator L, the non-zero auxiliary parameter , and the auxiliary function H(x).Obviously, when q = 0 and q = 1, both ψ(x; 0) = f 0 (x) and ψ(x; 1) = f (x), are hold.Thus as q increases from 0 to 1, the solution ψ(x; q) varies from the initial guess f 0 (x) to the solution f (x).Expanding ψ(x; q) in Taylor series with respect to q, one has where If the initial guess, the auxiliary linear operator, the non-zero auxiliary parameter , and the auxiliary function H(x) are properly chosen, then the series (4.17) converges at q = 1 that by using ψ(x; 1) = f (x) one has the so-called homotopy series solution, which must be one of solutions of the original nonlinear equation, as proved by Liao [23].Define the vector Differentiating eq. ( 4.16) n times with respect to the embedding parameter q and then setting q = 0 and finally dividing them by n!, we have the so-called nth-order deformation equation where and According to the definition (4.20), the right-hand side of Eq. (4.19) is only dependent upon f n−1 (x).Thus, we gain f 1 (x), f 2 (x). . .by mean of solving the linear high-order deformation equation (4.19) one after the other in order.

Application of the THAM on micropolar flow
According to Eqs. (2.1), (2.2) and the boundary conditions (2.3) and (2.4), solution can be expressed in the form ) where A n and B n (n=1,2,...) are vectors of order m to be determined.According to the rule of solution expression denoted by (4.22) and (4.23) and the boundary conditions (2.3) and (2.4), it is natural to choose Now from the Eq.(4.19) we have Against the standard HAM, there is no limitation in THAM to choose the auxiliary functions H 1 (η) and H 2 (η).Thus, we choose H 1 (η) = H 2 (η) = 1 for this problem.By using the operational matrices obtained in pervious section, we can define R 1,n (η) and R 2,n (η) as where Λ = D T A. Therefore, Eqs.(4.24) and (4.25) is simply obtained as where L 1 and L 2 are the matrices form of auxiliary linear operators L 1 and L 2 which are defined as Similar to the typical Tau method [29,7], one has subject to the boundary conditions 2 , where L 1 and L 2 are the modified matrices L 1 , L 2 , and Q 1 and Q 2 are modified vectors Q 1 and Q 2 , respectively, after implementing boundary conditions.Now by starting from A T 0 and B T 0 as initial approximation the above equations are ripe to be iteratively solved to obtain {A 1 , A 2 , ..., A l } and {B 1 , B 2 , ..., B l }, and consequently, it gives

Error analysis
We can check the accuracy of the method by defining the residual functions by using Eqs.
(2.1), (2.2), (4.30) and (4.31) as (5.33) Now by defining the norm of ∥ • ∥ ω(η) and knowing that Res 1 and Res 2 ∈ E 2m , where } , the norm of the residual functions can be determined by the well-known Gauss quadrature formulation as [31] ∥Res where η i are zero of ET m+1 (η) and obtained as It's clear that the minimum of ∥Res j (η, )∥ ω(x) might be determined conveniently by the auxiliary parameter .suppose that s j is an appropriate guessed value.Thus, we can obtain the valid interval of convergence that determines the optimal values of .The appropriate region of to get convergent solution of F ′′ (0) and H ′ (0) with ε = 0.1, ϕ = 0 and various values β is shown in Figs. 1 and 2. The residual errors in Eq. (5.34) is displayed in Fig. 3 which indicates the interval for the admissible values of .6 Concluding remarks

Results
A full description of the flow field could be obtained by considering suction, injection, varying coupling parameter and various values of β.For the purpose of discussing the results, the numerical calculations are presented in the form of F ′′ (0) and H ′ (0) representing the sheer stress at the surface and gradient of microrotation, respectively.The tabulated data in tables 1-6 display a comparison between results of and present methods.During the comparison, stability, high accuracy and efficiency of the applied approaches are found.

Present methods ϕ
Numerical [6] HAM [6] Shooting method THAM 1.0 A comparison between THAM results with m = 5 and r = 10 and those obtained by [6], of F ′′ (0) for various ε when ϕ = 0 and β = 6.A comparison between THAM results with m = 5 and r = 10 and those obtained by [6], of H ′ (0) for various ε when ϕ = 0 and β = 6.Also it can be clearly seen that the suction (ϕ > 0) increases the surface sheer stress as well as wall gradient of microrotation, whereas the injection (ϕ < 0) has the opposite effect.The influence of varying coupling parameter ε on F ′′ (0) and H ′ (0) is observed in tables 5 and 6.

Present methods
The effect of different values of β on microrotation field are illustrated in Figures 4  and 5. From these figures it can be clearly seen that microrotation increases up to a certain distance from the surface, reach a maximum value, then decreases gradually and approaches to zero.Figures 6 and 7 illustrate solutions for the fluid velocity profile and microrotation profile in the case of suction ϕ > 0 and injection ϕ < 0 when β = 6 and ε = 0.1.Figure 7 indicates that a rise in ϕ from −1 through −0.5, 0, 0.5 to 1 induces a decline in F ′ values a growth in boundary layer thickness.Figure 6 presents the effect of ϕ on microrotation.From the figure it is clear that microrotation field first increases near the boundary as ϕ increases, but away from the boundary a reverse pattern is observed.
Figures 9 and 8 display the behavior of fluid velocity and microrotation field for different values of coupling parameter.Figure 9 shows the variation of fluid velocity when continuously decreases with an increase in ε value from 0.1 to 0.9.In figure 8 we illustrated the profile of microrotation which decreases near the surface as ε increases.Additionally, by increasing the value of coupling parameter the rate of decreasing to zero decreases.

Conclusion
In this paper, we have considered a novel Tau-homotopy analysis method (THAM) to solve self-similar flow of a micropolar fluid driven by a linearly stretching sheet included suction and injection.In this way, by constructing the exponential Chebychev functions in [0, ∞), the operational matrices of derivative and product are obtained.The following advantages that make our proposed method distinct are noticeable: 1.In Tau-homotopy analysis method (THAM), the rule of solution expression and the rule of ergoticity are not workable unlike the standard HAM approach.
2. In using THAM one may use any form of initial guess (even complicated ones) as long as it satisfies the boundary conditions whereas in HAM one is restricted to choose an initial approximation that will make the integration of the higher order deformation equations possible.
3. Using the operational matrices and Tau method change the standard homotopy analysis method to an iterative matrix method which is capable of greatly reducing the size of calculation while maintaining the accuracy of the solution.Actually our described method is not computationally expensive.

. 29 )
The boundary conditions (4.28) and (4.29) are imposed on two last columns of the L 1 and L 2 matrices on left hand side of the Eqs.(4.26) and (4.27) and their values replaced in the right hand side vector of the mentioned equation.Obviously, it gives