The Implementation of Spectral Element Method in a CAE System for the Solution of Elasticity Problems on Hybrid Curvilinear Meshes

Konovalov, Dmitriy; Vershinin, Anatoly; Zingerman, Konstantin; Levin, Vladimir

doi:https://doi.org/10.1155/2017/1797561

Modelling and Simulation in Engineering

On this page

Abstract Introduction Materials and Methods Results Discussion Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2017 | Article ID 1797561 | https://doi.org/10.1155/2017/1797561

The Implementation of Spectral Element Method in a CAE System for the Solution of Elasticity Problems on Hybrid Curvilinear Meshes

Dmitriy Konovalov,¹Anatoly Vershinin,^1,2Konstantin Zingerman ,³and Vladimir Levin²

Academic Editor: Ricardo Perera

Received28 Aug 2017

Accepted13 Nov 2017

Published06 Dec 2017

Abstract

Modern high-performance computing systems allow us to explore and implement new technologies and mathematical modeling algorithms into industrial software systems of engineering analysis. For a long time the finite element method (FEM) was considered as the basic approach to mathematical simulation of elasticity theory problems; it provided the problems solution within an engineering error. However, modern high-tech equipment allows us to implement design solutions with a high enough accuracy, which requires more sophisticated approaches within the mathematical simulation of elasticity problems in industrial packages of engineering analysis. One of such approaches is the spectral element method (SEM). The implementation of SEM in a CAE system for the solution of elasticity problems is considered. An important feature of the proposed variant of SEM implementation is a support of hybrid curvilinear meshes. The main advantages of SEM over the FEM are discussed. The shape functions for different classes of spectral elements are written. Some results of computations are given for model problems that have analytical solutions. The results show the better accuracy of SEM in comparison with FEM for the same meshes.

1. Introduction

The finite element method (FEM) [1, 2] was considered as the main approach for solution of the problems of elasticity theory taking into account finite deformations. The desire to find the stress-strain state of structures with high accuracy is still a hot topic at the moment and it forces us to seek new and effective methods for solving engineering problems. One of such approaches is the spectral element method (SEM) [3–5]. The SEM was firstly applied for the modeling of liquid flows [6]. These problems require a high accuracy and high rate of computations. Later the SEM was successfully adopted for seismic problems [3, 7, 8]. A special quadrature formula was constructed for integration over space. This quadrature formula permits one to develop a fully explicit scheme for the integration over time, which is an important advantage of the SEM.

Among the main advantages of the SEM over the FEM is the high accuracy of approximation of the solution at a substantially smaller number of mesh elements required. The error of the numerical solution of SEM decreases exponentially with the order of elements [9]. When working with a model, the user does not need to rebuild and refine the mesh to verify the mesh convergence [9] of the obtained solution, as it has been when using the FEM, since with the use of the SEM the mesh can remain initial with only changing the order of elements. The possibility of effective paralleling of computing systems with shared and distributed memory using OpenMP and MPI technology makes SEM attractive for industrial applications in various software systems [10]. In particular, a spectral element finite element scheme that efficiently solves elliptic problems on unstructured hexahedral meshes is developed in [11]. It is demonstrated that problems with over 50 million degrees of freedom can be solved in a few seconds on an off-the-shelf GPU.

The industrial application of the SEM is hindered by the problem of mesh generation [12, 13]. Typically, for a multibody model geometry it is quite difficult to build a conformal finite element mesh consisting of hexahedral elements only for which the SEM [3] was initially developed. In general, the three-dimensional finite element mesh contains inmost dominantly hexahedral and tetrahedral elements. The prismatic and pyramidal elements are used in order to connect triangular and quadrilateral surfaces of elements. A two-dimensional mesh may contain rectangular and triangular elements. Such finite element meshes are called hybrid (mixed) meshes [14]; the SEM was adopted for them within the present work.

This article discusses the variant of implementation of spectral element method on hybrid curvilinear meshes for three-dimensional problems of elasticity theory and its industrial application in CAE Fidesys [15]. The shape functions for different classes of spectral elements are written. Some results of computations are given for model problems that have analytical solutions.

2. Materials and Methods

Let the domain be divided into elements . The classes of quadrilaterals and triangles serve as elements for two-dimensional case. The classes of hexahedrons, tetrahedrons, pyramids, and prisms serve as elements for three-dimensional case. Each element is defined by reference points. The number of reference points for a particular element is indicated in its name.

Each reference point is determined by the index varying from 1 to , where is a number of reference points of the finite element. A nondegenerate mapping of the base (reference) element in is built for each element. The coordinates of points in the global coordinate system and in the reference coordinate system are related as follows:where is the vector of global coordinates, , the global coordinates of the reference point of the element , the reference coordinates, and the th shape function of the finite element, .

Generally, the Jacobian [16] of the mapping is required for the calculation of the integral of an arbitrary function over the element through reference coordinates :where and is the Jacobian in the reference point . The Jacobian matrix can be calculated in a standard way: In the three-dimensional case:In the two-dimensional case, we assume that the finite elements lie in a plane and reference points lie in a plane :

2.1. QUAD: The Class of Quadrilaterals

The class of quadrilaterals includes the following types of finite elements, named by number of reference points: four-node, QUAD4 ; eight-node, QUAD8 ; nine-node, QUAD9 . The reference element for the class of quadrilaterals The nodes of a quadrilateral spectral element are the Gauss-Lobatto-Legendre (GLL) nodes. Let be an order of the spectral element; then the number of nodes in the element In one-dimensional case the GLL-nodes are calculated as roots of the derivative of the Legendre polynomial of an order , which can be defined as . In two-dimensional case the coordinates of GLL-nodes are defined as a direct product of one-dimensional coordinates.

The shape functions for the element are built based on the direct product of one-dimensional Lagrange polynomials: is the th shape function of the order . The one-dimensional Lagrange polynomials of the order are defined as follows [3]: Each polynomial is equal to 1 at its node and is equal to 0 at the remaining nodes of the element , where is the Kronecker symbol.

The approximation of solution of the problem on a quadrilateral spectral element will be as follows :In order to integrate an arbitrary function over the element , the Gauss-Lobatto-Legendre (GLL) quadrature formula is used:where are the GLL-weights and are the coordinates of GLL-nodes.

GLL-weights are calculated as follows:The coordinates of GLL-nodes and GLL-weights (indexes vary from to ) can be written as follows:One of the most important features of the SEM is that the computation of the Lagrange polynomials that are used for the approximation of solution is based on the same GLL-nodes that are necessary to calculate integrals over the region using the Gauss-Lobatto-Legendre quadrature formula.

2.2. TRI: The Class of Triangles

The class of triangles contains the following types of finite elements named by the number of reference points: three-node, TRI3; six-node, TRI6. The reference element for the class of trianglesThe nodes of a triangular spectral element are the points , which are obtained from the solution of the electrostatic problem described by Hesthaven and Teng [17] and Taylor et al. [18]. Let be an order of the spectral element; then the number of nodes in the element .

Shape functions for the element are built based on an orthogonal basis ( is the order of a spectral element):where are the Jacobi polynomials of the order on the interval , orthogonal with the weight function . The Vandermonde matrix elements are calculated in the nodes of the element ,. Then the shape functions will be as follows:The approximation of solution of the problem on a triangular spectral element will be as follows:To integrate an arbitrary function over the element , the symmetric quadrature formulas, described by Zhang et al. [19], are used:where are positive weights, are the coordinates of points, and is the number of points in the quadrature formula corresponding to the spectral element of the order .

2.3. HEX: The Class of Hexahedrons

The class of hexahedrons contains the following types of finite elements named by the number of reference points: eight-node, HEX8; twenty-node, HEX20; twenty-seven-node, HEX27. The reference element for the class of hexahedronsThe nodes of quadrilateral spectral element are the GLL-nodes. Let be an order of the spectral element; then the number of nodes in the element .

The shape functions for the element are built based on the direct product of one-dimensional Lagrange polynomials: , th shape function of the order .

The approximation of solution of the problem on hexahedral spectral element will be as follows :To integrate an arbitrary function over the element , the Gauss-Lobatto-Legendre (GLL) quadrature formula is used:where are the GLL-weights and are the coordinates of GLL-nodes.

Coordinates of the GLL-nodes and GLL-weights (indexes vary from to )

2.4. TETRA: The Class of Tetrahedrons

The class of tetrahedrons contains the following types of finite elements named by the number of reference points: four-node, TETRA4; ten-node, TETRA10. The reference element for the class of tetrahedronsThe nodes of a tetrahedron spectral element are the points , which are obtained from the solution of the electrostatic problem described by Hesthaven and Teng [17]. Let be an order of the spectral element; then the number of nodes in the element .

The shape functions for the element are built based on the following orthogonal basis ( is the order of a spectral element):where are the Jacobi polynomials of the order on the interval , orthogonal with the weight function . The Vandermonde matrix elements are calculated in the nodes of the element , . Then the shape functions will be as follows:The approximation of solution of the problem on a tetrahedral spectral element will be as follows:In order to integrate an arbitrary function over the element , the symmetric quadrature formulas, described by Zhang et al. [19], are used:where are positive weights, are coordinates of points, and is the number of points in the quadrature formula corresponding to the spectral element of the order .

2.5. PYRAMID: The Class of Pyramids

The class of pyramids contains the following types of finite elements named by the number of reference points: five-node, PYRAMID5; thirteen-node, PYRAMID13. The reference element for the class of pyramidsThe nodes of a pyramidal spectral element are the points , matching the GLL-nodes on a square pyramid base and nodes for triangular spectral element on the side surfaces of the pyramid. Internal points are located in the planes parallel to the square base of the pyramid in the scaled GLL-nodes. Let be an order of the spectral element; then the number of nodes in the elementThe shape functions for the element are built based on the following orthogonal basis ( is the order of a spectral element):where are the Jacobi polynomials of the order on the interval , orthogonal with the weight function . The Vandermonde matrix elements are calculated in the nodes of the element . Then the shape functions will be as follows:The approximation of solution of the problem on pyramidal spectral element will be as follows:To integrate an arbitrary function over the element , the symmetric conical quadrature formulas, described by Felippa [20], are used:where are positive weights, are coordinates of points, and is the number of points in the quadrature formula corresponding to the spectral element of the order .

2.6. WEDGE: The Class of Prisms

The class of prisms contains the following types of finite elements named by the number of reference points: six-node: WEDGE6; fifteen-node, WEDGE15. The reference element for the class of prisms

The nodes of a prismatic spectral element are the points arising from the direct product of the nodes for the triangular spectral element and one-dimensional GLL-points. Let be an order of the spectral element; then the number of nodes in the element .

The shape functions for the element are built based on the product of the shape functions for the triangular spectral element and the Lagrange polynomials:The approximation of solution of the problem on a triangular spectral element will be as follows :To integrate an arbitrary function over the element , the quadrature formulas are used, which are the unions of quadrature formulas for triangular spectral elements and the Gauss-Lobatto-Legendre quadrature formula:where are weights, are coordinates of points, and is the number of points in the quadrature formula corresponding to the spectral element of the order .

2.7. Features of Implementation

The main steps of application of the spectral element method to the problems of the theory of elasticity are similar to the steps of problems solving using the finite element method, such as discretization of equilibrium equations in the integral form; selection of quadrature for calculation of integrals; building of local matrices of stiffness, mass, and damping for each element; assembling global matrices of stiffness, mass, and damping. At present, most of mesh generators for geometric models build only the finite element meshes of the first and second order, which forces us to rebuild the mesh model for calculation using the spectral element method. The most long-lasting operation at this step is to construct a graph of mesh connectivity, so one can use the connectivity graph of the initial finite element mesh model in order to speed up the algorithm. Increasing the speed of the algorithm is associated with an increase in memory consumption. In particular, for high order elements it is necessary to store the points of integration, quadrature weights, and computed values of derivatives of functions at these points. An important feature of the algorithm is the ability to use the spectral element method on any hybrid finite element mesh of the first and second order. In fact, uncoated classes of reference mesh elements currently are beam and shell elements.

3. Results

This approach to the implementation of spectral element method on hybrid curvilinear meshes for elasticity problems was industrially implemented in CAE Fidesys. Let us give the spectral element solution of the problem of stress analysis for a structural element and the results of spectral element analysis of wave propagation [13] in the elastic medium in comparison with the finite element and analytical solutions.

Numerical experiments clearly demonstrate that the mesh convergence in the model is achieved much faster when using the spectral element method as compared with the finite element method.

3.1. Calculation of Wave Propagation in Unbounded Elastic Medium under the Action of a Point Source of Disturbances

The wave propagation in unbounded elastic medium was analysed using the spectral element method, and the results were compared with the solution obtained via analytical calculation of displacement vector components, depending on time at fixed points of the body [1]. In Figure 1, one can see the time dependencies of the displacement vector components , and the relative errors , of calculation results for these components obtained via the SEM against the analytical solution. The results are given for one of the receivers within the medium; denotes time, and is the order of approximating polynomials. As it can be seen in the graphs, the maximum error of the calculation via SEM on a mesh of about 10 thousands of elements does not exceed 2.0% at the 8th order of approximating polynomials and 0.4% at the 10th order of approximating polynomials.

3.2. Analysis of Stress-Strain State of a Bridge Span under the Action of Pressure on Its Base

The calculation was performed for the finite elements of the first and second orders, as well as for the spectral elements of the orders from 1 to 4. An error of vector norm of maximum displacement of the bridge span was estimated. As a reference value for the assessment of the accuracy of the obtained solution, a value corresponding to the model mesh convergence was taken. Results are provided for the following cases: case 1, 10 thousand elements (the character element size is 0.8); case 2, 28 thousand elements (the character element size is 0.4); case 3, 104 thousand elements (the character element size is 0.2); case 4, 849 thousand elements (the character element size is 0.1). The case number is laid off as abscissa.

As it can be seen in Figure 2 on the graphs, the spectral element method on the coarsest mesh just on the 4th order demonstrates the accuracy of the order of 1%, whereas the finite element method on a mesh of 849 thousand elements still gives an error of more than 10% for the finite elements of the first order and about 2% for the finite elements of the second order.

4. Discussion

This article discusses the variant of spectral element method implementation on the hybrid curvilinear meshes for problems of elasticity theory and its industrial application in CAE Fidesys. The comparison with the finite element method was conducted, which allowed us to draw a conclusion about the high accuracy of the method and the correctness of algorithms and the program developed. In the future, it is planned to expand the implementation of the method for the cases of shell and beam elements within the static and dynamic elasticity problems.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The research for this article was performed in FIDESYS LLC and was financially supported by Russian Ministry of Education and Science (Project no. 14.579.21.0112, Project ID RFMEFI57915X0112).

References

T. J. R. Hughes, The Finite Element Method—Linear Static and Dynamic Finite Element Analysis, Prentice Hall, Englewood Cliffs, NJ, USA, 1987.
O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method. Volume 1: The Basis, Butterworth-Heinemann, Oxford, UK, 2000.
D. Komatitsch and J.-P. Vilotte, “The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures,” Bulletin of the Seismological Society of America, vol. 88, no. 2, pp. 368–392, 1998.
View at: Google Scholar
B. Schubert, “The spectral-element method for seismic wave propagation: theory, implementation and comparison to finite difference methods,” Tech. Rep., University of Munich, Munich, Germany, 2003.
View at: Google Scholar
M. Stupazzini, A spectral element approach for 3D dynamic soil-structure interaction problems [Dottorato in Ingeneria Seismica], Milan, Italy, Politecnico di Milano, 2004.
A. T. Patera, “A spectral element method for fluid dynamics: laminar flow in a channel expansion,” Journal of Computational Physics, vol. 54, no. 3, pp. 468–488, 1984.
View at: Publisher Site | Google Scholar
D. Komatitsch, C. Barnes, and J. Tromp, “Simulation of anisotropic wave propagation based upon a spectral element method,” Geophysics, vol. 65, no. 4, pp. 1251–1260, 2000.
View at: Publisher Site | Google Scholar
D. Komatitsch, J. Ritsema, and J. Tromp, “Geophysics: The spectral-element method, beowulf computing, and global seismology,” Science, vol. 298, no. 5599, pp. 1737–1742, 2002.
View at: Publisher Site | Google Scholar
V. S. Karpenko, A. V. Vershinin, V. A. Levin, and Zingerman K. M., “Some results of mesh convergence estimation for the spectral element method of different orders in Fidesys industrial package,” IOP Conference Series: Materials Science and Engineering, vol. 158, Article ID 012049, 2016.
View at: Google Scholar
K. Kannan, M. Khoshlessan, M. Herrmann, and Y. Peet, “Detailed numerical study of flow and vortex dynamics in staggered pin-fin arrays within a channel,” in Proceedings of the ASME Turbo Expo 2016: Turbomachinery Technical Conference and Exposition, 123972, June 2016.
View at: Google Scholar
J.-F. Remacle, R. Gandham, and T. Warburton, “GPU accelerated spectral finite elements on all-hex meshes,” Journal of Computational Physics, vol. 324, pp. 246–257, 2016.
View at: Publisher Site | Google Scholar | MathSciNet
M. Charara, A. Vershinin, E. Deger, D. Sabitov, and G. Pekar, “3D spectral element method simulation of sonic logging in anisotropic viscoelastic media,” SEG Technical Program Expanded Abstracts, vol. 30, no. 1, pp. 432–437, 2011.
View at: Publisher Site | Google Scholar
M. Charara, A. Vershinin, D. Sabitov, and G. Pekar, “SEM wave propagation in complex media with tetrahedral to hexahedral mesh,” in Proceedings of the 73rd European Association of Geoscientists and Engineers Conference and Exhibition, pp. 41–45, Vienna, Austria, 2011.
View at: Google Scholar
M. Bergot, G. Cohen, and M. Durufle, “Higher-order finite elements for hybrid meshes using new nodal pyramidal elements,” SIAM Journal of Scientific Computing, vol. 42, no. 3, pp. 345–381, 2010.
View at: Publisher Site | Google Scholar | MathSciNet
Fidesys LLC official website, http://cae-fidesys.com/.
J. N. Bronstein, K. A. Semendjajew, G. Musiol, and H. Muchkig, Taschenbuch der Mathematik, Auflage 4, Harri Deutch, Frankfurt, Germany, 1999.
J. S. Hesthaven and C. H. Teng, “Stable spectral methods on tetrahedral elements,” SIAM Journal on Scientific Computing, vol. 21, no. 6, pp. 2352–2380, 2000.
View at: Publisher Site | Google Scholar | MathSciNet
M. A. Taylor, B. A. Wingate, and R. E. Vincent, “An algorithm for computing Fekete points in the triangle,” SIAM Journal on Numerical Analysis, vol. 38, no. 5, pp. 1707–1720, 2000.
View at: Publisher Site | Google Scholar | MathSciNet
L. Zhang, T. Cui, and H. Liu, “A set of symmetric quadrature rules on triangles and tetrahedra,” Journal of Computational Mathematics, vol. 27, no. 1, pp. 89–96, 2009.
View at: Google Scholar | MathSciNet
C. A. Felippa, “A compendium of FEM integration formulas for symbolic work,” Engineering Computations, vol. 21, no. 8, pp. 867–890, 2004.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2017 Dmitriy Konovalov et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1698

Downloads

991

Citations