# Influence of Uncertainties on the Dynamic Buckling Loads of Structures Liable to Asymmetric Postbuckling Behavior

Authors: Paulo B. Gonçalves; Donald Mark Santee
Journal: Mathematical Problems in Engineering
 (2008)

---

## Abstract

Structural systems liable to asymmetric bifurcation usually become unstable at static load levels lower than the linear buckling load of the perfect structure. This is mainly due to the imperfections present in real structures. The imperfection sensitivity of structures under static loading is well studied in literature, but little is know on the sensitivity of these structures under dynamic loads. The aim of the present work is to study the behavior of an archetypal model of a harmonically forced structure, which exhibits, under increasing static load, asymmetric bifurcation. First, the integrity of the system under static load is investigated in terms of the evolution of the safe basin of attraction. Then, the stability boundaries of the harmonically excited structure are obtained, considering different loading processes. The bifurcations connected with these boundaries are identified and their influence on the evolution of safe basins is investigated. Then, a parametric analysis is conducted to investigate the influence of uncertainties in system parameters and random perturbations of the forcing on the dynamic buckling load. Finally, a safe lower bound for the buckling load, obtained by the application of the Melnikov criterion, is proposed which compare well with the scatter of buckling loads obtained numerically.


---

## Body

## 1. Introduction

During the past few decades, a considerable effort within the
engineering sciences has been directed towards understanding the behavior of
structures that exhibit unstable postbuckling behavior [1–3]. The main
motivation for this comes from a notorious and persistent discrepancy between
theoretical and experimental results of the buckling loads of several slender
structures, being that the
experimental results are lower than the theoretical ones. A general
explanation for this upsetting behavior is given by Koiter in his pioneering
work on the general theory of buckling and postbuckling behavior of elastic
structures [4]. He showed that imperfections in the geometry or in the load might
decrease substantially the load carrying capacity of these structures under
slow variation of the applied load. This scenario becomes even worse if the
unavoidable uncertainties in system parameters are also taken into account. Elishakoff
[5] and Kounadis [6], among others, studied imperfection-sensitive structures under a step
load. Because the expressions for the critical load are developed from a static
equilibrium analysis, they actually calculate an upper bound for the load
carrying capacity of the real structure, since they do not take into account
the disturbances imposed upon the imperfect structure during its service life
[7]. The influence of these disturbances on the integrity of the structure can
be evaluated by analyzing the evolution of the basin of attraction of the
stable equilibrium configuration as a function of the system parameters. To
take indirectly into account these deleterious effects, lower bounds of
buckling loads have been proposed for design. They are usually based on the
scatter of experimental buckling loads [8, 9]. However, in the past decades
researchers have sought to deduce theoretically well-founded lower bounds for imperfection-sensitive
structures under static load. Croll [8] developed the so-called reduced
stiffness method based on the elimination in the potential energy of the
structure of the energy components mostly eroded by the imperfections. Based on
this idea, reliable lower bounds have been deduced for a series of structures [8, 9].The estimation of the dynamic buckling load of
structures with unstable postbuckling branches—the load
corresponding to escape from the safe prebuckling well—considering the
effects of uncertainties and imperfections is a much more difficult task.
Structures under dynamic loads may exhibit both local and global bifurcations
that affect in different ways the load carrying capacity and degree of safety
of the structure. Global bifurcations are particularly important since they
control, as shown by Thompson et al. [10–13], the evolution of the basins of
attraction of the solutions in phase space. In addition, compared with the
static case, the number of load control parameters is higher. Finally,
experimental results of dynamic buckling loads of slender structures are rather
scarce in literature [14, 15]. Therefore, little is known on the effects of
uncertainties on the load carrying capacity of structures liable to unstable
static buckling. Therefore, the aim of the present work is to shed some light
on this problem by analyzing the behavior of an archetypal model of a
harmonically forced structure liable to asymmetric bifurcation under increasing
static load.First, the evolution of the basin of attraction of the static
equilibrium configuration is studied. Then, the dynamic buckling load under
different loading conditions is evaluated and the different types of
bifurcation connected with the instability boundaries in force control space are
identified. Next, a detailed parametric analysis clarifies the influence of
uncertainties in load and system parameters on the dynamic buckling load. Based
on these results, a lower bound is proposed which compares favorably with the
scatter of buckling loads obtained in the analysis.The evolution of the basins of attraction of these systems is
governed in a large extent by the evolution of the stable and unstable
manifolds of the saddle connected with the hilltop that separates the pre- and postbuckling
wells. Thompson et al. have studied this connection in detail [10–13]. They
show that the erosion and stratification of the basin of attraction increases
significantly after the first crossing of the stable and unstable manifolds.
The load level associated with such event can be obtained by the application of
the Melnikov criterion, which measures the distance between the manifolds [16, 17]. It can be applied to lightly damped system, which is usually the case of
slender structures found in structural engineering. The present work shows that
the zeros of the Melnikov function can be used as a basis for the deduction of
safe lower bounds that can be used effectively in design. Rega and Lenci [18]
discussed recently the use of integrity measures in nonlinear
mechanical oscillators based on the evolution of basins of attraction.A classical example that illustrates asymmetric buckling
behavior in structural system is the plane frame studied by Roorda [1–3, 19].
Recently, Galvão et al. [20] published a detailed parametric analysis of
this frame structure showing the influence of the system parameters on its postbuckling
behavior and imperfection sensitivity. Another system is the perfect shallow
spherical cap under lateral pressure [21].Analyses of the dynamics of structures liable to asymmetric
bifurcations have been studied by, among others, Virgin [22] and Donescu et al.
[23]. General analyses of the static buckling behavior of systems with asymmetric
postbuckling behavior have been conducted recently by Ohaki [24] and Banchio
and Godoy [25].


## 2. Formulation of the Problem

