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

VTCR is extended to shallow shell theory.The uniqueness property of the shallow-shell-VTCR solution is demonstrated.p-Refinement is the most performing.VTCR outperforms FEM at mid-frequency. 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.


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][7][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 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.

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.

The equilibrium equations
The general reference example is presented in Fig. 1. The focus is on a generic subdomain X 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 @ Ã X i symbol refers to a generic boundary of X i where condition Ã is applied. In the particular case of a boundary shared among subdomains, C is used instead. In the same way, for the conditions applied on corners, a symbol @@ Ã X i is used. The generic corner shared among subdomains is indicated with C. The over-line symbol Ã indicates that quantity Ã is known (i.e. the value of the boundary constraint). The termn i is the normal unit vector of the boundary directed outward X 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 g i (Fig. 2d). The dis- in-plane v i and out-of-plane w i components. 1 In the same way, the load per unit length in-plane b i and out-of-plane q i components. The rotation condition w i;n i is imposed along @ w i;n i X i while a bending moment per unit length m i is applied along @ m i X i . The corners of the subdomain are subject to out-of-plane displacements constraints w Ci on @@ w Ci X i and punctual forces q Ci on @@ q Ci X i . Coupling conditions are applied on C and along C, 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 jxt where j ¼ ffiffiffiffiffiffiffi À1 p is the imaginary unit, x ¼ 2pf 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 basiŝ x i ;ŷ i ;ẑ i f g and the displacement field can be restricted to (Kirchhoff's kinematics assumptions) where u z i is the displacement thorough the thickness of the shell, u i ; v i and w i are respectively the total, the in-plane and the out-of-plane displacements of the middle surface and R i is the curvature matrix.
where D i H i is the Hooke's operator for plane stress, . i the density, g i the damping coefficient, h i the shell thickness, E 0i the Young modulus, m i the Poisson's ratio, Ã : Ã the inner product matrix operator, and N i and M i are the stress and stress moment resultants tensors respectively. The sub-space of D i associated with homogenized con- This definition will be useful in the next Sections.

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 (v i or b i ), 2. an out-of-plane condition, either a displacement constraint or a load per unit length (w i or q i ), 3. either a rotation or a bending moment per unit length (w i;n i or m i ), 4. an out-of-plane condition, on corners either a displacement constraint or a punctual load (w Ci or 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 wheret i is the tangent unit vector, j C is the index relative to the other subdomains sharing with X i the boundary C, and n C is their total number (X i excluded) as is shown in Fig. 3. The possible conditions on corners in the most general case are wheren 1i andn 2i are the outward normal unit vectors of the two boundaries of X i that share the corner C.t 1i andt 2i are their respective tangent unit vectors directed toward the corner. As for j C , the index j C is relative to other subdomains sharing with X i the corner C. n C is their total number (X i excluded) as is shown in Fig. 4.

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 where Ã H is the hermitian operator. The operator I Ã f g extracts the imaginary part of Ã. These will be extensively used in the following Sections.

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 where L i and 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 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 where B; a and u are a generic matrix, vector, and scalar respectively. Such identities can be coupled with the divergence theorem to obtain the following equalities Þds: ð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 where the following equality is implicitly used which holds since N i is symmetric by definition. The last term can be substituted using the in-plane equilibrium Eq. (3) 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 The last term can be further developed introducing the identity expressed in Eq. (28) Using Eq. (4) The first term of the right side relative to the integral along the boundary can be further developed. In fact, by definition The term À H The is a continuous function along boundaries that jumps on corners. Thus: Back-substituting all these developments into Eq. (24) and simplifying leads to where the left hand side of Eq. (38) are the boundary residues and L 0e is the virtual work of external forces without considering the distributed loads 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 where L ge is the virtual work of distributed loads g i .

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 ¼ u si ; N si ; M si f g 2 D i where i 2 ½1; . . . ; n is the index related to the subdomain X i such that where D sd0i is the test functions space being the VTCR a Galerkin method.

Fundamental properties
In this Subsection uniqueness and existence properties are demonstrated.

Uniqueness
This property is demonstrated by contraction. Assume that two different solutions ; M s2 f g exist over the whole domain of the shallow shell problem and consider their difference D sd ¼ u sd ; N sd ; M sd f g . 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 where i is the index of the generic subdivision X i of the domain X.
Since this must be true for every test function, the weak form can be rewritten, without loss of generality, such that , the weak form (42) can be rearranged such that 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 Since the dissipation energy is positive-definite, the deformations N sd and M sd are null. Back-substituting this result in Eqs. (3) and (4) (with g sd ¼ 0 because this is the difference solution) implies that also the total displacement of the difference solution u sd must be zero. Thus the two solutions coincide and the uniqueness property is proven.

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.

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 where k ¼ k iki is the wave vector,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,ĉ 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 2 of the element andẑ is directed along the normal surface. The unknown parameters areĉ 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 vectork 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 vectork i that, in case of the propagative waves, iŝ is the rotation matrix and h i 2 ½0; 2pÞ is the discretized angle of the plane wave direction. In case of evanescent waves the unit direction vector is defined such that where Þ is another discretized angle in the complex domain. Indeed out-of bound values of / 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 ffiffiffi 2 p is selected. Substituting the generic shape function defined in Eq. (45) with a given unit direction vectork i in Eqs. (3) and (4) leads to a linear set of equilibrium equations that in matrix form is where Z i ¼ Z i ðk i Þ is a 3 Â 3 matrix that depends on k i . In order get not trivial results (a i ¼ 0) the dispersion Equation must be imposed It provides the values of k i for the specified unit direction vectork i . When k i is extracted from Eq. (50), the termĉ i can be obtained by re-injecting k i in Eq. (50) and enforcing that

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.

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.
The curved shallow shell is partially simply supported and is subject to a time-dependent punctual load directed alongẑ. The reference solution approximates it by a series of sines   (b) VTCR solution. Fig. 6. Comparison of the out-of-plane displacements of the vibrational problem explained in Section 4.1.
C mn sinðk n xÞ sinðl m xÞ where x ¼ 2pf ; g z is the load alongẑ, and dðx À pÞ is the Dirac delta function. The out-of-plane displacement w of the analytic solution is approximated by a Fourier series B mn sinðk n xÞ sinðl m xÞ The series is stopped at fn max ¼ 100; m max ¼ 100g because convergence is reached. Its graphical representation is given in Fig. 6 alongside the VTCR solution.
In the VTCR the punctual force is approximated by a complete Fourier approximation: Injecting it in Eqs. (3) and (4) and erasing the in-plane inertia produces the particular solution of the VTCR. 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 VTCR is led by the approximation of the particular solution. Kirchhoff theory allows errors of the order of magnitude h Rx ¼ 10 À3 . Considering this error threshold, the series is stopped at fn VTCRmax ¼ 40; m VTCRmax ¼ 40g. The error indicator is based on the kinetic energy of the out-of-plane displacement where w ref ¼ w ref ðx; yÞ is the out-of-plane displacement of the reference solution and w VTCR ¼ w VTCR ðx; yÞ is the out-of-plane displacement of the VTCR solution.  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 where n rtot is the total number of rays, n e is the number of elements, n rp is the number of propagative rays, n Ce is the number of boundaries per element e, and n Ce is the number of corners per subdomain.
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.    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 chose FEM element size h 0 is 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 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.

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 : where u ABAQUS is the displacement calculated by ABAQUS Ò and u VTCR is the displacement calculated by the VTCR theory. The total number of rays is 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.

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 performaces 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. Table 4 DoFs, time, and memory consumption comparisons between VTCR and FEM of the numerical example described in Section 4.2.