A Simplified Finite Difference Method (SFDM) Solution via Tridiagonal Matrix Algorithm for MHD Radiating Nanofluid Flow over a Slippery Sheet Submerged in a Permeable Medium

In this paper, we turn our attention to the mathematical model to simulate steady, hydromagnetic, and radiating nanoﬂuid ﬂow past an exponentially stretching sheet. A numerical modeling technique, simpliﬁed ﬁnite diﬀerence method (SFDM), has been applied to the ﬂow model that is based on partial diﬀerential equations (PDEs) which is converted to nonlinear ordinary diﬀerential equations (ODEs) by using similarity variables. For the resultant algebraic system, the SFDM uses the tridiagonal matrix algorithm (TDMA) in computing the solution. The eﬀectiveness of numerical scheme is veriﬁed by comparing it with solution from the literature. However, where reference solution is not available, one can compare its numerical results with the results of MATLAB built-in package bvp4c . The velocity, temperature, and concentration proﬁles are graphed for a variety of parameters, i.e., Prandtl number, Grashof number, thermal radiation parameter, Darcy number, Eckert number, Lewis number, and Brownian and thermophoresis parameters. The signiﬁcant eﬀects of the associated emerging thermophysical parameters, i.e., skin friction coeﬃcient, local Nusselt number, and local Sherwood numbers are analyzed and discussed in detail. Numerical results are compared from the available literature and found a close agreement with each other. It is found that the Eckert number upsurges the velocity curve. However, the dimensionless temperature declines with the Grashof number. It is also shown that the SFDM gives good results when compared with the results obtained from bvp4c and results from the literature.


