ESDEP WG 6

APPLIED STABILITY

To explain the energy methods for assessing critical loads for cases where it is not possible to get a closed-form solution to differential equilibrium equations.

Beam theory

Lecture 6.1: Concepts of Stable and Unstable Elastic Equilibrium

Lecture 6.2: General Criteria for Elastic Stability

Lecture 6.3: Elastic Instability Modes

Lecture 6.5: Iterative Methods for Solving Stability Problems

When certain assumptions are made about the nature of the deformation of an elastic system during the change of configuration associated with the buckling mode, the elastic system may be approximated by one involving suitable and adjustable parameters which are determined in order that the neutral equilibrium conditions are fulfilled. Using this concept practical approximate methods can be derived which are very useful to the design engineer; some of the best known of these methods are presented in this lecture, i.e. the Rayleigh coefficient, the Rayleigh-Ritz method, and the Galerkin method. An outline of some numerical methods, such as the Euler finite difference method and the finite element method, is also given.

Critical stability loads may be determined using either of the following methods:

- by direct solution of the differential equilibrium equations with exact values of critical loads as a result.
- by using approximate methods, which are, for the most part, based on energy methods, and which lead to approximate solutions of buckling problems.

The solution of differential equilibrium equations to satisfy prescribed boundary conditions presents many difficulties and can only be achieved for simple buckling problems for structures with low degrees of freedom; such basic problems are presented and solved, in this way, in Lecture 6.3. This approach, however, will not be considered further in this lecture which instead concentrates on the alternative energy methods mentioned above. It should be noted that powerful iterative methods can also be used to solve stability problems; some of these methods are presented in Lecture 6.5.

When certain assumptions are made about the nature of the deformation of an elastic system during the change of configuration associated with neutral equilibrium (buckling mode), this elastic system may be approximated by one involving suitable and adjustable parameters or generalised coordinates which are to be determined in order that the neutral equilibrium conditions are fulfilled. This idea provides approximate methods which are very useful to the practical engineer; the best known of these methods are presented in this lecture, i.e. the Rayleigh coefficient, the Rayleigh-Ritz method, and the Galerkin method. An outline of some numerical methods, such as the Euler finite difference method and the finite element method, is also given.

If the adjustable parameters mentioned above are judiciously chosen and of adequate number (in the case of an approximation of a continuous system), these approximate methods will give results very close to the exact solution, at the expense of an increased design effort.

The energy methods for the solution of elastic stability problems are based on the principle of existence of a relative minimum of the total potential energy at neutral equilibrium (see Lecture 6.1 and Lecture 6.2).

It is stated that: considering the change DV, of total energy V of a system, from a starting equilibrium configuration to an adjoining configuration defined by an infinitesimally small and kinematically admissible virtual displacement, then the condition of neutral equilibrium is governed by the second variation d^{2}V of the increment DV, that is:

d^{2}V = 0 = minimum (1)

Generally, when considering steel structures, the stability of a structural system under a set of external loads is studied by considering a buckling deformation, f, from a starting configuration, and performing calculations, according to Equation (1), either to check the stability of this starting configuration, or to determine critical values of external loads leading to instability. The deformation f is expressed as a function of one or more independent variables (generally cartesian coordinates); for example, f(x) as the deflection of a beam or f(x,y) as the deflection of a plate. The change of deformation of a system at neutral equilibrium - the buckling deformation or buckling mode - will hereafter be denoted f(X), where X is the coordinate field of the reference system used (one, two or three dimensions). The function f can be continuous if the system is continuous, or can be defined by intervals if the system is discrete.

Solving Equation (1) requires evaluation of the potential function d^{2}V with respect to the function f, inside a domain of integration noted D (length for a beam, area for a plate, for example). Generally d^{2}V involves quadratic and homogeneous quantities with respect to the buckling deformation f(X) and its derivatives, and is a linear function of the external applied loads.

Introducing a common load multiplier, a, for all loading components and defining a reference loading S_{1} (corresponding to a = 1), the loads at any time of a proportional loading history are equal to:

S = a S_{1} (2)

Thus, generally, the function d ^{2}V can be evaluated by:

