**1. Introduction**

The heat transfer phenomenon plays a vital role when the temperature varies between different bodies or parts of the same body. Heat can be transferred by three methods: convection of fluids, conduction in solids, and radiation. To examine body heat transfer, the principle of heat transfer can be applied to the human body. Here, we can quote one example: the metal pan is used to transfer heat from the stove to food. Some applications of heat transfer are cooking food over metal pots, boiling milk in metal pots, and thermal treatment of pain by a hot water bag. Straughan [1] considered the Cattaneo–Christov (CC) model for heat flux and thermal convection over a Newtonian fluid. Khan et al. [2] analyzed it numerically by engaging the bvp4c MATLAB-based function on a Sisko fluid flow accompanied by generalized Fick's and Fourier's laws over a nonlinear stretched surface. Hayat et al. [3] examined analytically the magnetohydrodynamics (MHD) flow of Jeffrey fluid past a variable thick surface via the impacts of the chemical reaction and the CC model in a

stratified medium. Waqas et al. [4] investigated the CC heat flux model for energy equation formulation rather than Fourier's law of heat conduction. It was observed that a variable thermal conductivity remains inversely proportional to a temperature profile. The Soret–Dufour effects on walls with the second-grade fluid flow between inclined parallel plain walls were inspected by Khan et al. [5]. Heat transfer characteristics of the incompressible flow of the second-grade fluid flow produced by a stretching sheet were analyzed by Ghadikolaei et al. [6]. Khan et al. [7] investigated the behavior of homogenous–heterogeneous reactions against heat transfer flow due to a stretching sheet. They noticed that homogenous–heterogeneous reactions reduce fluid concentration.

In modern engineering processes, especially in metallurgical engineering and metalworking practices, the role of MHD is fundamental for electrically conducting fluids. The magnetic field function is crucial in cooling the hot plasma inside a nuclear reactor vessel. Similarly, the magnetic field is employed for the mixing of metals inside an electrical furnace [8]. Chamkha et. al. [9] analyzed the magnetic field effect on the mixed convection unsteady flow in an ambient fluid past a cone rotating with an unsteady angular velocity. Pullepu et al. [10] analyzed the free convection flow with variable surface temperature over a nonisothermal vertical cone. Akbar et al. [11] examined the two dimensional (2D) electrically conducting flow of the hyperbolic tangent fluid past a stretching surface. They observed that an increment in the Hartmann number decelerates the fluid velocity in the domain of the stretching sheet. Seini et al. [12] d the magnetic field impact over a stretching surface accompanied by the appearance of slip velocity near the stagnation point flow. They perceived that the impact of the magnetic field is more significant on the velocity profile. Ravindran et al. [13] considered the impact of a transverse magnetic field and heat generation and absorption on time-dependent mixed convection flow over a porous cone with a chemical reaction. Boland et al. [14] simulated MHD flow of viscous fluid over a circular cylinder covered with a permeable layer. They adopted the Darcy–Brinkman–Forchheimer model to study the flow inside a porous medium. Ellahi et al. [15] investigated the influence of a Hall current on MHD Jeffrey fluid flow over a nonuniform duct. Mishra et al. [16] explored heat transfer and mass in the appearance of a magnetic field of viscoelastic fluid flow. They determined that the behavior of the magnetic field against the velocity profile is opposite to temperature distribution and concentration profiles. Hussain et al. [17] analyzed numerically the influence of the applied magnetic field on a non-Newtonian fluid flow past a stretching surface.

Nanofluids are vital in many engineering applications, such as in biomedical engineering and many chemical processes. Nanofluid is composed of nanometer-sized particles with a diameter of less than 100 nm and some conventional fluid. The basic aim of using nanofluids is to upgrade the heat transfer and thermal conductivity to attain better cooling. Khan et al. [18] d nanofluid flow over a stretching sheet. Makinde [19] extended the work of Khan et al. [18] to convective boundary condition in nanofluid flow. Nadeem et al. [20] analyzed the second-grade nanofluid (nonorthogonal stagnation point) flow in the direction of the stretching surface. The impact of the variable magnetic field on the nanofluid flow between two disks was explored by Hayami et al. [21]. They found the analytical solution via the homotopy perturbation method and observed that the temperature of the boundary layer thickness decreases with the increase of the Brownian motion parameter and thermophoretic factor. Nanofluid flow in a permeable medium over a convectively heated permeable shrinking sheet was examined by Hayat et al. [22]. Sheikholeslami et al. [23,24] considered the behavior of the magnetic field on the free and forced convection flow of nanofluids respectively by making use of the two-phase model. Hassan et al. [25] elaborated on convective transport of heat transfer in a nanofluid through a porous medium. They concluded that convective heat transfer is improved by nanoparticle concentration, and the magnetic field impacts second-order slip flow. Nayak et al. [26] focused on the numerical solution of the three-dimensional (3D) nanofluid flow with nonlinear thermal radiation with convective conditions and slip. Hosseini et al. [27] observed the nanofluid MHD flow in a microchannel heat sink via the KKL (Koo–Kleinsteuer–Li) model. They noticed that the interaction between nanoparticles and the solid phase enhances the Nusselt number. In recent years, several scientists have used nanofluid heat transfer in their studies [28–34].

