Responses of a Strongly Forced Mathieu Equation—Part 1: Cyclic Loading

The forced Mathieu equation can model a variety of physical systems. One example is in the dynamics of large wind-turbine blades [1–5], for which the parametric excitation comes from the cyclic variation in effective stiffness while the blade rotates through its upright position with gravitational compression and softening, and downward orientation with gravitational tension and stiffening. The direct excitation to in-plane motion comes from gravitational loading toward the leading edge, and later the trailing edge, while rotating through horizontal positions, and the excitation to out-of-plane deflections comes from the cyclic variation in lift force as the blade rotates through large changes in altitude, and hence changes in wind speed and direction due to wind shear.

Equation (1) is linear and the solution can be expressed as a sum of homogeneous and particular solutions, as q = q_h + q_p, where q_h satisfies $q¨h+2εμq˙h+(ω2+εγcosΩt)qh=0$ and q_p satisfies $q¨p+2εμq˙p+(ω2+εγcosΩt)qp=F(t)$⁠.

The q_h term exhibits the behavior of the standard Mathieu equation. The regular Mathieu equation is well known to have instability regions in the parametric forcing amplitude-frequency parameter space, such that the fixed point at the origin can be stable or unstable depending on system parameters [23–25]. For lightly damped cases, these regions of instability are based near frequency ratios Ω/ω = 2/N, N being a positive integer. The dependence of the instability regions on the strength of damping has also been explored [26], as have transient response characteristics [27]. Through q_h, these properties apply to the system under study in this work.

The q_p part of the solution, to be approximated by perturbation techniques, can have resonance conditions as well as instabilities. Due to linearity, superposition holds for cases when F(t) has multiple terms. Also, if F(t) is scaled by a factor α, then q_p is scaled by α. In this paper, our work will focus on the case of cyclic excitation with the same frequency as the parametric excitation. Separate work, such as with constant loading F₀ or higher harmonics, can combine with results here by superposition.

This system under hard excitation is primed for study of secondary resonances. Simulations have revealed the existence of superharmonic resonances [1,5] and a primary resonance effect as well [5,28]. A first-order perturbation analysis uncovered the superharmonic of order two [5]. A second-order perturbation analysis detected additional resonances [6], and was applied to study the primary resonance parametric amplification in detail [28]. Second-order perturbation analyses have been useful in other systems with parametric excitation [17,19,23,29,30]. Parametric amplification has also been examined in the Mathieu equation with 2:1 forcing [12–14,31–34]. There have been numerous other studies on non-linear variations of forced and unforced Mathieu and parametric systems [4,5,35–48].

In this work, we analyze Eq. (1) under hard forcing by using second-order expansions of the method of multiple scales (MMS) [23,49], with the aim of capturing and characterizing the superharmonic of order three, and evaluating improvements on the description of the superharmonic of order two. It complements a second-order perturbation analysis of primary resonance and parametric amplification under weak excitation [28]. Elements of this work were introduced in Ref. [6,50], and are more thoroughly developed here, with numerical validations.

As a preview of behavior, Fig. 1 shows the simulated frequency sweep of Eq. (2) with $ε=0.1,μ=0.25,ω=1,F1=1$⁠, and γ = 1.0, 1.1, and 2.0. We see resonant activity at Ω = 1/3, 1/2, 1, and 2. The superharmonic of order three is represented by a small peak, which is indistinguishable between γ = 1.0 and 1.1, but is increased with γ = 2. The superharmonic of order two is an order of magnitude larger, and had previously been described by a first-order multiple-scales analysis [50]. The primary resonance, still another order larger, was explored in detail in Ref. [28]. The behavior of the subharmonic is in fact not a steady-state resonance, but a detected instability, which involves growing responses when parameters are within the instability tongue. In this section, we will analyze these secondary resonances with a second-order perturbation analysis, and discuss some details of their characteristics.

Here at first-order, Ω ≈ ω/2 and Ω ≈ 2ω lead to added secular terms. These cases were studied in the first-order expansions in Ref. [5], and are summarized later. There is also the non-resonant case, in which the leading-order response has the form q(t) = a cos (ωt + β) + ΛsinΩt, where a → 0 at steady-state, leaving us with a non-resonant periodic response q(t) = ΛsinΩt which resembles the standard oscillator.

In the case where $Ω=2ω+εσ$⁠, there is a subharmonic resonance condition. The resonant part of the forced response only has the zero steady-state solution. However, the zero solution can destabilize similarly to the unforced Mathieu equation. In other words, the homogeneous solution has the standard Mathieu equation subharmonic instability, and this governs the forced equation as well.

But Ω ≈ ω/3 is not yet revealed as a case that would lead to secular terms. This specific resonance condition can be reached if we start with the non-resonant case of Eq. (6), and carry forth to the second order.

