# Accurate Element of Compressive Bar considering the Effect of Displacement

Authors: Lifeng Tang; Jing Xu; Hongzhi Wang; Xinghua Chen
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

By constructing the compressive bar element and developing the stiffness matrix, most issues about the compressive bar can be solved. In this paper, based on second derivative to the equilibrium differential governing equations, the displacement shape functions are got. And then the finite element formula of compressive bar element is developed by using the potential energy principle and analytical shape function. Based on the total potential energy variation principle, the static and geometrical stiffness matrices are proposed, in which the large deformation of compressive bar is considered. To verify the accurate and validity of the analytical trial function element proposed in this paper, a number of the numerical examples are presented. Comparisons show that the proposed element has high calculation efficiency and rapid speed of convergence.


---

## Body

## 1. Introduction

Compressive bar element analysis is essential in structural engineering design. It can be used to calculate the static and stability problem. A large number of theories and methods for compressive bar element analysis have been suggested. The analysis methods include energy method and numerical method.At present, constructing the compressive bar element needs to calculate the stiffness. References [1–3] all agreed that the stiffness matrix of accurate element of static problem for compressive bar could be derived easily. But for the linear stability problem of compressive bar, the accurate solution cannot be obtained based on the common structural mechanics. A more exact solution can be got though the finite element method which was firstly proposed by Clough [4] in an article about plane elastic problem and was perfect in theory as a numerical method. The finite element method includes direct stiffness method, analytical trial function method, or interpolation trial function method. Direct stiffness method can accurately solve differential equation of static but has low efficiency affected by the number of meshes [5]. So Chen [5] proposed refined direct stiffness method evolved from the generalized hybrid element, which could be used in any displacement function (consistent or inconsistent) and had high precision. Interpolation shape function [6] is one of the main methods of numerical differentiation. This method based on interpolation polynomial can be used for the numerical differential formulas derivation. In practice, the element is needed to be divided into several nodes for interpolation, so interpolation shape function has low precision. Analytical trial function method [7–9] which had solved many traditional problems had widely been used in the elastic mechanics. The boundary element method [10], the boundary collocation method [11], and the hybrid element method [12] all utilized the analytical function for analysis. Fu [13] had carried on the system discussion based on the analytical trial function and proposed a series of analytical trial function (ATF) element, which solved many traditional finite element problem successfully. Cen et al. [14] introduced the variational principles of stress function into analytical trial function element and proposed a new idea of high precision and accurate element. The characteristic differential equation [15–17] which is based on the theory of operator matrix provided the basis method for finite element method to construct the analytical trial function.Nowadays the finite element method which performs very well only is for the case of small deformation problems and examples ranges from standard mixed elements [18–20] to enhanced strain elements [21, 22]. But for elements which seem to be ideal from a numerical and a theoretical perspective, it may fail in the large deformation range, due to the high compression states and many interesting methods [23–28] have been developed to solve the problem. However, a satisfactory analysis of finite element methodologies for large deformation problems is still not proposed imperfectly.In the present paper, using the differential equations of equilibrium, and considering the effect of large deformation, the displacement shape function is derived. And then based on the minimum potential energy principle, the elastic stiffness matrixK and geometrical stiffness matrix G are got. Utilizing the stiffness equilibrium equation, we can get the accurate ultimate load value and the buckling mode of the compressive bar subjected to different constraints. Several examples are presented in order to verify the accuracy and validity of the proposed analytical trial function method. Results show that the present method has many advantages including accuracy, efficiency, and simplicity compared to the direct stiffness method and the interpolation shape function element.


## 2. Characteristic of Large Displacement of Compressive Bar

For the static problem of compressive bar, the accurate solution can be obtained by the traditional finite element method. But for the stability issue, when compressive bar is subjected to the critical load, the large deformation is generated, which results in the interaction between the axial load and lateral displacement. Then the lateral displacement derived from the cubic polynomial Hermit interpolation trial function had a huge system error, while by using the direct Hermit interpolation trial function, the generated deformation of compressive bar element is an approximate and the error is also large.The stability problem of compressive bar should be ranked as a geometrical nonlinear problem; in which situation, the large deformation should be considered. Thusw′(x) is not a high order and small value, and the nodal lateral displacement ∫0l(1/2)w′2(x)dx generated by the longitudinal displacement cannot be ignored. Then the nodal lateral and longitudinal displacement are not two independent but interactive variables. Also, because of the influence of large deformation, the differential equilibrium equation between the shear force and bending moment V(x)=M′(x) [29] is not suitable anymore.Based on the characteristic of large displacement of compressive bar, the main key of generating the accurate element is searching for the accurate curve expression of lateral displacement.


## 3. Analytical Trial Function for Displacement of Compressive Bar Element

The steps of creating the displacement analytical shape function of the compressive bar element are as follows:Step 1.
Based on the differential equations of equilibrium, obtain the displacement equation of compressive bar, including the undetermined parameters;Step 2.
Based on the above equation, establish the nodal displacement equation of element.Step 3.
According to the equation of Step2, calculate the undermined parameters.Step 4.
On the basis of the above equations, create the deformation shape function.


### 3.1. Stiffness Balance Equation of Compressive Bar

Figure1 shows the elastic compressive bar, whose length is l, subjected to axial load P, bending moment M, and shear force Q at its two ends. The equilibrium differential governing equation [30] is(1)EId2wxdx2=Pwx-M1-Q1x,where E is elastic modulus; I is the moment of inertia; x denotes the coordinates; w(x) represents the corresponding displacement vector; P is axial load acting on the bar and is taken as negative when it is pressure; M1 is bending moment acting on end 1 of the compressive bar; and Q1 is shear force, as shown in Figure 1.Figure 1
Force diagram map of elastic compressive bar.Based on second derivative, (1) can be changed into(2)EId4wxdx4+λ2d2wxdx2=0,in which λ=P/EI.The solution of (2) is(3)wx=c1+c2x+c3sin⁡λx+c4cos⁡⁡λx=∑i=14fici.


### 3.2. Solution of Displacement Parameters of Element

Equation (3) can also be expressed as(4)wx=∑i=14fici=FC,in which(5)C=C1C2C3C4T,F=f1xf2xf3xf4x=1xsin⁡λxcos⁡⁡λx.By using differential calculation, we easily obtain the expression of the angle caused by the displacement of the compressive bar(6)θx=-dwxdx=-c2-c3λcos⁡⁡λx+c4λsin⁡λx.Define the vector equation of the nodal displacement and corner displacement as(7)δe=w1θ1w2θ2T.Combining (3) with (4) yields the expression of nodal displacement of compressive bar(8)δe=f10f20f30f40-f1′0-f2′0-f3′0-f4′0f1lf2lf3lf4l-f1′l-f2′l-f3′l-f4′lC=AC,in which(9)A=10010-1-λ01lsin⁡λlcos⁡⁡λl0-1-λcos⁡⁡λlλsin⁡λl.Based on (8), the expression of displacement function coefficient is (10)C=A-1δe.The inversion of matrixA is(11)A-1=a11a12a13a14a21a22a23a24a31a32a33a34a41a42a43a44.Each parameter of (11) can be expressed as(12)a11=λlsin⁡λl+cos⁡⁡λl-1R,a12=λlcos⁡⁡λl-sin⁡λlλR,a13=cos⁡⁡λl-1R,a14=--sin⁡λl+λlλR,a21=-λsin⁡λlR,a22=-cos⁡⁡λl-1R,a23=λsin⁡λlR,a24=-cos⁡⁡λl-1R,a31=sin⁡λlR,a32=-λlsin⁡λl+cos⁡⁡λl-1λR,a33=-sin⁡λlR,a34=cos⁡⁡λl-1λR,a41=cos⁡⁡λl-1R,a42=-λlcos⁡⁡λl-sin⁡λlλR,a43=-cos⁡⁡λl-1R,a44=-sin⁡λl+λlλR,where R=-2+λlsinλl+2cos⁡λl.


### 3.3. Analytical Trial Function

Substitution of (10) into (4) yields the element displacement expression of compressive bar(13)wx=FA-1δe=Nδe,in which(14)N=FA-1=N1N2N3N4,where N1, N2, N3, and N4 are shape functions for each displacement parameter, which can be derived by matrix calculation:(15)N1=λl-λx+sin⁡λxsin⁡λl+cos⁡⁡λx+1cos⁡⁡λl-1R,N2=--λlcos⁡⁡λl+λlcos⁡⁡λx-λl+sin⁡λx-sin⁡λlcos⁡⁡λl-1+λxλR,N3=λx-sin⁡λxsin⁡λl+cos⁡⁡λx-11-cos⁡⁡λlR,N4=-λx-sin⁡λxcos⁡⁡λl-1+cos⁡⁡λx-1sin⁡λl-λlλR.


## 3.1. Stiffness Balance Equation of Compressive Bar

Figure1 shows the elastic compressive bar, whose length is l, subjected to axial load P, bending moment M, and shear force Q at its two ends. The equilibrium differential governing equation [30] is(1)EId2wxdx2=Pwx-M1-Q1x,where E is elastic modulus; I is the moment of inertia; x denotes the coordinates; w(x) represents the corresponding displacement vector; P is axial load acting on the bar and is taken as negative when it is pressure; M1 is bending moment acting on end 1 of the compressive bar; and Q1 is shear force, as shown in Figure 1.Figure 1
Force diagram map of elastic compressive bar.Based on second derivative, (1) can be changed into(2)EId4wxdx4+λ2d2wxdx2=0,in which λ=P/EI.The solution of (2) is(3)wx=c1+c2x+c3sin⁡λx+c4cos⁡⁡λx=∑i=14fici.


## 3.2. Solution of Displacement Parameters of Element

Equation (3) can also be expressed as(4)wx=∑i=14fici=FC,in which(5)C=C1C2C3C4T,F=f1xf2xf3xf4x=1xsin⁡λxcos⁡⁡λx.By using differential calculation, we easily obtain the expression of the angle caused by the displacement of the compressive bar(6)θx=-dwxdx=-c2-c3λcos⁡⁡λx+c4λsin⁡λx.Define the vector equation of the nodal displacement and corner displacement as(7)δe=w1θ1w2θ2T.Combining (3) with (4) yields the expression of nodal displacement of compressive bar(8)δe=f10f20f30f40-f1′0-f2′0-f3′0-f4′0f1lf2lf3lf4l-f1′l-f2′l-f3′l-f4′lC=AC,in which(9)A=10010-1-λ01lsin⁡λlcos⁡⁡λl0-1-λcos⁡⁡λlλsin⁡λl.Based on (8), the expression of displacement function coefficient is (10)C=A-1δe.The inversion of matrixA is(11)A-1=a11a12a13a14a21a22a23a24a31a32a33a34a41a42a43a44.Each parameter of (11) can be expressed as(12)a11=λlsin⁡λl+cos⁡⁡λl-1R,a12=λlcos⁡⁡λl-sin⁡λlλR,a13=cos⁡⁡λl-1R,a14=--sin⁡λl+λlλR,a21=-λsin⁡λlR,a22=-cos⁡⁡λl-1R,a23=λsin⁡λlR,a24=-cos⁡⁡λl-1R,a31=sin⁡λlR,a32=-λlsin⁡λl+cos⁡⁡λl-1λR,a33=-sin⁡λlR,a34=cos⁡⁡λl-1λR,a41=cos⁡⁡λl-1R,a42=-λlcos⁡⁡λl-sin⁡λlλR,a43=-cos⁡⁡λl-1R,a44=-sin⁡λl+λlλR,where R=-2+λlsinλl+2cos⁡λl.


## 3.3. Analytical Trial Function

Substitution of (10) into (4) yields the element displacement expression of compressive bar(13)wx=FA-1δe=Nδe,in which(14)N=FA-1=N1N2N3N4,where N1, N2, N3, and N4 are shape functions for each displacement parameter, which can be derived by matrix calculation:(15)N1=λl-λx+sin⁡λxsin⁡λl+cos⁡⁡λx+1cos⁡⁡λl-1R,N2=--λlcos⁡⁡λl+λlcos⁡⁡λx-λl+sin⁡λx-sin⁡λlcos⁡⁡λl-1+λxλR,N3=λx-sin⁡λxsin⁡λl+cos⁡⁡λx-11-cos⁡⁡λlR,N4=-λx-sin⁡λxcos⁡⁡λl-1+cos⁡⁡λx-1sin⁡λl-λlλR.


## 4. Potential Energy and Stiffness Matrix of Compressive Bar

The finite element formula of compressive bar element is developed according to the potential energy principle and analytical shape function. The deformation of compressive bar subjected to axial load is geometrical nonlinear; the elastic and geometrical stiffness matrices are got by the variation of total potential energy.


### 4.1. The Calculation of Total Potential Energy of Element

By applying the potential energy principle, the function of total potential energy of system is(16)∏p=Ue+Up=12∫0lEIw′′2xdx-12P∫0lw′2xdx-∑i=12Qxiwxi-∑i=12Mxiw′xi,where Ue is strain energy of element; Up is potential energy of element; Q(xi) and M(xi) are external force and bending moment acting on bar, respectively.Substitution of (13) into (16) yields(17)∏p=12δeT∫0lEIN′′TN′′dxδe-12PδeT∫0lN′TN′dxδe-∑i=12Qiwi-∑i=12Miθi,where Qi and Mi are the shear force and bending moment.Based on the principle of minimum potential energy, (17) should meet the following condition:(18)∂∏p∂δe=0.According to (18), the equilibrium equation is derived:(19)Ke-PGeδe=V,where V=Q1M1Q2M2T; Ke and Ge are the elastic and geometrical stiffness matrices and the terms of these matrices are(20)Ke=∫0lEIN′′TN′′dx,(21)Ge=∫0lN′TN′dx.


### 4.2. Elastic Stiffness Matrix of Element

The equation of elastic stiffness matrix can be derived from (20)(22)K=EIk11k12k13k14k12k22k23k24k13k23k33k34k14k24k34k44,(23)kij=∫0lNi′′Nj′′dx.By using the integral calculation, we can easily obtain each of the stiffness parameters from (23); for example,(24)k11=λ3-sin⁡λl+λlS,in which S=cos⁡⁡λl+1λ2l2-4+8-4λlsin⁡λl.


### 4.3. Geometrical Stiffness Matrix of Element

Consider(25)G=g11g12g13g14g12g22g23g24g13g23g33g34g14g24g34g44,(26)gij=∫0lNi′Nj′dx.By using the integral calculation, we can easily obtain the geometrical stiffness parameters from (26); for example, (27)g11=λλlcos⁡⁡λl+2-3sin⁡λlS,where S has been given above.


## 4.1. The Calculation of Total Potential Energy of Element

By applying the potential energy principle, the function of total potential energy of system is(16)∏p=Ue+Up=12∫0lEIw′′2xdx-12P∫0lw′2xdx-∑i=12Qxiwxi-∑i=12Mxiw′xi,where Ue is strain energy of element; Up is potential energy of element; Q(xi) and M(xi) are external force and bending moment acting on bar, respectively.Substitution of (13) into (16) yields(17)∏p=12δeT∫0lEIN′′TN′′dxδe-12PδeT∫0lN′TN′dxδe-∑i=12Qiwi-∑i=12Miθi,where Qi and Mi are the shear force and bending moment.Based on the principle of minimum potential energy, (17) should meet the following condition:(18)∂∏p∂δe=0.According to (18), the equilibrium equation is derived:(19)Ke-PGeδe=V,where V=Q1M1Q2M2T; Ke and Ge are the elastic and geometrical stiffness matrices and the terms of these matrices are(20)Ke=∫0lEIN′′TN′′dx,(21)Ge=∫0lN′TN′dx.


## 4.2. Elastic Stiffness Matrix of Element

The equation of elastic stiffness matrix can be derived from (20)(22)K=EIk11k12k13k14k12k22k23k24k13k23k33k34k14k24k34k44,(23)kij=∫0lNi′′Nj′′dx.By using the integral calculation, we can easily obtain each of the stiffness parameters from (23); for example,(24)k11=λ3-sin⁡λl+λlS,in which S=cos⁡⁡λl+1λ2l2-4+8-4λlsin⁡λl.


## 4.3. Geometrical Stiffness Matrix of Element

Consider(25)G=g11g12g13g14g12g22g23g24g13g23g33g34g14g24g34g44,(26)gij=∫0lNi′Nj′dx.By using the integral calculation, we can easily obtain the geometrical stiffness parameters from (26); for example, (27)g11=λλlcos⁡⁡λl+2-3sin⁡λlS,where S has been given above.


## 5. Examples and Comparison

To verify the accuracy and validity of the analytical trial function proposed in this paper, a number of numerical examples are presented in this section. The results obtained from the analytical shape function are compared with those of interpolation shape function element and those of the direct stiffness method.


### 5.1. Calculation Model

In this section, four compressive bars with different constraints are analyzed to calculate their bifurcation buckling load. The dimensions of all those compressive bar are the same and are0.2m×0.3m, the lengths are all 3 m, the elastic modulus of material is E=210 GPa, and the Poisson’s ratio is υ=0.3.The types of those bars are Cantilever bar, simply supported bar, bar with one end pinned and the other end sliding supported, and bar with one end pinned and the other simply supported.


### 5.2. Theoretical Solution of Bifurcation Buckling Load

The solution of bifurcation buckling load of compressive bar with different constraints is calculated based on Euler’s formula [30]:(28)Pcr=π2EIμl2,where Pcr is bifurcation buckling load; μ is effective length coefficient of bar, whose value can be got from [29].


### 5.3. Numerical Solution of Bifurcation Buckling Load

#### 5.3.1. Step of Solving the Numerical Solution

(1)
The characteristic equation for stability problem of elastic compressive bar can be expressed asKδ+βP[G]δ=0. By transformation, the characteristic equation of stability problem is (1/P)[G]-1Kδ+βδ=0.
(2)
Letβ be the iteration variable, and assuming the initial value λ-0=1, Pcr0=P0=100N, and by using the iterative method, the characteristic value of the equation is obtained. Adopt the minimum characteristic value βi+1 as the final solution.
(3)
The value of characteristic loadsPi+1=βi+1Pi and λ-i=Pi+1/EI can be determined from the value of βi+1. When the relative error reaches β-1≤10-6, the iteration stops. Otherwise, repeat step (2). Then the bifurcation buckling load is developed. In which, β=Pcr/P, δ is the corresponding value of buckling mode.


#### 5.3.2. Direct Stiffness Method

By solving the fourth-order stiffness differential equations, the general solution of deformationwx is obtained. Let the nodal deformation δe be the unknown variable, and then the nodal displacement equations are established which contain the undetermined coefficients. By solving the undetermined coefficients and considering the axial and lateral deformation, the stiffness matrix equation [31] for stability problem of element is obtained. By calculating the equilibrium equation, the buckling load is obtained.


#### 5.3.3. Finite Element Method

(1) Interpolation Trial Function Element. In this method, based on the characteristic of deformation of compressive bar, the deformation trial function is created, in which the deformation trial functions generated by axial load and by bending moment are calculated by the Lagrange interpolation function and Hermite interpolation function, respectively. According to the above interpolation function, the potential energy function is established. Based on the principle of minimum potential energy and considering the large deformation of the compressive bar, by the variation of the potential energy, the elastic and geometrical stiffness matrices [32] are got and then the buckling load is obtained.(2) Analytical Trial Function Element. The process and method are the same as above.


### 5.4. Results and Discussion

#### 5.4.1. Element Precision

In this section, to calculate the buckling load of the elastic compressive bar under various working conditions, the analytical shape function method, the traditional finite element method, and the direct stiffness method without dividing element are all used. And those results and their accuracy are compared with the analytical solution derived from the Euler’s formula, as shown in Table1.Table 1
Comparison of calculated buckling loads of bar from different methods.
Types of bar
Theoretical   solution
Direct stiffnessmethod
Interpolation shapefunction method
Proposedmethod
BucklingLoad  (N)
RE (%)
Bucklingload (N)
RE (%)
Bucklingload (N)
RE (%)
Cantilever bar
1151.45
1151.41
4.13
E
-
3
1160.10
7.51
E
-
1
1151.45
1.80
E
-
8
Simply supported bar
4605.82
4605.62
4.13
E
-
3
5600.00
21.6
4605.82
3.80
E
-
9
Bar with one end pinned and the other end sliding supported
4605.82
4605.62
4.13
E
-
3
4666.70
1.32
4605.82
3.80
E
-
9
Bar with one end pinned  and the other  simply supported
9422.34
9421.88
4.94
E
-
3
14000.00
48.6
9422.34
6.37
E
-
8
Note: the sliding supported means that the bar can slide in longitudinal operation.The comparison shows that there is no need to divide the bar into multiple elements but only one element to get the accurate buckling load for the proposed finite element method, and only three iterations are needed to get this exact solution, which means that the iterative convergence speed is fast and the efficiency is high. The relative error is 10−10, which has a little difference with the real solution. This proves that this method has the analytical and accuracy result. Comparison between the theoretical analytical solution and the traditional finite element solution shows that the relative errors are 0.75% for cantilever bar and 1.3% for bar with one end pinned and the other end sliding supported, but for other restrained bars, the errors are large, which state that the precision between the polynomial displacement shape function and the actual displacement shape function is poor; the accuracy of the direct stiffness method is higher than that of the traditional finite element method but lower compared with the method proposed in this paper.


#### 5.4.2. Element Efficiency

For the compressive bar with one end pinned and the other end simply supported, in the traditional finite element method, by dividing the compressive bar into different number of elements, different buckling load values are got. And then those values are compared with the result calculated from the analytical trial function. The comparisons are shown in Table2.Table 2
Buckling load of bar with different number of elements (N).
Interpolation shape function element
Analytical trial function
Number of elements
1
2
3
4
5
8
20
Single element
Theory value
9422.3
—
Buckling load
14000
9664.1
9480.2
9441.7
9430.5
9422.4
9422.4
9422.34
Relative error (%)
48.58
2.56
6.14
E
-
1
2.05
E
-
1
8.66
E
-
2
1.34
E
-
2
6.38
E
-
04
6.37
E
-
8As shown in Table2, the accuracy is poor for traditional finite element without meshing; only when the compressive bar is divided into more than 20 elements, the calculation precision is close to theoretical solution. But for the method proposed in this paper, there is no need to divide the bar into elements and higher precision result for buckling load is got. That means that analytical trial function method has higher efficiency than the interpolation shape function element to create the compressive bar element.


### 5.5. Analysis of the Solutions

Based on above calculations and comparisons, the following results can be obtained:(1)
Due to the precision of displacement trial function method meeting the actual result, there is no need to mesh the bar when analyzing the stability of the bar, for the element created by analytical trial function is accurate.(2)
By applying analytical trial function to get an accuracy solution, the element does not need mesh; but for the interpolation shape function, to get a more accurate solution, the bar needs to be divided into more elements. Thus, the analytical trial function has higher efficiency than other functions to create compressive bar element.(3)
Infinite degrees of freedom of rigid-body are changed into limited degrees of freedom through the finite element method, which means that the constraints in the structure increased, resulting in the fact that the structure stiffness is increased and the buckling load is bigger than the actual value, as shown in Table2.(4)
For direct stiffness method, meshing has a direct effect on calculation precision when creating the compression bar element. With only one element, the results could not agree with those derived from the analytical solution. Thus, the precision of the direct stiffness method is lower than the analytical method.


## 5.1. Calculation Model

In this section, four compressive bars with different constraints are analyzed to calculate their bifurcation buckling load. The dimensions of all those compressive bar are the same and are0.2m×0.3m, the lengths are all 3 m, the elastic modulus of material is E=210 GPa, and the Poisson’s ratio is υ=0.3.The types of those bars are Cantilever bar, simply supported bar, bar with one end pinned and the other end sliding supported, and bar with one end pinned and the other simply supported.


## 5.2. Theoretical Solution of Bifurcation Buckling Load

The solution of bifurcation buckling load of compressive bar with different constraints is calculated based on Euler’s formula [30]:(28)Pcr=π2EIμl2,where Pcr is bifurcation buckling load; μ is effective length coefficient of bar, whose value can be got from [29].


## 5.3. Numerical Solution of Bifurcation Buckling Load

### 5.3.1. Step of Solving the Numerical Solution

(1)
The characteristic equation for stability problem of elastic compressive bar can be expressed asKδ+βP[G]δ=0. By transformation, the characteristic equation of stability problem is (1/P)[G]-1Kδ+βδ=0.
(2)
Letβ be the iteration variable, and assuming the initial value λ-0=1, Pcr0=P0=100N, and by using the iterative method, the characteristic value of the equation is obtained. Adopt the minimum characteristic value βi+1 as the final solution.
(3)
The value of characteristic loadsPi+1=βi+1Pi and λ-i=Pi+1/EI can be determined from the value of βi+1. When the relative error reaches β-1≤10-6, the iteration stops. Otherwise, repeat step (2). Then the bifurcation buckling load is developed. In which, β=Pcr/P, δ is the corresponding value of buckling mode.


### 5.3.2. Direct Stiffness Method

By solving the fourth-order stiffness differential equations, the general solution of deformationwx is obtained. Let the nodal deformation δe be the unknown variable, and then the nodal displacement equations are established which contain the undetermined coefficients. By solving the undetermined coefficients and considering the axial and lateral deformation, the stiffness matrix equation [31] for stability problem of element is obtained. By calculating the equilibrium equation, the buckling load is obtained.


### 5.3.3. Finite Element Method

(1) Interpolation Trial Function Element. In this method, based on the characteristic of deformation of compressive bar, the deformation trial function is created, in which the deformation trial functions generated by axial load and by bending moment are calculated by the Lagrange interpolation function and Hermite interpolation function, respectively. According to the above interpolation function, the potential energy function is established. Based on the principle of minimum potential energy and considering the large deformation of the compressive bar, by the variation of the potential energy, the elastic and geometrical stiffness matrices [32] are got and then the buckling load is obtained.(2) Analytical Trial Function Element. The process and method are the same as above.


## 5.3.1. Step of Solving the Numerical Solution

(1)
The characteristic equation for stability problem of elastic compressive bar can be expressed asKδ+βP[G]δ=0. By transformation, the characteristic equation of stability problem is (1/P)[G]-1Kδ+βδ=0.
(2)
Letβ be the iteration variable, and assuming the initial value λ-0=1, Pcr0=P0=100N, and by using the iterative method, the characteristic value of the equation is obtained. Adopt the minimum characteristic value βi+1 as the final solution.
(3)
The value of characteristic loadsPi+1=βi+1Pi and λ-i=Pi+1/EI can be determined from the value of βi+1. When the relative error reaches β-1≤10-6, the iteration stops. Otherwise, repeat step (2). Then the bifurcation buckling load is developed. In which, β=Pcr/P, δ is the corresponding value of buckling mode.


## 5.3.2. Direct Stiffness Method

By solving the fourth-order stiffness differential equations, the general solution of deformationwx is obtained. Let the nodal deformation δe be the unknown variable, and then the nodal displacement equations are established which contain the undetermined coefficients. By solving the undetermined coefficients and considering the axial and lateral deformation, the stiffness matrix equation [31] for stability problem of element is obtained. By calculating the equilibrium equation, the buckling load is obtained.


## 5.3.3. Finite Element Method

(1) Interpolation Trial Function Element. In this method, based on the characteristic of deformation of compressive bar, the deformation trial function is created, in which the deformation trial functions generated by axial load and by bending moment are calculated by the Lagrange interpolation function and Hermite interpolation function, respectively. According to the above interpolation function, the potential energy function is established. Based on the principle of minimum potential energy and considering the large deformation of the compressive bar, by the variation of the potential energy, the elastic and geometrical stiffness matrices [32] are got and then the buckling load is obtained.(2) Analytical Trial Function Element. The process and method are the same as above.


## 5.4. Results and Discussion

### 5.4.1. Element Precision

In this section, to calculate the buckling load of the elastic compressive bar under various working conditions, the analytical shape function method, the traditional finite element method, and the direct stiffness method without dividing element are all used. And those results and their accuracy are compared with the analytical solution derived from the Euler’s formula, as shown in Table1.Table 1
Comparison of calculated buckling loads of bar from different methods.
Types of bar
Theoretical   solution
Direct stiffnessmethod
Interpolation shapefunction method
Proposedmethod
BucklingLoad  (N)
RE (%)
Bucklingload (N)
RE (%)
Bucklingload (N)
RE (%)
Cantilever bar
1151.45
1151.41
4.13
E
-
3
1160.10
7.51
E
-
1
1151.45
1.80
E
-
8
Simply supported bar
4605.82
4605.62
4.13
E
-
3
5600.00
21.6
4605.82
3.80
E
-
9
Bar with one end pinned and the other end sliding supported
4605.82
4605.62
4.13
E
-
3
4666.70
1.32
4605.82
3.80
E
-
9
Bar with one end pinned  and the other  simply supported
9422.34
9421.88
4.94
E
-
3
14000.00
48.6
9422.34
6.37
E
-
8
Note: the sliding supported means that the bar can slide in longitudinal operation.The comparison shows that there is no need to divide the bar into multiple elements but only one element to get the accurate buckling load for the proposed finite element method, and only three iterations are needed to get this exact solution, which means that the iterative convergence speed is fast and the efficiency is high. The relative error is 10−10, which has a little difference with the real solution. This proves that this method has the analytical and accuracy result. Comparison between the theoretical analytical solution and the traditional finite element solution shows that the relative errors are 0.75% for cantilever bar and 1.3% for bar with one end pinned and the other end sliding supported, but for other restrained bars, the errors are large, which state that the precision between the polynomial displacement shape function and the actual displacement shape function is poor; the accuracy of the direct stiffness method is higher than that of the traditional finite element method but lower compared with the method proposed in this paper.


### 5.4.2. Element Efficiency

For the compressive bar with one end pinned and the other end simply supported, in the traditional finite element method, by dividing the compressive bar into different number of elements, different buckling load values are got. And then those values are compared with the result calculated from the analytical trial function. The comparisons are shown in Table2.Table 2
Buckling load of bar with different number of elements (N).
Interpolation shape function element
Analytical trial function
Number of elements
1
2
3
4
5
8
20
Single element
Theory value
9422.3
—
Buckling load
14000
9664.1
9480.2
9441.7
9430.5
9422.4
9422.4
9422.34
Relative error (%)
48.58
2.56
6.14
E
-
1
2.05
E
-
1
8.66
E
-
2
1.34
E
-
2
6.38
E
-
04
6.37
E
-
8As shown in Table2, the accuracy is poor for traditional finite element without meshing; only when the compressive bar is divided into more than 20 elements, the calculation precision is close to theoretical solution. But for the method proposed in this paper, there is no need to divide the bar into elements and higher precision result for buckling load is got. That means that analytical trial function method has higher efficiency than the interpolation shape function element to create the compressive bar element.


## 5.4.1. Element Precision

In this section, to calculate the buckling load of the elastic compressive bar under various working conditions, the analytical shape function method, the traditional finite element method, and the direct stiffness method without dividing element are all used. And those results and their accuracy are compared with the analytical solution derived from the Euler’s formula, as shown in Table1.Table 1
Comparison of calculated buckling loads of bar from different methods.
Types of bar
Theoretical   solution
Direct stiffnessmethod
Interpolation shapefunction method
Proposedmethod
BucklingLoad  (N)
RE (%)
Bucklingload (N)
RE (%)
Bucklingload (N)
RE (%)
Cantilever bar
1151.45
1151.41
4.13
E
-
3
1160.10
7.51
E
-
1
1151.45
1.80
E
-
8
Simply supported bar
4605.82
4605.62
4.13
E
-
3
5600.00
21.6
4605.82
3.80
E
-
9
Bar with one end pinned and the other end sliding supported
4605.82
4605.62
4.13
E
-
3
4666.70
1.32
4605.82
3.80
E
-
9
Bar with one end pinned  and the other  simply supported
9422.34
9421.88
4.94
E
-
3
14000.00
48.6
9422.34
6.37
E
-
8
Note: the sliding supported means that the bar can slide in longitudinal operation.The comparison shows that there is no need to divide the bar into multiple elements but only one element to get the accurate buckling load for the proposed finite element method, and only three iterations are needed to get this exact solution, which means that the iterative convergence speed is fast and the efficiency is high. The relative error is 10−10, which has a little difference with the real solution. This proves that this method has the analytical and accuracy result. Comparison between the theoretical analytical solution and the traditional finite element solution shows that the relative errors are 0.75% for cantilever bar and 1.3% for bar with one end pinned and the other end sliding supported, but for other restrained bars, the errors are large, which state that the precision between the polynomial displacement shape function and the actual displacement shape function is poor; the accuracy of the direct stiffness method is higher than that of the traditional finite element method but lower compared with the method proposed in this paper.


## 5.4.2. Element Efficiency

For the compressive bar with one end pinned and the other end simply supported, in the traditional finite element method, by dividing the compressive bar into different number of elements, different buckling load values are got. And then those values are compared with the result calculated from the analytical trial function. The comparisons are shown in Table2.Table 2
Buckling load of bar with different number of elements (N).
Interpolation shape function element
Analytical trial function
Number of elements
1
2
3
4
5
8
20
Single element
Theory value
9422.3
—
Buckling load
14000
9664.1
9480.2
9441.7
9430.5
9422.4
9422.4
9422.34
Relative error (%)
48.58
2.56
6.14
E
-
1
2.05
E
-
1
8.66
E
-
2
1.34
E
-
2
6.38
E
-
04
6.37
E
-
8As shown in Table2, the accuracy is poor for traditional finite element without meshing; only when the compressive bar is divided into more than 20 elements, the calculation precision is close to theoretical solution. But for the method proposed in this paper, there is no need to divide the bar into elements and higher precision result for buckling load is got. That means that analytical trial function method has higher efficiency than the interpolation shape function element to create the compressive bar element.


## 5.5. Analysis of the Solutions

Based on above calculations and comparisons, the following results can be obtained:(1)
Due to the precision of displacement trial function method meeting the actual result, there is no need to mesh the bar when analyzing the stability of the bar, for the element created by analytical trial function is accurate.(2)
By applying analytical trial function to get an accuracy solution, the element does not need mesh; but for the interpolation shape function, to get a more accurate solution, the bar needs to be divided into more elements. Thus, the analytical trial function has higher efficiency than other functions to create compressive bar element.(3)
Infinite degrees of freedom of rigid-body are changed into limited degrees of freedom through the finite element method, which means that the constraints in the structure increased, resulting in the fact that the structure stiffness is increased and the buckling load is bigger than the actual value, as shown in Table2.(4)
For direct stiffness method, meshing has a direct effect on calculation precision when creating the compression bar element. With only one element, the results could not agree with those derived from the analytical solution. Thus, the precision of the direct stiffness method is lower than the analytical method.


## 6. Conclusions

Based on the results of the investigation, the following conclusions can be obtained:(1)
The result of interpolation shape function element for calculating the buckling load is not in accordance with theoretical solution, and the relative error is larger. Only with element number increasing, the result is close to theoretical solution. It indicates that the precision requirement cannot be met by taking the interpolation shape function element.(2)
The solution of direct stiffness method for buckling load is close to theoretical solution. But in this method, meshing still has a direct effect on precision when creating the bar element, which reduces the efficiency.(3)
Elastic and geometrical stiffness matrices are obtained by applying analytical trial function in this paper. Buckling loads of four typical compressive bars are attained by analyzing eigenvalue and eigenvector and are compared with theoretical solution. The results show that the series of relative error is about1×10-9. The stiffness matrices in the light of analytical trial function proposed in this paper are high-precision element.



---

*Source: 762539-2015-05-19.xml*