1 Introduction

Recently,the study of flow and heat transfer of a finite thin film has attracted many researchers' attention because of its wide applications such as coating,polymer,metal extrusion,and drawing of plastic sheets. Wang^[1] first explored the motion behaviors of an unsteady thin finite Newtonian liquid film formed on an accelerating stretched sheet. Following his pioneering work,Andersson et al.^[2] and Liu and Andersson^[3] solved the nonlinear problem by using the multiple shooting subroutine method. The exact similarity solutions were obtained by using the homotopy analysis method (HAM)^[4] and the differential transform method (DTM)^[5]. A general surface temperature distribution of the wall-surface was taken into account,and feasible similarity transforms were obtained in Refs. ^[3] and ^[6]-[7]. Some of important references in the vast literature regarding to the thin finite liquid film can be found in Refs. ^[8]-[14].

The Newtonian fluid is important,and the fluid can display the non-Newtonian behaviors. Andersson et al.^[15] investigated the unsteady flow of a finite thin film within a non-Newtonian fluid following the power-law Ostwald-de Waele model. Chen^[16-18] presented heat transfer of power-law liquids in a finite thin film over an unsteady surface. The effect of viscous dissipation was considered by Chen^[17],and the Marangoni boundary condition was taken into account in Ref. ^[18]. Wang and Pop^[19] focused on the non-Newtonian behaviors of a power-law film and obtained the analytical expression of the critical value for the unsteadiness parameter. In addition,the effects of different constitutive models for non-Newtonian fluids,i.e.,the second grade fluid^[20],the biviscosity fluid^[21],and the power-law fluid^[22],were considered. Huang et al.^[23] solved the thin non-Newtonian power-law liquid film problem by using the Chebyshev finite difference technique. Vajravelu et al.^[24-25] also obtained a set of effective solutions by using the Keller-box technique. Recently,the topic of finite thin non-Newtonian fluid films was extended to pseudo-plastic nanoliquid films by Lin et al.^[26-27]. It should be noted that the base fluids,i.e.,sodium carboxymethyl cellulose (CMC) water,are non-Newtonian pseudo-plastic fluids in Refs. ^[26] and ^[27].

indent The study of the non-Newtonian rheological behaviors becomes more and more important because most liquids,such as pulps,varnishes,and multi-grade oils,do not obey the Newton inner friction law,and the molecular structures are different from the Newtonian fluids. Pop et al.^[28] and Gorla et al.^[29] developed a new constitutive model that the thermal diffusion coefficient of power-law liquids is a power-law function of the velocity gradient. The non-Newtonian heat transfer behaviors of power-law fluids were further explored in Refs. ^[30]-[36].

The aim of this research is to explore the influence of the internal heating source on momentum and heat transfer characteristic on a thin finite power-law liquid film,driven by an accelerating unsteady horizontal stretched sheet. Pop's thermal conductivity model and a general temperature distribution are taken into account. Numerical solutions to the similarity two-point boundary value differential equations are obtained by the shooting method coupled with the BVP4C,the Runge-Kutta method,and Newton's scheme.

2 Problem formulation and governing equations

Here,we consider an unsteady incompressible laminar flow and heat transfer within a finite thin non-Newtonian liquid film formed on a horizontal accelerating surface. The liquid film is issued from a narrow slot,and the thin film is smooth. The free interface wave,the shear viscous drag,and the thermal flux are in the adiabatic boundary interface. The schematic of the physical problem is shown in Fig. 1. Based on the above assumptions and the continuous medium hypothesis,the time-dependent governing partial differential equations (PDEs) for the mass conservation,momentum conservation,and thermal energy conservation can be written as

(1)

(2)

(3)

The corresponding boundary conditions can be expressed as

(4)

(5)

where $t$ is the time. $u$ and $v$ are the liquid velocity components in the $x$- and $y$-directions,respectively. $\tau_{xy}=K| \frac{\partial u}{\partial y}|^{n-1}\frac{\partial u}{\partial y}$ is the modified shear viscous drag of the power-law liquids,and $K$ is the viscosity coefficient,which is a positive constant. $n$ is the power-law index. The case 0<n< 1 represents the pseudo-plastic fluid,the case $n=1$ represents the Newtonian fluid,and the case $n>1$ represents the dilatant fluid. $T$ is the temperature of the liquids film,$T_0$ is the temperature at the origin of the narrow slot,$C_p$ is the specific heat capacity,and $\rho$ is the density of power-law fluids. In this research,the effects of a power-law velocity gradient and a viscosity coefficient on the thermal conductivity of power-law fluids are taken into account based on a new constitutive model proposed in Refs. ^[28]-[29],i.e.,$k=\omega|\frac{\partial u}{\partial y}|^{n-1}$,which $k$ is the effective thermal conductivity of non-Newtonian liquids,and $\omega$ is the consistency thermal coefficient. $Q(t)$ is the thermal absorption (< 0) or the thermal generation ($>0$),and it is proposed to have the form of ^[7]

