Modelling and Simulation in Engineering

Volume 2017, Article ID 1797561, 7 pages

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

^{1}Fidesys LLC, Office 402, 1 Bld. 77, MSU Science Park, Leninskie Gory, Moscow 119234, Russia^{2}Faculty of Mechanics and Mathematics, Lomonosov Moscow State University, GSP-1, 1 Leninskiye Gory, Main Building, Moscow 119991, Russia^{3}Tver State University, 33 Zhelyabov St., Tver 170100, Russia

Correspondence should be addressed to Konstantin Zingerman; ur.relbmar@namregniz

Received 28 August 2017; Accepted 13 November 2017; Published 6 December 2017

Academic Editor: Ricardo Perera

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.

#### 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.