High-Fidelity Fin–Actuator System Modeling and Aeroelastic Analysis Considering Friction Effect

Featured Application: A high ﬁdelity aeroelastic model of a typical ﬁn–actuator system is established, in which the friction model is more accurate. The inﬂuence of freeplay and friction on the system stability is analyzed using the time-domain method and frequency-domain method, and some suggestions for ﬂutter suppression design are put forward. Abstract: Both the dynamic characteristics and structural nonlinearities of an actuator will affect the ﬂutter boundary of a ﬁn–actuator system. The actuator models used in past research are not universal, the accuracy is difﬁcult to guarantee, and the consideration of nonlinearity is not adequate. Based on modularization, a high-ﬁdelity modeling method for an actuator is proposed in this paper. This model considers both freeplay and friction, which is easy to expand. It can be directly used to analyze actuator characteristics and perform aeroelastic analysis of ﬁn–actuator systems. Friction can improve the aeroelastic stability, but the mechanism of its inﬂuence on the aeroelastic characteristics of the system has not been reported. In this paper, the LuGre model, which can better reﬂect the friction characteristics, was integrated into the actuator. The inﬂuence of the initial condition, freeplay, and friction on the aeroelastic characteristics of the system was analyzed. The comparison of the results with the previous research shows that oversimpliﬁed friction models are not accurate enough to reﬂect the mechanism of friction’s inﬂuence. By changing the loads, material, and geometry of contact surfaces, ﬂutter can be effectively suppressed, and the power loss caused by friction can be minimized.


