A new fixed‐point algorithm to solve the blade element momentum equations with high robustness

The most common approach to solving blade element momentum (BEM) equations is through fixed point method. The fixed point method can provide reliable solutions with high precision, yet the robustness of the method has been challenged when infrequent failures of converging to a physical solution are found for some design space. Though the lack of robustness is alleviated by applying two improved algorithms, their shortcomings should not be of an understatement. Ning's method can result in a converged yet nonphysical solution, and Sun's method decreases the computational efficiency remarkably. To overcome these setbacks, a new algorithm has been proposed in this paper. A clear classification of a wind turbine operating states has been given first to correct the thrust relation for a > 1, followed by discussions of three failure cases encountered during solving BEM equations. Then, the new algorithm with three major modifications has been introduced and explained. The test of Section 4 reveals that the decreasing rf technique has a positive effect on improving the robustness. Besides, the first two tests in Section 5 prove that the new thrust equation can greatly enhance the robustness, and Aitken's squared process can significantly strengthen the efficiency. The results show that all three modifications contribute to offering a new FPA with high robustness and satisfactory computing efficiency, which serves as the best option for solving the BEM equations.


| INTRODUCTION
BEM (blade element momentum) theory is currently widely used in wind turbine-related research and industry due to its simplicity, efficiency, and reliability, which is named since it is a combination of the momentum theory and blade element theory (or strip theory). The most common approach to finding the solution is through an iterative process based on fixed point method (FPM). 1,2 Typically, the initial values of the axial induction ratio and tangential induction ratio are set to zero. Iterations then continue until the final convergence has been achieved or a message reporting failure is returned. The convergence criterion is whether the absolute differences errors of both axial induction ratio and tangential induction ratio between the adjacent ones are guaranteed within the error tolerance simultaneously. The specific numerical | 1735 JIN aNd YaNG algorithm developed based on FPM is named fixed-point algorithm (FPA) here for simplicity and consistency.
Though BEM was considered fast and accurate, Maniaci reported that the convergence problem did exist. 3 He claimed several factors were responsible for that, one of which was the singularity at the axial induction ratio of 1. The iteration routine might be trapped there due to a potentially severe oscillation of the thrust coefficient around the axial induction ratio of 1. Thus, a new smoothing technique had to be proposed to alleviate the convergence difficulty despite bringing more complexity alongside. However, the thrust coefficient deriving from the blade element theory could not be solely written as a function of the axial induction ratio. That is to say, all the blade element solution curves were not correctly given because they depended only on the axial induction ratio. Besides, the arbitrary statement that the axial induction ratio greater than 1 meant that the flow entered the vortex ring state was wrong.
Ning claimed that the convergence problems had not been fully resolved especially during optimization when extreme values were considered. 4 A new robust solution was put forward to describe the equations in a different way, whose core idea was to reduce the initial two-equation, two-variable problem to one-equation, and one-variable problem. Then, a robust root-finding method based on bisection could be used to find the solution after determining the searching span. The idea itself was highly appreciated, yet the relations between the rotor operating states and the local inflow angle were not fully understood.
The statement that two half planes corresponding to propeller brake state (a > 1) and momentum/empirical state were divided by the = 0 was not correct (see Figure 1). The propeller brake state can also be spotted with a negative if we traced back to the definition of (tan equals to the ratio of U ∞ (1 − a) to Ωr (1 + a�) ). Assuming that the incoming maintains the positive direction (U ∞ > 0), it can give > 0 for propeller brake state if the wind turbine rotational speed Ωr < 0. The case can happen during the idling analysis when the turbine can rotate in both clockwise and anticlockwise directions.
Secondly, the relation between the thrust coefficient and the axial induction factor was not correctly given by Ning for the propeller brake state. When a is larger than 1, the thrust coefficient C T would rather be negative than positive (see Figure 4). More detailed explanations about this would be given in Section 2. Thirdly, the robustness of Ning's method could still be challenged when it encountered a negative incoming wind speed or rotational speed. As indicated by Figure 2, a converged yet nonphysical solution was returned with Ning's one-equation method during the calculation of the static performance at some given operating condition (U ∞ = 9 m/s, Ω = −10.5 rpm, and pitch = 0 degree) for NREL 5-MW reference wind turbine. The detailed geometry and airfoil used to construct the blade were given by Jonkman et al. 5 As a contrast, the new method proposed in this paper would give a converged and physical solution at the same operating condition. It should be noted that the local speed ratio, solidity, twist, and unique airfoil characteristics would have a great influence on the actual operating state of the blade element. The pitch behavior, gusts, and changes in the wind direction would also have some effects on the operating state.
Studies on convergence issues had also been performed by McWilliam and Crawford by comparing the performance of FPM and the Newton-Raphson method in solving BEM equation. 6 Better results had been achieved with the latter one, yet the instabilities could not be fully eliminated. To argue that FPA can be competitive with appropriate initial values given, Sun et al 7 employed the one-variable concept by replacing the axial induction ratio and tangential induction ratio with the equation of relative flow angle to rewrite the equations in some way that could be used for one-dimensional FPA. Though there are obvious mistakes concerning the mathematical analyses of the singularity issue in that paper, the FPA exhibited rather good performance in finding the solution through the tests. Besides, it should not be ignored that the parameterization raised in Sun's paper was not that reasonable since singularities had been incorporated and the calculation of the zero-lift angle for each radial element brought more complexity.
A robust and efficient method is vital to obtain a reliable solution to the BEM equations, especially during optimal design. The integration of the one-equation concept to the FPM indeed offers a potential solution to tackle the convergence issues, but an appropriate parameterization of the relative flow angle should be adopted to get rid of the singularities if the FPA were to be employed. However, in this paper, the author will preserve the two-variable BEM equations and FPM, and present three major modifications to significantly improve the performance. The key to handle the multiple solutions is to ensure the "best" solution should keep consistent with the angles of attack along the blade.
The main sections are organized as follows. First, an introduction to the BEM theory will be given to facilitate the understanding along with a more comprehensive classification of the operating state of a wind turbine. Next, the causes of convergence issues will be analyzed with the aid of a simple test of solving BEM equations over given space with the traditional FPA. The specific modifications to the BEM equations and the underlying reason will be presented then. Last, comparisons of various methods will be held, and conclusions will be made.

