Analysis and suppression of a self-excitation vibration via internal stiffness and damping nonlinearity

Abstract

This work concerns the usage of the internal stiffness and damping nonlinearity for vibration suppression in a van der Pol–type mechanical self-excitation system. Changing the vertical or horizontal damping could limit the growth of the flow-induced instability. Applying harmonic balance method, the energy transfer for the self-excitation vibration is investigated through the limit circle. The steady amplitude of the limit circle is a key parameter that could be used to evaluate the effectiveness of the self-excitation oscillation suppressing. Increasing horizontal damping could reduce the rate of roll on for the steady amplitude curve of the limit circle, but the critical flow velocity for the limit circle occurring is minimally affected. Increasing vertical damping could increase the critical flow velocity, but the rate of the roll on is virtually unaffected when the parameters are properly chosen.

Keywords

Nonlinear vibration vibration suppression self-excitation oscillation stiffness and damping nonlinearities limit circle

Introduction

Self-excitation of a structure is a common nonlinear mechanism, which can be caused by steady wind flow¹ and dry friction² or artificially achieved by nonlinear feedback control.^3,4 The action of the self-excitation oscillation is two-fold: one is to operate utilizing self-excitation mechanism for the highest energy efficiency. Some examples are atomic force microscopy,^5,6 sensing,^7–10 and energy harvesting.¹¹ The self-excitation systems presented in these literatures always consist of constant source of energy, oscillator, and control device supplying power to the vibrating system that modifies the flow of energy and directs the vibrating energy.¹² The other considered by researchers is that the self-excitation vibration of the structure needs to be suppressed to enhance the performance.^13,14 This type of the self-excitation vibration has ultrasonic vibration of the tool in the process of cutting, oscillation of the rotating machineries, and flutter of the aerofoil.¹ The main objective of this contribution is that changing internal stiffness and damping nonlinearity suppresses the self-excitation vibration.

Appleton and van der Pol considered the self-sustained vibration of an electrical circuit involving a triode vacuum tube.^15,16 van der Pol described an alternative version of the equation that included in his now famous paper.¹⁵ The negative damping was of the linear viscous type, and the positive viscous damping was proportional to the product of the square of the displacement and the velocity. Self-excitation vibration is generated from the interaction between a negative near-equilibrium damping mechanism that destabilizes the static equilibrium of the system and a nonlinear damping operating far-from-equilibrium that limits the growth of the instability. As a result, the limit circle is obtained and the system is constrained to an autonomous vibration with the frequency very close to the natural frequency of the system. Van der Pol’s paper has attracted the attention of many researchers who have made countless attempts to deal with this equation analytically.^17,18 Recently, modifications have been made to van der Pol’s original equation in terms of the relationship between the damping and stiffness forces to velocity and displacement, respectively.^19–21 Krack et al.²² developed a novel reduced model for the analysis of self-excitation vibration and demonstrated the applicability of the proposed overall methodology that complemented by generalized Fourier–Galerkin method and numerical continuation. Chen¹ studied crosswind fatigue analysis of flow-excited flexible structures around vortex lock-in speed where the nonlinear aerodynamic damping effect is significant. He also derived analytical solutions of circle number, fatigue damage, and the validated accuracy through comparison with rain-flow circle counting method. Malas and Chatterjee³ proposed a nonlinear velocity feedback control method for generating artificial self-excitation vibration in a 2-degrees-of-freedom mechanical system. Kwaśniewki et al.²³ investigated the self-excitation acoustical system in elastic construction stress changing measurement.

Self-excitation vibration is usually undesirable. Literature on controlling the self-excitation vibration is vast.^24–30 Malas and Chatterjee⁴ proposed an acceleration feedback–based technique for the self-excitation vibration in a multi-degrees-of-freedom mechanical system. An adaptive control was proposed to keep the vibration amplitude at the desired level. Chatterjee² introduced a new method of controlling friction-driven self-excitation vibration. Das and Mallik²⁷ proposed proportional–derivative (PD) control via time-delayed velocity and displacement feedback for suppressing friction-induced self-excitation vibration. Chatterjee and Malas²⁸ reported stiffness-switching methods for generating artificial self-excitation vibration in a simple mechanical system. It was found that the on-off fashion can be effectively utilized for controlling self-excitation vibration in a simple mechanical system.

