Flow over a non-uniform sheet with non-uniform stretching (shrinking) and porous velocities

Abstract

Viscous flow over a porous and stretching (shrinking) surface of an arbitrary shape is investigated in this article. New dimensions of the modeled problem are explored through the existing mathematical analogies in such a way that it generalizes the classical simulations. The latest principles provide a framework for unification, and the consolidated approach modifies the classical formulations. A realistic model is presented with new features in order to explain variety of previous observations on the said problems. As a result, new and upgraded version of the problem is appeared for all such models. A set of new, unusual, and generalized transformations is formed for the velocity components and similarity variables. The modified transformations are equipped with generalized stretching (shrinking), porous velocities, and surface geometry. The boundary layer governing equations are reduced into a set of ordinary differential equations (ODEs) by using the unification procedure and technique. The set of ODEs has two unknown functions f and g. The modeled equations have five different parameters, which help us to reduce the problem into all previous formulations. The problem is solved analytically and numerically. The current simulation and its solutions are also compared with existing models for specific value of the parameters, and excellent agreement is found between the solutions.

Introduction

The present progress regarding fluid flow over a porous and stretching (shrinking) sheet is incremental and substantial. The physical models of fluid flow have profound applications in extrusion processes and metal industries.^1–3 Notable achievements and investigations have been made, and the researchers are recorded new observations regarding such models. Sakiadis ⁴ studied boundary layer flow of an incompressible fluid over a continuous stretching sheet moving with a constant speed. The next significant and analytical study was presented for the variable stretching velocity by Crane. ⁵ A new similarity solution is explored by Banks ⁶ for a boundary layer flow due to stretching surface, and he provided new analytical and numerical solutions to the problem. Moreover, he also presented some special non-similar solutions of the flow problem. Later on, the Eigen solutions of boundary layer flow over stretching sheets have been examined and investigated by Banks and Zaturska. ⁷ Flow and heat transfer over a continuous linearly stretched sheet with variable surface temperature have been discussed by Grubka and Bobba ⁸ and solved the problem by using special Kummer’s function. The reference papers are mainly concerned about similarity solutions of viscous fluid flows over a porous stretching (shrinking) sheet.^9–12 McLeod and Rajagopal ¹³ examined the unique solution of Navier–Stokes equation for a viscous flow over a stretching sheet. Chang et al. ¹⁴ and Lawrence et al. ¹⁵ investigated the multiple solutions of a non-Newtonian fluid flows over a stretching sheet. Recently, lot of works have been carried out by many researchers to model and solve the problems of fluid flow related to stretching (shrinking) over a sheet of an arbitrary shape in Newtonian and non-Newtonian fluids.^16–24 This section is strengthened by citing latest research work from the literature into it. In these research articles, nano and micropolar fluid models are investigated for both Newtonian and non-Newtonian fluids over a stretching sheet. The convective flow of nano fluids is analyzed in some of the studies,^25–31 and interesting results have been pursued for flow and heat transport in such special situations. A model problem of viscous flow and heat transfer is studied over a stretching sheet for passive control of nano particles on surface of a thin layer of fluid. ²⁵ Later on, this model is further extended for non-Newtonian fluids. ²⁶ Other researchers^27–31 investigated the flow, heat, and mass transfer in fluids maintained over a stretching sheet. The important application of boundary layer flow is further enhanced in some works,^32–36 and they discussed steady flow with magnetohydrodynamic (MHD) effects, unsteady flow and heat transfer with MHD effects, unsteady mixed convection flow with magnetic effect, and unsteady natural convection flow with MHD effects for micropolar fluid. The study of MHD thermosolutal Marangoni boundary layer flow without convection is discussed by Magyari, ³⁷ and the study of fluid flow with convection over a flat surface in the presence of heat generation or absorption effects is discussed by Chamkha. ³⁸ Moreover, MHD flows of chemically reactive species of micropolar fluid in the presence of heat generation or absorption are investigated in the works of Khedr et al. ³⁹ and Damseh et al. ⁴⁰ On the contrary, numerous studies have been carried out to improve the thermal conductivity of fluids by including nano/micro or large-sized particle materials into it. A comprehensive survey of all such cases has been completed by Chamkha.^41,42 Chamkha et al. ⁴¹ found the effects of stretching surfaces on three-dimensional micropolar fluids. The mathematical simulations in Chamkha et al. ⁴¹ offer the insights of wall-driven flows into the mechanisms studied so far. Chamkha et al. ⁴² investigated the effects of uniform lateral mass flux on non-Darcy natural convection boundary layer flow of a non-Newtonian fluid over a permeable vertical cone embedded in a porous medium filled with nano fluid. Note that this investigation is completed by considering the Brownian motion and the inputs of thermophoresis. Stretching (shrinking) problems are attempted by many researchers such that they transformed the governing non-linear boundary value problem into self-similar form, and later on they found either exact or the numerical solutions of the simplified problem. There are few stretching (shrinking) problems for which exact analytical solution are known.^43,44

