Next Article in Journal
Configuration Design and Dynamic Characteristics Analysis for Space Membrane Mechanism Based on Deployable Booms
Next Article in Special Issue
Implementation of a Holistic MCDM-Based Approach to Assess and Compare Aircraft, under the Prism of Sustainable Aviation
Previous Article in Journal
Analysis of the Effect of the Leading-Edge Vortex Structure on Unsteady Secondary Flow at the Endwall of a High-Lift Low-Pressure Turbine
Previous Article in Special Issue
Attitude Dynamics of Spinning Magnetic LEO/VLEO Satellites
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Analysis of Glauert Inflow Formula for Single-Rotor Helicopter in Steady-Level Flight below Stall-Flutter Limit

1
Section of Aerospace Engineering and Mechanics, Department of Military Mechanical Engineering, Military Academy, University of Defence in Belgrade, Veljka Lukica Kurjaka 33, 11042 Belgrade, Serbia
2
Department of Aerospace Engineering, Faculty of Mechanical Engineering, University of Belgrade, Kraljice Marije 16, 11120 Belgrade, Serbia
*
Author to whom correspondence should be addressed.
Aerospace 2023, 10(3), 238; https://doi.org/10.3390/aerospace10030238
Submission received: 7 January 2023 / Revised: 20 February 2023 / Accepted: 24 February 2023 / Published: 28 February 2023
(This article belongs to the Special Issue Advances in Aerospace Sciences and Technology III)

Abstract

:
This article addresses the numerical computation problem of induced inflow ratio based on the helicopter momentum theory in forward flight. The Glauert inflow formula (equation) is a nonlinear equation usually solved by the Newton–Raphson method in a relatively small number of iterations. However, many high-order convergence multipoint iterative methods have been developed over the last decade. The study examines several selected methods in terms of finding ones that provide a solution in only one iteration with acceptable accuracy. Furthermore, the influence of initial guesses on the accuracy of the obtained solutions has been investigated. In this regard, the practical range of parameters of the Glauert inflow equation for helicopters in forward flight is roughly determined by simplified modeling of a power and stall-flutter limitation. For these purposes, a basic low-fidelity longitudinal trim model of a single-rotor helicopter in steady-level flight is modified and numerically solved by a symbolic transformation of a system of 20+ nonlinear equations into a single nonlinear equation.

1. Introduction

Momentum theory is the first analytical tool for theoretically analyzing rotor thrust in forward flight conditions. It was developed back in 1926 by Hermann Glauert as a part of a broader, general theory for autogyro applications [1]. This basic theory represents a forerunner of the niches of today’s rotary-wing aerodynamic models [2]. Regardless of their complexity and primary purpose, all these models deal with a problem of induced velocity distribution at the rotor disc. In the early days, uniform distribution was routinely employed for rotorcraft performance estimation, real-time simulations, flight dynamics, and control applications. However, decades of strong computational development have enabled increasingly complex and more realistic aerodynamic models, including free-wake analysis, computational fluid dynamics, comprehensive tools for fluid-structure interaction, and dual-solver hybrid computational approaches [3]. Still, the current CPU/GPU performance level limits the use of CFD and CFD/CSD techniques in real-time applications. In this regard, when it comes to robust mathematical models applied in flight simulators, dynamic control, and aeroelastic stability evaluations, as well as preliminary performance estimations, nonuniform, both static and dynamic inflow models still play a dominant role [2,4]. This group of models includes linear ones whose induced velocity distribution is based on its main value at the rotor disc center, obtained by the momentum theory. Based on that, the problem of solving the Glauert inflow formula (GIF) is actual, and the momentum theory for vertical and forward flight is still a subject of improvement [5,6].
Momentum theory for forward flight is an analytical tool, and GIF, which relates forward flight speed, downwash flow, and thrust, has a complex analytical solution. Although the edgewise flight assumption leads to less accurate solutions for low advance ratios, it is still used in practice, especially for performance estimation during the preliminary design of the helicopter [7]. In the last more than a half-century, the numerical solution of GIF represents a standard approach in the rotorcraft community, and the classic Newton–Raphson (NR) method for solving nonlinear equations is commonly used. It is robust and has good local convergence. However, like all iterative one-point methods, it has a theoretical limitation established by Traub [8], who says that an iterative function ϕ x k that explicitly depends on the p 1 derivative of the function f x whose zero is required cannot attain an order of convergence higher than p . Consequently, the NR method has a convergence order p = 2 , which provides approximately two significant digits in one iteration.
Increasing the order of convergence of one-point methods is only possible by increasing the number of derivatives. However, it limits computational efficiency, which has resulted in decades of stagnation in this area [9]. The rapid development of computer hardware, symbolic algebra software, and high-precision arithmetic leads to increased interest in multipoint iterative methods because of their computational efficiency, which is much greater than the efficiency of one-point methods [10]. During the first decade of the 21st century, at least 200 new multipoint methods have been developed, and that number continues to increase [11].
On the other hand, Newton-based methods are also being developed and used to solve various engineering problems. Quasi-Newton methods are used in machine learning, neural network training, fluid-structure interaction, and other engineering areas [12,13,14,15]. Recently, the NR method has been used in coupling with other iterative algorithms as an alternative to the traditional finite element method for large multiphysics problems [16]. Jacobian-free NR-based methods for solving multiphysics nonlinear systems are also being improved [17].
Although powerful numerical methods have been generated in the last 30 years, the rotorcraft industry and science still dominantly use a method that originated more than 300 years ago. It is indisputable that the NR method has a quality that keeps it alone at the top of the list of used numerical methods for solving nonlinear equations. Still, there is a need to examine the utility value of newly developed ones. If the GIF is solved numerically by a novel multipoint method, that would primarily increase the computational efficiency. In addition, if an acceptable solution accuracy is reached in just one iteration step, that would increase the number of users preferring classical approximations rather than numerical methods. Finally, novel dual-solver hybrid approaches reopen doors for simple theories, such as momentum theory, to be used in a couple with high-fidelity models, which gives additional scientific importance to improving the GIF solution procedure.
It is essential to point out that a mathematical analysis of GIF cannot be found in the literature [18]. This paper aims to change that fact and use breakthroughs in symbolic and numerical mathematics to improve solving procedure of the low-fidelity model of single-rotor helicopter longitudinal force equilibrium in forward flight, which is necessary to use to determine the parameters of the GIF equation.

2. GIF Formulation and Solution Procedure

According to the momentum theory for forward flight (Figure 1), the GIF is a transcendent fixed-point equation for the induced inflow ratio:
λ i = C T 2 μ 2 + μ z + λ i 2 ,
which can be expressed in algebraic irrational zero-point form as
λ i C T 2 μ 2 + μ z + λ i 2 = 0 ,
and in rational (normal) form as follows
4 λ i 4 + 8 μ tan α D λ i 3 + 4 μ 2 tan 2 α D + 1 λ i 2 C T 2 = 0 .
Instead of the induced inflow ratio λ i , the rotor inflow ratio λ can be used as a variable in Equation (1):
λ = λ i + μ z = C T 2 μ 2 + λ 2 + μ tan α D
with the algebraic zero-point irrational form given as
λ C T 2 μ 2 + λ 2 μ tan α D = 0 ,
and rational form:
4 λ 4 8 μ tan α D λ 3 + 4 μ 2 tan 2 α D + 1 λ 2 - 8 μ 3 tan α D λ + 4 μ 4 tan 2 α D C T 2 = 0 .
Equations (3) and (6) are quartic algebraic equations that can be represented as
A 4 x 4 + A 3 x 3 + A 2 x 2 + A 1 x + A 0 = 0 .
Constants A 0 ,   A 1 ,   A 2 ,   A 3 ,   and   A 4 are given in Table 1.
Equations (2) and (5) are irrational because they have an unknown in the radicand. Equations (3) and (6) represent their normal (rational) forms, respectively. Although no recommendation can be found in the literature dealing with the question of which GIF form is the best for numerical solving, there is a clear reason for preferring an irrational form. Namely, in addition to the roots of the original (irrational) form, the normal form may contain some superfluous roots that are not solutions of the original form [19]. On the other hand, the rational form is more suitable for calculating derivatives, which can be a significant factor if a one-point method of a high order of convergence is engaged to solve the equation. There are four forms of the GIF, which can be solved numerically, and Equations (2), (3), (5), and (6) represent them.
Procedures for numerical solving of GIF can be found in [20] (pp. 97–98) and [21] (pp. 126, 579–580). Both authors use the NR method to solve an irrational form for λ (Equation (5)) and state that three-four iterations are sufficient in most cases. Leishman [20] suggests two numerical methods with practical engineering use: a simple fixed-point iteration (FPI) and NR method, where FPI typically requires 10–15 iterations, but in some cases, especially at lower advance ratios, a larger number of iterations may be necessary. For stopping criteria, the author uses an error estimator as follows:
ε = λ n + 1 λ n λ n + 1 ,
and states that when ε 0.0005 convergence is said to occur. According to Leishman, the initial value for λ 0 is usually the inflow ratio for hover:
λ h = C T / 2 .
W. Johnson [21,22] recommends the initial value for λ 0 as
λ w j = λ h 2 λ h + μ z 2 + μ 2 + μ z .
In addition, the author proposes using a relaxation factor r f = 0.5 for convergence improvement:
λ n + 1 = λ n r f f λ f λ .
The initial value given in Equation (10) can be modified to become even closer to the root as
λ m d = λ h 2 λ e + μ z 2 + μ 2 + μ z ,
where λ e is the main inflow ratio for edgewise flight ( α D = 0 ):
λ e = 1 2 μ 4 + 4 λ h 4 μ 2 .
Next to the numerical solution for λ , the GIF can also be solved for λ i . However, in the literature, no analysis compares the number of iterations and other parameters of the numerical solution for those two variables. Therefore, six initial guesses (values) for λ 0 are considered in this study, as it is shown in Table 2.
A few years ago, GIF’s analytical, closed-form solution was presented and offered to the rotorcraft community as a benchmark solution in helicopter aerodynamics for comparing numerical results [18]. It is a relatively new, easily programmable solution based on linear algebra and symbolic mathematics, with the potential to be used instead of numerical solving. Since the author examined only one case ( μ = 0.1 , C T = 0.003 , α D = 4 ° ) with the primary pedagogical goal of demonstrating the power of novel mathematical tools, the analytical solution is not taken as a benchmark solution. However, it is evaluated, along with numerical methods.

3. Models and Methods

GIF has only three parameters: advance ratio μ , thrust coefficient C T , and rotor disc incidence angle (angle of attack) α D . These parameters are mutually determined, and their ranges need to be defined. Only for this research, the union of all μ , C T , and α D combinations that may occur in steady-level flight conditions are named the Acceptable Parameter Range (APR). It is a fictive union of all steady-level flight regimes of an infinite number of single-rotor helicopters that may exist in practice.
Accordingly, two tasks are imposed:
  • Definition of a method for determining the APR;
  • Presentation of the APR with a finite number of cases.
Precise determination of APR boundaries is a complex and time-consuming process requiring sophisticated and computationally demanding tools. The main goal is to exclude from the numerical analysis as many cases as possible that have no practical significance and to ensure that a potential conclusion about the use of a particular numerical method in solving the Glauert inflow equation refers to the widest possible ranges of steady-level flight conditions, take-off masses, geometric and aerodynamic characteristics of helicopters. Therefore, it is enough to use a robust, fast, and sufficiently reliable tool capable of providing that the APR boundaries it sets encompass all possible cases. For that reason, the following steps are adopted:
  • Selection of a finite number of cases for APR boundaries determination purposes;
  • Establishing the APR boundaries as C T m i n = f 1 μ , C T m a x = f 2 μ , α D m i n = f 4 μ and α D m a x = f 4 μ ;
  • Establishing a more expansive APR zone with an acceptable safety margin;
  • Determination of the finite number of possible cases as APR zone representatives for numerical methods evaluation purposes.
The range of possible cases is very wide. There are many ultralight RW-UAV applications and a few single-rotor helicopters with a maximum take-off mass (MTOM) below 300 kg. The MTOM of the largest helicopter in the world is greater than 50,000 kg. The Official World Record in the category Rotorcraft–Altitude in horizontal flight set above ISA troposphere boundary (11,010 m) by sky-crane helicopter Sikorsky CH-54B on 4 November 1971. The fastest helicopter in the world is the British Lynx, with an achieved forward flight speed of 400.87 km/h. The Lynx held this title for 37 years.
Considering the high variability in rotorcraft take-off mass, flight speed, altitude, and other characteristics, around 7 million points are selected and divided into four working groups (WG), as shown in Table 3. The selection process is based on past and current design trends [21,22,23,24,25,26]. Next, the parameters a, b, and k are used to determine helicopter main and tail rotor diameter, angular velocity, and equivalent helicopter drag/lift area by statistics. Finally, the control group (CG) is used for grid density validation.
The first task in the APR establishing process is to develop an algorithm for GIF parameter determination. Considering the number of calculations that must be performed to determine the APR zone is measured in millions, a simplified approach to the problem of helicopter-trimmed flight is necessary. In this regard, only the longitudinal force equilibrium model (LFEM) is used, as shown in Figure 2.
LFEM is based on a simple low-fidelity longitudinal trim model given in [27] (pp. 116–118). The longitudinal force equilibrium equations for steady-level flight are:
T sin ( α D ) + L F sin ( Δ α D ) H c cos 2 ( α D ) D cos ( Δ α D ) H T R = 0 ,
and
T cos ( α D ) + H c sin ( α D ) cos ( α D ) W o g L f cos ( Δ α D ) D sin ( Δ α D ) = 0 .
The lack of complete trim equations (longitudinal moment balance and lateral force and moment balance) results in a lower quality of the obtained solutions, which is compensated by increasing the safety margin of the established APR boundaries. Because of demanding computational requirements, flapping effects, rotor imperfections, and nonlinearities are also neglected.
On the other hand, usual small angle assumption, typical for simplified low-fidelity models, cannot be applied, which significantly complicates the mathematical model and increases the number of nonlinear equations that must be solved simultaneously. In addition, airframe lift force and tail rotor drag force are addressed. The solution is obtained in two steps: LFEM transformation into a single longitudinal force equilibrium equation (LFE) and numerical solution of LFE using the NR method. LFEM and its numerical solution procedure are described in detail in Appendix A.

3.1. Available Power Limit

Many real limitations may be included in the LFEM model. The most important are available power and dynamic stall (stall-flutter) limit. Available power is strongly limited and is represented by Equation (16) which is obtained by statistical data from [23]:
P T O = 0.0764 W M T O 1.1455 ,   ε a v = 14 % ,     ε M A X = 37 % ,     R c = 0.989 ,     W M T O kg .
Since the maximum errors have always been on the side of overestimating the available power, the following expression is adopted for the maximum take-off power in [kW]:
P T O = 87.096 W M T O 1.1455 .
For altitudes up to 4000 m, the adopted model gives that the maximum available power is 95% of the output engine(s) power. The remaining 5% covers the power needs of the tail rotor, other consumers, and transmission losses. Furthermore, the novel hybrid and fully electrical power units were not observed because of their much lower gravimetric power density than hydrocarbon fuel-based power plants. For higher altitudes, available power is represented as a function of static pressure and temperature ratio:
P A h > 4000 m = 0.95 P T O δ θ .
On the other hand, the required power P R has been calculated using the following expression:
P R = C p ρ π R 2 V t 3 ,
where
C p = k i λ i C T + σ c d 8 F p + 1 2 f e π R 2 μ c 3 .
The empirical factor k i considers blade tip losses, nonuniform inflow, and other rotor non-perfections, and is given in Table 4. The blade profile drag empirical factor F p with an error of less than 1% for all speed ranges was developed in [28] as
F p = 1 + μ c 2 1 + 5 2 μ c 2 + 3 8 μ 2 4 + 7 μ c 2 + 4 μ c 4 1 + μ c 2 2 9 16 μ 4 1 + μ c 2 + 3 2 μ z 4 + 3 2 μ z 2 μ 2 + 9 16 μ 4 ln 1 + 1 + μ c 2 μ c ,
and it is used in this study. In addition, as the small angle assumption is not employed, μ c is used instead of μ in Equation (20).

3.2. Stall Flutter Limit

Dynamic stall is a phenomenon that allows the local angle of attack to grow beyond the static stall limit. Such conditions include high drag penalties, negative pitch damping, and amplified vibrations. Dynamic stall criteria usually represent a limit beyond which forward flight is unacceptable due to excessive pitch link vibrations, although the rotorcraft is usually still flyable. Many empirical dynamic stall models, such as the Boeing, ONERA, and the Leishman–Beddoes model, have been developed. These models are primarily used in complex rotary-wings modeling. On the other hand, F. Harris [29] developed simple stall criteria, which can be implemented in low-fidelity models:
C T σ = C l M A X 6 1 μ 2 + 9 4 μ 4 1 + 8 3 μ + 3 2 μ 2
In this study, a simple assumption is made that the dynamic stall limit on the retreating blade is 50% greater than the static stall limit, and a slightly modified Harris criterion for a retreating blade stall at 0.75 R is used as
C l r b < C l M A X r b , d y n a m i c ,
where
C l r b = C l a v 1 + 8 3 μ + 3 μ 2 + 4 μ 3 + 9 4 μ 4 1 μ 2 + 9 4 μ 4
The maximum lift coefficient at 0.75 R of the retreating blade depends on local Reynolds and Mach numbers. For this study, the blade airfoil maximum lift capabilities are based on NACA 0012 airfoil, and the influence of the local Reynolds number is modeled as
C l M A X r b , s t a t i c ,   M r b < 0.3 = 67.57 34.07 log R e l o c + 5.722 log 2 R e l o c 0.3141 log 3 R e l o c ,
where the local Reynolds number is defined as
R e l o c = 0.75 V t V cos α D · c ν
Equation (25) is obtained from [30] (Figure 13) by the GraphPad PRISM® software tool for nonlinear regression analysis and is used if 3 × 10 5 R e l o c 8 × 10 6 . The maximum lift coefficient on the retreating blade for lower Reynolds numbers is set to 1.0 and for higher ones to 1.75. Furthermore, the influence of the Mach number is also obtained from [30] (Figure 14) similarly:
C l M A X r b , s t a t i c ,     M r b > 0.3 = C l M A X r b ,     s t a t i c 0.8 M r + 1.1 C l M A X r b ,     s t a t i c 0.8
The maximum static lift coefficient C l M A X a v ,   s t a t , can be calculated by Equations (25) and (27) where M and R e are used instead of M r and R e l o c .

3.3. Numerical Methods and Procedures

This subsection presents two computational plans: a plan for APR boundary determination and a GIF numerical solution with certain numerical methods. In this regard, The Numerical Computation Plan 1 for APR determination is shown in Table 5. There are seven cycles of calculations: 6 cycles refer to 4 working groups (WG 1–4), and the last refers to 4 control groups (CG 1–4), as presented in Table 3. The control groups are introduced to check the influence of the grid density on the obtained results. In total, over 2.5 billion LFE calculations were performed. This demanding and computationally expensive task has been completed on a workstation with 20 processor cores. The computational code is written in the Julia® programming language because it is one of the best solutions for large-scale numerical work [31].
As can be seen in Table 5, there are three Power limits models. The Power 1 model is described in the previous subsection. Power 3 is practically the same, with the difference that the small angle assumption is adopted, which means that Equation (20) contains μ . Finally, Power 2 model uses the classic Gessow and Crim approach to the compressibility effect on the blade profile drag. The main task of this model is to serve as a benchmark for the Power 1 limit model.
In this paper, GIF is first numerically solved by a three-point eight-order Sharma-Sharma (Sharma) method using multi-precision numerical computing. The Sharma method is robust and helpful for single and multiple root-finding problems that require only four functional evaluations [32]. A high level of accuracy is obtained by the Arb.Numerics Julia® package for arbitrary precision, which is set to 200 of the 2500 available significant digits. Numerical computing of GIF in this research is performed with the:
  • Standard arithmetic precision;
  • Arbitrary arithmetic precisions (only for the benchmark solution);
  • Initial values are taken from Table (2);
  • Stopping criteria   E o < ε ( ε = 5 × 10 n ,   n = 4 ,   6 ,   8 ,   10 ,   12 ,   16 ) or “one iteration”;
  • Stopping criteria   E o < 5 × 10 150 (only for the benchmark solution);
  • Several one-point and multipoint iterative methods without memory for solving nonlinear equations.
The methods used are shown in Table 6 and thoroughly presented in [9,10,11]. As seen, 14 methods, 8 stopping criteria, 4 GIF forms, and 6 initial values were used, and 1344 numerical cycles were performed. Each cycle consists of about 3 million individual cases. For this work, GIF has been solved numerically about four billion times. The main goal of these numerical calculations is to obtain the following:
  • Maximum, minimum, and average iteration number (ITNumb) as a function of μ for the entire APR.
  • Maximum real relative error (Et) for all examined cases;
  • Percentage of cases of failed convergence, if any;
  • Presence of non-physical roots;
  • Quality assessment of used initial values;
  • Quality assessment of iterative methods.
In addition to the above, it is essential to point out that each numerical method used in this research is written in the Julia® programming language in such a way that the iterative sequence will be interrupted in the case the stopping criteria are satisfied or the number of iterations is greater than 1000 (failed convergence) or the solution is a NaN, zero or infinite. While the acceptable solution is practically impossible to obtain in the case of failed convergence, in other cases (NaN, zero, or infinitive), the value from the previous iteration may be taken as a solution. Additionally, for multipoint methods, a solution from the previous step (point) of the same iteration may be delivered as a final. Consequently, if some multipoint method is used, the number of iterations (ITNumb) does not necessarily have to be an integer.

4. Results and Discussion

Firstly, the results of the APR determination are shown and discussed. Then, the following analysis of the performance of the selected numerical methods for solving the GIF whose parameters are in the APR is presented.

4.1. APR Determination

Figure 3 (left) shows the maximum and minimum values of the angle of attack and thrust coefficient for the basic LFE model (No limits) and five additional models with limits. As can be seen, maximum α D values are restricted by Power 1 and Power 2 models. Furthermore, the matching of minimum and maximum values of C t and α D which are obtained by these two models is excellent. However, the results obtained by Power 3 model compared to the previous two models deviate significantly in obtaining very large C t and α D maximums for 0.05 μ 0.27 and 0.04 μ 0.36, respectively. Minimum values of these parameters are very similar for all three Power models, while in comparison to them, the other three models give lower C T m i n .
Regarding C T m a x , for an advance ratio less than 0.2, the limit is provided by Power 1 model. In contrast, the emergence of dynamic stalling represents the main limiting factor for higher advance ratios. According to the above, the results obtained using Dynamic stall and Power 1 limits served as the basis for establishing the APR boundaries, shown in Figure 3 (middle). Because the very simple longitudinal force equilibrium model is used as an APR solver, a more expansive APR zone is established, as shown in Figure 3 (right). Grid density is verified by the control group, with about 97 million points. As presented in Figure S1, there are no significant differences between the working and control groups.
The main goal of LFEM is to determine the APR zone, and in this regard, a dynamic stall (stall flutter) limit is fundamental, especially for an advance ratio higher than 0.2. Unlike fixed-wing aircraft, rotorcraft normally fly in mild and medium stall flutter conditions because the blades very quickly leave the retreating side of the rotor. The main consequence of stall flutter is a pitch link vibratory load often used for practical stall boundary determination. As already mentioned, a modified Harris criterion (MHC) is used in this study. Figure 4a shows MHC obtained by NCP 1 execution compared to standard Harris criteria.
The data from [33] are used and compared with the MHC line adapted to the given helicopter model and test flight conditions for validation. As shown in Figure 4b, there is a relatively good matching of MHC with the pitch link vibratory limit. The advantage of such a modification on the Harris criterion is that it allows variability of the maximum lift coefficient depending on the local (representative) Reynolds and Mach numbers on the retreating blade, which is the place of stall flutter occurrence.

4.2. GIF Numerical Analysis

The results of the NCP 2 execution are given in this subsection. As previously stated, 1344 cycles were conducted, and each cycle had N t o t . = 4 , 375 , 561 points. Given that one cycle requires solving the equation twice, GIF has been solved almost 11.8 billion times for regular and benchmark solutions. Due to the large number of obtained data, only the most important results are presented and discussed.

4.2.1. Percentage of Failed Convergence (PoF)

Iteration and failed convergence counters are built into the Julia® code for all used numerical methods. The number of iterations greater than 1000 is empirically chosen to stop the solver and to indicate the occurrence of failed convergence. Firstly, numerical solutions of Equation (6) are examined. According to the results shown in Table 7, the presence of failed convergence correlates with stopping criteria (required accuracy). Therefore, eliminating failed convergence presence by reducing the required accuracy is possible with all methods except the BWR II. In general, a small number of cases of failed convergence ( NoF ) concerning the total number of performed calculations N t o t . . The percentage of failed convergence ( P o F ) is calculated as
P o F = 100 N o F N t o t .   % ,
and its distribution over the advance ratio is shown in Figure 5. It can be seen that PoF increases with the advance ratio and reaches a value of around 50 percent for the Chun and Cordero methods. In addition, an initial value does not generally significantly influence the PoF. Finally, it is essential to note that using arbitrary precision in the Sharma method to calculate benchmark solutions has completely eradicated convergence failures. Since the double-precision floating-point arithmetic is set by default for real numbers in 64-bit system architecture (Float64), the maximum number of significant digits is about 15 to 17. Therefore, the accumulation of computational error is the probable cause of these failures.
To confirm this assumption, the case of the largest percentage of failed convergence from Table 7 ( PoF = 14.13 % ,   ε = 5 × 10 16 ,     λ o = λ w j , Cordero method) was taken and recalculated with appreciably increased arithmetic precision. The ArbNumerics Julia® package was used with 200 significant digits. As a result, the convergence was completely successful ( PoF = 0 % ). The same result was achieved using the BWR II method, where the failures were present for all stopping criteria.
Equation (3) has much better properties regarding convergence failures, as shown in Table 8. One-point methods are not presented in the table because no failures were detected for ε < 1 × 10 17 . When it comes to multipoint methods, there are certainly cases of convergence failures, but only if ε 1 × 10 16 . Equations (2) and (5) express excellent properties: no cases of failed convergence were found for declared precision.

4.2.2. Number of Iterations (ITNumb)

As in the previous subsection, Equation (6) was analyzed first. Figure 6 shows the number of iterations of one-point methods for an initial guess λ 0 = λ h and stopping criteria ε = 5 × 10 14 . The maximum number of iterations is relatively large and reaches over 150 for the E5 method and over 15 iterations for the Halley and Kiss methods. The minimum number of iterations ranges from two to three. The average I T N u m b is about four iterations for all one-point methods and is represented by dashed lines on two diagrams in Figure 6.
The total number of cases in which the required number of iterations is greater than 10 ranges from 0.02 to 4.15 percent, and three to five iterations are sufficient in 87.5 to 97.5% of cases, depending on the chosen one-point method. The results (ITNumb vs. advance ratio) for other stopping criteria ( ε = 5 × 10 n ,   n = 4 ,   6 ,   8 ,   10 ,   12 ,   16 ) and for multipoint methods (where λ 0 = λ h ) are very similar to the above (Figure 6).
The required number of iterations can be drastically reduced if λ w j or λ m d are used instead of λ h , as shown in Figure 7. The stopping criteria are essential because a small ε (high-required accuracy) causes a larger number of iterations. As can be seen in Figure S2, only one iteration in about 75% and up to two iterations in 99.8% of cases is sufficient if the Sharma method is used in conjunction with λ 0 = λ m d , where stopping criteria are set on ε = 5 × 10 14 .
Even better results are obtained if the irrational form is solved. Namely, under the same conditions, Equation (5) is solved with the Sharma method in two or fewer iterations in all cases. Figures S3–S8 show an average number of iterations and their percentage distribution of remaining GIF forms. As can be seen, the use of improved, i.e., modified initial value ( λ m d ) significantly reduces the required number of iterations, regardless of the stopping criteria, the type of method used, and the form of the Equation being solved.

4.2.3. Relative Errors ( E o , E t )

Let λ H P S be the benchmark solution of GIF obtained using high precision with 200 significant digits. Then, E t is a relative error with respect to the benchmark solution
E t = λ f λ H P S λ H P S ,
where λ f is a final numerical solution reached using a standard computational precision of up to approximately 16 significant digits. As a reminder, the relative input error that was used as a stopping criterion in the iterative cycle is marked with ε and defined by Equation (8). It represents the maximum acceptable relative distance between two adjacent solutions in the iterative sequence. The relative output error is defined as
E o = λ f λ f i λ f ,
where λ f i is a solution from the previous iteration ( i = 1 ) or the previous point of the same iteration step ( i 1 / 3 ,   1 / 2 ,   2 / 3 , applicable only to multipoint methods).
According to the above, all analyzed numerical methods come to the solution in two ways:
  • E o < ε —an iterative process is accomplished, convergence is achieved, the input criterion is fulfilled, and the output solution is λ f ;
  • E o ε —the stopping criterion is not activated, but the iterative process is interrupted due to the failed convergence or due to an inappropriate solution (NaN, infinite, or zero). Then, λ f i becomes a final numerical solution λ f .
  • The first case ( E o < ε ) can be seen in the following diagrams in Figure 8;
  • Equation (6), all methods, all initial values;
  • Equations (2), (3), and (5), one-point methods, all initial values;
  • Equation (3), two-point methods, all initial values.
