Nonlinear forced response of non-uniform beams carrying point masses using the harmonic balance method with an arc-length continuation scheme

Abstract

This work examines the forced vibrations of a non-uniform Euler-Bernoulli beam carrying point masses and subjected to a concentrated transverse harmonic force. By the means of Newton’s second law of motion, the integro-partial differential equation of motion is derived. The Galerkin approach is applied, and the normalized mode shapes of the non-uniform base beam are used to derive the nonlinear governing equation of motion that contains cubic and quintic nonlinearities. The method of harmonic balance method (HB) in conjunction with the pseudo arc-length continuation scheme are employed to obtain the amplitude-frequency curves. Non-uniform beams with two different boundary conditions are considered in this study; cantilever and simply supported beams. To examine the accuracy and validity of the proposed method and solution, some results obtained by the suggested approach have been compared with those found through numerical integration, where a good agreement was achieved. The influences of several factors such as the non-uniformity parameter, the magnitude of the point masses, and the boundary conditions on the steady state amplitude have been examined. It is observed that these parameters have significant effects on the dynamic behavior of the non-uniform beams. For generality and convenience, the results are displayed in dimensionless forms.

Introduction

Understanding the behavior of non-uniform structures is vital in studying the dynamic performance of several engineering applications such as helicopter rotor blades, wind turbine blades, offshore vertical risers, and others. ¹ Many researchers investigated the free vibrations of nonuniform beams. For example, Abrate ² used the Rayleigh-Ritz method to find the natural frequencies for beams with general shapes, where the variations of the cross-sectional area and the moment of inertia are represented by a polynomial of a certain degree. Wu and Chen ³ obtained closed form solutions for the natural frequencies and the mode shapes of beams with linearly tapered thickness. The modified Bessel functions of the first and second kinds were utilized, and the expansion theorem was incorporated to obtain the frequencies of the beam with point masses and different boundary conditions. Mahmoud ⁴ applied the Myklestad method to obtain the natural frequencies and the mode shapes of non-uniform and axially functionally graded cantilever Euler-Bernoulli beams carrying a tip mass. Lee and Lee ⁵ applied the transfer-matrix and the Frobenius methods to carry out the free vibration analysis of tapered Euler-Bernoulli beams. It was found that the natural frequencies and the mode shapes are affected by the taper ratio.

It is known that motors and engines are mounted on shafts and beams. Hence, their weights, sizes, and locations will have significant impact on these structures. As a result, it is important to analyze the dynamic behavior of structural elements that carry different masses at several locations. ⁶ In this regard, Hamdan and Jubran ⁷ analyzed the free and forced vibrations of a uniform Euler-Bernoulli beam with elastically restrained ends and carrying a mass with a rotary inertia. The beam was subjected to a concentrated harmonic force at the mass location. The eigen value problem was formulated, and the exact mode shapes were obtained. The Galerkin approach was utilized to derive the ordinary differential equation. A parametric study was conducted to investigate the impacts of the location of the mass and the rotary inertia, and the stiffnesses of the springs at both ends of the beam, on the resonance response of the beam. Saito et al. ⁸ carried out the nonlinear vibrations of a simply supported beam with a point mass and subjected to a concentrated harmonic force at an arbitrary location. The influence of gravity was considered, and the Galerkin approach was applied to derive the equation of motion. The harmonic balance method was utilized to find the dynamic behavior of the beam. The influence of the magnitude and the location of the mass on the natural frequency were discussed.

Özkaya et al. ⁹ considered the free and forced nonlinear vibrations of an Euler-Bernoulli beam carrying a point mass. The mode shapes and the transcendental equations for different set of boundary conditions were found. The method of multiple time scales was adopted to determine the steady state vibration amplitude as a function of the detuning parameter and the natural frequency. The effects of the different boundary conditions, the position and magnitude of the mass on the natural frequencies and the forced response were discussed. Fang et al. ¹⁰ used analytical solutions to deal with the free and forced oscillations of a double beam system having a tip mass with rotary inertia at the primary beam. The influences of several parameters such as the stiffness of the layer that connects the beams, the mass ratio, and the rigidity ratio of the beams on the system’s oscillation behavior were analyzed.

Alimoradzadeh et al. ¹¹ addressed the free and forced vibrations of Euler-Bernoulli beams lying on nonlinear viscoelastic layers. The beams are subjected to a moving mass and a thermal load. The governing equation of motion was established by the means of Hamilton’s principle, finite strain theory, and Galerkin approach. The method of multiple scales was introduced to obtain approximate solutions for the superharmonic and subharmonic responses of the beams. A parametric study was conducted to examine the impacts of the temperature rise, the coefficients of the viscoelastic layer, and the magnitude and velocity of the moving load on the dynamic behavior of the system. Karimi and Ziaei-Rad ¹² analyzed the impact of a moving mass on the nonlinear coupled vibrations of a beam with moving support. It was assumed that the mass may be moving with a constant or a variable speed. Newton’s second law of motion was employed to derive the coupled partial differential equations of motion, and the Galerkin method was utilized to obtain the governing coupled ordinary differential equations that were solved numerically. The effects of the motion of the support and the velocity of the mass on the response of the beam were reported.