| BEM THEORY
BEM theory is composed of momentum theory and blade element theory. A comprehensive classification of the operating states of a wind turbine is provided first, followed by the presentation of the key equations mentioned in the BEM theory to allow coherent reading and easier understanding. More in-depth discussions on the BEM theory can be found in the literature. 1,2,8,9

| Classification of the operating states of wind turbine
Since the fundamental theory used to study the wind turbine mainly comes from the helicopter theory, it is necessary to trace back to the primary derivation to present a reliable discussion on the operating state. Four regimes, that is, propeller state, windmill state, turbulent wake state, and propeller brake state, are usually employed to describe the wind turbine, whose counterparts can be easily found within the helicopter theory. A summary of the operating states for wind turbine and helicopter machines is given in Table 1.
The wind turbine normally operates in the windmill state where energy is extracted from the incoming flow and used to drive the rotor to generate electricity. The momentum theory breaks down when it enters the turbulent wake state and the high thrust correction needs to be applied. The axial induction ratio that divides the windmill state and turbulent wake state is called the critical axial induction ratio, and the corresponding thrust coefficient is called the critical thrust coefficient. In this article, this critical axial induction ratio is called the first critical axial induction ratio to distinguish it from the other one.
When the axial induction ratio a continues to grow and reaches 1, the energy extracted by the turbine blades from the moving air is totally dissipated into the wake, meaning the thrust alongside is equal to the rate of energy dissipation to the vortices. Yet, an equilibrium state can still be satisfied where no external power needs to be injected into the wind turbine, which can be called the ideal autorotation point.
As a continues to increase to move away from the ideal autorotation point, some external energy is demanded to drive the wind turbine to rotate. At the same time, the wake direction is altered gradually, and a new distinct slipstream boundary with a reversed direction begins to develop. Fundamentally, the momentum theory breaks down at the beginning of the vortex ring state. However, as suggested in Leishman's book, 10 there is a critical axial induction ratio beyond which the power curve conforms to that estimated by the momentum theory well. This argument can further be solidified by Wolkovitch, 11 who gave a guess for the critical axial induction ratio based on the theoretical derivation. In this paper, this critical axial induction ratio is called the second critical axial induction ratio. The detailed calculations of the second critical axial induction ratio are given in Appendix A.
It has to be admitted that there is a remarkable difference between two-second critical induction ratios computed  via two methods. However, this paper will not give a further derivation on the empirical correction for the propeller brake state and the discussion on the accuracy of the second induction ratio is not of our concern for the moment. A simple extension of the high thrust empirical correction of the turbulent wake state beyond the propeller brake state can greatly simplify the discussion on convergence issues. The reason why such assumption is applied can be ascribed to two vital factors: (1) the wind turbine normally operates within 0 < a < 1, which means that it is unlikely that the wind turbine can enter the propeller brake state; (2) the thrust relation derived based on BEM theory shall not be valid for the propeller brake state and a sudden drop of thrust coefficient from the turbulent wake state to the propeller brake state is detrimental for solving the BEM equations as it introduces the singular point. Always keep in mind that the first critical axial induction ratio is less than 1, and the second is greater than 1.
The curt presentations on the operating states can finally draw forth an important conclusion that the thrust equation derived from the momentum theory is not valid for the propeller brake state and an extension of the high thrust empirical correction beyond the propeller brake state in terms of improving the robustness of the solution algorithm. To get an intuitive understanding, the relationships between the induced velocity ratio and climb velocity ratio and between the power ratio and climb velocity ratio are provided in Figure 3. It can be easily seen that when the sum of induced velocity and climb velocity is less than 0, and the power ratio is negative and vice versa. In fact, the axial induction ratio can be obtained by the ratio between the induced velocity v i and −v c (negative sign here since the descent velocity is assumed negative for helicopter machine while the wind speed is assumed positive for wind turbine machine). Therefore, it should be convinced that beyond the ideal autorotation point (v i + v c = 0), external power is demanded to drive the machine to work normally for the vortex ring state (propeller brake state for a wind turbine).
Up to now, the paper has provided a thorough review of the classification of the operating states. Next, the authors feel it the right time to give the overall thrust curve describing the relationship between thrust coefficient and the axial induction ratio (see Figure 4).