The next three sections address the resonance cases in detail.

The resonance where Ω ≈ ω/2 can be captured and analyzed at the first-order of multiple-scales expansion [5]. Here we summarize this analysis, and then extend it to the second order and evaluate improvements.

The resonant response magnitude is therefore of $O(ε0)$⁠. The leading-order response from Eq. (5) can then be written as q(t) = acos (2ωt − ϕ) + Γsin (Ωt + θ). This particular solution is accompanied by the homogeneous solution (of the standard damped Mathieu equation), which can have a very slender instability wedge in the (ω, γ) space based at Ω = ω/2 which is quickly diminished with small but increasing damping. The forced superharmonic resonance occurs even if q_h is stable.

Numerical simulations show good agreement when γ is not too large. However, as γ increases, there are some discrepancies in the peak value and location. To this end, we continue with a second-order perturbation analysis to see if there are improvements.

Examination of the X and Y solutions, and a_p1 of Eq. (22), shows that the response amplitude a is order one in terms of the bookkeeping parameter $ε$⁠, and we have seen that a_p1 reduces to a_p0 when $ε$ is neglected in Eq. (22). The second-order analysis indicates that the peak occurs at a frequency slightly lower than the resonant frequency of ω/2, in contrast to the first-order expansion result for which σ_p0 = 0. This negative offset of the peak increases as γ increases. In terms of Ω, this is an $O(ε2)$ offset, which in practice may be small and negligible.

The response scales linearly with the external excitation amplitude, F₁ (through Γ in Eqs. (17) and (18)), but is independent of θ, and the peak frequency is independent of F₁ and θ. The amplitude a scales with a linear γ term and an $ε2γ2$ term (after canceling $ε2$ with the denominator). This latter term is neglected in the a_p1 approximation, and may become significant as γ gets large, but in such case the asymptotic approximation may also degrade.

Equation (19) shows that if μ > 0, then D > 0. Thus, if μ ≠ 0, the response amplitude is bounded, since D serves as the denominator of the amplitude a. For the case of primary resonance [28], also studied using a second-order MMS, the zero denominator coincided with the Mathieu instability boundary, such that forced responses became unbounded simultaneously as the Mathieu system destabilized. If the same mechanism occurs at the superharmonic, it is not captured by the second-order analysis. Furthermore, the trace of the matrix in Eq. (15) is negative (μ > 0), indicating that the approximated solution is stable. It is well known that there exists a slender instability wedge based at Ω ≈ ω/2, but is not captured with a second-order multiple scales. We would therefore expect large errors in the asymptotic analysis as the values of γ and μ approach an instability transition.

We can also see from Eqs. (18) and (22) that the resonating harmonic near the peak grows inversely with μ, typical of linear resonators. Parameter μ is not involved in the non-resonant response, so we would expect the resonance peaks to increase from the surrounding behavior as μ decreases, at least within the limits of the asymptotic perturbation analysis. The aforementioned instability wedge penetrates to lower values of γ as μ decreases. Approaching the unquantified subharmonic instability wedge would be one mechanism for deviating from the simple inverse dependence on μ indicated by the perturbation method.

The correction term $εq1(t)$ has harmonics that scale with $εa$⁠, and thus it should remain $O(ε)$ smaller than the q₀ term. However, q₁ contributes harmonics not expressed in q₀, and is seen in numerical comparisons to improve the approximation to q. As this is part of an asymptotic expansion, there could be limits on the benefit of including q₁ if $εγ$ cannot be regarded as “small.”

Figure 2(a) shows amplitudes of frequency response sweeps comparing numerical simulations (solid curves) and second-order MMS approximations (dashed). The horizontal axis is the excitation frequency Ω, and the parameters are $ε=0.1,μ=0.25,F1=1$⁠, and ω = 1, with γ = 1, 2, 3, and 4. The levels of the curves increase with γ. The simulated amplitudes at each frequency were obtained after transients were removed, and computed as the maximum of |q(t)| during a cycle. The analytical approximations were taken from $|q0+εq1|$ over one cycle, obtained from Eq. (36).

Figure 2(b) graphs the amplitude a of the resonant term as γ increases. The solid curve is the second-order MMS prediction of the maximum a from Eqs. (17) and (20) (which is visually indistinguishable from a plot based on Eq. (22)), while the dashed curve is the first-order MMS prediction from Eq. (9). They are nearly the same. However, as can be seen by Fig. 2(a), the resonance peaks occur at frequencies which deviate slightly from Ω = 0.5, as predicted for a in the second-order approximations but not in first-order approximations. While the generation of the curves involves other harmonics in both the simulations and from Eq. (23), we expect a to have the dominant effect.

