Mixed Convection Flow over an Elastic, Porous Surface with Viscous Dissipation: A Robust Spectral Computational Approach

Abstract

A novel computational approach is developed to investigate the mixed convection, boundary layer flow over a nonlinear elastic (stretching or shrinking) surface. The viscous fluid is electrically conducting, incompressible, and propagating through a porous medium. The consequences of viscous dissipation, Joule heating, and heat sink/source of the volumetric rate of heat generation are also included in the energy balance equation. In order to formulate the mathematical modeling, a similarity analysis is performed. The numerical solution of nonlinear differential equations is accomplished through the use of a robust computational approach, which is identified as the Spectral Local Linearization Method (SLLM). The computational findings reported in this study show that, in addition to being simple to establish and numerically implement, the proposed method is very reliable in that it converges rapidly to achieve a specified goal and is more effective in resolving very complex models of nonlinear boundary value problems. In order to ensure the convergence of the proposed SLLM method, the Gauss–Seidel approach is used. The SLLM’s reliability and numerical stability can be optimized even more using Gauss–Seidel approach. The computational results for different emerging parameters are computed to show the behavior of velocity profile, skin friction coefficient, temperature profile, and Nusselt number. To evaluate the accuracy and the convergence of the obtained results, a comparison between the proposed approach and the bvp4c (built-in command in Matlab) method is presented. The Matlab software, which is used to generate machine time for executing the SLLM code, is also displayed in a table.

1. Introduction

Heat transfer over an elastic surface has attracted the interest of researchers around the globe because of its numerous applications in manufacturing and commercial processes such as glass blowing, polymer extrusion, hot rolling, metal spinning, and wire drawing [1,2]. Many chemical engineering processes [3,4], such as metallurgical operations and polymer extrusion processes involve the cooling of a molten liquid that is stretched into a cooling system, have made use of the momentum and heat transfer of boundary layer flow across the stretching surface. In these applications, continuous strips of filaments are cooled by drawing them through a quiescent fluid. The process of drawing stretches the filaments. Rahmati et al. [5] investigated the non-Newtonian nanofluid flow under the impact of no-slip and slip conditions. Sene [6] discussed the effects of Newtonian heating using a second-grade fluid with Caputo fractional derivative.

The heat transfer process through porous media is applicable in a variety of fields ranging from Ceramic Engineering to Geophysiscs as it appears below in

Table 1. Applications of heat transfer through porous media.

Ceramic engineering [7]	Simultaneous mass and heat transfer in disordered media occurs during the burnout/drying of the binder system from green compacts during the colloidal process of ceramics.
Chemical engineering [8]	During the interim eventual and storage of nuclear waste, as well as in packet bed reactors.
Ground water hyrology [9]	The investigation of seepage water through river beds and underground water resources.
Industrial engineering [10]	The primary goal of filtration analysis is to examine the movement of fluid through a porous medium, leaving behind unwanted material. As a result of the mass deposition, the porous medium is constantly changing and altering the system’s pressure drop properties.
Mechanical engineering [11]	To achieve effective insulation, solid conduction must be minimized, porosity must be maximized to reduce effective thermal conductivity, and free convection must be suppressed. The same concept is useful when producing high-performance insulators for cryogenic containers.
Petroleum engineering [12]	For oil recovery mechanisms.
Geophysics [13]	In the analysis of geo-pressurized reservoirs, and extraction of geothermal energy.

.

All of these applications provide an incentive to further investigate the mixed convection flow in a porous medium with a non-linear elastic surface.

Among the previous studies on the subject, Sakiadis [14] investigated the behavior of the boundary layer on a continuous, stable subsurface, including both laminar and turbulent flow. Crane [15] further extended the work of Sakiadis and investigated a Newtonian fluid induced by a linearly stretched flat surface. He derived an exact solution for the problem in closed form. More recently, Rosca and Pop [16] examined the heat transfer across a stretching/shrinking sheet using a second order slip flow model. They used the Matlab’s bvp4c function for their numerical solutions. The effect of a heat source/sink on a viscous fluid passing through a porous medium was studied by Swain et al. [17]. Othman et al. [18] used the shooting method, which generates numerical solutions for the stagnation point of a steady 2-D mixed convection flow across a stretching/shrinking surface. Kumar et al. [19] also explored the role of Joule heating on the thermal boundary layer flow and heat transfer over a stretching surface. The effects of buoyancy on the two-dimensional mixed convection and Casson fluid along a linearly stretching sheet was studied by Gangadhar et al. [20]. Their numerical solutions were obtained by a spectral relaxation analysis. Prabha et al. [21] used the LTNE model to evaluate the 2-D mixed convection boundary layer flow through a porous medium. Badruddin et al. [22] provided a brief overview of heat transfer through porous materials. Ali et al. [23] examined the flow over a stretching sheet with a porous surface. Ilya et al. [24] used a Finite-Difference approach to investigate the occurrence of mixed convection flow along the surface of a cone immersed in a porous medium in the presence of a heat source/sink, a magnetic field, and density difference that induce buoyancy. Adeniyan et al. [25] studied the similarity solutions within the context of convective boundary layer flow and heat transfer around a stretching hollow cylinder, immersed and filled with a viscoelastic fluid. To obtain numerical results, they used the Runge–Kutta–Fehlberg integration method together with the shooting method.

Viscous dissipation is important in many areas: significant temperature rises have been experimentally observed in polymer processes, such as fast extrusion or injection molding. Aerodynamic heating in a thin boundary layer around high-speed rail and aircraft elevates their surface temperature. In a completely different application, the dissipation function is used to express the viscosity of dilute suspensions [26]. Viscous dissipation is also observed in convection mechanisms [27,28,29] for various devices that are operate at high relative velocities or at high rotational velocities. The impact of viscous dissipation can be also seen in planetary processes because it causes stronger gravitational fields with large masses of gas and larger planets.

