Coupled BE–FE based vibroacoustic modal analysis and frequency sweep using a generalized resolvent sampling method

Accepted Manuscript Coupled BE-FE based vibroacoustic modal analysis and frequency sweep using a generalized resolvent sampling method Tengfei Liang, Junpeng Wang, Jinyou Xiao, Lihua Wen

PII: DOI: Reference:

S0045-7825(18)30491-2 https://doi.org/10.1016/j.cma.2018.09.038 CMA 12098

To appear in:

Comput. Methods Appl. Mech. Engrg.

Received date : 2 April 2018 Revised date : 19 September 2018 Accepted date : 24 September 2018 Please cite this article as: T. Liang, J. Wang, J. Xiao et al., Coupled BE-FE based vibroacoustic modal analysis and frequency sweep using a generalized resolvent sampling method, Computer Methods in Applied Mechanics and Engineering (2018), https://doi.org/10.1016/j.cma.2018.09.038 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.

*Manuscript Click here to download Manuscript: paper2_revision_mark_blue.pdf

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Click here to view linked References

Coupled BE-FE based vibroacoustic modal analysis and frequency sweep using a generalized resolvent sampling method Tengfei Lianga , Junpeng Wanga , Jinyou Xiaoa , Lihua Wena,∗ a School

of Astronautics, Northwestern Polytechnical University, Xi’an 710072, China

Abstract For several decades, the coupled boundary element and finite element (BE-FE) approach has been recognized as a potential tool for modeling sound-structure interaction problems. However, till now, the frequency response computation, frequency sweep and modal analysis, which are key ingredients for the characterisation of a vibro-acoustic system, are still challenging due to the lack of effective numerical methods to deal with the complicated coupled matrices. We contribute to this subject in the following three aspects. Firstly, a fast coupled method is established using the acoustic BEM accelerated by the kernel-independent directional algorithm and the conventional structural FEM. The method uses quadratic elements and allows for non-matching meshes across the interfaces, and thus is accurate and versatile in a wide frequency range. Secondly, based on the Keldysh theorem of analytic matrix functions, a generalized resolvent sampling method is developed to construct the response subspaces for given excitations in a certain frequency range. By using the resulting response subspaces, a projection-based fast frequency sweep algorithm and a nonlinear eigensolver for modal analysis have been developed. Lastly, a numerical scheme for fast computation of the reduced system matrices in the fast sweep algorithm and the eigensolver is proposed, and several rules for practical applications of the proposed algorithms are discussed. The computational performance of the proposed methods has been tested and confirmed by academic benchmarks and a large-scale industrial application. Keywords: sound-structure interaction, coupled boundary and finite element method, frequency sweep, modal analysis

1. Introduction In aerospace engineering, the simulation of structural-acoustic interaction of rocket fairings has a high priority due to the fact that the vibroacoustic environment inside the fairings of large launch vehicles can pose a significant risk to the payload launch survivability. Such considerable vibroacoustic problems also arise in submarine, automobile and many other industrial sectors. The coupled boundary element and finite element (BE-FE) method is a favorable scheme for solving vibroacoustic problems, as it combines the advantage of the finite element method (FEM) in structural vibration analysis and the boundary element method (BEM) in acoustic field analysis. Started from the pioneering work of Everstine and Henderson in 1990 [1], the development and applications of coupled BE-FE methods for vibroacoustic analysis have now become a broad area of research; see e.g., [2, 3] and the references therein. In recent years, there has been a renewed interest in developing the coupled FE-BE method in the framework of isogeometric analysis [4, 5, 6]. Frequency sweep and modal analysis are two fundamental ways to acquire the dynamic characteristics and parameters of vibroacoustic systems. The former evaluates the frequency response functions in a frequency ∗ Corresponding

author Email addresses: [email protected] (Tengfei Liang), [email protected] (Jinyou Xiao), [email protected] (Lihua Wen) Preprint submitted to Elsevier

September 19, 2018

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

band, and the latter extracts the natural frequencies and modes by solving a certain type of eigenvalue problems. From an academic point of view, the frequency sweep problem (FSP) and the eigenvalue problem (EVP) are two different topics studied separately. In the practical design of vibroacoustic systems, however, the solutions of both problems are often conducted reciprocally to fully understand the intrinsic dynamic properties as well as the response profile to excitations. This practical demand inspires us to consider the two problems in a unified framework, which finally leads to the present work. Generally, the straightforward solution of large-scale FSPs with fine frequency increment is computationally expensive, if not prohibitive. On the other hand, the straightforward solution is unnecessary in most case as the responses typically lie in a small subspace. Motivated by this fact, model order reduction (MOR) approaches have been proposed to alleviate the computational burden; see [7] and the review paper [8]. These approaches generally fall into two classes: the projection based MOR and the interpolation based MOR. Given a certain projection subspace Q, the former searches for the solutions in Q by using a Galerkin or Petrov Galerkin projection or a least-squares minimization. The latter is also often referred to as moment matching; it constructs a certain type of expansion of the response functions, e.g., Taylor or Pad´e expansion. For a comparison of these two classes of approaches and their connections in dealing with structural dynamics and acoustics problems, the readers are referred to [9]. Here we notice that most moment matching approaches using high-order derivatives would be troublesome in treating coupled BE-FE systems, because the BE matrices are typically dense, nonsymmetric and the entries are functions of the frequency without explicit expressions. As a consequence, the accurate evaluation of the derivatives of the response function [9] is time and memory consuming. Therefore, this paper will follow the projection based MOR and construct the projection subspace Q using only the values of the response function at a series of sampling points. Finally, it is worth mentioning that MOR for coupled structural-acoustic BE-FE systems has been exploited by Peters et al. [10, 7]. In particular, in [7], the frequency independent FE matrices of the global system of equations for the fully coupled problem were condensed by using a Krylov subspace approach, while the BE matrices were left unchanged. The complicated structure of BE matrices also precludes the use of some well recognized sparse EVP solvers, such as the Lanczos, Arnoldi or subspace iteration method [11], to solve the coupled EVPs. Existing numerical methods for solving the coupled EVPs can be summarized according to how the BE domains are modeled and treated. In some fluid-structure interaction problems, the fluid is considered to be incompressible and modeled by the Laplace boundary integral equations (BIEs) that are frequency-independent [12]. This leads to nonsymmetric generalized EVPs whose matrices have dense and frequency-independent BE sub-matrices. For these EVPs, iterative solvers such as the Arnoldi and nonsymmetric Lanczos methods mentioned in [13] can be used. When the fluids are compressible, the BE domains are often modeled by frequency-dependent BIEs, leading to BE sub-matrices that are complicated matrix-valued functions of the frequency. As a result, the coupled EVPs become totally nonlinear, and their large-scale solution is sofar a challenging task in both mathematical and engineering communities [14, 15, 16]. Note that compressible fluids can also be modeled by boundary-volume integral equations by arranging the frequency z into the “inertial” terms [17]. This will lead to generalized EVPs as well. Whereas, following this clue, the evaluation of the volume integrals and the proper treatment of the volume variables are still bottlenecks restricting the application to complicated engineering problems, and of course, the dense matrices pose further difficulties. Numerical solution of the dense nonlinear EVPs in BE and coupled FE-BE analyses has become a research topic for several decades [18]. Early methods generally fall into two classes: the Newton methods [19] and the linearization methods [20]. However, these methods are computationally expensive for largescale problems, due to the need for repeated formulation and solution of dense linear systems and/or the storage of the dense coefficient matrices of the BE matrix in a polynomial basis. In recent years, the contour integral approach has emerged as an efficient method for solving nonlinear EVPs [21, 22, 23]. This approach possesses two salient features that are appealing in BEM applications: (1) its applicability is irrelevant to the structure of the system matrices, and the computationally intensive part can be easily parallelized; (2) it can simultaneously compute all the eigenvalues (and the associated eigenvectors) within a given contour in the complex plane. It has been applied to solve nonlinear EVPs in the BE and coupled FE-BE analyses by several authors, e.g., [24, 25, 26, 27, 28]. However, we notice that the involvement of high order contour moments of the resolvent matrix in the contour integral approach might lead to a loss of accuracy of the 2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

eigen-solutions [23, 29]. To overcome this drawback, the resolvent sampling based Rayleigh-Ritz method, denoted by RSRR, has been proposed by the present authors’ group [29, 30, 31]. The RSRR algorithm inherits the aforementioned advantages of the contour integral approach and enjoys much better stability and accuracy; see [32] for large-scale BEM applications. The merits of the RSRR algorithm stem mainly from the novel resolvent sampling scheme (abbreviated as RESS) for constructing eigenspaces. The theory of RESS can be found in [29, 31]. RESS has a connection with the resolvent moment scheme based on the contour integral approach in [23]: at best, the moment scheme can produce the same eigenspace as RESS, but very often the moment scheme produces a smaller and more ill-conditioned subspace than RESS due to the involvement of the high order contour moments. RESS also has a connection with the NLEIGS solver based on the rational Krylov method in [33]: both of them construct eigenspaces using the value of the resolvent matrices at a series of sampling points; however, NLEIGS is an iterative procedure, while RESS is a direct sampling procedure and is thus more suitable for parallelization [29]. Our latest study indicated that the theory of RESS can be more general than that is known yet. RESS can pave a solid way for solving some typical problems in structural dynamics and the related fields. The present paper is devoted to solving the frequency sweep and nonlinear eigenvalue problems based on RESS. We consider the coupled BE-FE based vibro-acoustic analysis as the model problem because it has an extensive industrial application background and exhibits the common bottlenecks in the numerical treatment, as mentioned before; nevertheless, the theory and methods presented in this paper can be used analogously to structural dynamics and acoustics. There are two main contributions: • The existing theory of RESS are generalized so that RESS can also be used to construct projection subspaces for FSPs. Some practical issues concerning the numerical implementation of RESS and the proper choice of the controlling parameters are discussed and properly resolved. • A fast projection based frequency sweep algorithm is developed to solve FSPs in coupled BE-FE, and the RSRR algorithm is applied to solve the nonlinear EVPs. Numerical strategies for fast computation of the reduced system matrices are outlined. The performance of the algorithms are confirmed by academic and large-scale industrial problems. The remainding content is organized as follows. In Section 2, a coupled BE-FE algorithm for large-scale vibro-acoustic analysis is described. The BE part uses the fast directional BEM developed in [34]. In Section 3 the projection based fast frequency sweep method and the nonlinear EVP solver used in this paper are described. An strategy for the fast computation of the projected system matrices is put forward as well. In Section 4, the generalized theory of RESS is derived and some practical issues are discussed. In Section 5, firstly two academic examples are provided to study the efficiency and accuracy of the fast frequency sweep method and EVP solver in this paper, and then an industrial example is provided to show the performance of the methods in treating real-world problems. In Section 6 the paper is summarized with a discussion on the future work in this direction. 2. Fast BE-FE coupling for vibro-acoustics 2.1. Coupling method Throughout this paper, all field variables are assumed to be harmonic in time with the behavior eiωt . The imaginary unit is denoted by i and ω = 2πf is the angular frequency with the excitation frequency f . The BEM and FEM are used to model the acoustic domain Ωa and the structural domain Ωs , respectively. The governing boundary integral and differential equations of the two subdomains and their discretizations can be found in many text books, e.g., [35, 36], and are thus not repeated here. The coupled system is illustrated in Figure 1. The unknown acoustic pressure on the interface in the BEM model acts as a loading on the structure in the FEM model and the unknown structure vibration in the FEM model acts as the excitation for the BEM model. Thus we obtain a fully coupled formulation.

3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 1: A coupled system with FEM and BEM sub-domains The finite element equation of the structural domain Ωs is obtained by introducing an FE discretization in the weak variational formulation of the equilibrium equations of Ωs (based on the principle of virtual work) and subsequently applying a Galerkin procedure. This provides the following set of equations: ( ) Ks + iωCs − ω 2 Ms u = fs + fsa (1)

where u is the nodal displacement vector of the structure; fs denotes the force vector resulting from the body force on Ωs and the tractions applied on the boundary of Ωs , while fsa is due to the acoustic pressure on the interface of the structural and acoustic domains; Ks , Cs and Ms are the stiffness, damping and mass matrices, respectively. The bracketed term on the left-hand side of (1) is identified as the dynamic stiffness matrix K(ω) of the structure, i.e., K(ω) = Ks + iωCs − ω 2 Ms . The acoustic field Ωa is modeled by the Nystr¨om BEM in [34], leading to the following system of equations: H(ω)p(ω) = G(ω)q(ω),

(2)

where H and G are complex square matrices obtained by the Nystr¨om discretization of the single-layer and double-layer integral operators in the boundary integral equation associated with the Helmholtz partial differential equation; p and q are vectors of the nodal pressure p and its normal derivative ∂p/∂n. In the implementation of the Nystr¨om method, the boundary of the acoustic domain Ωa is first partitioned into curved quadratic elements, and then 6-point Gaussian quadrature rules on each element are used to compute the boundary integrals. This implementation makes the present Nystr¨om method very similar to the collocation method with discontinuous elements. However, as Nystr¨om method is a point-based discretization, it is easier to be accelerated by fast techniques [34]. Further, because of this similarity, the coupling algorithms presented below apply equally to the collocation or Galerkin BEM. The vector q is related to the boundary normal velocity vector, denoted by va , via the continuity condition q = −iωρva . Therefore, equation (2) becomes H(ω)p(ω) = −iωρG(ω)va (ω),

(3)

Besides, there hold another two relations at the interface: the equilibrium of the acoustic pressure and normal traction fsa = Csa p (4) and the continuity of the velocities of the fluid and the structure in the normal direction across the interface va = iωCas u,

(5)

where, Csa and Cas are coupling matrices accounting for the data transfer between the two domains, which will be described in Subsection 2.2. 4

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Combining coupling conditions (4) and (5) with equations (1) and (3) leads to the following fully coupled system of equations [ ]{ } { } K(ω) −Csa u fs = , (6) −ρω 2 G(ω)Cas H(ω) p 0 | {z } | {z } | {z } T (ω)

y(ω)

f (ω)

or, simply expressed as

T (ω)y(ω) = f (ω).

(7)

2.2. Data transfer on non-matching meshes and numerical implementation To achieve a good modelling flexibility, the fluid and structural domains are discretized independently. This will result in non-matching meshes across the interface, thus a data transfer strategy is required. Most researches in this aspect consider linear shape functions at interfacial meshes, which typically cause gaps and overlaps in modelling curved interfaces and thus may deteriorate the accuracy of the coupling method [37]. To improve the accuracy of the data transfer, in this paper quadratic elements are used for both the structure and acoustic domains. Specifically, the structure is discretized using eight-node isoparametric quadrilateral shell finite elements (SHELL 281 in ANSYS), and the boundary of the acoustic domain is discretized using six-node triangular boundary elements with the variables represented by quadratic shape functions defined at the six Gauss points on triangular boundary elements [34]. Another advantage of using quadratic elements is that they can well resist the pollution effect [38]. Also note that [37] is the first one to combine quadratic elements for both structural finite and acoustical boundary elements. According to (4) and (5), the coupling matrix Csa is used to convert the acoustic pressure p at the Gauss points of acoustic elements to the equivalent nodal force fsa at the nodes of structure elements, while Cas is used to convert the displacement u at the nodes of the structure elements to the velocity va at the Gauss points of the acoustic elements. Here we notice that usually the Lagrange interpolation functions of the boundary elements are constructed at the vertices of the elements and the midpoints of the edges. Therefore, for boundary elements data transfer operators from Gauss points to interpolation nodes, denoted by Png , and conversely from interpolation nodes to Gauss points, denoted by Pgn , are introduced, leading to the following expressions Csa = Tsa · Png and Cas = Pgn · Tas , (8) where Tsa is the usual data transfer matrix from the interpolation nodes of boundary elements to the nodes of finite elements, and Tas is the data transfer matrix from the nodes of finite elements to the interpolation nodes of boundary elements. In this paper the Mortar method is employed to couple the two fields. The method is built upon the weak formulation of the conservation of loads or displacements over the interface, and the resulted transfer matrices are given by [39, 40] Tsa = Ds−1 Dsa and Tas = Da−1 Das , (9) where, Ds =

∫

Γs

NsT Ns dΓ,

Dsa =

∫

Γs

NsT Na dΓ,

Da =

∫

Γa

NaT Na dΓ,

Das =

∫

Γa

NaT Ns dΓ,

(10)

Ns and Na are the interpolation function matrices of the FEM and BEM, and Γa and Γs are the boundary element and finite element meshes over the interface, respectively. Ds and Da are the same as the mass matrix of the FEM; they are symmetric and positive, and thus their inverses can be easily computed by Cholesky factorization. The computations of Dsa and Das are somewhat troublesome, since they involve the integration of the product of two basis functions associated with different meshes. One common approach is to detect the potentially overlapping element pairs, and perform polygon clipping to find the overlapping region, and then perform a decomposition of the clip polygon to define easy-to-integrate subdomains, but this process is very complicated. In this paper, we use the relatively simple element-based integration in [40]. This method only needs to get the parametric coordinates of the projections of the Gauss points of 5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

acoustic elements on the structure meshes. The general locations of these projection points can be efficiently determined by aid of the oct-tree structure in the fast BEM technique, and their parametric coordinates are computed by using the Newton method; see [40] for the details. 2.3. Solution scheme for coupled linear systems The coupled system (6) can be recast into the following reduced system [ ] H(ω) − ρω 2 G(ω)Cas K(ω)−1 Csa p = f (ω) where

(11)

f (ω) = ρω 2 G(ω)Cas K(ω)−1 fs .

The coefficient matrix of (11) is the Schur complement of the FE matrix K(ω) in the coupled system matrix (6). Once the pressure p is solved from (11), the displacement u can be computed as u = K(ω)−1 (fs − Csa p) .

(12)

Alternatively, one can obtain a reduced system by using the Schur complement of the BE matrix H(ω), but it is noticed that, for large-scale coupled FE-BE analysis, solving the reduced system (11) by iterative solvers is often more cost effective in terms of computational time and memory consumptions because the FE systems involving K(ω)−1 can be efficiently solved once the Cholesky decomposition of K(ω) is computed and saved. We use the generalized minimal residual method (GMRES) to solve the linear system (11), and use the fast acoustic BEM in [34] to accelerate the matrix-vector product with the densely populated BE matrices H and G. This fast BEM is based on a kernel-independent wideband fast directional algorithm for fast summation of oscillatory kernels [41], and it is capable of solving large-scale acoustic problems in a wide frequency band with almost linear complexity. It had been used to solve acoustic problems with nondimensional wave number up to 1000 and degrees of freedom up to several millions on a computer with 24 GB memory [34]. This enables the present fast BE-FE coupling method to be used for problems with relatively higher frequencies. The ILU (incomplete LU factorization) preconditioner based on the near field matrix of H(ω) is used to speed up the convergence of GMRES in solving (6). 3. Efficient solution of coupled frequency sweep and nonlinear eigenvalue problems The nonlinear dependence of the coupled matrix on the frequency poses a restriction on the choice of solution methods for the FSP and EVP problems. Among the available methods, the projection method provides a common framework for solving these two problems. In the first two subsections, the projection procedures for these two problems will be described respectively, and then in Subsection 3.3 interpolation techniques for the fast computation of the projected system matrices will be introduced to effectively reduce the computational cost for solving the projected systems of equations. 3.1. Frequency sweep problem Frequency sweeping refers to the computation of the frequency response function in a frequency band. For the coupled system (7), the response function is of the form y(ω) = [T (ω)]

−1

f (ω),

ω ∈ [ωL , ωR ],

(13)

where [ωL , ωR ] is the angular frequency band of interest. To keep the notations consistent in the whole paper, the frequency band will also be denoted by D, i.e., D = [ωL , ωR ]. The brute force computation of y(ω) consisting of solving a large number of linear systems of equations becomes cost prohibitive as the angular frequency band [ωL , ωR ] is enlarged, or its sampling is refined. Several techniques have been developed to alleviate the computational burden; representative examples include the model order reduction (MOR) methods based on interpolation and projection [9]. In general, 6

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

the interpolatory techniques entail the computation of the response y and its high-order derivatives at a serial of points, which could be troublesome for the FE-BE coupled systems due to the involvement of the complicated BE matrices. This paper utilizes a projection based approach. The key is to construct a proper subspace spanned by the columns of a rectangular orthogonal matrix Q, which encloses the response y at any frequency; that is, y can be expressed as (14) y(ω) = Qyr (ω), ω ∈ [ωL , ωR ], where yr (ω) is of a reduced dimension, denoted by nr , compared with y and solved from the projected system [ H ] Q T (ω)Q yr (ω) = QH f (ω), (15) | {z } | {z } fr (ω)

Tr (ω)

or simply written as

Tr (ω)yr (ω) = fr (ω).

(16)

Throughout this paper, the superscript H denotes the complex conjugate transpose. The system (16) can be solved with a much lower cost since its dimension nr is much smaller than n, the dimension of the original system (13). In the rest of this paper, we will refer to the above projection based method as the fast frequency sweep method (FSM). The quality of the search space Q directly controls the accuracy and computational cost of the approximation. It is also a key characteristic for a model order reduction method. In Section 4, the resolvent sampling method (RESS) will be proposed for constructing high-quality search spaces. The constructed subspace retains the information of the entire frequency band, and thus can lead to a highly accurate approximation. 3.2. Coupled eigenvalue problem The coupled EVP is obtained by setting the right-hand side of (6) to zero, ]{ } { } [ K(λ) −Csa u 0 = , −ρλ2 G(λ)Cas H(λ) p 0

(17)

where the frequency ω is replaced with the eigenvalue (or, eigen-frequency) λ. The nonlinear EVP (17) has the general form T (λ)v = 0, (18) where the system T (ω) is a n × n matrix-valued function depending nonlinearly on the frequency ω, and v denotes the eigenvector. In this paper, the RSRR algorithm in [29] will be used to solve the nonlinear EVP (17). Note that instead of solving (17) directly, one can formulate smaller EVPs by using the Schur complements of the coupled system matrix. For example, in [12] the following EVP is obtained from (17) using the Schur complement of the BE matrix H: [ ] 0 = K(λ) − ρλ2 G(λ)Cas H(λ)−1 Csa u (19) { [ ]} = Ks + iλCs − λ2 Ms + ρG(λ)Cas H(λ)−1 Csa u.

This formulation indicates that the effect of the fluid is equivalent to adding an additional frequencydependent mass to the structure, but it might not be suitable for solving large-scale problems, because it requires to solve a linear system with H in each iteration, which is computationally too expensive. To reduce the computational cost, Fritze et al. [42] reformulated the original EVP using the Schur complement of the FE matrix K(λ), [ ] (20) H(λ) − ρλ2 G(λ)Cas K(λ)−1 Csa p = 0. The reason is that the FE matrix K(λ) is sparse, and its Cholesky factorization can be efficiently computed and saved, leading to much faster iterations. However, comparing with the original formulation (17), 7

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

there is a problem with the Schur complement formulations (19) and (20): the coefficient matrices become non-analytic due to the inclusion of H(λ)−1 or K(λ)−1 . This poses a difficulty to the approximation of the coefficient matrices using finite term expansions, which is crucial for the efficient reduction of the computational cost for large-scale problems; see Subsection 3.3. Therefore, in this paper we will focus on the direct solution of the original EVP (17). The RSRR algorithm uses the classic Rayleigh-Ritz projection to extract the eigenpairs. Given a search subspace containing all the interested eigenvectors with an orthogonal basis Q, the eigenpairs of the EVP of the general form (18) can be obtained by solving the following reduced EVP Tr (λ)g = 0,

(21)

where Tr (z) = QH T (z)Q, the same as (16). EVPs (18) and (21) have the same eigenvalues, and their eigenvectors are related via v = Qg. The previous Rayleigh-Ritz projection lies in the heart of several eigen-solvers, e.g., the Arnoldi, rational Krylov or Jacobi–Davidson method [43]. The main difference is how to construct the search spaces, denoted by Q. Most existing eigen-solvers, including the aforementioned ones, generate search spaces iteratively. In contrast, the recently proposed resolvent moment scheme and resolvent sampling scheme construct search spaces in a direct manner [44, 29]. A salient feature of these direct schemes is its potential parallelism: the computationally most expensive part is completely decoupled and can hence be assigned to different processing units on a parallel computer. Regarding the moment and sampling schemes, the latter is more robust and accurate due to the effective avoidance of the possible ill-conditionness in the basis vectors [29], and is thus further studied and used in the present paper; see Section 4. 3.3. Fast computation of the reduced matrix Tr (z) The solution of the reduced FSP (16) and EVP (21) requires the repeated computation of the reduced matrix Tr (z), typically several hundreds or even thousands of times. A naive computation from Tr (z) = QH T (z)Q entails the computation of the original matrix T (z), resulting in a heavy computational burden for large-scale problems. In the following, a matrix interpolation technique is introduced to circumvent this problem. To facilitate the use of the interpolation technique, a simple scheme for estimating the number of interpolating points is described as well. 3.3.1. Matrix interpolation. To accelerate the computation of Tr (z), the matrix T (z) is first approximated by an interpolation of the following form d ∑ T (z) ≈ Tj ϕj (z), (22) j=1

where zi denote the interpolation points, ϕj (z) are interpolation basis functions, Tj ≡ T (zj ) are the coefficient matrices, and d denotes the number of interpolation points. Since T (z) is analytic on D, the value of d is often not very large. As a result, the reduced matrix Tr (z) can be approximated as Tr (z) ≈

d ∑

QH Tj Q ϕj (z). | {z } j=1

(23)

Tj′

The small nr × nr coefficient matrices Tj′ can be explicitly computed by matrix-vector multiplications with the coupled matrix T (zj ) and saved in memory. Since the computation of each Tj′ is mutually independent, this process is suitable for parallelization. Two cases have to be considered in constructing the interpolation (22). For FSPs and real EVPs, the approximation domain is typically a real interval, and expression (22) can be obtained by Chebyshev interpolation, leading to an exponential rate of convergence. This idea has been used in [29] to solve the 8

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

EVPs in the BEM. In a more general class of coupled problems involving energy dissipation, the eigenvalues are complex, and the domains of interest are two dimensional, e.g., an ellipse or a rectangle in the complex plane. In this case, the Chebyshev interpolation technique becomes inadequate since its approximation quality quickly deteriorates apart from the major axis of the associated Bernstein ellipse. As a way out, we extend the Cauchy interpolation technique for scalar functions in [45] to interpolating matrix-valued functions. Let zj and wj , j = 1, · · · , d, be a set of numerical quadrature points and the associated weights on the closed boundary contour ∂D. Then the Cauchy interpolation formula for T (z) is given by / d d ∑ T (zj )wj ∑ wj T (z) ≈ , (24) zj − z z −z j=1 j=1 j which is of the form (22) with ϕj (z) being given by wj ϕj (z) = zj − z

/∑ d j=1

wj . zj − z

(25)

For analytic functions, the Cauchy interpolation has an exponential convergence rate [45]. 3.3.2. Evaluation of d. Using the minimal number of terms in interpolation (22) is also important for the effective control of the total computational cost. Here a strategy to estimate the smallest number d is proposed. The idea stems from the fact that the accuracy of the matrix approximation (22) is determined by the quality of the same type of approximation to each entry t(z) of the matrix, t(z) ≈ t˜(z; d) :=

d ∑

t(zk )ϕj (z).

(26)

k=1

Thus, the smallest d required for [a given accuracy is controlled by ]the most “irregular” entry. The coupled matrix T (z) = K(z), −Csa ; −ρz 2 G(z)Cas , H(z) consists of two parts. One is the FE matrix K(z), whose entries are quadratic polynomials of z which are easy to approximate by polynomials. As a result, d should mainly be determined by the other part, i.e., the BE matrices G(z) and H(z), due to the fact that the entries of these two matrices are complicated functions of z whose expressions are not easy to obtain. A close observation shows that the entries of the BE matrices are generally of the form t(z) = f (z)eizr ,

0 6 r 6 R,

(27)

where f (z) is analytic in D, r is the distance between the source and target points on the boundary of the fluid domain, and R is the maximal length of the fluid domain. Above observation suggests a way to estimate d using the function (27). Given a target accuracy ε of the matrix approximation, d is chosen to be the minimal integer d′ such that the error of the approximation (26) in the domain D is less than ε, i.e., ( ) |t(z) − t˜(z; d′ )| d = Argmin 6 ε, z ∈ D . (28) |t(z)| d′ =1,2,... In practice, the above condition can be attained by a trial and error process. In most situations, d < 50 is satisfactory, thus the computational cost of this process would be negligible. The quality of the estimation (28) is controlled by the function f (z) and the parameter r. Generally, f (z) can be chosen from the monomials f (z) = z m , m = 0, 1, ..., with the exponent m being determined with a reference to the underlying BE matrices. For example, to approximate the G and H matrices in the conventional collocation BEM, m = 0 and m = 1 respectively can often result in a good estimation. In light of this, if the collocation BEM is used in the coupled FE-BE analysis, f (z) = z 2 can be used. The parameter r controls the oscillating frequency of eizr , and the maximal value of r, i.e., the largest size R of the structure, should contribute the most to the final value of d. Therefore, r = R can be used in most cases. 9

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

4. A general resolvent sampling method The resolvent sampling scheme for constructing eigenspaces has been recently proposed and enriched by the present authors’ group [29, 30, 31]. In this section, the scheme is extended to a more general framework where it can also be used to construct search spaces for FSPs. This new version of resolvent sampling scheme will be abbreviated as RESS. 4.1. Theoretical derivation Inspired by the original version of RESS in [29], one can imagine that for a general n × L excitation matrix P (z), the response Y (z) = T (z)−1 P (z) for all z in a domain D should lie in a small subspace. Here, L is assumed to be a small integer, L ≪ n. For theoretical exploration, we also assume that both T (z) and P (z) are analytic functions and the determinant det T (z) does not vanish identically inside the domain D. In most engineering FE and BE analyses, these assumptions can be easily satisfied. The verification of the above idea is based on the Keldysh theorem for holomorphic matrix-valued functions [22]. In a nutshell, this theorem states that the resolvent matrix T (z)−1 can be expressed as T (z)−1 = VD Φ(z)WDH + RD (z).

(29)

The meanings of the notations in (29) are as follows: • Φ(z) is a diagonal block matrix with Jordan block structure    1  Φ1 (z) Φk (z)     .. .. Φ(z) =   , Φk (z) =  , . . ηk ΦnD (z) Φk (z)   −1 −2 −µl,k (z − λk ) (z − λk ) ··· (z − λk )   .. .. ..   . . . Φlk (z) =  . −1 −2   (z − λk ) (z − λk ) (z − λk )−1

(30)

• λk (k = 1, · · · , nD ) are the mutually different eigenvalues in D, and nD is their number. ηk and µl,k are the dimension the nullspace of T (λk ) and the lth partial multiplicity of T (z) at λk , respectively. ∑nD ∑of ηk Let n ¯ D = k=1 l=1 µl,k be the total number of eigenvalues counting the algebraic multiplicity. • VD and WD are n × n ¯ D matrices consisting of the canonical systems of generalized eigenvectors of T (z) and T (z)H , respectively; that is, ( ) VD = vjl,k , 0 ≤ j ≤ µl,k − 1, 1 ≤ l ≤ ηk , k = 1, · · · , nD (31) and

( WD = wjl,k ,

) 0 ≤ j ≤ µl,k − 1, 1 ≤ l ≤ ηk , k = 1, · · · , nD .

(32)

Note that here the vectors in VD and WD are scaled such that µt,k ( j ∑ ∑

l,k wj−α

α=0 β=1

where Tj,k =

)H

Tα+β,k vµt,k = δtl δ0j , t,k −β

0 ≤ j ≤ µl,k − 1, 1 ≤ l, t ≤ ηk ,

1 (j) (λk ). j! T

• RD is a n × n holomorphic matrix-valued function in D, accounting for the contribution of the eigenvalues outside the domain D.

For more detailed description of the theorem and notations, the readers are referred to references [22, 23]. 10

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

In light of expression (29), the response can be written as Y (ω) = VD Φ(z)WDH P (z) +RD (z)P (z) {z } | P ′ (z)

(33)

′

= VD P (z) + RD (z)P (z).

Note that the n ¯ D × L matrix P ′ (z) is singular at the eigenvalues inside D, and the response associated with it is restricted to a small subspace spanned by the columns of VD . The term RD (z)P (z) represents the response of the frequencies outside D. Since both RD (z) and P (z) are analytic functions of z, this term permits the following form of expansion, RD (z)P (z) =

J ∑

Rj′ gj (z),

(34)

j=1

where gj (z) are basis functions and Rj′ are n × L coefficient vectors, respectively. There always exists a suitable set of basis functions gj (z) such that the expansion converges rapidly. For example, exponential convergence can often be achieved when analytic functions are approximated by Chebyshev polynomials in an interval. By plugging expansion (34) into (33), one obtains, Y (z) = VD P ′ (z) + where

J ∑ j=1

Rj′ gj (z) = Sˆ · Z(z),

] · · · RJ′ ,   ′ P¯ (z)  D1 (z)      Z(z) =  D2 (z)  ,  ..   .  DJ (z)

