A Stability Analysis of the Direct Interpolation Boundary Element Method applied to acoustic wave propagation problems using the Modal Superposition Technique

Abstract

The formulations of the boundary element method for dynamic problems, whether acoustic or elastodynamic, that use a more straightforward fundamental solution and consequently require the use of Radial basis functions to transform domain integrals into boundary integrals present instability when reduced time increments are used. This problem does not arise in numerical techniques like finite element and finite difference methods. Despite the higher quality of results presented by the Direct Interpolation Method (DIBEM), a new radial basis formulation, instability problems persist for very small time steps, demanding mesh refinement. Several factors were examined, such as the robustness of the inertia DIBEM matrix concerning conditioning and positivity, densification of the number of internal interpolation points inside, the type of radial basis function, the composition of the time marching scheme, and others. This work evaluates if the cause of this stability is linked to the modal participation of non-real frequencies, which arise in the high spectrum due to numerical inaccuracy when dealing with non-symmetric matrices.

Keywords:
Direct Interpolation Boundary Element Method; Modal Superposition; Radial Basis Functions

Graphical Abstract

1 INTRODUCTION

Dynamics problems are among the most challenging for engineering and also for the powerful numerical methods currently used to solve them. Concerning BEM, its advantageous characteristics in reducing the dimension of the problem and other interesting features (Brebbia et al., 1984) are insufficient to meet specific demanding industrial sectors, which demand high computational performance solving large-size systems since the BEM matrices are full and asymmetric. Thus, many boundary element researches are directed to techniques for reducing computational time improving strategies such as the subdomain or sub-regions technique (Ramšak and Škerget, 2007; Chai et al., 2022). This technique becomes attractive since it considerably decreases the necessary amount of memory and makes the application of efficient solvers suitable. Economy in time processing impulses iterative solution schemes, such as Adaptive Cross Approximation (ACA) (Bebendorf, 2000; Bebendorf and Rjasanow, 2003; Bebendorf and Grzhibovskis, 2006), Fast Multipoles (Rokhlin, 1985; Greengard and Rokhlin, 1987; Liu, 2006), and others. Also significant in this field is the development of several strategies concerning out-of-core programming (Shiara and Paschoalini, 2023). However, developing more accurate and reliable models is essential, particularly for solving dynamic problems. In this sense, new models have been continually presented in the literature, regarding the effective transformation of domain integrals into boundary integrals with greater accuracy, such as the Dual Reciprocity Method (DRBEM) (Nardini and Brebbia, 1983a), General Radial Basis Approximation (Wen et al., 1998), Radial integration Method (RIM) (Gao,2002), Boundary Element–Rayleigh Integral Method (Kirkup et al., 2013), Meshless-based Integration Method (RBIM) (Narvaez and Useche, 2022), and Direct Interpolation Method (DIBEM) (Loeffler et al., 2015a).

It is well-known that the classical formulation of the Boundary Element Method (BEM) employs fundamental solutions strictly related to the partial differential equation one wants to solve. Such solutions are known as Green's functions and, particularly in the case of acoustic problems, they consist of solving the scalar wave equation considering an infinite medium subject to concentrated and impact loads (Mansur and Brebbia, 1985). However, this formulation is complex despite its elegance, excellent results, and mathematical consistency. Beyond this, it demands large computational effort beyond being restricted to homogeneous media unless one uses the mentioned subdomains. Thus, aiming to provide further simplicity and versatility to BEM, some formulations were developed in which more straightforward fundamental solutions are used, which implies additional approximations related to transforming domain integrals into boundary integrals.

Radial basis functions (Buhmann, 2003) are one of the most effective procedures for reaching this objective. In the case of elastodynamics, it was applied especially initially in free vibration problems and later in issues of response over time (Nardini and Brebbia, 1983b). Later, its extension to acoustic scalar problems and heat transmission problems was made (Loeffler and Mansur, 1987). However, particularly in dynamics, the DRBEM has presented extreme restrictions regarding the stability of the time advance scheme, based on the idea of discretization contained in the finite difference method, such as the Houbolt, Newmark, and Wilson-theta scheme (Castillo and Loeffler, 2003), used successfully with the Finite Element Method. For problems in which the response is dominated by the system's low and intermediate frequency components, all these incremental schemes produce satisfactory solutions. However, such integration schemes used by domain methods always induced uncertainties for impact problems in which more elevated modal content exists, motivated mainly due to the BEM discretization being limited exclusively to the boundary (Balista et al., 2023).

It is known from simulations with DRBEM that there is numerical instability when using minimal time intervals (Loeffler and Mansur, 1987; Castillo and Loeffler, 2003). In domain discrete methods, instability problems occur for huge time increments, and, therefore, unconditionally stable integration schemes are adopted. However, precisely what causes this instability with the BEM is still not known. One hypothesis is the inaccuracy in the constitution of the mass matrix, which is constructed by approximations with radial basis functions; another hypothesis is related to the stability criteria of incremental schemes studied on domain discrete methods, which may not be suitable for BEM. Yet another possibility is related to the fact that the mass matrix of DRBEM is not symmetric, generating complex eigenvalues in an eigenvalue analysis.

More recently, a new BEM formulation was developed based on Poisson's fundamental solution and radial basis functions called DIBEM (Loeffler et al., 2017a). The DIBEM technique, like DRBEM, also approximates non-self-adjoint kernels of domain integrals by sequences of radial basis functions. However, despite being similar to DRBEM in many respects, the DIBEM formulation does not require the construction of two auxiliary matrices that multiply the classical boundary elements matrices H and G since it directly addresses the integral domain, similar to what is done in an interpolation procedure using a primitive radial basis function. Therefore, the entire kernel of the domain integral is interpolated, including the fundamental solution. The DIBEM procedure was successfully applied to solve Poisson (Loeffler et al., 2017a), Helmholtz (Loeffler et al., 2017b), and Advection Diffusion (Pinheiro et al., 2022).

Regarding stationary wave problems, DIBEM's greater robustness and accuracy were demonstrated by comparison with DRBEM results (Balista et al., 2023), highlighting that both formulations use fundamental solutions independent of frequency. DIBEM model can also be successfully applied in consortium with RIM since this method requires an auxiliary technique for solving problems in which the kernel is composed of an unknown variable (Campos et al., 2020). However, in the dynamic response, the instability problem for reduced integration steps persisted for significantly smaller values. Aiming to understand this stability problem, the robustness of the inertia DIBEM matrix was analyzed concerning conditioning and positivity. Several conditioning norms were compared to several meshes constituted by linear boundary elements, used in the solution of classical problems of wave propagation in scalar fields. A positivity study was also undertaken, in addition to curves that relate the degree of mesh discretization as the minimum integration step capable of generating a stable solution.

As for positivity, a large percentage of negative and complex eigenvalues was observed from an analysis of free vibration problems in syncretized systems by DIBEM. Despite DIBEM's superior results in calculating natural frequencies, the high percentage of non-real eigenvalues drew attention to the need for a precise study of its effect on the dynamic response problem. For this purpose, a DIBEM model based on modal superposition was implemented. It should be noted that the use of this technique, considering that the BEM matrices are not symmetric, is very different from what is found in the literature on structural dynamics (Clough and Penzien, 1975), as it uses the Modified Modal Method (MMM), which calculates the eigenvalues of the adjoint problem (Prodonoff and Zepka, 1983). The pioneering work of Nardini and Brebbia (1983a) using DRBEM also uses the MMM and so did several authors (Ouisse and Foltete, 2011; Zhang and Lallement, 1987; Zhang et al., 1988). However, the solution of the adjoint problem doubles the computational time for solving the eigenvalue problem, which could be done using the classical model presented in the linear algebra literature (Friedberg et al., 1992), from now on named CMM, in which the inverse of the matrix composed of the eigenvectors is used.

