Surface and acoustic waves natural frequencies estimation by boundary element method

Grzegorz Waśniewski

doi:10.2478/sgem-2026-0004

1

Introduction

Fluid–structure interaction problems are often encountered in modern engineering practice. In particular, those must be taken into account in an analysis of tanks, pipelines, tunnels, dams, submarines, offshore objects, etc. The listed constructions were investigated by many authors widely, and a reader can find examples of considerations in cited papers. From the point of view of mechanics, the phenomenon under consideration is an evolutionary process, for which it is important to determine the dynamic characteristics of the system, such as eigenfrequencies and associated eigenforms. They allow for the risk assessment of the resonance or for solving equations of motion using the modal superposition method.

The boundary element method (BEM) [1] is very useful in numerical solutions of physical problems described by the Laplace, Poisson, or Helmholtz equations. Its advantage, in relation to the more commonly used finite element method (FEM) [2], is the smaller number of unknowns, which is particularly important in the case of problems, which analyze requires determination of the eigenvalues. In particular, significant benefits can be seen when an exterior wave problem is considered. Also, the well-known disadvantages of the boundary element method should be mentioned. In contrast to the finite element method, one obtains full-populated and nonsymmetrical coefficient matrices, which increases the computational burden.

In the 1970s, when the finite element method was being developed intensively, it seemed to be the natural choice for a fluid domain discretization. In the previous article [3] and the references therein, various methodologies can be found for the description of compressible and incompressible fluids and their interaction with an elastic structure using this method.

Currently, many authors use a hybrid approach: the boundary element method for the fluid physical model and the finite element method for the structure [4,5,6]: the first two articles focus on filled or submerged structure natural frequencies and restrict considerations to incompressible fluids and the third one considered acoustic waves. Ye et al. [7] estimated surface waves and natural frequencies using a specific discretization technique. Zhao et al. [8] analyzed the coupling between acoustic waves and submerged fluid-filled cylindrical shell, but the structure frequencies are only estimated. In addition, it is worth pointing out papers, in which the natural frequencies of surface [9] and acoustic [10] waves were determined experimentally. Despite the fact that the authors of the aforementioned articles presented analyses of various issues, it is possible to indicate specific difficulties, whose overcoming have been proposed in this contribution.

The purpose of this article is to present consistent numerical models for acoustic and surface waves natural frequencies estimation, which occur in the compressible and incompressible fluids, respectively. Although the frequency spectra for both wave types lay in separate ranges and their eigenfrequencies can be estimated independently, which is confirmed by the references enclosed, the discussed algorithm allows to include a coupling between them. Cases are considered when the fluid fills a perfectly rigid or flexible tank. To obtain algebraic equation systems, the boundary element method and the finite element method are employed for the fluid and the elastic structure domains, respectively.

The compressible fluid is described by the Helmholtz differential equation, for which the fundamental solution is complex and depends on the number of waves. This unknown quantity remains under integral operators in the equivalent boundary integral formulation. Disadvantages caused by it, described in Section 2.1.1, disappears when the fundamental solution of the Laplace equation is applied instead of the Helmholtz equation. In this case, the fluid domain integral appears, which could be transformed into a sum of boundary integrals using the multiple reciprocity method. This approach, in which higher order fundamental solutions of the Laplace equation are employed, makes it possible to exclude the wavenumber from the integrals. The result is a polynomial expression with real coefficient matrices, so that the natural frequencies could be determined more efficiently.

The incompressible fluid is described by the Laplace equation. It appears in the multiple reciprocity method as a particular case of the compressible fluid, for which the polynomial expansion includes only zero degree terms.

The eigenproblem solution for the compressible fluid does not finds proper eigenvalues alone. This means that the original equation is necessary but not a sufficient condition for estimating proper natural frequencies. The refined eigenproblem formulation, which allows one to identify the proper values, is presented in the article.

This article is a continuation of the author’s research related to the boundary element method applications in the continuous media mechanics and surface girders analysis [11,12,13].

2

Fluid in rigid tank

2.1

Compressible fluid

2.1.1

Problem formulation

A homogenous and inviscid fluid, characterized by a sound speed $a$ a and a mass density $ρ$ \rho , in the domain $Ω$ \Omega with the boundary $Γ = Γ_{1} \cup Γ_{2} \cup Γ_{3}$ \Gamma ={\Gamma }_{1}\cup {\Gamma }_{2}\cup {\Gamma }_{3} divided on separate fluid free surface $(Γ_{F})$ ({\Gamma }_{{\rm{F}}}) and rigid $(Γ_{R})$ ({\Gamma }_{{\rm{R}}}) parts, respectively, is considered (Figure 1). In addition, if there is a small disturbance of the pressure and density fields, exclusion of the volume loads, the irrotational flow and harmonic oscillations, then the conservation principles presents the Helmholtz equation for the velocity potential amplitude [14]: (1) $\nabla^{2} ϕ + k^{2} ϕ = 0 in Ω,$ {\nabla }^{2}\phi +{k}^{2}\phi =0\hspace{0.25em}\text{in}\hspace{0.25em}\Omega , subject to the Robin or Neumann boundary conditions in the separate Γ parts: (2) $\frac{\partial ϕ}{\partial n} = \frac{ω^{2}}{g} ϕ on Γ_{F}, \frac{\partial ϕ}{\partial n} = 0 on Γ_{R},$ \frac{\partial \phi }{\partial n}=\frac{{\omega }^{2}}{g}\phi \hspace{0.25em}\text{on}\hspace{0.25em}{\Gamma }_{{\rm{F}}},\text{}\frac{\partial \phi }{\partial n}=0\hspace{0.25em}\text{on}\hspace{0.25em}{\Gamma }_{{\rm{R}}}, where $k = ω / a$ k=\omega /a is the wavenumber, $ω$ \omega denotes the circular frequency, $n$ n is the outward normal direction, and $g$ g is the gravitational acceleration. The fluid velocity field is defined by a potential function $v = \nabla ϕ$ {\boldsymbol{v}}=\nabla \phi . If the effect of free surface waves is neglected, the first boundary condition becomes Dirichlet’s condition (3) $ϕ = 0 o n Γ_{F}, \frac{\partial ϕ}{\partial n} = 0 on Γ_{R} .$ \phi =0\hspace{0.25em}\text{o}\text{n}\hspace{0.25em}{\Gamma }_{{\rm{F}}},\text{}\frac{\partial \phi }{\partial n}=0\hspace{0.25em}\text{on}\hspace{0.25em}{\Gamma }_{{\rm{R}}}.

The equivalent integral formulation to equation (1), derived from the Green’s second identity [15], is given as follows: (4) $c ϕ = \int_{Γ} \frac{\partial ϕ}{\partial n} g^{*} - ϕ \frac{\partial g^{*}}{\partial n}) d Γ,$ c\phi =\mathop{\int }\limits_{\Gamma }\left(\frac{\partial \phi }{\partial n}{g}^{* }-\phi \frac{\partial {g}^{* }}{\partial n}\right){\rm{d}}\Gamma , where $g^{*}$ {g}^{* } is the fundamental solution of the Helmholtz equation. Coefficient $c$ c depends on location of the collocation point. It takes values 1.0 or 0.0 in or out of the domain, respectively, and 0.5 on the smooth boundary. As it will be pointed out in Section 2.1.3, its value is not to be determined.

The fundamental solutions of the Helmholtz equation for two- and three-dimensional problems include unknown wavenumber k and present the following forms [16,17,18]: (5) $g^{*} = \begin{matrix} \frac{- i}{4} H_{0}^{(2)} (kr) - 2 D problem \\ \frac{e^{- ikr}}{4 π r} - 3 D problem \end{matrix},$ {g}^{* }=\left\{\begin{array}{c}\frac{-i}{4}{H}_{0}^{\left(2)}({kr})-2\text{D problem}\\ \frac{{{\rm{e}}}^{-{ikr}}}{4{\rm{\pi }}r}\hspace{0.25em}-3\text{D problem}\end{array}\right., where r denotes the distance from the integration to the collocation point, $i = \sqrt{- 1}$ i=\sqrt{-1} , and $H_{0}^{(2)}$ {H}_{0}^{\left(2)} is the Hankel function of the second kind. In both cases, the fundamental solutions are complex and include the unknown wavenumber. This unknown quantity remains under integral operators in the boundary integral formulation (equation (4)). The resulting boundary element method discrete equations become complex, nonlinear, and transcendental, inconvenient for the numerical solution. In order to determine the eigenvalues in this case, one should calculate the coefficient matrices step by step, for small increments of the wavenumber, and look for zero values of the determinant of a homogeneous algebraic equation system. The inconveniencies listed earlier are documented in the previous studies [6,8], where the responses of the system subjected to excitations with successive frequencies increased by a small increment were calculated. Resonant peaks visible on graphs of real or imaginary parts and modules of the determined complex variables allowed for identification of natural frequencies. Zhao et al. [8] found that the eigenvalue problem is transcendental and proposed an iterative process for frequency determination.

To overcome the mentioned disadvantages, the wavenumber-independent, fundamental solution of the Laplace equation $g_{0}^{*}$ {g}_{0}^{* } will be employed. In that case, the Green’s second identity leads to the following equation, in which the domain integral appears: (6) $c ϕ + \int_{Γ} ϕ \frac{\partial g_{0}^{*}}{\partial n} - \frac{\partial ϕ}{\partial n} g_{0}^{*}) d Γ = k^{2} \int_{Ω} ϕ g_{0}^{*} d Ω .$ c\phi +\mathop{\int }\limits_{\Gamma }\left(\phi \frac{\partial {g}_{0}^{* }}{\partial n}-\frac{\partial \phi }{\partial n}{g}_{0}^{* }\right){\rm{d}}\Gamma ={k}^{2}\mathop{\int }\limits_{\Omega }\phi {g}_{0}^{* }{\rm{d}}\Omega .

In the aforementioned equation, the unknown wavenumber is excluded from both integrals: boundary and domain.

2.1.2

Multiple reciprocity method

The multiple reciprocity method alters the domain integral to boundary one [19,20], using higher order fundamental solutions (7) $\nabla^{2} g_{i + 1}^{*} = g_{i}^{*}, i = 0, 1, \dots$ {\nabla }^{2}{g}_{i+1}^{* }={g}_{i}^{* },{i}=0,\hspace{0.25em}1,\hspace{0.25em}\ldots

Substituting equation (7) into the domain integral of equation (6) and taking into account the Green theorem, one obtains the following: (8) $\int_{Ω} ϕ g_{0}^{*} d Ω = \int_{Ω} ϕ \nabla^{2} g_{1}^{*} d Ω = \int_{Γ} ϕ \frac{\partial g_{1}^{*}}{\partial n} - \frac{\partial ϕ}{\partial n} g_{1}^{*}) d Ω - k^{2} \int_{Ω} ϕ g_{1}^{*} d Ω .$ \mathop{\int }\limits_{\Omega }\phi {g}_{0}^{* }{\rm{d}}\Omega =\mathop{\int }\limits_{\Omega }\phi {\nabla }^{2}{g}_{1}^{* }{\rm{d}}\Omega =\mathop{\int }\limits_{\Gamma }\left(\phi \frac{\partial {g}_{1}^{* }}{\partial n}-\frac{\partial \phi }{\partial n}{g}_{1}^{* }\right){\rm{d}}\Omega -{k}^{2}\mathop{\int }\limits_{\Omega }\phi {g}_{1}^{* }{\rm{d}}\Omega .

After N substitutions, equation (6) takes the following form: (9) $c ϕ + \sum_{i = 0}^{N} {(- k^{2})}^{i} \int_{Γ} ϕ \frac{\partial g_{i}^{*}}{\partial n} d Γ - \sum_{i = 0}^{N} {(- k^{2})}^{i} \int_{Γ} \frac{\partial ϕ}{\partial n} g_{i}^{*} d Γ = {(- 1)}^{N} k^{2 (N + 1)} \int_{Ω} ϕ g_{N}^{*} d Ω .$ c\phi +\mathop{\sum }\limits_{i=0}^{N}{(-{k}^{2})}^{i}\mathop{\int }\limits_{\Gamma }\phi \frac{\partial {g}_{i}^{* }}{\partial n}{\rm{d}}\Gamma -\mathop{\sum }\limits_{i=0}^{N}{(-{k}^{2})}^{i}\mathop{\int }\limits_{\Gamma }\frac{\partial \phi }{\partial n}{g}_{i}^{* }{\rm{d}}\Gamma ={(-1)}^{N}{k}^{2(N+1)}\mathop{\int }\limits_{\Omega }\phi {g}_{N}^{* }{\rm{d}}\Omega .

For sufficiently large N, the domain integral is negligible.

The boundary element method discrete form of equation (9) can be written as follows: (10) $[F_{0} - k^{2} F_{1} + \dots + {(- k^{2})}^{N} F_{N}] Φ = [G_{0} - k^{2} G_{1} + \dots + {(- k^{2})}^{N} G_{N}] Φ_{n} or F (k) Φ = G (k) Φ_{n},$ {[}{{\boldsymbol{F}}}_{0}-{k}^{2}{{\boldsymbol{F}}}_{1}+\ldots +{(-{k}^{2})}^{N}{{\boldsymbol{F}}}_{N}]{\boldsymbol{\Phi }}\hspace{1em}={[}{{\boldsymbol{G}}}_{0}-{k}^{2}{{\boldsymbol{G}}}_{1}+\ldots +{(-{k}^{2})}^{N}{{\boldsymbol{G}}}_{N}]{{\boldsymbol{\Phi }}}_{{\boldsymbol{n}}}\text{or}\hspace{0.25em}{\boldsymbol{F}}(k){\boldsymbol{\Phi }}={\boldsymbol{G}}(k){{\boldsymbol{\Phi }}}_{{\boldsymbol{n}}}, where $Φ$ {\boldsymbol{\Phi }} and $Φ_{n}$ {{\boldsymbol{\Phi }}}_{{\boldsymbol{n}}} are the nodal unknown vectors of velocity potential and its normal derivatives, respectively.

2.1.3

Fundamental solutions

A sequence of the fundamental solutions $g_{i}^{*}$ {g}_{i}^{* } and its normal derivatives $\partial g_{i}^{*} / \partial n$ \partial {g}_{i}^{* }/\partial n , which define the matrices $G_{i}$ {{\boldsymbol{G}}}_{i} and $F_{i}$ {{\boldsymbol{F}}}_{i} in the equation (10), is defined by recurrence – equation (7), starting from the Laplace equation fundamental solution: (11) $\nabla^{2} g_{0}^{*} = δ,$ {\nabla }^{2}{g}_{0}^{* }=\delta , where symbol $δ$ \delta denotes the Dirac delta.

The fundamental solution of the Laplace equation of any order for the two dimensional problem can be written in the following form [21]: (12) $g_{i}^{*} = \frac{r^{2 i}}{2 π} C_{i} \ln \frac{r}{R} - D_{i}), R > 0, i = 0, 1, \dots,$ {g}_{i}^{* }=\frac{{r}^{2i}}{2{\rm{\pi }}}\left({C}_{i}\mathrm{ln}\frac{r}{R}-{D}_{i}\right),\hspace{0.25em}R\gt 0,\text{}i=0,\hspace{0.25em}1,\hspace{0.25em}\ldots , where r denotes the distance from the integration to the collocation point and coefficients $C_{i}$ {C}_{i} and $D_{i}$ {D}_{i} are defined by recurrence: (13) $C_{0} = 1, D_{0} = 0, C_{i} = \frac{C_{i - 1}}{4 i^{2}}, D_{i} = \frac{1}{4 i^{2}} \frac{C_{i - 1}}{i} + D_{i - 1}), i > 0 .$ {C}_{0}=1,\text{}{D}_{0}=0,{C}_{i}=\frac{{C}_{i-1}}{4{i}^{2}},\text{}{D}_{i}=\frac{1}{4{i}^{2}}\left(\frac{{C}_{i-1}}{i}+{D}_{i-1}\right),\hspace{0.25em}i\gt 0.

Fundamental solutions for normal derivatives of any order are derived from equation (12) and are as follows: (14) $\frac{\partial g_{i}^{*}}{\partial n} = \frac{r^{2 i - 2} (n_{x} r_{x} + n_{y} r_{y})}{- 2 π} 2 i C_{i} \ln \frac{r}{R} - D_{i}) + C_{i}],$ \frac{\partial {g}_{i}^{* }}{\partial n}=\frac{{r}^{2i-2}({n}_{x}{r}_{x}+{n}_{y}{r}_{y})}{-2{\rm{\pi }}}\left[2i\left({C}_{i}\mathrm{ln}\frac{r}{R}-{D}_{i}\right)+{C}_{i}\right], where $r_{x}, r_{y}$ {r}_{x},\text{}{r}_{y} and $n_{x}, n_{y}$ {n}_{x},\text{}{n}_{y} are $r$ r distance and unitary normal outward vector components, respectively.

For $i = 0$ i=0 , the fundamental solution takes the following form: (15) $\frac{\partial g_{0}^{*}}{\partial n} = \frac{- (n_{x} r_{x} + n_{y} r_{y})}{2 π r^{2}} .$ \frac{\partial {g}_{0}^{* }}{\partial n}=\frac{-({n}_{x}{r}_{x}+{n}_{y}{r}_{y})}{2{\rm{\pi }}{r}^{2}}.

For the collocation point placed on the boundary, in particular, the parameter $r = 0$ r=0 . In consequence, the diagonal coefficients of the $F_{0}$ {{\boldsymbol{F}}}_{0} matrix must not to be derived by direct integration. They are obtained in a different way.

For a constant potential field, the derivate $\partial ϕ / \partial n = 0$ \partial \phi /\partial n=0 holds on the whole $Γ$ \Gamma . Following that, equation (10) for $N = 0$ N=0 takes the following form: (16) $F_{0} 1 = 0,$ {{\boldsymbol{F}}}_{0}{\bf{1}}{\boldsymbol{=}}{\bf{0}}, where $1$ {\bf{1}} is a vector of all components equal to unity and $0$ {\bf{0}} denotes the null vector. According to equation (16), the diagonal coefficients of the $F_{0}$ {{\boldsymbol{F}}}_{0} matrix can be summed from others: (17) $\underset{n = 1, \dots, dof}{⋀} F_{0 nn} = - \sum_{\begin{matrix} m = 1 \\ m \neq n \end{matrix}}^{do f} F_{0 nm} .$ \mathop{\bigwedge }\limits_{n=1,\hspace{0.25em}\ldots ,\hspace{0.25em}{\rm{dof}}}{F}_{0{\rm{nn}}}=-\mathop{\sum }\limits_{\begin{array}{c}{m}=\hspace{0.25em}1\\ {m}\ne {n}\end{array}}^{{\rm{do}}f}{F}_{0{\rm{nm}}}.

