FDTD for Electric Fields

This article explains the concept behind the finite difference time domain technique, derives the equations and implements them in a Python script. The simulation result is shown in Figure 1.

Figure 1. E-field reflection at an air/glass interface

Maxwell’s Law tells us that a time-varying electric displacement field, D [C/m^2], induces a magnetic field strength, H [A/m]:

\[\begin{split}\frac{dD}{dt} &= \nabla × H \\\end{split}\]

In 1D, the curl operator reduces to a spatial derivative. The finite difference approximation is used to get the subsequent D value given the H value and the previous D value:

\[\begin{split}\frac{D\left(k,n+\frac{1}{2}\right)-D\left(k,n-\frac{1}{2}\right)}{\Delta t} &= \frac{H\left(k+\frac{1}{2},n\right)-H\left(k-\frac{1}{2},n\right)}{ℓ} \\ D\left(k,n+\frac{1}{2}\right) = \frac{\Delta t}{ℓ} &\left[ H\left(k+\frac{1}{2},n\right) - H\left(k-\frac{1}{2},n\right) \right] + D\left(k,n-\frac{1}{2}\right) \\ = \frac{1}{2c} &\left[ H\left(k+\frac{1}{2},n\right) - H\left(k-\frac{1}{2},n\right) \right] + D\left(k,n-\frac{1}{2}\right)\end{split}\]

where the sampling limit formula has been used, \(c \Delta t = ℓ/2\) (c being the speed of light in a vacuum). Note that the fields are interleaved; while D is aligned in space but misaligned in time, H is aligned in time but misaligned in space. A similar treatment is applied to the magnetic flux density, B [T], and the electric field, E [V/m]

\[\begin{split}\frac{dB}{dt} &= \nabla × E \\ \frac{B\left(k+\frac{1}{2},n+1\right)-B\left(k+\frac{1}{2},n\right)}{\Delta t} &= \frac{E\left(k+1,n+\frac{1}{2}\right)-E\left(k,n+\frac{1}{2}\right)}{ℓ} \\ B\left(k+\frac{1}{2},n+1\right) = \frac{\Delta t}{ℓ} &\left[ E\left(k+1,n+\frac{1}{2}\right)-E\left(k,n+\frac{1}{2}\right) \right] + B\left(k+\frac{1}{2},n\right) \\ = \frac{1}{2c} &\left[ E\left(k+1,n+\frac{1}{2}\right)-E\left(k,n+\frac{1}{2}\right) \right] + B\left(k+\frac{1}{2},n\right)\end{split}\]

FDTD uses a constitutive relation between two (or more) physical properties to iteratively solve for field data. At each time step, Property A will be a function of the previous state(s) of Property A and also of Property B. Similarly, Property B will be a function of its own previous state(s) and of Property A. With electric field simulations, these two properties are either the electric and magnetic fields, or the electric and displacement fields. This explanation uses the latter. The electric displacement field [C/m^2] equals the electrical permittivity [F/m] multiplied by the electric field [V/m]:

\[D = \epsilon E = \epsilon_r \epsilon_0 E\]

test

where \(\epsilon_r\) is the material-specific relative permittivity [unitless, 0 to 1] and \(\epsilon_0\) is the permittivity of free space. \(\epsilon_r\) is oftentimes a number from a lookup table, an estimate suitable for the frequency range of the simulation. We will use the Debye Relaxation Model,

\[\epsilon_{r_{Debye}} = \epsilon_{\infty} + \frac{\epsilon_{\Delta}}{1+j\omega\tau}\]

where \(\tau\) is the Debye relaxation time, \(\epsilon_{\infty}\) is the permittivity at the high frequency limit and \(\epsilon_{\Delta}\) is the DC permittivity minus the high-frequency value. These values are available in the literature for many different materials. The Debye Model is expanded to account for conductivity, \(\sigma\),

\[\epsilon_r = \epsilon_{\infty} + \frac{\epsilon_{\Delta}}{1+j\omega\tau} + \frac{\sigma}{j\omega\epsilon_0} = \epsilon_{\infty} + a(\omega) + b(\omega)\]

This expanded \(\epsilon_r\) serves as a useful general case for FDTD. If the Debye constants are not known, an operator can simply set \(\epsilon_{\Delta}\) to zero and \(\epsilon_{\infty}\) to an estimate suitable for the frequency range of the simulation. Likewise, the conductivity term may be skipped by using \(\sigma = 0\).

To go from the complex frequency domain to the time domain we use the Z-Transform, the discrete version of the Laplace Transform. The Z-Transform is often used in communications as a delay operator: if \(y(n) = x(n-1)\) then \(Y(z) = z^{-1}X(z)\). The following Z-Transform pair, taken from a lookup table, is needed:

\[\frac{1}{\alpha + j\omega} \leftrightarrow \frac{1}{1-z^{-1} e^{-\alpha \Delta t}}\]

\(a(\omega)\) is taken to the z-domain, multiplied by E(z) and rearranged to produce an displacement field contribution:

\[\begin{split}\frac{\epsilon_{\Delta}}{1+j\omega\tau} &= \frac{\epsilon_{\Delta} / \tau}{\frac{1}{\tau} + j \omega} \rightarrow A(z) = \left( \frac{\epsilon_{\Delta} / \tau}{1-z^{-1}e^{-\Delta t / \tau}} \right) \Delta t \\ A(z)E(z) &= E(z)\left(\frac{\epsilon_{\Delta}}{\tau}\right)\Delta t + z^{-1}A(z)e^{-\Delta t / \tau}\end{split}\]

\(b(\omega)\) is given a similar treatment:

\[\begin{split}\frac{\sigma}{j\omega\epsilon_0} &= \frac{\sigma / \epsilon_0}{0+j\omega} \rightarrow B(z) = \left( \frac{\sigma / \epsilon_0}{1-z^{-1}} \right) \Delta t \\ B(z)E(z) &= E(z)\left(\frac{\sigma}{\epsilon_0}\right)\Delta t + z^{-1}B(z)\end{split}\]

By substitution into the original displacement field equation,

\[\begin{split}\frac{D(z)}{\epsilon_0} &= \epsilon_{\infty}E(z) + A(z)E(z) + B(z)E(z) \\ &= \epsilon_{\infty}E(z) + E(z)\left(\frac{\epsilon_{\Delta}}{\tau}\right)\Delta t + z^{-1}A(z)e^{-\Delta t / \tau} + E(z)\left(\frac{\sigma}{\epsilon_0}\right)\Delta t + z^{-1}B(z)\end{split}\]

Rearranging to solve for E(z),

\[E(z) = \frac{D(z)/\epsilon_0 - z^{-1}A(z)e^{-\Delta t / \tau} - z^{-1}B(z)} {\epsilon_{\infty} + \left( \epsilon_{\Delta}/\tau + \sigma/\epsilon_0 \right) \Delta t}\]

We now have D(z) in terms of E(z) and previous values of itself. And, we have E(z) in terms of D(z) and previous values of itself. The criteria have been met for setting up an FDTD simulation. The equations are programmed into fdtd.py, below. This file contains the FDTD computations while ancillary functions, e.g. for plotting, are in util.py. The two files must be in the same directory before running fdtd.py.

import numpy as np
from util import scalarFields2Arrays, Medium, Source, render

EP0 = 8.8541878128e-12 # permittivity of free space, F/m
C0 = 299792458 # speed of light, m/s

def fdtd(Nts, dt, dz, mdms, srcs, ylim_E=[-1,1], ylim_H=[-1,1], dispNth=1):
    Ncells = np.sum([m.L for m in mdms])
    epR,sigma,chi1,t0 = scalarFields2Arrays(mdms,'epR','sigma','chi1','t0')
    axE,axH,pE,pH = None,None,None,None
    Ex,Hy,Dx,Ix,Sx = [np.zeros(Ncells) for _ in range(5)]
    Ex1_beforeLast, Ex1_last, ExN_beforeLast, ExN_last = 0,0,0,0
    
    for n in range(1, Nts+1): # for each time step
        t = n*dt
        
        # D from H
        for k in range(Ncells):
            Dx[k] += 0.5*(Hy[k-1] - Hy[k])
               
        # E from D
        for k in range(Ncells-1):
            Ex[k] = (Dx[k] - Ix[k] - np.exp(-dt/t0[k])*Sx[k])
            Ex[k] = Ex[k] / (epR[k] + sigma[k]*dt/EP0 + chi1[k]*dt/t0[k])
            Ix[k] = Ix[k] + sigma[k]*dt/EP0*Ex[k]
            Sx[k] = np.exp(-dt/t0[k])*Sx[k] + chi1[k]*dt/t0[k]*Ex[k]
        
        # Sources
        for src in srcs:
            i = int(round(src.z01*Ncells))
            Ex[i] += src.Ex(t)
            
        # Absorbing boundary condition
        if True:
            Ex[0] = Ex1_beforeLast
            Ex1_beforeLast = Ex1_last
            Ex1_last = Ex[1]
            
            Ex[Ncells-1] = ExN_beforeLast
            ExN_beforeLast = ExN_last
            ExN_last = Ex[Ncells-2]

        # H from E
        for k in range(Ncells-1):
            Hy[k] += 0.5*(Ex[k] - Ex[k+1])
        
        if not (n-1) % dispNth: # render on every dispNth step
            axE,axH,pE,pH = render(axE, axH, pE, pH, Ex, Hy, ylim_E, ylim_H, mdms)
#        plt.savefig('images/%d.png' % n)

if __name__ == '__main__':
    N = 160 # number of time steps
    dz = .01 # length of each cell
    dt = 2*dz/C0 # time step    
    ms = [Medium(50, epR=1, sigma=.002)] # air
    ms.append(Medium(25, epR=10, sigma=.001, chi1=2, t0=1e-9)) # glass    
    srcs = [Source(40*dt, 1/3, Source.gaussian, A=7, sigma=5*dt)]    
    fdtd(N, dt, dz, ms, srcs)