| Comprehensive derivations of BEM equations
As assumed by the BEM theory, a wind turbine rotor can be viewed as a multiple annulus element with no radial flow. The incremental thrust dT and torque dQ acting on the annulus at a distance r from the rotor axis can be obtained based on momentum theory and blade element theory, respectively: where is air density, v 0 the incoming wind speed, the angular speed, F the Prandtl tip loss factor, B the number of blade, c the chord length, C n the normal force coefficient, C t the tangential force coefficient, and the relative inflow angle.
The normal force coefficient and tangential force coefficient are defined as: After some simplifications, it yields the axial induction and tangential induction equations: where , the local solidity, is defined as = Bc 2 r . The following thrust coefficient deriving from the blade element theory is usually employed in the traditional FPA: Besides, from the flow geometry, it yields: where r represents the local speed ratio, is defined as r = wr v 0 . Then, local angle of attack can be given as: where is the sum of local pitch angle and local twist angle.
Concerning the tip loss correction, Shen's correction 12 seems to be the most advanced one. Additionally, researchers have proposed various empirical thrust corrections to correct the heavily

| CONVERGENCE ISSUES
Full knowledge of the causes of failure to converge when engaging the initial algorithm will help recognize the solution.
As a result, a simple test of solving BEM over given space will be conducted to aid the analyses. Remark that a = 1 and a′ = −1 will always be a trivial solution and makes no sense except for the representation of the ideal autorotation point, loaded wind turbine, such as Glauert,, 13 Spera 14 Anderson, 15 and Buhl. 16 As a well-known industrial code, GH Bladed 17 also puts forward its thrust correction formula. Only two corrections with the first critical induction ratio are tested in this paper to represent the linear and quadratic models, respectively.
Here, the value of a c in Equation (10) comes from the statement given by Hansen. 2 A brief description of the numerical algorithm of the traditional FPA with Spera's thrust correction model is presented below. If Bladed's thrust correction model was to be adopted, two lines would be modified: the 9th line would be modified from C T ≥ 0.64F to C T ≥ 0.9219F and 10th line from compute a with Equation (10) to with Equation (11).
Generally, a positive relaxation factor (rf) which is not larger than 1 will be adopted when updating the a and a′: where the subscript i is the number representing the i-th iteration. (10) The classification of helicopter operating state Compute the relative inflow angle using Eq. 8 and effective angle of attack using Eq. 9 5 Look up and 6 Calculate and with Eqs. 3 and 4 7 Compute the thrust coefficient using Eq. 7 8 Compute tip loss factor F 9 if then 10 Compute a with Eq. 10 11 else 12 Compute a with Eq. 5 13 end 14 Compute a' with Eq. 6 15 end 16 end | 1739 JIN aNd YaNG which will be treated as a failure though it is a converged solution. Therefore, three representative failure cases can be put forward: 1. a is returned as a trivial solution ( = 0, a = 1 and a′ = −1) 2. the BEM based code falls into an infinite loop 3. a is returned as a complex number

