**1. Introduction**

In a heat transfer mechanism, fluid is a main medium as a heat transfer carrier. Therefore, improving the thermal transfer efficiency of the fluid used is a vital challenge in the industry. Certain experiments have shown that the thermal conductivity of fluids containing metal and oxide particles is higher than that of traditional base liquids such as oil, water, and ethylene glycol [1–3]. For the sake of improving the heat transfer efficiency of the fluid, researchers have added metal and non-metallic nanoparticles into the traditional base liquid to form a new compound "nanofluid". Nanofluids are made up of base fluids and nanoparticles, but not a simple mixture, which are composed of nano-sized solid particle or tubes suspended in the base fluids, are solid–liquid composite materials. Nanoparticles have high surface-activity and tend to aggregate together with time. The idea was first proposed by Choi and Eastman [4]. Nanofluids are important in the fields of energy, chemical, microelectronics, and information. Recently, the flow and conduct heat of nanofluids have been studied by certain scholars. A quick overview is given here. Sheremet et al. [5] discussed natural convection of alumina-water nanofluid in an inclined wavy-walled cavity. Nanofluids flow in microchannels with heat conduction was discussed by Bowers et al. [6]. Hashim et al. [7] discussed the mixed convection and heat conduction of Williamson nanofluids under unsteady condition. Mahdy [8] presented the effects of magnetohydrodynamics (MHD) and variable wall temperature on non-Newtonian Casson nanofluid flow. Asadi et al. [9] presented the latest progress of preparation methods and thermophysical properties of oil-based nanofluids. Pourfattah et al. [10] simulated water/CuO nanofluid fluid flow and heat transfer inside a manifold microchannel. Alarifi et al. [11] investigated

the effects of solid concentration of nanoparticles, temperature, and shear rate on the rheological properties of nanofluid. For a traditional base fluid, there are two main types: Newtonian fluids and non-Newtonian fluids. In industry, non-Newtonian fluids play an important role, such as juices, starch solutions, egg whites, and apple pulp. To understand behavious of non-Newtonian fluids, certain models have been presented. Power law model is relatively simple, widely used among these models. Researchers have further investigated the flow and conduct heat of power law fluids. Javanbakht et al. [12] studyed the heat conduction on the surface of a power law fluid. Turan et al. [13] discussed mixed convection of power-law liquids in enclosures. The heat conduction of power law liquid in various section tubes was considered by Zhang et al. [14]. Ahmedet et al. [15] addressed MHD power law liquid flow in a Darcy–Brinkmann porous medium.

In this paper, the base fluid of a nanofluid is power law fluid. When nanoparticles are added into the traditional base liquid, local velocity slip may happen as an effect of high shear force between the fluid and the wall, and the slip condition is no longer negligible in the nanometer or micro scales. The velocity slip is a finite velocity boundary condition between the fluid and the solid [16]. Researchers have done certain studies on the slipping problems of nanofluids. Ramya et al. [17] studied the viscous flow and heat transfer of nanofluid through a stretched sheet with the effect of magnetic field, velocity, and thermal slip. Abbas et al. [18] discussed the stagnation flow of micropolar nanofluids through a cylinder with slip. The effect of heat and velocity slip on the flow of Carson nanofluids through a cylinder was discussed by Usman et al. [19]. Babu et al. [20] investigated the three-dimensional MHD nanofluid flow over a variable thickness slendering stretching sheet with the effect of thermophoresis, Brownian motion, and slip parameter. The above studies all discussed the first-order slip model, whereas higher-order slips should be considered when the velocity and temperature profiles of an average free path are nonlinear. It is now known that the inclusion of higher-order slip yields results closer to those by experiments [21]. Thus, various investigations on higher-order slip flows were published by Uddin et al. [22], Kamran et al. [23], Farooq et al. [24], and Yasin et al. [25]. These all sugges<sup>t</sup> that the power law constitutive equation should be considered on the basis of high order slip for a power law nanofluid.

In the aforementioned literature, there are few papers about the flow and heat transfer of magnetic nanofluids with higher-order slip parameters. Therefore, a new mathematical model is proposed. With the help of similarity transformation variables, governing equations are converted to ordinary differential equations, whose solution is solved using homotopy analysis method. The effects of nanofluid velocity, temperature, concentration, skin friction coefficient and Nusselt number on various physical parameters are simulated. In addition, the fluid flow situation is visualized by the computational fluid dynamics (CFD) software Ansys Fluent.

#### **2. Mathematical Modelling Formulation**

#### *2.1. Flow Behavior*

Consider a steady, two-dimensional, incompressible MHD fluid flow with copper through a stretching thin plate. All variables mentioned are presented in Tables 1 and 2 [26] gives some physical capabilities of the base liquid and nanoparticles. Meanwhile, a transverse magnetic field is utilized, where the strength is *Bx* and the presence of surface tension is also considered. Given the above hypotheses, the governing equations composed of continuity equation and momentum equation can be given as

$$
\frac{\partial \Pi}{\partial X} + \frac{\partial V}{\partial Y} = 0,\tag{1}
$$

$$
\hbar \frac{\partial \mathcal{U}}{\partial X} + V \frac{\partial \mathcal{U}}{\partial Y} = -\frac{1}{\rho\_{nf}} \frac{\partial P}{\partial X} + \frac{\partial S\_{XX}}{\partial X} + \frac{\partial S\_{XY}}{\partial Y} + \frac{\sigma B^2}{\rho\_{nf}} (\mathcal{U}\_t - \mathcal{U}), \tag{2}
$$

$$
\Delta U \frac{\partial V}{\partial X} + V \frac{\partial V}{\partial Y} = -\frac{1}{\rho\_{\rm nf}} \frac{\partial P}{\partial Y} + \frac{\partial S\_{YX}}{\partial X} + \frac{\partial S\_{YY}}{\partial Y}, \tag{3}
$$

$$\mathcal{S}\_{ij} = 2\mu\_{nf} \left( 2D\_{ml} D\_{ml} \right)^{\frac{n-1}{2}} D\_{ij}, \\ D\_{ij} = \frac{1}{2} \left( \frac{\partial \mathbf{U}\_i}{\partial \mathbf{X}\_j} + \frac{\partial \mathbf{U}\_j}{\partial \mathbf{X}\_i} \right). \tag{4}$$


