Journal Search Engine

Download PDF Export Citation Korean Bibliography
ISSN : 1226-525X(Print)
ISSN : 2234-1099(Online)
Journal of the Earthquake Engineering Society of Korea Vol.27 No.2 pp.77-82
DOI : https://doi.org/10.5000/EESK.2023.27.2.077

Automated Finite Element Mesh Generation for Integrated Structural Systems

Chongyul Yoon1)*
1)Professor, Department of Civil Engineering, Hongik University
*Corresponding author: Yoon, Chongyul E-mail: cyoon@hongik.ac.kr
December 14, 2022 December 26, 2022 December 27, 2022

Abstract


The structural analysis module is an essential part of any integrated structural system. Diverse integrated systems today require, from the analysis module, efficient real-time responses to real-time input such as earthquake signals, extreme weather-related forces, and man-made accidents. An integrated system may also be for the entire life span of a civil structure conceived during the initial conception, developed throughout various design stages, effectively used in construction, and utilized during usage and maintenance. All these integrated systems’ essential part is the structural analysis module, which must be automated and computationally efficient so that responses may be almost immediate. The finite element method is often used for structural analysis, and for automation, many effective finite element meshes must be automatically generated for a given analysis. A computationally efficient finite element mesh generation scheme based on the r-h method of mesh refinement using strain deviations from the values at the Gauss points as error estimates from the previous mesh is described. Shape factors are used to sort out overly distorted elements. A standard cantilever beam analyzed by four-node plane stress elements is used as an example to show the effectiveness of the automated algorithm for a time-domain dynamic analysis. Although recent developments in computer hardware and software have made many new applications in integrated structural systems possible, structural analysis still needs to be executed efficiently in real-time. The algorithm applies to diverse integrated systems, including nonlinear analyses and general dynamic problems in earthquake engineering.



통합 구조 시스템의 유한요소망 형성의 자동화

윤 종열1)*
1)홍익대학교 공과대학 건설환경공학과 교수

초록


    1. Integrated Structural Systems

    Structural analysis module is an essential part of any integrated structural system. Diverse integrated systems of today require from the analysis module, efficient real time responses to real time input such as earthquake signals, extreme weather related forces, and man-made accidents. An integrated system may also be for the entire life span of a civil structure conceived during the initial conception, developed throughout various design stages, effectively used in construction, and utilized during usage and maintenance. All these integrated systems’ essential part is the structural analysis module and the module must be automated and computationally efficient so that the responses may be almost immediate in terms of real time. The most widely used method of structural analysis is the finite element method and for automation of the method, many effective finite element meshes must be automatically generated for a given analysis.

    2. Automated Finite Element Structural Analysis Algorithm

    The finite element method with many variations is used diversely in applications of computational mechanics and structural analysis; and today’s modern integrated structural engineering system needs a reliable, robust and automated finite element analysis module as an essential component [1-4]. Complicated nonlinear and earthquake related structural analysis requires several finite element analyses of the same structural model. If the same mesh is used throughout the analysis, inefficiency in computation results because fine meshes are adopted where they are not needed and inaccuracy due to severely distorted elements formed during the analysis that may be processed undetected as only the starting and the finally deformed element shapes are checked in general practice. Thus, finite element meshes for complex and iterative analyses must be modified and generated many times as the analysis proceeds since the accuracy of results and computational efficiency depend on the meshes and the type of elements; generally in practice, an overly fine mesh for the entire domain is adopted but this is a poor scheme because of computational inefficiency caused by too many elements and inaccuracy caused by extremely distorted elements that are not detected [1, 5].

    The procedure for automated generation of mesh algorithm based on the r-h refinement using adaptive schemes is as follows: (1) start with a good initial mesh based on previous mesh, database, or an expert system to identify general areas that need relatively finer mesh; (2) compute error for elements; (3) refine meshes shifting from the r method and the h method; and (4) repeat until preset tolerance is reached. The scheme uses computationally simple error estimation based on strain deviations. The r-h refinement optimally combines the r method which moves an existing node and the h method which divides an element into many smaller elements of same shape using a dispersion parameter. The procedure has a check for distortion limit for each element using shape factor that is computed easily for a particular element and shape [6-8].

    2.1 Strain Deviations as Error Estimates

    Gauss point strain deviations are readily available as the strains at Gauss points have to be computed for numerical integration when forming the element stiffnesses; these values used as error estimates are computationally efficient estimates for mesh refinement purposes. For planar problems using quadrilateral elements, element i`th strain deviations may be expressed as :

    e i k = j = 1 n g k ( j k k * ) 2 n g k 1 , k = x , y , x y
    (1)

    e i = { e i x + e i y + e i x y } × A i A t
    (2)

    x a = i = 1 n a x c i ( e i ) i = 1 n a ( e i )
    (3)

    y a = i = 1 n a y c i ( e i ) i = 1 n a ( e i )
    (4)

    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, ngk is Gauss points in k direction, jk is the k direction strain of Gauss point j, 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 the element, Ai is the area of element i, and At is the total area. The error of element i is normalized with respect to the total area and the relative order of errors for all elements are calculated to identify elements to be changed adaptively. Past studies have shown that similar procedures are computationally efficient [6, 9].

    2.2 The r-h Method of Mesh Refinement and Shape Factors

    Mesh generation with strain deviations as relative error for elements is formulated by alternating the r method and the h method of mesh refinement. The r method moves node at coordinates (x,y) to a new adapted coordinates xa , ya :

    x a = i = 1 n a x c i ( e i ) i = 1 n a ( e i )
    (5)

    y a = i = 1 n a y c i ( e i ) i = 1 n a ( e i )
    (6)

    Here, xci , yci are the x, y coordinates of the center of element i, and na is the number of elements sharing the node. For boundary nodes, Eqs. (5) and (6) should yield x, y coordinates of the closest point on the boundary. Distortion of an element shape violating the tolerable limit is checked by the shape factor for the element. Shape factor of a quadrilateral element i with boundary length Li is defined as :

    S i = A i 0.25 L i
    (7)

    The above factor’s maximum value is 1 which is for a square. Shape factors for distorted quadrilateral becoming triangular are less than 0.8285(shape factor value for a right equilateral triangle). A reliable shape of a quadrilateral element has approximately equal side lengths and angles; the factors for these are close to 1. Meshes with all the shape factors near 1 throughout the analysis produce reliable results. Shape factors between 0.9 and 1.0 are appropriate but the r method should restrict node movements with a limitation on the shape factor to be above 0.95.

    The h method subdivides a shape into smaller similar shapes. The elements to be selected are based on the discretization parameter d defined as:

    d = α · m e a n [ e i n i t i a l ] P max
    (8)

    Here, α is constant value to be set, m e a n [ e i n i t i a l ] is the mean value of strain, and Pmax is the maximum value of the applied load or the inertia force. A parametric investigation suggests a value between 14.0 and 15.0 for α when analyzing dynamics or earthquake engineering analysis problems. Departure form this range for α causes computational inefficiency for no increase in accuracy [8].

    The power of the r-h method comes from effective combining strategies of the r method and the h method. Many means of combination strategies have been studied [8, 9]. An optimal combination scheme of the r-method and the h-method is obtained by using dispersion parameter D defined as the difference between the mean of the normalized strains and the mode of the distribution of normalized strains:

    D = | m e a n [ n o r m a l i z e d e i ] m o d e [ n o r m a l i z e d e i ] |
    (9)

    Value of D determines alteration ratio of the r method and the h method. A value around 19 sets the ratio of the r-method and the h-method to be about 3; a higher value reduces this ratio, and a lower value increases the ratio. The h method is much more efficient in reducing overall error as the element sizes are proportionately reduced and the number of element increases in single iteration from 1 element to 4, 16, or higher. However, the r method is needed to create new shapes and to sort out overly distorted elements. The preset tolerance, typically around 0.005%, ends mesh refinement; the tolerance is the change in the sum of the strain values.

    3. A Standard Dynamic Analysis Problem

    A cantilever beam loaded at the free end is analyzed using the finite element method in the time domain using the automated algorithm outlined in Section 2.

    3.1 Time Domain Analysis

    The iteration formulas for direct numerical integration in the time domain for a typical dynamic analysis based on the Newmark-β method are [10]:

    u i + 1 = u i + ( Δ t ) u ˙ i + ( Δ t ) 2 [ ( 1 2 β ) u ¨ i + β u ¨ i + 1 ]
    (10)

    and

    u ˙ i + 1 = u ˙ i + ( Δ t ) [ ( 1 γ ) u ¨ i + γ u ¨ i + 1 ]
    (11)

    Here, ui , u ˙ i , and u ¨ i are the ith step displacement, velocity, and acceleration vectors and ui+1 , u ˙ i + 1, and u ¨ i + 1 are the i+1th step corresponding quantities; Δt is the time step size; β and γ are parameters where for efficient convergence, 1/4 and 1/2 values are recommend.

    In matrix form, equilibrium equations for the i+1th step becomes:

    K u i + 1 = F i + 1
    (12)

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

    K = K + 1 β ( Δ t ) 2 M
    (13)

    F i + 1 = F i + 1 + M β ( Δ t ) 2 [ u i + ( Δ t ) 2 u ˙ i + ( 1 2 β ) ( Δ t ) 2 u ¨ i ]
    (14)

    Here, K, M are the stiffness and Mass matrices, and Fi+1 is i+1’th time step force vector.

    3.2 Response of a Cantilever Beam with Loading at the Free End

    Fig. 1 shows the dimension of the analyzed cantilever beam and the loading function P(t). The beam is a 100 cm deep rectangular section; Young's modulus E is 212×105 N/cm2; Poisson's ratio v is 0.33, and the unit mass is 7.86×10-3 kg/cm3. In the analysis, plane stress four node quadrilateral bilinear element is used. The loading P(t) is a sinusoidal load given by -500sin(2πt)N between t = 0 and 1 seconds.

    The time step Δt selected for the analysis is 0.0025 seconds, yielding 2000 steps for 5 seconds. Meshes are generated with dispersion parameter D = 19, discretization parameter d = 14, and limit on shape factor Si set to 0.98. The initial mesh is obtained from an expert system. The terminating iteration tolerance is set to 0.005%. A new mesh is generated after four time steps, i.e., at every 0.0100 seconds. Table 1 shows progression of sample strain deviation values in the initial steps of the dynamic adaptive mesh generation where the altercation of the r method and the h method are depicted. The adaptive mesh procedure satisfies the iteration tolerance limit after about nine steps.

    Fig. 2 shows a series of meshes automatically generated; Fig. 2(a) shows the initial mesh generated from an expert system, Fig. 2(b) shows the automatically generated mesh at t = 0.24 which is the time where the deflection is the maximum during forced vibration, and Fig. 2(c) shows the generated mesh at t = 4.70 which is the time where the deflection is the maximum during the free vibration. The generated meshes show the dependence of the mesh refinement on the loading (see Fig. 2(b)), deflection(see Fig. 2(c)), and highly stressed areas (fixed end part). The meshes also show the elimination of highly distorted elements.

    In common practice, a fine regular mesh is used throughout the analysis; a regular mesh of identical square shapes totaling 1600 elements are used to simulate this and the solution from this is the general solutions. For more accurate engineering solution without an adaptive scheme, a finer regular mesh with 6400 identical elements are generated by dividing each element used in the general solution into four identical square elements. The obtained solution using this refined mesh is the engineering solution. The solution from the automated adaptive meshes is the strategy solution where the number of elements are between 560 and 826. A 64 bit Personal Computer with Intel Core 17-6700 CPU, 32.0 GB RAM, and Windows 10Pro is used to run the analyses. Vertical displacement of the free end and normal stresses at the fixed end comparisons among the engineering, the general, and the strategy solutions are shown in Fig. 3. Close agreements among the three solutions are shown in the figure; however, numeric data depict that if the engineering solutions are assume to have no error, errors from the general solutions are much bigger than the errors from the strategy solutions. Table 2 shows the comparative computation times and errors. The total error is defined as the square root of the sum of errors at the selected points. The engineering solutions are assumed exact, i.e., no error. The data on the table show that with the automation, errors have been reduced remarkably (3.592% to 0.321% for total displacement and 2.822% to 0.222% for stress) with 62.96% decrease in real run time ; the finest mesh for the entire structure used in the engineering solution needed much longer real run time increase of 266.66%. Even as the computing power of computers is increasing continuously, automation still needs efficiency of computation algorithms in every step of the procedure in order for the method to be practical and to produce realistic real time responses.

    4. Conclusions

    An automated finite element mesh generation for finite element analysis of a structures is described. The algorithm is efficient in terms of real time without any apparent loss of accuracy in responses. The procedure optimally combines well known finite element related concepts such as the h and r method of mesh refinement, shape factors for robustness of element shapes, and strain deviations for error estimates; the computational requirements for these concepts are minimal. The cantilever dynamics problem considered as an example shows the specifics of the procedure. The procedure may be applied to the analysis module in any integrated system or for any part of nonlinear structural analysis and general dynamic analysis of a structure. The procedure is adaptable for real time numerical computation responses for finite element models of large complicated structures subjected to real time dependent loads such as earthquakes and extreme weather conditions as efficient automated analyses of these dynamic and nonlinear problems are an essential part of today’s integrated systems. Addition of a powerful expert system that generates more effective initial mesh should improve the performance of automation.

    / ACKNOWLEDGMENTS /

    This work was supported by 2020 Hongik University Research Fund.

    Figure

    EESK-27-2-77_F1.gif

    Cantilever Beam and Loading P(t)

    EESK-27-2-77_F2.gif

    Automatically Generated Finite Element Meshes

    EESK-27-2-77_F3.gif

    Analysis Responses

    Table

    A Progression of Strain Deviation Values

    Run Times and Error

    Reference

    1. Zienkiewicz OC, Taylor RL, Zhu JZ. The Finite Element Method: Its Basis and Fundamentals, 6th ed. Oxford: Butterworth-Heinemann; 2005. 689 p.
    2. Bathe KJ. Finite element procedures in engineering analysis. Englewood Cliffs: Prentice- Hall; c1982. 735 p.
    3. 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.
    4. 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.
    5. Belytschko T, Hughes JR, Bathe KJ. Finite element procedures. Englewood Cliffs: Prentice-Hall; c1996; 999 p.
    6. Yoon C. An Automated Adaptive Finite Element Mesh Generation for Dynamics. EESK J Earthquake Eng. 2019;23(1):83-88.
    7. Jeong YC, Yoon C. Representative Strain Value Based Adaptive Mesh Generation for Plane Stress. Hongik J Science and Tech. 2003 Dec;7:71-86.
    8. Yoon C. Adaptive mesh generation for dynamic finite element analysis. J Kor Soc Civil Eng. 2005 Jun;25(6A):989-998.
    9. Yoon C. Adaptive finite element mesh generation for dynamic planar problems. J Kor Soc Hazard Mitigation. 2014 Aug;14(4): 43-49.
    10. Newmark NM. A method of computation for structural dynamics. J Eng Mech Div. ASCE. 1959 Mar;85(EM3):67-94.
    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