1 INTRODUCTION

In the design and fabrication of many engineering structures and machines, vibrational analysis and dynamical response of structures like beams (or rods or shafts) are important factor to investigate in order to increase the performance of these structures. For the analysis of complex engineering structures like bridges, tall buildings, vehicle guide-ways, huge cranes, turbines and compressor blades, beams can be used as simple model. The dynamical response of beams is governed by linear as well as nonlinear differential equations both in space and time. To study their behavior it is therefore important to design methods for the numerical solutions of these partial differential equations. In this respect, recently some approximate analytical techniques (see for instance, ^{Barari et al., 2010}; ^{Bayat et al., 2010}; Baghari et al., 2014; ^{Jafari et al., 2014}; ^{Bayat et al., 2011};, ^{Ganji, 2012}; ^{He, 1999}; ^{He, 2006}; ^{Khan et al., 2012}; ^{Liu and Gurram, 2009}; ^{Mirgolbabaei et al., 2010}; ^{Sfahani et al., 2011}) as well as numerical methods by ^{Lai et al. (2008}) (see also, ^{Boukhalfa et al., 2010}; ^{Da Silva et al., 2009}; ^{Ganji et al., 2011}) and references therein, have been proposed to investigate the nonlinear vibrations of beams for designing and fabrication purpose.

The theory of Euler-Bernoulli beams which is based on the assumption that beam curvature is related to its bending moment provides a good explanation of long isotropic bar’s bending behavior. Several methods have been proposed for obtaining solution to the nonlinear equations of motions governing the Euler-Bernoulli beam’s response. Nonlinear vibrations of the Euler-Bernoulli beam has been investigated by ^{Barari et al. (2011}) using the Variational iteration and parameterized perturbation methods. Nikkar et al. (2014) studied the nonlinear vibration response of Euler-Bernoulli beam by approximate analytical techniques where they utilized Variational He’s approach and Laplace iteration technique in order to solve the respective nonlinear governing equations. ^{Pakar and Bayat (2013}) applied a max-min approach, also by He, to investigate the nonlinear vibrational response of Euler-Bernoulli beams subjected to an axial load. Nonlinear behavior of Euler-Bernoulli beam with different end conditions has also been studied by ^{Rafieipour et al. (2014}) using Laplace iteration method. Lindstedt-Poincare techniques for non linear vibrational analysis have been employed to double-clamped and simply-supported beams subjected to an axial load by ^{Ahmadian et al. (2009}). ^{Javanmard et al. (2013}) used He’s Energy method in order to get solution of nonlinear vibration problem of Euler-Bernoulli beam subjected to axial loads. A similar study has been performed earlier by Pirbodaghi et al. (2009) but with homotopy analysis method. The nonlinear equation of motion for large amplitude free vibrational analysis of Euler-Bernoulli beam resting on variable elastic foundation was solved by ^{Mirzabeigy and Madoliat (2016}) with the application of second order homotopy perturbation method. Moreover, the nonlinear vibrational response of Euler-Bernoulli beams with geometric nonlinarity and subjected to axial loads was computed by Johnson et al. S. (2014) where they used a differential transform and two auxiliary parameter based homotopy analysis techniques.

Here, in this article, we present a new approach to solve equations of motions governing the dynamics of Euler-Bernoulli clamped-clamped-beam which is fixed from one end. This article is organized as follows. In Section II, a mathematical framework governing the nonlinear vibrations of Euler-Bernolulli beam is presented. In Section III, continuous Galerkin-Petrov time-discretization method is developed. In Section IV, the developed method afterwards applied to the governing equations of Euler-Bernoulli beam presented in Section II. In Section V, simulations data and numerical results are discussed. In last Section, conclusions are drawn based on the obtained simulation results in Section V.

2 MATHEMATICAL FORMULATION

Let us consider a straight beam composed of homogeneous material having length *L*, placed on an elastic foundation with modulus of elasticity as *E*. The beam is experiencing an axial force *F* as shown in Figure 1, below. Let ‘A’ be the cross-sectional area of beam which is assumed to be uniform throughout its length. Further, suppose that there is no in plane deformation i.e. the planes of the cross section remains plane after deformation. Now by ignoring the transverse normal and shear strains one can write the equation of motion of the Euler-Bernoulli beam including the mid-plane stretching effect as follows

where, *d* is the coefficient of viscous damping, *K* is the stiffness of foundation and *q*(*x, t*) is the transverse directional load. In the absence of non-conservative forces above equation (1) reduced to

In order to undimensionalize the above equation we choose the following variables

where, _{
rG
} is the gyration radius of the cross-section of the beam. Using equation (3) and omitting the tilde sign from the variables equations (2) thus can be written as follows

Now, by assuming a product solution *Y* (*x, t*) = *φ*(*x*) *η*(*t*) where *φ*(*x*) is the beam’s eigen mode, and by applying the Galerkin method, an equation governing the beam dynamics is obtained as

Where ε and λ can be determined by the relations

Moreover, in order to complete the above problem description, the beam is subjected to following initial conditions

where *A* is the amplitude of oscillation. Equation (5) along with initial conditions in Equation (6) governs the nonlinear vibrational response of Euler-Bernoulli beams.

3 DESCRIPTION OF CONTINOUS GALERKIN-PETROV TIME-DISCRETIZATION METHOD

In this section, we present a time discretization method to solve nonlinear equations governing the dynamics of Euler-Bernoulli beams. It is easy to see that Equation (5) can be transformed into a system of nonlinear differential equations of the following form

Where, *U*(*t*) = (*η,*
*ξ*)^{T} and *f* denotes the vector of unknown-state variables and vector comprising non-linear functions in both time and space, respectively. Here, *ξ* is the first order time derivative of *η.* In order to develop a numerical method to solve such dynamical systems, here, we use the concept of a Galerkin-type formulation and time-discretization. Let *S* be a space of all possible solutions to set of nonlinear differential equations in Equation (7). We search a vector of unknown states *U*: [0, *τ*] → *S* such that

By choosing ℑ as a test space the problem in Equation (8) can be formulated as follows: Find *U* ∈ *S* such that

Equation (9) is the weak form of problem in Equation (8). Now, to solve it over the time interval *I* = (0, *τ*) let us discretized the interval *I* into further sub-intervals with each sub-interval _{
In
} = (_{
tn
}
_{-1}, _{
tn
} ) where *n* belongs to the set {1, ...,N} Now, each interval is having the following property

Moreover, it also satisfy

In order to compute the vector of state variables *U*(*t*) at time interval _{
In
} , we choose Ψ_{
n,j
} (*t*) ∈ P_{
k
} (_{
In
} ) as a set of basis functions. Here, **P**
_{
k
} (_{
In
} ) represents the *k*-th degree polynomials space defined on interval _{
In
} by

The vector of state variables *U*(*t*) in Equation (5) is thus approximated by

Where,
_{
SH
} Global vector of state variables *U*
_{σ}(*t*): _{
In
} → _{
SH
} In is in the discrete-solution space
*S* where,

with _{
Īn
}
*=*
_{
In
} ⋃ {_{
tn
}
_{-1}, _{
tn
} }. The symbol σ above represents a discretization parameter. The test function *V*(*t*) is taken from discrete test space

By denoting *V*
_{σ}(*t*) as a discrete approximation of the test function *V*(*t*), then the weak problem in Equation (9) takes the following Variational form which states:

Find

Above
_{
φn,i
} ∈ **P**
_{k-1}(_{
Īn
} ) on the time interval *I* zero everywhere in the set *I*, except at the time interval _{
Īn
} = [_{
tn
}
_{-1}, _{
tn
} ], and **v** ∊ _{
SH
} as an arbitrary scalar function then the variational problem on time interval _{
In
} in Equation (12) reads

Or, equivalently, we can write by using Equation (11)

To define Basis functions in Equation (14), a mapping
_{
Īn
} . Now, for every
*t* ∈ _{
Īn
} such that

The basis functions defined on physical time interval _{
Īn
} can now be calculated in natural time interval as

Where,

_{
δk,j
} Denote the Kronecker’s delta. To get a second-order approximation of discrete vector of state variables *U*
_{σ}(*t*) we require to find the coefficients
*j* ∈{0, 1, 2}. By using property of basis function from Equation (15) and with the application of initial condition the coefficients

Since, it is easy to deal with natural-time interval *Î*, thus Equation (14) is reconstructed into

where, *i* = 0, ..., *k* - 1 and

Now, with the application of (*k* + 1)-point Gauss-Lobatto formula in Equation (18) we write

Where, Ŵ_{
l
} represents the weights
**P**
_{k-1}(*t*). Using the abbreviation

Hence, we finally arrived at the following system of coupled equations

with initial conditions as described in Equation (16).

**3.1 Computation of Parameters α**
_{i,j}
**, β**
_{i,j}
**and η**
_{j,l}

To compute unknown parameters in Equation (21) let us select the following set of test functions

Where, each
_{
αi,j
} , _{
βi,j
} and _{
ηj,l
} are calculated and are given as under

4 APPLICATION TO GOVERNING EQUATIONS OF EULER-BERNOULLI BEAMS

The method presented in Section 3 is applied to motion equation in Equation (5) and the obtained discrete set of equations are solved using numerical algorithm outlined as under

**
Numerical Algorithm:
**

**
Step 1:
** Initialize the known parameters ε and

*λ*

**
Step 2:
** Choose a time step size σ

_{ n }

**
Step 3:
** Set the initial conditions at time

_{ tn }as

**
Step 4:
** First solve the following equations for unknowns

**
Step 5:
** Using the computed values for

*Step 4*, Solve the following nonlinear equations for

**
Step 6:
** Use the computed values of

*Step 5*to get the numerical values of

**
Step 7:
** Compute the solution at time step

_{ tn }by equation (11)

**
Step 8:
** Update the solution for the next time iteration by

**
Step 9:
** Update the time step

_{ tn }←

_{ tn }+ σ

_{ n }

**
Step 10:
** Go to Step 4.

5 RESULTS AND DISCUSSIONS

In this Section, nonlinear vibrational response of a clamped-clamped Euler-Bernoulli beam fixed at one end is presented by a new numerical method. The results computed by the presented method are compared with the other existing methods in literature (for instance, with the results of Barai et al. 2011, Nikkar et al. 2014 and Johnson et al. 2014). For the numerical computations a time step size of 0.05 is chosen for the discretization of time domain by the present method. The non-dimensional deflection *η*(*t*) is calculated by the presented method at the center of the beam and is shown in Figure 2, with the variation of non-dimensional maximum amplitude of oscillation over a wide range of time.

In Table 1 and Table 2, a comparison is made between the presented method and four different methods, with semi-analytical approach by ^{Barari et al. (2011}), with a semi-analytical method by Nikkar et al. (2014), with semi-analytical approach by Johnson et al. (2014) and with numerical method RK-4. It is seen that the computed values of deflection at the center of the beam for different values of maximum amplitude of oscillation are in strong agreement with the results by these methods. In Table 3, relative errors are provided which are computed by the formula

Where in, Error-1 is calculated by taking RK-4 solution as an existing solution. In Error-II the existing solution is taken from the semi-analytical approach by ^{Barari et al. (2011}). Error-III is obtained by taking the existing solution form the method by Nikkar et al. (2014). In all of these calculations the errors are shown for three different amplitude of oscillation of beam over a wide range of time. From Table 3, it can be clearly seen that the results by present method converges closely to other results in literature and are in strong agreement.

A | RK-4 |
Barari et al. (2011) |
Nikkar et al. (2014) | Present |
---|---|---|---|---|

0.01 | 0.008775827082332 | 0.008775710513259 | 0.008775710509477 | 0.008653113482869 |

0.1 | 0.087759719142174 | 0.087643206097855 | 0.087642827100616 | 0.086406655501930 |

0.2 | 0.175528204796203 | 0.174597455961528 | 0.174585253116176 | 0.172060842262657 |

0.3 | 0.263314167649526 | 0.260180487515907 | 0.260086892907199 | 0.263974531478741 |

0.4 | 0.351126207698248 | 0.343723299721559 | 0.343323578255470 | 0.349072483074928 |

0.5 | 0.438972761263577 | 0.424576467617775 | 0.423336627126820 | 0.417180214867134 |

0.6 | 0.526862050755839 | 0.502116086386647 | 0.498973411668122 | 0.492625565690737 |

0.7 | 0.614802037962082 | 0.575749233260282 | 0.568819355464365 | 0.563875657455006 |

0.8 | 0.702800381522887 | 0.644918846456226 | 0.631123805468134 | 0.630358772839251 |

0.9 | 0.790864399137227 | 0.709107932868351 | 0.683724780583061 | 0.691556439573670 |

1.0 | 0.879001034896931 | 0.767843030699686 | 0.723980197798699 | 0.766656852847480 |

A | |η_{
RK-4
} -η_{
present
} | |
|η_{
Barari et al. (2011
}_{
)
} -η_{
present
} | |
|η_{
Nikkar et al. (2014)
} -η_{
present
} | |
|η_{
Johnson et al. (2014)
} -η_{
present
} | |
---|---|---|---|---|

0.01 | 0.000122713599463 | 0.0001225970303900 | 0.00012259702660799 | 0.0003821333327880 |

0.1 | 0.001353063640244 | 0.0012365505959250 | 0.00123617159868600 | 0.0019824623866750 |

0.2 | 0.003467362533546 | 0.0025366136988710 | 0.00252441085351901 | 0.0008865790820170 |

0.3 | 0.000660363829215 | 0.0037940439628340 | 0.00388763857154201 | 0.0052723348189849 |

0.4 | 0.002053724623320 | 0.0053491833533690 | 0.00574890481945800 | 0.0086853097383650 |

0.5 | 0.021792546396443 | 0.0073962527506410 | 0.00615641225968599 | 0.0087110167343409 |

0.6 | 0.034236485065102 | 0.0094905206959099 | 0.00634784597738497 | 0.0002395324508129 |

0.7 | 0.050926380507076 | 0.0118735758052760 | 0.00494369800935901 | 0.0063966672702809 |

0.8 | 0.072441608683636 | 0.0145600736169750 | 0.00076503262888294 | 0.0161331435918892 |

0.9 | 0.099307959563557 | 0.0175514932946811 | 0.00783165899060900 | 0.0078649950295989 |

1.0 | 0.112344182049451 | 0.0011861778522060 | 0.04267665504878090 | 0.0153845773856570 |

A | time | Error-I | Error-II | Error-III |
---|---|---|---|---|

0.01 | 5 | 8.39406344436937 | 1.43006947021705 | 1.43006930029566 |

10 | 1.64295821637604 | 1.65335398475738 | 0.46126609193088 | |

25 | -0.40714155296337 | -0.30013048066676 | -0.30013048391456 | |

50 | 0.34564381067195 | 0.70609180050204 | 0.70609178786218 | |

100 | 0.89942811653737 | 1.32438720734438 | 1.32438716831965 | |

0.50 | 5 | -0.01611628030638 | -2.71860640853757 | -3.34003971902284 |

10 | -2.19636228576523 | 0.94834463220716 | -0.15087284130191 | |

25 | 1.33999271743426 | 0.18473004818300 | -0.63051410694749 | |

50 | 4.46407178480245 | -0.21897390156171 | -1.08460097418586 | |

100 | -5.91449682727209 | -0.55016306148146 | -0.57899624097327 | |

1.00 | 5 | -5.825758083900410 | 0.33425814206559 | -1.25501908713661 |

10 | -0.411754581578781 | -1.42949464499397 | -0.67810815791013 | |

25 | -0.614810964667922 | 3.07084443826624 | 0.52366254911929 | |

50 | 0.200159632081500 | -0.86954946868848 | 0.65475889510081 | |

100 | 1.993160046364360 | -0.27554338389560 | 0.96750105622486 |

