Effects of second-order slip and drag reduction in boundary layer flows

As we know, a moving flat plate in a fluid medium infuses a boundary layer. This kind of flow appears in several technological industries, such as extrusion process, wire and fiber coating, polymer processing, food-stuff processing, design of heat exchangers, and chemical processing equipment. The pioneering work of Sakiadis [1, 2] on the laminar boundary layer over a rigid surface moving in its own plane is quite different from the flow past a stationary surface (known as the classical Blasius [3] flow). Tsou et al. [4] examined the results of Sakiadis [2], analytically and experimentally, which includes both laminar and turbulent flow conditions. Vajravelu and Mohapatra [5] investigated the effects of injection, analytically, on the boundary layer flow past a moving sheet. Takhar et al. [6] analyzed the boundary layer flow due to a moving plate in the presence of variable fluid properties. Furthermore, Andersson and Aarseth [7] revisited the Sakiadis flow problem with the view to assess some crucial misconceptions of the flow behavior and provided numerical results for water at atmospheric pressure. Ahmad et al. [8] extended the Blasius and Sakiadis problems to nanofluids. Ding and Xin [9] used fixed-point method (FPM) and obtained semi-analytical solutions for classical Blasius flow. In recent years, several authors examined the boundary layer flow a past moving sheet or stretching sheet by considering Newtonian/non-Newtonian fluids with various constraints (Vajravelu et al. [10] and Prasad et al. ([11], [12], [13])).

The velocity slip is a phenomenon of non-adherence of the fluid to a solid boundary. More often “apparent wall slip” occur at solid boundaries for the fluids like slurries, gels, emulsions, and foam. Even for these fluids, equations of motion are governed by Navier-Stokes equations with no-slip condition. Beavers and Joseph [14], introduced the concept of a slip flow condition at the boundary. In view of this, Andersson [15] and Wang [16] considered Newtonian fluid flow past a linearly stretching sheet under partial slip (first order slip/general slip) condition. Many researchers such as Fang et al. [17], Sajid et al. [18], Matthews and Hill [19] considered partial slip in the absence of heat transfer, whereas Hayat et al. [20], Yazdi et al. [21], Sahoo [22] considered both partial slip and heat transfer character. Also, for futher applications see Zhu et al. [23], Mansur et al. [24], Zhu et al. [25], Sahoo et al. [26] and Hayat et al. [27]. The Maxwell [28] slip condition is widely accepted and implemented in current rarefied gas flow. However, the Maxwell model is only applicable for the gas flows where the rarefaction and roughness effects are not evident. Beskok and Karniadakis [29] proposed a classical second-order slip boundary condition and used it to solve the Navier-Stokes equations for confined fluids at the microscale and nanoscale. Wu [30] explored a new and advanced second-order slip velocity model. Fang et al. [31] used this model to obtain exact solutions of the governing Navier-Stokes equations.

In view of these applications, the problem studied here extends the work of Vajravelu and Mohapatra [5] to the second-order slip velocity at a moving plate. The coupled non-linear partial differential equations governing the problem are transformed into a system of coupled non-linear ordinary differential equations. The transformed equations are solved analytically via homotopy analysis method (HAM). Computed results for the flow characteristics are analyzed. The analysis of the results shows that the fluid flow is appreciably influenced by the sundry parameters. It is expected that the results obtained will not only provide useful information for industrial applications but also complement the earlier works.

Consider a flow of an incompressible viscous fluid, with constant velocity U_∞ past a parallel, porous, semi-infinite flat plate moving with a constant velocity U_w in the direction opposite to the main stream (see Fig. 1, not drawn to scale). Let the x-axis be taken along the plate, with positive x being the direction of the flow, and the y-axis is normal to it. Further, let the fluid properties be constant. Under these assumptions, the governing boundary-layer equations for the flow are

Fig. 1

∂u∂x+∂v∂y=0,$$\begin{array}{} \displaystyle \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0, \end{array}$$(1)

u∂u∂x+v∂u∂y=ν∂2u∂y2−σB02ρ2x(u−U∞),$$\begin{array}{} \displaystyle u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=\nu\frac{\partial^2 u}{\partial y^2}-\frac{\sigma B_0^2}{\rho 2x}(u-U_{\infty}), \end{array}$$(2)

where ν is the kinematic viscosity of the fluid and the other symbols have their usual meanings. The appropriate boundary conditions are

u=−(Uw+USlip),v=vw(x)=CνU∞2xaty=0,u=U∞asy→∞.$$\begin{array}{} \displaystyle u=-(U_w+U_{Slip}),\quad v=v_w(x)=C\sqrt{\frac{\nu U_{\infty}}{2x}} \quad \text{at} \quad y=0, \quad u=U_{\infty}\quad \text{as}\quad y\rightarrow \infty. \end{array}$$(3)

where u and v are the velocity components in the x-and y-directions, respectively, and U_Slip is the slip velocity at the wall. Wu’s [30] slip velocity model (valid for arbitrary Knudsen numbers K_n) is given by

USlip=233−αl3α−321−l2Knλ(x)∂u∂y−14l4+2Kn2(1−l2)λ2(x)∂2u∂y2=Aλ(x)∂u∂y+Bλ2(x)∂2u∂y2$$\begin{array}{} \displaystyle U_{Slip}=\frac{2}{3}\left(\frac{3-\alpha l^3}{\alpha}-\frac{3}{2}\frac{1-l^2}{K_n}\right)\lambda(x)\frac{\partial u}{\partial y}-\frac{1}{4}\left(l^4+\frac{2}{K_n^2}(1-l^2)\right)\lambda^2(x)\frac{\partial^2 u}{\partial y^2}=A\lambda(x)\frac{\partial u}{\partial y}+B\lambda^2(x)\frac{\partial^2 u}{\partial y^2} \end{array}$$(4)

where l = min[1/K_n, 1], α is the momentum accommodation coefficient with 0 ≤ α ≤ 1, and λ(x) = λx^1/2 is the molecular mean free path. The molecular mean free path is always positive. Thus, B < 0, and hence the second term in right hand side of Eq. (4) is positive. There are numerous applications that involve micro-scale devices including sensors, heat exchangers, and micro-power systems. Based on Knudsen number (K_n ), Beskok and Karniadakis [29] classified the gas flow in micro channels into four flow regimes: (a) continuum flow regime (K_n ≤ 0.001) (b) slip flow regime (0.001 ≤ K_n ≤ 0.1); (c) transition flow regime (0.1 ≤ K_n ≤ 10.0) and (d) free molecular flow regime (K_n > 10).

The stream function and the similarity variables can be written as

η=U∞2νxy,ψ(x,y)=2νxU∞f(η).$$\begin{array}{} \displaystyle \eta=\sqrt{\frac{U_{\infty}}{2\nu x}}y,\quad\psi(x,y)=\sqrt{2\nu x U_{\infty}}f(\eta). \end{array}$$(5)

Now, the velocity components can be expressed as

u=U∞f′(η),v=−νU∞2x[f(η)−ηf′(η)].$$\begin{array}{} \displaystyle u=U_{\infty}f'(\eta), v=-\sqrt{\frac{\nu U_{\infty}}{2x}}[f(\eta)-\eta f'(\eta)]. \end{array}$$(6)

Hence, the mass transfer velocity at the wall becomes

vw(x)=−νU∞2xf(0).$$\begin{array}{} \displaystyle v_w(x)=-\sqrt{\frac{\nu U_{\infty}}{2x}}f(0). \end{array}$$(7)

Using (5) and (6) in equations (1)-(3), we get

f‴+ff″−Mn(f′−1)=0$$\begin{array}{} \displaystyle f'''+ff''-\text{Mn}(f'-1)=0 \end{array}$$(8)

with the boundary conditions

f(0)=−C,f′(0)=−λ1−γf″(0)−δf‴(0),f′→1 as η→∞,$$\begin{array}{} \displaystyle f(0)=-C,\quad f'(0)=-\lambda_1-\gamma f''(0)-\delta f'''(0),\quad f' \to 1 \text{ as } \eta \to \infty, \end{array}$$(9)

