A Semianalytical Method for Nonlinear Vibration of Euler-Bernoulli Beams with General Boundary Conditions

Abstract

This paper presents a new semianalytical approach for geometrically nonlinear vibration analysis of Euler-Bernoulli beams with different boundary conditions. The method makes use of Linstedt-Poincaré perturbation technique to transform the nonlinear governing equations into a linear differential equation system, whose solutions are then sought through the use of differential quadrature approximation in space domain and an analytical series expansion in time domain. Validation of the present method is conducted in numerical examples through direct comparisons with existing solutions, showing that the proposed semianalytical method has excellent convergence and can give very accurate results at a long time interval.

1. Introduction

Geometrically nonlinear vibration of beams with different boundary conditions has long been a subject receiving numerous research efforts, as evidenced by many analytical and numerical studies reported in the open literature. Woinowsky-Krieger [1] used the elliptic integral function to solve the nonlinear equation of motion for simply supported beams with immovable ends. Lewandowski [2] applied the Rayleigh-Ritz technique to the nonlinear vibration of beams. A comprehensive review in this field was given by Sathyamoorthy [3]. Sze et al. [4] applied incremental harmonic balance method for nonlinear vibration of axially moving beams. Gadagi and Benaroya [5] studied the dynamic response of an axially loaded tendon of a tension leg platform. Ibrahim and Somnay [6] solved the nonlinear dynamic analysis of an elastic beam isolator sliding on frictional supports. Ozkaya and Tekin [7] analyzed the nonlinear vibrations of stepped beam system under different boundary conditions. Chen et al. [8] put forward the multidimensional Lindstedt-Poincare method for nonlinear vibration of axially moving beams. It is noted that many numerical investigations published so far employed step-by-step iterative time integration schemes to calculate the dynamic deflection response which requires a very small time step size in order to achieve the response with sufficient accuracy. This is inevitably computationally expensive. Another big concern associated with these iterative schemes is the error accumulation in each time step which may cause a huge loss of accuracy and sometimes numerical stability problem of the results in later time steps. To overcome this problem, several analytical and semianalytical approaches to the nonlinear vibration of beams which do not involve step-by-step time integration have been proposed; see, for example, those by Azrar et al. [9, 10], Leung and his coworkers [11–13], Lewandowski [14, 15], Nayfeh and his associates [16–18], and Ribeiro [19], to name just a few.

The differential quadrature method (DQM) is an efficient numerical approach developed by Bellman et al. [20] to solve linear and nonlinear differential equations. This method was later introduced to the structural analysis by Bert and his coworkers [21] and has been widely used in many structural engineering problems, including those in the nonlinear vibration of beams, plates, and shells; see, for example, the work by Zhong and Guo [22], Yang et al. [23, 24], Manaoach and Ribeiro [25], Hsu [26], and Tomasiello [27], among many others. Based on DQM approximation, this paper proposes a new semianalytic method for the geometrically nonlinear vibration of Euler-Bernoulli beams with different boundary conditions. The proposed method overcomes the disadvantage of error accumulation that occurs in conventional numerical methods involving iterative time integration and is therefore accurate and numerically stable even for a long time interval. The present paper is structured as follows: Section 2 briefly outlines the DQM rules. Section 3 gives a detailed description of mathematic formulations of the proposed semianalytical approach, starting with the transformation of the nonlinear governing equations into a group of linear differential equations by Linstedt-Poincare perturbation procedure which is followed by the semianalytical solution process for each perturbation equation by using DQM in space domain and an analytical series in the time domain. Section 4 presents some numerical results to validate the proposed method in both accuracy and convergence aspects through direct comparisons between the present results and existing solutions. Some concluding remarks are summarized in Section 4.

2. Mathematical Formulations

2.1. Perturbation Equations

Consider an isotropic, homogeneous slender beam of length L with uniform cross-section area subjected to a dynamic transverse load 𝑞(𝑥,𝑡). Let 𝑢 and 𝑤 be the displacements parallel to the 𝑥- and 𝑦-axes, 𝑡 the time, and 𝑚 denote the mass density. Based on classical Euler-Bernoulli beam theory and von-Karman nonlinear displacement-strain relationship, the geometrically nonlinear governing equations for flexural vibration of the beam can be derived as 𝜕𝜕𝑥[𝐸𝜕𝐴(𝑢+1𝜕𝑥2(𝜕𝑤)𝜕𝑥²𝜕)]=0,𝐸𝐼⁴𝑤𝜕𝑥⁴−𝐸𝜕𝐴[𝑢+1𝜕𝑥2(𝜕𝑤)𝜕𝑥²]𝜕²𝑤𝜕𝑥²𝜕+𝑚²𝑤𝜕𝑡²=𝑞(𝑥,𝑡), (2.1) where 𝐸 is the elastic modulus of the beam, and 𝐴 and 𝐼 are the area and the second moment of the cross section. The present analysis considers general boundary conditions, that is, the beam end may be either simply supported or clamped, with the following boundary conditions: 𝑢=0,𝑤=0,𝑀_𝑥=0(simplysupported),𝑢=0,𝑑𝑤=0,𝑤𝑑𝑥=0(clamped), (2.2) where 𝑀_𝑥 is the bending moment of the beam. To facilitate the solution process of (2.1) by using the differential quadrature method in space domain, the following quantities are introduced: 𝑢=𝑢𝜂,𝑤=𝑤𝜌𝑥,𝜉=𝐿,𝑡=𝜔𝐿𝑡,𝑞=⁴𝑞𝜔𝜌𝐸𝐼,𝑞=𝜀𝐹,𝑝=𝜔_𝐿,𝜔_𝐿1=(𝐿)²𝐸𝐼𝑚,𝜌²=𝐼𝐴, (2.3) where 𝜔_𝐿 and 𝜔 refer to the linear and nonlinear frequencies of the beam, respectively, and 𝜀 is the small perturbation parameter. Substituting these quantities into (2.1) leads to dimensionless nonlinear governing equations as below 𝜕²𝑢𝜕𝜉²+𝜕𝑤𝜕𝜕𝜉²𝑤𝜕𝜉²𝜕=0,₄𝑤𝜕𝜉⁴−[𝜕𝑢+1𝜕𝜉2(𝜕𝑤)𝜕𝜉²]𝜕²𝑤𝜕𝜉²+𝑝²𝜕²𝑤𝜕𝑡²=𝜀𝐹. (2.4)

Deflection 𝑤, axial displacement 𝑢, and frequency parameter 𝑝 can be expanded into series forms in terms of a small perturbation parameter 𝜀𝑤=𝑤₁𝜀+𝑤₃𝜀³+𝑤₅𝜀⁵+⋯,𝑢=𝑢₂𝜀²+𝑢₄𝜀⁴𝑝+⋯,₂=1+𝜀²𝜎₂+𝜀⁴𝜎₄+⋯. (2.5)(2.6)(2.7) Inserting (2.5)–(2.7) into (2.4), and equating the terms of the same power of 𝜀 yield a series of perturbation equations 𝜀¹:𝜕⁴𝑤₁𝜕𝜉⁴+𝜕²𝑤₁𝜕𝑡²𝜀=𝐹,₂:𝜕²𝑢₂𝜕𝜉²+𝜕𝑤₁𝜕𝜕𝜉²𝑤₁𝜕𝜉²𝜀=0,₃:𝜕⁴𝑤₃𝜕𝜉⁴−𝜕𝑢₂𝜕𝜕𝜉²𝑤₁𝜕𝜉²−12(𝜕𝑤₁)𝜕𝜉²𝜕²𝑤₁𝜕𝜉²+𝜕²𝑤₃𝜕𝑡²+𝜎₂𝜕²𝑤₁𝜕𝑡²𝜀=0,₄:𝜕²𝑢₄𝜕𝜉²+𝜕𝑤₁𝜕𝜕𝜉²𝑤₃𝜕𝜉²+𝜕𝑤₃𝜕𝜕𝜉²𝑤₁𝜕𝜉²=0. (2.8)(2.9)(2.10)(2.11)

The initial conditions considered in the present study are𝑤₁(𝜉,0)=𝑤₁₀,𝜕𝑤₁(𝜉,0)𝑤𝜕𝑡=0,₃(𝜉,0)=0,𝜕𝑤₃(𝜉,0)𝜕𝑡=0, (2.12)(2.13) where 𝑤₁₀ refers to initial vibration amplitude. By making use of the dimensionless quantities and perturbation series defined above, the dimensionless form of the associated boundary conditions in each perturbation can also be obtained.

2.2. Differential Quadrature Method

According to the DQM rule, an unknown function 𝑔(𝑥) and its 𝑛th order derivative with respect to x can be approximated as the weighted linear sums of its function values at a number of sampling points in the 𝑥-axis as 𝑔(𝑥)=^𝑁_𝑖=1𝑔(𝑥_𝑖)𝐿_𝑖𝑑(𝑥),_𝑛𝑔(𝑥)𝑑𝑥^𝑛|_{𝑥=𝑥𝑖}=^𝑁_𝑗=1𝐶^(𝑛)_𝑖𝑗𝑔(𝑥_𝑗)(𝑛=1,…,𝑁−1), (2.14)(2.15) where 𝑁 is the total number of sampling points, 𝐿_𝑖(𝑥) are Lagrange interpolation polynomials, and the weighting coefficients 𝐶^(𝑛)_𝑖𝑗 can be calculated from recursive formulae [20, 21].

2.3. Semianalytical Solution of the 1st-Order Perturbation Equation

This study is focused on the nonlinear vibration due to an initial deflection. In such a case, the right-hand-side term in (2.8) 𝐹=0. Hence, the first-order perturbation becomes𝜕⁴𝑤₁𝜕𝜉⁴+𝜕²𝑤₁𝜕𝑡²=0. (2.16) The solution of (2.16) can be readily obtained by separation of variables𝑤₁(𝜉,𝑡)=^𝑁𝑠_𝑚=0𝑇_1𝑚(𝑡)𝜙_1𝑚(𝜉), (2.17) where 𝑁_𝑠 is the number of terms in the truncated time series, 𝑇_1𝑚(𝑡)=cos(𝑝_𝑚𝑡) (2.18) is the solution in the time domain for the first-order perturbation equation, and 𝑝_𝑚 and 𝛿_1𝑚(𝜉) are the frequency parameter and the associated mode shape function.

The transient deflection at an arbitrary sampling point “𝑖’’ can be expressed as a cosine series with a supplementary term as 𝑤₁(𝜉_𝑖,𝑡)=𝛿^𝑖₀+^𝑁𝑠_𝑚=1𝑇_1𝑚(𝑡)𝛿^𝑖_1𝑚(𝑖=1,…,𝑁), (2.19) where 𝛿^𝑖_1𝑚=𝜙_1𝑚(𝜉_𝑖) needs to be determined. The frequency parameter 𝑝_𝑚 of beams with different boundary conditions can be written in a general form as 𝑝_𝑚=[(𝑚+𝑗)𝜋]² (2.20) in which 𝑗=0 for a simply supported beam, 𝑗=1/2 for a clamped beam, 𝑗=−1/2 for a beam clamped at one end and free at the other end, and 𝑗=1/4 for a beam clamped at one end but simply supported at the other end.

Since 𝑤₁(𝜉_𝑖,𝑡) needs to satisfy the initial conditions in (2.12), one has 𝛿^𝑖₀=𝑤^𝑖₁₀−^𝑁_𝑚=1𝛿^𝑖_1𝑚. (2.21) The time domain function in (2.19) can be readily obtained as𝑤₁(𝜉_𝑖,𝑡)=𝑤^𝑖₁₀+^𝑁𝑠_𝑚=1(cos𝑝_𝑚1𝑡−1)𝛿^𝑖_1𝑚(𝑖=1,…,𝑁). (2.22) Differentiation of (2.22) with respect to time twice gives𝑤̈₁(𝜉_𝑖,𝑡)=−^𝑁𝑠_𝑚=1𝑝²_𝑚(cos𝑝_𝑚1𝑡)𝛿^𝑖_1𝑚(𝑖=1,…,𝑁). (2.23)

The deflection response of the 1st-order perturbation equation can be determined by 𝑤₁(𝜉,𝑡)=^𝑁_𝑖=1𝑤₁(𝜉_𝑖,𝑡)𝐿_𝑖(𝜉). (2.24)

Substitution of (2.23) and application of DQM approximation (2.15) to (2.16) result in a semianalytical algebraic equation at an arbitrary sampling point^𝑁_𝑗=1𝐶^𝑁(4)_𝑖𝑗^𝑠_𝑚=1(cos𝑝_𝑚𝑡_𝑘−1)𝛿^𝑗_1𝑚−𝑝²₁^𝑁𝑠_𝑚=1𝑝²_𝑚cos(𝑝_𝑚𝑡_𝑘)𝛿^𝑖_1𝑚=−^𝑁_𝑗=1𝐶⁽⁴⁾_𝑖𝑗𝑤^𝑗₁₀(𝑖=1,…,𝑁;𝑘=1,…,𝑁_𝑠). (2.25) This linear equation system can be written in a matrix form as[𝐷₁]{𝛿₁}={𝑓₁}, (2.26) where [𝐷₁] and {𝛿₁} are the coefficient matrix and unknown vector to be solved, and {𝑓₁} is the generalized load vector due to initial deflection. Once {𝛿₁} is determined from (2.26), the deflections at all of the sampling points can be calculated from (2.22) and the 1st-order deflection solution of the beam is then obtained from (2.24).

2.4. Semianalytical Solution of the 2nd-Order Perturbation Equation

To solve the 2nd-order perturbation equation (2.9), the unknown axial displacement 𝑢₂ is expressed as𝑢₂(𝜉,𝑡)=^𝑁𝑠_𝑚=1𝑇_2𝑚(𝑡)𝜙_2𝑚(𝜉). (2.27) Substituting the 1st-order deflection solution and (2.27) into (2.9) yields ^𝑁𝑠_𝑚=1𝑇_2𝑚𝑑(𝑡)²𝜙_2𝑚𝑑𝜉²=−^𝑁𝑆_𝑚=1^𝑁𝑆_𝑛=1𝑇_1𝑚(𝑡)𝑇_1𝑛(𝑡)𝑑𝜙_1𝑚𝑑𝑑𝜉²𝜙_1𝑛𝑑𝜉². (2.28) This gives 𝑇_2𝑚∑=−^{𝑁𝑠𝑛=1}𝑇_1𝑚𝑇_1𝑛(𝑑𝜙_1𝑚/𝑑𝜉)(𝑑²𝜙_1𝑛/𝑑𝜉²)𝑑²𝜙_2𝑚/𝑑𝜉². (2.29) Therefore, (2.27) becomes𝑢₂(𝜉,𝑡)=^𝑁𝑠_𝑚=1^𝑁𝑠_𝑛=1𝑇_1𝑚𝑇_1𝑛𝜓_2𝑚𝑛(𝜉) (2.30) with𝜓_2𝑚𝑛(𝜉)=−(𝑑𝜙_1𝑚/𝑑𝜉)(𝑑²𝜙_1𝑛/𝑑𝜉²)𝑑²𝜙_2𝑚/𝑑𝜉²𝜙_2𝑚. (2.31) From (2.30), the axial displacement 𝑢₂ at an arbitrary sampling point “𝑖’’ is 𝑢₂(𝜉_𝑖,𝑡)=^𝑁𝑠𝑁_𝑚=1^𝑠_𝑛=1𝑇_1𝑚(𝑡)𝑇_1𝑛(𝑡)𝛿^𝑖_2𝑚𝑛. (2.32) Putting (2.30) into (2.9) and then applying DQM approximation, one has ^𝑁_𝑗=1𝐶^𝑁(2)_𝑖𝑗^𝑠𝑁_𝑚=1^𝑠_𝑛=1𝑇_1𝑚𝑇_1𝑛𝛿^𝑗_2𝑚𝑛=−^𝑁_𝑗=1𝐶^𝑁(1)_𝑖𝑗^𝑠_𝑚=1𝑇_1𝑚𝛿^𝑗𝑁_1𝑚𝑘=1𝐶^𝑁(2)_𝑖𝑘^𝑠_𝑛=1𝑇_1𝑛𝛿^𝑘_1𝑛(𝑖=1,…,𝑁) (2.33) in a matrix form [𝐷₂]⋅{𝛿₂}={𝑓₂}. (2.34) Note that the right-hand-side term {𝑓₂} is the pseudodynamic force vector totally dependent on the 1st-order deflection solution already obtained. After the unknown vector {𝛿₂} composed of 𝛿^𝑗_2𝑚𝑛 is determined from (2.34), 𝑢₂ can be calculated from the relationship 𝜓_2𝑚𝑛(𝜉)=^𝑁_𝑖=1𝛿^𝑖_2𝑚𝑛𝐿_𝑖(𝜉) (2.35) and (2.30).

2.5. Semianalytical Solution of the 3rd-Order Perturbation Equation

The 3rd-order perturbation equation (2.10) can be rewritten as 𝜕⁴𝑤₃𝜕𝜉⁴+𝜕²𝑤₃𝜕𝑡²=𝑓(𝜉,𝑡), (2.36) where the right-hand-side term𝑓(𝜉,𝑡)={[𝜕𝑢₂+1𝜕𝜉2(𝜕𝑤₁)𝜕𝜉²𝜕]}²𝑤₁𝜕𝜉²−𝜎₂𝜕²𝑤₁𝜕𝑡². (2.37)

Substituting (2.17) and (2.30) into (2.37) gives𝑓(𝜉,𝑡)=^𝑁𝑠_𝑙=1^𝑁𝑆_𝑚=1^𝑁𝑠_𝑛=1𝑇_1𝑙𝑇_1𝑚𝑇_1𝑛𝜑_𝑙𝑚𝑛+^𝑁𝑠_𝑛=1𝑝²_𝑛𝜎₂𝑇_1𝑛𝜙_1𝑛, (2.38) where𝜑_𝑙𝑚𝑛(𝜉)={𝑑𝜓_2𝑙𝑚+1𝑑𝜉2[𝑑𝜙_1𝑙𝑑𝜉𝑑𝜙_1𝑚𝑑𝑑𝜉]}²𝜙_1𝑛𝑑𝜉². (2.39) The solution of (2.36) can be readily obtained by the method of separation of variables 𝑤₃(𝜉,𝑡)=^𝑁𝑠_𝑟=1𝑇_3𝑟(𝑡)𝜙_3𝑟(𝜉). (2.40) It is evident from (2.16) and (2.36) that 𝜙_3𝑚(𝜉) is the same as 𝜙_1𝑚(𝜉). Therefore,𝑤₃(𝜉,𝑡)=^𝑁𝑠_𝑟=1𝑇_3𝑟(𝑡)𝜙_1𝑟(𝜉), (2.41) and (2.36) becomes^𝑁𝑠_𝑟=1𝑇_3𝑟𝜕(𝑡)⁴𝜙_1𝑟𝜕𝜉⁴+^𝑁𝑠_𝑟=1𝜕²𝑇_3𝑟𝜕𝑡²𝜙_1𝑟=𝑓(𝜉,𝑟). (2.42) From (2.16)–(2.18), one has 𝜕⁴𝜙_1𝑟𝜕𝜉⁴=𝑃²_𝑟𝜙_1𝑟. (2.43) Expanding 𝑓(𝜉,𝑡) into Fourier series in function 𝜙_1𝑚, (2.42) becomes ^𝑁𝑠_𝑟=1[𝑇^′′_3𝑟+𝑇_3𝑟𝑝²_𝑟]𝜙_1𝑟(𝜉)=^𝑁𝑠_𝑟=1𝑓_𝑟(𝜉,𝑡)𝜙_1𝑟(𝜉) (2.44) which yields 𝑇^′′_3𝑟+𝑇_3𝑟𝑝²_𝑟=𝑓_𝑟(𝜉,𝑡), (2.45) where right-hand-side term 𝑓_𝑟1(𝑡)=𝑁²_3𝑟¹₀𝑓(𝜉,𝑡)𝜙_1𝑟(𝜉)𝑑𝜉=^𝑁𝑠_𝑙=1^𝑁𝑠_𝑚=1^𝑁𝑠_𝑛=1𝑎_{𝑙𝑚𝑛𝑟}𝑇_1𝑙𝑇_1𝑚𝑇_1𝑛+^𝑁𝑠_𝑛=1𝜎₂𝑏_𝑛𝑟𝑝²_𝑛𝑇_1𝑛,𝑁²_3𝑟=¹₀𝜙²_1𝑟(𝜉)𝑑𝜉,𝑎_{𝑙𝑚𝑛𝑟}=1𝑁²_3𝑟¹₀𝜑_𝑙𝑚𝑛(𝜉)𝜙_1𝑟(𝜉)𝑑𝜉,𝑏_𝑛𝑟=1𝑁²_3𝑟¹₀𝜙_1𝑛(𝜉)𝜙_1𝑟(𝜉)𝑑𝜉. (2.46) Since 𝜙_1𝑟 are orthogonal functions, we have𝑏_𝑛𝑟=10,𝑛≠𝑟,𝑁²_3𝑟¹₀𝜙_1𝑟⋅𝜙_1𝑟(𝜉)𝑑𝜉=1,𝑛=𝑟. (2.47) To remove the secular terms in (2.43), 𝑓_3𝑚(𝑡) is expanded as𝑓_𝑟(𝑡)={^𝑁𝑠_𝑙=1^𝑁𝑠_𝑚=1^𝑁𝑠_𝑛=1cos(𝑝_𝑙+𝑝_𝑚+𝑝_𝑛)𝑡𝑎_{𝑙𝑚𝑛𝑟}+^𝑁𝑠_𝑙=1^𝑁𝑠_𝑚=1^𝑁𝑠_{𝑛=1𝑙≠𝑛,𝑛≠𝑚}cos(𝑝_𝑙+𝑝_𝑚−𝑝_𝑛)𝑡𝑎_{𝑙𝑚𝑛𝑟}+^𝑁𝑠_𝑙=1^𝑁𝑠_𝑚=1^𝑁𝑠_{𝑛=1𝑙≠𝑚,𝑚≠𝑛}cos(𝑝_𝑙−𝑝_𝑚+𝑝_𝑛)𝑡𝑎_{𝑙𝑚𝑛𝑟}+^𝑁𝑠_𝑙=1^𝑁𝑠_𝑚=1^𝑁𝑠_{𝑛=1𝑙≠𝑚,𝑙≠𝑛}cos(𝑝_𝑙−𝑝_𝑚−𝑝_𝑛)𝑡𝑎_{𝑙𝑚𝑛𝑟}}+{2^𝑁𝑠_𝑙=1^𝑁𝑠_{𝑚=1𝑙=𝑛}cos𝑝_𝑚𝑡𝑎_{𝑙𝑚𝑙𝑟}+2^𝑁𝑠_𝑙=1^𝑁𝑠_{𝑚=1𝑚=𝑛}cos𝑝_𝑙𝑡𝑎_{𝑙𝑚𝑚𝑟}+2^𝑁𝑠_𝑙=1^𝑁𝑠_{𝑛=1𝑙=𝑚}cos𝑝_𝑛𝑡𝑎_{𝑙𝑙𝑛𝑟}+3^𝑁𝑠_{𝑙=1𝑙=𝑚=𝑛}cos𝑝_𝑙𝑡𝑎_{𝑙𝑙𝑙𝑟}}+𝑝²_𝑟𝜎_2𝑟cos𝑝_𝑟𝑡. (2.48) The following condition must be satisfied in order to remove the last two terms in (2.48), {2^𝑁𝑠_𝑙=1^𝑁𝑠_{𝑚=1𝑙=𝑛}cos𝑝_𝑚𝑡𝑎_{𝑙𝑚𝑙𝑟}+2^𝑁𝑠_𝑙=1^𝑁𝑠_{𝑚=1𝑚=𝑛}cos𝑝_𝑙𝑡𝑎_{𝑙𝑚𝑚𝑟}+2^𝑁𝑠_𝑙=1^𝑁𝑠_{𝑛=1𝑙=𝑚}cos𝑝_𝑛𝑡𝑎_{𝑙𝑙𝑛𝑟}+3^𝑁𝑠_{𝑙=1𝑙=𝑚=𝑛}cos𝑝_𝑙𝑡𝑎_{𝑙𝑙𝑙𝑟}}+𝑝²_𝑟𝜎_2𝑟cos𝑝_𝑟𝑡=0. (2.49) The sum of the coefficients for cos(𝑝_𝑟𝑡) must be zero, which gives the expressions of 𝜎₂ for different values of 𝑚 as 𝜎_2𝑟=−1𝑝²_𝑟{2^𝑁𝑠_𝑚=1[𝑎_{𝑚𝑟𝑚𝑟}+𝑎_{𝑟𝑚𝑚𝑟}+𝑎_{𝑚𝑚𝑟𝑟}]+3𝑎_{𝑟𝑟𝑟𝑟}}. (2.50) Thus the term 𝑓_3𝑚(𝑡) without secular terms is𝑓_𝑟(𝑡)={^𝑁𝑠_𝑙=1^𝑁𝑠_𝑚=1^𝑁𝑠_𝑛=1cos(𝑝_𝑙+𝑝_𝑚+𝑝_𝑛)𝑡𝑎_{𝑙𝑚𝑛𝑟}+^𝑁𝑠_𝑙=1^𝑁𝑠_𝑚=1^𝑁𝑠_{𝑛=1𝑙≠𝑛,𝑛≠𝑚}cos(𝑝_𝑙+𝑝_𝑚−𝑝_𝑛)𝑡𝑎_{𝑙𝑚𝑛𝑟}+^𝑁𝑠_𝑙=1^𝑁𝑠_𝑚=1^𝑁𝑠_{𝑛=1𝑙≠𝑚,𝑚≠𝑛}cos(𝑝_𝑙−𝑝_𝑚+𝑝_𝑛)𝑡𝑎_{𝑙𝑚𝑛𝑟}+^𝑁𝑠_𝑙=1^𝑁𝑠_𝑚=1^𝑁𝑠_{𝑛=1𝑙≠𝑚,𝑙≠𝑛}cos(𝑝_𝑙−𝑝_𝑚−𝑝_𝑛)𝑡⋅𝑎_{𝑙𝑚𝑛𝑟}}. (2.51) Equation (2.45) can then be rewritten as𝑇^′′_3𝑟+𝑇_3𝑟𝑝²_𝑟=𝑓_𝑟(𝑡). (2.52) From (2.13), the initial condition for (2.52) can be derived by using separation of variable method 𝑇_3𝑟(0)=𝑑𝑇_3𝑟(0)𝑑𝑡=0. (2.53) The solution of (2.52) under the given initial condition (2.53) can be obtained through Laplace transform as𝑇_3𝑟(𝑡)=^𝑀_𝑙=1^𝑀_𝑚=1^𝑀_𝑛=1cos(𝑝_𝑟𝑡)−cos((𝑝_𝑙+𝑝_𝑚+𝑝_𝑛)𝑡)(𝑝_𝑙+𝑝_𝑚+𝑝_𝑛)²−𝑝²_𝑟𝛾_{𝑙𝑚𝑛𝑟}+^𝑀_𝑙=1^𝑀_𝑚=1^𝑀_{𝑛=1𝑙≠𝑛,𝑛≠𝑚}cos(𝑝_𝑟𝑡)−cos((𝑝_𝑙+𝑝_𝑚−𝑝_𝑛)𝑡)(𝑝_𝑙+𝑝_𝑚−𝑝_𝑛)²−𝑝²_𝑟𝛾_{𝑙𝑚𝑛𝑟}+^𝑀_𝑙=1^𝑀_𝑚=1^𝑀_{𝑛=1𝑙≠𝑚,𝑚≠𝑛}cos(𝑝_𝑟𝑡)−cos((𝑝_𝑙−𝑝_𝑚+𝑝_𝑛)𝑡)(𝑝_𝑙−𝑝_𝑚+𝑝_𝑛)²−𝑝²_𝑟𝛾_{𝑙𝑚𝑛𝑟}+^𝑀_𝑙=1^𝑀_𝑚=1^𝑀_{𝑛=1𝑙≠𝑚,𝑙≠𝑛}cos(𝑝_𝑟𝑡)−cos((𝑝_𝑙−𝑝_𝑚−𝑝_𝑛)𝑡)(𝑝_𝑙−𝑝_𝑚−𝑝_𝑛)²−𝑝²_𝑟𝛾_{𝑙𝑚𝑛𝑟}(where𝛾_{𝑙𝑚𝑛𝑟}=14𝑎_{𝑙𝑚𝑛𝑟}). (2.54) Finally, we have𝑤₃(𝜉,𝑡)=^𝑁𝑠_𝑟=1𝑇_3𝑟(𝑡)𝜙_1𝑟(𝜉). (2.55) The deflections response of nonlinear free vibration of the beam is then calculated from 𝑤=𝑤₁𝜀+𝑤₃𝜀³+⋯ (2.56) while the nonlinear frequency response, with higher-order terms being neglected, is given by (2.50) and (2.7) as𝑝²_𝑟=1+𝜀²𝜎_2𝑟, (2.57) and the nonlinear frequency response can be obtained from(𝜔_𝑟𝜔_𝐴)²=1−𝜀²1𝑝²_𝑟{2^𝑁𝑠_𝑚=1[𝑎_{𝑚𝑟𝑚𝑟}+𝑎_{𝑟𝑚𝑚𝑟}+𝑎_{𝑚𝑚𝑟𝑟}]+3𝑎_{𝑟𝑟𝑟𝑟}}. (2.58)

3. Numerical Results

To illustrate its efficiency and accuracy, the proposed semianalytical approach is used to study the nonlinear frequency and dynamic response of Euler-Bernoulli beams of rectangular cross section (𝑏×ℎ=0.01m×0.02m) under various boundary conditions with a given initial vibration amplitude and zero-valued initial velocity. The dimensionless vibration amplitude is defined as 𝐴=𝑊/𝜌 (or 𝐴=6𝑊/√3ℎ for a rectangular beam) where W is the vibration amplitude at the midpoint of the beam. The parameters used herein are 𝐿=0.5m, 𝐸=150 GPa, and 𝑚=8×10³kg/m³. The results are presented in both tabular and graphical forms in terms of the normalized frequency ratio 𝜔_𝑛/𝜔_𝐿 and dimensionless dynamic deflections at the midpoint of the beam.

Among the sampling point distribution schemes available, a nonuniform distribution is chosen in the present analysis due to its excellent convergence [21, 26], that is,𝑥_𝑘=12(1−cos𝑘𝜋𝑁−1)(𝑘=0,1,…,𝑁−1). (3.1)

In what follows, the total number of sampling points in the space domain is N = 11, unless stated otherwise.

In Figures 1–3, the present solutions and exact solutions are represented by solid lines and circular dots, respectively.

Figure 1 

Nonlinear fundamental frequency versus dimensionless central amplitude relationship of beams with different boundary conditions.

Figure 2 

Nondimensional central deflection response of a simply supported beam under different initial vibration amplitudes.

Figure 3 

Nondimensional central deflection response of a simply supported beam.

3.1. Nonlinear Vibration Frequency

Tables 1, 2, and 3 list, respectively, the normalized fundamental frequency ratios 𝜔_𝑛/𝜔_𝐿 for a simply supported beam, a clamped beam, and a beam simply supported at one end but clamped at the other end at varying dimensionless vibration amplitudes 𝐴=0.1∼2.0. All beams exhibit typical “hard-spring” behavior, that is, the frequency ratio increases with an increase in vibration amplitude. Comparisons between the results with varying total number of sampling points 𝑁=9,11,13 show that the present method converges very well to produce very accurate results when 𝑁≥11. Our results are compared with the existing results obtained by using elliptic function [1] and the finite element method [28]. The relative errors, defined as(PresentSolution−ExistingSolution)ExistingSolution×100%, (3.2) are also provided. Excellent agreement can be observed in these tables. The relationship between nonlinear fundamental frequency and dimensionless vibration amplitude (𝐴=0.0∼6.0) of these three beams are shown in Figure 1 where curves a, b, and c are for a simply supported beam, a beam simply supported at one end but clamped at the other end, and a clamped beam. The scatter points represent the elliptical function results using the formulae provided by Woinowsky-Krieger [1]. As can be seen, the present results and the elliptical function results are almost identical. The clamped beam has the highest fundamental frequency with the lowest deflection while the simply supported beam has the smallest fundamental frequency and the largest deflection. This is because the end support of the clamped beam is much stronger than its simply supported counterpart.

Table 1 

Normalized frequency ratio 𝜔_𝑛/𝜔_𝐿 for a simply supported beam.

Table 2 

Normalized frequency ratio 𝜔_𝑛/𝜔_𝐿 for a clamped beam.

Table 3 

Normalized frequency ratio 𝜔_𝑛/𝜔_𝐿 for a beam simply supported at one end and clamped at the other end.

3.2. Nonlinear Dynamic Response

Convergence study on the nonlinear dynamic response is conducted in both Tables 4 and 5 where central deflections of a simply supported beam at different time (in sec), calculated by using different total number of sampling points N in the space domain and different number of truncated series terms 𝑁_𝑠 in the time domain, are given and compared with the results calculated by using the formulae given by Woinowsky-Krieger [1]. The dimensionless vibration amplitude is 𝐴=1.0. Error 1, Error 2, and Error 3 in Table 4 indicate the relative error percentages between the exact solutions and the present results with 𝑁=9 and 𝑁_𝑠=1, 𝑁=9 and 𝑁_𝑠=3, 𝑁=9 and 𝑁_𝑠=5 while Error 4, Error 5, and Error 6 in Table 5 are the relative error percentages between the exact solutions and the present results with 𝑁=9 and 𝑁_𝑠=3, 𝑁=11 and 𝑁_𝑠=3, 𝑁=13 and 𝑁_𝑠=3. Our semianalytical solutions converge monotonically as 𝑁_𝑠 increases but nonmonotonically with an increasing N. However, the present method is capable of giving very accurate results with relative error less than 1.5% when 𝑁_𝑠≥3 and 𝑁≥11. We have calculated the nonlinear dynamic response of beams with other end supports as well and found that the results exhibit the same convergence characteristics as observed in Tables 4 and 5. These results are, therefore, not detailed here.

Table 4 

Nonlinear central deflection response of a simply supported beam: solutions with varying number of terms in time domain series 𝑁_𝑠.

Table 5 

Nonlinear central deflection response of a simply supported beam: solutions with varying total number of sampling points 𝑁.

Figure 2 depicts the nonlinear dynamic central deflection response of a simply supported beam with dimensionless initial deflection 𝑤₁₀=0.5sin(𝜋𝜉), 0.8sin(𝜋𝜉),1.0sin(𝜋𝜉), and 1.4sin(𝜋𝜉). The time interval in this example is relatively short, that is, 𝑡=0∼0.1 sec. The result within a much longer time interval 𝑡=0∼1 sec is displayed in Figure 3 where 𝑤₁₀=1.0sin(𝜋𝜉). 𝑁=11 and 𝑁_𝑠=3 are used in these two examples. Again, the present results agree very well with the elliptical function solutions by Woinowsky-Krieger [1]. It is also seen that a bigger initial deflection magnitude results in larger vibration amplitude.

It is worthy of noting that unlike the conventional numerical integration schemes commonly used to calculate dynamic response whose accuracy is quickly degraded at later time steps due to the accumulation of numerical errors, the proposed technique does not involve iterative time integration and accumulative errors and is therefore able to give results with excellent accuracy even for a very long time interval, as demonstrated in Figure 3.

4. Conclusions

A new semianalytical method for geometrically nonlinear vibration analysis of Euler-Bernoulli beams with different boundary conditions is presented in this paper. The proposed method is based on the perturbation technique, differential quadrature approximation, and an analytical series expansion in the time domain. Numerical results show that this method has excellent accuracy and convergence characteristics. Compared with other numerical approaches that use step-by-step time integration, the present method has a unique advantage of being capable of producing results with very good and stable accuracy at a long time interval because the error accumulation is avoided due to the use of the analytical time series.