Consider an SDOF system with quadratic nonlinearity that
exhibits under the variation of a static load parameter a transcritical
bifurcation point. The equation of motion of such a system can be written as(2.1)x¨+2ηω0x˙+ε+ω02x+βx2=Fcos(Ωt)where η is
the viscous damping parameter, ω0=(λcr−λ)/m is the natural frequency of the statically
loaded structure in which λcr is a critical load parameter, λ is an
applied load parameter and m is the
mass (this expression describes the load frequency, leading to ω02=0 at λ=λcr), ε is
an imperfection parameter, β is a nonlinearity parameter, F is the magnitude of the externally
applied load, and Ω is the excitation
frequency. The dots indicate derivation with respect to time t.Equation (2.1) is sometimes called Helmholtz
equation [26], meaning a single-well potential with one escape direction. It
applies also to other fields of mechanics such as the rolling of asymmetric vessels, in which case the
critical threshold corresponds to overturning [10]. In fact, (2.1) is the archetypal model
of an asymmetric bifurcation where the parameter ε is the perturbation
responsible for the unfolding [27, 28]. As an example, (2.1) is derived in the appendix
for a structural system liable to asymmetric bifurcation [1–3].


### 2.1. Analysis of the Autonomous System

For the autonomous undamped system, there are two fixed
points. They are(2.2a)x(eq1)=12−ω02+ω04−4βεβ,(2.2b)x(eq2)=12−ω02−ω04−4βεβ.For the equilibrium branch described by (2.2a), the
eigenvalues are(2.3)μ1=−ω04−4βε,μ2=−−ω04−4βε.For the equilibrium solution (2.2b), the eigenvalues are(2.4)μ1=ω04−4βε,μ2=−ω04−4βε.Figure1 shows, for β>0, the variation of ω02 as a function of the equilibrium position xeq.
Continuous lines correspond to stable equilibrium paths, and dashed lines correspond
to unstable paths. A similar figure, symmetric with respect to the ω02 axis, is obtained for β<0. For the perfect
system and for the imperfect system when β and ε have opposite signs, there are
for any load level two equilibrium positions, a center and a saddle. For the
imperfect system, when β and ε have the same sign, there is a region below and
above the critical value (ω02=0 or λ=λcr) where no solution
occurs. This region is bounded by two limit loads corresponding to saddle-node
bifurcations. The limit load defines thus the load carrying capacity of a real imperfect
system. The limit load parameter is given byEquilibrium paths of the perfect and imperfect
structures (β=1).
(a)
β and ε have opposite signs(b)
β and ε have the same sign(2.5)λlim=λcr−2βε.The limit load may be attained only in a slowly evolving,
quasistatic, system. Nonzero initial conditions will further decrease the
buckling load. In fact, the limit load is but an upper bound of the load
carrying capacity of the imperfect system under static loading. The area of the
basin of attraction at the critical point is zero. Therefore, any small
disturbance leads to buckling. A good measure of the integrity and safety of the
system is the area and topology of the basin of attraction [18]. For the
undamped autonomous system, bounded solution only occurs for initial condition
within the area defined by the homoclinic orbit of the saddle connected with
the unstable solution in Figure 1. The area enclosed by the homoclinic orbit is
given by(2.6)A=2∫xminxmax2(Esaddle−εx−12ω02x2−13βx3)dx,where Esaddle is the total energy of the system at
the saddle point, xmin is the coordinate corresponding to the saddle,
given by (2.2b), and xmax is given by(2.7)xmax=2ω04−βε−ω022β.The variation of this area withω02 is illustrated in Figure 2. By expanding (2.6)
in Taylor
series, one obtains as a first approximationFigure 2
Variation of the safe region defined by the
homoclinic orbit of the saddle corresponding to the unstable equilibrium point.ε=0, β=−1.(2.8)A=6ω055β2−6ω0εβwhich shows clearly the influence of
the nonlinearity, β, and imperfection, ε, on the safe area.The variation of the safe region with the load parameterλ
and the imperfection ε is illustrated in Figure 3. The curve on the λ-ε plane
is the so-called imperfection sensitivity curve [1–3].Figure 3
Variation of the safe region as a function of
the load parameterλ and the imperfection ε. β=−1.As the load level approaches the critical value, there is not
only a swift decrease of the safe area but also a decrease in the depth of the
safe potential well,h, which is
given by(2.9)h=(ω04−4βε)36β2.The degree of safety of a given autonomous undamped system
can be established by defining at the stage of design the magnitude of the safe
areaA. The maximum load, λmax, that can be applied to the structure with
a prescribed safe area A is given
approximately by [7](2.10)λmax=λcr−5βε−51/4Aββε12ε.The consideration of viscous damping changes the eigenvalues
but not the equilibrium solutions of (2.1).
For the damped case, the eigenvalues connected with the stable solution (2.2a)
become(2.11)μ1=−ηω0+η2ω02−ω04−4βε,μ2=−ηω0−η2ω02−ω04−4βε.So, for positive damping, the equilibrium is asymptotically
stable.For solutions (2.2b), the eigenvalues are(2.12)μ1=−ηω0+η2ω02+ω04−4βε,μ2=−ηω0−η2ω02+ω04−4βε,which is a saddle.The basin of attraction of the lightly damped system is
illustrated in Figure4, where gray corresponds to bounded solutions and white
to unbounded solutions. Disregarding the infinite tail that corresponds to
initially large amplitude motions, the area of this basin of attraction is only
slightly higher than the area enclosed by the homoclinic orbit. So, (2.6) can be used as a safe measure of the
basin area [7].Figure 4
Basin of attraction of the damped autonomous
system. Gray area: bounded solutions.Figure5 shows the variation of the basin area of the damped (η=0.05) autonomous system with the nonlinear
parameter β and the imperfection magnitude ε.Variation of the basin area as a function of
the nonlinearity and imperfections.
(a)
Nonlinearity.η=0.05, ε=0(b)
Imperfection.η=0.05, β=−1


### 2.2. The Forced System

