Nonlinear Vibration Analysis of Euler-Bernoulli Beams by Using Continuous Galerkin-Petrov Time-Discretization Method

Abstract

In this paper, we present a new numerical method for nonlinear vibrational analysis of Euler-Bernoulli beams. Our approach is based on the continuous Galerkin-Petrov time discretization method. The Euler-Bernoulli beam equation which governs its vibrations is transformed into set of ordinary differential equations and the presented method is employed in order to investigate the vibrational response. A comparison is made between present method and different other methods available in literature. It is observed that the obtained results are in strong agreement with other results in literature. We conclude that the present method has a great potential to deal with nonlinear vibration analysis problems of beams and related structures like rods and shafts.

Keywords:
Nonlinear vibration; Euler-Bernoulli beams; Time discretization; Numerical method

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

$EI\frac{{\partial }^{4}Y}{\partial x^4}+\mu \frac{{\partial }^{2}Y}{\partial t^2}+F\frac{{\partial }^{2}Y}{\partial x^2}+d\frac{\partial Y}{\partial t}+KY-\frac{EA}{2L}\frac{{\partial }^{2}Y}{\partial x^2}{\displaystyle \underset{0}{\overset{L}{\int }}{\left(\frac{\partial Y}{\partial x}\right)}^{2}dx}=q(x,t),$ (1)

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

$EI\frac{{\partial }^{4}Y}{\partial x^4}+\mu \frac{{\partial }^{2}Y}{\partial t^2}+F\frac{{\partial }^{2}Y}{\partial x^2}+KY-\frac{EA}{2L}\frac{{\partial }^{2}Y}{\partial x^2}{\displaystyle \underset{0}{\overset{L}{\int }}{\left(\frac{\partial Y}{\partial x}\right)}^{2}dx}=0.$ (2)

Figure 1
A schematic representation of an Euler-Bernoulli beam fixed at one end and subjected to axial load.

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

$\tilde{X}=\frac{X}{L},\tilde{Y}=\frac{Y}{r_G},\tilde{t}=t\sqrt{\frac{EI}{\mu L^4}},\tilde{F}=\frac{F{L}^2}{EI}\text{and}\tilde{K}=\frac{K{L}^4}{EI},$ (3)

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

$\frac{{\partial }^{4}Y}{\partial x^4}+\frac{{\partial }^{2}Y}{\partial t^2}+F\frac{{\partial }^{2}Y}{\partial x^2}+KY-\frac{1}{2}\frac{{\partial }^{2}Y}{\partial x^2}{\displaystyle \underset{0}{\overset{1}{\int }}{\left(\frac{\partial Y}{\partial x}\right)}^{2}dx}=0.$ (4)

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

$\frac{d^2\eta (t)}{d{t}^2}+\epsilon \eta (t)+\lambda {\eta }^3(t)=0$ (5)

Where ε and λ can be determined by the relations

$\epsilon =K+\frac{{\displaystyle \underset{0}{\overset{1}{\int }}{\varphi }^{(iv)}\varphi dx}}{{\displaystyle \underset{0}{\overset{1}{\int }}{\varphi }^{2}dx}}+F\frac{{\displaystyle \underset{0}{\overset{1}{\int }}{\varphi }^{″}\varphi dx}}{{\displaystyle \underset{0}{\overset{1}{\int }}{\varphi }^{2}dx}},\begin{array}{c}\end{array}and\begin{array}{c}\end{array}\lambda =\frac{-{\displaystyle \underset{0}{\overset{1}{\int }}\left({\varphi }^{″}{\displaystyle \underset{0}{\overset{1}{\int }}{{\varphi }^{\prime }}^{2}dx}\right)\varphi dx}}{2.{\displaystyle \underset{0}{\overset{1}{\int }}{\varphi }^{2}dx}}.$

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

$\eta (0)=A,\frac{d\eta }{dt}(0)=\mathrm{0,}$ (6)

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

$\frac{d}{dt}\left\{U(t)\right\}=f(t,U(t)).$ (7)

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

$\begin{array}{l}\frac{d}{dt}\left\{U(t)\right\}=f(t,U(t)),\forall t\in (\mathrm{0,}\tau ),\begin{array}{c}\end{array}where\\ U(0)=U_0.\end{array}$ (8)

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

${\displaystyle \underset{0}{\overset{\tau }{\int }}{V}^T(t)\frac{d}{dt}U(t)dt}={\displaystyle \underset{0}{\overset{\tau }{\int }}{V}^T\left(t\right)f\left(t,U\left(t\right)\right)dt,\forall V\in }\Im .\begin{array}{cc}& \end{array}With\begin{array}{c}\end{array}U(0)=U_0,$ (9)

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

${{\displaystyle \underset{n=1}{\overset{N}{\cup }}\left(t_{n-\mathrm{1,}}{t}_n\right)=\left(\mathrm{0,}\tau \right).}}_{}$

Moreover, it also satisfy

$\left(t_{n-\mathrm{1,}}{t}_n\right){\displaystyle \cap \left(t_{m-\mathrm{1,}}{t}_m\right)=\left\{0\right\},\text{for all}n}\ne mn,m=\text{1}\mathrm{,...,}N.$

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

${\rm P}\left({{\rm I}}_n,S_H\right)\text{:}=\left\{U\left(t\right):{{\rm I}}_n\to S_H\begin{array}{c};\begin{array}{c}\end{array}U\left(t\right)={\displaystyle \sum _{j=0}^k{u}^j{t}^j,\begin{array}{c}\end{array}\text{for all}t\in {{\rm I}}_n,u^j\in S_H}\end{array}\right\}.$ (10)

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

$U\left(t\right)\approx U_{\sigma }\text{(}t\text{)}={\displaystyle \sum _{j=0}^k{u}_n^j{\psi }_{n,j}\left(t\right)\text{, for all t}\in I_n}$ (11)

Where, $u_n^j$ are vector elements in Hilbert space _^SH Global vector of state variables U _σ(t): _^In → _^SH In is in the discrete-solution space $S_{\sigma }^k$ ⊂ S where,

$S_{\sigma }^k\text{:}=\{U\left(t\right)\in \text{C}\left({{\rm I}}_n,S_H\right)\text{:}U\left(t\right){\text{|}}_{{{\rm I}}_n}\in {{\rm P}}_{\text{k}}\left({\overline{{\rm I}}}_n\text{,}{S}_H\right);\text{for all}n=\text{1}\mathrm{,...,}N\},$

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

${\wp }_{\sigma }^k\text{:}=\{V\left(t\right)\in {\text{L}}^{\text{2}}\left({{\rm I}}_n,S_H\right)\text{:}V\left(t\right){|}_{{{\rm I}}_n}\in {{\rm P}}_{\text{k-1}}\left({{\rm I}}_n,S_H\right)\text{;}\forall n=\text{1}\mathrm{,...,}N\}.$

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 $U_{\sigma }\left(t\right)\in U_o+S_{\sigma \text{,0}}^k$ such that

${\displaystyle \underset{0}{\overset{\tau }{\int }}{V}_{\sigma }^{\text{T}}}\left(t\right)\frac{\text{d}}{\text{dt}}{U}_{\sigma }\left(t\right)dt={\displaystyle \underset{\text{0}}{\overset{\tau }{\int }}{V}_{\sigma }^{\text{T}}\left(t\right)}f\left(t\text{,}{U}_{\sigma }\left(t\right)\right)dt\text{;}{V}_{\sigma }\in {\wp }_{\sigma }^k.$ (12)

Above $s_{\sigma \mathrm{,0}}^k\text{:}=S_{\sigma }^k{\displaystyle \cap S_o}$ denotes the subspace of $S_{\sigma }^k$ with zero as initial condition. Define a polynomial_^φ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

${\displaystyle {\int }_{{{\rm I}}_n}\left(v^{\text{T}}\left(t\right)\frac{d}{dt}{U}_{\sigma }\left(t\right)\right)}{\varphi }_{n,i}\left(t\right)dt={\displaystyle {\int }_{{{\rm I}}_n}\left(v^{\text{T}}\left(t\right)f\left(t\text{,}{U}_{\sigma }\left(t\right)\right)\right)}{\varphi }_{n.i}\left(t\right)dt\text{;}v\in S_H.$ (13)

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

${\displaystyle \underset{I_n}{\int }\left(v^{\text{T}}\left(t\right){\displaystyle \sum _{\text{j}=\text{0}}^k{u}_n^j\frac{d}{dt}{\psi }_{n,j}\left(t\right)}\right)}\begin{array}{c}\end{array}{\varphi }_{n,i}\left(t\right)dt={\displaystyle \underset{I_n}{\int }\left(v^{\text{T}}\left(t\right)f\left(t\text{,}{\displaystyle \sum _{\text{j}=\text{0}}^k{u}_n^j{\psi }_{n,j}\left(t\right)}\right)\right)}\begin{array}{c}\end{array}{\varphi }_{n.i}\left(t\right)dt\text{;}v\in S_H.$ (14)

To define Basis functions in Equation (14), a mapping $M_n:\stackrel{⌢}{I}\to {\overline{{\rm I}}}_n$ is used which maps the time interval $\stackrel{⌢}{I}=[-\mathrm{1,1}]$ from natural time interval $\stackrel{⌢}{I}$ to physical time interval _^Īn . Now, for every $\stackrel{⌢}{t}\in \stackrel{⌢}{I}$ there exists t ∈ _^Īn such that

$t=M_n\left(\widehat{t}\right)\text{:}=\frac{t_{n\text{-1}}+t_n}{\text{2}}+\frac{{\sigma }_n}{2}\widehat{t}\in {\overline{{\rm I}}}_n\text{;}n\in \left\{\text{1}\mathrm{,...,}N\right\}.$

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

$\begin{array}{l}{\psi }_{n,j}\left(t\right)\text{:}={\stackrel{⌢}{\psi }}_j\left(M_n^{-1}\left(t\right)\right)\text{;}j=\text{0}\mathrm{,...,}k,\\ {\phi }_{n.i}\left(t\right)\text{:}={\widehat{\phi }}_i\left(M_n^{-1}\left(t\right)\right)\text{;}i=\text{0}\mathrm{,...}k\text{-1.}\end{array}$

Where, ${\stackrel{⌢}{\psi }}_j\in {{\rm P}}_{\text{k}}\left(\widehat{{\rm I}}\right)\text{,}$ and ${\widehat{\phi }}_i\in {{\rm P}}_{\text{k-1}}\left(\widehat{{\rm I}}\right)$ are satisfying the property

${\stackrel{⌢}{\psi }}_j\left(-1\right)={\delta }_{\text{0,}j}\text{,}{\widehat{\psi }}_j\left(\text{1}\right)={\delta }_{k,j}\text{,}$ (15)

_^δ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 $u_n^j$; j ∈{0, 1, 2}. By using property of basis function from Equation (15) and with the application of initial condition the coefficients $u_n^0$ can be determined from the following

$u_n^0=\{\begin{array}{l}{U}_{\text{0}},n=\text{1}\\ u_{n-1}^k\text{,}n>\text{1.}\end{array}$ (16)

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

${\displaystyle \sum _{j=0}^k{\alpha }_{i,j}\left(v^T{u}_n^j\right)=\frac{{\sigma }_{\text{n}}}{\text{2}}{\displaystyle \underset{\widehat{{\rm I}}}{\int }{v}^{T}f}}\left(M_n\left(\widehat{t}\right)\text{,}{\displaystyle \sum _{j=\text{0}}^k{u}_n^j}{\widehat{\Psi }}_j\left(\widehat{t}\right)\right){\widehat{\varphi }}_i\left(\widehat{t}\right)d\widehat{t}\text{;}v\in S_{\text{H}}$ (17)

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

${\alpha }_{i,j}\text{:}={\displaystyle \underset{\widehat{{\rm I}}}{\int }\frac{d}{d\widehat{t}}{\widehat{\psi }}_{\text{j}}\left(\widehat{t}\right){\widehat{\varphi }}_i\left(\widehat{t}\right)d\widehat{t}}.$ (18)

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

${\displaystyle \sum _{j=0}^k{\alpha }_{i,j}\left(v^{\text{T}}{u}_n^j\right)\approx {\displaystyle \sum _{l=\text{0}}^k{\widehat{w}}_l\{v^{T}f\left(M_n\left({\widehat{t}}_l\right)\text{,}{\displaystyle \sum _{j=0}^k{u}_n^j{\widehat{\psi }}_j\left({\widehat{t}}_l\right)}\right)}}\begin{array}{c}\end{array}{\widehat{\varphi }}_i\left({\widehat{t}}_l\right)\text{;}v\in S_{\text{H}}$ (19)

Where, Ŵ_l represents the weights ${\widehat{t}}_0=\text{-1 ,}{\widehat{t}}_k=\text{1}$, and ${\widehat{t}}_l\in \widehat{{\rm I}};\begin{array}{c}l=\mathrm{1,...,}k-1\end{array}$ are roots of the Legendre polynomial P _k-1(t). Using the abbreviation

$t_{n,l}\text{:}={\text{g}}_n\left({\widehat{t}}_l\right)\text{,}{\beta }_{i\text{,}l}\text{:}={\widehat{w}}_l\widehat{\varphi }\left({\widehat{t}}_l\right)\text{and}{\eta }_{j.l}\text{:}={\widehat{\psi }}_j\left({\widehat{t}}_l\right).$ (20)

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

${\displaystyle \sum _{j=0}^k{\alpha }_{i,j}\left(v^{\text{T}}{u}_n^j\right)=\frac{{\sigma }_{\text{n}}}{\text{2}}{\displaystyle \sum _{l=0}^k{\beta }_{i,l}\left(v^{\text{T}}f\left(t_{n,l}\text{,}{\displaystyle \sum _{j=\text{0}}^k{u}_n^j}{\eta }_{j,i}\right)\right)}}\begin{array}{c}\end{array};v\in S_{\text{H}}\text{,}i=\text{0}\mathrm{,...,}k\text{-1 and}j=\mathrm{0,...,}k,$ (21)

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

${\widehat{\psi }}_i\left({\widehat{t}}_j\right)={\delta }_{\text{i},j}\text{,}{\widehat{\phi }}_{\text{0}}\left({\widehat{t}}_j\right)=-\frac{\text{3}}{\text{4}}{\widehat{t}}_j\text{,}{\phi }_{\text{1}}\left({\widehat{t}}_j\right)=\text{1}\begin{array}{c}\end{array}\text{for all}i,j\in \left\{\text{0,1}\mathrm{,2}\right\},$

Where, each ${\widehat{\psi }}_{\text{j}}\left(\widehat{t}\right)\in {{\rm P}}_{\text{k}}\left(\widehat{{\rm I}}\right),\begin{array}{c}\end{array}{\widehat{\phi }}_i\left(\widehat{t}\right)\in {{\rm P}}_{\text{k-1}}\left(\widehat{{\rm I}}\right)\text{,}$ and ${\widehat{t}}_0=\text{-1 ,}{\widehat{t}}_{\text{1}}=\text{0 ,}$ and ${\widehat{t}}_2=1$. By this choice of test functions _^αi,j , _^βi,j and _^ηj,l are calculated and are given as under

$\{\begin{array}{l}{\alpha }_{\mathrm{0,0}}=-\frac{1}{2},{\alpha }_{\mathrm{0,1}}=1={\alpha }_{\mathrm{1,2}},{\alpha }_{\mathrm{0,2}}=-{\alpha }_{\mathrm{0,0}},{\alpha }_{\mathrm{1,0}}=-{\alpha }_{\mathrm{0,1}},{\alpha }_{\mathrm{1,1}}=\mathrm{0,}\\ {\beta }_{\mathrm{0,0}}=\frac{1}{4},{\beta }_{\mathrm{0,1}}=\mathrm{0,}{\beta }_{\mathrm{0,2}}=-{\beta }_{\mathrm{0,0}},{\beta }_{\mathrm{1,1}}=-3{\beta }_{\mathrm{0,0}},{\beta }_{\mathrm{1,0}}={\beta }_{\mathrm{1,2}}=-\frac{1}{4{\beta }_{\mathrm{1,1}}}\\ {\eta }_{i,j}={\delta }_{i,j},\text{for all}i,j\in \{\text{0,1}\mathrm{,2}\}.\end{array},$

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

${\eta }_n^{m-1}=A,{\xi }_n^{m-1}=0$

Step 4: First solve the following equations for unknowns ${\eta }_n^m,{\xi }_n^m$

$\begin{array}{l}{\eta }_n^m=\frac{1}{\text{2}}{\eta }_n^{m-1}+\frac{1}{\text{2}}{\eta }_n^{m+1}+\frac{{\sigma }_n}{\text{8}}\left({\xi }_n^{m-1}-{\xi }_n^{m+1}\right)\\ {\xi }_n^m=\frac{1}{\text{2}}{\xi }_n^{m-1}+\frac{1}{\text{2}}{\xi }_n^{m+1}+\frac{{\sigma }_n}{\text{8}}\left[-\epsilon \left({\eta }_n^{m-1}-{\eta }_n^{m+1}\right)-\lambda \left\{{\left({\eta }_n^{m-1}\right)}^3-{\left({\eta }_n^{m+1}\right)}^3\right\}\right]\end{array}$

Step 5: Using the computed values for ${\eta }_n^m$ and ${\xi }_n^m$ from Step 4, Solve the following nonlinear equations for ${\eta }_n^{m+1}$ and ${\xi }_n^{m+1}.$

$\begin{array}{l}{\eta }_n^{m+1}={\eta }_n^{m-1}+\frac{{\sigma }_n}{6}\left({\xi }_n^{m-1}+4{\xi }_n^m+{\xi }_n^{m+1}\right)\\ {\xi }_n^{m+1}={\xi }_n^{m-1}-\frac{{\sigma }_n}{6}\left\{\epsilon \left({\eta }_n^{m-1}+4{\eta }_n^m+{\eta }_n^{m+1}\right)+\lambda \left({\left({\eta }_n^{m-1}\right)}^3+4{\left({\eta }_n^m\right)}^3+{\left({\eta }_n^{m+1}\right)}^3\right)\right\}\end{array}$

Step 6: Use the computed values of ${\eta }_n^{m+1}$ and ${\xi }_n^{m+1}$ from Step 5 to get the numerical values of ${\eta }_n^m$ and ${\xi }_n^m$ from Step 4 by back substitution.

Step 7: Compute the solution at time step _^tn by equation (11)

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

${\eta }_{n-1}^{m-1}={\eta }_n^{m+1},{\xi }_{n-1}^{m-1}={\xi }_n^{m+1}.$

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.

Figure 2
Non-dimensional deflection with varying values of amplitude A and with parametric values ( = λ = 1.0.

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

$\frac{(\text{Solution by Present Method - Solution by Existing Method})}{\text{Solution by Existing Method}}\times 100\%,$

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.

Table 1
Comparison of time marching solution of equation (5) with parametric values ( = λ = 1.0.

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

Table 2
Differences between the time marching solution of equation (5) by presented method and different other schemes in literature. The parameters are chosen to be ( = λ = 1.0.

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

Table 3
Relative errors between the time marching solution of equation (5) by presented method with different other schemes in literature. The parameters are chosen to be ( = λ = 1.0.

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.

Figure 3
A comparison of time history of central deflection curve of beam. The time step used for this simulation is 0.05 and the parameters ( = λ = A - 1.0.

Figure 4
A comparison of time history of central deflection curve of beam. The time step used for this simulation is 0.05 and the parameters ( = λ = A - 1.0.

Figure 5
Time-history of central-deflection of beam with varying ( The parametric value of λ is chosen to be 1.0.

Figure 6
Time-history of central-deflection of beam with varying λ The parametric value of ( is chosen to be 1.0.

Table 4
A comparison of accumulated errors and accuracy for 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.

Acknowledgement

The authors are grateful to the anonymous reviewers for their careful reading and constructive suggestions in order to improve the quality of this research article.

References

• Ahmadian, M.T., Mojahedi, M., Moeenfard, H. (2009). Free vibration analysis of a nonlinear beam using homotopy and modified Lindstedt-Poincare Methods. Journal of Solid Mechanics. 29-36.
• Bagheri, S., Nikkar, A., Ghaffarzadeh, H. (2014). Study of nonlinear vibration of Euler-Bernoulli beams by using analytical approximate techniques. Lat. Am. J. Solids Struct. 11: 157-168.
• Barari, A., Ganjavi, B., Ghanbari, J M., Domairry, G. (2010). Assessment of two analytical methods in solving the linear and nonlinear elastic beam deformation problems. Journal of Engineering, Design and Technology, 8(2):127-145.
• Barari, A., Kaliji, H.D., Ghadimi, M., Domairry, G. (2011). Non-linear vibration of Euler-Bernoulli beams. Lat. Am. J. Solids Struct. 8: 139-148.
• Barari, B., Kimiaeifar, A., Domairry, G., Moghimi, M. (2010). Analytical evaluation of beam deformation problem using approximate methods. Songklanakarin Journal of Science and Technology, 32(3):207-326.
• Bayat, M., Pakar, I., Bayat, M., (2011). Analytical study on the vibration frequencies of tapered beams, Lat. Am. J. Solids Struct. 8(2):149-162.
• Bayat, M., Shahidi, M., Barari, A., Domairry, G. (2010). On the approximate analysis of nonlinear behavior of structure under harmonic loading. International Journal of Physical Sciences, 5(7):1074-1080.
• Boukhalfa, A., Hadjoui, A., (2010). Free vibration analysis of an embarked rotating composite shaft using the hp-version of the FEM, Lat. Am. J. Solids Struct. 7:105-141.
• Da Silva, J. C. R. Á., Beck, A. T., da Rosa, E., (2009). Solution of the stochastic beam bending problem by Galerkin method and the Askey-Wiener scheme, Lat. Am. J. Solids Struct. 6(1):51-72.
• Ganji, D. D. (2012). A semi-Analytical technique for non-linear settling particle equation of Motion, Journal of Hydro-environment Research 6(4): 323-327.
• Ganji, S.S., Barari, B., Ganji, D.D., (2011). Approximate analyses of two mass-spring systems and buckling of a column. Computers & Mathematics with Applications, 61(4):1088-1095.
• He, J. H. (2006). Some asymptotic methods for strongly nonlinear equations, International Journal of Modern Physics B 20(10):1141-1199.
• He., J.H. (1999). Variational iteration method - a kind of non-linear analytical technique: Some examples. International Journal of Nonlinear Mechanics, 34:699-708.
• Jafari, S.S., Rashidi, M.M., Johnson, S. (2014). Analytical solution of the nonlinear vibration of euler-bernoulli beams via the homotopy analysis method and differential transform method. Int. J. of Appl. Math and Mech. 10(9): 96-110.
• Javanmard, M., Bayat, M., Ardakani, A. (2013). Nonlinear vibration of Euler-Bernoulli beams resting on linear elastic foundation. Steel Comp. Struct. 15(4):
• Jian-She Peng, J., Yan Liu, Y., Yang, J. (2010). A semi-analytical method for nonlinear vibration of Euler-Bernoulli beams with general boundary conditions. Math. Probl. Eng.
• Khan, Y., Taghipour, R., Fallahian, M., Nikkar, A., (2012). A new approach to modified regularized long wave equation, Neural Computing and Applications doi:10.1007/s00521-012-1077-0
  » https://doi.org/10.1007/s00521-012-1077-0
• Lai, H.Y., Hsu, J.C., Chen, C.K., (2008). An innovative eigenvalue problem solver for free vibration of Euler-Bernoulli beam by using the Adomian Decomposition Method, Computers and Mathematics with Applications, 56: 3204-3220.
• Liu, Y., Gurram, C. S., (2009). The use of He’s Variational iteration method for obtaining the free vibration of an Euler-Beam beam, Mathematical and Computer Modeling 50:1545-1552.
• Mirgolbabaei, H., Barari, A., Ibsen, L.B., Sfahani, M.G. (2010). Numerical solution of boundary layer flow and convection heat transfer over a flat plate. Archives of Civil and Mechanical Engineering, 10:41-51.
• Mirzabeigy, A., Madoliat, R. (2016). Large amplitude free vibration of axially loaded beams resting on variable elastic foundation. Alexandria Engineering Journal 55, 1107-1114.
• Pakar, I., Bayat, M. (2013). An analytical study of nonlinear vibrations of buckled Euler Bernoulli beams, Acta Physica Polonica A 123: 11-15.
• Rafieipour, H., Tabatabaei, S.M., Abbaspour, M. (2014). A novel approximate analytical method for nonlinear vibration analysis of Euler-Bernoulli and Rayleigh Beams on the nonlinear elastic foundation. Arab J. Sci. Eng. 39: 3279-3287.
• Sfahani, M.G., Barari, A., Omidvar, M., Ganji, S.S., Domairry, G. (2011). Dynamic response of inextensible beams by improved Energy Balance Method. In Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamic, 225(1): 66-73.

Publication Dates

• Publication in this collection
  Sept 2017

History

• Received
  29 Aug 2016
• Reviewed
  20 Jan 2017
• Accepted
  22 June 2017