Automatic CHIEF Point Selection for Finite Element–Boundary Element Acoustic Backscattering

: Computing the backscattering of harmonic acoustic waves from underwater elastic targets of arbitrary shapes is a challenging problem of considerable practical signiﬁcance. The ﬁnite element method is well suited for the discretization of the target, while the boundary element method addresses the radiation boundary condition at inﬁnity. A disadvantage of the boundary integral method is that it yields non-unique solutions at certain wavenumbers. This failure is associated with the existence of eigensolutions of the Helmholtz equation in the interior of the complement of the ﬂuid domain (acoustic modes). The combined Helmholtz integral equation formulation (CHIEF) credited to Schenk is employed to combine the surface Helmholtz boundary integral with equations of the interior Helmholtz relation written down at selected points within the cavity of the scatterer (i.e., in the complement of the ﬂuid domain).The difﬁculty associated with this approach has always been the lack of guidance on the necessary number of interior points and on their locations. The solution to this problem proposed here is to compute the acoustic modes using the ﬁnite element method to complement of the ﬂuid domain and to identify locations of the peaks.This novel approach aids the decision as to how many points should be employed and where they should be located. Our numerical experiments demonstrate the robustness of the proposed automatic selection of the CHIEF points’ numbers and locations.


Introduction
Consider the classical problem of the scattering of a pressure wave by an arbitrary three-dimensional elastic structure submerged in an acoustic fluid [1,2]. The task of calculating the scattered sound-pressure field is addressed here with a coupled finite element and boundary element computational framework. The finite element method (FEM) is a suitable tool for the simulation of structural vibrations [3]. The FEM can also be employed for the acoustic part, however special schemes such as infinite elements or other types of treatment of the silent boundaries are required for unbounded media. On the other hand, this is where the boundary element method (BEM) can be employed since the radiation condition for exterior domains is automatically satisfied [4].
The integral equation approach results in the surface Helmholtz equation (SHE) [5]. It relates the acoustic pressure on the surface of the structure to the normal surface velocity. A disadvantage of this method is that, like any classical formulations of the Helmholtz equation in the exterior region, whether through the use of single-or double-layer virtual source distributions or via the use of Green's theorem, it breaks down due to the nonuniqueness of their solutions at certain wavenumbers, i.e., when the wavenumber k is an eigenvalue of the Dirichlet problem for the interior region [6][7][8]. This failure is completely non-physical in nature, being entirely due to the integral equation formulation, and is associated with the existence of eigensolutions of the Helmholtz equation in the interior of the wet surface of the scatterer. As cogently elucidated by [9][10][11], the problem is really the ill conditioning of the coefficient matrix in the vicinity of the characteristic frequencies of the interior modal problem.
It is possible to write the equivalent of the SHE at points interior to the outer surface of the scatterer, i.e., in the complement of the fluid domain: the so called Interior Helmholtz Relation (IHR). The solution of the exterior problem should vanish in the interior of the scatterer. The combined Helmholtz integral equation formulation (CHIEF) conceived by Schenck [7,12] uses the numerical stability of the SHE and the uniqueness property of the IHR to form an overdetermined set of linear equations. These equations can then be solved by a least-squares procedure. Schenck proved that at a critical wavenumber, only the required (unique) solution of the SHE will simultaneously satisfy the IHR, provided the IHR is not evaluated on the nodal surface of the relevant interior eigenfunction. This is the aspect of the CHIEF that resulted in a large subsequent volume of work: how to determine how many additional interior points are needed, and where they should be located. A recent search yielded more than 700 citations of the original paper [7] in papers and books [13]. See [14] for an early systematic study. To name a few examples, we refer to SuperCHIEF [15], Block (aka weighted residual) CHIEF [16], HyperCHIEF [17], and IHIEF [18]. A similar approach has recently been developed for the indirect formulation [19]. A related approach due to Panich was developed in [20]. Bartolozzi et al. have recently argued for a random selection as the optimal solution [21]. The CHIEF procedure has been recently adopted in [22]. For a thorough review, refer to [11,23].
Of course, the difficulties associated with CHIEF may be solved, at least to a significant extent, by switching to a different method. The Burton-Miller (BM) method appeared on the heels of CHIEF [24], and many variants followed. A major study of both the CHIEF and the BM methods is available in [23]. Alternative approaches have also appeared, such as the dual-surface method [25]. Refer also to a recent review [26].
As stated already in [14], ". . . The question of the appropriate selection of 'CHIEF points' is a potentially serious shortcoming to CHIEF, as neither the characteristic wavenumbers nor the nodal surfaces are known without solving the related interior problem first"; this is precisely where we take our inspiration. Since we are solving the backscattering problem with a combination of finite element and boundary element methods, we are naturally in possession of a volume mesh of the scatterer. Hence, the quote displayed above becomes a starting point for our solution: we solve the interior modal problem using a finite element mesh. Thus, the wavenumbers of the interior problem become known, and by identifying peaks of the associated mode shapes of the pressure (both positive and negative), we obtain an estimate of both the number and the locations of the CHIEF points.
The structure of this article is as follows. The method is developed in Section 2, where the boundary and the finite element treatments of the exterior acoustic problem and of the motion of the solid scatterer, respectively, are briefly reviewed. We address the coupling of the boundary element method used for the acoustic medium and the finite element method used for the elastic target, and we discuss the suppression of the spurious solutions with the CHIEF method. Then we introduce the main idea: that the finite element model of the interior can yield an approximation of the spectrum of the interior problem. We develop a method of identifying suitable points for the interior Helmholtz relation (IHR). In Section 3 we demonstrate that the proposed method can solve the backscattering problems without loss of accuracy due to ill conditioning. The additional cost of solving the eigenvalue problem is shown to be insignificant. Section 5 summarizes the key takeaways from this work.

