Abstract
In this paper, the dynamic stability analysis of a simply supported beam excited by a sequence of moving masses is investigated. All components of the mass acceleration including the centripetal, the Coriolis and the vertical one are considered. The periodical traverse of masses across the beam results to a linear timeperiodic problem. The Floquet theory and the Incremental Harmonic Balance (IHB) method are implemented to obtain the boundary between stable and unstable regions in the parameters plane. A new approach for identifying the conditions of resonance is investigated by presenting an intuitive definition of resonance for timevarying systems. This approach enables the IHB method to determine inherent curves of resonance conditions besides its ability to find the boundary curve separating the stable and unstable regions. Numerical simulations confirm the correctness of resulted curves.
1. Introduction
A substantial variety of practical systems in engineering can be represented as a flexible beam carrying a moving mass. Many applications such as motion of vehicles or trains on bridges, cranes transporting loads along its span, robotic arms, fluid transfer pipe systems, space structures and high speed machining processes are some of these systems, to enumerate a few. Therefore, determining dynamic behavior of a beam subjected to moving mass has been an interesting subject of investigation for a long time. Historically, the first known attempts to solve a moving load problem arose in the study of the collapse of Chester Railway Bridge by Willis [1] and Stokes [2]. After that, there have been a lot of efforts in this field, including the investigations done by Ayre and Jacobson [3] and two well known monographs by Inglis [4] and Hillerborg [5]. Recently, two books have been published by Fryba [6] and Yang et al. [7]. The growing usage of heavy and rapid truck vehicles besides the employment of lighter and more flexible structures have drawn engineers attention on terms which were not considered in the past [817].
Studies about the dynamic behavior of transiting objects across flexible beams can be classified into two major groups. The first group studies dynamic behavior of a beam under moving mass transition in time and/or frequency domain, focusing mainly on simulations and numerical methods. The second group mainly focuses on the analysis of system stability. This latter involves identification of system parameters for which system’s instability occurs or finding the resonance frequencies in presence of external excitation. These analyses usually conclude to the determination of stability boundary curves by semianalytical or numerical methods. Despite numerous published papers for finding dynamic response of such systems, the number of investigations on the conditions for dynamic stability is restricted. Nelson and Conover [18] utilized the Floquet theory to determine the parametric regions of stability of a simplysupported EulerBernoulli beam lying on a uniform elastic foundation excited by a continuous series of equallyspaced mass particles. Benedetti [19] pointed out that in Nelson and Conover’s study, multiple separate regions of instability may appear and just for certain combinations of particle intervals and foundation moduli, single instability region occurs. Katz et al. [20] studied the dynamic stability and transverse vibration of a simply supported beam subjected to a deflection dependent moving load caused by the cutting forces in machining operations. Aldraihem and Baz^{}[21] investigated dynamic stability of a stepped beam subjected to a moving mass. They used the impulsive parametric excitation theory to find the stable range for beam’s parameters when subjected to periodic parametric excitations. Verichev and Metrikine [22] studied the stability of vibrations induced by a moving mass along a beam lying on a periodically inhomogeneous continuous foundation using perturbation analysis. Mackertich [23] studied the dynamic stability of an elastically supported Timoshenko beam under travelling masses. He utilized the Floquet theory to determine regions of stability in the parameters plane.
It should be noted that all previous related studies on dynamic stability of the beammoving mass problem just introduced the boundary curve that separates stable and unstable regions in the parameters plane. Although some references [6, 2426] studied the problem of repetitive moving forces crossing the beam and identified those velocities of the force that conclude in resonance of the beam, to the best of the author's knowledge, the resonance analysis of a beam excited by moving masses has not been addressed before in open literature.
In this study, the mathematical model for a simply supported beam excited by a continuous sequence of equallyspaced identical masses is presented. The Galerkin method is then used to discretize the spatial domain and to obtain the governing ordinary differential equation of motion. Next, the IHB method is used to analyze the resulted equation and an iterative approach is implemented to determine the boundary of instability and the resonance curves in the plane of moving mass parameters. Numerical simulations are employed to confirm the correctness of the resulted curves.
2. Mathematical modeling and derivation of equation of motion
Fig. 1Schematic of a flexible beammoving mass problem
The system under consideration, which is composed of a traversing mass $m$ and a simply supported uniform EulerBernoulli beam of length $l$, mass density $\rho $, crosssectional area $A$ and flexural rigidity $EI$, is schematically shown in Fig. 1. The purpose of investigating a flexible system influenced by repetitive mass transition is to deal with a timevarying system on which stability and parametric resonance studies can be performed comprehensively. Due to timevarying identity, such systems can exhibit unexpected phenomena that may occur in practice.
It is assumed that the mass moves along the beam from left to right with a constant velocity $V$ and remains always in contact with the beam. The governing equation that describes small vertical vibration of the beam represented by $v(x,t)$, can be written as [25]:
where $g$ is the gravitational acceleration and $\delta $ is the Dirac delta function.
In order to apply Galerkin method to discretize Eq. (1), the solution is considered as being expressible by the following expansion:
where ${q}_{i}\left(t\right)$ is the corresponding influence coefficient for the $i$th shape function ${\phi}_{i}\left(x\right)$ of the undamped beam modes (without the moving mass), defined by:
Substituting Eq. (2) into Eq. (1), then multiplying the resultant equation by the $j$th shape function ${\phi}_{j}\left(x\right)$, next integrating over the beam length and considering the orthogonality condition between modes, the partial differential Eq. (1) is converted to a set of ordinary differential equations on the modal coordinates described by:
where $\mathbf{q}=\left[{q}_{1}\right(t),\mathrm{}{q}_{2}(t),...,{q}_{n}(t){]}^{T}$ is the array of the modal coordinates. The matrix components are expressed as:
where ${\delta}_{ij}$ and ${\omega}_{i}$ are the Kronecker delta and the $i$th natural frequency of a simply supported beam $\left({\omega}_{i}={\left(i\pi /l\right)}^{2}\sqrt{EI/\rho A}\right)$, respectively. Considering one mode, the equation governing the modal coordinate becomes:
To rewrite the governing equation in dimensionless form, nondimensional parameters are defined as:
Using the chain rule and the nondimensional parameters given in Eq. (7), the dimensionless equation results as:
where dot denotes derivation with respect to $\tau $.
Eq. (8) consists of a linear equation whose coefficients vary with time up to the moment that the moving mass is still on the beam span. As soon as it leaves the beam span, the mass will have no more influence on the beam and consequently the timevarying coefficients of the governing Eq. (8) will vanish ($\alpha =\text{0}$), resulting in the removal of any vibration amplification factor. By considering the following sequence that once a mass leaves the beam, the next one enters the beam span and consequently results in a repetitive excitation, coefficients of Eq. (8) become periodical with period ${T}_{p}=l/V$. Performing Fourier expansion of the coefficients of Eq. (8) leads to:
To study the dynamic stability of the system, the right hand side of Eq. (9) is momentarily disregarded, which follows as:
3. Incremental harmonic balance method
IHB method represents an ingenious technique for establishing dynamic stability, applicable to the most general form of differential systems. According to IHB method, $\text{(}{\alpha}^{\mathrm{*}}\text{,}\text{}{\mathrm{\Omega}}^{\mathrm{*}}\text{)}$ is considered as a point on the merge of instability corresponding to a periodic solution ${Q}^{\mathrm{*}}\left(\tau \right)$ of Eq. (10). A neighboring point on this boundary is denoted by:
Substituting Eq. (11) into Eq. (10) and maintaining linear terms of $\mathrm{\Delta}Q\left(\tau \right)$, $\mathrm{\Delta}\alpha $ and $\mathrm{\Delta}\mathrm{\Omega}$, the linear incremental equation is obtained as:
where:
$\left.\mathrm{}\mathrm{}\mathrm{}\mathrm{}+\left(1{\alpha}^{\mathrm{*}}{\mathrm{\Omega}}^{{\mathrm{*}}^{2}}\left(1\mathrm{c}\mathrm{o}\mathrm{s}\left(2\tau \right)\right)\right){Q}^{\mathrm{*}}\right]$
is the corrective term and will go to zero on instability boundary points [27]. The next step of the IHB method is to find an approximate solution by applying Galerkin's method. The unknown periodic solution ${Q}^{\mathrm{*}}$ can be expanded into a truncated Fourier series as:
Accordingly, $\mathrm{\Delta}Q$ is expressed as:
For conciseness, the functions ${Q}^{\mathrm{*}}$ and $\mathrm{\Delta}Q$ can be written as:
where:
$\mathbf{A}={\left[{a}_{0},\mathrm{}{a}_{1},\mathrm{}{a}_{2},\dots ,\mathrm{}{b}_{1},\mathrm{}{b}_{2},\dots \right]}^{T},$
$\mathrm{\Delta}\mathbf{A}={\left[\mathrm{\Delta}{a}_{0},\mathrm{}\mathrm{\Delta}{a}_{{\mathrm{}}_{1}},\mathrm{}\mathrm{\Delta}{a}_{{\mathrm{}}_{2}},\mathrm{}\dots \mathrm{},\mathrm{}\mathrm{\Delta}{b}_{{\mathrm{}}_{1}},\mathrm{}\mathrm{\Delta}{b}_{{\mathrm{}}_{2}},\mathrm{}\dots \right]}^{T}.$
Applying the Galerkin procedure to Eq. (12) yields:
$={\int}_{0}^{{T}_{p}}\delta {\left(\mathrm{\Delta}Q\right)}^{T}\left[R\left(\begin{array}{c}{\mathrm{\Omega}}^{{\mathrm{*}}^{2}}\left(1\mathrm{cos}\left(2\tau \right)\right){\ddot{Q}}^{\mathrm{*}}+2{\mathrm{\Omega}}^{{\mathrm{*}}^{2}}\mathrm{sin}\left(2\tau \right){\dot{Q}}^{\mathrm{*}}\\ {\mathrm{\Omega}}^{{\mathrm{*}}^{2}}\left(1\mathrm{cos}\left(2\tau \right)\right){Q}^{\mathrm{*}}\end{array}\right)\mathrm{\Delta}\alpha \right.$
$\left.\left(\begin{array}{c}2{\mathrm{\Omega}}^{\mathrm{*}}\left(1+{\alpha}^{\mathrm{*}}\left(1\mathrm{cos}\left(2\tau \right)\right)\right){\ddot{Q}}^{\mathrm{*}}+4{\mathrm{\Omega}}^{\mathrm{*}}{\alpha}^{\mathrm{*}}sin\left(2\tau \right){\dot{Q}}^{\mathrm{*}}\\ 2{\mathrm{\Omega}}^{\mathrm{*}}{\alpha}^{\mathrm{*}}\left(1\mathrm{c}\mathrm{o}\mathrm{s}\left(2\tau \right)\right){Q}^{\mathrm{*}}\end{array}\right)\mathrm{\Delta}\mathrm{\Omega}\right]d\tau $
Substituting the terse form of Eq. (16) into Eq. (18), a set of linear equations containing $\mathrm{\Delta}\mathbf{A}$, $\mathrm{\Delta}\alpha $ and $\mathrm{\Delta}\mathrm{\Omega}$ as main variables are obtained as:
where:
Eq. (19) is utilized to obtain the instability boundary. It is a set of nonlinear equations in terms of ${\alpha}^{\mathrm{*}}$, ${\mathrm{\Omega}}^{\mathrm{*}}$ and ${\mathbf{A}}^{\mathrm{*}}$. Instead of solving these nonlinear equations for ${\alpha}^{\mathrm{*}}$, ${\mathrm{\Omega}}^{\mathrm{*}}$ and ${\mathbf{A}}^{\mathrm{*}}$, a set of equations which are linear in terms of $\mathrm{\Delta}\mathbf{A}$, $\mathrm{\Delta}\alpha $ and $\mathrm{\Delta}\mathrm{\Omega}$ are solved in a recursive approach.
In Eq. (19) the number of unknown variables exceeds the number of equations by two. Yet, since in vector ${\mathbf{A}}^{\mathrm{*}}$ only the relative values of the coefficients are needed, it is possible to prescribe one of them as a unity reference constant with its corresponding increment in $\mathrm{\Delta}\mathbf{A}$ equal to zero. Also by selecting one parameter out of $\alpha $ or $\mathrm{\Omega}$ as the active parameter and consequently setting its increment equal to zero (i.e., $\mathrm{\Delta}\alpha =\text{0}$ or $\mathrm{\Delta}\mathrm{\Omega}=\text{0}$), two unknowns are already set as zero and the equations can be solved for the rest.
As explained, let us select $\alpha $ as the active parameter and set $\mathrm{\Delta}\alpha $ and the $i$th element of $\mathrm{\Delta}\mathbf{A}$ equal to zero. Then ${\mathrm{\Omega}}^{\mathrm{*}}$ and ${\mathbf{A}}^{\mathrm{*}}$ are found through the following recursive algorithm:
1) Select ${\alpha}^{\mathrm{*}}$.
2) Select initial values for ${\mathrm{\Omega}}^{\mathrm{*}}$ and ${\mathbf{A}}^{\mathrm{*}}$.
3) Solve the following linear equation:
where ${\stackrel{~}{\mathbf{S}}}_{\mathrm{\Delta}A}$ is the modified form of ${\mathbf{S}}_{\mathrm{\Delta}A}$ by eliminating its $i$th column and $\mathrm{\Delta}\stackrel{~}{\mathbf{A}}$ is the changed form of vector $\mathrm{\Delta}\mathbf{A}$ by removing its $i$th element.
4) Check for convergence of the solution by examining $\Vert \mathbf{R}\Vert \le \epsilon $, where $\epsilon $ is an enough small number.
5) Update ${\mathrm{\Omega}}^{\mathrm{*}}$ and ${\mathbf{A}}^{\mathrm{*}}$ as ${\mathrm{\Omega}}_{New}^{\mathrm{*}}={\mathrm{\Omega}}^{\mathrm{*}}+\mathrm{\Delta}\mathrm{\Omega}$ and ${A}_{New}^{\mathrm{*}}={A}^{\mathrm{*}}+\mathrm{\Delta}A$ if needed and repeat the procedure from step 2 if convergence is not reached.
6) Select a new value for ${\alpha}^{\mathrm{*}}$ and repeat steps 1 to 5 to determine the corresponding ${\mathrm{\Omega}}^{\mathrm{*}}$ and ${\mathbf{A}}^{\mathrm{*}}$.
By considering series expansion of different order for the solution, convergence of the instability boundary curve and other resonance curves can be expected with arbitrary accuracy. All these curves are determined through the abovementioned iterative procedure.
4. Results and discussion
4.1. Dynamic stability
Results of dynamic stability investigation for a simply supported EulerBernoulli beam excited by a continuous sequence of identical equally spaced moving masses are presented in this section. According to Floquet theory for $T$periodic systems, the solutions corresponding to the transition curves that separate stable and unstable regions in the parameters plane are of either $T$ or 2$T$ period. Therefore, the following series expansions are considered for ${Q}^{\mathrm{*}}\left(\tau \right)$ and $\mathrm{\Delta}Q\left(\tau \right)$ to represent the $T$periodic solution:
and also the following expansions for the 2$T$periodic solution:
The trend consists of finding those parameters of traversing mass in $(\mathrm{\alpha}\mathrm{\Omega})$ plane which result to a periodic response for Eq. (10), hence establishing the boundary curve iteratively. This curve, which separates stable and unstable regions in the system parameters plane, is depicted in Fig. 2. As shown, the regions above and below this curve display unstable and stable behavior, respectively. In order to verify the validation of the IHB method, this curve is compared to the one resulted by Floquet theory. Coincidence of these two curves proves the accuracy of the analysis performed by presented procedure.
Fig. 2Stable and unstable regions from IHB method (five term expansion) and Floquet theory
Fig. 3Vibration simulation for parameters selected in the vicinity of the a) stable and b) unstable sides
a)
b)
Fig. 3 depict the numerical simulations of the beam midspan deflection (i.e. $v\left(0.5l,t\right)={\phi}_{1}\left(0.5l\right){q}_{1}\left(t\right)$) for two sets of numerical values of $(\mathrm{\alpha}\mathrm{\Omega})$ from stable and unstable regions. Beam parameters are selected as $EI=$5×10^{5} N.m^{2}, $\rho A=\text{10}\text{}\text{kg/m}$ and $l=$20 m. The operational set points $\text{(}\mathrm{\alpha}=\text{0.3,}\mathrm{}\mathrm{\Omega}=\text{1.76)}$ and $\text{(}\mathrm{\alpha}=\text{0.3,}\mathrm{}\mathrm{\Omega}=\text{1.86)}$ belong to the stable and unstable regions, respectively. The corresponding stable and unstable behavior of the system can be seen in Fig. 3. It should be noticed that although there is an apparent damping term in the equation, nonetheless the solution in the stable region does not diminish to zero but remains bounded as predicted by Floquet theory.
4.2. Occurrence of resonance in stability region
In this section, the conditions under which the moving mass transition conducts to external resonance in the system are investigated. It is revealed that there are some curves lying in the stable region leading to situations where the beam experiences vibrations of growing amplitudes.
The superposition principle for linear timevarying systems permits us to consider separately each term on the righthand side of Eq. (27) as a plausible exciting term that may induce resonance in the system. It is claimed that if a solution can be synchronized with one of those excitation terms, similar to timeinvariant systems, resonance can occur. Indeed, the definition of natural frequency is extended here for a timevarying system as the frequency of its homogeneous response, like for timeinvariant systems. By this way, those parameters of the moving mass are being sought which generate a solution with principal frequency corresponding to one of the excitation frequencies. In order to fulfill this condition, let's consider a periodic solution with the same principal period as the kth excitation term:
As discussed, the IHB method permits to design a response which matches one of the excitation terms and consequently realizes the conditions for resonance. The expected resonance curve can be obtained by following the procedure explained in section 3. Through increasing $k$, different candidate solutions corresponding to other excitation terms can be generated, resulting to multiple resonance curves shown in Fig. 4. More concisely, those values of massspeed parameters in the $(\mathrm{\alpha}\mathrm{\Omega})$ plane lying on these curves will cause resonance in the system.
Fig. 4Resonance curves located in the stable region
Although these resonance curves have not been addressed in previous studies, their intersection with vertical axis, corresponding to moving load transition $\text{(}\mathrm{\alpha}=\text{0)}$, have been reported in [25] as:
as it is verified in Fig. 4.
In order to further validate the obtained results, numerical simulations of the beam midspan deflection (i.e. $v\left(0.5l,t\right)={\phi}_{1}\left(0.5l\right){q}_{1}\left(t\right)$) are provided using four sets of numerical values for $\mathrm{\alpha}$ and $\mathrm{\Omega}$ located on these resonance curves, namely $\text{(}\mathrm{\alpha}=\text{0.3,}\mathrm{\Omega}=\text{0.4302)}$, $\text{(}\mathrm{\alpha}=\text{0.3,}\text{}\mathrm{\Omega}=\text{0.2198)}$, $\text{(}\mathrm{\alpha}=\text{0.3,}\mathrm{}\mathrm{\Omega}=\text{0.1480})$ and $\text{(}\mathrm{\alpha}=\text{0.3,}\mathrm{}\mathrm{\Omega}=\text{0.1115)}$. Fig. 5 depicts the simulation results for beam midspan response. As seen from the figures, selected parameters for the moving mass cause resonant behavior in the beam. It is worth noting that such intermittent resonance conditions appearing for a determined mass ratio were not expected in practice.
Fig. 5Beam midpoint response for parameters selected from a) 1st, b) 2nd, c) 3rd, d) 4th resonance curve
a)
b)
c)
d)
5. Conclusions
In the present study, the problem of dynamic stability of a simply supported EulerBernoulli beam under the periodic loading of sequential moving masses is investigated by adopting the IHB method. The resulted instability curve, which separates the stable and unstable regions, is verified by comparing with the one derived by applying Floquet theory. In addition, by employing above method, curves describing resonance conditions in the stable parameter region are depicted. These curves, in fact, represent those moving mass parametric values that yield resonance in the beam response. Numerical simulations verify the occurrence of resonance for the parameters selected on these curves.
References

