# A Computational Method for Acoustic Interaction with Large Complicated Underwater Structures Based on the Physical Mechanism of Structural Acoustics

Authors: Yongzhuang Tang; Qidou Zhou; Xiaowei Wang; Zhiyong Xie
Journal: Advances in Materials Science and Engineering
 (2022)

---

## Abstract

A numerical coupling approach is proposed to fast predict the acoustic radiation from a vibrating large-complicated underwater structure. In this study, the physical mechanism of sound radiation from underwater large target is used for the first time to improve the efficiency and keep the accuracy of the numerical algorithm. Although the traditional coupled finite element method/boundary element method (FEM-BEM) is accurate, it contains a large number of boundary elements and thus requires a long computation time for large-complicated structures. The research on the physical mechanism of structural acoustics shows that when BEM is applied on the near-field artificial boundary at a proper distance away from the wet structural surface, large-size boundary elements are acceptable and the number of boundary elements and computation time are remarkably reduced. Thus, the fluid outside the structure is divided into the interior domain and the exterior domain by the artificial boundary. Then, the numerical method is realized by coupling structural finite element modelling with interior fluid finite element modelling and with exterior fluid boundary element modelling. Compared with the theoretical value, the experimental value and the results of the traditional FEM-BEM, the correctness of the proposed algorithm and its advantage of computational efficiency are verified. The computation time of the proposed method is over 99% shorter than that of FEM-BEM in the calculation example of a large-complicated structure. The proposed method can be further applied to multidomain acoustic and multibody acoustic calculations.


---

## Body

## 1. Introduction

Because of the concealment requirements of military targets such as submarines, the study on acoustic radiation from vibrating underwater structures has been a focus in structural acoustics for decades [1, 2]. Fast and accurate numerical simulation is an effective way to study the acoustic radiation characteristics and control strategies of large underwater structures, so the acoustic radiation algorithm attracts more research interests. Due to the strong acoustic-structure coupling, the elastic and acoustic wave problems must be solved simultaneously [3]. For structures with a simple geometry, analytical solutions are available [4]. Acoustic calculation in general needs to solve wave equation, which belongs to partial differential equation and boundary value problem. Alshehry et al.[5] recently proposed a new reliable technique for solving fractional partial differential equations. The residual power series method and the Laplace transform are combined in the suggested technique. Muhammad and Syed [6] applied the reduced differential transform method (RDTM) to solve the first order linear and nonlinear system of fractional PDEs. More methods for solving diffusion partial differential equations are widely studied in engineering problems such as heat transfer and flow [7, 8]. Generally, complex structural acoustic are investigated by numerical methods [9], such as the finite element method (FEM), infinite element method (IEM), and the coupled finite element method/boundary element method (FEM- BEM) [10–12]. BEM not only satisfies the Sommerfeld radiation condition by the fundamental solution but also establishes an independent coupling equation between the structure and fluid. Thus, FEM-BEM is usually employed to predict the noise radiated from underwater vibrating structures due to exciting forces acting on the structures [13–15].However, fully populated matrices resulting from the coupling of FEM and BEM make its computational efficiency lower than others. For large-complicated structures, even at low frequencies, the long computation time limits the application of FEM-BEM. The problem is brought by the extremely large number of boundary elements at the fluid-structure interface because the boundary element size at the interface has to be conformed [16] to the finite element size on the wet structural surface. The existing fast BEM [17–20] uses special iterative solvers to speed up matrix operations and to reduce memory consumption, which remains solving the problem on the structural surface and contains a large number of boundary elements. Besides, the FE-mixed FE-BE coupling method is in the line of using a virtual surface [21, 22]. Nevertheless, the existing numerical algorithms do not seem to improve the computational efficiency from the perspective of the sound radiation mechanism. Different from previous studies, the present study on the physical mechanism of acoustic radiation from large structures shows that it is feasible to reduce the computation time by using less boundary elements when their size becomes larger than in the existing method.In recent years, researchers have performed extensive simulations on acoustic radiation from large complicated underwater structures in order to find out the noise reduction mechanism of structures [23]. It has been found that for large, complex structures, as shown in Figure 1, the structural vibration is usually a summation of vibrations wave components of different structural wavelengths even if the excited force is harmonic. The far-field noise mainly comes from the structural vibration with a relatively long wavelength while short waves only exist in the near field around the structure which satisfies the kinematic continuity condition at the structure-fluid interface and form the evanescent waves from the near field. When the structural bending wavelength λb on the wet surface is smaller than the acoustic wavelength λ in the fluid domain, the evanescent waves [24], decayed exponentially with distance rs and normal to the fluid-structure interface, are formed in the near field around the wet surface. Especially, when λ >> λb, the evanescent waves decay as exp(−2πrs/λb). This means that the short waves only exist in the near field around the wet surface. Therefore, the FEM for fluids can be employed to well simulate the near-field acoustic waves in the interior region Ω1 as presented in Figure 2, as the BEM will lead to fully populated matrices, leading to unacceptable computation time.Figure 1
The acoustic radiation model of large, complicated structure vibrating at low frequencies.Figure 2
The geometry of the structural acoustic problem of the FE-FE-BE method.Based on the fact that the short acoustic waves only exist in the near field around the structure for the low-frequency structural vibrations, and that the exterior outgoing acoustic wavelength is long at the low frequency, the FE-FE-BE method is proposed in order to divide the acoustic fluid domain into an interior domain Ω1 and an exterior domain Ω0 by the artificial boundary surface ΓA as shown in Figure 2. To obtain the dynamic response of the structure, the FEM for structures is coupled with the FEM for fluids in the interior fluid domain Ω1, and the BEM for fluids is adopted to simulate the acoustic radiation in the exterior domain Ω0. The coupling effects of FEM in the domain Ω1 and BEM in the domain Ω0 is included by satisfying the pressure continuity condition on the artificial boundary surface ΓA. On ΓA, the large element size is used for BEM because of the long acoustic wavelength. Both the number of elements and the computation time are remarkably reduced.In addition, the author presented the initial idea of artificial boundary element in the 24th International Congress on Sound and Vibration (ICSV) and discussed the prospects for future research [23]. However, the formulations were not given and the code did not complete at that time. In the present work, the theoretical approach to the FE-FE-BE method and the complete simulation algorithm will be presented in Section 2. Furthermore, the fluid added stiffness matrix and mass matrix on the artificial boundary are discussed. In Section 3, the implementation steps of the algorithm are given. The accuracy of the proposed method is verified by two examples with theoretical results and experimental values, respectively, in Section 4. The efficiency is confirmed by the example of a large structure in Section 5, and the results of vibration and sound pressure are further discussed.


