# Analysis Of Structures With Stochastic Interfaces In The Medium-Frequency Range

Claude Blanzé, Philippe Rouch

► To cite this version:

Claude Blanzé, Philippe Rouch. Analysis Of Structures With Stochastic Interfaces In The Medium-Frequency Range. Journal of Computational Acoustics, 2005, 13 (04), pp.711-729. 10.1142/s0218396x05002876 . hal-03177883

HAL Id: hal-03177883

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

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

Distributed under a Creative Commons Attribution 4.0 International License
ANALYSIS OF STRUCTURES WITH STOCHASTIC INTERFACES
IN THE MEDIUM-FREQUENCY RANGE

CLAUDE BLANZÉ\* and PHILIPPE ROUCH†

LMT Cachan — ENS de Cachan/Paris VI/CNRS
61, av. du Pt Wilson — 94235 Cachan Cedex

\*blanze@lmt.ens-cachan.fr

†rouch@lmt.ens-cachan.fr

This paper proposes efficient techniques to obtain effective quantities when dealing with complex structures (including stochastic parameters, such as interface parameters) in medium-frequency vibrations. The first ingredient is the use of a dedicated approach — the Variational Theory of Complex Rays (VTCR) — to solve the medium-frequency problem. The VTCR, which uses two-scale shape functions verifying the dynamic equation and the constitutive relation, can be viewed as a means of expressing the power balance at the different interfaces between substructures. The second ingredient is the use of the Polynomial Chaos Expansion (PCE) to calculate the random response. Since the only uncertain parameters are those which appear in the interface equations (which, in this application, are the complex connection stiffness parameters), this approach leads to very low computation costs.

*Keywords:* Medium frequency vibrations; stochastic interface; polynomial chaos; VTCR.

# 1. Introduction

Today, the main numerical modeling techniques for the analysis of medium-frequency vibrations are all based on finite element (FEM) or boundary element (BEM) approaches. In order to represent small-wavelength phenomena in complex structures (such as car chassis, satellites or ships), these techniques require huge numbers of degrees of freedom, since at least seven elements per wavelength are required to represent oscillating solutions properly. In addition, the solutions obtained are highly sensitive to the material properties and boundary conditions. Moreover, damping is often localized at the structural connections; its measurement is very difficult and very error-prone. In order to obtain the stochastic response of the structure, one must carry out a large number of costly computations. The use of high-frequency approaches, such as Statistical Energy Analysis (SEA)⁹ or any of its improvements, appears unsuitable for medium-frequency vibrations: the vibrational behavior becomes too global and, in general, the coupling loss factor cannot be calculated in a predictive way. This paper proposes efficient techniques to obtain effective quantities when dealing with complex structures involving stochastic parameters, such as interface parameters. The first ingredient of these techniques is the use of a dedicated approach — the
<figure><img src="image_2.png" /><figcaption>Fig. 1. The reference problem.</figcaption></figure>

Variational Theory of Complex Rays (VTCR)⁷,¹¹ — to solve the medium-frequency problem. The VTCR, which uses two-scale shape functions verifying the dynamic equation and the constitutive relation, can be viewed as a means of expressing the power balance at the different interfaces between substructures. The second ingredient is the use of the Polynomial Chaos Expansion (PCE)⁶ to calculate the random response. Since the only uncertain parameters are those which appear in the interface equations (which, in this application, are the complex connection stiffness parameters), this approach leads to very low computation costs. The application presented here concerns a 3D assembly of plates connected by an uncertain elastically bonded joint.

## 2. Basic Aspects of the VTCR

### 2.1. 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 be easily generalized to an assembly of $n$ substructures. As an illustration, let us consider two thin, homogeneous, isotropic and elastic Kirchhoff-Love plates and study the steady-state vibrations of this assembly. Classically, all quantities are defined in the complex domain and depend only on the space quantities: the complex quantity $Q(X)$ is associated with $Q(X, t)$ by the relation $Q(X, t) = Q(X) \cdot \exp(i\omega t)$ where $\omega$ is the fixed angular frequency.

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 the following, with $l = (1, 2)$: prescribed deflection $w_l^d$ on part $\partial_{wd}S_l$ of $\partial S_l$, prescribed slope $w_{nl}^d$ on $\partial_{wd}S_l$, prescribed bending moment $M_l^d$ on $\partial_{Ml}S_l$ and, finally, prescribed 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 $S_{ad}^l$ of the displacement-stress pairs $s_l = (w_l, M_l)$ defined on $S_l$ such that:

$$ s_l \in S_{ad}^l \Leftrightarrow \left\{ \begin{array}{ll} w_l \in \mathcal{U} & (\text{set of finite-energy displacement field } H^1(S)) \\ M_l \in \mathcal{S} & (\text{set of finite-energy stress field } [L^2(S)]^3) \\ \Delta\Delta w_l - k_l^4 w_l = 0 & \text{on } S_l \\ M_l = \frac{2h_l^3}{3}(1 + i\eta_l) K_{PSl} \chi(w_l) & \text{for } l = 1, 2 \\ \text{with } k_l^4 = \frac{3\rho_l\omega^2(1 - \nu_l^2)}{(1 + i\eta_l)E_l h_l^2} & \end{array} \right. \quad (1) $$
where $k_l$ is the wave number and $\rho_l, \nu_l, \eta_l, E_l$ and $2h_l$ designate respectively the density, Poisson's ratio, structural damping coefficient, Young's modulus and thickness of plate $S_l$. $K_{PSl}$ is Hooke's tensor for plane stress and $\mathcal{X}$ is the curvature operator. The fields of $S_{ad}^l$ 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: $s_l \in S_{ad}^l$

$$ \left\{ \begin{array}{lll} w_l = w_l^d & \text{on } \partial_{wd}S_l & w_1 = w_2 \quad \text{on } \Gamma_{12} \\ w_{l,n_l} = w_{n_l}^d & \text{on } \partial_{wd}S_l & w_{1,n_1} = w_{2,n_2} \quad \text{on } \Gamma_{12} \\ M_l n_l = M_l^d & \text{on } \partial_{Md}S_l & M_{1n_1} = M_{2n_2} \quad \text{on } \Gamma_{12} \\ K_l n_l = K_l^d & \text{on } \partial_{Kd}S_l & K_{1n_1} = K_{2n_2} \quad \text{on } \Gamma_{12} \end{array} \right. \quad \text{for } l=1,2 \qquad (2) $$

with $M_l n_l = \underline{n}_l \mathbf{M}_l \underline{n}_l$ and $K_l n_l = \underline{n}_l \underline{div}[\mathbf{M}_l] + (\underline{t}_l \mathbf{M}_l \underline{n}_l)_{,l}$.

The reference problem has at most one solution if:

- Hooke's tensors are positive definite

- $\eta_l \ge 0$.

## 2.2. 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 S_{ad}^l \\ \displaystyle \sum_l (A_{\partial S_l}(\delta s_l, s) - L_{\partial S}(\delta s_l)) + \sum_{\substack{l,m \\ l<m}} C_{\Gamma_{lm}} (\delta s_l, s_l, \delta s_m, s_m) = 0 \qquad (3) \\ \forall \delta s_l \in S_{ad}^l \quad \text{for } l=1,2 \quad m=1,2 \end{array} \right. $$

with the boundary conditions:

$$ A_{\partial S_l}(\delta s_l, s_l) - L_{\partial S_l}(\delta s_l) = \mathcal{R}_e \left\{ \begin{array}{l} -i\Omega \left( - \int_{\partial_{wd}S_l} \delta K_{l n_l} (w_l - w_l^d)^* dL \right. \\ \qquad \left. + \int_{\partial_{wd}S_l} \delta \underline{n}_l \mathbf{M}_l \underline{n}_l (w_{l,n_l} - w_{l,n_l}^d)^* dL \right. \\ \qquad \left. + \int_{\partial_{Md}S_l} (\underline{n}_l \mathbf{M}_l \underline{n}_l - M_l^d) \delta w_{l,n_l}^* dL \right. \\ \qquad \left. - \int_{\partial_{Kd}S_l} (K_l n_l - K_l^d) \delta w_l^* dL \end{array} \right\} \qquad (4) $$
and with the transmission conditions on $\Gamma_{12}$ (for a perfect interface, continuity of the displacement and surfaces forces):

$$ C_{\Sigma_{lm}}(\delta s_l, s_l, \delta s_m, s_m) = \mathcal{R}_e \left\{ -i\Omega \left( \int_{\Gamma_{lm}} \frac{1}{2} [-\delta(K_l n_l - K_m n_m)(w_l - w_m)^* \right. \right. \\ \left. \left. + \delta(\underline{n_l} \mathbf{M}_l \underline{n_l} + \underline{n_m} \mathbf{M}_m \underline{n_m})(w_l, n_l + w_m, n_m)^* \right. \right. \\ \left. \left. + (\underline{n_l} \mathbf{M}_l \underline{n_l} - \underline{n_m} \mathbf{M}_m \underline{n_m})\delta(w_l, n_l - w_m, n_m)^* \right. \right. \\ \left. \left. - (K_l n_l + K_m n_m)\delta(w_l + w_m)^* ] dL \right) \right\}. \quad (5) $$

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

In practice, the VTCR procedure consists of two stages: first, one builds admissible fields; then, one introduces the boundary conditions through the discretization of the variational formulation associated with a set of admissible fields chosen among the elements of $S_{ad}$.

## 2.3. Construction of admissible fields

Let us define a subset of $S_{ad}^l$. The VTCR uses two-scale approximations. For example, in the vicinity of a point $\underline{X}$ of $S_l$, the solution is described as a superposition of modes which can be written as follows:

$$ W_l(\underline{X}, \underline{P}) = e^{i\sqrt{\omega P \cdot X}} U_l(\underline{X}, \underline{P}). \quad (6) $$

The moments are deduced from the constitutive relation:

$$ \mathbf{M}_l(\underline{X}, \underline{P}) = (1 + i\eta_l)\mathbf{K}_{PS}\mathcal{X}(W_l(\underline{X}, \underline{P})) \quad (7) $$

$e^{i\sqrt{\omega P \cdot X}}$ represents the “fast” part of the field, which is defined explicitly. $U(\underline{X}, \underline{P})$ is the “slow” part. $\underline{P}$ denotes a vector which characterizes the local vibration mode and is chosen such that the mode verifies the dynamic equations. There are many possible such choices, which depend on the zone one is particularly interested in; the VTCR distinguishes the interior zone, the edge zone and the corner zone.

For example, let us concentrate on the properties of $n$th-order complex rays for the interior zone. Thus, $U$ includes an $n$th-order polynomial in $\underline{X}$. If the damping factor is small, $U$ can be expressed as:

$$ U_l(\underline{X}, \underline{P}) = e^{\frac{\eta_l}{4}\sqrt{\omega P \cdot X}} \mathbf{V}(\underline{X}, \underline{P}) [a_l(\underline{P})]. \quad (8) $$

The coefficients $a_l(\underline{P})$ of the polynomial $\mathbf{V}(\underline{X}, \underline{P})$ [$a_l(\underline{P})$] are unknowns; they form a column [$a_l(\underline{P})$] of “generalized amplitudes”. Consequently, all unknowns are large-wavelength quantities. The damping of the wave being modeled is represented by the expression $e^{\frac{\eta_l}{4}\sqrt{\omega P \cdot X}}$. These fields can be viewed as vibration modes in an infinite domain with the same mechanical properties as substructure $S_l$; therefore, they belong to $S_{adl}$.
In the following sections, for the sake of simplicity, we will consider the particular case of a 0th-order complex ray and omit the index $l$ of the plate. In that case, the simplified expression of the complex ray is:

$$W(\underline{X}, \underline{P}) = a(\underline{P}) \exp\left(\frac{\eta}{4}\sqrt{\Omega\underline{P}} \cdot \underline{X}\right) \exp(i\sqrt{\Omega\underline{P}} \cdot \underline{X}) = a(\underline{P})w_e(\underline{X}, \underline{P}), \quad (9)$$

$a(\underline{P})$ being the complex magnitude of the elementary complex ray $w_e(\underline{X}, \underline{P})$.

The admissibility relation

$$\Delta\Delta W - k^4 W = 0 \quad \text{on } S \qquad (10)$$

requires:

$$(\underline{P} \cdot \underline{P})^2 = r^4 \quad \text{with} \quad r^4 = \frac{3\rho(1 - \nu^2)}{Eh^2}. \qquad (11)$$

Therefore, the locus of the end of the admissible wave vector $\underline{P}$ for the interior zone is a circle $C$ of radius $r$ (see Fig. 2). This curve is a characteristic of the material. Following the circular path, each direction of the plate is taken into account. Examples of such modes can be seen in Fig. 3. The definitions of corner and edge modes are very similar. For the interior zone, one takes into account all admissible ray directions. The solution $w^{sol}$ is written as:

$$w^{sol}(\underline{X}) = \int_C W(\underline{X}, \underline{P}) \, ds. \qquad (12)$$

<figure><img src="image_3.png" /><figcaption>Fig. 2. Admissible $\underline{P}$ for interior modes.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Admissible interior modes for a plate.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Admissible interior modes for a plate.</figcaption></figure>

The integrals are calculated along curve $C$. Clearly, $w^{sol}$ depends on the generalized amplitude $[a]$ relative to $P$ defined on $C$.

## 2.4. Discretized form of the VTCR

All that is needed in order to derive approximations from the VTCR is the definition of subspace $S_{ad}^{lh}$ from $S_{ad}^{l}$. In practice, one discretizes curve $C$ into $N$ elements which can be of different sizes (on Fig. 4, $C$ is discretized into 8 elements of the same size). The amplitude $a(P)$ is considered to be constant along angular sectors $a_{i}^{h}$. Then:

$$w^h(\underline{X}) = \sum_{i=1}^{n} a_i^h \int_{C_i} w_e(\underline{X}, \underline{P}) \, ds = \sum_{i=1}^{n} a_i^h w_i^h(\underline{X}). \quad (13)$$

The generalized amplitude $a_{i}^{h}$ is associated with the basis function $w_{i}^{h}(\underline{X})$, corresponding to the integral of an elementary ray over an angular sector related to $C_{i}$. Then, the subspace $S_{ad}^{h}$ of the approximations deduced from $S_{ad}$ is defined by:

$$s^h \in S_{ad}^h \Leftrightarrow \begin{cases} w^h(\underline{X}) = \mathbf{W}^h(\underline{X})\mathbf{a}^h \\ \mathbf{M}^h(\underline{X}) = \Lambda^h(\underline{X})\mathbf{a}^h. \end{cases} \quad (14)$$

The same procedure can be applied to the corner and edge zones.

In the case of an assembly of two substructures $S_1$ and $S_2$, the variational formulation and the subspace $S_{ad}^h$ defined by Eq. (14) lead to a system of equations of finite dimension in the complex domain:

$$\mathbf{K}\mathbf{a} = \mathbf{f} \quad (15)$$

$\mathbf{K} = \mathbf{E}_{\mathrm{D}}^{h} + \mathbf{Z}_{\mathrm{h}}$, where $\mathbf{E}_{\mathrm{D}}^{h}$ is the symmetric, positive definite damping matrix associated with the dissipated power, and $\mathbf{Z}_{\mathrm{h}}$ is the matrix associated with a bilinear form defined
<figure><img src="image_6.png" /><figcaption>Fig. 5. Generic assembly modeling.</figcaption></figure>

such that $Z_h^* = -Z_h$. $f$ is the vector associated with the linear form and $a$ is the vector of the unknown amplitudes.

## 2.5. Effective quantities

From a mechanical point of view, the spatial distribution of the solution has no physical meaning in the medium-frequency range. Indeed, the response at a discrete location and at a discrete frequency is extremely sensitive to perturbations. In order to interpret the response, spatially-averaged quantities are more useful. Thus, it is necessary to extract from the solution some effective quantities defined on a domain greater than one wavelength. On a given substructure $S$, the quantities involved are expressed as $q(X)$, which can represent, for example, the elastic energy, the kinetic energy or the dissipation work; in the vicinity of any test point, the associated effective quantities are defined on a domain $D$ by:

$$q_D = \frac{1}{D} \int_D |q(X)| d\underline{X}. \qquad (16)$$

## 3. Modeling of the Interfaces

Today, the modeling of beams, plates, shells and their combination no longer presents difficulties for the VTCR, but the connections among these structures can be improved in order to take into account more realistic types of behavior. For example, in real structures such as cars, satellites or launchers, most of the damping is localized in the different joints connecting the substructures. These joints, whether welded, riveted or bolted, modify the structure's response and can be used by engineers in order to tune the structure. Since the wavelengths of the vibrational phenomena in medium-frequency are relatively large compared to the distance between two rivets or two weld points, we decided to introduce these by using distributed interface parameters. These interface parameters can be mass, damping and stiffness coefficients, and they can be identified on generic assemblies.
Generally, identification tests are carried out on several specimens, leading to probabilistic models in order to take into account the variability of the observed parameters (e.g. localized damping). Since this variability can have a huge influence on the solution, one must be able to predict this influence.

### 3.1. The interface model

This local behavior can be modeled by interface relations which represent elastic and/or dissipative joint parameters. Each interface is a mechanical entity whose specific behavior depends on the type of connection. Many different connection types can be modeled by this approach:

$$ -n_1 \underline{\mathbf{M}}_1 \underline{n}_1 + n_2 \underline{\mathbf{M}}_2 \underline{n}_2 = 0 \quad \text{Equilibrium} \qquad (17) $$

$$ \frac{1}{2}(n_1 \underline{\mathbf{M}}_1 \underline{n}_1 + n_2 \underline{\mathbf{M}}_2 \underline{n}_2) = -k_c(1 + i\eta_c)(w_{1,n_1} + w_{2,n_2}) \quad \text{Behavior} \qquad (18) $$

where $k_c$ is the bending stiffness coefficient and $\eta_c$ the damping coefficient.

The part of the variational formulation associated with Eqs. (17) and (18) is done by:

$$ \mathcal{R}_e \left\{ -i\Omega \left( \int_{\Gamma} \left[ \frac{1}{2}(n_1 \underline{\mathbf{M}}_1 \underline{n}_1 + n_2 \underline{\mathbf{M}}_2 \underline{n}_2) + k_c(1 + i\eta_c)(w_{1,n_1} + w_{2,n_2}) \right] (\delta w_{1,n_1} + \delta w_{2,n_2})^* dL \right. \\ \left. + \int_{\Gamma} \frac{1}{2}[n_1 \underline{\mathbf{M}}_1 \underline{n}_1 - n_2 \underline{\mathbf{M}}_2 \underline{n}_2](\delta w_{1,n_1} - \delta w_{2,n_2})^* dL \right) \right\}. \qquad (19) $$

Here, we are dealing with problems of structural assemblies in which the component's behavior and the external loads are deterministic: the randomness of the response comes from the random behavior of the connections. The system's parameters, calibrated using experimental data, are modeled as random variables or processes which are assumed to be properly represented by a set of random variables {$\xi(\theta)$}, where $\theta$ belongs to the space of random events $\Omega$. More precisely, the bending stiffness coefficient $k_c$, which may be uncertain, is modeled by a random variable.

### 3.2. Discretization in random space: Polynomial Chaos Expansion (PCE)

The stochastic interface parameters lead to a discretized linear system similar to Eq. (15):

$$ \mathbf{K}(\theta) \mathbf{a}(\theta) = \mathbf{f}. \qquad (20) $$

The generalized variables $\mathbf{a}(\theta)$ can be expressed formally as a nonlinear functional of the set {$\xi_j(\theta)$} used to represent the material's stochastic property. It has been shown² that this functional dependence can be expanded in terms of polynomial chaos. Then, the truncated PCE of the response takes the form:

$$ \mathbf{a}(\theta) = \sum_{i=0}^{P} \mathbf{a}_i \Psi_i(\theta) \qquad (21) $$

where {$\Psi_i(\theta)$} are polynomials in the Gaussian random variables {$\xi_i$}, which can be shown to form a complete basis of the Hilbert space of second-order random variables. The number
of polynomials P depends on the order p of the PCE and on the number L of stochastic parameters. Let us assume that the material parameters are constant along an interface. If this were not the case, a Karhunen–Loeve expansion could easily be used to represent the spatial randomness of the interface’s characteristics.⁵

The random character of a joint’s parameter α(θ) is made explicit by its argument θ. The stochastic parameter α(θ) is represented by

$$ \alpha(\theta) = \bar{\alpha}(1 + \delta\xi(\theta)) \qquad (22) $$

where $\bar{\alpha}$ is the mathematical expectation of α(θ), δ the coefficient of variation (standard deviation divided by the expectation) and ξ(θ) a standard normal random variable: ξ(θ) ~ N(0, 1). For non-Gaussian material properties, the PCE is used to represent the stochastic material’s parameters: the case of a Gaussian process is a particular case which was chosen for the sake of simplicity.⁴

Following the assembly procedure, this leads to the corresponding expansion of the matrix $\mathbf{K}(\theta)$:

$$ \mathbf{K}(\theta) = \sum_{i=0}^{L} \xi_i(\theta) \mathbf{K}_i \qquad (23) $$

where $\mathbf{K}_0$ denotes the matrix corresponding to the mean material properties and the other terms are related to the random fluctuations about this mean. The number of stochastic parameters is L (by convention, $\xi_0 = 1$). Expanding the generalized parameters $\mathbf{a}(\theta)$ with respect to the polynomial chaos basis, one has:

$$ \mathbf{a}(\theta) = \sum_{j=0}^{P} \mathbf{a}_j \Psi_j(\theta). \qquad (24) $$

Then, substituting Eqs. (23) and (24) into Eq. (20) leads to:

$$ \sum_{j=0}^{P} \sum_{i=0}^{L} \xi_i \Psi_j(\theta) \mathbf{K}_i \mathbf{a}_j = \mathbf{f}. \qquad (25) $$

An equality (in a weak sense) can be derived by projecting Eq. (25) onto the subspace spanned by the polynomial chaos subset used in the approximation; this process results in the following equations where $\langle \cdot \rangle$ represent the mathematical expectation:

$$ \sum_{j=0}^{P} \sum_{i=0}^{L} \langle \xi_i \Psi_j(\theta) \Psi_k(\theta) \rangle \mathbf{K}_i \mathbf{a}_j = \langle \Psi_k(\theta) \rangle \mathbf{f} \quad k = 0, 1, \dots, P. \qquad (26) $$

This system can be rewritten as:

$$ \sum_{j=0}^{P} \sum_{i=0}^{L} c_{ijk} \mathbf{K}_i \mathbf{a}_j = \delta_{0k} \mathbf{f} \quad k = 0, 1, \dots, P, \qquad (27) $$
where the coefficients $c_{ijk}$ denote $\langle \xi_i \Psi_j(\theta) \Psi_k(\theta) \rangle$ and need to be calculated only once. This system of linear equations must be solved for the unknown $\mathbf{a}_j$ of the PCE. The details of this procedure were published in Ref. 6 and its implementation issues addressed in a number of other references.⁴,⁵ These equations can be assembled into a matrix of size $(P+1)n \times (P+1)n$ ($n$ being the number of degrees of freedom) of the form:

$$ \begin{pmatrix} \mathbf{K}^{(00)} & \mathbf{K}^{(0k)} & \mathbf{K}^{(0P)} \\ \mathbf{K}^{(j0)} & \mathbf{K}^{(jk)} & \mathbf{K}^{(jP)} \\ \mathbf{K}^{(P0)} & \mathbf{K}^{(Pk)} & \mathbf{K}^{(PP)} \end{pmatrix} \begin{pmatrix} \mathbf{a}_0 \\ \mathbf{a}_k \\ \mathbf{a}_P \end{pmatrix} = \begin{pmatrix} \mathbf{f} \\ \mathbf{0} \\ \mathbf{0} \end{pmatrix} \quad (28) $$

where

$$ \mathbf{K}^{(jk)} = \sum_{i=0}^{L} \langle \xi_i \Psi_j(\theta) \Psi_k(\theta) \rangle \mathbf{K}_i. \quad (29) $$

In summary, this approach consists in expanding the random response process about a basis of the Hilbert space of random variables and calculating the coefficients of this expansion. The result is a convergent expansion of the response in terms of multidimensional orthogonal polynomials. Although the methodology used is becoming widespread, serious obstacles have been encountered from a computational point of view in practical implementations. In large realistic problems, the methodology is either cumbersome or computationally intensive. Some numerical strategies, such as iterative algorithms, have been devised to overcome the numerical difficulties which arise in this context.¹⁰ Our proposed approach does not present such drawbacks: it uses very small elementary matrices, thus leading to very small problems, and since the evolution of the generalized amplitude is relatively smooth it can be easily represented by the Hermitian polynomials of the chaos expansion.

# 4. Numerical Examples

## 4.1. Deterministic problem: Comparison with NASTRAN

In order to verify the efficiency of the VTCR for medium frequency problems, we compared the results obtained with our approach to those obtained with a classical FE method. The FE results were obtained with the MSC NASTRAN program using ten elements per wavelength. The benchmark structure represents a highly simplified car door made of 3D steel plate assemblies. The reference problem is described in Fig. 6.

All edges of the structure are simply-supported, except for the three edges P1P3, P3P5 and P5P7 along which a distributed shear force of 1N/m at 400 Hz is applied. Figure 7 shows the two solutions (displacement amplitude).

The VTCR solution was calculated using about 700 DOFs, while the FE method required about 15 000 DOFs. The two solutions were plotted on the same scale and the maximum displacements were $2.52 \cdot 10^{-6}$ m for VTCR and $2.82 \cdot 10^{-6}$ m for NASTRAN. This difference
<figure><img src="image_7.png" /><figcaption>Fig. 6. The reference problem.</figcaption></figure>

<table><thead><tr><th colspan="3">Points location (dimensions in meters):</th></tr></thead><tbody><tr><td>p1</td><td>= 0.08</td><td>0.08</td><td>-0.24</td></tr><tr><td>p2</td><td>= 1.0</td><td>0.08</td><td>-0.24</td></tr><tr><td>p3</td><td>= 0.</td><td>0.</td><td>0.</td></tr><tr><td>p4</td><td>= 0.96</td><td>0.</td><td>0.</td></tr><tr><td>p5</td><td>= 0.</td><td>0.</td><td>0.2</td></tr><tr><td>p6</td><td>= 0.96</td><td>0.</td><td>0.2</td></tr><tr><td>p7</td><td>= 0.08</td><td>0.05</td><td>0.33</td></tr><tr><td>p8</td><td>= 0.2</td><td>0.05</td><td>0.33</td></tr><tr><td>p9</td><td>= 0.88</td><td>0.05</td><td>0.33</td></tr><tr><td>p10</td><td>= 0.96</td><td>0.05</td><td>0.33</td></tr><tr><td>p11</td><td>= 0.15</td><td>0.05</td><td>0.4</td></tr><tr><td>p12</td><td>= 0.2</td><td>0.05</td><td>0.4</td></tr><tr><td>p13</td><td>= 0.56</td><td>0.17093</td><td>0.8</td></tr><tr><td>p14</td><td>= 0.96</td><td>0.17093</td><td>0.823</td></tr><tr><td>p15</td><td>= 0.54</td><td>0.18</td><td>0.83</td></tr><tr><td>p16</td><td>= 1.</td><td>0.18</td><td>0.86</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Fig. 7. VTCR solution (left) NASTRAN solution (right).</figcaption></figure>

is insignificant: the two methods yield the same displacement distribution and the FE maximum level is quite sensitive to the number of DOFs in the calculation. As is the case of any Trefftz method, the VTCR allows one to use very small numbers of DOFs which remain constant as the frequency increases.¹¹ Conversely, in FE methods, $k^3h^2$ should remain constant,³ where k is the wave number and h the distance between two nodes. Gain factors of up to 100 can be easily obtained. But, as in all Trefftz methods, the elementary matrices of the VTCR are full, which reduces the CPU gains to between 5 and 10.⁸

## 4.2. Assembly of two plates with an uncertain interface

### 4.2.1. *Description*

Figure 8 describes the reference problem involving two plates connected by a complex flexural spring. The imaginary part of the stiffness is linked directly to the local loss factor of that joint (Eq. (17)).
<figure><img src="image_9.png" /><figcaption>Fig. 8. The reference problem and the associated solution for the parametric joint.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Influence of the joint's parameters.</figcaption></figure>

These local parameters can have a huge influence on the solution. In order to have an idea of this influence, we plotted the variation of the average displacement level of the two plates against the local stiffness or damping (in either case, the other parameter remains constant).

### 4.2.2. Probabilistic approach

The joint's parameters are now modeled using random variables in order to introduce the uncertainties observed in identification tests (the local damping factor, for example, may vary by 1000 percent depending on the tightening pressure).

$$k_R(\theta) = \bar{k}_R(1 + \delta\xi(\theta)) \quad (30)$$

For this problem with a single random variable, a polynomial chaos of order 3 was used. This means that $P = 3$ and, therefore, 4 functions were used in the expansion (Table 1).

Figure 10 shows the four coefficients of the polynomial chaos expansion of order 3. The first coefficient represents the mean displacement of the response $\mathbf{a}_0$ and the others ($\mathbf{a}_1, \mathbf{a}_2$ and $\mathbf{a}_3$) the main fluctuations about the mean.
<table><caption>Table 1. One-dimensional polynomial chaoses and their variances ($L = 1$).</caption><thead><tr><th>$i$</th><th>$\Psi_i$</th><th>$\langle\Psi_i^2\rangle$</th><th>$p$</th></tr></thead><tbody><tr><td>0</td><td>1</td><td>1</td><td>0</td></tr><tr><td>1</td><td>$\xi$</td><td>1</td><td>1</td></tr><tr><td>2</td><td>$-1 + \xi^2$</td><td>2</td><td>2</td></tr><tr><td>3</td><td>$\xi(-3 + \xi^2)$</td><td>6</td><td>3</td></tr><tr><td>4</td><td>$3 - 6\xi^2 + \xi^4$</td><td>24</td><td>4</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>Fig. 10. The four coefficients (0, 1, 2, 3 from left to right).</figcaption></figure>

In order to compare the parametric resolution to the VTCR method with polynomial chaos, we show in Fig. 11 the displacement at the end of the second plate as a function of $\xi$ for different orders of polynomial chaos. As expected, the PCE yields a very good approximation of the exact solution. The accuracy increases with the order of the PCE and the third order is sufficient in terms of smoothness of the solution.
<figure><img src="image_12.png" /><figcaption>Fig. 11. Comparison between the parametric and chaos solutions.</figcaption></figure>

## 4.3. Longeron: Two uncertain interfaces

### 4.3.1. Description

The next example represents an assembly of four plates including two uncertain interfaces on edges $P_5P_6$ and $P_7P_8$. A distributed moment of $1$Nm/m is applied on edge $P_5P_8$ and the geometric and mechanical properties are defined in Fig. 12. The mean of the random variable associated with the two joints is the same ($\bar{k}_R = 5000$ Nm/rad), even though these two joints are independent.

The deterministic solution is presented in Fig. 13.

### 4.3.2. Probabilistic approach

In order to test the capabilities of the VTCR, the polynomial decomposition was conducted using second-order and third-order polynomial functions. The second-order decomposition led to the computation of six coefficients associated with six Hermite's functions (Table 2) while third order decomposition will need ten coefficients (Table 2). In this specific case,
<figure><img src="image_13.png" /><figcaption>Fig. 12. Characteristics of the longeron.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Deterministic solution at 2200 Hz.</figcaption></figure>

<table><caption>Table 2. Two-dimensional polyno-
mial chaos and their variances
($L = 2$).</caption><thead><tr><th><i>i</i></th><th>&Psi;<sub><i>i</i></sub></th><th>&lang;&Psi;<sup>2</sup><sub><i>i</i></sub>&rang;</th><th><i>p</i></th></tr></thead><tbody><tr><td>0</td><td>1</td><td>1</td><td>0</td></tr><tr><td>1</td><td>&xi;<sub>1</sub></td><td>1</td><td>1</td></tr><tr><td>2</td><td>&xi;<sub>2</sub></td><td>1</td><td></td></tr><tr><td>3</td><td>-1 + &xi;<sub>1</sub><sup>2</sup></td><td>2</td><td></td></tr><tr><td>4</td><td>&xi;<sub>1</sub> &xi;<sub>2</sub></td><td>1</td><td>2</td></tr><tr><td>5</td><td>-1 + &xi;<sub>2</sub><sup>2</sup></td><td>2</td><td></td></tr><tr><td>6</td><td>&xi;<sub>1</sub> (-3 + &xi;<sub>1</sub><sup>2</sup>)</td><td>6</td><td></td></tr><tr><td>7</td><td>(-1 + &xi;<sub>1</sub><sup>2</sup>) &xi;<sub>2</sub></td><td>2</td><td></td></tr><tr><td>8</td><td>&xi;<sub>1</sub> (-1 + &xi;<sub>2</sub><sup>2</sup>)</td><td>2</td><td>3</td></tr><tr><td>9</td><td>&xi;<sub>2</sub> (-3 + &xi;<sub>2</sub><sup>2</sup>)</td><td>6</td><td></td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>Fig. 14. The six coefficients associated with the polynomial functions.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. The ten coefficients associated with the polynomial functions.</figcaption></figure>
there is no significant difference between the six coefficients of the second-order calculation and the first six coefficients of the third-order calculation (Figs. 14 and 15).

### 4.3.3. Comparison between the chaos and parametric solutions

The graph (Fig. 16) shows the comparison between the chaos and the parametric solutions for the third-order decomposition. The normalized value $\tilde{a}_1$ of the first element $a_1(\xi_1, \xi_2)$ is represented for an excitation frequency of 2200 Hz: $\tilde{a}_1(\xi_1, \xi_2) = (a_1(\xi_1, \xi_2) - a_1(0, 0))/a_1(0, 0)$.

Figure 17 shows the comparison between the exact parametric solution and the polynomial chaos representation for the normalized amplitude $\tilde{a}_1$, as functions of $\xi_1 = \xi_2$ and for different orders of polynomial chaos.

<figure><img src="image_17.png" /><figcaption>Fig. 16. Exact parametric solution (left) and chaos representation (right).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. Comparison between the parametric solution (dots) and the chaos solution (lines).</figcaption></figure>
**5. Conclusion**

The Variational Theory of Complex Rays is a general approach to medium-frequency calcu-
lations. Its main features are the use of shape functions with a sound mechanical meaning
and a dedicated variational formulation allowing the use of independent approximations in
the different substructures. In some cases, such as assemblies of beams, the VTCR leads
to the exact solution in tension as well as in bending. For complex structures, the VTCR
seems to be a very promising theory due to its ability to give predictive results at a very
low numerical cost. It can also take into account complex connections among substructures,
described by stochastic parameters, in order to represent the physical reality at the lowest
possible numerical cost. In this case, the benefit of two-scale shape functions is twofold: first,
the global system defined by the PCE is relatively small; second, the smooth variation of the
unknowns (generalized amplitudes) when the variation of the connection parameters can
be easily represented by a polynomial function. Extensions of this approach to structural
heterogeneities,¹ transient dynamics and frequency bandwidth⁸ are in progress.

**Appendix**

A. Basis Functions of the Polynomial Chaos Expansion

Table 1 shows the expressions of the basis functions of the PCE for the orders $p = 0$ to
$p = 4$ when only one random variable $\xi$ is used ($L = 1$). Table 2 shows the expressions of
the basis functions of the PCE for the orders $p = 0$ to $p = 3$ when two random variables $\xi_1$
and $\xi_2$ are used ($L = 2$).

**References**

1. L. Blanc, C. Blanzé, P. Ladevèze and P. Rouch, A multiscale and trefftz computational method for medium-frequency vibrations of assemblies of heterogeneous plates, *Computer-Assisted Mechanics and Engineering Sciences* **10** (2003) 375–384.

2. R. H. Cameron and W. T. Martin, The orthogonal development of nonlinear functionals in series of Fourier–Hermite functionals, *Ann. Math.* **48**(16) (1947) 385–392.

3. A. Deraemaeker, I. Babuska and P. Bouillard, Dispersion and pollution of the fem solution for the Helmholtz equation in one, two and three dimensions, *Int. J. Numer. Meth. Engng.* **46** (1999) 471–499.

4. R. Ghanem, Ingredients for a general purpose stochastic finite elements implementation, *Com-put. Methods Appl. Mech. Engrg.* **168** (1999) 19–34.

5. R. Ghanem and R. Kruger, Numerical solution of spectral stochastic finite element systems, *Comput. Methods Appl. Mech. Engrg.* **129** (1996) 289–303.

6. R. G. Ghanem and P. D. Spanos, *Stochastic Finite Elements: A Spectral Approach* (Springer, Berlin, 1991).

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

8. P. Ladevèze, P. Rouch, H. Riou and X. Bohineust, Analysis of medium-frequency vibrations in a frequency range, *J. Comp. Acoustic* **11**(2) (2003) 255–283.
9. R. H. Lyon and G. Maidanick, Power flow between linearly coupled oscillators, JASA 34(5) (1967) 623-639.

10. M. Pellissetti and R. Ghanem, Iterative solution of systems of linear equations arising in the context of the stochastic fem, *J. Advances in Engineering Software* **31** (2000) 607-616.

11. P. Rouch and P. Ladevèze, The variational theory of complex rays: A predictive tool for medium-frequency vibrations, *Comp. Meth. Appl. Mech. Engrg.* **192** (2003) 3301-3315.