d^{2}V(f) = F (a, X, f, f´, f² ) dD = 0
(3)

where the function F also involves geometrical and material characteristics of the domain.

Solving Equation (3) with respect to f will lead to determining the values of the loads (via the load multiplier a) at neutral equilibrium, i.e. critical values of loads above which the starting configuration becomes unstable.

The method of calculation of critical loads by Rayleigh coefficient is derived directly from Equation (1); the second variation d^{2}V of the total potential energy may be expressed as follows:

d^{2}V(f) = d^{2}U(f) + d^{2}W ^{ }(f,S)
(4)

where: d^{2}U is the second variation of strain energy (a quadratic and homogeneous function of f); it represents the strain energy corresponding to the buckling mode;

d^{2}W is the second variation of potential energy of external loads S (a linear function of S and a quadratic and homogeneous function of f). It represents the opposite to the work of the external loads corresponding to the buckling mode.

At neutral equilibrium, a (see Equation (2)) takes the particular value a_{cr} which is the critical load multiplier above which instability occurs. Equations (1), (2) and (4) yield:

d^{2}V(f) = d^{2}U(f) + a_{cr} . d^{2}W ^{ }(f,S_{1}) = 0 = minimum
(5)

If the buckling deformation f is known, the critical load multiplier may be obtained from Equation (4); that is:

a_{cr} = -
(6)

Considering now an approximation f_{1} of f (f_{1} different from f), a_{cr} being known. Then, because of the minimum condition:

d^{2}V(f_{1}) = d^{2}U(f_{1}) + a_{cr} . d^{2}W ^{ }(f_{1},S_{1}) > 0
(7)

Equations (5) and (7) yield:

a_{cr} = -
whatever f_{1} ¹ f (8)

that is a_{cr} < min
whatever f_{1} ¹ f (9)

This defines the Rayleigh Principle which states that the critical load multiplier a_{cr}, calculated using Equation (6) with any kinematically admissible displacement different to the true buckling deformation, will give a value of a_{cr} greater than the exact value.

Assuming f_{1} = f + e f_{2}, where f_{2} is any kinematically admissible displacement and e is a constant, this gives:

a_{cr1} = - + 0(e^{2}) = a_{cr} + 0(e^{2})
(10)

where 0(e^{2}) is a quantity in e^{2}. This implies that a first order error in the choice of f gives a second order error in the value of a_{cr}. If calculations are performed, using Equation (6) with a good approximating function f, simply required to satisfy the boundary conditions, a precise value of a_{cr} may be obtained, with an excess error.

The Rayleigh-Ritz method assumes that the exact solution f(X) of the variational problem described by Equation (3) can be approximated by a linear combination of suitably chosen "coordinate functions" f_{1}(X), f_{2}(X), ... f_{n}(X). That is:

f_{n}(X) = q_{1}f_{1}(X) + q_{2}f_{2}(X) + ... + q_{n} f_{n}(X)
(11)

where the q's are constants to be determined, and are to be regarded as generalised coordinates (degrees of freedom) of the system.

When f(X) is approximated by f_{n}(X) in this way, the function d^{2}V, which is to be evaluated, becomes a quadratic and homogeneous function of the q's; Equation (3) can, therefore, be written as:

d^{2}V = {q}^{t} [a] {q} = 0 (12)

where {q} is the vector of the q's and [a] is a matrix whose coefficients a_{ij} are:

a_{ij} = (13)

The coefficients a_{ij} are functions of the load multiplier a and the properties of the system.

Considering the case of non-trivial configurations, {q} 0, Equation (12) for neutral equilibrium requires [a] to be a singular matrix, that is to say that the determinant of [a] must be zero. This condition provides an equation in a, of degree n, whose smallest positive solution is to be regarded as the critical load multiplier a_{cr}.

The functions f_{i} are chosen in advance, depending on the knowledge and assumptions made about the nature of the deformation. They are not unknown and, provided they satisfy the forced (or geometric) boundary conditions for any value of the q's, the choice of "shapes" is arbitrary. It should be noted, however, that the efficiency of the method does depend on a judicious choice for the f's and that it is an advantage to satisfy all the boundary conditions: in practical applications, one will have some idea of the general nature of the true solution f(X) so that the question of using "outrageous" shapes for the f's rarely occurs. If the functions are judiciously chosen, very good accuracy can be attained with relatively few functions. The efficiency of the Rayleigh-Ritz process may be considerably enhanced if, in addition to the forced (or geometric) boundary conditions (concerning translations or rotations at supports, i.e. f and f´), the natural (or mechanical) boundary conditions (concerning curvature, i.e. f² ) are also satisfied.

To get an idea of the accuracy of the results, a more elaborate procedure is used to obtain a sequence of successive approximations; thus, the following may be taken as a first approximation:

f_{1}(X) = q_{1}f_{1}(X) (14)

and as a second approximation:

f_{2}(X) = q_{1}´ f_{1}(X) + q_{2}´ f_{2}(X)
(15)

and successive approximations in a similar way.

Comparison of successive solutions then gives some indication of how accurate the current solution is. It is worth noting that the solution f_{i + 1}(X) will always be better, or at least no worse, than the preceding solution f_{i}(X).

In contrast to the Rayleigh-Ritz method which gives a solution after forming the variational problem, the Galerkin method provides approximate solutions to the differential equations directly, and is applicable whether the transformation into a variational problem is possible or not. It therefore seems to have wider scope than the Rayleigh-Ritz technique and is more attractive in practice since there is no need to evaluate the potential function. It can be demonstrated, however, that the two methods are closely related. The Galerkin method proceeds as follows.

Generally, the governing differential equation of a buckling problem can be written as follows:

L[f(X)] = 0 (16)

where f(X) is the buckling deformation and L stands for a linear and homogeneous differential operator.

Suppose that the exact solution f(X) of Equation (16), is expressible in the form of a complete series of functions:

f(X) = q_{j} f_{j}(X)
(17)

satisfying all the required boundary conditions; the "exactness" of this solution can be expressed by the statement that the left hand side of Equation (16) is orthogonal to every term in the series of Equation (17); that is:

L[f(X)]f_{j}(X) dD = 0
j=1,2,... (18)

Suppose the series of Equation (17) is truncated to a finite number of terms, n, then using the above idea, n conditions of orthogonality may be imposed, as follows:

i=1,2,...n (19)

This can be written as follows because L is a linear operator:

L[f_{j}(X)]f_{i}(X) dD = 0
i=1,2,...n
(19´)

This provides an averaging basis for evaluating the n unknown q's such that:

f_{n}(X) = q_{j} f_{j}(X)
(20)

which constitutes an approximate solution to the differential equation.

The left hand side of Equation (19), which involves properties of the system and external loads via the load multiplier a, is quadratic and homogeneous in q's; this equation can be written in the form of Equation (11) and then treated in the same manner as for the Rayleigh-Ritz method to find the critical loads.

Numerical methods that require the use of a computer may be used to determine the critical loading. An outline of the Euler finite difference method and the Finite Element Method is now given.

Euler Finite Difference Method

In the Rayleigh-Ritz method, it is required that the admissible functions be continuously differentiable throughout the region of integration. The admissible range may be extended by admitting functions having piecewise continuous derivatives. Therefore, the basis of Euler's method of finite differences is to divide the region of integration into a certain number of sub-regions or intervals, assuming linear functions within each one. If f_{i} denotes the value of the function f at the frontier between intervals i and i+1, derivatives of f may be expressed as functions of the f's, and the sum of the second variation of energy over all the intervals is also a function of f's. Here, the f's are to be regarded as the q's in the Rayleigh-Ritz method; Figure 1a illustrates this approach.

Finite Element Method

This method is used especially for solving stability problems for plate and shell structures. The solution of plate-buckling problems by means of Finite Element theory has been increasing in popularity because it makes use of a matrix formulation suitable for computing. The Finite Element Method has the character of a "piecewise" Rayleigh-Ritz technique; the plate is 'cut' into a number of flat elements joined only at specified nodes, and continuity and equilibrium are established at these nodes. A large number of small elements gives a virtually continuous structure, the behaviour of which is similar to a complete plate. Mixed boundary conditions and varying flexural rigidity can be examined without difficulty. Complete structures may be analysed and post-critical behaviour may also be investigated. An example of a plate mesh is shown in Figure 1b.

