1 Introduction

Thin film fluid flow has become dependent on many theoretical and experimental studies in recent years due to its widespread applications in industry and engineering. These applications can be found in crystal production, continuous casting, wire and fiber coating, chemical processing equipment, and tinning of copper wires. Many authors have studied the thin film fluid flow and heat transfer under different cases^[1-6]. After these previous publications, a number of researchers have successfully applied numerical methods in this field. For example, the multiple shooting subroutine method was introduced to get the numerical solution for laminar liquid film moving on a horizontal stretching sheet by Andersson et al.^[7]. The method of adjoins was presented to evaluate the numerical solution for the fluid flow and heat transfer through a liquid film over a horizontal stretching sheet under different conditions by Dandapat et al.^[8-9]. In the same context, Dandapat and Maity^[10] obtained an analytic expression for the film thickness which governs the thin film flow past an unsteady stretching sheet by using the singular perturbation technique and the method of characteristics. Likewise, Noor and Hashim^[11] used the homotopy analysis method to get the numerical solution for the influence of thermo-capillarity and a magnetic field on the flow of a liquid film due to an unsteady stretching sheet. Among these numerical methods is the Legendre collocation method (LCM) which is a general approximate analytical method to solve the differential equations. The LCM has some advantages for handling this class of problems in which the Legendre coefficients for the solution can be obtained easily after using the numerical programs. For this reason, this process is much faster than other methods. Legendre polynomials have been widely used because of their good properties in the approximation of functions^[12-16]. This technique provides parity in functions that meet the exact solution of the problem, where the problem is solved without the need for discretization of variables. Therefore, in some problems, the round of errors is not affected by one error that you do not experience with large computer memory. The proposed scheme provides the solution to the problem in a closed form, but for example, the finite difference method provides the approximation at the grid points only. Orthogonal polynomials have a great variety and wealth of properties. Some of these properties take a very concise form in the case of the Legendre polynomials, making Legendre polynomials of leading importance among orthogonal polynomials^[13]. Also, the advantage of using Legendre polynomials as a tool for expansion functions is the good representation of smooth functions by finite Legendre expansion provided that the function y(x) is infinitely differentiable. The LCM has been used to solve many problems^[14-16].

In this work, we will use the properties of the Legendre polynomials to derive an approximate formula of the integer derivative D⁽ⁿ⁾y(x) and estimate an error upper bound of this formula, and use this formula to solve numerically the proposed problem. Due to high accuracy of this method, it is inevitable to use it for solving the nonlinear system of ordinary differential equations which describes flow and heat transfer of thin liquid film which is affected by the presence of thermal radiation and magnetic field.

2 Formulation of the problem

The basic ordinary differential equations, which describe the thin film liquid layer flow and heat transfer, can be summarized as^[3]

(1)

(2)

(3)

(4)

where the prime denotes differentiation with respect to η, S is the unsteadiness parameter, M is the magnetic parameter, R is the radiation parameter, δ is the slip velocity parameter, γ is the dimensionless film thickness, and Pr is the Prandtl number. Here, we note that this system of equations (1)-(4) is a generalization of the pioneering research of Wang^[3] and this system can be reduced to Wang's problem when we take M=R=δ=0.

3 Approximate formula for the derivatives of Legendre polynomial expansion

The well-known Legendre polynomials are defined on the interval [-1, 1] and can be determined with the aid of the following recurrence formula^[13]:

where L₀(z)=1 and L₁(z)=z. In order to use these polynomials on the interval [0, 1], we define the so-called shifted Legendre polynomials by introducing the change of variable z=(2η-1). Let the shifted Legendre polynomials L_k(2η-1) be denoted by . Then, can be obtained as follows:

where =1, and =2η-1. The analytic form of the shifted Legendre polynomials of degree k is given by

(5)

Note that (0)=(-1)^k, and (1)=1. The orthogonality condition is

The function u(η)∈ L₂[0, 1] may be expressed in terms of the shifted Legendre polynomials as

(6)

where the coefficients u_k (k=0, 1, ...) are given by

(7)

In practice, we use only the first (m+1)-terms of the shifted Legendre polynomials as follows:

(8)

The main approximate formula of the derivative of u_m(η) is given in the following theorem.

Theorem 1  Let u(η) be approximated by the shifted Legendre polynomials as in (8) and suppose that n is an integer number. Then,

(9)

Proof  The proof of this theorem can be done directly with the help of the formula (5) and some properties of the shifted Legendre polynomials.

In this section, special attention is given to study the convergence analysis and evaluate an upper bound of the error for the proposed approximate formula.

Theorem 2^[16]  The derivative of integer order n for the shifted Legendre polynomials can be expressed as a linear combination of the shifted Legendre polynomials in the following form:

(10)

where

Lemma 1^[16]  The error in approximating u⁽ⁿ⁾(η) by u_m⁽ⁿ⁾(η) is bounded by

4 Procedure solution

In this section, we will present the proposed method to numerically solve the system of ordinary differential equations of the form (1)-(2) coupled with the boundary conditions (3)-(4). The unknown functions f(η) and θ(η) may be expanded by finite series of the shifted Legendre polynomials as follows:

(11)

From Eqs. (1)-(2) and (11) and the formula (9), we can obtain

(12)

(13)

We now collocate Eqs. (12)-(13) at (m-n+1) points η_s (s=0, 1, ..., m-n) as

(14)

