Multi-level residue harmonic balance method for nonlinear vibration of the beam

Abstract

This study presents the nonlinear vibration and chaotic response of a beam subjected to harmonic excitation. The multi-level residue harmonic balance method is applied to solve the geometrically cubic nonlinear vibration of the simply supported beam. The obtained results agree well with those of the numerical integration method. The amplitude frequency response curves are presented to illustrate the nonlinear dynamic system response both for a damping and without damping model. Also, the chaotic response is examined for a simply supported beam with a nonlinear dynamic system.

Introduction

Dynamic responses such as vibration behavior are an important part of structural analysis. In recent years, extensive research has been carried out, for example, on geometrically nonlinear beams, ¹ nonlinear vibration of a curved beam with quadratic and cubic nonlinearities, ² large deflection of a simply supported beam due to pure bending moment, ³ nonlinear dynamics of an axially moving viscoelastic beam, ⁴ nonlinear dynamics of a buckled beam subjected to primary resonance, ⁵ theoretical study of and experiments on a buckle beam, ⁶ and nonlinear vibrations and stability of an axially moving Timoshenko beam. ⁷ Nonlinear vibration and stability of duffing oscillators have been enormously investigated.^8–10 Lin et al. ¹¹ studied the nonlinear dynamic of a cantilever beam excited by periodic force. In their study, they considered the combined impacts and nonlinear terms in relation to beam deflection and employed numerical simulations in order to investigate chaotic motions. In the literature, as mentioned above, the nonlinear vibrations of beams are formulated by nonlinear partial differential equation in space and time with the different boundary conditions. The partial differential equations are discretized into ordinary differential equations as a mathematical model containing cubic or quadratic and cubic nonlinearity. Analytical solutions of such kind of nonlinear forced vibrations are highly complicated. Owing to the presence of nonlinearity in the mathematical models formulated from the governing equation, explicit solutions are rarely obtained.

A general computational formulation for geometrically nonlinear structures excited by harmonic forces and executing periodic motion in a steady state was developed by Lewandowski ¹² in part 1 and provided the resulting matrix amplitude equation and corresponding tangential matrix in an explicit form. Then, the numerical strategy for solving a resulting set of nonlinear algebraic equations and an example application of general theory for beam vibrations were presented by Lewandowski ¹³ in part 2. Although the efficiency and accuracy of the predicted method were demonstrated, the multi-level residue and higher-order analytical approximations were not considered in that analysis.

In the past, an extensive study with an approximate analytical method was carried out for solving nonlinear vibration problems, for example, harmonic balance^14–16 and incremental harmonic balance method (IHB). ¹⁷ The incremental harmonic balance method eliminates all nonlinear terms during the variational process. On the other hand, the total HB method keeps all the nonlinear terms to produce the possible multiple solutions in a set of nonlinear algebraic equations. In this method, the solution of nonlinear differential equation is assumed by Fourier series expansion to produce a set of nonlinear algebraic equations in terms of the Fourier coefficients. Since a set of nonlinear algebraic equations is produced, it is very difficult to use HB and IHB to construct higher-order analytical approximation. ¹⁸ They introduced an iterative homotopy HB method by constructing homotopy parameter P, which is considered a small parameter in the HB method. The analytical solutions of algebraic equations are also complicated. Traditional numerical methods are unable to provide all branches of solution of such nonlinear algebraic equation since the solution process is defined as the starting point as an initial guess.

In recent years, He and El-Dib ¹⁹ developed a new approach called the homotopy perturbation method (HPM) with the reducing rank method introduced to analyze a nonlinear Klein–Gordon equation with a strong damping parameter. They derived the oscillating solutions individually and finally obtained a frequency–amplitude formulation. In the same year, He and El-Dib ²⁰ developed an accurate frequency–amplitude relationship for describing the properties of the oldest device Fangzhu. He and Sedighi ²¹ applied the simplest frequency–amplitude formulation to study the Fangzhu device. They showed the effectiveness of HPM for solving a singular Duffing-like oscillator. In order to obtain frequency–amplitude formulation, an extensive study was carried out by Anjum et al. ²² They employed hybridization of the enhanced perturbation method and parameter expansion technology called the Li-He approach. This method is able to provide highly accurate results for nonlinear oscillators. For accurate analysis of nonlinear vibration behavior, Leung et al. ²³ developed the multi-parameter homotropy HB method for steady-state multiple solutions. The total and tangential stiffness matrix with respect to Fourier components of polynomial nonlinearity has been provided explicitly. The present study was carried out by adopting the multi-level residue HB method.

The main objective of this study was to perform an analytical solution for determination of the vibration behavior of a harmonically excited simply supported beam. The cubic nonlinear response of a beam to a primary-resonance excitation of its first vibration mode is reported. The multi-level residue HB method of solution process was used to solve the nonlinear dynamic response of a simply supported beam. The obtained results from this method were compared with those of nonlinear dynamic analysis using the direct numerical integration (NI) method; second, a frequency and amplitude relationship was investigated to illustrate the nonlinear dynamic system response both for a damping and without damping model, and finally, chaotic response was also examined.

Description of the problem

The governing equation of motion of nonlinear vibration of the Euler–Bernoulli theory including the effect of mid-plane stretching ²⁴ is as follows

$$EIW{,}_{xxxx}+ \rho AW{,}_{tt}+\xi W{,}_{t }+PW{,}_{xx}+KW-\left(\frac{EA}{2L}{\displaystyle {\int }_0^L{\left(W{,}_x\right)}^{2}dx}\right)W{,}_{xx}=F\hspace{0.17em}cos\left(\omega t\right)$$

(1)

where x is the longitudinal coordinate of the beam, the subscripts following commas stand for partial differentiation, W is the transverse displacement by the mid-plane stretching, and t denotes time.

Considering a straight beam on an elastic foundation with length L, the cross-sectional area of the beam A = B × h, the beam width B, the beam thickness h, the properties of the beam including the Young’s modulus of the beam E, the material density ρ = 2700 kgm⁻³, the damping coefficient of the beam ξ, and foundation modulus K, there will be no axial force acting on the beam; therefore, p = 0, F is the excitation amplitude (F = ρA κg), $\omega$ is the excitation frequency, κ is the base excitation magnitude, and g is the acceleration due to gravity = 9.81 ms⁻². The associated boundary conditions for a simply supported beam are

$$W\left(x,t\right)=W″\left(x,t\right)=0 , at x=0, L$$

(2)

The governing equation is discretized by applying the Galerkin procedure assuming that the transverse displacement is expressed in terms of simply supported beam mode shapes

$$W\left(x,t\right)= {\displaystyle {\sum }_i^n{q}_i\left(t\right){\phi }_i\left(x\right)}$$

(3)

where q _i = the modal amplitude of the ith mode, φ _i = the ith structural mode shape, ${\phi }_i\left(\text{x}\right)=\text{ sin}\left(n\pi x/L\right)$ , i = the structural mode number, and n = the number of structural modes considered.

Hence, the residual can be defined by substituting equation (3) into equation (1) as follows

$$\Delta =EI{\displaystyle \sum q_i{\phi }_i^{\left(iv\right)}}+ \rho A{\displaystyle {\sum }_{i=1}^n{\ddot{q}}_i{\phi }_i}+\text{ξ}{\displaystyle {\sum }_{i=1}^n{\dot{q}}_i{\phi }_i}- \frac{EA}{2L}{\displaystyle {\sum }_{i=1}^n{\displaystyle {\sum }_{j=1}^n{\displaystyle {\sum }_{k=1}^n{q}_i{q}_j{q}_k}}}{\displaystyle {\int }_0^L\left(\phi ″\prime {\displaystyle {\int }_0^L\phi {{}^{\prime }}_i\phi {{}^{\prime }}_{j }dx}\right)dx}-F\hspace{0.17em}\text{cos}\left(\omega t\right)$$

(4)