where C is the wall mass transfer parameter, λ1=UwU∞$\begin{array}{} \displaystyle \lambda_1=\frac{U_w}{U_{\infty}} \end{array}$ is the velocity ratio parameter, γ is the first-order velocity slip parameter with 0<γ=AU∞2ν$\begin{array}{} \displaystyle 0<\gamma=A\sqrt{\frac{U_{\infty}}{2\nu}} \end{array}$ and δ is the second-order velocity slip parameter with 0>δ=BU∞2ν$\begin{array}{} \displaystyle 0>\delta=B\frac{U_{\infty}}{2\nu} \end{array}$ and Mn=σB02ρU∞$\begin{array}{} \displaystyle \text{Mn}=\frac{\sigma B_0^2}{\rho U_{\infty}} \end{array}$ is the magnetic parameter.

The physical quantity of interest is the skin friction coefficient C_f which takes the form

Cf=τwxρUw2,$$\begin{array}{} \displaystyle C_f=\frac{\tau_{wx}}{\rho U_w^2}, \end{array}$$(10)

where τ_w is the wall skin friction, given by

τw=μ∂u∂yy=0.$$\begin{array}{} \displaystyle \tau_{w}=\mu\left(\frac{\partial u}{\partial y}\right)_{y=0}. \end{array}$$(11)

From Eq. (5), the dimensionless skin friction coefficient becomes

Re2Cf=f″(0),$$\begin{array}{} \displaystyle \sqrt{\frac{Re}{2}}C_f=f''(0), \end{array}$$(12)

where Re=Uxν$\begin{array}{} \displaystyle Re=\frac{Ux}{\nu} \end{array}$ is the local Reynolds number.

For obtaining solution to the problem, we turn to homotopy analysis method (HAM). In 1992, Liao proposed HAM to serve as a general analytical method for solving nonlinear problems [32]. The accuracy and robustness of the HAM for solving nonlinear boundary value problems has been repeatedly confirmed in a wide range of papers, including: Lane-Emden equation [33], time-dependent Michaelis-Menton equation [34], non-local Whitham equation [35], and Zakharov system [36] to name a few. Here we outline the solution method and later discuss the convergence and accuracy of the HAM solution.

For the present problem, we choose the auxiliary linear operator 𝓛 as

L=∂3∂η3+β∂2∂η2,$$\begin{array}{} \displaystyle \mathscr{L}=\frac{\partial^3}{\partial \eta^3}+\beta\frac{\partial^2}{\partial\eta^2}, \end{array}$$(13)

with an initial approximation to f(η) as

f0(η)=η+λ1+1β−γβ2+δβ3e−βη−C+λ1+1β−γβ2+δβ3,$$\begin{array}{} \displaystyle f_0(\eta)=\eta+\frac{\lambda_1+1}{\beta-\gamma\beta^2+\delta\beta^3}e^{-\beta\eta}-\left(C+\frac{\lambda_1+1}{\beta-\gamma\beta^2+\delta\beta^3}\right), \end{array}$$(14)

where β > 0 is a convergence control parameter to be chosen later. From Eq. (8), we define the nonlinear operator 𝓝 as

N[f^(η,q,ℏ,β)]=∂3f^(η,q,ℏ,β)∂η3+f^(η,q,ℏ,β)∂2f^(η,q,ℏ,β)∂η2−Mn∂f^(η,q,ℏ,β)∂η−1.$$\begin{array}{} \displaystyle \mathscr{N}[\widehat{f}(\eta,q,\hbar,\beta)]=\frac{\partial^3\widehat{f}(\eta,q,\hbar,\beta)}{\partial\eta^3}+\widehat{f}(\eta,q,\hbar,\beta)\frac{\partial^2\widehat{f}(\eta,q,\hbar,\beta)}{\partial\eta^2}-\text{Mn}\left(\frac{\partial\widehat{f}(\eta,q,\hbar,\beta)}{\partial\eta}-1\right). \end{array}$$(15)

The so-called zeroth-order deformation equation is

