On the investigation of state space reconstruction of nonlinear aeroelastic response time series

Stall-induced aeroelastic motion may present severe non-linear behavior. Mathematical models for predicting such phenomena are still not available for practical applications and they are not enough reliable to capture physical effects. Experimental data can provide suitable information to help the understanding of typical non-linear aeroelastic phenomena. Dynamic systems techniques based on time series analysis can be adequately applied to non-linear aeroelasticity. When experimental data are available, the methods of state space reconstruction have been widely considered. This paper presents the state space reconstruction approach for the characterization of the stall-induced aeroelastic non-linear behavior. A wind tunnel scaled wing model has been tested. The wing model is subjected to different airspeeds and dynamic incidence angle variations. The method of delays is used to identify an embedded attractor in the state space from experimentally acquired aeroelastic response time series. To obtain an estimate of the time delay used in the state space reconstruction from time series, the autocorrelation function analyis is used. For the calculation of the embedding dimension the correlation integral approach is considered. The reconstructed attractors can reveal typical non-linear structures associated, for instance, to chaos or limit cycles.


Introduction
The treatment of aeroelastic phenomena with linear models have provided a reasonable amount of tools to the assessment and analysis of most of the adverse behaviour [1,2].Nonetheless, modern aviation has shown advances that lead to lighter and faster aircraft, thereby increasing the danger for severe aeroelastic behavior.For instance, transonic flight is surrounded by a complex mixture of flows experiencing different speeds, aggravated by shock waves appearance.Aeroelastic phenomena associated to those compressibility effects introduce a great deal of non-linear effects, and can not be predicted with linear models [3].Highly separated flows also lead to complex aeroelastic phenomena that are difficult to model [4].
Non-linear aeroelasticity research has recently become more relevant.Various approaches have been taken to model non-linear aeroelastic behaviour.In some cases the complex unsteady aerodynamics has been resolved with computational fluid dynamics (CFD) methods [3], or by other methodologies to reduce the computational effort [4,5].However, the majority of non-linear aeroelastic models still need to be validated or checked.In this case, experimental data is of great importance, but it is not common to find significant non-linear aeroelastic data available.
Among the possible behaviour that a non-linear system presents one can assess are many equilibrium points, bifurcations, limit cycles, chaos, etc. Bifurcations and limit cycles occur mainly in transonic aeroelastic responses or when introducing non-linear structural dynamics [6].Moreover, chaotic motion almost certainly appears in highly separated flows, such as the case of dynamic stall [7].
The analysis of non-linear dynamical systems can be based on data from either a mathematical model or an experiment.A variety of mathematical tools have been available to explore and analyse such possible non-linear features, which can also be applied to aeroelastic problems [8].Mathematical models for aeroelastic response associated to the dynamic stall behaviour or stall-induced motion are very hard to obtain.In this case, experimental or flight data seems to provide a more suitable basis for non-linear dynamical analysis.Experimental data furnishes a sequence of measurements that corresponds to a time series with the embedded system dynamics.
Typical dynamical systems responses can be assessed by means of reconstructing the state space from time series using the so-called method of delays (MOD).This technique has been shown to be robust enough to characterize non-linear dynamic systems, as well as to analyze chaotic behaviour.The fundamentals of this method have been introduced by Packard [9] and Takens [10].The MOD uses delayed values of the time series to build a new coordinate system.This leads to the reconstruction of the state space for the observed dynamical system and any embedding attractor of interest in the true state space can be reconstructed.The main task of the MOD is then to provide adequate values for the time delays and the so-called embedding dimension, that is, attractor dimension.Several approaches to obtain the MOD parameters have been investigated.
The aim of this paper is to explore techniques from the theory of time series analysis for the investigation of experimentally acquired non-linear aeroelastic responses.The characterization of the non-linear behaviour of stallinduced oscillations of an aeroelastic wing is achieved using the method of delays.An aeroelastic wing model has been constructed and tested in a wind tunnel.The wing model has been mounted on a turntable that allows variations in its incidence angle.Structural deformation is captured by means of strain gages, thereby providing information on the aeroelastic responses.The test cases correspond to aeroelastic response at specific strain gages points due to oscillatory motions of the turntable.A self-sustained oscillatory motion observed when the turntable is left free to move in the flow field, is also considered for analysis.The parameters for applying the method of delays are the time delay and the embedding (attractor) dimension.The autocorrelation function is used to estimate the time delays, while the correlation integral is used to obtain the embedding dimension.Evolutions of the reconstructed state spaces, with respect to flow and motion parameters, are presented and discussed.

