# A Multimode Approach to Geometrically Nonlinear Free and Forced Vibrations of Multistepped Beams

Authors: Issam El Hantati; Ahmed Adri; Hatim Fakhreddine; Said Rifai; Rhali Benamar
Journal: Shock and Vibration
 (2021)

---

## Abstract

The scope of this study is to present a contribution to the geometrically nonlinear free and forced vibration of multiple-stepped beams, based on the theories of Euler–Bernoulli and von Karman, in order to calculate their corresponding amplitude-dependent modes and frequencies. Discrete expressions of the strain energy and kinetic energies are derived, and Hamilton’s principle is applied to reduce the problem to a solution of a nonlinear algebraic system and then solved by an approximate method. The forced vibration is then studied based on a multimode approach. The effect of nonlinearity on the dynamic behaviour of multistepped beams in the free and forced vibration is demonstrated and discussed. The effect of varying some geometrical parameters of the stepped beams in the free and forced cases is investigated and illustrated, among which is the variation in the level of excitation.


---

## Body

## 1. Introduction

In different engineering fields, such as the automotive, aeronautical, civil, or mechanical engineering, many structural components that permit the reinforcement of the whole system, a good resistance to the working loads, or the transmission of motions, can be modeled as stepped beams. In their operation, this type of beams is often subjected to a high dynamic load that needs to be investigated in order to be able and reliably predict the dynamic response of the beam. Many research works have examined the linear behaviour of this type of beams, making the analysis very easy, but leading generally to unreliable results when large vibration amplitudes are encountered. Since, strictly speaking, no system behaves always linearly, the linear results are often misleading. As the nonlinear effects arise in the case of the large displacements induced in many dynamic operating conditions, this led us to include the influence of geometrical nonlinearity in the present analysis.A literature survey on the transverse linear vibration of stepped beams goes backwards to Taleb and Suppiger [1], who presented some of the results of the integral equation theory, applied the Cauchy function method to achieve an approximate estimation of the fundamental frequency and modal configuration of the transverse vibration of beams with one step, and compared it with the exact solution. The study of the performance of a four-degree-of-freedom-per-node element model for the vibration analysis of uniform and multistepped beams was carried out by Balasubramanian and Subramanian [2]. Laura et al. [3] used the finite element method and compared their results with the experimental ones. Popplewell and Daqing Chang [4] used the Rayleigh–Ritz method to determine the unknown constants of the transverse deflection of a stepped beam. Krishnan et al. [5] discussed the use of the finite difference method in the study of stepped beams. Naguleswaran [6–9] studied the vibration of stepped beams in different configurations involving different boundary conditions, various numbers of steps, and types of cross sections. The natural frequencies were determined by iterative procedures based on linear interpolations. The vibration of a stepped beam with elastic supports at the ends was studied by De Rosa [10]. Sato [11] presented the analytical results relative to the effect of an abrupt change of cross section on the natural frequency of a beam with a groove, which were very well correlated with the experimental results. Sairigul and Aksu [12] used a method based on variational principles in conjunction with the finite difference technique to investigate the free vibration analysis of stepped Timoshenko beams. The Adomian decomposition method (ADM) was applied by Mao in [13, 14] in order to examine the free vibration of stepped and multistepped beams. A new method for calculating the frequencies and mode shapes was presented by Lin and Ng [15], who considered this method to be well adapted to forced-vibration applications.Although most work on stepped beams has focused on linear vibration analysis, few papers have included the effect of geometric nonlinearity in their studies. Using the transfer matrix method, Sato [16] exposed a study on the nonlinear free vibration of stepped beams. Nuttawit Wattanasakulpong and Arisara Chaikittiratana [17] studied the problems of linear and nonlinear free vibrations of stepped beams using the Adomian modified decomposition method. None of the articles mentioned above have investigated the geometrically nonlinear dynamic response of forced stepped beams.This study presents an extension of the works previously initiated by Benamar et al. [18–20], with the aim to contribute to the modal analysis of the geometrically nonlinear vibrations of straight thin structures, such as beams, plates, and circular cylindrical shells. Azrar et al. [21] presented an analysis of the geometrically nonlinear forced vibrations of uniform beams based on the multimode model. El Kadiri et al. [22] developed a simplified multimode approach (MMA), which is, according to the authors, easy to use for the study of the geometrically nonlinear forced vibrations of clamped-clamped and simply supported-clamped beams at both sides. El Kadiri and Benamar [23] extended this method to investigate the case of rectangular plates. Merrimi et al. [24] studied a cracked beam excited by a harmonic concentrated force, using both multimode and single-mode approaches. Fakhreddine et al. [25] investigated using the same method the geometrically nonlinear vibrations of beams resting on multiple supports and subjected to concentrated and uniformly distributed harmonic forces. Fakhreddine et al. [26] also extended this method to the case of nonlinear forced vibration of beams carrying concentric masses of different magnitudes and for vibration amplitudes.The main objective of this work was to present a contribution to the geometrically nonlinear free and forced vibration of beams with multiple steps based on the Euler–Bernoulli and von Karman theories. First, the natural frequencies are calculated by the Newton–Raphson algorithm leading to the determination of linear modes. The discretized expressions for the bending strain energy, the axial strain energy due to the nonlinear stretching forces, and the kinetic energy are then derived. The problem is reduced to a nonlinear algebraic system using Hamilton’s principle and then solved by the so-called second formulation reported in [22]. The numerical method convergence study has been carried out in the free and forced vibrations. The effects of the changes in cross section, inertia, and length ratios are analyzed and then illustrated. Considering a uniformly distributed excitation, the forced vibration case is studied on the basis of a multimode approach. The nonlinear frequency response functions are calculated and presented in the neighborhood of the predominant mode. The effects of the variation in the excitation level, the type of cross section, and the length ratio of stepped beams are investigated and illustrated. The presented method can also be extended for other types of beams including nano-/microbeams, previously investigated in [27] using the modified strain gradient and hyperbolic shear deformation beam theories and in [28] using the Euler–Bernoulli beam theory via the enhanced Eringen differential model. It can also be extended for the case of free vibration nanotubes, earlier studied in [29], using the discrete singular convolution technique.


## 2. General Formulation

### 2.1. Linear Formulation

This study concerns nonuniform homogeneous beams of the Euler–Bernoulli type, corresponding to the three configurations shown in Figure1, characterized by n different steps at various points. The first case is a beam with a rectangular cross-section variation in the width (a), the second in the height (b), and the third in both (c). In all cases, the beams are vibrating transversely.Figure 1
Three representative cases of stepped beams.Figure2 shows that each of the multistepped beams considered is divided into n different steps of uniform cross section and is supported by linear and rotational springs at both ends.According to Taleb and Suppiger [1], the transverse free vibration of any segment i of the beam with a uniform cross section is governed by the following equation:(1)d4wxdx4−β4wx=0;xϵ0Ljforj=1,2,…,N+1,where(2)βij4=ωi2ρjSjEjIj,Sj=ejHj,Ij=ejHj312,fori=1,2,…,10.Each step of the stepped beams has a specified length, thickness, width, Young’s modulus, area of cross section, moment of inertia, mass per unit length, and eigenvalue parameter (Lj,Hj,ej,Ej,, Sj,Ij,ρj, and βij). The variation of each step is presented by the ratios of the cross sections, moments of inertia, step, length, and eigenvalue parameters (uj,vj,λj,Rj, and zj).(3)xj∗=xjL,Rj∗=LjL,zij=βijβi1,zij∗=βijβij−1,uj=SjS1,vj=IjI1,vj∗=IjIj−1.The ‘active’ dimension denoted bydj of beams (a), (b), and (c) is, respectively, the width, the height, and both of them, which allow one to write the ratio of steps and moments of inertia as follows.For beam (a),(4)λj=djd1,vj=djd1.For beam (b),(5)λj=djd1,vj=djd13.For beam (c),(6)λj=djd12,vj=djd14.The general solution of equation (1) leads to the linear mode shapes, which will be used as basic functions in the nonlinear theory detailed in the following. The transverse displacement function of the multistepped beam, shown in Figure 2, can be expressed in each span as given by [25](7)wix∗=wijx∗⟶0,x1∗wijx∗⟶xj∗,xj+1∗……wijx∗⟶xN∗,1,where xj−1∗≤x∗≤xj∗, and also,(8)wijx∗=AjcoshβijLx∗−xj−1∗+BjsinhβijLx∗−xj−1∗+CjcosβijLx∗−xj−1∗+DjsinβijLx∗−xj−1∗.Figure 2
The coordinate system for a multistepped beam.By varying indicesi from (1) to (7) and j from 1 to N + 1, the satisfaction of the boundary conditions together with the compatibility conditions and the transfer matrix leads to a homogeneous system. The boundary conditions for a beam segment supported by linear and rotational springs at both sides are given in [25].At the left side,(9)d3wi1x∗dx∗3|x∗=0=−kt1wi1x∗|x∗=0,d2wi1x∗dx∗2|x∗=0=kθ1dwi1x∗dx∗|x∗=0.At the right side,(10)d3wiN+1x∗dx∗3|x∗=1=kt2wiN+1x∗|x∗=1,d2wiN+1x∗dx∗2|x∗=1=−kθ1dwiN+1x∗dx∗|x∗=1.The stiffness values of the transverse and rotary springs at both ends of the beam are denoted askt1,kt2,kθ1,kθ2. The compatibility conditions for each beam step are expressed as shown in [13]:(11)wijx∗|x∗=xj∗=wij+1x∗|x∗=xj∗,dwijx∗dx∗|x∗=xj∗=dwij+1x∗dx∗|x∗=xj∗,EIjd2wijx∗dx∗2|x∗=xj∗=EIj+1d2wij+1x∗dx∗2|x∗=xj∗,EIjd3wijx∗dx∗3|x∗=xj∗=EIj+1d3wij+1x∗dx∗3|x∗=xj∗.


### 2.2. The Transfer Matrix Method

