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

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.

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𝜂𝜔₀𝑥̇+𝜀+𝜔²₀𝑥+𝛽𝑥²=𝐹cos(Ω𝑡) (2.1)where ? is the viscous damping parameter, 𝜔₀=(𝜆_cr−𝜆)/𝑚 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 𝜔²₀=0 at 𝜆=𝜆_cr), e is an imperfection parameter, ß is a nonlinearity parameter, F is the magnitude of the externally applied load, and O 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 e 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

𝑥_(eq1)=12−𝜔²₀+𝜔⁴₀−4𝛽𝜀𝛽,𝑥_(eq2)=12−𝜔²₀−𝜔⁴₀−4𝛽𝜀𝛽. (2.2a)(2.2b)

For the equilibrium branch described by (2.2a), the eigenvalues are

𝜇₁=−𝜔⁴₀𝜇−4𝛽𝜀,₂=−−𝜔⁴₀−4𝛽𝜀. (2.3)

For the equilibrium solution (2.2b), the eigenvalues are

𝜇₁=𝜔⁴₀𝜇−4𝛽𝜀,₂=−𝜔⁴₀−4𝛽𝜀. (2.4)

Figure 1 shows, for 𝛽>0, the variation of 𝜔²₀ as a function of the equilibrium position 𝑥_eq. Continuous lines correspond to stable equilibrium paths, and dashed lines correspond to unstable paths. A similar figure, symmetric with respect to the 𝜔²₀ axis, is obtained for 𝛽<0. For the perfect system and for the imperfect system when ß and e have opposite signs, there are for any load level two equilibrium positions, a center and a saddle. For the imperfect system, when ß and e have the same sign, there is a region below and above the critical value (𝜔²₀=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 by

(a) β and ε have opposite signs

(b) β and ε have the same sign

Figure 1 

Equilibrium paths of the perfect and imperfect structures (𝛽=1).

𝜆_lim=𝜆_cr−2𝛽𝜀. (2.5)

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^𝑥max_𝑥min2(𝐸_saddle1−𝜀𝑥−2𝜔²₀𝑥²−13𝛽𝑥³)𝑑𝑥, (2.6)where 𝐸_saddle is the total energy of the system at the saddle point, 𝑥_min is the coordinate corresponding to the saddle, given by (2.2b), and 𝑥_max is given by

𝑥_max=2𝜔⁴₀−𝛽𝜀−𝜔²₀2𝛽. (2.7)

The variation of this area with 𝜔²₀ is illustrated in Figure 2. By expanding (2.6) in Taylor series, one obtains as a first approximation

Figure 2 

Variation of the safe region defined by the homoclinic orbit of the saddle corresponding to the unstable equilibrium point. 𝜀=0, 𝛽=−1.

𝐴=6𝜔⁵₀5𝛽²−6𝜔₀𝜀𝛽 (2.8)which shows clearly the influence of the nonlinearity, ß, and imperfection, e, on the safe area.

The variation of the safe region with the load parameter ? and the imperfection e is illustrated in Figure 3. The curve on the ?-e 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

ℎ=(𝜔⁴₀−4𝛽𝜀)³6𝛽². (2.9)

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 area A. The maximum load, 𝜆_max, that can be applied to the structure with a prescribed safe area A is given approximately by [7]

𝜆_max=𝜆_cr−55𝛽𝜀−^1/4𝐴𝛽𝛽𝜀12𝜀. (2.10)

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

𝜇₁=−𝜂𝜔₀+𝜂²𝜔²₀−𝜔⁴₀𝜇−4𝛽𝜀,₂=−𝜂𝜔₀−𝜂²𝜔²₀−𝜔⁴₀−4𝛽𝜀. (2.11)

So, for positive damping, the equilibrium is asymptotically stable.

For solutions (2.2b), the eigenvalues are

𝜇₁=−𝜂𝜔₀+𝜂²𝜔²₀+𝜔⁴₀𝜇−4𝛽𝜀,₂=−𝜂𝜔₀−𝜂²𝜔²₀+𝜔⁴₀−4𝛽𝜀, (2.12)which is a saddle.

The basin of attraction of the lightly damped system is illustrated in Figure 4, 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.

Figure 5 shows the variation of the basin area of the damped (𝜂=0.05) autonomous system with the nonlinear parameter ß and the imperfection magnitude e.

(a) Nonlinearity. 𝜂=0.05, 𝜀=0

(b) Imperfection. 𝜂=0.05, 𝛽=−1

Figure 5 

Variation of the basin area as a function of the nonlinearity and imperfections.

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 Szemplinska-Stupnicka [30], one is interested in the values of F and O 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, 𝐹_𝑒. 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].

Figure 6 shows the escape load, 𝐹_𝑒, as a function of the forcing frequency, O, for ? = 0.05, e = 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 O 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 Figure 7. 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 O = 0.70. After this minimum, the initially stable period one solution undergoes a stable period doubling bifurcation (D₁) and escape occurs just after this solution becomes unstable, as shown in Figure 7(b), for O = 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 O = 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 (D₂). If the bifurcation load D₂ is higher than the bifurcation load D₁, the escape load of the slowly evolving system is then controlled by D₂, as illustrated in Figure 7(d), for O = 1.90. Here, after D₁, escape becomes unpredictable [10–13]. Figure 8 shows a summary of the bifurcation events connected with the escape boundary of the slowly evolving system.

(a) Ω = 0.70

(b) Ω = 1.00

(c) Ω = 1.50

(d) Ω = 1.90

Figure 7 

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; 𝐹𝑒: escape load. η = 0.05, ε = 0, ω 0 = 1, and β = −1.

Figure 8 

Summary of the bifurcations connected with the escape boundary. η = 0.05, ε = 0, ω 0 = 1, and β = −1.

The influence of the geometric imperfection parameter e on the escape load is illustrated in Figure 9 that depicts the stability boundaries for increasing values of e. The stability boundaries show a shift to the lower frequency range as e 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, the basin 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 limit 𝑡→∞. 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 𝑇=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 level F 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 𝐹=𝐹_𝑒. 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.

Figure 10 

Variation of the basin area of the damped forced system. Black: steady-state bounded solutions. White: unbounded solutions. η = 0.05, ε = 0, 𝜔20 = 1, β = −1, and Ω = 0.68.

Figure 11 

Variation of the basin area of the damped forced system. Black and gray: steady-state bounded solutions. White: unbounded solutions. η = 0.5, ε = 0, 𝜔20 = 1, β = −1, and Ω = 0.43.

Figure 12 shows the variation of the basin area parameterized by the basin area of the corresponding unloaded system, 𝐴/𝐴₀, as a function of the force ratio 𝐹/𝐹_𝑒, for three different values of O.

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, 𝐴₀, 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 𝐴/𝐴₀ 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

𝐹_𝑡(𝑡)=𝐹cos(Ω𝑡)+𝐺(𝑡;𝐹,Ω), (4.1)where the random term 𝐺(𝑡;𝐹,Ω) depends on the deterministic parameters F and O.

For the numerical simulation, the following hypotheses about G 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) An ergodic process is a process where the statistics of the random variable 𝐺(𝑡;𝐹,Ω) 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 𝐺(𝑡;𝐹,Ω) is an ergodic process and, consequently, stationary.

Another hypothesis is that G has expected value zero, that is,

𝐸[𝐺(𝑡;𝐹,Ω)]=0. (4.2)

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

Φ_𝐺𝐺𝜎(𝜔)=^2𝐺𝐺2𝜔_𝑙𝜔forΩ−^𝑙2𝜔<𝜔<Ω+^𝑙2, (4.3)where 𝜎²_𝐺𝐺 is the variance of the random force amplitude and 𝜔_𝑙 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)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 parameter a with respect to the deterministic force amplitude F, and the frequency bandwidth 𝜔_𝑙 around the forcing frequency O. 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.05, 𝜔²₀ = 1.00, and ß₀ = -1. For each excitation frequency, five to ten load samples are considered, depending on the dispersion, and the dynamic buckling load, 𝐹_𝑒, 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 𝜔_𝑙 (0.1 and 0.3).

Figure 14 

Influence of the standard deviation parameter of the radom force, 𝜎𝐺𝐺, 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, 𝜔𝑙, on the dynamic buckling load of the structure. Comparison with the dynamic buckling load of the system under deterministic harmonic forcing

The 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: e, ?, 𝜔²₀, and ß. For each control parameter, the following probability density function is assumed:

𝑓(𝛼)=1002𝛼₀𝑄,if𝛼₀−𝛼₀𝑄100<𝛼<𝛼₀+𝛼₀𝑄,1000,otherwise, (4.5)where a is the system parameter, a₀ 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, a₀.

The mean values (design values) considered in the analysis are e₀ = -0.05, ?₀ = 0.05, 𝜔²₀ = 1.00, ß₀ = -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 (𝜔₀=(𝜆_cr−𝜆)/𝑚), 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−𝜆)/𝑚, 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 Figures 14–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.

Figure 17 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 𝐹^𝑀, 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].

Figure 17 