Willis R. Appendix to the Report of the Commissioners Appointed to Inquire into the Application of Iron to Railway Structures. H. M. Stationary Office, London, 1849.

Stokes G. G. Discussion of a differential equation relating to the breaking of railway bridges. Transactions of the Cambridge Philosophical Society, Vol. 8, Issue 5, 1849, p. 707735.

Ayre R. S., Jacobson L. S. Transverse vibration of a twospan beam under the action of a moving alternating force. Journal of Applied Mechanics, Vol. 17, Issue 3, 1950, p. 283290.

Inglis C. E. A Mathematical Treatise on Vibrations in Railway Bridges. Cambridge University Press, London, 1934.

Hillerborg A. Dynamic influences of smoothly running loads on simply supported girders. Inst. of Structural Engineering and Bridge Building of the Royal Inst. of Technology, Stockholm, 1951.

Fryba L. Vibration of Solids and Structures Under Moving Loads. Thomas Telford Ltd., London, 1999.

Yang Y. B., Yau J. D., Wu Y. S. Vehicle Bridge Interaction Dynamics: with Applications To High Speed Railways. World Scientific Publishing, Singapore, 2004.

Gbadeyan J. A., Oni S. T. Dynamic behaviour of beams and rectangular plates under moving loads. Journal of Sound and Vibration, Vol. 182, Issue 5, 1995, p. 677695.