(6)

where $\phi$ is the heating source parameter. $\phi< 0$ is for the thermal absorption,and $\phi ＞0$ is for the thermal generation. $u_{\mathrm s}$ is the horizontal velocity of the stretched surface of the thin finite film for the continuously movement on the wall interface,and has the form of

(7)

where $a$ and $b$ are both positive constants with the same unit. Furthermore,the wall interface temperature of the stretched surface $T_{\mathrm s}$ is supposed to have a generalized power-law form of the time $t$ and the horizontal coordinate $x$ as

(8)

where $T_{\mathrm {ref}}$ is a constant,i.e.,it is regarded as a standard temperature,$d$ is a positive constant,and $r_1$ and $r_2$ are the power indices. Equation (8) is restricted to $r_1+2\ge0$ and $r_2\ge0$,and the relational expressions given by Eqs. (7) and (8) are limited to the condition $t ＜a^{-1}$.

To proceed,we introduce the following similarity transformation ($\psi$ satisfies $u=\frac{\partial \psi}{\partial y}$ and $v=-\frac{\partial \psi}{\partial x}$):

(9)

(10)

(11)

(12)

(13)

(14)

Using these new variables,Eqs. (1) -(3) and (6) -(7) can be written as

(15)

(16)

(17)

(18)

where $Pr$ represents the generalized Prandtl number,$S$ represents the unsteadiness parameter,$\beta$ represents the dimensionless liquid film thickness,and $Re_x$ represents the local Reynolds number.

The liquid film velocity component and the shear viscous drag $\tau_{xy}$ are

(19)

(20)

(21)

For practical applications,a quantity of major interest is the local Nusselt number $Nu_x$,which can be expressed as $Nu_x=xq_{\rm w}(x)k^{-1}T_{\mathrm {ref}}^{-1}$,where $q_{\mathrm w}(x)=-k \frac{\partial T}{\partial y}|_{y=0}$ is the thermal flux from the wall surface of the thin liquid film,and $k=\omega|\frac{\partial u}{\partial y}|^{n-1}$. $Nu_x$ can be obtained as

(22)

3 Numerical methods

In order to obtain effective solutions for Eqs. (15) -(18) ,we convert Eqs. (15) -(18) into the following equations by denoting $f$,$f'$,$f"$,$\theta$,and $\theta'$ using variables $f_1$,$f_2$,$f_3$,$\theta_1$,and $\theta_2$,respectively,

(23)

(24)

(25)

(26)

(27)

There are five first-order ordinary differential equations (ODEs) in the governing equations (23) -(25) ,while there are six relational expressions in the boundary value conditions (26) -(27) . Therefore,there exists a relationship between the parameters $\beta$ and $S$. The above boundary value problem (23) with the condition (26) is solved by using the shooting technique coupled with the Newtonian scheme for a given value of $S$. The initial guessed value of $\beta$ is given with the program BVP4C in MATLAB. And $\beta$ is adjusted so that the iteration of the loop program meets the condition $f_1(\beta)=(2-n)\beta S/(2n)$. This is done on the basis of a trial-error method. Furthermore,the effective numerical solutions of the ODEs (23) -(25) and (26) -(27) are achieved by using the shooting technique combined with Newton's scheme and the standard fourth-order Runge-Kutta iterative program for given values of $S$ and $\beta$.

The relationship between the parameters $\beta$ and $S$ is solved by using the program BVP4C in MATLAB. For example,when $n=0.8$,Eqs. (23) and (26) are reduced to

(28)

(29)

Equations (28) -(29) are a two-point boundary value problem with an unknown parameter $S$. We solve it by using the program BVP4C in MATLAB and get $S=0.918 9$ when $\beta=1$.

Let $\beta$ change as $0.01\rightarrow1.00$ and $N$ be fixed,i.e.,$N=100$. Therefore,we can get the relationship between the parameters $\beta$ and $S$ when $0.01\le \beta\le 1.00$. The initial value of $S$ is $S=1.3$,and the initial solutions of Eqs. (28) -(29) are (when $0.01\le \beta\le 1.00$)

