Nonlinear Aeroelastic Characteristics of a Supersonic All-Movable Fin with Single or Multiple Freeplay Nonlinearities

Abstract

Establishing the aeroelastic characteristic of all-movable fins with freeplay nonlinearities is one of the most common problems in the design of supersonic flight vehicles. In this context, this study provided novel points of view on the nonlinear aeroelastic characteristics of an all-movable fin with freeplay nonlinearities in its root. The unsteady aerodynamic model that was employed uses the second-order piston theory considering thickness effects. For a system with multiple freeplay nonlinearities, a discrete scanning method based on the describing function method was established to solve the limit-cycle oscillations (LCOs) and avoid the loss of solutions. Combining this with the time-domain integration method, the influences of the support stiffness at the root of the fin and the freeplay size ratio of the bending and torsional degrees of freedom on the dynamical response of the system were analyzed. The results demonstrate that systems with a single freeplay nonlinearity exhibit two completely different types of LCO, while systems with multiple freeplay nonlinearities exhibit complex dynamical behaviors such as LCO and quasi-periodic and chaotic motions. The path of a quasi-periodic torus breaking into chaos was observed. Furthermore, a harmonic initial condition for the time-domain integration is proposed; this can be used for a quick check of the frequency-domain calculation results.

1. Introduction

Control surfaces play a critical role in the flight of aircraft. With an all-movable fin, due to the looseness and excess tolerances in bolts, bearings, hinges, and actuators, freeplay is unavoidable. This will result in various nonlinear aeroelastic responses, such as limit-cycle oscillations (LCOs), a phenomenon involving structure-sustained limited-amplitude vibrations [1]. LCOs can occur below the flutter boundary, and these are referred to as “detrimental” LCOs; conversely, “benign” LCOs occur only above the flutter boundary [2]. Fatigue problems may occur in aircraft structures if the amplitude of LCOs is relatively small, while large amplitudes can directly lead to instantaneous structural failure [3]. Furthermore, there is a need for more frequent lubrication and maintenance checks, especially for aging aircraft, to ensure the freeplay magnitudes comply with relevant requirements, and this may be a notable financial burden [1]. Therefore, exploring the effect of control surface freeplay on aeroelastic characteristics has become an important subject.

The earliest work performed to understand the effects of concentrated structural nonlinearities on aeroelastic stability was initiated in the 1950s. Through analog computer and wind-tunnel tests, Woolston et al. [4] studied the effects of freeplay, hysteresis, and cubic spring stiffness on the flutter of a two-degree-of-freedom system; they noted that the flutter boundary may decrease due to the existence of structural nonlinearities. Later, Shen [5] applied the first-order harmonic balance method, also known as the describing function (DF) method, to nonlinear aeroelasticity problems, and the results using Woolston’s model were found to be consistent with the computer simulation results. In 1988, Yang and Zhao [6] carried out a wind tunnel test for a typical free-play airfoil section. The related results showed that the DF method can explain and predict most LCO behaviors with sufficient accuracy. With further research, a more comprehensive understanding of nonlinear aeroelasticity problems caused by structural nonlinearities has been gleaned. The jump and hysteresis phenomena [7], bifurcation, and chaotic behaviors [8] have been observed in the research. Tang and Dowell [9] considered a three-degree-of-freedom aeroelastic model with freeplay. They discussed the LCO behaviors and the sensitivity of nonlinear systems to initial conditions. In addition, the influence of the modal frequency ratio on LCO was studied. It was found that when the flap and pitch modal frequency ratio is lower than a certain value, the stable LCO branch will disappear, and only an unstable LCO branch will remain.

In recent years, some studies have been carried out from a more practical perspective. Padmanabhan and Dowell [10] considered the influence of gust load on LCO characteristics. Wayhs-Lopes et al. [11] discussed the influence of asymmetric freeplay and friction on the LCO characteristics of a typical airfoil section with a control surface. Iannelli et al. [12] considered the influence of model uncertainty on the Hopf bifurcation. Moreover, due to the sensitivity of nonlinear systems to the initial conditions, it is often necessary to explore a large number of initial conditions when performing time-domain simulations so as to completely analyze the nonlinear characteristics of the system. Wang et al. [13] proposed a state-space iterating scheme that can consider large numbers of combinations of initial conditions and effectively analyze the LCO stability with respect to initial conditions. Most of these studies took a typical airfoil section as the research object, and they established some classic analysis approaches for nonlinear aeroelastic problems: the frequency-domain approaches represented by the DF method [5] and the time-domain integration approaches represented by a combined scheme of the Hénon method and the fourth–fifth-order Runge–Kutta method (Hénon–RK45) [11,13,14,15,16]. The DF method has the common advantage of frequency-domain approaches, that is, high efficiency, while numerical integration methods have high accuracy but low efficiency. In 1999, Lee et al. [17] comprehensively summarized the previous research results on the nonlinear aeroelastic system of a typical airfoil section. In recent decades, Panchal et al. [1] and Bueno et al. [18] have also reviewed the research status of nonlinear problems of control surface structures.

For research using three-dimensional (3D) models, there are two main kinds of objects: wings with stores [19,20,21,22,23] and all-movable control surfaces [24,25,26,27]. A primary problem in such research is to find ways to establish reduced-order models of these systems. The Lagrange equation can be directly used for simple objects. For complex objects such as folding fins or folding wings, there are two main modeling methods: the modal synthesis method [26,27] and the fictitious mass method [23,28]. The main analysis approaches are also the DF method and the Hénon–RK45 numerical integration method. In addition to considering improvements and innovations in modeling and analysis methods, it is necessary to consider how research results can be applied to aircraft design. Padmanabhan et al. [20,21] found that friction nonlinearity will benefit a system in which there is no LCO occurrence below the linear flutter boundary; this may provide a solution for the suppression of LCOs. For an all-movable fin with freeplay nonlinearities, finding ways in which to delay or suppress LCOs through the reasonable design of structural parameters is a problem worth examining. So far, however, there has been little discussion of this problem.

Regarding multiple nonlinearities, Firouz-Abadi et al. [29] considered the flutter of a 3D rigid fin with freeplay in the flapping, plunging, and pitching degrees of freedom (DOFs). The results of the time-domain analysis showed that the freeplay nonlinearities lead to lower flutter-speed predictions compared to the linear system results. In addition to the time-domain integration approach, a frequency-domain iterative procedure based on the DF method has also been established and developed for systems with multiple structural nonlinearities. Laurenson and Trn [30] considered a missile control surface with roll and pitch freeplay nonlinearities with the DF method. However, in this study, the coupled responses between nonlinearities were not deeply considered; rather, the amplitude of one nonlinear element was fixed while the other was changed. In 1980, Breitbach [31] proposed an iterative procedure to achieve alignment of the oscillatory amplitude of each nonlinear element. Later, Lee [32] proposed a comprehensive frequency-domain iteration procedure for flutter analysis of large dynamical systems with multiple structural nonlinearities; the results of this frequency-domain approach were found to be consistent with those of the time-domain approach. Seo et al. [33] simplified the above solution process and emphasized the importance of the frequency ratio when dealing with a two-dimensional control-fin model. For all-movable fins, the current iterative procedure may lead to the loss of solutions. Furthermore, the iterative procedure may not converge when the selected iteration parameters are not appropriate.

In this study, we provided novel points of view on the nonlinear aeroelastic characteristics of a small aspect ratio all-movable fin with freeplay nonlinearities in its root. In the remainder of the paper, the basic modeling and analysis methods are introduced. Then, for a single freeplay nonlinearity, two different LCO characteristics are introduced, and a harmonic initial condition is proposed for time-domain integration to verify the frequency-domain LCO results. For multiple freeplay nonlinearities, a discrete scanning method based on the DF method is established to avoid the loss of solutions. Finally, the frequency- and time-domain approaches are respectively used to analyze the effects of the support stiffness and freeplay size ratio on the LCO characteristics.

2. Aeroelastic Model and Analysis Method

2.1. Structural Model

Figure 1 shows a diagram of the typical all-movable fin in a supersonic flow that is analyzed in this paper. A right-handed reference coordinate system is coincident with the fin root. The x-axis is along the free-stream direction, and the y-axis is along the fin span. The flexibility of the fin root support is modeled by two springs that constrain its bending (B) and torsional (T) motions, and $K_{\mathrm{B}}$ and $K_{\mathrm{T}}$ denote their linear spring stiffnesses.

Using the Branch Mode method [34], the linear fin can be divided into the following two branches: (1) a rigid fin with elastic supports at the root; (2) an elastic fin with rigid supports at the root. The bending and torsional modes of branch 1 and some elastic modes of branch 2 are retained for analysis. The displacement response in the z direction of the fin at position $(x,y)$ can be written as follows:

$$z(x,y,t)={\Phi }_{\mathrm{B}}\left(x,y\right)q_{\mathrm{B}}\left(t\right)+{\Phi }_{\mathrm{T}}\left(x,y\right)q_{\mathrm{T}}\left(t\right)+\sum _{j=1}^n{\Phi }_{\mathrm{E}j}\left(x,y\right)q_{\mathrm{E}j}\left(t\right)$$

(1)

