Explicit Solutions of MHD Flow and Heat Transfer of Casson Fluid over an Exponentially Shrinking Sheet with Suction

Abstract

In this study, the magnetohydrodynamic (MHD) flow and heat transfer of a Casson fluid over an exponentially shrinking sheet with suction is investigated using the homotopy analysis method (HAM). Different from previous numerical methods and analytical techniques, we have obtained an explicit formula solution to the presented nonlinear problem. The explicit solutions of $f\left(\eta \right)$ and $\theta \left(\eta \right)$ are obtained and are valid in the whole domain. The changes in velocity and temperature profiles are studied in cases of different Casson fluid parameter $\gamma$, magnetic interaction parameter M, suction parameter s, and Prandtl number $Pr$. The convergent solutions are verified by comparison with the numerical results. In addition, the skin friction coefficient $C_f$ and local Nusselt number $N{u}_x$ are analyzed using the analytic formulas of $f^{″}\left(0\right)$ and ${\theta }^{\prime }\left(0\right)$, respectively. The analytical formulas help us intuitively analyze the influence of various parameters at the theoretical level. The effects of different physical quantities on $C_f$ and $N{u}_x$ are thoroughly investigated.

1. Introduction

Magnetohydrodynamics (MHD) refers to the study of the properties of magnetic fields and electrically conducting fluids. It focuses on the mutual interaction between fluid flow and magnetic fields. Thus, the fluids are required to be electrically conductive and non-magnetic, e.g., saltwater, plasmas, liquid metals, and electrolytes. MHD flows have found application in many branches of fluid mechanics [1,2,3,4]. The appearance of three-dimensional objects in wall-bounded MHD flows are characterized experimentally [5]. Reddy et al. [6] obtained numerical solutions of MHD flows in an electrically conductive fluid driven by a traveling magnetic field imposed at the end caps of a cylindrical annulus. Camobreco et al. [7] examined base flow influence in the context of the transition to turbulence in a quasi-two-dimensional MHD flow. With the development of MHD flow, many researchers have begun paying attention to boundary layer problems in MHD flow.

Most boundary layer problems considered in MHD flow are nonlinear. The flows are usually governed by one- or multiple-coupled nonlinear ordinary differential equations (ODEs). Thus, it is of great significance to develop effective methods to solve these nonlinear problems. Many researchers have conducted investigations on MHD boundary layer flow. Mukhopadhyay et al. [8] numerically studied the MHD boundary layer flow over a heated stretching sheet with variable viscosity. We investigate an MHD free-convection boundary layer flow saturated in a Darcian–Forchheimer porous medium over a vertical flat plate in the presence of suction/injection effect. This is conducted numerically by employing the routine bvp4c of the symbolic computer algebra software MATLAB [9]. The MHD boundary layer flow over a nonlinear stretching sheet is studied by a direct collocation method based on rational Legendre functions [10]. In addition to the MHD boundary layer problem of Newtonian fluids, attention has also increasingly shifted toward the certain flow problems of non-Newtonian fluids acting on the MHD boundary layer.

Recently, MHD boundary layer flow combined with Casson fluid has attracted the interest of many researchers. As a non-Newtonian fluid, Casson fluid is very common in nature and man-made products, such as jelly, soap, honey, and human blood. Some theoretical analyses of the related MHD Casson fluid flow have been conducted by several researchers through various numerical techniques [11,12,13,14]. Although there exists an abundance of numerical methods to deal with the nonlinear MHD boundary layer flow of a Casson fluid, the obtainment of an explicit analytical approximation of these nonlinear problems remains a challenging task.

We note that Nadeema et al. [15] employed the Adomian Decomposition Method (ADM) and numerically studied the MHD boundary layer flow of a Casson fluid induced by an exponentially shrinking sheet. In their work, the heat transfer of Casson fluid over the exponentially shrinking sheet is not considered. However, heat transfer is a crucial aspect applied in different engineering areas such as materials, energy, machinery, chemical industry, medicine, and other fields. Casson fluid is a typical non-Newtonian fluid, and it is of great significance to study its heat transfer characteristics in MHD boundary layer flow in practice. To the best of our knowledge, no explicit analytical solutions have been presented for the MHD flow and heat transfer of a Casson fluid over an exponentially shrinking sheet. This provides the inspiration behind us conducting this study.

In this paper, we apply the homotopy analysis method (HAM) [16,17,18,19] to the nonlinear MHD flow and heat transfer of a Casson fluid over an exponentially shrinking sheet with suction. The HAM method allows great freedom in the selection of proper base functions, auxiliary linear operators, initial guesses of unknowns, and convergence-control parameters. This method provides a straightforward way to ensure the convergence of solution series in nonlinear problems even for strongly nonlinear equations. Nowadays, it is widely applied to solve many strongly nonlinear problems in different areas [20,21,22,23,24,25,26,27]. The success of these applications demonstrates the great potential of the HAM method. In particular, the HAM method has been used to deal with many boundary layer flows, such as Blasius’ viscous flow [28], MHD Falkner–Skan flow of nano-fluids [29], Casson fluid flow with stretching sheet [30], and fluid flow over an exponentially stretching porous sheet [31]. These encourage us to employ the HAM method to the present boundary layer flow and further improve it.

We emphasize that explicit analytic formulas for the velocity and temperature distribution in a considered system do not currently exist, and which should also be valid in the whole domain. This motivates us to solve the problem. In the next section, we provide the mathematical model and formulas of the present boundary layer flow. The explicit analytic solutions of the current system are given in Section 3. In Section 4, we test the obtained explicit analytic solutions. The convergence and accuracy of solutions are evaluated in detail. Discussions concerning the effects of various physical parameters are investigated in Section 5. Summaries are presented in the last Section 6.

2. Mathematical Model

Consider the steady two-dimensional MHD boundary layer flow of a Casson fluid over an exponentially shrinking sheet. As illustrated in Figure 1, we assume that $U_w\left(x\right)$ denotes the velocity of the shrinking sheet, $V_m\left(x\right)$ is the variable wall mass transfer velocity, $T_w$ is the uniformly distributed sheet temperature, $T_{\infty }$ is the free stream temperature assumed to be constant, where

