A 2.5D spectral approach based on the BEM and the FEM to study wave propagation in fluid acoustics and elastodynamics

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Máster Máster Universitario en Ingeniería de Caminos, Canales y Puertos A 2.5D spectral approach based on the BEM and the FEM to study wave propagation in fluid acoustics and elastodynamics Autor: Francisco Javier Cruz Muñoz Tutor: Antonio Romero Ordóñez Tutor: Pedro Galvín Barrera Dep. Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016 Trabajo Fin de Máster Máster Universitario en Ingeniería de Caminos, Canales y Puertos A 2.5D spectral approach based on the BEM and the FEM to study wave propagation in fluid acoustics and elastodynamics Author: Francisco Javier Cruz Muñoz Advisors: Antonio Romero Ordóñez Pedro Galvín Barrera Dep. Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016 Trabajo Fin de Máster: A 2.5D spectral approach based on the BEM and the FEM to study wave propagation in fluid acoustics and elastodynamics Autor: Francisco Javier Cruz Muñoz Tutor: Antonio Romero Ordóñez Pedro Galvín Barrera El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos He aquí mis mas sinceros agradecimientos a todas las personas sin las cuales no habría podido llegar hasta aquí. A Pedro Galvín y Antonio Romero, por haberme dado la oportunidad de aprender tanto durante el tiempo que hemos compartido y por su inestimable ayuda estos años. A Cristina Lugo, por ser como es y por estar ahí siempre. A mi familia y amigos, por el apoyo mostrado a lo largo de la carrera. A mis compañeros, en breve colegas de la profesión, por haberme hecho el camino hacia la meta más sencillo. Al resto de profesores de la Escuela Técnica Superior de Ingeniería que han participado de mi formación y propiciado que haya llegado hasta aquí. Al proyecto BIA2013-43085-P financiado por el Ministerio de Economía y Competitividad en el que está enmarcado este trabajo, y al Centro Informático Cient´ fico de Andalucia (CICA). I Resumen En este trabajo se presenta una formulación espectral en dos dimensiones y media (2.5D) basada en el método de los elementos finitos (MEF) y el método de los elementos de contorno (MEC) para estudiar la propagación de ondas en medios acústicos y elásticos en problemas tridimensionales (3D). El análisis se lleva a cabo mediante la superposición de problemas bidimensionales (2D) para diferentes números de onda a lo largo de la coordenada longitudinal. El método numérico está basado en la descomposición del dominio para representar el sistema fluido-estructura. Se presenta en este trabajo un elemento finito espectral para representar ondas guiadas en sólidos con sección transversal arbitraria. Además, el MEC se extiende a la formulación espectral para estudiar medios fluidos infinitos. Ambas aproximaciones utilizan polinomios de interpolación de Lagrage como funciones de forma en los puntos de Legendre-Gauss-Lobatto (LGL). Las técnicas propuestas se han verificado mediante dos problemas de referencia: la propagación de ondas en un medio fluido infinito en presencia de una cavidad rígida fija y la propagación de ondas en un carril libre, y los resultados obtenidos presentan un buen grado de acuerdo con las soluciones de referencia. Se realiza un breve análisis h−ppara evaluar la precisión de los métodos. III 2Chapter 1. Introduction likely to dominate the total numerical error, p -FEM elements proved to efficiently mitigate the pollution effect also for convected applications. Many authors have successfully applied the spectral element method in different fields. Degrande and De Roeck [ 7 ] developed a high-order method for studying the dynamic response of a porous saturated medium using spectral elements. They use as shape functions, in a displacement formulation, the analytical solution for Biot equations restricted to two-dimensional (2D) plane strain wave propagation. Bar-Yoseph et al. [ 3 ] presented a space-time spectral element method to solve a simply supported modified Euler-Bernoulli beam undergoing external forced lateral vibrations. They examined two Galerkin formulations using Hermitian and Lagrangian polynomials as interpolants for spatial and time discretisations. The space-time SEM allows the use of larger time steps while maintaining the accuracy of the solution. Kudela et al. [ 15 ] applied spectral finite elements for studying the wave propagation in an isotropic rod and in a Timoshenko beam. The authors proposed this methodology to detect small damage in structures. The computed results were compared to those obtained from the classical FEM and from measurements. This comparison highlighted the efficiency of the spectral elements and it confirms that the methodology reflects this phenomenon quite well. Zhu et al. [ 24 ] coupled Chebyshev spectral elements with the implicit Newmark time integral method to simulate acoustic field. The presented results indicated that the accuracy increased as the time step decreased. The result of this study also concluded that uniform distribution of node is better for the improvement of numerical accuracy. 1.2 Objetives The method proposed in this work regards with a two-and-a-half dimensional (2.5D) approach to represent fluid-solid interaction (FSI) problems. The study of three-dimensional (3D) scattered waves in fluid and solid media requires the use of highly computational demanding models. The solution becomes much simpler if the medium is homogeneous in one direction, even if the dynamic source remains 3D. Such situation is referred to as a two-and-a-half dimensional problems [ 6 ]. This work proposes a 2.5D spectral formulation based on the FEM and the BEM to study wave propagation in fluid-structure interaction problems. The proposed approach is useful for problems where the material and geometric properties are uniform along one direction, and the source exhibits 3D behaviour. 1.3 Original contributions The following points give an outline of what are perceived to be the significant contributions from the present work. • The boundary element method in 2.5D in fluid acoustics has been extended to its spectral formulation to allow for the computation of wave propagation problems in fluid media. This introduces a new accurate method, with lower computational efforts, to solve problems at a higher frequency range. •Moreover, a new 2.5D spectral finite element is presented to study solid waveguides. 1.4 Organization of the text The outline of this work is as follows. First, the numerical model is presented. The spectral finite elements formulation in elastodynamics is described. Then, the spectral boundary elementh method in fluid acoustics is discused. Finally, the proposed methods are then verified with two benchmark problems [20, 10]. 2 Numerical model The spectral formulation is based on a 2.5D coupled BEM-FEM formulation, where the solid ( Ωs ) is represented by the FEM, while the limiting interface ( Γf ) between the solid and the fluid ( Ωf∞ ) is modelled with the BEM (Figure 2.1). In this work both formulations have been developed separately and the coupling procedure is a line of a future research. The 2.5D formulation computes the problem solution as the superposition of two-dimensional (2D) problems with a different longitudinal wavenumber, kz , in the z direction. An inverse Fourier transform is used to compute the 3D solution: a(x,ω) = Z+∞ −∞b a(bx,kz,ω)e−ikzzdkz(2.1) where a(x,ω) is the unknown magnitude (e.g., displacement or pressure), b a(bx,kz,ω) is its representation in the frequency-wavenumber domain, bx=x(x,y,0) , ω is the angular frequency, and i=√−1 . The integral defined by Equation (2.1) is transformed into a summation as follows: a(x,ω) = +∞ ∑ −∞b a(bx,kz,ω)e−ikzz∆kz(2.2) with kz=m∆kz. This equation converges and can be aproximated by a finite sum of terms. 2.1 The 2.5D spectral boundary element formulation The boundary element formulation presented in this work considers an arbitrary boundary submerged in an unbounded fluid medium. The integral representation of the pressure pi for a point i located at the fluid subdomain Ωf∞, with zero body forces and zero initial conditions may be written as [9]: cipi(xi,ω) = ZΓf pi∗(x,ω;xi)ui(x,ω)dΓ−ZΓf ui∗(x,ω;xi)pi(x,ω)dΓ(2.3) Ωs Ωf∞ Γf x y z Figure 2.1 Domain decomposition. 3 4Chapter 2. Numerical model where ui(x,ω) and pi(x,ω) are respectively the normal displacement to boundary Γf and the nodal pressure. ui∗(x,ω;xi) and pi∗(x,ω,xi) are respectively the fluid full-space fundamental solution for normal displacement and pressure at point x due to a point load at xi . The integral-free term ci depends only on the boundary geometry at point i . The integration boundary Γf represents the boundary between the unbounded fluid medium (Ωf∞) and the solid subdomain (Ωs). Assuming that the unbounded medium is constant in the longitudinal direction z , Equation (2.3) is expressed in terms of integrals in this direction and over the cross-section boundary, Σf: cipi(x,ω) = Z+∞ −∞ZΣf pi∗(x,ω;xi)ui(x,ω)dSdz −Z+∞ −∞ZΣf ui∗(x,ω;xi)pi(x,ω)dSdz (2.4) Equation (2.4) is then transformed to the wavenumber domain as: cib pi(bx,ω,kz) = ZΣfb pi∗(bx,ω,kz;bxi)b ui(bx,ω,kz)dS −ZΣfb ui∗(bx,ω,kz;bxi)b pi(bx,ω,kz)dS (2.5) where a hat above a variable denotes its representation in the frequency-wavenumber domain. The problem is discretised into elements, leading to a boundary approximation of the normal displacement and pressure using the interpolation shape functions φ φ φj. Then, Equation (2.5) is written as: cib pi= Q ∑ j=1"(ZΣj fb pi∗φ φ φjdΣ)b ui−(ZΣj fb ui∗φ φ φjdΣ)b pi#= Q ∑ j=1hb Hi j b ui−b Gi j b pii(2.6) where Q is the number of boundary nodes at the boundary Σf and Σj f are the elements which contains the node j . After interpolating the boundary variables, the integral representation defined by Equation (2.6) yields a system of equations that is solved for each frequency. The spatial integrations in Equation (2.6) are numerically evaluated by a Lobatto-Gauss-Legendre (LGL) quadrature when the collocation point does not belong to the integration element. In a first approach, a LGL quadrature is also used to evaluate the singular integrals when the collocation point belongs to the integration element using an appropiate number of integration points. References [ 21 , 22 ] state that the use of a sufficient number of quadrature points leads to a good approximation of the singular integrals. The system of equations for all the boundary elements becomes: b H(ω,kz)b u(bx,ω,kz) = b G(ω,kz)b p(bx,ω,kz)(2.7) The proposed spectral boundary element method for the 2.5D fluid element uses Legendre polynomials of order pas interpolation shape functions. The shape interpolation functions φ φ φare given by: φi=∏ j6=i ξ−ξi ξj−ξi (2.8) where the local nodal coordinates ξare found at the LGL integration points: (1−ξ2)∂φ φ φ(ξ) ∂ξ =0(2.9) The element integration is carried out using the LGL points. Then, the computation of matrices b H y b G is simplified as follows: bhi j =b pi∗(bx,ω,kz;b xi)wjJ(ξj) b gi j =b ui∗(bx,ω,kz;b xi)wjJ(ξj)(2.10) where wj is the LGL integration weight at the natural nodal coordinate ξj , and J(ξj) is the Jacobian evaluated at ξj , b ui∗(bx,ω,kz;b xi) and b pi∗(bx,ω,kz;b xi) are respectively normal displacement to boundary Γf and nodal pressure, and both are derived from a potential: b ui∗(bx,ω,kz;b xi) = −i 4H(2) 1(kαr)∂r ∂n b pi∗(bx,ω,kz;b xi) = i 4H(2) 0(kαr)(2.11) 2.2 The 2.5D spectral finite element formulation 5 (a) (b) (c) x y x y ξ 1 2 345 6 7 Figure 2.2 (a) Element mesh seed representation, (b) element definition and (c) physical representation of spectral fluid boundary element of order p=6. ξ -1 -0.5 0 0.5 1 -0.2 0 0.2 0.4 0.6 0.8 1φ1 φ2 φ3 φ4 φ5 φ6 φ7 Figure 2.3 Element shape functions for an element with order p=6. where H(2) n are second Hankel funtions of order n , kα=qω2 α2−k2 z , with Im(kα)<0 , α is the sound propagation velocity and r=p(x−x0)2+(y−y0)2. The boundary Γf is discretized into elements that are generated from a base mesh. The element definition is done by using a polynomial interpolation from the element base mesh and its representation in the natural coordinates system. As an example, Figure 2.2.(a) shows an element base mesh defined by four nodes used to generate an element with order p=6 (Figure 2.2.(c)). The spectral element nodes are defined in natural coordinates in ascending order (Figure 2.2.(b)) with ξ∈[−1,1]. The shape function φ φ φ and its derivatives are symbolically computed in straightforward procedure for an arbitrary element order. Figure 2.3 shows the interpolation shape functions for the p=6 spectral element presented in Figure 2.2. 2.2 The 2.5D spectral finite element formulation The spectral formulation for solid subdomains is based on the virtual work principle. For any virtual displacement field δu imposed on the solid subdomain Ωs , the virtual work due to the internal and the inertial forces must be equal to the virtual work produced by the external loads [25]: −ω2ZΩs δu(x,ω)ρsu(x,ω)dΩ+ZΩs δε ε ε(x,ω)σ σ σ(x,ω)dΩ =ZΩs δu(x,ω)ρsb(x,ω)dΩ+ZΓs δu(x,ω)f(x,ω)dΓ (2.12) where u(x,ω) is the displacement vector, ε ε ε and σ σ σ are respectively the strain and stress tensors, ρsb(x,ω) is the body force in the domain Ωs , ρs is the solid density, and f(x,ω) is the nodal force vector. A variable preceded by δdenotes a virtual change of its magnitude. The Voigt notation is used to write both the symmetrical stress tensor σ σ σ and the strain tensor ε ε ε , σ σ σ= {σxx,σyy,σzz,σxy,σyz,σzx}T and ε ε ε={εxx,εyy,εzz,γxy,γyz,γzx}T , respectively. The use of engineering shear strains allows to write the internal work as an inner product εT bσb. 6Chapter 2. Numerical model The displacement vector u(x,ω)is discretized as: u(x,ω) = Nu(x,ω)(2.13) where Nare the globally defined shape functions and u(x,ω)is the nodal displacement vector. The strain tensor ε ε εis derived from the nodal displacement vector u(x,ω)as: ε ε ε=L1Nu +L2N∂u ∂z(2.14) where L1=          ∂ ∂x0 0 0∂ ∂y0 000 ∂ ∂y ∂ ∂x0 0 0 ∂ ∂y 0 0 ∂ ∂x          (2.15) and L2=         000 000 001 000 010 100         (2.16) The strain vector ε ε ε in Equation (2.14) is related to both the displacement vector u and its derivative ∂u/∂z with respect to the longitudinal coordinate z . This differs from the relationship ε ε ε=LNu for the 3D formulation, where the strain vector is written as a linear combination of the elements of the nodal displacement vector u . Equation (2.14) is alternatively written as: ε ε ε=B1u+B2 ∂u ∂z(2.17) where B1=L1Nand B2=L2N. The stress tensor is related to the strain vector through the constitutive relation: σ σ σ=Dε ε ε(2.18) where, in the case of a linear isotropic material, the constitutive matrix D depends on the Young’s modulus E and the Poisson’s ratio ν: D=E (1+ν)(1−2ν)         1−ν ν ν 0 0 0 ν1−ν ν 0 0 0 ν ν 1−ν0 0 0 0001−2ν 20 0 0 0 0 0 1−2ν 20 0 0 0 0 0 1−2ν 2         (2.19) A standard Galerkin procedure is followed, using the same approximation for the virtual displacement vector. Substituting the strain-displacement relation defined by Equation (2.17), and the constitutive Equation (2.18) into the virtual work Equation (2.12) yields: −ω2ZΩs δuTNTρsNudΩ+ZΩs δuTBT 1+δ∂u ∂zT BT 2!DB1u+B2 ∂u ∂zdΩ =ZΩb δuTNTρsbdΩ+ZΓf δuNTfdΓ (2.20) 2.2 The 2.5D spectral finite element formulation 7 Equation (2.20) is further elaborated, rewriting the volume integrals over the longitudinal coordinate z and the cross-section As: −ω2Z∞ −∞ δuTZAs NTρsNdAudz +Z∞ −∞ δuTZAs BT 1DB1dAudz +Z∞ −∞ δuTZAs BT 1DB2dA∂u ∂zdz +Z∞ −∞ δ∂u ∂zTZAs BT 2DB1dAudz +Z∞ −∞ δ∂u ∂zTZAs BT 2DB2dA∂u ∂zdz =Z∞ −∞ δuTZAs NTρsbdAdz +Z∞ −∞ δuTZΣf NTfdΣdz (2.21) where Σf is the intersection of the surface Γf with the plane z=0 (Figure 2.1). The discretized equation is obtained through elimination of the virtual displacement vector δu . This requires integration by parts on the terms containing derivatives δ(∂u/∂z)in Equation (2.21): −ω2Z∞ −∞ δuTZAs NTρsNdAudz +Z∞ −∞ δuTZAs BT 1DB1dAudz +Z∞ −∞ δuTZAs BT 1DB2dA −ZAs BT 2DB1dA∂u ∂zdz −Z∞ −∞ δuTZAs BT 2DB2dA∂2u ∂z2dz =Z∞ −∞ δuTZAs NTρsbdAdz +Z∞ −∞ δuTZΣf NTfdΣdz (2.22) Since Equation (2.22) holds for any virtual displacement δu , all integrals over the longitudinal coordinate zvanish and Equation (2.22) is written as follows: −ω2Mu +K0u+K1∂u ∂z+K2∂2u ∂z2=f(2.23) where the mass matrix Mis defined as: M=ZAs NTρsNdA (2.24) and stiffness matrices K0,K1and K2are defined as: K0=ZAs BT 1DB1dA K1=ZAsBT 1DB2−BT 2DB1dA K2=ZAs BT 2DB2dA (2.25) The external load vector: f=ZAs NTρsb(x,y,z,ω)dA +ZΣf NTf(x,y,z,ω)dΣ(2.26) contains contributions for both body forces and surface loads and it is evaluated for every point z along the longitudinal axis. The differential Equation (2.23) is solved by a Fourier transform of the longitudinal coordinate z to the horizontal wavenumber kz , where the Fourier transform is defined as Ff(z),kz=R∞ −∞exp(+ikzz)f(z)dz . Equation (2.23) is transformed to the wavenumber domain as: −ω2M+K0−ikzK1−k2 zK2b u(kz,ω) =bf(kz,ω)(2.27) 8Chapter 2. Numerical model (a) (b) (c) x y x y ξ η 1 8 15 22 29 36 43 2 9 16 23 30 37 44 3 10 17 24 31 38 45 4 11 18 25 32 39 46 5 12 19 26 33 40 47 6 13 20 27 34 41 48 7 14 21 28 35 42 49 Figure 2.4 (a) Element mesh seed representation, (b) element definition and (c) physical representation of a solid spectral element of order p=6. 1 0 ξ [-] -1 -1 η [-] 0 1 Figure 2.5 Element shape function of p=6N25 for solid finite element. where a hat above a variable denotes its representation in the frequency-wavenumber domain. Equation (2.27) is rewritten if an equivalent dynamic stiffness matrix b Kis considered: b K(kz,ω)b u(bx,kz,ω) =bf(bx,kz,ω)(2.28) where the dynamic stiffness matrix b Kis easily computed for each wavenumber and frequency step. The proposed spectral finite element method uses the 2D Legendre polynomials of order p as specific shape functions N , and the local nodal element coordinates ξ ξ ξ(ξ,η)∈[−1,1]x[−1,1] are defined as the Lobatto-Gauss-Legendre (LGL) points: ((1−ξ2)∂N(ξ,η) ∂ξ =0 (1−η2)∂N(ξ,η) ∂η =0(2.29) Then, the local nodal element coordinate matches the LGL integration points, and defines an orthogonal element shape function basis. Therefore, the solid element matrices defined by Equations (2.24) and (2.25) are easily computed as follows: mi j =δi jρwi j J(ξi,ηj) k2 i j =δi jLT 2DL2wi j J(ξi,ηj)(2.30) where δi j is the Kronecker delta, wi j is the LGL weight at the nodal coordinate system (ξi,ηj) , and J(ξi,ηj) is the Jacobian evaluated at (ξi,ηj). The solid subdomain Ωs is discretized into elements that are generated from a base element mesh. The element definition is done by using a polynomial interpolation from the starting element geometry and its representation in the natural coordinates system. As an example, Figure 2.4. shows an element with order p=6 . The element nodes are defined in natural coordinates in ascending order (Figure 2.4.(b)) with ξ ξ ξ(ξ,η)∈[−1,1]×[−1,1] . The spectral solid element has three degrees of freedom at each node, b u(bx,kz,ω) = {b ux(bx,kz,ω),b uy(bx,kz,ω),b uz(bx,kz,ω)}. The element shape functions are computed as Nk=φ φ φi(ξ)φ φ φj(η) , where k is the element node in which the shape function takes unity value, k= (i−1)(p+1) + j , and φ φ φi is the one-dimensional LGL shape function defined by Equation (2.8). The shape function Nk and its derivatives are symbolically computed in straightforward procedure for an arbitrary element order. Figure 2.5 shows the interpolation shape function N25 for an spectral element with order p=6(Figure 2.4.(b)). 3 Numerical verification The proposed spectral boundary and finite element methods were verified with two benchmark problems regarding to a fixed cylindrical inclusion in an unbounded fluid medium and the propagation of elastic waves in free rail. 3.1 Fixed cylindrical inclusion The BEM model was verified with a benchmark problem. The model was implemented and validated by applying it to a fixed cylindrical circular cavity, submerged in a homogeneous unboundend fluid medium. The cavity is subjected to a harmonic point pressure load. The analytical solution to this problem can be found in Reference [20]. The cavity has a radius r=5 m , and the unbounded fluid medium properties are pressure wave velocity α=1500 m/sand density ρ=1000 kg/m3. The problem was solved for a frequency range varying from 2.5 Hz to 320 Hz, with a frequency step ∆f=2.5 Hz . A constant value of kz=10−5rad/m was assumed. The problem solution was computed directly from Equation (2.7), and then the radiated wave field at receivers was assessed using the Somigliana identity [9]. The problem solution was computed for a dilatational point source placed at the fluid medium bx0= (x0,y0) = (0,15)15 m away from the cavity centre. This loads emits a harmonic incident field b pinc at a point bxdescribed by: b pinc = (bx,ω,kz) = −iA 2H(2) 0(kαq(x−x0)2+(y−y0)2)(3.1) where A is the source amplitude, H(2) n are the second Hankel functions of order n and kα=qω2 α2−k2 z is the effective wavenumber (Im{kp}=0). The element size and the approximation order was set to keep the nodal density per wavelength at least dλ=12 [ 4 ]. The characteristic problem wavelength was λmin =α/fmax =4.69 m , so to accomplish the above condition the minimum number of nodes must be 80. Three different boundary discretization were investigated for characteristic element distribution 1/h=2/πm−1 , 1/h=3/πm−1 and 1/h=4/πm−1 . In this study, the element order is p=5 . As an example, Figure 3.1 shows a boundary discretization using 1/h=2/πm−1and p=5. The scattered wave field was computed over a grid of receivers as shown in Figure 3.2. Figure 3.3 shows the amplitude of the analytical solution of the scattered pressure field when acts a harmonic pressure load of 200 Hz and the total incident fields. The maximum amplitude was found near to the source point, and a shadowed region appeared after the cavity. Figure 3.4 represents the difference between the analytical solution and computed results when the inclusion is modelled with a different number of boundary elements. Maximum errors occur in the shadowed zone and around the surface inclusion. The improvement of the computed solution as the distance to the surface of the cavity increases is noticeable. As expected, the boundary element method accuracy improves as element size decreases. The results showed that higher order elements or small element size must be used to compute the solution at points close to the cavity boundary, where the scattered field changes rapidly. 9 10 Chapter 3. Numerical verification Figure 3.1 Cylindrical cavity embedded in an unbounded fluid medium: problem discretization with element distribution 1/h=2/πand order p=6. x [m] -10 -5 0 5 10 y [m] -10 -5 0 5 10 Figure 3.2 Receivers grid. (a) (b) Figure 3.3 (a) Scattered pressure field and (b) total pressure field due to a harmonic pressure load acting at 200 Hz. 3.2 Propagative waves in a free rail The solid spectral element presented in Section 2.1 was verified through computing of propagative waves in a free rail UIC861-3. This problem was previously studied by L. Gavric [ 10 ]. The dispersion curves and the cross-section modes were computed for the wavenumber range 0−35 rad/m , with a wavenumber step ∆kz=0.05 rad/m. First, the problem was solved using linear finite elements with two different discretizations. The first mesh consists of 227 nodes and 176 finite elements, which results in a mesh of 681 degrees of freedom (Figure 3.5.(a)). The element size varies from h=5.7 mm to h=8.6 mm . This mesh is similar to the mesh used in 3.2 Propagative waves in a free rail 11 (a) (b) (c) (d) (e) (f) Figure 3.4 Absolute error obtained from discretizations (a-b) 1/h=2/πm−1, (c-d) 1/h=3/πm−1and (f-g) 1/h=4/πm−1. Reference [ 10 ]. The element size in the second case varies from h=2.9 mm to h=4.3 mm , which means 805 nodes and 704 elements with 2415 degrees of freedom (Figure 3.5.(b)). The computed results are compared with the reference solution [ 10 ]. Eight characteristic cross-sectional modes shapes were computed in the frequency range 0−6 kHz . Four modes had axial/vertical displacement componets (Figure 3.6), while the rest were characterized by dominant lateral in plane vibrations of the cross-section (Figure 3.7). For the low frequency range, the computed results show a good agreement with reference solution, regadless of elements size. Some differences appear between numerical results and reference results for higher frequencies. Furthermore, an analysis using spectral finite elements has been made. The cross-section has been represented by seven elements, three elements to model the railhead, three elements to model the foot and one element to model the neck. This is an advantage of the SEM, which allows to discretize the domain with fewer elements than the FEM. The problem was solved using different element order to assess the accuracy. Figure 3.8 shows the mesh with element order p=4and p=8. Computed results show a good agreement for most of the previously described waves. The lateral waves 4 Conclusions and further developement This work has proposed two spectral based formulations based on the FEM and the BEM to study fluid and solid wave propagation. Both models look at 3D problems whose materials and geometric properties remain homogeneous in one direction. Solid subdomain was modelled with the FEM, whereas the boundary of unbounded fluid media was modelled with the BEM. A spectral 2.5D FEM element for solid media and a spectral 2.5D BEM element for fluid-acoustics media were developed. The models were verified with two benchmark problems with known solutions. The BEM was verified with a problem consisting of a cavity embedded in an unbounded fluid medium subjected to an incident wave field. The FEM was verified by computation the dispersion curves of propagative waves in free rail. Numerical results were in good agreement with the reference results. Future researches should complete the work presented herein. The first step is improvement of the BEM model. This can be achieved by solving the singular integrals using more accurate methods, such as closed forms or a regularization procedure. The next step is to complete the coupling procedure in order to study fluid-solid interaction problems. Then, the radiated wavefield in an unbounded fluid medium caused by the deformation of a solid structure due to a point load can be studied. 19 List of Figures 2.1 Domain decomposition 3 2.2 (a) Element mesh seed representation, (b) element definition and (c) physical representation of spectral fluid boundary element of order p=65 2.3 Element shape functions for an element with order p=65 2.4 (a) Element mesh seed representation, (b) element definition and (c) physical representation of a solid spectral element of order p=68 2.5 Element shape function of p=6N25 for solid finite element 8 3.1 Cylindrical cavity embedded in an unbounded fluid medium: problem discretization with element distribution 1/h=2/πand order p=610 3.2 Receivers grid 10 3.3 (a) Scattered pressure field and (b) total pressure field due to a harmonic pressure load acting at 200 Hz 10 3.4 Absolute error obtained from discretizations (a-b) 1/h=2/πm−1, (c-d) 1/h=3/πm−1and (f-g) 1/h=4/πm−111 3.5 (a) Rail discretization for linear element size h=5.7−8.6 mm elements and (b) h=2.9−4.3 mm 12 3.6 Dispesion curves of propagative vertical/longitudinal modes for the frequency range 0-6000 Hz. , computed results for h=2.9−4.3 mm; , computed results for h=5.7−8.6 mm;4,,◦,+, reference results [10] 12 3.7 Dispesion curves of propagative lateral modes for the frequency range 0-6000 Hz. , computed results for h=2.9−4.3 mm; , computed results for h=5.7−8.6 mm;4,,◦,+, reference results [10] 12 3.8 Rail discretization with spectral elements. Element order (a) p=4and (b) p=813 3.9 Dispesion curves of propagative vertical/longitudinal modes for the frequency range 0-6000 Hz. , computed results for p=8; , computed results for p=4;4,,◦,+, reference results [10] 13 3.10 Dispesion curves of propagative lateral modes for the frequency range 0-6000 Hz. , computed results for p=8; , computed results for p=4;4,,◦,+, reference results [10] 14 3.11 Axionometric representation of the cross-section modes of three axial waves at excitation frequency 5600 Hz; (a) ax1wave, (b) ax2wave, (c) ax3wave 14 3.12 Front view of the cross-section mode of the vertical flexural wave f v at excitation frequencies (a) 450 Hz, (b) 1500 Hz and (c) 4200 Hz 15 3.13 Front view of the cross-section mode of the horizontal flexural wave f h at excitation frequencies (a) 450 Hz, (b) 1500 Hz and (c) 4200 Hz 15 3.14 Front view of the cross-section mode of the torsional wave to at excitation frequencies (a) 450 Hz, (b) 1500 Hz and (c) 4200 Hz 16 3.15 Front view of the cross-section mode of the propagative vawe c1at excitation frequencies (a) 1500 Hz, (b) 2800 Hz and (c) 4200 Hz 16 3.16 Front view of the cross-section mode of the propagative wave c2at excitation frequencies (a) 4200 Hz and (b) 5200 Hz 17 21 Bibliography [1] I. Babuška, F. Ihlenburg, E.T. Paik, and S.A. Sauter, A generalized finite element method for solving the helmholtz equation in two dimensions with minimal pollution, Computer Methods in Applied Mechanics and Engineering 128 (1995), no. 3-4, 325–359. [2] I. Babuška and M. Suri, P and H-P versions of the finite element method, basic principles and properties, SIAM Review 36 (1994), no. 4, 578–632. [3] P.Z. Bar-Yoseph, D. Fisher, and O. Gottlieb, Spectral element methods for nonlinear spatio-temporal dynamics of an euler-bernoulli beam, Computational Mechanics 19 (1996), no. 1, 136–151. [4] H. Bériot, G. Gabard, and E. Perrey-Debain, Analysis of high-order finite elements for convected wave propagation, International Journal for Numerical Methods in Engineering 96 (2013), no. 11, 665–688. [5] J. Biermann, O. von Estorff, S. Petersen, and C. Wenterodt, Higher order finite and infinite elements for the solution of helmholtz problems, Computer Methods in Applied Mechanics and Engineering 198 (2009), no. 13–14, 1171 – 1188. [6] M. Bouchon and K. Aki, Discrete wave-number representation of seismic source wavefields, Bulletin of the Seismological Society of America 67 (1977), no. 2, 259–277. [7] G. Degrande and G. De Roeck, A spectral element method for two-dimensional wave propagation in horizontally layered saturated porous media, Computers & Structures 44 (1992), no. 4, 717 – 728. [8] A. Deraemaeker, I. Babuška, and P. Bouillard, Dispersion and pollution of the fem solution for the helmholtz equation in one, two and three dimensions, International Journal for Numerical Methods in Engineering 46 (1999), no. 4, 471–499. [9] J. Domínguez, Boundary elements in dynamics, Computational Mechanics Publications and Elsevier Aplied Science, Southampton, 1993. [10] L. Gavrić, Computation of propagative waves in free rail using a finite element technique, Journal of Sound and Vibration 185 (1995), no. 3, 531–543. [11] F. Ihlenburg, The medium-frequency range in computational acoustics: Practical and numerical aspects, Journal of Computational Acoustics 11 (2003), no. 2, 175–193. [12] F. Ihlenburg and I. Babuška, Finite element solution of the helmholtz equation with high wave number Part I: The h-version of the FEM, Computers and Mathematics with Applications 30 (1995), no. 9, 9–37. [13] F. Ihlenburg and I. Babuška, Finite element solution of the helmholtz equation with high wave number part II: the h-p version of the FEM, SIAM Journal on Numerical Analysis 34 (1997), no. 1, 315–358. [14] F. Ihlenburg, I. Babuška, and S. Sauter, Reliability of finite element methods for the numerical computation of waves, Advances in Engineering Software 28 (1997), no. 7, 417–424. 23 [15] P. Kudela, M. Krawczuk, and W. Ostachowicz, Wave propagation modelling in 1d structures using spectral finite elements, Journal of Sound and Vibration 300 (2007), no. 1–2, 88 – 100. [16] O.Z. Mehdizadeh and M. Paraschivoiu, Investigation of a two-dimensional spectral element method for helmholtz’s equation, Journal of Computational Physics 189 (2003), no. 1, 111 – 129. [17] A.A. Oberai and P.M. Pinsky, A numerical comparison of finite element methods for the helmholtz equation, Journal of Computational Acoustics 08 (2000), no. 01, 211–221. [18] A.T. Patera, A spectral element method for fluid dynamics: Laminar flow in a channel expansion, Journal of Computational Physics 54 (1984), no. 3, 468–488. [19] S. Petersen, D. Dreyer, and O. von Estorff, Assessment of finite and spectral element shape functions for efficient iterative simulations of interior acoustics, Computer Methods in Applied Mechanics and Engineering 195 (2006), no. 44–47, 6463 – 6478. [20] A.J.B. Tadeu and L.M.C. Godinho, Three-dimensional wave scattering by a fixed cylindrical inclusion submerged in a fluid medium, Engineering Analysis with Boundary Elements 23 (1999), no. 9, 745–755. [21] A.J.B. Tadeu, P.F.A. Santos, and E. Kausel, Closed-form integration of singular terms for constant, linear and quadratic boundary elements. Part 1. SH wave propagation, Engineering Analysis with Boundary Elements 23 (1999), no. 8, 671–681. [22] A.J.B. Tadeu, P.F.A. Santos, and E. Kausel, Closed-form integration of singular terms for constant, linear and quadratic boundary elements. Part 2. SV-P wave propagation, Engineering Analysis with Boundary Elements 23 (1999), no. 9, 757–768. [23] P.E.J. Vos, S.J. Sherwin, and R.M. Kirby, From h to p efficiently: Implementing finite and spectral/hp element methods to achieve optimal performance for low- and high-order discretisations, Journal of Computational Physics 229 (2010), no. 13, 5161 – 5181. [24] C. Zhu, G. Qin, and J. Zhang, Implicit chebyshev spectral element method for acoustics wave equations, Finite Elements in Analysis and Design 47 (2011), no. 2, 184 – 194. [25] O.C. Zienkiewicz, The finite element method, third ed., McGraw-Hill, 1986. 24 Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Máster Máster Universitario en Ingeniería de Caminos, Canales y Puertos Un método espectral formulado en 2.5D basado en el MEC y el MEF para estudiar propagación de ondas en medios acústicos y elásticos. Resumen Autor: Francisco Javier Cruz Muñoz Tutor: Antonio Romero Ordóñez Tutor: Pedro Galvín Barrera Dep. Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016 Resumen En este trabajo se presenta una formulación espectral en dos dimensiones y media (2.5D) basada en el método de los elementos finitos (MEF) y el método de los elementos de contorno (MEC) para estudiar la propagación de ondas en los ámbitos de la acústica de fluidos y la elastodinámica tridimensionales (3D). El análisis se lleva a cabo mediante la superposición de problemas bidimensionales (2D) para diferentes números de onda a lo largo de la coordenada longitudinal. El método numérico está basado en la descomposición del dominio para representar el sistema fluido-estructura. Se presenta en este trabajo un elemento finito espectral para representar ondas guiadas en sólidos con sección transversal arbitraria. Además, el MEC se extiende a la formulación espectral para estudiar medios fluidos infinitos. Ambas aproximaciones utilizan polinomios de interpolación de Lagrage como funciones de forma en los puntos de Legendre-Gauss-Lobatto (LGL). Las técnicas propuestas se han verificado mediante dos problemas de referencia: la propagación de ondas en un medio fluido infinito en presencia de una cavidad rígida fija y la propagación de ondas en un carril libre, y los resultados obtenidos presentan un buen grado de acuerdo con las soluciones de referencia. Se realiza un breve análisis h−ppara evaluar la precisión de los métodos. I Índice Resumen I 1 Introducción 1 1.1 Esquema del problema y motivación 1 1.2 Objetivos 1 1.3 Contribuciones originales 2 1.4 Organización del texto 2 2 Modelo numérico 3 2.0.1 La formulación espectral 2.5D del método de los elementos de contorno 3 2.0.2 La formulación espectral 2.5D del método de los elementos finitos 4 3 Verificación del modelo numérico 9 3.0.1 Cavidad cilíndrica en un medio fluido 9 3.0.2 Propagación de ondas en un carril libre 10 4 Conclusiones y desarrollos futuros 15 Índice de Figuras 17 Bibliografía 19 III 7 (a) (b) (c) x y x y ξ η 1 8 15 22 29 36 43 2 9 16 23 30 37 44 3 10 17 24 31 38 45 4 11 18 25 32 39 46 5 12 19 26 33 40 47 6 13 20 27 34 41 48 7 14 21 28 35 42 49 Figura 2.3 (a) Representación de la malla base, (b) definición del elemento en coordenadas naturales y (c) representación en coordenadas físicas de un elemento espectral de orden p=6. 3 Verificación del modelo numérico Los métodos espectrales de elementos finitos y de contorno han sido verificados mediante la resolución de dos problemas de referencia, una cavidad cilíndria sumergida en un medio fluido infinito, y la propagación de ondas en un carril libre. 3.0.1 Cavidad cilíndrica en un medio fluido El modelo MEC se ha verificado con un problema de referencia. El modelo se ha validado mediante la solución de una cavidad cilíndrica, sumergida en un medio fluido homogéneo. La cavidad está sometida a una carga de presión puntual armónica. La solución analítica se puede encontrar en la Referencia [20]. La cavidad tiene un radio r=5 m , y las propiedades del fluido son la velocidad de propagación de las ondas α=1500 m/sy la densidad ρ=1000 kg/m3. El problema se ha resuelto para un rango de frecuencias de entre 2.5 Hz hasta 320 Hz, con un paso de ∆f=2.5 Hz . Se ha tomado un valor constante kz=10−5rad/m . El campo de presiones radiado se ha calculado mediante la identidad de Somigliana.[9]. La solución del problema se ha calculado para una fuente sutiada en el fluido en el punto bx0= (x0,y0) = (0,15)a15 m del centro de la cavidad. El tamaño de los elementos y el orden de aproximación se ha establecido para mantener una densidad de nodos por longitude de onda de al menos dλ=12 [4]. Se han considerado tres discretizaciones diferentes, con una longitud de elementos de 1/h=2/πm−1 , 1/h=3/πm−1y1/h=4/πm−1. En este trabajo se ha tomado un orden de elemento de p=5. Figura 3.1 Cavidad cilíndrica embebida en un medio fluido infinito: discretización con una distribución de elementos 1/h=2/πy orden p=6. El campo de ondas radiado se ha calculado en una malla de receptores como se muestra en la Figura ?? . La Figura 4.3 muestra la amplitud de la solución analítica del campo de presiones radiado y el campo total de presiones cuando actua una presión puntual de 200 Hz. La Figure 4.4 representa el error absoluto entre los resultados numéricos y la solución analítica para las diferentes discretizaciones estudiadas. Los máximos errores aparecen en la zona de sombra y alrededor del contorno de la inclusión. La solución numérica mejora cuanto mayor sea la distancia al contorno. Como era de esperar, la precisión del modelo 9 10 Capítulo 3. Verificación del modelo numérico x [m] -10 -5 0 5 10 y [m] -10 -5 0 5 10 Figura 3.2 Malla de receptores. (a) (b) Figura 3.3 (a) Campo de presiones radiado y (b) campo total de presiones debido a una carga armónica actuando a 200 Hz. mejora al disminuir el tamaño de los elementos. Los resultados muestran que un orden de aproximación mayor o un menor tamaño de elementos debe usarse para calcular la solución en puntos cercanos al contorno. 3.0.2 Propagación de ondas en un carril libre El elemento sólido espectral presentado en la Sección 2.0.1 se ha verificado mediante la el cálculo de las ondas propagativas en un carril libre UIC861-3. Este problema ya ha sido estudiado por L. Gavric [ 10 ]. Las curvas de dispersión y los modes de vibración se han calculado para un rango de número de onda 0−35 rad/m , con un paso de ∆kz=0.05 rad/m. La sección transversal se ha representado mediante siete elementos, tres para la cabeza del carril, tres para el patín y uno para el alma. Esta es una ventaja del metodo de los elementos finitos espectrales, que permite discretizar el dominio con un menor número de elementos. El problema ha sido resuelto con diferentes órdenes de aproximación para evaluar la precisión del método. Los resultados muestran que aparecen ocho modos de vibración característicos en el rango de frecuencias 0−6 kHz . Cuatro modos tienen componentes de desplazamiento axiales/verticales (Figure ?? ), mientras que el resto se caracterizan por un movimiento lateral de la sección transversal (Figure ??). Los resultados numéricos muestran un buen grado de acuerdo con los resultados de referencia para casi todas las ondas descritas. Las ondas de movimiento lateral son muy parecidas a la solución de referencia. Por otra parte, las ondas ax1 , ax2 y ax3 muestran diferencias apreciables con los resultados presentados en la Referencia [ 10 ] para el rango de frecuencias más elevado. Estas diferencias podrían ser explicadas por la malla gruesa utilizada en la Referencia [ 10 ] para frecuencias tan altas. Los resultados propuestos sugieren que la solución para mallas finas convergen a la solución del SEM que aparece en las Figuras 4.6 y 4.7. Estos resultados están en concordancia con los descritos en la Referencia [ 10 ]. De modo que la precisión del método propuesto queda verificadaa. 11 (a) (b) (c) (d) (e) (f) Figura 3.4 Error absoluto obtenido para las discretizaciones (a-b) 1/h=2/πm−1, (c-d) 1/h=3/πm−1y (f-g) 1/h=4/πm−1. x y z (a) (b) Figura 3.5 Discretización del carril con elementos espectrales. Orden de elementos (a) p=4y (b) p=8. 12 Capítulo 3. Verificación del modelo numérico Frequency [kHz] 0123456 Wavenumber [rad/m] 0 5 10 15 20 25 30 35 fv ax1 ax2 ax3 Figura 3.6 Curvas de dispersión de los modos axiales/verticales para el rango de frecuencias 0-6000 Hz. , resultados numéricos para p=8 ; , resultados numéricos para p=4 ; 4 ,  , ◦ , + , resultados de referencia [10]. Frequency [kHz] 0123456 Wavenumber [rad/m] 0 5 10 15 20 25 30 35 fh to c1 c2 Figura 3.7 Curvas de dispersión de los modos laterales para el rango de frecuencias 0-6000 Hz. , resultados numéricos para p=8 ; , resultados numéricos para p=4 ; 4 ,  , ◦ , + , resultados de referencia [10]. 13 (a) (b) (c) Figura 3.8 Vista frontal del modo de la onda axiales a la frecuencia de 5600 Hz; (a) ax1, (b) ax2, (c) ax3. (a) (b) (c) Figura 3.9 Vista frontal del modo de la onda f v a las frecuencias (a) 450 Hz, (b) 1500 Hz y (c) 4200 Hz. 14 Capítulo 3. Verificación del modelo numérico (a) (b) (c) Figura 3.10 Vista frontal del modo de la onda f h a las frecuencias (a) 450 Hz, (b) 1500 Hz y (c) 4200 Hz. (a) (b) (c) Figura 3.11 FVista frontal del modo de la onda to a las frecuencias (a) 450 Hz, (b) 1500 Hz y (c) 4200 Hz. (a) (b) (c) Figura 3.12 Vista frontal del modo de la onda c1a las frecuencias (a) 1500 Hz, (b) 2800 Hz y (c) 4200 Hz. 15 (a) (b) Figura 3.13 Vista frontal del modo de la onda c2a las frecuencias (a) 4200 Hz y (b) 5200 Hz. 4 Conclusiones y desarrollos futuros En este trabajo se han propuesto dos fomulaciones espectrales basadas en el MEC y en el MEF para estudiar propagación de ondas. Ambos modelos resuelven problemas tridimensionales con propiedades del material y geometría constante a lo largo de la coordenada longitudinal. Los subdominios sólidos se representan mediante el MEF, mientras que el contorno del medio fluido se modeliza con el MEC. Los elementos finitos y de contorno espectrales en 2.5D han sido desarrollados. Los modelos han sido verificados con dos problemas de referencia, de soluciones conocidas. El MEC ha sido validado mediante un problema consistente en una cavidad embebida en un medio fluido infinito sometido a un campo de presiones incidente. El MEF se ha verificado mediante el cálculo de las curvas de dispersión de ondas propagativas en un carril libre. Los resultados numéricos presentaron un buen grado de acuerdo con los resultados de referencia. Los desarrollos futuros debería completar el trabajo aquí presentado. El primer paso es la mejora del modelo MEC. Esto puede alcanzarse resolviendo las integrales singulares con métodos más precisos, como resolviéndolas analíticamente o mediante un proceso de regularización. El siguiente paso consiste en el proceso de acoplamiento de ambos métodos para estudiar problemas de interacción fluido-estructura. Se podrá estudiar el problema del campo radiado en un fluido infinito causado por la deformación de una estructura debido a una carga puntual. 17 Índice de Figuras 2.1 Descomposición del dominio 3 2.2 (a) Representación de la malla base, (b) definición del elemento en coordenadas naturales y (c) representación en coordenadas físicas del elemento de contorno espectral fluido de orden p=64 2.3 (a) Representación de la malla base, (b) definición del elemento en coordenadas naturales y (c) representación en coordenadas físicas de un elemento espectral de orden p=67 3.1 Cavidad cilíndrica embebida en un medio fluido infinito: discretización con una distribución de elementos 1/h=2/πy orden p=69 3.2 Malla de receptores 10 3.3 (a) Campo de presiones radiado y (b) campo total de presiones debido a una carga armónica actuando a 200 Hz 10 3.4 Error absoluto obtenido para las discretizaciones (a-b) 1/h=2/πm−1, (c-d) 1/h=3/πm−1y (f-g) 1/h=4/πm−111 3.5 Discretización del carril con elementos espectrales. Orden de elementos (a) p=4y (b) p=811 3.6 Curvas de dispersión de los modos axiales/verticales para el rango de frecuencias 0-6000 Hz. , resultados numéricos para p=8; , resultados numéricos para p=4;4,,◦,+, resultados de referencia [10] 12 3.7 Curvas de dispersión de los modos laterales para el rango de frecuencias 0-6000 Hz. , resultados numéricos para p=8; , resultados numéricos para p=4;4,,◦,+, resultados de referencia [10] 12 3.8 Vista frontal del modo de la onda axiales a la frecuencia de 5600 Hz; (a) ax1, (b) ax2, (c) ax312 3.9 Vista frontal del modo de la onda f v a las frecuencias (a) 450 Hz, (b) 1500 Hz y (c) 4200 Hz 13 3.10 Vista frontal del modo de la onda f h a las frecuencias (a) 450 Hz, (b) 1500 Hz y (c) 4200 Hz 13 3.11 FVista frontal del modo de la onda to a las frecuencias (a) 450 Hz, (b) 1500 Hz y (c) 4200 Hz 13 3.12 Vista frontal del modo de la onda c1a las frecuencias (a) 1500 Hz, (b) 2800 Hz y (c) 4200 Hz 14 3.13 Vista frontal del modo de la onda c2a las frecuencias (a) 4200 Hz y (b) 5200 Hz 14 19