# A Multidomain Spectral Method for Analysis of Interior Vibroacoustic Systems with Segmented Boundaries

Authors: Zhenguo Zhang; Haiting Yu; Ningyuan Duan; Hongxing Hua
Journal: Shock and Vibration
 (2020)

---

## Abstract

Spectral methods have previously been applied to analyze a multitude of vibration and acoustic problems due to their high computational efficiency. However, their application to interior structural acoustics systems has been limited to the analysis of a single plate coupled to a fluid-filled cavity. In this work, a general multidomain spectral approach is proposed for the eigenvalue and steady-state vibroacoustic analyses of interior structural-acoustic problems with discontinuous boundaries. The unified formulation is derived by means of a generalized variational principle in conjunction with the spectral discretization procedure. The established framework enables one to easily accommodate complex systems consisting of both a structure assembly and a built-up cavity with moderate geometric complexities and to effectively analyze vibroacoustic behaviors with sufficient accuracy at relatively high frequencies. Two practical examples are chosen to demonstrate the flexibility and efficiency of the proposed formulation: a built-up cavity backed by an assembly of multiple connected plates with arbitrary orientations and a thick irregular elastic solid coupled with a heavy acoustic medium. Comparison to finite element simulations and convergence studies for these two examples illustrate the considerable computational advantage of the method as compared to finite element procedures.


---

## Body

## 1. Introduction

Interior vibroacoustic interaction systems (i.e., systems in which an acoustic fluid is enclosed by a vibrating structure) with segmented or discontinuous boundaries provide fundamental components of many commonly used technological and scientific applications [1], including acoustic quality engineering, submerged turbomachinery, biomechanical systems, and aircraft cabins. Effective prediction or evaluation of vibroacoustic behaviors has long been regarded as one of the most significant issues for the design and control of such systems [2]. Due to increasingly complicated structural-acoustic systems and the potential mathematical complexities arising in the solution of fluid-structural interaction problems with segmented or discontinuous boundaries, development of generally applicable and efficient prediction techniques will remain a challenge and is anticipated to become even more significant.Considerable attention has been focused on vibroacoustic behaviors of acoustic fluid-structural interaction systems, and researchers have adopted different computational methods involving element-based methods and analytical methods. The finite element method (FEM) and boundary element method (BEM) are the most widely utilized numerical procedures [3–6] that have been used to analyze vibroacoustic problems. The superiority of element-based solutions lies in the excellent applicability and flexibility for the analysis of acoustic systems with complicated geometries and complex boundary conditions; they are probably the most practical methods for many real-life problems that are full of small details [7]. However, the trade-off between prediction accuracy and efficiency limits these methods to low-frequency ranges because of the rapidly growing model size and CPU requirement for increasing frequencies [8, 9]. Analytical and semianalytical methods, such as the modal coupling theorem [10, 11], the impedance method [12], the finite strip method [13], the Rayleigh–Ritz method [14–16], and the Trefftz methods [17, 18], have also been used to analyze structural acoustics and vibration behaviors. The attractive features of these models rely on their flexibility for extensive parametric studies and their ability to offer a clear physical insight into vibroacoustic phenomena as well as the reasonable CPU requirements. However, most analytical work has concentrated on a relatively simple system consisting of a flexible vibrating panel backed by a regular or irregular acoustic cavity, which is often not representative of actual vibroacoustic systems encountered in engineering applications. So there is still a strong need to increase the applicability of analytical methods for more complex vibroacoustic problems in order to fully take benefit of the strengths of analytical approaches. This study attempts to demonstrate the capability of analytical techniques to accommodate built-up structural-acoustic systems of moderate complexity.The spectral method (SM), which constructs approximate solutions from the space of admissible functions, can reasonably yield model of small size for the considered systems [19]. By employing orthogonal families of polynomials, which possess complete, linearly independent, and highly smooth properties, as admissible functions, SM can achieve rapid convergence and superior stability in numerical operations [20]. Moreover, in conjunction with blending function [21], the Lagrange multipliers technique [22], or the artificial spring technique [15], the SM can be further extended to deal with problems involving complex arbitrary geometries or boundary conditions. However, it cannot directly scale to the combined systems modeled with more than one component since mathematical difficulties arising from the explicit representation of discontinuous boundary conditions. Thus, it is a good approach to take advantage of the fundamental idea of dynamic substructuring [23, 24], which was developed for the dynamic simulation of complex structural systems. Separating the whole combined system properly into several subdomains according to geometries and boundaries, these techniques synthesize admissible functions (or modal vectors) of individual subdomains to form governing equations of the whole system by imposing certain compatibility on internal boundaries. Such multidomain approaches can dramatically simplify the analytical characterization, decrease the computation time, and also shows the potential to tackle dynamic problems at relatively high frequencies; that is why it is adopted in the present study.Although the multidomain approaches based on analytical techniques have been used in structural dynamic applications, few attempts have been made to use them for fluid-structure interaction problems. Desmet et al. [25] introduced the multidomain wave-based method for vibroacoustic problems based on an indirect Trefftz approach, where wave functions are adopted to describe the internal dynamic field of subdomains and the weighted residual approach is applied to enforce the continuity and boundary conditions. The method has been proven to be successful for vibroacoustic problems with relatively moderate geometry complexity in a considerable frequency range [17]. Nevertheless, the system matrices are always complex and frequency dependent, so that these matrices need to be reconstructed for every frequency of interest and a standard eigenvalue calculation cannot be efficiently carried out. Cheng [26] and Amabili [27] used the Rayleigh–Ritz based substructuring approach to study vibroacoustic behaviors of a plate-shell combined structure filled with fluids, where the artificial spring technique is employed to combine the substructures. However, the accuracy of this technique is sensitive to the choice of spring parameters, and the method becomes cumbersome as the number of substructures increases since an increasing number of spring systems will be needed. Xie et al. [16] employed a domain decomposition method to describe vibroacoustic responses of an irregular acoustic cavity enclosed by a single flat plate. However, their formulation of structural model is limited to transverse vibrations of a single thin plate based on Kirchoff bending theory; the analysis of more complex vibroacoustic systems that include an assembly of multiple plates (which requires to consider the close coupling of vibrations in different directions), or a thicker structure (for which thin plate theory is not valid), is not directly possible in their work.Motivated by overcoming the aforementioned limitations, the objective of this study is to propose a general multidomain spectral method for eigenvalue analysis and steady-state response analysis of interior vibroacoustic systems with moderately complex geometries and segmented boundaries. The novelty of this work stems from the unified formulation of such system models that include both the structure assembly and built-up cavity of different dimensionalities. The approach makes it possible to analyze structures that consist of an assembly of plates of arbitrary directions or thick structures that require using general elasticity theory. The development begins by partitioning a physical space, which is either relatively complicated or large to be described by a single domain, into several nonoverlapping substructures and fluid subdomains according to the distribution of boundaries. A generalized variational principle [28] together with the spectral discretization procedure [20] is then extended to describe interaction dynamics between the substructures or/and subdomains. The proposed approach combines the versatility of the spectral method and the flexibility of multidomain synthesis by inheriting from the dynamic substructuring for structural dynamics [24]. It is shown from calculation verification and convergence studies that the developed approach exhibits tremendous computational advantages and allows reaching higher frequency range than FEM.


## 2. Statement of the Problem

The general vibroacoustic problem studied here is schematically sketched in Figure1, where the fluid domain Ωf is partially in contact with the structural domain Ωs at the interface Γfs. The structure is supposed to be linear-elastic, and the structural boundaries with prescribed displacements and surface forces are denoted by Γus and Γσs, respectively. The compressible fluid is assumed to be homogeneous, irrotational, and inviscid, and Γf (Γf = Γpf ∪ Γvf ∪ Γzf) represents the fluid boundaries, where three acoustic boundaries may apply: the Dirichlet boundaries Γpf, the Neumann boundaries Γvf, and the Robin boundaries (where a normal acoustic impedance is specified), Γzf. nf and τs are the unit vectors normal to the fluid and structural surfaces, respectively.Figure 1
A principal sketch of the structural-acoustic system and the corresponding multidomain model. Variables are velocity potentialϕfs and structure displacement wsie.To achieve a more effective computational procedure and derive solutions describing vibroacoustic behaviors at relatively high frequencies, a multidomain spectral method is developed. Firstly, according to the distribution of the interface Γfs and Γf, the whole system is partitioned into Ns nonoverlapping substructures and Nf (Nf ≥ Ns) acoustic subdomains along the x-axis as shown in Figure 1. To generalize the model, the symbol e is used to mark both the substructures Ωse and acoustic subdomains Ωfe, and for the components without fluid-structure interaction Ωse = 0. Moreover, the responses of any subdomain in the global coordinate system (o-xyz) can be independently described by the associated local coordinate system (oe-xe·ye·ze), where o and oe are the origins. Secondly, the subdomains are synthesized using the identified Lagrange multipliers (λf for the fluid domains and λs for the structural domain) by enforcing a weak continuity of the solution at the connecting interfaces, while the fluid and the structure components are coupled by the distributed internal forces. Finally, the nonconditional variational principle [28] and the spectral discretization procedure [19, 20] are applied to each computational domain to derive the governing equation of the entire vibroacoustic system. Contrary to FEM, which partitions the problem domain into a large amount of small elements, the proposed approach divides it into several large subdomains.


## 3. Theoretical Formulations

### 3.1. Structural Vibration Model

