Computers and Structures 83 (2005) 2632–2649

www.elsevier.com/locate/compstruc

Buckling, flutter and vibration analyses of beams

by integral equation formulations
Z. Elfelsoufi, L. Azrar *

Equipe de Modélisation Mathématique de Problèmes Mécaniques, UFR SPI, Faculté des Sciences et Techniques de Tanger,
Université Abdelmalek Essaädi, BP. 416, Tanger, Morocco

Received 20 January 2004; accepted 14 March 2005

Abstract

This paper presents a model for the investigation of buckling, flutter and vibration analyses of beams using the inte-
gral equation formulation. A mathematical formulation based on Euler–Bernoulli beam theory is presented for beams
with variable sections on elastic foundations and subjected to lateral excitation, conservative and non-conservative
loads. Using the boundary element method and radial basis functions, the equation of motion is reduced to an algebraic
system related to internal and boundary unknowns. Eigenvalue problems related to buckling and vibrations are formu-
lated and numerically solved. Buckling loads, natural frequencies and associated eigenmodes are computed. The cor-
responding slope, bending and shear forces can be directly obtained. The load-frequency dependence is investigated
for various elastic foundations and the divergence critical loads are evidenced. Under non-conservative loads, a
dynamic stability analysis is illustrated numerically based on the coalescence of eigenfrequencies. The flutter load
and instability regions with respect to the elastic concentrated and distributed foundations are identified. Using the
eigenmodes, numerically computed, non-linear vibrations of beams are investigated based on one mode analysis.
The presented model is quite general and the obtained numerical results are in agreement with available data.

2005 Elsevier Ltd. All rights reserved.

Keywords: Beam; Buckling; Flutter; Vibration; Integral equation; BEM; Radial basis functions

1. Introduction on a better design of light weight structures which can be

safely used in the pre-buckling and post buckling ranges.
Slender structural components such as beams and Amongst the numerical methods available for thin struc-
columns constitute basic parts of many complex engi- ture problems, the finite element method is undoubtedly
neering structures. Buckling, flutter and vibrations are the most versatile. The only problem with this method is
the main forms of instability of these structures. The that its formulation is quite laborious and it takes a large
accurate prediction of static and dynamic critical loads amount of computer storage. A powerful alternative
and responses are very significant in practice and it relies method based on integral equations is the Boundary Ele-
ment Method. One of the main reason of the rapid
*
Corresponding author. Tel.: +212 39 39 39 54/55; fax: +212 development of the BEM is a possibility of reduced
39 39 39 53. dimensionality of the problem, which leads to a reduced
E-mail addresses: azrar@fstt.ac.ma, azrar@lpmm.univmetz. set of equations and a smaller amount of data required
fr, azrar@hotmail.com (L. Azrar). for the computation. Using the fundamental solution

0045-7949/$ - see front matter
2005 Elsevier Ltd. All rights reserved.
doi:10.1016/j.compstruc.2005.04.001
 Z. Elfelsoufi, L. Azrar / Computers and Structures 83 (2005) 2632–2649 2633

corresponding to the exact solution of a part of the method must be employed. The instability analysis of
problem, the inappropriate terms are moved to the beams under non-conservative force characterized by
right-hand side of the governing equation and consid- the flutter which occurs when two of the natural fre-
ered as a fictitious source density. For buckling and quencies coincide became complex conjugate can be
vibrations of beams under elastic foundations, domain done only dynamically. The presence of non-conserva-
integrals are necessary in the formulation. Thus, the tive loads makes the system of equations mathematically
main advantage of dimensionality reduction is elimi- non-self-adjoint and the corresponding eigenvalue prob-
nated. But, the use of dual reciprocity method (DRM), lem is dictated by a non-symmetric matrix. This problem
introduced by Nardini and Brebbia [1], permits to com- has been studied by many authors and various numeri-
bine the dimensionality reduction advantage with a sim- cal procedures have been set up for its solution. A com-
ple fundamental solution and to formulate the problem prehensive discussion of this subject based on analytical
on boundary unknowns only. A comprehensive litera- procedures can be found in the books [15–17]. More re-
ture review of the DRM and multiple reciprocity cently, the flutter and divergence instability of beams
method (MRM) as applied to elastodynamics can be and plates subjected to non-conservative loads are ana-
found in the review paper of Beskos [2]. Details and lyzed by Zuo and Schreyer [18] by solving the resulting
applications to various engineering problems are clearly characteristic equations. Similarly, the flutter and inter-
presented in the book of Partridge et al. [3]. The eigen- nal damping effects on the dynamic stability of rods with
value analyses of Helmhotz equation using the DRM intermediate spring support and with relocatable lamped
and MRM have been discussed by Kamiya et al. [4]. mass under follower loads have been largely investigated
Combining the MRM and singular value decomposition by Lee [19–21]. The influence of the tangency coefficient
method, the rod vibration problem has been analyzed by of follower load and the elastically restrained boundary
Chang et al. [5]. The solution of plate bending problems conditions on the elastic instability of beams has been
by MRM has been formulated by Sladek and Sladek [6]. discussed by Lee and Hsu [22]. Enhancing flutter and
Using DRM and differential quadrature method, the buckling capacity of beams by using piezoelectric layer
longitudinal vibrations of plates and membranes are is presented by Wang and Quek [23]. Based on the finite
investigated by Tanaka and Chen [7]. For bending prob- element method, the stability and instability of cantile-
lems of inhomogeneous Euler–Bernoulli beams, an ver elastic beams subjected to a follower force have been
investigation is carried out by Rong et al. [8]. Based investigated by Gasparini et al. [24] and by Ryu and
on Timoshenko beams theory and a quadrature method, Sugiyama [25]. Using the static approach, the divergence
the dynamic behaviors of beams have been analyzed by instability of thin walled beams in pre-buckling and
Schanz and Antes [9]. An extension to beams with arbi- post-buckling ranges has been recently analyzed by
trary cross-section has been developed by Sapountzakis Mohri et al. [26]. Based on the Ramm finite elements
[10] and to the non-linear dynamic analysis of beams and a perturbation method, the load-frequency depen-
with variable cross-section has been done by Katsikad- dence has been investigated for arches and shells with
elis and Tsiatas [11]. Vibrations of beams with variable large rotations by Boutyour et al. [27]. The critical loads
sections using BEM and radial basis functions are ana- and the stability and instability regions are evidenced for
lyzed in [12]. However, no general modeling based on largely deformed shells with the smaller eigenfrequen-
integral equation formulations concerning buckling, cies. The divergence and flutter instability are generally
vibration and flutter analyses of beams has been estab- analyzed by analytical method or by finite element meth-
lished in previous works. This paper intends to provide ods. To the best knowledge of the authors, there is no
a compact formulation for these behaviors and to available compact formulation and results based on
numerically investigate the static and dynamic instabili- the integral equation formulation for buckling, flutter
ties of beams. and vibration analyses of thin structures.
The force acting on a beam can be divided into con- In this paper, a mathematical modeling based on the
servative and non-conservative forces. Generally, the integral equations for buckling, flutter and vibration
instability of a beam under a conservative force is char- analyses of beams is presented. The Euler–Bernoulli
acterized by the divergence which occurs when one of beam theory is used and the governing equation is for-
the natural frequencies falls to zero. This critical solu- mulated for beams on elastic foundation and subjected
tion corresponds to buckling load and can be directly to conservative and non-conservative loads. The radial
investigated with static analysis. The buckling analysis basis functions are used and the required matrices are
of beams has been studied by many authors and is trea- explicitly presented for various boundary conditions
ted in almost any textbook on mechanics of solids [13]. and loads. The eigenvalue problems corresponding to
For a non-conservative system, it has been shown by buckling, vibration and flutter are explicitly formulated.
Ziegler [14] that the usual Euler method and minimum The displacement, slope, bending and shear forces can
potential energy methods (static methods) are inade- be directly obtained. The buckling and vibration modes
quate to predict their instability and that a dynamic and the load-frequency dependence are presented for
2634 Z. Elfelsoufi, L. Azrar / Computers and Structures 83 (2005) 2632–2649

various boundary conditions and elastic foundations. where K1 and K2 are functions of the axial coordinate z.
For axial tangential follower forces, the flutter load is Using non-dimensional parameters, Eq. (2) can be read
evidenced for uniform and concentrated elastic founda- as
tions based on the coalescence criterion. The flutter load

o2 o2 W o2 W
variations with respect to the position and amplitude of 2
K 1 ðxÞ 2  x2 K 2 ðxÞW þ k 2 þ j ðxÞW ¼ 0
the foundation and the flutter zone are investigated. The ox ox ox
non-linear vibration analyses of beams based on one ð4Þ
mode numerically computed are elaborated. Some where
benchmark tests are investigated ratifying the effective-
ness of the presented methodological approach. Sð0ÞL4 x2 L2 L3
x2 ¼ q ; k ¼ k ; j ¼ j ;
EIð0Þ EIð0Þ EIð0Þ
rffiffiffi
V I z
2. Basic beam equations W ¼ ; R¼ and x ¼
R S L

Let us consider a slender beam of length L with a in which R is the radius of gyration of the beam and
variable cross-section. The Euler–Bernoulli beam for- 0 6 x 6 1.
mulation based on the assumption that both shear Eq. (4) may be solved by the finite element method or
deformation and rotational inertia of the cross-section analytically for standard boundary conditions. The aim
are negligible is used. The axial displacement will be ne- of this paper is the development of an integral equation
glected and the equation of motion is formulated using formulation for numerical solution of (4) and the inves-
the transverse displacement only. Based on the integral tigation of the static and dynamic instability analyses of
equation formulation, buckling, flutter and transverse beams based on the resulting formulation.
vibrations of beams will be formulated. The governing
partial differential equation of motion of beams on elas-
tic foundations and subjected to axial compression and 3. Integral equation formulation
lateral excitation (Fig. 1) is formulated by

 The fundamental solution of (4) is hard to be explic-
o2 o2 V ðz; tÞ o2 V ðz; tÞ o2 V ðz; tÞ itly determined due to variable coefficients K1(x) and
2
EIðzÞ 2
þ qðzÞSðzÞ 2
þk
oz oz ot oz2 K2(x) even if only buckling (x = 0) or free vibration
þ jðzÞV ðz; tÞ ¼ pðz; tÞ ð1Þ (k = 0) problem is considered. For simplified cases, Bes-
sel functions may be used but will lead to some numer-
where V is the transverse displacement, E, I, S and q are ical difficulties at standard integral equation formulation
YoungÕs modulus, inertia, the area and the mass density of solving the resulting boundary value problem [3–5].
respectively. j(z) is the elastic foundation, k is the axial As the domain integrals are inevitable due to the excita-
compression, p(z, t) is the lateral excitation and z is the tion and load, a simple fundamental solution will be
axial coordinate. Assuming harmonic motion, the free used and the resulting domain integrals will be treated
vibration problem of axially loaded beam is given by by the dual reciprocity method. Let us denote W the

 fundamental solution of the following problem:
o2 o2 V o2 V
2
EIðzÞ 2
 qSðzÞx2 V þ k 2 þ jðzÞV ¼ 0 ð2Þ

oz oz oz o2 o2 W  ðx; sÞ
K 1 ðxÞ ¼ dðx; sÞ ð5Þ
For homogeneous beams with a variable section, the ox2 ox2
parameters E, I, S and q can be assumed in the following
form: where d is the Dirac function and s is the source point.
This fundamental solution will be used and the differen-
EIðzÞ ¼ Eð0ÞIð0ÞK 1 ðzÞ tial Eq. (4) will be transformed into an integral equation.
ð3Þ
qSðzÞ ¼ qð0ÞSð0ÞK 2 ðzÞ Following the boundary element method procedure

z r p(z,t)

λ
←
z : Field s : Source κ ( z)
point point
(Moving) (Marked)

Fig. 1. C–S beam subjected to axial force k, elastic foundation j and lateral excitation p(z,t).
 Z. Elfelsoufi, L. Azrar / Computers and Structures 83 (2005) 2632–2649 2635

[1–12], the resulting integral equation will be reduced to More details about W , fj, gj and hj, used in this anal-
algebraic equation. ysis, are given in the Appendix A. Making use of these
As well known in the bending problem of beams, the transformations, it is now possible to evaluate the inte-
following variables have physical meanings and may be gral formulation (7) using boundary values only. Based
also known at boundaries: on the decomposition (9) and Eq. (10), the three domain
integral in the r.h.s. of (7) are transformed into bound-
oW o2 W oM ary values as follows:
hðxÞ ¼ ; MðxÞ ¼ K 1 ðxÞ and QðxÞ ¼
ox ox2 ox 8Z 1
> Xnþ2
ð6Þ >
> K 2 ðxÞW ðxÞW  ðs; xÞdx ¼ aj Bj ðsÞ
>
>
>
> 0 j¼1
>
> 
where h(x) is the slope, M(x) is the bending moment and >
>
>
> d3 h
Q(x) is the shear force related to the derivatives of the < Bj ðsÞ ¼ hj ðsÞ þ W  ðs; xÞ 3j ðxÞ
dx
deflection W. Multiplying Eq. (4) by W  and integrating > 2
>
>
> oW 
d hj o2 W  dhj
from 0 to 1, one obtains: >  ðs; xÞ 2 ðxÞ þ ðs; xÞ ðxÞ
>
> ox dx ox2 dx
Z 1 2
 >
> 
o2 W >
> 1
o >
> o3 W 
2
K 1 ðxÞ ðxÞ W  ðs; xÞdx :  3 ðs; xÞhj ðxÞ
0 ox ox2 ox 0
Z 1 ð11aÞ
¼ x2 K 2 ðxÞW ðxÞW  ðs; xÞdx 8Z
0
>
>
1 2
oW Xnþ2
Z >
> ðxÞW 
ðs; xÞdx ¼ aj C j ðsÞ
1
o2 W >
> ox2
 k ðxÞW  ðs; xÞdx >
> 0 j¼1
ox2 >
> "
0 >
> d2 gj d5 gj
>
> 
Z 1
>
< C j ðsÞ ¼ ðsÞ þ W ðs; xÞ ðxÞ
dx2 dx5
 j ðxÞW ðxÞW  ðs; xÞdx ð7Þ
>
> oW  d4 gj o2 W  d3 g j
0 >
>
>
>  ðs; xÞ 4 ðxÞ þ ðs; xÞ 3 ðxÞ
Integrating by parts four times, the first term of (7) >
> ox dx ox 2 dx
>
> #1
>
> o3 W  d2 g j
becomes >
>
>
:  3 ðs; xÞ 2 ðxÞ
8Z 1 2
 ox dx
>
> o o2 W 0
>
>
> 2
K 1 ðxÞ ðxÞ W  ðs; xÞdx ¼ W ðsÞ þ AðsÞ ð11bÞ
>
> 0 ox ox2
>
>
>
> 
< oW  For a uniform elastic foundation j
AðsÞ ¼ W  ðs;xÞQðxÞ þ ðs; xÞMðxÞ 8Z 1
>
> ox X
nþ2
>
> >
>
>
> >
> j ðxÞW ðxÞW  ðs; xÞdx ¼ j aj Dj ðsÞ
>
>
 1 >
>
>
> o2 W  o o2 W  >
>
0 j¼1
>
: þK 1 ðs;xÞhðxÞ  K1 ðs; xÞW ðxÞ >
> "
ox2 ox ox2 >
> d3 gj
0 >
> 
< Dj ðsÞ ¼ gj ðsÞ þ W ðs; xÞ dx3 ðxÞ
>
ð8Þ
>
> oW  d2 gj o2 W  dgj
In the right-hand side of Eq. (7), three domain inte- >
>
>
>  ðs; xÞ ðxÞ þ ðs; xÞ ðxÞ
grals have to be evaluated. Making use of radial basis >
> ox dx 2 ox 2 dx
>
> #1
functions, thus avoiding the additional task of domain >
> o3 W 
>
>
>
:  3 ðs; xÞgj ðxÞ
integration, these domain integrals are transformed into ox
0
boundary values. Let us assume that for these integrals,
the function W(x) is assumed to be: ð11cÞ

X
nþ2 For a concentrated elastic foundation j at point
W ðxÞ ¼ aj fj ðxÞ ð9Þ 8Z
> 1 X
nþ2
j¼1 >
> j
ðxÞW ðxÞW 
ðs; xÞdx ¼ j
aj Dj ðsÞ
>
< 0
j¼1
where fj are radial basis functions, ÔnÕ is the number of ð11dÞ
interior points, aj are unknown coefficients [3,4]. Given >
> Dj ðsÞ ¼ 0 j 6¼ L;
>
>
fj defines two other functions gj and hj which satisfy :
DL ðsÞ ¼ W ðxL ÞW  ðs; xL Þ; aL ¼ 1
the following equations:
Based on Eqs. (8) and (11), the integral formulation (7)
d4 gj d4 hj is reduced to the following algebraic equation at interior
ðxÞ ¼ fj ðxÞ and ðxÞ ¼ K 2 ðxÞfj ðxÞ ð10Þ
dx4 dx4 points:
2636 Z. Elfelsoufi, L. Azrar / Computers and Structures 83 (2005) 2632–2649
8
X
nþ2 X
nþ2
>
>
nPþ2 nPþ2 nPþ2
W ðsÞ þ AðsÞ ¼ x2 aj Bj ðsÞ þ k aj C j ðsÞ >
> W i þ Ai ¼ x2 aj Bij þ k aj C ij þ j aj Dij
>
>
j¼1 j¼1 >
>
j¼1 j¼1 j¼1
>
>
X
nþ2 >
> P
nþ2 P
nþ2 P
nþ2
>
> h þA b i ¼ x2 b ij þ k
aj B b ij þ j
aj C b ij
aj D
þ j aj Dj ðsÞ ð12Þ >
< i j¼1 j¼1 j¼1
j¼1
>
> b nPþ2
b nPþ2
b nPþ2
b
>
> b i ¼ x2
M i þ A b ij þ k
aj B b ij þ j
aj C b ij
aj D
where A, B, C and D will be explicitly given later. >
>
>
> j¼1 j¼1 j¼1
The boundary conditions of the beam may be classi- >
>
>
> b nP b nP b nP b
fied as simply-supported, clamped, free or more general >
> b
b i ¼ x2
þ2
b
b ij þ k
þ2
b
b ij þ j
þ2
b
b ij
>
: Qi þ A aj B aj C aj D
as elastically supported edges. In order to present a gen- j¼1 j¼1 j¼1
eral formulation for various boundary conditions, more ð15a–dÞ
equations than (12) related to h, M and Q are needed.
They are obtained by derivatives of Eq. (12) according in which i = 1 and i = n + 2 correspond to beam ends
to the variable ÔsÕ. For a compact equation representa- and i = 2 to i = n + 1 correspond to interior points
tion, the following notations are introduced: which may correspond to uniform or non-uniform dis-
cretisation. Let us recall that we have (n + 4) unknowns,
oE b
b
EðsÞ ¼ ðsÞ; b ðsÞ ¼ K 1 ðsÞ o E ðsÞ and
b
E (n) interior and four unknowns related to the assumed
os os boundary conditions. Eq. (15a) leads to (n + 2) equa-
b b
b tions and two extra equations are then needed. Eqs.
b ðsÞ ¼ o E ðsÞ
b
E ð13Þ (15b) or (15c) may be used to complete the system for
os
S–S, S–C, C–C and other boundary conditions. A com-
where E may be A, B, C or E = D. bination of Eqs. (15a)–(15d) can also be used to solve
Finally, one obtains for a uniform elastic foundation the system for more general boundary conditions.
the following algebraic system: The present formulation is quite general and is given
8 nP þ2 nPþ2
for beams with variable cross-sections. As the aim of this
>
> W ðsÞ þ AðsÞ ¼ x2 aj Bj ðsÞ þ k aj C j ðsÞ paper is the development of a simple and general formu-
>
>
>
> j¼1 j¼1 lation leading to a unique approach for buckling, flutter,
>
>
>
> nPþ2 transverse vibrations and interaction with elastic foun-
>
>
>
> þj aj Dj ðsÞ dation, beams with constant sections are considered.
>
>
>
>
j¼1
For the dynamic analyses of beams with variable
>
> nP nP
>
> b
þ2 þ2
b j ðsÞ
> hðsÞ þ AðsÞ
> ¼ x2 b j ðsÞ þ k
aj B aj C cross-sections, more theoretical developments with re-
>
>
>
> j¼1 j¼1 spect to the associated fundamental solution and radials
>
>
>
> nP þ2 basis functions are needed. For the sake of clearness, let
>
> þj b j ðsÞ
aj D
>
< us rewrite (15) in a matrix form and give more details
j¼1
about introduced matrices for specified boundary
>
> b nP þ2
b nPþ2
b conditions.
>
> b ðsÞ ¼ x2
MðsÞ þ A b j ðsÞ þ k
aj B b j ðsÞ
aj C
>
> Assuming that the beam is clamped–simply-sup-
>
> j¼1 j¼1
>
>
>
> nPþ2
b ported as presented in Fig. 2 (W1 = 0, h1 = 0, Wn+2 = 0
>
> þj b j ðsÞ
aj D
>
> and Mn+2 = 0) and introducing the following notations:
>
> j¼1
>
>
>
> b b P b oW  j1  sj
b ðsÞ ¼ x2 P aj B
nþ2 nþ2
>
> b b b
>
> QðsÞ þ A b j ðsÞ þ k aj Cb j ðsÞ A1 ðsÞ ¼ 
ox
ðs; 0ÞM 1 ; A2 ðsÞ ¼
2
hnþ2 ;
>
>
>
>
j¼1 j¼1
>
>
> A3 ðsÞ ¼ W  ðs; 0ÞQ1 ; A4 ðsÞ ¼ W  ðs; 1ÞQnþ2
> nP þ2 b
b
>
: þj aj Db j ðsÞ fW g ¼ fW 2 ; W 3 ; . . . ; W n ; W nþ1 g and
j¼1

ð14a–dÞ fT g ¼ fM 1 ; hnþ2 ; Q1 ; Qnþ2 g

These equations give analytical solution representations The vector {W} represents the deflection unknowns
with respect to the interior variable s. For a numerical at interior point and {T} represents the boundary un-
solution, a discretization of (14) and the consideration knowns. For the considered boundary conditions and
of boundary conditions are needed. using previous notations, the algebraic system (15) can
be written in the following matrix form:
  
     
I A fW g B C D
4. Matrix formulations ¼ x2 þ k þ j F fW g
O A0 fT g B0 C0 D0
After discretization of Eq. (14), one can write ð16Þ
 Z. Elfelsoufi, L. Azrar / Computers and Structures 83 (2005) 2632–2649 2637
(
A0 fT g ¼ ðx2 B0 þ k C 0 þ j D0  OÞF fW g
ð17a; bÞ
AfT g ¼ ðx2 B þ k C þ j D  I ÞF fW g

Let us recall that the vectors {W} and {T} and the
1/(n+1)
previous matrices depend on the boundary conditions
1 1 3 n+1 n+2 considered. Matrices I and O may be also changed
according to boundary conditions and the same nota-
Fig. 2. Uniform discretization of a simply-supported and tions are kept for a general representation. In Appendix
clamped beam. W1 = Wn+2 = 0, h1 = 0 and Mn+2 = 0.
B, details about vectors and matrices for simply-sup-
ported, clamped–clamped and clamped–free are given.
The solutions of the considered problems are obtained
where I (n · n), and O (4 · n) are the identity and zero by numerically solving Eq. (16). This system is reduced
matrix, respectively. to an eigenvalue problem in deflection vector only as
follows:
A : matrix ðn  4Þ; Aik ¼ Ak ðsi Þ; for k ¼ 1; 2; 3 or 4
8
si ¼ ði  1Þ=ðn þ 1Þ and i ¼ 2 to n þ 1 < ½X  fW g ¼ ½Y  fW g
>
2
X ¼ ðI  AA1  1
0 OÞ  x ðAA0 B0  BÞF
>
:
A0 : matrix ð4  4Þ; A01k ¼ Ak ð0Þ; A02k ¼ Ak ð1Þ; 
Y ¼ ðk ðC  AA1  1
0 C 0 Þ þ j ðD  AA0 D0 ÞÞF
b b
b k ð0Þ; A04k ¼ A
b k ð1Þ; for k ¼ 1 to 4 ð18a–cÞ
A03k ¼ A
The unknowns at boundaries can be easily computed
B : matrix ðn  ðn þ 2ÞÞ; Bij ¼ Bj ðsi Þ;
by the following algebraic equation:
i ¼ 2 to n þ 1 and j ¼ 1 to n þ 2

fT g ¼ A1 2 
0 ½ðx B0 þ k C 0 þ j D0 ÞF  O fW g ð19Þ
B0 : matrix ð4  ðn þ 2ÞÞ; B01j ¼ Bj ð0Þ; B02j ¼ Bj ð1Þ; This permits to break down the matrix problem (16) into
b b
b j ð0Þ; B04j ¼ B
b j ð1Þ; j ¼ 1 to n þ 2 an eigenvalue problem and an algebraic one. Matrices
B03j ¼ B
[X] and [Y] can be easily formulated for each considered
C : matrix ðn  ðn þ 2ÞÞ; C ij ¼ C j ðsi Þ; problem.
For specified boundary conditions, load and founda-
i ¼ 2 to n þ 1 and j ¼ 1 to n þ 2
tions, matrices in (18) and (19) have first to be computed
following the developments presented in Appendix B.
C 0 : matrix ð4  ðn þ 2ÞÞ; C 01j ¼ Cð0Þ; C 02j ¼ C j ð1Þ; For numerical solutions, a computing program in MAT-
b b
b j ð0Þ; C 04j ¼ C
b j ð1Þ; j ¼ 1 to n þ 2 LAB has been developed. The MATLAB environment is
C 03j ¼ C
exploited for a standard use of the presented formula-
D : matrix ðn  ðn þ 2ÞÞ; Dij ¼ Dj ðsi Þ; tion. Formulations (18) and (19) are quite general and
allow one to investigate buckling, vibrations and combi-
i ¼ 2 to n þ 1 and j ¼ 1 to n þ 2
nation of them leading to load-frequency dependence
for beams on various types of elastic foundations.
D0 : matrix ð4  ðn þ 2ÞÞ; D01j ¼ Dj ð0Þ; D02j ¼ Dj ð1Þ;
b
b j ð0Þ; D04j ¼ Db
b j ð1Þ; j ¼ 1 to n þ 2 4.1. Buckling problem
D03j ¼ C

F 1 : matrix ððn þ 2Þ  ðn þ 2ÞÞ of radial function matrix The buckling problem may be formulated by omit-

 ting the frequency parameter in (18) and (19). The criti-
ði  1Þ cal buckling loads and associated eigenmodes can be
F 1ij ¼ fj ; i ¼ 1 to n þ 2 and j ¼ 1 to n þ 2
ðn þ 1Þ determined for various types of boundary conditions
and elastic foundations j by solving the following
F : matrix ððn þ 2Þ  ðnÞÞ; F ij ¼ F 11
iðjþ1Þ ; eigenvalue problem:
i ¼ 1 to n þ 2 and j ¼ 1 to n 8
>
> ½X fW g ¼ k1 fW g
because for clamped simply-supported boundary <  1
conditions W 1 ¼ W nþ2 ¼ 0 X ¼ ðI  AA1  1
0 OÞ  j ðD  AA0 D0 ÞF
ð20a; bÞ
>
>
: 1
ðC  AA0 C 0 ÞF
Eq. (16) presents an algebraic system on the deflec-
tion at unknown interior points {W} and unknowns at This allows the investigation of numerical critical
boundaries represented by {T}. This system is rewritten buckling loads and associated eigenmodes at interior
as: points. The corresponding slope, moment and shear
2638 Z. Elfelsoufi, L. Azrar / Computers and Structures 83 (2005) 2632–2649

force can be numerically computed using (15). The un- can be directly computed. The control of the linear and
knowns at boundaries are computed by the following non-linear flutter may be performed based on the
algebraic equation: numerically obtained modes and will be an extension

of this work.
fT g ¼ A1 
0 ½ðk C 0 þ j D0 ÞF  O fW g ð20cÞ
4.4. Nonlinear vibration problem
4.2. Linear vibration problem
One of the main objectives of the present work is to
The numerical solution of Eq. (18) permits, on one establish a multi-modal formulation based on boundary
hand, to study the linear vibration behaviors of beams element method for non-linear vibration and post-buck-
by omitting the load parameter. On the other hand, ling analyses of beams. The numerical solution of the
the linear vibration analysis of beams under an axial linear vibration problem (21) permits one to get the nat-
compression also can be investigated by numerically ural frequencies and associated eigenmodes. Using the
solving the following eigenvalue problem for each fixed obtained eigenmodes, a multimodal formulation can
load parameter: be developed for beams with various boundary condi-
8 tions and elastic foundations. Based on harmonic bal-
>
> ½X fW g ¼ x12 fW g ance method, a semi analytical method has been
< 
X ¼ I  AA1  1 presented for nonlinear free and forced vibrations of
0 O þ k ðC  AA0 C 0 ÞF ð21Þ
>
> 1 1 beams [28,29]. In that work, analytical beam modes,
:  1
þj ðD  AA0 D0 ÞF ðA0 B0  BÞF available for classical boundary conditions, are used.
In the present work, the numerically obtained modes
The unknowns at boundaries can be computed by the will be used and any boundary condition and elastic
algebraic Eq. (18). The load-frequency dependence can foundation may be inserted. The nonlinear harmonic re-
be investigated and the divergence stability may be ana- sponse may be computed in the following form:
lyzed for various types of elastic foundations.
X
n
W ðx; tÞ ¼ cosðxtÞ ai wi ðxÞ ð23Þ
4.3. Flutter problem i¼1

Let recall that the matrix X is load dependent. With where wi(x) are the computed eigenmodes and {A}t =
adjusted matrices, in accordance with the considered {a1, a2, . . . , an} is the corresponding amplitude vector.
boundary conditions, Eq. (21) can be used for conserva- Following the formulation presented in [28,29], the non-
tive and non-conservative loads. When the applied load linear forced vibration of beams can be analyzed using
is a non-conservative follower force, the frequencies can the following a dimensional and multidimensional
be either real or complex. Therefore, at divergence insta- formulation:
bility, the lowest frequency vanishes, as for the conserva-
tives system or two frequencies can approach each other, 3
ð½K   x2 ½M  ÞfAg þ ½B ðAÞ fAg ¼ fF  g ð24Þ
coalesce and then become complex conjugate. This cor- 2
responds to flutter instability and the load at the two fre- where [M ], [K ] and [B (A)] are the mass matrix, the lin-
quencies coincide is defined as the flutter load. This ear and nonlinear rigidity matrices respectively. {F } is
study is intended to extend the previous analysis based the column vector of generalized transverse excitations.
on boundary integral formulation to the stability of a The multimodal analysis will be investigated in the next
cantilevered beam subjected to a tangential follower work. In this work, we limit ourselves to one mode anal-
force at the free end (BeckÕs problem) under elastic foun- ysis and the response is assumed to be
dation. The assumed boundary conditions for a C–F
beam on an elastic foundation at the free end is: V ðz; tÞ ¼ Ra1 w1 ðxÞ cosðxtÞ
in which R is the radius of gyration and a1 the amplitude
W 1 ¼ 0; h1 ¼ 0; M nþ2 ¼ 0 and Qnþ2 ¼ j  W nþ2 corresponding to the first mode. The 1-D nonlinear fre-
ð22Þ quency response function is then given by [28]

 2
The dynamic stability analysis can be performed by x 3 b1111 2 1 f1
¼1þ  a1   ð25Þ
Eq. (21) with adjusted matrices taking into account the xL
2 k 11 k 11 a1
assumed boundary conditions as presented in Appendix
B. The load-frequency dependences and the flutter load where
corresponding to coalescence of two natural frequencies Z 1
can be investigated. The mode, moment and shear force xL2 ¼ k 11 =m11 ; m11 ¼ ðw1 ðxÞÞ2 dx
corresponding to critical frequency or to the flutter load 0
 Z. Elfelsoufi, L. Azrar / Computers and Structures 83 (2005) 2632–2649 2639

Z
2 Z 1 beam tests are analyzed and only some benchmark ones

1
d2 w1 ðxÞ
k 11 ¼ 2
dx þ j ðw1 ðxÞÞ2 dx are presented in order to demonstrate the effectiveness of
0 dx 0
the developed approach.
in the case of a uniform elastic foundation
Z 1
2 2
d w1 ðxÞ 5.1. Buckling
k 11 ¼ dx þ j ðw1 ðx0 ÞÞ2
0 dx2
The obtained results depend on the number of inter-
in the case of a concentrated elastic foundation
nal points considered. The convergence to the analytical
Z
2 !2
 1 1 dw1 ðxÞ solution is tested and results are presented. For the sake
b1111 ¼ dx of brevity, only few tested case are presented. The con-
2 0 ðdxÞ
vergence of the solution with respect to interior points
For a concentrated harmonic force of amplitude Fcc at ÔnÕ is presented in Tables 1 and 2. It can be seen clearly
x0, f1 is given by from these tables that 60 interior points is largely en-
L3 ough for accurate results. Results obtained with the
f1 ¼ F c w1 ðx0 Þ ¼ F c0r w1 ðx0 Þ ð26Þ present model for critical buckling loads are favorably
EI R
compared to analytical ones in case of constant sections.
The coefficients m11 ; k 11 ; b1111 and w1(x0) are evaluated Results are presented only for clamped–simply sup-
using the first mode shape numerically obtained by the ported and clamped–free beams. For buckling of beams
present analysis. The nonlinear frequency–amplitude on a concentrated elastic foundation, the critical load
dependence for beams under various boundary condi- changes according to the amplitude and position of the
tions and elastic foundations for free and forced vibra- foundation. Numerical results obtained by 60 interior
tion can be easily analyzed based on Eq. (25). points for S–S beams are presented in Table 3 and are
favorably compared to the available analytical ones
[13]. The effect of the concentrated foundation at
5. Numerical results Xc = 0.75 on the fifth first modes are presented in
Fig. 3 for j = 150 The same study is done for a C–F
The mathematical formulation is presented for beams beam (Xc = 1) leading to the same behavior as presented
with variable section under various elastic foundation in Fig. 4. The concentrated elastic foundation effect on
and boundary conditions and subjected to various types the first and higher buckling modes are clearly shown.
of loads. In this paper, the analysis is limited to isotropic
beams with constant sections. The extension to variable 5.2. Linear vibrations
sections, anisotropic beams and static and dynamic mul-
timodal analyses of beams will be investigated in the The vibration analyses are investigated by numeri-
next work. The critical buckling loads and the natural cally solving Eq. (20) and a large number of vibration
frequencies and the corresponding eigenmodes are deter- modes can be obtained. The fifth first free vibration
mined by solving the resulting eigenvalue problem at modes are presented in Figs. 5 and 6 for C–S and C–F
concatenation points. The deflection, the slope, the beams on concentrated elastic foundations at Xc = 0.75
bending moment and the shear force can be also investi- and at Xc = 1 respectively. The convergence of fifth five
gated at interior and boundaries of the beam. A large natural frequencies with respect to the number of inte-
number of natural frequencies or buckling loads and rior points is presented in Tables 4 and 5. Again, a few
associated eigenmodes can be numerically computed. number of interior points is enough for a good predic-
The accuracy of the result depends on the number of tion of eigenfrequencies. For beams submitted to an
interiors points chosen. Tests of convergence of the axial compression and on an elastic foundation, the
buckling and vibration analyses are investigated. Many natural frequency ratio x2 =x2
1 according to the axial

Table 1
The fifth first buckling loads of a C–S beam obtained by the present model for various numbers of internal points n
Order k Analytical solution [13]
n = 10 n = 20 n = 40 n = 60 n = 80 n = 100
1 20.47 20.27 20.21 20.20 20.20 20.19 20.19
2 62.17 60.36 59.86 59.76 59.72 59.71 59.68
3 128.90 121.59 119.60 119.22 119.08 119.02 118.90
4 225.91 205.36 199.81 198.74 198.36 198.18 197.90
5 360.00 313.51 300.94 298.53 297.67 297.27 296.60
2640 Z. Elfelsoufi, L. Azrar / Computers and Structures 83 (2005) 2632–2649

Table 2
The fifth first buckling loads of a C–F beam obtained by the present model for various numbers of internal points n
Order k Analytical solution [13]
n = 10 n = 20 n = 40 n = 60 n = 80 n = 100
1 2.46 2.47 2.47 2.47 2.47 2.47 2.47
2 22.85 22.38 22.25 22.23 22.22 22.21 22.21
3 62.99 62.04 61.78 61.73 61.71 61.70 61.69
4 135.37 124.73 121.89 121.35 121.16 121.07 120.90
5 220.07 205.33 201.28 200.50 200.22 200.09 199.86

Table 3
Comparison between critical buckling load ratios obtained by the present model and analytical ones of a simply-supported beam on
concentrated elastic foundations for different positions XL and amplitudes of j (n = 60)
Position XL k (j )/k (j = 0) k (j = 0) = p2

j = 20 j = 40 j = 60 j = 80 j = 100
Analytic Present Analytic Present Analytic Present Analytic Present Analytic Present
[13] model [13] model [13] model [13] model [13] model
0 1 1.0002 1 1.0002 1 1.0002 1 1.0002 1 1.0002
0.1 1.038 1.043 1.076 1.085 1.111 1.125 1.146 1.163 1.179 1.2
0.2 1.137 1.142 1.262 1.273 1.378 1.393 1.483 1.503 1.578 1.602
0.3 1.259 1.263 1.498 1.506 1.714 1.726 1.905 1.921 2.072 2.091
0.4 1.365 1.366 1.714 1.718 2.043 2.049 2.348 2.353 2.612 2.623
0.5 1.408 1.408 1.809 1.810 2.203 2.205 2.590 2.592 2.967 2.971

2.0 5
1 2
1.5 3

1.0

0.5
Wi(x)

0.0

-0.5

-1.0
4
-1.5

-2.0
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
x

Fig. 3. The fifth first buckling mode shapes, mass normalized, of a clamped–simply supported beam on a concentrated foundation
j = 150 at Xc = 0.75 (n = 60).

compression ratio k =k1 is presented in Figs. 7 and 8 any desired interior point in the pre-buckling region.
(x1 ¼ x ðk ¼ 0; j ¼ 0Þ; k1 ¼ k ðj ¼ 0ÞÞ. One can ob- For nonlinear pre-buckling and the post-buckling re-
serve that the first natural frequency vanishes at critical gions, a nonlinear analysis is needed [26,27].
buckling load (divergence). If one adds to axial compres-
sion the elastic foundation j the natural frequencies 5.3. Flutter analysis
change significantly as clearly shown in Figs. 7 and 8.
The load-frequency dependence can be easily obtained For a clamped–free beam loaded by a tangential fol-
for uniform and concentrated elastic foundations at lower force (BeckÕs problem), the flutter phenomenon is
 Z. Elfelsoufi, L. Azrar / Computers and Structures 83 (2005) 2632–2649 2641

2.0
3 4
1.5
5
1.0

Wi(x) 0.5

0.0

-0.5
2

-1.0
1

-1.5
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
x

Fig. 4. The fifth first buckling mode shapes, mass normalized, of a clamped–free beam on a concentrated foundation j = 150 at
Xc = 1 (n = 60).

1.5 4 1
5
2
1.2 3

0.9
0.6
0.3
Wi(x)

0.0
-0.3
-0.6
-0.9
-1.2
-1.5
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
x

Fig. 5. The fifth first vibration mode shapes, mass normalized, of a clamped–simply supported beam on a concentrated foundation at
j = 150 at Xc = 0.75 (n = 60).

2.5

2.0
2
5 3 1
1.5

1.0
Wi(x)

0.5

0.0

-0.5

-1.0 4

-1.5
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
x

Fig. 6. The fifth first vibration mode shapes, mass normalized, of a clamped–free beam on a concentrated foundation at j = 150 at
Xc = 1 (n = 60).
2642 Z. Elfelsoufi, L. Azrar / Computers and Structures 83 (2005) 2632–2649

Table 4
The fifth first natural frequencies of a C–S beam obtained by the present model for various numbers of internal points n
Order x Analytical solution [30]
n = 10 n = 20 n = 40 n = 60 n = 80 n = 100
1 239.61 238.23 237.86 237.78 237.76 237.74 237.81
2 2 571.20 2 516.80 2 501.80 2 498.89 2 497.85 2 497.36 2 497.07
3 11 593.84 11 062.42 10 918.36 10 890.49 10 890.49 10 875.93 10 866.83
4 35 599.54 32 786.83 32 041.05 31 897.70 31 897.70 31 822.94 31 780.09
5 88 198.82 77 668.23 74 944.78 74 425.67 74 425.67 74 155.50 74 000.84

Table 5
The fifth first natural frequencies of a C–F beam obtained by the present model for various numbers of internal points n
Order x Analytical solution [30]
n = 10 n = 20 n = 40 n = 60 n = 80 n = 100
1 12.355 12.360 12.362 12.362 12.362 12.362 12.362
2 489.95 486.74 485.83 485.66 485.60 485.57 485.52
3 3 927.73 3 839.64 3 815.21 3 810.46 3 808.77 3 807.97 3 806.62
4 15 636.15 14 892.46 14 689.09 14 649.68 14 635.65 14 629.09 14 617.45
5 44 875.56 41 257.42 40 284.90 40 097.59 40 030.97 39 999.85 39 943.81

1.8

1.6

1.4

1.2
*2
ω / ω1

1.0 Xc= 0.75

*2

0.8
Xc= 0.5 buckling load
0.6 Xc= 0.25

0.4 Xc= 0

0.2

0.0
0.0 0.5 1.0 1.5 2.0
* *
(λ / λ1)

 
Fig. 7. Frequency ratio (x2 =x2
1 ) with respect to axial compression ratio (k =k1 ) of C–S beam at various positions of a concentrated
elastic foundation j = 50 (n = 60).

investigated. In the present analysis, the coalescence cri- foundation with various amplitudes j . It can be seen
terion is used. The numerically obtained flutter load with that the flutter load is amplitude-foundation indepen-
60 internal points is k = 20.0625 (k = 20.05[19–21]). At dent and the BeckÕs solution [15–21] is obtained. The
this load, the first and the second eigenfrequencies coin- case of concentrated foundations is presented in Fig.
cide (x2 2
1 ¼ x2 ¼ 121:46) and become complex conju- 11 and shows a higher dependence between the flutter
gate after the flutter load. In Fig. 9 are presented the load and the position of the concentrated foundation.
fifth first vibration modes and shows a perfect coinci- For this test, some more information is supplied in order
dence between the first and the second mode. The third, to explain its particular behavior. The flutter load varia-
the fourth and the fifth eigenmodes are real. Fig. 10 tion with respect to the position and the amplitude of the
shows the load-frequency curves for a uniform elastic concentrated elastic foundations is presented in Fig. 12.
 Z. Elfelsoufi, L. Azrar / Computers and Structures 83 (2005) 2632–2649 2643

*2 3
ω / ω1 X c= 0.75
*2

2
buckling load

X c= 0.5
1 X c= 0.25

0
0.0 0.5 1.0 1.5 2.0
(λ* / λ*1)
 
Fig. 8. Frequency ratio (x2 =x2
1 ) with respect to axial compression ratio (k =k1 ) of C–F beam at various positions of a concentrated

elastic foundation j = 20 (n = 60).

2.5
2.0
4 3 5
1.5
1.0

0.5
1 and 2
Wi(x)

0.0
-0.5
-1.0
-1.5
-2.0
-2.5
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
x

Fig. 9. The fifth first vibration mode shapes at flutter load of a clamped–free beam subjected to a tangential follower force at the free
end, kðflutterÞ ¼ 20:0625 (n = 60).

It is clearly shown that for 0 < j < 35 the flutter may to divergence zone according to the amplitude and the
happen at every position Xc and the flutter load varies position of the concentrated elastic foundation is pre-
slowly from BeckÕs solution and the variation increases sented in Fig. 13.
for Xc > 0.5. In the present analysis, for Xc = 1, the flut-
ter limit correspond to jlim ¼ 35, which is favorably 5.4. Nonlinear vibrations
compared to the results given in [23] ðjlim  36Þ. For
j > 35, the flutter load is strongly position dependent The main extension of this work is the develop-
and there is no flutter (divergence) for some positions. ment of the multimodal analyses of post-buckling, non-
For j < 200, the position Xc  0.77 leads to the smallest linear vibration and nonlinear flutter of beams based on
value of the flutter load and the beam may flutter at a the numerically computed eigenmodes. For nonlinear
very low value. For j > 500, the flutter position zone vibration, this study is limited to the 1-mode analysis
is largely reduced and the flutter load increases highly and the nonlinear amplitude-frequency (24) is used.
from BeckÕs solution. The transition from flutter zone For the sake of brevity, only some benchmark cases
2644 Z. Elfelsoufi, L. Azrar / Computers and Structures 83 (2005) 2632–2649

65
60

∗2
ω2,κ /ω 1
55

∗ 2
50
45
40
35
30 *
κ =320
25 *
κ =240
∗2

20
ω 1,κ /ω 1

*
15 κ =160
∗ 2

*
10 κ =80
5 *
κ =0
0
0 4 8 12 16 20
λ*

Fig. 10. Frequency ratios (x2 2 2 2 

1;j =x1 ) and (x2;j =x1 ) with respect to axial tangential follower force k at the free end of a clamped–free
2
beam on various uniform elastic foundations j ðx1 ¼ 12:36Þ (n = 60).

50
∗2
ω2,κ /ω 1

45
∗ 2

40
Xc=0.5
35
30 Xc=0.25
25
20
Xc=1
∗2
ω1,κ /ω 1

15
Xc=0.75
∗ 2

10
Xc=0
5
0
0 5 10 15 20 25 30 35
λ*

Fig. 11. Frequency ratios (x2 2 2 2 

1;j =x1 ) and (x2;j =x1 ) with respect to axial tangential follower force k at the free end of a clamped–free
beam at various positions Xc of concentrated elastic foundations; j = 20 (n = 60).

*
60 κ =inf
55
*
50 * κ =36
κ =500
45 * *
κ =50 κ =40
40
λ*flutter

35 *
κ =300 *
30 κ =30
25 *
20 κ =0
*
κ =60
*
15
*
κ =70
10 κ =100
5 *
κ = 200
0
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
xc

Fig. 12. kflutter for a clamped–free beam according to the position of the concentrated elastic foundation Xc for different j , j = 0–500
and j = inf (n = 60).
 Z. Elfelsoufi, L. Azrar / Computers and Structures 83 (2005) 2632–2649 2645

800

700

ne
zo
e-
600

nc
ge
er
iv
500 Flutter-zone

D
400
κ*

300
Flutter limit
200
* *
100 λ =λ min=2.8
κ*=35
0
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
xc

Fig. 13. Flutter limit and zone with respect to the amplitude and the position of concentrated foundation (n = 60).

Table 6
Frequency ratio (x =xL ) at various amplitudes of the nonlinear free vibration of a C–C beam under various concentrated elastic
foundations at Xc = 0.75 (n = 60)
a1 Analytic [28] (x =xL ) Present model

j =0 j = 0 j = 300 j = 600 j = inf
0.2 1.000899 1.000897 1.000711 1.000675 1.000321
0.4 1.003590 1.003585 1.002843 1.002699 1.001282
1 1.022231 1.022197 1.017636 1.016749 1.007985
1.5 1.049357 1.049284 1.039259 1.037305 1.017879
2 1.086197 1.086070 1.068798 1.065417 1.031570
2.5 1.131801 1.131612 1.105617 1.100507 1.048911
3 1.185159 1.184899 1.149017 1.141930 1.069724
3.5 1.245275 1.244937 1.198283 1.189027 1.093810

Table 7
Frequency ratio (x =xL ) at various amplitudes of the nonlinear free vibration of a C–S beam under various concentrated elastic
foundations at Xc = 0.75 (n = 60)
a1 Analytic [28] (x =xL ) Present model

j =0 j = 0 j = 300 j = 600 j = inf
0.2 1.002001 1.001999 1.000685 1.000594 1.000252
0.4 1.007980 1.007975 1.002738 1.002375 1.001008
1 1.017868 1.017850 1.016989 1.014755 1.006283
1.5 1.048881 1.048849 1.037834 1.032902 1.014082
2 1.106951 1.106884 1.066333 1.057786 1.024902
2.5 1.183471 1.183360 1.101891 1.088944 1.038646
3 1.275125 1.274963 1.143852 1.125856 1.055202
3.5 1.378898 1.378682 1.191538 1.167976 1.074440

are presented such as C–C and C–S beams. The nonlin- resonance curves are presented in Fig. 14 for a C–S
ear frequencies with respect to mode-amplitudes for beam on concentrated elastic foundations at (Xc =
various concentrated elastic foundation amplitudes at 0.75) and subjected to a lateral concentrated har-
Xc = 0.75 are presented in Tables 6 and 7. The rigidity monic excitation at the beam center. Increasing the
effect added by the elastic foundation on the nonlin- foundation amplitude leads to reducing the nonlinear
ear free vibration is analyzed. The backbone and effect.
2646 Z. Elfelsoufi, L. Azrar / Computers and Structures 83 (2005) 2632–2649

4.0

0
=5
00

*
=1
3.5

κ
*

=0
κ

*
κ
3.0

2.5
a1
2.0

1.5

1.0

0.5

0.0 0.5 1.0 1.5 2.0 2.5

* * 2
(ω /ω L)

Fig. 14. Resonance curves of forced vibrations of a C–S beam under harmonic concentrated force at the center of the beam (for
c = 20) for different values of concentrated elastic foundations at Xc = 0.75 (n = 60).

6. Conclusion The fundamental solution W used in this analysis

4 
corresponds to ooxW4 ðx; sÞ ¼ dðx; sÞ and is
A methodological approach based on integral equa-
tion formulations for buckling, flutter and vibration j x  sj3
W  ðx; sÞ ¼
analyses of beams is presented in simple and compact 12
forms. Critical buckling loads, natural frequencies, flut-
ter loads and load-frequency dependences are investi- Several types of radial basis functions fj(x) are tested and
gated for beams on various elastic foundations and the general form is:
boundary conditions. Tests of convergence with respect
to interior points are carried out and showed that 60 fj ðxÞ ¼ 1 þ arj þ br2j þ cr3j where rj ¼j x  xj j
points are largely enough for good accuracies of eigen-
values and eigenmodes. Nonlinear free and forced vibra- The other functions gj and hj are calculate using the fol-
tions of beams based on one mode analysis are lowing expression:
investigated. The presented model is quite general and
all obtained results are in agreement with available data. d4 g j d4 hj
ðxÞ ¼ fj ðxÞ and ðxÞ ¼ K 2 ðxÞfj ðxÞ
Multimodal analyses based on the computed eigen- dx4 dx4
modes for post-buckling, nonlinear vibrations and non-
r4j r5j r6j r7j
linear flutter will be investigated in next work. gj ðxÞ ¼ hj ðxÞ ¼ þa þb þc
24 120 360 840

Acknowledgements where a, b and c are the chosen constants [3,4].

The authors wish to greatly thank the ÔMinistère de

lÕEnseignement Supérieur et de la Recherche Scientifi- Appendix B. Matrices corresponding to specified
queÕ and the center ÔCNRSTÕ of Morocco and the center boundary conditions
ÔCNRSÕ of France for financial support from both pro-
jects: ‘‘PROTARS III (D11/22)’’ and ‘‘Action Intégrée The system of Eq. (15) is put in the following matrix
(MA/05/17)’’. form vectors and matrices used will be specified for each
boundary condition considered
" #( )
Appendix A I A fW g

0 A0 fT g
In this study, we limit ourselves to isotropic elastic
" # " # " #!
beams with a constant section. In this case, the rigidity B C D
2  
and mass functions K1(x) and K2(x) are constant ¼ x þk þK F fW g
(K1(x) = K2(x) = 1). B0 C0 D0
 Z. Elfelsoufi, L. Azrar / Computers and Structures 83 (2005) 2632–2649 2647

Following the notations given in Fig. 2, let us put: C 0 : matrix ð3  ðn þ 2ÞÞ; C 01j ¼ Cð0Þ;
b b
b j ð0Þ; C 03j ¼ C
b j ð1Þ; j ¼ 1 to n þ 2
B-1 S–S beam: W1 = 0, M1 = 0, Wn+2 = 0 and Mn+2 = 0 C 02j ¼ C

jsj j1sj D : matrix ððn þ 1Þ  ðn þ 2ÞÞ; Dij ¼ Dj ðsi Þ;

A1 ðsÞ ¼  h1 ; A2 ðsÞ ¼ hnþ2 ;
2 2
i ¼ 2 to n þ 2 and j ¼ 1 to n þ 2
A3 ðsÞ ¼ W  ðs; 0ÞQ1 ; A4 ðsÞ ¼ W  ðs; 1ÞQnþ2
D0 : matrix ð3  ðn þ 2ÞÞ; D01j ¼ Dj ð0Þ;
fW g ¼ fW 2 ; W 3 ; . . . ; W n ; W nþ1 g and
b
b j ð0Þ; D03j ¼ Db
b j ð1Þ; j ¼ 1 to n þ 2
fT g ¼ fh1 ; hnþ2 ; Q1 ; Qnþ2 g D02j ¼ D

B-2 C–C beam: W1 = 0, h1 = 0, Wn+2 = 0 and hn+2 = 0 F 1 : matrix ððn þ 2Þ  ðn þ 2ÞÞ of radial function matrix


oW  oW  ði  1Þ
A1 ðsÞ ¼  ðs; 0ÞM 1 ; A2 ðsÞ ¼ ðs; 1ÞM nþ2 ; F 1ij ¼ fj ; i ¼ 1 to n þ 2 and j ¼ 1 to n þ 2
ox ox ðn þ 1Þ
A3 ðsÞ ¼ W  ðs; 0ÞQ1 ; A4 ðsÞ ¼ W  ðs; 1ÞQnþ2
F : matrix ððn þ 2Þ  ðn þ 1ÞÞF ij ¼ F 11
iðjþ1Þ ;
fW g ¼ fW 2 ; W 3 ; . . . ; W n ; W nþ1 g and
fT g ¼ fM 1 ; M nþ2 ; Q1 ; Qnþ2 g F iðnþ2Þ ¼ 0; i ¼ 1 to n þ 2 and j ¼ 1 to n þ 1
because for clamped–free boundary conditions
The other matrices are I, O, A, A0, B, B0, C, C0, D, D0,
W1 ¼ 0
F1, F as already defined for a C–S beam.
B-4 C–F beam with a conservative axial force. W1 = 0,
B-3 C–F beam with a follower tangential force: W1 = 0, h1 = 0, Mn+2 = 0 and Qn+2 = j*, Wn+2 + k*hn+2
h1 = 0, Mn+2 = 0 and Qn+2 = j Wn+2
oW 
oW  j1  sj A1 ðsÞ ¼  ðs; 0ÞM 1 ; A2 ðsÞ ¼ W  ðs; 0ÞQ1
A1 ðsÞ ¼  ðs; 0Þ; M 1 ; A2 ðsÞ ¼ hnþ2 ; ox
ox 2 fW g ¼ fW 2 ; W 3 ; . . . ; W nþ1 ; W nþ2 ; hnþ2 g and

A3 ðsÞ ¼ W ðs; 0ÞQ1 fT g ¼ fM 1 ; Q1 g
fW g ¼ fW 2 ; W 3 ; . . . ; W nþ1 ; W nþ2 g and
fT g ¼ fM 1 ; hnþ2 ; Q1 g I : matrixððn þ 2Þ  ðn þ 2ÞÞ;
8
> I i;i ¼ 1; I i;nþ1 ¼ j W  ðs; 1Þ; 0:5;
>
>
I : matrix ððn þ 1Þ  ðn þ 1ÞÞ; I i;i ¼ 1; >
> j1sj
  < I i;nþ2 ¼ 2 ; for i ¼ 1 to n
>
I i;nþ1 ¼ j W ðs; 1Þ  0:5; for i ¼ 1 to n; I nþ1;nþ1 ¼ 0:5
>
>
I nþ1;nþ1 ¼ j W  ð1; 1Þ þ 0:5 the other terms are null >
> I nþ2;nþ1 ¼ j W  ð0; 1Þ  0:5
>
>
:
I nþ2;nþ2 ¼ 0:5
O : matrix ð3  ðn þ 1ÞÞ; O1;nþ1 ¼ j W  ð0; 1Þ  0:5;
and the other terms are null.
o2 W  ðs; 1Þ
O2;nþ1 ¼ j ð0; 1Þ and the other terms are null
os2
o2 W  ðs; 1Þ
O : matrix ð2  ðn þ 2ÞÞ; O1;nþ1 ¼ j ð0; 1Þ
A : matrix ððn þ 1Þ  3Þ; Aik ¼ Ak ðsi Þ; for k ¼ 1–3; os2
i ¼ 2 to n þ 2 and si ¼ ði  1Þ=ðn þ 1Þ and the other terms are null

b A : matrix ððn þ 2Þ  2Þ; Aik ¼ Ak ðsi Þ; for k ¼ 1 to 2;

A0 : matrix ð3  3Þ; A01k ¼ Ak ð0Þ; b k ð0Þ;
A02k ¼ A
i ¼ 2 to n þ 2 and si ¼ ði  1Þ=ðn þ 1Þ
b
b ð1Þ; for k ¼ 1–3
A03k ¼ A
b
b k ð0Þ;
A0 : matrix ð2  2Þ; A01k ¼ A
B : matrix ððn þ 1Þ  ðn þ 2ÞÞ; Bij ¼ Bj ðsi Þ;
b
b k ð1Þ; for k ¼ 1–2
i ¼ 2 to n þ 2 and j ¼ 1 to n þ 2 A02k ¼ A

B0 : matrix ð3  ðn þ 2ÞÞ; B01j ¼ Bj ð0Þ; B : matrix ððn þ 2Þ  ðn þ 2ÞÞ; Bij ¼ Bj ðsi Þ;

b b i ¼ 2 to n þ 2 and j ¼ 1 to n þ 2
b j ð0Þ; B03j ¼ B
B02j ¼ B b j ð1Þ; j ¼ 1 to n þ 2
b
b j ð0Þ;
B0 : matrix ð2  ðn þ 2ÞÞ; B01j ¼ B
C : matrix ððn þ 1Þ  ðn þ 2ÞÞ; C ij ¼ C j ðsi Þ;
b
b j ð1Þ; j ¼ 1 to n þ 2
i ¼ 2 to n þ 2 and j ¼ 1 to n þ 2 B02j ¼ B
2648 Z. Elfelsoufi, L. Azrar / Computers and Structures 83 (2005) 2632–2649

C : matrix ððn þ 2Þ  ðn þ 2ÞÞ; C ij ¼ C j ðsi Þ; [9] Schanz M, Antes H. A boundary integral formulation for
the dynamic behavior of a Timoshenko beam. Adv Bound
i ¼ 2 to n þ 2 and j ¼ 1 to n þ 2
Elem Techn II, Hoggar 2001:475–82.
[10] Evangelos J. Sapountzakis solution of non-uniform torsion
b
b j ð0Þ;
C 0 : matrix ð2  ðn þ 2ÞÞ; C 01j ¼ C of bars by an integral equation method. Comput Struct
b 2000;77:659–67.
b j ð1Þ; j ¼ 1 to n þ 2
C 02j ¼ C [11] Katsikadelis JT, Tsiatas GC. Non-linear dynamic analysis
of beams with variable stiffness. J Sound Vibrat [in press].
D : matrix ððn þ 2Þ  ðn þ 2ÞÞ; Dij ¼ Dj ðsi Þ; [12] El Felsoufi Z., Azrar L. Modélisation des vibrations des
poutres à section variable par la méthode des équations
i ¼ 2 to n þ 2 and j ¼ 1 to n þ 2 intégrales. 6éme congrès de mécanique, Tanger, Maroc,
2003;(1):7–8.
b
b j ð0Þ; [13] Brush Don O, Almroth BO. Buckling of bar, plates and
D0 : matrix ð2  ðn þ 2ÞÞ; D01j ¼ D
shells. New York: McGraw-Hill, Inc.; 1975.
b
b j ð1Þ; [14] Ziegler H. Principle of structural stability. Waltham,
D02j ¼ D j ¼ 1 to n þ 2
MA: Blasdell Publishing Company; 1968.
[15] Bolotin VV. Non conservative problems of the theory of
F 1 : matrix ððn þ 2Þ  ðn þ 2ÞÞ
elastic stability. New York: Pergamon, Press; 1963.


ði  1Þ [16] Herrmann G, Nemat-Nasser S. Energy considerations in
of radial function matrix F 1ij ¼ fj ; the analysis of the stability of non conservative structural
ðn þ 1Þ
systems in dynamic stability of structures. Oxford: Per-
i ¼ 1 to n þ 2 and j ¼ 1 to n þ 2 gamon Press; 1967.
[17] Leipholz H. Stability of elastic systems. The Nether-
F : matrix ððn þ 2Þ  ðn þ 2ÞÞ; F ij ¼ F 11
iðjþ1Þ ;
lands: Sijthoff & Noordhoff; 1980.
[18] Zuo QH, Schreyer HL. Flutter and divergence instability
i ¼ 1 to n þ 2 and j ¼ 1 to n þ 1 of non conservative beams and plates. Int J Solids Struct
1996;33(9):1355–67.
because for clamped–free boundary conditions
[19] Lee HP. Divergence and flutter of a cantilever rod with an
W 1 ¼ 0. The other terms are null. intermediate spring support. Int J Solids Struct 1995;
32(10):1371–82.
[20] Lee HP. Damping effects on the dynamic stability of a rod
subjected to intermediate follower loads. Comput Methods
References Appl Mech Eng 1996;131(1–2):147–57.
[21] Lee HP. Flutter of a cantilever rod with a relocatable
[1] Nardini D, Brebbia CA. A new approach to free vibra- lumped mass. Comput Methods Appl Mech Eng
tion analysis using boundary elements. In: Brebbia CA, 1997;144(1–2):23–31.
editor. Boundary element method in engineering. [22] Lee SY, Hsu KC. Elastic instability of beams subjected to a
Berlin: Springer; 1982. p. 313–26. partially tangential force. J Sound Vibrat 1995;186(1):
[2] Beskos DE. Boundary element methods in dynamic anal- 111–23.
ysis. Part II (1986–1996). Appl Mech Rev ASME 1997; [23] Wang Q, Quek ST. Enhancing flutter and buckling
50:149–97. capacity of column by piezoelectric layers. Int J Solids
[3] Partridge PW, Brebbia CA, Wrobel LW. The dual Struct 2002;39(16):4167–80.
reciprocity boundary element method. Southamp- [24] Gasparini AM, Saetta AV, Vitaliani RV. On the stability
ton: Computational Mechanics Publications; 1992. and instability regions of non-conservative continuous
[4] Kamiya N, Andoh E, Nogae K. Eigenvalue analysis by the system under partially follower forces. Comput Methods
boundary element method: new developments. Eng Anal Appl Mech Eng 1995;124(1–2):63–78.
Elem 1993;2:151–62. [25] Si-Ung Ryu, Sugiyama Y. Computational dynamics
[5] Chang JR, Yeih W, Chen JT. Determination of the natural approach to the effect of damping on stability of a
frequencies and natural mode of a rod using the dual BEM cantilevered column subjected to a follower force. Com-
in conjunction with the domain partition technique. puters & Structures 2003;81(4):265–71.
Computat Mech 1999;24:29–40. [26] Mohri F, Azrar L, Potier-Ferry M. Vibration of post-
[6] Sladek V, Sladek J. Multiple reciprocity method in BEM buckled thin walled elements. Journal of Sound and
formulations for solution of plate bending problems. Eng vibration, In press, 2004.
Anal Bound Elem 1996;17:161–73. [27] Boutyour EH, Azrar L, Potier-Ferry M. Linear vibration
[7] Tanaka M, Chen W. Dual reciprocity BEM applied to of buckled elastic structures with large rotations by an
transient elstodynamic problems with differential quadra- asymptotic numerical method’’. Submitted to Computers
ture method in time. Comput Methods Appl Mech Eng and Structures.
2001;190:2331–47. [28] Azrar L, Benamar R, White RG. A Semi-analytical
[8] Rong G, Kisu H, Huang C. A new algorithm for bending approach to the nonlinear dynamic response problem of
problems of continuous and inhomogeneous beam by the elastic structures at large vibration amplitudes. Part I:
BEM. Adv Eng Softw 1999;30:339–46. General theory and application to the single mode
 Z. Elfelsoufi, L. Azrar / Computers and Structures 83 (2005) 2632–2649 2649

approach to the free and forced vibrations of beams. J the steady state forced periodic response. J Sound Vibrat
Sound Vibrat 1999;224:183–207. 2002;255(1):1–41.
[29] Azrar L, Benamar R, White RG. A Semi-analytical [30] Geradin M, Rixen D. Théorie des vibrations: Application à
approach to the nonlinear dynamic response problem of la dynamique des structures, Masson, 1993.
beams at large amplitudes Part II: Multimode approach to