# On wave based Trefftz and weak Trefftz discontinuous Galerkin approaches for medium-frequency heterogeneous Helmholtz problem

Hao Li, Pierre Ladevèze, Hervé Riou

► To cite this version:

Hao Li, Pierre Ladevèze, Hervé Riou. On wave based Trefftz and weak Trefftz discontinuous Galerkin approaches for medium-frequency heterogeneous Helmholtz problem. Internoise 2017, Aug 2017, Hong Kong, China. pp.1-9. hal-05431397

HAL Id: hal-05431397

https://hal.science/hal-05431397v1

Submitted on 24 Dec 2025

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
<figure><img src="image_2.png" /></figure>

inter.noise

27-30 AUGUST
HONG KONG

2017

# On wave based Trefftz and weak Trefftz discontinuous Galerkin approaches for medium-frequency heterogeneous Helmholtz problem

Hao LI<fnref n="1" />; Pierre LADEVÈZE<fnref n="2" />; Hervé RIOU<fnref n="3" />

¹Ecole Centrale de Lyon; Insa Lyon; Ecole Normale Supérieure Paris-Saclay, France

²Ecole Normale Supérieure Paris-Saclay, France

³Ecole Normale Supérieure Paris-Saclay, France

## ABSTRACT

Medium frequency heterogeneous Helmholtz problem gives rise to severe challenge to existing numerical approaches. Two methods are developed to solve this problem, namely extended Variational Theory of Complex Rays(VTCR) method and wave based weak Trefftz discontinuous Galerkin(WTDG) method. The VTCR belongs to Trefftz method, where general solutions of governing equation are used as shape functions to get approximate solutions for problem in vibroacoustics. When the square of wave number varies linearly in the media, one considers the extended VTCR with shape functions namely Airy wave functions. Then a general way to handle heterogeneous media by the wave based WTDG is proposed. There is no a priori restriction for the wave number. One locally develops general approximate solution of the governing equation, the gradient of the wave number being the small parameter. In this way, zero order and first order approximations are defined, namely zero order WTDG and first order WTDG. Their shape functions only satisfy the local governing equation in average sense. It proves that the extended VTCR and the wave based WTDG propose an efficient and flexible way to solve problem in vibroacoustics.

Keywords: medium frequency, heterogeneous, Helmholtz problem, Trefftz, weak Trefftz

## 1. INTRODUCTION

Nowadays, the numerical simulation has become indispensable to analyze and optimize problems in every part of engineering processes. Such as in automotive industry and aviation industry, abiding by the standards against pollution, the objective of company is to produce a lighter vehicle or airplane with improved comfort for passenger. However, decreasing the weight of vehicle or airplane often leads to the fact that it is more susceptible to vibrations, which are mainly generated by acoustic effect. Therefore the conception of reliable numerical approaches to solve the corresponding problems becomes a core interest.

Depending on characteristics of problem, a vibration in the mechanic field could be classified into three zones, correspondingly low-frequency range, mid-frequency range and high-frequency range. The low-frequency range is characterized by the local response. The resonance peaks are distinct from one to another. The behavior of vibration can be represented by the combination of several normal modes. The Finite Element Method (FEM) is no doubt most commonly used to solve vibration problem in low-frequency range. But this method suffers from pollution effect when frequency increases[1]. In higher frequency range, there exists much more modes and more overlapping resonance peaks. This phenomenon could be observed in mid-frequency and high-frequency range. Thus new approaches need to be proposed. In high-frequency range the dimension of object is much larger than the wave length. The system is extremely sensible to uncertainties. The Statistical Energy Analysis (SEA)[2] is the well-known method in this range. It

<footnote n="1">¹ olivierlihao@gmail.com</footnote>

<footnote n="2">² ladeveze@lmt.ens-cachan.fr</footnote>

<footnote n="3">³ riou@lmt.ens-cachan.fr</footnote>
neglects the local response and it studies the global energy by taking the averages and variances of dynamic field over large sub-systems. Despite its efficiency presented, SEA is limited into the high frequency range. The reason is that the problem in mid-frequency range contains both the characteristics of low-frequency and high-frequency problem. It contains intense modal density but the local response could not be totally neglected. Numerical methods proposed in latter sections are focusing on problems in mid-frequency range.

Facing to mid-frequency problem, approaches could be classified into the extensions of the standard FEM, the energy based methods and Trefftz methods. Among these categories, the Trefftz methods are based on the use of exact approximations of the governing equation. Such methods are, for example, the partition of unity method (PUM)[3], the ultra weak variational method(UWVF)[4], the least square method[5], the discontinuous Galerkin method (DEM)[6], the wave based method[7], the wave boundary element method[8]. The Variational Theory of Complex Rays (VTCR), first introduced in[9], belongs to this category of numerical strategies. Recently, a new approach called the Weak Trefftz Discontinuous Galerkin (WTDG) method is developed in[10,11] to enable the hybrid use of the FEM and VTCR. Most of these Trefftz methods consider the systems as piecewise homogeneous. It is not necessarily the case when the wave number varies in space. Hence this paper will introduce an extended TVRC and a wave based WTDG to solve mid-frequency heterogeneous vibration problem. These two methods are illustrated in the fluid dominated by Helmholtz equation. The variational formulation is introduced in Section 2. Then numerical examples are shown in Section 3. Conclusion is drawn in Section 4.

# 2. Reference Problem and Variational Formulation

## 2.1 Reference Problem

A 2-D Helmholtz problem serves 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$. The problem considered is: find $u \in H^1(\Omega)$ such that

$$ (1-i\eta)\Delta u + k^2(x,y)u + r_d = 0 \quad \text{over } \Omega \\ u = u_d \quad \text{over } \partial\Omega_1 \\ (1-i\eta)\partial_n u = g_d \quad \text{over } \partial\Omega_2 \tag{1} $$

where $\partial_n = \text{grad} \cdot \mathbf{n}$ and $\mathbf{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. $r_d, u_d, g_d$ are the source, prescribed Dirichlet and Neumann data. $k(x,y)$ is the wave number and its value changes with respect to location on $\Omega$.

## 2.2 Variational Formulation

### 2.2.1 Extended VTCR

Let the computational domain $\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 extended VTCR approach consists in searching solution $u$ in functional space $U$ such that

$$ U = \{u | u_{|\Omega_E} \in U_E\} \\ U_E = \{u_E | u_E \in V_E \subset H^1(\Omega_E); (1-i\eta)\Delta u_E + k^2(x,y)u_E + r_d = 0\} \tag{2} $$

when $r_d = 0$, the vector spaces associated with $U$ and $U_E$ are defined as $U_0$ and $U_{E,0}$. Especially when $k^2(x,y) = \alpha x + \beta y + \gamma$, where $\alpha, \beta, \gamma$ are constant parameters, it is proved in [12]
that $\psi(x, y, \mathbf{P}_j)$ namely Airy wave function is the shape function of (2). $\mathbf{P}_j$ denotes the propagation direction of Airy wave.

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)\mathrm{grad}u
\end{aligned} 
\quad (3) 
$$

The variational formulation of (1) can be written as: find $u \in U$ such that

$$ 
\mathrm{Re}\left(-ik\left(\sum_{E,E' \in \mathbb{E}} \int_{\Gamma_{EE'}} \left(\frac{1}{2}\{\mathbf{q}_u \cdot \mathbf{n}\}_{EE'}\{\bar{v}\}_{EE'} - \frac{1}{2}[\bar{\mathbf{q}}_v \cdot \mathbf{n}]_{EE'}[u]_{EE'}\right)dS - \sum_{E \in \mathbb{E}} \int_{\Gamma_{EE \cap \partial \Omega}} \bar{\mathbf{q}}_v \cdot \mathbf{n}(u-u_d)dS \right.\right. \\
\left. \left. + \sum_{E \in \mathbb{E}} \int_{\Gamma_{EE \cap \partial \Omega}} (\mathbf{q}_u \cdot \mathbf{n} - g_d)\bar{v}dS\right)\right) = 0 \quad \forall v \in U_0 
\quad (4) 
$$

### 2.2.2 Wave based WTDG

The wave based WTDG approach consists in searching solution $u$ in functional space $U$ such that

$$ 
\begin{aligned}
U &= \{u | u_{\Omega_E} \in U_E\} \\
U_E &= \{u_E | u_E \in V_E \subset H^1(\Omega_E); (1 - i\eta)\Delta u_E + k_E^2 u_E + r_d = 0\}
\end{aligned} 
\quad (5) 
$$

where $k_E$ is an approximation of $k(x, y)$ in subdomain $\Omega_E$. Developing $k(x, y)$ by Taylor series on the point $(x_e, y_e)$, which is the center of subdomain $\Omega_E$, one generally takes $k_E$ as its zero order approximation or first order approximation. Analogous to formulation proposed in (4), the variational formulation can be written as: find $u \in U$ such that

$$ 
\mathrm{Re}\left(-ik\left(\sum_{E,E' \in \mathbb{E}} \int_{\Gamma_{EE'}} \left(\frac{1}{2}\{\mathbf{q}_u \cdot \mathbf{n}\}_{EE'}\{\bar{v}\}_{EE'} - \frac{1}{2}[\bar{\mathbf{q}}_v \cdot \mathbf{n}]_{EE'}[u]_{EE'}\right)dS - \sum_{E \in \mathbb{E}} \int_{\Gamma_{EE \cap \partial \Omega}} \bar{\mathbf{q}}_v \cdot \mathbf{n}(u-u_d)dS \right.\right. \\
\left. \left. + \sum_{E \in \mathbb{E}} \int_{\Gamma_{EE \cap \partial \Omega}} (\mathbf{q}_u \cdot \mathbf{n} - g_d)\bar{v}dS - \sum_{E \in \mathbb{E}} \int_{\Omega_E} (\mathrm{div}\mathbf{q}_{u_E} + k_E^2 u + r_d)\bar{v}d\Omega\right)\right) = 0 \quad \forall v \in U_0 
\quad (6) 
$$

## 3. Numerical performance

### 3.1 Study of the performance of the extended VTCR on mid frequency heterogeneous Helmholtz problem

The domain $\Omega$ is a simple geometry of square $[0 \text{ m}, 1 \text{ m}] \times [0 \text{ m}, 1 \text{ m}]$. In this domain, $\eta = 0.01$ and $k^2(x, y) = \alpha x + \beta y + \gamma$, where $\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)$. The choice of geometry and boundary conditions allow one to calculate the relative error of extended VTCR method with exact solution. The definition of the problem and the discretization strategy can be seen on Figure 2.
<figure><img src="image_3.png" /><figcaption>Figure 2 – Example considered in Section 3.1. From left to right, First: definition of domain, Second: 1subdomain discretization, Third: 4 subdomains discretization, Fourth: 9 subdomains discretization.</figcaption></figure>

The relative error is defined as:

$$|u_{exact} - u| / |u_{exact}| = \sqrt{\int_{\Omega} |u_{exact} - u|^2 d\Omega / \int_{\Omega} |u_{exact}|^2 d\Omega} \quad (7)$$

As one could see from Figure 3, the extended VTCR only requires a small amount of degrees of freedom to make the simulation result converge. This phenomenon is consistent with the performance of VTCR when the medium is homogeneous[13]. There exists two 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 Figure 3 extended VTCR with only one computational domain converges the fastest.

<figure><img src="image_4.png" /><figcaption>Figure 3 – The convergence curves for the example of Section 3.1.</figcaption></figure>

## 3.2 Study of the performance of the wave based WTDG on mid frequency heterogeneous Helmholtz problem

The geometry definition keeps the same as the one in Section 3.1 and $\eta = 0.01$. However, in
this example one has $k(x,y) = 5x + 5y + 40$. In this case, the exact solution of governing equation in (1) is unknown, therefore one could not use extended VTCR to solve it. By smoothly approximating the governing equation, the wave based WTDG could handle this problem. The boundary condition on $\partial\Omega$ are Dirichlet type such that $u_d = 1$.

<figure><img src="image_5.png" /><figcaption>Figure 4 – Results of example considered in Section 3.2. From left to right: First: the zero order WTDG with 4 subdomains and 100 waves per subdomain. Second: the zero order WTDG with 25 subdomains and 80 waves per subdomain. Third: the zero order WTDG with 100 subdomains and 40 waves per subdomain.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5 – Results of example considered in Section 3.2. From left to right: First the first order WTDG with 1 subdomain and 160 waves per subdomain. Second: the first order WTDG with 4 subdomains and 120 waves per subdomain. Third: Referential solution calculated by the FEM with 625 elements of quadric mesh of order 3.</figcaption></figure>

One could see from Figure 4 and Figure 5 that both the zero order WTDG and the first order WTDG could lead to a convergent result. In addition, the first order WTDG converges much faster than the zero order WTDG. The reason is that the first order WTDG makes a higher order approximation of the governing equation. In Figure 5, there is a referential result calculated by the FEM. It could be observed that results obtained by the WTDG are consistent with the reference one.
## 3.3 Semi-unbounded harbor agitation problem

<figure><img src="image_7.png" /><figcaption>Figure 6 – Example considered in Section 3.3. Left: top view of Harbor. Right: profile of seafloor inside the harbor</figcaption></figure>

This example focuses on the water agitation of harbor[12,14]. In this problem the movement of waves yield to the Helmholtz equation. The definition of harbor can be seen in Figure 6. The incident wave comes with angle $\theta_0^+$. The origin of coordinate is O, located in the middle point of the harbor entrance. Outside the harbor the depth of seafloor is nearly constant. Inside the harbor the depth of seafloor varies as Figure 6 shows and it could be modelled as:

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

where a, b are constant parameters. All boundaries of the harbor are supposed to be reflecting boundaries, which are denoted by $\Gamma_R$:

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

The relation between wave frequency $\omega$ and water depth $h$ follows the non-linear dispersion relation:

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

where $g = 9.81 \text{ m/s}^2$ is the gravitational acceleration and $k$ is the wave number. In the case that the depth of water is much less than the length of wave, there is the shallow water approximation:

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

Associating (8) and (11), one could obtain that:

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

For the waves propagate towards far away field, they need to satisfy also Sommerfeld radiation condition.
<figure><img src="image_8.png" /><figcaption>Figure 7 – Discretization strategies in Section 3.3. Left: Strategy of the VTCR. Right: Strategy of the zero order WTDG</figcaption></figure>

As Figure 7 shows, both the extended VTCR and the zero order WTDG take the same computational strategies on $\Omega_1$ and $\Omega_2$:

$$U_{\Omega_1} = \{u(x,y): u(x,y) = A_0^+ e^{ik_0^+(\cos\theta_a x + \sin\theta_a y)}, \sum_{n=0}^{N_1} A_{1n} \cos n\theta H_n^{(1)}(kr), A_0^+, A_{1n} \in C, n = 0, \dots, N_1\} \quad (13)$$

$$U_{\Omega_2} = \{u(x,y): u(x,y) = \sum_{n=0}^{N_2} A_{2n} e^{ik(\cos\theta_a x + \sin\theta_a y)}, A_{2n} \in C, n = 0, \dots, N_2\} \quad (14)$$

where $A_0^+, A_{1n}$ are the amplitudes of waves on $\Omega_1$ and $A_{2n}$ is the amplitude of wave on $\Omega_2$. $\theta_a = \pi - \theta_0^+$ where $\theta_0^+$ is the angle of incident wave. $H_n^{(1)}$ is the first type of Hankel function.

Inside the harbor, the extended VTCR only has one subdomain and the working space is in form of:

$$U_{\Omega_3} = \{u(x,y): u(x,y) = \sum_{n=0}^{N_3} A_{3n} \psi(x,y, \mathbf{P}_{3n}), A_{3n} \in C, n = 0, \dots, N_3\} \quad (15)$$

where $A_{3n}$ is the amplitude of wave on $\Omega_3$. For the zero order WTDG, the domain inside the harbor is discretized into several subdomains and the working space is in form of:

$$U_{\Omega_3} = \{u(x,y): u(x,y) = \sum_{n=0}^{N_3} A_{jn} e^{ik_j(\cos\theta_a x + \sin\theta_a y)}, A_{jn} \in C, n = 0, \dots, N_j\} \quad (14)$$

where $k_j = k(x_j, y_j)$ and $(x_j, y_j)$ is the coordinate of center point on $\Omega_j$. $A_{jn}$ is the amplitude of wave on $\Omega_j$.

Parameters of model are chosen as $a = 4.8 \cdot 10^{-2} \text{ m}^{-1}$, $b = 4.8 \cdot 10^{-5} \text{ m}^{-2}$, $\eta = 0.03$. The amplitude of incoming wave corresponds to $A_0^+ = 2 \text{ m}$ and the angle of incoming wave $\theta_0^+ = 45^\circ$. Results could be seen in Figure 8 and Figure 9. Figure 8 shows the global result and Figure 9 shows only the result inside the harbor. It could be seen that results calculated by the extended VTCR and the zero order WTDG are consistent.
<figure><img src="image_9.png" /><figcaption>Figure 8: Example considered in Section 3.3. Up left: extended VTCR. Up right: zero order WTDG with five subdomains inside the harbor. Down left: zero order WTDG with ten subdomains inside the harbor. Down right: zero order WTDG with fifteen subdomains inside the harbor.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Example considered in Section 3.3. Up left: extended VTCR. Up right: zero order WTDG with five subdomains inside the harbor. Down left: zero order WTDG with ten subdomains inside the harbor. Down right: zero order WTDG with fifteen subdomains inside the harbor.</figcaption></figure>

## 4. CONCLUSIONS

This paper proposes an extended VTCR method and a wave based WTDG method, which are both able to solve mid-frequency heterogeneous Helmholtz problem. The equivalent variational
formulations are introduced. Academy study is done in this paper to show the convergence property of the proposed methods. Furthermore, these two approaches are applied to solve a harbor agitation problem. Results show that both methods have the potentiality to solve engineering problem in an efficient and flexible way.

## ACKNOWLEDGEMENTS

The author acknowledges the financial support by Labex CeLyA of Université de Lyon, operated by the French National Research Agency (ANR-10-LABX-0060/ANR-11-IDEX-0007). The author also acknowledges the support of work by Laboratoire de mécanique et technologie of Ecole Normale Supérieure Paris-Saclay. These supports are greatly appreciated.

## REFERENCES

1. Deraemaeker A, Babuška I, Bouillard P. Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions[J]. International journal for numerical methods in engineering, 1999; 46(4): 471-499.

2. Lyon R H, Maidanik G. Power flow between linearly coupled oscillators[J]. The journal of the Acoustical Society of America, 1962; 34(5): 623-639.

3. Strouboulis T, Hidajat R. Partition of unity method for Helmholtz equation: q-convergence for plane-wave and wave-band local bases[J]. Applications of Mathematics, 2006; 51(2): 181-204.

4. Cessenat O, Despres B. Application of an ultra weak variational formulation of elliptic PDEs to the two-dimensional Helmholtz problem[J]. SIAM journal on numerical analysis, 1998; 35(1): 255-299.

5. Monk P, Wang D Q. A least-squares method for the Helmholtz equation[J]. Computer Methods in Applied Mechanics and Engineering, 1999; 175(1-2): 121-136.

6. Gittelson C J, Hiptmair R, Perugia I. Plane wave discontinuous Galerkin methods: analysis of the h-version[J]. ESAIM: Mathematical Modelling and Numerical Analysis, 2009; 43(2): 297-331.

7. Desmet W, Sas P, Vandepitte D. An indirect Trefftz method for the steady-state dynamic analysis of coupled vibro-acoustic systems[J]. Computer Assisted Mechanics and Engineering Sciences, 2001; 8(2-3): 271-288.

8. Perrey-Debain E, Trevelyan J, Bettess P. Wave boundary elements: a theoretical overview presenting applications in scattering of short waves[J]. Engineering Analysis with Boundary Elements, 2004; 28(2): 131-141.

9. Ladevèze P. A new computational approach for structure vibrations in the medium frequency range[J]. Comptes Rendus de L Academie des Sciences Serie II Fascicule B-Mecanique Physique Chimie Astronomie, 1996; 322(12): 849-856.

10. Ladevèze P, Riou H. On Trefftz and weak Trefftz discontinuous Galerkin approaches for medium-frequency acoustics[J]. Computer Methods in Applied Mechanics and Engineering, 2014; 278: 729-743.

11. Li H, Ladevèze P, Riou H. Hybrid Finite Element Method and Variational Theory of Complex Rays for Helmholtz Problems[J]. Journal of Computational Acoustics, 2016; 24(04): 1650015.

12. Li H, Ladevèze P, Riou H. Extended variational theory of complex rays in heterogeneous Helmholtz problem[J]. Computational Mechanics, 2017; 6(59): 909-918.

13. Riou H, Ladevèze P, Sourcis B. The multiscale VTCR approach applied to acoustics problems[J]. Journal of Computational Acoustics, 2008; 16(04): 487-505.

14. Modesto D, Zlotnik S, Huerta A. Proper generalized decomposition for parameterized Helmholtz problems in heterogeneous and unbounded domains: Application to harbor agitation[J]. Computer Methods in Applied Mechanics and Engineering, 2015; 295: 127-149.