Fluids are basically divided into two groups: Newtonian [35] and non-Newtonian fluids [36]. Fluids that abide by Newton's law of viscosity are termed Newtonian fluids. However, a contradiction to Newton's law of viscosity is observed in non-Newtonian fluids. Applications of non-Newtonian fluids may be found in many industrial and engineering areas, glass fiber, hot rolling, casting, and paper production. Amongst many non-Newtonian fluids, Williamson fluid possesses shear thinning property (i.e., viscosity tends to decrease when shear stress increases). Abundant articles may be found in the literature that highlight the importance of Williamson fluid in numerous scenarios. Ramzan et al. [37] examined Williamson nanofluid flow over a Riga plate. They found that with an increase in Williamson fluid parameter velocity, distribution decreases. Ramzan et al. [38] also analyzed the numerical solution of the 2D MHD stagnation point of Williamson fluid flow under the effect of homogeneous–heterogeneous reactions over a linearly stretched surface and found opposite behavior of temperature and velocity distribution against the Williamson fluid parameter. Nadeem et al. [39,40] analyzed the 2D Williamson fluid flow over a stretching sheet considering the influence of nanosized particles, also characterized as Williamson nanofluid. They studied the Williamson nanofluid peristaltic flow in a curved channel, including compliant walls.

Homotopy analysis method (HAM) was suggested by Liao [41] in 1992 to solve highly nonlinear differential equations. This technique has an edge over the rest of the contemporary techniques. HAM is one of the best and simplest technique for obtaining the convergent series solution for weakly, as well as, highly nonlinear differential equations. This technique includes the concept of homotopy from topology. HAM is used for finding a convergent series solution with high nonlinearity. Homotopy discriminates itself from other methods in the following ways:


There have been many attempts in the literature to discuss the varied fluid problems utilizing the homotopy analysis method [42–49].

A literature review discloses that copious literature may be quoted in the case of 2D non-Newtonian flows. Less work is available on 3D geometries, and this group becomes narrower if we talk about 3D Williamson nanofluid flows. The subject matter of 3D MHD flow of Williamson nanofluid over a bidirectional stretched surface with second-order slip and double stratification is even more rarely discussed. The structure of this paper is as follows: In Section 2, we present the mathematical model. In Sections 3–5, we discuss the homotopic scheme in detail with zeroth and *m*th order solutions. In Section 6, we address the convergence analysis. In Section 7, we present the results and discuss their physical importance, and finally, we provide concluding remarks.

#### **2. Mathematical Modeling**

Here, we consider the steady 3D Williamson nano liquid flow with velocities of *Uw* = *ax* in the *x*-direction and *Vw* = *by* in the *y*-direction, respectively, over a bidirectional extended sheet. While *a* and *b* are constants (Figure 1), concentration buoyancy force and thermal are used by the fluid with double stratification phenomena to study heat and mass transfers.

**Figure 1.** Schematic flow diagram.

The following are the governing boundary layer equations:

$$
u\_x + v\_y + w\_z = 0\tag{1}$$

$$\tau u u\_x + v u\_y + w u\_z = \nu u\_{zz} + \sqrt{2} \Gamma \nu u\_z u\_{zz} - \frac{\sigma B\_0^2}{\rho} u + \begin{bmatrix} g \left[ a\_1 (T - T\_{\text{co}}) + a\_2 (T - T\_{\text{co}})^2 \right] + \\ g \left[ a\_3 (\text{C} - \text{C}\_{\text{co}}) + a\_4 (\text{C} - \text{C}\_{\text{co}})^2 \right] \end{bmatrix} \tag{2}$$

$$
\Delta v v\_x + \upsilon v\_y + \upsilon w\_z = \nu v\_{zz} + \sqrt{2} \Gamma \nu v\_z v\_{zz} - \frac{\sigma B\_0^2}{\rho} v \tag{3}
$$