Consider first the structure domain,Ωs, where the interface boundary between the substructures Ωse and Ωse+1 is denoted as Γse,e+1. The total boundary of the e-th substructure can be given as ∑se=Γfse∪Γuse∪Γσse∪Γse,e+1∪Γse−1,e, where Γfse, Γuse, and Γσse denote, respectively, the fluid-structure interface, prescribed displacement boundary, and surface force boundary with regard to the e-th substructure. The Lagrangian functional for the structure domain Ωs can be written using Cartesian tensor notation as(1)Js=∫t0t1∑e=1Ns∫Ωseρs2W˙sieW˙sie−12σijeεijedΩse+∫ΓseFsieWsiedΓsedt+∫t0t1∑e=1Ns−1Wλse,e+1dt−∫t0t1∑e=1NsWsfedt,where the first term gives the Lagrangian for the uncoupled substructure, in which repeated indices i and j (i, j = 1, 2, 3) imply the Einstein summation convention, Wsie is the i-th component of the displacement, ρs is the mass density of the structure, σije and εije denote the stress and strain tensor, respectively, and Fsie are prescribed surface forces and the superposed dot is the time derivative. The interface potential Wλse,e+1 is incorporated into Js to relax the constraint conditions between adjacent substructures.The work done by the acoustic back pressure excitationpfe on the fluid-structure interface Γfse is represented by Wsfe, which is expressed as [1](2)Wsfe=∫ΓfsewsiepfenfiedΓfse=∫Γfse−ρfwsieϕ˙fenfiedΓfse,where nfie is the unit normal external to Ωfs and ϕfe is the velocity potential.Performing variational calculations of equation (1) and applying the stationary condition of Js, yields δJs=∫t0t1δHsdt=0, in which δHs is derived as(3)δHs=∑e=1Ns∫Ωseρsw˙sieδw˙sie−σijeδεijedΩse+∑e=1Nsδwie∫ΓsieFsiedΓse+∫Γfsieρfϕ˙fenfedΓfse+∑e=1Ns−1∫Γse,e+1δWλse,e+1dΓse,e+1.For a general elastic structure, the strain-displacement relations and strain-stress relations can be given according to the linearized elasticity theory as [29](4)εije=12wsi,je+wsj,ie,σije=aijkleεkle,where wsi,je=∂wsie/∂χj and aijkle are the elastic constants.The interface potentialWλse,e+1 is expressed by defining a vector of Lagrange multipliers λsie,e+1 on the boundary interface Γse,e+1 of connecting substructures:(5)Wλse,e+1=∫Γse,e+1λsie,e+1wsie−wsie+1dΓse,e+1.Integrating equation (3) by parts, and employing Green’s formula, δHs can be rewritten with respect to small perturbations δλsie,e+1 and δwsie as follows:(6)δHs=∑e=1Nsδwsie∫Ωse−ρsw¨sie+σij,je︸1dΩse+∫ΓσseFsie−σijeτje︸2dΓσse+∑e=1Ns−1δwsie∫Γfsie−σijeτje+ρfϕ˙fenie︸3dΓfse+∑e=1Ns−1∫Γse,e+1δλsie,e+1wsie−wsie+1︸4dΓse,e+1+∑e=1Ns−1∫Γse,e+1δwsieλsie,e+1−σijeτje−δwsie+1σije+1τje+1+λsie,e+1︸5dΓse,e+1,where τje is the outward drawn normal to the substructure surface ∑se. Since δλsie,e+1 and δwsie are assumed to be arbitrary variations in Ωse, the following relations can be easily derived from equation (6): (1) the governing equations, (2) the prescribed force boundary conditions on Γσse, (3) the continuity conditions on the fluid-structure interface Γfse, (4) the continuity conditions on the substructures interfaces Γse,e+1, and (5) Lagrange multipliers derived in terms of original variables, namely,(7)λsie,e+1=σijeτjeorλsie,e+1=−σije+1τje+1.The Lagrange multipliersλsie,e+1 physically represent the additional constraint quantity needed to maintain the motion inside the interface, and here it is determined to be the intrinsic interface stress of the structure.It should be noted that, if thin plates or shell assemblies are used to enclose the acoustic cavity, the variational equation (6) and interface potential Wλse,e+1 cannot be directly employed due to assumptions concerning the stress state through the thickness in thin plate and shell theories. A slight modification is necessary to allow for the coupling effects between the in-plane and out-of-plane vibrations occur at the interfaces. Using the plate assemblies as an example, the solid stress and strain σije and εije in equations (1) and (3) should be replaced by stress resultants τ¯e and generalized strains ε¯e for the plate, and(8)ε¯eεseκseT,τ¯e=N¯seM¯seT,where the superscript T denotes the transpose, εse=ɛxeɛyeɛxye is the normal and shear strains of the neutral surface, κseκxeκyeκxye is curvature changes and twist of the neutral surface. N¯se=N¯xeN¯yeN¯xye and M¯se=M¯xeM¯yeM¯xye denote, respectively, the force and moment resultant vectors.Following the thin plate theories [30], the strain-stress relations and strain-displacement relations can be given as(9)τ¯e=Eρε¯e=Kεe00Kκeε¯e,ε¯e=Dρdpe,where Kεe and Kκe are the extensional and bending stiffness, Dρ is a differential operator derived from plate theory, and dpe=ws1ews2ews3e is midsurface displacements of the plate in the three local coordinate directions.Moreover, the interface potentialWλse,e+1 should be rewritten according to thin plate theories as(10)Wλse,e+1=∫Γse,e+1λse,e+1TTeCdpe−Te+1Cdpe+1dΓse,e+1.Here, the matrix expression within the bracket is the displacement constraint equations on the common interface, andλse,e+1=λsue,e+1λsve,e+1λswe,e+1λsϕe,e+1T is the vector of Lagrange multipliers defined on the interface between the e-th and (e + 1)-th plates. The differential operator C and the transformation matrix Te of the e-th plate which maps dpe into the global coordinate can be given as(11)Te=cosαe0−sinαe00100sinαe0cosαe00001,C=10001000100∂∂x,where αe is the angle between the xe-axis of the e-th plate and x-axis of the global coordinate system. Submitting equations (8)–(10) into equation (1) and performing the variational calculation in the same way as elastic solids, the Lagrange multipliers can be determined as(12)λsue,e+1=N¯xecosαe−Vxesinαe,λsve,e+1=N¯ye,λsϕe,e+1=−M¯xe,λswe,e+1=N¯xesinαe+Vxecosαe,where Vxe is the Kevlin–Kirchhoff edge reaction (see Ref. [30]). A similar procedure can work for the case of thin shell assemblies without any difficulties. The identified Lagrange multipliers can then be substituted back into equation (3). Compared to the common Lagrange multiplier technique, the present procedure yields the nonconditional variation principle in a weak form in terms of only the original primary variables without adding new variables. It should be emphasized that the present formulation is general; it can not only be directly combined with different simplified theory to analyze dynamics of structures, such as plates and shells, but also can be applied to structures that consist of an assembly of components or thick solids that require using general elasticity theory.


### 3.2. Acoustic Vibration Model

Consider next the fluid domain Ωf, where the total acoustic boundary ∑fe for the e-th subdomain is given as ∑feΓpfe∪Γvfe∪Γzfe∪Γfe,e+1∪Γfe−1,e∪Γfse (for the subdomains without fluid-structure interaction, Γfse=0), in which Γpfe, Γvfe, and Γvfe represent the Dirichlet boundaries, the Neumann boundaries, and the Robin boundaries related to the e-th subdomain, respectively.Dealing with the fluid in terms of the velocity potentialϕf, which assures that the flow is irrotational, the corresponding kinetic and potential energies for small-amplitude oscillations in e-th fluid subdomain can be given as [31](13)Tfe=ρf2∫Ωfeu˙fe2dΩfe=ρf2∫Ωfegradϕfe2dΩfe,(14)Ufe=1ρfc02∫Ωfepfe2dΩfe=ρf2c02∫Ωfeϕ˙fe2dΩfe,where ufe, ϕfe, and pfe denote, respectively, the particle displacement, velocity potential, and sound pressure. The constant co is the speed of sound in the fluid, ρf is the mean mass density of the fluid, and grad is the notation for the gradient.The work done over the boundary surfaces is equal to the acoustic pressure multiplied by the particle acoustic displacement. Considering only harmonic motion in this work, the corresponding energy expressions (Wzfe, Wvfe, and Wsfe) with regard to the impedance boundary, the Neumann boundary and the fluid-structure interface can then be, respectively, expressed as follows:(15)Wzfe=∫ΓZfepfeunfedΓzfe=−ρf2jωZfe∫ΓZfeϕ˙fe2dΓzfe,(16)Wvfe=−∫Γvfepfev¯nfejωdΓvfe=ρfjω∫Γvfeϕ˙fkev¯nfedΓvfe,(17)Wsfe=−∫Γsfepfewse⋅nfedΓsfe=ρf∫Γsfeϕ˙fewse⋅nfedΓsfe,where ω is the angular frequency, n in the subscripts denotes the outward normal to the fluid boundaries, and j=−1; wse is the displacement vector of the structure and nfe is unit outward normal. Zfe and v¯fe are the specific acoustic impedance and the externally imposed particle velocity, respectively.The adjacent subdomains can be synthesized by introducing a Lagrange multiplierλfe,e+1 defined along the interface λfe,e+1. The additional interface energy can be consequently obtained by multiplying the constraint relationship by the undetermined Lagrange multiplier as(18)Wλfe,e+1=∫Γλfeλfe,e+1ρfϕfe−ϕfe+1dΓλfe,Γλfe=Γpfe∪Γfe,e+1.The Dirichlet boundary herein is treated as a special interface as those between adjacent segments by replacingϕfe with −pnf−e/jρfω, in which pnf−e is the prescribed acoustic pressure. Moreover, the work done by the harmonic volume distributed acoustic source can be given as(19)WFfe=−∫ΩλfepfeQfejωdΩfe=pfjω∫Ωfeϕ˙feQfedΩfe,where Qfe is the distributed strength of the volume velocity source. The corresponding variational functional can then be written as(20)Jf=∫t0t1∑e=1NfTfe−Ufe+Wzfe+Wvfe+WFfedt+∫t0t1∑e=1NfWsfedt+∫t0t1∑e=1Nf−1Wλfe.e+1dt.Invoking the stationarity ofJf with respect to independent variables λfe,e+1 and ϕfe gives the principle of virtual potentials for the fluid. Therefore, the first variation of Jf with respect to small perturbation λfe,e+1 and δϕfe should be equal to zero. Substituting equations (13)–(19) into equation (20), performing variational calculations, and dividing both sides by ρf readily yield the variational equation as ∫t0t1Δdt=0, where(21)Δ=∑e=1Nfδϕfe∫Ωfe−divgradϕfe+1c02ϕ¨fe−Qfe︸1+∫Γsfe∂ϕfe∂nfe−w˙se⋅nse︸2dΓsfe+δϕfe∑e=1Nf∫Γsfe∂ϕfe∂nfe+ρfzfeϕ˙fe︸3dΓzfe+∫Γvfe∂ϕfe∂nfe−v¯nfe︸4dΓzfe+∑e=1Nf−1∫Γsfe∂λfe,e+1+ϕfe−ϕfe+1︸5+δϕfe∂ϕfe∂nfe+λfe,e+1+δϕfe∂ϕfe+1∂nfe+1−λfe,e+1︸6dΓλfe,and nfe denotes the unit normal external to the fluid boundary surface ∑fe. Since δλfe,e+1 and δϕfe, are arbitrary, expressions 1–5 in equation (21) under the integral intervals must be zero; this naturally leads to the following independent relations: (1) the Helmholtz equation, (2) fluid-structural conditions on Γsfke, (3) impedance boundary conditions on Γsfke, (4) Neumann boundary conditions on Γvfke, (5) continuity condition on common interfaces Γfke as well as the Dirichlet boundary Γpfke, and (6) the undetermined Lagrange multipliers identified in terms of original variables ϕfe, namely,(22)λke,e+1=−∂ϕfe∂nfe,orλke,e+1=−∂ϕfe+1∂nfe,e+1.Analogously, the constraint relations between the adjacent subdomains can be removed by substituting equation (22) into equation (20), which may result in the nonconditional variational principle by eliminating the appearance of the Lagrange multiplier λfe,e+1 and giving a weak form expressed entirely in terms of ϕfε, which is not shown here for brevity. Here, λfe,e+1 is physically identified as the interface normal velocity.


