NiHu  2.0
An introductory BEM example

Table of Contents

Problem definition

Let us consider the Helmholtz equation with a Neumann boundary condition:

\( \displaystyle \left(\nabla^2 + k^2 \right) p({\bf x}) = 0 \qquad {\bf x} \in V \\ \displaystyle \partial p({\bf x}) / \partial n = q({\bf x}) = \bar{q}({\bf x}) \qquad {\bf x} \in S \)

where \( V\subset\mathbb{R}^d \) is a bounded domain with a smooth boundary \( S \) , \( p \) denotes the acoustic pressure, \( q \) denotes its normal derivative (velocity), \( \bar{q} \) denotes the prescribed normal velocity on the boundary \( S \) , and \( k \) denotes the acoustic wave number. The pressure field \( p({\bf x}) \) is sought on the boundary and in internal points of the domain \( V \) .

Fundamental solution and integral operators

Solution of the boundary value problem is formulated using the fundamental solution \( G({\bf x}_0, {\bf x}) \) of the Helmholtz equation defined by

\( \displaystyle \left(\nabla^2 + k^2 \right) G({\bf x}_0, {\bf x}) = -\delta({\bf x}-{\bf x}_0), \qquad {\bf x} \in \mathbb{R}^d \setminus \left\{ {\bf x}_0 \right\} \)

where \( \delta \) denotes the Dirac delta function. The fundamental solution and its normal derivative are used as kernel functions to form the following single and double layer potential integral operators, respectively

\( \displaystyle \left(\mathcal{L}u\right)_S({\bf x}) = \int_S G({\bf x},{\bf y}) u({\bf y}) \mathrm{d} S_{{\bf y}} \\ \displaystyle \left(\mathcal{M}u\right)_S({\bf x}) = \int_S \frac{\partial G({\bf x}, {\bf y})}{\partial n_{{\bf y}}}({\bf x},{\bf y}) u({\bf y}) \mathrm{d} S_{{\bf y}} \)

The introduced operators are applied to pose an integral equation problem equivalent to the Helmholtz problem:

\( \displaystyle \left(\mathcal{M}p\right)_S({\bf x}) - \left(\mathcal{L}q\right)_S({\bf x}) = \frac{1}{2} p({\bf x}) \qquad {\bf x} \in S \\ \displaystyle \left(\mathcal{M}p\right)_S({\bf x}) - \left(\mathcal{L}q\right)_S({\bf x}) = p({\bf x}) \qquad {\bf x} \in V \setminus S \)

The first, boundary integral equation needs to be solved, using the prescribed velocity field \( \bar{q}({\bf x}) \) as excitation, to obtain the acoustic pressure \( p({\bf x}) \) on the boundary \( S \) . Consecutively, the acoustic pressure field \( p({\bf x}) \) in \( V \) can be determined by evaluating the second equation.

Weighted residuals

The solution of the boundary integral equation is based on the weighted residual method: The boundary integral equation is assumed to hold in a weak sense:

\( \displaystyle \left< v, \left(\mathcal{M}p\right)_S - \left(\mathcal{L}q\right)_S - \frac{1}{2}p \right>_S = 0 \)

where \( \left<\cdot,\cdot\right>_S \) denotes the inner product defined as

\( \displaystyle \left<v,f\right>_S = \int_{S} v({\bf x}) f({\bf x}) \mathrm{d} S_{{\bf x}} \)

and \( v({\bf x}) \) denotes an arbitrary test function chosen from the test function space \( T \subset S \to \mathbb{R} \) .

Discretisation

Numerical solution of the weak form is found by discretising the test space \( T \) and the domain spaces of the boundary integral operators \( \mathcal{L} \) and \( \mathcal{M} \) . Discretisation is performed by introducing their finite sets of basis functions \( t_i \) , \( d^\mathcal{L}_j \) and \( d^\mathcal{M}_j \) :

\( \displaystyle v \in T \to v({\bf x}) = \sum_{i} v_i t_i({\bf x}) \\ \displaystyle p \in D_{\mathcal{M}} \to p({\bf y}) = \sum_{j} p_j d^{\mathcal{M}}_j({\bf y}) \\ \displaystyle q \in D_{\mathcal{L}} \to q({\bf y}) = \sum_{l} q_l d^{\mathcal{L}}_l({\bf y}) \)

The discretised weak form of the boundary integral equation simplifies to a linear system of equations

\( \displaystyle \sum_j \left< t_i, \left(\mathcal{M}d^{\mathcal{M}}_j\right)_S \right>_S p_j - \sum_l \left< t_i, \left(\mathcal{L}d^{\mathcal{L}}_l\right)_S \right>_S q_l = \frac{1}{2} \sum_j \left< t_i, d^{\mathcal{M}}_j \right>_S p_j \)

or, with matrix-vector notations,

\( \displaystyle {\bf M} {\bf p} - {\bf L} {\bf q} = \frac{1}{2} {\bf B} {\bf p} \)

where \( {\bf M} \) is a square dense matrix, \( {\bf L} \) is a not necessarily square dense matrix and \( {\bf B} \) is a sparse square matrix. The Neumann problem is then solved by rearranging and inverting the linear equation

\( \displaystyle {\bf p} = \left( {\bf M} - \frac{1}{2} {\bf B} \right)^{-1} {\bf L} \bar{{\bf q}} \)

BEM formalisms

Different BEM formalisms can be defined by classifying the selection of the function spaces.

  • Galerkin BEM methods refer to the selection case, where the test function space coincides with the domain space of the solution: \( t_i = d_i^{\mathcal{M}} \) . In this selection case, the inner products are evaluated by double surface integrals, and, applying fem terminology, matrix \( \bf B \) simplifies to a ,,boundary mass matrix''.
  • Collocational BEM methods denote the case, where the test function space is generated by Dirac delta functions located at nodal locations \( {\bf x}_i \) of the solution: \( t_i({\bf x}) = \delta_{{\bf x}_i} = \delta\left({{\bf x}-{\bf x}_i}\right) \) . In the collocational case, the inner products boils down to single surface integrals, and matrix \( \bf B \) typically simplifies to the unity matrix.

After the solution has been performed by means of direct or iterative techniques, the radiated field \( p({\bf x}_i), {\bf x}_i \in V \) can be computed by evaluating the discretised version of the radiation integral

\( \displaystyle \sum_j \left< \delta_{{\bf x}_i}, \left(\mathcal{M}d^{\mathcal{M}}_j\right)_S \right>_F p_j - \sum_k \left< \delta_{{\bf x}_i}, \left(\mathcal{L}d^{\mathcal{L}}_l\right)_S \right>_F q_l = p({\bf x}_i) \)