# 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

Authors: M. Asif Farooq; A. Salahuddin; Asif Mushtaq; M. Razzaq
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, we turn our attention to the mathematical model to simulate steady, hydromagnetic, and radiating nanofluid flow past an exponentially stretching sheet. A numerical modeling technique, simplified finite difference method (SFDM), has been applied to the flow model that is based on partial differential equations (PDEs) which is converted to nonlinear ordinary differential equations (ODEs) by using similarity variables. For the resultant algebraic system, the SFDM uses the tridiagonal matrix algorithm (TDMA) in computing the solution. The effectiveness of numerical scheme is verified 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 packagebvp4c. The velocity, temperature, and concentration profiles 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 significant effects of the associated emerging thermophysical parameters, i.e., skin friction coefficient, 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.


---

## Body

## 1. Introduction

The 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 two-dimensional 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. The 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.The 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. These 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]). Therefore, the slip coefficients are the result of these assessments. Due to its simplicity, the slip flow phenomenon is always preferred to no-slip situations. The 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]. They also analyzed influence of buoyancy along with viscous dissipation, and the flow is governed by the mixed convection parameter Gr/Re2. The 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]. They 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. They 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–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. The 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. The 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. The 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. The 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.The paper is planned in the following order. Section2 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.


## 2. 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. The geometrical description of fluid flow over a sheet is shown in Figure1. In the figure, x-axis has been chosen along the sheet and y-axis normal to it.Figure 1
Schematic diagram of the problem.After making use of these assumptions, the set of continuity, momentum, energy, and concentration equations incorporating the Buongiorno model is written as follows [44]:(1)∂xu+∂yv=0,(2)uux+vuy=νuyy−σBo2uρ−νuK+gβT−T∞,(3)uTx+vTy=kρCpTyy+νCpuy2+σBo2u2ρCp+νu2CpK−1ρCpqr,y+QT−T∞ρCp+τDBCyTy+DTT∞Ty2,(4)uCx+vCy=DBCyy+DTT∞Tyy,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, Cp is the specific heat capacity at constant pressure, qr is the radioactive heat flux, Q is the heat source coefficient, C is the concentration, and τ=ρCp/ρCf, where ρCp and ρCf are heat capacities of the nanofluid and base fluid, respectively. Also, DB and DT are Brownian and thermophoretic diffusion coefficients, respectively. T∞ is the ambient fluid temperature, and C∞ is the ambient fluid concentration.


### 2.1. Boundary Conditions

The preceding mathematical model allows the following boundary condition:(5)ux,0=Uw+μL1uy,vx,0=0,Tx,0=Tw,Cx,0=Cw,u⟶0,T⟶T∞,C⟶C∞,asy⟶∞,where L1 is the slip length. Here, Uw=U0ex/L is the stretching velocity, where U0 is the reference velocity. And, Tw=T∞+T0ex/2L is the variable temperature at the sheet with T0 being a reference temperature. Also, Cw=C∞+C0ex/2L is the variable concentration at the sheet with C0 being a constant.


### 2.2. 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),(6)η=U02νLex/2Ly,ψη=2U0νLex/2Lfη,u=U0f′ex/L,v=−U0ν2Lex/2Lf′η+f,θη=T−T∞T0e−x/2L,ϕη=C−C∞C0e−x/2L,gives the following nonlinear ordinary differential equations:(7)f‴+ff″−2f′2−M+1Daf′+2Grθ=0,(8)Prθf′−fθ′−1+Nrθ″−PrEcf″2−PrEcM+1Daf′2−2PrSθ−PrNbϕ′θ′−PrNtθ′2=0,(9)ϕ″+NtNbθ″+Lefϕ′−f′ϕ=0,(10)fη=0,f′η=1+λf″η,θη=1,ϕη=1,asη⟶0,(11)f′η⟶0,θη⟶0,ϕη⟶0,asη⟶∞.In the above equations, various parameters appear which areM, 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. Their expressions are grouped as follows:(12)Pr=μCpk,M=2LσB02Uwρ,Nt=τDTTw−T∞T∞ν,Nb=τDBCw−C∞ν,Le=νDB,Gr=βgLTw−T∞Uw2,Da=KUw2νL,Ec=Uw2TwCp,S=QLUwρCp.