[ Sˆ = VD

R1′

R2′

(35)

(36)

(37)

Dj (z) = gj (z)IL and IL is the identity matrix of dimension L. Relation (35) has two aspects of implications. For EVPs, all the eigenvectors VD are enclosed in the ˆ denoted by span(S). ˆ Therefore, instead subspace spanned by the (¯ nD + J · L) columns of the matrix S, ˆ of using span(VD ), one can directly use the slightly larger subspace span(S) as an effective eigenspace in the Rayleigh-Ritz projection outlined in Section 3.2. Note that in order to correctly extract all the eigenvalues with right multiplicities, all the L columns of P are required to be linearly independent and the number L has to at least equal to the maximal algebraic multiplicity of eigenvalues of T in D, i.e., ∑ηbe k L ≥ maxk=1,··· ,nD ( l=1 µl,k ); see e.g., [21, 22]. In practical implementations, P can be chosen as a constant ˆ and thus span(S) ˆ can be used in random matrix. For FSPs, the response Y (z) lies in the subspace span(S), the projection-based MOR outlined in Section 3.1 to reduce the computational cost without loss of accuracy. ˆ becomes a common demand for solving both probTherefore, the construction of the subspace span(S) lems. 4.2. Numerical algorithms ˆ from the response T (z)−1 P (z) acThe RESS algorithm provides a way to extract a base of span(S) cording to (35). Given an excitation P (z), RESS consists of the following two basic steps: (1) Choose N sampling points zi (i = 1, . . . , N ) within D. 11

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(2) Compute the response Y (zi ) = T (zi )−1 P (zi ) (i = 1, . . . , N ) and form the n × (N · L) sampling matrix S = [Y (z1 ), Y (z2 ), · · · , Y (zN )] .

(38)

Then, it follows from (35) that ˆ S = Sˆ · [Z(z1 ), Z(z2 ), · · · , Z(zN )] = Sˆ · Z, | {z }

(39)

ˆ Z

ˆ is of the dimensions (¯ where the matrix Z nD + J · L) × (N · L). ˆ the rank of matrix Z ˆ has to be equal to n To retrieve span(S) ¯ D + J · L, which requires the values of N and L to fulfill the condition n ¯ D + J · L 6 N · L. (40) However, due to the fact that J is typically unknown, condition (40) can not be directly used in practice. Instead, we propose to use the following rule-of-thumb to determine the values of N and L: N · L = ⌈Cnl · n ¯D ⌉ ,

(41)

where ⌈·⌉ denotes the ceiling function, and Cnl is a user-defined number. In most cases, Cnl = 2.0 ∼ 3.0 can lead to a good balance between the overall accuracy and computational cost; this rule-of-thumb will be verified by numerical results in Section 5.2.1 and Section 5.3.2. Once a proper sampling matrix S is obtained, the projection matrix Q in (14) and (21) can then be computed by an orthogonalization of the columns of S. Before orthogonalization, it is important to normalize the column vectors of S so that they are of the same norm. This operation does not alter the resultant subspace but can greatly enhance the robustness of the method, because otherwise some vectors may have very large values when the sampling points zi are close to the eigenvalues of T (z) [31]. In this paper, the orthogonalization is accomplished by using the truncated SVD of S with tolerance δ, i.e., S ≈ QΣV ∗ . The reduced dimention nr is the numerical rank of S obtained by the truncated SVD. The tolerance δ controls the numerical rank of S. It cannot be too large, otherwise it will cause the loss of accuracy or even miss of eigenvalues. For computations with double precision, we suggest to use δ = 10−14 . 4.3. Practical issues As previously shown, the RESS algorithm can return a qualified eigenspace or response space if the parameters L, N and zi (i = 1, · · · , N ) are properly chosen. Now a practical question is how to determine these parameters. One way is to resort to some adaptive or iterative procedures, like the Jacobi-Davidson method [15] or the Krylov subspace method [16], however this will offset the parallelizability of the RSRR algorithm to a certain extent and thus will not be used in this paper. There are also problems concerning with the possible singularity of the coupled matrices and the scaling of the solution vectors (responses). In what follows, practical strategies will be put forward to properly deal with these issues so that the RESS algorithm will work with good robustness and efficiency. (1) Values for N and L. Here our discussion is confined to the use of the practical criterion (41). To this end, a rough estimation of the number of eigenvalues n ¯ D is required. There exist several ways to do this. Theoretically, n ¯ D of an analytic matrix T (z) can be estimated as [46]  ′  ∫ L ∑ 1  vlT T (z)−1 T ′ (z)vl  dz n ¯D ≈ 2πL′ i ∂D l=1 (42) ′ N L′ 1 ∑ ∑ T ≈ ωj vl T (zj )−1 T ′ (zj )vl , 2πL′ i j=1 l=1

12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

where, ∂D denotes the boundary contour of domain D, vl are sample vectors whose entries are 1 or −1 with equal probability, L′ is the number of sample vectors, T ′ (z) is the derivative of T (z) with respect to z, zj and ωj denote the numerical quadrature points and weights used to discretize the contour integral, and N ′ is the order of the quadrature rule. The above formula is fundamental and frequently used in literature. Besides, one can also estimate n ¯ D from the eigenvalues of a similar but easy-to-solve problem, such as a FE-BE discretization using a coarser mesh; or, for loosely coupled systems, one can directly estimate n ¯D from the natural frequencies of the structure. The estimation of n ¯ D for BEM using coarser meshes has been studied in our recent work [32]. For EVPs, the value of L has to be not less than the algebraic multiplicities of the eigenvalues. The latter can be estimated in several ways. A simple way is by judging from the symmetry property of the structure system. Given an estimation of the eigenvalues, one may also determine the algebraic multiplicities accordingly. In this term, there is a special situation where two computed eigenvalues are very close to one another; if so, it is suggested to consider them as an eigenvalue with multiplicity 2. Consequently, we suggest to set the value of L to be the estimated maximal algebraic multiplicity plus one. Of course, one can also use an even larger L, but this will increase the computational cost. Given the values of n ¯ C and L, N can be estimated from (41), e.g., N = ⌈Cnl · n ¯ C /L⌉. Usually, Cnl = 2 ∼ 3 can be used to achieve a good accuracy for both FSPs and EVPs. The theoretical and numerical studies in [31, 30] have shown that, compared with L, N has a more remarkable influence on the quality of the computed subspace Q. (2) Sampling points. Typically, the domain of a FSP or a real EVP is an interval. In this case, the sampling points of RESS can be chosen as a set of interpolation points in the interval, which corresponds to approximating RD (z)P (z) in (34) using an interpolation. For example, the Chebyshev points of the first kind are often a good choice; in the unit interval [−1, 1] they are given by zi = cos (2i−1)π (i = 1, . . . , N ). 2N For complex EVPs, the search domain D is usually a rectangle or an ellipse in the complex plane. Then, the sampling points are determined according to shape of the domain or the priori knowledge of the eigenvalues. In case of knowing an approximation of the eigenvalues, one can distribute the sampling points close to the eigenvalues. Otherwise, the sample points can either be chosen as a set of numerical quadrature points on the boundary of D or a special set of interpolation points inside the domain. An example of the latter case is the two-dimensional Chebyshev points in a rectangular domain. (3) Solvers for the coupled linear systems. When a sampling points zi is close to an (unknown) eigenvalue, the matrices T (zi ) will become nearly singular. This singularity may cause slow convergence or even failure of the GMRES solver. Specifically, the behavior of GMRES depends heavily on the orthogonalization method used to build up the Arnoldi basis. There are two commonly-used variants of GMRES, one conducts the orthogonalization using the classic Gram-Schmidt algorithm (CGS), denoted by GMRESCGS, and the other uses the modified Gram-Schmidt algorithm (MGS), denoted by GMRES-MGS; see e.g., [47]. The CGS is simpler and more suitable for parallelization, but may sometime lead to severe loss of orthogonality or even the loss of linear independence of the computed vectors; in comparison, the MGS has a better numerical stability. Our numerical experience indicates that for regular linear systems of equations these two variants of GMRES have almost the some convergence rate (number of iterations). However, when the sampling point |zi −λ| zi approaches an eigenvalue λ, e.g., max(|z < 10−2 , the convergence of GMRES-CGS tends to be slower, i |,|λ|) |zi −λ| while it has no considerable influence on the convergence of GMRES-MGS even when max(|z < 10−5 . i |,|λ|) Therefore, we suggest to use GMRES-MGS to solve the coupled linear systems in RESS. (4) Scaling of the sampling matrix. As the standard international units are used, the p and u components of the coupled system (6) are typically of totally different magnitudes. As a result, a suitable scaling of the p and u components of the solution vectors is required, otherwise either the structural part or the fluid part will dominate the basis vectors and the other domain will not be solved precisely. In this work we simply apply a scaling factor E, Young’s modulus of the structure material, to u. It turns out that this simple scaling can significantly improve the quality of the solutions for both FSPs and

13

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

EVPs. Therefore, the following scaled system is used in the projection operations ]{ } { } [ 1 −Csa Eu fs E2 K(ω) = . p 0 − ρω G(ω)C H(ω) as E

(43)

More specifically, the scaled system matrix in (43) will be used to compute the reduced matrices Tr in (16) and (21), and after solving the scaled FSP and EVP the scaling factor E is removed to obtain the true solutions. Finally, it is noticed that the vectors in S has to be scaled before the normalization operation mentioned in Section 4.2. 5. Numerical examples To show the correct implementation and computational performance of the present methods for solving coupled FSPs and EVPs, three numerical examples are considered herein. The first two examples are academic test cases, and the third one represents a common engineering application. All computations are carried out on a workstation with 12 cores (Intel Xeon E5-2643 v4 @ 3.40GHz) and 128 GB RAM. The linear systems of equations are solved using the GMRES solver with convergence tolerance 10−6 . The accuracy of the fast BEM and the accuracy controller ϵ of the matrix interpolation in (28) are set to be 10−6 as well. The accuracy of an eigenpair (λ, v) is measured by the relative residual ||T (λ)v||2 /||v||2 . For FSPs, the relative error of the response function computed by the fast FSM, denoted by y(ω), is evaluated by e(ω) =

||y(ω) − y ref (ω)||2 , ||y ref (ω)||2

ω ∈ [ωL , ωR ],

(44)

where y ref (ω) is the reference response obtained by directly solving the coupled linear system of the full size. 5.1. Spherical shell submerged in water The spherical shell depicted in Figure 2 is submerged in water. The geometrical and material parameters of this model are: radius 5.0 m, thichness 0.05 m, Young’s modulus 210 GPa, Poisson’s ratio 0.3, density of the structure 7900 kg/m3 , density of the water 1000 kg/m3 and speed of sound in water 1481 m/s. The structure is assumed to be undamped. The associated FSP and EVP will be solved over the range 0.5—200 Hz. To ensure enough accuracy, the structure is discretized into 2400 8-node quadrilateral elements SHELL 281 and the fluid boundary is divided into 4800 6-node triangular elements, leading to a coupled FE-BE model with 72012 DOFs. The mesh of the submerged spherical shell and the 2nd pressure mode resolved are shown in Figure 2.

Figure 2: Submerged spherical shell: Left: FEM model; Right: 2nd pressure mode extracted by RSRR

14

5.1.1. Frequency sweep The structure is assumed to be excited on its outer surface by a harmonic unit point force in the radial direction; see the left plot of Figure 2. The response over the aforementioned frequency band are computed by the brute force approach. In the brute force approach the frequency increment is 0.5 Hz, and totally 400 linear systems of the full size are solved. Analytical solutions of the coupled displacement and pressure can be found in, e.g., reference [2]. Figure 3 depicts the displacement at the location of the point load versus the analytical solution. Throughout the entire frequency band, the numerical results match accurately with the analytical solutions, even when the frequency is close to 0 where the displacement tends to infinite. In the frequency range between 70 and 200 Hz, the coupled resonances are dense, however, our numerical simulations predict these eigen-frequencies with high accuracy. We notice that this benchmark problem has been solved in Ref. [2] using a coupled fast BE-FE method. From the reported displacement responses in 0—100 Hz (Figure 3 in [2]) we see that the overall accuracy of our method is higher. The above results provide an evidence for the correctness of the numerical implementation and the good accuracy of our fast BE-FE coupled method.

2

displacement: mm

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

×10-5 analytical scheme numerical scheme

1.5 1 0.5 0

0

50

100

150

200

frequency: Hz Figure 3: Submerged spherical shell: displacement at the exciting point. 5.1.2. Eigenvalue problem Numerical modal analysis of the submerged spherical shell model is a challenging task [48], due to the fact that the natural frequencies are complex (with small imaginary part) and of high multiplicities. Here we solve this problem by using the RSRR algorithm. The frequency band below 200 Hz is considered. It can be identified from the response curves in Figure 3 that there are totally 18 mutually different eigenvalues. The blunt peak at the first resonance frequency indicates that the frequency is a complex number; for higher order natural frequencies, the peaks become sharper, indicating the smaller imaginary parts. This is in consistence with the theoretical prediction in [49]. However, from Figure 3 one can only evaluate the real parts of the natural frequencies, which are listed in the 3rd column of Table 1. To use the RSRR algorithm, one needs an estimation of the multiplicities of the eigenvalues. By the aid of the theoretical prediction in [49] and several trials, we found that the multiplicities of the eigenvalues increases regularly from 5 with an increment of 2; that is, the multiplicity of the first eigenvalue is 5, and that of the 18th becomes 39. So there are totally 396 eigenvalues below 200 Hz, which lie in the interval [50, 200] Hz; therefore, this interval is set as the search domain D in RSRR. The Chebyshev points with N = 30 are chosen to be the sampling points, and the same set of points are used to interpolate the coupled system matrix, i.e., we let d = 30. Besides, L = 40 is used. 15

Table 1: Submerged spherical shell: The 18 different resonance frequencies and their multiplicities Index 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

RSRR 55.813 + 1.181i 70.441 + 0.309i 80.560 + 0.044i 88.288 + 0.006i 94.656 + 0.004i 100.248 + 0.006i 105.441 + 0.009i 110.517 + 0.012i 115.709 + 0.015i 121.220 + 0.018i 127.228 + 0.022i 133.892 + 0.026i 141.349 + 0.031i 149.716 + 0.036i 159.086 + 0.042i 169.538 + 0.049i 181.129 + 0.057i 193.905 + 0.065i

Response curves 54.94 70.27 80.56 88.30 94.66 100.28 105.48 110.57 115.78 121.31 127.34 134.03 141.50 149.88 159.26 169.69 181.24 193.93

Multiplicity 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39

The RSRR algorithm correctly returns the 396 eigenpairs. The computed mutually different eigenvalues are listed in the 2nd column of Table 1; note that eigenvalues are recorded to three decimal places, and to this accuracy all the computed eigenvalues of the same order are regarded as the same one with high multiplicity. To measure the accuracy of RSRR, the relative residuals of the 396 eigenpairs are depicted by the curves in Figure 4. The effect of the scaling scheme discussed in Section 4.3 has also been tested. The results are illustrated in Figure 4. The curves with squares and circles show the residuals of the eigenpairs obtained with and without using the matrix scaling scheme, respectively. The simple matrix scaling in this paper reduces the residuals of the eigenpairs by more than two orders of magnitude.

100

relative error

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

non-scaling scaling

10-5

10-10

0

50

100

150

200

250

300

350

400

index of eigenpairs Figure 4: Residuals of the eigenpairs obtained with and without using the matrix scaling scheme In Ref. [48], the authors have solved the same eigenvalue problem by the contour integral method. They 16

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

reported the first two eigenvalues (Table 1 in [48]) to show the convergence of their method with respect to the refinement of the boundary mesh. The two converged eigenvalues therein are in good agreement with our results, i.e., the first two mutually different eigenvalues in Table 1. Also, we would like to emphasize that the contour integral method is inherently of lower accuracy and robustness than the Rayleigh-Ritz projection methods, including our RSRR algorithm and the SSRR algorithm in [23]; for theoretical and numerical comparisons, the readers are referred to Refs. [23, 29, 30]. 5.2. Rectangular panel backed by a closed cavity This example demonstrates the structural-acoustic coupling between an elastic rectangular panel backed by a closed rectangular cavity filled with water, as shown in Figure 5. The cavity has dimensions 1 m × 1 m × 1 m, and is backed by a steel panel on one side and by rigid walls on all other sides. The panel has dimensions 1 m × 1 m and thickness 0.01 m, and is simply supported on its edges. The physical parameters of the steel and water are the same as those in the submerged spherical shell model in Section 5.1. We consider the frequency sweep and eigenvalue problems in the frequency band 1—1000 Hz. The panel is partitioned into 900 8-node quadrilateral elements SHELL 281 and the boundary of the cavity is divided into 2580 6-node triangular elements, leading to a coupled FE-BE model with 37428 DOFs. The mesh of the model and the 2nd pressure mode resolved are shown in Figure 5.

Figure 5: Rectangular panel backed by a closed cavity: Left: FEM model; Right: the 2nd pressure mode extracted by RSRR

5.2.1. Frequency sweep Consider the scenario that a harmonic unit point force is applied at location (0.2,0.3) on the panel to excite the system in the frequency band 1—1000 Hz. We use this example to illustrate the performance of the proposed fast FSM, compared with the brute force approach, and at the same time, provide a guideline for choosing N (the number of sampling points) for the fast FSM. First, by comparing the response computed by the brute force approach and our fast BE-FE coupled method with the analytical solution developed in [50], we provide a further evidence for the correctness of the numerical implementation and the high accuracy of our fast BE-FE coupled method; see Figure 6a. In the brute force approach the frequency increment is 1 Hz, thus totally 1000 linear systems of the full size are solved. One can see that, throughout the entire frequency band, the numerical results show a good coincidence with the analytical solutions. The numerical simulations predict these eigen-frequencies with high accuracy. 17

Next, the performance of the fast FSM is tested. The sampling points are chosen to be the Chebyshev points of order N in [0, 1000] Hz. Given an estimate of the number of eigenvalues n ¯ D in the frequency band, N is obtained via relation (41), i.e., N = ⌈Cnl · n ¯ D ⌉. To show the influence of N on the overall accuracy of the fast FSM, we have run the computation under four cases, with Cnl = 1.0, 1.5, 2.0 and 4.0, respectively, and the error e(ω) defined in (44) is shown in Figure 6b; here n ¯ D = 37 is used, see the discussion about the EVP in Section 5.2.2. In all the computations, the coupled matrix T (z) is interpolated by the Chebyshev polynomials with d = 20 over the frequency band. One can see that the overall accuracy of the fast FSM improves when Cnl increases from 1 to 2, and after that, a larger Cnl does not improves the result remarkably but increases the computational cost. Therefore, Cnl around 2.0 is a good choice for this problem, which agrees with the experience rule as relation (41). When Cnl = 2.0, one only needs to solve N = 74 coupled systems, the total computational time is less than 1/13 of the brute force approach.

displacement: mm

100

analytical scheme numerical scheme

10-5

0

200

400

600

800

1000

frequency: Hz (a) Computed displacements versus the analytical results.

105

relative error

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

C n1=1.0 C n1=1.5 C n1=2.0

100

C n1=3.0

10-5

10-10

0

200

400

600

800

1000

frequency: Hz (b) Relative error of the fast frequency sweep method with different Cnl .

Figure 6: Rectangular panel backed by a closed cavity: performance of fast BE-FE coupled method and the fast frequency sweep method.

5.2.2. Eigenvalue problem The coupled eigenvalue problem is solved by the RSRR algorithm to obtain the natural frequencies and modes of the system within [0, 1000] Hz. A rough estimation using formula (42) indicates that there are 18

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

about 35 eigenvalues in this interval. The maximal multiplicity of the eigenvalues should be 2 due to the symmetry of the structure. In light of the criterions for determining N and L in Section 4.3, N = 60 and L = 2 are used. Since all the eigenvalues are real numbers, the Chebyshev points are chosen to be the sampling points. Chebyshev interpolation with d = 20 is used to approximate the coupled matrix. Totally 37 eigenpairs are extracted by the RSRR algorithm. The computed eigenvalues and the relative residuals of the eigenpairs are shown in Table 2. Several pairs of eigenvalues are very close to each other, i.e., the first two, the 6th and 7th; each of these pairs is likely to be an eigenvalue of multiplicity 2. Judging from the difference between the two-fold eigenvalues, the errors of the computed eigenvalues should be less than 1%. The much lower residuals of eigenpairs imply high accuracy of the eigenvectors. Table 2: Rectangular panel backed by a closed cavity: The 37 computed eigenvalues and the residuals k 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

λk 61.19 61.31 111.32 134.33 150.47 200.63 201.69 275.36 276.39 293.67 329.31 334.53 339.85 434.26 434.29 453.54 454.23 511.22 511.90

||T (λ)v||2 /||v||2 1.88 × 10−7 8.25 × 10−8 5.97 × 10−7 8.20 × 10−7 5.80 × 10−7 2.14 × 10−6 1.74 × 10−6 2.81 × 10−6 2.97 × 10−6 5.36 × 10−6 9.65 × 10−6 2.16 × 10−6 4.18 × 10−6 1.84 × 10−6 1.84 × 10−6 4.29 × 10−7 5.69 × 10−7 8.92 × 10−7 6.77 × 10−7

k 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37

λk 575.61 610.08 620.93 668.49 670.61 729.85 740.78 758.54 758.55 814.28 814.29 844.09 844.10 917.47 943.30 946.79 982.39 995.07 −−

||T (λ)v||2 /||v||2 1.93 × 10−7 3.60 × 10−6 4.80 × 10−6 3.21 × 10−6 3.50 × 10−6 2.75 × 10−7 1.51 × 10−7 1.51 × 10−6 1.51 × 10−6 1.51 × 10−7 1.51 × 10−7 1.51 × 10−6 1.51 × 10−6 1.51 × 10−7 1.51 × 10−7 1.51 × 10−7 1.51 × 10−7 1.51 × 10−7 −−

5.3. An industrial application Complex shell submerged in water is a prototypical structure in many navigation problems, such as the submarine, torpedo, etc. In the related applications, the vibroacoustic property of the entire system is one of the major concerns. Here we consider such a shell model, as shown in Figure 7. The geometrical and material parameters are: length 6.0 m, max radius 0.5 m, wall thickness 0.02 m, reinforcing ribs and partition thickness 0.05 m, Young’s modulus 210 GPa, Poisson’s ratio 0.3, density of the structure 7900 kg/m3 , density of the water 1000 kg/m3 and speed of sound in water 1481 m/s. The structure is assumed to be undamped. We will solve the frequency sweep and eigenvalue problems over the frequency band 1—500 Hz. To this end, the structure is discretized into 20702 8-node quadrilateral elements SHELL 281, where the number of elements on the coupling interface is 6900, and the fluid boundary is divided into 13800 6-node triangular elements, leading to a coupled FE-BE model with 420654 DOFs. 5.3.1. Eigenvalue problem The coupled eigenvalue problem is solved by the RSRR algorithm. A rough estimation using formula (42) indicates that there are around 10 eigenvalues in this interval. The maximal multiplicity of the eigenvalues should be 2 due to the symmetry of the structure, thus N = 20 and L = 2 are used. For this exterior problem, the eigenvalues should be complex numbers with small imaginary parts, hence the sampling points 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 7: FEM model of the cylindrical shell considered in Section 5.3. are still chosen to be the Chebyshev points in the real interval, and the coupled matrix is approximated by Chebyshev interpolation with d = 20. The RSRR algorithm returns 11 eigenpairs. The computed eigenvalues and the relative residuals of the eigenpairs are listed in Table 3. The relative residuals are all below 10−4 , implying a high accuracy of the results. All the 20 × 2 = 40 linear systems are solved with a single thread, the total computational time is 8.7 hours and the memory usage is around 14.1 GB. Figure 8 exhibits the displacement modal shapes of the first three mutually different eigenvalues. Table 3: cylindrical shell: The 11 computed eigenvalues and the residuals k 1 2 3 4 5 6

λk 135.037 + 0.0199i 141.785 + 0.0266i 263.579 + 0.0006i 299.379 + 0.8923i 320.978 + 1.4759i 343.394 + 0.9353i

(a) 135.037 + 0.0199i Hz

||T (λ)v||2 /||v||2 2.8 × 10−6 1.7 × 10−6 1.2 × 10−6 8.5 × 10−7 5.5 × 10−6 7.7 × 10−7

k 7 8 9 10 11

λk 365.974 + 0.6099i 415.664 + 2.4130i 417.560 + 0.4813i 461.466 + 3.2674i 498.289 + 0.5245i −−

(b) 141.785 + 0.0266i Hz

||T (λ)v||2 /||v||2 3.4 × 10−6 3.7 × 10−5 4.1 × 10−6 1.2 × 10−6 2.4 × 10−6 −−

(c) 263.579 + 0.0006i Hz

Figure 8: Displacement eigenmodes corresponding to the first three eigenvalues of the cylindrical shell model.

5.3.2. Frequency sweep As shown in Figure 7, a harmonic pressure with 1 N/m2 amplitude is applied at the rear face and the leftmost partition of the cylindrical shell along the X and Y directions, respectively, to excite the system. The FSP is solved by using the brute force approach and the proposed fast FSM respectively, and the responses at points 1, 2 and 3 are displayed for comparison. In the brute force approach, the frequency increment is 1 Hz, thus totally 500 linear systems of the full size are solved. The total computational time 20

is 101.4 hours, and the memory usage is around 13.8 GB. Figure 9a depicts the displacements at points 1,2 and 3, along the three directions, computed by the brute force approach. The fast FSM is run under three cases with Cnl = 2.0, 2.5 and 3; note that from the results in Section 5.3.1 n ¯ D = 11 is used to evaluate N = Cnl · n ¯ D . Figure 9b shows the error distributions corresponding to different Cnl . Clearly, Cnl = 2.5 or 3.0 is a good choice for this problem, which is again in consistence with the rule-of-thumb regarding (41). When Cnl = 3.0, only 33 linear systems of the full size are solved, and the total computational time is 9.7 hours, which is 10 times lower than that of the brute force approach.

displacement: mm

100

X-Position 1 X-Position 2 X-Position 3 Y-Position 1 Y-Position 2 Y-Position 3 Z-Position 1 Z-Position 2 Z-Position 3

10-5

10-10

10-15

0

100

200

300

400

500

frequency: Hz (a) Computed displacements at position 1,2 and 3.

C n1=2.0

relative error

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

C n1=2.5

10-4

C n1=3.0

10-6

0

100

200

300

400

500

frequency: Hz (b) Accuracy of the fast frequency sweep method.

Figure 9: Cylindrical shell: performance of fast BE-FE coupled method and the fast frequency sweep method.

6. Conclusions In this paper, firstly, a fast coupled BE-FE method has been described for solving the sound-structure interaction problems. The method allows for non-matching meshes across the interfaces. An efficient data transfer strategy for curved quadratic elements is employed to raise the accuracy. The BE acoustic analysis is accelerated by a kernel-independent wideband fast directional algorithm. All these features make the present coupled method fast and accurate in a wide frequency range. Secondly, a generalized resolvent sampling method (RESS) has been developed for solving modal analysis and frequency sweep problems in the coupled BE-FE based sound-structure interaction problems. In modal 21

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

analysis, the coupled nonlinear eigenvalue problem is solved by the recently proposed RSRR method, in which RESS is used to construct eigenspaces for the Rayleigh-Ritz projection procedure. In frequency sweeping analysis, a projection-based fast sweep method has been established, in which RESS serves to construct response spaces for the condensation of the original coupled equation using the Galerkin projection. A numerical scheme for the fast computation of the reduced system matrices in both the RSRR method and the fast sweep method is proposed, and several rules for practical applications of these two methods are also discussed. Here we mention that based on the proposed RESS method one can also construct modal analysis and frequency sweep methods for structural dynamics and acoustics. These two achievements have led to an efficient numerical tool for the frequency response computation, frequency sweeping and modal analysis of engineering sound-structure interaction problems. The computational performances of the proposed methods have been tested and confirmed by academic benchmarks and a large-scale industrial application. Acknowledgements The third author gratefully acknowledges the financial supports from the National Science Foundations of China under Grants 11102154 and 11472217, Fundamental Research Funds for the Central Universities in China and the Alexander von Humboldt Foundation (AvH) in Germany. References [1] G. C. Everstine, F. M. Henderson, Coupled finite element/boundary element approach for fluid–structure interaction, Journal of the Acoustical Society of America 87 (5) (1990) 1938–1947. [2] D. Brunner, M. Junge, L. Gaul, A comparison of FE–BE coupling schemes for large-scale problems with fluid–structure interaction, International Journal for Numerical Methods in Engineering 77 (5) (2009) 664–688. [3] T. van Opstal, E. van Brummelen, R. de Borst, M. Lewis, A finite-element/boundary-element method for largedisplacement fluid-structure interaction, Computational Mechanics (2012) 1–10. [4] L. Heltai, J. Kiendl, A. DeSimone, A. Reali, A natural framework for isogeometric fluid–structure interaction based on BEM–shell coupling, Computer Methods in Applied Mechanics and Engineering 316 (2017) 522–546. [5] J. Maestre, J. Pallares, I. Cuesta, M. A. Scott, A 3D isogeometric BE–FE analysis with dynamic remeshing for the simulation of a deformable particle in shear flows, Computer Methods in Applied Mechanics and Engineering 326 (2017) 70–101. [6] Z. Liu, M. Majeed, F. Cirak, R. N Simpson, Isogeometric FEM-BEM coupled structural-acoustic analysis of shells using subdivision surfaces, International Journal for Numerical Methods in Engineering 113 (9) (2018) 1507–1530. [7] H. Peters, N. Kessissoglou, S. Marburg, Modal decomposition of exterior acoustic-structure interaction problems with model order reduction, The Journal of the Acoustical Society of America 135 (5) (2014) 2706–2717. [8] B. Besselink, U. Tabak, A. Lutowska, N. V. D. Wouw, H. Nijmeijer, D. J. Rixen, M. E. Hochstenbach, W. H. A. Schilders, A comparison of model reduction techniques from structural dynamics, numerical mathematics and systems and control, Journal of Sound & Vibration 332 (19) (2013) 4403–4422. [9] U. Hetmaniuk, R. Tezaur, C. Farhat, Review and assessment of interpolatory model order reduction methods for frequency response structural dynamics and acoustics problems, International Journal for Numerical Methods in Engineering 90 (13) (2012) 1636–1662. [10] H. Peters, N. Kessissoglou, S. Marburg, Modal decomposition of exterior acoustic-structure interaction, The Journal of the Acoustical Society of America 133 (5) (2013) 2668–2677. [11] F. Tisseur, K. Meerbergen, The quadratic eigenvalue problem, SIAM Review 43 (2) (2001) 235–286. [12] M. Junge, D. Brunner, L. Gaul, Solution of FE-BE coupled eigenvalue problems for the prediction of the vibro-acoustic behavior of ship-like structures, International Journal for Numerical Methods in Engineering 87 (7) (2011) 664–676. [13] U. Tabak, D. J. Rixen, Vibro-Lanczos, a symmetric Lanczos solver for vibro-acoustic simulations, International Journal for Numerical Methods in Engineering 107 (4) (2016) 290–311. [14] V. Mehrmann, C. Schr¨ oder, Nonlinear eigenvalue and frequency response problems in industrial practice, Journal of Mathematics in Industry 1 (1) (2011) 1–18. ´ [15] C. Effenberger, Robust solution methods for nonlinear eigenvalue problems, Ph.D. thesis, Ecole polytechnique f´ ed´ erale de Lausanne (2013). [16] R. Van Beeumen, Rational Krylov methods for nonlinear eigenvalue problems, Ph.D. thesis, KU Leuven (2015). [17] A. Ali, C. Rajakumar, S. Yunus, Advances in acoustic eigenvalue analysis using boundary element method, Computers & Structures 56 (5) (1995) 837–847. [18] G. R. Tai, R. P. Shaw, Helmholtz-equation eigenvalues and eigenmodes for arbitrary domains, The Journal of the Acoustical Society of America 56 (3) (1974) 796–804. [19] G. De Mey, Calculation of eigenvalues of the Helmholtz equation by an integral equation, International Journal for Numerical Methods in Engineering 10 (1) (1976) 59–66.

22

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[20] S. M. Kirkup, S. Amini, Solution of the Helmholtz eigenvalue problem via the boundary element method, International Journal for Numerical Methods in Engineering 36 (2) (1993) 321–330. [21] J. Asakura, T. Sakurai, H. Tadano, T. Ikegami, K. Kimura, A numerical method for nonlinear eigenvalue problems using contour integrals, JSIAM Letters 1 (0) (2009) 52–55. [22] W.-J. Beyn, An integral method for solving nonlinear eigenvalue problems, Linear Algebra and Its Applications 436 (10) (2012) 3839–3863. [23] S. Yokota, T. Sakurai, A projection method for nonlinear eigenvalue problems using contour integrals, JSIAM Letters 5 (0) (2013) 41–44. [24] H. Gao, T. Matsumoto, T. Takahashi, H. Isakari, Eigenvalue analysis for acoustic problem in 3D by boundary element method with the block Sakurai–Sugiura method, Engineering Analysis with Boundary Elements 37 (6) (2013) 914–923. [25] A. Leblanc, A. Lavie, Solving acoustic nonlinear eigenvalue problems with a contour integral method, Engineering Analysis with Boundary Elements 37 (1) (2013) 162–166. [26] C.-J. Zheng, H.-F. Gao, L. Du, H.-B. Chen, C. Zhang, An accurate and efficient acoustic eigensolver based on a fast multipole bem and a contour integral method, Journal of Computational Physics 305 (2016) 677–699. [27] C. Zheng, C. Zhang, C. Bi, H. Gao, L. Du, H. Chen, Coupled FE–BE method for eigenvalue analysis of elastic structures submerged in an infinite fluid domain, International Journal for Numerical Methods in Engineering 110 (2) (2017) 163–185. [28] C.-J. Zheng, C.-X. Bi, C. Zhang, H.-F. Gao, H.-B. Chen, Free vibration analysis of elastic structures submerged in an infinite or semi-infinite fluid domain by means of a coupled FE–BE solver, Journal of Computational Physics 359 (2018) 183–198. [29] J. Xiao, S. Meng, C. Zhang, C. Zheng, Resolvent sampling based Rayleigh–Ritz method for large-scale nonlinear eigenvalue problems, Computer Methods in Applied Mechanics and Engineering 310 (2016) 33–57. [30] J. Xiao, C. Zhang, T.-M. Huang, T. Sakurai, Solving large-scale nonlinear eigenvalue problems by rational interpolation and resolvent sampling based Rayleigh–Ritz method, International Journal for Numerical Methods in Engineering 110 (8) (2017) 776–800. [31] J. Xiao, H. Zhou, C. Zhang, C. Xu, Solving large-scale finite element nonlinear eigenvalue problems by resolvent sampling based Rayleigh-Ritz method, Computational Mechanics 59 (2) (2017) 317–334. [32] J. Xiao, J. Wang, T. Liang, L. Wen, The RSRR method for solving large-scale nonlinear eigenvalue problems in boundary element method, Engineering Analysis with Boundary Elements 93 (2018) 150–160. [33] S. G¨ uttel, R. Van Beeumen, K. Meerbergen, W. Michiels, NLEIGS: A class of fully rational Krylov methods for nonlinear eigenvalue problems, SIAM Journal on Scientific Computing 36 (6) (2014) A2842–A2864. [34] Y. Cao, L. Wen, J. Xiao, Y. Liu, A fast directional BEM for large-scale acoustic problems based on the Burton–Miller formulation, Engineering Analysis with Boundary Elements 50 (2015) 47–58. [35] M. Kaltenbacher, Computational Acoustics, Springer, 2018. [36] T. J. Hughes, The finite element method: linear static and dynamic finite element analysis, Courier Corporation, 2012. [37] H. Peters, S. Marburg, N. Kessissoglou, Structural-acoustic coupling on non-conforming meshes with quadratic shape functions, International Journal for Numerical Methods in Engineering 91 (1) (2012) 27–38. [38] S. Marburg, A pollution effect in the boundary element method for acoustic problems, Journal of Theoretical and Computational Acoustics 26 (02) (2018) 1850018. [39] A. de Boer, A. H. van Zuijlen, H. Bijl, Comparison of conservative and consistent approaches for the coupling of nonmatching meshes, Computer Methods in Applied Mechanics and Engineering 197 (49) (2008) 4284–4297. [40] P. Farah, A.-T. Vuong, W. Wall, A. Popp, Volumetric coupling approaches for multiphysics simulations on non-matching meshes, International Journal for Numerical Methods in Engineering 108 (12) (2016) 1550–1576. [41] B. Engquist, L. Ying, Fast directional multilevel algorithms for oscillatory kernels, SIAM Journal on Scientific Computing 29 (4) (2007) 1710–1737. [42] D. Fritze, S. Marburg, H.-J. Hardtke, FEM–BEM-coupling and structural–acoustic sensitivity analysis for shell geometries, Computers & structures 83 (2-3) (2005) 143–154. [43] V. Mehrmann, H. Voss, Nonlinear eigenvalue problems: A challenge for modern eigenvalue methods, GAMM-Mitteilungen 27 (2) (2004) 121–152. [44] T. Sakurai, H. Sugiura, A projection method for generalized eigenvalue problems using numerical integration, Journal of Computational and Applied Mathematics 159 (1) (2003) 119–128. [45] A. P. Austin, P. Kravanja, L. N. Trefethen, Numerical algorithms based on analytic function values at roots of unity, SIAM Journal on Numerical Analysis 52 (4) (2014) 1795–1821. [46] Y. Maeda, Y. Futamura, T. Sakurai, Stochastic estimation method of eigenvalue density for nonlinear eigenvalue problem on the complex plane, JSIAM letters 3 (2011) 61–64. [47] L. Giraud, J. Langou, M. Rozloznik, The loss of orthogonality in the Gram-Schmidt orthogonalization process, Computers & Mathematics with Applications 50 (7) (2005) 1069–1075. [48] A. Kimeswenger, O. Steinbach, G. Unger, Coupled finite and boundary element methods for fluid-solid interaction eigenvalue problems, SIAM Journal on Numerical Analysis 52 (5) (2014) 2400–2414. [49] Y. K. Lou, T. C. Su, Free oscillations of submerged spherical shells, Journal of the Acoustical Society of America 63 (5) (1978) 1402–1408. [50] A. J. Pretlove, Forced vibration of a rectangular panel backed by a closed rectangular cavity, Journal of Sound & Vibration 2 (3) (1965) 197–209.

23