24. Helmholtz Equation

Inserting a time harmonic right hand side \(f(x,t) = f(x) e^{i \omega t}\) into the wave equation leads to a time-harmonic solution solution \(u(x,t) = u(x) e^{i \omega t}\) satisfying

\[ \frac{\partial^2}{\partial t^2} u(x) e^{i \omega t} - \Delta u(x) e^{i \omega t} = f(x) e^{i \omega t}. \]

This leads to the Helmholtz equation (also known as frequency domain wave equation)

\[ -\Delta u - \omega^2 u = f. \]

We consider Dirichlet (hard) boundary conditions

\[ u(x) = 0 \]

and Robin (absorbing) boundary conditions

\[ \frac{\partial u}{\partial n} - i \omega u = 0 \qquad \text{ on } \Gamma_R \]

The weak form is

\[ \int_\Omega \nabla u \nabla v - \omega^2 u v - i \omega \int_{\Gamma_R} u v = \int_\Omega f v \]

from ngsolve import *
from ngsolve.webgui import Draw

mesh = Mesh(unit_square.GenerateMesh(maxh=0.05))

fes = H1(mesh, order=5, complex=True)
u,v = fes.TnT()

omega = 2*pi*10

a = BilinearForm(fes)
a += (grad(u)*grad(v)-omega**2*u*v)*dx
a += -1j*omega*u*v*ds
a.Assemble()

source = exp(-50**2*((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)))
f = LinearForm(source*v*dx(bonus_intorder=5)).Assemble()

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
Draw (1e3*gfu, mesh, order=3, animate_complex=True, deformation=True);

24.1. Incoming beam from the left side

we model an incoming beam (e.g. laser beam) from the left boundary by means of the boundary condition

\[ \frac{\partial u}{\partial n} - i \omega u = g \]

from netgen.occ import *
shape = MoveTo(-1,-1).Rectangle (2,2).Face()
shape.edges.Min(X).name="left"
shape.edges.Min(Y).name="bottom"
shape.edges.Max(X).name="right"
shape.edges.Max(Y).name="top"
circ = Circle ( (0.6, 0.0), 0.07).Face()
rect = MoveTo(-0.2,-0.4).Rotate(45).Rectangle(1,0.1).Face()
# shape = shape-circ
# shape = shape-rect
mesh = Mesh (OCCGeometry( shape, dim=2).GenerateMesh(maxh=0.05)).Curve(3)
Draw (shape);

fes = H1(mesh, order=5, complex=True)
u,v = fes.TnT()

a = BilinearForm(fes)
f = LinearForm(fes)

omega = 20*pi
outer = "left|right|bottom|top"
a += grad(u)*grad(v)*dx-omega**2*u*v*dx - 1j*omega*u*v*ds(outer)
f = LinearForm(fes)

incoming = 100*exp(-5**2*y*y)
f += incoming*v * ds("left")

a.Assemble()
f.Assemble()

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
Draw (gfu, order=3, animate_complex=True);