Numerical Solutions for Laminar Boundary Layer Nanofluid Flow along with a Moving Cylinder with Heat Generation, Thermal Radiation, and Slip Parameter

The investigation of the numerical solution of the laminar boundary layer flow along with a moving cylinder with heat generation, thermal radiation, and surface slip effect is carried out. The fluid mathematical model developed from the Navier-Stokes equations resulted in a system of partial differential equations which were then solved by the multidomain bivariate spectral quasilinearization method (MD-BSQLM). The results show that increasing the velocity slip factor results in an enhanced increase in velocity and temperature profiles. Increasing the heat generation parameter increases temperature profiles; increasing the radiation parameter and the Eckert numbers both increase the temperature profiles. The concentration profiles decrease with increasing radial coordinate. Increasing the Brownian motion and the thermophoresis parameter both destabilizes the concentration profiles. Increasing the Schmidt number reduces temperature profiles. The effect of increasing selected parameters: the velocity slip, Brownian motion, and the radiation parameter on all residual errors show that these errors do not deteriorate. This shows that the MD-BSQLM is very accurate and robust. The method was compared with similar results in the literature and was found to be in excellent agreement.


Introduction
The boundary layer flow on heat and mass transfer has been overmoving, and stretching surfaces have been studied and remained an active area in the past decade. This is because it has numerous applications in areas such as hot rolling, processes of polymer extrusions, wire drawing, extrusions in aerodynamic plastic sheets, process of condensation in metallic plates during cooling, and many other applications. According to Poply et al. [1], the study of flows in cylinders is considered two-dimensional when the cylinder radius is much larger than the boundary layer thickness. For lean and thin cylinders, the two dimensions may be of the same order; in this case, the flow is referred to as axisymmetric rather than two-dimensional. These dimensions affect velocity, temperature, and concentration profiles which in turn affect the skin friction coefficient.
In light of the importance of laminar boundary layer fluid flow, many researchers have considered several flow geometries with different boundary conditions, using many different techniques to solve similar models. These include the research carried out by Shateyi and Marewo [2] who used the successive linearization method (SRM) to solve a problem on laminar boundary layer flow and heat transfer in nonlinear differential equations; they also considered stretching cylinder, porous media, and thermal conductivity. In their investigation, they observed that the curvature significantly affects temperature and velocity fields. Also, both the skin friction coefficient and local Nusselt number increase as the curvature increases. Rangi and Naseem [3] used the Keller-box technique to solve the equations describing boundary layer flow of heat transfer with nonconstant thermal conductivity along with a stretching cylinder. Their results also showed that the cylinder curvature affects temperature, velocity, and skin friction fields. In the results, thermal conductivity and curvature aid heat transfer and reduction in fluid viscosity aids the convectional heat transfer rate. Numerical solutions concerning boundary layer flow, heat transfer along a stretching cylinder in porous media were obtained by Mukhopadhyay [4] using the shooting method. The results of this study found that increasing permeability parameter results in a decrease in velocity profiles. The skin friction and rate of heat transfer are much less on the flat plate than on the cylinder surface. Some other laminar boundary layer fluid models by different researchers including Elbashbeshy et al. [5], Lin and Shih [6], and Ali and Alabdulkarem [7] among others investigated the Casson fluid flow on a stretching surface based on an exponential model; they also applied theoretical analysis using lie groups in MHD fluid flow. In these works, a fourth-order Runge-Kutta method was used.
The study of the fluid flow and the physical properties of nanofluids has been widely studied in the past few years due to the wide applications of these fluids in cancer therapy, fuel cells, electronics, and pharmaceutical processes, just to mention a few. The concept of nanofluids was introduced by Choi [8], where the introduction for the proposal of nanoparticles suspended in base fluids such as ethylene glycol, oil, and water. A nonhomogeneous equation with two components for nanofluids was introduced by Buongiorno [9], where the seven slip mechanisms were proposed. The results explained effects in thermophoresis Brownian motion in nanofluids. Alamri et al. [10] used the homotopy analysis method to investigate fluid flow in nanofluids in the presence of secondorder slip and Stefan blowing effects in a tube (Poiseuille). In this investigation, the results show that the retarding effects of Stefan blowing is observed for temperature and velocity profiles, the opposite effect was noticed for the case of particle concentrations. A more enhanced response in the field of velocity was observed in the slip of the second order than in the slip of the first order. Nadeem et al. [11] used the homotopy analysis method in the solutions for nonlinear differential equations in which the oblique stagnation point flow was considered for a Casson nanofluid flow with stretching surface and heat transfer. The Runge-Kutta method of the fourth order together with the shooting technique was used by Khan et al. [12] to find a numerical solution of the Maxwell nanofluid stagnation point fluid flow over a stretching sheet with slip conditions and chemical reaction effects. Their results showed that the skin friction coefficient is inversely proportional to slip parameter but an opposite is noticed in the case of fluid relaxation parameter.
Dhanai et al. [13] investigated MHD boundary layer fluid flow on multiple solutions of and heat transfer power-law nanofluids, viscous dissipation, and permeable shrinking/stretching effects in which the shooting method was used. In this study, viscous dissipation was significant and the Brownian motion was neglected in the case of heat transfer. The results of the analysis of the mass and heat transfer in nanofluid flow in Casson fluid flow between par-allel plates with Hall current effects were obtained by Shah et al. [14] using an optimal and numerical method. The results showed that Hall currents decrease conductivity, which then increase the velocity and temperature of the fluid. Nadeem and Khan [15] obtained dual solutions of inclined stagnation point nanofluid flow with MHD over an oscillatory shrinking/stretching sheet using the fourthorder Runge-Kutta method. The results obtained indicated that dual solutions occur in both cases that is the shrinking and stretching cases. Furthermore, the obtained lower solution branch shows the same behaviour for the coefficient of skin friction in the shrinking case. On the contrary, the solution in the upper branch has perfect behaviour in both the shrinking and stretching cases. More recent studies of nanofluid flow include among others Kamal et al. [16], Besthapu et al. [17], Sadiq et al. [18], Nayak et al. [19], Yousif et al. [20], Salawu and Ogunseye [21], Kumar and Srinivas [22], Sreedevi and Reddy [23,24], Gireesha et al. [25,26], Archana et al. [27], and Awais et al. [28].
To the best of the authors' knowledge, the study of laminar boundary layer nanofluid flow along a fixed or moving cylinder with heat generation, radiation parameter, and slip parameter has not been widely studied. This work is aimed at studying the effects of thermal radiation, heat generation, slip effects, and other important parameters discussed in Section 4. The model of Lin and Shih [6] is modified into a laminar boundary layer nanofluid model. The fluid flow includes the effects of thermophoretic forces, thermal radiation, heat generation, magnetic field, and Brownian motion. We will focus more on the nanoparticle slip boundary condition, thermophoretic force, and Brownian motion.
The solution of the partial differential equations is obtained by the multidomain bivariate spectral quasilinearization method (MD-BSQLM). This method uses more than one domain in its collocation process; it uses the quasilinearization method based on the Newton-Raphson technique previously used by Bellman and Kalaba [29]. The obtained system of equations is used together into multiple subintervals using the spectral method based on the Chebyshev and Lagrange interpolation polynomials. These are used on the linearized nonlinear equations of PDEs independently in both directions of time and space. This approach of the MD-BSQLM gives very accurate solutions. The method yields more accurate solutions when compared to other methods such as bivariate spectral quasilinearization [30], bivariate spectral homotopy analysis method [31], and bivariate spectral relaxation method (BSRM) [32], among others.
This paper mostly focuses on the robustness and accuracy of the MD-BSQLM and in finding solutions of problems of high complexity. The robustness and accuracy of the numerical method are obtained by analysis and computations of solutions of the momentum, heat, and mass transfer systems of equations. Obtained solutions are compared to those found in the literature, and an acceptable agreement is observed. The effect of different physical parameters on the temperature, velocity, and concentration fields is discussed in tabular and graphical forms in Section 4. The results show the different behaviours that occur when these parameters are changed. The paper will be arranged as follows: the mathematical formulation is in Section 2, the method of solution is in Section 3, and the results will be discussed in Section 4 and conclusions in Section 5.

