A modified harmonic balance method for solving forced vibration problems with strong nonlinearity

Abstract

In this paper, a modified harmonic balance method is presented to solve nonlinear forced vibration problems. A set of nonlinear algebraic equations appears among the unknown coefficients of harmonic terms and the frequency of the forcing term. Usually a numerical method is used to solve them. In this article, a set of linear algebraic equations is solved together with a nonlinear one. The solution obtained by the proposed method has been compared to those obtained by variational and numerical methods. The results show good agreement with the results obtained by both methods mentioned above.

Introduction

Nonlinear vibration is an important issue in science and engineering. Most of the differential equations involving physical phenomena are nonlinear. Therefore, the analysis of nonlinear vibration problems has attracted intensive research attention in the past few decades.^1–6 Generally, it is very hard to find the exact or closed-form solution of such a nonlinear problem. Therefore, many researchers have paid attention to both numerical and analytical methods. Numerical methods are comparatively easy but require heavy computational effort and a proper initial gauge value to obtain the desired results. Moreover, numerical methods cannot provide an overall view of the behavior of the nonlinear dynamical systems. Consequently, the quest for accurate behavior of the nonlinear dynamical systems led to the development of many analytical approximations. In the literature, several analytical approximate methods are found such as perturbation,^7,8 homotopy analysis,^9,10 homotopy perturbation,^11,12 variational iteration,^13,14 harmonic balance (HBM),^15–18 etc.

The perturbation methods^19–23 are a widely used technique for handling nonlinear problems which were originally developed for weak nonlinear problems. Furthermore, the classical perturbation methods have been extended or modified by several authors to investigate strong nonlinear problems. Jones ²⁴ introduced an approximate method to enhance the range of validity of the classical perturbation where the parameter is not small. Later, Cheung et al. ²⁵ modified the Lindstedt–Poincare method using the concept of Jones. ²⁴ Recently, Alam, Yeasmin and Ahamed ²⁶ have generalized the modified Lindstedt–Poincare method ²⁵ which covers a wide variety of nonlinear oscillators.

The HBM is another powerful technique for determining periodic solutions of nonlinear differential equations where a truncated Fourier series is chosen as the solution of the differential equations. In the classical HBM, a set of nonlinear algebraic equations is solved by a numerical method to determine the unknown coefficients. Furthermore, the classical HBM has been modified by some researchers. For example, Rahman et al. ²⁷ used the HBM to investigate the nonlinear Ven der Pol equation. Lau and Cheung ²⁸ introduced an incremental variational principle for nonlinear vibration of an elastic system. Azrar et al. ²⁹ presented a semi-analytical approach based on the concept of harmonic balance to analyze the nonlinear response of a large-amplitude beam. Rahman and Lee ³⁰ introduced a modified multi-level residue HBM for solving nonlinearly vibrating double-beam problems. Belendez et al. ³¹ applied a modified HBM to solve a class of strongly nonlinear oscillators including a rational nonlinear term. Wagner and Lentz ³² used a HBM to investigate a Duffing oscillator with excitation force. They also studied an extended Duffing oscillator.

Recently, variational and Hamiltonian principles are two interesting approaches to solve strong nonlinear problems. Earlier, He ³³ used the semi-inverse method to find the variational principle for handling some nonlinear partial differential equations with variable coefficients. Furthermore, He ³⁴ applied this method to solve some nonlinear oscillators with fractional nonlinear terms. He ³⁵ developed a new variational method for nonlinear oscillators using the Hamiltonian approach. Akbarzade and Kargar ³⁶ also applied the Hamiltonian approach to nonlinear free and forced vibrating systems without any damping effect. Then Yildirim et al. ³⁷ utilized the same method to solve nonlinear oscillators with rational and irrational elastic terms. Sadeghzadeh and Kabiri ³⁸ presented a higher-order Hamiltonian approach to the nonlinear vibration of micro electro-mechanical systems without any damping effect. From the above literature, it is seen that in most cases, variational and Hamiltonian principles are applied for nonlinear free vibration problems. Yet for the higher-order approximate solution, a set of nonlinear algebraic equations is raised and generally solved by the numerical method.