The change in the stable and unstable manifolds as the load F increases.

The prediction of the first crossing of the stable and unstable manifolds can be obtained by the Melnikov 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)where the vector 𝑔=𝑔(𝑥,𝑣,𝑡)=(𝑔₁,𝑔₂), ? is a small parameter, and 𝐸(𝑥,𝑣) 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)

The Melnikov function is given by

𝑀(𝑡₀)=^∞_−∞𝑔^∗⋅∇𝐸(𝑥^∗,𝑣^∗)𝑑𝑡 (5.3)where 𝑔^∗=𝑔(𝑥^∗,𝑣^∗,𝑡+𝑡₀). 𝑥^∗(𝑡) and 𝑣^∗(𝑡) 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,

𝑀(𝑡₀)=0. (5.4)

Equation (5.4) leads to an algebraic expression that can be used to calculate the load level 𝐹^𝑀 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 (e = 0), it is possible to obtain an analytic solution for the homoclinic orbit using the law of conservation of energy [7]. The solution is

𝑥_𝑝𝜔(𝑡)=−²⁰𝛽𝑒^2𝜔0𝑡−4𝑒^𝜔0𝑡+1𝑒^2𝜔0𝑡+2𝑒^𝜔0𝑡+1. (5.5)

The effect of the imperfection parameter e 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)

The coefficients A and L in (5.6) can be obtained by the restrictions

lim_𝑡→∞𝑥(𝑡)=𝑥_saddle,𝑥(0)=𝑥_max. (5.7)

This leads to

𝜀𝐴(𝜀)=−𝜔²₀,𝐿(𝜀)=1−2𝛽𝜀𝜔⁴₀. (5.8)

The time scale coefficient a(e) in (5.6) can be obtained by applying Galerkin method on the residue

1𝑅(𝑡)=2𝑥̇²+12𝜔²₀𝑥²1+𝜀𝑥+3𝛽𝑥³−𝐸_saddle (5.9)and by using the weight function 𝛿𝑥=𝑥(𝑡). This leads to

𝑎(𝜀)=1−𝛽𝜀𝜔⁴₀. (5.10)

Remembering that 𝑣=𝑑𝑥/𝑑𝑡, the Melnikov function becomes

𝑀(𝑡₀)=^∞_−∞𝐹𝑣cos[Ω(𝑡+𝑡₀)]𝑑𝑡−2𝜂𝜔₀^∞_−∞𝑣²𝑑𝑡. (5.11)

In (5.11), 𝐹=𝜉𝐹, 𝜂=𝜉𝜂, and ? is the small perturbation parameter. By substituting x(t) into (5.11) and calculating the integrals, one obtains

𝑀(𝑡₀)=𝐹sin(Ω𝑡₀)[−6𝐿𝜋Ω²𝑎²𝛽sinh(𝜋Ω/𝑎𝜔₀)]−2𝜂𝜔₀[6𝐿𝑎𝜔⁵₀5𝛽²]. (5.12)

Thus, the first crossing occurs when M(𝑡₀) = 0, which leads to

𝐹^𝑀=𝜂𝜋[2𝜔⁶₀5𝛽Ω²−𝜀𝜔²₀Ω²]sinh[𝜋Ω𝜔₀+𝜋Ω𝛽𝜀𝜔⁵₀]. (5.13)

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 𝜔𝑙=0.1 and two values of 𝜎𝐺𝐺.

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 Figure 19, 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 𝐴₀.

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.

Appendices

A. Simulation of the Random Force

In the following, the theoretical fundaments and methodology used to generate the random force in time domain are presented [7].

The idea of an algorithm to generate a stochastic process sample G(t) comes from the expression of the process variance in terms of the spectral density function

𝜎²_𝐺𝐺=^∞_−∞Φ_𝐺𝐺(𝜔)𝑑𝜔 (A.1)

Assuming that the process is ergodic, the variance can also be calculated in time domain as

𝜎²_𝐺𝐺=lim_𝑇0→∞1𝑇₀^𝑇0₀𝑔²1(𝑡)𝑑𝑡≅𝑇₀^𝑇0₀𝑔²(𝑡)𝑑𝑡, (A.2)where 𝑇₀ is the force duration and g(t) is a sample of the stochastic process G(t).

Based on (A.1) and (A.2), the following relation between the time function g(t) and the spectral density function is obtained:

1𝑇₀^𝑇0₀𝑔²(𝑡)𝑑𝑡≅^∞_−∞Φ_𝐺𝐺(𝜔)𝑑𝜔. (A.3)

Discretizing (A.3), one obtains

1𝑁^𝑁−1_𝑚=0𝑔²(𝑚Δ𝑡)≅2^𝑁/2_𝑘=1Φ_𝐺𝐺(𝑘𝜔₀)𝜔₀, (A.4)where Δ𝑡=𝑇₀/𝑁 and Δ𝜔=𝜔₀=2𝜋/𝑇₀.

Parceval theorem [31], which relates the amplitude of a stochastic process in time with the process amplitude on frequency domain, states that

1𝑁^𝑁−1_𝑚=0𝑔²(𝑚Δ𝑡)≅^𝑁−1_𝑘=0|𝐶_𝑔(𝑘Δ𝜔)|², (A.5)where 𝐶_𝑔(𝜔) is the discrete Fourier transform (DFT) coefficient of the process sample g(t).

Substituting (A.5) on the right-hand side of (A.4) and remembering that, for g(t) to be real, it is necessary that 𝐶_𝑔(𝑁/2+𝑖)=𝐶^∗_𝑔(𝑁/2−𝑖), (A.4) can be rewritten as

2^𝑁/2_𝑘=1|𝐶_𝑔(𝑘𝜔₀)|²≅2^𝑁/2_𝑘=1Φ_𝐺𝐺(𝑘𝜔₀)𝜔₀. (A.6)

The above expression is true if

|𝐶_𝑔(𝑘𝜔₀)|=Φ_𝐺𝐺(𝑘𝜔₀)𝜔₀𝑁,𝑘=1,…,2. (A.7)

This expression allows determining the modulus of the coefficients 𝐶_𝑔 of a discrete Fourier transform sample of the stochastic process G(t) in a way that it has a specified spectral density function. Finally, each DFT coefficient of g(t) can be calculated from

𝐶_𝑔(𝑘𝜔₀)=|𝐶_𝑔(𝑘𝜔₀)|cos(𝜃_𝑘)+𝑖|𝐶_𝑔(𝑘𝜔₀)|sin(𝜃_𝑘), (A.8)where the phase angles ?_k are random variables with constant distribution between 0 and 2p. Samples of the random variables can be obtained using a random number generator. In order to use expression (A.8), the following initial values are necessary:

(i)𝑇₀: random process duration,(ii)N: number of points analyzed on the process,(iii) F _GG(?): specified spectral density function.

B. Structural System Liable to Asymmetric Bifurcation

Consider the well-known SDOF structural system shown in Figure 20 comprising an inverted pendulum of length L and mass m, supported laterally by a linear spring of stiffness K in both tension and compression and inclined initially at 45 degrees. The structure is loaded by a vertical dead load of magnitude P (which includes the weight of the mass m). To generate a family of imperfect systems, a small perturbation moment M is applied to the system. This can be caused, for example, by a small load eccentricity, or any other moment generating disturbance. The effects of various imperfections on the response of the model are very similar. The rotation of the inverted pendulum is denoted by x.

: Laterally supported inverted pendulum submitted to an axial force P and a small disturbance M.

The potential energy of the system in terms of the rotation x is given by

𝑉=𝐾𝐿²[1+sin(𝑥)−1]²−𝑃𝐿[1−cos(𝑥)]−𝑀𝑥. (B.1)The kinetic energy of the pendulum is

𝑇=𝑚𝐿²2𝑥̇². (B.2)By expanding the potential energy in Taylor series and retaining all terms up to the third order, one obtains the following nonlinear equilibrium equation:

−𝑀+[𝐾𝐿²2−𝑃𝐿]𝑥−3𝐾𝐿²8𝑥²=0. (B.3)If the load imperfection M is not considered in the analysis, one obtains from the linearized equilibrium equation the following critical load:

𝑃_cr=𝐾𝐿²2. (B.4)The associated equation of motion is given by

𝑚𝐿²𝑥̈−𝑀+𝐿[𝑃_cr−𝑃]𝑥−3𝐾𝐿²8𝑥²=0. (B.5)By introducing the following auxiliary parameters (notice that we use the usual symbols found in literature for load, imperfection, and nonlinearity parameters):

𝜆_cr=𝑃_cr𝐿𝑃,𝜆=𝐿𝑀,𝜀=−𝑚𝐿²,𝛽=−3𝐾8𝑚, (B.6)the following equation of motion is obtained:

𝑥̈+𝜀+𝜔²₀𝑥+𝛽𝑥²=0, (B.7)where

𝜔₀=𝜆_cr−𝜆𝑚. (B.8)The term 𝐾_ef=(𝜆_cr−𝜆) is usually referred to as effective stiffness in the technical literature and the parameter ? is usually referred to as load parameter.