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

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.


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 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.

Perturbation Equations
Consider an isotropic, homogeneous slender beam of length L with uniform cross-section area subjected to a dynamic transverse load q x, t .Let u and w be the displacements parallel to the x-and y-axes, t the time, and m 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 where E is the elastic modulus of the beam, and A and I 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: where M x 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: where ω L 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

2.4
Deflection w, axial displacement u, and frequency parameter p can be expanded into series forms in terms of a small perturbation parameter ε w w Inserting 2.5 -2.7 into 2.4 , and equating the terms of the same power of ε yield a series of perturbation equations The initial conditions considered in the present study are where w 10 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.

Differential Quadrature Method
According to the DQM rule, an unknown function g x and its nth 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 x-axis as where N is the total number of sampling points, L i x are Lagrange interpolation polynomials, and the weighting coefficients C n ij can be calculated from recursive formulae 20, 21 .

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 F 0. Hence, the first-order perturbation becomes The solution of 2.16 can be readily obtained by separation of variables where N s is the number of terms in the truncated time series, is the solution in the time domain for the first-order perturbation equation, and p m and δ 1m ξ are the frequency parameter and the associated mode shape function.
The transient deflection at an arbitrary sampling point "i" can be expressed as a cosine series with a supplementary term as where δ i 1m φ 1m ξ i needs to be determined.The frequency parameter p m of beams with different boundary conditions can be written in a general form as in which j 0 for a simply supported beam, j 1/2 for a clamped beam, j −1/2 for a beam clamped at one end and free at the other end, and j 1/4 for a beam clamped at one end but simply supported at the other end.
Since w 1 ξ i , t needs to satisfy the initial conditions in 2.12 , one has

2.21
The time domain function in 2.19 can be readily obtained as

2.23
The deflection response of the 1st-order perturbation equation can be determined by 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

2.25
This linear equation system can be written in a matrix form as where D 1 and {δ 1 } are the coefficient matrix and unknown vector to be solved, and {f 1 } is the generalized load vector due to initial deflection.Once {δ 1 } 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 .

Semianalytical Solution of the 2nd-Order Perturbation Equation
To solve the 2nd-order perturbation equation 2.9 , the unknown axial displacement u 2 is expressed as Substituting the 1st-order deflection solution and 2.27 into 2.9 yields

2.28
This gives Therefore, 2.27 becomes with From 2.30 , the axial displacement u 2 at an arbitrary sampling point "i" is Putting 2.30 into 2.9 and then applying DQM approximation, one has

2.34
Note that the right-hand-side term {f 2 } is the pseudodynamic force vector totally dependent on the 1st-order deflection solution already obtained.After the unknown vector {δ 2 } composed of δ j 2mn is determined from 2.34 , u 2 can be calculated from the relationship and 2.30 .

Semianalytical Solution of the 3rd-Order Perturbation Equation
The 3rd-order perturbation equation 2.10 can be rewritten as where the right-hand-side term where right-hand-side term

2.46
Since φ 1r are orthogonal functions, we have

2.47
To remove the secular terms in 2.43 , f 3m t is expanded as l m n cos p l ta lllr p 2 r σ 2r cos p r t.

Mathematical Problems in Engineering
The following condition must be satisfied in order to remove the last two terms in 2.48 , l m n cos p l ta lllr p 2 r σ 2r cos p r t 0.

2.49
The sum of the coefficients for cos p r t must be zero, which gives the expressions of σ 2 for different values of m as a mrmr a rmmr a mmrr 3a rrrr .

2.50
Thus the term f 3m t without secular terms is l / m,l / n cos p l − p m − p n t • a lmnr .

2.51
Equation 2.45 can then be rewritten as

2.52
From 2.13 , the initial condition for 2.52 can be derived by using separation of variable method T 3r 0 dT 3r 0 dt 0.

2.53
The solution of 2.52 under the given initial condition 2.53 can be obtained through Laplace transform as where γ lmnr 1 4 a lmnr .

2.55
The deflections response of nonlinear free vibration of the beam is then calculated from while the nonlinear frequency response, with higher-order terms being neglected, is given by 2.50 and 2.7 as and the nonlinear frequency response can be obtained from a mrmr a rmmr a mmrr 3a rrrr . 2.58

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 b × h 0.01 m × 0.02 m under various boundary conditions with a given initial vibration amplitude and zero-valued initial velocity.The dimensionless vibration amplitude is defined as A W/ρ or A 6W/ √ 3h for a rectangular beam where W is the vibration amplitude at the midpoint of the beam.The parameters used herein are L 0.5 m, E 150 GPa, and m 8 × 10 3 kg/m 3 .The results are presented in both tabular and graphical forms in terms of the normalized frequency ratio ω n /ω L 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, 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.

Nonlinear Vibration Frequency
Tables 1, 2, and 3 list, respectively, the normalized fundamental frequency ratios ω n /ω L 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 A 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 Mathematical Problems in Engineering of sampling points N 9, 11, 13 show that the present method converges very well to produce very accurate results when N ≥ 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 Present Solution − Existing Solution Existing Solution × 100%, 3.2 are also provided.Excellent agreement can be observed in these tables.The relationship between nonlinear fundamental frequency and dimensionless vibration amplitude A 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.

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 N s in the time domain, are given and compared with the results calculated by using the formulae given by Woinowsky-Krieger 1 .The dimensionless  Our semianalytical solutions converge monotonically as N s 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 N s ≥ 3 and N ≥ 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.
Figure 2 depicts the nonlinear dynamic central deflection response of a simply supported beam with dimensionless initial deflection w 10 0.5 sin πξ , 0.8 sin πξ , 1.0 sin πξ , and 1.4 sin πξ .The time interval in this example is relatively short, that is, t 0 ∼ 0.1 sec.The result within a much longer time interval t 0 ∼ 1 sec is displayed in Figure 3 where w 10 1.0 sin πξ .N 11 and N s 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.

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.

1 l 1 m 1 n 1 l 1 m 1 n 1 l 2 l 1 m 1 l n cos p m ta lmlr 2 l 1 m 1 m n cos p l ta lmmr 2 l 1 n 1 l m cos p n ta llnr 3 l 1
/ n,n / m cos p l p m − p n ta lmnr Ns l Ns Ns / m,m / n cos p l − p m p n ta lmnr Ns l Ns Ns / m,l / n cos p l − p m − p n ta lmnr Ns Ns Ns Ns Ns Ns Ns

1 l n cos p m ta lmlr 2 l 1 m 1 m n cos p l ta lmmr 2 l 1 n 1 l m cos p n ta llnr 3 l 1
Ns NsNs Ns Ns

1 l 1 m 1 n 1 l
/ n,n / m cos p l p m − p n ta lmnr Ns l Ns Ns / m,m / n cos p l − p m p n ta lmnr

M n 1 cos p r t − cos p l p m p n t p l p m p n 2 1 l 1 l 1 l
/ n,n / m cos p r t − cos p l p m − p n t p l p m − p n 2 / m,m / n cos p r t − cos p l − p m p n t p l − p m p n 2 / m,l / n cos p r t − cos p l − p m − p n t p l − p m − p n 2 − p 2 r γ lmnr

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

Figure 2 :Figure 3 :
Figure 2: Nondimensional central deflection response of a simply supported beam under different initial vibration amplitudes.
It is evident from 2.16 and 2.36 that φ 3m ξ is the same as φ 1m ξ .Therefore,

Table 1 :
Normalized frequency ratio ω n /ω L for a simply supported beam.

Table 2 :
Normalized frequency ratio ω n /ω L for a clamped beam.

Table 3 :
Normalized frequency ratio ω n /ω L for a beam simply supported at one end and clamped at the other end.

Table 4 :
Nonlinear central deflection response of a simply supported beam: solutions with varying number of terms in time domain series N s .Error 2, and Error 3 in Table4indicate the relative error percentages between the exact solutions and the present results with N 9 and N s 1, N 9 and N s 3, N 9 and N s 5 while Error 4, Error 5, and Error 6 in Table5are the relative error percentages between the exact solutions and the present results with N 9 and

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