Also shown in Fig. 2(b) are values of a obtained from fast Fourier transforms (FFTs) of numerical solutions, at the peak frequencies as estimated from the perturbation analysis. In the numerical solutions, the sampling rate was set to have an integer number of samples (100) per cycle to remove leakage distortions on spectrum peaks. Agreement is good up through moderate values of γ. Since the peak frequency from the perturbation analysis was used for the numerical solutions, it is possible that the numerical solution peak is slightly underestimated, although a difference is probably imperceptible when γ is not large. At γ ≥ 7 there is large deviation between the simulations and amplitude curves. At γ = 8 simulations exhibit stable amplitudes well off the plotting scale. By 8.5, unstable responses are observed in simulations. As mentioned, the slender instability tongue at this superharmonic is not uncovered by the second-order perturbation analysis, and thus the analysis does not accommodate the unbounded amplitude for large γ. Interpretations must also heed the assumption that the asymptotic solution is “valid” for small $ε$ and order-one γ, a notion which is pushed with the larger values of γ in these examples.

Assembling the observations, the advantages of performing the second-order perturbation for this resonance are minimal. The quality of the first- and second-order approximations is very similar, and the small deviations that might be seen at larger γ, or in the frequency peaks, are likely to be negligible relative to other sources of error, such as asymptotic divergence, or in modeling or estimated parameter values.

Stability can be ascertained from the Jacobian A of Eqs. (28) and (29). We find that if μ > 0, then det(A) > 0 and tr(A) < 0 and so the response at this resonance is stable.

Based on Eq. (34), the superharmonic component is $O(ε)$ while the non-resonant term is O(1). The peak amplitude scales with γ², and inversely with damping μ. The amplitude is not affected by θ, although there may be a slight peak-to-peak effect by the phasing of harmonics.

Figure 3(a) shows example amplitude frequency sweeps through superharmonic resonances of orders four and three. The figure shows amplitudes defined as the maximum of the absolute value of the response in a cycle of steady-state oscillation at each forcing frequency in the sweep. The solid lines show the numerical simulations, and the dashed lines show the results from $q=q0+εq1$ as expressed in Eqs. (35) and (36). While the resonances of orders three and four (and even five for large enough γ) are observed in numerical simulations, only that of order three is captured by the second-order MMS analysis. The plots show agreement between the numerical and asymptotic solutions if γ is not too large, and indicates increasing error as γ becomes larger.

Order four and order five superharmonics may be captured by extending the MMS analysis to higher orders of $ε$⁠. This can be seen by considering an MMS expansion that is non-resonant until the nth-order. The parametric term q₀cosΩT₀ in the $O(ε)$ equation in Eq. (4) generates exponential terms with frequencies ±2Ω and ±Ω ± ω, which become frequency components of the q₁ solution. At $O(ε2)$⁠, the q₁cosΩT₀ term then generates frequencies ±2Ω ± Ω and ±Ω ± Ω ± ω, the first of which leads to the forced superharmonic resonance Ω ≈ ω/3 analyzed here. (There is also a primary resonance, which is not part of the non-resonant expansion and is analyzed separately with weak excitation [28].) At $O(ε3)$⁠, the q₂ cos ΩT₀ term would generate frequencies of ±2Ω ± Ω ± Ω and ±Ω ± Ω ± Ω ± ω, among which is a potential superharmonic resonance with Ω ≈ ω/4, as well as Ω ≈ 2ω/3 (an unforced Mathieu wedge). Continuing, at $O(εn)$⁠, the q_n-1cosΩT₀ term would generate frequencies of ±2Ω + (n − 1~occurrences of~ ± Ω) and (n~occurrences of~ ± Ω) ± ω, revealing a potential forced superharmonic resonance of Ω ≈ ω/(n + 1), and an unforced Mathieu wedge at Ω ≈ 2ω/n. A full analysis by nth-order MMS would generate analytical predictions of steady-state amplitudes and stability criteria. Evidence of this cascade of superharmonic resonances was present in simulations of the undamped case [1].

Figure 3(b) depicts the amplitude a of the resonant term as γ increases. The solid curve is the second-order MMS prediction from Eq. (34). Also shown in Fig. 3(b) are values of a obtained from FFTs of numerical solutions at the peak frequencies as estimated from the perturbation analysis. In the numerical solutions, an integer number of samples per cycle (100) were chosen to remove leakage distortions on spectrum peaks. Agreement is reasonable up through moderate values of γ. Since the peak frequency from the perturbation analysis was used for the numerical solutions, it is again possible that the numerical solution peak is slightly underestimated, although we expect this effect to be small when γ is not large, since peaks are locally quadratic. The trend suggests that the perturbation analysis slightly overestimates the response amplitudes. When γ gets larger than eight, the simulation value becomes large and off the scale. At about γ = 11, the numerical solution is unstable. It is known that there is a very slender instability tongue near this superharmonic for very small damping, but it is not uncovered by the second-order perturbation analysis, and thus the analysis does not accommodate the unbounded amplitude for large γ. Interpretations must consider that the perturbation study is asymptotic, “valid” for small $ε$ and order-one parameters such as γ, such that $εγ$ is “small.”