Introduction
e stretching sheet flow has several interesting engineering applications such as in a chemical engineering plant's polymer handling unit and in metallurgy for the metal working system. Crane [1] researched the continuous twodimensional boundary layer flow induced by stretching the sheet moving in its own plane at a velocity linearly varying from a fixed point on the sheet. Immediately after Crane [1], abundant work in this direction is reported and discussed.
Makinde and Aziz [2] explored the effect of boundary layer flow over linearly stretching nanofluid while Mustafa et al. [3] concentrated on boundary layer flow for an exponentially stretching sheet and solved the issue using the technique of the homotopy analysis method to calculate analytical solutions. Realistically, as discussed by Gupta and Gupta [4], stretching a plastic sheet may not necessarily be linear.
Since Choi and Eastman's pioneering research [5], surveys associated to nanofluid dynamics have risen significantly in contemporary times due to the low thermal conductivity of prevalent heat transfer liquids, which causes the device to function inefficiently and consume additional energy. A new method has been introduced to optimize machine operation by dispersing solid particles with a base fluid. Nanofluid defines the suspension in standard base liquids such as water, ethylene glycol, and motor oil of strong particles of a nanometer size. References from [2,3,5] give a thorough overview of the nanofluid literature. Innovative fluid types are needed in these days to achieve more effective output. Sheikholeslami and Bhatti [6] studied forced convection of nanofluid considering nanoparticles' shape impacts.
e Brownian motion impact on nanofluid flow within a porous cavity was recently regarded by Sheikholeslami [7]. He found that convective flow enhances with increase in Darcy number. e research of electrically conductive fluid flow has many applications in engineering issues such as MHD generators, plasma research, nuclear reactors, geothermal energy extraction, and aerodynamic boundary layer control [8]. Reddy et al. [9] studied the effects of frictional and irregular temperature on non-Newtonian MHD fluid flows owing to stretched surface.
Mishra and Singh [10] addressed dual solutions of mixed convection flow with momentum and heat slip over a permeable shrinking cylinder. When studying the flow models in nanoscales or microscales, the interaction of the fluid surface is mostly controlled by models of slip flow. ese models were checked from asymptotic solution using the Boltzmann technique where the internal kinetic solution matches the outer (i.e., bulk) Navier-Stokes solution, and matching is only achieved when the slip/jump coefficient is regarded at the border or surface (Hadjiconstantinou [11,12]). erefore, the slip coefficients are the result of these assessments. Due to its simplicity, the slip flow phenomenon is always preferred to no-slip situations. e Navier-Stokes equations are still valid here, and only the boundary conditions change in compliance with the slip flow model. With the newly suggested second-order slip flow model, Fang et al. [13] evaluated the slip flow over a permeable shrinking surface. Ullah et al. [14] examined the two-dimensional flow of Reiner-Philippoff fluid thin films over an unstable stretching sheet in the variable heat distribution and radiation.
Khan et al. [15] provided thermal radiation and viscous dissipation impacts on the unstable nanofluid boundary layer flow over a stretching sheet. In this research, they accounted for the viscous dissipation impact and discovered that the heat boundary layer thickness is increased by increasing the values of Eckert number. Ibrahim and Shankar [16] evaluated the impact of thermophoresis on Brownian fluid movement owing to stretching sheet.
Several technological systems depend on the impact of buoyancy. Makinde et al. [17] examined combined impacts of buoyancy force, convective warming, Brownian movement, thermophoresis, and magnetic field on stagnation point stream and heat exchange due to nanofluid stream towards an extending sheet. Ali and Yousef [18] analyzed laminar mixed convection heat transfer from continuously stretching vertical surface with energy functional form for wall temperature by considering the impact of buoyancy. Mixed convection heat transfer from an exponentially stretching sheet was explored by Partha et al. [19]. ey also analyzed influence of buoyancy along with viscous dissipation, and the flow is governed by the mixed convection parameter (Gr/Re 2 ). e effect of viscous dissipation in natural convection process has been investigated by Gebhart [20] and Gebhart and Mollendorf [21]. Magyari and Keller [22] analytically as well as numerically evaluated the continuous free fluid flow and thermal transfer from an exponentially stretching vertical surface with an exponential temperature distribution. Unsteady flow of thermally radiating nanofluid over nonlinearly stretching sheet was discussed by Seth et al. [23]. ey noted that the nanofluid's velocity curve depends on the unsteadiness, velocity slip, and stretching velocity nonlinearity. Makinde et al. [24] reported the two-dimensional unsteady MHD radiating electrically conducting fluid past a slippery stretching sheet embedded in a porous medium. Using the explicit finite difference scheme, they solved the system of higher-order nonlinear PDEs. Hamid and Khan [25] have discussed the thermophysical properties of the flow of Williamson nanofluid and solved their problem numerically. ey concluded that the stronger the magnetic field resulted in decreasing of boundary layer thickness. Some other references in this direction can be consulted in [26][27][28].
Qing et al. [29] researched the entropy generation of nanofluid owing to a magnetic field over a stretching surface. Hosseini et al. [30] discussed heat transfer of nanofluid flow in microchannel heat sink (MCHS) in the presence of a magnetic field. e influence of chemical reaction and heat generation/ absorption on mixed convective flow of nanofluid past an exponentially stretched surface has been examined by Eid [31], and numerical solutions have been obtained by utilizing the shooting technique along with the Runge-Kutta-Fehlberg method. Afify and Elgazery [32] investigated numerically the boundary layer flow of Maxwell nanofluid with convective boundary condition and heat absorption. e result showed that nanoparticle concentration reduces with higher chemical reaction parameter whereas a reverse pattern is noted for temperature. Reviews of viscous fluid flow problems for nonlinear stretching sheet have been presented by Prasad et al. [33], Afzal [34], and Nandeppanavar et al. [35]. Nadeem and Lee [36] studied analytically the problem of steady boundary layer flow of nanofluid over an exponentially stretching surface including the effects of Brownian motion and thermophoresis parameters. e influence of solar energy radiations in the time-dependent Hiemenz flow of nanofluid over a wedge was discussed by Mohamad et al. [37]. In [38], Sheikholeslami et al. discussed natural convection inside a sinusoidal annulus. Tripathi et al. [39] reported shape effects of nanoparticle on blood flow in a microvascular flow. Bhatti et al. [40] have discussed the movement of gyrotactic microorganism in a magnetized nanofluid over a plate. Ibrahim and Anbessa [41] discussed Casson nanofluid with Hall and Ion slip effects. Ibrahim and Negera [42] investigated Williamson nanofluid over a stretching cylinder with activation energy. For similar work in this direction, the reader referred to [43].
In all previous studies, a usual course is followed in one way or the other and discussion is intended towards linearly or nonlinearly stretching sheets in the absence of some important emerging parameters. e aim of this work is to add numerical methodology, SFDM, in the literature so that it can be applicable to many problems containing coupled ODEs. To the best of our knowledge, the current mathematical model along with numerical consideration has not been discussed before. e paper is planned in the following order. Section 2 commences by laying out the mathematical model of the physical problem. Numerical procedure is opted and discussed in Section 3. In the same section, the detailed description of the SFDM is given. As a consequence of numerical calculations, results and discussion are followed in Section 4. At the end of the paper, the conclusions are presented in Section 5.