Modal superposition is an expensive technique, as it involves the calculation of eigenvalues of the discrete system, which in dynamics must necessarily consist of numerous degrees of freedom. Furthermore, it becomes inadequate if there are physical or geometric non-linearities, as natural frequencies change under these conditions. This is why direct integration schemes, in which temporal discretization is done through finite difference-type approximations, are widely used in practice. In direct integration, the solution advances in time through successive increments, of small magnitude, to obtain precision, stability and convergence of the response.

Thus, in this work, an original application is made in which the modal superposition method is applied to DIBEM to verify whether the effect of non-real modes and frequencies in solving the direct problem is effectively harmful, being the factor responsible for the instability of the answer for small integration steps. By comparing results with direct analysis, in which the matrix system is coupled, you can check whether or not the integration of separate equations presents the same pattern of instability for very restricted time increments. This work also compares two modal superposition methods for non-symmetric systems, MMM and CMM, to verify the reason for the massive use of the first in specialized literature without any justification or presentation of comparative results.

2 GOVERNING INTEGRAL EQUATION

Consider a homogeneous domain Ω(X), representing a system, a body, or a control volume of a given physical problem. The differential equation that describes the propagation of an acoustic wave in an isotropic medium is given (Graff, 2012) as follows:

$\mathrm{u}{,}_{\mathrm{i}\mathrm{i}}\left(\mathbf{X}\right)=\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{${\mathrm{k}}^2$}\right.\right)\mathrm{u}\left(\mathbf{X}\right)$ (1)

The scalar k means the speed of propagation of the acoustic wave. Being a scalar equation, the effects of inertia and stiffness contribute to establishing a dynamic equilibrium given only by spherical contents of the stress state or pressure, which can be understood as a potential given u(X,t). The boundary conditions are provided by Dirichlet condition, in which the potentials are prescribed on Γ_u; and the Neumann conditions, given by the first derivative of the potential concerning the external normal direction to the boundary Γ_q, given by q(X,t). Thus, considering the well-known BEM mathematical procedures (Kythe, 1995), Eq. (2) represents the inverse integral form associated with the Laplacian operator on the left-hand side of Eq. (1):

$c\left(\mathit{\xi }\right)u\left(\mathit{\xi }\right)+{\int }_{\Gamma }u\left(\mathit{X}\right)q^{\mathrm{*}}\left(\mathit{\xi };\mathit{X}\right)d\Gamma -{\int }_{\Gamma }q\left(\mathit{X}\right)u^{\mathrm{*}}\left(\mathit{\xi };\mathit{X}\right)d\Gamma =-\frac{1}{k^2}{\int }_{\Omega }\ddot{u}\left(\mathit{X}\right)u^{\mathrm{*}}\left(\mathit{\xi };\mathit{X}\right)d\Omega$ (2)

In Eq. (2), the term X represents the field point, any point related to the domain Ω(X), limited by the boundary Γ(X). The basis point of the integrations is called the source point ξ. The term c(ξ) is a coefficient related to the location of the source point ξ relative to Ω(X) and, considering that it can be located on the boundary Γ(X), inside or outside it. Smoothness also influences the term c(ξ) (Brebbia et al., 1984). It should be noted that u(X) is the scalar potential and q(X) is its normal derivative; u*(ξ;X) is the fundamental solution correlated with Laplace's problem, and q*(ξ;X) is its normal derivative:

$u^{\mathrm{*}}\left(\mathit{\xi };\mathit{X}\right)=-\frac{\mathit{ln}\left[r\left(\mathit{\xi };\mathit{X}\right)\right]}{2\pi }$ (3)

$q^{\mathrm{*}}\left(\mathit{\xi };\mathit{X}\right)=-\frac{1}{2\pi r\left(\mathit{\xi };\mathit{X}\right)}{r}_i\left(\mathit{\xi };\mathit{X}\right){\eta }_i\left(\mathit{X}\right)$ (4)

In Eq. (3) and (4), r(ξ;X) is the Euclidean distance between the source point ξ and any field point X; and n_i(X) is the external normal on the boundary at field point X.

3 DIRECT INTERPOLATION TECHNIQUE

Using the DIBEM technique, all functions comprising the domain integral's kernel on the right-hand side of Eq. (2) are approximated by the sequence of radial functions. Thus, the fundamental solution is also included, which depends on the position of the source point. For this reason, DIBEM presents singularity problems if the source points coincide with the field point X, which is the usual procedure. Thus, to avoid this, the regularization procedure is used as a strategy (Loeffler and Mansur, 2017), adding two new integrals, as follows:

$c\left(\mathit{\xi }\right)u\left(\mathit{\xi }\right)+{\int }_{\mathrm{\Gamma }}u\left(\mathit{X}\right)q^{\mathrm{*}}\left(\mathit{\xi };\mathit{X}\right)d\mathrm{\Gamma }-{\int }_{\mathrm{\Gamma }}q\left(\mathit{X}\right)u^{\mathrm{*}}\left(\mathit{\xi };\mathit{X}\right)d\mathrm{\Gamma }=\frac{1}{k^2}\left[{\int }_{\mathrm{\Omega }}(\ddot{u}\left(\mathit{\xi }\right)-\ddot{u}\left(\mathit{X}\right))u^{\mathrm{*}}\left(\mathit{\xi };\mathit{X}\right)d\mathrm{\Omega }\right]-\frac{1}{k^2}\left[\ddot{u}\left(\mathit{\xi }\right){\int }_{\mathrm{\Omega }}{u}^{\mathrm{*}}\left(\mathit{\xi };\mathit{X}\right)d\mathrm{\Omega }\right]$ (5)

The singularity in the domain integral approximation is removed in equation (5), but an additional domain integral is created. However, it is easily transformed into a boundary integral, through the Galerkin Tensor concept (Wrobel and Aliabadi, 2002):

$\ddot{u}\left(\mathit{\xi }\right){\int }_{\mathrm{\Omega }}{u}^{\mathrm{*}}\left(\mathit{\xi };\mathit{X}\right)d\mathrm{\Omega }=\ddot{u}\left(\mathit{\xi }\right){\int }_{\mathrm{\Omega }}{G}_{,ii}^{\mathrm{*}}\left(\mathit{\xi };\mathit{X}\right)=\ddot{u}\left(\mathit{\xi }\right){\int }_{\mathrm{\Gamma }}{G}_{,i}^{\mathrm{*}}\left(\mathit{\xi };\mathit{X}\right)n_{i}d\mathrm{\Gamma }$ (6)

Where:

$G_{,i}^{\mathrm{*}}\left(\mathit{\xi };\mathit{X}\right)n_i=P_X^{\xi }\left(\mathit{\xi };\mathit{X}\right)=\frac{1}{4\pi }\{0.5-\ln \left[\mathrm{r}\left(\mathit{\xi };\mathit{X}\right)\right]\}{r}_i{n}_i$ (7)

Considering the additional term for regularization, the approximation via radial basis functions Fⁱ(Xⁱ;X) consists of the following expression:

$[\mathrm{}\ddot{u}\left(\mathit{X}\right)-{\ddot{u}\left(\mathit{\xi }\right)]u}^{\mathrm{*}}\left(\mathit{\xi };\mathit{X}\right)\cong {\ddot{\alpha }}^i{F}^i\left({\mathit{X}}^{\mathit{i}};\mathit{X}\right)$ (8)

Therefore, substituting Eq. (8) in the domain integral concerning the first domain integral on the right hand side of Eq. (5), one has:

$\left\{{\int }_{\mathrm{\Omega }}^{}\left[\ddot{u}\left(\mathit{X}\right)-\ddot{u}\left(\mathit{\xi }\right){]u}^{\mathrm{*}}\left(\xi ;\mathit{X}\right)\right]d\Omega \right\}\approx \left\{{\int }_{\mathrm{\Omega }}^{}{}^{\xi }{\ddot{\alpha }}^i{F}^i\left({\mathit{X}}^i;\mathit{X}\right)d\mathrm{\Omega }\right\}=\left\{{\int }_{\mathrm{\Omega }}^{}{}^{\xi }{\ddot{\alpha }}^i{{\psi }_{,kk}}^i\left({\mathit{X}}^i;\mathit{X}\right)d\mathrm{\Omega }\right\}={}^{\xi }{\ddot{\alpha }}^i{\int }_{\mathrm{\Gamma }}^{}{\eta }^i\left({\mathit{X}}^i;\mathit{X}\right)d\mathrm{\Gamma }$ (9)

DIBEM also uses an auxiliary interpolation function Ψⁱ(X^j;X), primitive of Fⁱ(X^j;X), transforming the domain integral into a boundary one. It is interesting to note at this point that the same radial basis functions (RBF) used in DRBEM can be employed. The function ${\eta }^i$ depends on the radial function used. In this article, the simple radial basis function is used, so that the ${\eta }^i$ function is given by:

${\eta }^i\left({\mathit{X}}^i;\mathit{X}\right)=r{(\mathit{X}}^i;\mathit{X})r_i{n}_i$ (10)

In Eq. (10) $r_i$ is the component of the radius vector and $n_i$ is the unit normal vector to the boundary. The surplus integral term can be easily rewritten as a boundary integral:

${\int }_{\mathrm{\Omega }}^{}\left[\ddot{u}\left(\mathit{\xi }\right)u^{\mathrm{*}}\left(\mathit{\xi };\mathit{X}\right)\right]d\mathrm{\Omega }=\ddot{u}\left(\mathit{\xi }\right){\int }_{\mathrm{\Omega }}^{}\left[{G\left(\mathit{\xi };\mathit{X}\right)}_{,ii}^{\mathrm{*}}\right]d\mathrm{\Omega }=\ddot{u}\left(\mathit{\xi }\right){\int }_{\mathrm{\Gamma }}^{}\left[{G\left(\mathit{\xi };\mathit{X}\right)}_{,i}^{\mathrm{*}}{n}_i\right]d\mathrm{\Gamma }=\ddot{u}\left(\mathit{\xi }\right)P_{\xi }$ (11)

Substituting Eq. (9) and (11) into Eq. (5) and using the standard BEM discretization procedures, where “n” source points are composed by the sum of internal values and boundary nodal points, the matrix system given by Eq. (12) is obtained. Details can be found in Loeffler et al.(2015b).

$\left(\begin{array}{ccc}{\mathrm{H}}_{11}& \dots & {\mathrm{H}}_{1\mathrm{n}}\\ ⋮& \ddots & ⋮\\ {\mathrm{H}}_{\mathrm{n}1}& \dots & {\mathrm{H}}_{\mathrm{n}\mathrm{n}}\end{array}\right)\left(\begin{array}{c}{\mathrm{u}}_1\\ ⋮\\ {\mathrm{u}}_{\mathrm{n}}\end{array}\right)-\left(\begin{array}{ccc}{\mathrm{G}}_{11}& \dots & {\mathrm{G}}_{1\mathrm{n}}\\ ⋮& \ddots & ⋮\\ {\mathrm{G}}_{\mathrm{n}1}& \dots & {\mathrm{G}}_{\mathrm{n}\mathrm{n}}\end{array}\right)\left(\begin{array}{c}{\mathrm{q}}_1\\ ⋮\\ {\mathrm{q}}_{\mathrm{n}}\end{array}\right)=\frac{1}{{\mathrm{k}}^2}\left(\begin{array}{ccc}{{}^1\mathrm{\alpha }}^1& \dots & {{}^1\mathrm{\alpha }}^{\mathrm{n}}\\ ⋮& \ddots & ⋮\\ {{}^{\mathrm{n}}\mathrm{\alpha }}^1& \dots & {{}^{\mathrm{n}}\mathrm{\alpha }}^{\mathrm{n}}\end{array}\right)\left(\begin{array}{c}{\mathrm{N}}_1\\ ⋮\\ {\mathrm{N}}_{\mathrm{n}}\end{array}\right)-\frac{1}{{\mathrm{k}}^2}\left(\begin{array}{c}{\mathrm{Z}}_1\\ ⋮\\ {\mathrm{Z}}_{\mathrm{n}}\end{array}\right)$ (12)

The vector Z is composed by the following integrals:

$$\left.\left(\begin{array}{l}{Z}_1\\ Z_2\\ \mathrm{...}\\ Z_n\end{array}\right.\right)=\left.\left(\begin{array}{l}[{\displaystyle \underset{1}{\int }{}^{1}P{}^{1}d{\Gamma }_1}+{\displaystyle \underset{1}{\int }{}^{1}P{}^{2}d{\Gamma }_2}+\mathrm{...}+{\displaystyle \underset{1}{\int }{}^{1}P{}^{n}d{\Gamma }_n}]{{\displaystyle }}^{}0{{\displaystyle }}^{}\mathrm{.....}{{\displaystyle }}^{}0\\ 0{{\displaystyle }}^{}[{\displaystyle \underset{1}{\int }{}^{2}P{}^{1}d{\Gamma }_1}+{\displaystyle \underset{1}{\int }{}^{2}P{}^{2}d{\Gamma }_2}+\mathrm{...}+{\displaystyle \underset{1}{\int }{}^{2}P{}^{n}d{\Gamma }_n}]{{\displaystyle }}^{}\mathrm{.....}{{\displaystyle }}^{}0\\ \mathrm{........}\\ 0{{\displaystyle }}^{}0{{\displaystyle }}^{}\mathrm{...}{{\displaystyle }}^{}[{\displaystyle \underset{1}{\int }{}^{n}P{}^{1}d{\Gamma }_1}+{\displaystyle \underset{1}{\int }{}^{n}P{}^{2}d{\Gamma }_2}+\mathrm{...}+{\displaystyle \underset{1}{\int }{}^{n}P{}^{n}d{\Gamma }_n}]\end{array}\right.\right)\left.\left(\begin{array}{l}{u}_1\\ u_2\\ \mathrm{...}\\ u_n\end{array}\right.\right)$$ (13)

To explicit the potential in the discretized governing equation (see Eq. (12)), it is necessary to make the following substitution, considering new coefficients A_ξ which is given by:

$$A_{\xi }=\left.\left(N_1{{\displaystyle }}^{}{N}_2\mathrm{....}{{\displaystyle }}^{}{N}_m\right.\right)\left.\left(\begin{array}{l}{}^{\xi }\alpha {}_1\\ \mathrm{...}\\ {}^{\xi }\alpha {}_m\end{array}\right.\right)$$ (14)

In turn, the ^ξα can be calculated using the basic interpolation formula, that is:

$$[{}^{\xi }\alpha ]={[F]}^{-1}[{}^{\xi }\Lambda ][F]\alpha ={[F]}^{-1}[{}^{\xi }\Lambda ][u]$$ (15)

Thus, in this case:

$$\begin{array}{l}{A}_{\xi }=\left.\left(N_1{{\displaystyle }}^{}{N}_2\mathrm{....}{{\displaystyle }}^{}{N}_n\right.\right){\left.\left(\begin{array}{l}{F}^{11}\mathrm{...}{{\displaystyle }}^{.}{F}^{1n}\\ \mathrm{.....}\\ F^{n1}{{\displaystyle }}^{}\mathrm{...}{F}^{nn}\end{array}\right.\right)}^{-1}\left\{\left.\left(\begin{array}{l}{}^{\xi }\Lambda {}^1\mathrm{...}{{\displaystyle }}^{}0\\ \mathrm{......}\\ 0{{\displaystyle }}^{}\mathrm{...}{}^{\xi }\Lambda {}^n\end{array}\right.\right)\left.\left(\begin{array}{l}{u}_1\\ \\ u_n\end{array}\right.\right)-\left.\left(\begin{array}{l}{}^{\xi }\Lambda {}^1.u_{\xi }\mathrm{..}{{\displaystyle }}^{}0\\ \mathrm{......}\\ 0{{\displaystyle }}^{}\mathrm{...}{}^{\xi }\Lambda {}^m{u}_{\xi }\end{array}\right.\right)\right\}=\\ \left.\left(S_1{{\displaystyle }}^{}{S}_2\mathrm{....}{{\displaystyle }}^{}{S}_n\right.\right)\left\{\left.\left(\begin{array}{l}{}^{\xi }\Lambda {}^1\mathrm{...}{{\displaystyle }}^{}0\\ \mathrm{......}\\ 0{{\displaystyle }}^{}\mathrm{...}{}^{\xi }\Lambda {}^n\end{array}\right.\right)\left.\left(\begin{array}{l}{u}_1\\ \\ u_n\end{array}\right.\right)-\left.\left(\begin{array}{l}{}^{\xi }\Lambda {}^1.u_{\xi }\mathrm{..}{{\displaystyle }}^{}0\\ \mathrm{......}\\ 0{{\displaystyle }}^{}\mathrm{...}{}^{\xi }\Lambda {}^n{u}_{\xi }\end{array}\right.\right)\right\}\end{array}$$ (16)

After applying this procedure, it is possible to write the domain integral of the term related to inertia only in terms of integral involving boundary variables. Therefore, the final system can be written as follows:

$\left[\mathrm{H}\right]\left\{\mathbf{u}\right\}-\left[\mathbf{G}\right]\left\{\mathbf{q}\right\}=\left[\mathrm{M}\right]\left\{\ddot{\mathbf{u}}\right\}$ (17)

In Eq. (17), the matrices [H] and [G] are standard matrices that appear from the discretization operations involving the integral of the fundamental solution and its normal derivative, respectively. The matrix [M] corresponds to the acoustic inertia property, while {u} and {q} are vectors that contain the values of the potential and its derivative at the nodes.

4 THE EIGENVALUE PROBLEM

The Helmholtz equation can be interpreted as a simplification of the scalar wave equation (see Eq. (1)), when searching for a response produced in the system by an excitation variable whose frequency ω is known. In this condition, the potential U(X,t) is equivalent to a periodic response, composed of specific temporal content:

$$U(X,t)=u(X)e^{-i_{}{\omega }_{}t}$$ (18)

In Eq. (18), "i" is the imaginary unit and u(X) is the space amplitude of the stationary system response to harmonic excitation, with frequency ω. Thus, considering the free vibration condition, one has:

$$\ddot{U}(X,t)=-{\omega }^{2}U(X,t)$$ (19)

However, eigenvalue BEM matrices involve stationary amplitudes of potential U(X) and normal derivatives Q(X), such as the equations shown in Eq. (17) need to be adequately manipulated to formulate the eigenvalue problem, in which all natural frequencies ω can be calculated. In this sense, the prescribed u and q nodal values, which now mean amplitudes, depending just on the position X, are highlighted below:

$$\left.\left(\begin{array}{l}{H}_{u\overline{u}}{{\displaystyle }}^{}{H}_{u\overline{q}}\\ H_{q\overline{u}}{{\displaystyle }}^{}{H}_{q\overline{q}}\end{array}\right.\right)\left.\left(\begin{array}{l}\overline{U}\\ U\end{array}\right.\right)-\left.\left(\begin{array}{l}{G}_{u\overline{u}}{{\displaystyle }}^{}{G}_{u\overline{q}}\\ G_{q\overline{u}}{{\displaystyle }}^{}{G}_{q\overline{q}}\end{array}\right.\right)\left.\left(\begin{array}{l}Q\\ \overline{Q}\end{array}\right.\right)={\omega }^2\left.\left(\begin{array}{l}{M}_{u\overline{u}}{{\displaystyle }}^{}{M}_{u\overline{q}}\\ M_{q\overline{u}}{{\displaystyle }}^{}{M}_{q\overline{q}}\end{array}\right.\right)\left.\left(\begin{array}{l}\overline{U}\\ U\end{array}\right.\right)$$ (20)

It must be highlighted that the improvement of the DIBEM model is given not only by the mesh refinement but also with the insertion of internal source points, usually called poles, which improve the consistency of the DIBEM model, introducing a better approximation of inertia properties inside. Considering that the prescribed nodal values of U(X) and Q(X) are zero null for free vibration analysis, one has:

$$H_{u\overline{q}}U-G_{u\overline{u}}Q={\omega }^2{M}_{u\overline{q}}U$$ (21)

$$H_{q\overline{q}}U-G_{q\overline{u}}Q={\omega }^2{M}_{q\overline{q}}U$$ (22)

Eliminating the normal derivative of potential Q from Eq. (21) and Eq. (22), a typical eigenvalue matrix system is obtained:

$$\left[\overline{\text{H}}\right]\left\{\text{U}\right\}={\omega }^2\left[\overline{\text{M}}\right]\left\{\text{U}\right\}$$ (23)

In previous equation new BEM matrices are comprised by:

$\begin{array}{l}\left[\overline{\text{H}}\right]=\left[{\text{H}}_{\text{qq}}\right]-\left[{\text{G}}_{\text{qu}}\right]{\left[{\text{G}}_{\text{uu}}\right]}^{-1}\left[{\text{H}}_{\text{uq}}\right]\\ \left[\overline{\text{M}}\right]=\left[{\text{M}}_{\text{qq}}\right]-\left[{\text{G}}_{\text{qu}}\right]{\left[{\text{G}}_{\text{uu}}\right]}^{-1}\left[{\text{M}}_{\text{uq}}\right]\end{array}$ (24)

New BEM matrices H and M can be interpreted as mass and stiffness matrices, which are asymmetric and fully populated. Therefore, the effectiveness of the BEM eigenvalue problem hinges on the availability of an efficient eigenvalue extraction method. Here, it was used the Hessenberg algorithm. Concerning the eigenvalues, the asymmetry of BEM matrices can produce negative and complex frequencies associated with inaccuracies in the numerical model. Such non-real eigenvalues, however, are related to the approximation of higher frequencies in the spectrum (Loeffler and Mansur, 1987). The improvement of the BEM model, through the mesh refinement and insertion of poles, minimizes this problem but does not eliminate it (Barbosa et al., 2019). Comparisons performed with DRBEM have shown that DIBEM model results are meaningful superior, improving the quality of the Helmholtz problem's direct solution and producing eigenvalues calculation with more accuracy.

It also must be pointed out that the mathematical steps of deduction of the integral BEM equation for the Helmholtz equation are different if the strictly correlated Green fundamental solution is employed. In this case, a transcendental characteristic equation is obtained, in which the matrix system appears nonlinear, depending on the frequency. In the case of eigenvalue problems, this naturally demands using iterative solution schemes.

Eigenvalues such as presented in Eq. (23). However, in order to apply the modal superposition method with the BEM, a dynamic [D] matrix is usually generated, in the following form:

$$\text{[D]}\left\{u\left(X\right)\right\}={\text{[M]}}^{\text{-1}}\text{[H]}\left\{u\left(X\right)\right\}={\omega }^2\text{[I]}\left\{u\left(X\right)\right\}$$ (25)

The eigenvalues are obtained relative to the dynamic matrix. If the stiffness and mass matrices were symmetric, such an operation would destroy the symmetry of the matrices, which is why it is not used with the Finite Element Method. However, the BEM matrices are asymmetric. Such an inversion, which implies an additional cost of the inversion, is not essential, but in this context it will greatly facilitate the understanding and comparison of the schemes shown below.

5 MODAL SUPERPOSITION

Modal superposition was a technique widely used in structural analysis in the 1970s. However, its use in standard form declined later due to computational cost issues and also because it was prohibitive in solving physically and geometrically non-linear problems. Several adaptations of the method were proposed in order to reduce the computation time by a selective approach to the load spectrum of frequencies or neglecting the participation of high modes. As mentioned, the modal superposition was performed in elastodynamics with DRBEM. Still, the results were unsatisfactory, corroborating that DRBEM does not work correctly in dynamics problems, whether in acoustics or elasticity, because its inertia matrix is not robust enough.

To present the method to dynamics, it is necessary to consider the matrix system referring to the acoustic wave propagation model, see Eq. (17), and then partition the matrix system, as shown in Eq. (23). Due to the absence of time derivatives of tractions it is possible to rewrite the BEM system in a similar way to other methods, that is:

$\left[H\mathrm{\text{'}}\right]\left\{u\left(X,t\right)\right\}+\left[M\mathrm{\text{'}}\right]\left\{\ddot{u}\left(X,t\right)\right\}=\left\{f\left(X,t\right)\right\}$ (26)

In the previous Equation, the vector {f} contains the prescribed values. From now on the arguments will be omitted for simplicity. Then, a coordinate transformation is performed using the matrix of normalized automatic vectors [φ(X)]:

$\left\{u\right\}=\left[\varphi \right]\left\{y\right\};\mathrm{}\ddot{u}\left\{X\right\}=\left[\varphi \right]\left\{\ddot{y}\right\}$ (27)

In the previous equation, {y} is the nodal vector of modal coordinates. Multiplication by [M’]^-1 generates the dynamic matrix:

$\left[D\right]\left[\varphi \right]\left\{y\right\}+\left[I\right]\left[\varphi \right]\left\{\ddot{y}\right\}={\left[M\mathrm{\text{'}}\right]}^{-1}\left\{f\right\}$ (28)

Pre-multiplication by the inverse of [φ] diagonalizes the system:

${\left[\varphi \right]}^{-1}\left[D\right]\left[\varphi \right]\left\{y\right\}=\left[\wedge \right]\left\{y\right\}$ (29)

The diagonal matrix [Λ] is composed of the various natural frequencies associated with the discrete system, now known. The new traction vector {p} is given by:

$\left\{p\right\}={\left[\varphi \right]}^{-1}{\left[M\mathrm{\text{'}}\right]}^{-1}\left\{f\right\}$ (30)

The equation of motion is decoupled, consisting of n ordinary differential equations, where n is the number of degree of freedom, as shown below:

$\left\{\ddot{y}\right\}+\left[\wedge \right]\left\{y\right\}=\left\{p\right\}$ (31)

Solution of Eq. (31) for each mode is given by:

$y_i\left(t\right)=\frac{p_i}{{\omega }_i^2}\mathrm{}[1-\mathrm{c}\mathrm{o}\mathrm{s}({\omega }_i\mathrm{}\mathrm{t}\left)\right]$ (32)

In the case of using the solution of the adjoint problem, the eigenvectors referring to the transpose of the dynamic matrix must be calculated. Its transpose has the same eigenvalues, but its eigenvectors [φ*] are different. It should be noted that the eigenvector matrix of the original problem [φ'] is not normalized in the usual sense, that is, it does not have a unitary module. Thus, the modal coordinates are different:

$\left\{u\right\}=\left[\varphi \mathrm{\text{'}}\right]\left\{y\mathrm{\text{'}}\right\}$ (33)

However, it is demonstrated that there is a cross-orthogonality relationship (Atkinson, J. H., 1965; Prodonoff and Zepka, 1983), in the sense that for each specific frequency the modes φ' and φ* are orthogonal. Based on this it can also be demonstrated that:

${\left[\varphi \mathrm{\text{'}}\right]}^T\left[{\varphi }^{\mathrm{*}}\right]=\left[I\right]$ (34)

For the previous identity matrix to be obtained, a particular metric must be used in the eigenvector matrix. Obeying this condition, one has:

${\left[\varphi \mathrm{\text{'}}\right]}^T\left[D\right]\left[{\varphi }^{\mathrm{*}}\right]\left\{y^{\mathrm{*}}\right\}=\left[\wedge \right]\left\{y^{\mathrm{*}}\right\}$ (35)

Thus, the matrix governing equation is given by:

${\left\{{\ddot{y}}^{\mathrm{*}}\right\}+\left[\wedge \right]\left\{y^{\mathrm{*}}\right\}=\left[\varphi \mathrm{\text{'}}\right]}^T{\left[M\mathrm{\text{'}}\right]}^{-1}\left\{f\right\}=\left\{g\right\}$ (36)

In this way, observing Eq. (33) and (36) the two methods can be compared. It can be seen that ordinary differential equations have a similar structure, but their modal coordinates are different, and, mainly, the forcing terms are completely distinct. In the CMM, this term is obtained by multiplying the inverse of the eigenvector matrix. In the MMM, the forcing term involves only the matrix product of the eigenvectors of the transposed matrix. The adjoint matrix method requires the solution of two eigenvalue problems and is therefore much more computationally expensive. However, the literature that applies the modal method uses it widely, without making any mention of its possible advantages compared to the classical method. Given this gap, the numerical simulations presented in this work will also include a comparison of these two methods, showing that they are radically different and the MMM is much more accurate, something not presented in the literature.

6 TRACTIONS CALCULATION

The BEM is a mixed numerical technique in which the primal variable and its normal derivative are calculated simultaneously. In the modal method, the matrix system is only arranged as a function of displacements. Thus, it is necessary to adapt the BEM matrices to use the displacement values already calculated. It is done through a new condensation of the BEM matrices to make the tractions explicit, considering the partitioning given by Eq. (37) as follows:

$\left[\mathrm{G}\mathrm{B}\right]\left\{{\mathrm{Q}}_u\right\}=\left[\mathrm{G}\mathrm{C}\right]\left\{{\mathrm{Q}}_q\right\}+\left[\mathrm{M}\mathrm{C}\right]\left\{{\mathrm{Q}}_u\right\}+\left[\mathrm{H}\mathrm{C}\right]\left\{{\mathrm{Q}}_u\right\}+\left[\mathrm{H}\mathrm{D}\right]\left\{{\mathrm{Q}}_u\right\}=\left\{\mathrm{Z}\left(\mathrm{t}\right)\right\}$ (37)

The matrices shown in Eq. (37) are composed of submatrices based on Eq. (20), that is:

$\left[\mathrm{G}\mathrm{B}\right]=\left[{\mathrm{G}}_{uu}\right]+\left[{\mathrm{M}}_{uq}\right]{\left[{\mathrm{M}}_{qq}\right]}^{-1}\left[{\mathrm{G}}_{qu}\right]$ (38)

$\left[\mathrm{G}\mathrm{C}\right]=\left[{\mathrm{M}}_{uq}\right]\left[{\mathrm{M}}_{qq}\right]\left[{\mathrm{G}}_{qu}\right]-\left[{\mathrm{G}}_{uu}\right]$ (39)

$\left[\mathrm{M}\mathrm{C}\right]=\left[{\mathrm{M}}_{uu}\right]+\left[{\mathrm{M}}_{uq}\right]{\left[{\mathrm{M}}_{qq}\right]}^{-1}\left[{\mathrm{M}}_{qu}\right]$ (40)

$\left[\mathrm{H}\mathrm{C}\right]=\left[{\mathrm{H}}_{uu}\right]+\left[{\mathrm{M}}_{uq}\right]{\left[{\mathrm{M}}_{qq}\right]}^{-1}\left[{\mathrm{H}}_{qu}\right]$ (41)

$\left[\mathrm{H}\mathrm{D}\right]=\left[{\mathrm{H}}_{uq}\right]+\left[{\mathrm{M}}_{uq}\right]{\left[{\mathrm{M}}_{qq}\right]}^{-1}\left[{\mathrm{H}}_{qq}\right]$ (42)

Then, one has:

$\left[\mathrm{G}\mathrm{B}\right]\left\{{\mathrm{Q}}_u\right\}=\left\{\mathrm{Z}\left(\mathrm{t}\right)\right\}\to \left\{{\mathrm{Q}}_u\right\}={\left[\mathrm{G}\mathrm{B}\right]}^{-1}\left\{\mathrm{Z}\left(\mathrm{t}\right)\right\}$ (43)

It may be possible to calculate the tractions where the u values are prescribed using another methodology, but the proposal presented here is original and has not been shown in any previous work. It must be considered that all applications of the modal superposition method were made using the FEM or the techniques of matrix analysis of structures, according to a displacement formulation, without dealing with the calculation of stresses. It should be noted that the Q_u values are obtained in the parts of the boundary where the potential was prescribed; therefore, the number of degrees of freedom involved in determining stresses is much smaller than that used in characterizing displacements.

The tractions are strongly influenced by the high natural frequencies (Clough and Penzien, 1975; Loeffler and Mansur, 1987), which are poorly represented by any discrete matrix system. This leads to inaccuracies in the response of the tractions that appear as spurious fluctuations of small amplitude along the time.

7 TIME MARCHING SCHEME

Among the first conclusions obtained in the dynamics simulations with DRBEM is the need for an unconditionally stable time advance scheme. However, BEM is a mixed method in which potentials and their normal derivatives are calculated simultaneously. In dynamics, it is known that the derivatives of the potential, which represent the tractions, are strongly affected by the higher frequencies, which in turn are poorly calculated. Therefore, it is necessary to eliminate the contribution of the high-frequency spectrum, which can be done through incremental schemes equipped with fictitious damping. In this context, the Houbolt scheme stands out for its substantial incidence of numerical dissipation and has been widely used with DRBEM.

On the other hand, as presented in the literature, the decoupled ordinary equations obtained by the modal superposition method can be solved analytically through the well-known Duhamell integral or numerically through an incremental time advance scheme. Aiming to compare the results of this article with those obtained by direct analysis, the Houbolt integration scheme is used to get an answer over time. Also, through the incremental scheme, the integration steps could be different for each equation, with the equations relating to high modes being imprecise and sensitive to the instability and precision of minimal steps, unlike the equations relating to more fundamental modes. However, a single value for the time step was used to integrate all equations.

Thus, the Houlbolt time stepping scheme approximates the acceleration as follows:

${\ddot{\mathit{u}}}_{\mathit{n}+1}=\mathrm{}\frac{2{\mathit{u}}_{\mathit{n}+1}-5{\mathit{u}}_{\mathit{n}}+4{\mathit{u}}_{\mathit{n}-1}-{\mathit{u}}_{\mathit{n}-2}}{{∆t}^2}$ (44)

In the last formula, Δt is the discretization time interval and n are the temporal instants. Houbolt scheme is unconditionally stable, but instability was observed with the reduction of the time interval in the computational DIBEM experiments using direct response with full matrices.

8 NUMERICAL SIMULATIONS

8.1 First test

Just to illustrate the differences between the methods, two non-symmetric matrices are taken as hypothetical stiffness and inertia matrices, whose composition is shown in Figure 1, with six degrees of freedom. The system is subject to a strong term capable of exciting all of its six natural modes of response.

Figure 1
Hipotheticous discrete system.

The CMM and MMM methods are used to calculate frequencies and obtain the response over time, which was obtained considering different modal contents, as seen in Figures 2, 3, and 4:

Figure 2
Time response using one mode.

Figure 3
Time response using three modes.

Figure 4
Time response using six modes.

It is noted that the final result was the same. Still, the CMM constructs the response based on modal contents referring to the highest frequencies, which is an unusual procedure and has profound implications. This distinctive pattern of behavior is due to the multiplication of the forcing term by the inverse of the matrix composed of the eigenvectors. Naturally, these results imply severe numerical issues, as shown in the following example, whose dimension is more significant and represents a physical problem. This occurs because, unlike the lowest modes, the highest modes and their respective frequencies are not well represented by any discrete numerical method.

8.2 Second test

To verify the effectiveness of the Modified Modal Method in solving dynamic problems with DIBEM, computational tests were carried out on the solution of the system commonly known as embedded bar, which has analytics and is widely explored in the specialized literature. Figure 5 illustrates the geometry of the physical problem and its boundary conditions, which have unitary dimensions and unitary wave propagation velocity.

Figure 5
Geometry of the physical problem and its boundary conditions.

The initial conditions are also null. 160 linear boundary elements of equal size and double nodes at the corners were used, with a displaced source point. 144 internal interpolating points were allocated regularly. The integration step used was equal to 0.01s, and the size of each boundary element is 0.0286 units long. According to the CFL criterion, used as step estimation classifications in numerical domain methods, such as the FEM and the Finite Difference Method, this value (or even smaller values), would be sufficient for a good integration without losing stability in a direct integration process.

In this problem, all odd modes are excited but composing an uneven response. As in the Variable Separation Method (Kreiszing, 2015), the modal superposition procedure allows evaluation of the importance of each vibrational mode in the dynamical response. Therefore, the contribution of two closely related parameters can be investigated: the frequency of each mode and the energy provided by all modes in the response. The amplitudes of the natural vibration modes are arbitrary. Still, for the set of initial conditions and typology of the forcing term, the ratio between the amplitudes is directly related to the level of influence of each mode. The most excited modes have greater vibration amplitudes, and their contribution to the final configuration of the system's motion is more significant. Therefore, use amplitude as a metric of the relevance of each mode in the composition of the system response, that is:

$A_i=\left|f_i/{\omega }_i^2\right|$ (45)

To better understand the hierarchy between the modes, all modes are normalized with respect to the fundamental mode. The result is shown in Figure 6.

Figure 6
Normalized modes with respect to the fundamental mode.

The results in Figure 6 clearly show the preponderance of the first modes over those associated with higher frequencies. Note that the amplitude value for the first mode is much greater than for all other modes.

Figure 7 shows the amplitudes of the first three modes in the displacements where the traction is applied. As the Houbolt integration scheme was used, the effect of fictitious damping can now be visually seen in the attenuation of the amplitudes of the third mode and the lengthening of its period. A reduction in the step value would reduce this effect.

Figure 7
Amplitudes of the first three modes.

Figures 8, 9, and 10 show the displacement response where the traction is applied as a function of an increasing number of vibration modes, that is, adding up the responses obtained for each decoupled in the dynamic solution. It is noted that with 40 modes, the numerical response curve practically coincides with the analytical response.

Figure 8
MMM time response for displacements and tractions considering one, third and seventh modes.

Figure 9
MMM time response for displacements and tractions considering twenty modes.

Figure 10
MMM time response for displacements and tractions considering eighty modes.

It is worth mentioning that the same problem was solved for various time increment values, reaching a value equal to 0.00001s without any issue of numerical instability. The graphs for these smaller time values were not shown because there is no noticeable change in what is shown in Figure 10. Obviously, the processing time in these extreme cases is prohibitive. Using the same mesh employing the direct integration method, in which modal decoupling is not performed; the minor possible integration step to obtain a stable response was 0.001s. This fact is robust evidence that the asymmetry of DIBEM dynamic matrices produces instability problems. This means that non-real frequencies, which implicitly are present in DIBEM dynamic matrix restrict the minimum integration time step for integration.

It should be noted that non-real frequencies only have physical meaning if the boundary element formulation using the complex fundamental solution is used. In this case, specific complex values may be associated with dispersive effects. However, the DIBEM uses the fundamental Poisson solution, which is real, so the complex frequencies are due exclusively to the non-symmetry of the BEM matrices. For example, the eigenvalue problem solved by the Finite Element with its symmetric matrices produces only real values, even for the high-frequency spectrum, which generally has low precision.

Regarding the traction behavior, spurious fluctuations are expected due to the influence of higher modes. However, the Gibss phenomenon is also present in the response, and its effect remains in the time response indefinitely unless that time marching scheme has fictitious damping. In addition to the factors mentioned previously, the traction design is also less precise because it has fewer degrees of freedom for its solution, as the modal method is a technique aimed at displacement calculations.

8.3 Third test

This is a solution to the same problem, but examining the answer obtained by the CMM. In the following figures, the dynamic response for displacements is presented as a function of the number of modes. The response with one mode, twenty modes, 200 modes and, 245 modes, are presented respectively in Figures 11, 12, 13 and 14.

Figure 11
CMM dynamic response for displacement calculated with one mode.

Figure 12
CMM dynamic response for displacement calculated with 20 modes.

Figure 13
CMM dynamic response for displacement calculated with 200 modes.

Figure 14
CMM dynamic response for displacement calculated with 245 modes.

The response becomes closer to analytical only when using the totality of natural modes present. The system has a total of 268 degrees of freedom. However, this response could not be improved from this value, probably due to numerical inaccuracy. It should be noted that the sensitivity of this method is very high. The same problem, if discretized with an irregular arrangement of internal points, presents a much lower quality of response.

A final experiment was carried out using a less refined discretization, with 40 linear boundary elements and 64 interpolating points inside. The free displacements at the end for two different modal contents are shown in Figures 15 and 16. The integration step was the same, equal to 0.01s. Although the amplitude obtained has improved, the curves have become noisier due to the lower quality of higher modes. Figures 15 and 16 show the response for 20 and 70 modes. This suggests an additional conclusion in the sense that CMM is very sensitive to the mesh parameters used.

Figure 15
CMM dynamic response for a simpler mesh with 20 modes.

Figure 16
CMM dynamic response for a simpler mesh with 70 modes.

9 CONCLUSION

This work is part of a line of research in which the DIBEM formulation is developed to solve acoustic wave propagation problems. One of the challenges in this field is identifying the factor responsible for the loss of stability of the direct response model for reduced time increments. This behavior is not presented in other numerical methods that discretize the domain but which, in turn, work with symmetric matrices. With DIBEM, using very limited steps direct time response model is subject to working with more refined meshes. Many factors have already been tested: the type of time advancement scheme, the effect of internal points, and the type of utilized radial function, among others. Despite the improved robustness of the DIBEM acoustic inertia matrix, a definitive conclusion has not yet been reached. Thus, the modal superposition technique was required since this method can easily identify and eliminate non-real frequencies in the response after solving the eigenvalue problem. A decoupled scheme of ordinary equations consisting of real frequencies is then created, which can be solved analytically or using incremental time advance schemes.

Based on the experiments developed here, a robust conclusion can be reached in the sense that the complex frequencies, implicitly represented in the DIBEM dynamic matrix, are responsible for the instability of the model for reduced time steps using the direct integration method. In simulations using a less refined discretization, the instability problem must be credited to the bad quality of modal content. When there is a refinement of the mesh, the composition of frequencies becomes more accurate, and thus, smaller steps can be used to improve the accuracy of the numerical solution; however, for small time steps, the effect of fictitious damping in eliminating spurious modes and frequencies is also reduced, causing instability.

Thus, using the modal superposition procedure non-real frequencies were eliminated and time marching schemes that have not fictitious damping could be used, such as the Newmark. The use of the Houbolt scheme was solely due to its ability to dampen the action of spurious high modes in the response and to be able to make a comparison with the results previously obtained by the DIBEM using the direct integration method.

On the other hand, in the studies that were carried out to use the modal method with BEM, two alternatives were presented: the possibility of using the classical method available in the literature and the technique that is used in the main references that deal with non-symmetric matrices called the modified modal method. The reason for choosing the latter, which is more computationally expensive, is not shown in any of the bibliographical references consulted. Both models are mathematically consistent. However, the reason for the choice was computationally proven: although the classical method is simple, it involves multiplying the vector of applied forces by the inverse of the eigenvector matrix. This results in physical behavior utterly opposite to what was expected. It differs entirely from the results obtained by the modal analysis of systems composed of symmetric matrices and from the results presented by the Variable Separation Method in solving partial differential equations of acoustic wave problems. Using a classical algebra strategy, there is a valorization of high vibrational modes, numerically poorly represented, which harm the dynamic response, as shown in this work.

Regarding future works, the results presented here show that it is suitable to perform a review of the techniques related to the symmetrization of mass matrices and also to avail the application of dynamic condensation schemes with BEM, aiming to eliminate the contribution of non-real frequencies within the modal content, if the direct integration scheme is chosen for the numerical solution of the wave equation.

References

• Atkinson, J. H. (1965). The Algebraic Eigenvalue Problem, Claredon Press, Oxford.
• Balista, T.G., Loeffler, C.F., Lara, L.O.C., Mansur, W.J. (2023). Comparisons Between Direct Interpolation and Reciprocity Techniques ofthe Boundary Element Method for Solving Two-Dimensional Helmholtz. Engineering Computations, vol. 40 (9/10), 2841-2861. https://doi.org/10.1108/EC-06-2023-0290
  » https://doi.org/10.1108/EC-06-2023-0290
• Barbosa, J.P., Loeffler, C.F. Lara, L.O.C. (2019). The Direct Interpolation Boundary Element Technique Applied to Three-Dimensional Scalar Free Vibration Problems. Engineering Analysis with Boundary Elements, 108, 295-300. https://doi.org/10.1016/j.enganabound.2019.09.002
  » https://doi.org/10.1016/j.enganabound.2019.09.002
• Bebendorf, M. (2000). Approximation of boundary element matrices. Numer math., 86 (4), 565-589. https://doi.org/10.1007/PL00005410
  » https://doi.org/10.1007/PL00005410
• Bebendorf, M., Grzhibovskis, R. (2006). Accelerating Galerkin BEM for Linear Elasticity Using Adaptive Cross Approximation. Mathematical Methods Applied Sciences, 29(14), 1721-1747. https://doi.org/10.1002/mma.759
  » https://doi.org/10.1002/mma.759
• Bebendorf, M., Rjasanow, S. (2003). Adaptive low-rank approximation of collocation matrices. Computing., 70 (1), 1-24. https://doi.org/10.1007/s00607-002-1469-6
  » https://doi.org/10.1007/s00607-002-1469-6
• Brebbia, C.A., Telles, J.C.F., Wrobel, L.C. (1984). Boundary Element Techniques, First Ed., Berlin Heidelberg: Springer-Verlag.
• Buhmann, M.D. (2003). Radial Basis Function: Theory and Implementations. United Kingdom: Cambridge University Press.
• Campos, L.S., Loeffler, C.F., Netto, F.O., dos Santos, A.J. (2020). Testing the accomplishment of the radial integration method with the direct interpolation boundary element technique for solving Helmholtz problems. Engineering Analysis with Boundary Elements, 110, 16–23.
• Castillo, G.A.V., Loeffler, C.F. (2003). Performance Evaluation of Some New Time Integration Methods in Elastodynamic Problems Formulated by Dual Reciprocity Boundary Element Method. WIT Transactions on Modelling and Simulation, 35, 329-338.
• Chai, P., Zhang, J., Xiao, R., He, R., Lin, W. (2022). A multi-domain BEM on dual interpolation boundary face method for 3D elasticity problem. Engineering Analysis with Boudary Elements, 143, 568-578.
• Clough, R., Penzien, J. (1975). Dynamics of structures. McGraw-Hill Kogakusha, NY, USA.
• Friedberg, S.H., Insel, A.J., Spence, L.E. (1992). Linear Algebra, Prentice Hall.
• Gao, X. (2002). The radial integration method for evaluation of domain integrals with boundary-only discretization. Engineering Analysis with Boundary Elements, 26, 905–916.
• Graff, K.K.F. (2012). Wave motion in elastic solids. Dover Publications.
• Greengard, L., Rokhlin, V. (1987). A fast algorithm for particle simulations. Journal of Computational Physics, 73 (2), 325-348. https://doi.org/10.1016/0021-9991(87)90140-9
  » https://doi.org/10.1016/0021-9991(87)90140-9
• Kirkup, S.M., Kolbrek, B., Yazdani, j., Thompson, A. (2013). Simulation of the Acoustic Field of a Horn Loudspeaker by the Boundary Element–Rayleigh Integral Method. Journal of Computational Acoustics, 21(1). https://doi.org/10.1142/S0218396X12500208
  » https://doi.org/10.1142/S0218396X12500208
• Kreiszing, E. (2015). Advanced Engineering Mathematics, Wiley.
• Kythe, P. (1995). An Introduction to Boundary Element Methods. [S.l.]: Taylor & Francis.
• Liu, Y. (2006). A new fast multipole boundary element method for solving large-scale two-dimensional elastostatic problems. International Journal for Numerical Methods in Engineering, 65(6), 863–881.
• Loeffler, C.F., Barcelos, H.M., Mansur, W.J., Bulcão A. (2015b). Solving Helmholtz Problems with the Boundary Element Method Using Direct Radial Basis Function Interpolation. Engineering analysis with Boundary Elements, 61, 218–25.
• Loeffler, C.F., Cruz, A.L., Bulcão A. (2015a). Direct Use of Radial Basis Interpolation Functions for Modelling Source Terms with the Boundary Element Method. Engineering Analysis with Boundary Elements, 50, 97-108. https://doi.org/10.1016/j.enganabound.2014.07.007
  » https://doi.org/ https://doi.org/10.1016/j.enganabound.2014.07.007
• Loeffler, C.F., Mansur, W.J. (1987). Analysis of Time Integration Schemes for Boundary Element Applications to Transient Wave Propagation Problems. Boundary Element Techniques: Applications in Stress Analysis and Heat Transfer, Computational Mechanics Pub., 105-124.
• Loeffler, C.F., Mansur, W.J. (2017). A regularization scheme applied to the direct interpolation boundary element technique with radial basis functions for solving eigenvalue problem. Engineering Analysis with Boundary Elements, 74, 14–18.
• Loeffler, C.F., Pereira, P.V., Lara, L.O.C., Mansur, W.J. (2017b). Comparison between the formulation of the boundary element method that uses fundamental solution dependent of frequency and the direct radial basis boundary element formulation for solution of Helmholtz problems. Engineering Analysis with Boundary Elements, 79, 81-87. doi: https://doi.org/10.1016/j.enganabound.2017.02.014
  » https://doi.org/10.1016/j.enganabound.2017.02.014
• Loeffler, C.F., Zamprogno, L., Mansur, W.J., Bulcão, A. (2017a). Performance of Compact Radial Basis Functions in the Direct Interpolation Boundary Element Method for Solving Potential Problems. Computer Modeling in Engineering & Sciences, 113(3),367-387.
• Mansur, W.J., Brebbia, C.A. (1985). Further Developments on the Solution of the transient scalar wave equation. Topics in Boundary Element Research, 87-123.
• Nardini, D., Brebbia, C.A. (1983a). A new approach to free vibration analysis using boundary elements. Applied mathematical modelling, 7(3), 157–162. https://doi.org/10.1016/0307-904X(83)90003-3
  » https://doi.org/10.1016/0307-904X(83)90003-3
• Nardini, D., Brebbia, C.A. (1983b). Transient Dynamic Analysis by the Boundary Element Method, Proceedings of the 5th International Conference (BEM V), C.A. Brebbia Editor, Springer Verlag, Berlin, 719-730.
• Narvaez, A, Useche, J. (2022). New Radial Basis Integration Method Applied to the Boundary Element Analysis of 2D Scalar Wave equations. Engineering Analysis with Boundary Elements, 136, 77-92. https://doi.org/10.1016/j.enganabound.2021.12.005
  » https://doi.org/10.1016/j.enganabound.2021.12.005
• Ouisse, M., Foltete, E. (2011). On the properness condition for modal analysis of non symmetric second order systems. Mechanical Systems and Signal Processing, 25 (2), 601-620. https://doi.org/10.1016/j.ymssp.2010.08.017
  » https://doi.org/10.1016/j.ymssp.2010.08.017
• Pinheiro, V.P., Loeffler, C.F., Mansur, W.J. (2022). Boundary element method solution of stationary advective–diffusive problems: A comparison between the direct interpolation and dual reciprocity technique. Engineering Analysis with Boundary Elements, 142, 39–51. https://doi.org/10.1016/j.enganabound.2022.05.003
  » https://doi.org/10.1016/j.enganabound.2022.05.003
• Prodonoff, V., Zepka, S. (1983).Análise modal de sistemas com matrizes não simétricas. In: Procedings of the VII Congresso Brasileiro de Engenharia Mecânica, In portuguese, pag. . [S.l.: s.n.].
• Ramšak, M., Škerget, L. (2007). 3D Multi-domain BEM for solving the Laplace equation. Engineering Analysis with Boundary Elements, 31(6), 528-538. https://doi.org/10.1016/j.enganabound.2006.10.006
  » https://doi.org/10.1016/j.enganabound.2006.10.006
• Rokhlin, V. (1985). Rapid solution of integral equations of classical potential theory. Journal of Computational Physics, 60 (2), 187-207. https://doi.org/10.1016/0021-9991(85)90002-6
  » https://doi.org/10.1016/0021-9991(85)90002-6
• Shiara, L.S., Paschoalini, A.T. (2023). A Detailed Implementation of Multithreading and Out-of-core Computation to the Convention Boundary Element Algorithm with Minimum Code Changes. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 45, 114. https://doi.org/10.1007/s40430-023-04034-y
  » https://doi.org/10.1007/s40430-023-04034-y
• Wen, P.H., Aliabadi, M.H., Rooke, D.P. (1998). A New Method for Transformation of Domain Integrals to Boundary Integrals in Boundary Element Method. Communication in Numercial Method in Engineeering, 14, 1055-1065. https://doi.org/10.1002/(SICI)1099-0887(199811)14:11<1055::AID-CNM209>3.0.CO;2-6
  » https://doi.org/ https://doi.org/10.1002/(SICI)1099-0887(199811)14:11<1055::AID-CNM209>3.0.CO;2-6
• Wrobel, L.C., Aliabadi, M.H. (2002). The boundary element method. Chichester, UK: Wiley.
• Zhang, Q., Lallement, G., (1987). Comparison of normal eigenmodes calculation methods based on identified complex eigenmodes. Journal of Spacecraft and Rockets, 24 (1), 69–73.
• Zhang, Q., Lallement, G., Fillod, R. (1988). Relations between the right and left eigenvectors of non-symmetric structural models applications to rotors. Mechanical Systems and Signal Processing 2 (1), 97–103.

Edited by

Editor: Marco L. Bittencourt

Publication Dates

• Publication in this collection
  02 Feb 2024
• Date of issue
  2024

History

• Received
  03 Oct 2023
• Reviewed
  29 Nov 2023
• Accepted
  04 Dec 2023
• Published
  04 Dec 2023