Viscous Vortex Particle Method Coupling with Computational Structural Dynamics for Rotor Comprehensive Analysis

: A panel/vortex particle hybrid method is coupled with a computational structure dynamics code to predict helicopter rotor loads. The rotor blade surfaces and near wakes are modeled by the panel method, while the far wake is modeled by resorting to the vortex particles method. A fast summation method is introduced to accelerate the evolution of particle–particle-induced velocity and its derivative as well as panel–particle interactions. The developed vortex particle method code is coupled with the multibody code MBDyn to predict the rotor airloads. Numerical validations are carried, out and the results are compared with the experiments and simulation results in the literature.


Introduction
Rotorcrafts are characterized by strong interactions between the slender blades' structural dynamics and the unsteady aerodynamics environment; among them, the blade vortex interaction (BVI) causes significant noise and vibration problems. These strong interactions present a substantial challenge for the prediction of airloads. A high fidelity aeroelastic simulation, coupling rotor aerodynamics with structure dynamics, is often required. Research efforts have been undertaken to develop sophisticated and comprehensive analysis methods for rotorcrafts coupling structural [1,2] and aerodynamics models [3,4]. The structural models, also referred to as "computational structural dynamics models" (CSD), perform sophisticated structural dynamics simulations by resorting to nonlinear beam models and different joints in a multibody dynamics framework and carries out the trim calculation. Several CSD codes have been developed during the past decades; among them are CAMRAD II [5], DYMORE [6], RCAS [7], UMARC [8], and MBDyn [9]. Usually, CSD codes are equipped with simplified aerodynamics models such as the lifting line model; these models are not suitable for an accurate prediction of unsteady airloads. In recent years, significant progress has been made in the use of coupled computational fluid mechanics (CFD) and CSD analysis to calculate the trim parameters and to predict the airloads and structural responses [4,[10][11][12][13]. CFD methods such as Reynolds-averaged Navier-Stokes (RANS) or large-eddy simulations (LES) for rotorcraft analysis usually make use of the overset mesh methodology; provide a high-fidelity simulation; and are capable of predicting the transonic flow (shock wave) effects, separated flow, and reverse flow phenomena. However, these grid-based methods may suffer from excessive dissipation due to the numerical discretization unless a fine grid system and advanced meshing techniques such as adaptive mesh refinement (AMR) that allows us to adaptively refine the discretization around key regions of high vorticity are adopted. This makes the CFD methods more computationally expensive. Hence, a rapid assessment of aerodynamics that, at the same time, provides an acceptable fidelity may be of interest. As such, an accurate wake model is essential for rotorcraft aerodynamic models. The Lagrangian vorticity-velocity form of the incompressible Navier-Stokes equations is a natural choice in order to solve the vortex wake transport problem. The viscous vortex particle method (VPM) solves the vorticity transport equation using a grid-free Lagrangian formulation. The VPM represents the wake with a cloud of vortex particles; each particle is associated with vector-valued vorticity, velocity, and position. Based on first principles, VPM captures, without introducing any artificial numerical dissipation, the fundamental fluid flow physics for both the vorticity stretching and the physical viscous diffusion. In VPM, the aerodynamic surfaces are treated as a vorticity source that emits vortex particles into the flow field. Opoku et al. [14] modeled the far-field wake with vortex blobs, coupled with a panel method, and incorporated Fowcs Williams-Hawkings (FW-H) equation to predict the blade-vortex interaction (BVI) and aeroacoustics of a rotorcraft; however, the viscosity of the flow was not accounted for in this model. Willis et al. [15] developed the FastAero code, a fast, unsteady panel method with vortex particle wakes, and used the pre-corrected Fast Fourier Transform (FFT) approach and the fast multipole method (FMM) to accelerate the computation. Similar works have been conducted by Kenyon and Brown [16], He and Zhao [17], Tan and Wang [18] and Wang et al. [19], who combined the potential flow panel method or the lifting line model with the viscous vortex particle method. Tugnoli et al. [20] developed a flexible medium fidelity aerodynamic computational tool, called DUST, which enables different aerodynamic elements such as lifting line, vortex lattice, and vortex particles, all available in a single model. In the above works, either the fast multipole method (FMM) [21] or the Treecode method [22] were utilized to accelerate the evaluation of particle convection velocity as well as its gradient. However, when complex configurations are analyzed, the number of vortex panels increases and the computational cost from the vortex panels to the particles becomes nontrivial. For this reason, a fast summation method is required also for the vortex panels. Alvarez and Ning [23] developed a parallel viscous vortex particle method called FLOWVPM, based on the open-source code ExaFMM developed by Yokota and Barba [24] and Hu et al. [25], to evaluate the multi-rotor configuration performance. In the work of Alvarez and Ning [23], the load distribution was calculated by resorting to the blade element method (BEM) and was used to derive the circulation along the blade in order to perform all the computations using FMM, to model the lifting surfaces using embedded vortex particles, and thus to obtain an O(N) simulation. However, obtaining the surface pressure of the lifting surfaces, which is of great importance for aeroacoustics analysis, is not straightforward.
In present study, a panel method is coupled with VPM. The panel method is used to model the rotor blade surfaces and near wakes, and VPM is used for far wake. Besides the classical FMM, a semi-analytical fast summation method is developed to speed up the computation of panels' mutual interactions and of the influence exherted from the panels to the particles. The developed VPM code is coupled with the multibody code MBDyn to predict the blade unsteady aerodynamic loads. Numerical validations are carried out, and the results are compared with the experiments and simulation results in the literature.

Panel/Viscous Vortex Particle Aerodynamic Model
Consider the unsteady incompressible flow around the rotor configuration. According to the Helmholtz decomposition theorem, the divergence-free velocity field u can be decomposed as u = u ∞ + u φ + u ω , where u ∞ denotes the free stream field; u φ is the irrotational part induced by lifting surfaces, is defined by scalar potential φ as u φ = ∇φ, and is the solution of a three dimensional boundary element method; and u ω is the rotational part due to the wake, defined by a vector potential ϕ as u ω = ∇ × ϕ.

Panel Method
The irrotational, inviscid, incompressible potential flow is characterized by the scalar velocity potential φ, which is governed by the Laplace equation The boundary condition is derived from the requirement that the normal velocity on the boundary surface should be zero. A regularity condition at infinity is also required: The potential flow is solved by resorting to a low-order source-doublet panel method. The body surface and the wake surface are divided into quadrilateral panels with constant strength. A constant source strength on each body surface is defined as where u B is the velocity of a point on a body surface and n denotes the surface outward normal vector. Applying the internal Dirichlet boundary condition stating that the perturbed potential just inside the body surface is zero yields the linear equation system for the unknown doublets distributed on the body surface: where A k , B k , and C k are the potential influence coefficients of the body doublet, body source, and wake doublet, respectively, and N B and N W denote the number of body surface panels and wake surface panels, respectively. The doublets of wake panels are known from the previous time step, with the exception of those adjacent to the trailing edge, for which the doublets are determined by applying the so-called Kutta condition. The pressure distribution on the body surface is computed from the unsteady Bernoulli equation. The nondimensionalized pressure coefficient can be calculated at any point as where u is the local fluid velocity and u re f is the local body surface velocity. The scalar velocity potential φ is, for the current method, equal to ∂φ ∂t = − ∂µ ∂t .

Vortex Particle Method
It is well established that the wake has a large influence on the aerodynamics of rotorcrafts [3]. As such, an accurate wake model is essential for the rotorcraft aerodynamic model. In the present work, the wakes are modeled using the viscous vortex particle method, which is a Lagrangian formulation of the incompressible Navier-Stokes equation: where ω is the vorticity, dω dt is the material derivative, and υ is the kinematic viscosity. The vorticity field ω is represented by N vortex particles with a vector-valued vorticity strength and by using moving Gaussian basis functions ξ σ as follows: where α j is the total vorticity inside the vortex particle, x j is the location of the j th particle, V j is the particle volume, and σ j is the particle core radius. Here, ξ σ (x) = ξ(|x|)/σ 3 , and ξ(|x|) is the Gaussian basis function that is used to smooth the vorticity field and deprecate singularity and is defined by With the discretized vortex particles, the viscous vorticity dynamic equation can be rewritten as follows: Equation (10) describes the vortex particles convection; in Equation (11), the first term on the right hand-side is the stretching-effect term, which stands for vortex stretching and rotation due to the velocity gradient, while the second term on the right-hand side is the viscous diffusion term, which represents the vorticity diffusion caused by viscous effects. The velocity of particles can be calculated by the well-known Biot-Savart law, which is the solution to the unbounded Poisson equation: where is the regularized Biot-Savart kernel for velocity evaluation, and G(ρ) is the Green's function for vector potential and can be expressed as follows: where is the error function. The vortex particle method was shown to convergence as the number of particles increases provided that the smoothing core radius is larger than the distance between neighbor particles. The stretching term is handled by using the so-called direct scheme, which is the dot product between the velocity gradient and the particle vorticity; the velocity gradient is obtained by analytically differentiating the velocity field, formulated as where α j × is the skew symmetric matrix built from vector α j . The viscous diffusion term is accounted for by using the Particle Strength Exchange (PSE) method, which approximates the Laplacian operator with an integral operator with a Gaussian kernel η σ ij , thus obtaining higher accuracy compared with numerical differentiation:

Hybrid Method
In the present method, almost all the wake panels that one would use in a normal panel method are replaced by vortex particles. At each time step, a new row of wake panels are emitted from the trailing edge, which is used to define the scalar potential jump, and their strength is determined to satisfy the Kutta condition. A second row of wake panels is generated by convecting the previous time step first row with the local velocity and by keeping their strength unchanged and is prepared to be converted into vortex particles at next time step, as depicted in Figure 1. Thus, there are always two rows of wake panels after the first time step, with the remaining portion of the wake represented by vortex particles. The conversion of doublet wake panels to vortex particles follows Hess's equivalence principle for which the velcocity induced at any point by a dipole with strength µ is equivalent to a surface vorticity γ = n × ∇µ and a boundary line vortex µ. For the loworder panel method, the dipole strength is constant within each panel; thus, the surface vorticity is zero and the dipole is equivalent to a vortex ring around the perimeter of the panel. The wake panel is transformed into a vortex particle in such a way that the zero and first moments of vorticity are conserved, as illustrated in Figure 1. It is noted that, in the present work, only three edges of the wake panel are integrated into the vortex particle while the fourth edge, which is adjacent and parallel to the second row wake panels' trailing edge, is treated as a vortex filament. The remaining vortex line can be regarded as modeling the circulation jump at the edges of the wake panel of the last row or as canceling out the circulations on the vortex particles so that Kelvin's circulation theorem is met [15,19]. When referring to the particle velocity, the contribution from panels must be accounted for. This can be evaluated by the theory of basic singularities. Meanwhile, the stretching effect from panels to particle can be considered by differentiating the velocity field. The particles also have an influence on each panel. In a traditional source-doublet combination panel method, the particles influence is added into the freestream conditions to solve for the source strength for each panel: where σ i and n i are the source strength and normal vector of the ith panel, respectively, and u i,j is the induced velocity of particle j to panel i.

Fast Summation Method
As shown in Equation (12), calculation of the convection velocity requires, for each particle, summation over all the particles within the computation domain as well as the velocity gradient, which renders the vortex particle method a typical N-body problem with a direct summation computation complexity of O(N 2 ) at each time step, where N is the number of particles. Two effective acceleration approaches are the TreeCode [22] and the Fast Multipole Method (FMM) [21,26], which reduce the computational complexity from O(N 2 ) to O(N log(N)) and O(N), respectively. Both approaches utilize an Octree data structure to hierarchically subdivide the whole computation space occupied by particles into three dimensional cells and to exactly calculate the near-filed interactions by particleto-particle (P2P) interactions while approximating the far-field interactions of the particles that are well separated. The separation threshold is defined by the Multipole Acceptance Criteria (MAC) r 1 +r 2 R < θ, where r 1 and r 2 are the radius of the source and target cell, respectively, and R is the distance between the two cells. However, FMM can approximate far-field interactions to a specified tolerance. TreeCode approximates the interaction between well-separated particles using the interactions between particle and cells by multipole expansion, which is a multipole-to-particle (M2P) interaction. In FMM, for two well-seperated cells, as depicted in Figure 2, the influence of the source cell on the center of the target cell is first calculated by a multipole-to-local (M2L) transformation, the influence can be propagated to offspring cells of the target cell by local-to-local (L2L) transformation down to the finest level, and finally the velocity of each particle in the target cell can be evaluated by accumulating the local-to-particle (L2P) interactions at the finest level. In the present work, an adaptive Octree is built, the dual tree traversal algorithm is implemented to identify the well-separated interaction list, and the classical FMM in spherical coordinate is adopted to speed up the calculations of the convection velocity and stretching term for the particle-particle interactions. Meanwhile, as in [17] and due to the rapid decay of the kernel function, only the particles close to the target particle contribute significantly to the PSE for dα i /dt; thus, the PSE can be considered only in the process of the particle-to-particle (P2P) interaction of FMM.
In the numerical experiments, as the rotor wake develops, the number of vortex particles increases rapidly and the vortex panel-particle interaction becomes a time-consuming computation step. A fast panel-particle evaluation is thus required as well. The velocity potential at x j generated by a constant dipole quadrilateral and source quadrilateral can be expressed as Compared with the potential generated by point sources and point dipoles, one needs to compute an integral over the panel. When the target point is outside the panel, the integral is regular and can be computed accurately by a numerical method, such as a Gauss quadrature rule. When the target point is inside the panel, the integral is singular. Thus, specific approaches are needed depending on the location of the target point. In the present work, a semi-analytical method is adopted, where the inside integral is computed by resorting to an analytical method and the outside integral is evaluated by Gauss quadrature. For the inside integral, analytical expressions are available; see, e.g., [27]. When the target point is located outside the panel, as depicted in Figure 3, the integral over the quadrilateral panel is evaluated numerically. By applying the Gauss quadrature on the panel surface, Equation (18) yields where x k is the Gauss quadrature point on the element; n is the element normal vector; σ and µ are panel source and dipole strength, respectively; and w k is a suitable integration weight. In order to compute the integration weight, the coordinates on the panel are interpolated by linear Lagrangian polynomials N j (ξ, η) by x(ξ, η) = ∑ j N j (ξ, η)x j , where ξ and η are the local coordinates in the local panel coordinate system and span the range [−1, 1]. An integral´A f (x)dA over a quadrilateral element can thus be approximated, as customary in finite element modeling, aŝ where w 0 k is the Gauss quadrature weight. Thus, w k = w 0 k ∂x ∂ξ × ∂x ∂η |ξ k ,η k . From Equation (19), it can be seen that the quadrilateral element can be represented by several particles located at the Gauss point and with a weighted strength. This makes it possible to implement the panel FMM into the classical particle-based fast multipole method. The multipole expansion for source and dipole quadrilateral elements can be obtained: whereM m n andM m n are multipole expansion coefficients for the source and dipole elements, respectively; ρ, α, and β are spherical coordinates; and Y m n is the spherical harmonics that are consistent with those defined in Cheng et al. [26]. Compared with a standard particle, the panel itself has a finite area that could not be ignored. When constructing the Octree, some panels may belong to several boxes, making it difficult to evaluate the MAC of a panel. In the present method, the panel is firstly represented by a point located at it center, and the Octree is built with the classical procedure. After this, following a post-traversal algorithm, the box center is kept unchanged and its radius is increased in order to occupy the entire panel for the leaf box or its children box for a non-leaf box. In the modified FMM, the multipole and local expansion coefficients have 5 components: 3 for the vector-valued vorticity, 1 for the weighted dipole strength, and 1 for the weighted source strength. A case with randomly distributed particles within a unit cubic was studied to demonstrate the convergence and acceleration of the FMM. The numerical experiments were performed on an Intel Xeon CPU E5-2640 0 @ 2.50 GHz desktop. Figure 4a shows the L 2 -norm of the relative difference between the induced velocity computed by the FMM and the result computed by the analytical method. It can be seen that, when the truncation order is larger than 5, a relative error of 10 −4 can be obtained. Figure 4b shows the time consumption of the two methods; the FMM shows its advantage for an increasing numer of particles, as expected, and when the particle number is 10 6 , an acceleration of over 100 times is obtained. It is clear that the computational complexity is proportional to the number of particles, O(N), as it should be. In order to verify the validity of the proposed approach for the panels, the interaction of n × n × n panels that are uniformly distributed in the unit cube was studied and the velocity measured at the center of each panel was compared with the analytical results computed with the formulae from [27]. Figure 5a shows the relative errors obtained with 30 × 30 × 30 panels. Figure 5b reports the computational time cost and shows that a nearly linear computation complexity was obtained.  The overall speedup allowed by applying the FMM also to panels is significant for real-word simulations. As an example, the fully coupled aeroelastic simulation of the rotor of Section 5.4 takes 224 min for ten rotor revolutions if the FMM is applied to the particles alone, but the overall computational time is reduced down to 85 min after applying it also to the panels.

CSD Model
MBDyn is a general purpose multibody dynamics code developed at Politecnico di Milano [9]. MBDyn solves initial value problems (IVP), be them either static or dynamic. The constrained dynamics is formulated as where M is the mass matrix, β is the momentum, f is a vector containing generic force and moment, and λ φ and λ ψ are the Lagrange multipliers associated with holonomic and nonholonomic constraint equations. The differential-algebraic equation is linearized and integrated in time by using an implicit A/L stable linear multi-step integration scheme and solved by a Newton-Rapson iterative procedure. Based on the Cosserat beam formulation, an original geometrically exact, three-node finite volume beam formulation that yields equilibrium of finite portions of beam was implemented in MBDyn as depicted in Figure 6 [28]. The internal force was computed from the nodes kinematics through the constitutive law and contributed to the equilibrium of the beam nodes. The beam cross sectional shear and warping were also accounted for by the constitutive properties of the beam cross section. The cross-sectional properties could be computed either by the variational asymptotic methods [29] or semi-analytical methods [30,31].

Aeroelastic Coupling
The data exchange between the aerodynamics module and structure module was accomplished by using local UNIX sockets. A loose coupling scheme was adopted. At each time step, the structure module predicts the kinematics by extrapolation and sends the predicted kinematics to the aerodynamics module; the aerodynamics module deforms the configuration based on the predicted kinematics, computes the aerodynamic forces, and sends them back to the structure module; the structure module receives the forces and iterates until convergence; and, after this, it advances to the next time step. The aeroelastic interface between the structural nodes of the beams and the aerodynamics surfaces was computed by an energy-conserving moving least-squared algorithm [32].

Numerical Examples
In order to validate the aerodynamic model, different helicopter rotors that have been tested by experiments under hover and forward flight condition were studied. Finally, the aerodynamic model was coupled with the CSD module and the HART II model was studied. In the present work, the rotor wake was cut off at a distance from the rotor equal to four times the rotor radius.

Caradonna-Tung Rotor
The two-bladed Caradonna-Tung rotor [33] in hover was studied here. The blade was rectangular, untwisted, and untapered and had a radius of 1.143 m and a chord of 0.1905 m. The sectional airfoil was NACA0012. The root cutoff used in simulation was equal to the chord length. The rotation speed was 1250 rpm, and the tip mach was 0.439. A total of 60 panels were arranged along the chordwise direction, with 20 panels in the spanwise direction on each blade. An azimuthal step of 5 • was adopted.
In Figure 7, the predicted thrust coefficient at a collective of 5 • , 8 • , and 12 • was compared with the experimental data. Figure 8 shows the predicted tip vortex location in both the radial direction (r/R) and the axis direction (z/R) variation as a function of the vortex age at a collective of 8 • . Figure 9 shows the pressure coefficient at span locations of 50%R, 68%R, 80%R, 87%R, and 96%R. The results obtained by the present method were compared with both the experimental data and the numerical results that were obtained with a code coupling CFD and free wakes. It can be seen that good correlations exist between the present method and both the hybrid CFD method and the experiment.

AH-1G Rotor in Forward Flight
The AH-1G rotor [34] in forward flight was simulated. The rotor had two blades with the optimal load survey/tip aerodynamics and acoustics test (OLS/TAAT) airfoil. The rotor radius was equal to 6.7 m, and the blade aspect ratio was 9.8 and had a linear twist from 0°to −10 • from root to tip. Each blade was modeled with 2400 panels, with 60 panels along the chordwise direction and 40 panels along the spanwise direction. The rotor was operated at an advanced ratio of 0.19, the tip mach number was equal to 0.65, and the thrust coefficient was equal to 0.00464. An azimuthal angle step of 3 • was used in the present simulation. Figures 10 and 11 show the pressure coefficient at a span location equal to 60%R and 96%R and at azimuth angles of 30 • , 90 • , 270 • , and 300 • . Figure 12 illustrates the variation in the sectional thrust coefficient at span positions of 75%R, 87%R, 91%R, and 97%R over one revolution. The predicted results agree reasonably well with both the experimental data and the results computed by a CFD simulation [35]. The BVI influence on the normal force was captured by the present method on the advancing side at azimuthal angle of 50-100 • and on the retreating side around an azimuthal angle of 270 • .

CSD
To validate the structural model used in the present work, a clamped flexible blade, shown in Figure 13 and taken from Cheng [36], was first simulated. The results computed by the present method were compared with Dymore and with the Geometrically Exact Beam Theory (GEBT) code [37]. The blade was 1 m long and subjected to a transversal periodic force f = 10 sin(50t) N at the free tip. Since the geometric and material properties were not provided in Cheng [36], the cross-sectional properties reported in Cheng [36] were used directly. The properties are reported in Table 1 for completeness. The blade was discretized by means of 20 three-node finite volume beam elements. Table 2 lists the fundamental frequency results, and Figure 14 plots the axis displacement, transversal displacement, and flap rotation history at the blade tip. It can be seen that the current result has a good coincidence with those in Dymore and GEBT. Note that the maximum flap rotation was about 47 • and that the maximum tip deflection reached about 0.5 m, 50% of the beam span. In a nutshell, the beam had definitely large rotations and deformations.

HART II Rotor Coupling with CSD
In this section, the aeroelastic simulation was validated against the experiment data from the well-known Second Higher Harmonic Control Aeroacoustic Rotor Tests (HART II) performed in the large low-speed facility of the German-Dutch Wind Tunnel (DNW) [38]. The HART-II test used a 40 percent mach scaled, four-bladed hingeless BO105 model rotor. The blades were rectangular, with a radius of 2 m, a chord of 0.121 m, a −8 • linear twist, and a precone of 2.5 • . The HART II tests were operated at descent conditions, where strong blade-vortex interactions (BVI) may occur. A detailed description of the experiment parameters can be found in van der Wall and Burley [39] and Smith et al. [40]; the blade properties are listed in Table 3 for completeness. The center of gravity (CG) and the tension axis (TA) were measured from the elastic axis (EA), while the elastic center was measured from the quarter chord; the positive direction was toward the leading edge. The inertial properties were measured with respect to the center of gravity.  Validations were performed for the baseline rotor at an advance ratio of 0.151 and a thrust coefficient equal to 0.00457. The blade was discretized by 40 finite volume beam elements, and 20 elements were placed in the blade inboard region from the hub center to r/R = 0.22, where a significant variation in blade cross-sectional properties took place. The region from r/R = 0.22 to tip had a constant mass distribution and structural properties and was discretized by means of 20 equally spaced beam elements. In order to model the pitching control system stiffness, as stated in [13], a torsion spring was modeled at the pitching hinge and the spring stiffness was tuned to match a non-dimensional first torsion frequency of 3.845 /rev, which was reported in the HART II document. At first, the fan plot was predicted by the current CSD method for the HART II rotor blade to verify the structure modeling, as shown in Figure 15. The dominant natural mode shapes, i.e., the first three flap shapes, the first two lead-lag shapes, and the first torsion shape, were also compared at 100% rotational speed with those from the HART II documents [39,41] and with those computed by UMARC [42]. As seen in Figure 16, the present results correlate well with those computed by UMARC; the difference in the torsion mode at the root was due to the different boundary condition applied at the pitching hinge: in the present simulation, a lumped torsion spring was used, and this led to a twist discontinuity at the pitching hinge.
In addition, as seen from Figure 15, the present model is able to predict transitions between modes of vibration. In general, the predicted mode shapes are in good agreement with those from the HART II document data. In the simulation, each blade was modeled by 40 × 20 panels and an azimuthal angle step of 3 • was used. The rotor was trimmed in order to reach the experiment thrust coefficient and zero longitudinal and lateral moments at an effective shaft tilt of 4.5 • using the Newton-Raphson iterative method [35,43]. The trimmed parameters are listed in Table 4.
The azimuth-wise variation for the sectional normal load C N M 2 at the 87% span-wise location is shown in Figure 17 and is compared with the experiment results and CFD/CSD coupling results from the work of Biedron and Lee-Rausch [44] (FUN3D/CAMRAD-II), Boyd [12] (OVERFLOW2/CAMRAD-II), and Kim and Kwon [45], which couples KFLOW and a nonlinear Euler-Bernoulli beam CSD solver. The airloads show good correlation with the experimental data and the CFD/CSD results, and from the azimuthal angle of 0 • to 75 • at the advancing side and from 270 • to 360 • at the retreating side, the magnitude and location of the unsteady loading (BVI) events are well captured. C N M 2 is computed by integrating the lower surface pressure and upper surface pressure along the chordwise direction, following the method described by Swanson et al. [46] and Johnson [47], that is where N is the sectional normal force, a is the air speed, c is the chord length, and P L and P U are upper and lower surface pressure, respectively.  Figure 17. HART II baseline case.

Conclusions
In this paper, a panel/viscous vortex particle hybrid method was developed and was coupled to the multibody dynamic code MBDyn to predict the blade unsteady aerodynamic loads. The lifting surfaces were modeled by vortex panels with constant dipole and source strength. The wake panels were shed from the trailing edge according to the Kutta conditon and were converted to viscous vortex particles as soon as they were convected downstream. The fast multipole method was utilized to accelerate the computation of the induced velocity and of its gradient. The proposed method was validated by comparison with the Caradonna-Tung rotor in hover and the AH-1G rotor in forward flight. Finally, the vortex particle code was coupled with MBDyn, and the HART II baseline model was simulated. The following conclusions can be drawn:

1.
When the vortex panel method was combined with the viscous vortex particle method, the predicted blade pressure, normal thrust, and wake geometry match well with experimental data and CFD results both in the hover and forward flight conditions. 2.
The classical particle-based fast multipole method was extended to a hybrid particlepanel computational domain by using a semi-analytical approach. The convergence of numerical tests show that, as the truncation order grows, a reasonable accuracy can be obtained; in addition, the numerical tests show that the presented algorithm has a computational complexity of O(N).

3.
When t vortex code was coupled with the computational structural dynamics code MBDyn, the results show that the VPM/MBDyn coupling approach can predict the unsteady airloads well and is capable of capturing the BVI phenomenon, a necessary requisite for a future, accurate evaluation of the rotor aeroacoustic response. Institutional Review Board Statement: Not applicable.