### 2.3. 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 coefficientCf is given by(13)Cf=νUw2∂u∂yy=0.Second, the local Nusselt numberNux is written as(14)Nux=−1+NrxTw−T∞∂T∂yy=0.Third, the local Sherwood numberShx is defined as(15)Shx=−xCw−C∞∂C∂yy=0.After substituting similarity variables in (13)–(15), this yields the expressions as follows:(16)Cf=12Rexf″0,Nux=−1+NrxRex2Lθ′0,Shx=−xRex2Lϕ′0.Here,Rex=Ux/ν is a local Reynolds number.


## 2.1. Boundary Conditions

The preceding mathematical model allows the following boundary condition:(5)ux,0=Uw+μL1uy,vx,0=0,Tx,0=Tw,Cx,0=Cw,u⟶0,T⟶T∞,C⟶C∞,asy⟶∞,where L1 is the slip length. Here, Uw=U0ex/L is the stretching velocity, where U0 is the reference velocity. And, Tw=T∞+T0ex/2L is the variable temperature at the sheet with T0 being a reference temperature. Also, Cw=C∞+C0ex/2L is the variable concentration at the sheet with C0 being a constant.


## 2.2. 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),(6)η=U02νLex/2Ly,ψη=2U0νLex/2Lfη,u=U0f′ex/L,v=−U0ν2Lex/2Lf′η+f,θη=T−T∞T0e−x/2L,ϕη=C−C∞C0e−x/2L,gives the following nonlinear ordinary differential equations:(7)f‴+ff″−2f′2−M+1Daf′+2Grθ=0,(8)Prθf′−fθ′−1+Nrθ″−PrEcf″2−PrEcM+1Daf′2−2PrSθ−PrNbϕ′θ′−PrNtθ′2=0,(9)ϕ″+NtNbθ″+Lefϕ′−f′ϕ=0,(10)fη=0,f′η=1+λf″η,θη=1,ϕη=1,asη⟶0,(11)f′η⟶0,θη⟶0,ϕη⟶0,asη⟶∞.In the above equations, various parameters appear which areM, 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. Their expressions are grouped as follows:(12)Pr=μCpk,M=2LσB02Uwρ,Nt=τDTTw−T∞T∞ν,Nb=τDBCw−C∞ν,Le=νDB,Gr=βgLTw−T∞Uw2,Da=KUw2νL,Ec=Uw2TwCp,S=QLUwρCp.


## 2.3. 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 coefficientCf is given by(13)Cf=νUw2∂u∂yy=0.Second, the local Nusselt numberNux is written as(14)Nux=−1+NrxTw−T∞∂T∂yy=0.Third, the local Sherwood numberShx is defined as(15)Shx=−xCw−C∞∂C∂yy=0.After substituting similarity variables in (13)–(15), this yields the expressions as follows:(16)Cf=12Rexf″0,Nux=−1+NrxRex2Lθ′0,Shx=−xRex2Lϕ′0.Here,Rex=Ux/ν is a local Reynolds number.


## 3. 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. The first numerical method we use is the SFDM, and the second one is the famous algorithm written in MATLAB and commonly known as bvp4c. Thus, due details on the SFDM will be presented first followed by brief description on bvp4c.


### 3.1. SFDM