In accordance to equation (4), the diagonal components of the $F_{0}$ {{\boldsymbol{F}}}_{0} matrix contain $c$ c coefficients. The application of Equation (17) ensures proper values for $F_{0 nn}$ {F}_{0{\rm{nn}}} without $c$ c coefficient determination. Other fundamental solutions and their derivatives are integrable functions.

In the case of three-dimensional problems, fundamental solutions of Laplace equations of any order take a simpler form than this defined by equation (12) and are expressed as follows [22]: (18) $g_{i}^{*} = \frac{r^{2 i - 1}}{4 π (2 i)!}, i = 0, 1, \dots$ {g}_{i}^{* }=\frac{{r}^{2i-1}}{4\pi (2i)\!},\text{}i=0,\hspace{0.25em}1,\hspace{0.25em}\ldots

2.1.4

Natural frequencies estimation

Taking into account the boundary conditions in the form of equation (2), equation (10) takes the following form: (19) $\begin{matrix} F_{Γ_{F}} & F_{Γ_{R}} \end{matrix}] \begin{matrix} Φ_{Γ_{F}} \\ Φ_{Γ_{R}} \end{matrix}\} = \begin{matrix} G_{Γ_{F}} & G_{Γ_{R}} \end{matrix}] \begin{matrix} \frac{ω^{2}}{g} Φ_{Γ_{F}} \\ 0 \end{matrix}\} \Rightarrow \begin{matrix} F_{Γ_{F}} (k) - \frac{ω^{2}}{g} G_{Γ_{F}} (k) & F_{Γ_{R}} \end{matrix} (k)] Φ = 0 .$ \left[\begin{array}{cc}{{\boldsymbol{F}}}_{{\Gamma }_{F}}& {{\boldsymbol{F}}}_{{\Gamma }_{R}}\end{array}\right]\left\{\begin{array}{c}{{\boldsymbol{\Phi }}}_{{\Gamma }_{F}}\\ {{\boldsymbol{\Phi }}}_{{\Gamma }_{R}}\end{array}\right\}=\left[\begin{array}{cc}{{\boldsymbol{G}}}_{{\Gamma }_{F}}& {{\bf{G}}}_{{\Gamma }_{R}}\end{array}\right]\left\{\begin{array}{c}\frac{{\omega }^{2}}{g}{{\boldsymbol{\Phi }}}_{{\Gamma }_{F}}\\ {\bf{0}}\end{array}\right\}\Rightarrow \left[\begin{array}{cc}{{\boldsymbol{F}}}_{{\Gamma }_{F}}(k)-\frac{{\omega }^{2}}{g}{{\boldsymbol{G}}}_{{\Gamma }_{F}}(k)& {{\boldsymbol{F}}}_{{\Gamma }_{R}}\end{array}(k)\right]{\boldsymbol{\Phi }}{\boldsymbol{=}}{\bf{0}}{\boldsymbol{.}}

The existence of a nontrivial solution of equation (19) gives the expression for natural frequencies estimation: (20) $\det \begin{matrix} F_{Γ_{F}} \frac{ω}{a}) - \frac{ω^{2}}{g} G_{Γ_{F}} \frac{ω}{a}) & F_{Γ_{R}} \end{matrix} \frac{ω}{a})] = 0 \to ω_{j} .$ \det \left[\begin{array}{cc}{{\boldsymbol{F}}}_{{\Gamma }_{F}}\left(\frac{\omega }{a}\right)-\frac{{\omega }^{2}}{g}{{\boldsymbol{G}}}_{{\Gamma }_{F}}\left(\frac{\omega }{a}\right)& {{\boldsymbol{F}}}_{{\Gamma }_{R}}\end{array}\left(\frac{\omega }{a}\right)\right]{\boldsymbol{=}}0\to {\omega }_{j}.

Equation (20) estimates the natural frequencies for both surface and acoustic waves and takes into account the coupling between them. For $N = 0$ N=0 , only free surface waves frequencies are obtained. For higher N values, acoustic wave frequencies appear.

For alternative boundary conditions in the form of equation (3), when surface waves are neglected, one obtains the following form: (21) $\begin{matrix} F_{Γ_{F}} & F_{Γ_{R}} \end{matrix}] \begin{matrix} 0 \\ Φ_{Γ_{R}} \end{matrix}\} = \begin{matrix} G_{Γ_{F}} & G_{Γ_{R}} \end{matrix}] \begin{matrix} Φ_{n Γ_{F}} \\ 0 \end{matrix}\} \Rightarrow \begin{matrix} - G_{Γ_{F}} (k) & F_{Γ_{R}} \end{matrix} (k)] \begin{matrix} Φ_{n Γ_{F}} \\ Φ_{Γ_{R}} \end{matrix}\} = 0,$ \left[\begin{array}{cc}{{\boldsymbol{F}}}_{{\Gamma }_{F}}& {{\boldsymbol{F}}}_{{\Gamma }_{R}}\end{array}\right]\left\{\begin{array}{c}{\bf{0}}\\ {{\boldsymbol{\Phi }}}_{{\Gamma }_{R}}\end{array}\right\}=\left[\begin{array}{cc}{{\boldsymbol{G}}}_{{\Gamma }_{F}}& {{\boldsymbol{G}}}_{{\Gamma }_{R}}\end{array}\right]\left\{\begin{array}{c}{{\boldsymbol{\Phi }}}_{{\boldsymbol{n}}{\Gamma }_{F}}\\ {\bf{0}}\end{array}\right\}\Rightarrow \left[\begin{array}{cc}-{{\boldsymbol{G}}}_{{\Gamma }_{F}}(k)& {{\boldsymbol{F}}}_{{\Gamma }_{R}}\end{array}(k)\right]\left\{\begin{array}{c}{{\boldsymbol{\Phi }}}_{{\boldsymbol{n}}{\Gamma }_{F}}\\ {{\boldsymbol{\Phi }}}_{{\Gamma }_{R}}\end{array}\right\}{\boldsymbol{=}}{\bf{0}}, so that the natural frequencies are estimated from the following equation: (22) $\det \begin{matrix} - G_{Γ_{F}} \frac{ω}{a}) & F_{Γ_{R}} \end{matrix} \frac{ω}{a})] = 0 \to ω_{j} .$ \det \left[\begin{array}{cc}-{{\boldsymbol{G}}}_{{\Gamma }_{F}}\left(\frac{\omega }{a}\right)& {{\boldsymbol{F}}}_{{\Gamma }_{R}}\end{array}\left(\frac{\omega }{a}\right)\right]{\boldsymbol{=}}0\to {\omega }_{j}.

In the aforementioned expressions, the surface waves are neglected, and only the acoustic wave frequencies are obtained from equation (22).

2.1.5

Proper eigenvalues identification

The solution of equation (20) or equation (22) gives some improper eigenvalues for acoustic waves. This means that those equations are only necessary but not sufficient conditions for proper eigenvalue estimation. The appearance of such fictitious natural frequencies, when a compressible fluid is discretized by boundary elements, is a common problem [6, 18].

The refined problem formulation, which estimates only the proper natural frequencies, was proposed by Kamiya and Andoh [18]. Suppose that for a given specific eigenvalue, an eigenvector x is to be determined. Because both equations (19) and (21) are homogenous, only the relative magnitudes of the eigenvector elements are obtained. Thus, with the $l th$ l{\rm{th}} of this vector thought to be unity and that of the right side null vector to be $η$ \eta , for equations (19) and (21), one obtains (23) $A (ω) \begin{matrix} \begin{matrix} x_{1} \\ ⋮ \\ x_{l - 1} \end{matrix} \\ 1 \\ \begin{matrix} x_{l + 1} \\ ⋮ \\ x_{dof} \end{matrix} \end{matrix}\} = B (ω) \begin{matrix} \begin{matrix} 0 \\ ⋮ \\ 0 \end{matrix} \\ η \\ \begin{matrix} 0 \\ ⋮ \\ 0 \end{matrix} \end{matrix}\} .$ {\boldsymbol{A}}(\omega )\left\{\begin{array}{c}\begin{array}{c}{x}_{1}\\ \vdots \\ {x}_{l-1}\end{array}\\ 1\\ \begin{array}{c}{x}_{l+1}\\ \vdots \\ {x}_{{dof}}\end{array}\end{array}\right\}{\boldsymbol{=}}{\boldsymbol{B}}(\omega )\left\{\begin{array}{c}\begin{array}{c}0\\ \vdots \\ 0\end{array}\\ \eta \\ \begin{array}{c}0\\ \vdots \\ 0\end{array}\end{array}\right\}.