The field of magnetohydrodynamics (MHD) flow and heat transfer has applications in a wide range of industries, ecological systems, and geophysical explorations. For example, paper manufacturing, glass-fiber manufacturing, and the fabrication of polymer sheets from dies, continuous casting and wire drawing are all examples of industrial processes related to this field. Because of their industrial applications, such flows have been addressed by a variety of researchers, including: Jafar et al. [30] evaluated the steady boundary layer flow of incompressible and viscous fluids under the effects of MHD. Dessie et al. [31] studied the relationship of viscous dissipation as a heat source/sink in MHD flow and heat transfer via porous media across a stretched sheet. The solutions of the governing equations in this study were calculated by employing Lie’s group transformations. Bibi et al. [32] used the shooting technique to numerically investigate the unsteady MHD flow of a Williamson fluid with heat transfer on a pervious stretchable sheet. Alarifi et al. [33] numerically assessed the MHD boundary layer flow over a vertical stretched surface, in the presence of a heat source/sink and a magnetic force. They used the RKFI approach to identify numerical solutions to the nonlinear differential equations. Also, Swain et al. [34] explored the effects of Joule heating and viscous dissipation on MHD flow and heat transfer via porous material. Megahed et al. [35] described an unsteady magnetized fluid flow with thermal radiation effects and heat flux conditions. To obtain numerical solutions, they combined the shooting method with a Runge–Kutta algorithm. Sarda et al. [36] considered a non-Newtonian fluid propagating over a stretching surface and non-equilibrium thermal conditions to investigate the effect of magnetic fields on heat transfer. They also used a combination of shooting technique and the fourth-fifth-order Runge–Kutta–Fehlberg (RKF) method. Zhou et al. [37] examined the role of slip on a stretchable surface with a non-uniform heat source using Casson fluid. They used the bvp4c command to develop solutions to the nonlinear problem. Ullah et al. [38] used a new stochastic method to investigate the MHD boundary layer flow. Recently, Bhatti et al. [39] used a Lie group analysis and a successive linearization technique to solve a nonlinear Jeffrey fluid model with mass transport.

The primary goal of this paper is to present a rigorous computational method for solving nonlinear boundary value problems. The proposed methodology is a simple, precise, and convergent algorithm and is known as the spectral local linearization method, abbreviated as SLLM. This method involves linearizing and decoupling equations using a univariate linearization approach and spectral collocation linearization. The proposed algorithm’s main distinguishing feature is that it converts a large coupled system of equations into a set of smaller systems that can be solved sequentially in a computationally efficient manner. We applied this methodology to a nonlinear problem, such as two coupled nonlinear differential equations of mixed convection heat transfer. The results show that, when compared to other similar methods, the computed methodology is simple to develop, fast converging, efficient, accurate, and reliable.

Not enough emphasis has been given to the mixed convection fluid flow with viscous dissipation. The majority of the literature reported above used the bvp4c, a shooting method, the Runge–Kutta approach or other numerical and analytical techniques to solve nonlinear boundary layer flow problems. Spectral numerical methods, like SLLM, were not used to solve such problems. The primary goal of this research is to use a robust computational technique to assess the mixed convection fluid flow including magnetic and porous effects. The energy equation takes into account the simultaneous effects of viscous dissipation, heat source/sink, and Joule heating. Following a similarity analysis, the numerical solutions to the governing equations were obtained using the computational approach of Spectral Local Linearization. The Gauss–Seidel technique is also used to improve the convergence of this method. The suggested methodology’s machine time is also given in tabular form utilizing the Matlab software. Full numerical and graphical comparisons with the bvp4c [40] and previously reported studies are also provided in this paper.

2. Problem Description and Modeling

We consider a two-dimensional mixed convective flow of a viscous incompressible fluid induced by a nonlinear elastic sheet, while stretching or shrinking, with a porous medium. The flow is modeled in a Cartesian coordinate system, with the $x^{\prime }$-coordinate along the moving surface and the $y^{\prime }$-coordinate in the normal to the surface direction. The elastic sheet is at rest on the plane $y^{\prime }=0$, which also incorporates the heat sink/source with volumetric heat absorption or generation $Q_0=x^{\prime }{q}_0$ (where $q_0$ is constant). The nonlinear elastic sheet is deformed with a wall velocity $w_s\left(x^{\prime }\right)=c{x}^{\prime 2}$, where the constant parameter c is positive during stretching $(c>0)$ and negative during shrinking $(c<0)$. The elastic surface is permeable. The process is an injection when the velocity is $v_s<0$ and suction when the velocity is $v_s>0$. The fluid is electricity conducting when exposed to the influence of a uniform variable magnetic field $B_0=b_0\sqrt{x^{\prime }}$ and traveling through a porous material with variable permeability $k=k_0/x^{\prime }$. Furthermore, $T_w$ signifies temperature at the wall, while $T_{inf}$ indicates the ambient temperature distant from the surface. The governing equations for this physical problem are formulated using the standard boundary layer approximations [41]:

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial u^{\prime }}{\partial x^{\prime }}}+{\displaystyle \frac{\partial v^{\prime }}{\partial y^{\prime }}}=0,\end{array}$$

(1)

$$\begin{array}{c}\hfill u^{\prime }{\displaystyle \frac{\partial u^{\prime }}{\partial x^{\prime }}}+v^{\prime }{\displaystyle \frac{\partial u^{\prime }}{\partial y^{\prime }}}=\nu {\displaystyle \frac{{\partial }^2{u}^{\prime }}{{\partial y^{\prime }}^2}}+dg\chi \left(T-T_{inf}\right)-{\displaystyle \frac{\sigma {B_0}^2}{\rho }}{u}^{\prime }-{\displaystyle \frac{\nu }{k}}{u}^{\prime },\end{array}$$

(2)