$T_w>T_{\infty }$ corresponds to the heated shrinking sheet. The fluid is electrically conductive with the uniform magnetic field applied normal to the shrinking sheet. Under the boundary layer approximations, the governing equations of continuity, motion, and energy equations are written as

$$\begin{array}{cc}\hfill & \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=\nu (1+\frac{1}{\gamma })\frac{{\partial }^{2}u}{\partial y^2}-\frac{\sigma B^2}{\rho }u,\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill & u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\alpha \frac{{\partial }^{2}T}{\partial y^2},\hfill \end{array}$$

(3)

where u and v are the corresponding velocities in x- and y-directions, respectively. $\nu =\mu /\rho$ is the kinematic fluid viscosity, $\rho$ is the fluid density, and $\mu$ is the Casson viscosity coefficient. $\gamma$ is the Casson fluid parameter. $\sigma$ is the electrical conductivity of the Casson fluid. $B=B_{0}exp(x/L)$ denotes the magnetic field, where $B_0$ is the constant magnetic field, and L is the characteristic length. $\alpha =\kappa /\left(\rho c_p\right)$ is the thermal diffusivity, where $\kappa$ is the fluid thermal conductivity and $c_p$ is the specific heat. The boundary conditions of governing Equations (1)–(3) are given by

$$\begin{array}{cc}\hfill & u=U_w=-U_{0}exp\left(\frac{x}{L}\right),\mathrm{at}y=0,\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill & v=V_w=V_{0}exp\left(\frac{x}{2L}\right),\mathrm{at}y=0,\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill & T=T_w=T_{\infty }+T_0,\mathrm{at}y=0,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & U=0,T=T_{\infty },\mathrm{as}y\to \infty \hfill \end{array}$$

(7)

where $U_0>0$ is shrinking constant, $V_0$ is a constant with $V_0<0$ for masss suction, $T_0$ is a constant measuring the uniform increase in temperature along the shrinking sheet.

The governing Equations (1)–(3) are nonlinear partial differential equations. In order to obtain the self-similar equations, we use the similarity transformation reported by Refs. [15,32], and thus Equations (1)–(3) are reduced to the nonlinear ordinary differential equations

$$\begin{array}{cc}\hfill & (1+1/\gamma )f^{‴}\left(\eta \right)-M{f}^{\prime }\left(\eta \right)+f\left(\eta \right)f^{″}\left(\eta \right)-2f^{\prime 2}\left(\eta \right)=0,\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & {\theta }^{″}\left(\eta \right)+Prf\left(\eta \right){\theta }^{\prime }\left(\eta \right)=0,\hfill \end{array}$$

(9)

with boundary conditions

$$\begin{array}{c}\hfill f\left(0\right)=s,f^{\prime }\left(0\right)=-1,f^{\prime }(+\infty )=0,\theta \left(0\right)=1,\theta (+\infty )=0,\end{array}$$

(10)

where the prime denotes the differentiation with respect to $\eta$, given by $\eta =y\sqrt{\frac{U_0}{2\nu L}}{e}^{\frac{x}{2L}}$, $f\left(\eta \right)$ is related to the stream function $\psi$ by $f\left(\eta \right)=\psi /\left(\sqrt{2\nu L{U}_0}{e}^{x/2L}\right)$, $\theta \left(\eta \right)=(T-T_{\infty })/(T_w-T_{\infty })$, $M=2L\sigma B^2/\rho U_0$ is the magnetic interaction parameter, $Pr=\mu c_p/\kappa$ is the Prandtl number, and $s=-V_0\sqrt{2L/\nu U_0}$ is the suction parameter.

The related skin friction coefficient $C_f$ and local Nusselt number $N{u}_x$ are given by using the similarity variables

$$\begin{array}{c}\hfill \sqrt{R{e}_x}{C}_f=(1+1/\gamma )f^{″}\left(0\right),\end{array}$$

(11)

$$\begin{array}{c}\hfill \frac{N{u}_x}{\sqrt{R{e}_x}}=-{\theta }^{\prime }\left(0\right),\end{array}$$

(12)

where $R{e}_x$ is the local Reynolds number.

3. Explicit Solutions

For the current problem, considering the convenience of calculation and obtaining the convergent solution [28], we first introduce a spatial scale factor $\lambda =3$ to transform Equations (8) and (9) into

$$\begin{array}{cc}\hfill & (1+\frac{1}{\gamma }){\lambda }^2{F}^{‴}\left(\xi \right)-M{F}^{\prime }\left(\xi \right)+\lambda F\left(\xi \right)F^{″}\left(\xi \right)-2\lambda F^{\prime 2}\left(\xi \right)=0,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill & {\mathsf{\Theta }}^{″}\left(\xi \right)+PrF\left(\xi \right){\mathsf{\Theta }}^{\prime }\left(\xi \right)=0,\hfill \end{array}$$

(14)

by the transformation $\xi =\lambda \eta$, $F\left(\xi \right)=f\left(\eta \right)$, and $\mathsf{\Theta }\left(\xi \right)=\theta \left(\eta \right)$. The corresponding boundary conditions (10) become

$$\begin{array}{c}\hfill F\left(0\right)=s,F^{\prime }\left(0\right)=-\frac{1}{\lambda },F^{\prime }(+\infty )=0,\mathsf{\Theta }\left(0\right)=1,\mathsf{\Theta }(+\infty )=0.\end{array}$$

(15)

In the framework of HAM, a nonlinear problem is transformed into an infinite number of linear sub-problems. According to Equations (13)–(15), $F\left(\xi \right)$ and $\mathsf{\Theta }\left(\xi \right)$ should hold the following “solution expressions” form

$$\begin{array}{cc}\hfill F\left(\xi \right)& =\sum _{m=0}^{+\infty }{b}_{m,0}^0+\sum _{m=0}^{+\infty }\sum _{n=1}^{+\infty }\sum _{k=0}^{+\infty }{b}_{m,n}^k{\xi }^{k}exp(-n\xi ),\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill \mathsf{\Theta }\left(\xi \right)& =\sum _{m=0}^{+\infty }\sum _{n=1}^{+\infty }\sum _{l=0}^{+\infty }{d}_{m,n}^l{\xi }^{l}exp(-n\xi ),\hfill \end{array}$$