In this article, the results of those classical models are recovered which describe fluid flow over a variable (uniform), stretched (shrunk), and porous (variable and uniform) sheet of variable thickness. The problem is formulated in terms of the continuity equation and boundary layer momentum equation with generalized boundary conditions. The boundary conditions are taken on stretched (shrunk) and porous (injection/suction) surface of variable thickness. A set of generalized and unusual similarity transformations is formed for the velocity components and similarity variable. By means of the transformation, the governing partial differential equations (PDEs) together with boundary conditions are converted into a set of coupled boundary value problems. The last system has two ordinary differential equations (ODEs) in two unknown functions, which are the representatives of velocity components. The new modeled equations are matched with some of the studies^4,10,16–24 for fixed value of parameters involved in the problem, and details of the comparison are shown in tables and graphs. New exact, numerical, and power series solutions of the modeled problem are attempted. The approximate solutions are compared with the exact and numerical solutions, and excellent agreement is found among them. Moreover, the solutions are also compared with the published results for fixed value of the parameters, and these observations are presented in different tables and figures. Velocity and shear stress profiles are graphed for different values of all these parameters, and interesting results have been produced and discussed. New profiles are also added for the sake of comparison with classical results. The comprehensive graphs of new and classical predications ensure the difference between the two.

Problem formulation

Here, we assumed two-dimensional, steady, and laminar flow of an incompressible viscous fluid, and the flow is maintained over a stretched (shrunk) and porous sheet. The sheet of variable thickness is porous, and the porous velocity is variable, whereas the stretching/shrinking velocity is also non-uniform. Moreover, the ambient fluid is at rest. The following governing equations are taken for the modeled problem, which include continuity and boundary layer momentum equation

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$$

(1)

$$u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=\nu \frac{{\partial }^{2}u}{\partial y^2}$$

(2)

It is customary to define the no-slip condition at fluid–solid interface. Therefore, the velocity components are considered on stretched (shrunk) and porous surface, $y=f(x)$ , and the following boundary conditions are taken on the solid–fluid interface

$$u(x,y)=U\left(x\right),v(x,y)=V\left(x\right)\mathrm{when}y=f\left(x\right)$$

(3)

$$u(x,y)=0\mathrm{as}y\to \infty$$

(4)

where $U(x)=U_{0}r(x{)}^m$ and $V(x)=V_0/(r(x))$ are representing variable stretching (shrinking if $U(x)<0$ ) and injection (suction if $V(x)<0$ ) velocities of the fluid at the surface of the sheet, respectively. Note that $f(x)=\delta r(x)$ is an arbitrary shape of the surface and $q(x)=q_0/(r(x))$ is the characteristic velocity. The characteristic length is defined by $r(x)=\gamma (a_0+\alpha x{)}^{1/\beta }$ , where $\beta =2+m$ and $\gamma ={\beta }^{1/\beta }$ . A set of following transformations is introduced for the velocity components and similarity variable

$$\begin{array}{c}u=U\left(x\right)g\left(\eta \right),v=q\left(x\right)h\left(\eta \right)+V\left(x\right),\hfill \\ \mathrm{where}\eta =\frac{y-f\left(x\right)}{r\left(x\right)}\hfill \end{array}$$

(5)

The new variables in equation (5) are formed on the basis of analogies used in the study of Kármán ⁴⁵ and substituted into the governing PDEs (1) and (2), and we get the following ODEs

$${\delta }_{2}h\prime +{\delta }_3\mathrm{mg}-{\delta }_3\delta g\prime -{\delta }_3\eta g\prime =0$$

(6)

$$g″+{\delta }_2(\mathrm{gh}\prime -\mathrm{hg}\prime )-{\delta }_{1}g\prime =0$$

(7)

The boundary conditions in equations (3) and (4) (see Figure 1) are reduced as

$$g\left(0\right)=1,h\left(0\right)=0,\mathrm{and}g\left(\eta \right)=0\mathrm{when}\eta \to \infty$$

(8)