Passing through the compatibility conditions, it is found that the constants (Aj,Bj,Cj, and Dj) in the jth step depend on those of the j−1th step (Aj−1,Bj−1,Cj−1, and Dj−1), which can be presented in a matrix form as in [30]:(12)Hjηj=Gj−1ηj−1,with(13)ηj=AjBjCjDj;ηj−1=Aj−1Bj−1Cj−1Dj−1,Gj−1=coshβij−1Lj−1sinhβij−1Lj−1cosβij−1Lj−1sinβij−1Lj−1sinhβij−1Lj−1coshβij−1Lj−1−sinβij−1Lj−1cosβij−1Lj−1coshβij−1Lj−1sinhβij−1Lj−1−cosβij−1Lj−1−sinβij−1Lj−1sinhβij−1Lj−1coshβij−1Lj−1sinβij−1Lj−1−cosβij−1Lj−1,Hj=10100zij∗0zij∗zij∗2vj∗0−zij∗2vj∗00zij∗3vj∗0−zij∗3vj∗.From equation (12), one obtains(14)ηj=Hj−1Gj−1ηj−1,where Hj−1;Gj−1 are 4×4 matrices, which form the transfer matrix Q expressed by(15)Qj−1=Hj−1Gj−1.Equations (17) and (18) lead to(16)ηj=Qj−1ηj−1.Under successive repetition of equation (12), the four constants of the first segment (A1,B1,C1, and D1) appear to be related with those of the final segment (AN+1,BN+1,CN+1, and DN+1), which decreases the number of independent constants to four.(17)AN+1BN+1CN+1DN+1=QN×⋯×Q2×Q1A1B1C1D1.The four constants (Aj,Bj,Cj, and Dj) are determined when the boundary conditions are fully satisfied and expressed as given by Chajdi et al. [31], but in this case, the eigenvalue parameter βij differs from one span to another. The boundary conditions at the right side lead to the following expressions:(18)AN+1βiN+13sinhβiN+1LN+1−kt2coshβiN+1LN+1+BN+1βiN+13coshβiN+1LN+1−kt2sinhβiN+1LN+1+CN+1βiN+13sinβiN+1LN+1−kt2cosβiN+1LN+1−DN+1βiN+13cosβiN+1LN+1+kt2sinβiN+1LN+1=0,(19)AN+1βiN+1coshβiN+1LN+1+kθ2sinhβiN+1LN+1+BN+1βiN+1sinhβiN+1LN+1+kθ2coshβiN+1LN+1−CN+1βiN+1cosβiN+1LN+1+kθ2sinβiN+1LN+1−DN+1βiN+1sinβiN+1LN+1−kθ2cosβiN+1LN+1=0.Equations (18) and (19) can be expressed in a matrix form as follows:(20)D×AN+1BN+1CN+1DN+1=D11D21D12D22D13D23D14D24T×AN+1BN+1CN+1DN+1=00such as(21)D11=βiN+13sinhβiN+1LN+1−kt2coshβiN+1LN+1,D12=βiN+13coshβiN+1LN+1−kt2sinhβiN+1LN+1,D13=βiN+13sinβiN+1LN+1−kt2cosβiN+1LN+1,D14=−βiN+13cosβiN+1LN+1−kt2sinβiN+1LN+1,D21=βiN+1coshβiN+1LN+1+kθ2sinhβiN+1LN+1,D22=βiN+1sinhβiN+1LN+1+kθ2coshβiN+1LN+1,D23=−βiN+1cosβiN+1LN+1−kθ2sinβiN+1LN+1,D24=−βiN+1sinβiN+1LN+1+kθ2cosβiN+1LN+1.Substitution of (17) into equation (20) leads to(22)D×AN+1BN+1CN+1DN+1=D×QN×⋯×Q2×Q1A1B1C1D1=00.The boundary conditions at the left side can be formulated in a matrix form as given in [31]:(23)kt1βi3kt1−βi3−βikθ1βikθ1T×A1B1C1D1=R×A1B1C1D1=00.Equations (22) and (23) lead to(24)D×QN×⋯×Q2×Q1R×A1B1C1D1=Y×A1B1C1D1=0000,where(25)Y=D×QN×⋯×Q2×Q1R.The homogeneous system is solved by the Newton–Raphson algorithm, where the determinant of matrixY is assigned to be equal to zero in order to allow nontrivial solutions to exist:(26)detY=Y=0.After determining the natural frequencies from equation (26), the constants (Aj, Bj, Cj and Dj) are calculated by the boundary and continuity conditions written above.


## 2.1. Linear Formulation

This study concerns nonuniform homogeneous beams of the Euler–Bernoulli type, corresponding to the three configurations shown in Figure1, characterized by n different steps at various points. The first case is a beam with a rectangular cross-section variation in the width (a), the second in the height (b), and the third in both (c). In all cases, the beams are vibrating transversely.Figure 1
Three representative cases of stepped beams.Figure2 shows that each of the multistepped beams considered is divided into n different steps of uniform cross section and is supported by linear and rotational springs at both ends.According to Taleb and Suppiger [1], the transverse free vibration of any segment i of the beam with a uniform cross section is governed by the following equation:(1)d4wxdx4−β4wx=0;xϵ0Ljforj=1,2,…,N+1,where(2)βij4=ωi2ρjSjEjIj,Sj=ejHj,Ij=ejHj312,fori=1,2,…,10.Each step of the stepped beams has a specified length, thickness, width, Young’s modulus, area of cross section, moment of inertia, mass per unit length, and eigenvalue parameter (Lj,Hj,ej,Ej,, Sj,Ij,ρj, and βij). The variation of each step is presented by the ratios of the cross sections, moments of inertia, step, length, and eigenvalue parameters (uj,vj,λj,Rj, and zj).(3)xj∗=xjL,Rj∗=LjL,zij=βijβi1,zij∗=βijβij−1,uj=SjS1,vj=IjI1,vj∗=IjIj−1.The ‘active’ dimension denoted bydj of beams (a), (b), and (c) is, respectively, the width, the height, and both of them, which allow one to write the ratio of steps and moments of inertia as follows.For beam (a),(4)λj=djd1,vj=djd1.For beam (b),(5)λj=djd1,vj=djd13.For beam (c),(6)λj=djd12,vj=djd14.The general solution of equation (1) leads to the linear mode shapes, which will be used as basic functions in the nonlinear theory detailed in the following. The transverse displacement function of the multistepped beam, shown in Figure 2, can be expressed in each span as given by [25](7)wix∗=wijx∗⟶0,x1∗wijx∗⟶xj∗,xj+1∗……wijx∗⟶xN∗,1,where xj−1∗≤x∗≤xj∗, and also,(8)wijx∗=AjcoshβijLx∗−xj−1∗+BjsinhβijLx∗−xj−1∗+CjcosβijLx∗−xj−1∗+DjsinβijLx∗−xj−1∗.Figure 2
The coordinate system for a multistepped beam.By varying indicesi from (1) to (7) and j from 1 to N + 1, the satisfaction of the boundary conditions together with the compatibility conditions and the transfer matrix leads to a homogeneous system. The boundary conditions for a beam segment supported by linear and rotational springs at both sides are given in [25].At the left side,(9)d3wi1x∗dx∗3|x∗=0=−kt1wi1x∗|x∗=0,d2wi1x∗dx∗2|x∗=0=kθ1dwi1x∗dx∗|x∗=0.At the right side,(10)d3wiN+1x∗dx∗3|x∗=1=kt2wiN+1x∗|x∗=1,d2wiN+1x∗dx∗2|x∗=1=−kθ1dwiN+1x∗dx∗|x∗=1.The stiffness values of the transverse and rotary springs at both ends of the beam are denoted askt1,kt2,kθ1,kθ2. The compatibility conditions for each beam step are expressed as shown in [13]:(11)wijx∗|x∗=xj∗=wij+1x∗|x∗=xj∗,dwijx∗dx∗|x∗=xj∗=dwij+1x∗dx∗|x∗=xj∗,EIjd2wijx∗dx∗2|x∗=xj∗=EIj+1d2wij+1x∗dx∗2|x∗=xj∗,EIjd3wijx∗dx∗3|x∗=xj∗=EIj+1d3wij+1x∗dx∗3|x∗=xj∗.


## 2.2. The Transfer Matrix Method

Passing through the compatibility conditions, it is found that the constants (Aj,Bj,Cj, and Dj) in the jth step depend on those of the j−1th step (Aj−1,Bj−1,Cj−1, and Dj−1), which can be presented in a matrix form as in [30]:(12)Hjηj=Gj−1ηj−1,with(13)ηj=AjBjCjDj;ηj−1=Aj−1Bj−1Cj−1Dj−1,Gj−1=coshβij−1Lj−1sinhβij−1Lj−1cosβij−1Lj−1sinβij−1Lj−1sinhβij−1Lj−1coshβij−1Lj−1−sinβij−1Lj−1cosβij−1Lj−1coshβij−1Lj−1sinhβij−1Lj−1−cosβij−1Lj−1−sinβij−1Lj−1sinhβij−1Lj−1coshβij−1Lj−1sinβij−1Lj−1−cosβij−1Lj−1,Hj=10100zij∗0zij∗zij∗2vj∗0−zij∗2vj∗00zij∗3vj∗0−zij∗3vj∗.From equation (12), one obtains(14)ηj=Hj−1Gj−1ηj−1,where Hj−1;Gj−1 are 4×4 matrices, which form the transfer matrix Q expressed by(15)Qj−1=Hj−1Gj−1.Equations (17) and (18) lead to(16)ηj=Qj−1ηj−1.Under successive repetition of equation (12), the four constants of the first segment (A1,B1,C1, and D1) appear to be related with those of the final segment (AN+1,BN+1,CN+1, and DN+1), which decreases the number of independent constants to four.(17)AN+1BN+1CN+1DN+1=QN×⋯×Q2×Q1A1B1C1D1.The four constants (Aj,Bj,Cj, and Dj) are determined when the boundary conditions are fully satisfied and expressed as given by Chajdi et al. [31], but in this case, the eigenvalue parameter βij differs from one span to another. The boundary conditions at the right side lead to the following expressions:(18)AN+1βiN+13sinhβiN+1LN+1−kt2coshβiN+1LN+1+BN+1βiN+13coshβiN+1LN+1−kt2sinhβiN+1LN+1+CN+1βiN+13sinβiN+1LN+1−kt2cosβiN+1LN+1−DN+1βiN+13cosβiN+1LN+1+kt2sinβiN+1LN+1=0,(19)AN+1βiN+1coshβiN+1LN+1+kθ2sinhβiN+1LN+1+BN+1βiN+1sinhβiN+1LN+1+kθ2coshβiN+1LN+1−CN+1βiN+1cosβiN+1LN+1+kθ2sinβiN+1LN+1−DN+1βiN+1sinβiN+1LN+1−kθ2cosβiN+1LN+1=0.Equations (18) and (19) can be expressed in a matrix form as follows:(20)D×AN+1BN+1CN+1DN+1=D11D21D12D22D13D23D14D24T×AN+1BN+1CN+1DN+1=00such as(21)D11=βiN+13sinhβiN+1LN+1−kt2coshβiN+1LN+1,D12=βiN+13coshβiN+1LN+1−kt2sinhβiN+1LN+1,D13=βiN+13sinβiN+1LN+1−kt2cosβiN+1LN+1,D14=−βiN+13cosβiN+1LN+1−kt2sinβiN+1LN+1,D21=βiN+1coshβiN+1LN+1+kθ2sinhβiN+1LN+1,D22=βiN+1sinhβiN+1LN+1+kθ2coshβiN+1LN+1,D23=−βiN+1cosβiN+1LN+1−kθ2sinβiN+1LN+1,D24=−βiN+1sinβiN+1LN+1+kθ2cosβiN+1LN+1.Substitution of (17) into equation (20) leads to(22)D×AN+1BN+1CN+1DN+1=D×QN×⋯×Q2×Q1A1B1C1D1=00.The boundary conditions at the left side can be formulated in a matrix form as given in [31]:(23)kt1βi3kt1−βi3−βikθ1βikθ1T×A1B1C1D1=R×A1B1C1D1=00.Equations (22) and (23) lead to(24)D×QN×⋯×Q2×Q1R×A1B1C1D1=Y×A1B1C1D1=0000,where(25)Y=D×QN×⋯×Q2×Q1R.The homogeneous system is solved by the Newton–Raphson algorithm, where the determinant of matrixY is assigned to be equal to zero in order to allow nontrivial solutions to exist:(26)detY=Y=0.After determining the natural frequencies from equation (26), the constants (Aj, Bj, Cj and Dj) are calculated by the boundary and continuity conditions written above.