(17)

where $b_{m,n}^k$ and $d_{m,n}^l$ are constant coefficients to be determined by HAM. The solution expressions (16)–(17) guide us to choose the following auxiliary linear operator

$$\begin{array}{c}\hfill {\pounds }_F\left[\mathsf{\Phi }(\xi ;q)\right]=\frac{d^3\mathsf{\Phi }}{d{\xi }^3}+\frac{d^2\mathsf{\Phi }}{d{\xi }^2},\end{array}$$

(18)

$$\begin{array}{c}\hfill {\pounds }_{\mathsf{\Theta }}\left[\mathsf{\Psi }(\xi ;q)\right]=\frac{d^2\mathsf{\Psi }}{d{\xi }^2}+\frac{d\mathsf{\Psi }}{d\xi },\end{array}$$

(19)

with the properties

$$\begin{array}{cc}\hfill & {\pounds }_F[C_1\xi +C_{2}exp(-\xi )+C_3]=0,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & {\pounds }_{\mathsf{\Theta }}[C_4+C_{5}exp(-\xi )]=0,\hfill \end{array}$$

(21)

where $C_1$, $C_2$, $C_3$, $C_4$ and $C_5$ are constants. Considering Equations (15)–(17), the following initial guesses are chosen

$$\begin{array}{ccc}\hfill F_0\left(\xi \right)& =& s-\frac{1}{\lambda }+\frac{1}{\lambda }{e}^{-\xi },\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill {\mathsf{\Theta }}_0\left(\xi \right)& =& e^{-\xi }.\hfill \end{array}$$

(23)

Then, according to Equations (13) and (14), we construct a family of zeroth-order deformation equations

$$\begin{array}{l}(1-q){\pounds }_F\left[\mathsf{\Phi }(\xi ;q)-F_0\left(\xi \right)\right]\\ =q{c}_0^f\left\{(1+\frac{1}{\gamma }){\lambda }^2\frac{{\partial }^3\mathsf{\Phi }(\xi ;q)}{\partial {\xi }^3}-M\frac{\partial \mathsf{\Phi }(\xi ;q)}{\partial \xi }+\lambda \mathsf{\Phi }(\xi ;q)\frac{{\partial }^2\mathsf{\Phi }(\xi ;q)}{\partial {\xi }^2}-2\lambda {\left[\frac{{\partial }^2\mathsf{\Phi }(\xi ;q)}{\partial {\xi }^2}\right]}^2\right\},\end{array}$$

(24)

$$\begin{array}{cc}\hfill & (1-q){\pounds }_{\mathsf{\Theta }}\left[\mathsf{\Psi }(\xi ;q)-{\mathsf{\Theta }}_0\left(\xi \right)\right]=q{c}_0^{\theta }\left\{\lambda \frac{{\partial }^2\mathsf{\Psi }(\xi ;q)}{\partial {\xi }^2}+Pr\mathsf{\Phi }(\xi ;q)\frac{\partial \mathsf{\Psi }(\xi ;q)}{\partial \xi }\right\},\hfill \end{array}$$

(25)

with boundary conditions

$$\begin{array}{cc}\hfill & \mathsf{\Phi }(0;q)=s,{\mathsf{\Phi }}^{\prime }(0;q)=-\frac{1}{\lambda },{\mathsf{\Phi }}^{\prime }(+\infty ;q)=0,\hfill \\ & \mathsf{\Psi }(0;q)=1,\mathsf{\Psi }(+\infty ;q)=0,\hfill \end{array}$$

(26)

where $q\in [0,1]$ is an embedding parameter, and $c_0^f$ and $c_0^{\theta }$ are two non-zero convergence-control parameters. When $q=0$, we have the solution

$$\mathsf{\Phi }(\xi ;0)=F_0\left(\xi \right),\mathsf{\Psi }(\xi ;0)={\mathsf{\Theta }}_0\left(\xi \right),$$

(27)

and at $q=1$, Equations (24)–(26) are the same as Equations (13)–(15), respectively, so that

$$\mathsf{\Phi }(\xi ;1)=F\left(\xi \right),\mathsf{\Psi }(\xi ;1)=\mathsf{\Theta }\left(\xi \right).$$

(28)

It is clear that the embedding parameter q plays the role of mapping. In other words, when q gradually increases from 0 to 1, this mapping ensures $\mathsf{\Phi }(\xi ;q)$ and $\mathsf{\Psi }(\xi ;q)$ continuously deform from the initial guesses $F_0\left(\xi \right)$ and ${\mathsf{\Theta }}_0\left(\xi \right)$ to the exact solutions $F\left(\xi \right)$ and $\mathsf{\Theta }\left(\xi \right)$, respectively.

For this direct mapping, we express $\mathsf{\Phi }(\xi ;q)$ and $\mathsf{\Psi }(\xi ;q)$ in Maclaurin series

$$\mathsf{\Phi }(\xi ;q)=\mathsf{\Phi }(\xi ;0)+\sum _{m=1}^{+\infty }\frac{q^m}{m!}\frac{{\partial }^m\mathsf{\Phi }(\xi ;q)}{\partial q^m}{|}_{q=0}=F_0\left(\xi \right)+\sum _{m=1}^{+\infty }{F}_m\left(\xi \right)q^m,$$

(29)

$$\mathsf{\Psi }(\xi ;q)=\mathsf{\Psi }(\xi ;0)+\sum _{m=1}^{+\infty }\frac{q^m}{m!}\frac{{\partial }^m\mathsf{\Psi }(\xi ;q)}{\partial q^m}{|}_{q=0}={\mathsf{\Theta }}_0\left(\xi \right)+\sum _{m=1}^{+\infty }{\mathsf{\Theta }}_m\left(\xi \right)q^m,$$

(30)

respectively. Note that the convergence of the above series is related to the two auxiliary linear operators ${\pounds }_F$, ${\pounds }_{\mathsf{\Theta }}$, and convergence-control parameters $c_0^f$, $c_0^{\theta }$. Fortunately, the HAM method allows us to choose them freely. This point is different from other analytic methods, due to the absence of physical meaning regarding the convergence-control parameter. If all of them are properly selected, we obtain via Equation (28) that

$$F\left(\xi \right)=F_0\left(\xi \right)+\sum _{m=1}^{+\infty }{F}_m\left(\xi \right)\approx F_0\left(\xi \right)+\sum _{m=1}^{\overline{M}}{F}_m\left(\xi \right),$$

(31)

$$\mathsf{\Theta }\left(\xi \right)={\mathsf{\Theta }}_0\left(\xi \right)+\sum _{m=1}^{+\infty }{\mathsf{\Theta }}_m\left(\xi \right)\approx {\mathsf{\Theta }}_0\left(\xi \right)+\sum _{m=1}^{\overline{M}}{\mathsf{\Theta }}_m\left(\xi \right),$$

(32)

where $\overline{M}$ is a sufficiently large truncation number and denotes $\overline{M}$th-order approximation. $F_m\left(\xi \right)$ and ${\mathsf{\Theta }}_m\left(\xi \right)$ ($m\ge 1$) are calculated from the corresponding governing equations

$$\begin{array}{cc}\hfill & {\pounds }_F[F_m\left(\xi \right)-{\chi }_n{F}_{m-1}\left(\xi \right)]\hfill \\ \hfill & =c_0^f\left\{(1+\frac{1}{\lambda }){\lambda }^2{F}_{m-1}^{‴}\left(\xi \right)-M{F}_{m-1}^{\prime }\left(\xi \right)+\sum _{i=0}^{m-1}{F}_i^{″}\left(\xi \right)F_{m-1-i}\left(\xi \right)-2\lambda \sum _{i=0}^{m-1}{F}_i^{\prime }\left(\xi \right)F_{m-1-i}^{\prime }\left(\xi \right)\right\}\hfill \\ \hfill & =c_0^f{\mathcal{R}}_{m-1}^F\left(\xi \right),\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill & {\pounds }_{\mathsf{\Theta }}[{\mathsf{\Theta }}_m\left(\xi \right)-{\chi }_m{\mathsf{\Theta }}_{m-1}\left(\xi \right)]=c_0^{\theta }\left\{\lambda {\mathsf{\Theta }}_{m-1}^{″}\left(\xi \right)+Pr\sum _{i=0}^{m-1}{F}_{m-1-i}\left(\xi \right){\mathsf{\Theta }}_i^{\prime }\left(\xi \right)\right\}=c_0^{\theta }{\mathcal{R}}_{m-1}^{\mathsf{\Theta }}\left(\xi \right),\hfill \end{array}$$

(34)

with the boundary conditions

$$\begin{array}{c}\hfill F_m\left(0\right)=F_m^{\prime }\left(0\right)=F_m^{\prime }(+\infty )=0,{\mathsf{\Theta }}_m\left(0\right)={\mathsf{\Theta }}_m(+\infty )=0,\end{array}$$

(35)

where the right-hand terms ${\mathcal{R}}_{m-1}^F\left(\xi \right)$ and ${\mathcal{R}}_{m-1}^{\mathsf{\Theta }}\left(\xi \right)$ are derived by substituting Equations (29) and (30) into Equations (24)–(26) and equalizing the like-powers of q. ${\chi }_m$ is defined by