$$q + \lambda\_E (q\_t + V\nabla q - q\_\cdot \nabla V + (\nabla\_\cdot V)\mathbf{q}) = -\nabla(\mathbf{k}\mathbf{T})\tag{4}$$

$$J + \lambda\_{\mathbb{C}} (\mathbb{J}\_t + V.\nabla f - f.\nabla V + (\nabla.\mathbb{V})\mathbb{J}) = -D\_{\mathbb{B}} \mathbb{V} \mathbb{C} \tag{5}$$

$$q + \lambda \mathbb{E} \left( V.\nabla q - q.\nabla V \right) = -\nabla (kT) \tag{6}$$

$$J + \lambda\_{\mathbb{C}}(V.\nabla I - f.\nabla V) = -D\_{\mathbb{B}}\nabla \mathbb{C} \tag{7}$$

under the supervision of above-mentioned consideration and the impression of thermophoresis and Brownian-motion, Equations (6) and (7) takes the form:

$$
\mu T\_x + vT\_y + wT\_z + \lambda\_E \phi\_E = \frac{1}{\rho c\_p} \frac{\partial}{\partial z} (kT\_z) + \tau [D\_B \mathbb{C}\_z T\_z + \frac{D\_T}{T\_\infty} (T\_z)^2 \tag{8}
$$

$$
\mu \mathbf{C}\_x + \nu \mathbf{C}\_y + \nu \mathbf{C}\_z + \lambda\_C \phi\_\mathbf{C} = D\_B \mathbf{C}\_{zz} + \frac{D\_T}{T\_\infty} T\_{zz} \tag{9}
$$

where

$$\phi\_E = u^2 T\_{xx} + v^2 T\_{yy} + w T\_{zz} + 2uwT\_{xy} + 2uwT\_{xz} + 2vwT\_{yz} + \begin{bmatrix} (uu\_x + vu\_y + wu\_z)T\_x + \\ (uv\_x + vv\_y + wv\_z)T\_y + \\ (uu\_x + vw\_y + ww\_z)T\_z \end{bmatrix} \tag{10}$$

$$\phi\_{\mathbb{C}} = u^2 \mathbb{C}\_{xx} + v^2 \mathbb{C}\_{yy} + w \mathbb{C}\_{zz} + 2uv \mathbb{C}\_{xy} + 2uw \mathbb{C}\_{xz} + 2vw \mathbb{C}\_{yz} + \begin{bmatrix} (uu\_x + vu\_y + wu\_z)\mathbb{C}\_x + \\ (uv\_x + vv\_y + wv\_z)\mathbb{C}\_y + \\ (uv\_x + vw\_y + ww\_z)\mathbb{C}\_z \end{bmatrix} \tag{11}$$

following boundary conditions supports the above-mentioned system of equations:

$$\begin{array}{llll} u = \mathsf{U}\_{\mathsf{w}} + \mathsf{U}\_{\mathrm{slip}} & v = V\_{\mathsf{w}} + V\_{\mathrm{slip}} & w = 0\\ T = T\_{\mathsf{w}} = T\_{0} + d\_{1} \mathsf{x} & \mathsf{C} = \mathsf{C}\_{\mathsf{w}} = \mathsf{C}\_{0} + d\_{2} \mathsf{x} & \text{at} & z = 0\\ u \to 0 \; \; v \to 0 & T \to T\_{\mathrm{co}} = T\_{0} + e\_{1} \mathsf{x} & \mathsf{C} \to \mathsf{C}\_{\mathsf{co}} = \mathsf{C}\_{0} + e\_{2} \mathsf{x} & \text{as} & z \to \infty \end{array} \tag{12}$$

here

$$
\Delta I\_{\rm slip} = \frac{2}{3} \left( \frac{3 - \alpha l^3}{\alpha} - \frac{3}{2} \frac{1 - l^2}{K\_n} \right) \Lambda u\_{\overline{z}} - \frac{1}{4} \Big( l^4 + \frac{2}{K\_n^2} (1 - l^2) \Big) \Lambda^2 u\_{\overline{z}\overline{z}} = A u\_{\overline{z}} + B u\_{\overline{z}\overline{z}} \tag{13}
$$

$$\Delta V\_{slip} = \frac{2}{3} \left( \frac{3 - \alpha l^3}{\alpha} - \frac{3}{2} \frac{1 - l^2}{K\_n} \right) \Lambda v\_z - \frac{1}{4} \left( l^4 + \frac{2}{K\_n^2} (1 - l^2) \right) \Lambda^2 v\_{zz} = \mathcal{C}v\_z + Dv\_{zz} \tag{14}$$

where *A* = <sup>2</sup> 3 - 3−α*l* 3 <sup>α</sup> <sup>−</sup> <sup>3</sup> 2 1−*l* 2 *Kn* <sup>Λ</sup>, *<sup>B</sup>* = <sup>−</sup><sup>1</sup> 4 *l* <sup>4</sup> + <sup>2</sup> *K*2 *n* (1 − *l* 2) Λ2, *C* = <sup>2</sup> 3 - 3−α*l* 3 <sup>α</sup> <sup>−</sup> <sup>3</sup> 2 1−*l* 2 *Kn* Λ, *D* = <sup>4</sup> + <sup>2</sup> 2) <sup>Λ</sup>2, *<sup>l</sup>* = min <sup>1</sup>