Mathematical Formulation of the Problem
An axisymmetric, steady boundary layer flow of an incompressible viscous fluid along with a cylinder in motion as shown in Figure 1 is considered. The cylinder is moving with a constant U w in a U ∞ . The T w , T ∞ , C w , and C ∞ are considered. The model will consider the motion of the cylinder in both directions. The mathematical model is presented as subject to the boundary conditions To solve this problem, we introduce the following nondimensional groups and variables: where ψ is the stream function which is defined as which satisfies the continuity equation (1). Substituting equation (9) in equations (2)- (5) gives subject to the boundary conditions: where prime defines derivatives with respect to η, magnetic parameter is M = ðσB 2 0 xÞ/ðρuÞ, Prandtl number is Pr = ðρc p Þ/ðkÞ, thermal radiation parameter is N r = ð16T 3 ∞ σ * Þ/ð3k * kÞ, Eckert number is Ec = ðU 2 Þ/ð4c p Ax m Þ, Ec = ð U 2 Þ/ð4c p ðT w − T ∞ ÞÞ, nondimensional heat generation parameter is He = ð4Q 0 xÞ/ðρc p UÞ, Brownian motion parameter is N b = ðτD B ðC w − C ∞ ÞÞ/ðνÞ, thermophoresis parameter is N t = ðτD T ðT w − T ∞ ÞÞ/ðνT ∞ Þ, Schmidt number is Sc = ðνÞ/ðD B Þ, and nondimensional velocity slip parameter is δ = ðβ 1 /2Þðνu/xÞ 1/2 .