Michaltsos G. T., Kounadis A. N. The effects of centripetal and coriolis forces on the dynamic response of light bridges under moving loads. Journal of Vibration and Control, Vol. 7, Issue 3, 2001, p. 315326.

Michaltsos G. T., Sophianopoulos D., Kounadis A. N. The effect of moving mass and other parameters on the dynamic response of a simply supported beam. Journal of Sound and Vibration, Vol. 191, Issue 3, 1996, p. 357362.

Rao V. G. Linear dynamics of an elastic beam under moving loads. Journal of Vibration and Acoustics, Vol. 122, Issue 3, 2000, p. 281289.

Mamandi A., Kargarnovin M. H., Farsi S. An investigation on effects of traveling mass with variable velocity on nonlinear dynamic response of an inclined Timoshenko beam with different boundary conditions. International Journal of Mechanical Sciences, Vol. 52, Issue 12, 2010, p. 16941708.

Nikkhoo A., Rofooei F. R., Shadnam M. R. Dynamic behavior and modal control of beams under moving mass. Journal of Sound and Vibration, Vol. 306, Issues 35, 2007, p. 712724.

Nayyeri Amiri S., Onyango M. Simply supported beam response on elastic foundation carrying repeated rolling concentrated loads. Journal of Engineering Science and Technology, Vol. 5, Issue 1, 2010, p. 5266.

