Finite Element Method for Non-Newtonian Radiative Maxwell Nanoﬂuid Flow under the Inﬂuence of Heat and Mass Transfer

: The recent study was concerned with employing the ﬁnite element method for heat and mass transfer of MHD Maxwell nanoﬂuid ﬂow over the stretching sheet under the effects of radiations and chemical reactions. Moreover, the effects of viscous dissipation and porous plate were considered. The mathematical model of the ﬂow was described in the form of a set of partial differential equations (PDEs). Further, these PDEs were transformed into a set of nonlinear ordinary differential equations (ODEs) using similarity transformations. Rather than analytical integrations, numerical integration was used to compute integrals obtained by applying the ﬁnite element method. The mesh-free analysis and comparison of the ﬁnite element method with the ﬁnite difference method are also provided to justify the calculated results. The effect of different parameters on velocity, temperature and concentration proﬁle is shown in graphs, and numerical values for physical quantities of interest are also given in a tabular form. In addition, simulations were carried out by employing software that applies the ﬁnite element method for solving PDEs. The calculated results are also portrayed in graphs with varying sheet velocities. The results show that the second-order ﬁnite difference method is more accurate than the ﬁnite element method with linear interpolation polynomial. However, the ﬁnite element method requires less number of iterations than the ﬁnite difference method in a considered particular case. We had high hopes that this work would act as a roadmap for future researchers entrusted with resolving outstanding challenges in the realm of enclosures utilized in industry and engineering.


Introduction
The extensive manufacturing and industrial applications have made studying boundary layer non-Newtonian fluid much more substantial for the research scholars. Several processes such as reparating plastic polymers, hot rolling, cooling metallic plates, drilling muds and assembling optical fibers have followed the principles of boundary layer non-Newtonian fluid. Scientists have suggested integral and rate models for boundary layer studies as a single model is not enough to cover the versatility of properties. The given study comprises the Maxwell fluid model, which is the subclass of the rate type model and is used to determine the effect of relaxation time.
The most stimulating area in applied sciences is the non-Newtonian fluid flow, which has attracted many scholars due to its extensive medical and engineering applications (Rivlin and Ericksen [1]). This fluid exhibits several alternating properties such as variable Details of MHD flow of Carreau fluid based on Cattaneo-Christov flux theory can easily be spotted in [31,32].
A new focus on mass heat flow has brought experts together, and they discovered that it could be used in a wide range of fields such as air conditioning, electronic device cooling and nuclear reactor cooling. It can also be used to desalinate water and in the food and pharmaceutical industries. Initially, the Fourier [33] and Fick [34] laws were used for interpreting mass and heat flow best, but later on, scientists observed some limitations of these laws, which could affect the studies. One of the most scavenging drawbacks is forming a parabolic type of equation. Javed and Nadeem [35] utilized two concentric cylinders to observe Casson fluid's mass and energy flow. Double stratified second-grade fluid flow was studied by Mallawi and colleagues [36] using the Riga plate heat flux model and thermal radiation. A review of the literature [37,38] revealed the focus of various scholars on heat transfer of turbulent nanofluid flow to minimize the energy consumption in a solar collector.
The study of Darcy-Forchheimer flow of Reiner-Philippoff nanofluid over the stretching sheet was investigated in [39] with the involvement of Motile microorganisms. The effects of heat source/sink and melting phenomenon were also considered. The Matlab solver was considered to solve the obtained dimensionless equations. The results confirmed that rising thermophoresis and Schmidt number values enhanced the heat transfer coefficient. By changing thermophysical characteristics, entropy reduction in the thermos and non-Newtonian nanofluid models was addressed [40]. Two types of nanofluids, namely Copper-Engine Oil and Zinconium Dioxide-Engine Oil, were considered for the study. The effects of different parameters on velocity, temperature and entropy distribution are shown in graphs. The Keller-Box scheme was implemented to solve differential equations obtained from the fluid phenomenon. A mixture model was considered to simulate a rotating tube bundle [41]. The microchannel heat exchanger's performance was assessed by applying various operating factors such as Reynolds number, rotation speed and concentration of the nanofluid. Instead of finding numerical solutions, exact expressions for dimensionless velocity and temperature were obtained in [42]. It was pointed out that the Newtonian fluid flow was faster than Maxwell fluid flow.
Since the finite element method is one of the numerical methods that can be used to find the solutions to linear and nonlinear ordinary and partial differential equations, the recent approach to applying the method is based on Galerkin weighted residuals, and integrals in this approach are computed using numerical Gauss-Legendre three points formula integrations. The model for non-Newtonian fluid flow over stretching sheet is given, and converted ODEs were solved by applying three numerical approaches. The second applied method is the finite difference method, which is second-order accurate at all the grid points except the last one. The numerical experiments conclude that the finite difference method (FDM) is more accurate than the finite element method (FEM) with linear interpolation polynomial for the considered problem.
Moreover, a Matlab solver bvp4c was also employed for solving nonlinear dimensionless ODEs. The solver can be used to find solutions to problems in finite domains. Even though the domain of the considered problem is infinite, a finite domain is chosen for numerical purposes. The solver, in most cases, converges, but it may give an error. As a result, an additional loop is utilized to make the procedure workable for some or all of the cases of the problem under consideration.

Problem Formulation
Consider incompressible, laminar, steady, electrically conductive two-dimensional non-Newtonian Maxwell nanofluid flow over a stretching sheet. Let the plate be stretched with stretching velocity U W . Let x-axis be taken along the sheet and y-axis be taken perpendicular to the sheet. The flow is generated by the sudden movement of the sheet towards the positive x-axis. The fluid is embedded in a porous medium, and the permeability of this porous medium is denoted by k. The magnetic field has strength B 0 applied perpendicular to the sheet. Let u and v be the velocity components in x and y directions, respectively. The fluid temperature is denoted by T, and concentration is denoted by C. Let T w and C w be temperature and concentration at the plate, whereas T ∞ and C ∞ is the temperature and concentration away from the plate. Under the boundary layer assumption(s), the governing equations of this phenomenon can be expressed as: Subject to the boundary conditions where ν represents the kinematic velocity, λ 0 represents the relaxation time, σ denotes an electrical conductivity, ρ denotes the density of the fluid, α represents the thermal diffusivity, τ denotes the ratio of heat capacities, D B denotes Brownian motion coefficients, D T denotes thermophoretic diffusion coefficient, µ is the dynamic viscosity, c p is the specific heat capacity, k 1 is the ratio of reaction and q r is the Rosseland radiative heat flux. The linearized Rosseland radiative heat flux is expressed as: By considering the flux q r , the corresponding term in energy Equation (3) was rewritten as: where σ * denotes the Stefan-Boltzmann constant and k * denotes mean absorption coefficient. For reducing governing Equations (1)-(5) into dimensionless forms, the following transformations are considered: Under the mentioned transformations, Equations (1)-(5) can be expressed into dimensionless ODEs as where Deborah number λ, porosity parameter k p , magnetic parameter M, Prandtl number P r , Brownian motion parameter N b , thermophoretic parameter N t , radiation parameter R d , Eckert number E c , Schmidt number S c and reaction rate parameter γ are denoted as The skin friction coefficient, Local Nusselt and Local Sherwood are defined as: The dimensionless form of skin friction coefficient, Local Nusselt and Sherwood numbers can be expressed as: where R e x = U w x ν .

Finite Element Method
Equations (8)- (11) are solved by employing the finite element method. For doing so, let the whole domain [0, η ∞ ] be divided into a finite number of sub-domains called elements. Each element is a line segment containing two nodes at both ends of each element. One node of two adjacent elements in common. The linear polynomial is interpolated in each element. Let the length of each element be h. The ith element is comprised of two nodes at η = η i and η = η i+1 . Linear interpolating polynomial for f , θ and φ in ith element can be expressed as: By using two points at η = η i and η = η i+1 , Equations (18)- (20) can be expressed as: The left node of ith element in the domain is assigned with nodal variables f i , θ i and φ i with the nodal coordinate value η i . Similarly, the right node of ith element is assigned with corresponding nodal variables and nodal coordinates. The variables f , θ and φ appeared in Equations (18)- (20) and are called trail functions, whereas ψ 1 and ψ 2 are called shape functions. The shape function satisfies two properties. The first property of the shape function is to satisfy the following conditions: where ψ 1 has value 1 at the left node of ith element and 0 at the right node also ψ 2 have values 0 and 1 at the left and right node of ith element, respectively. This first property is useful for obtaining continuous solutions over the whole domain.
The second property states that the sum of the shape functions for a single element is unity, i.e., The advantage of this property is to ensure the constant solution within each element, provided that the solution is the same at each node of that element.
For finding a numerical solution using Galerkin finite element method, weighted residuals of Equations (8)-(10) were constructed. Since Equation (8) is third-order ODE, to apply the method, it is converted into a system of two equations: The weighted residual of Equations (26), (27), (9) and (10) are given as: Since formulations (29)-(31) are strong formulations, one of the disadvantages of using strong forms is that these have vanished when linear interpolation is considered for each element. Thus, in this manner, weak formulations are constructed by integrating the second-order terms in Equations (29)-(31) as: Integrals Equations (32)- (35) are constructed on the whole domain, and these integrals are constructed on ith element as in the forms From Equations (36)-(39), the following matrix-vector equations can be achieved Let K i be the coefficient matrix in Equation (40) corresponding to ith element in the domain. Each entry in K i is a matrix of order 2 × 2, and each [R s ] is a vector of length 2 for s = 1, 2, 3, 4. The entries of the matrix K i are calculated from the following integrals and remaining matrices are zero matrices. Where the variables have " − " notation is kept fixed. The bar notation terms are treated as and ψ 1 and ψ 2 are shape functions or components of the test function, ψ 1 and ψ 2 denote derivatives of the shape function with respect to independent variable η. A numerical integration based on Gauss Legendre Quadrature three-point formula is adopted to carry out fast integration. Points in this numerical integration are the roots of third-degree polynomial and weights in this integration 5 9 , 8 9 and 5 9 . In order to find more accurate derivatives of solutions, a modified approach of the finite element method proposed in [43] was considered. This approach employs different formulas for calculating skin friction coefficients, Local Nusselt and Sherwood numbers. The discretization of these quantities is given as where 1 f , 1 f 1 , 1 θ and 1 φ denote the values of f , f 1 , θ and φ at the first grid point. The standard or classical finite element method using first-degree polynomial provides firstorder accuracy of derivatives of solutions. In contrast, the modified approach given in [43] provides at least second-order accuracy for derivatives of solutions subject to the condition that the solution is at least second-order accurate.

Validations
Since Equations (8)-(10) are nonlinear ordinary differential equations, to handle nonlinear differential equations, an iterative procedure was adopted by fixing some quantities. If the desired stopping criteria are fulfilled, iterations are stopped. The accuracy is checked by finding the difference of solutions of Equations (9), (10), (26) and (27) at two consecutive iterations and at each grid point. The stopping criteria can be expressed as where is a small number near zero and F k s,i denotes one of three dependent variables f , θ and φ at some iteration, say k 2 and at grid point i. The mesh-free study was also carried out by constructing Table 1, which shows that the results are independent of h. Table 1 shows the mesh-free study for the finite element method and Matlab solver bvp4c. The different numbers of elements are considered with results computed over variations in the number of elements, and changes in dependent variables can be seen in Table 1. Table 1. Meshfree analysis for f , θ and φ with varying amounts of elements using = 0.9, P r = 1.5, For the verification of Matlab code and computed results from the modified finite element method, a comparison was also made with the finite difference method. For applying the classical approach of the finite difference method for Equations (9), (10), (26) and (27), difference or discretization equations can be expressed in the forms:

No. of Elements/Nodes
Equations (65)-(68) were solved by employing an iterative method. The iterative method needs one initial guess and stops when the criteria are satisfied. Table 2 shows the comparison of the classical and modified finite element method given in [43] and the classical finite difference method for computing − f (0) with those results given in the past research [44,45]. The modified finite element method uses second-order difference formula for computing − f (0). The standard or classical finite element method with linear polynomial interpolation uses the first-order difference formula.

Results and Discussions
The finite element method employed in this study for solving ordinary differential equations is based on the Galerkin weighted residual method using weak formulations for equations with order two. The linear elements were chosen with equal lengths. Any numerical integration can be adopted for the computations of integrals due to quick calculations. The considered finite element method uses two nodes for each element. Each element produces a matrix-vector equation, and the assembly of each matrix equation corresponding to each element gives a global matrix called a stiffness matrix. By comparing the finite element method (F.E.M) with the finite difference method (F.D.M), it was concluded that F.E.M consumes less time due to utilizing fewer iterations than those utilized by F.D.M. One of the reasons behind this is the use of an iterative procedure with F.D.M. that slows the whole procedure, and the solutions' convergence was obtained by utilizing many iterations. The F.E.M uses the linearized forms of nonlinear equations, and the computations of results are also performed by solving matrix-vector equations from the Matlab solver. The Matlab solver finds the exact solutions for matrix-vector equations, whereas the F.D.M finds the solutions of nonlinear differential equations iteratively without linearizing nonlinear equations. Figure 1 shows the geometry of the considered problem. Figures 2-13 are drawn by using the standard finite element method for solving differential Equations (9), (10), (26) and (27) with boundary conditions (11). Figure 2 deliberates the impact of the magnetic parameter on the velocity profile. Figure 2 shows that velocity decays by increasing magnetic parameters. This growth of velocity profile is the consequence of resistance of Lorenz force that enhances by increment in the magnetic parameter. The impact of the Deborah number λ on the velocity profile is portrayed in Figure 3. The Deborah number contains the relaxation time parameter, which states the time when the fluid obtains an equilibrium state by applying stress. Larger values of the Deborah number give larger relaxation time, which gives more resistivity of stress in the fluid or enhances the viscosity of the fluid and, consequently, a slower velocity profile is obtained. The variation in porosity parameter on velocity profile is deliberated in Figure 4. The velocity profile de-escalates by enhancing the porosity parameter because the porosity parameter enhances the resistivity in the flow, leading to decays in the velocity profile. Figure 5 shows the variation in the Prandtl number on the temperature profile. The velocity profile de-escalates by rising values of the Prandtl number. The higher Prandtl number is responsible for smaller thermal diffusivity because both have an inverse relationship. The slower thermal diffusivity reduces thermal conductivity, which yields a slower temperature profile. Figure 6 shows the impact of thermophoretic parameters on the temperature profile. The temperature profile escalates by rising values of the thermophoretic parameter. This is due to the thermophoresis phenomenon that tends to move the hot particles from the immediate vicinity of the plate to its surroundings, and consequently, the temperature rises. Figure 7 deliberates the effect of the Brownian motion parameter on the temperature profile. The temperature profile grows by enhancing the Brownian motion parameter because the escalation in the Brownian motion parameter is responsible for spreading the hotter particles due to the increase in random movement of particles, which results in a temperature rise. The effect of the Eckert number on temperature profile is portrayed in Figure 8. The temperature profile rises by enhancing the values of the Eckert number. The increment in the Eckert number gives rise to the temperature difference between wall and ambient temperature, and therefore temperature rises. The effect of the radiation parameter on the temperature profile is shown in Figure 9. The temperature rises by enhancing the radiation parameter because the energy of the flow increases by incoming radiation. The impact of the Schmidt number on the concentration profile is portrayed in Figure 10. The concentration profile gets down by boosting the values of the Schmidt number. The de-escalation in the concentration profile is the decay of mass diffusivity due to an increase in Schmidt number. Figure 11 shows the thermophoretic parameter on the concentration profile. Rising values of thermophoretic parameters boost the concentration profile. This happened due to the thermophoresis phenomenon that tends to move nanoparticles from the immediate vicinity of the plate to its surroundings, and thus concentration profile grows. Figure 12 portrays the variation in the Brownian motion parameter on the concentration profile. The concentration profile decays by enhancing the Brownian motion parameter. The boosting up of the Brownian motion parameter gives the faster random movement of nanoparticles, resulting in decays in the concentration profile. The effect of the reaction rate parameter on the concentration profile is shown in Figure 13. From this Figure 13, it can be seen that the concentration profile de-escalates by rising values of the reaction rate parameter because higher values of the reaction rate parameter create impurity in the fluid responsible for decaying the concentration profile.                     grows by enhancing the Deborah number and magnetic parameter whereas decays by de-escalating the Brownian motion parameter and thermophoretic parameter.   x Nu x decays by de-escalating the Brownian motion parameter and thermophoretic parameter. Table 3. Numerical values and comparisons of finite element methods and finite difference methods using P r = 1.5, = S c , E c = 0.4, R d = 0.7, k p = 0.5, γ = 0.9.  Simulations were also carried out by employing a finite element method for solving governing equations of heat transfer of Newtonian fluid under the effects of radiations over moving surface towards x-axis. The geometry consists of a four-sided rectangular shape in which two sides are used as input and output, one side is a moving boundary and the rest have a no-free boundary condition. Figures 14-16 were drawn by solving partial differential equations with different velocities of the bottom surface. The momentum and thermal boundary layers can be seen in Figures 14-16.

Conclusions
The finite element method was utilized for heat and mass transfer of non-Newtonian fluid flow over the sheet. The linear interpolation polynomials were considered in the application procedure of the finite element method. The results were also compared with those given in past research and the classical finite difference method. In addition, the simulations for heat transfer of Newtonian fluid flow model over stretching sheet were also provided by using software that implements finite element method. By using the considered iterative method, the classical finite difference method consumed more iterations, and it was computationally expensive but gave better accuracy. The concluded points can be expressed as: • Velocity profile was de-escalated by escalating magnetic parameter, Deborah number and porosity parameter; • The temperature profile was grown by rising values of radiation parameters, Brownian motion and thermophoretic parameters; • The concentration profile was escalated by enhancing the Brownian motion and thermophoretic parameters and de-escalating by rising Schmidt number and reaction rate parameters.
Further, the finite element method considered in this work can be employed to solve efficiently nonlinear problems of similar type that arise in computational fluid dynamics with some extra effects. Following the completion of this work, it will be possible to propose other applications for the currently employed methodology, if desired [46][47][48][49]. In addition, the developed method is easy to use and can solve a wider range of differential equations in both practice and theory.

Data Availability Statement:
The manuscript included all required data and implementing information.