State space reconstruction
Consider a dynamical system where x and F are n-dimensional vectors.The embedding theorem attributed to Takens and Ma ñé [10,11], established that if one is able to observe a single scalar quantity, say h(•), of some function g(x(k)) then the geometric structure of the system dynamics can be unfolded from this set of scalar measurements h(g(x(k)) in a space made out of new vector which define motion in a d-dimensional space [12].If d is large enough, for smooth functions h(•) and g(•) it is shown that many important properties of x(k) are reproduced in the new space given by y k without ambiguity [13].
Let s(k) denote the actual measured variable.If one choose h(•)=s(k) and g i (x(k)) = x(k + T i ) = x(t 0 + (k + T i )τ s ) with τ s being the sampling time, the new vector for T i = iT takes the form The space constructed by using the vectors y k is called the reconstructed space, the parameter T is called time delay and d the embedded dimension.According to the theory of state space reconstruction, the geometric structure of an attractor can be observed in the d-dimensional reconstructed space if d 2d a + 1, with d a the dimension of the attractor of interest.The central issue in the reconstruction of the state space is the choice of the time lag T τ s along with the dimension d.
The time lag T τ s is usually chosen as the quarter of the period of the predominant frequency of the measured variables in the Fourier spectrum or equivalently T is found as the first zero of the linear autocorrelation function where s(k) with N 0 being the total number of sampled points.
There are different methodologies to estimative the embedding dimension.When dynamical system responses are obtained from experiments, noise contamination is practically inevitable.The determination of the MOD parameters must follow specific procedures in order to guarantee proper state space reconstruction.The methodology considered here uses the saturation of system invariants, that is, the invariance properties of an attractor calculated from the reconstructed trajectory does not change by increasing d.To estimate d, the average fraction of the number of points on the attractor with interdistances less than r are calculated from the correlation integral C(r) [8] where and H(t) is the Heaviside function, that is The |.| is taken as the Euclidean distance The correlation integral C d (r) is a function of r and the embedding dimension d.The slope of log 10 C d (r) versus log 10 r is calculated as a function of d over a sufficient range for small interdistances r and the embedding dimension d is thus obtained when the slope becomes independent of d.The slope tends to saturate into a value called the correlation dimension.

Lyapunov exponents
The Lyapunov exponents are essentially a measure of the average rates of expansion and contraction of points in trajectories in phase space.They are asymptotic quantities, defined locally in state space, describing the exponential rate at which a perturbation to a trajectory of a system grows or decays with time at a certain location in the state space [8].
Given a continuous dynamical system in an n-dimensional phase space, the Lyapunov exponents can be monitored in terms of the long-term evolution of an infinitesimal n-sphere of inicial conditions.The infinitesimal n-sphere will become an n-ellipsoid due to the locally deforming nature of the flow [14].The i th one-dimensional Lyapunov exponent is them defined in terms of the length of the ellipsoidal principal axis p i (t) where λ i are ordered from largest to smallest.The Lyapunov exponents signs provide a qualitative picture of a system dynamics.In a three-dimensional continuous dissipative dynamical system the only possible spectra, and the attractors can be described as: (+, 0, −), a strange attractor; (0, 0, −), a two-torus; (0, −, −), a limite cycle; and (−, −, −), a fixed point [14].

Experimental apparatus and database
The aeroelastic wing comprises a wind tunnel model of an arbitrary straight rectangular semi-span wing.The wind tunnel facility presents a testing chamber with about 2m 2 cross-section area.The maximum flow speed in the testing chamber is 50 m s with turbulence level of 0.3%.The wing model is fixed to a turntable that allows incidence variation to the wing.The wing semi-span is 800 mm and the chord is 290 mm.
The model main structure is constructed using fiber glass and epoxy resin in the shape of a tapered plate.The taper ratio is of 1:1.67,where the width at the wing root is 250 mm.To provide aerodynamic shape high density foam and wooden cover are used.The NACA0012 airfoil from wing's root to tip is used.In order to minimize as much as possible the effects of the skin to the wing structure stiffness, both foam and wooden shell are segmented at each 100 mm spanwise.
Figure 1 illustrates the experimental apparatus with indications of the strain gages locations inside the wing model.Incidence motion is achieved with an electrical motor mounted beneath the turntable.The motor actions are controlled by software integrated to the acquisition system.Strain gages are fixed to the plate surface to furnish proper measurement of the dynamic response of the wing main structure.The strain gages are distributed along three lines spanwise.The first and last lines present three strain gages each, all to capture bending motions.The intermediate line presents three strain gages for torsional motion.
Data acquisition and the motion control of the servo motor are achieved by using a dSPACE DS1103 PPC controller board and real-time interface for SIMULINK .The HBM KWS 3073 amplifier for strain gage bridge energizing is used to acquire and amplify the strain gage signals.The resulting signals are directly acquired by the dSPACE controller board, allowing subsequent data storage into a PC compatible.

Results and discussion
During experiments different turntable motions are executed and the respective aeroelastic responses at each strain gage point are acquired.Prescribed turntable motions correspond to oscillatory and random ones, carried out at different airspeeds (from 9 to 16.5 m s , approximately).Oscillatory motions are run at relatively low amplitude values (maximum of 5.5 • ), but such oscillations have been considered around two average angles, that is, zero and 9.5 • .For the cases where the average oscillatory angle is about 9.5 • , highly unsteady separated flow occurs.These cases provide an adequate database for non-linear phenomena investigation.Figure 2 presents a typical aeroelastic  time-history response measured during an oscillatory test case.This test case is obtained with airspeed 16.89 m s , reduced frequency k = 0.1374, average angle of incidence α m = 9.5 • , and motion amplitude equals to 5.5 • .One can observe the existance of complex aeroelastic response that can be infered as the result of highly separated flow.
The randomly generated motions follow the same strategy.Tests are also proceeded for static turntable at different incidences and for free turntable at a range of flow speeds.For the last one, a peculiar self-sustained oscillatory motion at higher angles of attack is observed.The aeroelastic responses associated to this free turntable condition is also examined in this paper.
The techniques for assessing the time delay and embedding dimension are used to provide the basic parameters for state space reconstruction.Figure 3 shows the time delay estimation by means of the autocorrelation function.The time delay is taking where the autocorrelation function becomes zero.Figure 4 illustrates the embedding dimension estimation via the analysis of the correlation integral and its saturation as the dimension increases.Saturation can be observed when the curves of log 10 C d (r) versus log 10 r for each dimension present a similarity in their slopes.
State space reconstruction are achieved for two different motion cases.Firstly, aeroelastic responses due to oscillatory turntable motion are considered, and, then, stall-induced aeroelastic responses when the turntable is free to move are also considered.

Oscillatory turntable
Aeroelastic responses associated to the strain gages at points 1 and 4 (cf.Fig. 1) are analysed.These strain gages are used to measure bending and torsional displacements, respectively, thereby allowing to access time series of important components of the wing aeroelastic behaviour.
Bending measurements have been acquired at three different flow speeds and for different oscillatory frequencies.The outcome of time delay determination for those cases are summarized in the Table 1, where U ∞ is the freestream velocity, ω is the oscillatory frequency, and k is the reduced frequency.The resulting time delays show a direct relation to the oscillatory frequencies, but no relation to the freestream velocities.The embedding dimension is  found to be 3, which leads to the reconstructed state space vectors as State space reconstructions for cases when the freestream velocities and reduced frequencies vary are presented in Figs 5 and 6, respectively.Here the evolution of the reconstructions are shown in terms for the attractor R 2 projections y k = [x(k) x(k + T )] T , in order to make it easy to observe the phenomena.Figure 5 presents the case where the reconstructed spaces evolve for encreasing freestream velocities and fixed oscillatory frequency.The turntable oscillatory frequency for this case is of 1.27 Hz (cf.Table 1).Reconstructed state spaces present similar pattern as the velocity encreases, differing only on the displacement amplitudes.However, in Fig. 6 the evolution of the reconstructed state spaces for a fixed freestream velocity shows complex patterns as reduced frequency increases.Such complexity in both shape and amplitude of the orbits may be connected to strong non-linear phenomena such as chaos or bifurcations.
For torsional measurements (from strain gage at point 4 in Fig. 1), the same procedure adopted in the case of bending measurements is considered.Table 2 presents the time delays that are obtained for the wing torsion  x(t+τ) Fig. 5. Reconstructed state spaces of bending measurements with freestream velocity variation (cases 2, 4 and 6 from Table 1).x(t+τ) Fig. 6.Reconstructed state spaces for bending measurements with reduced frequency variation (cases 3 to 5 from Table 1).measurement time series, where U ∞ is the freestream velocity, ω is the oscillatory frequency, and k is the reduced frequency.Here, the time delay values also show a coherence with respect to the oscillatory frequency.
Figures 7 and 8 present the reconstructed state spaces for cases when the freestream velocities and reduced frequencies vary, respectively.The evolution of the reconstructions are shown in terms of the projections y k = [x(k) x(k + T )] T , as in the bending measurement time series analysis.The results have shown strong influence from noisy data, however the structure of each trajectory can be observed.Differently of the bending measurement time series, in these cases both evolution in freestream velocity and reduced frequency have followed a characteristic pattern.The only parameter that effectively increases from the reconstructed spaces is the motion amplitude.

Free turntable
When the turntable is left free, the aerodynamic forces and other flow effects are the responsable for the wing motion.It has been observed a peculiar self-sustained oscillatory motion of the wing, typical of a limit cycle.The oscillatory motion has kept the amplitude confined in between 4.0 • to 14.0 • turntable incidence angle.In these cases the stall induced a deep break to the increasing pitching moment that builds up for the free wing motion immerse into the flow.Both incidence angle and frequencies increases as flow speed also increases.2).x(t+τ) Fig. 9. Evolution with respect to the freestream velocity of reconstructed state spaces for free turntable cases.
The torsional measurement time series acquired from the strain gage at the point 4 (cf.Fig. 1), has been used to reconstruct the state space.The experiments have been carried out for six different freestream velocities and the resulting time delays for these cases are summarized in the Table 3, where U ∞ is the freestream velocity.Similarly to the prescribed oscillatory turntable motion cases, here the embedding dimension is found to be 3, which also leads to the reconstructed vectors as Figure 9 shows the evolution of the reconstructed state spaces in terms of freestream velocity.Periodic motion is evident from the obtained trajectories.The bouncing phenomenon is also observed in these cases, which characterizes the existence of a resonance mode.The effect of stall inducing the breakdown in pitching moment may be the reason for such phenomenon.The reconstructed state spaces depict a closed orbit that corresponds to a limit cycle.

