Elsevier

Applied Mathematical Modelling

Volume 102, February 2022, Pages 115-136
Applied Mathematical Modelling

Frequency domain Bernstein-Bézier finite element solver for modelling short waves in elastodynamics

https://doi.org/10.1016/j.apm.2021.09.034Get rights and content

Highlights

  • BBFEM is extended to solve elastic wave problems.

  • Analytical rules are used to set up the local FE matrices for affine triangular elements.

  • The tensorial construction of Bernstein polynomials is exploited via sum factorisation for curved triangular elements.

  • The p BBFEM is more effective in coping with the pollution effect and resolving high order evanescent P and S wave modes.

  • Non uniform p refinement reduces the computational cost and enhances accuracy.

Abstract

This work presents a high-order Bernstein-Bézier finite element (FE) discretisation to accurately solve time harmonic elastic wave problems on unstructured triangular mesh grids. Although high-order FEs possess many advantages over standard FEs, the computational cost of matrix assembly is a major issue in high-order computations. A key ingredient to address this drawback is to resort to low complexity procedures in building the local high order FE matrices. This is achieved in this work by exploiting the tensorial property of Bernstein polynomials and applying the sum factorisation method for curved elements. An efficient implementation of the analytical rules for affine elements is also proposed. Furthermore, element-level static condensation of the interior degrees of freedom is performed to reduce the memory requirements. Additionally, the applicability of the method with a variable polynomial order, based on a simple a priori indicator, is investigated.

The computational complexities of sum factorisation, analytical rules and standard quadrature are first evaluated, in terms of the CPU time against the polynomial order. The analysis shows that the achieved numerical complexities compare favourably to those expected theoretically. A significant runtime saving is also obtained by using analytical rules and sum factorisation. The performance of the Bernstein-Bézier FEs is then assessed on various benchmark tests, over a wide range of frequencies. Results from the elastic wave scattering problem demonstrate the effectiveness of this method in coping with the pollution error, and its accuracy in resolving high order evanescent wave modes. Additionally, a wave transmission problem with high wave speeds contrast and a curved interior interface is considered, where a simple a priori indicator is proposed to assign the variable polynomial order. The study provides evidence of the great benefit of a non uniform p-refinement in reducing the computational cost and enhancing accuracy.

Introduction

Computer modelling of elastic wave propagation and scattering is an effective tool for predicting and testing in a wide variety of practical applications, including traffic-induced vibrations from roads and railways, seismic induced vibrations, seismic inversion, design of foundation elements, geophysical exploration, nondestructive evaluation and structural engineering.

The Finite Element Method (FEM) is a popular discretisation method commonly adopted for the numerical solution of wave problems. It is favoured over other methods, owing to its ability in accurately handling complex geometries, material heterogeneity and anisotropy. When typically piecewise linear FEs are used, around ten grid points per wavelength are needed to ensure a satisfactory resolution of the wave pattern. It is noteworthy that this rule of thumbnail of ten nodal points per wavelength enables to control the local approximation error [1]. However, due to the pollution effect in the mid or high frequency regimes [1], [2], [3], [4], even more than ten nodal points are required to achieve an acceptable level of accuracy, and thus the procedure becomes prohibitively expensive and less effective.

The past few decades have seen many attempts to design robust numerical schemes alleviating the pollution phenomenon. A substantial improvement was made in developing wave based, or Trefftz, methods enabling to achieve engineering accuracy with a significant reduction in the computational effort. The common idea in these approaches is to incorporate the oscillatory behaviour of the solution into the approximate space by including plane waves or Bessel functions. Examples of such methods, among others, are the Partition of Unity Finite Element Method (PUFEM) [5], [6], [7], [8], [9], [10], [11], [12], [13], [14], the Generalized Finite Element Method (GFEM) [15], [16], the Ultra Weak Variational Formulation (UWVF) [17], [18], [19], the least-squares method [20], the Discontinuous Enrichment Method (DEM) [21], [22] and the wave-based discontinuous Galerkin method [23], [24], [25]. It is worth noticing that all the aforementioned methods (except PUFEM, GFEM and DEM) share the same discontinuous Galerkin framework and differ only in the used numerical fluxes to ensure continuity between elements [26].

Wave based methods generally make use of the plane wave basis decomposition by pressure and shear waves in the approximation of time harmonic elastic wave problems. Perrey-Debain et al. [27] proposed the partition of unity boundary element method. Huttunen et al. [28], [29] developed the UWVF in two and three space dimensions. Also, Zhang et al. [30] extended the DEM to frequency domain elastic wave computations. Note that inter-element continuity in the UWVF and DEM is, respectively, enforced weakly by means of numerical fluxes and Lagrange multipliers. In References [31], [32], [33], [34], [35] and [36], PUFEM was successfully applied to time harmonic elastic wave problems, in two and three space dimensions. More recently, Yuan and Liu [37] proposed another approach in the discontinuous Galerkin setting based on plane wave basis and local spectral elements.

Higher-order polynomials such as those used in the hp-version of FEM and spectral elements are less vulnerable to the pollution effect [3], [4], [38], [39]. In the context of Helmholtz problems, Petersen et al. [40] assessed the efficiency of various p-FEM shape functions including Lagrange Gauss-Lobatto, integrated Legendre and Bernstein polynomials. They pointed out the advantage of high-order polynomials in reducing the pollution error and the good performance of Bernstein-Bézier FEs in conjunction with Krylov subspace solvers. Further, Lieu et al. [41] compared the performance of p-FEM, with Lobatto polynomials, and the wave-based discontinuous Galerkin method on various benchmarks. Similarly, El Kacimi et al. [42] compared the performance of the Bernstein-Bézier Finite Element method (BBFEM) and PUFEM. It was concluded from these comparative studies that high-order polynomial methods in combination with static condensation are able to deliver comparable or, in some cases, even superior performance. On the other hand, high order hp-discontinuous Galerkin methods have been demonstrated as competitive discretisation schemes (see, e.g., [43], [44] and the references cited therein), due to their flexibility in dealing with hp-adaptivity and general shaped elements, and their ability in achieving the exponential convergence rate of spectral techniques.

Another alternative to cope with the pollution error is the use of high-order isogeometric elements. The Isogeometric Analysis (IGA) [45] has gained popularity in recent years, due to a number of advantages afforded by spline basis functions, such as exact representation of geometries, higher order continuity and improved convergence properties. This methodology was successfully applied in a number of applications, as e.g. acoustic [46], [47], [48], [49], elastodynamic [50], [51], [52], [53] and electromagnetic [54], [55], [56] waves. Peake et al. [57], [58] extended the IGA concept within the framework of boundary elements for exterior acoustic problems, using the partition of unity in two and three space dimensions. They reported superior convergence compared to the partition of unity boundary element method. More recently, Ayala et al. [59] studied the performance of the enriched isogeometric collocation method in time-harmonic acoustics. They found improved convergence compared to the isogeometric collocation method. Willberg et al. [50] assessed the performance of isogeometric elements with non-uniform rational B-splines (NURBS) against high order SEM and p-FEM in Lamb wave propagation analysis. They concluded that higher order schemes deliver much improved accuracy, with a significant reduction in the total number of degrees of freedom (DoF), compared to conventional FEMs. Moreover, they pointed out that isogeometric elements enable to achieve high convergence rates.

In contrast to standard low order FEMs, where most of the CPU time is spent in the solution process, the computational burden of high order methods shifts to the evaluation of the element matrices and assembling (see, e.g., [60]). Therefore, the design of efficient quadrature rules is crucial to alleviate this major drawback. On tensor product elements, sum factorisation dating back to Orszag [61], has proven fundamental in the efficient implementation of SEMs. This technique relies mainly on a tensor-product structure of the shape functions and enables to achieve near optimal complexity. For simplicial elements, Karniadakis and Sherwin [62] designed special bases preserving a tensor product property. Recently, Ainsworth et al. [63] have shown that the Bernstein-Bézier basis naturally possesses the tensorial property needed for the sum factorisation method. They proposed an algorithm to set up the element mass and stiffness matrices with optimal complexity. Kirby and Thinh [64] developed low complexity matrix-free algorithm with Bernstein polynomials for applying the local finite element operators.

The main focus of this paper is to extend BBFEM to accurately solve time-harmonic elastic wave problems. As one of the key ingredients of the method, particular attention is paid to the computation of the element mass and stiffness matrices. Here, analytical rules are used for affine triangular elements and the sum factorisation method together with Stroud quadrature [76] is applied for elements with curved edges, where the geometry of such elements is interpolated via the linear blending map of Gordon and Hall [65], [66], [67]. The applicability of BBFEM with a variable polynomial order, using a simple a priori indicator, is also investigated throughout a benchmark dealing with the transmission of elastic waves.

The remainder of this paper is organized as follows. Section 2 introduces the model problem and its weak form. The Bernstein-Bézier FE approximation of the governing equations is presented in Section 3. Section 4 gives a description of the analytical and quadrature rules used for evaluating element integrals. Section 5 is devoted to numerical results. Finally, some concluding remarks are made in Section 6.

Notation

The following notation will be used throughout this paper. Each point x in R2 is identified by its components (x1,x2) relative to the Cartesian vector system denoted (e1,e2), i.e. x=x1e1+x2e2. The dot product of two vectors a and b in C2 is a scalar given by a·b=i=1,2aibi. The double dot product of two second-order tensors A and B in C2×2 is a scalar given by A:B=i,j=1,2AijBij. The double dot product of a fourth-order tensor CC2×2×C2×2 and a second-order tensor BC2×2 is a second order tensor D=C:B given by Dij=k,l=1,2CijklBkl. We denote the scalar product either in C2 or C2×2 by (·,·), that is, (a,b)=a·b¯ and (A,B)=A:B¯, where the notation ’¯’ refers to the complex conjugate. The induced norms in C2 or C2×2 will be denoted by ·.

We also denote the usual inner product on the complex-valued Sobolev space Hs(D), where s=0,1, by (·,·)s,D, with H0(D)=L2(D). We keep for simplicity the same notation for the inner products on the space of vector and tensor valued functions [Hs(D)]2 and [L2(D)]2×2, that is,(u,v)s,D=i=12(ui,vi)s,D,u,v[Hs(D)]2,(σ,τ)0,D=i,j=12(σij,τij)0,D,σ,τ[L2(D)]2×2.Likewise, for a given ΣD, the L2 inner products on L2(Σ) and [L2(Σ)]2 are denoted by (·,·)0,Σ. We further introduce the induced norm ·s,D in [Hs(D)]2 and the semi-norm |·|1,D in [H1(D)]2, where |v|1,D2=(v,v)0,D, =(1,2) is the gradient operator, and the superscript ’’ denotes the transpose.

Standard multi-index notation will be used. For αZ+3 and λR3, we set |α|=i=13αi, λα=i=13λiαi, α!=i=13αi! and (|α|α)=|α|!α!. If α,βZ+3 such that βα, i.e., βiαi, 1i3, (βα)=i=13(βiαi). We denote by eiZ+3 the multi-index whose the ith entry is unity and remaining entries are zero.

Section snippets

Mathematical model

Let Ω be a bounded Lipschitz domain in R2, consisting of a solid, homogeneous, isotropic and linear elastic material. We denote by Γ=Ω its boundary, n and t the outward unit normal and tangent vectors to Γ. In the absence of a source, the time harmonic Navier equation reads as [68]ρω2u·σ(u)=0inΩ.Herein, ω is the angular frequency, ρ is the material density, u is the displacement field and · is the divergence operator. The Cauchy stress tensor σ(u) is linearly related to the infinitesimal

Bernstein-Bézier FE approximation

Let Th be a conforming partition of the domain Ω into triangular elements such that Ω¯h=TThT, where Ωh is an approximation of Ω and h is the mesh size of Th given by h=maxTThhT, with hT=diam(T). Each TTh is the image of the triangular master element T^ defined byT^={ξ=(ξ1,ξ2):0ξ11,0ξ21ξ1},i.e. T=FT(T^), where FT is a suitable reference map. For p1, let Vh be the finite dimensional approximation space defined byVh={v[C0(Ω¯h)]2:vFT[Pp(T^)]2,TTh},where Pp(T^) is the space of

Computation of element integrals

This section is devoted to the computation strategies of the element matrices and right hand side vectors of the stated problem. For ease of presentation, we suppose the polynomial order p to be uniform.

Numerical results

In this section, numerical results are presented to assess the performance of the proposed approach on various two dimensional benchmark problems. The BBFEM is efficiently implemented via static condensation such that the internal DoF are eliminated at the element level. This leads to a condensed linear system involving only DoF attached the mesh skeleton. Once the solution related to the vertex and edge modes is computed, the internal DoF can be recovered at the post-processing stage by

Conclusions

In this paper, BBFEM has been extended to efficiently solve time-harmonic elastic wave problems, using unstructured triangular mesh grids. Key aspects of the method rely on the use of analytical rules to set up the local FE matrices for affine elements, while the tensorial construction of Bernstein polynomials amenable to sum factorisation is exploited for curved elements to speed up the assembly process. These yield a notable saving in terms of computational cost and hence afford high order

References (81)

  • G. Gabard

    Discontinuous galerkin methods with plane waves for time-harmonic problems

    Journal of Computational Physics.

    (2007)
  • A. El Kacimi et al.

    Numerical analysis of two plane wave finite element schemes based on the partition of unity method for elastic wave scattering

    Computers and Structures.

    (2010)
  • A. El Kacimi et al.

    Wavelet based ILU preconditioners for the numerical solution by PUFEM of high frequency elastic wave scattering

    Journal of Computational Physics.

    (2011)
  • K. Christodoulou et al.

    High-order finite elements for the solution of helmholtz problems

    Computers & Structures.

    (2017)
  • S. Petersen et al.

    Assessment of finite and spectral element shape functions for efficient iterative simulations of interior acoustics

    Computer Methods in Applied Mechanics and Engineering.

    (2006)
  • A. Lieu et al.

    A comparison of high-order polynomial and wave-based methods for helmholtz problems

    Journal of Computational Physics.

    (2016)
  • A. El Kacimi et al.

    Bernstein-Bézier based finite elements for efficient solution of short wave problems

    Computer Methods in Applied Mechanics and Engineering.

    (2019)
  • P.F. Antonietti et al.

    High-order discontinuous Galerkin methods for the elastodynamics equation on polygonal and polyhedral meshes

    Comput Methods Appl Mech Eng

    (2018)
  • T.J.R. Hughes et al.

    Isogeometric analysis: cad, finite elements, nurbs, exact geometry and mesh refinement

    Computer Methods in Applied Mechanics and Engineering.

    (2005)
  • H. Wu et al.

    Isogeometric finite element analysis of interior acoustic problems

    Applied Acoustics.

    (2015)
  • L. Coox et al.

    A performance study of NURBS-based isogeometric analysis for interior two-dimensional time-harmonic acoustics

    Computer Methods in Applied Mechanics and Engineering.

    (2016)
  • G.C. Diwan et al.

    Pollution studies for high order isogeometric analysis and finite element for acoustic problems

    Computer Methods in Applied Mechanics and Engineering.

    (2019)
  • C. Willberg et al.

    Comparison of different higher order finite element schemes for the simulation of lamb waves

    Computer Methods in Applied Mechanics and Engineering.

    (2012)
  • A. Idesman et al.

    Accurate solutions of wave propagation problems under impact loading by the standard, spectral and isogeometric high-order finite elements. comparative study of accuracy of different space-discretization techniques

    Finite Elements in Analysis and Design.

    (2014)
  • L. Dedè et al.

    Isogeometric numerical dispersion analysis for two-dimensional elastic wave propagation

    Comput Methods Appl Mech Eng

    (2015)
  • M.J. Peake et al.

    Extended isogeometric boundary element method (XIBEM) for two-dimensional helmholtz problems

    Computer Methods in Applied Mechanics and Engineering.

    (2013)
  • M.J. Peake et al.

    Extended isogeometric boundary element method (XIBEM) for three-dimensional medium-wave acoustic scattering problems

    Computer Methods in Applied Mechanics and Engineering.

    (2015)
  • T. Ayala et al.

    Enriched isogeometric collocation for two-dimensionaltime-harmonic acoustics

    Computer Methods in Applied Mechanics and Engineering.

    (2020)
  • J.C. Brigham et al.

    A spectral finite element approach to modeling soft solids excited with high-frequency harmonic loads

    Computer methods in applied mechanics and engineering.

    (2011)
  • S.A. Orszag

    Spectral methods for problems in complex geometries

    Journal of Computational Physics.

    (1980)
  • A. El Kacimi et al.

    Enhanced conformal perfectly matched layers for bernstein-bézier finite element modelling of short wave scattering

    Computer Methods in Applied Mechanics and Engineering.

    (2019)
  • P.R. Amestoy et al.

    Multifrontal parallel distributed symmetric and unsymmetric solvers

    Computer Methods in Applied Mechanics and Engineering.

    (2000)
  • F. Ihlenburg et al.

    Dispersion analysis and error estimation of Galerkin finite element methods for the Helmholtz equation

    Int. J. Numer. Methods Eng

    (1995)
  • F. Ihlenburg et al.

    Finite element solution of the Helmholtz equation with high wavenumber. part II: the hp-version of the FEM

    SIAM Journal on Numerical Analysis.

    (1997)
  • I. Babuška et al.

    Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave numbers?

    SIAM Review.

    (2000)
  • I. Babuška et al.

    The partition of unity method

    International Journal for Numerical Methods in Engineering.

    (1997)
  • O. Laghrouche et al.

    Short wave modelling using special finite elements

    Journal of Computational Acoustics.

    (2000)
  • P. Ortiz et al.

    An improved partition of unity finite element model for diffraction problems

    International Journal for Numerical Methods in Engineering.

    (2001)
  • O. Laghrouche et al.

    Modelling of short wave diffraction problems using approximating systems of plane waves

    International Journal for Numerical Methods in Engineering.

    (2002)
  • O. Laghrouche et al.

    Plane wave basis for wave scattering in three dimensions

    Communications in Numerical Methods in Engineering.

    (2003)
  • View full text