Three-dimensional dispersion analysis and stabilized finite element methods for acoustics

Accepted Manuscript Three-dimensional dispersion analysis and stabilised finite element methods for acoustics Fr´ed´eric Magoul`es, Hanyu Zhang

PII: DOI: Reference:

S0045-7825(18)30093-8 https://doi.org/10.1016/j.cma.2018.02.014 CMA 11789

To appear in:

Comput. Methods Appl. Mech. Engrg.

Received date : 25 September 2017 Revised date : 7 February 2018 Accepted date : 9 February 2018 Please cite this article as: F. Magoul`es, H. Zhang, Three-dimensional dispersion analysis and stabilised finite element methods for acoustics, Comput. Methods Appl. Mech. Engrg. (2018), https://doi.org/10.1016/j.cma.2018.02.014 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Three-dimensional dispersion analysis and stabilised finite element methods for acoustics Fr´ed´eric Magoul`es, Hanyu Zhang CentraleSup´elec, France and Providence University, Taiwan

Abstract Galerkin/least-squares and Galerkin gradient/least-squares stands out among several approaches designed to improve the numerical solution accuracy and counteract the pollution effect by adding terms to the standard Galerkin formulation. These added terms are multiplied by a ‘stability parameter’, which must be properly defined. In this paper, an original three-dimensional dispersion analysis is performed for the Helmholtz equation, together with the determination of the three-dimensional stability parameters for structured and unstructured meshes. Numerical experiments show the relative efficiency of the proposed methods for solving acoustic problems arising from industry. Keywords: stabilised finite element, Galerkin Least Squares, Galerkin Gradient Least Squares, dispersion, pollution effect, Helmholtz equation, acoustics

1. Introduction Boundary-value problems governed by the Helmholtz equation are important in a variety of applications involving time-harmonic wave propagation phenomena, as for example, computational acoustic. It is known and understood that the accuracy of the standard Galerkin finite element method deteriorates rapidly with increasing the wave number [17, 25]. This problem arises from the use of piecewise polynomial shape functions to approximate highly oscillatory wave propagation solutions and this effect is the so-called ‘pollution effect’ [2, 1, 15]. In fact, when applying numerical methods for the computation of solution to the Helmholtz equation, one obtains ‘numerical waves’ that is dispersive also in a non-dispersive media [24], which

Preprint submitted to Comput Methods Appl Mech Eng.

February 7, 2018

is called ‘dispersion’. For the Helmholtz equation, the pollution effect is directly related to the dispersion [30, 6]. In order to get an acceptable accuracy, more than ten elements per wavelength are required, as low-order piecewise polynomials may not capture the high frequency oscillation, which is known as the ‘rule of thumb’ [17]. However, in practice, this would lead to a cost prohibitively expensive, especially in three dimensions. Variety of enhancements are made to the standard finite element methods in order to circumvent the problem [20]. Stabilized methods and enrichment methods are two related categories of the most used strategies. Here we review some works that employs enrichment functions to the solution space to improve the accuracy. A possible approach is to employ residual-free bubbles (RBF) which is able to capture the singularity well [11], yet the bubbles should be an approximation to the element Greens function to remain effective. Later on a stabilized method of adjoint type that is related to RFB called the unusual stabilized finite element methods (US-FEMs) was proposed by Franca et al. [10]. Essentially this is the method based on piecewise linear shapes enriched by a space of bubble funcitons spanned with a single function on each element. In practical this method is form-identical to the Galerkin/least-squares (GLS) method [20]. More recently, Farhat propsed the discontinuous enrichment method [9] which assures high coarse-mesh accuracy. Enrichments in this method involves free-space solutions of the underlying partial differential equation, hence, they are obtainable, independent to element geometry and polynomial orders. When it is applied to the acoustic problem, a superior performance of the Discontinuous Galerkin Method (DGM) with Lagrange Multipliers derived from the discontinuous enrichment method over the GLS method with quadratic element on both structured and unstructured mesh is reported in the work of Harari et al. [23]. The similar techniques is also successfully applied to the unsteady advection-diffusion problems [3]. However, inter-element continuity enforced weakly by Lagrange multipliers in the DGM makes it more complex to implement compared to the GLS method, and we are more interested in the GLS method with linear cuboid elements as it eliminates completely the plane wave dispersion. It is also worthy to note that Boundary Element Method (BEM) is current widely used in the industry as well for some obvious advantages over the Finite Element Method. Both methods need to prepare an appropriate mesh of the domain under consideration for numerical analysis, however, the construction of the geometry is tedious, time consuming and creates inaccuracies [21], as for complex shape of geometry, the mesh can only be generated half-automatically. The problem is especially severe for the FEM as it needs to mesh the whole domain while the BEM only discretizes the surface of the domain for it’s theoretical foundation mainly relies on the Green’s formula that the integration over the volume can be described by integration 2

on the boundary. Therefore when the domain has low ratio of surface to volume, it is tended to choose the BEM. Besides, BEM is considered to be more precise generally. There are mainly two sources of errors in approximating the solution by FEM, both of them stem from the discretization, but are of different nature. One is the inaccuracy in geometry approximation by piecewise polynomial elements which may only be eliminated for certain shapes, another one is an unavoidable side effect of numerical analysis in transforming an infinite problem to a finite problem. Clearly the BEM suffers less from the second error after solving the problem on boundary, as one may calculate arbitrary interior region directly by the controlling equation. By the same reason, if one has to compute the solution only in a small region away from the domain, the better approach tends to be the BEM [8]. This is especially true for exterior acoustic problems, since there is no need to artificially set a DtN boundary [16] for finite element analysis. However, the BEM always generate a dense non-symmetric linear system, which may cause difficulties in solving it. There is an excellent work [27] done by Harari and Hughes in 1992 on the cost comparison between FEM and BEM. The situation of computational power is quite different now, still the work may provide some insights by the rigorous comparison, ”The work clearly demonstrates that finite element methods are economically competitive with boundary element methods for obtaining solutions to problems of time-harmonic acoustics.” Thus for the moderate size of problems in our work, the FEM possibly remains as an effective method to BEM. This paper investigates the performance of two stabilised finite element methods, namely the Galerkin/least-squares (GLS) [22, 18] and the Galerkin gradient/leastsquares (GGLS) [12, 26], which stands out among several approaches [14, 13] designed to improve the numerical solution accuracy. Both of them belongs to the class of residual based finite element method, which counteract the pollution effect by adding residuals of the Euler-Lagrange equation to the standard Galerkin formulation. In the case of GLS, the residuals are in least-squares form and in the case of GGLS, the added terms contain gradients. These added terms are multiplied by a meshdependent ‘stability parameter’, which must be properly defined. Harari and Hughes [17] determined in one dimension for a structured mesh, the stability parameter which leads to the corresponding GLS method for which the discrete solution is nodally exact for the plane-waves, regardless of the mesh refinement. Thompson and Pinsky [33] generalized in two dimensions for the Helmholtz equation, the procedure of finding the GLS stability parameter with a structured mesh, which yields to exact phase for the plane-waves directed along uniform mesh lines, thus nodally exact solution. Harari and Magoul`es [19] investigated numerically the per3

formance of both GLS and GGLS on unstructured meshes in two dimensions. In this paper, we extend the GLS and GGLS methods for the Helmholtz equation to three dimensions by generalizing the definition of the stability parameter. Despite some tentatives have been done in 3D, none of them proposed a detailled dispersion analysis for finite element and its extension to unstructured meshes. For instance, in [32], the authors studied separately the efficiency of applications of the Galerkin and Least Square methods on three-dimensional Poisson and Helmholtz equations, leaving out the possible combination of the two methods by adding a stability parameter. In [6], the authors presented a method to measure the dispersion related to the classical Galerkin method in 1D, 2D, 3D by two detailed examples of onedimensional and two-dimensional model problems respectively. They also compared the efficiency in 2D and studied the influence of the topological triangle elements, however, they only give a brief and incomplete dispersion result for three-dimensional problems. In [5], the authors proposed a twenty seven points finite difference method for three-dimensional Helmholtz equation by minimizing the numerical dispersion. They carried out the dispersion analysis as in [6], but opposite to our approach the authors changed the exact wavenumber and the numerical one, which may not harm the design of the finite difference method and even made the Taylor expansion relatively simpler, but did made the understanding of the numerical pollution indirect. Nonetheless there are also two drawbacks, first, by the nature of the finite difference method on structured meshes; second, there are up to seven parameters to determine instead of one in our finite element approach. In [28], the authors studied the GLS method for time harmonic Maxwell equation in 2D and 3D. They deduced the stability parameter by a dispersion analysis, however they only tested the parameter within uniform structured meshes. In the following we first derive, for the finite element method, the complete dispersion analysis and determine the GLS and GGLS stability parameter for the Helmholtz equation in three dimensions. Then we extend the computation of the stability parameter from node to element to improve the accuracy on unstructured mesh. We also present an asymptotic analysis followed by some error estimate related to the stability parameter. Finally, intensive numerical experiments arising from real problems illustrate and validate the robustness of the proposed stabilised finite element methods. The reminder of this paper is organized as follows. In Section 2 we first present the numerical methods, giving the strong and weak form of the boundary value problem. Then in Section 3, a full dispersion analysis is carried out in 3D with uniform hexahedron elements; the difference between the exact wavenumber and the numerical one is also given in an explicit form. Based on this dispersion analysis, 4

in Section 4 we deduce the formula of the stability parameter (defined on node) for GLS and GGLS, and present the process of passing the stability parameter to be used at the element level, upon with three element sizes’ definition. In Section 5 we report numerical tests that estimate the performance of GLS and GGLS on both structured and unstructured meshes. Section 6 concludes this paper. 2. Stabilised finite element methods Let us consider a d-dimensional open bounded domain Ω ∈ Rd with a smooth boundary Γ = ∂Ω, and the Helmholtz equation with homogeneous Dirichlet boundary conditions: Lu = f in Ω and u = 0 on Γ, with Lu = ∆u + k 2 u the Helmholtz operator, ∆u the Laplace operator, k the given wavenumber, f a prescribed function, and u the unknown acoustic pressure. Galerkin. The standard finite element method with Galerkin approximation consists to discretize the domain Ω into non-overlapping elements, and to write the Helmholtz equation in the weak form, i.e., involving a set of functions V h ⊂ H01 (Ω). This can be summarized as find uh ∈ V h such as: a(v h , uh ) = (v h , f ) ∀v h ∈ V h

(1)

where the weak operator a(., .) is defined as a(v, u) = (∇v, ∇u) − (v, k 2 u), with (., .) the L2 (Ω) inner product. Galerkin/least-squares. The Galerkin/least-squares (GLS) method, originally inspired from the streamline-upwind/Petrov-Galerkin (SUPG) [4] method for incompressible flows, is obtained by adding to equation (1) a residual term, equal to the Helmholtz operator applied to the test functions, and multiplied by a stability parameter τ : wh = v h + τ Lv h . The associated weak form contains now a residual-based term controlled by the stability parameter and can be written as: a(v h , uh ) + (Lv h , τ Luh )Ω˜ = (v h , f ) + (Lv h , τ f )Ω˜ ˜ is used to indicate that the integration is performed on the interior where subscript Ω domain. The definition of the stability parameter is the key ingredient of the GLS method. This parameter can be deduced by a dispersion analysis. In Section 4, for the three dimensional case we derive the parameter τ by eliminating spurious dispersion of plane waves of direction (θ, φ) in uniform structured meshes composed of elements of size h aligned in the direction of propagation. 5

Galerkin-gradient/least-squares. The Galerkin-gradient/least-squares (GGLS) method adds gradients to equation (1) as follows: a(v h , uh ) + (∇Lv h , τ G ∇Luh )Ω˜ = (v h , f ) + (∇Lv h , τ G ∇f )Ω˜ In the three dimensional case, eliminating spurious dispersion of plane waves of direction (θ, φ) in uniform structured meshes composed of elements of size h aligned in the direction of propagation, allows us to obtain the stability parameter. For GLS and GGLS methods, during the derivation of the stability parameter two points need to be outlined: the first one is the element size h and the second one is the direction (θ, φ). Since the stability parameter is obtained by dispersion analysis on uniform structured meshes, its natural element size is the element side length. This element size should be generalized for unstructured meshes and possible choice of element size for hexahedrons are investigated in Section 4. Regarding the direction (θ, φ), depending on the problem, the most suitable a priori direction is used to eliminate the spurious dispersion, as illustrated in Section 5. 3. Dispersion analysis It is known that for high wavenumbers of the Helmholtz equation, the Galerkin method is polluted by dispersion [6]. Lots of methods aim at reducing the dispersion to lower the so-called pollution effect. The dispersion analysis performed in this section has two roles: first it shows the dispersion in form of a series of powers of elements size and is an accurate indicator to the pollution error, second it actually sets the starting point for deducing the stability parameters. Exact solution. We consider the homogeneous Helmholtz equation ∆u + k 2 u = 0 in R3 , then the planar wave u(x, y, z) = cei(k1 x+k2 y+k3 z) with any constant c is a solution for any values (k1 , k2 , k3 ) lying on a sphere of radius k, k 2 = k12 + k22 + k32 . Numerical solution. We discretize the domain with a regular grid composed of cuboid elements of size hx × hy × hz , so the nodes form a set {(x, y, z) ∈ Ω = (rhx , shy , thz ) + (xref , yref , zref ), (r, s, t) ∈ Z3 }. where (xref , yref , zref ) is any given reference position, without loss of generality, it can be set to (0, 0, 0), thus a node’s coordinates can be simplified as (rhx , shy , thz ) 6

(r-1,s+1,t+1)

•

(r,s+1,t+1) (r+1,s+1,t+1)

(r-1,s,t+1)

•

•

(r,s,t+1)

(r-1,s-1,t+1)

•

•

•

•

(r,s-1,t+1)

•

•

•

(r+1,s,t+1)

•

(r+1,s-1,t+1)

•

•

•

(r-1,s+1,t) (r,s+1,t) (r+1,s+1,t) (r-1,s,t) • • (r+1,s,t) j • (r,s,t) (r-1,s-1,t) (r,s-1,t) (r+1,s-1,t)

•

(r-1,s+1,t-1)

•

•

•

•

(r,s+1,t-1)

(r-1,s,t-1)

•

(r-1,s-1,t-1)

•

(r+1,s+1,t-1)

•

(r,s,t-1)

•

•

•

•

(r+1,s,t-1)

(r,s-1,t-1)

•

•

(r+1,s-1,t-1)

• •

• • • •

•

• •

•

j

• •

• • •

•

•

• •

•

• • •

Figure 1: Stencil of an interior node with its eight surrounding elements in structured (left) and unstructured meshes (right)

only. Consider a stencil formed by an interior node j at coordinates (rhx , shy , thz ) with a set Pj = {pi , i = 1, ..., 26} of 26 surrounding nodes, as illustrated in Figure 1. Using linear interpolation leads to the linear system AU h = 0 where A is the acoustic impedance matrix and U h the vector P of nodal unknowns. To be more precise, we approximate the solution u by uh = nj=1 Ujh Nj where n is the number of nodes in the grid, Nj the linear shape function of node j and Ujh the j-th component of the nodal unknowns vector U h . Taking the weak form of the homogeneous Helmholtz equation leads after discretisation to the P coefficients of the impedance matrix obtained from ne e the assembly procedure Aij = e=1 Kij − k 2 Mije where ne denotes the number of elements e the node i and j belongs to, and Kije the elementary stiffness matrix and Mije the elementary mass matrix defined as Z Z e e Kij = (∇Ni , ∇Nj ) dω, Mij = Ni Nj dω Ωe

Ωe

The numerical solution is given by solving the system AU h = 0, which corresponds, for all the interior node j, to: ! X X i,j,k h h h (2) Ajp Up = Ar,s,t Ur+i,s+j,t+k = 0 Ajj Uj + p∈Pj

i,j,k=−1,0,1

where we use another notation Ai,j,k r,s,t to indicate non zero entries for the j-th row of matrix A. To be more specific, consider the interior node j, suppose its coordinate 7

is (rhx , shy , thz ) then all nonzero entries for the j-th row are Ajp where p ∈ Pj ∪ {j}. However, it would be difficult to map node p to its coordinate generally. Therefore we use another notation Ai,j,k r,s,t to directly relate the non zero entries to the coordinate ((r + i)hx , (s + j)hy , (t + k)hz ) of the corresponding nodes, notice that i, j, k ∈ {−1, 0, 1}. Then we can explicitly yields the 27 nonzero entries bound to the 27 nodes, 9 in three layers as shown on the left of Figure 1 by a standard assembly procedure [7], 0,1,1 1,1,1 −1,0,1 0,0,1 1,0,1 −1,−1,1 1,−1,1 Aj,top = {A−1,1,1 , A0,−1,1 r,s,t , Ar,s,t , Ar,s,t , Ar,s,t , Ar,s,t , Ar,s,t , Ar,s,t r,s,t , Ar,s,t }

0,1,0 1,1,0 −1,0,0 0,0,0 1,0,0 −1,−1,0 0,−1,0 1,−1,0 Aj,middle = {A−1,1,0 , Ar,s,t , Ar,s,t } r,s,t , Ar,s,t , Ar,s,t , Ar,s,t , Ar,s,t , Ar,s,t , Ar,s,t

1,1,−1 −1,0,−1 0,0,−1 1,0,−1 −1,−1,−1 0,−1,−1 Aj,bottom = {A−1,1,−1 , A0,1,−1 , Ar,s,t , Ar,s,t , Ar,s,t , Ar,s,t , A1,−1,−1 } r,s,t r,s,t , Ar,s,t , Ar,s,t r,s,t Ajp∈Pj ∪{j} = Aj,top ∪ Aj,middle ∪ Aj,bottom

Thanks to the symmetry, there are only four kinds of entries to calculate as follows: (i) the entry corresponding to the center node j of the reference element, A0,0,0 r,s,t = A0 =

8 h2x h2y + h2y h2z + h2x h2z 8k 2 hx hy hz − 9 hx hy hz 27

(ii) the six entries corresponding to the element’s faces’ center nodes, A0,0,±1 = Axy = r,s,t

2 h2x h2z + h2y h2z − 2h2x h2y 2k 2 hx hy hz − 9 hx hy hz 27

2 h2x h2y + h2y h2z − 2h2x h2z 2k 2 hx hy hz − 9 hx hy hz 27 2 2 2 2 2 2 2 hx hy + hx hz − 2hy hz 2k 2 hx hy hz A±1,0,0 = Ayz = − r,s,t 9 hx hy hz 27 (iii) the twelve entries corresponding to the element’s edges’ center nodes, A0,±1,0 = Axz = r,s,t

A0,±1,±1 r,s,t

1 h2y h2z − 2h2x h2y − 2h2x h2z k 2 hx hy hz = Ax = − 18 hx hy hz 54

1 h2x h2z − 2h2x h2y − 2h2y h2z k 2 hx hy hz = Ay = − 18 hx hy hz 54 2 2 2 2 2 2 2 1 hx hy − 2hx hz − 2hy hz k hx hy hz A±1,±1,0 = Az = − r,s,t 18 hx hy hz 54 (iv) and the eight entries corresponding to the reference corner nodes, A±1,0,±1 r,s,t

A±1,±1,±1 = Axyz = − r,s,t

1 h2y h2z + h2x h2y + h2x h2z k 2 hx hy hz − 36 hx hy hz 216 8

Error analysis between numerical and exact wavenumber. We now assume that the h h h numerical solution has a form uh (x, y, z) = cei(k1 x+k2 y+k3 z) with k1h = k h cos(φ)sin(θ), k2h = k h sin(φ)sin(θ), k3h = k h cos(θ) where φ ∈ [0, 2π), θ ∈ [0, π], which means h h h h = ei(k1 (r+i)hx +k2 (s+j)hy +k3 (t+k)hz ) . Substitution into equation (2) gives us Ur+i,s+j,t+k the relation: X i(k1h (r+i)hx +k2h (s+j)hy +k3h (t+k)hz ) Ai,j,k =0 (3) r,s,t e i,j,k=−1,0,1

which become after simplification

A0 + 8Axyz cos(k1h hx ) cos(k2h hy ) cos(k3h hz ) + 4Ax cos(k2h hy ) cos(k3h hz ) + 4Ay cos(k1h hx ) cos(k3h hz ) + 4Az cos(k1h hx ) cos(k2h hy ) + 2Axy cos(k3h hz ) + 2Ayz cos(k1h hx ) + 2Axz cos(k2h hy ) = 0 (4) We can now introduce a new theorem. Theorem 3.1. For the Galerkin method with hexahedron finite elements hx ×hy ×hz , where hx = β1 h, hy = β2 h and hz = β3 h with some constant β1 , β2 and β3 , then the numerical wavenumber k h is linked to the exact wavenumber k with the relation when kh → 0: k 3 h2 2 kh = k − (β cos4 (φ) sin4 (θ) + β22 sin4 (φ) sin4 (θ) + β32 cos4 (θ)) + O(k 5 h4 ) (5) 24 1 where θ ∈ [0, π] denotes the polar angle and φ ∈ [0, 2π) the azimuthal angle with the wave’s propagation direction. Proof. As the dispersion analysis shows, we have an implicit relation between the exact wavenumber k and the numerical wavenumber k h as the equation (4) suggests, from which we are going to simplify it by denoting fxh = cos(k1h hx ) fyh = cos(k2h hy ) fzh = cos(k3h hz ) Then we have γ0 + γx fxh + γy fyh + γz fzh + γxy fxh fyh + γyz fyh fzh + γxz fxh fzh + γxyz fxh fyh fzh k2 = (2 + fxh )(2 + fyh )(2 + fzh )

(6)

(7)

with the coefficients as follows −2 −2 −2 −2 −2 γ0 = 24(h−2 x + hy + hz ) γx = 12(hy + hz − 2hx )

−2 −2 −2 −2 −2 γy = 12(h−2 x + hz − 2hy ) γz = 12(hx + hy − 2hz )

−2 −2 −2 −2 −2 γyz = 6(h−2 x − 2hy − 2hz ) γxz = 6(hy − 2hx − 2hz )

−2 −2 −2 −2 −2 γxy = 6(h−2 z − 2hx − 2hy ) γxyz = −6(hx + hy + hz )

9

(8)

Now we replace fxh , fyh and fzh by their Taylor expansions using the expansion formula for cos(x) and using the expression of k1h , k2h and k3h upon k h cos2 (φ) sin2 (θ)(k h hx )2 cos4 (φ) sin4 (θ)(k h hx )4 + + O((k h hx )6 ) 2 24 sin2 (φ) sin2 (θ)(k h hy )2 sin4 (φ) sin4 (θ)(k h hy )4 h fy = 1 − + + O((k h hy )6 ) 2 24 2 h 2 4 h 4 cos (θ)(k h ) cos (θ)(k h ) z z fzh = 1 − + + O((k h hz )6 ) 2 24

fxh = 1 −

(9)

and inserting them into the equation (7) by dividing the common factor we get the numerator kh2 −

kh4 h2 2 (β1 cos2 (φ) sin2 (θ)(7 − cos(2φ) + 2 cos2 (φ) cos(2θ)) 48 + β22 sin2 (φ) sin2 (θ)(7 + cos(2φ) + 2 cos2 (φ) cos(2θ)) + 2β32 cos2 (θ)(3 − cos(2θ))) + O(kh6 h4 )

(10)

and the denominator 1 1 − kh2 h2 (β12 cos2 (φ) sin2 (θ) + β22 sin2 (φ) sin2 (θ) + β32 cos2 (θ)) + O(kh4 h4 ) 6

(11)

after evaluation of the division and multiplication of both side by h2 , we get: kh4 h4 2 (β1 cos4 (φ) sin4 (θ) + β22 sin4 (φ) sin4 (θ) + β32 cos4 (θ)) 12 + O(kh6 h6 )

k 2 h2 = kh2 h2 +

(12)

However what we need is the expansion of k h h as a series of powers of kh, thus we may write k h h as follows k h h = α0 + α1 kh + α2 (kh)2 + α3 (kh)3 + α4 (kh)4 + O((kh)5 ) Inserting it into the right side of the equation (12) then a comparison with its left side leads to α0 = 0, α1 = 1, α2 = 0, α4 = 0 and α3 = −

1 2 (β cos4 (φ) sin4 (θ) + β22 sin4 (φ) sin4 (θ) + β32 cos4 (θ)) 24 1

which gives us k h h = kh + α3 (kh)3 + O((kh)5 ) and concludes the proof. 10

(13)

(a) hx : hy : hz = 1 : 2 : 3

(b) hx : hy : hz = 3 : 1 : 2

(c) hx : hy : hz = 1 : 1 : 1

Figure 2: The −α3 for (θ, φ) ∈ [0, π2 ] × [0, π2 ] and different ratios among hx , hy and hz

It is quite clear now that no matter what the wave’s propagation direction is and no matter the ratios between hx , hy and hz are, there is always a phase lag for the numerical solution of the the Galerkin method. Here are three examples of −α3 in Figure 2. It is sufficient to show them in the range of (θ, φ) ∈ [0, π2 ]×[0, π2 ] since the rest parts are symmetric. It is obvious from the figure that the numerical wavenumber for the cuboid elements with different length of sides in the worst case could be even worse than the one for the cubic elements and it can be observed that the worst case happens when the propagation direction is aligned with the longest edge of the elements. However, for some particular directions, the numerical wavenumber could be closer to the exact one for the cuboid elements than that of the best case for cubic elements, which suggests that one may exploit the cuboid elements in some situations. In fact we have the following Lemma 3.2. For the Galerkin method with hexahedron finite elements as defined in Theorem 3.1, the difference between the exact wavenumber k and the numerical wavenumber k h depends on the wave’s propagation direction. In fact the coefficient α3 2 , − 24( 1 +11 + 1 ) ] with βmax = max{β1 , β2 , β3 }. Moreover belongs to the range [− βmax 24 2 β1

2

2 β2

2 β3

β2

β 2 (β 2 +β 2 )

α3 hits the upper bound while tan (φ) = β12 and tan2 (θ) = 3 β 21β 2 2 and α3 takes the 2 1 2 lower bound while the wave’s propagation is aligned with the longest edge.

11

Proof. Let us first prove the lower bound. β12 cos4 (φ) sin4 (θ) + β22 sin4 (φ) sin4 (θ) + β32 cos4 (θ) ≤ (β1 cos2 (φ) sin2 (θ) + β2 sin2 (φ) sin2 (θ) + β32 cos2 (θ))2 2 ≤ βmax (cos2 (φ) sin2 (θ) + sin2 (φ) sin2 (θ) + cos2 (θ))2 2 = βmax

(14)

2 /24. Now consider if βmax = β1 , then the equality holds which means α3 ≥ −βmax if and only if sin(φ) = cos(θ) = 0, similarly for βmax = β2 , then we need cos(φ) = cos(θ) = 0 and for βmax = β3 , only sin(θ) = 0 is required. In other words, the wave’s propagation direction is aligned with the longest edge.

In order to get the upper bound we use the Cauchy-Schwarz inequality (β12 cos4 (φ) sin4 (θ) + β22 sin4 (φ) sin4 (θ) + β32 cos4 (θ))(

1 1 1 + 2 + 2) 2 β1 β2 β3

≥ (cos2 (φ) sin2 (θ) + sin2 (φ) sin2 (θ) + cos2 (θ))2 =1

(15)

β2β2β2

1 2 3 which means α3 ≤ − 24(β 2 β 2 +β 2 β 2 +β 2 β 2 ) . The equality is attained when 1 2

1 3

2 3

β12 cos2 (φ) sin2 (θ) = β22 sin2 (φ) sin2 (θ) = β32 cos2 (θ) By the left equality, we have tan2 (φ) =

β12 β22

(16)

(17)

By the right equality and replacing sin2 (φ) by tan2 (φ)/(1 + tan2 (φ)), we get tan2 (θ) =

β32 (β12 + β22 ) β12 β22

(18)

which concludes the proof. Given the ratio β1 , β2 and β3 , It is possible to plot the surface representing the numerical wave number in the space. We are particularly interested in the case of β1 = β2 = β3 = 1 as shown in Figure 3 with k = 2 and h = 1, compared to the sphere of radius k = 2 representing the exact wave number. Furthermore we show in Figure 4 the evolutions (in percentage, i.e., 100(kh − k h h)/kh) of the phase error with k = 1 when h varies along some extreme directions. 12

(a) Front view

(b) Right view

(c) Left view

(d) Top view

Figure 3: Representation of the exact wave number (half sphere of radius) and numerical wave number (k h (φ, θ)) for solving the Helmholtz equation by the Galerkin method using linear shape functions with k = 2 and h = 1. The distance between the two surfaces varies with the angle φ and θ.

4. Stability parameter As we saw in the dispersion analysis for the Helmholtz equation, the classical Galerkin method inevitably introduced a numerical dispersion during the discretisation. Hence both GLS and GGLS methods add an additional term controlled by a stability parameter. The objective is now to find a suitable stability parameter to annihilate the dispersion with the homogeneous equation and uniform meshes. Stability parameter for GLS method. Theorem 4.1. For the homogeneous Helmholtz equation in three dimensions, with a structured mesh composed of hexahedron elements of hx × hy × hz , where hx , hy and hz denote the element size in x-axis, y-axis and z-axis respectively, choosing the stability parameter τ as: τ k2 = 1 −

γ0 + γx fx + γy fy + γz fz + γxy fx fy + γyz fy fz + γxz fx fz + γxyz fx fy fz (19) k 2 (2 + fx )(2 + fy )(2 + fz )

with the coefficients γ’s defined equation (8), and fx = cos(k1 hx ), fy = cos(k2 hy ), fz = cos(k3 hz ) with k1 = k cos(φ) sin(θ), k2 = k sin(φ) sin(θ), k3 = k cos(θ), the GLS method yields exact phase for the plane-waves directed along uniform mesh lines, and gives the exact solution. In particular, if hx = hy = hz = h, then the stability parameter can be reduced to τ k2 = 1 −

18 (4 − fx fy − fx fz − fy fz − fx fy fz ) h2 k 2 (fx + 2)(fy + 2)(fz + 2) 13

(20)

(a) Normal scale

(b) Log-Log scale

Figure 4: Percentage of error evolution with h for numerical wave number along different wave’s directions.

Proof. The so-called GLS method adds a residual term controlled by the stability parameter τ to the Galerkin weak form leading to Z Z Z h 2 h (∇u , ∇v) − k u v + τ (∆v + k 2 v)(∆uh + k 2 uh ) = 0 Ω

Ω

˜ Ω

Considering the linearity of the shape functions, and comparing to the weak form of equation (1), the only difference is the coefficient of the mass matrix ; all the procedure keeps the same, besides that we should replace k 2 by k 2 (1−τ k 2 ). Thus the left side of the equation (7) become k 2 (1 − τ k 2 ). In addition, the stability parameter τ is then by design to guarantee the nodal exact solution, which means we should also replace fxh , fyh , fzh by fx , fy , fz respectively in the equation (7). and taking its solution of τ , we can find the stability parameter τ as noted in equation (19). The reduction to equation (20) in the case of cubic elements is a direct calculation by taking hx = hy = hz = h within equation (19), which concludes the proof. Stability parameter for GGLS method. Theorem 4.2. For the homogeneous Helmholtz equation in three dimensions, with a structured mesh composed of hexahedron elements of hx × hy × hz , where hx , hy and hz denote the element size in x-axis, y-axis and z-axis respectively, choosing the stability parameter τ as: τ k4 =

k 2 (2 + fx )(2 + fy )(2 + fz ) − 1 (21) γ0 + γx fx + γy fy + γz fz + γxy fx fy + γyz fy fz + γxz fx fz + γxyz fx fy fz 14

where fx = cos(k1 h), fy = cos(k2 h), fz = cos(k3 h) and k1 = k cos(φ) sin(θ), k2 = k sin(φ) sin(θ), k3 = k cos(θ), the GGLS method yields exact phase for the planewaves directed along uniform mesh lines, and gives the exact solution. In particular, if hx = hy = hz = h, then the stability parameter can be reduced to τ k4 =

(fx + 2)(fy + 2)(fz + 2) h2 k 2 −1 18 (4 − fx fy − fx fz − fy fz − fx fy fz )

(22)

Proof. For the GGLS method, the added term controlled by the stability parameter τ contains gradients: Z Z Z h 2 h (∇u , ∇v) − k u v + τ (∇(∆v + k 2 v), ∇(∆uh + k 2 uh )) = 0 Ω

Ω

˜ Ω

By the same technique we used to prove Theorem 4.1, considering linear shape function and comparing to the weak form of equation (1), the only difference is the coefficient for the stiffness matrix, but it is equivalent by dividing the coefficient in case it is not zero since the equation is homogenous, so the difference reduces to the coefficient of the mass term, which should be k 2 /(1 + τ k 4 ) in place of k 2 . So we replace the left side of the equation (7) by k 2 /(1 + τ k 4 ) and with the same object – eliminating completely the dispersion – in mind, we also replace fxh , fyh , fzh by fx , fy , fz respectively in the equation (7). Solving it for the stability parameter leads us to the formula (21). Thus the nodal exactness is inherent from the construction. Again the reduction to the formula (22) in the case of cubic elements is a simple calculation by taking hx = hy = hz = h, which concludes the proof. Remark 4.3. This 3D stability parameter could be reduced for 2D’s and 1D’s usages, leading for instance in 2D for quadrangular elements to the same expression than in [33]. With the propagation direction fixed as in Figure 5 for the stability parameter τGLS , we could have a general sense of the variation of the stability parameter in the space of (hx , hy , hz ). It is not hard to notice that the level face is perpendicular to the propagation direction at their intersection while the total surface bends to the origin when the direction is not aligned to any of the axis, in addition their curvatures increase as they get closer to the origin. The figure for the stability parameter τGGLS for the GGLS method is quite similar to the one above. Asymptotic analysis. Figure 6 shows how the stability parameter evolves as the element’s size decreases. These curves can actually be extracted from Figure 5 along the 15

(a) θ = 0, φ = 0

(b) θ =

π 4,

(c) θ =

φ=0

π 4,

φ=

π 4

(d) θ =

π 6,

φ=

π 8

Figure 5: Stability parameter τGLS varies with hx , hy and hz and k = 1

(a) τGLS

k=1

(b) τGGLS

k=1

Figure 6: Stability parameter with h decreasing to 0

diagonal. Each curve corresponds to a diagonal for some fixed wave’s propagation direction. We observe that the curve first falls faster and faster as the element’s size increase and then rises even faster once the elements size h gets larger than certain critical point, which suggests that the curve is relatively stable to the change of element’s size for small values. Thus we may expect that a relative small perturbation in element’s size might cause a significant change when the element’s size is large, and hence requires more accuracy for calculating the element’s size. In the case of uniform structured meshes, it wouldn’t be a big problem in general, but while the stability parameter is used in unstructured meshes, the definition of element’s size becomes more crucial when it is large as the same level of difference to the potential optimal element’s size may cause larger difference in stability parameters, hence we may expect a slight performance improvements of GLS and GGLS methods with refinements of unstructured meshes.

16

Stability parameter from node to element. We saw that the stability parameter is defined for one particular interior node. The relative simple form of the node-based stability parameter takes great advantage of the uniform element size h, which leads to simplification. It will no longer be true for an unstructured mesh. In order to use GLS and GGLS on unstructured meshes, we need somehow to pass the node-based stability parameter to be used element-wise, and thus we generalize the definition of the element’s size on unstructured meshes. Consider an interior node j ; there are eight elements defining a set Ej = {ei , i = 1, . . . , 8} surrounding this node. The vertices of the elements of Ej define a set of nodes Pj ∪{j}. The discretised Helmholtz equation with GLS stabilisation writes as (K − (1 − τ k 2 )k 2 M )U h = 0 and the j-th row of the linear system writes X h (Kjm − (1 − τ k 2 )k 2 Mjm )Um =0 m∈Pj ∪{j}

e e If we note Kjm and Mjm the elementary contribution from element e to the entries Kjm and Mjm , then we have X X X e e h Kjm − (1 − τ k 2 )k 2 Mjm )Um =0 ( m∈Pj ∪{j} e∈Ej

e∈Ej

In this equation, if τ is the stability parameter of node j, i.e., τ = τj , then by h h,n construction of τj we have Um (= Um ) = Um , where here n indicates that the solution is obtained with node-based stability parameter, and we get P P e h,n m∈M e∈Ej Kjm Um (1 − τj k 2 )k 2 = P P h,n e m∈M e∈Ej Mjm Um

e e e e As Kjm = Kjm (he ) and Mjm = Mjm (he ), then τj is function of the mesh size of the surrounding elements, i.e., τj = τj (he1 , he2 , he3 , he4 , he5 , he6 , he7 , he8 ). If we have a structured mesh, i.e., he = h, then the formula reduces to τj = τj (h). If elements are unstructured, it is much more complicated to calculate, but replacing τ by τ e , the equation becomes X X X e h e h Um )=0 Kjm Um − (1 − τ e k 2 )k 2 Mjm ( m∈Pj ∪{j} e∈Ej

e∈Ej

17

h h,e If we take τ e = τj , then still we have Um (= Um ) = Um , where the superscript e indicates that the solution is obtained with element-based stability parameters. But there are also other possible values for τ e to validate the equality as well since X X e e ( Kjm Um − (1 − τ e k 2 )k 2 Mjm Um ) = 0 e∈Ej m∈Pj ∪{j}

As a consequence if we take τ e = τ e (he ) by the same formula for τj = τj (h), then h,e we may expect that the solution Um will be close to Um (as prove later within Theorem 4.4) . The advantage is that τ e only depends on he , therefore it is much simpler to compute, however the performance depends on the definition of he . Here ¯ cubic we consider three possible definition of element’s size, average side length (h), 1 root of volume (V 3 ) and the ratio of volume to the square of the average side length ¯ 2 ). (V /h We now introduces a definition of the irregularity of a mesh. Definition 4.1. Given a direction of wave propagation and a given definition of element’s size h, consider an interior node j with its surrounding nodes Pj = {pk | k = 1, . . . 26} belonging to Ej , each element having its own size hi according to the element’s size definition. We can construct a structured stencil composed of eight cubes of size hi along the direction of wave propagation, which gives 26 shifted nodes Pji = {pik | k = 1, . . . 26}. We define the irregularity of the interior node {j} along the given direction and the underlying element’s size as Ij = maxi,k ||pik − pk ||, where || · || denotes any vector norm. With this definition it is clear that a mesh is structured if and only if the irregularity of all interior nodes is equal to zero. We can now announce a new theorem: Theorem 4.4. Suppose that the solution U h,n (resp U h,e ) is obtained with GLS nodebased stability parameter τj (resp element-based stability parameter τ e ), then for any  > 0, there exists ξ > 0 such that ||U h,e − U h,n || <  as long as Ij < ξ for all interior nodes j. Proof. The system derived from the stabilised method with element-based stability parameters τ e , denoted as Ae , is perturbed by δA from the system A that is obtained with node-based stability parameters τj . The solutions U h,e and U h,n are different and the difference can be written as δU h = U h,e − U h,n . Thus we have the equation Ae U h,e = (A + δA)(U h,n + δU h ) = AU h,n

We first prove that ||δU h || is bounded by ||δA||. By simplification, we have δU h = −(A + δA)−1 δA U h,n 18

and considering the vector norm, we have ||δU h || ≤ ||(A + δA)−1 ||||δA|| = ||(I + A−1 δA)−1 ||||A−1 ||||δA|| h,n ||U || Since we have ||(I + A−1 δA)−1 || = ||I − A−1 δA + (A−1 δA)2 − (A−1 δA)3 + . . . || ≤

1 1 − ||A−1 δA||

as long as ||A−1 δA|| < 1, we obtain: ||δU h || ||A−1 ||||δA|| ≤ ||U h,n || 1 − ||A−1 δA|| Second we bound ||δA|| which is caused by the difference of two kinds of stability parameter. Consider an interior node j, its stability parameter τj is actually a continuous function of the coordinates of its 26 surrounding nodes. This means that we can compare it to the element-based stability parameter of each surrounding element e of node j. In fact, the element-based stability parameter τ e could also be regarded as the node-based stability parameter τje for the interior node j with identical cubic elements of size he surrounding it, thus according to the continuity of τj and the definition of irregularity, for any 0 > 0, there exists some δj > 0 so that |τ e − τj | = |τje − τj | < 0 as long as Ij < δj for all surrounding elements. Now consider the entries of j-th row of the system. The assembly procedure leads to |δAji | <

8 X e=1

|(τj − τ e )k 4 Mjie | < 0 |k 4 Mji | −

1−||A 1δA|| Finally, to have ||δU h || < , we need ||δA|| < η, with η = ||A − 1||||U h,n || . This will be satisfied if |δAji | < ξ, ∀j, i for some positive number ξ. From the inequality above, ξ we need to have 0 < |k4 M , ∀j, i, which is also satisfied as soon as δ is small enough. ji | This last condition could always be determined since there is a finite number of nodes.

Remark 4.5. The irregularity depends on the definition of element’s size, thus a good choice of the definition of element’s size makes a lot of difference since it has a small irregularity, which finally influences the perturbation of the system.

19

5. Numerical results Computations are carried out on structured and unstructured meshes to evaluate the numerical performance of the GLS and GGLS methods for several configurations with different kinds of boundary conditions. Special interests are given to the definitions of the elements size and wave’s propagation direction, on the efficiency of the stability parameters. We report the results of some of these numerical tests in the following section. We mainly consider the error caused by the dispersion which is introduced by the discretisation of the problem domain. When the wavelength is small, the mesh must be refined for the Galerkin method, otherwise, the dispersion is non-negligible and is the principle source of error. The stabilized method takes this phenomena into account and generates exact solution by design no matter how coarse is the mesh for plane wave on structured mesh. 5.1. Plane wave in a cuboid domain First we consider a cuboid domain, discretised with cuboid elements whose edges’ length are hx = 0.15, hy = 0.1 and hz = 0.25. The base discretisation consists of 8 × 8 × 8 elements, while the other two meshes as shown in Figure 7 are two successive refinements of the base one, which respectively consists of 16 × 16 × 16 and 32 × 32 × 32 elements.

Figure 7: Structured meshes (successive refinements) for the cuboid domain

The physical problem is defined by the Helmholtz equation with a wave number k = 5, without any distributed source. Inhomogeneous Dirichlet boundary conditions, i.e., u = ei(k1 x+k2 y+k3 z) with k1 , k2 , k3 as defined in Theorem 4.1 and Theorem 4.2, are imposed on the six faces of the domain, so the exact solution is known. We compared the numerical solutions of GLS and GGLS method and the Galerkin method to the exact solution with five different wave’s propagation directions, three 20

of them are chosen to be aligned with x, y and z axis, one is chosenp to minimise the phase error according to the form of the element, i.e., θ = arctan(β3 β12 + β22 /β1 β2 ) and φ = arctan(β1 /β2 ) as obtained in Lemma 3.2, and the other is chosen P to distinguish any of the previous. The error is measured by the norm ||V || = (( ni=1 Vi2 )/n)1/2 applied on U h − U . The results shows that both the GLS and GGLS method provides almost identical solution to the exact one for all five propagation directions with errors (not shown in the figure) of an order 10−8 . This superior performance of both GLS and GGLS along all directions confirms that the stability parameter defined equation (19) and (21) are well established. However, we would like to say a little more about the numerical solution by Galerkin method. We observe in Figure 8a that the smallest error doesn’t occur while the propagation direction is chosen to minimise the error as shown in Figure 8c, nor are some others cases. In other words, the order of errors of the numerical solution for the Galerkin method is not preserved as the order of relative phase error for chosen directions. This is because of the Dirichlet boundary condition. In fact, the dispersion analysis is carried out only for interior nodes of the mesh, hence the phase error is a precise indicator for the numerical solution only if we can solve the problem inside the mesh without any boundary condition. This can be simulated by imposing another inhomogeneous Dirichlet bounday condition, which use the numerical wavenumber k h instead. The corresponding result of error is shown in Figure 8b, it is clearly easy to notice the correspondence between Figure 8b and Figure 8c now. One more interesting observation is the following paradox, an inexact Dirichlet boundary condition may actually generate far better numerical solutions for the Galerkin method. In our experiments, the exactness of the Dirichlet boundary condition is obtained at the cost of more errors for interior nodes, and even worse, the errors accumulate. This phenomena seems to be inevitable for the Galerkin method since a proper modification in boundary condition is generally unknown, in contrast, both GLS and GGLS methods provides exact solution only when Dirichlet boundary conditions are exact, which is another advantages comparing to the Galerkin method for the paradox is automatically avoided. In the rest tests, experiments like Figure 8b are not considered since they are not known and hereafter we only test the reduced form of the stability parameter and its element based version in structured and unstructured meshes. Various boundary conditions and directions are investigated.

21

(a) Error with wavenumber = k (b) Error with wavenumber = k h

(c) Relative phase error

Figure 8: Plane wave in a cuboid domain along different wave’s propagation directions (θ, φ). Left: Errors ||U h − U || with Dirichlet boundary condition with ei(k1 x+k2 y+k3 z) , Center: Errors ||U h − U || h h h with Dirichlet boundary condition with ei(k1 x+k2 y+k3 z) , Right: Numerical wavenumber’s relative error ||k h − k||/||k||

5.2. Plane wave in a cube domain This domain consists of a cube of dimension a × a × a, meshed with two different kinds of meshes, each includes a set of three increasingly refined meshes. One set consists of nested structured meshes of 8×8×8, 16×16×16 and 32×32×32 elements, which means the element’s edges are halved from one level of refinement to the next. Figure 9 represents these meshes. The final set consists of three nested

Figure 9: Structured meshes (successive refinements) for the cube domain

unstructured meshes, i.e., the discretisation elements are irregular, but each has 22

a mean mesh size comparable to the uniform mesh. In the coarsest unstructured mesh, there are 2,744 nodes with 13 × 13 × 13 hexahedrons. The resolution in this mesh varies from 6.5 to 12.13 (with a mean of 9.4) points per wavelength. In the intermediate unstructured mesh, there are 19,683 nodes of 26 × 26 × 26 hexahedrons. The resolution in this mesh varies from 11.8 to 26.3 (with a mean of 19.0) points per wavelength. In the most refined unstructured mesh, there are 52 × 52 × 52 hexahedrons consisting of 148,877 node. The resolution in this mesh varies from 21.9 to 56.5 (with a mean of 38.3) points per wavelength. These meshes contains highly distorted elements with large variation in mesh size, aiming at testing the computational performance of the stability parameters under extreme conditions. Figure 10 represents these meshes.

Figure 10: Unstructured meshes (successive refinements) for the cube domain

The physical problem is defined by the Helmholtz equation with a wave number satisfying ka = 8 and a null distributed source. Various combinations of boundary conditions are imposed on the faces of the cube, so that the exact solution is a plane wave in the x-direction, i.e., uh = ei(k1 x+k2 y+k3 z) with φ = 0, ϕ = 0 and k1 , k2 , k3 as defined in Theorem 4.1 and Theorem 4.2. Appropriate inhomogeneous Dirichlet conditions are imposed on the six faces of the domain for the first set of tests. The results for the two kinds of meshes with Q-0-0 stability parameter are shown in Figure 11. The error is measured by the norm ||U h − U ||. The extremely good performance of the stabilised methods on structured meshes is a direct result of the design of stability parameters, which is devoted to eliminate 23

(a) structured mesh

(b) average size

(c) ratio size

(d) volume size

Figure 11: Plane wave in a cube: error for propagation at (0, 0) with Dirichlet-Dirichlet boundary condition and Q-0-0 stability parameter on structured (top, left) and unstructured meshes, with element size defined as the average side (top, right), the ratio of the volume to the square of the average side (bottom, left) and cubic root of the volume (bottom, right)

the dispersion phenomena in the wave’s propagation direction. On unstructured meshes, though the stabilised methods don’t perform as good as in the structured meshes, it still provides superior improvement comparing to the standard Galerkin method. This deterioration of performance mainly results from the highly distorted elements in these meshes. As Theorem 4.4 indicates the performance of elementbased stability parameter only tends to be as good as the exact solution when the mesh tends to be a structured one. Here GGLS performs better than GLS slightly, both of them improve the Galerkin method with the alternative three definitions of the element size, among which the average definition gives the best performance, typically reducing the error by over 80%. A comparison of the stability parameters with different φ and θ is reported in Figure 12. The results of the Q-0-0 parameter are uniformly the best whereas those of the Q- π4 - π4 parameter are invariably the worst. Results of the Q- π8 - π8 lie between the two extremes. It naturally matches our expectation since that closer the chosen direction to the wave’s propagation direction is, better the performance of the stabilised methods is. Once again, the GGLS method is a little better than the GLS method. And again the average definition of element size provides the best performance overall three directions. Typically the Q- π4 - π4 stability parameter reduces the error by almost 50% while the Q- π8 - π8 one reduces the error by about 65%. In all above cases, the unstructured meshes deteriorate the performance of the stabilised methods for the same reasons discussed earlier. Now we replace the Dirichlet boundary condition by an inhomogeneous Robin condition on the face whose normal vector is parallel to the x-axis and points to the positive direction, Neumman conditions are also imposed to the faces that are perpendicular to y-axis or z-axis, while the Dirichlet condition only retains on the last 24

(a) structured mesh

(b) average size

(c) ratio size

(d) volume size

Figure 12: Plane wave in a cube: error for propagation at (0, 0) with Dirichlet-Dirichlet boundary condition and various stability parameters on the coarsest structured (top, left) and unstructured meshes, with element size defined as the average side (top, right), the ratio of the volume to the square of the average side (bottom, left) and cubic root of the volume (bottom, right)

(a) structured mesh

(b) average size

(c) ratio size

(d) volume size

Figure 13: Plane wave in a cube: error for propagation at (0, 0) with Dirichlet-Neumann-Robin boundary condition and Q-0-0 stability parameter on structured (top, left) and unstructured meshes, with element size defined as the average side (top, right), the ratio of the volume to the square of the average side (bottom, left) and cubic root of the volume (bottom, right)

face. Results for the structured and unstructured meshes with the Q-0-0 stability parameter are reported in Figure 13 The presence of the Robin and Neumann boundary conditions degrades the performance of all methods, however, the stabilised methods still keep good performance on structured meshes, particularly the GLS method. Once again, both stabilised methods on the unstructured meshes, improve the performance with the three alternative definitions of the element size, and again the average element size provides the best improvement over the Galerkin method, typically reducing the error by over 60%.

25

Figure 14: Amphitheater - Computer Aided Model

5.3. Room acoustics Finite element analysis is an invaluable tool for investigating room acoustics. The usual approach is for the 3D-geometry to be obtained by Computer Aided Design (CAD) which is processed to produce smooth interpolated triangular surfaces. Then a volume mesh of solid tetrahedral or hexahedron finite elements is generated. An alternative approach is to use a voxel based finite element mesh where each voxel is converted directly into a finite element. Apart from this simplicity of model creation, working directly with voxel data offers a number of other advantages when modelling room acoustics, mainly: the voxels are converted into almost identical cuboid elements, with no meshing, distortion or scaling issues; resampling of the model is straightforward, allowing the resolution of the geometry and the model to be readily increased. Recently isogeometric analysis (IGA), which is more geometrically based and Computer Aided Design (CAD) friendly method, was proposed by Hughes et al [21]. The method is analogical to the FEM with some fundamental ideas. It is so called for it requires the solution space for unknows is represented in terms of the same functions which represent the geometry. One of its primary goal is to eliminate completely the geometrical error introduced by the discretization of the domain even when the discretization is very coarse [21]. As we stated in the section of introduction, the geometric approximation inherent in the mesh can induce accuracy problems, like in thin shell analysis, which is strongly sensitive to geometric imperfections; Hence it would be a better choice in practical use and it is actually the current tendance in the industry. NURBS (Non-Uniform Rational B-Splines) [31] are generally employed 26

to construct the exact geometric model. For the acoustic problem in the amphitheater and in the cavity of cars, the IGA would probably provide an improvement on the accuracy. However we are satisfied with relative easy implementation of voxel meshes and yet the IGA method still needs to discretize the domain, while the dispersion is inherently caused by the discretization, this kind of pollution is observed carried over to the NURBS-based approach [29], hence the stabilization technology of GLS method developed for FEM may also be transplanted for working around the problem in IGA; an issue under current investigation by the authors inspired by the work of Hughes et al. [21] of the application of SUPG within IGA for the advection–diffusion equation. In this section, the geometry consists of an amphitheater, see Figure 5.3 with a plateform and benchs. The problem is defined with a Dirac source located at the speaker position, absorbing material, modeled by Robin boundary conditions involving complex number, located on the roof of the amphiteater, and the floor, modeled with Neumann boundary conditions.

(a) average size

(b) ratio size

(c) volume size

Figure 15: Amphitheater: error with a Dirac source and Dirichlet-Neumann-Robin boundary condition and Q- π4 -0 stability parameter with element size defined as the average side (left), the ratio of the volume to the square of the average side (center) and cubic root of the volume (right)

The coarsest mesh contains 530 elements composed of 770 nodes, while the first level of refinement makes it 4,240 elements and 5,145 nodes, and the next level of refinement bring the number to 33,920 elements and 37,433 nodes, while the finest mesh which is used for calculating a Galerkin reference solution and consists of 271,360 elements and 285,201 nodes. With a wavenumber of ka = 10, the resolution for each mesh are respectively 3.1, 6.2, 12.4 and 24.9 points per wavelength, their elements are quite regular cuboid though not exact cubes. Results with stability parameter Q- π4 -0 can be found in Figure 15. Performance among the three alternative definitions of element’s size distinguish themselves very 27

slightly, this phenomena relies on the fact that the voxel element is so regular that their edge differentiate only by about 2%, leading to nearly same element size for all three size definitions. The stabilised methods continue to exhibit superior performance on these unstructured meshes. The GGLS method usually performs worse on coarser meshes. On finer meshes, the GLS method tends to be the same as the GGLS method, yet provide slightly better results. 5.4. Acoustic cavity Consider the interior of a car compartment, as illustrated in Figure 16. We report the results obtained with a wavenumber of ka = 2.8. The vibration of the firewall are modelled by an inhomogeneous Dirichlet boundary condition, and the absorbing roof is modelled by inhomogeneous Robin conditions involving complex number, while metal part of the structure are modelled by Neuman boundary conditions.

Figure 16: Meshes (three refinements) for a domain of an automotive interior

The domain is discretised using a set of increasingly refined unstructured meshes. In the coarsest mesh (not shown), there are 294 nodes within 151 hexahedrons, the resolution of this mesh ranges from 6.08 to 16.13 (with a mean of 8.3) points per wavelength. The lower limit of resolution in this mesh is far less than adequate for resolving the high wave number problems. After the first level of refinement, the mesh contains 1,727 nodes with 1,208 hexahedrons, the resolution in this mesh ranges from 11.65 to 38.125 (with a mean of 16.6) points per wavelength. The next refinement gives a mesh containing 11,637 nodes with 9,664 hexahedrons, the resolution in this mesh ranges from 22.875 to 81.5 (with a mean of 33.2) points per wavelength. The last mesh is used for compute the Galerkin reference solution, and consists of 85,001 nodes in 77,312 hexahedrons, the resolution ranging from 45.25 to 169.12 (with a mean of 67.1) points per wavelength.

28

(a) average

(b) ratio

(c) volume

Figure 17: Automotive: error for radiation in an automotive with Q-0-0 stability parameter and element size defined as the average side (left), the ratio of the volume to the square of the average side (center) and cubic root of the volume (right)

Results obtained with various definitions of the stability parameter of the three size definitions are reported in Figure 17. The same pattern observed previously appears again, i.e., that the average size definition performs best while the ratio of volume to square of average size definition works worst among the three. The stabilised methods consistently exhibit considerable improvement over the Galerkin method. The GGLS method provides better results over the GLS method with average size definition, while the GLS method prevails on the other two size definitions for the Q-0-0 stability parameter. An unexpected phenomena appears in Figure 17a, i.e., both curves of the GLS method and the GGLS methods bends to down side while the other two bends to up side as the curve of the Galerkin solution does. This is because that the stabilised methods with the Q-0-0 stability parameter defined with average size definition actually performs so good that they even beat the Galerkin solution on the next refined mesh, we need a more accurate reference solution to eliminate this phenomena. Opposite, the results with stability parameter Q- π4 -0 shown in Figure 18 don’t have such a problem which confirms the theory. The Q- π6 -0 and Q- π8 -0 stability parameters defined with the three element’s size definitions are also considered. Their results are similar to those in Figure 18, except that they perform better. A comparison of the four stability parameters on the intermediate refined mesh is reported in Figure 19. The results of the Q-0-0 parameter are uniformly the best (over 80% improvement), while the other three stability parameters ameliorate the performance over the Galerkin solution respectively by about 27%, 35% and 44% for ratio size definition, 35%, 46% and 58% for cubic root of volume size definition and 40%, 53% and 67% for average size definition. We observe that smaller the φ value is, better the performance is, which may convince 29

(a) average

(b) ratio

(c) volume

Figure 18: Automotive: error for radiation in an automotive with Q- π4 -0 stability parameter and element size defined as the average side (left), the ratio of the volume to the square of the average side (center) and cubic root of the volume (right)

(a) average

(b) ratio

(c) volume

Figure 19: Automotive: error with Dirichlet-Neumann-Robin boundary condition and various stability parameters on the intermediate coarse meshes with element size defined as the average side (left), the ratio of the volume to the square of the average side (center) and cubic root of the volume (right)

us that the φ value of the wave’s propagation direction is 0 or very close to it. This pattern is quite similar to the cube tests. 6. Conclusion This paper investigates the derivation of the stability parameters involved within the Galerkin/least-squares (GLS) and the Galerkin gradient/least-squares (GGLS) for the three-dimensional Helmholtz equation. After a detailed dispersion analysis, we have shown than GLS and GGLS have quite similar performance, and both of them are superior to the standard Galerkin method with a significant improvement 30

for uniform meshes. Though the performance of both methods slightly degrade on unstructured meshes, they still considerably make the solution better to varying extents with different element’s size definition. Among the three element’s size definitions we have proposed, the average edge length works best in general, the cubic root of the volume is also a good alternative, but the ratio of volume to the square of average size should probably be avoided. The extension of the node-based stability parameter to the element-based stability parameter we have proposed is a good compromise between the speed of calculation with easy integration of the stabilised method and the numerical performance.

31

[1] I. Babuˇska, F. Ihlenburg, E.T. Paik, and S.A. Sauter. A generalized finite element method for solving the Helmholtz equation in two dimensions with minimal pollution. Computer Methods in Applied Mechanics and Engineering, 128(3):325–359, 1995. [2] I. Babuska and S.A. Sauter. Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave numbers? SIAM review, 42(3):451– 484, 2000. [3] Raunak Borker, Charbel Farhat, and Radek Tezaur. A high-order discontinuous Galerkin method for unsteady advection–diffusion problems. Journal of Computational Physics, 332:520–537, March 2017. [4] A. Brooks and T.J.R. Hughes. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Computer Methods in Applied Mechanics and Engineering, 32(1):199–259, 1982. [5] Zhongying Chen, Dongsheng Cheng, and Tingting Wu. A dispersion minimizing finite difference scheme and preconditioned solver for the 3D Helmholtz equation. Journal of Computational Physics, 231(24):8152–8175, 2012. [6] A. Deraemaeker, I. Babuˇska, and P. Bouillard. Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions. International Journal for Numerical Methods in Engineering, 46(4):471–499, 1999. [7] G. Dhatt, G. Touzot, and E. Lefrancois. M´ethode des ´el´ements finis. Lavoisier, pages 143,269–282, 2005. [8] S. Falletta, G. Monegato, and L. Scuderi. On the discretization and application of two space-time boundary integral equations for 3d wave propagation problems in unbounded domains. Applied Numerical Mathematics, 124, 2017. [9] Charbel Farhat, Isaac Harari, and Leopoldo P. Franca. The discontinuous enrichment method. Computer Methods in Applied Mechanics and Engineering, 190(48):6455–6479, September 2001. [10] Leopoldo P. Franca, Charbel Farhat, Michel Lesoinne, and Alessandro Russo. Unusual stabilized finite element methods and residual free bubbles. International Journal for Numerical Methods in Fluids, 27(1-4):159–168, January 1998. 32

[11] Leopoldo P. Franca, Charbel Farhat, Antonini P. Macedo, and Michel Lesoinne. Residual-free bubbles for the Helmholtz equation. International Journal for Numerical Methods in Engineering, 40(21):4003–4009, November 1997. [12] L.P. Franca and E.G. Dutra Do Carmo. The Galerkin gradient least-squares method. Computer Methods in Applied Mechanics and Engineering, 74(1):41– 54, 1989. [13] L.P. Franca, G. Hauke, and A. Masud. Revisiting stabilized finite element methods for the advective–diffusive equation. Computer Methods in Applied Mechanics and Engineering, 195(13):1560–1572, 2006. [14] T.P. Fries and H.G. Matthies. A review of Petrov-Galerkin stabilization approaches and an extension to meshfree methods. Technical report, Institute of Scientific Computing, Univ. of Braunschweig–Inst. of Technology, Brunswick, Germany., 2004. [15] K. Gerdes and F. Ihlenburg. On the pollution effect in FE solutions of the 3DHelmholtz equation. Computer Methods in Applied Mechanics and Engineering, 170(1):155–172, 1999. [16] Dan Givoli and Joseph B. Keller. A finite element method for large domains. Computer Methods in Applied Mechanics and Engineering, 76(1):41–66, November 1989. [17] I. Harari and T.J.R. Hughes. Finite element methods for the Helmholtz equation in an exterior domain: model problems. Computer Methods in Applied Mechanics and Engineering, 87:59–96, 1991. [18] I. Harari and T.J.R. Hughes. Galerkin/least-squares finite element methods for the reduced wave equation with non-reflecting boundary conditions in unbounded domains. Computer Methods in Applied Mechanics and Engineering, 98(3):411–454, 1992. [19] I. Harari and F. Magoul`es. Numerical investigations of stabilized finite element computations for acoustics. Wave Motion, 39(4):339–349, 2004. [20] Isaac Harari. A survey of finite element methods for time-harmonic acoustics. Computer Methods in Applied Mechanics and Engineering, 195(13):1594–1607, February 2006.

