Two-domain and three-domain limit cycles in a typical aeroelastic system with freeplay in pitch

TL;DR

This paper investigates the existence and characteristics of two-domain and three-domain limit cycle oscillations in a typical aeroelastic system with freeplay in pitch, combining theoretical analysis and wind tunnel experiments.

Contribution

It demonstrates the existence of two- and three-domain limit cycles in aeroelastic systems with freeplay through theoretical derivation and experimental validation.

Findings

Two-domain limit cycles can result from grazing bifurcation.

Limit cycle amplitudes and frequencies can be predicted by equivalent linearization.

Experimental results agree with theoretical predictions.

Abstract

Freeplay is a significant source of nonlinearity in aeroelastic systems and is strictly regulated by airworthiness authorities. It splits the phase plane of such systems into three piecewise linear subdomains. Depending on the location of the freeplay, limit cycle oscillations can result that span either two or three of these subdomains. The purpose of this work is to demonstrate the existence of two-domain cycles both theoretically and experimentally. A simple aeroelastic system with pitch, plunge and control deflection degrees of freedom is investigated in the presence of freeplay in pitch. It is shown that two-domain and three-domain cycles can result from a grazing bifurcation and propagate in the decreasing airspeed direction. Close to the bifurcation, the two limit cycle branches interact with each other and aperiodic oscillations ensue. Equivalent linearization is used to derive the conditions of existence of each type of limit cycle and to predict their amplitudes and frequencies. Comparisons with measurements from wind tunnel experiments demonstrate that the theory describes these phenomena with accuracy.

Tables

Table 1. Table 1: Nominal freeplay of the system

Nominal freeplay (deg)	Measured freeplay (deg)	$\delta$ (deg)
1	0.6 - 0.9	0.3 - 0.45
2	1.8 - 2.0	0.9 - 1
3	3.1 - 3.2	1.55 - 1.6
8	7.4 - 7.6	3.7 - 3.8

Table 2. Table 2: Structural parameters of the experimental system

	Wing dimensions
Chord (with flap)	$c$	25.4 cm
Span	$s$	52 cm
Flexural axis	$x_f$	0.25 c
	Flap dimensions
Chord (flap alone)	$-$	6.25 cm
Span	$s_2$	52 cm
Hinge axis	$x_h$	0.75c
	Inertial parameters
Plunge mass	M	2.562 kg
Pitch inertia	$I_{\alpha }$	0.0181 m kg
Control inertia	$I_{\beta }$	${2.6610}^{-4}$ m kg
Pitch static imbalance	S	0.0943 m.kg
Pitch-Flap inertia product	$I_{\alpha ,\beta }$	0.0013 mkg
Flap static imbalance	$S_{\beta }$	0.0084 mkg
	Stiffness parameters
Plunge stiffness	$K_h$	850.7 N/m
Pitch stiffness	$K_{\alpha }$	34 Nm/rad
Flap stiffness	$K_{\beta }$	1.512 Nm/rad
	Modal parameters
Plunge mode frequency	$f_1$	2.9 Hz
Pitch mode frequency	$f_2$	7.1 Hz
Control mode frequency	$f_3$	17.0 Hz
Plunge mode damping	${\zeta }_1$	0.87 %
Pitch mode damping	${\zeta }_2$	1.39%
Flap mode damping	${\zeta }_3$	0.6 %

Equations