(1−q)L[f^(η,q,ℏ,β)−f0(η)]−qℏN[f^(η,q,ℏ,β)]=0,$$\begin{array}{} \displaystyle (1-q)\mathscr{L}[\widehat{f}(\eta,q,\hbar,\beta)-f_0(\eta)]-q\hbar\mathscr{N}[\widehat{f}(\eta,q,\hbar,\beta)]=0, \end{array}$$(16)

with boundary conditions

f^(0,q,ℏ,β)=−C,f^′(0,q,ℏ,β)=−λ1−γf^″(0,q,ℏ,β)−δf^‴(0,q,ℏ,β),f^′(∞,q,ℏ,β)=1,$$\begin{array}{} \displaystyle \widehat{f}(0,q,\hbar,\beta)=-C,\quad \widehat{f}'(0,q,\hbar,\beta)=-\lambda_1-\gamma\widehat{f}''(0,q,\hbar,\beta)-\delta\widehat{f}'''(0,q,\hbar,\beta),\quad \widehat{f}'(\infty,q,\hbar,\beta)=1, \end{array}$$(17)

where q ∈ [0,1] is the embedding parameter, ℏ ≠ 0 is a convergence control parameter, and the prime denotes differentiation with respect to η. At q = 0 and q = 1, we get respectively

f^(η,0,ℏ,β)=f0(η),andf^(η,1,ℏ,β)=f(η).$$\begin{array}{} \displaystyle \widehat{f}(\eta,0,\hbar,\beta)=f_0(\eta),\quad \text{and}\quad \widehat{f}(\eta,1,\hbar,\beta)=f(\eta). \end{array}$$(18)

Therefore, as q increases from 0 to 1, f̂(η, q, ℏ, β) varies continuously from the initial approximation f₀(η) to the solution of interest f(η). By defining

fm(η,ℏ,β)=1m!dmf(η,q,ℏ,β)dηmq=0$$\begin{array}{} \displaystyle f_m(\eta,\hbar,\beta)=\left.\frac{1}{m!}\frac{d^mf(\eta,q,\hbar,\beta)}{d\eta^m}\right|_{q=0} \end{array}$$(19)

we expand f̂(η, q, ℏ, β) about q to obtain

f^(η,q,ℏ,β)=f0(η)+∑m=1∞fm(η,ℏ,β)qm.$$\begin{array}{} \displaystyle \widehat{f}(\eta,q,\hbar,\beta)=f_0(\eta)+\sum_{m=1}^{\infty}f_m(\eta,\hbar,\beta)q^m. \end{array}$$(20)

With prudent choice of the auxiliary linear operator, initial approximation, and convergence control parameters, the series in (20) can be made convergent at q = 1. The HAM series solution now becomes

f(η)=f0(η)+∑m=1∞fm(η,ℏ,β)$$\begin{array}{} \displaystyle f(\eta)=f_0(\eta)+\sum_{m=1}^{\infty}f_m(\eta,\hbar,\beta) \end{array}$$(21)

In order to obtain f_m(η, ℏ, β), we differentiate Eq. (16)m times with respect to q, divide by m!, and finally evaluate at q = 0 to arrive at the so-called m^th-order deformation equation