### 3.3. Spectral Discretization of the Interaction Model

To derive the discretized equations of motion for vibroacoustic coupled systems, the spectral method is employed to discretize the energy functionals (1) and (20). The i-th (i = x, y, z) displacement wie for the e-th substructure and the velocity potential ϕfe for the e-th subdomain can be, respectively, approximated in terms of truncated boundary characteristic polynomial series [14, 20]:(23)wiexe,ye,t=fi,we∑k=1Nx∑l=1NyPkxePyedi,k,let=Wiexe,yedsieejωt,(24)ϕexe,ye,ze,t=fϕe∑r=1Mx∑s=1My∑p=1MzPrxePsyePpzeqrspet=Φexe,ye,zeqfeejωt,where Nx, Ny, Mx, My, and Mz are the truncation orders of the polynomial series, Pαβ1,2,3⋯;β=xe,ye,ze are the α order polynomials, di,klet and qrspet are generalized coordinates for the e-th substructure and fluid subdomain, Wie and Φie represent the admissible function vectors, while qfe and dsie are generalized coordinate vectors. The functions fi,we=f¯i,wxexef¯i,wyeye and fϕxe=f¯ϕxexef¯ϕyeyef¯ϕzeze are boundary functions enforcing boundary conditions, where f¯i,wxe and f¯ϕxeδ=x,y,z are components in different coordinate directions and one can refer to [32] for detailed examples. It should also be noted that only the 2D structures and 3D fluid are used here to illustrate the spectral approximation process for example. However, polynomial expansions for problems in other dimensions (such as 2-D fluid, 2-D or 3-D elastic solid) can be analogously performed using the present formulation without any difficulty.Substituting equations (7) and (22) as well as matrix expressions of equations (23) and (24) into equations (1) and (20) and performing variational operation with regard to qfe and dsie readily result in the following assembled system of equations:(25)M˜q˜¨+C˜q˜˙+K˜q˜=F˜,where M˜=Mf00Ms, K˜=Kf+KLf00Ks+KLs, C˜=ZfCfs−Csf0, and q˜=q¯fTd¯sTT, F˜=FfTFsTT. The subscripts f and s will be generated by relevant submatrices with regard to the fluid and structure, respectively. The matrices Mδ and kδδ=f,s are equivalent mass and stiffness matrices. The global acoustic-structural coupling matrix is denoted as Cfs, and Zf is the matrix related to the impedance boundary of the fluid. The matrices KLf and KLS are the interface stiffness matrices introduced by the generalized variational principle, Fδ (δ=f,s) are generalized force vectors, and q¯f and d¯s are global generalized coordinate vectors for the fluid and the structure, respectively, which can be expressed as(26)d¯s=dsx1Tdsy1Tdsz1T︸1−stsubstructure⋯dsxNsTdsyNsTdszNsT︸NS−thsubstructureT,q¯f=qf1T⋯qfeT⋯qfNfTT.Submatrices of the system matrices are given by integrals of various combinations of admissible functions and their spatial derivatives, as shown in Appendix A for completeness.The steady-state response ofq can be readily given by assuming the harmonic excitation as(27)q˜=−ω2M˜+jωC˜+K˜−1F˜.The substitution of equation (27) into equations (23) and (24) leads to steady-state responses of the structural displacements and the fluid velocity potential. Furthermore, the acoustic pressure for the e-th fluid subdomain can be given as(28)piexe,ye,t=jωρfΦxe,ye,Zeqfeejωt.Obviously, using the state vectorX˜T=q˜Tq˜˙T, equation (25) can also be easily modified to the generalized eigenvalue problem of the form AX˜=λ˜X˜ by setting F˜ to be zero, in which the matrix A is given by(29)A=0I−M−1K˜−M−1C˜,and the i-th eigenvalue λ¯i can be written in the following form:(30)λ¯i,λ¯i∗=−ςiωi±jωi1−ςi2,where the notation ∗ denotes complex conjugate, ζi and ωi are the i-th modal damping ratio and frequency, which can be expressed as(31)ωi=Reλ˜i2+Imλ˜2,ζi=iReλ˜i−ωi.


## 3.1. Structural Vibration Model

Consider first the structure domain,Ωs, where the interface boundary between the substructures Ωse and Ωse+1 is denoted as Γse,e+1. The total boundary of the e-th substructure can be given as ∑se=Γfse∪Γuse∪Γσse∪Γse,e+1∪Γse−1,e, where Γfse, Γuse, and Γσse denote, respectively, the fluid-structure interface, prescribed displacement boundary, and surface force boundary with regard to the e-th substructure. The Lagrangian functional for the structure domain Ωs can be written using Cartesian tensor notation as(1)Js=∫t0t1∑e=1Ns∫Ωseρs2W˙sieW˙sie−12σijeεijedΩse+∫ΓseFsieWsiedΓsedt+∫t0t1∑e=1Ns−1Wλse,e+1dt−∫t0t1∑e=1NsWsfedt,where the first term gives the Lagrangian for the uncoupled substructure, in which repeated indices i and j (i, j = 1, 2, 3) imply the Einstein summation convention, Wsie is the i-th component of the displacement, ρs is the mass density of the structure, σije and εije denote the stress and strain tensor, respectively, and Fsie are prescribed surface forces and the superposed dot is the time derivative. The interface potential Wλse,e+1 is incorporated into Js to relax the constraint conditions between adjacent substructures.The work done by the acoustic back pressure excitationpfe on the fluid-structure interface Γfse is represented by Wsfe, which is expressed as [1](2)Wsfe=∫ΓfsewsiepfenfiedΓfse=∫Γfse−ρfwsieϕ˙fenfiedΓfse,where nfie is the unit normal external to Ωfs and ϕfe is the velocity potential.Performing variational calculations of equation (1) and applying the stationary condition of Js, yields δJs=∫t0t1δHsdt=0, in which δHs is derived as(3)δHs=∑e=1Ns∫Ωseρsw˙sieδw˙sie−σijeδεijedΩse+∑e=1Nsδwie∫ΓsieFsiedΓse+∫Γfsieρfϕ˙fenfedΓfse+∑e=1Ns−1∫Γse,e+1δWλse,e+1dΓse,e+1.For a general elastic structure, the strain-displacement relations and strain-stress relations can be given according to the linearized elasticity theory as [29](4)εije=12wsi,je+wsj,ie,σije=aijkleεkle,where wsi,je=∂wsie/∂χj and aijkle are the elastic constants.The interface potentialWλse,e+1 is expressed by defining a vector of Lagrange multipliers λsie,e+1 on the boundary interface Γse,e+1 of connecting substructures:(5)Wλse,e+1=∫Γse,e+1λsie,e+1wsie−wsie+1dΓse,e+1.Integrating equation (3) by parts, and employing Green’s formula, δHs can be rewritten with respect to small perturbations δλsie,e+1 and δwsie as follows:(6)δHs=∑e=1Nsδwsie∫Ωse−ρsw¨sie+σij,je︸1dΩse+∫ΓσseFsie−σijeτje︸2dΓσse+∑e=1Ns−1δwsie∫Γfsie−σijeτje+ρfϕ˙fenie︸3dΓfse+∑e=1Ns−1∫Γse,e+1δλsie,e+1wsie−wsie+1︸4dΓse,e+1+∑e=1Ns−1∫Γse,e+1δwsieλsie,e+1−σijeτje−δwsie+1σije+1τje+1+λsie,e+1︸5dΓse,e+1,where τje is the outward drawn normal to the substructure surface ∑se. Since δλsie,e+1 and δwsie are assumed to be arbitrary variations in Ωse, the following relations can be easily derived from equation (6): (1) the governing equations, (2) the prescribed force boundary conditions on Γσse, (3) the continuity conditions on the fluid-structure interface Γfse, (4) the continuity conditions on the substructures interfaces Γse,e+1, and (5) Lagrange multipliers derived in terms of original variables, namely,(7)λsie,e+1=σijeτjeorλsie,e+1=−σije+1τje+1.The Lagrange multipliersλsie,e+1 physically represent the additional constraint quantity needed to maintain the motion inside the interface, and here it is determined to be the intrinsic interface stress of the structure.It should be noted that, if thin plates or shell assemblies are used to enclose the acoustic cavity, the variational equation (6) and interface potential Wλse,e+1 cannot be directly employed due to assumptions concerning the stress state through the thickness in thin plate and shell theories. A slight modification is necessary to allow for the coupling effects between the in-plane and out-of-plane vibrations occur at the interfaces. Using the plate assemblies as an example, the solid stress and strain σije and εije in equations (1) and (3) should be replaced by stress resultants τ¯e and generalized strains ε¯e for the plate, and(8)ε¯eεseκseT,τ¯e=N¯seM¯seT,where the superscript T denotes the transpose, εse=ɛxeɛyeɛxye is the normal and shear strains of the neutral surface, κseκxeκyeκxye is curvature changes and twist of the neutral surface. N¯se=N¯xeN¯yeN¯xye and M¯se=M¯xeM¯yeM¯xye denote, respectively, the force and moment resultant vectors.Following the thin plate theories [30], the strain-stress relations and strain-displacement relations can be given as(9)τ¯e=Eρε¯e=Kεe00Kκeε¯e,ε¯e=Dρdpe,where Kεe and Kκe are the extensional and bending stiffness, Dρ is a differential operator derived from plate theory, and dpe=ws1ews2ews3e is midsurface displacements of the plate in the three local coordinate directions.Moreover, the interface potentialWλse,e+1 should be rewritten according to thin plate theories as(10)Wλse,e+1=∫Γse,e+1λse,e+1TTeCdpe−Te+1Cdpe+1dΓse,e+1.Here, the matrix expression within the bracket is the displacement constraint equations on the common interface, andλse,e+1=λsue,e+1λsve,e+1λswe,e+1λsϕe,e+1T is the vector of Lagrange multipliers defined on the interface between the e-th and (e + 1)-th plates. The differential operator C and the transformation matrix Te of the e-th plate which maps dpe into the global coordinate can be given as(11)Te=cosαe0−sinαe00100sinαe0cosαe00001,C=10001000100∂∂x,where αe is the angle between the xe-axis of the e-th plate and x-axis of the global coordinate system. Submitting equations (8)–(10) into equation (1) and performing the variational calculation in the same way as elastic solids, the Lagrange multipliers can be determined as(12)λsue,e+1=N¯xecosαe−Vxesinαe,λsve,e+1=N¯ye,λsϕe,e+1=−M¯xe,λswe,e+1=N¯xesinαe+Vxecosαe,where Vxe is the Kevlin–Kirchhoff edge reaction (see Ref. [30]). A similar procedure can work for the case of thin shell assemblies without any difficulties. The identified Lagrange multipliers can then be substituted back into equation (3). Compared to the common Lagrange multiplier technique, the present procedure yields the nonconditional variation principle in a weak form in terms of only the original primary variables without adding new variables. It should be emphasized that the present formulation is general; it can not only be directly combined with different simplified theory to analyze dynamics of structures, such as plates and shells, but also can be applied to structures that consist of an assembly of components or thick solids that require using general elasticity theory.


