The block Gauss-Seidel method in sound transmission problems

THE BLOCK GAUSS-SEIDEL METHOD IN SOUND TRANSMISSION PROBLEMS J. Poblet-Puig∗ , A. Rodr´ıguez-Ferran† Laboratori de C`alcul Num`eric E.T.S. d’Enginyers de Camins, Canals i Ports de Barcelona Universitat Polit`ecnica de Catalunya May 17, 2011 Abstract Sound transmission through partitions can be modelled as an acoustic fluidelastic structure interaction problem. The block Gauss-Seidel iterative method is used in order to solve the finite element linear system of equations. The blocks are defined, respecting the fluid and structural domains. The convergence criterion is analysed and interpreted in physical terms by means of simple onedimensional problems. This analysis highlights the negative influence on the convergence of a strong degree of coupling between the acoustic domains and the structure. A selective coupling strategy has been developed and applied to problems with strong coupling (e.g. double walls). 1 Introduction Fluid-structure interaction is the key aspect of many acoustic problems of practical interest. This is the case, for instance, of sound propagation through partitions in buildings. In this field, vibroacoustic equations are solved in the frequency domain and multiple frequencies have to be considered in order to obtain the frequency response spectrum of the system. Several formulations of the vibroacoustic equations are available in the literature [1–3]. The structure is usually described by means of the displacement formulation while several options have been used for the fluid. When the fluid is described by means of displacement or a mixed pressure-displacement formulations the coupled problem is globally symmetric. However there exist spurious modes. Modifications ∗correspondence: UPC, Campus Nord B1, Jordi Girona 1, E-08034 Barcelona, Spain, e-mail: [email protected] †e-mail: [email protected] 1 Poblet-Puig, J. and Rodríguez-Ferran, A, The block Gauss-Seidel method in sound transmission problems, Journal of Computational Acoustics, Vol. 18, Issue 1, pp. 13-30, 2010 have to be introduced in order to avoid this drawback, see for example [4]. This is not the case of the velocity potential [3] or pressure formulations [5–7]. In the first case the systems of equations to be solved are symmetrical but an extra postprocessing is required if the final output is the acoustic pressure (which is often the case in building acoustics). In the pressure formulation the systems of equations are unsymmetrical but smaller than in displacement or mixed formulations. In all of these cases the finite element discretisation leads to systems of linear equations with a block structure. The diagonal blocks in the global matrix are typically symmetric and indefinite, but the off-diagonal blocks (which represent the coupling between the acoustic fluid and the elastic structure) break the symmetry of the global matrix. For this reason, a monolithic solution approach requires the use of general solvers for unsymmetrical and indefinite matrices, such as Crout factorisation or GMRES iterations with an appropriate preconditioner [8,9]. Alternatively, block iterative solvers can be used. By doing so, the symmetry of the diagonal blocks can be exploited, and the storage requirements are decreased. The block Gauss-Seidel iterative solver is considered here. The well-known convergence condition (spectral radius of iteration matrix smaller than one) is interpreted from a physical viewpoint, by considered simple, one-dimensional vibroacoustic models. This analysis shows the detrimental effect on the convergence of the iterative solver of i) the excitation frequency being close to an acoustic or structural eigenfrequency and ii) a large level of coupling between the acoustic fluid and the structure. The coupling in a fluid-structure problem can be weak or strong depending on different factors such as the physical properties of the media, the geometry of the structure or the imposed boundary conditions. This has important consequences in the performance of the solving strategies and has been widely studied when the problem is solved in the time domain (see for example [10–14]). Computational costs can be drastically reduced if the coupling is weak and the fluid and structure problems can be solved in a staggered way. In the frequency domain, it is not easy to predict when the coupling is strong or weak. In the first case a monolithic solver is required while in the second case a one-way coupling strategy is enough. The block Gauss- Seidel method considered here is between these two strategies, taking advantage of weak coupling when it is possible and reducing computational costs. The block Gauss-Seidel method can be understood as a domain decomposition method. These methods base their efficiency in the splitting of the physical domain of the problem into smaller subdomains. The system of equations is then solved at two different levels. On the one hand each subdomain is solved as an individual problem and on the other hand the continuity between them is imposed. These techniques have been mainly designed to be used in parallel computing machines. Each CPU deals with a single smaller domain using the more adequate solver for each region. In [15–17] domain decomposition techniques have been used in order to solve scattering problems governed by the Helmholtz equations in big physical domains. The continuity of the pressure field (and its normal derivative) in the interface between regions is imposed by means of Lagrange multipliers. Moreover each subdomain has to be regularised by means of fictitious boundary conditions in order to avoid problems caused by artificial eigenfrequencies. Domain decomposition techniques has also 2 been used for vibroacoustic problems. In [18] the partitions have been done in both the acoustic domains and the structure. Finally, in [19–21], block Jacobi and block Gauss-Seidel algorithms have been used for vibroacoustic problems where the decomposition of the domain strictly respects the physical regions (fluid and structure). The only interface between subdomains is the fluid-structure boundary. Since the goal in [19–21] is to propose a general solver (also for strongly coupled problems) their discussion is focused on the convergence of the methods. It seems clear that using the physical interface conditions in order to transfer information between fluid and structural subdomains leads to divergence in a large number of situations. They propose relaxed coupling conditions that cause the block Gauss-Seidel algorithm to have fast convergence for all the analysed situations. In the application examples the performance of the modified algorithms around the eigenfrequencies of the problem has not been analysed. Moreover, the use of the modified interface conditions require some modifications at finite element level. An outline of the paper follows. The convergence condition and its physical interpretation are covered in Section 3. The application examples of Section 4corroborate this interpretation, and motivate the selective coupling strategy presented in Section 5, which is applied to the problem of sound propagation through double walls. The concluding remarks of Section 6close the paper. 2 The block Gauss-Seidel algorithm The block Gauss-Seidel algorithm will be presented in matrix form. The coupled system of linear equations is F CF S CSF SxF xS=fF fS(1) The pressure-displacement formulation is taken as reference here. Fis the flexibility matrix governing the fluid domain with nodal unknowns xF(typically pressures) and Sis the stiffness matrix governing the structural domain with nodal unknowns xS (typically displacements and rotations). If FEM is used, Fand Sare typically sparse, symmetric and indefinite matrices. fFand fSare the forces acting in the fluid and structural domains. The coupling is taken into account by means of matrices CF S and CSF . The forces acting on the structure due to the acoustic pressures in the fluid are fSF =CSF xF(2) and the acoustic forces in the fluid contour caused by the structural vibrations are fF S =CF SxS(3) The global matrix of Equation (1) is non-symmetric for the more widely used formulations. The block Gauss-Seidel algorithm is summarised in Table 1. The initial guess can be chosen as the solution of the uncoupled problems. The convergence is checked by means of the relative errors in the solution e(i) F=x(i) F−x(i+1) F x(i+1) F ;e(i) S=x(i) S−x(i+1) S x(i+1) S (4) 3 Table 1: The block Gauss-Seidel method Choose an initial guess x(0) S,x(0) F for i= 0,1,2,... Fx(i+1) F=fF−CF Sx(i) S Sx(i+1) S=fS−CSF x(i+1) F check convergence; continue if necessary end and the relative residual r(i) F=Fx(i) F+CF Sx(i) S−fF fF (5) The two systems of equations in Table 1have to be solved several times with different force vectors but constant matrices Fand S. This has to be exploited for maximum efficiency. A first option is to use a direct solver (for small matrix dimensions) and save the factorisation of the matrices. Another possibility is to use the adequate iterative solver (GMRES, MINRES,... see [8,9] for more details and [22–24] for robust implementations) and save the preconditioner, which is calculated only once for i= 0 and can be reused for the successive iterations. Wave problems often require to perform calculations for successive frequencies or different types of force terms. Matrices are then very similar. The possibility of using the same preconditioner for several successive frequencies has also to be considered. 3 Analysis of the block Gauss-Seidel method 3.1 The convergence condition As other stationary iterative methods, the block Gauss-Seidel algorithm converges if the spectral radius ρ(i.e. the maximum modulus of the eigenvalues) of the iteration matrix Gis less than one, see [9]. The algorithm in Table 1can be rewritten as (x(i+1) F x(i+1) S)=0−F−1CF S 0 S−1CSF F−1CF S (x(i) F x(i) S)+F−1fF S−1(fS−CSF F−1fF)(6) The iteration matrix Gis the matrix in Equation (6), so the convergence condition is ρS−1CSF F−1CF S<1 (7) 3.2 Physical interpretation of the convergence condition The simplified model of Figure 1will be used in order to understand and illustrate the phenomena of vibroacoustic coupling and the performance of the block Gauss- Seidel algorithm. A vibrating mass is coupled with an acoustic domain. Both can be 4 x K ΩF M u(t) vn ρ a, c F(t) l Figure 1: Simple one-dimensional coupled system with two degrees of freedom. excited: the mass by means of an exterior force F(t) = Re ϕeiωtand the acoustic fluid cavity by an exterior imposed velocity vn n n. Note that the model is formulated for a unit surface. Thereby ϕis the phasor of force per unit surface, and Mand Kthe mass and stiffness per unit surface. The interaction between the acoustic fluid and the single mass can be characterised by the pressure applied by the fluid on the mass and the displacement imposed by the mass at the acoustic contour. The governing equation and boundary conditions for the fluid domain are d2p(x) d x2+k2p(x) = 0 xin ΩF(8) d p(x) d x x=0 =ρFω2u(9) d p(x) d x x=ℓ=−ρFiωvn n n(10) and if the frequency of the problem is a real value, the pressure field is p(x) = C1cos(kx) + C2sin(kx) (11) where C1and C2are unknown complex constants. Taking into account the dynamic equilibrium of the single mass, a linear system with three equations results:   sin(kℓ)−cos(kℓ)0 0 1 −ρFωc 1 0 K−ω2M    C1 C2 u   =   ρFicvn n n 0 ϕ   (12) This system is the particularisation for this simple one-dimensional example of Equation (1). The convergence condition (7) leads to ρ(G) = cos (kℓ)ρFωc sin (kℓ) (K−ω2M)<1 (13) A similar analysis can be done for the one-dimensional model with two fluid domains shown in Figure 2and [25], which represents two rooms separated by a partition. The expression of the spectral radius is ρ(G) = ρFωc K−ω2Mcos (kℓ1) sin (kℓ1)+cos (kℓ2) sin (kℓ2)(14) 5 x x K M Ω1Ω2 l1l2 1 2 Figure 2: Simple one-dimensional coupled system with two acoustic domains. Several conclusions can be obtained from Equations (13) and (14). The method is less efficient for denser fluids (i.e. larger density ρF) or fluids with higher wave speed c, because the coupling between the structure and the fluid increases. The geometry of the problem is also important. For this one-dimensional case, the geometry is represented by terms cos (kℓ) and sin (kℓ) (ℓ=ℓ1or ℓ2in Equation (14)). If sin(kℓ)≈0 the method will not converge. This happens for the eigenfrequencies of the acoustic cavity, kn=nπ/ℓ, but also when ℓis very small (small fluid domains). If cos(kℓ) = 0 the method converges in one iteration. This is a very specific situation of the one-dimensional model and cannot be generalised to higher dimensions. Note that the method also diverges for frequencies close to the structural eigenfrequency pK/M. Finally, Equations (13) and (14) also show that the performance of the iterative solver improves with the frequency. A similar parameter (λ=ρFc/ρStω) has been defined by [26]. t is the typical thickness of the structure and ρSits density. However, λdoes not take into account the influence of the geometry nor the stiffness. By condensing out the unknown C2and noting from Equation (11) that C1is p(x= 0), system (12) can be recast as sin(kℓ)−ρFωc cos(kℓ) 1K−ω2Mp(x= 0) u=ρFicvn n n ϕ(15) and ρ(G) can then be viewed as the ratio of stiffness of the fluid and the structure, including the effect of coupling: ρ(G) = (d p(x= 0)/d u)F (d p(x= 0)/d u)S (16) The subscript F(S) means here derivative from the point of view of the fluid (structure). The conceptual behaviour of the fluid-structure system has been plotted in Figure 3. The harmonic equilibrium is reached at the pair x∗ S−x∗ F. The acoustic pressure caused by the structural displacement is xF S =−F−1CF Sx∗ S(17) and the displacement caused by the pressure is xSF =−S−1CSF x∗ F(18) 6 xF * xS xF unc. xS unc. xS xFS xSF * 1 xF1 FLUIDSTRUCTURE −1 −F CFS −S C −1 SF Figure 3: Conceptual behaviour of a coupled fluid-structure system. Figure 4shows a sketch of the convergence and divergence of the algorithm depending on the spectral radius. G ρ( ) <1 S x(0) S x(2) F x(1) F x(2) xS xF S x(1) F x * S x * STRUCTURE FLUID (a) S x(0) S x(1) F x(1) G ρ( ) >1 xS xF F x * S x * FLUID STRUCTURE (b) Figure 4: Convergence of the block Gauss-Seidel algorithm: (a) Convergence for ρ(G)<1; (b) Divergence for ρ(G)>1. 7 4 Application examples The performance of the block Gauss-Seidel method is illustrated here with various 2D and 3D vibroacoustic problems (sound transmission through a single wall). The method has already been used to study the sound transmission through single and double walls in [25,27]. The goals are 1) to show the influence on the convergence of the iterative solver of i) the damping, ii) the coupling between fluid and structure and iii) acoustic and structural eigenfrequencies, and 2) to demonstrate the use of the solver in practical simulations. A FEM-FEM approach (i.e. finite elements for the fluid and acoustic domains) is used here. 4.1 Influence of damping Two acoustic domains are separated by a single wall (represented in this two-dimensional setting by a concrete beam), see Figure 5. The acoustic excitation is a punctual sound source placed in the left bottom corner of the first domain, at a distance of 0.5 m to the contours. The room dimensions are 3 ×3 m2and 4 ×3 m2. The material and geometrical parameters are summarised in Table 2. Note that we are dealing with air, which is a very light fluid. This is the typical situation where the method will have a very good behaviour. A relative tolerance of 10−9is used in the stopping criteria defined in Equations (4) and (5). 01 0011 1 m Sound source Robin or pure reflecting boundary Structural element (Euler beam) Figure 5: Sound transmission through a single wall Two different situations have been analysed. On the one hand, an undamped problem (no acoustic absorption and no structural damping). On the other hand, the same problem with an acoustic absorption of 30 % at the boundaries (introduced by means of a Robin boundary condition) and hysteretic structural damping (5 %). These are reasonable values. The results (number of iterations required) have been plotted in Figure 6. Note that damping considerably decreases the number of iterations required, especially near eigenfrequencies. 8 Table 2: Material and geometrical data for the acoustic domains (air), the heavy single wall and the lightweight leaves of the double wall. STRUCTURE Meaning Symbol Heavy Lightweight Young’s modulus E2.943 ·1010 N/m24.5·109N/m2 Poisson’s ratio ν0.25 0.25 Wall density ρS2500 kg/m3913 kg/m3 Wall thickness t 0.10 m 0.013 m Loss factor η0−5 % 0 −5 % FLUID Meaning Symbol Value Speed of sound c 340 m/s Density of fluid ρF1.18 kg/m3 Source strength Q0.005i m3/s Acoustic absorption α0−30 % 0 10 20 30 40 50 20 40 60 80 100 120 140 160 180 200 f (Hz) Iterations Eig. ΩF(1) Eig. ΩF(2) Eig. Structure (a) 0 10 20 30 40 50 20 40 60 80 100 120 140 160 180 200 f (Hz) Iterations Eig. ΩF(1) Eig. ΩF(2) Eig. Structure (b) Figure 6: Iterations of the block Gauss-Seidel solver: (a) undamped problem; (b) damped problem (30 % acoustic absorption and 5 % structural damping). Eigenfrequencies of the sending and receiving domains, Ω(1) Fand Ω(2) F, and the structure are also shown. 4.2 Influence of particular eigenfrequencies The performance of the block Gauss-Seidel solver for three particular frequencies has been analysed. The example of Section 4.1 has been considered (damped situation). The aim of the analysis is to show differences in the efficiency of the method depending on the type of eigenfrequencies that are close to the excitation frequency. The studied frequencies are: i) 70 Hz, which is close to uncoupled eigenfrequencies of the structure (70.83 Hz) and the receiving room (69.27 Hz); ii) 90 Hz, which is not close to any of the eigenfrequencies of the problem; iii) 156 Hz, which is close to an uncoupled eigenfrequency of the structure (156.22 Hz). 9 0 20 40 60 80 100 50 100 150 200 250 f (Hz) Selective Coupling Standard BGS Eig. ΩF(1) Eig. ΩF(2) Eig. Leaves Nonconvergence (a) 0 20 40 60 80 100 50 100 150 200 250 f (Hz) Selective Coupling Standard BGS Eig. ΩF(1) Eig. ΩF(2) Eig. Leaves Nonconvergence (b) Figure 12: Application of the selective coupling to a problem of sound transmission through a lightweight double wall. Iterations required for each algorithm: (a) air cavity; (b) absorbing material. Tolerance: = 10−9; maximum number of iterations: 100. Seidel iteration with the appropriate ordering of blocks, along sound trajectory) is not accurate enough, and an iterative procedure to reach convergence is required. For larger degrees of coupling, these iterations may fail to converge. This is the case, for instance, in the numerical simulation of sound transmission through double walls: the two leaves and the small acoustic cavity between them are strongly coupled. This observation has suggested a selective coupling strategy, where the structure and the problematic acoustic domains (e.g. the cavities in a double wall) are treated 16 together, in the same block. The convergence analysis also shows the negative effect of an excitation frequency close to an (acoustic or structural) eigenfrequency for undamped problems. However, the convergence improves significantly if one uses realistic values of structural damping and acoustic absorption. One-, two- and three-dimensional numerical examples have been used to illustrate the capabilities of the block Gauss-Seidel solver. The systematic use of this solver for the numerical simulation of sound transmission through single and double has been reported in [25,27]. Acknowledgements The financial support of the Fons Social Europeu (2003 FI 00652), Ministerio de Educaci´on y Ciencia (BIA2007-66965, DPI2007-62395) and the Research Fund for Coal and Steel (RFSR-CT-2003-00025) is gratefully acknowledged. Free software has been used [22,23,28,30]. References [1] J.P. Morand and R. Ohayon. Fluid-structure interaction. Wiley, 1995. [2] J. Wandinger. A symmetric Craig-Bampton method of coupled fluid-structure systems. Eng. Comput., 15(4):450–461, 1998. [3] G.C. Everstine. Finite element formulations of structural acoustics problems. Comput. Struct., 65(3):307–321, 1997. [4] K.J. Bathe, C. Nitikitpaiboon, and X. Wang. A mixed displacement-based finite element formulation for acoustic fluid-structure interaction. Comput. Struct., 56(2):225–237, 1995. [5] F. Ihlenburg. Finite element analysis of acoustic scattering. Springer, 1998. [6] P. Bouillard, V. Lacroix, and E. De Bel. A wave-oriented meshless formulation for acoustical and vibro-acoustical applications. Wave Motion, 39(4):295–305, 2004. [7] J.P. Tuck-Lee and P.M. Pinsky. Adaptive frequency windowing for multifrequency solutions in structural acoustics based on the matrix Pad´e-via-Lanczos algorithm. Int. J. Numer. Methods Engng., 73(5):728–746, 2008. [8] R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, and H. Van der Vorst. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, 2nd Edition. SIAM, Philadelphia, PA, 1994. [9] Y. Saad. Iterative methods for sparse linear systems. http://wwwusers.cs.umn.edu/˜saad/books.html, 2000. 17 [10] C.A. Felippa, K.C. Park, and Ch. Farhat. Partitioned analysis of coupled mechanical systems. Comput. Meth. Appl. Mech. Eng., 190(24–25):3247–3270, 2001. [11] Q. Zhang and T. Hisada. Studies of the strong coupling and weak coupling methods in FSI analysis. Int. J. Numer. Methods Eng., 60(12):2013–2029, 2004. [12] H.G. Matthies and J. Steindorf. Partitioned strong coupling algorithms for fluidstructure interaction. Comput. Struct., 81(8–11):805–812, 2003. [13] S. Piperno, C. Farhat, and B. Larrouturou. Partitioned procedures for the transient solution of coupled aroelastic problems Part I: Model problem, theory and two-dimensional application. Comput. Meth. Appl. Mech. Eng., 124(1–2):79–112, 1995. [14] S. Piperno and C. Farhat. Partitioned procedures for the transient solution of coupled aeroelastic problems – Part II: energy transfer analysis and three-dimensional applications. Comput. Meth. Appl. Mech. Eng., 190(24–25):3147–3170, 2001. [15] C. Farhat, A. Macedo, and M. Lesoinne. A two-level domain decomposition method for the iterative solution of high frequency exterior Helmholtz problems. Numer. Math., 85(2):283–308, 2000. [16] C. Farhat, A. Macedo, M. Lesoinne, F Roux, F. Magouls, and A. de La Bourdonnaie. Two-level domain decomposition methods with Lagrange multipliers for the fast iterative solution of acoustic scattering problems. Comput. Meth. Appl. Mech. Eng., 184(2–4):213–239, 2000. [17] R. Tezaur, A. Macedo, and Ch. Farhat. Iterative solution of large-scale acoustic scattering problems with multiple right hand-sides by a domain decomposition method with Lagrange multipliers. Int. J. Numer. Methods Eng., 51(10):1175– 1193, 2001. [18] J. Mandel. An iterative substructuring method for coupled fluid-solid acoustic problems. J. Comput. Phys., 177(1):95–116, 2002. [19] X. Feng. Interface conditions and non-overlapping domain decomposition methods for a fluid-solid interaction problem. Contemporary Mathematics, 218:417–424, 1998. [20] P. Cummings and X. Feng. Domain decomposition methods for a system of coupled and elastic Helmholtz equations. In Ch. Lai, P.E. Bjorstad, M. Cross, and O. Widlund, editors, 11th International conference on domain decomposition methods (www.ddm.org), Greenwich (England), 1998. [21] X. Feng and Z. Xie. Domain decomposition methods for a coupled vibration between an acoustic field and a plate. In T. Chan, T. Kako, H. Kawarada, and O. Pironneau, editors, 12th International conference on domain decomposition methods (www.ddm.org), Chiba (Japan), 1999. 18 [22] Satish Balay, Kris Buschelman, William D. Gropp, Dinesh Kaushik, Matthew G. Knepley, Lois Curfman McInnes, Barry F. Smith, and Hong Zhang. PETSc Web page, 2001. http://www.mcs.anl.gov/petsc. [23] Satish Balay, Kris Buschelman, Victor Eijkhout, William D. Gropp, Dinesh Kaushik, Matthew G. Knepley, Lois Curfman McInnes, Barry F. Smith, and Hong Zhang. PETSc users manual. Technical Report ANL-95/11 - Revision 2.1.5, Argonne National Laboratory, 2004. [24] Satish Balay, William D. Gropp, Lois Curfman McInnes, and Barry F. Smith. Efficient management of parallelism in object oriented numerical software libraries. In E. Arge, A. M. Bruaset, and H. P. Langtangen, editors, Modern Software Tools in Scientific Computing, pages 163–202. Birkh¨auser Press, 1997. [25] J. Poblet-Puig. Numerical modelling of sound transmission in lightweight structures. PhD thesis, Universitat Polit`ecnica de Catalunya (UPC), January 2008. [26] N. Atalla and R.J. Bernhard. Review of numerical solutions for low-frequency structural-acoustic problems. Appl. Acoust., 43(3):271–294, 1994. [27] J. Poblet-Puig, A. Rodr´ıguez-Ferran, C. Guigou-Carter, and M. Villot. The role of studs in the sound transmission of double walls. Acta Acust. United Acust., 95(3):555–567, 2009. [28] C. Geuzaine and J. Remacle. Gmsh reference manual. http://www.geuz.org/gmsh/, 2009. [29] F. Fahy. Sound and structural vibration. Academic Press, London, 1989. [30] EDF. Code-aster home page. http://www.code-aster.org, 2009. 19