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

Authors: Si Yuan; Yongliang Wang; Kangsheng Ye
Journal: Mathematical Problems in Engineering
 (2013)

---

## 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.


---

## Body

## 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:(1)Lv≡(EI(x)v′′(x))′′=-λ(P(x)v′(x))′,0<x<l,
subject to the default BCs
(2)v(0)=0,v′(0)=0,v(l)=0,v′(l)=0,
where prime denotes differentiation with respect to the longitudinal coordinate x; l is the member length; EI(x) is the member’s flexural rigidity; P(x) is the distribution function of the axially compressive force; λ is the critical buckling load coefficient (eigenvalue); v(x) is the buckling mode (eigenfunction); and L 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)KD=λKgD,
where D is the so-called mode vector, K is the conventional static stiffness matrix, and Kg 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 (λh,vh), 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 (λh,vh) 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(λ,v) 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 (λh,vh) satisfy
(4)|λ-λh|≤Tol·(1+|λ|)max|v-vh|≤Tol,  with  max|v|=1.


## 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 -λ(P(x)v′)′ being replaced by a given “load” term f(x).It is well known that for elements of degreem (>3), the FE solution vh gains, for sufficiently smooth problems, a super-convergence order of h2(m-1) at element end nodes [9], while the FE solutions vh and its first three derivatives (vh)(k) (k=1,2,3) on element interior points only gain convergence orders of hm+1 and hm+1-k, 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 C1 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 x-a on the element (x-1,x-2) under consideration. For simplicity and conciseness, only the formula for super-convergent displacement [5] is given as follows:
(5)EIava*=EIavah+J20(x-1-x-a)36+J21(x-1-x-a)22-J10(x-2-x-a)36-J11(x-2-x-a)22(6)J1α=∫x-1x-a(f-Lvh)N-2αdx,J2α=∫x-ax-2(f-Lvh)N-1αdx,α=0,1,
where ()* and ()a represent, respectively, super-convergent values and values at the interior point x-a∈(x-1,x-2) and N-iα (i=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 leasthm+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 displacementv* is at least one order higher than that in the ordinary FE solution vh, for elements of degree m>3, v* will be much more accurate than vh when the mesh is fine enough. Thus, a very simple and reasonable strategy for adaptive process is to use v* instead of the exact solution v to estimate the errors in vh 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 vh on π satisfies the maximum norm criterion [12]:
(7)max|v*-vh|≤Tol.
A detailed description and algorithm of the adaptive procedure for linear BVP are given in [12].


## 2.1. EEP Super-Convergence Formula for Linear BVP

Consider the fourth-order BVP of (1) and (2) with -λ(P(x)v′)′ being replaced by a given “load” term f(x).It is well known that for elements of degreem (>3), the FE solution vh gains, for sufficiently smooth problems, a super-convergence order of h2(m-1) at element end nodes [9], while the FE solutions vh and its first three derivatives (vh)(k) (k=1,2,3) on element interior points only gain convergence orders of hm+1 and hm+1-k, 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 C1 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 x-a on the element (x-1,x-2) under consideration. For simplicity and conciseness, only the formula for super-convergent displacement [5] is given as follows:
(5)EIava*=EIavah+J20(x-1-x-a)36+J21(x-1-x-a)22-J10(x-2-x-a)36-J11(x-2-x-a)22(6)J1α=∫x-1x-a(f-Lvh)N-2αdx,J2α=∫x-ax-2(f-Lvh)N-1αdx,α=0,1,
where ()* and ()a represent, respectively, super-convergent values and values at the interior point x-a∈(x-1,x-2) and N-iα (i=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 leasthm+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 displacementv* is at least one order higher than that in the ordinary FE solution vh, for elements of degree m>3, v* will be much more accurate than vh when the mesh is fine enough. Thus, a very simple and reasonable strategy for adaptive process is to use v* instead of the exact solution v to estimate the errors in vh 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 vh on π satisfies the maximum norm criterion [12]:
(7)max|v*-vh|≤Tol.
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 Section2.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 λ(P(x)v′)′. However, it is observed that after the conventional FE solution (λh,vh) has been converged by inverse iteration on a given mesh, further inverse iteration by using (λh,vh) as the load input will reproduce (λh,vh) 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 λ(P(x)v′)′ in (1) with the FE solution (λh,vh) as follows:
(8)Lv=f(x),0<x<l,withf(x)=-λh(P(x)vh′)′v(0)=0,v′(0)=0,v(l)=0,v′(l)=0.
It is evident that the FE solution of the linear BVP (8) on the current mesh is exactly the same as vh itself. In other words, the FE solution vh 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 v* can be calculated by using (5) with f(x) replaced by -λh(P(x)vh′)′, 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
(9)max|v*-vh|≤Tol,withmax|vh|=1.
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 (λh,vh) 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 vh 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 (λh,vh) 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 (λh,vh), the EEP solution v* is calculated using (5) with f(x) replaced by -λh(P(x)vh′)′, and the maximum error on each element maxe|v*-vh| 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 to10-9 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 Figure1 with length l and
(10)EI(x)=EI0(1+xl)4,P(x)=1.
The critical buckling load and corresponding mode are(11a)λ=4π2EI0l2(11b)v(x)=(1+xl)sin2π(1+x/l).To check the accuracy of the mode function, the domain is divided into np=1000 uniform subintervals with the dividing point being xi=il/np (i=0,1,…,np). Then, the error between the calculated vh(x) and the exact v(x) is computed by
(12)εv=maxi=0,1,…,np|v(xi)v(xj)-vh(xi)vh(xj)|with|v(xj)|=maxi=0,1,…,np|v(xi)|.
Noting the introduction of |v(xj)| implies that the normalization is based on the discrete values v(xi) of the exact solution.
The computed mode and the final mesh are shown in Figure2. The computed critical load is λh=39.478417604EI0/l2 which has the relative error ελ=1.39E-14 compared with the exact one, and the corresponding computed mode error is εv=1.09E-12. It is evident that the prespecified error tolerance Tol=10-9 for both the critical load and mode is well satisfied.Figure 1
Simply supported nonuniform beam of Example1.Figure 2
Buckling mode of Example1.Example 2 (buckling of a spindle-shaped beam).
Consider buckling of a simply supported spindle-shaped beam with(13)I(x)=I0[1+4(xl)-4(xl)2],P(x)=1
as shown in Figure 3. The computed critical load is λ=1.85193832π2EI0/l2 which agrees very well with the computed result λ=1.85373576π2EI0/l2 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 Example2.Figure 4
Buckling mode of Example2.Example 3 (buckling of a vertical column bearing uniform load).
Consider a vertical uniform column as shown in Figure5, which is fixed at the bottom and free at the top. The column is subjected to vertical uniform load q. Consider
(14)EI=const,P(x)=1-xH.
The present method gives the critical buckling load qcr=7.83734744EI/H3 which is consistent with the exact solution qcr=7.837EI/H3 [14]. The buckling mode and final mesh are shown in Figure 6.Figure 5
Vertical column bearing uniform load of Example3.Figure 6
Buckling mode of Example3.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 inertiaI is related to the cross-sectional area A by I=γA2 where the constant γ reflects the property of the cross-section shape. The column has volume V, length l, and Young’s modulus E and is subjected to an axially compressive force, the value of which is λ at the critical buckling load. The dimensionless cross-sectional area α(x) and critical buckling load λ- are defined by
(15)α(x)=A(x)lV,λ-=λl4(EγV2),
where the coordinate x (0≤x≤1) is nondimensional coordinate relative to l.
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-10 as plotted in Figure 7.Figure 7
Optimized clamped-clamped column of Example4.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 Table1, 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.
Mode type
Present method
Reference [2]
ODE solver method [3]
DSM [4]
Symmetrical
52.35625427
52.3563
52.35625427
52.35625427
Antisymmetrical
52.35625424
52.3563
52.35625427
52.35625428Figure 8
Symmetrical buckling mode of Example4.Figure 9
Antisymmetrical buckling mode of Example4.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.



---

*Source: 461832-2013-10-30.xml*