where $u^{\prime }$ and $v^{\prime }$ represent the velocity components in $x^{\prime }$- and $y^{\prime }$- coordinates, $\rho$ denotes the density, $\nu$ denotes the kinematic viscosity, g is the gravitational acceleration, $\chi$ denotes the coefficient of the thermal expansion, $\mu$ is the fluid viscosity, $\sigma$ denotes the electrical conductivity, and d is a dummy parameter with values $0,+1$ or $-1$.

$$\begin{array}{c}\hfill u^{\prime }{\displaystyle \frac{\partial T}{\partial x^{\prime }}}+v^{\prime }{\displaystyle \frac{\partial T}{\partial y^{\prime }}}=\alpha {\displaystyle \frac{{\partial }^{2}T}{{\partial y^{\prime }}^2}}+{\displaystyle \frac{Q_0}{\rho c_p}}\left(T-T_{inf}\right)+\mathrm{V}.\mathrm{D}.,\end{array}$$

(3)

where $\alpha$ denotes the thermal diffusivity, $c_p$ denotes the specific heat capacity, and V.D. is the viscous dissipation, defined as follows:

$$\begin{array}{c}\hfill \mathrm{V}.\mathrm{D}.={\displaystyle \frac{1}{\rho c_p}}\left[\mu {\left({\displaystyle \frac{\partial u^{\prime }}{\partial y^{\prime }}}\right)}^2+\sigma {B_0}^2{u}^{\prime 2}+{\displaystyle \frac{\mu }{k}}{u}^{\prime 2}\right].\end{array}$$

(4)

The boundary conditions at $y^{\prime }=0$ are defined as:

$$\begin{array}{c}\hfill u^{\prime }=w_s\left(x^{\prime }\right)=c{x^{\prime }}^2,v^{\prime }={x^{\prime }}^{\frac{1}{2}}{v}_s,T=T_w\left(x^{\prime }\right)=T_{inf}+b{x}^{\prime 3},\end{array}$$

(5)

where b is a constant that indicates a hot plate when $b>0$ and a cool plate when $b<0$.

The boundary conditions at $y^{\prime }\to \infty$ are defined as:

$$\begin{array}{c}\hfill u^{\prime }\to 0,{\displaystyle \frac{\partial u^{\prime }}{\partial y^{\prime }}}\to 0,T\to T_{inf}.\end{array}$$

(6)

Similarity Analysis

The following transformations are introduced for the solution of the governing equations:

$$\begin{array}{cc}\hfill & \eta =y^{\prime }\sqrt{{\displaystyle \frac{3\left|c\right|x^{\prime }}{2\nu }}},u^{\prime }=\left|c\right|{x^{\prime }}^2{f}^{\prime }\left(\eta \right),\hfill \\ & v^{\prime }=-{\displaystyle \frac{1}{2}}\sqrt{{\displaystyle \frac{2}{3}}\nu \left|c\right|x^{\prime }}\left[3f\left(\eta \right)+\eta f^{\prime }\left(\eta \right)\right],T=b{x^{\prime }}^3\theta +T_{inf},\hfill \end{array}$$

(7)

and the velocity of the wall mass transfer becomes $v_s=-\sqrt{{\displaystyle \frac{3\nu \left|c\right|}{2}}}\varsigma$. Applying Equation (7) to the governing Equations (2)–(6), we derive the following set of differential equations

$$\begin{array}{c}\hfill f^{‴}+f{f}^{″}-{\displaystyle \frac{4}{3}}{f^{\prime }}^2-(\beta +\gamma )f^{\prime }+\Lambda \theta =0,\end{array}$$

(8)

$$\begin{array}{c}\hfill {\theta }^{″}+\Gamma \left[f{\theta }^{\prime }-2f^{\prime }\theta +\lambda {f^{″}}^2+\lambda (\beta +\gamma ){f^{\prime }}^2+\delta \theta \right]=0,\end{array}$$

(9)

where $\beta ={\displaystyle \frac{2\sigma {b_0}^2}{3\left|c\right|\rho }}$ denotes the magnetic parameter, $\gamma ={\displaystyle \frac{2\upsilon }{3k_0\left|c\right|}}$ denotes the permeability parameter, $\lambda ={\displaystyle \frac{w_s^2}{c_p(T_w-T_{inf})}}$ denotes the Eckert number, $\Gamma ={\displaystyle \frac{\nu }{\alpha }}$ denotes the Prandtl number, $\delta ={\displaystyle \frac{2}{3}}{\displaystyle \frac{q_0}{\rho c_p\left|c\right|}}$ denotes the local heat generation/absorption coefficient, $\Lambda =d{\displaystyle \frac{2}{3}}{\displaystyle \frac{G_r}{{\Re }^2}}$ denotes the mixed convention parameter, $G_r={\displaystyle \frac{g\chi (T_w-T_{inf}){x^{\prime }}^3}{{\upsilon }^2}}$ denotes the Grashof number, and $\Re ={\displaystyle \frac{w_s{x}^{\prime }}{\upsilon }}$ denotes the Reynolds number. Mixed convection occurs when the mechanisms of forced and natural convection work in tandem to transfer heat. The Grashof number (for natural convection) and the Reynolds number (for forced convection) are two dimensionless numbers that are frequently used to define the strength of each part of the convection [42,43].

The boundary conditions reduce to the following form:

$$\begin{array}{cc}\hfill & f\left(0\right)=\varsigma ,f^{\prime }\left(0\right)=\varphi ,\theta \left(0\right)=1,\hfill \\ & f^{\prime }\left(\infty \right)=0,\theta \left(\infty \right)=0,\hfill \end{array}$$

(10)

where $\varsigma$ denotes the wall transpiration, and $\varphi$ is the stretching/shrinking parameter. If $\varphi >0$ the surface is stretching and if $\varphi <0$ the surface shrinks.

The skin friction coefficient and the Nusselt number are expressed as follows:

$$\begin{array}{c}\hfill C_f={\displaystyle \frac{{\tau }_s}{\rho w_s^2}}{\left({\displaystyle \frac{\partial u^{\prime }}{\partial y^{\prime }}}\right)}_{y=0},N_u={\displaystyle \frac{x^{\prime }{q}_f}{\kappa \left(T_w-T_{inf}\right)}},\end{array}$$

(11)

where $\kappa$ represents the thermal conductivity. The wall shear stress ${\tau }_s$, and the surface heat flux $q_f$ are given by the equations:

$$\begin{array}{c}\hfill {\tau }_s=\mu {\left({\displaystyle \frac{\partial u^{\prime }}{\partial y^{\prime }}}\right)}_{y^{\prime }=0},q_f=-\kappa {\left({\displaystyle \frac{\partial T}{\partial y^{\prime }}}\right)}_{y^{\prime }=0}.\end{array}$$

(12)

Using Equations (7) and (11), we obtain the following expressions:

$$\begin{array}{c}\hfill {\Re }^{\frac{1}{2}}{C}_f=\sqrt{{\displaystyle \frac{3}{2}}}{f}^{″}\left(0\right),{\Re }^{-\frac{1}{2}}{N}_u=-\sqrt{{\displaystyle \frac{3}{2}}}{\theta }^{\prime }\left(0\right).\end{array}$$

(13)

It can be seen that the resulting Equations (8) and (9) are nonlinear and coupled. Exact and closed-form solutions to these differential equations are impossible to derive. Therefore, in the following section, we will present a robust computational method to solve the system of equations.

3. Implementation of Numerical Method

3.1. Basic Steps to Apply Spectral Local Linearization Method

Let us consider a system of differential equations $\mathcal{F}=\left[f_1\left(\eta \right),f_2\left(\eta \right),\dots ,f_m\left(\eta \right)\right]$ that fulfill the condition:

$$\begin{array}{c}\hfill {\mathcal{L}}_n+{\mathcal{N}}_n={\mathcal{H}}_n,n=1,2,\dots ,m,\end{array}$$

(14)

where m denotes the number of differential equations. Each ${\mathcal{H}}_n$ is a function of $\eta \in \left(\mathcal{A},\mathcal{B}\right)$, and ${\mathcal{L}}_n,{\mathcal{N}}_n$ represent the system’s linear and nonlinear components, respectively.

To continue the iteration process, the local linearization of ${\mathcal{N}}_n$ at approximately ${\mathcal{F}}_{n,t}$ (latest iteration) is employed for the nth nonlinear equation by assuming that all other ${\mathcal{F}}_{l,t}(l\ne n)$ values are known.

$$\begin{array}{cc}\hfill {\mathcal{N}}_n\left[f_1\left(\eta \right),f_2\left(\eta \right),\dots ,f_m\left(\eta \right)\right]& ={\mathcal{N}}_n\left[f_{1,t},f_{2,t},\dots ,f_{m,t}\right]\hfill \\ & +\nabla {\mathcal{N}}_n\left[f_{1,t},f_{2,t},\dots ,f_{m,t}\right](f_n-f_{n,t}).\hfill \end{array}$$

(15)

As a result, for the current iteration with $f_n=f_{n,t+1}$, Equation (14) becomes

$$\begin{array}{cc}\hfill & {\mathcal{L}}_n\left[f_{1,t},f_{2,t},\dots f_{m,t+1}\right]+\nabla {\mathcal{N}}_n\left[f_{1,t},f_{2,t},\dots ,f_{m,t}\right]f_{1,t+1}\hfill \\ & ={\mathcal{H}}_n+\nabla {\mathcal{N}}_n\left[f_{1,t},f_{2,t},\dots ,f_{m,t}\right]f_{1,t}-{\mathcal{N}}_n\left[f_{1,t},f_{2,t},\dots ,f_{m,t}\right],\hfill \end{array}$$

(16)

$$\begin{array}{c}\hfill ⋮\\ \hfill ⋮\end{array}$$

$$\begin{array}{cc}\hfill & {\mathcal{L}}_m\left[f_{1,t},f_{2,t},\dots ,f_{m,t+1}\right]+\nabla {\mathcal{N}}_m\left[f_{1,t},f_{2,t},\dots ,f_{m,t}\right]f_{m,t+1}\hfill \\ & ={\mathcal{H}}_m+\nabla {\mathcal{N}}_m\left[f_{1,t},f_{2,t},\dots ,f_{m,t}\right]f_{m,t}-{\mathcal{N}}_m\left[f_{1,t},f_{2,t},\dots ,f_{m,t}\right],\hfill \end{array}$$

(17)

where the nth equation $\left[f_{1,t},f_{2,t},\dots ,f_{m,t}\right]\equiv \left[f_{1,t+1},f_{2,t+1},\dots ,f_{n-1,t+1},f_{n,t},\dots ,f_{m,t}\right]$.

Traditionally, the SLLM is an iterative procedure for solving differential equations that begins with an initial approximation $f_0$ and then is implemented iteratively, generating new approximations $f_1,f_2,f_3,\dots ,$ where ${\mathcal{F}}_t=[f_{1,t},f_{2,t},\dots f_{m,t}]$ for each $t=0,1,2,\dots$. Once the nonlinear components ${\mathcal{N}}_n$ have been linearized, Equations (16) and (17) can be numerically solved utilizing the Chebyshev spectral collocation approach.

The nth differential Equation (14) after the first $t+1$ iterations can be expressed in this manner:

$$\begin{array}{c}\hfill {\mathcal{L}}_{n+1}+{\mathcal{N}}_{n+1}={\mathcal{H}}_n.\end{array}$$

(18)

Taylor series can be used to linearize the nonlinear components. For instance

$$\begin{array}{c}\hfill {\mathcal{N}}_n{|}_{t+1}={\mathcal{N}}_n{{|}_t+\nabla {\mathcal{N}}_n|}_t\left[{\mathcal{V}}_{t+1}-{\mathcal{V}}_t\right],\end{array}$$

(19)

