NOMENCLATURE
- b
-
semi–span [m]
-
$\hat{b}$
-
bi-coherence
- B
-
bi-spectrum
- f
-
frequency [Hz]
- G
-
pitching freeplay nonlinearity term
-
h,
$\dot{h}$
,
$\ddot{h}$
-
plunging displacement [m], velocity [m/s] and acceleration [m/s2]
- I α
-
aerofoil moment of inertia about the elastic axis [kg·m2]
- K α
-
pitching stiffness [N/rad]
- K h
-
plunging stiffness [N/m]
- L
-
time varying aerodynamic lift/unit span[N]
- m
-
aerofoil mass [kg]
- M
-
higher-order spectra number of data segments
- M α
-
time-varying aerodynamic pitching moment about the elastic axis [N·m]
- M ∞
-
freestream Mach number
- N
-
total number of data points
- p ∞
-
freestream static pressure [Pa]
- q ∞
-
freestream dynamic pressure [Pa]
- R
-
correlation coefficient
- S α
-
static unbalance moment about the elastic axis [kg·m]
- t
-
time [s]
- T
-
tri-spectrum
- T ∞
-
freestream static temperature [K]
- V*
-
velocity index
$V^* = V_\infty /b\omega _\alpha \sqrt{\mu }$
- V ∞
-
freestream velocity [m/s]
- y(t)
-
time-series
- Y(f)
-
Fourier transform of y(t)
- Y*(f)
-
complex conjugate of Y(f)
- α,
$\dot{\alpha }$
,
$\ddot{\alpha }$
-
pitching displacement [°], velocity [°/s] and acceleration [°/s2]
- α s
-
pitching freeplay margin [°]
- α0
-
angle-of-attack
- μ
-
structural-to-fluid mass ratio μ = m/πρ∞ b 2
- ρ∞
-
freestream air density [kg/m3]
- τ
-
tri-coherence
- Φ
-
phase
- ω h , ωα
-
plunging and pitching uncoupled natural frequencies [Hz]
- ω h, aero, ωα, aero
-
wetted plunging and pitching un-coupled natural frequencies [Hz]
Greek Symbol
1.0 INTRODUCTION
The development of reliable and efficient modelling for nonlinear aeroelastic phenomena is currently a critical field of study. Augmenting the fidelity of aeroelastic modelling and ensuring that the computational cost of simulation remains feasible has significant benefits in the aeroelastic community, in particular understanding and diagnosing types of nonlinear behaviour which induce nonlinear aeroelastic phenomena, such as Limit Cycle Oscillation (LCO). LCOs are a result of structural and/or aerodynamic nonlinearities within an aeroelastic system(Reference Dowell1). More specifically, an LCO is a bounded periodic or quasi-periodic dynamic phenomenon which, in the context of aeroelasticity, is bound by structural or aerodynamic nonlinear phenomena. The resultant oscillatory (potentially high-amplitude) motion can cause significant structural damage on the airframe including fatigue failure due to LCO or catastrophic failure if divergence is encountered.
Aerodynamic nonlinearities generally result from transonic flow phenomena, i.e. shockwave–boundary-layer interactions and shock-induced separation, high angle-of-attack dynamic stall operations or, within critical regions of the flight regime, transonic buffet. Structural nonlinearity can be either distributed or concentrated. Distributed nonlinearities will affect the entire wing (i.e. nonlinear material properties or stiffening due to large deformation), whereas concentrated nonlinearities are localised to regions of the airframe, e.g. at locations where there are structural loads due to externally mounted pylons with scientific payload and ordnance, freeplay, backlash providing nonlinear stiffness and friction(Reference Dowell1). For a given aeroelastic system exhibiting a nonlinear response, nonlinear system identification techniques are employed with the aim of detecting, characterising and quantifying the type of nonlinearity. Identifying such systems can become an arduous task if the system contains multiple types of nonlinearity, the type of nonlinearity is not distinguishable due to high-amplitude dynamic oscillations or contaminated by a high signal-to-noise ratio. Considering this, there is a requirement for the development of a robust system identification framework capable of detecting and characterising single or multiple types of nonlinearity in complex noisy aeroelastic signals. Furthermore, there is a requirement to improve our comprehension of how different types of nonlinearity drive nonlinear aeroelastic phenomena, in particular, for this study, LCO.
Generally, freeplay is considered to be inherited from loosened mechanical linkages between the main wing and control surfaces, main wing and store link system or within an all-movable horizontal tail. Current design criteria for freeplay-induced LCO are outdated and based on an elementary comprehension of the phenomena(Reference Ni, Hu, Zhao and Dowell2); hence, it is necessary to improve our understanding of freeplay-induced LCO characteristics. Control surface freeplay is the most common form of freeplay and has been found to induce LCO in various civilian- and defense-based assets. However, these occurrences are not well documented in the public domain. Loosened mechanical linkages can also lead to cubic stiffening within a linkage (rather than a freeplay, which is representative of a dead-zone in the control actuation). Authoritative studies on freeplay nonlinearity for typical sections are provided for subsonic flow in References(Reference Price, Lee and Alighanbari3–Reference Tang, Dowell and Virgin5) and transonic flow in References(Reference Kim and Lee6,Reference Dowell, Thomas and Hall7) . In transonic flow regimes, complex unsteady shock motion may occur as a result of the wing dynamic structural motion. The relationship between the motion of the shock and the dynamic structural motion of the wing is generally linear. However, as the motion increases in amplitude it can shift from Tijdeman Type-A to Type-B, shock motion(Reference Tijdeman8), which is characterised by disappearance and re-appearance of a shock wave and is inherently an inviscid phenomenon. This discontinuity of the pressure distribution across the lifting surface can introduce a dynamic nonlinearity(Reference Thomas, Dowell and Hall9). Nonlinear transonic aerodynamic phenomena for typical sections have been studied extensively, with some examples included in Refs 10–14.
Higher-Order Spectra (HOS) analysis is a valuable tool when analysing nonlinear aeroelastic systems. Their superiority when comparing to traditional linear methods, such as the power spectrum, comes from the ability of the higher-order statistics to predict the presence of nonlinearities. The power spectrum is only able to define second-order statistics and, therefore, can only rigorously unveil physics associated with linear processes. HOS methods are a powerful tool in nonlinear aeroelastic analysis as they demonstrate quadratic and cubic phase coupling between frequencies which result from nonlinearities; hence, HOS can be used to identify the presence and form of nonlinearities, and the transition from linear to nonlinear behaviour within an aeroelastic system(Reference Nikias and Petropulu15–Reference Kerschen, Worden, Vakakis and Golinval18).
HOS methods have been utilised to analyse various nonlinear aeroelastic studies including the analysis of wind-tunnel and flight test experiments, and fundamental research(Reference Silva19). Experimental flutter wind-tunnel results of a pitch/plunge apparatus were analysed using HOS in a study conducted by Silva et al(Reference Silva, Stragnac and Muhammad20). The analysis proved beneficial in identifying regions of linear and nonlinear behaviour and also in identifying the transition from linear to nonlinear aeroelastic behaviour. Further, HOS has been applied to the results of a wind-tunnel flutter experiment for the flexible High-Speed Civil Transport (HSCT) semi-span model in a study conducted by Hajj and Silva(Reference Hajj and Silva21). Nonlinearly coupled frequencies within the aerodynamic forces on the model were identified. HOS is applied to F/A-18 flight flutter test data in a study conducted by Silva and Dunn(Reference Silva, Dunn and Muhammad22). It was found that HOS were valuable in detecting nonlinearities within the experimental data. Hajj and Beran(Reference Hajj and Beran23) performed HOS analysis on F-16 flight flutter test data. The analysis demonstrated areas in which nonlinearities were most prominent and relations between nonlinearities and vibration modes of various components. Finally, the Volterra series has been applied in fundamental analytical studies of nonlinear aeroelasticity by Marzocca et al(Reference Marzocca, Librescu and Silva24). Whilst the aforementioned studies involving aircraft wings or full aircraft proved to be beneficial in identifying nonlinear behaviour, the paradox remains that in all cases the studies are conducted via the analysis of flight test or wind-tunnel data, which is expensive in both cost and time. What has not been considered extensively are comprehensive studies which couple HOS methods with nonlinear aeroelastic predictions. The development of such a method would alleviate the requirement of costly experiments when identifying and analysing nonlinear aeroelastic systems.
In the present study, HOS are utilised to analyse the nonlinear aspects of a two-dimensional pitch/plunge aerofoil system in the presence of pitching freeplay and viscous/inviscid aerodynamic nonlinearities, and under various initial setting angles and speeds. The magnitude of pitching freeplay, initial setting angle and speed are varied, investigating the individual and combined effect of the three parameters on the amplitude and frequency of the limit cycle response. Furthermore, the types of nonlinearity present at various conditions are identified and correlations are made with the identified nonlinear form (quadratic or cubic phase coupling). This provides an initial platform by which different types of aerodynamic and structural nonlinearity can be diagnosed in situations where LCO is observed. A high-fidelity numerical model of the Isogai benchmark configuration(25,Reference Isogai26) is considered modified to include freeplay nonlinearities. The methodology applied in the present paper provides a basis for the development of a robust nonlinear aeroelastic analysis capability which, when LCO or other nonlinear aeroelastic phenomena are encountered, is able to diagnose and quantify the type of nonlinearity present.
2.0 METHODS
In this study, two-way Fluid-Structure Interaction (FSI) simulations are performed; this is a numerical technique which combines transient rigid body pitch/plunge aerofoil motion with high-fidelity Computational Fluid Dynamics (CFD).
2.1 Fluid structure interaction
The FSI simulations are computed using the commercial CFD solver ANSYS Fluent 16.2(27). The nonlinear pitch/plunge structural system is embedded within the ANSYS Fluent solver via a user defined function (UDF). Explicit coupling between the aerodynamic and structural models is achieved via UDF to prescribe generalised motion of the aerofoil. The time–step is chosen to be 0.001 s according to a time-step convergence study to be presented. The transient simulations are initialised from a converged steady-state solution.
2.2 Computational fluid Dynamics
The Navier-Stokes equations for transient flowfields are solved via a coupled pressure-based solver with second-order upwind spatial accuracy. The convergence criteria are set to 10−3 for the scaled residuals at each time step. The one-equation Spalart-Allmaras(Reference Spalart and Allmaras28) RANS turbulence model is chosen to model the viscous flowfield. In order to facilitate the movement of the aerofoil, the C-type grid (Fig. 1) is designed in the interest of preserving cell quality. The grid consists of approximately 9,000 cells with a y-plus value of less than 3. A diffusion smoothing method is implemented with a diffusion parameter of 1.5, aiming to preserve cells closest to the wall and apply the dynamic motion to the cells in the interior of the far-field. The diffusion smoothing method is chosen as it generally results in a better-quality mesh when compared to other smoothing methods such as spring-based smoothing(27).
Figure 1. Computation grid.
2.3 Validation
In order to validate the fluid structure interaction methodology which has been implemented in this study, the linear flutter boundary is computed at a range of Mach numbers and compared to the results of Hall et al(Reference Hall, Thomas and Dowell29) and Timme et al(Reference Timme and Badcock30), as presented in Fig. 2.
Figure 2. Validation of the numerical model.
A comparison is made between results from different numerical methods illustrating the instability boundary as a critical flutter speed index V* f versus freestream Mach number M ∞. The decrease in flutter speed index in the deep transonic region, followed by a second stable branch for higher values of the flutter speed index, is typical of the well-known transonic dip and is distinct for each approach. In the present work, the stability boundary is obtained by increasing V* incrementally until the oscillations are observed to no longer converge (this is repeated for each Mach number). The discrepancy in the comparison of the studies illustrated is generally attributed to the treatment of the boundary layer, as the results of Hall et al(Reference Hall, Thomas and Dowell29) are using the inviscid Euler equations, while the study of Timme et al(Reference Timme and Badcock30) employs viscous simulations with a RANS methodology. For the conditions used in this validation study, see Refs 29 and 30.
2.4 Structural model
The NACA 64A010 aerofoil is used with a structural model based on the Isogai case(25,Reference Isogai26) . Nonlinearities are introduced via a freeplay dead-zone and cubic stiffness, which are incorporated within the pitch/plunge system depicted in Fig. 3. The stiffness values are also varied from the Isogai case(25,Reference Isogai26) , as presented in Table 3.
Figure 3. (a) Two-degree-of-freedom NACA64A010 aerofoil geometry and (b) stiffness as a function of the pitching displacement.
The aeroelastic governing equations for the system are define in matrix form as
$$\begin{eqnarray*}
\left[ \begin{array}{cc}m & S_\alpha \\
S_\alpha & I_\alpha\end{array}\right] \left\lbrace \begin{array}{c}\ddot{h} \\
\ddot{\alpha } \end{array}\right\rbrace + \left[ \begin{array}{cc}k_h & 0\\
0 & k_\alpha\end{array}\right] \left\lbrace \begin{array}{c}h \\
{\alpha }\end{array}\right\rbrace + \left[ \begin{array}{c}0 \\
{G(\alpha )}\end{array}\right] = \left\lbrace \begin{array}{c}L \\
{M_{\alpha }} \end{array}\right\rbrace ,
\end{eqnarray*}$$
where h represents plunging motion and α represents pitching motion and K h and K α are the linear stiffness coefficients in plunge and pitch, respectively. S α is the static unbalance moment about the elastic axis per unit span and I α is the cross-section mass moment of inertia about the elastic axis per unit span. Finally, L is the time-varying lift per unit span, and M EA is the time-varying aerodynamic moment about the elastic axis per unit span. G(α) represents the pitching freeplay nonlinearity as depicted in Fig. 3 and can be defined as
$$\begin{eqnarray}
\begin{array}{rc@{\quad}ll}
G(\alpha ) &=& K_{\alpha }\alpha &\quad \textrm {for} \quad \alpha >\alpha _s,\\
G(\alpha ) &=& 0 &\quad \textrm {for} \quad {-}\alpha _s \ge \alpha \ge \alpha _s,\\
G(\alpha ) &=& K_{\alpha }\alpha &\quad \textrm {for} \quad \alpha <\alpha _s \end{array}
\end{eqnarray}$$
Finally, this can be represented by a system of first-order equations and solved via the fourth-order Runge-Kutta method.
2.5 Operating conditions
The present study considers 20 cases in which the pitching freeplay range of motion α s and initial setting angle α0 are varied. This provides understanding of the individual effect of each parameter as well as the effect of the parameters on each other. All simulations are conducted at a fixed freestream Mach number of M ∞ = 0.8 with the two velocity index values, V* = 0.85 and V* = 1.0. This provides the freestream conditions as presented in Table 1. Table 2 summarises the freeplay and initial setting angle variables.
Table 1 Operating conditions
Table 2 Operating conditions
2.6 Configuration
A benchmark case often used for numerical comparison is the Isogai Case A, a two-dimensional pitch/plunge aeroelastic model with a NACA 64A010 aerofoil section(25,Reference Isogai26) . The structural parameters of this configuration are chosen to represent the dynamics of the outer section of a swept-back wing, with centre of gravity, x CG = 0.4, with offset to centre of rotation, x α = −1.8, and radius of gyration about the centre of rotation, r α = 1.865. This places the pivot point forward of the leading edge. The ratio of uncoupled natural frequencies ω h /ωα = 0.63. No structural damping is considered. The remainder of the structural variables are summarised in Table 3.
Table 3 Structural parameters
3.0 HIGHER–ORDER SPECTRA
HOS are an increasingly used and valuable tool in identifying higher-order statistics that can predict the presence of nonlinearities. The present section provides a brief description of the HOS and their workings, for a more detailed description and derivation the reader is referred Refs 15–18. In this respect, HOS methods can unveil phase coupling between frequencies within an aeroelastic system which originates from nonlinearity, hence, HOS can be used to identify the presence of nonlinearity and the transition from linear to nonlinear behaviour.
Bi-spectral and tri-spectral density analysis can be used to identify quadratic and cubic processes respectfully within the aeroelastic system. The bi-spectrum determines nonlinear interactions within the aeroelastic system by estimating third-order moments in the frequency-domain. The bi-spectral density or bi-spectrum can be defined as
$$\begin{equation}
B(f_1, f_2) = \frac{1}{M}\sum _{i=1}^{M}Y_i(f_1)Y_i(f_2)Y_i^*(f_1 + f_2),
\end{equation}$$
where the block length M is the number of data segments to be considered. The bi-spectrum B(f
1, f
2) can be displayed against two frequency variables f
1 and f
2 in a three-dimensional plot. Each point on the plot describes the bi-spectral energy density of the signal at the bi-frequency (f
1, f
2). The bi-spectrum at any bi-frequency (f
1, f
2) measures the level of interaction between the two frequencies f
1 and f
2. The interactions are a result of quadratic phase coupling between the frequency components, and hence the bi-spectrum detects the presence of quadratic nonlinearity. Y(f) is the Fourier transform of the time series
$\mathcal {F}(y(t))$
. To give a better understanding of the nature of the identification of phase coupling which occurs via the bi-spectrum, one can consider the definition of the bi-spectrum as
$$\begin{equation}
Y(f_1)Y(f_2)Y^*(f_1 + f_2) = ||Y(f_1)||||Y(f_2)||||Y^*(f_1 + f_2)||e^{j\Phi (f_1,f_2)},
\end{equation}$$
where ||Y(f 1)||||Y(f 2)||||Y*(f 1 + f 2)|| measures the magnitude of the interaction between f 1 and f 2, and e jΦ(f 1, f 2) measures the phase. It is also worth noting that
$$\begin{equation}
\Phi (f_1,f_2) = \arctan {\frac{Im[B(f_1,f_2)]}{Re[B(f_1,f_2)]}} = \Phi (f_1) + \Phi (f_2) -\Phi (f_1 + f_2)
\end{equation}$$
While the bi-spectrum detects quadratic relationships between frequency components, it does not provide a basis for the comparison of levels of nonlinearity; hence, as a general practice, the bi-spectrum is conveniently normalised to give the bi-coherence which is bounded between 0 and 1, where spikes in the magnitude of the bi-coherence function represent levels of nonlinearity.
The bi-coherence function
$\hat{b}$
is defined as
$$\begin{equation}
\hat{b}(f_1, f_2) = \frac{|\frac{1}{M}\sum _{i=1}^{M}Y_i(f_1)Y_i(f_2)Y_i^*(f_1 + f_2)|}{\frac{1}{M}\sum _{i=1}^{M}|Y_i(f_1)Y_i(f_2)|^2\cdot \frac{1}{M}\sum _{i=1}^{M}|Y_i^*(f_1 + f_2)|^2}
\end{equation}$$
Similarly, the tri-spectrum identifies cubic nonlinearity within the aeroelastic system by estimating fourth-order moments in the frequency domain, the tri-spectrum is formulated in the same manner as the bi-spectrum; however, it takes into consideration phase coupling between three frequency components (as opposed to two), which can be defined as
$$\begin{equation}
T(f_1, f_2, f_3) = \frac{1}{M}\sum _{i=1}^{M}Y_i(f_1)Y_i(f_2)Y_i(f_3)Y_i^*(f_1 + f_2 + f_3)
\end{equation}$$
Here the tri-spectrum at any tri-frequency (f 1, f 2, f 3) measures the level of cubic interaction between the three frequencies f 1, f 2 and f 3 as a result of cubic nonlinearity within the system. The tri-spectrum is also conveniently normalised, producing the tri-coherence, which is defined as
$$\begin{equation}
\tau (f_1, f_2, f_3) = \frac{|\frac{1}{M}\sum _{i=1}^{M}Y_i(f_1)Y_i(f_2)Y_i(f_3)Y_i^*(f_1 + f_2 + f_3)|}{\frac{1}{M}\sum _{i=1}^{M}|Y_i(f_1)Y_i(f_2)Y_i(f_3)|^2\cdot \frac{1}{M}\sum _{i=1}^{M}|Y_i^*(f_1 + f_2 + f_3)|^2}
\end{equation}$$
Essentially, the bi-spectrum and tri-spectrum indicate whether a phase relationship exists at f 1, f 2 and f 1 + f 2 (bi-spectrum) or f 1, f 2, f 3 and f 1 + f 2 + f 3 (tri-spectrum). If no phase relationship exists, the bi- or tri-coherence will be at or near zero, indicating no quadratic or cubic coupling. Conversely, if a phase relationship does exist, and hence quadratic or cubic coupling, the value will be near unity. Values between zero and unity indicate partial quadratic or cubic coupling.
3.1 Time-series and block length
In estimating the HOS, the data is split into blocks which are evaluated individually and averaged. It is essential that particular attention be given to the block length M in comparison to the total data length N. A larger block size will provide a finer resolution, but this comes with greater variance(Reference Subba and Gabr31). It is suggested by Dalle Molle and Hinch(Reference Dalle Molle and Hinich32) that when identifying n th-order cumulants, the block length should be (n – 1)th root of the sample size to have confidence in the spectral estimate. The present work uses a higher-order Blackman window(Reference Oppenheim and Schafer33) to smooth the time-domain and reduce spectral spillage.
Although particular attention must be paid to the number data points N and block length M in estimating the HOS content of a system, the recommendation by Dalle Molle and Hinch(Reference Dalle Molle and Hinich32) is highly impractical in most cases; for example, to estimate the tri-spectrum using a 512-point FFT (block length of 512 points) the required number of data points would be 134,217,728. As this number of data points is formidable, this section presents a convergence study where the number of data points is increased incrementally until the bi-coherence and tri-coherence values converge to a value sufficient for this research, providing a realistic minimum number of data points to be used.
As the objective of this paper is to characterise various types of nonlinearity via the higher-order frequency content of the system, the requirement of the bi-coherence and tri-coherence analysis to provide insight as to whether each nonlinear form has an underlying quadratic, cubic or coupled quadratic/cubic relationship. Hence, the convergence criteria for the minimum required time-series length are that the bi-coherence and tri-coherence values are bounded within a range of ±0.1. To test the time–series length requirements, the single degree of freedom (s-DOF) Helmholtz-Duffing equation(Reference Hickey, Worden, Platten, Wright and Cooper34) is considered
$$\begin{equation}
m\ddot{y} + c\dot{y} + k_{1}y + k_{2}y^2 + k_{3}y^3 = x(t)
\end{equation}$$
The parameters are chosen to be equivalent to those in Ref. 34, where m = 1 kg, c = 2 N/m/s, k
1 = 104 N/m, k
2 = 107 N/m2 and k
3 = 5 · 109 N/m3 resulting in an undamped natural frequency of 15.92 Hz. The input function x(t) is Gaussian white noise, low-pass pass filtered through 100 Hz. The bi-coherence and tri-coherence are estimated from the velocity
$\dot{y}(t)$
where a clear peak can be seen at 15.9 Hz via self-interaction of the natural frequency as presented in Fig. 4.
Figure 4. Bi-coherence and tri-coherence for the Helmholtz-Duffing equation.
By modifying Equation (8) such that (i) just the cubic term is removed and (ii) just the quadratic term is removed, it would be expected that (i) the bi-coherence presents a peak via self-interaction of the natural frequency while the tri-coherence is zero and (ii) the tri-coherence presents a peak via self-interaction of the natural frequency while the bi-coherence is zero. The bi-coherence and tri-coherence estimates in the present research are based on a block length of a 1,024-point FFT, respectively. To examine the convergence of the bi-coherence and tri-coherence the total number of data points N = 2 n are incremented from n = 10 to 20 for the bi-coherence and n = 9 to 22 for the tri-coherence, as presented in Fig. 5.
Figure 5. Convergence of the bi-coherence and tri-coherence with respect to the length of the data set N.
The HOS estimates are highly sensitive to total number of data points, number of blocks (averages) and block length; thus, from a simulation perspective, HOS analysis is limited by the total number of data which are able to be simulated whilst minimising computational cost. Furthermore, from an experimental perspective HOS analysis is limited due to nonlinear phenomena, such as LCO being dangerous for an asset or wind-tunnel model to endure for long periods of time, again limiting the total number of data points which can be collected. To give an indication as to whether the system is linear, quadratic or cubic (in normalised form) whilst minimising computational cost a minimum of 1017 data points have been found to be sufficient for a SDOF MSD system, no less that 1015 can be used as a minimum. This is to give an indication of the bi-/tri-coherence value to within 0.1 where the entire bi-/tri-coherence spectrum is defined between 0 and 1. For more precise HOS estimations, the recommendations of Dalle Molle and Hinich should be considered. When considering real-world applications for this analysis method rigorous convergence studies should be completed prior to gaining full confidence in the HOS estimates.
3.2 Mesh and time-step convergence
The mesh density and time-step Δt convergence studies are conducted using the operating conditions as presented in Table 4. The convergence of mesh densities and time-steps is defined in terms of the correlation coefficient R, which quantifies the statistical relationship between two or more observed data values as a percentage. The nominal mesh density consists of 9,000 cells, fine of 36,000 cells and super-fine of 144,000 cells. As presented in Table 5, the nominal fine and super-fine mesh densities differentiate statistically by less than 3%; therefore, for the sake of computational efficiency, the nominal mesh density was selected. Table 6 presents the time-step convergence. It can be seen that Δt = 0.01 s and 0.005 s have less than 5% statistical relations to the smaller time-steps. By decreasing the time-step, it can be seen that Δt = 0.001 s and 0.0005 s differ statistically by less than 2%. Hence, Δt = 0.001 s is chosen to be a sufficient time-step.
Table 4 Operating conditions for mesh and time-step convergence
Table 5 Correlation coefficient values for the mesh density convergence study
Table 6 Correlation coefficient values for the time-step convergence study
4.0 RESULTS
The results are presented in the form of time-histories, Power Spectral Density (PSD), bi-coherence plots, tri-coherence plots, pressure distributions and Mach number flowfield contours. Examples of the Mach number flowfields for both setting angles can be observed in Fig. 6. It can be seen that the steady-state flow properties differ significantly between the two setting angles. At α0 = 0°, there is an even distribution of sub- and supersonic flow on both the upper and lower surfaces. Whereas for α0 = 2°, a strong shockwave is observed on the upper surface whilst the flow remains subsonic on the lower surface.
Figure 6. Example static Mach number flowfields at (a) α0 = 0° and (b) α0 = 2°.
4.1 Case ID 1–5
Observation of the phase portraits presented in Fig. 7 indicates that the aeroelastic system is in limit cycle. In plunge, the motion can be seen to be quasi-periodic, and in pitch, it tends towards periodicity. Both modes indicate that the amplitude of the oscillations is low, and in the pitching mode, the deflection of the torsional spring past the freeplay margin is low (not indifferent to the aerofoil re-bounding back and forth between the freeplay boundaries). Mapping the lift and pitching moment against pitching spatial displacement provides insight as to the nature of the aerodynamic forces on the aerofoil as presented in Fig. 8. A clear inflection can be observed in the pitch phase plane, pitching moment and lift diagrams just past the point at which the aerofoil induces maximum deflection of the torsional spring, this indicates that there is some oscillation occurring at the peak deflection as opposed to a clean reversal of the pitching velocity. This is likely due to the impact/release from the freeplay boundary and the interaction with the flow. The power spectral densities (linear frequency contributions) in plunge and pitch are presented in Fig. 9. It is clear that all modes oscillate about a driving frequency of approximately 11.23 Hz (the wetted plunging natural frequency ω h, aero) and its super-harmonics. The strong presence of super-harmonics is an expected finding considering the structural freeplay nonlinearity. The difference between the two DOFs here is that the plunge mode is oscillating at 3ω h, aero whilst the pitch mode oscillates about ω h, aero. The super-harmonics are increasing in odd integers according to f = ω h, aero, 3ω h, aero, 5ω h, aero. . .; this is due to the symmetry of the system about α0 = 0° and hence super-harmonics are generated twice per cycle. The even integers ( f = 2ω h, aero, 4ω h, aero, 4ω h, aero. . .) are non-existent at α s = 0.2° and become weak at α s = 1.0° and, hence, are related to the size of the freeplay margin. Finally, close inspection of Fig. 9 indicates a slight decrease in frequency as the freeplay margin increases; this is to be expected as a greater fraction of the cycle is spent in freeplay and only being acted upon by aerodynamic forces.
Figure 7. Phase portraits in (a) plunge and (b) pitch with V* = 0.85 and α0 = 0° (Case ID 1–5). Vertical lines indicate freeplay dead-zone region.
Figure 8. Change in (a) lift and (b) moment with respect to pitch for V* = 0.85 and α0 = 0° (Case ID 1–5). Vertical lines indicate freeplay dead-zone region.
Figure 9. PSD plots for (a) plunge and (b) pitch V* = 0.85 and α0 = 0° (Case ID 1–5).
The bi-coherence estimates at this speed are not shown as in both DoFs the bi-coherence estimates are null. Contrastingly, the tri-coherence contributions are strong as presented in Fig. 10. These plots indicate strong cubic phase coupling in both modes via the tri-interaction of ω h, aero. There is also very weak cubic phase coupling in the plunge mode where the self-interacting ω h, aero couples 3ω h, aero. In both modes, only α s = 1.0° is shown as the others do not differ. The low amplitude of the oscillations (absence of aerodynamic nonlinearity), absence of quadratic phase coupling and strong cubic phase coupling suggests that the freeplay nonlinearity (in isolation from aerodynamic nonlinearity) is a cubic phenomenon; this is to be further investigated in the coming sections.
Figure 10. Tri-coherence estimates in (a), (b) plunge and (c), (d) pitch for V* = 0.85 and α0 = 0° with α s = 1.0° (Case ID 5). Right-hand image represents a 2D planar interpretation of the left-hand figure.
4.2 Case ID 6–10
With a setting angle of α0 = 2° the nature of the flowfield changes as indicated in Fig. 6 such that the static flowfield is characterised by a shock on the upper surface of the aerofoil and subsonic flow on the lower surface. As a result, the system is no longer symmetrical, and hence the phase portraits and spectral content of the system changes significantly. The phase portraits presented in Fig. 11 indicate that at α s = 0.2° and α s = 0.4° the momentum gained when the aerofoil rebounds off the upper freeplay boundary is not great enough for the aerofoil to strike the lower boundary at the trough of the cycle; hence, the trough of the cycle is characterised by aerodynamics. Furthermore, due to the positive setting angle the aerofoil never plunges below the zero-axis. Figure 12 indicates how the aerodynamic forces are changing with the pitching displacement. It can be seen that below α s = 0.8°, the lift and pitching moment are lift are constantly positive. With close examination of the pitching moment cycle, it can be seen that, at α s = 0.2° and α s = 0.4°, the motion is linear, this is due to the linear nature of the torsional spring and the small margin of freeplay dead-zone that aerofoil is exposed to. The change in lift with pitching angle (Fig. 12) indicates that the trajectory for α s = 0.2° and α s = 0.4° is elliptical; however, above this, it moves to being a lemniscate trajectory. This indicates that the lift when the aerofoil motion is positive is greater than when it is negative for approximately half of the cycle and vice versa for the other half; this is due to the aerodynamic phenomenon of Type-B shock such that shock waves are appearing and disappearing on both surfaces in an asymmetrical fashion. Figure 13 presents the pitching and plunging PSD plots. It can be seen that for α s = 0.2° and α s = 0.4°, the system is characterised by ωα, aero, whereas for α s = 0.6°–1.0°, the system is characterised by ω h, aero. In both cases, odd and even super-harmonics are present. The presence of ωα, aero at low freeplay margins can be attributed to the range of motion in pitch taking precedence over the range of motion in plunge (Fig. 11), and hence the pitching mode is activated.
Figure 11. Phase portraits in (a) plunge and (b) pitch with V* = 0.85 and α0 = 2° (Case ID 6–10). Vertical lines indicate freeplay dead-zone region.
Figure 12. Change in (a) lift and (b) moment with respect to pitch for V* = 0.85 and α0 = 2° (Case ID 6–10). Vertical lines indicate freeplay dead-zone region.
Figure 13. PSD plots for (a) plunge and (b) pitch V* = 0.85 and α0 = 2° (Case ID 6–10).
Similar to α0 = 0.0° the bi-coherence estimates in plunge are zero for all setting angles aside from α s = 1.0° where very weak quadratic phase coupling can be observed via the self-interaction of ω h, aero as per Fig. 14. On the other hand, in pitch there is no quadratic phase coupling for α s = 0.2° and α s = 0.4°; however, above this, the quadratic phase coupling is at unity via the self-interaction of ω h, aero. A physical explanation for this finding is provided in Figs 15 and 16, where it can be seen that at α s = 1.0°, the aerodynamics is highly nonlinear, characterised by strong Type-B shock motion (the shockwave on the upper surface disappears and re-appears within the cycle) and moderate-scale shock-induced flow separation on the upper aerofoil surface. The flowfields are similar for α s = 0.6° and 0.8°; however, for α s = 0.2° and α s = 0.4°, no nonlinear aerodynamic phenomena are found to occur. This suggests that the nonlinear aerodynamics can be characterised by quadratic phase coupling.
Figure 14. Bi-coherence estimates in (a) plunge and (b) pitch for V* = 0.85 and α0 = 2° with α s = 1.0° (Case ID 10).
Figure 15. Pressure coefficient distributions at the (a) peak and (b) trough of the cycle for V* = 0.85 and α0 = 2° with α s = 1.0° (Case ID 10).
Figure 16. Mach number distributions at the (a) peak and (b) trough of the cycle for V* = 0.85 and α0 = 2° with α s = 1.0° (Case ID 10) (0 < M ∞ < 1.3).
The cubic phase coupling is negligible in plunge and hence is not presented. On the other hand, in pitch (Fig. 17), the magnitude of the cubic phase coupling can be seen to be dependent on the nature of the oscillations. At α s = 0.2° and α s = 0.4°, it is found that the oscillations are low magnitude and not impacting with both freeplay margins within the cycle (Figs 11 and 12). As a result, the cubic form of the freeplay nonlinearity is not complete and hence difficult to detect, for these freeplay margins a very low-magnitude peak is observed via the tri-interaction of ω h, aero. Above α s = 0.4°, the system is completing a full cycle, rebounding off both freeplay boundaries, and hence strong cubic phase coupling is observed via the tri-interaction of ω h, aero.
Figure 17. Tri-coherence estimates in pitch for V* = 0.85 and α0 = 2° with (a,b) α s = 0.2° and (c,d) α s = 1.0° (Case ID 10). Right-hand image represents a 2D planar interpretation of the left-hand figure.
4.3 Case ID 11–15
By increasing the velocity index to V* = 1.0, the system is well above the linear flutter boundary and ordered limit cycle behaviour is observed. Figure 18 indicates that above α s = 0.4°, the plunge mode contains a symmetrical knot; however, at α s = 0.2° and α s = 0.4°, the knot does not exist, and hence the super-harmonics in the system are subdued. In the pitching mode, it can now be seen that the oscillations extend far beyond the freeplay margins. In plunge, as the freeplay margin increases, it can be seen that the inflection in the dead-zone region increases; however, the pitching inflection in the dead-zone remains unchanged. This suggests that the amplitude of the pitching oscillations is not being driven solely by the freeplay nonlinearity and that the aerodynamics is taking effect. Figure 19 indicates that the change in lift and pitching moment with respect to pitching displacement remains constant for all freeplay margins, providing further evidence to suggest that the aerodynamics is now a key factor in defining the amplitude of the limit cycle. Another interesting finding is that the amplitude of the plunging motion has decreased when comparing to V* = 0.85; however, the amplitude of the lift has increased. This indicates that the aerodynamics is suppressing the plunging mode of the aerofoil and essentially it is rotating about its pitch axis at high amplitude with low amplitude oscillations in plunge.
Figure 18. Phase portraits in (a) plunge and (b) pitch with V* = 1.0 and α0 = 0° (Case ID 11–15). Vertical lines indicate freeplay dead-zone region.
Figure 19. Change in (a) lift and (b) moment with respect to pitch for V* = 1.0 and α0 = 0° (Case ID 11–15). Vertical lines indicate freeplay dead-zone region.
The power spectral density plots presented in Fig. 20 indicate that the driving frequency in both modes has shifted to approximately 15.5–16.5 Hz (depending on the freeplay margin); this is denoted by ωα, aero. Similar to V* = 0.85 the odd super-harmonics are prominent in all modes; however, in plunge, 2ωα, aero is also prominent. Furthermore, in the plunge mode it can be seen that there is a low-magnitude peak at approximately 38 Hz for α s = 0.2° and α s = 0.4°. This is the a result of the modular frequency interaction ωα, aero + 2ω h, aero.
Figure 20. PSD plots for (a) plunge and (b) pitch V* = 1.0 and α0 = 0° (Case ID 11–15).
The bi-coherence estimates are now strong (near unity); as presented in Figs 21 and 22, this coincides with high-amplitude pitching oscillations and, as a result, the development of aerodynamic nonlinearity. As the size of the freeplay margin increases, it can also be seen that the magnitude of the super-harmonic 2ω h, aero increases, as to be expected for a freeplay-type nonlinearity. To give some insight into the type of aerodynamic nonlinearity, Figs 23 and 24 present the pressure distributions and flowfield at the peak and trough of the cycle; this is symmetrical about α0 = 0.0°. It is clear that the Type-B shock motion is strong as the shock disappears and reappears on both surfaces as the aerofoil moves through the cycle. Some separation can be observed; however, this is not large scale, and hence the strong quadratic phase coupling is likely attributed to the Type-B shock motion. This is further investigated for the next case (V* = 1.0, α0 = 2.0°).
Figure 21. Bi-coherence estimates in plunge for V* = 1.0 and α0 = 0° with (a) α s = 0.2°, (b) α s = 0.4°, (c) α s = 0.6°, (d) α s = 0.8°, (e) α s = 1.0° (Case ID 11–15).
Figure 22. Bi-coherence estimates in pitch for V* = 1.0 and α0 = 0° with α s = 1.0° (Case ID 15).
Figure 23. Pressure coefficient distributions at the (a) peak and (b) trough of the cycle for V* = 1.0 and α0 = 0° with α s = 1.0° (Case ID 15).
Figure 24. Mach number distributions at the (a) peak and (b) trough of the cycle for V* = 1.0 and α0 = 0° with α s = 1.0° (Case ID 15) (0 < M ∞ < 1.3).
The tri-coherence estimates are at unity for all cases via the tri-interaction of ωα, aero due to the freeplay. This is expected due to the freeplay nonlinearity inducing cubic phase coupling. As these plots do not differ in either mode or for aerodynamic forces, they are not shown.
4.4 Case ID 16–20
Figure 25 presents the phase portraits for the system at V* = 1.0 with a setting angle of α0 = 2°. It is demonstrated that in plunge a knot begins to develop as the freeplay margin increases; this is evident for α s = 0.8° at h = 0.08 m and for α s = 1.0° at h = 0.068 m. In pitch, the system exhibits typical high-amplitude limit cycle behaviour. The magnitude of the pitching moment (Fig. 26) changes significantly with freeplay amplitude; however, the form appears to remain constant – this is contrary to α0 = 0.0° where the pitch and pitching moment remain unchanged for all freeplay margin values. This indicates that the amplitude of the oscillations at α0 = 2.0° are not limited solely by aerodynamic phenomena, and the size of the freeplay margin is a significant factor in determining the amplitude of the oscillations. The PSD plots in Fig. 27 demonstrate that the system (in both modes) is driven by ωα, aero and both its odd and even super-harmonics. Interestingly, this is the only case at which 2ωα, aero increases as the freeplay margin increases. The bi-coherence estimates presented in Fig. 28 corroborate the observations made in the aforementioned sections with respect to the effect of the nonlinear aerodynamics and quadratic form that the nonlinear aerodynamics is defined by. In plunge, it can be seen that very weak quadratic phase coupling is up to α s = 0.8 which becomes strong at α s = 1.0, this coincides with the presence of Type-B shock motion and large scale separation (Figs 29 and 30). Thus, it is strongly suggested that the presence of nonlinear aerodynamics can be diagnosed by quadratic phase coupling within the system. Although not shown here, similar observations are made in the pitch DOF.
Figure 25. Phase portraits in (a) plunge and (b) pitch with V* = 1.0 and α0 = 2° (Case ID 16–20). Vertical lines indicate freeplay dead-zone region.
Figure 26. Change in (a) lift and (b) moment with respect to pitch for V* = 1.0 and α0 = 2° (Case ID 16–20). Vertical lines indicate freeplay dead-zone region.
Figure 27. PSD plots for (a) plunge and (b) pitch V* = 1.0 and α0 = 2° (Case ID 16–20).
Figure 28. Bi-coherence estimates in plunge for V* = 1.0 and α0 = 2° with (a) α s = 0.2°, (b) α s = 0.4°, (c) α s = 0.6°, (d) α s = 0.8°, (e) α s = 1.0° (Case ID 16–20).
Figure 29. Pressure coefficient distributions at the (a) peak and (b) trough of the cycle for V* = 1.0 and α0 = 2° with α s = 1.0° (Case ID 20).
Figure 30. Mach number distributions at four different evenly spaced stages between the peak and trough of the cycle where (a) represents the peak and (d) represents the trough for V* = 1.0 and α0 = 2° with α s = 1.0° (Case ID 20) (0 < M ∞ < 1.3).
Similar to V* = 0.85 the cubic phase coupling in plunge is negligible and hence is not shown. However, in pitch (Fig. 31), it is found that there is a strong cubic interaction for all freeplay margins via the tri-interaction of ωα, aero. From the perspective of structural health monitoring and nonlinearity diagnosis, this is a promising finding such that the cubic form freeplay nonlinearity can be detected even when the oscillations are characterised by high-amplitude limit cycle behaviour.
Figure 31. Tri-coherence estimates in pitch for V* = 1.0 and α0 = 2° with (a,b) α s = 0.2° and (c,d) α s = 1.0° (Case ID 16 and 20). Right-hand image represents a 2D planar interpretation of the left-hand figure.
5.0 CONCLUDING REMARKS
HOS are utilised to analyse high-fidelity numerical simulations of a two-dimensional pitch/plunge aeroelastic aerofoil system with the inclusion of a pitching freeplay nonlinearity in pitch. The limit cycle condition is examined for various combinations of speed, initial setting angle and freeplay margin. The significant effect of freeplay and viscous/inviscid aerodynamic nonlinearity on the structure of the limit cycle and nonlinear phase coupling is presented. At a lower speed, the response is dominated by the plunging natural frequency and its harmonics, whereas, at a higher speed, the response is dominated by the pitching natural frequency and its harmonics.
It is found that at a lower speed (V* = 0.85), the zero initial setting angle condition is characterised by low-amplitude oscillations in both modes, which rebound back and forth between the freeplay margins. The nonlinear spectra for this case demonstrate zero quadratic phase coupling and strong cubic phase coupling in all modes via the tri-interaction of the plunging natural frequency. Weak cubic phase coupling is also evident via the interaction of super-harmonic at double the plunging natural frequency. No aerodynamic nonlinearity is clearly evident. At the same speed with a 2 degree initial setting, angle oscillations at low freeplay margins fail to make impact with both freeplay margins, but rather, move into a low-amplitude limit cycle which rebound back and forth off the upper margin, no aerodynamic nonlinearity is clearly evident. The bi-coherence here is null in both modes, and the tri-coherence presents weak cubic phase coupling. At higher freeplay margins, the aerofoil moves between both boundaries of the dead-zone and the amplitude of the limit cycle increases significantly in pitch. The magnitude of the nonlinear spectra here are at unity in both modes via the self-interaction of the plunging natural frequency. This coincides with the presence of aerodynamic nonlinearity in the form of strong Type-B shock motion and some low-scale separation on the upper surface.
It is found that at a higher speed (V* = 1.0) with zero initial setting angle, the pitch mode is characterised by high-amplitude symmetrical limit cycle behaviour for all freeplay margins, while in plunge, the motion is characterised by a very low-amplitude limit cycle for all freeplay margins. The flowfield is by strong Type-B shock motion. The nonlinear spectra here are at unity in all modes.
Finally, at the same speed with a 2 degree initial setting angle, the system is characterised by a high-amplitude asymmetrical limit cycle in pitch and a very low-amplitude asymmetrical limit cycle in plunge. A high gradient drop in lift is observed in the plunge mode. The flowfields are characterised by strong Type-B shock motion, large-scale flow separation and dynamic stall for large freeplay margins and just moderate to strong Type-B shock motion for lower freeplay margins. Here the nonlinear spectra are negligible in plunge. However, in pitch, both nonlinear spectra display low-moderate magnitude nonlinear phase coupling for low freeplay margins and high-magnitude for higher freeplay margins.
The results of the present study suggest that freeplay is a cubic phenomenon and can be diagnosed by the presence of cubic phase coupling whilst inviscid and viscous aerodynamic nonlinearities can be diagnosed by the presence of quadratic phase coupling. Interestingly, if a freeplay dead-zone is present, however, asymmetrical aerodynamic loading leads to limit cycle behaviour which is characterised by the aerofoil only rebounding back and forth off one of the freeplay margins, the cubic phase coupling is weak and the freeplay nonlinearity may be difficult to diagnose.
Future work will consider the inclusion of a linearised reduced-order aerodynamic model to further isolate the sources of nonlinearity. Further, this analysis is to be extended to high-fidelity aeroelastic simulations of three-dimensional wing and full aircraft models or experimental/flight test data with known types of nonlinearity. The aim is to develop a rigorous system identification framework to diagnose the type of nonlinearity present when an asset encounters limit cycle or chaotic response.
ACKNOWLEDGEMENTS
The authors are grateful for the financial support provided by the Defence Science Institute (DSI) for: High-Fidelity Modelling of Wing Flutter and Nonlinear Aeroelastic Predictions. WBS: RE-02290. The authors are also appreciative of the ongoing support provided by the Defence Science and Technology Group (DSTG).