where the prime denotes differentiation with respect to $\eta$ and ${\delta }_1=V_0/\nu ,{\delta }_2=q_0/\nu ,\mathrm{and}{\delta }_3=U_0\alpha /\nu$ .

Figure 1.

Geometry of the physical model under consideration and its coordinates.

Solutions of the problem

In this section, we present the exact solution for specific value of parameters involved in the problem. A series solution coupled with Padé approximation is also found for general value of parameters. The series solution is compared with the exact and numerical solutions of the problem. The series solution is exactly matched with the other two solutions, and accuracy of the solutions is established. The exact solutions of the highly non-linear Navier–Stokes equations are few; therefore, in some cases, the closed form solutions of the governing equations are found and these recommended solutions are cited by many researchers. Moreover, the numerical solutions or computer’s generated solutions are usually compared with the closed form solutions (exact solutions) in order to determine their accuracy and establish their correctness. So these exact solutions are meaningful and applicable, and work as a baseline solution by which one can validate the numerical results. This is the most reasonable argument to compare the computer-generated solutions with the closed form solutions. The exact solutions cannot be formed for all situations, and only special cases are known (came into being) for which such solutions can be found. On the other hand, the classical exact solutions of the popular modeled problems Crane ⁵ , Banks ⁶ and Grubka & Bobba ⁸ are highly cited by many researchers for the comparison of the numerical/approximate solution of their problems. The main objective of the exact solution is to facilitate the computer-generated approximate solutions in the form of comparison. Note that a set of classical solutions is also retrieved from these solutions. In this section, we have presented four different closed form solutions of boundary value problem in equations (6)–(8). These solutions are in the form of rational, hyperbolic, and exponential functions.

Exact solution

• Case I (when ${\delta }_3<0$ and $m>-1$ )

The closed form solution, given below, is in the form of rational functions and clearly asymptotic in nature. This solution is not valid for all choices of the value of parameters. Therefore, some constraints have been imposed on the governing parameters. The denominator of these rational functions cannot be zero, so we need to avoid all those values of the parameters for which the denominator vanishes. The boundary value problem in equations (6)–(8) has a closed form solution for ${\delta }_3<0$ and $m>-1$

$$g\left(\eta \right)=\frac{s_1}{{(s_2+s_3\eta )}^2}$$

(9)

$$h\left(\eta \right)=\frac{{\delta }_3{s}_1\{(1+m)s_2+s_3(\delta +(2+m)\eta )\}}{{\delta }_2{s}_3{(s_2+s_3\eta )}^2}-\frac{{\delta }_1}{{\delta }_2}$$

(10)

where $s_1=-6s_3^2/({\delta }_3(2+m))$ , $s_3=\sqrt{(-{\delta }_3(2+m)s_2^2)/6}$ and $s_2=-\beta \delta -\sqrt{(-6{\delta }_3(1+m))/2+m}$ .

• Case II (when ${\delta }_3>0$ and ${\delta }_2\ne 0$ )

Another closed form solution of the boundary value problem in equations (6)–(8) is found by hit and trial method, and it has the following form

$$g\left(\eta \right)=\sec h{\left(\frac{\sqrt{{\delta }_3}\eta }{2}\right)}^2$$

(11)

$$h\left(\eta \right)=\frac{-\delta {\delta }_3+{\delta }_3(\delta +\eta )\sec h{\left(\frac{\sqrt{{\delta }_3}\eta }{2}\right)}^2-\sqrt{{\delta }_3}\tanh \left(\frac{\sqrt{{\delta }_3}\eta }{2}\right)}{{\delta }_2}$$

(12)

• Case III

For ${\delta }_1=-f_w$ , $m=2m_1/(1-m_1)$ , ${\delta }_3=1-m_1/1+m_1$ , $\delta =0$ , and $m_1=1$ , an exponential-type closed form solution of equations (6)–(8) is

$$h\left(\eta \right)=\frac{e^{-\frac{1}{2}{\overline{f}}_w\eta }}{{\delta }_2{\overline{f}}_w{(1+m_1)}^2(2+2m_1-2e^{\frac{1}{2}{\overline{f}}_w\eta })+{\overline{f}}_w\eta -{\overline{f}}_w{m}_1\eta }$$

(13)

$$g\left(\eta \right)=e^{\frac{-1}{2}{\overline{f}}_w\eta }$$

(14)

where ${\overline{f}}_w=(f_w+\sqrt{4+f_w^2})$ and this solution is exactly matched with the solution of modeled problem published in Magyari and Keller. ¹⁰

• Case IV