where $\phi {{}^{\prime }}_i,\phi {\mathit{"}}_i,\phi \mathit{"}{\mathit{"}}_i$ are the first, second, and fourth derivatives of the ith mode shape, respectively, and i and j are the mode numbers. Equation (4) contains cubic nonlinear term. According to the Galerkin approach, the weighted residual in equation (4) is set to zero. Multiplying ${\phi }_m$ by each term on the right-hand side of equation (4) and taking integration over the beam length, we have

$$\rho A{\displaystyle {\sum }_{i=1}^n{\ddot{q}}_i{\alpha }_{i,m}^{\mathrm{0,0}}}+EI{\displaystyle {\sum }_{i=1}^n{q}_i{\alpha }_{i,m}^{\mathrm{4,0}}}-\frac{EA}{2L}{\displaystyle {\sum }_{i=1}^n{\displaystyle {\sum }_{j=1}^n{\displaystyle {\sum }_{k=1}^n{q}_i{q}_j{q}_k{\alpha }_{i,m}^{\mathrm{2,0}}{\alpha }_{i,j}^{\mathrm{1,1}}}}}-F_m=0$$

(5)

where ${\alpha }_{i,m}^{\mathrm{0,0}}={\displaystyle {\int }_0^L{\phi }_i{\phi }_{m}dx}$ , ${\alpha }_{i,m}^{\mathrm{1,1}}={\displaystyle {\int }_0^L\phi {{}^{\prime }}_i\phi {{}^{\prime }}_{m}dx}$ , ${\alpha }_{i,m}^{\mathrm{2,0}}={\displaystyle {\int }_0^L\phi {{}^{″}}_i{\phi }_{m}dx}$ , ${\alpha }_{i,m}^{\mathrm{1,1}}={\displaystyle {\int }_0^L\phi {{}^{\prime }}_i\phi {{}^{\prime }}_{m}dx}$ , and ${\alpha }_{i,m}^{\mathrm{4,0}}={\displaystyle {\int }_0^L\phi ″{{}^{″}}_i{\phi }_{m}dx}$ .

For the single mode approach, equation (5) can be written as

$$\ddot{q}\left(t\right)+\xi \dot{q}\left(t\right)+\alpha q\left(t\right)+\beta q^3\left(t\right)= F\hspace{0.17em}\cos \left(\omega t\right)$$

(6)

where $\alpha =\left(EI{{\displaystyle \int }}_0^L\phi ″{}^{″}{\phi }_{m}dx\right)/\left(\rho A{{\displaystyle \int }}_0^L\phi {\phi }_{m}dx\right)$ and $\beta =\left(-EA{{\displaystyle \int }}_0^L\phi \prime \phi {{}^{\prime }}_{m}dx{{\displaystyle \int }}_0^L\phi ″\phi {{}^{\prime }}_{m}dx\right)/\left(2\rho AL{{\displaystyle \int }}_0^L\phi {\phi }_{m}dx\right)$ .

Residue harmonic balance method

Governing equation (6) can be written as

$$\psi \left(\ddot{q},\dot{q},q,\lambda ,t\right)\triangleq \ddot{q}+\xi \dot{q}+\alpha q+\beta q^3- F\hspace{0.17em}\cos \left(\omega t\right)=0$$

(7)

where $F\hspace{0.17em}\cos \left(\omega t\right)$ is an external periodic excitation term with period $T=2\pi /\omega$ . Since the periodic solutions are of interest, the steady-state response of equation (7) can be expanded in Fourier series

$$q_0\left(t\right)={\displaystyle {\sum }_{i=1}^n{a}_i\hspace{0.17em}\text{cos}\hspace{0.17em}\frac{i\omega t}{m}} +{\displaystyle {\int }_0^L{b}_i\hspace{0.17em}\text{sin}\hspace{0.17em}\frac{i\omega t}{m}}$$

(8)

where $a_i,b_i$ are the unknown Fourier coefficients, i is an integer representing super-harmonic responses, and $m$ is an integer corresponding to the ratio of the period of the external forcing to be period of the response representing sub-harmonic responses. It allows periodic solutions other than those of the fundamental period to be found. $n$ is also an integer representing the maximum harmonic order. It is retained in the Fourier series in equation (8) and assumed that $n\ge m$ . For convenience, a new time scale $\tau =\omega t/m$ of period $2\pi$ is introduced to eliminate $\omega$ ; then equation (7) becomes

$$\psi \left({\omega }^2,q″,q\prime ,q,\lambda ,\tau \right)\triangleq {\omega }^{2}q″+\xi \omega mq\prime +\alpha m^{2}q+\beta m^2{q}^3- F{m}^2\hspace{0.17em}\cos \left(m\tau \right)=0$$

(9)

The approximate harmonic balance solution of equation (8) is

$$q_0\left(\tau \right)={\displaystyle {\sum }_{i=1}^n{a}_i\hspace{0.17em}\text{cos}\left(i\tau \right)}+{\displaystyle {\sum }_{i=1}^n{b}_i\hspace{0.17em}\text{sin}\left(i\tau \right)}$$

(10)

Substituting equation (10) into equation (9) and applying the Galerkin procedure, we have the following harmonic balance equations

$$R_i^c\left(a_0,a_1,\dots ,a_n,b_1,\dots ,b_n\right)=\frac{2}{\pi }{\displaystyle {\int }_0^{\pi }\psi \left({\omega }^{2}q{{}^{″}}_0,\omega q{{}^{\prime }}_0,q_0,\lambda ,\tau \right)\cos \left(i\tau \right)d\tau }=0, \hspace{1em} i=0,\dots n,$$

$$R_i^s\left(a_0,a_1,\dots ,a_n,b_1,\dots ,b_n\right)=\frac{2}{\pi }{\displaystyle {\int }_0^{\pi }\psi \left({\omega }^{2}q{{}^{″}}_0,\omega q{{}^{\prime }}_0,q_0,\lambda ,\tau \right)\sin \left(i\tau \right)d\tau }=0, \hspace{1em} i=1,\dots n$$

(11)

Carrying out the integration yields the nonlinear algebraic equation

$$\left[K\left(L,\lambda \right)\right]\left\{L\right\}-\left\{F\right\}=\left\{R\right\}=0$$

(12)

where $\left\{R\right\}={\left[R_0^c,R_1^c,R_2^c,\dots R_n^c,R_1^s,R_2^s,\dots R_n^s\right]}^T$ is the vector of residues to be annihilated, $\left\{L\right\}={\left[2a_0,a_1,\dots ,a_n,b_1,\dots ,b_n\right]}^T$ and $\left\{F\right\}$ are the Fourier coefficient vectors of the response and excitation, respectively, and $\left[K\left(L,\lambda \right)\right]$ is the total stiffness matrix. In this study, the total stiffness matrices were programmed using the formula given by Fergusson and Leung ²⁵

$$\left[K\left(L,\lambda \right)\right]=\left({\mathrm{\Omega }}^2+\xi \mathrm{\Omega }+\alpha I+\beta {\mathit{H}}^2\right)\mathit{D}$$

where

$$\mathit{H}=\frac{1}{2}\left[\begin{array}{cc}{a}_0& p^T\\ 2p& Q+a_{0}I\end{array}\right], \hspace{1em}\mathit{D}=\left[\begin{array}{cc}1/2& 0\\ 0& I\end{array}\right],\hspace{1em}\mathrm{\Omega }=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& \mathrm{\Lambda }\\ 0& -\mathrm{\Lambda }& 0\end{array}\right]$$

$$Q=Q_0+Q_0^T+Q_1, Q_0=\left[\begin{array}{cc}{A}_0& B_0\\ -B_0& A_0\end{array}\right], Q_1=\left[\begin{array}{cc}{A}_1& B_1\\ -B_1& A_1\end{array}\right]$$

$$A_0=\left[\begin{array}{cccccc}0& a_1& a_2& a_3& a_4& \dots \\ 0& 0& a_1& a_2& a_3& \dots \\ 0& 0& 0& a_1& a_2& \dots \\ 0& 0& 0& 0& a_1& \dots \\ 0& 0& 0& 0& 0& \dots \\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right], A_1=\left[\begin{array}{cccccc}{a}_2& a_3& a_4& a_5& a_6& \dots \\ a_3& a_4& a_5& a_6& a_7& \dots \\ a_4& a_5& a_6& a_7& a_8& \dots \\ a_5& a_6& a_7& a_8& a_9& \dots \\ a_6& a_7& a_8& a_9& a_{10}& \dots \\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right]$$

$$B_0=\left[\begin{array}{cccccc}0& b_1& b_2& b_3& b_4& \dots \\ 0& 0& b_1& b_2& b_3& \dots \\ 0& 0& 0& b_1& b_2& \dots \\ 0& 0& 0& 0& b_1& \dots \\ 0& 0& 0& 0& 0& \dots \\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right], B_1=\left[\begin{array}{cccccc}{b}_2& b_3& b_4& b_5& b_6& \dots \\ b_3& b_4& b_5& b_6& b_7& \dots \\ b_4& b_5& b_6& b_7& b_8& \dots \\ b_5& b_6& b_7& b_8& b_9& \dots \\ b_6& b_7& b_8& b_9& b_{10}& \dots \\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right]$$

The unknown $u\left(t\right)=\left(1/2\right)a_0+{\mathit{p}}^T\mathit{z},\mathit{ }\mathit{z}=\left\{\cos \hspace{0.17em}i\omega t|\sin \hspace{0.17em}i\omega t\right\}, \cos \hspace{0.17em}i\omega t=\left\{\begin{array}{ccc}\text{cos}\hspace{0.17em}\omega t& \text{cos}\hspace{0.17em}2\omega t& \cdots \end{array}\right\}$ , etc. is presented by one of the following equivalent column vectors

$$\left\{\begin{array}{ccccccccc}{a}_0& a_1& a_2& a_3& \dots & b_1& b_2& b_3& \dots \end{array}\right\}=\left\{\begin{array}{ccc}{a}_0& {\mathit{a}}^T& {\mathit{b}}^T\end{array}\right\}=\left\{\begin{array}{cc}{a}_0& {\mathit{p}}^T\end{array}\right\}={\mathit{v}}^T$$

Then $u^n\left(t\right)=\left(1/2\right)a_0^{\left(n\right)}+{\mathit{p}}^{\left(n\right)T}\mathit{z}={\overline{\mathit{p}}}^{\left(n\right)T}\left\{1,\mathit{z}\right\}$ for positive integer n, where

$${\overline{\mathit{p}}}^{\left(n\right)}={\mathit{H}}^{\mathit{n}-1}\mathit{Dv}={\mathit{K}}^{\left(n\right)}\mathit{v}$$

Finally, the Jacobian associated with $u^n\left(t\right)$ is given by

$$\mathrm{J}=\mathrm{k}=\nabla \left(K^{\left(n\right)}\mathit{v}\right)=n{\mathit{K}}^{\left(n\right)}$$

The tangential stiffness is $\mathrm{k}=\left({\mathrm{\Omega }}^2+\xi \mathrm{\Omega }+\alpha I+3\beta {\mathit{H}}^2\right)\mathit{D}$ .

Hamiltonian-based frequency–amplitude formulation

We consider a nonlinear oscillator in the form

$$u″+f\left(u\right)=0, \hspace{1em}u\left(0\right)=A,\hspace{1em}u\prime \left(0\right)=0$$

(13)

The variational principle of equation (13) can be written as ²⁶

$$J\left(u\right)={\displaystyle \int \left(\frac{1}{2}u{{}^{\prime }}^2-F\left(u\right)du\right)}$$

(14)

where $\left(1/2\right)u{{}^{\prime }}^2$ is the kinetic energy, and $F\left(u\right)={\displaystyle \int f\left(u\right)du}$ is the potential energy.

The total energy remains unchanged during the oscillation

$$\frac{1}{2}u{{}^{\prime }}^2+F\left(u\right)=H$$

(15)

where $H$ is the Hamiltonian constant, which can be determined by the initial conditions given in equation (13). Finally, we obtain the following first-order differential equation

$$\frac{1}{2}u{{}^{\prime }}^2+F\left(u\right)-F\left(A\right)=0$$

(16)

We use the following trial function to determine the angular frequency

$$u\left(t\right)=A \hspace{0.17em}\cos \hspace{0.17em}\omega t$$

(17)

Substituting (17) into (16), we obtain the following residual equation

$$R\left(t\right)={\omega }^2{A}^2\hspace{0.17em}{\text{sin}}^2\left(\omega t\right)+F\left(A \cos \omega t\right)-F\left(A\right)$$

(18)

Similarly, we define two average residuals

$${\tilde{R}}_1=\frac{4}{T_1}{\displaystyle {\int }_0^{T_1/4}{R}_1\hspace{0.17em}\cos \hspace{0.17em}{\omega }_{1}tdt}$$

(19)

$${\tilde{R}}_2=\frac{4}{T_2}{\displaystyle {\int }_0^{T_2/4}{R}_2\hspace{0.17em}\cos \hspace{0.17em}{\omega }_{2}tdt}$$

(20)

A modification of the frequency–amplitude formulation is given as follows

$${\omega }^2=\frac{{\omega }_2^2 {\tilde{R}}_1-{\omega }_1^2 {\tilde{R}}_2}{{\tilde{R}}_1-{\tilde{R}}_2}$$

(21)

Homotopy continuation method

The analytical solutions of the nonlinear algebraic equations are rather difficult to obtain accurate and all kinds of possible solutions. The initial starting values for iteration in the numerical method cannot be defined properly; therefore, this method is unable to provide all solutions. The homotopy continuation method is a global numerical method which is capable of solving nonlinear equations with all solutions and useful for a wide variety of problems. Hence, this method is applied to peak up all solutions of equation (12). We defined a target equation for solution in a linear-product homotopy and homogenization. Homogenization is recommended because path endpoints approaching infinity are very common and the available projective transformation after homogenization keeps both the magnitudes of the coordinates and the arc-lengths of the homotopy paths finite. ²⁷ For the nonlinear algebraic equation (12)

$$\left[K\left(L,\lambda \right)\right]\left\{L\right\}-m^{2}F\left\{e_{m+1}\right\}=0$$

(22)

where $\left\{e_{m+1}\right\}$ denotes a unit column vector having unity at the ${\left(m+1\right)}^{th}$ position and zeros elsewhere. In a unit-homogeneous treatment on a continuation parameter $\gamma$ using the coordinates $\left[a_0,a_1, \mathrm{...},a_n,b_1,b_2, \mathrm{...} , b_n,\gamma \right]\in P^{2n+1}$ for projection on the 2n+1 dimension plane, equation (22) becomes

$$\left[\left({\omega }^2{\mathrm{\Omega }}^2\gamma +\xi \omega \mathrm{\Omega }\gamma +\alpha I{m}^2\gamma +m^2\beta {\mathit{H}}^2\right)\mathit{D}\right]\left\{L\right\}-m^{2}F\left\{e_{m+1}\right\}{\gamma }^2=0$$

(23)

and we specify a compatible homogeneous structure to append an inhomogeneous linear equation ${\alpha }_0{a}_0+{\alpha }_1{a}_1+\cdots +{\alpha }_n{a}_n+{\beta }_1{b}_1+{\beta }_n{b}_n+\beta \gamma =1$ for some random, complex vector $\{{\alpha }_0, {\alpha }_1,\dots ,{\alpha }_n,{\beta }_n,\gamma \}$ . By choosing a random path on $P^{2n+1}$ to the material equation, we specify the linear-product structure. The steady-state harmonic solutions may be obtained by solving nonlinear algebraic equation (22) using the above method. However, this is rather difficult for large $n$ . Therefore, we often select $n$ as a small integer that satisfies $n\ge m$ in solving the steady-state responses for reasonable accuracy. In this study, we selected $m=1, n=1$ in equation (6) for studying the steady-state super-harmonic response.

Residue

Substituting the obtained harmonic solutions (8) into governing equation (7), the residue vector ${\tilde{R}}_0$ can be expressed as follows

$$\Vert {\tilde{R}}_0\Vert =\text{sup}\left\{\left|{\displaystyle {\sum }_{i=n+1}^{2n}{R}_i^c\cos \left(i\tau \right)+R_i^c\sin \left(i\tau \right)}\right|, \tau \in \left[\mathrm{0,2}\pi \right]\right\}\ne 0$$

(24)

where

$$R_i^c\left(a_0,a_1,\dots ,a_n,b_1,\dots ,b_n\right)=\frac{2}{\pi }{\displaystyle {\int }_0^{\pi }\psi \left({\omega }^{2}q{{}^{″}}_0,\omega q{{}^{\prime }}_0,q_0,\lambda ,\tau \right)\cos \left(i\tau \right)d\tau }=0,\hspace{1em}i=n+1,\dots 2n$$

$$R_i^s\left(a_0,a_1,\dots ,a_n,b_1,\dots ,b_n\right)=\frac{2}{\pi }{\displaystyle {\int }_0^{\pi }\psi \left({\omega }^{2}q{{}^{″}}_0,\omega q{{}^{\prime }}_0,q_0,\lambda ,\tau \right)\sin \left(i\tau \right)d\tau }=0,\hspace{1em}i=1+1,\dots 2n$$

Note that the residue vector $\left\{R\right\}$ in equation (12) is associated with $i=\mathrm{0,1},\dots ,n$ while $\left\{{\tilde{R}}_0\right\}$ is associated with $i=n+1,\dots , 2n$ . The residues ${\tilde{R}}_0$ are convergent to zero when the representation (10) is sufficiently accurate. If they are far from zero, the representation (10) is not physically realistic. Therefore, the residuals ${\tilde{R}}_0$ are the measures of the accuracy about the steady-state harmonic solutions.

Assume that the accurate harmonic solutions of equation (9) are in the form

$$q\left(\tau \right)=q_0\left(\tau \right)+\epsilon q_1\left(\tau \right)+{\epsilon }^2{q}_2\left(\tau \right)+\cdots$$

(25)

where $\epsilon$ is an embedding parameter and $q_j\left(\tau \right)={\displaystyle {\sum }_{i=0}^{\left(j+1\right)n}{a}_{i,k}\hspace{0.17em}\text{cos}\left(i\tau \right)}+{\displaystyle {\sum }_{i=1}^{\left(j+1\right)n}{b}_{i,k}\hspace{0.17em}\text{sin}\left(i\tau \right)}$ $j=1, 2\dots$ . When $q_j\left(\tau \right)$ is included, the number of terms is virtually j times more than what was represented in equation (10). We shall find $q_j\left(\tau \right)$ successively in the following sub-sections. To find $q_0\left(\tau \right)$ , one needs to solve the 2n + 1 nonlinear algebraic equation (12). However, when $q_0\left(\tau \right)$ qualitatively describes the solution $q\left(\tau \right)$ sufficiently, the solution can be improved quantitatively by including $q_j\left(\tau \right)$ which requires solving $2\left(j+1\right)n+1$ linear equations sequentially in the following manner.

Zero level solution ${\mathit{\epsilon }}^0$

Substituting equation (25) into equation (9) and equating the coefficients of ${\epsilon }^0$ we have

$$\begin{array}{l}\psi \left({\omega }^2,q{{}^{″}}_0,q{{}^{\prime }}_0,q_0,\lambda ,\tau \right)\triangleq {\omega }^{2}q{{}^{″}}_0+\xi \omega mq{{}^{\prime }}_0+\alpha m^2{q}_0+\beta m^2{q}_0^3- F{m}^2\cos \left(m\tau \right)\\ ={\displaystyle {\sum }_{i=0}^{2n}{R}_i^c\cos \left(i\tau \right)}+{\displaystyle {\sum }_{i=1}^{2n}{R}_i^s\sin \left(i\tau \right)}\end{array}$$

Then, after separating the summation into two parts

$$\begin{array}{l}\psi \left({\omega }^2,q{{}^{″}}_0,q{{}^{\prime }}_0,q_0,\lambda ,\tau \right)\\ ={\displaystyle {\sum }_{i=0}^n{R}_i^c\cos \left(i\tau \right)}+{\displaystyle {\sum }_{i=1}^n{R}_i^s\sin \left(i\tau \right)}+ \epsilon \left({\displaystyle {\sum }_{i=n+1}^{2n}{R}_i^c\cos \left(i\tau \right)}+{\displaystyle {\sum }_{i=n+1}^{2n}{R}_i^s\sin \left(i\tau \right)}\right)=\epsilon {\tilde{R}}_0\end{array}$$

(26)

The zero solution ${\epsilon }^0$ can be obtained in equation (10) for the initial $2n+1$ unknowns $a_0,a_1,\dots ,a_n,b_1,\dots ,b_n$ by solving nonlinear algebraic equation (12). The following sub-sections are for the ${\epsilon }^j$ solution by solving $2\left(j+1\right)n+1$ linear equations.

First-level solution ${\mathit{\epsilon }}^1$

Substituting equations (25) and (26) into equation (9) and equating the coefficients of ${\epsilon }^1$ terms, one has the following linear ordinary differential equation in $q_1\left(\tau \right)$

$${\omega }^{2}q{{}^{″}}_1+\xi \omega mq{{}^{\prime }}_1+\alpha m^2{q}_1+3\beta m^2{q}_0^2{q}_1+{\tilde{R}}_0={\displaystyle {\sum }_{i=0}^{3n}{R}_{i,1}^c\cos \left(i\tau \right)}+{\displaystyle {\sum }_{i=1}^{3n}{R}_{i,1}^s\sin \left(i\tau \right)}$$

Then, after separating the summation into two parts

$${\omega }^{2}q{{}^{\prime }}_1+\xi \omega mq{{}^{\prime }}_1+\alpha m^2{q}_1+3\beta m^2{q}_0^2{q}_1+{\tilde{R}}_0={\displaystyle {\sum }_{i=0}^{2n}{R}_{i,1}^c\cos \left(i\tau \right)}+{\displaystyle {\sum }_{i=1}^{2n}{R}_{i,1}^s\sin \left(i\tau \right)}+\epsilon {\tilde{R}}_1$$

(27)

Based on harmonic balance, we eliminate the terms $\text{cos}\left(i\tau \right)\left(i=\mathrm{0,1},\dots ,2n,\sin \left(i\tau \right)\right)\left(i=1,\dots ,2n\right)$ by equating their coefficients to zero

$$R_{i,1}^c\left(a_{\mathrm{0,1}},a_{\mathrm{1,1}},\dots ,a_{2n,1},b_{\mathrm{1,1}},\dots ,b_{2n,1}\right)=0,\hspace{1em}i=\mathrm{0,1}, \dots ,2n,$$

$$R_{i,1}^s\left(a_{\mathrm{0,1}},a_{\mathrm{1,1}},\dots ,a_{2n,1},b_{\mathrm{1,1}},\dots ,b_{2n,1}\right)=0, \hspace{1em}i=1, \dots ,2n,$$

(28)

There are $4n+1$ independent unknowns $a_{\mathrm{0,1}},a_{\mathrm{1,1}},\dots ,a_{2n,1},b_{\mathrm{1,1}},\dots ,b_{2n,1}$ to be solved by linear equation (28) which contains $4n+1$ equations. Once they are determined, according to equations (10), (25), and (28), the first-level solution can be obtained as

$$q\left(\tau \right)={\displaystyle {\sum }_{i=1}^n\left(a_i+a_{i,1}\right)\cos \left(i\tau \right)}+{\displaystyle {\sum }_{i=1}^n\left(b_i+b_{i,1}\right)\sin \left(i\tau \right)}+{\displaystyle {\sum }_{i=n+1}^{2n}\left(a_{i,1}\cos \left(i\tau \right)+b_{i,1}\text{sin}\left(i\tau \right)\right)}$$

(29)

Note that the zero-level ${\epsilon }^0$ solution is obtained by solving nonlinear algebraic equation (12) while the first-level ${\epsilon }^1$ solution is determined by solving linear equation (28).

Second-level solution ${\mathit{\epsilon }}^2$

Substituting equations (25) and (27) into equation (9) and equating the coefficients of the ${\epsilon }^2$ terms yields the following linear ordinary differential equation in $q_2\left(\tau \right)$

$$\begin{array}{l}{\omega }^{2}q{{}^{″}}_2+\xi \omega mq{{}^{\prime }}_2+\alpha m^2{q}_2+3\beta m^2{q}_0^2{q}_2+3\beta m^2{q}_0{q}_1^2+{\tilde{R}}_1\\ ={\displaystyle {\sum }_{i=0}^{4n}{R}_{i,2}^c\cos \left(i\tau \right)}+{\displaystyle {\sum }_{i=1}^{4n}{R}_{i,2}^s\sin \left(i\tau \right)}\end{array}$$

Then, after separating the summation into two parts

$$\begin{array}{l}{\omega }^{2}q{{}^{″}}_2+\xi \omega mq{{}^{\prime }}_2+\alpha m^2{q}_2+3\beta m^2{q}_0^2{q}_2+3\beta m^2{q}_0{q}_1^2+{\tilde{R}}_1\\ ={\displaystyle {\sum }_{i=0}^{3n}{R}_{i,2}^c\cos \left(i\tau \right)}+{\displaystyle {\sum }_{i=1}^{3n}{R}_{i,2}^s\sin \left(i\tau \right)}+\epsilon {\tilde{R}}_2\end{array}$$

(30)

Based on the harmonic balance procedure, we eliminate the terms $\text{cos}\left(i\tau \right)\left(i=\mathrm{0,1},\dots ,3n,\sin \left(i\tau \right)\right)\left(i=1,\dots ,3n\right)$ by letting their coefficients to be zero

$$R_{i,2}^c\left(a_{\mathrm{0,2}},a_{\mathrm{1,2}},\dots ,a_{3n,2},b_{\mathrm{1,2}},\dots ,b_{3n,2}\right)=0, \hspace{1em}i=\mathrm{0,1}, \dots ,3n,$$

$$R_{i,2}^s\left(a_{\mathrm{0,2}},a_{\mathrm{1,2}},\dots ,a_{3n,2},b_{\mathrm{1,2}},\dots ,b_{3n,2}\right)=0,\hspace{1em}i=1, \dots ,3n,$$

(31)

There are $6n+1$ independent unknowns $a_{\mathrm{0,2}},a_{\mathrm{1,2}},\dots ,a_{3n,2},b_{\mathrm{1,2}},\dots ,b_{3n,2}$ to be determined by the $6n+1$ linear equation (31). Thus, the second-level solution can be obtained from the following equation

$$q\left(\tau \right)={\displaystyle {\sum }_{i=0}^n\left(a_i+a_{i,1}+a_{i,2}\right)\cos \left(i\tau \right)}+{\displaystyle {\sum }_{i=1}^n\left(b_i+b_{i,1}+b_{i,2}\right)\sin \left(i\tau \right)}+{\displaystyle {\sum }_{i=n+1}^{2n}\left[\left(a_{i,1}+a_{i,2}\right)\cos \left(i\tau \right)+\left(b_{i,1}+b_{i,2}\right)\text{sin}\left(i\tau \right)\right]}+{\displaystyle {\sum }_{i=2n+1}^{3n}\left(a_{i,2}\cos \left(i\tau \right)+b_{i,2}\sin \left(i\tau \right)\right)}$$

(32)

Modified homotopy perturbation method

Homotopy perturbation method is an impotent method for solving nonlinear vibration problems. Several authors modified this method to improve the accuracy of the method. Li and He proposed a modified homotopy perturbation technique called Li-He method combining homotopy perturbation and enhanced perturbation methods. In this section, a modified homotopy perturbation method (Li-He method) to solve the nonlinear forced vibration problem is presented as follows ²⁰ .

Consider the equation

$$\ddot{q}+\xi \dot{q}+\alpha q+\beta q^3= F\hspace{0.17em}\cos \left(\omega t\right)$$

(33)

Operating on equation (30) by $\left(D^2+{\omega }^2\right)$ , we obtain

$$\left(D^2+{\omega }^2\right)\left(D^2+\xi D+\alpha +\beta q^2\right)q= 0$$

(34)

Equation (34) can be written as follows

$$\left(D^2+{\omega }^2\right)\left(D^2+\alpha \right)q=-\left(D^2+{\omega }^2\right)\left(\xi \dot{q}+\beta q^3\right)=-\left(\xi \dot{q}+\xi {\omega }^2\dot{q}+3\beta q^2\ddot{q}+6\beta q{\dot{q}}^2+\beta {\omega }^2{q}^3\right)$$

(35)

Substituting $\ddot{q}$ from equation (33) into equation (35) and constructing the corresponding homotopy equation has the following form

$$\begin{array}{l}\left(D^2+{\omega }^2\right)q=\rho \left[-\frac{1}{D^2+\alpha }\left(\xi \dot{q}+\xi {\omega }^2\dot{q}+3\beta q^2\left(-\xi \dot{q}-\alpha q-\beta q^3+ F\cos \left(\omega t\right)\right)+6\beta q{\dot{q}}^2+\beta {\omega }^2{q}^3\right)\right]\\ =\rho \left[-\frac{1}{D^2+\alpha }\left(\xi \dot{q}+\xi {\omega }^2\dot{q}-3\beta \xi \dot{q}{q}^2+6\beta q{\dot{q}}^2+\left({\omega }^2-3\alpha \right)\beta q^3-3{\beta }^2{q}^5+ 3\beta F{q}^2\cos \left(\omega t\right)\right)\right]\end{array}$$

(36)

The solution of equation (36) is chosen as

$$q=q_0+\rho q_1+{\rho }^2{q}_2+\cdots$$

(37)

Substituting equation (37) into equation (36) and equating the identical power of $\rho$ to zero, we have the system of linear algebraic equation as follows

$${\ddot{q}}_0+{\omega }^2{q}_0=0$$

(38)

$${\ddot{q}}_1+{\omega }^2{q}_1=-\frac{1}{D^2+\alpha }\left(\xi {\stackrel{⃛}{q}}_0+\xi {\omega }^2{\dot{q}}_0-3\beta \xi {\dot{q}}_0{q}_0^2+6\beta q_0{\dot{q}}_0^2+\left({\omega }^2-3\alpha \right)\beta q_0^3-3{\beta }^2{q}_0^5+3\beta F{q}_0^2\cos \left(\omega t\right)\right)$$

(39)

………………………………………………………………………………………………….

The solution of equation (38) is as follows

$$q_0=a\hspace{0.17em}\cos \left(\omega t\right)+b\hspace{0.17em}\sin \left(\omega t\right)$$

(40)

Substituting equation (40) into equation (39) and eliminating secular terms, we obtain

$$\frac{9}{4}{a}^{2}F\beta +\frac{3}{4}{b}^{2}F\beta -\frac{9}{4}{a}^3\alpha \beta -\frac{9}{4}a{b}^2\alpha \beta -\frac{15a^5{\beta }^2}{8}-\frac{15}{4}{a}^3{b}^2{\beta }^2-\frac{15}{8}a{b}^4{\beta }^2-\frac{3}{4}{a}^{2}b\beta \xi \omega -\frac{3}{4}{b}^3\beta \xi \omega +\frac{9}{4}{a}^3\beta {\omega }^2+\frac{9}{4}a{b}^2\beta {\omega }^2=0$$

(41)

$$\frac{3}{2}abF\beta -\frac{9}{4}{a}^{2}b\alpha \beta -\frac{9}{4}{b}^3\alpha \beta -\frac{15}{8}{a}^{4}b{\beta }^2-\frac{15}{4}{a}^2{b}^3{\beta }^2-\frac{15b^5{\beta }^2}{8}+\frac{3}{4}{a}^3\beta \xi \omega +\frac{3}{4}a{b}^2\beta \xi \omega +\frac{9}{4}{a}^{2}b\beta {\omega }^2+\frac{9}{4}{b}^3\beta {\omega }^2=0$$

(42)

Solving equations (41) and (42), the values of $a$ and $b$ are calculated for different values of $\omega$ .

Results and discussion

The comparison of displacement amplitude between zero- and second-level solutions both for damping and without damping of the cubic nonlinear beam is shown in Figure 1. This figure shows that zero-level and second-level solutions are much closer to each other. A little difference of displacement amplitude corresponding to normalized frequency of 1.5 is about 0.55% and the maximum difference is about 1.2% corresponding to normalized frequency of 2.053.

Figure 1.

Comparison of displacement amplitude between zero- and second-level (a) without damping (b) with damping of the nonlinear cubic beam.

On the other hand, Figure 2(b) represents a difference of displacement amplitude occurring at normalized frequency $\left(\omega /{\omega }_0\right)=1.095$ of about 1.97% and the maximum difference is about 4.16% corresponding to normalized frequency of 2.053. Therefore, it is important for an accurate analysis of the cubic nonlinear beam with damping to continue upto the higher-level solution rather than zero-level solution. However, the maximum difference between the zero- and second-level solutions is observed within 5%. A representative result corresponding to normalized frequency of 0.864 is summarized in Table 1.

Figure 2.

(a) Comparison of time history and (b) deformation shape plot of the existing stable responses when the normalized frequency = 0.684.

Table 1.

Zero- and second-level solutions when $\left(\omega /{\omega }_0\right)=0.684$ .

	$q_0\left(t\right)$	q2(t)
Constant	$0$	$0$
$\text{cos}\left(\omega t\right)$	$1.16834\times {10}^{-5}$	$1.16745\times {10}^{-5}$
$\text{sin}\left(\omega t\right)$	$2.05427\times {10}^{-7}$	$2.04985\times {10}^{-7}$
$\text{cos}\left(3\omega t\right)$	—	$1.17051\times {10}^{-7}$
$\text{sin}\left(3\omega t\right)$	—	$2.28473\times {10}^{-8}$

Comparison of results

The results obtained from the RHB method are $q\left(\tau \right)=1.16745\times {10}^{-5}\cos \left(\tau \right)+2.04985\times {10}^{-7}\sin \left(\tau \right)+1.17051\times {10}^{-7}\cos \left(3\tau \right)+2.28473\times {10}^{-8}\sin \left(3\tau \right)$ corresponding to the normalized driving frequency of 0.684. A numerical integration method (NI) was also employed to compare the results obtained from the RHB method. Note that numerical simulation was performed assuming that the initial conditions of both the displacement and velocity are zero. In Figures 2(a) and 3(b), the solid and dashed lines represent the solution obtained from the RHB and NI methods for time histories and phase trajectories, respectively. The solution from two methods agrees well both for time histories and phase trajectories. A representative deformation shape diagram corresponding to the periodic response obtained by the RHB method is shown to illustrate the deformation of the beam under the external harmonic excitation in Figure 2(b). It shows that the maximum deformation is at the center. Figure 3(a) shows the Poincare section of the nonlinear continuous dynamic system. The single point attributes in this Poincare section that means the solution of the nonlinear dynamic system is periodic with period-1. One can say in other words there is no period doubling occurred in that system. Figure 4, with representative phase trajectories at the different values of driving frequency of 30, 40, and 50, shows that velocity increases with constant displacement.

Figure 3.

(a) Poincare map of the existing stable responses when ω = 0.684 and (b) comparison of phase trajectories.

Figure 4.

Phase trajectories varying with different frequencies.

Backbone curves

In Figure 5, a representative amplitude response curve is shown to illustrate the jump phenomenon. The backbone curve is plotted considering no damping. There is exactly one solution branch for $\left(120/2\pi \right)<\left(370/2\pi \right)$ and three coexisting solutions $\left(370/2\pi \right)<\left(370/2\pi \right)$ to $\propto$ . Therefore, the steady-state solutions are $\left(120/2\pi \right)<\left(370/2\pi \right)$ and the lower displacement amplitude regarding the frequency $\left(370/2\pi \right)<\left(370/2\pi \right)$ to $\propto$ .

Figure 5.

Backbone curve of the beam equation (12) without damping.

In Figure 6, a comparison between amplitude–frequency response curves is shown to illustrate considering the damping model when the base excitation $\kappa =0.1$ . The overall displacement amplitudes are calculated by using the formula $\left|\sqrt{a_i^2+b_i^2}\right|$ . From the figure, it is seen that the result obtained by the proposed method agrees reasonably well with the corresponding results obtained by the classical harmonic balance method. However, the results obtained by the homotopy perturbation method slightly deviated. An exact one solution branch was found for $\left(\omega /{\omega }_0\right)<1.095$ and three co-existing solutions for $1.095<\left(\omega /{\omega }_0\right)<2.053$ , which is termed “interval of bistability”. At $\left(\omega /{\omega }_0\right)=1.095$ , only two solutions were found; thereafter, at this location two solutions merged to form more solutions with the increase of driving frequency. Therefore, the periodic response of the cubic beam with external harmonic excitation loses stability leading to a jump in the response. To elucidate the importance of several coexisting solutions in Figure 6, we classified the three branches of solutions as (i) upper, (ii) lower, and (iii) middle branches of solutions. Now the driving frequency which is a control parameter gradually increased from $\left(\omega /{\omega }_0\right)>1.095$ , and the response amplitude will provide upper branch or the large amplitude solution branch of the response diagram. When $\left(\omega /{\omega }_0\right)=2.053$ , the large amplitude forced vibration decreases and a transition takes place to lower branches consisting of small amplitude solutions. Thus, a jump forms from the large amplitude solution branch to the small amplitude solution branch.

Figure 6.

Comparison between the amplitude versus frequency response curve obtained by the proposed method, harmonic balance method, and He’s homotopy perturbation method.

If the driving frequency is slowly decreased from the normalized frequency of 2.053, the displacement amplitude of the steady-state forced response increases. At $\left(\omega /{\omega }_0\right)=1.095$ , a transition occurs to a solution on the upper branch and lower branch, and vice versa at different values of the driving frequency. Thus, hysteresis arises. The middle branch solutions are unstable since they are perturbed and will not return to the original solution, so normally they merge to another solution.

Chaotic response

It was found that the transition occurs when the normalized frequencies are 1.095 and 2.053. The chaotic responses are summarized in Figure 7, through the time histories, phase trajectories, and the Poincare map.

Figure 7.

(a) Time history, (b) phase trajectories, and (c) Poincare map of irregular attractor of equation (12) when the normalized frequency = 2.053.

Conclusion

In this study, the RHB method was applied to examine nonlinear vibration of a straight Euler–Bernoulli’s beam with cubic nonlinearities subjected to harmonic excitation. The Galerkin method was employed to discretize the governing partial differential equation to a single-mode nonlinear ordinary differential equation with mid-plane stretching effect. This method can be employed for the analysis of multi-modes vibration of beam and plate. Also, this method has enabled to solve nonlinear couple equations, and hence it can be used for sub-harmonic analysis of beam and plate. From this study it can be concluded that: (i) no period doubling was observed, (ii) the comparison shows good agreement between the present results and the numerical results along with the time histories, phase trajectories, and backbone curves, (iii) the quasi-periodic response was observed in the present study, and the small magnitude of external harmonic excitation is the main reason for a quasi-periodic response, and (iv) chaotic responses through time histories, phase trajectories, and Poincare sections were investigated numerically. The backbone curves prove that the chaotic region started after the primary resonance at the normalized frequency of 1.095.