Integration of Freeplay-Induced Limit Cycles Based On a State Space Iterating Scheme

Featured Application: To derive all possible limit cycles of a freeplay nonlinear dynamic system using time integrations with massive combinations of initial conditions efﬁciently and to analyze the limit cycle oscillation (LCO) stability of the initial conditions. Abstract: Time integration is commonly used to obtain accurate system responses, such as the limit cycle oscillations (LCOs) for an aeroelastic system with freeplay. However, the integrations that start with various initial conditions (I.C.s) are usually studied case by case, so only a few system states can possibly be focused on. This paper proposes a state space iterating (SSI) scheme to ﬁnd LCO solutions using time integration by using another method. First, a large number of arbitrary I.C. cases are used for time integrations, but only a very short integration time is required for each I.C. case. Second, system behaviors are depicted visually through a method that combines a modiﬁed Poincar é map and Lorenz map, in which the LCO solutions are found as ﬁxed points via visual inspections. To verify the SSI scheme’s ability to ﬁnd LCOs, a typical plunge–pitch wing section is established numerically. Time integrations with both the classic scheme and the proposed SSI scheme are carried out. The LCO results of the SSI scheme are well-aligned with those from the classic scheme. The SSI scheme visualizes the patterns of system responses using arbitrary I.C. cases and analyzes the LCO stability, which provides more mathematical insights into an aeroelastic system with freeplay.


Introduction
Nonlinearities inevitably occur in most real dynamic systems and can induce abundant nonlinear behaviors such as limit cycle oscillations (LCOs), quasi-periodic motions, and chaos. Aeroelasticity is a discipline that usually studies aircraft dynamic systems that involve elastic structures subject to aerodynamic forces, inertia forces, and/or control systems [1]. Over the decades, scholars have been aware that when freeplay-one of the most ubiquitous and typical nonlinearities on the connecting parts of an aeroelastic system-occurs, persistent and complex nonlinear behaviors will arise and affect the system's safety [2,3]. Therefore, the development of effective and efficient methods to analyze the nonlinear behaviors in an aeroelastic system with freeplay has become a major subject of concern.
Analyses of nonlinear behaviors involve at least two important aspects: (1) obtaining the nonlinear responses in the form of time histories or some parameters that can reconstruct the responses (e.g., the amplitudes and frequencies for an LCO) to understand "what is happening for the system", and (2) obtaining as many and as detailed physical and mathematical insights as possible to suggest "why it happens or how important it is"-for example, analyzing the stability of an LCO or determining the type of an LCO bifurcation.
Various methodologies are available for obtaining nonlinear responses, such as the harmonic balance (HB) method, continuation, and time integration, among which time of the graph of state patterns, which is very efficient. The LCO stability could also be confirmed in the same way.
Because this new scheme requires an iteration procedure in state space, we named it the state space iterating (SSI) scheme and applied it to a typical plunge-pitch wing section subject to unsteady aerodynamic loads. The proposed SSI scheme yields two LCO branches for the wing section, which are well aligned with those derived from the classic time integration using the Hénon-RK45 method. The highlights of the SSI scheme are as follows: (1) The spatial patterns of system states are clearly pictured, so one can easily find LCO solutions and confirm the stabilities directly and visually; (2) the integration time for each I.C. case is reduced to a few LCO periods, so a large number of I.C. cases would not decrease the calculation efficiency; and (3) the proposed method can be easily extended to all other kinds of structural nonlinearities and applied to an aeroelastic system with higher dimensions. This paper is organized as follows: Section 2 presents the state space equations for a general aeroelastic system with freeplay and introduces the classic Hénon-RK45 method. Section 3 introduces the proposed SSI method and explains how the Poincaré map and Lorenz map are modified, applied, and finally form the SSI scheme together. A plungepitch aeroelastic model is established in Section 4. Then, an LCO analysis based on the Hénon-RK45 method is implemented in Section 5. Time integrations based on the SSI scheme are carried out in Section 6, along with the discussions and a comparison to the Hénon-RK45 method.

State Space Equations of an Aeroelastic System with Freeplay
The equations of the motion of an aeroelastic system with freeplay can generally be expressed in the form of a set of piecewise linear state space equations (e.g., Equation (5) in [3], Equation (8) in [4], and Equation (1) in [13]): where x is n x system states; A 0 is an n x × n x constant matrix; b is an n x × 1 constant vector, and f FP is the freeplay function, reflecting the nonlinear connection between a displacement and a structurally restoring force. Assume that only one freeplay involved in the present work takes the k-th state x k as its input; then, the freeplay nonlinearity can be expressed as: (2) where K lin is the underlying linear stiffness of the freeplay and 2δ is the amount of freeplay. In the present paper, the k-th degree of freedom (DoF) of the system is called the "freeplay DoF" and x k is called the "freeplay state". Equation (2) defines two non-smooth boundaries, which are known as the two "freeplay boundaries": Σ 1 : x k = −δ and Σ 2 : x k = δ. These two boundaries separate the state space into three subdomains: (1) S 0 , where |x k | ≤ δ; (2) S 1 , where x k < −δ; and (3) S 2 , where x k > δ. Within each subdomain, the system operates as a linear system, so the system is actually governed by two linear subsystems: (1) the underlying linear system (ULS), when the state x is in S 0 , and (2) the overlying linear system (OLS), when x belongs to S 1 ∪ S 2 . Figure 1 illustrates freeplay and the three subdomains S 1,2,3 , separated by two freeplay boundaries Σ 1,2 in state space. 2 S , where k x   . Within each subdomain, the system operates as a linear system, so the system is actually governed by two linear subsystems: (1) the underlying linear system (ULS), when the state x is in 0 S , and (2) the overlying linear system (OLS), when x belongs to 1 2 S S  . Figure 1 illustrates freeplay and the three subdomains 1,2,3 S , separated by two freeplay boundaries 1,2  in state space.

Time Integrations Based on the Runge-Kutta and Hénon Methods
Time integration using a combined fourth-fifth-order Runge-Kutta and the Hénon method (Hénon-RK45) was proposed by Hénon [7] and is widely adopted when solving a set of state space equations with non-smooth nonlinearities. For simplicity, the state space Equation (1)  x to the about-crossing boundary precisely via the following steps: (1) Change the integrating variable from time t to the freeplay state k x ; that is, Thus, a new set of state space equations can be constructed as below:

Time Integrations Based on the Runge-Kutta and Hénon Methods
Time integration using a combined fourth-fifth-order Runge-Kutta and the Hénon method (Hénon-RK45) was proposed by Hénon [7] and is widely adopted when solving a set of state space equations with non-smooth nonlinearities. For simplicity, the state space Equation (1) are expressed as . x = f(x). The Hénon-RK45 method integrates the equations using an ordinary fourth-fifth-order Runge-Kutta scheme (RK45) initially until any freeplay boundary is crossed by x k (t). Assume that the latest time instance before crossing is t break , and the states at that time are x break . Then, the Hénon method moves the freeplay state x k to the about-crossing boundary precisely via the following steps: (1) Change the integrating variable from time t to the freeplay state x k ; that is, Thus, a new set of state space equations can be constructed as below: where y = [x T , t] T is the new vector of the states. (2) Integrate Equation (4) using RK45 T as x k is moved from x k,break to the about-crossing boundary, where y becomes y e = [x T e , t e ] T . (3) Resume the integration of x(t) using the RK45 scheme from x e (t e ) until the next crossing event occurs, and then repeat steps (1) to (3) until the system responses are obtained within a sufficiently long time range.

Preliminary Concepts
In this section, four important concepts are introduced before we propose the state space iterating (SSI) scheme: (1) the Poincaré map, (2) the Lorenz map, (3) the attractors, and (4) the basins of attraction. The first two concepts are introduced in Section 3.1.1, while the last two are provided in Section 3.1.2.

