This notebook illustrates the numerical solution of the wave equation for an harmonic excitation and Robin boundary conditions using the 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.
The inhomogeneous 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$.
Starting from the variational formulation of the Helmholtz equation (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.
The numerical solution of the variational problem is based on FEniCS, 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.
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);
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.
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.
# 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.
# 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$).
# 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. Feel free to use the notebook for your own purposes. The text is licensed under Creative Commons Attribution 4.0, the code of the IPython examples under the MIT license.