Helmholtz Equation

25. Helmholtz Equation#

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

\frac{\partial^{2}}{\partial t^{2}} u (x) e^{i ω t} - Δ u (x) e^{i ω t} = f (x) e^{i ω t} .

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

- Δ u - ω^{2} u = f .

We consider Dirichlet (hard) boundary conditions

u (x) = 0

and Robin (absorbing) boundary conditions

\frac{\partial u}{\partial n} - i ω u = 0 on Γ_{R}

The weak form is

\int_{Ω} \nabla u \nabla v - ω^{2} u v - i ω \int_{Γ_{R}} u v = \int_{Ω} 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);

25.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 ω 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);

Helmholtz Equation

Contents

25. Helmholtz Equation#

25.1. Incoming beam from the left side#