## 3. Nonlinear Formulation

### 3.1. Free Vibration

The dynamic behaviour of the structure is determined by Hamilton’s principle which is formally expressed as follows:(27)δ∫02ω/tV−Tdt=0.T is the beam kinetic energy. V is the total strain energy of the multistep beam, expressed as the sum of the energy due to the bending Vf and the energy Va due to the nonlinear stretching forces induced by the large deflections.SettingW·=∂W/∂t, the kinetic energy can be expressed as in [32]:(28)T=12∫0LρSxW•2dx.SettingW″=d2W/dx2, the deformation energy due to bending can be expressed as in [32]:(29)Vf=12∫0LEIxW″2dx.Equations (28) and (29) can be developed as the following form:(30)T=ρ2∫0x1S1W•12dx+∫x2x1S2W•22dx+⋯+∫xNLSN+1W•N+12dx,(31)Vf=E2∫0x1I1W1″2dx+∫x1x2I2W2″2dx+⋯+∫xNLIN+1WN+1″2dx.The multistep beam has a variable cross section that leads to a variable moment of inertia, which changes from step to step. The axial deformation energy due to the nonlinear stretching forces is calculated as follows:(32)Va=12∫0LεsNxdx.The normal loadNx induced by the nonlinear stretching force of the beam is assumed, when the axial inertia forces are negligible, to be constant along the whole beam span as in [33].SettingW′=dW/dx, the nonlinear stretching force can be given as(33)Nx=E2.∫0LW′2dx∫0Ldx/Sx,where(34)∫0LdxSx=LS1C,C=R1u1+R2u2+⋯+RN+1uN+1.The nonlinear axial strainεS is expressed as(35)εs=NxES.Substituting equations (33) and (35) in equation (32) gives the following equation:(36)Va=ES18L.∫0LW′2dx2C.By assuming a harmonic motion and expanding the transverse displacement as a finite series of basic spatial functions, it is possible to write(37)Wx,t=aiwixsinωt,where wix and ai are the multistepped beam ith linear mode and basic function contribution coefficient, respectively. The replacement of the new w form in equations (30), (31), and (36) and the discretization lead to the following expressions for the kinetic energy, the axial strain energy due to nonlinear stretching forces, and the strain energy due to bending as in [22]:(38)T=12ω2aiaicosωt2mij,(39)Va=12aiajakalsinωt4bijkl,(40)Vf=12aiaisinωt2kij,in which mij, kij, and bijkl are the mass matrix due to T, the rigidity matrix due to Vf, and the nonlinearity tensor due to Va, expressed by(41)mij=ρS12∫0x1u1wi1wj1dx+⋯+∫xNxN+1uN+1wiN+1wjN+1dx,kij=EI12∫0x1v1wi1″wj1″dx+⋯+∫xNxN+1vN+1wiN+1″wjN+1″dx,bijkl=ES14.C∫0Lwi′wj′dx∫0Lwk′wl′dx.Using equations (38)–(40) and applying Hamilton’s principle lead to(42)2aikir+3aiajakbijkr−2ω2aimir=0,r=1,…,n.According to El Kadiri et al. [22], the above equation can be expressed in a matrix form as follows:(43)2KA+3BAA−2ω2MA=0,where kir,mir,ajakbijkl, and ai are the general terms of K,M,BA, and A. The column vector of the contribution coefficients A and the frequency ω are the unknowns. Before they are determined, equation (43) is put in a dimensionless form by replacing the dimensional parameters by the corresponding nondimensional ones:(44)wix∗wi∗x∗=H1,(45)ω2ω∗2=EI1ρS1L4,(46)MijMij∗=ρS1H12L,(47)Kij=EI1H12L3Kij∗,(48)Bijkl=EI1H12L3Bijkl∗.Replacing equations (44)–(48) in equation (42) or (43) results in the following equation:(49)K∗A+32B∗AA−ω∗2M∗A=0,which can also be written as(50)2aikir∗+3aiajakbijkr∗−2ω∗2aimir∗=0.The determination of the unknowns is based on the so-called second formulation developed previously in [22], where it was shown that the contribution of the fundamental mode (a1) remains predominant for the full range of the vibration amplitudes considered. The other contributions, which are small in comparison with (a1), are denoted by ε2,…,εn. Consequently, following [22], one can write aiajakbijkr∗=a13b111r+a12εib11ir, which allows equation (50) to become, in the vicinity of the rth mode,(51)Kr∗R−ω∗2Mr∗RArR+32αr∗ArR=−32ar3biiir∗,where(52)αr∗=ar3bijrr∗forr=1,…,n.This approximate linear system, given in [22], makes it very easy to calculate the contribution coefficients.


### 3.2. Forced Vibration

The study of nonlinear forced vibrations is made in the three configurations, which differ by the type of the cross section. In all cases, the beam is subjected to a distributed force, which excites the modes of the structure through a series of generalized forces. The distributed harmonic force is denoted byFid and depends on the expression for Fd, the portion of the beam in which the excitation is applied, and the mode considered. The expressions for the distributed harmonic force are extracted from [20] and written as follows:(53)Fid=Fdsinωt∫0Lwix0dx=fidsinωt.To investigate the dynamic behaviour of the conservative system, Hamilton principle, in which the forcing term has been introduced, has been applied as in [26]:(54)δ∫02ω/tV−T+WFdt=0.The kinetic energy, the total strain energy, and the work done by the external loads of the multistep beam are denoted byT, V, and WF. After calculations, one gets(55)2KA+3BAA−2ω2MA=F.The dimensionless generalized distributed forces are denoted byfid∗:(56)fid∗=L4EI1H1Fd∫01wi∗x∗dx∗.Equation (57) is made dimensionless by using the dimensionless parameters in equation (55) to get(57)K∗−ω∗2M∗A+32B∗AA=F∗.According to El Kadiri et al. [22], an approximate linear system has been obtained as in the nonlinear free vibration case developed above, which permits equation (57) to be expressed as follows:(58)Kr∗R−ω∗2Mr∗RArR+32αr∗ArR=F∗−32ar3biiir∗.The necessary use of the multimode approach in the neighborhood of the predominant mode makes it very easy to calculate the contribution coefficients, as in [22].


## 3.1. Free Vibration