Rearranging the terms yields $(d \equiv dof)$ (d\equiv {\rm{dof}}) (24) $[a_{1} \dots - b_{l} \dots] \begin{matrix} \begin{matrix} x_{1} \\ ⋮ \end{matrix} \\ η \\ \begin{matrix} ⋮ \\ x_{d} \end{matrix} \end{matrix}\} = [b_{1} \dots - a_{l} \dots] \begin{matrix} \begin{matrix} 0 \\ ⋮ \end{matrix} \\ 1 \\ \begin{matrix} ⋮ \\ 0 \end{matrix} \end{matrix}\} = - a_{l} .$ {[}{{\boldsymbol{a}}}_{1}{\boldsymbol{\ldots }}{\boldsymbol{-}}{{\boldsymbol{b}}}_{l}{\boldsymbol{\ldots }}]\left\{\begin{array}{c}\begin{array}{c}{x}_{1}\\ \vdots \end{array}\\ \eta \\ \begin{array}{c}\vdots \\ {x}_{d}\end{array}\end{array}\right\}{\boldsymbol{=}}{[}{{\boldsymbol{b}}}_{1}{\boldsymbol{\ldots }}{\boldsymbol{-}}{{\boldsymbol{a}}}_{l}{\boldsymbol{\ldots }}]\left\{\begin{array}{c}\begin{array}{c}0\\ \vdots \end{array}\\ 1\\ \begin{array}{c}\vdots \\ 0\end{array}\end{array}\right\}{\boldsymbol{=}}{\boldsymbol{-}}{{\boldsymbol{a}}}_{l}.

Now $η$ \eta is determined by the Cramer formula: (25) $η = \frac{\det A}{\det \begin{matrix} a_{1} \dots a_{l - 1} & - b_{l} & a_{l + 1} \dots a_{dof} \end{matrix}]} .$ \eta =\frac{\det {\boldsymbol{A}}}{\det \left[\begin{array}{ccc}{{\boldsymbol{a}}}_{1}{\boldsymbol{\ldots }}{{\boldsymbol{a}}}_{l-1}& {\boldsymbol{-}}{{\boldsymbol{b}}}_{l}& {{\boldsymbol{a}}}_{l+1}{\boldsymbol{\ldots }}{{\boldsymbol{a}}}_{{\rm{dof}}}\end{array}\right]}.

The condition $η = 0$ \eta =0 holds, so that the following equation is possible (26) $\frac{\det A}{\det \begin{matrix} a_{1} \dots a_{l - 1} & - b_{l} & a_{l + 1} \dots a_{dof} \end{matrix}]} = 0 \to ω_{j} .$ \frac{\det {\boldsymbol{A}}}{\det \left[\begin{array}{ccc}{{\boldsymbol{a}}}_{1}{\boldsymbol{\ldots }}{{\boldsymbol{a}}}_{l-1}& {\boldsymbol{-}}{{\boldsymbol{b}}}_{l}& {{\boldsymbol{a}}}_{l+1}{\boldsymbol{\ldots }}{{\boldsymbol{a}}}_{{\rm{dof}}}\end{array}\right]}=0\to {\omega }_{j}.

It can be identified as the original eigenproblem for equations (19) and (21), which simultaneously satisfies equations (20) and (22). Equation (26) estimates only the proper eigenvalues.

An alternative approach to identifying improper values is worth mentioning. When eigenmodes are determined from the original eigenproblems defined by equation (19) or equation (21), their shape can help to distinguish between the proper and fictitious frequencies associated with them. However, this approach is subjective and requires the experience of the evaluator, while the procedure proposed earlier is an objective analytical method.

2.2

Incompressible fluid

Following the assumptions specified earlier, the continuity equation for the incompressible fluid becomes the Laplace one [14]: (27) $\nabla^{2} ϕ = 0 in Ω$ {\nabla }^{2}\phi =0\hspace{0.25em}\text{in}\hspace{0.25em}\Omega with boundary conditions given by equation (2). The integral formulation of equation (27) takes the following form: (28) $c ϕ = \int_{Γ} \frac{\partial ϕ}{\partial n} g_{0}^{*} - ϕ \frac{\partial g_{0}^{*}}{\partial n}) d Γ .$ c\phi =\mathop{\int }\limits_{\Gamma }\left(\frac{\partial \phi }{\partial n}{g}_{0}^{* }-\phi \frac{\partial {g}_{0}^{* }}{\partial n}\right){\rm{d}}\Gamma .

The discrete form of the boundary element method of equation (28) is equivalent to equation (10) for $N = 0 :$ N=0: (29) $F_{0} Φ = G_{0} Φ_{n},$ {{\boldsymbol{F}}}_{0}{\boldsymbol{\Phi }}={{\boldsymbol{G}}}_{0}{{\boldsymbol{\Phi }}}_{{\boldsymbol{n}}},

and natural frequencies could be calculated from equation (20) after omitting the higher-order fundamental solutions: (30) $\det \begin{matrix} F_{0 Γ_{F}} - \frac{ω^{2}}{g} G_{0 Γ_{F}} & F_{0 Γ_{R}} \end{matrix}] = 0 \to ω_{j} .$ \det \left[\begin{array}{cc}{{\boldsymbol{F}}}_{0{\Gamma }_{F}}-\frac{{\omega }^{2}}{g}{{\boldsymbol{G}}}_{0{\Gamma }_{F}}& {{\boldsymbol{F}}}_{0{\Gamma }_{R}}\end{array}\right]{\boldsymbol{=}}0\to {\omega }_{j}.

For the incompressible fluid, there exist only free surface waves, so that, for the boundary conditions in the form of equation (3), for $N = 0$ N=0 , equation (22) is frequency independent.

2.3

Numerical examples

2.3.1

General remarks

In all the examples included in the article, a fluid that occupies a rectangular two-dimensional domain is considered. In the BEM discretization, constant elements with central node were used. Figure 2 shows the geometry of the problem. The parameter $R$ R in the fundamental solutions was assumed to be equal to $1.0 m$ 1.0\text{m} .

Example 1

Model validation for surface waves

In the example, the convergence of the surface waves natural frequencies for incompressible fluid, due to mesh refinement, is analyzed. The gravitational acceleration $g$ g was assumed to be equal to $9.81 m / s^{2}$ 9.81\hspace{0.25em}\text{m}/{\text{s}}^{2} . The results for the lowest natural frequencies $f_{j} = ω_{j} / 2 π$ {f}_{j}={\omega }_{j}/2{\rm{\pi }} estimated on the base of equation (30) are shown in Table 1. A convenient way to identify the unknowns is to trace the determinant sign (Figure 3). In the considered case, the Laplace equation solution is expressed as follows: (31) $ϕ = (A \sin λ x + B \cos λ x) (C e^{- λ y} + D e^{λ y}) .$ \phi =(A\sin \lambda x+B\cos \lambda x)(C{{\rm{e}}}^{-\lambda y}+D{{\rm{e}}}^{\lambda y}).

Table 1

The lowest natural frequencies for surface waves (Hz).

No of dof	$f_{1}$ {f}_{1}	$f_{2}$ {f}_{2}	$f_{3}$ {f}_{3}	$f_{4}$ {f}_{4}	$f_{5}$ {f}_{5}	$f_{6}$ {f}_{6}	$f_{7}$ {f}_{7}
36	0.8496	1.226	1.532	1.782	2.009	2.202	2.352
60	0.8438	1.210	1.493	1.740	1.962	2.167	2.357
84	0.8417	1.203	1.481	1.721	1.935	2.133	2.319
Exact	0.8375	1.191	1.459	1.685	1.884	2.064	2.229

Source: Author’s contribution.

Taking into account the boundary conditions defined in equation (2), it is possible to determine constants and exact values of natural frequencies: $A = 0, C = D, λ = j π / L_{x}$ A=0,\hspace{0.25em}C=D,\hspace{0.25em}\lambda =j{\rm{\pi }}/{L}_{x} : (32) $f_{j} = \sqrt{\frac{jg}{4 π L_{x}} \frac{1 - e^{\frac{- 2 j π L_{y}}{L_{x}}}}{1 + e^{\frac{- 2 j π L_{y}}{L_{x}}}}}, j = 0, 1, 2, \dots$ {f}_{j}=\sqrt{\frac{{jg}}{4{\rm{\pi }}{L}_{x}}\frac{1-{e}^{\frac{-2j{\rm{\pi }}{L}_{y}}{{L}_{x}}}}{1+{e}^{\frac{-2j{\rm{\pi }}{L}_{y}}{{L}_{x}}}}},\hspace{1em}j=0,1,2,\ldots

Although $j = 0$ j=0 is allowable, the results for it will not be presented.

The natural frequencies in the following examples will be determined for discrete models with 84 BEM degrees of freedom (Figure 2).

Example 2

Model validation for acoustic waves

The example is focused on the natural frequencies of acoustic waves in the absence of surface waves. The convergence of solutions is investigated when the number of fundamental solutions employed in the multiple reciprocity method increases. The sound speed, which characterizes acoustic fluid, was assumed to be $a = 1,450 m / s$ a=\mathrm{1,450}\text{m}/\text{s} . The natural frequencies, including the improper ones, obtained from equation (22) are presented in Table 2. The proper frequencies are shown in bold and are accompanied by an analytical solution in the same column. In the considered case, the solution of the Helmholtz equation is expressed as follows: (33) $ϕ = (A \cos μ_{1} x + B \sin μ_{1} x) (C \cos μ_{2} y + D \sin μ_{2} y) .$ \phi =(A\cos {\mu }_{1}{\rm{x}}+B\sin {\mu }_{1}x)(C\cos {{\rm{\mu }}}_{2}y+D\sin {\mu }_{2}y).

Table 2

The lowest natural frequencies for acoustic waves (Hz).

$N$ N	$f_{1}$ {f}_{1}	$f_{2}$ {f}_{2}	$f_{3}$ {f}_{3}	$f_{4}$ {f}_{4}	$f_{5}$ {f}_{5}	$f_{6}$ {f}_{6}
1		300.2
3	393.7	838.1	589.7		665.4
5	401.0	502.3	785.2	1,030	970.3	1,461
7	401.0	502.0	772.7	1,115	1,226	1,577
9	401.0	502.0	772.6	1,137	1,209	1,311
14	401.0	502.0	772.6	1,138	1,206	1,326
15	401.0	502.0	772.6	1,138	1,206	1,326
Exact	402.8		772.4		1,208

Source: Author’s contribution.

Taking into account the boundary conditions defined in equation (3), it is possible to determine constants and exact values of natural frequencies: $B = D = 0, k = \sqrt{μ_{1}^{2} + μ_{2}^{2}}$ B=D=0,\hspace{0.25em}k=\sqrt{{\mu }_{1}^{2}+{\mu }_{2}^{2}} , (34) $f_{j} = \frac{a}{2} \sqrt{{\frac{n}{L_{x}})}^{2} + {\frac{m + 0.5}{L_{y}})}^{2}}, m, n = 0, 1, 2, \dots$ {f}_{j}=\frac{a}{2}\sqrt{{\left(\frac{n}{{L}_{x}}\right)}^{2}+{\left(\frac{m+0.5}{{L}_{y}}\right)}^{2}},\hspace{1em}m,{n}=0,\hspace{1em}1,\hspace{1em}2,\hspace{1em}\ldots

The determinant sign obtained from equation (22), for the number of fundamental solutions $N = 9$ N=9 , is shown in Figure 4. To identify proper frequencies, the left side ratio from equation (26) is presented in Figures 5 and 6. As it turns out, the effectiveness of the procedure depends on the position of the $l th$ l{\rm{th}} component highlighted in equation (23) (Figure 2). When it is located in the center of the edge $Γ_{2}$ {\Gamma }_{2} (node – $l_{V}$ {l}_{{\rm{V}}} ), valid values can be identified (Figure 5). On the other hand, for the position in the center of $Γ_{3}$ {\Gamma }_{3} (node – $l_{H}$ {l}_{{\rm{H}}} ), the frequency $f_{3}$ {f}_{3} is lost (Figure 6).

For a few numbers of the fundamental solutions, the process converges for a few first natural frequencies. For the number of fundamental solutions $N \geq 9$ N\ge 9 , the proper values are close to the analytical solution. In the following examples, $N = 9$ \hspace{0.25em}N=9 is assumed to be enough in the acoustic fluid approximation.

Example 3

Surface and acoustic waves coupling

Equation (20) determines the natural frequencies with the simultaneous presence of the surface and the acoustic waves and takes into account the coupling between them. The sound speed occurring in it is the only physical parameter that distinguishes the fluid analyzed. In the example, the value of $a$ a is investigated for which, the coupling between both types of waves has a significant impact on the natural frequencies estimation.

The surface waves frequencies for different sound speed values are presented in Table 3. The results are compared with those derived for incompressible fluid in Example 1.

Table 3

Surface waves natural frequencies vs sound speed (Hz).

$a (m / s)$ {a}(\text{m}/\text{s})	$f_{1}$ {f}_{1}	$f_{2}$ {f}_{2}	$f_{3}$ {f}_{3}	$f_{4}$ {f}_{4}	$f_{5}$ {f}_{5}	$f_{6}$ {f}_{6}
5	0.8113	1.183	1.464	1.708	1.921	2.119
10	0.8340	1.198	1.477	1.717	1.932	2.130
30	0.8408	1.203	1.481	1.720	1.935	2.133
100	0.8416	1.203	1.481	1.721	1.935	2.133
500	0.8417	1.203	1.481	1.721	1.935	2.133
1,450	0.8417	1.203	1.481	1.721	1.935	2.133
Example 1	0.8417	1.203	1.481	1.721	1.935	2.133

Source: Author’s contribution.

Table 4 presents the natural frequencies of the acoustic waves, including the improper ones, estimated for the coupled problem (Equation (20)). The results are compared with values derived under the condition of surface waves absence (Equation (22)). The proper frequencies are marked by bold.

Table 4

Acoustic waves natural frequencies vs sound speed (Hz).

$a (m / s)$ {a}(\text{m}/\text{s})	Case	$f_{1}$ {f}_{1}	$f_{2}$ {f}_{2}	$f_{3}$ {f}_{3}	$f_{4}$ {f}_{4}	$f_{5}$ {f}_{5}	$f_{6}$ {f}_{6}
5	Equation (20)	1.553	1.749	2.693	3.921	4.238	4.515
5	Equation (22)	1.383	1.731	2.664	3.922	4.170	4.520
10	Equation (20)	2.859	3.464	5.342	7.844	8.374	9.038
10	Equation (22)	2.766	3.462	5.328	7.844	8.340	9.040
30	Equation (20)	8.330	10.39	15.99	23.53	25.03	27.12
30	Equation (22)	8.297	10.39	15.98	23.53	25.02	27.12
100	Equation (20)	27.67	34.62	53.28	78.44	83.40	90.40
100	Equation (22)	27.66	34.62	53.28	78.44	83.40	90.40
500	Equation (20)	138.3	173.1	266.4	392.2	417.0	452.0
500	Equation (22)	138.3	173.1	266.4	392.2	417.0	452.0
1,450	Equation (20)	401.0	502.0	772.6	1,137	1,209	1,311
1,450	Equation (22)	401.0	502.0	772.6	1,137	1,209	1,311

Source: Author’s contribution.

Sound speed for real liquids does not exceed $1,000 m / s$ \mathrm{1,000}\hspace{0.25em}\text{m}/\text{s} , so for considered problem geometry, a coupling effect is not important. Moreover, for a hypothetical fluid, characterized by extremely low sound speed values, even when frequency spectra overlap, results modification occurs to be small.

3

Fluid in flexible tank

3.1

Fluid structure interaction

One divides the tank walls on separate wet $(Γ_{W})$ ({\Gamma }_{{\rm{W}}}) and dry $(Γ_{D})$ ({\Gamma }_{{\rm{D}}}) parts (Figure 1). On the fluid–structure interface $(Γ_{W})$ ({\Gamma }_{{\rm{W}}}) , kinetic and kinematic coupling must be considered.

The kinetic coupling will be based on the Navier–Stokes equation [14], which under assumptions specified in Section 2.1.1 takes a simple form [3]: (35) $\nabla p = - ρ_{f} \dot{v},$ \nabla p=-{\rho }_{{\rm{f}}}\dot{{\boldsymbol{v}}}, which joins fluid acceleration $\dot{v}$ \dot{{\boldsymbol{v}}} with hydrodynamic pressure $p$ p and fluid mass density $ρ_{f}$ {\rho }_{{\rm{f}}} . Differentiating equation (35) with respect to time and taking into account that $v = \nabla ϕ$ {\boldsymbol{v}}=\nabla \phi and BEM approximation for fluid velocity potential function $ϕ = N^{B} ϕ$ \phi ={{\boldsymbol{N}}}^{{\rm{B}}}{\boldsymbol{\phi }} , one obtains expression for pressure velocity on the fluid–structure interface, for the considered problem of natural vibrations: (36) $\nabla \dot{p} = - ρ_{f} \ddot{v} = - ρ_{f} \nabla \ddot{ϕ} ⟹ {\dot{p}}_{Γ_{W}} = ω^{2} ρ_{f} N_{Γ_{W}}^{B} ϕ_{Γ_{W}} .$ \nabla \dot{p}=-{\rho }_{{\rm{f}}}\ddot{{\boldsymbol{v}}}=-{\rho }_{{\rm{f}}}\nabla \ddot{\phi }\Longrightarrow {\dot{p}}_{{\Gamma }_{{\rm{W}}}}{\boldsymbol{=}}{\omega }^{2}{\rho }_{{\rm{f}}}{{\boldsymbol{N}}}_{{\Gamma }_{{\rm{W}}}}^{{\rm{B}}}{{\boldsymbol{\phi }}}_{{\Gamma }_{{\rm{W}}}}.