Anomalies can be observed on other charts ( E o ε ). The occurrence of failed convergence has already been discussed, and a comparative analysis of data shown in Table 7 and Table 8 with the charts in Figure 8 leads to the two groups of anomalies whose cause is different from convergence failure. The first group refers to Equation (6) solved by:
  • Chun and KT Diff, which gives E o m a x 1 for all stopping criteria where λ 0 = λ h ;
  • KT Diff with λ 0 λ w j , λ m d , which gives E o m a x 2.5 × 10 8   ε < 10 8 ;
  • Sharma with λ 0 λ h , λ w j , λ m d , which gives E o m a x 2.5 × 10 8 for ε < 10 8 .
The second group includes Equations (3) and (5) solved by the Sharma method with λ 0 λ h ,     λ w j ,   λ m d , which also give E o m a x > ε , as shown in Figure 8. The most common cause of anomalies is NaN, which does not necessarily represent an unfavorable outcome. For example, when finding a benchmark solution, NaN occurs in 100% of cases, but λ f i is already a sufficiently high-quality solution after 1 + i iterations, and in that case, E o m a x 5 × 10 51 . The total maximum varies with the change in advance ratio, as shown in Figure 9. The same applies to the spectrum of E o m a x values given in Figure 8, while their dependence on the advance ratio is shown in Figures S9–S12.
Given that the subject of the study is the examination of numerical methods with standard accuracy (up to 16–17 digits), it is much more essential to analyze maximum relative errors compared to the benchmark solution ( E t m a x ), which are given in Figure 10 and Figure 11. Figure 10 shows E t m a x for all methods and initial values used in solving Equations (2), (3), and (5). As can be seen, the actual minimal accuracy (the maximal error) E t m a x is far better compared to the stopping criteria ε and the maximal output relative error E o m a x . However, when it comes to Equations (3) and (5), it is noted that E t ε but only for the highest required precision ( ε = 5 × 10 16 ). Furthermore, the lowest accuracy is detected for the Sharma method used with λ 0 = λ m d in all three equations, which is explained by the small number of iterations and the high order of convergence of the method.
Figure 11 shows E t m a x for Equation (6). These results are presented separately from Figure 10 because they require a more detailed analysis. Namely, all one-point methods used in conjunction with λ h indicate that E t m a x 0.5 regardless of stopping criteria. The same applies to all two-point methods except the Chun and the three-point Sharma and BWR methods. With omitted methods, the results are even worse (extremely high values of E t m a x ), but they are easily explained. It is enough to look at the corresponding diagrams in Figure 8 and see that anomalies ( E o m a x ε ) have already been detected there.
Additionally, BWR II methods have a high maximal actual relative error ( 4 × 10 2 < E o m a x < 5 × 10 8 ) for all initial values, but it is clear that this phenomenon is a consequence of failed convergence. Completely satisfied stopping criteria with excellent output relative error on the one side and disappointingly large actual relative error on the other can only be explained only by a successful convergence to superfluous roots, which is the case.
Further analysis revealed superfluous roots are always at a distance   Δ λ 0.0035 from the physical solution. This fact enabled the calculation of the probability of convergence towards the “wrong” solution during solving Equation (6). Figure 12 shows that probability increases with the advance ratio and reaches over 60 percent. It should be noted again that this applies only to λ 0 = λ h , which is essential considering that using λ h as an initial value is standard practice. Equation (3) also represents a rational form of the inflow equation, but the “superfluous roots” problem was not detected.

4.2.4. One Iteration Solution

A solution obtained in just one iteration can be considered noniterative in the practical sense because stopping criteria are not used. Figure 13 and Figure 14 show maximal relative error with respect to the benchmark solution, depending on the advanced ratio μ for Equations (5), (6) and (2), (3), respectively. After analyzing these charts, several findings can be highlighted:
  • Convergence towards an extraneous root may occur if λ h is used as an initial value to solve Equation (6), regardless of chosen method;
  • The solutions of Equation (3) obtained in one iteration with λ o = λ h has no practical significance regardless of the chosen method due to the very large maximal relative error concerning the benchmark solution ( E t m a x );
  • Failed convergence may occur if the BWR II method and λ m d initial value are used to solve Equation (6) with standard arithmetic precision;
  • Some anomalies in E t m a x magnitude may appear if multipoint methods and λ m d are used together to solve all forms of GIF except Equation (2), where no anomalies;
  • E t   M a x < 5.3 × 10 16 is certainly reached over the entire range of advance ratio if the Sharma or KT Diff method is used in conjunction with λ m d to solve Equation (2) which applies to standard arithmetic precision architecture. Figure S13 shows that the above conclusion applies only to the APR zone, i.e., that the numerical methods have worse performance outside the above zone;
  • E5 one-point method used with λ o = λ m d to solve Equations (2) and (3) provides high accuracy of obtained solution where E t   M a x < 9.7 × 10 10 for low advance ratios and E t   M a x < 2.5 × 10 14 for μ > 0.23 .
The previously mentioned failed convergence, which occurs if the BWR II/ λ m d pair is engaged in solving Equation (6) in standard precision, can be eliminated by using arbitrary precision with 200 significant digits, as shown in Figure 15. The same applies to the anomalies in E t m a x magnitude that arise when some analyzed multipoint methods are used with λ o = λ m d for numerical solving of all GIF forms except Equation (2), shown in Figure 15. Furthermore, it can be concluded that the accumulation of computational error at standard precision affects the final accuracy of the obtained solution even after one iteration. The situation is even more unfavorable if multipoint methods are used in several iterations, which could already be seen in the previous subsection. The reason lies in the large difference between E o and E t after the first iteration (Figure 16a), which practically unnecessarily activates the second one. As a result, calculation error accumulates, and the final solution has less precision than one obtained after the first iteration.
For example, the numerical solving of Equation (5) by the Sharma method with an initial value λ o = λ h and stopping criteria ε 5 × 10 8 always gives ITNumb = 1.67 and occurrence of NaN in 100% of cases, which means that the stopping criteria are not fulfilled, as is seen in Figure 8.
Computational error accumulates even during the first iteration, especially if the high order of convergence methods is used with the most precise initial value ( λ m d ). In that regard, Figure 16b gives a comparative view of the maximal relative errors after the first iteration and after 2/3 of the first iteration for Sharma/ λ m d duo used in solving Equations (3), (5), and (6). As can be seen, better precision is obtained after the first two steps compared to the whole iteration.

4.2.5. Newton–Raphson Method

The NR method, as the most commonly used tool for solving nonlinear equations, is analyzed separately due to its outstanding historical significance. It has been applied to solve all GIF forms with three initial values ( λ h ,   λ w j ,   λ m d ), with or without the use of the convergence acceleration factor ( f = 0.5 ) recommended by W. Johnson in [21]. Regarding failed convergence, the NR method behaves similarly to the previously discussed one-point methods: it only occurs with Equation (6) in about two percent of examined cases. The relaxation factor reduces the number of convergence failures by about 7.5 times, as seen in Table S1. On the other side, using the relaxation factor f = 0.5 significantly increases the number of iterations, which is presented in Figure 17. The maximum values of the required number of iterations are higher than the average and are shown in Figure S14. The maximum number of iterations of the NR method is at the level of examined one-point methods or even slightly higher. At the same time, the introduction of the relaxation factor significantly increases the number of required iterations. It does not apply to the case of solving Equation (6) with λ h , where the maximum number of iterations is large for all methods.
Figure 18 shows the change in E t M a x depending on the stopping criteria ε for all GIF forms. As in previous cases, convergence towards extraneous roots is noticeable for Equation (6) and initial value λ h . The positive influence of the relaxation factor f = 0.5 in increasing the accuracy of the numerical solution is not detected. The relative error is higher or at least the same if the mentioned factor is used. In general, the relative accuracy of the numerical solution achieved by the NR method is at the level of the previously analyzed one-point methods.
The solution obtained in one iteration is shown in Figure 19. If the results are compared with the diagrams given in Figure 13a and Figure 14a, it can be concluded that the NR method does not perform better than other one-point methods in case of use in one iteration. In addition, Figure 19 shows the negative influence of the relaxation factor on the achieved accuracy.

4.3. GIF Analytical Solution

According to [18], an analytical solution gives four roots: one complex conjugate pair and two real roots. Unlike [18], where only one case (point) was considered, in this study, a closed-form solution is tested over the entire APR zone in 436,067,402 points, and the observations are as follows:
  • The first two roots form a complex conjugate pair in 100% of tested cases and can be completely rejected;
  • The third and fourth roots form a complex conjugate pair in 2.8% of tested points;
  • The third and fourth roots are a NaN in 1.25% of cases;
  • The third root is a negative real number in 49.6% of cases;
  • If it is not a NaN or complex number, the fourth root represents a physical solution;
  • The analytical solution is found in about 95% of the tested cases.
The accuracy of the analytical solution is shown in Figure 20. As can be seen, it has a much larger relative error, minimum and maximum, compared to very precise, previously analyzed numerical solutions obtained in only one iteration using the Sharma and E5 methods with an initial value λ 0 = λ m d for solving Equation (2). The above, as well as the fact that a solution cannot be found in about five percent of cases, leads to the conclusion that closed-form solution does not show sufficient quality in comparison with numerical methods.

4.4. GIF Computational Cost