A fourth solution of equations (6)–(8) is found for special values of the parameters, that is, ${\delta }_1=-f_w$ , $m=2m_1/(1-m_1)$ , ${\delta }_3=1-m_1/1+m_1$ , $\delta =0$ , and $m_1=-(1/3)$ . The solution is given below

$$g\left(\eta \right)=\frac{{\tilde{f}}_w}{1+({\tilde{f}}_w-1)\cosh \left(\sqrt{{\tilde{f}}_w}\eta \right)+f_w\sqrt{{\tilde{f}}_w}\sinh \left(\sqrt{{\tilde{f}}_w}\eta \right)}$$

(15)

$$h\left(\eta \right)=\frac{-{\stackrel{⌣}{f}}_w+{\tilde{f}}_w\left(-1+m_1\right)\eta +{\stackrel{⌣}{f}}_w\cosh \left(\sqrt{{\tilde{f}}_w}\eta \right)+\sqrt{{\tilde{f}}_w}(1+m_1)\sinh \left(\sqrt{{\tilde{f}}_w}\eta \right)}{{\delta }_2\left(1+m_1\right){\left(\sqrt{{\tilde{f}}_w}\cosh \left(\frac{\sqrt{{\tilde{f}}_w}\eta }{2}\right)+f_w\sinh \left(\frac{\sqrt{{\tilde{f}}_w}\eta }{2}\right)\right)}^2}$$

(16)

where ${\tilde{f}}_w=2+f_w^2$ and ${\stackrel{⌣}{f}}_w=f_w(1+m_1)$ . The solutions in equations (15) and (16) are exactly matched with the solutions of modeled equations, reported in the study of Magyari and Keller ¹⁰ for fixed value of the parameters.

Power series solution

The laborious power series solution of equations (6)–(8) is presented here. Notice that the unknown functions $g(\eta )$ and $h(\eta )$ are expressed in power series form as

$$g\left(\eta \right)=\sum _{n=0}^{\infty }{a}_n{\eta }^n$$

(17)

$$h\left(\eta \right)=\sum _{n=0}^{\infty }{b}_n{\eta }^n$$

(18)

where the coefficients $a_n$ and $b_n$ are need to be determined by substituting equations (17) and (18) into equations (6) and (7). We meet the following recurrence relations for the coefficients $a_n$ and $b_n$

$$\begin{array}{c}{a}_{k+2}=\frac{1}{k+2}\hfill \\ ({\delta }_1{a}_{k+1}-\frac{{\delta }_2}{k+1}\sum _{i=0}^k(k+1-i)(a_i{b}_{k+1-i}-b_i{a}_{k+1-i}))\hfill \end{array}$$

(19)

$$b_{k+1}=\frac{{\delta }_3}{{\delta }_2}(\delta a_{k+1}-\frac{m-k}{k+1})$$

(20)

where $k\in {\mathbb{Z}}^{+}$ . The coefficients $a_0$ , $b_0$ are determined by utilizing the boundary conditions $g(0)=1$ and $h(0)=0$ , which are given in equation (8). Note that the two coefficients are reduced to $a_0=1$ and $b_0=0$ . Moreover, the arbitrary constant $a_1$ in equations (19) and (20) can be determined by applying the condition $g(\infty )=0$ and Padé approximation. The missing constant $a_1$ is also calculated by the Shooting method for different value of parameters. Four different values of $a_1$ are calculated and compared with the initial guess $g\prime (0)$ , and the results are presented in Table 1. Besides that, the two solutions (power series and numerical) for $g(\eta )$ are compared in graphs and tables for different sets of value of parameters, and excellent agreement is found between the two.

Table 1.

The arbitrary coefficient $(a_1)$ in equation (17) is compared with the missing condition $g\prime (0)$ in equation (8), when $m=-3$ , $\delta =2$ , ${\delta }_2=-1$ , and ${\delta }_3=-1/2$ .

${\delta }_1$	$a_1$	$g\prime (0)$	Absolute error
1	−0.555556	−0.5556	4e–5
−1	−1.14778	−1.1486	8e–4
1/2	−0.650488	−0.6505	1e–5
1/5	−0.720663	−0.7207	4e–5

Comparison of the power series solution with a numerical solution

The power series solution in equations (17) and (18) and numerical solution of equations (6)–(8) are compared in Table 2 for different values of ${\delta }_1$ . The absolute error of the two solutions is calculated in this table. It is noticed that the two solutions are matched and very close to each other. Therefore, the error between solutions is very small. In this table, we have taken only those values of the unknown function at some points of the similarity variable for which the error between solutions is of significant order. A comprehensive analysis of the comparison between solutions is shown in Figure 2. In this figure, the power series solution in equation (17) is compared with the numerical solution of equations (6)–(8), which is obtained by Shooting method. It is observed that the two solutions are exactly matched to each other for this set of value of parameters.