Methods
We begin by providing a rudimentary synopsis of the requisite equations. This material is not new, as it already has been the subject of many articles and even textbooks [27].

Exterior Acoustic Problem
For the application of the boundary element technique we employ an integral equation formulation of the exterior Helmholtz problem. (The wide-hat symbol • is used to indicate harmonic amplitude.) For an acoustic wave, p inc (x), disturbs the structure, resulting in a scattered acoustic wave, p scat (x), the integral equation for the total acoustic wave, The middle equations in (1) are known as the surface Helmholtz equation (SHE) and the bottom is known as the interior Helmholtz relation (IHR) [8]. In this work, G(x, y) is the free-space Green's function. There are two unknown quantities in the resulting equation, the pressure and the surface normal displacement.
The second equation of the coupled system is the governing equation of the solid scatterer, which will be now written directly for the steady state of harmonic vibration. We assume that body forces are absent and the vibrating solid is free of displacement constraints: that is, the solid is free-floating in the acoustic medium. We employ here the so called "Voigt" vector notation. Here B is the strain-displacement operator, n is the outer normal, and η is the test function. The boundary conditions assumed on the surfaces of the voids within the solid scatterer (if present) are the absence of tractions (i.e., traction-free surfaces). The weak formulation of the governing equations results in the well-known Galerkin weighted residual equation The numerical solution can now be sought by discretizing the fluid equations using the boundary element method, and by discretizing the equations of the solid using the finite element method. The finite element model is based on the nodally integrated four-node tetrahedra with stabilization, as detailed in Reference [28]. Consequently, the surface of the solid scatterer (i.e., the surface S, and all the surfaces of the internal voids) is tiled with triangular faces. We shall call these faces panels here.
The middle equation in (1) can be written in discrete form using collocation, where each of the triangular panels corresponds to a face of the tetrahedral element employed for the solid scatterer that is exposed to the acoustic medium where where i, j = 1, . . . , N. The objects A and B are often referred to as the double-and singlelayer potential matrices. The vector V n consists of the normal components of the velocity at the centroids of the panels. Both integrands in (4) possess integrable singularities when y → x i , which can be handled by a variety of methods [29]. The finite element model is based on the nodally-integrated four-node tetrahedra with stabilization as detailed in Reference [28]. Energy-sampling stabilization is applied to the nodally integrated elements, effectively removing spurious (unphysical) modes at both lower and at higher ends of the frequency spectrum. The formulation is shown by the numerical inf-sup tests that it is coercive and locking-free, and hence suitable for nearly incompressible materials such as rubber. In order to increase the fixed-mesh accuracy, the two stabilization parameters introduced in [28] have been adjusted to a = 2 and b = 3. The goal is to derive additional accuracy from the nearly equilateral tetrahedral elements, typical of the meshes used in the acoustics simulations. The finite element model discrete equations are derived from a variant of (2), as described in considerable detail in Reference [28], resulting in the well known matrix expression of undamped harmonic forced vibration The mass matrix M and the stiffness matrix K result from the first two terms in (2), and the vector of the nodal loads L is due to the acoustic pressure p on the boundary S, which is the last term in (2). The vector L represents the loading of the scatterer by the total pressure on its surface. Approximating the distribution of the total pressure on the surface allows us to write where G of size (D, N) is the rectangular coupling matrix defined in components as with the notation (K, j) indicating a map from node number K and coordinate direction j to a degree of freedom number (K, j); D is the total number of displacement degrees of freedom; m is the pressure degree of freedom in the boundary element method, and S m is the surface panel m; N is the total number of pressure degrees of freedom (i.e., number of triangle panels); and P is the vector of pressure values on the surface panels. Analogously, the normal velocity of the surface panel, piecewise constant along each surface panel, can be expressed from the vector of the nodal velocities as the surfaceweighted mean of the nodal velocities [12] where D is a diagonal matrix of surface panel areas, Therefore, using V = −iω U, we can express the normal velocity on the panels as Consequently, combining (10) with the equation of motion of the scatterer (5), with the definition of the pressure loads (6), and with the BEM Equation (3), yields or which may be solved for the panel values of total pressure P. It bears emphasis that the matrices A and B, and the vector P inc depend on ω.
The matrix −ω 2 M + K is known as the dynamic stiffness. Note that it depends on the frequency of the excitation, ω. In the applications considered here, i.e., in the socalled "acoustic color" computation of the backscattered target strength as a function of aspect angle and frequency, the scattered pressures need to be calculated for many values of the frequency ω. In order to avoid the large cost of inverting (factorizing, really) the large dynamic stiffness matrix, a modal expansion may be introduced [12]. Using the truncated rectangular matrix of eigenvectors Φ r , where each of the r columns represents an eigenvector, we can write approximately Here, the matrix subscript r indicates that only r eigenvectors (columns) have been retained in the matrix Φ r . Therefore (12) may be approximated with the computationally much more tractable expression Equation (14) needs to be solved for the vector of the pressures, P.