In order to complete the comparison of multipoint with single-point methods, it is necessary to evaluate the computational cost of the applied methods. Furthermore, verifying whether the multipoint method will achieve its computational efficiency in practice is necessary. In a limited numerical experiment applied to the rational forms of GIF, by comparing the time required for identical calculations on a sample of over 21.4 million different combinations of GIF parameters, the comparative advantage of the three-point methods was confirmed. Introducing an additional four-point method (order of convergence equal to 16) showed that these methods have the weakest performance of all analyzed methods. Therefore, three-point methods represent the optimal choice. In addition, introducing the relaxation factor in the NR method increases the required calculation time compared to the standard method by about 50 percent. The selected initial guess and stopping criteria have an impact on the required calculation time, but to a much lesser extent than the choice of the method itself. The KT Diff method stands out for its calculation speed. The obtained results are shown in Table S2 and given in the Supplementary Material (https://github.com/marjanvtsl/GIF-Numerical-Analysis.git, accessed on 23 February 2023).

5. Conclusions and Recommendations

Symbolic algebra, arbitrary-precision arithmetics, and multipoint numerical methods are proficient mathematical tools used in engineering practice. They allow old, entrenched methods to be replaced by more modern, faster, and efficient ones. The low-fidelity longitudinal force equilibrium model of a single-rotor helicopter in forward flight is modified (no small angle assumption) and represented by one nonlinear variable equation instead of a system of nonlinear equations, which significantly simplifies the numerical solution procedure. The used algorithm allows changing the variant of the equation during the iterative cycle. Further practical implementation of the LFEM requires complete trim equations to be included in the model.
Regarding the numerical analysis of the GIF solution, the conclusions are as follows:
  • The use of λ w j for the numerical solution of the Glauert inflow equation results in a smaller number of iterations and a higher accuracy of the obtained solution compared to the case when the induced ratio in hover λ h is used as an initial guess;
  • Modifying the previously mentioned initial value ( λ m d ) gives an even better result: the smallest average number of iterations and the highest precision of the solution regardless of the used method and stopping criteria;
  • The chosen form of the Glauert inflow equation (rational or irrational, with unknown λ or λ i ) directly affects the performance of the numerical solution where the best convergence failure-free performance gives Equation (2);
  • Optimal three-point methods used along with λ m d and high accuracy requested often lead to anomalies that significantly degrade the performance of the used method, and in this regard, it is much better to use them in only one iteration procedure for standard precision arithmetic;
  • Multi-precision arithmetic is computationally expensive, but it completely frees the solutions from the accompanying anomalies;
  • Analytical solutions in the standard arithmetic precision environment perform much worse than numerical solutions;
  • The Sharma–Sharma, KT Diff, and optionally E5 methods, used in only one iteration, with λ m d as an initial value for induced inflow ratio in solving Equation (2), have all attributes to become methods of choice for solving the GIF for a single-rotor helicopter in steady-level flight conditions below dynamic stall limit.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/aerospace10030238/s1, Figure S1: Working Group and Control Group results of No limit model; Figure S2: 𝐼𝑇𝑁𝑢𝑚𝑏 percentage distribution of Equation (6): (a) One-point methods Halley and Kiss; (b) One-point methods Chebishev and E5; (c) Two-point methods Ostrowski and Jarrat; (d) Two-point methods Chun and Cordero; (e) Three-point methods Sharma and BWR II; (f) Three-point methods KT Diff and BWR.; Figure S3: Average ITNumb distribution of Equation (5): (a) One-point methods; (b) Two-point methods; (c) Three-point methods; Figure S4: 𝐼𝑇𝑁𝑢𝑚𝑏 percentage distribution of Equation (5): (a) One-point methods; (b) Two-point methods; (c) Three-point methods; Figure S5: Average ITNumb distribution of Equation (3): (a) One-point methods; (b) Two-point methods; (c) Three-point methods; Figure S6:𝐼𝑇𝑁𝑢𝑚𝑏 percentage distribution of Equation (3): (a) One-point methods Halley and Kiss; (b) One-point methods Chebishev and E5; (c) Two-point methods Ostrowski and Jarrat; (d) Two-point methods Chun and Cordero; (e) Three-point methods.; Figure S7: Average ITNumb distribution of Equation (2): (a) One-point methods; (b) Two-point methods; (c) Three-point methods; Figure S8:𝐼𝑇𝑁𝑢𝑚𝑏 percentage distribution of Equation (2): (a) One-point methods; (b) Two-point methods; (c) Three-point methods; Figure S9: The maximum output errors Eo max vs advance ratio for Equation (6): (a) One-point methods; (b) Two-point methods; (c) Three-point methods; Figure S10: The maximum output errors Eo max vs advance ratio for Equation (5): (a) One-point methods; (b) Two-point methods; (c) Three-point methods; Figure S11: The maximum output errors Eo max vs advance ratio for Equation (3): (a) One-point methods; (b) Two-point methods; (c) Three-point methods; Figure S12: The maximum output errors Eo max vs advance ratio for Equation (2): (a) One-point methods; (b) Two-point methods; (c) Three-point methods; Figure S13: Comparative view of Et max for solutions of Equation (2) obtained in one iteration by three-point methods with 𝜆𝑚𝑑 initial value in the APR zone and out of the APR zone (0.15 ≤ 𝐶𝑇 ≤ 0.2; 20° ≤ ∝𝐷 ≤ 35°). Figure S14: Maximum number of iterations for 4 GIF forms with 3 initial values: Newton-Raphson method (NR) and Newton-Raphson method with relaxation factor (NRf); Table S1: Percentage of failed convergence, PoF (%) generated by Newton-Raphson method used with (NRf) and without (NR) relaxation factor in solving of Equation (6), ε = 5 Δ 10 n . Symbolic Equations (S1)–(S16) represent 16 variants of the LFE, and they are given in the HTML file: Code No.1, which is uploaded along with the article. The complete codes with obtained results developed and used for NCP_1 execution can be downloaded at: https://github.com/marjanvtsl/LFEM.git, accessed on 30 December 2022. The code for NCP_2 with computational cost tests is given at: https://github.com/marjanvtsl/GIF-Numerical-Analysis.git, accessed on 23 February 2023.

Author Contributions

All authors discussed and agreed upon the idea, and made scientific contributions: Conceptualization, M.D. (Marjan Dodic), B.K., M.D. (Mirko Dinulovic) and A.B.; Data curation, M.D. (Marjan Dodic) and B.K.; Formal analysis, M.D. (Marjan Dodic), B.K., B.R., M.D. (Mirko Dinulovic) and A.B.; Investigation, M.D. (Marjan Dodic) and B.K.; Methodology, M.D. (Marjan Dodic), B.K., B.R. and A.B.; Resources, M.D. (Marjan Dodic) and B.K.; Software, M.D. (Marjan Dodic), B.K. and M.D. (Mirko Dinulovic); Supervision, B.K.; Validation, M.D. (Marjan Dodic), B.K. and B.R.; Visualization, M.D. (Marjan Dodic) and B.K.; Writing—original draft, M.D. (Marjan Dodic) and B.K.; Writing—review and editing, M.D. (Marjan Dodic) and B.K. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Not applicable.

Acknowledgments

This research was supported by the University of Defence in Belgrade.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

a rotor tip speed statistical factor
a m i n minimum value of rotor tip speed statistical factor
a m a x maximum value of rotor tip speed statistical factor
a l blade lift-curve slope with compressibility effect
a l 0 blade lift-curve slope without compressibility effect
a s speed of sound
a 1 s first harmonic coefficient of longitudinal blade flapping with respect to shaft
b rotor diameter statistical factor
b m i n minimum value of rotor diameter statistical factor
b m a x maximum value of rotor diameter statistical factor
b 1 s first harmonic coefficient of lateral blade flapping with respect to shaft
c blade chord
c d main drag coefficient of rotor blades
c d T R main drag coefficient of tail rotor blades
f e equivalent helicopter drag area, including fuselage, gear, rotor hub and tail, D / 0.5 ρ V 2
f e c equivalent helicopter drag area at α F = 0
f e l equivalent helicopter lift area
h altitude
h m i n minimum altitude
h m a x maximum altitude
Δ h altitude increment
k equivalent helicopter drag area empirical factor
k m i n minimum value of equivalent helicopter drag area empirical factor
k m a x maximum value of equivalent helicopter drag area empirical factor
k i induced power empirical factor
Δ m mass increment
r f convergence relaxation factor
C l a v mean lift coefficient, 6 C t σ 1 1 + 3 2 μ 2
C p power coefficient, P R / ρ A Ω 3 R 3
Δ C p profile power increase due to the compressibility effect
C T thrust coefficient, T / ρ A Ω 2 R 2
C T m i n minimum thrust coefficient
C T m a x maximum thrust coefficient
D parasite drag of rotorcraft
D f fuselage parasite drag force
E o local relative error during iterative process
E o m a x maximum local relative error during iterative process
E t real relative error
E t m a x maximum real relative error
F p profile power factor
H rotor drag force or H-force
H T R tail rotor drag force
ITNumbnumber of iterations
L f fuselage lift force,
M Mach number, 0.75 V t / a s
M r blade advancing Mach number at 0.75 R , Vcos α D + 0.75 V t / a s
M a t blade advancing-tip Mach number, Vcos α D + V t / a s
M t i p blade advancing-tip Mach number at hover, V t / a s
P R required power for forward flight
P T O maximum (available) take-off power
R rotor radius
R c correlation coefficient
R e Reynolds number, 0.75 V t Δ c / ν
R e l o c blade retreating Reynolds number at 0.75 R
T main rotor thrust
v i rotor induced velocity
V , V rotor velocity with respect to the air or flight speed or true airspeed
V m i n minimum flight speed
V m a x maximum flight speed
Δ V flight speed increment
V f local air velocity with respect to the fuselage, V 2 + 2 v i V sin α D + v i 2
V t rotor tip speed, Ω R
W M T O maximum take-off mass
W T O take-off mass
α l airfoil angle of attack, C l a v a l + α n
α n airfoil zero-lift angle of attack
α D rotor disc plane incidence angle or angle of attack, positive down for forward tilt
α D m i n minimum rotor disc plane incidence angle
α D m a x maximum rotor disc plane incidence angle
α F fuselage angle of attack, positive down for forward tilt
α M mast tilt, positive forward
Δ α D downwash angle produced at fuselage by the main rotor
δ atmospheric (static) pressure ratio, p / p o , p o = 101325   P a
δ o blade average profile drag coefficient without compressibility effect
δ 0 blade average profile drag coefficient with compressibility effect
θ atmospheric (static) temperature ratio, T / T o , T o = 288.15   K
θ F P flight path angle, positive down
ε relative error, stopping criteria
ε a v average error [%]
ε M A X maximum error [%]
λ rotor inflow ratio, Vsin α D + v i / Ω R or λ i + μ z , positive down through disc
λ 0 initial value of rotor (induced) inflow ratio for numerical solution
λ e analytical solution of rotor inflow ratio λ for edgewise flight ( α D = 0 )
λ w j initial value of rotor inflow ratio proposed by Wayne Johnson
λ m d initial value of rotor inflow ratio proposed by Wayne Johnson—modified
λ i induced inflow ratio, v i / Ω R
λ h rotor inflow ratio in hover
λ f rotor inflow ratio or induced inflow ratio numerically obtained by standard precision
λ f 1 value from previous iteration
λ H P S high precision solution of rotor inflow ratio (benchmark solution)
μ rotor advance ratio, Vcos α D / Ω R
μ z rotor axial velocity ratio, Vsin α D / Ω R or μ tan α D
μ c rotor velocity (with respect to the air) ratio, V / Ω R or μ 2 + μ z 2
ρ air density
ρ I S A air density according to the International Standard Atmosphere
ν kinematic viscosity
σ rotor solidity, N c / π R
σ m i n minimum rotor solidity
σ m a x maximum rotor solidity
σ m e a n mean rotor solidity
Ω rotor rotational speed
Ω n nominal rotor rotational speed
APRacceptable parameter range
BERPBritish Experimental Rotor Programme
CPUcentral processing unit
CFDcomputational fluid dynamics
CSDcomputational structural dynamics
CGcontrol group
GIFGlauert inflow formula
GPUgraphics processing unit
ISAInternational Standard Atmosphere
MTOMmaximum take-off mass
MTOMminminimum value of maximum take-off mass
MTOMmaxmaximum value of maximum take-off mass
NaNnot a number
NCPnumerical computational plan
NRNewton–Raphson
NoFnumber of failed convergence
PoFpercentage of failed convergence
RWrotary wing
TOMtake-off mass
UAVunmanned aerial vehicle
WGworking group
1P2Rone-point two-order of convergence
1P3Rone-point three-order of convergence
1P4Rone-point four-order of convergence
1P4Rone-point five-order of convergence
2P4Rtwo-point five-order of convergence
3P8Rthree-point eight-order of convergence

Appendix A

Longitudinal force equilibrium model (LFEM)
The main two equations of the model are previously given in the study, which is Equations (14) and (15). Impossibility of using a small angle assumption leads to the more complex longitudinal equilibrium model, given below.
Airframe drag is usually written in terms of an equivalent parasite drag area as:
D = 0.5 ρ V f 2 f e ,
where f e   is equivalent parasite drag area, which depends on the fuselage angle of attack. For this study, the following model is adopted:
f e = f e c A + B sin n α F ,  
where A = 1 ,   B = 10 , and n = 3 . Constant f e c is a value of the equivalent parasite drag area at α F = 0 , usually represented as a function of maximum take-off weight as:
f e c = k W M T O / 1000 2 / 3
According to [24,34], the constant k is based on historical data and varies from k max = 1.3533 for high airframe drag to   k min = 0.3934 for the very best industrial solutions. The average value is k av = 0.66 . Fuselage lift force is calculated as
L f = 0.5 ρ V f 2 f e l .
Equivalent lift area f e l   is also scalable, based on the helicopter model from [25] (p. 679), as in Table A1.
Table A1. Equivalent lift area, f e l .
Table A1. Equivalent lift area, f e l .
α F f e l   1
α F 15 ° 0.113 f e c m α F 0.282 f e c m
15 ° < α F 22.5 ° 0.0133 f e c m α F 1.7 f e c m
22.5 ° < α F 30 ° 2 f e c m
α F > 30 ° 0
1  f e c m is f e c calculated for k av = 0.66 .
As is shown in Figure 2b, the fuselage angle of attack can be evaluated as
α F = 180 π ( α D α M + Δ α D ) ,
where Δ α D is a downwash angle caused by the main rotor stream and, according to Figure 2b, calculated as
Δ α D = a r c t a n λ i c o s α D μ c + λ i s i n α D ,  
where λ i as modeled with an assumption that C T C w , as
λ i = C w 2 λ e + μ tan α D 2 + μ 2 ,
The main rotor drag force H is roughly estimated as follows:
H = 3.65 8 π ρ σ c d V t V R 2 cos α D ,  
and the tail rotor drag force H T R as
H T R = 3.65 8 π ρ σ T R c d T R V t T R V R T R 2 cos α D ,  
where c d T R = 0.01 and σ T R = 0.1 . V t T R and R T R   are obtained by statistics, given in [23].
The most used expression for c d in low-fidelity rotor models is the quadratic equation for NACA 23012 airfoil at R e = 2 × 10 6 , developed by Bailey [35]:
c d = 0.0087 0.0216 α l + 0.645 α l 2 .
The biggest disadvantage of Bailey’s expression is a significant drop in accuracy for α l > 11.8 ° [21]. Therefore, based on experimental test data for 0 α l 20 ° for NPL9618 airfoil given in [36] and with the use of the GraphPad PRISM® regression analysis tool, a cubic expression is developed:
c d = δ 0 R e ,   M + 0.02274 α l 0.6453 α l 2 + 5.001 α l 3 ,
where δ 0 R e ,   M represents a blade profile drag, and is calculated as given in Table A2.
Table A2. Rotor blade profile drag estimation.
Table A2. Rotor blade profile drag estimation.
M r   10 6 R e 10 7   R e < 10 6 ,   R e > 10 7
M r 0.8 δ 0 = 0.015058 0.0010138 log R e δ 0 = 0.00867 2 × 10 6 R e 1 8
0.8 < M r 1 δ 0 = δ M r 0.8 + 0.4 ( M r 0.8 )
M r > 1 1
For R e = 2 × 10 6 and M r < 0.8 , the Equation (A11) becomes:
c d = 0.00867 + 0.02274 α l 0.6453 α l 2 + 5.001 α l 3
NPL9618 airfoil has a superior drag characteristic over NACA 23012 and is used for the Lynx helicopter main rotor blades. The modeling of Reynolds and Mach numbers effects on the blade average profile drag coefficient is based on [37,38] and is shown in Table A2. For the determination of an average blade angle of attack, well-known relations for α l , C l a v and C t have been used (see nomenclature).
Generally, an airfoil lift-curve slope a l primarily depends on the Mach number and significantly less on the Reynolds number [30,37,39,40]. Therefore, the Reynolds number influence is neglected. As shown in Table A3, the lift-curve slope modeling is based on a simple Prandtl–Glauert correction for M r 0.9 . For 0.9 < M r 1.1 , the McCroskey expression has been used (Equation (10) from [37], where t / c = 0.12 ). For higher M r , the lift-curve slope is practically equal to zero.
Table A3. Main rotor blade lift-curve slope determination.
Table A3. Main rotor blade lift-curve slope determination.
M r a l
M r 0.9 a l o 1 0.75 M r 2
0.9 < M r 1.1 3.1515 0.26 M r 2 1 3
The system of nonlinear equations within the longitudinal force equilibrium model (LFEM) is transformed into a single nonlinear algebraic longitudinal force equilibrium (LFE) equation with α D as a variable using symbolic algebra. Furthermore, depending on the values of M r and α F , there are 16 variants of LFE, represented by Equations (S1)–(S16). LFE is successively solved by the NR method about three times faster than the classic Newton-Kantorovich method for solving systems of nonlinear equations.

References

  1. Glauert, H. The Theory of the Autogyro. J. R. Aeronaut. Soc. 1927, 31, 483–508. [Google Scholar] [CrossRef]
  2. Peters, D.A. How Dynamic Inflow Survives in the Competitive World of Rotorcraft Aerodynamics. J. Am. Helicopter Soc. 2009, 54, 011001. [Google Scholar] [CrossRef] [Green Version]
  3. Smith, M.J.; Moushegian, A. Dual-Solver Hybrid Computational Approaches for Design and Analysis of Vertical Lift Vehicles. Aeronaut. J. 2022, 126, 187–208. [Google Scholar] [CrossRef]
  4. Chen, R.T. A survey of nonuniform inflow models for rotorcraft flight dynamics and control applications. In Proceedings of the Fifteenth European Rotorcraft Forum, Amsterdam, The Netherland, 12–15 September 1989. [Google Scholar]
  5. Peters, D.A.; Modarres, M.; Howard, A.B.; Rahming, B. A Third Approximation to Glauert’s Momentum Theory. J. Am. Helicopter Soc. 2016, 61, 042007. [Google Scholar] [CrossRef]
  6. Rubin, R.L.; Zhao, D. New Development of Classical Actuator Disk Model for Propellers at Incidence. AIAA J. 2021, 59, 1040–1054. [Google Scholar] [CrossRef]
  7. Petrovic, Z.; Stupar, S.; Kostic, I.; Simonovic, A. Determination of a Light Helicopter Flight Performance at the Preliminary Design Stage. Stroj. Vestn. J. Mech. Eng. 2010, 56, 535–543. [Google Scholar]
  8. Traub, J.F. Iterative Methods for the Solution of Equations, 2nd ed.; Chelsea Publishing Company: New York, NY, USA, 1982. [Google Scholar]
  9. Petkovic, M.S.; Neta, B.; Petkovic, L.; Dzunic, J. Multipoint Methods for Solving Nonlinear Equations, 1st ed.; Elsevier Science Publishing Co., Inc.: San Diego, CA, USA, 2014; ISBN 9780123972989. [Google Scholar]
  10. Dzunic, J. Multipoint Methods for Solving Nonlinear Equations. Ph.D. Thesis, University of Nis, Nis, Serbia, 2012. [Google Scholar]
  11. Petkovic, M.S.; Neta, B.; Petkovic, L.; Dzunic, J. Multipoint methods for solving nonlinear equations: A survey. Appl. Math. Comput. 2014, 226, 635–660. [Google Scholar] [CrossRef] [Green Version]
  12. Berahas, A.S.; Jahani, M.; Richtárik, P.; Takáč, M. Quasi-Newton methods for machine learning: Forget the past, just sample. Optim. Methods Softw. 2022, 37, 1668–1704. [Google Scholar] [CrossRef]
  13. Goldfarb, D.; Ren, Y.; Bahamou, A. Practical quasi-newton methods for training deep neural networks. Adv. Neural Inf. Process. Syst. 2020, 33, 2386–2396. [Google Scholar] [CrossRef]
  14. Davis, K.; Schulte, M.; Uekermann, B. Enhancing Quasi-Newton Acceleration for Fluid-Structure Interaction. Math. Comput. Appl. 2022, 27, 40. [Google Scholar] [CrossRef]
  15. Bogaers, A.E.J.; Kok, S.; Reddy, B.D.; Franz, T. Quasi-Newton methods for implicit black-box FSI coupling. Comput. Methods Appl. Mech. Eng. 2014, 279, 113–132. [Google Scholar] [CrossRef] [Green Version]
  16. Akbari, A.; Giannacopoulos, D. An efficient multi-threaded Newton–Raphson algorithm for strong coupling modeling of multi-physics problems. Comput. Phys. Commun. 2021, 258, 107563. [Google Scholar] [CrossRef]
  17. Liu, L.; Zhang, H.; Wu, Y.; Liu, B.; Guo, J.; Li, F. A modified JFNK with line search method for solving k-eigenvalue neutronics problems with thermal-hydraulics feedback. Nucl. Eng. Technol. 2023, 55, 310–323. [Google Scholar] [CrossRef]
  18. Ganguli, R.A. Pedagogical Example for STEM Using the Glauert Inflow Equation, Mathematica and Python. In Proceedings of the AIAA SciTech Forum, San Diego, CA, USA, 7–11 January 2019. [Google Scholar] [CrossRef]
  19. Bronshtein, I.N.; Semendyayev, K.A.; Musiol, G.; Muehlig, H. Handbook of Mathematics, 5th ed.; Springer: Berlin/Heidelberg, Germany, 2007; ISBN 978-3-540-72121-5. [Google Scholar]
  20. Leishman, G.J. Principles of Helicopter Aerodynamics, 2nd ed.; Cambridge University Press: Cambridge, UK, 2006. [Google Scholar]
  21. Johnson, W. Rotorcraft Aeromechanics; Cambridge University Press: New York, NY, USA, 2013. [Google Scholar]
  22. Johnson, W. NDARC—NASA Design and Analysis of Rotorcraft: Theory; Appendix 3, Release 1.11; TP–2015-218751; National Aeronautics and Space Administration: Washington, DC, USA, 2016.
  23. Rand, O.; Khromov, V. Helicopter Sizing by Statistics. J. Am. Helicopter Soc. 2004, 49, 300–317. [Google Scholar] [CrossRef]
  24. Harris, F.D. Introduction to Autogyros, Helicopters, and Other V/STOL Aircraft, Volume II: Helicopters; National Aeronautics and Space Administration: Washington, DC, USA, 2012.
  25. Prouty, R. Helicopter Performance, Stability and Control, 2nd ed.; Krieger Publishing Company: Malabar, FL, USA, 2002. [Google Scholar]
  26. Keys, C.N.; Wiesner, R. Guidelines for Reducing Helicopter Parasite Drag; Boeing Vertol Company: Ridley Park, PA, USA, 1974. [Google Scholar]
  27. Bramwell, A.R.S.; Done, G.; Balmford, D. Bramwell’s Helicopter Dynamics, 2nd ed.; Butterworth-Heinemann: Oxford, UK, 2000. [Google Scholar]
  28. Harris, F.D. Rotor Performance at High Advance Ratio: Theory versus Test; NASA/CR-2008-215370; National Aeronautics and Space Administration: Washington, DC, USA, 2008.
  29. Harris, F.D. Rotary Wing Aerodynamics—Historical Perspective and Important Issues. In Proceedings of the American Helicopter Society National Specialist’s Meeting on Aerodynamics and Aeroacoustics, Arlington, TX, USA, 25–27 February 1987. [Google Scholar]
  30. McCroskey, W.J. A Critical Assessment of Wind Tunnel Results for the NACA 0012 Airfoil; NASA-TM-100019; National Aeronautics and Space Administration: Washington, DC, USA, 1987.
  31. Bezanson, J.; Edelman, A.; Karpinski, S.; Shah, V.B. Julia: A Fresh Approach to Numerical Computing. SIAM Rev. 2017, 59, 65–98. [Google Scholar] [CrossRef] [Green Version]
  32. Sharma, J.R.; Sharma, R. A new family of modified Ostrowski’s methods with accelerated eighth order convergence. Numer. Algor. 2010, 54, 445–458. [Google Scholar] [CrossRef]
  33. Lau, B.H.; Louie, A.W.; Griffiths, N.; Sotiriou, C.P. Performance and Rotor Loads Measurements of the Lynx XZ170 Helicopter with Rectangular Blades; NASA-TM-104000; National Aeronautics and Space Administration: Washington, DC, USA, 1993.
  34. Dodic, M.; Krstic, B.B.; Rasuo, B.P. Single Rotor Helicopter Parasite Drag Estimation in the Preliminary Design Stage. Tehnika 2018, 73, 65–77. [Google Scholar] [CrossRef] [Green Version]
  35. Bailey, F.J., Jr. A Simplified Theoretical Method of Determining the Characteristics of a Lifting Rotor in Forward Flight; NACA Report 716; US Government Printing Office: Washington, DC, USA, 1941.
  36. Lau, B.H.; Louie, A.W.; Griffiths, N.; Sotiriou, C.P. Correlation of the Lynx-XZ170 Flight-Test Results Up to and Beyond the Stall Boundary. In Proceedings of the AHS Forum 49, St. Louis, MO, USA, 19–21 May 1993. [Google Scholar]
  37. Totah, J. A Critical Assessment of UH-60 Main Rotor Blade Airfoil Data; NASA-TM-103985; National Aeronautics and Space Administration: Washington, DC, USA, 1993.
  38. Yamauchi, G.K.; Johnson, W. Trends of Reynolds Number Effects on Two-Dimensional Airfoil Characteristics for Helicopter Rotor Analyses; NASA-TM-84363; National Aeronautics and Space Administration: Washington, DC, USA, 1983.
  39. Rasuo, B. The influence of Reynolds and Mach numbers on two-dimensional wind-tunnel testing: An experience. Aeronaut. J. 2011, 115, 1166. [Google Scholar] [CrossRef]
  40. Rasuo, B. Scaling between Wind Tunnels–Results Accuracy in Two-Dimensional Testing. Trans. Jpn. Soc. Aeronaut. Space Sci. 2012, 55, 109–115. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Rotor disc in forward flight.
Figure 1. Rotor disc in forward flight.
Aerospace 10 00238 g001
Figure 2. (a) A simplified view of determining the GIF parameters within the longitudinal force equilibrium model (LFEM); (b) Single-rotor helicopter longitudinal force model in forward flight.
Figure 2. (a) A simplified view of determining the GIF parameters within the longitudinal force equilibrium model (LFEM); (b) Single-rotor helicopter longitudinal force model in forward flight.
Aerospace 10 00238 g002
Figure 3. Graphical representation of the APR boundaries determination: NCP 1 results (left), Power 1 and Dynamic stall limits (middle), APR with safety margin (right).
Figure 3. Graphical representation of the APR boundaries determination: NCP 1 results (left), Power 1 and Dynamic stall limits (middle), APR with safety margin (right).
Aerospace 10 00238 g003
Figure 4. (a) Modified Harris criterion for Dynamic stall limit obtained by WG; (b) Dynamic stall boundaries in accordance to Lynx XZ170 flight test data.
Figure 4. (a) Modified Harris criterion for Dynamic stall limit obtained by WG; (b) Dynamic stall boundaries in accordance to Lynx XZ170 flight test data.
Aerospace 10 00238 g004
Figure 5. PoF distribution over advance ratio, Equation (6), ε = 5 × 10 16 (a) One-point methods; (b) Two-point methods; (c) Three-point methods.
Figure 5. PoF distribution over advance ratio, Equation (6), ε = 5 × 10 16 (a) One-point methods; (b) Two-point methods; (c) Three-point methods.
Aerospace 10 00238 g005
Figure 6. ITNumb distribution over advance ratio, Equation (6), λ 0 = λ h ,   ε = 5 × 10 14 , one-point methods: (a) Halley and Kiss methods; (b) Chebyshev and E5 methods.
Figure 6. ITNumb distribution over advance ratio, Equation (6), λ 0 = λ h ,   ε = 5 × 10 14 , one-point methods: (a) Halley and Kiss methods; (b) Chebyshev and E5 methods.
Aerospace 10 00238 g006
Figure 7. Average I T N u m b distribution of Equation (6): (a) One-point methods; (b) Two-point methods; (c) Three-point methods.
Figure 7. Average I T N u m b distribution of Equation (6): (a) One-point methods; (b) Two-point methods; (c) Three-point methods.
Aerospace 10 00238 g007
Figure 8. E o   M a x (maximum relative error between the last two iterations or points) for all GIF forms solved by: (a) one-point methods, (b) two-point methods; (c) three-point methods.
Figure 8. E o   M a x (maximum relative error between the last two iterations or points) for all GIF forms solved by: (a) one-point methods, (b) two-point methods; (c) three-point methods.
Aerospace 10 00238 g008
Figure 9. Maximum output relative error of the benchmark solution.
Figure 9. Maximum output relative error of the benchmark solution.
Aerospace 10 00238 g009
Figure 10. E t   M a x (maximum relative error related to the benchmark solution) for all GIF forms solved by: (a) one-point methods, (b) two-point methods; (c) three-point methods.
Figure 10. E t   M a x (maximum relative error related to the benchmark solution) for all GIF forms solved by: (a) one-point methods, (b) two-point methods; (c) three-point methods.
Aerospace 10 00238 g010
Figure 11. E t   M a x for all methods and initial values—Equation (6).
Figure 11. E t   M a x for all methods and initial values—Equation (6).
Aerospace 10 00238 g011
Figure 12. Percentage of convergence towards superfluous roots for all examined methods.
Figure 12. Percentage of convergence towards superfluous roots for all examined methods.
Aerospace 10 00238 g012
Figure 13. E t   M a x (maximum relative error related to the benchmark solution) for rational GIF forms solved by: (a) one-point methods, (b) two-point methods; (c) three-point methods.
Figure 13. E t   M a x (maximum relative error related to the benchmark solution) for rational GIF forms solved by: (a) one-point methods, (b) two-point methods; (c) three-point methods.
Aerospace 10 00238 g013
Figure 14. E t   M a x (maximum relative error related to the benchmark solution) for irrational GIF forms solved by: (a) one-point methods, (b) two-point methods; (c) three-point methods.
Figure 14. E t   M a x (maximum relative error related to the benchmark solution) for irrational GIF forms solved by: (a) one-point methods, (b) two-point methods; (c) three-point methods.
Aerospace 10 00238 g014
Figure 15. The effects of an arbitrary precision on maximum relative error reduction for three-point methods used in solving of Equations (3), (5) and (6) with an initial value λ m d .
Figure 15. The effects of an arbitrary precision on maximum relative error reduction for three-point methods used in solving of Equations (3), (5) and (6) with an initial value λ m d .
Aerospace 10 00238 g015
Figure 16. (a) The example of discrepancy between output relative error and actual relative error after one iteration (three-point methods, λ o = λ m d , Equation (2)); (b) Three-point Sharma–Sharma numerical method—the error accumulation in the third step of the first iteration.
Figure 16. (a) The example of discrepancy between output relative error and actual relative error after one iteration (three-point methods, λ o = λ m d , Equation (2)); (b) Three-point Sharma–Sharma numerical method—the error accumulation in the third step of the first iteration.
Aerospace 10 00238 g016
Figure 17. Average number of iterations for 4 GIF forms solved by Newton–Raphson method (NR) and Newton–Raphson method with relaxation factor (NRf).
Figure 17. Average number of iterations for 4 GIF forms solved by Newton–Raphson method (NR) and Newton–Raphson method with relaxation factor (NRf).
Aerospace 10 00238 g017
Figure 18. Maximum relative error related to the benchmark solution ( E t M a x ) for NR and NRf methods in solving: (a) Equations (5) and (6); (b) Equations (2) and (3).
Figure 18. Maximum relative error related to the benchmark solution ( E t M a x ) for NR and NRf methods in solving: (a) Equations (5) and (6); (b) Equations (2) and (3).
Aerospace 10 00238 g018
Figure 19. Maximum relative error related to the benchmark solution ( E t M a x ) for NR and NRf methods obtained in one iteration.
Figure 19. Maximum relative error related to the benchmark solution ( E t M a x ) for NR and NRf methods obtained in one iteration.
Aerospace 10 00238 g019
Figure 20. Maximum relative error ( E t m a x ) related to the benchmark solution for analytical (closed form) solution and the numerical solutions obtained in one iteration by a high order of convergence methods with an initial value λ 0 = λ m d .
Figure 20. Maximum relative error ( E t m a x ) related to the benchmark solution for analytical (closed form) solution and the numerical solutions obtained in one iteration by a high order of convergence methods with an initial value λ 0 = λ m d .
Aerospace 10 00238 g020
Table 1. Coefficients of the GIF rational form.
Table 1. Coefficients of the GIF rational form.
Constant Variable   x = λ i Variable   x = λ
A 4 4 4
A 3 8 μ tan α D -   8 μ   tan α D
A 2 4 μ 2 tan 2 α D + 1 4 μ 2 ( tan 2 α D + 1 )
A 1 0   8 μ 3 tan   α D
A 0 -   C T 2 4 μ 4 tan 2 α D C T 2
Table 2. Initial guesses λ 0 for GIF numerical solutions.
Table 2. Initial guesses λ 0 for GIF numerical solutions.
Main Inflow Ratio λ
Equations (5) and (6)
Induced Inflow Ratio λ i ,
Equations (2) and (3)
Reference
λ h λ h Leishman [20]
λ h 2 λ h + μ z 2 + μ 2 + μ z λ h 2 λ h + μ z 2 + μ 2 W. Johnson [15,16]
λ h 2 λ e + μ z 2 + μ 2 + μ z λ h 2 λ e + μ z 2 + μ 2 Current study
Table 3. Finite number of cases of Working (WG) and Control Groups (CG).
Table 3. Finite number of cases of Working (WG) and Control Groups (CG).
CharacteristicsGroup 1Group 2Group 3Group 4
MTOMmin/Δm/MTOMmax 1WG1000/500/85009000/2000/19,00020,000/5000/35,00040,000/10,000/80,000
CG1000/200/88009000/400/19,60020,000/800/40,00041,600/1600/80,000
hmin/Δh/hmax 0/500/11,500
fe,max/fe,mean/fe,minCalculated based on [21,22,23,24,25,26]
Vmin/ΔV/Vmax1/1/131
σminmeanmax0.04/0.01/0.2
amin/amax130/150
bmin/bmax0.85/1.15
kmin/kmax0.39/1.35
No. of WG points3,621,8881,358,208905,4721,131,840
No. of CG points32,823,36022,155,76821,335,18420,514,600
Total No. of WG points:7,017,408Total No. of CG points:96,828,912
1 Three different TOM are used for each MTOM: empty, mean, and max (see the Julia® LFEM code).
Table 4. The rotor induced power empirical factor k i .
Table 4. The rotor induced power empirical factor k i .
μ k i   Reference
0 < μ 1 1.15 ÷ 1.2 Vast literature
μ 0.5 1.075 cosh 6.76 μ 2 [28]
0.5 < μ 1 1.0 29.332 μ + 92.439 μ 2 51.746 μ 3
μ 0.2 1.2Current study
0.2 < μ 0.4 2 μ + 0.8
μ > 0.4 14 μ 4
Table 5. Numerical Computational Plan 1 (NCP 1).
Table 5. Numerical Computational Plan 1 (NCP 1).
Cycle No.GroupModelCycle No.GroupModel
1.WG 1–4No limit5.WG 1–4Power 2 limit
2.WG 1–4Power limit 6.WG 1–4Power 3 limit
3.WG 1–4Dynamic stall7.CG 1–4No limit
4.WG 1–4Power 1 limit
Table 6. GIF Numerical Computing Plan (NCP 2).
Table 6. GIF Numerical Computing Plan (NCP 2).
MethodsTypeStopping
Criteria
Equations Initial   Guesses   λ o
Halley1P3R ε < 5 × 10 4 Equation (6) λ h
Chebyshev1P3R ε < 5 × 10 6 Equation (5) λ w j
Kiss1P4R ε < 5 × 10 8 Equation (3) λ m d
E51P5R   ε < 5 × 10 10 Equation (2)
NR1P2R   ε < 5 × 10 12
NR with f = 0.5-   ε < 5 × 10 14
Ostrowski2P4R   ε < 5 × 10 16
Chun2P4R1 iteration
Cordero2P4R
Jarratt2P4R
Sharma–Sharma3P8R
BWR 3P8R
BWR II3P8R
KT 3P8R
Table 7. Percentage of Failed Convergence, PoF (%), ε = 5 × 10 n , Equation (6).
Table 7. Percentage of Failed Convergence, PoF (%), ε = 5 × 10 n , Equation (6).
λ 0 nOne-Point MethodsTwo-Point MethodsThree-Point Methods
HalleyChebyshevKissE5OstrowskiChunCorderoJarratSharmaBWRBWR IIKT Diff
λ h 162.162.152.152.150.7612.3912.842.151.622.681.281.75
1400001 × 10−46 × 10−36 × 10−300.01500.0370.04
12000000002 × 10−400.0240.002
1000000000000.0170
800000000001.2 × 10−40
600000000005.9 × 10−50
400000000009.8 × 10−60
λ w j 162.122.082.102.110.6313.6914.132.11.682.711.202.37
1400005 × 10−55 × 10−35 × 10−300.01200.050.032
12000000005 × 10−500.0240.002
1000000000000.072 × 10−5
800000000000.00150
600000000007.3 × 10−60
400000000003 × 10−60
λ m 162.112.112.112.110.6413.6013.952.11.772.691.212.55
1400005 × 10−54 × 10−33 × 10−300.00800.040.033
12000000007 × 10−500.0320.001
1000000000000.0272 × 10−5
800000000000.02230
600000000001.3 × 10−40
4000000000010−50
Table 8. Percentage of Failed Convergence, PoF (%), ε = 5 × 10 n , Equation (3).
Table 8. Percentage of Failed Convergence, PoF (%), ε = 5 × 10 n , Equation (3).
λ 0 nTwo-Point MethodsThree-Point Methods
Ostr.ChunCord.Jarr.Sh.Sh.BWRBWR IIKT
λ h 160.0030.0950.6280000.00150.0015
1400000000
λ w j 160.0060.1180.7710000.00210.001
1400000000
λ m d 160.0060.0870.5710000.00280.0015
1400000000
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Dodic, M.; Krstic, B.; Rasuo, B.; Dinulovic, M.; Bengin, A. Numerical Analysis of Glauert Inflow Formula for Single-Rotor Helicopter in Steady-Level Flight below Stall-Flutter Limit. Aerospace 2023, 10, 238. https://doi.org/10.3390/aerospace10030238

AMA Style

Dodic M, Krstic B, Rasuo B, Dinulovic M, Bengin A. Numerical Analysis of Glauert Inflow Formula for Single-Rotor Helicopter in Steady-Level Flight below Stall-Flutter Limit. Aerospace. 2023; 10(3):238. https://doi.org/10.3390/aerospace10030238

Chicago/Turabian Style

Dodic, Marjan, Branimir Krstic, Bosko Rasuo, Mirko Dinulovic, and Aleksandar Bengin. 2023. "Numerical Analysis of Glauert Inflow Formula for Single-Rotor Helicopter in Steady-Level Flight below Stall-Flutter Limit" Aerospace 10, no. 3: 238. https://doi.org/10.3390/aerospace10030238

APA Style

Dodic, M., Krstic, B., Rasuo, B., Dinulovic, M., & Bengin, A. (2023). Numerical Analysis of Glauert Inflow Formula for Single-Rotor Helicopter in Steady-Level Flight below Stall-Flutter Limit. Aerospace, 10(3), 238. https://doi.org/10.3390/aerospace10030238

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop