Abstract
We propose a method of defining absorbing boundary conditions for use in finite element modeling of mechanoacoustic systems. The finite element model of an elastic body is expanded by a water domain, bounded by an axisymmetric surface, to which boundary conditions are applied. The boundary conditions are formed using the method of equivalent sources, so the integral relations for the pressure and its gradient are derived. The method retains its validity for the modeling of radiation in the low-frequency range, as it provides not only the absorption of outgoing waves but also the modeling of the added mass of vibrating elastic bodies.
1. Introduction
Numerical modeling of wave processes by means of the finite element method (FEM) inevitably requires dealing with the boundaries of the computational domain. The question is quite general and arises not only in underwater acoustics simulations but also, for example, when modeling electromagnetic waves. The so-called absorbing boundary conditions (BCs) should ensure the absorption of waves that come out of the computational domain and correct simulation of the near field. In the case of underwater acoustics, which will be discussed in this paper, that means that BC should be capable to simulate an added mass of liquid formed by the volume of water outside the computational domain. In other words, a search is made for such BC that would act similarly to the Sommerfeld condition but lead to a smaller error at finite distances from the sound source.
A literature review helps one to identify several competing approaches [1]. The first one to be mentioned is such a BC that is derived from the Helmholtz equation with certain approximations. Namely, they relate the pressure and its derivatives along the border to the normal derivative of the pressure. An example of BC, used in commercial software, is [2], and some overview of known methods is given in [3]. The applicability of such BC is limited either by the wave incidence angles or by other parameters. However, the mathematical expressions that define the BC turn out to be relatively simple, so they can be easily generalized to the time domain or applied in those problems which require more complicated research than a harmonic analysis [4].
The next approach is PML (perfectly matched layer), which consists of creating an “exotic” medium layer around the edges of the model [5]. This layer provides an increased absorption coefficient for propagating waves and a sufficient attenuation of the near field due to a complex transformation of coordinates, so the waves do not reach the new boundaries of the computational domain. In terms of the development of PML algorithms, topical issues are dealing with arbitrary shape boundaries [6], the accuracy of near-field modeling [7], and application in layered and dispersive media [8, 9].
Next follows the integral (i.e., nonlocal or global) BCs, that is, when one mathematical expression can relate the pressure values in far-away nodes. Acoustic infinite elements [10] are used in some cases. Their formulation differs from the above-mentioned BCs in several key respects. The outer region is subdivided into elements, each of which has a “base” on the surface of the finite element model. The weighted residual method is applied, similar to the conventional finite element method. The degrees of freedom corresponding to the interpolation functions in the infinite direction are added to the overall matrix of the system.
Mathematically, BCs for radiation in the finite element model are governed by expressing the normal derivative of the pressure at the boundary of the region through the values of pressure at the nodes of this boundary (except for PML, that is, a layer by itself). For boundaries of a regular shape, one can use the outgoing wave expansion to explicitly derive the so-called Dirichlet-to-Neumann (DtN) operator [11]. This operator transforms the vector of values at all border nodes into a vector of values at the same nodes.
For an arbitrary shape of the boundary, the relationship between and can be found by means of the boundary element method (BEM), that is, by numerically solving the Helmholtz integral equation [12, 13]. BEM is an independent and actively used [14] method for solving the problems of radiation and sound scattering, but the types of scattering bodies under consideration are limited to absolutely rigid and absolutely soft bodies and those bodies that can be described in terms of local impedance, for example, strongly absorbing bodies. As for radiation problems, BEM allows to calculate the field of panels with a given vibration velocity on the surfaces of such bodies. Modeling of elastic bodies is achieved by jointly solving the system of equations formed by the FEM and BEM [15], and in this case, the BEM mesh directly coincides with the surface mesh of the elastic body modeled by the FEM. This leads, firstly, to a long BEM computation time for the matrix coefficients connecting in pairs all nodes on the model surface, and secondly, it is expressed in further computational difficulties when operating with a large fully populated matrix. Recall that the FEM matrix is usually sparse, and the BEM matrix is dense. Some researchers try to optimize operations with matrices consisting of sparse and filled parts, starting from a proper rearrangement of columns and rows [11] and ending with iterative procedures based on the fast multipole method [14, 16].
A number of papers compare the performance of various methods for specifying absorbing BC [3, 17, 18]. In particular, according to [17], BC given by the BEM provides the best convergence in terms of the fastest rate of error reduction when refining the mesh. However, in general, which type of BC is preferable depends on the particular application. The result in performance may be different, depending on the implementation of the rest part of the software, such as the formation of a finite element model and its matrix, subsequent matrix operations, postprocessing, etc.
The algorithm for generating absorbing BC, proposed in this article, is based on the idea of the combined use of FEM- and BEM-like approaches. The problem of sound radiation and scattering by elastic bodies in water is studied. The difference from the known BEM implementations is the usage of the equivalent source method [13, 19–21] instead of calculating singular surface integrals. An original optimization of the BEM matrix with its presentation in the Toeplitz form is also proposed. Finite element models suitable for applying the algorithm are prepared in a special way. The region of the external medium, modeled by the BEM, is not directly connected to the complex contact surface of the liquid and the body. The object under study is expanded with a water domain, i.e., finite elements that simulate a certain volume of the water medium, so that the shape of the resulting outer boundary is a rotating body of revolution, e.g., a cylinder with end caps. The axisymmetric geometry of the outer region contributes to a significantly faster formation of the BC matrix in comparison with the general BEM. This extra water domain allows the evanescent waves to decay, so the mesh on the outer boundary is coarse, and its step is determined by the length of the acoustic wave, not by the flexural waves in the body. However, at the same time, the extra water domain is not too large in order not to greatly increase the dimension of the system. Elongated bodies are usually investigated in underwater acoustics, and it is convenient to expand them up to a cylinder (compare, for example, if the task is to inscribe an elongated body into a ball; this will lead to the addition of a large number of extra degrees of freedom to the model).
2. Methods
Finite element modeling of the studied harmonic vibroacoustic processes is performed in the space of two variables, namely, the displacement field , which oscillates with a cyclic frequency ω, and the displacement potential θ, which describes the fluid medium dynamics. One should note that usually, acoustic waves in water are described by a pressure field , while finite element modeling is carried out. Here, we exploit the displacement potential that is related to pressure as in the case of harmonic process. This change of variables helps to make some of the further expressions symmetric. Thus, the studied problem of finding a response to give force can be described in a discrete form by a system of linear equations with symmetric matrix elements, which is similar to the one presented by the authors in [4, 22]: where and are the matrices of stiffness and mass of the deformable solid body, and are the matrices describing the transfer and compression of a fluid, and MFSI and Z are the matrices determining the boundary conditions of no fluid loss at the fluid-structure contact surface and nonreflection of acoustic waves from the boundaries of the computational domain, respectively. Algorithms for the formation of matrix components , , , and can be found in the classical papers [23, 24]. A procedure for efficient formation of the matrix MFSI, which establishes the relationship between and θ, can be found in [22].
The matrix is frequency-dependent and has nonzero elements only for nodes at the outer boundary of the computational domain. It defines the connection between the value of the θ potential and its derivative θn with respect to the normal to the boundary of the computational domain. This paper is just dedicated to the method of calculation of this matrix. Formally, this matrix can be obtained by differentiating the following function with respect to the values of the generalized degrees of freedom: where ρ is the fluid density, ω is the angular frequency, is the fluid boundary, Θ is the displacement potential (as a continuous distribution), and is the outer normal to . The quantity δW determines the energy flux increment through at a given increment of the values of the generalized degrees of freedom.
A link between and Θ is an essential point in derivation of BC. Assume that the whole domain inside is filled with water medium of the same kind as outside. It is known [13, 19–21] that an arbitrary distribution of pressure (or the potential) on the surface can be very closely approximated by placing in some way the sources of variable mass (i.e., monopole acoustic sources) inside (see Figure 1). As soon as the equivalent sources, creating the required field on the boundary, are determined, the value of the normal derivative will be determined in a straightforward way.
At this point, we must require that the mesh on the boundary is structured (regular). It must coincide with itself when rotated by some angle Δφ. So it should have elements along the circle (see figures below and mind that the notation M/2 is used due to so-called finite elements of the second order; i.e., there is one intermediate point on each edge). The nodes of these elements are going to be organized in rings, and there will be two poles, one node per each pole. To ensure that the matrix is conditioned, the number of degrees of freedom, provided by equivalent sources, should be the same as the number of degrees in the nodes of boundary . A group of equivalent sources is a scaled copy of the boundary mesh nodes. According to this, equivalent sources are defined as continuous coaxial rings (Figure 1) with numbers in a short distance away from the surface (proper offset was found using additional numerical experiments). Herein, it is important to specify options for the distribution of along the angular coordinate φ for each location option. That is, each ring as a sound source is modulated with a series of functions with , and each term of the series is treated as a unique degree of freedom. So the volumetric velocity is expressed in the form , where is the amplitude of the mth harmonic of the kth ring source.
Matrix is substituted in (1) so that it is multiplied by unknown variables, which are the values of the potential in the ith nodes of the FEM. Before we derive the matrix in that form, a supplementary problem is going to be considered. Let us derive expression for δW as if were the unknown variables in the problem. Later, we will derive a transform from the values θi in FEM nodes to the amplitudes of sources . Consider the fundamental solution of the Helmholtz equation in the form where is the longitudinal coordinate of the ring source location. Denote a distance as
Expression (2) can be reduced to the matrix form where the elements of the matrix are determined by the expression
Recall that the submerged structure that is supposed to be studied can have an arbiter shape, and it is going to be expanded by water elements to make an axisymmetric outer boundary that we derive the BC for. Due to these facts, the integration with respect to in Equation (2) can be reduced to integration with respect to the generatrix of the axisymmetric FEM boundary. Then, the elements of the matrix with allowance for the azimuthal angle orthogonality can be determined as follows: where is equal to
Note that both integrals contain symmetric functions of at . So, one may compute only cases.
In theory, according to Green’s second formula, the matrix itself should be complex symmetric (not Hermitian) in terms of swap and indices. It is essential to obtain this property in the result of the numerical calculation routine, carried out according to (8). For this purpose, the Gauss-Lobatto quadrature is applied and the integration is performed with an adaptive mesh step. Integral (8) is replaced with a finite sum over line elements that coincide with the edges of the FEM elements in the first iteration. Then, each line element is reduced to a small size until the difference in the values of the integral over the surface element before and after the next reduction exceeds the specified permissible error.
Expression (5) relates the flux increment of acoustic energy transmitted through the surface from equivalent acoustic sources. To use this expression in FEM calculations, it is necessary to pass from the generalized degrees of freedom (the amplitudes of ring sources) to the degrees of freedom θ (the potential values at the nodes of the FE mesh at the boundary ). In this case, it is expedient to perform this transformation not directly, but in two stages, in the first of which the generalized degrees of freedom can be expressed in terms of the values of the Fourier image of the potential θ at the generatrix points of the surface . Namely, if we apply the Fourier transform in the azimuthal coordinate to expression (3) for fixed coordinates of the generatrix points and , we can obtain in matrix form where the elements of the matrix are equal to
The key point is the possibility for us to compute the inverse matrix . The authors did not perform a strict mathematical proof of the convergence of this method in the problem being solved. Instead, wide experience in the practical use of the equivalent source method was taken into account, and the accuracy in model problems with a known solution was explored. Experience has shown that there are the following requirements for using the equivalent source technique in this geometry of the problem: the displacement of ring sources deep from the surface should be from a quarter to a half of the step of the FE mesh on the surface .
Thus, in view of (9), the above functional reduces to the following form:
Passing from the Fourier transform to the original, the potential values in their final form can, instead of (11), be written as where the elements of the matrix of coupling between the mth harmonic of the Fourier image and the original are determined by the relationship where is the degree of freedom corresponding to the ring cross section of the surface ; φs+ and φs- are the angular coordinates of the integration sector, which are determined by the FEM mesh; and Ψ is the shape function representing the interpolation of the pressure field between neighboring nodes of the ring cross section in the angle sector . In this case, according to whether the ring cross section lies in the plane of the FE angular nodes or passes in the plane of only the middle nodes (as in Figure 2), the shape functions combine the potential of either three or eight nodes, respectively. Namely, in the first case (for the ring cross section passing through nodes 1, 4, and 8 in Figure 2), and for the second case (for the ring cross section passing through nodes 5 and 7 in Figure 2), where is a parametric coordinate that within the limits of each sector of angles from varies from for to for .
If we replace the variables in Equation (13), where then Equation (13) will be reduced to the form
The final analysis matrix of the environment can be obtained by differentiating expression (12) by the vector of increments of the degree of freedom {δθ} and will represent a multiplier before the vector {θ}.
It is crucial that the matrix (19) is determined by a complex frequency dependence in accordance with expressions (8) and (10), which precludes identifying it as a matrix of mass loss or stiffness of a mechanical system. Therefore, this matrix is further called the frequency-dependent impedance matrix of the acoustic environment or the impedance matrix for short.
Once again, the impedance matrix by its origin relates the potential (pressure) on the boundary with the weighted normal derivative of the potential (of the pressure) so that we are able to compute (2), knowing only Θ:
One may note that since we defined that converts the field into the equivalent sources’ amplitudes and we are able to compute the gradient of the field, produced by the equivalent sources, it is possible to derive in a simpler way [25]. However, the procedure, described in this paper, is more preferable due to the fact that the matrix numericaly computed only in this way turns out to be very close to the symmetric matrix. It was declared above that the symmetry of is a property of physical processes that lie behind the mathematics, but an inefficient algorithm may contain such series of floating point operations that turns the result unsymmetrical.
The next remarkable property of matrix is its Toeplitz-like structure that can be explained in the following way. Consider that the whole model is rotated by around the cylinder axis. Boundary nodes are going to take places of other nodes; i.e., in the case of proper numbering, a node goes to the position of a node j+Ngroup, if it is not a too large number, and otherwise, the node number is shifted in a circular way. That is true except for two nodes, which are located right on the rotational axis. An element of , denoted as , defines a reaction of external homogeneous medium that appears in the jth node, when a load is applied in the kth node. The reaction must be the same in cases like . This leads to the following structure of the whole matrix: where , , … are the matrix blocks of size; is a matrix block of size (the number of rows goes first); and is size.
Authors’ software implementation of the described method functions so that it saves the computed matrices on a hard drive, before it is used in the FEM subprogram. This turns out to be a practical approach when the FEM model is adjusted, but its border remains the same. The above-noted Toeplitz-like structure of allows to store only the set of unique blocks: , , , …, and . This speeds up input-output operations and saves the disk space.
3. Results and Discussion
3.1. Numerical Validation
As a verification of the above method for calculating the impedance matrix , we assume that there is a certain test acoustic source inside . In Figure 3, we show the result of the test calculation of the projection of the pressure gradient of the test source field on the outer normal to the cylindrical surface ). Markers (dots and crosses) are obtained by the given numerical algorithm. That is, first, we compute the pressure filed, produced by the test source at the positions of nodes laying on the surface . And second, we treat the set of these pressure values as a vector and multiply them by the matrix . Lines are obtained using a reference analytical equation; i.e., we do the straightforward computation of the normal projection of the pressure gradient in those points. The markers and lines coincide that prove the high precision of the discrete computation routine.
(a)
(b)
(c)
(d)
Two lines on each plot show the real and imaginary parts. The algorithm has been tested using the following model: mesh step 0.01 m, cylinder diameter 1.6 m, and length 14 m; the test source was a monopole at the cylinder axis (Figures 3(a) and 3(b)) and a quadrupole source (Figures 3(c) and 3(d)). The frequency was 500 Hz (Figures 3(a) and 3(c)) and 1500 Hz (Figures 3(b) and 3(d)).
The development of this type of boundary conditions is a part of the general line for the development of the super-element technique for harmonic problems of noise emission [25], which is a kind of a domain decomposition technique. Within the framework of such a statement, it is possible to represent the computational domain containing the emitting body as a set of conjugate finite element models, so-called super-elements, and the environment outside the computational domain as a separate “boundary” super-element, whose matrix is calculated according to the above algorithm.
The validity of this approach can be demonstrated by the example of a simple system consisting of a super-element of a cylindrical volume of the medium, which is simulated by finite elements (Figure 4(a)) and a “boundary” super-element. The cylinder sizes correspond to the above boundary and are as follows: the diameter is 1.6 m and the length is 14 m. Figures 4(b) and 4(c) show the results of calculating the field of a point source placed in a cylindrical volume.
(a)
(b)
(c)
3.2. Postprocessing
Additional calculations are used to determine the field characteristics at any other points of the environment. To find the complex amplitude of the potential field and the respective sound pressure at an arbitrary point outside , the Helmholtz-Kirchhoff integral is employed in the calculations: where is Green’s function of free space. The quantity is the result of calculating the near field in terms of FEM, while its normal derivative is computed using again. After the comparison of (2) and (20), one can follow the same rule and write out
An infinitely distant point should be considered as to find the radiation pattern. Its radius vector can conveniently be expressed via a unit directing vector and the module . Under the conditions and , Green’s function can be replaced by an approximate expression:
Figure 5 shows an example of calculation of the radiation pattern. The model again consists of a cylindrical volume of water enclosed within a boundary that simulates the presence of the same acoustic medium in the outer region (Figure 5(a)). Test sources, namely, a dipole or a quadrupole, are placed inside this volume. First of all, the pressure field inside the cylindrical volume was calculated (Figures 5(b) and 5(c)). Then, the radiation pattern in the far-field region was calculated (Figures 5(d) and 5(e)). The calculation result complies with the theoretical model of the chosen source.
(a)
(b)
(c)
(d)
(e)
3.3. Estimation of Computation Time
Finally, we provide the results of computation time measurement according to the mesh size. For testing, we use models consisting of a cylinder, made of a “finite element water” and a “boundary” super-element, which in fact is a set of BC coefficients, governed by the impedance matrix computed by the described method. Recall that the boundary mesh contains approximately KM nodes, i.e., K nodes in line and M nodes along the circle, and “approximately” has been written due to second-order elements. Corresponding finite element models were meshed in a similar way. First, we measured , which is a time to compute (see Table 1). A set of the model has been chosen to do that, and the one used for another test above and depicted on Figure 4(a) is no. 4 in the present list.
Tests have shown that the most time-consuming operation is computation of coefficients according to expression (8). And the most complicated part of it is integration in the vicinity of the sources, i.e., when or are small enough. Integration step is automatically refined in those places. So the complexity of computation of the total integral in (8) is , and it appeared that in the tested models, . The number of coefficients is , so the total complexity is . The time values, given by Table 1, approximately match this function.
However, in the case of larger models, one may expect that the computation of will take . The other computation-intensive operations, while computing , are the matrix inversion for each and matrix multiplications (see equation (19)). So we have to invert a dense matrix and multiply several dense matrices of that size same amount of times. Each inversion and each multiplication have a complexity of , so total complexity of (19) is estimated as .
Usually, a surface mesh is scaled with the decrease of an acoustic wavelength as . So theoretical computation complexity of the proposed algorithm for matrix might be between and . Hopefully, while dealing with models, that could be studied on an “advanced table computer”; without large computational servers, computation complexity is found.
Next, we studied the solution of the whole problem that could be done in two ways. Both of them start with the following: (i)Compute submatrices for each finite element for time (ii)Compute for time (see Table 1)
Then, the first way continues with the following: (i)Assemble a total matrix, consisting of a sparse part originated by the finite element model and a dense part, originated by BC, namely, matrix (negligible time for assembly)(ii)Factorize that sparse matrix with a sufficiently large dense part and solve the system, the sparse solver is used, and time is
The second way is as follows: (i)Compute a Schur complement for the sparse finite element matrix to obtain a dense matrix of a smaller size, each column and row of this dense matrix corresponds to the nodes on the external boundary, the sparse solver is used, and time is (ii)Assemble a dense total matrix, consisting of the Schur complement from finite elements and (negligible time for assembly)(iii)Factorize that dense matrix and solve the system, with time
The time for listed operations is given in Table 2. The model numbers correspond to parameters, given in Table 1. The Schur complement as an intermediate step helps to reduce the duration of the last step; however, the total computational time is approximately the same for both ways, including or excluding the Schur step.
All tests were carried out with the use of the authors’ software on a computer with a 4-core 3.4 GHz CPU (might be automatically boosted up to 3.9 GHz) and 32 GB of RAM. All the computation times are valid for an analysis at a single frequency. For analysis in a frequency band, the values of time should be multiplied by the number of frequencies.
4. Concluding Remarks
We have presented a method to define the boundary conditions of sound radiation in the problems of numerical modeling of mechanoacoustic systems. The method does not require significant computational resources and simulates the acoustic processes with high accuracy for inhomogeneous elastic elongated bodies in a wide frequency range. The method of super-element formulation is successfully used in the acoustic design of marine facilities.
Data Availability
The data used to support the findings of this study is the numerical model geometry. We believe that it can be easily reproduced basing on the description in the text. If one faces problem with that, the numerical model is available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
The authors are grateful to Eugeny M. Sokov and the members of Software development sector of the Institute of Applied Physics RAS for their contribution in the realization of the author’s methods. This work was fulfilled under the State Contracts with the Ministry of Education and Science of the Russian Federation [0030-2021-0017 and 0030-2022-0003].