A subspace iteration eigensolver based on Cauchy integrals for vibroacoustic problems in unbounded domains

Despite the potential and the increasing popularity of the boundary element method (BEM), modal analyses based on BEM are not yet put into engineering practice, mainly due to the lack of efficient solvers for the underlying nonlinear eigenvalue problem (EVP). In this article, we review a subspace iteration method based on FEAST for the solution of vibroacoustic EVPs involving the finite element method (FEM) and BEM. The subspace is obtained by applying a spectral projector and is computed by contour integration, whereas the contour is also used to subsequently solve the projected EVP by rational approximation. The computation of the projection matrices is addressed by two approaches. In the case of heavy fluid loading, we solve the underlying coupled linear systems by an iterative block Krylov method. In the case of light fluid loading, we exploit the fact that the coupled system admits accurate model order reduction solely based on the structural subsystem. Applications to a spherical shell and to a musical bell indicate that only a few contour points are required for an accurate solution without inducing spurious eigenvalues. The results are compared with those of a contour integral method and illustrate the efficiency of the proposed eigensolver.


INTRODUCTION
Vibroacoustics is a branch of physics dealing with oscillations of solid structures (vibrations) and mechanical waves in gases or liquids (acoustics). The study of vibroacoustics spans over various engineering disciplines often aiming at improving insulation and radiation characteristics of structures. For example, the mutual interaction between acoustics and elasticity is an essential aspect in the design process of phononic crystals 1 that are used to attenuate wave propagation in certain frequency bands. Our recent work has shown that acoustic loading can also have a significant effect on the dynamic properties of lightweight sandwich structures. 2 Moreover, it is crucial to consider these two physical phenomena simultaneously in numerical analysis, in order to enable meaningful comparisons to experiments, 3 which of course are conducted in a real multiphysics world.
Fully coupled numerical formulations describing the mutual interaction between structural vibrations and acoustic wave propagation are well established in the scientific literature 4,5 and appearing in commercial practice as well. The elastodynamic equations underlying the vibrating structure are typically discretized by the finite element method (FEM), 6 whereas several methods are available for addressing the surrounding, often unbounded acoustic domain. Here, we will consider the direct boundary element method (BEM) with a collocation discretization. 7 Compared with other strategies for modeling the acoustic field, BEM is particularly convenient in the case of unbounded domains. The discretization is restricted to the sound radiating surface and BEM does not require special boundary conditions for truncating the far-field sound radiation. Those advantages however come at the cost of an implicitly frequency-dependent kernel, which poses significant challenges when performing frequency sweep analyses or modal analyses.
In the case of frequency sweep analyses with BEM, the computational effort mainly stems from the assembly of coefficient matrices at each frequency point of interest along with the solution of the resulting sequence of linear systems. Early works on this topic focused on interpolating either the coefficient matrix, 8 the underlying kernel function, 9 or even the resulting response itself. 10 More recent work is concerned with efficient strategies for solving the sequence of linear systems by means of model order reduction [11][12][13] and greedy sampling. 14 Along with frequency sweeps, modal analysis accounts for another fundamental tool to study the dynamical behavior and provides system-inherent properties such as eigenfrequencies and modal damping values. In the case of bounded acoustic domains, modal analyses based on FEM 15 or BEM with frequency-independent coefficient matrices 16,17 are well established. However, in the context of coupled structural exterior acoustic systems involving BEM, modal analyses require the solution of the underlying eigenvalue problem (EVP) that depends nonlinearly on the eigenvalue parameter. The solution of such nonlinear EVPs arising from frequency-dependent boundary element (BE) matrices dates back to the 1970s 18 and has since then drawn the attention of researchers from both mathematical and engineering communities. Despite the presence of this topic for nearly half a century and the growing number of boundary element software as well as eigensolver packages, 19 modal analyses based on BEM are not yet put into engineering practice. This hesitation may be largely attributed to the complicated choice of the adequate solver for the specific application at hand. For example in our case, an important distinction concerns the extent of fluid loading. Here, most papers focus on applications with heavy fluid loading 11,20,21 and address them by forming the Schur complement of the finite element (FE) submatrix with respect to the coupled system. 11,21 But in the case of light fluid loading, it may be more efficient and intuitive to form the Schur complement of the BE submatrix and treat the acoustic field as additional mass and damping contributions. Furthermore, most researchers rely on direct eigensolvers in order to avoid repeated assembly of frequency-dependent BE matrices. In the case of light fluid loading, iterative eigensolvers may actually be beneficial, since the in vacuo modes of such structures serve as a good initial guess and facilitate finding the wet (shifted) eigenfrequencies and modal damping values. The in vacuo modes may even be used a posteriori to check the occurrence of spurious eigenfrequencies.
The different approaches for solving nonlinear EVPs can be categorized into algorithms based on Newton's method, those based on polynomial and rational approximation 22 and methods based on contour integration. 23 The use of Newton's method is typically accompanied by deflation and hence requires subsequent computation of individual eigenpairs. Polynomial approximations are typically based on frequency sampling at Chebychev nodes 24 or truncated Taylor series. 20,25 When dealing with complex eigenvalues, the quality of these polynomial approximations quickly deteriorate when the eigenvalue lies apart from the real axis. Moreover, linerization of higher order polynomials induces numerical instabilities and thus precludes accurate computation of eigenpairs in large frequency ranges. Contour integral methods 21,26,27 account for another family of eigensolvers and are particularly appealing because of the low memory requirements and since the main computations can be executed in parallel. They essentially work by projecting the nonlinear EVP into a generalized EVP. However, the projection often leads to ill-conditioning 28 and results in spurious eigenvalues. Moreover, accurate projection via contour integration often requires a large number of integration points at which the BE matrices have to be explicitly assembled. Finally more recently, resolvent sampling based model order reduction has been employed for a Rayleigh-Ritz procedure. 11,29 In this article, we extend two recent approaches-the nonlinear FEAST algorithm 30 and the linearization of rational approximations 31,32 -in order to solve structural acoustic EVPs with high accuracy and efficiency. The nonlinear FEAST (NLFEAST) is a generalization of the popular FEAST algorithm, 33 for which an open source library package 34 exists, and works by iteratively refining a projection matrix that is used for a Rayleigh-Ritz projection. However, the library package as well as the applications in the aforementioned paper 30 are limited to polynomial EVPs, leaving problems out of consideration that depend nonlinearly on the eigenvalue parameter. Structural acoustic EVPs fall into the latter category and-after projection by NLFEAST-still require solution of a nonlinear EVP. We address this issue by employing a rational approximation based on the Cauchy integral equation and subsequent linearization. 32 This enables us to simply reuse the information for computing the projection matrix. Hence, the actual solution of the EVP comes at a negligible computational cost. Moreover, the same rational approximation is also used to avoid explicit assembly of the original system matrix at intermediate eigenvalues over the course of NLFEAST iterations. Finally, in order to enable fast computation of the projection matrix, we present two tailored approaches for solving the underlying coupled linear systems: In the case of heavy fluid loading, we form the Schur complement of the structural subsystem and employ an iterative block Krylov solver. In the case of light fluid loading, we form the Schur complement of the acoustic subsystem and employ an unilateral, frequency-independent model order reduction, which is computed solely based on the structural subsystem. The presented numerical framework enables efficient and accurate modal analysis of structural acoustic interaction without inducing spurious eigenfrequencies. Applications to an academic example of a spherical shell in water as well as a musical bell in air illustrate the performance of the method.

COUPLED FINITE AND BOUNDARY ELEMENT METHODS FOR VIBRATION ANALYSIS IN UNBOUNDED ACOUSTIC DOMAINS
We consider time-harmonic problems involving a mutual interaction between structural vibrations and the surrounding acoustic field. The underlying equation of linear elasticity is addressed by FEM, resulting in the linear system Therein, u ∈ C n s is a vector containing the n s unknown displacement degrees of freedom (DOF). The angular frequency is defined as = 2 f , where f is the frequency. The stiffness and mass matrices of the structure are denoted with K, M ∈ R n s ×n s . Structure-inherent damping is not considered in this work, but a damping matrix can be added to Equation (1) without affecting the further derivations. The structure is excited by external and fluid forces f s , f f ∈ C n s on the nodes. The latter act by virtue of the acoustic field which can be described by the Helmholtz equation. After a collocation discretization with BEM, the linear system of equations for the acoustic subdomain reads where p ∈ C n f is the vector containing the n f unknown sound pressure DOFs. The BE matrices H( ), G( ) ∈ C n f ×n f relate the structural particle velocity v s ∈ C n f to the sound pressure. Further, the acoustic field may be excited by a sound source with an incident pressure field p i ∈ C n f . Quadrilateral boundary elements with discontinuous bilinear sound pressure interpolation are used in this work. Hence, neighboring boundary elements do not share their interpolation nodes.
Equations (1) and (2) are mutually coupled on the submerged surface of the structure. There, the structure is subject to normal tractions due to the acoustic sound pressure, and the particle velocity v s equals the time derivative of the normal displacement on the boundary. Due to the use of discontinuous boundary elements, the interpolation nodes of the acoustic and structural subdomains do not coincide. Hence, the coupling conditions are reformulated in the weak sense. After discretization by the Bubnov-Galerkin approach, the coupling conditions can be expressed as 35 where C sf ∈ R n s ×n f and C fs ∈ R n f ×n s are the mesh coupling matrices that relate the displacement and pressure DOFs, and i denotes the imaginary unit. Combining Equations (1) to (3) yields the fully coupled linear system By forming the Schur complement of H( ) with respect to L( ) and thereby eliminating the pressure DOFs from Equation (4), we obtain a structural equation that includes the effect of acoustic loading, 20 that is, in which i C sf H −1 ( )G( )C fs can be interpreted as the additional mass and damping contributions due to the acoustic field. The excitation on the right-hand side of Equation (5) comprises structural loading f s as well as acoustic loading by virtue of the incident pressure field p i . Alternatively, the displacement DOFs can be eliminated from Equation (4) by forming the Schur complement of K − 2 M with respect to L( ), which yields Equation (6) can be interpreted as an acoustic equation incorporating the elasticity of the structure as admittance boundary condition. 36 By setting the right-hand sides of Equations (5) and (6) to zero, we obtain the structural acoustic EVPs respectively. The nonzero vectors and are the displacement and pressure modes, respectively. The EVPs in Equations (7) and (8) are nonlinear, since the BE matrices H(̃) and G(̃) implicitly depend on the complex eigenvalue parameter̃. The real part of̃corresponds to the eigenfrequency and the imaginary part is associated with modal radiation damping. It is well known that if the matrix H(̃) is invertible, theñis an eigenvalue of L(̃) if and only if̃is an eigenvalue of Equation (7). Similarly, if K −̃2M is invertible, theñis an eigenvalue of L(̃) if and only if̃is an eigenvalue of Equation (8). Furthermore, the nonzero vectors with L(̃) = 0 can be related with the eigenvectors and , respectively. 37 Both EVP formulations (7) and (8) have been considered in the literature, 2,11,20,21 but the proper choice between them has not yet been discussed in a comprehensive manner. Most researchers focus on applications with heavy fluid loading and address them by Equation (8). The linear systems that arise over the course of solving Equation (8) can be addressed by an iterative method once a factorization of the sparse and symmetric matrix K −̃2M is computed. On the other hand, in the case of light fluid loading, it is more intuitive to form the Schur complement of the BE submatrix and treat the acoustic field as additional mass and damping contributions. Moreover, we will later see that the use of Equation (7) admits efficient unilateral model order reduction in the case of light fluid loading. We will consider both formulations (7) and (8) in the following and provide tailored solution schemes for both of them.