$${\chi }_m=\left\{\begin{array}{cc}0,\hfill & \mathrm{when} m\le 1,\hfill \\ 1,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(36)

It should be emphasized that Equations (33) and (34) are linear. Thus, first considering the initial guesses (22)–(23), then $F_1\left(\xi \right)$, $F_2\left(\xi \right)$, …, $F_m\left(\xi \right)$ and ${\mathsf{\Theta }}_1\left(\xi \right)$, ${\mathsf{\Theta }}_2\left(\xi \right)$, …, ${\mathsf{\Theta }}_m\left(\xi \right)$ can be calculated by solving these linear equations. Interestingly, we can now find the structures of $F_m\left(\xi \right)$ and ${\mathsf{\Theta }}_m\left(\xi \right)$ under each mth-order approximation. Then the recurrence formulas regarding the coefficients $b_{m,n}^k$ of Equation (16) and $d_{m,n}^l$ of Equation (17) are deduced through a lengthy derivation process. For details, please refer to the process in [28]. Here, we list the explicit recurrence formulas of $b_{m,n}^k$ in Equation (16) and $d_{m,n}^l$ in Equation (17).

3.1. $b_{m,n}^k$ and $f\left(\eta \right)$

The coefficients $b_{m,n}^k$ in Equation (16) are as follows:

$$b_{0,0}^0=s-\frac{1}{\lambda },b_{0,1}^0=\frac{1}{\lambda },$$

(37)

and for $m\ge 1$, $0\le n\le m+1$ and $0\le k\le m-n+1$,

$$\begin{array}{cc}\hfill b_{m,0}^0=& {\chi }_m{b}_{m-1,0}^0-\sum _{q=0}^{m-1}{\mathsf{\Gamma }}_{m,1}^q{\mu }_{1,1}^q-\sum _{n=2}^{m+1}(n-1){\mathsf{\Gamma }}_{m,n}^0{\mu }_{n,0}^0-\sum _{n=2}^m\sum _{q=1}^{m-n+1}{\mathsf{\Gamma }}_{m,n}^q((n-1){\mu }_{n,0}^q-{\mu }_{n,1}^q),\hfill \end{array}$$

(38)

$$b_{m,0}^1=0,$$

(39)

$$\begin{array}{cc}\hfill b_{m,1}^0=& {\chi }_m{b}_{m-1,1}^0+\sum _{q=0}^{m-1}{\mathsf{\Gamma }}_{m,1}^q{\mu }_{1,1}^q+\sum _{n=2}^{m+1}n{\mathsf{\Gamma }}_{m,n}^0{\mu }_{n,0}^0+\sum _{n=2}^m\sum _{q=1}^{m-n+1}{\mathsf{\Gamma }}_{m,n}^q(n{\mu }_{n,0}^q-{\mu }_{n,1}^q),\hfill \end{array}$$

(40)

$$b_{m,1}^k={\chi }_m{b}_{m-1,1}^k+\sum _{q=k-1}^{m-1}{\mathsf{\Gamma }}_{m,1}^q{\mu }_{1,k}^q,1\le k\le m-1,$$

(41)

$$b_{m,1}^m={\mathsf{\Gamma }}_{m,1}^{m-1}{\mu }_{1,m}^{m-1},$$

(42)

$$\begin{array}{cc}\hfill & b_{m,n}^k={\chi }_m{b}_{m-1,n}^k-\sum _{q=k}^{m-n+1}{\mathsf{\Gamma }}_{m,n}^q{\mu }_{n,k}^q,0\le k\le m-n,2\le n\le m,\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill & b_{m,n}^{m-n+1}=-{\mathsf{\Gamma }}_{m,n}^{m-n+1}{\mu }_{n,m-n+1}^{m-n+1},2\le n\le m,\hfill \end{array}$$

(44)

$$b_{m,m+1}^0=-{\mathsf{\Gamma }}_{m,m+1}^0{\mu }_{m+1,0}^0,$$

(45)

for $m\ge 1$, please refer to Appendix A for ${\mathsf{\Gamma }}_{m,n}^q$ and ${\mu }_{n,k}^q$ mentioned in above formulas.

Therefore, owing to the transformation $\xi =\lambda \eta$, $f\left(\eta \right)=F\left(\xi \right)$, we obtain the explicit analytic solution of $f\left(\eta \right)$

$$\begin{array}{cc}\hfill f\left(\eta \right)=& \underset{\overline{M}\to +\infty }{lim}\sum _{m=0}^{\overline{M}}{f}_m\left(\eta \right)\hfill \\ \hfill =& \underset{\overline{M}\to +\infty }{lim}\sum _{m=0}^{\overline{M}}{b}_{m,0}^0+\underset{\overline{M}\to +\infty }{lim}\sum _{n=1}^{\overline{M}+1}exp(-n\lambda \eta )\times \left(\sum _{m=n-1}^{\overline{M}}\sum _{k=0}^{m-n+1}{b}_{m,n}^k{\left(\lambda \eta \right)}^k\right),\hfill \end{array}$$

(46)

3.2. $d_{m,n}^l$ and $\theta \left(\eta \right)$

The coefficients $d_{m,n}^l$ in Equation (17) are as follows:

$$\begin{array}{cc}\hfill & d_{0,1}^0=1,d_{m,0}^0=0,m\ge 1,\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill & d_{m,1}^0={\chi }_m{d}_{m-1,1}^0-\sum _{n=2}^{m+1}\sum _{q=0}^{m-n+1}{\mathsf{\Omega }}_{m,n}^q{w}_{n,0}^q,\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill & d_{m,1}^l={\chi }_m{d}_{m-1,1}^l+\sum _{q=l}^m{\mathsf{\Omega }}_{m,1}^{q-1}{w}_{1,l}^{q-1},1\le l\le m-1,\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill & d_{m,1}^m={\mathsf{\Omega }}_{m,1}^{m-1}{w}_{1,m}^{m-1},\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill & d_{m,n}^l={\chi }_m{d}_{m-1,n}^l+\sum _{q=l}^{m-n+1}{\mathsf{\Omega }}_{m,n}^q{w}_{n,l}^q,2\le n\le m,0\le l\le m-n,\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill & d_{m,n}^l=\sum _{q=l}^{m-n+1}{\mathsf{\Omega }}_{m,n}^q{w}_{n,l}^q,2\le n\le m,l=m-n+1,\hfill \end{array}$$

(52)

$$\begin{array}{cc}\hfill & d_{m,m+1}^0={\mathsf{\Omega }}_{m,m+1}^0{w}_{m+1,0}^0,\hfill \end{array}$$

(53)

for $m\ge 1$, the details of ${\mathsf{\Omega }}_{m,n}^q$ and $w_{n,l}^q$ are listed in Appendix B.

Similarly, the explicit expression of $\theta \left(\eta \right)$ is gained

$$\begin{array}{c}\hfill \theta \left(\eta \right)=\underset{\overline{M}\to +\infty }{lim}\sum _{m=1}^{\overline{M}}{\theta }_m\left(\eta \right)=\underset{\overline{M}\to +\infty }{lim}\sum _{n=1}^{\overline{M}+1}exp(-n\lambda \eta )\left(\sum _{m=n-1}^{\overline{M}}\sum _{l=0}^{m-n+1}{d}_{m,n}^l{\left(\lambda \eta \right)}^l\right).\end{array}$$

(54)

Using all of the above recurrence formulas, the coefficients $b_{m,n}^k$ and $d_{m,n}^k$ are easily calculated in turn. It is worth emphasizing that the series solutions (46) and (54) are convergent as long as the convergence-control parameters are properly selected. Thus, one can gain accurate results under different values of $\gamma$, M, and s.

In particular, the total averaged value of the squared residual error in the governing equations [19] is evaluated by substituting the $\overline{M}$th-order approximations (46) and (54) into the original governing Equations (8) and (9), respectively. The squared residual error clearly indicates the accuracy of the analytic approximations (46)–(54). It is crucial to guarantee the convergence of an approximation series. As Liao [19] reported, the series approximations contain the convergence-control parameters $c_0^f$ and $c_0^{\theta }$. Thus, the squared residual errors also contain $c_0^f$ and $c_0^{\theta }$. The proper values of $c_0^f$ and $c_0^{\theta }$ can always be found to guarantee the convergence of the homotopy series owing to the great freedom. Obviously, at the given order of approximation $\overline{M}$, the optimal approximation is defined by minimizing the squared residual error with the corresponding optimal convergence-control parameters $c_0^{f^{*}}$ and $c_0^{{\theta }^{*}}$, respectively.

4. Convergence Test

Physically, the values of $f^{″}\left(0\right)$ and ${\theta }^{\prime }\left(0\right)$ are of great significance. They are related to the skin friction coefficient $C_f$ (11) and the local Nusselt number $N{u}_x$ (12), respectively. In this section, we first give the analytical formulas of $f^{″}\left(0\right)$ and ${\theta }^{\prime }\left(0\right)$ and then test their convergence.

Through the use of Equations (46) and (54), the $\overline{M}$th-order approximation of $f^{″}\left(0\right)$ and ${\theta }^{\prime }\left(0\right)$ are derived

$$\begin{array}{cc}\hfill & f^{″}\left(0\right)=\sum _{n=1}^{\overline{M}+1}\sum _{m=n-1}^{\overline{M}}{\lambda }^2(n^2{b}_{m,n}^0-2n{\beta }_{m,n}^1{b}_{m,n}^1+2{\beta }_{m,n}^2{b}_{m.n}^2),\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill & {\theta }^{\prime }\left(0\right)=\sum _{n=1}^{\overline{M}+1}\sum _{m=n-1}^{\overline{M}}\lambda ({\Lambda }_{m,n}^1{d}_{m,n}^1-n{d}_{m.n}^0),\hfill \end{array}$$

(56)

where ${\beta }_{m,n}^k$ is defined by Equation (A14), and ${\Lambda }_{m,n}^l$ reads

$${\Lambda }_{m,n}^l=\left\{\begin{array}{cc}0,\hfill & m=n=0,l\ge 0,\hfill \\ 0,\hfill & m>0,n=0,l\ge 0,\hfill \\ 0,\hfill & n>m+1,\hfill \\ 0,\hfill & l>m-n+1\hfill \\ 1,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(57)

To reveal the accuracy and superiority of the series solution, let us first consider the case $\gamma =M=1$ with different values of the suction parameter s, which has the convergent series solution of $f^{″}\left(0\right)$ and ${\theta }^{\prime }\left(0\right)$. As shown in Figure 2, the convergent values of $f^{″}\left(0\right)$ and ${\theta }^{\prime }\left(0\right)$ obtained by HAM agree quite well with the numerical results [33]. Without the loss of generality, we further investigate the changes of skin friction and Nusselt number with various values of Casson fluid parameter $\gamma$, as shown in Figure 3. The convergent values of $f^{″}\left(0\right)$ and ${\theta }^{\prime }\left(0\right)$ are in good agreement with the corresponding numerical results in a region of $\gamma$. All of these cases indicate that the analytic formulas of $f^{″}\left(0\right)$ and ${\theta }^{\prime }\left(0\right)$ (55)–(56) are valid, so that one can obtain a sufficiently accurate approximation by means of the optimal convergence-control parameters. Moreover, according to Figure 2, it is seen that both the skin friction and the Nusselt number increase linearly with an increase in s. However, the influence of $\gamma$ is relatively weak, $C_f$, $f^{″}\left(0\right)$, $N{u}_x$ and $-{\theta }^{\prime }\left(0\right)$ increase slowly and nonlinearly with respect to $\gamma$.

Furthermore, we would like to emphasize two points. It is known that the HAM method bears superiority over other analytical/semi-analytical methods. Firstly, HAM can ensure the convergence accuracy of nonlinear problems, even in those with strong nonlinearity. A comparison of $f^{″}\left(0\right)$ and ${\theta }^{\prime }\left(0\right)$ with the corresponding numerical values for various values of M, $\gamma$, s is illustrated in

Table 1. Analytic approximations (HAM) of $f^{″}\left(0\right)$, ${\theta }^{\prime }\left(0\right)$ compared with numerical results for various M, $\gamma$ and s in case of $Pr=0.71$.

Parameters		${\mathit{f}}^{″}\left(0\right)$		${\mathit{\theta }}^{\prime }\left(0\right)$
		HAM	Num	HAM	Num
M	2	1.18455279	1.184553	−1.10536254	−1.105533
	3	1.47878566	1.478786	−1.14720967	−1.147279
	5	1.89590282	1.895903	−1.18877394	−1.188803
	10	2.60700267	2.607003	−1.23492995	−1.234942
$\gamma$	1	0.61356181	0.613026	−0.94909557	−0.944217
	3	1.18865714	1.188657	−1.11302087	−1.113159
	5	1.36433170	1.364332	−1.13918925	−1.139268
	10	1.52194504	1.521945	−1.15851866	−1.158571
s	2.5	1.00391561	1.003915	−1.48397701	−1.483987
	3	1.31012112	1.310121	−1.90669518	−1.906695
	5	2.39570505	2.395705	−3.42764833	−3.427648
	10	4.94949146	4.949491	−7.04073641	−7.040736

. Note that our convergent results are sufficiently accurate with high decimal precision, especially for large values of M, $\gamma$, s. HAM remains independent of small/large physical parameters and provides a convenient way to control the convergence of homotopy series solutions even for large disturbances, which distinguishes it from all other analytic techniques. Secondly, the homotopy approximation quickly converges with the optimal convergence-control parameters chosen. As shown in Figure 4, the total average residual error decreases sharply for each case as the order of approximation increases. Notice that, using the optimal convergence-control parameters $c_0^{f^{*}}=-0.058$, $c_0^{{\theta }^{*}}=-0.458$, results in a convergence speed faster than that using non-optimal values, such as $c_0^f=-0.075$, $c_0^{\theta }=-0.504$ for $\gamma =1$. The proper convergence-control parameters guarantee the convergence of the homotopy series solution. Furthermore, it is worth noting that other analytic techniques cannot guarantee the same. Thus, this demonstrates an obvious advantage of using HAM. In practice, the sufficiently accurate approximations are obtained in far fewer terms by the optimal convergence-control parameters.

5. Discussions

The explicit solutions of the velocity and temperature distribution (46)–(54) are obtained in Section 3, and are valid in the whole domain $\eta \ge 0$. Therefore, according to the explicit solutions, one can investigate the influence of physical quantities, such as Casson fluid parameter $\gamma$, magnetic interaction parameter M, suction parameter s and Prandtl number $Pr$ on the velocity and temperature profiles, skin friction coefficient $C_f$ and local Nusselt number $N{u}_x$.

5.1. Effect of $\gamma$

First, we investigate the effects of Casson fluid parameter $\gamma$ on the velocity profile and temperature profile. As shown in Figure 5, the magnitude of $f^{\prime }\left(\eta \right)$ decreases as $\gamma$ increases. The thickness of the velocity boundary layer decreases with increases in $\gamma$. This is because the yield stress decreases as $\gamma$ increases, which results in the velocity being suppressed. It is observed that when $\gamma$ approaches infinity, the problem in the given case reduces to a Newtonian case. We emphasize that the effect of $\gamma$ on the velocity profile in this study is consistent with the Casson fluid flow along an exponentially stretching surface [34]. Whether in regard to a stretching or shrinking surface, the magnitude of velocity is found to decrease with increasing $\gamma$. The decreasing nature of the momentum boundary layer thickness with increasing $\gamma$ appears accordingly. This relationship is reasonable in non-Newtonian fluids because an increase in yield stress suppresses the velocity in the boundary layer.

Figure 6 shows the effect of $\gamma$ on the temperature profile $\theta \left(\eta \right)$ in the cases of $M=2$, $s=3$ and $Pr=0.2$. It can be seen that the temperature decreases slightly with the increasing values of $\gamma$. Hence, the thermal boundary layer thickness decreases as the $\gamma$ increases. Here, we would like to illustrate a point regarding the wall temperature $Tw$. Note that $Tw$ given by Equation (6) in this study is a uniformly distributed constant. The temperature field given by constant $Tw$ is much more suppressed in the same value of $\gamma$ than that of the wall temperature condition increasing exponentially with x given by $Tw=T_{\infty }+T_0{e}^{\frac{x}{2L}}$. Figure 7 exhibits the temperature profiles for constant $Tw$ and exponentially increased $Tw$ in the cases of $\gamma =0.3$, $Pr=0.2$, $M=2$ and $s=3$, respectively. It is observed that the temperature increases for $Tw=T_{\infty }+T_0{e}^{\frac{x}{2L}}$ under the same parameters. Physically, as the wall temperature increases, the temperature of flow within the boundary layer increases. This causes an increase in the thermal boundary layer thickness.

5.2. Effects of M and s

The curves of $f^{\prime }\left(\eta \right)$ versus M and s are shown in Figure 8 and Figure 9, respectively. It is noticed that the magnitude of velocity profiles shows an appreciable decrease for large values of M and s. From Figure 8, the magnitude of velocity in the boundary layer is suppressed as M increases because the force of the magnetic field opposes the motion of the fluid. As shown in Figure 9, the magnitude of velocity decreases significantly with increasing mass suction, which causes a decrease in the boundary layer thickness. This phenomenon can be explained physically. The heated fluid is sucked closer to the wall as the mass suction becomes stronger, where the flow is slowed down due to the greater influence of viscosity. This effect suppresses the maximum velocity in the boundary layer. Therefore, an increased s leads to a faster reduction in the magnitude of velocity.

Figure 10 and Figure 11 show the effects of M and s on the temperature profile $\theta \left(\eta \right)$. Figure 10 depicts the features of temperature profile as a function of $\eta$ for various M. As M increases, the thermal temperature thickness becomes slightly thinner, and is not so sensitive to M. The convergent analytic approximation $\theta \left(\eta \right)$ for $s=2.5$, 3, 5 and 10 is shown in Figure 11. The increase in s obviously reduces the temperature profile. This change is quite significant due to the increase of mass suction. Compared with M, s has a greater impact on thermal distribution. Meanwhile, $\theta \left(\eta \right)$ obtained by HAM matches the numerical ones quite well for each case, as shown in Figure 11.

M and s show similar effects on temperature and velocity. Due to the applied magnetic field and suction, the velocity and temperature distributions become more uniform within the boundary layer. The presence of a magnetic field force opposite to the velocity direction and suction tends to reduce the momentum and thermal thickness of the boundary layer. This shows the effect of decreasing both the velocity and temperature within the boundary layer.

Additionally, similar to the analysis in $\gamma$, we compare the temperature profile for different wall temperature conditions in the cases of $M=2$ and $s=3$ (see Figure 12). It is observed that the temperature profile caused by $Tw=T_{\infty }+T_0$ decreases faster than that caused by $Tw=T_{\infty }+T_0{e}^{\frac{x}{2L}}$ under the same parameters. The exponentially increasing wall temperature with x raises the temperature of the fluid within the boundary layer and increases the thermal boundary layer thickness.

5.3. Effect of $Pr$

Notice that since the Prandtl number $Pr$ is closely related to the temperature of the boundary layer, we study the influences of $Pr$ on the temperature profile $\theta \left(\eta \right)$. As shown in Figure 13, the temperature profile decreases with the increasing values of $Pr$. The thermal boundary layer thickness is reduced with an increase in $Pr$. Since $Pr$ signifies the ratio of momentum diffusivity to thermal diffusivity, the thickness of momentum and thermal boundary layers can be controlled using $Pr$. The heat diffuses faster than the momentum for a small value of $Pr$. Therefore, the thermal boundary layer is thicker than that of the momentum boundary layer. In other words, higher thermal conductivity corresponds to a thicker thermal boundary layer. Therefore, it is observed from Figure 13 that the heat diffuses slowly and the thermal boundary layer becomes thinner as $Pr$ increases.

5.4. Analysis of $C_f$ and $N{u}_x$

The effects of s and $\gamma$ on skin friction $C_f$ and Nusselt number $N{u}_x$ are investigated by the explicit solutions of $f^{″}\left(0\right)$ and ${\theta }^{\prime }\left(0\right)$ in Section 4. Figure 2 illustrates that an increase in s leads to a linear increase in $f^{″}\left(0\right)$ and linear decrease in ${\theta }^{\prime }\left(0\right)$. According to Equations (11) and (12), $C_f$ and $N{u}_x$ exhibit a linearly increasing trend as s increases. However, $f^{″}\left(0\right)$ and ${\theta }^{\prime }\left(0\right)$ are nonlinear functions of $\gamma$ (see Figure 3). Note that $C_f\propto (1+1/\gamma )f^{″}\left(0\right)$, $C_f$ and $N{u}_x$ increase slowly with an increase in $\gamma$. As $\gamma$ increases to infinity, a Newtonian case appears.

In addition, the comparison of the present results corresponding to $f^{″}\left(0\right)$ and ${\theta }^{\prime }\left(0\right)$ with numerical solutions for various values of M, $\gamma$ and s is presented in Table 1. It shows the high accuracy of the applied scheme and verifies the effectiveness of the HAM approach. From Table 1, it is seen that the influence of M on $f^{″}\left(0\right)$ and ${\theta }^{\prime }\left(0\right)$ is similar to $\gamma$. An increment in M causes $f^{″}\left(0\right)$ and $-{\theta }^{\prime }\left(0\right)$ to increase slowly. The same is true for the trend of $C_f$ and $N{u}_x$ due to Equations (11) and (12).

Figure 14 displays the wall temperature gradient $-{\theta }^{\prime }\left(0\right)$ against $Pr$. It is clearly shown that the increment in $Pr$ leads to a linear increase in the wall temperature gradient. The wall temperature gradient is proportional to the heat transfer rate or $N{u}_x$. As mentioned earlier, $Pr$ represents the ratio of momentum diffusivity to thermal diffusivity. The momentum diffusivity increases, whereas thermal diffusivity decreases as $Pr$ increases, so the heat transfer rate increases. This reflects that the Nusselt number $N{u}_x$ increases with the increasing values of $Pr$. Significantly, the heat transfer rate under the condition of an exponentially increasing wall temperature is lower than that under a constant wall temperature. As can be observed from

Table 2. Analytic approximations (HAM) of $-{\theta }^{\prime }\left(0\right)$ under different wall temperature conditions for various $Pr$ with $M=\gamma =s=3$.

Parameters		$-{\mathit{\theta }}^{\prime }\left(0\right)$
		$\mathit{T}\mathit{w}={\mathit{T}}_{\mathbf{\infty }}+{\mathit{T}}_{\mathbf{0}}$	$\mathit{T}\mathit{w}={\mathit{T}}_{\mathbf{\infty }}+{\mathit{T}}_{\mathbf{0}}{\mathit{e}}^{\frac{\mathit{x}}{\mathbf{2}\mathit{L}}}$
$Pr$	0.1	0.2664	0.23258
	0.2	0.53845	0.47562
	0.5	1.37692	1.24919
	1	2.81671	2.62471
	4	11.71837	11.42673

, the convergent solution $-{\theta }^{\prime }\left(0\right)$ is reduced by increasing the wall temperature along the x-direction. This means that the temperature profile within the boundary layer enhances, which causes an increase in thermal boundary layer thickness.

Finally, it should be emphasized that our series solutions are sufficiently accurate by comparison with the numerical results. In addition, as observed in Figure 4, the squared residual error of the homotopy approximations decreases exponentially as the order of approximation increases. Thus, it is in fact unnecessary to compare our convergent results with numerical ones. The variation of the squared residual error evaluates the accuracy of the homotopy approximation.

6. Conclusions

The MHD flow and heat transfer of a Casson fluid over an exponentially shrinking sheet with suction is investigated using the HAM approach. First, the governing boundary layer equations are transformed into nonlinear ordinary differential equations using similarity transformations. Then, the nonlinear ordinary differential governing Equations (8) and (9) are replaced utilizing an infinite number of linear sub-equations, which are solved analytically in the whole domain. Due to the freedom in constructing the zeroth-order deformation equations, we can choose the appropriate initial guesses and the auxiliary linear operators so that the explicit solutions are derived further. The generosity in freedom is based on the concept of convergence-control by means of convergence-control parameters. The optimal values of the convergence-control parameters are strongly suggested for use in practice. In general, it is enough to obtain an accurate homotopy approximation by using optimal convergence-control parameters determined by the minimum of the squared residual error corresponding to the governing equations.

In this study, by solving the coupled nonlinear differential equations in the MHD flow and heat transfer of a Casson fluid over an exponentially shrinking sheet, we arrive at the following main points:

• The explicit analytic solutions of $f\left(\eta \right)$ and $\theta \left(\eta \right)$ are obtained and valid in the whole region $\eta =[0,+\infty )$.
• The important quantities $f^{″}\left(0\right)$ and ${\theta }^{\prime }\left(0\right)$ related to the skin friction coefficient $C_f$ and local Nusselt number $N{u}_x$ are derived in an explicit form.
• The convergent analytic solutions are in good agreement with the numerical solutions. The rapid decrease in squared residual error ensures the accuracy of the homotopy approximation.
• An increase in the Casson fluid parameter $\gamma$ suppresses the magnitude of velocity profile $f^{\prime }\left(\eta \right)$ due to the reduced yield stress as $\gamma$ increases. This leads to a thinner momentum boundary layer thickness. The velocity profile magnitude $f^{\prime }\left(\eta \right)$ is found to decrease with increasing $\gamma$ for both stretching and shrinking surfaces.
• The temperature profile $\theta \left(\eta \right)$ decreases slightly with increasing values of $\gamma$ in the current case, which decreases the thermal boundary layer thickness.
• The magnitudes of $f^{\prime }\left(\eta \right)$ and $\theta \left(\eta \right)$ decrease significantly with increases in the magnetic interaction parameter M and suction parameter s.
• The velocity and thermal boundary layer thicknesses decrease as M and s increase. The presence of a magnetic field force opposite to the velocity and suction reduces the momentum and thermal thickness of the boundary layer.
• The temperature profile and thermal boundary layer thickness decrease with increasing values of Prandtl number $Pr$. The heat diffuses faster corresponding to the higher thermal conductivity for a small value of $Pr$.
• $C_f$ and $N{u}_x$ exhibit a linearly increasing trend as s becomes stronger.
• $C_f$ and $N{u}_x$ increase nonlinearly with increases in $\gamma$ and M.
• The wall heat transfer rate $-{\theta }^{\prime }\left(0\right)$ increases linearly as $Pr$ increases, as a result, $N{u}_x$ also increases linearly.
• Compared with the constant wall temperature condition, the exponentially increasing wall temperature with x raises the temperature of the fluid within the boundary layer and leads to increased thickness of the thermal boundary layer.