Error estimation in vibroacoustic problems solved by means of finite elements

1 Error estimation in vibroacoustic problems solved by means of finite elements Poblet-Puig, J.; Rodríguez-Ferran, A. Affiliation: {Laboratori de Càlcul Numèric, E.T.S. d'Enginyers de Camins, Canals i Ports de Barcelona, Universitat Politècnica de Catalunya} e-mail: {[email protected]; antonio.rodriguez-ferr[email protected]u} Abstract The vibroacoustic equations can be solved by means of the finite element method. A discretisation of the structure and the acoustic domains is required and highly influences the quality of the numerical solution. There exist meshing criteria (a priori error estimators) for the case of the Helmholtz equation but these studies have not focused their attention in the case of the vibroacoustic problem. The fluid structure interaction represents a new source of numerical errors and meshes in the interaction zone should be designed by not only taking into account the physical properties of the acoustic medium but also the mechanical properties of the structure. The goal of the work is to obtain an a priori error estimation criterion for the vibroacoustic problem and Illustrate its efficiency by means of numerical experiments. Keywords: finite element, boundary element, error estimation. 1 Introduction The solution of physical problems where the expected solution is a wave can be still considered an open problem of the numerical techniques such as the finite element method (FEM) or the boundary element method (BEM) [17]. If the wave length of the solution field is too short, the discretisation of the physical domain must be fine enough in order to describe this wave and the computational costs of the problem increase till they become unaffordable. Acoustics and structural vibration are among those physical problems whose solution is a wave. They are governed by the vibroacoustic equations. Several formulations are available [5,14]. However, the discussion and the examples shown here are based on the steady harmonic pressure-displacement equations: INTERNOISE 2010 │ JUNE 13-16 │ LISBON │ PORTUGAL 2 (1) Where p is the phasor of acoustic pressure, k is the air wave number, the acoustic force term is due to point sources, v n is an imposed velocity, ρ 0 is the air density, n are outward normals and linear elasticity is assumed for the solid part of the problem. The numerical formulation of this problem is reviewed in [1]. Examples of the numerical resolution of these equations in order to model physical situations of industrial interest can be found in the literature. The sound transmission through walls has been studied by means of the FEM in [11] or BEM in [3]. Studies on the numerical error for the uncoupled (acoustic or structural) problem have been done. In the field of numerical acoustics, a usual rule of thumb is that six linear finite elements per wave length are enough in order to obtain accurate results. More detailed analyses [9,10]) predict that the interpolation error for linear elements can be estimated as (2) where and are dimensionless wave number and element size, is a characteristic length of the problem, C local is a constant that depends on each particular situation but it is independent of and , p is the exact solution and p I is an interpolation of p. The six-elements-per-wave-length rule of thumb only takes into account this local interpolation error. Nevertheless, two additional phenomena make this criterion insufficient: the dispersion effect (k-singularity) and the existence of eigenfrequencies of the equation (λsingularity) [2]. The wave number of the numerical solution has been proved to be different from the exact wave number in the Helmholtz equation (for a linear one-dimensional finite element solution of Helmholtz equation we have . This causes an increasing phase shift of the discrete solution. The error affects all the domain. Due to dispersion, a numerical solution obtained by keeping constant would have more error for high INTERNOISE 2010 │ JUNE 13-16 │ LISBON │ PORTUGAL 3 frequencies. The increase of error due to the increase of dimensionless wave number is known as pollution effect [4]. Introducing the numerical wave number (affected by dispersion error) in the analysis, an expression for the total error can be obtained as (3) p h is the numerical solution obtained by FEM. C 1 and C 2 are constants that should be calculated for each particular situation (but they are independent of and ) and q the degree of polynomial interpolation. This is an important drawback because these a priori error estimates give a tendency but the mesh should be correctly designed for each problem. Note that in equation (3) we can associate the first term (multiplied by C 1 ) with the interpolation or local error (e I = (p-p I /p)) and the second term (multiplied by C 2 ) with the pollution error (e POL = (p I -p h )/p h ). Numerical experiments that show the evolution of numerical errors in the BEM for the Helmholtz equation can be found in [12,13]. In [16] an a posteriori error estimator for the case of structural dynamics in the time-domain has been developed. The work in the present contribution is focused in the analysis of possible sources of numerical error due to the fluid-structure coupling in the vibroacoustic problem. In all the previous references, the study is done in uncoupled acoustic or solid problems and uniform domains. Thus, the waves of the problem (pressure or displacement waves) depends only on the physical properties of the medium and the governing equation. The wavelength of the solution fields is provided by the dispersion relation. In a bounded problem, the solution tends to be similar to the eigenmode whose eigenfrequency is closer to the problem pulsation. However, the solution of the problem can also be composed of forced waves. They appear when the solution is required to fit irregularities at the boundary of the system, or satisfy some imposed force. In general, they can be different of any characteristic mode of oscillation of the system. This possibility is not considered in uncoupled system error analysis because for usual forces acting on the problem (point forces, uniform loads,...) forced waves are not generated. This is not the case of coupled vibroacoustic problems where an spatial oscillatory pressure is often the excitation force on a structure or a sinusoidal velocity is imposed on an acoustic contour. If the discretisation of a medium is designed taking as reference the wave length of resonant waves and the solution is composed of forced waves of shorter wave length the numerical errors would be larger than expected. A similar idea is found when modal analysis is used for sound transmission problems. In [6] both resonant and forced modes must be considered in order to accurately describe the pressure fields in the rooms. 2 Structural problem The dynamic response of the two structural systems in Figure 1 is analysed. The first one, Figure 1(a), is a simply supported beam with a point force placed at from the support ( is the length of the beam). It is a reference case where the solution will be only composed of INTERNOISE 2010 │ JUNE 13-16 │ LISBON │ PORTUGAL 4 resonant waves. On the contrary, in the second case (Figure 1(b)) a sinusoidal pressure is imposed ( ). The vibration field can be composed of resonant and forced waves, depending on the relationship between the geometric and mechanical properties of the beam, the pulsation of the problem and the force parameters K f and θ . Figure 1 – Sketches of the two structural problems analysed: (a)Beam loaded with a point force at the position (b)Beam loaded with a sinusoidal load . In all the examples shown in the contribution, the coupling and forcing terms are calculated with enough accuracy. As shown in [8] for static beams loaded with sinusoidal forces, the incorrect integration of the force term can be a source of error, especially if the force is integrated using the same accuracy that has been used for the mass matrices (typically linear or quadratic polynomials for structural analysis FEM). In any case, this possible source of error in vibroacoustic coupled problems will not be analysed in the contribution. The error measure in this structural problem is (4) u h is the phasor of the numerical value of the structural displacement for a mesh of size h, while u is the reference value. It is calculated by means of the modal analysis solution. The material properties of the structure are shown in Table 1. The calculated numerical errors can be seen in Figures 2 and 3. The point force example should be understood as the reference for the convergence behaviour. The interpolation of the normal displacement in an Euler Beam element is of degree q = 3 and the slope of the log 10 (||e||) - log 10 (h) curve is q +1 = 4. The error is clearly related with the structural wave number and increases with frequency as shown in Figure 1. The expression of the structural wave number for an Euler beam is (5) Table 1 – Geometrical and mechanical properties of the structure. Meaning Symbol Value Thickness 5 mm Young’s modulus 2,061·10 11 N/m 2 Density 7500 kg/m 3 Damping 2 % Length 4 m INTERNOISE 2010 │ JUNE 13-16 │ LISBON │ PORTUGAL 5 Figure 2 – Error analysis for a point loaded beam: on the left, dependence on the number of elements per wave length; on the right dependence on the frequency. In the sinusoidal load case of Figure 3 the value of the `force wave number' K f = 100 π / , has been chosen in order to perfectly match the excitation force with the mode that has 50 wave lengths inside the beam. In that case the exact solution is composed by only one forced wave. Figure 3 – Error analysis for a beam loaded by means of an oscillatory force where K f = 100 π / (exactly 50 load waves per beam length) : (a)relative error depending on the number of nodes per structural wave length; (b)relative error depending on the number of nodes per pressure load wave length. In Figure 3 (a) the error has been plotted depending on the structural wave number and in Figure 3 (b) depending on the force wave number K f (number of structural elements per imposed pressure wave length). The results take sense in the second case because the INTERNOISE 2010 │ JUNE 13-16 │ LISBON │ PORTUGAL 6 displacement field wave length is K f . This means that the mesh must be designed by considering the excitation wave number K f instead of the structural wave number. 3 Acoustic problem The effect of spatial oscillatory excitations can be also found in the acoustic problem. The acoustic radiation in a semi-infinite semi-space has been calculated by means of boundary elements (a similar situation has been analysed in [15]). An imposed velocity v n = K v y = 2 π n/L y has been imposed along a 4 m length vibrating boundary. The error measure in this acoustic problem is (6) where p h is the phasor of the numerical value of the acoustic pressure for a mesh of size h, while p ref is the reference value. This value is now taken as the pressure in a mesh with smaller elements. The integration of the outputs is performed along the 4 m length contour with boundary elements (radiating part of the contour). Figure 4 – Convergence of the relative error in acoustic pressure for the acoustic radiation of sound inside a semi-infinite space. Constant frequency of the problem (250 Hz), with increasing value of the spatial wave number of the imposed velocity at the contour. The obtained results are shown in Figure 4. We can see how the convergence of the numerical error depends again on the wave number of the imposed velocity. For the case of n = 15 waves inside the computational domain the expected convergence slope is obtained. INTERNOISE 2010 │ JUNE 13-16 │ LISBON │ PORTUGAL 7 However, for larger number of waves inside the computational domain (n = 30 and n = 50), the convergence with the expected slope starts for a smaller element size. At less, more than one element per excitation wave length is required in order to reach the convergence of order 2. In these cases the excitation wave length is shorter than the length of propagation waves in the air. The wave number of the excitation oscillation is smaller than the wave number of air (k = ω /c). In this situation the pressure field is mainly composed of evanescent waves. This waves decrease with distance and vanish away from the source location (vibrating contour in this case). 4 Conclusions It has been shown that for both structural and acoustic problems the convergence of numerical methods such as FEM or BEM can depend on the excitation action. In these methods, the spatial discretisation of the physical domain must be designed according to the expected waves in the solution. In general the length of these waves can be predicted by considering only the physical properties of the medium (resonant or propagating waves). However for some excitation actions forced waves can be induced. This causes the final solution to have contributions related to different frequencies and consequently modifies the expected performance of the considered numerical methods. This behaviour in uncoupled problems should be studied for vibroacoustic problems in the following cases: • Acoustic domains that are in contact with structures vibrating with a wavelength smaller than the length of the waves generated in the air. It is the case of acoustic domains coupled with lightweight structures and frequencies below the critical frequency. The discretisation in the acoustic part of the problem should be done by taking into account the structural wave length. • Structures surrounded by acoustic fluids and frequencies above the critical frequency. The discretisation in the structural part of the problem should be done by taking into account the acoustic wave length. According to the presented preliminary results, it can be concluded that the element size should be decided not only by taking into account the medium properties but also the excitation forces of the problem. In addition, it will be checked if for coupled vibroacoustic problems this effect is also relevant in the coupling interface and the shortest wave length of any of the media involved in the problem should be considered. Acknowledgments The financial support of the Ministerio de Educación y Ciencia (BIA2007-66965, DPI200762395) and the Facultat de Matemàtiques i Estadística de la UPC (FME, CERMET) is gratefully acknowledged. Free software has been used [7]. INTERNOISE 2010 │ JUNE 13-16 │ LISBON │ PORTUGAL 8 References [1] Atalla,N.; Bernhard, R.J. Review of numerical solutions for low-frequency structuralacoustic problems, Appl. Acoust., Vol 43, 1994, pp 271-294. [2] Bouillard,Ph.; Ihlenburg,F. Error estimation and adaptivity for the finite element method in acoustics: 2D and 3D applications, Comput. Methods Appl. Mech. Engrg., Vol 176(1-4), 1999, pp 147-163. [3] Coyette, J.P. The use of finite-element and boundary-element models for pre-dicting the vibro-acoustic behaviour of layered structures, Adv. Eng. Softw., Vol 30, 1999, pp 133139. [4] Deraemaeker, A.; Babuska,I.; Bouillard,Ph. Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions, Int.J. Numer. Meth. Engng., Vol 46, 1999, pp 471-499. [5] Everstine, G.C. Finite element formulations of structural acoustics problems, Comput. Struct., Vol 65(3), 1997, pp 307-321. [6] Gagliardini, L.; Roland, J.; Guyader,J.L. The use of a functional basis to calculate acoustic transmission between rooms, J. Sound Vibr., Vol 145(3), 1991, pp 457-478. [7] Geuzaine, C.; Remacle, J.-F. Gmsh: a three-dimensional finite element mesh generator with built-in preand post-processing facilities, Int. J. Numer. Meth. Engng., Vol 11(79), 2009, pp 1309-1331. [8] Gudla, P.K.; Ganguli,R. Error estimates for inconsistent load lumping ap-proach in finite element solution of differential equations, Applied Mathematics and Computation, Vol 194, 2007, pp 21-37. [9] Ihlenburg,F. Finite element analysis of acoustic scattering. Springer, 1998. [10] Ihlenburg,F. ; Babuska,I. Finite element solution of the Helmholtz equationwith high wave number part II: the h-p version of the FEM, SIAM J. Numer.Anal., Vol 34(1), 1997, pp 315-358. [11] Maluski,S.; Gibbs, B.M. Application of a finite-element model to low-frequency sound insulation in dwellings, J. Acoust. Soc. Am., Vol 108(4), 2000, pp 1741–1751. [12] Marburg, S. Six elements per wavelength. is that enough?, J. Comput. Acoust., Vol 10, 2002, pp 25–51. [13] Marburg, S.; Schneider, S. Influence of element types on numeric error for acoustic boundary elements, J. Comput. Acoust., Vol 11, 2003, pp 363-386. [14] J.P. Morand and R. Ohayon. Interactions fluides-structures. Masson, Paris, 1992. [15] Poblet-Puig,J.; Rodrguez-Ferran,A.; Guigou-Carter,C.; Villot,M. Numerical modelling of the radiation efficiency of asymmetrical structures, Appl. Acoust., Vol 70(5), 2009, pp 777–780. [16] Wiberg, N.; Zeng,L.; Li,X. Error estimation and adaptivity in elastodynamics, Comput. Methods Appl. Mech. Engrg., Vol 101, 1992, pp 369-395. [17] Zienkiewicz, O.C. Achievements and some unsolved problems of the finite element method, Int. J. Numer. Meth. Engng., Vol 47, 2000, pp 9-28.