Table 2.

Power series solution in equation (17) is compared with the numerical solution of equations (6)–(8) for $m=-3$ , $\delta =2$ , ${\delta }_2=-1$ , and ${\delta }_3=-1/2$ .

$\eta$	Solution in equation (17)		Numerical method		Absolute error
	${\delta }_1=1/2$	${\delta }_1=1$	${\delta }_1=1/2$	${\delta }_1=1$	${\delta }_1=1/2$	${\delta }_1=1$
0	1	1	1	1	0	0
0.6531	0.66694	0.706856	0.667	0.7069	4e–4	4e–5
1.6327	0.385866	0.441696	0.3859	0.4417	4e–5	4e–6
2.9388	0.201227	0.253048	0.2013	0.253	8e–5	5e–5
3.5918	0.148721	0.195925	0.1489	0.1959	2e–4	2e–5
4.5714	0.0962875	0.1364	0.0967	0.1364	4e–4	0
5.2245	0.0726422	0.108406	0.0733	0.1083	7e–4	1e–4
6.5306	0.0415744	0.0699509	0.0428	0.0698	1e–3	2e–4
7.5102	0.0271059	0.0510427	0.0288	0.0507	2e–3	3e–4

Figure 2.

Axial velocity (normalized) is graphed against $\eta$ , and the power series solution in equation (17) is compared with the numerical solution of equations (6)–(8).

Comparison of the current model and its solutions with literature

The current simulation is compared with the published models for proper and fixed value of the parameters. The unknown functions of the current model are replaced by the respective functions used in the published problems. In addition, we recovered the previous models of the earlier studies^4,10,16–24 such that the parameters $(m,\delta ,{\delta }_1,{\delta }_2,{\delta }_3)$ of current model are adjusted accordingly. Moreover, if we reconsider the classical stream function formulation then we can present the velocity components $(h(\eta )\&g(\eta ))$ in the classical form as mentioned in the above stated papers. In such a situation, the comparison of the new approach and the published work is established in Table 3. Conveniently, the parameters in Table 3 have the meanings, that is, ${\delta }_4=(-c+\mathrm{AB}(d_1+d_2))/d_2$ and ${\delta }_1$ and ${\delta }_2$ are replaced by ${\lambda }_1$ and ${\lambda }_2$ , respectively, in the work of Ali et al., ¹⁷ and $m$ is replaced by $m_1$ in the earlier studies.^10,19,20,23 Note that 10 different research papers are taken from literature, and these articles are simultaneously described the fluid flow over a uniform (or non-uniform) stretching (or shrinking) and porous (uniform or non-uniform injection/suction) sheets (flat or non-flat). These papers are equipped with coupled effects of boundary sources. In all these models, the governing equations are transformed into a single ODE by using appropriate similarity transformations. The model presented in this article is the generalization of all such models, which are relevant to a specific mechanism. The governing equations are transformed with the help of unusual and generalized similarity transformation to a set of ODEs. In Table 4, $g\prime (0)$ is computed from equations (6)–(8) numerically by Shooting method and compared with $f″(0)$ of earlier studies^{10,16,17,20,22} for fixed value of parameters. For these choices of value of parameters, the current results are exactly matched with classical data, which are taken from different papers. It is known fact that the plate suction is a very practical means of delaying boundary layer transition to turbulence. Note that for $m=-3$ , ${\delta }_3=-1/2$ , and different ${\delta }_1$ , we are failed to calculate $g\prime (0)$ for different values of $\delta$ . It is confirmed that for $m=-3$ , ${\delta }_1=-1$ , and for any value of $\delta$ , $g\prime (0)$ cannot be evaluated for a wide range of ${\delta }_3\succ 0$ (injection), so the values recorded in Table 5 for $\delta$ and ${\delta }_3$ are beyond the critical values. The observation in Figure 6 enhances the facts about the existence of critical values. It is worth mentioning that the critical values are difficult to determine exactly since they are found by successive integration of the ODEs in equations (6)–(8). The velocity profiles in Figure 3 change shape gradually from a near-Sakiadis flow at the leading edge to the asymptotic exponential, hyperbolic, and general rational function profiles of equations (6)–(8) for downstream (Table 6).

Table 3.

The current model in equations (6)–(8) is compared with published papers.

S. No.	$g(\eta )$	$h(\eta )$	$m$	$\delta$	${\delta }_1$	${\delta }_3$	Published models
$1.$	$\frac{1}{B}f\prime (\eta )$	$\frac{-1}{{\delta }_2}({\delta }_1+{\lambda }_{2}f(\eta )-({\lambda }_1+{\lambda }_2)\eta f\prime (\eta ))$	$\frac{-{\lambda }_1}{{\lambda }_1+{\lambda }_2}$	$0$	$-{\lambda }_2{\delta }_4$	$B({\lambda }_1+{\lambda }_2)$	Ali et al. 17
$2.$	$f\prime (\eta )$	$\frac{-1}{{\delta }_2}({\delta }_1+f(\eta )+\frac{m_1-1}{m_1+1}\eta f\prime (\eta ))$	$\frac{2m_1}{1-m_1}$	$0$	$-f_w$	$\frac{1-m_1}{1+m_1}$	Magyari and Keller 10
$3.$	$f\prime (\eta )$	$\frac{-1}{2{\delta }_2}(2{\delta }_1+f(\eta )+(-1+2\beta )\eta f\prime (\eta ))$	$\frac{2\beta }{1-2\beta }$	$0$	$-\frac{\alpha }{2}$	$\frac{1-2\beta }{2}$	Liao 22
$4.$	$-f\prime (\eta )$	$\frac{-1}{{\delta }_2}({\delta }_1+f(\eta )+(m_1-1)\eta f\prime (\eta ))$	$\frac{m_1}{1-m_1}$	$0$	$-s$	$m_1-1$	Fang 19
$5.$	$f\prime (\eta )$	$\frac{-1}{{\delta }_2(1+m_1)}((1+m_1)({\delta }_1+f(\eta ))+(m_1-1)\eta f\prime (\eta ))$	$\frac{2m_1}{1-m_1}$	$0$	$-\alpha \left(\frac{1-m_1}{1+m_1}\right)$	$\frac{1-m_1}{1+m_1}$	Fang et al. 20
$6.$	$-f\prime (\eta )$	$\frac{1}{2{\delta }_2{m}_1}(-2{\delta }_1{m}_1-(1+m_1)f(\eta )+(1-m_1)\eta f\prime (\eta ))$	$\frac{2m_1}{1-m_1}$	$0$	$-\frac{1+m_1}{2m_1}{f}_w$	$\frac{-1+m_1}{2m_1}$	Prasad et al. 23
$7.$	$-f\prime (\eta )$	$\frac{-1}{{\delta }_2}({\delta }_1+f(\eta )-2\eta f\prime (\eta ))$	$\frac{-1}{2}$	$0$	$2\alpha$	$-2$	Asghar et al. 16
$8.$	$-f\prime (\eta )$	$\frac{-1}{{\delta }_2(1+n)}((1+n)({\delta }_1+f(\eta ))+(n-1)\eta f\prime (\eta ))$	$\frac{2n}{1-n}$	$0$	$S$	$\frac{1-n}{1+n}$	Ghosh et al. 21
$9.$	$f\prime (\eta )$	$\frac{-1}{{\delta }_2}(f(\eta )+\eta f\prime (\eta ))$	$-2$	$0$	$0$	$-1$	Al-Mdallal et al. 18 and Sajjid and Hayat 24
$10.$	$f\prime (\eta )$	$\frac{-1}{{\delta }_2}(f(\eta )-\eta f\prime (\eta ))$	$0$	$0$	$0$	$\frac{1}{2}$	Sakiadis 4

Table 4.

$g\prime (0)$ is computed from equations (6) to (8) numerically when $m=-1/2$ , $\delta =0$ , and ${\delta }_2=-1$ and compared with $f″(0)$ of the previous works.^{10,16,17,20,22}

${\delta }_1,{\delta }_3$	$1,2$	$-1/2,1$	$-2,2$	$-2,-2$	$1,2$
$f″(0)$ from equations (16) to (18) in the work of Magyari and Keller 10	$1$	–	–	–	–
$f″(0)$ from equation (14) in the work of Liao 22	–	$-0.5$	–	–	–
$f″(0)$ from equations (9) and (13) in the work of Fang et al. 20	–	–	$-2$	–	–
$f″(0)$ from equation (10) in the work of Asghar et al. 16	–	–	–	$2$	–
$f″(0)$ from equations (10) to (12) in the work of Ali et al. 17	–	–	–	–	$1$
$g\prime (0)$ from equations (6) to (8) (Current results)	$1$	$-0.5$	$-2$	$2$	$1$

Table 5.

Critical values for different ranges of parameters when $m=-3,{\delta }_2=-1$ .

${\delta }_3=-1/2\&{\delta }_1\downarrow$	Critical values for $\delta$	${\delta }_1=-1\&\delta \downarrow$	Critical values for ${\delta }_3$
$4$	No critical value	$-3$	${\delta }_3\succ 0$
$3$	$\delta \succ 92.6$	$0$	${\delta }_3\succ 0$
$2$	$\delta \succ 94.5$	$1$	${\delta }_3\succ 0.1$
$1$	$\delta \succ 96.6$	$2$	${\delta }_3\succ 0.2$
$0$	$\delta \succ 98.6$	$3$	${\delta }_3\succ 4.2$

Figure 3.

Axial velocity (normalized) is graphed against $\eta$ and new solutions of equations (6)–(8) and their comparison with previous plots where the parameters are given in Table 6.

Table 6.

Different values of parameters are fixed for each profile in Figure 3.

S. No.	$m$	$\delta$	${\delta }_1$	${\delta }_2$	${\delta }_3$	Resulting graph
$1.$	$-1/2$	$0$	$-4.5$	$-1$	$-10$	Ali et al. 17
$2.$	$-1/2$	$0$	$1.5$	$-1$	$2$	Magyari and Keller 10
$3.$	$-20/19$	$0$	$-1/2$	$-1$	$-19/2$	Liao 22
$4.$	$-3/8$	$0$	$-2.5$	$-1$	$-8/5$	Fang 19
$5.$	$-5/2$	$0$	$0.3333$	$-1$	$-2/3$	Fang et al. 20
$6.$	$-1/2$	$0$	$-5.4$	$-1$	$-2$	Asghar et al. 16
$7.$	$-6$	$0$	$2.5$	$-1$	$-0.2$	Ghosh et al. 21
$8.$	$-2$	$0$	$0$	$-1$	$-1$	Al-Mdallal et al. 18 and Sajjid and Hayat 24
$9.$	$0$	$0$	$0$	$-1$	$1/2$	Sakiadis 4
$10.$	$-1/2$	$1/2$	$1$	$-1$	$0.4$	Current results
$11.$	$-1/2$	$-1$	$1$	$-1$	$0.4$	Current results

Discussion

In this section, the field properties are elaborated for fixed value of parameters appeared in modeled equations (6)–(8). The effects of all these parameters are seen on the axial and normal velocity components $(h(\eta )\&g(\eta ))$ , and the results are presented in different graphs and tables. Moreover, shear stress is found from Newton’s law of viscosity. The skin friction at the walls is computed against the parameters ${\delta }_1\&\delta$ and graphed in two different figures.

Shear stress

The shear stress due to moving fluid can be obtained from the following Newton’s law of viscosity

$${\tau }_{\mathrm{xy}}=\mu (\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x})$$

(21)

In view of equation (5), equation (21) becomes

$$\begin{array}{c}\multicolumn{1}{c}{{\tau }_{\mathrm{xy}}=\mu [\frac{U\left(x\right)}{r\left(x\right)}g\prime \left(\eta \right)-\frac{\alpha }{\beta (a_0+\alpha x)}(V\left(x\right)+q\left(x\right)h\left(\eta \right)}\\ \multicolumn{1}{c}{+\frac{U\left(x\right)\alpha (\delta +\eta )}{r{\left(x\right)}^{m+1}}((\delta +\eta )g\prime \left(\eta \right)-\mathrm{mg}\left(\eta \right)))]}\end{array}$$

(22)

The shear stress at the surface of the sheet (skin friction) is modified and non-dimensionalized, which is presented by the following relation

$$\hat{\tau }=\frac{\overline{\tau }}{\theta }=g\prime \left(0\right)$$

(23)

where $\overline{\tau }={\tau }_{\mathrm{xy}}+(\alpha \mu /(\beta (a_0+\alpha x)))(V(x)-(m\alpha \delta U(x)/(r(x{)}^{m+1})))$ and $\theta =(\mu U(x)/(r(x)))(1-({\alpha }^2{\delta }^2)/(\beta (a_0+\alpha x)r{(x)}^m))$ . The modified shear stress is furnished with stretching (shrinking), injection (suction), and boundary layer inputs. Note that this simulation is independent of stream function formulation, and the skin friction coefficient is dependent upon the gradient of axial velocity.

Figure 4 shows the profiles of skin friction coefficients $(g\prime (0))$ against ${\delta }_1$ which is growing thicker as $\delta$ decreases from positive to negative. However, the skin friction is strongly dependent on ${\delta }_1$ . It is approaching to zero at the separation point and breaks down completely for non-similar flows in the region where skin friction is zero. Similar behavior of skin friction coefficient is noted for different values of ${\delta }_1$ and $\delta$ in Figures 5 and 6, respectively. In Figure 7, the effects of ${\delta }_3$ are seen on the profiles of axial velocity for different values of ${\delta }_1$ . There are three sets of profiles in this figure: one set is representing by solid lines, the other one by dotted lines, and the last one is denoted by dashed lines. In each set of profiles, the velocity and boundary layer thickness are increased with the increasing of ${\delta }_3(<0)$ . A thin boundary layer is seen for different values of ${\delta }_3$ and ${\delta }_1=-1$ . In this figure, the maximum axial velocity and a thick boundary layer are observed for ${\delta }_3=-1/2$ and ${\delta }_1=2$ , and all the profiles are concaved up. In Figure 8, effects of ${\delta }_3$ are observed on profiles of axial velocity and the profile has highest peak for smaller value of ${\delta }_3$ with maximum boundary layer thickness. The velocity is increased with the decreasing of ${\delta }_3$ and profiles are concaved down. Besides that, each profile is increased smoothly against the similarity variable $\eta$ . It starts decreasing which is approaching to zero asymptotically outside the boundary layer. Moreover, the axial velocity is increased in the vicinity of the sheet surface and then starts decreasing as we are getting away from the plate surface. Similar behavior of the axial velocity is noted in Figure 9 for different values of ${\delta }_3$ . The boundary layer thickness is decreased with the increasing of ${\delta }_3$ . In Figure 9, the peaks of the profiles are higher than the peaks of the profiles in Figure 8. From Figures 8 and 9, it is observed that the axial velocity is changed rapidly with ${\delta }_1$ and increased triple times as ${\delta }_1$ increased from 1 to 2.

Figure 4.

Skin fraction coefficient is plotted against ${\delta }_1$ .

Figure 5.

Skin fraction coefficient is plotted against $\delta$ .

Figure 6.

Skin fraction coefficient is plotted against ${\delta }_3$ .

Figure 7.

Axial velocity (normalized) is graphed against $\eta$ , and the numerical solution of equations (6)–(8) is plotted for different values of parameters.

Figure 8.

Axial velocity (normalized) is graphed against $\eta$ , and the numerical solution of equations (6)–(8) is plotted for different values of ${\delta }_3$ .

Figure 9.

Axial velocity (normalized) is graphed against $\eta$ , and the numerical solution of equations (6)–(8) is plotted for different values of ${\delta }_3$ .

In Figure 10, two sets of profiles are presented for the normal velocity component. One set of profiles is the reflection of the other about the similarity line $\eta =0$ . In the first set, positive velocity profiles are decreased with the increasing of ${\delta }_2(<0)$ . Meanwhile, the other set consists of negative velocity profiles and increased with increasing ${\delta }_2(>0)$ . The normal velocity profiles reached different values of constants outside the boundary layer. The constants are dependent upon the values of ${\delta }_2$ .

Figure 10.

Normal velocity (normalized and translated) is graphed against $\eta$ , and the numerical solution of equations (6)–(8) is plotted for different values of ${\delta }_2$ .

Conclusion

In this article, we have studied the flow of viscous fluid over a porous stretching/shrinking sheet of a variable thickness. Closed form solutions of the problem are found for specific ranges of the value of parameters. The power series solutions of the general problem are obtained with the combination of Padé approximation. This solution is in excellent agreement with the numerical solution. Furthermore, a detailed comparative study of the current model is presented, and remarkable results are retrieved that are available in literature on the topic. Following are the key findings of this research work:

• Direct simulation of the governing PDEs is taken and pursued their similarity solutions. A set of transformations are defined for the velocity components and similarity variable. Note that the stream function formulation is not considered in the modeled problem. As a result, all the previous models of stream function formulation are recovered easily.

• Multiple solutions of the final governing ODEs are found. The new solutions of the current problem are compared with classical simulations, and excellent agreement is found between them. The comparison is established in different tables and graphs.

• Different and multiple solutions (closed form and series) of the model equations are presented and exactly matched with the numerical solution.

• A new set of solutions is found, and these solutions are missing in the classical simulations. Ultimately, a set of new profiles in Figures 8–10 are presented and such solutions are not seen in the literature. It is one of the reasons that the flow behavior in these graphs is totally different.