Introduction
Flutter is a critical problem in modern flight vehicle design. To pursue lighter weight, wings and control surfaces tend to be more flexible, which may lead to flutter. Flutter instability will seriously affect aircraft performance and flight safety [1]. The conventional aeroelastic analysis of the control surface often focused only on the surface itself, and the actuator supporting the surface was regarded as a linear torsion spring with constant stiffness [2], as shown in Figure 1a. In fact, many experiments have proved that the actuator itself has dynamic characteristics, especially electromechanical servo actuators. Electromechanical servo actuators have good maintainability, low static power consumption, and storage convenience, so they have been widely used in missile control fin position servo systems [3]. Thus, they are the modeling objects of this paper. An electromechanical actuator consists of multiple components that can be classified as the moment of inertia (MOI), elastic spring, and damper, as shown in Figure 1b, which means that the actuator should not be over-simplified as a linear spring. Moreover, the stiffness of the actuator will change according to the change in the external load, which is called dynamic stiffness. In addition, the actuator has control algorithms and electrical and mechanical nonlinearities such as current/voltage saturation, freeplay, and friction [4]. The analysis of the fin-actuator Both the dynamics and structural nonlinearity of the actuator will greatly affect the aeroelastic characteristics of the system [6,7]. In the field of aeroelasticity, there are some studies that model the actuator and perform flutter analysis while considering the actuator dynamics and structural nonlinearities. Yehezkely et al. [5] performed an aeroelastic analysis of the missile fin with a nonlinear pneumatic actuator and proposed a flutter suppression method. Paek and the team [8] analyzed the flutter characteristics of the control fin of a rocket considering dynamic actuator characteristics. Shin et al. [6] established a basic framework for obtaining the aeroelastic stability results of a fin-actuator system considering the actuator's dynamic characteristics. They found that the flutter boundary could be increased in this way. In the same year, Shin et al. [9] investigated the same system with the consideration of actuator nonlinearities, such as freeplay, backlash, and transmission error. The results show that both the freeplay and backlash may influence the actuator dynamic stiffness. They also found that the nonlinearity and gear reductions can greatly affect the limit cycle oscillation (LCO) characteristics. Yang et al. [7] obtained the aeroelastic characteristics of the fin-actuator system containing preload freeplay and found that an appropriate preload angle can suppress LCOs. Zhang and his team [10] analyzed the influence of actuator parameters on the dynamic stiffness. The results show that factors such as inertia of the motor rotor, connection stiffness between different levels of reducers, reduction ratio of each level of reducer, and the damping at the actuator fin shaft have significant effects. It can be seen that the studies above only considered freeplay or other nonlinearities in the actuator and rarely analyzed the friction. Friction actually has a great influence on aeroelastic properties, which exists widely and can be used to improve aeroelastic properties. Therefore, it was innovative and necessary to add a friction nonlinearity module to the model in this paper. Moreover, the previous work was only for a specific actuator, which is often too simplified, and the reliability of the model was not easy to guarantee. Furthermore, they often used the dynamic differential equations of the actuator to carry out modeling and analysis. The coupling between these equations often makes the modeling process cumbersome and the solving process complex. Therefore, a new actuator modeling method is needed. Based on modularization, a high-fidelity and systematic modeling method for an actuator is formed that can ensure accuracy, ease of use, and consider freeplay and friction at the same time. It can be easily extended to different configurations of actuators. Using this model, the characteristics of the actuator can be obtained directly for the analysis of aeroservoelasticity.
In addition to the natural friction at the joints of the system, engineers, based on experience, take measures to increase friction in order to suppress flutter. For example, friction sheets [11] or friction rings [12] were added to some actuators to increase friction. Friction plays a positive role in the aeroelastic system, but friction in servomechanism systems dissipates power [13] and may cause position signal steady-state errors or oscillations. The added friction sheet is often set empirically. So far, the mechanism of friction on the aeroelastic characteristics of the system is still unclear, and relevant research or analysis has not been reported. How can we reduce friction while ensuring the effect of friction on flutter stability? How can we achieve a balance between ensuring aeroelastic stability and reducing power loss? These are the issues that this article explores in the study of the friction influence mechanism.
Friction is a nonlinear and complex phenomenon that is determined by the relative motion between the two contact surfaces. Friction is related to interface material and operating conditions, such as loads acting on the contact [14]. Introducing nonlinearity significantly complicates all aspects of already complicated aeroelastic equations of motion, to the point that friction is typically omitted from a flutter-related analytical workflow or approximated by a constant value. Some researchers use simple friction models to analyze the influence of friction on aeroelastic characteristics. Khalak [15] researched stability criteria with consideration of mechanical damping, which was assumed constant. Griffin et al. [16] considered friction as a hysteretic spring with damping to study the blade response. The blade was simplified as a single degree of freedom system, and then the maximum response amplitude of it was analyzed using the Ritz method. Whiteman and the team [17] studied the possibility of dry friction for flutter suppression. The results showed that the aeroelastic stability boundary can be enhanced by changing the position of friction and the geometry of the interfaces. Mignolet et al. [18] validated for the first time that friction may play a stabilizing role in LCOs that occur after flutter. Tan et al. [19] studied the stability of an aeroelastic system with 1.5 degrees of freedom. The results showed that different friction stiffness has an influence on the critical aerodynamic damping ratio and the stability boundary of the system, but it has little influence on the stability boundary. Lu et al. [20] studied the effect of the Coulomb friction on airfoil aeroelastic stability by a numerical method and the harmonic balance method. The results showed that Coulomb friction can increase the flutter velocity, and it is affected by the initial condition. Wayhs-Lopes et al. [21] regarded Coulomb friction as a linear damping ratio and analyzed the influence of the Coulomb friction on LCOs caused by symmetric and asymmetric freeplay. The results show that friction is dissipative, but it is not enough to restrain LCOs. The friction model used in the above work is relatively simple and is a preliminary exploration of the friction effect mechanism and cannot fully reflect the friction characteristics. Therefore, the prediction of the stability boundary of a fin-actuator aeroelastic system is not enough. In fact, in engineering, the friction can be changed by adjusting the preload, the roughness of the friction surface, lubrication, material, temperature, and other factors that cannot be represented fully by simple models. Therefore, a model that can more accurately reflect the friction characteristics is introduced into the aeroelastic analysis model to improve the model's accuracy, study further the influence mechanism of the above factors on the aeroelastic stability, and give some guidance suggestions.
The friction models that are currently used are mainly separated into two main types: the static model and the dynamic model [22]. The static model may include different components such as viscous friction, Coulomb friction, the Stribeck effect, and static friction. However, the static friction model cannot describe the influence of friction when the contact surfaces are relatively static. Some experiments show that some phenomena of friction can only be reproduced by dynamic models with memory, such as hysteresis, rate-dependence, and pre-displacement. The dynamic friction model can reflect the friction phenomenon more accurately, and therefore has more application value. Many researchers have developed a variety of dynamic models, mainly including the Dahl model [23], the LuGre model [24,25], the Leuven model [26,27], the Generalized Maxwell-Slip (GMS) model [28,29], and two-state friction model (2SEP) [30]. The Dahl model is not accurate enough when the velocity is low, and it cannot capture stick-slip motion. The LuGre friction model extends the Dahl [24,25,[31][32][33][34]. It can reflect the complex processes of friction static and dynamic characteristics, including viscous friction, Stribeck effect, variable static friction, memorial friction, and pre-sliding displacement. The disadvantage of the LuGre model is that the hysteresis behavior does not have nonlocal memory characteristics in the presliding regime, which has attracted some attention in the field of precise control. The Leuven model can reflect this characteristic, but it is difficult to implement. The same author also proposed the GMS model. Compared with the Leuven model, the GMS model can more accurately reflect the physical mechanism of friction, but the model is also more complex, which leads to difficulties in implementation or identification. The 2SEP model considers to a greater extent the micro-scale elastoplasticity, which is not the focus of this paper. A comparison of these models is shown in Table 1. It is worth noting that when choosing a model, a proper trade-off must be made between model fidelity and implementation. The accuracy of the LuGre model is sufficient for the problem we want to study. The parameters' physical meaning of the LuGre model is clear, and thus can easily match experimental data. It is a first-order model, which is more convenient to use and analyze than other second-order models [35]. Various benefits make this model widely used in control systems [36][37][38][39][40], which means it is easy to integrate into the actuator model based on SIMULINK, launched by MathWorks company, and it is also convenient for subsequent work with which to perform active flutter suppression. Using more complex models to consider the effect of friction in a more micro and accurate way can be an open question for future investigations.
The LuGre model contains only six parameters, which are σ 0 , σ 1 , σ 2 , F c , F s , and V s . Many researchers have proposed effective LuGre parameters identification methods [31,[41][42][43]. Based on interpreting the asperities contact, the LuGre model uses elastic spring-like bristles with damping to accurately represent the asperities between two contact surfaces. σ 0 is the stiffness of the assumed elastic spring, and σ 1 is the damping that is associated with micro-displacement. The bristles deform because of the action of external force, and the sum of the force generated by bristles' deformation plus the viscous friction represents the friction force [44]. The average deformation of the bristles reflects the characteristics of dynamic friction. Once the bristles reach their maximum deformation, which also means friction force has been larger than the maximum static friction force F s , slip starts, as illustrated by the schematic representation in Figure 2.
In this article, a high-fidelity simulation model of a typical electromechanical actuator was established, considering its structural characteristics, details of servomechanism, freeplay, and friction. The actuator was divided into modules to solve the problem that the previous model was not universal and extensible. Then the fin-actuator system was established by connecting the actuator with the fin using the branch mode method. The unsteady aerodynamic force acted on the control fin. The influence on the system of freeplay and friction was studied using numerical simulation. Compared with previous literature, the friction model used in this paper could better reflect the real physical phenomenon, and the mechanism of friction effect on aeroelastic system was studied. Based on the results, some guidelines applicable to engineering projects are put forward. This paper is organized as follows. Section 2 shows the study object of this article, which is a typical fin-actuator system with a friction sheet. Section 3 presents the high-fidelity modeling method of a fin-actuator system considering aerodynamics. The frequency domain and time domain flutter analysis methods are established in Section 4. Section 5 introduces the results and discussion of the stability analysis of the fin-actuator system, considering freeplay and friction nonlinearity of the actuator, and proposes some rules for the design reference of flutter suppression.

Study Object
The electromechanical actuator used in this paper was composed of a DC motor, a primary gear reducer, a lead screw-nut pair, a fork, a controller, and a sensor. Figure 3a shows the whole structure of the fin-actuator system. The fork between the screw and the fin shaft drove the fin. A friction sheet was installed under the fork as shown in Figure 3b. The mutual extrusion of the sheet and the fork produced the friction force, which caused the friction torque relative to the center of the fin shaft. One bolt fastened one section of the friction sheet, and the other bolt, which can provide a pressing force changed artificially, pressed the other section. Adjusting the pressure applied by the screw to the friction sheet, changing the materials in contact, or changing the roughness of the contact surface between the sheet and the fork could regulate friction. The MOI of the fin was 0.0105 kg·m 2 . A branch mode analysis of the fin was performed, which was adopted for aeroelastic analysis, and the first two elastic natural frequencies and the corresponding mode shapes of branch 1 are listed in Table 2.  The MOI of the fin was 0.0105 kg·m 2 . A branch mode analysis of the fin was performed, which was adopted for aeroelastic analysis, and the first two elastic natural frequencies and the corresponding mode shapes of branch 1 are listed in Table 2.

Modeling of the System
The electromechanical actuator installed on a missile provides sufficient torque to rotate the control fin against aerodynamic loads and allows the fin to follow guidance commands. The actuator must have appropriate stiffness to offset the disturbance to meet the stability and accuracy requirements. Different support stiffnesses will change the natural mode of the fin. Therefore, the normal mode cannot effectively represent the fin. To model the motion of the system, the substructure technique was introduced.

Modeling of the Fin Structure
Here, the Gladwell branch modal synthesis method was applied and used to model the fin structure. This method divides the modes into two branches-namely, the elastic branch and the rotational rigid branch [10,45,46], as illustrated by the schematic diagram in Figure 4. Note that the elastic modal branch, which is branch 1, was obtained by adding a spring to the root of the fin shaft in the direction of the x axis and the fin was clamped in the y and z directions. The MOI of the fin was 0.0105 kg·m 2 . A branch mode analysis of the fin was performed, which was adopted for aeroelastic analysis, and the first two elastic natural frequencies and the corresponding mode shapes of branch 1 are listed in Table 2.

Modeling of the System
The electromechanical actuator installed on a missile provides sufficient torque to rotate the control fin against aerodynamic loads and allows the fin to follow guidance commands. The actuator must have appropriate stiffness to offset the disturbance to meet the stability and accuracy requirements. Different support stiffnesses will change the natural mode of the fin. Therefore, the normal mode cannot effectively represent the fin. To model the motion of the system, the substructure technique was introduced.

Modeling of the Fin Structure
Here, the Gladwell branch modal synthesis method was applied and used to model the fin structure. This method divides the modes into two branches-namely, the elastic branch and the rotational rigid branch [10,45,46], as illustrated by the schematic diagram in Figure 4. Note that the elastic modal branch, which is branch 1, was obtained by adding a spring to the root of the fin shaft in the direction of the x axis and the fin was clamped in the y and z directions.

Modeling of the System
The electromechanical actuator installed on a missile provides sufficient torque to rotate the control fin against aerodynamic loads and allows the fin to follow guidance commands. The actuator must have appropriate stiffness to offset the disturbance to meet the stability and accuracy requirements. Different support stiffnesses will change the natural mode of the fin. Therefore, the normal mode cannot effectively represent the fin. To model the motion of the system, the substructure technique was introduced.

Modeling of the Fin Structure
Here, the Gladwell branch modal synthesis method was applied and used to model the fin structure. This method divides the modes into two branches-namely, the elastic branch and the rotational rigid branch [10,45,46], as illustrated by the schematic diagram in Figure 4. Note that the elastic modal branch, which is branch 1, was obtained by adding a spring to the root of the fin shaft in the direction of the x axis and the fin was clamped in the y and z directions.

x+Kx=F
(1) for which the details of matrix M, C, K, and F are shown in Appendix A. Moreover, the generalized coordinates vector is x T = [q e , δ , where q e is the displacement vector in generalized coordinates of elastic modes, and δ is the rotation angle of the actuator output shaft.

Modeling of the Unsteady Aerodynamics
The aerodynamic forces were calculated by the panel method. The reference Mach number was 2.5, and the air density was 1.225 kg/m 3 The aerodynamic influence coefficient (AIC) matrix was generated using ZONA7 [48]. Figure 5 shows the aerodynamic model of the fin. The generalized aerodynamic forces in frequency domain can be written as follows: The details of q ∞ , A, and x are shown in Appendix A.

State-Space Form of the Fin Model
Combining Equations (1) and (2), yields: The fin model in state-space form is: where the state vector is x, x a , x a is the vector of m augmented states. Input vector is u s = M t , and output vector is y s T = [δ, . δ]. The size of x is (r + 1 + m) × 1, where r represents the number of elastic natural modes of branch 1. Other matrices can be found in Appendix A.

Modeling of the Electromechanical Actuator
To model the actuator, some generally accepted simplifications are used concerning power supplies, electromagnetic circuits, machine structures, and physical properties. The motor model directly relates the current and the output torque. The gearing mechanism and the lead screw-nut pair can be treated as mass-damping-spring systems to study their properties.

Model of DC Motor
Modeling of a DC motor ignores the details of the three-phase electromagnetic field and establishes the relationship between current and torque directly. The Lorentz force of the motor is given by: where K m is the torque coefficient. The motor winding voltage balance equation is given by: where L is the inductance, R is the resistance, C e is the back electromotive force (EMF) coefficient, u m is the driving voltage, K i is the current feedback coefficient, θ 0 is the angle of motor shaft. The dynamic equation of motor rotor is given by where J m is the MOI of the rotor, b m is the viscous damping coefficient of the rotor, T 0 is the load torque acting on the motor shaft, and its sign is defined as positive in the same direction as θ 0 .

Model of the Gear Pair
What is connected to the DC motor is the gear pair. It should be noticed that the connection stiffness k m exists between the motor and the gear pair.
The force diagram of gears is shown in Figure 6. The dynamic equations of gears are given by: where θ 1 and θ 2 are the rotation angles of the two gears, respectively; J 1 and J 2 are the MOIs of the two gears, respectively; b 1 and b 2 are the viscous damping coefficients of the two gears, respectively; T 1 and T 2 are the external moments acting on gears 1 and 2, respectively; r 1 and r 2 are the radii of the two gears, respectively; and F is the force between gears 1 and 2. The meshing stiffness between gears is k g . The equation reflecting the relationship is given by:

Model of the Screw-Nut Pair and Fork
The structure of the screw-nut pair which is connected to the gear pair is shown in Figure 7. The force of the screw-nut pair and fork model was complicated. Freeplay existed in the entire actuator. The screw-nut pair and fork model was the section closest to the actuator output shaft, with dynamic characteristics that will have a huge impact on the dynamic stiffness of the actuator. Therefore, in this paper freeplay was considered here as the comprehensive freeplay of the whole actuator. The friction produced by the friction sheet is described in this section, compared with the friction in other parts of the actuator, which was negligible.
where θ 3 is the screw rotation angular displacement, b 3 is the viscous damping coefficient of the screw, η is the efficiency of the screw-nut pair, λ is the helix angle of the lead screw, r sg is the radius of the screw, k c is the comprehensive stiffness between the screw-nut pair, x n is the displacement of the nut, x f is the displacement of the fork, e is half the width of symmetrical freeplay at the fin shaft, L bc is the length of the fork, δ e is the effective rotation angle of the actuator output shaft, and T d is the output torque of the fork.
Notice that the connection stiffness between the gear pair and the screw-nut pair and fork is k z . Therefore T 3 is given by As described in Section 1 (Introduction), a friction sheet is installed on the fork. Before the actuator starts to work, a certain pressure is applied to the friction sheet. The friction torque of the LuGre model is [49]: where T f is the friction torque; σ 0 and σ 1 have been explained in Section 1. σ 2 is the viscous damping coefficient, and z is an internal state variable that is introduced to represent the average deflection of all the bristles. The positive function g( . δ) represents the Stribeck effect and Coulomb friction, which is related to several factors, such as temperature and material. F c is Coulomb friction, and F s corresponds to the maximum static friction. V s is the Stribeck velocity, which reflects how fast friction approaches F c .
The torque of the actuator to the fin is given by

Model of the Controller and Sensor
The fin deflection angle is measured with an angular displacement sensor that gives position feedback to the controller. The transfer function of the sensor can be considered as a second-order element, which is given by: The controller is used to adjust the input signal so that the output signal of the actuator can keep the track of the command signal. The most commonly used control law is to control the error between the command signal and the measured signal using proportional (P), integral (I), and derivative (D) terms [50]. The conventional actuator is a closed-loop system, and its controller usually includes amplification, filtering, and a PID control. The principle diagram of a PID control is shown in Figure 8, where δ cmd is the command signal.

Model of the System
The actuator drives the fin to rotate and provides support stiffness for the fin in the meantime. The output of the actuator is the input of the fin and vice versa. Therefore, M t = T fin . A schematic diagram of the fin-actuator system is presented in Figure 9.

Dynamic Stiffness
The describing function (DF) method is widely used in control engineering. Nikolay Mitrofanovich Krylov and Bogoliubov developed the DF method in the 1930s [51,52]. They used the average method and introduced the concept of equivalent linearization. Ralph Kochenburger extended the DF method [53]. When using the DF method, the Fourier coefficients of the nonlinear term need to be calculated and divided into two terms, which are respectively interpreted as the equivalent stiffness and damping terms related to the excitation signal. In order to ensure the quasi-linearity of the actuator, the excitation signal amplitude needs to be small enough. The actuator dynamic stiffness is essentially the transfer function of the quasi-linear actuator model obtained by the DF method when the input amplitude is small. Dynamic stiffness performance reflects the robustness of the actuator to disturbance. It refers to the dynamic torque required to be applied to produce a unit angle at the output shaft of the actuator. Dynamic stiffness reflects the characteristic that the support stiffness provided by the actuator to the fin will change with the frequency of the external load. As shown in Figure 10, the dynamic stiffness of the actuator is given by: where T load (ω) = T load e iωt is the external load on the actuator output shaft, and θ out (ω) = θ out e i(ωt−ϕ) is the rotation angle of the actuator output shaft. The dynamic stiffness can be expressed by the amplitude phase curve varying with the input frequency ω; k δδ is the amplitude and ϕ is the phase. According to the relationship between the output and input power spectral densities (PSDs), the actuator dynamic stiffness can be calculated [6]. Taking the external load as the input signal x(t) and the deflection angle as the output signal y(t), the transfer function is: where P xy (ω) is cross PSD of x(t) and y(t), and P xx (ω) is auto PSD of x(t).
In numerical simulation, the excitation signal is a constant-amplitude sweep frequency signal, and the sweep frequency range should include the frequency range of interest. Another method to obtain dynamic stiffness is to combine a step sine sweep (SSS) signal with a least square (LS) frequency response function (FRF) estimation algorithm. The SSS signal contains considerable energy to excite out all the characteristics of the actuator, and the LS algorithm is an averaging process in a sense and essentially plays an effective filtering role, but this method is time-consuming.

Flutter Analysis with V-g Method
Using dynamic stiffness to conduct the frequency-domain flutter analysis. The force in Equation (1) is the aerodynamic force in Equation (2), and ignoring C yields: The generalized stiffness matrix K is a complex matrix due to the dynamic stiffness. Therefore, the equation above using the V-g method [54] is: where ig is the added artificial complex structural damping. Let the eigenvalue of Equation (19) be λ = ω 2 /(1 + ig), and substituting λ and k = ωb/V into Equation (19) yields: In Equation (20), both k and K(ω) are related to the oscillation frequency ω. Therefore, an iterative solution is required for this eigenvalue problem (see Figure 11).

Time Domain Method
The whole time domain mathematical fin-actuator model was integrated into SIMULINK, and the in-built "ode23s" was used for numerical simulation. When the time domain simulation starts, release the fin from an initial angular position, then change the flight speed, and the angular position response of the fin can be seen. Flutter has occurred when the oscillation is diverging. The fin may maintain a certain angle of oscillation with constant frequency before large-angle divergence occurs. This phenomenon is called LCO. Although the time domain method is time-consuming for nonlinear problems, the results are more accurate than the frequency domain method. Therefore, the subsequent analysis mainly focuses on the time domain. The frequency domain is used as a comparative verification.

Aeroelastic Analysis Results and Discussion
In this article, the impact of the actuator structural nonlinearities such as freeplay and friction on the fin-actuator aeroelastic characteristics were studied, and other nonlinearities were excluded.

Identification of Actuator Parameters
The actuator parameters are listed in Table 3. The parameters of the friction model are identified according to the method in Reference [41]. Several constant velocity experiments were carried out to obtain the static parameters of the model. The dynamic parameters are obtained during the pre-sliding process.

Preliminary Flutter Results for No-Freeplay Gap (e = 0) and No-Friction
The fin was released from the initial fin shaft deflection angle of 1 • . Figure 12

Preliminary Flutter Results for No-Freeplay Gap ( = 0 e ) and No-Friction
The fin was released from the initial fin shaft deflection angle of 1°. Figure 12 shows the time domain response of the fin under aerodynamic force. Comparison has been made between the responses of the fin shaft angle at two flight velocities. Flutter occurred at 765 m/s, and when the flight velocity was lower than 765 m/s, the fin shaft angle response was convergent.  Figure 13. Eighty-two frequency points were selected from 40 Hz to 120 Hz. The corresponding dynamic stiffness was adopted to obtain the flutter frequency and velocity of the system at each frequency point. Then a line was drawn with slope 1 that indicated the frequencies of abscissa and ordinate were equal. The frequency corresponding to the intersection of the flutter frequency line and the line with slope 1 was the real flutter frequency, which was 55.1077 Hz. The velocity corresponding to the real flutter frequency was the real flutter velocity, which was 769.2596 m/s. The frequency-domain and time-domain results were in good agreement; the discrepancy between which was of the order of 0.5%. This result shows that the calculation method of dynamic stiffness was correct and the flutter result of fin-actuator system calculated using dynamic stiffness was accurate enough. The initial conditions will affect the stability of nonlinear systems. The half width of freeplay at the fin shaft was 0.2 • , and the initial fin shaft deflection angles δ 0 were 0.5 • , 1 • , and 2 • , respectively. The corresponding actuator dynamic stiffness under different δ 0 was calculated, and the results are shown in Figure 14.  Figure 13. Eighty-two frequency points were selected from 40Hz to 120Hz. The corresponding dynamic stiffness was adopted to obtain the flutter frequency and velocity of the system at each frequency point. Then a line was drawn with slope 1 that indicated the frequencies of abscissa and ordinate were equal. The frequency corresponding to the intersection of the flutter frequency line and the line with slope 1 was the real flutter frequency, which was 55.1077 Hz. The velocity corresponding to the real flutter frequency was the real flutter velocity, which was 769.2596 m/s. The frequency-domain and time-domain results were in good agreement; the discrepancy between which was of the order of 0.5%. This result shows that the calculation method of dynamic stiffness was correct and the flutter result of finactuator system calculated using dynamic stiffness was accurate enough. The initial conditions will affect the stability of nonlinear systems. The half width of freeplay at the fin shaft was 0.2°, and the initial fin shaft deflection angles 0  were 0.5°, 1°, and 2°, respectively. The corresponding actuator dynamic stiffness under different 0  was calculated, and the results are shown in Figure 14. It can be discovered from Figure 14a that as δ 0 increases, the gain of the dynamic stiffness also increased greatly, but it was not proportional to the increase in δ 0 , and the increase in phase was very small. In Reference [55], the actuator dynamic stiffness was obtained by experiments and simulations in a series of cases including changing the initial preload. The conclusion was consistent with the change of the initial angle in this paper. However, flutter analysis was not carried out in Reference [55]. Actually, δ 0 will also affect the aeroelastic characteristics of the system, which can be seen in Figure 15b. δ 0 will affect the speed at which the fin begins to oscillate. As δ 0 increases, the oscillation occurs at smaller flight speed, and the oscillation amplitude becomes larger. However, different δ 0 have little effect on A/e-namely, the ratio of oscillation amplitude of the fin to the width of freeplay at the same flight speed. It only has some influence on A/e at the beginning of oscillation. When δ 0 is 1 • , the comparison of the frequency domain and time domain results is shown in Figure 15, which shows that these results were highly consistent and the maximum error did not exceed 1.2%. As the flight speed increased, the oscillation amplitude suddenly increased at a certain speed, which means that flutter occurred. The oscillation frequency increased slowly within a small range with the increase in flight speed.

Influence of Different Freeplay with
Keep δ 0 at 1 • and change the freeplay. When the half-widths of freeplay are 0.1 • , 0.2 • , and 0.5 • , calculate the corresponding actuator dynamic stiffness as shown in Figure 16.
The experimental results in Reference [55] showed that the dynamic stiffness of the actuator with a smaller freeplay gap will reach a larger value. What is not mentioned in the conclusion of Reference [55] is that the phase was not affected much by the freeplay gap, although it can be seen from the experimental results. Freeplay is symmetrical nonlinearity; thus, the phase of its describing function was zero. Figure 16a also confirms that the smaller the freeplay is, the larger the amplitude is, while the phase changes little. Figure 16b shows the time domain results of aeroelastic analysis. The torsional response amplitude increases as the flight velocity increases, and the oscillation amplitude increases as the width of freeplay increases, which are mentioned in Reference [56] too. Moreover, different widths of freeplay also have an influence on the velocity at which the fin begins to oscillate. It should be noted that the ratio of oscillation amplitude of the fin to e was almost not affected by e just after reaching the oscillating velocity. However, as the flight speed increased, the smaller the freeplay width, the larger the ratio.

Aeroelastic Characteristics of the System with Freeplay and Friction
The LuGre model is represented by several nonlinear functional equations, including six parameters in which σ 0 and σ 1 are two dynamic parameters, and σ 2 , F c , F s , and V s are four static parameters. Analyzing the effect of friction on the flutter characteristics is challenging since coupling exists among the six parameters [57,58]. However, it is undoubtedly an effective way to analyze the aeroelastic characteristics with the help of the concept of dynamic stiffness and the physical meaning of the parameters. σ 0 and σ 1 are related to material properties. If a lubricant exists, σ 2 indicates the viscous properties of the lubricant. For dry friction, σ 2 is pretty small and close to zero [26]. The pressure influences on F s , F c , and V s [59]. The variable-parameter analysis is conducted to evaluate the influence of each parameter. Table 4 summarizes the range of friction parameters and the factors that parameter depends principally upon. As can be seen, lubricant is an important factor, but there is no lubricant added to the friction sheet of the actuator generally. Since the LuGre parameters depend on the material properties and working conditions, when the parameters are changed for stability analysis, they cannot be divorced from reality. Here, the upper and lower bounds of the parameters were derived from the literature [26,41,[60][61][62]. For comparison, the standard value of each friction parameter in Table 3 was multiplied by factors of 0.1 and 10, after which the parameter values were still in the range of Table 4.  Figure 17 shows the impact of the variation of six friction parameters on the actuator dynamic stiffness when δ 0 is 1 • and e is 0.2 • . It was discovered that changing σ 2 , F s , or V s had little effect on the dynamic stiffness, while changes in σ 0 , σ 1 , and F c had a greater effect. It can be observed that, generally, friction had a huge influence on both gain and phase of actuator dynamic stiffness. This is because the friction model used in this paper is complex and highly nonlinear. Friction has been oversimplified in some literatures. The actuator in Reference [64] was equivalent to a spring with friction nonlinearity, and the relationship between displacement and restoring force is given by: The equivalent stiffness of this nonlinear spring is calculated by the describing function method. The results in Reference [64] show that friction only affects the amplitude of actuator stiffness, which has been modified in this paper according to Figure 17. Figure 18 shows the comparison of the frequency-domain and time-domain results. The friction parameters are standard values in Table 3 and keep e at 0.2 • . Due to the nonlinearity of the friction model, as the oscillation amplitude increased, the difference between the time-domain and frequency-domain results became larger, but it did not exceed 4.4%. The oscillation frequency increased slightly as the flight speed increased. It is mentioned in Reference [21] that adding friction makes the frequency of LCO decrease slightly. Figure 19 compares time-domain frequencies in Figures 15b and 18b, which confirms this conclusion.   Figure 20 shows the time-domain results, which also confirm that the six parameters had different effects on the system's stability. In Figure 20, we can draw a conclusion that F s and V s have little influence on the aeroelastic characteristics. This corresponds to the fact that they have little effect on dynamic stiffness. As F s increases, the flight speed at which the fin begins to oscillate first increases and then decreases slightly, while the oscillation amplitude at the same flight speed remains almost unchanged. As V s increases, the flight speed at which the fin begins to oscillate also increases slightly, while the oscillation amplitude at the same flight speed is almost unchanged. σ 1 seems to be the parameter that has the greatest influence on the aeroelastic characteristics in the variable-parameter analysis. It can be seen that increasing σ 1 can greatly improve the stability of the system. When σ 1 is 10 times the standard value, the fin deflection angle always converges within the flight speed range of interest. This parameter is difficult to quantify in actual engineering. It can be increased by changing the contact surface geometry, that is, increasing the roughness of the contact surface. More work is needed to quantify and standardize σ 1 in the future.
F c , which stands for Coulomb friction, can be studied as a separate friction model. Compared with what is mentioned in Reference [20]-that increasing the maximum friction torque of Coulomb friction can increase the flutter speed and reduce the amplitude of the pitch response-we discover that the influence of F c on system stability is complicated. As F c increases, the speed at which the fin starts to oscillate first increases and then decreases, while the oscillation amplitude at the same flight speed increases. Therefore, F c is not the larger the better. An optimal F c can be obtained by appropriately designing the lubrication, pressure, and the geometry of the contact surface between the friction sheet and fork, so that the system is stable and friction is reduced.
In addition to σ 1 and F c , σ 0 and σ 2 also have a non-negligible effect on the aeroelastic characteristics. The smaller the σ 0 , the greater the flight speed at which the fin starts to oscillate, which is the desired effect. However, once the oscillation starts in this case, the amplitude will increase sharply. As mentioned earlier, σ 2 is related to lubricants, so it is not the focus of this study. It can be seen that larger σ 2 will slightly increase the speed at which the fin starts to oscillate, but it will also cause the amplitude of oscillation to increase at the same flight speed.
As mentioned earlier, a large number of studies have used the Coulomb model. Some researchers have discovered the limitations of Coulomb friction. They have noticed that with the increase of friction, the stable region increased, but the stable boundary changed little [19]. Wayhs-Lopes et al. studied the effects of both freeplay and Coulomb friction last year, pointed out that although friction is dissipative, it is not enough to restrain LCO [21]. However, the analysis in this article is sufficient to show that a single factor may have limitations. The synergistic effect of multiple parameters can suppress the LCO in the low-speed range, and it is also effective in increasing the flutter velocity and reducing the oscillation amplitude.

Conclusions
In this article, a high-fidelity model of a fin-actuator system with a typical configuration was established based on modularization, in which a friction sheet was modeled utilizing the LuGre model, which can better reflect the friction characteristics. The whole modeling method can be used for dynamic analysis and aeroelastic analysis.
The aeroelastic characteristics of the system were analyzed using frequency-domain and time-domain methods. The results of the two methods were highly consistent. The time-domain method can deal with highly nonlinear situations, while the frequencydomain method based on the DF method is suitable for weak nonlinear systems.
Initial conditions, freeplay, and friction will all affect the stability of the nonlinear system. As the initial deflection angle δ 0 increases, the gain of the dynamic stiffness will be larger, but the increase was not proportional, and the increase of the phase was minor. The oscillation occurred at smaller flight speed, and the oscillation amplitude A became larger. Different δ 0 have little effect on the ratio of A to the half width of freeplay e at the same flight speed. Dynamic stiffness of the actuator with a smaller freeplay gap will reach a larger value, while the phase was not affected much by freeplay gap. The oscillation amplitude increased as the width of freeplay increased.
Friction affects the amplitude and phase of actuator dynamic stiffness greatly. Adding friction makes the frequency of LCO decrease slightly. The LuGre friction model has six parameters. Friction sheets that do not use lubricant, σ 2 , which reflects the viscosity of the lubricant, were not the focus of this study. F s and V s have little influence on the aeroelastic characteristics, which corresponds to the fact that they have little effect on dynamic stiffness. σ 1 is the parameter that has the greatest influence on the aeroelastic characteristics. Increasing σ 1 can greatly improve the stability of the system. However, σ 1 is difficult to quantify in actual engineering. The influence of F c on system stability is complicated. As F c increased, the speed at which the fin started to oscillate first increased and then decreased, while the oscillation amplitude at the same flight speed increased. Therefore, F c was not better when larger. An optimal F c can be obtained by appropriately designing the lubrication, pressure, and the geometry of the contact surface between the friction sheet and fork, so that the system is stable and friction is reduced. The smaller the σ 0 , the greater the flight speed at which the fin starts to oscillate, which was the desired effect. However, once the oscillation started in this case, the amplitude increased sharply. A single factor may have limitations. The synergistic effect of multiple parameters can suppress the LCO in the low-speed range, and it is also effective in increasing the flutter velocity as well as reducing the oscillation amplitude.
Author Contributions: Data curation, formal analysis, investigation, software, validation, visualization, and writing (original draft) were conducted by J.L. Conceptualization, methodology, and project administration were conducted by Z.W. Writing (review and editing) was conducted by both J.L. and Z.W. Resources and supervision were conducted by both Z.W. and C.Y. All authors have read and agreed to the published version of the manuscript.