where ${\mathcal{V}}_t$ denotes the n-tuples of $F_{n,t}$ and its differentials.

It is now possible to use Equations (18) and (19) into Equation (14) to obtain

$$\begin{array}{c}\hfill {\mathcal{L}}_n{|}_{t+1}+\nabla {\mathcal{N}}_n{|}_t{\mathcal{V}}_{t+1}={\mathcal{H}}_n+\nabla {\mathcal{N}}_n{{|}_t{\mathcal{V}}_t-{\mathcal{N}}_n|}_t.\end{array}$$

(20)

3.2. Implementation of the Spectral Local Linearization Method

To develop the SLLM, the current system of differential equations must be reduced in order, hence, we assume that $f^{\prime }=h$, to obtain the transformed Equations (8) and (9) in the form:

$$\begin{array}{c}\hfill h^{″}+f{h}^{\prime }-\left(\beta +\gamma \right)h-{\displaystyle \frac{4}{3}}{h}^2+\Lambda \theta =0,\end{array}$$

(21)

$$\begin{array}{c}\hfill {\theta }^{″}+\Gamma \left[f{\theta }^{\prime }-2h\theta +\lambda {h^{\prime }}^2+\lambda \left(\beta +\gamma \right)h^2+\delta \theta \right].\end{array}$$

(22)

We linearized the non-linear term $h^2$ using a Taylor series expansion, which is described in more detail below:

$$\begin{array}{c}\hfill h_{t+1}^2=2h_t\left(h_{t+1}-h_t\right)+h_t^2=2h_{t+1}{h}_t-h_t^2,\end{array}$$

(23)

where the component with the subscript $t+1$ represents the current approximated value and the component with the subscript t represents the previous value. Equation (23) is substituted in Equation (21), and the nonlinear system is decoupled by applying the Gauss–Seidel relaxation method to obtain the desired result. Thus, we have:

$$\begin{array}{c}\hfill {f^{\prime }}_{t+1}=h_t,\end{array}$$

(24)

$$\begin{array}{c}\hfill {h^{″}}_{t+1}+f_t{h^{\prime }}_{t+1}-\left(\beta +\gamma \right)h_{t+1}-{\displaystyle \frac{8}{3}}{h}_t{h}_{t+1}=-{\displaystyle \frac{4}{3}}{h_t}^2-\Lambda {\theta }_t,\end{array}$$

(25)

$$\begin{array}{c}\hfill {{\theta }^{″}}_{t+1}+\Gamma [f_t{{\theta }^{\prime }}_{t+1}-2h_{t+1}{\theta }_{t+1}+\delta {\theta }_{t+1}]=-\lambda \left[{h^{\prime }}_{t+1}^2+\left(\beta +\gamma \right)h_{t+1}^2\right].\end{array}$$

(26)

The pertinent boundary conditions are transformed as follows:

$$\begin{array}{c}\hfill f_{t+1}\left(0\right)=\varsigma ,h_{t+1}\left(0\right)=\varphi ,{\theta }_{t+1}\left(0\right)=1,\end{array}$$

(27)

$$\begin{array}{c}\hfill h_{t+1}\left(\infty \right)=0,{\theta }_{t+1}\left(\infty \right)=0.\end{array}$$

(28)

We write Equations (24)–(26) in a concise form as:

$$\begin{array}{c}\hfill {f^{\prime }}_{t+1}={\mathcal{E}}_{00},\end{array}$$

(29)

$$\begin{array}{c}\hfill {h^{″}}_{t+1}+{\mathcal{E}}_{11}{h^{\prime }}_{t+1}-{\mathcal{E}}_{13}{h}_{t+1}-{\mathcal{E}}_{12}{h}_{t+1}={\mathcal{E}}_{1,t},\end{array}$$

(30)

$$\begin{array}{c}\hfill {{\theta }^{″}}_{t+1}+\Gamma \left({\mathcal{E}}_{11}{{\theta }^{\prime }}_{t+1}-2h_{t+1}{\theta }_{t+1}+\delta {\theta }_{t+1}\right)={\mathcal{E}}_{2,t},\end{array}$$

(31)

where

$$\begin{array}{cc}\hfill & {\mathcal{E}}_{00}=h_t,{\mathcal{E}}_{11}=f_t,{\mathcal{E}}_{12}={\displaystyle \frac{8}{3}}{h}_t,{\mathcal{E}}_{13}=\left(\beta +\gamma \right),\hfill \\ & {\mathcal{E}}_{1,t}=-{\displaystyle \frac{4}{3}}{h_t}^2-\Lambda {\theta }_t,{\mathcal{E}}_{2,t}=-\lambda \left[{h^{\prime }}_{t+1}^2+\left(\beta +\gamma \right)h_{t+1}^2\right].\hfill \end{array}$$

(32)

Now, we apply the Chebyshev spectral collocation method to the system of Equations (24)–(26), where the differentiation matrix [44] $\mathcal{D}={\displaystyle \frac{2}{\mathfrak{l}}}D$ is used to approximate the derivatives of the unknown variables in the above equations. The physical region of the problem $[0,\infty )$ is truncated to $[0,\mathfrak{l}]$ to allow the numerical implementation using the spectral technique, where $\mathfrak{l}$ is stipulated to be sufficiently high. Our new system is written as:

$$\begin{array}{c}\hfill {\mathcal{D}f}_{t+1}=h_t,\end{array}$$

(33)

$$\begin{array}{c}\hfill \{{\mathcal{D}}^2+\mathrm{diag}\left[{\mathcal{E}}_{11}\right]\mathcal{D}-{\mathcal{E}}_{13}I-{\mathcal{E}}_{12}I\}{h}_{t+1}={\mathcal{E}}_{1,t},\end{array}$$

(34)

