Abstract
This paper demonstrates the applicability of the large parameter spectral perturbation method (LSPM) to a coupled system of partial differential equations that cannot be solved exactly. The LSPM is a numerical method that employs the Chebyshev spectral collocation method in the solution of a sequence of ordinary differential equations (ODEs) that are derived from decomposing coupled systems of nonlinear partial differential equations (PDEs) using series expansion about a large parameter. The validity of the LSPM is applied to the problem of boundary layer flow and heat transfer in a micropolar fluid past a permeable flat plate in the presence of heat generation and thermal radiation. The coupled nature of the PDEs that define the problem under investigation precludes the option of using series-based methods that seek to generate analytical solutions even in the presence of small or large parameters. The present study demonstrates that the LSPM can easily overcome this limitation while giving very accurate results in a computationally efficient manner. For qualitative validation of the results and the numerical method used, calculations were carried out to graphically obtain the velocity, microrotation, and temperature profiles for selected physical parameter values. The results obtained were found to correlate with the results from a published literature. For quantitative confirmation of our findings, the LSPM numerical solutions were again validated against known results from the literature and against results obtained using the multidomain bivariate spectral quasilinearisation method (MD-BSQLM), and the results were observed to be in perfect agreement. Further accuracy validation is displayed by using residual error and solution error analysis on the governing PDEs and their underlying solutions. This study’s findings indicate that the heat generation and thermal radiation parameters have related effects on the temperature profile, enhancing both the fluid temperature and the thermal boundary layer thickness.
1. Introduction
The theory of micropolar fluids, as first introduced by Eringen [1, 2], can be used to explain fluids with microstructure and fluids which contain a suspension of rigid bodies such as dirty water, liquid crystals, colloidal fluids, paint, polymers, and blood. Physically, micropolar fluids represent non-Newtonian fluids containing randomly oriented molecules suspended in a viscous medium. An extensive collection of articles on this subject and its applications can be found in the books by Eringen [3, 4] and Lukaszewicz [5] and articles by Ariman et al. [6, 7]. Micropolar fluids have recently received the attention of many researchers because of their essential applications in some processes in the industry. Such applications include cooling metallic plates in a bath, extrusion of polymer fluids, solidification of liquid crystals, exotic lubricants, suspension solutions, magnetic fluids, and dust clouds, among others. Inspired by the prior work of Eringen [1, 2], several researchers have investigated the characteristics of micropolar fluids under different conditions. These include the investigation of Kumari and Nath [8] who used the quasilinear finite difference scheme to gain numerical solutions to the equations governing the flow of unsteady incompressible boundary-layer flow of a micropolar fluid near a stagnation point. Rees and Andrew [9] investigated the Blasius boundary layer flow of a micropolar fluid over a flat plate. The MHD boundary-layer flow of a micropolar fluid past a wedge with variable wall temperature has been studied by Ishak et al. [10]. Nazar et al. [11] investigated the unsteady boundary layer flow over a stretching sheet in a micropolar fluid. In the study of [9–11], numerical solutions were generated using a finite-difference scheme known as the Keller-box method. Lok et al. [12] studied the unsteady mixed convection flow of a micropolar fluid near the stagnation point on a vertical surface using the Keller-box method. Their study discovered that the skin friction increases with time and a smooth transition occurs from the small-time solution to the large-time solution. Mixed convection flow of micropolar fluid over an isothermal plate with variable spin gradient viscosity was studied by Hossain and Chowdhury [13] using an implicit finite difference scheme. Dual solutions for unsteady mixed convection flow of MHD micropolar fluid over a stretching/shrinking sheet with a nonuniform heat source/sink were presented by Sandeep and Sulochana [14]. Their study generated numerical solutions to the governing equations using the shooting method technique. The shooting method with the fifth-order Runge-Kutta-Fehlberg integration scheme was utilized by El-Aziz [15] to examine the influence of viscous dissipation on mixed convection of a micropolar fluid from an unsteady stretching surface. Rashidi et al. [16] generated approximate analytic solutions for the heat of a micropolar fluid through a porous medium with radiation using the homotopy analysis method. Hayat et al. [17] carried out a numerical study using the homotopy analysis to investigate the mixed convection flow of a micropolar fluid with radiation and chemical reaction. Hossain et al. [18] used the series method to investigate the boundary layer flow and heat transfer in a micropolar fluid past a porous flat plate. Other micropolar fluid flow studies have been conducted by, among others, Shit et al. [19], Govardhan and Kishan [20], Ahmad et al. [21], Bhattacharyya et al. [22], El-Dabe et al. [23], Singh et al. [24–29], Mabood et al. [30], Khan et al. [31], and Pattnaik et al. [32].
The mathematical equations modeling boundary layer and heat transfer flow of a micropolar fluid past a permeable flat plate with heat generation and thermal radiation effects are defined by a coupled system of three nonlinear partial differential equations. The problem investigated in this study extends the analysis of Hossain et al. [18] by putting into account the effects of heat generation and thermal radiation. The problem without the effects of heat generation and thermal radiation has earlier been solved by Hossain et al. [18] analytically using the perturbation method and numerically using an implicit finite difference method. The perturbation method and implicit finite method have been used extensively by many researchers on related problems defined on large parameter intervals (see, for example, [18, 33, 34]). The perturbation method is an analytic technique for finding approximate solutions to differential equations. The observation made in earlier studies that used the perturbation method to obtain series solutions about a large perturbation parameter is that higher-order perturbation equations may be impossible or difficult to solve analytically beyond a specific order of approximation. This leads to less accurate results if one, two, or three perturbation approximations are used to generate approximate solutions. Hence, there is a need to obtain higher-order perturbation approximations. Besides, there are limits to how far the perturbation series analytic approach can be used in nonlinear systems of partial differential equations involving many coupled equations. This is because nonlinear systems of partial differential equations involving many coupled equations are difficult to solve exactly. The implicit finite difference method is a known numerical method used by many researchers for solving nonlinear PDEs defined on large parameter intervals. Cebeci and Bradshaw introduced the technique [35]. The finite difference schemes require many grid points to achieve reasonable accuracy; hence, a lot of computer memory and computational time may be needed to obtain accurate results.
This study is aimed at demonstrating the applicability of the large parameter spectral perturbation method (LSPM) on coupled systems of nonlinear partial differential equations defined over a large parameter interval, which cannot be solved analytically even with methods that look for series solutions. To the best of the authors’ knowledge, series approaches similar to the LSPM have not been reported in the literature for solving coupled nonlinear PDEs describing micropolar fluid flow. In this work, we apply, for the first time, the LSPM on coupled systems of nonlinear PDE modeling boundary layer flow and heat transfer in a micropolar fluid past a porous flat plate with heat generation and thermal radiation effects. Also, this current study is a numerical investigation of boundary layer flow and heat transfer in a micropolar fluid past a porous flat plate with heat generation and thermal radiation effects. The LSPM combines both analytical and numerical solution techniques. In the LSPM, series expansions are constructed about a large perturbation parameter, and the Chebyshev spectral collocation method is used to solve the resulting series equations numerically. In this present study, we observe that the sequence of ordinary differential equations (ODEs) derived from the series expansion is coupled and cannot be solved exactly, hence the introduction of the Chebyshev spectral collocation method to resolve the resulting series equations numerically. With the use of the Chebyshev spectral collocation method, it becomes possible to gain the approximate numerical solution of the higher-order perturbation equations, which are impossible to obtain analytically due to the nature of the problem investigated in this study. The LSPM approximate solutions are validated using the multidomain bivariate spectral quasilinearisation method (MD-BSQLM) and published results in the literature. Our results are found to be in excellent agreement. The MD-BSQLM is a numerical method that blends the quasilinearisation idea with the Chebyshev spectral collocation and bivariate Lagrange interpolation. In the MD-BSQLM, the nonlinear systems of PDEs are linearized using the quasilinearization method of Bellman and Kalaba [36]. The resulting equations are then integrated into multiple subintervals using the Chebyshev spectral collocation method. The numerical results from this study show that the LSPM is accurate and can be utilized as an alternative numerical tool for solving coupled nonlinear systems of partial differential equations defined over large parameter domain. The advantage of the LSPM is that, unlike other numerical methods, it applies discretization only in the space -direction when solving a partial differential equation. This feature enables the LSPM to compute approximate numerical solutions in a very efficient and computationally fast manner.
2. Mathematical Formulation
We consider a steady, two-dimensional laminar boundary layer flow of a micropolar fluid along a permeable vertical flat surface in the presence of heat generation and thermal radiation. The external flow comprises a uniform free stream with temperature and velocity . The ambient fluid temperature is assumed to be uniform at and that of the surface at . The flow configuration and coordinate system are displayed in Figure 1. Under the usual boundary layer approximation, the flow is governed by the following equations (see [18]). We consider a steady, two-dimensional laminar boundary layer flow of a micropolar fluid along a permeable vertical flat surface in the presence of heat generation and thermal radiation. The external flow comprises a uniform free stream with temperature and velocity . The ambient fluid temperature is assumed to be uniform at and that of the surface at . The flow configuration and coordinate system are displayed in Figure 1. Under the usual boundary layer approximation, the flow is governed by the following equations (see [18]):
In the above equations, and are the velocity components along - and -axis, is the vortex viscosity parameter, is the kinematic viscosity, is the porous medium permeability, is the microrotation component normal to the plane, is the thermal diffusivity, is the density of the fluid, is the specific heat at constant pressure, is the heat generation constant, is the fluid temperature in the boundary layer, and is the radiative heat flux defined using the Rosseland approximation as
The boundary conditions are
In equation (6), , is the transpiration parameter, where is the suction velocity of the fluid through the surface of the plate, is the reference length, and is the Reynolds number. When is positive, it denotes suction or withdrawal and a negative stands for injection or blowing of fluid through the surface of the plate. In this present investigation, the case of the suction will only be considered rather than the blowing case, and therefore, is taken to be positive throughout in this study. Furthermore, is a constant which ranges between . It should be noted that the case indicates the strong concentration of micro-elements, corresponds to the vanishing of anti-symmetric part of the stress tensor and represents weak concentration, and may be utilized for the modeling of turbulent boundary layer flows. The following transformations are then introduced: where is the pseudosimilarity variable, is the nonsimilar parameter, and is the stream function defined as
Substituting the transformations given in equation (7) into equations (2)–(4) gives the following set of partial differential equations:
The boundary conditions (6) become
In the above equations, prime denotes differentiation with respect to , is the porous medium parameter, is the Prandtl number, is the thermal radiation parameter, and is the heat generation parameter. These parameters are mathematically defined as where is the local Darcy number and is the local Reynolds number.
The shear stress, , the couple-stress, , and the rate of heat transfer, , are calculated as
The shear stress, , defines the tangential force per unit area parallel to a surface. The presence of shearing stress in a fluid enables us to observe the fluid’s viscosity. Due to the shearing stress, a viscous fluid produced resistance to the body moving through it and between the particles of the fluid itself. The couple-stress, , of fluids defines the rotation component in terms of the velocity component. The rate of heat transfer, , also known as the Nusselt number, is a nondimensional heat transfer coefficient. It gives a comparison between the conductive and convection heat transfer rates.
3. Large Parameter Spectral Perturbation Method (LSPM)
In this section, we derive the asymptotic solution of equations (9)–(11) subject to the boundary conditions (12)–(14) when the nonsimilar parameter is large. When is large, the dominant terms in equation (9) are and , and in (10), and and in (11). It is sufficient to balance these terms in magnitude. Therefore, given that as , it is necessary to find appropriate scaling for and . Balancing the order of magnitude of and in equation (9), and in (10) and and in (11) using the method of dominant balance, it is found that and , and as . Accordingly, we introduce the following transformations that are valid for large values of :
Substituting equation (17) into equations (9)–(11) reduces the equations to where prime denotes differentiation with respect to . The corresponding boundary conditions are
Since is large, we seek solutions to equations (18)–(20) using the perturbation series approach. The functions , , and are expanded in powers of as
Substituting equation (22) into equations (18)–(20) and then equating the coefficients of like powers of , we obtain the equation for as
subject to the following boundary conditions:
The sequences of equations for are given as
subject to the following boundary conditions:
The initial solution at used to begin the LSPM algorithm can be found by solving equations (23)–(26) subject to the boundary conditions (26). We remark that equations (23) and (24) are interlinked and cannot be solved independently of each other or solved exactly and therefore must be solved numerically as a coupled system. Solving equation (25) analytically, we get the exact solution as
The Chebyshev spectral collocation method is very accurate and convenient to use in solving differential equations of the type (23)–(25) and (29)–(31) which are linear with constant coefficients. For brevity, details of the spectral method have been omitted in this paper. Interested readers can refer to the works of (see [37–42]) on how the spectral collocation method has been used to solve related PDEs. Applying the spectral method on equations (23) and (24) gives where where , with being the Chebyshev spectral differentiation matrix of size whose entries are defined in [43, 44], is a finite value selected to be large to approximate the conditions at infinity, is an identity matrix of size , and is a zero vector of size . The boundary conditions are imposed on the zero vectors, starting from the initial solution from , , and which are known by solving equations (23)–(24) numerically and (25) exactly. Solutions to the higher order approximations , , and for given by equations (29)–(31) can then be obtained using the spectral method. It is for this reason the method is referred to as the large parameter spectral perturbation method. Similarly, applying the spectral collocation method on equations (29)–(31) gives where where , , and are defined as
Hence, starting from already known , , and , the solutions are obtained as where
The shear stress, , the couple-stress, , and rate of heat transfer, , can be calculated as
4. Results and Discussion
In this section, we present the approximate numerical solutions of the governing system of nonlinear partial differential equations (9)–(11), subject to the boundary conditions (12) obtained using the large parameter spectral perturbation method (LSPM). To verify the accuracy the LSPM, the present numerical results are confirmed against the multidomain bivariate spectral quasilinearisation method (MD-BSQLM). In the MD-BSQLM which is a numerical method, we first linearise the governing nonlinear equation using the Newton-Raphson based quasilinearisation approach developed by Bellman and Kalaba [36] and then solve the resulting sequence of linearised equations on multiple intervals using the Chebyshev spectral collocation method. We remark that the multidomain approach was applied only in the direction. From our numerical simulations, it was observed that collocation points in the space variable and collocation points for the nonsimilar variable were sufficient to give accurate results in all the spectral method discretization done. Furthermore, the finite value selected to be large used to approximate the boundary condition at infinity was set to be in the LSPM and in the MD-BSQLM. For qualitative validation of our methods, calculations were carried out to graphically obtain the shear stress, , the couple-stress, , the rate of heat transfer, , and the velocity , microrotation , and temperature profiles for different physical parameters, and our results were found to correlate with results from published literature. To quantitatively verify our results, the LSPM and the MD-BSQLM numerical solutions were also confirmed against published results reported by Hossain et al. [18] for different values of the vortex viscosity parameter, . We remark that results that are comparable with those of Hossain et al. [18] were achieved in the LSPM with while in the MD-BSQLM with at . The comparison is displayed in Table 1, and the results were observed to be in very good agreement. This indicates that the numerical methods used in this study are reliable.
Table 2 shows the computed values of the shear stress, couple-stress, and rate of heat transfer for different values of the non-similar parameter, , and the transpiration parameter, . It is noticed from the table that there is a good agreement between the LSPM and the MD-BSQLM numerical solutions. The table also gives the computational run time required to obtain accurate results consistent with up to six decimal places. It can be seen from Table 2 that the computational run time for the LSPM is significantly faster than that of the MD-BSQLM. The LSPM computational speed can be attributed to the fact that in the LSPM, discretization is done only in the domains, while discretization is done in both the and domains in the MD-BSQLM. Another observation from the table is that when is small, many terms of the LSPM series approximation are required to give converged results comparable with the MD-BSQLM, and only a few order of approximation is needed when is large.
To validate the accuracy of the LSPM, we define the following solution error:
The errors defined by equation (42) are considered the solution-based error where and denote the solution with the first , terms of the series approximation, respectively. We remark that if the iterative scheme converges, the errors are expected to decrease with an increase in the order of series approximation. To further assess the accuracy of the LSPM numerical method, we have also considered the residual error, which is a representation of the extent to which the solution of the governing partial differential equations (18)–(20) is approximated by the numerical solution. Accordingly, we define the maximum error of the residual as where and are the governing nonlinear PDEs defined as
and are the LSPM approximate solutions.
Tables 3–5 present the variation of the solution error of the approximate numerical solutions of , , and against increasing order of the LSPM series approximation for different values of . We observe the error reduces as the order of series approximation increases. Also, as becomes larger, the error becomes smaller because the series expansion was done inversely proportional to . The method’s success is linked to the fact that series expansion was done inversely about a large parameter, and the Chebyshev spectral collocation method was used to solve the higher-order perturbation equations in the space direction.
Figure 2 shows the velocity profile graph for the vortex viscosity parameter . We notice from the figure that the velocity profile decreases with an increase in the vortex viscosity parameter . This suggests that as the value of increases, the coupling between Equations (9) and (10) increases, thereby causing a reduction in the fluid velocity. This is due to the extra mixing of fluid layers due to the enhancement of shear stress and the shear stress coefficient increase with an increase in the values of . A similar trend is observed in the study carried out by Ishak et al. [10] and El-Dabe et al. [23].
The effect of vortex viscosity parameter on the microrotation profile is shown in Figure 3. We observe from the figure that initially the microrotation profile is an increasing function of but later decreases for larger values of . As the values of increase, the fluid boundary layer increases near the boundary, but the opposite is observed after a certain value of .
The variation of the constant parameter on the velocity and microrotation profiles is sketched in Figures 4 and 5. We observe from Figure 4 that the velocity profile is enhanced for higher values , while, in Figure 5, it can be seen from the graph that there is a decrease in the microrotation profile when increasing the constant parameter. Here, the different values corresponding to indicate different fluid natures. For instance, when , an intense concentration that is, near the wall, represents concentrated particle flows in which micro-elements close to the wall surface cannot rotate. Also, indicates the vanishing of the antisymmetric part of the stress tensor and denotes weak concentrations. Again, is used to model turbulent boundary layer flows. This assumption is invoked to allow the field of equations to predict the correct behavior of the limiting case when the microstructure effects become negligible, and the microrotation reduces the angular velocity. These findings are consistent with those of Ibrahim et al. [45].
The effects of transpiration parameter on the velocity, microrotation, and temperature profiles are displayed in Figures 6–8, respectively. Figure 6 indicates that as the transpiration parameter increases, the velocity profile increases. It can be observed from Figure 7 that the microrotation profile decreases at the initial stage as increases but later increases with an increase in . In Figure 8, we notice that an increase in transpiration parameter results in a decrease in the temperature profile. The thermal boundary layer thickness also decreases as the transpiration parameter increases.
Figure 9 shows the effect of the porous medium parameter on the velocity profile. It can be noted that increasing the porous medium parameter increases the velocity profile. This suggests that the result could be used in determining the relevance of enhanced oil recovery in engineering. This observation is consistent with a recent study by Jimoh et al. [46].
Figure 10 illustrates the effect of the Prandtl number on the temperature profile. The temperature profile, as well as the thermal boundary layer thickness, decreases with increasing . The effect is more pronounced with smaller values of because smaller values of increase the thermal conductivity of the fluid, and hence, heat diffuses away more rapidly from the heated surface than for higher values of . These results are in agreement with those obtained by El-Dabe et al. [23], Ibrahim et al. [45], and Hayat and Qasim [47].
Figure 11 depicts the temperature profile for the heat generation parameter . The figure shows that the temperature profile is enhanced by increasing the heat generation parameter. Physically, the presence of within the boundary layer causes energy levels to increase. As a result, the increasing value of the heat generation parameter increases the thermal boundary layer thickness. A similar trend has been reported in the investigation carried out by Gnaneswara Reddy [48].
The effect of radiation parameter on the temperature profile is shown in Figure 12. We notice from the figure that an increase in the radiation parameter results in an increasing temperature profile within the boundary layer. This is due to a rise in thermal boundary layer thickness. Similar results have been observed in Hayat and Qasim [47], and Gnaneswara Reddy [48].
Figures 13 and 14 represent the wall shear stress for different values of vortex viscosity parameter when the constant parameter and , respectively. We observe from Figure 13 that the wall shear stress is decreased by increasing the vortex viscosity parameter. Also, in Figure 14, it is observed that the wall shear stress is enhanced by increasing the vortex viscosity parameter. A similar trend was observed by Rees et al. [9].
Figure 15 illustrates the influence of the constant parameter on the wall shear stress. From the graph, we observe that an increase in the constant parameter leads to an increase in the wall shear stress values. This flow pattern with a similar parameter effect has been reported by Rees et al. [9].
Figures 16 and 17 are prepared for the effects of vortex viscosity parameter on the gyration component at the wall when the constant parameter and . It is evident from Figure 16 that the gyration component at the wall is decreased for higher values of the vortex viscosity parameter. In addition, it is seen from Figure 17 that the gyration component at the wall is an increasing function of the vortex viscosity parameter.
Figure 18 depicts the influence of heat generation parameter on the local heat transfer. It is noticed that the values of the local heat transfer decrease with an increase in the heat generation parameter . This is because increasing values of increases the thermal boundary layer but decreases the temperature gradient at the boundary.
Figures 19 and 20 shows how the local heat transfer varies with the radiation parameter and the Prandtl number . Figure 19 shows a decrease in the values of the local heat transfer with an increase in radiation parameter. In contrast, an increase in the values of the local heat transfer is observed in Figure 20 with increasing values of Prandtl number. This can be explained by how as increases, the thermal boundary layer thickness decreases, and the wall temperature gradient increases. It is worth mentioning that the graphs displayed in Figures 13–20 were obtained using the MD-BSQLM.
The variation of the residual errors for , , and , respectively, against increasing order of the LSPM series approximation is depicted in Figures 21–23. It can be seen from the figures that there is a decrease in the residual errors with an increase in the order of LSPM series approximation. For different values of considered, we notice that the residual curve is nearly uniform, and the rate of reduction of the residual error is observed to be linear. This indicates that the LSPM will converge up to a certain maximum level, equivalent to the level at which the curve levels off. In Figure 21 which is the graph for equation, it can be seen that the maximum level is at least and about in Figures 22 and 23 which corresponds to the graphs for and equations, respectively. A further observation from the graphs is that as becomes larger, the accuracy of the LSPM improves. The accuracy of the LSPM not deteriorating as becomes larger can be accorded to the fact that series expansion was done about a large parameter. This result indicates that the LSPM gives accurate results that are valid for large parameters or intervals.
Figures 24–26 shows the variation of the solution error of , and against increasing order of the LSPM series approximation. The decrease in the solution error as the order of approximation increases is an indication that the numerical scheme converges. It can be seen that the error appears to be very small after ten orders of approximation. This is because the series expansion was inversely proportional to . This observation accords with Tables 3–5.
The graphs of the LSPM solution for , , and against the order of series approximation is shown in Figures 27–29. We note a decrease in LSPM solution graphs as the order of approximation increases. Full convergence is said to have been attained when the solution graphs start to plateau off. This further confirms that the LSPM numerical scheme is accurate and converges.
5. Conclusion
We have studied the boundary layer flow and heat transfer in a micropolar fluid past a permeable flat plate with heat generation and thermal radiation effects. Approximate numerical solutions of the nonlinear system of partial differential equations were obtained using a large parameter spectral perturbation method (LSPM). The investigation is mainly aimed at demonstrating the applicability of the large parameter spectral perturbation for nonlinear systems of coupled partial differential equations related to the fluid model considered in this study that cannot be solved exactly even with methods that look for series solutions. In addition, the model studied the effects of heat generation and thermal radiation on boundary layer flow and heat transfer in a micropolar fluid past a permeable flat plate. The results from this study demonstrate that the LSPM gives accurate results valid for higher values of the non-similar parameter . Applying the spectral collocation method discretization only in the space domain in the LSPM led to significant computation time savings. The LSPM numerical method validation was established by comparing with the multidomain bivariate spectral quasilinearization method (MD-BSQLM), and a good agreement was achieved between the two sets of results. The LSPM and MD-BSQLM approximate numerical solutions were also found to be in agreement with those existing in the literature. Numerical simulations were conducted to show the accuracy of the LSPM, and it was observed that the accuracy of the method improves as becomes larger with the residual and solution errors approaching zero. The heat generation and thermal radiation parameters were found to increase both the temperature and thermal boundary layer thickness. These observations are consistent with earlier findings in the literature. Also, increasing the heat generation and thermal radiation parameters reduces the values of the local heat transfer. The methods used in this study proved to be reliable and easy to use, thereby making it an excellent numerical tool for solving coupled nonlinear systems of PDEs related to the one investigated in this study. In the future, we intend to establish if the iterative scheme can be suitable for solving coupled nonlinear PDEs in three dimensions whose solutions cannot be sought exactly even with series solution methods.
Data Availability
No data was used in this study; calculations were performed from code simulations in Matlab.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this article.
Acknowledgments
This work is supported by the National Research Foundation (NRF) of South Africa (Grant Number 129490), the DST-NRF Centre of Excellence in Mathematical and Statistical Sciences (CoE-MaSS), South Africa, and the Central University of Technology, South Africa.