**5.** *m***th-Order Deformation Problems**

Here, we have

$$\mathbb{E}L\_f[f\_m(\eta) - \chi\_m f\_{m-1}(\eta)] = \hbar\_f \mathbb{R}\_f^n(\eta) \tag{36}$$

$$L\_{\mathfrak{F}}[\mathbf{g}\_m(\eta) - \chi\_m \mathbf{g}\_{m-1}(\eta)] = \hbar\_{\mathfrak{F}} R\_{\mathfrak{F}}''(\eta) \tag{37}$$

$$L\_0[\theta\_m(\eta) - \chi\_m \theta\_{m-1}(\eta)] = \hbar\_0 R\_0^n(\eta) \tag{38}$$

$$L\_{\phi}[\phi\_{\mathfrak{m}}(\eta) - \chi\_{\mathfrak{m}}\phi\_{\mathfrak{m}-1}(\eta)] = \hbar\_{\phi}R\_{\phi}^{\mathfrak{n}}(\eta) \tag{39}$$

$$\begin{aligned} R\_f^m(\eta) &= f\_{m-1}^{\prime\prime} - f\_{m-1}^{\prime2} + \sum\_{k=0}^{m-1} \left( f\_{m-1-k} + g\_{m-1-k} \right) f\_k^{\prime\prime} + \mathcal{W}e \sum\_{k=0}^{m-1} f\_{m-1-k}^{\prime\prime} f\_k^{\prime\prime} + \lambda \theta\_{m-1} + \\ \beta\_2 \lambda \sum\_{k=0}^{m-1} \theta\_{m-1-k} \theta\_k &+ \lambda Nr \phi\_{m-1} + \lambda Nr \beta\_3 \sum\_{k=0}^{m-1} \phi\_{m-1-k} \phi\_k - Haf\_{m-1}^{\prime} \end{aligned} \tag{40}$$

$$R\_{\mathcal{g}}^{m}(\eta) = \mathcal{g}\_{m-1}^{\prime} - \mathcal{g}\_{m-1}^{\prime 2} + \sum\_{k=0}^{m-1} (f\_{m-1-k} + g\_{m-1-k}) \mathcal{g}\_{k}^{\prime \prime} + \mathcal{N}\varepsilon \sum\_{k=0}^{m-1} \mathcal{g}\_{m-1-k}^{\prime \prime} \mathcal{g}\_{k}^{\prime \prime \prime} - H \mathcal{g}\_{m-1}^{\prime} \tag{41}$$

*R<sup>m</sup>* <sup>θ</sup> (η) = <sup>θ</sup>-- *<sup>m</sup>*−<sup>1</sup> <sup>+</sup> <sup>ε</sup> *m*(−1 *k*=0 θ*m*−1−*k*θ-- *<sup>k</sup>* + <sup>ε</sup> *m*(−1 *k*=0 θ- *m*−1−*k* θ- *<sup>k</sup>* + Pr*Nb m*(−1 *k*=0 θ- *m*−1−*k* φ- *<sup>k</sup>* + Pr*Nt m*(−1 *k*=0 θ- *m*−1−*k* θ- *k* − Pr(*S*<sup>1</sup> + θ)*f*- *<sup>m</sup>*−<sup>1</sup> <sup>+</sup> Pr*m*(−<sup>1</sup> *k*=0 (*fm*−1−*k*θ- *<sup>k</sup>* + *gm*−1−*k*θ- *k* ) <sup>−</sup> <sup>δ</sup>*t*Pr*m*(−<sup>1</sup> *k*=0 *fm*−1−*<sup>k</sup>* ( *k l*=0 *fk*−1θ-- <sup>1</sup> − <sup>δ</sup>*t*Pr*m*(−<sup>1</sup> *k*=0 *gm*−1−*<sup>k</sup>* ( *k l*=0 *gk*−1θ-- <sup>1</sup> <sup>−</sup> <sup>2</sup>δ*t*Pr*m*(−<sup>1</sup> *k*=0 *fm*−1−*<sup>k</sup>* ( *k l*=0 *gk*−1θ-- <sup>1</sup> − −2δ*t*Pr*m*(−<sup>1</sup> *k*=0 *fm*−1−*<sup>k</sup>* ( *k l*=0 *f*- *k*−1 θ- 1− <sup>2</sup>δ*t*Pr*m*(−<sup>1</sup> *k*=0 *gm*−1−*<sup>k</sup>* ( *k l*=0 *f*- *k*−1 θ- <sup>1</sup> <sup>−</sup> *<sup>S</sup>*1δ*t*Pr*m*(−<sup>1</sup> *k*=0 *f*- *<sup>m</sup>*−*k*−<sup>1</sup> *<sup>f</sup>*- *<sup>k</sup>*−<sup>1</sup> <sup>+</sup> *<sup>S</sup>*1δ*t*Pr*m*(−<sup>1</sup> *k*=0 *fm*−1−*<sup>k</sup> f* -- *<sup>k</sup>* + *<sup>S</sup>*1δ*t*Pr*m*(−<sup>1</sup> *k*=0 *gm*−1−*<sup>k</sup> f* -- *<sup>k</sup>* <sup>−</sup> <sup>δ</sup>*t*Pr*m*(−<sup>1</sup> *k*=0 *f*- *m*−*k*−1 ( *k l*=0 *f*- *k*−1 <sup>θ</sup><sup>1</sup> <sup>+</sup> <sup>δ</sup>*t*Pr*m*(−<sup>1</sup> *k*=0 *fm*−1−*<sup>k</sup>* ( *k l*=0 *f* -- *k*−1 θ1+ <sup>δ</sup>*t*Pr*m*(−<sup>1</sup> *k*=0 *gm*−1−*<sup>k</sup>* ( *k l*=0 *f* -- *k*−1 <sup>θ</sup><sup>1</sup> <sup>+</sup> <sup>δ</sup>*t*Pr*m*(−<sup>1</sup> *k*=0 *fm*−1−*<sup>k</sup>* ( *k l*=0 *f*- *k*−1 θ- <sup>1</sup> <sup>+</sup> <sup>δ</sup>*t*Pr*m*(−<sup>1</sup> *k*=0 *fm*−1−*<sup>k</sup>* ( *k l*=0 *g*- *k*−1 θ- 1 <sup>+</sup>δ*t*Pr*m*(−<sup>1</sup> *k*=0 *gm*−1−*<sup>k</sup>* ( *k l*=0 *f*- *k*−1 θ- <sup>1</sup> <sup>+</sup> <sup>δ</sup>*t*Pr*m*(−<sup>1</sup> *k*=0 *gm*−1−*<sup>k</sup>* ( *k l*=0 *g*- *k*−1 θ- 1 (42)

*Rm* <sup>φ</sup> (η) = <sup>φ</sup>-- *<sup>m</sup>*−<sup>1</sup> <sup>+</sup> *Nt Nb m*(−1 *k*=0 θ- *m*−1−*k* θ- *<sup>k</sup>* − PrLeS2 *f*- *<sup>m</sup>*−<sup>1</sup> <sup>−</sup> PrLe*m*(−<sup>1</sup> *k*=0 *f*- *m*−1−*k* θ*<sup>k</sup> f*- *<sup>m</sup>*−<sup>1</sup> <sup>+</sup> PrLe*m*(−<sup>1</sup> *k*=0 *fm*−1−*k*φ- *k* + PrLe*m*(−<sup>1</sup> *k*=0 *gm*−1−*k*φ- *<sup>k</sup>* − PrLeδ*<sup>c</sup> m*(−1 *k*=0 *fm*−1−*<sup>k</sup>* ( *k l*=0 *fk*−1φ-- <sup>1</sup> − PrLeδ*<sup>c</sup> m*(−1 *k*=0 *gm*−1−*<sup>k</sup>* ( *k l*=0 *gk*−1φ-- <sup>1</sup> − <sup>2</sup>δ*c*PrLe*m*(−<sup>1</sup> *k*=0 *fm*−1−*<sup>k</sup>* ( *k l*=0 *gk*−1φ-- <sup>1</sup> <sup>−</sup> <sup>2</sup>δ*c*PrLe*m*(−<sup>1</sup> *k*=0 *fm*−1−*<sup>k</sup>* ( *k l*=0 *f*- *k*−1 φ- 1− <sup>2</sup>δ*c*PrLe*m*(−<sup>1</sup> *k*=0 *gm*−1−*<sup>k</sup>* ( *k l*=0 *f*- *k*−1 φ- <sup>1</sup> <sup>+</sup> *<sup>S</sup>*2δ*c*PrLe*m*(−<sup>1</sup> *k*=0 *f*- *<sup>m</sup>*−*k*−<sup>1</sup> *<sup>f</sup>*- *<sup>k</sup>*−<sup>1</sup> <sup>+</sup> *<sup>S</sup>*2δ*c*PrLe*m*(−<sup>1</sup> *k*=0 *fm*−1−*<sup>k</sup> f* -- *<sup>k</sup>* + *<sup>S</sup>*2δ*c*PrLe*m*(−<sup>1</sup> *k*=0 *gm*−1−*<sup>k</sup> f* -- *<sup>k</sup>* − PrLeδ*<sup>c</sup> m*(−1 *k*=0 *f*- *m*−*k*−1 ( *k l*=0 *f*- *k*−1 φ1+ <sup>δ</sup>*c*PrLe*m*(−<sup>1</sup> *k*=0 *fm*−1−*<sup>k</sup>* ( *k l*=0 *f* -- *k*−1 <sup>φ</sup><sup>1</sup> <sup>+</sup> <sup>δ</sup>*c*PrLe*m*(−<sup>1</sup> *k*=0 *gm*−1−*<sup>k</sup>* ( *k l*=0 *f* -- *k*−1 <sup>φ</sup><sup>1</sup> <sup>+</sup> <sup>δ</sup>*c*PrLe*m*(−<sup>1</sup> *k*=0 *fm*−1−*<sup>k</sup>* ( *k l*=0 *f*- *<sup>k</sup>*−1φ- 1+ <sup>δ</sup>*c*PrLe*m*(−<sup>1</sup> *k*=0 *fm*−1−*<sup>k</sup>* ( *k l*=0 *g*- *k*−1 φ- <sup>1</sup> <sup>+</sup> <sup>δ</sup>*c*PrLe*m*(−<sup>1</sup> *k*=0 *gm*−1−*<sup>k</sup>* ( *k l*=0 *f*- *k*−1 φ- <sup>1</sup> <sup>+</sup> <sup>δ</sup>*c*PrLe*m*(−<sup>1</sup> *k*=0 *gm*−1−*<sup>k</sup>* ( *k l*=0 *g*- *k*−1 φ- 1 (43) χ*<sup>m</sup>* = 0, *<sup>m</sup>* <sup>≤</sup> <sup>1</sup> 1, *<sup>m</sup>* <sup>&</sup>gt; <sup>1</sup> (44)

the final solutions can be transcribed in the subsequent forms:

$$\begin{cases} f\_m(\eta) = f\_m^\*(\eta) + D\_1 + D\_2 \mathbf{e}^{\eta} + D\_3 \mathbf{e}^{-\eta} \\ g\_m(\eta) = g\_m^\*(\eta) + D\_4 + D\_5 \mathbf{e}^{\eta} + D\_6 \mathbf{e}^{-\eta} \\ \theta\_m(\eta) = \theta\_m^\*(\eta) + D\_7 \mathbf{e}^{\eta} + D\_8 \mathbf{e}^{-\eta} \\ \phi\_m(\eta) = \phi\_m^\*(\eta) + D\_9 \mathbf{e}^{\eta} + D\_{10} \mathbf{e}^{-\eta} \end{cases} \tag{45}$$

where *fm*, *gm*, θ*m*, and φ*<sup>m</sup>* symbolize the special solutions.

#### **6. Convergence Analysis**

HAM is used to obtain the solution of higher order nonlinear problems or those in series form. It gives several choices to control and modify the convergence region for the series solutions. Figure 2 represents the <sup>−</sup> curves behavior of all distributions. Characteristic parameters *<sup>f</sup>* , *g*, <sup>θ</sup> and <sup>φ</sup> have permissible ranges <sup>−</sup>1.6 <sup>≤</sup> *<sup>f</sup>* ≤ −0.4, <sup>−</sup>2.15 <sup>≤</sup> *<sup>g</sup>* ≤ −0.2, <sup>−</sup>2.75 <sup>≤</sup> <sup>θ</sup> ≤ −0.8 and <sup>−</sup>2.6 <sup>≤</sup> <sup>φ</sup> ≤ −0.6 when γ = 0.2, Le = 1, *Nt* = 0.2, *Nb* = 0.3, *P*r = 1.0, ε = 0.3, λ = 0.002, β1= β<sup>3</sup> = 0.2, β = 0.1 and M = 0.2 Table 1 represents the numerical results obtained for series solutions depicting the convergence of approximations up to the 25th order of approximations, that is, enough for series solution convergence. It can be verified that the graphical depiction in Figure 2 and the tabular results in Table 1 are in total consensus.

**Figure 2.** <sup>−</sup> curves for *<sup>f</sup>*, *<sup>g</sup>*, <sup>θ</sup> and <sup>φ</sup>.


**Table 1.** Convergence analysis of series solution using different order approximations when M = 0.2, γ<sup>1</sup> =0.2, γ<sup>2</sup> = 0.2, *Nb* = 0.3, *Nt* = 0.2, λ = 0.002, Le = 1.0, Pr = 1.0 δ<sup>t</sup> = 0.2, δ<sup>c</sup> = 0.2, ε =0.1, β = 0.1.

Table 2 was developed to validate the results obtained in the current model for skin friction in both directions by comparison with Malik et al. [50], and excellent harmony in both outcomes is achieved.

**Table 2.** Comparative estimates of *We* with Malik et al. [50] for skin friction along both directions in the limiting case, i.e., by considering the Hartmann number, second-order slip, temperature, and concentration profiles to zero.


#### **7. Results and Discussion**

In this section, we analyze the impact of appearing factors on particular distributions in Figures 3–21. The behavior on velocity profiles of β (ratio parameter) is described in Figures 3 and 4. It is noticed that contradictory behavior shown by both velocities (*f*-, *g*-) for an increasing rate of β. As β = *<sup>b</sup> <sup>a</sup>* , *a* was smaller for higher values of β, which specified a decreasing velocity rate along the *x*-direction, or *b* with higher values specified an increasing rate along the *y*-direction. In Figures 5 and 6, δ*<sup>t</sup>* and δ*<sup>c</sup>* illustrate the influence of thermal relaxation and the concentration relaxation factor on temperature concentration and distributions. We found that both concentration and temperature fields associated with the thickness of the boundary layers were the functions of decreasing δ*<sup>c</sup>* and δ*t*, respectively. Furthermore, δ*<sup>c</sup>* = 0 and δ*<sup>t</sup>* = 0 existing model will transform into classical laws of Fick's and Fourier's respectively. The influence of thermal conductivity ε on the temperature distribution is described in Figure 7. For higher values of ε, an increasing rate for the thermal boundary layer is found, which in result increases the temperature distribution. In Figure 8 the impact of Lewis number Le on concentration field is described. The strength of Lewis number depends on smaller estimations of mass diffusivity than the thermal diffusivity, which shows that exhausted Brownian motion coefficient decreases nanoparticle concentration profile. Figure 9 d the influence of mixed convective factor λ on the velocity field (*g*-). Higher estimations of λ produce stronger buoyancy force, which indicates an increasing rate in the velocity field (*g*-). The behavior of Prandtl number Pr against temperature distribution is presented in Figure 10. It is inspected that heat diffusion is very slow from the heated surface for higher estimates of Pr than smaller estimations of Pr. Therefore, temperature decreases with increasing values of Pr. Figures 11 and 12 show the influence of *H*a (Hartmann number) on both velocity profiles (*f*-,*g*-). Retardation in the fluid motion is seen due to resistance effered by strong Lorentz force. This act finally points out the decreasing rate on both velocity distributions. Figures 13 and 14 demonstrate the impact of Brownian motion factor *Nb* on concentration and temperature distribution. For the larger estimates of *Nb*, fluid temperature increases and rapidly reduces the deposition of particles far away from the fluid on the stretched surface. Due to which it increases and decreases concentration. The influence of the *Nt* on the concentration distribution is described in Figure 15. When the estimates of *Nt* are high and they are directly proportional to the temperature distribution, as a result it enhances the concentration distribution and its concentration associated with thickness of boundary layer. Figure 16 illustrates the influence of *Nt* on the temperature profile. When the values of *Nt* are higher, nanoparticles move from ambient fluid with higher temperature to the ambient fluid with lower temperature, and as a results temperature is higher in the boundary layer region. Finally, we identified the thickness of augmented thermal boundary layer. Figures 17 and 18 are drawn to display the influence of Williamson fluid parameter *We* on both velocity profiles. Increasing values of We decrease both velocities profiles. By increasing Williamson factor, relaxation time enhances. It causes to increase liquid viscosity, which results in decrease in the velocity profile. In Figure 19 the impact of stretching ratio factor β and mixed convective parameter λ on coefficient of Skin friction along *x*-direction is displayed. Which d that the coefficient of Skin friction shpws increasing behavior versus β and λ. Similar tendency can be observed in Hartmann number *Ha* and mixed convective parameter λ, against coefficient of Skin friction in *x*-direction, as displayed in Figure 20. Analysis of the impact of λ and *Nr* on Skin friction is described in Figure 21. It is noted that a thinner boundary layer is associated with larger λ, which result in higher velocity gradient near the wall. That's why Skin friction reduces against λ.

**Figure 4.** Impact of β against *g'*.

**Figure 8.** Impact of *L*e against φ.

**Figure 15.** Impact of *Nt* against φ.

**Figure 16.** Impact of *Nt* against θ.

**Figure 17.** Impact of *We* against *f'*.

**Figure 18.** Impact of We against *g'*.

**Figure 19.** Impact of λ and β against Skin friction.

**Figure 20.** Impact of λ and *Ha* against Skin friction.

**Figure 21.** Impact of *Nr* and λ against Skin friction.

#### **8. Concluding Remarks**

Three-dimensional Williamson nanofluid flow was investigated considering the Cattaneo–Christov heat flux model. The originality of the envisioned mathematical model was boosted by considering the influence of double stratification and second-order slip. HAM was applied to obtain the problem solution in series form. The salient outcomes of the problem are as follows:

• The stretching ratio parameter had an opposite impact on both velocities.


**Author Contributions:** Conceptualization, M.R. and A.L.; methodology, M.R.; software, S.K.; validation, S.Y., S.N., D.L. and Y.N; funding acquisition, Y.N.

**Funding:** This research was supported by Basic Science Research Program through the National Research Foundation of Korea(NRF) funded by the Ministry of Education(NRF-2017R1D1A3B03028309) and also supported by the Soonchunhyang University Research Fund.

**Conflicts of Interest:** The authors declare no conflicts of interest regarding this publication.
