Low-Frequency Acoustic-Structure Analysis Using Coupled FEM-BEM Method

Abstract

A numerical algorithm based on finite element method (FEM) and boundary element method (BEM) is proposed for the analyses of acoustic-structure coupled problems. By this algorithm, the structural domain and the acoustic domain are modelled by FEM and BEM, respectively, which are coupled with each other through the consideration of the appropriate compatibility and equilibrium conditions on the interface of the two domains. To improve the computational efficiency, the adaptive cross approximation (ACA) approach is incorporated into the proposed algorithm to deal with the nonsymmetric and fully populated matrices resulting from the coupling of the FEM and BEM. The validity and the high efficiency of the present approach are demonstrated by two examples.

1. Introduction

Acoustic-structure coupling problems can be found in various engineering fields. For example, in space exploration, spacecrafts are subject to heavy acoustic load, particularly during launching, which can impose severe and adverse effect on the astronauts and the structure of the spacecrafts. Therefore, it is necessary to develop numerical methods to predict the structural behaviour in the acoustic-structure coupling problems since these numerical techniques are essential to the improvement of practical engineering design.

The finite element method (FEM) and boundary element method (BEM) are both widely used in various engineering problems. However, FEM and BEM are generally applied in different areas. For instance, the former is especially well suited to the analyses of the problems which involve inhomogeneity or nonlinearity [1, 2], whilst the latter is one of the most powerful and suitable numerical techniques for predicting the noise around a vibrating system [3]. Compared to FEM [4], BEM reduces the mathematical dimension of the problem under analysis by one. Furthermore, in BEM, the Sommerfeld condition is automatically satisfied so that the external domain does not need to be bounded. As a result, BEM is regarded as the most suitable tool to deal with the problems in exterior unbounded domains.

Therefore, for acoustic-structure problems involving submerged bodies, it should be a natural strategy to combine FEM and BEM in a computational model to make use of the respective advantages of the two methods. The combination of FEM and BEM can be implemented by decomposing the concerned domain into several subdomains, each of which, according to its specific physics, is modelled by either FEM or BEM [5–7]. Because of the appropriate compatibility and equilibrium conditions on the respective interface boundaries, the subdomains are coupled to each other.

It should be mentioned that the boundary integral equation has a major defect for exterior problems; that is, it has nonunique solution at a set of fictitious eigenfrequencies associated with the resonant frequencies. In the past, several methods have been proposed to overcome this difficulty. The two most popular methods are the CHIEF method and Burton-Miller approach. The CHIEF method [8] is based on the Helmholtz integral equation on the surface of the radiating body combined with Helmholtz equations for some interior points. The resulting system equation is usually solved by the least-square approach. The Burton-Miller approach [9], inspired by Panich formulation [10], forms a linear combination of the Helmholtz boundary integral equation and its normal derivative, providing also a valid solution at any frequency but leading to hypersingular integrals which can be handled by Guiggiani algorithm [11, 12]. In this study, the Burton-Miller approach is adopted.

Similar to BEM, the coupled FEM-BEM approach leads to a fully populated system matrix. This is thought to be a major drawback of the coupled approach as compared with FEM which normally results in a sparse symmetrical system matrix. If iterative methods like the generalized minimal residual (GMRES) method [13] in combination with suitable preconditioners are used to solve the FEM-BEM system with 𝑁 degrees of freedom, the computational expense is in the order of 𝑂(𝑁²). If a direct solver is applied, the computational cost is even in the order of 𝑂(𝑁³). Fortunately, over the last few years, several fast solution methods including fast multipole method (FMM) [14–16] and adaptive cross approximation approach [17, 18] have been developed. These methods reduce the memory requirements and the computational time significantly when dealing with the matrix-vector product. In this study, the ACA approach is adopted since it is much easier to operate compared with FMM which is based on a series expansion of kernel-shape function products.

