Passive Aeroelastic Control of a Near-Ground Airfoil with a Nonlinear Vibration Absorber

1. Introduction

Aeroelasticity, the interaction between aerodynamic loads and structural dynamics, presents a major challenge in the design of modern flying vehicles. This interaction can induce undesirable vibrations, potentially leading to structural fatigue or even failure. The severity of aeroelastic effects, influenced by the complex dynamics of fluid-structure coupling, is more pronounced in innovative flexible wings [1]. Such wings are increasingly employed in unmanned aerial vehicles, high-altitude long-endurance aircraft [2], and unconventional aircraft [3] designed for enhanced mission performance. When flying close to the ground or during critical phases like takeoff and landing, these aircraft are subjected to aerodynamic changes due to the ground effect. While the ground effect enhances aerodynamic lift and reduces drag [4], research has indicated that it can also adversely affect dynamic aeroelastic stability [5,6]. This imposes significant constraints on performance and operational capabilities in low-altitude environments, emphasizing the importance of aeroelastic mitigation strategies.

Various control strategies offer promising solutions for mitigating aeroelastic instability. Both active and passive control methods have been explored over the years. Active control, which incorporates sensors, actuators, and feedback systems, allows real-time adjustment of aerodynamic loads and vibration modes, making it effective for varying flight conditions [7,8]. Semi-active devices, on the other hand, adapt their stiffness or damping properties dynamically in response to vibrations, offering adaptability with lower energy requirements than active systems [9]. In contrast, passive control methods employ energy-absorbing or dissipating devices to mitigate vibration. Passive systems are energy-efficient, less prone to failure, and may be fully reliable over the aircraft’s lifespan, as well as being simpler and more cost-effective than active and semi-active systems [10]. A tuned mass damper (TMD) is a classical example of a passive control device, consisting of a small mass attached to the primary structure via a linear spring and damper. When TMD parameters are tuned correctly, it can significantly delay the onset of flutter [11]. However, because a TMD operates linearly, it is not efficient on post-flutter instability behavior, such as limit cycle oscillations (LCO). Instead, researchers have proposed the use of nonlinear absorbers, based for instance on frictional effects [12,13] or a cubic restoring force – then called ”nonlinear energy sink” (NES) [14,15,16,17,18], and have reported a reduction in post-flutter LCO amplitude. When a TMD and a NES are combined – referred to as a nonlinear tuned vibration absorber (NLTVA) – the resulting system can offer even more robust control through optimized design. The linear component of the NLTVA can substantially increase the onset of flutter, while the nonlinear components can be tuned to reduce the LCO amplitude and promote supercritical behavior. Both numerical [19,20,21], and experimental studies [22,23] have demonstrated the effectiveness of NLTVA in controlling the flutter and post-flutter behavior. Most studies have primarily focused on structural effects, often using low-fidelity aerodynamics models. However, for aircrafts operating outside standard flight conditions, where aerodynamic interactions may play a significant role, assessing the NLTVA’s effectiveness from the aerodynamic perspective is essential.

Both quasi-steady and Theodorsen-based theories are often employed as complementary aerodynamic models in flutter mitigation. Theodorsen’s theory [24] is regarded as more accurate than quasi-steady theory, as it accounts for flow unsteadiness, including the phase lag between airfoil motion and aerodynamic loads, which is crucial in flutter dynamics. This approach is generally sufficient for predicting flutter onset and performing aeroelastic analysis, particularly when structural nonlinearities are significant and nonlinearities in the flow are minimal. However, when considering the ground effect, a high-fidelity model such as computational fluid dynamics (CFD) is often preferred in aerodynamics studies [25,26]. In aeroelastic studies, CFD is yet rarely considered due to its high computational demands. Instead, the unsteady discrete vortex method (UDVM), based on a potential flow theory, will be used in this study. This method is better suited for modeling interference effects, such as the ground effect [27,28], while providing a computationally efficient approach.

The objective of this study is to numerically assess how an NLTVA can improve aeroelastic performance under ground effect conditions, addressing both the critical flutter speed and the resulting post-flutter behavior. The problem is analyzed using a typical section model supported elastically by nonlinear springs, which provides a practical approximation of a high-aspect-ratio wing and is widely used in flutter calculations [20,29,30]. The aeroelastic equations of motion are formulated using Lagrange’s equations. Aerodynamic loads are computed using the UDVM, while the ground effect is simulated using the mirror-image method. The transient aeroelastic behavior at, before, and after flutter onset is obtained by numerical time integration. The influence of key design parameters, including the linear stiffness and damping of the NLTVA and the airfoil’s height from the ground, on flutter dynamics is analyzed. Additionally, the effect of the NLTVA’s nonlinear component on LCO amplitude and bifurcation behavior is investigated. The modeling and numerical solution procedures are presented in Section 2, with results and analysis in Section 3.