Tian et al. ¹³ considered the nonlinear vibrations of a beam bridge under a moving point mass. A smooth and discontinuous absorber was placed at the center of the beam to reduce the vibrations of the system. The governing equation was introduced, and the Galerkin method was utilized to obtain the ordinary differential equation of motion. The harmonic balance method, the Fourier expansion, the complete elliptic integrals, and the power flow approach were employed to obtain approximate solutions for the amplitude of beam’s vibrations. The influences of the speed of the moving mass, the mass ratio between the mass and the beam, and the smoothness parameter of the absorber on the beam’s amplitude were illustrated.

Based on the above literature review, it is observed that carrying out the nonlinear forced response of non-uniform beams attached to concentrated masses is not yet fully addressed. Thus, the purpose of this research is to fill this gap in literature. The main objective of this article is to investigate the dynamic behavior of non-uniform beams carrying concentrated masses and subjected to transverse harmonic forces. To analyze the proposed nonlinear system, we will use the harmonic balance method (HB), the Galerkin approach, and the arc-length continuation method to obtain the amplitude-frequency responses in the primary resonance state. The mathematical model of the non-uniform beam is formulated by utilizing Newton’s second law of motion and the Euler-Bernoulli beam theory. The effects of the magnitudes of the point masses, the boundary conditions, and the non-uniformity parameter on the dynamic characteristics of the non-uniform beam are presented in frequency response curves, phase portrait, and time response diagrams. The reliability of the HB method was illustrated by comparing the results with those obtained through direct numerical integration and good agreement is achieved. From the physical viewpoint, the non-uniform beam with the attached concentrated masses can be used as a model of structural elements used in aerospace and mechanical engineering applications carrying motors or engines. Moreover, the suggested model and solution procedure are appropriate for the establishment of the amplitude-frequency diagrams of strongly nonlinear systems.

Mathematical modeling

As previously mentioned in the introduction, the nonuniform beams can be used to represent tapered wings, fuselages, and helicopter rotor blades.^14,15 Furthermore, a point mass can be used to model a wide range of examples in engineering where concentrated loads or localized masses significantly influence the structure’s behavior.^16,17 For example, a point mass can be used to represent a motor or a gripper at the end of a robotic arm. Additionally, in aircraft wings, the engine can be represented as a point mass. In crane beams supporting a trolley and payload, the moving trolley and payload may be represented as point masses. In spacecraft booms supporting solar panels with attached instruments or equipment, the instruments can be modeled as point masses.

Following Rincón-Casado et al., ¹⁸ the governing partial differential equation of motion for free transverse vibrations of an Euler-Bernoulli beam is expressed as

$${\left[\frac{\mathrm{EI}{\left(\frac{w″}{1-\frac{{w\prime }^2}{2}}\right)}^{\prime }}{1-\frac{{w\prime }^2}{2}}\right]}^{\prime }+\rho A\stackrel{··}{w}-\left[\frac{\mathrm{Hw}\prime }{1-\frac{{w\prime }^2}{2}}\right]\prime =0$$

(1)

where E is the modulus of elasticity, I is the second moment of area, w is the transverse displacement, A is the cross-sectional area of the beam, $\rho$ is the density of the beam’s material, H is the horizontal projection of the internal forces (normal and shear forces) onto the global x and y plane. The prime indicates the partial derivative with respect to $x$ , whereas the dot indicates the derivative with respect to time. In the current study, it is assumed that the beam is subjected to a concentrated harmonic force. Besides, point masses are attached to the beam, and the internal structural damping of the beam’s material is considered. Applying Newton’s second law of motion by summing the forces in the transverse direction (the inertia force due to the point masses, the damping force, and the concentrated harmonic force), considering equation (1) and expanding the first term by finding its full derivative, the equation of motion is modified as

$$\begin{array}{l}E(I^{\prime \prime }{w}^{\prime \prime }+I^{\prime \prime }{w^{\prime }}^2{w}^{\prime \prime }+\frac{1}{4}{I}^{\prime \prime }{w^{\prime }}^4{w}^{\prime \prime }+3I^{\prime }{w}^{\prime }{w}^{\prime \prime 2}+\frac{3}{2}{I}^{\prime }{w^{\prime }}^3{w}^{\prime \prime 2}+I{w}^{\prime \prime 3}+\frac{3}{2}I{w^{\prime }}^2{w}^{\prime \prime 3}+2I^{\prime }{w}^{\prime \prime \prime }+2I^{\prime }{w^{\prime }}^2{w}^{\prime \prime \prime }\\ +\frac{1}{2}{I}^{\prime }{w^{\prime }}^4{w}^{\prime \prime \prime }+4I{w}^{\prime }{w}^{\prime \prime }{w}^{\prime \prime \prime }+2I{w^{\prime }}^3{w}^{\prime \prime }{w}^{\prime \prime \prime }+I{w}^{\prime \prime \prime \prime }+I{w^{\prime }}^2{w}^{\prime \prime \prime \prime }+\frac{1}{4}I{w^{\prime }}^4{w}^{\prime \prime \prime \prime })\\ +\rho A\stackrel{··}{w}-{\left(H{w}^{\prime }\right)}^{\prime }+c_{v}EI{\stackrel{.}{w}}^{\prime \prime \prime \prime }+\sum _{i=1}^n{m}_i\delta \left(x-x_i\right)\stackrel{··}{w}=F_0\cos \left(\mathrm{\Omega }t\right)\delta \left(x-x_f\right)\end{array}$$

(2)

where n is the number of the attached point masses, $x_i$ is the location of the mass $m_i$ , $x_f$ is the point of application of the harmonic force, $\delta ()$ is the Dirac delta function, $F_0$ is magnitude of the transverse force, $\mathrm{\Omega }$ is the excitation frequency, t is time, and

$$(Hw\prime )\prime =-\rho (w\prime {\int }_l^{s}A{\int }_0^s({\stackrel{··}{w}}^{″}w\prime +{\stackrel{·}{w}\prime }^2)\mathrm{ds}\mathrm{ds})\prime$$

(3)

For more details regarding the derivation of the integro-differential equation (equation (2)), interested readers are encouraged to refer to Rincón-Casado et al. ¹⁸ where a step-by-step formulation was shown for the case of a uniform beam without point masses. For brevity, the formulation of equation (2) was omitted in the current study.

Figure 1 shows a non-uniform cantilever beam with attached point masses and acted upon by a tip harmonic force. In this study, it is assumed that the thickness of the beam is kept constant, and the width of the beam varies exponentially along the length. Thus, the cross-sectional area and the moment of inertia of the beam are given as ¹⁹

$$A=A\left(x\right)=A_0{e}^{\Lambda x}$$

(4)

$$I=I\left(x\right)=I_0{e}^{\Lambda x}$$

(5)

where $A_0$ and $I_0$ are the area and the moment of inertia at the left end of the beam ( $x=0$ ), respectively, and $\mathrm{\Lambda }$ is the non-uniformity parameter. In Ece et al., ¹⁹ the non-uniformity parameter was taken in the range of $-2<\mathrm{\Lambda }<2$ . For convenience and better analysis of the dynamic behavior and numerical computations, the equation of motion is transformed into a nondimensional form utilizing the following dimensionless parameters

$$x^{*}=\frac{x}{L},w^{*}=\frac{w}{L},t^{*}=t\sqrt{\frac{E{I}_0}{\rho A_0{L}^4}},{\mathrm{\Omega }}^{*}=\mathrm{\Omega }\sqrt{\frac{\rho A_0{L}^4}{E{I}_0}},F_0^{*}=\frac{F_0{L}^3}{E{I}_0},c=c_v\sqrt{\frac{E{I}_0}{\rho A_{0}L},}{\mu }_i=\frac{m_i}{M}$$

(6)

where M is the total mass of the beam given as

$$M={\int }_0^1\rho A(x)\mathrm{dx}$$

(7)

Figure 1.

Schematic diagram of a non-uniform cantilever beam carrying n point masses and subjected to a reciprocating tip force.

These dimensionless parameters are substituted into equation (2) to convert the governing equation of motion into a non-dimensional form. It should be stated that the asterisk notations are removed for simplicity. The governing equation of motion is obtained as

$$\begin{array}{c}(I^{″}{w}^{″}+I^{″}{w\prime }^2{w}^{″}+\frac{1}{4}{I}^{″}{w\prime }^4{w}^{″}+3I\prime w\prime {w^{″}}^2+\frac{3}{2}I\prime {w\prime }^3{w^{″}}^2+I{w^{″}}^3+\frac{3}{2}I{w\prime }^2{w^{″}}^3+2I\prime w^{″\prime }+2I\prime {w\prime }^2{w}^{″\prime }\hfill \\ +\frac{1}{2}I\prime {w\prime }^4{w}^{″\prime }+4Iw\prime w^{″}{w}^{″\prime }+2I{w\prime }^3{w}^{″}{w}^{″\prime }+I{w}^{″″}+I{w\prime }^2{w}^{″″}+\frac{1}{4}I{w\prime }^4{w}^{″″})\hfill \\ +A\stackrel{··}{w}-(Hw\prime )\prime +\mathrm{cI}{\stackrel{·}{w}}^{″″}+\sum _{i=1}^n{\mu }_i\delta (x-x_i)\stackrel{··}{w}=F_0\cos \left(\mathrm{\Omega }t\right)\delta (x-x_f)\hfill \end{array}$$

(8)

Solution procedure

The nonlinear partial differential equation (equation (8)) can be discretized in space and time utilizing the Galerkin’s approach as ²⁰

$$w(x,t)=\sum _{j=1}^N{\varphi }_j\left(x\right)q_j\left(t\right)$$

(9)

where N denotes the number of retained modes, ${\varphi }_j\left(x\right)$ are the linear mode shapes for a cantilever or a simply supported base beam, and $q_j\left(t\right)$ are the generalized coordinates. In the current study, it is assumed that the concentrated harmonic force will only excite the fundamental mode.^13,21 Moreover, it is supposed that this is the dominant mode for the beam, and higher modes can be ignored. Furthermore, it is presumed that this mode is not involved in an internal resonance with any other higher modes. ²² Therefore, other modes not being directly or indirectly excited, will decay to zero with time due to the existence of damping. Thus, the beam bending is given as

$$w(x,t)\approx {\varphi }_1\left(x\right)q\left(t\right)$$

(10)

where the first bending mode shape of the base beam (without the point masses) is given as ¹⁹

$${\varphi }_1\left(x\right)=B_2{e}^{\frac{-\Lambda }{2}x}[b_1\cos \left({\lambda }_{1}x\right)+\sin \left({\lambda }_{1}x\right)-b_1\cosh \left({\lambda }_{2}x\right)+b_4\sinh \left({\lambda }_{2}x\right)]$$

(11)

The relations that govern $b_1$ , $b_4$ , ${\lambda }_1$ , ${\lambda }_2$ , and the natural frequency $\omega$ are found in Ece et al. ¹⁹ The constant $B_2$ is determined by normalizing the mode shape of the base beam such that

$${\int }_0^1{\varphi }_1^2\left(x\right)\mathrm{dx}=1$$