$$\begin{array}{c}\hfill \{{\mathcal{D}}^2+\Gamma \left(\mathrm{diag}\left[{\mathcal{E}}_{11}\right]\mathcal{D}-2\mathrm{diag}\left[h_{t+1}\right]+\delta I\right)\}{\theta }_{t+1}={\mathcal{E}}_{2,t}.\end{array}$$

(35)

The corresponding boundary conditions are:

$$\begin{array}{c}\hfill f_{t+1}\left({\eta }_N\right)=\varsigma ,h_{t+1}\left({\eta }_N\right)=\varphi ,{\theta }_{t+1}\left({\eta }_N\right)=1,\end{array}$$

(36)

$$\begin{array}{c}\hfill h_{t+1}\left({\eta }_0\right)=0,{\theta }_{t+1}\left({\eta }_0\right)=0.\end{array}$$

(37)

This system of equations may be expressed in a more concise form as follows:

$$\begin{array}{c}\hfill {\mathcal{B}}_1{f}_{t+1}={\mathcal{C}}_1,\end{array}$$

(38)

$$\begin{array}{c}\hfill {\mathcal{B}}_2{h}_{t+1}={\mathcal{C}}_2,\end{array}$$

(39)

$$\begin{array}{c}\hfill {\mathcal{B}}_3{\theta }_{t+1}={\mathcal{C}}_3,\end{array}$$

(40)

where

$$\begin{array}{c}\hfill {\mathcal{B}}_1=\mathcal{D},{\mathcal{C}}_1=h_t,\end{array}$$

(41)

$$\begin{array}{c}\hfill {\mathcal{B}}_2={\mathcal{D}}^2+\mathrm{diag}\left[{\mathcal{E}}_{11}\right]\mathcal{D}-{\mathcal{E}}_{13}I-{\mathcal{E}}_{12}I,{\mathcal{C}}_2={\mathcal{E}}_{1,t},\end{array}$$

(42)

$$\begin{array}{c}\hfill {\mathcal{B}}_3={\mathcal{D}}^2+\Gamma \left(\mathrm{diag}\left[{\mathcal{E}}_{11}\right]\mathcal{D}-2\mathrm{diag}\left[h_{t+1}\right]+\delta I\right),{\mathcal{C}}_3={\mathcal{D}}_{2,t},\end{array}$$

(43)

where the diagonal matrices “diag” are defined as follows:

$$\begin{array}{c}\hfill \mathrm{diag}\left[{\mathcal{E}}_{11}\right]=\left[\begin{array}{ccc}{\mathcal{E}}_{11}\left({\eta }_0\right)& \dots & .\\ ⋮& \ddots & ⋮\\ .& \dots & {\mathcal{E}}_{11}\left({\eta }_N\right)\end{array}\right],\mathrm{diag}\left[{\mathcal{E}}_{12}\right]=\left[\begin{array}{ccc}{\mathcal{E}}_{12}\left({\eta }_0\right)& \dots & .\\ ⋮& \ddots & ⋮\\ .& \dots & {\mathcal{E}}_{12}\left({\eta }_N\right)\end{array}\right],\end{array}$$

(44)

$$\begin{array}{c}\hfill \mathrm{diag}\left[{\mathcal{E}}_{1,t}\right]=\left[\begin{array}{c}{\mathcal{E}}_{1,t}\left({\eta }_0\right)\\ ⋮\\ {\mathcal{E}}_{1,t}\left({\eta }_N\right)\end{array}\right],\mathrm{diag}\left[{\mathcal{E}}_{2,t}\right]=\left[\begin{array}{c}{\mathcal{E}}_{2,t}\left({\eta }_0\right)\\ ⋮\\ {\mathcal{E}}_{2,t}\left({\eta }_N\right)\end{array}\right],\end{array}$$

(45)

and

$$\begin{array}{c}\hfill f_{t+1}={[f\left({\eta }_0\right),f\left({\eta }_1\right),\dots ,f\left({\eta }_N\right)]}^T,\end{array}$$

(46)

$$\begin{array}{c}\hfill h_{t+1}={[h\left({\eta }_0\right),h\left({\eta }_1\right),\dots ,h\left({\eta }_N\right)]}^T,\end{array}$$

(47)

$$\begin{array}{c}\hfill {\theta }_{t+1}={[\theta \left({\eta }_0\right),\theta \left({\eta }_1\right),\dots ,\left({\eta }_N\right)]}^T,\end{array}$$

(48)

where I represents the identity matrix of order $\left(N+1\right)\times \left(N+1\right)$. Equations (46)–(48) represents the vectors of size $\left(N\times 1\right)\times 1$.

In light of Equation (10), the following initial guesses are chosen for the implementation of the method:

$$\begin{array}{c}\hfill f_0\left({\eta }_0\right)=\varsigma +\varphi \left(1-exp[-\eta ]\right),h_0\left({\eta }_0\right)=\varphi exp[-\eta ],{\theta }_0\left({\eta }_0\right)=exp[-\eta ].\end{array}$$

(49)

These initial stipulations are used to obtain further approximations of $f_t,h_t,{\theta }_t$ for each $t=1,2,3,\dots$ by using the suggested technique.

3.3. Convergence Analysis

The Gauss–Seidel technique with successive over relaxation parameter is often used to accelerate the convergence of the linear system of equations. In this case, an equivalent approach is utilized to accelerate the convergence rate for SLLM. At the $(t+1)$th iteration, the SLLM approach is used to solve the function $\mathcal{Z}$.

$$\begin{array}{c}\hfill {\mathcal{B}}_1{\mathcal{Z}}_{t+1}={\mathcal{C}}_1.\end{array}$$

(50)

The new mode of the SLLM approach is then represented by rewriting the equation as:

$$\begin{array}{c}\hfill {\mathcal{B}}_1{\mathcal{Z}}_{t+1}=\left(1-\Omega \right){\mathcal{B}}_1{\mathcal{Z}}_t+\Omega {\mathcal{C}}_1.\end{array}$$

(51)