−1 4 *l K*2 *n* (1 − *l Kn* , 1 , α describes momentum accommodation coefficient and varies from 0 ≤ α ≤ 1 *Kn* denotes Knudsen number and Λ denotes molecular mean free path. On the basis of definition of *l*, we found for any particular estimates of *Kn* we own 1 ≥ *l* ≥ 0. The molecular mean free path is always positive. Therefore, we know that B, D < 0 and C and A are positive.

To solve Equations (1),(3) and (8),(9), following similarity transformations are introduced:

$$\begin{array}{ll} u = a \text{x} f \, (\eta) & v = a \text{y} \mathbf{y} \, (\eta) & w = -\sqrt{a \nu} (f(\eta) + \mathbf{g}(\eta)) \\ \theta(\eta) = \frac{T - T\_{\text{cv}}}{T\_w - T\_{\text{cv}}} & \phi(\eta) = \frac{\mathbb{C} - \mathbb{C}\_{\text{cv}}}{\mathbb{C}\_w - \mathbb{C}\_{\text{cv}}} & \eta = \sqrt{\frac{\mathbb{A}}{\nu}} z \end{array} \tag{15}$$

here, *f*, *g*, θ and φ are the non-dimensional form for both velocities, temperature and the concentration.

Condition for incompressibility is self-satisfied and Equations (2),(3) and (8),(9) reduce to:

$$f^{\prime\prime\prime} - f^2 + (f+g)f^{\prime\prime} + \mathcal{We}f^{\prime\prime}f^{\prime\prime\prime} + \lambda (1 + \beta\_2 \theta)\theta + \lambda \mathcal{N}r(1 + \beta\_3 \phi)\phi - \mathcal{H}af \prime = 0 \tag{16}$$

$$\mathbf{g}^{\prime\prime\prime} - \mathbf{g}^{\prime} + (f + \mathbf{g})\mathbf{g}^{\prime\prime} + \mathcal{W}\mathbf{c} \mathbf{g}^{\prime\prime} \mathbf{g}^{\prime\prime\prime} - \mathcal{H}\mathbf{a} \mathbf{g} = \mathbf{0} \tag{17}$$

$$\begin{aligned} &(1+\varepsilon\theta)\theta''+\varepsilon\theta'^2+\text{Pr}N\_b\theta\nu\phi'+\text{Pr}N\_l\theta'^2-\text{Pr}f(\mathbb{S}\_1+\theta)+\\ &\text{Pr}(f+\mathfrak{g})\theta\nu-\delta\mathfrak{p}\text{Pr}\left(\begin{array}{c}(f+\mathfrak{g})^2\theta''-2f\nu\theta\nu(f+\mathfrak{g})+\left(f\tau^2-f''(f+\mathfrak{g})\right)(\mathbb{S}\_1+\theta)\\+(f+\mathfrak{g})(f\nu+\mathfrak{g}\nu)\theta\nu\end{array}\right)=0\end{aligned} \tag{18}$$

$$\begin{aligned} &\phi''' + \frac{N\_l}{N\_b} \theta'' - \text{PrL} \text{f} \prime (\text{S}\_2 + \phi) + \text{PrL} \varepsilon (f + \text{g}) \phi \prime - \\ &\text{PrL} \text{c} \delta\_b \bigg( \begin{array}{l} (f + \text{g})^2 \phi'' - 2f \tau (f + \text{g}) \phi \prime + \Big( f \tau^2 - f^{\prime \prime} (f + \text{g}) \Big) (\text{S}\_2 + \phi) + \\ (f + \text{g}) (f \prime + \text{g} \prime) \phi \prime \end{array} \tag{19}$$

and boundary conditions hold the form

$$\begin{cases} f(0) = 0 \text{ } f\prime(0) = 1 + \gamma\_1 f\prime\prime(0) + \gamma\_2 f\prime\prime(0) \text{ } g(0) = 0\\ \text{g}\prime(0) = 1 + \gamma\_3 \text{g}\prime\prime(0) + \gamma\_4 \text{g}\prime\prime(0) \text{ } \theta(0) = 1 - \text{S}\_1 \text{ } \phi(0) = 1 - \text{S}\_2\\ \text{f}\prime(\infty) \to 0 \text{ } \text{g}\prime(\infty) \to 0 \text{ } \theta(\infty) = 0 \text{ } \phi(\infty) = 0 & \text{as} \quad z \to \infty \end{cases} \tag{20}$$

where the parameters given above are defined as follows:

$$\begin{cases} \lambda = \frac{Gr\_{\overline{v}}}{Rv\_{x}^{2}} \quad Gr\_{x} = \frac{\xi \beta \Gamma(T\_{w} - T\_{w}) \mathbf{x}^{3}}{v^{2}} \operatorname{Re}\_{x} = \frac{\mu\_{0} \mathbf{x}}{v} \operatorname{N} \mathbf{v} = \frac{a\_{3} (\mathsf{C}\_{w} - \mathsf{C}\_{0})}{a\_{1} (\mathsf{T}\_{w} - \mathsf{T}\_{0})} \operatorname{S}\_{1} = \frac{c\_{1}}{d\_{1}}\\ \boldsymbol{\gamma}\_{1} = A \sqrt{\frac{a}{v}} \quad \boldsymbol{\gamma}\_{2} = B \frac{a}{v} \ \boldsymbol{\gamma}\_{3} = \mathbf{C} \sqrt{\frac{a}{v}} \ \boldsymbol{\gamma}\_{4} = D \frac{a}{v} \ \boldsymbol{\gamma}\_{2} = \frac{a\_{2}}{a\_{1}} (T\_{w} - T\_{0}) \ \boldsymbol{\beta} = \frac{b}{a}\\ \boldsymbol{\beta}\_{3} = \frac{a\_{3}}{a\_{3}} (\mathsf{C}\_{w} - \mathsf{C}\_{0}) \quad a\_{1} = \boldsymbol{\beta}\_{T} \ \mathsf{W} \boldsymbol{e} = \mathsf{U}\_{w} \boldsymbol{\Gamma} \sqrt{\frac{a\_{0}}{v}} \ \operatorname{Pr} = \frac{\mu\_{T}}{\bar{k}} \ \operatorname{Ha} = \frac{\sigma B\_{0}}{\rho a} \ \mathbf{S}\_{2} = \frac{c\_{2}}{d\_{2}}\\ \boldsymbol{L} = \frac{a}{B} \ \operatorname{N} = \frac{\tau D\_{T} (T\_{w} - T\_{0})}{T\_{w} \underline{v}} \ \operatorname{N}b = \frac{\tau D\_{\overline{b}} (\mathsf{C}\_{w} - \mathsf{C}\_{0})}{T\_{w} \underline{v}} \ \delta\_{\overline$$

*Cf x* is the coefficients of Skin friction in x- and *Cf y* in the y-direction are represented as follows:

$$\mathbb{C}\_{fx} = \frac{\tau\_{wx}}{\rho l I\_w^2} \quad \mathbb{C}\_{fy} = \frac{\tau\_{wy}}{\rho l I\_w^2} \tag{22}$$

$$\text{where } \tau\_{wx}|\_{z=0} = u\_z + \frac{\Gamma}{\sqrt{2}} (\mathbf{u}\_z)^2 \text{ and } \tau\_{wy}|\_{z=0} = v\_z + \frac{\Gamma}{\sqrt{2}} (\mathbf{v}\_z)^2 \tag{23}$$

Coefficients of Skin friction in dimensionless forms are:

$$\begin{array}{l} \mathbb{C}\_{f\mathbf{x}} \mathrm{Re}^{\frac{1}{2}} = \left[ f^{\prime\prime} + \frac{\mathrm{Wc}}{2} (f^{\prime\prime})^2 \right]\_{\eta=0} \\ \mathbb{C}\_{f\mathbf{y}} \mathrm{Re}^{\frac{1}{2}} = \left[ \mathbf{g}^{\prime\prime} + \frac{\mathrm{Wc}}{2} (\mathbf{g}^{\prime\prime})^2 \right]\_{\eta=0} \end{array} \tag{24}$$