(12)

Equations (9) and (10) are substituted into equation (8), then the result is multiplied by ${\varphi }_1\left(x\right)$ and integrated over the length of the beam to obtain the nonlinear ordinary differential equation that is expressed as

$$m_{\mathrm{eq}}\stackrel{··}{q}+c\beta \stackrel{·}{q}+{\alpha }_{1}q+{\alpha }_2(q{\stackrel{·}{q}}^2+q^2\stackrel{··}{q})+{\alpha }_3{q}^3+{\alpha }_5{q}^5=F_0{\varphi }_1\left(x_f\right)\cos \left(\mathrm{\Omega }t\right)$$

(13)

where

$$m_{\mathrm{eq}}={\int }_0^1{\varphi }_1^2\left(x\right)e^{\Lambda x}\mathrm{dx}+\sum _{i=1}^n\left({\mu }_i{\varphi }_1^2\left(x_i\right)\right){\int }_0^1{e}^{\Lambda x}\mathrm{dx}$$

(14)

$$\beta ={\int }_0^1{\varphi }_1\left(x\right)e^{\Lambda x}\left(x\right)\varphi ″{\prime }_1\left(x\right)\mathrm{dx}$$

(15)

$${\alpha }_1={\int }_0^1{\varphi }_1\left(x\right)(e^{\mathrm{\Lambda }x}{\mathrm{\Lambda }}^2{\varphi }_1^{″}\left(x\right)+2e^{\mathrm{\Lambda }x}\mathrm{\Lambda }{\varphi }_1^{″\prime }\left(x\right)+e^{\mathrm{\Lambda }x}{\varphi }_1^{″″}\left(x\right))\mathrm{d}x$$

(16)

$${\alpha }_2={\int }_0^1{\varphi }_1\left(s\right)[{\varphi }_1\prime \left(s\right){\int }_1^{s}A\left(s\right){\int }_0^{s}2{{\varphi }_1\prime }^2\left(s\right)\mathrm{ds}\mathrm{ds}]\prime \mathrm{ds}$$

(17)

$$\begin{array}{l}{\alpha }_3={\int }_0^1{\varphi }_1\left(x\right)[e^{\Lambda x}{\Lambda }^2{\varphi }_1^{\prime 2}\left(x\right){\varphi }_1^{\prime \prime }\left(x\right)+3e^{\Lambda x}\Lambda {\varphi }_1^{\prime }\left(x\right){\varphi }_1^{\prime \prime 2}\left(x\right)+e^{\Lambda x}{\varphi }_1^{\prime \prime 3}\left(x\right)+2e^{\Lambda x}\Lambda {\varphi }_1^{\prime 2}\left(x\right){\varphi }_1^{\prime \prime \prime }\left(x\right)\\ +4e^{\Lambda x}{\varphi }_1^{\prime }\left(x\right){\varphi }_1^{\prime \prime }\left(x\right){\varphi }_1^{\prime \prime \prime }\left(x\right)+e^{\Lambda x}{\varphi }_1^{\prime 2}{\varphi }_1^{\prime \prime \prime \prime }\left(x\right)]\mathrm{d}x\end{array}$$

(18)

$$\begin{array}{c}{\alpha }_5={\int }_0^1{\varphi }_1\left(x\right)(\frac{1}{4}{e}^{\mathrm{\Lambda }x}{\mathrm{\Lambda }}^2{\varphi }_1{\prime }^4\left(x\right){\varphi }_1^{″}\left(x\right)+\frac{3}{2}{e}^{\mathrm{\Lambda }x}\mathrm{\Lambda }{\varphi }_1{\prime }^3{\varphi }_1^2+\frac{3}{2}{e}^{\mathrm{\Lambda }x}{\varphi }_1{\prime }^2{\varphi }_1^3+\hfill \\ \frac{1}{2}{e}^{\mathrm{\Lambda }x}\mathrm{\Lambda }{\varphi }_1{\prime }^4\left(x\right){\varphi }_1^{″\prime }\left(x\right)+2e^{\mathrm{\Lambda }x}{\varphi }_1{\prime }^3\left(x\right){\varphi }_1^{″}\left(x\right){\varphi }_1^{″\prime }\left(x\right)+\frac{1}{4}{e}^{\mathrm{\Lambda }x}{\varphi }_1{\prime }^4\left(x\right){\varphi }_1^{″″}\left(x\right))\mathrm{d}x\hfill \end{array}$$

(19)

$x_f=1$ for the cantilever beam that is subjected to a harmonic tip force, whereas $x_f=0.5$ for the simply supported and the clamped-clamped beams that are acted upon by harmonic center forces. The natural frequency is given as

$${\omega }_n=\sqrt{\frac{{\alpha }_1}{m_{\mathrm{eq}}}}$$

(20)

The cantilever, simply supported, and clamped-clamped beams are considered in the present study. The boundary conditions for these beams are given in terms of the mode shapes as ¹⁹

Cantilever beam:

$$\varphi \left(0\right)=\frac{\partial \varphi \left(0\right)}{\partial x}=\frac{{\partial }^2\varphi \left(1\right)}{\partial x^2}=\frac{{\partial }^3\varphi \left(1\right)}{\partial x^3}=0$$

(21)

Simply supported beam:

$$\varphi \left(0\right)=\frac{{\partial }^2\varphi \left(0\right)}{\partial x^2}=\varphi \left(1\right)=\frac{{\partial }^2\varphi \left(1\right)}{\partial x^2}=0$$