| Failure case 1
Failure case 1 is inherent in the BEM algorithm derived above. The fact that the trivial solution always exists can be easily seen that a and a′ will be equal to 1 and −1, respectively, if and only if = 0. Besides, When tends to 0, (1a) term will always converge to 0 faster than (1 + a′), which can be proved as follows: From Equations (5) and (6), it can be obtained that: Thus, it can be proved that: The explanations why the solution converges to the trivial solution can be given in two aspects. The one is that if the relaxation factor is set to 1, the algorithm may not be able to find the solution other than the trivial solution, which is attributed to the way how FPA works. A relaxation factor smaller than 1 may help alleviate this scenario.
The other one is that the derived Equation (5) is associated with some drawbacks that should not be neglected, which cannot be removed with changing relaxation factor. More specifically, two scenarios discovered can be used to illustrate the drawbacks. One scenario will be that a continuous shift of axial induction ratio from negative to positive and back to negative, which pushes the algorithm to fall into failure case 2. The other one is that the axial induction ratio becomes greater than 1 at some iteration and the algorithm will then gradually converge to the trivial solution. It should be mentioned that the existence of the axial induction ratio greater than 1 during the iterative process is a prerequisite for FPA to converge to the trivial solution, which can be easily seen in how FPA works.
Through analyzing the test results obtained, it can be found Equation (5)  is not that reasonable.
When C n is negative, there is a great chance for a to be greater than 1. As mentioned earlier, a > 1 is not desired during the iteration process as it is a prerequisite for the algorithm to converge to the trivial solution. Also, the trivial solution is not a singular point in Equation (5), which means that the iterations will be ended and considered as a failure once the trivial solution is obtained.

| Failure case 2
Apart from the potential failure case 2 mentioned in the above segment, another factor leading to failure case 2 is a combination of the inappropriate setting of the relaxation factor and the randomness of airfoil characteristics. By changing either the relaxation factor or the airfoil, the FPA could reach a converged solution under the same working condition. But, remember better convergence performance cannot be guaranteed with a smaller relaxation factor. This is to say, a smaller relaxation factor may even bring negative effects on convergence performance. In addition, the iteration may be trapped to random regions, in contrast to the statement given by Maniaci 3 that the infinite loop is limited to the regions close to axial induction ratio of 1. For example, when twist = 8.333333, tip speed ratio = 10.72222, solidity = 0.07889, and rf = 1, the FPA will be trapped into an infinite cycle for FFA-W3-241 where the axial induction ratio deviates far from 1.