L[fm(η,ℏ,β)−χmfm−1(η,ℏ,β)]=ℏfm−1‴+∑n=0m−1fm−1−nfn″−Mn(fm−1′)+Mn(1−χm),$$\begin{array}{} \displaystyle \mathscr{L}[f_m(\eta,\hbar,\beta)-\chi_mf_{m-1}(\eta,\hbar,\beta)]=\hbar\left[f_{m-1}'''+\sum_{n=0}^{m-1}f_{m-1-n}f''_n-\text{Mn}(f'_{m-1})+\text{Mn}(1-\chi_m)\right], \end{array}$$(22)

with corresponding boundary conditions

fm(0,ℏ,β)=0,fm′(0,ℏ,β)=0,andfm′(∞,ℏ,β)=0,$$\begin{array}{} \displaystyle f_m(0,\hbar,\beta)=0,\quad f_m'(0,\hbar,\beta)=0,\quad \text{and}\quad f'_m(\infty,\hbar,\beta)=0, \end{array}$$(23)

where χ_m if defined as

χm=0,m≤1,1,m>1.$$\begin{array}{} \displaystyle \chi_m = \begin{cases} 0, \quad m \leq 1, \\ 1, \quad m \gt 1. \end{cases} \end{array}$$(24)

It should be noted that f_m(η, ℏ, β) for m ≥ 1 is governed by the linear Eq. (22) with boundary conditions (23). Thus, the HAM has efficaciously converted the nonlinear Eq. (8) with relevant boundary conditions into an infinite series of linear sub-problems, which can be solved.

As pointed out by Liao [37], the convergence rate and the convergence region of the HAM series solution is strongly tied to the convergence control parameter. Therefore, the optimal solution is obtained vis-á-vis optimal choice of the convergence control parameter. In 2007, Yabushita et al. [38] suggested squared residual error as a means to determine the optimal convergence-control parameter in the framework of the HAM. The exact squared residual error for the m^th-order solution is defined as

Em(ℏ,β)=∫0∞(N[fm(η,ℏ,β)])2dη.$$\begin{array}{} \displaystyle \text{E}_m(\hbar,\beta)=\int_0^{\infty}(\mathscr{N}[f_m(\eta,\hbar,\beta)])^2d\eta. \end{array}$$(25)

In practice, it is often too computationally demanding to evaluate the integral in (25), even for lower order approximations. Alternatively, we make use of an average squared residual error

Em(ℏ,β)=1N+1∑i=0N(N[∑j=0mfj(ηi,ℏ,β)])2,$$\begin{array}{} \displaystyle E_m(\hbar,\beta)=\frac{1}{N+1}\sum_{i=0}^{N}(\mathscr{N}[\sum_{j=0}^{m}f_j(\eta_i,\hbar,\beta)])^2, \end{array}$$(26)

where η_i = iΔη, and N is a positive integer.

Notice that the series in (21) gives a family of solutions in two parameters, ℏ and β.We may further reduce computation time by selecting an appropriate value of β and minimizing the error function (26) over one parameter, ℏ. It should be reiterated that the HAM provides considerable freedom in selecting the values of ℏ and β, as to ensure convergence of the series solution. On selecting an appropriate value of β, the reader is referred to Liao [39]. In that work, Liao presents a HAM solution to the Blasius flow problem and finds admissible values of β to be β > β_c, where β ≈ 2.5. For computational purposes, we set β = 3 and minimize the error function (26) to obtain the optimal value of ℏ.

We consider the following three sets of values for the parameters:

Mn = 15, C = 0.001, λ₁ = 0.3, γ = 1.0, and δ = –2.0,

and calculate the 10^th-order HAM solution with β = 3. The optimal values for the convergence control parameter ℏ are found to be, I: ℏ = –0.17916, II: ℏ = –0.24102, and III: ℏ = –0.30055. In Fig. 2, we presented the average squared residual error against order of approximation for the above three sets of parameters. For the 2^nd-order solution, we observed errors of E₂ = 2.8 × 10^–3, E₂ = 5.0 × 10^–5, and E₂ = 3.9 × 10^–6 for the parametric schemes I, II, and III, respectively. The error decreases as the order of approximation increases, reaching errors of E₁₀ = 1.9 × 10^–6, E₁₀ = 1.1 × 10^–7, and E₁₀ = 1.5 × 10^–9 for schemes I, II, and III, respectively.

Fig. 2

Solutions to the governing Eq. (8) with associated boundary conditions (9) are obtained using the HAM. The mathematical computations were executed using Mathematica 9.0. To confirm the accuracy of the present results, we compare values of the skin friction coefficient f″(0) with numerical (obtained by Runge-Kutta method with shooting technique) and approximate results reported by Vajravelu and Mohapatra [11], for the special case where the magnetic parameter Mn, first-order velocity slip parameter γ, and second-order velocity slip parameter δ are neglected. The results, reported in Table 1, are in excellent agreement. Hence, the outlined solution method in the paper can be applied to the present problem with confidence.

Table 1

Comparison of skin friction coefficient for different values of C and λ₁ in the absence of magnetic parameter Mn, first-order slip parameter γ, and second-order slip parameter δ.

Here we illustrate the effects of the physical parameters on the dimensionless velocity components f(η), f′(η), and dimensionless shear stress f″(η). That is, the effects of magnetic parameter Mn, injection parameter C, velocity ratio parameter λ₁, first-order velocity slip parameter γ > 0, and second-order velocity slip parameter δ < 0, are analyzed graphically. Fig. 3 shows the f versus η for different values of Mn and C. It can be seen that an increase in Mn or C results in a decrease in f. Additionally, f(η) is observed to increase over the entire domain.

Fig. 3

The dimensionless velocity and shear stress profiles for different combinations of parameters are displayed in Figs. 4-7. From Figs. 4 and 5, it can be seen that increasing the magnetic parameter Mn has the effect of decreasing the velocity profile, and decreasing the magnitude of the shear stress. Due to an enhanced magnetic field, a Lorentz force is produced opposite to the flow, causing the velocity profile to decrease. The effects of mass injection parameter C on the velocity and shear stress profiles are displayed in Figs. 6 and 7. From these figures we see that higher velocity and higher magnitude of shear stress occur for larger values of C. Physically, the introduction of mass in to the flow allows the flow to penetrate deeper into the fluid, causing an increase in velocity.

Fig. 4

Fig. 5

Fig. 6

Fig. 7

The influence of velocity ratio parameter λ₁ on the velocity and shear stress profiles are shown Figs. 4 and 7. We see that an increase in λ₁ results in an increase in the velocity field and increase in magnitude of the shear stress. From the velocity ratio parameter λ1=UwU∞$\begin{array}{} \displaystyle \lambda_1=\frac{U_w}{U_{\infty}} \end{array}$ , it is interesting to note that the velocity within the boundary layer increases as the free stream velocity decreases (i.e. with increasing values of λ₁ for a constant U_w ), thereby causing an increase in the velocity gradient at the surface and increase in skin friction.

The effects of first-order velocity slip parameter γ on the dimensionless velocity and shear stress profiles are displayed in Fig. 4. From this figure, the velocity field and magnitude of the shear stress decrease for increasing values of γ. Fig. 6 depicts the effects of the second-order velocity slip parameter δ on the dimensionless velocity and shear stress profiles. From the figure, we see that the velocity field and magnitude of shear stress decreases with an increase in |δ|.

To further illustrate the effects of the physical parameters on the surface drag force, the dimensionless skin friction coefficient –f″(0) as a function of different pertinent parameters is displayed in Figs. 8-10. As one can see in Fig. 8, the magnitude of the skin friction coefficient increases as λ₁ increases, but decreases with increasing Mn. Fig. 9 shows the skin friction coefficient as a function of γ for two values of C. From Fig. 10, we see that the magnitude of the skin friction coefficient decrease as |δ| increases. The functional dependence of the skin friction coefficient on C is illustrated in Figs. 9 and 10. It is noticed that the magnitude of skin friction coefficients is lower when there is no fluid injection in to the boundary layer (C = 0). Additional values for the dimensionless skin friction coefficient –f″(0) are presented in Table 2. This table further illustrates the effects of the physical parameters on the shear stress at the surface.

Fig. 8

Fig. 9

Fig. 10

Table 2

Values of dimensionless skin friction coefficient – f ″(0) for different values of Mn, C, λ₁, γ, and δ.

The viscous flow with second-order velocity slip over a porous moving flat plate subjected to a perpendicular magnetic field is investigated. The governing nonlinear partial differential equation is transformed in to a nonlinear ordinary differential equation by an appropriate similarity transformation. We solved the resulting equation analytically by the homotopy analysis method. The influence of various physical parameters on the velocity field and shear stress are analyzed and discussed. Some of the conclusions are as follows:

Velocity and shear stress values are lower in the presence of the two slip parameters;

Thus, the inclusion of mass injection and two slip parameters have significant influence on the fluid flow and shear stress at the surface of the moving plate, as well as in the fluid medium.