33

[21] T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics & Engineering, 194(39):4135–4195, 2005. [22] T.J.R. Hughes, L.P. Franca, and G.M. Hulbert. A new finite element formulation for computational fluid dynamics: Viii. the Galerkin least squares method for advective-diffusive equations. Computer Methods in Applied Mechanics and Engineering, 73:173–189, 1989. [23] R. Tezaur I. Harari and C. Farhat. A study of higher-order discontinuous galerkin and quadratic least-squares stabilized finite element computations for acoustics. Journal of Computational Acoustics, 14:1–19, 2006. [24] F. Ihlenburg and I. Babuska. Dispersion analysis and error estimation of Galerkin finite element methods for the Helmholtz equation. International Journal for Numerical Methods in Engineering, 38(22):3745–3774, 1995. [25] F. Ihlenburg and I. Babuˇska. Finite element solution of the Helmholtz equation with high wave number part i: The h-version of the FEM. Computers and Mathematics with Applications, 30(9):9–37, 1995. [26] F. Ilinca and J.F. H´etu. Galerkin gradient least-squares formulations for transient conduction heat transfer. Computer Methods in Applied Mechanics and Engineering, pages 3073–3097, 2002. [27] Isaac Harari and Thomas J. R. Hughes. A cost comparison of boundary element and finite element methods for problems of time-harmonic acoustics. Computer Methods in Applied Mechanics and Engineering, 97(1):77–102, 1992. [28] J. Jagalur-Mohan, G. Feij´oo, and A. Oberai. A Galerkin least squares method for time harmonic Maxwell equations using N´ed´elec elements. Journal of Computational Physics, 235:67–81, 2013. [29] T. Khajah, X. Antoine, and S. P. A. Bordas. Isogeometric finite element analysis of time-harmonic exterior acoustic scattering problems. ArXiv e-prints, October 2016. [30] K.C. Patidar. Dispersion induced by the pollution for the wave equation. Applied Mathematics and Computation, 160(2):329–341, 2005. [31] Les Piegl and Wayne Tiller. The NURBS Book. Monographs in Visual Communication. Springer-Verlag, Berlin Heidelberg, 2 edition, 1997. 34

[32] EC Rom˜ao, MD Campos, and LFM Moura. Application of the Galerkin and least-squares finite element methods in the solution of 3D Poisson and Helmholtz equations. Computers and Mathematics with Applications, 62(11):4288–4299, 2011. [33] L.L. Thompson and P.M. Pinsky. A Galerkin least squares finite element method for the two-dimensional Helmholtz equation. International Journal for Numerical Methods in Engineering, 38(3):371–397, 1995.

35

The highlights of this paper are the following: • Complete dispersion analysis for 3D finite element • Determination of relation between exact and numerical wavenumber for 3D finite element • Evaluation of bounds for the error between exact and numerical wavenumber • Derivation for GLS method of 3D stabilized parameter in structured/unstructured mesh • Derivation for GGLS method of 3D stabilized parameter in structured/unstructured mesh • Justification for the computation of stabilized parameter within elements, i.e., not at nodes, in 3D • Asymptotic analysis upon mesh size • Numerical investigation and analysis upon several boundary conditions • Numerical experiments for solving acoustic problems arising from industry

1