In this paper, a coupling algorithm combining FEM and BEM is developed for the analysis of acoustic-structure response. This paper is organized as follows. Section 1 is a short introduction to the background and corresponding research status. In Section 2, the FE formulations of the structure domain are constructed. In Section 3, the application of the BEM to Helmholtz problem and its implementation is presented. Then, the coupled equations are built in Section 4. In Section 5, the solving algorithms of coupled equations including the GMRES iterative solver and ACA approach are introduced. In Section 6, several numerical examples are presented.

2. FE Formulation of Structural Domain

A coupling system of acoustic-structure is shown in Figure 1. The structure is assumed to be fully submerged in interior and exterior acoustic fields or one of them. In Figure 1, the structure domain is denoted by Ω^𝑆. The interior and the exterior acoustic fields are denoted by Ω^𝐼 and Ω^𝐸, respectively. The inner interface and the outer interface are named as 𝑆^𝐼 and 𝑆^𝐸, respectively. The structure is assumed to be thin and flexible, which is susceptible to acoustic-structure interaction.

Figure 1 

Illustration of a coupling system of acoustic-structure.

FEM is predominantly chosen for simulating linear elastodynamic systems. When the structure Ω^𝑆 is discretized by FEM, the resulting system of linear equations reads as 𝐌𝐮̈+𝐂𝐮̇+𝐊𝐮=𝐟_𝑠+𝐓_𝑠, (1) where 𝐌, 𝐂, and 𝐊 denote the global mass matrix, the damping matrix, and the stiffness matrix, respectively. The vector 𝐟_𝑠 is load vector, while the vector 𝐓_𝑠 resulted from the tractions on the coupling interface, and it will be further discussed in the following section. The global displacement vector is denoted by 𝐮.

In this paper, Rayleigh damping is considered for the damping matrix, which means that 𝐂=𝜉𝐌+𝜂𝐊, (2) where 𝜉 and 𝜂 are the damping parameters.

If the load vector 𝐟_𝑠 and the acoustic pressure vector 𝐓_𝑠 are both assumed to be harmonic, the displacement vector 𝐮 should be also harmonic. This means that 𝐟_𝑠=𝐟_𝑠𝑒^𝑖𝜔𝑡,𝐓_𝑠=𝐓_𝑠𝑒^𝑖𝜔𝑡,𝐮=𝐮𝑒^𝑖𝜔𝑡, (3) where the imaginary unit is denoted by 𝑖 and 𝜔=2𝜋𝑓 is the circular frequency with the excitation frequency 𝑓. For the relationship between 𝐮 and 𝐮̇, one obtains 𝐮̇(𝑡)=𝑖𝜔𝐮𝑒^𝑖𝜔𝑡=𝐮̇(𝜔)𝑒^𝑖𝜔𝑡,𝐮̈(𝑡)=−𝜔²𝐮𝑒^𝑖𝜔𝑡=𝑖𝜔𝐮̇(𝜔)𝑒^𝑖𝜔𝑡. (4)

Substituting (4) into (1) yields 𝑖𝜔𝐌+𝐂−𝑖𝐊𝜔𝐮̇(𝜔)=𝐟_𝑠+𝐓_𝑠. (5)

NASTRAN, a finite element package, is employed to set up 𝐌, 𝐂, 𝐊, and the right-hand side vector 𝐟_𝑠. Since 𝐌, 𝐂, and 𝐊 are frequency independent, these matrices only have to be calculated once for a given model.

3. The Boundary Element Method for the Helmholtz Equation

Acoustic problems in frequency domain can be described by the Helmholtz equation which has a three-dimensional form given by ∇²𝑝(𝐱)+𝑘²𝑝(𝐱)=0, (6) where ∇² is the Laplace operator, 𝑝(𝐱) is the sound pressure at point 𝐱 in the fluid,𝑘=𝜔/𝑐 is the acoustic wavenumber, and 𝑐 is the speed of sound in the fluid.

The boundary condition can be expressed as 𝑝(𝐱)=𝑝(𝐱),𝐱∈𝑆^𝐷,𝜕𝑝(𝐱)=𝜕𝑛𝑞(𝐱),𝐱∈𝑆^𝑁. (7)

Additionally, the three-dimensional Sommerfeld radiation condition [19] is 𝜕𝑝≤𝑐𝜕𝑟−𝑖𝑘𝑝𝑟²,𝑟⟶∞. (8) For the exterior problem, (8) has to be satisfied to ensure that the wave is purely outgoing.

A three-dimensional fundamental solution for the Helmholtz equation is given by 𝐺^∗()=𝑒𝐱,𝐲^𝑖𝑘𝑟,4𝜋𝑟 (9) where 𝐱 is a field point, 𝐲 is a source point, and 𝑟=|𝐱−𝐲| is the distance from 𝐱 to 𝐲. The key idea is to use Green’s second identity in combination with the property of the Dirac distribution: _{Ω𝐸∪Ω𝐼}𝑝∇²𝐺^∗−𝐺^∗∇²𝑝=d𝑉_{𝑆𝐸∪𝑆𝐼}𝑝𝜕𝐺^∗𝜕𝑛−𝐺^∗𝜕𝑝𝜕𝑛d𝑆. (10)

By this way, the pressure 𝑝 at an arbitrary point 𝐱 within the acoustic domain Ω^𝐸∪Ω^𝐼 is given by the integral representation 𝑝(𝐱)=−_{𝑆𝐸∪𝑆𝐼}𝐺^∗(𝐱,𝐲)𝜕𝑝(𝐲)𝜕𝑛d𝑆+_{𝑆𝐸∪𝑆𝐼}𝜕𝐺^∗(𝐱,𝐲)𝜕𝑛_𝑦𝑝(𝐲)d𝑆,𝐱∈Ω^𝐸∪Ω^𝐼. (11)

For scattering problems (11) can be written as 𝑝(𝐱)=−_{𝑆𝐸∪𝑆𝐼}𝐺^∗(𝐱,𝐲)𝜕𝑝(𝐲)𝜕𝑛d𝑆+_{𝑆𝐸∪𝑆𝐼}𝜕𝐺^∗(𝐱,𝐲)𝜕𝑛_𝑦𝑝(𝐲)d𝑆+𝑝^𝐼(𝐱),𝐱∈Ω^𝐸∪Ω^𝐼. (12)

Let the point 𝐱 approach the boundary; then the following conventional boundary integral equation (CBIE) [20] is obtained: 𝐶(𝐱)𝑝(𝐱)=_𝑆𝐺^∗𝜕𝑝𝜕𝑛−𝑝𝜕𝐺^∗𝜕𝑛d𝑆, (13) where 𝐶(𝐱)=1/2  if 𝑆 is smooth around 𝐱. This CBIE can be employed to solve the unknown 𝑝 and 𝜕𝑝/𝜕𝑛 on 𝑆^𝐸∪𝑆^𝐼.

It is well known that the CBIE has a major defect for exterior problems; that is, it has nonunique solution at a set of fictitious eigenfrequencies associated with the resonant frequencies [21]. A remedy to this problem is to use the normal derivative BIE in conjunction with the CBIE. Taking the derivative of the integral representation given in (13) with respect to the normal at a typical point 𝐱 on 𝑆^𝐸∪𝑆^𝐼 and letting 𝐱 approach 𝑆^𝐸∪𝑆^𝐼, one obtains the following hypersingular boundary integral equation (HBIE): 𝐶^′(𝐱)𝑞(𝐱)=−_{𝑆𝐸∪𝑆𝐼}𝜕𝐺^∗(𝐱,𝐲)𝜕𝑛_𝑥𝜕𝑝(𝐲)+𝜕𝑛d𝑆_{𝑆𝐸∪𝑆𝐼}𝜕²𝐺^∗(𝐱,𝐲)𝜕𝑛_𝑥𝜕𝑛_𝑦𝑝(𝐲)d𝑆,𝐱∈𝑆^𝐸∪𝑆^𝐼, (14) where 𝐶^′(𝐱)=1/2 if 𝑆^𝐸∪𝑆^𝐼 is smooth around 𝐱. It should be noted that (14) suffers from the same kind of defect as (13). However, if a linear combination of CBIE and HBIE is used, the uniqueness of the results can be ensured for exterior acoustic wave problems. Therefore, _{𝑆𝐸∪𝑆𝐼}𝜕𝐺^∗𝜕𝑛_𝑦𝑝(𝐲)d𝑆+𝐶(𝐱)𝑝(𝐱)+𝛼_{𝑆𝐸∪𝑆𝐼}𝜕²𝐺^∗(𝐱)(𝐲)𝑝(𝐲)𝜕𝑛𝜕𝑛d𝑆=_{𝑆𝐸∪𝑆𝐼}𝐺^∗𝜕𝑝𝜕𝑛d𝑆+𝛼_{𝑆𝐸∪𝑆𝐼}𝜕𝐺^∗𝜕𝑛_𝑥𝜕𝑝𝜕𝑛d𝑆−𝐶^′(𝐱)𝜕𝑝𝜕𝑛_𝑥(𝐱), (15) where 𝛼 is the coupling constant. This formulation is called Burton-Miller formulation [9] for acoustic wave problems and has been shown to yield unique solutions at all frequencies, if 𝛼 is a complex number which, for example, can be chosen as 𝛼=𝑖/𝜅 with 𝑖=√−1 [21]. In (15), 𝜕𝑝/𝜕𝑛=−𝑖𝜔𝜌𝐯, and 𝜌 and 𝐯 denote the density of the fluid and the normal velocity on the interface, respectively.

When the boundary of the acoustic domain is discretized into boundary elements, the resulting linear system equations can be expressed as [𝐇][𝐩]=[𝐆][𝐯], (16) where the matrices [𝐇] and [𝐆] are obtained by integrating the fundamental solutions over each boundary element.

Equation (16) is the boundary element equation used to solve the acoustical problem. When there is incident sound, a free term is needed, and (16) becomes [𝐇][𝐩]=[𝐆][𝐯]+𝐩_𝐼. (17)

4. Coupled FE-BE Formulations

By taking into account the appropriate compatibility and equilibrium conditions at the respective interface boundaries, the fully coupled FE-BE formulations can be derived. Hereby, the matching grids for the FE and the BE parts are required. The compatibility condition at the respective interface boundaries links the acoustic pressure 𝑝 and the tractions 𝑇_𝑛 on the structure, which can be written in the following matrix form: 𝑇_𝑛=𝐓⋅𝐧=−𝑝. (18)

For an element, one has 𝐓^𝑒_𝑠=_𝑆𝜎𝐍^𝑇𝐓dΓ=_𝑆𝜎𝐍^𝑇𝑇_𝑛𝐧dΓ=−_𝑆𝜎𝐍^𝑇𝐧𝐍^𝑏𝑇dΓ𝐩=−𝐀^s_𝑒𝐩, (19) where 𝐀^s_𝑒=∫_𝑆𝜎𝐍^𝑇𝐧[𝐍^𝑏]^𝑇dΓ.

Substituting (19) into (5) and using uniform symbol for the matrices and vectors yield [𝐌]+[𝐂]−𝑖[𝐊]𝑖𝜔𝜔𝐮̇(𝜔)=𝐟_𝑠−𝐀^𝑠[𝐩]. (20)

From (20), one obtains 𝐮̇(𝜔)=𝑖𝜔[𝐌]+[𝐂]−𝑖[𝐊]𝜔⁻¹𝐟_𝑠−𝐀^𝑠[𝐩]. (21)

Please note that [𝐮̇(𝜔)] in (21) incorporates the velocity vector of all of the nodes in the structure, and [𝐯] in (17) just denotes the velocity vector of the nodes on the surface of the structures. The relationship between them reads as [𝐯]=[𝐑]𝐮̇(𝜔), (22) where [𝐑] is a transformation matrix from the velocity vector for all nodes in the structure to the velocity vector for the surface nodes [22].

Substituting (21) into (22) yields [𝐯]=[𝐑[𝐌]+[𝐂]−𝑖[𝐊]]𝑖𝜔𝜔⁻¹𝐟_𝑠−𝐀^𝑠[𝐩]. (23)

