A Floquet-Based Analysis of Parametric Excitation Through the Damping Coefficient

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 [1–4]. 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 [5]. In ecological systems, the population dynamics in periodic environments causes periodic variability in systems [6]. In mechanical systems, a pendulum with periodic vertical base excitation [7,8], gear whine with cyclic contact points [9], and horizontal-axis wind-turbine blades with cyclic gravitational forces [10] 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 [11]. 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. [12] 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 [13] 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. [14] applied a perturbation method to a parametrically driven pendulum to determine the stability of the harmonic solution. Hsu [15] 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 [16] and perturbation methods [19–29]. Ward [30] determined the region of stability of the Mathieu equation by assuming the solution to be consist of only a periodic part. Acar and Feeny [10] 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 [31] 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. [10], 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.

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 c₀. If c₀ ≥ 0, the circle implies stable decay with (typically) quasiperiodic oscillation. If c₀ ≥ 0, the system stability depends on $λ$ as follows:

• $|λi|≤1$ stable

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

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

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)$⁠.

Hartono and Burgh [34] 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.

Equation (30) specifies the relationship between the two principle Floquet exponents, μ₁ and μ₂, with the mean damping term c₀, 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 c₀ and $ω$ are selected to satisfy Eq. (30).

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

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

As explained in Sec. 3.2, when the real parts of μ₁ and μ₂ 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 c₁ 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 c₁ 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 [34] studied a cyclically damped system (with c₀ = 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. [34] and those in our model are $ϵ=2c1/ω$ and $Λ=4/ω2$⁠. The results in Ref. [34] for cyclic damping show instability wedges based at Λ = 1 and 9 (and generally values of m² for odd values of m). The wedges based at $Λ=4,16,…$ (generally m² 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.

The classical Floquet-based analysis describes the properties of Eq. (1) such as stability boundaries [35], 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 c₀ = 0, c₁ = 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 c₀ = 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.

In the approximate solution, the real part of the Floquet exponent specifies the exponential growth. Figure 5 shows Re(μ) as a function of $ω$ when c₀ = 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 c₀ = 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.

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$⁠.

For c₀ ≠ 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 c₀ ≠ 0.

Figure 7 shows the stability diagrams when c₀ = 0.05 (star symbols, slightly lifted in the wedges) for three different truncation orders and compares with the stability boundaries when c₀ = 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.

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

Figure 9 illustrates that the added mean damping shifts the growth/decay diagram downward, and the shift is c₀/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 c₀ = 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(μ₁) < Re(μ₂) < 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 (c₁), 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(μ₁) < 0 and Re(μ₂) < 0)) and the system is stable (in this example). The response has exponential decay that modulates an oscillation that is otherwise of period T.

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 c₁ 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 c₀/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.

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 c₀ = 0 and c₀ = 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$ (T₂ 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 (c₀ ≠ 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.

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.

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.

There are no conflicts of interest.