# The Variational Theory of Complex Rays: An answer to the resolution of mid-frequency 3D engineering problems

Hervé Riou, Pierre Ladevèze, Louis Kovalevsky

► To cite this version:

Hervé Riou, Pierre Ladevèze, Louis Kovalevsky. The Variational Theory of Complex Rays: An answer to the resolution of mid-frequency 3D engineering problems. Journal of Sound and Vibration, 2013, 332 (8), pp.1947-1960. 10.1016/j.jsv.2012.05.037. hal-01647868

HAL Id: hal-01647868

https://hal.science/hal-01647868

Submitted on 13 Jan 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The Variational Theory of Complex Rays: An answer to the resolution of mid-frequency 3D engineering problems

H. Riou, P. Ladevèze*,<fnref n="1" />, L. Kovalevsky

LMT-Cachan (ENS Cachan/CNRS/Paris 6 University, PRES UniverSud Paris) 61 avenue du Président Wilson, F-94230 Cachan, France

The Variational Theory of Complex Rays (VTRC) is an approach to the simulation of mid-frequency phenomena whose wavelengths are relatively small compared to the dimensions of the domain. This is a wave-based computational technique which involves a nonclassical variational formulation. This paper focuses on the development of the approach for 3D engineering problems and shows that it is a mature technique. Illustrations are given to show the capabilities of this method.

## 1. Introduction

Due to regulatory requirements and growing customer expectations, the vibrational and acoustical behavior of structures has become an important issue in engineering design. Virtual prototypes are being increasingly used in order to reduce the number of (expensive) experimental tests. This requires efficient and accurate numerical prediction techniques.

Today, most engineering problems are solved using the Finite Element Method (FEM). However, because the FEM seeks an approximate solution by using continuous, piecewise polynomial functions, it is incapable of handling problems involving rapidly oscillating functions other than by using highly refined meshes. High-order finite elements can be used as an alternative to a very fine mesh. However, due to the pollution effect, the number of degrees of freedom must still be appropriately increased for the standard FEM to maintain the desired level of accuracy. Thus, in practice, the FEM is limited to low-frequency applications. Although various techniques have been proposed in order to overcome this limitation (see [1-7]). Indeed, these techniques have been shown to reduce computation costs compared to standard finite elements, but the frequency ranges which they are capable of addressing are still lower than that would be desirable.

Another often used prediction technique is the Boundary Element Method (BEM). Since the BEM is based on a boundary integral formulation of the problem, only the boundary of the domain being considered needs to be discretized. Within the boundary elements which compose the model, the variables are expressed in terms of simple polynomial shape functions. Unlike the FEM, the enforcement of the boundary conditions results in only a small system of equations to be solved for the nodal values at the discretized boundary. Once these nodal values have been determined, the field variables within the domain can be reconstructed in a postprocessing step using the boundary integral formulation. Since only the boundary needs to be discretized, the BEM is exempt from some of the limitations of the FEM mentioned previously. Besides, many works have been dedicated to the rapid evaluation of the necessary integrals (see [8,9]). However, these methods are

<footnote>* Corresponding author. Tel.: +33 147402238; fax: +33 147402240.</footnote>

<footnote>E-mail addresses: riou@lmt.ens-cachan.fr (H. Riou), ladeveze@lmt.ens-cachan.fr (P. Ladevèze), kovalevsky@lmt.ens-cachan.fr (L. Kovalevsky).</footnote>

<footnote n="1">¹ EADS Foundation chair, Advances Computational Structural Mechanics.</footnote>
penalized by system matrices which are fully populated, frequency-dependent, complex and sometimes nonsymmetric, leading to calculations which are computationally intensive.

Another approach called Statistical Energy Analysis (SEA) (see [11,10]) involves the description of the energy exchanged among various systems and leads to the dynamic response of each system averaged over time and space. However, because it is characterized by a single energy level, SEA cannot account for the spatial variation of the response within each system and is subject to some restrictions (see [12]). Extensions of SEA have been developed (see [13–17]), but still produce information which is less local than would be desirable in the mid-frequency regime.

Today, there are also computational strategies dedicated to the resolution of mid-frequency problems, known as Trefftz methods (see [18]), which differ from the FEM in that the shape functions are often exact solutions of the governing differential equations. These approaches include, for example, a particular use of the partition of unity method [19–21], the ultra-weak variational method [22], the least-squares method [23], the discontinuous enrichment method [24], the element-free Galerkin method [25], the wave boundary element method [26], and the wave-based method [27]. All of these numerical techniques use oscillating functions to solve the problem. They differ primarily in the treatment of the transmission conditions among the substructures and boundary conditions, and in the types of shape functions used. The Variational Theory of Complex Rays (VTCR), which is the subject of this paper and was first introduced in Ref. [28], belongs to that category.

In [29,30], the VTCR was used to predict the vibrational response of a 3D plate assembly. In [31], plates with heterogeneities were taken into account. In [32], the approach was extended to shell structures. Vibration analysis over a range of frequencies was presented in [33]. The use of the VTCR for transient dynamic problems was addressed in [34]. The extension to linear acoustic problems was introduced in [35] and an adaptive version was developed in [36]. Many examples have shown that this approach is capable of producing accurate solutions using only limited numbers of DOFs. More recently, a new discretization of the repartition of the waves based on Fourier series has led to the development of an efficient semi-analytical method for calculating the operator [37].

After a general presentation of the basic concepts of the VTCR in Section 2, the paper focuses on the development of the VTCR for 3D engineering problems and shows that it is a competitive numerical prediction tool for the analysis of 3D structures made of plates (see Section 3) or shells (see Section 4), as well as for 3D acoustic problems (see Section 5). Conclusions are drawn in Section 6.

## 2. Basic concepts of the VTCR

Let us consider a general steady-state dynamic problem in $\Omega$, a $d$-dimensional domain with boundary $\partial\Omega$. The problem to be solved can be expressed as: find $u$ such that

$$ \mathcal{L}(u) = 0 \quad \text{in } \Omega \qquad (1) $$

$$ \mathcal{B}(u) = b \quad \text{over } \partial\Omega \qquad (2) $$

where $b$ denotes a prescribed function associated with prescribed sources located at the boundary of $\Omega$; $\mathcal{L}$ is a differential operator associated with the governing equation; and $\mathcal{B}$ is a differential operator associated with the boundary conditions (and also, if $\Omega$ is partitioned into subregions, with the continuity across the interfaces between these subregions).

The first characteristic feature of the VTCR is that it uses a global formulation of boundary conditions (2) in terms of both primal and dual quantities. Thus, Problems (1) and (2) are expressed as: find $u \in S_{ad}$ such that

$$ (u, \delta u)_{\mathcal{B}} = (b, \delta u) \quad \forall \delta u \in S_{ad} \qquad (3) $$

where $(\bullet, \bullet)_{\mathcal{B}}$ and $(b, \bullet)$ are respectively a bilinear form and a linear form equivalent to (2), and $S_{ad}$ is the space of the functions which satisfy (1). Formulation (3) is written in such a way that approximations which are a priori independent of one another can be used in the different subregions of $\Omega$.

Then, all that is required to build a VTCR approximation is the definition of a subspace $S_{ad}^a$ of $S_{ad}$. The approximate formulation can be written as: find $u^a \in S_{ad}^a$ such that

$$ (u^a, \delta u^a)_{\mathcal{B}} = (b, \delta u^a) \quad \forall \delta u^a \in S_{ad}^a \qquad (4) $$

The second characteristic feature of the VTCR is the description of the space $S_{ad}^a$ using two-scale approximations with a strong mechanical content: the solution of (1) is considered to be the superposition of an infinite number of plane waves which satisfy the governing equation exactly. Then, the construction of space $S_{ad}^a$ is the result of an approximation of the distribution of the amplitudes of the waves, which can be chosen among different options in order to benefit from physical and numerical advantages.

Once the discretization of each subregion of $\Omega$ has been chosen, the VTCR leads to a linear system of equations

$$ \mathbf{KA} = \mathbf{F} \qquad (5) $$

where matrix $\mathbf{K}$ corresponds to the discretization of the bilinear form of weak formulation (4), $\mathbf{A}$ is the vector of the unknown quantities which approximate the distribution of the wave amplitudes, and $\mathbf{F}$ corresponds to the discretization of the linear form of (4) projected onto $S_{ad}^a$. The approximate VTCR solution of Problem (1) and (2) is obtained by selecting from $S_{ad}^a$ the functions whose amplitudes are $\mathbf{A}$.
Several remarks can be made:

* **Weak formulation:** The VTCR is based on a specific weak formulation of the problem which was developed in order to enable the approximations within the subregions to be a priori independent of one another and which takes into account both the boundary conditions and the continuity across the interfaces, thus eliminating the need for a special treatment to guarantee interelement continuity. Any type of shape function can be used within each subregion, which gives this approach great flexibility and robustness.

* **Injected waves:** Space $S_{ad}^d$ is spanned by propagative and evanescent waves which are viewed as two-scale approximations. Only the slowly varying scale (which corresponds to the amplitudes of the waves) is discretized. The rapidly varying scale (which corresponds to the spatial shape of the waves) is taken into account analytically in Formulation (4). These two scales give the approximation a strong physical content and several computational advantages.

* **Dimension of subspace $S_{ad}^u$:** The number of waves to be used in $S_{ad}^u$ was studied in [30] for plates, in [32] for shells and in [37] for acoustics. Each of these studies showed that the VTCR, along with all other Trefftz methods, converges rapidly when the number of wave DOFs increases. Heuristic criteria and energy criteria for determining a priori the number of propagative waves required for an accurate numerical solution were presented in the last two references. Each of these criteria involves a relation between that number of waves and the ratio between the wavelength and the characteristic length of the domain being studied. Concerning evanescent waves, [30] showed that these are essential, but that not many of them are needed. If they are related to higher wavelengths than the propagative waves (e.g. in the case of pure acoustics problems), they are not even introduced at all.

* **Properties of matrix $K$:** Matrix $K$ in (5) is complex, nonsymmetric and frequency-dependent. However, since the shape functions are defined over acoustic subregions $\Omega_e$, $K$ is populated by blocks. Each block $K_{(i,i)}$ associated with subregion $\Omega_i$ is fully populated, but the off-diagonal submatrix $K_{(i,j)}$ between subregions $\Omega_i$ and $\Omega_j$ is zero if these subregions are not connected.

* **Accuracy of derived variables:** In the FEM, the primary response variables are usually approximated by simple polynomial shape functions. Therefore, the quantities formed by derivation of these variables (such as stresses and strains, acoustic velocities, etc.) are less accurate than the variables themselves. In the VTCR, there is no loss of accuracy because the derivatives of the wave functions are also wave functions.

* **Treatment of problems with complex geometries:** The VTCR can handle complex geometries, whose effects are included in (3), but sometimes a proper definition of $S_{ad}$ (and, thus, of $S_{ad}^u$) requires a specific discretization of $\Omega$ in the subregions $\Omega_e$.

* **Computational performance:** The calculation of the matrix coefficients involves the integration of highly oscillatory functions. Special care must be taken in carrying out these integrations in order to ensure good convergence. Semi-analytical methods can be used in the case of meshes composed of straight lines. Numerical integration schemes can be used for curved boundaries. Graphics Processing Units (GPUs) or parallel computer systems with multiple servers can also be used to accelerate many of the computations without resorting to low-level programming.

* **Handling of internal sources:** The presence of an internal source in (1) can be taken into account simply by adding a particular solution in approximation spaces $S_{ad}$ (see [29] for further details and some numerical examples).

# 3. 3D assemblies of plates

## 3.1. Problem description

In order to simplify the presentation of the VTCR for 3D plate assemblies, the problem will be formulated for an assembly of two plates, but this can be easily generalized. Let us consider two isotropic homogeneous plates whose reference surfaces are $\Omega$ and $\Omega'$ with boundaries $\partial\Omega$ and $\partial\Omega'$ respectively. In each plate, the quantities of interest are the normal displacements $w$ and $w'$ and the moment operators $M$ and $M'$. Fig. 1 presents the geometry of the problem along with the prescribed displacements $w_d$, rotations $w_{nd}$, bending moments $M_d$ and shear forces $K_d$ of $\Omega$ and the corresponding quantities of $\Omega'$. Using Kirchhoff-Love's plate theory, the steady-state out-of-plane displacements $w$ in $\Omega$ at circular

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Definition of the reference problem for an assembly of two plates.</figcaption></figure>
frequency $\omega$ are governed by the following partial differential equations:

$$ \Delta\Delta w \frac{(1+i\eta)Eh^3}{12(1-v^2)} = \rho h\omega^2 w \quad \text{in } \Omega \quad (6) $$

$$ \mathbf{M} = \frac{h^3}{12}(1+i\eta)\mathbf{K}_{\text{PS}}\mathbf{X}(w) \quad \text{in } \Omega \quad (7) $$

where $E, v, h, K_{\text{PS}}, \rho, \eta$ and $\mathbf{X}$ denote respectively Young's modulus, Poisson's coefficient, the plate's thickness, Hooke's plane stress tensor, the density, the damping coefficient and the curvature operator. Eqs. (6) and (7) are the governing equation and the constitutive relation of plate $\Omega$. Similar equations exist for plate $\Omega'$.

Let $S_{ad}$ be the space defined by

$$ S_{ad} = \{(w, \mathbf{M}) \in W \times \mathbf{M} / (w, \mathbf{M}) \text{ satisfy (6) and (7)}\} \quad (8) $$

where $\mathbf{W}$ and $\mathbf{M}$ are the sets of the finite-energy displacements and moments. The problem of Fig. 1 can be expressed as: find $(w, \mathbf{M}) \in S_{ad}$ and $(w', \mathbf{M}') \in S'_{ad}$ which satisfy the boundary conditions:

$$ w = w_d \quad (\text{resp. } w' = w'_d) \quad \text{along } \partial w_d \Omega \quad (\text{resp. } \partial w'_d \Omega') $$

$$ w_{,n} = w_{nd} \quad (\text{resp. } w'_{,n} = w'_{nd}) \quad \text{along } \partial w_{nd} \Omega \quad (\text{resp. } \partial w'_{nd} \Omega') $$

$$ \mathbf{nMn} = M_d \quad (\text{resp. } \mathbf{n'M'n'} = M'_d) \quad \text{along } \partial M_d \Omega \quad (\text{resp. } \partial M'_d \Omega') $$

$$ K_n = \mathbf{n} \operatorname{div}(\mathbf{M}) + (\mathbf{nMn})_{,t} = K_d \quad \text{along } \partial K_d \Omega $$

$$ (\text{resp. } K'_n = K'_d) \quad (\text{resp. } \partial K_d \Omega') $$

$$ [[ \mathbf{tMn} ]] = 0 \quad (\text{resp. } [[ \mathbf{t'M'n'} ]] = 0) \quad \text{at the corner of } \partial\Omega \quad (\text{resp. } (\partial\Omega') $$

$$ w = w' \quad \text{along } \partial\Omega \cap \partial\Omega' $$

$$ w_{,n} = w'_{,n} \quad \text{along } \partial\Omega \cap \partial\Omega' $$

$$ \mathbf{nMn} + \mathbf{n'M'n'} = 0 \quad \text{along } \partial\Omega \cap \partial\Omega' $$

$$ K_n + K'_n = 0 \quad \text{along } \partial\Omega \cap \partial\Omega' \quad (9) $$

Problem (9) can also be written as: find $(w, \mathbf{M}) \in S_{ad}$ and $(w', \mathbf{M}') \in S'_{ad}$ such that

$$ \begin{aligned} \operatorname{Re} \left\{ -i\omega \left( \int_{\partial w_{nd} \Omega} \mathbf{n} \delta \mathbf{Mn} (w_{,n} - w_{nd})^* \, d\partial\Omega - \int_{\partial w_d \Omega} \delta K_n (w - w_d)^* \, d\partial\Omega \right. \right. \\ \left. + \int_{\partial M_d \Omega} (\mathbf{nMn} - M_d) \delta w_{,n}^* \, d\partial\Omega - \int_{\partial K_d \Omega} (K_n - K_d) \delta w^* \, d\partial\Omega \right. \\ \left. + \int_{\partial w_{nd} \Omega'} \mathbf{n}' \delta \mathbf{M'n'} (w'_{,n} - w'_{nd})^* \, d\partial\Omega' - \int_{\partial w_d \Omega'} \delta K'_n (w' - w'_d)^* \, d\partial\Omega' \right. \\ \left. + \int_{\partial M_d \Omega'} (\mathbf{n'M'n'} - M'_d) \delta w'_{,n'}^* \, d\partial\Omega' - \int_{\partial K_d \Omega'} (K'_n - K'_d) \delta w^* \, d\partial\Omega' \right. \\ \left. - \sum_{\partial \Omega \text{ corners}} [[ \mathbf{tMn} ]] \delta w^* - \sum_{\partial \Omega \text{ corners}} [[ \mathbf{t'M'n'} ]] \delta w^* \right. \\ \left. + \int_{\partial \Omega \cap \partial \Omega'} \frac{1}{2} [(\mathbf{nMn} - \mathbf{n'M'n'})(w_{,n} - w'_{,n})^* - \delta(K_n - K'_n)(w - w')^*] \, d\partial\Omega \right. \\ \left. + \int_{\partial \Omega \cap \partial \Omega'} \frac{1}{2} [(\mathbf{nMn} + \mathbf{n'M'n'})(\delta w_{,n} + \delta w'_{,n})^* - (K_n + K'_n)(\delta w + \delta w')^*] \, d\partial\Omega \right) \end{aligned} \quad (10) $$

where $\bullet^*$ and $\operatorname{Re}(\bullet)$ denote the conjugate part and the real part of the complex quantity $\bullet$. Eqs. (8) and (9) correspond to (1) and (2) for an assembly of two plates, and Formulation (10) corresponds to (3). Complete details about Formulation (10), especially its equivalence with (8) and (9), can be found in Ref. [29].

## 3.2. Approximate solution using plane waves

Among the several possible ways to generate the space $S_{ad}^a$ in which to seek an approximate solution (see (4)), the expansion proposed in Ref. [29] seems promising. One assumes that the solution is correctly represented by a specific set of propagative and evanescent waves. Thus, $S_{ad}^a$ is generated by a set of functions of the form $e^{ik \cdot x}$, where $i$ is the complex
<figure><img src="image_3.png" /><figcaption>Fig. 2. Interior, edge and corner waves which satisfy (11) for a homogeneous rectangular plate. The nature of the wave depends on the real and imaginary parts of the coordinates of **k** in function e**ik**x.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Definition of the example problem of Section 3.3: an assembly of 12 plates subjected to a shear force **F** distributed linearly along an edge.</figcaption></figure>

imaginary unit, **k** is the wave vector and **x** is the spatial location. This type of function satisfies (6) and (7) if, and only if

$$ \mathbf{k}^4 = \frac{12(1-v^2)\rho\omega^2}{1+i\eta} \quad (11) $$

See Ref. [29] for more details. Depending on the real and imaginary parts of the coordinates of **k**, different types of waves satisfy (11): interior waves, edge waves or corner waves (see [29]). Fig. 2 gives a representation of such waves for a homogeneous rectangular plate. Then, $S_{ad}^a$ is spanned by the entire set of the interior, edge and corner waves

$$ S_{ad}^a = \operatorname{span}\{e^{i\mathbf{k}_{\text{int},j}\cdot\mathbf{x}}, e^{i\mathbf{k}_{\text{edge},h}\cdot\mathbf{x}}, e^{i\mathbf{k}_{\text{corner},l}\cdot\mathbf{x}} \\ j \in [1; n_{\text{int}}], h \in [1; n_{\text{edge}}], l \in [1; n_{\text{corner}}]\} \quad (12) $$

where $n_{\text{int}}$, $n_{\text{edge}}$ and $n_{\text{corner}}$ are the numbers of waves of each type which are being considered in the approximate solution. The choice of these numbers depends on the problem and on the level of accuracy desired.

### 3.3. Numerical example

Fig. 3 shows the geometry and the mechanical properties of the structure being studied, which is a typical stringer as can be found in a car chassis. The reference solution was calculated using the FEM program MSC Nastran with about seven CQUA4 elements per wavelength. In each plate, the VTCR solution was calculated with 32 interior waves (distributed regularly in all the possible 2-D propagation directions of the reference surface), nine edge waves per edge and no corner waves. The results obtained with the two methods are shown in Fig. 4.

One can see that the VTCR solution and the FEM solution are very similar. The spatial distributions and the amplitudes of the vibrational behavior are nearly the same. Moreover, the maximum relative differences between the VTCR and FEM effective displacements (1/mes($\Omega_i$)) $\int_{\Omega_i} |w|$ d$\Omega$ in each plate $\Omega_i$ ($i=1..12$) are 0.1, 0.03 and 0.07 at 400 Hz, 800 Hz and 1500 Hz respectively. However, significantly fewer degrees of freedom were necessary with the VTCR, especially at higher frequencies (about two hundred times fewer at 1500 Hz). This example clearly demonstrates the benefits of the VTCR in dealing with 3D plate assemblies.

## 4. 3D shell assemblies

### 4.1. Problem description

Again, in order to simplify the presentation, the problem will be formulated for an assembly of two shells, but could be easily generalized to a larger number. Let us consider two isotropic homogeneous shells whose reference surfaces are $\Omega$ and $\Omega'$ with boundaries $\partial\Omega$ and $\partial\Omega'$ respectively. In each shell, the quantities of interest are the displacement vector
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Comparison of the normal displacements (in m) obtained using the VTCTR (left) and the FEM program MSC Nastran (right) for the problem of Fig. 3. The numbers of DOFs used and the circular frequencies considered are given next to each plot.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Problem definition of the 3D assembly of two shells.</figcaption></figure>

**u** = (**v**, **w**) (**v** and **w** being the tangential and normal displacements), the moments (associated with operator **M**) and the resultants (associated with operator **N**). Both shells are assumed to be slightly curved. The action of the environment on **Ω** consists of a prescribed displacement **u**<sub>d</sub> along ∂**u**<sub>d</sub>**Ω** and a force density **F**<sub>d</sub> along ∂**F**<sub>d</sub>**Ω** (see Fig. 5). Similar quantities are defined for **Ω'**.

The formulation uses Koiter's linear shell theory (see Refs. [38,39]). The displacement class is restricted to **u** = **v**+**we**<sub>3</sub>+**z**β with **β** = −∇(**w**)-**Bv**, **z** being the thickness, **e**<sub>3</sub> the normal to the reference surface and **B** the curvature tensor (Kirchhoff's kinematic assumption). The transverse strain energy is neglected.
For structure $\Omega$, let us define the space $S_{ad} = \{(v, w, N, M) \in V \times W \times N \times M\}$ ($V, W, N$ and $M$ being the finite-energy sets of the corresponding physical quantities) such that

$$
\left.
\begin{aligned}
\mathbf{div} \, \mathbf{N} - \mathbf{B}(\mathbf{div} \, \mathbf{M}) &= -\rho \omega^2 h \mathbf{v} && \text{in } \Omega \\
\mathbf{div}(\mathbf{div} \, \mathbf{M}) + \text{Tr}(\mathbf{N}\mathbf{B}) &= -\rho \omega^2 h w && \text{in } \Omega \\
\mathbf{M} &= \frac{h^3}{12}(1+i\eta)\mathbf{K}_{CP}\mathbf{X}(\mathbf{u}) && \text{in } \Omega \\
\mathbf{N} &= h(1+i\eta)\mathbf{K}_{CP}\gamma(\mathbf{u}) && \text{in } \Omega \\
\mathbf{X}(\mathbf{u}) &= \epsilon(\beta) - [\mathbf{B}\epsilon(\mathbf{v} + w\mathbf{e}_3)]_{\text{sym}} \\
\gamma(\mathbf{u}) &= \epsilon(\mathbf{v} + w\mathbf{e}_3)
\end{aligned}
\right\} \quad (13)
$$

where $K_{CP}$ is Hooke's plane stress operator, $\epsilon$ is the gradient operator, $\rho$ is the density, $\eta$ is the damping coefficient and $h$ is the thickness of the shell. $S'_{ad}$ is defined in the same way. Thus, the problem to be solved can be expressed as: find $(v, w, N, M) \in S_{ad}$ and $(v', w', N', M') \in S'_{ad}$ which verify the boundary conditions

$$
\left.
\begin{aligned}
\mathbf{v} &= \mathbf{v}_d \quad \text{along } \partial\mathbf{v}_d \Omega \\
w &= w_d \quad \text{along } \partial w_d \Omega \\
w_n &= w_{nd} \quad \text{along } \partial w_{nd} \Omega \\
\mathbf{N}_n &= \mathbf{N}_n - \mathbf{B}\mathbf{M}\mathbf{n} = \mathbf{N}_d \\
&\quad \text{along } \partial \mathbf{N}_d \Omega \\
K_n &= (\mathbf{div} \, \mathbf{M})\mathbf{n} + (\mathbf{n}\mathbf{M}\mathbf{t})_{,t} = K_d \\
&\quad \text{in } \partial K_d \Omega \\
M_n &= \mathbf{n}\mathbf{M}\mathbf{n} = M_d \quad \text{in } \partial M_d \Omega \\
[\mathbf{n}\mathbf{M}\mathbf{t}] &= 0 \\
&\quad \text{at the corners of } \partial\Omega \\
\mathbf{v} &= \mathbf{v}' \quad \text{along } \partial\Omega \cap \partial\Omega' \\
w &= w' \quad \text{along } \partial\Omega \cap \partial\Omega' \\
w'_n &= w'_n \quad \text{along } \partial\Omega \cap \partial\Omega' \\
\mathbf{N}_n - \mathbf{B}\mathbf{M}\mathbf{n} &= \mathbf{N}'\mathbf{n}' - \mathbf{B}'\mathbf{M}'\mathbf{n}' \quad \text{along } \partial\Omega \cap \partial\Omega' \\
(\mathbf{div} \, \mathbf{M})\mathbf{n} + (\mathbf{n}\mathbf{M}\mathbf{t})_{,t} \\
&= (\mathbf{div} \, \mathbf{M}')\mathbf{n}' + (\mathbf{n}'\mathbf{M}'\mathbf{t}')_{,t} \quad \text{along } \partial\Omega \cap \partial\Omega' \\
\mathbf{n}\mathbf{M}\mathbf{n} &= \mathbf{n}'\mathbf{M}'\mathbf{n}' \quad \text{along } \partial\Omega \cap \partial\Omega'
\end{aligned}
\right\} \quad (14)
$$

Problem (14) can also be written as: find $(v, w, N, M) \in S_{ad}$ and $(v', w', N', M') \in S'_{ad}$ such that

$$
\begin{aligned}
& \text{Re} \left\{ -i\omega - \int_{\partial w_d \Omega} \delta \mathbf{N}_n \cdot (\mathbf{v} - \mathbf{v}_d)^* \, d\partial\Omega - \int_{\partial \mathbf{N}_d \Omega} (\mathbf{N}_n - \mathbf{N}_d) \cdot \delta \mathbf{v}^* \, d\partial\Omega \right. \\
& \qquad \left. + \int_{\partial w_{nd} \Omega} \delta M_n \cdot (w_n - w_{nd})^* \, d\partial\Omega + \int_{\partial \mathbf{M}_d \Omega} (M_n - M_d) \cdot \delta w_n^* \, d\partial\Omega \right. \\
& \qquad \left. - \int_{\partial w_d \Omega} \delta K_n (w - w_d)^* \, d\partial\Omega - \int_{\partial K_d \Omega} (K_n - K_d) \delta w^* \, d\partial\Omega \right. \\
& \qquad \left. - \int_{\partial w'_d \Omega} \delta \mathbf{N}'_n \cdot (\mathbf{v}' - \mathbf{v}_d)^* \, d\partial\Omega' - \int_{\partial \mathbf{N}'_d \Omega} (\mathbf{N}'_n - \mathbf{N}'_d) \cdot \delta \mathbf{v}'^* \, d\partial\Omega' \right. \\
& \qquad \left. + \int_{\partial w_{nd}' \Omega} \delta M_n' \cdot (w'_n - w_{nd})^* \, d\partial\Omega' + \int_{\partial \mathbf{M}'_d \Omega} (M_n' - M_d') \cdot \delta w_n'^* \, d\partial\Omega' \right. \\
& \qquad \left. - \int_{\partial w'_d \Omega} \delta K_n' (w' - w_d)^* \, d\partial\Omega' - \int_{\partial K'_d \Omega} (K_n' - K_d') \delta w'^* \, d\partial\Omega' \right. \\
& \qquad \left. - \sum_{\partial \Omega \text{ corners}} [\mathbf{n}\mathbf{M}\mathbf{t}] \delta w^* - \sum_{\partial \Omega' \text{ corners}} [\mathbf{n}'\mathbf{M}'\mathbf{t}'] \delta w'^* \right. \\
& \qquad \left. + \int_{\Gamma} \frac{1}{2} [-(\delta \mathbf{N}_n - \delta \mathbf{N}'_n)(\mathbf{v} - \mathbf{v}')^* + (\delta M_n - \delta M_n')(w_n - w_n')^* \right. \\
& \qquad \left. -(\delta K_n - \delta K_n')(w - w')^*] d\partial\Omega + \int_{\Gamma} \frac{1}{2} [-(\mathbf{N}_n + \mathbf{N}'_n)(\delta \mathbf{v} + \delta \mathbf{v}')^* \right. \\
& \qquad \left. +(M_n + M_n')(\delta w_n + \delta w_n')^* - (K_n + K_n')(\delta w + \delta w')^*] d\partial\Omega \right\} \\
&= 0 \quad \forall (\delta \mathbf{v}, \delta w, \delta \mathbf{N}, \delta \mathbf{M}) \in S_{ad} \text{ and } (\delta \mathbf{v}', \delta w', \delta \mathbf{N}', \delta \mathbf{M}') \in S'_{ad}
\end{aligned}
\quad (15)
$$
<figure><img src="image_7.png" /><figcaption>Fig. 6. Interior, edge and corner waves which satisfy (16) for a homogeneous cylinder. The nature of the wave depends on the real and imaginary parts of the coordinates of $k$ in function $e^{ik \cdot x}$.</figcaption></figure>

Eqs. (13) and (14) correspond to (1) and (2) for an assembly of two shells, and Formulation (15) corresponds to (3). Complete details about Formulation (15), especially its equivalence with (13) and (14), can be found in Ref. [32].

## 4.2. Approximate solution using plane waves

Again, among the several possible ways to generate the space $S_{ad}^a$ in which to seek an approximate solution (see (4)), the expansion proposed in Ref. [32] seems promising. One assumes that $S_{ad}^a$ is correctly represented by a specific set of propagative and evanescent waves expressed as $e^{ik \cdot x}$. In order to find the equation that the wave vector $k$ must verify, one can perform an asymptotic expansion of (13) with respect to the small parameter $\sqrt{\hbar/R}$, where $R$ is the smallest radius of curvature of the shell being considered (which is assumed to be only slightly curved). Then, one can easily prove that the wave vector must verify

$$ \mathbf{k}^8 = \frac{1}{1+i\eta} \frac{12(1-v^2)\rho\omega^2}{Eh^2} \mathbf{k}^4 - \frac{12(1-v^2)}{h^2} (\mathbf{k}R\mathbf{B}\mathbf{R}\mathbf{k})^2 \quad (16) $$

where $R$ is a rotation matrix in the coordinate system of the shell's reference surface. Of course, in the case of a plate, i.e. if $B=0$, (16) boils down to (11). Otherwise, the additional portion of Eq. (16) is necessary in order to calculate waves in shells because it is used to represent some well-known characteristics of propagative waves in certain directions at certain frequencies (see [32] for details).

Again, (16) can be satisfied by interior waves, edge waves or corner waves depending on the real and imaginary parts of the coordinates of $k$. Fig. 6 illustrates these different types of waves in the case of a homogeneous cylinder. Then, $S_{ad}^a$ is spanned by the entire set of the interior, edge and corner waves

$$ S_{ad}^a = \operatorname{span}\{e^{ik_{\text{int},j} \cdot x}, e^{ik_{\text{edge},h} \cdot x}, e^{ik_{\text{corner},l} \cdot x} \\ j \in [1; n_{\text{int}}], h \in [1; n_{\text{edge}}], l \in [1; n_{\text{corner}}]\} \quad (17) $$

where $n_{\text{int}}$, $n_{\text{edge}}$ and $n_{\text{corner}}$ are the numbers of waves of each type which are being considered in the approximate solution. Again, the choice of these numbers depends on the problem and on the level of accuracy desired.

## 4.3. Numerical example

Fig. 7 shows the geometry of a 3D assembly of two cylindrical shells and one plate. The mechanical properties are $E = 75$ GPa, $\eta = 0.0001$, $v = 0.3$, $\rho = 2750$ kg/m$^3$, $\omega = 2\pi \cdot 1000$ rad/s, $F_d = 1$ N/m. The cylinder's height $L_z$ and its curvature are both equal to 1 m and its thickness is $h = 0.01$ m.

We used the solution obtained with the FEM program MSC Nastran as the reference solution. The mesh seed used in the program was set to create 10 QUADR shell elements per wavelength and the solution was obtained with about 1,200,000 DOFs. For the VTCR solution, the structure was divided into three parts (shell $\Omega_1$, shell $\Omega_3$, and plate $\Omega_2$, see Fig. 7). The model involved 264 DOFs (corresponding in each substructure to 24 interior waves, five edge waves per edge and no corner waves). One can verify in Fig. 8 that the result given by the VTCR is quite good and very similar to that given by the FEM regarding the distribution and amplitudes of the displacements).

In order to compare the solutions in terms of an effective quantity, we calculated the effective displacement $(1/\text{mes})(\Omega_i) \int_{\Omega_i} \|v\| d\Omega$ in each substructure $\Omega_i$, $i=1..3$. The comparison is shown in Table 1. One can see that the VTCR gave accurate results at a very low cost. The VTCR can calculate any number of homogeneous substructures (plates, shells, etc.) without difficulty.
<figure><img src="image_8.png" /><figcaption>Fig. 7. Definition of the example problem of Section 4.3: an assembly of two cylindrical shells and one plate.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Comparison of the displacements (in m) obtained using the VTCR (left with 264 DOFs) and the FEM program MSC Nastran (right with about 1,200,000 DOFs) for the problem of Section 4.3 illustrated in Fig. 7.</figcaption></figure>

<caption>**Table 1**</caption>

Comparison of the effective displacements in substructures $\Omega_i$, $i=1..3$ of the problem of Fig. 7, whose solutions are given in Fig. 8.

<table><thead><tr><th>DOFs</th><th>VTCR<br>264</th><th>FEM (MSC Nastran)<br>1,200,000</th></tr></thead><tbody><tr><td>Effective displacement (m)</td><td></td><td></td></tr><tr><td>Shell Ω<sub>1</sub></td><td>1.78 × 10<sup>-8</sup></td><td>1.62 × 10<sup>-8</sup></td></tr><tr><td>Plate Ω<sub>2</sub></td><td>2.57 × 10<sup>-8</sup></td><td>2.52 × 10<sup>-8</sup></td></tr><tr><td>Shell Ω<sub>3</sub></td><td>1.01 × 10<sup>-8</sup></td><td>1.09 × 10<sup>-8</sup></td></tr></tbody></table>

# 5. 3D acoustic problems

## 5.1. Problem description

Let us consider a general steady-state interior dynamic problem in a 3D acoustic cavity $\Omega$ filled with a fluid characterized by its sound velocity $c$, its density $\rho$ and its damping coefficient $\eta$. The problem of the steady-state dynamic
<figure><img src="image_10.png" /><figcaption>Fig. 9. The general 3D acoustic problem.</figcaption></figure>

behavior of $\Omega$ at a fixed circular frequency $\omega$ is: find the pressure $p \in H^1(\Omega)$ such that

$$ egin{cases} 
\Delta p + k^2 p = 0 & \text{in } \Omega \\ 
p - Z \cdot L_v(p) = h_d & \text{over } \partial Z \\ 
p = p_d & \text{over } \partial p \\ 
L_v(p) = v_d & \text{over } \partial v 
\end{cases} 
\quad (18) $$

where $k = \omega/c$ is the wavenumber, $Z$ is an impedance coefficient and $L_v(\square)$ an operator defined as $L_v(\square) = (i/\rho\omega)(\partial\square/\partial\mathbf{n}) = (i/\rho\omega)\mathbf{n} \cdot \nabla(\square)$, $\mathbf{n}$ being the outward normal to $\partial\Omega$, $h_d$, $p_d$ and $v_d$ denote respectively a prescribed excitation over $\partial Z\Omega$, a prescribed pressure over $\partial p\Omega$ and a prescribed velocity over $\partial v\Omega$.

Let us partition $\Omega$ into $n_{el}$ non-overlapping subcavities $\Omega_E$ and use the notations $\Gamma_E = \partial\Omega_E$ and $\Gamma_{E,E'} = \partial\Omega_E \cap \partial\Omega_{E'}$. The continuity conditions over $\Gamma_{E,E'}$ are

$$ \left| \begin{array}{l} p_E = p_E \\ L_v(p_E) = -L_v(p_E) \end{array} \right. \qquad (19) $$

The problem defined by (18) and (19) is illustrated in Fig. 9.

The VTCR formulation is obtained by rewriting boundary value problem (18) and (19) in a weak form. Let $S_{ad}^E$ denote the space of the functions which satisfy the Helmholtz equation (the first equation of Problem (18)) everywhere in subcavity $\Omega_E$:

$$ S_{ad}^E = \{p_E \in H^1(\Omega_E) | \Delta p_E + k^2 p_E = 0\} \quad (20) $$

The VTCR formulation of Problem (18) and (19) is

Find $((p_1, p_2, \dots, p_{n_{el}}) \in S_{ad}^1 \times S_{ad}^E \times \dots \times S_{ad}^{n_{el}})$ such that

$$ \begin{aligned} \operatorname{Re} \left\{ \frac{1}{2} \sum_{E=1}^{n_{el}} \int_{\partial Z \Omega_E} (p_E - Z L_v(p_E) - h_{dE}) \cdot L_v(\delta p_E)^* + \left(\frac{p_E}{Z} - L_v(p_E) - \frac{h_{dE}}{Z}\right)^* \cdot \delta p_E \, d\partial\Omega_E \right. \\ + \sum_{E=1}^{n_{el}} \int_{\partial p \Omega_E} (p_E - p_{dE}) \cdot L_v(\delta p_E)^* \, d\partial\Omega_E \\ + \sum_{E=1}^{n_{el}} \int_{\partial v \Omega_E} (L_v(p_E) - v_{dE})^* \cdot \delta p_E \, d\partial\Omega_E \\ + \sum_{E,E' < E} \frac{1}{2} \int_{\Gamma_{E,E'}} \left( (p_E - p_{E'}) \cdot L_v(\delta p_E - \delta p_{E'})^* + L_v(p_E + p_{E'})^* \cdot (\delta p_E + \delta p_{E'}) \right) d\partial\Omega_E \right\} = 0 \\ \forall \{ (\delta p_1, \delta p_2, \dots, \delta p_{n_{el}}) \in S_{ad}^1 \times S_{ad}^E \times \dots \times S_{ad}^{n_{el}} \} \end{aligned} \quad (21) $$

Eqs. (18) and (19) correspond to (1) and (2) for acoustic problems, and Formulation (21) corresponds to (3). Complete details about Formulation (21), especially its equivalence with (18) and (19), can be found in Ref. [37].

## 5.2. Approximate solution using plane waves

Again, there are several possible ways to generate the space $S_{ad}^{E,a}$ in which to seek an approximate solution (see (4)). In the cases of plates and shells of Sections 3 and 4, although the structure was defined in 3D, the waves propagated in 2D, i.e. in the plate's or the shell's reference surface. Here, the difficulty comes from the fact that waves can propagate in all 3D directions. Taking into account all possible directions can lead to a huge problem. Therefore, the discretization of these directions must be carried out efficiently.

In order to do that, a direct and straightforward extension of what was proposed for 2D acoustics in Ref. [37] seems promising. The pressure field is sought as an integral distribution of plane waves, also known as Herglotz wave functions,
<figure><img src="image_11.png" /><figcaption>Fig. 10. Illustration of the shape functions $\Phi_l^m(\mathbf{x})$ used in (23) to define $S_{ad}^{E,a}$. For each index $l$ and each momentum $m$, the real and imaginary parts are plotted on top and bottom.</figcaption></figure>

which satisfy the homogeneous Helmholtz equation. In the 3D case, this plane wave distribution is defined over the unit sphere:

$$p_E(\mathbf{x}) = \int_{\theta = -\pi}^{\pi} \int_{\varphi = 0}^{\pi} A_E(\theta, \varphi) \cdot e^{i\mathbf{k}(\theta, \varphi) \cdot \mathbf{x}} \, d\theta \sin \varphi \, d\varphi \quad (22)$$

where $A_E$ denotes the amplitude of the plane waves in the 3D spherical direction $(\theta, \varphi)$. One can easily prove that (22) satisfies (20) and, therefore, that $p_E(\mathbf{x})$ belongs to $S_{ad}^E$.

In order to preserve all possible wave directions, avoid too refined a quadrature of the unit sphere and benefit from the advantages of the Fourier decomposition introduced in Ref. [37], function $A_E(\theta, \varphi)$ is described using a truncated Laplace series, which is a 3D extension of a Fourier series. Then, the discrete space $S_{ad}^{E,a}$ of (21) can be expressed as

$$S_{ad}^{E,a} = \operatorname{span}\left\{ \int_{\theta = -\pi}^{\pi} \int_{\varphi = 0}^{\pi} Y_l^m(\theta, \varphi) e^{i\mathbf{k}(\theta, \varphi) \cdot \mathbf{x}} \, d\theta \sin \varphi \, d\varphi, m = -l, \dots, l, l = 0, \dots, L_E \right\} \\ = \operatorname{span}\{\Phi_l^m(\mathbf{x}), m = -l, \dots, l, l = 0, \dots, L_E\} \quad (23)$$

where $\mathbf{k}(\theta, \varphi)$ is the wave vector of the plane wave propagating in the spherical direction $(\theta, \varphi)$ and $Y_l^m(\theta, \varphi) = \sqrt{2} \cdot (l-m)!/(l+m)! \cdot P_l^m(\cos \theta) \cdot e^{im\varphi}$ is the spherical harmonic of nonnegative index $l$ and momentum $m$, and $L_E$ the maximum index retained in the discrete space $S_{ad}^{E,a}$ of dimension $(L_E+1)^2$. In the last expression, $P_l^m(X) = ((-1)^m/2^l \cdot l!) (1-X^2)^{m/2} (\hat{\chi}^{m+l}(\chi^2-1)^{l/2} \hat{\chi}^{m+l})$ denotes the Legendre polynomial and parameter $m$ varies from $-l$ to $l$. The set $\{Y_l^m(\theta, \varphi)\}_{|m| \le l \le \infty}$ constitutes a complete orthogonal coordinate system for the unit sphere with regard to the classical $L^2$ scalar product.

Fig. 10 illustrates the corresponding shape functions $\Phi_l^m(\mathbf{x})$ in a 3D cubic zone. The choice of the number of these shape functions depends on the problem and on the level of accuracy desired.

## 5.3. Numerical example

The example considered (see Fig. 11) is a Z-shaped cavity filled with air ($\rho = 1.25 \text{ kg/m}^3$, $c = 340 \text{ m/s}$, $\eta = 0.0001$). The boundary conditions consisted of a unit external pressure, a prescribed impedance $Z = 425 \text{ Pa s/m}$ or a rigid-wall condition as indicated in the figure. The circular frequencies considered were $\omega = 2\pi \cdot 500 \text{ Hz}$, $1000 \text{ Hz}$ and $1500 \text{ Hz}$. For each frequency, a reference solution (in the form of a pressure distribution $p_{ref}$) was calculated using the FEM program COMSOL. In order to achieve an appropriate quality of results, the mesh seed used in the program was set to create 14 3D TETRA elements per wavelength, leading to three models with 7, 160, 57, 230 and 192, 420 DOFs respectively for the three frequencies.

For the VTCR solution, the structure was divided into three subcavities $\Omega_1, \Omega_2$ and $\Omega_3$ (see Fig. 11). Fig. 12 shows for each frequency the real part of the VTCR pressure, the number of DOFs used and the corresponding relative error $\int_{\Omega} |p_{VTCR} - p_{ref}| \, d\Omega / \int_{\Omega} |p_{ref}| \, d\Omega$ where $p_{ref}$ is the solution given by the FEM. One can see that in all three cases the solutions
<figure><img src="image_12.png" /><figcaption>Fig. 11. Problem definition for the Z-shaped cavity filled with air of Section 5.3.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. The real part of the VTCR pressure for the problem of Section 5.3 illustrated in Fig. 11. The circular frequency considered, the number of DOFs and the relative error are indicated in each case.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Comparison of the FRFs given by the FEM and the VTCR for the problem of Section 5.3 illustrated in Fig. 11. The corresponding numbers of DOFs are indicated on the vertical scale on the right. For the FEM, 14 elements per wavelength were used. For the VTCR, the number of shape functions was chosen so that the relative error compared to the FEM would be less than $10^{-3}$ for each frequency.</figcaption></figure>

given by the VTCR and the FEM are very similar, but with the VTCR the solution was obtained using far fewer DOFs (over 100 times fewer at $\omega = 2\pi \cdot 150$ Hz).

Fig. 13 shows another interesting comparison concerning the global Frequency Response Functions (FRFs) in energy for the FEM and the VTCR. In addition, the number of DOFs used is indicated on the vertical scale on the right. Again, in the case of the FEM, 14 elements per wavelength were used in order to ensure the accuracy of the solution (especially for the 3D calculations). In the case of the VTCR, the number of shape functions was chosen so that the relative error compared to the FEM would be less than $10^{-3}$ for each frequency. One can see that the two FRFs thus predicted are very similar, which is normal since the number of DOFs for the VTCR was chosen specifically so the error between the two curves would be very small. However, as mentioned before, the number of DOFs needed for the VTCR was considerably smaller than for the FEM. Furthermore, it is clear that the need for large numbers of DOFs increases faster for the FEM than for the VTCR. Again, this leads to the conclusion that the VTCR is a competitive numerical technique for the analysis of 3D acoustic problems.
# 6. Conclusion

This paper presents an overview of the development of the VTCR for the resolution of mid-frequency 3D engineering problems. First, the paper reviews the basic concepts and main characteristics of the VTCR, a new variational formulation in which approximations within subregions are defined a priori independently of one another, and a two-scale approximation with a strong mechanical content is used to derive the approximation spaces. Following that introduction, the details of the resolution of problems involving 3D assemblies of plates or shells as well as acoustic problems are given.

For all these problems, the discussion of the VTCR approach follows the same outline: the problem is defined, the expansion of the field variables is presented and a numerical example is given. Such a systematic presentation was possible because the VTCR is a general-purpose predictive numerical technique for mid-frequency problems. Throughout these numerical examples, we show that the VTCR is capable of finding an accurate solution with great numerical benefits compared to classical deterministic numerical techniques, such as the FEM.

Some possible subjects for additional work are the investigation of vibroacoustic problems and the development of an efficient and robust numerical strategy for large calculations in a frequency range.

# Acknowledgments

The authors gratefully acknowledge the support of ITN Marie Curie project GA-214909 "MID-FREQUENCY - CAE Methodologies for Mid-Frequency Analysis in Vibration and Acoustics".

# References

[1] I. Babuska, F. Ihlenburg, E.T. Paik, S.A. Sauter, A generalized finite element method for solving the Helmholtz equation in two dimensions with minimal pollution, *Computer Methods in Applied Mechanics and Engineering* 128 (1995) 325-359.

[2] L.P. Franca, C. Farhat, A.P. Macedo, M. Lesoinne, Residual-free bubbles for the Helmholtz equation, *International Journal for Numerical Methods in Engineering* 40 (1997) 4003-4009.

[3] I. Harari, T.J.R. Hughes, Galerkin/least-squares finite element methods for the reduced wave equation with non-reflecting boundary conditions in unbounded domains, *Computer Methods in Applied Mechanics and Engineering* 98 (3) (1992) 411-454.

[4] T.J.R. Hughes, Multiscale phenomena: Greens functions and the Dirichlet-to-Neumann formulation and subgrid scale models and bubbles and the origins of stabilized methods, *Computer Methods in Applied Mechanics and Engineering* 127 (1995) 387-401.

[5] J.M. Melenk, I. Babuska, The partition of unity finite element method: basic theory and applications, *Computer Methods in Applied Mechanics and Engineering* 139 (1996) 289-314.

[6] C. Soize, Reduced models in the medium frequency range for the general dissipative structural dynamic systems, *European Journal of Mechanics and A/Solids* 17 (1998) 657-685.

[7] T. Strouboulis, K. Copps, I. Babuska, The generalized finite element method: an example of its implementation and illustration of its performance, *International Journal for Numerical Methods in Engineering* 47 (2000) 1401-1417.

[8] L.F. Greengard, V. Rokhlin, A fast algorithm for particle simulations, *Journal of Computational Physics* 73 (1987) 325-348.

[9] M. Bonnet, S. Chaillat, J.-F. Semblat, A multi-level fast multipole BEM for 3-D elastodynamics in the frequency domain, *Computer Methods in Applied Mechanics and Engineering* 197 (2008) 4233-4249.

[10] R.S. Langley, V. Cotoni, Response variance prediction in the statistical energy analysis of built-up systems, *Journal of the Acoustical Society of America* 115 (2004) 706-718.

[11] R.H. Lyon, G. Maidanik, Power flow between linearly coupled oscillators, *Journal of the Acoustical Society of America* 34 (1962) 623-639.

[12] B.R. Mace, On the statistical energy analysis hypothesis of coupling power proportionality and some implications of its failure, *Journal of Sound and Vibration* 178 (1994) 95-112.

[13] M.N. Ichchou, A. Le Bot, L. Jézéquel, Energy model of one-dimensional, multipropulsive systems, *Journal of Sound and Vibration* 201 (1997) 534-554.

[14] R.S. Langley, A wave intensity technique for the analysis of high frequency vibrations, *Journal of Sound and Vibration* 159 (1992) 483-502.

[15] B. Mace, Statistical energy analysis: coupling loss factor indirect coupling and system modes, *Journal of Sound and Vibration* 279 (2005) 141-170.

[16] L. Maxit, J.L. Guyader, Estimation of SEA coupling loss factors using a dual formulation and FEM modal information. Part I. Theory, *Journal of Sound and Vibration* 230 (2001) 907-930.

[17] V.D. Belov, S.A. Ryback, Applicability of the transport equation in the one-dimensional wave propagation problem, *Journal of Soviet Physics Acoustics* 21 (1975) 173-180.

[18] E. Trefftz, Ein gegenstück zum ritzschen verfahren (A counterpart to the Ritz's method), in: *Second International Congress on Applied Mechanics*, 1926, pp. 131-137.

[19] O. Laghrouche, P. Bettess, Short wave modelling using special finite elements, *Journal of Computational Acoustics* 8 (2000) 189-210.

[20] T. Strouboulis, I. Babuska, R. Hidajat, The generalized finite element method for Helmholtz equation: theory and computation and open problems, *Computer Methods in Applied Mechanics and Engineering* 195 (2006) 4711-4731.

[21] T. Strouboulis, R. Hidajat, Partition of unity method for Helmholtz equation: q-convergence for plane-wave and wave-band local bases, *Applications of Mathematics* 51 (2006) 181-204.

[22] O. Cessenat, B. Despres, Application of an ultra weak variational formulation of elliptic PDES to the two-dimensional Helmholtz problem, *SIAM Journal on Numerical Analysis* 35 (1998) 255-299.

[23] P. Monk, D.Q. Wang, A least-squares method for the Helmholtz equation, *Computer Methods in Applied Mechanics and Engineering* 175 (1999) 121-136.

[24] C. Farhat, I. Harari, L.P. Franca, The discontinuous enrichment method, *Computer Methods in Applied Mechanics and Engineering* 190 (2001) 6455-6479.

[25] P. Bouillard, S. Suleau, Element-free Galerkin solutions for Helmholtz problems: formulation and numerical assessment of the pollution effect, *Computer Methods in Applied Mechanics and Engineering* 162 (1998) 317-335.

[26] E. Perrey-Debain, J. Trevelyan, P. Bettess, Wave boundary elements: a theoretical overview presenting applications in scattering of short waves, *Engineering Analysis with Boundary Elements* 28 (2004) 131-141.

[27] W. Desmet, P. Sas, D. Vandepitte, An indirect Trefftz method for the steady-state dynamic analysis of coupled vibro-acoustic systems, *Computer Assisted Mechanics and Engineering Sciences* 8 (2001) 271-288.
[28] P. Ladevèze, A new computational approach for structure vibrations in the medium frequency range, *Comptes Rendus de l'Académie des Sciences Paris* 322 (IIb) (1996) 849-856.

[29] P. Ladevèze, L. Arnaud, P. Rouch, C. Blanzé, The variational theory of complex rays for the calculation of medium frequency vibrations, *Engineering Computations* 18 (1-2) (2001) 193-214.

[30] R. Rouch, P. Ladevèze, The variational theory of complex rays: a predictive tool for medium-frequency vibrations, *Computer Methods in Applied Mechanics and Engineering* 192 (2003) 3301-3315.

[31] P. Ladevèze, L. Blanc, P. Rouch, C. Blanzé, A multiscale computational method for medium-frequency vibrations of assemblies of heterogeneous plates, *Computers and Structures* 81 (2003) 1267-1276.

[32] H. Riou, P. Ladevèze, P. Rouch, Extension of the variational theory of complex rays to shells for medium-frequency vibrations, *Journal of Sound and Vibration* 272 (1-2) (2004) 341-360.

[33] P. Ladevèze, P. Rouch, H. Riou, X. Bohineust, Analysis of medium-frequency vibrations in a frequency range, *Journal of Computational Acoustics* 11 (2003) 255-284.

[34] P. Ladevèze, M. Chevreuil, A new computational method for transient dynamics including the low- and the medium-frequency ranges, *International Journal for Numerical Methods in Engineering* 64 (2005) 503-527.

[35] H. Riou, P. Ladevèze, B. Sourcis, The multiscale VTCR approach applied to acoustics problems, *Journal of Computational Acoustics* 16 (2008) 487-505.

[36] H. Riou, P. Ladevèze, B. Sourcis, B. Faverjon, L. Kovalevsky, An adaptive numerical strategy for the medium-frequency analysis of Helmholtz's problem, *Journal of Computational Acoustics*, in press. http://dx.doi.org/10.1142/S0218396X11004481.

[37] L. Kovalevsky, P. Ladevèze, H. Riou, The Fourier version of the Variational Theory of Complex Rays for medium-frequency acoustics, *Computer Assisted Mechanics and Engineering Sciences*, in press. http://dx.doi.org/10.1016/j.cma.2012.03.009.

[38] W.T. Koiter, A consistent first approximation in the general theory of thin elastic shells, *Proceedings of the IUTM Symposium on Theory of Thin Elastic Shells*, Delft, August 1959, Amsterdam, North-Holland, 1960, pp. 12-33.

[39] A.W. Leissa, Vibration of shells, *Acoustical Society of America* (1993). ISBN:1563962934.