Hello FEniCS community!
I am trying to model a 2D acoustic plane wave with FEniCS. The domain is a rectangle whose left boundary is the source of the wave, and the other 3 boundaries should absorb the wave with no reflections.
After reading some tutorials and practicing with some examples, I am still unable to meet my requirements. So I thought I’d share my code with you in the hope that you can help me.
The wave equation
\partial u^2 / \partial t^2 = c^2 \nabla^2 u + f
The weak form
a(u,v) = \int_{\Omega} (vu^{n+1} - c^2 \Delta t^2 \nabla v \nabla u^{n+1})dx
L_{n+1}(v) = \int_{\Omega}(2u^{n} - u^{n-1} + \Delta t^2 f^{n+1})vdx - \int_{\Gamma^N} (gv)ds
Left boundary condition
u_D^L (t) = Asin(\omega t)
Other boundary conditions (I think the problem is here)
\partial u / \partial n = g = 0
The code
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
T = 10.0
num_steps = 250
dt = T / num_steps
tol = 1E-14
c = 2
# Create mesh and define function space
nx = 40
ny = 80
mesh = RectangleMesh(Point(-2,-1), Point(2,1),nx,ny)
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)
#Boundary expression
u_L = Expression('50*sin(10*t)', degree=2, t=0)
#Initial pressure
init = "0"
g_ini = Expression(init,t=0, degree=0)
u_0 = interpolate(g_ini,V)
u_1 = interpolate(g_ini,V)
#Boundary Dirichlet condition on the Left side
def boundary_L(x, on_boundary):
if on_boundary:
if near(x[0],-2,tol):
return True
else:
return False
else:
return False
bc_L = DirichletBC(V, u_L, boundary_L)
bc = [bc_L]
#Null source term
f = Constant(0)
#Neumann boundary condition
g = Expression('0', degree=1)
#Variational Form
a = u*v*dx - c*c*dt*dt*dot(grad(u), grad(v))*dx
L = (2*u_1 - u_0 + dt*dt*f)*v*dx - g*v*ds
# Time-stepping
u = Function(V)
t = 0
file = File("wave/wave.pvd")
for n in range(num_steps):
# Update current time
t += dt
u_L.t = t
# Compute solution
solve(a == L, u, bc)
# Plot solution
plot(u)
# Update previous solution
u_0.assign(u_1)
u_1.assign(u)
file<<u
This code produces the following result in the very first time step
which seems correct, but then the waves don’t propagate well. I have changed the parameters c, \omega, A,u_D and g, but it didn’t help at all.
I believe I should replace the boundary conditions g on \Gamma^N with some function to account for the wave absorption, but I can’t figure out how.
Can someone help with that?
Thanks a lot!