Extended variational theory of complex rays in heterogeneous Helmholtz problem

In the past years, a numerical technique method called Variational Theory of Complex Rays (VTCR) has been developed for vibration problems in medium frequency. It is a Trefftz Discontinuous Galerkin method which uses plane wave functions as shape functions. However this method is only well developed in homogeneous case. In this paper, VTCR is extended to the heterogeneous Helmholtz problem by creating a new base of shape functions. Numerical examples give a scope of the performances of such an extension of VTCR.


Introduction
With the development of numerical simulation, many partial differential equations PDEs in engineering problem can by solved by computer-aided methods. Among these PDEs based problems, Helmholtz equation models the time harmonics wave equation, which appears in a large range of applications such as acoustics, electromagnetism and quantum mechanics. Relying upon the wave number k, Finite Element Method (FEM) could only be used in small wave number cases since when the wave number increases, the meshes need drastically to be refined. This leads to ineffi- In order to have a good resolution for medium frequency problems, a great number of Trefftz method based approaches are developed, which makes use of a priori knowledge of exact solution to build shape functions. These methods are, for example, the partition of unity method [1], the ultra weak variational method (UWVF) [2], the least square method [3], the plane wave discontinuous Galerkin methods [4], the method of fundamental solutions [5], the discontinuous enrichment method (DEM) [6], the wave based method (WBM) [7]. In the slowly varying wave number Helmholtz problem, the UWVF method makes use of plane wave functions to approximate the solution of the one dimension (1D) problem. The DEM method develops a series of Airy function based enrichment shape functions to study 2D scattering problem [8].
The Variational Theory of Complex Rays (VTCR), first introduced in [9], belongs to this category of numerical strategies. VTCR uses wave functions as shape functions to get approximated solutions for vibration problems. It has been developed for 3-D plate assemblies in [10], for plates with heterogeneities in [11] and for shells in [12]. Its extensions to acoustics problems can be seen in [13,14]. However, all these cases are limited to homogeneous Helmholtz problems. Facing to slowly varying wave number Helmholtz problem where the square of wave number depends linearly on the coordinates, extended VTCR generates a serie of new shape functions, which are composed by Airy functions and satisfy a priori the dominant equation. Comparing to the enrichment shape functions proposed in [8], these new shape functions are created in a different way in this paper.
Heterogeneous Helmholtz problems exist in problems such as underwater acoustics, wave propagation in geophysics, electromagnetism. Sometimes these problems need to be resolved in unbounded domain or semi-unbounded domain. Additional techniques are required to transform the unbounded domain into bounded computational domain for numerical calculation. In this paper, a semi-unbounded heterogeneous Helmholtz problem is solved by extended VTCR method in a simple way without applying additional techniques to semi-unbounded domain.
As the objectives mentioned above, the paper is organised as follows: an extended VTCR is introduced by adding a new base of shape functions for slowly varying wave number Helmholtz problem in Sect. 2 [15]. The following problem is considered: find u ∈ H 1 (Ω) such that where ∂ n u = gradu · 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 and i is the imaginary unit. Being different from Helmholtz problem resolved by VTCR in previous work, k is not a constant and its value changes with respect to location (x, y) on Ω. u d and g d are the prescribed Dirichlet and Neumann data.

The variational formulation of the reference problem
Let Ω be partitioned into N non overlapping subdomains The VTCR approach consists in searching solution u in functional space U such that Denoting the variational formulation of (1) can be written as: find u ∈ U such that where˜ represents the conjugation of . The existence and uniqueness of solution in this kind of variational formulation have been proved in [15].

The new shape function in VTCR
As mentioned in Sect. 1, in this paper the wave number is in the form that k 2 = αx + βy + γ , where α, β, γ are constant parameters. For simplicity, in this section we denote that k 2 respectively. In VTCR method, exact solution needs to be known a priori to serve as shape function. Therefore exact solution of heterogeneous Helmholtz equation in (1) is required to be found. Separation of variable is considered here. By introducing u(x, y) = F(x)G(y) into (1), it can be obtained that: where δ is a free constant parameter. The analytic solutions of (5) are: where Ai and Bi are Airy functions [16]. C 1 , C 2 , D 1 , D 2 are constant coefficients. In the interval [0, +∞] function Ai tends towards 0 and function Bi tends towards infinity (see Fig. 1). Moreover when a variable named −z tends to −∞, the asymptotic expression of function Ai and Bi are: In the interval [0, +∞], Bi goes to infinity and it has no physical meaning. Thus choosing Airy function in the [−∞, 0] to compose new shape function is taken into account. By this way, a new shape function ψ(x, y, P) is defined as follow procedures: Defining where k 2 m represents the minimum value of k 2 on Ω and (x m , y m ) is the coordinate which enables k 2 to take its minimum value k 2 m . Denoting P = [P 1 , P 2 ] = [cos(θ ), sin(θ )], where θ represents an angle parameter ranging from 0 to 2π . By such definition, k 2 can be expressed in form that: As the similar procedure to get (6) and (7), functions F and G can be composed by: wherex andỹ are defined as follows: By such a way, −x and −ỹ always locate in [−∞, 0] on the domain Ω. The new shape function ψ(x, y, P) is built as: Asymptotically, when α tends to 0 Asymptotically, when β tends to 0 It can be observed that ψ(x, y, P) function is the general solution of Helmholtz equation in (1). Especially when α = 0 and β = 0, ψ(x, y, P) function becomes plane wave function. The angle parameter θ in P describes the propagation direction of plane wave. Analogous to plane wave case, when α = 0 and β = 0, ψ function still represents a wave propagates in 2-D plane. P decides its propagation direction. In order to be distinct from plane wave, this wave is named Airy wave. An example of Airy wave and plane wave can be seen in Fig. 2.

Implementation of the extended VTCR
In order to implement VTCR method, it is required to take into account a finite number of propagation directions of Airy wave.
where M E is the number of Airy wave directions selected in the subdomain Ω E and A m j is the amplitude of wave. Then VTCR leads to resolve linear system of equations: K corresponds to the discretization of the bilinear form of weak formulation. Inside K there are N 2 partitioning of blocks K E E . When Γ E E = ∅, the blocks corresponding to K E E are non zero fully populated. Otherwise K E E are zero blocks. The vector A corresponds to the amplitudes of waves, which is the degree of freedom in VTCR. F is the linear form of weak formulation and corresponds to the loading. ψ(x, y, P j ), where ψ(x, y, P j ) is the Airy wave solution of heterogeneous Helmholtz equation in domain Ω. θ 1 = 10 • , θ 2 = 55 • , θ 3 = 70 • correspond to propagation angle in P 1 , P 2 , P 3 respectively. Such kind of geometry and boundary conditions allow to calculate the relative error of extended VTCR method with exact solution. Thus it allows to see the performance of this approach. The definition of the problem and the discretization strategy can be seen on Fig. 3. In order to capture the relative error, the following definition is used in this paper: As one can see from Fig. 4, the convergence curves of this extended VTCR method in heterogeneous problem behaves in the same way as the convergence curves of VTCR applied in k 2 constant case with plane wave [13]. Merely a small amount of degrees of freedom is sufficient to attain the convergence of numerical result, which is under a small relative error. The geometrical heuristic criterion of convergence for VTCR with plane waves is that N e = τ k R e /(2π) where N e is the number of directions of waves, τ a parameter to be chosen, k the wave number and R e is the characteristic radius of domain [14]. Since k is not constant in this example, its maximum value on the domain is used in the heuristic criterion.
Here we take τ = 10. It can be seen that to obtain the result with same precision, refinement of subdomains results in the need of more degrees of freedom. This phenomena relates to VTCR convergence property. For VTCR there exists two  The convergence curves for the example of Sect. 3.1. The three convergence curves of extended VTCR calculated with Airy wave correspond to the three discretization strategies shown in Fig. 3 convergence strategies. The first one named p-convergence is to keep the number of subdomains and to increase the number of wave directions. The second one named h-convergence is to keep the number of wave directions and to increase the number of subdomains. These two strategies all lead to convergent results but p-convergence performs in a far more efficient way. This is the reason why in Fig. 4 VTCR with only one computational domain converges the fastest.

Study of the extended VTCR on semi-unbounded harbor agitation problem
This example corresponds to a study of water agitation of a harbor. The movement of waves is dominated by Helmholtz equation. Incoming wave from far away field gives rise to reflected wave inside the harbor [17]. The water wave length is much smaller than the geometry size of harbor. It is a medium frequency Helmholtz problem since there exist many periods of wave in the harbor. Definition of the harbor is shown in Fig. 5. It is obvious that the agitation of harbor depends on the direction of incoming wave. In later part of this section, one can see three different numerical results calculated with different incoming waves. All boundaries of the harbor are supposed to be reflecting boundaries, which is denoted by Γ R : u + 0 represents incoming wave from far away onto the harbor. It can be expressed as u + 0 = A + 0 exp ik + 0 (cosθ + 0 x+sinθ + 0 y) , where A + 0 is the amplitude of wave and θ + 0 is the angle of wave propagation direction. The origin of coordinate is O, located in the middle point of the harbor entrance. As Fig. 6 shows, the sea bottom of the region outside the harbor varies slowly and the depth of water is considered as constant there. The depth of water inside the harbor decreases when it is closer to the land. Thus it makes the length of wave vary inside the harbor. An assumption is proposed in this example that the depth of water h complies with the following relation: where a, b are constant parameters. This relation could describe the variation of the water depth with respect to y. The relation between wave frequency ω and water depth h follows the non linear dispersion relation: where g = 9.81 m/s 2 is the gravitational acceleration and k is the wave number. In the case h λ, when the depth of water is far more less than the length of wave, there is the following shallow water approximation: This approximation is valid in the underwater field near seashore. The numerical result of this section will further approve the validation of this approximation. Thus it can be obtained that: Incoming waves cause two kinds of reflection, which include the wave reflected by the boundary inside the harbor and the wave reflected by the boundary locating outside the harbor. Part of these reflected waves propagate from the harbor to far away field. This phenomenon leads to a semiunbounded problem. In physics these waves need to satisfy Sommerfeld radiation condition. In our 2-D model it is represented by: where r is the radial direction in polar coordinate. Many methods have been proposed to solve unbounded problem such as perfectly matched layer (PMLs) [17,18], Non-reflecting artificial boundary conditions (NRBC) [19], Bayliss, Gunzburger and Turkel Local non-reflecting boundary conditions (BGT-like ABC) [20,21] and Dirichlet to Neumann non-local operators [22]. PMLs creates an artificial boundary and a layer outside the region of interest in order to absorb the outgoing waves. NRBC, ABC and Dirichlet to Neumann non-local operators introduce a far away artificial boundary which leads to minimize spurious reflections. VTCR method can combine these artificial boundary techniques to solve the semi-unbounded harbor problem without difficulty. But here analytic solution is taken into account to solve the problem. This choice allows us to take great advantage of VTCR method. Since analytic solution verifies Helmholtz equation and Sommerfeld radiation condition, it can be used as shape functions in VTCR. Compared with artificial boundary techniques, this approach leads to a simpler strategy of calculation.
The idea of seeking for analytic solution on the domain outside the harbor can be illustrated by two steps. As Fig. 7 shows, in the first step a relatively simple problem is considered. Without the region inside the harbor, incoming wave u + 0 agitates on a straight boundary which is infinitely long. The boundary condition here is same as (20). The reflected wave is denoted by u a . It is evident that for such a problem, when u + 0 = A + 0 exp ik + 0 (cosθ + 0 x+sinθ + 0 y) , it can be obtained that u a = A + 0 exp ik + 0 (cosθ a x+sinθ a y) , where θ a = π − θ + 0 . For the second step as Fig. 8 shows, it is exactly the original harbor agitation problem in this Section. If u a of the first step is taken as exact solution here, it will create the residual value because the governing equation inside the harbor and boundary conditions are not satisfied. It is logical to add a complementary solution outside the harbor to offset the residual value. In this point of view, the origin O is chosen to develop the expansion of this complementary solution, which is denoted by u b . Here u b is required to satisfy governing equation outside the harbor, where the wave number is constant. Furthermore u b is required to satisfy the boundary condition on Γ O and Sommerfeld radiation condition.
In previous work of VTCR [14], it is shown that for 2-D acoustic domain exterior to a circular boundary surface, the analytic solution of reflected wave U s of scattering problem in polar coordinate is in form of [23]: (A n sin(nθ) + B n cos(nθ)) H (1) n (kr) where H (1) n (kr) is the first type of Hankel function, which represents outgoing waves. A n and B n are constant coefficients. Solution U s satisfies Sommerfeld radiation condition automatically. H (1) n (kr) is singular on the origin. It should be noticed that if (26) can be modified to satisfy the boundary condition on Γ O , it will be the u b we search. Therefore the problem can be simplified and abstracted into Fig. 9. Γ O are straight and infinitely long boundaries with the same boundary conditions as (20). d is set to be a arbitrary distance from singular point O of (26). The main purpose of Fig. 9 is to find the analytic solution u b outside the semicircular domain. u b can be expressed in the form that: It can be verified that (27) satisfies boundary conditions on Γ O . Therefore u b is found. Except on the origin point, the analytic solution on the domain outside the harbor equals to the sum of u + 0 , u a and u b . As mentioned before, our computational strategies are shown as Fig. 10. The domain outside the harbor is divided into two computational subdomains Ω 1 and Ω 2 . The subdomain Ω 2 is a semicircular domain, whose center locates at the origin point. The subdomain Ω 1 ranges from the boundary of Ω 2 to infinity. On this domain the analytic solution presented before is used. Computational domain Ω 2 is created to separate origin point from Ω 1 . Since k is considered as constant value of the region outside the harbor, plane wave function is used as shape function on subdomain Ω 2 .
Inside the harbor two strategies of discretization are chosen as Fig. 10 shows. The first strategy is that the domain inside the harbor is divided into one computational subdomain. The second strategy is that the domain inside the harbor is divided into four computational subdomains. By comparing the numerical results calculated by these two strategies, When the subdivision of computational domain is done, one needs to choose shape functions used on each subdomain. As mentioned before, u on domaine Ω 1 contains u + 0 , u a and u b . This relation can be represented by u| Ω 1 = u + 0 + u a + u b . The unknown value u b can be expanded in the series written as (27). To achieve a discrete version of the VTCR, finite-dimensional space is required. Thus (27) needs to be truncated into finite series. The working space of u b denoted by U b Ω 1 is defined as: where A 1n is the unknown degree of freedom. N 1 is the number of degree of freedom on Ω 1 . Working space of Ω 2 is defined as follows: A 2n exp ik(cosθ n x+sinθ n y) , where A 2n is the unknown amplitude of plane wave. N 2 is the number of degree of freedom on Ω 2 . On the computational domain of inside harbor, ψ(x, y, P j ) function forms the working space in the form of (18). Here ω = 0.5 rad/s, a = 4.8 × 10 −2 m −1 , b = 4.8 × 10 −5 m −2 , η = 0.03 are the chosen as parameters. Therefore the depth of water ranges from −20.83 to −8.33 m, which corresponds to slow variation of water depth near the seashore. The relation between k 2 and y follows (24). Taking into account the parameters, it can be derived that: Inside the harbor k 2 ∈ [1.2 × 10 −3 m −2 , 3.0 × 10 −3 m −2 ] and λ ∈ [104.72 m, 181.38 m]. The shallow water approximation (23) is approved to be valid since λ h. Let the amplitude of incoming wave corresponds to A + 0 = 2 m and the angle of incoming wave corresponds to θ + 0 = 45 • . Following the computational strategies mentioned above, numerical results are shown in Fig. 11. For the first strategy, N 1 = 20, N 2 = 100, N 3 = 160 are the degrees of freedom on Ω 1 , Ω 2 and Ω 3 respectively. For the second strategy, N 1 = 20, N 2 = 100, N 3 = 160, N 4 = 160, N 5 = 160, N 6 = 160 are the degrees of freedom on Ω 1 ,  5 , Ω 6 respectively. The geometrical heuristic criterion of convergence for Ω 1 and Ω 2 can refer to [14]. For the Airy wave, it is the same criterion presented in Sect. 3.1.
Since Ω 1 is the semi-unbounded domain, here the numerical result only shows a truncated part with r ∈ [1000 m, 2000 m] in polar coordinate. Figure 12 is the comparison of results inside the harbor calculated by the first strategy and the second strategy. It shows that the two different computational strategies of extended VTCR method lead to the same result. One point interesting to be noticed is that the performance of the first strategy is slightly better than the second strategy. But the first strategy has less computational subdomains than the second strategy inside the harbor. The reason is explained in Sect. 3.1 that VTCR always performs more efficient when less subdomains are used. It should also be noticed that only 280 degrees of freedom in all are sufficient to solve this medium frequency heterogeneous Helmholtz problem. Such a low requirement of domain subdivisions and of degrees of freedom embodies the advantage of VTCR method. It also can be seen from Fig.  11 that the numerical solution has a good continuity between adjacent subdomains. Combined with the same parameters and with the first computational strategy mentioned before, two other results are calculated by changing the angle of incoming wave to θ + 0 = 35 • and θ + 0 = 65 • (see Fig. 13).

Conclusion
This paper proposes an extended VTCR method, which is able to solve heterogeneous Helmholtz problem. In this extended VTCR method a new type of shape function is created. In the context of Trefftz Discontinuous Galerkin method, this new shape function satisfies governing equation a priori. Thus this extended VTCR method is only required to meet the continuity conditions between subdomains and boundary conditions. All these conditions are included in the variational formulation, which is equivalent to the reference problem. Academic study is done in this paper to show the convergence property of extended VTCR method. And it shows that this approach converges in the same way as classical VTCR method. A study of water agitation of harbor with engineering application background is made. By applying VTCR method, this problem with complex geometry is solved with simple computational domain division and with a small amount of degrees of freedom. It successfully illustrates that VTCR has a significant potential to solve true engineering problem in an efficient and flexible way.
An extension to vibro-acoustic problems would present no particular difficulty and will be addressed in future works.