{'fraction_non_alphanumeric': 0.028189700099493058, 'fraction_numerical': 0.016108400056853177, 'mean_word_length': 4.920897615708275, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multi-level residue harmonic balance method for nonlinear vibration of the beam

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, 1 nonlinear vibration of a curved beam with quadratic and cubic nonlinearities, 2 large deflection of a simply supported beam due to pure bending moment, 3 nonlinear dynamics of an axially moving viscoelastic beam, 4 nonlinear dynamics of a buckled beam subjected to primary resonance, 5 theoretical study of and experiments on a buckle beam, 6 and nonlinear vibrations and stability of an axially moving Timoshenko beam. 7 Nonlinear vibration and stability of duffing oscillators have been enormously investigated. [8][9][10] Lin et al. 11 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 12 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 13 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][15][16] and incremental harmonic balance method (IHB). 17 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. 18 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 19 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 20 developed an accurate frequency-amplitude relationship for describing the properties of the oldest device Fangzhu. He and Sedighi 21 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. 22 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. 23 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 24 is as follows 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 À3 , 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), ω is the excitation frequency, κ is the base excitation magnitude, and g is the acceleration due to gravity = 9.81 ms À2 . The associated boundary conditions for a simply supported beam are 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 where q i = the modal amplitude of the ith mode, φ i = the ith structural mode shape, φ i ðxÞ ¼ sinðnπx=LÞ, 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 where φ 0 i ,φ} i ,φ}} 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 φ m by each term on the right-hand side of equation (4) and taking integration over the beam length, we have For the single mode approach, equation (5) can be written as

Residue harmonic balance method
Governing equation (6) can be written as where F cosðωtÞ is an external periodic excitation term with period T ¼ 2π=ω. Since the periodic solutions are of interest, the steady-state response of equation (7) can be expanded in Fourier series 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 ≥ m. For convenience, a new time scale τ ¼ ωt=m of period 2π is introduced to eliminate ω; then equation (7) becomes The approximate harmonic balance solution of equation (8) is Substituting equation (10) into equation (9) and applying the Galerkin procedure, we have the following harmonic balance equations Carrying out the integration yields the nonlinear algebraic equation where fRg ¼ ½R c 0 ,R c 1 ,R c 2 ,…R c n ,R s 1 ,R s 2 ,…R s n T is the vector of residues to be annihilated, fLg ¼ ½2a 0 ,a 1 ,…,a n ,b 1 ,…,b n T and fFg are the Fourier coefficient vectors of the response and excitation, respectively, and ½KðL,λÞ is the total stiffness matrix. In this study, the total stiffness matrices were programmed using the formula given by Fergusson and Leung 25 ½KðL,λÞ ¼ À The unknown uðtÞ ¼ ð1=2Þa 0 þ p T z, z ¼ fcos iωtjsin iωtg, cos iωt ¼ f cos ωt cos 2ωt / g, etc. is presented by one of the following equivalent column vectors Finally, the Jacobian associated with u n ðtÞ is given by

Hamiltonian-based frequency-amplitude formulation
We consider a nonlinear oscillator in the form The variational principle of equation (13) can be written as 26 where ð1=2Þu 0 2 is the kinetic energy, and FðuÞ ¼ R f ðuÞdu is the potential energy. The total energy remains unchanged during the oscillation 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 We use the following trial function to determine the angular frequency Substituting (17) into (16), we obtain the following residual equation Similarly, we define two average residualsR A modification of the frequency-amplitude formulation is given as follows

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. 27 For the nonlinear algebraic equation (12) ½KðL,λÞfLg À m 2 Ffe mþ1 g ¼ 0 where fe mþ1 g denotes a unit column vector having unity at the ðm þ 1Þ th position and zeros elsewhere. In a unithomogeneous treatment on a continuation parameter γ using the coordinates ½a 0 ,a 1 , :::,a n ,b 1 ,b 2 , ::: , b n ,γ 2 P 2nþ1 for projection on the 2n+1 dimension plane, equation (22) becomes and we specify a compatible homogeneous structure to append an inhomogeneous linear equation α 0 a 0 þ α 1 a 1 þ / þ α n a n þ β 1 b 1 þ β n b n þ βγ ¼ 1 for some random, complex vector fα 0 , α 1 ,…,α n ,β n ,γg. 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 ≥ 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 vectorR 0 can be expressed as follows where R c i ða 0 ,a 1 ,…,a n ,b 1 ,…,b n Þ ¼ Note that the residue vector fRg in equation (12) is associated with i ¼ 0; 1,…,n while fR 0 g is associated with i ¼ n þ 1,…, 2n. The residuesR 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 residualsR 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 where ε is an embedding parameter and q j ðτÞ ¼ P ðjþ1Þn i¼0 a i,k cosðiτÞ þ P ðjþ1Þn i¼1 b i,k sinðiτÞ j ¼ 1, 2… . When q j ðτÞ is included, the number of terms is virtually j times more than what was represented in equation (10). We shall find q j ðτÞ successively in the following sub-sections. To find q 0 ðτÞ, one needs to solve the 2n + 1 nonlinear algebraic equation (12). However, when q 0 ðτÞ qualitatively describes the solution qðτÞ sufficiently, the solution can be improved quantitatively by including q j ðτÞ which requires solving 2ðj þ 1Þn þ 1 linear equations sequentially in the following manner.

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 20 . Consider the equation Operating on equation (30) by ðD 2 þ ω 2 Þ, we obtain À Equation (34) can be written as follows Substituting € q from equation (33) into equation (35) and constructing the corresponding homotopy equation has the following form The solution of equation (36) is chosen as Substituting equation (37) into equation (36) and equating the identical power of ρ to zero, we have the system of linear algebraic equation as follows The solution of equation (38)   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.

Comparison of results
The results obtained from the RHB method are qðτÞ ¼ 1:16745 × 10 À5 cosðτÞ þ 2:04985 × 10 À7 sinðτÞþ 1:17051 × 10 À7 cosð3τÞ þ 2:28473 × 10 À8 sinð3τÞ 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.

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 ð120=2πÞ < ð370=2πÞ and three coexisting solutions ð370=2πÞ < ð370=2πÞ to }. Therefore, the steady-state solutions are ð120=2πÞ < ð370=2πÞ and the lower displacement amplitude regarding the frequency ð370=2πÞ < ð370=2πÞ to }. In Figure 6, a comparison between amplitude-frequency response curves is shown to illustrate considering the damping model when the base excitation κ ¼ 0:1. The overall displacement amplitudes are calculated by using the formula j ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 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 ðω=ω 0 Þ < 1:095 and three co-existing solutions for 1:095 < ðω=ω 0 Þ < 2:053, which is termed "interval of bistability". At ðω=ω 0 Þ ¼ 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 ðω=ω 0 Þ > 1:095, and the response amplitude will provide upper branch or the large amplitude solution branch of the response diagram. When ðω=ω 0 Þ ¼ 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.
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 ðω=ω 0 Þ ¼ 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.

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.

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) received no financial support for the research, authorship, and/or publication of this article.