It is significant to suppress self-excitation vibration passively. An additional degree of freedom that couples the tangential force was provided by an inclined vibration absorber. Chatterjee et al.²⁴ investigated the usage of high-frequency tangential excitation in controlling self-excitation vibration. Dohnal²⁹ studied on a flow-induced self-excitation oscillation and presented stability analysis on vibration suppression employing variable stiffness and damping.

Recognizing the significance of the problem of control of self-excitation vibration in mechanical systems, this manuscript concerns a passive suppression method for self-excitation vibration. Self-excitation vibration suppression is achieved by internal stiffness and damping nonlinearities that are provided by appending horizontal springs and horizontal dampers, respectively.

Mechanical model

Figure 1 shows a lumped parameter model of flow-induced self-excitation system. The horizontal dampers and springs provide internal nonlinearities. The mass $m$ is a rigid body with length $l_{m}$ , width $w_{m}$ , and thickness $h_{m}$ . The vibration of the mass is confined to the x-direction. The horizontal stiffness $k_{h}$ and horizontal damping $c_{h}$ are placed onto the mass, besides vertical stiffness $k_{v}$ and vertical damping $c_{v}$ . The self-excitation vibration is induced by flow with constant velocity $u_{f}$ .

Figure 1.

Schematic representation of a lumped parameter model of flow-induced self-excitation system that has both stiffness and damping nonlinearities.

The internal force matrix $f$ can be written as

f = {\begin{matrix} f_{k} \\ f_{c} \end{matrix}} = {\begin{matrix} k_{v} x + 2 k_{h} (\frac{x}{{(x^{2} + l^{2})}^{1 / 2}}) ({(x^{2} + l^{2})}^{1 / 2} - l_{0}) \\ c_{v} x' + 2 c_{h} (\frac{x}{{(x^{2} + l^{2})}^{1 / 2}}) (\frac{x}{{(x^{2} + l^{2})}^{1 / 2}}) x' \end{matrix}}

(1)

where $f_{c}$ is the internal damping force, $f_{k}$ is the internal stiffness force, $l_{0}$ is the initial length of the lateral springs, and $l$ is their length when they are in the horizontal position.

Bernoulli principle is applied. The external force induced by the flow can be obtained¹

f_{d} = - 4 (l_{m} + w_{m}) ρ u_{f} h_{m} x'

(2)

where $ρ$ is the density of the flow.

The equivalent viscous damping, which is a commonly used concept in nonlinear systems, is exploited and this equivalence is based on the equal energy dissipation per cycle

c_{d} = - 4 (l_{m} + w_{m}) ρ u_{f} h_{m}

(3)

The equations of the motion for the system as shown in Figure 1 can be written as

\begin{matrix} mx ″ + (c_{v} + 2 c_{h} \frac{x^{2}}{x^{2} + l^{2}}) x' - c_{d} x' \\ + k_{v} x + 2 k_{h} (1 - \frac{l_{0}}{{(x^{2} + l^{2})}^{1 / 2}}) x = 0 \end{matrix}

(4)

The damping type can be approximated to the type of the damping referred to by the Van der Pol, provided that the level of vibration is not too large.

Limit circle

Calculation of limit circle

Equation (4) can be written in non-dimensional form as

\begin{matrix} \hat{x} ″ + (2 ζ_{v} + 4 ζ_{h} \frac{(1 - {\hat{l}}^{2}) {\hat{x}}^{2}}{(1 - {\hat{l}}^{2}) {\hat{x}}^{2} + {\hat{l}}^{2}} - 2 ζ_{d}) \hat{x}' \\ + (1 + 2 \hat{k} (1 - \frac{1}{{({\hat{x}}^{2} (1 - {\hat{l}}^{2}) + {\hat{l}}^{2})}^{1 / 2}})) \hat{x} = 0 \end{matrix}

(5)

where $ζ_{v} = c_{v} / 2 m ω_{n}$ , $ζ_{h} = c_{h} / 2 m ω_{n}$ , $ζ_{d} = c_{d} / 2 m ω_{n}$ , $ω_{n} = \sqrt{k_{v} / m}$ , $\hat{x} = x / x_{s}$ , $x_{s} = (l_{0}^{2} - l^{2})^{1 / 2}$ , $\hat{l} = l / l_{0}$ , $\hat{k} = k_{h} / k_{v}$ , and $(\cdot)' = d (\cdot) / d τ$ .