Strain energy expressions are needed to perform calculations using the various energy methods. Given below are some useful typical strain energy expressions for structural elements, such as members and plates. These expressions denote the change of strain energy corresponding to the buckling deformation.

Members

Notation

L member length

E Young's modulus

G shear modulus

A cross-section area

A_{v} cross-section shear area

I second moment of area

I_{w} warping constant

I_{t} torsion constant

x abscissa along the member (origin at the beginning of the member)

u(x) axial elongation in member at x }

w(x) lateral deflection in member at x } components of

q(x) slope due to curvature alone in member at x } buckling

y(x) shearing angle (dw/dx - q(x)) in member at x } deformation

f(x) torsion angle in member at x }

Strain energy

- Elongation..... d
^{2}U = (21) - Bending......... d
^{2}U = (22) - Torsion......... d
^{2}U = (23) - Shear........... d
^{2}U = (24) - Warping........ d
^{2}U = (25)

Thin Plates

Notation

a plate dimension along x axis

b plate dimension along y axis

t plate thickness

x,y cartesian coordinates of any point (origin at a corner of the plate)

w(x,y) deflection

D plate stiffness = Et^{3}/(12(1-u^{2}))

u Poisson coefficient

Strain Energy

Bending:

d^{2}U =
(26)

Flexural buckling of a compression member is investigated; the critical load is determined using Rayleigh coefficient, Rayleigh-Ritz and Galerkin methods.

The studied compression member is shown in Figure 2: it is pin-ended and the ends are laterally restrained; L, I and E denote the member length, the second moment of area of the section and Young's modulus respectively; the load P acts vertically downwards; the critical value P_{cr} of P is to be determined; the buckling deformation of the member is shown in Figure 3.

Rayleigh Coefficient Method

The following expression is chosen as an approximation of the buckling deflection w(x):

w(x) = a (x^{2} - xL) where a = any non zero constant (27)

which satisfies the boundary conditions w = 0 at x = 0 and x = L. The derivatives are:

dw/dx = 2ax d^{2}w/dx^{2} = 2a (28)

Performing the integration according to Equation (22), the change of strain energy for the buckling deformation is:

d^{2}U = 2a^{2}EIL (29)

The change in the potential energy of P is the opposite of the work done by P for the buckling deformation. The vertical displacement of the point of application of P due to bending deformation is expressed by:

e = (30)

and the change of potential energy of external load yields after integration of Equation (30):

d^{2}W ^{ }= - P e = - P a^{2}L^{3}/6
(31)

The critical load multiplier is obtained by Equation (6), that is:

a_{cr} =
(32)

and, from Equation (2), the critical value P_{cr} is as follows:

P_{cr} = 12 (33)

This value is to be compared with the exact value obtained from the exact buckling deflection:

w(x) = a sin px/L

that is:

P_{cr} = p^{2}
= 9,8696 (34)

This shows that a parabolic buckling deflection is not a very good approximation of the exact buckling mode. If as an approximation the deflection of a simply supported beam under a uniformly distributed load is chosen, that is:

w(a) = a (x^{4} - 2x^{3}L + xL^{3}) (35)

the above calculations give:

P_{cr} = 9,88 (36)

which is very close to the exact value given in Equation (34).

Rayleigh-Ritz Method

To simplify calculations, the origin of the abscissa is set at mid-length of the member (see Figure 4). A buckling deflection is chosen which is a linear combination of the two following coordinate functions:

f_{1}(x) = x^{2} - L^{2}/4
(37)

f_{2}(x) = x^{4} - L^{4}/16
(38)

which both satisfy the boundary conditions w = 0 at x = -L/2 and x = L/2.

Therefore the expression for the buckling deflection is:

w(x) = a f_{1}(x) + bf_{2}(x) = a (x^{2} - L^{2}/4) + b(x^{4} - L^{4}/16)
(39)

The derivatives are:

dw/dx = 2ax + 4bx^{3} d^{2}w/dx^{2} = 2a + 12bx^{2} (40)

The change in strain energy is expressed by:

d^{2}U = 2. dx = EI (2 a^{2}L + b^{2}L^{5} + 2 abL^{3})
(41)

The change in the potential energy of the compression load is expressed by:

d^{2}W ^{ }= -2 . dx = - P
(42)

and Equations (4), (41) and (42) give:

d^{2}V = a^{2} (2EIL - PL^{3}/6) + b^{2} (9EIL^{5}/10 - PL^{7}/56) +
ab (2EIL^{3} - PL^{5}/10) (43)

The required derivatives are:

(44)

and the matrix [a] of Equation (22) is:

[a] = (45)

Its determinant is given by:

Det [a] = (4EIL - PL^{3}/3) (9EIL^{5}/5 - PL^{7}/28) - (2EIL^{3} - PL^{5}/10)^{2} (46)

The smallest positive solution of Det [a] = 0 is:

P_{cr} = 9,875 (47)

which is to be compared to the exact value given in Equation (34). Although the coordinate functions (37) and (38) are individually not really good approximations of the exact buckling mode, their combination (2 degrees of freedom) gives satisfactory results.

Galerkin method

The same approximation of the buckling deformation and co-ordinate sign convention as for the Rayleigh-Ritz method is chosen:

w(x) = af_{1}(x) + bf_{2}(x)
(39)

where f_{1}(x) = x^{2} - L^{2}/4
(37)

f_{2}(x) = x^{4} - L^{4}/16
(38)

which both satisfy the required boundary conditions which are w = 0 at x = -L/2 and x = L/2 (no imposed condition for the end rotations).

It has been seen in Lecture 6.3 that the differential equation governing flexural buckling of a compression member is:

EI + Pw = 0 (48)

The set of equations obtained from Equation (19´) is:

(49)

The required derivatives are:

d^{2}f_{1}/dx^{2} = 2
d^{2}f_{2}/dx^{2} = 12x^{2}

After integration, the following set of equations is obtained:

a.[PL^{5}/30 - EIL^{3}/3] + b.[PL^{7}/105 - EIL^{5}/10]
= 0
(50)

a.[PL^{7}/105 - EIL^{5}/10] + b.[PL^{9}/360 - EIL^{7}/28] = 0

A non-trivial solution exists if the determinant of Equation (50) is equal to zero; that is:

(PL^{5}/30 - EIL^{3}/3) (PL^{9}/360 - EIL^{7}/28) - (PL^{7}/105 - EIL^{5}/10)^{2} = 0 (51)

whose smallest solution is:

P_{cr} = 9,8697 (52)

which is nearly equal to the exact value given in Equation (34).

- Approximate energy methods provide the engineer with practical means to determine critical loadings for most engineering stability problems; these methods make assumptions about the nature of the buckling deformation of the elastic system, involving adjustable parameters which are determined in order to fulfil the neutral equilibrium conditions.
- The Rayleigh coefficient, the Rayleigh-Ritz method and the Galerkin method presented herein are well known methods which generally can be manually applied to simple buckling problems of isolated structural elements under basic loadings.
- As the number of degrees of freedom increases these methods generally require computer analysis as do the Euler finite difference and finite element methods.
- Several other methods of analysis are given in the technical literature; one of these, involving iterative procedures, is described in Lecture 6.5.

- Richards T.H., "Energy Methods in Stress Analysis", Rainbow Bridge Book Company, 1977.
- Mason J., "Variational, Incremental and Energy Methods in Solid Mechanics and Shell Theory", Elsevier Scientific Publishing Company, Amsterdam, Oxford, New York, 1980.
- Langhaar H.L., "Energy Methods in Applied Mechanics", JohnWiley & Sons, New York, London, 1962.
- Timoshenko S., "Theory of Elastic Stability", McGraw Hill Book Company, New York, 1960.
- Massonnet C., "Résistance des Matériaux" - Volume 2, Dunod, Paris, 1963.
- Bleich, F., "Buckling Strength of Metal Structures", McGraw Hill Book Company, New York, 1952.