## 2. Theoretical Approaches

### 2.1. BE Equations for Exterior Domain Fluid

In the exterior domain, the acoustic pressurep in the fluid satisfies the wave equation, the Sommerfeld condition and the boundedness condition at infinity. Besides, acoustic pressure or normal velocity continuity is required at the artificial boundary surface ΓA. The following equation constitutes the governing equation for fluid in the exterior domain:(1a)∂2p∂t2=c02∇2p,(1b)limr⟶∞r∂p∂r+kω∂p∂t=0,(1c)limrp<∞,r⟶∞(1d)p=p1or∂p∂n=∂p1∂nonΓA,where ω is the circular disturbance frequency, c0 is the speed of sound in the fluid, the wave number is k = ω/c0, r is the distance from the field point x (x1, x2, x3) to the source point y (y1, y2, y3), p1 is the interior sound pressure, ρ is the fluid density, and n is the outward normal to the interior fluid. The pressure p is the solution of the Helmholtz integral equation:(2)∫ΓApy∂Gx,y∂ndS−∫ΓA∂py∂nGx,ydS=px2,xonΓA,px,xinΩ0,0,xinΩ1,(3)Gx,y=e−ikr4πr,where G (x and y) is the free-space Green’s function. When the boundary condition at the free water surface is required, Green’s function is as follows:(4)Gx,y=e−ikr4πr−e−ikr14πr1.And the normal derivative of Green’s functionG can be calculated as follows:(5a)∂Gx,y∂n=∂Gx,y∂r∂r∂n=e−ikr4πrik+1rcosα,nofreewatersurface,(5b)∂Gx,y∂n=e−ikr4πrik+1rcosα−e−ikr14πr1ik+1r1cosα1,withfreewatersurface,where i = (−1)1/2, r = x − y, r = |r| = [(x1 − y1)2 + (x2 − y2)2 + (x3 − y3)2]1/2, r1 = x − y1, and y1(y1, y2, −y3) are the mirror points of y about free water surface; r1 = |r1| = [(x1 − y1)2 + (x2 − y2)2 + (x3 + y3)2]1/2, and y3 denotes the axis pointing upward to the free water surface; α is the angle between the normal n and the vector r; and α1 is the angle between the normal n and the vector r1.In the first formula of equation (2), the field point x and the source point y are both on the surface ΓA. To operate numerical computation, the surface integral over ΓA is discretized as the sum of integral over each small surface element:(6)∑i=1Ne∫Sipy∂Gx,y∂ndS−∑i=1Ne∫Si∂py∂nGx,ydS=px2,where Ne is the number of small surface elements (boundary elements). When the element size is smaller than the acoustic wavelength, the integral for Si could be replaced by the multiplication of element area and medium integral (i.e., average):(7)∑i=1Nep¯i∂Gx,y∂n∆Si−p¯j2=∑i=1Ne∂p¯i∂nGx,y∆Si,where ΔSi is the area of element i. The source point y is on the element i, the field point x is on the element j. Equation (7) means that the pressure of element j could be represented by the information of elements numbered from 1 to Ne. Then, we obtain the matrix equation for the artificial surface ΓA:(8)Hp¯=Q∂p¯∂n.Using the equation (1d), the equation (8) is written as follows:(9)Hp¯1=Q∂p¯∂n,where {} represents row-by-row or column-by-column matrix assembly, the bold symbols are matrices. The element of matrix H is as follows:(10)Hij=∂Gij∂n∆Si,i≠j,∂Gij∂n∆Si−12,i=j.Using equation (5), though without free water surface, equation (10) can be written as follows:(11a)Hij=∆Sicosα4πr2coskr−krsinkr+isinkr+krcoskr,i≠j,−12,i=j.Considering free water surface, equation (10) can be written as follows:(11b)Hij=∆Sicosα4πr2coskr−krsinkr+isinkr+krcoskr−∆Sicosα14πr12coskr1−kr1sinkr1+isinkr1+kr1coskr1,i≠j,−∆Sicosα14πr12coskr1−kr1sinkr1+isinkr1+kr1coskr1−12,i=j.The element of matrixQ is Qij = G∆Si. When without free water surface, it can be written as follows:(12a)Qij=coskr−isinkr4πr∆Si,i≠j,∫02π∫0bicoskr−isinkr4πrrdrdθ,i=j.In equation (12a), when i = j, the source point and field point are on the same element. The element is equivalent to a circular surface and bi is the equivalent radius, so bi=(∆Si/π)1/2. Assume that the source point is at the element center and the field point is anywhere on the element. After integral, Qii = [sin(kbi) + i(cos(kbi) − 1)]/2k. Considering free water surface, Qij can be written as follows:(12b)Qij=coskr−isinkr4πr∆Si−coskr1−isinkr14πr1∆Si,i≠j,sinkbi+icoskbi−12k−coskr1−isinkr14πr1∆Si,i=j.


### 2.2. FE Equations for Fluid in the Interior Domain

