# Extended variational theory of complex rays in heterogeneous Helmholtz problem

Hao Li, Pierre Ladevèze, Hervé Riou

► To cite this version:

Hao Li, Pierre Ladevèze, Hervé Riou. Extended variational theory of complex rays in heterogeneous Helmholtz problem. Computational Mechanics, 2017, 59 (6), pp.909-918. 10.1007/s00466-017-1385-4. http://dx.doi.org/10.1007/s00466-017-1385-4

HAL Id: hal-01647645

https://hal.science/hal-01647645

Submitted on 17 Dec 2019

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

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Extended variational theory of complex rays in heterogeneous Helmholtz problem

Hao Li<fnref n="1" />, ∗ · Pierre Ladeveze<fnref n="1" /> · Hervé Riou<fnref n="1" />

**Abstract** 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.

**Keywords** VTCR · Trefftz methods · Discontinuous Galerkin methods · Heterogeneous Helmholtz problem

## 1 Introduction

With the development of numerical simulation, many partial differential equations PDEs in engineering problem can be 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-

cient calculations for problems where the medium and the high frequencies are considered.

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 series 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

<figure><img src="image_2.png" /></figure>

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

Pierre Ladeveze
pierre.ladeveze@lmt.ens-cachan.fr

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

<footnote n="1">¹ LMT, ENS Cachan, CNRS, Université Paris-Saclay,
61 Avenue de Président Wilson, 94235 Cachan, France</footnote>
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. Then its properties can be seen in the academic examples in Sect. 3. And a semi-unbounded harbor agitation example solved by extended VTCR is presented in this section. Conclusions are drawn in Sect. 4.

## 2 Extended VTCR for slowly varying wave number Helmholtz problem

### 2.1 The reference problem

A 2-D Helmholtz problem is taken as reference problem. Let $\Omega$ be the computational domain and $\partial\Omega = \partial\Omega_1 \cup \partial\Omega_2$ be the boundary. Without losing generality, Dirichlet and Neumann conditions are prescribed on $\partial\Omega_1$ and $\partial\Omega_2$ to illustrate this extended VTCR method. Treatment of other different boundary conditions can be seen in [15]. The following problem is considered: find $u \in H^1(\Omega)$ such that

$$ 
\begin{array}{l@{\qquad}l}
(1 - i\eta)\Delta u + k^2 u = 0 & \text{over } \Omega \\
u = u_d & \text{over } \partial\Omega_1 \\
(1 - i\eta)\partial_n u = g_d & \text{over } \partial\Omega_2
\end{array}
\qquad (1) $$

where $\partial_n u = \text{grad} u \cdot n$ and $n$ is the outward normal. $u$ is the physical variable studied such as the pressure in acoustics. $\eta$ 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 $\Omega$. $u_d$ and $g_d$ are the prescribed Dirichlet and Neumann data.

### 2.2 The variational formulation of the reference problem

Let $\Omega$ be partitioned into $N$ non overlapping subdomains $\Omega = \cup_{E=1}^N \Omega_E$. Denoting $\partial\Omega_E$ the boundary of $\Omega_E$, we define $\Gamma_{EE} = \partial\Omega_E \cap \partial\Omega$ and $\Gamma_{EE'} = \partial\Omega_E \cap \partial\Omega_{E'}$. The VTCR approach consists in searching solution $u$ in functional space $\mathcal{U}$ such that

$$ 
\begin{aligned}
\mathcal{U} &= \{u \mid u|_{\Omega_E} \in \mathcal{U}_E\} \\
\mathcal{U}_E &= \{u_E \mid u_E \in \mathcal{Y}_E \subset H^1(\Omega_E) | (1 - i\eta)\Delta u_E \\
&\qquad + k^2 u_E = 0\}
\end{aligned}
\qquad (2) $$

Denoting

$$ 
\begin{aligned}
\{u\}_{EE'} &= (u_E + u_{E'})|_{\Gamma_{EE'}} \\
[u]_{EE'} &= (u_E - u_{E'})|_{\Gamma_{EE'}} \\
\mathbf{q}_u &= (1 - i\eta)\text{grad} u
\end{aligned}
\qquad (3) $$

the variational formulation of (1) can be written as: find $u \in \mathcal{U}$ such that

$$ 
\begin{aligned}
&\mathrm{Re} \left( -ik \left( \sum_{E,E' \in \mathcal{E}} \int_{\Gamma_{EE'}} \left( \frac{1}{2} \{\mathbf{q}_u \cdot \mathbf{n}\}_{EE'} \{\tilde{\mathbf{v}}\}_{EE'} - \frac{1}{2} [\tilde{\mathbf{q}}_v \cdot \mathbf{n}]_{EE'} [u]_{EE'} \right) dS \right. \right. \\
&\qquad \left. - \sum_{E \in \mathcal{E}} \int_{\Gamma_{EE} \cap \partial_1 \Omega} \tilde{\mathbf{q}}_v \cdot \mathbf{n} (u - u_d) dS \right. \\
&\qquad \left. + \sum_{E \in \mathcal{E}} \int_{\Gamma_{EE} \cap \partial_2 \Omega} (\mathbf{q}_u \cdot \mathbf{n} - g_d) \tilde{v} dS \right) \\
&= 0 \quad \forall v \in \mathcal{U}
\end{aligned}
\qquad (4) $$

where $\tilde{\mathcal{D}}$ represents the conjugation of $\mathcal{D}$. The existence and uniqueness of solution in this kind of variational formulation have been proved in [15].

### 2.3 The new shape function in VTCR

As mentioned in Sect. 1, in this paper the wave number is in the form that $k^2 = \alpha x + \beta y + \gamma$, where $\alpha, \beta, \gamma$ are constant parameters. For simplicity, in this section we denote that $k_{\uparrow}^2 = k^2/(1 - i\eta) = \alpha_{\uparrow}x + \beta_{\uparrow}y + \gamma_{\uparrow}$, where $\alpha_{\uparrow} = \alpha/(1 - i\eta)$, $\beta_{\uparrow} = \beta/(1 - i\eta)$, $\gamma_{\uparrow} = \gamma/(1 - i\eta)$ 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:

$$ \left( \frac{F''}{F} + \alpha_{\uparrow}x + \gamma_{\uparrow} \right) = - \left( \frac{G''}{G} + \beta_{\uparrow}y \right) \equiv \delta \qquad (5) $$

where $\delta$ is a free constant parameter. The analytic solutions of (5) are:

$$ F(x) = \begin{cases} C_1 \mathrm{Ai} \left( \frac{-\alpha_{\uparrow}x - \gamma_{\uparrow} + \delta}{\alpha_{\uparrow}^{2/3}} \right) + C_2 \mathrm{Bi} \left( \frac{-\alpha_{\uparrow}x - \gamma_{\uparrow} + \delta}{\alpha_{\uparrow}^{2/3}} \right) & |\alpha_{\uparrow}| \neq 0 \\ C_1 \cos \left( \sqrt{\gamma_{\uparrow} - \delta x} \right) + C_2 \sin \left( \sqrt{\gamma_{\uparrow} - \delta x} \right) & |\alpha_{\uparrow}| = 0 \end{cases} \qquad (6) $$

$$ G(y) = \begin{cases} D_1 \mathrm{Ai} \left( \frac{-\beta_{\uparrow}y - \delta}{\beta_{\uparrow}^{2/3}} \right) + D_2 \mathrm{Bi} \left( \frac{-\beta_{\uparrow}y - \delta}{\beta_{\uparrow}^{2/3}} \right) & |\beta_{\uparrow}| \neq 0 \\ D_1 \cos \left( \sqrt{\delta y} \right) + D_2 \sin \left( \sqrt{\delta y} \right) & |\beta_{\uparrow}| = 0 \end{cases} \qquad (7) $$
<figure><img src="image_3.png" /><figcaption>Fig. 1 Behaviors of Airy functions</figcaption></figure>

where Ai and Bi are Airy functions [16]. C₁, C₂, D₁, D₂ 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:

$$
\left\{
\begin{aligned}
Bi(-z) &\sim \frac{\cos(\frac{2}{3}z^{\frac{3}{2}} + \frac{\pi}{4})}{\sqrt{\pi z^{\frac{1}{4}}}} & |\arg(z)| < 2\pi/3 \\
Ai(-z) &\sim \frac{\sin(\frac{2}{3}z^{\frac{3}{2}} + \frac{\pi}{4})}{\sqrt{\pi z^{\frac{1}{4}}}} & |\arg(z)| < 2\pi/3
\end{aligned}
\right.
\qquad (8)
$$

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

$$
k_m^2 = \inf_{x \in \Omega} \inf_{y \in \Omega} \alpha x + \beta y + \gamma = \alpha x_m + \beta y_m + \gamma \quad (9)
$$

where $k_m^2$ represents the minimum value of $k^2$ on $\Omega$ and $(x_m, y_m)$ is the coordinate which enables $k^2$ to take its minimum value $k_m^2$. Denoting $\mathbf{P} = [P_1, P_2] = [\cos(\theta), \sin(\theta)]$, where $\theta$ represents an angle parameter ranging from 0 to $2\pi$. By such definition, $k^2$ can be expressed in form that:

$$
\begin{align}
k^2 &= k_m^2 + \alpha(x - x_m) + \beta(y - y_m) \tag{10} \\
    &= k_m^2 P_1^2 + k_m^2 P_2^2 + \alpha(x - x_m) + \beta(y - y_m) \nonumber
\end{align}
$$

As the similar procedure to get (6) and (7), functions F
and G can be composed by:

$$
F(\tilde{x}) = \mathrm{Bi}(-\tilde{x}) + i * \mathrm{Ai}(-\tilde{x}) \quad (11)
$$

$$
G(\tilde{y}) = \mathrm{Bi}(-\tilde{y}) + i * \mathrm{Ai}(-\tilde{y}) \quad (12)
$$

where $\tilde{x}$ and $\tilde{y}$ are defined as follows:

$$
\tilde{x} = \frac{k_m^2 * P_1^2 + \alpha(x - x_m)}{\alpha^{2/3}(1 - i\eta)^{1/3}} = \frac{k_1^2}{\alpha^{2/3}(1 - i\eta)^{1/3}} \quad (13)
$$

$$
\tilde{y} = \frac{k_m^2 * P_2^2 + \beta(y - y_m)}{\beta^{2/3}(1 - i\eta)^{1/3}} = \frac{k_2^2}{\beta^{2/3}(1 - i\eta)^{1/3}} \quad (14)
$$

By such a way, − ˜ x and − ˜ y always locate in [−∞, 0] on
the domain Ω. The new shape function ψ ( x, y, **P** ) is built
as:

$$
\psi(x, y, \mathbf{P}) = F(\tilde{x}) * G(\tilde{y}) \quad (15)
$$

Asymptotically, when $\alpha$ tends to 0

$$
F(\tilde{x}) \rightarrow \cos(k_1 \cdot x) + i \cdot \sin(k_1 \cdot x) \quad (16)
$$

Asymptotically, when $\beta$ tends to 0

$$
G(\tilde{y}) \rightarrow \cos(k_2 \cdot y) + i \cdot \sin(k_2 \cdot y) \quad (17)
$$

It can be observed that $\psi(x, y, \mathbf{P})$ function is the general solution of Helmholtz equation in (1). Especially when $\alpha = 0$ and $\beta = 0$, $\psi(x, y, \mathbf{P})$ function becomes plane wave function. The angle parameter $\theta$ in $\mathbf{P}$ describes the propagation direction of plane wave. Analogous to plane wave case, when $\alpha \neq 0$ and $\beta \neq 0$, $\psi$ function still represents a wave propagates in 2-D plane. $\mathbf{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.

**2.4 Implementation of the extended VTCTR**

In order to implement VTCTR method, it is required to take
into account a finite number of propagation directions of Airy
wave.

$$
\mathcal{U} = \left\{
\begin{array}{l}
u \in L^2(\Omega) : u(x, y)|_{\Omega_E} = \sum_{m_E=1}^{M_E} A_{m_E} \psi(x, y, \mathbf{P}_{m_E}), \\
A_{m_E} \in \mathbb{C}, E = 1, \dots, N
\end{array}
\right\} \qquad (18)
$$

where $M_E$ is the number of Airy wave directions selected in
the subdomain $\Omega_E$ and $A_{m_j}$ is the amplitude of wave. Then
VTCTR leads to resolve linear system of equations:

$$
KA = F \tag{19}
$$

$K$ corresponds to the discretization of the bilinear form of weak formulation. Inside $K$ there are $N^2$ partitioning of blocks $K_{EE'}$. When $\Gamma_{EE'} \neq \emptyset$, the blocks corresponding to $K_{EE'}$ are non zero fully populated. Otherwise $K_{EE'}$ are zero blocks. The vector $A$ corresponds to the amplitudes of waves, which is the degree of freedom in VTCTR. $F$ is the linear form of weak formulation and corresponds to the loading.
<figure><img src="image_4.png" /><figcaption>Fig. 2 Example of Airy wave and plane wave. Left Airy wave with $\eta = 0.001$, $\alpha = 300 \text{ m}^{-3}$, $\beta = 300 \text{ m}^{-3}$, $\gamma = 600 \text{ m}^{-2}$, $\mathbf{P} = [\cos(\pi/6), \sin(\pi/6)]$. Right plane wave with $\eta = 0.001$, $\alpha = 0 \text{ m}^{-3}$, $\beta = 0 \text{ m}^{-3}$, $\gamma = 600 \text{ m}^{-2}$, $\mathbf{P} = [\cos(\pi/6), \sin(\pi/6)]$</figcaption></figure>

## 3 Numerical example

### 3.1 Study of the performance of the extended VTCTR on medium frequency heterogeneous Helmholtz problem

A simple geometry of square [0 m; 1 m] x [0 m; 1 m] is considered for domain $\Omega$. In this domain, $\eta = 0.01$, $\alpha = 150 \text{ m}^{-3}$, $\beta = 150 \text{ m}^{-3}$, $\gamma = 1000 \text{ m}^{-2}$. Boundary conditions on $\partial\Omega$ are Dirichlet type such that $u_d = \sum_{j=1}^3 \psi(x, y, \mathbf{P}_j)$, where $\psi(x, y, \mathbf{P}_j)$ is the Airy wave solution of heterogeneous Helmholtz equation in domain $\Omega$. $\theta_1 = 10^\circ$, $\theta_2 = 55^\circ$, $\theta_3 = 70^\circ$ correspond to propagation angle in $\mathbf{P}_1, \mathbf{P}_2, \mathbf{P}_3$ respectively. Such kind of geometry and boundary conditions allow to calculate the relative error of extended VTCTR 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:

$$ \frac{|u_{exact} - u|}{|u_{exact}|} = \sqrt{\int_{\Omega} |u_{exact} - u|^2 d\Omega / \int_{\Omega} |u_{exact}|^2 d\Omega} $$

As one can see from Fig. 4, the convergence curves of this extended VTCTR method in heterogeneous problem behaves in the same way as the convergence curves of VTCTR 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 VTCTR with plane waves is that $N_e = \tau k R_e / (2\pi)$ where $N_e$ is the number of directions of waves, $\tau$ 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 $\tau = 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 VTCTR convergence property. For VTCTR there exists two
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Example considered in Sect. 3.1. From left to right—*First* definition of domain. *Second* 1 subdomain discretisation. *Third* 4 subdomains discretisation. *Fourth* 9 subdomains discretisation</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4** The convergence curves for the example of Sect. 3.1. The three convergence curves of extended VTCTR calculated with Airy wave correspond to the three discretization strategies shown in Fig. 3</figcaption></figure>

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 VTCTR with only one computational domain converges the fastest.

## 3.2 Study of the extended VTCTR 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

<figure><img src="image_7.png" /><figcaption>**Fig. 5** Example considered in Sect. 3.2. Top view of Harbor. $\theta_0^+$ represents the direction of incident wave</figcaption></figure>

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 $\Gamma_R$:

$$ (1 - i\eta)\partial_n u = 0 \quad \text{over } \Gamma_R \qquad (20) $$

$u_0^+$ represents incoming wave from far away onto the harbor. It can be expressed as $u_0^+ = A_0^+ \exp ik_0^+ (\cos\theta_0^+ x + \sin\theta_0^+ y)$, where $A_0^+$ is the amplitude of wave and $\theta_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
<figure><img src="image_8.png" /><figcaption>Fig. 6 Example considered in Sect. 3.2. Side view of Harbor. Variable h represents depth of water from sea surface to the bottom. The depth h increases when it points from harbor inside to harbor outside</figcaption></figure>

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:

$$h = \frac{1}{a + by} \qquad (21)$$

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:

$$\omega^2 = kg\tanh(kh) \qquad (22)$$

where g = 9.81 m/s² 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:

$$\tanh(kh) \approx kh \qquad (23)$$

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:

$$k^2 = g^{-1}\omega^2(a + by) \qquad (24)$$

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 semi-unbounded problem. In physics these waves need to satisfy Sommerfeld radiation condition. In our 2-D model it is represented by:

$$\lim_{r \to +\infty} \sqrt{r} \left( \frac{\partial u(r)}{r} - iku(r) \right) = 0 \qquad (25)$$

where r is the radial direction in polar coordinate.

<figure><img src="image_9.png" /><figcaption>Fig. 7 Example considered in Sect. 3.2. First step for seeking analytic solution outside the harbor</figcaption></figure>

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\theta_0^+ x + \sin\theta_0^+ y)$, it can be obtained that $u_a = A_0^+ \exp ik_0^+ (\cos\theta_a x + \sin\theta_a y)$, where $\theta_a = \pi - \theta_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 $\Gamma_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]:
$$U_s = \sum_{n=0}^{\infty} (A_n \sin(n\theta) + B_n \cos(n\theta)) H_n^{(1)}(kr) \quad (26)$$

where $H_n^{(1)}(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_n^{(1)}(kr)$ is singular on the origin. It should be

<figure><img src="image_10.png" /><figcaption>Fig. 8 Example considered in Sect. 3.2. Second step for seeking analytic solution outside the harbor</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9 Half plane problem with boundary $\Gamma_O$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10 Example considered in Sect. 3.2. Left The first strategy: Domain inside the harbor divided into one computational subdomain. Right The second strategy: Domain inside harbor divided into four computational subdomains</figcaption></figure>

noticed that if (26) can be modified to satisfy the boundary condition on $\Gamma_O$, it will be the $u_b$ we search. Therefore the problem can be simplified and abstracted into Fig. 9. $\Gamma_O$ are straight and infinitely long boundaries with the same boundary conditions as (20). $d$ is set to be an 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:

$$u_b = \sum_{n=0}^{\infty} B_n \cos(n\theta) H_n^{(1)}(kr) \quad (27)$$

It can be verified that (27) satisfies boundary conditions on $\Gamma_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 $\Omega_1$ and $\Omega_2$. The subdomain $\Omega_2$ is a semicircular domain, whose center locates at the origin point. The subdomain $\Omega_1$ ranges from the boundary of $\Omega_2$ to infinity. On this domain the analytic solution presented before is used. Computational domain $\Omega_2$ is created to separate origin point from $\Omega_1$. Since $k$ is considered as constant value of the region outside the harbor, plane wave function is used as shape function on subdomain $\Omega_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,

<figure><img src="image_13.png" /></figure>
one could see VTCR leads to a unique and convergent simulation result in this problem.

When the subdivision of computational domain is done, one needs to choose shape functions used on each subdomain. As mentioned before, $u$ on domain $\Omega_1$ contains $u_0^+$, $u_a$ and $u_b$. This relation can be represented by $u|_{\Omega_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 $\mathcal{U}_{\Omega_1}^b$ is defined as:

$$ \mathcal{U}_{\Omega_1}^b = \left\{ u_b \in L^2(\Omega_1): u_b(x, y) = \sum_{n=0}^{N_1} A_{1n} \cos n\theta H_n^{(1)}(kr), \right. \\ \left. A_{1n} \in \mathbb{C}, n = 0, \dots, N_1 \right\} \quad (28) $$

where $A_{1n}$ is the unknown degree of freedom. $N_1$ is the number of degree of freedom on $\Omega_1$. Working space of $\Omega_2$ is defined as follows:

$$ \mathcal{U}_{\Omega_2} = \left\{ u \in L^2(\Omega_2): u(x, y) = \sum_{n=0}^{N_2} A_{2n} \exp ik(\cos\theta_n x + \sin\theta_n y), \right. \\ \left. A_{2n} \in \mathbb{C}, n = 1, \dots, N_2 \right\} \quad (29) $$

where $A_{2n}$ is the unknown amplitude of plane wave. $N_2$ is the number of degree of freedom on $\Omega_2$.

On the computational domain of inside harbor, $\psi(x, y, P_j)$ function forms the working space in the form of (18). Here $\omega = 0.5$ rad/s, $a = 4.8 \times 10^{-2}$ m$^{-1}$, $b = 4.8 \times 10^{-5}$ m$^{-2}$, $\eta = 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:

$$ k^2 = 1.2 \times 10^{-3} - 1.2 \times 10^{-6}y \quad (30) $$

Inside the harbor $k^2 \in [1.2 \times 10^{-3}$ m$^{-2}$, $3.0 \times 10^{-3}$ m$^{-2}$] and $\lambda \in [104.72$ m, $181.38$ m]. The shallow water approximation (23) is approved to be valid since $\lambda \gg h$.

Let the amplitude of incoming wave corresponds to $A_0^+ = 2$ m and the angle of incoming wave corresponds to $\theta_0^+ = 45^\circ$. 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 $\Omega_1$, $\Omega_2$ and $\Omega_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 $\Omega_1$,

<figure><img src="image_14.png" /><figcaption>Fig. 11 Example considered in Sect. 3.2. Up numerical result calculated by the first strategy with $\theta_0^+ = 45^\circ$. Down Numerical result calculated by the second strategy with $\theta_0^+ = 45^\circ$. Results of semi-unbounded domain $\Omega_1$ are shown in a truncated part with $r \in [1000$ m, $2000$ m] in polar coordinate</figcaption></figure>

$\Omega_2$, $\Omega_3$, $\Omega_4$, $\Omega_5$, $\Omega_6$ respectively. The geometrical heuristic criterion of convergence for $\Omega_1$ and $\Omega_2$ can refer to [14]. For the Airy wave, it is the same criterion presented in Sect. 3.1. Since $\Omega_1$ is the semi-unbounded domain, here the numerical result only shows a truncated part with $r \in [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
<figure><img src="image_15.png" /><figcaption>**Fig. 12** Example considered in Sect. 3.2. *Up* numerical result inside the harbor calculated by the first strategy with $\theta_0^+ = 45^\circ$. *Down* Numerical result inside the harbor calculated by the second strategy with $\theta_0^+ = 45^\circ$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 13** Example considered in Sect. 3.2. *Up* numerical result calculated by the first strategy with $\theta_0^+ = 35^\circ$. *Down* Numerical result calculated by the first strategy with $\theta_0^+ = 65^\circ$. Results of semi-unbounded domain $\Omega_1$ are shown in a truncated part with $r \in [1000 \text{ m, } 2000 \text{ m}]$ in polar coordinate</figcaption></figure>

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 $\theta_0^+ = 35^\circ$ and $\theta_0^+ = 65^\circ$ (see Fig. 13).

## 4 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 conver-
gence 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.

References

1. Strouboulis T, Hidajat R (2006) Partition of unity method for Helmholtz equation: q-convergence for plane-wave and wave-band local bases. Appl Math 51:181–204

2. Cessenat O, Despres B (1998) Application of an ultra weak variational formulation of elliptic PDEs to the two-dimensional Helmholtz problem. SIAM J Numer Anal 35:255-299

3. Gabard G, Gamallo P, Huttunen T (2011) A comparison of wave-based discontinuous Galerkin, ultra-weak and least-square methods for wave problems. Int J Numer Methods Eng 85:380-402

4. Gittelson CJ, Hiptmair R, Perugia I (2009) Plane wave discontinuous Galerkin methods: analysis of the h-version. ESAIM Math Model Numer Anal 43:297-331

5. Barnett A, Betcke T (2008) Stability and convergence of the method of fundamental solutions for Helmholtz problems on analytic domains. J Comput Phys 227:7003-7026

6. Farhat C, Harari I, Franca L (2001) The discontinuous enrichment method. Comput Methods Appl Mech Eng 190:6455-6479

7. Desmet W, Sas P, Vandepitte D (2001) An indirect Trefftz method for the steady-state dynamic analysis of coupled vibro-acoustic systems. Comput Assist Mech Eng Sci 8:271-288

8. Tezaur R, Kalashnikova I, Farhat C (2014) The discontinuous enrichment method for medium-frequency Helmholtz problems with a spatially variable wavenumber. Comput Methods Appl Mech Eng 268:126-140

9. Ladevèze P (1996) A new computational approach for structure vibrations in the medium frequency range. C R Acad Sci Paris 332:849–856

10. Rouch P, Ladevèze P (2003) The variational theory of complex rays: a predictive tool for medium-frequency vibrations. Comput Methods Appl Mech Eng 192:3301–3315

11. Ladevèze P, Blanc L, Rouch P, Blanzé C (2003) A multiscale computational method for medium-frequency vibrations of assemblies of heterogeneous plates. Comput Struct 81:1267–1276

12. Riou H, Ladevèze P, Rouch P (2004) Extension of the variational theory of complex rays to shells for medium-frequency vibrations. Sound Vib 272:341–360

13. Riou H, Ladevèze P, Sourcis B (2008) The multiscale VTCR approach applied to acoustics problems. Comput Acoust 16:487–505

14. Kovalevsky L, Riou H, Ladevèze P (2013) On the use of the variational theory of complex rays for the analysis of 2-D exterior Helmholtz problem in an unbounded domain. J Wave Motion 50:428–436

15. Ladevèze P, Riou H (2014) On Trefftz and weak Trefftz discontinuous Galerkin approaches for medium-frequency acoustics. Comput Methods Appl Mech Eng 278:729-743

16. Polyanin AD, Zaitsev VF (2002) Handbook of exact solutions for ordinary differential equations. CRC Press

17. Modesto D, Zlotnik S, Huerta A (2015) Proper generalized decomposition for parameterized Helmholtz problems in heterogeneous and unbounded domains: application to harbor agitation. Comput Methods Appl Mech Eng 295:127–149

18. Berenger JP (1994) A perfectly matched layer for the absorption of electromagnetic waves. J Comput Phys 114:185-200

19. Givoli D (2004) High-order local non-reflecting boundary conditions: a review. J Wave Motion 39:319–326

20. Bayliss A, Turkel E (1980) Radiation boundary conditions for wave like equations. Commun Pure Appl Math 33:707-725

21. Antoine X, Barucq H, Bendali A (1999) Bayliss Turkel-like radiation conditions on surfaces of arbitrary shape. J Math Anal Appl 229:184–211

22. Givoli D (1999) Recent advances in the DtN FE method. Arch. Comput. Methods Eng. 6:71–116

23. Herrera I (1984) Boundary methods: an algebraic theory. Pitman Advanced Publishing Program, Boston