This 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. The 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 TDMAThe results are shown forN=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(17)d2F¯dη2=−fdF¯dη+2F¯2+M+1DaF¯−2Grθ,Define a new variable as(18)ξ1η,F¯,F¯′=−fdF¯dη+2F¯2+M+1DaF¯−2Grθand approximate dF¯/dη by forward difference approximation with constant width h(19)ξ1η,F¯,F¯′=−fiF¯i+1−F¯ih+2F¯i2+M+1DaF¯i−2Grθi.The coefficients are written as(20)A¯n=−∂f∂F¯′=−−f=f=fi,B¯n=−∂f∂F¯=−4F¯−M+1Da,B¯n=−4F¯i−M+1Da,D¯n=ξ1η,F¯,F¯′+B¯nF¯i+A¯nF¯i+1−F¯ih.Simplifying the above, we reach at(21)a¯iF¯i−1+b¯iF¯i+c¯iF¯i+1=r¯i,i=1,2,3,…,N,where(22)a¯i=2−hA¯n,b¯i=2h2B¯n−4,c¯i=2+hA¯n,r¯i=2h2D¯n.In the matrix-vector form, it is written in compact as(23)A¯F¯=s¯,where(24)A¯=b¯1c¯1a¯2b¯2c¯2….a¯N−2b¯N−2c¯N−2a¯N−1b¯N−1,F¯=F¯1F¯2⋅⋅F¯N−1,s¯=s¯1s¯2⋅⋅s¯N−1.The matrixA¯ is a tridiagonal matrix and is written in LU factorization as(25)A¯=L¯U¯,where(26)L¯=β¯1a¯2β¯2….a¯N−2β¯N−2a¯N−1β¯N−1,U¯=1γ¯11γ¯2….1γ¯N−21,where L¯ and U¯ are the lower and upper triangular matrices, respectively. Here, the unknowns β¯i,γ¯i,i=1,2,…,N−1 are to be related as(27)β¯1=−1−λ¯h,γ¯1=λ¯β¯1h,β¯i=b¯i−a¯iγ¯i−1,i=2,3,…,N−1,β¯iγ¯i=c¯i,i=2,3,…,N−2.After defining these relations, (23) becomes(28)L¯U¯F¯=s¯,U¯F¯=z¯,L¯z¯=s¯,and we have(29)β¯1a¯2β¯2….a¯N−2β¯N−2a¯N−1z¯1z¯2z¯3⋅⋅⋅z¯N−2z¯N−1=s¯1s¯2s¯3⋅⋅⋅s¯N−2s¯N−1.The unknown elements ofs¯ are written as(30)z¯1=s¯1β¯1,s¯i=s¯i−a¯iz¯i−1β¯i,i=2,3,…,N−1,1γ¯11γ¯2….1γ¯N−21F¯1F¯2⋅⋅⋅F¯N−2F¯N−1=z¯1z¯2⋅⋅⋅z¯N−2z¯N−1.We get(31)F¯i−1=z¯i−1,F¯i=z¯i−γ¯iF¯i+1,i=N−2,N−3,…,3,2,1,which is a solution of (17). We can easily find f from f′=F¯ which is in the discretization form written as follows:(32)fi+1−fih=F¯i,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:(33)d2θdη2=11+NrPrθF¯−fdθdη−PrEcdF¯dη2−PrEcM+1DaF¯2−2PrSθ−PrNbdϕdηdθdη−PrNtdθdη2,ξ2η,θ,θ′=11+NrPrθiF¯i−fiθi−θi−1h−PrEcF¯i−F¯i−1dη2−PrEcM+1DaF¯i2−2PrSθi−11+NrPrNbϕi−ϕi−1hθi−θi−1h−PrNtθi−θi−1dη2,Ann=−∂f∂θ′=−11+Nr−Prf−PrNbdϕdη−2PrNtdθdη,Ann=−∂f∂θ′=11+NrPrfi+PrNbdϕi−ϕi−1h+2PrNtdθi−θi−1h,Bnn=−∂f∂θ=−11+NrPrF¯−2PrS,Bnn=−∂f∂θ=−11+NrPrF¯i−2PrS,d2ϕdη2=−NtNbd2θdη2−Lefdϕdη−F¯ϕ,ξ3η,ϕ,ϕ′=−NtNbθi−1−2θi+θi+1dη2−Lefiϕi−ϕi−1h−F¯iϕi.Similarly, the coefficients for (9) are written as(34)Annn=Lefi,Bnnn=−LeF¯i.Boundary conditions are discretized as(35)F¯1=1+λF¯2−F¯1h.


### 3.2.bvp4c

This 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(36)y1=f,y2=f′,y3=f″,y4=θ,y5=θ′,y6=ϕ,y7=ϕ′.The system of first-order equations is given as follows:(37)y1′=f′=y2,y2′=f″=y3,y3′=f‴=−y1y3+2y22+M+1Day2−2Gry4,y4′=θ′=y5,y5′=θ″=11+NrPry4y2−Pry1y5−PrEcy32−M+1DaEcPry22−2PrSy4−PrNby7y5−NtPry52,y6′=ϕ′=y7,y7′=ϕ″=Ley6y2−Ley1y7−NtNby5′,and boundary conditions are given as follows:(38)y01=0,y02−1−λy03=0,y04−1=0,y06−1=0,yinf2=0,yinf4=0,yinf6=0.


## 3.1. SFDM

