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

Authors: Dmitriy Konovalov; Anatoly Vershinin; Konstantin Zingerman; Vladimir Levin
Journal: Modelling and Simulation in Engineering
 (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.


---

## Body

## 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Ω∈R2R3 be divided into ne elements Ωe. 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 Ωe 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 indexa varying from 1 to d, where d is a number of reference points of the finite element. A nondegenerate mapping Fe:Λ→Ωe of the base (reference) element Λ in Ωe is built for each element. The coordinates of points in the global coordinate system xe=xe,ye,ze and in the reference coordinate system ξ=ξ,η,θ are related as follows:(1)xeξ=Feξ=∑a=1dNaξxae,where xe=xe,ye,ze is the vector of global coordinates, xe∈Ωe, xae the global coordinates of the reference point a of the element Ωe, ξ=ξ,η,θ the reference coordinates, and Naξ the ath shape function of the finite element, ξ∈Λ.Generally, the Jacobian [16] of the mapping Feξ is required for the calculation of the integral of an arbitrary function hx over the element Ωe through reference coordinates ξ:(2)∫Ωehxdx=∫Λhex·dxdξdξ=∫ΛheξJeξdξ,where hex=hx∣Ωe and Jeξ=det⁡Jeξ is the Jacobian in the reference point ξ. The Jacobian matrix Jeξ can be calculated in a standard way: (3)Jeξ=dxeξdξ=∑a=1ddNadξdξxae.In the three-dimensional case:(4)Jeξ=JxξJxηJxθJyξJyηJyθJzξJzηJzθ=∂x∂ξ∂x∂η∂x∂θ∂y∂ξ∂y∂η∂y∂θ∂z∂ξ∂z∂η∂z∂θ,Je-1ξ=∂ξ∂x∂ξ∂y∂ξ∂z∂η∂x∂η∂y∂η∂z∂ζ∂x∂ζ∂y∂ζ∂z.In the two-dimensional case, we assume that the finite elements lie in a plane x,y and reference points lie in a plane ξ,η:(5)Jeξ=JxξJxηJyξJyη=∂x∂ξ∂x∂η∂y∂ξ∂y∂η,Je-1ξ=∂ξ∂x∂ξ∂y∂η∂x∂η∂y.


### 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, QUAD4d=4; eight-node, QUAD8 d=8; nine-node, QUAD9 d=9. The reference element Λ for the class of quadrilaterals (6)Λ=ξ=ξ,η:ξ,η∈-1,1.The nodes of a quadrilateral spectral element are the Gauss-Lobatto-Legendre (GLL) nodes. Let N be an order of the spectral element; then the number of nodes in the element NφQUAD=N+12. In one-dimensional case the GLL-nodes are calculated as roots ξa of the derivative of the Legendre polynomial PN of an order N, which can be defined as PNξ=1/2NN!(dN/dξN)ξ2-1N. 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Ωe are built based on the direct product of one-dimensional Lagrange polynomials: φij,NQUADξ=liNξljNη is the ijth shape function of the order N. The one-dimensional Lagrange polynomials laN of the order N are defined as follows [3]: (7)laNξ=∏b=0d≠aNξ-ξbξa-ξb.Each polynomial laN is equal to 1 at its node a and is equal to 0 at the remaining nodes of the element Ωe:laNξb=δab, where δab is the Kronecker symbol.The approximation of solution of the problemue=u∣Ωe on a quadrilateral spectral element Ωe will be as follows m=iN+j:(8)ueξ≈∑m=0NφQUADφmQUADξueξm=∑i,j=0Nφij,NQUADξueξij=∑i,j=0Nφij,NQUADξUije.In order to integrate an arbitrary function hx over the element Ωe, the Gauss-Lobatto-Legendre (GLL) quadrature formula is used:(9)∫Ωehxdx=∫ΛheξJeξdξ=∑l=1N+12wlheξlJeξl,where wl are the GLL-weights and ξl are the coordinates of GLL-nodes.GLL-weights are calculated as follows:(10)wa=2NN+1PNξa2,ξa≠±1,wa=2NN+1,ξa=±1.The coordinates of GLL-nodes and GLL-weights (indexes i,j vary from 0 to N) can be written as follows:(11)ξl=ξij=ξi,ηj,wl=wij=wi·wj,l=iN+j.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 Ωe 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 triangles(12)Λ=ξ=ξ,η:ξ,η∈0,1,ξ+η≤1.The nodes of a triangular spectral element are the points ξj, which are obtained from the solution of the electrostatic problem described by Hesthaven and Teng [17] and Taylor et al. [18]. Let N be an order of the spectral element; then the number of nodes in the element NφTRI=(N+1)(N+2)/2.Shape functions for the elementΩe are built based on an orthogonal basis (N is the order of a spectral element):(13)ψabTRI=Pa0,02ξ1-η-11-ηaPb2a+1,02η-1,0≤a+b≤N,where Pγα,βx are the Jacobi polynomials of the order γ on the interval [-1,1], orthogonal with the weight function 1-xα1+xβ. The Vandermonde matrix elements are calculated in the nodes ξj of the element Vm,jTRI=ψmTRIξj,1≤m,j≤NφTRI. Then the shape functions will be as follows:(14)φmTRIξ=∑j=1NφTRIVm,jTRI-1ψjTRIξ.The approximation of solution of the problem ue=u∣Ωe on a triangular spectral element Ωe will be as follows:(15)ueξ≈∑m=0NφTRIφmTRIξueξm=∑m=0NφTRIφmTRIξUme.To integrate an arbitrary function hx over the element Ωe, the symmetric quadrature formulas, described by Zhang et al. [19], are used:(16)∫Ωehxdx=∫ΛheξJeξdξ=∑l=1NmaxwlheξlJeξl,where wl are positive weights, ξl are the coordinates of points, and Nmax is the number of points in the quadrature formula corresponding to the spectral element of the order N.


### 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 hexahedrons(17)Λ=ξ=ξ,η,θ:ξ,η,θ∈-1,1.The nodes of quadrilateral spectral element are the GLL-nodes. Let N be an order of the spectral element; then the number of nodes in the element NφHEX=N+13.The shape functions for the elementΩe are built based on the direct product of one-dimensional Lagrange polynomials: φijk,NHEXξ=liNξljNηlkNθ, ijkth shape function of the order N.The approximation of solution of the problemue=u∣Ωe on hexahedral spectral element Ωe will be as follows m=iN2+jN+k:(18)ueξ≈∑m=0NφHEXφm,NHEXξueξm=∑i,j,k=0Nφijk,NHEXξueξijk=∑i,j,k=0Nφijk,NHEXξUijke.To integrate an arbitrary function hx over the element Ωe, the Gauss-Lobatto-Legendre (GLL) quadrature formula is used:(19)∫Ωehxdx=∫ΛheξJeξdξ=∑l=1N+13wlheξlJeξl,where wl are the GLL-weights and ξl are the coordinates of GLL-nodes.Coordinates of the GLL-nodes and GLL-weights (indexesi,j,k vary from 0 to N)(20)ξl=ξijk=ξi,ηj,θk,wl=wijk=wi·wj·wk,l=iN2+jN+k.


### 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 tetrahedrons(21)Λ=ξ=ξ,η,θ:ξ,η,θ∈0,1,ξ+η+θ≤1.The nodes of a tetrahedron spectral element are the points ξj, which are obtained from the solution of the electrostatic problem described by Hesthaven and Teng [17]. Let N be an order of the spectral element; then the number of nodes in the element NφTET=(N+1)(N+2)(N+3)/6.The shape functions for the elementΩe are built based on the following orthogonal basis (N is the order of a spectral element):(22)ψabcTET=Pa0,02ξ1-η-θ-11-η-θaPb2a+1,02η1-θ-11-θbPc2a+2b+2,02θ-10≤a+b+c≤N,where Pγα,βx are the Jacobi polynomials of the order γ on the interval [-1,1], orthogonal with the weight function 1-xα1+xβ. The Vandermonde matrix elements are calculated in the nodes ξj of the element Vm,jTET=ψmTETξj, 1≤m,j≤NφTET. Then the shape functions will be as follows:(23)φmTETξ=∑j=1NφTETVm,jTET-1ψjTETξ.The approximation of solution of the problem ue=u∣Ωe on a tetrahedral spectral element Ωe will be as follows:(24)ueξ≈∑m=0NφTETφmTETξueξm=∑m=0NφTETφmTETξUme,In order to integrate an arbitrary function hx over the element Ωe, the symmetric quadrature formulas, described by Zhang et al. [19], are used:(25)∫Ωehxdx=∫ΛheξJeξdξ=∑l=1NmaxwlheξlJeξl,where wl are positive weights, ξl are coordinates of points, and Nmax is the number of points in the quadrature formula corresponding to the spectral element of the order N.


### 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 pyramids(26)Λ=ξ=ξ,η,θ:ξ,η,θ∈-1,1,2ξ+θ≤1,2η+θ≤1.The nodes of a pyramidal spectral element are the points ξj, 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 N be an order of the spectral element; then the number of nodes in the element(27)NφPYR=N+1N+22N+36.The shape functions for the element Ωe are built based on the following orthogonal basis (N is the order of a spectral element):(28)ψabcPYR=Pa0,02ξ1-θPb0,02η1-θ1-θ2bPc2max⁡a,b⁡+2,0θ,0≤c≤N,0≤a+b≤N-c,where Pγα,βx are the Jacobi polynomials of the order γ on the interval [-1,1], orthogonal with the weight function 1-xα1+xβ. The Vandermonde matrix elements are calculated in the nodes ξj of the element Vm,jPYR=ψmPYRξj1≤m,j≤NφPYR. Then the shape functions will be as follows:(29)φmPYRξ=∑j=1NφPYRVm,jPYR-1ψjPYRξ.The approximation of solution of the problem ue=u∣Ωe on pyramidal spectral element Ωe will be as follows:(30)ueξ≈∑m=0NφPYRφmPYRξueξm=∑m=0NφPYRφmPYRξUme.To integrate an arbitrary function hx over the element Ωe, the symmetric conical quadrature formulas, described by Felippa [20], are used:(31)∫Ωehxdx=∫ΛheξJeξdξ=∑l=1NmaxwlheξlJeξl,where wl are positive weights, ξl are coordinates of points, and Nmax is the number of points in the quadrature formula corresponding to the spectral element of the order N.


### 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(32)Λ=ξ=ξ,η,θ:θ∈-1,1,ξ,η∈0,1,ξ+η≤1.The nodes of a prismatic spectral element are the pointsξj arising from the direct product of the nodes for the triangular spectral element and one-dimensional GLL-points. Let N be an order of the spectral element; then the number of nodes in the element NφWEDGE=NφTRI(N+1).The shape functions for the elementΩe are built based on the product of the shape functions for the triangular spectral element and the Lagrange polynomials:(33)φikWEDGEξ=φiTRIξ,ηlkNθ.The approximation of solution of the problem ue=u∣Ωe on a triangular spectral element Ωe will be as follows m=kNφTRI+i:(34)ueξ≈∑m=0NφWEDGEφmWEDGEξueξm=∑m=0NφWEDGEφmWEDGEξUme.To integrate an arbitrary function hx over the element Ωe, the quadrature formulas are used, which are the unions of quadrature formulas for triangular spectral elements and the Gauss-Lobatto-Legendre quadrature formula:(35)∫Ωehxdx=∫ΛheξJeξdξ=∑l=1NmaxwlheξlJeξl,where wl are weights, ξl are coordinates of points, and Nmax is the number of points in the quadrature formula corresponding to the spectral element of the order N.


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


## 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, QUAD4d=4; eight-node, QUAD8 d=8; nine-node, QUAD9 d=9. The reference element Λ for the class of quadrilaterals (6)Λ=ξ=ξ,η:ξ,η∈-1,1.The nodes of a quadrilateral spectral element are the Gauss-Lobatto-Legendre (GLL) nodes. Let N be an order of the spectral element; then the number of nodes in the element NφQUAD=N+12. In one-dimensional case the GLL-nodes are calculated as roots ξa of the derivative of the Legendre polynomial PN of an order N, which can be defined as PNξ=1/2NN!(dN/dξN)ξ2-1N. 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Ωe are built based on the direct product of one-dimensional Lagrange polynomials: φij,NQUADξ=liNξljNη is the ijth shape function of the order N. The one-dimensional Lagrange polynomials laN of the order N are defined as follows [3]: (7)laNξ=∏b=0d≠aNξ-ξbξa-ξb.Each polynomial laN is equal to 1 at its node a and is equal to 0 at the remaining nodes of the element Ωe:laNξb=δab, where δab is the Kronecker symbol.The approximation of solution of the problemue=u∣Ωe on a quadrilateral spectral element Ωe will be as follows m=iN+j:(8)ueξ≈∑m=0NφQUADφmQUADξueξm=∑i,j=0Nφij,NQUADξueξij=∑i,j=0Nφij,NQUADξUije.In order to integrate an arbitrary function hx over the element Ωe, the Gauss-Lobatto-Legendre (GLL) quadrature formula is used:(9)∫Ωehxdx=∫ΛheξJeξdξ=∑l=1N+12wlheξlJeξl,where wl are the GLL-weights and ξl are the coordinates of GLL-nodes.GLL-weights are calculated as follows:(10)wa=2NN+1PNξa2,ξa≠±1,wa=2NN+1,ξa=±1.The coordinates of GLL-nodes and GLL-weights (indexes i,j vary from 0 to N) can be written as follows:(11)ξl=ξij=ξi,ηj,wl=wij=wi·wj,l=iN+j.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 Ωe 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 triangles(12)Λ=ξ=ξ,η:ξ,η∈0,1,ξ+η≤1.The nodes of a triangular spectral element are the points ξj, which are obtained from the solution of the electrostatic problem described by Hesthaven and Teng [17] and Taylor et al. [18]. Let N be an order of the spectral element; then the number of nodes in the element NφTRI=(N+1)(N+2)/2.Shape functions for the elementΩe are built based on an orthogonal basis (N is the order of a spectral element):(13)ψabTRI=Pa0,02ξ1-η-11-ηaPb2a+1,02η-1,0≤a+b≤N,where Pγα,βx are the Jacobi polynomials of the order γ on the interval [-1,1], orthogonal with the weight function 1-xα1+xβ. The Vandermonde matrix elements are calculated in the nodes ξj of the element Vm,jTRI=ψmTRIξj,1≤m,j≤NφTRI. Then the shape functions will be as follows:(14)φmTRIξ=∑j=1NφTRIVm,jTRI-1ψjTRIξ.The approximation of solution of the problem ue=u∣Ωe on a triangular spectral element Ωe will be as follows:(15)ueξ≈∑m=0NφTRIφmTRIξueξm=∑m=0NφTRIφmTRIξUme.To integrate an arbitrary function hx over the element Ωe, the symmetric quadrature formulas, described by Zhang et al. [19], are used:(16)∫Ωehxdx=∫ΛheξJeξdξ=∑l=1NmaxwlheξlJeξl,where wl are positive weights, ξl are the coordinates of points, and Nmax is the number of points in the quadrature formula corresponding to the spectral element of the order N.


## 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 hexahedrons(17)Λ=ξ=ξ,η,θ:ξ,η,θ∈-1,1.The nodes of quadrilateral spectral element are the GLL-nodes. Let N be an order of the spectral element; then the number of nodes in the element NφHEX=N+13.The shape functions for the elementΩe are built based on the direct product of one-dimensional Lagrange polynomials: φijk,NHEXξ=liNξljNηlkNθ, ijkth shape function of the order N.The approximation of solution of the problemue=u∣Ωe on hexahedral spectral element Ωe will be as follows m=iN2+jN+k:(18)ueξ≈∑m=0NφHEXφm,NHEXξueξm=∑i,j,k=0Nφijk,NHEXξueξijk=∑i,j,k=0Nφijk,NHEXξUijke.To integrate an arbitrary function hx over the element Ωe, the Gauss-Lobatto-Legendre (GLL) quadrature formula is used:(19)∫Ωehxdx=∫ΛheξJeξdξ=∑l=1N+13wlheξlJeξl,where wl are the GLL-weights and ξl are the coordinates of GLL-nodes.Coordinates of the GLL-nodes and GLL-weights (indexesi,j,k vary from 0 to N)(20)ξl=ξijk=ξi,ηj,θk,wl=wijk=wi·wj·wk,l=iN2+jN+k.


## 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 tetrahedrons(21)Λ=ξ=ξ,η,θ:ξ,η,θ∈0,1,ξ+η+θ≤1.The nodes of a tetrahedron spectral element are the points ξj, which are obtained from the solution of the electrostatic problem described by Hesthaven and Teng [17]. Let N be an order of the spectral element; then the number of nodes in the element NφTET=(N+1)(N+2)(N+3)/6.The shape functions for the elementΩe are built based on the following orthogonal basis (N is the order of a spectral element):(22)ψabcTET=Pa0,02ξ1-η-θ-11-η-θaPb2a+1,02η1-θ-11-θbPc2a+2b+2,02θ-10≤a+b+c≤N,where Pγα,βx are the Jacobi polynomials of the order γ on the interval [-1,1], orthogonal with the weight function 1-xα1+xβ. The Vandermonde matrix elements are calculated in the nodes ξj of the element Vm,jTET=ψmTETξj, 1≤m,j≤NφTET. Then the shape functions will be as follows:(23)φmTETξ=∑j=1NφTETVm,jTET-1ψjTETξ.The approximation of solution of the problem ue=u∣Ωe on a tetrahedral spectral element Ωe will be as follows:(24)ueξ≈∑m=0NφTETφmTETξueξm=∑m=0NφTETφmTETξUme,In order to integrate an arbitrary function hx over the element Ωe, the symmetric quadrature formulas, described by Zhang et al. [19], are used:(25)∫Ωehxdx=∫ΛheξJeξdξ=∑l=1NmaxwlheξlJeξl,where wl are positive weights, ξl are coordinates of points, and Nmax is the number of points in the quadrature formula corresponding to the spectral element of the order N.


## 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 pyramids(26)Λ=ξ=ξ,η,θ:ξ,η,θ∈-1,1,2ξ+θ≤1,2η+θ≤1.The nodes of a pyramidal spectral element are the points ξj, 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 N be an order of the spectral element; then the number of nodes in the element(27)NφPYR=N+1N+22N+36.The shape functions for the element Ωe are built based on the following orthogonal basis (N is the order of a spectral element):(28)ψabcPYR=Pa0,02ξ1-θPb0,02η1-θ1-θ2bPc2max⁡a,b⁡+2,0θ,0≤c≤N,0≤a+b≤N-c,where Pγα,βx are the Jacobi polynomials of the order γ on the interval [-1,1], orthogonal with the weight function 1-xα1+xβ. The Vandermonde matrix elements are calculated in the nodes ξj of the element Vm,jPYR=ψmPYRξj1≤m,j≤NφPYR. Then the shape functions will be as follows:(29)φmPYRξ=∑j=1NφPYRVm,jPYR-1ψjPYRξ.The approximation of solution of the problem ue=u∣Ωe on pyramidal spectral element Ωe will be as follows:(30)ueξ≈∑m=0NφPYRφmPYRξueξm=∑m=0NφPYRφmPYRξUme.To integrate an arbitrary function hx over the element Ωe, the symmetric conical quadrature formulas, described by Felippa [20], are used:(31)∫Ωehxdx=∫ΛheξJeξdξ=∑l=1NmaxwlheξlJeξl,where wl are positive weights, ξl are coordinates of points, and Nmax is the number of points in the quadrature formula corresponding to the spectral element of the order N.


## 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(32)Λ=ξ=ξ,η,θ:θ∈-1,1,ξ,η∈0,1,ξ+η≤1.The nodes of a prismatic spectral element are the pointsξj arising from the direct product of the nodes for the triangular spectral element and one-dimensional GLL-points. Let N be an order of the spectral element; then the number of nodes in the element NφWEDGE=NφTRI(N+1).The shape functions for the elementΩe are built based on the product of the shape functions for the triangular spectral element and the Lagrange polynomials:(33)φikWEDGEξ=φiTRIξ,ηlkNθ.The approximation of solution of the problem ue=u∣Ωe on a triangular spectral element Ωe will be as follows m=kNφTRI+i:(34)ueξ≈∑m=0NφWEDGEφmWEDGEξueξm=∑m=0NφWEDGEφmWEDGEξUme.To integrate an arbitrary function hx over the element Ωe, the quadrature formulas are used, which are the unions of quadrature formulas for triangular spectral elements and the Gauss-Lobatto-Legendre quadrature formula:(35)∫Ωehxdx=∫ΛheξJeξdξ=∑l=1NmaxwlheξlJeξl,where wl are weights, ξl are coordinates of points, and Nmax is the number of points in the quadrature formula corresponding to the spectral element of the order N.


## 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 Ux, Uy and the relative errors εxSEM, εySEM 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; t denotes time, and N 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.Figure 1


### 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 Figure2 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.Figure 2


## 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 Ux, Uy and the relative errors εxSEM, εySEM 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; t denotes time, and N 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.Figure 1


## 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 Figure2 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.Figure 2


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



---

*Source: 1797561-2017-12-06.xml*