(30)

where $i=\beta^{-1}(2.25S-3) $ and $j=\beta^{-2}(3-2.25S)$. In the same way,we can obtain the relationship between $\beta$ and $S$ when $\beta\ge1.00$.

4 Results and discussion

For the hydrodynamic thin liquid film,there is a critical value of $S_0$ for the unsteadiness parameter,where no solution could be achieved^[16-19]. Results for the critical value are compared with those obtained from Refs. ^[16]-[19] in Table 1. It can be seen that the results of this research agree well with the other related results for different values of $n$.

Figure 2 illustrates the influence of the power-law index $n$ on the liquid film thickness $\beta$ and the unsteadiness parameter $S$. The $\beta-S$ diagram displays that the liquid film thickness $\beta$ decreases as the power-law index $n$ decreases at a specified value of the unsteadiness parameter $S$. In general,for the fixed values of the power-law index $n$ and other involved parameters,the dimensionless thickness decreases monotonically with the increase of $S$ from zero to the critical value. The dimensionless thickness first decreases very quickly with the increase of the unsteadiness parameter $S$,i.e.,the $\beta-S$ profile is close to the $y$-axis,then decreases more gradually,i.e.,the $\beta-S$ profile in the middle zone is closely paralleled to the $x$-axis,and finally drops to be zero rapidly as the unsteadiness parameter $S$ comes close to $S_0$ (these results are similar to those of Ref. ^[19]).

Figures 3 and 4 illustrate the influence of the unsteadiness parameter $S$ on the horizontal velocity component and temperature distribution. For particular values of $n$ ($n=0.8$) and other parameters ($Pr=2.0$,$r_1=0$,$r_2=1$,and $\phi=0$),the liquid film thickness $\beta$ decreases while the dimensionless horizontal wall interface velocity and temperature profiles increase as the unsteadiness parameter increases. For example,Fig. 3 displays that the horizontal velocity component at $n=0.8$ varies by $14.33%$ through the non-Newtonian liquid film for $S=1.2$,while by as much as $55.00%$ for $S=0.8$. Figure 4 shows that the temperature varies by $33.27%$ across the film for $S=1.2$ and by as much as $81.52%$ for $S=0.8$.

Figures 5-7 illustrate the effects of the generalized Prandtl number and the power indices ($r_1$ or $r_2$) of the wall surface temperature on the dimensionless temperature. The distribution curves display the dimensionless temperature decreases when $\eta$ changes from zero to $\beta$,i.e.,from the wall interface to the free interface of the thin liquid film,for both the power indices and the generalized Prandtl number. It should be emphasized that the dimensionless temperature also decreases as the generalized Prandtl number and the two power indices increase. Furthermore,the temperature decreases from one to zero very quickly as $Pr\to +\infty$ or $r_1\rightarrow +\infty$,$r_2\rightarrow +\infty$,and $\phi\to +\infty$. The generalized Prandtl number and the two power indices of the power-law fluid film have a significant effect on the results of the temperature distribution curves. In addition,a temperature distribution of $\theta(\eta)\approx1.00$ or $T=T_{\mathrm s}$,i.e.,the temperature of the fluid field is a constant,is obtained to prevail in the thin finite liquid film as $Pr\to 0$.

Figure 8 illustrates the effects of the heating source parameter $\phi$ on the dimensionless temperature distribution. In general,the thermal absorption ($\phi< 0$) has a tendency to cool down the liquid film,while the thermal generation ($\phi>0$) has a tendency to warm it up. The present research shows that the dimensionless temperature distribution curve declines as the heating source parameter ($\phi$,both including the positive case and the negative case) decreases,and the heating source parameter has a significant effect on the results of the temperature distribution.

5 Conclusions

In this article,we present the research for flow and heat transfer of the non-Newtonian power-law fluids within a finite thin liquid film over an unsteady stretching sheet. The Pop's thermal conductivity model and a general surface temperature distribution are taken into account. Some of the interesting results are listed as follows:

(a) The liquid film thickness increases while the free-surface horizontal velocity component and temperature decrease as the unsteadiness parameter declines.

(b) The temperature distribution curves go down as the generalized Prandtl number and the power indices increase,while go up as the heating source parameter increases.

Acknowledgements The authors want to thank the anonymous referees for their valuable sugges-tions and comments.