This section addresses the problem of calculating the dynamic
buckling load of imperfection-sensitive structures under a harmonic load.
Actually, because of the resonance phenomenon this is one of the worst possible
types of dynamic load.The solution set of (2.1) can be classified in bounded and unbounded
solutions. Unbounded solutions indicate ruin of the structure, as its
displacements become increasingly large and incompatible with the structure’s
use and hypotheses embodied into the mathematical modeling. Unbounded solutions
are also called escapesolutions, or simply escape. In this work, as in, for
example, the works of Malasoma et al. [29], Thompson [10], and
Szemplińska-Stupnicka [30], one is interested in the values of F and Ω that
lead to escape from a given potential well. The minimum value of F at which escape occurs, when all other
parameters are maintained fixed, is called the escape load, Fe. The underlying dynamics that ultimately
leads to escape can be very complex. Consequently, the escape boundary, which is the set of escape loads in the parameter
space, is rather involved and can even be of fractal nature [7].Figure6 shows the escape load, Fe, as a function of the forcing frequency, Ω, for η = 0.05, ε = 0, ω0 = 1, and β = −1. Two loading processes are considered:
a suddenly applied harmonic load (dashed line) and a gradually increasing
harmonic load (continuous line). For the suddenly applied harmonic load,
after each load increment, (2.1) is integrated numerically considering zero
initial position and velocity, that is, after each increment the system starts
from rest. For the gradually increasing harmonic load, for a given excitation frequency,
the final position and velocity of the previous load level are taken as the
initial conditions for the current load level (here, a load step of 0.001 is
considered). For comparison purposes, Figure 6 also shows the escape load for a structure under a step
load of infinite duration as well as the static
critical load. As the value of the forcing frequency Ω varies,
there is a series of valleys associated to super harmonics of various orders
culminating with a deep valley around the natural frequency. For higher
excitation frequencies, the escape load increases and can be, due to the
appearance of new attractors, many times larger than the corresponding static
critical load.Figure 6
Stability boundaries in force control space. System
under harmonic and constant load.η = 0.05, ε = 0, ω0 = 1, and β = −1.The escape is connected with a series of local bifurcations,
as illustrated in the bifurcation diagrams depicted in Figure7. In the main resonance
region, for excitation frequencies smaller than those corresponding to the
minimum escape load, escape occurs due to a saddle-node (S-N) bifurcation, as
illustrated in Figure 7(a) for Ω = 0.70. After this minimum, the initially stable
period one
solution undergoes a stable period doubling bifurcation (D1) and
escape occurs just after this solution becomes unstable, as shown in Figure 7(b),
for Ω = 1.00. As the forcing frequency increases, the period doubling bifurcation
becomes unstable and initially escape occurs at this point, as illustrated in
Figure 7(c), for Ω = 1.50. Finally, as the forcing frequency increases even
further, a secondary stable branch appears along the bifurcated path after a
saddle-node bifurcation. This solution also becomes unstable (D2).
If the bifurcation load D2 is higher than the bifurcation load D1,
the escape load of the slowly evolving system is then controlled by D2,
as illustrated in Figure 7(d), for Ω = 1.90. Here, after D1, escape
becomes unpredictable [10–13]. Figure 8 shows a summary of the bifurcation
events connected with the escape boundary of the slowly evolving system.Bifurcation diagrams of the slowly evolving
system for selected values of the excitation frequency, Ω. Bifurcations
connected with the stability boundaries in force control space: S-N:
saddle-node bifurcation;D1: first period doubling bifurcation; D2:
second period doubling bifurcation; Fe:
escape load. η = 0.05, ε = 0, ω0 = 1, and β = −1.
(a)
Ω = 0.70(b)
Ω = 1.00(c)
Ω = 1.50(d)
Ω = 1.90Figure 8
Summary of the bifurcations connected with the
escape boundary.η = 0.05, ε = 0, ω0 = 1, and β = −1.The influence of the geometric imperfection parameterε on
the escape load is illustrated in Figure 9 that depicts the stability
boundaries for increasing values of ε. The stability boundaries show a shift to
the lower frequency range as ε increases. This is due to the decrease in the
natural frequency caused by the imperfection. For comparison, the static
buckling load of the imperfect system is also shown in Figure 9. While in some
frequency ranges the imperfection sensitivity of the escape load is of the same
order or even higher than that of the static case, in other regions the escape
load is almost insensitive to imperfections. However, the escape load of the
slowly evolving system represents only an upper bound of the actual load
bearing capacity of the structure under harmonic loading. Because of dynamic
perturbations, an imperfection-sensitive structure can escape at load levels
much lower than at the
escape load, as will be shown herein through the analysis of the basins of
attraction.Figure 9
Influence of the imperfection parameterε on
the static and dynamic buckling loads.


## 2.1. Analysis of the Autonomous System

