Waveguides

17.5. Waveguides#

Waveguides are one-dimensional structures to propagate waves. Typical examples are optical fibres. The simple idea is to have a fibre with refrection index larger than one, such that the wave stays inside the fibre due to total reflection.

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw
from time import sleep

omega = 50
epsr = 5
thick = 0.04

rect = Rectangle(1,1).Face()
rect.edges.name="outer"

f = HalfSpace(Pnt(0,0,0), Vec(0,0,1)).faces[0]
curve = MoveTo(-0.1,0.8).Line(0.7).Arc(0.2,-90).Line(0.5) \
    .Finish().Wire()
thickcurve = curve.Offset(f, thick/2, 'arc', False)

strip = Face(Wire(thickcurve.edges))
board = rect-strip
strip = rect*strip

strip.faces.col=[1,0,0]
board.faces.col=[0,0,1]
strip.faces.name="strip"
board.faces.name="board"

shape = Glue([strip, board])
Draw (shape);

mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.01)).Curve(3)
Draw (mesh);

17.5.1. waveguide modes#

We introduce local coordinates \((\xi, \eta)\) in longitudinal and transveral direction of a straight waveguide of thickness \(t\). We assume that

\[\begin{split} n^2(\xi, \eta) = \left\{ \begin{array}{cl} \varepsilon_r & | \eta | \leq \frac{t}{2} \\ 1 & \text{else} \end{array} \right. \end{split}\]

Look for solutions

\[ u(\xi, \eta, t) = u_y(\eta) e^{i k_\xi \xi - i \omega t} \]

In each region \(\eta < -\frac{t}{2}\), \(-\frac{t}{2} < \eta < \frac{t}{2}\), and \(\eta > \frac{t}{2}\), the cross-section function \(u_y(\eta) = e^{i k_\eta \eta}\), with a real or imaginary \(k_\eta\).

Inside the waveguide we look for a

\[ k_\xi^2 + k_{\eta,1}^2 = n^2 \omega^2, \]

and in the air domain we want an exponentially decaying solution due to an imaginary \(k\), i.e.

\[ k_\xi^2 + (i k_{\eta,2})^2 = \omega^2 \]
\[\begin{split} u_y = \left\{ \begin{array}{ll} \cos(k_{\eta,1} \eta) & \text{for } | \eta | \leq \tfrac{t}{2} \\ \cos(k_{\eta,1} \tfrac{t}{2} ) \exp ( k_{\eta,2} \big(|\eta|-\tfrac{t}{2} ) \big) & \text{else} \end{array} \right. \end{split}\]

This \(u_y\) is continuous. Demanding \(C^1\) continuity leads to the equation

\[ k_{\eta,1} \tan (k_{\eta,1} \tfrac{t}{2}) = k_{\eta,2} \]

These are now 3 equations for the components of the wave vectors \(k\).

import sympy as sym

ky1 = sym.symbols('ky1')
ky1 = float(sym.nsolve( ky1**2 * (1+sym.tan(ky1*thick/2)**2) - (epsr-1)*omega**2, ky1, sqrt(epsr-1)*omega))

kx = sqrt( epsr*omega**2 - ky1**2 )
ky2 = sqrt (kx**2-omega**2)

print ("kx = ", kx, "ky1 =", ky1, ", ky2 = ", ky2)

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[4], line 1
----> 1 import sympy as sym
      3 ky1 = sym.symbols('ky1')
      4 ky1 = float(sym.nsolve( ky1**2 * (1+sym.tan(ky1*thick/2)**2) - (epsr-1)*omega**2, ky1, sqrt(epsr-1)*omega))

ModuleNotFoundError: No module named 'sympy'

def uy(eta):
    eta = IfPos(eta, eta, -eta)
    return IfPos( thick/2-eta, 
                   cos(ky1 * eta),
                   cos(ky1 * thick/2) * exp( -ky2 * (eta-thick/2) ))

order=2   # order=1 or order=2 is supported

tau = 5e-4
tend = 2

u0 = uy(eta=y-0.8) * sin(kx*x) * exp(-10**2 * (x-0.3)**2)

v0 = 0

fes = H1LumpingFESpace(mesh, order=order)  
u,v = fes.TnT()

cfepsr = mesh.MaterialCF( {"strip" : 5, "board" : 1} )
mform = cfepsr*u*v*dx(intrules=fes.GetIntegrationRules())
aform = grad(u)*grad(v)*dx
rform = sqrt(BoundaryFromVolumeCF(epsr))*u*v*ds("outer", intrules=fes.GetIntegrationRules())

m = BilinearForm(mform, diagonal=True).Assemble()
a = BilinearForm(aform).Assemble()
r = BilinearForm(rform, diagonal=True, check_unused=False).Assemble()

mstar = BilinearForm(mform+tau/2*rform, diagonal=True).Assemble()
mstarinv = mstar.mat.Inverse(fes.FreeDofs())  

gfu = GridFunction(fes)
gfu.Set(u0)

scene = Draw(gfu, order=2, deformation=True, scale=1, euler_angles=[-60,0,-40])

unew = gfu.vec.CreateVector()
uold = gfu.vec.CreateVector()
uold.data = gfu.vec

rhs = gfu.vec.CreateVector()

with TaskManager(): 
    for n in range(int(tend/tau)):
        rhs.data = 2*m.mat*gfu.vec - m.mat*uold - tau**2*a.mat*gfu.vec + tau/2*r.mat*uold
        unew.data = mstarinv*rhs.data
        uold.data = gfu.vec
        gfu.vec.data = unew.data
        if n % 10 == 0:
            scene.Redraw()

scene.Redraw()

17.5.2. Exercises:#

  • extend the waveguide

  • try crossings and junctions

  • let a second waveguide come close, the wave packet should partially tunnel over

Contents