The dynamic behaviour of the structure is determined by Hamilton’s principle which is formally expressed as follows:(27)δ∫02ω/tV−Tdt=0.T is the beam kinetic energy. V is the total strain energy of the multistep beam, expressed as the sum of the energy due to the bending Vf and the energy Va due to the nonlinear stretching forces induced by the large deflections.SettingW·=∂W/∂t, the kinetic energy can be expressed as in [32]:(28)T=12∫0LρSxW•2dx.SettingW″=d2W/dx2, the deformation energy due to bending can be expressed as in [32]:(29)Vf=12∫0LEIxW″2dx.Equations (28) and (29) can be developed as the following form:(30)T=ρ2∫0x1S1W•12dx+∫x2x1S2W•22dx+⋯+∫xNLSN+1W•N+12dx,(31)Vf=E2∫0x1I1W1″2dx+∫x1x2I2W2″2dx+⋯+∫xNLIN+1WN+1″2dx.The multistep beam has a variable cross section that leads to a variable moment of inertia, which changes from step to step. The axial deformation energy due to the nonlinear stretching forces is calculated as follows:(32)Va=12∫0LεsNxdx.The normal loadNx induced by the nonlinear stretching force of the beam is assumed, when the axial inertia forces are negligible, to be constant along the whole beam span as in [33].SettingW′=dW/dx, the nonlinear stretching force can be given as(33)Nx=E2.∫0LW′2dx∫0Ldx/Sx,where(34)∫0LdxSx=LS1C,C=R1u1+R2u2+⋯+RN+1uN+1.The nonlinear axial strainεS is expressed as(35)εs=NxES.Substituting equations (33) and (35) in equation (32) gives the following equation:(36)Va=ES18L.∫0LW′2dx2C.By assuming a harmonic motion and expanding the transverse displacement as a finite series of basic spatial functions, it is possible to write(37)Wx,t=aiwixsinωt,where wix and ai are the multistepped beam ith linear mode and basic function contribution coefficient, respectively. The replacement of the new w form in equations (30), (31), and (36) and the discretization lead to the following expressions for the kinetic energy, the axial strain energy due to nonlinear stretching forces, and the strain energy due to bending as in [22]:(38)T=12ω2aiaicosωt2mij,(39)Va=12aiajakalsinωt4bijkl,(40)Vf=12aiaisinωt2kij,in which mij, kij, and bijkl are the mass matrix due to T, the rigidity matrix due to Vf, and the nonlinearity tensor due to Va, expressed by(41)mij=ρS12∫0x1u1wi1wj1dx+⋯+∫xNxN+1uN+1wiN+1wjN+1dx,kij=EI12∫0x1v1wi1″wj1″dx+⋯+∫xNxN+1vN+1wiN+1″wjN+1″dx,bijkl=ES14.C∫0Lwi′wj′dx∫0Lwk′wl′dx.Using equations (38)–(40) and applying Hamilton’s principle lead to(42)2aikir+3aiajakbijkr−2ω2aimir=0,r=1,…,n.According to El Kadiri et al. [22], the above equation can be expressed in a matrix form as follows:(43)2KA+3BAA−2ω2MA=0,where kir,mir,ajakbijkl, and ai are the general terms of K,M,BA, and A. The column vector of the contribution coefficients A and the frequency ω are the unknowns. Before they are determined, equation (43) is put in a dimensionless form by replacing the dimensional parameters by the corresponding nondimensional ones:(44)wix∗wi∗x∗=H1,(45)ω2ω∗2=EI1ρS1L4,(46)MijMij∗=ρS1H12L,(47)Kij=EI1H12L3Kij∗,(48)Bijkl=EI1H12L3Bijkl∗.Replacing equations (44)–(48) in equation (42) or (43) results in the following equation:(49)K∗A+32B∗AA−ω∗2M∗A=0,which can also be written as(50)2aikir∗+3aiajakbijkr∗−2ω∗2aimir∗=0.The determination of the unknowns is based on the so-called second formulation developed previously in [22], where it was shown that the contribution of the fundamental mode (a1) remains predominant for the full range of the vibration amplitudes considered. The other contributions, which are small in comparison with (a1), are denoted by ε2,…,εn. Consequently, following [22], one can write aiajakbijkr∗=a13b111r+a12εib11ir, which allows equation (50) to become, in the vicinity of the rth mode,(51)Kr∗R−ω∗2Mr∗RArR+32αr∗ArR=−32ar3biiir∗,where(52)αr∗=ar3bijrr∗forr=1,…,n.This approximate linear system, given in [22], makes it very easy to calculate the contribution coefficients.


## 3.2. Forced Vibration

The study of nonlinear forced vibrations is made in the three configurations, which differ by the type of the cross section. In all cases, the beam is subjected to a distributed force, which excites the modes of the structure through a series of generalized forces. The distributed harmonic force is denoted byFid and depends on the expression for Fd, the portion of the beam in which the excitation is applied, and the mode considered. The expressions for the distributed harmonic force are extracted from [20] and written as follows:(53)Fid=Fdsinωt∫0Lwix0dx=fidsinωt.To investigate the dynamic behaviour of the conservative system, Hamilton principle, in which the forcing term has been introduced, has been applied as in [26]:(54)δ∫02ω/tV−T+WFdt=0.The kinetic energy, the total strain energy, and the work done by the external loads of the multistep beam are denoted byT, V, and WF. After calculations, one gets(55)2KA+3BAA−2ω2MA=F.The dimensionless generalized distributed forces are denoted byfid∗:(56)fid∗=L4EI1H1Fd∫01wi∗x∗dx∗.Equation (57) is made dimensionless by using the dimensionless parameters in equation (55) to get(57)K∗−ω∗2M∗A+32B∗AA=F∗.According to El Kadiri et al. [22], an approximate linear system has been obtained as in the nonlinear free vibration case developed above, which permits equation (57) to be expressed as follows:(58)Kr∗R−ω∗2Mr∗RArR+32αr∗ArR=F∗−32ar3biiir∗.The necessary use of the multimode approach in the neighborhood of the predominant mode makes it very easy to calculate the contribution coefficients, as in [22].


## 4. Convergence Study

The purpose of this section is to investigate the convergence of the series expansion used to solve the nonlinear vibration problem via the amplitude equation obtained in the numerical approach exposed above. These analyses deal with a two-step fully clamped beam of type (b), whereR1∗=R3∗=0.5×R2∗, u1=u3, and u2=0.6×u1, subjected to a free and forced vibration.


### 4.1. Free Vibration

In order to solve equation (51), the assumptions stated by El Kadiri et al. were used to simplify the nonlinear algebraic system obtained. The contribution coefficients have been calculated as well as the associated nonlinear frequency. The normalized first nonlinear mode Wnl∗x,a1, for a given value of a1, is obtained as a series of basic functions involving the beam modal parameters depending on w1∗,…,w10∗ and given as follows:(59)Wnl∗x,a1=a1w1∗+⋯+ε10w10∗.The predominant term, relative to the first linear mode shape, isa1w1∗, while the other terms proportional to the symmetric and antisymmetric higher modes ε2w2∗,ε3w3∗,…,ε10w10∗ present the corrections due to the geometrical nonlinearity.Table1 shows the result of contribution coefficients obtained using equation (51); these contribution coefficients are included in equation (59) in order to calculate the normalized first nonlinear mode Wnl∗x,a1. Its maximal value is denoted by Wmax and is located at the middle of the beam. From Table 1, some observations can be clearly made:(i)
The contribution of the fundamental modea1 remains predominant for the full range of the vibration amplitudes considered. The other contributions denoted by ε2,…,εn remain small with respect to a1.(ii)
The contribution coefficients of the modes decrease asi indices increase.(iii)
The contribution coefficients increased according to the maximum nondimensional amplitudeWmax∗.Table 1
The linear and nonlinear frequency ratio and the contribution coefficients of stepped beam (b) subjected to a free vibration.
Wmaxa1ωnl∗/ωl∗ε2ε3ε4ε5ε6ε7ε8ε9ε100.9680.51.233.43E−03−4.21E031.42E−04−4.52E−04−2.912E−05−1.04E−04−2.77E−05−2.96E−051.39E−051.1660.61.325.64E−03−7.01E−032.69E−04−7.74E−04−4.449E−05−1.79E−04−4.61E−05−5.07E−052.36E−051.3650.71.428.48E−2−1.06E−024.65E−03−1.21E−03−6.058E−05−2.83E−04−7.03E−05−7.98E−053.66E−051.5660.81.531.19E−21.51E−027.48E−04−1.78E−03−7.46E−05−4.19E−03−1.00E−041.17E−045.32E−05From these observations, it is quite interesting to verify the convergence of the nonlinear beam deflection equation, which is obtained as a series of basic functions. The increase of these basic functions added in equation (59) enhanced the correction due to geometrical nonlinearity. The effect of varying the number of basic function series on the convergence of the numerical method has been investigated by comparing the results of a nonlinear frequency estimated ωeN∗. This latter can be expressed as in [26]:(60)ωeN∗=ATK∗A+2/3ATB∗AAATM∗A,forA=a1,ε2w2∗,…,εNwN∗.The nonlinear frequency estimatedωeN∗ can be extracted from equation (60), where A presents the contribution coefficient vector already obtained from the second formulation, using N basic functions. For a1=0.5, the nonlinear frequency ωeN∗ is calculated for different cases, in which the number of basic functions N increased from 7 to 10. The same test is repeated for different values of a1, varying from 0.5 to 0.8. The differences between ωeN∗ and ωeN-1∗ are calculated, as well as the residual of each test, and summarized hereafter in Table 2.Table 2
The frequency differences and residual for each value ofa1 and at different numbers of basic functions.
Basic function numbersTesta1ωnl∗/ωl∗78910Residual10.51.231.45E−051.80E−072.90E−071.00E−081.452E−0520.61.321.08E−044.50E−067.00E−064.00E−071.079E−0430.71.422.30E−049.60E−061.50E−059.00E−072.305E−0440.81.534.42E−041.90E−052.90E−051.00E−064.433−04(61)DiffN=ωN−1∗−ωN∗,residual=∑i=710DiffN2.The numerical results presented in Table2 show that the differences DiffN for the four tests decrease by increasing the number of basic functions, while the residual increases by increasing the a1 value. The difference Diff7 of the first test is equal to 1.45×10−05, and after increasing the number of basic functions to 10, this difference decreases to a low value of 1×10−08, meaning that(i)
The use of 6 basic functions induces an error of1.5×10−05 with respect to ωe10∗(ii)
The use of 7 basic functions induces an error of4.8×10−07 with respect to ωe10∗(iii)
The use of 8 basic functions induces an error of3×10−07 with respect to ωe10∗(iv)
The use of 9 basic functions induces an error of1×10−08 with respect to ωe10∗Also, ata1 = 0.8, the difference for 7 basic functions increases from 1.45×10−05 to 4.42×10−04, and by adding 3 other terms, the difference decreased to Diff10=1×10−06. In other terms,(i)
The use of 6 basic functions induces an error of4.91×10−04 with respect to ωe10∗(ii)
The use of 7 basic functions induces an error of4.90×10−05 with respect to ωe10∗(iii)
The use of 8 basic functions induces an error of3.00×10−05 with respect to ωe10∗(iv)
The use of 9 basic functions induces an error of1×10−06 with respect to ωe10∗When the basic functions’ number increases and exceeds ten, the difference between the final iteration ofωeN∗ and its precedent decreases by value less than 1×10−08 for the first test and less than 1×10−06 for the last one, whereas the residual will be increased by less than these values. These corrections are almost negligible compared to ωe∗. These results are considered to be tightly convergent in the nonlinear free vibration case using ten or more basic functions.


### 4.2. Forced Vibration