In Figure 3, time marching deflection response over a wide range of time is shown where the result obtained from present method is graphically compared with the result computed from semi-analytical approach by Nikkar et al. (2014). In Figure 4, the time marching response of deflection by present method is compared with the solution of deflection response of center of the beam by ^{Barari et al. (2011}). It can be observed from these two Figures that the center of the beam shows same deflection behavior for a large time period as depicted by the semi-analytical approaches in literature. A parametric study is carried out to observe time marching response of the central deflection of beam. In Figure 5 and Figure 6, the central deflection of beam is observed with varying the parameters *(* and *λ* The initial maximum amplitude of vibration was chosen to be 1.0 and a time step size of 0.005 for the calculation by the presented method. By varying the parameters *(* and *λ* an increase in the maximum amplitude of vibration is seen as the time marches. In Table 4, a comparison is made between the accumulated errors by classical numerical scheme RK-4 and present method. Also the accuracy of solution by both of the numerical methods is shown for two different time step sizes. In all the computations of accumulated errors in Table 4 the reference solution by Johnson et al. (2014) is taken into consideration. It is observed that our present numerical method is more accurate in predicting the solution over the classical RK-4 method. This is highly due to the fact that present method does not involve iterative time integration and therefore have less accumulated errors as shown in Table 4. Moreover, it is observed that that the accumulated error decreases with an increase in the time step size for both the classical RK-4 and present method.

Time- step size | Number of time steps | Accumulated error by RK-4 | Accuracy (in percent) by RK-4 | Accumulated error by present method | Accuracy (in percent) by present method |
---|---|---|---|---|---|

0.05 | 5 | 7.1629284856E-04 | 7.838643342877E+01 | 4.7825078312E-04 | 8.556916328361E+01 |

10 | 4.1601192713E-04 | 9.991054140377E+01 | 2.8166615852E-04 | 9.996668334601E+01 | |

25 | 8.9254489341E-05 | 9.910446043034E+01 | 6.3862965951E-05 | 9.935922760337E+01 | |

50 | 1.8633928689E-04 | 9.809327174198E+01 | 1.4120960638E-04 | 9.855506398418E+01 | |

100 | 3.6919714161E-04 | 9.584366665341E+01 | 2.9192878467E-04 | 9.671353538314E+01 | |

250 | 7.2383261971E-04 | 7.572278868091E+01 | 5.9430126403E-04 | 8.006724623713E+01 | |

500 | 3.7148542418E-04 | 9.975956260686E+01 | 1.8132325720E-04 | 1.000734824724E+02 | |

0.005 | 5 | 5.04555298800E-05 | 9.82223873076E+01 | 2.39457037950E-05 | 9.91563623037E+01 |

10 | 2.64336613610E-05 | 1.00163169326E+02 | 1.36876932130E-05 | 1.00315112462E+02 | |

25 | 4.97951263280E-05 | 1.00561378909E+02 | 2.28701839830E-05 | 1.01222287224E+02 | |

50 | 1.70105932140E-05 | 9.98284057000E+01 | 3.19282474000E-06 | 9.99677923915E+01 | |

100 | 3.59407763030E-05 | 9.96277324234E+01 | 6.42406862100E-06 | 9.99334607456E+01 | |

250 | 1.61756733776E-04 | 9.35316240962E+01 | 2.41848699480E-05 | 9.90328882986E+01 | |

500 | 4.74408476010E-05 | 9.94577701222E+01 | 1.21912724690E-05 | 1.00139341359E+02 |

6 CONCLUSIONS

A new numerical method is presented for nonlinear vibrational analysis of Euler-Bernoulli beams. The non-dimensional central deflection response of a clamped-clamped-Euler-Bernoulli beam (fixed from one end side) is investigated under different initial vibration amplitudes using this numerical technique. A comparison of computed results by present method is made with the results by other semi-analytical schemes in literature and the numerical values of central deflection of beam are tabulated. Relative errors between the present method and other semi-analytical approaches in literature have been shown to observe the convergence of solution obtained by present method. Moreover, a parametric study is carried out to investigate effect of different parameters on beam’s central deflection. It is worth mentioning that the obtained results by presented method are in strong agreement with results by semi-analytical approaches (for instance, ^{Barari et al. 2011}, Nikkar et al. 2014 and Johnson et al. 2014) in literature and also with classical numerical scheme RK-4. It is worthy of noting that the present scheme does not involve iterative time integration like other methods in literature and thus no associated accumulative errors. Therefore, it is capable to provide accurate results with high accuracy and convergence characteristics as can be seen by the presented results and thus can be utilized to compute dynamical response of structures like rods and shafts. We conclude that the present method has great potential in dealing with nonlinear vibrational response of beams and similar structures like rods and shafts.