- Email: [email protected]
A computationally efficient prediction technique for the steady-state dynamic analysis of coupled vibro-acoustic systems
Advances in Engineering Software 33 (2002) 527–540 www.elsevier.com/locate/advengsoft
A computationally efficient prediction technique for the steady-state dynamic analysis of coupled vibro-acoustic systems W. Desmet*, B. van Hal, P. Sas, D. Vandepitte K.U. Leuven, Department of Mechanical Engineering, Division PMA, Celestijnenlaan 300 B, B-3001 Leuven (Heverlee), Belgium Received 14 November 2000; accepted 1 July 2002
Abstract A new prediction technique has been developed for the steady-state dynamic analysis of coupled vibro-acoustic systems. In contrast with the finite element method, in which the dynamic field variables within each element are expanded in terms of local, non-exact shape functions, the dynamic field variables are expressed as global wave function expansions, which exactly satisfy the governing dynamic equations. The contributions of the wave functions to the coupled vibro-acoustic response result from an integral formulation of the boundary conditions. This paper describes the basic concept of the new technique for the modelling of the vibro-acoustic coupling between the pressure field in an acoustic cavity with arbitrary shape and the out-of-plane displacement of a flat plate with arbitrary shape. It is illustrated through a threedimensional validation example that the new prediction technique yields a high accuracy with a substantially smaller computational effort than the finite element method, so that the new prediction technique can be applied up to much higher frequencies. q 2002 Civil-Comp Ltd and Elsevier Science Ltd. All rights reserved. Keywords: Numerical prediction technique; Coupled vibro-acoustic systems; Trefftz-method; Computational efficiency; Steady-state dynamic behaviour; Mid-frequency technique
1. Introduction At present, the finite element and boundary element methods are the most commonly used numerical prediction techniques for solving steady-state dynamic problems, defined in structural and acoustic continuum domains. An important implication of the conventional element concept, using non-exact shape functions, is that the large model sizes involved practically restrict the applicability of these deterministic prediction techniques to the low-frequency range. Above a certain frequency limit, these methods would require, even with the nowadays-available highperformance computer resources, a prohibitively large amount of computational effort and memory resources to get an acceptable level of prediction accuracy. For coupled vibro-acoustic problems, the conflict between computational effort and accuracy is even more pronounced than for uncoupled structural or uncoupled acoustic problems. First of all, the size of coupled prediction * Corresponding author. Tel.: þ 32-16-322480; fax: þ 32-16-322987. E-mail address: [email protected] (W. Desmet).
models is substantially larger, since a structural and an acoustic problem must be solved simultaneously to incorporate the acoustic pressure loading on the elastic structure and the continuity of the normal fluid and structural displacements along the fluid –structure coupling interface. Secondly, the numerical solution procedure for coupled vibro-acoustic models is less efficient than for uncoupled structural or uncoupled acoustic models, since coupled models, at least in their most commonly used acoustic pressure/structural displacement formulation [1], are no longer symmetric. Finally, the efficiency of the modal expansion method in reducing the model size is significantly reduced for coupled vibro-acoustic problems. The most appropriate mode selection for a modal expansion is the use of the modes of the coupled system. However, since coupled models are no longer symmetric, the calculation of coupled modes requires a non-symmetric eigenvalue calculation, which is very time consuming and which makes it a practically impossible procedure for many real-life vibroacoustic problems. The most commonly used alternative is a modal expansion in terms of uncoupled structural and uncoupled acoustic modes, which result from symmetric
0965-9978/02/$ - see front matter q 2002 Civil-Comp Ltd and Elsevier Science Ltd. All rights reserved. PII: S 0 9 6 5 - 9 9 7 8 ( 0 2 ) 0 0 0 6 2 - 5
528
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
and computationally efficient eigenvalue problems. However, the fact that uncoupled acoustic modes have a zero displacement component, normal to the fluid – structure coupling interface, implies that a large number of high-order uncoupled acoustic modes is required to accurately represent the normal displacement continuity along the fluid –structure interface. Hence, the benefit of a computationally efficient construction of the modal base is significantly reduced by the smaller model size reduction, obtained with an uncoupled modal base. As a result, there exists still a substantial discrepancy between the limited frequency range, in which conventional element based models can provide accurate coupled vibroacoustic predictions, and the significantly larger frequency range, in which accurate, deterministic predictions are needed for many industrial engineering problems. In view of narrowing this currently existing frequency discrepancy, a new deterministic prediction technique has been developed [2], which is based on the indirect Trefftz method [3,4]. The steady-state dynamic field variables in the entire—or at least in large subdomains of the—acoustic and structural domains of a coupled vibro-acoustic system are approximated in terms of a set of acoustic and structural wave functions, which are exact solutions of the homogeneous parts of the governing dynamic equations, and some particular solution functions of the inhomogeneous equations. In this way, the governing dynamic equations are exactly satisfied, irrespective of the contributions of the wave functions to the field variable expansions. Since the proposed field variable expansions are exact solutions of the governing dynamic equations, the contributions of the wave functions to the field variable expansions are merely determined by the acoustic and structural boundary conditions. Since only finite sized prediction models are amenable to numerical implementation, the boundary conditions can only be satisfied in an approximate way. Therefore, the wave function contributions result from a weighted residual formulation of the boundary conditions. A key issue in the use of the indirect Trefftz method is the definition of a T-complete function set, which ensures the convergence of the subsequent field variable expansions towards the exact solutions. Several T-complete function sets have already been defined for solving steady-state acoustic problems [5,6] and steady-state dynamic plate bending problems [7]. Although the theoretical convergence for these function sets has been proven, their practical convergence is, however, seriously disturbed or even prevented, due to the ill-conditioning of the involved model matrices. These numerical condition problems may be circumvented, to some extent, by subdividing the considered continuum domain into small elemental subdomains. This has led to the development of the T(refftz)-element approach [8], which allows the introduction of the Trefftz-
idea into a standard finite element scheme: the internal field variables within the T-elements are approximated in terms of a suitably truncated non-conforming T-complete set of trial functions, satisfying the governing equations a priori, while the boundary conditions and interelement continuity are enforced in an average integral sense. Two main approaches may be distinguished. In the hybrid T-element approach, an auxiliary frame is defined at the element boundaries. A conforming frame field distribution in terms of conventional (polynomial) trial functions is defined on the boundary frame and is used to enforce the boundary conditions and interelement conformity. The least-squares T-element approach is based only on the contributions of the T-complete functions to the field variable distributions within the element, without the introduction of auxiliary frame degrees of freedom. Zielinski and Zienkiewicz [9] indicated how the interelement conformity and the boundary conditions may be enforced directly in a least-squares sense. Jirousek and Wroblewski [10] translated this approach into an element formulation, which is suitable for implementation into a standard finite element code. One of the advantages of the T-element approach is its ability to include special-purpose functions to the internal element solution expansions, so that singularity problems can be handled efficiently without troublesome local mesh refinement [11]. However, in its current stage of development, the efficiency and accuracy of the T-element approach are still strongly problem-dependent, so that some further research efforts are needed for the approach to become a generally applicable alternative to conventional element based techniques. Within the considered field of vibroacoustic analysis, a least-squares T-element approach has been applied recently for acoustic analysis [12] and a hybrid T-element approach for elasto-dynamic analysis [13]. The innovative character of the new prediction technique consists mainly of using a new class of T-complete wave functions, which allow the use of the indirect Trefftz method for coupled vibro-acoustic prediction models, whose poor condition is no longer preventing the numerical results from converging towards the exact solution or does not require anymore to divide the problem domain into small elemental subdomains. This paper presents the basic principles of the new prediction technique and applies it for the modelling of the vibro-acoustic coupling between the pressure field in an acoustic cavity with arbitrary shape and the out-of-plane displacement of a flat plate with arbitrary shape. This application illustrates that a high accuracy is obtained from substantially smaller models, compared to the finite element method, since an approximation is induced only in the representation of the boundary conditions. Moreover, a comparison of the computational efforts, involved with the new technique and the finite element method, indicates the beneficial convergence rate of the new prediction technique. In this way, the new technique can be applied up to much higher frequencies than the existing finite element method.
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
529
governed by the inhomogeneous Kirchhoff equation for plate bending motion 74 wðj; hÞ 2 kb4 wðj; hÞ ¼
pðux ðj; hÞ; uy ðj; hÞ; uz ðj; hÞÞ F dðrsF ðjF ; hF ÞÞ þ ; D D
ð3Þ
where 74 ¼ Fig. 1. Three-dimensional coupled vibro-acoustic problem.
2. Basic concept of the new prediction technique 2.1. Problem definition A detailed discussion of the new prediction technique is given in Ref. [2]. This paper briefly presents its basic principles, when applied for a three-dimensional coupled vibro-acoustic problem, as shown in Fig. 1. The boundary surface Va of an acoustic cavity V consists of four parts. On parts Vp, Vv and VZ, prescribed pressure, normal velocity and normal impedance distributions are imposed, respectively, while part Vs consists of a flat, flexible plate. The cavity domain, whose geometry is defined in an orthogonal co-ordinate system ðx; y; zÞ; is filled with a fluid with an ambient density r0 and speed of sound c. The thin, flat plate in domain Vs has a thickness t and its material has a density rs, an elastic modulus E, a material loss factor h and a Poisson coefficient n. The geometry of the plate middle surface Vs may be described in an orthogonal co-ordinate system ðj; hÞ; which is linked to the co-ordinate system ðx; y; zÞ through a parameterisation of the form ;rs ðj; hÞ ¼ rs ðx; y; zÞ [ Vs 8 x ¼ ux ðj; hÞ ¼ ux;0 þ ux;j j þ ux;h h > > < : y ¼ uy ðj; hÞ ¼ uy;0 þ uy;j j þ uy;h h ; > > : z ¼ uz ðj; hÞ ¼ uz;0 þ uz;j j þ uz;h h
ð1Þ
where ui;0 ; ui;j and ui;h ði ¼ x; y; zÞ are constants. The coupled vibro-acoustic system is dynamically excited by a normal point force excitation F at location rsF on the flexible plate. The dynamic force excitation has a harmonic time dependence with circular frequency v. The steady-state acoustic pressure p at any position rðx; y; zÞ in the cavity domain V is governed by the homogeneous Helmholtz equation 2
2
7 pðrÞ þ k pðrÞ ¼ 0;
ð2Þ
where kð¼ v=cÞ is the acoustic wave number. The steady-state out-of-plane plate displacement w, with positive orientation away from the acoustic cavity, is
›4 ›4 ›4 þ2 2 2 þ ; 4 › h4 ›j ›j ›h
ð4Þ
and where d is a two-dimensional Dirac delta function. The plate bending wave number kb and the plate bending stiffness D are defined as sffiffiffiffiffiffiffiffi 2 4 r tv Et3 ð1 þ jhÞ s kb ¼ ; D¼ ; ð5Þ D 12ð1 2 n 2 Þ pffiffiffiffi with j ¼ 21 the imaginary unit. The acoustic boundary conditions for the coupled vibroacoustic problem are pðrÞ ¼ p ðrÞ;
;r [ Vp ;
ð6Þ
j ›pðrÞ ¼ v n ðrÞ; r0 v ›n
;r [ Vv ;
ð7Þ
j ›pðrÞ pðrÞ ; ¼ ZðrÞ r0 v ›n
;r [ VZ ;
ð8Þ
j ›pðrÞ ¼ jvwðrÞ; r0 v ›n
;r [ Vs ;
ð9Þ
where n represents the normal direction of the boundary surface Va, with positive orientation away from the cavity, and where p ; v n and Z are prescribed pressure, normal velocity and normal impedance functions, respectively. The boundary curve Gs of the plate domain Vs may consist of three parts, as shown in Fig. 2. On the first part Gwu , the out-of-plane translational and the normal rotational displacements are imposed. On the second part Gwm ; the outof-plane translational displacements and the normal bending moments are imposed. On the third part GmQ ; the normal bending moments and the generalised normal shear forces are imposed. By defining the differential operators Lu ¼ 2
› ; ›gn
ð10Þ
! ›2 ›2 Lm ¼ 2D þn 2 ; ›g2n ›gs
› LQ ¼ 2D ›gn
! ›2 ›2 þ ð2 2 nÞ 2 ; ›g2n ›gs
ð11Þ
ð12Þ
where gn and gs are, respectively, the normal and tangential directions of the plate boundary curve (Fig. 2), the structural
530
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
^ approximated as a solution expansion w ^ j ; hÞ wðj; hÞ < wð ¼
ns X
w s C s ð j ; hÞ þ
s¼1
na X
^ F ðj; hÞ ^ a ðj; hÞ þ w pa w
a¼1
^ a {pa } þ w^ F : ¼ ½Cs {ws } þ ½w
ð19Þ
Each function in the ð1 £ ns Þ vector ½Cs is a structural wave function, which satisfies the homogeneous part of the dynamic plate equation (3) Fig. 2. Flexible plate domain Vs :
Cs ðj; hÞ ¼ e2jðkjs ðj2fjs Lj Þþkhs ðh2fhs Lh ÞÞ ;
boundary conditions become ( wðrs Þ ¼ wðr s Þ; ;rs [ Gwu ; Lu ½wðrs Þ ¼ un ðrs Þ; ( wðrs Þ ¼ wðr s Þ; ;rs [ Gwm ; Lm ½wðrs Þ ¼ m n ðrs Þ; ( Lm ½wðrs Þ ¼ m n ðrs Þ; ;rs [ GmQ ; n ðrs Þ; LQ ½wðrs Þ ¼ Q
with ð13Þ
ð14Þ
ð15Þ
n are prescribed functions for, where w; un ; m n and Q respectively, the out-of-plane translational displacement, the normal rotational displacement, the normal bending moment, and the generalised normal shear force. 2.2. Field variable expansions The steady-state cavity pressure p is approximated as a solution expansion p^ pðx; y; zÞ < p^ ðx; y; zÞ ¼
na X
pa Fa ðx; y; zÞ ¼ ½Fa {pa }:
ð16Þ
a¼1
Each function in the ð1 £ na Þ vector [Fa is an acoustic wave function, which satisfies the Helmholtz equation (2)
Fa ðx; y; zÞ ¼ e
2jðkxa ðx2fxa Lx Þþkya ðy2fya Ly Þþkza ðz2fza Lz ÞÞ
;
with 2 2 2 kxa þ kya þ kza ¼ k2 :
ð17Þ
Lx, Ly and Lz are the dimensions of the (smallest) rectangular prism that encloses the cavity domain V. In order to ensure that the amplitudes of the wave functions are not greater than 1 within the cavity domain V, which is beneficial for the numerical condition of the resulting model, the scaling factors fxa, fya and fza are defined as follows ( 1; if Imðka Þ . 0; ð18Þ ða ¼ xa; ya; zaÞ: fa ¼ 0; if Imðka Þ # 0; The contributions of the acoustic wave functions to the solution expansion are comprised in the ðna £ 1Þ vector {pa }: The steady-state out-of-plane plate displacement w is
ðkj2s þ kh2 s Þ2 ¼ kb4 :
ð20Þ
Lj and Lh are the dimensions of the (smallest) rectangle that encloses the plate domain Vs. Again, scaling factors fjs and fhs are introduced to keep the amplitudes of the wave functions smaller than 1 within the plate domain Vs ( 1; if Imðks Þ . 0; fs ¼ ð21Þ ðs ¼ js; hsÞ: 0; if Imðks Þ # 0; The contributions of the structural wave functions to the solution expansion are comprised in the ðns £ 1Þ vector {ws }: ^ F is a particular solution for the mechanical Function w force term in the right-hand side of Eq. (3). Several mathematical expressions may serve as particular solution. It is advantageous, however, to select an expression, which is already close to the physical response of the plate. Therefore, the displacement of an infinite plate, excited by a normal point force F is selected jF ð22Þ ½H ð2Þ ðkb rÞ 2 H0ð2Þ ð2jkb rÞ; 8kb2 D 0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi with r ¼ ðj 2 jF Þ2 þ ðh 2 hF Þ2 and where H0ð2Þ is the zero-order Hankel function of the second kind. ^ a ðj; hÞ in the ð1 £ na Þ vector ½w^ a is a Each function w particular solution for the part of the cavity pressure loading term in the inhomogeneous right-hand side of Eq. (3), resulting from one of the acoustic wave functions Fa in the cavity pressure expansion (16). The displacement of an infinite plate, excited by this distributed pressure, may be used as particular solution function. However, due to the flat shape of the plate, a particular solution function may be defined, which is proportional to the pressure loading function
w^ F ðj; hÞ ¼ 2
w^ a ðj; hÞ ¼
e2jðkxa ðux 2fxa Lx Þþkya ðuy 2fya Ly Þþkza ðuz 2fza Lz ÞÞ ; D½A4j þ 2A2j A2h þ A4h 2 kb4
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
531
with Aj ¼ kxa ux;j þ kya uy;j þ kza uz;j ;
ð23Þ
Ah ¼ kxa ux;h þ kya uy;h þ kza uz;h :
2.3. Coupled vibro-acoustic wave model Due to the particular choice of the field variable expansions (16) and (19), the governing dynamic Eqs. (2) and (3) are exactly satisfied, irrespective of the values of the unknown wave function contributions pa and ws. These contributions are merely determined by the acoustical and structural boundary conditions. To obtain a finite sized wave model, which is amenable to numerical implementation, these boundary conditions are transformed into a weighted residual formulation. 2.3.1. Acoustic boundary conditions Depending on the wave function contributions pa and ws in the field variable expansions (16) and (19), some approximation errors are induced in the representation of the acoustic boundary conditions (6) – (9). For each of these boundary conditions, an associated residual error function may be defined as follows Rp ðrÞ ¼ p^ ðrÞ 2 p ðrÞ;
r [ Vp ; r [ Vv ;
j ›p^ ðrÞ p^ ðrÞ ; 2 RZ ðrÞ ¼ ZðrÞ r0 v ›n
r [ VZ ;
j ›p^ ðrÞ ^ 2 jvwðrÞ; r0 v ›n
2
VZ
r [ Vs :
ð Vp
Vs
2.3.2. Structural boundary conditions In the similar way as for the acoustic boundary conditions, some residual error functions for the structural boundary conditions (13) – (15) may be defined as follows ^ s Þ 2 wðr Rw ðrs Þ ¼ wðr s Þ;
rs [ ðGwu < Gwm Þ;
ð30Þ
r s [ G wu ;
ð31Þ
^ s Þ 2 un ðrs Þ; Ru ðrs Þ ¼ Lu ½wðr ^ s Þ 2 m Rm ðrs Þ ¼ Lm ½wðr n ðrs Þ;
rs [ ðGmQ < Gwm Þ; ð32Þ
n ðrs Þ; ^ s Þ 2 Q RQ ðrs Þ ¼ LQ ½wðr
rs [ GmQ :
ð33Þ
ð25Þ ð26Þ ð27Þ
A weighted residual formulation of the acoustic boundary conditions is defined as ð ð ð p~ Rv dV þ p~ RZ dV þ p~ Rs dV Vv
matrices ½Aaa and ½Caa and the ðna £ ns Þ matrix ½Cas ; the reader is referred to Ref. [2].
ð24Þ
j ›p^ ðrÞ 2 v n ðrÞ; Rv ðrÞ ¼ r0 v ›n
Rs ðrÞ ¼
Fig. 3. Normal and tangential directions at a corner point.
ð28Þ
j ›p~ R dV ¼ 0; r0 v ›n p
where p~ is a weighting function. In a similar way as in the Galerkin weighting procedure, used in the finite element method, each of the acoustic wave functions Fa in the cavity pressure expansion (16) is used as a weighting function p~ in Eq. (28). This yields a set of na algebraic equations in the ðna þ ns Þ unknown wave function contributions pa and ws ( )
pa ¼ {fa }: Aaa þ Caa Cas ð29Þ ws For a full derivation of the coefficients in the ðna £ na Þ
For the sake of model symmetry, some additional residuals are defined at the corner points of the plate boundary curve Gs, at which the normal and tangential directions are not uniquely defined (Fig. 3). For the nw corner points, which belong to the boundary part Gwu < Gwm ; the additional residual is the approximation error for the out-of-plane plate displacement at the corner point location rsc ^ sc Þ 2 wðr Rcw ðrsc Þ ¼ wðr sc Þ;
ðc ¼ 1::nw Þ:
ð34Þ
For the nF corner points, which belong to the boundary part GmQ, the additional residual is the approximation error for the concentrated corner point force, which is associated with the discontinuity of the torsional moment ^ sc Þ 2 ðm RcF ðrsc Þ ¼ LF ½wðr ns ðrþ ns ðr2 sc Þ 2 m sc ÞÞ ^ sc Þ ^ sc Þ ›2 wðr ›2 wðr ¼ 2Dð1 2 nÞ 2 þ þ 2 ›gn ›gs ›gn ›g2 s 2
ðm ns ðrþ sc Þ
2
! ð35Þ
m ns ðr2 sc ÞÞ;
ðc ¼ 1::nF Þ: A weighted residual formulation of the structural boundary
532
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
conditions is defined as ð Gwu
2
LQ ½wR ~ w dG þ
ð Gwm
þ
nw X c¼1
associated uncoupled plate problem without acoustic pressure loading from the cavity.
ð Gwu
Lm ½wR ~ u dG
Lu ½wR ~ m dG 2
LF ½wðr ~ sc ÞRcw ðrsc Þ 2
2.4. Convergence ð GmQ nF X
wR ~ Q dG
wðr ~ sc ÞRcF ðrsc Þ ¼ 0;
ð36Þ
c¼1
where w ~ is a weighting function. Again, by using each of the structural wave functions Cs in the displacement expansion (19) as a weighting function w~ in Eq. (36), a set of ns algebraic equations is obtained in the ðna þ ns Þ unknown wave function contributions pa and ws ( )
pa ¼ {fs }: Csa Ass ð37Þ ws For a full derivation of the coefficients in the ðns £ ns Þ matrix ½Ass and the ðns £ na Þ matrix ½Csa ; the reader is referred to Ref. [2]. 2.3.3. Model properties The combination of the weighted acoustic boundary condition (29) and the weighted structural boundary condition (37) yields the coupled vibro-acoustic model " #( ) ( ) Aaa þ Caa Cas fa pa ¼ : ð38Þ Csa Ass ws fs Solving model (38) for the unknown wave function contributions pa and ws and back-substituting the results in expansions (16) and (19) yield the prediction of the steadystate dynamic response of the coupled vibro-acoustic problem. Several model properties may be identified. † Since the acoustic and structural wave functions Fa and Cs are defined in, respectively, the entire acoustic domain V and the entire structural domain Vs, the model matrices are fully populated. † Since the acoustic and structural wave functions Fa and Cs are complex functions and since they are implicitly dependent on frequency by virtue of the frequency dependence of their wave number components, the model matrices are complex and frequency dependent and cannot be partitioned into frequency independent submatrices. † The model matrix in Eq. (38) is not symmetric. However, it is proved in Ref. [2] that the submatrices ½Aaa and ½Ass are symmetric. Note that these submatrices would become, respectively, the model matrix for the associated uncoupled acoustic problem, having no elastic boundary surface Vs, and the model matrix for the
Since Eqs. (17) and (20) have an infinite number of real and complex solutions, the key question is whether a set of acoustic and structural wave functions can be selected from the infinite number of possible wave functions, such that the subsequent pressure and displacement expansions (16) and (19) converge towards the exact solution. 2.4.1. Convex problem domain The acoustic wave functions Fa with the following wave number components are proposed ffi! m1x p m1y p qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ðkxa ; kya ; kza Þ ¼ ; ; ; ^ k2 2 kxa 2 kya Lx Ly ffi m2x p qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 k2 ; m1z p ; ; ^ k2 2 kxa and za Lx Lz qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m p m p ! 2y 2 2 k2 ; ; ; 2z ^ k2 2 kya za Ly Lz
ð39Þ
with mix ; miy ; miz ¼ 0; ^1; ^2; … ði ¼ 1; 2Þ: The structural wave functions Cs with the following wave number components are proposed ! sx p qffiffiffiffiffiffiffiffiffiffi 2 2 ðkjs ; khs Þ ¼ ; ^ kb 2 kjs ; Lj ! qffiffiffiffiffiffiffiffiffiffi sx p 2 2 ; ^j kb þ kjs ; and Lj qffiffiffiffiffiffiffiffiffiffiffi s p ! qffiffiffiffiffiffiffiffiffiffi s p ! y y 2 2 ; ^j kb2 þ kh2 s ; ; ^ kb 2 khs ; Lh Lh
ð40Þ
with sx, sy ¼ 0; ^1; ^2; … Recall that Lx, Ly and Lz are the dimensions of the (smallest) rectangular prism, enclosing the acoustic cavity domain V, and that Lj and Lh are the dimensions of the (smallest) rectangle, enclosing the plate domain Vs. It is proved in Ref. [2] that the convexity of the acoustic domain V and the plate domain Vs is a sufficient condition for the solution expansions (16) and (19), using the proposed wave function selections (39) and (40), to converge, in the limit for mix, miy, miz, sx, sy : 0 ! 1; towards the exact coupled vibro-acoustic response. Moreover, it can be proved that the proposed acoustic wave
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
533
Fig. 4. Two-dimensional acoustic validation example (convex domain).
functions may be combined into !
m1y p 2jkza ðz2fza Lz Þ m p e Fa ðx; y; zÞ ¼ cos 1x cos ; Lx Ly m2x p 2jkya ðy2fya Ly Þ m1z p cos e ; and cos Lx Lz ! m2y p m2z p 2jkxa ðx2fxa Lx Þ cos cos e ; Lz Ly
ð41Þ
with mix, miy, miz ¼ 0; 1; 2; …; ði ¼ 1; 2Þ: Note that in this ^ a in Eq. (23) must be case, the particular solution functions w modified, as described in Ref. [2]. To illustrate the convergence of the proposed wave function selection, based on the definition of an enclosing rectangular domain, a two-dimensional acoustic validation problem is considered, as shown in Fig. 4. At the right boundary of a slender convex acoustic cavity (length of 5 m and height of 1 m), a real normal impedance of 2000 kg/(m2 s) is prescribed. The upper and lower parts of the boundary surface are assumed to be rigid. The dynamic excitation of the air-filled cavity consists of a unit prescribed normal velocity, imposed at the left boundary. The air in the cavity has a density r0 ¼ 1.225 kg/m3 and a speed of sound c ¼ 340 m/s. For this two-dimensional example, the two-dimensional counterpart of the acoustic wave function selection, proposed in Eq. (39), should be used, i.e. qffiffiffiffiffiffiffiffiffiffi mx p 2 ; ðkxa ; kya Þ ¼ and ; ^ k2 2 kxa Lx ð42Þ qffiffiffiffiffiffiffiffiffiffi m p ! y 2 2 ; ^ k 2 kya ; Ly with Lx and Ly the dimensions of an enclosing rectangle ðmx ; my ¼ 0; ^1; ^2; …Þ: Fig. 5 shows the acoustic frequency response function (FRF), relating the pressure at the response location (indicated in Fig. 4) to the unit velocity input, obtained from a wave model with a total of 120 acoustic wave functions. The enclosing rectangle, used for the wave function selection, is indicated in dashed lines in Fig. 4. To compare the performance of the new prediction technique with the finite element method, three FE models are
constructed for the considered validation example. In these models, the two-dimensional acoustic cavity is discretised into 4-noded linear elements, resulting in FE model sizes of, respectively, 105, 369 and 1292 nodal pressure degrees of freedom. The corresponding acoustic FRFs, obtained from these three FE models, are also plotted (in dashed lines) in Fig. 5. This figure clearly illustrates that, by increasing the number of elements, the prediction accuracy increases and the finite element results gradually converge towards the results, obtained from the small wave model. This figure illustrates also the typical approximation errors, involved with FE discretisations for dynamic analysis. The use of simple, but approximating (polynomial)
Fig. 5. Acoustic FRF at response point (solid: wave model—120 dof, dashed: FE models (upper: 105 dof, middle: 369 dof, lower: 1292 dof)).
534
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
Fig. 6. Acoustic convergence curve at 100 Hz (accuracy versus dof) (wave models (solid)—FE models (dashed)).
shape functions, describing the dynamic response within each element, involves not only an interpolation error, but also a numerical dispersion error, in that the wave number of the numerical solution differs from the wave number of the exact problem solution [14]. This error results usually in an overestimation of the resonance frequencies. Since the dispersion error increases for increasing frequency, the mesh density and subsequent computational effort, required for getting reasonable prediction accuracy, increase also with frequency. Another indication of the convergence of the proposed technique is given by the convergence curve, shown in Fig. 6. This curve plots the relative prediction error on the pressure at the considered response point against the number of degrees of freedom. The relative error is defined as the difference between the calculated and
Fig. 7. Two-dimensional vibro-acoustic validation example.
exact solution, divided by the exact solution. However, since the exact solution for the considered problem is not known analytically, the result, obtained from a wave model with a very large number of wave functions, is used instead to define the relative prediction error. The convergence curves, associated with both the novel wave based prediction technique (solid line) and with the finite element method (dashed line), are plotted in Fig. 6 for the case of a dynamic excitation at 100 Hz. This figure illustrates not only that the novel technique provides a higher accuracy with a smaller number of degrees of freedom, but also that its convergence rate, i.e. the slope of the convergence curve, is substantially higher. The size of the enclosing rectangle has some influence on the overall accuracy, as is illustrated for the twodimensional vibro-acoustic validation problem, shown in Fig. 7. One part of the boundary surface of an air-filled (r0 ¼ 1.225 kg/m3, c ¼ 340 m/s) cavity (length of 1.5 m and height of 1 m) consists of a 0.75 m long, flexible plate (E ¼ 70 £ 109 N/m2 , rs ¼ 2790 kg/m3 , n ¼ 0.3, t ¼ 0.002 m) with clamped boundaries, while the remaining part of the cavity boundary surface is perfectly rigid. The coupled vibro-acoustic system is excited by a unit normal line force F at j F ¼ 0.25 m with a harmonic time dependence at 200 Hz. As shown in Fig. 8, two different enclosing rectangles have been used for the definition of the wave number components. This figure shows the resulting convergence curves for the out-of-plane plate displacement at the excitation point. It illustrates that it is beneficial for both the accuracy and the convergence rate to select the wave functions, associated with the smallest rectangular domain, enclosing the considered problem domain. 2.4.2. Substructuring So far, the convergence of the novel prediction technique is discussed for convex problem domains. Although a firm mathematical proof is still under construction, the authors have not yet found an example, counterproving the statement that any type of non-convex cavity or plate domain can be decomposed into some subdomains and that a convergent pressure or displacement expansion for each subsequent subdomain is obtained by using a wave function selection of type (39) or (40), which is based on the dimensions of a rectangular domain, enclosing the subdomain. When the acoustic domain V must be decomposed for convergence reasons, a weighted residual formulation of type (28) must be specified for each acoustic subdomain, resulting in a wave model of type (29) for each subdomain. In this respect, some additional boundary conditions must be specified to ensure the continuity of the pressure and normal fluid velocity along the common boundary interfaces between the various subdomains. Hence, for each pair of adjacent acoustic subdomains, a pressure residual function of type (24), specifying the numerical pressure difference along the common interface, must be added to the
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
535
Fig. 10. Substructuring of the concave domain.
Fig. 8. Influence of the size of the enclosing rectangular domain on the convergence.
weighted residual formulation, associated with one subdomain, and a normal velocity residual function of type (25), specifying the numerical normal velocity difference along the common interface, must be added to the weighted residual formulation, associated with the other subdomain. In the same way, when the plate domain Vs must be
Fig. 9. Two-dimensional acoustic validation example (concave domain).
decomposed, a weighted residual formulation of type (36) must be specified for each structural subdomain, resulting in a wave model of type (37) for each subdomain. For each pair of adjacent structural subdomains, a residual function of type (30), specifying the numerical out-of-plane displacement difference along the common interface, and a residual function of type (31), specifying the numerical difference in normal rotational displacement along the interface, must be added to the weighted residual formulation, associated with one subdomain. In addition, residual functions of type (32) and (33) must be added to the weighted residual formulation, associated with the other subdomain, to ensure the balances of the generalised normal shear forces and the normal bending moments along the common interface. The total coupled vibro-acoustic model for a non-convex problem is then obtained from combining the weighted residual formulations of all acoustic and all structural subdomains. The motivation for the proposed substructuring procedure for concave problems is similar to the one in the T-element approach, discussed in Section 1, in that convergence is quested by dividing the problem domain. However, due to the advantageous properties of the proposed wave function selections (39) and (40), the problem domain must only be divided in a few, large subdomains, rather than in many small elemental subdomains as is the case in the T-element approach. To illustrate the convergence of the proposed substructuring procedure, a two-dimensional acoustic validation problem is considered, as shown in Fig. 9. At the right boundary of an L-shaped acoustic cavity (outer length of 4 m and height of 2 m), a real normal impedance of 2000 kg/(m2 s) is prescribed. The dynamic excitation of the airfilled cavity consists of a unit prescribed normal velocity, imposed at the upper part of the boundary surface. The remaining boundary parts are assumed to be rigid. The air in the cavity has a density r0 ¼ 1.225 kg/m3 and a speed of sound c ¼ 340 m/s. For this concave problem, the L-shaped cavity is divided in two subdomains, as shown in Fig. 10. For both subdomains, a wave function selection of type (42) is used. Fig. 11 shows the acoustic FRF, relating the pressure at the response location (indicated in Fig. 9) to the unit velocity input, obtained from a wave model with a total of 108 wave functions. To compare the performance of the prediction technique with the finite element method, three
536
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
acoustic wavelengths at the considered excitation frequency. When the excitation frequency v coincides with an undamped structural eigenfrequency of the circumscribing rectangular plate domain with simply supported boundaries, i.e. !sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 q22 Et2 2 q1 v¼p þ 2 ; 2 Lj Lh 12ð1 2 n2 Þrs ð44Þ ðq1 ; q2 ¼ 1; 2; …Þ; or when it coincides with an undamped acoustic eigenfrequency of the circumscribing rectangular cavity with rigid side walls, i.e. vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !2 u u a1 p 2 a2 p a3 p 2 t ; v¼c þ þ Lx Ly Lz ð45Þ ða1 ; a2 ; a3 ¼ 0; 1; …Þ; it follows from Eqs. (39) and (40) that not all wave functions are linearly independent and that some wave functions must be eliminated to restore the linear independence without reducing the span of the wave function sets.
3. Comparison with the finite element method Fig. 11. Acoustic FRF at response point (solid: wave model—108 dof, dashed: FE models (upper: 105 dof, middle: 217 dof, lower: 793 dof)).
FE models are constructed for this validation problem using 4-noded linear elements, resulting in model sizes of, respectively, 105, 217 and 793 nodal pressure degrees of freedom. The corresponding acoustic FRFs, obtained from these three FE models, are also plotted (in dashed lines) in Fig. 11. Again, this figure illustrates that, by increasing the number of elements, the prediction accuracy increases and the finite element results gradually converge towards the results, obtained from the small wave model. 2.4.3. Truncation rule For the numerical implementation of the coupled wave model (38), the infinite wave function series (mix, miy, miz, sx, sy : 0 ! 1) must be truncated at some finite integer values m ix ; m iy ; m iz ; s x ; sy : In this respect, the following truncation rule is proposed m sy iy m m s 2maxðk; kb Þ ix : < < iz < x < $ p Lx Ly Lz Lj Lh
ð43Þ
This truncation rule expresses that the smallest wavelength components in the structural and acoustic wave function sets may not be larger than half the structural and half the
In comparison with the finite element method, the new prediction technique has several advantageous properties. † Since the field variable expansions satisfy a priori the governing dynamic equations, approximation errors are only involved with the representation of the boundary conditions. This yields substantially smaller models than corresponding FE models. † Predictions of dynamic quantities such as structural stresses and strains or fluid velocities are obtained from deriving the prediction results of the primary field variables, i.e. the structural displacement components or the acoustic pressure. In the FE method, the latter are usually expressed in terms of polynomial expansions, which are complete up to a certain order. Hence, the derived secondary variables are also expressed in terms of polynomial expansions, which are, however, only complete up to a lower-order. As a result, the lower-order (polynomial) expansions of the derived secondary variables represent a smaller spatial variation and are less accurate than the expansions of the corresponding primary variables. In the proposed prediction technique, however, the spatial gradients of the wave functions are wave functions, which have the same spatial variation as the ones, used in the primary field variable expansions. This is advantageous for the convergence rate, especially in the case of coupled vibro-acoustic problems, for which
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
the effect of the fluid on the structure is pressure controlled, but for which the effect of the structure on the fluid is velocity controlled. † A model refinement is easily obtained. Firstly, the accuracy of a model is easily checked, since only the approximation errors of the boundary conditions must be verified. Secondly, increasing the number of degrees of freedom, i.e. the number of wave functions in the field variable expansions, requires only the calculation of the additional matrix coefficients, associated with the additional wave functions, while the original model coefficients remain unchanged. For FE models, increasing the number of degrees of freedom requires the calculation of a completely new model, at least when a global refinement of the element discretisation is needed. The new prediction technique has also some drawbacks, compared with the FE method. † Since the wave functions are complex functions, defined within the entire continuum domain or at least within large subdomains, the wave model matrices are fully populated and complex. † Due to the implicit frequency dependence of the wave functions, the wave model matrices cannot be decomposed into frequency independent submatrices, which implies that a wave model must be recalculated for each frequency of interest and that natural frequencies and mode shapes cannot be obtained from a standard eigenvalue calculation. † Due to the global and oscillatory nature of the wave functions and due to the poor condition of a wave model, the involved numerical integrations are computationally more demanding. † The FE method is a widely applicable technique, since the convergence of a FE model is obtained with simple element shape functions, whereas the convergence of a wave model requires the definition of a complete wave function set and appropriate particular solution functions.
537
Fig. 12. Three-dimensional vibro-acoustic validation example.
The acoustic cavity is comprised in the volume of an enclosing rectangular prism with dimensions Lx ¼ 1:5 m; Ly ¼ 0:5 m; Lz ¼ 1 m: The vibro-acoustic system is excited by a unit mechanical point force F, applied at location ðjF ; hF Þ ¼ ð2Lj =3; Lh =3Þ on the plate, in the direction normal to the plate and with a harmonic time dependence at a certain frequency v. Fig. 13 shows the real part of the out-of-plane plate displacement at 200 Hz, obtained from a wave model with 314 (na ¼ 258 acoustic and ns ¼ 56 structural) wave functions, and the corresponding normal fluid displacement along the cavity – plate coupling interface. This figure indicates that the wave model enables an accurate representation of the clamped plate boundary conditions of type (13) and the normal displacement continuity condition of type (9). Fig. 14 shows the contour plot of
4. Three-dimensional validation example To illustrate the accuracy of the proposed prediction technique for three-dimensional coupled vibro-acoustic analysis and to assess its computational performance in comparison with the finite element method, the numerical validation example is considered, as shown in Fig. 12. One boundary surface of an acoustic cavity consists of a flexible rectangular ðLj ¼ 0:75 m; Lh ¼ 0:5 mÞ plate with clamped boundaries, while all other cavity boundary surfaces are perfectly rigid. The air in the cavity has a density r0 ¼ 1:225 kg=m3 and a speed of sound c ¼ 340 m=s: The undamped aluminium plate has a thickness t ¼ 2 £ 1023 m; a density rs ¼ 2790 kg=m3 ; a Poisson coefficient n ¼ 0:3 and an elasticity modulus E ¼ 70 £ 109 N=m2 :
Fig. 13. Real part of the out-of-plane plate displacement at 200 Hz (a) and the corresponding normal fluid displacement along the cavity – plate coupling interface (b).
538
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
Table 1 Properties of the FE models
FE model 1 FE model 2 FE model 3
No. of structural elements
No. of acoustic elements
dof
54 216 486
1296 10 368 34 992
1849 12 586 39 991
the real part of the acoustic pressure in the cavity plane y ¼ Ly =3: Since the pressure contour lines are perpendicular to the rigid cavity boundary surfaces, the zero normal fluid velocity boundary conditions of type (7) at these boundary surfaces are accurately represented. Since the solution expansions (16) and (19) inherently satisfy the governing dynamic equations, and since all boundary conditions of the considered problem are accurately represented, this validation example illustrates that the proposed prediction technique is practically converging and enables an accurate approximation of the exact three-dimensional vibro-acoustic response. To compare the performance of the new prediction technique with the finite element method, three FE models are constructed for the considered validation example. In these models, the plate is discretised into four-noded quadrilateral shell elements and the acoustic cavity is discretised into eight-noded hexahedral fluid elements. Table 1 lists the number of elements in each model and their corresponding number of unconstrained degrees-offreedom (dof). Fig. 15 compares the direct frequency response functions (with a frequency resolution of 1 Hz), obtained from the aforementioned wave model with 314 wave functions and from the three FE models. Fig. 16 compares the associated FRFs of the pressure at location ðx; y; zÞ ¼ ð2Lx =3; 2Ly =3; Lz =4Þ: Again, these figures clearly shown that, by increasing the number of elements, the prediction accuracy increases and the finite element results gradually
converge towards the results, obtained from the small wave model. In order to assess the convergence rate, several wave models have been constructed, as listed in Table 2. The solid lines on Fig. 17 indicate the resulting relative prediction errors for the out-of-plane plate displacement at the excitation point ðjF ; hF Þ ¼ ð2Lj =3; Lh =3Þ and for the cavity pressure at ðx; y; zÞ ¼ ð2Lx =3; 2Ly =3; Lz =4Þ for two excitation frequencies, i.e. 60 and 180 Hz. The dashed lines indicate the corresponding convergence curves for the finite element models, listed in Table 1. This figure shows that the new prediction technique yields accurate prediction results with substantially smaller models than the finite element method and with a higher convergence rate.
Fig. 14. Contour plot of the real part of the acoustic pressure in the cavity plane y ¼ Ly =3 at 200 Hz.
Fig. 15. Direct FRF (solid: wave model, dashed: FE model 1 (a), FE model 2 (b) and FE model 3 (c)).
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
539
Table 2 Properties of the wave models
Wave model 1 Wave model 2 Wave model 3
No. of structural wave functions
No. of acoustic wave functions
dof
32 56 80
78 258 542
110 314 622
However, as indicated in Section 3, the beneficial effect of a smaller model size on the computational effort is partly annihilated by the fact that wave models are fully populated and cannot be decomposed into frequency-independent submatrices. Therefore, it is more appropriate to consider the prediction accuracy in terms of required CPU time rather than model size, as shown in Fig. 18. The CPU times in this figure denote the times, needed for a direct response calculation at one frequency. The indicated CPU times for the wave models include both the times for constructing the
Fig. 17. Structural (left) and acoustic (right) convergence curves (accuracy versus dof) (wave models (solid)—FE models (dashed)) (þ : 60 Hz, L: 180 Hz).
model as well as for solving the resulting matrix equation. Since FE models consist of frequency-independent submatrices, these submatrices must only be calculated once and the construction of the models at each frequency requires only a simple submatrix assembly. As a result, the computational effort for constructing the models is negligible, compared with the efforts for solving the (large) models, at least when the dynamic response is calculated in a large frequency range. Therefore, the indicated CPU times for the FE models comprise only the direct solution times. All calculations are performed on a HP-C180 workstation (SPECfp95 ¼ 18.7, SPECint95 ¼ 11.8). Fig. 18 clearly shows the beneficial convergence rate of the new prediction technique. In comparison with the finite element method, the technique provides highly accurate predictions of the coupled vibro-acoustic response with a
Fig. 16. Acoustic FRF at cavity location ðx; y; zÞ ¼ ð2Lx =3; 2Ly =3; Lz =4Þ (solid: wave model, dashed: FE model 1 (a), FE model 2 (b) and FE model 3 (c)).
Fig. 18. Structural (left) and acoustic (right) convergence curves (accuracy versus CPU time) (wave models (solid)—FE models (dashed)) (þ: 60 Hz, L: 180 Hz).
540
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
substantially smaller computational effort. Note also that the high computational efficiency will most likely become even more apparent, when the technique is implemented in an efficient programming environment instead of the currently used MATLAB environment.
5. Conclusion A new prediction technique for coupled vibro-acoustic analysis has been developed. The technique is based on the indirect Trefftz method and approximates the structural and acoustic field variables of a coupled vibro-acoustic system in terms of solution expansions, which exactly satisfy the governing dynamic equations. The contribution factors in these expansions are obtained from a weighted residual formulation of the boundary conditions. The paper describes the basic concept of the new technique for the modelling of the vibro-acoustic coupling between the pressure field in an acoustic cavity with arbitrary shape and the out-of-plane displacement of a flat plate with arbitrary shape. It is illustrated through a numerical validation example that, in comparison with the finite element method, the new prediction method provides accurate prediction results with a substantially smaller computational effort. Due to this beneficial convergence rate, the proposed prediction technique may be a suitable method for obtaining accurate coupled vibro-acoustic predictions up to much higher frequencies than the conventional finite element method.
Acknowledgements Wim Desmet is a Postdoctoral Fellow of the Fund for Scientific Research, Flanders (Belgium) (FWO). The research work of Bas van Hal is financed by a scholarship
of the Institute for the Promotion of Innovation by Science and Technology in Flanders (IWT).
References [1] Craggs A. An acoustic finite element approach for studying boundary flexibility and sound transmission between irregular enclosures. J Sound Vibr 1973;30:343–57. [2] Desmet W. A wave based prediction technique for coupled vibroacoustic analysis. PhD Dissertation. K.U. Leuven, Leuven; 1998. [3] Trefftz E. Ein Gegenstu¨ck zum Ritzschen Verfahren. Proc 2nd Int Cong Appl Mech, Zurich 1926;131–7. [4] Kita E, Kamiya N. Trefftz method: an overview. Adv Engng Software 1995;24:3–12. [5] Herrera I. Boundary methods: an algebraic theory. Boston: Pitman Advanced Publishing Program; 1984. [6] Masson P, Redon E, Priou JP, Gervais Y. The application of the Trefftz method to acoustics. In: Crocker M, editor. Proceedings of the Third International Congress on Air- and Structure-borne Sound and Vibration, Montreal; 1994. p. 1809–16. [7] Langley RS. Some perspectives on wave-mode duality in SEA. In: Fahy FJ, Price WG, editors. Solid mechanics and its applications. IUTAM Symposium on Statistical Energy Analysis, Dordrecht: Kluwer; 1997. p. 1– 12. [8] Jirousek J, Wroblewski A. T-elements: state-of-the-art and future trends. Arch Comput Meth Engng 1998;41(5):831– 49. [9] Zielinski AP, Zienkiewicz OC. Generalized finite element analysis with T-complete boundary solution functions. Int J Numer Meth Engng 1985;21:509 –28. [10] Jirousek J, Wroblewski A. Least-squares T-elements: equivalent FE and BE forms of a substructure-oriented boundary solution approach. Commun Numer Meth Engng 1994;10:21–32. [11] Piltner R. Special finite elements with holes and internal cracks. Int J Numer Meth Engng 1985;21:1471–85. [12] Stojek M. Least-squares Trefftz-type elements for the Helmholtz equation. Int J Numer Meth Engng 1998;41:831–49. [13] Texeira de Freitas JA. Hybrid finite element formulations for elastodynamic analysis in the frequency domain. Int J Solids Struct 1998;36:1883–923. [14] 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 1999;46:471 –99.
A computationally efficient prediction technique for the steady-state dynamic analysis of coupled vibro-acoustic systems W. Desmet*, B. van Hal, P. Sas, D. Vandepitte K.U. Leuven, Department of Mechanical Engineering, Division PMA, Celestijnenlaan 300 B, B-3001 Leuven (Heverlee), Belgium Received 14 November 2000; accepted 1 July 2002
Abstract A new prediction technique has been developed for the steady-state dynamic analysis of coupled vibro-acoustic systems. In contrast with the finite element method, in which the dynamic field variables within each element are expanded in terms of local, non-exact shape functions, the dynamic field variables are expressed as global wave function expansions, which exactly satisfy the governing dynamic equations. The contributions of the wave functions to the coupled vibro-acoustic response result from an integral formulation of the boundary conditions. This paper describes the basic concept of the new technique for the modelling of the vibro-acoustic coupling between the pressure field in an acoustic cavity with arbitrary shape and the out-of-plane displacement of a flat plate with arbitrary shape. It is illustrated through a threedimensional validation example that the new prediction technique yields a high accuracy with a substantially smaller computational effort than the finite element method, so that the new prediction technique can be applied up to much higher frequencies. q 2002 Civil-Comp Ltd and Elsevier Science Ltd. All rights reserved. Keywords: Numerical prediction technique; Coupled vibro-acoustic systems; Trefftz-method; Computational efficiency; Steady-state dynamic behaviour; Mid-frequency technique
1. Introduction At present, the finite element and boundary element methods are the most commonly used numerical prediction techniques for solving steady-state dynamic problems, defined in structural and acoustic continuum domains. An important implication of the conventional element concept, using non-exact shape functions, is that the large model sizes involved practically restrict the applicability of these deterministic prediction techniques to the low-frequency range. Above a certain frequency limit, these methods would require, even with the nowadays-available highperformance computer resources, a prohibitively large amount of computational effort and memory resources to get an acceptable level of prediction accuracy. For coupled vibro-acoustic problems, the conflict between computational effort and accuracy is even more pronounced than for uncoupled structural or uncoupled acoustic problems. First of all, the size of coupled prediction * Corresponding author. Tel.: þ 32-16-322480; fax: þ 32-16-322987. E-mail address: [email protected] (W. Desmet).
models is substantially larger, since a structural and an acoustic problem must be solved simultaneously to incorporate the acoustic pressure loading on the elastic structure and the continuity of the normal fluid and structural displacements along the fluid –structure coupling interface. Secondly, the numerical solution procedure for coupled vibro-acoustic models is less efficient than for uncoupled structural or uncoupled acoustic models, since coupled models, at least in their most commonly used acoustic pressure/structural displacement formulation [1], are no longer symmetric. Finally, the efficiency of the modal expansion method in reducing the model size is significantly reduced for coupled vibro-acoustic problems. The most appropriate mode selection for a modal expansion is the use of the modes of the coupled system. However, since coupled models are no longer symmetric, the calculation of coupled modes requires a non-symmetric eigenvalue calculation, which is very time consuming and which makes it a practically impossible procedure for many real-life vibroacoustic problems. The most commonly used alternative is a modal expansion in terms of uncoupled structural and uncoupled acoustic modes, which result from symmetric
0965-9978/02/$ - see front matter q 2002 Civil-Comp Ltd and Elsevier Science Ltd. All rights reserved. PII: S 0 9 6 5 - 9 9 7 8 ( 0 2 ) 0 0 0 6 2 - 5
528
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
and computationally efficient eigenvalue problems. However, the fact that uncoupled acoustic modes have a zero displacement component, normal to the fluid – structure coupling interface, implies that a large number of high-order uncoupled acoustic modes is required to accurately represent the normal displacement continuity along the fluid –structure interface. Hence, the benefit of a computationally efficient construction of the modal base is significantly reduced by the smaller model size reduction, obtained with an uncoupled modal base. As a result, there exists still a substantial discrepancy between the limited frequency range, in which conventional element based models can provide accurate coupled vibroacoustic predictions, and the significantly larger frequency range, in which accurate, deterministic predictions are needed for many industrial engineering problems. In view of narrowing this currently existing frequency discrepancy, a new deterministic prediction technique has been developed [2], which is based on the indirect Trefftz method [3,4]. The steady-state dynamic field variables in the entire—or at least in large subdomains of the—acoustic and structural domains of a coupled vibro-acoustic system are approximated in terms of a set of acoustic and structural wave functions, which are exact solutions of the homogeneous parts of the governing dynamic equations, and some particular solution functions of the inhomogeneous equations. In this way, the governing dynamic equations are exactly satisfied, irrespective of the contributions of the wave functions to the field variable expansions. Since the proposed field variable expansions are exact solutions of the governing dynamic equations, the contributions of the wave functions to the field variable expansions are merely determined by the acoustic and structural boundary conditions. Since only finite sized prediction models are amenable to numerical implementation, the boundary conditions can only be satisfied in an approximate way. Therefore, the wave function contributions result from a weighted residual formulation of the boundary conditions. A key issue in the use of the indirect Trefftz method is the definition of a T-complete function set, which ensures the convergence of the subsequent field variable expansions towards the exact solutions. Several T-complete function sets have already been defined for solving steady-state acoustic problems [5,6] and steady-state dynamic plate bending problems [7]. Although the theoretical convergence for these function sets has been proven, their practical convergence is, however, seriously disturbed or even prevented, due to the ill-conditioning of the involved model matrices. These numerical condition problems may be circumvented, to some extent, by subdividing the considered continuum domain into small elemental subdomains. This has led to the development of the T(refftz)-element approach [8], which allows the introduction of the Trefftz-
idea into a standard finite element scheme: the internal field variables within the T-elements are approximated in terms of a suitably truncated non-conforming T-complete set of trial functions, satisfying the governing equations a priori, while the boundary conditions and interelement continuity are enforced in an average integral sense. Two main approaches may be distinguished. In the hybrid T-element approach, an auxiliary frame is defined at the element boundaries. A conforming frame field distribution in terms of conventional (polynomial) trial functions is defined on the boundary frame and is used to enforce the boundary conditions and interelement conformity. The least-squares T-element approach is based only on the contributions of the T-complete functions to the field variable distributions within the element, without the introduction of auxiliary frame degrees of freedom. Zielinski and Zienkiewicz [9] indicated how the interelement conformity and the boundary conditions may be enforced directly in a least-squares sense. Jirousek and Wroblewski [10] translated this approach into an element formulation, which is suitable for implementation into a standard finite element code. One of the advantages of the T-element approach is its ability to include special-purpose functions to the internal element solution expansions, so that singularity problems can be handled efficiently without troublesome local mesh refinement [11]. However, in its current stage of development, the efficiency and accuracy of the T-element approach are still strongly problem-dependent, so that some further research efforts are needed for the approach to become a generally applicable alternative to conventional element based techniques. Within the considered field of vibroacoustic analysis, a least-squares T-element approach has been applied recently for acoustic analysis [12] and a hybrid T-element approach for elasto-dynamic analysis [13]. The innovative character of the new prediction technique consists mainly of using a new class of T-complete wave functions, which allow the use of the indirect Trefftz method for coupled vibro-acoustic prediction models, whose poor condition is no longer preventing the numerical results from converging towards the exact solution or does not require anymore to divide the problem domain into small elemental subdomains. This paper presents the basic principles of the new prediction technique and applies it for the modelling of the vibro-acoustic coupling between the pressure field in an acoustic cavity with arbitrary shape and the out-of-plane displacement of a flat plate with arbitrary shape. This application illustrates that a high accuracy is obtained from substantially smaller models, compared to the finite element method, since an approximation is induced only in the representation of the boundary conditions. Moreover, a comparison of the computational efforts, involved with the new technique and the finite element method, indicates the beneficial convergence rate of the new prediction technique. In this way, the new technique can be applied up to much higher frequencies than the existing finite element method.
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
529
governed by the inhomogeneous Kirchhoff equation for plate bending motion 74 wðj; hÞ 2 kb4 wðj; hÞ ¼
pðux ðj; hÞ; uy ðj; hÞ; uz ðj; hÞÞ F dðrsF ðjF ; hF ÞÞ þ ; D D
ð3Þ
where 74 ¼ Fig. 1. Three-dimensional coupled vibro-acoustic problem.
2. Basic concept of the new prediction technique 2.1. Problem definition A detailed discussion of the new prediction technique is given in Ref. [2]. This paper briefly presents its basic principles, when applied for a three-dimensional coupled vibro-acoustic problem, as shown in Fig. 1. The boundary surface Va of an acoustic cavity V consists of four parts. On parts Vp, Vv and VZ, prescribed pressure, normal velocity and normal impedance distributions are imposed, respectively, while part Vs consists of a flat, flexible plate. The cavity domain, whose geometry is defined in an orthogonal co-ordinate system ðx; y; zÞ; is filled with a fluid with an ambient density r0 and speed of sound c. The thin, flat plate in domain Vs has a thickness t and its material has a density rs, an elastic modulus E, a material loss factor h and a Poisson coefficient n. The geometry of the plate middle surface Vs may be described in an orthogonal co-ordinate system ðj; hÞ; which is linked to the co-ordinate system ðx; y; zÞ through a parameterisation of the form ;rs ðj; hÞ ¼ rs ðx; y; zÞ [ Vs 8 x ¼ ux ðj; hÞ ¼ ux;0 þ ux;j j þ ux;h h > > < : y ¼ uy ðj; hÞ ¼ uy;0 þ uy;j j þ uy;h h ; > > : z ¼ uz ðj; hÞ ¼ uz;0 þ uz;j j þ uz;h h
ð1Þ
where ui;0 ; ui;j and ui;h ði ¼ x; y; zÞ are constants. The coupled vibro-acoustic system is dynamically excited by a normal point force excitation F at location rsF on the flexible plate. The dynamic force excitation has a harmonic time dependence with circular frequency v. The steady-state acoustic pressure p at any position rðx; y; zÞ in the cavity domain V is governed by the homogeneous Helmholtz equation 2
2
7 pðrÞ þ k pðrÞ ¼ 0;
ð2Þ
where kð¼ v=cÞ is the acoustic wave number. The steady-state out-of-plane plate displacement w, with positive orientation away from the acoustic cavity, is
›4 ›4 ›4 þ2 2 2 þ ; 4 › h4 ›j ›j ›h
ð4Þ
and where d is a two-dimensional Dirac delta function. The plate bending wave number kb and the plate bending stiffness D are defined as sffiffiffiffiffiffiffiffi 2 4 r tv Et3 ð1 þ jhÞ s kb ¼ ; D¼ ; ð5Þ D 12ð1 2 n 2 Þ pffiffiffiffi with j ¼ 21 the imaginary unit. The acoustic boundary conditions for the coupled vibroacoustic problem are pðrÞ ¼ p ðrÞ;
;r [ Vp ;
ð6Þ
j ›pðrÞ ¼ v n ðrÞ; r0 v ›n
;r [ Vv ;
ð7Þ
j ›pðrÞ pðrÞ ; ¼ ZðrÞ r0 v ›n
;r [ VZ ;
ð8Þ
j ›pðrÞ ¼ jvwðrÞ; r0 v ›n
;r [ Vs ;
ð9Þ
where n represents the normal direction of the boundary surface Va, with positive orientation away from the cavity, and where p ; v n and Z are prescribed pressure, normal velocity and normal impedance functions, respectively. The boundary curve Gs of the plate domain Vs may consist of three parts, as shown in Fig. 2. On the first part Gwu , the out-of-plane translational and the normal rotational displacements are imposed. On the second part Gwm ; the outof-plane translational displacements and the normal bending moments are imposed. On the third part GmQ ; the normal bending moments and the generalised normal shear forces are imposed. By defining the differential operators Lu ¼ 2
› ; ›gn
ð10Þ
! ›2 ›2 Lm ¼ 2D þn 2 ; ›g2n ›gs
› LQ ¼ 2D ›gn
! ›2 ›2 þ ð2 2 nÞ 2 ; ›g2n ›gs
ð11Þ
ð12Þ
where gn and gs are, respectively, the normal and tangential directions of the plate boundary curve (Fig. 2), the structural
530
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
^ approximated as a solution expansion w ^ j ; hÞ wðj; hÞ < wð ¼
ns X
w s C s ð j ; hÞ þ
s¼1
na X
^ F ðj; hÞ ^ a ðj; hÞ þ w pa w
a¼1
^ a {pa } þ w^ F : ¼ ½Cs {ws } þ ½w
ð19Þ
Each function in the ð1 £ ns Þ vector ½Cs is a structural wave function, which satisfies the homogeneous part of the dynamic plate equation (3) Fig. 2. Flexible plate domain Vs :
Cs ðj; hÞ ¼ e2jðkjs ðj2fjs Lj Þþkhs ðh2fhs Lh ÞÞ ;
boundary conditions become ( wðrs Þ ¼ wðr s Þ; ;rs [ Gwu ; Lu ½wðrs Þ ¼ un ðrs Þ; ( wðrs Þ ¼ wðr s Þ; ;rs [ Gwm ; Lm ½wðrs Þ ¼ m n ðrs Þ; ( Lm ½wðrs Þ ¼ m n ðrs Þ; ;rs [ GmQ ; n ðrs Þ; LQ ½wðrs Þ ¼ Q
with ð13Þ
ð14Þ
ð15Þ
n are prescribed functions for, where w; un ; m n and Q respectively, the out-of-plane translational displacement, the normal rotational displacement, the normal bending moment, and the generalised normal shear force. 2.2. Field variable expansions The steady-state cavity pressure p is approximated as a solution expansion p^ pðx; y; zÞ < p^ ðx; y; zÞ ¼
na X
pa Fa ðx; y; zÞ ¼ ½Fa {pa }:
ð16Þ
a¼1
Each function in the ð1 £ na Þ vector [Fa is an acoustic wave function, which satisfies the Helmholtz equation (2)
Fa ðx; y; zÞ ¼ e
2jðkxa ðx2fxa Lx Þþkya ðy2fya Ly Þþkza ðz2fza Lz ÞÞ
;
with 2 2 2 kxa þ kya þ kza ¼ k2 :
ð17Þ
Lx, Ly and Lz are the dimensions of the (smallest) rectangular prism that encloses the cavity domain V. In order to ensure that the amplitudes of the wave functions are not greater than 1 within the cavity domain V, which is beneficial for the numerical condition of the resulting model, the scaling factors fxa, fya and fza are defined as follows ( 1; if Imðka Þ . 0; ð18Þ ða ¼ xa; ya; zaÞ: fa ¼ 0; if Imðka Þ # 0; The contributions of the acoustic wave functions to the solution expansion are comprised in the ðna £ 1Þ vector {pa }: The steady-state out-of-plane plate displacement w is
ðkj2s þ kh2 s Þ2 ¼ kb4 :
ð20Þ
Lj and Lh are the dimensions of the (smallest) rectangle that encloses the plate domain Vs. Again, scaling factors fjs and fhs are introduced to keep the amplitudes of the wave functions smaller than 1 within the plate domain Vs ( 1; if Imðks Þ . 0; fs ¼ ð21Þ ðs ¼ js; hsÞ: 0; if Imðks Þ # 0; The contributions of the structural wave functions to the solution expansion are comprised in the ðns £ 1Þ vector {ws }: ^ F is a particular solution for the mechanical Function w force term in the right-hand side of Eq. (3). Several mathematical expressions may serve as particular solution. It is advantageous, however, to select an expression, which is already close to the physical response of the plate. Therefore, the displacement of an infinite plate, excited by a normal point force F is selected jF ð22Þ ½H ð2Þ ðkb rÞ 2 H0ð2Þ ð2jkb rÞ; 8kb2 D 0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi with r ¼ ðj 2 jF Þ2 þ ðh 2 hF Þ2 and where H0ð2Þ is the zero-order Hankel function of the second kind. ^ a ðj; hÞ in the ð1 £ na Þ vector ½w^ a is a Each function w particular solution for the part of the cavity pressure loading term in the inhomogeneous right-hand side of Eq. (3), resulting from one of the acoustic wave functions Fa in the cavity pressure expansion (16). The displacement of an infinite plate, excited by this distributed pressure, may be used as particular solution function. However, due to the flat shape of the plate, a particular solution function may be defined, which is proportional to the pressure loading function
w^ F ðj; hÞ ¼ 2
w^ a ðj; hÞ ¼
e2jðkxa ðux 2fxa Lx Þþkya ðuy 2fya Ly Þþkza ðuz 2fza Lz ÞÞ ; D½A4j þ 2A2j A2h þ A4h 2 kb4
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
531
with Aj ¼ kxa ux;j þ kya uy;j þ kza uz;j ;
ð23Þ
Ah ¼ kxa ux;h þ kya uy;h þ kza uz;h :
2.3. Coupled vibro-acoustic wave model Due to the particular choice of the field variable expansions (16) and (19), the governing dynamic Eqs. (2) and (3) are exactly satisfied, irrespective of the values of the unknown wave function contributions pa and ws. These contributions are merely determined by the acoustical and structural boundary conditions. To obtain a finite sized wave model, which is amenable to numerical implementation, these boundary conditions are transformed into a weighted residual formulation. 2.3.1. Acoustic boundary conditions Depending on the wave function contributions pa and ws in the field variable expansions (16) and (19), some approximation errors are induced in the representation of the acoustic boundary conditions (6) – (9). For each of these boundary conditions, an associated residual error function may be defined as follows Rp ðrÞ ¼ p^ ðrÞ 2 p ðrÞ;
r [ Vp ; r [ Vv ;
j ›p^ ðrÞ p^ ðrÞ ; 2 RZ ðrÞ ¼ ZðrÞ r0 v ›n
r [ VZ ;
j ›p^ ðrÞ ^ 2 jvwðrÞ; r0 v ›n
2
VZ
r [ Vs :
ð Vp
Vs
2.3.2. Structural boundary conditions In the similar way as for the acoustic boundary conditions, some residual error functions for the structural boundary conditions (13) – (15) may be defined as follows ^ s Þ 2 wðr Rw ðrs Þ ¼ wðr s Þ;
rs [ ðGwu < Gwm Þ;
ð30Þ
r s [ G wu ;
ð31Þ
^ s Þ 2 un ðrs Þ; Ru ðrs Þ ¼ Lu ½wðr ^ s Þ 2 m Rm ðrs Þ ¼ Lm ½wðr n ðrs Þ;
rs [ ðGmQ < Gwm Þ; ð32Þ
n ðrs Þ; ^ s Þ 2 Q RQ ðrs Þ ¼ LQ ½wðr
rs [ GmQ :
ð33Þ
ð25Þ ð26Þ ð27Þ
A weighted residual formulation of the acoustic boundary conditions is defined as ð ð ð p~ Rv dV þ p~ RZ dV þ p~ Rs dV Vv
matrices ½Aaa and ½Caa and the ðna £ ns Þ matrix ½Cas ; the reader is referred to Ref. [2].
ð24Þ
j ›p^ ðrÞ 2 v n ðrÞ; Rv ðrÞ ¼ r0 v ›n
Rs ðrÞ ¼
Fig. 3. Normal and tangential directions at a corner point.
ð28Þ
j ›p~ R dV ¼ 0; r0 v ›n p
where p~ is a weighting function. In a similar way as in the Galerkin weighting procedure, used in the finite element method, each of the acoustic wave functions Fa in the cavity pressure expansion (16) is used as a weighting function p~ in Eq. (28). This yields a set of na algebraic equations in the ðna þ ns Þ unknown wave function contributions pa and ws ( )
pa ¼ {fa }: Aaa þ Caa Cas ð29Þ ws For a full derivation of the coefficients in the ðna £ na Þ
For the sake of model symmetry, some additional residuals are defined at the corner points of the plate boundary curve Gs, at which the normal and tangential directions are not uniquely defined (Fig. 3). For the nw corner points, which belong to the boundary part Gwu < Gwm ; the additional residual is the approximation error for the out-of-plane plate displacement at the corner point location rsc ^ sc Þ 2 wðr Rcw ðrsc Þ ¼ wðr sc Þ;
ðc ¼ 1::nw Þ:
ð34Þ
For the nF corner points, which belong to the boundary part GmQ, the additional residual is the approximation error for the concentrated corner point force, which is associated with the discontinuity of the torsional moment ^ sc Þ 2 ðm RcF ðrsc Þ ¼ LF ½wðr ns ðrþ ns ðr2 sc Þ 2 m sc ÞÞ ^ sc Þ ^ sc Þ ›2 wðr ›2 wðr ¼ 2Dð1 2 nÞ 2 þ þ 2 ›gn ›gs ›gn ›g2 s 2
ðm ns ðrþ sc Þ
2
! ð35Þ
m ns ðr2 sc ÞÞ;
ðc ¼ 1::nF Þ: A weighted residual formulation of the structural boundary
532
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
conditions is defined as ð Gwu
2
LQ ½wR ~ w dG þ
ð Gwm
þ
nw X c¼1
associated uncoupled plate problem without acoustic pressure loading from the cavity.
ð Gwu
Lm ½wR ~ u dG
Lu ½wR ~ m dG 2
LF ½wðr ~ sc ÞRcw ðrsc Þ 2
2.4. Convergence ð GmQ nF X
wR ~ Q dG
wðr ~ sc ÞRcF ðrsc Þ ¼ 0;
ð36Þ
c¼1
where w ~ is a weighting function. Again, by using each of the structural wave functions Cs in the displacement expansion (19) as a weighting function w~ in Eq. (36), a set of ns algebraic equations is obtained in the ðna þ ns Þ unknown wave function contributions pa and ws ( )
pa ¼ {fs }: Csa Ass ð37Þ ws For a full derivation of the coefficients in the ðns £ ns Þ matrix ½Ass and the ðns £ na Þ matrix ½Csa ; the reader is referred to Ref. [2]. 2.3.3. Model properties The combination of the weighted acoustic boundary condition (29) and the weighted structural boundary condition (37) yields the coupled vibro-acoustic model " #( ) ( ) Aaa þ Caa Cas fa pa ¼ : ð38Þ Csa Ass ws fs Solving model (38) for the unknown wave function contributions pa and ws and back-substituting the results in expansions (16) and (19) yield the prediction of the steadystate dynamic response of the coupled vibro-acoustic problem. Several model properties may be identified. † Since the acoustic and structural wave functions Fa and Cs are defined in, respectively, the entire acoustic domain V and the entire structural domain Vs, the model matrices are fully populated. † Since the acoustic and structural wave functions Fa and Cs are complex functions and since they are implicitly dependent on frequency by virtue of the frequency dependence of their wave number components, the model matrices are complex and frequency dependent and cannot be partitioned into frequency independent submatrices. † The model matrix in Eq. (38) is not symmetric. However, it is proved in Ref. [2] that the submatrices ½Aaa and ½Ass are symmetric. Note that these submatrices would become, respectively, the model matrix for the associated uncoupled acoustic problem, having no elastic boundary surface Vs, and the model matrix for the
Since Eqs. (17) and (20) have an infinite number of real and complex solutions, the key question is whether a set of acoustic and structural wave functions can be selected from the infinite number of possible wave functions, such that the subsequent pressure and displacement expansions (16) and (19) converge towards the exact solution. 2.4.1. Convex problem domain The acoustic wave functions Fa with the following wave number components are proposed ffi! m1x p m1y p qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ðkxa ; kya ; kza Þ ¼ ; ; ; ^ k2 2 kxa 2 kya Lx Ly ffi m2x p qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 k2 ; m1z p ; ; ^ k2 2 kxa and za Lx Lz qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m p m p ! 2y 2 2 k2 ; ; ; 2z ^ k2 2 kya za Ly Lz
ð39Þ
with mix ; miy ; miz ¼ 0; ^1; ^2; … ði ¼ 1; 2Þ: The structural wave functions Cs with the following wave number components are proposed ! sx p qffiffiffiffiffiffiffiffiffiffi 2 2 ðkjs ; khs Þ ¼ ; ^ kb 2 kjs ; Lj ! qffiffiffiffiffiffiffiffiffiffi sx p 2 2 ; ^j kb þ kjs ; and Lj qffiffiffiffiffiffiffiffiffiffiffi s p ! qffiffiffiffiffiffiffiffiffiffi s p ! y y 2 2 ; ^j kb2 þ kh2 s ; ; ^ kb 2 khs ; Lh Lh
ð40Þ
with sx, sy ¼ 0; ^1; ^2; … Recall that Lx, Ly and Lz are the dimensions of the (smallest) rectangular prism, enclosing the acoustic cavity domain V, and that Lj and Lh are the dimensions of the (smallest) rectangle, enclosing the plate domain Vs. It is proved in Ref. [2] that the convexity of the acoustic domain V and the plate domain Vs is a sufficient condition for the solution expansions (16) and (19), using the proposed wave function selections (39) and (40), to converge, in the limit for mix, miy, miz, sx, sy : 0 ! 1; towards the exact coupled vibro-acoustic response. Moreover, it can be proved that the proposed acoustic wave
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
533
Fig. 4. Two-dimensional acoustic validation example (convex domain).
functions may be combined into !
m1y p 2jkza ðz2fza Lz Þ m p e Fa ðx; y; zÞ ¼ cos 1x cos ; Lx Ly m2x p 2jkya ðy2fya Ly Þ m1z p cos e ; and cos Lx Lz ! m2y p m2z p 2jkxa ðx2fxa Lx Þ cos cos e ; Lz Ly
ð41Þ
with mix, miy, miz ¼ 0; 1; 2; …; ði ¼ 1; 2Þ: Note that in this ^ a in Eq. (23) must be case, the particular solution functions w modified, as described in Ref. [2]. To illustrate the convergence of the proposed wave function selection, based on the definition of an enclosing rectangular domain, a two-dimensional acoustic validation problem is considered, as shown in Fig. 4. At the right boundary of a slender convex acoustic cavity (length of 5 m and height of 1 m), a real normal impedance of 2000 kg/(m2 s) is prescribed. The upper and lower parts of the boundary surface are assumed to be rigid. The dynamic excitation of the air-filled cavity consists of a unit prescribed normal velocity, imposed at the left boundary. The air in the cavity has a density r0 ¼ 1.225 kg/m3 and a speed of sound c ¼ 340 m/s. For this two-dimensional example, the two-dimensional counterpart of the acoustic wave function selection, proposed in Eq. (39), should be used, i.e. qffiffiffiffiffiffiffiffiffiffi mx p 2 ; ðkxa ; kya Þ ¼ and ; ^ k2 2 kxa Lx ð42Þ qffiffiffiffiffiffiffiffiffiffi m p ! y 2 2 ; ^ k 2 kya ; Ly with Lx and Ly the dimensions of an enclosing rectangle ðmx ; my ¼ 0; ^1; ^2; …Þ: Fig. 5 shows the acoustic frequency response function (FRF), relating the pressure at the response location (indicated in Fig. 4) to the unit velocity input, obtained from a wave model with a total of 120 acoustic wave functions. The enclosing rectangle, used for the wave function selection, is indicated in dashed lines in Fig. 4. To compare the performance of the new prediction technique with the finite element method, three FE models are
constructed for the considered validation example. In these models, the two-dimensional acoustic cavity is discretised into 4-noded linear elements, resulting in FE model sizes of, respectively, 105, 369 and 1292 nodal pressure degrees of freedom. The corresponding acoustic FRFs, obtained from these three FE models, are also plotted (in dashed lines) in Fig. 5. This figure clearly illustrates that, by increasing the number of elements, the prediction accuracy increases and the finite element results gradually converge towards the results, obtained from the small wave model. This figure illustrates also the typical approximation errors, involved with FE discretisations for dynamic analysis. The use of simple, but approximating (polynomial)
Fig. 5. Acoustic FRF at response point (solid: wave model—120 dof, dashed: FE models (upper: 105 dof, middle: 369 dof, lower: 1292 dof)).
534
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
Fig. 6. Acoustic convergence curve at 100 Hz (accuracy versus dof) (wave models (solid)—FE models (dashed)).
shape functions, describing the dynamic response within each element, involves not only an interpolation error, but also a numerical dispersion error, in that the wave number of the numerical solution differs from the wave number of the exact problem solution [14]. This error results usually in an overestimation of the resonance frequencies. Since the dispersion error increases for increasing frequency, the mesh density and subsequent computational effort, required for getting reasonable prediction accuracy, increase also with frequency. Another indication of the convergence of the proposed technique is given by the convergence curve, shown in Fig. 6. This curve plots the relative prediction error on the pressure at the considered response point against the number of degrees of freedom. The relative error is defined as the difference between the calculated and
Fig. 7. Two-dimensional vibro-acoustic validation example.
exact solution, divided by the exact solution. However, since the exact solution for the considered problem is not known analytically, the result, obtained from a wave model with a very large number of wave functions, is used instead to define the relative prediction error. The convergence curves, associated with both the novel wave based prediction technique (solid line) and with the finite element method (dashed line), are plotted in Fig. 6 for the case of a dynamic excitation at 100 Hz. This figure illustrates not only that the novel technique provides a higher accuracy with a smaller number of degrees of freedom, but also that its convergence rate, i.e. the slope of the convergence curve, is substantially higher. The size of the enclosing rectangle has some influence on the overall accuracy, as is illustrated for the twodimensional vibro-acoustic validation problem, shown in Fig. 7. One part of the boundary surface of an air-filled (r0 ¼ 1.225 kg/m3, c ¼ 340 m/s) cavity (length of 1.5 m and height of 1 m) consists of a 0.75 m long, flexible plate (E ¼ 70 £ 109 N/m2 , rs ¼ 2790 kg/m3 , n ¼ 0.3, t ¼ 0.002 m) with clamped boundaries, while the remaining part of the cavity boundary surface is perfectly rigid. The coupled vibro-acoustic system is excited by a unit normal line force F at j F ¼ 0.25 m with a harmonic time dependence at 200 Hz. As shown in Fig. 8, two different enclosing rectangles have been used for the definition of the wave number components. This figure shows the resulting convergence curves for the out-of-plane plate displacement at the excitation point. It illustrates that it is beneficial for both the accuracy and the convergence rate to select the wave functions, associated with the smallest rectangular domain, enclosing the considered problem domain. 2.4.2. Substructuring So far, the convergence of the novel prediction technique is discussed for convex problem domains. Although a firm mathematical proof is still under construction, the authors have not yet found an example, counterproving the statement that any type of non-convex cavity or plate domain can be decomposed into some subdomains and that a convergent pressure or displacement expansion for each subsequent subdomain is obtained by using a wave function selection of type (39) or (40), which is based on the dimensions of a rectangular domain, enclosing the subdomain. When the acoustic domain V must be decomposed for convergence reasons, a weighted residual formulation of type (28) must be specified for each acoustic subdomain, resulting in a wave model of type (29) for each subdomain. In this respect, some additional boundary conditions must be specified to ensure the continuity of the pressure and normal fluid velocity along the common boundary interfaces between the various subdomains. Hence, for each pair of adjacent acoustic subdomains, a pressure residual function of type (24), specifying the numerical pressure difference along the common interface, must be added to the
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
535
Fig. 10. Substructuring of the concave domain.
Fig. 8. Influence of the size of the enclosing rectangular domain on the convergence.
weighted residual formulation, associated with one subdomain, and a normal velocity residual function of type (25), specifying the numerical normal velocity difference along the common interface, must be added to the weighted residual formulation, associated with the other subdomain. In the same way, when the plate domain Vs must be
Fig. 9. Two-dimensional acoustic validation example (concave domain).
decomposed, a weighted residual formulation of type (36) must be specified for each structural subdomain, resulting in a wave model of type (37) for each subdomain. For each pair of adjacent structural subdomains, a residual function of type (30), specifying the numerical out-of-plane displacement difference along the common interface, and a residual function of type (31), specifying the numerical difference in normal rotational displacement along the interface, must be added to the weighted residual formulation, associated with one subdomain. In addition, residual functions of type (32) and (33) must be added to the weighted residual formulation, associated with the other subdomain, to ensure the balances of the generalised normal shear forces and the normal bending moments along the common interface. The total coupled vibro-acoustic model for a non-convex problem is then obtained from combining the weighted residual formulations of all acoustic and all structural subdomains. The motivation for the proposed substructuring procedure for concave problems is similar to the one in the T-element approach, discussed in Section 1, in that convergence is quested by dividing the problem domain. However, due to the advantageous properties of the proposed wave function selections (39) and (40), the problem domain must only be divided in a few, large subdomains, rather than in many small elemental subdomains as is the case in the T-element approach. To illustrate the convergence of the proposed substructuring procedure, a two-dimensional acoustic validation problem is considered, as shown in Fig. 9. At the right boundary of an L-shaped acoustic cavity (outer length of 4 m and height of 2 m), a real normal impedance of 2000 kg/(m2 s) is prescribed. The dynamic excitation of the airfilled cavity consists of a unit prescribed normal velocity, imposed at the upper part of the boundary surface. The remaining boundary parts are assumed to be rigid. The air in the cavity has a density r0 ¼ 1.225 kg/m3 and a speed of sound c ¼ 340 m/s. For this concave problem, the L-shaped cavity is divided in two subdomains, as shown in Fig. 10. For both subdomains, a wave function selection of type (42) is used. Fig. 11 shows the acoustic FRF, relating the pressure at the response location (indicated in Fig. 9) to the unit velocity input, obtained from a wave model with a total of 108 wave functions. To compare the performance of the prediction technique with the finite element method, three
536
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
acoustic wavelengths at the considered excitation frequency. When the excitation frequency v coincides with an undamped structural eigenfrequency of the circumscribing rectangular plate domain with simply supported boundaries, i.e. !sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 q22 Et2 2 q1 v¼p þ 2 ; 2 Lj Lh 12ð1 2 n2 Þrs ð44Þ ðq1 ; q2 ¼ 1; 2; …Þ; or when it coincides with an undamped acoustic eigenfrequency of the circumscribing rectangular cavity with rigid side walls, i.e. vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !2 u u a1 p 2 a2 p a3 p 2 t ; v¼c þ þ Lx Ly Lz ð45Þ ða1 ; a2 ; a3 ¼ 0; 1; …Þ; it follows from Eqs. (39) and (40) that not all wave functions are linearly independent and that some wave functions must be eliminated to restore the linear independence without reducing the span of the wave function sets.
3. Comparison with the finite element method Fig. 11. Acoustic FRF at response point (solid: wave model—108 dof, dashed: FE models (upper: 105 dof, middle: 217 dof, lower: 793 dof)).
FE models are constructed for this validation problem using 4-noded linear elements, resulting in model sizes of, respectively, 105, 217 and 793 nodal pressure degrees of freedom. The corresponding acoustic FRFs, obtained from these three FE models, are also plotted (in dashed lines) in Fig. 11. Again, this figure illustrates that, by increasing the number of elements, the prediction accuracy increases and the finite element results gradually converge towards the results, obtained from the small wave model. 2.4.3. Truncation rule For the numerical implementation of the coupled wave model (38), the infinite wave function series (mix, miy, miz, sx, sy : 0 ! 1) must be truncated at some finite integer values m ix ; m iy ; m iz ; s x ; sy : In this respect, the following truncation rule is proposed m sy iy m m s 2maxðk; kb Þ ix : < < iz < x < $ p Lx Ly Lz Lj Lh
ð43Þ
This truncation rule expresses that the smallest wavelength components in the structural and acoustic wave function sets may not be larger than half the structural and half the
In comparison with the finite element method, the new prediction technique has several advantageous properties. † Since the field variable expansions satisfy a priori the governing dynamic equations, approximation errors are only involved with the representation of the boundary conditions. This yields substantially smaller models than corresponding FE models. † Predictions of dynamic quantities such as structural stresses and strains or fluid velocities are obtained from deriving the prediction results of the primary field variables, i.e. the structural displacement components or the acoustic pressure. In the FE method, the latter are usually expressed in terms of polynomial expansions, which are complete up to a certain order. Hence, the derived secondary variables are also expressed in terms of polynomial expansions, which are, however, only complete up to a lower-order. As a result, the lower-order (polynomial) expansions of the derived secondary variables represent a smaller spatial variation and are less accurate than the expansions of the corresponding primary variables. In the proposed prediction technique, however, the spatial gradients of the wave functions are wave functions, which have the same spatial variation as the ones, used in the primary field variable expansions. This is advantageous for the convergence rate, especially in the case of coupled vibro-acoustic problems, for which
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
the effect of the fluid on the structure is pressure controlled, but for which the effect of the structure on the fluid is velocity controlled. † A model refinement is easily obtained. Firstly, the accuracy of a model is easily checked, since only the approximation errors of the boundary conditions must be verified. Secondly, increasing the number of degrees of freedom, i.e. the number of wave functions in the field variable expansions, requires only the calculation of the additional matrix coefficients, associated with the additional wave functions, while the original model coefficients remain unchanged. For FE models, increasing the number of degrees of freedom requires the calculation of a completely new model, at least when a global refinement of the element discretisation is needed. The new prediction technique has also some drawbacks, compared with the FE method. † Since the wave functions are complex functions, defined within the entire continuum domain or at least within large subdomains, the wave model matrices are fully populated and complex. † Due to the implicit frequency dependence of the wave functions, the wave model matrices cannot be decomposed into frequency independent submatrices, which implies that a wave model must be recalculated for each frequency of interest and that natural frequencies and mode shapes cannot be obtained from a standard eigenvalue calculation. † Due to the global and oscillatory nature of the wave functions and due to the poor condition of a wave model, the involved numerical integrations are computationally more demanding. † The FE method is a widely applicable technique, since the convergence of a FE model is obtained with simple element shape functions, whereas the convergence of a wave model requires the definition of a complete wave function set and appropriate particular solution functions.
537
Fig. 12. Three-dimensional vibro-acoustic validation example.
The acoustic cavity is comprised in the volume of an enclosing rectangular prism with dimensions Lx ¼ 1:5 m; Ly ¼ 0:5 m; Lz ¼ 1 m: The vibro-acoustic system is excited by a unit mechanical point force F, applied at location ðjF ; hF Þ ¼ ð2Lj =3; Lh =3Þ on the plate, in the direction normal to the plate and with a harmonic time dependence at a certain frequency v. Fig. 13 shows the real part of the out-of-plane plate displacement at 200 Hz, obtained from a wave model with 314 (na ¼ 258 acoustic and ns ¼ 56 structural) wave functions, and the corresponding normal fluid displacement along the cavity – plate coupling interface. This figure indicates that the wave model enables an accurate representation of the clamped plate boundary conditions of type (13) and the normal displacement continuity condition of type (9). Fig. 14 shows the contour plot of
4. Three-dimensional validation example To illustrate the accuracy of the proposed prediction technique for three-dimensional coupled vibro-acoustic analysis and to assess its computational performance in comparison with the finite element method, the numerical validation example is considered, as shown in Fig. 12. One boundary surface of an acoustic cavity consists of a flexible rectangular ðLj ¼ 0:75 m; Lh ¼ 0:5 mÞ plate with clamped boundaries, while all other cavity boundary surfaces are perfectly rigid. The air in the cavity has a density r0 ¼ 1:225 kg=m3 and a speed of sound c ¼ 340 m=s: The undamped aluminium plate has a thickness t ¼ 2 £ 1023 m; a density rs ¼ 2790 kg=m3 ; a Poisson coefficient n ¼ 0:3 and an elasticity modulus E ¼ 70 £ 109 N=m2 :
Fig. 13. Real part of the out-of-plane plate displacement at 200 Hz (a) and the corresponding normal fluid displacement along the cavity – plate coupling interface (b).
538
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
Table 1 Properties of the FE models
FE model 1 FE model 2 FE model 3
No. of structural elements
No. of acoustic elements
dof
54 216 486
1296 10 368 34 992
1849 12 586 39 991
the real part of the acoustic pressure in the cavity plane y ¼ Ly =3: Since the pressure contour lines are perpendicular to the rigid cavity boundary surfaces, the zero normal fluid velocity boundary conditions of type (7) at these boundary surfaces are accurately represented. Since the solution expansions (16) and (19) inherently satisfy the governing dynamic equations, and since all boundary conditions of the considered problem are accurately represented, this validation example illustrates that the proposed prediction technique is practically converging and enables an accurate approximation of the exact three-dimensional vibro-acoustic response. To compare the performance of the new prediction technique with the finite element method, three FE models are constructed for the considered validation example. In these models, the plate is discretised into four-noded quadrilateral shell elements and the acoustic cavity is discretised into eight-noded hexahedral fluid elements. Table 1 lists the number of elements in each model and their corresponding number of unconstrained degrees-offreedom (dof). Fig. 15 compares the direct frequency response functions (with a frequency resolution of 1 Hz), obtained from the aforementioned wave model with 314 wave functions and from the three FE models. Fig. 16 compares the associated FRFs of the pressure at location ðx; y; zÞ ¼ ð2Lx =3; 2Ly =3; Lz =4Þ: Again, these figures clearly shown that, by increasing the number of elements, the prediction accuracy increases and the finite element results gradually
converge towards the results, obtained from the small wave model. In order to assess the convergence rate, several wave models have been constructed, as listed in Table 2. The solid lines on Fig. 17 indicate the resulting relative prediction errors for the out-of-plane plate displacement at the excitation point ðjF ; hF Þ ¼ ð2Lj =3; Lh =3Þ and for the cavity pressure at ðx; y; zÞ ¼ ð2Lx =3; 2Ly =3; Lz =4Þ for two excitation frequencies, i.e. 60 and 180 Hz. The dashed lines indicate the corresponding convergence curves for the finite element models, listed in Table 1. This figure shows that the new prediction technique yields accurate prediction results with substantially smaller models than the finite element method and with a higher convergence rate.
Fig. 14. Contour plot of the real part of the acoustic pressure in the cavity plane y ¼ Ly =3 at 200 Hz.
Fig. 15. Direct FRF (solid: wave model, dashed: FE model 1 (a), FE model 2 (b) and FE model 3 (c)).
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
539
Table 2 Properties of the wave models
Wave model 1 Wave model 2 Wave model 3
No. of structural wave functions
No. of acoustic wave functions
dof
32 56 80
78 258 542
110 314 622
However, as indicated in Section 3, the beneficial effect of a smaller model size on the computational effort is partly annihilated by the fact that wave models are fully populated and cannot be decomposed into frequency-independent submatrices. Therefore, it is more appropriate to consider the prediction accuracy in terms of required CPU time rather than model size, as shown in Fig. 18. The CPU times in this figure denote the times, needed for a direct response calculation at one frequency. The indicated CPU times for the wave models include both the times for constructing the
Fig. 17. Structural (left) and acoustic (right) convergence curves (accuracy versus dof) (wave models (solid)—FE models (dashed)) (þ : 60 Hz, L: 180 Hz).
model as well as for solving the resulting matrix equation. Since FE models consist of frequency-independent submatrices, these submatrices must only be calculated once and the construction of the models at each frequency requires only a simple submatrix assembly. As a result, the computational effort for constructing the models is negligible, compared with the efforts for solving the (large) models, at least when the dynamic response is calculated in a large frequency range. Therefore, the indicated CPU times for the FE models comprise only the direct solution times. All calculations are performed on a HP-C180 workstation (SPECfp95 ¼ 18.7, SPECint95 ¼ 11.8). Fig. 18 clearly shows the beneficial convergence rate of the new prediction technique. In comparison with the finite element method, the technique provides highly accurate predictions of the coupled vibro-acoustic response with a
Fig. 16. Acoustic FRF at cavity location ðx; y; zÞ ¼ ð2Lx =3; 2Ly =3; Lz =4Þ (solid: wave model, dashed: FE model 1 (a), FE model 2 (b) and FE model 3 (c)).
Fig. 18. Structural (left) and acoustic (right) convergence curves (accuracy versus CPU time) (wave models (solid)—FE models (dashed)) (þ: 60 Hz, L: 180 Hz).
540
W. Desmet et al. / Advances in Engineering Software 33 (2002) 527–540
substantially smaller computational effort. Note also that the high computational efficiency will most likely become even more apparent, when the technique is implemented in an efficient programming environment instead of the currently used MATLAB environment.
5. Conclusion A new prediction technique for coupled vibro-acoustic analysis has been developed. The technique is based on the indirect Trefftz method and approximates the structural and acoustic field variables of a coupled vibro-acoustic system in terms of solution expansions, which exactly satisfy the governing dynamic equations. The contribution factors in these expansions are obtained from a weighted residual formulation of the boundary conditions. The paper describes the basic concept of the new technique for the modelling of the vibro-acoustic coupling between the pressure field in an acoustic cavity with arbitrary shape and the out-of-plane displacement of a flat plate with arbitrary shape. It is illustrated through a numerical validation example that, in comparison with the finite element method, the new prediction method provides accurate prediction results with a substantially smaller computational effort. Due to this beneficial convergence rate, the proposed prediction technique may be a suitable method for obtaining accurate coupled vibro-acoustic predictions up to much higher frequencies than the conventional finite element method.
Acknowledgements Wim Desmet is a Postdoctoral Fellow of the Fund for Scientific Research, Flanders (Belgium) (FWO). The research work of Bas van Hal is financed by a scholarship
of the Institute for the Promotion of Innovation by Science and Technology in Flanders (IWT).
References [1] Craggs A. An acoustic finite element approach for studying boundary flexibility and sound transmission between irregular enclosures. J Sound Vibr 1973;30:343–57. [2] Desmet W. A wave based prediction technique for coupled vibroacoustic analysis. PhD Dissertation. K.U. Leuven, Leuven; 1998. [3] Trefftz E. Ein Gegenstu¨ck zum Ritzschen Verfahren. Proc 2nd Int Cong Appl Mech, Zurich 1926;131–7. [4] Kita E, Kamiya N. Trefftz method: an overview. Adv Engng Software 1995;24:3–12. [5] Herrera I. Boundary methods: an algebraic theory. Boston: Pitman Advanced Publishing Program; 1984. [6] Masson P, Redon E, Priou JP, Gervais Y. The application of the Trefftz method to acoustics. In: Crocker M, editor. Proceedings of the Third International Congress on Air- and Structure-borne Sound and Vibration, Montreal; 1994. p. 1809–16. [7] Langley RS. Some perspectives on wave-mode duality in SEA. In: Fahy FJ, Price WG, editors. Solid mechanics and its applications. IUTAM Symposium on Statistical Energy Analysis, Dordrecht: Kluwer; 1997. p. 1– 12. [8] Jirousek J, Wroblewski A. T-elements: state-of-the-art and future trends. Arch Comput Meth Engng 1998;41(5):831– 49. [9] Zielinski AP, Zienkiewicz OC. Generalized finite element analysis with T-complete boundary solution functions. Int J Numer Meth Engng 1985;21:509 –28. [10] Jirousek J, Wroblewski A. Least-squares T-elements: equivalent FE and BE forms of a substructure-oriented boundary solution approach. Commun Numer Meth Engng 1994;10:21–32. [11] Piltner R. Special finite elements with holes and internal cracks. Int J Numer Meth Engng 1985;21:1471–85. [12] Stojek M. Least-squares Trefftz-type elements for the Helmholtz equation. Int J Numer Meth Engng 1998;41:831–49. [13] Texeira de Freitas JA. Hybrid finite element formulations for elastodynamic analysis in the frequency domain. Int J Solids Struct 1998;36:1883–923. [14] 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 1999;46:471 –99.