**Table 1.** Nomenclature.

In the above, *X* and *Y* are the Cartesian coordinates along and normal to the extension sheet, respectively. **U** is the velocity field. *U* and *V* are the x and y components of **U**. *P* is the pressure, *σ* the electric conductivity, *Bx* the magnetic field along the forward direction of *Y*-axis, *Ue* the free stream speed, *Sij* the deviatoric part of the stress tensor *ςij* = −*Pδij* + *Sij*, *δij* the unit tensor, and *Dij* the rate-of-strain tensor. *ρnf*the effective density and *μnf*the effective dynamic viscosity given by [27]

$$
\rho\_{nf} = (1 - \varrho)\rho\_f + \varrho\rho\_{s\prime} \quad \mu\_{nf} = \frac{\mu\_f}{(1 - \varrho)^{2.5}}.\tag{5}
$$

The other parameters of nanofluid (*ρCp*)*n f* , *αn f* , *kn f* are given [27]

$$(\rho \mathbb{C}\_p)\_{nf} = (1 - \rho)(\rho \mathbb{C}\_p)\_f + \varphi(\rho \mathbb{C}\_p)\_{sf}, \quad a\_{nf} = \frac{k\_{nf}}{(\rho \mathbb{C}\_p)\_{nf}},\tag{6}$$

$$\frac{k\_{nf}}{k\_f} = \frac{k\_s + 2k\_f - 2\varphi(k\_f - k\_s)}{k\_s + 2k\_f + \varphi(k\_f - k\_s)}\text{ }\tag{7}$$

where subscripts *s*, *f* , and *n f* represent the solid particle, base liquid, and the thermophysical properties of nanofluid, respectively. *ϕ* is the solid volume fraction of nanoparticles, (*ρCp*)*n f* the effective heat capacity. The thermal conductivity is *kn f* and the thermal diffusivity is *αn f* .

For the sake of analyzing the boundary layer in a better way, the following nondimensional variables are introduced,

$$\text{tr} = \frac{X}{L}, y = \frac{Y}{\delta}, u = \frac{\mathcal{U}}{\mathcal{U}\_w}, v = \frac{\mathcal{U}V}{\delta \mathcal{U}\_w}, p = \frac{P}{\rho\_f \mathcal{U}\_w^2}, \tau\_{\text{ij}} = \frac{S\_{\text{ij}}}{\rho\_f \mathcal{U}\_w^2},\tag{8}$$

where *L* and *δ* represent the characteristic length in the *X* and *Y* direction, respectively. *Uw* denotes the velocity in the *X*-direction.

Thus, Equations (1)–(4) become

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

$$u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = -\frac{\rho\_f}{\rho\_{nf}}\frac{\partial p}{\partial x} + \frac{\partial}{\partial y}(\frac{\mu\_{nf}\rho\_f}{\rho\_{nf}}\left|\frac{\partial u}{\partial y}\right|^{n-1}\frac{\partial u}{\partial y}) + \frac{\rho\_f \sigma B^2}{\rho\_{nf}}(u\_\varepsilon - u),\tag{10}$$

$$\frac{\partial p}{\partial y} = 0.\tag{11}$$

From Equation (11), it can be concluded that the pressure *p* is identical with the pressure of mainstream flow.

$$-\frac{\partial p}{\partial \mathbf{x}} = \mu\_{\mathbf{e}} \frac{\partial \mu\_{\mathbf{e}}}{\partial \mathbf{x}}.\tag{12}$$

For a power law nanofluid, velocity slip effect need be considered. In many investigations, the first-order model is adopted widely. The model is suitable under the assumption that temperature and velocity profiles are linear through a average free path. However, when temperature and velocity profiles are nonlinear through a average free path, higher-order slip would become possible. Mitsuya [28] has obtained a second-order slip model from a physical phenomenon by considering the accommodation coefficient:

$$F = u f\_1 m\_1 \left[ (\frac{2}{3}\lambda) \frac{\partial u}{\partial y} + \frac{1}{2} (\frac{2}{3}\lambda)^2 \frac{\partial^2 u}{\partial y^2} + u\_{slip} \right]|\_{y=0\nu} \tag{13}$$

where *F* is the shear stress, *α* an accommodation coefficient relative to momentum, *f*1 the frequency of molecular bombardment, *m*1 the molecular mass density, and *λ* the local molecular average free path.

In this paper, as the base fluid is a power flow fluid, namely, the shear stress *F* = *μn f* --- *∂u∂y* --*<sup>n</sup>*−1 *∂u∂y* , the constitutive equation of a power flow fluid with a higher-order slip is considered. The enhanced slip model is written as

$$u(\mathbf{x},0) = \mathcal{U}\_w + \left(A\_1 \frac{\partial u}{\partial y} + A\_2 \frac{\partial^2 u}{\partial y^2} + A\_3 \left|\frac{\partial u}{\partial y}\right|^{n-1} \frac{\partial u}{\partial y}\right)|\_{y=0\prime} \tag{14}$$

$$w(\mathbf{x},0) = 0, \mu(\mathbf{x}, \infty) = \mu\_{\mathfrak{e}} = a\mathbf{x}^{\mathfrak{m}},\tag{15}$$

where *A*1, *A*2, and *A*3 denote the velocity slip coefficients; *Uw* is the the speed of the stretch plate; and *Uw* = *cxm*.

For the sake of deriving a simplified model by converting governing equations into ordinary differential equations, a stream function *ψ*(*<sup>x</sup>*, *y*) is introduced in this paper such that *u* = *∂ψ∂y* , *v* = <sup>−</sup>*∂ψ∂x* . Then Lie-group transformationsis also introduced to obtain a new set of similar variables.

$$\Gamma: \mathfrak{x}^\* = \mathfrak{x} \mathfrak{e}^{\mathfrak{x}\_1}, \quad y^\* = y \mathfrak{e}^{\mathfrak{x}\_2}, \quad \psi^\* = \psi \mathfrak{e}^{\mathfrak{x}\_3}, \quad \mathfrak{u}^\* = \mathfrak{u} \mathfrak{e}^{\mathfrak{x}\_4}, \quad \upsilon^\* = \upsilon \mathfrak{e}^{\mathfrak{x}\_5}, \quad \mathfrak{u}\_\varepsilon^\* = \mathfrak{u}\_\varepsilon \mathfrak{e}^{\mathfrak{x}\_6}.\tag{16}$$