Further more, substituting (23) into (17) yields 𝐇]+[𝐆𝐑[𝐌]+[𝐂]−𝑖[𝐊][][]𝑖𝜔𝜔⁻¹𝐀^𝑠𝐩]=[𝐆𝐑[𝐌]+[𝐂]−𝑖[𝐊][][]𝑖𝜔𝜔⁻¹𝐟_𝑠+𝐩_𝐼. (24)

Equation (24) is the coupled FEM-BEM equation, from which the sound pressure of the nodes on the interfaces can be obtained.

5. Solving Algorithm of the Coupled FEM-BEM Equation

Because of the nonsymmetric and fully populated matrices in the coupled FEM-BEM equations, the ACA approach is used in this study to reduce the memory requirements and the computation cost. During iteration process, the generalized minimal residual (GMRES) method iterative solver is adopted. In the following, the ACA and the GMRES algorithms are introduced briefly. More details can be found in relevant literatures [17, 18, 23].

5.1. Adaptive Cross Approximation

The ACA algorithm produced by M. Bebendorf and S. Rjasanow is an effective technique for solving nonsymmetric and fully populated matrices. The main process is as follows.(1)Index Octree. The low-rank approximation of the matrix of FEM-BEM is based on asymptotically smooth functions, which happens in a well-separated domain. Thus, a tree structure for all boundary points should be constructed at first. For a three-dimensional domain, an octree is necessary to be built to describe the geometry relationship of sections on surface. The information that a tree should record includes its father, its children, its geometry information, and the points contained. The least nodes are called leaves whose numbers of boundary points contained are less than a given threshold value used to judge whether a new division is necessary.(2)Partitioning of Matrix. Based on the tree, the matrix in the FEM-BEM equation can be divided into some blocks. At first, the relationship between nodes should be established. For some node NP in an octree, its neighbour nodes are defined as the nodes which have the same geometry size and one common point at least. The interaction nodes are those nodes owing the same geometry size, whose fathers are neighbour nodes, but not for them. The relationship is shown as given Figure 2 in a two-dimensional domain. Then for each node of the tree, the elements contained are taken as source points and the elements of its neighbour as field points, and the submatrices coming from the source points and field points are defined as 𝐴_𝑛. The submatrices whose entries come from elements of the node and its interaction nodes are defined as 𝐴_𝑓. All of the submatrices 𝐴_𝑛 and 𝐴_𝑓 form the whole efficient matrix of FEM-BEM equation.(3)Low-Rank Approximation. The low-rank approximation in ACA is defined as 𝐀^𝑚×𝑛=𝐔^𝑚×𝑟𝐕^𝑟×𝑛=^𝑟_𝑖=1𝐮_𝑖^𝑚×1𝐯_𝑖^1×𝑛. (25) The target of ACA is to achieve 𝐑^𝑚×𝑛=𝐀^𝑚×𝑛−𝐀^𝑚×𝑛≤𝜀𝐀^𝑚×𝑛, (26) where 𝜀 is the tolerance of ACA.

Figure 2 

Relationship between nodes.

The ACA algorithm has been applied to the low-rank blocks achieving approximately 𝑂(𝑁) for both storage and matrix-vector multiplication [24]. It must be noted that the fully pivoted approach is well known to be much slower than the partially pivoted approach. The main reason behind this is that the fully pivoted approach requires the knowledge of the full matrix, whereas the partially pivoted approach would only require generation of individual matrix entries.

5.2. Generalized Minimal Residual Method

The generalized minimal residual (GMRES) method is one of the most popular iterative solvers for nonsymmetric linear systems. It was proposed by Saad and Schultz [25] and further developed by other researchers [26–28]. It has the property of minimizing at every step the norm of the residual vector over a Krylov subspace. The algorithm is derived from the Arnoldi process for constructing an orthogonal basis of Krylov subspace. Because of existence of error, the gained vectors lose orthogonality gradually. Therefore, preconditioner is needed. In the computer codes for this study, the solution of the diagonal blocks is taken as the initial solution, and the maximum number of iteration is set to 500.

6. Numerical Examples

In this section, the proposed coupling approach is applied to two examples for the simulation of acoustic-structure coupling problems.

6.1. Benchmark Example

The spherical test structure is an elastic sphere shell with a radius of 5 m and thickness of 0.15 m, full of some kind of liquid. A point sound source is located in the center of the sphere emitting harmonic spherical wave. The outer surface of the sphere is free, and the inner surface of the sphere is coupled with the acoustic field. The detailed geometrical and material parameters for the sphere and the fluid are shown in Table 1. In this example, the damping is ignored.

Table 1 

Characteristics of the fluid and the sphere.

This problem is simulated by using FEM and BEM, respectively. The BEM results are obtained using the proposed algorithm. The FEM results are obtained using NASTRAN software. The sound pressure amplitudes obtained by the two methods are plotted in Figure 3. From the two curves, we can see that the results using BEM agree well with the results obtained from FEM. However, from Figure 4 in which computational efficiencies of different solving methods are compared, it can be seen that the ACA approach requires the least computing time.

Figure 3 

Sound pressure curves obtained by BEM and FEM.

Figure 4 

Comparison of computing time.

6.2. Analysis of the Frequency Response of a Rectangular Box

Further more, the frequency response analysis of a square box is also presented. The rectangular box is a closure chamber made of thin plates. Its inner and outer surfaces of the box are coupled with the interior and the exterior acoustic fields, respectively. In this example, the exterior acoustic field is assumed to be infinite. A point acoustic source located outside the box initiates the coupled response of the structure and its surrounding acoustic fields. The mesh of the box is shown in Figure 5.

Figure 5 

Mesh of the square box.

Parameters of the box and the fluid are listed in Table 2.

Table 2 

Characteristics of the fluid and the structure.

To model this problem, the origin of the Cartesian coordinate system is set to coincide with the center of the bottom surface of the box, while the coordinate axes are parallel with the edges of the box. The box’s bottom surface is restrained. The point acoustic source generating harmonic spherical wave (𝑝₀=1Pa) is located at the point (4.0 m, 0.0 m, 4.0 m). Damping is ignored.

Two cases are considered. In the first case, the interior acoustic field is neglected, whilst in the second case the interior and the exterior acoustic fields are both considered. In the two cases, the sound pressures at a selected point (1.18 m, 0.0 m, 2.057 m) are calculated. They are compared in Figure 6.

Figure 6 

Sound pressure under different frequencies.

From Figure 6, it can be seen that the trends of the two pressure-frequency curves are similar, but the pressures are a bit different at some specific frequencies; for instance, the frequency is around 20. Figure 7 demonstrates the acoustic pressure on the outer surface of the box. According to the graphs, it is found that, when the frequency varies from 20 to 100, the pressure distributions on the outer surface of the box are almost the same for the two cases. Therefore, it can be concluded that, for this example, the neglection of the inner acoustic field in the modelling does not lead to ruinous result.

(a)

(b)

Figure 7 

(a) Sound pressure distribution for the coupled system without considering the inner field. (b) Sound pressure distribution for the coupled system considering the inner field.

7. Conclusions

In this study, an algorithm based on the coupling of FEM and BEM is developed for the analysis of acoustic-structure response. The FEM is employed to model the structure part, while the BEM is used to discretize the acoustic domain. The two domains are coupled with each other through the consideration of the appropriate compatibility and equilibrium conditions on the interface of the two domains. To improve the computational efficiency, the ACA approach is incorporated into the proposed algorithm to deal with the nonsymmetric and fully populated matrices resulting from the coupling of the FEM and BEM. By the proposed algorithm, the fictitious frequency problem can also be avoided. The validity and the high accuracy of the present algorithm are demonstrated by a benchmark example. The modeling for the frequency response of a rectangular box is also presented. The numerical examples show that the coupling method has great potential to deal with large-scale complex acoustic-structure problems.

Acknowledgment

This work is jointly supported by the National Natural Science Foundation of China (nos. 11272182 and 11072128).