| Failure case 3
The failure case 3 only occurs for the particular thrust correction model, which results from a negative square root calculation when updating axial induction ratio. A general linear and quadratic correction can be given as: Equating Equation (17) to Equation (7), it gives: As the term C n F(sin ) 2 will always be a positive term since Equation (19) is adopted for C T larger than a positive number, the discriminant will, therefore, be positive if the coefficients of the linear correction are positive. In other words, failure case 3 can be avoided with a linear correction with all the positive coefficients.
However, it is not the case for quadratic correction. Following the similar procedures, the discriminant for the quadratic correction can be given as: Take the GH correction model for example, by substituting the coefficients into Equation (21), it gives: If C n F(sin ) 2 is smaller than 0.1905, it will result in a negative discriminant and failure case 3. This can occur when the axial induction ratio used to calculate C T in Equation (7) is smaller than −1.1914 or greater than 3.1914. Therefore, there is a possibility in encountering the failure case 3 with the GH correction. Similar analyses can be performed for other quadratic corrections; it turns out that the failure case 3 is associated with the quadratic correction.

| Verification of various cases
Next, the three failure cases will be verified with a test of solving BEM equations over 10 × 10 × 10 grids covering the twist from −5 to 25 degree, the local speed ratio from 0.5 to 12, and the local solidity from 0.005 to 0.1 for a single airfoil, FFA-W3-241. More airfoils will be analyzed in the later sections. The intervals of the above parameters have also been adopted by Ning 4 to test the performance of his bisection method. A linear distribution of the parameters is used here. The FFA-W3-241 airfoil is used to design the 10 MW reference wind turbine by Bak et al. 18 The airfoil characteristics are first obtained in 2D (EllipSys2D) and corrected using Bak's correction model. 19 Although the 3D airfoil characteristics are obtained with the CFD method and extended to full 360 degrees with flat plate theory, it is not the concern to worry about whether the data are completely accurate as the aim was to identify the reasons leading to failure to converge at this stage.
Two correction models shown in Equations (10) and (11), and three relaxation factors, 1, 0.5 and 0.25, have been employed and tested separately. The nondimensional radius of the blade element r/R is chosen as 0.5 to assure that the tip loss factor can be set to 1 for the moment. The convergence criterion is that the variations of the axial and tangential induction ratios should both be within the error tolerance = 1 × 10 − 8 simultaneously. If the deviations of the a from 1 and the a′ from −1 are within 1 × 10 − 6 , the solution will be considered as the trivial solution even though it converges.
As is shown in Table 2, as rf decreases from 1 to 0.5, the total number of failures drops, mainly due to the elimination of the failure case 2. This is not the case when rf is further reduced. On the contrary, when rf is reduced from 0.5 to 0.25 using both correction models, the number of failure case 2 increases rather than decreases, which confirms the viewpoint that a smaller relaxation factor may bring negative effects on the convergence performance.
Neither of the two thrust correction models is superior to the other one. Fewer number of failure case 1 occurs with Bladed's model but more of the failure case 3. As a consequence, the choice of the thrust correction model actually will have a minor influence on the overall failure rate. To put forward a new induction factor equation, a revisit to Equation (5) is essential. If a more reasonable equation relating C T , a i , and a i+1 can be built, the performance of the algorithm ought to be improved. By equating the thrust coefficients deriving from the blade element theory and momentum theory, it gives: By eliminating the (1a) term on both sides, Equation (23) can be organized as: Instead of moving all terms including a to the left-hand side and obtain Equation (5), Equation (24) is kept unchanged and embedded in the improved fixed-point algorithm. It is interesting to find that the denominator can only be composed of the tip loss factor and the axial induction ratio. In particular, during the choice of the new axial induction ratio equation, the singularities are confined to known locations containing = 0, = and = − .
Rewrite Equation (24) as: Note, Equation (25) is only valid for the lightly loaded turbine. With Equation (25), the trivial solution becomes a singular point, which will not be returned as a solution anymore. Besides, special measures have been taken to deal with the singular point in the improved algorithm, including the automatic relaxation factor decreasing technique. The automatic relaxation factor decreasing technique means that whenever the iteration step reaches some integers times the given iteration limit and the algorithm still does not converge, the relaxation factor will automatically decrease to half of the relaxation factor at that time. Generally, the initial rf is set 1 and converged solutions can be obtained before rf = 1/8. Also, to accelerate the convergence speed, Aitken's squared process is introduced to compute the axial induction factor only. It should be noted that the switch controlling the use of Aitken's squared process will be turned off whenever a singular value is encountered within that process.
The improved algorithm for solving the blade element momentum equations is outlined below with the modified part decorated with rectangular borders. Note that Spera's correction model is employed. Compute the relative inflow angle using Eq. 8 and effective angle of attack using Eq. 9 5 Look up and 6 Calculate and using Eqs. 3 and 4 7 Compute the thrust coefficient using Eq. 7 8 Compute tip loss factor F 9 if then 10 Compute a with Eq. 10 11 else 12 Compute a with Eq. 25 Compute a' using Eq. 9 15 end 16 end

