## Abstract

The present study deals with the response of a forced Mathieu equation with damping, with weak harmonic direct excitation at the same frequency as the parametric excitation. A second-order perturbation analysis using the method of multiple scales unfolds parametric amplification at primary resonance. The parametric effect on the primary resonance behavior occurs with a slow time scale of second-order, although the effect on the steady-state response is of order 1. As the parametric excitation level increases, the response at primary resonance stretches before becoming unbounded and unstable. Analytical expressions for predicting the response amplitudes are presented and compared with numerical results for a specific set of system parameters. Dependence of the amplification behavior, and indeed possible deamplification, on parameters is examined. The effect of parametric excitation on the response phase behavior is also presented.

## 1 Introduction

Our studies on the forced Mathieu equation were originally motivated by our interest in the dynamics of large wind-turbine blades. Horizontal-axis wind-turbine blades are seen as slender beams rotating with the hub about a nearly horizontal axis [1]. As the blade rotates through the gravitational field, it cycles from phases of tension (downward orientation), which stiffens the blade, through phases of compression (upright orientation), which softens the blade, together resulting in parametric stiffness. The blade cycles through phases of in-plane transverse gravitational loading against the leading edge and then the trailing edge, contributing to external cyclic forcing for in-plane deflections. Also, because of wind shear, the wind velocity is stronger at higher altitudes, such that the lift loading on out-of-plane deflections is also cyclic. When linearized, the single-mode model with in-plane or out-of-plane deflection for a rotor at steady speed $Ω$ reduces to a forced Mathieu equation,
$q¨+2ϵμq˙+(ω2+ϵγcosΩt)q=F(t)$
(1)
where q represents a modal displacement, ε is a small parameter, εμ/ω is the damping ratio, ω is the undamped modal frequency, γ is the scaled strength of variation of the effective stiffness due to gravity, and F(t) is the cyclic aerodynamic and direct gravitational (in-plane case) loading on the blade. The cyclic stiffening, εγ, may be significant in the wind-turbine application as the wind turbines become extremely large [24]. Other examples of modeling wind turbines with parametric excitations include [57]. Forced Mathieu systems also show up in other applications [811].

Since Eq. (1) is linear, the forced system is nonhomogeneous and its solution can be expressed as a sum of homogeneous and particular solutions, given as q = qh + qp, where qh satisfies the damped Mathieu equation, $q¨h+2ϵμq˙h+(ω2+ϵγcosΩt)qh=0$, and qp satisfies Eq. (1). There is a large body of work aimed at understanding the unforced Mathieu equation, and hence qh. It 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, even if the damping and time-varying stiffness are positive. These stability wedges can be studied by applying the Floquet theory [12] with harmonic balance solutions of periodic transitional motions [13] and also by a higher-order perturbation expansion [14]. For lightly damped cases, these regions of instability are based near frequency ratios $Ω/ω=2/N$, N being an integer. The dependence of the instability regions on the strength of damping has been further explored in detail in Ref. [15]. Free responses, whether stable or unstable, have also been analyzed based on the Floquet theory [16].

By linearity, if F(t) in Eq. (1) is scaled by a factor α, then the forced or particular solution qp is also scaled by α. (However, scaling the level γ of parametric excitation does not simply scale the response.) The solution to the forced Mathieu equation can be approximated by perturbation techniques to reveal resonance conditions as well as instabilities. Our previous work [17] dealt with the superharmonic and subharmonic resonances for the forced Mathieu equation with cubic nonlinearity. A first-order perturbation analysis using the method of multiple scales revealed the existence of superharmonics at one-half and one-third the natural frequency and subharmonics at twice and thrice the natural frequency. The first-order analysis showed that the superharmonic at order 2 persists for the linear system. Numerical simulations of the linear system showed a noticeable stretch in the primary resonance (Fig. 1), suggesting a parametric amplification, as well as a superharmonic of order 3, which were not uncovered by the first-order perturbation multiple-scales analysis [17]. The aim of this study is to apply second-order multiple scales to unfold the parametric amplification at primary resonance.

Fig. 1
Fig. 1
Close modal

Resonance conditions in parametrically excited systems with forcing and also nonlinearity have been examined [18], and there have been other studies on dynamics, stability, and bifurcations of nonlinear generalizations of forced Mathieu-type equations [1924].

Although our interest began with wind turbines, such structures are unlikely to operate near primary resonance, unless there is an uncontrolled runaway incident. However, in this work, our interest is in the forced Mathieu equation as a fundamental dynamical system. An example area of study has been subharmonic parametric amplification (with primary direct excitation), which has undergone numerous analytical and experimental studies [810,25,26] in microelectromechanical systems (MEMS) and nanoelectromechanical systems with applications [2730], including with arrays of resonators [3133]. These examples involve a 2:1 ratio of parametric to direct excitation frequencies, and as such, the subharmonic parametric resonances are exploited. In contrast, Eq. (1) has a 1:1 ratio between the frequencies of the parametric and direct excitation terms. Conceivably, other configurations of MEMS resonators might involve 1:1 excitations and be conducive to primary parametric amplification.

We aim to apply higher-order perturbation techniques to unfold the existence and features of primary parametric amplification in the forced Mathieu equation. Higher-order perturbation techniques have been used extensively across various sciences to get more accurate expressions of the solutions to problems, to derive an understanding of the underlying physics [34,35]. Dynamical systems have been described by higher-order perturbation expansions [14,36,37] as well as other scaling techniques [38].

Hence, in the current work, we analyze Eq. (1) under hard forcing with two orders of multiple-scale expansion as formulated by Nayfeh [36] to capture the primary resonance and its features. The primary parametric amplification was introduced in Refs. [3,39] and is further developed here to describe the amplification and possible suppression, phase effects, and the effects of parameters.

## 2 Weak Excitation—Primary Resonance Amplification

