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

HAL
open science

A multiscale "Trefftz" computational method for
medium-frequency vibrations of assemblies of
heterogeneous plates with uncertainties

Laurent Blanc, Claude Blanzé, Philippe Rouch

► To cite this version:

Laurent Blanc, Claude Blanzé, Philippe Rouch. A multiscale "Trefftz" computational method for medium-frequency vibrations of assemblies of heterogeneous plates with uncertainties. Computers & Structures, 2007, 85 (10), pp.595-605. 10.1016/j.compstruc.2006.08.078. hal-03177687

HAL Id: hal-03177687

https://hal.science/hal-03177687v1

Submitted on 15 May 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A multiscale “Trefftz” computational method for medium-frequency vibrations of assemblies of heterogeneous plates with uncertainties

L. Blanc *, C. Blanzé, P. Rouch

LMT-Cachan, ENS Cachan/CNRS/Paris 6 University, 61 avenue du Président Wilson, 94235 Cachan Cedex, France

A new approach, called the “Variational Theory of Complex Rays”, is being developed for the calculation of the vibrations of slightly damped elastic plates in the medium-frequency range. In this approach, the resolution of a small system of equations which does not result from a refined spatial discretization of the structure leads to the evaluation of effective quantities (deformation energy, vibration amplitude, ...). This theory has already been validated for assemblies of homogeneous substructures and extended to the case of heterogeneous substructures. In this paper, we propose some improvements to deal with sensitivity analysis. The objective is to evaluate the influence of discontinuities with uncertain positions, sizes or directions on the structural response. A specific partial inversion of the constrained system of equations enables multiple resolutions for combinations of the parameters defining the possible shapes of the discontinuity to be carried out at low cost.

*Keywords: Vibrations; Medium-frequency range; Complex rays; Heterogeneous structures; Uncertainties

## 1. Introduction

The modeling and analysis of the vibrational response of elastic structures are undoubtedly among the key issues which arise in the design of satellites or car chassis. Today, at least concerning modeling and analysis, there remain no major difficulties in the low frequency range, even for complex structures [6]. Concerning high frequencies, computational tools quite distinct from those used for low frequencies do exist, particularly the SEA method [26,27].

Conversely, the modeling and analysis of medium-frequency vibrations, which constitute the topic of this paper, continue to cause certain problems. Difficulties are experienced in extending the SEA method, which is appropriate mostly for high frequencies. In particular, the scale on which space is described is too coarse for medium-

frequency analysis since the spatial aspect disappears almost entirely. The difficulty in attempting to extend low frequency methods to this context is that the lengths of variation of the phenomena being studied are very small compared to the characteristic dimension of the structure. Therefore, the finite element calculations involved would require an unreasonable number of degrees of freedom. Apart from these serious numerical difficulties, the quantities involved in the calculations remain associated with small variation lengths which are not very significant. These are not “effective” quantities, i.e., they are not representative of the vibratory response of the structure on the time and space scales being considered. Therefore, they lead to results which are very sensitive to inaccuracies in the data. Nevertheless, various improved finite element approaches have been studied, in particular in [1,14,5,2, 9,15,16,12,13,34–36,24,25,30].

There have been only few works dedicated specifically to the calculation of medium-frequency vibrations. For rods

<footnote>* Corresponding author. Tel.: +33 1 47 40 22 41; fax: +33 1 47 40 27 85.
E-mail address: blanc@lmt.ens-cachan.fr (L. Blanc).</footnote>
and beams, the problem was solved in particular in [7,11,29]. For more complex structures, such as plates or shells, boundary element methods were proposed in [8] and methods based on the use of specific reduced bases can be found in [28,32,33]. Most of these works are closely related to finite element methods and the quantities calculated are still not effective quantities. Therefore, in our opinion, these are not “true” medium-frequency methods. On the contrary, the theory initiated by Belov and Ryback [3] is built upon the concepts of effective energy density and effective power flow. However, despite the improvements proposed in [17], this theory still faces theoretical obstacles.

The approach followed here is the Variational Theory of Complex Rays (VTCR), a predictive tool designed specifically to deal with medium-frequency problems which was first introduced by Ladevèze [19]. This is a true medium-frequency method in the sense that the calculations are performed on effective quantities. Previous works have already validated this strategy for two-dimensional and three-dimensional assemblies of elastic homogeneous plates with low damping [20–22].

However, most industrial structures present some kind of structural heterogeneity, whether intentional (such as portholes or equipment connections) or unintentional (such as cracks). Their vibratory response in the medium-frequency range is very sensitive to these structural discontinuities. Although the VTCR was originally based on the assumption that the structure is an assembly of homogeneous substructures, an extension to the VTCR enabling the designer to take heterogeneity into account has already been proposed [23]. This paper goes one step further and deals with three-dimensional assemblies of homogeneous and heterogeneous substructures in an uncertain context.

## 2. The basic aspects of the VTCR

### 2.1. Principle

In the general case of a structure composed of several substructures, the first step of the VTCR consists in associating each substructure with a superelement described by degrees of freedom which correspond to local basic modes. These basic modes, which are defined on two scales and verify the dynamic equations (local equilibrium and constitutive relation) exactly, have a strong mechanical meaning and are called complex rays. They can be divided into families related to interior, edge or corner zones. In the vicinity of a point, the solution is assumed to be properly described locally as the superimposition of an infinite number of such local vibration modes. The unknowns are large-wavelength quantities, each representing the amplitude of the basic mode with which it is associated.

Next, an appropriate variational formulation enables the boundary and transmission conditions to be verified in an average sense. Our choice of reasoning on a superelement's level requires the formulation to allow approximations within the substructures which are a priori

independent of one another. These approximations must not necessarily verify the transmission conditions a priori. In practice, the variational formulation associates each superelement with an elementary matrix which represents the interaction among the fields and a right-hand side which represents the interaction between the fields and the boundary conditions. The transmission conditions at the interface are taken into account automatically during the assembly stage.

The last characteristic of the VTCR is that from the calculated discretized amplitudes it retains only effective quantities related to the elastic energy, the kinetic energy, the dissipation work, etc.

### 2.2. The reference problem

Here, in order to simplify the presentation, the problem will be formulated for an assembly of only two substructures, but this can easily be generalized to an assembly of n substructures. For the purpose of illustration, let us consider two thin, isotropic and elastic Kirchhoff–Love plates. First, we assume that they are homogeneous.

Let us study the steady-state vibrations of this assembly. Classically, all quantities are defined in the complex domain: an amplitude $Q(X)$ is associated with $Q(X) \cdot \exp(i\omega t)$.

Let $S_1$ and $S_2$ be two plates, $\partial S_1$ and $\partial S_2$ their respective boundaries solicited harmonically at a fixed angular frequency $\omega$ and $\Gamma_{12}$ the perfect interface between $S_1$ and $S_2$. The boundary conditions are as follows, with $l = (1, 2)$: deflection $w_l^d$ on part $\partial_{wd}S_l$ of $\partial S_l$; slope $w_{nl}^d$ on $\partial_{wd}S_l$; bending moment $M_l^d$ on $\partial_{Ml}S_l$; finally, Kirchhoff shear $K_l^d$ on $\partial_{Kl}S_l$ (Fig. 1). The quantities of interest are the deflections $w_l$ and the moments $M_l$.

For $S_l$, let us introduce the space $\mathcal{S}_{adl}$ of the displacement-stress pairs $s_l = (w_l, M_l)$ defined on $S_l$ such that:

$$ s_l \in \mathcal{S}_{adl} \iff \left\{ \begin{array}{l} s_l \in \mathcal{U} \times \mathcal{S} \quad (\text{set of finite-energy fields}) \\ \Delta \cdot \Delta w_l - k_l^4 w_l = 0 \quad \text{on } S_l \\ \mathbf{M}_l = (1 + i\eta_l) \mathbf{K}_{PSl} \mathcal{X}(w_l) \\ \text{for } l = 1, 2 \\ \\ \text{with } k_l^4 = \frac{3\rho_l \omega^2 (1-nu_l^2)}{(1+\eta_l)E_l h_l^2} \end{array} \right. \quad (1) $$

where $k_l$ is the wave number and $\rho_l$, $v_l$, $\eta_l$, $E_l$ and $h_l$ designate respectively the density, Poisson's ratio, structural damping coefficient, Young's modulus and thickness of plate $S_l$. $\mathbf{K}_{PSl}$ is Hooke's tensor for plane stress and $\mathcal{X}$ is

<figure><img src="image_3.png" /><figcaption>Fig. 1. The reference problem.</figcaption></figure>
the curvature operator. Let us express $\mathbf{M}_l$ and $\mathcal{X}$ in cartesian coordinates:

$$[\mathbf{M}_l] = \begin{bmatrix} M_{lx} \\ M_{ly} \\ M_{lxy} \end{bmatrix} \qquad (2)$$

$$[\mathcal{X}] = - \begin{bmatrix} \frac{\partial^2 w}{\partial x^2} \\ \frac{\partial^2 w}{\partial y^2} \\ 2 \frac{\partial^2 w}{\partial x \partial y} \end{bmatrix} \qquad (3)$$

then, $\mathbf{K}_{PSl}$ is defined by:

$$[\mathbf{K}_{PSl}] = \frac{E_l h_l^3}{12(1 - v_l^2)} \begin{bmatrix} 1 & v_l & 0 \\ v_l & 1 & 0 \\ 0 & 0 & \frac{1-v_l}{2} \end{bmatrix} \qquad (4)$$

The fields of $\mathcal{S}_{adl}$ are admissible in the sense that they verify both the local equilibrium and the constitutive relation exactly. The reference problem is:

**Problem 1.** Find $s_l$ such that:

$$\left\{
\begin{array}{l}
s_l \in \mathcal{S}_{adl} \\
w_l = w_l^d \quad \text{on } \partial_{wd} S_l \\
w_{l,n_l} = w_{n_l}^d \quad \text{on } \partial_{wdn} S_l \\
M_{lnl} = M_l^d \quad \text{on } \partial_{Md} S_l \\
K_{lnl} = K_l^d \quad \text{on } \partial_{Kd} S_l \\
\text{with } l = 1, 2 \\
w_1 = w_2 \quad \text{on } \Gamma_{12} \\
w_{1,n_1} = w_{2,n_2} \quad \text{on } \Gamma_{12} \\
M_{1n_1} + M_{2n_2} = 0 \quad \text{on } \Gamma_{12} \\
K_{1n_1} + K_{2n_2} = 0 \quad \text{on } \Gamma_{12}
\end{array}
\right.
\qquad (5)$$

with $M_{lnl} = n_l \mathbf{M}_l n_l$

and $K_{lnl} = \underline{n_l} \text{div}[\mathbf{M}_l] + (\underline{t_l} \mathbf{M}_l \underline{n_l})_{,l}$.

The reference problem has, at the most, one solution if:

* Hooke's tensors are positive definite.

* $\eta_l > 0$.

## 2.3. The variational formulation associated with the VTCR

The VTCR is primarily a global formulation of the boundary conditions and transmission conditions in terms of both displacements and forces across $\Gamma_{12}$. With this understanding, the problem becomes:

**Problem 2.** Find $s_l$ such that:

$$\left\{
\begin{array}{l}
s_l \in \mathcal{S}_{adl} \\
\sum_l (A_{\hat{c}S_l}(\delta s_l, s) - L_{\hat{c}S}(\delta s_l)) + \sum_{\substack{l,m \\ l<m}} C_{\Gamma_{lm}}(\delta s_l, s_l, \delta s_m, s_m) = 0 \\
\forall \delta s_l \in \mathcal{S}_{adl} \\
\text{for } l=1,2 \quad m=1,2
\end{array}
\right.
\qquad (6)$$

with the boundary conditions:

$$\begin{align}
A_{\hat{c}S_l}(\delta s_l, s_l) - L_{\hat{c}S_l}(\delta s_l) & \nonumber \\
&= \operatorname{Re} \left\{ -i\omega \left( -\int_{\partial_{wd}S_l} \delta K_{lnl} (w_l - w_l^d)^* dL \right. \right. \nonumber \\
& \qquad \left. + \int_{\partial_{wdn}S_l} \delta n_l \mathbf{M}_l n_l (w_{l,n_l} - w_{l,n_l}^d)^* dL \right. \nonumber \\
& \qquad \left. + \int_{\partial_{Md}S_l} (n_l \mathbf{M}_l n_l - M_l^d) \delta w_{l,n_l}^* dL \right. \nonumber \\
& \qquad \left. - \int_{\partial_{Kd}S_l} (K_{lnl} - K_l^d) \delta w_l^* dL \right) \right\} \tag{7}
\end{align}$$

and with the transmission conditions:

$$C_{\Sigma_{lm}}(\delta s_l, s_l, \delta s_m, s_m) = \operatorname{Re} \left\{ -i\omega \left( \int_{\Gamma_{lm}} \frac{1}{2} \left[ \begin{aligned}[t]
& -\delta(K_{lnl} - K_{mnm})(w_l - w_m)^* \\
& + \delta(n_l \mathbf{M}_l n_l - \underline{n}_m \mathbf{M}_m \underline{n}_m)(w_{l,n_l} - w_{m,n_m})^* \\
& + (\underline{n}_l \mathbf{M}_l n_l + \underline{n}_m \mathbf{M}_m \underline{n}_m)\delta(w_{l,n_l} + w_{m,n_m})^*
\end{aligned} \right] - (K_{lnl} + K_{mnm})\delta(w_l + w_m)^* \right] dL \right\} \qquad (8)$$

Re[$Q$] designates the real part of $Q$ and $Q^*$ designates the conjugate of $Q$. Note that the terms which appear in this formulation are powers.

In practice, the VTCR procedure consists of two stages: first, admissible fields are built; then, the boundary conditions are taken into account through the discretization of the variational formulation associated with a set of admissible fields chosen among the elements of $\mathcal{S}_{adl}$.

### 2.3.1. Construction of admissible fields

In practice, the first stage of the VTCR consists in building admissible fields. Let us define a subset of $\mathcal{S}_{adl}$ which can contain local interior vibration modes described using complex rays of the nth order. If the damping coefficient $\eta_l$ ($\eta_l \ll 1$) is small, the basic mode associated with the wave vector $\underline{P}_l$ can be written as:

$$\left\{
\begin{array}{l}
w_l(\underline{X}, \underline{P}_l) = e^{i\sqrt{\omega} \underline{P}_l \cdot \underline{X}} e^{\frac{\eta_l}{4}\sqrt{\omega} \underline{P}_l \cdot \underline{X}} [\underline{V}(\underline{X}, \underline{P}_l)][\underline{a}_l(\underline{P}_l)] \\
= [\underline{W}(\underline{X}, \underline{P}_l)][\underline{a}_l(\underline{P}_l)], \\
\mathbf{M}_l(\underline{X}, \underline{P}_l) = (1 + i\eta_l)\mathbf{K}_{PSl}\mathcal{X}(w(\underline{X}, \underline{P}_l)), \\
\underline{P}_l \in \mathcal{C}_{adl}.
\end{array}
\right.
$$

The admissibility relation requires that the locus of the end of $\underline{P}_l$ be a curve $\mathcal{C}_{adl}$ characteristic of the material. $[\underline{V}(\underline{X}, \underline{P}_l)][\underline{a}_l(\underline{P}_l)]$ is a polynomial expression of degree $n$ in $\underline{X}$. The small-length part appears explicitly. The unknown generalized amplitudes $\underline{a}_l(\underline{P}_l)$ are independent of $\underline{X}$ and, therefore, are large-length quantities. The edge and corner modes are built in the same way.

For the interior zone, we take into account all the admissible directions of the rays. The solution $w_l^{\text{sol}}$ is written as:
$$w_l^{\text{sol}}(\underline{X}) = \int_{\mathcal{C}_{\text{adl}}} w_l(\underline{X}, P_l) \mathrm{d}\mathcal{C}. \qquad (9)$$

The integrals are calculated along the curve $\mathcal{C}_{\text{adl}}$. Clearly, $w_l^{\text{sol}}$ depends on the generalized amplitude $a_l$ relative to $P_l$ defined on $\mathcal{C}_{\text{adl}}$.

### 2.3.2. The discretized form of the VTCR

In order to derive approximations from the VTCR, one discretizes the curve $\mathcal{C}_{\text{adl}}$ using finite elements; for example, $P_l \to a_l(P_l)$ is assumed to be constant within each element. In other terms, $a_l^h$ is associated with $\mathcal{C}_{\text{adl}}^{r_l}$ (i.e., $\mathcal{C}_{\text{adl}}$ discretized with $r_l$ elements). The subspace $\mathcal{S}_{\text{adl}}^h$ of the approximations deduced from $\mathcal{S}_{\text{adl}}$ yields a discretized variational formulation leading to a system of equations of dimension $\sum_l r_l$ in the complex domain:

$$[\mathbf{A}^h][a^h] = [\underline{L}_d^h], \qquad (10)$$

$[a^h]$ designates the column of all the amplitudes $a_l^h$.

### 2.4. Effective quantities

For relatively high frequencies, only effective quantities associated with the generalized amplitudes are meaningful. Indeed, the precise spatial distribution of $s_l^h$ is meaningless from a mechanical point of view. On a given substructure $S$ with a characteristic dimension $L$, the quantities involved are expressed as $q(\underline{X})$, which can represent, for example, the elastic energy, the kinetic energy or the dissipation work; in the vicinity of any test point $\underline{X}_0$, the associated effective quantities are defined on a square patch $S_{\underline{X}_0}$ with side $L$ by:

$$q_{\text{eff}}(\underline{X}_0) = \frac{1}{L^2} \int_{S_{\underline{X}_0}} q(\underline{X}) \mathrm{d}S. \qquad (11)$$

These effective quantities are large-wavelength quantities. Compared to classical quantities such as the displacements of points, which can vary considerably due to slight variations of input data, they are much less sensitive. Effective quantities can also be defined on the edge of $S$.

### 2.5. Example

A three-dimensional assembly of plates is used to model a part of the hull of a ship (Fig. 2). The 13 plates are made of sheet steel. Their mechanical properties are: $E = 210$ GPa; $\eta = 2\%$; $v = 0.3$; $\rho = 7800$ kg/m$^3$; $h = 2$ mm. The harmonic excitation is a distributed momentum $M_d$ on the edge $P_{17}P_{18}$: $M_d = 1$ Nm/m; $f = 500$ Hz. All boundaries are simply supported. We used the program CORAY-MF (COMplex RAYs for Medium Frequencies) developed by Rouch at LMT-Cachan to compute the VTCR solution. This solution and the solution obtained with finite elements are very similar (Fig. 3).

## 3. Extension of the VTCR to heterogeneous structures

### 3.1. The reference problem

Here, we focus on the construction of the superelement associated with an isolated heterogenous substructure (Fig. 4). Let us consider the steady-state vibrations of a thin, isotropic and elastic Kirchhoff-Love plate with local heterogeneity. Each type of heterogeneity imposes specific conditions on its boundary. In order to test our strategy, let us consider the case of a plate with an arbitrarily-shaped hole. Its boundary $\Gamma_0$ is assumed to be free. Therefore, the boundary conditions on $\Gamma_0$ are $M^d = 0$ and $K^d = 0$. The notations are the same as before, except for dropping the suffix $l$ in order to simplify the equations. The space $\mathcal{S}_{\text{ad}}$ is defined exactly like $\mathcal{S}_{\text{adl}}$ in the first part. The reference problem becomes:

**Problem 3.** Find $s$ such that:

$$ \left\{ \begin{array}{l} s \in \mathcal{S}_{\text{ad}} \\ w = w^d \quad \text{on } \partial_{\text{wd}} S \\ w_{,n} = w_n^d \quad \text{on } \partial_{\text{wd}} S \\ M_n = M^d \quad \text{on } \partial_{\text{Md}} S \\ K_n = K^d \quad \text{on } \partial_{\text{Kd}} S \\ M_n = M^d = 0 \quad \text{on } \Gamma_0 \\ K_n = K^d = 0 \quad \text{on } \Gamma_0 \\ \text{with } M_n = \underline{n} \mathbf{M} \underline{n} \\ \text{and } K_n = \underline{n} \operatorname{div}[\mathbf{M}] + (\underline{t} \mathbf{M} \underline{n})_{,t} \end{array} \right. \qquad (12) $$

### 3.2. Direct application of the VTCR

#### 3.2.1. The variational formulation associated with the VTCR

The first option is to treat $\Gamma_0$ like $\partial S$. This leads to Problem 4:

**Problem 4.** Find $s$ such that:

$$ \left\{ \begin{array}{l} s \in \mathcal{S}_{\text{ad}} \\ (A_{\partial S}(\delta s, s) + A_{\Gamma_0}(\delta s, s)) - (L_{\partial S}(\delta s) + L_{\Gamma_0}(\delta s)) = 0 \quad (13) \\ \forall \delta s \in \mathcal{S}_{\text{ad}} \end{array} \right. $$

with:

$$ A_{\Gamma_0}(\delta s, s) - L_{\Gamma_0}(\delta s) = \operatorname{Re} \left\{ -i\omega \left( \int_{\Gamma_0} \underline{n} \mathbf{M} \underline{n} \delta w_{,n}^* dL - \int_{\Gamma_0} K_n \delta w^* dL \right) \right\}. $$

The term $(A_{\partial S}(\delta s, s) - L_{\partial S}(\delta s))$ was developed while defining Problem 2 in Eq. (7).

#### 3.2.2. Limitations of the direct formalism

While the influence of heterogeneity is taken into account through the terms integrated over $\Gamma_0$, we would like the treatment of local heterogeneity to be distinguished from that of the master structure. Indeed, from a mechanical
<figure><img src="image_4.png" /><figcaption>Fig. 2. Definition of the geometry.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. VTCR solution (left) and FE solution (right): $|w|$ (m).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. The reference problem.</figcaption></figure>

point of view, the hole creates only a perturbation of the vibratory response. One should note that this interpretation is not valid for all types of heterogeneity: for example, a long stiffener would behave differently. The case of stiffeners will be developed in a subsequent paper.

### 3.3. Modification of the VTCR in order to deal with local heterogeneity

#### 3.3.1. The new modified formulation

The modification consists in taking the perturbation due to local heterogeneity into account as soon as the admissible fields of $\mathcal{S}_{ad}$ have been built. The idea is to develop a priori, once and for all, families of admissible fields which verify the boundary conditions resulting from each type of heterogeneity. Therefore, let us define the set $\mathcal{S}_{ad}^c$ (i.e., corrected $\mathcal{S}_{ad}$) of the pairs $s$ such that:

$$ s \in \mathcal{S}_{ad}^c \iff \begin{cases} s \in \mathcal{S}_{ad} \\ M_n = M^d = 0 \quad \text{on } \Gamma_0 \\ K_n = K^d = 0 \quad \text{on } \Gamma_0 \end{cases} \qquad (14) $$

This leads to Problem 5:
**Problem 5.** Find $s$ such that:

$$ \left\{ \begin{array}{l} s \in \mathcal{S}_{ad}^{c} \\ A_{\delta S}(\delta s, s) - L_{\delta S}(\delta s) = 0 \\ \forall \delta s \in \mathcal{S}_{ad}^{c} \end{array} \right. \qquad (15) $$

### 3.3.2. Principle of construction of the corrected admissible fields

We want to build fields which correspond to the new definition of admissibility. Thus, we seek $s^c \in \mathcal{S}_{ad}^c$ among the admissible pairs $\bar{s}$ in the sense of $\mathcal{S}_{ad}$ which already exist. $\bar{s}$ is a nth complex ray which results from the linear combination of elementary complex rays $\bar{s}_e = (\bar{w}_e, \mathbf{M}_e)$, associated with the amplitudes $a_e$:

$$ \bar{w}(\underline{X}, \underline{P}) = \sum_{e=0}^{n} \bar{w}_e(\underline{X}, \underline{P}) a_e(\underline{P}). \qquad (16) $$

Each $\bar{s}_e$ has a unitary amplitude. The corresponding correction field $\tilde{s}_e$ is calculated so that the sum $s_e^c$ of $\bar{s}_e$ and $\tilde{s}_e$ verifies the boundary conditions resulting from heterogeneity. Of course, this correction must belong to $\mathcal{S}_{ad}$, too. One can view each $s^c$ as the superimposition of an incident field $\bar{s}$ and a diffracted field $\tilde{s}$ due to heterogeneity:

$$ (w_e^c, \mathbf{M}_e^c) = (\bar{w}_e, \mathbf{M}_e) + (\tilde{w}_e, \tilde{\mathbf{M}}_e). \qquad (17) $$

### 3.3.3. Practical construction of the corrected admissible fields

For our reference problem, the corrections $\tilde{s}_e(\underline{X}, \underline{P})$ associated with an incident field of unit amplitude $\bar{s}_e(\underline{X}, \underline{P})$ are defined by Problems 6 and 7:

**Problem 6.** Find $\tilde{s}_e$ such that:

$$ \left\{ \begin{array}{l} \tilde{s}_e \in \mathcal{S}_{ad}, \\ \tilde{M}_{en} = -\tilde{M}_{en}|_{\Gamma_0}, \\ \tilde{K}_{en} = -\tilde{K}_{en}|_{\Gamma_0}. \end{array} \right. \qquad (18) $$

In order to take into account the radiating aspect of the wave diffracted due to heterogeneity, we wish to calculate $\tilde{s}_e$ in polar coordinates with an origin located at the “center of gravity” of the hole. In order to do that, we define an “equivalent circular hole” with appropriate boundary conditions: this is an imaginary circular contour $\Gamma_{0e}$ with radius $R$, centered at the origin of the polar coordinates system and containing the heterogeneous property (Fig. 5). Then, $\tilde{w}_e$ is sought in the classical form:

$$ \tilde{w}_e(r, \theta) = \sum_{m \in \mathbb{Z}} C_{fm}(\bar{w}) H_{fm}(r, \theta) \qquad (19) $$

$$ \text{with } H_{fm} = H_m^{(f)}(\alpha_f kr) e^{im\theta} \quad \text{and} \quad \alpha_1 = i, \alpha_2 = 1, \qquad (20) $$

where $H_m^{(1)}$ and $H_m^{(2)}$ denote the Hankel functions of the first and second kind respectively. The solution is a combination of an outbound propagating wave and a boundary wave which decreases rapidly. Both waves are emitted by

<figure><img src="image_7.png" /><figcaption>Fig. 5. Definition of the imaginary circular contour.</figcaption></figure>

$\Gamma_0$ and neither reflects at infinity. The coefficients $C_{fm}$ are determined through the boundary conditions on $\Gamma_0$, which are relocated on $\Gamma_{0c}$. These relocated boundary conditions are the restrictions to the contour $\Gamma_{0c}$ of the solution to Problem 6. In fact, outside $\Gamma_{0c}$, the solution to Problem 6 is the same as the solution to Problem 7.

**Problem 7.** Find $\tilde{s}_e$ such that:

$$ \left\{ \begin{array}{l} \tilde{s}_e \in \mathcal{S}_{ad} \\ \tilde{w}_e = \tilde{w}_e^d|_{\Gamma_{0c}} \\ \tilde{w}_{e,n} = \tilde{w}_{en}^d|_{\Gamma_{0c}} \end{array} \right. \qquad (21) $$

Problem 6 can be solved by a boundary element method (e.g., [18]) or an extension of the VTCR to unbounded media which is currently under development. Then, $\tilde{w}_e^d|_{\Gamma_{0c}}$ and $\tilde{w}_{en}^d|_{\Gamma_{0c}}$ are expanded into Fourier series:

$$ \left\{ \begin{array}{l} \tilde{w}_e^d(R, \theta) = \sum_{m \in \mathbb{Z}} \tilde{w}_e^{dm}(R) e^{im\theta} \\ \tilde{w}_{e,n}^d(R, \theta) = \sum_{m \in \mathbb{Z}} \tilde{w}_{e,n}^{dm}(R) e^{im\theta} \end{array} \right. \qquad (22) $$

$$ \text{with } \tilde{w}_e^{dm}(R) = \frac{1}{2\pi} \int_0^{2\pi} e^{-im\theta} \tilde{w}_e^d(R, \theta) d\theta, \qquad (23) $$
$$ \text{and } \tilde{w}_{e,n}^{dm}(R) = \frac{1}{2\pi} \int_0^{2\pi} e^{-im\theta} \tilde{w}_{e,n}^d(R, \theta) d\theta. $$

The coefficients $C_{fm}$ are calculated by inverting linear systems of dimension 2:

$$ \left\{ \begin{array}{l} C_{1m} H_m^{(1)}(ikr) + C_{2m} H_m^{(2)}(kr) = \tilde{w}_e^{dm}(R), \\ C_{1m} H_m^{(1)'}(ikr) + C_{2m} H_m^{(2)'}(kr) = \tilde{w}_{e,n}^{dm}(R). \end{array} \right. \qquad (24) $$

This can be expressed in matrix form by introducing the operator $C_{\Gamma_0/\Gamma_{0c}}$: it can be considered as a projection operator of the correction on the basis of the functions $H_{\Gamma_{0c}}$ which are associated with the contour $\Gamma_{0c}$:

$$ \tilde{w}_e(r, \theta) = [\mathbf{H}_{\Gamma_0/\Gamma_{0c}}(r, \theta)][\mathbf{C}_{\Gamma_0/\Gamma_{0c}}(\bar{w}_e)]. \qquad (25) $$

The expression of the field $w_e^c$, adapted from the elementary complex ray $\bar{w}_e$, is:

$$ w_e^c(\underline{X}, \underline{P}) = \bar{w}_e(\underline{X}, \underline{P}) + [\mathbf{H}_{\Gamma_{0c}}(r, \theta)][\mathbf{C}_{\Gamma_0/\Gamma_{0c}}(\bar{w}_e)], \qquad (26) $$

and the complete corrected complex ray can be expressed as:

$$ w^c(\underline{X}, \underline{P}) = ([\bar{\mathbf{W}}(\underline{X}, \underline{P})] + [\mathbf{H}_{\Gamma_{0c}}(r, \theta)][\mathbf{C}_{\Gamma_0/\Gamma_{0c}}(\bar{\mathbf{W}})])[\underline{a}(\underline{P})]. \qquad (27) $$

In short, the superelement corresponding to a heterogeneous substructure includes a family of additional basic modes, weighted in order to verify conditions on the
boundary of the discontinuity a priori. The weights corresponding to ranges of non-dimensional parameters associated with each type of local heterogeneity are processed and stored in advance in sorts of “libraries”. This approach can be generalized to two or more zones containing local heterogeneities.

### 3.3.4. Discretized expression of the modified VTCR

The expression of the solution $w^{\text{sol}}$ for the heterogeneous structure is:

$$w^{\text{sol}}(\underline{X}) = \int_{\mathcal{C}_{ad}} w^c(\underline{X}, \underline{P}) d\mathcal{C}. \qquad (28)$$

The integral is evaluated over the admissibility domain $\mathcal{C}_{ad}$ of the field $\bar{w}$. Before discretizing the variational formulation, all the contributions of each correction function to the correction of all the elementary complex rays defined over the admissibility domain are put together. Then, the generalized amplitudes $\hat{a}$ associated with the correction functions are defined. These generalized amplitudes, considered as the new unknowns, are linked with the $\underline{a}(\underline{P})$ by the relation:

$$[\hat{a}] = \int_{\mathcal{C}_{ad}} [\mathbf{C}_{\Gamma_0/\Gamma_{0C}}(\bar{\mathbf{W}}(\underline{X}, \underline{P}))][\underline{a}(\underline{P})] d\mathcal{C}. \qquad (29)$$

The discretization is performed as follows: an approximate expression of $\bar{w}$ is obtained by discretizing $\mathcal{C}_{ad}$ according to the process described in Section 2.3.2; an approximate expression of $\tilde{w}$ is obtained by truncating the series of correction functions to $2 \times (2M + 1)$ terms, $M$ being the maximum order of the harmonics kept; finally, the subspace $\mathcal{S}_{ad}^{ch}$ of the approximations deduced from $\mathcal{S}_{ad}^{c}$ is defined by:

$$s^{ch} \in \mathcal{S}_{ad}^{ch} \iff \left\{ \begin{array}{l} w^{ch}(\underline{X}) = [\bar{\mathbf{W}}^h(\underline{X})][\underline{a}^h] + [\mathbf{H}_{\Gamma_{0C}}^h(r, \theta)][\underline{a}^h] \\ \mathbf{M}^{ch}(\underline{X}) = [\bar{\underline{\Lambda}}^h(\underline{X})][\underline{a}^h] + [\bar{\underline{\Xi}}_{\Gamma_{0C}}^h(r, \theta)][\underline{a}^h] \\ [\underline{a}^h] = [\mathbf{D}_{\Gamma_0/\Gamma_{0C}}^h][\underline{a}^h] \end{array} \right. \qquad (30)$$

$\mathbf{D}_{\Gamma_0/\Gamma_{0C}}^h$, linking the generalized amplitudes of uncorrected fields to the ones of correction functions, designates the correction operator. The modified variational formulation defined by Eq. (15) and the subspace $\mathcal{S}_{ad}^{ch}$ defined by Eq. (30) lead to a constrained linear system of equations of dimension $r + r^c$ in the complex domain, where $r^c$ denotes the number of correction fields:

$$\begin{bmatrix} \bar{\underline{\mathbf{A}}}^{ch} & \bar{\underline{\mathbf{A}}}^{ch} \\ \bar{\underline{\tilde{\mathbf{A}}}}^{ch} & \bar{\underline{\tilde{\mathbf{A}}}}^{ch} \end{bmatrix} \begin{bmatrix} \underline{a}^h \\ \underline{a}^h \end{bmatrix} = \begin{bmatrix} \underline{L}_d^{ch} \\ \underline{L}_d^{ch} \end{bmatrix} \qquad (31)$$
with $[\underline{a}^h] = [\mathbf{D}_{\Gamma_0/\Gamma_{0C}}^h][\underline{a}^h]$

This formalism enables us to associate each heterogeneous plate superelement with a stiffness matrix and a right-hand side relative to uncorrected fields and to correction functions which enrich the space of the basis functions. The correction functions are adjusted to the discontinuity through the correction operator $\mathbf{D}_{\Gamma_0/\Gamma_{0C}}$ which links the generalized

amplitudes of the uncorrected fields to those of the correction functions.

The constraint on the generalized amplitudes can be taken into account during the resolution, for example by introducing $r^c$ Lagrange multipliers which form the vector $\mu$. The new linear system to be solved has the dimension $r + 2 \times r^c$ in the complex domain, and its expression is:

$$\begin{bmatrix} \bar{\underline{\mathbf{A}}}^{ch} & \bar{\underline{\mathbf{A}}}^{ch} & -\mathbf{D}_{\Gamma_0/\Gamma_{0C}}^{hT} \\ \bar{\underline{\tilde{\mathbf{A}}}}^{ch} & \bar{\underline{\tilde{\mathbf{A}}}}^{ch} & \mathbf{I} \\ -\mathbf{D}_{\Gamma_0/\Gamma_{0C}}^{h} & \mathbf{I} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \underline{a}^h \\ \underline{a}^h \\ \mu \end{bmatrix} = \begin{bmatrix} \underline{L}_d^{ch} \\ \underline{L}_d^{ch} \\ \underline{0} \end{bmatrix}. \qquad (32)$$

### 3.4. Example

This method can be applied, for example, to the problem of a rectangular steel plate with length $L = 1$ m, width $l = 0.8$ m and thickness $h = 0.8$ mm containing a circular hole 0.11 m in diameter at the center. The material's properties are $E = 210$ GPa, $\eta = 0.001$, $v = 0.3$ and $\rho = 7800$ kg/m$^3$. The harmonic excitation is a distributed force $F_d$ on one side: $f = 330$ Hz; $F_d = 1$ N/m; the plate is clamped on the opposite side. The modified VTCR solution was calculated using 84 interior modes, 48 edge modes and eight harmonics of the correction functions (Fig. 6). This solution and the solution obtained using finite elements (Fig. 7) are very similar.

<figure><img src="image_8.png" /><figcaption>Fig. 6. VTCR solution: $|w|$ (μm).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Finite element solution: $|w|$ (μm).</figcaption></figure>
# 4. Sensitivity analysis

## 4.1. Definition of the problem

The objective is to evaluate the influence of discontinuities with uncertain positions, sizes or directions on the structural response. The point of departure is the library of corrected fields defined in Section 3.3.3. This library includes the elements, stored in non-dimensional form for each family of discontinuities (e.g., circular holes), which are necessary to rebuild the specific corrections dedicated to particular members of the family (for example, the hole whose diameter is twice the wavelength of the flexural waves in the plate). These elements, defined in the local coordinate system associated with the equivalent discontinuity, constitute the basis of the correction functions and the projection operator of the correction.

In the case of a perfectly known structure (i.e., in the absence of uncertainties), using the library to calculate the superelement requires the positioning of the equivalent discontinuity's imaginary circular contour in the global coordinate system of the structure: then, the polar coordinate system associated with the correction functions is centered and oriented with respect to the boundaries of the structure.

When the parameters of the discontinuity vary about their nominal values, the imaginary circular contour also varies in the global coordinate system: this may affect its size (e.g., when the diameter of a hole changes), its position (when the center of the hole moves) or its orientation (when the axis of symmetry of an elliptical hole tilts). The expressions of the correction functions in the global coordinate system depend on the evolution of the imaginary contour; consequently, the matrix of the interactions among the basis functions and the right-hand side vector of the interactions between the fields and the boundary conditions also vary. Then, each set of parameters requires a full calculation. Because many very similar problems are considered, the sensitivity analysis is very expensive. It is of paramount importance to choose a well suited and efficient calculation method.

## 4.2. Principle of the method: use of a single imaginary circular contour

In the framework of structural mechanics, a multiresolution strategy has been developed: it is dedicated to the LATIN method which deals with non-linear evolution problems [4]. The fundamental idea is to exploit the same basis of functions to describe the solution of similar problems.

An analogous principle is introduced here. It consists in using a single imaginary circular contour for all the discontinuities corresponding to the range of uncertain parameters we wish to consider (Fig. 8). The basis functions which are associated with this contour do not change when the parameters vary. The interaction matrix and the right-hand side remain unchanged, which drastically reduces the

<figure><img src="image_10.png" /><figcaption>Fig. 8. Definition of the single imaginary contour.</figcaption></figure>

time needed for calculating the matrices defining the super-element of the heterogeneous substructure.

The only change concerns the constraint relations between the generalized amplitudes of the corrected fields and those of the correction functions: these are deduced from the constraint relations which were established when the equivalent discontinuity was defined in order to build the library of corrected fields.

## 4.3. Practical aspects involved in building a model of uncertain local heterogeneity

The single imaginary contour $\Gamma_U$ is chosen such that it includes all the imaginary contours associated with the equivalent discontinuities envisaged. The set of parameters relative to each discontinuity is evaluated using, for example, the Monte Carlo method. The origin of the common coordinate system can be chosen such that it coincides with the “center of gravity” of the discontinuity with nominal parameters.

The field $s^c$ (an adjusted version of field $\bar{s}$ for the particular discontinuity marked i), whose contour is $\Gamma_0^i$ and whose equivalent discontinuity has the contour $\Gamma_{0c}^i$, is projected onto the basis functions associated with the common imaginary contour through an operator $C_{\Gamma_0^i/\Gamma_U}$ defined as follows:

$$ w^c(\underline{X}, \underline{P}) = \left( [\bar{W}(\underline{X}, \underline{P})] + [\mathbf{H}_{\Gamma_U}(r, \theta)] \left[ \mathbf{C}_{\Gamma_0^i/\Gamma_U}(\mathbf{H}_{\Gamma_0^i}) \right] \times \left[ \mathbf{C}_{\Gamma_0^i/\Gamma_{0c}}(\bar{\mathbf{W}}) \right] \right) [\underline{a}(\underline{P})]. \quad (33) $$

$C_{\Gamma_0^i/\Gamma_U}$ can be viewed as an operator which projects the basis functions relative to contour $\Gamma_{0c}$ onto the basis functions relative to the common contour $\Gamma_U$. The process of constructing $C_{\Gamma_0^i/\Gamma_U}$ is analogous in every respect to the process described in Section 3.3.3 to express the projection of $\tilde{w}$ onto the basis of the functions $\mathbf{H}_{\Gamma_{0c}}$.

When one discretizes the variational formulation, the formalism introduced in Section 3.3.4 leads to the use of generalized amplitudes associated with common correction functions $\mathbf{H}_{\Gamma_U}$ defined by:

$$ \begin{aligned} [\hat{a}] &= \int_{\mathcal{E}_{n,d}} \left[ \mathbf{C}_{\Gamma_0^i/\Gamma_U}(\mathbf{H}_{\Gamma_0^i}) \right] \left[ \mathbf{C}_{\Gamma_0^i/\Gamma_{0c}}(\bar{\mathbf{W}}(\underline{X}, \underline{P})) \right] [\underline{a}(\underline{P})] \mathrm{d}s \\ &= \left[ \mathbf{D}_{\Gamma_0^i/\Gamma_U} \right] [\underline{a}(\underline{P})] \end{aligned} \quad (34) $$
The specificity of the corrections relative to discontinuity i are taken into account through the operator $D_{Γ_0^i/Γ_U}$. System (32) can be partially inverted so that $\underline{\mu}$ and $D_{Γ_0^i/Γ_U}$ remain parameters (Schur method):

$$ \left[ \frac{\underline{a}^h}{\tilde{\underline{a}}^h} \right] = \left[ \begin{array}{cc} \bar{\underline{\mathbf{A}}}^{ch} & \tilde{\underline{\mathbf{A}}}^{ch} \\ \tilde{\underline{\mathbf{A}}}^{ch} & \tilde{\underline{\mathbf{A}}}^{ch} \end{array} \right]^{-1} \left( \left[ \frac{\bar{\underline{L}}_d^{ch}}{\tilde{\underline{L}}_d^{ch}} \right] - \left[ \begin{array}{c} -\mathbf{D}_{\Gamma_0^i/\Gamma_{0C}}^{hT} \\ \mathbf{I} \end{array} \right] \left[ \underline{\mu} \right] \right) \quad (35) $$

with

$$ [\underline{\mu}] = [\mathbf{M}]^{-1} [\mathbf{V}] \quad (36) $$

where

$$ [\mathbf{M}] = \left[ -\mathbf{D}_{\Gamma_0^i/\Gamma_{0C}}^h \quad \mathbf{I} \right] \left[ \begin{array}{cc} \bar{\underline{\mathbf{A}}}^{ch} & \tilde{\underline{\mathbf{A}}}^{ch} \\ \tilde{\underline{\mathbf{A}}}^{ch} & \tilde{\underline{\mathbf{A}}}^{ch} \end{array} \right]^{-1} \left[ \begin{array}{c} -\mathbf{D}_{\Gamma_0^i/\Gamma_{0C}}^{hT} \\ \mathbf{I} \end{array} \right] \quad (37) $$

and

$$ [\mathbf{V}] = \left[ -\mathbf{D}_{\Gamma_0^i/\Gamma_{0C}}^h \quad \mathbf{I} \right] \left[ \begin{array}{cc} \bar{\underline{\mathbf{A}}}^{ch} & \tilde{\underline{\mathbf{A}}}^{ch} \\ \tilde{\underline{\mathbf{A}}}^{ch} & \tilde{\underline{\mathbf{A}}}^{ch} \end{array} \right]^{-1} \left[ \frac{\bar{\underline{L}}_d^{ch}}{\tilde{\underline{L}}_d^{ch}} \right] \quad (38) $$

From here on, for each discontinuity, the only calculation consists in building the operator $D_{Γ_0^i/Γ_U}$ using the data previously stored in the library: thus, the generalized amplitudes are obtained almost instantaneously through the inversion of matrix $\mathbf{M}$ which is cheap compared to the inversion of the whole system (32).

The formalism associated with the parametric model of local heterogeneity lays the foundations of a multiresolution strategy which is fitted for the VTCTR. Coupled with a response surface technique [10,31], it makes it possible to deal with many kinds of problems linked to the study of the sensitivity of the response to a variation of heterogeneity parameters: parametric analysis, Monte Carlo simulation, reliability analysis.

### 4.4. Example

Let us consider a plate with a hole, using the same dimensions and material characteristics as in the example of Section 3.4. We wish to know the influence of a variation in the position of the hole (due, e.g., to manufacturing inaccuracies) on the vibrational response (Fig. 9).

For the sake of simplicity, let us consider only the possibility of an offset $y$ parallel to the width of the plate. The frequency of the harmonic excitation is $f = 277$ Hz; the diameter of the hole is $D = 0.11$ m (the ratio between $D$ and the wavelength $\lambda$ is 0.58). Let us assume that $y$ varies about its nominal value according to a normal centered Gaussian distribution with standard deviation $\sigma = 0.03$ m. A sampling of 50000 possible values of $y$ is obtained through a Monte Carlo draw ($\tilde{y}$ designates the reduced variable associated with $y$). Distribution tails are truncated at $\pm 3\sigma$.

First, the response curve giving the global effective displacement $|w|$ as a function of $y$ is calculated using the

<figure><img src="image_11.png" /><figcaption>Fig. 9. Variation in the position of the hole.</figcaption></figure>

multiple resolution process we just defined (Fig. 10). It is an interpolation function built on a small number of values of the diameter of the hole. This number, of the order of a few tens, is considered to be enough when the precision of the interpolated global effective displacement does not improve significantly by further refining.

Fig. 11 shows the probability density function $f_{|w|}$ associated with $|w|$. It reveals three peaks of values of $|w|$ which are more probable than others around 2.5 µm, 3.75 µm and 8.75 µm. These values correspond to local extrema of $|w|$, hence a possible accumulation of the samples of $f_{\tilde{y}}$ associated with these values (the maximum of $|w(y)|$ at 17.5 µm is reached only for values situated in the tail of $f_{\tilde{y}}$, that is why the accumulation effect does not happen there). In other respects, there is no value below 2.5 µm and there is very little chance of getting a global effective displacement above 8.75 µm.

<figure><img src="image_12.png" /><figcaption>Fig. 10. Response curve of $|w|$ to an offset of the hole.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Approximate probability density function associated with $|w|$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 12. Approximate cumulative density function associated with $|w|$.</figcaption></figure>

The cumulative distribution function $F_{|w|}$ associated with $|w|$ (Fig. 12) makes it possible to evaluate the probability of exceeding a critical threshold: then only 2% of the values of $|w|$ are above 8.75 µm.

## 5. Conclusion

The VTCR initially introduced to calculate the vibrations of slightly damped elastic plates in the medium-frequency range is a very general approach with a strong mechanical basis. This tool is not only efficient for assemblies of homogeneous plates, but also flexible enough to deal with locally heterogeneous substructures economically. Moreover, using an appropriate formalism, the designer can even take uncertainty into account in the context of a sensitivity analysis.

In the near future, the enriched superelement and the sensitivity analysis procedure will be implemented in the program CORAY-MF for more general applications. Afterwards, the capabilities of the VTCR will be extended from local heterogeneity to other complex cases, such as multi-stiffened substructures.

## References

[1] Babuska I, Ihlenburg F, Paik E, Sauter S. A generalized finite element method for solving the Helmholtz equation in the two dimensions with minimal pollution. Comput Meth Appl Mech Eng 1995;128: 325-59.

[2] Barbone PE, Montgomery JM, Harari I. Scattering by a hybrid asymptotic/finite element method. Comput Meth Appl Mech Eng 1998;164:141-56.

[3] Belov VD, Ryback SA. Applicability of the transport equation in the one-dimensional wave propagation problem. J Soviet Phys Acoust 1975;21(2):110-4.

[4] Boucard P-A, Ladevèze P. A multiple solution method for non-linear structural mechanics. Mech Eng 1999;50(5):317-28.

[5] Buvailo LE, Ionov AV. Application of the finite element method to the investigation of the vibroacoustical characteristics of structures at high audio frequencies. J Soviet Phys Acoust 1980;26(4): 277-9.

[6] Craig RR. Substructure Methods in Vibration. J Vib Acoust 1995;117:207-13 [50th Anniversary].

[7] Cuschieri JM. Vibration transmission through periodic structures using a mobility power flow approach. J Sound Vib 1990;143(1): 65-74.

[8] De Langre E. Fonctions de transfert de plaques en flexion par équation intégrale. Tests de validation et de performance. CEA: DMT/90/395, 1991.

[9] Demkowicz L, Karafiat A, Oden JI. Solution of elastic scattering problems in linear acoustics using h-p boundary element method. Comput Meth Appl Mech Eng 1992;101:251-82.

[10] Faravelli L. Response-surface approach for reliability analysis. J Eng Mech ASCE 1989;115(12):317-28.

[11] Girard A, Defosse H. Frequency response smoothing and structural path analysis: application to beam trusses. J Sound Vib 1993;165(1): 165-170.

[12] Greenstadt J. Solution of wave propagation problems by the cell discretisation method. Comput Meth Appl Mech Eng 1999;174:1-21.

[13] Grosh K, Pinsky PM. Galerkin generalized least square finite element methods for time harmonic structural acoustics. Comput Meth Appl Mech Eng 1998;154:299-318.

[14] Harari I, Grosh K, Hughes TJR, Malhotra M, Pinsky PM, Stewart JR, et al. Recent developments in finite element methods for structural acoustics. Arch Comput Meth. Eng. - State of the art reviews 1996;3:131-309.

[15] Harari I, Haham S. Improved finite element method for elastic waves. Comput Meth Appl Mech Eng 1998;166:143-64.

[16] Harari I, Hughes TJR. Galerkin least-squares finite element method for the reduced wave equation with non reflecting boundary condition in unbounded domains. Comput Meth Appl Mech Eng 1992;98: 411-54.

[17] Ichchou MN, Le Bot A, Jézéquel L. Energy model of one dimensional multi-propagative systems. J Sound Vibration 1997;201(5):535-54.

[18] Kitahara M. Boundary integral equation methods in eigenvalue problems of elastodynamics and thin plates. Elsevier; 1985.

[19] Ladevèze P. A new computational approach for structure vibrations in the medium frequency range. CR Acad Sci Paris, Sér IIb 1996; 322(12):849-56.

[20] Ladevèze P, Arnaud L. A new computational method for structural vibrations in the medium-frequency range. Comput Assist Mech Engineering Sciences 2000;7:219-26.

[21] Ladevèze P, Arnaud L, Rouch P, Blanzé C. The variational theory of complex rays for the calculation of medium-frequency vibrations. Revue Européenne des Eléments Finis 1999;9:67-88.

[22] Ladevèze P, Arnaud L, Rouch P, Blanzé C. The variational theory of complex rays for the calculation of medium-frequency vibrations. Eng Comput 2001;18(1/2):193-214.

[23] Ladevèze P, Blanc L, Rouch P, Blanzé C. A multiscale and computational method for medium-frequency vibrations of assemblies of heterogeneous plates. Comput Struct (Special Issue: Advanced computational models and techniques in dynamics) 2003;81(12): 1267-76.

[24] Leung AYT, Chan JKW. Fourier p-elements for the analysis of beams and plates. J Sound Vib 1998;212(1):179-85.

[25] Liu WK, Zhang Y, Ramirez MR. Multiple scale finite element methods. Int J Numer Meth Eng 1991;32:969-90.

[26] Lyon RH, Maidanik G. Power flow between linearly coupled oscillators. J Acoust Soc of Am 1962;34(5):623-39.

[27] Lyon RH, Richard H, Richard G. Statistical Energy Analysis. Heinemann: Butterworth; 1995.

[28] Morand JPH. A modal hybridization method for the reduction of dynamic models. In: Ladevèze P, Zienkiewicz OC, editors. New Advances in Computational Structural Mechanics. Elsevier; 1992. p. 347-65.

[29] Nefske DJ, Sung SH. Power flow finite element analysis of dynamic systems: basic theory and application to beams. J Vib Stress Reliab Des, Transa of the ASME 1989;111:94-100.

[30] Oberai AA, Pinsky PM. A multiscale finite element method for the Helmholtz equation. Comput Meth Appl Mech Eng 1998;154: 281-98.

[31] Schueller GI, Bucher CG, Pradlwarter HJ. The response surface method, an efficient tool to determine the failure probability of large structural systems. Proc of the Int Conf on Spacecraft Structural and Mechanical Testing, SP-321. Netherlands: ESA, Noordwijk; 1991. p. 247-51.
[32] Soize C. The local effects in the linear dynamic analysis structures in the medium frequency range. In: Ladevèze P, editor. Local effects in the analysis of structures. Barking and Amsterdam: Elsevier; 1985. p. 253–78.

[33] Soize C. Reduced models in the medium frequency range for general dissipative structural-dynamics systems. Eur J Mech, A/Solids 1998;17/4:657–85.

[34] Wu K, Ginsberg JH. Mid frequency range acoustic radiation from slender elastic bodies using the surface variational principle. J Vib Acoust 1998;120:392–400.

[35] Yoon Young K, Jeaong Hoon K. Free vibration analysis of membrane using wave type base functions. JASA 1996;99(5):2938–46.

[36] Zielinski AP, Herrera I. Trefftz method : fitting boundary conditions. Int J Numer Meth Eng 1987;24(5):871–91.