Introduce the automatic relaxation factor decreasing technique Introduce Aitken's delta-squared process to compute a
In Table 3, the improved FPA will be applied to the previous simple test and the results are compared with traditional FPA with Spera's thrust correction. It is exciting to find the robustness of the FPA is greatly improved, especially for rf = 0.5 and 0.25. None of the failures case except for the second is found, which is removed later when rf is reduced. Considering the convergence speed and failure rate, the rf = 0.5 will be the best choice up to now. Also, it can be inferred that the automatic rf decreasing technique can enhance the robustness remarkably.

DIFFERENT METHODS
In this section, various performance tests have been performed to prove that the new thrust equation can greatly enhance the robustness, and Aitken's squared process can significantly strengthen the efficiency. Firstly, the same working space but with denser grids (20 × 20 × 20) will be used to examine the performance of three FPAs for 14 airfoils while keeping the rest parameters including F the same for the moment. Three FPAs are Meng's FPA, traditional FPA, and one-variable FPA proposed by Sun et al. 6 Note that it has already been shown in Section 4 that the under-relaxation factor can improve the robustness and the decreasing rf technique will be equipped for all the three FPAs. To indicate the superiority of the new thrust equation, Aitken's squared process will be deactivated. Then, the same test will be carried out for the improved FPM with or without Aitken's squared process equipped to reveal the benefit of Aitken's squared process. Lastly, thanks to the first two tests assume a unit F, the influence of tip loss factor F on the convergence performance will be investigated for the new FPA with all three modifications equipped.
A summary of the airfoils to be used is listed in Table 4. The Re is an abbreviation of the Reynolds number, which is a dimensionless quantity to measure the ratio of inertia forces to viscous forces. The Reynolds number is defined as the characteristic flow speed times the characteristic length over the kinematic viscosity. The chord length of the airfoil is adopted as the characteristic length. The incoming air velocity is adopted as the characteristic velocity.
The anterior three airfoils are S809 and S815, whose 2D airfoil characteristics between −20 degree and 40 degrees comes from. 20,21 The terms "clean" and "rough" represent the airfoil is tested with or without leading edge grit roughness, respectively. S809 is used for the wind turbine in the NREL Phase VI experiment. 22 The 3D airfoil data with rotational corrections are obtained and extended to full 360 degrees with the help of AirfoilPrep, 23 a tool developed by the National Renewable Energy Laboratory to prepare 3D airfoil data for AeroDyn. The posterior 5 FFA series airfoils together make up the noncylindrical part of DTU 10-MW reference wind turbine blade model. 18 Descriptions of how to obtain their 3D airfoil data are presented in Section 4. The last 6 airfoils, including one NACA airfoil and five DU series airfoils, comes from the NREL 5-MW reference wind turbine model 5 and the 3D airfoil data are given by Timmer. 24 The numbers of different failure cases for each airfoil are shown in Table 5. The new FPA and Sun' FPA outperform the traditional FPA in the success rates, which both give 0 failure rate for every airfoil tested. It can be found that failure case 1 is the exclusive failure type for the traditional FPA.