For the autonomous undamped system, there are two fixed
points. They are(2.2a)x(eq1)=12−ω02+ω04−4βεβ,(2.2b)x(eq2)=12−ω02−ω04−4βεβ.For the equilibrium branch described by (2.2a), the
eigenvalues are(2.3)μ1=−ω04−4βε,μ2=−−ω04−4βε.For the equilibrium solution (2.2b), the eigenvalues are(2.4)μ1=ω04−4βε,μ2=−ω04−4βε.Figure1 shows, for β>0, the variation of ω02 as a function of the equilibrium position xeq.
Continuous lines correspond to stable equilibrium paths, and dashed lines correspond
to unstable paths. A similar figure, symmetric with respect to the ω02 axis, is obtained for β<0. For the perfect
system and for the imperfect system when β and ε have opposite signs, there are
for any load level two equilibrium positions, a center and a saddle. For the
imperfect system, when β and ε have the same sign, there is a region below and
above the critical value (ω02=0 or λ=λcr) where no solution
occurs. This region is bounded by two limit loads corresponding to saddle-node
bifurcations. The limit load defines thus the load carrying capacity of a real imperfect
system. The limit load parameter is given byEquilibrium paths of the perfect and imperfect
structures (β=1).
(a)
β and ε have opposite signs(b)
β and ε have the same sign(2.5)λlim=λcr−2βε.The limit load may be attained only in a slowly evolving,
quasistatic, system. Nonzero initial conditions will further decrease the
buckling load. In fact, the limit load is but an upper bound of the load
carrying capacity of the imperfect system under static loading. The area of the
basin of attraction at the critical point is zero. Therefore, any small
disturbance leads to buckling. A good measure of the integrity and safety of the
system is the area and topology of the basin of attraction [18]. For the
undamped autonomous system, bounded solution only occurs for initial condition
within the area defined by the homoclinic orbit of the saddle connected with
the unstable solution in Figure 1. The area enclosed by the homoclinic orbit is
given by(2.6)A=2∫xminxmax2(Esaddle−εx−12ω02x2−13βx3)dx,where Esaddle is the total energy of the system at
the saddle point, xmin is the coordinate corresponding to the saddle,
given by (2.2b), and xmax is given by(2.7)xmax=2ω04−βε−ω022β.The variation of this area withω02 is illustrated in Figure 2. By expanding (2.6)
in Taylor
series, one obtains as a first approximationFigure 2
Variation of the safe region defined by the
homoclinic orbit of the saddle corresponding to the unstable equilibrium point.ε=0, β=−1.(2.8)A=6ω055β2−6ω0εβwhich shows clearly the influence of
the nonlinearity, β, and imperfection, ε, on the safe area.The variation of the safe region with the load parameterλ
and the imperfection ε is illustrated in Figure 3. The curve on the λ-ε plane
is the so-called imperfection sensitivity curve [1–3].Figure 3
Variation of the safe region as a function of
the load parameterλ and the imperfection ε. β=−1.As the load level approaches the critical value, there is not
only a swift decrease of the safe area but also a decrease in the depth of the
safe potential well,h, which is
given by(2.9)h=(ω04−4βε)36β2.The degree of safety of a given autonomous undamped system
can be established by defining at the stage of design the magnitude of the safe
areaA. The maximum load, λmax, that can be applied to the structure with
a prescribed safe area A is given
approximately by [7](2.10)λmax=λcr−5βε−51/4Aββε12ε.The consideration of viscous damping changes the eigenvalues
but not the equilibrium solutions of (2.1).
For the damped case, the eigenvalues connected with the stable solution (2.2a)
become(2.11)μ1=−ηω0+η2ω02−ω04−4βε,μ2=−ηω0−η2ω02−ω04−4βε.So, for positive damping, the equilibrium is asymptotically
stable.For solutions (2.2b), the eigenvalues are(2.12)μ1=−ηω0+η2ω02+ω04−4βε,μ2=−ηω0−η2ω02+ω04−4βε,which is a saddle.The basin of attraction of the lightly damped system is
illustrated in Figure4, where gray corresponds to bounded solutions and white
to unbounded solutions. Disregarding the infinite tail that corresponds to
initially large amplitude motions, the area of this basin of attraction is only
slightly higher than the area enclosed by the homoclinic orbit. So, (2.6) can be used as a safe measure of the
basin area [7].Figure 4
Basin of attraction of the damped autonomous
system. Gray area: bounded solutions.Figure5 shows the variation of the basin area of the damped (η=0.05) autonomous system with the nonlinear
parameter β and the imperfection magnitude ε.Variation of the basin area as a function of
the nonlinearity and imperfections.
(a)
Nonlinearity.η=0.05, ε=0(b)
Imperfection.η=0.05, β=−1


## 2.2. The Forced System

This section addresses the problem of calculating the dynamic
buckling load of imperfection-sensitive structures under a harmonic load.
Actually, because of the resonance phenomenon this is one of the worst possible
types of dynamic load.The solution set of (2.1) can be classified in bounded and unbounded
solutions. Unbounded solutions indicate ruin of the structure, as its
displacements become increasingly large and incompatible with the structure’s
use and hypotheses embodied into the mathematical modeling. Unbounded solutions
are also called escapesolutions, or simply escape. In this work, as in, for
example, the works of Malasoma et al. [29], Thompson [10], and
Szemplińska-Stupnicka [30], one is interested in the values of F and Ω that
lead to escape from a given potential well. The minimum value of F at which escape occurs, when all other
parameters are maintained fixed, is called the escape load, Fe. The underlying dynamics that ultimately
leads to escape can be very complex. Consequently, the escape boundary, which is the set of escape loads in the parameter
space, is rather involved and can even be of fractal nature [7].Figure6 shows the escape load, Fe, as a function of the forcing frequency, Ω, for η = 0.05, ε = 0, ω0 = 1, and β = −1. Two loading processes are considered:
a suddenly applied harmonic load (dashed line) and a gradually increasing
harmonic load (continuous line). For the suddenly applied harmonic load,
after each load increment, (2.1) is integrated numerically considering zero
initial position and velocity, that is, after each increment the system starts
from rest. For the gradually increasing harmonic load, for a given excitation frequency,
the final position and velocity of the previous load level are taken as the
initial conditions for the current load level (here, a load step of 0.001 is
considered). For comparison purposes, Figure 6 also shows the escape load for a structure under a step
load of infinite duration as well as the static
critical load. As the value of the forcing frequency Ω varies,
there is a series of valleys associated to super harmonics of various orders
culminating with a deep valley around the natural frequency. For higher
excitation frequencies, the escape load increases and can be, due to the
appearance of new attractors, many times larger than the corresponding static
critical load.Figure 6
Stability boundaries in force control space. System
under harmonic and constant load.η = 0.05, ε = 0, ω0 = 1, and β = −1.The escape is connected with a series of local bifurcations,
as illustrated in the bifurcation diagrams depicted in Figure7. In the main resonance
region, for excitation frequencies smaller than those corresponding to the
minimum escape load, escape occurs due to a saddle-node (S-N) bifurcation, as
illustrated in Figure 7(a) for Ω = 0.70. After this minimum, the initially stable
period one
solution undergoes a stable period doubling bifurcation (D1) and
escape occurs just after this solution becomes unstable, as shown in Figure 7(b),
for Ω = 1.00. As the forcing frequency increases, the period doubling bifurcation
becomes unstable and initially escape occurs at this point, as illustrated in
Figure 7(c), for Ω = 1.50. Finally, as the forcing frequency increases even
further, a secondary stable branch appears along the bifurcated path after a
saddle-node bifurcation. This solution also becomes unstable (D2).
If the bifurcation load D2 is higher than the bifurcation load D1,
the escape load of the slowly evolving system is then controlled by D2,
as illustrated in Figure 7(d), for Ω = 1.90. Here, after D1, escape
becomes unpredictable [10–13]. Figure 8 shows a summary of the bifurcation
events connected with the escape boundary of the slowly evolving system.Bifurcation diagrams of the slowly evolving
system for selected values of the excitation frequency, Ω. Bifurcations
connected with the stability boundaries in force control space: S-N:
saddle-node bifurcation;D1: first period doubling bifurcation; D2:
second period doubling bifurcation; Fe:
escape load. η = 0.05, ε = 0, ω0 = 1, and β = −1.
(a)
Ω = 0.70(b)
Ω = 1.00(c)
Ω = 1.50(d)
Ω = 1.90Figure 8
Summary of the bifurcations connected with the
escape boundary.η = 0.05, ε = 0, ω0 = 1, and β = −1.The influence of the geometric imperfection parameterε on
the escape load is illustrated in Figure 9 that depicts the stability
boundaries for increasing values of ε. The stability boundaries show a shift to
the lower frequency range as ε increases. This is due to the decrease in the
natural frequency caused by the imperfection. For comparison, the static
buckling load of the imperfect system is also shown in Figure 9. While in some
frequency ranges the imperfection sensitivity of the escape load is of the same
order or even higher than that of the static case, in other regions the escape
load is almost insensitive to imperfections. However, the escape load of the
slowly evolving system represents only an upper bound of the actual load
bearing capacity of the structure under harmonic loading. Because of dynamic
perturbations, an imperfection-sensitive structure can escape at load levels
much lower than at the
escape load, as will be shown herein through the analysis of the basins of
attraction.Figure 9
Influence of the imperfection parameterε on
the static and dynamic buckling loads.