## 3.2. Acoustic Vibration Model

Consider next the fluid domain Ωf, where the total acoustic boundary ∑fe for the e-th subdomain is given as ∑feΓpfe∪Γvfe∪Γzfe∪Γfe,e+1∪Γfe−1,e∪Γfse (for the subdomains without fluid-structure interaction, Γfse=0), in which Γpfe, Γvfe, and Γvfe represent the Dirichlet boundaries, the Neumann boundaries, and the Robin boundaries related to the e-th subdomain, respectively.Dealing with the fluid in terms of the velocity potentialϕf, which assures that the flow is irrotational, the corresponding kinetic and potential energies for small-amplitude oscillations in e-th fluid subdomain can be given as [31](13)Tfe=ρf2∫Ωfeu˙fe2dΩfe=ρf2∫Ωfegradϕfe2dΩfe,(14)Ufe=1ρfc02∫Ωfepfe2dΩfe=ρf2c02∫Ωfeϕ˙fe2dΩfe,where ufe, ϕfe, and pfe denote, respectively, the particle displacement, velocity potential, and sound pressure. The constant co is the speed of sound in the fluid, ρf is the mean mass density of the fluid, and grad is the notation for the gradient.The work done over the boundary surfaces is equal to the acoustic pressure multiplied by the particle acoustic displacement. Considering only harmonic motion in this work, the corresponding energy expressions (Wzfe, Wvfe, and Wsfe) with regard to the impedance boundary, the Neumann boundary and the fluid-structure interface can then be, respectively, expressed as follows:(15)Wzfe=∫ΓZfepfeunfedΓzfe=−ρf2jωZfe∫ΓZfeϕ˙fe2dΓzfe,(16)Wvfe=−∫Γvfepfev¯nfejωdΓvfe=ρfjω∫Γvfeϕ˙fkev¯nfedΓvfe,(17)Wsfe=−∫Γsfepfewse⋅nfedΓsfe=ρf∫Γsfeϕ˙fewse⋅nfedΓsfe,where ω is the angular frequency, n in the subscripts denotes the outward normal to the fluid boundaries, and j=−1; wse is the displacement vector of the structure and nfe is unit outward normal. Zfe and v¯fe are the specific acoustic impedance and the externally imposed particle velocity, respectively.The adjacent subdomains can be synthesized by introducing a Lagrange multiplierλfe,e+1 defined along the interface λfe,e+1. The additional interface energy can be consequently obtained by multiplying the constraint relationship by the undetermined Lagrange multiplier as(18)Wλfe,e+1=∫Γλfeλfe,e+1ρfϕfe−ϕfe+1dΓλfe,Γλfe=Γpfe∪Γfe,e+1.The Dirichlet boundary herein is treated as a special interface as those between adjacent segments by replacingϕfe with −pnf−e/jρfω, in which pnf−e is the prescribed acoustic pressure. Moreover, the work done by the harmonic volume distributed acoustic source can be given as(19)WFfe=−∫ΩλfepfeQfejωdΩfe=pfjω∫Ωfeϕ˙feQfedΩfe,where Qfe is the distributed strength of the volume velocity source. The corresponding variational functional can then be written as(20)Jf=∫t0t1∑e=1NfTfe−Ufe+Wzfe+Wvfe+WFfedt+∫t0t1∑e=1NfWsfedt+∫t0t1∑e=1Nf−1Wλfe.e+1dt.Invoking the stationarity ofJf with respect to independent variables λfe,e+1 and ϕfe gives the principle of virtual potentials for the fluid. Therefore, the first variation of Jf with respect to small perturbation λfe,e+1 and δϕfe should be equal to zero. Substituting equations (13)–(19) into equation (20), performing variational calculations, and dividing both sides by ρf readily yield the variational equation as ∫t0t1Δdt=0, where(21)Δ=∑e=1Nfδϕfe∫Ωfe−divgradϕfe+1c02ϕ¨fe−Qfe︸1+∫Γsfe∂ϕfe∂nfe−w˙se⋅nse︸2dΓsfe+δϕfe∑e=1Nf∫Γsfe∂ϕfe∂nfe+ρfzfeϕ˙fe︸3dΓzfe+∫Γvfe∂ϕfe∂nfe−v¯nfe︸4dΓzfe+∑e=1Nf−1∫Γsfe∂λfe,e+1+ϕfe−ϕfe+1︸5+δϕfe∂ϕfe∂nfe+λfe,e+1+δϕfe∂ϕfe+1∂nfe+1−λfe,e+1︸6dΓλfe,and nfe denotes the unit normal external to the fluid boundary surface ∑fe. Since δλfe,e+1 and δϕfe, are arbitrary, expressions 1–5 in equation (21) under the integral intervals must be zero; this naturally leads to the following independent relations: (1) the Helmholtz equation, (2) fluid-structural conditions on Γsfke, (3) impedance boundary conditions on Γsfke, (4) Neumann boundary conditions on Γvfke, (5) continuity condition on common interfaces Γfke as well as the Dirichlet boundary Γpfke, and (6) the undetermined Lagrange multipliers identified in terms of original variables ϕfe, namely,(22)λke,e+1=−∂ϕfe∂nfe,orλke,e+1=−∂ϕfe+1∂nfe,e+1.Analogously, the constraint relations between the adjacent subdomains can be removed by substituting equation (22) into equation (20), which may result in the nonconditional variational principle by eliminating the appearance of the Lagrange multiplier λfe,e+1 and giving a weak form expressed entirely in terms of ϕfε, which is not shown here for brevity. Here, λfe,e+1 is physically identified as the interface normal velocity.


## 3.3. Spectral Discretization of the Interaction Model

To derive the discretized equations of motion for vibroacoustic coupled systems, the spectral method is employed to discretize the energy functionals (1) and (20). The i-th (i = x, y, z) displacement wie for the e-th substructure and the velocity potential ϕfe for the e-th subdomain can be, respectively, approximated in terms of truncated boundary characteristic polynomial series [14, 20]:(23)wiexe,ye,t=fi,we∑k=1Nx∑l=1NyPkxePyedi,k,let=Wiexe,yedsieejωt,(24)ϕexe,ye,ze,t=fϕe∑r=1Mx∑s=1My∑p=1MzPrxePsyePpzeqrspet=Φexe,ye,zeqfeejωt,where Nx, Ny, Mx, My, and Mz are the truncation orders of the polynomial series, Pαβ1,2,3⋯;β=xe,ye,ze are the α order polynomials, di,klet and qrspet are generalized coordinates for the e-th substructure and fluid subdomain, Wie and Φie represent the admissible function vectors, while qfe and dsie are generalized coordinate vectors. The functions fi,we=f¯i,wxexef¯i,wyeye and fϕxe=f¯ϕxexef¯ϕyeyef¯ϕzeze are boundary functions enforcing boundary conditions, where f¯i,wxe and f¯ϕxeδ=x,y,z are components in different coordinate directions and one can refer to [32] for detailed examples. It should also be noted that only the 2D structures and 3D fluid are used here to illustrate the spectral approximation process for example. However, polynomial expansions for problems in other dimensions (such as 2-D fluid, 2-D or 3-D elastic solid) can be analogously performed using the present formulation without any difficulty.Substituting equations (7) and (22) as well as matrix expressions of equations (23) and (24) into equations (1) and (20) and performing variational operation with regard to qfe and dsie readily result in the following assembled system of equations:(25)M˜q˜¨+C˜q˜˙+K˜q˜=F˜,where M˜=Mf00Ms, K˜=Kf+KLf00Ks+KLs, C˜=ZfCfs−Csf0, and q˜=q¯fTd¯sTT, F˜=FfTFsTT. The subscripts f and s will be generated by relevant submatrices with regard to the fluid and structure, respectively. The matrices Mδ and kδδ=f,s are equivalent mass and stiffness matrices. The global acoustic-structural coupling matrix is denoted as Cfs, and Zf is the matrix related to the impedance boundary of the fluid. The matrices KLf and KLS are the interface stiffness matrices introduced by the generalized variational principle, Fδ (δ=f,s) are generalized force vectors, and q¯f and d¯s are global generalized coordinate vectors for the fluid and the structure, respectively, which can be expressed as(26)d¯s=dsx1Tdsy1Tdsz1T︸1−stsubstructure⋯dsxNsTdsyNsTdszNsT︸NS−thsubstructureT,q¯f=qf1T⋯qfeT⋯qfNfTT.Submatrices of the system matrices are given by integrals of various combinations of admissible functions and their spatial derivatives, as shown in Appendix A for completeness.The steady-state response ofq can be readily given by assuming the harmonic excitation as(27)q˜=−ω2M˜+jωC˜+K˜−1F˜.The substitution of equation (27) into equations (23) and (24) leads to steady-state responses of the structural displacements and the fluid velocity potential. Furthermore, the acoustic pressure for the e-th fluid subdomain can be given as(28)piexe,ye,t=jωρfΦxe,ye,Zeqfeejωt.Obviously, using the state vectorX˜T=q˜Tq˜˙T, equation (25) can also be easily modified to the generalized eigenvalue problem of the form AX˜=λ˜X˜ by setting F˜ to be zero, in which the matrix A is given by(29)A=0I−M−1K˜−M−1C˜,and the i-th eigenvalue λ¯i can be written in the following form:(30)λ¯i,λ¯i∗=−ςiωi±jωi1−ςi2,where the notation ∗ denotes complex conjugate, ζi and ωi are the i-th modal damping ratio and frequency, which can be expressed as(31)ωi=Reλ˜i2+Imλ˜2,ζi=iReλ˜i−ωi.


## 4. Theoretical Formulations

In this section, two practical examples are presented to illustrate the applicability of the method and its advantages for describing the vibroacoustic interaction between elastic structures and the internal acoustic cavities. The Gauss elimination algorithm is utilized to get the numerical solutions of equation (25), while the QZ algorithm is implemented to obtain eigenvalues and eigenvectors from equation (29). Corresponding calculation procedures are worked out in Matlab 7.0 and run on an Intel Xeon X5570 2.93 GHz computer system. The results are validated by comparison with those from the commercial software ANSYS 17.1. In order to compare the calculation efficiency, the computational time for both methods is recorded. The time for the present approach includes efforts for the construction of system matrices and the solution of equation (25) as well as the postprocessing of response data, while only the time required for the sparse solution in the harmonic analysis is taken into account for the FEM. Although the proposed formulation mathematically permits the use of any spectral orthogonal basis functions as appropriate admissible functions, the first kind of Chebyshev orthogonal polynomials, which have been rationally verified to have rapid convergence and superior stability in previous numerical computations for vibration and acoustic problems [33], are adopted in this work.


### 4.1. A Simplified Vehicle Cabin Backed by Plate Assembly

To ensure the applicability and computational advantages of the present formulation for analyzing vibroacoustic problems with relatively complex geometries, a simplified vehicle cabin in Figure2 is employed; it consists of a built-up cavity backed by a plate assembly. The structure is assembled by three plate components oriented arbitrarily. Unlike the pure plate bending problems as in [15, 16], the coupling effects of the out-of-plane modes and in-plane modes of the plates should be considered here, especially for vibrations in higher frequencies [34]. Moreover, besides the transformation of the physical coordinate to the dimensionless local coordinate, the mapping of the irregular geometries into the regular ones should also be included:(32)xe=Lxe2ξe+Lxe2,ye=Lye2ηe+Lye2,ze=Lze2ζe+Lze2,ye=Lxetanαi41+ξe+1+ζe,where Lxe, Lye, and Lze are dimensions of subdomains in three directions (Lze is the dimension of the left surface of each acoustic subdomain in z). xe, ye, and ze are the local physical coordinates for each fluid subdomain. ξe, ηe, and ζe are the dimensionless local coordinates defined in the interval [−1, 1]. αi (i = 1, 2) is the leaning angle as in Figure 2. Following the standard isoparametric mapping procedure as used in FEM [7], the Jacobian matrix for the transformation can be readily achieved, and thus the functional in equations (1) and (20) can be mapped from the physical domain to computational domain. Substituting equation (32) into equation (25) finally yields the system of equations in the local coordinates.Figure 2
Interior vibroacoustic interaction of a simplified vehicle cabin. Point coordinates: A(0.8, 0.6, 1), B(0.5, 0.45, 0.5), C(1.0, 0.2, 1), and D(0.75, 0.45, 0.5).The uncoupled edges of all plates are taken to be clamped, while all uncoupled walls of the cavity are assumed to be rigid. Taking benefit of the generality of the present formulation, in addition to representing the out-of-plane motion of the plate by Kirchhoff’s plate bending theory, the theory of elasticity for two-dimensional solids [29] is also employed to formulate the in-plane motion of the plates. The interaction between the in-plane and out-of-plane vibrations occurs only at the coupling interface. The physical parameters for the system are the same as those used in Ref. [35], the plates are 1.5 mm thick and made of steel with Young’s modulus E = 210 GPa, the mass density ρs=7850kg/m3 and the Poisson ratio vs = 0.3, while the cavity is filled with air of fluid density ρf=1.225kg/m3 and the speed of sound cf = 340 m/s. During the calculation, the plate and the cavity are divided into three components in the x direction according to the plate junction positions. Both the in-plane and out-of-plane displacements should be approximated in terms of Chebyshev polynomials, and thus the total number of degrees of freedom (DoFs) is given by Nsum=3Nx⋅Ny⋅NL+Mx⋅Mx⋅Mx⋅NL.The ability of the present method to characterize well free vibration behaviors of the coupled system is first demonstrated. The prediction accuracy of the present methodology depends mainly on the truncation order (Nx,Ny,Mx,My,Mz) of polynomial functions; thus, a convergence analysis needs to be conducted. The simple way is to recalculate the concerned quantities with an increasing number of polynomial functions until the quantities cease to change. The 15-th and the 85-th natural frequencies (f15 = 80.53 Hz and f85 = 293.11 Hz), which are, respectively, located in the low-frequency and mid-frequency bands, are arbitrarily chosen for example. The convergence curves for the present method and FEM are shown depending on the number of degrees of freedom of the problem in Figure 3. A very fine FE analysis with the element sizes he = 0.013 m and the total number of DoFs 668,920, which ensures the required accuracy for comparison purposes, serves as reference (one can find the detailed information about the FEM model in Appendix B). The formula for the relative error is used to quantify the quality of the convergence:(33)εi=ωi−ωi,refωi,ref×100%,where ωi,ref is the i-th frequency obtained from the FE model with the fine mesh. The convergence curves for the present model are obtained by increasing the truncation order (Mx,My,Mz) of Chebyshev polynomial functions for the fluid from 4 to 14 while keeping Nx=Ny=14, while those for the FEM are achieved by refining element sizes from 0.06 m for both the fluid and structure domains until the FEM predictions accuracy at a similar accuracy as the present converged results. It is shown in Figure 3 that the present model can provide a stable and rapid convergence for both the lower- and higher-order frequencies with a rather small model size as compared to the FEM, where further increments of DoFs produce relatively little change in the relative error. Moreover, a good agreement, with relative errors less than 3×10−4 between the present converged results and FE predictions, can be guaranteed with a lower number of DoFs than the reference FEM, resulting in a considerably shorter computation time. In addition, the present method is able to represent well the mode shapes. For example, comparison of the mode shapes for the 85-th mode between the present and the reference FEM solutions in Figure 4 demonstrates that the spatial distribution of both the acoustic pressure and the displacement obtained by the present method agrees well with the FEM results. To highlight the benefits of the accuracy and efficiency of the present method, several natural frequencies for different combinations of leaning angles α1 and α2 are further presented in Table 1. It shows clearly that the present results are in good agreement with those obtained from the FE model, which may be utilized as benchmark solutions for the similar numerical studies.Figure 3
Convergence curves for the present method and FEM with regard to the 15-th and 85-th natural frequencies. A FE model with element sizeshe = 0.013 m is used to provide a base of comparison (f15 = 80.53 Hz, f85 = 293.11 Hz, and NL = 3).Figure 4
The mode shapes of the structure displacements (top panels) and the acoustic pressure (bottom panels) of the cavity with regard to the 85-th natural frequencies. (a) Present solutions. (b) FEM solutions.
(a)(b)Table 1
Comparison of natural frequencies in Hz for the built-up system in Figure2 with different leaning angles α1 (Nx=Ny=14, Mx=My=Mz=8, and the leaning angle α=π/4).
α1 (rad)MethodMode no.12345203040α1 = atan (6/5)Present16.1622.4434.1535.5135.7597.80130.7166.6FEM16.1622.4434.1435.5035.7497.75130.5166.4α1 = π/3Present14.7721.3126.5631.9234.9289.98123.1152.2FEM14.7621.3126.5531.9134.9189.91123.0151.9α1 = π/4Present16.8723.4834.5337.0037.1999.82143.8174.5FEM16.8723.4734.5236.9937.1899.77143.5174.1α1 = π/6Present18.0127.5134.7338.4242.23117.6149.6189.1FEM18.0127.5134.7238.4042.21117.5149.3189.1The reliability and efficiency for numerical calculations of steady-state vibroacoustic responses are then checked, in which a unit point forceFs in the z direction is applied at point A (see Figure 2). The sound pressure responses at point B and the displacement responses at point C are collected. Comparisons of the response amplitude spectrum between the present results and the reference FEM solutions with the element sizes he = 0.013 m are plotted in Figure 5 between 1 and 400 Hz at intervals of 1 Hz. Note that both the flexural wave-number of plates and the acoustic wave-number in the fluid are considerably larger than one, and it is not a stretch to infer that both the low-frequency and mid-frequency behaviors of the coupled system can be exhibited within the considered frequency range. A good agreement between the results demonstrates the reliability of the present method over the frequency range of interest. The discrepancies in magnitude of the resonant peaks are mainly because (1) no damping is introduced in both simulations and (2) a small difference in natural frequency may cause a large difference in amplitude close to the resonant frequency. It should be mentioned that the observation point A in Figure 5(b) is located at the interface between the fluid subdomains ωf1 and ωf2; thus, the consistency between the present results and the converged FE solutions also reveals that the constraint conditions between adjacent subdomains and substructures can be incorporated correctly by using the Lagrange multiplier technology. To illustrate that the present method is capable of providing accurate spatial distribution of response fields in an explicit fashion, detailed comparisons of pressure amplitude on several lines inside the cavity between the present method and the reference FEM are plotted in Figure 6. Additionally, Figure 7 shows a color plot of the acoustic pressure and the structural deformation fields at 225 Hz. Again a good match between the results of both modeling approaches is observed. For the current calculations, the truncation orders Nx=Ny=14 of polynomial functions for each substructure and Mx=My=Mz=8 for each fluid subdomain are used; this leads to a model with a total of 3300 DoFs and the entire solution time for a single frequency of the present method is about 62 s. Moreover, following the convergence analysis in Figure 3, a FE model with 334,476 DoFs (and corresponding element sizes he = 0.016 m) is needed to gain the similar accuracy as the present method. With the same hardware, the calculation time of this FE model is approximately 1324 s, revealing an obvious computational cost saving of the proposed procedure under the similar level of accuracy.Figure 5
(a) Displacement amplitude responses at point C and (b) acoustic pressure amplitude responses of point B at the interface between the fluid subdomainsΩf1 and Ωf2Nx=Ny=14,My=Mz=8,NL=3.
(a)(b)Figure 6
Comparison of pressure amplitude on Line 1 (y = 0.3 m, z = 0.45 m) and Line 2 (y = 0.5 m, z = 0.45 m) obtained from the present method and FEM due to a point force excitation of 225 Hz Nx=Ny=14,My=Mz=8,NL=3.Figure 7
The spatial distribution of the amplitude of the structure displacements and the acoustic pressure of the cavity due to a point force excitation of 225 Hz applied on the center plate at point ANx=Ny=14,My=Mz=8,NL=3. (a) Present solutions. (b) FEM solutions.
(a)(b)The present procedure is able to introduce acoustic impedance boundaries. A simulation is subsequently followed to calculate the dynamic responses of the system with the wall impedance. Although the method can be easily utilized to take into account different impedance boundaries, only the nondissipative boundary defined by pure imaginary impedance is employed here for demonstration purposes. Considering the right surface of the cavity in Figure2 as the nondissipative boundary with uniform impedance Zf=102j, comparisons of the plate responses and sound pressure between the present and FEM solutions from 1 Hz to 400 Hz are, respectively, plotted in Figure 8 at the intervals 1 Hz. The plate response point is also selected at point C, while the pressure response point is set at point D inside the second fluid subdomain. It is obvious that both the acoustic pressure and the displacement responses obtained by the present solution well match the reference FEM results. In this case, using the FE model (he = 0.016 m) under the similar level of accuracy as the present results, required time for a sparse solution is about 1587 s per each calculation step. However, only an average of 0.32 s is used for each frequency in the harmonic analysis within the present procedure, further illustrating a significant saving of the calculation effort of the present procedure for solutions of vibroacoustic problems with acoustic impedance boundaries.Figure 8
(a) Normal displacement at (1.0 m, 0.2 m, 1 m) and (b) acoustic pressure at (0.75 m, 0.45 m, 0.5 m). Nondissipative boundary with a uniform impedanceZf=102j is applied to the right surface of the cavity Nx=Ny=14,My=Mz=8,NL=3.
(a)(b)


### 4.2. Irregular Elastic Solid with Heavy Acoustic Medium

An irregular 2D elastic structure backed by heavy acoustic medium is finally considered to further illustrate that the present formulation can expediently combine with different theories of elastic structures and easily accommodate problems in different dimensionalities. The geometric configuration as well as the dimensions of the system is shown in Figure9. The right and left edges of the solid are taken to be clamped, and all the uncoupled walls of the cavity are assumed to be rigid. The cavity is filled with water of density ρf=1000kg/m3 and acoustic speed cf = 1480 m/s. The structure is made of rubber with Young’s modulus E = 23 MPa, mass density ρs=1600kg/m3, Poisson ratio vs = 0.47, and the inclined angle α=π/12 rad. For comparison purposes, the vibroacoustic system is modeled by both the FEM and the proposed spectral method. The description of the FEM model can be found in Appendix B. For the present model, because the displacement depends on both x and y, a plate model would not be appropriate, so a 2D elastic model which is possible due to the generality of the present formulation is necessary here. And the exact small-strain 2D elasticity theory [29] is employed to model the structure. The irregular structure and the cavity are partitioned into three components NL=3 in the x direction as shown in Figure 9, and the irregular substructure can be easily mapped into the regular one similarly to the previous case. Moreover, approximating the velocity potential of the fluid and the displacements of the solid in the x and y directions in terms of spectral Chebyshev polynomials, the total number of DoFs for the present model can be given by Nsum=2Nx⋅Ny⋅NL+Mx⋅My⋅NL, where NL is the number of subdomains and Nx, Ny, Mx, and My are the truncation orders of the polynomial series.Figure 9
An irregular 2D elastic solid coupled with heavy medium.Considering modal characteristics first, convergence curves of lower-order (35-th,f35 = 490.09 Hz) and higher-order (65-th, f36 = 714.73 Hz) natural frequencies are shown in Figure 10 as examples. The global curve for the present model is performed by increasing the polynomial order (Mx=My) of the velocity potential from 5 to 20 and fixing the number of the polynomial (Nx=Ny=13) with regard to structure displacements, while the curves for the FEM are obtained by gradually refining the mesh sizes from 0.025 m until a similar accuracy as the present converged results reaches. The FE model with a very fine mesh (element sizes he = 0.0025 m and the total DoFs 175,861) is utilized in Figure 10 to provide a base of comparison, and the definition of the relative error is the same as that in equation (33). It can be seen that the accuracy of the FEM model increases steadily with an increasing number of DoFs, while the relative errors for the present model slump at the beginning and converge rapidly towards a situation where the increase of polynomial terms does not produce further improvements of prediction accuracy. An excellent agreement with relative errors on the order of 10−4 between the present converged results and FE predictions for both the lower- and higher-order frequencies can be observed, which shows that even with small models, the proposed method can provide accurate solutions within a fairly large frequency range. Note that, although the convergence analysis is only carried out for several selected frequencies, it is expected the convergence property observed is valid for other frequencies.Figure 10
Convergence curves for the present method and FEM with regard to the 35-th and the 65-th natural frequencies (f35=490.09Hz,f65=714.73Hz,NL=3).Finally, the efficiency for calculations of steady-state responses is verified, in which a unit harmonic forceFs in the y direction is applied on the top surface of the elastic structure at the global position (x,y) = (0.3 m, 0.6 m). Figures 11(a) and 11(b) present displacement responses in the y direction at (0.6 m, 0.54 m) and sound pressure responses at (0.7 m, 0.3 m) between 1 and 700 Hz at intervals of 1 Hz for the present method, respectively. The FEM reference solutions with element sizes he = 0.0025 m are also presented in Figure 11 for method verification. Good correlation illustrates that the developed method is able to accurately predict the vibroacoustic behaviors at relatively high frequencies. Additionally, comparisons of spatial distribution of the structural deformation field in the x direction between the present and FEM predictions at a higher frequency 610 Hz are further plotted in Figure 12. As can be observed, the present result corresponds well with the reference result. The contour plots of the pressure as well as displacements in the y direction are also in good agreement, which are not shown for brevity. For the present analysis, truncation orders Nx=Ny=13, Mx = 13, and My = 7 of polynomial functions for subdomains are adopted in the calculations, and the dimension of resulting system matrices is only 1287. It is much smaller than that of FEM analysis under the similar level of accuracy, where 123,292 DoFs is needed for the FE model as can be seen from Figure 10. Furthermore, the present results in Figure 11 are calculated about 15.2 times faster than the FEM sparse solution with the similar accuracy. The comparisons not only show the reliability and modeling efficiency of the present method even if both the irregular elastic solid and the heavy acoustic medium are introduced but also indicate that the proposed multidomain spectral procedure is applicable for vibroacoustic analyses in a broad frequency range.Figure 11
(a) Displacement iny direction at (0.6 m, 0.54 m) and (b) acoustic pressure at (0.7 m, 0.3 m) inside the cavity Nx=Ny=13,My=Mz=13,NL=3.
(a)(b)Figure 12
The spatial distribution of the structure displacements in thex direction due to a point force of 610 Hz applied on the structure at position (0.3 m, 0.6 m) Nx=Ny=13,My=Mz=13,NL=3.


### 4.3. Discussion

It may be concluded from the above analyses that the present method may offer a very powerful tool for providing an efficient modeling approach for vibroacoustic problems with relatively complicated geometries. By inheriting the flexibility from the multidomain synthesis as used in structural mechanics, the approach divides the full continuous system into a few large continuous subdomains instead of a large number of small elements as in FEM, which results in substantially smaller system matrices and higher computational efficiency under the similar level of prediction accuracy as compared to FEM. Moreover, the application of multidomain synthesis makes the analytical techniques possible to analyze vibroacoustic systems consisting of an assembly of structures of arbitrary directions or thick structures requiring general elasticity theory.Benefiting from the versatility of the spectral method, the present formulation mathematically permits the use of any orthogonal polynomials (Fourier, Legendre, Chebyshev, Hermite polynomials, etc.), which possess complete, linearly independent, and highly smooth properties, as admissible functions to approximate field variables. Although these polynomial functions can usually bring about rapid convergence and superior stability in numerical operations, they are required to satisfy the desired boundary conditions, such as prescribed displacement boundaries. The selection of these functions depends partly on boundary treatment methods including the boundary function method [32], the Lagrange multiplier method [22], and the artificial spring method [26, 27]. The former procedure directly ensures the satisfaction of natural boundary constrains by embedding the boundary equations into the admissible functions. The latter ones permit the application of the trial functions for unconstrained substructures and subdomains due to the potential satisfaction of desired boundary conditions in the variational principle, but an additional stiffness matrix needs to be derived and appended. Although only the boundary function method is adopted to deal with the boundary conditions here, a suitable combination of the above approaches according to the actual situations can be flexibly utilized to treat analytically the structure and acoustic subsystems constrained or supported in arbitrary manners, which reasonably reflects the high flexibility and expandability of the present procedure.With the present method, the explicit expression of derivative quantities, such as mechanical strain or stress, can be directly obtained by applying differential operators to the orthogonal polynomials. Since orthogonal polynomials are generally multiorder differentiable, derivative quantities can be obtained with the similar spatial resolution as the primary variables. It is also worthwhile to notice that, when the point loads are applied at the nodes in the FEM, the artificial stress/strain concentration usually appears at the exciting position due to the point load singularity. Using the second case as an example, the comparisons of spatial distribution ofy component of the mechanical strain between the present and FEM predictions at 610 Hz are plotted in Figure 13. The artificial strain concentration can be clearly observed in the contour plots of FEM model. However, it seems like that the present method has effectively attenuated this phenomenon. In fact, except the small areas of artificial strain concentration in FE analysis, the distribution of the mechanical strain is in good agreement between the proposed method and FEM. Considering an arbitrarily chosen point (0.365, 0.55), εyy is calculated to be 1.209×10−5 for the FEM, while εyy=1.207×10−5 for the proposed method.Figure 13
The spatial distribution ofy component of mechanical strain due to a point force of 610 Hz applied on (0.3 m, 0.6 m) Nx=Ny=13,Mx=13,My=13,NL=3.It should also be noted that most of the computation time within the present procedure is consumed in integration calculations in equation (25). The time cost is mainly affected by the number of quadrature nodes, the truncation orders of polynomial series, and the number NL of subdomains. Generally, the choice of truncation orders and quadrature points should satisfy some specified truncation rules in order to ensure a good convergence rate and numerical stability [19]. Taking these restrictions into account, the smaller number NL of subdomains is preferred in the case of a guaranteed convergence rate, and thus, the present method may be most efficient for industrial applications with moderately complex geometries.


## 4.1. A Simplified Vehicle Cabin Backed by Plate Assembly

To ensure the applicability and computational advantages of the present formulation for analyzing vibroacoustic problems with relatively complex geometries, a simplified vehicle cabin in Figure2 is employed; it consists of a built-up cavity backed by a plate assembly. The structure is assembled by three plate components oriented arbitrarily. Unlike the pure plate bending problems as in [15, 16], the coupling effects of the out-of-plane modes and in-plane modes of the plates should be considered here, especially for vibrations in higher frequencies [34]. Moreover, besides the transformation of the physical coordinate to the dimensionless local coordinate, the mapping of the irregular geometries into the regular ones should also be included:(32)xe=Lxe2ξe+Lxe2,ye=Lye2ηe+Lye2,ze=Lze2ζe+Lze2,ye=Lxetanαi41+ξe+1+ζe,where Lxe, Lye, and Lze are dimensions of subdomains in three directions (Lze is the dimension of the left surface of each acoustic subdomain in z). xe, ye, and ze are the local physical coordinates for each fluid subdomain. ξe, ηe, and ζe are the dimensionless local coordinates defined in the interval [−1, 1]. αi (i = 1, 2) is the leaning angle as in Figure 2. Following the standard isoparametric mapping procedure as used in FEM [7], the Jacobian matrix for the transformation can be readily achieved, and thus the functional in equations (1) and (20) can be mapped from the physical domain to computational domain. Substituting equation (32) into equation (25) finally yields the system of equations in the local coordinates.Figure 2
Interior vibroacoustic interaction of a simplified vehicle cabin. Point coordinates: A(0.8, 0.6, 1), B(0.5, 0.45, 0.5), C(1.0, 0.2, 1), and D(0.75, 0.45, 0.5).The uncoupled edges of all plates are taken to be clamped, while all uncoupled walls of the cavity are assumed to be rigid. Taking benefit of the generality of the present formulation, in addition to representing the out-of-plane motion of the plate by Kirchhoff’s plate bending theory, the theory of elasticity for two-dimensional solids [29] is also employed to formulate the in-plane motion of the plates. The interaction between the in-plane and out-of-plane vibrations occurs only at the coupling interface. The physical parameters for the system are the same as those used in Ref. [35], the plates are 1.5 mm thick and made of steel with Young’s modulus E = 210 GPa, the mass density ρs=7850kg/m3 and the Poisson ratio vs = 0.3, while the cavity is filled with air of fluid density ρf=1.225kg/m3 and the speed of sound cf = 340 m/s. During the calculation, the plate and the cavity are divided into three components in the x direction according to the plate junction positions. Both the in-plane and out-of-plane displacements should be approximated in terms of Chebyshev polynomials, and thus the total number of degrees of freedom (DoFs) is given by Nsum=3Nx⋅Ny⋅NL+Mx⋅Mx⋅Mx⋅NL.The ability of the present method to characterize well free vibration behaviors of the coupled system is first demonstrated. The prediction accuracy of the present methodology depends mainly on the truncation order (Nx,Ny,Mx,My,Mz) of polynomial functions; thus, a convergence analysis needs to be conducted. The simple way is to recalculate the concerned quantities with an increasing number of polynomial functions until the quantities cease to change. The 15-th and the 85-th natural frequencies (f15 = 80.53 Hz and f85 = 293.11 Hz), which are, respectively, located in the low-frequency and mid-frequency bands, are arbitrarily chosen for example. The convergence curves for the present method and FEM are shown depending on the number of degrees of freedom of the problem in Figure 3. A very fine FE analysis with the element sizes he = 0.013 m and the total number of DoFs 668,920, which ensures the required accuracy for comparison purposes, serves as reference (one can find the detailed information about the FEM model in Appendix B). The formula for the relative error is used to quantify the quality of the convergence:(33)εi=ωi−ωi,refωi,ref×100%,where ωi,ref is the i-th frequency obtained from the FE model with the fine mesh. The convergence curves for the present model are obtained by increasing the truncation order (Mx,My,Mz) of Chebyshev polynomial functions for the fluid from 4 to 14 while keeping Nx=Ny=14, while those for the FEM are achieved by refining element sizes from 0.06 m for both the fluid and structure domains until the FEM predictions accuracy at a similar accuracy as the present converged results. It is shown in Figure 3 that the present model can provide a stable and rapid convergence for both the lower- and higher-order frequencies with a rather small model size as compared to the FEM, where further increments of DoFs produce relatively little change in the relative error. Moreover, a good agreement, with relative errors less than 3×10−4 between the present converged results and FE predictions, can be guaranteed with a lower number of DoFs than the reference FEM, resulting in a considerably shorter computation time. In addition, the present method is able to represent well the mode shapes. For example, comparison of the mode shapes for the 85-th mode between the present and the reference FEM solutions in Figure 4 demonstrates that the spatial distribution of both the acoustic pressure and the displacement obtained by the present method agrees well with the FEM results. To highlight the benefits of the accuracy and efficiency of the present method, several natural frequencies for different combinations of leaning angles α1 and α2 are further presented in Table 1. It shows clearly that the present results are in good agreement with those obtained from the FE model, which may be utilized as benchmark solutions for the similar numerical studies.Figure 3
Convergence curves for the present method and FEM with regard to the 15-th and 85-th natural frequencies. A FE model with element sizeshe = 0.013 m is used to provide a base of comparison (f15 = 80.53 Hz, f85 = 293.11 Hz, and NL = 3).Figure 4
The mode shapes of the structure displacements (top panels) and the acoustic pressure (bottom panels) of the cavity with regard to the 85-th natural frequencies. (a) Present solutions. (b) FEM solutions.
(a)(b)Table 1
Comparison of natural frequencies in Hz for the built-up system in Figure2 with different leaning angles α1 (Nx=Ny=14, Mx=My=Mz=8, and the leaning angle α=π/4).
α1 (rad)MethodMode no.12345203040α1 = atan (6/5)Present16.1622.4434.1535.5135.7597.80130.7166.6FEM16.1622.4434.1435.5035.7497.75130.5166.4α1 = π/3Present14.7721.3126.5631.9234.9289.98123.1152.2FEM14.7621.3126.5531.9134.9189.91123.0151.9α1 = π/4Present16.8723.4834.5337.0037.1999.82143.8174.5FEM16.8723.4734.5236.9937.1899.77143.5174.1α1 = π/6Present18.0127.5134.7338.4242.23117.6149.6189.1FEM18.0127.5134.7238.4042.21117.5149.3189.1The reliability and efficiency for numerical calculations of steady-state vibroacoustic responses are then checked, in which a unit point forceFs in the z direction is applied at point A (see Figure 2). The sound pressure responses at point B and the displacement responses at point C are collected. Comparisons of the response amplitude spectrum between the present results and the reference FEM solutions with the element sizes he = 0.013 m are plotted in Figure 5 between 1 and 400 Hz at intervals of 1 Hz. Note that both the flexural wave-number of plates and the acoustic wave-number in the fluid are considerably larger than one, and it is not a stretch to infer that both the low-frequency and mid-frequency behaviors of the coupled system can be exhibited within the considered frequency range. A good agreement between the results demonstrates the reliability of the present method over the frequency range of interest. The discrepancies in magnitude of the resonant peaks are mainly because (1) no damping is introduced in both simulations and (2) a small difference in natural frequency may cause a large difference in amplitude close to the resonant frequency. It should be mentioned that the observation point A in Figure 5(b) is located at the interface between the fluid subdomains ωf1 and ωf2; thus, the consistency between the present results and the converged FE solutions also reveals that the constraint conditions between adjacent subdomains and substructures can be incorporated correctly by using the Lagrange multiplier technology. To illustrate that the present method is capable of providing accurate spatial distribution of response fields in an explicit fashion, detailed comparisons of pressure amplitude on several lines inside the cavity between the present method and the reference FEM are plotted in Figure 6. Additionally, Figure 7 shows a color plot of the acoustic pressure and the structural deformation fields at 225 Hz. Again a good match between the results of both modeling approaches is observed. For the current calculations, the truncation orders Nx=Ny=14 of polynomial functions for each substructure and Mx=My=Mz=8 for each fluid subdomain are used; this leads to a model with a total of 3300 DoFs and the entire solution time for a single frequency of the present method is about 62 s. Moreover, following the convergence analysis in Figure 3, a FE model with 334,476 DoFs (and corresponding element sizes he = 0.016 m) is needed to gain the similar accuracy as the present method. With the same hardware, the calculation time of this FE model is approximately 1324 s, revealing an obvious computational cost saving of the proposed procedure under the similar level of accuracy.Figure 5
(a) Displacement amplitude responses at point C and (b) acoustic pressure amplitude responses of point B at the interface between the fluid subdomainsΩf1 and Ωf2Nx=Ny=14,My=Mz=8,NL=3.
(a)(b)Figure 6
Comparison of pressure amplitude on Line 1 (y = 0.3 m, z = 0.45 m) and Line 2 (y = 0.5 m, z = 0.45 m) obtained from the present method and FEM due to a point force excitation of 225 Hz Nx=Ny=14,My=Mz=8,NL=3.Figure 7
The spatial distribution of the amplitude of the structure displacements and the acoustic pressure of the cavity due to a point force excitation of 225 Hz applied on the center plate at point ANx=Ny=14,My=Mz=8,NL=3. (a) Present solutions. (b) FEM solutions.
(a)(b)The present procedure is able to introduce acoustic impedance boundaries. A simulation is subsequently followed to calculate the dynamic responses of the system with the wall impedance. Although the method can be easily utilized to take into account different impedance boundaries, only the nondissipative boundary defined by pure imaginary impedance is employed here for demonstration purposes. Considering the right surface of the cavity in Figure2 as the nondissipative boundary with uniform impedance Zf=102j, comparisons of the plate responses and sound pressure between the present and FEM solutions from 1 Hz to 400 Hz are, respectively, plotted in Figure 8 at the intervals 1 Hz. The plate response point is also selected at point C, while the pressure response point is set at point D inside the second fluid subdomain. It is obvious that both the acoustic pressure and the displacement responses obtained by the present solution well match the reference FEM results. In this case, using the FE model (he = 0.016 m) under the similar level of accuracy as the present results, required time for a sparse solution is about 1587 s per each calculation step. However, only an average of 0.32 s is used for each frequency in the harmonic analysis within the present procedure, further illustrating a significant saving of the calculation effort of the present procedure for solutions of vibroacoustic problems with acoustic impedance boundaries.Figure 8
(a) Normal displacement at (1.0 m, 0.2 m, 1 m) and (b) acoustic pressure at (0.75 m, 0.45 m, 0.5 m). Nondissipative boundary with a uniform impedanceZf=102j is applied to the right surface of the cavity Nx=Ny=14,My=Mz=8,NL=3.
(a)(b)


## 4.2. Irregular Elastic Solid with Heavy Acoustic Medium

An irregular 2D elastic structure backed by heavy acoustic medium is finally considered to further illustrate that the present formulation can expediently combine with different theories of elastic structures and easily accommodate problems in different dimensionalities. The geometric configuration as well as the dimensions of the system is shown in Figure9. The right and left edges of the solid are taken to be clamped, and all the uncoupled walls of the cavity are assumed to be rigid. The cavity is filled with water of density ρf=1000kg/m3 and acoustic speed cf = 1480 m/s. The structure is made of rubber with Young’s modulus E = 23 MPa, mass density ρs=1600kg/m3, Poisson ratio vs = 0.47, and the inclined angle α=π/12 rad. For comparison purposes, the vibroacoustic system is modeled by both the FEM and the proposed spectral method. The description of the FEM model can be found in Appendix B. For the present model, because the displacement depends on both x and y, a plate model would not be appropriate, so a 2D elastic model which is possible due to the generality of the present formulation is necessary here. And the exact small-strain 2D elasticity theory [29] is employed to model the structure. The irregular structure and the cavity are partitioned into three components NL=3 in the x direction as shown in Figure 9, and the irregular substructure can be easily mapped into the regular one similarly to the previous case. Moreover, approximating the velocity potential of the fluid and the displacements of the solid in the x and y directions in terms of spectral Chebyshev polynomials, the total number of DoFs for the present model can be given by Nsum=2Nx⋅Ny⋅NL+Mx⋅My⋅NL, where NL is the number of subdomains and Nx, Ny, Mx, and My are the truncation orders of the polynomial series.Figure 9
An irregular 2D elastic solid coupled with heavy medium.Considering modal characteristics first, convergence curves of lower-order (35-th,f35 = 490.09 Hz) and higher-order (65-th, f36 = 714.73 Hz) natural frequencies are shown in Figure 10 as examples. The global curve for the present model is performed by increasing the polynomial order (Mx=My) of the velocity potential from 5 to 20 and fixing the number of the polynomial (Nx=Ny=13) with regard to structure displacements, while the curves for the FEM are obtained by gradually refining the mesh sizes from 0.025 m until a similar accuracy as the present converged results reaches. The FE model with a very fine mesh (element sizes he = 0.0025 m and the total DoFs 175,861) is utilized in Figure 10 to provide a base of comparison, and the definition of the relative error is the same as that in equation (33). It can be seen that the accuracy of the FEM model increases steadily with an increasing number of DoFs, while the relative errors for the present model slump at the beginning and converge rapidly towards a situation where the increase of polynomial terms does not produce further improvements of prediction accuracy. An excellent agreement with relative errors on the order of 10−4 between the present converged results and FE predictions for both the lower- and higher-order frequencies can be observed, which shows that even with small models, the proposed method can provide accurate solutions within a fairly large frequency range. Note that, although the convergence analysis is only carried out for several selected frequencies, it is expected the convergence property observed is valid for other frequencies.Figure 10
Convergence curves for the present method and FEM with regard to the 35-th and the 65-th natural frequencies (f35=490.09Hz,f65=714.73Hz,NL=3).Finally, the efficiency for calculations of steady-state responses is verified, in which a unit harmonic forceFs in the y direction is applied on the top surface of the elastic structure at the global position (x,y) = (0.3 m, 0.6 m). Figures 11(a) and 11(b) present displacement responses in the y direction at (0.6 m, 0.54 m) and sound pressure responses at (0.7 m, 0.3 m) between 1 and 700 Hz at intervals of 1 Hz for the present method, respectively. The FEM reference solutions with element sizes he = 0.0025 m are also presented in Figure 11 for method verification. Good correlation illustrates that the developed method is able to accurately predict the vibroacoustic behaviors at relatively high frequencies. Additionally, comparisons of spatial distribution of the structural deformation field in the x direction between the present and FEM predictions at a higher frequency 610 Hz are further plotted in Figure 12. As can be observed, the present result corresponds well with the reference result. The contour plots of the pressure as well as displacements in the y direction are also in good agreement, which are not shown for brevity. For the present analysis, truncation orders Nx=Ny=13, Mx = 13, and My = 7 of polynomial functions for subdomains are adopted in the calculations, and the dimension of resulting system matrices is only 1287. It is much smaller than that of FEM analysis under the similar level of accuracy, where 123,292 DoFs is needed for the FE model as can be seen from Figure 10. Furthermore, the present results in Figure 11 are calculated about 15.2 times faster than the FEM sparse solution with the similar accuracy. The comparisons not only show the reliability and modeling efficiency of the present method even if both the irregular elastic solid and the heavy acoustic medium are introduced but also indicate that the proposed multidomain spectral procedure is applicable for vibroacoustic analyses in a broad frequency range.Figure 11
(a) Displacement iny direction at (0.6 m, 0.54 m) and (b) acoustic pressure at (0.7 m, 0.3 m) inside the cavity Nx=Ny=13,My=Mz=13,NL=3.
(a)(b)Figure 12
The spatial distribution of the structure displacements in thex direction due to a point force of 610 Hz applied on the structure at position (0.3 m, 0.6 m) Nx=Ny=13,My=Mz=13,NL=3.


## 4.3. Discussion

It may be concluded from the above analyses that the present method may offer a very powerful tool for providing an efficient modeling approach for vibroacoustic problems with relatively complicated geometries. By inheriting the flexibility from the multidomain synthesis as used in structural mechanics, the approach divides the full continuous system into a few large continuous subdomains instead of a large number of small elements as in FEM, which results in substantially smaller system matrices and higher computational efficiency under the similar level of prediction accuracy as compared to FEM. Moreover, the application of multidomain synthesis makes the analytical techniques possible to analyze vibroacoustic systems consisting of an assembly of structures of arbitrary directions or thick structures requiring general elasticity theory.Benefiting from the versatility of the spectral method, the present formulation mathematically permits the use of any orthogonal polynomials (Fourier, Legendre, Chebyshev, Hermite polynomials, etc.), which possess complete, linearly independent, and highly smooth properties, as admissible functions to approximate field variables. Although these polynomial functions can usually bring about rapid convergence and superior stability in numerical operations, they are required to satisfy the desired boundary conditions, such as prescribed displacement boundaries. The selection of these functions depends partly on boundary treatment methods including the boundary function method [32], the Lagrange multiplier method [22], and the artificial spring method [26, 27]. The former procedure directly ensures the satisfaction of natural boundary constrains by embedding the boundary equations into the admissible functions. The latter ones permit the application of the trial functions for unconstrained substructures and subdomains due to the potential satisfaction of desired boundary conditions in the variational principle, but an additional stiffness matrix needs to be derived and appended. Although only the boundary function method is adopted to deal with the boundary conditions here, a suitable combination of the above approaches according to the actual situations can be flexibly utilized to treat analytically the structure and acoustic subsystems constrained or supported in arbitrary manners, which reasonably reflects the high flexibility and expandability of the present procedure.With the present method, the explicit expression of derivative quantities, such as mechanical strain or stress, can be directly obtained by applying differential operators to the orthogonal polynomials. Since orthogonal polynomials are generally multiorder differentiable, derivative quantities can be obtained with the similar spatial resolution as the primary variables. It is also worthwhile to notice that, when the point loads are applied at the nodes in the FEM, the artificial stress/strain concentration usually appears at the exciting position due to the point load singularity. Using the second case as an example, the comparisons of spatial distribution ofy component of the mechanical strain between the present and FEM predictions at 610 Hz are plotted in Figure 13. The artificial strain concentration can be clearly observed in the contour plots of FEM model. However, it seems like that the present method has effectively attenuated this phenomenon. In fact, except the small areas of artificial strain concentration in FE analysis, the distribution of the mechanical strain is in good agreement between the proposed method and FEM. Considering an arbitrarily chosen point (0.365, 0.55), εyy is calculated to be 1.209×10−5 for the FEM, while εyy=1.207×10−5 for the proposed method.Figure 13
The spatial distribution ofy component of mechanical strain due to a point force of 610 Hz applied on (0.3 m, 0.6 m) Nx=Ny=13,Mx=13,My=13,NL=3.It should also be noted that most of the computation time within the present procedure is consumed in integration calculations in equation (25). The time cost is mainly affected by the number of quadrature nodes, the truncation orders of polynomial series, and the number NL of subdomains. Generally, the choice of truncation orders and quadrature points should satisfy some specified truncation rules in order to ensure a good convergence rate and numerical stability [19]. Taking these restrictions into account, the smaller number NL of subdomains is preferred in the case of a guaranteed convergence rate, and thus, the present method may be most efficient for industrial applications with moderately complex geometries.


## 5. Conclusion

The principal contribution of this paper is the development of a unified multidomain spectral formulation to characterize vibroacoustic properties of interior structural-acoustic systems with moderate geometrical complexity. Combining the versatility of the spectral method and the flexibility of generalized variational principle, the established framework enables one (1) to expediently accommodate complex fluid-structural systems consisting of both a structure assembly and a built-up cavity, (2) to conveniently tackle arbitrary combinations of boundary conditions and problems in different dimensionalities, and (3) to effectively analyze vibroacoustic characteristics with reliable accuracy at relatively high frequencies.The second contribution is the detailed calculation validation and convergence studies of the given procedure. Comparisons between the present method and FEM are performed for two practical problems of moderate complexity, namely, a simplified vehicle cabin backed by plate assembly and an irregular solid coupled with heavy acoustic medium. The good agreements for both the eigenvalues and steady-state responses demonstrate the reliability and efficiency of the present formulation. With relatively small model size, the developed procedure is found to provide a rapid convergence for a considerable frequency range and a significant computational advantage as compared to FEM. In fact, the computational advantage can be exploited to significative parametric studies for the design and control of vibroacoustic systems in engineering practices.



---

*Source: 6910174-2020-05-29.xml*