The Non-Uniqueness of Solutions: CHIEF
The original boundary value problem has a unique solution; on the other hand, the linear algebraic Equation (14) corresponding to the discrete integral equation may not be uniquely solvable at certain values of the wavenumber corresponding to the adjoint interior modal problem (i.e., solutions of acoustic modes in the complement to the fluid domain with a homogeneous boundary condition).While the continuous integral equation yields non-unique solution only at a discrete set of wavenumbers, the approximating linear algebraic equations become ill-conditioned when the wavenumber is merely in the neighborhood of a critical value [11]. The two main techniques that have been applied to exterior acoustic problems in order to resolve the non-uniqueness problem are CHIEF and Burton-Miller methods, as discussed at length in the Introduction.
In CHIEF, the SHEs are augmented by imposing some IHRs (refer to (1)) at properly selected interior points. The uniqueness requires that the interior points (CHIEF points) do not lie on internal nodal surfaces of the interior modes. The resulting overdetermined system can then be solved using the least-squares [7] or Lagrange multipliers [14] techniques. Various modifications and improvements have been proposed, mainly to ensure uniqueness (references have been provided in the Introduction).
The crucial aspect of CHIEF is how to select both the number and the locations of the interior points for which the IHRs are added to the system of equations. Usually the hope is that by generating some points within the interior, at least one of those points will manage not to lie on a zero isosurface of the interior acoustic modes; such points would not be effective in constraining the solution. Points which manage to make the solution unique have been labeled "good" in [14]. The more interior points are generated, preferably at random locations, the more likely it is that some of them will be "good." Denoting the coefficient matrix of Equation (14) as H, we can formulate the classical CHIEF equations as [12,23] where the subscript • int denotes the interior points of the IHR. This equation with a rectangular coefficient matrix is best solved as a classical least squares problem, perhaps with the QR decomposition. A symmetric formulation can be also derived [14,23,30,31], for instance, as the "unconjugated" form [23] H, Here Λ int is a vector of Lagrange multipliers, whose values measure the violation of (1) at the interior points. The symmetric formulation leads to a faster solution, given that only an LU decomposition is required. On the other hand, there is some evidence [23] that the least-squares solution of (15) is more successful at eliminating the irregular frequencies than (16). Some evidence supporting this contention will be provided in the section on examples.