Relaxation factor
The decreasing rf methodology helps the traditional FPA and Sun's FPA to get rid of failure case 2 reported in Sections 3 and 4.
Despite the rather satisfactory robustness, it is disappointed to find that the time consumed by Sun's method is the most expensive. The overall failure rate and average times of iterations required to converge successfully for each airfoil are shown in Table 6. The failure rate is defined as the ratio of the number of failures to the total number of tests for a single airfoil, 8000. Only those which converge successfully  are counted in the computation of the number of average iteration. Note that the average numbers of iterations consumed are normalized by the corresponding ones obtained via Sun's method, which are bracketed in Table 6. Note that the "Avg." and "Iters." in Table 6 are abbreviations for "average" and "iterations." The convergence speed of Meng's FPA is only slightly boosted to the traditional one, which is acceptable considering the enormous improvement in the robustness. It is quite interesting that the convergence speed of the traditional FPA is faster than Sun's method, which seems to be contradictory to the results in the literature. 7 This discrepancy may be attributed to the choice of the convergence criteria. Remark that the is set to 1 × 10 − 6 by Sun et al. 7 The performance of Meng's FPA can further be enhanced by coupling with Aitken's delta-squared process, which is the alternative to Steffensen's method. Steffensen's method is a root-finding scheme with quadratic convergence, which avoids the difficulty to calculate the derivative in contrast to Newton's method. Comparisons of the convergence speed and failure rate for Meng's FPA with and without Aitken's delta-squared process under the same working space is given in Table 7. It can be easily seen that the average number of iterations has been cut down significantly without losing any robustness. Meng's FPA with Aitken's delta-squared process is even 5 times faster and can certainly save lots of computational costs in the BEM-related optimization in the future. It is worth mentioning that the new method is about 2 times faster than Ning's one-equation method, which is 11.3. 4 After coupled with Aitken's squared process, Meng's FPA turns to be quite promising. To fully understand the performance of the new FPA, the influence of tip loss factor F on convergence ability should be studied. Overall, five different F will be used. The convergence results will also be compared with those obtained by the traditional FPA, which are shown in Table 8. Here, the failure rates and numbers of iterations are obtained by averaging the results from all 14 airfoils. It is inspiring that Meng's FPA still gives 0 failure rate. As contrast, it becomes harder for the traditional FPA to reach a converged solution with decreasing F. The largest failure rate is found at F = 0.2. Besides, one common characteristic of the two FPAs is the average number of iterations goes up when F decreases. This can be partly explained that a smaller F tends to return a larger axial induction ratio when the rest parameters in Equations (5) and (25) are kept the same, which is detrimental to the convergence ability.

| CONCLUSIONS
The comprehensive classification of the operating state of the wind turbine gives a clear indication of the demarcation between each state. The thrust curve beyond the a > 1 is replaced by an extension of the high thrust correction of the turbulent wake state instead. It is straightforward that the resultant elimination of the singular point will have a positive effect on improving the robustness of the solution algorithm. In addition, three failure types of the traditional FPA are addressed and discussed with verification by a simple test. Three modifications are then applied to introduce the new FPA, including a new thrust equation, automatic relaxation factor decreasing technique, and Aitken's squared process. It has been shown the decreasing rf technique has a positive effect on improving the robustness via the test in Section 4. The superiority of the rest two modifications has been

Comput ations of the second critical induction ratio
There are two methods for computing the second critical induction ratio. The first one is based on the theoretical derivation, and the other one based on experimental observation. For the first method, the second critical induction ratio is considered to be the upper boundary presented by Wolkovitch 11 during analyses of the vortex ring state for vertical descent without tip loss. The descent velocity (V c ) is approximately −0.707 times the hover velocity (V h ), which is quite close to the ONERA VRS test result presented in the reference. 25 Taking advantage of the quartic fitting formulation of the induced velocity given by Leishman 16 and assuming that the induced power factor, which controls the effects of tip loss and nonuniform inflow, equals 1.25 based on work done by Leshiman, 16 it gives Then, the second induction ratio a crit2 based on theoretical derivation can be obtained: The minus sign is used here to in line with the axial induction ratio definition used in the wind turbine. For the other method, different from the theoretical derivation, Leishman 10 observed that the induced velocity curve obtained from flight tests and wind tunnel measurements agreed well with the one predicted based on the momentum theory for − 1.5 ≤