For a flexible tank, in addition to stiffness, its inertia must also be taken into account. In a natural vibration problem, external loads do not appear, so only hydrodynamic pressure acts on the structure part $Γ_{W}$ {\Gamma }_{{\rm{W}}} . As mentioned earlier, the equilibrium equation of the finite element method, written for the assumed global vector of nodal parameters $q$ {\boldsymbol{q}} , takes the following form: (37) $\begin{matrix} K_{Γ_{W}} & K_{Γ_{D}} \end{matrix}] - ω^{2} \begin{matrix} B_{Γ_{W}} & B_{Γ_{D}} \end{matrix}] \begin{matrix} q_{Γ_{W}} \\ q_{Γ_{D}} \end{matrix}\} = \begin{matrix} f_{Γ_{W}} (p) \\ 0 \end{matrix}\},$ \left[\begin{array}{cc}{{\boldsymbol{K}}}_{{\Gamma }_{{\rm{W}}}}& {{\boldsymbol{K}}}_{{\Gamma }_{{\rm{D}}}}\end{array}\right]-{\omega }^{2}\left[\begin{array}{cc}{{\boldsymbol{B}}}_{{\Gamma }_{{\rm{W}}}}& {{\boldsymbol{B}}}_{{\Gamma }_{{\rm{D}}}}\end{array}\right]\left\{\begin{array}{c}{{\boldsymbol{q}}}_{{\Gamma }_{{\rm{W}}}}\\ {{\boldsymbol{q}}}_{{\Gamma }_{{\rm{D}}}}\end{array}\right\}{\boldsymbol{=}}\left\{\begin{array}{c}{{\boldsymbol{f}}}_{{\Gamma }_{{\rm{W}}}}(p)\\ {\boldsymbol{0}}\end{array}\right\}{\boldsymbol{,}} where the coefficient matrices do not depend on time and its parts, representing wet and dry nodes, will be denoted as follows: (38) $\begin{matrix} K_{Γ_{W}} & K_{Γ_{D}} \end{matrix}] - ω^{2} \begin{matrix} B_{Γ_{W}} & B_{Γ_{D}} \end{matrix}] = \begin{matrix} S_{Γ_{W}}^{Γ_{W}} (ω) & S_{Γ_{D}}^{Γ_{W}} (ω) \\ S_{Γ_{W}}^{Γ_{D}} (ω) & S_{Γ_{D}}^{Γ_{D}} (ω) \end{matrix}] .$ \left[\begin{array}{cc}{{\boldsymbol{K}}}_{{\Gamma }_{{\rm{W}}}}& {{\boldsymbol{K}}}_{{\Gamma }_{{\rm{D}}}}\end{array}\right]-{\omega }^{2}\left[\begin{array}{cc}{{\boldsymbol{B}}}_{{\Gamma }_{{\rm{W}}}}& {{\boldsymbol{B}}}_{{\Gamma }_{{\rm{D}}}}\end{array}\right]{\boldsymbol{=}}\left[\begin{array}{cc}{{\boldsymbol{S}}}_{{\Gamma }_{{\rm{W}}}}^{{\Gamma }_{{\rm{W}}}}(\omega )& {{\bf{S}}}_{{\Gamma }_{{\rm{D}}}}^{{\Gamma }_{{\rm{W}}}}(\omega )\\ {{\boldsymbol{S}}}_{{\Gamma }_{{\rm{W}}}}^{{\Gamma }_{{\rm{D}}}}(\omega )& {{\bf{S}}}_{{\Gamma }_{{\rm{D}}}}^{{\Gamma }_{{\rm{D}}}}(\omega )\end{array}\right]{\boldsymbol{.}}

When equation (37) is differentiated with respect to time, then the nonzero elements of the load vector can be determined from equation (36) and a resulting relation will be expressed by the BEM nodal unknowns vector $ϕ :$ {\boldsymbol{\phi }}: (39) $\begin{matrix} S_{Γ_{W}}^{Γ_{W}} & S_{Γ_{D}}^{Γ_{W}} \\ S_{Γ_{W}}^{Γ_{D}} & S_{Γ_{D}}^{Γ_{D}} \end{matrix}] \begin{matrix} {\dot{q}}_{Γ_{W}} \\ {\dot{q}}_{Γ_{D}} \end{matrix}\} = \begin{matrix} f_{Γ_{W}} (\dot{p}) \\ 0 \end{matrix}\} = \begin{matrix} T (ω) ϕ_{Γ_{W}} \\ 0 \end{matrix}\} .$ \left[\begin{array}{cc}{{\boldsymbol{S}}}_{{\Gamma }_{{\rm{W}}}}^{{\Gamma }_{{\rm{W}}}}& {{\boldsymbol{S}}}_{{\Gamma }_{{\rm{D}}}}^{{\Gamma }_{{\rm{W}}}}\\ {{\boldsymbol{S}}}_{{\Gamma }_{{\rm{W}}}}^{{\Gamma }_{{\rm{D}}}}& {{\boldsymbol{S}}}_{{\Gamma }_{{\rm{D}}}}^{{\Gamma }_{{\rm{D}}}}\end{array}\right]\left\{\begin{array}{c}{\dot{{\boldsymbol{q}}}}_{{\Gamma }_{{\rm{W}}}}\\ {\dot{{\boldsymbol{q}}}}_{{\Gamma }_{{\rm{D}}}}\end{array}\right\}{\boldsymbol{=}}\left\{\begin{array}{c}{{\bf{f}}}_{{\Gamma }_{{\rm{W}}}}(\dot{p})\\ {\bf{0}}\end{array}\right\}{\boldsymbol{=}}\left\{\begin{array}{c}{\boldsymbol{T}}(\omega ){{\boldsymbol{\phi }}}_{{\Gamma }_{{\rm{W}}}}\\ {\bf{0}}\end{array}\right\}{\boldsymbol{.}}

The $T (ω)$ {\boldsymbol{T}}(\omega ) matrix definition results from the FEM assembling procedure and the tank wall transverse displacement approximation $w = N^{F} q .$ w={{\boldsymbol{N}}}^{{\rm{F}}}{\boldsymbol{q}}{\boldsymbol{.}} (40) $T (ω) = - ω^{2} ρ_{f} \int_{Γ_{W}} {(N_{Γ_{W}}^{F})}^{T} N_{Γ_{W}}^{B} d Γ .$ {\boldsymbol{T}}(\omega ){\boldsymbol{=}}{\boldsymbol{-}}{\omega }^{2}{\rho }_{{\rm{f}}}\mathop{\int }\limits_{{\Gamma }_{{\rm{W}}}}{({{\boldsymbol{N}}}_{{\Gamma }_{{\rm{W}}}}^{{\rm{F}}})}^{{\rm{T}}}{{\boldsymbol{N}}}_{{\Gamma }_{{\rm{W}}}}^{{\rm{B}}}{\rm{d}}\Gamma {\boldsymbol{.}}

The negative sign in the aforementioned equation results from the orientations of the finite element local coordinate systems axes, which were assumed for the normal to the walls ones, were assumed inside the tank.

The kinematic coupling follows from the equality of fluid and structure velocities in the normal to wall direction and leads to the relation between the BEM and FEM nodal parameters: (41) $\frac{\partial ϕ}{\partial n} = - \dot{w} ⟹ Φ_{n Γ_{W}} = - N_{Γ_{W}}^{F} {\dot{q}}_{Γ_{W}} .$ \frac{\partial \phi }{\partial n}=-\dot{w}\Longrightarrow {{\boldsymbol{\Phi }}}_{{\boldsymbol{n}}{\Gamma }_{{\rm{W}}}}{\boldsymbol{=}}{\boldsymbol{-}}{{\boldsymbol{N}}}_{{\Gamma }_{{\rm{W}}}}^{{\rm{F}}}{\dot{{\boldsymbol{q}}}}_{{\Gamma }_{{\rm{W}}}}.

The reason for the negative sign is the same as shown in equation (40).

3.2

Acoustic waves in compressible fluid

The results presented in Example 3 recommend that, in further considerations, the surface waves can be analyzed for the incompressible fluid model, and for the acoustic waves examination, surfaces may be neglected.

By introducing the boundary part $Γ_{W}$ {\Gamma }_{{\rm{W}}} and substituting equation (41) into equation (21), one obtains (42) $\begin{matrix} - G_{Γ_{F}} & F_{Γ_{W}} & F_{Γ_{R}} \end{matrix}] \begin{matrix} Φ_{n Γ_{F}} \\ Φ_{Γ_{W}} \\ Φ_{Γ_{R}} \end{matrix}\} = - G_{Γ_{W}} N_{Γ_{W}}^{F} {\dot{q}}_{Γ_{W}} .$ \left[\begin{array}{ccc}-{{\boldsymbol{G}}}_{{\Gamma }_{{\rm{F}}}}& {{\boldsymbol{F}}}_{{\Gamma }_{{\rm{W}}}}& {{\boldsymbol{F}}}_{{\Gamma }_{{\rm{R}}}}\end{array}\right]\left\{\begin{array}{c}{{\boldsymbol{\Phi }}}_{{\boldsymbol{n}}{\Gamma }_{{\rm{F}}}}\\ {{\boldsymbol{\Phi }}}_{{\Gamma }_{{\rm{W}}}}\\ {{\boldsymbol{\Phi }}}_{{\Gamma }_{{\rm{R}}}}\end{array}\right\}{\boldsymbol{=}}{\boldsymbol{-}}{{\boldsymbol{G}}}_{{\Gamma }_{{\rm{W}}}}{{\boldsymbol{N}}}_{{\Gamma }_{{\rm{W}}}}^{{\rm{F}}}{\dot{{\boldsymbol{q}}}}_{{\Gamma }_{{\rm{W}}}}{\boldsymbol{.}}