In this paper, a modified harmonic balance is presented. The advantage of the proposed method is that a set of linear algebraic equations together with a nonlinear one are solved, thereby reducing the computational effort and requiring less computational time than other HBMs. The method is applied to solve nonlinear forced vibration problems. The results have been compared to those obtained by numerical and a variational method to verify the accuracy of the method.

The method

Consider a nonlinear forced vibration equation of the following form

$$\ddot{x}+{\omega }_0^{2}x+\mu \dot{x}+\epsilon f\left(x\right)=\mathit{E}\cos (\omega t)$$

(1)where over dots denote differentiation with respect to $t$ , ${\omega }_0$ is the natural frequency, $\mu$ is the damping coefficient, $f\left(x\right)$ is a given nonlinear function of $x$ , $\epsilon$ is a parameter, $E$ is the excitation force amplitude, and $\omega$ is the excitation frequency.

The solution of equation (1) is chosen as follows

$$x=\text{a}\cos \left(\omega t\right)+\mathit{b}\sin \left(\omega t\right)+a_3\cos \left(3\omega t\right)+b_3\sin \left(3\omega t\right)+\cdots$$

(2)

Substituting equation (2) into equation (1) and expanding $f\left(x\right)$ in Fourier series and equating the coefficients of equal harmonics, the following equations are obtained

$$a\left(-{\omega }^2+{\omega }_0^2\right)+b\mu \omega +\epsilon C_1\left(a,b,a_3,b_3,\cdots \right)=E$$

(3a)

$$b\left(-{\omega }^2+{\omega }_0^2\right)-a\mu \omega +\epsilon S_1\left(a,b,a_3,b_3,\cdots \right)=0$$

(3b)

$$a_3\left(-9{\omega }^2+{\omega }_0^2\right)+3b_3\mu \omega +\epsilon C_3\left(a,b,a_3,b_3,\cdots \right)=0$$

(3c)

$$b_3\left(-9{\omega }^2+{\omega }_0^2\right)-3a_3\mu \omega +\epsilon S_3\left(a,b,a_3,b_3,\cdots \right)=0$$

(3d)

Eliminating ${\omega }^2$ from equations (3b–d) with the help of equation (3a), and neglecting the terms whose responses are negligible, the above equations can be re-written as follows

$${\omega }^2={\omega }_0^2+\epsilon C_1\left(a,b,a_3,b_3,\cdots \right)-E/a$$

(4a)

$$-a\mu \omega +\epsilon C_1\left(a,b,a_3,b_3,\cdots \right)+\epsilon S_1\left(a,b,a_3,b_3,\cdots \right)-bE/a=0$$

(4b)

$$-8{\omega }_0^2{a}_3+a_3\epsilon C_1\left(a,b,a_3,b_3,\cdots \right)+\epsilon C_3\left(a,b,a_3,b_3,\cdots \right)-a_{3}E/a=0$$

(4c)

$$-8{\omega }_0^2{b}_3+b_3\epsilon C_1\left(a,b,a_3,b_3,\cdots \right)+\epsilon S_3\left(a,b,a_3,b_3,\cdots \right)-3a_3\mu \omega -b_{3}E/a=0$$

(4d)

Now using equation (4b), eliminating $\omega$ from equations (4c and d) and considering only linear terms of $a_3$ , $b_3$ and neglecting the terms whose responses are negligible, a system of linear algebraic equations of $a_3$ , $b_3$ is obtained. Solving these two equations, $a_3$ , $b_3$ are determined in terms of $a,b$ . Then the values of $a_3$ , $b_3$ are substituted into equation (4b) and $b$ is expressed in powers of small parameter $\lambda (\mu , \omega ,E)$ as follows

$$b=l_0+l_1\lambda +l_2{\lambda }^2+l_3{\lambda }^3+\cdots$$

(5)where l₀,l₁,l₂,… are functions of $a$ . Finally, substituting $a_3$ , $b_3$ , and $b$ into equation (4a) and solving, the value of $a$ is obtained. Consequently, $b$ , $a_3$ , and $b_3$ are obtained.

Solution obtained by proposed method

Example 1

Consider a damped forced vibration equation of the form

$$\ddot{x}+x+\mu \dot{x}+\epsilon x^3=\mathit{E}\cos (\omega t)$$

(6)

The solution of equation (6) is considered as follows

$$x=\text{a}\cos \left(\omega t\right)+b\sin \left(\omega t\right)+a_3\cos \left(3\omega t\right)+b_3\sin \left(3\omega t\right)$$

(7)

Substituting equation (7) into equation (6) and equating the coefficients of equal harmonics from both sides and ignoring the terms whose responses are negligible, the following equations were found

$$a\left(-{\omega }^2+1\right)+b\mu \omega +\frac{3\epsilon a}{4}\left(a^2+a{a}_3+b^2+2a_3^2+2b{b}_3-\frac{a_3{b}^2}{a}\right)=E$$

(8a)

$$b\left(-{\omega }^2+1\right)-a\mu \omega +\frac{3\epsilon b}{4}\left(a^2-2a{a}_3+b^2-b{b}_3+\frac{a^2{b}_3}{b}\right)=0$$

(8b)

$$a_3\left(-9{\omega }^2+1\right)+3b_3\mu \omega +\frac{\epsilon a}{4}\left(a^2+6a{a}_3-3b^2+6\frac{b^2{a}_3}{a}\right)=0$$

(8c)

$$b_3\left(-9{\omega }^2+1\right)-3a_3\mu \omega +\frac{\epsilon b}{4}\left(3a^2-b^2+\frac{6a^2{b}_3}{b}+6b{b}_3\right)=0$$

(8d)

Eliminating ${\omega }^2$ from equations (8b–d) with the help of equation (8a), and neglecting the terms whose responses are negligible, the above equations can be written as follows

$$4bE-9\epsilon a^2{a}_{3}b+3\epsilon a_3{b}^3+3\epsilon a^3{b}_3-9\epsilon a{b}^2{b}_3-4\left(a^2+b^2\right)\mu \omega =0$$

(9a)

$$-32a{a}_3+36a_{3}E+\epsilon a^4-21\epsilon a^3{a}_3-3\epsilon a^2{b}^2-21\epsilon a{a}_3{b}^2+(3a{b}_3-9a_{3}b)\mu \omega =0$$

(9b)

$$-32a{b}_3+36b_{3}E+3\epsilon a^{3}b-\epsilon a{b}^3-21\epsilon a^3{b}_3-21\epsilon a{b}^2{b}_3-(3a{a}_3+9b{b}_3)\mu \omega =0$$

(9c)

Now using equation (9a), eliminating $\omega$ from equations (9c and d) and considering only linear terms of $a_3$ , $b_3$ , and ignoring the terms whose responses are negligible, the following equations are obtained

$$-32a{a}_3+36a_{3}E+\epsilon a^4-21\epsilon a^3{a}_3-2\epsilon a^2{b}^2-42\epsilon a{a}_3{b}^2-3\epsilon b^4=0$$

(10a)

$$-32a{b}_3+36b_{3}E+3\epsilon a^{3}b+2\epsilon a{b}^3-21\epsilon a^3{b}_3-42\epsilon a{b}^2{b}_3-\epsilon b^5/a=0$$

(10b)

Solving equations (10a and b), $a_3$ and $b_3$ are obtained as follows

$$a_3=\frac{\epsilon (a^4-2a^2{b}^2-3b^4)}{32a-36E+21\epsilon a^3+42\epsilon a{b}^2},b_3=\frac{\epsilon (a^{4}b+2a^2{b}^3-b^5)}{a(32a-36E+21\epsilon a^3+42\epsilon a{b}^2)}$$

(11)

Substituting the values of $a_3$ and $b_3$ into equation (9a), and then expanding in powers of the small parameter $\lambda$ as follows

$$b=l_0+l_1\lambda +l_2{\lambda }^2+l_3{\lambda }^3+\cdots$$

(12)where

$$\lambda =\frac{\mu \omega }{E},l_0=\frac{a^2\mu \omega }{E},l_1=\frac{a^4{\mu }^2{\omega }^2}{E^2},l_2=\frac{2a^6{\mu }^3{\omega }^3}{E^3},l_3=\frac{5a^8{\mu }^4{\omega }^4}{E^4},\cdots$$

Finally, after substituting $a_3$ , $b_3$ , and $b$ into equation (8a) and solving, the value of $a$ is determined.

Example 2

Consider a damped forced vibration equation of the form

$$\ddot{x}+x+\mu \dot{x}+\epsilon x^2=\mathit{E}\cos (\omega t)$$

(13)

The solution of equation (13) is taken as follows

$$x=a_0+a_1\cos \left(\omega t\right)+b_1\sin \left(\omega t\right)+a_2\cos \left(2\omega t\right)+b_2\sin \left(2\omega t\right)$$

(14)

Substituting equation (14) into equation (13) and equating the coefficients of equal harmonics from both sides and neglecting the terms whose responses are negligible, the following equations are found

$$\epsilon a_0^2+a_0+\frac{1}{2}\epsilon (a_1^2+b_1^2)=0$$

(15a)

$$a_1\left(-{\omega }^2+1\right)+b_1\mu \omega +\epsilon (2a_0{a}_1+a_1{a}_2+b_1{b}_2)=E$$

(15b)

$$b_1\left(-{\omega }^2+1\right)-a_1\mu \omega +\epsilon \left(2a_0{b}_1-a_2{b}_1+a_1{b}_2\right)=0$$

(15c)

$$a_2\left(-4{\omega }^2+1\right)+2b_2\mu \omega +\frac{1}{2}\epsilon \left(a_1^2+2a_0{a}_2-b_1^2\right)=0$$

(15d)

$$b_2\left(-4{\omega }^2+1\right)-2a_2\mu \omega +\epsilon \left(a_1{b}_1+2a_0{b}_2\right)=0$$

(15e)

Eliminating $\omega$ from equations (15c–e) and neglecting the terms whose responses are negligible, from equations (15a, c, and d) one can obtain the following

$$\begin{array}{c}{a}_0=-\frac{\epsilon a_1^2}{4}(\left(2+{\epsilon }^2{a}_1^2+{\epsilon }^4{a}_1^4\right)+a_1^2\left(2+2{\epsilon }^2{a}_1^2+3{\epsilon }^4{a}_1^4\right){\lambda }^2+a_1^4\left(4+5{\epsilon }^2{a}_1^2+9{\epsilon }^4{a}_1^4\right){\lambda }^4\\ +a_1^6\left(5+7{\epsilon }^2{a}_1^2+14{\epsilon }^4{a}_1^4\right){\lambda }^6)\end{array}$$

(16a)

$$b_1=a_1^2\lambda +a_1^4{\lambda }^3+2a_1^6{\lambda }^5$$

(16b)

$$a_2=\frac{\epsilon (a_1^4-b_1^4)}{2a_1(3a_1-4E+6\epsilon a_0{a}_1)}$$

(16c)

$$b_2=\frac{\epsilon (a_1^2{b}_1+b_1^3)}{3a_1-4E+6\epsilon a_0{a}_1}$$

(16d)where

$$\lambda =\frac{\mu \omega }{E}$$

Finally, after substituting $a_0$ , $a_2$ , $b_1$ , and $b_2$ into equation (15b) and solving, the value of $a_1$ is determined.

Solution obtained by the variational approach

Pierre de Fermat (1601–1665) ³⁹ developed the variational principles in physics which played a key role in the variational principle. Furthermore, Gottfried Leibniz (1646–1716) made significant contributions to the development of variational principles in classical mechanics. In the past few decades, variational principles have been used for various differential equations. Several methods are used to obtain the variational principles such as Lagrange multiplier,^40,41 Semi-inverse, ³³ Noether’s theorem, ⁴² etc. Earlier, the Lagrange multiplier was the most convenient way to develop generalized variational principles. In 1997, He ⁴³ proposed a more effective way to establish the generalized variational principles using the semi-inverse method. Later, the method became popular to the researcher and was used to solve various types of differential equations. In this section, a variational approach based on the semi-inverse method ³³ to solve the nonlinear forced vibration problem with damping is presented as follows.

Consider a damped forced vibration equation of the form

$$\ddot{x}+x+\mu \dot{x}+\epsilon x^3=\mathit{E}\cos (\omega t)$$

(17)

The variational parameter of equation (17) can be written as^33,36

$$J\left(x\right)={\displaystyle {\int }_{\hspace{0.17em}}^T{e}^{\mu t}(-\frac{1}{2}{\dot{x}}^2+\frac{1}{2}{x}^2+\frac{1}{4}\epsilon x^4-\mathit{xE}\cos (\omega t))\text{d}t}$$