Here $\Omega$ denotes the convergence parameter. If the value of $\Omega$ is between 1 and 2, the current value has a higher weight than the prior value. The current value must also be close to the old value. This procedure will help to improve the convergent system even more. This form of change is known as over relaxation, and the process is known as the SOR approach (successive over relaxation). The improved SLLM approach significantly contributes to the improvement of the accuracy and numerical efficiency.

4. Numerical Results and Discussion

The numerical solution for the velocity profile, the skin friction coefficient, the Nusselt number, and the temperature profile are discussed in this section. The Matlab software is used to plot all of the numerical results graphically. We chose the appropriate parametric values for the graphical results such as: $N=100,\beta =0.2,\gamma =0.2,\delta =0.1,\Gamma =1,\Lambda =0.1,\lambda =0.5,\varsigma =0,\varphi =1$.

Table 2. Machine time to execute the Matlab code on intel core i7.

Function Name	Calls	Total Time (s)	Self Time (s)
SLLM code	1	0.284	0.202
Newplot	2	0.033	0.013
xlabel	2	0.029	0.023
newplot > ObserveAxesNextPlot	2	0.017	0.004
cla	2	0.013	0.007
hold	2	0.011	0.009
ylabel	2	0.010	0.008
convertStringToCharArgs	4	0.007	0.005
graphics/private/clo	2	0.005	0.005
gobjects	4	0.003	0.003
convertStringsToChars	4	0.002	0.001
markfigure	2	0.001	0.001
graphics/private/claNotify	2	0.001	0.001
graph2d/private/labelcheck	4	0.001	0.001
convertStringsToChars > convertStrings	16	0.001	0.001
axescheck	2	0.000	0.000
newplot > ObserveFigureNextPlo	2	0.000	0.000

shows the machine time (s), which indicates the execution time of numerical coding in Matlab.

Table 3. Impact of varying N on skin friction coefficient and Nusselt number.

Iterations	N	Skin Friction Coefficient	Nusselt Number
4	5	$-1.15627$	$1.17653$
4	10	$-1.29290$	$1.43602$
4	100	$-1.30162$	$1.44072$
−	bvp4c	$-1.30162$	$1.44072$

shows the effect of N on the skin friction coefficient and the Nusselt number. This table indicates that the desired level of accuracy (up to five decimal places) can be obtained when $N=100$ with four iterations. In addition,

Table 4. Numerical comparison between the present numerical method and bvp4c for Skin friction against different values of emerging parameters.

$\mathit{\beta }$	$\mathit{\gamma }$	$\mathbf{\Lambda }$	$\mathbf{\Gamma }$	$\mathit{\delta }$	$\mathit{\lambda }$	$\mathit{\varsigma }$	$\mathit{\varphi }$	SLLM	bvp4c
0.2	0.2	0.5	1	0.1	0.2	0	1	$-1.30162$	$-1.30162$
1								$-1.65671$	$-1.65671$
1.5								$-1.85064$	$-1.85064$
	0.5							$-1.44264$	$-1.44264$
	0.8							$-1.57383$	$-1.57383$
	1.2							$-1.73640$	$-1.73640$
		0.4						$-1.34991$	$-1.34991$
		1.0						$-1.07167$	$-1.07167$
		1.5						$-0.85505$	$-0.85505$
			0.71					$-1.27273$	$-1.27273$
			2					$-1.40238$	$-1.40238$
			3					$-1.44594$	$-1.44594$
				$-0.1$				$-1.31463$	$-1.31463$
				0				$-1.30860$	$-1.30860$
				0.12				$-1.30008$	$-1.30008$
					0.3			$-1.29782$	$-1.29782$
					0.4			$-1.29404$	$-1.29404$
					0.5			$-1.29029$	$-1.29029$
						$-0.5$		$-1.00143$	$-1.00143$
						0		$-1.301625$	$-1.301625$
						0.5		$-1.70488$	$-1.70488$
							$-0.6$	$0.942750$	$0.942750$
							0	$0.45307$	$0.45307$
							0.6	$-0.48097$	$-0.48097$

and

Table 5. Numerical comparison between present method and bvp4c for Nusselt number against different values of the other parameters.

$\mathit{\beta }$	$\mathit{\gamma }$	$\mathbf{\Lambda }$	$\mathbf{\Gamma }$	$\mathit{\delta }$	$\mathit{\lambda }$	SLLM	bvp4c
0.2	0.2	0.5	1	0.1	0.2	$1.44072$	$1.44072$
1						$1.25957$	$1.25957$
1.5						$1.16052$	$1.16052$
	0.5					$1.36888$	$1.36888$
	0.8					$1.30192$	$1.30192$
	1.2					$1.21886$	$1.21886$
		0.4				$1.42396$	$1.42396$
		1				$1.50962$	$1.50962$
		1.5				$1.56319$	$1.56319$
			0.71			$1.18551$	$1.18551$
			2			$2.15499$	$2.15499$
			3			$2.69120$	$2.69120$
				$-0.1$		$1.55523$	$1.55523$
				0		$1.49965$	$1.49965$
				0.12		$1.42846$	$1.42846$
					0.3	$1.38320$	$1.38320$
					0.4	$1.32612$	$1.32612$
					0.5	$1.26948$	$1.26948$

illustrate the numerical results for the Skin friction and Nusselt number profile as functions of the adjustable parameters. We have also included a comparison with the bvp4c (built-in command in Matlab) in both tables (Table 4 and Table 5). We observe in these Tables that the results obtained by the two methods agree very well, indicating that the current approach, the SLLM, has been appropriately implemented. The data shows that the magnetic parameter $\beta$, the permeability parameter $\gamma$, the wall transpiration parameter $\varsigma$, the shrinking or stretching parameter $\varphi$, and the Prandtl number $\Gamma$ cause the skin friction profile to reduce. On the contrary, the mixed convection parameter $\Lambda$, the heat source or sink parameter $\delta$, and the Eckert number $\lambda$ augment the skin friction profile. The Nusselt number increases significantly because of the effect of the mixed convective parameter $\Lambda$ and the Prandlt number $\Gamma$, whereas the magnetic parameter $\beta$, the permeability parameter $\gamma$, the Eckert number $\lambda$, and the heat source/sink $\delta$ tend to reduce the Nusselt number. Furthermore,