Figure 4 shows the amplitude a (log scale) and phase ϕ behavior of the resonant term in the response expressed in Eq. (35). Indeed, based on Eqs. (32) and (31), the negative amplitude with a phase near π suggests that the response is nearly in phase with the excitation near the resonance peak, which is offset slightly below the superharmonic frequency (Ω = ω/3) for sizable γ. The phase change and the peak amplitude occur at decreasing frequencies as γ increases.

While a superharmonic resonance of order 1/3 is apparent in Eq. (25), we note that the analysis does not reveal a subharmonic resonance of Ω ≈ 3ω in this linear forced Mathieu equation.

Figure 1 shows an example of how the instability may be observed in a simulation. For γ = 1 (solid line), the parameters come very close to the instability boundary. The peak seen at Ω = 2 has a value of approximately 1.75. This is not a steady-state resonant amplitude, but a nearly neutrally stable response amplitude essentially defined by the initial conditions. Changing the time duration for removing the transients has very little effect on this value, in this example. When γ = 1.1 (dash-dot line), at Ω = 2 the parameters just inside the instability boundary, and dictate very slow unstable growth. In this simulation, the amplitude has reached a value of about 5.2. Longer simulation times for removal of transients will lead to a higher peak. For the case of γ = 2 (dashed line), the parameters are well inside the unstable zone near Ω = 2. When varying Ω increases into the instability wedge (at about Ω = 1.85) the instability achieves an unstable exponent with quick growth. Thus, while the stable steady-state resonances near Ω = 1/3, 1/2, and 1 show smooth peak transitions from the surrounding non-resonant state, the subharmonic is signified by abrupt peaks departing from the non-resonant curve, which will grow with integration time. If we choose γ = 0.9 (not shown), then Ω = 2 is stable, but very slowly, and there is a shrinking peak as the response decays from the initial state and the transient time is increased.

Motivated by reduced-order models of large wind-turbine blades under steady conditions, in which gravity and aerodynamics provide cyclic stiffening and direct loading, we have looked at the linear Mathieu equation with combined parametric and direct excitations of the same frequency. The forced linear Mathieu equation has homogeneous and particular solutions relative to consideration of the direct excitation term. The homogeneous solution has all of the characteristics shown in classical studies of the Mathieu equation. Thus, our attention was on the forced responses, which we studied for hard cyclic excitation using a second-order multiple-scales analysis.

Following an earlier study which described the superharmonic resonance of order two using a first-order perturbation analysis, the second-order perturbation analysis was applied to see if there was any refinement in the analytical quality. The response amplitudes as predicted by the first- and second-order analyses were nearly the same. However, the second-order analysis predicted the peak to be slightly shifted from the superharmonic frequency, while first-order analysis predicted the peak to be at the superharmonic frequency.

In earlier studies, numerical simulations revealed a superharmonic resonance of order three which was not uncovered by a first-order multiple-scales analysis. Here, the second-order multiple-scales analysis with hard excitation revealed the existence of a stable superharmonic resonance of order three. This resonance peak is of order $ε$ in magnitude, which is small compared to the superharmonic resonance of order two. The possible existence of such a resonance may be important to wind turbines, as they are designed to operate at frequencies well below the first natural frequency, and a small resonance can induce a response that is significantly larger than that of the non-resonant oscillator with direct excitation only. Furthermore, a recursive examination indicated that an nth-order MMS analysis will reveal resonance conditions of order n + 1, although the amplitudes of such superharmonics are expected to diminish with n.

The analysis showed that cyclic external excitation does not incur a steady-state subharmonic resonance, but consistently described the subharmonic instability of the Mathieu equation. Numerical simulations show how this might be observed in a frequency sweep.

Together with the previous work on primary resonance [28], the second-order perturbation analysis unfolds superharmonic resonances of orders three and two, primary resonance amplification, and subharmonic instability. Numerical simulations suggest that superharmonic resonances of higher orders can occur. The detection of such resonances with higher orders of perturbation expansions has been mapped out, but the full analyses would be needed to describe them.

This material is based on work supported by National Science Foundation (CBET-0933292, CMMI-1335177). Any opinions, findings, and conclusions or recommendations expressed are those of the authors and do not necessarily reflect the views of the NSF.

There are no conflicts of interest.

The authors attest that all data for this study are included in the paper.