Table of Contents
Introduction
This tutorial demonstrates the usage of NiHu for solving the Helmholtz equation in 3D by means of a Galerkin type boundary element method (BEM). In this tutorial the Helmholtz equation will be applied for the solution of an acoustic wave propagation problem. In this tutorial we will use NiHu's Matlab interface for the mesh definition and the solution of the resulting matrix equations. This tutorial is based on the theory and notations introduced previously in the theoretical introduction BEM example.
Theory
In this example the Helmholtz integrals are utilised in order to solve an external radiation problem. Our radiating surface is a closed surface in the three-dimensional space \( S \subset \mathbb{R}^3 \). Neumann boundary conditions are specified on the entire surface, given as \( q(\mathbf{x}) = \bar{q}(\mathbf{x}) \) if \( \mathbf{x} \in S \). We need to determine the radiated acoustic field \( p(\mathbf{x}) \) on the surface \( S \) and in an other set of points \( F \subset \mathbb{R}^3 \setminus S \) referred to as field points.
The Helmholtz integrals
The sound pressure radiated by the closed vibrating surface \( S \) is obtained by means of the Helmholtz integrals
\( \displaystyle \frac{1}{2} p({\bf x}) = \int_S H({\bf x}, {\bf y}) p({\bf y}) dS_y - \int_S G({\bf x}, {\bf y}) q({\bf y}) dS_y \qquad \) if \( \quad x \in S \)
\( \phantom{\displaystyle \frac{1}{2}} p({\bf x}) = \displaystyle \int_S H({\bf x}, {\bf y}) p({\bf y}) dS_y - \int_S G({\bf x}, {\bf y}) q({\bf y}) dS_y \qquad \) if \( \quad x \in F , \)
where \( G({\bf x}, {\bf y}) \) and \( H({\bf x}, {\bf y}) \) denote the fundamental solution of the inhomogeneous Helmholtz equation and its normal derivative with respect to the source point. The kernel functions are given as
\( G({\bf x}, {\bf y}) = \displaystyle \frac{\mathrm{e}^{-\mathrm{i}kr}}{4\pi r} \quad \) and \( \quad H({\bf x}) = \displaystyle \frac{\partial G(\mathbf{x}, \mathbf{y})}{\partial n_y} = -G({\bf x}, {\bf y}) \left(\displaystyle \frac{1 + \mathrm{i} k r}{r} \frac{{\bf r} \cdot {\bf n}_y }{r} \right) , \)
with the distance vector \( \bf r \) and the distance \( r \) defined as \( \mathbf{r} = \mathbf{y} - \mathbf{x} \) and \( r = \left| \mathbf{r} \right| . \) The normal vector \( \mathbf{n}_y \) is pointing outward from the surface \( S \) at the point \( \mathbf{y} \).
Galerkin BEM formalism
In order to solve the Helmholtz integral equations numerically, the radiating surface \( S \) is divided into a number of non-overlapping elements. The sound pressure function \( p(\mathbf{x}) \) and its normal derivative \( q(\mathbf{x}) \) are approximated by means of weighting functions \( w(\mathbf{x}) \), also called trial functions, such that
\( p(\mathbf{x}) \approx \displaystyle \sum_{j=1}^n w_j (\mathbf{x}) p_j \qquad \) and \( \qquad q(\mathbf{x}) \approx \displaystyle \sum_{j=1}^n w_j (\mathbf{x}) q_j , \)
with \( p_j \) and \( q_j \) denoting the corresponding weights and \( n \) representing the total number of weighting functions, also known as the number of degrees of freedom.
As a next step, the weak form of the Helmholtz integral equation is obtained by multiplying the equations with a set of test functions \( t_j(\mathbf{x}) \) and integrating the result over the complete surface \( S \) with respect to the variable \( \mathbf{x} \). In the Galerkin formulation, the test function are equivalent to the weighting functions, i.e. \( t_j(\mathbf{x}) \equiv w_j(\mathbf{x}) \). Hence, the Galerkin formalism is a special case of weighted residuals with the weighting functions being equivalent to the trial functions.
Applying the operator notation introduced for weighted residuals , the double integral for the surface is expressed in the form
\( \left< w_i , \left(\left[ \mathcal{M} - \frac{1}{2} \mathcal{I} \right] w_j \right)_S \right>_S p^{(s)}_j = \left< w_i, \left(\mathcal{L} w_j \right)_S \right>_S q^{(s)}_j . \)
The upper inidices \( \cdot^{(s)} \) and \( \cdot^{(f)} \) symbolise the variables on the surface \( S \) and the field \( F \) meshes, respectively. For the calculation of the acoustic field of the radiated surface a Dirac-delta test function is applied herein, which gives
\( p^{(f)}_i = \left< \delta_i , \left( \mathcal{M} w_j \right)_S \right>_F p^{(s)}_j - \left< \delta_i, \left(\mathcal{L} w_j \right)_S \right>_F q^{(s)}_j . \)
The equations can be rewritten by utilizing matrix notations
\( M^{(s)}_{ij} = \left< w_i , \left(\left[ \mathcal{M} - \frac{1}{2} \mathcal{I} \right] w_j \right)_S \right>_S \), \( \qquad L^{(s)}_{ij} = \left< w_i, \left(\mathcal{L} w_j \right)_S \right>_S \),
\( M^{(f)}_{ij} = \left< \delta_i , \left( \mathcal{M} w_j \right)_S \right>_F \qquad \) and \( \qquad L^{(f)}_{ij} = \left< \delta_i , \left( \mathcal{L} w_j \right)_S \right>_F \)
in the forms
\( \mathbf{M}^{(s)} \mathbf{p}^{(s)} = \mathbf{L}^{(s)} \mathbf{q}^{(s)} \qquad \) and
\( \phantom{\mathbf{M}^{(s)}} \mathbf{p}^{(f)} = \mathbf{M}^{(f)} \mathbf{p}^{(s)} - \mathbf{L}^{(f)} \mathbf{q}^{(s)} .\)
The vectors \( \mathbf{p} \) and \( \mathbf{q} \) denote the column vectors of the weights for the degrees of freedom of the pressure and its normal derivative, respectively.
In the sequel, the above Matrix equations are determined and solved using NiHu in the following steps
- Assembling the system matrices
- Solution of the surface equation
- Evaluation of the radiation equation for the field points
- Displaying the results and quantifying the error
Step (1) is performed in C++, by means of a mex function, whereas steps (2-4) are carried out utilising the Matlab interface of NiHu.
The C++ code
The objective of our C++ code is to assemble the systems matrices in a function callable from Matlab by means of mex.
Header and mesh creation
As it was done in the Rayleigh integral tutorial the header of our C++ source file includes the necessary header files and defines some basic types for convenience.
Since we have to calculate four different matrices, that are \( \mathbf{M}^{(s)} \), \( \mathbf{L}^{(s)} \), \( \mathbf{M}^{(f)} \) and \( \mathbf{L}^{(f)} \), our mex function will return four left hand side parameters.
Similar to the Rayleigh integral tutorial, two meshes will be created. Thus, our mex function header and the mesh generation commands will read as
Definition of function spaces
Next, we create the discretised function spaces of our mesh. In this example we will use homogeneous function spaces, which are easily created by means of function space views of the meshes. On the radiating surface, a piecewise constant approximation of trial and test functions are applied, whereas on the field point mesh, a piecewise constant approximation of the trial function with Dirac test functions is applied. The Dirac delta test functions are created by using the simple dirac view approach.
The number of degrees of freedom on the radiating surface \( S \) and in the field mesh \( F \) are denoted by \( n \) and \( m \), respectively. After the function spaces are built, the number of DOFs are known, and the memory for the system matrices can be preallocated. The surface system matrices \( \mathbf{M}^{(s)} \) and \( \mathbf{L}^{(s)} \) are of size \( n \times n \), whereas the field system matrices \( \mathbf{M}^{(f)} \) and \( \mathbf{L}^{(f)} \) are of size \( m \times n \). The four system matrices are preallocated by means of the following lines of code.
Integral operators and their evaluation
In the next steps, the thre integral operators \( \mathcal{M} \), \( \mathcal{L} \) and \( \mathcal{I} \) are defined. The operators \( \mathcal{M} \) and \( \mathcal{L} \) are created as the double and single layer potential kernels of the three-dimensional Helmholtz equations, while operator \( \mathcal{I} \) is simply the identity operator. Since the kernel functions depend on the wave number \( k \), which is passed as the fifth right hand side parameter of the mex function, the integral operators are created as follows.
Finally, the system matrices are obtained by the evaluation of the integral operators on the function spaces defined above.
It should be realised that the second line of the above code snippet contains the evaluation two integral operators, \( \mathcal{M} \) and \( \mathcal{I} \) and stores the summed result in the matrix \( \mathbf{M}^{(s)} \).
The resulting mex function is capable of computing the four system matrices for arbitatry radiating surface meshes consisting of quadrilateral and triangular elements and field meshes consisting of quadrilateral elements only at a give wave number \( k \).
The Matlab code
The example creates a sphere surface of radius \( R = 1\,\mathrm{m} \) centered at the origin, consisting of linear quadrilateral and triangular elements. The field mesh is defined as a disc segment, parametrised as \( 1.125\,\mathrm{m} \le r \le 3.625\,\mathrm{m} \), \( 0 \le \varphi \le \pi / 2 \) and \( \vartheta = 0 \) .
The Neumann boundary conditions are defined by a point source located at the point \( \mathbf{x}_0 \). By this definition the geometry is considered transparent, and the resulting acoustic fields both on \( S \) and \( F \) are obtained as the field of the point source. Therefore it is easy to compare the calculated numerical results to the analytical ones.
The compiled mex file can be called from Matlab, as demonstrated in the example below:
The generated plot looks like
And the resulting errors read as
Surface log10 error: -2.079379 Field log10 error: -2.448697
The full codes of this tutorial are available here:
- C++ code: tutorial/helmholtz_bem_3d.mex.cpp
- Matlab code: tutorial/helmholtz_bem_3d_test.m
