Abstract
The near-space hypersonic aerodynamic glider has strong maneuverability in wide flight envelope. The glide is generally achieved in a smooth manner with no or weak altitude oscillations in the altitude. The maximum lift-to-drag ratio glide is a typical trajectory that can approximate the maximum glide range, which is a crucial indicator for a glider. However, another important indicator, the minimum glide range, which is used in some time-sensitive missions and is expected to reduce the velocity to a specified threshold in a short longitudinal range, is difficult to be realized. In practice, the excessive velocity or energy is usually dissipated during lateral manipulation, wherein either the entire glide range or the glide time is not shortened. An innovative guidance strategy is proposed for achieving the minimum glide range based on a typical maximum glide scheme and bang-bang control scheme only on the longitudinal plane, and the flight time can be reduced considerably. Then, a practical extended state observer based pitch control is utilized to efficiently track the bang-bang command within a wide velocity envelope to achieve the guidance objective. Extensive simulation results demonstrate the effectiveness of the proposed methods.
1. Introduction
In the past several decades, hypersonic vehicles have attracted worldwide attention in aerospace industry due to their high flight speeds. Among them, powered cruise vehicles and unpowered aerodynamic gliders have different trajectories. The powered cruise hypersonic vehicle employs a cruising mode with a fixed altitude and velocity as the major phase to achieve long flight range. The cruising altitude and velocity can be optimized simultaneously to maximize the fuel efficiency in order to obtain the maximum flight range. Note that the minimum flight range is also an important indicator which reflects the lower utility range of a flight vehicle, especially for a weaponry system. Based on these two indicators, the entire utility range of a flight vehicle can then be determined. Since the scramjet can be switched off at any time, the minimum cruising range can be very small, and then the utility range of a cruising hypersonic vehicle is only determined by the maximum cruising range. However, it is different for the unpowered hypersonic gliders, which generally have lift-to-drag ratios of not less than 2 in near space. In fact, the maximum glide range can almost be determined based on the angle of attack corresponding to the maximum lift-to-drag ratio at each operating point. Unlike the powered cruising hypersonic vehicles, the minimum glide range of a hypersonic glider is more critical when the velocity must be reduced to a specified threshold as soon as possible such that the subsequent actions can be conducted to deal with time-sensitive missions. Otherwise, some subsystems, for example, the terminal seeker, the thermal protection system, and/or the airframe, will be affected severely. Nowadays, the dominant means of dissipating the excessive energy is based on lateral manipulation while the vertical motion should not be changed since the airframe of the hypersonic glider is only suitable for a narrow altitude envelope in near space. It should be noted that the total glide range or the glide time is still not minimized for this strategy, which just makes a trade-off between the glide time and the longitudinal projection of the total glide range. In other words, the minimum longitudinal glide range can only be achieved with an extended glide time. There are two problems with this scheme: (1) this extended glide time is not favourable for the time-sensitive missions accompanying the minimum glide range requirement, and (2) there might be no lateral maneuvering space to avoid metropolitan areas and adversary environments. Therefore, a minimum time glide guidance strategy on the independent longitudinal plane is worth investigating in order to enhance the flexibility of hypersonic gliders.
In academia, considerable progress has been achieved for reentry glide guidance, especially for the space shuttle and reusable launch vehicle missions which require accurate landing at specific areas. Many autonomous, adaptive, and robust reentry guidance laws have been proposed, and the typical methods include baseline guidance, numerical predictor-corrector method [1], evolved acceleration guidance logic for entry (EAGLE) [2], and quasi-equilibrium glide method [3]. All these methods are comprehensive combinations of longitudinal and lateral guidance laws to ensure terminal accuracy, and the bank angle is used to perform lateral manipulation to consume excessive energy. The angle of attack sequence is usually predefined offline according to certain criteria. The fundamental principle is fully utilizing the large lift-to-drag ratio of a hypersonic glider. For example, the angles of attack of the maximum glide range approximately correspond to the largest lift-to-drag ratios to make the glide as far as possible with quite small elevation angles that guarantee the altitude to vary slowly. This type of trajectory exhibits a certain degree of oscillation, which is good for evading hostile interception but unfavourable for the manned flight vehicle. Thus, a quasi-equilibrium glide method [3] was proposed to tackle this problem by determining the angle of attack based on the principle of dynamic balance of the elevation angle. In such a way, the trajectory is smooth without any oscillation at the cost of losing a certain glide range. The objective of these reentry designs is to land a flight vehicle at a specified location, which has sufficient manipulation capability margin for this kind of missions, and the extremal capability does not need to be investigated thoroughly. However, the limit ranges of a weaponry system are determined by its crucial indicators, including the minimum flight range. To the best of our knowledge, there are no open reports on the minimum glide range design on an independent longitudinal plane.
During the deceleration process of a hypersonic glider, the Mach number will be reduced considerably by almost 10. This will cause remarkable time-varying dynamic characteristics in a wide velocity envelope, and therefore, an effective attitude controller is necessary to meet the opposite extremal glide range requirements such that the original guidance objectives can be achieved. There are many research results on the attitude control for hypersonic vehicles [4], mostly in the longitudinal plane with strong time-varying characteristics, and there are few lateral control methods, which are basically concerning coupling dynamics. In the literature, the most popularly used approaches are dynamic inversion [5], back-stepping control [6], sliding mode control [7–9], neural network control [10–14], and disturbance observer-based control [15–17]. In recent years, prescribed performance control is also applied to hypersonic vehicles [18]. For most of the above advanced control methods, sufficient knowledge about the aerodynamic characteristics is required and complicated modelling must be included in the controller design, which is not preferred in practice. Moreover, most of these methods completely differ from the traditional PID-type controller, limiting their acceptances by practitioners. Nowadays, the model free and weak model control methods for controlling hypersonic vehicles have attracted a lot of attention worldwide. The characteristic model-based all-coefficient adaptive control [19, 20] and reinforcement learning adaptive critic control are two typical representatives [21–26]. The former one supplements necessary coefficient constraint in the identification process when using a first/second-order characteristic model to describe the dynamics in order to ensure a numerically stable result. The adaptive critic control directly establishes the control value function without using a dynamic model, which is appealing to academia. However, these adaptive control strategies are far beyond the traditional PID controller framework and substantial expertise is needed in their implementations. An effective and weak model-dependent control strategy, which has a similar framework of the PID controller, is desired to regulate the attitude in reality. Recently, active disturbance rejection control (ADRC) [27, 28], as an extension of traditional PID controller, has been widely used in industrial applications [29–38], and flight control is possibly the most active arena of ADRC [15, 38–42]. In the design, the ADRC only needs an open-loop control gain, which can be obtained relatively accurately than other characteristic parameters and can significantly reduce difficulty of modelling.
In this paper, a unified guidance and control strategy is proposed for an unmanned hypersonic glider to realize the maximum and minimum glide ranges completely on an independent longitudinal plane. The maximum glide range trajectory is selected as the baseline. The aerodynamic characteristics of the drag and lift are sufficiently utilized to produce bang-bang angles of attack on the basis of the maximum glide trajectory to quickly reduce the velocity while the vertical trajectory shape can be kept almost similar. The guidance sampling interval serves as a compromising parameter among the realizable minimum glide range, the realistically generating trajectory shape, and the tracking capability. Then, an extended state observer based attitude controller is designed to track the fast time-varying bang-bang commands. The graphical tuning rules are also presented. Extensive simulation results demonstrate the effectiveness of the proposed method.
The major contribution of this paper can be summarized as follows. (1)A unified energy dissipation strategy solely on the longitudinal plane is proposed according to the aerodynamic characteristics of the near space glider(2)A stability margin tester-based attitude control is provided to achieve satisfactory control performance within a large velocity and altitude envelope by only using a linear time-invariant (LTI) controller
The remaining parts of paper are organized as follows. The mathematical model and the problem formulation are presented in Section 2. A brief investigation of aerodynamic characteristics is conducted in Section 3. The preliminary guidance design is provided in Section 4 as the foundation for further extension. The main deceleration guidance strategy is proposed in Section 5. The extended state observer-based attitude control and its tuning method are offered in Section 6. The numerical design and extensive simulation results are presented in Section 7. Finally, the conclusions are given in Section 8.
2. Mathematical Model with Problem Formulation
2.1. Mathematical Model
In this study, a hypersonic glider model is used and described as where is the velocity of the hypersonic glider; is the drag; is the lift, is the elevation angle, is the angle of attack, is the pitch angle, is the pitch angular rate, is the glide range, is the altitude, is the mass of the glider, is the pitching moment, is the pitching moment of inertia, is the radius of the Earth, is the gravitational constant, is the drive voltage of the elevator , and is a control law to be designed.
Here, a combined aerodynamic model of [3, 40, 42, 43] is employed, which is slightly modified so it is suitable for a wider range of angle of attack. The aerodynamic forces can be calculated as where m2 is the aerodynamic reference area, is the dynamic pressure represented by and and are the drag and lift coefficients, respectively, and they can be fitted by
Here, is the angle of attack in the unit of degree, and the Mach number is defined as where is the normal sound speed. The pitching moment can be calculated as where m is the aerodynamic reference length and where
Here, is the elevator in the unit of degree.
The mass of the glider is . The pitching moment of inertia is .
The elevator dynamics can be described as where .
2.2. Problem Formulation
In this paper, a practical unified guidance strategy will be sought to realize the maximum and minimum glide ranges for the proposed unmanned hypersonic glider model only on the longitudinal plane without lateral manipulation. In addition, a reliable control scheme is also needed to efficiently achieve this guidance law within a wide Mach number interval.
The process constraints are where is the normal acceleration whose range is limited by the structural strength of the airframe and all embedded components. Moreover, we expect to be less than and to be less than 3 after a short period of trajectory jump in the initial phase.
We specify that the initial Mach number is and the terminal Mach number is . The initial altitude is , and the terminal altitude is near with an allowable miss distribution of . Because a terminal seeker is on-board to perform the subsequent terminal guidance, a reasonable operation condition range with certain width instead of a group of terminal equality constraints is sufficient for a weaponry system. In a word, an inequality terminal constraint of altitude rather than an equality constraint is practically useful. This is an optimization problem for the midcourse guidance, unlike the case of space shuttle.
3. Aerodynamic Analysis
For the hypersonic glider guidance design, the first step is investigating the aerodynamic characteristics. According to (1) and (2), the drag and lift coefficients can be illustrated in Figures 1 and 2, respectively. It is clear that both the drag and the lift coefficients monotonically increase with respect to for each Mach number. In addition, is almost proportional to and is nearly linear in terms of based on (1) and (2) or Figures 1 and 2, respectively. This observation matches our empirical experience. We will make full use of these properties later in the unified guidance design.
Because the maximum glide range can be approximately achieved at corresponding to the maximum lift-to-drag ratio [44], we explicitly plot these ratios in Figure 3. It can be observed that the largest lift-to-drag ratio can be obtained at around for each Mach number, and these largest ratios are all about 2, which is the foundation for the subsequent guidance design.
4. Preliminary Maximum Glide Range Guidance Design
In this section, the baseline guidance design for the maximum glide range will be determined. Assuming an ideal attitude controller and ignoring the control error, the mass point model of the hypersonic glider is considered in the guidance design.
According to the above aerodynamic characteristic observation, we can present a guidance strategy by fixing . The corresponding important states are shown in Figures 4–6, respectively. The glide range is 2710 km with a glide time of 788 s. The altitude exhibits obvious oscillations and exceeds the upper limit of 40 km sometimes while the elevation angle varies between . This obvious oscillation results from the fact that no damping effect exists in the pure guidance. This strategy is not preferred in practice due to its persistent trajectory oscillation, which is harmful to the thermal protection system.
To alleviate the altitude persistent oscillation in the pure guidance, we can implement pitch angle-based guidance to enhance the altitude smoothness because can be regarded as a damping variable for the altitude. Considering a small during the glide, we can straightforwardly let to approximately realize the maximum glide range, wherein is the value of corresponding to the maximum lift-to-drag ratio. In this paradigm, . The corresponding states are shown in Figures 7–9, respectively. It can be seen that the altitude trajectory is smoother with smaller oscillations comparing with the scheme of fixed angle of attack. In this case, the glide range is 2697 km with a glide time of 781 s. This range loss is relatively tiny comparing with the constant strategy. In addition, the requirements for and are also satisfactorily met. Therefore, this strategy can be used to replace the previous one.
To completely eliminate the altitude oscillation, we further investigate the quasi-equilibrium glide method [3]. In this approach, the dynamic balance condition is used to derive the altitude-velocity profile, which can generate a proper attitude command according to the current velocity. According to (1), is equivalent to
By utilizing numerical methods, for example, the fsolve function in Matlab, we can obtain the altitude profile with respect to the velocity as shown in Figure 10, which can be fitted by a fourth-order polynomial. We design a guidance law as where is the initial pitch angle, is the altitude reference in Figure 11, and is a constant gain. The related states are illustrated in Figures 12–14, respectively. It can be seen that the altitude trajectory is monotonously decreasing without any oscillation and the altitude can track its command perfectly even with a simple guidance law of (15). In this case, the glide range is 2565 km with a glide time of 743 s, and the range loss is relatively large and the cost for the altitude smoothness is too high because deviates from to some degree only to ensure the smoothness of altitude trajectory according to Figure 13.
In summary, the second scheme, , is a preferred scheme to balance the oscillation and the maximum glide range. Therefore, we concentrate on this approach to design a minimum glide range guidance in a unified way.
5. Minimum Glide Range Guidance Design
Here, we attempt to provide the guidance law for achieving the minimum glide range based on the maximum glide range trajectory to demonstrate the effectiveness of fast deceleration without any lateral maneuver. In the quasi-equilibrium glide method, is always less than zero or is monotonously decreasing; therefore, the high altitude advantage during the glide will soon disappear. Thereafter, we hope to keep the altitude trajectory as much as possible to take the full advantage of high lift-to-drag ratio characteristics in near space while enhancing drag considerably. The minimum glide range guidance law can be designed as
The central idea is to maintain the glide mode by using a type of bang-bang control strategy: the maximum lift is produced when the altitude is decreasing and vice versa, so that the altitude is almost constant. This is quite similar to the variable structure control. Considering the linearity of with respect to , we have to ensure that the average can approximately be as an effective glide mode. On the other hand, the effective range is . As is almost proportional to , we can select and to achieve the maximum drag while making the altitude trajectory shape almost unchanged. It is obvious that this bang-bang control strategy will lead to a somewhat high-frequency oscillation in the pitch command; therefore, the guidance sampling interval plays a crucial role in this design. Let the guidance sampling intervals be 50 ms, 500 ms, and 5000 ms; the glide ranges are 1495 km, 1494 km, and 1484 km, with glide times of 443 s, 443 s, and 440 s, respectively. Comparing with the original glide range of 2697 km, it is obvious that this attitude self-oscillation strategy can reduce the glide range considerably. Next, the important variables are illustrated in Figures 15–17, respectively. Due to the page limitation, we only show the case of 5 s guidance sampling interval since the shorter one will make the attitude controller rather difficult to track the command considering the reasonable bandwidth of the attitude controller and the elevator constraints such that the guidance command is exactly ineffective. Comparing Figure 15 with Figure 7, it can be observed that the altitude has tiny oscillations around the baseline trajectory from the constant pitch scheme. In Figure 16, the elevation angle oscillates with small amplitudes around zero, satisfying the specified constraints.
6. Pitch Attitude Controller Design and Tuning
6.1. Attitude Controller Design Based on an Extended State Observer
In the abovementioned generic guidance design for both the maximum and minimum glide ranges, the Mach number will be reduced by 10. This substantial variation imposes a challenge for the attitude controller, which should have robust dynamic performance against time-varying environment and certain aerodynamic uncertainties. Moreover, the pitch attitude controller must track the bang-bang command in the minimum glide range mode in a fast way; otherwise, the longitudinal range cannot be minimized efficiently. As the proposed guidance design is practically implementable, we will subsequently derive a useful attitude controller.
Because the elevator dynamics (12) is usually much faster than that of the airframe, it can be ignored in the controller design and it is only included in the tuning process as well as simulations. Therefore, we do not distinguish between and in this subsection although they have different units of degree and voltage, respectively. Consider the dynamic equation where is
Define where is defined as an extended state [27], which includes the total right-hand-side dynamics of (18) except the explicit effect of . Here, the extended state includes the major uncertain aerodynamics, and the time-consuming modelling procedures can be avoided. This is the foundation for the following controller design. Equation (18) can then be reformulated as where is the total disturbance, which makes the real plant (18) deviate from an integrator. Equation (20) treats as a whole regardless of its physical meaning and origin. We can design a state observer for (21) as where is the estimation of , wherein approximates . Thus, this observer is called extended state observer (ESO). To simplify the tuning process, Gao [28] proposed a linear ESO (LESO) that can be parameterized as where the observer bandwidth, , is the only tunable parameter. When the control signal employs the form of where is a virtual control variable, the original plant (18) can be approximated as and then
Because this is a second-order integrator, we can design a PD controller to achieve no steady-state error, wherein and are the proportional and derivative gains, respectively.
In summary, the pitch attitude controller is
It can be seen that the control law (28) is an approximate PID-type controller in form. The effect of is an extension of the conventional integral feedback. This kind of controller is likely to be understood by practitioners who are quite familiar with the traditional PID control.
6.2. Attitude Controller Tuning Based on a Stability Margin Tester
The proposed LESO-based attitude controller is a PID-type controller. It is well known that there is no rigorous guarantee for the nominal closed-loop stability when tuning a PID controller. This is because (1) the PID controller is not established based on the Lyapunov theory which is the theoretical basis for closed-loop stability and (2) the order of a PID controller is generally lower than that of the plant, which causes difficulty in analyzing nominal closed-loop stability. However, there are mature procedures for tuning PID controller from the perspective of robustness in the presence of uncertainties. An explicit tuning rule is necessary for a control algorithm if it can be used in practice. In this section, we present an explicit tuning methodology for the proposed attitude controller, and the gain and phase margins (GPMs) are used to evaluate the robustness of the designed control system. In [38, 40, 42], we developed a graphical method for fast parameter tuning of linear ADRC based on the stability margin tester. This explicit approach utilizes the D-decomposition philosophy to partition the parameter space. Due to the page limitation, this approach is not repeated here and the readers can refer to [38, 40, 42]. The only requirement for this tuning is the linearized model below at the trimming point. where the aerodynamic coefficients are defined as
Here, denotes the partial derivative of with respect to , and represents the partial derivative of with respect to . Because or each part of is a polynomial of , , and in this paper, all these partial derivatives can be explicitly calculated. Note that the linearized small-perturbation model around a trimming point is traditionally only valid around small the angle of attack, which does not hold true for the large , for example . However, the tuning result is helpful for the scenario of large angles of attack.
7. Numerical Example and Simulations
In this section, the effectiveness of the proposed method will be evaluated by numerical simulations.
7.1. Operating Point Analysis and Tuning
At first, the aerodynamic coefficients can be obtained at several operating points. The principles of operating point selection are covering the major dynamic variations of the nominal trajectory and selecting as few operating points as possible. In the glide phase, the altitude variation is not quite remarkable and the velocity variation is more distinct. In addition, the altitude matches the velocity on the nominal trajectory to some degree. Therefore, we select the operating point every 1 Mach number. For the constant pitch case of the mass point trajectory shown in Figure 7, we can obtain the aerodynamic coefficients as shown in Table 1.
Next, the stability margin tester is used to tune the controller. We specify based on our empirical experience. In fact, can be the first choice for various flight vehicles. According to Table 1, has a tight variation range, which can facilitate the controller design because this is the only parameter related to the characteristics of the glider. We select and take all operating points into account. The stability margin testers for all 11 operating points are shown in Figure 18. Here, each curve corresponds to a specific GPM combination. On one of the two sides bounded by the curve, the GPMs more than the given specification can be achieved, and vice versa. We use and to denote the gain and phase margin specifications, respectively. In Figure 18, four specifications are as follows: (1) or the nominal stability with red colour, (2) or the double gain perturbation stability with green colour, (3) or the half gain perturbation stability with cyan colour, and (4) or the phase perturbation stability with magenta colour. Unlike the abovementioned curves related to robust stability, we add a nominal dynamic performance curve corresponding to the loop transfer function with unity feedback and phase margin with blue colour. We recommend the points on this curve in the common robust stability area to have high priority as tuning results, which can reduce the 2-dimensional search to a linear one. For all these operating points, we can employ an LTI controller when there is a common intersection area. Fortunately, there is a common area in this scenario and we can specify as illustrated in Figure 18. The corresponding step responses for all these operation points are shown in Figure 19. It can be seen that satisfactory dynamic performances can be achieved.
7.2. Nominal Nonlinear Simulation
Then, the tuned parameters are used in the nonlinear simulation. Three cases are considered when the pitch error limit in (28) is included to avoid elevator rate saturation. Three error limits are specified as , , and , respectively. In addition, we impose a hard constraint on the elevator deflection rate at when taking the practical implementation into account, which is a fairly conservative specification in practice. The glide ranges are 1867 km, 1680 km, and 1649 km, respectively. The deceleration effect is not remarkable as that in the mass point trajectory scenario because of the control error. It is evident that the deceleration effect is more remarkable when the control action is aggressive and there is a risk of elevator rate saturation. This is because the pitch command has the unique property of large step and short steady time simultaneously, which is a huge challenge for controller realization. Without loss of generality, we only investigate the case of in detail. The corresponding states are shown in Figures 20–26, respectively. Comparing Figure 20 with Figure 15, it can be found that the altitude shape is kept well. According to Figure 22, the speediness of the command can be almost guaranteed with a certain accuracy. The elevator deflection rate and the pitch angular rate are illustrated in Figures 25 and 26, respectively. The maximum pitch angular rate is less than , and the peak angular rate is produced in an impulsive manner, and they are within the measurement range of conventional gyro with considerable margin. It is clear that the control performance is still satisfactory even considering all the practical constraints. In fact, the difference between and , , can be selected to make trade-off between the vibration amplitude of and the achievable minimum glide range. Smaller corresponds to smoother attitude and weak deceleration effect, and vice versa. This variable can be regulated to fit a concrete scenario and enhance the flexibility of this strategy. For example, the glide range will increase from 1680 km to 2271 km when is changed from to or and , and the deceleration effect is weaker.
Furthermore, the maximum glide range case is used to assess the tracking accuracy of the controller for the slow time-varying command. Due to the page limitation, only the pitch response is shown in Figure 27. It can be observed that the attitude accuracy can be ensured.
7.3. Monte Carlo Simulations
Finally, the robustness of the attitude controller should be evaluated because we employ an LTI controller to deal with nonlinear dynamics within a wide velocity envelope. As there are two kinds of aerodynamic uncertainties such as aerodynamic forces and moments, they should be considered separately. In practice, the relative errors in the aerodynamic forces are much lower than those in the aerodynamic moments. According to existing experience, there are nearly 10% relative errors in the aerodynamic forces and 30% relative errors in the aerodynamic moments. As a matter of fact, the aerodynamic force uncertainties mainly affect the mass point trajectory, and the aerodynamic moment uncertainties fundamentally influence the control performance. To evaluate the robustness of the controller, we concentrate on the pitching moment uncertainties. From the aerodynamic point of view, the pitching moment generally includes three parts and they are the stability moment in terms of , the damping moment in terms of , and the manipulation moment in terms of , which are corresponding to , , and , respectively. For the sake of brevity, we take the influence of and as a whole with multiplicative uncertain factor uniformly distributed within the interval of . The uncertain factor related to , , is also uniformly distributed in . We conducted 500 Monte Carlo simulations, and the related variables are shown in Figures 28–33, respectively. It can be observed that all variables are well within their specified ranges. The long period variables, the altitude and Mach number, are almost insensitive to the aerodynamic moment uncertainties. According to Figure 30, the small oscillatory amplitudes of the elevation angles can be maintained during the whole period. In summary, the proposed LTI controller can track the fast and large bang-bang commands within a wide range of velocity quite well subject to diverse constraints.
8. Conclusion
An innovative guidance strategy was proposed for the minimum glide range realization based on a typical maximum glide range scheme by properly using the bang-bang control scheme only on the longitudinal plane. Then, a practical extended state observer-based pitch control was utilized to track the fast bang-bang command rather efficiently within a wide velocity envelope to achieve the guidance objective. Both the control speediness for the minimum glide range and the control accuracy for the maximum glide range can be realized simultaneously by using an LTI controller, and the flight time could be reduced considerably. Extensive simulation results demonstrated the effectiveness of the proposed methods. This work can provide helpful guidelines for the near space hypersonic glider design.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors disclosed receipt of the following financial support for the research, authorship, and/or publication: National Natural Science Foundation of China (Nos. 62073177, 61973175, and 51777013), South African National Research Foundation Grants (Nos. 112108 and 112142), South African National Research Foundation Incentive Grant (No. 114911), and Eskom Tertiary Education Support Programme Grant of South Africa.