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

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.*e unified formulation is derived by means of a generalized variational principle in conjunction with the spectral discretization procedure. *e 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.


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. e finite element method (FEM) and boundary element method (BEM) are the most widely utilized numerical procedures [3][4][5][6] that have been used to analyze vibroacoustic problems. e superiority of elementbased 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][15][16], and the Trefftz methods [17,18], have also been used to analyze structural acoustics and vibration behaviors. e 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. is study attempts to demonstrate the capability of analytical techniques to accommodate built-up structural-acoustic systems of moderate complexity. e 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. us, 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. e 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. e 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. e 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. e 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. e 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.

Statement of the Problem
e general vibroacoustic problem studied here is schematically sketched in Figure 1, where the fluid domain Ω f is partially in contact with the structural domain Ω s at the interface Γ fs . e structure is supposed to be linear-elastic, and the structural boundaries with prescribed displacements and surface forces are denoted by Γ us and Γ σs , respectively. e compressible fluid is assumed to be homogeneous, irrotational, and inviscid, and 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 . n f and τ s are the unit vectors normal to the fluid and structural surfaces, respectively.
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 N s nonoverlapping substructures and N f (N f ≥ N s ) 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 Ω e s and acoustic subdomains Ω e f , and for the components without fluid-structure interaction Ω e s � 0. Moreover, the responses of any subdomain in the global coordinate system (o-xyz) can be independently described by the associated local coordinate system (o e -x e ·y e ·z e ), where o and o e 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.

Theoretical Formulations
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, W e si is the i-th component of the displacement, ρ s is the mass density of the structure, σ e ij and ε e ij denote the stress and strain tensor, respectively, and F e si are prescribed surface forces and the _ ϕ e f n e fi dΓ e fs , (2) where n e fi is the unit normal external to Ω fs and ϕ e f is the velocity potential.
Performing variational calculations of equation (1) and applying the stationary condition of J s , yields δJ s � For a general elastic structure, the strain-displacement relations and strain-stress relations can be given according to the linearized elasticity theory as [29] ε e ij � where τ e j is the outward drawn normal to the substructure surface e s . Since δλ e,e+1 si and δw e si are assumed to be arbitrary variations in Ω e s , the following relations can be easily derived from equation (6): (1) the governing equations, (2) the prescribed force boundary conditions on Γ e σs , (3) the continuity conditions on the fluid-structure interface Γ e fs , (4) the continuity conditions on the substructures interfaces Γ e,e+1 s , and (5) Lagrange multipliers derived in terms of original variables, namely, e Lagrange multipliers λ e,e+1 si 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 e,e+1 λs 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 σ e ij and ε e ij in equations (1) and (3) should be replaced by stress resultants τ e and generalized strains ε e for the plate, and where the superscript T denotes the transpose, ε e s � [ε e x ε e y ε e xy ] is the normal and shear strains of the neutral surface, Following the thin plate theories [30], the strain-stress relations and strain-displacement relations can be given as where K e ε and K e κ are the extensional and bending stiffness, D ρ is a differential operator derived from plate theory, and d e p � [w e s1 w e s2 w e s3 ] is midsurface displacements of the plate in the three local coordinate directions.
Moreover, the interface potential W e,e+1 λs should be rewritten according to thin plate theories as Here, the matrix expression within the bracket is the displacement constraint equations on the common interface, and λ e,e+1 su λ e,e+1 sv λ e,e+1 sw λ e,e+1 sϕ ] T is the vector of Lagrange multipliers defined on the interface between the e-th and (e + 1)-th plates.
e differential operator C and the where α e is the angle between the x e -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 where V e x is the Kevlin-Kirchhoff edge reaction (see Ref. [30]). A similar procedure can work for the case of thin shell assemblies without any difficulties. e 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. 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 smallamplitude oscillations in e-th fluid subdomain can be given as [31]

Acoustic Vibration
where u e f , ϕ e f , and p e f denote, respectively, the particle displacement, velocity potential, and sound pressure. e constant c o 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.
e 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 (W e zf , W e vf , and W e sf ) with regard to the impedance boundary, the Neumann boundary and the fluid-structure interface can then be, respectively, expressed as follows: _ ϕ e f w e s · n e f dΓ e sf , (17) where ω is the angular frequency, n in the subscripts denotes the outward normal to the fluid boundaries, and j � �� � −1 √ ; w e s is the displacement vector of the structure and n e f is unit outward normal. Z e f and v e f are the specific acoustic impedance and the externally imposed particle velocity, respectively. e adjacent subdomains can be synthesized by introducing a Lagrange multiplier λ e,e+1 f defined along the interface λ e,e+1 f . e additional interface energy can be consequently obtained by multiplying the constraint relationship by the undetermined Lagrange multiplier as e Dirichlet boundary herein is treated as a special interface as those between adjacent segments by replacing ϕ e f with −p −e nf /jρ f ω, in which p −e nf is the prescribed acoustic pressure. Moreover, the work done by the harmonic volume distributed acoustic source can be given as where Q e f is the distributed strength of the volume velocity source. e corresponding variational functional can then be written as Invoking the stationarity of J f with respect to independent variables λ e,e+1 f and ϕ e f gives the principle of virtual potentials for the fluid. erefore, the first variation of J f with respect to small perturbation λ e,e+1 f and δϕ e f 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 and n e f denotes the unit normal external to the fluid boundary surface e f . Since δλ e,e+1 f and δϕ e, f are arbitrary, expressions 1-5 in equation (21) 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 λ e,e+1 f and giving a weak form expressed entirely in terms of ϕ ε f , which is not shown here for brevity. Here, λ e,e+1 f is physically identified as the interface normal velocity.

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). e i-th (i � x, y, z) displacement w e i for the e-th substructure and the velocity potential ϕ e f for the e-th subdomain can be, respectively, approximated in terms of truncated boundary characteristic polynomial series [14,20]: � Φ e x e , y e , z e q e f e jωt , (24) where N x , N y , M x , M y , and M z are the truncation orders of the polynomial series, P α (β)(1, 2, 3 · · · ; β � x e , y e , z e ) are the α order polynomials, d e i,kl (t) and q e rsp (t) are generalized coordinates for the e-th substructure and fluid subdomain,  [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 q e f and d e si readily result in the following assembled system of equations: where . e subscripts f and s will be generated by relevant submatrices with regard to the fluid and structure, respectively. e matrices M δ and k δ (δ � f, s) are equivalent mass and stiffness matrices. e global acoustic-structural coupling matrix is denoted as C fs , and Z f is the matrix related to the impedance boundary of the fluid. e matrices K Lf and K LS 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 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. e steady-state response of q can be readily given by assuming the harmonic excitation as e 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 p e i x e , y e , t � jωρ f Φ x e , y e , Z e q e f e jωt . (28) Obviously, using the state vector 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 and the i-th eigenvalue λ i can be written in the following form: where the notation * denotes complex conjugate, ζ i and ω i are the i-th modal damping ratio and frequency, which can be expressed as

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. e 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. e 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. e 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.

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 Figure 2 is employed; it consists of a built-up cavity backed by a plate assembly. e 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: e 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.
e interaction between the in-plane and out-ofplane vibrations occurs only at the coupling interface. e 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 � 7850 kg/m 3 and the Poisson ratio v s � 0.3, while the cavity is filled with air of fluid density ρ f � 1.225 kg/m 3 and the speed of sound c f � 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 e ability of the present method to characterize well free vibration behaviors of the coupled system is first demonstrated. e prediction accuracy of the present methodology depends mainly on the truncation order (N x , N y , M x , M y , M z ) of polynomial functions; thus, a convergence analysis needs to be conducted. e simple way is to recalculate the concerned quantities with an increasing number of polynomial functions until the quantities cease to change. e 15-th and the 85-th natural frequencies (f 15 � 80.53 Hz and f 85 � 293.11 Hz), which are, respectively, located in the low-frequency and mid-frequency bands, are arbitrarily chosen for example. e 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 h e � 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). e formula for the relative error is used to quantify the quality of the convergence: where  Shock and Vibration 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. e reliability and efficiency for numerical calculations of steady-state vibroacoustic responses are then checked, in which a unit point force F s in the z direction is applied at point A (see Figure 2). e 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 h e � 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. e 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 N x � N y � 14 of polynomial functions for each substructure and M x � M y � M z � 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 e 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 Figure 2 as the nondissipative boundary with uniform impedance Z f � 10 2 j, 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. e 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 (h e � 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.

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.
e geometric configuration as well as the dimensions of the system is shown in Figure 9. e 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. e cavity is filled with water of density ρ f � 1000 kg/m 3 and acoustic speed c f � 1480 m/s. e structure is made of rubber with Young's modulus E � 23 MPa, mass density ρ s � 1600 kg/m 3 , Poisson ratio v s � 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. e 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 smallstrain 2D elasticity theory [29] is employed to model the structure. e irregular structure and the cavity are partitioned into three components (N L � 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 N sum � 2N x · N y · N L + M x · M y · N L , where N L is the number of subdomains and N x , N y , M x , and M y are the truncation orders of the polynomial series.  Considering modal characteristics first, convergence curves of lower-order (35-th, f 35 � 490.09 Hz) and higherorder (65-th, f 36 � 714.73 Hz) natural frequencies are shown in Figure 10 as examples. e global curve for the present model is performed by increasing the polynomial order (M x � M y ) of the velocity potential from 5 to 20 and fixing the number of the polynomial (N x � N y � 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. e FE model with a very fine mesh (element sizes h e � 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.
Finally, the efficiency for calculations of steady-state responses is verified, in which a unit harmonic force F s in the   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. e 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 N x � N y � 13, M x � 13, and M y � 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. e 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.

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    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]. e former procedure directly ensures the satisfaction of natural boundary constrains by embedding the boundary equations into the admissible functions. e 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 of y component of the mechanical strain between the present and FEM predictions at 610 Hz are plotted in Figure 13. e 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.
It should also be noted that most of the computation time within the present procedure is consumed in integration calculations in equation (25). e time cost is mainly affected by the number of quadrature nodes, the truncation orders of polynomial series, and the number N L 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 N L 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.

Conclusion
e 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 fluidstructural 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.
e 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. e 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.
where M e f , K e f , and C e fs are the mass matrix, the stiffness matrix, and the coupled fluid-structural matrix, respectively. F e f is the distributed acoustic source vector and K e,e+1 Lf is the generalized interface stiffness matrix.
Finally, the assembly of derived system matrices for all substructures and fluid subdomains leads to the equations of motion equation (25) for the entire vibroacoustic interaction system. Moreover, the Gauss-Legendre integration method is used to calculate the system matrices in this work and the number of the quadrature nodes is set to be 15, which is sufficient for reaching good solution accuracy for the considered numerical cases.

B. Derivation of System Matrices
e finite element models needed for the validation of the present method were built and analyzed using the commercial software ANSYS 17.1. For the first case, the acoustic subdomains and substructures are discretized by using 8node acoustic fluid elements (FLUID30) and 4-node elastic shell elements (SHELL63), respectively. e acoustic fluid elements only have the pressure DoFs, while the interface elements (in purple) in contact with the plate include both pressure and displacement DoFs. Both effects of the bending and membrane stiffness of the plates are included in FEM. For the second case, the acoustic fluid and substructures are discretized by using 4-node acoustic elements (FLUID29) and 8-node quadrilateral plane stress elements (PLANE183), respectively. Beside the solid, the displacement DoFs are only kept for the interface acoustic elements in contact with the solid.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.