Poincaré Map and Lorenz Map
For a given initial condition (I.C.) x 0 , a start time t 0 , and a time range T, Equation (1) determines a solution x(t), where t 0 ≤ t ≤ t 0 + T. Such a solution can be depicted as an orbit Γ in state space, as shown in Figure 2a, where the closed orbit Γ C indicates a periodic solution of the system. The Poincaré map was defined by Poincaré in 1881 [15]. Perko provided good interpretations of the Poincaré map in Section 3.4 of his book [16]. Suppose that Σ is a hyperplane perpendicular to a periodic orbit Γ C at x 0 ; then, for any point x ∈ Σ sufficiently near x 0 , the orbit Γ through x at t = 0 will cross Σ again at another point P(x). The mapping x → P(x) is called a Poincaré map, or the first return map. Figure 2a illustrates the Poincaré map in state space.
If we could find a fixed point x L on the Poincaré map-that is, x L − P(x L ) = 0-then a periodic solution would be found at x L and an LCO could be determined. However, only a few researchers have tried to find the fixed points of a Poincaré map. An example can be found in Monfared's paper [17], in which he derived the expression of a Poincaré map P(x) numerically and calculated the fixed points of P(x). Nevertheless, the results presented by Monfared were dependent on an artificial parameter λ that has no physical meaning, as the author stressed. The sophisticated procedures for finding fixed points on a Poincaré map could be difficult to apply to a general aeroelastic system.
Other applications of Poincaré maps are actually based on phase diagrams, upon which a set of points at the Poincaré section Σ are projected. This phase diagram consists of many discrete points and is usually simply called a "Poincaré section". Figure 2b shows a Poincaré section derived by Dimitriadis [13], based on which the author explained the mechanism of the associated bifurcations and confirmed that the system was likely undergoing chaotic motion. Many examples can be found in Figure 12 in [10], Figure 10 in [18], Figure 13 in [19], and Figure 11 in [20]. The Poincaré section provides evidence for the type of nonlinear response (e.g., LCO, quasi-periodic motion, or chaotic motion). However, it cannot aid in time integration for finding LCOs. Instead, it usually plays the role of confirming what has been derived from time integration. The Lorenz map was proposed by Lorenz [21] while studying his famous three-DoF system from a simplified model of convection rolls in the atmosphere. An excellent introduction on Lorenz's system and the Lorenz map is presented by Strongatz in Section 9 of his book [22]. In Lorenz's system, three system states, x , y , and z , vary with time. To grasp the features of the system, Lorenz focused only on the connection between the nth local maximum n z and the next one,  The Lorenz map was proposed by Lorenz [21] while studying his famous three-DoF system from a simplified model of convection rolls in the atmosphere. An excellent introduction on Lorenz's system and the Lorenz map is presented by Strongatz in Section 9 of his book [22]. In Lorenz's system, three system states, x, y, and z, vary with time. To grasp the features of the system, Lorenz focused only on the connection between the n-th local maximum z n and the next one, z n+1 , as shown in Figure 3. Then, a map z n+1 = L(z n ), known as the Lorenz map, was obtained through numerous combinations of z n and z n+1 , as presented in Figure 4. Both Figures 3 and 4 were originally presented by Strogatz [22].
The Lorenz map was proposed by Lorenz [21] while studying his famous three-DoF system from a simplified model of convection rolls in the atmosphere. An excellent introduction on Lorenz's system and the Lorenz map is presented by Strongatz in Section 9 of his book [22]. In Lorenz's system, three system states, x , y , and z , vary with time. To grasp the features of the system, Lorenz focused only on the connection between the nth local maximum n z and the next one,    The Lorenz map ( Figure 4) is significant because it reflects the nonlinear pattern of ( ) z t . Note that there are almost no "thicknesses" to the graph in Figure 4, so the pattern of ( ) z t is clear. For any given 0 z , we can predict 1 z by 1 z  0 ( ) L z , followed by   Figure 4). The Lorenz map also provides a visual approach to determine the stability of ( ) z t . Lorenz noted that ( ) z t is always unstable since ( ) 1 n L z  can be seen at every point in Figure 4.
Since Figure 4 was obtained at a series of discrete points, the "worst-case" scenario of 0 z could disobey the pattern of L shown in Figure 4. Tucker [23,24] proved that the Lorenz map does accurately reflect the real features of the system and that the Lorenz The Lorenz map (Figure 4) is significant because it reflects the nonlinear pattern of z(t). Note that there are almost no "thicknesses" to the graph in Figure 4, so the pattern of is clear. For any given z 0 , we can predict z 1 by z 1 =L(z 0 ), followed by z 2 = L(z 1 ) . . . ; thus, an iteration procedure can be constructed and allows us to predict z n simply by using a Lorenz map n times [22]: The fixed points of the Lorenz map, z L : z L − L(z L ) = 0, also determine the periodic solutions, that is, the LCOs, and one can find these fixed points via a visual inspection (e.g., finding the intersections of z n+1 = L(z n ) and z n+1 = z n in Figure 4). The Lorenz map also provides a visual approach to determine the stability of z(t). Lorenz noted that z(t) is always unstable since |L(z n )| > 1 can be seen at every point in Figure 4.
Since Figure 4 was obtained at a series of discrete points, the "worst-case" scenario of z 0 could disobey the pattern of L shown in Figure 4. Tucker [23,24] proved that the Lorenz map does accurately reflect the real features of the system and that the Lorenz attractor exists. Stewart [25] and Viana [26] provided some interpretations of Lorenz's work and stressed the significance of attractors in the Lorenz map. These attractors will be introduced in the next section.