(15)

For suitable collocation points, we use the roots of the shifted Legendre polynomial . Also, by substituting (11) into the boundary conditions (3)-(4), we can obtain six equations as follows:

(16)

where , and . Equations (14)-(15), together with six equations of the boundary conditions (16) give an algebraic system of (2m+2) equations which can be solved for the unknowns f_i and θ_i (i=0, 1, ..., m), using the Newton iteration method. In our numerical study, we take m=5, i.e., five terms of the truncated series solution (11) at η=1.

5 Results and discussion

Since the problems which arise in industrial and technological situations are highly non-linear, it is difficult to get pure exact solutions to this type of problems. One of these problems occurs in the tinning of copper wires, polymer industry, and crystal growth involving the thin flow of a Newtonian fluid past a stretching sheet. A relatively new method called the shifted LCM which is beneficial for developing solutions that are valid even for moderate to large values of parameters. This method gives the solution in a series form. One of the chief factors that influences the convergence of the solution series is the type of the base functions used to express the solution. After applying the shifted LCM, the solution for the non-linear system of ordinary differential equations (1)-(4) is presented through Figs. 1-9.

Firstly in this section, we mentioned in our problem that the system under (1)-(4) can be reduced to Wang's problem when we take M=R=δ=0. We must report that there exist some quantities which are very useful in engineering industry. These physical quantities are called the skin-friction coefficient and the Nusselt number, respectively. The first quantity is proportional to the second derivative for the stream function. In other words, it is proportional to -f"(0)^[3]. The second quantity is proportional to the first derivative of the dimensionless temperature, i.e., it is proportional to θ'(0)^[3]. Therefore, to validate our numerical solution, we make the following comparison which proves thoroughness and trust for our proposed method. Our comparison is illustrated in Table 1. An excellent agreement is clearly noted between our data and the previous published work.

Figure 1 highlights the effect of S on the dimensionless velocity f'(η). These plots reveal the fact that increasing the parameter S results in enhancing the velocity distribution along the thin film region, but a reverse is observed for the film thickness. The dimensionless temperature for different values of the unsteadiness parameter S is displayed in Fig. 2. It is noticed that both the free surface temperature θ(γ) and the dimensionless temperature increase with the increase in S.

The influence of the magnetic parameter M on the dimensionless velocity is presented in Fig. 3. It is observed that, along the sheet, the dimensionless velocity increases with an increase in the parameter M, but a reverse trend is noted away from the sheet. Also, it is further found that the increasing value of parameter M is to decrease the film thickness.

In Fig. 4, we depict the effects of the same parameter M on the dimensionless temperature. It is noticed that, with the increasing values of the magnetic parameter M, both the free surface temperature θ(γ) and the dimensionless temperature increase inside the film region.

Figure 5 represents the variation of velocity profiles for different values of the slip velocity parameter δ. From this plot, it is evident that, along the sheet, the effect of increasing values of the parameter δ is responsible for thinning the film thickness and decreasing the dimensionless velocity.

Figure 6 displays the temperature θ(η) profiles versus η for various values of the parameter δ. It elucidates that both the dimensionless temperature distribution and the free surface temperature θ(γ) increase with an increase in the parameter δ.

Figure 7 depicts the effects of the radiation parameter R on the temperature profile. It is interesting to note that the dimensionless temperature distribution increases when the radiation parameter R increases. Likewise, some qualitative behaviors for the thin film flow and heat transfer characteristics are demonstrated in the same figure in which the thickness for the film is fixed for different values of the radiation parameter R.

Finally, at this step, we can illustrate the effect of the magnetic parameter M on both the skin-friction coefficient -f"(0) and the Nusselt number in Fig. 8. It is very clear that increasing the magnetic parameter M results in increasing the skin-friction coefficient. However, from the same figure, it is found that the Nusselt number -θ'(0) decreases with the increasing values of the magnetic parameter M.

At the same context, both the skin-friction coefficient -f"(0) and the Nusselt number -θ'(0) decrease with enhancing the values of the slip parameter δ as illustrated from Fig. 9.

6 Conclusions

The proposed numerical study highlights the lineaments of the radiation parameter, velocity slip parameter, magnetic number, and unsteadiness parameter in flow and heat transfer for a Newtonian thin liquid film. The LCM is used here to investigate the approximate solution of the resulting non-linear system of ordinary differential equations for the proposed problem. Numerical evaluations of the closed form results are performed, and graphical results are obtained to illustrate the details of the kthin film flow and heat transfer characteristics and their dependence on some physical parameters. On the basis of the obtained results, we can confirm the following observations:

(ⅰ) Increasing both the values of the unsteadiness parameter and the magnetic parameter causes an increase in both the dimensionless velocity and the dimensionless temperature throughout the film layer.

(ⅱ) The temperature distribution can be affected by changing the values of the velocity slip parameter and the thermal radiation parameter.

(ⅲ) The presence of slip velocity parameter may result in a lower velocity distribution near the stretching sheet and also thin the film thickness.

(ⅳ) The illustrated numerical results show the reliability and the efficiency of the suggested approach.

(ⅴ) The results obtained from the LCM are compared with those obtained from the homotopy analysis method in the special case of the given problem.

(ⅵ) It is seen from the presented results in tables and figures that the results obtained from the suggested method are more precise than those obtained from the existing methods.

(ⅶ) The results confirm the ability of the LCM to solve the problem under study.