Table 6. Numerical comparison between present results with previously report results [45] for skin friction profile.

$\mathit{\beta }$	$\mathit{\gamma }$	$\mathbf{\Lambda }$	$\mathit{\varsigma }$	$\mathit{\varphi }$	Present Results	Fang et al. [45]
0	0	0	$-\frac{1}{6}$	1	$-1.0234$	$-1.0234$

shows a numerical comparison of the skin friction profile with previously published data for $\beta =0,\gamma =0,\Lambda =0,\varsigma =-1/6,\varphi =1$. It is clear from this table that the current results are compatible with the previous studies.

Figure 1 and Figure 2 provide a graphical comparison of the velocity and temperature profiles. The solid line depicts the results of SLLM, whereas the red circle depicts the results of bvp4c. Matlab software is used to implement both solutions. The bvp4c (built-in command) is based on the shooting method. We can observe from these data that the two graphical outcomes agree well.

Figure 3 demonstrates the impact of the magnetic parameter $\beta$ on the velocity profile. It is evident in this Figure that the velocity profile diminishes as the magnetic parameter increases. Because electromagnetic forces are proportional to magnetic fields, boosting the magnetic parameter increases the electromagnetic force, and this restricts the fluid motion over the whole surface. The temperature profile was significantly improved by increasing the parameter $\beta$, as shown in Figure 4. The Lorentz force opposes the fluid motion and this causes the temperature distribution and thermal boundary layer to increase. Since $\beta =0$ signifies the consequence of non-magnetic mixed convection flow, these findings agree with the case of a linear stretching surface [46].

Figure 5 and Figure 6 show the influence of the permeability parameter $\gamma$ on the velocity and temperature profiles. It is observed that, as the permeability parameter is increased, the velocity profile decreases. In the momentum Equation (8), the permeability parameter includes the drag force component, i.e., $-\gamma f^{\prime }$. Increasing the permeability parameter results in the restriction of the flow. However, as shown in Figure 6, the permeability parameter has the opposite effect on the temperature distribution. The restriction of flow impedes the forced convection and results in the increase of temperature.

Figure 7 demonstrates the effect of the heat generation parameter $\delta$ on the temperature profile and shows that an increase of the heat generation parameter $\delta$ causes an increase in the temperature profile. When $\delta >0$, the heat source effect is to increase the fluid temperature while when $\delta <0$, (heat absorption by a heat sink) the fluid temperature profile decreases.

Figure 8 reveals the impact of the Prandtl number $\Gamma$, which is the ratio of momentum diffusion rate to thermal diffusion rate, on the temperature profile. The temperature profile is seen to decrease as the Prandtl number increases. It is obvious that the momentum diffusion rate is more influential and the increased fluid flow reduces the developed temperature profile.

We can see in Figure 9 that the Grashof number has a significant and positive effect on the velocity profile, because the buoyancy force, $+\Lambda \theta$ in Equation (8), increases as the Grashof number increases. The velocity profile increases because the natural convection velocity intensifies. When $\Lambda \to 0$, natural (free) convection ceases and only forced convection takes place.This coincides to the lowest calculated velocity. On the other contrary, increasing the natural convection parameter $\Lambda$ in Figure 10 reduces the temperature distribution and the thermal boundary layer thickness. However, the impact of this parameter is minor.

Figure 11 depicts the effect of the Eckert number $\lambda$ on the temperature profile. It is apparent that increasing the Eckert number $\lambda$ slightly increases the temperature distribution and the thickness of the thermal boundary layer. Because enhancing the Eckert number actually augments the advective transport and weakens the heat dissipation, the temperature profile tends to rise.

5. Summary and Conclusions

This study focuses on the mixed convective viscous fluid flow across an elastic surface with a porous medium and the role of viscous dissipation. The viscous fluid is electrically conducting, irrotational and incompressible. A similarity analysis was conducted for the momentum and energy equations. The governing equations for the system are numerically solved using a robust numerical method, known as the Spectral Local Linearization Method. The Gauss–Seidel technique was employed to improve the convergence rate of the derived solutions. The machine time required to execute the suggested approach was reported. The main conclusion of the study are shown below:

(i)

The magnetic and permeability parameters augment the temperature profile, while decreasing the velocity profile.

(ii)

The mixed convection parameter increases the velocity profile, while reducing the temperature profile.

(iii)

The Eckert number and the heat generation parameter have a significant effect to increase the temperature profile and the thickness of the thermal boundary layer.

(iv)

Increasing the Prandtl number reduces the temperature profile, because the momentum diffusivity is more dominant and the increased flow dampens the development of the temperature profile.

(v)

The numerical results obtained with the Spectral local linearization method are consistent with the results obtained with other related methodologies such as bvp4c.

(vi)

When the skin friction profile is compared to previously published data, it is found that the current numerical results show a good agreement and this validates the proposed method.

(vii)

The SLLM does not lose accuracy as the number of collocation points increases.

(viii)

For all of the model parameters investigated in this work, the method has been found to rapidly converge to the respective solutions.

(ix)

When compared to other approaches, the suggested methodology is computationally more efficient and shows better performance with fewer collocation points $\left(N\right)$ and four iterations. Because this method is more reliable, simple, and efficient, the SLLM methodology is considered more suitable for the solution of nonlinear boundary value problems.

In general, the present analysis shows the excellent performance of the Spectral Local Linearization Method for the simulation of nonlinear boundary layer problems with porosity and viscous dissipation effects. Since our study is restricted to a Newtonian fluid model, future studies can be extended to non-Newtonian fluids with slip and no-slip effects.