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 computations 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 elements 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 proportionally 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 memory 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 conjugate 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 computational 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 presented 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 artificial 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 boundary 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-overlapping 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 Steklov–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 minimization 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 iteration 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 algorithm, by 3.1–3.7 requires one matrix-vector product. The quantity xH denotes the complex conjugate transpose 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 operations 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 component 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 composed 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 different 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 infinite element shape function i.e. the functions Nm. When two infinite elements share a common edge, the integral 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 radiation, 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 compartments, 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 structure for an exterior acoustic domain. The method also solves coupled fluid-structure problems, and is well suited 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 module. The implemented domain decomposition technique (DDT) aims at partitioning large model, using a geometric 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. Subdomains 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 system, 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 subdomains are then handled by different processes which can run on different computer. The available non-overlapping 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 interface 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 (smaller 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 solution of the interface problem. The BiCGstab is based on the conjugate gradient method, the GMRES (generalized 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. Nonblocking communication was used for the computation of the Lagrange multipliers. The communication consists 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 (DOMINOS). 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 nonoverlapping 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 nonoverlapping 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-overlapping 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-overlapping 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 procedure (TO0). It is important to point out that one iteration of TO0 involves exactly the same number of operations 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 comparison to vehicle testing it is more convenient to examine the acoustic radiation within free field boundary conditions. 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 partitioning 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 convergence 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 deterioration in the number of iterations needed to reach convergence since all iteration vectors from the preceding window 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 supplementary 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 purpose. This method is interesting for solving acoustical radiation problems in unbounded domain. A finite elements 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.