2. Mathematical Modeling

2.1. Typical Airfoil Section with Nonlinear Absorber

The structural model of the wing used in the following developments is illustrated in Figure 1. It incorporates a typical airfoil section (TSM) in a uniform flow $\left(\mathrm{U}\right)$. The airfoil is supported by torsional and flexural springs at the elastic axis. These springs provide linear stiffness, denoted by $k_{\alpha }$ and $k_h$, and include damping, denoted by $c_{\alpha }$ and $c_h$, which influence airfoil pitch ($\alpha$) and heave (h) motions, respectively. Moreover, the structural model accounts for geometric nonlinearities. The impact of the nonlinear effects in the heave mode has been shown to be negligible [20], therefore, the nonlinearity considered is exclusively in the pitch mode, modeled by a cubic stiffness denoted $k_{{\alpha }_3}$. The nonlinear absorber, connected to the TSM at a distance l from the elastic axis, consists of mass ($m_z$) supported by linear stiffness ($k_z$), cubic stiffness ($k_{z_3}$), and damping ($c_z$). The cubic stiffness term in the absorber primarily affects the amplitude of the nonlinear response after flutter; without it, the absorber functions as a traditional tuned mass damper. The absorber relative transverse motion to the airfoil is denoted by z. The governing equations of the TSM coupled with the nonlinear absorber are written as [20]:

$$\begin{array}{c}{I}_{\alpha }\ddot{\alpha }-S_{\alpha }\ddot{h}+c_{\alpha }\dot{\alpha }+c_{z}l(l\dot{\alpha }+\dot{h}+\dot{z})+k_{\alpha }\alpha +k_{z}l(l\alpha +h+z)\hfill \\ \hfill +k_{{\alpha }_3}{\alpha }^3+k_{z_3}l{(l\alpha +h+z)}^3=M_{\alpha },\end{array}$$

(1a)

$$m\ddot{h}-S_{\alpha }\ddot{\alpha }+c_h\dot{h}+c_z(l\dot{\alpha }+\dot{h}+\dot{z})+k_{h}h+k_z(l\alpha +h+z)+k_{z_3}{(l\alpha +h+z)}^3=L,$$

(1b)

$$m_z\ddot{z}+c_z(l\dot{\alpha }+\dot{h}+\dot{z})+k_z(l\alpha +h+z)+k_{z_3}{(l\alpha +h+z)}^3=0,$$

(1c)

where m is the mass of the TSM, and $I_{\alpha }$ and $S_{\alpha }$ are its moment of inertia and its static moment of inertia about its elastic axis, respectively. L and $M_{\alpha }$ represent the aerodynamic lift and moment about the elastic axis, which will be detailed in the following section. Equation (1) can be rearranged into standard matrix form as:

$$\mathit{M}\ddot{\mathit{q}}+\mathit{C}\dot{\mathit{q}}+\mathit{K}\mathit{q}+{\mathit{K}}_{nl}\left(\mathit{q}\right)=\mathit{f}\left(t\right)$$

(2)

where $\mathit{M}$, $\mathit{C}$, and $\mathit{K}$ are the mass, damping, and stiffness matrices of the system. ${\mathit{K}}_{nl}$ represents the nonlinear restoring forces, including geometric nonlinearities and nonlinear effects from the absorber, $\mathit{f}\left(t\right)$ accounts for the aerodynamic loads, and $\mathit{q}={[\alpha ,h,z]}^T$ contains the generalized coordinates.