(22)

Clamped-clamped beam:

$$\varphi \left(0\right)=\frac{\partial \varphi \left(0\right)}{\partial x}=\varphi \left(1\right)=\frac{\partial \varphi \left(1\right)}{\partial x}=0$$

(23)

The normalized mode shapes for non-uniform cantilever, simply supported, and clamped-clamped beams are expressed as

Cantilever beam:

• a. For the non-uniformity parameter $\mathrm{\Lambda }$  = 1

$$\varphi \left(x\right)=0.868282e^{-x/2}(-1.8636\cos (1.5216x)+1.8636\cosh (1.6779x)+\sin (1.5216x)-0.90686\sinh (1.6779x))$$

(24)

• b. For the non-uniformity parameter $\mathrm{\Lambda }$  = 2

$$\varphi \left(x\right)=0.765143e^{-x}(-3.43074\cos (0.916842x)+3.43074\cosh (1.68541x)+\sin (0.916842x)-0.5439884\sinh (1.68541x))$$

(25)

Simply Supported beam:

• a. For the non-uniformity parameter $\mathrm{\Lambda }$  = 1

$$\varphi \left(x\right)=1.83986e^{-x/2}(-0.133581\cos (3.08592x)+0.133581\cosh (3.165898x)+\sin (3.08592x)-0.150029\sinh (3.165898x))$$

(26)

• b. For the non-uniformity parameter $\mathrm{\Lambda }$  = 2

$$\varphi \left(x\right)=2.4967e^{-x}(-0.219949\cos (2.913289x)+0.219949\cosh (3.2384024x)+\sin (2.913289x)-0.2552436\sinh (3.2384024x))$$

(27)

Clamped-clamped beam:

• a. For the non-uniformity parameter $\mathrm{\Lambda }$  = 1

$$\varphi \left(x\right)=1.265647e^{-x/2}(-1.005857\cos (4.718227x)+1.00585684\cosh (4.7709192x)+\sin (4.718227x)-0.98896\sinh (4.7709192x))$$

(28)

• b. For the non-uniformity parameter $\mathrm{\Lambda }$  = 2

$$\varphi \left(x\right)=1.63878e^{-x}(-0.97179\cos (4.68377x)+0.97179\cosh (4.89262x)+\sin (4.68377x)-0.95731\sinh (4.89262x))$$

(29)

The normalized mode shapes for the different nonuniform Euler-Bernoulli beams are displayed in Figure 2.

Figure 2.

The first mode shape for nonuniform: (a) cantilever, (b) simply supported, and (c) clamped Euler-Bernoulli beams.

Harmonic balance

The periodic solutions of the equation of motion can be obtained using the harmonic balance method. The periodic solution is assumed as

$$q\left(t\right)=a_0+\sum _{i=1}^N(a_i\cos \left(i\mathrm{\Omega }t\right)+b_i\sin \left(i\mathrm{\Omega }t\right))$$

(30)

where $a_0$ is a constant that represents the solution mean of the generalized coordinate $q\left(t\right)$ , N is the number of harmonics retained in the approximation, $a_i$ are the coefficients of the cosine terms, and $b_i$ are the coefficients of the sine terms. The assumed solution shown in equation (30) is substituted into the equation of motion, multiplied by each harmonic, then the result is integrated with respect to the time t from 0 to $2\pi /\mathrm{\Omega }$ , leading to a system of $2N+1$ nonlinear algebraic equations, given as

$$f(X,\mathrm{\Omega })=0$$

(31)

where

$$X={\left\{\begin{array}{cc}\begin{array}{ccc}{a}_0\hfill & a_1\hfill & \dots \hfill \end{array}\hfill & \begin{array}{ccc}{a}_k\hfill & b_1\hfill & \begin{array}{cc}\dots \hfill & b_k\hfill \end{array}\hfill \end{array}\hfill \end{array}\right\}}^T$$

(32)

Arc length continuation method

The Pseudo arc-length continuation method is utilized in this research to track all possible stable and unstable solutions of X (at a certain value of $\mathrm{\Omega }$ ).^23,24 Introducing the arc-length parameter $\eta$ , the augmented equation is expressed as

$$f(X,\mathrm{\Omega })-\eta =0$$

(33)

Using a quadratic form of the pseudo arc length, equation (33) is rewritten as

$$f\left(Y\right)-\eta ={Y\prime }^T(Y-Y_{k-1})-\eta =0$$

(34)

where $Y={[X,\mathrm{\Omega }]}^T$ .

The slope of the equilibrium path ( $Y^{\prime }$ ) is defined in terms of two nearby previous points $Y_{k-2}$ and $Y_{k-1}$ as

$$Y\prime =\frac{Y_{k-1}-Y_{k-2}}{\Vert Y_{k-1}-Y_{k-2}\Vert }$$

(35)

The initial prediction of the next point is determined as

$$Y_u=Y_{k-1}+\mathrm{\Delta }\eta Y\prime$$

(36)

To obtain the periodic solution, the prediction obtained in equation (36) is corrected iteratively using the Newton-Raphson method until the exact solution on the equilibrium path is found. As the governing equation of motion (equation (13)) includes odd nonlinearities only, the solution given in equation (31) can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{q\left(t\right)=a_0+a_1\cos \left(\mathrm{\Omega }t\right)+b_1\sin \left(\mathrm{\Omega }t\right)+a_3\cos \left(3\mathrm{\Omega }t\right)}\\ +b_3\sin \left(3\mathrm{\Omega }t\right)+a_5\cos \left(5\mathrm{\Omega }t\right)+b_5\sin \left(5\mathrm{\Omega }t\right)\hfill \end{array}$$

