On the Numerical Analysis of Unsteady MHD Boundary Layer Flow of Williamson Fluid Over a Stretching Sheet and Heat and Mass Transfers

Abstract

A thorough and detailed investigation of an unsteady free convection boundary layer flow of an incompressible electrically conducting Williamson fluid over a stretching sheet saturated with a porous medium has been numerically carried out. The partial governing equations are transferred into a system of non-linear dimensionless ordinary differential equations by employing suitable similarity transformations. The resultant equations are then numerically solved using the spectral quasi-linearization method. Numerical solutions are obtained in terms of the velocity, temperature and concentration profiles, as well as the skin friction, heat and mass transfers. These numerical results are presented graphically and in tabular forms. From the results, it is found out that the Weissenberg number, local electric parameter, the unsteadiness parameter, the magnetic, porosity and the buoyancy parameters have significant effects on the flow properties.

1. Introduction

The study of non-Newtonian fluids has attracted many researchers owing to its enormous applications in industrial and engineering sectors, especially in manufacturing and processing industries. Unlike Newtonian fluids, Kahshan et al. [1], non-Newtonian fluids are more complicated because there is no single construction relation that can be used to explain them all, Hussanan et al. [2]. The relationship between the shear stress and rate of strain is non-linear at a given temperature and pressure in a non-Newtonian fluid. To this end, these fluids cannot be modelled using the classical Navier-Stokes equations. Non-Newtonian fluids have the ability to: (i) shear-thin or shear thicken, (ii) exhibit thixotropy, (iii) allow stress relaxation, (iv) creep in a nonlinear manner, (v) develop normal stress differences, and (vi) exhibit a threshold for the shear stress before it starts to flow. Many non-Newtonian fluids possess one or more of these characteristics.

Non-Newtonian fluids are generally classified into three main categories, namely: (i) the differential type, (ii) the rate type and (iii) the integral type. The detailed descriptions of each category can be found in Cioranescu et al. [3]. In non-Newtonian fluids, pseudoplastic fluids are the most frequently encountered fluids. But as expected, Navier-stokes equations alone are insufficient to describe the rheological properties of these fluids. Therefore, to overcome this challenge, many rheological models like Maxwell model, Jeffrey model, Ellis model, Power law model, Carreau model, among others, have been developed. The Williamson fluid model is another powerful model used to explain the rheological properties of pseudoplastic fluids. A pseudoplastic is a shear-thinning fluid and has less resistance at high strain rates like polymer solution, paint, blood, plasma. In rheology, shear thinning is the pseudoplastic fluid whose viscosity decreases under shear strain. The shear thinning fluid behaviour is usually exhibited by inks which are used in inkjet printing as discussed by Miccichè et al. [4] and Dybowska-Sarapuk et al. [5].

A lot of authors have investigated the flow of non-Newtonian fluids due to their vast applications especially when with suspensions of nano-sized particle. Some work on nanoparticles was done by Mozaffari et al. [6], Mozaffari et al. [7], Darjani et al. [8] and Xing et al. [9]. Considering a flow underlying spreading surface through a non-Darcian porous medium, Elgazery [10] anlayzed the effects of internal heat generation/absorption of a non-Newtonian Casson fluid with suspension of gold and alumina nanoparticles. Hsiao [11] presented a study for thermal energy extrusion system conversion problem with electric hydromagnetic heat and mass mixed convection of a viscoelastic non-Newtonian Carreau-Nanofluid with radiation and viscous dissipation effects. Khan et al. [12] studied the influence of the nanoparticles and uniform magnetic field on the slip flows in arterial vessels with blood conveyed through hollow arterial tubes described as a third grade non-Newtonian fluid. An investigation on the multislip effects on the magnetohydrodynamic mixed convection unsteady flow of microploar nanofluids over a stretching/shrinking sheet along with radiation in the presence of a heat source was done by Abdal et al. [13]. Adesanya et al. [14] investigated the steady flow of the non-Newtonian fluid via the inclined channel heated isothermally at the boundaries.

Williamson [15] pioneered the discussion of the flow of pseudoplastic materials and proposed a model equation to describe the flow of pseudoplastic fluids and experimentally verified the results. Nadeem et al. [16] presented a paper modelling a two-dimensional Williamson fluid flow over a stretching sheet. Nadeem and Hussain [17] explored the effects of heat transfers on the Williamson fluid over a porous exponentially stretching sheet surface. The study also considered the two cases of heat transfer, namely, the prescribed exponential order surface temperature case and the prescribed exponential order heat flux case. Khan et al. [18] described the effect of thermal radiation on the thin film nanofluid flow of a Williamson fluid over an unsteady stretching surface with variable fluid properties. Ijaz Khan et al. [19] developed a model for a boundary layer stagnation point flow of an electrically conducting Williamson fluid in the presence of a constant magnetic field.

Monica et al. [20] presented an analysis that dealt with the study of stagnation point flow of a Williamson fluid over nonlinearly stretching sheet with thermal radiation. Malik et al. [21] studied the Williamson fluid past a stretching cylinder with combined effects of variable thermal conductivity and heat generation/absorption. Mabood et al. [22] performed an analysis for MHD Williamson nanofluid flow over a continuously moving heated surface with thermal radiation and heat source. Dawar et al. [23] did an analysis of the flow of a Williamson fluid taken over a linear porous stretching sheet under the influence of thermal radiation.

Kurmar et al. [24] carried out a mathematical analysis of two-phase boundary layer flow and heat transfer of a Williamson fluid with particle suspension over a stretching sheet. Shateyi et al. [25] investigated a Casson fluid flow in the presence of free convection of combined heat and mass transfer towards an unsteady permeable stretching sheet with thermal radiation, viscous dissipation and chemical reaction. Hayat et al. [26] studied the unsteady two-dimensional boundary layer flow of an incompressible Williamson fluid over an unsteady permeable stretching surface with thermal radiation. Khan et al. [27] examined the influence of chemically reactive species and mixed convection on the magnetohydrodynamic Williamson nanofluid induced by a nonisothermal cone and plate in a porus medium.