Generally, the FEM for fluids is different from the FEM for structures as the fluid pressure is a scalar but the structural displacement is a vector. Traditionally, the FE equations for fluids could be derived by variational principle and weight residual method [25]. The fluid finite elements can alternatively be created from elastic finite elements by the use of analogy between the wave equation and the equations of elasticity [26]. In this section, the interactive “forces” between the interior fluid, the structure, and the exterior fluid are developed by analogy, which is possibly concise and easy to code. Based on the variational principle or analogy, the finite element equation for fluid has a general form:(13)−ω2Mf+Kfp1=0,where Mf and Kf are the fluid mass matrix and stiffness matrix. The wave equation for fluid could also be written as follows:(14)∇2p1=1c02∂2p1∂t2.Thex-component of Navier equations of elasticity is as follows:(15)λ+2GeGe∂2u∂x2+∂2u∂y2+∂2u∂z2+λ+GeGe∂2v∂x∂y+∂2w∂x∂z+1Gefx=ρsGe∂2u∂t2,(16)Ge=E21+ν,λ=Eν1+ν1−2ν,where u, v, and w are the components of displacement in the Cartesian coordinate system, λ is the Lame constants, Ge is the shear modulus, E is the modulus of elasticity, ν is the Poisson’s ratio, fx is the body force per unit volume of x-component, ρs is the density of elastomer. Comparing equation (14) with equation (15), it is not difficult to discover that if the constants of equation (17) are used in equation (15), equation (15) reduces the same form of equation (14) and u is equivalent with p1.(17a)λ+2GeGe=1,(17b)λ+GeGe=0,(17c)1Gefx=0,(17d)ρsGe=1c02.Now, the fluid could be analyzed in the same way as the structure. According to the physical equations of elasticity, the normal stress on the interface is described by the following:(18)Tn=σijn,where n is the outward normal to the interface at the node, and σij is the stress tensor. Only x-component exists in the analogy, so equation (18) gives the following:(19)Tn=σxxnx+σxyny+σxznz.Based on the physical and geometric equations of elasticity, the relationship between stress tensorσ and strain tensor ε is written as follows:(20a)σ=Dε,(20b)ε=∂u∂x,∂v∂y,∂w∂z,∂u∂y+∂v∂x,∂v∂z+∂w∂y,∂w∂x+∂u∂zT,(20c)D=E1+ν1−2ν1−νννν1−νννν1−ν1−2ν21−2ν21−2ν2.In equation (17b), λ/Ge = −1 and then equation (16) are written as (1 − 2ν)/2ν = −1, where ν is an infinite value. Hence, the elastic constant D could be written as follows:(21)D=E21+ν1−1−1−11−1−1−11111=Ge1−1−1−11−1−1−11111.In equation (20b), only u exists for fluid. Then σxx, σxy, and σxz are written as follows:(22)σxx=Ge∂u∂x,σxy=Ge∂u∂y,σxz=Ge∂u∂z.For fluid in NASTRAN, the shear modulusGe = 1 [27]. Using equation (22) and the gradient of u, equation (19) yields(23)Tn=∂u∂xnx+∂u∂yny+∂u∂znz=∂u∂n.Equation (23) gives the stress vector acting on a surface with the unit outward normal. The surface could be the fluid-structure interface or the artificial boundary between the interior and exterior fluids. If the surface is discretized with finite elements, the force Fn applied to the node can be replaced by Tn × A, where A is the area assigned to the node. And p1 is equivalent to u. Hence, the final expression of the nodal force is as follows:(24)Fn=A∂p1∂n.The exterior fluid load and the structural load are determined by the boundary conditions on the artificial surface and wet structural surface. On the boundary ΓS, the relationship is given by the following:(25)∂p1∂n=iωρvn.Thus, the structural load acted on the interior fluid becomes(26)Fs=A∂p1∂n=ρAsu¨onΓS,where As is the area matrix on Γs. On the artificial boundary ΓA, the loading of the exterior fluid acted on the interior fluid can be written as follows:(27)Fa=Aa∂p1∂nonΓA,where Aa is the area matrix on ΓA and {∂p1/∂n} is dependent on the pressure on the boundary ΓA. Considering the structural and exterior fluid loads acted on the interior fluid, equation (13) becomes(28)−ω2Mf+Kfp1=Fs+Fa.


### 2.3. FE Dynamic Equation for the Structure

Considering an arbitrary elastic structure modelled with finite elements, including time-dependent loads and pressure from interior fluid Ω1, the resulting dynamic equation of FEM for structure is as follows:(29)−ω2Ms+iωBs+Ksu=F−AsTp1,where Ms, Bs, and Ks represent the structure mass, viscous damping, and stiffness matrix, u is the structural displacement vector, the second term at the right of the equation is the nodal force translated from interior fluid pressure p1, and F is the load on structure.


### 2.4. Coupled Fluid-Structure Equations

There are two coupling relationships in the fluid-structure system, the structure is coupled with the interior fluid, and the interior fluid is coupled with the exterior fluid. Combining equation (29) for structural dynamics with equation (28) for the interior fluid yields the coupled fluid-structure equations:(30)Ms0−ρAsMfu¨p¨1+Bs000u˙p˙1+KsAsT0Kfup1=FFa,where there are three groups of unknowns ({u}, {p1}, and {Fa}) in the equation. To solve it, equation (9) for boundary elements has to be included. However, ∂p1/∂n in Fa is the nodal value, and p1¯ is average value for the element. It is necessary to obtain the relationship between nodal value and average value for the element. The transformation matrices G and L are defined to build the connections.(31a)G∂p¯1∂n=∂p¯1∂n,(31b)p¯1=Lp1.On the artificial boundary, triangular elements are modelled. Thus, the average pressure of any element could be solved by interpolating the pressure of three element nodes. As shown in Figure3, for element j, the average sound pressure pj¯ can be calculated by pj¯ = Lj{pj1, pj2, pj3}, where Lj is the special case in( [15], equation (16)) and Lj = [1/3, 1/3, 1/3]. The matrix L (Ne × Ng) is assembled by Lj(j = 1, 2, 3, …, Ne) and G (Ng × Ne) is the transposition of L, where Ng is the number of boundary element nodes.Figure 3
The acoustic pressure on the element j.Combing equation (31) with equation (9) gives the final form of {Fa}:(32)Fa=AaGQ−1HLp1.Obviously, {Fa} is the matrix containing complex numbers which can be written as the sum of real and imaginary parts:(33)Fa=−ω2Mfa−iωBfap1,where(34a)Mfa=ReAaGQ−1HLω2,(34b)Bfa=−ImAaGQ−1HLω.The final coupled fluid-structure FE equation considering the exterior fluid is modified as follows:(35)Ms0−ρAsMf+Mfau¨p¨1+Bs00Bfau˙p˙1+KsAsT0Kfup1=F0.The solution of this equation gives the structural dynamic response {u} and the acoustic pressure of the interior fluid {p1}, where {p1} contains the nodal acoustic pressure on the artificial boundary. The matrices defined by equation (34) may be interpreted as the added mass and damping of the exterior fluid for the interior fluid, which is developed with the concept of fluid element mass and stiffness in FEM. Comparing the added mass and damping is to obtain the effects of fluid on the structure in FEM-BEM in [15], and in this paper is to see the effects of the exterior fluid Ω0 on the interior fluid Ω1.


