Studies of an infinite element method for acoustical radiation

Jean-Christophe Autrique a

, Fre´de´ric Magoule`s b,*

a LMS International, Interleuvenlaan 70, Researchpark Haasrode Z1, 3001 Leuven, Belgium b Institut Elie Cartan de Nancy, Universite´ Henri Poincare´, BP 239, 54506 Vandoeuvre-les-Nancy Cedex, France

Received 11 July 2005; received in revised form 2 August 2005; accepted 8 August 2005

Available online 22 December 2005

Abstract

Infinite element computations are very efficient for predicting the vibro-acoustic response and sensitivities of a vibrating

structure for an exterior acoustic domain. In addition, domain decomposition methods are very powerful algorithms for

solving large linear systems in parallel. In this paper, an infinite element method is proposed and analyzed for parallel com￾putations purpose. An original formulation of this method with Lagrange multipliers defined on (semi-)infinite space is

presented. The implementation aspects of this method in an industrial acoustic software (SYSNOISE) are discussed.

New numerical results illustrate the efficiency of the proposed method for realistic acoustical radiation problems.

  2005 Elsevier Inc. All rights reserved.

Keywords: Infinite element; Finite element; Acoustic radiation; SYSNOISE; Domain decomposition method; Parallel computing

1. Introduction

The infinite element method [1–3] is an elegant extension of the finite element method [4–6] which allows for

the modelling of exterior acoustic problems. Such exterior problems involve unbounded media and require an

appropriate treatment of the Sommerfeld radiation condition. The formulation of infinite elements relies on a

truncated multipole expansion [7,8] of the acoustic field outside a regular convex surface enclosing acoustic

sources. Such a multipole expansion can be expressed in various separable coordinate systems as spherical,

spheroidal or ellipsoidal coordinates, for instance. In contrast with conventional finite elements [4–6], infinite

elements [3,9] incorporate frequency dependent interpolation functions along the radial (outward) direction.

Interpolation (or trial) functions are supplemented by appropriate test functions: the particular infinite ele￾ments available in SYSNOISE software rely on a conjugated formulation where test functions are basically

complex conjugates of trial functions [10]. This particular choice avoids highly oscillating functions in element

integrals, and so far standard Gauss–Legendre quadratures can be used for integrating element contributions.

An additional benefit of such conjugated infinite elements is related to the particular frequency dependence of

element contributions which remains polynomial as for conventional acoustic elements. These features allow

for the use of efficient solution techniques as iterative Krylov methods for instance [11]. The linear system

obtained after the discretization consists of a sparse matrix. The dimension of this matrix increases proportion￾ally with the number of finite and infinite elements, and with the order of the infinite elements too. When dealing

with large size problems and high frequencies, the iterative methods meet some difficulties to converge [11].

Generally speaking, the solution of acoustic problems, involving large meshes and high order elements,

requires a large amount of memory and a high computation time. The domain decomposition methods allow

to solve a large problem in a more reasonable amount of time [12–16]. Indeed, the discrete linear system

obtained from a finite or an infinite element method over a large mesh leads to a matrix with a large bandwidth.

If the physical mesh is partitioned into several sub-domains and if each sub-domain is allocated to one process,

the local bandwidths are smaller than the whole model bandwidth [12,17,18]. This creates a gain both in mem￾ory and computing time. A reformulation of the mathematical problem with Lagrange multipliers leads to a

linear system defined on the interface between the sub-domains. In SYSNOISE [19], two iterative solvers are

available to solve this linear system defined on the interface. The BiCGstab [20] algorithm is based on the con￾jugate gradient method, and the GMRES [20] (generalized minimum residual) algorithm is designed to solve

efficiently non-symmetric linear system. Communications between the different processes, involved during

the iterative algorithm, are handled with a message passing interface library (MPI). The speed-up obtained

by parallelization depends on the size of the interface between the sub-domains and the speed of the network.

In this paper, an infinite element method is proposed and analyzed for parallel computations purpose. The

implementation aspects in an industrial software (SYSNOISE) are then discussed. The scope of this paper is as

follows. Section 2 describes the general radiation problem analyzed in the following. Then, in Section 3 the

conjugate infinite element method is reminded. Section 4 presents a domain decomposition method well suited

for parallel acoustics computations. An original formulation of this method with Lagrange multipliers defined

on (semi-)infinite space is introduced. In Section 5, implementation aspects of this method in an industrial

software (SYSNOISE) are detailed. In Section 6, new numerical experiments are presented on large compu￾tational acoustics problems which demonstrate the performance and robustness of the proposed domain

decomposition method. This analysis investigates the dependency of the domain decomposition method upon

different parameters for general mesh partitioning obtained with the METIS software. Three-dimensional

analysis performed on academic and industrial test cases are presented. The conclusions of this paper are pre￾sented in Section 7.

Acoustical radiation problem

A model radiation problem is considered in an unbounded domain. The main motivation for this analysis is

to determine the frequency response functions arising from the vibrations of a structure. In the following the

radiation of an object delimited by a boundary CN immersed in an unbounded domain Xe is considered. This

model problem can be expressed as the following system of equations:

where g 2 L2

ðCNÞ is the prescribed Neumann boundary conditions and k 2 Rþ the wave number. The normal

unitary vector along the boundary CN is denoted by n, and r represents the radius in the spherical coordinates.

3. Infinite element method

3.1. Principles

In order to apply the infinite element method to the analysis of an acoustical radiation problem, involving a

non-convex object like an engine for instance, a convex envelope surrounding this object should first be

defined. The volume between the object and the convex envelope is then meshed with finite elements [4–6] and

infinite elements [9,10] are defined on the surface of the convex envelope. An example of such a mesh is shown

in Fig. 1 for a two-dimensional case. In the general case, the first step consists of defining a truncation of the

unbounded domain Xe called Xi;e

c as

Xi;e

c ¼ Xe \ fx 2 R3

; jxj < cg;

where the artificial boundary Sc (here, the sphere of radius c > 1) has been introduced. The bounded domain

Xi;e

c is meshed with finite elements. The exterior of the domain Xi;e

c is then defined by

Xo;e

c ¼ fx 2 R3

; c < jxjg

and is meshed with infinite elements ðXe ¼ Xi;e

c [ Xo;e

c Þ. The optimal distance between the object and the arti￾ficial boundary Sc depends of the order of the infinite elements considered [1–3].

3.2. Variational formulation

The variational formulation of the problem differs in the domain Xi;e

c and in the domain Xo;e

c .

In the domain Xi;e

c , the Helmholtz equation is first multiplied by the complex conjugate of the test function v

(noted vv).

 The integration in the domain Xi;e

c is then performed, and the Green formula is applied. The solution

u belongs to the space

H1

ðXi;e

c Þ¼fu : kuk1 < 1g

with kuk1 the norm associated to the scalar product

4. Parallel computing

4.1. Non-overlapping Schwarz algorithm

In order to solve the previous linear system, the non-overlapping Schwarz algorithm with absorbing bound￾ary conditions defined on the interface between the sub-domains is considered [21,22]. The case of a general

domain X partitioned into two sub-domains X(1) and X(2) with an interface C is now considered, as shown in

Fig. 2. The degrees of freedom located inside sub-domain X(s)

, s = 1, 2 and on the interface C are denoted by

subscripts i and p. With this notation the contribution of sub-domain X(s)

, s = 1, 2 to the impedance matrix

and to the right-hand side can be written as in [12,13,17]

where SðqÞ ¼ ZðqÞ

pp   ZðqÞ

pi ½ZðqÞ

ii 

11

ZðqÞ

ip is the condensed matrix and cðqÞ

p ¼ bðqÞ

p   ZðqÞ

pi ½ZðqÞ

ii 

11

b

ðqÞ

i is the condensed

right-hand side, for q = 1, 2. As already discussed, this linear system is solved with an iterative method, and

each iteration involves a solution with a direct method of an Helmholtz sub-problem in each sub-domain.

The choice of the matrices A(1) and A(2) has a strong influence on the convergence speed of the non-over￾lapping Schwarz algorithm. Different choice of these matrices has been investigated in [22,23,25,26]. In the

following the matrices A(1) and A(2) are obtained from a Taylor zeroth order approximation of the Stek￾lov–Poincare´ operator and from an optimized zeroth order approximation of the Steklov–Poincare´ operator,

as introduced in [21]. These matrices are equal to

Að1Þ ¼ aMC; Að2Þ ¼ aMC;

where a is equal to ik for a Taylor zeroth order approximation and obtained from the solution of a minimi￾zation problem for an optimized zeroth order approximation [21,18]. The matrix MC is a surface mass matrix

defined on the interface between the sub-domains.

4.2. Iterative solution of the interface problem

An iterative algorithm is considered for the solution of the interface problem, namely the GMRES(m). At

each iteration of this algorithm, only one matrix-vector product by the matrix F is required.

The conjugate gradient method is designed for Hermitian positive definite matrices and the work per iter￾ation is usually dominated by the matrix-vector product with F. The storage of five vectors of the dimension of

k is required. The method is optimal since it minimizes the error in the F-norm and has a smooth convergence

behavior. For non-Hermitian matrices, other methods must be used. The GMRES(m) method, proposed by

Saad [20], keeps the property of optimal and smooth convergence behavior, but the memory requirements can

be large, since the basis vectors of a large Krylov space should be stored. The GMRES algorithm have shown

robust convergence for acoustics problem [27,28]. One step of GMRES(m), as defined in the following algo￾rithm, by 3.1–3.7 requires one matrix-vector product. The quantity xH denotes the complex conjugate trans￾pose of x. The quantity TOL is the residual tolerance used for the stopping criterion.

GMRES algorithm for the solution of the interface problem

1. Compute the residual r = d   Fk0, and compute b = krk.

2. First basis vector: v1 = r/b.

3. For j = 1,...,m.

3.1. Compute wj = Fvj.

3.2. Comput0e Gram–Schmidt coefficients: hij ¼ vH

i wj; i ¼ 1; ... ;j.

3.3. Gram–Schmidt orthogonalization: wj ¼ wj   Pj

i¼1hijvi.

3.4. Normalization: vj+1 = wj/hj+1,j with hj+1,j = kwjk.

3.5. Let Hj ¼ ½hil

ðjþ1;jÞ

ði;jÞ¼ð1;1Þ with hil = 0 for i > l + 1.

3.6. Solve the least squares problem Hjzj = be1.

3.7. Update the solution k = k0 + [v1,..., vj]zj

4. Compute r = d   Fk.

5. If krk > TOL, set k0 = k, and go to step 1.

The parallel solution of the linear system (Fk = d) takes place as follows. The vectors k, d as well as the

iteration vectors are distributed in the same way among the processors of the distributed computer. The oper￾ations on the vectors of large dimension are all carried out in parallel. This does not require communication

between processors for vector updates (y y + ax) since this operation can be carried out for each compo￾nent of y independently, but communication is required for the inner product f = xHy. The other operations

are carried out simultaneously without communication on each processor.

4.3. Case of (semi-)infinite interface

The block form of the linear system (4.1) is very similar to the linear system (3.2), and can be assimilated to

the case of an interface between a sub-domain composed with finite elements only, with a sub-domain com￾posed with infinite elements only. When a general mesh partitioning of the global domain is performed, the

interface joins some (finite or infinite) elements sharing a common edge on the interface and belonging to dif￾ferent sub-domains. Three possibilities may appear: two finite elements sharing an edge on the interface, or

one finite element and one infinite element sharing an edge on the interface, or two infinite elements sharing

an edge on the interface. In this last case the length of the interface is infinite.

If some Lagrange finite elements are considered, for example P1-finite elements, the degrees of freedom of

an element corresponds to the nodes of the triangle. Defining the Lagrange multipliers at the nodes of the

finite element helps to apply the sub-structuring methodology described in (4.1). Fig. 3 shows the definition

of the degrees of freedom and of the Lagrange multipliers for two finite elements sharing one edge on the

interface.

In the second case, the Lagrange multipliers should be defined at the element nodes as shown in Fig. 4.

Indeed in this case, the restriction on the edge of the angular basis functions of the infinite element is similar

to the restriction of the P1-finite element basis functions. As a consequence, there is no difference between this

case and the previous one.

In the third case, the Lagrange multipliers should be defined at the element nodes and at the Gauss points

of the infinite elements as shown in Fig. 5. These Gauss points correspond to the degree of freedom of the

infinite element and are used to compute the integrals. Increasing the order of the infinite element implies

red in the non-overlapping Schwarz algorithm, a surface mass matrix should

be computed on the interface between the sub-domain. This matrix is of the form

MC ¼

Z

C

uvdS.

In the case of an interface between two finite elements, the coefficients of the matrix MC are computed as

½MC

lm ¼

Z

Xð1Þ

\Xð2Þ

NlN m dS;

where Nl and Nm are the finite element shape functions associated with node l and node m on the common

edge on the interface between sub-domains X(1) and X(2). In the case of an interface between one finite element

and one infinite element, the finite element shape functions on the common edge is similar to the angular infi￾nite element shape function i.e. the functions Nm. When two infinite elements share a common edge, the inte￾gral along this infinite edge only involves the radial shape functions i.e. the functions Nl, and the integral is

computed using the Gauss points along the infinite edge.

5. Implementation in SYSNOISE

SYSNOISE is an acoustic software developed by LMS International [19]. SYSNOISE predicts the radia￾tion, scattering and transmission of sound waves and the structural vibrations induced by the loading effects

of an acoustic fluid onto a structure. The program calculates a wide variety of results such as sound pressure

and radiated sound power, acoustic velocities and intensities, contributions of panel groups of the sounds,

energy densities, vibro-acoustic sensitivities, normal modes and structural deflections.

SYSNOISE utilizes state-of-the-art numerical methods based on direct and indirect boundary element

method (DBEM and IBEM) and a pressure formulation for acoustic finite and infinite element modelling

(FEM and I-FEM). Finite element methods are well adapted for enclosed regions, such as passenger compart￾ments, ducts or covers [27,28]. They are used to calculate the vibro-acoustic response of the cavity due to a

defined excitation, in the time or frequency domain, with the possibility to take into account flow effects. A

finite element model represents the elasticity of the fluid-loaded structure. SYSNOISE uses a complementary

infinite element method (I-FEM) to predict the vibro-acoustic response and sensitivities of a vibrating struc￾ture for an exterior acoustic domain. The method also solves coupled fluid-structure problems, and is well sui￾ted to multifluid applications.

The SYSNOISE DDT, implemented in the commercial version of SYSNOISE Revision 5.4 [19], is related

to the implementation of the optimized Schwarz algorithm into the harmonic acoustic FEM and I-FEM mod￾ule. The implemented domain decomposition technique (DDT) aims at partitioning large model, using a geo￾metric or an automatic mesh partitioning like METIS [29,30]. The elements of the mesh are connected via their

Fig. 5. Degree of freedom (white and black bullets) of the elements and of the Lagrange multiplier (black bullets) between two infinite

faces, which are surfaces in 3D. The elements of the global domain are decomposed into several sub-domains

by numbering or coloring all elements. All elements with the same number or color form a sub-domain. Sub￾domains are sets of elements, so the interface consists of a set of faces. The interface is defined by the faces

connected to elements with a different sub-domain number. The mesh partitioning should be done in such

a way to obtain a good load balancing among the sub-domains, and so that the number of interface nodes

is small in order to have a small interface problem. In many cases, this decomposition can be done by hand,

but for practical applications, it is often a very difficult and tedious task.

The non-overlapping Schwarz algorithm is by nature parallel: the solution can be computed independently

for each sub-domain, while the interface equation connects all the independent sub-domains into a global

problem. Because of this independence, each domain can be allocated to a single processor of the parallel sys￾tem, for example. Operations related to a single sub-domain take place without communication. The solution

of the interface problem, however, requires communication since it connects all the sub-domains. The sub￾domains are then handled by different processes which can run on different computer. The available non-over￾lapping Schwarz algorithm involves the direct solution method in each sub-domain, and an iterative solution

method for solving the interface problem. Each process uses the standard SYSNOISE procedure to solve its

own acoustic problem and communications between the process are ensured with the message passing inter￾face library (MPI). These communications are mandatory to solve the interface problem. The purpose of a

good load balancing among the sub-domains allows to provide a load balance between the various processors

allocated to the numerical treatment. The numerical benefits derive from the local bandwidth reduction (smal￾ler than the whole model bandwidth) that creates a gain both in memory space and computing time.

Since each iteration of the algorithm requires the solution of a linear system with a local impedance matrix,

it is advantageous to factorize once and perform just the forward and backward substitutions. The local

assembled finite element matrix was factorized and the local linear system was solved by the SYSNOISE

built-in direct solver, which is based on an LDLt factorization. Two iterative solvers are available for the solu￾tion of the interface problem. The BiCGstab is based on the conjugate gradient method, the GMRES (general￾ized minimum residual) is designed to solve non-symmetric linear system. In these algorithms, the

computation of vector inner products is performed by global communication commands from the MPI library

(MPI_Allreduce). The current implementation does not overlap communication and computation. Non￾blocking communication was used for the computation of the Lagrange multipliers. The communication con￾sists of a sequence of two-processor exchanges, that represent two neighboring sub-domains. These are the

MPI commands MPI_Isend, and MPI_Irecv. The other operations of the algorithms are fully parallel without

communication.

6. Numerical experiments

In this section the performance of the non-overlapping Schwarz method is evaluated for the solution of

acoustical problems arising from infinite and finite element discretization in the frequency domain. Two main

test cases have been studied for the evaluation of the Schwarz method. The first model considered here is a

simplified engine that has been used within the frame of a European project (PIANO). The second model

is an engine in free field conditions which has been used within the framework of a European project (DOM￾INOS). The objective of these two test cases is to evaluate the sound radiation generated by a normal velocity

boundary condition defined along the surface of the object. This normal velocity boundary condition is

obtained from measured vibration or calculated from a finite element analysis (FEA). From these boundary

conditions, SYSNOISE determines the radiated sound on the body surface as well as anywhere in the acoustic

field.

6.1. Simplified engine

This first example is related to a simplified engine. The objective is to evaluate the sound radiation with a

normal velocity boundary condition along the engine surface. Radiation problems are exterior problems that

can be handled quite efficiently using an infinite element model [11]. The current implementation of the non￾overlapping Schwarz algorithm has been done in such a way that infinite elements can be handled within each

6.1. Simplified engine

This first example is related to a simplified engine. The objective is to evaluate the sound radiation with a

normal velocity boundary condition along the engine surface. Radiation problems are exterior problems that

can be handled quite efficiently using an infinite element model [11]. The current implementation of the non￾overlapping Schwarz algorithm has been done in such a way that infinite elements can be handled within each

sub-domain together with conventional finite elements. The conjugate infinite element formulation of first

order is here considered [3,10].

The global mesh and the mesh partitioning into three sub-domains are illustrated in Fig. 6. The non-over￾lapping Schwarz method with absorbing boundary conditions is considered, and the interface problem is

solved iteratively with the GMRES algorithm. The stopping criterion is krnk2 < 1066

kr0k2, where rn and r0 are

the nth and the initial global residuals, and where kÆk2 denotes the L2

-norm. All simulations were performed

on an SGI ORIGIN 2000 machine with eight CPUs. Each sub-domain is allocated to a different processor,

and data exchange between the processors are performed with the MPI library.

Table 1 presents the convergence results of the Schwarz algorithm for different frequencies and different

number of sub-domains for the simplified engine. This table clearly outlines the robustness of the non-over￾lapping Schwarz method with absorbing boundary conditions. The improvement of the Schwarz method is

clearly illustrated when using an optimized procedure (OO0) on the interface rather than an expansion pro￾cedure (TO0). It is important to point out that one iteration of TO0 involves exactly the same number of oper￾ations and exactly the same computational than one iteration of the OO0 method [22]. For this reason, the

CPU time is not indicated in Table 1.

6.2. Engine in free field conditions

Large applications in automotive acoustic simulations are concerned with cavity problem, like passenger

compartments [27,28] with upper frequencies higher than 200 Hz. On the part of the engine compartment

The order of the infinite element is equal to m = 3.

exterior acoustic radiation of engine and transmission up to 2.5 kHz are of interest. The simulation of the

acoustic radiation of an engine is necessary to predict the impact on the passenger compartment. For compar￾ison to vehicle testing it is more convenient to examine the acoustic radiation within free field boundary con￾ditions. In FEM analysis, a free field boundary condition can be taken into account by means of infinite

elements. For this purpose a single layer of infinite elements of variable order is matched onto an ellipsoidal

surface of a conventional finite element mesh. The interior boundary of the acoustic model is represented by

the surface of the engine. The acoustic excitation, caused by the engine vibrations is calculated in a separate

FE model by means of MSC/NASTRAN and is imported into the acoustic model as frequency dependant

boundary condition on the engine surface.

The ellipsoidal air volume surrounding the engine has been modelled by 250 000 FEM-elements and 30 000

IFEM-elements, and 54 000 nodes. Figs. 7 and 8 illustrate the geometry of the engine and the mesh partition￾ing of the ellipsoidal air volume into four sub-domains.

The non-overlapping Schwarz algorithm is considered for the solution of the acoustical radiation problem.

The GMRES algorithm is considered in the following analysis. As already discussed, the speed of the conver￾gence of the GMRES algorithm is linked with the maximum number of stored descent direction vectors. This

number is associated to a window of iterations. Within one window each descent direction vector used to

search the solution is kept in memory. Unfortunately, when the number of iterations becomes larger than

the maximum number of iteration vectors, the iteration method is restarted. Any restart induces a deteriora￾tion in the number of iterations needed to reach convergence since all iteration vectors from the preceding win￾dow are lost. To achieve fast convergence, the maximum number of iteration vectors should be chosen large

enough to perform all iterations within one window.

Table 2 presents the number of iterations for different frequencies and different number of sub-domains for

the engine in free field conditions. Once again, the non-overlapping Schwarz method with OO0 conditions

converges faster than with TO0 conditions. In Table 2, it can be noticed that increasing the frequency increases

the number of iterations. This property is due to the fact that at higher frequencies the acoustic solution

becomes more complex and therefore deteriorate the conditioning number of the linear system defined on

the interface. As already noticed for internal acoustic problems [27], increasing the number of sub-domains

increases the number of iterations. Anyway, this does not implies additional CPU time, since the supplemen￾tary number of iteration is overcome by smaller sub-domains which are solved by a direct solution method.

7. Conclusion

In this paper, an infinite element method has been presented and analyzed for parallel computations pur￾pose. This method is interesting for solving acoustical radiation problems in unbounded domain. A finite ele￾ments domain decomposition method has then been reminded. It has been shown, that an original

formulation of this domain decomposition method with Lagrange multipliers defined on (semi-)infinite space

allows to extend this method to infinite elements. This technique has been implemented in an industrial acoustic software (SYSNOISE). The numerical performance of this software have been illustrated on two

problems in the field of the engine development process. Some enhancements of the current implementation

should be done with respect to the scalability upon the number of sub-domains and with the frequency

parameter.

Acknowledgements

The authors acknowledge partial support by the European Commission under Grant ESPRIT-25009, and

are grateful to DaimlerChrysler and to LMS International, partners in this project, for providing the test cases

and equipment for running the numerical examples. The authors would like to acknowledge A. Kongeter, }

J.-P. Coyette, C. Meir and J.-L. Migeot for the useful discussions, comments and remarks.