In this section, we perform second-order analysis to uncover resonance conditions under soft excitation (O(ε) forcing). Our primary interest here is in describing primary resonance amplification observed in simulations. However, we will take note of other resonances that may arise. We focus our attention on the original equation (1) with weak excitation, restated as follows:
$q¨+2ϵμq˙+(ω2+ϵγcosΩt)q=ϵFsin(Ωt+θ)$
(2)
We follow a second-order multiple-scales expansion analysis. To this end, we incorporate three time scales, (T0, T1, T2) and allow for a dominant solution q0 and small variations q1 and q2, such that
$q=q0(T0,T1,T2)+ϵq1(T0,T1,T2)+ϵ2q2(T0,T1,T2)+⋯$
(3)
where Ti = εit. Then, d/dt = D0 + εD1 + ε2D2, where Di = ∂/∂Ti. We substitute these into Eq. (2) and then simplify and balance the coefficients of ε0, ε1, and ε2 as follows:
$O(ϵ0):D02q0+ω2q0=0O(ϵ1):D02q1+ω2q1=−2μD0q0−2D0D1q0−γq0cosΩT0+Fsin(ΩT0+θ)O(ϵ2):D02q2+ω2q2=−2D0D1q1−(D12+2D0D2)q0−2μ(D0q1+D1q0)−γq1cosΩT0$
(4)
The solution of the O(ε0) equation is expressed as follows:
$q0=A(T1,T2)eiωT0+c.c.$
(5)
where $A=12aeiβ$ and c.c. means complex conjugate of the preceding terms.
The differential equation for q1 is then rewritten by substituting the solution for q0 into the O(ε) equation as follows:
$D02q1+ω2q1=−2iωD1AeiωT0−2μiωAeiωT0−γ(A2ei(ω+Ω)T0+A¯2ei(Ω−ω)T0)−iF2ei(ΩT0+θ)+c.c.$
(6)
where $A¯$ is the complex conjugate of A.

We have three cases for Eq. (6), which can contribute secular terms:

1. No specific relationship between $Ω$ and ω

2. $Ω≈2ω$

3. $Ω≈ω$

The treatment of these cases will be laid out in the following sections. In each case, the second-order expansions will give us the slow time scale variations of the amplitude A of our fast time scale solution q0. A second-order expansion will involve solvability conditions obtained by eliminating secular terms at both O(ε) and O(ε2).

### 2.1 Case 1: No Resonance at O(ε).

When there is no specific relationship between the forcing frequency $Ω$ and the natural frequency ω, we remove secular terms by letting −2iωD1A − 2μiωA = 0 and solve for the particular solution using the remaining terms in differential equation (6), noting that we treat A as a constant with respect to the independent variable T0. As such,
$q1=KAei(Ω+ω)T0+LA¯ei(Ω−ω)T0+NeiΩT0+c.c.$
(7)
where
$K=γ2(Ω2+2ωΩ),L=γ2(Ω2−2ωΩ),andN=iFeiθ2(Ω2−ω2)$
(8)
We substitute the solutions for q0 and q1 into Eq. (4) at O(ε2). From the resulting q2 differential equation, we seek to identify terms that can contribute to a resonance condition. Although there were no resonances specified at O(ε), new resonances can be identified at O(ε2). When there is not a specified relationship between $Ω$ and ω at either O(ε) or O(ε2), i.e., no resonances are specified, secular terms are removed and the resulting solvability conditions from both O(ε) and O(ε2) are expressed as follows:
$O(ϵ):−2iωD1A−2μiωA=0O(ϵ2):−D12A−2D2Aiω−2μD1A−γ2(KA+LA¯)=0$
(9)
Following Refs. [13,14], these equations can be analyzed for their fixed points and stability.
The removal of the secular terms in Eq. (4) at O(ε2) leaves a differential equation for q2 as follows:
$D02q2+ω2q2=−2(i(Ω+ω)KD1Aei(Ω+ω)T0+i(Ω−ω)LD1A¯ei(Ω−ω)T0)−(D12AeiωT0+2iωD2AeiωT0)−2μ(i(Ω+ω)KAei(Ω+ω)T0+i(Ω−ω)LA¯ei(Ω−ω)T0+iΩNeiΩT0+D1AeiωT0)−γ2(KAei(2Ω+ω)T0+KAeiωT0+LA¯ei(2Ω−ω)T0+LA¯e−iωT0+Nei2ΩT0+N)+c.c.$
(10)
for the nonresonant case, where K, L, and N are defined in Eq. (8).

But there are also specific combinations of $Ω$ and ω appearing at the second-order that could lead to resonance conditions in the system, and thereby add terms to the solvability condition at O(ε2). These can be identified by examining the remaining terms of the O(ε2) equation given by Eq. (10). We arrive at the following candidates:

• $Ω≈ω/2$

• $Ω≈2ω$

• $Ω≈ω$

The first case was treated with a first-order multiple-scales analysis in Ref. [17], in which hard excitation-induced secular terms at the O(ε) expansion. Similarly, in this case, the excitation, this time at O(ε), induced secular terms one order down at O(ε2). Hence, similar results are expected, and so this case is not emphasized in this article. The latter two cases actually produce secular terms at O(ε) and are therefore cannot be treated within case 1, since it is a nonresonant case at O(ε).

### 2.2 Case 2: Subharmonic Resonance of Order 2.

This case was examined in the previous first-order multiple-scales analysis [17], but we continue here to verify consistency in the second-order analysis. When $Ω≈2ω$, i.e., $Ω=2ω+ϵσ$, the solvability condition for Eq. (6) is expressed as follows:
$O(ϵ):−2iωD1A−2μiωA−γA¯2eiσT1=0$
(11)
The particular solution of Eq. (6) for q1 then becomes
$q1=KAei(Ω+ω)T0+NeiΩT0+c.c.$
(12)
Applying to Eq. (4) at O(ε2), and removing secular terms, leaves us with a solvability condition as
$O(ϵ2):−D12A−2D2Aiω−2μD1A−γ2KA=0$
(13)
We recombine Eqs. (11) and (13) by letting dA/dt = εD1A + ε2D2A [36], and apply $Ω≈2ω$ to K, to get an ordinary differential equation:
$−2iωdAdt+ϵ(−2μiωA−γA¯2eiσϵt)+ϵ2(μ2A+γσA¯4ωeiσϵt−γ2A32ω2)=0$
(14)

It is noticed that F does not appear in the solvability conditions and thus does not associate with a forced subharmonic response. There is a small linear forced-response term in Eq. (12), representing the forced component of the response. The solvability conditions thus provide the dynamic characteristics of the homogeneous solution and its dependence on the parametric excitation.

To examine this, we follow the steps to be detailed in the next section by letting A = (X + i Y)eiσεt in Eq. (14) to arrive at $X˙$ and $Y˙$ equations, which have the matrix form of $x˙=Ax$. The expressions for the trace (T) and determinant (D) of the resulting matrix A provide information about stability. The trace and determinant are given by
$T=−4ϵωμD=(2ϵωμ)2+(2ϵσω+ϵγ2+ϵ2μ2−ϵ2σγ4ω−3ϵ2γ232ω2)(2ϵσω−ϵγ2+ϵ2μ2+ϵ2σγ4ω−3ϵ2γ232ω2)=(2ϵωμ)2+(2ϵσω+ϵ2μ2−3ϵ2γ232ω2)2−(ϵγ2−ϵ2σγ4ω)2$
(15)

The expression for D is quadratic in γ2, which when equated to zero gives the stability boundaries of the order 2 subharmonic under hard excitation. These agree with analyses in the existing literature [13,14], and we do not present further details of this analysis. A first-order analysis, in fact, is enough to approximate the stability boundaries (from Eq. (11)), and the second-order analysis refines the expressions of the stability boundaries.

### 2.3 Case 3: Primary Resonance.

This case is the focus of this article. When $Ω≈ω$, i.e., letting $Ω=ω+ϵσ$ in Eq. (6), the solvability condition is expressed as follows:
$O(ϵ):−2iωD1A−2μiωA−iFeiθ2eiσT1=0$
(16)
The secular terms at this order do not involve parametric excitation, γ, and therefore do not uncover the parametric effect on primary resonance. The particular solution for q1 becomes
$q1=KAei(Ω+ω)T0+LA¯ei(Ω−ω)T0+c.c.$
(17)
where K and L are defined in Eq. (8).

We substitute the solution for q0 and q1 into the Eq. (4) at O(ε2). After expanding the right-hand side of the resulting q2 differential equation, we seek to identify secular terms at this resonance condition.

For the case of primary resonance, we let $Ω=ω+ϵσ$ in the last of Eqs. (4) and eliminate secular terms to get
$O(ϵ2):−D12A−2D2Aiω−2μD1A−γ2((K+L)A+LA¯e2iσT1)=0$
(18)
These second-order secular terms do indeed involve the parametric excitation γ.
Rewriting the O(ε) equation in Eq. (16), we have
$D1A=−μA−Feiθ4ωeiσT1$
We compute $D12A$ from the aforementioned expression as follows:
$D12A=μ2A+μFeiθ4ωeiσT1−iσFeiθ4ωeiσT1$
We substitute these into the O(ε2) expression in Eq. (18) to write
$−2D2Aiω−μ2A−μFeiθ4ωeiσT1+iσFeiθ4ωeiσT1−2μ(−μA−Feiθ4ωeiσT1)−γ2((K+L)A+LA¯e2iσT1)=0$
(19)
Applying reconstitution as in Refs. [14,36] by letting dA/dt = εD1A + ε2D2A, and substituting $Ω≈ω$ into K and L, the O(ε) solvability condition from Eq. (18) and the O(ε2) solvability condition from Eq. (19) are seen to result from a multiple-scales expansion of
$−2iωdAdt+ϵ(−2μiωA−iFeiθ2eiσϵt)+ϵ2(μ2A+γ2A¯4ω2e2iσϵt+γ2A6ω2+iσFeiθ4ωeiσϵt+μFeiθ4ωeiσϵt)=0$
(20)

The second-order analysis indicates that the parametric excitation affects the resonant amplitude at a very slow time scale, i.e., in the ε2 terms in the reconstituted differential equation for A.

## 3 Behavior of Primary Resonance: Parametric Amplification

In this section, we study the amplitude and the phase of responses at primary resonance. Transforming Eq. (20) to polar coordinates via $A=12aeiβ$ and letting β = σT1ϕ leads to slow-flow equations in a and ϕ. The steady-state equations include sin (ϕ), cos (ϕ) and sin (2ϕ), andcos (2ϕ) terms, which make them more difficult to solve than other familiar examples.

Instead, following Ref. [14], we let A = (X + i Y)eiσεt, with real X and Y, in Eq. (20), separate imaginary and real parts, and cancel the common exponential term to obtain
$2ωdXdt+2ϵμωX−(2ϵωσ+ϵ2μ2−ϵ2γ212ω2)Y=ϵα$
(21)
$2ωdYdt+(2ϵωσ+ϵ2μ2+5ϵ2γ212ω2)X+2ϵμωY=ϵδ$
(22)
where εα = −(εF/2)cosθ + (ε2/4ω)cosθ + (ε2/4ω)sinθ and εδ = −(εF/2)sinθ + (ε2/4ω)sinθ − (ε2/4ω)cosθ.
At steady-state $X˙=Y˙=0$. Hence, from Eqs. (22) and (21), we get
$2ϵμωX+−(2ϵωσ+ϵ2μ2−ϵ2γ212ω2)Y=ϵα(2ϵωσ+ϵ2μ2+5ϵ2γ212ω2)X+2ϵμωY=ϵδ$
(23)

### 3.1 Parametric Amplification.

We solve (23) for X and Y and then write the square of the total amplitude r2 = X2 + Y2 as follows:
$r2=[2ωμα+δ(2ωσ+ϵμ2−ϵγ212ω2)]2[(2ωμ)2+(2ωσ+ϵμ2−ϵγ212ω2)(2ωσ+ϵμ2+5ϵγ212ω2)]2+[2ωμδ−α(2ωσ+ϵμ2+5ϵγ212ω2)]2[(2ωμ)2+(2ωσ+ϵμ2−ϵγ212ω2)(2ωσ+ϵμ2+5ϵγ212ω2)]2$
(24)

We have noted that γ affects the response amplitude at the second-order slow time scale. However, the steady-state response amplitude, via the expression for r2, is at order O(ε0).

Equation (24) provides information related to the amplitude of the primary resonance as a function of parameters, including the detuning parameter σ and the phase angle θ between the parametric and direct excitation terms. Since α and δ depend on θ, the response amplitude is also dependent on θ.

Note that we had A = (X + iY)eiεσt = Beiεσt and $Ω=ω+ϵσ$. Letting $B=12ae−iϕ$, the leading-order solution (5) takes the form
$q0=BeiωT0+c.c.=BeiΩT0+c.c.=acos(ΩT0−ϕ)$
(25)
Thus, the response is at the excitation frequency, and ϕ is the phase lag relative to the parametric excitation. Also, since |X + i Y| = r, the response amplitude is a = 2r.
Looking, for example, at the special case of θ = 0, we have
$r02=[(2ωσ+ϵμ2−ϵγ212ω2)(−ϵFμ4ω)−(F2−ϵFσ4ω)(2ωμ)]2[(2ωμ)2+(2ωσ+ϵμ2−ϵγ212ω2)(2ωσ+ϵμ2+5ϵγ212ω2)]2+[(2ωμ)(−ϵFμ4ω)+(F2−ϵFσ4ω)(2ωσ+ϵμ2+5ϵγ212ω2)]2[(2ωμ)2+(2ωσ+ϵμ2−ϵγ212ω2)(2ωσ+ϵμ2+5ϵγ212ω2)]2$
(26)
Similarly, the case of θ = π/2 yields
$rπ/22=[2ωμ(ϵFμ4ω)+(−F2+ϵFσ4ω)(2ωσ+ϵμ2−ϵγ212ω2)]2[(2ωμ)2+(2ωσ+ϵμ2−ϵγ212ω2)(2ωσ+ϵμ2+5ϵγ212ω2)]2+[2ωμ(−F2+ϵFσ4ω)−(ϵFμ4ω)(2ωσ+ϵμ2+5ϵγ212ω2)]2[(2ωμ)2+(2ωσ+ϵμ2−ϵγ212ω2)(2ωσ+ϵμ2+5ϵγ212ω2)]2$
(27)
In the limit as ε → 0, σ → 0 (and using a = 2r), Eqs. (24) and (26) are consistent with the common forced harmonic oscillator.
Equation (24) (or special cases of Eq. (26) or (27)) provides a predicted response amplitude as a function of system parameters for the 1:1 primary resonance of the linear forced Mathieu equation with damping. The response phase can also be obtained from X and Y. The steady-state solution (25) can be recast in the form of
$q0=asin(ΩT0+θ−ψ)$
(28)
where the phase lag ψ, relative to the direct excitation, can be determined from X and Y. To gain insight into the behavior of the response amplitude and phase, let us look at some numerical plots and simulations.

Figure 2(a) plots the amplitude a = 2r, where r comes from Eq. (26) for θ = 0 (meaning the direct and parametric excitations are out of phase by π/2, since they are present as sine and cosine in Eq. (2)), for nondimensional values of parametric excitation level γ ranging from γ = 0 to γ = 4. The parameter values are ε = 0.1, ω = 1, μ = 0.25, and F = 5, and also nondimensional. The plot shows a clear resonant peak. This peak occurs at a value of σp (equivalently $Ωp=ω+ϵσp$), which is very close to σ = 0 when γ is small, as γ increases, the peak moves to lower values of $Ω$. When γ < 1.5, the resonance curve is indistinguishable from the curve of the regular forced linear oscillator. When γ = 1.5, the curves have a slight perceptible separation, which grows quadratically as γ increases further, indicating a parametric amplification of primary resonance. The case of γ = 4 is cropped for better visualization of the other curves, but remains bounded for all $Ω$. For slightly larger values of γ, the primary resonance becomes unbounded for some values of $Ω$. We will further discuss this shortly. For the cases of γ = 3 and γ = 3.5, simulation values are plotted. These were obtained by numerically solving the differential equation (using matlab ode45) and examining the fast Fourier transform (FFT) of 25 cycles of the response at the steady state. The sampling interval was 1/100 of the excitation period, such that leakage issues did not distort the FFT. The FFT indicated a strong peak at the primary resonance excitation frequency (close to the natural frequency) and a minor peak at its first harmonic. The amplitude of the primary term only is plotted in Fig. 2 and subsequent plots. The next harmonic (at twice the natural frequency) can be seen in the q1 expression in Eq. (17). For example, when ε = 0.1, γ = 3, and $Ω=1$, this term should be about 1/20 of the primary resonance amplitude; this is very consistent with the detected first and second harmonic amplitudes of 13.97 and 0.6988. The simulations shown have excellent agreement. Not shown are simulations samples for γ = 4. At this value of γ, some error can be detected. For example, when $Ω=ω$, there is a 4.6% error between the simulated response amplitude and that predicted by the second-order perturbation expansion.

Fig. 2
Fig. 2
Close modal

Figure 2(b) shows the associated phase lag ψ behind the direct excitation for the same parameter values. The plot magnifies the phase transition region near the frequency ratio of $Ω/ω=1$. For γ = 1, the phase lag overlaps with the nonparametric case (γ = 0). Increasing γ pushes the phase transition to a slightly lower frequency ratio, along with the trend of the response amplitude peak. Otherwise, the forced and parametrically excited system exhibits a phase lag, which approaches 0 for very low frequencies, approximately π/2 radians near resonance, and approaches π radians for very large frequencies.

Figure 3(a) plots the amplitude a = 2r, from Eq. (26) for θ = π/4 for nondimensional values of parametric excitation level γ ranging from γ = 0 to γ = 4, and for the same values of parameters as the previous figure. In this case, the plotted values of γ of 2 or less do not show a deviation from the regular forced linear oscillator, which is represented by the solid curve with a single local maximum near $Ω=1$. As γ increases beyond a value of about 2, the resonance peak first decreases. Thus, this value of θ shows parametric deamplification. At the value of γ = 3.5, the resonance peak begins to split. The two peaks that form increase as γ continues to increase. Further slight increases in γ are accompanied by large growth in the amplitude before becoming unbounded. Simulations for the cases of γ = 3.5 and γ = 4 show strong agreement with the perturbation results. Figure 3(b) shows the associated phase lag ψ for θ = π/4 and the same parameter values. For γ = 1, the phase lag overlaps with the nonparametric case (γ = 0). Increasing γ causes the phase transition to become steeper, i.e., more abrupt with frequency changes near resonance. This trend, with increasing γ, is reminiscent of the trend seen in a regular (nonparametric) forced oscillator as damping decreases. Curiously, at this value of θ, the effect on the response amplitude (Fig. 3(a)) when increasing γ has a trend opposite that caused by decreasing damping. The phase steepening occurs with a crossover frequency, below which the phase is decreased and above which the phase is increased, with increasing γ. At the crossover point, the phase is independent of γ.

Fig. 3
Fig. 3
Close modal

Insight can also be obtained from the plot of the complex amplitude in the complex plane. Figure 4 shows plots of X+iY, representing half complex amplitudes, phased relative to the direct excitation, for the cases of (a) θ = 0 and (b) θ = π/4, with various parametric levels γ. (In this case, the various levels of γ are plotted with solid curves to maintain a cleaner image.) The curves are parameterized by −1 < σ < 1 (equivalently $Ω$). If the range σ were increased, the plots would approach closed curves. Figure 4(a) shows that the effect of γ is to extend the amplitudes from the γ = 0 reference, which is the smaller circle, to the oval curves that increase as γ increases. The plot shows that the maximum amplitude occurs with a phase lag ψ > π/2 radians (π/2 being the phase at resonance of the reference γ = 0 system). Figure 4(b) shows that as γ increases, the reference circle distorts into the oval curves that extend to both sides. From this, it can be seen that if γ is large enough, there can be two amplitude peaks associated with phase lags that are closer to 0 and π radians than to the reference phase of π/2. This feature can be cross checked in Fig. 3.

Fig. 4
Fig. 4
Close modal

In both examples of Fig. 4, there are two crossover points, which indicate, for any γ, there are two frequencies that produce a common amplitude and phase. We have observed that these events occur, however, at different frequencies $Ω$ for different values of γ (not shown here).

Figure 5 delves more into the effect of θ on the resonance behavior, for a fixed value of γ = 4. The top left shows how the half complex amplitude (X+iY) varies with 0 < θ < π at increments of π/8 radians. θ = 0 is represented by the dashed curve and also shows up in Fig. 4(a). The oval rotates clockwise with increase in θ. At θ = π/4, the oval is small and horizontal, represented by the tightly hatched curve. Continued increase in θ rotates the oval and increases it size. The dotted curve represents θ = 3π/4 and leads to the maximum resonance amplification of the cases depicted. Continued increase to θ = π brings the ovals back to the dashed case, and the process repeats. The solid black circular curve shows the case of γ = 0 for reference.

Fig. 5
Fig. 5
Close modal

The right column of Fig. 5 separates the amplitude and the phase of the response. The dotted, hatched, and dashed curves are consistent with those in the top left image, and the other curves show the transitions with varying θ. The solid black curve (lowest single maximum) shows γ = 0 for reference. The top right curves show that the resonance occurs at a slightly lower frequency than the regular (γ = 0) case. θ has a strong effect on the character of the peak. The case of θ = π/4 is the hatched curve with twin peaks. It can be cross checked in Fig. 3(a). The lower right image shows the associated phase lags. The case of θ = π/4 is the hatched (looks dotted in this image, however) curve that is steepest. It can be cross checked in Fig. 3(b). The γ = 0 reference is the solid black curve that transits π/2 radians at $Ω=1$. The general, but not strict trend, is that for most θ, the phase transition occurs at slightly lower frequencies than $Ω=1$, consistent with the amplitude graph. Also, θ can be associated with steeper transitions (e.g., θ = π/4) or softer transitions (e.g., the dotted curve of θ = 3π/4).

Steepening the phase curve, and amplifying the amplitude peak, might be suspected to be “equivalent damping reduction” effects. When θ = π/4, the amplitude peak is suppressed (like increased damping), while the phase curve is more steep (like decreased damping). When θ = 3π/4, both effects are reversed. Hence, if one tries to look for “equivalent damping effects,” the parametric resonance has conflicting patterns.

Thus, θ can have a significant effect on the parametric amplification [8,10,26]. In this case, depending on the value of θ, parametric-plus-direct excitation can lead to parametric amplification or suppression. We can define a gain as the ratio of the response amplitude a = 2r to that of the system without parametric excitation (γ = 0) [8]. For the latter, we use the resonant ($Ω=ω$) value for the simple harmonic oscillator, namely, $ah=F/2μω2$. Figure 6 shows the dependence of the maximum gain on the phase θ, for a family of values of γ ranging from γ = 0 at the horizontal line (unit gain) to γ = 4 at the highest curve, with the same parameters as mentioned earlier. This figure was produced by choosing a value of θ and then sampling Eq. (24) to determine the maximum amplitude. The maximum amplitude is then normalized by ah to produce the gain and plotted versus θ. The figure shows that the maximum amplification occurs near θ ≈ 3π/4, and the minimum amplification, indeed suppression, occurs near θπ/4. However, for most values of θ, the parametric effect amplifies the resonance. The curve is periodic in θ, with one period depicted. As γ increases, the amplification becomes larger (and in cases of deamplification, the suppression becomes less pronounced until it eventually shows slight amplification). At some point, γ becomes large enough that the maximum response becomes uniformly unbounded. (For these parameters, this occurs around γ =4.5.)

Fig. 6
Fig. 6
Close modal

Figure 6(b) shows that as damping decreases, the maximum system gain uniformly increases. Much like when γ achieved a high enough value, it turns out that when μ reaches a small enough value, the system becomes unbounded.

Figure 7 shows the amplification as a function of parametric excitation level γ, for various values of θ. The graph indicates that for γ near zero, the gain has a zero slope, and toward a value of γ slightly larger than 4, the gain becomes unbounded (as it enters a parametric instability). Similar γ dependence has been observed for 2:1 parametric pumps [8,26,28]. For values of θ for which the system is parametrically de-amplified, with large enough γ, there is a parametric amplification up to a point at which the system becomes unbounded. In the figure, the lowest curve corresponds to θ = π/4, and the highest curve corresponds to θ = 3π/4, with equal increments between. Similar phase dependence has been observed for 2:1 parametric pumps [8,26,28]. For small γ, there is little amplification effect. As γ increases, the amplification increases sharply (or suppression develops gently) before becoming unbounded. There may be an opportunity to introduce nonlinearity to broaden the range of γ for significant amplification and provide robustness in applications [10].

Fig. 7
Fig. 7
Close modal

The features of Figs. 2 and 3 show that the resonance may have a single peak or two peaks or may be unbounded. With damping present, we might in principle seek the presence of one or more peaks by finding zeros of d(r2)/. However, this is a complicated calculation, and thus, we do not pursue it here. Furthermore, since the peaks can stray from $Ω=ω$ (σ = 0), it may not be useful to approximate the peak value with an $Ω=ω$ (or σ = 0), as is equivalently done to estimate the resonance peak of a lightly damped regular linear forced vibration system. The upshot is that it is not easy to find a simple expression for approximating the peak response under parametric amplification.

### 3.2 Stability of Primary Resonance.

Equations (22) and (21) can be rewritten in the following form:
$2ωx˙=Ax+F$
(29)
where, for example, in the case of θ = 0,
$x=[XY],F=[−ϵF2+ϵ2Fσ4ω−ϵ2Fμ4ω],andA=[−2ϵμω2ϵωσ+ϵ2μ2−ϵ2γ212ω2−2ϵωσ−ϵ2μ2−5ϵ2γ212ω2−2ϵμω]$
Note that while F is written in the case of θ = 0, A is independent of θ.

The equilibrium x0 = A−1F was found previously, and its amplitude r was studied. Since Eq. (29) is linear, dynamics of X and Y, and hence the leading-order slowly varying response amplitude $A=12a(T1,T2)eiβ(T1,T2)$, about the equilibrium, and hence the stability, are governed by matrix A. In addition, the homogeneous solution has its own equilibrium at the origin. The homogeneous solution is that for which F = 0 in Eq. (2) and consequently F = 0 in Eq. (29). Thus, matrix A simultaneously governs the stability of the homogeneous solution and the characteristics of the particular solutions of Eq. (2).

The stability is determined by the trace (T) and determinant (D) of A, which are as follows:
$T=−4ϵμωD=4ϵ2ω2(μ2+σ2)+4ϵ3μ2ωσ+2ϵ3γ2σ3ω+ϵ4μ4+ϵ4μ2γ23ω2−5ϵ4γ4144ω4$
(30)
Since ε > 0, μ > 0 and ω > 0, the trace T < 0. Therefore, the stability is governed by the determinant D. The instability transition is given by parameters for which D = 0. Close examination of Eq. (26) shows that the denominator of the r2 expression is in fact D2, consistent with the denominator squared of A−1. Thus, the boundedness of the forced response coincides with the approach of r2 to infinity and with the instability of xh of the unforced Mathieu equation. Solving D = 0 for γ yields the critical parametric excitation as follows:
$γcr=24ω25ϵ(2ωσ+ϵμ2)±12ω25ϵ9(2ωσ+ϵμ2)2+5(2μω)2$
(31)
The second term is larger than the first term, and hence, the positive inner root leads to real values for γcr.

The calculated value of γcr is the value of γ in Fig. 7 for which the maximum amplitude becomes unbounded for specific values of detuning parameter. The increase in amplitude until that point is a result of parametric amplification, analytically approximated by Eq. (26). The analytically predicted values of γcr are plotted versus σ in Fig. 8. Numerical simulations were also checked for stability in this parameter space, and stable (dots) and unstable (×’s) events are plotted along with the predicted stability boundary from Eq. (31) with good agreement. For the parameters of the plot, the minimum γ for which the system destabilizes is near γ = 4.5, which is consistent with the implied asymptote in Fig. 7.

Fig. 8
Fig. 8
Close modal

Looking at Eq. (26), for the lightly damped case with sufficiently large γ, there are two values of σ (equivalently $Ω$) for which the denominator is zero and the amplitude is unbounded. Thus, there are two frequencies, $Ω1$ and $Ω2$, which meet the stability-wedge boundary. For $Ω1<Ω<Ω2$, the response is unstable and unbounded. Fixing γ and increasing $Ω→Ω1$, the response amplitude plot approaches a vertical asymptote, while decreasing $Ω→Ω2$, the response amplitude plot approaches another vertical asymptote. The sides of the unbounded resonance peak surround an interval of values of $Ω$, in contrast to a single frequency of unbounded responses in the regular (γ = 0) oscillator.

### 3.3 Summary of Primary Parametric Resonance.

We have seen that the linear forced Mathieu equation has a homogeneous and particular solution. The homogeneous solution has the known behavior of the Mathieu equation. The particular solution is a forced response with the potential for a primary resonance, which was analyzed in this section by multiple scales. A second-order analysis shows that the primary resonance is stretched, or amplified (for most values of θ), by the parametric excitation. Given the parameter set, as the level of parametric excitation γ grows, the resonance amplitude is at first hardly affected by small γ, but then begins to grow sharply for larger γ, and approaches infinity at a critical γcr. The bounded resonance response is stable. As the response becomes infinite at γcr, the response destabilizes. At the same time, the homogeneous part of the response also becomes unstable as the parameters cross into the unstable part of the Arnold tongue.

### 3.4 Cautionary Note on Second-Order Expansions in ε.

We have done an expansion based on powers of ε and derived a predictive amplitude equation (24), which is a function of parameters. Equation (24) has a denominator which is a function of ε. There may be a question as to whether this expression for the amplitude should be expanded into powers of ε, retaining up to ε2, for consistency with the perturbation expansion. It turns out, in our example, that such a treatment reduces the accuracy of the amplitude expression.

To check this, we have expanded Eq. (24) and omitted terms of O(ε3) and higher (not shown). For example, for μ = 0.25, ε = 0.1, F = 5, γ = 3 and σ = 0, the amplitude obtained from Fig. 1 at primary resonance ($Ω≈ω$) is a = 13.95. The numerical simulation of Eq. (2), as described earlier, led to the resonant harmonic amplitude of as = 13.97. However, the quadratic expansion of Eq. (24) produced an estimation of ax = 13.3.

Thus, it seems that the reconstituted equation (20) has information of a second-order expansion, and the approximate steady-state solution (26) uses all of this information. A truncated expansion of (26) loses information and deviates from the true solution.

## 4 Conclusion

We have looked at the linear Mathieu equation with combined parametric and direct excitation at 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. In an earlier study, numerical simulations revealed resonance phenomena that were not uncovered by a first-order multiple-scales analysis, including a parametric amplification at primary resonance, noticeable to a trained eye.

In this study, we performed a second-order multiple-scales analysis with weak excitation to uncover the parametric amplification effects on primary resonance. The work complements parametric amplification studies exploiting subharmonic resonance. With a parametric excitation of $γcosΩt$ and a direct excitation of $Fsin(Ωt+θ)$, the parametric amplification shows a stretching of the primary resonance peak with the increasing strength γ of parametric excitation, until a critical value γcr at which the resonance becomes unbounded. Amplitudes of the numerical solutions were compared with the analytical expression to very high agreement. The response phase lag relative to the direct excitation varies with the excitation frequency from zero to π radians as with a standard linear oscillator. The parametric term tends to push the resonance peak and phase-lag transition to a frequency slightly less than the natural frequency.

We looked at the role of a phase θ in the direct excitation term. Depending on parameters, it may be possible to choose a phase such that the parametric term suppresses the resonance. Under most phases, the parametric term has an amplifying effect. Variation in the response phase lag with respect to the direct-parametric phase difference was also examined. It is possible to have an excitation phase such that there is a frequency for which the response phase lag is independent of the parametric excitation level γ.

The amplified primary resonance goes unstable coincidentally as it becomes unbounded. This instability coincides with that of the 1:1 wedge in the unforced Mathieu equation. Thus, at primary resonance, the forced response and the homogeneous solution destabilize together.

Parametric amplification could be applied by design to accentuate (or suppress) resonance quality. However, care should be taken regarding the possibility of instability associated with the slender instability wedge associated with the unforced Mathieu equation near the natural frequency. The next study could address effects of stiffness nonlinearity on bounding the unstable behavior or bending the resonance, as done for the case when the parametric excitation has twice the frequency of the direct excitation [10].

The problem studied here shows an example where the second-order perturbation analysis uncovers behavior of nominally O(ε0) (or greater based on numerical results), but is not uncovered in the first-order perturbation analysis. The parametric effect at primary resonance occurs with a second-order slow time scale.

## Acknowledgment

This material is relevant to projects supported by the National Science Foundation (CBET-0933292 and 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.

## Conflicts of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The data and information that support the findings of this article are freely available. The authors attest that all data for this study are included in the paper.

## References

1.
Ramakrishnan
,
V.
, and
Feeny
,
B. F.
,
2011
, “
In-Plane Nonlinear Dynamics of Wind Turbine Blades
,”
ASME International Design Engineering Technical Conferences, 23th Biennial Conference on Vibration and Noise
,
Washington, DC
,
Aug. 28–31
, Paper No. DETC2011–48219.
2.
Allen
,
M. S.
,
Sracic
,
M. W.
,
Chauhan
,
S.
, and
Hansen
,
M. H.
,
2011
, “
Output-Only Modal Analysis of Linear Time-Periodic Systems With Application to Wind Turbine Simulation Data
,”
Mech. Syst. Signal Process.
,
25
(
4
), pp.
1174
1191
.
3.
Ramakrishnan
,
V.
,
2017
, “
,” Ph.D. thesis,
Michigan State University
,
East Lansing
.
4.
Acar
,
G.
, and
Feeny
,
B. F.
,
2018
, “
Bend-Bend-Twist Vibrations of a Wind Turbine Blade
,”
Wind Energy
,
21
(
1
), pp.
15
28
.
5.
Inoue
,
T.
,
Ishida
,
Y.
, and
Kiyohara
,
T.
,
2012
, “
Nonlinear Vibration Analysis of the Wind Turbine Blade (Occurrence of the Superharmonic Resonance in the Out-of-Plane Vibration of the Elastic Blade)
,”
ASME J. Vib. Acoust.
,
134
(
3
), p.
031009
.
6.
Ikeda
,
T.
,
Harata
,
Y.
, and
Ishida
,
Y.
,
2018
, “
Parametric Instability and Localization of Vibrations in Three-Blade Wind Turbines
,”
ASME J. Comput. Nonlinear Dyn.
,
13
(
7
), p.
071001
.
7.
Acar
,
G. D.
,
Acar
,
M. A.
, and
Feeny
,
B. F.
,
2020
, “
Parametric Resonances of a Three-Blade-Rotor System With Reference to Wind Turbines
,”
ASME J. Vib. Acoust.
,
142
(
2
), p.
021013
.
8.
Rugar
,
D.
, and
Grutter
,
P.
,
1991
, “
Mechanical Parametric Amplification and Thermomechanical Noise Squeezing
,”
Phys. Rev. Lett.
,
67
(
6
), pp.
699
702
.
9.
,
J. F.
,
Shaw
,
S. W.
,
Turner
,
K. L.
,
Moehlis
,
J.
,
DeMartini
,
B. E.
, and
Zhang
,
W.
,
2006
, “
Generalized Parametric Resonance in Electrostatically Actuated Microelectromechanical Oscillators
,”
J. Sound Vib.
,
296
(
4–5
), pp.
797
829
.
10.
,
J. F.
, and
Shaw
,
S. W.
,
2010
, “
The Impact of Nonlinearity on Degenerate Parametric Amplifiers
,”
Appl. Phys. Lett.
,
96
(
23
), p.
234101
.
11.
,
M. A.
, and
Sapsis
,
T.
,
2016
, “
Probabilistic Response and Rare Events in Mathieu’s Equation Under Correlated Parametric Excitation
,”
Ocean Eng.
,
120
, pp.
289
297
.
12.
McLachlan
,
N.
,
1964
,
Theory and Application of Mathieu Functions
,
Dover Publications
,
New York
.
13.
14.
Nayfeh
,
A. H.
, and
Mook
,
D. T.
,
1979
,
Nonlinear Oscillations
,
Wiley Interscience Publications. John Wiley and Sons
,
New York
.
15.
Susskind
,
H. J.
,
1991
, “
A Stability Analysis of the Mathieu Equation With Order-One Damping
,” Master’s thesis,
Cornell University
,
Ithaca, NY
.
16.
Acar
,
G.
, and
Feeny
,
B. F.
,
2016
, “
Floquet-Based Analysis of General Responses of the Mathieu Equation
,”
ASME J. Vib. Acoust.
,
138
(
4
), p.
041017
.
17.
Ramakrishnan
,
V.
, and
Feeny
,
B. F.
,
2012
, “
Resonances of the Forced Mathieu Equation With Reference to Wind Turbine Blades
,”
ASME J. Vib. Acoust.
,
134
(
6
), p.
064501
.
18.
Sharma
,
A.
,
2020
, “
A Re-Examination of Various Resonances in Parametrically Excited Systems
,”
ASME J. Vib. Acoust.
,
142
(
3
), p.
031010
.
19.
Belhaq
,
M.
, and
Houssni
,
M.
,
1999
, “
Quasi-Periodic Oscillations, Chaos and Suppression of Chaos in a Nonlinear Oscillator Driven by Parametric and External Excitations
,”
Nonlinear Dyn.
,
18
(
1
), pp.
1
24
.
20.
Pandey
,
M.
,
Rand
,
R.
, and
Zehnder
,
A. T.
,
2007
, “
Frequency Locking in a Forced Mathieu-van der Pol-Duffing System
,”
Nonlinear Dyn.
,
54
(
1–2
), pp.
3
12
.
21.
Newman
,
W. I.
,
Rand
,
R.
, and
Newman
,
A. L.
,
1999
, “
Dynamics of a Nonlinear Parametrically Excited Partial Differential Equation
,”
Chaos
,
9
(
1
), pp.
242
253
.
22.
Ng
,
L.
, and
Rand
,
R.
,
2002
, “
Bifurcations in a Mathieu Equation With Cubic Nonlinearities
,”
Chaos Solitions Fractals
,
14
(
2
), pp.
173
181
.
23.
Marghitu
,
D. B.
,
Sinha
,
S. C.
, and
Boghiu
,
D.
,
1998
, “
Stability and Control of a Parametrically Excited Rotating System. Part 1: Stability Analysis
,”
Dyn. Control
,
8
(
1
), pp.
7
20
.
24.
Tondl
,
A.
, and
Ecker
,
H.
,
2003
, “
On the Problem of Self-Excited Vibration Quenching by Means of Parametric Excitation
,”
Appl. Mech.
,
72
(
11–12
), pp.
923
932
.
25.
Carr
,
D. W.
,
Evoy
,
S.
,
Sekaric
,
L.
,
,
H. G.
, and
Parpia
,
J. M.
,
2000
, “
Parametric Amplification in a Torsional Microresonator
,”
Appl. Phys. Lett.
,
7
(
10
), pp.
1545
1547
.
26.
,
J.
,
Miller
,
N.
,
Shaw
,
S.
, and
Feeny
,
B.
,
2008
, “
Mechanical Domain Parametric Amplification
,”
ASME J. Vib. Acoust.
,
130
(
6
), p.
061006
.
27.
Dana
,
A.
,
Ho
,
F.
, and
Yamamoto
,
Y.
,
1998
, “
Mechanical Parametric Amplification in Piezoresistive Gallium Arsenide Microcantilevers
,”
Appl. Phys. Lett.
,
72
(
10
), pp.
1152
1154
.
28.
Zalalutdinov
,
M.
,
Olkhovets
,
A.
,
Zehnder
,
A.
,
Ilic
,
B.
,
Czaplewski
,
D.
,
,
H. G.
, and
Parpia
,
J. M.
,
2001
, “
Optically Pumped Parametric Amplification for Micromechanical Oscillators
,”
Appl. Phys. Lett.
,
78
(
20
), pp.
3142
3144
.
29.
Ouisse
,
T.
,
Stark
,
M.
,
Rodrigues-Martins
,
F.
,
Bercu
,
B.
,
Huant
,
S.
, and
Chevrier
,
J.
,
2005
, “
Theory of Electric Force Microscopy in the Parametric Amplification Regime
,”
Phys. Rev. B
,
71
(
20
), p.
205404
.
30.
Ono
,
T.
,
Wakamatsu
,
H.
, and
Esashi
,
M.
,
2005
, “
Parametrically Amplified Thermal Resonant Sensor With Pseudo-Cooling Effect
,”
J. Micromech. Microeng.
,
15
(
12
), pp.
2282
2288
.
31.
,
J.-P.
,
Brown
,
A. R.
,
Khuri-Yakub
,
B. T.
, and
Rebeiz
,
G. M.
,
2000
, “
A Novel Parametric-Effect MEMS Amplifier
,”
J. Microelectromech. Syst.
,
9
(
4
), pp.
528
537
.
32.
,
R.
, and
Turner
,
K. L.
,
2003
, “
Mechanical Domain Coupled Mode Parametric Resonance and Amplification in a Torsional Mode Micro Electro Mechanical Oscillator
,”
J. Micromech. Microeng.
,
13
(
5
), pp.
701
707
.
33.
Wallin
,
C. B.
,
De Alba
,
R.
,
Westy
,
D.
,
Holland
,
G.
,
Grutzik
,
S.
,
Rand
,
R. H.
,
Zehnder
,
A. T.
,
Aksyuk
,
V. A.
,
Krylov
,
S.
, and
Ilic
,
B. R.
,
2018
, “
Nondegenerate Parametric Resonance in Large Ensembles of Coupled Micromechanical Cantilevers With Varying Natural Frequencies
,”
Phys. Rev. Lett.
,
121
(
26
), p.
264301
.
34.
Holmes
,
M. H.
,
1995
,
Introduction to Perturbation Methods
,
Springer-Verlag
,
New York
.
35.
Murdock
,
J. A.
,
1991
,
Perturbations: Theory and Methods
,
A Wiley-Interscience Publication
,
New York
.
36.
Nayfeh
,
A. H.
,
1986
, “Perturbation Methods in Nonlinear Dynamics” (
Lecture Notes in Physics
), Vol.
247
, M. Jowett, M. Month, and S. Turner, eds.,
Springer-Verlag
,
Berlin
, pp.
238
314
.
37.
Sayed
,
M.
, and
Hamed
,
Y. S.
,
2011
, “
Stability and Response of a Nonlinear Coupled Pitch-Roll Shi Model Under Parametric and Harmonic Exitations
,”
Nonlinear Dyn.
,
64
(
3
), pp.
207
220
.
38.
Romero
,
L. A.
,
Torczynski
,
J. R.
, and
Kraynik
,
A. M.
,
2011
, “
A Scaling Law Near the Primary Resonance of the Quasiperiodic Mathieu Equation
,”
Nonlinear Dyn.
,
64
(
4
), pp.
395
408
.
39.
Ramakrishnan
,
V.
, and
Feeny
,
B. F.
,
2012
, “
Second-Order Multiple-Scales Analysis of the Nonlinear Forced Mathieu Equation
,”
ASME International Design Engineering Technical Conferences, 24th Conference on Vibration and Noise
, Chicago, Paper No. DETC2012-71532.