The Variational Theory of Complex Rays: An answer to the resolution of mid-frequency 3D engineering problems

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.


Introduction
Due to regulatory requirements and growing customer expectations, the vibrational and acoustical behavior of structures has become an important issue in engineering design. Virtual prototypes are being increasingly used in order to reduce the number of (expensive) experimental tests. This requires efficient and accurate numerical prediction techniques.
Today, most engineering problems are solved using the Finite Element Method (FEM). However, because the FEM seeks an approximate solution by using continuous, piecewise polynomial functions, it is incapable of handling problems involving rapidly oscillating functions other than by using highly refined meshes. High-order finite elements can be used as an alternative to a very fine mesh. However, due to the pollution effect, the number of degrees of freedom must still be appropriately increased for the standard FEM to maintain the desired level of accuracy. Thus, in practice, the FEM is limited to low-frequency applications. Although various techniques have been proposed in order to overcome this limitation (see [1][2][3][4][5][6][7]). Indeed, these techniques have been shown to reduce computation costs compared to standard finite elements, but the frequency ranges which they are capable of addressing are still lower than that would be desirable.
Another often used prediction technique is the Boundary Element Method (BEM). Since the BEM is based on a boundary integral formulation of the problem, only the boundary of the domain being considered needs to be discretized. Within the boundary elements which compose the model, the variables are expressed in terms of simple polynomial shape functions. Unlike the FEM, the enforcement of the boundary conditions results in only a small system of equations to be solved for the nodal values at the discretized boundary. Once these nodal values have been determined, the field variables within the domain can be reconstructed in a postprocessing step using the boundary integral formulation. Since only the boundary needs to be discretized, the BEM is exempt from some of the limitations of the FEM mentioned previously. Besides, many works have been dedicated to the rapid evaluation of the necessary integrals (see [8,9]). However, these methods are

Basic concepts of the VTCR
Let us consider a general steady-state dynamic problem in O, a d-dimensional domain with boundary qO. The problem to be solved can be expressed as: find u such that BðuÞ ¼ b over qO (2) where b denotes a prescribed function associated with prescribed sources located at the boundary of O; L is a differential operator associated with the governing equation; and B is a differential operator associated with the boundary conditions (and also, if O is partitioned into subregions, with the continuity across the interfaces between these subregions).
The first characteristic feature of the VTCR is that it uses a global formulation of boundary conditions (2) in terms of both primal and dual quantities. Thus, Problems (1) and (2) are expressed as: find u 2 S ad such that ðu,duÞ B ¼ ðb,duÞ 8du 2 S ad (3) where ð,Þ B and ðb,Þ are respectively a bilinear form and a linear form equivalent to (2), and S ad is the space of the functions which satisfy (1). Formulation (3) is written in such a way that approximations which are a priori independent of one another can be used in the different subregions of O.
Then, all that is required to build a VTCR approximation is the definition of a subspace S a ad of S ad . The approximate formulation can be written as: find u a 2 S a ad such that The second characteristic feature of the VTCR is the description of the space S a ad using two-scale approximations with a strong mechanical content: the solution of (1) is considered to be the superposition of an infinite number of plane waves which satisfy the governing equation exactly. Then, the construction of space S a ad is the result of an approximation of the distribution of the amplitudes of the waves, which can be chosen among different options in order to benefit from physical and numerical advantages.
Once the discretization of each subregion of O has been chosen, the VTCR leads to a linear system of equations where matrix K corresponds to the discretization of the bilinear form of weak formulation (4), A is the vector of the unknown quantities which approximate the distribution of the wave amplitudes, and F corresponds to the discretization of the linear form of (4) projected onto S a ad . The approximate VTCR solution of Problem (1) and (2) is obtained by selecting from S a ad the functions whose amplitudes are A.
Several remarks can be made: Weak formulation: The VTCR is based on a specific weak formulation of the problem which was developed in order to enable the approximations within the subregions to be a priori independent of one another and which takes into account both the boundary conditions and the continuity across the interfaces, thus eliminating the need for a special treatment to guarantee interelement continuity. Any type of shape function can be used within each subregion, which gives this approach great flexibility and robustness.
Injected waves: Space S a ad is spanned by propagative and evanescent waves which are viewed as two-scale approximations. Only the slowly varying scale (which corresponds to the amplitudes of the waves) is discretized. The rapidly varying scale (which corresponds to the spatial shape of the waves) is taken into account analytically in Formulation (4). These two scales give the approximation a strong physical content and several computational advantages.
Dimension of subspace S a ad : The number of waves to be used in S a ad was studied in [30] for plates, in [32] for shells and in [37] for acoustics. Each of these studies showed that the VTCR, along with all other Trefftz methods, converges rapidly when the number of wave DOFs increases. Heuristic criteria and energy criteria for determining a priori the number of propagative waves required for an accurate numerical solution were presented in the last two references. Each of these criteria involves a relation between that number of waves and the ratio between the wavelength and the characteristic length of the domain being studied. Concerning evanescent waves, [30] showed that these are essential, but that not many of them are needed. If they are related to higher wavelengths than the propagative waves (e.g. in the case of pure acoustics problems), they are not even introduced at all.
Properties of matrix K: Matrix K in (5) is complex, nonsymmetric and frequency-dependent. However, since the shape functions are defined over acoustic subregions O e , K is populated by blocks. Each block K ½i,i associated with subregion O i is fully populated, but the off-diagonal submatrix K ½i,j between subregions O i and O j is zero if these subregions are not connected.
Accuracy of derived variables: In the FEM, the primary response variables are usually approximated by simple polynomial shape functions. Therefore, the quantities formed by derivation of these variables (such as stresses and strains, acoustic velocities, etc.) are less accurate than the variables themselves. In the VTCR, there is no loss of accuracy because the derivatives of the wave functions are also wave functions.
Treatment of problems with complex geometries: The VTCR can handle complex geometries, whose effects are included in (3), but sometimes a proper definition of S ad (and, thus, of S a ad ) requires a specific discretization of O in the subregions O e .
Computational performance: The calculation of the matrix coefficients involves the integration of highly oscillatory functions. Special care must be taken in carrying out these integrations in order to ensure good convergence. Semianalytical methods can be used in the case of meshes composed of straight lines. Numerical integration schemes can be used for curved boundaries. Graphics Processing Units (GPUs) or parallel computer systems with multiple servers can also be used to accelerate many of the computations without resorting to low-level programming.
Handling of internal sources: The presence of an internal source in (1) can be taken into account simply by adding a particular solution in approximation spaces S ad (see [29] for further details and some numerical examples).

Problem description
In order to simplify the presentation of the VTCR for 3D plate assemblies, the problem will be formulated for an assembly of two plates, but this can be easily generalized. Let us consider two isotropic homogeneous plates whose 3 frequency o are governed by the following partial differential equations: where E, n, h, K PS , r, Z and X denote respectively Young's modulus, Poisson's coefficient, the plate's thickness, Hooke's plane stress tensor, the density, the damping coefficient and the curvature operator. Eqs. (6) and (7) where W and M are the sets of the finite-energy displacements and moments. The problem of Fig. 1 can be expressed as: find ðw,MÞ 2 S ad and ðw 0 ,M 0 Þ 2 S 0 ad which satisfy the boundary conditions: Problem (9) can also be written as: find ðw,MÞ 2 S ad and ðw 0 ,M 0 Þ 2 S 0 ad such that where n and Refg denote the conjugate part and the real part of the complex quantity . Eqs. (8) and (9) correspond to (1) and (2) for an assembly of two plates, and Formulation (10) corresponds to (3). Complete details about Formulation (10), especially its equivalence with (8) and (9), can be found in Ref. [29].

Approximate solution using plane waves
Among the several possible ways to generate the space S a ad in which to seek an approximate solution (see (4)), the expansion proposed in Ref. [29] seems promising. One assumes that the solution is correctly represented by a specific set of propagative and evanescent waves. Thus, S a ad is generated by a set of functions of the form e ik:x , where i is the complex imaginary unit, k is the wave vector and x is the spatial location. This type of function satisfies (6) and (7) if, and only if See Ref. [29] for more details. Depending on the real and imaginary parts of the coordinates of k, different types of waves satisfy (11): interior waves, edge waves or corner waves (see [29]). Fig. 2 gives a representation of such waves for a homogeneous rectangular plate. Then, S a ad is spanned by the entire set of the interior, edge and corner waves S a ad ¼ span e ik int,j :x ,e ik edge,h :x ,e ik corner,l :x È j 2 ½1; n int ,h 2 ½1; n edge ,l 2 ½1; n corner É where n int , n edge and n corner are the numbers of waves of each type which are being considered in the approximate solution.
The choice of these numbers depends on the problem and on the level of accuracy desired. Fig. 3 shows the geometry and the mechanical properties of the structure being studied, which is a typical stringer as can be found in a car chassis. The reference solution was calculated using the FEM program MSC Nastran with about seven CQUA4 elements per wavelength. In each plate, the VTCR solution was calculated with 32 interior waves (distributed regularly in all the possible 2-D propagation directions of the reference surface), nine edge waves per edge and no corner waves. The results obtained with the two methods are shown in Fig. 4.

Numerical example
One can see that the VTCR solution and the FEM solution are very similar. The spatial distributions and the amplitudes of the vibrational behavior are nearly the same. Moreover, the maximum relative differences between the VTCR and FEM effective displacements ð1=mesðO i ÞÞ R O i 9w9 dO in each plate O i (i¼1: :12) are 0.1, 0.03 and 0.07 at 400 Hz, 800 Hz and 1500 Hz respectively. However, significantly fewer degrees of freedom were necessary with the VTCR, especially at higher frequencies (about two hundred times fewer at 1500 Hz). This example clearly demonstrates the benefits of the VTCR in dealing with 3D plate assemblies.

Problem description
Again, in order to simplify the presentation, the problem will be formulated for an assembly of two shells, but could be easily generalized to a larger number. Let us consider two isotropic homogeneous shells whose reference surfaces are O and O 0 with boundaries qO and qO 0 respectively. In each shell, the quantities of interest are the displacement vector  The formulation uses Koiter's linear shell theory (see Refs. [38,39]). The displacement class is restricted to u ¼ v þ we 3 þ zb with b ¼ ÀrðwÞÀBv, z being the thickness, e 3 the normal to the reference surface and B the curvature tensor (Kirchhoff's kinematic assumption). The transverse strain energy is neglected.  For structure O, let us define the space S ad ¼ fðv,w,N,MÞ 2 V Â W Â N Â Mg (V, W, N and M being the finite-energy sets of the corresponding physical quantities) such that div NÀBðdiv MÞ ¼ Àro 2 hv where K CP is Hooke's plane stress operator, e is the gradient operator, r is the density, Z is the damping coefficient and h is the thickness of the shell. S 0 ad is defined in the same way. Thus, the problem to be solved can be expressed as: find ðv,w,N,MÞ 2 S ad and ðv 0 ,w 0 ,N 0 ,M 0 Þ 2 S 0 ad which verify the boundary conditions Problem (14) can also be written as: find ðv,w,N,MÞ 2 S ad and ðv 0 ,w 0 ,N 0 ,M 0 Þ 2 S 0 ad such that Eqs. (13) and (14) correspond to (1) and (2) for an assembly of two shells, and Formulation (15) corresponds to (3). Complete details about Formulation (15), especially its equivalence with (13) and (14), can be found in Ref. [32].

Approximate solution using plane waves
Again, among the several possible ways to generate the space S a ad in which to seek an approximate solution (see (4)), the expansion proposed in Ref. [32] seems promising. One assumes that S a ad is correctly represented by a specific set of propagative and evanescent waves expressed as e ik:x . In order to find the equation that the wave vector k must verify, one can perform an asymptotic expansion of (13) with respect to the small parameter ffiffiffiffiffiffiffiffi h=R p , where R is the smallest radius of curvature of the shell being considered (which is assumed to be only slightly curved). Then, one can easily prove that the wave vector must verify where R is a rotation matrix in the coordinate system of the shell's reference surface. Of course, in the case of a plate, i.e. if B ¼ 0, (16) boils down to (11). Otherwise, the additional portion of Eq. (16) is necessary in order to calculate waves in shells because it is used to represent some well-known characteristics of propagative waves in certain directions at certain frequencies (see [32] for details). Again, (16) can be satisfied by interior waves, edge waves or corner waves depending on the real and imaginary parts of the coordinates of k. Fig. 6 illustrates these different types of waves in the case of a homogeneous cylinder. Then, S a ad is spanned by the entire set of the interior, edge and corner waves S a ad ¼ spanfe ik int,j :x ,e ik edge,h :x ,e ik corner,l :x j 2 ½1; n int ,h 2 ½1; n edge ,l 2 ½1; n corner g (17) where n int , n edge and n corner are the numbers of waves of each type which are being considered in the approximate solution. Again, the choice of these numbers depends on the problem and on the level of accuracy desired. Fig. 7 shows the geometry of a 3D assembly of two cylindrical shells and one plate. The mechanical properties are

Numerical example
The cylinder's height L z and its curvature are both equal to 1 m and its thickness is h ¼ 0:01 m. We used the solution obtained with the FEM program MSC Nastran as the reference solution. The mesh seed used in the program was set to create 10 QUADR shell elements per wavelength and the solution was obtained with about 1; 200,000 DOFs. For the VTCR solution, the structure was divided into three parts (shell O 1 , shell O 3 , and plate O 2 , see Fig. 7).
The model involved 264 DOFs (corresponding in each substructure to 24 interior waves, five edge waves per edge and no corner waves). One can verify in Fig. 8 that the result given by the VTCR is quite good and very similar to that given by the FEM regarding the distribution and amplitudes of the displacements).
In order to compare the solutions in terms of an effective quantity, we calculated the effective displacement

Problem description
Let us consider a general steady-state interior dynamic problem in a 3D acoustic cavity O filled with a fluid characterized by its sound velocity c, its density r and its damping coefficient Z. The problem of the steady-state dynamic

Table 1
Comparison of the effective displacements in substructures O i , i ¼ 1: :3 of the problem of Fig. 7, whose solutions are given in Fig. 8 where k ¼ o=c is the wavenumber, Z is an impedance coefficient and L v ð&Þ an operator defined as L v ð&Þ ¼ ði=roÞðq&=qnÞ ¼ ði=roÞn Á rð&Þ, n being the outward normal to qO. h d , p d and v d denote respectively a prescribed excitation over q Z O, a prescribed pressure over q p O and a prescribed velocity over q v O. Let us partition O into n el non-overlapping subcavities O E and use the notations The problem defined by (18) and (19) is illustrated in Fig. 9.

Approximate solution using plane waves
Again, there are several possible ways to generate the space S E,a ad in which to seek an approximate solution (see (4)). In the cases of plates and shells of Sections 3 and 4, although the structure was defined in 3D, the waves propagated in 2D, i.e. in the plate's or the shell's reference surface. Here, the difficulty comes from the fact that waves can propagate in all 3D directions. Taking into account all possible directions can lead to a huge problem. Therefore, the discretization of these directions must be carried out efficiently.
In order to do that, a direct and straightforward extension of what was proposed for 2D acoustics in Ref. [37] seems promising. The pressure field is sought as an integral distribution of plane waves, also known as Herglotz wave functions, which satisfy the homogeneous Helmholtz equation. In the 3D case, this plane wave distribution is defined over the unit sphere: jÞ Á e ikðy,jÞ:x dy sin j dj (22) where A E denotes the amplitude of the plane waves in the 3D spherical direction ðy,jÞ. One can easily prove that (22) satisfies (20) and, therefore, that p E ðxÞ belongs to S E ad . In order to preserve all possible wave directions, avoid too refined a quadrature of the unit sphere and benefit from the advantages of the Fourier decomposition introduced in Ref. [37], function A E ðy,jÞ is described using a truncated Laplace series, which is a 3D extension of a Fourier series. Then, the discrete space S E,a ad of (21) can be expressed as S E,a ad ¼ span where kðy,jÞ is the wave vector of the plane wave propagating in the spherical direction ðy,jÞ and Y m l ðy,jÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 Á ðlÀmÞ!=ðl þmÞ! p Á P m l ðcos yÞ Á e imj is the spherical harmonic of nonnegative index l and momentum m, and L E the maximum index retained in the discrete space S E,a ad of dimension ðL E þ 1Þ 2 . In the last expression, P m l ðXÞ ¼ ððÀ1Þ m =2 l Á l!Þð1ÀX 2 Þ m=2 ðq m þ l ðX 2 À1Þ l =qX m þ l Þ denotes the Legendre polynomial and parameter m varies from Àl to l. The set fY m l ðy,jÞg 9m9 r l r 1 constitutes a complete orthogonal coordinate system for the unit sphere with regard to the classical L 2 scalar product.
The boundary conditions consisted of a unit external pressure, a prescribed impedance Z ¼ 425 Pa s=m or a rigid-wall condition as indicated in the figure. The circular frequencies considered were o ¼ 2p Á 500 Hz, 1000 Hz and 1500 Hz. For each frequency, a reference solution (in the form of a pressure distribution p ref ) was calculated using the FEM program COMSOL. In order to achieve an appropriate quality of results, the mesh seed used in the program was set to create 14 3D TETRA elements per wavelength, leading to three models with 7; 160, 57; 230 and 192; 420 DOFs respectively for the three frequencies.
For the VTCR solution, the structure was divided into three subcavities O 1 , O 2 and O 3 (see Fig. 11). Fig. 12 shows for each frequency the real part of the VTCR pressure, the number of DOFs used and the corresponding relative error  given by the VTCR and the FEM are very similar, but with the VTCR the solution was obtained using far fewer DOFs (over 100 times fewer at o ¼ 2p Á 150 Hz). Fig. 13 shows another interesting comparison concerning the global Frequency Response Functions (FRFs) in energy for the FEM and the VTCR. In addition, the number of DOFs used is indicated on the vertical scale on the right. Again, in the case of the FEM, 14 elements per wavelength were used in order to ensure the accuracy of the solution (especially for the 3D calculations). In the case of the VTCR, the number of shape functions was chosen so that the relative error compared to the FEM would be less than 10 À3 for each frequency. One can see that the two FRFs thus predicted are very similar, which is normal since the number of DOFs for the VTCR was chosen specifically so the error between the two curves would be very small. However, as mentioned before, the number of DOFs needed for the VTCR was considerably smaller than for the FEM. Furthermore, it is clear that the need for large numbers of DOFs increases faster for the FEM than for the VTCR. Again, this leads to the conclusion that the VTCR is a competitive numerical technique for the analysis of 3D acoustic problems.   Fig. 11. The circular frequency considered, the number of DOFs and the relative error are indicated in each case.  Fig. 11. The corresponding numbers of DOFs are indicated on the vertical scale on the right. For the FEM, 14 elements per wavelength were used. For the VTCR, the number of shape functions was chosen so that the relative error compared to the FEM would be less than 10 À3 for each frequency.

Conclusion
This paper presents an overview of the development of the VTCR for the resolution of mid-frequency 3D engineering problems. First, the paper reviews the basic concepts and main characteristics of the VTCR, a new variational formulation in which approximations within subregions are defined a priori independently of one another, and a two-scale approximation with a strong mechanical content is used to derive the approximation spaces. Following that introduction, the details of the resolution of problems involving 3D assemblies of plates or shells as well as acoustic problems are given.
For all these problems, the discussion of the VTCR approach follows the same outline: the problem is defined, the expansion of the field variables is presented and a numerical example is given. Such a systematic presentation was possible because the VTCR is a general-purpose predictive numerical technique for mid-frequency problems. Throughout these numerical examples, we show that the VTCR is capable of finding an accurate solution with great numerical benefits compared to classical deterministic numerical techniques, such as the FEM.
Some possible subjects for additional work are the investigation of vibroacoustic problems and the development of an efficient and robust numerical strategy for large calculations in a frequency range.