Attractor and Basin of Attraction
The attractor is a special set where mapping F shows the attracting features and is helpful for analyzing nonlinear behaviors. For a system with n x states x(t), Strogatz defined an attractor as a set U min ⊂ R n x of the mapping F : x n+1 =F(x n ) as follows (see Section 9.3 in [22]): (1) Any x(t) that starts in U min stays in U min for all time; (2) there exists an open set U containing U min such that if x 0 ∈ U, then the distance from x(t) to U tends toward zero as t → ∞ ; and (3) there is no proper subset of U min that satisfies both (1) and (2).
In summary, a set U min is an attractor of a mapping F : 3. there is no U sub ⊂ U min satisfies 1. and 2.
where U is called the "basin of attraction" for the attractor U min . An LCO is actually a special attractor that contains only one point, a fixed point x L , which allows Note that the mapping F could be a Poincaré map, Lorenz map, or any other carefully designed type of mapping which reflects the periodicity of the system when any LCO occurs.

Basic Ideas of SSI
As explained in Section 3.1.2, one can find an LCO solution for a system by finding the fixed point of a mapping F. Suppose that the system has n L LCOs. Since each LCO x L,j is an attractor and is associated with a basin of attraction U L,j , according to Equation (6), we could define a set X L that collects all LCOs and a set U L that contains all basins of attraction-that is, X L ={x L | x L − F(x L ) = 0} and U L = ∪ n L j=1 U L,j . Our final goal is to find X L and distinguish each separate point x L in X L .
Two basic conceptions thus emerge: (1) If we use the Poincaré map P as F, there is no method for finding fixed points easily; if we use a Lorenz map L as F, although it is easy to find fixed points via a visual inspection of Figure 4, the results highly depend on the accuracy of the local maxima. Since we have already used the Hénon-RK45 method to move system states accurately into discontinuous boundaries, we expect to find a way to use the states that we already have instead of spending additional and considerable time to find local maxima. Therefore, we need to construct a new mapping F instead of using P or L directly.
(2) If the set X L is extremely difficult to obtain, we can initially find the basins of attraction U L because we know that X L ⊆ U L according to Equation (6). Moreover, if U L is still too difficult obtain, then we can try to find an iteration in state space: We can start with any X {j} 0 that could possibly be obtained, confirm that U L ⊇ X {j} 0 ⊇ X L , and then try to contract X {j} 0 by finding another set X {j+1} 0 which would satisfy U L ⊇ in state space needs to be constructed, which may need some additional results from time integration. Finally, when we find X L and the system patterns are clear enough, we can read all x L values easily from a figure which is similar to Figure 4.
Thus, we have two goals: (1) design a mapping F and (2) construct the iteration procedure shown in Equation (8), which will be implemented in Sections 3.3 and 3.4, respectively.

Modification of the Poincaré Map and Iteration Plots
First, an N T -th time return mapping F = F N T is proposed based on the Poincaré map. We specify a Poincaré section Σ at the upper freeplay boundary, Σ : x = δ; then, let the system evolve from a point x 0 (t 0 ) at Σ to another point x * 0 (t * 0 ) at Σ after N T times crossing Σ, where N T is the "crossing number" and can be specified as any positive integer. The orbit of x(t), which represents the evolution, is obtained by time integration with the Hénon-RK45 method and is illustrated in Figure 5a. The proposition of N T is based on the consideration that one may encounter a lengthy transient phase of the system responses due to an inappropriate selection of x 0 and the complexity of the aeroelastic system with freeplay. Figure 5b provides a general sense of the freeplay state x k when the system is evolving from Thus, we have two goals: (1) design a mapping F and (2) construct the iteration procedure shown in Equation (8), which will be implemented in Sections 3.3 and 3.4, respectively.