In the numerical simulations presented in Section 3, non-dimensional values will be used. The airfoil model is then characterized by the following parameters: the non-dimensional locations of the elastic axis $\left(a\right)$ and center of mass $\left(x_a\right)$; the uncoupled heave-to-pitch frequency $(\sigma ={\omega }_h/{\omega }_{\alpha })$; the radius of gyration $\left(r_{\alpha }\right)$; and the non-dimensional airfoil mass $\left(\kappa \right)$. The structural damping coefficients are denoted by ${\zeta }_{\alpha }$ and ${\zeta }_h$ for pitch and heave, while the cubic restoring forces coefficient for pitch is denoted by ${\xi }_{\alpha }$. Furthermore, the quantity ${\mathrm{U}}^{*}$ and $\rho$, referring to the reduced airspeed and air density, are used. For a given pitch frequency $\left({\omega }_{\alpha }\right)$ and semi-chord length $\left(b\right)$, the following relationships are applied in terms of physical quantities:

$$\begin{array}{cc}\hfill & m=\pi \rho b^2\kappa ,S_{\alpha }=mb{x}_{\alpha },I_{\alpha }=m{b}^2{r}_{\alpha }^2,\hfill \\ \hfill & k_{\alpha }=I_{\alpha }{\omega }_{\alpha }^2,k_h=m{\omega }_h^2,c_{\alpha }={\zeta }_{\alpha }m{b}^2{\omega }_{\alpha },c_h={\zeta }_{h}m{\omega }_{\alpha },\hfill \\ \hfill & k_{{\alpha }_3}={\xi }_{\alpha }m{b}^2{\omega }_{\alpha }^2,U^{*}=U/\left({\omega }_{\alpha }b\right).\hfill \end{array}$$

Concerning the NLTVA, the ratio between its mass and that of the airfoil $\left(ϵ\right)$ is set to $5\%$, following [20], with the NLTVA positioned at the leading edge of the airfoil. The other parameters that define the NLTVA are considered design parameters. These are the ratio of NLTVA-to-pitch frequency $(\gamma ={\omega }_z^2/{\omega }_{\alpha }^2)$, the damping coefficient of the NLTVA $\left({\zeta }_z\right)$, and the coefficient of the NLTVA’s nonlinear restoring force $\left({\xi }_z\right)$. The following relations are applied in terms of physical quantities:

$$\begin{array}{cc}\hfill & m_z=ϵm,k_z=m_z{\omega }_z^2,c_z={\zeta }_z{m}_z{\omega }_{\alpha },\hfill \\ \hfill & k_{z_3}={\xi }_z{m}_z{\omega }_{\alpha }^2/b^2.\hfill \end{array}$$

2.2. Ground Effect Aerodynamics

The aerodynamics of the airfoil in proximity to the ground are modeled using the UDVM, with the ground interference effect simulated using the mirror image method [31]. UDVM is a computational approach based on potential flow theory. It is inherently nonlinear, applies to two-dimensional thin bodies in a uniform flow $\left(\mathrm{U}\right)$, and accounts for flow interference. However, it does not incorporate flow separation or thickness effects. Figure 2 illustrates the discrete vortex representation of the airfoil, its wake, and their mirror images. Here, the airfoil is discretized into N panels, in which each panel is represented by a lumped vortex strength (noted $\mathsf{\Gamma }$) placed at the quarter chord of the panel. The mirror image method simulates the ground effect by considering the mirror image of the primary model, positioned symmetrically relative to the ground plane. Thus, the influence of the ground is modeled by the image vortices, which have the same strength but are opposite in direction to those of the primary model. The image vortices are located at a distance of $2d$ from the primary model, where d is the height of the structure above the ground. The numerical solution proceeds in a time-marching process formulated by satisfying the continuity equation complemented with boundary conditions. The aerodynamic model applied in this study is well-established; for a comprehensive mathematical foundation, readers are encouraged to consult Katz and Plotkin’s textbook [31]. The main constitutive equations are recalled next.

The basic condition of physical flow over an airfoil (i.e. the flow should not penetrate the surface) is satisfied by enforcing the zero normal flow boundary condition. This condition is applied at the control point of each panel to solve for the unknown bound and wake vortex strengths. The zero normal flow boundary condition near the ground is written as:

$$(\nabla {\mathsf{\Phi }}_B+\nabla {\mathsf{\Phi }}_{I_B}+\nabla {\mathsf{\Phi }}_W+\nabla {\mathsf{\Phi }}_{I_W}-{\mathit{V}}_0-{\mathit{v}}_{\mathrm{rel}}-\mathsf{\Omega }\times \mathit{r}).\mathit{n}=0$$

(3)

where $(\nabla {\mathsf{\Phi }}_B).\mathit{n}$ and $(\nabla {\mathsf{\Phi }}_{I_B}).\mathit{n}$ represent the self-induced component of the normal velocity due to the primary and image vortices, respectively. Likewise, $(\nabla {\mathsf{\Phi }}_W).\mathit{n}$ and $(\nabla {\mathsf{\Phi }}_{I_W}).\mathit{n}$ accounts for the induced normal velocity from the wake vortices of the primary and image airfoils. The remaining terms account for the kinematic contribution to the normal velocity from the airfoil’s motion, where ${\mathit{V}}_0$ is the instantaneous velocity of the airfoil origin, ${\mathit{v}}_{\mathrm{rel}}$ is the relative motion of the chord line within the body coordinate system, $\mathit{r}$ is the position vector, and $\mathsf{\Omega }$ is the instantaneous rotation of the body coordinate system.

The flow leaving the trailing edge is governed by the Kutta condition, ensuring smoothing vorticity shedding. In addition, Kelvin’s theorem is applied to determine the strength of the shed wake according to the conservation of total circulation in the flowfield, written as:

$$\mathsf{\Gamma }\left(t\right)-\mathsf{\Gamma }\left(t-\Delta t\right)+{\mathsf{\Gamma }}_{W_t}=0$$

(4)

where $\mathsf{\Gamma }\left(t\right)$ and $\mathsf{\Gamma }\left(t-\Delta t\right)$ are the sum of all vortex strengths on the airfoil at the current and previous time step, respectively. The most recent wake vortex, with strength ${\mathsf{\Gamma }}_{W_t}$, is positioned within the interval $[0.2,0.3]\mathrm{U}\Delta t$ behind the trailing edge to correct for the wake-discretization error [31]. At each time step, Equations (3) and (4) are solved simultaneously to determine the strength of the airfoil vortices and the latest wake vortex. Once these strengths are known, the pressure along the airfoil is computed from the unsteady Bernoulli’s equation. The pressure difference at the ith control panel is given by:

$$\Delta p_i={\displaystyle \frac{{\rho }_{\infty }}{2}}\left({\left(V_{rel}^2\right)}_i-{\left(V_i^{\tau }\right)}^2-2{\displaystyle \frac{{\sum }_{k=1}^i\left({\mathsf{\Gamma }}_k\left(t\right)-{\mathsf{\Gamma }}_k(t-\Delta t)\right)}{\Delta t}}\right)$$

(5)

where ${\rho }_{\infty }$ is the fluid density, $V_{rel}$ is the magnitude of the velocity $(-{\mathit{V}}_{\mathbf{0}}-{\mathit{v}}_{\mathit{rel}}-\mathsf{\Omega }\times \mathit{r})$ which is the kinematic velocity due to the motion of an airfoil, and $V^{\tau }$ is the tangential speed. The aerodynamic lift $\left(L\right)$ and moment $\left(M_{\alpha }\right)$ at the elastic axis are then calculated by integrating the pressure along the chord, and subsequently used in Equations (1)-(2).

The wake vortices shed from the airfoil are free to deform and move downstream with the local flow velocity. The circulation strength around each vortex remains unchanged after detaching from the trailing edge, according to Helmholtz’s theorem. As the vortices are convected downstream, their new position is determined by calculating the induced velocity components $(u_w,w_w)$, which are influenced by the singularities from the primary and image airfoils as well as their wake. Consequently, the wake vortex points are transported along the distance $\Delta x$ and $\Delta z$ in the plane, achieving a vortex wake rollup, as given by:

$$(\Delta x,\Delta z)=(u_w,w_w)\Delta t.$$

(6)

A numerical problem is often encountered while estimating the induced velocity between two closely spaced vortices, which often results in large values that are unrealistic compared to the real flow. To address this, the Lamb-Oseen vortex model in conjunction with vortex-core growth mode is employed [32,33]. Vortex-core growth model incorporates the effect of eddy viscosity and allows the core radius of a vortex to grow with time, which leads to more accurate vortex dynamics.

2.3. Aeroelastic Solution

The second-order differential equation (Equation (2)) is first transformed into a set of first-order equations and expressed in state-space form. The state-space aeroelastic equation is solved iteratively in the time domain using a fifth-order Runge-Kutta numerical integrator to obtain the dynamic response of the system. The integrator calculates the airfoil’s position, which is then provided to the aerodynamic model. In return, the aerodynamic model computes the aerodynamic loads and feeds them back to the integrator. To maintain stability and prevent the response from diverging, an integral time step is employed within the integrator, during which the aerodynamic load is assumed to remain unchanged. Figure 3 illustrates the flowchart for performing the time-domain aeroelastic analysis using the UDVM.

3. Results and Discussion

3.1. Model Validation

The present method is first validated against the results presented in [5] – using the same numerical values – for near and far from the ground cases, without the absorber. Far from the ground, the present method predicts flutter instability at a reduced speed of 2.091 as compared to 2.073 reported in [5]. The difference in the reduced speed is attributed to the placement of the latest wake vortex: in the present approach, the latest wake vortex is positioned at $0.2\mathrm{U}\Delta t$ behind the trailing edge, whereas in [5], it is positioned at the trailing edge. The aeroelastic behavior is characterized by examining the evolution of the transient response over time. In Figure 4 transient responses are shown at three different reduced speeds, each starting with an initial pitch angle of 5 degrees. At reduced flutter speed of $2.091$, the response is self-sustaining. Below this speed, the response dampens out, while above it, the response grows to a larger amplitude and stabilizes into a LCO. The development of LCO above the flutter speed, in the absence of structural nonlinearities, is attributed solely to aerodynamic nonlinearities, which are captured in this study using the UDVM.

As the airfoil moves closer to the ground, a significant change occurs in the wake dynamics, resulting in increased aerodynamic loads. Consequently, the closer the airfoil is to the ground, the stronger the ground effect becomes, destabilizing the airfoil at lower speeds. This is illustrated in Figure 5, showing a significant decrease in the flutter speed as the airfoil height above the ground is decreased. In the figure, $L{\mathrm{w}}_1$ relates to the case where the most recent wake vortex is placed at the trailing edge whereas it is placed at $0.2\mathrm{U}\Delta t$ behind the trailing edge for $L{\mathrm{w}}_2$. The latter case shows a slight increase in flutter speed. Nevertheless, all results are in good agreement with the references [5,6].

3.2. Linear Flutter Analysis

In this section, neither the geometric nonlinearities of the airfoil nor the nonlinear stiffness of the absorber are considered. The absorber is then similar to a simple TMD. The structure is treated as purely linear, though simulations with the UDVM inherently include aerodynamic nonlinearity effects. The results obtained from the UDVM are first compared with those derived using the Wagner approach (a linear, time-domain approximation of Theodorsen’s theory) as aerodynamic model [34].

The numerical values reported in [20], as recalled in Table 1, are used for the simulations. Without the absorber and far from the ground, the present approach predicts a reduced flutter speed of $1.350$, while the Wagner approach predicts a reduced flutter speed of $1.376$. For the same system, Malher et al. [20] reported the flutter speed at $0.934$. Their analysis utilized quasi-steady aerodynamics, which ignores the effects of flow unsteadiness, thus underestimating the flutter speed. Similar to the approach taken in [20], an optimization of the linear parameters of the NLTVA ($\gamma$ and ${\zeta }_z$) is next conducted, but incorporating the unsteady effects, in order to identify the optimal pair $({\gamma }_{opt},{\zeta }_{z_{opt}})$ that maximizes flutter speed. Using the Wagner method, eigenvalue analyses are performed over a range of $\gamma$ and ${\zeta }_z$. As shown in the results illustrated in Figure 6, the maximum reduced flutter speed of $1.762$ is achieved at ${\gamma }_{opt}=0.4523$ and ${\zeta }_{z_{opt}}=0.1492$.

Using the enhanced aerodynamic model of the UDVM, the flutter speed is calculated only as a function of $\gamma$, with simulations run at the optimal value of ${\zeta }_z$ due to the high computational demands of the model. Despite some differences attributed to the linear and nonlinear aerodynamic models, the results from both methods are consistent, as shown in Figure 7, with a maximum flutter speed of $1.734$ obtained using the UDVM, indicating that the addition of the NLTVA results in an approximate $29\%$ increase in flutter speed as compared to without absorber. It is also observed that increasing $\gamma$ has a stabilizing effect on the airfoil up to the optimal value, beyond which it destabilizes. Eigenvalue analyses for two cases of $\gamma$ using Wagner approach, one before and one after the optimal value, are shown in Figure 8. The results indicate that the unstable flutter mode for $\gamma =0.3$ is mainly associated with the pitch mode, while the instability for $\gamma =0.6$ is mainly linked to the heave mode. The sudden change in behavior illustrated in Figure 7 is then attributed to the transition between unstable flutter modes. The NLTVA absorbs energy from the pitch mode while dissipating energy to the heave mode, resulting in stabilization of the pitch mode but destabilization of the heave mode as $\gamma$ increases. The optimal value marks the point at which the switching of the flutter unstable mode occurs. In subsequent analyses, the optimal values of $\gamma$ and ${\zeta }_z$ are utilized.

For the aeroelastic system approaching the ground, the influence of ground proximity on flutter instability is computed as a function of the normalized height above the ground. The results, shown in Figure 9, compare systems with and without the absorber. In both cases, the ground destabilizes the airfoil and lowers the flutter speed as the height between the model and the ground is reduced. The inclusion of an absorber improved the flutter stability in the presence of the ground effect as compared to the system without it. Moreover, since the results in Figure 9 correspond to the optimal values determined for the system far from the ground, it is possible that different optimal values for $\gamma$ and ${\zeta }_z$ may exist under the influence of ground effects, which could further mitigate flutter instability.

3.3. Post-Flutter Analysis with a Nonlinear Structural Model

The post-flutter solution follows the process outlined in Figure 3. The coefficient of nonlinear restoring force for pitch $\left({\xi }_{\alpha }\right)$ is set to unity. In the numerical time integration simulation, the airfoil motion is initiated with an initial disturbance of pitch angle, while all other degrees of freedom are set to zero. The flow speed is initially set below the critical speed and is then gradually increased and decreased back to the initial value within a single computation. The resulting post-flutter responses, for both increasing and decreasing flow speed, are examined to determine the system’s criticality that describes the nature of the Hopf bifurcation related to flutter. The Hopf bifurcations of supercritical nature are characterized by the gradual change in the onset of LCOs whereas subcritical ones are characterized by jump phenomenon and the appearance of a bistable region, see [20] for instance.

First, the nonlinear analysis of an aeroelastic system without an absorber, and positioned far from the ground, is realized. The simulation is initiated with a pitch angle of 2 degrees, and reduced speed $\left(U^{*}\right)$ of $1.348$, which is below the reduced flutter speed of the system without NLTVA. The value of $U^{*}$ is gradually increased up to $1.84$ and subsequently decreased back to its initial value over time. Figure 10 shows the resulting pitch and heave responses. In the figure, the solid blue curve traces the maximum amplitude of the response during the increase in $U^{*}$, while the dashed red curve shows the maximum amplitude as $U^{*}$ decreases. Together, these curves effectively give insight on the LCO amplitude as a function of $U^{*}$. The post-flutter response would be supercritical if these two curves were to coincide. However, achieving this requires a much slower variation in $U^{*}$ in order to give the system enough time to reach a dynamic behavior close to a steady state at a given $U^{*}$ value. This problematic is even more pronounced when the disturbances are relatively small, leading to higher computation costs. Nonetheless, it is inquired that the post-flutter instability in the Figure 10 of the system without NLTVA is supercritical.

Similarly, the nonlinear analysis of the aeroelastic system with a NLTVA far from the ground is carried out with an initial pitch angle of 2 degrees and reduced speed $\left(U^{*}\right)$ of $1.71$ . The starting $U^{*}$ is below the critical speed of the system with NLTVA, but higher than that of the system without NLTVA. Figure 11 illustrates the pitch and heave responses as $U^{*}$ is gradually increased and then decreased for three different values of the nonlinear stiffness ${\xi }_z$. For ${\xi }_z=0$, which corresponds to a traditional TMD, the amplitude of the response rapidly increases to a large value as $U^{*}$ rises beyond the critical speed and exhibits a bistable region when $U^{*}$ decreases. This behavior is related to a subcritical bifurcation and it shows the negative impact of TMD on shifting the original supercritical bifurcation (as illustrated in Figure 10) to a subcritical one. In contrast, the introduction of a nonlinear absorber, as observed with ${\xi }_z=0.1$, promotes supercritical bifurcation leading to more controlled response amplitude. However, increasing ${\xi }_z$ may not guarantee a smooth behavior. This parameter significantly impacts the nonlinear response, potentially leading to the emergence of additional bifurcation points. As exhibited in [20], relatively high values of ${\xi }_z$ can result in the creation of limit points that correspond to distinct jumps: a transition to high-amplitude LCO as $U^{*}$ increases, and a shift to low-amplitude LCO as $U^{*}$ decreases. This phenomenon is depicted in Figure 11 with ${\xi }_z=0.22$.

Overall, the results show the NLTVA’s effectiveness in mitigating LCO amplitudes compared to TMD (${\xi }_z=0)$ and underscore the importance of optimally tuning ${\xi }_z$ to avoid subcritical instabilities and jumps, which pose a risk of structural failure. Similar conclusions were drawn in [20], though the numerical values differ, likely due to the use of a more simplified aerodynamic model.

In order to support the distinction between supercritical and subcritical bifurcations, the pitch response from Figure 11 is plotted on a logarithmic scale in Figure 12. Two cases are examined, ${\xi }_z=0$ and ${\xi }_z=0.1$. In Figure 12a, pitch envelopes are shown for increasing and decreasing $U^{*}$, while in Figure 12b, pitch envelopes are shown as a function of time for a constant speed of $U^{*}=1.737$ just beyond the critical speed. The envelope in Figure 12b initially follows the same path in both cases, but over time, they settle at different steady-state values. The time required to reach a steady state would increase with a smaller initial disturbance. The steady-state values in Figure 12b, marked as A and B, correspond to the markers A and B in Figure 12a at $U^{*}=1.737$. Marker A represents a pitch LCO amplitude of approximately 11 degrees with the TMD, while marker B corresponds to a pitch LCO amplitude of approximately 3 degrees reached with the NLTVA, clearly demonstrating the NLTVA’s effectiveness in suppressing oscillation. Furthermore, for the envelopes to pass through their respective markers A and B during both increasing and decreasing $U^{*}$, a slower variation in $U^{*}$ would be needed to allow the response to stabilize. This required variation in $U^{*}$ must be especially slow when the disturbance is small, as a larger initial disturbance drives the response to a steady state more rapidly. This behavior is observed for decreasing $U^{*}$, where the envelopes effectively pass through markers A and B. Furthermore, the primary change in slope in Figure 12a, observed with decreasing $U^{*}$ values, occurs approximately at the critical speed for ${\xi }_z=0.1$. However, this change is more pronounced and occurs at a lower critical speed for ${\xi }_z=0$. The latter case suggests a jump phenomenon, where the solution branch extends below the critical speed as $U^{*}$ decreases, further indicating subcritical instability.

A detailed bifurcation analysis was conducted in [12,20] through the use of continuation procedures to identify multiple solution branches as the flow speed U varies. In our current approach, the weakly coupled model between structural dynamics and aerodynamics, combined with time-delayed terms in the aerodynamic model, makes the application of continuation techniques particularly challenging and merits further consideration.

Considering the ground effect, Figure 13 shows pitch and heave responses of a system with a NLTVA in the post-flutter regime, analyzed for increasing and decreasing $U^{*}$. In this analysis, the airfoil is positioned closer to the ground at a height equal to the airfoil chord ($d=c$), and motion is initiated with 2 degrees pitch disturbance at a reduced speed of $1.656$, below the critical speed. Similarly to the simulations presented far from the ground, the responses are evaluated for three different values of ${\xi }_z$. When ${\xi }_z=0$, equivalent to a TMD, a subcritical bifurcation is observed even in the presence of the ground effect. Introducing a nonlinear component encourages supercritical bifurcation, as observed for ${\xi }_z=0.1$, with potential limit points bifurcation (and associated jumps) occurring for larger ${\xi }_z$, as observed for ${\xi }_z=0.22$. Comparing, with ground (Figure 13) and without ground effects (Figure 11) results, show no significant changes in the nature of the bifurcations and overall dynamics.

Figure 14 retrieves time histories and phase plots of the system near the ground at the reduced speed of $1.7$ and presents comparisons between various cases: without the absorber, with a TMD, and with a NLTVA configured at ${\xi }_z=0.1$ or at ${\xi }_z=0.22$. The reduced speed chosen voluntarily exceeds the critical speed for all systems. Simulations are initiated with a small disturbance of pitch angle $(0.01$ deg.) and are allowed to have prolonged transient for the system to develop LCO. The phase in Figure 14 is plotted for part of the transient response after the response has reached a steady state. Due to the ground influence, the equilibrium position of the responses is no longer symmetric but has shifted upward. Among the cases, the system without an absorber attains LCO most quickly and exhibits the largest amplitude. In contrast, the systems with absorber exhibit significantly reduced LCO amplitudes, as the absorber particularly affects the pitch, which is considered as most energetic degree of freedom of the system. Furthermore, the mitigation of LCO is especially effective when the NLTVA is employed.

