Free vibration and postbuckling of laminated composite Timoshenko beams

1 Introduction

Nowadays, laminated composite structures have several applications in different engineering fields such as aeronautical, automobile, and aerospace industries. Among different laminated composite structures, composite beams are extensively used in civil and mechanical engineering structures because of their high strength-to-weight and high stiffness-to-weight ratios over the conventional isotropic beams. Accordingly, the research topic on buckling [1–10] and vibration [11–20] of these structures has attracted much attention in the literature.

Aydogdu [3] studied the thermal buckling of cross-ply laminated composite beams under different types of boundary conditions based on a three-degree-of-freedom shear deformable beam theory. He solved the problem using the Ritz method with the polynomial series and showed that some cross-ply beams buckle upon cooling rather than heating and some of them do not buckle regardless of whether they are heated or cooled. Gupta et al. [5] applied the Rayleigh-Ritz technique to obtain closed-form expressions for the postbuckling of composite beams subject to different boundary conditions with axially immovable ends. Baghani et al. [6] used the variational iteration method to investigate the postbuckling of unsymmetrically laminated composite beams on a nonlinear elastic foundation. Their study revealed that the beam response can be controlled by a suitable lay-up configuration of the laminate. Emam [7] investigated the significance of the shear deformation on the static postbuckling response of composite beams. It was indicated that classical and first-order theories tend to underestimate the amplitude of buckling load. Vo and Thai [9] employed a refined shear deformation theory to analyze the buckling behavior of composite beams. They developed a two-node C¹ beam element of five degrees of freedom per node to carry out the analysis. Dong and co-workers [13] studied the free vibration of a stepped laminated composite Timoshenko beam considering the effect of shear deformation and rotary inertia. Based on the Timoshenko beam theory (TBT), Çalım [14] investigated the free and forced vibrations of non-uniform composite beams in the Laplace domain. The problem was numerically solved by the complementary function method, and it was concluded that the non-uniformity parameters and angle of fiber orientation have significant effect on the dynamic response of non-uniform composite beams. Ghayesh et al. [16] analytically studied the vibrations and stability of an axially traveling laminated composite beam based on the method of multiple scales. Kim and Wang [17] applied a finite element-based formal asymptotic expansion method to study the free vibration of composite beams. Ghayesh [19] conducted a parametric study on the dynamic response of an axially traveling laminated composite beam with special consideration to natural frequencies, complex mode functions, and critical speeds of the system. Recently, Jafari-Talookolaei and his associates [20] studied the free vibration of laminated composite beams using the TBT and the method of Lagrange multipliers.

Also, it is of great technical importance to study the vibrations of buckled beams. There are several papers in the literature in which the vibration of beams around the buckled configurations has been studied (e.g., [21–23]). In 2008, Nayfeh and Emam [21] published a paper on the postbuckling and vibrations of isotropic Euler-Bernoulli beams. They conducted a closed-form solution for the postbuckling configurations of beams under different boundary conditions in terms of the applied axial load. The vibration around buckled configurations was also analyzed in that work. Later, Emam and Nayfeh [22] extended their previous work on the isotropic Euler-Bernoulli beams to the composite Euler-Bernoulli beams.

Unlike the Euler-Bernoulli beam theory (EBT), the TBT considers the effects of transverse shear deformation and rotary inertia. The difference between the results of these theories becomes more prominent for short beams, which might be encountered in some engineering applications. Motivated by this, the aim of the present work is to extend the study of Emam and Nayfeh [22] based on the EBT to the TBT. By developing a new numerical solution approach, the postbuckling and vibrations around the buckled configurations of laminated composite Timoshenko beams are studied in this paper. It should be noted that the vibration problem around the first three buckled configurations is solved here, while the solution given in [22] is only for vibration around the first buckled configuration. The effect of in-plane inertia is taken into account, and the analysis is carried out for beams under simply supported-simply supported (SS), clamped-clamped (CC), and simply supported-clamped (SC) boundary conditions. Also, the solution is presented for composite beams with both symmetric and unsymmetric lay-up configurations of the laminate. First, the governing equations and boundary conditions are derived via the Hamilton principle. Based on the generalized differential quadrature (GDQ) technique, differential matrix operators are introduced so as to discretize the governing equations and boundary conditions. The pseudo-arc-length continuation method is employed to determine the postbuckling configurations. To study the vibrations around the buckled equilibrium positions, the corresponding eigenvalue problem is numerically solved using the postbuckling configuration modes which are obtained by the pseudo-arc-length method.

2 Timoshenko beam formulation

The stored strain energy in a continuum constructed by a linear elastic material occupying region Ω undergoing infinitesimal deformations is expressed as

where ε_ij and σ_ij are the strain and stress tensors given by

in which u_i denotes the components of the displacement vector u and δ is the Kronecker delta. λ and u can be obtained by [24]

Consider a beam of length L, cross-section A, and thickness h subjected to an axial compressive load N₀. According to the TBT, the displacement components can be written as

where U(t, x), W(t, x), and ψ(t, x) are the axial displacement of the center of sections, lateral deflection of the beam, and the rotation angle of the cross-section with respect to the vertical direction, respectively. The von Karman nonlinear strain-displacement relations are given by

The main components of the symmetric section of the stress tensor can be written in the following forms:

where k_s is a correction factor which shows the non-uniformity of shear strain over the beam cross-section. The strain energy with respect to the initial configuration, Π_s, the kinetic energy of beam, Π_τ, and work done by the axial force, Π_w, can be expressed as

where the normal resultant force N_x, shear force Q_x, and bending moment M_x in a section are defined as

For a composite beam, the stiffnesses are defined as

where n denotes the number of layers, b is the width of the beam, and Z_k and Z_k-1 are the distances from the beam reference surface to the outer and inner surfaces of the kth layer, respectively. Q¯ij are the reduced-transformed stiffnesses of the cross-section of the beam defined as

in which m=cosθ and n=sinθ (θ is the angular orientation of the fibers), and Q is the stiffness matrix. Also, the inertia terms can be given by

where ρ_k is the density of the kth layer.

By applying the Hamilton principle, the governing equilibrium equations and corresponding boundary conditions of a composite beam are obtained as

By introducing the non-dimensional parameters

the governing equations can be rewritten as

where

and boundary conditions are

for CC beams,

for SS beams, and

for SC beams.

3 Differential matrix operators

Consider the column vector F as

in which f(x_i) shows the value of f(x) at each grid point x_i and N symbolizes the number of grid points. Similarly, the values of the rth derivative of f(x) at each point can be defined through the following column vector:

Based on the GDQ method [25], one can arrive at

where the differential matrix operator Dx(r) is computed at grid point values by means of the following formula:

where i, j=1, 2, …, N and δ denotes the Kronecker delta. Also, P_k is expressed as

Using the shifted Chebyshev-Gauss-Lobatto grid points, the mesh generation in the x direction can be given by

4 Pseudo-arc-length method

The pseudo-arc-length continuation [26] is a technique to approximate the solution set of a system of nonlinear equations such as F(V, p)=0, for different values of the system parameter p, in which F: ℝ^N+1→ ℝ^N:(V, p)→F(V, p) and V∈ ℝ^N. Xi=[ViT,  pi]T are defined for the subsequent points on the solution curve {X_i∣i=0, 1, 2, …, k-1 and F(X_i)=0}, where k is the number of solution sets.

First, it is necessary to select an initial point on the solution path. Each point in the algorithm consists of two distinct stages known as the prediction stage and the correction stage. In the prediction stage, the next point on the path, Xi+1p, is predicted using the unit length tangent vector X˙i at X_i and the step size Δs>0,

In the correction stage, in order to improve the predicted point Xi+1p and find the next solution point X_i+1, Newton’s iteration method is applied on the following augmented system:

The main problem in solving this augmented system is to satisfy the constraint F(X_i+1)=0. Furthermore, it is required that X_i+1 be on the hyperplane which passes through Xi+1p and has a normal vector in the direction of the tangent vector X˙i. The augmented system, which is now considered as a square system of linear equations, maps ℝ^N+1 to ℝ^N+1 in which its solution can be found by employing Newton’s method starting from an appropriate initial guess Xi+1p. The tangent vector of the next step can be determined by solving the following system of linear equations:

and after that normalizing ||X˙i+1||=1. The right-hand side of Eq. (29) is a column vector in which all elements are zero except those appeared in the last row. It is worth mentioning that in order to determine the tangent direction X˙i+1 and also the solution of the augmented system of Eqs. (28) in the correction stage, it is necessary to compute Jacobian F_X (X_i+1) which can be calculated numerically using the central finite difference or analytically by symbolic manipulation.

5 Buckling problem

For the buckling problem, in order to obtain the static configuration of the beam, the inertia terms I¯1 and I¯3 are dropped in Eqs. (16). Also, u, w, and ψ are denoted by u_s, w_s, and ψ_s, respectively. The governing equations are discretized together with boundary conditions using the differential operators introduced in Section 3. To this end, grid points in the ξ_i direction are generated using Eq. (26) and column vectors U_s , W_s, and Ψ_s are introduced as

with N entries equal to the number of grid points. In Eq. (30), Usi=us(ξi),  Wsi=ws(ξi), and Ψsi=ψs(ξi). The discretized form of Eqs. (16) for the buckling problem is

where ∘ denotes the Hadamard product (if A=[A_ij]_N×M and B=[B_ij]_N×M, then the Hadamard product of these matrices takes the form A∘B=[AijBij]N×M). In the same way, the boundary conditions can be written as

for CC beams,

for SS beams, and

for SC beams. Note that M=η2(Dξ(1)Us+1α(Dξ(1)Ws)∘(Dξ(1)Ws))+η4Dξ(1)Ψs. Equations (31) can be considered as a nonlinear parameterized boundary value problem of the form

which is solved by the pseudo-arc-length continuation method. It should be noted that with the aim of satisfying the boundary conditions, during the nonlinear solution procedure, the residual of equations of boundaries is inserted into the residual vector of domain that is obtained from the set of nonlinear equations of (35). To accomplish this aim, the elements of the residual vector of domain equivalent to the grid points at the boundaries are replaced with the residual values of boundaries which are obtained from one of the equations of (32), (33), or (34).

6 Vibration around the buckled configurations

To study the free vibration of a buckled beam, a small disturbance around the buckled configurations is considered. Accordingly, the field variables of governing equations can be split into two static and dynamic parts as

where subscript s and d denote the static configuration of a beam in the buckling problem and small dynamic disturbance around the buckled configuration, respectively. By substituting Eq. (36) into (16), one can arrive at

and the boundary conditions become

ud=wd=ψd=0 at ξ=0,1,

for CC beams,

ud=wd=η2(∂ud∂ξ+12α(∂wd∂ξ)2+1αdwsdξ∂wd∂ξ)+η4∂Ψd∂ξ=0=0 at ξ=0,1,

for SS beams, and

for SC beams. In order to investigate the linear vibration around the buckled configurations of beams, the nonlinear terms are ignored and the following relations are substituted into Eqs. (37):

The final form of governing equations for the vibration is obtained as

η3(d2w˜ddξ2+αdψ˜ddξ)+η1α(dusdξd2w˜ddξ2+d2wsdξ2d2u˜ddξ)+η1α(d2usdξ2dw˜ddξ+dwsdξd2u˜ddξ2)+3η12α2((dwsdξ)2d2w˜ddξ2+2dwsdξd2wsdξ2dw˜ddξ)+η2α(dψsdξd2w˜ddξ2+d2wsdξ2dψ˜ddξ)+η2α(d2ψsdξ2dw˜ddξ+dwsdξd2ψ˜ddξ2)+N¯0d2w˜ddξ2=-ω2w˜d(42b)

η4d2ψ˜ddξ2-η3α(d2w˜ddξ+αψ˜d)+η2d2u˜ddξ2+η2α(dwsdξd2w˜ddξ2+d2wsdξ2dw˜ddξ)=-ω2(ψ˜d+I¯1u˜d)(42c)

where ω is the natural frequency and [u˜d w˜d ψ˜d]T stands for the corresponding mode shape. Using the differential matrix operators, Eqs. (42) are discretized as follows:

and the boundary conditions can be written as

for CC beams,

for SS beams, and

for SC beams. Note that Md=η2Dξ(1)Ud+η4Dξ(1)Ψd. Equations (43) can be shown in the following form:

where