This 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. The 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 TDMAThe results are shown forN=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(17)d2F¯dη2=−fdF¯dη+2F¯2+M+1DaF¯−2Grθ,Define a new variable as(18)ξ1η,F¯,F¯′=−fdF¯dη+2F¯2+M+1DaF¯−2Grθand approximate dF¯/dη by forward difference approximation with constant width h(19)ξ1η,F¯,F¯′=−fiF¯i+1−F¯ih+2F¯i2+M+1DaF¯i−2Grθi.The coefficients are written as(20)A¯n=−∂f∂F¯′=−−f=f=fi,B¯n=−∂f∂F¯=−4F¯−M+1Da,B¯n=−4F¯i−M+1Da,D¯n=ξ1η,F¯,F¯′+B¯nF¯i+A¯nF¯i+1−F¯ih.Simplifying the above, we reach at(21)a¯iF¯i−1+b¯iF¯i+c¯iF¯i+1=r¯i,i=1,2,3,…,N,where(22)a¯i=2−hA¯n,b¯i=2h2B¯n−4,c¯i=2+hA¯n,r¯i=2h2D¯n.In the matrix-vector form, it is written in compact as(23)A¯F¯=s¯,where(24)A¯=b¯1c¯1a¯2b¯2c¯2….a¯N−2b¯N−2c¯N−2a¯N−1b¯N−1,F¯=F¯1F¯2⋅⋅F¯N−1,s¯=s¯1s¯2⋅⋅s¯N−1.The matrixA¯ is a tridiagonal matrix and is written in LU factorization as(25)A¯=L¯U¯,where(26)L¯=β¯1a¯2β¯2….a¯N−2β¯N−2a¯N−1β¯N−1,U¯=1γ¯11γ¯2….1γ¯N−21,where L¯ and U¯ are the lower and upper triangular matrices, respectively. Here, the unknowns β¯i,γ¯i,i=1,2,…,N−1 are to be related as(27)β¯1=−1−λ¯h,γ¯1=λ¯β¯1h,β¯i=b¯i−a¯iγ¯i−1,i=2,3,…,N−1,β¯iγ¯i=c¯i,i=2,3,…,N−2.After defining these relations, (23) becomes(28)L¯U¯F¯=s¯,U¯F¯=z¯,L¯z¯=s¯,and we have(29)β¯1a¯2β¯2….a¯N−2β¯N−2a¯N−1z¯1z¯2z¯3⋅⋅⋅z¯N−2z¯N−1=s¯1s¯2s¯3⋅⋅⋅s¯N−2s¯N−1.The unknown elements ofs¯ are written as(30)z¯1=s¯1β¯1,s¯i=s¯i−a¯iz¯i−1β¯i,i=2,3,…,N−1,1γ¯11γ¯2….1γ¯N−21F¯1F¯2⋅⋅⋅F¯N−2F¯N−1=z¯1z¯2⋅⋅⋅z¯N−2z¯N−1.We get(31)F¯i−1=z¯i−1,F¯i=z¯i−γ¯iF¯i+1,i=N−2,N−3,…,3,2,1,which is a solution of (17). We can easily find f from f′=F¯ which is in the discretization form written as follows:(32)fi+1−fih=F¯i,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:(33)d2θdη2=11+NrPrθF¯−fdθdη−PrEcdF¯dη2−PrEcM+1DaF¯2−2PrSθ−PrNbdϕdηdθdη−PrNtdθdη2,ξ2η,θ,θ′=11+NrPrθiF¯i−fiθi−θi−1h−PrEcF¯i−F¯i−1dη2−PrEcM+1DaF¯i2−2PrSθi−11+NrPrNbϕi−ϕi−1hθi−θi−1h−PrNtθi−θi−1dη2,Ann=−∂f∂θ′=−11+Nr−Prf−PrNbdϕdη−2PrNtdθdη,Ann=−∂f∂θ′=11+NrPrfi+PrNbdϕi−ϕi−1h+2PrNtdθi−θi−1h,Bnn=−∂f∂θ=−11+NrPrF¯−2PrS,Bnn=−∂f∂θ=−11+NrPrF¯i−2PrS,d2ϕdη2=−NtNbd2θdη2−Lefdϕdη−F¯ϕ,ξ3η,ϕ,ϕ′=−NtNbθi−1−2θi+θi+1dη2−Lefiϕi−ϕi−1h−F¯iϕi.Similarly, the coefficients for (9) are written as(34)Annn=Lefi,Bnnn=−LeF¯i.Boundary conditions are discretized as(35)F¯1=1+λF¯2−F¯1h.


## 3.2.bvp4c

