Helmholtz Scattering Problem

1

Hello!

I’m trying to deal with wave scattering in frequency domain using FEniCS.

I need to solve the exterior Helmholtz problem (p is the deviation from ambient pressure and k the wave number):
\nabla ^2 p + k^2 p = 0 on \Omega
\partial p / \partial \vec{n} = g on \Gamma_{obstacle}
p = 0 on \Gamma_{Dirichlet}

using the Weak Formulation (q is a test funcion):
\int_{\Omega} (\nabla p \cdot \nabla q - k^2 pq)d \Omega = \int_{\Gamma} (gq)ds

For that end, I’d like to address two issues

  1. How to apply the boundary condition on \Gamma_{obstacle}?
    At first, the obstacle may be a circle of radius r in the middle of the domain, but at some point I’d like to use more complicated geometries. I have seen the Subdomains and boundary conditions tutorial, but I couldn’t make much sense of it.

  2. How to implement an incident plane wave moving in x+ direction?
    I’m not sure if I have to account for an internal source f or if this is done with boundary conditions.

I’ve come this far with the code (I’m a newbie so any comments will help me a lot):

from fenics import *
import matplotlib.pyplot as plt

nx = 80
ny = 80
k = 1

#mesh and function space
mesh = RectangleMesh(Point(-1, -1), Point(1, 1), nx, ny)
V = FunctionSpace(mesh, "Lagrange", 1)

#boundary conditions
p0 = Constant(0.0)
bc = DirichletBC(V,p0, "on_boundary")

#variational problem
p = TrialFunction(V)
q = TestFunction(V)
#g = ... Issue 1
#f = ... Issue 2
a = inner(grad(p),grad(q))*dx - k*k*p*q*dx
L = g*q*ds

#compute solution
p = Function(V)
solve(a == L, p, bc)

plot(p)
plt.show()

I’d be glad if someone could show me how to proceed! :slightly_smiling_face: