Abstract
This study presents a detailed numerical analysis of nonlinear aeroelastic behavior in a two degree of freedom (DOF) model, focusing on plunge and pitch motions and employing the continuation method (CM) with an adaptive step size control algorithm. The research incorporates free-play nonlinearity at the plunge hinge, a common structural nonlinearity in aeronautics that can induce detrimental limit cycle oscillations (LCOs) during flight. By examining three scenarios—linear response, unhindered plunge motion, and nonlinear stiffness behavior—the study assesses the effects of free play on flutter and LCO phenomena, including discontinuity-induced bifurcations like grazing bifurcation. Additionally, the study explores parameter variation for nonlinear flutter analysis, revealing the dynamics of grazing bifurcation and its impact on LCO behavior. The research also demonstrates the method’s superior accuracy in flutter speed estimation and mode-switching identification, despite higher computational demands. The findings underscore the diminishing influence of nonlinear free-play behavior on LCO amplitude, providing insights with significant implications for aeroelastic design and aircraft safety.
Keywords: adaptive step size; aeroelasticity; continuation method; flutter; free play; limit cycle oscillation
1. Introduction
The investigation of aeroelasticity, a critical aspect of military and commercial aircraft design, encompasses the examination of interactions among aerodynamic, structural elastic, and inertial forces [1]. Within this context, both flutter and limit cycle oscillation (LCO) emerge as critical considerations for ensuring flight safety, gathering escalating attention and concern from the scientific and practitioner communities [2–5]. The nominal flutter is characterized as an infinitesimally small oscillation poised at the threshold of dynamic instability [6], while LCO represents a postflutter phenomenon distinguished by a relatively substantial and consistent amplitude and frequency, governed by the nonlinear characteristics of the system [7]. LCO denotes a self-sustained, steady-state fluid-induced vibration, typically triggered by either gradually surpassing the critical flutter speed or by a sudden external disturbance, such as a gust force. The majority of LCOs can be categorized into two distinct types [8]: benign LCO and detrimental LCO, both diverging from linear oscillations due to the amplitude of oscillatory behavior being independent of initial conditions. The benign LCO, or “good” LCO as depicted in Figure 1(a), consistently exhibits stability or attraction. This signifies that if an aeroelastic system experiences a benign LCO, any sudden augmentation or reduction in oscillatory amplitudes tends to converge toward the benign LCO state, irrespective of the initial condition and forcing function. In contrast, the detrimental LCO, or “bad” LCO, as illustrated in Figure 1(b), invariably demonstrates instability or repulsion, indicating that trajectories in the proximity of the limit cycle will diverge from the detrimental LCO state.
(a)
(b)
The prevalence of flutter and LCO phenomena in actual aircraft flight dynamics is predominantly instigated by various structural nonlinearities. These include, but are not limited to, the presence of external stores, free play, hysteresis, and cubic stiffness, as substantiated in seminal works [9–12]. Prior investigations have established, through both numerical [13] and empirical methodologies [14], the coexistence of both detrimental and benign LCO manifestations at velocities below the conventional flutter threshold in aircraft exhibiting free play nonlinearity within a singular control surface. Free play nonlinearity, a significant contributor to structural irregularities, commonly originates from factors such as spring fatigue, loose linkage, or joint slippage. This phenomenon is characterized by a range of structural motion devoid of restorative forces or torque as depicted in Figure 2. Consequently, such conditions may precipitate LCOs through Hopf bifurcation mechanisms [15]. According to the directives outlined in Guidance Military Specification [16], aircraft featuring movable control surfaces are mandated to ensure that the free play gap angle remains within specified limits to mitigate the risk of LCOs and potential structural failure. However, the required limits in practice are usually too conservative to meet and thus increase the cost of manufacturing and additional costs associated with inspection [17]. Consequently, there is a notable interest in the development of an effective methodology for assessing the impact of control surface free play on LCO. Such an approach is of high interest during the preliminary stages of design, manufacturing, and inspection, offering potential cost savings and efficiency improvements [18].
The literature is replete with both numerical and experimental inquiries into the aeroelastic characteristics of systems exhibiting free-play nonlinearity. Focused attention has been given to flutter and LCO phenomena. Investigations have encompassed various forms of LCOs on rigid wing models with free play in pitch motion, as demonstrated in studies by Yang and Zhao [19]; Lee, Price, and Wong [20]; Liu, Wong, and Lee [21]; and Vasconcellos et al. [22]. Additionally, comprehensive studies have been conducted on a flexible, fully movable wing model with a free play gap in the pitch shaft, with corresponding flutter responses detailed in Tang and Dowell [23, 24]. The LCO responses of typical wing sections with varying degree of freedom (DOF)—including plunge, pitch, and control surface motion as explored by Block and Strganac [25]; Ko, Strganac, and Kurdila [26]; and Singh and Yim [27], or an expanded four DOFs incorporating an additional leading control surface motion as in Panta et al. [28]; Tang, Li, and Dowell [29]; and Mannarino and Dowell [30]—have also been thoroughly investigated. Complex dynamical behaviors of a conceptual two-dimensional airfoil with both structural and aerodynamic nonlinearities excited by different types of external loads are thoroughly investigated by Liu, Xu, and Li [31] and Guo et al. [32], and an overview of recent developments related to this field is provided by Liu et al. [33].
At its core, the study of both flutter and LCO phenomena is an exploration of fluid-structure interaction (FSI) problems. A direct approach to analyzing these phenomena involves coupling fluid dynamics (via Euler/Navier–Stokes equations) with structural modal equations within the time domain, a method effective in monitoring the onset of dynamic instabilities. This coupling, as applied in computational structure dynamics (CSD)–computational fluid dynamics (CFD) solvers, has yielded accurate time-domain solutions for LCO/flutter predictions, as evidenced in the work of Robinson, Batina, and Yang [34]. However, this time-domain approach is notably time-intensive, particularly in predicting LCO/flutter onset, due to the requirement for multiple iterative and ad hoc simulations to accurately delineate the flutter boundary within a given aerostructural system, as noted by Schuster, Liu, and Huttsell [35]. To circumvent this drawback, numerous unsteady aerodynamic reduced-order models (ROMs) have been developed. For instance, the harmonic motion method, a lower-order linear aerodynamic model, offers a more pragmatic and efficient estimation of LCO, as outlined by Greco et al. [36] and Thomas, Dowell, and Hall [37]. In recent studies, Liu et al. [38] have introduced a novel linear matrix inequality (LMI) approach using convex optimization techniques, that demonstrates the numerical solvability of all uncertainties encountered during maneuvering processes. This method is not only pertinent to spacecraft attitude control but also applicable to second-order system applications.
Therefore, the crux of addressing the flutter challenge lies in the development of an efficient, accurate, nonlinear aeroelastic analysis tool, balancing computational expediency with robust physical interpretation. In evaluating the impact of multiple structural nonlinearities on flutter and LCO phenomena, translating linearized unsteady aerodynamic responses into the frequency domain emerges as a more viable alternative. This approach allows for the decoupling of structural dynamic equations from the repeated execution of costly CFD solvers, facilitating a swifter analysis of aeroelastic systems. In frequency-domain aeroelastic analyses, it becomes imperative to track aeroelastic modes as functions of dynamic pressure (or airspeed) to assess critical growth rates. However, this task grows increasingly complex in aircraft models exhibiting single or multiple structural nonlinearities. Traditional eigenvalue analysis methods often struggle with closely clustered modes in the distribution of eigenvalues and eigenmodes, leading to heightened sensitivity and potential discontinuities or misidentification of modes at mode-switching points. Such numerical challenges have been observed in commercial software relying on classical eigenvalue methods [39].
In this investigation, we conduct a comprehensive numerical analysis of the nonlinear aeroelastic behavior in a two DOF aeroelastic model, specifically focusing on plunge and pitch motions. This analysis is facilitated by an innovative and efficient methodology, namely, the continuation method (CM) augmented with an adaptive step size control algorithm. In this model, the aerodynamic response to oscillatory frequencies is presumed linear for minimal deflections, a common assumption in aeroelastic studies. A critical aspect of our study is the incorporation of free-play nonlinearity at the plunge hinge—a typical structural nonlinearity observed in practical aeronautical contexts. This nonlinearity is known to potentially induce LCOs that can be detrimental during flight. Our research specifically delves into the influence of this nonlinearity on the aeroelastic response, particularly in terms of LCO manifestation. Our analysis stratifies the problem into three distinct scenarios: (i) a linear structural response in the absence of free play motion; (ii) an unhindered plunge motion with no translational stiffness, as the motion amplitude remains within the free play gap angle; and (iii) a nonlinear stiffness behavior in the plunge motion when the amplitude exceeds the gap angle. This nonlinear stiffness is attributed to the presence of free play, introducing a complex behavior in the structure. Moreover, the study extends to evaluate the impact of this free-play nonlinearity on flutter and LCO phenomena, encompassing an exploration of discontinuity-induced bifurcations, such as grazing bifurcation. These investigations are pivotal in enhancing our understanding of the complex interplay between structural nonlinearities and aeroelastic responses, with significant implications for the design and safety of aerial vehicles.
2. Two DOF Airfoil Sections With Structural Nonlinearity
2.1. Mathematical Description
Consider a two-dimensional typical section of an aircraft wing with two DOFs, that is, plunge motion and pitch motion as presented in Figure 3. Both plunge and pitch motions are measured relative to the elastic axis located at from the middle point of the chord. The direction of the airspeed is assumed to be the positive -direction. The quantity is the dimensionless length scaled by the half-chord length between the center of gravity and the elastic axis of the wing section. For a two DOF airfoil as depicted in Figure 3, the linearized governing equations can be expressed into the following Equation (1) of the matrix form under the linearization assumption: where where , , and are the mass matrix, damping, and stiffness matrix, respectively. is the total mass of the wing, is the static mass moment, is the mass moment of inertia, is the plunge stiffness, and is the rotational stiffness. and are the structural damping which is modelled by viscous damping proportional to the velocity of plunge and rotation, respectively. The aerodynamic forces consists of the aerodynamic lift force and the moment around the elastic axis are given in the following Equations (2) and (3) based on the Theodorsen theory: where is the Theodorsen coefficient which is the function of the reduced frequency . The structural nonlinearity considered in the present study is the typical free-play behavior. It is assumed to occur in the translational motion, that is, the plunge motion . It is allowed to move freely within a specific range where the restoring force drops to zero. The relationship between the restoring force and displacement can be expressed in the following Equation (4) as a piece-wise function. where is the gap range for the plunge motion. Due to the existence of the free-play phenomenon, the aeroelastic system will consequently experience a nonlinear behavior when the amplitude of the plunge motion exceeds . Equation (4) in the form of a piece-wise function cannot be used directly in the flutter analysis conducted in the frequency domain. An alternative nonlinear function in the continuous form is hence needed to approximate the restoring force based on the described function method as outlined in Danowsky, Thompson, and Kukreja [17]. Assume that the plunge motion follows a harmonic way presented in Equation (5): where is the oscillatory plunge amplitude. The restoring force in Equation (4) can be expended by the Fourier’s series in the following Equation (6): where
The higher order of harmonic components is considered to contribute little to the final resulting force; hence, the approximate takes the following form in Equation (7) by keeping only the first fundamental components:
It is necessary to calculate the equivalent stiffness, , as a continuous function of , which is the amplitude of the harmonic plunge motion. Let , , then is calculated as follows:
Now, rewrite the as follows in Equation (8) and the equivalent translational stiffness can be obtained:
Noted that the equivalent stiffness in Equation (9) is expressed as the continuous nonlinear function of the amplitude . Consequently, can deviate substantially from its nominal values with varying amplitudes and hence could affect the aeroelastic behavior significantly. The structural parameters chosen for the aeroelastic model by Tang, Dowell, and Virgin [40] are listed in Table 1.
3. Numerical Algorithm for Aeroelastic Analysis
Assuming the solution of Equation (1) is of the form and substituting Equation (10) into Equation (1) leads to where is the dynamic pressure and is the density of the air. are the aerodynamic influence coefficient matrices which are functions of reduced frequency at a certain Mach number. In the present study, it is a matrix that consists of lift coefficient and moment coefficient which can be calculated as follows:
In the context of aeroelastic systems, the state variables and are defined as and , respectively, where denotes the complex eigenvector. The parameter represents the growth rate, which serves as a crucial indicator of the stability of the aeroelastic system, while signifies the oscillation frequency. It is noteworthy that the critical state of the system, characterized by , marks the initiation of the flutter phenomenon, denoted as . Additionally, it is imperative to acknowledge that the stiffness matrix is inherently a nonlinear function of the plunge amplitude . For a comprehensive examination of the flutter characteristics, it becomes imperative to solve the nonlinear flutter equations defined in Equation (11). This equation encompasses nonlinear structural variables and necessitates the application of an appropriate numerical algorithm for resolution. In the current study, we employ the modified CM coupled with an adaptive step size control algorithm. The foundational principles underlying this approach are succinctly summarized as follows:
The system of nonlinear equations represented by necessitates a solution, whereas functions are either already known or relatively straightforward to solve. To facilitate the incremental tracking of solutions from to , a parameter is introduced and systematically increased from 0 to 1. This enables the continuous tracing of solutions along the continuation function in the above Equation (12), where the path spans from to . Depending on the different selection of the continuation parameter , the CMs can be divided into two categories, that is, natural parameter CM or artificial parameter CM. In the natural parameter CM, is selected as a system parameter that possesses a physical interpretation. Conversely, in the artificial parameter CM, is constructed artificially without any physical meanings. In engineering practices, the CM predominantly utilizes a natural parameter for , integrating both and within the governing equations. For instance, in our research, the airspeed is selected as the continuation parameter for the mode tracking task, where represents the governing Equation (1). In this context, the specific formulations of and become secondary, as the primary focus shifts to resolving the solutions of , that is, the governing equations of the airfoil model under varying airspeeds . Throughout the tracking process, solutions are acquired using a predictor–corrector scheme at each incremental step. During the prediction phase, the Jacobian matrix , which encapsulates the partial derivatives of the nonlinear function with respect to as presented in Equation (13), is computed. Additionally, the tangent vector is determined as part of the predictive computation shown in Equation (14). This sets the stage for subsequent correction phases in the iterative scheme.
The predicted solution, for example, at step , that is, is computed following Equations (15) and (16) from the corrected solution at the step :
In the corrector phase, Newton–Raphson’s method is employed to compute the corrected solution as shown in Equation (17), wherein incremental adjustments are made to improve the accuracy of the solution obtained in the prediction phase. This iterative process is reiterated until either the absolute tolerance , expressed as , is met or the maximum iteration number is reached. The parameter represents the step size applied in the iterative correction process, guiding the incremental adjustments toward the accurate determination of the solution along the tracking curve.
The careful selection of the step size along the tracking curve is paramount, as it significantly influences the convergence and stability of the algorithm. Inappropriately large step sizes may lead to convergence to erroneous solutions or divergence, compromising computational efficiency. Conversely, setting the step size to an unnecessarily small value may sacrifice computational efficiency without providing substantial benefits. To address this critical aspect, a local control algorithm is employed in this study for adaptive step size adjustment. This algorithm dynamically adjusts the step size computed by Equations (18) and (19) during the iterative process, optimizing its magnitude to ensure a balance between computational efficiency and solution accuracy. where is the closeness index to assess the proximity between two solution curves, both and are small positive constants, and is the initial step size. and are the growth rate of oscillation for the th and th aeroelastic mode, respectively. In the present study, the tracking results for aeroelastic modes (growth rate vs. airspeed ) are illustrated. When using smaller step sizes ( and 0.5 m/s), the continuation algorithm successfully tracks plunge and pitch modes, revealing the development of flutter phenomena. However, with a larger step size (), the algorithm fails to distinguish between the two aeroelastic modes correctly. It struggles to converge the solution of the plunge mode to the correct path, and at higher airspeeds (), the two aeroelastic modes, and , become indistinguishable. This numerical failure is attributed to the selection of an improper large step size (), emphasizing the importance of careful step size tuning in continuation algorithms to ensure accurate and reliable tracking of solution curves (Figure 4).
(a)
(b)
4. Results and Discussions
In this section, the study presents and discusses results to highlight the key features of the proposed methods. The numerical algorithm, detailed in Section 3, is applied to solve the linearized aeroelastic equations of the test model depicted in Figure 3. The aeroelastic responses obtained from this numerical solution are then presented and subjected to discussion. Furthermore, flutter analyses are conducted for both linear and nonlinear structural systems using the modes tracking approach. The outcomes of these analyses are presented and thoroughly discussed. This includes an examination of the aeroelastic responses and an exploration of the findings related to flutter phenomena in both linear and nonlinear structural systems. The intention is to provide a comprehensive understanding of the performance and characteristics of the proposed methods, drawing insights from the results obtained through the applied numerical algorithm.
The behavior of the equivalent stiffness, , is dictated by three scenarios, based on the relationship between the gap angle and the oscillation amplitude of the plunge motion: (a) When , the absence of free play motion results in a linear structure, implying a constant equivalent stiffness, that is, . (b) For and , the plunge oscillates freely within the gap range, leading to . (c) When and , the equivalent stiffness becomes a nonlinear function of the oscillation amplitude according to Equation (9). This nonlinearity induces a nonlinear behavior in the structure due to free play. Realistically, all the above three situations could occur. The study investigates various cases, as outlined in Table 2, to comprehensively analyze the behavior of the system under different conditions. This approach allows for a thorough exploration of the effects of free play on the equivalent stiffness and the resulting structural behavior.
4.1. Case 1: Flutter Analysis of Aeroelastic Model With Linear Structural Properties
The provided passage describes a flutter analysis of the aeroelastic model shown in Figure 3, assuming linear structural properties (case no. 1). This analysis is compared with the results of an airfoil having a trailing-edge control surface (3 DOF system addressed in Yu, Damodaran, and Khoo [3]). The nominal stiffness for the plunge motion is restored by setting the free-play gap () to eliminate nonlinear free-play behavior. The governing equations (Equation (11)) become independent of the LCO amplitude, that is, . The results are presented in Figure 5, showing the modes tracking (variation of the growth rate of all structural modes vs. airspeed) in Figure 5(a) and frequency tracking (variation of the frequencies of all structural modes vs. airspeed) in Figure 5(b). The adaptive step size in the CM is utilized. The plunge mode goes into flutter first at the estimated flutter speed , which is determined to be 23.9 m/s. Near the flutter speed, the oscillating frequencies of the two associated aeroelastic modes (arising from the plunge and pitch mode) show a phenomenon called “frequency coalescence,” indicating energy transfer between the two modes, which is commonly considered to be a potential cause of flutter. Comparisons are made with a 3 DOF aeroelastic model (indicated in blue), which exhibits a slightly lower flutter speed at 23.4 m/s. For verification, the flutter analysis results for the linear structural system (gapless) through modes tracking are also compared with published results in Conner et al. [14]. The comparison is tabulated in Table 3. This comprehensive analysis provides insights into the flutter characteristics of the aeroelastic model under consideration and establishes comparisons with existing models and published results for validation purposes.
(a)
(b)
Table 4 provides a comparison between the estimated flutter speed and the total number of iterations needed to complete the modes tracking analysis using the CM algorithm with adaptive step sizes against several fixed step sizes. The adaptive step size approach demonstrates either enhanced efficiency or increased accuracy compared to fixed step sizes. For instance, when comparing the adaptive step sizes of 0.5 m/s and of 0.1 m/s with the reference case using the smallest fixed step size of 0.05 m/s), the adaptive approach achieves the same estimated flutter speed of 23.9 m/s while reducing the total iterations required by more than 59%, from 504 to 85. This highlights the improved efficiency of the adaptive step size methodology.
Moreover, it is important to note that using a fixed step size of might lead to algorithm failure due to the misidentification of mode-switching phenomena. Figure 6 illustrates an example of erroneous results due to misidentified aeroelastic mode switching when using a fixed step size of of 2 m/s. In this case, the tracking of two aeroelastic modes arising from the pitch mode and the plunge mode is not properly executed, resulting in a single tracking path after an airspeed of 13 m/s. This comparison underscores the advantages of employing adaptive step sizes in the modes tracking analysis, emphasizing either enhanced efficiency or increased accuracy over fixed step sizes while mitigating the risk of algorithm failure due to mode misidentification.
In order to validate the findings, a comprehensive CFD simulation is carried out, wherein the airfoil is directly integrated with the CFD solver in the time domain. The primary objective of this simulation is to determine the flutter speed for comparison with the data presented in Table 3. The outcomes are illustrated in Figure 7, revealing a flutter speed close to m/s. This observation is attributed to the stabilization of the airfoil displacement at an airspeed of m/s, as depicted in Figure 7(a), and its divergence at an airspeed of m/s, as shown in Figure 7(b). The flutter speed obtained from the full-order CFD simulation is found to be in excellent agreement with the current result () derived from the aerodynamic ROM and the CM. However, it is important to note that the time domain approach is computationally demanding for predicting the onset of LCO and flutter phenomena, as it necessitates numerous exhaustive and ad hoc runs for a given coupled aerostructural system. In contrast, the proposed methods in this study enable the rapid determination of the flutter speed within a few minutes. It is worth highlighting that both cases, as depicted in Figure 7, required several hours per simulation to gather sufficient data for assessing the flutter speed, in contrast to the efficiency of the proposed methods.
(a)
(b)
In the context of nonlinear flutter analysis, a specific test case featuring is selected to facilitate a comparative assessment with the adaptive step size control algorithm. The CM detailed in Section 3 is employed for this purpose. The LCO amplitude which exceeds the gap angle by a factor of 1.5 is deliberately chosen to induce a pronounced nonlinear effect on the structure. The real part of for plunge and pitch modes, which serve as indicators of aeroelastic stability, is systematically monitored as the flight speed is incrementally increased as part of the continuation process. The corresponding equivalent stiffness is computed using Equation (9), yielding . Through the application of modes tracking analysis procedures, the findings are presented in Figure 8.
In Figure 8, the employment of an adaptive step size within the CM is evident. Specifically, when the modes crossing region, delineated by the red ellipse, is identified within the airspeed range of 5 m/s to 15 m/s, the step size is automatically adjusted to 0.1 m/s. This adaptive modification is implemented to preemptively mitigate the potential divergence of the solver or misidentification of aeroelastic modes. In contrast, conventional point solution methods may necessitate manual intervention by the analyst, introducing a potential source of confusion and difficulty, especially in scenarios where two aeroelastic modes exhibit close proximity without intersecting. Outside the modes crossing region, a pragmatic restoration of the step size to 0.5 m/s ensues, thereby optimizing computational efficiency. In instances related to the neutral points regions, denoting the manifestation of instabilities in the plunge mode at an airspeed of , the step size undergoes another reduction. However, in this case, the reduction aims to enhance the accuracy in estimating the critical (lowest) flutter speed. Additionally, Figure 8 offers a detailed examination of the growth rate versus airspeed in proximity to these regions through zoomed variations presented as inset figures. This detailed visualization contributes insights into the aeroelastic behavior in the vicinity of critical regions.
4.2. Case 2: Flutter Analysis of Aeroelastic Model With Complete Free Play
In this examined case, the rotational displacement, denoted as pitch , is assumed to remain linear. While the translational motion, represented by plunge , is considered to undergo a complete free-play motion, devoid of any rotational stiffness, that is, . This assumption is grounded in the oscillation amplitude being smaller than the gap angle (). The outcomes of the two-mode tracking analysis for this configuration, employing the CM with adaptive step sizes ( of 0.5 m/s and of 0.1 m/s) and a fixed step size ( of 1 m/s), are presented in Figures 9(a) and 9(b). The onset of the flutter phenomenon is identified at an airspeed of , attributed to the pitch mode. Furthermore, Figure 9(a) delineates modes switching phenomena occurring around an airspeed of 27 m/s. The adept handling of the modes switching phenomenon is evident as the adaptive reduction of the step size from 0.5 m/s to 0.1 m/s is implemented when approaching the proximity of the switching region, as depicted within the red box. Conversely, Figure 9(b) illustrates a typical failure within the numerical scheme of the CM when a fixed step size is employed for the same case. This failure arises from the incorrect identification of the two switched aeroelastic modes, particularly evident in the red box of Figure 9(b), at an airspeed of 30 m/s, attributable to the use of a fixed and excessively large step size. This discrepancy underscores the significance of adaptive step size strategies in capturing nuanced aeroelastic phenomena accurately.
(a)
(b)
4.3. Case 3: Parameter Variation in Nonlinear Flutter Analysis
Parameter variation serves as an alternative methodology for nonlinear flutter analysis, obviating the necessity for repeated modes tracking. This approach involves systematically varying the LCO amplitude within the CM algorithm while directly monitoring the flutter speed for diverse amplitudes of the displacements of structural modes. Specifically, the investigation focuses on maintaining linearity in the rotational displacement , while concurrently varying the normalized amplitude of the plunge motion, denoted as , across a spectrum ranging from 1.0 to 3.0. This exploration accommodates nonzero free-play gap angles, that is, . Given the airfoil’s oscillation at freestream velocities below the linear flutter speed, Figure 10(a) presents the normalized LCO amplitude variation with increasing airspeed. Notably, a continuous decrease in normalized amplitude is observed, transitioning from 3.0 to 1.3 as the airspeed varies from 27 m/s to 34 m/s. Subsequently, at , a subtle decrement in amplitude from 1.3 to 1.15 is observed. This abrupt change, termed “grazing bifurcation” following the nomenclature of Kowalczyk et al. [41], signifies a pivotal event in the LCO behavior. After the occurrence of the grazing bifurcation, the LCO amplitude undergoes a continual decrease in airspeed, extending from 34 m/s to 35.5 m/s. To gain deeper insights into the dynamics associated with grazing bifurcation, scrutiny of Figure 10(b) is imperative. This figure, derived from the parameter variation approach within the CM solver, delineates the stable part of the detrimental LCO curve, that is, the curve and the unstable part of the detrimental LCO curve, that is, the curve within the red dashed square. As articulated by Van Rooij [42], the unstable LCO behaves like a repeller, while the stable LCO functions as an attractor. Consequently, as the normalized LCO amplitude decreases with increasing airspeed, reaching point D, any incremental airspeed rise triggers a precipitous drop in amplitude from point D to point C, aligning with the observations in Figure 10(a). Similarly, if the LCO amplitude diminishes with airspeed along the unstable LCO curve to the saddle-node bifurcation point A, a sudden jump from A to B transpires. This dichotomy implies that if the oscillation amplitude remains below the unstable LCO amplitude, it tends to decay to zero. Conversely, if the postdisturbance oscillation amplitude surpasses that of the unstable LCO, a stable LCO ensues. As the LCO amplitude increases, the influence of nonlinear free-play behavior gradually diminishes, becoming negligible as the LCO speed approaches the threshold for the gapless flutter curve beyond a normalized plunge amplitude of of 24.
(a)
(b)
5. Conclusions
In this study, we elucidate the efficacy of an advanced adaptive step size control algorithm, developed within the framework of the CM, and its applicability to two typical aeroelastic analysis tasks: modes tracking and parameter variation. To benchmark its performance in these areas, particularly in modes tracking and modes switching, the results derived from the CM are cross-checked with extant numerical and experimental data. A notable outcome of this research is the enhanced accuracy in flutter speed estimation achieved by the CM with adaptive step size, in comparison to conventional eigenvalue analysis or the - method. This increased precision is attributed to a denser concentration of interpolation points in the vicinity of the flutter onset, resulting from the algorithm’s step size reduction. Compared to the traditional eigenvalue methods and the fixed step size CM, the proposed algorithm shows an obvious advantage in managing the challenge of mode switching. This is primarily due to its innovative step size control mechanism, which allows for adaptive adjustments of the step size in proximity to switching points. This adaptability ensures a delicate balance between accuracy, stability, and computational efficiency, effectively reducing the risk of numerical inaccuracies or misidentification of aeroelastic modes. It is, however, important to acknowledge that, relative to traditional methods, the algorithm proposed in this paper demands greater computational resources. Nonetheless, the trade-off is justified, particularly when considering two practical aspects: (1) In preliminary design phases where repeated tracking of multiple modes with varying system parameters is essential, the risk of algorithmic failure and subsequent mode misidentification can lead to substantial time consumption and possibly necessitate manual correction by researchers. (2) With the ongoing enhancement of computer processing capabilities and the adoption of parallel computing across multiple central processing units (CPUs), the gap in computational efficiency between classical methods and our proposed algorithm is expected to narrow. Additionally, this study observes a marginal increase in flutter speed (from 23.4 m/s to 23.9 m/s) in the current model compared to the aeroelastic responses of a 3 DOF airfoil section. However, under the influence of complete free-play behavior in plunge motion, the flutter speed escalates by a significant 52.7%. In the realm of nonlinear flutter analysis, a consistent reduction in flutter speed is noted (from 35.5 m/s to 27.2 m/s) as LCO amplitudes increase. This study highlights that as the LCO amplitude escalates, the impact of nonlinear free-play behavior progressively wanes, becoming negligible when the LCO speed converges with the threshold of the gapless flutter curve.
Nomenclature
| : | damping matrix |
| : | Jacobian matrix |
| : | stiffness matrix |
| : | mass matrix |
| : | state vector for airfoil movement |
| : | continuation method equations |
| : | aerodynamic influence coefficient matrices |
| : | tangent vector in the predictor of continuation algorithm |
| : | location of the airfoil rotation axis |
| : | dimensionless length between the center of gravity and the elastic axis |
| : | structural damping to the plunge motion |
| : | structural damping to the pitch motion |
| : | dynamic pressure |
| : | the restoring force to the plunge motion |
| : | spring stiffness for plunge motion (per span) |
| : | equivalent translational stiffness for plunge motion (per span) |
| : | rotation stiffness for pitch motion (per span) |
| : | the lift force |
| : | the moment around the elastic axis |
| : | first moment of inertia of airfoil about the elastic axis |
| : | second moment of inertial of airfoil about elastic axis |
| : | flight speed |
| : | flutter velocity |
| : | flutter frequency |
| : | airfoil chord length |
| : | plunge motion of airfoil |
| : | total mass of the wing model |
| : | air density |
| : | pitch motion of airfoil |
| : | continuation parameter |
| : | step size of continuation parameter in th iteration |
| : | gap angle for the plunge free-play motion |
| : | oscillating frequency |
Data Availability Statement
Some data, models, or codes generated or used during this study are available from the corresponding author by request.
Conflicts of Interest
The authors declare no conflicts of interest.
Funding
The research results presented in this paper are funded by the general research project under the Beijing Office for Education Sciences the “14th Five-Year” Planning (the development and application of critical algorithms for virtual simulation experiments in the course “Fundamentals of Mechanical Design” aimed at vocational education (project number: CDDB23201)).
Acknowledgments
The authors have nothing to report.