Review of vortex methods for rotor aerodynamics and wake dynamics

Electric vertical take-off and landing (eVTOL) aircraft with multiple lifting rotors or prop-rotors have received significant attention in recent years due to their great potential for next-generation urban air mobility (UAM). Numerical models have been developed and validated as predictive tools to analyze rotor aerodynamics and wake dynamics. Among various numerical approaches, the vortex method is one of the most suitable because it can provide accurate solutions with an affordable computational cost and can represent vorticity fields downstream without numerical dissipation error. This paper presents a brief review of the progress of vortex methods, along with their principles, advantages, and shortcomings. Applications of the vortex methods for modeling the rotor aerodynamics and wake dynamics are also described. However, the vortex methods suffer from the problem that it cannot deal with the nonlinear aerodynamic characteristics associated with the viscous effects and the flow behaviors in the post-stall regime. To overcome the intrinsic drawbacks of the vortex methods, recent progress in a numerical method proposed by the authors is introduced, and model validation against experimental data is discussed in detail. The validation works show that nonlinear vortex lattice method (NVLM) coupled with vortex particle method (VPM) can predict the unsteady aerodynamic forces and complex evolution of the rotor wake.

Rotorcraft undergo complex aerodynamic phenomena due to sharp gradients of velocity and pressure near the blade tips, strong wake vortices, compressible dynamic stall, large fluctuation amplitude, and unsteady flow reversal. Rotorcraft's complex aerodynamics create many challenges for engineers. Progress in the aviation industry has come from initial experimental rotorcraft development, conducted in wind tunnel facilities, along with simple mathematical calculations [1][2][3][4]. However, the limitation of experimental studies is that they can consider only a small number of aerodynamic rotorcraft layouts. Wind-tunnel tests also cannot be used to carry out multi-objective, multiparameter aerodynamic optimizations that address interactional effects. Nonetheless, wind-tunnel testing is required as the last step in the manufacturing process to validate numerical simulations [5][6][7][8]. A valuable list of empirical research done to study rotor flow fields is provided in the references [9][10][11].
With improvements in computing power and resources, engineers are now able to use numerical methods to analyze complex rotorcraft configurations and improve their understanding of areas of rotor aerodynamics that had been difficult to investigate via wind-tunnel tests [12][13][14]. A variety of numerical methods, ranging from engineeringlevel models to intensive high-fidelity approaches, have been developed and validated as predictive tools to analyze rotor aerodynamics and wake dynamics. The blade element momentum (BEM) theory was derived by combining the basic principles of the onedimensional momentum theory and blade element theory. The BEM theory has a simple structure and can rapidly calculate the rotor performance with very low computational cost. Therefore, the BEM theory is an engineering-level model, and has been widely used for design purposes. However, owing to its inherent assumptions, it can only describe quasi-steady solutions and ignores mutual aerodynamic interactions between adjacent blade elements. Computational Fluid Dynamics (CFD) methods were used to simulate rotor flow fields and had an undeniable impact on rotorcraft design developments. CFD methods are the high-fidelity, expensive approaches used to predict unsteady and transient phenomena by directly solving entire flow-fields containing rotor blades and downstream regions. As computing performance and resources have begun to grow dramatically, industrial and academic research has been attempting to simulate the isolated rotor, multirotor, and complete rotorcraft configurations using various CFD methods. In general, they can be classified into the Reynolds-Averaged Navier-Stokes (RANS), Large Eddy Simulation (LES), and Detached Eddy Simulation (DES) techniques depending on turbulence models with the range of length and time scales. However, CFD methods inherently suffer from grid-induced dissipation errors because of the numerical discretization over the flow field. Excessive dissipation error caused by the numerical discretization may result in rapid decay of the intensity of the rotor tip vortex downstream and a breakdown of the helical rotor wake structure. Moreover, although many modifications have been made to make CFD methods more efficient [15][16][17][18], high computational costs and complicated equation setup still make these methods unfeasible for designers aiming for rapid feedback to optimize their models [19,20].
Accordingly, there has been increasing demand for a fast and reliable method to analyze new eVTOL aircraft. Vortex methods are becoming useful tools to predict the unsteady flow physics and to study rotorcraft's complex wake flow. Although vortex methods are basically used for irrotational and incompressible flow, they can produce acceptable and accurate results for complex configuration designs in less time than the CFD methods [21][22][23][24]. For these reasons, most current comprehensive rotorcraft analysis codes have adopted vortex methods, in which the lifting surface (e.g., rotor blade) can be modeled in various ways, including the lifting line theory (LLT), Weissinger's LLT (extended lifting-line theory), the lifting surface method, and the source-doublet panel method. LLT or Weissinger's LLT use the simplest representation of the blade model in terms of vortex singularities, called horseshoe vortex filaments, along the spanwise direction. However, these models do not capture the three-dimensional effects on the rotor blade because the lifting surface is represented by only a single chordwise vortex element. The lifting surface method, referred to as the vortex lattice method (VLM), represents the curvature of the rotor blade using both the chordwise and the spanwise distribution of vortex ring elements on the actual camber surface. Therefore, the VLM has been shown to give much better representation of the three-dimensionality of the flow on the blade. The source-doublet panel method, in which the airfoil involved in the rotor blade is divided into upper and lower panels, can consider the thickness effect of the rotor blade. This is a major advantage compared to VLM. The rotor wake model is also used to represent the wake structure and to describe vorticity fields in the wake through the use of straight/curved vortex filaments or vortex particles. A vortex filament is the concentrated vortices along a segment with a singularity at the center, whereas a vortex particle is the concentrated vortices within a certain volume. These vortices are generated from the trailing-edge of each blade and convected downstream. Tables 1 and 2 present a simple overview of the blade and wake models, respectively, for the rotor simulation reported in this paper.
The most challenging problems of rotorcraft can be divided into three fields of engineering: aerodynamics, aeroacoustics, and aeroelastics. Researchers have used different vortex methods to tackle these challenges. The main topics in rotor aerodynamics are the influences of interactions among components [41,[45][46][47] and with nearby infrastructure (e.g., brownout and shipboard operation) [48,49], on rotorcraft performance, stability, control, and safety. There have also been several investigations on the design of rotorcraft configuration to improve aerodynamic performance [50][51][52]. The next Hess [33], Crispin [34] Wachspress et al. [35] Hybrid method Horseshoe vortex Only volume grid for flow field (CFD domain) Slow (Faster than only CFD simulation) Rajmohan and He [36] Zhao et al. [37] Bae and He [38]  rotorcraft design challenge is to reduce their noise, especially in urban environments. Blade-vortex interaction (BVI) is identified as the dominant source of noise in rotorcraft. Vortex methods are powerful tools for researchers to study interactions of blade wake vortices with each other as well as with other rotorcraft components [23,[53][54][55]. Rotorcraft are strongly affected by unsteady aerodynamic loads (e.g., BVI), which significantly contribute to vibration and structural deformation. Aeroelastic analysis has always been a crucial step in designing rotorcraft. With the new generation of compound rotorcraft, which generally have hingeless and bearingless rotors, rotor aeroelastic and response problems are attracting great attention in the aviation industry. Vortex methods are reliable and low-cost aerodynamic approaches for studying aeroelasticity in rotary wing aircraft [56][57][58].
The main objectives of this paper are to provide an overall review of vortex methods including principles, advantages, and drawbacks, as well as development and applications (specifically in the rotorcraft industry). In addition, the authors' recent progress in their numerical method of modeling rotor aerodynamics and wake dynamics is introduced, and its applications are discussed. The paper is organized as follows: Section 2 gives an overview of the principles of vortex methods and their fundamental formulation. Section 3 explains the three most common vortex methods used to predict rotor aerodynamics and introduces hybrid approaches. Applications of the vortex method for modeling rotor wake dynamics are presented in Section 4, which reviews the timemarching free-wake method and viscous vortex particle method in detail. Section 5 introduces our proposed nonlinear vortex lattice methods and compares our numerical results with experimental data. Finally, concluding remarks are made in Section 6.

Vortex theory
The vortex method is based on the assumption of inviscid, incompressible, and irrotational flow over the entire domain surrounding the body surface and wake region. The velocity field is obtained by solving Laplace's equation using appropriate velocity boundary conditions on the body surface and far-field. The continuity equation for incompressible potential flow, also represented by Laplace's equation for the velocity potential, is given by Eq. (1). Here, φ is the velocity potential; the velocity of an irrotational flow can be defined as in Eq. (2): Leishman and Bhagwat [39] Lee and Na [40] Vortex particle method Particles No grid Jang et al. [41], He and Zhao [42] Winckelmans and Leonard [43] Singh and Friedmann [44] Green's second identity can be used to find the general solution to Laplace's equation. The wake surface is supposed to be thin; it can be represented as the sum of only the doublet distribution on the wake surface S w , and the sum of source (σ) and doublet (μ) distribution on the body surfaces S B . The integral form used to generate the solution is derived as in Eq. (3) and Eq. (4), in which surface integration should be made over all boundaries containing singularity elements. Here, G is Green's function; β is a solid angle, of which the value is 4π when the point x is located at the outside of the body surface boundary for three-dimensional flow; n is the outward unit normal vector of the surface. σ and μ are the strength of the source and doublet singular elements, respectively.
The difference between the values of external and internal potential on the solid boundary defines the strength of the doublet elements in Eq. (5); discontinuity in the normal derivative of the velocity potential on the solid boundary can be referred to as the strength of the source element in Eq. (6).
Using the definitions of the source and doublet strength, the general solution is rearranged as in Eq. (7). Equation (8) shows the resulting velocity induced by the source and doublet distribution on the body surface and wake.

Lifting line method
The lifting line method dates back to the development of Prandtl's lifting-line theory (LLT) [59]. It was the first analytical method used to predict lift and induced drag on lifting surfaces. In this theory, each spanwise section of a finite wing has a section lift equivalent to that acting on a similar section of an infinite wing having the same section circulation. The local circulation is related to the local aerodynamic force of a three-dimensional wing using the two-dimensional vortex lifting law of Kutta and Joukowski. LLT assumes an incompressible and inviscid fluid for which compressibility and viscous effects are negligible for application areas of interest. However, for Mach numbers less than 0.6, the effect of lowspeed compressibility can be introduced by the Prandtl-Glauert rule [60]. The classical theory assumes an infinitesimally thin and uncambered rigid flat plate with zero spanwise twist and zero sweep. Wing twist about a spanwise axis can be included as an alteration to the wing geometric angle of attack. Lifting surfaces, such as wings or rotor blades, are modeled as bound vortices with the strength of Ŵ(r) at the aerodynamic center. The goal of LLT is to determine Ŵ(r) as a function of the wing geometric properties; then, the Kutta-Joukowski theorem is used to obtain the lift per unit span [25,61]: where ρ is air density, V(r) is sectional flow speed along wing length, r is wing radial position measured from the wing root, and c is the chord. C lα,2d is the airfoil lift curve slope, α g and α i are the wing geometric angle and the induced angle of attack, respectively. Thus, Γ(r) can be computed as shown in Eq. (10), where w(r) is the induced downwash velocity distribution along the wing length.
LLT has wide application in evaluating the aerodynamics of wings with prescribed rotary and flapping motions. Goldstein [62] applied the original LLT of Prandtl for propellers. Lerbs [63] extended the formulations for moderately loaded propellers of arbitrary circulation distribution using the induction concept proposed by Kawada [64]. Conlisk [65] discussed the use of LLT for rotary wings in hover, emphasizing the importance of accounting for the influence of linear velocity variation along the blade on the bound circulation distribution. Johnson [26] emphasized the significance of applying modifications to LLT to handle specific rotary-wing aerodynamic phenomena such as wake periodicity, whereas Leishman [25] gave a generic formulation of LLT for rotary-wing motions.
LLT can be applied to a rotating frame with constant angular velocity [66][67][68]. The rotor encounters a uniform velocity field ( V ∞ ) aligned perpendicular to the rotation axis, rotating with constant angular velocity (Ω). The relative velocity field is calculated as in Eq. (11).
The Kutta-Joukowski theorem connects the lift force per unit length of the lifting line (L) to the total velocity field (V), as in Eq. (12). Here, ρ is the fluid density. The flow is assumed to be nearly two-dimensional at each radial position along the lifting lines.
The lifting line approach accurately predicts the hover performance of a wide range of conventional and advanced rotor designs. Landgrebe et al. [27,28] showed that the lifting line approach was adequate for predicting the aerodynamic performance of rotor blades in hover and forward flight conditions. Miller [69] explored the aerodynamics and dynamics of a rotor and the dynamics and control characteristics of a vehicle by modeling the rotor with LLT. Analytical modeling and design of the Apache helicopter were performed by Jones and Kunz [70,71] using CAMRAD II, which uses LLT for blade modeling. Yeo et al. [72,73] studied tail rotor flutter, exploring a wide range of design parameters and examining their effects on whirl flutter speed. Jain et al. [74] studied rotor performance in hover and forward flight and compared their results with experimental data. Interference of coaxial rotors in hover was studied by Ho et al. [75]. Wachspress et al. [76] also studied coaxial rotor performance with wake geometry and aeroacoustics by changing key rotor design parameters. Conceptual design of a slowedrotor compound helicopter emphasizing aerodynamic efficiency was performed by Moodie et al. [77].
Although classical LLT does not consider compressibility or viscosity and has limitations in its consideration of low aspect ratio wings and the thickness effect, it is a powerful tool to model lifting devices, offering a simple approach with low computational cost for the preliminary design of wings and for aerodynamic predictions. Though the method has some inherent limitations, researchers have proposed variations and adaptations over the years. Guermond [78] and Phillips [79] modified the classical LLT to account for curvature of swept wings. Several researchers also modified LLT to include unsteady effects [80][81][82].

Vortex lattice method
The vortex lattice method (VLM) is a branch of computational fluid dynamics that mathematically stems from finite-difference concepts [29]. In the VLM, lifting surfaces are approximated as thin surfaces represented by vortex sheets of unknown circulation. Because it allows the rotor blade to be discretized into both chordwise and spanwise directions, VLM can geometrically represent the curvature of a rotor blade surface, camber, and swept shape, a significant advantage of VLM compared to LLT [83].
Vortex lattices (or vortex rings) with circulation (Γ) are placed on the mean surface of a wing (neglecting the thickness of the wing) and the wake. The normal vector (n) is defined at the center of the vortex ring element, and is referred to as the collocation point (see Fig. 1). The leading edge of the vortex ring element and collocation point are located at the quarter and the three-quarter chord lines of each vortex panel, respectively. Because there is no flow separation in potential flow, the velocity field around the surface will be tangential. Moreover, no flow should cross the solid surface, in what is usually known as the no-penetration boundary condition, shown in Eqs. (13) and (14). The zero normal flow condition indicates that the sum of normal velocity components at each collocation point on the camber line induced by the rotor blade, wake, and timedependent kinematic velocity due to unsteady motion of the body must be zero.
Here, the V term is the time-dependent kinematic velocity, which is the sum of the system's velocity (V ∞ ), the relative velocity of the body (V rel ), and the body's frame (13) rotational velocity (Ω × r) at each collocation point. V ind,bound, and V ind,wake indicate the velocity components induced by the rotor blade and wake vortices.
Applying the no penetration boundary condition to the vortex sheet reduces the problem of finding the vortex strengths to a system of linear algebraic equations, as in Eq. (15). The bound-vortex circulations (Γ i ), and thus the self-induced velocity by the rotor blade can be represented by a combination of influence coefficients (a ij ). These influence coefficients are defined as the velocity component normal to the surface induced by the j-th vortex ring element with unit strength to the collocation point of the i-th vortex ring element. The right-hand side of the equation includes the contributions of wake-induced velocity and body kinematic velocities. Katz and Plotkin's book [60] provides a detailed description of the derivation of VLM. The VLM mathematical model has seen much improvement and its application has been extended to compressible fluids [84], unsteady flow [30], and many other complex aerodynamic phenomena [29,85]. Because the simplicity and low computational cost of VLMs increase their applicability, they are still being widely used in many aviation projects instead of CFD methods, which are far more demanding [31]. VLMs have helped researchers investigate the influences of rotorcraft configuration and flying environment on its aerodynamic performance [52], noise generation [55,86], flight control [87], and aeroelasticity response [88].

Source-doublet panel method
Although VLM is a useful tool for analyzing the aerodynamics of many rotorcraft configurations, it is not suitable for simulating flow-field over structures like fuselage or empennage. Engineers usually use the source-doublet panel method (also known as the panel method) to evaluate the aerodynamic performance of complex planforms [41,50,[89][90][91].
Like VLM, the source-doublet panel methods assume irrotational and incompressible flow. Therefore, they can also be used to solve Laplace's equations for potential function in the vortex field. The general solution of Laplace's equation for each element in an inertial system can be derived following Green's second identity. Considering the source (σ) and doublet (μ) distribution on the lifting surface and its wakes (see Fig. 2), the general solution of Eq. (1) can be expressed as: The integral equation can be made into a more straightforward form by applying the Dirichlet boundary condition to each of the collocation points in Eq. (19), where the total inner potential can be set equal to the freestream potential. To obtain a numerical solution of the integral equation, the body and wake surfaces are divided into a number of rectilinear panels consisting of N body surface panels and N W additional wake panels with constant-strength singularities. For the constant-strength source (σ) and doublet (μ) elements, the influences of body panel k and wake panel j at point P can be computed using Eqs. (21,22,23).

Fig. 2 Surface and wake panels arrangements on wing
The integral equation for finding the unknown doublet distribution (μ) on the body surface can be established as a linear algebraic equation, as in Eq. (24). Here, μ W is the strength of wake doublets, which can be expressed in terms of the unknown surface doublets (μ k ) by enforcing the Kutta condition as in Eq. (25). The Kutta condition implies that the wake doublets are related to the difference between the doublet strengths of the upper and lower panels of the trailing edge (See Fig. 2). In addition, for unsteady rotor simulation, the source strength (σ) can be determined by applying the zero normal flow condition on the surface, defined as in Eq. (26) with the local kinematic velocity.
Note that the strengths of the source and wake doublet are known at each time step. Moving the source term of Eq. (24) to the right-hand side of the equation results in a linear system of N equations with the unknown μ k , in the following form [60]: Once the linear system of equations is solved, the local velocity components of the surface panel coordinates and the pressure coefficients can be evaluated. In addition, the aerodynamic loads are then obtained by integrating the pressure over the lifting surface.
The literature refers to the articles of Hess and Rubbert [32,33] as pioneering implementations of the panel method in the aviation industry [34,92]. Because of their special ability in simulating realistic, complex aircraft geometries, and their accuracy, these source-doublet panel methods are well-known in studying aerodynamic research, including area of optimizing all rotor design parameters [50,93] and investigating component interaction effects [41,94]. Rotor performance in the vicinity of other rotorcraft components, such as fuselage and wings (for tilt rotors), is also important to ensure complete rotorcraft analysis. Not only is a rotor affected by other vehicle components, but it also causes many changes in the components' structural dynamic and control behavior. Therefore, many studies have been conducted to better understand interactive effects between rotor and airframe and to better predict rotorcraft characteristics.
Recent developments in vortex methods have enabled researchers to analyze full rotorcraft. Because rotor blades usually have low thickness sections, it is advantageous to use the VLM to simulate various blade shapes. The source-doublet panel method can adequately represent thick components, such as fuselage. One well-known example of using these methods for comprehensive analysis of full-configuration rotorcraft is CHARM, which stands for Comprehensive Hierarchical Aeromechanics Rotorcraft Model. Continuum Dynamics, Inc. [35,95] developed this comprehensive analysis model by combining several extended versions of their previous models; they later employed CHARM to analyze various rotorcraft phenomena such as rotor/airframe flow interaction [35], vibration [96], and noise propagation [97]. Other noteworthy combinations of VLM and the panel method for comprehensive analysis of rotorcraft can be found in Jang et al. [41] and Cao et al. [45], which investigate pressure fluctuation on fuselage surfaces and aerodynamic interaction among rotorcraft components, respectively. For more on coupled VLM and panel analysis, we encourage you to read the NASA Ames report by Wayne Johnson [98] on developing a comprehensive analysis.

Hybrid vortex methods
Although vortex methods, including LLT, VLM, and the source-doublet panel method, are useful approaches for analyzing rotor aerodynamics and simulating a complete rotorcraft configuration, they do not consider viscous and compressibility effects. These effects become more important than the influence of rotor wake when the rotorcraft operates in high-speed forward flight. CFD methods allow us to simulate a wide range of flow regimes and to accurately capture the complex flow physics occurring in nearfields around rotor blades. However, CFD simulations suffer from excessive numerical dissipation on coarse grids; hence, wake structure and vorticity tend to dissipate rapidly after shedding from rotating blades. Therefore, hybrid approaches combining the advantages of the vortex and CFD methods have been suggested to overcome these inherent shortcomings.
Advanced Rotorcraft Technology, Inc. developed a coupled vortex particle method (VPM) and CFD analysis model in which VPM and CFD are employed to resolve the wake vorticity dynamics and to compute the near-body flow solution, respectively [36,37]. VPM is used to solve the dynamics of wakes generated by rotors, wings, and other rotorcraft components, convect the wake vorticity under the combined influence of the freestream, and account for self-induced wake velocities and flow field perturbations caused by bodies like the fuselage. The CFD solver is used to compute the flow field around the bodies, including rotor wake-induced velocities. Yang et al. [99] developed another hybrid method, which combines Navier-Stokes equations near the blade to consider the viscous effects and near wake, and potential flow analysis in the farfield to model inviscid isentropic flow. Wie et al. [100,101] suggested an approach that reduces computational time by minimization of the computational domain: the strength and motions of the wake are modeled using the time-marching free-wake approach, and the rotating blades and flow field around the rotor are estimated using CFD. A moving overset grid approach was used to consider rotor motions during hovering and forward flight. The time-marching free-wake method provides inflow and outflow conditions in the CFD domain from an induced flow at each time step. A similar approach can be observed in the works of Shi et al. [102] and Zhao et al. [103]. In the present paper, a coupled VPM/CFD methodology is presented briefly, as follows.
Flow field velocity from the CFD solution, termed u CFD , is treated as an interference velocity in the VPM simulation. CFD flow field solution directly affects the rotor wake dynamics, represented by the vortex particles, as in Eqs. (29) and (30). The grid points and flow properties on the grid are derived from the CFD solution.
The governing equation of the near-body CFD solver can be changed to account for the influence of VPM wake induced velocity, as follows [104]: In the equation above, Q is the flow state variable; F i and F v are the inviscid and viscous fluxes. V is the considered fluid domain, S is the boundary of the fluid domain, and n is the normal vector pointing into the fluid domain. A flux correction term, ∆F i , is provided to account for the VPM driven velocity. Due to the coupled approach, the flux correction is dependent on both the VPM-induced velocity field and the CFD solver velocity field. The flow state variables for compressible flow are shown in Eq. (32).
In the above equation, u CFD = (u CFD , v CFD , w CFD ) is the flow field; ρ and E are density and total energy obtained from CFD simulation. In CFD solvers, for both structured and unstructured frameworks, the Flux Correction Approach can be used to account for the VPM-induced velocity. The implementation of VPM-induced velocity contributes to the CFD solver's correction of the convective flux term. The VPM-induced velocity correction terms can be formatted as follows: Here, u VPM and p are the VPM-induced velocity field and the pressure, respectively. The boundary condition on the body surfaces must be changed in the tightly coupled VPM/CFD solver. Due to the viscosity effect, the flow velocity on the body surfaces will be the same as the body velocity. The no-slip boundary condition was enforced in VPM/ CFD coupling in such a way that the combined velocity field, u CFD + u VPM , is the same as the body velocity, u body , as follows: According to the no-slip boundary condition, for a stationary body, u body = 0. Therefore, the boundary condition becomes u CFD = -u VPM in VPM/CFD coupling analysis. Similarly, the far-field boundary condition must be applied to the outer boundary of the CFD grid domain. Riemann invariants are commonly used to impose the far-field boundary condition. The far-field Riemann boundary condition can be modified to include VPM/CFD interaction and to account for the effect of the VPM rotor wake on the CFD far-field boundary. When the flow at the local grid possesses outgoing characteristics, the flow state at the boundary can be extrapolated from the inner CFD grid nodes. When the flow at the local grid has incoming characteristics, the flow state at the boundary can take the combined values from the free stream and the VPM-induced velocity data. For an incompressible flow, the freestream state is written as shown in Eq. (35), whereas for compressible flow it is a function of the local speed of sound.

Vortex methods for rotor wake dynamics
The wake field plays an important role in rotorcraft performance with significant impacts on rotorcraft's performance, vibration [105], maneuverability [106], and noise [54,105]. The dominant vortical structures of the rotor wake are the inboard vortex sheet and the tip vortices [25]. The vortices generated by a rotorcraft form different structures according to the flight condition. Compared to other wake components, the concentrated vortices evolving from the blade tips have the highest contributions to the flow field and blade aerodynamics [25]. Therefore, it is critical for the wake model to simulate tip vortices accurately.
Many experiments have been conducted to investigate disturbances and instabilities in rotor wake flow fields; resulting data have become the basis for several semi-empirical rotor wake models [107][108][109][110][111]. Wake models aim to predict the strengths and structures of rotorcraft wake vortex fields [112]. Choosing a proper wake model is an important responsibility for designers of rotorcraft. Therefore, researchers are continuously trying to improve the ability of wake models to accurately predict the aerodynamic loads of certain configurations. A common problem with the CFD method is that its inherent numerical dissipation causes the wake vorticity to diffuse too quickly. This causes problems in applying far-wake boundary conditions, especially for hovering flight, which are critical in wake modeling [112][113][114]. These difficulties of CFD, alongside its computational burden, reduce its application for practical problems. Therefore, vortex methods are still being implemented to model the rotor wake flow field [115] since they can offer accurate solutions with affordable computational burden and model the wake geometry without numerical dissipation error. References [9,34,116] provide exhaustive information on the history of wake models, as well as their development and application in rotorcraft studies. In this section, two wake models based on vortex theory, the freewake method (FWM) and the viscous vortex particle method (VPM), are introduced.

Time-marching free-wake method
In the free-vortex or free-wake methods, the wake vortex sheet can be divided into trailing and shed vortices, as shown in Fig. 3. The strength of a trailed vortex is the radial change between the bound vortices, whereas the strength of a shed vortex is related to the azimuthal change in the bound vortices with time [60,83]. The time-marching freewake method, in contrast to the prescribed wake model, allows the vortex sheet consisting of vortex filaments behind the rotor to move in free motion as the wake propagates downstream. In addition, the vortex sheet will deform at each time step based on the induced velocities by other vortex elements during unsteady rotor simulation [117]. According to the Kutta condition, which acts as a bridge between the wake field and the lifting surface vortices in the vortex methods, the strength of a nascent shed wake vortex element is equal to the strength of a blade vortex element placed at the trailing edge in Because the vortex surface has no force exerted on it, a number of discrete vortex elements are allowed to convect freely with the local stream velocity, which is induced by other vortex elements in the wake region. Biot-Savart's law in Eq. (36) provides the induced velocity components of a Lagrangian marker, which is connected to adjacent markers by straight or curved vortex filaments [60,118].
Here, V is the velocity of the marker induced by other markers, dl is the vortex line element, r is the distance between the Lagrangian markers, and Γ is the strength of the vortex filaments.
The positions of the Lagrangian markers are defined as a function of the blade azimuth angle (ψ) and wake age (ζ). Thus, the initial helical wake structure will change as the markers convect downstream. The rate of change of the Lagrangian markers' position can be expressed in the partial differential equation (PDE) form, which should be transformed into a finite number of finite difference equation (FDE) to obtain the wake solution by the time-accurate numerical means. Various time-marching approaches, such as Euler explicit, Adams-Bashforth, and two predictor-corrector schemes, have been employed for the time integration of the rotor wake equations. However, the free-wake method is sensitive to numerical instabilities, particularly in hovering and low-speed forward flight conditions [155] since the rotor wake is inherently unsteady and unstable [156]. The time-marching methods affect the overall stability of the free-wake analysis and the numerical errors could induce nonphysical disturbances. Bhagwat and Leishman [39] suggested a predictor-corrector central (PCC) and a predictor-corrector with second-order backward (PC2B) difference schemes. The PC2B utilizes the solutions at the three previous time steps for approximating the temporal derivative. This algorithm implicitly introduces additional effective damping terms, making it more stable with second-order accuracy. Kini and Conlisk [157] used a fourth-order Adams-Moulton method, which is implicit and computationally expensive for a free-wake analysis. Gupta and Leishman [158] performed both linear and nonlinear stability analyses of timemarching methods used in the free-wake analysis. Bagai and Leishman [159] developed a pseudo-implicit predictor-corrector method with a five-point central differencing approach for improving the convergence characteristics of the free-wake method. The numerical results obtained from free-wake method with a relaxation implementation were compared with the measurements of the tip vortex locations and flow field in hover and forward flight. The numerical instability of time marching free-wake modeling is also caused by the impulsive rotation at the beginning of computation, which creates a strong starting vortex. Initial treatments for this problem used a helicoidal spiral wake at the start [119] or assumed a uniform axial velocity [30]. Lee and Na [40,120] suggested that, for this problem, it is better to increase the rotating speed of the blade slowly from zero to the desired speed after some revolution. Their model achieved a close agreement with experiments [114]. Leishman and Bhagwat [39] provided more detailed information on this matter. Abedi et al. [160] suggested the VLM with the free-wake method to (36) predict aerodynamic loads of wind turbine rotor blades operating the unsteady flow field by employing tabulated airfoil data and the dynamic stall model. Wachspress and Yu [31] also proposed the lifting surface blade model with free-vortex wake model for comprehensive rotorcraft analysis, and Wachspress et al. applied to the problems of rotor/wake/ body interaction [35] and rotor/airframe noise of eVTOL aircraft [86] by incorporating fast panel methods and aeroacoustic software. Yeo [161] explored the performance potential of advanced compound helicopter configurations with a wide range of sizes, gross weight, and rotor systems using an aeromechanics analysis tool based on lifting line theory with the free-wake method.

Viscous vortex particle method
The potential flow assumption in the rotor wake limits the free-wake method's ability to forecast the rotor wake dynamics. The effects of turbulence on the diffusive characteristics of the vortex and the vortex stretching are considered empirically. The vortex decay factor or vortex core growth are modeled empirical formulations with the parameters, which are often derived experimentally [121]. Therefore, the numerical results of the free-wake analysis strongly depend on empirical formulations and the values of parameters. The viscous vortex particle method (VPM) solves the vorticity-velocity form of the incompressible Navier-Stokes equations with a Lagrangian description for obtaining the wake vorticity field. Therefore, consideration of viscous effect utilizing the VPM, which directly simulates the viscous flow and avoids artificial numerical dissipation, is crucial for solution accuracy. This method, which was developed for 2-D [122, 123] and 3-D viscous flow [123][124][125][126], has been used in rotor flow field analysis [42]. A rotor wake generated from the trailing edge of the rotor blade is modeled by a number of vortex particles, which influence each other and induce a velocity in the flow-field in the wake evolution process, as shown in Fig. 4. Compared to the vortex filament method, VPM has the advantage that the particles do not necessarily need to maintain connectivity with adjacent particles because particles propagate independently during the time-marching step. This property is especially useful for investigating wake interaction phenomena occurring with rotors or rotor-body configuration. In addition, VPM can be used to easily treat the wake vortex penetration problem. The Navier-Stokes equation can be represented in vorticity-velocity form with a Lagrangian description for incompressible flow. The Lagrangian description of the vorticity-velocity form is used in the VPM approach. The vorticity dynamics equation can be expressed as in Eq. (37), where ω is the vorticity and u is the velocity.
A natural way to solve the vortex transport problem is to use the vorticity-velocity form with a Lagrangian description. Grid generation is not required for the VPM simulation. Moreover, the convection term in Eq. (37) is not explicitly treated, resulting in a dissipation-free approach for the vortex particle model. To solve the vorticity-velocity governing equation, the vorticity field can be represented by a set of S Lagrangian vector-valued particles, as in Eq. (38). (

37) Dω Dt
= ω · ∇u + ν∇ 2 ω where x i and α i are the position and the vector-valued total vorticity inside the vortex particle with volume V i , respectively. The three-dimensional regularization function or smoothing function (ξ σ ) can be defined as in Eq. (39), where σ is the smoothing parameter [42]. A Gaussian distribution is one of the distribution functions that can be utilized, as expressed in Eq. (40) [43,127].
The velocity component of the i-th vortex particle (x i ) induced by other particles can be computed using Eq. (41), where σ ij is a symmetrized smoothing parameter used to conserve the linear and angular vortex impulses.
Here, ρ is a non-dimensional distance parameter, Κ(ρ) is the regularized Biot-Savart kernel used for the velocity calculation, and G(ρ) is Green's function, used for the stream function evaluation, as follows:

Fig. 4
Rotor wake modeling using vortex particle method (VPM) During the time-marching step for unsteady rotor simulation, the location of particles x i (t) will be updated using the local convection velocity, which is the sum of the freestream velocity, self-induced velocity, and wake-induced velocity. Then, the convection equation governs particle positions x i as follows: The governing equation for the vorticity field is defined in Eq. (37). The left-hand side of the equation is the material derivative of the vorticity. The first term on the right-hand side of the governing equation for the vortex dynamics represents the stretching effect, which describes vortex stretching and rotation owing to the velocity gradient. The viscous diffusion term is the second term on the right-hand side of Eq. (37); it describes the vorticity diffusion due to viscous effects. In the direct scheme, the vortex stretching effect is accounted for by directly multiplying the velocity gradient matrix by the particle vorticity [42].
The particle strength exchange (PSE) method [122,124,[128][129][130][131] can be utilized to consider the viscous diffusion effect, for which the second term is on the right-hand side of the vorticity governing equation in Eq. (37). The fundamental idea of the PSE algorithm is to approximate the Laplacian operator with an integral operator, avoiding numerical differentiation, which has lower numerical precision than the integral operation. The approximated Laplacian can be written as in Eq. (48).
The kernel η σ in the above equation is considered to have a Gaussian distribution function. The integral in Eq. (48) can be discretized across all particles by using the midpoint quadrature, resulting in the following equation: Koumoutsakos et al. [132] implemented the Neumann type vorticity boundary conditions to consider the no-slip state, which is expressed in terms of vorticity flux. Singh and Friedmann [133] used viscous VPM to simulate coaxial rotors in hover flight, which yields complex unsteady aerodynamic interaction effects. The flow separation during the dynamic stall is an important source of vibrations on a rotor at high advance ratios. They (44) also used VPM simulation to model the shedding of concentrated vorticity from the airfoil's leading edge and study the wake evolution of the coaxial rotors [44]. Su et al. [134] applied the VPM to an electrically controlled swashplateless rotor to investigate aerodynamic characteristics and wake structure. He et al. [42] investigated the effects of modeling parameters such as wake cut-off distance, time-step size, and other model tuning parameters for hovering and forward flight. Helicopter rotor loads were predicted using a source-doublet panel/VPM hybrid approach combined with a computational structural dynamics (CSD) algorithm. In a hybrid approach, the panel method was employed to simulate the rotor blade surfaces and close wakes, and VPM was used to model the far wake [91]. Alvarez and Ning [162] developed a viscous VPM code to study the unsteady wake dynamics of individual propellers of distributed propulsion electric aircraft, including near-far field transition and vortex breakdown. They presented its capacity to model rotor-on-rotor aerodynamic interactions in a side-by-side configuration [163]. Tan et al. [164] suggested coupling VPM and a discrete element method to simulate the helicopter brownout phenomenon and investigated the particle-surface interactions and the flow field of a helicopter in ground effect. Huberson et al. [165] employed the panel method with VPM and Farassat's formulation 1A to predict the vortex-solid interactions, such as helicopter blade-vortex interaction noise. While flow is passing over the rotor blade, the wake elements keep emitting from it and the number of vortex particles increases at every time step. Although the number of wake elements determines the prediction accuracy of the wake region, any increase in the number of vortex particles extends the computational time of the numerical model [114,120]. Therefore, an efficient summation rule needs to be incorporated for VPM simulation. In the literature, there are two broad groups of fast summing algorithms for resolving the N-body problem. The first is the Tree-Code approach [129][130][131]; the second is the fast multipole method (FMM) [135,136]. The detailed description is out of the scope of this paper. Berdowski et al. [166] developed the efficient framework coupled with actuator disc method, VPM, and open-source FMM library for handling the efficient data-parallelism on a CPU, and Willis et al. [167] suggested a pFFT-Fast Multipole Tree algorithm to accelerate the calculation of particle-induced velocity.

Nonlinear vortex lattice method
A comprehensive simulation tool capable of accurately and efficiently predicting rotor aerodynamics and wake dynamics is needed to design more advanced rotor systems. Among the various numerical approaches, VLM is one of the most suitable because it can provide reasonably accurate solutions with affordable computational cost and can represent the wake vorticity fields without numerical dissipation error. However, VLM simulation cannot consider nonlinear aerodynamic characteristics, which are mainly associated with viscosity, flow separation, and low-Reynolds number flow. To overcome these drawbacks, the authors of this paper have suggested a nonlinear vortex lattice method (NVLM) that combines VLM with airfoil lookup tables, semi-empirical models, and vortex strength correction. Moreover, NVLM is tightly coupled with VPM to simulate unsteady wake dynamics. Details of the numerical strategies used in NVLM/ VPM are elucidated in previous studies [83,137]. In the following sections, NVLM/VPM methodology is briefly presented, and its applications are discussed.

Airfoil look-up table
VLM simulation yields accurate results with an affordable computational cost for subsonic attached flow. However, it fails to predict the nonlinear aerodynamic behavior of rotor blades that occurs at high angles of attack (above stalling or critical angle of attack) or low-Reynolds number flow. For these flow regimes, an abrupt drop in the lift coefficient and nonlinear variation in lift coefficient with respect to the angle of attack become important in determining the rotor aerodynamics. The airfoil look-up table is an efficient and practical way of obtaining aerodynamic information about the airfoils of the rotor blades operating in both attached and stalled flow conditions. The sectional lookup table containing pre-computed airfoil aerodynamic data for a wide range of Reynolds numbers and angles of attack is used in the NVLM/VPM simulation to compute the lift and drag coefficients at each blade section's control point. The location of the control point should be explicitly specified to conduct an airfoil look-up table. The mathematical equation is derived by applying zero normal boundary conditions, and it is numerically solved using the root-finding approach. The result shows that the most suitable control point location is at half of the chord for NLVM [138]. A detailed derivation is outside the scope of this paper.
The aerodynamic coefficients of the individual blade sections rely largely on the local inflow velocity and effective angle of attack, which are calculated using Eqs. (50) and (51), based on VLM solutions.
Here, V inflow is the local inflow velocity. θ twist and θ pitch are the local twist angle and collective pitch angle at each blade section, respectively. a 1 and a 3 are unit vectors along directions tangential and normal to the rotating plane. Once the strengths of the bound vortices on the rotor blade are determined by solving the linear system of equations with an instantaneous boundary condition at each time step, the self-induced velocity (V ind,bound ) and wake-induced velocity (V ind,wake ) can be calculated by Biot-Savart's law. Finally, each blade section's sectional lift and drag forces can be calculated using Eqs. (52) and (53).
According to the force conventions for the rotor blade, all other aerodynamic coefficients, including normal (C N ), tangential (C A ), thrust (C T ), and torque (C Q ) can be evaluated using the lift (C L ) and drag (C D ) coefficients in conjunction with their reference angles. The normal and tangential force coefficients are determined relative to the chord line plane, while the thrust and torque coefficients are determined relative to the rotation plane.

Semi-empirical models for airfoil aerodynamics
The centrifugal and Coriolis forces caused by three-dimensional effects and rotational augmentation affect the stability of the boundary layer of the rotor blades, resulting in significantly higher lift coefficients than in the two-dimensional or three-dimensional non-rotating cases. These events are called as stall delay phenomenon, which is especially noticeable at the inboard section of the rotating blade, where they strongly influence the onset of flow separation. The boundary layer can be stabilized, minor separation bubbles can return to attached flow, and flow separation can be postponed above stall onset angle corresponding to the two-dimensional situation, despite the comparatively high angle of attack due to the low rotational speed at the inboard of the blade. As a result, the rotor blade in the inboard region experiences different post-stall airfoil characteristics than in the two-dimensional scenario. Several stall delay models have been developed to account for the three-dimensional rotating effect by adjusting the twodimensional aerodynamic coefficient data. Among several other models, Raj and Selig's model [139] can be utilized to account for the influence of stall delay on the aerodynamic coefficients at the inboard portion of the rotor blade; their model is an improvement of Du and Selig's model [140]. The corrected lift and drag coefficients, c l,3D and c d,3D , can be calculated in the following way: where c l,2D and c d,2D are two-dimensional lift and drag coefficients. ∆c l and ∆c d are the differences between aerodynamic coefficients obtained from the potential theory and the two-dimensional configuration. Incremental factors for the lift and drag coefficients, g c l and g c d , are determined as follows: where c/r and λ are the ratio of a local chord to the local radius and a local speed ratio. a, b, d, and n are empirical constants for the correction formula; their values were determined from the measurements. The dynamic stall phenomenon is also most likely to occur in the inboard regions of the rotor blade because the effective angle of attack is much higher than in the outboard regions due to the low rotational speed. This creates significantly unstable aerodynamic behavior in the blade section, resulting in a situation in which the aerodynamic coefficients enter a hysteresis loop. As a result of the presence of shedding vortices from the leading edge of the airfoil, nonlinear variations in lift, drag, and pitching moment coefficients as a function of angle of attack occur, and their values in the stall and post-stall regions are completely different from the static aerodynamic coefficients. Dynamic stall (54) generally causes unsteady and high loads on the rotor blades, which could cause structural damage to the rotor blade and other components. To forecast the aerodynamic loads acting on a rotor blade subjected to unsteady inflow circumstances, it is crucial to include the dynamic stall effect and the nonlinear aerodynamic features of the airfoils in the stall or post-stall regions. Various dynamic stall models have been proposed, including the ONERA model [141], the Leishman and Beddoes (L-B) model [142,143], and the Stig Øye model [144]. The L-B model, devised by Leishman and Beddoes, is the most widely used semi-empirical modeling method for the hysteresis loop of the aerodynamic coefficients. This model consists of three parts that were formulated to represent flow behaviors occurring around rotor blades: unsteady attached flow, trailing-edge separated flow, and dynamic stall flow. The details of the L-B model are omitted as it is rather complicated and lengthy.

Vortex strength correction
The sectional aerodynamic forces acting on the rotor blades for separated flow or low-Reynolds number flow can be evaluated through the look-up table. However, the bound circulation strength of the rotor blades, derived by solving a linear system of equations may be significantly over-predicted. Vortex elements placed on the trailing edge of the rotor blades will shed into the wake during time-marching simulation, and their strength will remain constant due to Helmholtz's theorem. As a result, the over-predicted strengths of bound and wake vortices can cause inaccurate estimation of the rotor wake geometry and evolution process, leading to errors in induced velocity predictions. Therefore, the influence of nonlinear aerodynamics on the bound circulation strength should be included to obtain a more accurate solution. In NVLM, sectional lift forces from the look-up table are employed to correct the bound vortex strength derived from VLM simulation. By equating the formulas for sectional lift forces based on airfoil theory in Eq. (52) and the Kutta-Joukowski theorem, the representative circulation strength at the control point, represented as Γ table , can be newly determined. A correction process with an under-relaxation factor is iteratively conducted until a convergence criterion is satisfied. If the difference between the current and updated circulation strengths reaches a value under 0.001%, then the iterative correction is stopped, and the corrected circulation strength is assigned to both the chordwise and spanwise vortex elements depending on the ratio of strength of each chordwise element to the average strength. A full description of the circulation strength correction technique can be found in the references [138].
After the bound circulation strength on the rotor blade converges, the corrected bound vortices located at the trailing edge will shed into the wake according to the Kutta condition. The rotor blade is rotated at each time step, and wake vortex elements are generated from the rotor blades' trailing edges. Once the strength of wake vortices is determined, the wake structure can be modeled through two approaches, vortex filaments and vortex particles, as mentioned above. In NVLM/VPM, the nascent wake that was recently shed from the full span of blades is represented using curved vortex filaments during approximately 4 ~ 5 discretized time steps. After that, the vortex filaments are split into a finite number of vortex particles, except for the nascent wake panels, to avoid wake instability problems. Vortex particles mutually interact with each other, allowing them to distort and propagate freely downstream with local convection velocity. A large number of vortex particles are generated during the time-marching step for wake evolution. The computing time to evaluate each vortex particle's local convection velocity increases exponentially as the number of vortex particles increases. To alleviate the computational burden, parallel computing using the Message Passing Interface (MPI) library on a multi-core processor is applied to evaluate the induced velocities at each vortex particle.

Applications
As previously mentioned, NVLM/VPM has been suggested to overcome the intrinsic drawbacks of the existing VLM, which is impractical for many applications. NVLM/ VPM has been applied to investigate the rotor aerodynamics, wake dynamics, and acoustics of various types of rotors, such as propellers [55,[145][146][147][148] and wind turbine blades [137,138,149,150]; its predictive capability has been compared with those of other numerical predictions and experimental results. The validation results show that NVLM/VPM can consider the nonlinear aerodynamic characteristics, which are mainly introduced by viscous effects and low Reynolds number flow. In this paper, specific application examples are introduced to assess the model accuracy of NVLM/VPM.

Caradonna-Tung rotor
The experiment on the Caradonna-Tung rotor model was conducted in 1981 [151], and provided extensive measurements of rotor aerodynamics and wake dynamics that have been widely used in the rotorcraft field to validate the accuracy of numerical methods. The model rotor is made of two blades installed on a tall column with a drive shaft. The Caradonna-Tung rotor is a rectangular blade without twist or swept angles, and its sectional shape is that of a NACA 0012 airfoil. The aspect ratio of the rotor blade is 6, the chord length is 0.1905 m (0.625 ft), and the rotor diameter is 2.286 m (7.5 ft), as listed in Table 3. A wide range of test parameters were used under ambient conditions, with a tip Mach number ranging from 0.226 to 0.890; the collective pitch setting varied from 0 to 12 degrees. The tip vortex trajectory was retrieved using a hotwire approach after pressure distributions were recorded at five cross-sections of the blade. The specific flow conditions chosen for the validation work are listed in Table 4.
The thrust coefficient is a non-dimensional parameter representing the aerodynamic load acting on the rotor blades in a direction normal to the rotating plane. Figure 5 provides a comparison between measurements and numerical predictions for the thrust coefficient distribution along the radial direction for various blade pitch angles. The results of VLM (blue dashed line with triangle symbols) and NVLM (red solid line with square symbols) are compared with experimental data to validate the prediction capability of the numerical models. Results obtained from both VLM and NVLM simulations matched well with the measurements, although VLM tended to slightly over-predict the sectional thrust force, especially in the blade tip region. The variation in the integrated thrust force coefficient with respect to collective pitch angle is shown in Fig. 6. Exact agreement between the experiment and NVLM was observed, whereas slight overprediction was observed in the case of VLM.
The rotor wake is defined as an unsteady fluctuation flow that generates unsteady rotor blade aerodynamics and complicates the flow field. The wake vortex particles have different sizes depending on the circulation strength. The stronger the wake vortex strength, the larger the particle size. The color of the vortex particles also varies with their circulation strength, just as the size does. The evolution of wake geometries of the Caradonna Tung rotor with respect to revolutions is shown in Fig. 7. The development of near and far wakes is easily predicted as rotor revolutions increase; tip vortex descent and wake contraction can be clearly noticed. Figure 8 shows a comparison of vorticity magnitude contours on a vertical cross-section (x-z plane) through the center of the rotating axis with an increase in revolutions. It can be clearly seen that the periodic shedding of wake vortices behind the rotor plane gives rise to a symmetric wake structure. A comparison of the tip vortex trajectory locus between experiment and computation is provided in Fig. 9, which shows an excellent agreement with the experiment.

NREL Phase VI rotor
Vortex methods also have broad applicability in simulating wind turbines, which operate at lower Mach numbers than rotorcraft. Here, we again validate the proposed   NVLM against extensive and high-quality measurements of an NREL Phase VI wind turbine operating in both axial and yawed flow conditions. A well-known experimental study was conducted in the NASA-Ames wind tunnel facility under controlled wind conditions [152]. The NREL Phase VI wind turbine consists of two-bladed rotors without hub tilt, coning, and prebend angles. The rotor blade is an S809 airfoil with a tapered-twisted configuration. Details of the NREL Phase VI model and its operating conditions for the validation work are provided in Tables 5 and 6, respectively. Because the NREL Phase VI is a stall-regulated wind turbine, the rotor blades are designed so that flow begins to separate from their upper surfaces at high wind speeds. The aerodynamic thrust and power output are controlled by stall effects occurring on the wind turbine blades. Below a wind speed of 10 m/s, the flow remains fully attached over the rotor blade. However, with increases in wind speed, separated flow starts to cover the rotor blade from the inboard to outboard regions. Finally, the flow is completely separated from the rotor blade at a wind speed of 20 m/s, and massive flow separation occurs at a wind speed of 25 m/s. Once flow separation occurs, the flow  field around the rotor blades becomes highly unsteady and transient. Figures 10 and  11 respectively provide comparisons of the normal and tangential force distributions along the radial direction, depending on the wind speed. Above a wind speed of 10 m/s, there is a distinct difference between the results of VLM and NVLM simulations; significant improvements associated with flow separation and the stall effect are evident. NVLM predictions are seen to be quite close to the measurements, even if there are minor discrepancies, whereas the VLM simulations show significant overprediction of the normal and tangential force coefficients due to the neglecting of the nonlinear aerodynamic behaviors caused by the flow separation. In Fig. 12, the overall aerodynamic performances of wind turbines in terms of low-speed shaft torque and  thrust force predicted by VLM and NVLM simulations are compared with experimentally measured data. Above a wind speed of 10 m/s, VLM tends to overestimate the aerodynamic load on the rotor blades, while the NVLM results are in reasonably good agreement with experimental data and the CFD results, even if flow separation occurs at high wind speed. Figures 13 and 14 show NVLM predictions of unsteady aerodynamic loads on a wind turbine blade subjected to wind speed of 7 m/s and yaw angle of 30°. Under yawed flow conditions, the wind turbine blades suffer from cycle-to-cycle variation in aerodynamic load. This is mainly attributed to the advancing and retreating blade effect and the skewed wake effect. It can be observed that azimuthal variations in the sectional aerodynamic loads are much more pronounced at the inboard section because, there, the rotor blade experiences significant variation in angle of attack due to asymmetric inflow distribution. The time-averaged normal and tangential force coefficients of the rotor blade are shown in Fig. 15. Comparing Fig. 15 (a) and (b), we can confirm that NVLM also accurately predicts wind turbine aerodynamic loads with yaw angle. Figures 16,17,18 compare the results of NVLM and LLM [153] with NREL measurement data obtained under wind speed of 10 m/s and yaw angle of 30°. At higher wind speeds, advancing and retreating blade effects become more dominant. This can induce a dynamic stall of rotor blades; the blades experience periodic variation in angle of attack with large amplitude, particularly at the inboard section. These results indicate that NVLM can provide more accurate predictions of aerodynamic load from the inboard to the outboard regions than LLM, which tends to overestimate the aerodynamic loads. Based on the authors' experience, it is difficult to use the LLM approach to accurately predict tangential forces of rotor blades because a single vortex element in the chordwise direction cannot sufficiently represent the three-dimensional geometry and various planform shapes of the rotor blades [31,138,154].

Conclusion
This paper provides a brief introduction to the most common vortex methods for analyzing rotor aerodynamics and wake dynamics. The purpose of this paper is not  to present an exhaustive review but to present state-of-the-art vortex methods and address well-known uses of these methods for simulating flow over rotorcraft. The vortex methods are still appealing today due to their negligible numerical dissipation, conservation of flow invariants, relaxed stability condition at time steps, and ability to capture high-resolution wake structure. The vortex methods are coupled with grid-free wake modeling methods, such as the time-marching free-wake method and viscous vortex particle method, which do not require any grid generation effort and minimize the dissipation of vorticities over long distance traveled. Moreover, the vortex methods are useful for preliminary designs and parametric studies because they produce numerical results much more quickly than grid-based CFD simulation. However, the numerical methods based on vortex theory fail to provide accurate representations of the viscous boundary layer and lead to underestimation of the drag force due to the assumption of potential flow. They are also unable to consider nonlinear aerodynamic characteristics of airfoils involved in the rotor blade, including viscosity, flow separation, and low-Reynolds number flow. To overcome the intrinsic shortcomings of VLM, the authors have proposed NVLM; its applications to rotor aerodynamics and wake dynamics were discussed. Simulations indicated that NVLM/VPM has great capability to assess aerodynamic loads acting on rotor blades for a wide range of operating conditions and to simulate the generation and evolution of rotor wake, allowing for higher resolution simulation of the helical wake structure.
Recent increases in efficiency of electric propulsion, particularly in areas of motor and battery technology, have been driving the development of electric vertical takeoff and landing (eVTOL) aircraft, which use electric power to hover, take off, and land vertically in highly populated urban areas. Among various configurations, eVTOL aircraft designed with multiple lifting rotors or prop-rotors are popular in the UAM market because distributed electric propulsion (DEP) systems using multiple rotors can improve safety and reduce noise. With the increasing number of rotor systems, vortex methods have emerged as useful tools for comprehensive analysis of eVTOL aircraft; these methods can provide a comprehensive solution with an affordable computational cost by combining structural, flight dynamic, and acoustic analysis solvers. Furthermore, in the presence of fuselage or any other such body, a hybrid method that combines vortex methods with CFD or a generalized treatment of boundary conditions on solid walls can be used, an efficient and accurate way to compute flow fields around bodies while considering the effects of the rotor wake.