To investigate the influence of the NLTVA on the LCO amplitude under ground influence, numerical simulations are performed at a constant speed of $U^{*}=1.75$. This speed is above the critical speed for cases far from the ground, thus, ensuring the development of LCO across all ground heights. In a single computation, the airfoil height relative to the ground is gradually varied, starting from twice the chord length and 2 degrees pitch disturbance, reducing to $40\%$ of the chord length, and then returning to the initial height. Figure 15 illustrates the pitch and heave responses for three different values of ${\xi }_z$. The results clearly demonstrate the effectiveness of the NLTVA in suppressing the amplitude of post-flutter amplitude even in the presence of strong ground effects. Notably, with the NLTVA $({\xi }_z=0.1)$, a reduction in amplitude approximately from $220\%$ to $40\%$ is observed as the airfoil height is varied from twice the chord length to the chord length compared to the TMD case (${\xi }_z=0$). However, as the ground height $\left(d\right)$ decreases, the critical flutter speed also drops (as it was shown in Figure 5 in the linear case), and thus the global effectiveness of the NLTVA diminishes at lower ground heights. For ${\xi }_z=0.22$, the results in Figure 15 indicate the coexistence of several stable branches within the height range of approximately $1.5$ to 2 times the chord length. A relatively small initial condition in pitch allows the system to follow a low-amplitude branch while decreasing the height. Conversely, if the initial condition were larger, the system could transition to a different basin of attraction, reaching the red stable branch -— this behavior is observed when increasing the height.

4. Conclusions

The study investigated the aeroelastic behavior of an airfoil equipped with a nonlinear passive absorber in uniform flow near the ground. This absorber acted as both a tuned mass damper and a nonlinear energy sink, effectively dissipating energy from the unstable system. The airfoil was modeled as a typical section with geometric nonlinearity, subjected to nonlinear aerodynamic loads using the unsteady discrete vortex method, while the ground effect was simulated via the mirror image method. The primary goal was to assess the absorber’s effectiveness in mitigating aeroelastic instability, particularly by delaying flutter onset and reducing limit cycle oscillation amplitude. Key design parameters of the nonlinear tuned vibration absorber (NLTVA) were analyzed, including linear stiffness, damping, cubic stiffness, and airfoil height from the ground. The approach also provided insights into modeling instability behaviors, such as supercritical and subcritical bifurcations, influenced by various factors including geometric nonlinearities and wake interactions. However, the use of continuation techniques to a nonlinear aeroelastic model with time-domain nonlinear aerodynamics is particularly challenging and merits further consideration.

The nonlinear tuned vibration absorber was designed to increase flutter onset speed. A linear stability analysis showed a $28.67\%$ increase in flutter speed with an optimal combination of linear stiffness and damping. Both the UDVM and Wagner’s approach demonstrated consistent results in the absence of ground effects. However, introducing ground effects reduced the flutter speed as the airfoil approached the ground, due to increased aerodynamic loads from altered wake dynamics. Importantly, the presence of the absorber significantly enhanced the flutter onset speed, improving it by $3.8\%$ to $28.6\%$ as the airfoil’s height above the ground varied from $30\%$ to $200\%$ of the chord length, compared to cases without the absorber.

Post-flutter analyses revealed that the nonlinear tuned vibration absorber can induce subcritical or supercritical behavior, and also jumps associated with limit points. Using only the linear component of the absorber resulted in subcritical bifurcation, while incorporating the nonlinear component shifted the behavior to supercritical bifurcation, significantly reducing the LCO amplitude under the ground influence. The global effectiveness of the NLTVA in reducing LCO amplitudes decreases at lower ground heights due to stronger ground influence. The system’s bifurcation behavior was sensitive to the nonlinear stiffness of the NLTVA, with higher stiffness values leading to jumps to higher amplitude LCO. To avoid subcritical bifurcation, which can cause sudden and explosive motion, it is crucial to determine an optimal nonlinear parameter for the NLTVA that works effectively in both near-ground and far-from-ground scenarios. This study underscores the effectiveness of NLTVA in mitigating the aeroelastic instabilities in a near-ground scenario, which could be advantageous in the dynamic stability control of aero-elastically vulnerable aircraft with highly flexible wings.