Modification of the Poincaré Map and Iteration Plots
where T N is the "crossing number" and can be specified as any positive integer. The orbit of ( ) t x , which represents the evolution, is obtained by time integration with the Hénon-RK45 method and is illustrated in Figure 5a. The proposition of T N is based on the consideration that one may encounter a lengthy transient phase of the system responses due to an inappropriate selection of 0 x and the complexity of the aeroelastic system with freeplay. Figure 5b provides a general sense of the freeplay state k x when the system is evolving following the orbit  . Next, we will observe the orbit  from another perspective and focus on the jth state j x . When the integration starts with 0 0 x t x  , while when the integration is terminated at . Inspired by the Lorenz map, we next construct a two-DoF plane, as shown in Figure 6, where point A represents the orbit  by means of its coordinates x and a point on the two-DoF plane. If we select S N I.C. cases, then we will have Next, we will observe the orbit Γ from another perspective and focus on the j-th state x j . When the integration starts with x 0 (t 0 ), we have x j (t 0 ) = x 0,j , while when the integration is terminated at Inspired by the Lorenz map, we next construct a two-DoF plane, as shown in Figure 6, where point A represents the orbit Γ by means of its coordinates (x 0,j , x * 0,j ). This is similar to what Lorenz did in Figure 4, where Lorenz used a point with the coordinates (z n+1 , z n ) to represent a segment of system behaviors. By means of this method, each I.C. case x 0 is associated with x * 0 and a point on the two-DoF plane. If we select N S I.C. cases, then we will have N S points on the plane. When N S is a large number, sufficient points on the plane could illustrate a spatial pattern, indicating system behaviors. Based on this pattern, we can construct the mapping . We call this two-DoF plane an "iteration plot". Note that we need to construct n x − 1 iteration plots (where n x is the number of system states), with each one corresponding to a state except for the freeplay state x k , because x k always equals δ since both x 0 and x * 0 are at Σ.
illustrate a spatial pattern, indicating system behaviors. Based on this pattern, we can construct the mapping x . We call this two-DoF plane an "iteration plot".
Note that we need to construct 1 x n  iteration plots (where x n is the number of system states), with each one corresponding to a state except for the freeplay state k x , because k x always equals  since both 0 x and * 0 x are at  . Figure 6. Illustration of an iteration plot, in which point A represents how a state behaves when the system evolves from the upper freeplay boundary and returns to that boundary after crossing it several times.

LCO Anlaysis Using the SSI Scheme
As explained in Section 3.2, an LCO is a fixed point L x of the mapping x , from which we can find the fixed point simply by visual inspection. Add a line equ l in the iteration plot for j x , which has a unit slope and indicates that equ l : . If a fixed point can be found at all iteration plots, then at those fixed points, all system states will behave as periodic motions, and an LCO is determined. Figure 6. Illustration of an iteration plot, in which point A represents how a state behaves when the system evolves from the upper freeplay boundary and returns to that boundary after crossing it several times.

LCO Anlaysis Using the SSI Scheme
As explained in Section 3.2, an LCO is a fixed point x L of the mapping F N T : , from which we can find the fixed point simply by visual inspection. Add a line l equ in the iteration plot for x j , which has a unit slope and indicates that l equ : x * 0,j = x 0,j , as shown in Figure 7a,b. Assume that N S points on the iteration plot form a perfect line l 0 , which intersects l equ at Point A; then, A provides the j-th component of . If a fixed point can be found at all iteration plots, then at those fixed points, all system states will behave as periodic motions, and an LCO is determined. Stability can also be determined visually. For simplicity, we use F to represent the mapping instead of    is given as the input of the map- Stability can also be determined visually. For simplicity, we use F to represent the mapping instead of F N T . We define the slope F j as We can prove this by linearizing F j near a fixed point. Suppose that the fixed point is at x 0,j = x * 0,j = x. A small perturbation δx 0,j = x 0,j − x is given as the input of the mapping. Then, which means that or As time progresses, the mapping will be repeated over and over again a ), finally yielding: which indicates that the LCO is unstable. Conversely, if F j < 1, the LCO is stable. Figure 7 illustrates situations of both a stable and an unstable LCO. For simplicity, only one l equ representing F j = 1 is plotted. However, to find fixed points may not be easy, since the patterns shown in the iteration plots may be complicated. As shown in Figure 8, the N S points form a range R 0 instead of the l 0 in Figure 7; therefore, there is no obvious fixed point that can be found. Presumably, a fixed point is hidden inside of the intersection line L 0 : L 0 = R 0 ∩ l equ . Here, we introduce an iterating procedure to gradually find the fixed point hidden in L 0 , which is the SSI scheme proposed in this paper. Here, we introduce an iterating procedure to gradually find the fixed point hidden in 0 L , which is the SSI scheme proposed in this paper.
x values still remain leads to 0 . Second, construct iteration plots using x 0 and x * 0 as explained in Section 3.4, where N S points form a region R 0 and intersect with l equ at L 0 , as shown in Figure 8. The key point is that L 0 normally occupies a smaller range X After that, if R 0 no longer has a "thickness" to the graph, as shown in Figures 4 and 7, the L 0 in Figure 8 will degenerate to a fixed point x L , where an LCO solution is found; otherwise, quasiperiodic motion or chaotic motion should be suspected in the set {x|x ∈ L 0 = R 0 ∩ l equ . The present paper only focuses on LCOs, so any other kind of motion will be estimated as an LCO with a perturbation on the orbit. Studies on LCOs together with all other kinds of nonlinear responses are expected in our future works.
Note that when N T = 1, the proposed mapping is a Poincaré map. However, we found that using N T > 1 could accelerate the determination of clear patterns on the iteration plots and make the SSI scheme more robust. We tried using N T = 3, 5, 9, and 11, all of which yielded the same LCO results but required different numbers of SSI iterations to contract the set X {n} 0 .

Numerical Model and Results of the Hénon-RK45 Method
In this section, a plunge-pitch wing section is introduced, as shown in Figure 9. Time integration with the Hénon-RK45 method is used to obtain the nonlinear LCO behaviors when symmetrical and non-preloaded freeplay nonlinearity is added to the pitch DoF of the wing section. The amount of freeplay is 2δ = 0.1 • . The results in this section will serve as the benchmark for the next section, where the SSI scheme is applied to the same numerical model. The parameters of the wing section are listed in Table A1 in Appendix A.

Numerical Model and Results of the Hénon-RK45 Method
In this section, a plunge-pitch wing section is introduced, as shown in Figure 9. Time integration with the Hénon-RK45 method is used to obtain the nonlinear LCO behaviors when symmetrical and non-preloaded freeplay nonlinearity is added to the pitch DoF of the wing section. The amount of freeplay is 2 = 0.1°. The results in this section will serve as the benchmark for the next section, where the SSI scheme is applied to the same numerical model. The parameters of the wing section are listed in Table A1 in Appendix A. Theodorsen's method [27] and can be expressed as  An unsteady aerodynamic model is established in the frequency domain following Theodorsen's method [27] and can be expressed as f a (k) = 1 2 ρU 2 Qq(k), where ρ = 1.225 kg/m 3 is the air density, k is the reduced frequency, U is the flow velocity (m/s), Q is the aerodynamic influence coefficient (AIC) matrix, and q = [h, α] T are the generalized modal coordinates, where h and α denote the plunge (m) and pitch angle (rad), respectively. Then, the rational function approximation (RFA) introduced by Roger [28] and developed by Karpel [29] and the minimal state (MS) method proposed by Sherwood and Karpel [30] are used to approximate the aerodynamic loads in the time domain, f a (t). The equations of motion for the wing section model are as follows: where M S , D S , and K S1 represent the structural mass, damping, and stiffness, respectively; A 0 , A 1 , A 2 , A D , and A E are matrices of the aerodynamic coefficients; and R is a diagonal matrix that contains two lag terms, r 1 = −0.08 and r 2 = −0.60. The accuracy of f a (t) approximating f a (k) in Equation (13) is examined by converting f a (t) in the Laplace domain, where L( f a (t)) = 1 2 ρU 2 Q(s)q(s). Next, let the Laplace operator be s = ik, where i is an imaginary unit. Figure 10    When freeplay nonlinearity is added to the pitch DoF, the linear structurally restoring forces S1 K q will be replaced by S0 K q . Equations of motion for the wing section with freeplay can be expressed in the form of a set of piecewise linear state space equations: When freeplay nonlinearity is added to the pitch DoF, the linear structurally restoring forces K S1 q will be replaced by K S0 q. Equations of motion for the wing section with freeplay can be expressed in the form of a set of piecewise linear state space equations: . x(t) = A lin x(t) + bK α f NL (α) (14) where x = [q T , . q T , r T ] T represents the system states; K α is the linear stiffness of pitch DoF; A lin represents the linear coefficient matrices for the ULS; b is a constant vector; and the eigenvalues of A lin are shown in Figure 11. For a clear comparison, the eigen-values of matrix A NL of the OLS are presented in Figure 12. The parameters involved in Equations (13) and (14)

LCO Results of Time Integrations with the Hénon-RK45 Method
An I.C. case 0

LCO Results of Time Integrations with the Hénon-RK45 Method
An I.C. case 0

LCO Results of Time Integrations with the Hénon-RK45 Method
An I.C. case x 0 = x(t) is specified as a case where only the pitch angle α(t) starts with a non-zero value α 0 , while the rest of the system states are zero-that is, x 0 = [0, α 0 , 0, 0, 0, 0] T . The amount of freeplay is 2δ = 0.1 • . We first carried out time integrations with the Hénon-RK45 method at an airspeed of U = 20 m/s and various α 0 values. The results show that the system behaviors depend highly on α 0 . For example, when the initial non-dimensional pitch angle, α 0 /δ, equals 5.00, a three-domain LCO (3-D LCO) is found in which the phase trajectory of α(t) between 13 and 15 s passes all three subdomains defined in Section 2.1, as illustrated in Figure 13b. However, when α 0 /δ = 1.00, a two-domain LCO (2-D LCO) emerges, and its phase trajectory passes only two subdomains, S 0 and S 2 , as presented in Figure 14. To determine the effect of 0   on LCO behaviors, 100 cases of 0   that vary from 0.1 to 100 subjects with a 1-cos function and 73 cases of airspeed U that increases from 12 to 20 m/s following a piecewise linear function were considered, as illustrated in Figure 15. Therefore, 7300 cases featuring a combination of 0   and U were specified, and the results derived from the time integrations are presented in Figure 16. To determine the effect of 0   on LCO behaviors, 100 cases of 0   that vary from 0.1 to 100 subjects with a 1-cos function and 73 cases of airspeed U that increases from 12 to 20 m/s following a piecewise linear function were considered, as illustrated in Figure 15. Therefore, 7300 cases featuring a combination of 0   and U were specified, and the results derived from the time integrations are presented in Figure 16. To determine the effect of α 0 /δ on LCO behaviors, 100 cases of α 0 /δ that vary from 0.1 to 100 subjects with a 1-cos function and 73 cases of airspeed U that increases from 12 to 20 m/s following a piecewise linear function were considered, as illustrated in Figure 15. Therefore, 7300 cases featuring a combination of α 0 /δ and U were specified, and the results derived from the time integrations are presented in Figure 16.

LCO Results under an Airspeed of 20 m/s
In this section, we seek to find all possible LCO solutions using time integration with the proposed SSI scheme. The SSI scheme is introduced in Section 3. First, an airspeed of The time integrations for the 7300 cases consume a total central processing unit (CPU) time of 4687 s (about 1.3 h); however, only the effect of α 0 is investigated among all system states. All calculations in this paper were carried out via the software package MATLAB r2014 running on a computer with a quad-core and an eight-thread processor (Intel Core i5-10210U CPU @ 1.60 GHz and 2.11 GHz) and 16.0 GB of random access memory. ⊂ R n x −1 is specified as follows:  Figures 18 and 19 provides an interpretation of Figure 18 to visualize the connection between the mapping F in state space and the iteration plot for x 1 . x , is presented in Figures 18 and 19 provides an interpretation of Figure 18 to visualize the connection between the mapping F in state space and the iteration plot for 1

LCO Analysis of Time Integrations with the SSI Scheme
x .
. Similar situations are seen in the other system states. Therefore, we confirm that for any Thus, we know that there must be some fixed points of F hidden inside region X  Figure 18, and find the intersection region L 0 = R 0 ∩ l equ , where region R 0 consists of N S points.
Then, let X      After Iteration-4, clear spatial patterns emerge on each of the five iteration plots, as shown in the rightmost column of Figure 21. For example, we can focus on the iteration plot for 0,4 x , which is located at the third row and fourth column in Figure 21, and replot it in Figure 22. Three fixed points can be easily distinguished: one is at 0,4 x  10.6 deg/s, one is at 0,4 x  0.45 deg/s, and the last one is at 0,4 x  0. The first two fixed After Iteration-4, clear spatial patterns emerge on each of the five iteration plots, as shown in the rightmost column of Figure 21. For example, we can focus on the iteration plot for x 0,4 , which is located at the third row and fourth column in Figure 21, and re-plot it in Figure 22. Three fixed points can be easily distinguished: one is at x 0,4 ≈ 10.6 deg/s, one is at x 0,4 ≈ 0.45 deg/s, and the last one is at x 0,4 = 0. The first two fixed points represent two LCO solutions, while the last fixed point indicates a convergence of the system responses. The amplitudes and frequencies are derived from the time integrations of Iteration-4; then, a plot of the LCO amplitude/LCO frequency versus x 0,4 is given in Figure 23. It is found that the amplitude of the LCO arising at x 0,4 ≈ 10.6 deg/s is larger than 1.0, which indicates that this is a 3-D LCO, and the other LCO is determined as a 2-D LCO, as its relative amplitudes are less than or around 1.0.
As a result, the system under an airspeed of 20 m/s has two LCOs: one is a 3-D LCO that has a non-dimensional amplitude fluctuating within 7.22-7.27 and a frequency of about 4.98 Hz, and the other one is a 2-D LCO whose non-dimensional amplitude is about 0.46 and whose frequency is 4.26 Hz. Moreover, both LCOs are stable since the slopes of the mapping F near the fixed points satisfy |F | < 1, as explained in Section 3.3.    The SSI scheme also confirms that all LCOs are stable with respect to any system state.
The LCO results derived from time integrations with the Hénon-RK45 method and the SSI scheme are compared in Figures 24 and 25. The Hénon-RK45 method uses 7300 I.C.s (see Figure 15), and the LCO results are analyzed based on a simulation time of 15-20 s, which means that the frequency resolution is 0.25 Hz according to Fast Fourier Transformation (FFT), leading to the discontinuity of frequency shown in Figure 25. In the case of airspeed U > 24 m/s, the system response can be a stable LCO, quasi-periodic motion, or chaotic motion, so the points on Figure 24 after U > 24 m/s have a somewhat irregular distribution. However, with respect to the stable 3-D and 2-D LCOs, the results from the SSI scheme agree well with those from the Hénon-RK45 method. 20 s, which means that the frequency resolution is 0.25 Hz according to Fast Fourier Trans formation (FFT), leading to the discontinuity of frequency shown in Figure 25. In the case of airspeed U  24 m/s, the system response can be a stable LCO, quasi-periodic motion or chaotic motion, so the points on Figure 24 after U  24 m/s have a somewhat irregular distribution. However, with respect to the stable 3-D and 2-D LCOs, the results from the SSI scheme agree well with those from the Hénon-RK45 method.  or chaotic motion, so the points on Figure 24 after U  24 m/s have a somewhat ir distribution. However, with respect to the stable 3-D and 2-D LCOs, the results f SSI scheme agree well with those from the Hénon-RK45 method.

Conclusions
This paper proposed a state space iterating scheme (SSI) to find LCO solutions via a visualization procedure. Inspired by the Poincaré map and the Lorenz map, an N T -time return mapping F was defined and visualized as a set of iteration plots, in which the LCO solutions were represented as the fixed points x L of F and were found within a range X 0 that continuously contracts through a series of iterations: X plots are clear enough, the fixed points can be easily distinguished via a visual inspection, and thus, LCOs can be found.
To verify the SSI scheme, a typical plunge-pitch wing section with pitching freeplay was established. Time integrations with a combined fourth-fifth-order Rung-Kutta method and the Hénon method (Hénon-RK45), as well as the proposed SSI scheme, were applied to the model. Both methods determined a three-domain LCO (3-D LCO) and a two-domain LCO (2-D LCO) under a wide range of flow velocities. The SSI scheme confirmed the minimal X {n} 0 at the fourth iteration and ultimately found the two LCOs x L ∈ X {4} 0 simply through a visual inspection of the iteration plots. The LCO results obtained by SSI were found to be well aligned with those of the Hénon-RK45 method, as shown in Figures 24 and 25.
The highlights of the study are as follows: (1) To obtain the LCO results, Hénon-RK45 checked 100 initial pitch angles for each of the 73 airspeed cases, without guaranteeing that any other arbitrary initial conditions (I.C.s) would lead to the same LCO results, while SSI used 1000 arbitrary I.C.s for each of the 18 airspeed cases, uncovered clear spatial patterns of the system responses, and determined the basins of attraction of the system, which makes the LCO results suitable for any I.C.; (2) the spatial patterns depicted by the SSI scheme could be easily utilized in future studies on other nonlinear responses, such as quasiperiodic motions and chaotic motions, which are usually analyzed via Poincaré sections; (3) the SSI scheme can be easily extended to any other kind of structural nonlinearities and applied to an aeroelastic system with higher dimensions.
However, further studies on SSI are needed, including its application to points (2) and (3) above. We must also determine a better method for selecting arbitrary I.C.s and the parameter N T , which is currently part of our ongoing research.
Author Contributions: Conceptualization, data curation, formal analysis, investigation, methodology, software, validation, visualization and writing (original draft) were conducted by X.W. Project administration was conducted by Z.W. Funding acquisition was conducted by C.Y. Writing (review & editing) was conducted by both X.W. 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.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A Parameters of the Wing Section Model
This appendix provides the parameters and matrices of the numeric model in Section 3.1. The structural modeling procedures are similar to those in Edwards' work [31], which deals with a three-DoF typical wing section: