## Abstract

The Floquet theory has been classically used to study the stability characteristics of linear dynamic systems with periodic coefficients and is commonly applied to Mathieu’s equation, which has parametric stiffness. The focus of this article is to study the response characteristics of a linear oscillator for which the damping coefficient varies periodically in time. The Floquet theory is used to determine the effects of mean plus cyclic damping on the Floquet multipliers. An approximate Floquet solution, which includes an exponential part and a periodic part that is represented by a truncated Fourier series, is then applied to the oscillator. Based on the periodic part, the harmonic balance method is used to obtain the Fourier coefficients and Floquet exponents, which are then used to generate the response to the initial conditions, the boundaries of instability, and the characteristics of the free response solution of the system. The coexistence phenomenon, in which the instability wedges disappear and the transition curves overlap, is recovered by this approach, and its features and robustness are examined.

## 1 Introduction

Parametric excitation occurs when a system has time-varying coefficients that drive the oscillation. Dynamical systems that are exposed to parametric excitation experience notable behavior . Parametric excitation happens in various types of systems. For example, in electrical systems, an inductor-capacitor (LC) circuit with periodic time varying capacitance represents Hill’s equation . In ecological systems, the population dynamics in periodic environments causes periodic variability in systems . In mechanical systems, a pendulum with periodic vertical base excitation [7,8], gear whine with cyclic contact points , and horizontal-axis wind-turbine blades with cyclic gravitational forces  have parametric stiffness.

Our previous study on vertical-axis wind-turbine blades (VAWTs) suggested the existence of periodic damping in the equations of motion due to the aeroelastic forces involving a cyclic angle of attack . A simplified mechanism for this can be seen in the context of an airfoil of chord length c. The lift force per unit length F is proportional to the speed v squared through a coefficient, which is assumed to be linear in the angle of attack α, such that $F=12Cd(α)ρcv2=cdαv2$. In a VAWT under ideal steady wind conditions with wind speed u, the angle of attack varies with the cyclic angle $α0(t)≅a0+a1cosωt$ of ambient wind and also with the velocity $x˙$ of the blade deflection (Fig. 1).

If $α≅α0(t)−x˙/u$ and $v≅u−bx˙$, then the force takes the form $F≅cda0u2+cdu2a1cosωt$$c0x˙−c1cosωtx˙$. The first terms are direct excitation, and the latter terms involve cyclic damping. (In an ideal VAWT, the magnitude of the relative flow velocity also varies with the VAWT rotation. Other aeroelastic effects might also affect the mean damping.)

Systems with parametric excitation have been studied broadly with various methods. Cantero et al.  looked into the parametric excitation of mooring lines. They used the method of harmonic balance to obtain an analytical expression of the unstable layout. Lilien and Da Costa  studied the vibration amplitudes caused by parametric excitation of nonlinear cable-stayed structures. They defined the nondimensional parameters for the cable-stayed bridge and then used the harmonic balance method to find the limit-cycle amplitude. Lenci et al.  applied a perturbation method to a parametrically driven pendulum to determine the stability of the harmonic solution. Hsu  examined parametric excitation of a dynamic system with multiple degrees-of-freedom. He combined the method of variation of parameters and the series expansion of the perturbation method to derive the instability criteria for the instability regions.

Besides all the declared methods, the Floquet theory applies to the linear differential equations with parametric coefficients, and it is a very useful tool for studying the stability of dynamic systems [16,17]. However, Hsu explains that in the Floquet theory, approximation is often required to discover the stability characteristics. Furthermore, most of the approximate methods are only valid when the magnitude of the excitation is small [4,18].

The Mathieu equation is one of the most well-known differential equations with parametric excitation. There have been extensive studies on the damped and undamped Mathieu equation using harmonic balance  and perturbation methods . Ward  determined the region of stability of the Mathieu equation by assuming the solution to be consist of only a periodic part. Acar and Feeny  studied the Floquet-based analysis of the Mathieu equation’s responses and assumed a Floquet solution including both exponential and periodic parts. They combined the Floquet theory with the harmonic balance method and investigated the time response and the stability of the equations of motion.

In some parametric systems, a phenomenon called coexistence happens. Coexistence occurs when the boundaries of an unstable region overlap and the instability wedge disappears. Recktenwald and Rand  studied coexistence and obtained the conditions for coexistence to occur in a generalized Ince’s equation.

In this study, we study the response characteristics of a linear oscillator with parametric damping. We first use the Floquet theory to draw conclusions about the effects of cyclic damping on the Floquet multipliers. We then use the original Floquet solution with a non-zero exponential part and a periodic part, which is approximated as a truncated Fourier series. As in Ref. , the harmonic balance method is applied to obtain the Floquet exponents and Fourier coefficients. This analysis provides not only the stability information but also the general response of the system in each stability situation. We uncover the coexistence phenomenon and study its characteristics and then observe the effects of a specific perturbation in the damping.

## 2 Floquet Theory

In this section, we follow a standard development of the Floquet theory  and make adaptions for our system. The Floquet theory is applicable to a linear differential system of the form
$x˙=A(t)x(1)$
(1)
where A(t) = A(t + T) is an n × n periodic matrix of period T and xRn is an n × 1 column vector . Suppose
$x(1)(t):isthesolutionvectorwithI.C.sx0(1)=[10…0]T⋮x(n)(t):isthesolutionvectorwithI.C.sx0(n)=[00…1]T$
where I.C. refers to the initial condition. For a general initial condition, x0 = [c1c2cn]T, the solution becomes
$x(t)=[x(1)(t)x(2)(t)…x(n)(t)](c1c2⋮cn)=X(t)x0$
(2)
where X(t) is called the fundamental solution matrix (FSM). The FSM is the solution vectors of Eq. (1) when initial conditions are unity, such that X(0) = I. Any other solution matrix Z(t) may be written in the following form
$Z(t)=X(t)C$
(3)
where C is a nonsingular n × n matrix. In that case, each of the columns of Z(t) may be written as a linear combination of the columns of X(t). Note that Z(0) = X(0)C = IC = C. Let Z(t) = X(t + T). Then, Z(0) = C = X(T) ≠ I, and from Eq. (3), we have
$X(t+T)=X(t)C$
(4)
This implies that X(t) and X(t + T) are related by the constant matrix C = X(T). By successively iterating, we get
$X(t+mT)=X(t)Cm$
(5)
Equation (4) is a difference equation. If $λi$, the eigenvalues of C are distinct, then Eq. (4) can be diagonalized via the transformation X = YP−1, where the columns of matrix P are eigenvectors of C, and J = P−1CP is a diagonal matrix of eigenvalues. The diagonalized difference equation becomes
$yi(t+T)=λiyi(t),i=1,…,n$
(6)
In this setting, yi are diagonal elements of Y.
Analysis Based on Floquet Theory. A second-order linear differential equation with time-varying coefficients can be expressed in a general form as follows:
$x¨+g(t)x˙+h(t)x=f(t)$
(7)
where g(t) and h(t) are periodic functions. If g(t) = 0 and f(t) = 0, the equation reduces to the Hill’s equation , $x¨+h(t)x=0$. If $h(t)=δ+ϵcosωt$, where $ω$ is the frequency ratio (dimensionless frequency), and $δ$ and $ϵ$ are constants, Hill’s equation takes the form of Mathieu’s equation as follows:
$x¨+(δ+ϵcosωt)x=0$
(8)
Here, we study the application of the Floquet theory to the generalized homogeneous differential equation with parametric damping and stiffness, i.e., Eq. (7) when f(t) = 0. In the state-space form, the equation with cyclic damping and cyclic stiffness appears as follows:
$x˙=(x˙x¨)=[01−h(t)−g(t)](xx˙)=A(t)x$
(9)
where A(t) is a periodic matrix with period $T=2π/ω$. The fundamental matrix solution in Eq. (2) takes the form
$X(t)=[x1(1)(t)x1(2)(t)x2(1)(t)x2(2)(t)]$
(10)
where x1 = x and $x2=x˙$. $W(t)=detX(t)$ = $x1(1)(t)x2(2)(t)−x1(2)(t)x2(1)(t)$ is defined as the Wronskian. Taking the time derivative of W and using Eq. (9), $W˙(t)$ in terms of x1(t) and x2(t) is (for simplicity we drop (t))
$W˙=x˙1(1)x2(2)+x1(1)x˙2(2)−x˙1(2)x2(1)−x1(2)x˙2(1)=−x2(1)x2(2)+x1(1)(−hx1(2)−gx2(2))+x2(2)x2(1)−x1(2)(−hx1(1)−gx2(1))=−h(x1(1)x1(2)−x1(2)x1(1))−g(x1(1)x2(2)−x1(2)x2(1))$
(11)
Unlike with Hill’s equation [16,30], Eq. (12) reduces to
$W˙(t)=−g(t)W(t)$
(12)
By decomposing g(t) into a mean and a time-varying part as $g(t)=c0+g~(t)$ and then solving Eq. (11), we obtain $W(t)=e−c0t−G(t)+k$, where $G(t)=∫g~(t)dt$ and k is an integration constant. At t = 0, $W(0)=e−G(0)+k$ = $detX(0)=detI=1$. Therefore, k = G(0). By periodicity, G(T) = G(0), and therefore,
$W(T)=e−2πc0/ω$
(13)
We have observed that X(T) = X(0)C = C. $λi$, the eigenvalues of C, are the Floquet multipliers, and they satisfy $|C−λI|=0$, i.e., the characteristic equation
$λ2−bλ+d=0$
(14)
where b = trace(C) and $d=det(C)=W(T)$. Different values of b and d lead to real and complex eigenvalues. For the case that $λ$ is real,
$λ=b±b2−4d2$
(15)
and for the case that $λ$ is complex,
$λ=b±i4d−b22|λ|=12(b2+4d−b2)1/2=d$
(16)
Using Eq. (13),
$|λ|=W(T)=e−2πc0/ω=e−πc0/ω$
(17)

Thus, for complex $λ$, Eq. (16) describes a circle of radius $r=|λ|=e−πc0/ω$ (Fig. 2). This holds for any system of the form of Eq. (7), including the damped Mathieu equation, where the mean damping is c0. If c0 ≥ 0, the circle implies stable decay with (typically) quasiperiodic oscillation. If c0 ≥ 0, the system stability depends on $λ$ as follows:

• $|λi|≤1$ stable

• $|λ1|≤1,λ2=±1$ stability transition

• $|λ1|<1$, $|λ2|>1$ unstable

When $λ=1$, the system has an underlying periodicity of period T, since
$y(i)(t+T)=λiy(i)(t)=y(i)(t)$
(18)
for one of the y(i). In such case, Eq. (14) implies b = 1 + d, and from Eqs. (16) and (17),
$b=1+e−2πc0/ω$
(19)
When $λ=−1$, the underlying periodic solution is associated with period 2T, since one of the y(i) obeys
$y(i)(t+T)=−y(i)(t)$
(20)
Also $b=−(1+e−2πc0/ω)$. When b = 0, $λ=±id$.

Note that the transition between stable and unstable solutions of the differential equation (Eq. (7)) involves a solution with underlying periodicity T$(λ=1)$ or 2T$(λ=−1)$.

Equation (6) is a first-order linear difference equation, and if $λi$ are distinct, we seek a solution of the form
$y(i)(t)=y0(i)λist$
(21)
where s is an unknown constant. Substituting Eq. (21) into Eq. (6)
$y(i)(t+T)=y0(i)λis(t+T)=y0(i)λistλisT=λiy(i)(t)$
(22)
Equation (22) is satisfied if s = 1/T . Then,
$y(i)(t)=y0(i)λit/T=y0(i)eμit$
(23)
where $μi=lnλi/T$. The $λi$ were defined earlier as the Floquet multipliers, and μi are the associated Floquet exponents. When all $|λi|≤1$, then all Re(μi) ≤ 0, and as discussed earlier, the system is stable. $W(T)=detC$ indicates the growth of the state-space volume associated with iterations of the Poincare map . When c0 = 0, W(T) = 1, and the map is area preserving. When c0 > 0, the map is contracting and dissipative. When c0 < 0 and both multipliers are complex with $|λi|>1$, the system is unstable. Furthermore, when c0 < 0, and the multipliers are real, at least one multiplier will be such that $|λi|>1$ and the system is unstable. Thus, parametric excitation cannot stabilize a system with negative mean damping.
Special Case Whenc0 = 0. When c0 = 0, the circle has radius 1. So when $λ=±1$, the eigenvalues are repeated and the Jordan form is not diagonal. However, the conclusions about stability transitions hold. Based on Eq. (13), when c0 = 0, W(T) = 1, and accordingly, the characteristic equation for the differential equation with parametric damping found by Floquet theory is $λ2−trace(C)λ+1=0$. In such circumstances, the Floquet circle of periodic damping appears to have the same form as that of Hill’s equation when stiffness is periodic, namely,
$λ=a±a2−42,a=trace(C)$
(24)
although a = trace(C) is not the same as in Hill’s equation.

## 3 Oscillator With Parametric Damping

The concern of this study is to inquire into a linear differential equation for which, unlike Mathieu’s equation, stiffness is constant but damping is parametric, such that
$x″+(d0+d1cosω^τ)x′+ωn2x=0$
(25)
where d0 and d1 are constants and $()′=d()/dτ$. By performing a change of variables, $t=ωnτ$, the alternative equation is
$x¨+(c0+c1cosωt)x˙+x=0$
(26)
where $ω=ω^/ωn$ is a dimensionless frequency ratio, ($˙$)=d()/d t, and c0 and c1 are constants. Defining a state vector $x=[xx˙]T$, Eq. (26) can be expressed in the state-space form as $x˙=A(t)x$. Therefore, the Floquet theory can be applied to the intended differential equation [16,22].

Hartono and Burgh  studied the equation with a time-varying damping coefficient with zero mean damping and investigated the stability. They transformed variables and reduced the equation to a standard Hill’s equation $y¨+[Λ+Q(t)]y=0$, where Q(t) is a periodic function with period $π$. They approximated Λs for which the stability transitions occur. Later, we will discuss Hartono and Burgh’s results and compare with our observed results.

### 3.1 Approximate Solution.

The Floquet theory is combined with harmonic balance to find an approximate solution to Eq. (26). Following the Floquet theory, the solution consists of an exponential term and periodic term, x(t) = eμtp(t), where p(t) is periodic with period T. A truncated series solution is assumed such that
$x(t)=eμt∑j=−n+najeijωt$
(27)
where μ is the Floquet exponent and ajs are Fourier coefficients. The assumed solution is plugged into Eq. (26) and harmonic balance is applied to find the coefficients of $eijωt$. The harmonic balance leads to an equation of the form D(μ)a = 0, where D(μ) is a (2n + 1) × (2n + 1) coefficient matrix and a is a (2n + 1) × (1) vector of the Fourier coefficients (aj’s) for the assumed solution. To have nontrivial and nonzero aj’s, the determinant of the coefficient matrix, D(μ), is required to be zero. Setting the determinant to zero gives the characteristic equation with unknown μ.
The characteristic equation of |D(μ)| = 0 is a (4n + 2) degree polynomial in terms of μ, which yields (4n + 2) roots for μ. From the Floquet theory, the number of μ’s should be equal to the dimension of the state space, N = 2, so it turns out that some of the μ’s of the characteristic equation are spurious. Thus, N = 2 of the roots are independent and are called the principle roots. According to the Floquet theory , the characteristic multipliers satisfy
$λ1λ2⋯λN=exp(∫0TTr(A(s))ds)$
(28)
where $λk=eμkT$. For this problem, two distinct Floquet exponents (μ’s) satisfy Eq. (28) such that
$e(μ1+μ2)T=exp(∫0TTr(A(s))ds)$
(29)
In Eq. (9), $Tr(A)=−g(t)=−(c0+c1cosωt),$ and so, $∫0TTr(A(s))ds=−c0T$. Therefore, we can write $e(μ1+μ2)T=e(−c0T±2πik)$, since any multiple of 2π i can be added to an exponential argument. Using $T=2π/ω$, this implies
$μ1+μ2=−c0±kωi$
(30)

Equation (30) specifies the relationship between the two principle Floquet exponents, μ1 and μ2, with the mean damping term c0, the time period, $T=2π/ω$, and an integer k. The solution to |D(μ)| = 0 results in 4n + 2 roots for μ and the two individual principal roots for any pair of c0 and $ω$ are selected to satisfy Eq. (30).

In the assumed solution (Eq. (27)), the imaginary part of the Floquet exponent, Im(μ), combines with the frequencies of the harmonics of the periodic term to define the frequencies of the response. The general response is a linear combination of principal responses such that
$x(t)=A1eμ1t∑j=−n+najeijωt+A2eμ2t∑j=−n+najeijωt$
(31)
where A1 and A2 are determined by initial conditions. In Sec. 3.2, the truncated solution of Eq. (26) is explained.

### 3.2 Truncated Series Solution.

The truncated series solution of the parametric excited motion can be used to quantify the stability domain and approximate the response to the initial conditions and response frequencies. We consider the truncated solution in Eq. (27) for n = 2 as follows:
$x(t)=eμt∑j=−2+2ajeijωt$
(32)

Previously, Acar and Feeny  showed that for n = 2, the truncated solution of the Mathieu equation converges well except when $ω$ was small.

Substituting the approximate truncated solution into the governing equation (Eq. (26)) and balancing coefficients of $eijωt$ lead to the following equation:
$[B−2c12000c12B−1c12000c12B0c12000c12B1c12000c12B2](a−2a−1a0a1a2)=(00000)$
(33)
where
$B−2=ic0μ−2ic0ω−μ2+4μω−4ω2+1B−1=ic0μ−ic0ω−μ2+2μω−ω2+1B0=12(2ic0μ−2μ2+2)B1=ic0μ+ic0ω−μ2−2μω−ω2+1B2=ic0μ+2ic0ω−μ2−4μω−4ω2+1$
(34)
Equation (33) specifies an eigenvalue problem. The characteristic equation is obtained by setting the determinant of the coefficient matrix to be zero. For a set of parameters, when n = 2, there will be ten μs, yet there exist only two principle roots. The four μs with the smallest real parts are identified among all of the μs. We choose one of these as μ1 and then apply Eq. (30) with k = 0 to select the second principle μ, such that μ2 = −c0μ1. Letting $μ=α+iβ$,
$x1(t)=e(α1+iβ1)t∑j=−2+2ajeijωt=eα1t∑j=−2+2ajei(jω+β1)tx2(t)=e(α2+iβ2)t∑j=−2+2ajeijωt=eα2t∑j=−2+2ajei(jω+β2)t$
(35)
Then, the total response is x(t) = A1x1(t) + A2x2(t). Real parts of μs govern the exponential growth or decay of the solution. If any Re(μ) > 0, the solution grows exponentially (unstable). The truncated response has terms of frequencies $jω+β1$ and $jω+β2$, $j=−n,…,n$. If $β1$ and $β2$ are accurate, the full set of response frequencies will cover $j=−∞,…,∞$. Equation (30) with k = 0 implies $β2=−β1$ for our chosen principal roots.

## 4 Results

The analysis of the parametric excitation with periodic damping without and with the mean damping term (when c0 = 0 and c0 ≠ 0, respectively) are presented in this section.

### 4.1 Zero Mean Damping, $c0=0$

#### 4.1.1 Stability Boundaries.

As explained in Sec. 3.2, when the real parts of μ1 and μ2 are negative or zero the solution is stable, whereas for any μ with a positive real part, the solution grows exponentially. Figure 3 demonstrates the stable and unstable regions for three different truncation orders, n = 2, n = 6 and n = 10 in the $(ω,c1)$ plane, where $ω$ and c1 are the frequency and amplitude of the harmonic damping. In this range, increasing the number of terms in the Fourier expansion provides a more precise solution. Increasing n = 2 to n = 10 in the truncated solution enhances the transition curves. Inside the tongues of instability (white zone), the response is unstable with underlying periodicity of period 2T, whereas outside of the tongues (shaded zone), the response is quasi-periodic (to be discussed in more detail later).

In Fig. 3, the truncated solution captures the subharmonic instability at $ω=2$ and superharmonic instability wedges based at $ω=2/3$ and $ω=2/5$. The existence of the wedge based at $ω=2/3$ in Fig. 3 persists with the increasing n. The wedge is not observed to reach the $ω$ axis because of limited resolution of the plot, which was generated by plotting pixels. Other superharmonic tongues are not detected, and a wedge based at $ω≈0$ and its connection to the subharmonic wedge is erroneous. An improved description for low $ω$ and large c1 would require higher values of n. As shown in Fig. 3, when n = 10, the higher truncation order reveals more and finer wedges of superharmonic instability.

Hartono and Burgh  studied a cyclically damped system (with c0 = 0) based on a cyclic damping equation of the form $x″+ϵcos2tx′+Λx=0$. They plotted transition curves as $ϵ$ versus Λ, the square of the nondimensionalized natural frequency, when the excitation frequency was fixed at 2. The relationships between the parameters in Ref.  and those in our model are $ϵ=2c1/ω$ and $Λ=4/ω2$. The results in Ref.  for cyclic damping show instability wedges based at Λ = 1 and 9 (and generally values of m2 for odd values of m). The wedges based at $Λ=4,16,…$ (generally m2 for even values of m), which are known to occur in the Mathieu equation, are collapsed, such that the boundaries overlap into a single curve and instability does not occur (associated with coexistence, discussed later). All of these wedge base points correspond to $ω=2,1,2/3,$ and 1/2 in Fig. 3, for which wedges based at 2 and 2/3 are not coexistent and are observed in convergent plots of Fig. 3. Thus, the two studies are consistent.

#### 4.1.2 Response Frequency Analysis.

The classical Floquet-based analysis describes the properties of Eq. (1) such as stability boundaries , while the truncated solution method provides the characteristics of the response in addition to the tongues of instability. The approximate solution contains an exponential and a periodic part. A combination of Floquet exponents and excitation harmonics comprises the frequencies of the response. The response frequencies in the truncated solution are $Im(μ)+jω$, where j = ±1, ± 2, …, ± n and n is the number of assumed harmonics.

Figure 4 shows the response frequencies of the truncated solution as functions of $ω$ for c0 = 0, c1 = 1, and the truncation orders n = 2 and n = 4. At the instability wedges based at $ω=2/m$, m = 1, 3, 5, …, the Fourier expansion of the periodic part of the response has harmonics of $mω/2$, m = 1, 3, 5, …. In instability wedges based at $ω=2/m$, m = 2, 4, 6, …, the Fourier expansion of the periodic part of the response has harmonics of $mω/2$, m = 2, 4, 6, …. For example, near $ω=2,$ the response frequencies include $ω/2,3ω/2$, and $5ω/2$, and at the crossing, near $ω=1,$ the response frequencies include $ω,2ω,$ and $3ω$.

With varying $ω$, pairs of branches may merge into single branches or cross one another at isolated points. Intervals where two frequency branches have merged into one branch are associated with instability zones in the stability diagram and with the stability transition points in the Floquet circle (Fig. 2, case of c0 = 0). In the intervals of instability Re(μ1,2) = ±α and $Im(μ1,2)=±β$. Merging of two branches into one branch indicates a stability transition, and the interval of merging represents the instability tongues. Branches crossing at one point indicate that there is no frequency-merging interval, such that would-be wedge boundaries overlap and the periodic responses on the boundaries coexist, and therefore, the coexistence phenomenon is happening.

#### 4.1.3 Growth/Decay Factor.

In the approximate solution, the real part of the Floquet exponent specifies the exponential growth. Figure 5 shows Re(μ) as a function of $ω$ when c0 = 0. The graph shows two intervals where the Floquet exponents have two different values for the real parts, forming “bubbles” in the graph (a tiny bubble near $ω≈0.61$, and a large bubble centered at $ω=2$). Re(μ) > 0 denotes growth in the truncated solution and therefore yields instability. In the case of stability, the least negative value governs the decay rate. Therefore, the diagram indicates the growth and decay factors of the response. While Fig. 5(a) demonstrates an instability near $ω=2$, $ω=0.61,$ and $ω<0.38$ for c0 = 0, numerical studies suggest that the latter case is a truncation error, and there are indeed mostly stable solutions for $ω<0.38$. For low values of $ω,$ the lower truncation is not sufficient to describe the vibration solution or stability. This erroneous range of $ω$ is reduced with the increasing truncation order n. In Fig. 5, two different orders of truncation illustrate the benefit of higher values of n in approximating the solution for lower frequencies.

#### 4.1.4 Coexistence When c0 = 0.

Coexistence phenomenon occurs when stability transitions overlap and the tongues of instability disappear (black lines in Fig. 6(a)). Figure 6(a) shows that stability tongues collapse into lines based at $ω≈0.91$ and $ω≈0.46$. The thin instability tongue based at $ω≈0.61$ is associated with the small bubble in the growth/decay factor diagram (Fig. 6(a)) and the small frequency-merging interval in the response frequency plot (Fig. 6(b)). Figure 6(a) shows the properties of the stable and unstable responses, transition boundaries, and coexistence curves, and whether the response is periodic or quasi-periodic. “Unstable periodic” means unstable exponential growth modulates an oscillation that is otherwise periodic (not quasi-periodic). In the wedges, based at $ω=2/m$, $m=1,3,5,…$, this periodic part has period 2T. At coexistence, the response is periodic with period T.

In the growth/decay factor diagram of Fig. 6(a), the bubbles, normally associated with the instability zones, collapse, and thus no bubble appears at coexistence points. In Fig. 6(b), one can see that coexistence corresponds to the crossing points in the response frequency plot, where periodic solutions cross one another at $ω≈0.91$ and $ω≈0.46$.

### 4.2 Non-Zero Mean Damping $c0≠0$⁠.

For c0 ≠ 0, the stability boundaries and each of the response features are affected by mean damping. To determine the effect of damping, the previous analysis is repeated for the case c0 ≠ 0.

#### 4.2.1 Stability Boundaries.

Figure 7 shows the stability diagrams when c0 = 0.05 (star symbols, slightly lifted in the wedges) for three different truncation orders and compares with the stability boundaries when c0 = 0 (dots). As noted earlier, a higher truncation number produces more accurate stability curves within lower frequency ranges. It can be perceived how the stability boundaries are affected by damping. A small addition of mean damping smoothens the wedge tips and shrinks the instability zones.

#### 4.2.2 Response Frequency Analysis.

Figure 8 zooms into the response frequency plot when c0 = 0.05. In Fig. 6(b), a frequency crossing at $ω≈0.91$ indicated the coexistence phenomenon for c0 = 0. However, for c0 ≠ 0, the crossing turns into an interval and the coexistence line disappears. For the lower frequency, $ω≈0.46$, at which we had coexistence with c0 = 0, the truncation order and resolution are not sufficient to clearly distinguish whether c0 ≠ 0 produces an interval.

#### 4.2.3 Growth/Decay Factor.

Figure 9 illustrates that the added mean damping shifts the growth/decay diagram downward, and the shift is c0/2, which is consistent with Eq. (30). In Fig. 9, at $ω≈0.91,$ a small bubble pops up at the location where coexistence happens (no bubble for c0 = 0). One might expect the coexistence line in Fig. 6(a) to separate and form an instability wedge when the merging interval of Im(μ) appears. But since the bubble occurs with values of Re(μ1) < Re(μ2) < 0, the response stays stable. Thus, no instability tongues appear near $ω≈0.91$ in Fig. 7, although the crossing point has turned into an interval in Fig. 8.

By increasing the amplitude of the parametric damping (c1), the bubble in Re(μ) values in Fig. 9 grows but does not cross zero and the system stays stable. As $ω$ increases and the small bubble is reached in Fig. 9, two complex multipliers $λ$ (recall $λ=eμT$) on the circle in Fig. 2 collide on the positive real axis. As $ω$ continues to increase, $λ1$ and $λ2$ split on the real axis and then come back together and collide again in Fig. 2, as the far end of the bubble is reached in Fig. 9. When $ω$ increases beyond the bubble, $λ1$ and $λ2$ are on the circle and are complex. For values of $ω$ in the bubble, $λ1$ and $λ2$ are real with $|λ1|<1$ and $|λ2|<1$ (Re(μ1) < 0 and Re(μ2) < 0)) and the system is stable (in this example). The response has exponential decay that modulates an oscillation that is otherwise of period T.

#### 4.2.4 No Coexistence When c0 ≠ 0.

The previous discussion indicates that by adding a mean damping to the system, the coexistence is disrupted. The frequency-branch crossings are perturbed into frequency-branch mergers over small intervals. With the frequency-branch interval, a bubble forms in the graph of the Re(μ), which was not a characteristic of coexistence. However, the Re(μ) is pushed to be negative (in the example shown), and so, the bubble may not generate an instability although it is conceivable that instability could happen if c1 is increased enough. The overlapping instability transition boundaries, which are seen in coexistence, separate with the emergence of the frequency-branch interval, but are no longer serving as stability boundaries since the Re(μ) has been reduced by c0/2.

While coexistence is disrupted by the addition of slight damping (or probably other perturbations), it is important as a limiting case, since it points to a near description of behavior for the case of very small mean damping.

The importance of studying coexistence comes into sight when the system parameters are positioned on the coexistence line. It is possible that arbitrary perturbations may open up the coexistence line and generate instabilities and the response may fall into unstable zone. The shaded area in the stability diagrams, where the system parameters are located in the “safe” stability zone, indicate the best possible parameters for design. As such, parameter values near the coexistence line may be avoided if instability is undesirable.

### 4.3 Response of Undamped and Damped Systems.

Depending on the parameters, the system may be stable, neutrally stable, or unstable. Three sets of parameters, associated with three types of stability, are chosen from the stability diagram to obtain the truncated solution of Eq. (26). Figures 10 and 11 illustrate both the numerical and the truncated solutions for the sets of parameters when c0 = 0 and c0 = 0.05, respectively. The truncated solution with n = 3 truncation order agrees with the numerical analysis for the chosen sets of parameters.

Figure 10(a) represents a quasi-periodic bounded response. In this case, the quasi-periodic drift occurs on a long time scale and is only slightly detectable in the plot. Figure 10(b), for which $ω=ω1=0.91,$ is very close to coexistence and therefore looks periodic with period $T1=2π/ω1≈6.904$. Figure 10(c) shows a subharmonic unstable solution at a periodicity of $2T2=4π/ω2≈6.28$ for which $ω=ω2=2$ (T2 is the period of excitation).

Figures 11(a) and 11(b) show exponentially stable solutions with frequency characteristics similar to those in Figs. 10(a) and 10(b) and additional mean damping. Figure 11(c) is unstable despite the mean damping. A very slight decrease in the growth rate compared to Fig. 10(c) can be detected on careful inspection. Putting the mean damping into effect for the same system parameters (c0 ≠ 0) dampens the response with little effect on frequency content.

The fast Fourier transforms (FFTs) of the truncated solutions are compared with those of the numerical solutions in Figs. 12 and 13. The analytical truncation order is n = 3. The frequency content of analytical and numerical solutions for the chosen excitation frequencies are consistent. However, for small values of $ω,$ the truncated and numerical solutions do not agree. Increasing the truncation order improves the accuracy of the results compared with the numerical solution.

## 5 Conclusion

We studied the dynamics of an oscillator with cyclic damping, which was motivated by vertical-axis wind-turbine blade dynamics. Specifically, we sought an approximate analytical solution to a linear second-order differential equation with a periodic damping coefficient and observed the characteristics of the response.

We sought a general Floquet solution in which the periodic part was approximated with a truncated Fourier series. Applying harmonic balance led to a relationship among the Floquet exponents, Fourier coefficients, and the parameters.

The exponents were then used to formulate the response to the initial conditions, the response frequency content, and the stability characteristics. The theoretical responses of the system and their FFTs were compared with the numerical solutions, and the consistency of the responses for different sets of parameters was shown. The analysis is good if $ω$ is not too small and can be improved by increasing the analytical truncation order n. Limits on the truncation order were not investigated.

For the case when there was zero mean damping, we captured a phenomenon called “coexistence.” Coexistence has been recognized to occur when stability boundaries overlap such that the instability wedges disappear. This is accompanied by a crossing of response frequencies as the frequency of excitation changes and an absence of a “bubble” in the real part of the Floquet exponents (meaning the real part of the Floquet exponents is repeated). We found that the addition of mean damping disrupted the features of coexistence such that the response frequencies merged over a parameter interval, in which the real parts of the Floquet exponents became distinct (forming a bubble in the plot of Re(μ) versus $ω$). In such case, an unstable wedge may not appear in the parameter range at hand since the mean damping can stabilize the responses in very slender wedges.

The results of this study may have relevance to VAWTs. The presence of cyclic damping in VAWTs, especially if aeroelastic effects reduce the effective mean damping, suggests the potential for instabilities of the same pattern as is familiar with cyclic stiffness, namely, at excitation frequencies near $2ωn/k,k=1,2,…$. In the case of cyclic damping, the even values of k represent coexistence in the ideal model, and the system does not actually destabilize. However, perturbations of the model could destabilize the coexistence events, and thus, designs should avoid these frequencies.

## Acknowledgment

This work is based on a project supported by the National Science Foundation, under grant CMMI-1435126. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

## Conflict of Interest

There are no conflicts of interest.

## References