Automatic Selection of CHIEF Points
Many (if not all) of the papers and book chapters citing the original 1968 paper [7] (and there are hundreds) make a statement to the effect that ". . . the appropriate selection of 'CHIEF points' is a potentially serious shortcoming to CHIEF, as neither the characteristic wavenumbers, nor the nodal surfaces, are known without solving the related interior problem first". However, we are solving the backscattering problem with a combination of finite element and boundary element methods, and so a volume mesh of the scatterer is available to us. Hence, we simply propose to solve the interior modal problem using a finite element mesh, and to use the extracted information to place the IHR points in the "good" locations, away from the nodal surfaces.
A straightforward Galerkin discretization of the acoustic equation using the traditional isoparametric formulation on four-node tetrahedra [32] leads to a moderately sized eigenvalue problem that can typically be solved in seconds: A total of 70 modes within a sphere with almost 20,000 nodes can be calculated in under 2 seconds of wall-clock time on a modest laptop; for a larger model of the sphere with 80,000 nodes, the interior-mode problem can be solved for 70 modes under 25 seconds, while the solution of the boundary-element problem for a single frequency takes 1020 seconds. With the solution of the modal problem in hand, we locate the propitious locations of the interior points for the IHR equations using the Algorithm 1. Note that the number of modes n m can be tailored to the problem at hand: if we know the lower and upper limit of the frequency sweep, we can work only with modes within that range.

Algorithm 1 Algorithm Find-CHIEF-points
Require: Solution of the interior modal problem for n m number of modes Make the set of interior points, I, initially empty for i = 1, . . . , n m do For all pressure modes Compute the magnitude of the pressure mode: p a,i = |p i |. Compute threshold p thr = α max k p a,k . Collect finite elements connected to nodes j with nodal dofs p a,j > p thr into a set E. while E is not empty do Mark the subset S of E such that all such elements are connected together Find all nodes connected by elements from the set S, adding them into a set N. Add the node with the largest value p a,j within the set N to the set I. Remove the elements of S from E end while end for return I The goal is to identify peaks of the magnitude of the pressure mode: these are definitely points away from the nodal isosurfaces of the pressure mode, and hence "good" point locations. The threshold pressure is determined by a heuristic. Either: α = min(0.8 max k p a,k , 2 mean p a,k ) max k p a,k , which enables selection of multiple points for each mode, or in which case only a single point per mode (the point with the largest mode amplitude) is selected. The task of locating the "peaks" is reminiscent of finding local maxima of an arbitrary function. Because of the potential expense of such an operation, we choose to employ the above described heuristic. This typically indicates multiple "good" points per mode. Compare with Figure 1: the blue and red surfaces show regions in which the acoustic mode results in negative or positive pressures of significant amplitude. The "good" interior points should be located inside these surfaces (ideally at the peaks).
There is one paper where a distantly related idea was employed [25]: ill-conditioning of the coefficient matrix was monitored, and, when detected, the internal acoustic field was evaluated along the axis of symmetry (obviously, the technique was only applicable to solids of revolution) and the "good " point was located at the maximum of the internal acoustic field.

Results
The computational experiments described below were implemented in the Julia programming language [33,34] in the framework of the FinEtools.jl Julia package [35]. The experiments were executed on a Windows 10 machine with 2× Intel (R) Xeon(R) E5-2670 at 2.6 GHz, each with 8 cores, and with a total of 256 GB of memory, L1 cache of 2 × 8 × 32 KiB, L2 cache of 2 × 8 × 256 KiB, and L3 cache of 20,480 KiB.
The performance of the proposed point selection algorithm is illustrated here on an example of acoustic backscattering: solid steel sphere of 0.5 m radius in water. The properties of the scatterer were: Young's modulus E = 205,000 MPa, Poisson ratio ν = 0.28, mass density ρ = 7850 kg/m 3 ; those of the acoustic medium were: mass density ρ w = 1000 kg/m 3 and sound speed c w = 1500 m/s. The sweep in the frequency domain consisted of 801 frequencies, from 1000 Hz to 5000 Hz by 5 Hz steps. The mesh of the sphere consisted of 4-node tetrahedra with a prescribed edge length of 0.025 m and 19,693 nodes. The modal matrix of (14) consisted of 300 modes, computed to default accuracy with an Arnoldi-iteration-based method [36]. The largest natural frequency of the sphere was approximately 10.8 kHz. The surface mesh consisted of 10,228 triangles. Each of the solutions of Equation (15) or Equation (16) typically took on the order of one minute of wall-clock time or more, and hence the outlay of a few seconds to solve the interior acoustic-mode problem was negligible.
The spherical volume has interior Dirichlet modal frequencies as shown in Table 1: the frequencies are reported according to their multiplicity. In the range of interest there are eleven distinct frequencies, for a total of 53 modes, when all multiplicities are accounted for. Refer also to Figure 1 for a graphical representation of four particular interior modes. Mode 4 has a multiplicity of 3, mode 9 has a multiplicity of 5, mode 16 has a multiplicity of 7, and mode 25 has a multiplicity of 9. Of course, the modes computed on the finite element mesh displayed some discretization errors: for instance, the computed frequency of mode 25 was higher by approximately 2.5% than the analytical value of 3.90688 kHz of Table 1. 3.33649 × 10 3 4.97380 × 10 3 6.54032 × 10 3 9 3.90688 × 10 3 5.58868 × 10 3 7.18091 × 10 3 11 4.46707 × 10 3 6.19106 × 10 3 7.80880 × 10 3 In the graphs below, the backscattering amplitude as a function of frequency is presented such that the solid black line is the reference partial wave solution [37] and the thin colored line represents the coupled finite and boundary element solution.
According to Table 1, we should expect the method based solely on the SHE equation to experience loss of accuracy due to ill-conditioning around frequencies lower than 5 kHz. As explained above, in the usual setup of the CHIEF, we do not know how many points we should use and where they should be located. Therefore, we would randomly choose a subset of the centroids of the tetrahedral elements that compose the solid domain as the IHR points. Given that we know in this special case of analytically known interior resonance frequencies that there are 53 of those in the range of interest, we pick 53 random points in the interior of the sphere. Figure 2 shows the less-than-ideal outcome where the interior modes cause ill conditioning manifested by the pronounced departures from the reference solution. We can identify a very minor blip at 1.5 kHz, another rather larger error around 2.15 kHz, and another at 2.75 kHz, and so on.
However, the poor performance of the 53 randomly chosen points was in fact due to the use of the symmetric formulation: the problem was also solved with a rectangular coefficient matrix using the QR decomposition (i.e., the least-square solution). Figure 3 demonstrates that the ill-conditioning errors vanished. Therefore, at least in this case, the formulation of the square system of equations is much less effective than the least-squares formulation. Moreover, there is other corroborating information in the open literature [23].
In Figure 3 we have evidence that the 53 random points were sufficient to ameliorate the ill-conditioning issue. Nevertheless, that does nothing to address the major deficiency of the original CHIEF: it is not known how many points will be needed and in which locations.    Next we direct our attention to the automatic point-selection procedure based on the computed interior modes. From the computed acoustic interior modes we can extract 450 interior points using the threshold (17). The results are presented in Figure 4 for the square system of equations, and in Figure 5 for the least-squares solution. There is not a big difference between these two solutions, but the least-squares solution seems to fit the reference curve a little bit better.
We can also demonstrate the use of a single interior point selected for each mode, using the zero threshold (18). For 53 acoustic modes we select 53 points at the locations of the maximum amplitude for each of the modes. With a square system, the loss of accuracy again appears at the locations of the interior eigenfrequencies, as shown in Figure 6. On the other hand, with the least-squares method a good solution is obtained (Figure 7).      Results computed with 53 interior points collected from 53 interior acoustic modes. Formulation (15) (i.e., a rectangular system with a QR solver) was used.

