# The Variational Theory of Complex Rays applied to the shallow shell theory

Alessandro Cattabiani, Hervé Riou, Andrea Barbarulo, Pierre Ladevèze, Guillaume Bézier, Bernard Troclet

► To cite this version:

Alessandro Cattabiani, Hervé Riou, Andrea Barbarulo, Pierre Ladevèze, Guillaume Bézier, et al. The Variational Theory of Complex Rays applied to the shallow shell theory. Computers & Structures, 2015, 158, pp.98-107. 10.1016/j.compstruc.2015.05.021. hal-01647848

HAL Id: hal-01647848

https://hal.science/hal-01647848

Submitted on 8 Oct 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 applied to the shallow shell theory

Alessandro Cattabiani<sup>a</sup>, Hervé Riou<sup>a</sup>, Andrea Barbarulo<sup>a</sup>, Pierre Ladevèze<sup>a,*</sup>, Guillaume Bézier<sup>b</sup>,
Bernard Troclet<sup>a,c</sup>

In the last few decades the interest of aerospace and automotive industries towards the study of the medium-frequency response of complex shell structure frames has grown. Recently some dedicated "wave" computational approaches have been developed. Among them, a Trefftz technique called Variational Theory of Complex Rays (VTCR) is catching on as an *ad hoc* method to deal with such vibration problems. This work presents the development of the VTCR in the shallow shell theory to increase its effectiveness and flexibility. First, general theory is given and two key properties of the solution demonstrated. After that, two numerical examples are deeply analyzed.

## 1. Introduction

In recent years, the interest of aerospace and automotive industries has been focused on efficient virtual testing of the vibration response. Shallow shell structures are widely used in these industrial contests due to their high resistance and light weight. The equilibrium equations of shallow shells are quite complex and in almost every real case an analytic solution cannot be obtained. Thus, an effective method to predict vibrational behavior in shallow shell structures is needed. The Modal Overlap Factor [1] defines three zones: low, mid and high frequency range. The low-frequency range has been extensively studied by the Finite Element Method (FEM) [2] and the Boundary Element Method (BEM) [3]. On the other side, the high-frequency range can be addressed by the Statistical Energy Analysis (SEA) [1]. This technique neglects almost entirely spatial quantities to focus on global energy. This effective approach is based on some key assumptions assured in the high-frequency range. The medium-frequency range is still an open question. On one hand the FEM and BEM are not indicated in this frequency domain since the phenomena variation length is very small if compared to characteristic dimensions of the structure. For this reason the required number of Degrees of

Freedom (DoFs) explodes [4]. On the other hand the SEA is not suggested because the key assumptions of the theory might be unsatisfied [5]. Notwithstanding a lot of work has been done to extend such theories to the medium-frequency range [6–8]. There are also methods developed specifically for the medium-frequency range such as the partition of unity method [9], the ultra-weak variational method [10], the asymptotic scaled modal analysis [11], the energy operator eigenmodes [12], Galerkin method [13], the wave boundary element method [14] or the wave-based method [15,16]. One of them is the Variational Theory of Complex Rays (VTCR). It approximates the vibrational problem solution with a sum of shape functions that identically satisfies the equilibrium equations while addressing the boundary conditions in weak form. This approach allows a priori independent approximations among subdomains. Thus, different (in number and type) shape functions can be chosen for each subdomain giving great flexibility to the method. It has already been applied to plate theory [17], to general shell theory [18], to transient dynamics [19], to 3D acoustic [20] and, on a wide frequency band [21,22]. Nevertheless the shell version of the VTCR can still be improved. Yet the in-plane inertia was not taken into account in previous works, the weak variational formulation must be customized for the specific geometry, and the VTCR formulation does not address the general case of a boundary (or a corner) shared by multiple subdomains. Such problems are analyzed and solved in this work; the in-plane inertia assumption is relaxed and two propagative waves that lead in-plane stresses and displacements are introduced. The customization phase of

<footnote><sup>a</sup> LMT-Cachan (ENS Cachan/CNRS/Paris-Saclay-University), 61 avenue du Président Wilson, F-94230 Cachan, France</footnote>

<footnote><sup>b</sup> CNES Direction des Lanceurs, 52 rue Jacques Hillairet, 75612 Paris Cedex, France</footnote>

<footnote><sup>c</sup> Airbus Defence and Space, 66 Route de Verneuil, BP3002, 78133 Les Mureaux Cedex, France</footnote>

<footnote>* Corresponding author. Tel.: +33 (0) 147402402; fax: +33 (0) 147402240.

E-mail addresses: cattabiani@lmt.ens-cachan.fr (A. Cattabiani), riou@lmt.ens-cachan.fr (H. Riou), barbarulo@lmt.ens-cachan.fr (A. Barbarulo), ladeveze@lmt.ens-cachan.fr (P. Ladevèze), guillaume.bezier@cnes.fr (G. Bézier), troclet@lmt.ens-cachan.fr (B. Troclet).</footnote>
the weak variational formulation is avoided using the shallow shell approximation providing effectiveness and flexibility to the method. Since in this theory the surface geometry is projected to the underlying area, the tuning phase of the weak variational formulation is no more needed. Finally, a more general version of the VTCR is presented.

The present work is structured as follows: first, the general theory is proposed providing some useful properties. After, two numerical examples are presented. The first one is an academic case where the analytic solution is known. Convergence tests are performed and performances are compared with a FEM reference. The second one is a complex structure frame.

## 2. Theory

In this Section the equilibrium and the boundary equations are examined using the standard shallow shell approximations. The theory is akin to the one provided in [23,24]. After, some useful energy quantities are derived and the virtual work theorem is adapted for this specific case.

### 2.1. The equilibrium equations

The general reference example is presented in Fig. 1. The focus is on a generic subdomain $\Omega_i$ of the frame structure in Fig. 1. For the sake of clarity various boundary, corner, coupling, and surface conditions are split in Fig. 2. The $\partial\square\Omega_i$ symbol refers to a generic boundary of $\Omega_i$ where condition $\square$ is applied. In the particular case of a boundary shared among subdomains, $\Gamma$ is used instead. In the same way, for the conditions applied on corners, a symbol $\partial\square\Omega_i$ is used. The generic corner shared among subdomains is indicated with $C$. The over-line symbol $\square$ indicates that quantity $\square$ is known (i.e. the value of the boundary constraint). The term $\hat{n}_i$ is the normal unit vector of the boundary directed outward $\Omega_i$. A subdomain is subject to loads, displacements constraints, and continuity conditions along the boundaries (Fig. 2a) and on the corners (Fig. 2b) as well as a distributed load per unit surface $\bar{g}_i$ (Fig. 2d). The displacement constraint $\bar{u}_i = [\bar{v}'_i, \bar{w}_i]'$ along $\partial\bar{u}, \Omega_i$ can be divided in in-plane $\bar{v}_i$ and out-of-plane $\bar{w}_i$ components.<fnref n="1" /> In the same way, the load per unit length $\bar{p}_i = [\bar{b}'_i, \bar{q}_i]'$ along $\partial\bar{p}, \Omega_i$ can be divided in in-plane $\bar{b}_i$ and out-of-plane $\bar{q}_i$ components. The rotation condition $\bar{w}_{i,n_i}$ is imposed along $\partial\bar{w}_{i,n_i}, \Omega_i$ while a bending moment per unit length $\bar{m}_i$ is applied along $\partial\bar{m}_i, \Omega_i$. The corners of the subdomain are subject to out-of-plane displacements constraints $\bar{w}_{ci}$ on $\partial\partial\bar{w}_{ci}, \Omega_i$ and punctual forces $\bar{q}_{ci}$ on $\partial\partial\bar{q}_{ci}, \Omega_i$. Coupling conditions are applied on $C$ and along $\Gamma$, in order to ensure continuity of stresses and displacements among subdomains (Fig. 2c).

All quantities of interest are defined in the complex domain. Each one is considered multiplied by $e^{j\omega t}$ where $j = \sqrt{-1}$ is the imaginary unit, $\omega = 2\pi f$ is the angular frequency and $t$ is the time.

The geometry of the subdomain can be approximated by its projection on the local plane defined by the orthonormal basis $\{\hat{x}_i, \hat{y}_i, \hat{z}_i\}$ and the displacement field can be restricted to (Kirchhoff's kinematics assumptions)

$$ 
\begin{aligned}
\mathbf{u}_i^z &= \mathbf{u}_i - z_i \phi_i \\
\phi_i &= \nabla w_i - \mathbf{R}_i \cdot \mathbf{v}_i \\
\mathbf{R}_i &= \begin{bmatrix} \frac{1}{R_{xi}} & 0 \\ 0 & \frac{1}{R_{yi}} \end{bmatrix}
\end{aligned}
$$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Generic frame structure described in Section 2.1.</figcaption></figure>

where $\mathbf{u}_i^z$ is the displacement thorough the thickness of the shell, $\mathbf{u}_i, \mathbf{v}_i$ and $w_i$ are respectively the total, the in-plane and the out-of-plane displacements of the middle surface and $\mathbf{R}_i$ is the curvature matrix. $D_i = \{\mathbf{u}_i, \mathbf{N}_i, \mathbf{M}_i\}$ is the set of fields that satisfies the equilibrium equations

$$ \mathbf{u}_i^z \in \mathcal{U}_i^z \quad \text{finite energy displacement set,} \quad (1) $$

$$ \{\mathbf{N}_i, \mathbf{M}_i\} \in S_i \quad \text{finite energy generalized stress set,} \quad (2) $$

$$ \nabla \cdot \mathbf{N}_i + \bar{g}_{xyi} + \varrho_i \hbar_i \omega^2 \mathbf{v}_i = 0 \quad \text{over } \Omega_i, \quad (3) $$

$$ \nabla \cdot (\nabla \cdot \mathbf{M}_i) + \mathbf{R}_i : \mathbf{N}_i + \bar{g}_{zi} + \varrho_i \hbar_i \omega^2 \mathbf{w}_i = 0 \quad \text{over } \Omega_i, \quad (4) $$

$$ \bar{g}_i = [\bar{g}_{xi}, \bar{g}_{yi}, \bar{g}_{zi}]', \quad (5) $$

$$ \bar{g}_{xyi} = [\bar{g}_{xi}, \bar{g}_{yi}]', \quad (6) $$

$$ \mathbf{M}_i = -D_i \mathbf{H}_i : \nabla \nabla w_i, \quad (7) $$

$$ \mathbf{N}_i = \frac{12}{h_i^2} D_i \mathbf{H}_i : (\mathbf{E}_i - \mathbf{R}_i w_i), \quad (8) $$

$$ \mathbf{E}_i = [\nabla \mathbf{v}_i]_{\text{sym}} = \frac{1}{2} (\nabla \mathbf{v}_i + \nabla \mathbf{v}_i'), \quad (9) $$

$$ D_i = \frac{E_i h_i^3}{12(1 - v_i^2)}, \quad (10) $$

$$ E_i = E_{0i} (1 + j \eta_i), \quad (11) $$

where $D_i \mathbf{H}_i$ is the Hooke's operator for plane stress, $\varrho_i$ the density, $\eta_i$ the damping coefficient, $h_i$ the shell thickness, $E_{0i}$ the Young modulus, $v_i$ the Poisson's ratio, $\square: \square$ the inner product matrix operator, and $\mathbf{N}_i$ and $\mathbf{M}_i$ are the stress and stress moment resultants tensors respectively. The sub-space of $D_i$ associated with homogenized conditions ($\bar{g}_i = 0$) is denoted as $D_{0i} = \{\delta\mathbf{u}_i, \delta\mathbf{N}_i, \delta\mathbf{M}_i\}$. This definition will be useful in the next Sections.

### 2.2. The boundary conditions

In order to present a well-posed problem, three conditions imposed along each boundary and one on each corner are needed. The boundary and corner conditions presented in Fig. 1 can be classified in this way:

1. an in-plane condition, either a displacement constraint or a load per unit length ($\bar{v}_i$ or $\bar{b}_i$),

2. an out-of-plane condition, either a displacement constraint or a load per unit length ($\bar{w}_i$ or $\bar{q}_i$),

3. either a rotation or a bending moment per unit length ($\bar{w}_{i,n_i}$ or $\bar{m}_i$),

4. an out-of-plane condition, on corners either a displacement constraint or a punctual load ($\bar{w}_{ci}$ or $\bar{q}_{ci}$).

These boundary conditions can be imposed along edges and on corners shared among subdomains. On one hand, if the constraint is a displacement or a rotation and there is more than one subdomain involved there is no need of the relative coupling condition. In other words the subdomains can be considered decoupled for what concerns that particular constraint. On the other hand, if the condition is a load and there is more than one subdomain

<footnote n="1">¹ $\square'$ is the transpose operator.</footnote>
<figure><img src="image_3.png" /></figure>

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

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

<figure><img src="image_6.png" /><figcaption>Fig. 2. Boundary, corner, coupling and surface constraints of the generic structure frame of Fig. 1 subdivided and highlighted.</figcaption></figure>

involved, the coupling conditions hold. The boundary conditions in the most general case are

$$ \mathbf{v}_i = \begin{cases} \bar{\mathbf{v}}_i & \nexists \bar{\mathbf{v}}_i \\ \mathbf{v}_{j_r} & \nexists \bar{\mathbf{v}}_i \end{cases}, \qquad (12) $$

$$ w_i = \begin{cases} \bar{w}_i & \exists \bar{w}_i \\ w_{j_r} & \nexists \bar{w}_i \end{cases}, \qquad (13) $$

$$ w_{i, \hat{n}_i} = (\nabla w_i) \hat{n}_i = \begin{cases} \bar{w}_{i, \hat{n}_i} & \exists \bar{w}_{i, \hat{n}_i} \\ w_{j_r, \hat{n}_i} & \nexists \bar{w}_{i, \hat{n}_i} \end{cases}, \qquad (14) $$

$$ \mathbf{b}_i = \mathbf{N}_i \hat{n}_i = \bar{\mathbf{b}}_i - \sum_{j_r=1}^{n_r} \mathbf{b}_{j_r}, \qquad (15) $$

$$ \begin{aligned} q_i &= (\nabla \cdot \mathbf{M}_i) \hat{n}_i + \nabla (\hat{\mathbf{t}}_i \mathbf{M}_i \hat{n}_i) \cdot \hat{\mathbf{t}}_i \\ &= (\nabla \cdot \mathbf{M}_i) \hat{n}_i + (\hat{\mathbf{t}}_i \mathbf{M}_i \hat{n}_i)_{\hat{\mathbf{t}}_i} = \bar{q}_i - \sum_{j_r=1}^{n_r} q_{j_r}, \end{aligned} \qquad (16) $$

$$ m_i = \hat{\mathbf{n}}_i \mathbf{M}_i \hat{\mathbf{n}}_i = \bar{m}_i - \sum_{j_r=1}^{n_r} m_{j_r}, \qquad (17) $$

where $\hat{\mathbf{t}}_i$ is the tangent unit vector, $j_r$ is the index relative to the other subdomains sharing with $\Omega_i$ the boundary $\Gamma$, and $n_r$ is their total number ($\Omega_i$ excluded) as is shown in Fig. 3.

The possible conditions on corners in the most general case are

$$ w_{ci} = \begin{cases} \bar{w}_{ci} & \exists \bar{w}_{ci} \\ w_{cj_c} & \nexists \bar{w}_{ci} \end{cases}, \qquad (18) $$

$$ q_{ci} = \hat{\mathbf{t}}_{1i} \mathbf{M}_i \hat{\mathbf{n}}_{1i} + \hat{\mathbf{t}}_{2i} \mathbf{M}_i \hat{\mathbf{n}}_{2i} = \bar{q}_{ci} - \sum_{j_c=1}^{n_c} q_{j_c}, \qquad (19) $$

<figure><img src="image_7.png" /><figcaption>Fig. 3. $n_r + 1$ subdomains ($\Omega_i$ included) sharing the same boundary $\Gamma$.</figcaption></figure>

where $\hat{\mathbf{n}}_{1i}$ and $\hat{\mathbf{n}}_{2i}$ are the outward normal unit vectors of the two boundaries of $\Omega_i$ that share the corner $C$. $\hat{\mathbf{t}}_{1i}$ and $\hat{\mathbf{t}}_{2i}$ are their respective tangent unit vectors directed toward the corner. As for $j_r$, the index $j_c$ is relative to other subdomains sharing with $\Omega_i$ the corner $C$. $n_c$ is their total number ($\Omega_i$ excluded) as is shown in Fig. 4.

## 2.3. Energies

In the shallow shell theory and in case of forced vibrations some useful energy quantities can be defined. These are the strain energy $E_S$, the kinetic energy $E_K$, and the dissipation energy $E_D$.

$$ E_S = \frac{1}{2} \sum_{i=1}^{n} \left\{ \int_{\Omega_i} \mathbf{N}_i^H : (\mathbf{E}_i - \mathbf{R}_i w_i) - \mathbf{M}_i^H : \nabla \nabla w_i \, d\mathbf{s} \right\}, \qquad (20) $$

$$ E_K = \frac{1}{2} \sum_{i=1}^{n} \left\{ \int_{\Omega_i} \rho_i h_i \omega^2 \mathbf{u}_i^H \cdot \mathbf{u}_i \, d\mathbf{s} \right\}, \qquad (21) $$

$$ E_D = \mathfrak{F}\{E_S\}. \qquad (22) $$

where $\mathfrak{F}$ is the hermitian operator. The operator $\mathfrak{F}\{\square\}$ extracts the imaginary part of $\square$. These will be extensively used in the following Sections.

## 2.4. Virtual work theorem

This Section starts enunciating the virtual work theorem for shallow shells subject to forced vibrations. After some rearrangements a useful version of the theorem is illustrated. In general, the standard virtual work theorem is

<figure><img src="image_8.png" /><figcaption>Fig. 4. $n_c + 1$ subdomains ($\Omega_i$ included) that share the same corner $C$.</figcaption></figure>
$$ \mathcal{L}_i = \mathcal{L}_e \qquad (23) $$

where $\mathcal{L}_i$ and $\mathcal{L}_e$ are the internal and the external virtual works respectively. In case of forced vibrations in the shallow shell theory it can be developed by means of the previously defined energies

$$ 2(E_s - E_K) = \mathcal{L}_e. \qquad (24) $$

The quantity $E_s - E_K$ can be rearranged through residues along boundaries and on corners. This equivalence lies on a the divergence vector calculus identity. It can be specialized in two slightly different forms

$$ \mathbf{V} \cdot (\mathbf{B}\mathbf{a}) = (\mathbf{V} \cdot \mathbf{B})\mathbf{a} + \mathbf{B} : \mathbf{V}\mathbf{a}, \qquad (25) $$

$$ \mathbf{V} \cdot (\varphi\mathbf{a}) = \varphi(\mathbf{V} \cdot \mathbf{a}) + \mathbf{a}\mathbf{V}\varphi, \qquad (26) $$

where $\mathbf{B}, \mathbf{a}$ and $\varphi$ are a generic matrix, vector, and scalar respectively. Such identities can be coupled with the divergence theorem to obtain the following equalities

$$ \int_{\Omega_i} \mathbf{B} : \mathbf{V}\mathbf{a}\mathbf{d}s = \oint_{\partial\Omega_i} (\mathbf{B}\mathbf{a})\hat{\mathbf{n}}\mathbf{d}s - \int_{\Omega_i} (\mathbf{V} \cdot \mathbf{B})\mathbf{a}\mathbf{d}s, \qquad (27) $$

$$ \int_{\Omega_i} \mathbf{a}\mathbf{V}\varphi\mathbf{d}s = \oint_{\partial\Omega_i} (\varphi\mathbf{a})\hat{\mathbf{n}}\mathbf{d}s - \int_{\Omega_i} \varphi(\mathbf{V} \cdot \mathbf{a})\mathbf{d}s. \qquad (28) $$

the strain energy $E_s$ can be developed using these two equations. In particular, the part related to in-plane stresses and displacements becomes

$$ \int_{\Omega_i} \mathbf{N}_i^H : \mathbf{E}_i \mathbf{d}s = \oint_{\partial\Omega_i} \mathbf{b}_i^H \mathbf{v}_i \mathbf{d}s - \int_{\Omega_i} (\mathbf{V} \cdot \mathbf{N}_i)^H \mathbf{v}_i \mathbf{d}s, \qquad (29) $$

where the following equality is implicitly used

$$ \mathbf{N}_i^H : \mathbf{E}_i = \mathbf{N}_i^H : \mathbf{V}\mathbf{v}_i \qquad (30) $$

which holds since $\mathbf{N}_i$ is symmetric by definition. The last term can be substituted using the in-plane equilibrium Eq. (3)

$$ \int_{\Omega_i} \mathbf{N}_i^H : \mathbf{E}_i \mathbf{d}s = \oint_{\partial\Omega_i} \mathbf{b}_i^H \mathbf{v}_i \mathbf{d}s + \int_{\Omega_i} \varphi_i h_i \omega^2 \mathbf{v}_i^H \mathbf{v}_i \mathbf{d}s + \int_{\partial\Omega_i} \bar{\mathbf{g}}_{\mathbf{x}y_i}^H \mathbf{v}_i \mathbf{d}s. \qquad (31) $$

In the same way the part related to the out-of-plane stresses and displacements of the strain energy $E_s$ using Eq. (27) becomes

$$ - \int_{\Omega_i} \mathbf{M}_i^H : \mathbf{V}\mathbf{V}\mathbf{w}_i \mathbf{d}s = - \oint_{\partial\Omega_i} (\mathbf{M}_i \hat{\mathbf{n}}_i)^H \mathbf{V}\mathbf{w}_i \mathbf{d}s + \int_{\Omega_i} (\mathbf{V} \cdot \mathbf{M}_i)^H \mathbf{V}\mathbf{w}_i \mathbf{d}s. \qquad (32) $$

The last term can be further developed introducing the identity expressed in Eq. (28)

$$ \int_{\Omega_i} (\mathbf{V} \cdot \mathbf{M}_i)^H \mathbf{V}\mathbf{w}_i \mathbf{d}s = \oint_{\partial\Omega_i} [(\mathbf{V} \cdot \mathbf{M}_i)\hat{\mathbf{n}}_i]^H \mathbf{w}_i \mathbf{d}s - \int_{\Omega_i} [\mathbf{V} \cdot (\mathbf{V} \cdot \mathbf{M}_i)]^H \mathbf{w}_i \mathbf{d}s. \qquad (33) $$

Using Eq. (4)

$$ - \int_{\Omega_i} \mathbf{M}_i^H : \mathbf{V}\mathbf{V}\mathbf{w}_i \mathbf{d}s = \oint_{\partial\Omega_i} [(\mathbf{V} \cdot \mathbf{M}_i)\hat{\mathbf{n}}_i]^H \mathbf{w}_i - (\mathbf{M}_i\hat{\mathbf{n}}_i)^H \mathbf{V}\mathbf{w}_i \mathbf{d}s \\ + \int_{\Omega_i} [\mathbf{R}_i : \mathbf{N}_i]^H \mathbf{w}_i + \varphi_i h_i \omega^2 \mathbf{w}_i^H \mathbf{w}_i + \bar{\mathbf{g}}_{zi}^H \mathbf{w}_i \mathbf{d}s. \qquad (34) $$

The first term of the right side relative to the integral along the boundary can be further developed. In fact, by definition

$$ \oint_{\partial\Omega_i} [(\mathbf{V} \cdot \mathbf{M}_i)\hat{\mathbf{n}}_i]^H \mathbf{w}_i - (\mathbf{M}_i\hat{\mathbf{n}}_i)^H \mathbf{V}\mathbf{w}_i \mathbf{d}s = - \oint_{\partial\Omega_i} m_{sdi}^H w_{sdi,\hat{\mathbf{n}}_i} \mathbf{d}s + \oint_{\partial\Omega_i} [(\mathbf{V} \cdot \mathbf{M}_{sdi})\hat{\mathbf{n}}_i]^H w_{sdi} - (\hat{\mathbf{t}}_i \mathbf{M}_{sdi} \hat{\mathbf{n}}_i)^H w_{sdi,\hat{\mathbf{n}}_i} \mathbf{d}s. \qquad (35) $$

The term $-\oint_{\partial\Omega_i} (\hat{\mathbf{t}}_i \mathbf{M}_{sdi} \hat{\mathbf{n}}_i)^H w_{sdi,\hat{\mathbf{n}}_i} \mathbf{d}s$ can be integrated by parts

$$ - \oint_{\partial\Omega_i} (\hat{\mathbf{t}}_i \mathbf{M}_{sdi} \hat{\mathbf{n}}_i)^H w_{sdi,\hat{\mathbf{n}}_i} \mathbf{d}s = - \left[ (\hat{\mathbf{t}}_i \mathbf{M}_{sdi} \hat{\mathbf{n}}_i)^H w_{sdi} \right]_{\hat{\mathbf{n}}_i}^{\mathbf{s}=0} + \oint_{\partial\Omega_i} (\hat{\mathbf{t}}_i \mathbf{M}_{sdi} \hat{\mathbf{n}}_i)^H w_{sdi} \mathbf{d}s. \qquad (36) $$

The term $[(\hat{\mathbf{t}}_i/\mathbf{M}_{sdi} \hat{\mathbf{n}}_i)^H w_{sdi}]_{\hat{\mathbf{n}}_i}$ is a continuous function along boundaries that jumps on corners. Thus:

$$ - \left[ (\hat{\mathbf{t}}_i \mathbf{M}_{sdi} \hat{\mathbf{n}}_i)^H w_{sdi} \right]_{\hat{\mathbf{n}}_i}^{\mathbf{s}=0} = -q_{cADI}^H w_{cADI} \Big|_{\partial\partial\Omega_i}. \qquad (37) $$

Back-substituting all these developments into Eq. (24) and simplifying leads to

$$ \mathfrak{Z} \left\{ \sum_{i=1}^{n} \left\{ \oint_{\partial\Omega_i} \mathbf{b}_i^H \mathbf{v}_i + q_i^H w_i - m_i^H w_{i,\hat{\mathbf{n}}_i} \mathbf{d}s - q_{cCi}^H w_{cCi} \Big|_{\partial\partial\Omega_i} \right\} \right\} = \mathcal{L}_{0e} \qquad (38) $$

where the left hand side of Eq. (38) are the boundary residues and $\mathcal{L}_{0e}$ is the virtual work of external forces without considering the distributed loads $\bar{\mathbf{g}}_i$. In fact, it is erased by the simplification with the right hand side of the equation. This development of the virtual work theorem leads to an equation that links boundary residues to dissipation energy $E_D$. Since $E_K$ is real valued by definition, the following equation holds

$$ \mathfrak{Z} \left\{ \sum_{i=1}^{n} \left\{ \oint_{\partial\Omega_i} \mathbf{b}_i^H \mathbf{v}_i + q_i^H w_i - m_i^H w_{i,\hat{\mathbf{n}}_i} \mathbf{d}s - q_{cCi}^H w_{cCi} \Big|_{\partial\partial\Omega_i} \right\} \right\} = 2E_D - \mathfrak{Z}\{\mathcal{L}_{ge}\} \qquad (39) $$

where $\mathcal{L}_{ge}$ is the virtual work of distributed loads $\bar{\mathbf{g}}_i$.

### 3. Shallow shell VTCR

In this Section the shallow shell version of the VTCR technique is illustrated. Since the VTCR is a Trefftz method, the solution is researched in the set of functions that satisfy the interior equilibrium equations. Boundary and corner residuals are addressed using Eq. (39). The weak variational problem definition is:

find the solution set $D_{si} = \{\mathbf{u}_{si}, \mathbf{N}_{si}, \mathbf{M}_{si}\} \in \mathcal{D}$ where $i \in [1, \dots, n]$ is the index related to the subdomain $\Omega_i$ such that

$$ \sum_{i=1}^{n} \left\{ \begin{aligned} & \int_{\partial \bar{\mathbf{v}}_i \Omega_i} \delta \mathbf{b}_i^H (\mathbf{v}_{si} - \bar{\mathbf{v}}_i) \mathbf{d}s - \int_{\partial \bar{\mathbf{b}}_i \Omega_i} \delta \mathbf{v}_i^H (\mathbf{b}_{si} - \bar{\mathbf{b}}_i) \mathbf{d}s \\ & + \int_{\partial \bar{\mathbf{w}}_i \Omega_i} \delta q_i^H (w_{si} - \bar{w}_i) \mathbf{d}s - \int_{\partial \bar{\mathbf{q}}_i \Omega_i} \delta w_i^H (q_{si} - \bar{q}_i) \mathbf{d}s \\ & - \int_{\partial \bar{\mathbf{w}}_{i, \hat{\mathbf{n}}_i} \Omega_i} \delta m_i^H (w_{si, \hat{\mathbf{n}}_i} - \bar{w}_{i, \hat{\mathbf{n}}_i}) \mathbf{d}s \\ & + \int_{\partial \bar{\mathbf{m}}_i \Omega_i} \delta w_{i, \hat{\mathbf{n}}_i}^H (m_{si} - \bar{m}_i) \mathbf{d}s + \delta w_{cCi}^H (q_{csi} - \bar{q}_{cCi}) \Big|_{\partial \bar{\mathbf{q}}_{cCi} \Omega_i} \\ & - \delta q_{cCi}^H (w_{cCi} - \bar{w}_{cCi}) \Big|_{\partial \bar{\mathbf{q}}_{cCi} \Omega_i} + \frac{n_{\Gamma}}{n_{\Gamma} + 1} \int_{\Gamma} \delta \mathbf{b}_i^H \mathbf{v}_{sj} + \delta q_i^H w_{sj} + \delta m_i^H w_{sj, \hat{\mathbf{n}}_j} + \delta \mathbf{v}_i^H \mathbf{b}_{sj} \\ & - \frac{n_c}{n_c + 1} \delta q_{cCi}^H w_{cCi} \Big|_c - \frac{1}{n_{\Gamma} + 1} \int_{\Gamma} \delta \mathbf{v}_i^H \mathbf{b}_{si} + \delta w_i^H q_{si} \\ & - \delta w_{i, \hat{\mathbf{n}}_i}^H m_{si} \mathbf{d}s + \frac{1}{n_c + 1} \delta w_{cCi}^H q_{cCi} \Big|_c \\ & - \sum_{j=1, j \neq i}^{n} \left\{ \frac{1}{n_{\Gamma} + 1} \int_{\Gamma} \delta \mathbf{b}_i^H \mathbf{v}_{sj} + \delta q_i^H w_{sj} + \delta m_i^H w_{sj, \hat{\mathbf{n}}_j} + \delta \mathbf{v}_i^H \mathbf{b}_{sj} + \delta w_i^H q_{sj} + \delta w_{i, \hat{\mathbf{n}}_i}^H m_{sj} \mathbf{d}s - \frac{1}{n_c + 1} (\delta q_{cCi}^H w_{cCi} + \delta w_{cCi}^H q_{cCi}) \Big|_c \right\} \end{aligned} \right\} = 0 \forall \{\delta \mathbf{u}_i, \delta \mathbf{N}_i, \delta \mathbf{M}_i\} \in \mathcal{D}_{0i}. \qquad (40) $$
where $D_{sd0}$ is the test functions space being the VTCR a Galerkin method.

## 3.1. Fundamental properties

In this Subsection uniqueness and existence properties are demonstrated.

### 3.1.1. Uniqueness

This property is demonstrated by contraction. Assume that two different solutions $D_{s1} = \{\mathbf{u}_{s1}, \mathbf{N}_{s1}, \mathbf{M}_{s1}\}$ and $D_{s2} = \{\mathbf{u}_{s2}, \mathbf{N}_{s2}, \mathbf{M}_{s2}\}$ exist over the whole domain of the shallow shell problem and consider their difference $D_{sd} = \{\mathbf{u}_{sd}, \mathbf{N}_{sd}, \mathbf{M}_{sd}\}$. Since the equilibrium equations are all linear, $D_{sd}$ is solution of the difference problem. This means that it satisfies the homogeneous equilibrium equations, the constraints along the boundaries and on the corners are all zero and $D_{sd0} = D_{sd}$. In this case

$$ \sum_{i=1}^{n} \left\{ \begin{aligned} & \int_{\partial \bar{\omega}_{sdi}} \Omega_i \delta \mathbf{b}_{sdi}^H \mathbf{v}_{sdi} \mathrm{d}\mathbf{s} - \int_{\partial \bar{\omega}_{sdi}} \Omega_i \delta \mathbf{v}_{sdi}^H \mathbf{b}_{sdi} \mathrm{d}\mathbf{s} + \int_{\partial \bar{\omega}_{sdi}} \Omega_i \delta q_{sdi}^H \mathbf{w}_{sdi} \mathrm{d}\mathbf{s} \\ & - \int_{\partial \bar{\omega}_{sdi}} \Omega_i \delta w_{sdi}^H q_{sdi} \mathrm{d}\mathbf{s} - \int_{\partial \bar{\omega}_{sdi}} \Omega_i \delta \dot{m}_{sdi}^H \mathbf{w}_{sdi.\mathbf{n}_i} \mathrm{d}\mathbf{s} + \int_{\partial \bar{\omega}_{sdi}} \Omega_i \delta w_{sdi.\mathbf{n}_i}^H m_{sdi} \mathrm{d}\mathbf{s} \\ & + \delta w_{Csd}^H q_{Csd} \Big|_{\partial \Omega_{Csd}} \Omega_i - \delta q_{Csd}^H w_{Csd} \Big|_{\partial \Omega_{Csd}} \Omega_i + \frac{n_F}{n_F + 1} \int_F \delta \mathbf{b}_{sdi}^H \mathbf{v}_{sdi} \\ & + \delta q_{sdi}^H \mathbf{w}_{sdi} - \delta m_{sdi}^H \mathbf{w}_{sdi.\mathbf{n}_i} \mathrm{d}\mathbf{s} - \frac{n_C}{n_C + 1} \delta q_{sdi}^H \mathbf{w}_{Csd} \Big|_C \\ & - \frac{1}{n_F + 1} \int_F \delta \mathbf{v}_{sdi}^H \mathbf{b}_{sdi} + \delta w_{sdi}^H q_{sdi} - \delta w_{sdi.\mathbf{n}_i}^H m_{sdi} \mathrm{d}\mathbf{s} + \frac{1}{n_C + 1} \delta w_{Csd}^H q_{Csd} \Big|_C \end{aligned} \right\} \\ = 0 \{\delta \mathbf{u}_{sdi}, \delta \mathbf{N}_{sdi}, \delta \mathbf{M}_{sdi}\} \in D_{sdi} \quad (41) $$

where $i$ is the index of the generic subdivision $\Omega_i$ of the domain $\Omega$. Since this must be true for every test function, the weak form can be rewritten, without loss of generality, such that

$$ \begin{aligned} \Re \left\{ \sum_{i=1}^{n} \left\{ \begin{aligned} & \int_{\partial \bar{\omega}_{sdi}} \Omega_i \delta \mathbf{b}_{sdi}^H \mathbf{v}_{sdi} \mathrm{d}\mathbf{s} - \int_{\partial \bar{\omega}_{sdi}} \Omega_i \delta \mathbf{v}_{sdi}^H \mathbf{b}_{sdi} \mathrm{d}\mathbf{s} + \int_{\partial \bar{\omega}_{sdi}} \Omega_i \delta q_{sdi}^H \mathbf{w}_{sdi} \mathrm{d}\mathbf{s} \\ & - \int_{\partial \bar{\omega}_{sdi}} \Omega_i \delta w_{sdi}^H q_{sdi} \mathrm{d}\mathbf{s} - \int_{\partial \bar{\omega}_{sdi}} \Omega_i \delta \dot{m}_{sdi}^H \mathbf{w}_{sdi.\mathbf{n}_i} \mathrm{d}\mathbf{s} + \int_{\partial \bar{\omega}_{sdi}} \Omega_i \delta w_{sdi.\mathbf{n}_i}^H m_{sdi} \mathrm{d}\mathbf{s} \\ & + \delta w_{Csd}^H q_{Csd} \Big|_{\partial \Omega_{Csd}} \Omega_i - \delta q_{Csd}^H w_{Csd} \Big|_{\partial \Omega_{Csd}} \Omega_i + \frac{n_F}{n_F + 1} \int_F \delta \mathbf{b}_{sdi}^H \mathbf{v}_{sdi} \\ & + \delta q_{sdi}^H \mathbf{w}_{sdi} - \delta m_{sdi}^H \mathbf{w}_{sdi.\mathbf{n}_i} \mathrm{d}\mathbf{s} - \frac{n_C}{n_C + 1} \delta q_{Csd}^H w_{Csd} \Big|_C \\ & + \delta w_{sdi}^H q_{sdi} - \delta w_{sdi.\mathbf{n}_i}^H m_{sdi} \mathrm{d}\mathbf{s} + \frac{1}{n_C + 1} \delta w_{Csd}^H q_{Csd} \Big|_C \end{aligned} \right\} \right\} \\ = 0 \{\delta \mathbf{u}_{sdi}, \delta \mathbf{N}_{sdi}, \delta \mathbf{M}_{sdi}\} \in D_{sdi}. \quad (42) $$

Considering $\Re\{\Box_1^H \Box_2\} = \Re\{-\Box_2^H \Box_1\}$, the weak form (42) can be rearranged such that

$$ \Re \left\{ \sum_{i=1}^{n} \left\{ \oint_{\partial \Omega_i} \delta \mathbf{b}_{sdi}^H \mathbf{v}_{sdi} + \delta q_{sdi}^H \mathbf{w}_{sdi} - \delta m_{sdi}^H \mathbf{w}_{sdi.\mathbf{n}_i} \mathrm{d}\mathbf{s} - \delta q_{Csd}^H w_{Csd} \Big|_{\partial \Omega_i} \right\} \right\} \\ = 0 \{\delta \mathbf{u}_{sdi}, \delta \mathbf{N}_{sdi}, \delta \mathbf{M}_{sdi}\} \in D_{sdi}. \quad (43) $$

Substituting Eq. (39) and keeping in mind that in this particular case $L_{ge} = 0$, leads to the conclusion that the weak variational formulation implies that

$$ E_{sdD} = 0. \quad (44) $$

Since the dissipation energy is positive-definite, the deformations $\mathbf{N}_{sd}$ and $\mathbf{M}_{sd}$ are null. Back-substituting this result in Eqs. (3) and

(4) (with $\bar{\mathbf{g}}_{sd} = \mathbf{0}$ because this is the difference solution) implies that also the total displacement of the difference solution $\mathbf{u}_{sd}$ must be zero. Thus the two solutions coincide and the uniqueness property is proven.

### 3.1.2. Existence

If the solution of the forced vibration problem respects the regularity conditions in Eqs. (1) and (2), it is also solution of the weak variational formulation. If $D_{si}$ is a discrete space approximation, the uniqueness properties leads also to the existence property.

## 3.2. Discretization

Different sets of shape functions can be chosen. In [18] plane waves have been chosen, meanwhile in [20] Fourier functions has been selected. In this work plane waves are used since integrals over straight boundaries can be calculated analytically reducing computational costs. Their generic form is

$$ f_{SF} = a_i \hat{c}_i e^{\mathbf{k}_i \cdot \mathbf{x}_{rel}} \quad (45) $$

where $\mathbf{k} = k_i \hat{\mathbf{k}}_i$ is the wave vector, $\hat{\mathbf{k}}_i$ is the unit direction vector of the wave vector and $k_i$ is a complex scalar term, $a_i$ is the amplitude of the plane wave, $\hat{c}_i$ is the unit direction of the plane wave, and $\mathbf{x}_{rel}$ is the relative position vector. The relative coordinate system is located in the geometric center<fnref n="2" /> of the element and $\hat{\mathbf{z}}$ is directed along the normal surface. The unknown parameters are $\hat{c}_i, k_i, a_i$. The first two are chosen so that equilibrium equations are identically satisfied meanwhile the amplitude is set by the weak variational formulation. The unit direction vector $\hat{\mathbf{k}}_i$ is set a priori. It is the only discretized parameter of the VTCR. This class of shape functions can be further divided in propagative and evanescent waves. These are both needed as explained in [18]. The difference between them lies in the unit direction vector $\hat{\mathbf{k}}_i$ that, in case of the propagative waves, is

$$ \hat{\mathbf{k}}_i = \mathbf{T} \begin{bmatrix} 1 \\ 0 \end{bmatrix} \quad (46) $$

where

$$ \mathbf{T} = \begin{bmatrix} \cos(\theta_i) & -\sin(\theta_i) \\ \sin(\theta_i) & \cos(\theta_i) \end{bmatrix} \quad (47) $$

is the rotation matrix and $\theta_i \in [0; 2\pi)$ is the discretized angle of the plane wave direction. In case of evanescent waves the unit direction vector is defined such that

$$ \hat{\mathbf{k}}_i = \mathbf{T} \cdot \begin{bmatrix} \cosh(\phi_i) \\ j \sinh(\phi_i) \end{bmatrix} \quad (48) $$

where $\phi_i \in [-\log(1 + \sqrt{2}); \log(1 + \sqrt{2}))$ is another discretized angle in the complex domain. Indeed out-of bound values of $\phi_i$ can be chosen. Nevertheless a discretized range must be chosen so, as criteria, a span where the imaginary part ranges from $-1$ to $1$ and the real part from $0$ to $\sqrt{2}$ is selected.

Substituting the generic shape function defined in Eq. (45) with a given unit direction vector $\hat{\mathbf{k}}_i$ in Eqs. (3) and (4) leads to a linear set of equilibrium equations that in matrix form is

$$ \mathbf{Z}_i \mathbf{c}_i a_i = 0 \quad (49) $$

where $\mathbf{Z}_i = \mathbf{Z}_i(k_i)$ is a $3 \times 3$ matrix that depends on $k_i$. In order get not trivial results ($a_i = 0$) the dispersion Equation must be imposed

$$ \det[\mathbf{Z}_i] = 0 \quad (50) $$

<footnote n="2">² The geometric center is chosen just for the sake of convenience. Any other point could fulfill the same role.</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 5. Geometric dimensions and boundary conditions of the curved shallow shell described in Section 4.1.</figcaption></figure>

It provides the values of $k_i$ for the specified unit direction vector $\hat{k}_i$. When $k_i$ is extracted from Eq. (50), the term $\hat{c}_i$ can be obtained by re-injecting $k_i$ in Eq. (50) and enforcing that

$$ \|\hat{c}_i\| = 1. \qquad (51) $$

## 4. Numerical examples

Two tests are presented in this Section. The first one studies a very simple vibrational problem where the analytic solution is known. The second is a complex frame structure where the VTCR solution is compared with a FEM one.

### 4.1. First numerical example: curved shallow shell

This academic example is conducted to demonstrate the effectiveness of the method in the medium-frequency range and to analyze convergence rates of the VTCR. In this particular case the “quasi-analytic” solution is known. This reference solution [23] neglects the in-plane inertia. This approximation is acceptable as illustrated in [23].

The considered structure is a curved shallow shell. Geometric dimensions and boundary conditions are reported in Fig. 5. Other important quantities are shown in Table 1.

<table><caption>**Table 1**

Quantities of interest of the example in Section 4.1.</caption><tr><td><b>x</b></td><td>Generic relative position</td><td>[x<br>y]</td></tr><tr><td><i>h</i></td><td>Thickness</td><td>3 mm</td></tr><tr><td><i>R</i><sub>x</sub></td><td>Radius of curvature along <b>x</b></td><td>2.0304 m</td></tr><tr><td><i>R</i><sub>y</sub></td><td>Radius of curvature along <b>y</b></td><td>+∞</td></tr><tr><td><i>f</i></td><td>Frequency</td><td>2000 Hz</td></tr><tr><td><i>E</i></td><td>Young modulus</td><td>200 GPa</td></tr><tr><td><i>ν</i></td><td>Poisson's ratio</td><td>0.3</td></tr><tr><td><i>ρ</i></td><td>Density</td><td>7800 kg/m<sup>3</sup></td></tr><tr><td><i>η</i></td><td>Damping factor</td><td>0.01</td></tr><tr><td><b>p</b></td><td>Relative position of the punctual load</td><td>[0.4<br>0.4]<br>m</td></tr><tr><td><b>g</b></td><td>Punctual load</td><td>[0<br>0<br>1]<br>δ(x - p)e<sup>j2πf</sup></td></tr></table>

The curved shallow shell is partially simply supported and is subject to a time-dependent punctual load directed along $\hat{z}$. The reference solution approximates it by a series of sines

<figure><img src="image_10.png" /><figcaption>(a) Reference solution.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) VTCR solution.</figcaption></figure>

<caption>**Fig. 6.** Comparison of the out-of-plane displacements of the vibrational problem explained in Section 4.1.</caption>
$$ 
\begin{aligned}
g_z &= e^{j\omega t} \sum_{m=1}^{\infty} \sum_{n=1}^{\infty} C_{mn} \sin(\lambda_n x) \sin(\mu_m x) \\
C_{mn} &= \frac{4}{L_x L_y} \int_{L_x} \int_{L_y} \delta(\mathbf{x} - \mathbf{p}) \sin(\lambda_n x) \sin(\mu_m x) dy dx
\end{aligned} 
\quad (52) $$

$$ \lambda_n = \frac{n\pi}{L_x} $$

$$ \mu_m = \frac{m\pi}{L_y} $$

where $\omega = 2\pi f$, $g_z$ is the load along $\hat{\mathbf{z}}$, and $\delta(\mathbf{x} - \mathbf{p})$ is the Dirac delta function. The out-of-plane displacement $w$ of the analytic solution is approximated by a Fourier series

$$ 
\begin{aligned}
g_z &= e^{j\omega t} \sum_{m=1}^{\infty} \sum_{n=1}^{\infty} B_{mn} \sin(\lambda_n x) \sin(\mu_m x) \\
B_{mn} &= \left[ \frac{Eh(\frac{\lambda_n^2}{R_x} + \frac{\mu_m^2}{R_y})^2}{(\lambda_n^2 + \mu_m^2)^2} + D(\lambda_n^2 + \mu_m^2)^2 - \rho h \omega^2 \right]^{-1} C_{mn}.
\end{aligned} 
\quad (53) $$

<figure><img src="image_12.png" /><figcaption>Fig. 7. p-iteration. Convergence rate of the VTMR with the propagative rays number.</figcaption></figure>

The series is stopped at {$n_{max} = 100, m_{max} = 100$} because convergence is reached. Its graphical representation is given in Fig. 6 alongside the VTMR solution.

In the VTMR the punctual force is approximated by a complete Fourier approximation:

$$ 
\begin{aligned}
g_z &= e^{j\omega t} \sum_{m=-m_{VTCR_{max}}}^{m_{VTCR_{max}}} \sum_{n=-n_{VTCR_{max}}}^{n_{VTCR_{max}}} C_{VTCR_{mn}} e^{j(\lambda_{VTCR_{n}}x + \mu_{VTCR_{m}}x)} \\
C_{VTCR_{mn}} &= \frac{1}{L^2} \int_{L}^{L} \int_{L}^{L} \delta(\mathbf{x} - \mathbf{p}) e^{j(\lambda_{VTCR_{n}}x + \mu_{VTCR_{m}}x)} dy dx \\
\lambda_{VTCR_{n}} &= \frac{n\pi}{L} \\
\mu_{VTCR_{m}} &= \frac{m\pi}{L} \\
L &= \sqrt{L_x^2 + L_y^2}.
\end{aligned} 
\quad (54) $$

Injecting it in Eqs. (3) and (4) and erasing the in-plane inertia produces the particular solution of the VTMR. Since the Fourier approximation can tackle almost every possible distributed load, the construction of the particular solution is very general and can be applied straightforwardly to any other load.

The error of the converged VTMR is led by the approximation of the particular solution. Kirchhoff theory allows errors of the order of magnitude $\frac{h}{R_x} = 10^{-3}$. Considering this error threshold, the series is stopped at {$n_{VTCR_{max}} = 40, m_{VTCR_{max}} = 40$}. The error indicator is based on the kinetic energy of the out-of-plane displacement

$$ E_{Kz}(w) = \frac{1}{2} \rho h \omega^2 \int_{S} w^H \cdot w ds \quad (55) $$

$$ \text{err} = \frac{E_{Kz}(w_{VTMR} - w_{ref})}{E_{Kz}(w_{ref})} \quad (56) $$

where $w_{ref} = w_{ref}(x, y)$ is the out-of-plane displacement of the reference solution and $w_{VTMR} = w_{VTMR}(x, y)$ is the out-of-plane displacement of the VTMR solution.

<figure><img src="image_13.png" /><figcaption>Fig. 8. p-iteration and p-h-iteration. This set of figures summarizes the results of the p-h iteration. The domain is divided in 1; 2; 4; 8 subdomains and the number of rays is increased up to convergence. With the increasing of subdomains, convergence is reached with less rays per subdomain but more total rays.</figcaption></figure>
For the sake of simplicity the number of evanescent waves per boundary and per element are kept equal to half the number of propagative rays. This is not a binding choice. If the convergence rate is not satisfactory this ratio can be changed. Nevertheless, in order to present results in a simple and clean way (without comparing them with too many variables) a characteristic parameter has to be chosen. A number of propagative waves is selected for this purpose and the number of evanescent waves is linked to it. The total number of rays is

$$n_{r_{tot}} = n_e \left( n_{rp} + n_r \frac{n_{rp}}{2} + n_r + n_c \right) \qquad (57)$$

where $n_{r_{tot}}$ is the total number of rays, $n_e$ is the number of elements, $n_{rp}$ is the number of propagative rays, $n_r$ is the number of boundaries per element $e$, and $n_c$ is the number of corners per subdomain.

<table><caption>**Table 2**

Characteristics of workstation used for numerical examples described in Section 4.1.</caption><tr><td>Cores</td><td>16</td><td></td></tr><tr><td>Clock frequency</td><td>2.4</td><td>GHz</td></tr><tr><td>RAM</td><td>50</td><td>Gb</td></tr></table>

The subsequent analysis aims to study the influence of the number of rays and the number of elements on the convergence rate of the VTCR using the p-iteration and the p-h-iteration. The first one increases the number of rays up to 100 without dividing the domain as illustrated in Fig. 7. The p-h-iteration instead is a two-steps procedure. First, the domain is divided in *n* subdomains. Four configurations are chosen *n* = [1; 2; 4; 8]. Second, a p-iteration in each set-up is performed. Results are illustrated in Fig. 8.

It is interesting to note that the needed number of rays per subdomain to reach convergence decreases as the number of elements grows. Qualitatively this behavior can be explained observing that decreasing the characteristic dimension of a subdomain is equivalent to decreasing the frequency because there are less wave lengths per subdomain. In fact, if the domain becomes smaller, the frequency needed to present a medium-frequency behavior increases. Conversely, the total number of needed rays increases with the number of subdivisions. Therefore the effect of decreasing the characteristic dimension of each subdomain is not able to counteract the term $n_e$ in Eq. (57). This implies that there is a serious advantage in convergence rate if the number of subdomains is kept as low as possible. This is consistent with the study [25] developed in acoustic domain.

<figure><img src="image_14.png" /><figcaption>**Fig. 9.** Performance study of the numerical example described in Section 4.1.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 10.** Boundary conditions and geometric dimensions of the complex frame structure described in Section 4.2.</figcaption></figure>
**Performance analysis** Once the p-iteration proved to be the most advantageous choice, the focus passes on VTCR performances. In this Section the previous academic case is investigated over a frequency band and VTCR results are compared with a FEM reference calculated by ABAQUS®. These tests are performed on a workstation. Table 2 reports its principal characteristics.

According to [4,26] at low-frequency the error is dominated by the local error. The “a priori” rule-of-thumb to choose FEM element size $h_0$ is

$$h_0 k = \frac{2\pi}{10}, \qquad (58)$$

<table><caption>**Table 3**

Quantities of interest of the example in Section 4.2.</caption><tr><td>h</td><td>Thickness</td><td>1</td><td>mm</td></tr><tr><td>R</td><td>Radius of curvature</td><td>2.5</td><td>m</td></tr><tr><td>f</td><td>Frequency</td><td>500</td><td>Hz</td></tr><tr><td>E</td><td>Young modulus</td><td>200</td><td>GPa</td></tr><tr><td>ν</td><td>Poisson's ratio</td><td>0.33</td><td></td></tr><tr><td>ρ</td><td>Density</td><td>7850</td><td>kg/m³</td></tr><tr><td>η</td><td>Damping factor</td><td>0.01</td><td></td></tr></table>

<figure><img src="image_16.png" /><figcaption>(a) ABAQUS® solution.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 11.** Magnitude of the displacement of the vibrational problem of Section 4.2.</figcaption></figure>

where $k$ is the dimensionless wavenumber. Conversely, at mid-frequency the error is led by the pollution effect. In this case the rule-of-thumb suggested by Deraemaeker and Ihlenburg in [4,26] respectively is

$$h_0^2 k^3 = \frac{4\pi^2}{100L} \qquad (59)$$

where $L$ is a characteristic dimension of the numerical example considered. In this Section frequency varies over a frequency band [2000; 8000] Hz. In order to measure performances degrees of freedom (DoFs) cannot be directly compared since methods are intrinsically different. Nevertheless, for the sake of completeness, they are reported with time and memory consumption in Fig. 9.

At mid-frequency FEM suffers of pollution effect as proved in [4,26]. DoFs, time, and memory required increase faster as frequency grows. Conversely, VTCR is unaffected by pollution effect as illustrated by Fig. 9. Since VTCR remains stable, computational cost differences grow with frequency. At 2000 Hz consumption differences are already of one order of magnitude and increase with frequency. Due to low memory requirements VTCR could have been run on a laptop while FEM needs a workstation. Moreover, the difference in computational time is so massive that VTCR on a laptop would still have outperformed FEM on a workstation. The fact that these results are obtained comparing an unoptimized MATLAB® program against a FEM implemented in the commercial code ABAQUS® confirms even further VTCR effectiveness.

## 4.2. Second numerical example: complex frame structure

The present Section studies a complex frame structure. A triple-joint structure is investigated and VTCR solution compared with a FEM reference.

Geometric dimensions and boundary conditions are reported in Fig. 10 while the other important quantities are stated in Table 3.

The test is performed at low-frequency since FEM computational costs are prohibitive at mid-frequency. In this frequency range the rule-of-thumb reported in Eq. (58) is valid.

The error is based on the kinetic energy $E_K$:

$$err = \frac{E_K(\mathbf{u}_{ABAQUS} - \mathbf{u}_{VTCR})}{E_K(\mathbf{u}_{ABAQUS})} \qquad (60)$$

where $\mathbf{u}_{ABAQUS}$ is the displacement calculated by ABAQUS® and $\mathbf{u}_{VTCR}$ is the displacement calculated by the VTCR theory. The total number of rays is

$$n_{r_{tot}} = n_e \left( 3n_{r_p} + n_I \frac{n_{r_p}}{2} + n_I + n_c \right). \qquad (61)$$

<figure><img src="image_18.png" /><figcaption>**Fig. 12.** Error of the complex frame structure example presented in Section 4.2. The definition of the error is given in Eq. (60). Convergence is achieved around 80 propagative rays per subdomain (1224 total DoFs). The remaining error (around 10%) stays stable afterwards and it is due small differences in theories. The reference solution is computed by ABAQUS® using 372084 DoFs.</figcaption></figure>
**Table 4**

<table><caption>DoFs, time, and memory consumption comparisons between VTCR and FEM of the numerical example described in Section 4.2.</caption><thead><tr><th></th><th>DoFs</th><th>Wallclock time (s)</th><th>RAM (Mb)</th></tr></thead><tbody><tr><td>VTCR</td><td>1224</td><td>91</td><td>85</td></tr><tr><td>FEM</td><td>372084</td><td>24</td><td>202</td></tr></tbody></table>

In this case the in-plane component is no more negligible therefore two new types of propagative rays must be considered. For the sake of simplicity the number of in-plane propagative rays is kept equal to the number of out-of-plane waves.

The displacement magnitude of the ABAQUS® and VTCR solutions are compared in Fig. 11. The general behavior is very similar. This denotes that the VTCR correctly addresses the three-shell-joint. The error with the number of propagative rays is reported in Fig. 12. Error decreases until 80 propagative rays. After that, a there is a plateau at ≈10% due to small differences in theories. For example, since FEM approximates the equilibrium equations while VTCR does not, the wavelengths calculated are slightly different. This leads to a small discrepancy between solutions. Table 4 compares DoFs and performances between VTCR and FEM. As was expected, at this frequency FEM outperforms VTCR in terms of computational time because the pollution effect is negligible. Nevertheless VTCR still outrun FEM in terms of memory consumption.

## 5. Conclusions

In order to increase its flexibility and effectiveness, VTCR was extended to shallow shell theory. The uniqueness property of the shallow-shell-VTCR solution was provided and two numerical examples were presented to validate the theory. The first one was an academic example where the analytic solution is known. Some convergence tests were run to analyze p- and h-refinements. p-Refinement proved to be the most effective in terms of total DoFs required as well as memory and time consumptions. After that, VTCR performs were compared with a FEM reference over a frequency band. At mid-frequency VTCR outperforms FEM of some order of magnitude in terms of time and memory consumption. This happens because FEM is affected by pollution effect while VTCR is not. The second example focused on coupling condition testing. Due to FEM computational costs, this test was performed at low-frequency. In this case, as expected, FEM outperformed VTCR in terms of computational time. Nevertheless VTCR consumed less memory. For this reason VTCR is a more suited method when mid-frequency is involved or when memory consumption is the principal concern.

In the near future, in order to extend the theory, orthotropic materials will be investigated. Moreover, VTCR will be coupled with a Reduced Order Model to extend its applicability to a frequency band.

## Acknowledgements

The authors gratefully acknowledge the “Centre National d'Études Spatiales (CNES)” and “Airbus Defense and Space”

## References

[1] Lyon RH. Statistical energy analysis of dynamical systems: theory and applications. Cambridge, MA: MIT press; 1975.

[2] Hughes TJR. The finite element method: linear static and dynamic finite element analysis. Courier Dover Publications; 2012.

[3] Hall WS. Boundary element method. Springer; 1994.

[4] Deraemaeke A. Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions. Int J Numer Methods Eng 1999;46(4):471-99. http://dx.doi.org/10.1002/(SICI)1097-0207(19991010)46:4<471::AID-NME684>3.0.CO;2-6.

[5] Mace B. Statistical energy analysis, energy distribution models and system modes. J Sound Vib 2003;264(2):391-409. http://dx.doi.org/10.1016/S0022-460X(02)01201-4.

[6] Farhat C, Roux FX. A method of finite element tearing and interconnecting and its parallel solution algorithm. Int J Numer Methods Eng 1991;32(6):1205-27.

[7] Liu Y. Fast multipole boundary element method: theory and applications in engineering. Cambridge University Press; 2009.

[8] Tezaur R, Kalashnikova I, Farhat C. The discontinuous enrichment method for medium-frequency Helmholtz problems with a spatially variable wavenumber. Comput Methods Appl Mech Eng 2014;268:126-40. http://dx.doi.org/10.1016/j.cma.2013.08.017.

[9] Strouboulis T, Hidajat R. Partition of unity method for Helmholtz equation: q-convergence for plane-wave and wave-band local bases. Appl Math 2006;51(2):181-204.

[10] Cessenat O, Despres B. Application of an ultra weak variational formulation of elliptic PDEs to the two-dimensional Helmholtz problem. SIAM J Numer Anal 1998;35(1):255-99.

[11] De Rosa S, Franco F. On the use of the asymptotic scaled modal analysis for time-harmonic structural analysis and for the prediction of coupling loss factors for similar systems. Mech Syst Signal Process 2010;24(2):455-80. http://dx.doi.org/10.1016/j.vmssp.2009.07.008.

[12] Soize C. Reduced models in the medium frequency range for general external structural-acoustics systems. J Acoust Soc Am 1998;103(6):3393-406. http://dx.doi.org/10.1121/1.423052.

[13] Bouillard P, Suleaub S, Suleau S. Element-free Galerkin solutions for Helmholtz problems: formulation and numerical assessment of the pollution effect. Comput Methods Appl Mech Eng 1998;162(14):317-35.

[14] Perrey-Debain E, Trevelyan J, Bettess P. Wave boundary elements: a theoretical overview presenting applications in scattering of short waves. Eng Anal Boundary Elements 2004;28:131-41.

[15] Desmet W, Van Hal B, Sas P, Vandeputte D, Hal BV. A computationally efficient prediction technique for the steady-state dynamic analysis of coupled vibro-acoustic systems. Adv Eng Softw 2002;33(7):527-40. http://dx.doi.org/10.1016/S0965-9978(02)00062-5.

[16] Genechten BV, Vandeputte D, Desmet W. A direct hybrid finite element-wave based modelling technique for efficient coupled vibro-acoustic analysis. Comput Methods Appl Mech Eng 2011;200(5):742-64. http://dx.doi.org/10.1016/j.cma.2010.09.017.

[17] 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. http://dx.doi.org/10.1108/02644400110365879.

[18] Riou H, Ladevèze P, Rouch P. Extension of the variational theory of complex rays to shells for medium-frequency vibrations. J Sound Vib 2004;272(1):341-60. http://dx.doi.org/10.1016/S0022-460X(03)00775-2.

[19] Chevreuil M, Ladevèze P, Rouch P. Transient analysis including the low- and the medium-frequency ranges of engineering structures. Comput Struct 2007;85(17-18):1431-44.

[20] Kovalevsky L, Ladevèze P, Riou H, Bonnet M. The variational theory of complex rays for three-dimensional helmholtz problems. J Comput Acoust 2012;1-25.

[21] Ladevèze P, Riou H. Calculation of medium-frequency vibrations over a wide frequency range. Comput Methods Appl Mech Eng 2005;194(27-29):3167-91. http://dx.doi.org/10.1016/j.cma.2004.08.009.

[22] Barbarulo A, Ladevèze P, Riou H, Kovalevsky L. Proper Generalized Decomposition applied to linear acoustic: a new tool for broad band calculation. J Sound Vib 2014;333(11):2422-31. http://dx.doi.org/10.1016/j.jsv.2014.01.014.

[23] Ventsel E, Krauthammer T. Thin plates and shells: theory: analysis, and applications. Basel: CRC Press; 2001.

[24] van der Heijden AMA. WT Koiter's elastic stability of solids and structures. Cambridge: Cambridge University Press; 2009.

[25] Riou H, Ladevèze P, Sourcis B. The multiscale VTCR approach applied to acoustics problems. J Comput Acoust 2008;16(4):487-505.

[26] Ihlenburg F. Finite element analysis of acoustic scattering. Springer; 1998.