where ${\Phi }_{\mathrm{B}}$ and ${\Phi }_{\mathrm{T}}$ represent the mode shape functions of the bending and torsional modes, ${\Phi }_{\mathrm{E}j}$ represents the jth elastic mode shape function, $q_{\mathrm{B}}$, $q_{\mathrm{T}}$ and $q_{\mathrm{E}j}$ are the corresponding generalized coordinates, n is the number of modes retained in branch 2. By means of the Lagrange equations, the motion equations of the linear fin can be derived as

$$\left[\begin{array}{cc}\begin{array}{cc}{I}_{xx}& -I_{xy}\\ -I_{xy}& I_{yy}\end{array}& \mathbf{RE}\\ {\mathbf{RE}}^{\mathrm{T}}& {\mathbf{M}}_{\mathrm{E}}\end{array}\right]\left[\begin{array}{c}{\ddot{q}}_{\mathrm{B}}\\ {\ddot{q}}_{\mathrm{T}}\\ {\ddot{\mathbf{q}}}_{\mathrm{E}}\end{array}\right]+\left[\begin{array}{ccc}{K}_{\mathrm{B}}& & \\ & K_{\mathrm{T}}& \\ & & {\mathbf{K}}_{\mathrm{E}}\end{array}\right]\left[\begin{array}{c}{q}_{\mathrm{B}}\\ q_{\mathrm{T}}\\ {\mathbf{q}}_{\mathrm{E}}\end{array}\right]=\left[\begin{array}{c}{Q}_{\mathrm{B}}\\ Q_{\mathrm{T}}\\ {\mathbf{Q}}_{\mathrm{E}}\end{array}\right]$$

(2)

where $I_{xx}$, $I_{yy}$, and $I_{xy}$ denote the mass moments of inertia about the corresponding axes in the reference coordinates, $\mathbf{RE}$ is the coupling mass term of branches 1 and 2, ${\mathbf{M}}_{\mathrm{E}}$ and ${\mathbf{K}}_{\mathrm{E}}$ are the generalized mass and stiffness matrix of branch 2, ${\mathbf{q}}_{\mathrm{E}}$ is the corresponding generalized modal coordinate vector, $Q_{\mathrm{B}}$ and $Q_{\mathrm{T}}$ are the generalized aerodynamic forces, ${\mathbf{Q}}_{\mathrm{E}}$ is the generalized aerodynamic force vector.

Generally, for a small aspect ratio supersonic fin, the structure of the fin is more rigid than the attachment part [35], which means the modal frequency of branch 2 is much higher than that of branch 1. So the consideration of branch 2 has little influence on the aeroelastic characteristics of the fin [36], and it can be approximated as a rigid body. Equation (2) can be reduced to the following form:

$$\left[\begin{array}{cc}{I}_{xx}& -I_{xy}\\ -I_{xy}& I_{yy}\end{array}\right]\left[\begin{array}{c}{\ddot{q}}_{\mathrm{B}}\\ {\ddot{q}}_{\mathrm{T}}\end{array}\right]+\left[\begin{array}{cc}{K}_{\mathrm{B}}& \\ & K_{\mathrm{T}}\end{array}\right]\left[\begin{array}{c}{q}_{\mathrm{B}}\\ q_{\mathrm{T}}\end{array}\right]=\left[\begin{array}{c}{Q}_{\mathrm{B}}\\ Q_{\mathrm{T}}\end{array}\right]$$

(3)

For an all-movable fin with symmetric freeplay in the bending and torsional directions, the nonlinear elastic restoring moments $M_{\mathrm{B}}$ and $M_{\mathrm{T}}$ can be described by

$$\begin{array}{c}{M}_{\mathrm{B}}\left(q_{\mathrm{B}}\right)=\left\{\begin{array}{cc}{K}_{\mathrm{B}}\left(q_{\mathrm{B}}-{\delta }_{\mathrm{B}}\right)\hfill & q_{\mathrm{B}}>{\delta }_{\mathrm{B}}\hfill \\ 0\hfill & -{\delta }_{\mathrm{B}}\le q_{\mathrm{B}}\le {\delta }_{\mathrm{B}}\hfill \\ K_{\mathrm{B}}\left(q_{\mathrm{B}}+{\delta }_{\mathrm{B}}\right)\hfill & q_{\mathrm{B}}<-{\delta }_{\mathrm{B}}\hfill \end{array}\right.\\ M_{\mathrm{T}}\left(q_{\mathrm{T}}\right)=\left\{\begin{array}{cc}{K}_{\mathrm{T}}\left(q_{\mathrm{T}}-{\delta }_{\mathrm{T}}\right)\hfill & q_{\mathrm{T}}>{\delta }_{\mathrm{T}}\hfill \\ 0\hfill & -{\delta }_{\mathrm{T}}\le q_{\mathrm{T}}\le {\delta }_{\mathrm{T}}\hfill \\ K_{\mathrm{T}}\left(q_{\mathrm{T}}+{\delta }_{\mathrm{T}}\right)\hfill & q_{\mathrm{T}}<-{\delta }_{\mathrm{T}}\hfill \end{array}\right.\end{array}$$

(4)

where ${\delta }_{\mathrm{B}}$ and ${\delta }_{\mathrm{T}}$ represent the unilateral sizes of the freeplay in the bending and torsional directions, respectively. Thus, the motion equations of the nonlinear fin can be written in the following form.

$$\left[\begin{array}{cc}{I}_{xx}& -I_{xy}\\ -I_{xy}& I_{yy}\end{array}\right]\left[\begin{array}{c}{\ddot{q}}_{\mathrm{B}}\\ {\ddot{q}}_{\mathrm{T}}\end{array}\right]+\left[\begin{array}{c}{M}_{\mathrm{B}}\\ M_{\mathrm{T}}\end{array}\right]=\left[\begin{array}{c}{Q}_{\mathrm{B}}\\ Q_{\mathrm{T}}\end{array}\right]$$

(5)

2.2. Aerodynamic Model

Hayes’ [37] and Lighthill’s [38] piston theory, a simple and practical method, has been applied to many supersonic aeroelastic problems [39]. In this study, Van Dyke’s second-order piston theory [40] considering thickness effects was selected as the aerodynamic model. The pressure coefficients can be expressed as

$$C_p=\frac{2}{M^2}\left[\frac{M}{\beta }\left(\frac{w}{a}\right)+\frac{G}{2\beta }{\left(\frac{w}{a}\right)}^2\right]$$

(6)

where

$$\beta =\sqrt{M^2-1},G=\frac{M^4\left(\gamma +1\right)-4{\beta }^2}{2{\beta }^3}$$

M is the Mach number, a is the local speed of sound, w represents the piston upwash, and $\gamma$ is the ratio of specific heat. The total pressure differential $\Delta p$ can be expressed as

$$\Delta p\left(x,y,t\right)=-\left(\frac{1}{2}\rho V^2\right)\left(\frac{4}{\beta }\right)\sum _{j=1}^n\left[1+G\frac{\partial H\left(x,y\right)}{\partial x}\right]\left[\frac{\partial {\Phi }_j\left(x,y\right)}{\partial x}{q}_j\left(t\right)+\frac{1}{V}{\Phi }_j\left(x,y\right){\dot{q}}_j\left(t\right)\right]$$

(7)

where ${\Phi }_j\left(x,y\right)$ is the jth modal shape, $q_j\left(t\right)$ is the jth generalized modal coordinate, n is the number of modes, $H\left(x,y\right)$ is the thickness function measured form the middle profile of the fin, $\rho$ is the air density, and V is the velocity of the free stream. The ith generalized aerodynamic forces can be expressed as

$$Q_i=\underset{S}{\iint }\Delta p\left(x,y,t\right){\Phi }_i\left(x,y\right)\mathrm{d}x\mathrm{d}y$$

(8)

Substituting Equation (7) into Equation (8) yields

$$Q_i=-(\frac{1}{2}\rho V^2)\left(\frac{4}{\beta }\right){\sum }_{j=1}^n\left[\left(E_{ij}+G{C}_{ij}\right)q_j+\frac{1}{V}\left(B_{ij}+G{D}_{ij}\right){\dot{q}}_j\right]$$

(9)

where

$$E_{ij}=\underset{S}{\iint }\frac{\partial {\Phi }_j\left(x,y\right)}{\partial x}{\Phi }_i\left(x,y\right)\mathrm{d}x\mathrm{d}y$$

$$B_{ij}=\underset{S}{\iint }{\Phi }_j\left(x,y\right){\Phi }_i\left(x,y\right)\mathrm{d}x\mathrm{d}y$$

$$C_{ij}=\underset{S}{\iint }\frac{\partial H\left(x,y\right)}{\partial x}\frac{\partial {\Phi }_j\left(x,y\right)}{\partial x}{\Phi }_i\left(x,y\right)\mathrm{d}x\mathrm{d}y$$

$$D_{ij}=\underset{S}{\iint }\frac{\partial H\left(x,y\right)}{\partial x}{\Phi }_j\left(x,y\right){\Phi }_i\left(x,y\right)\mathrm{d}x\mathrm{d}y$$

Equation (9) can be written in matrix form:

$$\mathbf{Q}=-(\frac{1}{2}\rho V^2)\left(\frac{4}{\beta }\right)\left[\left(\mathbf{E}+G\mathbf{C}\right)\mathbf{q}+\frac{1}{V}\left(\mathbf{B}+G\mathbf{D}\right)\dot{\mathbf{q}}\right]$$

(10)

Thus, the generalized aerodynamic forces in the time domain can be expressed as:

$$\mathbf{Q}=\frac{1}{2}\rho V^2{\mathbf{A}}_0\mathbf{q}+\frac{1}{2}\rho V{\mathbf{A}}_1\dot{\mathbf{q}}$$

(11)

where

$${\mathbf{A}}_0=-\left(\frac{4}{\beta }\right)\left(\mathbf{E}+G\mathbf{C}\right)$$

$${\mathbf{A}}_1=-\left(\frac{4}{\beta }\right)\left(\mathbf{B}+G\mathbf{D}\right)$$

2.3. Analysis Methods for Aeroelastic Model with Freeplay

2.3.1. Describing Function Method

Common linear flutter solution methods, such as the V-g method, p-k method, and root-locus method, have high computational efficiency. However, when freeplay nonlinearity is involved, these methods are no longer applicable. In this instance, an equivalent linearization method, the DF method, has solved this problem. Taking into account only the first-order harmonic component of the output of the nonlinear element under the harmonic input, the relationship between the equivalent stiffness and the LCO amplitude is established as follows [5]:

$$\frac{\overline{K}}{K}=\left\{\begin{array}{cc}0& A\le \delta \\ 1-\frac{2t_1}{\pi }-\frac{sin2t_1}{\pi }& A>\delta \end{array}\right.$$

(12)

where $t_1=\arcsin \left(\delta /A\right)$, in which $\delta$ is the unilateral freeplay size and A is the LCO amplitude in the nonlinear element. $A\le \delta$ indicates that the system has been vibrating inside the freeplay in the nonlinear element, that is, the elastic restoring moment and the equivalent stiffness are 0. When $A>\delta$, the equivalent stiffness increases from 0 to 1 with increasing oscillatory amplitude. Here, we define the amplitude ratio $A/\delta$ as the nondimensional LCO amplitude.

According to Equations (3) and (11), the state-space equations of the linear system can be written in the following form:

$$\left[\begin{array}{c}\dot{\mathbf{q}}\\ \ddot{\mathbf{q}}\end{array}\right]=\left[\begin{array}{cc}\mathbf{0}& \mathbf{I}\\ {\mathbf{M}}^{-1}\left(\frac{1}{2}\rho V^2{\mathbf{A}}_0-\mathbf{K}\right)& {\mathbf{M}}^{-1}\left(\frac{1}{2}\rho V{\mathbf{A}}_1\right)\end{array}\right]\left[\begin{array}{c}\mathbf{q}\\ \dot{\mathbf{q}}\end{array}\right]$$

(13)

where

$$\mathbf{M}=\left[\begin{array}{cc}{I}_{xx}& -I_{xy}\\ -I_{xy}& I_{yy}\end{array}\right],\mathbf{K}=\left[\begin{array}{cc}{K}_{\mathrm{B}}& \\ & K_{\mathrm{T}}\end{array}\right],\mathbf{q}=\left[\begin{array}{c}{q}_{\mathrm{B}}\\ q_{\mathrm{T}}\end{array}\right]$$

Here, we define the matrix $\mathbf{A}$ as follows:

$$\mathbf{A}=\left[\begin{array}{cc}\mathbf{0}& \mathbf{I}\\ {\mathbf{M}}^{-1}\left(\frac{1}{2}\rho V^2{\mathbf{A}}_0-\mathbf{K}\right)& {\mathbf{M}}^{-1}\left(\frac{1}{2}\rho V{\mathbf{A}}_1\right)\end{array}\right]$$

For systems with a single freeplay nonlinearity, when the nondimensional LCO amplitude in the nonlinear element is given, an equivalent stiffness $\overline{K}$ can be obtained from Equation (12). By replacing $K_{\mathrm{B}}$ or $K_{\mathrm{T}}$ in Equation (13), a linear flutter velocity can be calculated through the root-locus method. The specific process is as follows:

(1)

Assume a nondimensional LCO amplitude A in the nonlinear element and obtain the equivalent stiffness $\overline{K}$;

(2)

Replace $K_{\mathrm{B}}$ or $K_{\mathrm{T}}$ with $\overline{K}$ in Equation (13) and obtain a new matrix $\mathbf{A}$;

(3)

Calculate the eigenvalues of $\mathbf{A}$ with the variation of the free-stream velocity V;

(4)

As V increases, find the eigenvalue first crossing the imaginary axis from negative to positive, and record the free-stream velocity V, the crossing eigenvalue and corresponding eigenvector at this time;

(5)

Change the assumed nondimensional LCO amplitude and repeat steps (1) to (5) until the end.

Finally, we can obtain the relationship between the free-stream velocity and the nondimensional LCO amplitude; this reflects the nonlinear aeroelastic characteristics of the system. The frequency and phase information of the LCO will be described in more detail in Section 2.3.3. For systems with multiple freeplay nonlinearities, a discrete scanning method will be introduced in Section 3.

2.3.2. Numerical Integration Approach

According to Equations (5) and (11), the state-space equations of the nonlinear system can be written in the following form:

$$\left[\begin{array}{c}\dot{\mathbf{q}}\\ \ddot{\mathbf{q}}\end{array}\right]=\left[\begin{array}{cc}\mathbf{0}& \mathbf{I}\\ {\mathbf{M}}^{-1}\left(\frac{1}{2}\rho V^2{\mathbf{A}}_0\right)& {\mathbf{M}}^{-1}\left(\frac{1}{2}\rho V{\mathbf{A}}_1\right)\end{array}\right]\left[\begin{array}{c}\mathbf{q}\\ \dot{\mathbf{q}}\end{array}\right]+\left[\begin{array}{c}\mathbf{0}\\ -{\mathbf{M}}^{-1}{\mathbf{F}}_{\mathrm{NL}}\end{array}\right]$$

(14)

where

$${\mathbf{F}}_{\mathrm{NL}}=\left[\begin{array}{c}{M}_{\mathrm{B}}\\ M_{\mathrm{T}}\end{array}\right]$$

In this study, the Hénon–RK45 method [13,14] is used to integrate Equation (14). For simplicity, Equations (14) can be expressed as $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x})$. The fourth–fifth-order Runge–Kutta scheme (RK45) is used for numerical integration until the monitored variable ${\mathbf{x}}^{(k)}$ passes through the upper or lower boundary of freeplay. Assume that the value of the monitored variable ${\mathbf{x}}^{(k)}$ is ${{\mathbf{x}}_c}^{(k)}$ at the crossing point (${{\mathbf{x}}_c}^{(k)}$ is $\delta$ or $-\delta$), the last time instance before crossing is $t_0$, and the corresponding states are ${\mathbf{x}}_0$. Then, the Hénon method is used to make the monitored variable $x_k$ reach the discontinuity point ${{\mathbf{x}}_c}^{(k)}$ precisely via the following steps:

(1)

change the integral variable from time t to the monitored variable ${\mathbf{x}}^{(k)}$, and assume a new vector of the states $\mathbf{y}={[{\mathbf{x}}^T,t]}^T$; that is,

$$\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}{\mathbf{x}}^{(k)}}=\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}t}·\frac{\mathrm{d}t}{\mathrm{d}{\mathbf{x}}^{(k)}}=\left[\begin{array}{c}\mathbf{f}(\mathbf{x})/{\mathbf{f}}^{(k)}(\mathbf{x})\\ 1/{\mathbf{f}}^{(k)}(\mathbf{x})\end{array}\right]$$

(2)

integrate the above equation using the RK45 method from ${\mathbf{y}}_0={[{{\mathbf{x}}_0}^T,t_0]}^T$ to the crossing point ${\mathbf{y}}_c={[{{\mathbf{x}}_c}^T,t_c]}^T$ as the integral variable ${\mathbf{x}}^{(k)}$ is moved from ${{\mathbf{x}}_0}^{(k)}$ to ${{\mathbf{x}}_c}^{(k)}$;

(3)

restore the integral variable to the time t, and integrate Equation (14) from $\mathbf{x}(t_c)={\mathbf{x}}_c$ using the RK45 method until the next crossing point. Then repeat steps (1) to (3) until the end.

For systems with a single freeplay nonlinearity, the Hénon–RK45 method can effectively identify the discontinuities and achieve the purpose of reducing the error. For systems with multiple freeplay nonlinearities, there may be two or more integral variables that simultaneously cross the upper or lower boundary of the freeplay in the same time step. In this case, it is necessary to determine the order in which the relevant integral variables arrive at the boundary, and then the integral variable and time step can be controlled so that the state of the system can reach these boundaries successively.

Regarding the integration time step, it is recommended to set it around $1/(20{\omega }_{\max })$, where ${\omega }_{\max }$ is the maximum modal frequency of the model modes considered.

2.3.3. Harmonic Initial Conditions Derived from Frequency-Domain LCO Solutions

The relationship between the free-stream velocity and the nondimensional LCO amplitude can be finally obtained by using the DF method. However, in this process, most studies ignore the information of the corresponding flutter eigenvector, which contains the amplitude and phase information of LCOs.

According to Equation (13), here, we consider an n-dimensional linear system $\dot{\mathit{x}}(t)=\mathbf{A}\mathit{x}(t)$. The system response can be expressed as follows [41]:

$$\mathit{x}(t)=\mathbf{V}{e}^{\mathbf{L}t}{\mathbf{V}}^{-1}\mathit{x}(0)$$

(15)

where: $\mathbf{L}$ is the diagonal matrix of eigenvalues, the jth element of which is the jth eigenvalue ${\lambda }_j$ of $\mathbf{A}$; and $\mathbf{V}$ is the matrix of the eigenvectors, and its jth column is the jth eigenvector ${\mathit{\phi }}_j$, which satisfies $\mathbf{A}{\mathit{\phi }}_j={\lambda }_j{\mathit{\phi }}_j$. The response of the kth state variable ${\mathit{x}}^{(k)}(t)$ can be written as a superposition of all eigenvalue exponential responses:

$${\mathit{x}}^{\left(k\right)}(t)={\sum }_{j=1}^n\left(c_j{\mathrm{e}}^{{\lambda }_{j}t}\right)$$

(16)

where: the weight $c_j={\mathit{\phi }}_j^{\left(k\right)}{\sum }_{i=1}^n{\mathit{\xi }}_i^{(j)}{\mathit{x}}^{(i)}(0)$; ${\mathit{\phi }}_j^{\left(k\right)}$ represents the kth element of the column vector ${\mathit{\phi }}_j$; ${\mathit{\xi }}_i^{(j)}$ represents the jth element of ${\mathit{\xi }}_i$, which can be obtained by inverting the matrix $\mathbf{V}$, i.e., ${\mathbf{V}}^{-1}=\left[{\mathit{\xi }}_1,\dots ,{\mathit{\xi }}_n\right]$; and ${\mathit{x}}^{(i)}(0)$ is the ith element of $\mathit{x}(0)$.

For a real matrix $\mathbf{A}$, its eigenvalues ${\lambda }_j$ are real or complex. If ${\lambda }_j$ is complex, its complex conjugate ${\tilde{\lambda }}_j$ is also an eigenvalue of $\mathbf{A}$, and their eigenvectors are conjugate. Here, we denote the real and imaginary parts of ${\lambda }_j$ as $\mathbf{\Re }({\lambda }_j)$ and $\mathbf{\Im }({\lambda }_j)$, respectively. If $\mathbf{\Im }({\lambda }_j)\ne 0$, the superimposed exponential response of ${\lambda }_j$ and ${\tilde{\lambda }}_j$ can be written as follows:

$$c_j{\mathrm{e}}^{{\lambda }_{j}t}+{\tilde{c}}_j{\mathrm{e}}^{{\tilde{\lambda }}_{j}t}=2\left|c_j\right|{\mathrm{e}}^{\mathbf{\Re }({\lambda }_j)t}cos\left(\mathbf{\Im }({\lambda }_j)t+\angle c_j\right)$$

(17)

where $\left|c_j\right|$ is the magnitude of $c_j$, and $\angle c_j$ is the phase of $c_j$. Equation (17) represents an oscillatory response with angular frequency $\mathbf{\Im }({\lambda }_j)$. If $\mathbf{\Re }({\lambda }_j)<0$, the oscillatory amplitude decreases; if $\mathbf{\Re }({\lambda }_j)>0$, it will increase. If ${\mathbf{\Im }(\lambda }_j)=0$, the exponential response of ${\lambda }_j$ is $c_j{\mathrm{e}}^{\mathbf{\Re }({\lambda }_j)t}$, which diverges or decays exponentially over time.

For a linear aeroelastic system, when the free-stream velocity does not reach the flutter boundary, the real parts of the eigenvalues of the system matrix are all negative, and the system response decays. When the flutter boundary is reached, the real parts of one or two eigenvalues will change from negative to zero, and the exponential responses of the other eigenvalues with negative real parts will decay over time, eventually leaving only the response of the eigenvalues with zero real parts. If the imaginary parts of the remaining eigenvalues are zero, this indicates a static divergence; otherwise, it is a typical dynamic aeroelastic flutter. Here, we are mainly concerned with the latter.

Linear aeroelastic systems tend to be in two states: convergence or divergence, and there is almost no case of a critical constant-amplitude vibration; there is only an approximate constant-amplitude vibration that can be maintained for a longer time near the flutter boundary. However, for nonlinear systems, amplitude-stabilized LCO is relatively common. This can be approximated as the critical flutter state of the equivalent linear system from the point of view of the DF method. It should be noted that the state-space equation Equation (13) is only applicable to the stage of amplitude-stabilized LCO, and it cannot be applied to the previous transition stage. Here, we assume that the nonlinear system already has amplitude-stabilized LCO at time $t=0$, and we refer to $\mathit{x}(0)$ as a harmonic initial condition. Therefore, the response of the kth state variable ${\mathit{x}}^{(k)}(t)$ can be written as:

$${\mathit{x}}^{(k)}\left(t\right)=c_s{\mathrm{e}}^{{\lambda }_{s}t}+c_{s+1}{\mathrm{e}}^{{\lambda }_{s+1}t}$$

(18)

where: $c_s={\mathit{\phi }}_s^{\left(k\right)}{\sum }_{i=1}^n{\mathit{\xi }}_i^{(s)}{\mathit{x}}^{(i)}(0)$; $c_{s+1}={\mathit{\phi }}_{s+1}^{\left(k\right)}{\sum }_{i=1}^n{\mathit{\xi }}_i^{(s+1)}{\mathit{x}}^{(i)}(0)$; ${\lambda }_s=\mathrm{i}\omega$ and ${\lambda }_{s+1}=-\mathrm{i}\omega$ are the eigenvalues with zero real part in the equivalent linear system; $\omega$ is the angular frequency of the LCO; i is the imaginary unit.

Equation (18) indicates that the output response does not contain the exponential responses of the eigenvalue whose real part is not zero. Referring to Equation (16), $\mathit{x}(0)$ should satisfy the following conditions:

(1)

${\sum }_{i=1}^n{\mathit{\xi }}_i^{(j)}{\mathit{x}}^{(i)}(0)=0$ when $j\ne s$ and $j\ne s+1$;

(2)

${\mathit{\phi }}_s{\sum }_{i=1}^n{\mathit{\xi }}_i^{(s)}{\mathit{x}}^{(i)}(0)+{\mathit{\phi }}_{s+1}{\sum }_{i=1}^n{\mathit{\xi }}_i^{(s+1)}{\mathit{x}}^{(i)}(0)=\mathit{x}(0)$.

A general solution of the harmonic initial condition $\mathit{x}(0)$ is as follows:

$$\begin{array}{cc}\hfill \mathit{x}(0)& =a\left({\mathit{\phi }}_s{\mathrm{e}}^{\mathrm{i}{\varphi }_c}+{\mathit{\phi }}_{s+1}{\mathrm{e}}^{-\mathrm{i}{\varphi }_c}\right)\hfill \\ \hfill & =2a\left[\mathbf{\Re }({\mathit{\phi }}_s)\cos ({\varphi }_c)+\mathbf{\Im }({\mathit{\phi }}_s)\sin ({\varphi }_c)\right]\hfill \end{array}$$

(19)

where $a\in {\mathit{R}}^{+}$ and ${\varphi }_c\in [0,2\pi )$. The calculation of a is determined by the given LCO amplitude. For example, the given LCO amplitude in the bending DOF is $A_{\mathrm{B}}$, which means that $A_{\mathrm{B}}=2a\left|{\mathit{\phi }}_s^{(0)}\right|$, and we can calculate $a=A_{\mathrm{B}}/\left(2\left|{\mathit{\phi }}_s^{(0)}\right|\right)$. Here we call ${\varphi }_c$ the phase parameter of the harmonic initial condition. Its value is optional and different values correspond to different phases of the system at time $t=0$. In this paper, the value of ${\varphi }_c$ is set to 0 or $\pi /2$.

Substituting Equation (19) into Equation (18) yields

$${\mathit{x}}^{\left(k\right)}\left(t\right)=2a\left|{\mathit{\phi }}_s^{(k)}\right|cos\left(\omega t+{\varphi }_c+\angle {\mathit{\phi }}_s^{(k)}\right)$$

(20)

Therefore, for the all-movable fin in this paper, $\mathit{x}={[q_{\mathrm{B}},q_{\mathrm{T}},{\dot{q}}_{\mathrm{B}},{\dot{q}}_{\mathrm{T}}]}^T$, and the LCO amplitude ratio R in the bending and torsional DOFs can be written as:

$$R={\displaystyle \frac{2a\left|{\mathit{\phi }}_s^{(0)}\right|}{2a\left|{\mathit{\phi }}_s^{(1)}\right|}}={\displaystyle \frac{\left|{\mathit{\phi }}_s^{(0)}\right|}{\left|{\mathit{\phi }}_s^{(1)}\right|}}$$

(21)

which also represents the modal amplitude ratio of the equivalent linear system in the flutter state. The LCO phase difference $\Delta \varphi$ between the bending and torsional DOFs is as follows.

$$\Delta \varphi ={\varphi }_{\mathrm{B}}-{\varphi }_{\mathrm{T}}=\angle {\mathit{\phi }}_s^{(0)}-\angle {\mathit{\phi }}_s^{(1)}$$

(22)

The above analysis is based on the harmonic oscillation assumption of the DF method. The essence of the harmonic initial condition is to select a state point on the LCO trajectory as the initial condition of time-domain integration according to the amplitude and phase information of the generalized coordinates of the system under the assumed harmonic oscillation. However, this harmonic oscillation assumption is approximate for actual nonlinear systems. Therefore, the harmonic initial condition places the state of the system in a small neighborhood near the theoretical LCO trajectory; thus, the system response can be quickly attracted to the theoretical LCO trajectory if the LCO is stable. If the LCO is unstable, the system trajectory may gradually move away from the LCO trajectory [42], but there may be a short period of approximately constant-amplitude oscillation response before that. The introduction of the harmonic initial condition is convenient for quickly verifying the existence and stability of the frequency-domain LCO solutions, especially for LCO with a small attraction domain.

3. Results and Discussion

A typical small aspect ratio all-movable fin is studied in this paper; a diagram of this is shown in Figure 2, and the relevant parameters are shown in

Table 1. Structure parameters of the all-movable fin.

Parameter	Description	Value
m	Mass	$5.49\mathrm{kg}$
$I_{xx}$	Mass moment of inertia	$9.28\times {10}^{-2}\mathrm{kg}·{\mathrm{m}}^2$
$I_{yy}$	Mass moment of inertia	$3.39\times {10}^{-2}\mathrm{kg}·{\mathrm{m}}^2$
$I_{xy}$	Mass moment of inertia	$1.18\times {10}^{-2}\mathrm{kg}·{\mathrm{m}}^2$
$K_{\mathrm{B}}$	Bending spring stiffness	$5.0\times {10}^3\mathrm{N}·\mathrm{m}/\mathrm{rad}$
$K_{\mathrm{T}}$	Torsional spring stiffness	$2.0\times {10}^3\mathrm{N}·\mathrm{m}/\mathrm{rad}$
$\left[\overline{x},\overline{y}\right]$	Center-of-gravity position	$[0.013\mathrm{m},0.111\mathrm{m}]$

. Unless otherwise specified, ${\delta }_{\mathrm{B}}={0.1}^{\circ }$, ${\delta }_{\mathrm{T}}={0.1}^{\circ }$, the integration time step is $0.001$ s, and the initial conditions $\mathit{x}(0)$ of the time-domain integration are set to ${[10{\delta }_{\mathrm{B}},10{\delta }_{\mathrm{T}},0.0,0.0]}^T$. The air density and Mach are $0.122\mathrm{kg}/{\mathrm{m}}^3$ and $3.0$, respectively. The mode shapes are shown in Figure 3, and the linear flutter results obtained using the root-locus method are shown in

Table 2. Linear flutter results of the all-movable fin.

Parameter	Description	Value
${\omega }_1$	Modal frequency of mode 1	$34.3\mathrm{Hz}$
${\omega }_2$	Modal frequency of mode 2	$42.6\mathrm{Hz}$
$V_f$	Flutter velocity	$921.3\mathrm{m}/\mathrm{s}$
$Q_f$	Flutter dynamic pressure	$51.8\mathrm{kPa}$
${\omega }_f$	Flutter frequency	$37.9\mathrm{Hz}$

.

3.1. Single Freeplay Nonlinearity

According to the DF method, different LCO amplitudes correspond to different spring stiffnesses. Therefore, the effect of different torsional or bending stiffness on the linear flutter dynamic pressure was investigated first, as shown in Figure 4a and Figure 4b, respectively.

According to Figure 4 and Equation (12), the nonlinear aeroelastic characteristics of an all-movable fin with a single freeplay can be quickly and qualitatively judged [21]. The stiffness values given in Table 1 above correspond to state (i) in Figure 4. From Figure 4a we can see that when there is a freeplay nonlinearity in the torsional DOF, as the given nondimensional LCO amplitude $A_{\mathrm{T}}/{\delta }_{\mathrm{T}}$ increases, the bending stiffness $K_{\mathrm{B}}$ does not change, the equivalent torsional stiffness ${\overline{K}}_{\mathrm{T}}$ increases continuously and tends to the torsional stiffness of state (i), and the corresponding flutter dynamic pressure first decreases and then increases to the linear flutter dynamic pressure. Similarly, from Figure 4b, we can see that when there is a freeplay nonlinearity in the bending DOF, the flutter dynamic pressure decreases continuously to the linear flutter dynamic pressure with an increase in the given nondimensional LCO amplitude $A_{\mathrm{B}}/{\delta }_{\mathrm{B}}$.

Figure 5 shows the variation of nondimensional LCO amplitudes $A_{\mathrm{B}}/{\delta }_{\mathrm{B}}$ and $A_{\mathrm{T}}/{\delta }_{\mathrm{T}}$ with dynamic pressure in state (i); this is consistent with the above qualitative analysis results. In terms of the morphological characteristics of the LCO, for state (i), there are two different types of LCO when the freeplay nonlinearity exists in the torsional or bending DOF, and these are recorded as LCO-i and LCO-ii. It can be seen from Figure 5a that LCO-i contains both stable and unstable LCOs. The system response decays below $Q_{\mathrm{L}0}$, which is here referred to as the LCO occurrence boundary. When the dynamic pressure Q reaches $Q_{\mathrm{L}0}$, stable and unstable LCOs exist at the same time, and the system exhibits stable LCO motions in most initial conditions, while in some small initial conditions the system response decays. Then, the amplitude of the stable LCO motions in the nonlinear DOF will increase with increasing dynamic pressure, and this tends to infinity when the dynamic pressure Q tends to the linear-flutter dynamic pressure $Q_{\mathrm{LF}}$. When $Q>Q_{\mathrm{LF}}$, the system response will diverge in most initial conditions, while in some small initial conditions the system response decays. The discussion of initial conditions and the unstable LCOs will be expanded later. In summary, for a system with LCO-i characteristics, LCO may appear below the linear flutter boundary, bringing detrimental effects.

Unlike LCO-i, as shown in Figure 5b, LCO-ii only contains unstable LCO. The system response decays when $Q<Q_{\mathrm{LF}}$. After Q reaches $Q_{\mathrm{LF}}$, the response of the system decays or diverges depending on the magnitude of the initial disturbance. The minimum initial condition $q_{\mathrm{B}0}/{\delta }_{\mathrm{B}}$ required for divergence is also shown in Figure 5b with all other initial states at zero, and its value decreases as the dynamic pressure increases. In summary, the divergence boundary of a nonlinear system with LCO-ii characteristics is not lower than the flutter boundary of the linear system, indicating that the freeplay nonlinearity will instead bring favorable effects.

In fact, the minimum point of the dynamic pressure in Figure 4 divides the entire calculation region into two parts, that is, I and II. The position of the support stiffness in the nonlinear DOF determines the LCO characteristics of the fin. In region I, the system exhibits LCO-i characteristics, while in region II, the system exhibits LCO-ii characteristics.

In addition, we can calculate the variation of the linear flutter dynamic pressure with different bending and torsional stiffnesses, as shown in Figure 6. In Figure 6, along the horizontal direction, the bending stiffness of the fin changes, while along the vertical direction, the torsional stiffness changes. The blue arrow indicates that the linear flutter dynamic pressure decreases with increasing stiffness, while the magenta arrow indicates that it increases. The results shown in Figure 4a,b finally map to the vertical and horizontal lines through point (i) in Figure 6, respectively. For example, the vertical line through point (i) (marked as (i)-V) reflects that the fin with single torsion freeplay nonlinearity (marked as Ndof = T) exhibits LCO-i characteristics, corresponding to the results of Figure 4a and Figure 5a. Moreover, the horizontal line through the point (i) (marked as (i)-H) reflects that the fin with single bending freeplay nonlinearity (marked as Ndof = B) exhibits LCO-ii characteristics, corresponding to the results of Figure 4b and Figure 5b. Similarly, it can be found that the LCO characteristics of the fin in state (ii) are completely opposite to those of state (i).

Furthermore, the tangent points of the pressure contour lines and the vertical or horizontal lines can be connected into two orange straight lines; the slopes of these lines are recorded here as $\xi ={\xi }_{\mathrm{T}}$ and $\xi ={\xi }_{\mathrm{B}}$, where $\xi$ represents the ratio of torsional and bending stiffness. Another understanding of ${\xi }_{\mathrm{B}}$ and ${\xi }_{\mathrm{T}}$ is that if the bending stiffness $K_{\mathrm{B}}$ is constant, when $K_{\mathrm{T}}$ satisfies $K_{\mathrm{T}}/K_{\mathrm{B}}={\xi }_{\mathrm{T}}$, the flutter dynamic pressure $Q_f$ will achieve its minimum value; similarly, if the torsional stiffness $K_{\mathrm{T}}$ is constant, $Q_f$ achieves the minimum value when $K_{\mathrm{B}}$ satisfies $K_{\mathrm{T}}/K_{\mathrm{B}}={\xi }_{\mathrm{B}}$. When there is a torsion freeplay nonlinearity, $\xi ={\xi }_{\mathrm{T}}$ divides the entire computational domain into two parts, where $\xi \le {\xi }_{\mathrm{T}}$ corresponds to LCO-ii, and $\xi >{\xi }_{\mathrm{T}}$ corresponds to LCO-i. When there is a bending freeplay nonlinearity, $\xi ={\xi }_{\mathrm{B}}$ divides the entire computational domain into two parts, where $\xi <{\xi }_{\mathrm{B}}$ corresponds to LCO-i, and $\xi \ge {\xi }_{\mathrm{B}}$ corresponds to LCO-ii. Therefore, if the bending and torsional stiffnesses of the system are located between the lines $\xi ={\xi }_{\mathrm{T}}$ and $\xi ={\xi }_{\mathrm{B}}$, that is, ${\xi }_{\mathrm{T}}<\xi <{\xi }_{\mathrm{B}}$, the system exhibits LCO-i characteristics no matter whether there is a bending or torsion freeplay nonlinearity. Furthermore, from the results of the calculation, it can be seen that changing the bending and torsional stiffness does not change the values of ${\xi }_{\mathrm{B}}$ and ${\xi }_{\mathrm{T}}$; therefore, the LCO characteristics of the system can be changed by changing the ratio of bending and torsional stiffnesses.

Figure 5 also shows the results of stable and unstable LCOs obtained from time-domain simulations. The verification and calculation of unstable LCO in the time domain mainly use harmonic initial conditions to find the process of an approximately constant-amplitude vibration response. Figure 7a shows the time-domain response of the nondimensional LCO amplitude with only $q_{\mathrm{T}0}$ input. The results show that when $q_{\mathrm{T}0}/{\delta }_{\mathrm{T}}$ changes from $4.53$ to $4.5305$, the system response changes from convergence to LCO. This transition occurs around $q_{\mathrm{T}0}/{\delta }_{\mathrm{T}}=7.5$, which is exactly the same as the nondimensional amplitude of the unstable LCO calculated by the DF method (Figure 5a, $Q=48.0\mathrm{kPa}$). When only the $q_{\mathrm{T}0}$ input is non-zero, the system response first transitions to an intermediate state. If the amplitude of this intermediate state is greater than the amplitude of the unstable LCO, the system response decays, otherwise, it increases to the stable LCO. Figure 7b shows the time-domain response with a harmonic initial condition corresponding to the unstable LCO at 48 kPa. The phase parameter ${\varphi }_c$ of the harmonic initial condition is set to 0. Since the harmonic initial condition causes the system state to fall within a small neighborhood near the theoretical LCO trajectory, it is necessary to slightly adjust the harmonic initial condition (the whole is multiplied by a number around 1.0) to observe the unstable LCO for a long enough time. Unlike Figure 7a, the system response initially enters an approximately constant-amplitude vibration process, which corresponds to the unstable LCO. Then as time goes on, the system response will decay or enter the stable LCO. The LCO amplitude is consistent with the result in Figure 5a.

In general, for both LCO-i and LCO-ii, the initial conditions have a certain impact on the response of the system. With the introduction of harmonic initial conditions, the effects of initial conditions and unstable LCOs on the system response can be studied more clearly. Furthermore, we find that the unstable LCO is more like an obstacle. When the initial conditions are so small that the system cannot cross this obstacle, a system that is expected to enter the LCO or diverge would decay.

So far, we have provided a comprehensive understanding of the LCO characteristics of an all-movable fin with a single freeplay nonlinearity and different bending and torsional stiffnesses. Two different types of LCO are proposed, and the influence of parameter variation on the LCO characteristics is investigated. From this perspective, one possible LCO suppression approach is to change the LCO characteristics from LCO-i to LCO-ii so that no LCO occurs below the linear flutter boundary. On the contrary, the appearance of freeplay nonlinearity will bring favorable influence. The specific implementation method may be to adjust the supported stiffness radio or add counterweights.

3.2. Multiple Freeplay Nonlinearities

For a system with both bending and torsion freeplay nonlinearities, when continuing to use the DF method, it is necessary to assume two LCO amplitudes $A_{\mathrm{B}}$ and $A_{\mathrm{T}}$ to obtain the equivalent stiffnesses ${\overline{K}}_{\mathrm{B}}$ and ${\overline{K}}_{\mathrm{T}}$. Using the root-locus method, the corresponding flutter dynamic pressure, frequency, and eigenvector of the equivalent linear system can be obtained. The LCO amplitude ratio R in the bending and torsional DOFs can be calculated using Equation (21). Then, it is necessary to solve the problem that R is inconsistent with the assumed amplitude. In this section, a discrete scanning method is first introduced, and then the influence of the freeplay size ratio ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$ on the LCO characteristics in two different stiffness cases is analyzed.

3.2.1. Discrete Scanning Method

For the convenience of the expression, the function $f_{\mathrm{DF}}()$ is used here to represent the relationship between the equivalent stiffness and the nondimensional LCO amplitude, and the function $f_{\mathrm{R}}()$ is used to represent the relationship between the flutter modal amplitude ratio R and the stiffnesses $K_{\mathrm{B}}$ and $K_{\mathrm{T}}$ in a linear system. The relevant expressions are as follows:

$$\left\{\begin{array}{cc}\hfill {\overline{K}}_{\mathrm{B}}& =f_{\mathrm{DF}}(A_{\mathrm{B}}/{\delta }_{\mathrm{B}})\hfill \\ \hfill {\overline{K}}_{\mathrm{T}}& =f_{\mathrm{DF}}(A_{\mathrm{T}}/{\delta }_{\mathrm{T}})\hfill \\ \hfill R& =f_{\mathrm{R}}(K_{\mathrm{B}},K_{\mathrm{T}})\hfill \end{array}\right.$$

(23)

The LCO solutions should satisfy the following equation to achieve alignment of the amplitude of each nonlinear element.

$$A_{\mathrm{B}}/A_{\mathrm{T}}=f_{\mathrm{R}}\left[f_{\mathrm{DF}}(A_{\mathrm{B}}/{\delta }_{\mathrm{B}}),f_{\mathrm{DF}}(A_{\mathrm{T}}/{\delta }_{\mathrm{T}})\right]$$

(24)

The left side of Equation (24) represents the assumed LCO amplitude ratio, which is a monotonically decreasing function of $A_{\mathrm{T}}$; the right side is the flutter modal amplitude ratio R of the equivalent linear system, which is a composite function composed of $f_{\mathrm{R}}()$ and $f_{\mathrm{DF}}()$. Solving the LCO characteristics of all-movable fins with both bending and torsion freeplay nonlinearities can be transformed into the following two equivalent problems:

• given ${\delta }_{\mathrm{B}}$ and ${\delta }_{\mathrm{T}}$, assuming $A_{\mathrm{B}}$, find a suitable $A_{\mathrm{T}}$ to satisfy Equation (24);
• given ${\delta }_{\mathrm{B}}$ and ${\delta }_{\mathrm{T}}$, assuming $A_{\mathrm{T}}$, find a suitable $A_{\mathrm{B}}$ to satisfy Equation (24).

Here, we choose the former. According to Equation (14), $f_{\mathrm{DF}}()$ is a monotonically increasing function of $A_{\mathrm{T}}$. For the subject of this study, $f_{\mathrm{R}}()$ is shown in Figure 8. Given different $A_{\mathrm{B}}$ values that correspond to different ${\overline{K}}_{\mathrm{B}}$ values, the minimum point of R is roughly on the yellow dotted line. When $A_{\mathrm{B}}$ is fixed, if $K_{\mathrm{T}}$ is smaller than the value on the yellow dotted line, the right side of Equation (24) is a monotonically decreasing function of $A_{\mathrm{T}}$; otherwise, the right side of Equation (24) first decreases and then increases with increasing $A_{\mathrm{T}}$. For the latter case, the LCO amplitude $A_{\mathrm{T}}$ corresponding to the minimum R is recorded as $A_c$ here.

In contrast to previous studies, the ranges of the bending and torsional stiffnesses of the all-movable fin considered in this paper are very large so that its nonlinear aeroelastic characteristics can be comprehensively studied. The torsional stiffness may be relatively low, but no static divergence occurs. In this case, because the function on the right side of Equation (24) may not be monotonous, the previous iterative methods may cause a loss of solutions, and the iteration may not even converge. A more general method is to scan the entire solution interval to find the solutions of Equation (24). However, in the discrete scanning method proposed in this paper, a critical point $A_{\mathrm{c}}$ is used to divide the solution interval into two parts: $A_{\mathrm{T}}\in [A_{\min },A_{\mathrm{c}}]\bigcup [A_{\mathrm{c}},A_{\max }]$, where $A_{\min }$ and $A_{\max }$ are the specified minimum and maximum LCO amplitudes in the torsional DOF, respectively. As shown in Figure 9, when $A_{\mathrm{T}}\in [A_{\mathrm{c}},A_{\max }]$, the dichotomy method can be directly used because the left and right sides of the equation have different monotonicity. When $A_{\mathrm{T}}\in [A_{\min },A_{\mathrm{c}}]$, the scanning method is used. Compared with scanning the whole solution interval, this discrete scanning method improves the calculation efficiency while avoiding the loss of solutions.

Equation (24) can be expressed as:

$$\frac{A_{\mathrm{B}}/{\delta }_{\mathrm{B}}}{A_{\mathrm{T}}/{\delta }_{\mathrm{T}}}\frac{{\delta }_{\mathrm{B}}}{{\delta }_{\mathrm{T}}}=f_{\mathrm{R}}(f_{\mathrm{DF}}(A_{\mathrm{B}}/{\delta }_{\mathrm{B}}),f_{\mathrm{DF}}(A_{\mathrm{T}}/{\delta }_{\mathrm{T}}))$$

(25)

where the nondimensional LCO amplitudes $A_{\mathrm{B}}/{\delta }_{\mathrm{B}}$ and $A_{\mathrm{T}}/{\delta }_{\mathrm{T}}$ can be considered as two variables, and the freeplay size ratio ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$ is a parameter. Therefore, in this paper ${\delta }_{\mathrm{T}}$ is set to ${0.1}^{\circ }$, and the freeplay size ratio ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$ is changed by using different ${\delta }_{\mathrm{B}}$ values. The variation of the nondimensional LCO amplitudes with the dynamic pressure is obtained with different freeplay size ratios. It should be noted that although ${\delta }_{\mathrm{T}}$ here is a definite value, no matter how it changes, as long as ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$ remains unchanged, the solution of the nondimensional LCO amplitudes $A_{\mathrm{B}}/{\delta }_{\mathrm{B}}$ and $A_{\mathrm{T}}/{\delta }_{\mathrm{T}}$ corresponding to the equation will not change; the only thing that changes is the corresponding LCO amplitudes $A_{\mathrm{B}}$ and $A_{\mathrm{T}}$. Therefore, the nonlinear flutter characteristics of the system under different freeplay sizes ${\delta }_{\mathrm{B}}$ and ${\delta }_{\mathrm{T}}$ can be obtained by establishing the variation of the nondimensional LCO amplitudes of the system under different freeplay size ratios ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$.

Here, $A_{\min }/{\delta }_{\mathrm{T}}$ is set to $1.2$ to avoid static aeroelastic divergence, and $A_{\max }/{\delta }_{\mathrm{T}}$ is set to 1000. The calculation error of $A_{\mathrm{T}}/{\delta }_{\mathrm{T}}$ in the dichotomy method is required to be less than ${10}^{-6}$. The scanning interval of $A_{\mathrm{T}}{\delta }_{\mathrm{T}}$ is $0.01$. When the scanning method is used, the computational domain is divided into n intervals: $[A_{\mathrm{T}}^i,A_{\mathrm{T}}^{i+1}]$ ($i=1,...,n$). Then, we successively calculate the values $F(A_{\mathrm{T}}^i)$. If $F(A_{\mathrm{T}}^i)·F(A_{\mathrm{T}}^{i+1})\le 0$, then we use the dichotomy method to find the solution of $F(A_{\mathrm{T}})=0$ while $A_{\mathrm{T}}\in [A_{\mathrm{T}}^i,A_{\mathrm{T}}^{i+1}]$.

3.2.2. Case 1: $K_{\mathrm{B}}=5.0\times {10}^3\mathrm{N}·\mathrm{m}/\mathrm{rad}$, $K_{\mathrm{T}}=2.0\times {10}^3\mathrm{N}·\mathrm{m}/\mathrm{rad}$

In Case 1, a system with only a bending freeplay nonlinearity will exhibit LCO-ii characteristics; with only a torsion freeplay nonlinearity, it will exhibit LCO-i characteristics. Figure 10 shows the LCO solutions obtained by the discrete scanning method at different freeplay size ratios ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$ and the corresponding time-domain bifurcation diagrams. The dynamic pressure range for the time-domain analysis was from 15 kPa to the divergence dynamic pressure. Here, different LCO branches are defined according to the solving interval and the continuity of the LCO solutions with the parameters.

From the frequency-domain calculation results, we can find that there are two main LCO branches for each ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$. Combined with the time-domain results, it can be found that when the freeplay size ratio is small, stable LCO motions appear below the linear flutter dynamic pressure. In addition, irregular motions with limited amplitude are observed in a certain region of dynamic pressure. With the increase of the freeplay size ratio, the dynamic pressure region of the LCO motions decreases, and the region of the irregular motions increases, but the amplitude tends to decrease. When the freeplay size ratio increases to $0.5$, the nondimensional amplitude of the vibration is very small and almost entirely located near the freeplay boundary.

After the LCO solutions are obtained, the equivalent stiffness corresponding to the nondimensional LCO amplitude in the bending and torsional DOFs can be calculated. Therefore, each set of LCO nondimensional amplitudes $[A_{\mathrm{B}}/{\delta }_{\mathrm{B}},A_{\mathrm{T}}/{\delta }_{\mathrm{T}}]$ will correspond to a point $[K_{\mathrm{B}},K_{\mathrm{T}}]$ in Figure 6. Finally, complete LCO solutions at different ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$ can be drawn in the form of Figure 11. The yellow dotted line in Figure 11 is the same as in Figure 8. For each given $A_{\mathrm{B}}$, it divides the solution domain into upper and lower parts, which correspond to the solution intervals $[A_{\mathrm{c}},A_{\max }]$ and $[A_{\min },A_{\mathrm{c}}]$. Finally, the results of branches A and B located on either side of the yellow line are obtained. The LCO solutions of a system with a single torsion or bending freeplay nonlinearity can be seen in the Figure 5 map to the right and the upper dotted lines in Figure 11, respectively.

When ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$ decreases, branch A gradually approaches the upper dotted line, and the LCO characteristics of a system with multiple freeplay nonlinearities are more similar to a system with only a torsion freeplay nonlinearity, that is, LCO-i. Similarly, when ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$ increases, the LCO characteristics tend to LCO-ii. Comparing the frequency-domain results at ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}=0.01$ in Figure 5a and Figure 10, it can be inferred that when ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$ decreases to 0, the LCO of branch A finally will become the stable LCO branch in Figure 5a, while branch B becomes the unstable LCO branch. As ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$ increases to infinity, the LCO solutions in the bending DOF of branch A will finally become the unstable LCO branch in Figure 5b, while the LCO amplitude of branch B will be reduced to become negligible.

In summary, when the difference between the two freeplay sizes is large, the nonlinear aeroelastic characteristics of the system are mainly determined by the freeplay with the larger value. In addition, expressing the LCO solutions in the form of Figure 11 allows more concise and clear observation of the variation rule of the LCO characteristics of the system with the freeplay size ratio. This also contains the variation of the nondimensional LCO amplitude with dynamic pressure.

Further, we construct the Poincaré section by plotting ${\dot{q}}_{\mathrm{B}}$ against $q_{\mathrm{B}}$ at time instances when $q_{\mathrm{T}}={\delta }_{\mathrm{T}}$ and ${\dot{q}}_{\mathrm{T}}>0$. For the Poincaré sections to be meaningful, the total simulation time was set to 200 s at each dynamic pressure state, and the last 100 s of data was selected for analysis. In Figure 12, dynamic pressures between $44.0$ kPa and $47.6$ kPa were chosen for the analysis at ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}=0.01$. Because the system response decays, the response trajectories do not intersect with the Poincaré section when the dynamic pressure is less than $44.2$ kPa. A closed curve appears suddenly at $44.4$ kPa, and at about $46.2$ kPa, this becomes a single point. Figure 13 shows a plot of the response trajectory of the system in the $q_{\mathrm{T}}$–${\dot{q}}_{\mathrm{T}}$ phase plane at $45.0$ kPa. The response trajectory of the system is always within a torus. It can be clearly seen that a Hopf bifurcation occurs before $44.4$ kPa. The system response changes from decaying motions to quasi-periodic motions somewhere between $44.2$ and $44.4$ kPa, and it transforms into LCO motions at about $46.2$ kPa. From the frequency-domain LCO solutions in Figure 10, we can find two LCO branches that grow from a point around 38 kPa, which have a similar amplitude. Until approximately 46 kPa, the difference in the torsional LCO amplitude increases sharply due to the turning point of branch B. It can also be seen in Figure 11 that the two LCO branches start from one point and then gradually move away from each other when ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}=0.01$. Therefore, the appearances of such quasi-periodic motions are probably caused by the interaction between the two LCO branches, which should have similar amplitudes and irrational frequency ratios [42]. Once the two LCO branches move far away from each other, these quasi-periodic motions will be destroyed, and the response trajectories will be attracted by one of the LCO branches.

Figure 14 shows plots of the variations of the LCO frequency $\omega$ and the phase difference $\Delta \varphi$ with dynamic pressure. The frequency of the quasi-periodic motion refers to the frequency of the component with the highest proportion. The frequency-domain calculation results show that the frequencies of the two LCO branches are relatively close until the turning point of the torsional LCO frequency at about $46.0$ kPa; however, the phase differences between the responses in the bending and torsional DOFs of the two branches have completely opposing trends. After the system enters the LCO at $46.2$ kPa, the frequency and phase of the time-domain LCO response are consistent with the LCO solutions of branch A. Figure 15 shows the spectral-analysis results of the time-domain response in the yellow area. Not surprisingly, we find that the LCO response at 47 kPa mainly contains the components of one frequency and its multiples, while the quasi-periodic motion at 45 kPa contains some irreducible frequencies. The frequency of the component with the highest proportion is consistent with the frequencies of the two LCO branches (Figure 14). However, some other frequency components cannot be obtained by the method in this paper; these may be generated by the coupling effects of LCO branches with similar amplitudes and frequencies but different phase differences. Another interesting phenomenon is that when the two LCO branches are very close, this quasi-periodic motion disappears, such as when the dynamic pressure is less than 44 kPa at ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}=0.01$ or less than 28 kPa at ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}=0.1$.

To further study the transformation process of this quasi-periodic motion to the LCO, the time-domain response of systems with different dynamic pressures and initial conditions at ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}=0.1$ was calculated, and the results are shown in Figure 16. When the default initial conditions ${[10{\delta }_{\mathrm{B}},10{\delta }_{\mathrm{T}},0.0,0.0]}^T$ are given, a quasi-periodic motion occurs at $42.27$ and $43.53$ kPa. If a harmonic initial condition corresponding to branch A is given, the response trajectory is immediately attracted to the LCO at $43.53$ kPa. However, at $42.27$ kPa, even given similar harmonic initial conditions, the amplitude of the system fluctuates increasingly as time goes on, and it eventually enters quasi-periodic motion. This phenomenon shows that with increasing dynamic pressure, the torsional LCO amplitude difference between the two branches expands, and this leads to the weakening of the coupling effect between the LCO branches. At $45.26$ kPa, when a harmonic initial condition corresponding to branch B is given, the system response decays. However, if the amplitude of the initial condition becomes $1.02$ times its original value, the system will eventually enter a stable LCO. This comparison indicates that the LCO of branch B is unstable at this dynamic pressure, and with sufficient initial amplitude, the system will be attracted to the stable LCO of branch A. In addition, we can find that the system at $45.26$ kPa with larger initial conditions has strong amplitude fluctuations at the beginning, but as time goes on, the fluctuations become progressively smaller, which is the complete opposite of the phenomenon at $42.27$ kPa. The above interesting phenomena reflect the weakening of the coupling effect of the two LCO branches and the expansion of the attraction domain of the stable LCO during the transformation process of the system from quasi-periodic motion to the LCO.

Figure 17 shows variations of the Poincaré section with dynamic pressure at ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}=0.1$ and ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}=0.5$. Figure 18 shows the variation of the largest Lyapunov exponent with time [43,44], and it finally converges to a positive value with increasing time step N, indicating that the response of the system is a chaotic motion. Combined with Figure 12, it can be clearly seen that with increasing ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$, the torus is gradually broken, and the system response gradually transitions from LCO and quasi-periodic motion to chaotic motion.

3.2.3. Case 2: $K_{\mathrm{B}}=5.0\times {10}^3\mathrm{N}·\mathrm{m}/\mathrm{rad}$, $K_{\mathrm{T}}=1.5\times {10}^3\mathrm{N}·\mathrm{m}/\mathrm{rad}$

For the LCO characteristics of a system with a single freeplay nonlinearity, cases 2 and 1 are completely opposite. A system with only a bending freeplay nonlinearity will exhibit LCO-i characteristics, while one with only a torsion freeplay nonlinearity will exhibit LCO-ii characteristics. Figure 19 plots the LCO solutions obtained by the discrete scanning method and the corresponding time-domain bifurcation diagrams at different freeplay size ratios ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$. Figure 20 plots the distribution of LCO solutions.

According to the solving interval and the continuity of the LCO solutions with the parameters, three LCO branches are defined here. In Figure 20, when ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}\le 0.25$, LCO branches A and B are located on either side of the yellow dotted line (two solution intervals) and converge at the top right. In contrast, branch C is relatively independent. With the decrease of ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$, Branch C moves toward the right dotted line, and branches A and B move toward the yellow dotted line. The intersection point of branches A and B moves down and left, which means the right boundary of the dynamic pressure region of LCO branches A and B in Figure 19 is shifted to the left. If ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$ decreases to 0, LCO branches A and B would disappear and the system would exhibit LCO-ii characteristics similar to Figure 5b, as determined by the torsion freeplay. From the frequency-domain LCO solutions in Figure 19, there is no matching LCO solution in a certain dynamic pressure region when ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}\le 0.25$. The corresponding time-marching results show that the system response decays in this dynamic pressure region. Branch C is in the right of the region. This dynamic pressure region keeps shrinking with increasing ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$, and it disappears completely at ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}=0.5$. Through further computational analysis, it was found that branches B and C meet and undergo recombination when ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$ is approximately $0.2685$. The relevant results are presented in Figure 20.

Later, with the increase of ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$, branch B moves downward and eventually disappears in Figure 20, while branches A and C move toward the upper dotted line. The system exhibits LCO-i characteristics similar to Figure 5a, determined by the bending freeplay. Comparing the trend of $A_{\mathrm{B}}/{\delta }_{\mathrm{B}}$ at ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}=5.0$ in Figure 20 with the trend of LCO-i in Figure 5a, it can be seen that branch A corresponds to the unstable LCO, and branch C corresponds to the stable LCO. The corresponding time-domain bifurcation diagram in Figure 20 also agrees well with the calculation results of branch C.

Next, we specifically look at the time-domain bifurcations of the system under different freeplay size ratios in Figure 19. At ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}=0.1$, the system response decays and eventually reaches a static equilibrium position when the dynamic pressure is lower than the linear flutter boundary. Due to the existence of the unstable LCO branch C, when the dynamic pressure is slightly higher than the linear flutter dynamic pressure, the system response does not diverge. This is because the initial conditions at this time are not of great enough amplitude to excite the system into a divergent state. The time-response results given in Figure 21 further verify the calculation results of the LCO branch C. However, the other two LCO branches A and B in the frequency-domain calculation results are not reflected in the corresponding time-domain bifurcation diagrams.

Considering that the current initial conditions might not excite the system to enter the LCO, we increased the amplitude of the initial conditions. The corresponding time-domain bifurcation diagrams are shown in Figure 22. When the initial conditions increase from 10 times the freeplay size to 100 times, quasi-periodic or chaotic motions appear between dynamic pressures of $25.4$ and $31.3$ kPa. However, the existence of such vibrations with a non-zero amplitude is uncertain because the system enters a decaying motion under some dynamic pressures between $25.4$ and $31.3$ kPa. In addition, a quasi-periodic motion occurs alone at $35.7$ kPa. In general, the existence of LCO branches A and B theoretically causes the system to exhibit more quasi-periodic or chaotic motions, but such vibrations are uncertain and have great dependencies on the initial conditions.

For the case of ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}\ge 0.25$, the system mainly exhibits quasi-periodic or chaotic motions when the dynamic pressure is less than about 38 kPa. The nondimensional vibration amplitude also becomes larger as the dynamic pressure increases. As ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}$ increases, the nondimensional vibration amplitude becomes smaller, but it changes more dramatically at around 38 kPa. Finally, the jumping phenomenon of the LCO appears at ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}=5.0$. When the dynamic pressure is higher than about 38 kPa but below the linear flutter boundary, the system exhibits LCO motions except for within the dynamic pressure region in which the response decays at ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}=0.25$.

Considering that ${\delta }_{\mathrm{B}}/{\delta }_{\mathrm{T}}\ge 0.5$, the nondimensional LCO amplitude increases with increasing dynamic pressure in the bending DOF, but it first decreases and then increases in the torsional DOF. There is a maximum value of the nondimensional LCO amplitude in the torsional DOF at around $38.0$ kPa; this should be carefully considered because structural damage can occur here.

4. Conclusions

This paper focused on the nonlinear aeroelastic characteristics of a small aspect ratio supersonic all-movable fin with single or multiple freeplay nonlinearities. The conclusions that can be drawn from the numerical results can be summarized as follows.

(1)

The DF method has shown good applicability to the analysis and solution of the LCO characteristics of the system with single or multiple freeplay nonlinearities. In this paper, a discrete scanning method based on the DF method was established to solve the problem of amplitude alignment, and this can avoid missing solutions and give consideration to efficiency.

(2)

Systems with a single freeplay nonlinearity exhibit two different types of LCO. A system with LCO-i characteristics has stable LCO motions below the linear flutter boundary, while a system with LCO-ii characteristics only has unstable LCO motions above the linear flutter boundary, and the divergence boundary is not lower than the linear flutter boundary. From a flutter suppression perspective, LCO-ii is more beneficial.

(3)

Systems with bending and torsion freeplay nonlinearities exhibit complex dynamical behaviors, such as LCO and quasi-periodic and chaotic motions. The initial conditions have a great influence on the system response. The introduction of harmonic initial conditions will effectively improve the efficiency of using time-domain integration to verify frequency-domain LCO solutions, especially for those with a small attraction domain.

(4)

The support stiffness at the root of the fin and freeplay size ratio have impressive effects on the LCO characteristics. The ratio of bending stiffness to torsional stiffness determines the LCO characteristics of a fin with a single freeplay nonlinearity. For a fin with multiple freeplay nonlinearities, when the difference between the two freeplay sizes is large, its nonlinear aeroelastic characteristics are mainly determined by the freeplay with the larger value.

The small aspect ratio fin in this paper is treated as a rigid fin because its elastic modes have higher frequencies and have little influence on the flutter characteristics of the fin. However, for the flexible fin with a high aspect ratio, the elastic modes cannot be ignored. Further research is needed to investigate the influence of the elastic modes. In addition, damping is not involved in this paper, and the suppression effect of damping on the vibrations of the fin with multiple freeplay nonlinearities needs further study.