Equation (16) can be considered as a point-transformation of coordinates (*<sup>x</sup>*, *y*, *ψ*, *u*, *v*, *ue*) into coordinates (*x*<sup>∗</sup>, *y*<sup>∗</sup>, *ψ*<sup>∗</sup>, *<sup>u</sup>*<sup>∗</sup>, *<sup>v</sup>*<sup>∗</sup>, *u*∗*e* ). Substituting Equation (16) in Equation (10), we ge<sup>t</sup>

$$\begin{split} & \quad \epsilon^{\varepsilon(a\_{1}+2a\_{2}-2a\_{3})} (\frac{\partial \Psi^{\*}}{\partial y^{\*}}, \frac{\partial^{2} \Psi^{\*}}{\partial x^{\*} \partial y^{\*}} - \frac{\partial \Psi^{\*}}{\partial x^{\*}}, \frac{\partial^{2} \Psi^{\*}}{\partial y^{\*2}}) \\ &= & \quad -\epsilon^{\varepsilon(a\_{1}-2a\_{6})} \frac{\rho\_{f}}{\rho\_{nf}} u^{\*}\_{\varepsilon} \frac{d u^{\*}\_{\varepsilon}}{dx^{\*}} + \epsilon^{\varepsilon(3a\_{2}-a\_{3})} \frac{\mu\_{nf} \rho\_{f}}{\rho\_{nf}} \frac{\partial^{3} \Psi^{\*}}{\partial y^{\*3}} \epsilon^{\varepsilon(n-1)(2a\_{2}-a\_{3})} (-\frac{\partial^{2} \Psi^{\*}}{\partial y^{\*2}})^{n-1} \\ &+ \quad \frac{\rho\_{f} \sigma B^{2}}{\rho\_{nf}} (\epsilon^{-a\_{6} \varepsilon} u^{\*}\_{\varepsilon} - \epsilon^{(a\_{2}-a\_{3}) \varepsilon} \frac{\partial \Psi^{\*}}{\partial y^{\*}}). \end{split} \tag{17}$$

The boundary condition Equations (14) and (15) become

$$\begin{split} \frac{\partial \psi^\*}{\partial y^\*} (\mathbf{x}^\*, 0) &= \mathfrak{e}^{\mathfrak{e}(\mathfrak{a}y - a\_2 - \mathfrak{a}\mathfrak{a}\_1)} c \mathfrak{x}^{\*m} + A\_1 \mathfrak{e}^{\mathfrak{a}\mathfrak{a}\_2} \frac{\partial^2 \psi^\*}{\partial y^{\*2}} + A\_2 \mathfrak{e}^{\mathfrak{a}\mathfrak{a}\_2} \frac{\partial^3 \psi^\*}{\partial y^{\*3}} \\ &+ A\_3 \mathfrak{e}^{(\mathfrak{n}\mathfrak{e}(2a\_2 - \mathfrak{a}\_3) + \mathfrak{a}\mathfrak{y} - \mathfrak{a}\_2)} (-\frac{\partial^2 \psi^\*}{\partial y^{\*2}})^{n-1} (-\frac{\partial^2 \psi^\*}{\partial y^{\*2}}), \quad \text{at} \quad y^\* = 0; \end{split} \tag{18}$$

$$\frac{\partial \Psi^\*}{\partial \mathbf{x}^\*} (\mathbf{x}^\*, 0) = 0, \quad \text{at} \quad y^\* = 0; \tag{19}$$

$$\frac{\partial \psi^\*}{\partial y^\*} (x^\*, \infty) = \epsilon^{x(a\_3 - a\_2 - ma\_1)} \alpha x^{\*m}, \quad \text{at} \quad y^\* \to \infty. \tag{20}$$

The system will remain unaltered under the group of transformations Γ, so the parameters have the following relations, namely,

$$
\pi\_2 + \pi\_4 - \pi\_3 \quad = \quad \pi\_1 + \pi\_5 - \pi\_3 = \pi\_3 - \pi\_1 - m \\
\pi\_2 = \pi\_3 - \pi\_2 - m \\
\pi\_1 = 0,\tag{21}
$$

$$2a\_2 - 2a\_3 + a\_1 = \quad n(2a\_2 - a\_3) + a\_2 = (n+1)a\_2 - na\_4 = a\_2 - a\_4 - a\_5. \tag{22}$$

Thus, Equation (16) becomes

$$\begin{array}{ll} \Gamma: & \mathfrak{x}^\* = \mathfrak{x} \mathfrak{e}^{\mathfrak{su}\_1}, & \mathfrak{y}^\* = \mathfrak{y} \mathfrak{e}^{\frac{m n - 2n + 1}{n + 1} a\_1 \mathfrak{e}}, & \mathfrak{y}^\* = \mathfrak{y} \mathfrak{e}^{\frac{2m \mathfrak{w} - m + 1}{n + 1} a\_1 \mathfrak{e}}, \\ & \mathfrak{u}^\* = \mathfrak{u} \mathfrak{e}^{\mathfrak{m} \mathfrak{a}\_1 \mathfrak{e}}, & \mathfrak{v}^\* = \mathfrak{v} \mathfrak{e}^{\frac{2m \mathfrak{w} - m - \mathfrak{a}}{n + 1} a\_1 \mathfrak{e}}. \end{array} \tag{23}$$

Based on the above Lie-group transformations, the stream function and similar parameter can be prescribed as follows,

$$
\eta = \left(\frac{c^{2-n}}{\mu\_f}\right)^{\frac{1}{n+1}} \ge^{\frac{2n-mv-1}{n+1}} y\_\prime \psi = \left(\frac{\mu\_f}{c^{1-2n}}\right)^{\frac{1}{n+1}} \ge^{\frac{2mu+1-m}{n+1}} f(\eta). \tag{24}
$$

After further similarity transformations, a nonlinear ordinary differential equation is obtained.

$$m f^{\prime\prime} |f^{\prime}|^{n-1} + m \varphi\_1 (d^2 - f^{\prime} f^{\prime}) + \varphi\_1 \varphi\_2 \frac{2mn - m + 1}{n + 1} f f^{\prime\prime} + \varphi\_1 M(d - f^{\prime}) = 0. \tag{25}$$

The boundary condition Equations (14) and (15) now develop into

$$f(0) = 0,\\ f'(\infty) = d,\tag{26}$$

$$f'(0) = 1 + \lambda\_1 f''(0) + \lambda\_2 f'''(0) + \lambda\_3 \left| f''(0) \right|^{n-1} f''(0),\tag{27}$$

where *d* = *ac* , *M* is the Hartmann number with *M* = *σ<sup>B</sup>*20 *c* , *λ*1, *λ*2, and *λ*3 are velocity slip parameters; these parameters and *ϕ*1, *ϕ*2 [27] can now be written as

$$
\lambda\_1 = A\_1 (\frac{\mathfrak{c}^{2-n}}{\mu\_f})^{\frac{1}{n+1}} \mathbf{x}^{\frac{2m-nn-1}{n+1}}, \\
\lambda\_2 = A\_2 (\frac{\mathfrak{c}^{2-n}}{\mu\_f})^{\frac{2}{n+1}} \mathbf{x}^{\frac{2(2m-nn-1)}{n+1}}, \tag{28}
$$

$$\lambda\_3 = A\_3 \left( c x^m (\frac{c^{2-n}}{\mu\_f})^{\frac{1}{n+1}} x^{\frac{2m-nu-1}{n+1}} \right)^n,\tag{29}$$

$$
\varphi\_1 = (1 - \varphi)^{2.5}, \quad \varphi\_2 = 1 - \varphi + \varphi \frac{\rho\_s}{\rho\_f}, \tag{30}
$$

where *A*1, *A*2, and *A*3 are arbitrary positive constants.

#### *2.2. Heat and Mass Transfer Behavior*

The heat and mass equations can now be formulated as follows,

$$\begin{split} \mathcal{U}\frac{\partial T}{\partial X} + \mathcal{V}\frac{\partial T}{\partial Y} &= \frac{\partial}{\partial Y} \left( k(T) \frac{\partial T}{\partial Y} \right) \\ &+ \frac{\tau}{\mu\_f} \mathcal{C}\_f \frac{3}{n+1} (\mathcal{C}^3 x^{3m-1}) \frac{n-1}{n+1} \left( D\_B \frac{\partial \mathcal{C}}{\partial Y} \frac{\partial T}{\partial Y} + \frac{D\_T}{T\_\infty} \left( \frac{\partial T}{\partial Y} \right)^2 \right), \end{split} \tag{31}$$

$$
\lambda \mathcal{U} \frac{\partial \mathcal{C}}{\partial X} + V \frac{\partial \mathcal{C}}{\partial Y} = \mu\_f \, ^{\frac{2}{n+1}} (\mathcal{C}^3 \mathbf{x}^{3m-1}) ^{\frac{n-1}{n+1}} \left( D\_B \frac{\partial^2 \mathcal{C}}{\partial Y^2} + \frac{D\_T}{T\_\infty} \frac{\partial^2 T}{\partial Y^2} \right), \tag{32}
$$

$$k(T) = \frac{k\_{nf}}{(\rho \mathbb{C}\_p)\_{nf}} (T\_w - T\_\infty)^{1-n} \mathcal{U}\_w^{n-1} \left| \frac{\partial T}{\partial Y} \right|^{n-1} \,. \tag{33}$$

The boundary conditions are as follows,

$$T(X,0) = T\_w + k\_{\rm nf}(T\_w - T\_\infty)^{1-n} \left| \frac{\partial T}{\partial Y} \right|^{n-1} \frac{\partial T}{\partial Y}|\_{y=0\prime} \tag{34}$$

$$\mathbb{C}(X,0) = \mathbb{C}\_{\mathfrak{w}\_{\prime}} T(X,\infty) = T\_{\mathfrak{w}\_{\prime}} \mathbb{C}(X,\infty) = \mathbb{C}\_{\mathfrak{w}\_{\prime}} \tag{35}$$

where *T* shows temperature in the boundary layer, *T*∞ denotes the temperature away from the sheet and is a constant, and *Tw* indicates the unified temperature of the fluid. *C* is the concentration of the fluid, *C*∞ is the fluid concentration in the free stream, and *Cw* the unified concentration of the fluid.

For the sake of gaining the similarity solutions of equations, the following similarity variables are introduced,

$$\theta(\eta) = \frac{T - T\_{\infty}}{T\_{\text{av}} - T\_{\infty}},\\\phi(\eta) = \frac{\mathbb{C} - \mathbb{C}\_{\infty}}{\mathbb{C}\_{\text{av}} - \mathbb{C}\_{\infty}}. \tag{36}$$

Then, Equations (31)–(33) become

$$n\rho\_4 \theta'^\prime |\theta'|^{n-1} + \frac{2mn - m + 1}{n + 1} Pr\rho\_3 f \theta' + PrNb\rho \gamma \phi' \theta' + PrNt\rho \gamma \theta'^2 = 0,\tag{37}$$

$$
\phi'' + \frac{2mn+1-m}{n+1} \text{Scf} \phi' + \frac{Nt}{Nb} \theta'' = 0. \tag{38}
$$

The boundary conditions Equations (34) and (35) are converted to

$$\theta(0) = 1 + \beta \theta'(0) |\theta'(0)|^{n-1}, \quad \theta(\infty) = 0,\tag{39}$$

$$
\phi(0) = 1, \quad \phi(\infty) = 0,\tag{40}
$$

where *Pr* denotes Prandtl number, *Nt* represents thermophoresis parameter, *Sc* is Schmidt number, and *Nb* is Brownian motion parameter. The above parameters, *ϕ*3, *ϕ*4, and *β*, are defined as

$$Pr = \frac{\mu\_f}{a\_f}, Nb = \frac{\tau D\_B(\mathbb{C}\_{\text{av}} - \mathbb{C}\_{\text{os}})}{\mu\_f}, Nt = \frac{\tau D\_T(T\_{\text{av}} - T\_{\text{os}})}{\mu\_f T\_{\text{os}}},\tag{41}$$

$$\varphi\_{\mathcal{B}} = 1 - \varphi + \varphi \frac{\left(\rho \mathbb{C}\_{\mathcal{P}}\right)\_{\mathfrak{s}}}{\left(\rho \mathbb{C}\_{\mathcal{P}}\right)\_{f}},\\\varphi\_{\mathfrak{k}} = \frac{k\_{\mathfrak{s}} + 2k\_{f} - 2\varphi(k\_{f} - k\_{\mathfrak{s}})}{k\_{\mathfrak{s}} + 2k\_{f} + \varphi(k\_{f} - k\_{\mathfrak{s}})},\tag{42}$$

$$\beta = \frac{k\_{nf}}{\mu\_{nf}^{\frac{n}{n+1}}} (\mathbb{C}^{2n-1} X^{2mn-n-m})^{\frac{1}{n+1}},\\ \mathcal{Sc} = \frac{\mu\_f}{D\_B}. \tag{43}$$

Momentous physical parameters are expressible as follows,

$$\mathcal{C}\_f = \frac{\mu\_{nf} \left| \frac{\partial u}{\partial y} \right|^{n-1} \frac{\partial u}{\partial y} \Big|\_{y=0}}{\frac{1}{2} \rho\_f u\_w^2} = \frac{|f''(0)|^{n-1} f''(0)}{(1-q)^{2.5}} \text{Re}\_x^{-\frac{1}{n+1}},\tag{44}$$

$$\mathcal{C}\_f \mathcal{R} e\_x^{-\frac{1}{n+1}} = \frac{|f''(0)|^{n-1} f''(0)}{(1-q)^{\frac{n}{2.5}}},\tag{45}$$

$$Nu\_x = -\frac{xk\_{nf}\frac{\partial T}{\partial y}|\_{y=0}}{k\_f(T\_w - T\_\infty)} = -\frac{k\_s + 2k\_f - 2\varrho(k\_f - k\_s)}{k\_s + 2k\_f + \varrho(k\_f - k\_s)} \text{Re}\_x \overset{\perp}{\tau \to 1} \theta'(0),\tag{46}$$

$$Nu\_{\lambda} Re\_{\lambda}^{-\frac{1}{\pi+1}} = -\frac{k\_{\sf s} + 2k\_f - 2\varrho(k\_f - k\_s)}{k\_{\sf s} + 2k\_f + \varrho(k\_f - k\_s)}\theta'(0),\tag{47}$$

$$Sh\_{\mathcal{X}} = -\frac{x D\_B \frac{\partial C}{\partial y}|\_{y=0}}{D\_B(C\_{\mathcal{w}} - C\_{\infty})} = -R c\_{\mathcal{X}} \frac{1}{r^{\frac{1}{n+1}}} \phi'(0),\tag{48}$$

$$S\hbar\_{\mathbf{x}}\mathcal{R}\mathbf{e}\_{\mathbf{x}}^{-\frac{\mathbf{l}}{\mathbf{n}+\mathbf{l}}}=-\phi'(0).\tag{49}$$

#### **3. Solution Procedures**

In this section, the homotopy analysis method (HAM) [29] is used to solve this problem. The initial guess solutions of velocity, temperature, and concentration, based on boundary conditions, are, respectively,

$$f\_0 = B\_1 + B\_2 e^{-\eta} + B\_3 \eta,\\ \theta\_0 = B e^{-\eta},\\ \phi\_0 = e^{-\eta}. \tag{50}$$

Three linear operators are selected as

$$\mathcal{L}\_f = f^{\prime\prime} + f^{\prime\prime},\\\mathcal{L}\_\theta = \theta^{\prime\prime} + \theta^{\prime},\\\mathcal{L}\_\phi = \phi^{\prime\prime} - \phi. \tag{51}$$

These operators satisfy some properties:

$$L\_f(\mathbb{C}\_1 + \mathbb{C}\_2 \varepsilon^{-\eta} + \mathbb{C}\_3 \eta) = 0,\\ L\_\theta(\mathbb{C}\_4 \varepsilon^{-\eta} + \mathbb{C}\_5) = 0,\\ L\_\phi(\mathbb{C}\_6 \varepsilon^{-\eta} + \mathbb{C}\_7 \varepsilon^{\eta}) = 0 \tag{52}$$

where *Ci*(*i* = 1, 2, ··· , 7) are arbitrary constants.

The 0-th order deformation equations and its boundary conditions are derived and the expressions are written as

$$(1-p)L[F(\eta,p)-f\_0(\eta)] = p\hbar\_f H\_f(\eta)N\_f[F(\eta,p)],\tag{53}$$

$$\mathbb{E}\left[(1-p)L[\Theta(\eta,p)-\theta\_0(\eta)]\right] = p h\_{\theta} H\_{\theta}(\eta) N\_{\theta}[F(\eta,p), \Theta(\eta,p), \Phi(\eta,p)],\tag{54}$$

$$(1-p)L[\Phi(\eta,p)-\phi\_0(\eta)] = ph\_\theta H\_\theta(\eta)N\rho[F(\eta,p),\Theta(\eta,p),\Phi(\eta,p)];\tag{55}$$

$$F(0, p) = 0,\\ F'(\infty, p) = d,\\ \Theta(\infty, p) = 0,\\ \Phi(0, p) = 1,\\ \Phi(\infty, p) = 0,\tag{56}$$

$$F'(0, p) = 1 + \lambda\_1 F''(0, p) + \lambda\_2 F'''(0, p) + \lambda\_3 |F''(0, p)|^{n-1} F''(0, p),\tag{57}$$

$$
\Theta(0, p) = 1 + \beta \Theta\_0^{\prime}(0, p) |\Theta\_0^{\prime}(0, p)|^{n-1}. \tag{58}
$$

In the above equations, *p* ∈ [0, 1] is the embedding parameter; *hf* , *hθ*, and *hφ* are auxiliary non-zero parameters; and *Hf*(*η*), *Hθ* (*η*), and *<sup>H</sup>ϕ*(*η*) are nonzero auxiliary functions [30]. Obviously, for *p* = 0 and *p* = 1, we have

$$\begin{array}{ll} F(\eta,0) = f\_0(\eta), & F(\eta,1) = f(\eta), \\ \Theta(\eta,0) = \theta\_0(\eta), & \Theta(\eta,1) = \theta(\eta), \\ \Phi(\eta,0) = \phi\_0(\eta), & \Phi(\eta,1) = \phi(\eta). \end{array} \tag{59}$$

As *p* increases from 0 to 1, *<sup>F</sup>*(*η*, *p*) is from the initial guess *f*0(*η*) to the exact solution *f*(*η*), <sup>Θ</sup>(*η*, *p*) is from the initial guess *<sup>θ</sup>*0(*η*) to the exact solution *<sup>θ</sup>*(*η*), and <sup>Φ</sup>(*η*, *p*) is from the initial guess *φ*0(*η*) to the exact solution *φ*(*η*) [30]. With Taylor's theorem, they can write

$$F(\eta, p) \quad = \left. F(\eta, 0) + \sum\_{k=1}^{+\infty} f\_k(\eta) p^k, f\_k(\eta) \right. = \frac{1}{k!} \frac{\partial^k F(\eta, p)}{\partial p^k}|\_{p=0\prime} \tag{60}$$

$$\Theta(\eta, p) \quad = \quad \Theta(\eta, 0) + \sum\_{k=1}^{+\infty} \theta\_k(\eta) p^k, \quad \theta\_k(\eta) = \frac{1}{k!} \frac{\partial^k \Theta(\eta, p)}{\partial p^k}|\_{p=0}. \tag{61}$$

$$\Phi(\eta, p) \quad = \quad \Phi(\eta, 0) + \sum\_{k=1}^{+\infty} \Phi\_k(\eta) p^k, \quad \Phi\_k(\eta) = \frac{1}{k!} \frac{\partial^k \Phi(\eta, p)}{\partial p^k}|\_{p=0}. \tag{62}$$

Assuming that the auxiliary parameters *hf* , *hθ*, and *hφ* are appropriate chosen, we can obtain convergen<sup>t</sup> solutions in the following form.

$$f(\eta) = f\_0(\eta) + \sum\_{k=1}^{\infty} f\_k(\eta), \theta(\eta) = \theta\_0(\eta) + \sum\_{k=1}^{\infty} \theta\_k(\eta), \phi(\eta) = \phi\_0(\eta) + \sum\_{k=1}^{\infty} \phi\_k(\eta). \tag{63}$$

For the sake of getting the higher order deformation equation, differentiating the 0-th order deformation Equations (53)–(55) *k* times with regard to *p*, set *p* = 0 and divide by *k*!, to attain

$$L\_f(f\_k(\eta) - \chi\_k f\_{k-1}(\eta)) = h\_f H\_f(\eta) R\_{f,k}(\eta),\tag{64}$$

$$L\_{\theta}(f\_{\theta}(\eta) - \chi\_{\theta}f\_{\theta-1}(\eta)) = h\_{\theta}H\_{\theta}(\eta)R\_{\theta,k}(\eta),\tag{65}$$

$$L\_{\Phi}(f\_{\Phi}(\eta) - \chi\_{\Phi}f\_{\Phi^{-1}}(\eta)) = h\_{\Phi}H\_{\Phi}(\eta)R\_{\Phi,k}(\eta),\tag{66}$$

where *Rf* ,*<sup>k</sup>*(*η*), *<sup>R</sup>θ*,*<sup>k</sup>*(*η*), and *<sup>R</sup>φ*,*<sup>k</sup>*(*η*) are, respectively,

$$\begin{split} & \quad R\_{f,k}(\eta) \\ &= \chi\_k \sum\_{l=0}^{k-2} f\_l^{\prime\prime\prime} \sum\_{j=2}^{k-1} \sum\_{\substack{i\_1 + 2 \cdots + k \\ i\_1 + 2 \cdots + (i-1)k\_1 - k^{-1} = l \\ 1 + 2j + \cdots + (1-j)k\_1 - k^{-1} = l \end{split}} |f\_0^{\prime\prime}|^{n-j} \prod\_{q=1}^{k-1} |f\_q^{\prime\prime\prime}|^{i\_q} \\ & \quad + n f\_{k-1}^{\prime} |f\_0^{\prime\prime}|^{n-1} - m \varrho\_1 \varrho\_2 \sum\_{i=0}^{k-1} f\_i^{\prime\prime} f\_{k-1-i}^{\prime} \\ & \quad + \varrho\_1 \varrho\_2 \frac{2nm - m + 1}{n+1} \sum\_{i=0}^{k-1} f\_i f\_{k-1-i}^{\prime\prime} - \varrho\_1 M f\_{k-1}^{\prime} \,^{i} \,^{i} \,^{i} \end{split} \tag{67}$$

$$\begin{split} & \quad R\_{\theta,k}(\eta) \\ &= \chi\_k \sum\_{l=0}^{k-2} \theta\_l^{\prime\prime} \sum\_{j=2}^{k-1} \sum\_{\begin{subarray}{c} i\_1 + i\_2 + \cdots + i\_{k-1} = j-1 \\ i\_1 + 2i\_2 + \cdots + (k-1)i\_{k-1} = k-l \end{subarray}} \frac{n(n-1)\cdots(n-j+1)\varrho\_4}{i\_1!i\_2!\cdots i\_{k-1}!} |\theta\_0^{\prime}|^{n-j} \prod\_{q=1}^{k-1} |\theta\_q^{\prime}|^{i\_q} \\ & \left. + n\varrho\_4 \theta\_{k-1}{}^{\prime\prime} |\theta\_0^{\prime}|^{n-1} + \frac{2nn-n+1}{n+1} \Pr\_{q>3} \frac{k-1}{j\_1} \theta\_1 \theta\_{k-1-i}{}^{\prime} \right. \\ & \left. + \Pr N b \, q \varphi\_3 \sum\_{i=0}^{k-1} \phi\_i \theta\_{k-1-i}{}^{\prime} + \Pr N t \, q \varphi\_3 \sum\_{i=0}^{k-1} \theta\_i^{\prime} \theta\_{k-1-i}{}^{\prime} \right. \end{split} \tag{68}$$

$$R\_{\phi,k}(\eta) = \phi\_{k-1}{}^{\eta} + \frac{2mn+1-m}{n+1} \text{Sc} \sum\_{i=0}^{k-1} f\_i \phi\_{k-1-i}{}^{\prime} + \frac{Nt}{Nb} \theta\_{k-1}{}^{\prime\prime} \tag{69}$$

$$\chi\_k = \begin{cases} 0 & k \le 1, \\ 1 & k > 1. \end{cases} \tag{70}$$

Boundary conditions Equations (56)–(58) become

$$f\_k(0) = 0,\\ f\_k^{\prime}(\infty) = 0,\\ \theta\_k(\infty) = 0,\\ \phi\_k(0) = 0,\\ \phi\_k(\infty) = 0,\tag{71}$$

$$\begin{aligned} &\quad f\_{k}^{\prime\prime}(0) \\ = \sum\_{l=0}^{k-1} f\_{l}^{\prime\prime}(0) \sum\_{j=2}^{k+1-l} \sum\_{i\_{1}2,\cdots,i\_{k}=0 \\ \quad \atop i\_{1}+2j\_{2}+\cdots+i\_{k}=j-1 \\ \dots+\lambda\_{5}f\_{k}^{\prime\prime}(0)|f\_{0}^{\prime\prime}(0)|^{n-1}+\lambda\_{1}f\_{k}^{\prime\prime}(0)+\lambda\_{2}f\_{k}^{\prime\prime\prime}(0) \end{aligned} \tag{72}$$

$$\begin{array}{ll} \theta\_{k}(0) \\ = \sum\_{l=0}^{k-1} \theta\_{l}^{\prime}(0) \sum\_{j=2}^{k+1-l} \sum\_{\substack{i\_{1}i\_{2}, \dots, i\_{k}=0 \\ i\_{1}+2, \dots, i\_{k}=j-1 \\ i\_{1}+2, \dots, i\_{k}=k-l \end{array}} \frac{\theta(n-1)(n-2)\cdots(n-j+1)}{i\_{1}!i\_{2}!\cdots i\_{k}!} |\theta\_{0}{}^{\prime}(0)|^{n-j} \prod\_{q=1}^{k} |\theta\_{q}{}^{\prime}(0)|^{i\_{q}} \\ \tag{73}$$

#### **4. Results and Discussion**

In homotopy analysis, the *h*-curves are ploted to select the effective region of parameter *h*. For the sake of obtaining the convergen<sup>t</sup> parameters *hf* , *hθ*, and *<sup>h</sup>φ*, Figures 1–3 plot the *h*-curves of various orders for *f* (0), *θ*(0) and *φ*(0). Ranges of *h*-curves are [−0.4, 0], [−0.5, −0.3], [−0.5, 0.3], that is, the horizontal segmen<sup>t</sup> of the curves, which is called the effective region, so *hf* = *hθ* = *hφ* = *h* = −0.35 is selected in the paper.

For the sake of proving the accuracy and effectiveness of homotopy analysis after determining values of *hf* , *hθ*, and *<sup>h</sup>φ*, Figure 4 plots the error curves of various power law index by the "BVPh2.0" procedure software package. As can be seen from Figure 4, the errors have reached 10−<sup>4</sup> in the second order, meeting the standards of engineering calculation. The larger the order, the smaller the error becomes. Further, surface friction coefficients are compared with the literature [31] for various first-order slip parameter *λ*1 in Table 3.

**Table 2.** The physical capabilities of base fluid and nanoparticles [26].


**Table 3.** Comparisons of *Cf R n*+1 *e* for various *λ*1 as *n* = 1, *m* = 1, *d* = 1.5, *λ*2 = *λ*3 = 0, *ϕ* = 0.

1


**Figure 2.** *hθ* -curves.

After attesting the accuracy and effectiveness of homotopy analysis, the impacts of various physical parameters are analyzed, such as nondimensional velocity *f* (*η*), temperature *<sup>θ</sup>*(*η*), etc.

Meanwhile, the flow of power law nanofluid is numerically simulated by the widely used software Ansys Fluent to further explore the flow properties.

**Figure 4.** Total error of approximation for various powers *n*.

#### *4.1. Behavior of Velocity Profiles*

Figures 5 and 6 demonstrate effects of power law exponential of the plate *m* and Hartmann number *M* on nondimensional velocity *f* (*η*). The velocity distribution for various *m* is showed in Figure 5. By increasing the power exponent of the plate *m*, the tensile speed of the plate increases. Greater deformation is effected in the fluid, leading to the increase of *f* (*η*). As pointed out in [32], the effects of *M* on *f* (*η*) are visible in Figure 5. Recall that Hartmann number *M* expresses the ratio of electromagnetic force to viscous force. Due to the fact that greater Hartmann number corresponds to larger Lorenz force, the velocity *f* (*η*) increases.

When the fluid is pseudoplastic and expansive, impacts of *d* on *f* (*η*) are illustrated in Figure 7. In Figure 7, the velocity of the fluid has upward tendency for various *d*. Whereas, the velocity of expansive fluid increases slower than that of pseudoplastic fluid due to the increase of the fluid viscosity.

**Figure 7.** Impacts of *d* on *f* (*η*) for *n* < 1 and *n* > 1.

Figure 8 clearly presents the impacts of various power law index *n* on *f* (*η*). As seen in Figure 8, thebuoyancybecomeslargerasthepowerlawindex *n*increases,whichcausestheincreaseofvelocity.

Influences of different velocity slip parameters *λ*1, *λ*2, and *λ*3 on *f* (*η*) are illustrated in Figures 9–11, respectively. Velocity slip mainly affects slip loss and, in a cascade, fluid velocity. With the increases of the second-order slip parameter *λ*2, velocity *f* (*η*) also increases; however, the results are contradictory when the first-order linear slip parameter *λ*1 and nonlinear slip parameter *λ*3 increase.

**Figure 8.** Impacts of *n* on *f* (*η*).

**Figure 10.** Effects of *λ*2 on *f* (*η*).

**Figure 11.** Effects of *λ*3 on *f* (*η*)).

#### *4.2. Behavior of Temperature Profiles*

Figures 12 and 13 indicate various temperature behavior for different *Nb* and *Nt*. Figure 12 displays the effects of *Nb* on temperature. Fluid particles generate more heat through random motions when *Nb* increases, which causes the rise in temperature. Figure 13 clearly shows temperature distribution for various thermophoresis parameter *Nt*. Thermophoresis indicates that particles move from a high temperature part to a low temperature one in a fluid with temperature gradient. Thus, the temperature increases with the enhancement of the parameter *Nt*.

**Figure 13.** Impacts of *Nt* on *<sup>θ</sup>*(*η*).

Figures 14 and 15 show temperature distribution for diverse temperature jump parameter *β* and power law index *n*. Figure 14 plots the temperature curves for diverse *β*. Increasing temperature jump parameter *β* leads to a rise in the thickness of temperature boundary layer. Thus, the temperature has an upward tendency. Figure 15 demonstrates the temperature distribution for various *n*. The temperature diminish when the power law index rises. In other words, temperature boundary layer becomes thinner with the enhancement of *n*.

**Figure 15.** Impacts of *n* on *<sup>θ</sup>*(*η*).

#### *4.3. Behavior of Concentration Profiles*

Figures 16 and 17 show the concentration distribution for diverse values of the Brownian motion parameter *Nb* and the thermophoresis parameter *Nt*. From Figure 16, the collision of fluid particles rises with the stronger Brown motion, which leads to the reduction of fluid concentration. Figure 17 indicates the concentration field for various thermophoresis parameter *Nt*. The magnitude of concentration variation is greater under the influence of thermophoresis parameter.

**Figure 17.** Impacts of *Nt* on *φ*(*η*).

#### *4.4. Analysis of Skin Friction and Nusselt Number*

In the study of fluids, vital physical parameters, such as skin friction coefficient and local Nusselt number, are discussed. In this paper, the impacts of various parameters on these two parameters are demonstrated in Table 4. Skin friction coefficients have ascending behavior with the increase of *ϕ*, *λ*1 and *λ*3. On the contrary, the downward trend is seen with the raise of *λ*2. For local Nusselt number, when *ϕ* and *λ*2 rise, the local Nusselt numbers have an upward trend, whereas the local Nusselt numbers diminish with the rise of *λ*1, *λ*3 and *β*.

#### *4.5. Simulated Behavior*

In this subsection,the laminar model is used to solved governing equations. Ansys Fluent uses the Gauss-–Siedel point-by-point iterative method combined with the algebraic multigrid (AMG) method to solve the algebraic equations. The effects of various parameters on the flow of power-law nanofluid over a stretched thin sheet are simulated. The computational results obtained by using CFD solver are compared with the available results of Chen [33] for some limiting conditions. The present results are proved to be in good agreemen<sup>t</sup> as shown in Table 5. The effects of various parameters such as power law exponential of the plate *m*, nanoparticle volume fraction *ϕ*, and power law index *n* on Nusselt number *Nu* and skin friction coefficient are shown in Figures 18–22.


**Table 4.** Effects of *ϕ*, *λ*1, *λ*2, *λ*3, and *β* on *Cf Re* 1 *n*+1 *x* and *NuxRe* − 1 *n*+1 *x* for *n* = 1/2, *m* = 0, *M* = 1, *d* = 1, *Pr* = 1, *Nb* = 1, *Nt* = 1, and *Sc* = 1.

**Table 5.** Comparisons of *Cf R* 1 *n*+1 *e* for various *n* with *m* = 0.5.


The velocity contours for nonlinear slip are simulated in Figure 18. From these diagrams, the flow produces velocity boundary layer near the entrance. Besides, the velocity boundary layer of the pseudoplastic fluid is thicker than that for a Newton and expansive fluid.

Figures 19 and 20 present the effect of nanoparticle volume fraction on Nusselt number *Nu* and skin friction coefficient *Cf* at fixed values of inlet velocity, power law index. From Figure 19, the local Nusselt number increases at any x-location When nanoparticles are added to the base fluid. This is because a lower local temperature difference between the sheet walls and fluid can be achieved. Therefore, the high thermal conductivity of Cu nanoparticles enhances the thermal performance of the fluid. As the viscosity of the liquid can be increased by adding Cu nanoparticles into the base fluid, the *Cf* along the thin sheet increases when using higher concentrations of nanoparticles, as shown in Figure 20.

**Figure 20.** Effect of *ϕ* on *Cf* .

Figure 21 shows the effect of power law index *n* on skin friction coefficient *Cf* . The skin friction coefficient decreases with the increase of x-location for a given power law index. However, for a constant x-location, the skin friction coefficient have an upward tendency as the power law index increases.

Figure 22 demonstrates the skin friction coefficient distribution for various *ϕ*. The skin friction coefficient increases as the fluid behavior changes from shear-thinning to shear-thickening for a certain *ϕ*. As the *ϕ* increases the skin friction coefficient increases for a constant power law index.

**Figure 21.** Effect of *n* on skin friction coefficient.

**Figure 22.** Variation of the skin friction coefficient at different *ϕ*.

#### **5. Conclusions**

The flow and heat transfer of magnetic nanofluid through a stretched thin sheet with higher-order slip parameters are discussed in the paper. The model contains the influences of Brown motion and thermophoresis impacts. Simplified ODEs are obtained by a series of similarity transformations. The similar solutions are solved through homotopy analysis theory and the stability of the solutions is analyzed. Moreover, the current results are shown to be in good agreemen<sup>t</sup> with the literature results, the error of Nusselt number and skin friction coefficient is less than 0.1%. The key conclusions follow.


**Author Contributions:** J.Z. conducted the original research, modified the model and contributed analysis tools. X.H. analyzed the data, simulated the modified model and prepared original draft. Y.X. made numerical simulation with ANSYS software. J.Z. and Y.X. revised the manuscript.

**Funding:** This paper was supported by the National Natural Science Foundation of China (No. 11772046; No. 81870345).

**Acknowledgments:** The authors would like to express their gratitude to the reviewers for their suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest.