Discussion
The selection of the interior points can be made adaptively: when performing the sweep through the frequencies, only the points for modes in the vicinity of the current frequency need to be added for the IHR. This, however, may not lead to a huge payoff, since adding even hundreds of interior points does not affect the solution times greatly.
We should state clearly that we cannot provide a guarantee that the automatically selected points as described above will be sufficient for the elimination of all inaccuracies due to ill conditioning. However, we are confident that the selected points will be in general better than any random collection of the same size.
In any case, solving the interior acoustic mode problem using a finite element model of the interior cavity generates information on the counts and locations of the interior frequencies. That should be valuable information irrespective of the chosen CHIEF variant.

Conclusions
Computing the backscattering of harmonic acoustic waves from underwater elastic targets of an arbitrary shape is a challenging problem of considerable practical significance. The finite element method applied to the discretization of the target was here combined with the surface Helmholtz boundary element method used for the surrounding fluid medium.
As is well known, the boundary integral method yields non-unique solutions at certain wavenumbers. This failure is associated with the existence of eigensolutions of the Helmholtz equation in interior of the surface of interest (acoustic modes). The combined Helmholtz integral equation formulation (CHIEF) suggested by Schenk is commonly employed to combine the surface Helmholtz boundary integral with equations of the interior Helmholtz relation written down at selected points within the cavity of the scatterer (i.e., in the complement of the fluid domain).The difficulty associated with this approach is the lack of guidance on the necessary number of interior points and on their locations.
In the method proposed here, the acoustic modes in the complement of the fluid domain (i.e., in the interior of the scatterer) are computed using the finite element method, and the locations of the peaks of the acoustic pressure are identified. This information is used to guide the selection of the interior points at which to add the interior Helmholtz relation; the user can thus determine how many points should be employed and where they should be located.
Our numerical experiments demonstrated the utility of the proposed automatic selection of the CHIEF points. With the least-squares procedure, robust elimination of the loss of accuracy in the vicinity of the interior frequencies could be achieved. Both the random selection of CHIEF points and our algorithm were less successful with the square (Lagrange multiplier) formulation. Possible extensions of our algorithm may include adaptive choice of the CHIEF equations based on the knowledge of the spectrum of the interior problem: the number and locations of the CHIEF points may be chosen only from the interior frequencies close to the frequency of interest, while those interior frequencies are removed from the frequency of interest may be ignored. In any case, information about the interior modes obtained by solving the relatively inexpensive finite element model for the complement of the fluid domain is likely to be valuable when constructing the acoustic response curves.