$$M_{\alpha}(\alpha)=\left\{\begin{array}[]{cc}K_{\alpha}(\alpha+\delta)&\mbox{if \hskip 8.61108pt}\alpha<-\delta\\ 0&\mbox{if \hskip 8.61108pt}|\alpha|\leq\delta\\ K_{\alpha}(\alpha-\delta)&\mbox{if \hskip 8.61108pt}\alpha>\delta\end{array}\right.$$

$$\dot{\mathbf{x}}={\mathbf{Q}}_1\mathbf{x}+{\mathbf{q}}_n{M}_{\alpha }(\alpha )+{\mathbf{q}}_p{\alpha }_p$$

$${\mathbf{Q}}_1$$

$${\mathbf{q}}_n$$

$${\mathbf{q}}_p$$

$$\mathbf{E}_{1}=\left(\begin{array}[]{ccc}K_{h}&0&0\\ 0&0&0\\ 0&0&K_{\beta}\end{array}\right)$$

$$\dot{\mathbf{x}}=\left\{\begin{array}[]{llc}\mathbf{Q}_{1}\mathbf{x}+\mathbf{q}_{p}\alpha_{p}&\mbox{\hskip 8.61108pt if \hskip 8.61108pt}|\alpha|\leq\delta&\mbox{\hskip 8.61108pt (a)}\\ \mathbf{Q}_{2}\mathbf{x}-\mathbf{q}_{n}K_{\alpha}\mathop{\mathrm{sgn}}(\alpha)\delta+\mathbf{q}_{p}\alpha_{p}&\mbox{\hskip 8.61108pt if \hskip 8.61108pt}|\alpha|>\delta&\mbox{\hskip 8.61108pt (b)}\end{array}\right.$$

$${\mathbf{x}}_{F_1}$$

$${\mathbf{x}}_{F_2}$$

$${\mathbf{x}}_{F_3}$$

$$T_{grav}=\sin (\theta )\cos (\alpha )Sg\approx \sin (\theta )(1-\frac{{\theta }^2}{2})S\approx \sin (\theta )Sg$$

$$\dot{\mathbf{x}}=\left\{\begin{array}[]{ll}\mathbf{Q}_{1}\mathbf{x}+\mathbf{q}_{p}\alpha_{p}+\mathbf{q}_{n}T_{grav}&\mbox{\hskip 8.61108pt if \hskip 8.61108pt}|\alpha|\leq\delta\\ \mathbf{Q}_{2}\mathbf{x}-\mathbf{q}_{n}K_{\alpha}\mathop{\mathrm{sgn}}(\alpha)\delta+\mathbf{q}_{p}\alpha_{p}+\mathbf{q}_{n}T_{grav}&\mbox{\hskip 8.61108pt if \hskip 8.61108pt}|\alpha|>\delta\end{array}\right.$$

$$\alpha (t)=A\sin \omega t+{\alpha }_0$$

$$M_{\alpha }(\alpha )=a_0+a_1\cos \omega t+b_1\sin \omega t$$

$$a_0$$

$$a_1$$

$$b_1$$

$${\sigma }_1$$

$${\sigma }_2$$

$$b_1=\frac{1}{\pi }{\int }_{-\pi +{\sigma }_2}^{-{\sigma }_2}{K}_{\alpha }(A\sin t+\delta )\sin t\mathrm{d}t+\frac{1}{\pi }{\int }_{{\sigma }_1}^{\pi -{\sigma }_1}{K}_{\alpha }(A\sin t-\delta )\sin t\mathrm{d}t$$

$$a_0$$

$$a_1$$

$$b_1$$

$$\sin \omega t=\frac{\alpha -{\alpha }_0}{A}$$

$$M_{\alpha }(\alpha )=a_0+b_1\frac{\alpha -{\alpha }_0}{A}=a_0-\frac{b_1{\alpha }_0}{A}+\frac{b_1}{A}\alpha$$

$$K_{eq}=\frac{K_{\alpha }}{2\pi }\left(2\pi -(\sin 2{\sigma }_1+\sin 2{\sigma }_2)-2({\sigma }_1+{\sigma }_2)\right)$$

$$\dot{\mathbf{x}}={\mathbf{Q}}_{eq}\mathbf{x}+{\mathbf{q}}_n\left(a_0-K_{eq}{\alpha }_0\right)+{\mathbf{q}}_p{a}_p$$

$${\mathbf{x}}_{eq}=-{\mathbf{Q}}_{eq}^{-1}\left({\mathbf{q}}_n\left(a_0-K_{eq}{\alpha }_0\right)+{\mathbf{q}}_p{a}_p\right)$$

$$\begin{array}[]{cc}|\delta-\alpha_{0}|&\leq A\\ |\delta+\alpha_{0}|&\leq A\end{array}$$

$$K_{eq}-\frac{K_{\alpha }}{2\pi }\left(2\pi -(\sin 2{\sigma }_1+\sin 2{\sigma }_2)-2({\sigma }_1+{\sigma }_2)\right)$$

$${\alpha }_{F_{eq}}-{\alpha }_0$$

$$a_0$$

$$a_1$$

$$b_1$$

$$K_{eq}=\frac{K_2}{2}-\frac{K_2}{2\pi }\left(2{\sigma }_1+\sin 2{\sigma }_1\right)$$

$$\mathrm{if}{\alpha }_0+A\ge \delta$$

Full text

††journal: Journal of Fluids and Structuresedrefedreffootnotetext: PhD Candidate, University of Liège, Allée de la Découverte 9, Belgiumgregrefgregreffootnotetext: Associate Professor, University of Liège, Allée de la Découverte 9, Belgium, Senior Member AIAA.gustavorefgustavoreffootnotetext: Aeronautics Engineer, Technology Development Department, Av. Brigadeiro Faria Lima, 2170, São José dos Campos, São Paulo, 12227-901, Brazilearlrefearlreffootnotetext: William Holland Hall Professor, Duke University, Department of Mechanical Engineering and Materials Science, Box 90300, Hudson Hall, Durham, NC, USA, Grade Honorary Fellow AIAA

1 Introduction

Freeplay in actuators and bearings is a significant source of nonlinearity in aeroelastic systems. Airworthiness authorities place very strict limits on the amount of freeplay allowed in aircraft control surfaces (see [1] for example). Assessing the impact of freeplay on aeroelastic responses is therefore an important aspect of nonlinear aeroelastic research.

Numerous experimental and numerical studies have been published in the literature on the aeroelastic behaviour of simple systems with freeplay nonlinearity, paying particular attention to Limit Cycle Oscillation (LCO) phenomena. The effect of freeplay was investigated on rigid wings with pitch and plunge degrees-of-freedom (DOFs) with freeplay in the pitch DOF [2, 3, 4, 5, 6, 7, 8, 9] but also in the plunge DOF [3] or with an external store added [10]. Tang et al. also investigated a similar system with a flap instead of a plunge DOF [11]. A flexible aerodynamic surface attached to a pitch shaft with freeplay was also considered [12, 13, 14] and Lee et al. even added a plunge DOF [15]. Another well known aeroelastic apparatus is the typical aeroelastic section (a rigid wing with pitch, plunge and control DOFs) with freeplay in the control surface [16, 17, 18, 19, 20, 21, 22, 23]. Most of these systems were shown to undergo several different types of LCO, including asymmetric and aperiodic oscillations at airspeeds lower than their linear flutter speed. Nevertheless, all the LCOs orbit the single fixed point lying at the origin.

In this paper, the typical aeroelastic section is investigated but with freeplay in the pitch DOF instead of the control surface. It is shown that this location of the freeplay can cause anti-symmetric fixed points to come into existence and that limit cycles can orbit these fixed points. This phenomenon is interesting because the system is completely symmetric and yet gives rise to asymmetric periodic responses and fixed points. The emergence of symmetric and asymmetric limit cycles is investigated for both symmetric freeplay and freeplay with preload. The phenomena are demonstrated mathematically and experimentally by means of an aeroelastic wind tunnel model.

2 General equations of motion of an aeroelastic system with freeplay and preload

The pitch-plunge-control aeroelastic system is a 2D symmetric flat plate wing with a control surface. The entire wing is suspended by an extension spring with stiffness $K_{h}$ and a rotational spring of stiffness $K_{\alpha}$ from its pitch axis $x_{f}$. These two springs provide restoring forces in the plunge, $h$, and pitch, $\alpha$, DOFs respectively. The control surface deflection angle $\beta$ is an additional DOF, restrained by a rotational spring with stiffness $K_{\beta}$. The control surface hinge lies at $x_{h}$ and the total chord of the wing is denoted by $c$. The complete system is shown in figure 1

It is assumed that there is freeplay in the pitch DOF, such that the restoring force in the corresponding spring is zero while $|\alpha|<\delta$, $2\delta$ being the width of the freeplay region. Figure 2 shows a typical restoring force diagram for freeplay, whereby the stiffness is $K$ if $|\alpha|>\delta$ and zero otherwise. Note that the freeplay region is centred around the origin.

In the case of the pitch-plunge-control wing with freeplay in the pitch DOF, the stiffness outside the freeplay region is given by $K_{\alpha}$, while the stiffness inside the freeplay region is zero. The restoring moment equation is

[TABLE]

where $M_{\alpha}$ is the pitching moment provided by the freeplay spring.

In addition to the freeplay, an aerodynamic preload angle $\alpha_{p}$ is considered. It models the fact that when the wing is perfectly centred in the freeplay region, it is not perfectly aligned with the airflow and vice versa. As a result, the structural moment, $M_{\alpha}$, depends on $\alpha(t)$ while the aerodynamic moment, $M_{aero}$, depends on $\alpha_{tot}=\alpha_{p}+\alpha(t)$, which introduces an asymmetry in the system.

The equations of motion of the system flying with airspeed $U$ in air of density $\rho$ can be developed using linear unsteady attached flow aerodynamic assumptions; a time-domain model can be written by means of Wagner function analysis [16, 24]. The structural displacements are denoted by the vector $\mathbf{y}=[h\mbox{\hskip 8.61108pt}\alpha\mbox{\hskip 8.61108pt}\beta]$ while the six aerodynamic states are denoted by the vector $\mathbf{w}=[w_{1}\mbox{\hskip 8.61108pt}\ldots\mbox{\hskip 8.61108pt}w_{6}]$. Then the complete state vector of the system is given by $\mathbf{x}=[\dot{\mathbf{y}}\mbox{\hskip 8.61108pt}\mathbf{y}\mbox{\hskip 8.61108pt}\mathbf{w}]^{T}$ and has dimensions $12\times 1$. The equations of motion of the system with freeplay and aerodynamic preload in the pitch DOF are given by

[TABLE]

where

[TABLE]

and $\mathbf{E}_{1}$, is the value of the structural stiffness matrix inside the freeplay region $\pm\delta$, given by

[TABLE]

Matrix $\mathbf{C}$ is the structural damping matrix, $\rho U\mathbf{D}$ is the aerodynamic damping matrix, $\rho U^{2}\mathbf{F}$ is the aerodynamic stiffness matrix, $\mathbf{W}$ is the aerodynamic state matrix, $\mathbf{W}_{1}$ and $\mathbf{W}_{2}$ are the aerodynamic state equation matrices, $\mathbf{M}=\mathbf{A}+\rho\mathbf{B}$, $\mathbf{A}$ is the structural mass matrix and $\mathbf{B}$ is the aerodynamic mass matrix. The matrix $\mathbf{P}$ is an aerodynamic stiffness vector that takes into account the effect of the preload angle $\alpha_{p}$ on the loads acting on the system. The notation $\mathbf{I}_{3\times 3}$ denotes a unit matrix of size $3\times 3$. The values of all the matrices are given in Appendix A. Equation 2 can be written as

[TABLE]

where $\mathbf{Q}_{2}\mathbf{x}=\mathbf{Q}_{1}\mathbf{x}+\mathbf{q}_{n}K_{\alpha}\alpha$.

In this work, we will define two linear sub-systems that are relevant to freeplay:

Underlying linear system: the system without structural stiffness that is only valid inside the freeplay region (equation 17(a)).

• 2.

Overlying linear system: the nominal system without freeplay and with full stiffness (equation 17(b) with $\delta$ = 0).

3 Fixed points

The freeplay function of figure 2 splits the phase plane of the system responses into three piecewise linear subdomains, $S_{1}$ for $|\alpha|\leq\delta$, $S_{2}$ for $\alpha>\delta$ and $S_{3}$ for $\alpha<-\delta$. Response trajectories can span one, two or all three of the subdomains. Furthermore, equation 2 has three fixed points given by

[TABLE]

i.e. they depend on the aerodynamic preload $\alpha_{p}$, the freeplay gap $\delta$ and the airspeed $U$. These fixed points do not coexist; only one of them is an attractor at any instance in time, depending on which subdomain the response trajectory lies in. Dividing equations 18 throughout by $\delta$ we obtain non-dimensional fixed points $\bar{\mathbf{x}}_{F_{i}}=\mathbf{x}_{F_{i}}/\delta$ that only depend on the airspeed and the ratio of the aerodynamic preload divided by the freeplay gap, $\alpha_{p}/\delta$.

Figure 3(a) plots the pitch component of the three fixed points for airspeeds between 0 and 40 m/s and for $\alpha_{p}/\delta$ ratios between 0 and 2. It can be seen that, as the aerodynamic preload ratio increases, fixed point $\mathbf{\overline{x}}_{F_{1}}$ moves from 0 to -1 and eventually exits the freeplay region. Similarly, $\mathbf{\overline{x}}_{F_{3}}=1$ for all airspeeds when $\alpha_{p}/\delta=1$ and enters the freeplay region for all $\alpha_{p}/\delta>1$. This means that only $\mathbf{\overline{x}}_{F_{2}}$ exists for $\alpha_{p}/\delta>1$, since $\mathbf{\overline{x}}_{F_{1}}$ and $\mathbf{\overline{x}}_{F_{3}}$ violate the conditions for existence given in equations 18.

As the aerodynamic preload increases, the system bifurcates from a system with three fixed points to a system with 1 fixed point. This bifurcation is known as a boundary-equilibrium bifurcation and occurs when $\alpha_{p}=\delta$ and $\mathbf{x}_{F_{1}}=\mathbf{x}_{F_{3}}=-\delta$ for all airspeeds. The bifurcation can be more easily visualised in the phase plane plot of figure 3(b). The system’s fixed points are plotted in the $\alpha$-$\dot{\alpha}$ plane for $\alpha_{p}=0$ and $U\neq 0$. The arrows denote the motion of the fixed points as $\alpha_{p}$ is increased. As mentioned earlier, the freeplay region divides the phase plane into three subdomains: $S_{1}$ inside the freeplay boundaries at $\pm 1$, $S_{2}$ and $S_{3}$ outside. Each fixed point is only defined inside its respective subdomain. The arrows show that, as $\alpha_{p}$ increases, fixed points $\mathbf{x}_{F_{1}}$ and $\mathbf{x}_{F_{3}}$ collide with the $-\delta$ boundary and disappear. Another boundary-equilibrium bifurcation occurs at $U=0$, where $\mathbf{x}_{F_{2}}=\delta$, $\mathbf{x}_{F_{3}}=-\delta$ for all values of the aerodynamic preload ratio.

The fixed points of systems with piecewise linear stiffness display a transient characteristic. In the present case there are three piecewise linear subdomains and three piecewise linear systems. The fixed point of each system exists and attracts the response trajectory while the latter lies in the corresponding subdomain. There are two types of fixed point:

Fixed point of system in $S_{i}$ that lies in subdomain $S_{i}$.

• 2.

Fixed point of system in $S_{i}$ that lies in another subdomain.

The first type of fixed point can attract static solutions, i.e. the system response can decay towards it (or, if the response trajectory starts on the fixed point it will stay on it forever). The second type of fixed point cannot attract static solutions and therefore is not a fixed point in the classic sense. However, it can still attract dynamic solutions while the response trajectory travels through subdomain $S_{i}$. The present discussion of the boundary-equilibrium bifurcation only concerns static solutions. Indeed, when a fixed point crosses a discontinuity boundary into the wrong subdomain it cannot attract static solutions anymore. It disappears in a static sense but it still affects the response intermittently every time a trajectory enters the relevant subdomain.

4 Two-domain and three-domain limit cycles

It is reasonable to make the assumption that the overlying (i.e. nominal) linear system is flutter-free and divergence-free inside its flight envelope and that any aeroelastic instabilities are due to the freeplay. Furthermore, it will be assumed that the flutter speed of the underlying linear system $U_{F_{1}}$ is lower than that of the overlying linear system, $U_{F_{2}}$. Two types of periodic solution are then possible:

Circles: These exist entirely in the $S_{1}$ subdomain and can only occur at the flutter point of the underlying linear system, i.e. when $U=U_{F_{1}}$.

• 2.

Limit cycles: These must span at least two subdomains as they can only exist if the system response is nonlinear. They can exist at a range of airspeeds.

The circles and limit cycles are related; the circles bifurcate into limit cycles when their amplitude becomes equal to the width of the freeplay boundary, as displayed in figure 4(a). This bifurcation is known as a grazing bifurcation (see for example [25]). Limit cycles that span two domains, i.e. $S_{1}$ and $S_{2}$ or $S_{1}$ and $S_{3}$ are referred to as two-domain cycles. Limit cycles that span all three domains are referred to as three-domain cycles. Figure 4(b) plots both types of limit cycle. In the absence of aerodynamic preload, it can be seen that a three-domain cycle will orbit $\mathbf{x}_{1}$ and $\mathbf{x}_{2,3}$ if they exist. In contrast, a two-domain cycle can only orbit either $\mathbf{x}_{2}$ or $\mathbf{x}_{3}$. It follows that two-domain cycles can only exist if the fixed points $\mathbf{x}_{2,3}$ also exist.

5 Experimental setup

The experimental system consists of a wing with pitch, plunge and control DOFs mounted vertically in the low-speed wind tunnel of Duke University. The two tips of the wing lie very close to the tunnel wall, so that the flow can be considered quasi-2D. For this study, the freeplay is placed in the pitch DOF. The wing is a NACA 0012 with span 52cm and chord 19cm with a flap of chord 6.35cm mounted to the wing’s trailing edge with micro-bearings and a pin. The model is externally similar to the one used by Conner et al. [16], however the support system and internal structures are different so the dynamic behaviours of the two models are different. The restoring force in the control DOF is provided by a thin piano wire glued to the pin and clamped on the main wing. The wing is attached to the mobile support structure by means of a single spar located at $x_{f}=0.25c$ and two bearings placed on the support. Tightening screws allow the user to change the alignment of the wing in the test section to modify the aerodynamic preload angle $\alpha_{p}$. In this paper, angles of approximately 0 and 5 degrees are considered. Another piano wire is used to provide the necessary restoring torque in pitch. Finally, the mobile support is clamped to the wind tunnel by means of two leaf springs that provide restoring force in the plunge DOF. Five configurations of pitch restoring forces are considered. Firstly the wire was clamped without any freeplay in order to study the overlying linear system. Then it was inserted into plates with holes of different dimensions in order to introduce different amounts of freeplay. The nominal freeplay ranges, defined as the nominal angular distance between the two ends of the freeplay range, considered in this study are given in table 1, along with the measured freeplay ranges and the corresponding values of $\delta$. The freeplay range measurements were obtained using a pitch angular sensor, by moving the piano wire manually between the two ends of the freeplay range. This type of measurement is prone to experimental error therefore a range of values is given in the table.

The structural parameters of the system, given in appendix B, were identified using different techniques. The inertia and stiffness of the different components were measured statically while dynamic tests were performed for validation and for damping computation. Modal analysis performed using the least square complex exponential method at several airspeeds on the overlying linear system highlighted the following modal characteristics:

The plunge-dominated mode has a frequency of 2.9 Hz and a damping of 0.87 % at wind-off conditions. Increasing the airspeed has a hardening effect that makes it interact with the pitch mode.

• 2.

The pitch-dominated mode’s wind-off frequency is 7.1 Hz and its damping is equal to 1.39 %. The airflow has a softening effect on this mode that helps the interaction with the plunge mode.

• 3.

The flap-dominated mode lies at 17 Hz at wind-off with a damping of 0.6 %. The airspeed has a softening effect on this mode, however its frequency is too high to allow interactions with any of the other modes in the airspeed range considered here.

A flutter test performed on the overlying linear system showed that hard flutter occurs at 27-28 m/s due to the interaction of the pitch and plunge modes. All the experiments with freeplay are performed at airspeeds between 8 and 20 m/s, which lie far below the flutter speed of the overlying linear system and are therefore safe.

The system is instrumented in pitch using a meas-spec R30D rotation sensor with a sensitivity of 0.125 mV/deg. A meas-spec R30A angular sensor with sensitivity 0.02 mV/deg is used to measure the flap angle because of its smaller size. Finally an ultrasound sensor with a sensitivity of 10 V/m measures the plunge displacement. The airspeed in the wind-tunnel is measured in real time using a hot wire probe. The data from all these instruments are acquired simultaneously on a NI CompactDAQ with a sampling frequency of 1 kHz. All the signals are 20 seconds long and low-pass filtered at 45 Hz in order to remove electrical noise.

6 Experimental results

Figures 6(a) and 6(b) plot the pitch amplitude and frequency of the limit cycles obtained without aerodynamic preload ($\alpha_{p}\approx 0$) for all airspeeds and freeplay cases. The largest three freeplay values (2, 3 and 8 deg) exhibit similar behaviour as all the amplitude and frequency results fall on a single curve for each of the three freeplay values, as already highlighted by numerous previous studies on freeplay. For the smallest freeplay gap (stars), the results are slightly different. In this case, the amplitude ratio and the LCO onset speed are higher than in the other cases and the LCO frequency is slightly lower. These differences are attributed to the nonlinear friction in the bearings and to the geometry of the freeplay gap. Marsden and Price observed a similar effect on a pitch-plunge system with freeplay and bearings in the pitch [7]. In all four cases, the system undergoes a slow and almost linear amplitude increase with airspeed once the oscillations have started. The LCO frequency features a main branch that starts at about 3.3 Hz and increases to up to 4.2 Hz. The points with frequencies under 3 Hz arise from a secondary peak observed in the Fast Fourier Transforms (FFT) of quasi-periodic oscillations.

The amplitude bifurcation diagram presented in fig. 6(a) plots the peak-to-peak amplitude of the measured response signals and misses an important part of the dynamics: the existence of two-domain and three-domain LCOs. Figure 7 plots pitch time history responses of the system with a freeplay gap of 3 deg and without aerodynamic preload. At 11.8 m/s, the system undergoes quasi-periodic limit cycle oscillations. The responses switch from two-domain to three-domain oscillations and vice-versa in a quasi-periodic fashion because the system is attracted by both solutions. Increasing the airspeed to 12.5 m/s leads to a motion dominated by three-domain oscillations with infrequent occurrences of two-domain oscillations. Then at 12.9 m/s, the two-domain oscillations completely disappear and only three-domain quasi-periodic oscillations are observed up to 13.1 m/s where the system undergoes mono-harmonic limit cycles. Figure 7 clearly shows that two-domain cycles become less and less frequent as the airspeed is increased until they totally disappear and only three-domain mono-harmonic LCOs are observed.

The bifurcation diagrams of the system with an aerodynamic preload angle $\alpha_{p}$ of 5 deg are shown in figure 8. Once again, the largest 3 freeplay values lead to similar results while the measurements from the case with the smallest freeplay are slightly different. This time, two-domain LCOs are clearly observed. They appear at around 11m/s and occur continuously, not intermittently as was the case for $\alpha_{p}=0$ deg. Their amplitude smoothly increases until approximately $U=12$ m/s, where they are replaced by three-domain LCOs. The frequency diagram is similar to the one obtained without preload, the main LCO branch undergoes a frequency increase with airspeed while another branch, corresponding to the secondary peak of the FFT of quasi-periodic oscillations, features a decreasing frequency variation with airspeed.

Figure 9 highlights the transition from two-domain LCOs to three-domain LCOs with a freeplay of 8 deg and a preload angle of 5 deg. At 10.7 m/s, a clear two-domain LCO is observed. Its positive peak amplitude is constant in time while its negative peak amplitude varies unpredictably from cycle to cycle. At 11.3 m/s, two responses can be observed, depending on the initial condition. The first one, plotted in figure 9(b), is another two-domain limit cycle with a larger negative amplitude than in the lower-speed case. The second one, depicted in figure 9(c), is a quasi-periodic LCO where all the periods span the three domains. Finally, at 11.8 m/s only mono-harmonic three-domain LCOs are observed.

Even though no clear two-domain LCOs were observed without aerodynamic preload, the system exhibited a smooth transition from a motion dominated (or influenced) by the two-domain LCOs to a purely three-domain LCOs motion where a single dominant harmonic was observed.

7 Mathematical model of the experiment

The primary mathematical model of the experiment was based on equations 17. While care was taken to make the experiment as symmetric as possible (aside from the aerodynamic preload), neither the wing nor its placement in the tunnel were perfectly symmetric. In particular, it was impossible to ensure that the wing was perfectly vertical, which means that gravity played a small role in introducing additional asymmetry. The mathematical model of the experiment represented all the sources of asymmetry as a gravity effect due to imperfect verticality. Therefore, a gravity term was added to equation 17; for a constant roll angle $\theta$, the pitch torque due to gravity is written as

[TABLE]

where $S=(x_{cg}-x_{f})m$ is the static imbalance of the wing and $g$ is the acceleration due to gravity. Assuming small pitch angles and neglecting second order terms yields a constant torque that depends on the roll angle of the setup. The constant moment of force is then multiplied by $\mathbf{q}_{n}$ and added to equation 17 to obtain the full equations of motion of the model.

[TABLE]

The values of the structural parameters of the equations of motion are those given in appendix B. An equivalent mathematical model representing the reference analytical solution was created in MSC/NASTRAN [26] for comparison. Pure plunge, pitch and control modes were enforced following the reasoning of the ZONATECH/ZAERO model given in [27]. The aerodynamic model consists of doublet lattice panels distributed along a very long span wing. Figure 10 plots the variation of the natural frequencies and damping ratios of the overlying linear system with airspeed. Modal estimates from the wind tunnel tests are compared to the predictions of the three models. Clearly the predictions of the models are in good agreement with each other and with the experimental measurements.

8 Bifurcation analysis using equivalent linearisation

According to the previous discussion, both two-domain and three-domain limit cycles will be asymmetric for $\alpha_{p}>0$. In order to apply equivalent linearisation [28] (also known as the describing function or Krylov and Bogoliubov method) we must look for limit cycles of the form

[TABLE]

where $A$ is the amplitude, $\omega$ is the frequency and $\alpha_{0}$ is the centre of the limit cycle. Three types of limit cycle can occur:

Three-domain cycles, spanning $S_{1}$, $S_{2}$ and $S_{3}$

• 2.

Two-domain cycles spanning $S_{1}$ and $S_{2}$

• 3.

Two-domain cycle spanning $S_{1}$ and $S_{3}$

In all cases, the nonlinear restoring moment of equation 1 is approximated as a Fourier series of the form

[TABLE]

where

[TABLE]

For a three-domain cycle, $A$ must be large enough that $\alpha(t)$ spans all three domains. Figure 11 plots the pitch variation in time and the corresponding freeplay load for the case $\alpha_{0}=0.3$, $\delta=0.5$, $A=1$ and $\omega=1$. The integrals of equation 23 must be carried out in the intervals $[-\pi,t_{1}]$, $[t_{1},t_{2}]$, $[t_{2},t_{3}]$, $[t_{3},t_{4}]$ and $[t_{4},\pi]$. These time instances are given by $t_{1}=-\pi+\sin^{-1}(\delta+\alpha_{0})/A$, $t_{2}=-\sin^{-1}(\delta+\alpha_{0})/A$, $t_{3}=\sin^{-1}(\delta-\alpha_{0})/A$, $t_{4}=\pi-\sin^{-1}(\delta-\alpha_{0})/A$.

As an example, after defining

[TABLE]

the equation for $b_{1}$ becomes

[TABLE]

After working out the integrals and repeating for $a_{0}$ and $a_{1}$, equations 23 become

[TABLE]

Looking back at equation 21, it can be re-arranged as

[TABLE]

so that, after setting $a_{1}=0$, equation 22 becomes

[TABLE]

In this expression, there is a constant term and a term proportional to $\alpha$, i.e. a linear stiffness term. An equivalent linear stiffness can be defined as $K_{eq}=b_{1}/A$, or

[TABLE]

Replacing the nonlinear function in equations 17 by the equivalent linear function of equation 22, the complete equivalent linear system for three-domain cycles is obtained as

[TABLE]

for different values of $A$ and $\alpha_{0}$, where $\mathbf{Q}_{eq}\mathbf{x}=\mathbf{Q}_{1}\mathbf{x}+\mathbf{q}_{n}K_{eq}(A,\alpha_{0})\alpha$. Note that the fixed point of the equivalent linearised system is given by

[TABLE]

The equivalent linearised system can only exist if both $\sigma_{1}$ and $\sigma_{2}$ are real. This means that

[TABLE]

simultaneously. If $\alpha_{p}$ is not equal to zero then $\alpha_{0}$ is also not equal to zero and is it not possible to obtain a symmetric response. Under these circumstances, three-domain limit cycles with amplitude $A=\delta$ cannot exist because they would violate one of the conditions of existence.

The case $K_{eq}=0$ reflects the bifurcation condition for the symmetric system, i.e. $\alpha_{p}=0$. A limit cycle with the lowest possible amplitude $A=\delta$ appears at $U_{F_{1}}$, the flutter speed of the underlying linear system. However, for an asymmetric system this limit cycle cannot occur, as explained earlier. Therefore, $K_{eq}=0$ is not the bifurcation condition for asymmetric systems and limit cycles will start appearing at airspeeds higher than $U_{F_{1}}$.

The equivalent linearisation problem consists in determining the values of $A$ and $\alpha_{0}$ that result in periodic solutions at each airspeed value of interest. The solution process is identical to the one developed in [29]. First, the flutter speed of the equivalent linearized system, $U_{F_{eq}}(K_{eq})$ is calculated for all values of $K_{eq}$ from 0 to $K_{\alpha}$. Then, for each $U_{F_{eq}}(K_{eq})$ we calculate the values of $A$ and $\alpha_{0}$ that give the correct value of the equivalent linearised stiffness when substituted into equation 30 and for which $\alpha_{0}$ is equal to the pitch component of the fixed point of the equivalent linearised system. The nonlinear algebraic system

[TABLE]

can be set up and solved using Newton iterations, where $\alpha_{F_{eq}}$ is the pitch component of $\mathbf{x}_{eq}$ from equation 32. The starting point of the limit cycle branch is $K_{eq}=0$, $U=U_{F_{1}}$. Initial guesses are chosen as $A=\delta$ and $\alpha_{0}=0$.

Once a converged limit cycle is obtained, its stability can be investigated by evaluating a new equivalent linearized system at the same airspeed but with slightly higher amplitude, i.e. $A+\delta A$, where $\delta A<<1$. We calculate the new value of $\alpha_{0}$ that corresponds to this amplitude and the resulting $K_{eq}$. Finally, we set up the new equivalent linearised system using equations 31. If all of its eigenvalues have negative real part, then the limit cycle at $A$ is stable. If any of the eigenvalues have a positive real part then the limit cycle is unstable.

Two-domain limit cycles can be approximated using the same equivalent linearisation scheme. The only difference lies in the values of $K_{eq}$ and $\alpha_{0}$, which are different for two-domain cycles. Applying the Fourier series procedure to two-domain cycles gives the following expressions for $a_{0}$, $a_{1}$ and $b_{1}$ [29]

[TABLE]

so that $K_{eq}$ becomes

[TABLE]

Two-domain limit cycles can exist as long as $|\delta-\alpha_{0}|\leq A$, i.e. $\sigma_{1}$ is real. Furthermore, they must only span two domains so that,

[TABLE]

This means that neither the upper nor the lower bound of the cycle can cross a freeplay boundary. For example, a two-domain cycle spanning $S_{1}$ and $S_{2}$ will disappear if either of the two bounds crosses $+\delta$. Note that $\sigma_{1}$ takes values between $-\pi/2\leq\sigma_{1}\leq\pi/2$, while $\delta-\alpha_{0}$ spans $-A\leq\delta-\alpha_{0}\leq A$.

For the right limit $\sigma_{1}=\pi/2$, $\delta-\alpha_{0}=A$, substituting into equations 35 and 36 yields

[TABLE]

while for the left limit $\sigma_{1}=-\pi/2$, $\delta-\alpha_{0}=-A$ we obtain

[TABLE]

In other words, two-domain limit cycles appear when the equivalent stiffness is equal to the stiffness of the underlying linear system and disappear when $K_{eq}$ is equal to the stiffness of the overlying linear system. Again, these conclusions are only true if $\alpha_{p}=0$; in the presence of aerodynamic preload the $K_{eq}=0$ and $K_{eq}=K_{\alpha}$ appearance and disappearance bounds are modified.

9 Comparison of the mathematical predictions and experimental observations

A comparison of the theoretical and experimental limit cycle amplitude variation with airspeed for $\alpha_{p}=0$ is plotted in figure 12. For all the freeplay gaps considered, the mathematical model predicts two-domain LCOs, which were not observed in practice. In contrast, the amplitudes of the three-domain cycles are predicted with satisfactory accuracy. In all four cases, the predicted two-domain limit cycle branch lies in the airspeed range where quasi-periodic motion was observed in the experiment. That is consistent with the fact that the quasi-periodicity is due to the co-existence of nearby limit cycles at the same airspeed. In all four cases, the model predicts limit cycles at airspeeds lower than the experimental LCO onset speed. Again, it is believed that this disparity is due to the nonlinear damping present in the bearings of the experimental system, which suppressed small amplitude oscillations.

The frequency bifurcation diagram for the system with $\alpha_{p}=0$ is depicted in figure 13. The model accurately predicts the three-domain limit cycle frequency for all the freeplay values. However, the predicted two-domain frequencies (i.e. the part of the branch below 3 Hz) are slightly lower than those observed experimentally.

Figure 14 displays the experimental and mathematical pitch amplitude bifurcation diagrams of the system with an aerodynamic preload angle of 5 deg. The model estimates well both the two-domain and the three-domain LCO amplitudes for freeplay gaps of 2 and 3 deg, slightly underestimates the amplitude when the freeplay gap is 1 degree and slightly overestimates the top branch for 8 degree of freeplay gap. In all four cases, the model predicts limit cycles at airspeeds lower than those observed in the experiment. These small amplitude vibrations are once again damped by the nonlinear friction in the bearings. This is especially noticeable in the 1 deg freeplay cases because the gap is so small that all LCOs have low amplitude. It is important to note that, for freeplay gaps of 2, 3 and 8 deg, the model predicts with satisfactory accuracy the amplitudes of both the two-domain and three-domain limit cycles.

The observed and predicted frequency variations with airspeed for $\alpha_{p}=5$ deg are compared in figure 15. In this case, the frequencies of both the two-domain and three-domain limit cycles are predicted with satisfactory accuracy, although the lowest airspeed at which LCOs can occur is still under-predicted.

The analytical and MSC/NASTRAN models gave nearly identical limit cycle amplitude and frequency predictions for all test cases. For the sake of clarity, sample MSC/NASTRAN predictions are plotted in figure 16 but only the estimates of the analytical model are shown in the other figures.

10 Conclusions

This paper demonstrated mathematically and experimentally the co-existence of two-domain and three-domain limit cycles in aeroelastic systems with freeplay in pitch. The two-domain cycles have small amplitudes, occur at low airspeeds and can be quite irregular or quasi-periodic. The mathematical treatment predicts such cycles in cases with and without aerodynamic preload. However, the limit cycles occurring without preload have even smaller amplitudes and are even less periodic; as a consequence, they were observed intermittently in the experimental work. In contrast, two-domain cycles were clearly observed in cases with preload, both experimentally and mathematically.

From a mathematical point of view, the two-domain and three-domain cycles can be seen as lying on the same limit cycle branch. The former are the result of a grazing bifurcation of the latter. When one of the bounds of a three-domain limit cycle grazes a discontinuity boundary and then enters the freeplay range, the cycle is transformed into a two-domain cycle and undergoes a fold bifurcation. The amplitude of the two-domain cycle keeps decreasing until its other bound also grazes a discontinuity boundary and the cycle disappears completely.

Acknowledgements

The authors would like to acknowledge the financial support of the European Union (ERC Starting Grant NoVib 307265 and Erasmus+ Programme).

Appendix A: Aerodynamic and structural matrices of the model

The matrices appearing in equation 15 are given by

[TABLE]

[TABLE]

where $a=x_{f}/b-1$, $b=c/2$, $c_{h}=x_{h}/b-1$, $I_{\alpha\beta}=I_{\beta}+b(c_{h}-a)S_{\beta}$ and the other quantities are given below. The total aerodynamic damping matrix is given by $\mathbf{D}=\mathbf{D}_{1}+\Phi(0)\mathbf{D}_{2}$ where $\Phi(0)=1-\Psi_{1}-\Psi_{2}$, $\Phi(t)=1-\Psi_{1}e^{-\varepsilon_{1}Ut/b}-\Psi_{2}e^{-\varepsilon_{2}Ut/b}$ is Wagner’s function and

[TABLE]

[TABLE]

The total aerodynamic stiffness is given by $\mathbf{F}=\mathbf{F}_{1}+\Phi(0)\mathbf{F}_{2}+\Xi\mathbf{F}_{3}$ where $\Xi=\Psi_{1}\varepsilon_{1}/b+\Psi_{2}\varepsilon_{2}/b$ and

[TABLE]

[TABLE]

[TABLE]

The aerodynamic state influence matrix is given by $\mathbf{W}=[2\pi b\mathbf{W}_{0}\mbox{\hskip 8.61108pt}-2\pi b^{2}(a+1/2)\mathbf{W}_{0}\mbox{\hskip 8.61108pt}b^{2}T_{12}\mathbf{W}_{0}]^{T}$ where

[TABLE]

The $T_{1}$-$T_{14}$ coefficients are defined in Theodorsen [30] and many other classic aeroelasticity texts.

The structural damping matrix is given by ${\bf D}={{\bf V}^{-1}}^{T}{\bf B}_{mod}{\bf V}^{-1}$, where ${\bf V}$ are the eigenvectors of the matrix ${\bf A}^{-1}{\bf E}$ and ${\bf B}_{mod}$ is given by

[TABLE]

In this latest expression, $\bar{m}_{i}$ are the diagonal elements of the matrix ${\bf V}^{T}{\bf A}{\bf V}$, $\omega_{i}$ are the square roots of the eigenvalues of the matrix ${\bf A}^{-1}{\bf E}$.

Finally, the preload vector $P$ is given by

[TABLE]

Appendix B: Characteristics of the experimental system