An Adaptive FEM for Buckling Analysis of Nonuniform Bernoulli-Euler Members via the Element Energy Projection Technique

Abstract

This paper presents a new adaptive finite element (FE) procedure for accurate, efficient, and reliable computation of the critical buckling loads and the associated modes of nonuniform Bernoulli-Euler members. After the conventional FE solution on a given mesh has been obtained, a novel conceptual and practical strategy is introduced, in which the FE solution of the eigenproblem is equivalently viewed as the FE solution of an associated linear problem. This strategy allows the recently developed element energy projection (EEP) technique to be readily used to calculate super-convergent FE solutions for buckling modes, which are subsequently used to estimate the FE solution errors and to further guide mesh refinements, until the accuracy of the FE solution matches the user-preset error tolerance. Numerical examples are given to show the effectiveness, efficiency, accuracy, and reliability of the proposed method, which also paves the way for development of an adaptive FE solver for more general and complicated eigenproblems.

1. Introduction

It is important to study buckling problems of slender or thin-walled structures in engineering practice. Although all these structures are practically three-dimensional (3D) elastic bodies, it is convenient and reasonable in many cases to simplify these problems by theories of beams, plates, and shells rather than the theory of elasticity. The procedure presented in this paper is the conventional finite element (FE) method with effective adaptivity for accurate and reliable computation of critical buckling loads (first eigenvalues) and the associated buckling modes (first eigenfunctions) of nonuniform Bernoulli-Euler members, which is an eigenproblem in the fourth-order ordinary differential equation (ODE). Because of the adaptivity capability, the final FE solution from the method is numerically exact in the sense that both the eigenvalue and eigenfunction satisfy the user-specified error tolerances in the strictest manner, that is, maximum norm for eigenfunctions. The paper does not aim at a complete solution of general eigenproblems and, instead, concentrates on the adaptive solution of the first eigenpair of the buckling problems, which is of prime importance in structural engineering. Also for simplicity, fixed-end boundary conditions (BCs) are taken as the default in the following. However, other BCs can be handled by common practice without difficulties.

The buckling loads and modes for flexure of an axially compressed nonuniform Bernoulli-Euler member are calculated by solving the fourth-order ODE eigenproblem: 𝐿𝑣≡𝐸𝐼(𝑥)𝑣^′′(𝑥)^′′=−𝜆𝑃(𝑥)𝑣^′(𝑥)^′,0<𝑥<𝑙, (1) subject to the default BCs 𝑣(0)=0,𝑣^′(0)=0,𝑣(𝑙)=0,𝑣^′(𝑙)=0, (2) where prime denotes differentiation with respect to the longitudinal coordinate 𝑥; 𝑙 is the member length; 𝐸𝐼(𝑥) is the member’s flexural rigidity; 𝑃(𝑥) is the distribution function of the axially compressive force; 𝜆 is the critical buckling load coefficient (eigenvalue); 𝑣(𝑥) is the buckling mode (eigenfunction); and 𝐿 is the associated fourth-order self-adjoint operator.

In conventional FEM, for a given mesh, the standard FE procedure leads to the following linear eigenproblem: 𝐊𝐃=𝜆𝐊_𝑔𝐃, (3) where 𝐃 is the so-called mode vector, 𝐊 is the conventional static stiffness matrix, and 𝐊_𝑔 is the geometric stiffness matrix. By using many available algorithms in FEM, the critical buckling loads and modes can be obtained. In this paper, to obtain a conventional FE solution for the first eigenpair (𝜆^ℎ,𝑣^ℎ), the inverse iteration is used for mode vector computation, and then, the Rayleigh quotient is adopted for eigenvalue estimation [1]. However, after the conventional FE solution (𝜆^ℎ,𝑣^ℎ) has been obtained, its quality and accuracy are unknown and left to be evaluated by the user.

This paper presents an error estimation scheme for buckling modes by using the recently developed element energy projection (EEP) technique [5–8]. Based on the estimated mode error, the FE mesh is successively refined until the accuracy of the solution matches the prespecified error tolerance.

Suppose that the exact solution for the first eigenpair is (𝜆,𝑣) and that the user-preset error tolerance for both buckling load and mode is Tol. The aim of the present method is to find a sufficiently fine mesh 𝜋 on which the FE solutions (𝜆^ℎ,𝑣^ℎ) satisfy 𝜆−𝜆^ℎ≤Tol⋅(1+|𝜆|)max𝑣−𝑣^ℎ≤Tol,withmax|𝑣|=1. (4)

2. Element Energy Projection Technique

This section introduces a recently developed element energy projection (EEP) technique for extracting super-convergent solutions from conventional FE solutions of linear boundary value problems (BVPs).

2.1. EEP Super-Convergence Formula for Linear BVP

Consider the fourth-order BVP of (1) and (2) with −𝜆(𝑃(𝑥)𝑣^′)^′ being replaced by a given “load” term 𝑓(𝑥).

It is well known that for elements of degree 𝑚 (>3), the FE solution 𝑣^ℎ gains, for sufficiently smooth problems, a super-convergence order of ℎ^2(𝑚−1) at element end nodes [9], while the FE solutions 𝑣^ℎ and its first three derivatives (𝑣^ℎ)^(𝑘) (𝑘=1,2,3) on element interior points only gain convergence orders of ℎ^𝑚+1 and ℎ^{𝑚+1−𝑘}, respectively [9].

In earlier study led by the first author [6–8], it has been found that the well-known projection theorem [10] in FE mathematical theory almost holds true for an individual element, with the error introduced being higher order. Based on this and other related studies, the EEP method has been developed for super-convergence computation at the postprocessing stage of FEM. The one-dimensional 𝐶¹ FE analysis, [5], also led by the first author, developed a simplified version of the EEP super-convergence formulae for computation of both super-convergent displacements and derivatives at any interior point 𝑥_𝑎 on the element (𝑥₁,𝑥₂) under consideration. For simplicity and conciseness, only the formula for super-convergent displacement [5] is given as follows: 𝐸𝐼_𝑎𝑣^∗_𝑎=𝐸𝐼_𝑎𝑣^ℎ_𝑎+𝐽⁰₂𝑥₁−𝑥_𝑎³6+𝐽¹₂𝑥₁−𝑥_𝑎²2−𝐽⁰₁𝑥₂−𝑥_𝑎³6−𝐽¹₁𝑥₂−𝑥_𝑎²2𝐽^𝛼₁=^𝑥𝑎_𝑥1𝑓−𝐿𝑣^ℎ𝑁^𝛼₂𝐽d𝑥,_𝛼2=^𝑥2_𝑥𝑎𝑓−𝐿𝑣^ℎ𝑁^𝛼₁d𝑥,𝛼=0,1, (5)(6) where ()^∗ and ()_𝑎 represent, respectively, super-convergent values and values at the interior point 𝑥_𝑎∈(𝑥₁,𝑥₂) and 𝑁^𝛼_𝑖 (𝑖=1,2, 𝛼=0,1) are conventional cubic shape functions of Hermite type.

A large number of numerical experiments showed that the quantity calculated from the above formula gained super-convergence order of at least ℎ^𝑚+2 and the correctness of the formula is also verified by the excellence of all of the results in Section 4. It is also worth mentioning that the calculation of (5) is performed at the postprocessing stage and involves only some definite integration which can easily be achieved numerically, as is commonly done in FE practice.

2.2. Adaptive Scheme for Linear Problems

Since the accuracy of the super-convergent displacement 𝑣^∗ is at least one order higher than that in the ordinary FE solution 𝑣^ℎ, for elements of degree 𝑚>3, 𝑣^∗ will be much more accurate than 𝑣^ℎ when the mesh is fine enough. Thus, a very simple and reasonable strategy for adaptive process is to use 𝑣^∗ instead of the exact solution 𝑣 to estimate the errors in 𝑣^ℎ and further to guide mesh refinements. This estimation tends to become more and more accurate as the mesh becomes finer and finer. The strategy therefore allows ambitious development [11, 12] of a very efficient and reliable adaptive procedure with the goal being, for a user-specified tolerance Tol, to find a proper mesh 𝜋 such that the FE solution 𝑣^ℎ on 𝜋 satisfies the maximum norm criterion [12]: max𝑣^∗−𝑣^ℎ≤Tol. (7) A detailed description and algorithm of the adaptive procedure for linear BVP are given in [12].

3. Adaptive Strategy

Although the adaptive scheme for linear BVPs presented in Section 2.2 is very powerful, versatile, and efficient, it is prohibitive to be directly applied to the buckling problems, which essentially belong to nonlinear problems due to the obviously nonlinear term 𝜆(𝑃(𝑥)𝑣^′)^′. However, it is observed that after the conventional FE solution (𝜆^ℎ,𝑣^ℎ) has been converged by inverse iteration on a given mesh, further inverse iteration by using (𝜆^ℎ,𝑣^ℎ) as the load input will reproduce (𝜆^ℎ,𝑣^ℎ) itself and no further improvement in accuracy will be gained unless the mesh is further refined. Thus, the key technology introduced in this paper is to reduce the eigenproblem in (1) and (2) to a linear problem by replacing the nonlinear term 𝜆(𝑃(𝑥)𝑣^′)^′ in (1) with the FE solution (𝜆^ℎ,𝑣^ℎ) as follows: 𝐿𝑣=𝑓(𝑥),0<𝑥<𝑙,with𝑓(𝑥)=−𝜆^ℎ𝑃(𝑥)𝑣^ℎ′′𝑣(0)=0,𝑣^′(0)=0,𝑣(𝑙)=0,𝑣^′(𝑙)=0. (8) It is evident that the FE solution of the linear BVP (8) on the current mesh is exactly the same as 𝑣^ℎ itself. In other words, the FE solution 𝑣^ℎ of (1) and (2) can also be viewed as an FE solution of the linear BVP (8). Therefore, it is at this stage that the EEP technique can be applied, and hence, the corresponding super-convergent solution 𝑣^∗ can be calculated by using (5) with 𝑓(𝑥) replaced by −𝜆^ℎ(𝑃(𝑥)𝑣^ℎ′)^′, after which error checking and mesh refinement can be implemented as is done for a linear BVP. The only difference is that in eigenproblem the mode needs to be normalized; that is, the error check is conducted as max𝑣^∗−𝑣^ℎ≤Tol,withmax𝑣^ℎ=1. (9) If the error check fails, which means the current mesh is not fine enough, a new refined mesh is generated. Then another round of FE procedure is conducted on the new mesh with a new FE solution (𝜆^ℎ,𝑣^ℎ) obtained. Then, the EEP technique is used again to estimate the error of the new solution followed by adaptive mesh refinement. This process is repeated until a sufficiently fine mesh 𝜋 is found on which the FE solution 𝑣^ℎ fully satisfies the imposed error tolerance Tol as defined in (9).

The adaptive strategy can be summarized in the following three steps.(1)FE Solution. On the current mesh (the initial one is provided by the user), by using the inverse iteration with Rayleigh quotient, the conventional FE (𝜆^ℎ,𝑣^ℎ) is obtained, which, on convergence, can also be viewed as a linear FE solution to the linear problem (8).(2)EEP Solution. For the linear problem (8) and its FE solution (𝜆^ℎ,𝑣^ℎ), the EEP solution 𝑣^∗ is calculated using (5) with 𝑓(𝑥) replaced by −𝜆^ℎ(𝑃(𝑥)𝑣^ℎ′)^′, and the maximum error on each element max_𝑒|𝑣^∗−𝑣^ℎ| is calculated.(3)Mesh Refinement. For those elements in which (9) is not satisfied, the error-averaging method [11] is used to subdivide each into two elements, and consequently, a new refined mesh is obtained, and then return to the first step until all elements satisfy (9).

The above adaptive procedure has been found to work very well, and its performance is illustrated by the practical numerical examples given in the next section.

4. Numerical Examples

The proposed adaptive strategy has been coded into a Fortran 90 program. This section presents four numerical examples showing the excellent performance of the procedure. Throughout, about 14 decimal digits are used for floating point number calculations, and the program is run on a DELL personal computer.

In all the examples, the tolerance Tol is set to 10⁻⁹ and the element degree of polynomials used is 5 throughout. To illustrate the feasibility and accuracy of the proposed procedure, the first example is chosen to be the buckling of a simply supported nonuniform beam whose analytic solution is available, so that the results obtained from the present method are compared with the analytic solution. In the other three examples, calculated results can only be compared, since no analytic solutions are available, with the related research results of other scholars.

Example 1 (buckling of a simply supported nonuniform beam). Consider the simply supported beam in Figure 1 with length 𝑙 and 𝐸𝐼(𝑥)=𝐸𝐼₀𝑥1+𝑙⁴,𝑃(𝑥)=1. (10)
The critical buckling load and corresponding mode are𝜆=4𝜋²𝐸𝐼₀𝑙²𝑣(𝑥)=𝑥1+𝑙sin2𝜋().1+𝑥/𝑙 (11a)(11b)To check the accuracy of the mode function, the domain is divided into 𝑛_𝑝=1000 uniform subintervals with the dividing point being 𝑥_𝑖=𝑖𝑙/𝑛_𝑝 (𝑖=0,1,…,𝑛_𝑝). Then, the error between the calculated 𝑣^ℎ(𝑥) and the exact 𝑣(𝑥) is computed by 𝜀_𝑣=max_{𝑖=0,1,…,𝑛𝑝}𝑣𝑥_𝑖𝑣𝑥_𝑗−𝑣^ℎ𝑥_𝑖𝑣^ℎ𝑥_𝑗with𝑣𝑥_𝑗=max_{𝑖=0,1,…,𝑛𝑝}𝑣𝑥_𝑖. (12) Noting the introduction of |𝑣(𝑥_𝑗)| implies that the normalization is based on the discrete values 𝑣(𝑥_𝑖) of the exact solution.
The computed mode and the final mesh are shown in Figure 2. The computed critical load is 𝜆^ℎ=39.478417604𝐸𝐼₀/𝑙² which has the relative error 𝜀_𝜆=1.39𝐸−14 compared with the exact one, and the corresponding computed mode error is 𝜀_𝑣=1.09𝐸−12. It is evident that the prespecified error tolerance Tol=10⁻⁹ for both the critical load and mode is well satisfied.

Figure 1 

Simply supported nonuniform beam of Example 1.

Figure 2 

Buckling mode of Example 1.

Example 2 (buckling of a spindle-shaped beam). Consider buckling of a simply supported spindle-shaped beam with 𝐼(𝑥)=𝐼₀𝑥1+4𝑙𝑥−4𝑙²,𝑃(𝑥)=1 (13) as shown in Figure 3. The computed critical load is 𝜆=1.85193832𝜋²𝐸𝐼₀/𝑙² which agrees very well with the computed result 𝜆=1.85373576𝜋²𝐸𝐼₀/𝑙² by the Rayleigh-Ritz method in [13]. The corresponding computed mode and final mesh are shown in Figure 4.

Figure 3 

Spindle-shaped beam of Example 2.

Figure 4 

Buckling mode of Example 2.

Example 3 (buckling of a vertical column bearing uniform load). Consider a vertical uniform column as shown in Figure 5, which is fixed at the bottom and free at the top. The column is subjected to vertical uniform load 𝑞. Consider (𝑥)𝑥𝐸𝐼=const,𝑃=1−𝐻. (14) The present method gives the critical buckling load 𝑞_cr=7.83734744𝐸𝐼/𝐻³ which is consistent with the exact solution 𝑞_cr=7.837𝐸𝐼/𝐻³ [14]. The buckling mode and final mesh are shown in Figure 6.

Figure 5 

Vertical column bearing uniform load of Example 3.

Figure 6 

Buckling mode of Example 3.

Example 4 (buckling of an optimized clamped-clamped column). For a clamped-clamped column, the problem of determination of the optimal shape was first considered analytically by Tadjbakhsh and Keller [15], using a single mode formulation with no cross-section constraints. However, Olhoff and Rasmussen [2] discovered that the single mode result was not generally optimal and obtained a better solution based on a bimodal formulation. Their results were confirmed by the first author of the present paper using an ODE-solver method [3]. This buckling problem was further studied by using the exact dynamic stiffness method (DSM) in [4].
For this problem, the variation of the cross section along the column is geometrically similar and similarly oriented so that the moment of inertia 𝐼 is related to the cross-sectional area 𝐴 by 𝐼=𝛾𝐴² where the constant 𝛾 reflects the property of the cross-section shape. The column has volume 𝑉, length 𝑙, and Young’s modulus 𝐸 and is subjected to an axially compressive force, the value of which is 𝜆 at the critical buckling load. The dimensionless cross-sectional area 𝛼(𝑥) and critical buckling load 𝜆 are defined by 𝛼(𝑥)=𝐴(𝑥)𝑙𝑉,𝜆=𝜆𝑙⁴𝐸𝛾𝑉², (15) where the coordinate 𝑥 (0≤𝑥≤1) is nondimensional coordinate relative to 𝑙.
In this example, this problem is revisited by the present method. The optimized shape is calculated in advance by using the ODE-solver method [3] with an error tolerance [16] of 0.5×10⁻¹⁰ as plotted in Figure 7.

Figure 7 

Optimized clamped-clamped column of Example 4.

The first two eigenvalues are coincident in theory with symmetrical and Antisymmetrical modes. To avoid dealing with coincident eigenvalues, the half-structure approach is used to solve the symmetric and antisymmetric modes, respectively. The computed results are listed in Table 1, from which it can be seen that the solutions agree very well with the other results. The first two modal shapes and corresponding final meshes are shown in Figures 8 and 9.

Table 1 

Critical buckling loads (𝜆) of clamped-clamped column of Example 4.

Figure 8 

Symmetrical buckling mode of Example 4.

Figure 9 

Antisymmetrical buckling mode of Example 4.

It is worth mentioning that in the theory of 3D elasticity there is severe stress concentration at the necked cross sections, which reveals the difference from the current beam theory.

5. Conclusions

A new adaptive strategy of finite element (FE) procedure for accurate, efficient, and reliable computation of the critical buckling loads and the associated buckling modes of nonuniform Bernoulli-Euler members is presented. Novel utilization of the element energy projection (EEP) technique yields a simple, efficient, and reliable adaptive FE procedure to find sufficiently fine meshes matching the user-preset error tolerance. The method is simple, general, and reliable and has the potential to be extended to other general eigenproblems.

Acknowledgments

The authors gratefully acknowledge financial support from the National Natural Science Foundation of China (Grant nos. 50678093, 51078199, and 51078198) the Chinese Ministry of Education (Grant no. IRT00736).