Equations (39) and (42), written together, lead to homogeneous equations set: (43) $\begin{matrix} 0 & - T & 0 & S_{Γ_{W}}^{Γ_{W}} & S_{Γ_{D}}^{Γ_{W}} \\ 0 & 0 & 0 & S_{Γ_{W}}^{Γ_{D}} & S_{Γ_{D}}^{Γ_{D}} \\ - G_{Γ_{F}} & F_{Γ_{W}} & F_{Γ_{R}} & G_{Γ_{W}} N_{Γ_{W}}^{F} & 0 \end{matrix}] \begin{matrix} Φ_{n Γ_{F}} \\ \begin{matrix} Φ_{Γ_{W}} \\ Φ_{Γ_{R}} \\ {\dot{q}}_{Γ_{W}} \end{matrix} \\ {\dot{q}}_{Γ_{D}} \end{matrix}\} = 0,$ \left[\begin{array}{ccccc}{\bf{0}}& -{\boldsymbol{T}}& {\bf{0}}& {{\boldsymbol{S}}}_{{\Gamma }_{{\rm{W}}}}^{{\Gamma }_{{\rm{W}}}}& {{\boldsymbol{S}}}_{{\Gamma }_{{\rm{D}}}}^{{\Gamma }_{{\rm{W}}}}\\ {\bf{0}}& {\bf{0}}& {\bf{0}}& {{\boldsymbol{S}}}_{{\Gamma }_{{\rm{W}}}}^{{\Gamma }_{{\rm{D}}}}& {{\boldsymbol{S}}}_{{\Gamma }_{{\rm{D}}}}^{{\Gamma }_{{\rm{D}}}}\\ -{{\boldsymbol{G}}}_{{\Gamma }_{{\rm{F}}}}& {{\boldsymbol{F}}}_{{\Gamma }_{{\rm{W}}}}& {{\boldsymbol{F}}}_{{\Gamma }_{{\rm{R}}}}& {{\boldsymbol{G}}}_{{\Gamma }_{{\rm{W}}}}{{\bf{N}}}_{{\Gamma }_{{\rm{W}}}}^{{\rm{F}}}& {\bf{0}}\end{array}\right]\left\{\begin{array}{c}{{\boldsymbol{\Phi }}}_{{\boldsymbol{n}}{\Gamma }_{{\rm{F}}}}\\ \begin{array}{c}{{\boldsymbol{\Phi }}}_{{\Gamma }_{{\rm{W}}}}\\ {{\boldsymbol{\Phi }}}_{{\Gamma }_{{\rm{R}}}}\\ {\dot{{\boldsymbol{q}}}}_{{\Gamma }_{{\rm{W}}}}\end{array}\\ {\dot{{\boldsymbol{q}}}}_{{\Gamma }_{{\rm{D}}}}\end{array}\right\}{\boldsymbol{=}}{\bf{0}}{\boldsymbol{,}} which allows to determine acoustic waves natural frequencies.

3.3

Surface waves for incompressible fluid

A coefficients matrix of the homogeneous equations set for surface waves natural frequencies estimation follows from equations (19), (29), (39), and (41) and takes the following form: (44) $\begin{matrix} 0 & - T & 0 & S_{Γ_{W}}^{Γ_{W}} & S_{Γ_{D}}^{Γ_{W}} \\ 0 & 0 & 0 & S_{Γ_{W}}^{Γ_{D}} & S_{Γ_{D}}^{Γ_{D}} \\ F_{0 Γ_{F}} - \frac{ω^{2}}{g} G_{0 Γ_{F}} & F_{0 Γ_{W}} & F_{0 Γ_{R}} & G_{0 Γ_{W}} N_{Γ_{W}}^{F} & 0 \end{matrix}] .$ \left[\begin{array}{ccccc}{\bf{0}}& -{\boldsymbol{T}}& {\bf{0}}& {{\boldsymbol{S}}}_{{\Gamma }_{{\rm{W}}}}^{{\Gamma }_{{\rm{W}}}}& {{\boldsymbol{S}}}_{{\Gamma }_{{\rm{D}}}}^{{\Gamma }_{{\rm{W}}}}\\ {\bf{0}}& {\bf{0}}& {\bf{0}}& {{\boldsymbol{S}}}_{{\Gamma }_{{\rm{W}}}}^{{\Gamma }_{{\rm{D}}}}& {{\boldsymbol{S}}}_{{\Gamma }_{{\rm{D}}}}^{{\Gamma }_{{\rm{D}}}}\\ {{\boldsymbol{F}}}_{0{\Gamma }_{{\rm{F}}}}-\frac{{\omega }^{2}}{g}{{\boldsymbol{G}}}_{0{\Gamma }_{{\rm{F}}}}& {{\boldsymbol{F}}}_{0{\Gamma }_{{\rm{W}}}}& {{\boldsymbol{F}}}_{0{\Gamma }_{{\rm{R}}}}& {{\boldsymbol{G}}}_{0{\Gamma }_{{\rm{W}}}}{{\boldsymbol{N}}}_{{\Gamma }_{{\rm{W}}}}^{{\rm{F}}}& {\bf{0}}\end{array}\right].

3.4

Numerical examples

3.4.1

General remarks

For the considered two-dimensional problem, which results in the cylindrical bending of the tank wall only, Kirchhoff’s theory is assumed to be a physical model of the structure. The conforming finite elements with third-degree polynomials as base functions are applied [2]. The inertia matrix is determined as a consistent one.

The steel tank, characterized by material constants: Young modulus $E = 205 GPa$ E=205\hspace{0.25em}\text{GPa} , Poisson’s ratio $ν = 0.3$ {\rm{\nu }}=0.3 , mas density $ρ_{s} = 7,850 kg / m^{3}$ {\rho }_{{\rm{s}}}=\mathrm{7,850}\text{kg}/{\text{m}}^{3} , will be analyzed. For fluid one assumes: sound speed $a = 1,450 m / s$ a=\mathrm{1,450}\text{m}/\text{s} and mass density $ρ_{f} = 1,000 kg / m^{3}$ {\rho }_{{\rm{f}}}=\mathrm{1,000}\text{kg}/{\text{m}}^{3} . The considered tank geometry and discrete model mesh are displayed in Figure 2.

Example 4

Acoustic waves in compressible fluid

To estimate the natural frequencies of acoustic waves, the sign of the coefficient matrix determinant defined in equation (43) is traced. The different stiffnesses and inertia of the tank will be expressed by the wall thickness $(t)$ (t) variation. The results of the analysis are presented in Table 5. The natural frequencies determined for the acoustic fluid in the rigid tank in Example 2 are displayed in the last row for comparison. Reference values for $t = 0$ t=0 , shown in the first two rows, were determined in two ways: analytically and numerically for the boundary condition $ϕ = 0 on Γ_{F} = Γ_{1} \cup Γ_{2}$ \phi =0\hspace{0.25em}\text{on}{\Gamma }_{{\rm{F}}}={\Gamma }_{1}\cup {\Gamma }_{2} (Figure 1). A rigorous exact solution follows from equation (34) with the important restriction that the counter $n$ n starts from 1. For $n = 0$ n=0 , on a part of $Γ_{2}$ {\Gamma }_{2} , Neumann and Dirichlet boundary conditions are satisfied simultaneously, what must not appear in a physical problem. The numerical solution follows from equation (22) with the boundary division mentioned earlier. All the proper frequencies are highlighted in bold. In addition, the presented algorithm allows to examine impact of the tank inertia on the analysis results. Therefore, the calculations for each wall thickness were carried out in two variants: when inertia is neglected (italic font in the table) and is taken into account (normal font).

Table 5

Acoustic waves natural frequencies vs wall thickness (Hz).

$t (m)$ {t}\left(\text{m})	Variant	$f_{1}$ {f}_{1}	$f_{2}$ {f}_{2}	$f_{3}$ {f}_{3}	$f_{4}$ {f}_{4}	$f_{5}$ {f}_{5}	$f_{6}$ {f}_{6}
0.000	Exact		772.4			1,376	1,378
0.000	equation (22)	497.4	774.7	1,135	1,333	1,355	1,368
0.001	$B = 0$ {\boldsymbol{B}}={\bf{0}}	497.1	767.4	1,132	1,330	1,358	1,371
0.001	$B \neq 0$ {\boldsymbol{B}}\ne {\bf{0}}	498.3	763.0	1,135	1,351	1,381	1,382
0.002	$B = 0$ {\boldsymbol{B}}={\bf{0}}	498.9	774.7	1,136	1,328	1,362	1,374
0.002	$B \neq 0$ {\boldsymbol{B}}\ne {\bf{0}}	497.2	757.8	1,135	1,338	1,494	1,534
0.005	$B = 0$ {\boldsymbol{B}}={\bf{0}}	499.0	770.1	1,137	1,332	1,357	1,370
0.005	$B \neq 0$ {\boldsymbol{B}}\ne {\bf{0}}	498.3	710.4	1,135	1,288	1,488	1,527
0.010	$B = 0$ {\boldsymbol{B}}={\bf{0}}	499.0	785.4	1,138	1,336	1,349	1,365
0.010	$B \neq 0$ {\boldsymbol{B}}\ne {\bf{0}}	497.9	710.2	1,134	1,179	1,460	1,498
0.020	$B = 0$ {\boldsymbol{B}}={\bf{0}}	498.1	743.4	1,137	1,352	1,482	1,497
0.020	$B \neq 0$ {\boldsymbol{B}}\ne {\bf{0}}	491.4	653.0	1,107	1,142	1,465	1,503
0.050	$B = 0$ {\bf{B}}={\bf{0}}	499.5	796.3	1,107	1,135	1,218	1,311
0.050	$B \neq 0$ {\boldsymbol{B}}\ne {\bf{0}}	491.7	518.8	970.5	1,138	1,256	1,301
0.100	$B = 0$ {\boldsymbol{B}}={\bf{0}}	496.9	542.2	836.6	1,137	1,237	1,310
0.100	$B \neq 0$ {\boldsymbol{B}}\ne {\bf{0}}	379.5	500.3	873.6	1,137	1,210	1,307
0.200	$B = 0$ {\boldsymbol{B}}={\bf{0}}	345.8	500.7	751.0	1,137	1,491	1,497
0.200	$B \neq 0$ {\boldsymbol{B}}\ne {\bf{0}}	410.4	502.1	762.1	1,137	1,231	1,306
0.500	$B = 0$ {\boldsymbol{B}}={\bf{0}}	397.6	501.9	758.6	1,137	1,181	1,313
0.500	$B \neq 0$ {\boldsymbol{B}}\ne {\bf{0}}	391.7	503.4	779.2	1,137	1,213	1,309
	Rigid	401.0	502.0	772.6	1,137	1,209	1,311

Source: Author’s contribution.

The results of the calculations indicate that the tank stiffness has a noticeable effect on the natural frequency proper values of the acoustic waves. In addition, the inertia of the tank wall cannot be neglected when estimated. In contrast to that mentioned earlier, the tank mass does not modify improper values significantly and those frequencies are stable due to the wall thickness variation. In addition, the coupling between the two media is evident in a wide range of variations in the parameter $t$ t . Therefore, to confirm the convergence of the results with the limit solutions: determined in the previous examples or analytically, it was necessary to perform calculations for large wall thicknesses. Although the Kirchhoff’s theory is not adequate for the tank physical model when $t > 0.1 L_{y} = 0.09 m$ t\gt 0.1{L}_{y}=0.09\text{m} , it was utilized consequently for the sake of results consistence.

The coupling between the vibrations of a tank filled with liquid and acoustic waves makes it difficult to distinguish the natural frequencies of different genesis (Figure 7). Equation (44), which estimates at once surface waves and filled tank natural frequencies, is helpful in it (Figure 8). The comparison of figures allows to distinguish between them because the surface wave spectrum is under 4 Hz and the tank frequencies manifest clearly on the diagram. Also, as mentioned earlier, the stability of improper values is very helpful in the recognition process. Moreover, significant coupling appears only for lower tank frequencies. When those are below the acoustic wave spectrum, the coupling effect is weak – see Figures 9 and 10.

Example 5

Surface waves for incompressible fluid

The natural frequency estimation results of the surface waves obtained on the basis of the determinant sign analysis of the coefficient matrix defined in equation (44) are displayed in Table 6 in the same manner as in the previous example. The corresponding values for the incompressible fluid in the rigid tank, determined in Example 1, are displayed in the last row. The natural frequencies for $t = 0$ t=0 , designated analytically, are shown in the first row and take the same values as in the case of the rigid tank wall. Counter $j$ j in equation (32) must start from 1. The reason is the same as in the previous example. Numerical solution for that thickness requires other boundary surface assignments then those, which led to equation (30): $Γ_{F} = Γ_{1}$ {\Gamma }_{{\rm{F}}}={\Gamma }_{1} , $Γ_{R} = Γ_{3}$ {\Gamma }_{{\rm{R}}}={\Gamma }_{3} and different condition on the $Γ_{2}$ {\Gamma }_{2} : $ϕ = 0$ \phi =0 (Figure 1). It results an equation for natural frequencies estimation in the considered case: (45) $\det \begin{matrix} F_{0 Γ_{F}} - \frac{ω^{2}}{g} G_{0 Γ_{F}} & - G_{0 Γ_{2}} & F_{0 Γ_{R}} \end{matrix}] = 0 \to ω_{j} .$ \det \left[\begin{array}{ccc}{{\boldsymbol{F}}}_{0{\Gamma }_{F}}-\frac{{\omega }^{2}}{g}{{\boldsymbol{G}}}_{0{\Gamma }_{F}}& -{{\boldsymbol{G}}}_{0{\Gamma }_{2}}& {{\boldsymbol{F}}}_{0{\Gamma }_{R}}\end{array}\right]{\boldsymbol{=}}0\to {\omega }_{j}.

Table 6

Surface waves natural frequencies vs wall thickness (Hz).

$t (m)$ {t}\left(\text{m})	Variant		$f_{1}$ {f}_{1}	$f_{2}$ {f}_{2}	$f_{3}$ {f}_{3}	$f_{4}$ {f}_{4}	$f_{5}$ {f}_{5}	$f_{6}$ {f}_{6}
0.0000	Exact		0.8375	1.191	1.459	1.685	1.884	2.064
0.0000	Equation (45)		0.8417	1.203	1.481	1.720	1.934	2.131
0.0005	$B = 0$ {\boldsymbol{B}}={\bf{0}}	0.0505	0.8361	1.195	1.484	1.773	1.962	2.141
0.0005	$B \neq 0$ {\boldsymbol{B}}\ne {\bf{0}}	0.0502	0.8322	1.189	1.470	1.758	1.946	2.123
0.0010	$B = 0$ {\boldsymbol{B}}={\bf{0}}	0.1418	0.8416	1.292	1.525	1.757	1.921	2.095
0.0010	$B \neq 0$ {\boldsymbol{B}}\ne {\bf{0}}	0.1403	0.8317	1.277	1.508	1.735	1.897	2.070
0.0020	$B = 0$ {\boldsymbol{B}}={\bf{0}}	0.3793	0.9165	1.220	1.391	1.589	1.744	1.955
0.0020	$B \neq 0$ {\boldsymbol{B}}\ne {\bf{0}}	0.3732	0.9006	1.201	1.376	1.575	1.741	1.952
0.0050	$B = 0$ {\boldsymbol{B}}={\bf{0}}	0.7803	1.164	1.452	1.694	1.913	2.101	2.293
0.0050	$B \neq 0$ {\boldsymbol{B}}\ne {\bf{0}}	0.7787	1.161	1.449	1.688	1.906	2.078	2.268
0.0100	$B = 0$ {\boldsymbol{B}}={\bf{0}}	0.8341	1.199	1.478	1.718	1.933	2.132	2.317
0.0100	$B \neq 0$ {\boldsymbol{B}}\ne {\bf{0}}	0.8340	1.199	1.478	1.718	1.933	2.132	2.317
0.0500	$B = 0$ {\boldsymbol{B}}={\bf{0}}	0.8416	1.203	1.481	1.721	1.935	2.133	2.319
0.0500	$B \neq 0$ {\boldsymbol{B}}\ne {\bf{0}}	0.8416	1.203	1.481	1.721	1.935	2.133	2.319
	Rigid – Equation (30)	0.8417	1.203	1.481	1.721	1.935	2.133	2.319
		$f_{1}$ {f}_{1}	$f_{2}$ {f}_{2}	$f_{3}$ {f}_{3}	$f_{4}$ {f}_{4}	$f_{5}$ {f}_{5}	$f_{6}$ {f}_{6}	$f_{7}$ {f}_{7}

Source: Author’s contribution.

The results of the calculations indicate that the stiffness of the tank has a noticeable effect on the natural frequency values of the surface waves, whereas its inertia can be neglected. This is because for small wall thicknesses, their mass is negligible in relation to the mass of the fluid, while as the thickness increases, the vibration spectra of the filled tank and surface waves move away from each other and the coupling effect disappears.

4

Conclusions

The article presents a consistent method for determining the natural frequencies of surface and acoustic waves. For the cases of compressible and incompressible fluids, numerical model equations were formulated with boundary conditions describing perfectly rigid or flexible tank walls. The fluid domain discrete model equations were obtained using the boundary element method, which, compared to other methods, significantly reduces the number of degrees of freedom in numerical models.

The results obtained in the examples allow to formulate useful advices in future research of the discussed issue.

Acceptable values of natural frequencies for surface and acoustic waves are obtained even with a small number of degrees of freedom.
Identification of improper natural frequencies for acoustic waves must be done carefully, as it does not always give satisfactory results. Determination of the eigenmodes can be useful in distinguishing between them.
The coupling effect is not important for real fluids, so that different wave types can be analyzed separately: the surface waves for incompressible fluid model, the acoustic waves without surface ones presence.
The stiffness of the tank walls significantly affects the correct natural frequency values for both types of waves, while their inertia is important only for acoustic waves.
The improper values of natural frequencies determined in the acoustic wave analysis are insensitive to changes in the stiffness and inertia of the tank’s walls.
In some specific cases, it will be worth considering applying more accurate tank physical models than those used in the examples presented in the article.
When coupling occurs between the vibrations of the tank and the waves in the fluid, the simultaneous analysis of spectra involving surface and acoustic waves makes it possible to distinguish between natural frequencies of different genesis.

The consistency and simplicity of the discrete equations for different physical models: compressible or incompressible fluid, rigid or flexible tank, and the quality of the results in each case confirm the versatility of the proposed method. In the author’s opinion, the aforementioned algorithm features and practical tips formulated on the basis of the numerical examples provided are a contribution to the analysis methods development of the issue under consideration.

Finally, it will be worth mentioning that the presented technique can be applied to the three-dimensional problems as well. In that case, any order fundamental solutions of the Laplace equation take the form displayed in equation (18).

Funding information

Author states no funding involved.

Author contributions

The author confirms the sole responsibility for the conception of the study, presented results and manuscript preparation.

Conflict of interest statement

Author states no conflict of interest.

Surface and acoustic waves natural frequencies estimation by boundary element method

Full Article

Paradigm

My account