ABSTRACT

The quasi-static explicit finite element method (FEM) and element free Galerkin (EFG) method are applied to trace the post-buckling equilibrium path of thin-walled members in this paper. The factors that primarily control the explicit buckling solutions, such as the computation time, loading function and dynamic relaxation, are investigated and suggested for the buckling analysis of thin-walled members. Three examples of different buckling modes, namely snap-through, overall and local buckling, are studied based on the implicit FEM, quasi-static explicit FEM and EFG method via the commercial software LS-DYNA. The convergence rate and accuracy of the explicit methods are compared with the conventional implicit arc-length method. It is drawn that EFG quasi-static explicit buckling analysis presents the same accurate results as implicit finite element solution, but is without convergence problem and of less-consumption of computing time than FEM.

Keywords:

Quasi-Static Explicit Solution, Buckling Analysis, Thin-Walled Members, FEM, EFG

1. Introduction

Thin-walled members of various shapes have been widely used in civil and mechanical engineering. Under many conditions, when these members are subjected to compressive loads, buckling and post-buckling behaviors usually dictate the design considerations. Analytical solutions of buckling of thin-walled members can only obtain for the simple cases of elastic linear/nonlinear buckling. Analysis of nonlinear buckling problems of plastic and large deformations primarily relies on the numerical methodology.

In the nonlinear buckling analysis of thin-walled structural member, the implicit arc-length algorithm is generally accepted as an effective numerical method for tracing the post-buckling path [1] [2] [3] . It is well known that the implicit method is established on the operation of stiffness matrix, where a large amount of computing resource is required for processing highly incremental iteration. For complex nonlinear problems, the disadvantage of non-conver- gence of this method is obvious and usually difficult to be overcome, for example, the singularity of the stiffness matrix near the critical point. Therefore it is necessary to study the explicit method for solving highly nonlinear post-buck- ling problems.

Quasi-static analysis is a simulation of static problem with motion analysis which restricts the load velocity so that the outcome of this analysis can only have a little inertia influence that can be neglected. As an explicit algorithm, the advantages of quasi-static buckling analysis lie in the lower computing cost and no convergence consideration. However, structural dynamic responses caused by loading speed and inertia force significantly influence the quasi-static explicit results. Usually very small loading step is needed to approach the static equilibrium state at each loading moment, which inversely decreases the convergence rate. The efficiency of quasi-static explicit method mainly depends on the problems being solved. The key factors that can reduce the dynamic responses, like computing time, loading function and damping relaxation, must be specified in order to keeping the ratio of dynamic energy to internal energy within a low level. Zhuang [4] presented the conventional method of quasi-static analysis and compared the differences between implicit and explicit methods. Ji [5] used the quasi-static analysis for solving stable problems of stiffened plates under axial pressure, and obtained the structure responses which static analysis could not reach with a little time cost. Fan [6] did research on the effect of the duration and functional form of the time histories of loads by the analysis of a linear spring mass oscillator. Li [7] interpreted the principle of static calculation by using dynamic relaxation method and put forward a new method on value adoption of mass damp and rigidity damp. Lee [8] applied the explicit arc-length method to trace the post-buckling equilibrium path of structures on the basis of dynamic relaxation method with kinetic damping.

Finite element method (FEM) is a stable and reliable computing method through meshing the continuum into discrete units. When structures undergo large deformations, the computing accuracy is significantly influenced by the distortion of discrete units. In explicit method, a stable time step must be very small if the distortion of discrete units occurs, which greatly adds up the computing cost. Element-free method (EFM) is studied by many researchers for avoiding the effects of discrete units on numerical consequence. Solution with EFM depends on the discrete nodes setting up within or at the edge of a domain. Shape function is constructed on local nodes, so there is no mesh-dependence problem. The primary advantage of this approach is that there is no singularity of stiffness matrix induced by distortion of discrete units in the solution of large deformation and discontinuity problems. Element free Galerkin (EFG) method is based on the global Galerkin weak forms and the integration of background grids. The displacement approximation functions are generated by using the least squares approximation constructed via nodes in local fields. The computational accuracy and convergence rate of EFG methods have been demonstrated to be the same as FEM. The stability of this method is not affected by the irregular nodes, and furthermore, it can be combined with FEM and BEM (boundary element method) to improve the computing efficiency.

EFG method has been well used in the buckling analysis of thin-walled members. Liu [9] developed an EFG formulation to calculate the buckling loads of symmetrically laminated composite plates based on the principle of minimum potential energy, and found that solving the eigenvalue problem is much more computationally efficient compared to the FEM. Chinnaboon [10] developed a BEM-based meshless method for buckling analysis of elastic plates with various boundary conditions that include elastic supports and restraints. Liew [11] used an EFG method to study the elastic buckling behavior of stiffened and un-stif- fened folded plates under partial in-plane edge loads. Tamijani [12] employed the EFG method for buckling and static analysis of plates with arbitrary curvilinear stiffeners. Peng [13] obtained the critical buckling load of ribbed plates using the mesh-free method based on the first-order shear deformation theory. Xiang [14] predicted buckling behavior of microtubules based on an atomistic- continuum model. Lu [15] developed an adaptive enrichment mesh-free method to capture wrinkling and post-buckling behavior in sheet metal forming. Li [16] used mesh-free method for numerical simulations of large deformation of thin shell structures, which showed simplicity in both formulation and implementation as compared to shell theory approach. Lin [17] used a non-linear dynamic explicit scheme for the post-buckling analysis of thin-walled structure based on the meshless shell formulation. Compared with the finite element method, the mesh-free method possesses the same accuracy and can save some computing time as well as work out the problems that can’t be solved by the traditional FEM.

The FEM and EFG quasi-static explicit methodologies are applied to trace the post-buckling path of thin-walled members in this paper. The key factors that control the convergence rate and dynamic responses, such as the computation time, loading function and damping relaxation, are discussed and suggested in the numerical buckling analysis. Three examples of thin-walled members occurred snap-through, overall and locally buckling are studied in detail by quasi-static explicit FEM and EFG method, and the efficiency and accuracy of the applied methods are demonstrated through the comparison with the conventional solution of implicit arc-length method.

2. Formulation of Quasi-Static Explicit FEM and EFG Algorithms

In explicit formulation, the basic dynamic equation of a volume element at time t is written in the form of equation of equilibrium as follow:

${\sigma }_{ij,j}+f_i-\rho u_{i,tt}-\mu u_{i,t}=0\text{in}V$ (1)

The constitutive equation and the boundary equations are given by

${\sigma }_{ij}=D_{ijkl}{\epsilon }_{kl}\text{in}V$ (2)

$u_i={\overline{u}}_i\text{on}{S}_u$ (3)

$s_{ij}{n}_j={\overline{T}}_i\text{on}{S}_{\sigma }$ (4)

where σ is the Gauss stress, f_i is body force, ρ is density, μ is damping coefficient, u_i,t and u_i,tt are the first and second derivatives of displacement u_i, D_ijkl is stiffness tensor, ε_kl is strain tensor, n_j is normal vector, and T_i is surface force.

By applying the Galerkin method to the Equation (1), Equation (3) and Equation (4), the corresponding weak form is expressed as

${\displaystyle \underset{V}{\int }\delta u_i\left({\sigma }_{ij,j}+f_i-\rho u_{i,tt}-\mu u_{i,t}\right)\text{d}V}-{\displaystyle \underset{S_{\sigma }}{\int }\delta u_i\left({\sigma }_{ij}{n}_j-{\overline{T}}_i\right)\text{d}s}=0$ (5)

Substituting Equation (2) to Equation (5), the weak form is then transformed into Equation (6).

${\displaystyle \underset{V}{\int }\left(\delta {\epsilon }_{ij}{D}_{ijkl}{\epsilon }_{kl}+\delta u_i\rho u_{i,tt}+\delta u_i\mu u_{i,t}\right)\text{d}V}={\displaystyle \underset{V}{\int }\delta u_i{f}_i}\text{d}V+{\displaystyle \underset{S_{\sigma }}{\int }\delta u_i{\overline{T}}_i}\text{d}s$ (6)

${\displaystyle \underset{V}{\int }\left(\delta {\epsilon }^{\text{T}}D\epsilon +\delta u^{\text{T}}\rho u+\delta u^{\text{T}}\mu u\right)\text{d}V}={\displaystyle \underset{V}{\int }\delta u^{\text{T}}f}\text{d}V+{\displaystyle \underset{S_{\sigma }}{\int }\delta u^{\text{T}}T}\text{d}s$ (7)

In the FEM, the shape function, N(x), is created by interpolation in elements of a set of fixed nodes. The displacement function is given by

$u\left(x\right)={\displaystyle \sum _{i=1}^n{N}_i\left(x\right)u_i}=N\left(x\right)u$ (8)

where $N=\left[N_1^{\text{T}},N_2^{\text{T}},\cdots ,N_n^{\text{T}}\right]$ , $u=\left[u_1^{\text{T}},u_2^{\text{T}},\cdots ,u_n^{\text{T}}\right]$ . Substituting Equation (8) into Equation (7) leads to the equation as below:

#Math_12# (9)

The above dynamic equilibrium equation can be reduced to the general form as follow.

$M\ddot{a}\left(t\right)+C\dot{a}\left(t\right)+Ka\left(t\right)=Q\left(t\right)$ (10)

in which

$M={\displaystyle \sum _e{M}^e}={\displaystyle \sum _e{\displaystyle {\int }_{V_e}\rho N^{\text{T}}N\text{d}V}},C={\displaystyle \sum _e{C}^e={\displaystyle \sum _e{\displaystyle {\int }_{V_e}\mu N^{\text{T}}N\text{d}V}}}$

$K={\displaystyle \sum _e{K}^e={\displaystyle \sum _e{\displaystyle {\int }_{V_e}{B}^{\text{T}}DB\text{d}V}}},Q={\displaystyle \sum _e{Q}^e={\displaystyle \sum _e\left({\displaystyle {\int }_{V_e}{N}^{\text{T}}f\text{d}V}+{\displaystyle {\int }_{S_{\sigma }^e}{N}^{\text{T}}T\text{d}S}\right)}}$

Different from FEM, the numerical discretization in EFG method is based on the moving least-squares (MLS) approximation [18] . The displacement function, u(x), is approximated by using the polynomial basis as follows:

$u\left(x\right)={\displaystyle \sum _{i=1}^m{p}_i\left(x\right)a_i\left(x\right)}=p^{\text{T}}\left(x\right)a\left(x\right)$ (11)

where m is the order of completeness in this approximation, the monomial p_i(x) is the basis function, and a_i(x) is the coefficient of the approximation. a_i(x) depends on the sampling point x_i that is collected by a weighting function $w\left(x\right)=w\left(x-x_i\right)$ , which is nonzero in a small domain called influence domain.

The weighted residual can be written as L₂-Norm, namely

$J={\displaystyle \sum _{I=1}^N{w}_I\left(x\right){\left[u\left(x_I\right)-u_I\right]}^2}={\displaystyle \sum _{I=1}^N{w}_I\left(x\right){\left[{\displaystyle \sum _{i=1}^m{p}_i}\left(x_I\right)a_i\left(x\right)-u_I\right]}^2}$ (12)

In MLS approximation, at an arbitrary point x, a(x) is chosen on the basis of minimizing the weighted residual, then we have:

$a\left(x\right)=A^{-1}\left(x\right)B\left(x\right)u$ (13)

Substitute Equation (13) to Equation (11), it gives

#Math_20# (14)

in which:

$A\left(x\right)={\displaystyle \sum _{I=1}^N{w}_I\left(x\right)}p\left(x_I\right)p^{\text{T}}\left(x_I\right)$

$B\left(x\right)=\left[\begin{array}{cccc}{w}_1\left(x\right)p\left(x_1\right)& w_2\left(x\right)p\left(x_2\right)& \cdots & w_N\left(x\right)p\left(x_N\right)\end{array}\right]$

By applying Equation (14) to Equation (10), the quasi-static explicit EFG formulation can be obtained.

3. Factors in Quasi-Static Explicit Buckling Analysis

In quasi-static explicit buckling analysis, to achieve the real dynamic loading process and the progress of unbalanced stress waves between elements, a large amount of time increment steps are usually required to obtain stable solutions. The kinetic energy usually increases rapidly with the enlargement of deformation after peak load point. For improving computing efficiency and having a stable solution, measures have to be employed to accelerate the computing process, and at the same time, the dynamic responses expressed by the ratio of kinetic energy to internal energy must be kept within 5% - 10%. The key factors in this slow-dynamic technique are how to choose the load duration and minimize the undesired dynamic effect originated from the inertia force of the governing equation. The response of structures is mainly controlled by the first mode, so the computation time in quasi-static explicit buckling analysis is usually set up more than ten times of the first mode period T.

Rapid movement can generate stress wave, which results in shock or inaccurate numerical solutions. Therefore, the curve of loading function must be smooth. The commonly used loading functions are shown as below:

Linear function curve: $F\left(t\right)=F_{0}t/T$

Parabolic function curve: $F\left(t\right)=F_0{\left(t/T\right)}^2$

Versin function curve: $F\left(t\right)=F_0\left(1-\cos \left(\text{π}t/T\right)\right)/2$

Cycloid function curve: $F\left(t\right)=F_0\left(t/T-\sin \left(2\text{π}t/T\right)/2\text{π}\right)$

Dynamic relaxation is originated in the steady-state solution of single degree of freedom damping system. The basic idea lies in keeping a system in an over damping state by setting the Rayleigh damping to a large value to weaken the dynamic effect on the system.

The damping of an actual structure can be expressed by the Rayleigh damping, namely $C=\alpha M+\beta K$ , where M is the mass matrix, K is the stiff matrix, α and β are the mass damping and stiff damping coefficients, respectively. The damping ratio ξ_i is written as

${\xi }_i=\left(\alpha /{\omega }_i+\beta {\omega }_i\right)/2$ (15)

where ω_i is the ith order circular frequency. The curves of unit step responses to different damping ratios of a single degree of freedom system are shown in Figure 1. It is found that the responses can easily approach to the static solution if ξ_i ≥ 1, preferably ξ_i = 1. As long as αβ ≥ 1, the damping ratio of the system ξ_i is more than 1 from Equation (15).

Refer to the research by Li et al. [7] , the value of mass damping and rigidity damping is taken as α = ω_min, β = 1/ω_min, where ω_min is the first order circular frequency. By this way, the shocks induced by the higher order frequencies are restrained, which can make the explicit solution of dynamic relaxation close to the static solution.

4. Application of the Numerical Methods

In this section, by using the commercial software LS-DYNA, the FEM and EFG quasi-static explicit method are employed to simulate the post-buckling behaviors of three typical buckling problems: snap-through buckling, overall buck-

ling and local buckling. The loading-displacement curves solved by explicit quasi-static FEM, implicit FEM and explicit quasi-static EFG method are compared and the reliability of explicit quasi-static buckling analysis with the suggested computation time, loading function and dynamic relaxation is demonstrated.

4.1. Example 1: Snap―Through Buckling of a Cylindrical Shell

A simply supported cylindrical shallow shell subjected to central loading is shown in Figure 2. The shell is assumed to be made of steel with uniform elastic modulus E = 210 GPa, Poisson’s ratio υ = 0.3 and density ρ = 7850 kg/m³. Length, radius, thickness and angle of the shell are L = 100 mm, R = 1000 mm, t = 4 mm and θ = 6˚, respectively. A quarter of the shell with symmetric constrain is modeled, and the displacement of 0-13mm is applied at point A. This is a typical example of snap-through buckling; in which the traditional load-controlled computing method is inapplicable due to the singularity of tangential stiffness matrix at the extreme point. The equilibrium path can be traced by implicit arc-length method, and being as a comparison, the accuracy of quasi-static explicit buckling analysis is studied from the following aspects.

As a versine displacement is applied and the computation time is specified as 10T, 20T, 40T, 80T, 160T and 320T, the curves of load-vertical displacement at point A are constructed in Figures 3-5. It is shown that the curves from explicit calculation are very close to those from implicit calculation before the critical load points. There are fluctuations about the implicit results in the rising segments of post-buckling in the explicit calculation. The shorter is the computation time, the lower is the fluctuating frequency and the higher is the fluctuating amplitude. From the ratio of kinetic energy to internal energy in Figure 6, it can be seen that the percentage decreases with the increase of computation time, but hardly decreases to 10% just by means of extending computation time.

When different loading functions such as linear, parabolic, cycloidal and versine curves are applied, and the computation time is taken as 20T, the load-dis- placement relationships are constructed in Figure 7. It is shown that the critical load of each curve is very close. The result from the applied linear displacement is the most fluctuant among all curves; and the one from parabolic displacement-loading is well coincident with that from the versine displacement-loading. The displacement from cycloidal loading is bigger than other consequences. Therefore, the versine or parabolic displacement-loading is more appropriate to be used in the quasi-static explicit buckling analysis of thin-walled members.

In order to keep the structural responses being quasi-statically damped and converged to the static equilibrium state, the dynamic relaxation (DR) method has to be used in the explicit buckling analysis. Based on the minimum frequency of free vibration of the calculated shell, 752.28 Hz determined from model analysis, coefficients α and β are calculated as α = 4726.7 and β = 2.12e⁻⁴, respectively. The load-displacement curves based on versine load and DR in the case of different computation time are plotted in Figure 8. Compared with the curves in Figures 3-5, the fluctuation on the curves almost disappears after DR is applied. The consequence of quasi-static explicit analysis is close to the implicit solution if long enough computation time is adopted. The post-bucking performance of displacement reduction obtained by implicit arc-length method cannot be identified by the explicit analysis since the monotonically increasing displacement is used in the explicit solution. However, the load-displacement curve of DR and 320T is demonstrated to be well enough to express the post-buckling path. The kinetic energy is remarkably reduced after the DR is applied as shown in Figure 9, even if the computation time equals 10T.

The computation time, loading function and dynamic relaxation are also important for quasi-static explicit post-buckling analysis with element-free Galerkin (EFG) method. The nodes in EFG are appointed based on the grid nodes in

FEM in the package LS-DYNA. Parameters and solution settings of shell are defined by the key words in *SECTION_SHELL_EFG and *CONTROL_EFG. According to the previous explicit FEM, versine load and DR are employed in the explicit EFG buckling analysis. The load-displacement curve of 160T solved by EFG is contrasted with those from the implicit arc-length method and explicit FEM as shown in Figure 10. It can be seen that the load-displacement curves from EFG and FEM explicit analysis are almost coincident, and the post-buck- ling path is well indicated as implicit solution. Furthermore, the curves from EFG are smoother than those from FEM.

4.2. Example 2: Overall Buckling of Thin-Walled Steel Angle

Dimensions of a thin-walled steel angle subjected to central loading are shown in Figure 11. The steel angle is with uniform elastic modulus E = 206 GPa, Poisson’s ratio υ = 0.3 and density ρ = 7850 kg/m³. Thickness, width and length are t = 0.7 mm, b = 15.85 mm and L = 180 mm, respectively. Yield strength is 360 MPa. Based on the previous study, axial versine displacement and DR are applied in the quasi-static explicit buckling analysis. Post-buckling paths of the thin-walled angle clamped at both ends under different computation time are contrasted with the results obtained by implicit nonlinear analysis with modified Crisfield arc-length method (MC-ALM).

Based on model analysis, the first frequency of the angle equals 752.28 Hz. The damping coefficients α = 2667.4 and β = 3.749e⁻⁴. Curves of axial load-ver- tical displacement in the case of different computation time 10T, 20T, 40T, 80T, 160T are plotted in Figure 12. It is shown that the curves from explicit solution approach to that from implicit solution with the extending of computation time.

Before the critical load, the curves obtained by explicit FEM and implicit MC- ALM are well coherent, while, after the critical load, the critical load determined by the explicit FEM is higher than that from implicit MC-ALM, as shown in Figure 12(a). The explicit consequence is close to the implicit result when the displacement extends to 0.3 mm and the computation time exceeds 40T, as shown in Figure 12(b). The kinetic energy shown in Figure 13 indicates that it upgrades rapidly when the axial load exceeds the ultimate point and the sudden enlargement of displacement occurs. This is the main factor that causes the different results by explicit method relative to implicit method.

At the point of critical load, stress and deformation contours obtained by implicit method and explicit solutions in the case of various computation time are shown in Figure 14 and Figure 15. Larger deformation can be revealed by explicit FEM when the computation time exceeds 160T. If the dynamic responses induced by rapid buckling are well controlled, the post-buckling path of large deformation can be precisely indicated by the explicit FEM. The convergence time and buckling load based on various solutions are listed in Table 1, where

the buckling load obtained by the explicit solution under 160T is close to the implicit result. However, the convergence time for the determination of buckling loads in explicit solution is far less than that in the implicit solution.

The overall buckling path of the thin-walled steel angle solved by EFG explicit method is similar to the results by FEM explicit method, as shown in Figure 16. Initially the buckling load determined by EFG is a little smaller than that from FEM as shown in Figure 16(a), but with the development of deformation, the post-buckling load obtained by EFG explicit method is higher than other results, as shown in Figure 16(b). From Table 2, it can be seen that the buckling load from EFG explicit method is smaller than other results, and the convergence time is higher than the FEM explicit solution.

4.3. Example 3: Locally Buckling of Thin-Walled Cylindrical Shell

A thin-walled cylindrical shell subjected to symmetrically concentrated force at midpoint is shown in Figure 17. The shell is assumed to be made of steel with yield stress σ = 360 MPa, elastic modulus E = 206 GPa, Poisson’s ratio υ = 0.3 and density ρ = 7850 kg/m³. Length, radius and thickness are L = 200 mm, R = 100 mm and t = 1 mm, respectively. Top and bottom are pinned-supported. The node arrangement in EFG and mesh in FEM are shown in Figure 18(a). Displacement loads, 0 - 100 mm, are applied on nine nodes as shown in Figure 18(b).

As a contrast, implicit MC-ALM has been conducted first. The natural frequency, ω = 1339.5Hz, is known by model study. Quasi-static explicit buckling analysis with versine load, DR and different computation time of 10T, 20T, 40T, 80T and 160T, is carried out. The curves of applied load with respect to radial displacement are drawn in Figure 19. It can be seen that the curve solved by implicit MC-ALM has no descending stage since non-convergence is generated for locally large deformation at the extreme point. But the post-buckling path can be traced by quasi-static explicit FEM. With the increase of computation time, the curve from explicit solution with 160T is almost coherent with implicit solution before the peak load. The kinetic energy stays in a very low level, as shown in Figure 20, which exhibits very weak dynamic response in the quasi-static solution.

The locally post-buckling contours in the case of different computation time are shown in Figure 21. It manifests that the post-buckling behavior of locally large deformation can be clearly revealed by the quasi-static explicit solution with sufficient computation time. The comparisons of convergence rate and buckling loads listed in Table 3 demonstrate the efficiency and accuracy of the quasi-static explicit method in bucking analysis.

Load-displacement curves solved by EFG and FEM explicit solutions are compared with the implicit solution, as shown in Figure 22. Before the peak point, there is little difference among the three curves. Critical load determined by EFG is closer to the implicit result except earlier descending segment occurred. The load-displacement curve from explicit EFG method is also smoother than that from the implicit FEM, but longer time is required to converge to determine buckling loads, as listed in Table 4.

5. Conclusion

Quasi-static explicit FEM and EFG method for tracing post-buckling path of thin-walled members are discussed in this study. Three examples, including snap-through, overall and local buckling, are analyzed by implicit, explicit FEM and EFG method, separately. It is found that the explicit methods with versine displacement loading, dynamic relaxation and sufficient computation time, such as more than 160T, offer almost the same consequences as the implicit FEM. The consumption of convergent time in explicit solutions is greatly less than that in implicit FEM, and there is no convergence problem. The post-buckling behavior can be easily uncovered by explicit methods. Since explicit methods are greatly affected by dynamic responses, the ratio of kinetic energy to internal energy should be controlled within 10%. Applied force loading in the analysis generates larger deformation, which can cause the dynamic response out of control. EFG quasi-static explicit buckling analysis presents as the same precision as FEM, and the load-displacement curves are smoother than those from FEM. However, it takes longer computation time than explicit FEM.

Acknowledgements

The authors are grateful for the sponsorships by State Key Laboratory of Structural Analysis for Industrial Equipment (Grant S14204), Liaoning Provincial Program for Science and Technology (Grant 2014028004), and the State Key Development Program for Basic Research of China (Grant 2015CB057300).

Cite this paper

Huang, L.H., Li, B. and Wang, Y.F. (2017) FEM and EFG Quasi-Static Explicit Buckling Analysis for Thin- Walled Members. Open Journal of Civil Engineering, 7, 432-452. https://doi.org/10.4236/ojce.2017.73030

References