Mathematical Formulation
We deliberate a two dimensional, steady, incompressible, laminar, and MHD flow of an electrically conducting nanofluid occupied over a slippery stretching sheet submerged in a porous medium. e geometrical description of fluid flow over a sheet is shown in Figure 1. In the figure, x-axis has been chosen along the sheet and y-axis normal to it.
After making use of these assumptions, the set of continuity, momentum, energy, and concentration equations incorporating the Buongiorno model is written as follows [44]: here the velocity components (u, v) are considered along and normal of the sheet. μ is the coefficient of viscosity, ρ is the density of the fluid, σ is the electrical conductivity of the fluid, T is fluid's temperature, K is the permeability, β is the thermal expansion coefficient, k is the thermal conductivity, C p is the specific heat capacity at constant pressure, q r is the radioactive heat flux, Q is the heat source coefficient, C is the concentration, and τ � (ρC) p /(ρC) f , where (ρC) p and (ρC) f are heat capacities of the nanofluid and base fluid, respectively. Also, D B and D T are Brownian and thermophoretic diffusion coefficients, respectively. T ∞ is the ambient fluid temperature, and C ∞ is the ambient fluid concentration.

Boundary
Conditions. e preceding mathematical model allows the following boundary condition: where L 1 is the slip length. Here, U w � U 0 e x/L is the stretching velocity, where U 0 is the reference velocity. And, T w � T ∞ + T 0 e x/(2L) is the variable temperature at the sheet with T 0 being a reference temperature. Also, C w � C ∞ + C 0 e x/(2L) is the variable concentration at the sheet with C 0 being a constant.

Method of Solution.
By introducing similarity variables η, ψ(η), θ(η), and ϕ(η) as dimensionless independent variable, stream function, temperature and concentration for the momentum, energy and concentration equations (1)-(4), and in the boundary conditions (5), gives the following nonlinear ordinary differential equations: Mathematical Problems in Engineering In the above equations, various parameters appear which are M, Pr, Nt, Nb, Gr, Da, Nr, Ec, S, Le, and λ. In order, these are the magnetic parameter, Prandtl number, thermophoresis parameter, Brownian parameter, Grashof number, Darcy's number, thermal radiation effect, heat source or sink, Lewis number, and the slip parameter. eir expressions are grouped as follows:

Physical
Quantities. Now that the flow equations are known, the physical quantities that measures roughness, heat transfer rate, and concentration rate at the sheet can be obtained. First, the skin friction coefficient C f is given by Second, the local Nusselt number Nu x is written as ird, the local Sherwood number Sh x is defined as After substituting similarity variables in (13)- (15), this yields the expressions as follows: Here, Re x � Ux/] is a local Reynolds number.

Numerical Procedures
In search of solution for the above problem given in equations (7)- (9), the only plausible way to compute solution is numerically. We find numerical solutions by using two numerical techniques. e first numerical method we use is the SFDM, and the second one is the famous algorithm written in MATLAB and commonly known as bvp4c. us, due details on the SFDM will be presented first followed by brief description on bvp4c.