FEAST FOR NONLINEAR STRUCTURAL ACOUSTIC EIGENVALUE PROBLEMS
Introduced in the context of linear EVPs, the original FEAST algorithm 33,38 essentially is a subspace iteration method which uses a Rayleigh-Ritz procedure in each iteration, that is, where A is the matrix of interest, Q is the complex-valued projection matrix, and (⋅) H denotes the Hermitian transpose. The reduced matrix A r contains some approximate eigenvalues of A. Those eigenvalues are obtained by solving the EVP A r X r = X r , and the approximate eigenvectors are retrieved from X = QX r . Different from, for example, the Arnoldi algorithm which uses Krylov subspaces for Q, FEAST constructs a subspace of fixed dimension by applying a spectral projector (A) to the desired eigenvectors. The spectral projection is accomplished by a Cauchy integral representation of the resolvent of A, yielding the complex contour integration where I is the identity matrix. The shift z is defined along the contour -a user-defined Jordan curve in the complex plane. Hence, FEAST yields the eigenvalues and corresponding eigenvectors lying in the domain  that is enclosed by , see Figure 1. The numerical effort of FEAST is determined by the accuracy of the integration in Equation (10), that is, the number of discrete integration points. The computation of the integrand at the individual integration points can be executed in parallel. Gavin et al. 30 have recently provided another interpretation of the original FEAST algorithm by mentioning that it is essentially a generalized shift-and-invert iteration. Standard shift-and-invert iteration uses a single shift z, whereas the contour integration in Equation (10) corresponds to using multiple shifts simultaneously. Inspired by this, Gavin et al. 30 have proposed a nonlinear version of FEAST (NLFEAST) which in turn is a generalization of the residual inverse iteration. 39 Again, contour integration is used in NLFEAST to handle multiple shifts efficiently. However, the authors limit their applications to quadratic EVPs leaving open the application to analytic matrix functions that depend nonlinearly on the eigenvalue parameter. The next section briefly reviews the NLFEAST algorithm, whereafter Sections 3.2 to 3.4 present our contribution toward an accurate and efficient application to structural acoustic problems that depend nonlinearly on the eigenfrequency.

General outline of the nonlinear FEAST algorithm
The nonlinear EVPs (7) and (8) are replaced by T(̃)x = 0, T(̃) ∈ C n×n , x ∈ C n for the sake of readability. We are interested in a given number of eigenpairs (̃, x) whose eigenvalues̃lie in a predefined region  that is enclosed by a contour  in the complex plane. Since the imaginary parts of the complex eigenvalues-which correspond to radiation damping-are typically small, we use elliptic contours that have their major axis aligned with the real axis. The intersections of the ellipse with the real axis are the lower and upper bounds of the frequency range of interest. An exemplary ellipse is shown in Figure 1 and is expressed by where Similar to its linear counterpart in Equation (9), NLFEAST uses a Rayleigh-Ritz projection to reduce the system dimension in each iteration. The resulting projected EVP is still nonlinear but of significantly reduced dimension m ≪ n. The subspace dimension m is fixed throughout the iterations and should be at least as large as the number of expected eigenvalues inside  to ensure that all desired eigenvalues are captured. In our case, the number of eigenvalues can be either estimated a priori by an empirical formula, 40 or by a preceding in vacuo analysis-that is, by solving the linear EVP (12) is solved, the m eigenvalues lying closest to the center of  are determined by evaluating the scaled distance Those eigenvalues̃1, … ,̃m and associated eigenvectors y 1 , … , y m are stored in = diag (̃1, … ,̃m), ∈ C m×m and X = [ , X ∈ C n×m , respectively. They are used to update the search subspace Q via F I G U R E 1 Elliptic contour  enclosing a complex domain , in which the desired eigenvalues lie in which the residual matrix B ∈ C n×m is obtained by The integrand in Equation (14) corresponds to the update rule for eigenvectors in the residual inverse iteration method with multiple shifts. 30,39 Its Cauchy integral representation (14) then yields an approximate spectral projection onto the eigenspace corresponding to the eigenvalues in .
With the definition of the contour (11) at hand, the integral can be evaluated using the trapezoidal rule, that is, where N denotes the number of integration points on the contour (cf. Figure 1). For an elliptic contour, the integration points and weights can, for example, be chosen as The general procedure for solving structural acoustic EVPs with NLFEAST is summarized in Algorithm 1. Convergence of the algorithm is monitored by evaluating the maximum residual among the eigenpairs whose eigenvalues lie inside the contour, that is, Alternatively, the accuracy of an eigenpair can be assessed by computing the relative residual by, for example, which requires estimating the 1-norm of the system matrix. 41 Other matrix norms are also possible such as the ∞− or Frobenius norm. Once the residual stagnates, the algorithm is terminated. If the accuracy is then deemed to be insufficient, the whole procedure can be repeated with the computed eigenvectors as new initial guess and by incorporating a more accurate solution strategy for the projected EVP in Equation (12). The latter is discussed in Section 3.3, where we will see that the accuracy can be increased while retaining information from the previous run. The rate of convergence of the algorithm is mainly driven by the accuracy of the numerical integration in line 15, and thus by the number of contour points N. On the other hand, N also determines the number of linear systems T(z j )Y = B that need to be solved in each iteration in line 15 accounting for the main computational effort. When the computations for each contour point are executed in parallel, N could be chosen according to the number of available computing nodes. The BE matrices associated with T(z j ) are stored in memory in order to avoid re-assembly in each iteration. In this respect, the use of data-sparse formats 42 or model order reduction of BE matrices 12 can alleviate the memory requirements. The initial guess of eigenvectors X (0) ∈ C n×m can either be random, or in the case of moderately strong structural acoustic interaction, X (0) can be chosen as the in vacuo modes of the structure. If one is then interested in the complex eigenfrequencies corresponding to those particular modes, the eigenvalue selection in lines 10 and 18 could be made based on a modal assurance criterion. 43 Regarding the computational efficiency, several issues need to addressed when applying Algorithm 1 to structural acoustic EVPs (7) and (8): • Firstly, the repeated evaluation of the residuals in Equation (15) requires assembling the frequency-dependent BE matrices at m eigenvalues in each iteration.
• Secondly, the projected EVP in Equation (12) is still nonlinear and requires a tailored solution scheme.

Rational approximation of frequency-dependent boundary element matrices
The repeated evaluation of the residuals in Equation (15) deteriorates the numerical efficiency of NLFEAST when the frequency-dependent BE matrices H(̃) and G(̃) are assembled explicitly for all intermediate (not yet converged) eigenvalues in each iteration. Several approaches for accelerating the setup of BE matrices have been proposed in the literature. 8 Typically, polynomial approximations based on frequency sampling at Chebychev nodes 24 or truncated Taylor series 20 are employed in this regard. However when evaluating at complex eigenvalues̃that lie apart from the real axis, the quality of those polynomial approximations quickly deteriorates. Therefore, we will use a Cauchy integral representation of the BE matrices which-after discretization with the trapezoidal rule-yields a rational approximation. For a detailed discussion on different approximations inside complex regions, we refer to the review by Austin et al. 44 The approximate BE matrices read where the scalar-valued rational functions v j are given by The error introduced by these rational approximations depends on the fluctuation of the entries of H(̃) and G(̃) with respect to the frequency. For example, the relative error in H(̃) can be assessed by where || ⋅ || F denotes the Frobenius norm. However, an evaluation of Equation (23) is generally not feasible and an alternative way is required for assessing the approximations given in Equations (21) and (22). For this purpose, we first notice that the frequency dependence of a BE matrix entry can be generally characterized by a product of a monomial and an exponential function, 45 that is, with the speed of sound c and the nonnegative integer exponent p. Given the frequency dependence of the single and double layer potential in the collocation BE formulation, p = 1 is a suitable choice. Further, the spatial distance r should be chosen as the largest distance between two points on the discretized boundary. Then, an error estimate can be obtained from In practice, h (̃) could be probed at a few points inside the complex domain  and the number of contour points N could be increased until an acceptable accuracy is achieved. We will see in Section 4.1 that an error estimate based on Equation (25) is in reasonable accordance with the relative error given by Equation (23).
For the convenient choice of using the same elliptic contour for the rational approximation as for the computation of the search subspace in Equation (16), the setup of the rational approximation does not require any additional numerical effort. The BE matrices H(z j ) and G(z j ) are readily available and the contour points z j and associated weights w j are identical to those in Equations (17) and (18). In this case, Equation (16) is a trapezoidal rule approximation of the Cauchy integral representation of eigenfunctions, which in turn includes BE matrices with the same underlying trapezoidal rule approximation. Therefore, we expect that the approximation of BE matrices do not further impair the convergence rate of NLFEAST nor the accuracy of the final eigenpairs. Finally, we note that some intermediate eigenvalues may also lie outside of the contour, particularly in the first couple of iterations of NLFEAST. Although the rational approximations (21) yield exponential convergence only inside the contour, 46 they are still used in these cases to evaluate the residual (15).

Solution of the projected nonlinear eigenvalue problem
A tailored procedure for solving the projected (but still nonlinear) EVP (12) is crucial, since it determines the final accuracy of the eigenpairs obtained by NLFEAST. The dimension of the EVP (12) is rather small, and hence it would admit the application of various nonlinear eigensolvers for dense systems that are available in the mathematical literature. 47,48 However, with regard to computational efficiency, we would favor an eigensolver that requires minimal additional effort in each NLFEAST iteration and in particular, avoids repeated assembly of the original system matrix. In this respect, an obvious choice would be to simply employ a rational approximation as has already been done for the acceleration of the residual computation in Equation (15). However, instead of only approximating the BE matrices as given in Equation (21), an approximation of the entire coupled system matrix T(̃) is necessary in order to obtain a rational EVP that can be subsequently linearized. The rational approximation of the projected EVP (12) is obtained by quadrature of the Cauchy integral representation of Q H T(̃)Q and reads where N p denotes the number of quadrature points. Again, for the convenient choice of the same underlying elliptic contour and trapezoidal rule as for the computation of the projection matrix in NLFEAST, N p = N holds. In that case, the quadrature points for the rational approximation of the EVP are identical to the contour points in Equation (17) (i.e., z j = z j ), and the BE matrices at the contour points can be reused. For solving the rational EVP (26), we follow the approach of El-Guide et al. 32 and reformulate the problem as a generalized EVP w = w (27) with identical eigenvalues̃= . The block matrices in Equation (27) are defined as The main issue with the linearization of EVPs is the significant increase of the system dimension. El-Guide et al. 32 address this issue by subspace iteration in which the columns of the projection matrix are individually computed by inverse power iteration. In our case, the inflation of the EVP associated with the linearization is mitigated by the projection (12), where the projection matrix (16) and the rational approximation of the EVP (26) can be based on the same contour points. The computational effort for solving the linearized EVP is reduced to an order of  ( (mN p ) 3 ) . We will later see that a small number N p is sufficient to achieve a high accuracy and hence, the solution of the reduced nonlinear EVP in Equation (12) only marginally contributes to the overall computational time for the NLFEAST algorithm. Lastly, we note that the final accuracy of the eigenpairs could be improved by adaptively increasing N p in subsequent runs of NLFEAST by adding new quadrature pointsẑ j in-between the previous ones.

Efficient solution of coupled systems with multiple right-hand sides
The main computational effort of the NLFEAST algorithm stems from the solution of the N linear systems in each iteration in order to update the projection matrix via Equation (16). For a discussion on efficient solution strategies for Equation (30), we need to distinguish between the two structural acoustic formulations given by Equations (7) and (8).
When forming the Schur complement of K − z 2 j M with respect to coupled block matrix in (4), the system matrices read whereas using the Schur complement of H(z j ) yields The system in Equation (31) admits efficient application of an iterative solver to Equation (30) once a factorization of the sparse and symmetric matrix K − z 2 j M is stored. Linear systems with multiple right-hand sides such as Equation (30) are usually solved by successively applying standard iterative schemes to obtain the solutions to each forcing vector individually. More efficient approaches such as subspace recycling 49 or block Krylov solvers 50 enable to carry over information among the solutions to different forcing vectors. In this work, we will employ a block variant of the generalized minimal residual (GMRes) method. 51,52 While block Krylov solvers have been successfully applied to acoustic BE equations, 53 to the best of our knowledge, the present contribution is the first application to coupled structural acoustic systems with multiple right-hand sides.
Block Krylov solvers generally work by refining the solution Y k in each iteration k such that holds, where Y 0 and R 0 = B − T(z j )Y 0 denote the initial guess and the corresponding initial residual, respectively. The block Krylov subspace  k is given by In the case of block GMRes, the orthonormal basis vectors spanning the subspace  k are computed via a block Arnoldi process. The resulting block subspace is of much larger dimension compared with the subspace in standard GMRes, and hence a faster convergence can be expected.
Having discussed the solution when using the formulation (31), we now turn our attention to the Schur complement of the acoustic subsystem in Equation (32). Here, an iterative solution would require evaluation of H −1 (z j )G(z j )C fs Y (k) j in each iteration k. In order to avoid two nested iterations, we employ a model order reduction 54 to the frequency-independent structural subsystem and thereby enable a direct solution of Equation (30). The reduced system (32) reads in which V ∈ R n s ×q could be a Krylov subspace or spanned by the in vacuo modes of the structure. In this work, we choose a Krylov basis, which is computed once for the whole frequency range of interest. 55 Once the Krylov basis is established, the reduced system T r (z j ) can be explicitly formed by evaluating H −1 (z j )G(z j )C fs V using the above described block GMRes. The N reduced matrices T r (z j ) can be saved and reused in all subsequent NLFEAST iterations. Hence, the computational effort of solving the high-fidelity system in Equation (30) in each iteration boils down to the effort of computing V and reducing the coupled system matrices at each contour point only once via Equation (35). The actual solution in each iteration requires negligible computational effort. The memory requirements to store the system matrices are reduced as well.

FIRST NUMERICAL EXAMPLE: MODAL ANALYSIS OF A SUBMERGED SPHERICAL SHELL
As a first example, we consider a spherical shell submerged in water. Its geometrical and material parameters are listed in Table 1. A finite element mesh with 864 eight-noded quadrilateral finite elements based on the Reissner Mindlin theory is used to discretize the sphere. This corresponds to twelve elements on a ∕2 arc and 15,564 displacement DOFs. The acoustic field is represented by a conforming BE mesh with discontinuous bilinear interpolation functions and 3456 pressure DOFs. We first illustrate the application of NLFEAST to perform a modal analysis in Section 4.1, and then compare the performance to a contour integral method in Section 4.2.

Eigenvalue solution using NLFEAST
In this section, we use NLFEAST to compute the eigenpairs of the submerged sphere in the frequency range from f min = 34 Hz to f max = 82 Hz. As illustrated in Figure 2, these frequency bounds are used to define the elliptic contour. A treatment for the nonuniqueness problem in BEM is not required in this frequency range, but would be needed above the first interior Dirichlet eigenfrequency of approximately 148 Hz. The Burton-Miller formulation 56 would resolve this issue, however care has to be taken since the spurious eigenfrequencies are not eliminated altogether but only moved away from the real axis. 57 They could be distinguished from the actual ones based on the sign of the imaginary parts, given a proper choice of the coupling parameter in the Burton-Miller formulation. 58 Based on the analytical solution, 59 we expect three distinct eigenvalues in the considered frequency range with geometric multiplicities of five, seven, and nine, respectively. The aspect ratio of the elliptic contour is set to = 0.05. The dimension of the search space in NLFEAST is set to m = 21, which is sufficient to capture all eigenvalues. The initial eigenvectors X (0) are chosen randomly. The convergence criterion for block GMRes is defined such that relative residuals smaller than gmres = 10 −12 are achieved for all right-hand sides. In this example, we use the Schur complement of the structural subsystem resulting in the EVP formulation given by Equation (8). Figure 3 displays the convergence behavior of NLFEAST for different numbers of contour points N, and a fixed number of N p = 16 for the solution of the projected EVP (26) by rational approximation. The accuracy is assessed by the maximum residual defined in Equation (19), in which T(̃i) is evaluated explicitly without the approximations in Equations (21) and (22). In this example, all m = 21 eigenvalues were found regardless of the number of contour points in the first iteration after the initialization step. No spurious eigenvalues occurred in any of the test cases with NLFEAST. As expected, the number of contour points N (and thus the quality of the projection matrix Q) determines the convergence rate of NLFEAST, whereas the final accuracy of the eigenpairs is unaffected by this choice. Even when using only N = 6 contour points, a maximum residual of max = 4 ⋅ 10 −10 is achieved after a couple of iterations. Note that in this example, the 1-norm of the system matrix is of order (10 0 ) in the whole frequency range on interest and hence, the absolute and relative residuals are roughly in the same order of magnitude, cf. Equations (19) and (20). For example, the above mentioned maximum residual of max = 4 ⋅ 10 −10 corresponds to a relative residual of rel = 6 ⋅ 10 −11 . The final accuracy of the eigenpairs is mainly controlled by the accuracy with which the projected EVP is solved. The latter is in turn dependent on the number of quadrature points N p for the rational approximation in Equation (26). This is illustrated in Figure 4, where now max is displayed for different values of N p and a fixed number of contour points N = 12. We observe that the final residuals stagnate at different values depending on the number of N p . For example, when using N p = 6, the maximum residual stagnates at max = 7 ⋅ 10 −4 . Rational approximations with more than N p = 16 did not further improve the results beyond max = 4 ⋅ 10 −10 .
In Section 3.2, we have proposed rational approximations of the BE matrices (21) in order to avoid their explicit assembly in each iteration. Figure 5 displays thereby introduced relative error H (̃) given by Equation (23) when using N = 12 contour points. The plot is generated by 1683 uniformly distributed points in the complex plane and linear interpolation in-between them. We notice that N = 12 contour points are sufficient to achieve errors H (̃) of order (10 −10 ) inside the elliptic contour. The error estimator for an matrix entry given in Equation (25) yields a maximum value of h (̃) = 6 ⋅ 10 −9 , which is in reasonable accordance to the actual error shown in Figure 5. The influence of the rational approximations of the BE matrices on the convergence of NLFEAST is examined in Figure 6. The residuals with and without approximation are compared for N = 8 contour points and a varying number of points N p for the rational approximation. We notice that the rational approximations only have a negligible impact on the convergence rate as well as on the final accuracy-at least up to N p = 12. When using N p = 16, the rational approximations prevent the eigenpair from reaching the same final accuracy as is reached when the BE matrices are explicitly formed. It stagnates at max = 3 ⋅ 10 −9 while without the rational approximations, max = 4 ⋅ 10 −10 is achieved. However, in view of the significant saving of computational effort, this minor deterioration of accuracy is deemed acceptable-also because other sources of inaccuracies such as discretization errors typically have a larger impact.
As discussed in Section 3.4, the computational effort is mainly associated with the solution of linear systems (30) with multiple right-hand sides for updating the projection matrix via Equation (16). The performances of conventional and block GMRes are compared with each other in Figure 7 by counting the equivalent matrix vector multiplications. When using conventional GMRes, the linear systems are solved individually for each right-hand side. The corresponding relative residuals are plotted one after the other in Figure 7 yielding the apparent sawtooth shape (i.e., each tooth corresponds to one right-hand side). The complete solution requires 822 matrix vector multiplications within GMRes corresponding to an average of 39 multiplications for the solution to each right-hand side.
In contrast, block GMRes solves the linear system (30) for all right-hand sides simultaneously by performing one matrix multiplication in each iteration. For the sake of comparability with conventional GMRes, each of those operations are counted as m = 21 equivalent matrix vector multiplications. The dashed red curve represents the relative residual associated with the first right-hand side. It exhibits a significantly steeper slope than conventional GMRes due to the larger dimension of the underlying (block) Krylov subspace. This also leads to a faster overall convergence which is illustrated by the solid red curve. The latter displays the maximum relative residual among all right-hand sides. In total, block GMRes requires 15 iterations and a total of 315 equivalent matrix vector multiplications. In terms of computational F I G U R E 4 Maximum residual among the eigenpairs whose eigenvalues lie inside the contour for different numbers N p for solving the projected EVP. The number of contour points is set to N = 16 time, block GMRes additionally benefits from efficient memory usage due to less frequent access to the stored matrices. In this example, block GMRes yields a speed up of more than eight times compared with conventional GMRes.

Comparison to contour integral methods
In the following, we briefly compare the performance of NLFEAST to two other nonlinear eigensolvers.

Comparison to the block Sakurai Sugiura method
The block Sakurai Sugiura method (block SS) 26 is a contour integral method based on resolvent moments and essentially transforms the nonlinear EVP to a generalized EVP involving block Hankel matrices. The block Hankel matrices contain moments of the resolvent of the system matrix up to an user-defined degree, whereas the moment matrices are computed with respect to a number L of random source vectors via complex contour integration. The algorithmic details of block SS can be found in the literature. 21,26,40 The computational efforts of both NLFEAST and block SS mainly comprise the assembly of BE coefficient matrices and the solution of linear systems at the contour points. The latter aspect is studied in Figure 8, which shows the maximum residual over the number of solved linear systems. NLFEAST is executed with N p = 16 and different numbers N of contour points. As we have already seen in the previous section, the residual in NLFEAST decreases monotonically until stagnation. In the case of block SS, we also observe that the residuals generally decrease with an increasing N (and thus increasing number of linear systems). However, there are also some apparent fluctuations in Figure 8, even though the random source vectors are kept the same among all runs of block SS for the sake of comparability. In other words, the accuracy of the eigenpairs computed by block SS may actually deteriorate in some cases even when more integration points are used for the projection. Moreover, we note that spurious eigenvalues occurred in all our test cases with block SS although we perform a singular value decomposition and truncation of the Hankel matrices as suggested by Asakura et al. 26 These spurious eigenvalues associated with large residuals are not included in Figure 8. Figure 8 also indicates the number of linear systems that need to be solved in order to reach the final (i.e., highest possible) accuracy. For the sake of comparability, we define a threshold of max = 10 −8 . When using block SS with L = 21, a total of N = 32 contour points (i.e., the solution of 32 linear systems) are required in order reach that threshold. Compared with that, NLFEAST is an iterative scheme that solves N linear systems in each iteration. For example with N = 8, the threshold is reached after five iterations, which corresponds to a total number of 40 linear systems. Similarly, when using N = N p = 16 contour points, only two iterations are required summing up to a total of 32 linear systems. Table 2 lists the associated computational times for the solution of those systems as well as the times that are required for the assembly of the BE matrices. Regarding the solution of linear systems, NLFEAST benefits from the reuse of the LU factorizations of K − z 2 j M over the course of the iterations. Hence, given a certain number of linear systems and right-hand sides, the solution time of NLFEAST is slightly shorter. Note that we also use block GMRes for solving the F I G U R E 8 Number of linear systems that need to be solved to reach a threshold of ε max = 10 −8 . All computations with NLFEAST are conducted using N p = 16. The Hankel matrices in block SS contain moments up to a degree of 3. The associated computational times are listed in Note: Associated convergence behavior is plotted in Figure 8. a The total time does not refer to the total wall clock time for the execution of the algorithms but only includes matrix assembly, rational approximation, and solution of linear systems. The other (marginal) contributions, which are mostly identical in all test cases, are neglected.

TA B L E 2 Computational effort that
is required to solve the benchmark EVP with ε max ≤ 10 −8 linear systems arising in block SS. But actually, the major saving does not stem from the decrease of solution time but from the fact that NLFEAST requires significantly less contour points in order to achieve the highest possible accuracy (with a minimum number of N = N p = 13). As a result of this, fewer BE matrices need to be explicitly assembled. In this example, the assembly of H(z j ) and G(z j ) at a single contour point took 57 s while the time for computing the rational approximation via Equation (21) only marginally contributed to the total time. Summing up, depending on the respective solver parameters, NLFEAST required 22 % to 61 % less computational time then block SS to achieve the final accuracy. Of course these numbers may change depending on the dimension of the problem and in the case of parallel computing, also on the number of available computing nodes.

Comparison to the Sakurai Sugiura method with Rayleigh Ritz projection
The ill-conditioning of the Hankel matrices in block SS has been addressed in a paper by Yokota and Sakurai,28 in which the authors suggest a Rayleigh Ritz procedure to extract the eigenvalues from the resolvent moments based subspace. This procedure is known as the Sakurai Sugiura method with Rayleigh Ritz projection (SSRR). While SSRR is expected to give more accurate results than block SS, it needs to be accompanied by a nonlinear eigensolver for the projected (but still nonlinear) EVP. This nonlinear EVP of small size can be conveniently solved by block SS with a large number L of source vectors without impairing the computational efficiency. A similar strategy has been employed for the solution of acoustic EVPs based on BEM, 29 however with the difference that the Rayleigh Ritz projector is built by sampling the resolvent matrix.
Our experience with SSRR is in line with those reported in previous works 28,29 and indicates that the method is computationally more efficient and robust than block SS. For example, with L = 21 and a degree of moments of 3, SSRR only requires N = 18 contour points to reach the threshold of max = 10 −8 (recall that block SS required N = 32 and NLFEAST N = N p = 13). Compared with NLFEAST, fewer linear systems need to be solved in SSRR, but still, NLFEAST requires less computational time in this example since fewer BE matrices need to be explicitly assembled. Further, we note that although the occurrence of spurious eigenvalues is effectively mitigated in SSRR, we still encountered them when solving the projected EVP via block SS, even with a very small degree of moments and large number L of source vectors. In contrast, no spurious eigenvalues occur when combining SSRR with the rational approximation method presented in Section 3.3.

SECOND NUMERICAL EXAMPLE: QUANTIFICATION OF MODAL RADIATION DAMPING IN A MUSICAL BELL
A bell excited by a clapper radiates sound by an ensemble of its fundamental modes, where each of them is characterized by an initial amplitude, eigenfrequency and a decay rate. While the eigenfrequencies can be anticipated by an in vacuo modal analysis of the bell, the decay characteristics, which make up for the desired sound of the bell, depend on the acoustic radiation damping of each individual mode. 60 It is the yearlong experience of bell founders that enables them to manually tune the bell after casting.
In this second numerical example, we apply the NLFEAST algorithm in order to determine radiation damping values associated with modes of a bronze musical bell. The main geometrical and material parameters of the bell are listed in Table 3. The bell is discretized using 20-noded hexahedral solid finite elements resulting in a total of 26,769 displacement DOFs. The respective stiffness and mass matrices are extracted from ANSYS. Figure 9 shows the mesh of the bell. The outer surfaces of the solid elements are coupled to a conforming acoustic BE mesh consisting of 2560 bilinear discontinuous elements with a total of 10,240 pressure DOFs.
The modal analysis of the bell in interaction with surrounding air is performed using an elliptic contour with f min = 80 Hz, f max = 390 Hz, = 0.05, and N = N p = 48. A treatment for the nonuniqueness problem in BEM is not required. From a preceding in vacuo analysis, we expect m = 10 eigenpairs in the considered frequency range. The in vacuo modes are also used as the initial guess of eigenvectors X (0) in NLFEAST. Due to the light fluid loading, the coupled system (4) admits accurate model order reduction solely based on the structural subsystem. Hence, we form the Schur complement of the acoustic subsystem resulting in the EVP formulation given by Equation (7). The system matrices at the contour points are reduced by a Krylov subspace model order reduction as described in Section 3.4. The Krylov basis is computed using a random initial vector at a single expansion point of 250 Hz. A dimension of q = 40 for the Krylov basis is required in order to obtain all eigenpairs in this example. Hence, the computational effort of the modal analysis basically comprises the setup of the BE matrices and solution of a single linear system H(z)X = G(z)C fs V with q = 40 right-hand sides for each contour point. Compared with that, using the EVP formulation (8) based on the Schur complement of the structural subsystem would necessitate solution of a linear system with multiple right-hand sides for each contour point every time when the projection matrix needs to be updated. The memory requirements for storing the system matrices of reduced dimension q is negligible. In contrast, using the EVP formulation (8) would require storage of several BE matrices of original dimension. Moreover, using the latter approach, we were not able to find accurate eigenpairs, which may be attributed to ill-conditioning of the FE matrices K − 2 M. In this example, NLFEAST requires two iterations to reach the final accuracy.
The first five distinct modes of the bell and the associated eigenfrequencies are shown in Figure 10. Note that due to rotational symmetry, all modes occur with a geometric multiplicity of two. While the eigenfrequencies corresponding to the real parts of the eigenvalues remain almost unaffected by air loading, the imaginary parts of the eigenvalues characterize the extent of acoustic radiation damping. Radiation damping corresponding to an individual mode can be expressed where Im (̃i) is negative due to the harmonic time dependency e −i t . Beside modal loss factors, harmonic loss factors are an alternative way to quantify radiation damping and can be computed by a frequencywise response analysis. 2 They can serve us to verify the modal loss factors and thus the accuracy of the computed eigenvalues. Figure 11 shows the modal loss factors obtained by NLFEAST as well as harmonic loss factors of the bell subject to point force excitation. The latter are obtained by computing the radiated sound power in a standard frequency sweep analysis, which can be performed either based on the Schur complement of the acoustic subsystem in conjunction with the above described Krylov subspace model order reduction or based on the Schur complement of the structural subsystem without model order reduction. The maximum relative difference in radiated sound power between these two approaches is 2.5 ⋅ 10 −5 which undermines that the coupled system admits accurate model order reduction solely based on the structural subsystem. Figure 11 also indicates that although the imaginary parts of the eigenvalues are very small, the resulting modal loss factors agree well with the harmonic ones at the respective eigenfrequencies.

SUMMARY AND CONCLUSION
We have proposed a subspace iteration method for the solution of nonlinear structural acoustic EVPs. At its core, the method is based on the nonlinear FEAST algorithm and essentially works by iteratively refining a projection matrix that is obtained by applying a spectral projector. We have shown that in the context of structural acoustic problems involving complex eigenvalues, the information for computing the projection matrix can be recycled in two ways and thus giving rise to a tailored eigensolver. Firstly, the integration points are reused for a rational approximation of the projected but still nonlinear EVP. This enables an accurate computation of eigenpairs at a negligible additional cost without inducing spurious eigenvalues. Secondly, the integration points are also reused to approximate the BE matrices at intermediate eigenvalues thereby avoiding their assembly in each FEAST iteration. Only a few points are sufficient to achieve small errors in the BE matrices, and hence to retain the convergence rate of FEAST while significantly reducing the numerical effort.
In consequence, the actual computation of the projection matrix remains the only computationally expensive task required for solving the EVP. In this respect, we provide two tailored approaches for updating the projection matrix. In the case of heavy fluid loading, we form the Schur complement of the structural subsystem and employ block GMRes for solving the underlying coupled linear systems with multiple right-hand sides. Compared with standard GMRes, we have achieved a speed-up of eight times in our test cases. In the case of light fluid loading, we exploit the fact that the coupled system admits accurate model order reduction solely based on the structural subsystem. By forming the Schur complement of the acoustic subsystem and computing the reduced coupled matrices again by block GMRes, we are able to significantly reduce both the memory requirements and the effort for updating the projection matrix.
Applications to a submerged spherical shell and to a musical bell have verified the proposed method. Comparison to resolvent moments based contour integral methods indicates two major advantages of the proposed eigensolver. While both methods yield similar final (i.e., highest possible) accuracies, the proposed method requires a significantly smaller number of contour points. Hence, less BE matrices need to be explicitly assembled. Moreover, in all our test cases, the proposed method yielded monotonous convergence without inducing spurious eigenvalues inside the contour. Finally, the application to a bell in interaction with air has shown that accurate modal radiation damping values are obtained despite the very weak acoustic loading.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.