# Numerical Solution of the Helmholtz Equation with Robin Boundary Conditions using the Finite Element Method

This notebook illustrates the numerical solution of the wave equation for an harmonic excitation and Robin boundary conditions using the [Finite Element Method](https://en.wikipedia.org/wiki/Finite_element_method) (FEM). The method aims at an approximate solution by subdividing the area of interest into smaller parts with simpler geometry, linking these parts together and applying methods from the calculus of variations to solve the problem numerically. The FEM is a well established method for the numerical approximation of the solution of partial differential equations (PDEs). The solutions of PDEs are often known analytically only for rather simple geometries. FEM based simulations allow to gain insights into other more complex cases.

## Problem Statement

The inhomogeneous [Helmholtz equation](https://en.wikipedia.org/wiki/Helmholtz_equation) is given as

\begin{equation}
\Delta P(\mathbf{x}, \omega) + \frac{\omega^2}{c^2} P(\mathbf{x}, \omega) = - Q(\mathbf{x}, \omega) .
\end{equation}

We aim for a numerical solution of the Helmholtz equation on the domain $V$ with respect to the homogeneous Robin boundary condition

\begin{equation}
V_n(\mathbf{x}, \omega) + \frac{1}{Z(\mathbf{x}, \omega)} P(\mathbf{x}, \omega) = 0 \qquad \text{for } x \in \partial V ,
\end{equation}

where $V_n(\mathbf{x}, \omega)$ denotes the particle velocity in inward normal direction to the boundary $\partial V$ of $V$ and $Z(\mathbf{x}, \omega)$ the acoustic impedance of the boundary.
The particle velocity can be linked to the pressure using the Euler equation

\begin{equation}
-\mathrm{j} \omega \rho_0 V_n(\mathbf{x}, \omega) = \frac{\partial}{\partial n} P(\mathbf{x}, \omega) ,
\end{equation}

where $\rho_0$ denotes the static density of air.
Introducing this into the Robin boundary equation above yields

\begin{equation}
\frac{\partial}{\partial n} P(\mathbf{x}, \omega) - \mathrm{j} \underbrace{\frac{\omega \rho_0}{Z}}_{\sigma} P(\mathbf{x}, \omega) = 0 \qquad \text{for } x \in \partial V .
\end{equation}

The medium impedance of air for free-field propagation is $Z_0 = \rho_0 c$, hence $\sigma_0 = \frac{\omega}{c}$ in this case. Free-field conditions can be simulated by matching the impedance of the boundary to $Z_0$.

## Variational Formulation

Starting from the [variational formulation of the Helmholtz equation](FEM_Helmholtz_equation_2D.ipynb#Variational-Formulation) (before application of Green's first theorem)

\begin{equation}
{-} \int_V \nabla P(\mathbf{x}, \omega) \cdot \nabla V(\mathbf{x}, \omega) \mathrm{d}x  + \int_{\partial V} V(\mathbf{x}, \omega) \frac{\partial}{\partial n}  P(\mathbf{x}, \omega) \mathrm{d}s + \frac{\omega^2}{c^2} \int_V P(\mathbf{x}, \omega) V(\mathbf{x}, \omega) \mathrm{d}x = -\int_V Q(\mathbf{x}, \omega) V(\mathbf{x}, \omega) \mathrm{d}x
\end{equation}

and introducing the Robin boundary condition into the second integral yields

\begin{equation}
{-} \int_V \nabla P(\mathbf{x}, \omega) \cdot \nabla V(\mathbf{x}, \omega) \mathrm{d}x  + \mathrm{j} \sigma \int_{\partial V} V(\mathbf{x}, \omega) P(\mathbf{x}, \omega) \mathrm{d}s + \frac{\omega^2}{c^2} \int_V P(\mathbf{x}, \omega) V(\mathbf{x}, \omega) \mathrm{d}x = -\int_V Q(\mathbf{x}, \omega) V(\mathbf{x}, \omega) \mathrm{d}x .
\end{equation}

It is common to express this integral equation in terms of the bilinear $a(P, V)$ and linear $L(V)$ forms 

\begin{equation}
a(P, V) = \frac{\omega^2}{c^2} \int_V P(\mathbf{x}, \omega) V(\mathbf{x}, \omega) \mathrm{d}x - \int_V \nabla P(\mathbf{x}, \omega) \cdot \nabla V(\mathbf{x}, \omega) \mathrm{d}x + \mathrm{j} \sigma \int_{\partial V} V(\mathbf{x}, \omega) P(\mathbf{x}, \omega) \mathrm{d}s , \end{equation}

\begin{equation}
L(V) = -\int_V Q(\mathbf{x}, \omega) V(\mathbf{x}, \omega) \mathrm{d}x ,
\end{equation}

where

\begin{equation}
a(P, V) = L(V) .
\end{equation}

Computational implementations of the FEM (like FEniCS) may not be able to handle complex numbers.
In this case the problem can be split into two coupled problems with respect to the real and imaginary part. By introducing $P(\mathbf{x}, \omega) = P_r(\mathbf{x}, \omega) + \mathrm{j} P_i(\mathbf{x}, \omega)$ and $V(\mathbf{x}, \omega) = V_r(\mathbf{x}, \omega) + \mathrm{j} V_i(\mathbf{x}, \omega)$ and identifying the real and imaginary parts of the bilinear and linear forms, we get

\begin{equation}
a_r = \int_V \left( \frac{\omega^2}{c^2} V_r(\mathbf{x}, \omega) P_r(\mathbf{x}, \omega) -  \frac{\omega^2}{c^2} V_i(\mathbf{x}, \omega) P_i(\mathbf{x}, \omega) - \nabla P_r(\mathbf{x}, \omega) \cdot \nabla V_r(\mathbf{x}, \omega) + \nabla P_i(\mathbf{x}, \omega) \cdot \nabla V_i(\mathbf{x}, \omega) \right) \mathrm{d}x -
\sigma \cdot \int_{\partial V} \left( V_i(\mathbf{x}, \omega) P_r(\mathbf{x}, \omega) + V_r(\mathbf{x}, \omega) P_i(\mathbf{x}, \omega) \right) \mathrm{d}s, 
\end{equation}

\begin{equation}
a_i = \int_V \left( \frac{\omega^2}{c^2} V_r(\mathbf{x}, \omega) P_i(\mathbf{x}, \omega) +  \frac{\omega^2}{c^2} V_i(\mathbf{x}, \omega) P_r(\mathbf{x}, \omega) - \nabla P_i(\mathbf{x}, \omega) \cdot \nabla V_r(\mathbf{x}, \omega) + \nabla P_r(\mathbf{x}, \omega) \cdot \nabla V_i(\mathbf{x}, \omega) \right) \mathrm{d}x +
\sigma \cdot \int_{\partial V} \left( V_r(\mathbf{x}, \omega) P_r(\mathbf{x}, \omega) - V_i(\mathbf{x}, \omega) P_i(\mathbf{x}, \omega) \right) \mathrm{d}s
\end{equation}

for the bilinear form.

## Numerical Solution

The numerical solution of the variational problem is based on [FEniCS](https://fenicsproject.org/), an open-source framework for numerical solution of PDEs.
Its high-level Python interface `dolfin` is used in the following to define the problem and compute the solution.
The implementation is based on the variational formulation derived above.
It is common in the FEM to denote the solution of the problem by $u$ and the test function by $v$.
The definition of the problem in FEniCS is very close to the mathematical formulation of the problem.

For the subsequent examples the solution of inhomogeneous wave equation for a point source $Q(\mathbf{x}) = \delta(\mathbf{x}-\mathbf{x_s})$ at position $\mathbf{x_s}$ is computed using the FEM.
A function is defined for this purpose, as well as for the plotting of the resulting sound field.


```python
import dolfin
import mshr
import matplotlib.pyplot as plt
%matplotlib inline


def Helmholtz_Robin(mesh, frequency, xs, sigma=dolfin.Constant(0), c=343):
    
    # squared wavenumber
    k2 = dolfin.Constant(2*dolfin.pi*frequency/c)**2
    
    # define function space
    V = dolfin.VectorFunctionSpace(mesh, "CG", 1, dim=2)

    # define variational problem
    (u_r, u_i) = dolfin.TrialFunction(V)
    (v_r, v_i) = dolfin.TestFunction(V)

    a_r = ( k2 * dolfin.inner(u_r,v_r) - k2 * dolfin.inner(u_i,v_i) - dolfin.inner(dolfin.grad(u_r), dolfin.grad(v_r)) + dolfin.inner(dolfin.grad(u_i), dolfin.grad(v_i)) ) * dolfin.dx  - sigma*dolfin.inner(u_r, v_i) * dolfin.ds - sigma*dolfin.inner(u_i, v_r) * dolfin.ds
    a_i = ( k2 * dolfin.inner(u_r,v_i) + k2 * dolfin.inner(u_i,v_r) - dolfin.inner(dolfin.grad(u_r), dolfin.grad(v_i)) - dolfin.inner(dolfin.grad(u_i), dolfin.grad(v_r)) ) * dolfin.dx  + sigma*dolfin.inner(u_r, v_r) * dolfin.ds - sigma*dolfin.inner(u_i, v_i) * dolfin.ds
    L_r = dolfin.Constant(0) * v_r * dolfin.dx
    L_i = dolfin.Constant(0) * v_i * dolfin.dx

    a = a_r + a_i
    L = L_r + L_i
    A, b = dolfin.assemble_system(a, L)

    # define inhomogenity
    delta = dolfin.PointSource(V.sub(0), xs, -1)  # negative amplitude accounts for -Q(x,w) in inhomogeneous wave equation
    delta.apply(b)

    # compute solution
    u = dolfin.Function(V)
    dolfin.solve(A, u.vector(), b)
    (u_r, u_i) = dolfin.split(u)
    
    return u_r


def plot_soundfield(u):
    '''plots solution of FEM-based simulation'''
    fig = plt.figure(figsize=(10,10))
    fig = dolfin.plot(u)
    plt.title(r'$P(\mathbf{x}, \omega)$')
    plt.xlabel(r'$x$ in m')
    plt.ylabel(r'$y$ in m')
    plt.colorbar(fig, fraction=0.038, pad=0.04);
```

### Sound Field in a Rectangular Room

The two-dimensional sound field in a rectangular room (rectangular plate) with homogeneous Robin boundary conditions is computed for the frequency $f=1000$ Hz and source position $x_s = (1.2,3.2)$ m.


```python
f = 1000  # frequency
xs = dolfin.Point(1.2, 3.2)  # source position

# define geometry and mesh
mesh = dolfin.RectangleMesh(dolfin.Point(0, 0), dolfin.Point(5, 4), 200, 200, "right/left")
```

First, the case of sound-hard (Neumann) boundary conditions with $\sigma = 0$ ($Z \to \infty$) is considered.


```python
# compute solution for sigma=0
u = Helmholtz_Robin(mesh, f, xs, sigma=dolfin.Constant(0))

# plot sound field
plot_soundfield(u)
plot_soundfield(abs(u))
plt.title(r'$|P(\mathbf{x}, \omega)|$');
```

Now the case of a matched boundary (free-field propagation) with $\sigma = \frac{\omega}{c}$ ($Z = Z_0$) is considered.


```python
# compute solution for free-field propogation
u = Helmholtz_Robin(mesh, f, xs, sigma=dolfin.Constant(2*dolfin.pi*f/343))

# plot sound field
plot_soundfield(u)
plot_soundfield(abs(u))
plt.title(r'$|P(\mathbf{x}, \omega)|$');
```

The last example shows the simulation result for sound-soft (Dirichlet) boundary conditions with $\sigma \to \infty$ ($Z=0$).


```python
# compute solution for very large sigma
u = Helmholtz_Robin(mesh, f, xs, sigma=dolfin.Constant(1e15))

# plot sound field
plot_soundfield(u)
plot_soundfield(abs(u))
plt.title(r'$|P(\mathbf{x}, \omega)|$');
```

The effect of the homogeneous Dirichlet boundary condition (zero pressure at the boundary) can be observed conveniently by inspecting the magnitude of the sound field close to the boundary.

**Copyright**

This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the IPython examples under the [MIT license](https://opensource.org/licenses/MIT).

755d3df79e9e8bc2e51ac6c9b50b6486effddf07
ipynb
Jupyter Notebook
FEM_Helmholtz_equation_Robin.ipynb
wgddd/computational_acoustics
6e1f249b562a73b3bb295981bea9d80d218a213f
MIT
2018-10-17T09:51:06.000Z
2022-02-09T02:09:48.000Z
FEM_Helmholtz_equation_Robin.ipynb
Zhengyu-Huang/computational_acoustics
6e1f249b562a73b3bb295981bea9d80d218a213f
MIT
FEM_Helmholtz_equation_Robin.ipynb
Zhengyu-Huang/computational_acoustics
6e1f249b562a73b3bb295981bea9d80d218a213f
MIT
2018-10-19T11:55:53.000Z
2022-02-10T15:24:20.000Z
Qwen/Qwen-72B
1. YES
2. YES