This 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(36)y1=f,y2=f′,y3=f″,y4=θ,y5=θ′,y6=ϕ,y7=ϕ′.The system of first-order equations is given as follows:(37)y1′=f′=y2,y2′=f″=y3,y3′=f‴=−y1y3+2y22+M+1Day2−2Gry4,y4′=θ′=y5,y5′=θ″=11+NrPry4y2−Pry1y5−PrEcy32−M+1DaEcPry22−2PrSy4−PrNby7y5−NtPry52,y6′=ϕ′=y7,y7′=ϕ″=Ley6y2−Ley1y7−NtNby5′,and boundary conditions are given as follows:(38)y01=0,y02−1−λy03=0,y04−1=0,y06−1=0,yinf2=0,yinf4=0,yinf6=0.


## 4. 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. The results demonstrate that the numerical values of SFDM are accurate and closely agreed with one another.Table 1
The comparison of skin friction coefficient to previous data forλ=S=Ec=Nr=Gr=0 and Da=∞ and for various values Pr, M, Nb, Nt, and Le.
Present resultPrMNbNtLeSharif et al. [47]bvp4cSFDM0.700.50.511.281831.28180891.2646694—0.1———1.321041.32101481.3030810—0.2———1.358951.35895751.3402296—0.3———1.395811.39577451.3762525In Table2, 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.Table 2
Results for−f″0, −θ′0, and −ϕ′0 obtained by fixing values of parameters Pr=6.2, Nt=2, Ec=0.2, Nb=8, S=0.1, Nr=5, and λ=3.
bvp4cSFDMMDaGrLe−f″0−θ′0−ϕ′0−f″0−θ′0−ϕ′0140.380.19450.01002.20700.19451750.00992192.1967401.1———0.19770.00942.17690.19770450.00934812.16699091.2———0.20080.00882.14770.20075730.00878792.13809161.3———0.20370.00832.11930.20368180.00824402.1100172150.380.19290.01032.22230.19287130.01022302.2119414—6——0.19180.01052.23270.19175380.0104232.2221985—7——0.19100.01062.24010.19094580.01056742.2295847—8——0.19030.01072.24580.19033400.01067632.2351598150.480.17540.01212.38050.17532670.01210912.3689312−2—0.5—0.15920.01362.51610.15915680.01360772.5033457——0.6—0.14410.01492.63540.14405710.01484172.6216478——0.7—0.12990.01592.74260.12982390.01588282.7277798150.790.13020.01592.90790.13017680.01496052.8912397———130.13130.01243.49000.13129320.01241133.4660959———170.13210.01083.98630.13210690.01083053.9552333———210.13280.00974.42630.13273880.00972924.3880280It is also evident from Table2 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. The 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. This also causes a surge in local Sherwood number whereas its effect on the skin friction coefficient and local Nusselt number is opposite.Table 3
Results for−f″0, −θ′0, and −ϕ′0 obtained by various values of parameters Pr=6.2, M=2, Da=5, Gr=0.7,Nr=5, Le=8, and λ=3.
bvp4cSimplified FDMNbNtEc−f″0−θ′0−ϕ′0−f″0−θ′0−ϕ′0220.20.16600.02702.48760.16821580.07464422.49738306——0.16480.01482.48000.16475640.01482022.46796918——0.16410.01042.47770.16402380.01040382.465646110——0.16360.00812.47740.16352270.00813782.4652570100.50.10.16560.01062.44770.16559620.01059042.4350691—1—0.16500.01062.45800.16494190.01060632.4455417—1.5—0.16440.01062.46680.16436950.01062052.4546071—2—0.16390.01062.47450.1638720.01063262.462460515200.16340.00862.47390.16329970.00862512.4618319——0.10.16310.00702.47630.16302720.00698072.4641318——0.20.16280.00532.47870.16274430.00532962.4665070——0.30.16250.00372.48120.16245000.00367142.4689664


### 4.1. Effect of Magnetic Parameter M

Figure2 shows a decreasing trend in velocity profiles against M2≤M≤8to 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. The 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.Figure 2
Velocity profiles for differentM.Figure 3
Temperature profiles for differentM.Figure 4
Concentration profiles for differentM.


### 4.2. Effects of Darcy Number Da

Darcy number Da characterizes the strength of permeability of the porous medium. Figures5 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. Thus, the thickness of the thermal boundary layer decreases as shown in Figure 7.Figure 5
Velocity profiles for differentDa.Figure 6
Concentration profiles for differentDa.Figure 7
Temperature profiles for differentDa.


### 4.3. Effects of Lewis Number Le

Figure8 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. This 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.Figure 8
Velocity profiles for differentLe.Figure 9
Temperature profiles for differentLe.Figure 10
Concentration profiles for differentLe.


### 4.4. Effects of Grashof Number Gr

The 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 Figures11–13. It can be observed that temperature and concentration decrease with the Grashof number Gr3≤Gr≤6.5, but there is an abrupt change in a velocity profile.Figure 11
Velocity profiles for differentGr.Figure 12
Temperature profiles for differentGr.Figure 13
Concentration profiles for differentGr.


### 4.5. Effects ofNt

Figure14 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.Figure 14
Velocity profiles for differentNt.Figure 15
Temperature profiles for differentNt.Figure 16
Concentration profiles for differentNt.


### 4.6. Effects ofNb

Incremental Brownian parameterNb5≤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.Figure 17
Velocity profiles for differentNb.Figure 18
Concentration profiles for differentNb.Figure 19
Temperature profiles for differentNb.


### 4.7. Effects of Eckert Number Ec

Eckert number plays an important role in high speed flows for which viscous dissipation is significant. It gives relative importance of the kinetic energy in heat transfer flows. ForEc≪1, the energy equation gives the balance between conduction and convection. From Figures 20–22, the effects of this dissipation on velocity, 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. The 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.Figure 20
Velocity profiles for differentEc.Figure 21
Temperature profiles for differentEc.Figure 22
Concentration profiles for differentEc.


## 4.1. Effect of Magnetic Parameter M

Figure2 shows a decreasing trend in velocity profiles against M2≤M≤8to 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. The 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.Figure 2
Velocity profiles for differentM.Figure 3
Temperature profiles for differentM.Figure 4
Concentration profiles for differentM.


## 4.2. Effects of Darcy Number Da

Darcy number Da characterizes the strength of permeability of the porous medium. Figures5 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. Thus, the thickness of the thermal boundary layer decreases as shown in Figure 7.Figure 5
Velocity profiles for differentDa.Figure 6
Concentration profiles for differentDa.Figure 7
Temperature profiles for differentDa.


## 4.3. Effects of Lewis Number Le

Figure8 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. This 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.Figure 8
Velocity profiles for differentLe.Figure 9
Temperature profiles for differentLe.Figure 10
Concentration profiles for differentLe.


## 4.4. Effects of Grashof Number Gr

The 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 Figures11–13. It can be observed that temperature and concentration decrease with the Grashof number Gr3≤Gr≤6.5, but there is an abrupt change in a velocity profile.Figure 11
Velocity profiles for differentGr.Figure 12
Temperature profiles for differentGr.Figure 13
Concentration profiles for differentGr.


## 4.5. Effects ofNt

Figure14 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.Figure 14
Velocity profiles for differentNt.Figure 15
Temperature profiles for differentNt.Figure 16
Concentration profiles for differentNt.


## 4.6. Effects ofNb

Incremental Brownian parameterNb5≤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.Figure 17
Velocity profiles for differentNb.Figure 18
Concentration profiles for differentNb.Figure 19
Temperature profiles for differentNb.


## 4.7. Effects of Eckert Number Ec

Eckert number plays an important role in high speed flows for which viscous dissipation is significant. It gives relative importance of the kinetic energy in heat transfer flows. ForEc≪1, the energy equation gives the balance between conduction and convection. From Figures 20–22, the effects of this dissipation on velocity, 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. The 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.Figure 20
Velocity profiles for differentEc.Figure 21
Temperature profiles for differentEc.Figure 22
Concentration profiles for differentEc.


## 5. Conclusions

This 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. The key observations of present work are as follows:(i)
The skin friction coefficient grows with a rise of the values of magnetic parameter but local Nusselt and Sherwood numbers reduces The 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)
The Darcy number causes the thermal boundary layer and solute concentration to reduce whereas it enhances the momentum boundary layer. The skin friction coefficient decreases yet the wall temperature gradient and nanoparticle concentration increase with an increase in Darcy number.(iii)
The Grashof number enhances local Nusselt number and local Sherwood number whereas this reduces the skin friction coefficient.(iv)
The 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)
The thermophoresis parameter causes both the thermal and momentum boundary layer to increase while its effect on nanoparticle volume fraction is insignificant. The skin friction coefficient and local Sherwood number increase while local Nusselt number remains unaltered.(vi)
The 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)
The Eckert numberEc 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)
The 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 withbvp4c and those with the literature which gives a good account of agreement with each other.



---

*Source: 6628009-2021-01-27.xml*