Harmonic balance method (HBM) is used. To improve the accuracy of the solution when relative low-order harmonic truncation is applied, the new time transformation is introduced³⁰

τ = σ + \sum_{n = 1}^{N} (\frac{α_{n}}{n} \cos n σ - \frac{β_{n}}{n} \sin n σ) + \sum_{n = 1}^{N} \frac{β_{n}}{n}

(6)

where $σ$ is a non-dimensional time transformation and $α_{n}$ and $β_{n}$ are constant transformation coefficients. Through this time transformation, the spectrum represented by the new time scale would be changed, so that the effect of higher order harmonic components could be reduced or even disappeared. For that more accuracy solution could be obtained when lower order truncation is used.

Derivation of equation (6) to $σ$ can be written as

\frac{d σ}{d τ} = g (σ) = \frac{1}{[1 + \sum_{n = 1}^{N} (α_{n} \cos n σ + β_{n} \sin n σ)]}

(7)

So, equation (5) has the form

F (\hat{x}, \hat{x}', \hat{x} ″, \hat{k}, ζ_{v}, ζ_{h}, ζ_{d}) = 0

(8)

where $(\cdot)' = d (\cdot) / d σ$ and $(\cdot) ″ = d^{2} (\cdot) / d σ^{2}$

Using HBM, the solution to equation (8) can be assumed as

\hat{x} = \sum_{n = 1}^{N} (a_{n} \cos n σ + b_{n} \sin n σ)

(9)

Truncation order N is the parameter that determines the accuracy of the solutions. Substituting equation (9) into equation (8), the coefficients of $\sin n τ$ and $\cos n τ$ equal 0; $2 N$ equations for $2 N + 1$ unknown parameters are obtained and can be solved using arc-length continuation.

Figure 2(a) shows the third, fifth, and seventh truncation HBM solutions. The parameters are $ζ_{h} = 0.25$ , $ζ_{v} = 0.25$ , $ζ_{d} = 0.75$ , $\hat{k} = 0$ , and $\hat{l} = 0.7$ . It can be seen that the difference between the limit circles is decreased as truncation order is increased. Seventh truncation solution is very similar to fifth truncation solution, and seventh truncation solution is reasonable and convincing. Convergence and divergence of the solutions depend on the truncation order. The results thus far have been as a result of using the HBM. When such an approximate method is used, it is necessary to check the results numerically, for example, by direct numerical integration of the equation of motion. To this end, the limit circle for the parameters used in Figure 2(b) with a horizontal length ratio $\hat{l} = 0.7$ is plotted for two cases (1) the analytical solution calculated using the seventh truncation HBM and (2) direct numerical integration of the matrix equation of motion at each frequency of excitation using Runge–Kutta method. The numerical solution is shown as “+”. It can be observed that there is a good agreement between the analytical and numerical results for the parameters chosen.

Figure 2.

The limit circle of the self-excitation system; the parameters are $ζ_{h} = 0.25$ , $ζ_{v} = 0.25$ , $ζ_{d} = 0.75$ , $\hat{k} = 1$ , and $\hat{l} = 0.7$ :(a) comparison of the limit circle between third, fifth, and seventh HBM solution. Red line: third HBM solution; blue line: fifth HBM solution; and black line: seventh HBM solution. (b) Comparison of the limit circle using seventh HBM and numerical method. Blue line: seventh HBM; black “+”: numerical method.

Table 1 shows comparison of the solutions between analytical and numerical methods. The points A, B, C, and D on the limit circle are chosen. It is found that the deviation of the harmonic balance solution is decreased as the truncation order is increased. The seventh truncation harmonic balance solution is very accurate to 0.05%–0.6%.

Table 1.

Comparison of the solutions between analytical and numerical methods.

	Numerical solution	HBM
		Third truncation	Fifth truncation	Seventh truncation
Point A	(1.7781, –0.6251)	(1.7683, –0.7986)	(1.7736, –0.5713)	(1.7805, –0.6251)
Error		2.9430%	1.1385%	0.1190%
Point B	(0.0728, –2.2147)	(0.0345, –2.3979)	(0.0751, –2.2955)	(0.0749, –2.2246)
Error		8.2246%	3.6476%	0.4495%
Point C	(–0.9677, –2.7640)	(–0.9754, –2.4502)	(–0.9728, –2.6796)	(–0.9692, –2.7468)
Error		9.9467%	2.6560%	0.5370%
Point D	(–1.9935, 0.0559)	(2.0280, –0.1010)	(–2.0040, 0.0916)	(–1.9940, 0.0767)
Error		1.8158%	0.5913%	0.0589%

HBM: harmonic balance method.

Characteristics of limit circle

The mechanical energy of the self-excitation system is the sum of kinetic energy, and elastic potential energy can be written as

E_{mech} = \frac{1}{2} \hat{x}'^{2} + \int_{0}^{τ_{p}} (1 + 2 \hat{k} (1 - \frac{1}{{({\hat{x}}^{2} (1 - {\hat{l}}^{2}) + {\hat{l}}^{2})}^{1 / 2}})) \hat{x} d \hat{x}

(10)

where $τ_{p}$ is the time interval through the limit circle.

The energy acting by the damping force can be written as

W_{d} = \int_{0}^{τ_{p}} (2 ζ_{v} + 4 ζ_{h} \frac{(1 - {\hat{l}}^{2}) {\hat{x}}^{2}}{(1 - {\hat{l}}^{2}) {\hat{x}}^{2} + {\hat{l}}^{2}} - 2 ζ_{d}) \hat{x}' d \hat{x}

(11)

The change in the mechanical energy of the self-excitation system $Δ E_{mech}$ depends on the energy acting by the damping force $W_{d}$ ; note that these quantities are the same for a self-excitation system $Δ E_{mech} = W_{d}$ .

Figure 3 shows the three-dimensional (3D) plots of the energy for the self-excitation vibration through the limit circle. Red line shows that energy is being supplied by the damper, and the mechanical energy of the system is increased. The blue line shows the energy is being dissipated. The reason for this is that the nonlinear damping performs significantly far from the equilibrium point but minimal near the equilibrium point. The sign of the energy only depends on the displacement. The energy is negative, until to the extreme level, and then, the energy is positive, as the displacement response is increased. Analysis of the power input and output characteristics through the limit circle could provide the theoretical basis for controlling the size of the limit circle, so that the self-excitation vibration could be suppressed.

Figure 3.

Normalized energy acting by damping for various non-dimensional velocity and non-dimensional displacement. The parameters are $ζ_{h} = 0.25$ , $ζ_{v} = 0.25$ , $ζ_{d} = 0.75$ , $\hat{k} = 1$ , and $\hat{l} = 0.7$ . Red line: the energy is being supplied; blue line: the energy is being dissipated; pink line: the projection of the limit circle; black line: the relationship between damping power and non-dimensional displacement through a limit circle; and green line: the relationship between damping power and non-dimensional velocity through a limit circle.

Vibration suppression

This section concerns the suppression of self-excitation vibration by regulating internal stiffness and damping nonlinearities. The steady amplitude of the stable limit circle is a key parameter that could be used to evaluate the suppressing effectiveness. It could be achieved by root mean square (RMS) of the time response

x_{max} = RMS (x (t))

(12)

Figure 4 illustrates the steady amplitude of the limit circle $x_{max}$ for variation flow velocity $u_{f}$ , showing the effects of changing the horizontal damping, vertical damping, horizontal stiffness, and vertical stiffness ( $c_{h}$ , $c_{v}$ , $k_{h}$ , $k_{v}$ ). Figure 4(a) shows the effects on the steady amplitude curve when horizontal damping $c_{h}$ is changed. For $u_{f} < 50 m / s$ , the system vibrates around the equilibrium position. For $u_{f} > 50 m / s$ , the system vibrates in the limit circle motion. The steady amplitude of the limit circle is increased as flow velocity $u_{f}$ is increased. The rate of roll on for the steady amplitude curve of the limit circle is decreased as horizontal damping $c_{h}$ is increased. Critical velocity $u_{critical}$ for the limit circle occurring is minimally affected. The reason for this is that as flow velocity $u_{f}$ is increased, the liner part of the damping is decreased and the displacement is increased, so the cubic nonlinear damping is enlarged. Figure 4(b) shows the effects on the steady amplitude of the limit circle when vertical damping $c_{v}$ is changed. The result shows that the critical velocity $u_{critical}$ is increased as the vertical damping $c_{v}$ is increased, but the rate of the roll on for the steady amplitude curve is minimally affected. Figure 4(c) shows the effects of changing the horizontal stiffness $k_{h}$ . It is found that the horizontal stiffness $k_{h}$ has minimal effect on the steady amplitude of the limit circle. Figure 4(d) shows the effects of changing the vertical stiffness $k_{v}$ . Simulations give the same conclusion as that of changing the horizontal stiffness. The limit circle is virtually unaffected by changing stiffness force.

Figure 4.

Illustration of the steady amplitude $x_{max}$ of the flow-induced self-excitation system with m = 1 kg, l = 0.35 m, l₀ = 0.5 m, l_m = 0.1 m, w_m = 0.1 m, h_m = 0.1 m, and ρ = 1.29 kg/m³, showing the effects of changing the horizontal damping, vertical damping, horizontal stiffness, and vertical stiffness ( $c_{h}$ , $c_{v}$ , $k_{h}$ , and $k_{v}$ ). (a) Effects of setting c_v = 20 N s/m, k_h = 10,000 N/m, and k_v = 1000 N/m and adjusting c_h: blue solid line, c_h = 20 N s/m; red dashed line, c_h = 40 N s/m; and black dotted line, c_h = 60 N s/m. (b) Effects of setting c_h = 20 N s/m, k_h = 10,000 N/m, and k_v = 1000 N/m and adjusting c_v: blue solid line, c_v = 20 N s/m; red dashed line, c_v = 40 N s/m; and black dotted line, c_v = 60 N s/m. (c) Effects of setting c_h = 20 N s/m, c_v = 20 N s/m, and k_v = 1000 N/m and adjusting $k_{h}$ : black solid line, k_h = 10,000 N/m and red “o,”k_h = 50,000 N/m. (d) Effects of setting c_h = 20 N s/m, c_v = 20 N s/m, and k_h = 10,000 N/m and adjusting $k_{v}$ : black solid points, $k_{v} = 1000 N / m$ , and blue “+” $k_{v} = 3000 N / m$ .

Conclusion

This article was concerned about the generation and suppression of flow-induced self-excitation vibrations in a type of spring-damper vertical mechanical system. Suppression of the self-excitation vibration was achieved by passive means. Internal stiffness and damping nonlinearities that coupled the vertical force were provided by appending horizontal springs and dampers. HBM was used for analysis of the limit circle of the system. The harmonic balance solutions with third, fifth, and seventh truncation were compared. The seventh truncation solution was also validated by numerical method. The energy transfer for the self-excitation vibration through the limit circle was analyzed. The amplitude of the limit circle was used to evaluate the effectiveness of the suppression of the self-excitation vibration, and the conclusions were as follows:

The rate of roll on for the steady amplitude curve of the limit circle was decreased as horizontal damping was increased, but critical velocity for the limit circle occurring was minimally affected.

The critical velocity was increased as the vertical damping was increased, but the rate of the roll on was unaffected.

Both the horizontal and vertical stiffness had minimal effect on the steady amplitude and critical velocity.

Footnotes

Acknowledgements

The authors of this article would like to acknowledge the help provided by Professor Michael J. Brennan from UNESP.

Handling Editor: Crinela Pislaru

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the State Key Program of National Natural Science of China (No.11232009) and National Natural Science Foundation of China (Nos 11502135 and 11572182).

References

Chen

. Analysis of crosswind fatigue of wind-excited structures with nonlinear aerodynamic damping. Eng Struct 2014; 74: 145–156.

Chatterjee

. Non-linear control of friction-induced self-excited vibration. Int J Nonlin Mech 2007; 42: 459–469.

Malas

Chatterjee

. Analysis and synthesis of modal and non-modal self-excited oscillations in a class of mechanical systems with nonlinear velocity feedback. J Sound Vib 2015; 334: 296–318.

Malas

Chatterjee

. Modal self-excitation by nonlinear acceleration feedback in a class of mechanical systems. J Sound Vib 2016; 376: 1–17.

Liu

et al . Selective photothermal self-excitation of mechanical modes of a micro-cantilever for force microscopy. Appl Phys Lett 2011; 99: 173501.

Yabuno

Kaneko

. Van der Pol type self-excited micro-cantilever probe of atomic force microscopy. Nonlinear Dynam 2008; 54: 137–149.

Lee

Lim

Moon

. A piezoelectric micro-cantilever bio-sensor using the mass-micro-balancing technique with self-excitation. Microsyst Technol 2007; 13: 563–567.

Kwasniewski

Dominik

Lalik

. Application of self-oscillating system for stress measurement in metal. J Vibroeng 2012; 14: 1392–8716.

Yabuno

Yasuhiro

Kuroda

. Self-excited coupled cantilevers for mass sensing in viscous measurement environments. Appl Phys Lett 2013; 103: 063104.

10.

Mizuno

Yaoita

Takeuchi

et al . Mass measurement using self-excited oscillations via two-relay controller. IEEE T Automat Contr 2009; 54: 416–420.

11.

Salem

SMS

. Study of wind turbine based self-excited induction generator under nonlinear resistive loads as a step to solve the Egypt electricity crisis. Comput Electr Eng 2016; 51: 1–11.

12.

Dominik

Iwaniec

Lech

. Low frequency damage analysis of electric Pylon model by fuzzy logic application. J Low Freq Noise V A 2013; 32: 239–251.

13.

Ding

. Steady-state responses of pulley-belt systems with a one-way clutch and belt bending stiffness. ASME J Vib Acoust 2014; 136: 041006.

14.

Ding

. Effect of one-way clutch on the nonlinear vibration of belt-drive systems with a continuous belt model. J Sound Vib 2013; 332: 6472–6478.

15.

Van der Pol . LXXXVIII. On relaxation oscillations. Lond Edinb Dublin Philos Mag J Sci 1926; 7: 978–992.

16.

Appleton

Van der Pol . On a type of oscillation-hysteresis in a simple triode generator. Lond Edinb Dublin Philos Mag J Sci 1922; 6: 177–193.

17.

Stoker

. Nonlinear vibrations in mechanical and electrical systems. New York: Interscience publishers, 1950.

18.

Raphael

. A radial solution of the van der Pol-Rayleigh limit circle oscillator for small ε. J Acoust Soc Am 1993; 94: 428–436.

19.

Andrianova

van Horssen

. Analytical approximations of the period of a generalized nonlinear van der Pol oscillator. J Sound Vib 2006; 295: 1099–1104.

20.

Waluya

van Horssen

. On the periodic solutions of a generalized nonlinear van der Pol oscillator. J Sound Vib 2003; 268: 209–215.

21.

Kovacic

. On the motion of a generalized van der Pol oscillator. Commun Nonlinear Sci 2011; 16: 1640–1649.

22.

Krack

Scheidt

Wallaschek

. A method for nonlinear model analysis and synthesis: application to harmonically forced and self-excited mechanical systems. J Sound Vib 2013; 332: 6798–6814.

23.

Kwaśniewki

Dominik

Lalik

. Influence of acoustoelastic coefficient on wave time of flight in stress measurement in piezoelectric self-excited system. Mech Syst Signal Pr 2016; 78: 143–155.

24.

Chatterjee

Singha

Karmakar

. Effect of high-frequency excitation on a class of mechanical systems with dynamic friction. J Sound Vib 2004; 269: 61–89.

25.

Malas

Chatterjee

. Generating self-excited oscillation in a class of mechanical systems by relay feedback. Nonlinear Dynam 2014; 76: 1253–1269.

26.

Chatterjee

. Self-excited oscillation under nonlinear feedback with time delay. J Sound Vib 2011; 330: 1860–1876.

27.

Das

Mallik

. Control of friction driven oscillation by time-delayed state feedback. J Sound Vib 2006; 297: 578–594.

28.

Chatterjee

Malas

. On the stiffness-switching methods for generating self-excited oscillations in simple mechanical systems. J Sound Vib 2012; 331: 1742–1758.

29.

Dohnal

. Suppressing self-excited vibrations by synchronous and time-periodic stiffness and damping variation. J Sound Vib 2007; 306: 136–152.

30.

Lau

Yuen

. The Hopf bifurcation and limit cycle by the incremental harmonic balance method. Comput Method Appl M 1991; 91: 1109–1121.