SFDM.
is work is influenced by Na [45] in which some numerical results are displayed for linear ODEs. For coupled nonlinear ODEs, we expand these ideas theoretically and execute them in MATLAB. e algorithm with necessary details for the SFDM is as follows: (1) Reduce third-order ODEs to a pair of ODEs of the first and second order (2) Use Taylor series to linearize the system of nonlinear ODEs (3) Substitute finite difference formulas in the derivatives (4) Finally, solve the algebraic system by TDMA e results are shown for N � 1000 grid points. Generally, the domain length varies with different parameters. However, the domain value η � 7 seems enough to show steady state results. To initiate the SFDM procedure, we assume f′ � F in equation (7) and we get Define a new variable as (18) and approximate (dF/dη) by forward difference approximation with constant width h e coefficients are written as Simplifying the above, we reach at where In the matrix-vector form, it is written in compact as where e matrix A is a tridiagonal matrix and is written in LU factorization as where Mathematical Problems in Engineering where L and U are the lower and upper triangular matrices, respectively. Here, the unknowns (β i , c i ), i � 1, 2, . . . , N − 1 are to be related as After defining these relations, (23) becomes and we have e unknown elements of s are written as We get which is a solution of (17). We can easily find f from f ′ � F which is in the discretization form written as follows: which gives a required solution of (7). A similar procedure can also be opted for solutions θ and ϕ. For the sake of brevity, we only present coefficients for these ODEs and leave the details which follows on the same line as presented above. For example, we have the energy and concentration equation as follows: 6 Mathematical Problems in Engineering Similarly, the coefficients for (9) are written as Boundary conditions are discretized as

bvp4c.
is section presents the second numerical method of the studied problem given in (7)-(9) which subject to the boundary conditions (10) and (11). We use MATLAB built-in function bvp4c for this purpose. For description and details of this method, one can refer to [46]. Let us define the variables as Mathematical Problems in Engineering e system of first-order equations is given as follows: and boundary conditions are given as follows:

Results and Discussion
In this section, the focus is to analyze the role of embedded parameters on the velocity, temperature, and concentration. Results of the current study are displayed in the tabular as well as graphical form. For − f ″ (0), the results are compared with the solutions published in the literature, and this comparison is listed in Table 1. e results demonstrate that the numerical values of SFDM are accurate and closely agreed with one another.
In Table 2, when the admissible values of the magnetic parameter increase resultantly, the skin friction coefficient also increases. However, reduction in both temperature and concentration gradients is observed. One can also observe that the magnitude of the local Nusselt and the local Sherwood numbers increases and the skin friction coefficient decreases with the rise of values of Darcy's number. Grashof number enhances local Nusselt number and local Sherwood number whereas this reduces the skin friction coefficient. However, Lewis number Le causes slight change in skin friction coefficient while concentration gradient and wall temperature gradient reduce.
It is also evident from Table 2 that local Sherwood number increases by increasing Nt, but the effect is seen to be reverse on skin friction coefficient while local Nusselt number remains constant. e skin friction coefficient, local Nusselt number, and local Sherwood number decrease with respect to thermophoretic parameter Nb (see Table 3). From Table 3, one can observe an increase in local Sherwood number along the range of Ec. is also causes a surge in local Sherwood number whereas its effect on the skin friction coefficient and local Nusselt number is opposite. Figure 2 shows a decreasing trend in velocity profiles against M(2 ≤ M ≤ 8)to the point where η ≈ 2.50. After this point, the boundary layer thickness demonstrates the opposite behaviour. Figure 3 illustrates an increase in thermal boundary layer thickness due to an increase in a magnetic parameter. However, minor increase in concentration profile is presented in Figure 4. e reduction of the momentum boundary layer is strongly influenced by the magnetic parameter strength which produces Lorentz force and that offer resistance to the flow.

Effects of Darcy Number Da.
Darcy number Da characterizes the strength of permeability of the porous medium. Figures 5 and 6 depict increasing values of Darcy number (0.5 ≤ Da ≤ 15.5) that increase the velocity profile while concentration profile decreases. However, the temperature decreases in the boundary layer region. us, the thickness of the thermal boundary layer decreases as shown in Figure 7. Figure 8 displays the variations of velocity profiles due to the variations in the values of Lewis number (0.5 ≤ Le ≤ 2). It is observed that the velocity profile decreases with an increase in Le. For Le ≥ 1, the mass transport is dominant that resists the flow. One can also observe that in Figure 9, the temperature profile as well as the thickness of the boundary layer initially increases and then decreases with an increase due to Lewis number. is implies that the momentum boundary layer thickness decreases when a ratio of thermal diffusivity to a mass diffusivity increases. Figure 10 demonstrates the nanoparticle volume fraction for several values of Lewis number Le accompanying reduction in concentration boundary layer thickness.

Effects of Grashof Number Gr.
e Gr approximates the ratio of buoyancy to viscous forces and represents how dominant is buoyancy force which is responsible for the conviction comparing to viscous forces. Either convection or viscous forces are dominant, and the results are displayed in Figures 11-13. It can be observed that temperature and concentration decrease with the Grashof number Gr(3 ≤ Gr ≤ 6.5), but there is an abrupt change in a velocity profile.  bvp4c   bvp4c  Figure 14 shows that the velocity profile increases with Nt in the range (0.5 ≤ Nt ≤ 2). Figure 15 illustrates the variations of thermophoretic parameter on temperature profile. It validated the fact that thermophoretic parameter enhances the temperature profile. Since the thermophoretic phenomenon transferred nanoparticles from hot surface to the cold region, it resulted in increasing the temperature of the fluid. Figure 16 suggests that a stronger thermophoretic parameter produces minor change in nanoparticle volume fraction.

Effects of Nb. Incremental Brownian parameter
Nb (5 ≤ Nb ≤ 20) causes slight change in nanoparticles volume fraction which increases the velocity profile as presented in Figures 17 and 18. Figure 19 suggests that a stronger Brownian motion is responsible for an increase in thermal boundary layer thickness.  temperature, and concentration profile have been shown. It depicts that, in the absence of Ec, the dimensionless velocity is lowest at the surface and then increases with increasing Ec. e dimensionless temperature is lowest inside the thermal boundary layer and increases with Ec while the effect of a viscous dissipation is insignificant on concentration profile.

Conclusions
is study focuses on two objectives as follows: (1) to produce a mathematical model for nanofluid flow considering MHD, radiation, porosity, and slippery exponentially stretching sheet which is immersed in a porous medium and (2) to develop a new numerical scheme for the solution of a general nonlinear ODEs which is applicable not only for the current problem but many others. e key observations of present work are as follows: (i) e skin friction coefficient grows with a rise of the values of magnetic parameter but local Nusselt and Sherwood numbers reduces e rise of magnetic parameter causes decrease in velocity initially after that the reverse effects can be seen. It enhances temperature and nanoparticle concentration of the boundary layer regime. (ii) e Darcy number causes the thermal boundary layer and solute concentration to reduce whereas it enhances the momentum boundary layer. e skin friction coefficient decreases yet the wall temperature gradient and nanoparticle concentration increase with an increase in Darcy number. (iii) e Grashof number enhances local Nusselt number and local Sherwood number whereas this reduces the skin friction coefficient. (iv) e Lewis number increases momentum and concentration boundary layers. Although it causes a slight change in the skin friction coefficient, wall temperature gradient and boosts nanoparticle concentration are reduced.
(v) e thermophoresis parameter causes both the thermal and momentum boundary layer to increase while its effect on nanoparticle volume fraction is insignificant. e skin friction coefficient and local Sherwood number increase while local Nusselt number remains unaltered. (vi) e Brownian parameter increases thickness of the thermal boundary layer and momentum boundary layer while there is a slight change in the concentration boundary layer. It causes a reduction in wall temperature gradient, the skin friction coefficient, the local Nusselt number, and the local Sherwood number. (vii) e Eckert number Ec reduces the skin friction coefficient, local Nusselt number, and concentration boundary layer while momentum boundary layer, thermal boundary layer, and local Sherwood number increase. (viii) e SFDM has been successfully developed and applied in the current problem. One can show that the SFDM is easy to implement and converges quickly. (ix) To validate, one compares the SFDM results with bvp4c and those with the literature which gives a good account of agreement with each other.

Data Availability
ere are no experiment data involved in this study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.