(37)

Results

In the current section, frequency response curves of non-uniform cantilever and simply supported beams with attached masses and subjected to concentrated harmonic forces are presented. Regarding the damping properties, it is known that aerospace systems often have very low damping ratios (ζ < 0.01) to ensure minimal energy loss and high precision. In machinery, the industrial equipment may have moderate damping (ζ ≈ 0.05–0.2) to reduce noise and vibrations. In the current study, the dimensionless damping coefficient c = 0.1. It is necessary to conduct a convergence test to determine the sufficient number of harmonics used in the analysis. Therefore, the steady state amplitudes are determined and plotted versus the excitation frequency for different numbers of harmonics for a non-uniform cantilever beam carrying three masses at x = 0.3, 0.6, and 0.9 for $\mathrm{\Lambda }=2.0,\mu =6.0,F_0=60$ . Figure 3 shows that three harmonics is an adequate number for convergence, so that three harmonics are used in the calculations. To examine the accuracy of the proposed solution obtained using the HB, a comparison study has been carried out. In Figure 4, the frequency amplitude curve obtained using the HB is compared to that achieved with the numerical integration approach (NI) for a non-uniform cantilever beam carrying three masses at x = 0.3, 0.6, and 0.9 for $\mathrm{\Lambda }=2.0,\mu =6.0,F_0=60$ . The solid line corresponds to the solution obtained from the HB, whereas the circles correspond to the results generated using the NI method. The accuracy of the suggested method can be validated by noting that the solutions of both approaches are in good agreement.

Figure 3.

The frequency response of q obtained using the harmonic balance method (HB), for different numbers of harmonics, of a non-uniform cantilever beam carrying three masses at x = 0.3, 0.6, and 0.9 for $\mathrm{\Lambda }=2.0,\mu =6.0,F_0=60$ .

Figure 4.

The frequency response of q obtained using the harmonic balance method (HB) compared with that obtained using direct numerical integration (NI) of a non-uniform cantilever beam carrying three masses at x = 0.3, 0.6, and 0.9 for $\mathrm{\Lambda }=2.0,\mu =6.0,F_0=60$ .

In Figure 5, the frequency amplitude curves are obtained for a non-uniform cantilever beam carrying three equal masses at x = 0.3, 0.6, and 0.9 for $\mathrm{\Lambda }=1.0,F_0=60$ at several values of the mass magnitudes. One can observe that on the upper stable branches, the maximum displacement of the beam rises as the magnitudes of the attached masses grow. It is known that on the upper stable branch where the driving frequency Ω is increased starting from a small value, the steady state amplitude is increased until Ω reaches the value where the periodic solution loses its stability due to the appearance of a saddle-node bifurcation. ²⁵ An additional rise in Ω leads to a sharp drop in the magnitudes of the steady state amplitude. In general, as Ω gets closer to the natural frequency ${\omega }_n$ of the system, the displacement grows, and as it moves away from ${\omega }_n$ , the displacement declines. It should be noted that adding masses to the beam results in decreasing ${\omega }_n$ of the system, whereas increasing the non-uniformity parameter may lead to a grow or a decline in ${\omega }_n$ depending on the boundary conditions of the beam.

Figure 5.

The influence of the magnitudes of the point masses on the frequency-response curves of a non-uniform cantilever beam carrying three masses at x = 0.3, 0.6, and 0.9 for $\mathrm{\Lambda }=1.0,F_0=60$ .

In contrary, on the lower stable branches, the steady state peak amplitude declines as the magnitudes of the attached points rise. Thus, the influence of the mass magnitudes on the steady state amplitudes is different as the forward and backward sweeping are conducted. Besides, at lower excitation frequencies, it is noted that loops are created where multiple periodic solutions are found. Moreover, the hardening nonlinearity performance becomes stronger as the magnitudes of the point masses increase. Similar observations are made for Figure 6. The hardening nonlinearity refers to a type of nonlinear behavior in a system where the stiffness of the system increases as the displacement rises, causing an upward bend in the response curve at higher frequencies. This means the restoring force grows more rapidly than linearly with displacement, leading to an increase in the natural frequency of the system at larger amplitudes.

Figure 6.

The influence of the magnitudes of the point masses on the frequency-response curves of a non-uniform cantilever beam carrying three masses at x = 0.3, 0.6, and 0.9 for $\mathrm{\Lambda }=2.0,F_0=60$ .

Based on Figure 4, for $\mathrm{\Omega }\ge 5.5$ , the cantilever beam will have a low and stable steady state amplitude regardless of the magnitudes of the three added point masses.

Figure 7 displays the steady state amplitude of the three harmonics considered in this work. These amplitudes are generated for a non-uniform cantilever beam carrying three masses at x = 0.3, 0.6, and 0.9 with $\mu =0.5,\mathrm{\Lambda }=1.0,\mathrm{and}{F}_0=60$ . The magnitudes of these amplitudes are given as

$$A_i=\sqrt{a_i^2+b_i^2}$$

(35)

Figure 7.

The amplitude-frequency response curves $A_1$ , $A_3$ , and $A_5$ of a non-uniform cantilever beam carrying three masses at x = 0.3, 0.6, and 0.9 for $\mu =0.5,\mathrm{\Lambda }=1.0,F_0=60$ .

It is clear that $A_1>A_3>A_5$ ; so, the contribution of the first harmonic is higher than that for the higher harmonics. Additionally, the amplitude-frequency response curves $A_3$ and $A_5$ are similar to each other, whereas the amplitude-frequency response curve $A_1$ looks similar to the total amplitude A that is shown in the previous figures.

Figure 8 illustrates the impact of the non-uniformity parameter $\mathrm{\Lambda }$ on the steady state amplitude of a cantilever beam carrying three masses at x = 0.3, 0.6, and 0.9 for $\mu =0.5$ , and $F_0=60$ . It can be seen that as the non-uniformity parameter advances, the maximum amplitude of the beam (the peaks) declines. Furthermore, on the upper stable branches and for $\mathrm{\Omega }\le 5$ , it is noted that the steady state amplitude increases as the non-uniformity parameter grows. In contrary, on the lower stable branches, the steady-state amplitude decreases as $\mathrm{\Lambda }$ rises. Both the frequency jump-up and jump down phenomena are evident in this figure and in the previous ones. It is worthwhile to mention that when $\mathrm{\Lambda }=0$ , the uniform beam is recovered. Based on Figure 6, to get stable and very low steady state amplitudes, it is preferred to have $\mathrm{\Omega }\ge 16$ for a cantilever beam with $\mathrm{\Lambda }=0$ , $\mathrm{\Omega }\ge 9$ for a cantilever beam with $\mathrm{\Lambda }=1.0$ , and $\mathrm{\Omega }\ge 5$ for a cantilever beam with $\mathrm{\Lambda }=2.0$ . Thus, for $\mathrm{\Omega }\ge 16$ , the cantilever beam will have a very low steady state amplitude regardless of the value of $\mathrm{\Lambda }$ .

Figure 8.

The Influence of the non-uniformity parameter on the frequency-response curves of a cantilever beam carrying three masses at x = 0.3, 0.6, and 0.9 for $\mu =0.5,F_0=60$ .

The impact of the magnitudes of the attached masses on the amplitude of a non-uniform simply supported beam carrying two masses at x = 0.25 and 0.75, for $\mathrm{\Lambda }=1.0,$ and $F_0=60$ is depicted in Figure 9. In this case, the beam is subjected to a concentrated harmonic force applied at the center of the beam (at $x=\frac{1}{2}$ ). It is noticed that the maximum amplitudes of the non-uniform beam (the peaks) grow as the magnitudes of the attached masses enhance. This behavior is evident on the upper stable branches, whereas it is reversed on the lower stable branches, and the steady state amplitudes are reduced as the magnitudes of these masses advance. The hardening nonlinearity gets stronger at larger magnitudes of the masses. On the upper stable branches, the frequency keeps increasing until it reaches a point where any further development in its value will lead to a jump-down, and there will be a sharp decrease in the response amplitude. For the lower stable branches, one will observe that as the frequency decreases, it will reach a value where the periodic solution loses its stability, and where the jump-up phenomenon will occur. This action of forward and backward frequency sweeping leads to the creation of the “hysteresis” phenomenon. It is well-known that the jump-up and jump-down in the response amplitude happens at turning points where a qualitative change in the system’s behavior occurs, which is known as “bifurcation.” Jump up and jump down refer to sudden transitions in the response of a system as the driving frequency is varied. These phenomena are associated with nonlinear resonance and arise due to the system’s nonlinear stiffness or damping characteristics. The jump up happens when the system’s steady state amplitude abruptly rises to a higher state as the driving frequency is gradually increased. On the other hand, the jump down occurs when the system’s steady state amplitude suddenly decreases to a lower state as the driving frequency is gradually decreased.^23,25 It represents a transition from a higher-energy state to a lower-energy state. In general, the jumps occur due to the presence of multiple stable equilibrium points or dynamic attractors in the system. As the driving frequency changes the system follows one branch of the equilibrium curve until it reaches a critical point (e.g. saddle-node bifurcation), where it loses stability. At this critical point, the system jumps to another stable branch of the response curve. Considering Figure 9, for $\mathrm{\Omega }\ge 7.5$ , the non-uniform simply supported beam will have a low and stable steady state amplitude regardless of the magnitudes of the two added point masses. Almost identical statements are made for Figure 10.

Figure 9.

The effect of the point masses magnitudes on the frequency-response curves of a non-uniform simply supported beam carrying two masses at x = 0.25, and 0.75, for $\mathrm{\Lambda }=1.0,F_0=60$ .

Figure 10.

The impact of the magnitudes of the point masses on the frequency-response curves of a non-uniform simply supported beam carrying two masses at x = 0.25, and 0.75, for $\mathrm{\Lambda }=2.0,F_0=60$ .

The influence of the force amplitude on the frequency-response diagrams of a non-uniform simply supported beam carrying two masses at x = 0.25 and 0.75, for $\mu =0.5,$ and $\mathrm{\Lambda }=1.0$ is displayed in Figure 11. As expected, the response amplitude rises as the force amplitude increases. Besides, the growth in the force amplitude results in a change in the frequency-amplitude response and leads to the appearance of additional local periodic solutions and branches at lower values of the excitation frequency. This may be attributed due to the contributions of the super harmonics of the non-uniform beam.

Figure 11.

The impact of the force amplitude on the frequency-response curves of a simply supported non-uniform beam carrying two masses at x = 0.25, and 0.75, for $\mu =0.5,\mathrm{\Lambda }=1.0$ .

Figure 12 shows the impact of the non-uniformity parameter $\mathrm{\Lambda }$ on the steady state amplitude of a simply supported beam carrying two masses at x = 0.25 and x = 0.75 for $\mu =2.0$ , and $F_0=60$ . It can be noticed that as the non-uniformity parameter increases, the maximum amplitude of the beam (the peaks) drops. Moreover, on the upper stable branches for $\mathrm{\Omega }\le 4.2$ and on the lower stable branches, it is observed that the steady state amplitude decreases as the non-uniformity parameter grows. However, for certain values of the excitation frequency on the upper stable branches, the steady-state amplitude is larger at lower values of $\mathrm{\Lambda }$ . Besides, as $\mathrm{\Lambda }$ advances, the hardening nonlinearity becomes weaker. The frequency jump-up and jump-down take place, where the dramatic increase and decrease of the response amplitude are reduced as $\mathrm{\Lambda }$ progresses. In light of Figure 12, for $\mathrm{\Omega }\ge 9.0$ , the non-uniform simply supported beam will have a low and stable steady state amplitude despite the magnitudes of the non-uniformity parameter.

Figure 12.

The influence of the non-uniformity parameter on the frequency-response curves of a simply supported beam carrying two masses at x = 0.25, and 0.75 for $\mu =2.0,F_0=60$ .

The influences of the magnitudes of the attached masses on the amplitude of a non-uniform clamped-clamped beam carrying two masses at x = 0.25 and 0.75, $F_0=200$ , for $\mathrm{\Lambda }=1.0$ and 2.0 are presented in Figures 13 and 14, respectively. As in the simply supported beam, the concentrated harmonic force is applied at the center of the beam (at $x=\frac{1}{2}$ ). Observations like those that were made for Figures 11 and 12 are noticed. Moreover, the curves look smoother than those obtained for the cantilever and simply supported beams, and there are no dramatic declines in the magnitudes of the steady state amplitudes especially for µ = 0.5 and 2.0, where it can be stated that the jump-down and jump-up phenomena are not present.

Figure 13.

The influence of the magnitudes of the point masses on the frequency-response curves of a non-uniform clamped-clamped beam carrying two masses at x = 0.25 and 0.75 for $\mathrm{\Lambda }=1.0,F_0=200$ .

Figure 14.

The influence of the magnitudes of the point masses on the frequency-response curves of a non-uniform clamped-clamped beam carrying two masses at x = 0.25 and 0.75 for $\mathrm{\Lambda }=2.0,F_0=200$ .

Figure 15 represents the time response and phase portrait for a non-uniform cantilever beam carrying three masses at x = 0.3, 0.6, and 0.9 with $\mathrm{\Lambda }=2,\mu =6.0,F_0=60$ , and $\mathrm{\Omega }=1$ . In Figure 15(a), it is clear that the beam vibrates around $q=0$ with a periodic response. In Figure 15(b), the plot resembles a single ellipse which implies that the vibration response contains higher order harmonic components. In a similar manner, in Figures 16(a) and 17(a), the simply supported and clamped-clamped beams, respectively, oscillate around $q=0$ with a bounded periodic amplitude. As illustrated in Figures 16(b) and 17(b), the phase portraits correspond to closed loop trajectories.

Figure 15.

(a) The steady state forced response $q\left(t\right)$ and (b) the phase portrait obtained for a non-uniform cantilever beam carrying three masses at x = 0.3, 0.6, and 0.9 with $\mathrm{\Lambda }=2,\mu =6.0,F_0=60$ , and $\mathrm{\Omega }=1$ .

Figure 16.

(a) The steady state forced response $q\left(t\right)$ and (b) the phase portrait obtained for a non-uniform simply supported beam carrying two masses at x = 0.25, and 0.75 with $\mu =2.0,\mathrm{\Lambda }=1,F_0=60$ , and $\mathrm{\Omega }=5$ .

Figure 17.

(a) The steady state forced response $q\left(t\right)$ and (b) the phase portrait obtained for a non-uniform clamped-clamped beam carrying two masses at x = 0.25, and 0.75 with $\mu =2.0,\mathrm{\Lambda }=1,F_0=200$ , and $\mathrm{\Omega }=8$ .

Conclusions

The forced nonlinear response of cantilever and simply supported non-uniform Euler Bernoulli beams carrying point masses and subjected to concentrated harmonic forces have been examined. The beam’s height is kept constant, and its width is assumed to alter exponentially along the length of the beam. The governing integro-partial differential equation of motion was formulated using Newton’s second law of motion, and the Galerkin approach was utilized to convert it into an ordinary differential equation. The HB method in conjunction with the pseudo-arclength continuation scheme were applied to generate frequency-amplitude curves. Some results were validated by comparing them with those found using the direct numerical integration of the equation of motion. The effects of the non-uniformity parameter, the magnitudes of the attached masses, and the boundary conditions on the behavior of the non-uniform beams were carried out. The results reveal that these parameters have hardening impacts on the frequency-amplitude diagrams, and the jump-down and jump-up phenomena were observed. The impact of the non-uniformity parameter and the magnitudes of the attached masses on the response amplitude depend on whether the frequencies are swept on the upper or lower stable branches. It is believed that the outcomes may be useful for researchers working on examining the forced behavior of thin structures with non-uniform cross sections and carrying point masses, as it may have direct applications in designing resilient, innovative, and efficient infrastructure. For future work, it is recommended to consider other cross-sections, boundary conditions, and other effects such as thermal and aerodynamic forces.