Kargarnovin M. H., Ahmadian M. T., JafariTalookolaei R. A. Dynamics of a delaminated Timoshenko beam subjected to a moving oscillatory mass. Mechanics Based Design of Structures and Machines, Vol. 40, Issue 2, 2012, p. 218240.

Ariaei A., ZiaeiRad S., Ghayour M. Vibration analysis of beams with open and breathing cracks subjected to moving masses. Journal of Sound and Vibration, Vol. 326, Issues 35, 2009, p. 709724.

Zou J. H., Feng W. X., Jiang H. B. Dynamic response of an embedded railway track subjected to a moving load. Journl of Vibroengineering, Vol. 13, Issue 3, 2011, p. 544551.

Nelson H. D., Conover R. A. Dynamic stability of a beam carrying moving masses. Journal of Applied Mechanics, Vol. 38, Issue 4, 1971, p. 10031006.

Benedetti G. A. Dynamic stability of a beam loaded by a sequence of moving mass particles. Journal of Applied Mechanics, Vol. 41, Issue 4, 1974, p. 10691071.

Katz R., Lee, W. U., Ulsoy A. G., Scott R. A. Dynamic stability and response of a beam subjected to a deflection dependent moving load. Journal of Vibration and Acoustics, Vol. 109, Issue 4, 1987, p. 361365.

Aldraihem O. J., Baz A. Dynamic stability of stepped beams under moving loads. Journal of Sound and Vibration, Vol. 250, Issue 5, 2002, p. 835848.

Verichev S. N., Metrikine A. V. Instability of vibrations of mass that moves uniformly along a beam on a periodically inhomogeneous foundation. Journal of Sound and Vibration, Vol. 260, Issue 5, 2003, p. 901925.

Mackertich S. Dynamic stability of a beam excited by a sequence of moving mass particles. Acoustical Society of America, Vol. 115, Issue 4, 2004, p. 14161419.

Dahlberg T. Vehiclebridge interaction. Vehicle System Dynamics, Vol. 13, Issue 4, 1984, p. 187206.

Ouyang H. Movingload dynamic problems: a tutorial (with a brief overview). Mechanical Systems and Signal Processing, Vol. 25, Issue 6, 2011, p. 20392060.

Yau J. D. Vibration of simply supported compound beams to moving loads. Journal of Marine Science and Technology, Vol. 12, Issue 4, 2004, p. 319328.

Lau S. L., Yuen S. W. Solution diagram of nonlinear dynamic systems by the IHB method. Journal of Sound and Vibration, Vol. 167, Issue 2, 1993, p. 303316.