# Extension of the variational theory of complex rays to orthotropic shallow shell structures

Alessandro Cattabiani, Andrea Barbarulo, Hervé Riou, Pierre Ladevèze

► To cite this version:

Alessandro Cattabiani, Andrea Barbarulo, Hervé Riou, Pierre Ladevèze. Extension of the variational theory of complex rays to orthotropic shallow shell structures. Advances in Aircraft and Spacecraft Science, 2016, 3 (3), pp.317-330. 10.12989/aas.2016.3.3.317. hal-01647841

HAL Id: hal-01647841

https://hal.science/hal-01647841

Submitted on 23 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
Extension of the variational theory of complex rays
to orthotropic shallow shell structures

Alessandro Cattabiani*¹, Andrea Barbarulo²ª, Hervé Riou¹ᵇ and Pierre Ladevèze¹ᶜ

¹LMT-Cachan, ENS Cachan/CNRS/Paris-Saclay-University, 61 avenue du Président Wilson,
F-94230 Cachan France

²Laboratoire MSSMat, Centrale-Supélec/CNRS UMR8579, Grande Voie des Vignes,
Cedex, Chatenay-Malabry, 92290, France

**Abstract.** Nowadays, the interest of aerospace and automotive industries on virtual testing of medium-frequency vibrational behavior of shallow shell structures is growing. The development of software capable of predicting the vibrational response in such frequency range is still an open question because classical methods (i.e., FEM, SEA) are not fully suitable for the medium-frequency bandwidth. In this context the Variational Theory of Complex Rays (VTCR) is taking place as an ad-hoc technique to address medium-frequency problems. It is a Trefftz method based on a weak variational formulation. It allows great flexibility because any shape function that satisfies the governing equations can be used. This work further develops such theory. In particular, orthotropic materials are introduced in the VTCR formulation for shallow shell structures. A significant numerical example is proposed to show the strategy.

**Keywords:** VTCR; medium-frequency; orthotropic materials; shallow shells

**1. Introduction**

Recently, aerospace and automotive industries have shown interest towards efficient virtual testing of the vibration response. Shallow shell structures are widely used in air- and spacecraft due to their high resistance and lightweight. Aircraft wings, fuselages, tails and spacecraft launchers are all composite shallow shell structures. The equilibrium Eqs. are quite complex and in almost every real case an analytic solution cannot be obtained. Therefore, an effective method to predict vibrational behavior in shallow shell structures is needed.

The Modal Overlap Factor Lyon (1975) defines three zones: low, mid and high frequency range. The low-frequency range has been extensively studied by deterministic techniques such as the Finite Element Method (FEM) Hughes (2012) and the Boundary Element Method (BEM) Hall (1994). On the other side, the high-frequency range can be addressed by averaging methods as the Statistical Energy Analysis (SEA) Lyon (1975). This technique neglects almost entirely spatial

<footnote>*Corresponding author, Ph.D. Student, E-mail: alessandro.cattabiani@gmail.com</footnote>

<footnote>ªProfessor, E-mail: andrea.barbarulo@ecp.fr</footnote>

<footnote>ᵇProfessor, E-mail: riou@lmt.ens-cachan.fr</footnote>

<footnote>ᶜProfessor Emeritus, E-mail: ladeveze@lmt.ens-cachan.fr</footnote>
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 suited in this frequency domain since the phenomena variation length is very small if compared to characteristic dimensions of the structure. For this reason computational costs explode Deraemaeker (1999). On the other hand the classical SEA is not indicated because the key assumptions of the theory might be unsatisfied Mace (2003). However a lot of work has been done to extend such theories to the medium-frequency range Farhat and Roux (1991), Liu (2009), Tezaur *et al.* (2014), De Rosa and Franco (2010), Soize (1998).

There is also a class of methods developed specifically for the medium-frequency range known as Trefftz methods such as the partition of unity method Strouboulist and hidajat (2006), the ultra-weak variational method Cessenatu and Despres (1998), the least-squares method Monk and Wang (1999), the discontinuous enrichment method Farhat *et al.* (2001), the element-free Galerkin method Bouillard *et al.* (1998), the wave boundary element method Perrey-Debain *et al.* (2004) or the wave-based method Desmet *et al.* (2002), Genechten *et al.* (2011). 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 Eqs. 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 providing great flexibility. It has already been applied to plate theory Ladevèze *et al.* (2001), shallow shell theory Cattabiani *et al.* (2015), to general shell theory Riou *et al.* (2004), to transient dynamics Chevreuil *et al.* (2007), to orthotropic plates Kovalevsky *et al.* (2014), to 3D acoustic Kovalevsky *et al.* (2012) and, on a wide frequency band Ladeveze and Riou (2005), Barbarulo *et al.* (2014).

In the present work the VTCR is applied to orthotropic shallow shell structures. Section 2 illustrates the general shallow shell theory and the notation. Section 3 presents the VTCR adapted for the shallow shell theory and its improvements. In the end, Section 4 validates the method on a relevant numerical example comparing the VTCR solution with a FEM reference and analyzing performances.

## 2. General shallow shell theory

In this Section the equilibrium and boundary equations of the general kirchhoff-Love theory for orthotropic shallow shells are examined. The present development is akin to the one provided in Ventsel and Krauthammer (2001) and van der Heijden (2009).

### 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_{\perp i}$ symbol refers to a generic boundary of $i$ where condition $\blacksquare^1$ 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_{\perp i}$ is used.

<footnote n="1">¹The symbol $\blacksquare$ stands for a generic constraint i.e. $\blacksquare = u$ is the displacement field.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1 Generic frame structure described in Section 2.1</figcaption></figure>

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

The generic corner shared among subdomains is indicated with C. The over-line symbol $\bar{\}$ indicates that quantity $\hat{\}$ 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. 2(a)) and on the corners (Fig. 2(b)) as well as a distributed load per unit surface $\bar{g}_i$ (Fig. 2(d)). The displacement constraint $\bar{u}_i = [\bar{v}_i', \bar{w}_i]'$ along $\partial\bar{u}_i\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}_i\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,\hat{n}_i}$ is imposed along $\partial\bar{w}_{i,\hat{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. 2(c)).

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.

<footnote n="1">¹$\hat{\}$ is the transpose operator.</footnote>
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 \boldsymbol{\phi}_i \\
\boldsymbol{\phi}_i &= \nabla w_i - \mathbf{R}_i \cdot \mathbf{v}_i \\
\mathbf{R}_i &= \begin{bmatrix} 1/R_{xi} & 0 \\ 0 & 1/R_{yi} \end{bmatrix}
\end{aligned}
$$

where $\mathbf{u}_i^z$ is the displacement through 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. $\mathbf{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 \text{ finite energy displacement set,} \quad (1) $$

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

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

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

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

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

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

$$ \mathbf{N}_i = \mathbf{B}_i : (\mathbf{E}_i - \mathbf{R}_i w_i), \quad (8) $$

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

$$ \mathbf{B}_i = \begin{bmatrix} B_{xi} & \nu_{xyi} B_{yi} & 0 \\ \nu_{yxi} B_{xi} & B_{yi} & 0 \\ 0 & 0 & B_{Gi} \end{bmatrix}, \quad (10) $$

$$ B_{xi} = \frac{E_{xi}h_{i}}{1 - \nu_{xyi}\nu_{yxi}}, \quad B_{yi} = \frac{E_{yi}h_{i}}{1 - \nu_{xyi}\nu_{yxi}}, \quad B_{Gi} = G_{xyi}h_{i} \quad (11) $$

$$ \mathbf{D}_i = \begin{bmatrix} D_{xi} & \nu_{xyi} D_{yi} & 0 \\ \nu_{yxi} D_{xi} & D_{yi} & 0 \\ 0 & 0 & D_{Gi} \end{bmatrix}, \quad (12) $$

$$ D_{xi} = \frac{E_{xi}h_i^3}{12(1 - \nu_{xyi}\nu_{yxi})}, \quad D_{yi} = \frac{E_{yi}h_i^3}{12(1 - \nu_{xyi}\nu_{yxi})}, \quad D_{Gi} = \frac{G_{xyi}h_i}{6} \quad (13) $$

$$ E_{xi} = E_{x0i}(1 + j\eta_{xi}), \quad (14) $$

$$ E_{yi} = E_{y0i}(1 + j\eta_{yi}), \quad (15) $$

$$ G_i = G_{0i}(1 + j\eta_{Gi}), \quad (16) $$

where $\mathbf{B}_i$ and $\mathbf{D}_i$ are Hooke's plane stress operators concerning in-plane and out-of-plane stresses respectively, $\rho_i$ is the density, $h_i$ is the shell thickness, $E_{x0i}$ and $E_{y0i}$ are the Young moduli along directions $x$ and $y$ respectively, $\eta_{xi}$ and $\eta_{yi}$ are the relative damping coefficients of the Young
moduli, $v_{xyi}$ and $v_{yxi}$ are the Poisson's ratios $(x, y)$ and $(y, x)$ respectively, $G_{xyi}$ is the in-plane shear modulus, $\eta_{Gi}$ is its specific damping coefficient, ■: ■ is the inner matrix product operator, $[\cdot]_{sym} = \frac{1}{2}(\cdot + \cdot')$ is the symmetric part operator, $e_{i(x,y)}$ is the component $(x, y)$ of the matrix $E_i$, and $N_i$ and $M_i$ are the stress and stress moment resultants tensors respectively. As suggested in Cattabiani et al. (2015), in-plane inertia is not neglected. The sub-space of $D_i$ associated with homogenized conditions ($\bar{g}_i = 0$) is denoted as $D_{0i} = \{\delta u_i, \delta N_i, \delta M_i\}$. This definition will be useful in the next Sections.

## 2.2 Boundary conditions

In order to present a well-posed problem, three conditions must be imposed along each boundary and one on each corner. 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,\hat{n}_i}$ or $\bar{m}_i$),

4. an out-of-plane condition, on corners either a displacement constraint or a point load ($\bar{w}_{Ci}$ or $\bar{q}_{Ci}$). Boundary conditions in the most general case are

$$ v_i = \begin{cases} \bar{v}_i \quad \exists \quad \bar{v}_i \\ v_{j\Gamma} \quad \notin \quad \bar{v}_i \end{cases}, \qquad (17) $$

$$ w_i = \begin{cases} \bar{w}_i \quad \exists \quad \bar{w}_i \\ w_{j\Gamma} \quad \notin \quad \bar{w}_i \end{cases}, \qquad (18) $$

$$ w_{i,\hat{n}_i} = (\nabla w_i)\hat{n}_i = \begin{cases} \bar{w}_{i,\hat{n}_i} \quad \exists \quad \bar{w}_{i,\hat{n}_i} \\ \bar{w}_{j\Gamma,\hat{n}_j} \quad \notin \quad \bar{w}_{i,\hat{n}_i} \end{cases}, \qquad (19) $$

$$ b_i = N_i \hat{n}_i = \bar{b}_i - \sum_{j\Gamma=1}^{n_\Gamma} b_{j\Gamma}, \qquad (20) $$

$$ q_i = (\nabla \cdot \mathbf{M}_i)\hat{n}_i + \nabla(\hat{t}'_i \mathbf{M}_i \hat{n}_i) \cdot \hat{t}_i =, $$

$$ = (\nabla \cdot \mathbf{M}_i)\hat{n}_i + (\hat{t}'_i \mathbf{M}_i \hat{n}_i)_{,\hat{t}_i} = \bar{q}_i - \sum_{j\Gamma=1}^{n_\Gamma} q_{j\Gamma}, \qquad (21) $$

$$ m_i = \hat{n}'_i \mathbf{M}_i \hat{n}_i = \bar{m}_i - \sum_{j\Gamma=1}^{n_\Gamma} m_{j\Gamma} \qquad (22) $$

<figure><img src="image_4.png" /><figcaption>Fig. 3 $n_\Gamma + 1$ subdomains ($\Omega_i$ included) sharing the same boundary $\Gamma$</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4 $n_c + 1$ subdomains ($\Omega_i$ included) sharing the same corner C</figcaption></figure>

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

Corner conditions in the most general case are

$$ w_{Ci} = \begin{cases} \bar{w}_{Ci} & \exists \bar{w}_{Ci} \\ w_{Cj\Gamma} & \nexists \bar{w}_{Ci} \end{cases}, \qquad (23) $$

$$ 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 (24) $$

where $\hat{\mathbf{n}}_{1i}$ and $\hat{\mathbf{n}}_{2i}$ are outward normal unit vectors of the two boundaries of $\Omega_i$ sharing the corner C. $\hat{\mathbf{t}}_{1i}$ and $\hat{\mathbf{t}}_{2i}$ are their respective tangent unit vectors directed towards the corner. As for $j_\Gamma$, 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.

### 3. Shallow shell-VTCR theory

In this Section the VTCR developed for shallow shell structures is illustrated. Since the VTCR is a Trefftz method, the solution is searched in a function set that satisfy equilibrium equations. Boundary and corner residuals are addressed in weak form. The weak variational problem is:

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

$$ \sum_{i=1}^{n} \left\{ \begin{aligned} & \int_{\partial \bar{v}_i \Omega_i} \delta \mathbf{b}_i^H (\mathbf{v}_{si} - \bar{\mathbf{v}}_i) \, ds - \int_{\partial \bar{b}_i \Omega_i} \delta \mathbf{v}_i^H (\mathbf{b}_{si} - \bar{\mathbf{b}}_i) \, ds \\ & + \int_{\partial \bar{w}_i \Omega_i} \delta q_i^H (w_{si} - \bar{w}_i) \, ds - \int_{\partial \bar{q}_i \Omega_i} \delta w_i^H (q_{si} - \bar{q}_i) \, ds \\ & - \int_{\partial \bar{w}_{i, \hat{n}_i} \Omega_i} \delta m_i^H (w_{si, \hat{n}_i} - \bar{w}_{i, \hat{n}_i}) \, ds + \int_{\partial \bar{m}_i \Omega_i} \delta w_{i, \hat{n}_i}^H (m_{si} - \bar{m}_i) \, ds \end{aligned} \right. $$
$$
\begin{align}
&+\delta w_{Ci}^{H}(q_{Csi} - \bar{q}_{Ci})|_{\partial\partial\bar{q}_{Ci}\Omega_i} - \delta q_{Ci}^{H}(w_{Csi} - \bar{w}_{Ci})|_{\partial\partial\bar{w}_{Ci}\Omega_i} \nonumber \\
&+\frac{n_{\Gamma}}{n_{\Gamma}+1} \int_{\Gamma} \delta \mathbf{b}_{i}^{H} \mathbf{v}_{si} + \delta q_{i}^{H} w_{si} - \delta m_{i}^{H} w_{si, \hat{n}_{i}} \, ds - \frac{n_{C}}{n_{C}+1} \delta q_{Ci}^{H} w_{Csi}|_{C} \nonumber \\
&-\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{n}_{i}}^{H} m_{si} \, ds + \frac{1}{n_{C}+1} \delta w_{Ci}^{H} q_{Csi}|_{C} \nonumber \\
&-\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{n}_{j}} + \delta \mathbf{v}_{i}^{H} \mathbf{b}_{sj} + \delta w_{i}^{H} q_{sj} - \delta w_{i, \hat{n}_{i}}^{H} m_{sj} \, ds \right. \nonumber \\
&\qquad \left. -\frac{1}{n_{\Gamma}+1} (\delta q_{Ci}^{H} w_{Csi} + w_{Ci}^{H} q_{Csi})|_{C} \right\} = 0 \tag{25} \\
&\forall \{\mathbf{u}_i, \mathbf{N}_i, \mathbf{M}_i\} \in \mathcal{D}_{0i} \nonumber
\end{align}
$$

where $\mathcal{D}_{0i}$ is the test functions space being VTCR a Galerkin method. Ladevèze in Ladevèze et al. (2003) proved uniqueness and existence properties in general elastic theory. Since shell theory is a particularization of such theory, those demonstrations (properly adapted to meet shell theory approximations) hold in our specific case.

## 3.1 Discretization

Different sets of shape functions can be chosen. In Riou et al. (2004) plane waves were chosen, meanwhile in Kovalevsky (2012) Fourier functions were selected. In this work plane waves are used since integrals over straight boundaries can be calculated analytically, drastically reducing computational costs. Their generic form is

$$f_{SF} = a_i \hat{c}_i e_i^{k_i \cdot x_{rel}} \qquad (26)$$

where $k_i = k_i \hat{k}_i$ is the wave vector, $\hat{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 $x_{rel}$ is the relative position vector. The relative coordinate system is located in the geometric center<fnref n="1" /> of the element and $\hat{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{k}_i$ is set a priori. This is the only discretized parameter of the VTCR. This class of shape functions can be divided into propagative and evanescent waves. These are both needed as explained in Riou et al. (2004). The difference between them lies in the unit direction vector $\hat{k}_i$ that, in case of propagative waves, is

$$\hat{k}_i = T \cdot \begin{bmatrix} 1 \\ 0 \end{bmatrix} \qquad (27)$$

<footnote n="1">¹The geometric center is chosen just for the sake of convenience.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 5 Qualitative behavior of a propagative wave described in Section 3.1</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6 Qualitative behavior of an evanescent wave described in Section 3.1</figcaption></figure>

where

$$ T = \begin{bmatrix} \cos(\theta_i) & -\sin(\theta_i) \\ \sin(\theta_i) & \cos(\theta_i) \end{bmatrix} \qquad (28) $$

is the rotation matrix and $\theta_i \in [0; 2\pi)$ is the discretized angle of the plane wave direction. Their qualitative behavior is pictured in Fig. 5.

The case of evanescent waves is slightly different. By definition these waves are evanescent along one direction and propagative along the orthogonal one as depicted in Fig. 6.

In Cattabiani et al. (2015) the unit direction vector of these rays is

$$ \hat{k}_i = T \cdot \begin{bmatrix} \cosh(\phi_i) \\ \sinh(\phi_i) \end{bmatrix} \qquad (29) $$

where $\phi_i$ is a different discretized angle in the complex domain. Substituting the generic shape function defined in Eq. (26) with a given unit direction vector $\hat{k}_i$ in Eqs. (3) and (4) leads to a linear set of equilibrium equations that in matrix form is

$$ Z_i c_i a_i = 0 \qquad (30) $$

where $Z_i = Z_i(k_i)$ is a 3 × 3 matrix that depends on $k_i$. In order to get untrivial results ($a_i = 0$) the dispersion equation must be imposed

$$ \det[Z_i] = 0. \qquad (31) $$

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

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

As explained in Cattabiani et al. (2015), since in-plane inertia is not neglected, Eq. (31) provides four wavenumbers for each wavevector. Three of them are propagative rays and one corresponds to an evanescent wave. Moreover, one of the three propagative rays controls the out-of-plane behavior while the other two carries in-plane stresses and displacements.

## 3.2 Orthotropic materials

Eqs. (27) and (29) are defined for isotropic materials in a cartesian coordinate system. This Section
illustrates necessary corrections for orthotropic materials. Kovalevsky et al. (2014) suggested to add a correction matrix $\mathbf{O}_i$ to address orthotropic materials in plates

$$ \mathbf{k}_i = \mathbf{O}_i \hat{\mathbf{k}}_i k_i, \qquad (33) $$

where

$$ \mathbf{O}_i = \sqrt[4]{\rho_i h_i \omega^2} \begin{bmatrix} \sqrt[4]{D_{xi}} & 0 \\ 0 & -\sqrt[4]{D_{yi}} \end{bmatrix}. \qquad (34) $$

In this formulation $\mathbf{O}_i = \mathbf{O}_i(\rho_i h_i \omega^2)$. This dependency can lead to numerical difficulties during calculation of $k_i$ and $\hat{\mathbf{c}}_i$ because $\mathbf{O}_i \hat{\mathbf{k}}_i$ is no more a unit vector. For example, let us consider the case where $\mathbf{k}_i = [a, 0]'$. If the direction vector corrected for orthotropic materials injected in the equilibrium equations is $\mathbf{O}_i \hat{\mathbf{k}}_i = [a + b, 0]'$ where $b \gg a$, the wavenumber provided by these equations should be $k_i = a/(b + a)$. Since this value is computed numerically, the round-off error could approximate the wavenumber as $k_i \approx 1/b$ degrading results. Hence, the direction vector should be kept unitary to avoid numerical difficulties. In the present approach we use a slight different formulation that solves such problem

$$ \mathbf{O}_i = \sqrt[8]{D_{xi} D_{yi}} \begin{bmatrix} \sqrt[4]{D_{xi}} & 0 \\ 0 & -\sqrt[4]{D_{yi}} \end{bmatrix}. \qquad (35) $$

In this case the wavevector is dimensionless suppressing numerical difficulties that can arise when high frequencies are involved. In other words, this modification increases the computational precision of $k_i$ and $\hat{\mathbf{c}}_i$.

## 4. Numerical example: wing

This Section investigates the vibrational response of a composite aircraft wing frame structure. The VTCR solution is compared with a FEM reference, investigating accuracy and performances.

<figure><img src="image_8.png" /><figcaption>Fig. 7 Boundary conditions and geometric dimensions of a typical composite wing described in Section 4.1</figcaption></figure>
<table><caption>Table 1 Studied frequencies, thicknesses, and material properties of Section 4.1</caption><thead><tr><th>f</th><th>Frequency</th><th>500 Hz</th></tr></thead><tbody><tr><td>h<sub>1</sub>, ..., h<sub>12</sub></td><td>Thicknesses</td><td>3 mm</td></tr><tr><td>E<sub>x1</sub>, ..., E<sub>x12</sub></td><td>Young moduli</td><td>125 GPa</td></tr><tr><td>E<sub>y1</sub>, ..., E<sub>y12</sub></td><td>Young moduli</td><td>60 GPa</td></tr><tr><td>G<sub>xy1</sub>, ..., G<sub>xy12</sub></td><td>Shear moduli</td><td>18 GPa</td></tr><tr><td>ν<sub>xy1</sub>, ..., ν<sub>xy12</sub></td><td>Poisson's ratios</td><td>0.3</td></tr><tr><td>ρ<sub>1</sub>, ..., ρ<sub>12</sub></td><td>Densities</td><td>2000 kg/m<sup>3</sup></td></tr><tr><td>h<sub>13</sub>, h<sub>14</sub></td><td>Thicknesses</td><td>1 cm</td></tr><tr><td>E<sub>13</sub>, E<sub>14</sub></td><td>Young moduli</td><td>72 GPa</td></tr><tr><td>ν<sub>13</sub>, ν<sub>14</sub></td><td>Poisson's ratios</td><td>0.33</td></tr><tr><td>ρ<sub>13</sub>, ρ<sub>14</sub></td><td>Densities</td><td>2800 kg/m<sup>3</sup></td></tr></tbody></table>

## 4.1 Description of the vibrational problem

Geometry and boundary conditions are illustrated in Fig. 7. It is divided in fourteen subdomains $\Omega_i$. Spars [$\Omega_{13}$, $\Omega_{14}$ are isotropic aluminum plates. The upper and lower skins as well as the leading edge are orthotropic composite shallow shells. For the sake of simplicity, every structural damping factor is $\eta = 0.01$. The left wing cross-section is clamped. The remaining boundaries are free but the four dotted boundaries where is applied a distributed and oscillatory load $\mathbf{p} = [0,0,1]^{\prime}$N/m. Thicknesses and material properties are reported in Table 1.

## 4.2 VTCR and FEM discretizations and error indicator

This Subsection describes how VTCR and FEM solutions are discretized to ensure convergence. The common FEM "rule-of-thumb" prescribes to keep the mesh characteristic dimension smaller than one out of ten times the wavelength. Ihlenburg in acoustic Ihlenburg (1998) and Deraemaeker for general Helmholtz problems Deraemaeker (1999) proved that this rule is no more valid at mid-frequency due to high-scattering behavior. In particular, they affirmed that pollution error becomes predominant as the wavenumber increases. They suggested a corrected version

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

where $L$ is a characteristic dimension (0.5 m in this numerical example) of the considered problem, $h_0$ is the characteristic mesh element length, and $k$ is the greatest possible wavenumber.

In VTCR the complex frame structure is divided in fourteen subdomains. The ray number increases until convergence. Displacements unequivocally define the solution. Since at this frequency a slight difference in theories can lead to very different frequency responses, the error indicator is based on the total kinetic energy

$$E_K(\mathbf{u}) = \frac{1}{2} \rho h \omega^2 \int_S \mathbf{u}^H \cdot \mathbf{u} \, ds, \qquad (37)$$

$$err = \frac{E_K(\mathbf{u}_{ABAQUS}) - E_K(\mathbf{u}_{VTCR})}{E_K(\mathbf{u}_{ABAQUS})}, \qquad (38)$$
<table><caption>Table 2 Characteristics of the workstation used for the numerical example described in 4.1</caption><tr><td>Cores</td><td>16</td></tr><tr><td>Clock frequency</td><td>2.4 GHz</td></tr><tr><td>RAM</td><td>50 Gb</td></tr></table>

<figure><img src="image_9.png" /><figcaption>Fig. 8 Displacement magnitude portraits comparison at 500 Hz of the wing described in Section 4.1</figcaption></figure>

<table><caption>Table 3 DoFs and performances comparison at 500 Hz of the numerical example described in Section 4.1</caption><thead><tr><td></td><td>DoFs</td><td>Wallclock Time h:m:s</td><td>RAM consumption</td></tr></thead><tbody><tr><td>VTCR</td><td>5936</td><td>00:00:47</td><td>470 Mb</td></tr><tr><td>FEM</td><td>18 520 944</td><td>02:36:21</td><td>9.676 Gb</td></tr></tbody></table>

where $u_{ABAQUS}$ is the displacement calculated by Abaqus and $u_{VTCR}$ is the displacement calculated by VTCR theory.

Table 2 reports principal characteristics of the used workstation.

### 4.3 VTCR-FEM comparison at 500 Hz

The FEM and VTCR displacement magnitudes portraits at 500 Hz are confronted in Fig. 8. DoFs as well as performances are compared in Table 3.
## 4.4 Remarks

The discrepancy between VTCR and FEM is $E_K(u) = 10\%$. Even if FEM and VTCR solutions do not perfectly match, the energy level is similar. The amplitude magnitude portrait difference is due to small theory differences that at mid-frequency can be determinant. As explained in Riou et al. (2004), at this frequency the relatively high leading edge curvature reflect vibrational waves. FEM and VTCR solutions confirm such behavior demonstrating that even the shallow shell-VTCR addresses this effect. Since VTCR is a Trefftz method, the number of required subdomains is very small (just fourteen in this case). Thus, the VTCR DoFs number is very small too. In order to study performances DoFs numbers cannot be compared directly since FEM and VTCR are intrinsically different. Nevertheless, Table 3 illustrates the great saving of memory and time consumption of the VTCR method even at 500 Hz.

# 5. Conclusions

The present work extended the VTCR to orthotropic shallow shells. Composite wings, fuselages, and rocket tanks are some examples of such structures in air- and spacecraft. Composite materials were addressed introducing a corrective term in the wave direction vector. Such correction is necessary to address different sound speeds along different principal directions. In this way, the vibrational waves (which are VTCR shape functions) propagate at different speeds along principal directions. In contrast to previous works, the corrective term is dimensionless avoiding numerical difficulties that can arise at high frequency.

A relevant complex numerical example of a composite aircraft wing was presented to illustrate the method. The VTCR solution calculated at 500 Hz was compared with a FEM reference and performances were studied. The DoFs comparison between VTCR and FEM cannot be directly used to investigate performances due to the different DoFs meanings. Nevertheless, VTCR greatly outperformed FEM in terms of time and memory consumption proving to be a very promising method.

## Acknowledgments

The research described in this paper was financially supported by the Centre National d'Études Spatiales (CNES) and Airbus Defense and Space.

## References

Barbarulo, A., Ladevèze, P., Riou, H. and Kovalevsky, L. (2014), “Proper generalized decomposition applied to linear acoustic: a new tool for broad band calculation”, *J. Sound Vib.*, **333**(11), 2422-2431.

Bouillard, P., Suleaub, S. and Suleau, S. (1998), “Element-free galerkin solutions for helmholtz problems: formulation and numerical assessment of the pollution effect”, *Comput. Meth. Appl. Mech. Eng.*, **162**(14), 317-335.

Cattabiani, A., Riou, H., Barbarulo, A., Ladevèze, P., Bézier, G. and Troclet, B. (2015), “The Variational Theory of Complex Rays applied to the shallow shell theory”, *Comput. Struct.*, **158**, 98-107.

Cessenat, O. and Despres, B. (1998), “Application of an ultra weak variational formulation of elliptic PDEs to the two-dimensional Helmholtz problem”, *SIAM J. Numer. Anal.*, **35**(1), 255-299.
Chevreuil, M., Ladevèze, P. and Rouch, P. (2007), "Transient analysis including the low- and the medium-frequency ranges of engineering structures", *Comput. Struct.*, **85**(17-18), 1431-1444.

De Rosa, S. and Franco, F. (2010), “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 Pr.*, **24**(2), 455-480.

Deraemaeker, A. (1999), “Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions”, *Int. J. Numer. Meth. Eng.*, **46**(4), 471-499.

Desmet, W., Van Hal, B., Sas, P., Vandepitte, D. and Hal, B.V. (2002), “A computationally efficient prediction technique for the steady-state dynamic analysis of coupled vibro-acoustic systems”, *Adv. Eng. Softw.*, **33**(7), 527-540.

Farhat, C. and Roux, F.X. (1991), “A method of finite element tearing and interconnecting and its parallel solution algorithm”, *Int. J. Numer. Meth. Eng.*, **32**(6), 1205-1227.

Farhat, C., Harari, I. and Franca, L.P. (2001), “The discontinuous enrichment method”, *Comput. Meth. Appl. Mech. Eng.*, **190**, 6455-6479.

Genechten, B.V., Vandepitte, D. and Desmet, W. (2011), “A direct hybrid finite element wave based modelling technique for efficient coupled vibro-acoustic analysis”, *Comput. Meth. Appl. Mech. Eng.*, **200**(5), 742-764.

Hall, W. S. (1994), *Boundary Element Method*, Springer.

Hughes, T.J.R. (2012), *The finite element method: linear static and dynamic finite element analysis*, Courier Dover Publications.

Ihlenburg, F. (1998), *Finite Element Analysis of Acoustic Scattering*, Springer.

Kovalevsky, L., Ladevèze, P., Riou, H. and Bonnet, M. (2012), “The variational theory of complex rays for three-dimensional helmholtz problems”, *J. Comput. Acoust.*, **20**(4), 1-25.

Kovalevsky, L., Riou, H. and Ladevèze, P. (2014), “A Trefftz approach for medium-frequency vibrations of orthotropic structures”, *Comput. Struct.*, **143**, 85-90.

Ladevèze, P. and Riou, H. (2005), “Calculation of medium-frequency vibrations over a wide frequency range”, *Comput. Meth. Appl. Mech. Eng.*, **194**(27-29), 3167-3191.

Ladevèze, P., Arnaud, L., Rouch, P. and Blanzé, C. (2001), “The variational theory of complex rays for the calculation of medium-frequency vibrations”, *Eng. Comput.*, **18**(1-2), 193-214.

Ladevèze, P., Rouch, P., Riou, H. and Bohineust, X. (2003), “Analysis of medium-frequency vibrations in a frequency range”, *J. Comput. Acoust.*, **11**(2), 255-283.

Liu, Y. (2009), *Fast Multipole Boundary Element Method: Theory and Applications in Engineering*, Cambridge, University Press.

Lyon, R.H. (1975), *Statistical Energy Analysis of Dynamical Systems: Theory and Applications*, MIT Press Cambridge, MA.

Mace, B. (2003), “Statistical energy analysis, energy distribution models and system modes”, *J. Sound Vib.*, **264**(2), 391-409.

Monk, P. and Wang, D.Q.Q. (1999), “A least-squares method for the helmholtz equation”, *Comput. Meth. Appl. Mech. Eng.*, **175**(1), 121-136.

Perrey-Debain, E., Trevelyan, J. and Bettess, P. (2004), “Wave boundary elements: a theoretical overview presenting applications in scattering of short waves”, *Eng. Anal. Bound. Elem.*, **28**, 131-141.

Riou, H., Ladevèze, P. and Rouch, P. (2004), “Extension of the variational theory of complex rays to shells for medium-frequency vibrations”, *J. Sound Vib.*, **272**(1), 341-360.

Soize, C. (1998), “Reduced models in the medium frequency range for general external structural-acoustics Systems”, *J. Acoust. Soc. Am.*, **103**(6), 3393-3406.

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

Tezaur, R., Kalashnikova, I. and Farhat, C. (2014), “The discontinuous enrichment method for medium-frequency Helmholtz problems with a spatially variable wavenumber”, *Comput. Meth. Appl. Mech. Eng.*, **268**, 126-140.

van der Heijden, A.M.A. (2009), *WT Koiter's Elastic Stability of Solids and Structures*, Cambridge University Press Cambridge.

Ventsel, E. and Krauthammer, T. (2001), *Thin Plates and Shells: Theory, Analysis, and Applications*, CRC press, Basel.