Lyapunov exponents analysis
Now, to analyse the oscillatory behavior the non-negative Lyapunov exponents of the reconstructed attractor for the fre turntable case is obtained.
Figure 10 presents the tridimensional reconstructed space state without low pass filtering and Fig. 11 presents the reconstructed space state with low pass filtering for case 5 in free turntable (cf.Table 3).Table 4 presents the positive Lyapunov exponents obtained.The choice of the time delay and replacement steps is very important to calculate of the Lyapunov exponents, and it is governed by the necessity to avoid catastrophes [15].A check of the stationarity of the exponent estimates with τ is presented in Fig. 12. Figure 13 shows the same results in Fig. 12, but with low pass filtering.
The choice of the replacement steps depend on additional parameters.Accurate exponent calculation therefore requires the consideration of the following interrelated points: the desirability of maximizing evolution times, the tradeoff between minimizing replacement vector size and minimizing the concomitant orientation error, and the manner in which orientation errors can be expected to accumulate.A check on the stationarity of exponent estimates with the replacement steps is presented in Fig. 14, and in Fig. 15 where a low pass filtering has been considered.The non-negative Lyapunov exponent found for the stall induced oscillations indicates the occurrence of chaos.

Concluding remarks
Techniques from the theory of time series analysis, in the context of non-linear systems, is used in this paper to investigate experimentally acquired non-linear aeroelastic responses.An aeroelastic wing model is tested in wind tunnel.Aeroelastic responses are obtained from strain gages outputs.The wing is mounted on a turntable that allows variations to the incidence angle.The aeroelastic time series are related to oscillatory and free turntable motions.The aeroelastic responses studied in this paper are influenced by highly separated flow fields.The method of delays is employed and for the oscillatory turntable motion a complex evolution of the reconstructed state space is   observed when reduced frequency varies for the same freestream velocity.It is possible to infer that bifurcation and chaotic behaviour are related to the aeroelastic system.It has been observed that when the turntable is left free in the aerodynamic flow, a self-sustained oscillation appears.In this case, the existence of the bouncing phenomenon, characterizing a resonance mode, is observed.The reconstructed attractors depict a closed orbits corresponding to irregular movements.

Fig. 12 .
Fig. 12. Evolution of the Lyapunov exponent with respect to time delay.

Fig. 13 .
Fig. 13.Evolution of the Lyapunov exponent with respect to time delay with low pass filtering.

Fig. 14 .
Fig. 14.Evolution of the Lyapunov exponent with respect to replacement steps.

Fig. 15 .
Fig. 15.Evolution of the Lyapunov exponent with respect to replacement steps with low pass filtering.