International Journal of Analysis and Applications Volume 19, Number 5 (2021), 725-742 URL: https://doi.org/10.28924/2291-8639 DOI: 10.28924/2291-8639-19-2021-725 APPLICATION OF SUCCESSIVE LINEARISATION METHOD ON THE BOUNDARY LAYER FLOW PROBLEM OF HEAT AND MASS TRANSFER WITH RADIATION EFFECT AHMED A. KHIDIR, SALIHAH L. ALSHARARI∗ Department of Mathematics, Faculty of Science, University of Tabuk, P.O. Box 741, Tabuk 71491, Saudi Arabia ∗Corresponding author: salihah1972@yahoo.com Abstract. In this paper, we applied the successive linearization method (SLM) in solving highly system of nonlinear boundary value problem. The method is presented in detail by solving the problem of free convective heat and mass transfer of an incompressible fluid past a moving vertical plate in the presence of radiation effect. The governing partial differential equations are converted into system of non linear ordinary differential equations by a similarity transformation, which are converted into system of linear ordinary differential equations using SLM. The linear system is solved using the Chebyshev spectral method to find solutions that are accurate and converge rapidly to the full numerical solution. Comparison with previously published works are performed to test the validity of the obtained results with focus on the accuracy and convergence of the solution. The effects of selected fluid parameters on the velocity as well as the temperature and concentration distribution are determined and discussed. 1. Introduction Most problems that arise in engineering are nonlinear with no analytic solutions and developing new methods that give rapid convergence, are robust and easy to use is a core function of numerical analysis. In the last few decades, a great deal of interest has been generated in the area of heat and mass transfer on a continuously stretching surface with a given temperature or heat flux distribution and they have Received June 20th, 2021; accepted July 26th, 2021; published August 12th, 2021. 2010 Mathematics Subject Classification. 80A20. Key words and phrases. successive linearisation method; boundary layer flow; heat and mass transfer. ©2021 Authors retain the copyrights of their papers, and all open access articles are distributed under the terms of the Creative Commons Attribution License. 725 https://doi.org/10.28924/2291-8639 https://doi.org/10.28924/2291-8639-19-2021-725 Int. J. Anal. Appl. 19 (5) (2021) 726 included several different physical models. Transport of heat is being extensively studied, as understanding the associated transport processes becomes increasingly important. This interest stems from the variety of cases which can be modeled or approximated as transport through porous media, such as packed sphere beds, high performance insulation for buildings, chemical catalytic reactors, grain storage, and many others. Literature concerning convective flow in porous media is abundant, and many published studies such as [3, 6, 7, 11, 14, 16, 17]. Finding exact solutions to steady hydromagnetic flow and heat transfer could be of great benefit to polymer technology. In particular, there are many applications involve the cooling of continuous strips or filaments by drawing them though a quiescent fluid. By drawing such strips in an electrically conducting fluid subjected to a magnetic field, the rate of cooling can be controlled and final products of desired characteristics are obtainable [4]. Numerical study for the effect of thermal-diffusion and diffusion-thermo on combined heat and mass transfer of flows induced a rotating disk has been investigated by Emmanuel et al. [13]. The effect of MHD coupled heat and mass transfer of free convection from a moving permeable vertical was studied by surface was investigated by Yih [18]. A similar problem, including natural convection about a vertical impermeable flat plate, was investigated by Sparrow and Cess [15]. The magneto-hydrodynamics (MHD) convection in porous medium have been studied by Makinde [9, 10]. Alan and Rahman [2] examined Dufour and Soret effects on mixed hydrogenair convective flow past a vertical porous flat plate embedded in a porous medium. Numerical study of natural convection of water in a partially heated enclosure with Soret and Dufour effects were discussed by Nithyadevi and Yang [12]. Recently, Ibrahim and Makinde [5] studied the combined effects of wall suction and magnetic fields on boundary layer flow with heat and mass transfer over an accelerating vertical plate. In this paper, we aims to extend the work of Olanrewaju et. all [1] to include the radiation effect of MHD boundary layer flow of heat and mass transfer. The successive linearisation method (SLM) has been used to convert the governing non linear equations into a system of linear differential equations. The Chebyshev pseudospectral method willbe used to solve the higher order deformation on linear differential equations. The auxiliary linear operator is defined in terms of the Chebyshev spectral collocation differentiation matrix described in [21]. The SLM has been used in a number of recent studies (see [21–25]. They showed that successive linearisation method is accurate and converges rapidly to the numerical results when compared to other recent semi-analytical methods such as the Adomian decomposition method. The SLM method can further be used in place of traditional numerical methods such as finite differences and Runge-Kutta methods to solve highly nonlinear systems of boundary value problems. Int. J. Anal. Appl. 19 (5) (2021) 727 2. Problem formulation We consider the steady free convective heat and mass transfer flow of a viscous, incompressible, electrically conducting fluid past a moving vertical plate in the presence of thermal diffusion (Soret) and diffusion-thermo (Dufour) effects. The non uniform transverse magnetic field B0 is imposed along the y-axis. The induced magnetic field is neglected as the magnetic Reynolds number of the flow is assumed very small. It is further assumed that the external electric field is zero and that the electric field due to charge polarization is negligible. The temperature and the concentration of the ambient fluid are T∞ and C∞ respectively, while those at the surface are, respectively, Tw(x) and Cw(x). The pressure gradient, viscous and electrical dissipation are also neglected. The fluid properties are assumed constant, apart from the density in the buoyancy terms of the linear momentum equation, which is estimated using Boussinesq’s approximation. Under the above assumptions, the boundary layer form of the governing equation can be written as [8]: ∂u ∂x + ∂v ∂y = 0 (2.1) u ∂u ∂x + v ∂u ∂y = v ∂2u ∂2y + gβT (T −T∞) + gβC(C −C∞) − σB20 ρ u (2.2) u ∂T ∂x + v ∂T ∂y = 1 cp ( k + 16σ∗T3∞ 3k∗ ) ∂2T ∂y2 + α ∂2T ∂y2 + DmKT cscp ∂2C ∂y2 (2.3) u ∂C ∂x + v ∂C ∂y = Dm ∂2C ∂y2 + DmKT Tm ∂2T ∂y2 (2.4) The boundary conditions for equations (2.1)–(2.4) are expressed as v = V,u = Bx,T = Tw = T∞ + ax,C = Cw = C∞ + bx at y = 0 u → 0, T → T∞,C → C∞ as y →∞ where B is a constant, a and b denote the stratification rate of the gradient of ambient temperature and concentration profiles and (u,v) are the fluid velocity components in the x and y directions, respectively regarding the plate, T is the temperature,βT is the volumetric coefficient of thermal expansion, α is the thermal diffusivity and g is the acceleration due to gravity. Fluid parameters are ν, the kinematic viscosity, Dm the coefficient of diffusion in the mixture, C the species concentration, σ the electrical conductivity, kT is the themal diffusion ratio, cs is the concentration susceptibility, cp is the specific heat at constant pressure and Tm the mean fluid temperature, k ∗ is the mean absorption coefficient and σ∗ is the Stefan–Boltzmann constant. B0 is the externally imposed magnetic field in the y direction. We introduce the following non- dimensional variables: η = √ B v y F(η) = ψ x √ Bv , θ(η) = T −T∞ Tw −T∞ , φ(η) = C −C∞ Cw −C∞ (2.5) where F(η) is a dimensionless stream function, θ(η) is a dimensionless temperature of the fluid in the boundary layer region, φ(η) is a dimensionless species concentration of the fluid in the boundary layer region Int. J. Anal. Appl. 19 (5) (2021) 728 and η is the similarity variable. The velocity components u and v are respectively obtained as follows u = ∂ψ ∂y = xBF ′, v = − ∂ψ ∂x = − √ BvF, (2.6) where Fw = V √ Bv is the dimensionless suction velocity. Following equation (2.6), the partial differential equations (2.2)–(2.4) are transformed into local similarity equations as follows: F ′′′(η) + F(η)F ′′(η) − (F ′(η) + M)F ′(η) + Grθ(η) + Gcφ(η) = 0 (2.7)( 1 + 4 3 Rd ) θ′′(η) + PrFθ ′(η) −PrF ′(η)θ(η) + PrDuφ′′(η) = 0 (2.8) φ′′(η) + ScFφ ′(η) −ScF ′(η)φ(η) + ScSrθ′′(η) = 0 (2.9) The boundary conditions are also transformed into the form F ′ = 1 F = −Fw,θ = 1,φ = 1 at η = 0 (2.10) F ′ = 0 θ = 0 φ = 0 as η →∞ (2.11) where M = σB20 ρB is the magnetic parameter, Pr = v α is the Prandtl number, Sc = v Dm is the Schmidt number, Gr = gβT (Tw −T∞) xB2 is the local temperature Grashof number, Gc = gβc(Cw −C∞) xB2 is the local concentration Grashof number, Du = DmkT (Cw−C∞) cscp(Tw−T∞) is the Dufour number, Sr = DmkT (Tw−T∞) Tm(Cw−C∞)v is the Soret number and Rd = 4σ ∗T3∞/k ∗k is the radiation number. 3. Successive Linearisation Method (SLM) The SLM procedure linearises the governing non linear equation. In order to fully describe of the SLM algorithm, let us consider the following boundary value problem of order n in the form: L[u(x),u′(x),u′′(x), ...,un + N[u(x),u′(x),u′′(x), ...,un] = g(x), (3.1) where L and N are linear and non linear operators ,u(x) is an unknown function to be determined, g(x) is a known function. We assume that equation (3.1) is to be solve for x ∈ [a,b] subject to the boundary conditions u(a) = a0, u(b) = b0 (3.2) We represent the vertical difference between the function u(x) and the initial guess u0(x) by a function U1(x) as (see Figure 1) U1(x) = u(x) −u0(x), (3.3) where U1(x) is an unknown functions and u0(x) is the initial guess which is chosen to satisfy boundary conditions (3.2). It is reasonable to assume, for example, that the initial approximation u0(x) is a linear Int. J. Anal. Appl. 19 (5) (2021) 729 Figure 1. Geometric representation of the function ui(x) function in case of second order problems defined on a finite domain and an exponential function for problems defined on an infinite or a semi-infinite domain. Substituting (3.3) into (3.1), yields L[U1(x),U ′ 1(x),U ′′ 1 (x), ....,U (n) 1 ] + L[u0(x),u ′ 0(x),u ′′ 0 (x), ....,u (n) 0 ] + N[U1(x) + u0(x),U ′ 1(x) + u ′ 0(x),U ′′ 1 (x) + u ′′ 0 (x), ...,U (n) 1 (x) + u (n) 0 (x)] = g(x), (3.4) This equation is non-linear in U1(x), so it may not be possible to find an exact solution. We therefore look for a solution which is obtained by solving the linear part of the equation and neglecting the non-linear terms containing U1(x) and its derivatives. We further assume that U1(x) and its derivatives are very small and denote the solution of the linearised equation (3.4) by U1(x), that is U1(x) ≈ u1(x). Equation (3.4) can be written as L[u1(x),u ′ 1(x),u ′′ 1 (x), ....,u (n) 1 ] + f0,0u1(x) + f1,0u ′ 1(x) + f2,0u ′′ 1 (x) + .... + fn,0u (n) 1 (x) = R1(x) (3.5) where f0,0 = ∂N ∂u1(x) (u0(x),u ′ 0(x),u ′′ 0 (x), ....,u (n) 0 ), f1,0 = ∂N ∂u′1(x) (u0(x),u ′ 0(x),u ′′ 0 (x), ....,u (n) 0 ), f2,0 = ∂N ∂u′′1 (x) (u0(x),u ′ 0(x),u ′′ 0 (x), ....,u (n) 0 ), fn,0 = ∂N ∂u (n) 1 (x) (u0(x),u ′ 0(x),u ′′ 0 (x), ....,u (n) 0 ), R1(x) = g(x) −L[u0(x),u′0(x),u ′′ 0 (x), ....,u (n) 0 ] −N[u0(x),u ′ 0(x),u ′′ 0 (x), ....,u (n) 0 ]. Int. J. Anal. Appl. 19 (5) (2021) 730 Since the left hand side of equation (3.5) is linear and the right hand side is known, the equation can be solved for u1(x) subject to the boundary conditions u(a) = 0,u(b) = 0 (3.6) Assuming that the solution of the linear equation (3.5) is close to the solution of the non linear equation (3.4), then the first approximation of the solution (order 1) is u(x) ≈ u0(x) + u1(x) (3.7) To improve this solution we define the vertical difference between the function U1(x) and u1(x) by the function U2(x) as U2(x) = U1(x) −u1(x) (3.8) Substitute (3.8) into equation (3.1) to give L[U2(x),U ′ 2(x),U ′′ 2 (x), ....,U (n) 2 ] + L[u0(x) + u1(x),u ′ 0(x) + u ′ 1(x),u ′′ 0 (x) + u ′′ 1 (x), ....,u (n) 0 + u (n) 1 (x)] + N[U2(x) + u0(x) + u1(x),U ′ 2(x) + u ′ 0(x) + u ′ 1(x) ′U′′2 (x) + u ′′ 0 (x) + u ′′ 1 (x), ...,U (n) 2 (x) + u (n) 0 (x) + u (n) 1 (x)] = g(x) (3.9) Since u0(x) and u1(x) are known and this equation is non-linear in U2(x), we solve the linearized equation after neglecting the non-linear terms containing U2(x) and its derivatives. We further assume that U2(x) and its derivatives are very small that is U2(x) ≈ u2(x). Equation(3.9)can be written as L[u2(x),u ′ 2(x),u ′′ 2 (x), ....,u (n) 2 (x)] + f0,0u2(x) + f1,0u ′ 2(x) + f2,0u ′′ 2 (x) + ..... + fn,0u (n) 2 (x) = R2(x), (3.10) where f0,0 = ∂N ∂u2(x) (u0(x) + u1(x),u ′ 0(x) + u ′ 1(x),u ′′ 0 (x) + u ′′ 1 (x), ....,u (n) 0 (x) + u n 1 (x)), f1,0 = ∂N ∂u′2(x) (u0(x) + u1(x),u ′ 0(x) + u ′ 1(x),u ′′ 0 (x) + u ′′ 1 (x), ....,u (n) 0 (x) + u n 1 (x)), f2,0 = ∂N ∂u′′1 (x) (u0(x) + u1(x),u ′ 0(x) + u ′ 1(x),u ′′ 0 (x) + u ′′ 1 (x), ....,u (n) 0 (x) + u n 1 (x)), fn,0 = ∂N ∂u (n) 1 (x) (u0(x) + u1(x),u ′ 0(x) + u ′ 1(x),u ′′ 0 (x) + u ′′ 1 (x), ....,u (n) 0 (x) + u n 1 (x)), R2(x) = g(x) −L[u0(x) + u1(x),u′0(x) + u ′ 1(x), ....,u (n) 0 (x) + u (n) 1 (x)] −N[u0(x) + u1(x),u′0(x) + u ′ 1(x), ....,u (n) 0 (x)]. After solving the equation (3.10), the second order approximation of u(x) is given by u(x) ≈ u0(x) + u1(x) + u2(x). (3.11) Int. J. Anal. Appl. 19 (5) (2021) 731 This process is repeated for m = 2, 3, ..., i. Thus, u(x) is given by u(x) = U1(x) + u0(x), = U2(x) + u0(x) + u1(x), = U3(x) + u0(x) + u1(x) + u2(x), ... = Ui+1(x) + u0(x) + u1(x) + u2(x) + ... + ui(x), = Ui+1(x) + i∑ m=0 um(x). Thus, for large i, we can approximate the ith order solution u(x) by u(x) = i∑ m=0 um(x). (3.12) The solution ui(x) can be determined from the linearized original equation (3.1) starting from the initial guess u0(x) and solving the linear equations for ui(x). In general, the form of the linearized equation for ui(x) is given by L[ui(x),u ′ i(x),u ′′ i (x), ....,u (n) i ] + f0,i−1ui(x) + f1,i−1u ′ i(x) + f2,i−1u ′′ i (x) + ..... + fn,i−1u (n) i (x) = Ri−1(x), i = 1, 2, ...,M, (3.13) subject to the boundary conditions ui(a) = 0,ui(b) = 0, (3.14) where M is termed the order of the SLM. f0,i−1 = ∂N ∂ui(x) ( i−1∑ m=0 um(x), i−1∑ m=0 u′m(x), i−1∑ m=0 u′′m(x), ...., i−1∑ m=0 u(n)m (x)), f1,i−1 = ∂N ∂u′i(x) ( i−1∑ m=0 um(x), i−1∑ m=0 u′m(x), i−1∑ m=0 u′′m(x), ...., i−1∑ m=0 u(n)m (x)), f2,i−1 = ∂N ∂u′′i (x) ( i−1∑ m=0 um(x), i−1∑ m=0 u′m(x), i−1∑ m=0 u′′m(x), ...., i−1∑ m=0 u(n)m (x)), Int. J. Anal. Appl. 19 (5) (2021) 732 ... fn,i−1 = ∂N ∂u (n) i (x) ( i−1∑ m=0 um(x), i−1∑ m=0 u′m(x), i−1∑ m=0 u′′m(x), ...., i−1∑ m=0 u(n)m (x)), Ri−1(x) = g(x) −L( i−1∑ m=0 um(x), i−1∑ m=0 u′m(x), i−1∑ m=0 u′′m(x), ...., i−1∑ m=0 u(n)m (x)) −N( i−1∑ m=0 um(x), i−1∑ m=0 u′m(x), i−1∑ m=0 u′′m(x), ...., i−1∑ m=0 u(n)m (x)). Then ordinary differential equation is linear and can easily be solved using any analytical or numerical method. 4. Method of solution The system of non linear equation to be solved usin SLM are F ′′′ + FF ′′ − (F ′ + M)F ′ + Grθ + Gcφ = 0 (4.1) λθ′′ + Fθ′ −F ′θ + Duφ′′ = 0 (4.2) φ′′ + ScFφ ′ −ScF ′φ + ScSrθ′′ = 0 (4.3) subject to the boundary conditions F ′ = 1 F = −Fw,θ = 1,φ = 1 at η = 0 (4.4) F ′ = 0 θ = 0 φ = 0 as η →∞ (4.5) where λ = 1 Pr ( 1 + 4 3 Rd ) . The SLM algorithm starts with the assumption that the independent variables F(η),θ(η) and φ(η) can be expressed as F(η) = Fi(η) + i−1∑ m=0 Fm(η), θ(η) = θi(η) + i−1∑ m=0 θm(η), φ(η) = φ(η) + i−1∑ m=0 φm(η) (4.6) where Fi,θi,φi (i = 1, 2, 3, ...) are unknown functions and Fm,θm and φm(m < i) are known functions which are obtained by recursively solving the linear part of the equation system that results from substituting in the governing equations. The main assumption of the SLM is that fi,θi,φi become increasingly smaller when i becomes large that is lim i→∞ Fi = lim i→∞ θi = lim i→∞ φ = 0. (4.7) Thus, starting from the initial guesses F0(η) = 1 − Fw − e−η,θ0(η) = e−η,φ0(η) = e−η, which are chosen to satisfy boundary conditions (4.4) and (4.5), the subsequent solutions for Fi,φi,θi, i ≥ 1 are obtained by successively solving the linearised form of the equations which are obtained by substituting equation (1)in the governing equations and considering only the linear terms. The linearised equations to be solved are Int. J. Anal. Appl. 19 (5) (2021) 733 given by substituting the assumptions (4.6) in the system (4.1) - (4.3) and neglecting the non linear terms, gives F ′′′i + i−1∑ m=0 FmF ′′ i + i−1∑ m=0 F ′′mFi − 2 i−1∑ m=0 F ′mF ′ i −MF ′ i + Crθi + Gcrφi = r1 (4.8) λθ′′i + i−1∑ m=0 θ′m(η)Fi + i−1∑ m=0 Fm(η)θ ′ i − i−1∑ m=0 θmF ′ i − i−1∑ m=0 F ′mθi + Duφ ′′ i = r2 (4.9) φ′′i + Sc i−1∑ m=0 φ′mFi + Sc i−1∑ m=0 Fmφ ′ i −Sc i−1∑ m=0 φmF ′ i + Sc i−1∑ m=0 F ′mφi + Scθ ′′ i = r3 (4.10) Subject to the boundary conditions Fi(0) = F ′ i (0) = F ′ i (∞) = θi(0) = θi(∞) = φ(0) = φi(∞) = 0 (4.11) where r1 = − ( i−1∑ m=0 F ′′′m + i−1∑ m=0 Fm i−1∑ m=0 F ′′m − i−1∑ m=0 F ′m i−1∑ m=0 F ′m −M i−1∑ m=0 F ′m + Gr i−1∑ m=0 θm + Gc i−1∑ m=0 φm ) r2 = − ( λ i−1∑ m=0 θ′′m + i−1∑ m=0 Fm(η) i−1∑ m=0 θ′m − i−1∑ m=0 F ′m t−1∑ m=0 θm + Du i−1∑ m=0 φ′′m ) r3 = − ( i−1∑ m=0 φ′′m + Sc i−1∑ m=0 Fm i−1∑ m=0 φ′m −Sc i−1∑ m=0 F ′m i−1∑ m=0 φm + ScSr i−1∑ m=0 θ′′m ) once each solution for fi,θi,φi(i ≥ 1) has been found from iteratively solving equations (3- 6 ), the approxi- mate solutions for f(η),θ(η) and φ(η) are obtained as f(η) ≈ M∑ m=0 fm(η),θ(η) ≈ M∑ m=0 θm(η),φ(η) ≈ M∑ m=0 φm(η), (4.12) Where M is the order of SLM approximation. Since the coefficient parameters and the right hand side of equations (4.8) - (4.10), for i = 1, 2, 3, ... are known (from previous iterations). The equation system can easily be solved analytically (wherever possible) or using any numerical method such as finite differences, finite elements, Runge-Kutta based shooting methods or collocation methods’ in this work, we used the Chebyshev spectral collocation method. This method is based on approximating the unknown functions by the Chebyshev interpolating polynomials in such a way that they are collocated at the Gauss-Lobatto points defined as (see [19, 20]) ξj = cos πj N , j = 0, 1, . . . ,N where N + 1 is the number of collocation points used. In order to implement the method, the physical region [ 0,∞) is transformed into the region [−1, 1] using the any domain truncation technique in which the problem is solved on the interval instead of [0,∞). The unknown functions fi,θi and φi are approximated Int. J. Anal. Appl. 19 (5) (2021) 734 at the collocation points by fi(ξ) ≈ N∑ k=0 fi (ξ z k) Tk (ξj) ,θi(ξ) ≈ N∑ k=0 θi (ξk) Tk (ξj) ,φ(ξ) ≈ N∑ k=0 φ (ξk) Tk (ξj) ,j = 0, 1, . . .N (4.13) Where Tkis the kth Chebyshev polynomial defined as Tk(ξ) = cos[k cos −1 (ξ)] (4.14) The derivatives of the variables at the collocation points are represented as drFi dηr = N∑ k=0 Drkj ·Fi (ξk) , drθi dηr = N∑ k=0 Drkjθi (ξk) , drφi dηr = N∑ k=0 Drkjφi (ξk) ,j = 0, 1, . . . ,N (4.15) Where r is the order of differentiation and D being the Chebyshev spectral differentiation matrix hose entries are defined as [19, 20] D00 = 2 N 2+1 6 Djk = cj(−1)j+k ckξj−ξk j = /k; j,k = 0, 1, . . . ,N Dkk = − ξk 2(1−ξ2k) ,k = 1, 2, . . . ,N − 1, DNN = −2 N 2+1 6   (4.16) Substituting equations (4.15) into the system (4.8) - (4.10) gives the following system of algebraic equations A11Fi + A12θi + A13φi = r1,i−1 A21Fi + A22θi + A23φi = r2,i−1 A31Fi + A32θi + A33φi = r2,i−1 which leads to the matrix equation given as Ai−1Xi = Ri−1, (4.17) in which Ai−1 is a (3N + 3)×(3N + 3) square matrix and Xi and Ri−1 are (3N + 1)1 column vectors defined by Ai−1 =   A11 A12 A31 A21 A22 A32 A31 A32 A33   , Xi =   Fi θi φi   , Ri−1 =   r1,i−1 r2,i−1 r3,i−1   (4.18) where Fi = [ Fi (η0) ,Fi (η1) , . . . ,Fi (ηN−1) ,Fi (ηN )] T θi = [ θi (η0) ,θi (η1) , . . . ,θi (ηN−1) ,θi (ηN )] T φi = [ φi (η0) ,φi (η1) , . . . ,φi (ηN−1) ,φi (ηN )] T r1,i−1 = [r1,i−1 (η0) ,r1,i−1 (η1) , . . . ,r1,i−1 (ηN−1) ,r1,i−1(ηN )] T Int. J. Anal. Appl. 19 (5) (2021) 735 r2,i−1 = [r2,i−1 (η0) ,r2,i−1 (η1) , . . . ,r2,i−1 (ηN−1) ,r2,i−1(ηN )] T r3,i−1 = [r3,i−1 (η0) ,r3,i−1 (η1) , . . . ,r3,i−1 (ηN−1) ,r3,i−1(ηN )] T A11 = D 3 + i−1∑ m=0 FmD 2 − ( 2 i−1∑ m=0 F ′m + M ) D + [ i−1∑ m=0 F ′′m ] , A12 = [Cr], A13 = [Gr], A21 = − i−1∑ m=0 θmD + [ i−1∑ m=0 θ′m(η) ] ,A22 = λD 2 + i−1∑ m=0 FmD − [ i−1∑ m=0 F ′m ] ,A23 = DuD 2, A31 = −Sc i−1∑ m=0 φmD + Sc [ i−1∑ m=0 φ′m(η) ] ,A32 = SrScD 2,A33 = D 2 and [..] stands to a diagonal matrix of size (3N + 3)×(3N + 3). The boundary conditions (4.15) are imposed to the system (4.17) as displayed in figure 2. After modifying the matrix system (4.17) to incorporate boundary conditions, the solution is obtained as Xi = A −1 i−1Ri−1 (4.19) Figure 2. Imposing the boundary conditions (4.15) into the system (4.19). Int. J. Anal. Appl. 19 (5) (2021) 736 5. Results and discussion The system of non-linear ordinary differential equations (2.7) - (2.9) together with the boundary conditions (2.10) and (2.11) are solved numerically using successive linearization method. In this chapter we give the obtained results for the main parameters affecting the flow. We use MATLAB software in all obtained results in this research. To check the accuracy of the proposed successive linearisation method (SLM), comparison was made with those obtained in literature. The graphs and tables presented in this work, unless otherwise specified, were generated using N = 30 and η∞ ' L = 15 gave sufficient accuracy for the SLM results. 6. Convergence of the solution In this section, comparisons with previously published works are performed to test the validity and the convergence of the obtained results. Also, the effects of some physical parameters on the Velocity, temperature and concentration Profiles are obtained and discussed. Table 1 showed the SLM results of F ′′(0), −θ′(0) and −φ′(0) at various values of Rd for different iterations. It is cleared form the Table that the results obtained by SLM are in excellent agreement with a few order of SLM series starting from the third iteration and giving accuracy up to eight decimal places.. In Table 2, the SLM results were compared for −F ′′(0), −θ′(0) and −φ′(0) for different values of Fw with those obtained by Olanrewaju et. all [1] in the absence of radiation effect. They used the shooting iteration technique together with a sixth order Runge–Kutta integration scheme. It can be seen from the table that they are in good agreement between the results.. 7. Velocity Profiles Figure 3 illustrates the effect of magnetic parameter M parameter on Velocity as a function of the similarity variable η. It is observed from the Figure that an increase in the magnetic parameter M leads to decreases in the velocity profile. The same result was obtained by Olanrewaju et. all [1]. Figure 4 represents the the effect of the radiation parameter Rd on dimensionless velocity distributions. It shows that the velocity enhanced as the radiation parameter increased. 8. Temperature Profiles The effect of magnetic filed parameter M on temperature distribution profile θ(η) as a function of variable η in displayed on Figure 5. From the Figure, we see that an increase in the magnetic parameter M leads to increases in the temperature distribution. Figure 6 shows the effect of increasing the radiation parameter Rd on dimensionless temperature distribu- tions. We observe that increasing the radiation parameter leads to enhanced the temperature distributions θ(η) Int. J. Anal. Appl. 19 (5) (2021) 737 Table 1. The SLM results for F ′′(0), −θ′(0) and −φ′(0) to different values of the radiation number Rd when Fw = 1,M = 0.1,Gr = 0.1,Gc = 0.1,Pr = 0.72,Sc = 0.62,Sr = 0.1,Du = 0.03 Rd First iteration Second iteration Third iteration Fourth iteration F ′′(0) 0.1 0.56946888 0.56947404 0.56947405 0.56947405 0.2 0.56880563 0.56881060 0.56881061 0.56881061 0.5 0.56717266 0.56717699 0.56717699 0.56717699 1.0 0.56523532 0.56523914 0.56523914 0.56523914 −θ′(0) 0.1 0.52976296 0.52974957 0.52974955 0.52974955 0.2 0.50638146 0.50636189 0.50636185 0.50636185 0.5 0.45023485 0.45019930 0.45019924 0.45019924 1.0 0.38601330 0.38596450 0.38596442 0.38596442 −φ′(0) 0.1 0.51484762 0.51483002 0.51482998 0.51482998 0.2 0.51622774 0.51621125 0.51621122 0.51621122 0.5 0.51958417 0.51957245 0.51957243 0.51957243 1.0 0.52347723 0.52347206 0.52347205 0.52347205 9. Concentration Profiles Figure 7 showed the effect of magnetic parameter M parameter on concentration Profile as a function of the similarity variable η. It is observed from the Figure that an increase in the magnetic parameter M leads to increases in the Concentration Profile φ(η). The effect of radiation parameter Rd is decreased concentration distribution via increasing Rd . Int. J. Anal. Appl. 19 (5) (2021) 738 Table 2. Comparison between the SLM, and Ref. [1] for −F ′′(0), −θ′(0) and −φ′(0) for different values of Fw when M = 0.1,Gr = 0.1,Gc = 0.1,Pr = 0.72,Rd = 0,Sc = 0.62,Sr = 0.1,Du = 0.03 Ref. [1] SLM Fw −F ′′(0) −θ′(0) −φ′(0) −F ′′(0) −θ′(0) −φ′(0) 0.5 0.724431 0.673004 0.605696 0.724428 0.673003 0.605696 0.3 0.801102 0.728095 0.648537 0.801097 0.728093 0.648536 0.0 0.934398 0.820892 0.720212 0.934389 0.820888 0.720212 −0.3 1.090854 0.926778 0.801866 1.090838 0.926772 0.801864 −0.5 1.208097 1.004942 0.862315 1.208074 1.004934 0.862312 0 5 10 15 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 η F ′ ( η ) M = 0 M = 0.3 M = 0.5 M = 1 Figure 3. Effect of M on the Velocity profile when Fw = 0.1,Gr = 0.1,Gc = 0.1,Pr = 0.72,Rd = 0.1,Sc = 0.62,Sr = 0.1,Du = 0.03 . 10. Conclusion In this paper, we applied the successive linearization method in solving highly system of nonlinear bound- ary value problem. The method is applied on the MHD free convective heat and mass transfer with radiation effect. The set of governing equations and the boundary conditions are reduced to ordinary differential equa- tions with appropriate boundary conditions. The results were compared with other methods in the literature Int. J. Anal. Appl. 19 (5) (2021) 739 0 2 4 6 8 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 η F ′ ( η ) Rd = 0.0 Rd = 0.5 Rd = 1.0 Rd = 1.5 Figure 4. Effect of M on the concentration profile when Fw = 0.1,Gr = 0.1,Gc = 0.1,Pr = 0.72,M = 0.1,Sc = 0.62,Sr = 0.1,Du = 0.03 . 0 5 10 15 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 η θ (η ) M = 0 M = 0.3 M = 0.5 M = 1 Figure 5. Effect of M on the temperature profile when Fw = 0.1,Gr = 0.1,Gc = 0.1,Pr = 0.72,Rd = 0.1,Sc = 0.62,Sr = 0.1,Du = 0.03 . such as the shooting iteration technique together with a sixth order Runge–Kutta integration scheme with focus on the accuracy and convergence of the results. Graphs were presented showing the effects of various physical parameters on the fluid properties. The main conclusions emerging from this paper are as follows: • The SLM suggested a standard method for choosing the linear operators and initial approximations by using any form of initial guess as long as it satisfies the boundary conditions, while the initial Int. J. Anal. Appl. 19 (5) (2021) 740 0 5 10 15 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 η θ (η ) Rd = 0.0 Rd = 0.5 Rd = 1.0 Rd = 1.5 Figure 6. Effect of M on the concentration profile when Fw = 0.1,Gr = 0.1,Gc = 0.1,Pr = 0.72,M = 0.1,Sc = 0.62,Sr = 0.1,Du = 0.03 . 0 5 10 15 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 η φ (η ) M = 0 M = 0.3 M = 0.5 M = 1 Figure 7. Effect of M on the concentration profile when Fw = 0.1,Gr = 0.1,Gc = 0.1,Pr = 0.72,Rd = 0.1,Sc = 0.62,Sr = 0.1,Du = 0.03 . guess in the other method such as HPM and HAM can be selected that will make the integration of the higher order deformation equations possible. • The SLM is highly accurate, efficient and converges rapidly with a few iterations required to achieve the accuracy of the numerical results, in this study it was found that for a few iterations of SLM was sufficient to give good agreement with the exact solution. Int. J. Anal. Appl. 19 (5) (2021) 741 0 2 4 6 8 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 η φ (η ) Rd = 0.0 Rd = 0.5 Rd = 1.0 Rd = 1.5 Figure 8. Effect of M on the concentration profile when Fw = 0.1,Gr = 0.1,Gc = 0.1,Pr = 0.72,M = 0.1,Sc = 0.62,Sr = 0.1,Du = 0.03 . • When the magnetic field increases, the velocity profile is decreased while the temperature and con- centration components are enhanced. • When the radiation parameter increases, the velocity profile temperature distribution are increased while the concentration component is reduced. Finally, the successive linearisation method has high accuracy and simple for solving nonlinear boundary value problems compared with the Runge Kutta, finite difference, finite element and Keller-Box methods. Because of its efficiency and easy of use. The extension to systems of nonlinear BVPs allows the method to be used as alternative to the traditional of those methods. Conflicts of Interest: The author(s) declare that there are no conflicts of interest regarding the publication of this paper. References [1] P. O. Olanrewaju, O. D. Makinde, Effects of thermal diffusion and diffusion thermo on chemically reacting MHD boundary layer flow of heat and mass transfer past a moving vertical plate with suction/injection, Arab. J. Sci. Eng. 36 (2011), 1607-1619. [2] M. S. Alam, M. M Rahman, Dufour and Soret effects on mixed convection flow past a vertical porous flat plate with variablesuction, Nonlinear Anal.: Model. Control, l1 (1) (2006), 3–12. [3] A. Bejan, I. Dincer, S. Lorente, A. F. Miguel, A. H. Reis, Porous and complex flow structures in modern technologies. Springer, New York, (2004). [4] A. Chakrabarti, A. S. Gupta, Hydromagnetic flow and heat transfer over a stretching sheet. Quart. Appl. Math. 37 (1979), 73–78. Int. J. Anal. Appl. 19 (5) (2021) 742 [5] S. Y. Ibrahim, O. D. Makinde, Chemically reacting MHD boundary layer flow of heat and mass transfer past a moving vertical plate with suction. Sci. Res. Essay. 5 (19) (2010), 2875–2882. [6] D. B. Ingham, A. Bejan, E. Mamut, I. Pop, Emerging technologies and techniques in porous media, Kluwer, Dordrecht, (2004). [7] D. B. Ingham, I. Pop, Transport phenomena in porous media, vol. III. Pergamon, Oxford, (2005). [8] N. G. Kafoussias, E. W. Williams, Thermal-diffusion and diffusion-thermo effects on mixed free convective and mass transfer boundary layer flow with temperature dependent viscosity. Int. J. Eng. Sci. 33 (1995), 1369–1384. [9] O. M. Makinde, MHD steady flow and heat transfer on the sliding plate. AMSE Model. Meas. Control B. 70(1) (2001), 61–70. [10] O. D. Makinde, On MHD boundary-layer flow and mass transfer past a vertical plate in a porous medium with constant heat flux. Int. J. Numer. Meth. Heat Fluid Flow, 19 (3/4) (2009), 546–554. [11] D. A. Nield, A. Bejan, Convection in porous media, 3rd edn. Springer, New York, (2006). [12] N. Nithyadevi, R. J. Yang, Double diffusive natural convection in a partially heated enclosure with Soret and Dufour effects. Int. J. Heat Fluid Flow, 30 (5) (2009), 902-910. [13] E. Osalusi, J. Side, R. Harris, Thermal-diffusion and diffusion-thermo effects on combined heat and mass transfer of a steady MHD convective and slip flow due to a rotating disk with viscous dissipation and Ohmic heating. Int. Commun. Heat Mass Transfer. 35 (2008), 908–915. [14] I. Pop, D. B. Ingham, Convective heat transfer: mathematical and computational modeling of viscous fluids and porous media. Pergamon, Oxford, (2001). [15] E. M. Sparrow, R. D. Cess, The effect of a magnetic field on free convection heat transfer, Int. J. Heat Mass Transfer. 3 (1961), 267–274. [16] P. Vadasz, Emerging topics in heat and mass transfer in porous media. Springer, New York (2008). [17] K. Vafai, Handbook of porous media. Taylor & Francis, New York (2005). [18] K. A. Yih, Free convection effect on MHD coupled heat and mass transfer of a moving permeable vertical surface. Int. Commun. Heat Mass Transf. 26(1), 95–104 (1999). [19] C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Zang, Spectral Methods in Fluid Dynamics, Springer-Verlag, Berlin (1998). [20] L. N. Trefethen, Spectral Methods in MATLAB, SIAM, 2000. [21] Z. G. Makukula, P. Sibanda and S. S. Motsa, A note on the solution of the Von Karman equations using series and chebyshev spectral methods, Bound. Value Probl. 2010 (2010), 471793. [22] Z. Makukula and S. S. Motsa, On new solutions for heat transfer in a viscolastic fluid between parallel plates, Int. J. Math. Models Meth. Appl. Sci. 4 (2010), 221-230. [23] S. S. Motsa, P. Sibanda and S. Shateyi, On a new quasi-linearization method for systems of nonlinear boundary value problems. Math. Meth. Appl. Sci. 34 (2011), 1406–1413. [24] F. G. Awad, P. Sibanda, S. S. Motsa and O. D. Makinde, Convection from an inverted cone in a porous medium with cross-diffusion effects, Computers Math. Appl. 61 (2011), 1431-1441. [25] F. G. Awad, P. Sibanda, M. Narayana and S. S. Motsa, Convection from a semi-finite plate in a fluid saturated porous medium with cross-diffusion and radiative heat transfer, Int. J. Phys. Sci. 6 (21) (2011), 4910–4923. 1. Introduction 2. Problem formulation 3. Successive Linearisation Method (SLM) 4. Method of solution 5. Results and discussion 6. Convergence of the solution 7. Velocity Profiles 8. Temperature Profiles 9. Concentration Profiles 10. Conclusion References