1.
Mandel’shtam
,
L.
, and
Papaleksi
,
N.
,
1934
, “
On the Parametric Excitation of Electric Oscillations
,”
Zhurnal teknicheskoy fiziki
,
4
(
1
), pp.
1
47
.
2.
Dugundji
,
J.
, and
,
V.
,
1973
, “
Lateral Bending-Torsion Vibrations of a Thin Beam Under Parametric Excitation
,”
ASME J. Appl. Mech.
,
40
(
3
), pp.
693
698
. 10.1115/1.3423075
3.
Zhang
,
W.-M.
, and
Meng
,
G.
,
2007
, “
Nonlinear Dynamic Analysis of Electrostatically Actuated Resonant MEMS Sensors Under Parametric Excitation
,”
IEEE Sens. J.
,
7
(
3
), pp.
370
380
. 10.1109/JSEN.2006.890158
4.
Hsu
,
C. S.
,
1972
, “
Impulsive Parametric Excitation: Theory
,”
ASME J. Appl. Mech.
,
39
(
2
), pp.
551
558
. 10.1115/1.3422715
5.
Kniffka
,
T. J.
, and
Ecker
,
H.
,
2015
, “
Parametererregte Mikroelektromechanische Systeme (MEMS)
,”
e & i Elektrotechnik und Informationstechnik
,
132
(
8
), pp.
456
461
. 10.1007/s00502-015-0372-8
6.
Klausmeier
,
C. A.
,
2008
, “
Floquet Theory: A Useful Tool for Understanding Nonequilibrium Dynamics
,”
Theoretical Ecology
,
1
(
3
), pp.
153
161
. 10.1007/s12080-008-0016-2
7.
Liang
,
Y.
, and
Feeny
,
B.
,
2008
, “
Parametric Identification of a Chaotic Base-Excited Double Pendulum Experiment
,”
Nonlinear Dyn.
,
52
(
1–2
), pp.
181
197
. 10.1007/s11071-007-9270-x
8.
Sartorelli
,
J. C.
, and
Lacarbonara
,
W.
,
2012
, “
Parametric Resonances in a Base-Excited Double Pendulum
,”
Nonlinear Dyn.
,
69
(
4
), pp.
1679
1692
. 10.1007/s11071-012-0378-2
9.
Perret-Liaudet
,
J.
,
Carbonelli
,
A.
,
Rigaud
,
E.
,
Nelain
,
B.
,
Bouvet
,
P.
, and
Vialonga
,
C. J.
,
2014
, “
Modeling of Gearbox Whining Noise
,”
Technical Report, SAE Technical Paper, Society of Automotive Engineers
.
10.
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
. 10.1115/1.4033341
11.
Afzali
,
F.
,
Kapucu
,
O.
, and
Feeny
,
B. F.
,
2016
, “
Vibrational Analysis of Vertical-Axis Wind-Turbine Blades
,”
Proceedings of the ASME 2016 International Design Engineering Technical Conferences
,
Charlotte, NC
,
Aug. 21–24
,
Paper No. IDETC2016-60374
.
12.
Cantero
,
D.
,
Rønnquist
,
A.
, and
Naess
,
A.
,
2016
, “
Recent Studies of Parametrically Excited Mooring Cables for Submerged Floating Tunnels
,”
Procedia. Eng.
,
166
, pp.
99
106
. 10.1016/j.proeng.2016.11.571
13.
Lilien
,
J.-L.
, and
Da Costa
,
A. P.
,
1994
, “
Vibration Amplitudes Caused by Parametric Excitation of Cable Stayed Structures
,”
J. Sound. Vib.
,
174
(
1
), pp.
69
90
. 10.1006/jsvi.1994.1261
14.
Lenci
,
S.
,
Pavlovskaia
,
E.
,
Rega
,
G.
, and
Wiercigroch
,
M.
,
2008
, “
Rotating Solutions and Stability of Parametric Pendulum by Perturbation Method
,”
J. Sound. Vib.
,
310
(
1–2
), pp.
243
259
. 10.1016/j.jsv.2007.07.069
15.
Hsu
,
C.
,
1963
, “
On the Parametric Excitation of a Dynamic System Having Multiple Degrees of Freedom
,”
ASME J. Appl. Mech.
,
30
(
3
), pp.
367
372
. 10.1115/1.3636563
16.
Rand
,
R. H.
,
2012
, “
Lecture Notes on Nonlinear Vibrations
,” https://ecommons.cornell.edu/handle/1813/28989, Accessed 30 September, 2020.
17.
Cesari
,
L.
,
2012
,
Asymptotic Behavior and Stability Problems in Ordinary Differential Equations
, Vol.
16
,
,
New York
.
18.
Hsu
,
C. S.
,
1965
, “
Further Results on Parametric Excitation of a Dynamic System
,”
ASME J. Appl. Mech.
,
32
(
2
), pp.
373
377
. 10.1115/1.3625809
19.
Ramakrishnan
,
V.
, and
Feeny
,
B. F.
,
2012
, “
Resonances of a Forced Mathieu Equation With Reference to Wind Turbine Blades
,”
ASME J. Vib. Acoust.
,
134
(
6
), p.
064501
. 10.1115/1.4006183
20.
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
. 10.1115/1.4005829
21.
,
J. F.
,
Miller
,
N. J.
,
Shaw
,
S. W.
, and
Feeny
,
B. F.
,
2008
, “
Mechanical Domain Parametric Amplification
,”
ASME J. Vib. Acoust.
,
130
(
6
), p.
061006
. 10.1115/1.2980382
22.
Nayfeh
,
A. H.
, and
Mook
,
D. T.
,
2008
,
Nonlinear Oscillations
,
John Wiley & Sons
,
New York
.
23.
Taylor
,
J. H.
, and
Narendra
,
K. S.
,
1969
, “
Stability Regions for the Damped Mathieu Equation
,”
SIAM J. Appl. Math.
,
17
(
2
), pp.
343
352
. 10.1137/0117033
24.
Younesian
,
D.
,
,
E.
, and
Sedaghati
,
R.
,
2005
, “
Existence of Periodic Solutions for the Generalized Form of Mathieu Equation
,”
Nonlinear Dyn.
,
39
(
4
), pp.
335
348
. 10.1007/s11071-005-4338-y
25.
Thomson
,
W.
,
1996
,
Theory of Vibration With Applications
,
CRC Press
,
Englewood Cliffs
.
26.
Nayfeh
,
A. H.
,
2008
,
Perturbation Methods
,
John Wiley & Sons
,
New York
.
27.
Benaroya
,
H.
, and
Nagurka
,
M. L.
,
2011
,
Mechanical Vibration: Analysis, Uncertainties, and Control
,
CRC Press
,
Englewood Cliffs
.
28.
Turrittin
,
H.
,
1952
, “Asymptotic Expansions of Solutions of Systems of Ordinary Linear Differential Equations Containing a Parameter,”
Contributions to the Theory of Nonlinear Oscillations
, Vol.
II
,
Princeton University Press
,
Princeton, NJ
, pp.
81
116
.
29.
Rand
,
R. H.
,
1969
, “
On the Stability of Hill’s Equation With Four Independent Parameters
,”
ASME J. Appl. Mech.
,
36
(
4
), pp.
885
886
. 10.1115/1.3564793
30.
Ward
,
M.
,
2010
, “
Lecture Notes on Basic Floquet Theory
,” http://www.math.ubc.ca/~ward/teaching/m605/every2_floquet1.pdf,
Accessed 30 September, 2020
.
31.
Recktenwald
,
G.
, and
Rand
,
R.
,
2005
, “
Coexistence Phenomenon in Autoparametric Excitation of Two Degree of Freedom Systems
,”
Int. J. Non Lin. Mech.
,
40
(
9
), pp.
1160
1170
. 10.1016/j.ijnonlinmec.2005.05.001
32.
Magnus
,
W.
, and
Winkler
,
S.
,
1979
,
Hill’s Equation
,
Dover
,
New York
.
33.
Arnold
,
V. I.
,
1978
,
Mathematical Methods of Classical Mechanics
,
Springer-Verlag
,
New York
.
34.
Hartono
, and
Burgh
,
A. H. P.
,
2002
,
An Equation Time-Periodic Damping Coefficient: Stability Diagram and an Application
,
Delft University of Technology
,
Delft, The Netherlands
.
35.
Kuchment
,
P. A.
,
1993
,
Floquet Theory for Partial Differential Equations
, Vol.
60
,
Birkhäuser
,
Basel
.