(18)where $T=2\pi /\omega$ is the period of oscillation.

The second approximation solution of equation (17) is chosen as follows

$$x=\text{A}\cos \left(\omega t\right)+\mathit{B}\sin \left(\omega t\right)+C\cos \left(3\omega t\right)+\mathit{D}\sin \left(3\omega t\right)$$

(19)where $\omega$ is frequency, and A, B, C, and D are constants to be determined.

Substituting equation (19) into equation (18) results in

$$\begin{array}{l}J(A,B,C,D)={\displaystyle {\int }_{\hspace{0.17em}}^T{e}^{\mu t}(-\frac{1}{2}{\omega }^2{(-\text{A}\sin \left(\omega t\right)+B\cos \left(\omega t\right)-3C\sin \left(3\omega t\right)+3D\cos \left(3\omega t\right))}^2}\\ +\frac{1}{2}{\left(\text{A}\cos \left(\omega t\right)+\mathit{B}\sin \left(\omega t\right)+C\cos \left(3\omega t\right)+\mathit{D}\sin \left(3\omega t\right)\right)}^2\\ +\frac{1}{4}\epsilon {\left(\text{A}\cos \left(\omega t\right)+\mathit{B}\sin \left(\omega t\right)+C\cos \left(3\omega t\right)+\mathit{D}\sin \left(3\omega t\right)\right)}^4\\ -(\text{A}\cos \left(\omega t\right)+\mathit{B}\sin \left(\omega t\right)+C\cos \left(3\omega t\right)+\mathit{D}\sin \left(3\omega t\right))\mathit{E}\cos (\omega t))\text{d}t\end{array}$$

(20)The values of A, B, C, and D are obtained by solving the following four equations

$$\frac{\partial J(A,B,C,D)}{\partial A}=0$$

(21a)

$$\frac{\partial J(A,B,C,D)}{\partial B}=0$$

(21b)

$$\frac{\partial J(A,B,C,D)}{\partial C}=0$$

(21c)

$$\frac{\partial J(A,B,C,D)}{\partial D}=0$$

(21d)

Results and discussion

In this section, a comparison among the results obtained by the proposed method, numerical method, and a variational approach has been presented and graphically presented in Figure 1(a) to (d) and Figure 2(a) and (b). From the figures, it is observed that our result agrees reasonably well with those obtained by the variational approach and numerical method.

Figure 1.

(a) Comparison between the time versus displacement results obtained by the proposed method, numerical method, and variational approach when $\omega =5, \epsilon =.5, \mu =.05, \text{ and } E=5$ . (b) Comparison between the time versus displacement results obtained by the proposed method, numerical method, and variational approach when $\omega =5, \epsilon =1, \mu =0.1, \text{and }E=10$ . (c) Comparison between the time versus displacement results obtained by the proposed method, numerical method, and variational approach when $\omega =5, \epsilon =1, \mu =0.2, \text{and }E=20$ . (d) Comparison between the time versus displacement results obtained by the proposed method, numerical method, and variational approach when $\omega =5, \epsilon =2, \mu =0.2, \text{and }E=20$ .

Figure 2.

(a) Comparison between the time versus displacement results obtained by the proposed method, numerical method, and variational approach when $\omega =3, \epsilon =.5, \mu =.1, \text{and }E=5$ . (b) Comparison between the time versus displacement results obtained by the proposed method, numerical method, and variational approach when $\omega =3, \epsilon =.5, \mu =.1, \text{and }E=10$ .

Moreover, both solutions have been compared in phase plane together with the numerical solution in Figure 3(a) and (b). These figures also indicate that our solution shows good agreement with those calculated by the variational approach and numerical method.

Figure 3.

(a) Comparison between the phase plane obtained by the proposed method and numerical method when $\omega =5, \epsilon =.5, \mu =0.05, \text{and }E=5$ . (b) Comparison between the phase plane obtained by the proposed method and numerical method when $\epsilon =1, \mu =0.1, \text{and }E=10.$

Conclusion

In this paper, a modified HBM has been presented. The main advantage of the proposed method is that only one nonlinear algebraic equation together with a set of linear algebraic equations have been solved with less effort. However, in the classic HBM and variational approach, a set of nonlinear algebraic equations is solved by the numerical method which needs more computational effort than the proposed method.