Recently, Panezai et al. [28] examined the influence of thermal radiation on two-dimensional incompressible MHD mixed convective heat transfer flow of Williamson fluid flowing over a porous wedge. Hsu et al. [29] numerically analysed the heat and mass transfer characteristics of the influence of uniform blowing/suction and MHD (magnetohydrodynamic) on the free convection of non-newtonian fluids over a vertical plate in porous media with internal heat generation and soret/dufour effects. Kebede [30] presented an analytic approximation to the heat and mass transfer characteristics of Williamson nanofluid flow. Reddy et al. [31] studied MHD flow and heat transfer characteristics of Williamson nanofluid over a stretching sheet with variable thickness and variable thermal conductivity.

Lastly, Megahed [32] researched on Williamson boundary layer fluid flow and heat transfer due to a nonlinearly stretching sheet. Many researchers have been attracted to study viscous flows due to a stretching sheet. This has been necessitated by their several applications in polymer processing industries, environmental pollution, biological processes and aerodynamic extrusion of plastic sheets, among other applications. Sakiadis [33] pioneered the discussion of the fluid flow due to a stretching surface. Crane [34] extended the work of Sakiadis to study the problem of fluid flow of Blasius type due to a stretching sheet. The literature on the stretching sheet topic is quite extensive and hence cannot be listed here in detail. The most recent works of great researchers regarding the flow over stretching sheet can be found in Hsiao [35,36], Shateyi [37], Sharma et al. [38], Shamashuddin et al. [39], Nagalakshm and Vijaya [40].

Motivated by the above mentioned studies as well as the vast applications of the different types of non-Newtonian fluids, the current study seeks to investigate unsteady free convection boundary layer flow of an electrically conducting Williamson fluid over a stretching sheet. The sheet is saturated with a porous medium under the combined effects of viscous dissipation, chemical reaction, thermal radiation as well as a uniform magnetic field. It is well known that fluid dynamics problems are analysed through material and virtual experimentation, Vedovoto et al. [41]. Material modeling entails the construction of benches and models of the physical problems, as well as instrumentation for gathering information. However, the virtual experimentation requires mathematical modeling which leads to models with differential equations, integral equations and/or intro-differential equations. It is remarked that valuable information about the problem is provided by virtual experimentation. Numeric algorithms and quantitative information about the problem are generated through computational modelling. The data is then used to perform the visualization of the simulation results ad the statistical analysis of the numerical experiment. To that end, the current study also seeks to use the virtual experimentation approach. We are also employing recently developed numerical technique known as the spectral quasi-linearization method. Therefore, novelty of the current paper lies in the application of the numerical technique to solve this two dimensional, unsteady free convection boundary layer flow of an incompressible electrically conducting Williamson fluid over a stretching sheet saturated with a porous medium. The innovation of this study lies in putting together more factors into the basic Williamson fluid model to give a through analysis of the model.

2. Mathematical Formulation

We consider a two dimensional, unsteady free convection boundary layer flow of an incompressible electrically conducting Williamson fluid over a stretching sheet saturated with a porous medium. The x-axis is taken parallel to the stretching sheet and the y-axis is in the vertical direction. We assume that the sheet is moving with a velocity $U_w=\frac{ax}{1-ct}$, where $a>0$ is a stretching rate along the x-axis. If $a<0$, then it becomes a shrinking velocity constraint. The the constant $ct$ is such that $ct<1$. A uniform transverse magnetic field $\overrightarrow{B}=(0,B_0,0)$ and uniform electric field $\overrightarrow{E}=(0,0,-E_0)$ are applied to the flow region as shown in Figure 1. The magnetic Reynolds number is assumed to be much less than unity and hence the induced magnetic field is negligible in comparison to the applied magnetic field. This is made possible by the fact that the fluid is assumed to be slightly conducting. In this study we also neglect the Hall effect/currents. We note that the magnetic field is weaker than the electric field and also that the magnetic field obeys Ohm’s law $\check{J}=\sigma (\overrightarrow{E}+\overline{V}\times \overrightarrow{B}),$ where $\check{J}$ is the Joule current, $\sigma$ is the electrical conductivity and $\overline{V}$ is the velocity. We take into account frictional heating in the from of viscous dissipation. Thermal radiation and chemical reaction are also considered in this study. We neglect the induced magnetic field and the hall effect.

The relevant governing flow equations through the conservation laws of mass, linear momentum, energy and mass by using the boundary layer approximations give the following Hayat et al. [26]

$$\begin{array}{c}\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0,\\ \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=\nu \frac{{\partial }^{2}u}{\partial y^2}+2\nu \mathsf{\Gamma }\frac{\partial u}{\partial y}\frac{{\partial }^{2}u}{\partial y^2}+g\beta (T-T_{\infty })+g{\beta }_c(C-C_{\infty })\end{array}$$

(1)

$$\begin{array}{c}+\frac{\sigma }{\rho }(E_0{B}_0-B_0^{2}u),\end{array}$$

(2)

$$\begin{array}{c}\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\frac{K}{\rho c_p}\frac{{\partial }^{2}T}{\partial y^2}+\frac{\mu }{\rho c_p}{\left(\frac{\partial u}{\partial y}\right)}^2+\frac{\mu \mathsf{\Gamma }}{\rho c_p}{\left(\frac{\partial u}{\partial y}\right)}^3+\frac{\sigma }{\rho c_p}{({\beta }_{0}u-E_0)}^2-\frac{1}{\rho c_p}\frac{\partial q_r}{\partial y},\end{array}$$

(3)

$$\begin{array}{c}\frac{\partial C}{\partial t}+u\frac{\partial C}{\partial x}+v\frac{\partial C}{\partial y}=D\frac{{\partial }^{2}C}{\partial y^2}-K_r(C-C_{\infty }),\end{array}$$

(4)

subject to boundary conditions