## 3. Basin of Attraction of the Forced System and Structural Stability

Mathematically, thebasin
of attraction of a periodic solution is the set of all initial conditions
that lead to a solution (attractor) as time goes to infinity. This means that
if a periodic solution has a large compact basin of attraction, it will be
stable under finite perturbations. On the other hand, if it has a small, or
fragmented, basin of attraction, small finite perturbations can lead the
solution to escape even if the solution is stable. Thus, a measure of the
stability of the structure, in particular its safety, has to be based on a
global view of the behavior of the structure. This global view can be expressed
mathematically by the characteristics of the basin of attraction and its
boundary.The concept of basin of attraction is based on the limitt→∞.
Because of limitations in the numerical integration, a practical concept for
basin of attraction is used. This practical concept is the basin of r-attraction [10–13]. A basin of r-attraction is the set of all
initial conditions that lead to the neighborhood of the respective periodic
solution in r times the forcing
period T=2π/Ω. As the integration time increases, the basin
of r-attraction tends asymptotically
to the basin of attraction. Our experience has shown that a basin of
32-attraction is a reasonable approximation.Numerical explorations have shown that the way the basin of
attraction changes as the load levelF increases can be classified into one of two groups: (a) it gradually decreases
until it vanishes completely; (b) its shape remains the same as the load level
increases, until it suddenly becomes fractal [7]. In these two types of basin
of attraction evolution, the area decreases as the load increases and becomes
zero at F=Fe. The two
types of behavior are illustrated in Figures 10 and 11. One important fact to
note is that, even at the eminence of escape, when the basin of attraction is
very small, the periodic solution represented by the fixed point of the Poincaré
map is a stable solution. This shows that when one uses the stability of the
periodic solution as a measure of the structure’s stability, this value
furnishes only an upper bound to the true stable load.Variation of the basin area of the damped
forced system. Black: steady-state bounded solutions. White: unbounded
solutions.η = 0.05, ε = 0, ω02 = 1, β = −1, and Ω = 0.68.Variation of the basin area of the damped
forced system. Black and gray: steady-state bounded solutions. White: unbounded
solutions.η = 0.5, ε = 0, ω02 = 1, β = −1, and Ω = 0.43.Figure12 shows the variation of the basin area parameterized
by the basin area of the corresponding unloaded system, A/A0, as a function of the force ratio F/Fe, for three different
values of Ω.Figure 12
Variation of the basin area as a function of
the excitation magnitude for selected values of the forcing frequency.For higher and lower values of the excitation frequency, the
variation of the basin area becomes smoother, decreasing the relative length of
the initial plateau. The variation of the basin area, as will be shown in the
next item, is closely related to the sensitivity of the dynamic buckling loads
to perturbations and uncertainties.As a safety
measure, the designer can specify a maximum erosion level with respect to the
initial safe area of the structure,A0,
and determine the corresponding maximum load that can be applied during the
service life of the structure. In the autonomous case, this load level is given
by (2.8). For the system under harmonic excitation, curves of constant A/A0 ratio, obtained numerically,
are depicted in Figure 13 and compared with the previously deduced escape
boundaries. From Figures 12 and 13, one can conclude that the erosion of the
safe basin area varies with the value the forcing frequency. For the excitation
frequency corresponding to the lowest escape load (Ω≅0.82), the safe area remains practically constant
up to the critical value and then drops suddenly to zero.Figure 13
Variation of the basin area as a function of
the load control parameters. Curves of constant parameterized area.


## 4. Influence of Uncertainties in System Parameters on the Dynamic Buckling Load

### 4.1. Nondeterministic Force