These analyses deal with a two-step fully clamped beam of type (b), whereR1∗=R3∗=0.5×R2∗, u1=u3, and u2=0.6×u1), subjected to a uniformly distributed force over an interval of [0.4, 0.6]. Firstly, the contribution coefficients are obtained using multimode approach (58) and then summarized in Table 3.Table 3
The linear and nonlinear frequency ratio and the contribution coefficients of stepped beam (b) subjected to a forced vibration.
Wmaxωnl∗/ωl∗a1ε2ε3ε4ε5ε6ε7ε8ε9ε100.0960.2990.05173.986E−185.258E−043.117E−194.381E−052.534E−191.175E−051.781E−184.455E−062.281E−180.7481.0990.39981.049E−153.6292E−3−1.027E−16−3.855E−044.186E−17−5.581E−05−7.892E−17−1.487E−054.443E−171.131.30.61313.271E−151.092E−2−3.642e−16−1.479E−31.369E−16−2.332E−04−2.810E−16−6.657E−051.459E−16Afterwards, these contribution coefficients are included in equation (59) in order to calculate the normalized first nonlinear mode Wnl∗x,a1. Its maximal value is denoted by Wmax and is located at the middle of the beam. The same observations can be concluded from Table 3. Furthermore, it can be clearly seen that, in the forced case, the contribution coefficients of the even modes have very low values. However, a test has been carried out to predict the effect of neglecting the antisymmetric terms in equation (59) on the normalized first nonlinear mode. Two first normalized nonlinear modes have been calculated: the first one is obtained as shown in equation (59), and the second is also obtained as a series of basic functions involving only the symmetrical terms ε3w3∗,…,ε9w9∗ and denoted by WnlSymmetric∗. The results are summarized in Table 4.Table 4
Comparison between normalized first nonlinear modesWnl∗ and WnlSymmetric∗ for ωnl∗/ωl∗=1.8.
W∗∫01W∗d2W∗/dx∗2|x∗=0.5∫01d2W∗/dx∗2dxWnl∗1.953070559600.99504143259212831.92929395580.005546040320137WnlSymmetric∗1.953070559600.995041432592128431.92929395580.005546040320142Difference—4.4408E−16—5.381E−15The numerical results presented in Table4 prove that, by neglecting the antisymmetric modes, the use of equation (59) leads to the same results. The normalized first nonlinear mode can also be written as(62)Wnl∗x,a1=a1w1∗+ε3w3∗+ε5w5∗+ε7w7∗+ε9w9∗.The effect of varying the basic functions’ number on the convergence of equation (62) has been investigated by comparing the normalized first nonlinear mode obtained for five basic functions Wnl∗ for three scenarios.Scenario 1.
The nonlinear beam deflection assumes the use of one basic function:(63)Wnl3∗x,a1=a1w1∗+ε3w3∗.Scenario 2.
The nonlinear beam deflection assumes the use of two basic functions:(64)Wnl5∗x,a1=a1w1∗+ε3w3∗+ε5w5∗.Scenario 3.
The nonlinear beam deflection assumes the use of three basic functions:(65)Wnl7∗x,a1=a1w1∗+ε3w3∗+ε5w5∗+ε7w7∗.
The convergence betweenWnl∗ and Wnl3∗,Wnl5∗,andWnl7∗ has been examined by calculating the integral of these nonlinear modes over an interval of [0, 1] using the Simpson numerical integration method. The results are obtained for different values of ωnl∗/ωl∗ in the case where Fd is equal to 100, 200, and 300. The comparison was performed and is summarized in Table 5.
In Table5, the final normalized first nonlinear mode Wnl∗ (62) was compared with the three other scenarios of nonlinear modes as shown in equations (63)–(65). The comparison was performed by varying the frequency ratio for each excitation. The results obtained for the scenario where Fd=100 and ωnl∗/ωl∗=1 show a small difference of Diff3=1.724×10−5 for the first scenario, and it decreases to Diff7=4.5×10−7 for the last one, in other terms:(i)
The use ofScenario 1 induces an error of 1.98×10−05 with respect to ∫01Wnl∗(ii)
The use ofScenario 2 induces an error of 2.59×10−06 with respect to ∫01Wnl∗(iii)
The use ofScenario 3 induces an error of 4.5023×10−07 with respect to ∫01Wnl∗
Also, for the same excitation, ifωnl∗/ωl∗=1.8, Diff3 increases until reaching 1.11×10−3 for the first scenario, and Diff7 increases to 6.7609×10−5 for the last one, meaning that(i)
The use ofScenario 1 induces an error of 1.39×10−03 with respect to ∫01Wnl∗(ii)
The use ofScenario 2 induces an error of 2.78×10−04 with respect to ∫01Wnl∗(iii)
The use ofScenario 3 induces an error of 6.76×10−05 with respect to ∫01Wnl∗
It can be noticed that the change in the excitation level from 100 to 400 withωnl∗/ωl∗=1.8 increases the error involved between Wnl∗ and Wnl7∗ by a small variation. For Fd=100, this error is equal to 6.7609×10−5, and for Fd=400, this error increases to 7.148×10−5. Indeed, the effect of the frequency ratio variation is more remarkable than the force variation effect on the convergence of the normalized first nonlinear mode. The differences between Wnl∗ and Wnl7∗ for the value of the frequency ratio and for the three excitations considered remain small. However, the increase in the number of basic functions exceeding five symmetric terms in equation (62) provides a correction, due to the geometric nonlinearity, less than 7.148×10−5. This correction is almost negligible compared to Wnlmax∗; hence, the numerical method proposed is well converged.Table 5
The convergence test of the normalized first nonlinear modes with different numbers of basic functions in the forced vibration case.
FdWnlmax∗ωnl∗/ωl∗∫01Wnl∗∫01Wnl3∗∫01Wnl5∗∫01Wnl7∗Residual1000.53610.25807410.2580940.2580760.25807461.738E−05Diffi1.724E−052.135E−064.5023E−071.1351.30.5566270.556900.5566790.5566392.302E−04Diffi2.265E−043.934E−051.201E−061.461.50.7291680.729780.72928690.72919655.09E−04Diffi5.010E−049.033E−052.826E−051.91.80.967220.96861680.96750640.9672961.132E−03Diffi1.11E−032.10E−046.7609E−052000.67310.3255040.325544260.32550960.3255053.489E−05Diffi3.46E−054.321E−069.229E−071.1841.30.58206890.582378150.58212530.582082032.568E−04Diffi2527E−044.333E−051.308E−051.4961.50.7451240.74577580.74524820.7451535.367E−04Diffi5.275E−049.460E−052.947E−051.9211.80.9767480.9781670.9770320.97681791.157E−03Diffi1.13E−032.143E−046.891E−054000.84310.41077320.41085340.41078390.41077517.01E−05Diffi6.952E−058.795E−061.913E−071.2661.30.6255550.6259230.6256210.6255703.071E−04Diffi30251E−045/087E−051.511E−051.5491.50.7744680.7751810.7746030.7745005.886–04Diffi5.787E−041.028E−043.179E−051.9531.80.9950410.9965200.9953350.9951121.207E−03Diffi1.184E−032.229E−047.14E−05(66)Diffi=∫01Wnli∗−∫01Wnli+2∗,residual=∑i=39Diffi2.


## 4.1. Free Vibration

In order to solve equation (51), the assumptions stated by El Kadiri et al. were used to simplify the nonlinear algebraic system obtained. The contribution coefficients have been calculated as well as the associated nonlinear frequency. The normalized first nonlinear mode Wnl∗x,a1, for a given value of a1, is obtained as a series of basic functions involving the beam modal parameters depending on w1∗,…,w10∗ and given as follows:(59)Wnl∗x,a1=a1w1∗+⋯+ε10w10∗.The predominant term, relative to the first linear mode shape, isa1w1∗, while the other terms proportional to the symmetric and antisymmetric higher modes ε2w2∗,ε3w3∗,…,ε10w10∗ present the corrections due to the geometrical nonlinearity.Table1 shows the result of contribution coefficients obtained using equation (51); these contribution coefficients are included in equation (59) in order to calculate the normalized first nonlinear mode Wnl∗x,a1. Its maximal value is denoted by Wmax and is located at the middle of the beam. From Table 1, some observations can be clearly made:(i)
The contribution of the fundamental modea1 remains predominant for the full range of the vibration amplitudes considered. The other contributions denoted by ε2,…,εn remain small with respect to a1.(ii)
The contribution coefficients of the modes decrease asi indices increase.(iii)
The contribution coefficients increased according to the maximum nondimensional amplitudeWmax∗.Table 1
The linear and nonlinear frequency ratio and the contribution coefficients of stepped beam (b) subjected to a free vibration.
Wmaxa1ωnl∗/ωl∗ε2ε3ε4ε5ε6ε7ε8ε9ε100.9680.51.233.43E−03−4.21E031.42E−04−4.52E−04−2.912E−05−1.04E−04−2.77E−05−2.96E−051.39E−051.1660.61.325.64E−03−7.01E−032.69E−04−7.74E−04−4.449E−05−1.79E−04−4.61E−05−5.07E−052.36E−051.3650.71.428.48E−2−1.06E−024.65E−03−1.21E−03−6.058E−05−2.83E−04−7.03E−05−7.98E−053.66E−051.5660.81.531.19E−21.51E−027.48E−04−1.78E−03−7.46E−05−4.19E−03−1.00E−041.17E−045.32E−05From these observations, it is quite interesting to verify the convergence of the nonlinear beam deflection equation, which is obtained as a series of basic functions. The increase of these basic functions added in equation (59) enhanced the correction due to geometrical nonlinearity. The effect of varying the number of basic function series on the convergence of the numerical method has been investigated by comparing the results of a nonlinear frequency estimated ωeN∗. This latter can be expressed as in [26]:(60)ωeN∗=ATK∗A+2/3ATB∗AAATM∗A,forA=a1,ε2w2∗,…,εNwN∗.The nonlinear frequency estimatedωeN∗ can be extracted from equation (60), where A presents the contribution coefficient vector already obtained from the second formulation, using N basic functions. For a1=0.5, the nonlinear frequency ωeN∗ is calculated for different cases, in which the number of basic functions N increased from 7 to 10. The same test is repeated for different values of a1, varying from 0.5 to 0.8. The differences between ωeN∗ and ωeN-1∗ are calculated, as well as the residual of each test, and summarized hereafter in Table 2.Table 2
The frequency differences and residual for each value ofa1 and at different numbers of basic functions.
Basic function numbersTesta1ωnl∗/ωl∗78910Residual10.51.231.45E−051.80E−072.90E−071.00E−081.452E−0520.61.321.08E−044.50E−067.00E−064.00E−071.079E−0430.71.422.30E−049.60E−061.50E−059.00E−072.305E−0440.81.534.42E−041.90E−052.90E−051.00E−064.433−04(61)DiffN=ωN−1∗−ωN∗,residual=∑i=710DiffN2.The numerical results presented in Table2 show that the differences DiffN for the four tests decrease by increasing the number of basic functions, while the residual increases by increasing the a1 value. The difference Diff7 of the first test is equal to 1.45×10−05, and after increasing the number of basic functions to 10, this difference decreases to a low value of 1×10−08, meaning that(i)
The use of 6 basic functions induces an error of1.5×10−05 with respect to ωe10∗(ii)
The use of 7 basic functions induces an error of4.8×10−07 with respect to ωe10∗(iii)
The use of 8 basic functions induces an error of3×10−07 with respect to ωe10∗(iv)
The use of 9 basic functions induces an error of1×10−08 with respect to ωe10∗Also, ata1 = 0.8, the difference for 7 basic functions increases from 1.45×10−05 to 4.42×10−04, and by adding 3 other terms, the difference decreased to Diff10=1×10−06. In other terms,(i)
The use of 6 basic functions induces an error of4.91×10−04 with respect to ωe10∗(ii)
The use of 7 basic functions induces an error of4.90×10−05 with respect to ωe10∗(iii)
The use of 8 basic functions induces an error of3.00×10−05 with respect to ωe10∗(iv)
The use of 9 basic functions induces an error of1×10−06 with respect to ωe10∗When the basic functions’ number increases and exceeds ten, the difference between the final iteration ofωeN∗ and its precedent decreases by value less than 1×10−08 for the first test and less than 1×10−06 for the last one, whereas the residual will be increased by less than these values. These corrections are almost negligible compared to ωe∗. These results are considered to be tightly convergent in the nonlinear free vibration case using ten or more basic functions.


## 4.2. Forced Vibration

These analyses deal with a two-step fully clamped beam of type (b), whereR1∗=R3∗=0.5×R2∗, u1=u3, and u2=0.6×u1), subjected to a uniformly distributed force over an interval of [0.4, 0.6]. Firstly, the contribution coefficients are obtained using multimode approach (58) and then summarized in Table 3.Table 3
The linear and nonlinear frequency ratio and the contribution coefficients of stepped beam (b) subjected to a forced vibration.
Wmaxωnl∗/ωl∗a1ε2ε3ε4ε5ε6ε7ε8ε9ε100.0960.2990.05173.986E−185.258E−043.117E−194.381E−052.534E−191.175E−051.781E−184.455E−062.281E−180.7481.0990.39981.049E−153.6292E−3−1.027E−16−3.855E−044.186E−17−5.581E−05−7.892E−17−1.487E−054.443E−171.131.30.61313.271E−151.092E−2−3.642e−16−1.479E−31.369E−16−2.332E−04−2.810E−16−6.657E−051.459E−16Afterwards, these contribution coefficients are included in equation (59) in order to calculate the normalized first nonlinear mode Wnl∗x,a1. Its maximal value is denoted by Wmax and is located at the middle of the beam. The same observations can be concluded from Table 3. Furthermore, it can be clearly seen that, in the forced case, the contribution coefficients of the even modes have very low values. However, a test has been carried out to predict the effect of neglecting the antisymmetric terms in equation (59) on the normalized first nonlinear mode. Two first normalized nonlinear modes have been calculated: the first one is obtained as shown in equation (59), and the second is also obtained as a series of basic functions involving only the symmetrical terms ε3w3∗,…,ε9w9∗ and denoted by WnlSymmetric∗. The results are summarized in Table 4.Table 4
Comparison between normalized first nonlinear modesWnl∗ and WnlSymmetric∗ for ωnl∗/ωl∗=1.8.
W∗∫01W∗d2W∗/dx∗2|x∗=0.5∫01d2W∗/dx∗2dxWnl∗1.953070559600.99504143259212831.92929395580.005546040320137WnlSymmetric∗1.953070559600.995041432592128431.92929395580.005546040320142Difference—4.4408E−16—5.381E−15The numerical results presented in Table4 prove that, by neglecting the antisymmetric modes, the use of equation (59) leads to the same results. The normalized first nonlinear mode can also be written as(62)Wnl∗x,a1=a1w1∗+ε3w3∗+ε5w5∗+ε7w7∗+ε9w9∗.The effect of varying the basic functions’ number on the convergence of equation (62) has been investigated by comparing the normalized first nonlinear mode obtained for five basic functions Wnl∗ for three scenarios.Scenario 1.
The nonlinear beam deflection assumes the use of one basic function:(63)Wnl3∗x,a1=a1w1∗+ε3w3∗.Scenario 2.
The nonlinear beam deflection assumes the use of two basic functions:(64)Wnl5∗x,a1=a1w1∗+ε3w3∗+ε5w5∗.Scenario 3.
The nonlinear beam deflection assumes the use of three basic functions:(65)Wnl7∗x,a1=a1w1∗+ε3w3∗+ε5w5∗+ε7w7∗.
The convergence betweenWnl∗ and Wnl3∗,Wnl5∗,andWnl7∗ has been examined by calculating the integral of these nonlinear modes over an interval of [0, 1] using the Simpson numerical integration method. The results are obtained for different values of ωnl∗/ωl∗ in the case where Fd is equal to 100, 200, and 300. The comparison was performed and is summarized in Table 5.
In Table5, the final normalized first nonlinear mode Wnl∗ (62) was compared with the three other scenarios of nonlinear modes as shown in equations (63)–(65). The comparison was performed by varying the frequency ratio for each excitation. The results obtained for the scenario where Fd=100 and ωnl∗/ωl∗=1 show a small difference of Diff3=1.724×10−5 for the first scenario, and it decreases to Diff7=4.5×10−7 for the last one, in other terms:(i)
The use ofScenario 1 induces an error of 1.98×10−05 with respect to ∫01Wnl∗(ii)
The use ofScenario 2 induces an error of 2.59×10−06 with respect to ∫01Wnl∗(iii)
The use ofScenario 3 induces an error of 4.5023×10−07 with respect to ∫01Wnl∗
Also, for the same excitation, ifωnl∗/ωl∗=1.8, Diff3 increases until reaching 1.11×10−3 for the first scenario, and Diff7 increases to 6.7609×10−5 for the last one, meaning that(i)
The use ofScenario 1 induces an error of 1.39×10−03 with respect to ∫01Wnl∗(ii)
The use ofScenario 2 induces an error of 2.78×10−04 with respect to ∫01Wnl∗(iii)
The use ofScenario 3 induces an error of 6.76×10−05 with respect to ∫01Wnl∗
It can be noticed that the change in the excitation level from 100 to 400 withωnl∗/ωl∗=1.8 increases the error involved between Wnl∗ and Wnl7∗ by a small variation. For Fd=100, this error is equal to 6.7609×10−5, and for Fd=400, this error increases to 7.148×10−5. Indeed, the effect of the frequency ratio variation is more remarkable than the force variation effect on the convergence of the normalized first nonlinear mode. The differences between Wnl∗ and Wnl7∗ for the value of the frequency ratio and for the three excitations considered remain small. However, the increase in the number of basic functions exceeding five symmetric terms in equation (62) provides a correction, due to the geometric nonlinearity, less than 7.148×10−5. This correction is almost negligible compared to Wnlmax∗; hence, the numerical method proposed is well converged.Table 5
The convergence test of the normalized first nonlinear modes with different numbers of basic functions in the forced vibration case.
FdWnlmax∗ωnl∗/ωl∗∫01Wnl∗∫01Wnl3∗∫01Wnl5∗∫01Wnl7∗Residual1000.53610.25807410.2580940.2580760.25807461.738E−05Diffi1.724E−052.135E−064.5023E−071.1351.30.5566270.556900.5566790.5566392.302E−04Diffi2.265E−043.934E−051.201E−061.461.50.7291680.729780.72928690.72919655.09E−04Diffi5.010E−049.033E−052.826E−051.91.80.967220.96861680.96750640.9672961.132E−03Diffi1.11E−032.10E−046.7609E−052000.67310.3255040.325544260.32550960.3255053.489E−05Diffi3.46E−054.321E−069.229E−071.1841.30.58206890.582378150.58212530.582082032.568E−04Diffi2527E−044.333E−051.308E−051.4961.50.7451240.74577580.74524820.7451535.367E−04Diffi5.275E−049.460E−052.947E−051.9211.80.9767480.9781670.9770320.97681791.157E−03Diffi1.13E−032.143E−046.891E−054000.84310.41077320.41085340.41078390.41077517.01E−05Diffi6.952E−058.795E−061.913E−071.2661.30.6255550.6259230.6256210.6255703.071E−04Diffi30251E−045/087E−051.511E−051.5491.50.7744680.7751810.7746030.7745005.886–04Diffi5.787E−041.028E−043.179E−051.9531.80.9950410.9965200.9953350.9951121.207E−03Diffi1.184E−032.229E−047.14E−05(66)Diffi=∫01Wnli∗−∫01Wnli+2∗,residual=∑i=39Diffi2.


## 5. Numerical Results and Discussion

### 5.1. Linear Free Vibration

A computer program has been performed using MATLAB software, which analyses the effects of the number and position of each step of the multiple-stepped beam, as well as the change in the ratios of the step cross sections and moments of inertia on the vibration frequencies and mode shapes of the beam. To verify the validity and precision of this analysis, three cases are examined first, named beam (a), beam (b), and beam (c), in which the beam cross sections vary in width, in depth, and in both. Secondly, the three beams are assumed to be fixed on both sides; for this purpose, an infinite value is adopted for the value of the stiffness of the linear and rotary springs at both ends (kt1=∞,kt2=∞,kθ1=∞, and kθ2=∞). Thirdly, the values of the vibration frequencies for the first four vibration modes are calculated and compared with those available in [6, 13, 14] as shown in Tables 6 and 7.Table 6
The first four eigenvalues for the three clamped-clamped beams (a, b, and c) where (λ1;λ2;λ3;λ4) = (1, 0; 0, 80; 0, 65; 0, 25) and (R1; R2; R3; R4) = (0, 25; 0, 30; 0, 25; 0, 20).
Mode number1234Beam (a)Present work4.54057.660310.808814.0643[6]4.54057.660310.8088—Beam (b)Present work3.54195.81168.481811.2834[14]3.54195.81168.481811.2834[6]3.54195.81168.4818—Beam (c)Present work3.31285.73168.584111.4491[6]3.31285.73168.5841—Table 7
The first four eigenvalues for the three clamped-clamped beams (a, b, and c) where (λ1;λ2) = (1, 0; 0, 50) and (R1; R2) = (0, 25; 0, 75).
Mode number1234Beam (a)Present work4.95797.904210.970814.1707[6]4.95797.904310.9708—[13]4.95797.904310.970814.1717Beam (b)Present work3.99066.16178.305310.7570[6]3.99066.16188.3054—[13]3.99066.16188.305410.7577Beam (c)Present work4.17666.38238.215410.6443[6]4.17666.38248.2155—[13]4.17666.38248.215510.6448Finally, the program developed for the case of tapered beams is validated by modelling the tapered beam with a 12-step beam as shown in Figure3, maintaining the same taper ratio and the same angle of inclination.Figure 3
Clamped-clamped tapered beams with a taper ratio equal to 0.2.The results obtained are compared in Table8 with those of [34] and show an excellent agreement, with a relative difference not exceeding 0.0094% for the whole range of frequencies considered. The increase in the number of steps is expected to further reduce this difference.Table 8
The first four eigenvalues for clamped-clamped tapered beams varying in depth with a taper ratio denoted byα: α=0.2 and Sx/S1=1−αx/L.
Mode index1234Present work4.48127.439710.416413.3922[34]4.48087.439010.415513.3911Rel diff (%)0.00940.0090.0080.007


### 5.2. Nonlinear Free Vibration

The nonlinear vibration of tapered beams similar to that of [33–35] is now investigated in the same way, and a comparison is made between the results obtained here, based on the second formulation and summarized in Table 9, and those of [33–35].Table 9
ωNL/ω for clamped-tapered beams in which beam (a) is varying in breadth and beam (b) is varying in depth, with a taper ratio α=0.4 and r=I/S.
Wmax/rPresent work[33]Rel diff [33] (%)[35]Rel diff [35] (%)Beam (a)0.21.00091.000901.000900.41.00361.00370.009——0.61.0081.00820.011.00790.0090.81.01431.01460.02——11.02231.02270.031.02180.0421.08591.08790.181.08350.22Beam (b)0.21.00141.00170.021.00150.0080.41.00561.00660.09——0.61.01261.01490.221.01340.070.81.02231.02630.38——11.03471.04080.581.03660.1721.13361.15451.81.13660.25An excellent agreement can be noticed since the relative difference does not exceed 0.2% for beam (a) and 1.8% for beam (b). Increasing the number of steps is expected to further reduce this error.


### 5.3. The Effect of the Section Ratio

The effect of the cross-section ratio and the vibration amplitude on the beam nonlinear behaviour is illustrated in Figures4–6 by means of the backbone curves associated to the first nonlinear mode shape of a fully fixed beam with two steps (R1∗=R3∗=0.6×R2∗, u1=u3, and u2=m×u1), where the value of m is increasing for the three beams (a, b, and c).Figure 4
Backbone curves of beam (a) in the vicinity of the first mode for various values ofm.Figure 5
Backbone curves of beam (b) in the vicinity of the first mode for various values ofm.Figure 6
Backbone curves of beam (c) in the vicinity of the first mode for various values ofm.It is clear that the nonlinearity remains of the hardening type and that the increase in the section ratio (u2) leads to an increase in the frequency ratio.


### 5.4. The Effect of the Moment of Inertia Ratio

A comparison between the three beam configurations with the same section ratio (u1=u3 and u2=0.5×u1) and length ratio (R1∗=R3∗=0.5×R2∗) but with a different inertia ratio is performed and illustrated in Figures 7–9.Figure 7
The curvature distribution of beam (a). Maximum nondimensional amplitudes (Wmax∗=1.5).Figure 8
The curvature distribution of beam (b). Maximum nondimensional amplitudes (Wmax∗=1.5).Figure 9
The curvature distribution of beam (c). Maximum nondimensional amplitudes (Wmax∗=1.5).Figures7–9 show the stress distributions along the beam for each configuration, derived from the linear theory and the nonlinear model. The effect of geometric nonlinearity is clearly shown in the figures. The calculation of the percentage correction introduced by the nonlinear theory compared to the linear theory shows that it can reach 36% for beam (a), 23% for beam (b), and 59% for beam (c) as summarized in Table 10.Table 10
The percentage correction between the curvatures estimated via the linear and nonlinear theories. Maximum nondimensional amplitudes (Wmax∗=1.5).
BeamCurvaturex∗NonlinearLinearPercentage correction (%)(a)045.55538.6678.880.24.4997.07436.390.5−26.518−29.47810.04(b)030.59427.53611.1070.2527.95636.77823.9880.5−42.4−44.3144.3196(c)037.34533.30312.1370.255.93114.60659.3910.5−32.104−35.2889.022In addition, it can be noticed that the maximum value of the curvature for the linear and nonlinear cases for beam (a) is located near the clamped ends, while it is located in the middle for beam (b). For beam (c), the location of the maximum curvature is not the same in the linear and nonlinear cases: in the linear case, it is located in the middle of the beam, while in the nonlinear case, it is located at its end.


### 5.5. The Effect of the Step Length Ratio

The effect of the step length ratio on the nonlinear behaviour of beam (b) is illustrated in Figure10 which gives the backbone curves associated to the first nonlinear mode shape for four cases: R1∗=R3∗≠R2∗, u1=u3, and u2=0.5×u1.Figure 10
Backbone curves of beam (b) in the vicinity of the first mode for various values ofR2.Figure10 shows the effect of the length ratios on the geometrically nonlinear beam behaviour. The modification of these ratios is obtained by increasing R2∗, resulting in a decrease in R1∗ and R3∗. It should be noted that, by increasing the length ratio R2∗, the frequency ratio increases accordingly.


### 5.6. Nonlinear Forced Vibration

The forced nonlinear dynamic behaviour is investigated for the three types of beams (a, b, and c) withR1∗=R3∗=0.5×R2∗, u1=u3, and u2=0.6×u1 and excited by a uniformly distributed force over an interval of [0.4, 0.6]. In order to examine the forced case, and as mentioned in the theoretical section, it is required to determine the predominant mode corresponding to each level of excitation in order to perform the study in its vicinity. The generalized forces have thereby been calculated and listed in Table 11.Table 11
Percentage of generalized forces exciting the first five symmetric modes of three types of beams excited by a uniformly distributed force over an interval of [0.4, 0.6].
BeamModes13579∫0.40.6wi∗/∑110∫0.40.6wi∗(a)53.216.68.4211.410.1(b)52.518.18.929.8210.5(c)52.717.28.7510.310.9The calculation of the generalized forces led to the conclusion that the first mode is predominant. Figures11–14 show the nonlinear frequency response functions in the vicinity of the first vibration mode.Figure 11
Linear and nonlinear frequency response functions’ comparison for beam (b).Figure 12
Nonlinear frequency response functions, based on a multimode approach, for the three different beams, which are excited by a uniformly distributed harmonic force.Figure 13
Nonlinear frequency response functions, based on a multimode approach, for beam (b) with differentR2∗.Figure 14
Nonlinear frequency response functions, based on a multimode approach, for beam (b) excited by different levels of distributed forces.Figure11 shows a comparison between the linear and nonlinear frequency response functions of configuration 2, where the beam is excited by a distributed force Fd=100 uniformly distributed over an interval of [0.4, 0.6].Figure12 shows a comparison between the frequency response curves obtained for the three beams considered, having different inertias. The comparison is conducted using a multimodal approach for three beams subjected to a uniformly distributed force over an interval of [0.4, 0.6], with Fd=100.Figure13 shows a comparison between the response curves obtained for the three configurations: R1∗=R3∗≠R2∗, u1=u3, and u2=0.6×u1. The beam is subjected to a uniformly distributed force Fd=100 in the interval of [0.4, 0.6].Figure14 presents a comparison between the frequency response curves of beam (b), excited by different levels of uniformly distributed forces in the range [0.4, 0.6]. The effect of geometrical nonlinearity can be seen with the multivalued region in which the jump phenomenon commonly found in nonlinear vibration systems can occur. The hardening behaviour is also clearly visible, as shown in Figure 11. As the vibration amplitude increases, the frequency ratio increases accordingly. In addition, Figure 12 shows how the nonlinear effect changes with the type of cross section when one moves from one beam to another. It appears from this figure that the hardening effect in beam (b) is stronger than in the other two. Figure 13 shows the effect of the change in length ratio on the geometrically nonlinear behaviour. The increase in the length ratio R2∗ makes the hardening behaviour more accentuated. Figure 14 illustrates the disproportionate increase in the frequency response amplitude relative to the increase in the intensity of the distributed excitation force.


## 5.1. Linear Free Vibration

A computer program has been performed using MATLAB software, which analyses the effects of the number and position of each step of the multiple-stepped beam, as well as the change in the ratios of the step cross sections and moments of inertia on the vibration frequencies and mode shapes of the beam. To verify the validity and precision of this analysis, three cases are examined first, named beam (a), beam (b), and beam (c), in which the beam cross sections vary in width, in depth, and in both. Secondly, the three beams are assumed to be fixed on both sides; for this purpose, an infinite value is adopted for the value of the stiffness of the linear and rotary springs at both ends (kt1=∞,kt2=∞,kθ1=∞, and kθ2=∞). Thirdly, the values of the vibration frequencies for the first four vibration modes are calculated and compared with those available in [6, 13, 14] as shown in Tables 6 and 7.Table 6
The first four eigenvalues for the three clamped-clamped beams (a, b, and c) where (λ1;λ2;λ3;λ4) = (1, 0; 0, 80; 0, 65; 0, 25) and (R1; R2; R3; R4) = (0, 25; 0, 30; 0, 25; 0, 20).
Mode number1234Beam (a)Present work4.54057.660310.808814.0643[6]4.54057.660310.8088—Beam (b)Present work3.54195.81168.481811.2834[14]3.54195.81168.481811.2834[6]3.54195.81168.4818—Beam (c)Present work3.31285.73168.584111.4491[6]3.31285.73168.5841—Table 7
The first four eigenvalues for the three clamped-clamped beams (a, b, and c) where (λ1;λ2) = (1, 0; 0, 50) and (R1; R2) = (0, 25; 0, 75).
Mode number1234Beam (a)Present work4.95797.904210.970814.1707[6]4.95797.904310.9708—[13]4.95797.904310.970814.1717Beam (b)Present work3.99066.16178.305310.7570[6]3.99066.16188.3054—[13]3.99066.16188.305410.7577Beam (c)Present work4.17666.38238.215410.6443[6]4.17666.38248.2155—[13]4.17666.38248.215510.6448Finally, the program developed for the case of tapered beams is validated by modelling the tapered beam with a 12-step beam as shown in Figure3, maintaining the same taper ratio and the same angle of inclination.Figure 3
Clamped-clamped tapered beams with a taper ratio equal to 0.2.The results obtained are compared in Table8 with those of [34] and show an excellent agreement, with a relative difference not exceeding 0.0094% for the whole range of frequencies considered. The increase in the number of steps is expected to further reduce this difference.Table 8
The first four eigenvalues for clamped-clamped tapered beams varying in depth with a taper ratio denoted byα: α=0.2 and Sx/S1=1−αx/L.
Mode index1234Present work4.48127.439710.416413.3922[34]4.48087.439010.415513.3911Rel diff (%)0.00940.0090.0080.007


## 5.2. Nonlinear Free Vibration

The nonlinear vibration of tapered beams similar to that of [33–35] is now investigated in the same way, and a comparison is made between the results obtained here, based on the second formulation and summarized in Table 9, and those of [33–35].Table 9
ωNL/ω for clamped-tapered beams in which beam (a) is varying in breadth and beam (b) is varying in depth, with a taper ratio α=0.4 and r=I/S.
Wmax/rPresent work[33]Rel diff [33] (%)[35]Rel diff [35] (%)Beam (a)0.21.00091.000901.000900.41.00361.00370.009——0.61.0081.00820.011.00790.0090.81.01431.01460.02——11.02231.02270.031.02180.0421.08591.08790.181.08350.22Beam (b)0.21.00141.00170.021.00150.0080.41.00561.00660.09——0.61.01261.01490.221.01340.070.81.02231.02630.38——11.03471.04080.581.03660.1721.13361.15451.81.13660.25An excellent agreement can be noticed since the relative difference does not exceed 0.2% for beam (a) and 1.8% for beam (b). Increasing the number of steps is expected to further reduce this error.


## 5.3. The Effect of the Section Ratio

The effect of the cross-section ratio and the vibration amplitude on the beam nonlinear behaviour is illustrated in Figures4–6 by means of the backbone curves associated to the first nonlinear mode shape of a fully fixed beam with two steps (R1∗=R3∗=0.6×R2∗, u1=u3, and u2=m×u1), where the value of m is increasing for the three beams (a, b, and c).Figure 4
Backbone curves of beam (a) in the vicinity of the first mode for various values ofm.Figure 5
Backbone curves of beam (b) in the vicinity of the first mode for various values ofm.Figure 6
Backbone curves of beam (c) in the vicinity of the first mode for various values ofm.It is clear that the nonlinearity remains of the hardening type and that the increase in the section ratio (u2) leads to an increase in the frequency ratio.


## 5.4. The Effect of the Moment of Inertia Ratio

A comparison between the three beam configurations with the same section ratio (u1=u3 and u2=0.5×u1) and length ratio (R1∗=R3∗=0.5×R2∗) but with a different inertia ratio is performed and illustrated in Figures 7–9.Figure 7
The curvature distribution of beam (a). Maximum nondimensional amplitudes (Wmax∗=1.5).Figure 8
The curvature distribution of beam (b). Maximum nondimensional amplitudes (Wmax∗=1.5).Figure 9
The curvature distribution of beam (c). Maximum nondimensional amplitudes (Wmax∗=1.5).Figures7–9 show the stress distributions along the beam for each configuration, derived from the linear theory and the nonlinear model. The effect of geometric nonlinearity is clearly shown in the figures. The calculation of the percentage correction introduced by the nonlinear theory compared to the linear theory shows that it can reach 36% for beam (a), 23% for beam (b), and 59% for beam (c) as summarized in Table 10.Table 10
The percentage correction between the curvatures estimated via the linear and nonlinear theories. Maximum nondimensional amplitudes (Wmax∗=1.5).
BeamCurvaturex∗NonlinearLinearPercentage correction (%)(a)045.55538.6678.880.24.4997.07436.390.5−26.518−29.47810.04(b)030.59427.53611.1070.2527.95636.77823.9880.5−42.4−44.3144.3196(c)037.34533.30312.1370.255.93114.60659.3910.5−32.104−35.2889.022In addition, it can be noticed that the maximum value of the curvature for the linear and nonlinear cases for beam (a) is located near the clamped ends, while it is located in the middle for beam (b). For beam (c), the location of the maximum curvature is not the same in the linear and nonlinear cases: in the linear case, it is located in the middle of the beam, while in the nonlinear case, it is located at its end.


## 5.5. The Effect of the Step Length Ratio

The effect of the step length ratio on the nonlinear behaviour of beam (b) is illustrated in Figure10 which gives the backbone curves associated to the first nonlinear mode shape for four cases: R1∗=R3∗≠R2∗, u1=u3, and u2=0.5×u1.Figure 10
Backbone curves of beam (b) in the vicinity of the first mode for various values ofR2.Figure10 shows the effect of the length ratios on the geometrically nonlinear beam behaviour. The modification of these ratios is obtained by increasing R2∗, resulting in a decrease in R1∗ and R3∗. It should be noted that, by increasing the length ratio R2∗, the frequency ratio increases accordingly.


## 5.6. Nonlinear Forced Vibration

The forced nonlinear dynamic behaviour is investigated for the three types of beams (a, b, and c) withR1∗=R3∗=0.5×R2∗, u1=u3, and u2=0.6×u1 and excited by a uniformly distributed force over an interval of [0.4, 0.6]. In order to examine the forced case, and as mentioned in the theoretical section, it is required to determine the predominant mode corresponding to each level of excitation in order to perform the study in its vicinity. The generalized forces have thereby been calculated and listed in Table 11.Table 11
Percentage of generalized forces exciting the first five symmetric modes of three types of beams excited by a uniformly distributed force over an interval of [0.4, 0.6].
BeamModes13579∫0.40.6wi∗/∑110∫0.40.6wi∗(a)53.216.68.4211.410.1(b)52.518.18.929.8210.5(c)52.717.28.7510.310.9The calculation of the generalized forces led to the conclusion that the first mode is predominant. Figures11–14 show the nonlinear frequency response functions in the vicinity of the first vibration mode.Figure 11
Linear and nonlinear frequency response functions’ comparison for beam (b).Figure 12
Nonlinear frequency response functions, based on a multimode approach, for the three different beams, which are excited by a uniformly distributed harmonic force.Figure 13
Nonlinear frequency response functions, based on a multimode approach, for beam (b) with differentR2∗.Figure 14
Nonlinear frequency response functions, based on a multimode approach, for beam (b) excited by different levels of distributed forces.Figure11 shows a comparison between the linear and nonlinear frequency response functions of configuration 2, where the beam is excited by a distributed force Fd=100 uniformly distributed over an interval of [0.4, 0.6].Figure12 shows a comparison between the frequency response curves obtained for the three beams considered, having different inertias. The comparison is conducted using a multimodal approach for three beams subjected to a uniformly distributed force over an interval of [0.4, 0.6], with Fd=100.Figure13 shows a comparison between the response curves obtained for the three configurations: R1∗=R3∗≠R2∗, u1=u3, and u2=0.6×u1. The beam is subjected to a uniformly distributed force Fd=100 in the interval of [0.4, 0.6].Figure14 presents a comparison between the frequency response curves of beam (b), excited by different levels of uniformly distributed forces in the range [0.4, 0.6]. The effect of geometrical nonlinearity can be seen with the multivalued region in which the jump phenomenon commonly found in nonlinear vibration systems can occur. The hardening behaviour is also clearly visible, as shown in Figure 11. As the vibration amplitude increases, the frequency ratio increases accordingly. In addition, Figure 12 shows how the nonlinear effect changes with the type of cross section when one moves from one beam to another. It appears from this figure that the hardening effect in beam (b) is stronger than in the other two. Figure 13 shows the effect of the change in length ratio on the geometrically nonlinear behaviour. The increase in the length ratio R2∗ makes the hardening behaviour more accentuated. Figure 14 illustrates the disproportionate increase in the frequency response amplitude relative to the increase in the intensity of the distributed excitation force.


## 6. Conclusion

The geometrically nonlinear free and forced vibrations of multistepped beams were examined, based on Euler–Bernoulli’s beam theory and von Karman’s nonlinearity assumptions. The solution of the linear problem was obtained for different configurations of multistepped beams. The calculation of the frequencies was performed using the Newton–Raphson algorithm, and the mode shapes were determined afterwards. Discrete expressions of the strain and kinetic energies were derived, and the problem was reduced using Hamilton’s principle to a nonlinear algebraic system solved by the second formulation. The numerical method convergence study was carried out in the free and forced vibrations, in which the results showed a good convergence. The results obtained were compared with those obtained using the Galerkin method in order to verify the validity and precision of the proposed approach and showed a good agreement since the relative difference does not exceed 1.8% forWmax∗=2. The effect of the geometrical nonlinearity on the dynamic behaviour of multistepped beams was illustrated, as well as the effects of the change in the beam cross section, inertia, and length ratios. It was clearly noticed that the nonlinearity remains of the hardening type and that the increase in the section or the length ratios leads to the increase of the frequency ratio. Furthermore, it appears that the location of the maximum curvature can differ from the linear to the nonlinear cases. The forced vibration was analyzed using the multimode approach. The effect of the excitation level, the type of section, and the length ratio on the nonlinear forced behaviour of the beam was illustrated and discussed. Besides, it can be clearly seen that the hardening effect is greatest when the cross section of the beams increases in height or when the second length ratio increases. Furthermore, a convergence study was also presented to show how the results of the numerical method are well convergent.



---

*Source: 6697344-2021-04-16.xml*