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

Nowadays, the interest of aerospace and automotive industries on virtual testing of mediumfrequency 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 mediumfrequency 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.


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 Corresponding author, Ph.D. Student, E-mail: alessandro.cattabiani@gmail.com a Professor, E-mail: andrea.barbarulo@ecp.fr b Professor, E-mail: riou@lmt.ens-cachan.fr c Professor Emeritus, E-mail: ladeveze@lmt.ens-cachan.fr 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 ultraweak 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.

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).

The equilibrium equations
The general reference example is presented in Fig. 1. The focus is on a generic subdomain 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 symbol refers to a generic boundary of where condition 1 is applied. In the particular case of a boundary shared among subdomains, is used instead. In the same way, for the conditions applied on corners, a symbol is used. The generic corner shared among subdomains is indicated with . The over-line symbol indicates that quantity is known (i.e. the value of the boundary constraint). The term ̂ is the normal unit vector of the boundary directed outward . 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 ( Fig. 2(d)). The displacement constraint , along can be divided in inplane and out-of-plane components 1 . In the same way, the load per unit length [ ] along can be divided in in-plane and out-of-plane components. The rotation condition ̂ is imposed along ̂ while a bending moment per unit length is applied along . The corners of the subdomain are subject to out-of-plane displacements constraints on and punctual forces on . Coupling conditions are applied on and along , 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 where √ is the imaginary unit, is the angular frequency and is the time.
The geometry of the subdomain can be approximated by its projection on the local plane defined by the orthonormal basis * ̂ ̂ ̂ + and the displacement field can be restricted to (Kirchhoff's kinematics assumptions)

[ ]
where is the displacement thorough the thickness of the shell, , and are respectively the total, the in-plane and the out-of-plane displacements of the middle surface and is the curvature matrix.
* + is the set of fields that satisfies the equilibrium equations.
finite energy displacement set, * + finite energy generalized stress set, , -, , -, , -( ), , , where and are Hooke's plane stress operators concerning in-plane and out-of-plane stresses respectively, is the density, is the shell thickness, and are the Young moduli along directions and respectively, and are the relative damping coefficients of the Young moduli, and are the Poisson's ratios ( ) and ( ) respectively, is the in-plane shear modulus, is its specific damping coefficient, is the inner matrix product operator, , -( ) is the symmetric part operator, ( ) is the component ( ) of the matrix , and and 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 associated with homogenized conditions ( ) is denoted as * +. This definition will be useful in the next Sections.

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 ( or ), 2. an out-of-plane condition, either a displacement constraint or a load per unit length ( or ), 3. either a rotation or a bending moment per unit length ( ̂ or ), 4. an out-of-plane condition, on corners either a displacement constraint or a point load ( or ). Boundary conditions in the most general case are  where ̂ is the tangent unit vector, is the index relative to the other subdomains sharing with the boundary , and is their total number excluded) as is shown in Fig. 3. Corner conditions in the most general case are where ̂ and ̂ are outward normal unit vectors of the two boundaries of sharing the corner . ̂ and ̂ are their respective tangent unit vectors directed towards the corner. As for , the index is relative to other subdomains sharing with the corner . is their total number ( excluded) as is shown in Fig. 4.

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 * + where , is the index related to the subdomain such that where 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.

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 where ̂ is the wave vector, ̂ is the unit direction vector of the wave vector and is a complex scalar term, is the amplitude of the plane wave, ̂ is the unit direction of the plane wave, and is the relative position vector. The relative coordinate system is located in the geometric center 1 of the element and ̂ is directed along the normal surface. The unknown parameters are * ̂ +. 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 ̂ 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 ̂ that, in case of propagative waves, is is the rotation matrix and , ) 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 where is a different discretized angle in the complex domain. Substituting the generic shape function defined in Eq. (26) with a given unit direction vector ̂ in Eqs. (3) and (4) leads to a linear set of equilibrium equations that in matrix form is (30) where ( ) is a matrix that depends on . In order to get untrivial results ( ) the dispersion equation must be imposed , -.
It provides the values of for the specified unit direction vector ̂ . When is extracted from Eq. (31) the term ̂ can be obtained by re-injecting in Eq. (31) and enforcing that 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.

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 to address orthotropic materials in plates ̂ , where In this formulation ( ). This dependency can lead to numerical difficulties during calculation of and ̂ because ̂ is no more a unit vector. For example, let us consider the case where , -. If the direction vector corrected for orthotropic materials injected in the equilibrium equations is ̂ , where , the wavenumber provided by these equations should be ( ). Since this value is computed numerically, the round-off error could approximate the wavenumber as 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 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 and ̂ .

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.

Description of the vibrational problem
Geometry and boundary conditions are illustrated in Fig. 7. It is divided in fourteen subdomains . Spars , 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 . 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 , -. Thicknesses and material properties are reported in Table 1.

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 , where is a characteristic dimension (0.5 m in this numerical example) of the considered problem, is the characteristic mesh element length, and 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  Table 2 reports principal characteristics of the used workstation.

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.

Remarks
The discrepancy between VTCR and FEM is ( ) . 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.

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.