# Integral Absorbing Boundary Conditions Optimized for Modelling of Acoustic Radiation of Elongated Bodies

Authors: Mikhail B. Salin; Sergey A. Smirnov; Anatoliy S. Suvorov; Irina A. Usacheva; Irina A. V’yushkina
Journal: Journal of Applied Mathematics
 (2022)

---

## 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.


---

## Body

## 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∂p/∂n at the boundary of the region through the values of pressure p 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 p values at all border nodes into a vector of ∂p/∂n values at the same nodes.For an arbitrary shape of the boundary, the relationship between∂p/∂n and p 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 fieldu, 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 p, while finite element modeling is carried out. Here, we exploit the displacement potential that is related to pressure as θ=p/ρω2 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 b 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]:
(1)KS000−ω2MSMFSIMFSIT−MF−ω4000GF−i⋅000Zuθ=b0,where Ks and Ms are the matrices of stiffness and mass of the deformable solid body, MF and GF 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 Ks, Ms, MF, and GF can be found in the classical papers [23, 24]. A procedure for efficient formation of the matrix MFSI, which establishes the relationship between u and θ, can be found in [22].The matrixZ 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:
(2)δW=ρω2∫S∂Θ∂nδΘdS,where ρ is the fluid density, ω is the angular frequency, S is the fluid boundary, Θ is the displacement potential (as a continuous distribution), and n is the outer normal to S. The quantity δW determines the energy flux increment through S at a given increment of the values of the generalized degrees of freedom.A link between∂Θ/∂n and Θ is an essential point in derivation of BC. Assume that the whole domain inside S 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 S can be very closely approximated by placing in some way the sources of variable mass q (i.e., monopole acoustic sources) inside S (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.Figure 1
The geometry of the fluid boundary and equivalent sources.At this point, we must require that the mesh on the boundaryS is structured (regular). It must coincide with itself when rotated by some angle Δφ. So it should have M/2=2π/Δφ 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 K 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 S. 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 k=1..K in a short distance away from the surface S (proper offset was found using additional numerical experiments). Herein, it is important to specify M options for the distribution of q along the angular coordinate φ for each k location option. That is, each ring as a sound source is modulated with a series of eimφ functions with m=−M/2..M/2−1, and each term of the series is treated as a unique degree of freedom. So the volumetric velocity is expressed in the form qkm=akmeimφ, where akm is the amplitude of the mth harmonic of the kth ring source.MatrixZ is substituted in (1) so that it is multiplied by unknown variables, which are the values of the θi potential in the ith nodes of the FEM. Before we derive the matrix Z in that form, a supplementary problem is going to be considered. Let us derive expression for δW as if akm 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 akm. Consider the fundamental solution of the Helmholtz equation in the form
(3)Θ=14π∑k,makmeimφrk′∫02πe−imφ′eiω/cr2+rk′2−2rrk′cosφ′+z−zk′2r2+rk′2−2rrk′cosφ′+z−zk′2dφ′,where zk′ is the longitudinal coordinate of the ring source location. Denote a distance Rk as
(4)Rk=r2+rk′2−2rrk′cosφ′+z−zk′2.Expression (2) can be reduced to the matrix form
(5)δW=ρω2δaplTGakm,where the elements of the matrix G are determined by the expression
(6)Gk,pm,l=116π2∫Srk′rp′eim+lφ∫02πe−imφ′eiω/cRkRkdφ′⋅∂∂n∫02πe−ilφ′eiω/cRpRpdφ′dS.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 boundaryS that we derive the BC for. Due to these facts, the integration with respect to S in Equation (2) can be reduced to integration with respect to the generatrix L of the axisymmetric FEM boundary. Then, the elements of the matrix G with allowance for the azimuthal angle orthogonality can be determined as follows:
(7)Gk,pm,l=Gk,pm,0,m≠−l,where Gk,pm is equal to
(8)Gk,pm=18π∫Lr⋅rk′⋅rp′⋅∫02πeiω/cRk−mφ′Rkdφ′⋅∫02πiωc−1Rpeiω/cRp+mφ′Rp∂Rp∂ndφ′dL.Note that both integrals contain symmetric functions ofφ at e−imφ′. So, one may compute only m≥0 cases.In theory, according to Green’s second formula, the matrixGmk,p itself should be complex symmetric (not Hermitian) in terms of swap k and p 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 S from equivalent acoustic sources. To use this expression in FEM calculations, it is necessary to pass from the generalized degrees of freedom a (the amplitudes of ring sources) to the degrees of freedom θ (the potential values at the nodes of the FE mesh at the boundary S). 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 S. Namely, if we apply the Fourier transform in the azimuthal coordinate φ to expression (3) for fixed coordinates of the generatrix points rp and zp, we can obtain in matrix form
(9)θpm=Tmakm,where the elements of the matrix Tm are equal to
(10)Tp,km=rk′2∫02πe−imφ′eiω/crp2+r′k2−2rprk′cosφ′+zp−zk′2rp2+r′k2−2rprk′cosφ′+zp−zk′2dφ′.The key point is the possibility for us to compute the inverse matrixTm−1. 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 S should be from a quarter to a half of the step of the FE mesh on the surface S.Thus, in view of (9), the above functional reduces to the following form:
(11)δW=ρω2∑mδθ~mTTm−TGmTm−1θ~m.Passing from the Fourier transform to the original, the potential values in their final form can, instead of (11), be written as
(12)δW=ρω2δθT∑mNmTTm−TGmTm−1Nmθ,where the elements of the matrix of coupling between the mth harmonic of the Fourier image and the original Nm are determined by the relationship
(13)Nj,pm=∑s∫φs−φs+Ψj,p,se−imφdφ,where j is the degree of freedom corresponding to the ring cross section p of the surface S; φs+ and φs- are the angular coordinates s 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 p in the angle sector s. 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),
(14)Ψ1=−η21−η,Ψ4=η21+η,Ψ8,s=1−η2,Ψ2=Ψ3=Ψ5=Ψ6=Ψ7=0,and for the second case (for the ring cross section passing through nodes 5 and 7 in Figure 2),
(15)Ψ1=Ψ2=Ψ3=Ψ4=−141−η2,Ψ6=Ψ8=121−η2,Ψ5=121−η,Ψ7=121+η,where η is a parametric coordinate that within the limits of each sector of angles from φs−доφs+ varies from η=−1 for φs− to η=1 for φs+.Figure 2
Layout of the boundary condition formation model.If we replace the variables in Equation (13),
(16)φ=Δφ2η+φs0,where
(17)φs0=φs++φs−2,Δφ=φs+−φs−, then Equation (13) will be reduced to the form
(18)Nj,pm=Δφ2∑se−imφs0∫−11Ψj,p,se−imΔφη/2dη.The final analysis matrix of the environmentZ 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 {θ}.
(19)Z=ρω2∑mNmTTm−TGmTm−1Nm.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 boundaryS with the weighted normal derivative of the potential (of the pressure) so that we are able to compute (2), knowing only Θ:
(20)δW=δθTZθ.One may note that since we definedTm−1 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 Z in a simpler way [25]. However, the procedure, described in this paper, is more preferable due to the fact that the matrix Z numericaly computed only in this way turns out to be very close to the symmetric matrix. It was declared above that the symmetry of Z 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 matrixZ 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 j 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 Z, denoted as Zjk, 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 Zjk=Zj+Ngroup,k+Ngroup. This leads to the following structure of the whole matrix:
(21)Z=DAAA⋯AATB1B2B3⋯B2TATB2TB1B2⋯B3TATB3TB2B1B4T⋯⋯⋯,where B1, B2, … are the matrix blocks of Ngroup×Ngroup size; A is a matrix block of 2×Ngroup size (the number of rows goes first); and D is 2×2 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 ofZ allows to store only the set of unique blocks: A, B1, B2, …, and D. 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 matrixZ, we assume that there is a certain test acoustic source inside S. 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 (∂P/∂n). 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 S. And second, we treat the set of these pressure values as a vector and multiply them by the matrix Z. 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.Figure 3
The result of testing the impedance matrix in the model problem. The diagrams show the projection of the pressure gradient of the test source field on the outer normal to the cylindrical surface: real and imaginary parts (see the text for explanations).
(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 (Figures3(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 Z 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 (Figure4(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.Figure 4
The result of calculations based on a model consisting of the SE of the cylindrical volume of the medium, the modeled FE (a), and the boundary super-element. The point source is located in the medium. The pressure modulus in dB at the points on the cylinder axis (b) and the instantaneous pressure value in the YZ plane (c) are shown. The frequency is 500 Hz and the mesh step is 0.1 m; cylindrical volume dimensions: diameter 1.6 m and length 14 m.
(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 pointr′ outside S, the Helmholtz-Kirchhoff integral is employed in the calculations:
(22)θr′=∫Sθr∂g∂nr′,r+∂θr∂ngr′,rdS,where
(23)gr′,r=14πe−iω/cr′−rr′−r,is Green’s function of free space. The quantity θr is the result of calculating the near field in terms of FEM, while its normal derivative is computed using Z again. After the comparison of (2) and (20), one can follow the same rule and write out
(24)∫S∂θr∂ngr′,r≈Zgr′,rj.An infinitely distant point should be considered asr′ to find the radiation pattern. Its radius vector can conveniently be expressed via a unit directing vector e′ and the module r′:r′=r′e′. Under the conditions ωr′>>c and r′>>r, Green’s function can be replaced by an approximate expression:
(25)gr′,r=14πe−iω/cre′r′.Figure5 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.Figure 5
Example of a test calculation. The frequency is 1000 Hz, the cylinder length is 8 m, and the diameter is 1.6 m. (a) The calculation model geometry: the green points are the nodes of the outer boundaryS, and the red and blue points are the nodes where the sources can be specified; (b) view of the pressure field of a dipole source in the cross section, in the near field: the real part is constructed, and the white rectangle covers the source area; (c) view of the pressure field of a quadrupole source in the cross section, in the near field; (d) radiation pattern of a dipole source, ∣P∣; and (e) radiation pattern of a quadrupole source, ∣P∣.
(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 matrixZ 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 tZ, which is a time to compute Z (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.Table 1
Estimation of computation time for impedance matrixZ, which defines BC coefficients in the model.
Model no.Nodes along,KNodes around,MComputation time,tZ (s)13416327036143804823414072945200120370Tests have shown that the most time-consuming operation is computation ofGk,pm coefficients according to expression (8). And the most complicated part of it is integration in the vicinity of the sources, i.e., when Rk or Rp are small enough. Integration step is automatically refined in those places. So the complexity of computation of the total integral in (8) is OKM+const, and it appeared that in the tested models, KM≪const. The number of Gk,pm coefficients is K2M, so the total complexity is OK2M. The time values, given by Table 1, approximately match this function.However, in the case of larger models, one may expect that the computation ofGk,pm will take OK3M2. The other computation-intensive operations, while computing Z, are the matrix inversion for each Tm−1 and matrix multiplications Tm−TGmTm−1 (see equation (19)). So we have to invert a dense K×K matrix M and multiply several dense matrices of that size same amount of times. Each inversion and each multiplication have a complexity of OK3, so total complexity of (19) is estimated as OK3M.Usually, a surface mesh is scaled with the decrease of an acoustic wavelength asKM~λ−2. So theoretical computation complexity of the proposed algorithm for matrix Z might be between Oλ−3 and Oλ−5. Hopefully, while dealing with models, that could be studied on an “advanced table computer”; without large computational servers, computation complexity Oλ−3 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 timetl(ii)
ComputeZ for time tZ (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, matrixZ (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 istiThe 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 ists(ii)
Assemble a dense total matrix, consisting of the Schur complement from finite elements andZ (negligible time for assembly)(iii)
Factorize that dense matrix and solve the system, with timetgThe time for listed operations is given in Table2. 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.Table 2
Estimation of computation time for the whole problem solution.
Model no.tl (s)ti (s)Total time for way 1tl+tZ+ti (s)ts (s)tg (s)Total time for way 2tl+tZ+ts+tg (s)1126116224201.53.62132.37324.57.13745441439.63514455.629467058.7222656All 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.


## 3.1. Numerical Validation

As a verification of the above method for calculating the impedance matrixZ, we assume that there is a certain test acoustic source inside S. 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 (∂P/∂n). 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 S. And second, we treat the set of these pressure values as a vector and multiply them by the matrix Z. 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.Figure 3
The result of testing the impedance matrix in the model problem. The diagrams show the projection of the pressure gradient of the test source field on the outer normal to the cylindrical surface: real and imaginary parts (see the text for explanations).
(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 (Figures3(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 Z 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 (Figure4(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.Figure 4
The result of calculations based on a model consisting of the SE of the cylindrical volume of the medium, the modeled FE (a), and the boundary super-element. The point source is located in the medium. The pressure modulus in dB at the points on the cylinder axis (b) and the instantaneous pressure value in the YZ plane (c) are shown. The frequency is 500 Hz and the mesh step is 0.1 m; cylindrical volume dimensions: diameter 1.6 m and length 14 m.
(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 pointr′ outside S, the Helmholtz-Kirchhoff integral is employed in the calculations:
(22)θr′=∫Sθr∂g∂nr′,r+∂θr∂ngr′,rdS,where
(23)gr′,r=14πe−iω/cr′−rr′−r,is Green’s function of free space. The quantity θr is the result of calculating the near field in terms of FEM, while its normal derivative is computed using Z again. After the comparison of (2) and (20), one can follow the same rule and write out
(24)∫S∂θr∂ngr′,r≈Zgr′,rj.An infinitely distant point should be considered asr′ to find the radiation pattern. Its radius vector can conveniently be expressed via a unit directing vector e′ and the module r′:r′=r′e′. Under the conditions ωr′>>c and r′>>r, Green’s function can be replaced by an approximate expression:
(25)gr′,r=14πe−iω/cre′r′.Figure5 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.Figure 5
Example of a test calculation. The frequency is 1000 Hz, the cylinder length is 8 m, and the diameter is 1.6 m. (a) The calculation model geometry: the green points are the nodes of the outer boundaryS, and the red and blue points are the nodes where the sources can be specified; (b) view of the pressure field of a dipole source in the cross section, in the near field: the real part is constructed, and the white rectangle covers the source area; (c) view of the pressure field of a quadrupole source in the cross section, in the near field; (d) radiation pattern of a dipole source, ∣P∣; and (e) radiation pattern of a quadrupole source, ∣P∣.
(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 matrixZ 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 tZ, which is a time to compute Z (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.Table 1
Estimation of computation time for impedance matrixZ, which defines BC coefficients in the model.
Model no.Nodes along,KNodes around,MComputation time,tZ (s)13416327036143804823414072945200120370Tests have shown that the most time-consuming operation is computation ofGk,pm coefficients according to expression (8). And the most complicated part of it is integration in the vicinity of the sources, i.e., when Rk or Rp are small enough. Integration step is automatically refined in those places. So the complexity of computation of the total integral in (8) is OKM+const, and it appeared that in the tested models, KM≪const. The number of Gk,pm coefficients is K2M, so the total complexity is OK2M. The time values, given by Table 1, approximately match this function.However, in the case of larger models, one may expect that the computation ofGk,pm will take OK3M2. The other computation-intensive operations, while computing Z, are the matrix inversion for each Tm−1 and matrix multiplications Tm−TGmTm−1 (see equation (19)). So we have to invert a dense K×K matrix M and multiply several dense matrices of that size same amount of times. Each inversion and each multiplication have a complexity of OK3, so total complexity of (19) is estimated as OK3M.Usually, a surface mesh is scaled with the decrease of an acoustic wavelength asKM~λ−2. So theoretical computation complexity of the proposed algorithm for matrix Z might be between Oλ−3 and Oλ−5. Hopefully, while dealing with models, that could be studied on an “advanced table computer”; without large computational servers, computation complexity Oλ−3 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 timetl(ii)
ComputeZ for time tZ (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, matrixZ (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 istiThe 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 ists(ii)
Assemble a dense total matrix, consisting of the Schur complement from finite elements andZ (negligible time for assembly)(iii)
Factorize that dense matrix and solve the system, with timetgThe time for listed operations is given in Table2. 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.Table 2
Estimation of computation time for the whole problem solution.
Model no.tl (s)ti (s)Total time for way 1tl+tZ+ti (s)ts (s)tg (s)Total time for way 2tl+tZ+ts+tg (s)1126116224201.53.62132.37324.57.13745441439.63514455.629467058.7222656All 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.



---

*Source: 9524376-2022-02-22.xml*