A 2D BEM-FEM approach for time harmonic fluid-structure interaction analysis of thin elastic bodies

A 2D BEM-FEM approach for time harmonic fluid-structure interaction analysis of thin elastic bodies∗ J.D.R. Bordón, J.J. Aznárez, O. Maeso Instituto Universitario de Sistemas Inteligentes y Aplicaciones Numéricas en Ingeniería, Universidad de Las Palmas de Gran Canaria, Edificio Central del Parque Científico y Tecnológico del Campus Universitario de Tafira, 35017 Las Palmas de Gran Canaria, Spain {jdrodriguez,jaznarez,omaeso}@iusiani.ulpgc.es Abstract This paper deals with two-dimensional time harmonic fluidstructure interaction problems when the fluid is at rest, and the elastic bodies have small thicknesses. A BEM-FEM numerical approach is used, where the BEM is applied to the fluid, and the structural FEM is applied to the thin elastic bodies. From the fluid point of view, the thin elastic bodies are considered of null thickness. This assumption is treated using simultaneously the Singular Boundary Integral Equation and the Hypersingular Boundary Integral Equation. It is assumed that the thin elastic bodies are under the Euler-Bernoulli hypotheses with added rotational inertia. The BEM equations (fluid) and the FEM equations (thin bodies) are coupled using appropriate equilibrium and compatibility conditions. The developed BEM-FEM model requires a simple discretization and leads to a small number of degrees of freedom, although it has some limitations that are studied in some depth. This approach is validated with existing results in the field of sound barriers, and new results using complex barrier shapes are presented. Also, a parametric study about a straight wall immersed in a fluid is done, which provides results of practical usage. Keywords: SBIE/HBIE dual boundary formulation, BEMFEM coupling, thin bodies, fluid-structure interaction, wave propagation, flexible sound barriers 1 Introduction The Boundary Element Method (BEM) and the Finite Element Method (FEM) can handle problems such as heat conduction, electrostatics, elastostatics and elastodynamics, just to name a few. Nevertheless, each method has its own strengths and weaknesses [4]. The combination of both methods comes up when neither the FEM nor the BEM is adequate to face a problem. This is the case of the Fluid-Structure Interaction (FSI) problem posed here, in which there are thin elastic bodies surrounded by a fluid where wave propagation phenomena take place. The BEM is widely used for time harmonic wave propagation in fluids, viscoelastic solids, poroelastic solids, and when ∗Draft of the paper originally published in Engineering Analysis with Boundary Elements 43 (2014) 19–29 http://dx.doi.org/10.1016/j.enganabound.2014.03.004 regions of any of these types are interacting with each other. When each region is treated by the BEM, the approach is called BEM-BEM. A lot of work has been done about it when applied to dynamic Fluid-Soil-Structure Interaction (FSSI) and its particular cases: FSI and Soil-Structure Interaction (SSI). In this field, the work of Domínguez and co-workers [41, 16, 3] must be highlighted. There are very complete reviews [5, 6, 21]. Also, the BEM-BEM approach for FSI problems in the field of sound barriers has been studied [50]. Taking this into account, the problem posed here can be solved by the BEM-BEM approach. However, when thin elastic solids appear in a problem, some interrelated difficulties emerge in the BEM-BEM approach: discretization needs, quasi-singular integration accuracy, and the Linear System of Equations (LSE) degeneracy. A thin body is characterized by having faces very close each other. Depending on the element size, the relative distance between an element and any node not belonging to it can be very small. This relative distance is the most relevant factor when evaluating BEM quasi-singular integrals. Thus, the maximum size of the elements heavily depends on the quasisingular integration capabilities. Several quasi-singular integration strategies have been developed through the years. It is worth mentioning two simple but powerful classical strategies: adaptative subdivision with selection of the quadrature order [32, 26] and adaptative cubic transformation [52, 53]; a mix of both is used in this work. When body thickness is too small, more elaborated strategies are needed. Among others, Liu and co-workers worked on it in 2D [36] and in 3D [35, 34, 12], and they showed that very thin bodies can be efficiently treated. Recent works contain brief updated reviews of quasi-singular integration strategies [55, 57, 56]. Given an exterior region with a thin body inside it, Krishnasamy et al. [27] showed that if only Singular Boundary Integral Equations (SBIE) are applied to build the BEM final LSE, then its condition number get worse as thickness decreases, becoming completely degenerated if the thickness is null. If the thin body is not surrounded by an exterior region, i.e. the thin body is considered alone, Liu et al. [34] demonstrated that, if the primary variables are not constrained at all boundaries, then the LSE does not degenerate. In both cases, the discretization must be carefully done, and a capable quasisingular strategy is mandatory. If the thin body region is a viscoelastic region, the well known structural beam/shell hypotheses are applicable to 1 model it. Doing so, its dimensional space is reduced to 1D for a beam, and to 2D for a shell. The FEM is appropriate to discretize these structural elements, which reduces heavily the discretization effort and the number of degrees of freedom when compared to the BEM. However, from the point of view of the region that contains the thin body, the FEM discretization is seen as a degenerated geometry, i.e. a null thickness geometry, that can not be directly handled by the conventional BEM. Two ways of solving this difficulty for our problem are: the multiregion approach [7], or employing the Hypersingular BIE (HBIE) in combination with the SBIE [27]. The multiregion approach needs the definition of some artificial boundaries, which can be hard to do. However, the SBIE/HBIE dual boundary formulation is applied directly to the null thickness geometry. The SBIE/HBIE dual boundary formulation emerged to solve fracture mechanics problems [15, 47], but it has been applied to other problems like sound propagation [30, 31, 48, 49]. Its main drawback is handling with the HBIE, which is more difficult to treat than the SBIE. The HBIE has received much attention, particularly about the continuity requirements [29, 39, 40] and regularization techniques [28, 20, 51, 10]. The Cauchy Principal Value (CPV) and Hadamard Finite Part (HFP) definitions are often used to deal with it, see for example [23]. Nevertheless, in line with Guiggiani [20], we prefer making explicit the whole limiting process in order to see clearly how the unbounded terms cancel out, and only regular or weakly singular integrals remain. We propose coupling directly the SBIE/HBIE dual boundary formulation with the structural FEM. The idea of BEM-FEM coupling arose since the BEM beginnings [4]. It is used to overcome difficulties such as nonlinearities, or to reduce the discretization and computational effort, that is our case. An example of its usefulness could be seen in the work of Padrón et al. [42, 43, 44], who coupled FEM beams (piles) with 3D BEM viscoelastic regions (stratified soils). Most works about BEM-FEM applied to FSI deal with closed thin structures like boxes, cylinders, spheres, ships or submarines [8, 54], where a dual boundary formulation is not needed. A much smaller number of works deal with open thin structures, and their methodologies differ from ours. Jean [25] used a variational approach discretized with 2D boundary elements for the fluid and 2D finite elements for the structure. Z.S. Chen et al. [13] used a Symmetric Galerkin BEM for 3D FSI problems. Recently, L.L. Chen et al. [11] used the FEM in combination with the Wideband Fast Multipole Method to handle 2D FSI problems. The proposed direct BEM-FEM approach is presented as follows. The fluid is treated by the BEM. From the fluid region point of view, the structure is considered as a null thickness body, and the simultaneous application of the SBIE and the HBIE is used to handle it. The fluid basic formulation is shown at 2.1, and some details of the HBIE regularization process are given at A. The thin bodies are discretized using structural straight FEM elements with Euler-Bernoulli hypotheses with added rotational inertia, which is shown at 2.2. At 2.3, the BEM equations (fluid) and the FEM equations (thin bodies) are coupled using equilibrium and compatibility conditions. Two limitations exist in this approach: the null thickness assumption of the thin elastic body from the fluid point of view, and the Euler-Bernoulli hypotheses of the thin elastic body; they are studied at 2.4. The proposed BEM-FEM approach is validated at 3.1. In order to demonstrate its potential, complex sound barrier shapes are studied at 3.2, and a parametric study about a straight wall is done at 3.3. 2 Methodology The problem consists in the harmonic analysis of a twodimensional domain composed by a fluid region and many viscoelastic thin regions (thin bodies). Onwards, as usual in the harmonic analysis, the frequency is denoted as f, and the angular frequency is ω=2πf. 2.1 Fluid (BEM) The fluid is considered homogeneous, inviscid, at rest, its body forces are neglected, and the excitations are low enough to admit small disturbances (linear behaviour). As it is well known, under these hypotheses the governing equation is the Helmholtz PDE. Let Ω⊂R2be a fluid region, ˜ ρits density, and ˜ cits wave propagation speed, the pressure pat any point x∈Ωobeys: ∇2p+k2p=0 in Ω(1) where kis the wave number (k=ω/˜ c), and sources are neglected for the sake of brevity. The boundary of the region Ωis denoted as Γ=∂Ω, and the normal vector nis defined outwards. The pressure pacts as the primary variable, while the secondary variable could be the pressure flux qor the displacement in the normal direction un: q=∂p ∂n,un=1 ˜ ρω2 ∂p ∂n(2) The latter is physically more meaningful than the former, and is used when establishing compatibility conditions. However, qis chosen as the secondary variable, which is more common in the literature. 2.1.1 Singular BIE The pressure BIE of (1) applied at a point xi(collocation point) is: cpi+Z Γ q∗pdΓ=Z Γ p∗qdΓ(3) where each term of the equation is: c=   0, xi∉Ω∪Γ 1, xi∈Ω ]0,1[,xi∈Γ p∗=1 2πK0(ikr ) q∗= − ik 2πK1(ikr )∂r ∂n (4) where iis the imaginary unit, r=¯¯x−xi¯¯is the distance between observation and collocation points, and Kn(z)is the 2 modified Bessel function of the second kind, order n, and argument z. Kn(z)properties and expansions can be found in [1, Chapter 9]. When xiis taken to Γ, the integrals in (3) contain a singularity, but they are integrable if a limiting process from inside or outside Ωis followed. The integration domain Γis partitioned as Γ=limǫ→0{(Γ−eǫ i)+Γǫ i}, where eǫ iis the exclusion zone of Γ, and Γǫ iis an arc of radius ǫthat surrounds xi. The integration over Γǫ iproduces the free-term c, which is in the interval ]0,1[, being 1/2 if Γis smooth at the collocation point, i.e. Γ(xi)∈C1. The integrals over Γ−eǫ iare at most weakly singular, as is well known. In this work, the collocation points are placed at smooth boundary points, so, in the following, the equations are written under this hypothesis. If the boundary Γis partitioned in Neelements, Γ= ∪Ne 1Υj, and geometry, p, and qare interpolated over each element Υjusing Lagrange elements, then the discretized SBIE can be written as: 1 2φ˜  i·p˜ + j=Ne X j=1 hj i·pj= j=Ne X j=1 gj i·qj,xi½∈Υ˜  ∉∂Υ˜ (5) where hj iand gj iare the integral kernels of the element jwhen the SBIE is applied at xi. The element ˜ is the one that contains xi, and φ˜  iis the vector of shape functions of the element ˜  evaluated at xi. 2.1.2 Hypersingular BIE In order to obtain the pressure flux BIE, the derivative of the pressure BIE with respect to a direction dis taken: cµ∂p ∂d¶i +Z Γ ∂q∗ ∂dpdΓ=Z Γ ∂p∗ ∂dqdΓ,xi∉Γ(6) where each term of the equation is: c=½0, xi∉Ω∪Γ 1, xi∈Ω ∂p∗ ∂d= − ik 2πK1(ikr )∂r ∂d ∂q∗ ∂d=ik 2π·ikK2(ikr )∂r ∂d ∂r ∂n+1 rK1(ikr )(d·n)¸ (7) When xiis taken to Γ,d=niis used, where niis the normal vector at the collocation point. The left-hand side integral of (6) contains a singularity, which is stronger than the right-hand side integral of (6) and the integrals of (3). Thus, when xiis taken to Γ, the pressure flux BIE is called the Hypersingular BIE. Given a hypersingular integral I=RB AF(x)/(x− xi)2dx,A<xi<B, if Fhas certain continuity properties, then Iexists. Fmust belong to the Hölder function space C1,α[40]. To do so, the pressure must be differentiable at the collocation point, i.e. p¡xi¢∈C1. Similarly to the SBIE, a limiting process where Γ= limǫ→0{(Γ−eǫ i)+Γǫ i} is needed. The integration over Γǫ inot only produces a free-term, but also produces an unbounded term: 1 2µ∂p ∂ni¶i −pi πlim ǫ→0µ1 ǫ¶+lim ǫ→0Z Γ−eǫ i ∂q∗ ∂ni pdΓ=lim ǫ→0Z Γ−eǫ i ∂p∗ ∂ni qdΓ w Γ Ω Γ Ω w→0 Figure 1: Approximation when using the null thickness assumption (8) In the following, the right-hand and the left-hand side integrals of (8) are called Land M, respectively. Lis regular since has the same kind of singularity as the left-hand side integral of (3). Taking into account that K1(ikr )=1/(ikr )+KR 1(ikr ), where KR 1(ikr )=O(rlnr), Mcan be decomposed into: M=(ik)2 2πlim ǫ→0Z Γ−eǫ i K2(ikr )∂r ∂ni ∂r ∂npdΓ+ +1 2πlim ǫ→0Z Γ−eǫ i 1 r2¡ni·n¢pdΓ +ik 2πlim ǫ→0Z Γ−eǫ i 1 rKR 1(ikr )¡ni·n¢pdΓ=M1+M2+M3 (9) where M1is regular, M2is hypersingular and M3is weakly singular. A regularization process is required for M2(see A). It is based on Sáez et al. contributions [45], who applied it to elastostatics. Unlike Sáez et al., M2is regularized before discretization, which gives some interesting insights into this integral. Through the regularization process of M2emerges an unbounded term that cancels out the one that appears in (8), which leads to the regularized HBIE. As did with the SBIE, the discretized HBIE can be written as: j=Ne X j=1 mj i·pj= − 1 2φ˜  i·q˜ + j=Ne X j=1 lj i·qj,xi½∈Υ˜  ∉∂Υ˜ (10) where mj iand lj iare the integral kernels of the element jwhen the HBIE is applied at xiwith a normal ni. Only the integral kernel vector m˜  iof the element ˜ uses the regularized M. 2.1.3 BIEs for coincident boundaries When nearly coincident boundaries (nearly coplanar boundaries) belong to the same region (i.e. a crack, a thin void, a thin inclusion, or a thin scatterer) the LSE is nearly-singular. In the limit when boundaries are coincident (null thickness discontinuity), the LSE becomes singular. The simultaneous application of the SBIE and the HBIE can solve this difficulty [27]. The null thickness assumption is very interesting because it can greatly reduce the discretization and the computational effort at the expense of an approximation of the field around the discontinuity (see Figure 1). The computational cost reduction can be >60% [30, Table 1], depending on the problem, the analysed frequencies, and the implementation. 3 Let Γbe the boundary of a region Ω, resulting from the approaching of two identical boundaries, whose normals are pointing at each other, until they are coincident. One of the boundaries is chosen as the positive face Γ+of Γ, which is used as the reference face for the whole Γ. Thus the normal vectors nand niare defined on it. Each face has two variables, the pressure and the pressure flux,so there are four variables at the collocation point i:p+ i,q+ i,p− iand q− i. The limiting process can be done using the integration domain depicted in Figure 2: Γ=lim ǫ→0nhΓ+−¡eǫ i¢+i+¡Γǫ i¢++£Γ−−¡eǫ i¢−¤+¡Γǫ i¢−o(11) The singularity is avoided twice: when integrating over Γ+, and when integrating over Γ−. The resulting SBIE is: 1 2p+ i+lim ǫ→0Z Γ+−³eǫ i´+ q∗pdΓi+1 2p− i+lim ǫ→0Z Γ−−³eǫ i´− q∗pdΓi= =lim ǫ→0Z Γ+−³eǫ i´+ p∗qdΓ+lim ǫ→0Z Γ−−³eǫ i´− p∗qdΓ,xi∈Γ (12) and the resulting HBIE is: 1 2q+ i−p+ i πlim ǫ→0µ1 ǫ¶+lim ǫ→0Z Γ+−³eǫ i´+ ∂q∗ ∂ni pdΓ− −1 2q− i−p− i πlim ǫ→0µ1 ǫ¶+lim ǫ→0Z Γ−−³eǫ i´− ∂q∗ ∂ni pdΓ= =lim ǫ→0Z Γ+−³eǫ i´+ ∂p∗ ∂ni qdΓ+lim ǫ→0Z Γ−−³eǫ i´− ∂p∗ ∂ni qdΓ,xi∈Γ (13) where q+ i=(∂p+/∂ni)iand q− i= −(∂p−/∂ni)i. Similarly to 2.1.2, two unbounded terms have been produced when solving the integrals over (Γǫ i)+and (Γǫ i)−. Likewise, the regularization process of the left hand side integrals of (13) produces two unbounded terms that cancel out the previous ones. The discretization process is similar to that followed for the SBIE and the HBIE, with the condition that the faces Γ+and Γ− have the same discretization. Because of that, it is possible to build a new type of element Υjthat is composed by two subelements Υ+ jand Υ− j. The variables of Υjcan be written as: pj=n³pj´+³pj´−o,qj=n³qj´+³qj´−o(14) Taking advantage of n=n+= −n−, the integral kernels of Υj can be written only in terms of the integral kernels of the subelement Υ+ j. The discretized SBIE and the discretized HBIE for a collocation point xibelonging to coincident boundaries are: 1 2nφ˜  iφ˜  io·p˜ + j=Ne X j=1 hj i·pj= j=Ne X j=1 gj i·qj(15) Ω Γ+ Ω Γ− w→0 n=n+ ni=n+ i −ni=n− i ³eǫ i´+ ³eǫ i´− −n=n− face + face − ³Γǫ i´+ ³Γǫ i´− Figure 2: Singularity treatment for coincident boundaries j=Ne X j=1 mj i·pj=1 2n−φ˜  iφ˜  io·q˜ + j=Ne X j=1 lj i·qj(16) If the fluid is uncoupled, then these equations are handled in the usual way, but having Γ+and Γ−independent boundary conditions. If the fluid is coupled with a thin elastic body, which is the main input of this paper, these equations together with those presented in 2.2 and 2.3 are used. 2.1.4 Discretization and collocation procedure The discretized equations written above have been developed under Γ¡xi¢∈C1and p¡xi¢∈C1hypotheses at the collocation point. By doing so, the collocation procedure described here can be applied simultaneously to the HBIE and the SBIE at coincident boundaries, giving a uniform approach to build the BEM equations. The C1requirement can be fulfilled by many ways, among others: cubic splines [33], an interpolation algorithm [19], Overhauser elements [9] or discontinuous Lagrange elements [45]. The way we deal with it is using continuous isoparametric Lagrange elements with non-nodal collocation at vertex nodes, and adding up the BIEs associated with each vertex node. This strategy is known as the Multiple Collocation Approach (MCA), and it was introduced by Gallego et al. [18, 17, 2]. It is very simple, gives accurate results, and makes BEM-FEM coupling relatively easy. In this work, quadratic elements are used. Given a vertex node iand its elements Υ1and Υ2, two SBIEs are added up to build the equation associated with the node: one SBIE is collocated inside the element 1 at x1 i′, and the other SBIE is collocated inside the element 2 at x2 i′; the same is done with HBIEs (see Figure 3). Let δbe defined as the displacement of the collocation point towards the inside of the element. Given an element with a local system of coordinates −1≤ξ≤1, if ξiis the nodal position of the node i, then the local coordinate ξi′of the displaced collocation point is: ξi′=ξi(1−δ),0 <δ<1⇒xi′=x¡ξi′¢(17) 4 ix1 i′ x2 i′Υ2Υ1 ξ1 i ξ2 i ξ1 i′=ξ1 i(1−δ) ξ2 i′=ξ2 i(1−δ) Ω Γ n Collocation pointNode Vertex node + ++ + + Figure 3: Multiple Collocation Approach Note that δ=0 gives a collocation point at the nodal position, while δ=1 gives a collocation point at the element centre. A question that arises is how much the collocation point should be displaced from the nodal position. Ariza et al. [2] used a value of δ=0.25, although it was not explained why. To the authors’ knowledge, there is no published work about an optimum value of δfor the MCA. Nevertheless, the MCA can be related to discontinuous elements because the set-up of the collocation points is the same. Marburg [38] studied the optimum position of nodes of discontinuous elements for a sound propagation problem. He found that nodes located at the zeros of the Legendre polynomials gives optimum results. He also stated that the hypersingular formulation may have other optimal locations. Thus, we use δ=0.2254. 2.2 Thin elastic bodies (FEM) In a plane deformation problem, a thin body has infinite width along x3, and finite thickness wand length Lin the x1−x2 plane. Under these conditions, the thin body could be considered as a beam with a cross section A=w·1, a length L, and a modified Young’s modulus E=Em/(1−ν2), where Emis the Young’s modulus of the material and νits Poisson’s ratio. Onwards, when the term “beam” is used, it must be understood this way. Let Ωsbe a thin elastic body. It can be split into straight beam FEM elements Υj, which are under the Euler-Bernoulli hypotheses with added rotational inertia. In order to have a node-to-node correspondence with a quadratic BEM element, a beam FEM element with three nodes and eight degrees of freedom is considered [42] (see Figure 4). The vertex nodes i=1,2 have translation u(i) 1,u(i) 2and rotation θ(i), while the central node i=3 has only translation u(3) 1,u(3) 2. Each element has a density ρ, a modified elastic modulus E, a thickness w, an inertia I=(1/12)·w3·1, and a length L. Damping of hysteretic type is introduced by defining a complex Young’s modulus E=Re(E)(1+i2ξ), where ξis the damping coefficient. Because axial behaviour and lateral behaviour are decoupled, axial and lateral elemental matrices can be obtained separately in the local system of coordinates. From now on, variables carrying an apostrophe are variables expressed in the local system of coordinates. The axial displacement u′ 1is interpolated using a Lagrange 231 u′ 2 (1) u′ 1 (1) u′ 1 (3) u′ 1 (2) θ(2) θ(1) u′ 2 (3) L 2 L 2 x′ 1 x′ 2 u′ 2 (2) s′ 2 (1) s′ 2 (3) s′ 2 (2) Figure 4: Three nodes / eight degrees of freedom FEM beam element quadratic element: u′ 1(ξ)=©φ1φ2φ3ª·nu′ 1 (1) u′ 1 (2) u′ 1 (3) oT=φT·u′a (18) The axial stiffness matrix and the axial translation mass matrix are obtained by using the Principle of Virtual Displacements, respectively: K′ i j a=2 LE AZ1 −1 dφi dξ dφj dξdξ M′ i j ta =L 2ρAZ1 −1φiφjdξ (19) The lateral displacement u′ 2is taken as a fourth degree polynomial in −1≤ξ≤1. The lateral displacement and the rotation are: u′ 2(ξ)=ϕT·u′l,θ(ξ)=ϑT·u′l(20) where: u′l=nu′ 2 (1) θ(1) u′ 2 (2) θ(2) u′ 2 (3) oT(21) ϕ=                    ϕ1 ϕ2 ϕ3 ϕ4 ϕ5                    =                    1 4ξ¡−3+4ξ+ξ2−2ξ3¢ L 8ξ¡−1+ξ+ξ2−ξ3¢ 1 4ξ¡3+4ξ−ξ2−2ξ3¢ L 8ξ¡−1−ξ+ξ2+ξ3¢ 1−2ξ2+ξ4                    ,ϑ=2 L dϕ dξ(22) The lateral stiffness matrix, the lateral translation mass matrix and the lateral rotation mass matrix are obtained by using the Principle of Virtual Displacements, respectively: K′ i j l=µ2 L¶3 EI Z1 −1 d2ϕi dξ2 d2ϕj dξ2dξ M′ i j tl =L 2ρAZ1 −1ϕiϕjdξ M′ i j r=L 2ρIZ1 −1ϑiϑjdξ (23) The lateral distributed load s′ 2along the beam is interpolated using a Lagrange quadratic element: s′ 2(ξ)=©φ1φ2φ3ª·ns′ 2 (1) s′ 2 (2) s′ 2 (3) oT=φT·s′ 2 5 (24) The load s′ 2can be transformed into equivalent nodal forces and moments by using the Principle of Virtual Work: S′ i j l=L 2Z1 −1ϑiφjdξ,i=1,...,5, j=1,...,3 (25) The axial and lateral kinematic variables can be gathered together in u′, and the lateral distributed load vector can be reordered as s′: u′=nu′ 1 (1) u′ 2 (1) θ(1) u′ 1 (2) u′ 2 (2) θ(2) u′ 1 (3) u′ 2 (3) oT s′=n∅s′ 2 (1) ∅ ∅ s′ 2 (2) ∅ ∅ s′ 2 (3) oT (26) so that a stiffness matrix K′is obtained by combining K′aand K′l, a mass matrix M′is obtained by combining M′ta and M′tl +M′r, and a distributed load matrix S′is obtained by reordering the matrix S′l. In the harmonic regime, the dynamic equilibrium equation in global coordinates for a given element is: hL·hK′−ω2M′i·LTi·u=£L·S′¤·s′ Kh·u=Q·s′(27) where Lis the coordinate transformation matrix of the element. This FEM equation is assembled considering all vertex nodes as rigid joints. It must be noticed that s′(distributed lateral loads) is unknown when coupled with the fluid. The element matrices can be easily obtained from (19), (23) and (25), or seen in [42] (except M′r). 2.3 Fluid-structure coupling (BEM-FEM) Once the fluid equations (BEM) and the thin elastic bodies equations (FEM) have been posed, it is possible to combine both by using coupling equations. Let Υjbe a BEM-FEM fluidstructure element composed by three sub-elements: Υ+ j,Υ− j and Υs j; being Υ+ jand Υ− jthe sub-elements associated with both faces of the coincident boundaries of the fluid, and Υs j the sub-element associated with the thin elastic body (see Figure 5). Since the fluid is inviscid, it interacts only laterally with the thin elastic body. Therefore, only lateral compatibility and equilibrium have to be established. The normal displacements of the fluid at the boundary (2) must coincide with the beam lateral displacements u′ 2, for a given node i: u′ 2 (i)= − 1 ˜ ρω2q+ i,u′ 2 (i)=1 ˜ ρω2q− i(28) that leads to q+ i= −q− i. In (27), the displacements are expressed in the global system of coordinates, so their projection onto the x′ 2axis gives the lateral displacements: u(i)·x′ 2= − 1 ˜ ρω2q+ i,u(i)·x′ 2=1 ˜ ρω2q− i(29) Ωf Γ+ Ωf Γ− 1 1 1 3 3 3 2 2 2 Ωs p+ 1 p− 1 p+ 3 p− 3 p+ 2 p− 2 q+ 1 q− 1 u′ 2 (1) q+ 3 q− 3 u′ 2 (3) q+ 2 q− 2 u′ 2 (2) Υ+ j Υ− j Υs ju′ 1 (1) u′ 1 (3) u′ 1 (2) θ(2) θ(1) Figure 5: Coupling between sub-elements Υ+ j,Υ− jand Υs j(local numbering) which are the compatibility equations for a node i. Note that these equations relate a primary variable of the structure (u) with a secondary variable of the fluid (pressure flux q). Thus, if the node iis a vertex node shared by two non-collinear elements, then the pressure flux qis undefined there (corner problem), being defined only just before and after the vertex. If both elements are almost collinear, and we are not interested in local effects, then it is acceptable assuming that the pressure flux is continuous. For this case, (29) is posed for each element and added up to build a unique compatibility condition. If both elements are far from collinear, the BEM variables of the vertex node are doubled, so that two sets of compatibility equations like (29) are posed. The pressure difference between both faces is equal to the lateral distributed load at each node of the beam: s′ 2 (i)=p− i−p+ i(30) which is the equilibrium equation for a node i. For each vertex node there are eight variables: p+ i,q+ i,p− i, q− i,u(i) 1,u(i) 2,θ(i)and s′ 2 (i); and eight equations: SBIE for node i(15), HBIE for node i(16), 2 compatibility equations (29), 1 equilibrium (30), and 3 FEM equations from (27). For each central node, the situation is similar to the vertex node, except that the rotation and its associated FEM equation does not exist. Although the number of unknowns are equal to the number of equations, more conditions are required. The structure needs the necessary kinematic boundary conditions in order to avoid any rigid body motion. The number of unknowns and equations for each node can be easily reduced from 8 to 6 (7 to 5 for the central node). It can be done by substituting (30) in (27), by using only the first equation of (29), and by substituting q+= −q−in (15) and (16) for every node of a BEM-FEM element. This reduction considerably decreases the computational effort. 6 L/w=1000 L/w=500 L/w=200 L/w=100 L/w=50 L/w=20 L/w=10 a0 Average relative error (front face) [%] 1001010.10.01 10 1 0.1 0.01 0.001 0.0001 Figure 6: Average relative error of the null thickness assumption 2.4 Limitations There are two relevant limitations in the proposed model: the null thickness assumption of the thin elastic body from the fluid point of view, and the Euler-Bernoulli hypotheses of the thin elastic body. For practical reasons, it is necessary establishing a validity range. Since studying the limitations using the complete FSI model needs many parameters, it seems to be more efficient studying each limitation in an uncoupled way. The null thickness assumption can be studied considering the thin body as a rigid obstacle. Lacerda et al. [30] worked about this problem in the sound barriers field. They made a study comparing results from real geometries and their null thickness geometries at certain points and frequencies. It is interesting to expand and generalize this topic by using a dimensionless problem. The experiment consists of a rectangular obstacle of length Land thickness wwithin a fluid ( ˜ ρ,˜ c), where a plane wave is propagating perpendicularly to the length with an angular frequency ω. The dimensionless frequency a0=(ωL)/˜ cis used. Eight cases are solved: seven different geometrical slendernesses L/w={10, 20, 50, 100, 200, 500, 1000}, and the case with null thickness. A conservative discretization of six quadratic elements per wavelength is used. Figure 6 shows the average relative error of the pressure field over the front face of the obstacle versus the dimensionless frequency. It can be seen that the error decreases as geometrical slenderness increases, which is an obvious result. For a0<2, the error decreases as a0decreases, being possible to define a frequency limit which ensures an error level. The maximum average error occurs around a0=2, being: 10% for L/w=10, 2% for L/w=100 and 0.3% for L/w=1000. For a0>2, the error slowly increases if L/w>200, and slowly decreases if L/w<200. Figure 7 shows the relative error of the pressure at some selected nodes of the front face versus the dimensionless frequency, for L/w=100. For a0<1 the error is approximately the same at all points. For a0>1, the error near the tip is at 1/100 from the tip at 1/4 from the tip at center a0 Relative error (points at the front face) [%] 1001010.10.01 10 1 0.1 0.01 0.001 0.0001 Figure 7: Relative error of the null thickness assumption (L/w=100) around twice the error at points far from the tip. This behaviour also occurs for other slendernesses. The Euler-Bernoulli hypotheses can be studied ignoring the fluid. A remarkable paper by Han et al. [22] studies the most widespread beam theories in dynamics, including our EulerBernoulli with added rotational inertia (called Rayleigh theory in that paper). Based on the study, Han et al. recommend using the Euler-Bernoulli theory when L/w>29. Nevertheless, from [22, Figure 22], where the first natural frequency versus the mechanical slenderness is studied, it can be seen that the Euler-Bernoulli theory is appropriate even when L/w>10. Thus, the studied limitations have compatible validity ranges. The methodology is valid for geometrical slendernesses greater than 10. However, it must be taken into account that the error produced by the approximations depends on the dimensionless frequency and geometrical slenderness. 3 Results and discussion 3.1 Validation The proposed model has been validated with results published by Jean [25], where a simple noise barrier problem is studied. The problem description is outlined in Figure 8. The fluid Ωf is air with ˜ ρ=1.3 kg/m3and ˜ c=340 m/s. The thin elastic body Ωsis a simple noise barrier 3 m high and 0.01 m thick, and it is clamped to the ground. Three different materials are considered for the barrier Ωf(see Table 1). The ground is a perfectly reflecting surface, i.e. the fluid displacement at the ground is null. A point source located at xs=(−2.3,0.5)is used. The point source is easily added to the BEM equations, as shown in [37]. A comparison between results from Jean [25] and results from the proposed model is shown in Figure 9. The results from [25] are shown as a coloured background image from the original paper. The figure shows three graphs, one for each material. The yaxis of each plot is the difference between pressures absolute values at a point xwhen using a rigid barrier 7 xs=(−2.3,0.5) 0.01 m 3.00 m x1 x2 Reflecting ground Real body Null thickness assumption Point source Ωf: air Ωs: wood, glass, paraglass Figure 8: Noise barrier problem studied by Jean [25] (thickness not to scale) Ωsρ£kg/m3¤E£MN/m2¤ν ξ Wood 650 12.0 0.01 0.0100 Glass 2400 87.0 0.24 0.0005 Paraglass 1190 3.3 0.40 0.0150 Table 1: Materials for the barrier considered by Jean [25] (q=0) and when using a flexible barrier. The natural frequencies fnof each case are plotted as vertical lines, and they are calculated using the cantilever beam equations [14]. The model used in [25] takes into account the real geometry of the barrier, while the proposed model uses a null thickness barrier. The slenderness is L/w=333, so from the barrier behaviour point of view, the Euler-Bernoulli hypotheses are valid. From the fluid behaviour point of view, the null thickness assumption is also valid, see section 2.4. Thus, the proposed model should be able to reproduce the results from [25]. Figure 9 shows excellent agreement between Jean’s model and the proposed model. Peak frequencies and amplitudes are very well reproduced, although some small discrepancies appear in the wood case at frequencies around 850 Hz. 3.2 Complex sound barrier shapes Jean [25] made a broad study comparing results between flexible and rigid simple sound barriers when varying material, thickness, damping coefficient, receiver and source position, and barrier height. In this section, the proposed model is used to study some complex barrier shapes. The layout of the numerical experiments is depicted in Figure 10. Two simple screen barriers (simple barrier and double simple barrier) together with three multi-edge barrier shapes (Y barrier, U barrier and E barrier) are considered. For each shape, all materials from the Table 1 are used, the thickness is w=0.01 m for all pieces, and the effective height is 3 m. The point source is located at ground level and 10 m ahead the barrier [46, 24, 37]. A grid of 3×11 receivers covering 6 ×60 m2is x1 x2 Point source Receivers area 6 m 60 m10 m Barrier 10 m Figure 10: Layout for studying complex sound barrier shapes Simple barrier Y barrier Double simple barrier U barrier E barrier SIL [dB] WoodParaglassGlassRigid 14.0 13.5 13.0 12.5 12.0 11.5 11.0 10.5 10.0 9.5 Figure 11: SIL for different barrier shapes and materials considered. A thousand frequencies uniformly distributed in log10(f) space from fmin =20 Hz to fmax =4000 Hz are used. Instead of taking the pressure as the variable of interest, the Insertion Loss IL is used [37]. The IL is the difference between pressures (in dB) when there is no barrier and when the barrier is placed, so it measures the efficiency of the barrier. We also consider the average Spectral Insertion Loss SIL, which is simply the average IL in the spectrum, leading to a frequencyindependent indicator. The IL and the SIL are averaged values over all receivers. In the literature, it is often assumed that noise barriers are rigid, so it is interesting finding when this hypothesis is valid or not. A first step is using the SIL, Figure 11 shows the SIL for all considered barrier shapes and materials, including the rigid case. It is seen that the rigid case is not conservative when using the SIL as an indicator. However, the maximum difference between the rigid case and any case is below 2 dB, being 1 dB for the simple barrier and double simple barrier, and 2 dB for the Y barrier. Thus, when a global indicator such as the SIL is going to be studied, the rigid assumption seems to be valid. As Jean [25] showed for the simple barrier, when considering the elastic nature of the barrier there is a widespread pressure increment at low frequencies. Although this behaviour seems reasonable, it is interesting to analyse what happens when barriers more complex than the simple one are used. Figure 12 shows the IL spectrum for all studied barrier shapes and materials, including the rigid case. For low frequencies (f<200Hz) appreciable differences be8 Barrier fn At x=(20,2) At x=(5,2) Figure 7 [38] Paraglass f[Hz] ¯¯p¯¯paraglass −¯¯p¯¯rigid [dB] 10009008007006005004003002001000 20 10 0 -10 -20 Barrier fn At x=(20,2) At x=(5,2) Figure 6 [38] Glass ¯¯p¯¯glass −¯¯p¯¯rigid [dB] 20 10 0 -10 -20 Barrier fn At x=(20,2) At x=(5,2) Figure 5 [38] Wood ¯¯p¯¯wood −¯¯p¯¯rigid [dB] 20 10 0 -10 -20 Figure 9: Comparison between results from Jean [25] and results from the proposed model tween rigid and flexible barriers are obtained. The simple barrier behaves as Jean described, with increments of pressure below 5dB, i.e. IL decrements below 5dB. The other barrier shapes have IL decrements below 10dB. For very low frequencies (f<80Hz) there is virtually no noise attenuation. The considered complex barrier shapes strongly influence the IL spectrum, especially at low frequencies. For mid-high frequencies (f>500Hz) the IL spectrum is very similar to a rigid barrier. For simple and double simple barriers, the differences are very small. For Y, U and E barriers, the differences are more noticeable, reaching up to 5dB at some frequencies. Nevertheless, these differences seem to be irrelevant for noise propagation problems. The human ear is less sensitive at low frequencies than at high frequencies, so, at first, this behaviour at low frequencies could be neglected. However, high frequencies are attenuated by losses in the air and on the absorbing surfaces, while low frequencies are not. Furthermore, when a building with windows closed is near the noise barrier, low frequency noises may be amplified inside the building. Therefore, depending on the context, the elasticity of a barrier similar to those studied should be considered. 3.3 Parametric study about a straight wall In order to provide results of practical usage from this BEMFEM approach, a simple but useful problem is studied. The problem consists of a straight wall (beam) (2L,w,ρ,Em,ν,ξ) with its centre clamped, surrounded by a fluid ( ˜ ρ,˜ c), where a pressure plane wave is propagating with unity amplitude, perpendicular direction, and angular frequency ω, see Figure 13. The problem parameters can be reduced to six dimensionless ones: • Wave propagation speeds ratio: ˜ c/c, where c=pEm/ρis the beam axial wave propagation speed. • Densities ratio: ˜ ρ/ρ. • Geometrical slenderness: L/w. • Dimensionless frequency: a0=(ωL)/˜ c. • Damping coefficient: ξ • Poisson’s ratio: ν Table 2 shows the studied values of the dimensionless parameters. The wave propagation speeds ratio and the densities ratio have ranges that include the most extreme fluidstructure combinations. The geometrical slenderness starts from L/w=10 to L/w=1000, which are within the validity interval. The dimensionless frequency range has been chosen so that at least the first natural frequency is clearly captured in all cases. This parametric study is oriented to know the FSI coupling degree. It seems obvious that a decoupled model could be used for extreme cases, e.g. a thick steel wall in air. In these extreme cases, the pressure field in the air is calculated considering a rigid obstacle, and if needed, the pressure field can be 9