The analysis conducted up to this point considered a harmonic
excitation. This is rarely true in practical situations where loads do not lend
themselves to explicit time description, being random or including at least
some kind of noise. So, it is important to know how departures from an ideally
perfect harmonic excitation may affect the performance of the system. Consider
that the applied load is composed of a harmonic deterministic portion plus a
random term such that(4.1)Ft(t)=Fcos(Ωt)+G(t;F,Ω),where the random term G(t;F,Ω) depends on the deterministic parameters F and Ω.For the numerical simulation, the following hypotheses aboutG are adopted in the present work [7, 31].(i) A force that varies randomly in time is mathematically a
stochastic process. A stochastic process is a random variable where the
probability distribution depends on a parameter. If the parameter is
continuous, the process is called continuous. In the present case, this
parameter is time. If the statistics of the process (mean and variance) are
time-independent, the process is called stationary.(ii) Anergodic process is a process where the
statistics of the random variable G(t;F,Ω) are the same as the statistics of only a
sample of the random process taken along time. An ergodic process is always
stationary, but a stationary process may not be ergodic. In this work, it is
assumed that the random term G(t;F,Ω) is an ergodic process and, consequently,
stationary.Another hypothesis is thatG has expected value zero, that is,(4.2)E[G(t;F,Ω)]=0.The description of a stochastic process is usually made in
the frequency domain. Here, it is assumed that the random term has a spectral
density function given by(4.3)ΦGG(ω)=σGG22ωlforΩ−ωl2<ω<Ω+ωl2,where σGG2 is the variance of the random force amplitude
and ωl is the frequency bandwidth.Additionally, it is considered that the standard deviation of
the random force amplitude is proportional to the deterministic force
amplitude, thus(4.4)σGG=aF,where a is the standard deviation parameter. Here, the random force
depends on the frequency and amplitude of the deterministic term.Physically, the random term is a noise that increases with an
increase in the applied force. Another point to be emphasized is that the
random term depends on two prescribed parameters: the standard deviation
parametera with respect to the
deterministic force amplitude F, and
the frequency bandwidth ωl around the forcing frequency Ω. The numerical methodology used to generate
the random force in time domain is presented in the appendix.The influence of the random noise on the dynamic buckling
load is studied considering the following system parameters:η0 = 0.05, ω02 = 1.00, and β0 = −1. For each excitation
frequency, five to ten load samples are considered, depending on the
dispersion, and the dynamic buckling load, Fe,
is computed numerically. The results considering two values of the standard
deviation parameter a (0.1F
and 0.3F) are shown in Figure 14 and compared
with the results obtained for the deterministic load. Figure 15 shows the
results considering two values of the bandwidth ωl (0.1 and 0.3).Figure 14
Influence of the standard deviation parameter
of the radom
force,σGG,
on the dynamic buckling load of the structure. Comparison with the dynamic
buckling load of the system under deterministic harmonic forcing.Figure 15
Influence of the bandwidth parameter of the radom  force,ωl,
on the dynamic buckling load of the structure. Comparison with the dynamic
buckling load of the system under deterministic harmonic forcingThe real values of the systems parameters, such as mass,
damping, and stiffness, are dependent on the quality of the fabrication
process. They can be, and usually are, different from the value assumed at the
design stage. In order to quantify the influence of variations on system
parameters in the vicinity of the design values on the dynamic buckling load, a
parametric analysis is carried out herein.Uncertainties in the following system parameters are
considered:ε, η, ω02, and β. For each control parameter, the following
probability density function is assumed:(4.5)f(α)={1002α0Q,ifα0−α0Q100<α<α0+α0Q100,0,otherwise,where α is
the system parameter, α0 is the mean value of the
chosen parameter, and Q is a parameter
which expresses the quality of the fabrication process as a percentage of the
mean value, α0.The mean values (design values) considered in the analysis
areε0 = −0.05, η0 = 0.05, ω02 = 1.00, β0 = −1, and Q = 10 (10%). In the
parametric analysis, for each excitation frequency, ten samples of the
perturbed parameter are considered and the escape load is computed considering
both a gradually applied harmonic load and a suddenly applied harmonic load, as
in the previous deterministic analysis. An example of such analysis is shown in
Figure 16 where the influence of small variations on the system stiffness
parameter, and, consequently, natural frequency (ω0=(λcr−λ)/m), is considered. Similar distribution is observed when the other
parameters are considered. The scatter of results shows the strong influence of
the stiffness value on the dynamic buckling load and the sensitivity of the
load carrying capacity of the structure to small variations in system
parameters.Figure 16
Influence of uncertainties in the stiffness
parameter,ω0=(λcr−λ)/m, on the dynamic buckling
load. GA load: gradually applied harmonic load. SA load: suddenly applied harmonic
load.It is interesting to notice that the scatter of results
presented in Figures14–16 follows the
pattern of variation of the basin area shown in Figure 13. So, it becomes clear
that there is a close relation between the variation of the basin of attraction
and the scatter of dynamic buckling loads. It is also clear that, in the present
case, the scatter is a function of the excitation frequency. The results show
that the escape load of the perturbed system is mostly lower than the escape
load under deterministic harmonic forcing. One can also observe that the
scatter of buckling loads is highly dependent on the forcing frequency.


## 4.1. Nondeterministic Force

The analysis conducted up to this point considered a harmonic
excitation. This is rarely true in practical situations where loads do not lend
themselves to explicit time description, being random or including at least
some kind of noise. So, it is important to know how departures from an ideally
perfect harmonic excitation may affect the performance of the system. Consider
that the applied load is composed of a harmonic deterministic portion plus a
random term such that(4.1)Ft(t)=Fcos(Ωt)+G(t;F,Ω),where the random term G(t;F,Ω) depends on the deterministic parameters F and Ω.For the numerical simulation, the following hypotheses aboutG are adopted in the present work [7, 31].(i) A force that varies randomly in time is mathematically a
stochastic process. A stochastic process is a random variable where the
probability distribution depends on a parameter. If the parameter is
continuous, the process is called continuous. In the present case, this
parameter is time. If the statistics of the process (mean and variance) are
time-independent, the process is called stationary.(ii) Anergodic process is a process where the
statistics of the random variable G(t;F,Ω) are the same as the statistics of only a
sample of the random process taken along time. An ergodic process is always
stationary, but a stationary process may not be ergodic. In this work, it is
assumed that the random term G(t;F,Ω) is an ergodic process and, consequently,
stationary.Another hypothesis is thatG has expected value zero, that is,(4.2)E[G(t;F,Ω)]=0.The description of a stochastic process is usually made in
the frequency domain. Here, it is assumed that the random term has a spectral
density function given by(4.3)ΦGG(ω)=σGG22ωlforΩ−ωl2<ω<Ω+ωl2,where σGG2 is the variance of the random force amplitude
and ωl is the frequency bandwidth.Additionally, it is considered that the standard deviation of
the random force amplitude is proportional to the deterministic force
amplitude, thus(4.4)σGG=aF,where a is the standard deviation parameter. Here, the random force
depends on the frequency and amplitude of the deterministic term.Physically, the random term is a noise that increases with an
increase in the applied force. Another point to be emphasized is that the
random term depends on two prescribed parameters: the standard deviation
parametera with respect to the
deterministic force amplitude F, and
the frequency bandwidth ωl around the forcing frequency Ω. The numerical methodology used to generate
the random force in time domain is presented in the appendix.The influence of the random noise on the dynamic buckling
load is studied considering the following system parameters:η0 = 0.05, ω02 = 1.00, and β0 = −1. For each excitation
frequency, five to ten load samples are considered, depending on the
dispersion, and the dynamic buckling load, Fe,
is computed numerically. The results considering two values of the standard
deviation parameter a (0.1F
and 0.3F) are shown in Figure 14 and compared
with the results obtained for the deterministic load. Figure 15 shows the
results considering two values of the bandwidth ωl (0.1 and 0.3).Figure 14
Influence of the standard deviation parameter
of the radom
force,σGG,
on the dynamic buckling load of the structure. Comparison with the dynamic
buckling load of the system under deterministic harmonic forcing.Figure 15
Influence of the bandwidth parameter of the radom  force,ωl,
on the dynamic buckling load of the structure. Comparison with the dynamic
buckling load of the system under deterministic harmonic forcingThe real values of the systems parameters, such as mass,
damping, and stiffness, are dependent on the quality of the fabrication
process. They can be, and usually are, different from the value assumed at the
design stage. In order to quantify the influence of variations on system
parameters in the vicinity of the design values on the dynamic buckling load, a
parametric analysis is carried out herein.Uncertainties in the following system parameters are
considered:ε, η, ω02, and β. For each control parameter, the following
probability density function is assumed:(4.5)f(α)={1002α0Q,ifα0−α0Q100<α<α0+α0Q100,0,otherwise,where α is
the system parameter, α0 is the mean value of the
chosen parameter, and Q is a parameter
which expresses the quality of the fabrication process as a percentage of the
mean value, α0.The mean values (design values) considered in the analysis
areε0 = −0.05, η0 = 0.05, ω02 = 1.00, β0 = −1, and Q = 10 (10%). In the
parametric analysis, for each excitation frequency, ten samples of the
perturbed parameter are considered and the escape load is computed considering
both a gradually applied harmonic load and a suddenly applied harmonic load, as
in the previous deterministic analysis. An example of such analysis is shown in
Figure 16 where the influence of small variations on the system stiffness
parameter, and, consequently, natural frequency (ω0=(λcr−λ)/m), is considered. Similar distribution is observed when the other
parameters are considered. The scatter of results shows the strong influence of
the stiffness value on the dynamic buckling load and the sensitivity of the
load carrying capacity of the structure to small variations in system
parameters.Figure 16
Influence of uncertainties in the stiffness
parameter,ω0=(λcr−λ)/m, on the dynamic buckling
load. GA load: gradually applied harmonic load. SA load: suddenly applied harmonic
load.It is interesting to notice that the scatter of results
presented in Figures14–16 follows the
pattern of variation of the basin area shown in Figure 13. So, it becomes clear
that there is a close relation between the variation of the basin of attraction
and the scatter of dynamic buckling loads. It is also clear that, in the present
case, the scatter is a function of the excitation frequency. The results show
that the escape load of the perturbed system is mostly lower than the escape
load under deterministic harmonic forcing. One can also observe that the
scatter of buckling loads is highly dependent on the forcing frequency.


## 5. Stable and Unstable Manifolds: Melnikov Criterion

When an imperfection-sensitive damped structure described by
(2.1) is unloaded (F = 0), the system has
only one stable equilibrium point and one saddle point. In this situation, the
saddle point’s stable manifold goes smoothly around the stable equilibrium
point defining the boundary of the basin of attraction of this point. One of
the branches of the unstable manifold lies inside the safe basin of attraction,
converging to the stable equilibrium solution. The other branch lies outside
the basin of attraction and tends to infinity.Figure17 shows how the stable and unstable manifolds change
as a function of the load. As the load increases, they approach each other and,
at a certain load level FM,
the stable and unstable manifolds cross transversally. When the stable and
unstable manifolds cross transversally at one point, they cross at an infinite
number of points, thus this crossing indicates the beginning of the erosion of
the basin of attraction [16, 17].The change in the stable and unstable manifolds as
the loadF increases.The prediction of the first crossing of the stable and
unstable manifolds can be obtained by theMelnikov
function. This function gives a measure of the distance of the stable and
unstable manifolds, when this distance is small [16, 17]. It applies to problems
where the damping is small, which is usually the case of slender structures, and
when the algebraic expressions for the stable and unstable manifolds for zero
damping are known.When both the damping and the externally applied forces are
small, the vector field of the system can be expressed generically as [16](5.1)x˙=∂E∂v+ξg1,v˙=−∂E∂x+ξg2,where the vector g=g(x,v,t)=(g1,g2), ξ is a
small parameter, and E(x,v) is the total energy of
the unforced, undamped system (ξ = 0). Also admit that g(t) is periodic, that is, it satisfies
the relation(5.2)g(t+T)=g(t).The Melnikov function is given by(5.3)M(t0)=∫−∞∞g*⋅∇E(x*,v*)dtwhere g*=g(x*,v*,t+t0). x*(t) and  v*(t) are the algebraic
expressions of position and velocity of the stable and unstable manifolds of
the conservative system.The two manifolds cross when this distance is zero, that is,(5.4)M(t0)=0.Equation (5.4) leads to an algebraic expression that can be
used to calculate the load level FM at which the tangling of the stable and unstable manifolds first occurs. Next,
we apply the Melnikov method to imperfection-sensitive structures whose motion
is described by (2.1).In the case of structures liable to asymmetric bifurcation
with no imperfection (ε = 0), it is possible to obtain an analytic
solution for the homoclinic orbit using the law of conservation of energy [7].
The solution is(5.5)xp(t)=−ω02βe2ω0t−4eω0t+1e2ω0t+2eω0t+1.The effect of the imperfection parameterε can be
introduced by observing that it does not change qualitatively the solution (the
homoclinic orbit continues to be a homoclinic orbit), but changes only the
position of the center and saddle points. Thus, the approximate solution can be
expressed by(5.6)x(t)=A(ε)+L(ε)xp[a(ε)t].The coefficientsA and L in (5.6) can be obtained by the
restrictions(5.7)limt→∞x(t)=xsaddle,x(0)=xmax.This leads to(5.8)A(ε)=−εω02,L(ε)=1−2βεω04.The time scale coefficienta(ε) in (5.6) can be obtained by applying Galerkin
method on the residue(5.9)R(t)=12x˙2+12ω02x2+εx+13βx3−Esaddleand by using the weight function δx=x(t). This leads to(5.10)a(ε)=1−βεω04.Remembering thatv=dx/dt, the Melnikov function becomes(5.11)M(t0)=∫−∞∞F¯vcos[Ω(t+t0)]dt−2η¯ω0∫−∞∞v2dt.In (5.11), F=ξF¯, η=ξη¯, and ξ is the small perturbation parameter. By
substituting x(t) into (5.11) and calculating the integrals, one obtains(5.12)M(t0)=F¯sin(Ωt0)[−6LπΩ2a2βsinh(πΩ/aω0)]−2η¯ω0[6Laω055β2].Thus, the first crossing occurs whenM(t0) = 0, which
leads to(5.13)FM=ηπ[2ω065βΩ2−εω02Ω2]sinh[πΩω0+πΩβεω05].This expression can thus be used to calculate the load level
above which the structure’s basin of attraction starts to loose its integrity.


## 6. Melnikov Lower Bound

The algebraic expressions for the Melnikov load given by (5.13)
allows the prediction of the load above which the stable solution’s basin of
attraction begins to loose its integrity by the tangling of the stable and
unstable manifolds of the respective saddle point. Figure 18 compares the
Melnikov load with the scatter of the dynamic buckling results shown in Figure
14. As one can observe, all results lie above the Melnikov load, indicating
that (5.13) can be considered as a safe lower bound in the whole range of
excitation frequencies considered in the present work. Reference [7] presents
similar comparisons considering uncertainties in all parameters. In all cases,
the scatter of results, considering reliable deviations from design values, lies
above the present lower bound.Figure 18
Comparison of the scatter of buckling load with
the Melnikov loadωl=0.1  and two values of σGG.The uncertainties in the structure’s parameters (stiffness,
nonlinearity, natural frequency, etc.) and random deviations of the applied
load make real structures become nondeterministic. In this sense, the random
perturbation of the parameters is a numerical simulation of a real structure.
Note that the random perturbation of 10% of the actual design value generates a
large scatter in the values of the escape loads. Despite the large variation of
the escape loads, the Melnikov load is always smaller indicating that it is a
safe lower bound.Finally, in Figure19, the lower bound is compared with the
curves of constant basin area, already shown in Figure 13, and the scatter of
buckling loads obtained considering uncertainties in all system parameters,
except the external load (Q = 10). The results corroborate the lower bound
character of the Melnikov load. However, if a good quality control is
considered at the fabrication stage, the designer may use a less conservative
estimate of the dynamic buckling load based on the safe basin area. In fact,
one can observe in Figure 19 that almost all results in this numerical
experiment are above the curve corresponding to a safe basin with an area equal
to 40% of the reference basin of the unloaded system A0.Figure 19
Comparing Melnikov load with the buckling load of
a realistic imperfection-sensitive structure.


## 7. Conclusions

For a structure liable to asymmetric bifurcation, the
critical load of the perfect or imperfect structure is an upper bound of its
buckling load, since it corresponds to a safe basin with null area. So, any
disturbance, however small, leads to buckling. To preserve the integrity of the
structure, the designer should prescribe a nonzero compact basin surrounding
the fixed point of the desired solution. In this paper, initially, the
integrity of the structure under static load is investigated by the variation
of the safe basin of attraction as a function of the system parameters,
including initial imperfections. The results show that the safe basin decreases
exponentially as one approaches the critical value. Next, the behavior of the
harmonically excited structure is analyzed and the stability boundaries in
force control space are obtained considering different loading histories. The
results show that uncertainties in system parameters or small random perturbations
of the applied load lead to dynamic buckling loads that are mostly lower than
the load of the unperturbed ideal system. The scatter of results varies with
the forcing frequency and is governed by the variation of the safe basin of
attraction. The variation of the safe basin is dictated by the evolution of the
stable and unstable manifolds of the saddle connected with the safe basin
boundary. Melnikov developed a procedure to determine an approximation for the
first crossing of the stable and unstable manifolds of the saddle-point related
to the fundamental stable solution. When the stable and unstable manifolds
cross transversally at one point, they cross transversally at an infinite
number of discrete points. Since the unstable manifold is the fundamental
solution’s basin of attraction boundary, this indicates that the basin of
attractions becomes, at least partially, fractal. Thus, the load level at which
the tangling of the stable and unstable manifolds first occurs can be taken as
the load that marks the beginning of the loss of stability of the structure,
consequently a lower bound for the structure load carrying capacity. The
proposed lower bound, based on a mathematical reasoning that accounts for the
effects of imperfection and dynamical perturbations on the structure, compares
well with the scatter of dynamic buckling loads and can be used as a safe
design recommendation for imperfection-sensitive structures under periodic loads.
Finally, the proposed procedure can be applied to a variety of
imperfection-sensitive structures, in particular structural systems liable to
unstable symmetric or asymmetric bifurcation.



---

*Source: 490137-2008-06-23.xml*