### 2.5. The Calculation of Acoustic Pressure in the Exterior Domain

Once the interior fluid pressure {p1} is obtained from the solution of equation (35), the acoustic pressure at any field point x in the exterior domain can be computed by BEM. The second equation of Helmholtz integral equation (equation (2)) is discretized for numerical computation according to the process of equations (6)–(9):(36)px=H′p¯y−Q′∂p¯y∂n,xinΩ0yonΓΑ.With the continuity condition of pressure on ΓA, and based on equations (9) and (31b), equation (36) becomes(37)px=H′Lp1y−Q′Q−1HLp1y,where {p1(y)} can be extracted from the interior fluid pressure {p1}, the far-field point x does not coincide with the source point y.(38)Hij′=∂Gx,y∂n∆Si,(39)Qij′=Gx,y∆Si,where i = 1, 2, 3, …, Ne; j = 1, 2, 3, …, Ng0; Ne is the number of boundary elements; and Ng0 is the number of acoustic field points.


## 2.1. BE Equations for Exterior Domain Fluid

In the exterior domain, the acoustic pressurep in the fluid satisfies the wave equation, the Sommerfeld condition and the boundedness condition at infinity. Besides, acoustic pressure or normal velocity continuity is required at the artificial boundary surface ΓA. The following equation constitutes the governing equation for fluid in the exterior domain:(1a)∂2p∂t2=c02∇2p,(1b)limr⟶∞r∂p∂r+kω∂p∂t=0,(1c)limrp<∞,r⟶∞(1d)p=p1or∂p∂n=∂p1∂nonΓA,where ω is the circular disturbance frequency, c0 is the speed of sound in the fluid, the wave number is k = ω/c0, r is the distance from the field point x (x1, x2, x3) to the source point y (y1, y2, y3), p1 is the interior sound pressure, ρ is the fluid density, and n is the outward normal to the interior fluid. The pressure p is the solution of the Helmholtz integral equation:(2)∫ΓApy∂Gx,y∂ndS−∫ΓA∂py∂nGx,ydS=px2,xonΓA,px,xinΩ0,0,xinΩ1,(3)Gx,y=e−ikr4πr,where G (x and y) is the free-space Green’s function. When the boundary condition at the free water surface is required, Green’s function is as follows:(4)Gx,y=e−ikr4πr−e−ikr14πr1.And the normal derivative of Green’s functionG can be calculated as follows:(5a)∂Gx,y∂n=∂Gx,y∂r∂r∂n=e−ikr4πrik+1rcosα,nofreewatersurface,(5b)∂Gx,y∂n=e−ikr4πrik+1rcosα−e−ikr14πr1ik+1r1cosα1,withfreewatersurface,where i = (−1)1/2, r = x − y, r = |r| = [(x1 − y1)2 + (x2 − y2)2 + (x3 − y3)2]1/2, r1 = x − y1, and y1(y1, y2, −y3) are the mirror points of y about free water surface; r1 = |r1| = [(x1 − y1)2 + (x2 − y2)2 + (x3 + y3)2]1/2, and y3 denotes the axis pointing upward to the free water surface; α is the angle between the normal n and the vector r; and α1 is the angle between the normal n and the vector r1.In the first formula of equation (2), the field point x and the source point y are both on the surface ΓA. To operate numerical computation, the surface integral over ΓA is discretized as the sum of integral over each small surface element:(6)∑i=1Ne∫Sipy∂Gx,y∂ndS−∑i=1Ne∫Si∂py∂nGx,ydS=px2,where Ne is the number of small surface elements (boundary elements). When the element size is smaller than the acoustic wavelength, the integral for Si could be replaced by the multiplication of element area and medium integral (i.e., average):(7)∑i=1Nep¯i∂Gx,y∂n∆Si−p¯j2=∑i=1Ne∂p¯i∂nGx,y∆Si,where ΔSi is the area of element i. The source point y is on the element i, the field point x is on the element j. Equation (7) means that the pressure of element j could be represented by the information of elements numbered from 1 to Ne. Then, we obtain the matrix equation for the artificial surface ΓA:(8)Hp¯=Q∂p¯∂n.Using the equation (1d), the equation (8) is written as follows:(9)Hp¯1=Q∂p¯∂n,where {} represents row-by-row or column-by-column matrix assembly, the bold symbols are matrices. The element of matrix H is as follows:(10)Hij=∂Gij∂n∆Si,i≠j,∂Gij∂n∆Si−12,i=j.Using equation (5), though without free water surface, equation (10) can be written as follows:(11a)Hij=∆Sicosα4πr2coskr−krsinkr+isinkr+krcoskr,i≠j,−12,i=j.Considering free water surface, equation (10) can be written as follows:(11b)Hij=∆Sicosα4πr2coskr−krsinkr+isinkr+krcoskr−∆Sicosα14πr12coskr1−kr1sinkr1+isinkr1+kr1coskr1,i≠j,−∆Sicosα14πr12coskr1−kr1sinkr1+isinkr1+kr1coskr1−12,i=j.The element of matrixQ is Qij = G∆Si. When without free water surface, it can be written as follows:(12a)Qij=coskr−isinkr4πr∆Si,i≠j,∫02π∫0bicoskr−isinkr4πrrdrdθ,i=j.In equation (12a), when i = j, the source point and field point are on the same element. The element is equivalent to a circular surface and bi is the equivalent radius, so bi=(∆Si/π)1/2. Assume that the source point is at the element center and the field point is anywhere on the element. After integral, Qii = [sin(kbi) + i(cos(kbi) − 1)]/2k. Considering free water surface, Qij can be written as follows:(12b)Qij=coskr−isinkr4πr∆Si−coskr1−isinkr14πr1∆Si,i≠j,sinkbi+icoskbi−12k−coskr1−isinkr14πr1∆Si,i=j.


## 2.2. FE Equations for Fluid in the Interior Domain

