Incident wave boundary condition for acoustic scattering problem

1

Hi everyone,

I am very new to Fenics and I’ve very slow learning curve at the moment. I am trying to solve an acoustic scattering problem but I’ve not been able to implement Dirichlet boundary condition for my incident plane wave. I’ve tried understanding and merging multiple codes and tutorials available on acoustic scattering and multiple BC problems. What is wrong with my code?

Define domain with cylinder and mesh

channel = Rectangle(Point(-0.75, -0.8), Point(0.75, 0.4))
cylinder = Circle(Point(0, 0), 0.05, 32)
domain = channel - cylinder
mesh = generate_mesh(domain, 64)
plot(mesh)

def FEM_Helmholtz(mesh, frequency, inc, c):
‘’‘numerical solution of the Helmholtz equation using the FEM’’’

# squared wavenumber
k2 = dolfin.Constant(2*dolfin.pi*frequency/c)**2

# define function space
V = dolfin.VectorFunctionSpace(mesh, "CG", 1, dim=2)

# define outer boundary conditions

tol = 1E-14

def boundary_west(x, on_boundary):
    return on_boundary and near(x[0],-0.75, tol) 
bc_w = dolfin.DirichletBC(V, dolfin.Constant(0), boundary_west)


def boundary_east(x, on_boundary):
    return on_boundary and near(x[0],0.75, tol)
bc_e = dolfin.DirichletBC(V, dolfin.Constant(0), boundary_east)


def boundary_south(x, on_boundary):
    return on_boundary and near(x[1],-0.8, tol)
bc_s = dolfin.DirichletBC(V, dolfin.Constant(0), boundary_south)


def boundary_north(x, on_boundary):
    return on_boundary and near(x[1],0.4, tol)
bc_n = dolfin.DirichletBC(V, inc, boundary_north, 0)


def boundary_cylinder(x, on_boundary):
    return on_boundary and near(x[0],-0.1, tol) and near(x[0],0.1, tol) and near(x[1],-0.1, tol) and near(x[1],0.1, tol)
bc_c = dolfin.DirichletBC(V, dolfin.Constant(0), boundary_cylinder, 0)

bcs = [bc_w, bc_e, bc_n, bc_s, bc_c]


# 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
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
L_r = dolfin.Constant(0) * v_r * dolfin.dx
L_i = dolfin.Constant(0) * v_i * dolfin.dx
  
# assemble the system
a = a_r + a_i
L = L_r + L_i
A, b = dolfin.assemble_system(a, L, bcs)

# 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):
‘’‘plot 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 / m’)
plt.ylabel(r’y / m’)
plt.colorbar(fig, fraction=0.038, pad=0.04);

compute solution

u = FEM_Helmholtz(mesh, 1000, dolfin.Constant((0, -10)), 343)

plot sound field

plot_soundfield(u)
plot_soundfield(abs(u))
plt.title(r’|P(\mathbf{x}, \omega)|’);