Journal Search Engine

Volume/Issue :

Year(s) : to

Search :

Title :

Author :

Keyword :

Abstract :

Figure :

Table :

Reference :

Search Advanced Search

Adode Reader(link)

Download PDF Export Citation Korean Bibliography

ISSN : 1226-525X(Print)
ISSN : 2234-1099(Online)

Journal of the Earthquake Engineering Society of Korea Vol.23 No.1 pp.83-88
DOI : https://doi.org/10.5000/EESK.2019.23.1.083

An Automated Adaptive Finite Element Mesh Generation for Dynamics

Chongyul Yoon¹⁾*

¹⁾Professor, Department of Civil Engineering, Hongik University

Corresponding author: Yoon, Chongyul E-mail: cyoon@hongik.ac.kr

Received December 13, 2018 Review December 17, 2018 Accepted December 18, 2018

Abstract

Structural analysis remains as an essential part of any integrated civil engineering system in today’s rapidly changing computing environment. Even with enormous advancements in capabilities of computers and mobile tools, enhancing computational efficiency of algorithms is necessary to meet the changing demands for quick real time response systems. The finite element method is still the most widely used method of computational structural analysis; a robust, reliable and automated finite element structural analysis module is essential in a modern integrated structural engineering system. To be a part of an automated finite element structural analysis, an efficient adaptive mesh generation scheme based on R-H refinement for the mesh and error estimates from representative strain values at Gauss points is described. A coefficient that depends on the shape of element is used to correct overly distorted elements. Two simple case studies show the validity and computational efficiency. The scheme is appropriate for nonlinear and dynamic problems in earthquake engineering which generally require a huge number of iterative computations.

Key Words : Automated finite element analyses , Adaptive mesh generations , Gauss point error estimates , Finite element dynamics

초록

키워드 :

This article has been cited by 0 article in crossref

Cited-By

Funding:

Hongik University

1. Introduction

Complex dynamic and nonlinear analyses of structural models often requires many finite element analyses of the same structural model. If the same mesh is used throughout an analysis, computational inefficiency results due to fine mesh where they are not needed and inaccuracy due to overly distorted elements formed during the analysis that may proceed undetected as only the initial and the final element shapes are generally checked in practice. Thus, finite element meshes for these types of analyses must be dynamically adaptive and computationally efficient as results depend on the mesh and the element types used [5-7].

An efficient adaptive mesh generation scheme based on R-H refinement [7] is described. The scheme uses computationally simple error estimation based on representative strains [8]. A coefficient depending on the element shape [9] is used to correct overly distorted elements. The R-H refinement optimally combines the r-method (moving an existing node) and the h-method (dividing an element into smaller elements) using a dispersion parameter [8, 9]. The scheme includes a check for limiting distortion for each element using shape factor that is easily computed for a particular shape [8-10]. The dynamic analysis considered is in time domain, based on direct integration. The results are based on the basic procedure outlined in reference [11] but the data are from diverse refinements and improvements on the software and tuning of parameters. Two simple case studies for planar problems using 4-node quadrilateral elements show the validity and computational efficiency. The scheme is appropriate for nonlinear and dynamic problems in earthquake engineering which generally require a huge number of iterative computations.

2. Time Domain Finite Element Dynamic Analysis

The direct numerical integration iterative formulas in the time domain for a typical dynamic structural analysis based on the widely used Newmark-β method may be summarized as follows [5, 12]:

u_{i + 1} = u_{i} + (Δ t) {\dot{u}}_{i} + {(Δ t)}^{2} [(\frac{1}{2} - β) {\ddot{u}}_{i} + β {\ddot{u}}_{i + 1}]

(1)

{\dot{u}}_{i + 1} = {\dot{u}}_{i} + (Δ t) [(1 - γ) {\ddot{u}}_{i} + γ {\ddot{u}}_{i + 1}]

(2)

Here, $u_{i}, {\dot{u}}_{i}, and {\ddot{u}}_{i}$ are the displacement, velocity, and acceleration vectors in the ith step and $u_{i + 1}, {\dot{u}}_{i + 1}, are {\ddot{u}}_{i + 1}$ are the similar quantities in the i+1th step; Δt is the time step duration; β and γ are parameters selected by the user and the recommended values of β=1/4 and γ=1/2 are commonly used for efficient convergence.

The matrix equilibrium equations for the i+1th step may be represented as follows:

K^{'} u_{i + 1} = {F^{'}}_{i + 1}

(3)

The matrices K′ and F′_i+1 represent the following:

K^{'} = K + \frac{1}{β {(Δ t)}^{2}} M

(4)

{F^{'}}_{i + 1} = F_{i + 1} + \frac{M}{β {(Δ t)}^{2}} [u_{i} + {(Δ t)}^{2} {\dot{u}}_{i} + (\frac{1}{2} - β) {(Δ t)}^{2} {\ddot{u}}_{i}]

(5)

Here, K is the stiffness matrix, M is the mass matrix, and F_i+1 is the force vector for the i+1th time step.

3. Mesh Errors and Refinements

3.1 Relative Gauss Point strains as Error Estimates

An adaptive mesh generation scheme automates mesh generation by repeated improvement of the previous mesh based on error estimations of the previous mesh and improvement strategies based on the errors [10,13]. Relative strains at Gauss points which are easily computed since the strains here need to be computed for the computation of the element stiffnesses; these relative and representative error estimates are proven to be computationally efficient error estimates for mesh refinement purposes [10].

For planar problems using quadrilateral elements, the representative strain value of element i may be represented by the following equations:

{‖ e ‖}_{i k} = \sqrt{\frac{\sum_{j = 1}^{n_{g k}} {(\in_{j k} - \in_{k}^{*})}^{2}}{n_{g k} - 1}}, k = x, y, x y

(6)

{‖ e ‖}_{i} = {{‖ e ‖}_{i x} + {‖ e ‖}_{i y} + {‖ e ‖}_{i x y}} \times \frac{A_{i}}{A_{t}}

(7)

Here, ${‖ e ‖}_{i k}$ is the norm of the k directional (x, y and xy for planar problems where xy is the shear component) standard deviation of the strain, n_gk is the number of Gauss points in the k direction, є_jk is the k directional strain of Gauss point j, $\in_{k}^{*}$ is the k directional strain, and ${‖ e ‖}_{i}$ is the norm of the representative strain of element i which is used as error for element I. In addition, A_i is the area of element i and A_t is the total area of the entire mesh. In Eq. (7), the representative error of element i is weighted by the element's area with respect to the total area and thus the relative order of errors for all elements are computed to identify elements to be adaptively changed. Previous studies have shown that this is computationally very efficient for achieving this [11].

3.2 The R-H Method of Mesh Refinement with Dispersion and Shape Factors

The adaptive mesh generation based on the error estimation with representative strains is formulated combining the r-method and the h-method. The r-method moves existing node at coordinates (x,y) to the new adapted coordinates x_a, y_a :

x_{a} = \frac{\sum_{i = 1}^{n_{a}} x_{c i} ({‖ e ‖}_{i})}{\sum_{i = 1}^{n_{a}} ({‖ e ‖}_{i})}

(8)

y_{a} = \frac{\sum_{i = 1}^{n_{a}} y_{c i} ({‖ e ‖}_{i})}{\sum_{i = 1}^{n_{a}} ({‖ e ‖}_{i})}

(9)

Here, x_ci, y_ci are the coordinates of the centroid of element i, and n_a is the number of elements sharing the considered node. For nodes on the boundary, Eqs. (8) and (9) are moved to the closest point on the boundary. Possible distortion of an element shape beyond the tolerable limit is checked by limiting the shape factor. Shape factor S_i of a quadrilateral element i with total boundary length L_i is defined as follows [8]:

S_{i} = \frac{\sqrt{A_{i}}}{0.25 L_{i}}

(10)

The shape factor is defined so that the maximum value is 1 which is for a square. Shape factors for quadrilaterals approaching a triangular shape are less than around 0.8285 (shape factor for a right equilateral triangle). A robust shape of a quadrilateral element has approximately equal side lengths and angles. The shape factors of these are close to 1. An attempt to keep shape factors close to 1 throughout the analysis produces reliable results. Generally, shape factors between 0.9 and 1.0 are acceptable and thus the r-method restricts node movements with a limitation on the shape factor to be above 0.95.

The h-method divides an element into smaller elements of the same shape. The elements to be subdivided are based on the discretization parameter d defined as follows:

d = \frac{α \cdot m e a n [{‖ e ‖}_{i n i t i a l}]}{P_{max}}

(11)

Here, α is a constant, $m e a n [{‖ e ‖}_{i n i t i a l}]$ is the mean representative strain value of the mesh and P_max is the maximum value of the applied load. A parametric study recommends a value between 14.0 and 15.0 for optimal value of α for general dynamics problems. Significant deviations of α from this range result in computational inefficiency for no apparent increase in accuracy [8, 11].

The effectiveness of the R-H method depends on strategies that combine the r-method and the h-method. Various means of combining have been studied [10, 13]. To obtain an optimal combination of the r-method and the h-method, first the representative strain values are normalized at each step so that the minimum and maximum values are from 0 to 100. A dispersion parameter D is defined as the absolute value of the difference between the mean of the normalized representative strains and the mode of the distribution of normalized representative strains:

D = | m e a n [n o r m a l i z e d {‖ e ‖}_{i}] - \mod e [n o r m a l i z e d {‖ e ‖}_{i}] |

(12)

A constant need to be set for D to shift between the r-method and the h-method. A reasonable value for this constant is around 18-20 where if D is larger than this value, the r-method is used, and in other cases, the h-method is used. A value around 19 sets the ratio of the r-method and the h-method to be about 3. A higher value for this constant will reduce this ratio. Generally, the h-method is much more efficient in reducing overall error as the element sizes are proportionately reduced (e.g. 1 element to 4 or 16 elements of the same shape reduced to 1/4 or 1/16 of the original size). However, the r-method is essential in creating new shapes and orientations that are not overly distorted by using the shape factors given by Eq. (10). The R-H method of adaptive refinement of the mesh is terminated when the change in the sum of the representative strain values is less than the preset tolerance, typically set at 0.01%.

4. Planar Problems using Four Node Quadrilateral Elements

4.1 A Portal Example

Fig. 1 shows the sample two dimensional 600 cm wide and 600 cm high portal steel portal frame. The beam and columns are all rectangular 100 cm deep steel sections where Young's modulus is 210×10⁵ N/cm², Poisson's ratio is 0.3, and the unit mass is 7.85×10^-3 kg/cm³. Element type is the four node isoparametric quadrilateral bilinear plane strain element. The only applied dynamic load is a concentrated load P at the center of the frame (point B in Fig. 1) where between time t=0 and t=1 second, one period of a sinusoidal load given by -500sin(2πt)Newton.

Free vibration response continues after the steady state response (1 second), and the total response time considered is 5 seconds. The time step Δt selected for analysis is 0.005 seconds, yielding 1000 steps for 5 seconds. The adaptive meshes are generated with dispersion parameter D = 19, discretization parameter d = 14, and the limit on shape factor S_i is set to 0.95. The initial mesh to start the adaptive algorithm is a regular mesh composed of 144 identical square elements. The tolerance set to terminate the adaptive refinements is set at 0.01% for change in the sum of the representative strain from the previous step. A new mesh is generated at every four time steps, i.e., at every 0.02 seconds. Table 1 shows a set of sample representative strain values in the initial steps of the dynamic adaptive mesh where the transitions of the r-method and the h-method are demonstrated. Overall, the adaptive mesh scheme satisfies the preset tolerance after about ten adaptive steps.

Fig. 2 shows the adaptive mesh at 0.275 seconds when the vertical deflection of point A is at the maximum which is during the steady state response; the number of elements generated here is the maximum at 840. Fig. 2 also shows the mesh when the number of elements is the minimum at 760, occurring at t=3.760.

The solution obtained by the dynamic adaptive meshes is named the strategic solution. In general practice, a regularly discretized mesh is used throughout the analysis; a regular mesh of square shapes with 576 identical elements are used to simulate this and the solution from this is called the general solutions. To estimate a more accurate engineering solution without an adaptive scheme, a finer regular mesh with 2304 elements is generated by dividing each element used in the general solution into four identical square elements. The solution obtained by using this fine mesh is the engineering solution. The finite element program is run on 64 bit Personal Computer with Intel (R) Pentium (R) CPU3240 @3.10GHz, 4.0 GB RAM, and Windows 10 Home Version 1803. Figs. 3, 4 and 5 show comparisons of the vertical displacement of point A, the horizontal (x directional) normal stress of point C, and the vertical (y directional) normal stress of point C of the engineering, the general, and the strategy solutions. The figures show close agreement among the three solutions. However numeric values show that if the engineering solutions are assume to be the most accurate, the errors from the general solutions are much larger than the errors from the strategic solutions. Table 2 shows the comparative computation times and errors among the engineering, the general, and the strategic solutions. The total error is computed as the square root of the sum of the errors at selected critical points (A, B, C in Fig. 1) where the engineering solutions are assumed to have no error. The data on the table show that with the dynamic adaptive scheme, errors have been reduced significantly (3.59% to 0.38% for total displacement and 2.82% to 0.24% for stress) with a 259.1% increase in real run time whereas the fine mesh for the entire structure for the engineering solutions required enormous increase in real run time of 562.93%. Even as the computing speed of computers is continuously increasing in general, the enormous amount of computing needed in the more accurate and complex analysis of structures requires computing efficiency in every aspect of the algorithm in order for the method to be practical [22, 23].

4.2 A cantilever beam example

Fig. 6 shows a series of meshes generated for a cantilever beam loaded with a concentrated load at the free end with time variation equal to P(t) in the portal sample in Section 4.1 (i.e. one period of a sinusoidal load given by -500sin(2πt) Newton). The figure shows overlap of undeformed previous with deformed refined meshes to depict how the generated meshes need to correct overly distorted elements as the structure deforms and to generate effective meshes at each time step as the analysis proceeds.

5. Conclusions

An efficient adaptive mesh generation scheme for dynamic analyses of structures in time domain is described. The scheme uses representative strain values from each element computed from the previous time step for estimation of error. Based on a dispersion parameter, the r-method and the h-method are combined optimally for mesh refinement. Shape factors used to correct overly distorted elements. Analyses of the applications of the scheme for two case studies show that the representative strain values which efficiently compute relative errors among elements successfully identify elements to be refined and that the algorithm is efficient computationally and have reasonable confidence level for accuracy for dynamic problems in structural analysis. The algorithm dynamically produces robust, effective (small errors) and computationally efficient meshes for finite element analyses. The scheme is applicable to other general nonlinear analyses where requirements for generated meshes are similar to those required by the case studies. The scheme is especially appropriate for real time computations of large complex structures under erratic time dependent loads such as earthquakes and turbulent winds. Efficient automated analyses of these dynamic and nonlinear problems are an essential part of today’s integrated structural engineering system. Addition of a powerful expert system that generates a reasonable initial mesh that starts the proposed adaptive mesh generation scheme will make the algorithm computationally more efficient.

/ ACKNOWLEDGMENTS /

This work was supported by 2014 Hongik University Research Fund.

Figure

Fig. 1.

Geometry of the portal example

Fig. 2.

Maximum and minimum number of elements

Fig. 3.

Vertical displacement of bottom mid point A of the frame

Fig. 4.

Horizontal normal stress at mid point C of the frame

Fig. 5.

Vertical normal stress at mid point C of the frame

Fig. 6.

Deformed and undeformed meshes in the catilever beam example (Left fixed, Concentrated load at the free end)

Table

Table 1.

A sample set of progression of representative strain values

Table 2.

Comparative computation times and error

Reference

Rabizuk T, Borda S, Zi GA. Three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction in statics and dynamics. Comp Mech. 2007 Jan;40:473-495.
Zi G, Song JH, Lee SH. A new method for growing multiple cracks without remeshing and its application to fatique crack growth. J Kor Soc Civil Eng. 2005 Jan;25:183-190.
Alfarra B, Lopez-Almansa F, Oller S. New methodology for calculating damage variable evolution in plastic damage model for RC structure. Eng Struct. 2017 Mar;132:70-86.
Snozzi L, Molinari JF. A cohesive element model for mixed mode loading with frictional mixed mode loading with frictional contactcapability. IJNME. 2013 Jan;93(5):510-526.
Bathe KJ. Finite element procedures in engineering analysis. Englewood Cliffs: Prentice- Hall; c1982. 735 p.
Belytschko T, Hughes JR, Bathe KJ. Finite element procedures. Englewood Cliffs: Prentice-Hall; c1996. 999 p.
Choi C, Jung H. Adaptive mesh generation for dynamic finite element analysis. J Kor Soc of Civil Eng. 1998 Mar;18(I-2):203-804.
Jeong YC, Yoon C. Representative Strain Value Based Adaptive Mesh Generation for Plane Stress. Hongik J Science and Tech. 2003 Dec;7:71-86.
Yoon C. Adap tive mesh generation f or dynamic f inite element analysis. J Kor Soc Civil Eng. 2005 Jun;25(6A):989-998.
de Las Casas EB. R-H mesh improvement algorithms for the finite element Method. Ph.D. Dissertation. West Lafayette: Purdue University. c1988. 389p.
Yoon C. Adaptive finite element mesh generation for dynamic planar problems. J Kor Soc Hazard Mitigation. 2014 Aug;14(4):43-49.
Newmark NM. A method of computation for structural dynamics. J Eng Mech Div. ASCE. 1959 Mar;85(EM3):67-94.
Zienkiewicz OC, Taylor RL, Zhu JZ. The Finite Element Method: Its Basis and Fundamentals, 6th ed. Oxford: Butterworth-Heinemann; c2005. 689 p.

Journal Abbreviation J. Earthq. Eng. Soc. Korea
Frequency Bimonthly
Doi Prefix 10.5000/EESK
Year of Launching 1997
Publisher Earthquake Engineering Society of Korea

Indexed/Tracked/Covered By