Spectral Element Method for Transient Wave Equation

Authors: Ernesto Jose Quiroz Vazquez Juan Angelo Vargas Fajardo

Supervisors: Jonas Hillenbrand

1. Introduction

In this Wiki, the supporting theory and a couple of examples for the implementation of Spectral Elements are presented. In the context of this project, spectral elements are finite elements of polynomial degree greater than one, that have the property of yielding a diagonal mass matrix. The high polynomial degree has the advantage of describing curvature better than linear elements, reducing the spatial discretization error when using them for Finite Element Analysis. The diagonal matrix has the advantage of being easily invertible, which reduces the computational effort to solve the linear system of equations that results from FEM.

In this document, the method to generate spectral elements is discussed. Afterwards, some examples that showcase the properties of these elements are presented, and the results are briefly discussed. In these examples, the spectral elements are used for spatial discretization while the Newmark method is chosen for the time discretization, which allows for transient wave analysis. First, a 1D Truss element (that can be placed in 2D space) example, and later, a 2D plane stress example of a beam.

Example simulation is done using spectral elements.

2. Spectral Element Method

2.1. Transient Wave Equation

This project's main topic is the solution through finite elements of the transient wave equation through spectral elements. Initially, it is defined for 1D and later expanded to 2D along the equations of motion.

2.1.1. 1D Formulation

The 1D wave equation reads:

$\begin{array}{l}\displaystyle \frac{\partial^2 u\left(x,t\right)}{\partial t^2}-c^2\frac{\partial^2 u\left(x,t\right)}{\partial x^2}=0\end{array}$

$\begin{array}{l}\displaystyle x \in \Omega \subseteq \mathbb{R}\end{array}$

$\begin{array}{l}\displaystyle t \in ]0,T[\end{array}$

where $\begin{array}{l}u(x,t)\end{array}$ is the solution of the previous equation with some boundary conditions.

This can be solved in a weak form as follows:

$\begin{array}{l}\displaystyle \frac{d^2}{dt^2}{\int_{\Omega}u v dx}+c^2\int_{\Omega}{\frac{\partial u}{\partial x}\frac{\partial v}{\partial x}}=0\end{array}$

Discretizing our spatial domain $\begin{array}{l}\Omega\end{array}$ in $\begin{array}{l}k\end{array}$ elements, each element having $\begin{array}{l}n\end{array}$ nodes, our solution at each element $\begin{array}{l}u^i\end{array}$ can be approximated by:

$\begin{array}{l}\displaystyle u^{(i)}(x,t)\approx\sum_{j=1}^{n}N_j(x) u_j^(i)(t)=[N]^T[u]\end{array}$

Where

$\begin{array}{l}\displaystyle [N]=[N1\quad N2\quad ... \quad N_n]^T\end{array}$

$\begin{array}{l}\displaystyle [u]=[u1\quad u2\quad ... \quad u_n]^T\end{array}$

$\begin{array}{l}\displaystyle \left[\frac{dN}{dx}\right]=\left[\frac{dN_1}{dx} \quad \frac{dN_2}{dx}\quad ... \quad \frac{dN_m}{dx}\right]^T\end{array}$

Using the Bubnov-Galerkin approach, we calculate the local stiffness matrix $\begin{array}{l}K_{ij}\end{array}$ and local mass matrix $\begin{array}{l}M_{ij}\end{array}$ as follows:

$\begin{array}{l}\displaystyle K_{ij}=\int_{\Omega^e}{\left[\frac{dN}{dx}\right]\left[\frac{dN}{dx}\right]^T dx}\end{array}$

$\begin{array}{l}\displaystyle M_{ij}=\int_{\Omega^e}{\left[N\right]\left[N\right]^T dx}\end{array}$

With this result, we can assemble our global stiffness matrix and global mass matrix and get the following equation, where we spatial discretized our equation:

$\begin{array}{l}\displaystyle M\ddot{u}+Ku=0\end{array}$

Which can be solve with a temporal scheme, such as the Newmark scheme.

2.1.2. 2D Formulation

Similarly, a dynamic 2D Linear Elastic Problem, assuming Plane Stress can be also spatial discretize as:

$\begin{array}{l}\displaystyle M\ddot{u}+Ku=F\end{array}$

Being the same result as a 1D wave equation. The local stiffness matrix $\begin{array}{l}K_{ij}\end{array}$ and local mass matrix $\begin{array}{l}M_{ij}\end{array}$ can be written as follows:

$\begin{array}{l}\displaystyle K_{ij}=\int_{\Omega^e}{\left(\mathbf{L}\left[N\right]\right)^T{C}\left(\mathbf{L}\left[N\right]\right)\ d\Omega^e}\end{array}$

$\begin{array}{l}\displaystyle M_{ij}=\mu\int_{\Omega^e}{\left[N\right]^{T}\left[N\right]d\Omega^e}\end{array}$

Where $\begin{array}{l}\mathbf{L}\end{array}$ is the linear operator for plane stress, $\begin{array}{l}C\end{array}$ , the material matrix and $\begin{array}{l}\mu\end{array}$ is the mass per unit area.

2.2. Spectral Elements

Spectral Elements is the name given to such elements that allow for lumping in the Mass Matrix obtained from the space discretization of a physical domain through the Finite Elements Method (FEM). A lumped mass matrix has the advantage of being diagonal and, therefore, simple to invert, reducing the time necessary for solving the system of equations. Spectral elements also have the property of using polynomial shape functions with degrees higher than 1.

To generate spectral elements with these properties, we need special shape functions that have the following characteristics:

Be high-order polynomial shape functions.
Yield a diagonal mass matrix (Orthogonality between shape functions).
Be valid for interpolation of the state variables.
Have a corresponding efficient and accurate integration rule.

The shape functions that fulfill these desired properties are the Gauss Lobatto shape functions described in the following section.

2.3. Gauss-Lobatto Shape Functions

The Gauss-Lobatto-Legendre quadrature, including the determination of the associated Gauss-Lobatto nodes (quadrature points) and shape functions, is an extension of the Gauss-Legendre quadrature which includes the endpoints of the integration interval as quadrature points.

These shape functions have the desired characteristics to form spectral elements, meaning they are orthogonal to each other, allowing for a diagonal mass matrix. Additionally, they have a corresponding quadrature similar to the Gauss-Legendre Quadrature, but with adjusted node coordinates and weights. Finally, the Gauss Lobatto shape functions are generated as Lagrange polynomials with not equidistant coordinate points. These concepts are further elaborated in the subsections below.

2.3.1. Orthogonality

In this context, the inner product of two functions is defined as the integral over the domain of the product of these functions.

$\begin{array}{l}\displaystyle < u,v>:=\int_{\Omega}{u\cdot v\ d\Omega}\end{array}$

we note that the inner product is always positive definite

$\begin{array}{l}< u,u>>0\leftrightarrow u\neq0 \\\end{array}$

, and that two functions are orthogonal with respect to an inner product if and only if their inner product is zero.

$\begin{array}{l}{Orthogonality}\leftrightarrow\ \ < u,v>=0 \\\end{array}$

Recalling the mass matrix from the spatially discretized equations of motion,

$\begin{array}{l}\mathbf{M_e}=\mu\int_{\Omega^e}\mathbf{N^TN}d\Omega^e \\\end{array}$

each of its entries can be written as the inner product of two shape functions, scaled my $\begin{array}{l}\mu\end{array}$

$\begin{array}{l}M_{ij}=\mu\int_{\mathrm{\Omega}^e}{N_i\left(x\right)N_j(x)dx}=\mu< N_i,N_j> \\\end{array}$

One of our desired outcomes is to have a diagonal mass matrix M, which means that

$\begin{array}{l}\displaystyle M_{ij} = 0 \ \ i\neq j\end{array}$

recognizing that $\begin{array}{l}\mu\end{array}$ is not zero, then the inner product of the shape functions must be zero

$\begin{array}{l}\displaystyle \mu< N_i,N_j>=0,\ \mu\neq0\end{array}$

therefore, we are looking for shape functions that are orthogonal with respect to this inner product

$\begin{array}{l}\displaystyle < N_i,N_j>=0,\ i\neq0 j\end{array}$

The Gauss-Lobatto shape functions fulfill this orthogonality condition, making them appropriate for use along with the Spectral Elements.

2.3.2. Gauss-Lobatto Quadrature

The general numerical approximation of the integral of a function reads

$\begin{array}{l}\displaystyle \int_{-1}^{1}f\left(\xi\right)d\xi\approx\sum_{i=1}^{n}{w_if(\xi_i)}\end{array}$

applied to the discretization of the mass matrix, we obtain

$\begin{array}{l}\displaystyle \int_{-1}^{1}{N_i\left(\xi\right)N_j(\xi)d\xi}\approx\sum_{i=1}^{n} {w_iN_i\left(\xi_i\right)N_j(\xi_i)}\end{array}$

for which

$\begin{array}{l}\xi_i=nodes\ Coordinates \\ w_i=Quadrature\ Weights\end{array}$

and for two dimensions, the quadrature rule becomes

$\begin{array}{l}\displaystyle \int_{-1}^{1}{\int_{-1}^{1}f\left(\xi,\eta\right)d\xi d\eta}=\sum_{j=1}^{n}\sum_{i=1}^{n}{w_iw_jf(\xi_i,\eta_j)}\end{array}$

Until here, there is not much difference compared to the frequently used Gauss-Legendre quadrature; through the definition of the polynomials and its respective roots, the differences begin to appear.

The necessary coordinates and weights for selected polynomial degrees can be found in tables published in various literature on the topic. For the implementation in the BM-MFEM software, an arbitrary polynomial degree can be chosen. Therefore, the way to find the coordinates and weights was programmed. Without going into detail, the steps taken are outlined here.

The Legendre Polynomials of the same order are found. These polynomials are recursive and are, in a general manner described as $\begin{array}{l}P_{k+1}\left(x\right)=\alpha\left(k\right)P_k(x)-\beta\left(k\right)P_{k-1}(x)\end{array}$
The derivative of the Legendre $\begin{array}{l}P^\prime(x)\end{array}$ polynomials is obtained.
The roots of this derivative are the roots for the Gauss Lobatto quadrature $\begin{array}{l}P^\prime(x) \longrightarrow \xi_i\end{array}$ . A Newton-Raphson algorithm was employed to find the roots, using the Chebyshev polynomials as the initial guess for the roots to avoid converging to an undesired root.
Finally, the weights can be deduced from the coordinates using the following equation $\begin{array}{l}w_i = \frac{2}{n(n-1)[P_{n-1}(x_i)]^2}, \space x_i \neq \pm1\end{array}$ .

2.3.3. Lagrange Polynomials

The interpolatory polynomials used for this implementation of the spectral elements are Lagrange Polynomials. The Basis polynomials used to construct the Lagrange polynomials use the position $\begin{array}{l}x\end{array}$ of the roots from the Gauss-Lobatto quadrature, shown in the previous section.

In general, a Lagrange polynomial $\begin{array}{l}P(x)\end{array}$ is a polynomial of degree $\begin{array}{l}n\end{array}$ that passes through $\begin{array}{l}n+1\end{array}$ given data points $\begin{array}{l}(x_0 ,y_0 ), (x_1,y_1),…,(x_n,y_n)\end{array}$ . It is defined by the formula:

$\begin{array}{l}\displaystyle P(x) = L_0(x)y_0 + L_1(x)y_1 + \ldots + L_n(x)y_n\end{array}$

where:

$\begin{array}{l}L_i(x)\end{array}$ is the i-th Lagrange basis polynomial,
$\begin{array}{l}y_i\end{array}$ is the corresponding function value at $\begin{array}{l}x_i\end{array}$ .

The Lagrange Basis Polynomial $\begin{array}{l}L_i(x)\end{array}$ for a given data point $\begin{array}{l}(x_i,y_i)\end{array}$ is given by the formula:

$\begin{array}{l}\displaystyle L_i(x) = \prod_{j=0, j\neq i}^{n} \frac{x - x_j}{x_i - x_j}\end{array}$

Using the Gauss-Lobatto roots instead of equidistant points has a visible impact on the resulting shape functions, which can be visualized in the following image for a polynomial of degree three:

Lagrange Polynomials using equidistant points (blue) and Gauss Lobatto roots (red).

It is observed that the Lagrange polynomials have now a maximum value larger than 1, and the zero crossings do not coincide. A comparison of the root points can be seen in the following image:

Equidistant points in local coordinates compared to Gauss-Lobatto quadrature roots.

2.3.4. Effect of Gauss Lobatto Shape Functions on Mass Matrix

The properties of the Gauss-Lobatto Shape functions, which include orthogonality, have an effect on the distribution of values on the Mass Matrix. This is observed as a diagonalization of the Mass Matrix; through the use of Gauss-Lobatto shape functions, a "natural" mass lumping is achieved, therefore we can say that our continuous system has been discretized in such a way that no mass coupling occurs.

In the following image, we can observe the effects of using Gauss-Lobatto shape functions for the discretization of a given continuous system, on the resulting Mass Matrix:


Sparcity pattern of the Mass Matrix for a given system using (blue) Lagrange Polynomials with equidistant points vs Lagrange Polynomials with Gauss-Lobatto points (red).

2.4. Newmark Scheme

The numerical time step procedure used along with the spectral elements for the simulation of dynamic systems in this project is the Newmark scheme. This scheme was already implemented into the BM-MFEM software and only a brief description of it is given in this report.

The Newmark Scheme is a numerical method used for the time-step discretization of dynamic systems, often used in the context of structural dynamics and finite element analysis. It is an implicit integration scheme that can be interpreted as a prediction-correction method.

The method is based on the following approximations for position, velocity, and acceleration

$\begin{array}{l}\displaystyle u_{k+1} = u_k + \Delta t \cdot \dot{u}_k + \frac{\Delta t^2}{2}[(1-2\beta)\ddot{u}_k + 2\beta\ddot{u}_{k+1}]\end{array}$

$\begin{array}{l}\displaystyle \dot{u}_{k+1} = \dot{u}_k + \Delta t[(1-\gamma)\ddot{u}_k + \gamma\ddot{u}_{k+1}]\end{array}$

$\begin{array}{l}\displaystyle \ddot{u}_{k+1} = \frac{1}{\beta\Delta t^2}(u_{k+1} - u_k) - \frac{1}{\beta\Delta t}\dot{u}_k - \frac{(1-2\beta)}{2\beta}\ddot{u}_k\end{array}$

where $\begin{array}{l}\Delta t\end{array}$ represents the chosen time step, and the parameters $\begin{array}{l}\beta\end{array}$ and $\begin{array}{l}\gamma\end{array}$ are introduced to control the numerical integration (1/4 and 1/2 respectively corresponding to trapezoidal rule).

Regarding stability and accuracy, the Newmark scheme is conditionally stable, and its accuracy depends on the choices of $\begin{array}{l}\beta\end{array}$ and $\begin{array}{l}\gamma\end{array}$ . The standard choice of $\begin{array}{l}\beta = \frac{1}{4}\end{array}$ and $\begin{array}{l}\gamma = \frac{1}{2}\end{array}$ imply unconditional stability and second-order accuracy.

It is noteworthy to mention that during this time discretization procedure, the expression

$\begin{array}{l}\displaystyle \mathbf{M}{\ddot{u}}_{k+1}+\mathbf{C}\left({\dot{\widetilde{u}}}_{k+1}+{\ddot{u}}_{k+\mathbf{1}}\gamma\Delta t\right)+\mathbf{K}({\widetilde{u}}_{k+\mathbf{1}}+{\ddot{u}}_{k+\mathbf{1}}\beta\mathbf{\Delta}t^2)=f(t_{k+1})\end{array}$

is used to approximate the equation of motion in the next time step. From this expression, the procedure solves for the acceleration $\begin{array}{l}\ddot u_{k+1}\end{array}$ . The right hand contains a term which requires the inversion of the Mass matrix.

$\begin{array}{l}\displaystyle {\ddot{u}}_{k+1}=\left(\mathbf{M}+\mathbf{C}\gamma\mathbf{\Delta}t+\mathbf{K}\beta\mathbf{\Delta}t^\mathbf{2}\right)^{-1}\left(f\left(t_{k+\mathbf{1}}\right)-{\mathbf{C}\dot{\widetilde{u}}}_{k+\mathbf{1}}-\mathbf{K}{\widetilde{u}}_{k+\mathbf{1}}\right)\end{array}$

Remembering that the choice of doing the spatial discretization using Gauss-Lobatto shape functions had as a consequence the diagonalization of the Mass matrix, which is generally computationally easier to invert, it can be speculated that the computational effort for each iteration of the Newmark Scheme might be reduced. The computational time is further explored in the following sections in the context of some examples.

3. 1D Truss Spectral Elements in 2D Space examples.

In the context of this project, the spectral elements were implemented inside the framework of the BM-MFEM software. This section presents and briefly discusses two examples of simulations using one-dimensional truss elements defined in 2D space. The capabilities of the implemented code are summarized in the following image:

Code implemented within the framework of BM-MFEM in MATLAB.
Spectral Element versions for Trusses in 2D space and Plane Stress Elements in 2D.
Shape function options using Equidistant or Gauss-Lobatto node coordinates.
Definition of material properties (Poisson Ratio, Youngs Modulus, density, etc.).
Option to refine the mesh in terms of elements (h-refinement) or polynomial degree (p-refinement). p-refinement is implemented for an arbitrary polynomial degree.
Implementation of the spectral elements has been tested for the static "Simple" solver and the dynamic time step procedure "Newmark Scheme".

3.1. Linear Initial Displacement

The first example is a truss element supported by pins on both sides, subject to initial displacement conditions under no initial velocity. The initial displacement has the shape of a triangular function, or triangle wave, which is described by the function

$\begin{array}{l}\displaystyle f(t) = \frac{y_0}{T} \left( T - \left| t - \frac{T}{2} \right| \right)\end{array}$

where $\begin{array}{l}y_0\end{array}$ is the initial amplitude and $\begin{array}{l}T\end{array}$ is the period of the wave.

The following model can be interchangeably interpreted for truss elements, strings with a small sag and torsional rods, by interpreting the resulting displacement field as axial displacement, perpendicular displacement or rotation respectively.

	Poisson Ratio: $\begin{array}{l}\nu = 0.3\end{array}$ Density: $\begin{array}{l}\rho = 1000 \frac{kg}{m^3}\end{array}$ Young Modulus: $\begin{array}{l}E = 4000 MPa\end{array}$ Max initial displacement: $\begin{array}{l}y_0 = 1m\end{array}$ Cross Section: $\begin{array}{l}A = 0.1 m^2\end{array}$ Length: $\begin{array}{l}l = 2m\end{array}$ Time step: $\begin{array}{l}\Delta t =0.1 s\end{array}$
String fixed on both sides with a prescribed initial displacement

This test case is chosen because there is a known analytical solution, coming from the superposition of opposing traveling waves, of half the amplitude of the initial displacement

$\begin{array}{l}\displaystyle y(x,t) = \frac{y_0(x-ct)+y_0(x+ct)}{2}\end{array}$

where $\begin{array}{l}y_0\end{array}$ is the description of the initial displacement and $\begin{array}{l}c\end{array}$ is the speed of the travelling wave. This analytical solution is discussed in depth in the Structural Dynamics script by Prof. Dr-Ing. Gerhard Müller.

The solution can be observed in the following figure for two different spatial discretization and polynomial degree:


p=1, 15 elements, 16 nodes (linear elements)	p=5, 3 elements, 16 nodes (spectral elements)

The graphics above compare the same test case when using linear elements against spectral elements of polynomial degree 5. Both scenarios contain the same amount of degrees of freedom, for a fair comparison. It is observed that the spectral elements can describe the displacement field with a smaller deviation from the analytical solution. In the following image, this error is quantified as the root mean square error of the whole displacement field, as it evolves with each time step to have a more quantifiable difference. In this graph, we have 50 linear elements, compared to just 5 spectral elements of polynomial degree 10, again maintaining the same nodes (analogously, DOFs).


a) Root Mean Square Error evolving with each time step (iterations) for the same amount of nodes.	b) RMSE for 1000 linear elements compared to 5 spectral elements

In graphic a) above, we observe an oscillating error increasing with time. This behavior is also present within the spectral elements but with a much smaller magnitude. This is showcased in graphic b), where a h-refinement is done to increase the amount of linear elements to reach a similar RMSE value compared to the spectral elements. We observe that both errors have the same behavior, and even though they are of a much lower magnitude than for graph a), their tendency is of continuous growth. It is deduced that the total error for this spatial discretization level is dominated by the time discretization error, which shows numerical dissipation evolving with time observable as a phase shift in the traveling waves (given enough time steps).

It is then possible to find, for each polynomial degree, the number of nodes (analogously DOFs) necessary to reach the point where the spatial discretization error is negligible compared to the time discretization error domination. With other words, how many nodes are needed for the spatial discretization error introduced by finite elements to become irrelevant. In the following graph, this concept is plotted for various polynomial degrees.

Nodes needed for the time discretization error to dominate

From the figure above, we observe a maximum value of approximately 100 nodes for polynomial degree 1 and a strong decrement when using $\begin{array}{l}p=2\end{array}$ and then similar values for the following polynomial degrees for this specific test case. We obtain the following graph if we plot the corresponding CPU time required for a simulation with dominating time discretization error.

CPU time required to simulate with a dominating time discretization error.

A trend comparable to that observed in the amount of nodes necessary is repeated in this graph. Again, a decrement from 41 seconds to 18 seconds is achieved by increasing the polynomial degree from 1 to 2. The speed-up flattens at the 4th polynomial degree, after which the CPU time remains in a range between 7 and 13 seconds for this example. It is natural that for a lower amount of nodes, consequently degrees of freedom, the time required to solve the corresponding system of equations is reduced. This example showcases one advantage of the spectral elements, mainly due to their possibility of being p-refined. The optimum polynomial degree will vary across applications depending on the behavior of the displacement field to be represented, which means that p-refinement is more accurate and efficient than h-refinement.

3.2. Wave Initial Displacement

The second example consists again of initial displacement conditions, in this instance with a summation of 4 sine functions with increasing frequency and diminishing amplitudes. The recursive function describing this sine sum is

$\begin{array}{l}\displaystyle f(x,n) = \frac{1}{a}sin(\frac{2a\pi (x+1)}{2}) + f(x,n-1) \ \ while \ \ n>0\end{array}$

where $\begin{array}{l}a = 2(n-1)\end{array}$ for $\begin{array}{l}n>1\end{array}$ and $\begin{array}{l}a = n\end{array}$ for $\begin{array}{l}n=1\end{array}$ , $\begin{array}{l}n \in \mathbb{N}\end{array}$

	Poisson Ratio: $\begin{array}{l}\nu = 0.3\end{array}$ Density: $\begin{array}{l}\rho = 1000 \frac{kg}{m^3}\end{array}$ Young Modulus: $\begin{array}{l}E = 4000 MPa\end{array}$ Cross Section: $\begin{array}{l}A = 0.1 m^2\end{array}$ Length: $\begin{array}{l}l = 2m\end{array}$ Time step: $\begin{array}{l}\Delta t =0.1 s\end{array}$
String fixed with a pin support on both sides with an initial displacement given by the sum of sine functions.

In the following graphs, we observe the results obtained from simulating using a polynomial degree of 2 and a polynomial degree of 4. In this example the displacement field cannot be described by 4 elements of polynomial degree 2. However, employing a higher polynomial for the same amount of elements allows the analytical solution to be closely approximated.


p = 2, 4 elements, 9 Nodes	p = 4, 4 elements, 17 Nodes

In the same manner as the analysis done for the triangular wave example, we are looking to compare the performance of the spectral elements with linear elements for similar spatial discretization errors. This error is chosen to be as small as to let it be dominated by the time discretization error coming from the Newmark Scheme time step procedure. In the following image the amount of nodes required to reach this relationship between the spatial and time error is depicted.

Nodes needed for the time discretization error to dominate the spatial error.

We observe the same pattern as in the previous example, with 500 nodes required to achieve this error with p=1, and a reduction to just 80 nodes needed with spectral elements of p=2. Further decrements can be seen, with varying amount of nodes needed with a lowest value at 29 nodes for p=14. Now, we look at the required CPU time in the following plots.


CPU time required to simulate with a dominating time discretization error.	CPU time per degree of freedom for dominating time discretization error.

Continuing the observed trend, the CPU time is reduced with increasing polynomial degree, showing a change from 250 seconds to 25 seconds corresponding to p=1 and p=2. Again, diminishing decrements are observed for higher polynomial degrees that plateau between 5 and 10 seconds. When the CPU time is considered for each degree of freedom, the reduction of time with increasing polynomial degree is still observed, but it is no longer an order of magnitudes of change. For example, the time needed for each degree of freedom with a polynomial degree of 1 is approximately 3 seconds for 10000 iterations. In comparison, the time required for a polynomial degree of 2 for the same amount of iterations is about 2 seconds per degree of freedom.

4. 2D Plane Stress Spectral Elements Examples

In the following examples, the second type of spectral elements implemented into the BM-MFEM software are presented. This elements are developed for linear, isotropic, plane stress conditions. Similar to the Truss elements, they can be refined with an arbitrary polynomial degree, and can be chosen to have either Gauss-Lobatto local coordinates, or equidistant ones. These elements are in the framework of BM-MFEM compatible with the Simple Solver and the Newmark Solver.

Depiction of functionalities implemented along with the spectral elements within the BM-MFEM framework

Two examples will be introduced in this section, the first corresponds to the modal analysis of a beam, and the second one looks at a structure subjected to an initial displacement.

4.1. Beam Modal Analysis

In this example a modal analysis of a beam with a fixed side is done using the implemented spectral elements. This configuration is frequently used as a benchmark to evaluate methods of computation and it was therefore chosen as a viable example. The purpose of this analysis is to obtain the eigenfrequencies of the structure and to visualize the corresponding eigenmodes.

	$\begin{array}{l}E = 2.10 \times 10^{11}\end{array}$ $\begin{array}{l}t=1\end{array}$ $\begin{array}{l}h=10\end{array}$ $\begin{array}{l}L = 100\end{array}$ $\begin{array}{l}\nu = 0.3\end{array}$ $\begin{array}{l}\rho = 0.3\end{array}$
Beam fixed on the left side, with length L and height h

For the modal analysis we look at the undamped equation of motion in matrix form and solve the inherent eigenvalue problem. For this example we solved the discretized equation of motion resulting from the use of the Finite Element Method with linear elements and spectral elements.

$\begin{array}{l}\displaystyle Mu ̈+Ku=0\end{array}$

$\begin{array}{l}\displaystyle (\mathbf{K} - \omega^2 \mathbf{M}) \mathbf{\phi} = \mathbf{0}\end{array}$

$\begin{array}{l}\displaystyle det(\mathbf{K} - \lambda \mathbf{M}) = 0\end{array}$

In the following graphics we can observe the obtained eigenmodes with their respective eigenfrequency. It is noteworthy that the mode shape for the 9th and 10th mode do not coincide, given that the frequencies of these two modes are close to one another, depending on the error of the calculation the 10th mode can become the 9th mode and vice-versa.

40 elements p = 1 DOF = 63
12 elements p = 2 DOF = 65

The corresponding relative error for the circular eigenfrequencies for this simulation are collected in the following table. The relative error is obtained by comparing to the analytical solution obtained from with the Bernoulli Beam Theory. We observe that the error when using a polynomial degree of 1 is generally greater that when using a polynomial degree of second order. A notable exception arises in the 3rd eigenfrequency, that corresponds with a longitudinal mode, which can be well described using linear finite elements.

	Relative Error
$\begin{array}{l}\omega\end{array}$	p=1	p=2
1	35.9%	0.5%
2	32.8%	1.6%
3	0.2%	0.2%
4	29.3%	3.1%
5	25.8%	5.2%
6	6.7%	6.9%
7	14.9%	0.9%
8	19.7%	13.0%
9	0.2%	0.2%
10	16.9%	0.1%

Further refinement using more elements with higher polynomial degree is possible with the implemented spectral elements. In the following graphic, the first 10 eigenmodes with the corresponding circular frequency and relative error are shown, for a polynomial degree of 4 and 250 elements. In the higher modes, it can be seen that the cross section does not remain perpendicular to the midplane of the beam, which is an assumption for the Euler Bernoulli Beam theory, indicating that higher modes might not be within the limits of its validity.

Eigenmodes for a beam with a fixed end, using 250 spectral elements of 4th polynomial degree.

The subsequent images show the relative error for the 1st and 10th eigenfrequency for polynomial degrees 1, 2 and 3. The error observed for p=2 is considerably lower than that of p=1, while further increase in the polynomial degree yields diminishing improvements.


a) Relative error of the 1st eigenfrequency for p=1,2,3	b) Relative error of the 10th eigenfrequency for p=1,2,3

The relative error is now compared between different shape function configurations, without varying the polynomial degree. The effect that using the Gauss-Lobatto coordinate points for the shape functions has on the accuracy of the results is evaluated. In the plot below, we observe that there is not a significant difference in relative error between the 2 approaches. From this comparison it can be stated that for this example, using Gauss-Lobatto does not improve or worsen the accuracy of the simulation. In the next example the computational time required for each of these approaches is explored.


Relative error of the 1st and 10th eigenfrequency when using Equidistant or Gauss-Lobatto node coordinates, p=4.

4.2. Initial Displacement

The second example with 2D plane stress elements of a structure fixed on one side, that is subject to an initial displacement that contracts the structure and to be later released and allowed to vibrate. The resulting displacement field for 3 chosen Poisson ratios are evaluated. An initial case with a Poisson ratio of 0.3, typical of steel materials. A second Poisson ratio with value of 0, implying a disconnection from the displacements in the axial and transversal direction. Finally, a 3 scenario of Poisson ratio approaching 0.5, that implies incompressibility .

	$\begin{array}{l}E = 30 \times 10^{6}\end{array}$ $\begin{array}{l}t=1\end{array}$ $\begin{array}{l}h=2\end{array}$ $\begin{array}{l}L = 10\end{array}$ $\begin{array}{l}\nu = 0.3, 0.0, 0.499\end{array}$ $\begin{array}{l}\rho = 7200\end{array}$ $\begin{array}{l}\Delta t = 0.01\end{array}$ initial displacement: -0.25L
Plane stress structure fixed on left side, with varying Poisson ratio

In the graphics below, the results for the 3 test cases are collected, each of them with a bilinear elements version and a biquadratic spectral elements version. It can be observed in all cases fewer elements are needed to describe the displacement field when using a quadratic polynomial. Furthermore, in figure b) we observe that linear elements are able to capture the behavior since the deformation is linear, the quadratic elements are also able to depict this linear behavior. Finally, in case c) it is observed that the incompressibility patter is better shown with a higher polynomial, where the volume (area) of each element is maintained.

p=1
p=2
	a) plane stress elements with $\begin{array}{l}\nu = 0.3\end{array}$	b) plane stress elements with $\begin{array}{l}\nu = 0.3\end{array}$	c) plane stress elements with $\begin{array}{l}\nu = 0.499\end{array}$

In the previous example with the modal analysis of the beam, we stated that the accuracy is not significantly affected by the choice of either equidistant or Gauss-Lobatto coordinates to build the shape functions. In the following image we now compare the Computational Time in terms of the speed up obtained when using spectral elements, compared to classical p-refined elements, for various polynomial degrees and a Poisson Ratio = 0.3. It is observed that the speed up is always in a range between 1 and 1.5, for all the test cases analyzed. It can be therefore stated that for this test cases, spectral elements are always the better option to reduce computational time without sacrificing accuracy.

$\begin{array}{l}\displaystyle Speed Up=\frac{CPU \ \ Time original}{CPU\ \ Time improved}\end{array}$

Speed up for between classical and spectral elements for p=2 to p=5

5. Conclusions

1D and 2D spectral elements of arbitrary polynomial degree were implemented within the framework of BM-MFEM, specifically to work with the Newmark and Simple solver.
The spectral elements of higher polynomial provided a higher accuracy in the representation of higher mode shapes and eigenfrequencies.
Non-Equidistant Gauss Lobatto Coordinates and its respective quadrature provided more efficiency in solving, by providing a diagonal mass matrix.
- Speed up in the solving of the system of equations in the range of 1 to 1.5 for a given accuracy.

6. References

Cohen. "Higher-Order Numerical Methods for Transient Wave Equations," 2002.
A. Eisinberg, G. Fedele. "Discrete orthogonal polynomials on Gauss–Lobatto Chebyshev nodes," in Journal of Approximation Theory, vol. 144, no. 2, pp. 238–246, 2007.
Li, S. Dang, W. Li, and Y. Chai, “Free and forced vibration analysis of two-dimensional linear elastic solids using the finite element methods enriched by interpolation cover functions,” Mathematics, vol. 10, no. 3, p. 456, Jan. 2022.
Müller. Structural Dynamics Script, 2023.

Seitenhierarchie