Multidomain Bivariate Spectral Quasilinearization Method
In this section, we describe how to apply the multidomain bivariate spectral quasilinearization method (MD-BSQLM) for solving the system of coupled PDEs that are highly coupled as shown in equations (13)- (15). The solution is decomposed into a large interval with small subdivisions. The solution is evaluated at each time step at the end of each subinterval. The multidomain approach is applied in the ξ direction. The system of equations is first linearized using the quasilinearization (QLM) as indicated in [29,33]. The QLM approach makes use of the Taylor series approximation to linearize the system of equations. In this method, we assume that the value of the function on the current and previous iterations is small. Applying the QLM on (13)-(15) gives the following: where a 0,s = 1 + ξη, Now, let ξ ∈ Ω, where Ω ∈ ½0, T and the domain Ω is decomposed into p nonoverlapping intervals as The PDEs are solved separately at each of the subintervals p. Solutions of the mth interval are obtained by using solutions from the m − 1th interval. This process is implemented as follows: Abstract and Applied Analysis the boundary conditions are imposed An appropriate initial guess is imposed to start the multidomain iteration. This gives the solution in the first interval and satisfies the boundary condition (16). The initial conditions at the next interval are prescribed by the conditions that ensure continuity as follows: The domains in the η and ξ directions have transformed the domain ðx, tÞ ∈ ½−1, 1 × ½−1, 1 for computation, at each subinterval using the linear transformation where Lx is a number that approximates infinity in η. The Chebyshev-Gauss-Lobatto nodes defined in [34,35] [36] for different values of λ and ξ when M = 0, Pr = 0:72, and N r = δ = He = m = r = 0 and in the absence of N t , N b , and Sc.
where ðN x + 1Þ and ðN t + 1Þ are all the points in ηand ξ -directions used for collocation, respectively.
If the solutions f , θ, and ϕ are approximated at each subinterval by using the bivariate Lagrange interpolation poly-nomial of the form where the Lagrange cardinal polynomial functions L p ðxÞ and L q ðtÞ are defined as At the Chebyshev-Gauss-Lobatto points ðx i , t j Þ for i = 0 , 1, 2, ⋯, N x , the derivative of f , θ, and ϕ with respect to η is evaluated as

Abstract and Applied Analysis
where T denotes transpose of the matrix. The derivative of order n for f , θ, and ϕ w.r.t η are approximated as follows: The derivatives of f , θ, and ϕ are evaluated at each point ðx i , t j Þ for j = 0, 1, 2, ⋯, N t as The solution at j = N t of each subinterval is obtained from the solution at the previous level, redenoting i and j indices; equations (42)-(44) can be written as

Abstract and Applied Analysis
For brevity, equations (45)-(47) can be written as

Results and Discussion
A numerical solution for the boundary layer nanofluid flow along a moving cylinder with thermal radiation, heat generation, and slip effects is studied. The solutions for the system of partial differential equations are presented. The bivariate spectral quasilinearization method (MD-BSQLM) is used to solve the coupled system. To validate the numerical method, a comparison with results in [36] was made, and the results are in acceptable agreement. The finite difference-based method used in [36] shows a slightly more rigorous approach than the one applied in this study. We present important and interesting results only. The results obtained are presented in Table 1 and  14 Abstract and Applied Analysis and ∂f /∂ξ = ∂f ′/∂ξ = 0 give the set of equations (13)-(15) to those of [36]. Figures 2 and 3 show the effect of varying the velocity slip parameter δ on the velocity and temperature profiles. Increasing δ result in increasing both the velocity and temperature profiles. Increasing the velocity slip parameter causes the fluid layer in contact with the cylinder surface to slip thereby causing the velocity at the surface to increase. The velocity decrease with increasing boundary layer. When fluid particles move faster on the cylinder surface, this causes more fluid particle interaction causing the temperature to increase with increasing δ.
Figures 4 and 5 depict the effect of changing the magnetic parameter M on velocity profiles and changing the heat generation parameter on the temperature profiles. Increasing the magnetic parameter result in the decrease in velocity profiles, this result agrees with many other results displayed in the literature in [20,37] and [21]. If the heat generation parameter is increased, this then results in the increase in the temperature profiles; this result has been shown in many other different studies. This shows a strong agreement in the results when the MD-BSQLM is used. Figure 6 shows the effect of increasing the radiation parameter N r on the temperature profiles. If the radiation parameter is increased, this has an effect of increasing temperature profiles. A point of "inflexion" is noticed close to the surface of the cylinder. When radiation is applied, it causes a rise in temperature causing the mass particles close to the surface with lower temperature to rise. They then become heated as they move further from the surface causing the temperature to increase further away from this point. Figures 7 and 8 show the effect of varying the Eckert number Ec and the radial coordinate r. Increasing the Eckert number results in the increase in temperature profiles. The trends in the increase in temperature fields are more pro-nounced at close to the cylinder surface. This is caused by the presence of velocity slip and radiation which increase temperature at the surface. The gradient of the velocity profile at the surface −θ′ð0Þ (heat transfer coefficient) can be seen increasing from negative to positive. Increasing the radial coordinate shows a decrease in the concentration profiles. This shows that as we move outward from the cylinder surface, the mass particle concentration decreases. This is consistent with the general settling of mass particles to the cylinder surface due to gravity. The settling is less enhanced due to the fluid flow which tends to carry the particles along with the flow regime. Figures 9 and 10 show the effect of varying the thermophoresis parameter N t and the the Brownian motion parameter N b and on the concentration profiles. If the Brownian motion parameter is increased, this results in decreasing the concentration profile at close to the surface. A reverse effect is noticed further away from the cylinder surface. Increasing the Brownian motion parameter causes the mass particles to move more rapidly in random motion this destabilizes the particle motion. A more enhanced point of inflexion is noticed. Increasing the thermophoresis parameter N t results in decreasing the concentration profiles at close to the cylinder surface.
In Figure 11, it is shown that if the Schmidt number Sc is increased, this result in the decrease in the concentration profiles. The result is consistent with those in the literature.
The effect of changing parameters of the coupled partial differential equations affects the stability of the numerical solution. Sometimes when parameters exceed certain ranges, the numerical method fails. In this section, we study the effect of changing various parameters on the stability of the numerical solution. This is shown in the graphs from  15 Abstract and Applied Analysis N b on the residual error of f , θ, and ϕ, and N r on f , θ, and ϕ. Figure 12 shows the effect of increasing the velocity slip δ on the residual error of the solution of f ; the error does not deteriorate with increasing iterations. The sudden drop in residual error in less than five iterations shows the accuracy of the MD BSQLM. Figure 13 shows the effect of increasing the Brownian motion parameter N b on the residual error of f ; the residual error becomes very small with larger values of the Brownian motion parameter. The error becomes more stable with increasing iterations. In Figures 14 and 15, the errors E ϕ and E θ improve with increasing Brownian motion parameter. At low values of the parameter convergence is delayed, happens after five iterations. In both cases, the error does not deteriorate with increasing iterations.
In Figures 16-18, it is shown that increasing the radiation parameter does not make the errors E f , E ϕ , and E θ deteriorate with increasing iterations. Increasing the radiation parameter N r result in improving the stability of the errors of the MD-BSQLM.

Conclusion
In this paper, we have studied the numerical solution of laminar boundary layer flow along with a moving cylinder with thermal radiation, heat generation, and surface slip effect. The coupled partial differential equations of the mass and heat transfer along with a cylinder in motion are solved by the multidomain bivariate spectral quasilinearization method (MD BSQLM). The most important results show that increasing the velocity slip factor results in an enhanced increase in velocity and temperature profiles. Increasing the heat generation parameter increases temperature profiles; if the radiation parameter and the Eckert number are both increased, this would then result in the increase temperature trends. The concentration fields decrease with increasing radial coordinate. Increasing the thermophoresis parameter and Brownian motion both destabilizes the concentration profiles. Increasing the Schmidt number reduce temperature fields. The effect of increasing selected parameters: the velocity slip, Brownian motion, and radiation parameter on all residual errors show that these errors do not deteriorate. This shows that the multidomain bivariate quasilinearization method (MD BSQLM) is very accurate and robust. The method is also easier to implement than the traditional finite difference methods as shown in Section 3. The method also shows very low residual errors only after five iterations. The method adds value to the way in which differential equations arising from fluid dynamics are solved.