Generally, the FEM for fluids is different from the FEM for structures as the fluid pressure is a scalar but the structural displacement is a vector. Traditionally, the FE equations for fluids could be derived by variational principle and weight residual method [25]. The fluid finite elements can alternatively be created from elastic finite elements by the use of analogy between the wave equation and the equations of elasticity [26]. In this section, the interactive “forces” between the interior fluid, the structure, and the exterior fluid are developed by analogy, which is possibly concise and easy to code. Based on the variational principle or analogy, the finite element equation for fluid has a general form:(13)−ω2Mf+Kfp1=0,where Mf and Kf are the fluid mass matrix and stiffness matrix. The wave equation for fluid could also be written as follows:(14)∇2p1=1c02∂2p1∂t2.Thex-component of Navier equations of elasticity is as follows:(15)λ+2GeGe∂2u∂x2+∂2u∂y2+∂2u∂z2+λ+GeGe∂2v∂x∂y+∂2w∂x∂z+1Gefx=ρsGe∂2u∂t2,(16)Ge=E21+ν,λ=Eν1+ν1−2ν,where u, v, and w are the components of displacement in the Cartesian coordinate system, λ is the Lame constants, Ge is the shear modulus, E is the modulus of elasticity, ν is the Poisson’s ratio, fx is the body force per unit volume of x-component, ρs is the density of elastomer. Comparing equation (14) with equation (15), it is not difficult to discover that if the constants of equation (17) are used in equation (15), equation (15) reduces the same form of equation (14) and u is equivalent with p1.(17a)λ+2GeGe=1,(17b)λ+GeGe=0,(17c)1Gefx=0,(17d)ρsGe=1c02.Now, the fluid could be analyzed in the same way as the structure. According to the physical equations of elasticity, the normal stress on the interface is described by the following:(18)Tn=σijn,where n is the outward normal to the interface at the node, and σij is the stress tensor. Only x-component exists in the analogy, so equation (18) gives the following:(19)Tn=σxxnx+σxyny+σxznz.Based on the physical and geometric equations of elasticity, the relationship between stress tensorσ and strain tensor ε is written as follows:(20a)σ=Dε,(20b)ε=∂u∂x,∂v∂y,∂w∂z,∂u∂y+∂v∂x,∂v∂z+∂w∂y,∂w∂x+∂u∂zT,(20c)D=E1+ν1−2ν1−νννν1−νννν1−ν1−2ν21−2ν21−2ν2.In equation (17b), λ/Ge = −1 and then equation (16) are written as (1 − 2ν)/2ν = −1, where ν is an infinite value. Hence, the elastic constant D could be written as follows:(21)D=E21+ν1−1−1−11−1−1−11111=Ge1−1−1−11−1−1−11111.In equation (20b), only u exists for fluid. Then σxx, σxy, and σxz are written as follows:(22)σxx=Ge∂u∂x,σxy=Ge∂u∂y,σxz=Ge∂u∂z.For fluid in NASTRAN, the shear modulusGe = 1 [27]. Using equation (22) and the gradient of u, equation (19) yields(23)Tn=∂u∂xnx+∂u∂yny+∂u∂znz=∂u∂n.Equation (23) gives the stress vector acting on a surface with the unit outward normal. The surface could be the fluid-structure interface or the artificial boundary between the interior and exterior fluids. If the surface is discretized with finite elements, the force Fn applied to the node can be replaced by Tn × A, where A is the area assigned to the node. And p1 is equivalent to u. Hence, the final expression of the nodal force is as follows:(24)Fn=A∂p1∂n.The exterior fluid load and the structural load are determined by the boundary conditions on the artificial surface and wet structural surface. On the boundary ΓS, the relationship is given by the following:(25)∂p1∂n=iωρvn.Thus, the structural load acted on the interior fluid becomes(26)Fs=A∂p1∂n=ρAsu¨onΓS,where As is the area matrix on Γs. On the artificial boundary ΓA, the loading of the exterior fluid acted on the interior fluid can be written as follows:(27)Fa=Aa∂p1∂nonΓA,where Aa is the area matrix on ΓA and {∂p1/∂n} is dependent on the pressure on the boundary ΓA. Considering the structural and exterior fluid loads acted on the interior fluid, equation (13) becomes(28)−ω2Mf+Kfp1=Fs+Fa.


## 2.3. FE Dynamic Equation for the Structure

Considering an arbitrary elastic structure modelled with finite elements, including time-dependent loads and pressure from interior fluid Ω1, the resulting dynamic equation of FEM for structure is as follows:(29)−ω2Ms+iωBs+Ksu=F−AsTp1,where Ms, Bs, and Ks represent the structure mass, viscous damping, and stiffness matrix, u is the structural displacement vector, the second term at the right of the equation is the nodal force translated from interior fluid pressure p1, and F is the load on structure.


## 2.4. Coupled Fluid-Structure Equations

There are two coupling relationships in the fluid-structure system, the structure is coupled with the interior fluid, and the interior fluid is coupled with the exterior fluid. Combining equation (29) for structural dynamics with equation (28) for the interior fluid yields the coupled fluid-structure equations:(30)Ms0−ρAsMfu¨p¨1+Bs000u˙p˙1+KsAsT0Kfup1=FFa,where there are three groups of unknowns ({u}, {p1}, and {Fa}) in the equation. To solve it, equation (9) for boundary elements has to be included. However, ∂p1/∂n in Fa is the nodal value, and p1¯ is average value for the element. It is necessary to obtain the relationship between nodal value and average value for the element. The transformation matrices G and L are defined to build the connections.(31a)G∂p¯1∂n=∂p¯1∂n,(31b)p¯1=Lp1.On the artificial boundary, triangular elements are modelled. Thus, the average pressure of any element could be solved by interpolating the pressure of three element nodes. As shown in Figure3, for element j, the average sound pressure pj¯ can be calculated by pj¯ = Lj{pj1, pj2, pj3}, where Lj is the special case in( [15], equation (16)) and Lj = [1/3, 1/3, 1/3]. The matrix L (Ne × Ng) is assembled by Lj(j = 1, 2, 3, …, Ne) and G (Ng × Ne) is the transposition of L, where Ng is the number of boundary element nodes.Figure 3
The acoustic pressure on the element j.Combing equation (31) with equation (9) gives the final form of {Fa}:(32)Fa=AaGQ−1HLp1.Obviously, {Fa} is the matrix containing complex numbers which can be written as the sum of real and imaginary parts:(33)Fa=−ω2Mfa−iωBfap1,where(34a)Mfa=ReAaGQ−1HLω2,(34b)Bfa=−ImAaGQ−1HLω.The final coupled fluid-structure FE equation considering the exterior fluid is modified as follows:(35)Ms0−ρAsMf+Mfau¨p¨1+Bs00Bfau˙p˙1+KsAsT0Kfup1=F0.The solution of this equation gives the structural dynamic response {u} and the acoustic pressure of the interior fluid {p1}, where {p1} contains the nodal acoustic pressure on the artificial boundary. The matrices defined by equation (34) may be interpreted as the added mass and damping of the exterior fluid for the interior fluid, which is developed with the concept of fluid element mass and stiffness in FEM. Comparing the added mass and damping is to obtain the effects of fluid on the structure in FEM-BEM in [15], and in this paper is to see the effects of the exterior fluid Ω0 on the interior fluid Ω1.


## 2.5. The Calculation of Acoustic Pressure in the Exterior Domain

Once the interior fluid pressure {p1} is obtained from the solution of equation (35), the acoustic pressure at any field point x in the exterior domain can be computed by BEM. The second equation of Helmholtz integral equation (equation (2)) is discretized for numerical computation according to the process of equations (6)–(9):(36)px=H′p¯y−Q′∂p¯y∂n,xinΩ0yonΓΑ.With the continuity condition of pressure on ΓA, and based on equations (9) and (31b), equation (36) becomes(37)px=H′Lp1y−Q′Q−1HLp1y,where {p1(y)} can be extracted from the interior fluid pressure {p1}, the far-field point x does not coincide with the source point y.(38)Hij′=∂Gx,y∂n∆Si,(39)Qij′=Gx,y∆Si,where i = 1, 2, 3, …, Ne; j = 1, 2, 3, …, Ng0; Ne is the number of boundary elements; and Ng0 is the number of acoustic field points.


## 3. Summary of the Theory and Numerical Program

For the solution procedure, the FEM software NASTRAN is employed to generate the matricesMs, Bs, Ks, Mf, Kf, As, and F. If the coupled effect from the exterior fluid (Fa) is not considered, equation (30) can be solved by NASTRAN completely. To solve the coupled fluid-structure equation (equations (35)–(37)) for the exterior sound, the information of boundary elements and field points are generated by preprocessing software PATRAN so that the matrices Q, H, G, L, Aa, H′, and Q′ can be calculated by a separate self-written code. The matrix operations are solved by direct matrix abstraction program (DMAP), a programming language of NASTRAN. The procedure of the FE-FE-BE method is as follows:(I)
Building the finite element model with the preprocessor PATRAN, and meanwhile generating the initial calculation file and separate information files of boundary elements and acoustic field points.(II)
Generating the matricesQ, H, G, L, Aa, H′, and Q′ with the self-written code.(III)
Solving equation (34) to obtain Mfa and Bfa with the DMAP language.(IV)
Adding the added mass and damping matrices of the exterior fluid to the mass and damping matrices of the interior fluid in the initial calculation file with the DMAP language and getting a new calculation file.(V)
Running the SOL108 of NASTRAN to solve the new calculation file and outputting the nodal pressure on the artificial boundary.(VI)
Solving equation (37) with the DMAP language to give the acoustic pressure of field points. The results can be visually displayed in PATRAN.


## 4. Validation with Theoretical and Experimental Data

### 4.1. A Uniformly Driven Spherical Shell

The sound pressure radiated from an underwater spherical shell driven by a uniform pressure unit is calculated, and the results are compared with the theoretical data [13]. The pressure load is p0 = cos (2πft). The structural and fluid characteristics are shown in Table 1, where a is the radius of the spherical shell, t is the thickness of the shell, ν is Poisson’s ratio, E is Young’s modulus, ρs is the density of structure, ρf is the density of water, and c0 is the speed of sound in water.Table 1
Properties of materials.
a (m)t (m)νE (Pa)ρs (kg·m−3)ρf (kg·m−3)c0 (m·s−1)0.50.010.32.07e11766910001524The spherical shell is meshed with 3200 triangular elements and 1602 nodes, as shown in Figure4(a). The artificial boundary is a spherical surface with radius R = 2 m, which is also meshed with 3200 triangular elements as shown in Figure 4(b). The interior fluid is meshed with CTETRA elements between the artificial boundary elements and wet interface elements as shown in Figure 4(b). Figure 4(c) depicts the 360 acoustic field points on a circle with a radius of 100 m.Figure 4
(a) Structural FE meshes; (b) interior fluid FE meshes; and (c) acoustic field points.
(a)(b)(c)The numerical model calculates in the range from 100 Hz to 2000 Hz. Figure5 shows the comparison between the numerical and theoretical results, where the former is the average acoustic pressure on 360 field points. The numerical results agree well with the theoretical results in the calculation range.Figure 5
Results of the acoustic pressure (10−4 Pa).


### 4.2. A Ribbed Cylindrical Shell Driven by Force

The sound pressure radiated from an underwater ribbed cylindrical shell is computed by the FE-FE-BE method and FEM-BEM respectively, and then the numerical results of the FE-FE-BE method are compared with the experimental data [28] and the FEM-BEM results. The FEM-BEM has been verified by Reference [15]. The parameters of the model can be found in [28] and is shown in Figure 6. The driving force is applied to a point on the middle rib, where F = 0.454cos (2πft) Newton and f = 258 Hz. The material properties of the cylinder are as follows: Young’s modulus E = 2.06 × 1011 N·m−2, Poisson’s ratio ν = 0.3, structural damping ratio 2ξ = 0.06, and density ρs = 7850 kg·m−3. The density of the fluid is ρ = 1000 kg·m−3 and the speed of sound in water is c = 1461 m·s−1.Figure 6
Geometry of the ribbed cylindrical shell and underwater arrangement.In the FEM-BEM, the cylinder shell is modelled with 3600 triangular elements and the ribs are meshed with 280 CQUAD4 elements as shown in Figure7(a). The boundary elements are coincided with 3600 triangular elements on the wet surface. As shown in Figure 6, the 360 acoustic field points are evenly arranged on a circle which is 6096 mm away from the cylinder axis.Figure 7
(a) FE elements of FEM-BEM; (b) FE elements of FE-FE-BE method.
(a)(b)In the FE-FE-BE method, the element size of the wet structural surface is the same as in the FEM-BEM. The artificial boundary surface is a cylindrical surface with lengthL = 2145 mm and radius R = 1000 mm, which has the common axis and centroid with the structure. The artificial boundary surface is modelled with 1628 triangular elements and 816 nodes, where the maximum element size is about λ/28, λ is the acoustic wavelength at 258 Hz. Figure 7(b) shows the FE elements of the FE-FE-BE method, the interior fluid is meshed with 26317 tetrahederal elements fitting the meshes of the wet surface and artificial boundary surface.The numerical results include the effect of reflections from the free water surface. The comparison of numerical results and experimental data is plotted in Figure8. The results of the FE-FE-BE method are in good agreement with the results of FEM-BEM. Both numerical results match well with the experimental data, besides at certain angles (around 60° and 120°) where the experimental data may contain a small error due to the low signal-to-noise ratio. The above comparisons indicate that the FE-FE-BE method is valid.Figure 8
Results of the FE-FE-BE method, FEM-BEM, and experiment.


## 4.1. A Uniformly Driven Spherical Shell

The sound pressure radiated from an underwater spherical shell driven by a uniform pressure unit is calculated, and the results are compared with the theoretical data [13]. The pressure load is p0 = cos (2πft). The structural and fluid characteristics are shown in Table 1, where a is the radius of the spherical shell, t is the thickness of the shell, ν is Poisson’s ratio, E is Young’s modulus, ρs is the density of structure, ρf is the density of water, and c0 is the speed of sound in water.Table 1
Properties of materials.
a (m)t (m)νE (Pa)ρs (kg·m−3)ρf (kg·m−3)c0 (m·s−1)0.50.010.32.07e11766910001524The spherical shell is meshed with 3200 triangular elements and 1602 nodes, as shown in Figure4(a). The artificial boundary is a spherical surface with radius R = 2 m, which is also meshed with 3200 triangular elements as shown in Figure 4(b). The interior fluid is meshed with CTETRA elements between the artificial boundary elements and wet interface elements as shown in Figure 4(b). Figure 4(c) depicts the 360 acoustic field points on a circle with a radius of 100 m.Figure 4
(a) Structural FE meshes; (b) interior fluid FE meshes; and (c) acoustic field points.
(a)(b)(c)The numerical model calculates in the range from 100 Hz to 2000 Hz. Figure5 shows the comparison between the numerical and theoretical results, where the former is the average acoustic pressure on 360 field points. The numerical results agree well with the theoretical results in the calculation range.Figure 5
Results of the acoustic pressure (10−4 Pa).


## 4.2. A Ribbed Cylindrical Shell Driven by Force

The sound pressure radiated from an underwater ribbed cylindrical shell is computed by the FE-FE-BE method and FEM-BEM respectively, and then the numerical results of the FE-FE-BE method are compared with the experimental data [28] and the FEM-BEM results. The FEM-BEM has been verified by Reference [15]. The parameters of the model can be found in [28] and is shown in Figure 6. The driving force is applied to a point on the middle rib, where F = 0.454cos (2πft) Newton and f = 258 Hz. The material properties of the cylinder are as follows: Young’s modulus E = 2.06 × 1011 N·m−2, Poisson’s ratio ν = 0.3, structural damping ratio 2ξ = 0.06, and density ρs = 7850 kg·m−3. The density of the fluid is ρ = 1000 kg·m−3 and the speed of sound in water is c = 1461 m·s−1.Figure 6
Geometry of the ribbed cylindrical shell and underwater arrangement.In the FEM-BEM, the cylinder shell is modelled with 3600 triangular elements and the ribs are meshed with 280 CQUAD4 elements as shown in Figure7(a). The boundary elements are coincided with 3600 triangular elements on the wet surface. As shown in Figure 6, the 360 acoustic field points are evenly arranged on a circle which is 6096 mm away from the cylinder axis.Figure 7
(a) FE elements of FEM-BEM; (b) FE elements of FE-FE-BE method.
(a)(b)In the FE-FE-BE method, the element size of the wet structural surface is the same as in the FEM-BEM. The artificial boundary surface is a cylindrical surface with lengthL = 2145 mm and radius R = 1000 mm, which has the common axis and centroid with the structure. The artificial boundary surface is modelled with 1628 triangular elements and 816 nodes, where the maximum element size is about λ/28, λ is the acoustic wavelength at 258 Hz. Figure 7(b) shows the FE elements of the FE-FE-BE method, the interior fluid is meshed with 26317 tetrahederal elements fitting the meshes of the wet surface and artificial boundary surface.The numerical results include the effect of reflections from the free water surface. The comparison of numerical results and experimental data is plotted in Figure8. The results of the FE-FE-BE method are in good agreement with the results of FEM-BEM. Both numerical results match well with the experimental data, besides at certain angles (around 60° and 120°) where the experimental data may contain a small error due to the low signal-to-noise ratio. The above comparisons indicate that the FE-FE-BE method is valid.Figure 8
Results of the FE-FE-BE method, FEM-BEM, and experiment.


## 5. Comparison between the FE-FE-BE Method and FEM-BEM for Large Complicated Structures

