Hao Li
email: hao-li@lmt.ens-cachan.fr

Pierre Ladevèze
email: pierre.ladeveze@lmt.ens-cachan.fr

Hervé Riou
email: herve.riou@lmt.ens-cachan.fr

Hybrid Finite Element Method and Variational Theory of Complex Rays for Helmholtz Problems

Keywords: Hybrid methods, Trefftz methods, discontinuous Galerkin methods, VTCR, FEM

In this paper, we consider the Weak Trefftz Discontinuous Galerkin (WTDG) method, which enables one to use at the same time the Finite Element Method (FEM) or Variational Theory of Complex Rays (VTCR) discretizations (polynoms and waves), for vibration problems. It has already been developed such that the FEM and the VTCR can be used in different adjacent subdomains in the same problem. Here, it is revisited and extended in order to allow one to use the two discretizations in the same subdomain, at the same time. Numerical examples illustrate the performances of such an approach.

Introduction

Today, the numerical simulation has become indispensable to analyze and optimize the problems in every part of engineering processes. Based on the piecewise, continuous polynomial approximation, the standard Galerkin Finite Element Method (FEM) 1 is a wellestablished method which is commonly used for the analysis of time-harmonic dynamic problems. However, it suffers from the pollution effect. [START_REF] Deraemaeker | Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions[END_REF] As a consequence, when the wave number is increased, the meshes need drastically to be refined. This leads to inefficient calculations for problems when the medium and the high frequencies are considered. As the FEM is based on polynomial approximations, the problem cannot fundamentally be solved.

Other approaches have been developed in order to alleviate this problem, namely the Trefftz approaches. [START_REF] Trefftz | Ein gegenstück zum ritzschen verfahren[END_REF] They are based on the use of exact approximations of the governing equation. Such methods are, for example, the partition of unity method, [START_REF] Strouboulis | Partition of unity method for Helmholtz equation: q-convergence for plane-wave and wave-band local bases[END_REF] the ultra weak variational method, [START_REF] Cessenat | Application of an ultra weak variational formulation of elliptic PDEs to the two-dimensional Helmholtz problem[END_REF][START_REF] Huttunen | An ultra-weak method for acoustic fluid-solid interaction[END_REF] the least square method, [START_REF] Monk | A least-squares method for the Helmholtz equation[END_REF][START_REF] Gabard | A comparison of wave-based discontinuous Galerkin, ultra-weak and least-square methods for wave problems[END_REF] the plane wave discontinuous Galerkin methods, [START_REF] Gittelson | Plane wave discontinuous Galerkin methods: Analysis of the h-version[END_REF] the method of fundamental solutions [START_REF] Fairweather | The method of fundamental solutions for elliptic boundary value problems[END_REF][START_REF] Barnett | Stability and convergence of the method of fundamental solutions for Helmholtz problems on analytic domains[END_REF] the discontinuous enrichment method, [START_REF] Farhat | The discontinuous enrichment method[END_REF][START_REF] Farhat | A domain decomposition method for discontinuous Galerkin discretizations of Helmholtz problems with plane waves and Lagrange multipliers[END_REF] the element free Galerkin method, [START_REF] Bouillard | Element-free Galerkin solutions for Helmholtz problems: Fomulation and numerical assessment of the pollution effect[END_REF] the wave boundary element method [START_REF] Perrey-Debain | Wave boundary elements: A theoretical overview presenting applications in scattering of short waves[END_REF][START_REF] Bériot | Plane wave basis in Galerkin bem for bidimensional wave scattering[END_REF] and the wave based method. [START_REF] Desmet | An indirect Trefftz method for the steady-state dynamic analysis of coupled vibro-acoustic systems[END_REF][START_REF] Van Genechten | An efficient wave based method for solving Helmholtz problems in three-dimensional bounded domains[END_REF] The Variational Theory of Complex Rays (VTCR), first introduced in Ref. 19, belongs to this category of numerical strategies which uses waves in order to get some approximations for vibration problems. It has been developed for 3D plate assemblies in Ref. 20, for plates with heterogeneities in Ref. 21 and for shells in Ref. 22. Its extensions to acoustics problems can be seen in Refs. 23 and 24. In opposition to FEM, the VTCR has good performances for medium frequency applications, but is less efficient for very low frequency problems.

Recently, a new approach called the Weak Trefftz Discontinuous Galerkin (WTDG) method has been introduced. [START_REF] Ladevèze | On Trefftz and weak Trefftz discontinuous Galerkin approaches for medium-frequency acoustics[END_REF] It differs from the pure Trefftz methods, because the necessity to use exact solution of the governing equations can be weaken. Then it allows one to couple the FEM (based on polynoms) and the VTCR (based on waves) approximations at the same time. In the illustrations proposed Ref. 25, the polynoms and the waves approximations where used in different adjacent subdomains of a problem. Here, we propose to mix the approximations in the same subdomains, at the same time. This allow us to give to the strategy a strong flexibility for engineering problems. Through this paper, it will be shown that such a mix approach (which uses at the same time the polynomial and the wave descriptions) presents better performances than a pure FEM approach (which uses only a polynomial description) or a pure VTCR approach (which uses only a wave description).

This work is presented as follows: in Sec. 2 we present the reference problem and develop some fundamental properties of WTDG. In Sec. 3, we propose some numerical examples to illustrate the strategy. Conclusions are drawn in Sec. 4.

The Weak Trefftz Discontinuous Galerkin Method

The reference problem

The reference problem to solve is illustrated in the context of Helmholtz problem, which can be applied into acoustic vibration problems. Let us consider the domain Ω with the boundaries ∂Ω = ∂ 1 Ω ∪ ∂ 2 Ω, where ∂ 1 Ω corresponds to the boundary where Dirichlet conditions are prescribed, and ∂ 2 Ω corresponds to the boundary where Neumann conditions are prescribed (see Fig. 1). The following problem is considered: such that

find u ∈ H 1 (Ω) Ω Ω E Γ EE r d Ω u d ∂ 1 Ω ∂ 2 Ω g d Ω E
(1 + iζ)∆u + k 2 u + r d = 0 over Ω, u = u d over ∂ 1 Ω, (1 + iζ)∂ n u = g d over ∂ 2 Ω, (1)
where ∂ n u = grad u • n and n is the outward normal. u is the physical variable studied such as the pressure in acoustics. ζ is the damping coefficient, which is positive or equals to zero. k is the wave number, which is a real, positive and constant value. u d and g d are the prescribed Dirichlet and Neumann data. The first equation of ( 1) is the standard Helmholtz equation. The second and third equations correspond to the boundary conditions.

The WTDG variational formulation of the reference problem

In order to get an equivalent variational formulation of (1), the domain is divided into subdomains Ω E with E ∈ E (see Fig. 1). Γ EE denotes the interface between two subdomains E and E . Γ EE denotes the interface between subdomain Ω E and boundary ∂Ω. The WTDG method consists in using the working space U ⊂ H 1 (Ω):

U = {u | u |Ω E ∈ U E }, U E = {u E | u E ∈ V E ⊂ H 1 (Ω E )}. (2)
The vector spaces associated with U and U E where r d = 0 are denoted by U 0 and U E,0 . Denoting

{u} EE = (u E + u E ) |Γ EE , [u] EE = (u E -u E ) |Γ EE , q u = (1 + iζ)grad u (3)
the WTDG formulation can be written [START_REF] Ladevèze | On Trefftz and weak Trefftz discontinuous Galerkin approaches for medium-frequency acoustics[END_REF]

: find u ∈ U such that Re   -ik   E,E ∈E Γ EE 1 2 {q u • n} EE {ṽ} EE - 1 2 [q v ] EE [u] EE dS - E∈E Γ EE ∩∂ 1 Ω qv • n (u -u d )dS + E∈E Γ EE ∩∂ 1 Ω α • i • ṽ (u -u d )dS + E∈E Γ EE ∩∂ 2 Ω (q u • n -g d )ṽdS - E∈E Ω E (div q u + k 2 u + r d )ṽdΩ = 0 ∀ v ∈ U 0 , ( 4
)
where α is a parameter strictly positive to enforce the boundary Dirichlet condition.

As one can see, there is no a priori constraint on the choice of the spaces U and U 0 . As a consequence, one can select polynoms, like in FEM, or waves, like in VTCR, or even both. This gives to the strategy a very strong flexibility.

Properties of the WTDG formulation

Let us note that (4) can be written as:

find u ∈ U such that b(u, v) = l(v), ∀ v ∈ U 0 , ( 5
)
where b and l, are u-linear in U and v-(anti

)linear in U 0 . b is such that b(u, u) is real. Property 1. For u ∈ U 0 , we have b(u, u) = E∈E kζ Ω E grad ũ • grad udΩ + E∈E Γ EE ∩∂ 1 Ω kαuũdS ≥ 0. (6) Proof. b(u, u) = Re -ik E∈E ∂Ω E (q u • n)ũdS - E∈E Ω E div q u ũdΩ + E∈E Γ EE ∩∂ 1 Ω αiuũdS . (7)
Consequently,

b(u, u) = E∈E kζ Ω E grad ũ • grad udΩ + E∈E Γ EE ∩∂ 1 Ω kαuũdS. ( 8
)
From Property 1 it can be deduced that if b(u, u) = 0, then u is equal to zero over ∂Ω E ∩ ∂ 1 Ω. It is a piecewise constant within subdomains Ω E , E ∈ E and is not equal to zero. To keep the uniqueness of the solution, condition (P) is introduced to be satisfied by the shape functions which belong to U 0 .

Condition (P). Let a

E ∈ U E be a piecewise constant function within subdomains E ∈ E. a E satisfies condition (P) if    ∀ v ∈ U 0 , ∀ E ∈ E, Re   -ik   E,E ∈E ∂Ω E (q v • n)ã E dS     = 0    ⇒ a E = ±a, ( 9
)
where E is a subdomain sharing a common boundary with E. And let us take the convention

a E = -a E over ∂Ω E ∩ ∂Ω.
Property 2. If U 0 satisfies condition (P) and if ζ is positive, the WTDG formulation (4) has a unique solution.

Proof. Existence of solution is evident. And suppose (4) has two solutions u 1 and u

2 . v = u 1 -u 2 ∈ U 0 and b(v, v) = E∈E kζ Ω E grad ṽ • grad vdΩ + E∈E Γ EE ∩∂ 1 Ω kαvṽdS = 0. ( 10
)
It can be observed that v E = a E with E ∈ E, where a E is piecewise constant within the subdomains and a E = 0 in the subdomains sharing a common boundary with ∂ 1 Ω. Back-substituting this result into (4), one also finds b(v, v * ) = 0 ∀ v * ∈ U 0 , which leads to

∀ v * ∈ U 0 , Re   -ik   E,E ∈E ∂Ω E (q v • n) ãE dS     = 0. ( 11
)
Equation ( 11) corresponds to the condition (P), where E represents a subdomain sharing a common boundary with E, with the convention

a E = -a E over ∂Ω E ∩ ∂Ω. Consequently, a E = ±a ∀ E ∈ E.
Moreover, given that a E = 0 over ∂ 1 Ω, we have a = 0.

Property 3.

If U E is the combination of solution spaces of FEM and VTCR, then the condition (P) is satisfied, and

u 2 U 0 = b(u, u) + γ 2 (u) (12)
is a norm over U 0 , with the definition

U E,0 = {u|u ∈ V E , u = C E • X E }, ( 13
)
where C E is constant over Ω E and X E is the position vector relative to the center of inertia of element E. U 0 denotes the associated space defined over Ω of U E,0 . And for u ∈ U 0 the definition of quantity γ is introduced as

γ(u) = sup v∈U 0 b(u, v)/||C v || L 2 (Ω) , (14)
where C v corresponds to the vector C E of v according to (13).

Proof. To demonstrate that condition (P) is satisfied, let us take a E ∈ U E , E ∈ E a piecewise constant. Since U E is the combination of FEM and VTCR, u could be any linear combination of polynoms and wave functions. Therefore when u = a E a piecewise constant, u could only be the polynomial function of order 0. Let us note that

∀ E ∈ E a E = β E (β E constant over Ω E ) for any subdomain E sharing a common boundary with E. Let us introduce z E = β E + a E . z E is continuous because z E|Γ EE = a E + a E = z E |Γ EE .
It follows that z is constant over Ω. Since z is zero over ∂Ω, z = 0 over Ω and β E = -a E . Consequently, a E can be only the values of +a or -a, a being a constant over Ω.

To demonstrate that ||u|| 2 U 0 is a norm over U 0 , let us consider that u 2 U 0 = b(u, u) + γ 2 (u) = 0. It follows that b(u, u) = 0 and γ(u) = 0. From (8) it can be obtained that

u |Ω E = a E is constant over Ω E and that u = 0 over ∂ 1 Ω. Then γ(u) is equal to γ(u) = sup v∈U 0 1 ||C v || L 2 (Ω) Re   -ik   E,E ∈E ∂Ω E (q v • n)ã E dS     = 0. ( 15
)
Since condition (P) is satisfied, it can be derived that u E = ±a, a being a constant over Ω.

Finally from u = 0 on ∂ 1 Ω, one gets u = 0 over Ω.

Implementation of the WTDG method

As said before, we here seek the approximation solution as a sum of traditional FEM shape functions and traditional VTCR shape functions. As a consequence, we can write

U =    u ∈ L 2 (Ω) : u(x) |Ω j = N j n j =1 λ n j ϕ n j (x) + M j m j =1 A m j e ikdm j •x , λ n j ∈ C, A m j ∈ C, j = 1, . . . , n s    , ( 16
)
where ϕ n correspond to the standard, continuous piecewise polynomial FEM shape function, associated to the unknown degree of freedom λ n , and e ikdm•x is the standard wave VTCR shape function, associated to the unknown degree of freedom A m (d m is the unitary vector defining a direction of propagation of wave). As a consequence, after back-substituting ( 16) into (4), we get the matrix problem

K FF K FV K VF K VV • λ A = L F L V , ( 17
)
where K FF and K VV are the matrix blocks of pure FEM part and pure VTCR part, respectively. K FV and K VF denote the matrix blocks of couple terms of FEM and VTCR. λ and A are the unknown degrees of freedom to be found. The terms of K FF are calculated by the Gauss integration. The K VV , K FV and K VF terms are calculated through analytical calculations. The exponential form of the VTCR shape functions allow one to do this, at least on straight line boundaries. On the seldom cases where the boundaries are not straight, these terms are computed through numerical computations.

Numerical Examples

Homogeneous Helmholtz problem with variation of scales in the solution

The domain being considered is the square Ω = [0 m; 0.5 m] × [0 m; 0.5 m] The prescribed boundary conditions are such that the exact solution is ζ = 0.0001. We satisfy Eq. ( 1) with k which varies from 5 m -1 to 72 m -1 . This example is interesting, because it covers different scales (from slow varying scale with k = 5 m -1 to fast varying scale with k = 72 m -1 ) for the same model of problem, in the solution. It can be seen from this example how the WTDG method works from lower frequency problem to medium frequency problem. This definition of the problem and the discretization strategy can be seen in Fig. 2.

u exact = 5 j=1 exp(ik(1 + iζ) -1/2 (x • cos(θ j ) + y • sin(θ j ))) with θ 1 = 5.
The WTDG strategy is used to find an approximation of u exact . The WTDG formulation ( 4) is applied with Dirichlet boundary condition with α = 0.0001. For the FEM part, we use a regular squared mesh with a polynomial discretization of degree 1. The discretization uses 10 elements per wave length. For the VTCR, we use a 2D regular polar repartition of the directions of the waves in Ω. The discretization strategy complies with a geometrical heuristic criterion: N e = δkR e /(2π), where N e is the number of directions of waves, δ a parameter to be chosen [START_REF] Kovalevsky | On the use of the variational theory of complex rays for the analysis of 2-D exterior Helmholtz problem in an unbounded domain[END_REF] and R e is the characteristic radius of Ω. Here we take δ = 10. The convergence of the WTDG strategy is assessed by computing the relative error |u exact -u|/|u exact | = Ω |u exact -u| 2 dΩ/ Ω |u exact | 2 dΩ. We compare the pure FEM approach (which uses only a polynomial description), the pure VTCR approach (which uses only a wave description) and the WTDG approach (which uses at the same time the polynomial and the wave descriptions). Each time, the comparisons are done with a WTDG which is based on the same discretization than the pure FEM or the pure VTCR. The convergence curve is represented in Fig. 3.

As one can see, the WTDG presents a better behavior than the pure FEM or the pure VTCR. The pure FEM suffers from a lack of accuracy when the frequency becomes to be too high. The pure VTCR is not so efficient in the low frequency domain. This shows the benefits of using the WTDG method for finding the solution for low and medium frequency problems with the same descriptions, at the same time.

The convergence of the WTDG method depends on discretization strategies of both polynomial and wave shape functions. A study has been made to see how the discretization strategy of wave shape functions affects the performance of the WTDG method. With the model of the same computational domain and the same boundary condition defined in this section, we take k = 25 m -1 . Seven different discretization strategies of wave shape functions have been used to draw the convergence curves of the WTDG method as Fig. 4 shows. For each discretization strategy of wave functions, increasing of degree of freedoms of polynomial functions leads to convergence of the WTDG method. As one can see, for a fixed number of waves, the curves converge when the number of polynomial shape functions is increased. It can be seen that when using the same degree of freedoms of polynomial functions, enrichment of wave functions used in WTDG leads to more precise result. For the VTCR method, 32 wave shape functions are sufficient to have a precise simulation result in this problem. So here the convergence curve of FEM with 32 VTCR waves locates always in a zone of good precision.

Nonhomogeneous Helmholtz problem with two scales in the solution

We now consider the Helmholtz problem defined on Ω = [-0.5 m; 0.5 m] × [-0.5 m; 0.5 m], with k = 100 m -1 , θ 1 = 30

• , θ 2 = 82 • , r d = (k 2 -k 2 e ) exp(ik e (1 + iζ) -1/2 (cos θ 2 • x + sin θ 2 • y))
, k e = 10 m -1 and ζ = 0.0001. We prescribe the boundary conditions such that the exact solution is u exact = exp(ik(1 + iζ) -1/2 (cos θ 1 • x + sin θ 1 • y)) + exp(ik e (1 + iζ) -1/2 (cos θ 2 • x + sin θ 2 • y)). This example is again interesting, because it corresponds to a nonhomogeneous Helmholtz problem with two scales in the solution (slow varying scale with k e and fast varying scale with k). The exact solution u exact is represented in Fig. 5.

We here consider the WTDG method, used with Dirichlet boundary condition with α = 0.0001. We mix the polynomial and wave shape functions. As we expect here that the polynomial shape functions help us to find the slow varying scale, and that the waves recover the fast varying scale, we consider a regular squared mesh with a polynomial discretization of degree 2, and 2 waves, propagating in the 30 • and 210 • directions. The FEM descriptions uses 10 elements per wave length to simulate the slow varying scale. The exact solution of fast varying scale is taken directly as shape function of VTCR part. There is no need to add more shape function to simulate the fast varying scale. With such a choice, the solution given by the wave description, denoted u VTCR and the polynomial description, denoted u FEM is depicted in Fig. 5.

The comparison between the exact and the WTDG solution is represented in Fig. 6. According to the definition of the error in Sec. 3.1, the error is here 4.48 × 10 -5 . As one can see, the WTDG method is able to recover the exact approximation with a very good accuracy and adopts the scale of the approximation to the physics of the problem.

The WTDG method applied with different types of approximations

We finally consider the computational domain defined in Fig. 7. This L shape domain is filled with a fluid with k = 30 m -1 and ζ = 0.0001. The boundary conditions are such elements of degree 1 per wave length for the FEM description, and δ = 14 for the VTCR description. The WTDG uses δ = 17 for the wave approximations and 6 elements of degree 1 in the subdomain. Again, Dirichlet boundary conditions with α = 0.0001 allow us to prescribe u exact on the boundaries of Ω.

The exact solution and the WTDG solution are depicted in Fig. 8. As one can see, the solutions are very closed. This is because the formulation of weak Trefftz Galerkin method allows the couple using FEM and VTCR and their hybrid. According to the definition of the error in Sec. 3.1, the error is here 2.187 × 10 -2 . It can be deduced from this example that all combinations of methods such as pure FEM, pure VTCR, hybrid of FEM and VTCR can be integrated together in one complex geometry problem. In each subdomain the concrete method can be chosen depend on specific requirement of engineering problems.

Conclusion

This paper proposed an hybrid use of FEM and VTCR thanks to the Weak Trefftz Discontinuous Galerkin method. It is illustrated on the Helmholtz problem. The WTDG method allows one to use a combination of traditional FEM shape functions and VTCR shape functions. It is based on a variational formulation which is equivalent to the reference problem. All the conditions such as the governing equation, the transmission continuity and the boundary conditions are included in the formulation. No a priori constraint is needed for the definition of the shape functions that are used. As a consequence, any shape function can be used, with no difficulty. It gives then to the WTDG method a great flexibility, as the used can select polynomials or wave shape functions (or a combination of them) very easily in the approximated space, with no restriction. It is successfully illustrated on different examples of different complexity, homogeneous or not, with sometimes two scales in the solution. An extension to vibro-acoustic problems would present no particular difficulty and will be addressed in future works.

Fig. 1 .

 1 Fig. 1. Left: Definition of the computational domain Ω. Right: Definition of the subdomains of Ω.

Fig. 2 .

 2 Fig. 2. Example considered in Sec. 3.1. Left: Definition of domain. Middle: VTCR wave directions discretization. Right: FEM mesh refinement.

Fig. 3 .

 3 Fig.3. The convergence curves for the example of Sec. 3.1. The FEM curve corresponds to the solution obtained with a pure FEM discretization explained in Sec. 3.1. The VTCR curve corresponds to the solution obtained with a pure VTCR discretization explained in Sec. 3.1. The WTDG curve corresponds to the solution obtained with an enrichment of the FEM shape functions with waves, according to the WTDG strategy explained in Sec. 2.4.

Fig. 4 .

 4 Fig.4. The convergence curves for the example of Sec. 3.1. For each convergence curve, a fixed number of wave directions of VTCR part is chosen in WTDG strategy. The degree of freedoms of FEM part is varied in order to attain the convergence.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Example considered in Sec. 3.2. Up left: Definition of the computational domain. Up right: Exact solution u exact . Down left: Representation of the fast varying scale result simulated by VTCR part u VTCR . Down right: Representation of the slow varying scale result simulated by FEM part u FEM .

that the exact solution is

• cos( π 3 )y). On this example, we use different kinds of approximations: even a pure FEM approximation, or a pure VTCR approximation, or a mix between polynomial and wave shape functions (see Fig. 7). The formulation of Discontinuous Galerkin method allows this possibility. This example differs from the one given in Ref. 25 by the fact that here, in one subdomain of Ω, we mix the approximations. In order to be sure to have a good approximation, we select 20