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

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.


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 userspecified 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: subject to the default BCs 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); V() 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: where D is the so-called mode vector, K is the conventional static stiffness matrix, and K  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 ( ℎ , V ℎ ), 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 ( ℎ , V ℎ ) 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][6][7][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

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).
In earlier study led by the first author [6][7][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  1 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 ( 1 ,  2 ) under consideration.For simplicity and conciseness, only the formula for super-convergent displacement [5] is given as follows: (5) where () * and ()  represent, respectively, super-convergent values and values at the interior point   ∈ ( 1 ,  2 ) 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 superconvergence 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.

Adaptive Scheme for Linear Problems.
Since the accuracy of the super-convergent displacement V * is at least one order higher than that in the ordinary FE solution V ℎ , for elements of degree  > 3, V * will be much more accurate than V ℎ 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 V ℎ 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 V ℎ on  satisfies the maximum norm criterion [12]: A detailed description and algorithm of the adaptive procedure for linear BVP are given in [12].

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 (()V  )  .However, it is observed that after the conventional FE solution ( ℎ , V ℎ ) has been converged by inverse iteration on a given mesh, further inverse iteration by using ( ℎ , V ℎ ) as the load input will reproduce ( ℎ , V ℎ ) 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 (()V  )  in (1) with the FE solution ( ℎ , V ℎ ) as follows: It is evident that the FE solution of the linear BVP (8) on the current mesh is exactly the same as V ℎ itself.In other words, the FE solution V ℎ 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 () replaced by − ℎ (()V ℎ  )  , 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 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 ( ℎ , V ℎ ) 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 V ℎ 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 ( ℎ , V ℎ ) 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 ( ℎ , V ℎ ), the EEP solution V * is calculated using (5) with () replaced by − ℎ (()V ℎ  )

󸀠
, and the maximum error on each element max  |V * − V ℎ | 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.

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 − 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 Figure 1 with length  and The critical buckling load and corresponding mode are 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 V ℎ () and the exact V() is computed by Noting the introduction of |V(  )| implies that the normalization is based on the discrete values V(  ) of the exact solution.
The computed mode and the final mesh are shown in Figure 2. The computed critical load is  ℎ = 39.4784176040 / 2 which has the relative error Mathematical Problems in Engineering   = 1.39 − 14 compared with the exact one, and the corresponding computed mode error is  V = 1.09 − 12.It is evident that the prespecified error tolerance Tol = 10 − 9 for both the critical load and mode is well satisfied.
Example 2 (buckling of a spindle-shaped beam).Consider buckling of a simply supported spindle-shaped beam with as shown in Figure 3.The computed critical load is  = 1.85193832 2  0 / 2 which agrees very well with the computed result  = 1.85373576 2  0 / 2 by the Rayleigh-Ritz method in [13].The corresponding computed mode and final mesh are shown in Figure 4.
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 The present method gives the critical buckling load  cr = 7.83734744/ 3 which is consistent with the exact solution  cr = 7.837/ 3 [14].The buckling mode and final mesh are shown in Figure 6.
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  =  2 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 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 −10 as plotted in Figure 7.
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.
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.

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.