$$\begin{array}{c}\hfill \left\{\begin{array}{c}u(x,0)=U_w,v(x,0)=V_w,T(x,0)=T_w,C(x,0)=C_w,\hfill \\ u\to 0,T\to T_{\infty },C\to C_{\infty },\mathrm{as}y\to \infty .\hfill \end{array}\right.\end{array}$$

(5)

Here u and v are the velocity components in the x and y directions, respectively, $\nu$ is the kinematic viscosity, $\mathsf{\Gamma }$ is the relaxation time, g is the acceleration due to gravity, T and C are the fluid temperature and concentration, respectively, $T_{\infty }$ and $C_{\infty }$ are the respective ambient temperature and concentration, $\beta$ and ${\beta }_c$ are the respective coefficients of thermal and concentration expansion, $\sigma$ is the electrical conductivity, $\rho$ is the fluid density, $\mu$ is the dynamic viscosity, $\kappa$ is the permeability parameter, $c_p$ is the specific heat at constant pressure, $q_r$ is the radiative heat flux, $K_r$ is the chemical reaction rate, K is thermal conductivity and lastly, D is the mass diffusivity. $V_w$ is expressed by $V_w=-\frac{V_0}{{(1-ct)}^{\frac{1}{2}}}$.

As now known, $V_w$ express the mass transfer at the surface with $V_w<0$ for injection and $V_w>0$ being for suction. The surface temperature $T_w(x,t)$ and the surface concentration $C_w(x,t)$ are given respectively as:

$$\begin{array}{c}\hfill T_w(x,t)=T_{\infty }+\frac{ax}{2\nu {(1-ct)}^2}{T}_0\mathrm{and}{C}_w(x,t)=C_{\infty }+\frac{ax}{2\nu {(1-ct)}^2}{C}_0,\end{array}$$

(6)

where $c\ge 0$. $T_0(0\le T_0\le T_w)$ and $C_0(0\le C_0\le C_w)$ are the reference temperature and concentration respectively.

By employing the Rosseland approximations, we have

$$q_r=-\frac{4{\sigma }^{*}}{3k_1}\frac{\partial T^4}{\partial y},$$

(7)

where ${\sigma }^{*}$ is the Stefan-Boltzmann constant and $k_1$ is the mean absorption coefficient. Applying Taylor’s series, we have $T^4\simeq 4T_{\infty }^{3}T-3T_{\infty }^4,$ where $T_{\infty },$ is the ambient temperature and the energy equation now becomes:

$$\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\frac{1}{\rho c_p}\left(\frac{16{\sigma }^{*}{T}_{\infty }^3}{3k_1}+K\right)\frac{{\partial }^{2}T}{\partial y^2}+\frac{\mu }{\rho c_p}{\left(\frac{\partial u}{\partial y}\right)}^2+\frac{\mu \mathsf{\Gamma }}{\rho c_p}{\left(\frac{\partial u}{\partial y}\right)}^3+\frac{\sigma }{\rho c_p}{({\beta }_{0}u-E_0)}^2.$$

(8)

Similarity Transformation

Following Hayat et al. [26], among others, we introduce the following dimensionless variables.

$$\begin{array}{c}\hfill \eta =\sqrt{\frac{U_w}{x\nu }}y,\psi =\sqrt{\nu x{U}_w}f\left(\eta \right),\theta \left(\eta \right)=\frac{T-T_{\infty }}{T_w-T_{\infty }},\varphi \left(\eta \right)=\frac{C-C_{\infty }}{C_w-C_{\infty }},\end{array}$$

(9)

with the velocity components given by $u=\frac{\partial \psi }{\partial y}$ and $v=-\frac{\partial \psi }{\partial x}$.

The continuity equation is satisfied and the resulting governing partial differential equations together with the boundary conditions become:

$$\begin{array}{c}(1+2We{f}^{″})f^{‴}+f{f}^{″}-{\left(f^{\prime }\right)}^2-S\left(f^{\prime }+\frac{\eta }{2}{f}^{″}\right)+M^2(E_1-f^{\prime })+{\lambda }_1\theta +{\lambda }_2\varphi =0,\end{array}$$

(10)

$$\begin{array}{c}\frac{1}{Pr}\left(1+\frac{4R}{3}\right){\theta }^{″}+f{\theta }^{\prime }-\frac{S\eta }{2}{\theta }^{\prime }-(2S+f^{\prime })\theta +Ec({\left(f^{″}\right)}^2+We{\left(f^{″}\right)}^3)+M^{2}Ec{(f^{\prime }-E_1)}^2=0,\end{array}$$

(11)

$$\begin{array}{c}\frac{1}{Sc}{\varphi }^{″}+\left(f-\frac{S\eta }{2}\right){\varphi }^{\prime }+(2S-f^{\prime }-\gamma )\varphi =0.\end{array}$$

(12)

The corresponding boundary conditions become:

$$\begin{array}{c}f\left(0\right)=f_w,f^{\prime }\left(0\right)=1,\theta \left(0\right)=1,\varphi \left(0\right)=1,\hfill \end{array}$$

(13)

$$\begin{array}{c}{f}^{\prime }(\infty )\to 0,f^{″}(\infty )\to 0\theta (\infty )\to 0,\varphi (\infty )\to 0,\hfill \end{array}$$

(14)

where $S\left(=\frac{c}{a}\right)$ is the unsteadiness parameter, $f_w\left(=\frac{V_0}{\sqrt{a\nu }}\right)$ is the suction/injection parameter, $We\left(=\mathsf{\Gamma }{U}_w\sqrt{\frac{a}{\nu (1-ct)}}\right)$ is the Weissenberg number, $M^2\left(=\frac{\sigma B_0^2(1-ct)}{\rho a}\right)$ and M is the magnetic field parameter, $E_1\left(=\frac{E_0(1-ct)}{B_{0}ax}\right)$ is the local electric parameter, $R\left(=\frac{4{\sigma }^{*}{T}_{\infty }^3}{\kappa K}\right)$ is the thermal radiation parameter, $Sc\left(=\frac{\nu }{D}\right)$ is the Schmidt number, $Pr\left(=\frac{\mu c_p}{K}\right)$ is the Prandtl number, $\gamma \left(=\frac{K_r(1-ct)}{a}\right)$ is the chemical reaction parameter, ${\lambda }_1\left(=\frac{g\beta T_0}{2a\nu }\right)$ is the thermal buoyancy parameter, ${\lambda }_2\left(=\frac{g{\beta }_c{C}_0}{2a\nu }\right)$ is the solutal buoyancy parameter, $Ec\left(=\frac{U_w^2}{c_p(T_w-T_0)}\right)$ is the Eckert number.

The flow characteristics which are of engineering significance are the skin friction coefficient, the local Nusselt number and the Sherwood number. These are respectively, defined as: (Hayat et al. [26])

$$\begin{array}{c}\hfill C_f=\frac{{\tau }_w}{\rho U_w^2},N{u}_x=\frac{x{q}_w}{K(T_w-T_{\infty })},S{h}_x=\frac{x{j}_w}{k_{\infty }(C_w-C_{\infty })}.\end{array}$$

(15)

Upon applying the necessary expressions for ${\tau }_w,q_w,\mathrm{and}{j}_w$, we get the following:

$$\begin{array}{c}\hfill C_f=R{e}_x^{-1/2}[f^{″}\left(0\right)+We{f^{″}}^2\left(0\right)],N{u}_x=-R{e}_x^{1/2}\left(1+\frac{4R}{3}\right){\theta }^{\prime }\left(0\right),S{h}_x=-R{e}_x^{1/2}{\varphi }^{\prime }\left(0\right).\end{array}$$

(16)

3. Method of Solution

In this section we present a spectral method based on a quasi-linearization method called the spectral quasi-linearization method (SQLM), Motsa et al. [42]. The SQLM combines two methods, the quasi-linearization method (QLM) and a Chebyshev spectral collocation method (CSCM). The QLM, a Newton-Raphson based quasi-linearzation method which was introduced by Bellman and Kalaba [43], is used to linearize the governing non-linear equations. The CSCM is used to integrate the resulting iterative sequence of linear differential equations.

3.1. Quasi-Linearization

With reference to Equations (10)–(12), we derive the quasi-linearization formula by considering a system of differential equations

$$\begin{array}{c}\hfill {\mathsf{\Lambda }}_p[{\zeta }^1,{\zeta }^2,{\zeta }^3]=0,p=1,2,3,\end{array}$$

(17)

where ${\mathsf{\Lambda }}_p(p=1,2,3)$ are non-linear operators, ${\zeta }^1=\{f,f^{\prime },f^{″},f^{‴}\}$, ${\zeta }^2=\{\theta ,{\theta }^{\prime },{\theta }^{″}\}$ and ${\zeta }^3=\{\varphi ,{\varphi }^{\prime },{\varphi }^{″}\}$. Denote $f_s,{\theta }_s,{\varphi }_s$ and $f_{s+1},{\theta }_{s+1},{\varphi }_{s+1}$ be the solutions of Equation (17) at the previous and the current iteration levels, respectively. Assuming the difference between the solution at the previous and the current iteration levels and its derivatives is small enough, linear Taylor series expansion of Equation (17) about some previous solution, upon simplification yields

$$\begin{array}{c}{\zeta }_{s+1}^1·{\nabla }_f{\mathsf{\Lambda }}_p[f_s,{\theta }_s,{\varphi }_s]+{\zeta }_{s+1}^2·{\nabla }_{\theta }{\mathsf{\Lambda }}_p[f_s,{\theta }_s,{\varphi }_s]+{\zeta }_{s+1}^3·{\nabla }_{\varphi }{\mathsf{\Lambda }}_p[f_s,{\theta }_s,{\varphi }_s]=\hfill \\ {\zeta }_s^1·{\nabla }_f{\mathsf{\Lambda }}_p[f_s,{\theta }_s,{\varphi }_s]+{\zeta }_s^2·{\nabla }_{\theta }{\mathsf{\Lambda }}_p[f_s,{\theta }_s,{\varphi }_s]+{\zeta }_s^3·{\nabla }_{\varphi }{\mathsf{\Lambda }}_p[f_s,{\theta }_s,{\varphi }_s]\hfill \\ -{\mathsf{\Lambda }}_p[f_s,{\theta }_s,{\varphi }_s],s=0,1,2,\cdots ,\hfill \end{array}$$

(18)

where ${\nabla }_f=\{\frac{\partial }{\partial f},\frac{\partial }{\partial f^{\prime }},\frac{\partial }{\partial f^{″}},\frac{\partial }{\partial f^{‴}}\}$, ${\nabla }_{\theta }=\{\frac{\partial }{\partial \theta },\frac{\partial }{\partial {\theta }^{\prime }},\frac{\partial }{\partial {\theta }^{″}}\}$ and ${\nabla }_{\varphi }=\{\frac{\partial }{\partial \varphi },\frac{\partial }{\partial {\varphi }^{\prime }},\frac{\partial }{\partial {\varphi }^{″}}\}$. Applying the formula Equation (18) to the system Equations (10)–(12) results in an iterative sequence of linear equations

$$\begin{array}{ccc}\hfill a_{0,s}{f}_{s+1}^{‴}+a_{1,s}{f}_{s+1}^{″}+a_{2,s}{f}_{s+1}^{\prime }+a_{3,s}{f}_{s+1}+a_{4,s}{\theta }_{s+1}+a_{5,s}{\varphi }_{s+1}& =& R_{1,s},\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill b_{0,s}{\theta }_{s+1}^{″}+b_{1,s}{\theta }_{s+1}^{\prime }+b_{2,s}{\theta }_{s+1}+b_{3,s}{f}_{s+1}^{″}+b_{4,s}{f}_{s+1}^{\prime }+b_{5,s}{f}_{s+1}& =& R_{2,s},\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill c_{0,s}{\varphi }_{s+1}^{″}+c_{1,s}{\varphi }_{s+1}^{\prime }+c_{2,s}{\varphi }_{s+1}+c_{3,s}{f}_{s+1}^{\prime }+c_{4,s}{f}_{s+1}& =& R_{3,s},\hfill \end{array}$$

(21)

where the variable coefficients are given by

$$\begin{array}{ccc}\hfill a_{0,s}& =& 1+2We{f}_s^{″},a_{1,s}=2We{f}_s^{‴}+f_s-\frac{\eta }{2}S,a_{2,s}=-2f_s^{\prime }-S-M^2,a_{3,s}=f_s^{″},a_{4,s}={\lambda }_1,\hfill \\ \hfill a_{5,s}& =& {\lambda }_2,b_{0,s}=\frac{1}{Pr}(1+\frac{4}{3}R),b_{1,s}=f_s-\frac{\eta }{2}S,b_{2,s}=-2S-f_s^{\prime },b_{3,s}=2E_c{f}_s^{″}+3E_{c}We{\left(f_s^{″}\right)}^2,\hfill \\ \hfill b_{4,s}& =& 2M^2{E}_c(f_s^{\prime }-E_1)-{\theta }_s,b_{5,s}={\theta }_s^{\prime },c_{0,s}=\frac{1}{Sc},c_{1,s}=f_s-\frac{\eta }{2}S,c_{2,s}=2S-f_s^{\prime }-\gamma ,\hfill \\ \hfill c_{3,s}& =& -{\varphi }_s,c_{4,s}={\varphi }_s^{\prime }.\hfill \end{array}$$

The terms on the right are defined as follows:

$$\begin{array}{ccc}\hfill R_{1,s}& =& 2We{f}_s^{″}{f}_s^{‴}+f_s{f}_s^{″}-{\left(f_s^{\prime }\right)}^2-M^2{E}_1,\hfill \\ \hfill R_{2,s}& =& f_s{\theta }_s^{\prime }-f_s^{\prime }{\theta }_s+E_c{\left(f_s^{″}\right)}^2+2E_{c}We{\left(f_s^{″}\right)}^3+M^2{E}_c[{\left(f_s^{\prime }\right)}^2-E_1^2],\hfill \\ \hfill R_{3,s}& =& f_s{\varphi }_s^{\prime }-{\varphi }_s{f}_s^{\prime }.\hfill \end{array}$$

3.2. Chebyshev Differentiation

We approximate the solutions of Equations (10)–(12) by functions of the form

$$\begin{array}{ccc}\hfill f\left(\eta \right)& =& \sum _{k=0}^{N}f\left({\eta }_k\right)L_k\left(\eta \right),\theta \left(\eta \right)=\sum _{k=0}^N\theta \left({\eta }_k\right)L_k\left(\eta \right),\hfill \\ \hfill \varphi \left(\eta \right)& =& \sum _{k=0}^N\varphi \left({\eta }_k\right)L_k\left(\eta \right),\hfill \end{array}$$

(22)

where

$$L_k\left(\eta \right)=\prod _{\begin{array}{c}j=0\\ k\ne j\end{array}}^N\frac{\eta -{\eta }_j}{{\eta }_k-{\eta }_j},$$

is the Lagrange interpolating polynomial. For convenience in the application of the collocation method, we transform the semi-infinite physical domain $[0,\infty )$ in the $\eta$-direction to $[-1,1]$ in the x-direction using the linear transformation $\eta =\frac{l_{\infty }}{2}(1+x)$. ($l_{\infty }$ is a number chosen to be large enough that boundary conditions at ∞ hold). Evaluating the approximating functions (22) and their derivatives at Chebyshev-Gauss-Lobatto points

$$\begin{array}{c}\hfill \left\{{\eta }_i\right\}=cos\left(\frac{\pi i}{N}\right),i=0,1,\cdots ,N,\end{array}$$

gives

$$\begin{array}{ccc}\hfill \frac{df}{d\eta }{|}_{\eta ={\eta }_i}& =& \sum _{k=0}^{N}f\left({\eta }_k\right)\frac{d}{d\eta }{L}_k\left({\eta }_i\right)=\mathbf{D}\mathbf{F},\hfill \\ \hfill \Rightarrow \frac{d^{n}f}{d{\eta }^n}{|}_{\eta ={\eta }_i}& =& {\mathbf{D}}^n\mathbf{F},\hfill \end{array}$$

(23)

where $\mathbf{D}=\left(\frac{2}{l_{\infty }}\right)\mathcal{D}$, $\mathcal{D}$ is an $(N+1)\times (N+1)$ Chebyshev differentiation matrix, Trefethen [44], $\mathbf{F}={[f\left({\eta }_0\right),f\left({\eta }_1\right),\cdots ,f\left({\eta }_{N-1}\right),f\left({\eta }_N\right)]}^T$ and T denotes matrix transpose. $\theta$ and $\varphi$ have similar expressions for derivatives. Evaluating Equations (19)–(21) at the collocation points and approximating derivatives with Chebyshev derivatives gives, in vector matrix form, the system

$$\begin{array}{c}\hfill \left[\begin{array}{ccc}{\mathbf{A}}_{11}& {\mathbf{A}}_{12}& {\mathbf{A}}_{31}\\ {\mathbf{A}}_{21}& {\mathbf{A}}_{22}& {\mathbf{A}}_{23}\\ {\mathbf{A}}_{31}& {\mathbf{A}}_{32}& {\mathbf{A}}_{33}\end{array}\right]\left[\begin{array}{c}{\mathbf{F}}_{s+1}\\ {\mathsf{\Phi }}_{s+1}\\ {\mathsf{\Theta }}_{s+1}\end{array}\right]=\left[\begin{array}{c}{\mathbf{R}}_{1,s}\\ {\mathbf{R}}_{2,s}\\ {\mathbf{R}}_{3,s}\end{array}\right],\end{array}$$

(24)

where

$$\begin{array}{c}\hfill {\mathbf{A}}_{11}={\mathbf{a}}_{0,s}{\mathbf{D}}^3+{\mathbf{a}}_{1,s}{\mathbf{D}}^2+{\mathbf{a}}_{2,s}\mathbf{D}+{\mathbf{a}}_{3,s},{\mathbf{A}}_{12}={\mathbf{a}}_{4,s},{\mathbf{A}}_{13}={\mathbf{a}}_{5,s}\\ \hfill {\mathbf{A}}_{21}={\mathbf{b}}_{3,s}{\mathbf{D}}^2+{\mathbf{b}}_{4,s}\mathbf{D}+{\mathbf{b}}_{5,s},{\mathbf{A}}_{22}={\mathbf{b}}_{0,s}{\mathbf{D}}^2+{\mathbf{b}}_{1,s}\mathbf{D}+{\mathbf{b}}_{2,s}.{\mathbf{A}}_{23}=\mathbf{0},\\ \hfill {\mathbf{A}}_{31}={\mathbf{c}}_{3,s}\mathbf{D}+{\mathbf{c}}_{4,s},{\mathbf{A}}_{32}=\mathbf{0},{\mathbf{A}}_{33}={\mathbf{c}}_{0,s}{\mathbf{D}}^2+{\mathbf{c}}_{1,s}\mathbf{D}+{\mathbf{c}}_{2,s},\\ \hfill {\mathbf{R}}_{k,s}={\{R_{k,s}\left({\eta }_0\right),R_{k,s}\left({\eta }_1\right),\cdots ,R_{k,s}\left({\eta }_{N-1}\right),R_{k,s}\left({\eta }_N\right)\}}^T,(k=1,2,3).\end{array}$$

4. Results and Discussions

The spectral quasi-linearization method was used to generate the results to be discussed in this section. All the numerical results were obtained using MATLAB 2016. For all the computations, the default parameters considered, unless otherwise stated, are $N=60$, $Pr=1.0$, $We=0.1$, $E_1=0.3$, $E_c=0.5$, $R=0.3$, ${\lambda }_1=0.1$, ${\lambda }_2=0.1$, $Sc=0.1$, $\gamma =1.0$ and $f_w=0.3$.

The error infinity norms defined by

$$\begin{array}{ccc}\hfill \mathrm{Error}\left(\mathrm{f}\right)& =& \underset{0\le i\le N}{max}\left|\right|{\mathbf{F}}_{s+1,i}-{\mathbf{F}}_{s,i}|{|}_{\infty },\mathrm{Error}\left(\theta \right)=\underset{0\le i\le N}{max}||{\mathsf{\Theta }}_{s+1,i}-{\mathsf{\Theta }}_{s,i}|{|}_{\infty },\hfill \\ \hfill \mathrm{Error}\left(\varphi \right)& =& \underset{0\le i\le N}{max}\left|\right|{\mathsf{\Phi }}_{s+1,i}-{\mathsf{\Phi }}_{s,i}|{|}_{\infty },\hfill \end{array}$$

are used to confirm the convergence of the SQLM solution. The accuracy of the SQLM is tested using the residual error infinity norms defined as follows:

$$\begin{array}{ccc}\hfill {\left|\right|Res\left(f\right)\left|\right|}_{\infty }& =& \left|\right|(1+2We{f}^{″})f^{‴}+f{f}^{″}-{\left(f^{\prime }\right)}^2-S\left(f^{\prime }+\frac{\eta }{2}{f}^{″}\right)+M^2(E_1-f^{\prime })+{\lambda }_1\theta +{\lambda }_2\varphi {\left|\right|}_{\infty },\hfill \\ \hfill {\left|\right|Res\left(\theta \right)\left|\right|}_{\infty }& =& \left|\right|\frac{1}{Pr}\left(1+\frac{4R}{3}\right){\theta }^{″}+f{\theta }^{\prime }-\frac{S\eta }{2}{\theta }^{\prime }-(2S+f^{\prime })\theta +Ec({\left(f^{″}\right)}^2+We{\left(f^{″}\right)}^3)+M^{2}Ec{(f^{\prime }-E_1)}^2{\left|\right|}_{\infty },\hfill \\ \hfill {\left|\right|Res\left(\varphi \right)\left|\right|}_{\infty }& =& \left|\right|\frac{1}{Sc}{\varphi }^{″}+\left(f-\frac{S\eta }{2}\right){\varphi }^{\prime }+(2S-f^{\prime }-\gamma )\varphi {\left|\right|}_{\infty }.\hfill \end{array}$$

Figure 2 shows the convergence and accuracy of the SQLM. We note that an increase in the number of iterations results in a decrease of the error infinity norm and method converges after six iterations. Also, the decrease in residual error infinity norms against the number of iterations confirms the accuracy of the numerical method. The SQLM achieves an accuracy of order ${10}^{-15}$ after five iterations, showing that the method is highly accurate.

Figure 3 shows the sensitivity of temperature and concentration profiles with the unsteadiness parameter S and the chemical reaction parameter $\gamma$, respectively. It is shown that for values of $S\ge 0$, the temperature profiles show the expected trend, satisfying boundary conditions (13) and (14). For values of $S<0$, the profiles behave unexpectedly. Similarly, the concentration profiles show the desired pattern for values of $\gamma >0$, satisfying the underlying boundary conditions (13) and (14). We remark that these chaotic and high sensitivity behaviour happen when the method has not yet reached the convergence stage. Once convergence has been reached, the behaviour becomes normal as expected. This can clearly be observed in subsequent figures.

The effects of the emerging physical parameters of importance in this study on the fluid velocity, temperature and concentration fields are shown in Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12. The effects of the Weissenberg number $We$ on the fluid velocity $f^{\prime }\left(\eta \right)$ and the fluid temperature $\theta \left(\eta \right)$ is shown in Figure 4. It is noted that the velocity profiles are decreasing functions of $We$. The Weissenberg number gives the ratio of relaxation time to specific process time. Increasing $We$ translates to a decrease in the specific process time which in turn will cause a decrease in the velocity component as well as the boundary layer thickness. Increasing the values of $We$ results in a decrease in the fluid temperature profiles.

Figure 5 shows the graphs of the dimensionless concentration profiles for different values of the Schmidt number $Sc$ and the chemical reaction parameter $\gamma$. We can observe that increasing $Sc$ results in a decrease in the fluid concentration. Physically, since $Sc$ is dependent on fluid mass diffusivity, high values of $Sc$ correspond to a decrease in mass diffusion hence the concentration profile is reduced. The concentration profiles are also decreased by increasing the chemical reaction. This is due to the fact that a higher chemical reaction results in a decrease in chemical molecular diffusivity leading to a reduction in mass diffusion. The decrease in concentration of the diffusing species will result in the thinning of the concentration boundary layer.

The influence of the magnetic field parameter on the fluid velocity and temperature is presented in Figure 6. We notice that increasing the values of the magnetic field parameter has an effect of slowing down the fluid motion. Physically, the presence of a magnetic field creates a drag-like force called the Lorentz force that opposes the motion of the fluid and resists the velocity of the fluid. The variation of the magnetic parameter M on $\theta \left(\eta \right)$ shows an opposite trend. We observe that the temperature profiles increase as the magnetic field parameter increases. This is because, as mentioned earlier on, an increase in M reduces the magnitude of the velocity profiles in the boundary layer, hence will cause a rise in the temperature in the boundary layer.

Figure 7 shows that an increase in the magnetic parameter causes an increase in the concentration field. We noted earlier that increasing the magnetic parameter reduces the velocity profiles in the boundary which in turn induces the diffusion of the particles in the boundary layer. The increment of the thermal buoyancy parameter ${\lambda }_1$ causes a decline in the fluid concentration.

The effects of the thermal radiation parameter R and the Prandtl number $Pr$ on the temperature profiles are depicted on Figure 8. We can see that the fluid temperature increases with increasing values of the thermal radiation parameter R. $Pr$ is inversely proportional to thermal diffusivity. Increasing the values of $Pr$ causes a decrease in the temperature field. Physically, high values of $Pr$ are associated with low thermal conductivity which reduces conduction and hence the thermal boundary layer resulting in a decrease in fluid temperature.

As shown in Figure 9, the temperature profiles decrease with the increase in the values of the buoyancy parameters ${\lambda }_1$ and ${\lambda }_2$. It is clearly noted that an increase in thermal buoyancy parameter causes a decrease in the thermal boundary layer thickness and consequently the fluid temperature decreases due to the buoyancy effect.

From Figure 10 it can be seen that the velocity and temperature profiles decrease as the unsteadiness parameter S increases. The effect of $E_1$ on the fluid velocity is the opposite of that of S. Figure 11 depicts that increasing the values of the local electric parameter $E_1$ results in an increase in both the velocity and concentration fields. Figure 12 shows the effects of increasing the values of Eckert number $E_c$ on the temperature profiles and ${\lambda }_1$ on velocity profiles. The Eckert number characterize the influence of self heating of a fluid as a consequence of dissipation effects. Viscous dissipation due to internal friction of the fluid will cause an increase in the fluid temperature. Increasing the value of the thermal buoyancy parameter ${\lambda }_1$ leads to an increase in the velocity profiles.

The examination of influence of the emerging parameters on the physical quantities of engineering importance, namely the skin friction coefficient, Nusselt number and Sherwood number is depicted in

Table 1. The numerical values of $C_{f}R{e}_x^{\frac{1}{2}}$ for different values of the parameters $We$, S, M, $E_1$, ${\lambda }_1$, ${\lambda }_2$, $Pr$, R, $E_c$, $\gamma$ and $Sc$ when $\eta =0.5$.

S	M	$E_1$	${\lambda }_1$	${\lambda }_2$	Pr	$\gamma$	Sc	R	$E_c$	We	$-C_f{Re}_x^{\frac{1}{2}}$
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	0.1	0.1	1.078951883204829
0.2											1.104506706276102
0.3											1.129271687273391
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	0.1	0.1	1.078951883204829
	0.8										1.100265152109772
	0.9										1.120409610249936
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	0.1	0.1	1.078951883204829
		0.2									1.033110225329453
		0.3									0.987166455664011
0.1	0.7	0.1	0.2	0.1	1.0	0.5	0.2	0.1	0.1	0.1	1.046544533993518
			0.4								0.981494769194006
			0.6								0.916019062597979
0.1	0.7	0.1	0.1	0.2	1.0	0.5	0.2	0.1	0.1	0.1	1.025760860799882
				0.4							0.921685018418928
				0.6							0.819005181818752
0.1	0.7	0.1	0.1	0.1	2.0	0.5	0.2	0.1	0.1	0.1	1.091693867285284
					3.0						1.103983070497612
					4.0						1.112467387854574
0.1	0.7	0.1	0.1	0.1	1.0	1.0	0.2	0.1	0.1	0.1	1.085984779187155
						2.0					1.090279447513868
						3.0					1.093114307914497
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.1	0.1	0.1	0.1	1.078951883204829
							0.15				1.081939366249072
							0.2				1.084395556441141
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	0.1	0.1	1.078951883204829
								0.4			1.076840518443974
								0.8			1.074763680141721
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	1.0	0.1	1.076389805314664
									2.0		1.071224585593461
									3.0		1.066003195448882
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	0.1	0.2	0.948884379906569
										0.4	0.620322356942535
										0.6	0.194330797445648

,

Table 2. The numerical values of $N{u}_{x}R{e}_x^{\frac{1}{2}}$ for different values of the parameters $We$, S, M, $E_1$, ${\lambda }_1$, ${\lambda }_2$, $Pr$, R, $E_c$, $\gamma$ and $Sc$ when $\eta =0.5$.

S	M	$E_1$	${\lambda }_1$	${\lambda }_2$	Pr	$\gamma$	Sc	R	$E_c$	We	${Nu}_x{Re}_x^{\frac{1}{2}}$
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	0.1	0.1	0.915860341811270
0.2											1.017901652122212
0.3											1.109121647214719
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	0.1	0.1	0.915860341811270
	0.8										0.895915260713988
	0.9										0.876756392212799
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	0.1	0.1	0.915860341811270
		0.2									0.977894447068897
		0.3									1.030403768637055
0.1	0.7	0.1	0.2	0.1	1.0	0.5	0.2	0.1	0.1	0.1	0.929041372189117
			0.4								0.953439839063926
			0.6								0.975582164474820
0.1	0.7	0.1	0.1	0.2	1.0	0.5	0.2	0.1	0.1	0.1	0.948908798318697
				0.4							1.001281793736204
				0.6							1.041796255284524
0.1	0.7	0.1	0.1	0.1	2.0	0.5	0.2	0.1	0.1	0.1	1.304303729032015
					3.0						1.597280137353339
					4.0						1.845769077311545
0.1	0.7	0.1	0.1	0.1	1.0	1.0	0.2	0.1	0.1	0.1	0.907214610540086
						2.0					0.902204998128932
						3.0					0.899080115679210
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.1	0.1	0.1	0.1	0.907214610540086
							0.15				0.903342484386425
							0.2				0.900478136478512
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	0.1	0.1	0.907214610540086
								0.4			0.992636450654252
								0.8			1.093677046062223
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	1.0	0.1	0.636725832004835
									1.5		0.367594034183057
									3.0		0.099808347796939
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	0.1	0.2	0.914843548830036
										0.4	0.936472243301972
										0.6	0.967724950819449

and

Table 3. The numerical values of $S{h}_{x}R{e}_x^{\frac{1}{2}}$ for different values of the parameters $We$, S, M, $E_1$, ${\lambda }_1$, ${\lambda }_2$, $Pr$, R, $E_c$, $\gamma$ and $Sc$ when $\eta =0.5$.

S	M	$E_1$	${\lambda }_1$	${\lambda }_2$	Pr	$\gamma$	Sc	R	$E_c$	We	${Sh}_x{Re}_x^{\frac{1}{2}}$
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	0.1	0.1	0.223598982783657
0.2											0.137162963897435
0.3											0.062147655092748
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	0.1	0.1	0.223598982783657
	0.8										0.221909175374051
	0.9										0.220302263370217
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	0.1	0.1	0.220302263370217
		0.2									0.244356594673201
		0.3									0.262718334286973
0.1	0.7	0.1	0.2	0.1	1.0	0.5	0.2	0.1	0.1	0.1	0.222825140750917
			0.4								0.227598540998689
			0.6								0.232071072563289
0.1	0.7	0.1	0.1	0.2	1.0	0.5	0.2	0.1	0.1	0.1	0.232796360927152
				0.4							0.252429377039816
				0.6							0.267886177341174
0.1	0.7	0.1	0.1	0.1	2.0	0.5	0.2	0.1	0.1	0.1	0.218996596175997
					3.0						0.218382226977623
					4.0						0.218017429115231
0.1	0.7	0.1	0.1	0.1	1.0	1.0	0.2	0.1	0.1	0.1	0.220049536513613
						2.0					1.090279447513868
						3.0					1.093114307914497
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.1	0.1	0.1	0.1	0.220049536513613
							0.15				0.286846328785100
							0.2				0.351830465481003
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	0.1	0.1	0.351830465481003
								0.4			0.352601074008495
								0.8			0.353544944638484
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	1.0	0.1	0.352768392930626
									2.0		0.353807745690575
									3.0		0.354844653039502
0.1	0.7	0.1	0.1	0.1	1.0	0.5	0.2	0.1	0.1	0.2	0.349345275097305
										0.4	0.343567754574247
										0.6	0.337874422565249

. Table 1 shows that the skin friction increases as the values of S, M, $Pr$, $\gamma$, $Sc$ and $E_c$ are increased. The opposite trend is observed for the parameters $E_1$, ${\lambda }_1$, R and $We$. In Table 2, we can see that the heat transfer rate increased by the increasing values of S, $E_1$, ${\lambda }_1$, ${\lambda }_2$, $Pr$ and R. The heat transfer rate is depressed as the the values of M, $\gamma$, $Sc$ and $E_c$. Table 3 shows that the coefficient of the Sherwood number is an increasing function of $E_1$, ${\lambda }_1$, ${\lambda }_2$, $\gamma$, $Sc$, R and $E_c$ and a decreasing function of S, M, $Pr$ and $We$.

5. Conclusions

A Williamson fluid is characteristic of a non-Newtonian fluid model with shear thinning property. The model constitute a coupled system of nonlinear partial differential equations. The transformed coupled system of dimensionless non-linear ordinary differential equations was successfully solved using the spectral quas-linearization method. The use of error infinity norms and residual error infinity norms confirmed that the numerical technique is convergent and accurate, respectively. In this present study, we numerically analyzed the effects of viscous dissipation, thermal radiation, chemical reaction and uniform magnetic field on the unsteady boundary layer flow of an electrically conducting Williamson fluid over a stretching sheet. We conclude the most important findings of the study as follows:

• The SQLM is a very efficient and accurate method.
• The fluid velocity and the momentum boundary layer decrease with respective, increases in the Williamson parameter, unsteadiness parameter, magnetic parameter, Eckert number as well as the Prandtl and Schmidt numbers.
• The fluid velocity and the momentum boundary layer increase with increasing values of the electric parameter, buoyancy parameters, thermal radiation and chemical reaction parameter.
• The fluid temperature increases as the values of the magnetic parameter, thermal radiation parameter, electrical parameter and Eckert number increase.
• The fluid temperature is a decreasing function of the buoyancy parameter, Prandtl number, unsteadiness parameter as well as the Williamson number.
• The stretching parameter, chemical reaction parameter, suction, Schmidt number, buoyancy parameters and the Williamson number were found to reduce the concentration profiles.
• The concentration was observed to be increasing as the values of the magnetic parameter, injection and Eckert number increase.
• The skin friction increases with the increase of the unsteadiness parameter, magnetic parameter, Prandtl number, Schmidt number, chemical reaction parameter, and thermal radiation parameter.
• However, the skin friction decreases with increasing values of the Eckert number, buoyancy parameters, thermal radiation and the Williamson number.
• The wall temperature gradient decreases with the increasing values of the Williamson number, suction, magnetic parameter, chemical reaction parameter, Schmidt number and Eckert number.
• The study observed that the Nusselt number increases with the increase of the unsteadiness parameter, electric parameter, buoyancy parameters, Prandtl number, thermal radiation parameter, and the Williamson number.
• The unsteadiness parameter, magnetic parameter, the Prandtl number and the Williamson number cause the wall concentration gradient to decrease.
• Lastly, the Sherwood number increases as the values of the electric parameter, buoyancy parameters, chemical reaction, Schmidt number, thermal radiation and Eckert number increase.

The joy of the current paper lies in the accuracy and fast convergence of the used numerical technique. Therefore, the method may also be valid for other complex nonlinear boundary value problems even with chaotic behaviour. We recommend that these results can be used as a standard example for other applications in engineering and applied sciences.