The vibration response and acoustic radiation from a large, complicated structure at low frequencies are calculated by the FE-FE-BE method and FEM-BEM, then their results and the computation time are compared. To simulate the large underwater structure, a simplified cylinder with ribs and bulkheads is modelled. As shown in Figure9, the parameters of the cylinder are as follows: length L = 36 m, diameter d = 3.6 m, length-diameter ratio is 10, the ring-ribs are evenly spaced with a gap of 500 mm, and the cylinder is divided into six compartments by five bulkheads. The one Newton of force is applied to the cylinder middle at frequency f = 100 Hz. The material properties of every part are the same as follows: Young’s modulus E = 210000 MPa, Poisson’s ratio ν = 0.3, shear modulus G = 80769 MPa, structural damping ratio 2ξ = 0.06, and density ρs = 7800 kg·m−3. The thickness of the shell and bulkhead are 30 mm and 15 mm. The ring-rib section is rectangular and the size is shown in Figure 9.Figure 9
The geometry of the cylinder.As shown in Figure10, the wet structural surface is meshed with 31267 triangular elements with three elements on one rib space. To make the number of boundary elements acceptable by the FEM-BEM computation, only three elements are modelled on one rib space but at least five elements are needed in practice. The bulkhead is meshed with CQUAD4 elements and the rib is meshed with bar elements. In Figure 11, the boundary elements are coincided with the wet surface elements and the acoustic field points are on the circle with a radius of 100 m in FEM-BEM.Figure 10
Finite elements of the structure.Figure 11
Boundary elements and acoustic field points in FEM-BEM.In the FE-FE-BE method, the finite elements of the structure are the same as in the FEM-BEM. The artificial boundary surface is a cylindrical shell with a length of 50 m and a diameter of 9 m, and the distance to the wet surface is approximately half the acoustic wavelength. To ensure the boundary element size is less thanλ/6, the artificial boundary surface is meshed into 1392 triangulate elements as shown in Figure 12. The interior fluid is meshed with 312173 tetrahedral elements compatible with the boundary elements and wet surface elements. It can be seen in Figure 13 that the tetrahedral element size increases gradually from the wet surface elements to the boundary elements. The acoustic field points are the same as in FEM-BEM in Figure 12.Figure 12
Artificial boundary elements and acoustic field points in the FE-FE-BE method.Figure 13
Finite elements of the structure and interior fluids in the FE-FE-BE method.To analyze the influencing factors of computation time, five cases are calculated, namely, Case 1, the structural vibration response in vacuum; Case 2, the response of structure only coupled with the interior fluid; Case 3, the acoustic radiation and response of structure coupled with all fluid domains using FEM-BEM; Case 4, the acoustic radiation and response of the structure coupled with all fluid domains using the FE-FE-BE method, whose model is shown in Figures12 and 13; Case 5, the number of artificial boundary elements is increased to 4768 compared with Case 4. Note that Case 2 does not satisfy the boundary conditions in practice and it is only used to obtain the computation time. To ensure comparability of computation time, the computer configurations are given in Table 2.Table 2
The configurations of computer used in this paper.
ItemIndexCPUIntel(R)XEON(R)CPU E5-2690 v2 @3.00 GHzOSWindows 7RAM128 GBFigure14 shows a comparison between the results of the two numerical methods. The polar axis indicates the sound pressure level (Ref. 1 μPa). The results of the FE-FE-FE method are in good agreement with those of FEM-BEM. Since the FEM-BEM has been verified by Reference [13], the results of the FE-FE-BE method are credible enough. The results of Case 4 agree well with those of Case 5. It indicates that the artificial boundary element size less than λ/6 is enough to predict the far-field SPL correctly. Figure 15 shows the vibration response in vacuum. Figure 16 is the response of the structure coupled with all fluid domains using the FE-FE-BE method. It can be seen that the vibration response in vacuum differs significantly from that in fluid. Hence, the effect of the fluid cannot be ignored in underwater structural vibration and acoustic radiation. Comparing Figure 16 with Figure 17, the shape response of the structure in fluid in the FE-FE-BE method is coincident with that in FEM-BEM. The results in Figures 15–17 show that the vibration range of the structure in the fluid is larger than that in vacuum. It may be that the fluid action makes the structure more easily excited at this frequency.Figure 14
The comparison results of the FE-FE-BE method and FEM-BEM.Figure 15
The response of the structure in vacuum (case 1).Figure 16
The response of the structure coupled with the fluid (case 4).Figure 17
The response of the structure coupled with the fluid (case 3).Table3 lists the computation time of all calculation cases. Compared with Case 3, the computation time of Case 1 and Case 2 are extremely short even though there is a huge number of finite elements. The computation time of the FE-FE-BE method decreases over 99% compared with FEM-BEM. For Case 1 and Case 2, only sparse matrices are involved in the FEM. However, it is time-consuming to solve the dense matrices in BEM in Case 3. Hence, the main factor affecting the computation time is the dimension of dense matrices that increases with the number of boundary elements. Although it is difficult to avoid dense matrices in the FEM-BEM and FE-FE-BE methods, as the BE element size in the FE-FE-BE method can model larger than FEM-BEM at low frequencies, the number of boundary elements and computational time decreases significantly.Table 3
Computational time.
CaseMethodNumber of BEDOF of fluid FEDOF of structural FEComputation timeCase 1FE0017101 × 69 sCase 2FE0312173 × 117101 × 633 sCase 3FE-BE31267017101 × 65.5 hCase 4FE-FE-BE1392312173 × 117101 × 679 sCase 5FE-FE-BE4768516909 × 117101 × 66 min 48 sFigure18 shows the interior fluid pressure in the FE-FE-BE method, indicating that the acoustic pressure near the structural surface is complicated and changeable even in a small region. This means that the near-field evanescent waves exert effects and using enough small elements (compared with acoustic wavelength) is necessary. At the same time, the sound pressure in the far field and the sound pressure on the artificial boundary change slowly in space. Therefore, larger-scale units can be used on the artificial boundary. These phenomena and analysis conclusions are consistent with the physical mechanism of sound radiation from large complex structures.Figure 18
Interior fluid acoustic pressure (case 4).


## 6. Conclusions

Based on the physical mechanism of structural acoustics, a new approach combining the multidomain finite element method and artificial boundary element method was proposed to simulate the structure-acoustic interaction for large, complicated structures at low frequencies. Theoretical analysis shows that the effects of the exterior fluid on the interior fluid were the added mass and damping of the exterior fluid, which were added to the interior fluid mass and damping matrices with the DMAP language in NASTRAN. Theoretical analysis shows that FE-FE-BE results were in good agreement with the theoretical results and experimental data. Compared with FEM-BEM, the computation time of the proposed method decreased over 99% in the acoustic radiation calculation from a large, complicated vibrating structure. Furthermore, the study indicates that the unit scale on the artificial boundary can be larger than that on the structural surface of the FE-BE method by using the physical mechanism of the sound radiation from large complex structures while maintaining the same calculation accuracy. The FE-FE-BE method in this paper can be used for rapid prediction of radiation noise of large underwater structures (such as submarines) and further used for radiation noise control and structure optimization. Besides, the method is applicable for the multidomain or multibody structural acoustic calculation, including some equipment with complex surfaces (acoustic vibration analysis of automobile shell, muffler design, etc.).



---

*Source: 3631241-2022-12-31.xml*