15. Computational Electromagnetics
15. Computational Electromagnetics¶
Learning Objectives¶
- Review the physical meaning of Maxwell's equations
- Understand numerical discretization of electromagnetic fields
- Introduction to the FDTD (Finite-Difference Time-Domain) method
- Understand the Yee lattice structure
- Learn the Courant condition for electromagnetic waves
1. Review of Maxwell's Equations¶
1.1 Differential Form¶
Maxwell's Equations (vacuum, SI units):
1. Gauss's Law (Electric):
βΒ·E = Ο/Ξ΅β
2. Gauss's Law (Magnetic):
βΒ·B = 0
3. Faraday's Law:
βΓE = -βB/βt
4. Ampère-Maxwell Law:
βΓB = ΞΌβJ + ΞΌβΞ΅β βE/βt
Where:
- E: Electric field [V/m]
- B: Magnetic field (magnetic flux density) [T]
- Ο: Charge density [C/mΒ³]
- J: Current density [A/mΒ²]
- Ξ΅β = 8.854Γ10β»ΒΉΒ² F/m (vacuum permittivity)
- ΞΌβ = 4ΟΓ10β»β· H/m (vacuum permeability)
- c = 1/β(ΞΌβΞ΅β) β 3Γ10βΈ m/s (speed of light)
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Physical constants
eps0 = 8.854e-12 # Vacuum permittivity [F/m]
mu0 = 4 * np.pi * 1e-7 # Vacuum permeability [H/m]
c0 = 1 / np.sqrt(mu0 * eps0) # Speed of light [m/s]
print(f"Vacuum permittivity Ξ΅β = {eps0:.3e} F/m")
print(f"Vacuum permeability ΞΌβ = {mu0:.3e} H/m")
print(f"Speed of light cβ = {c0:.3e} m/s")
def maxwell_equations_overview():
"""Maxwell's equations overview"""
print("=" * 60)
print("Maxwell's Equations and Physical Meaning")
print("=" * 60)
overview = """
βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
β Maxwell's Equations System β
βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ€
β β
β Gauss (Electric): βΒ·E = Ο/Ξ΅β β
β β Charges are sources of electric field divergence β
β β
β Gauss (Magnetic): βΒ·B = 0 β
β β No magnetic monopoles (always N-S pairs) β
β β
β Faraday: βΓE = -βB/βt β
β β Time-varying magnetic field induces electric curl β
β β
β AmpΓ¨re-Maxwell: βΓB = ΞΌβJ + ΞΌβΞ΅ββE/βt β
β β Current and time-varying E induce magnetic curl β
β β
βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
Derivation of wave equation:
Combining Faraday and Ampère laws:
βΓ(βΓE) = -β(βΓB)/βt = -ΞΌββJ/βt - ΞΌβΞ΅ββΒ²E/βtΒ²
Applying vector identity βΓ(βΓE) = β(βΒ·E) - βΒ²E:
In vacuum without charges/currents:
βΒ²E = ΞΌβΞ΅β βΒ²E/βtΒ² = (1/cΒ²) βΒ²E/βtΒ²
β Wave equation with velocity c = 1/β(ΞΌβΞ΅β)!
"""
print(overview)
maxwell_equations_overview()
1.2 1D Wave Equation¶
def electromagnetic_wave_1d():
"""Visualization of 1D electromagnetic wave analytical solution"""
# 1D TEM wave (propagating in x, y polarization)
# Only Ey and Hz components exist
# Spatial/temporal setup
L = 10.0 # Domain length [m]
T = 3 * L / c0 # Simulation time
x = np.linspace(0, L, 500)
times = [0, L/(3*c0), 2*L/(3*c0), L/c0]
# Initial Gaussian pulse
x0 = L / 4
sigma = L / 20
wavelength = L / 5
k = 2 * np.pi / wavelength
omega = c0 * k
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
for idx, t in enumerate(times):
ax = axes[idx // 2, idx % 2]
# Electric field Ey (Gaussian-modulated sine wave)
envelope = np.exp(-((x - x0 - c0*t) / sigma)**2)
Ey = envelope * np.sin(k*(x - x0) - omega*t)
# Magnetic field Hz (proportional to Ey, same phase)
Hz = Ey / (c0 * mu0) # Hz = Ey/(c*ΞΌβ) for plane wave
ax.plot(x, Ey, 'b-', linewidth=2, label=r'$E_y$ (Electric field)')
ax.plot(x, Hz * c0 * mu0, 'r--', linewidth=2, label=r'$\mu_0 c H_z$ (Magnetic field)')
ax.axhline(y=0, color='gray', linestyle='-', alpha=0.3)
ax.set_xlabel('x [m]')
ax.set_ylabel('Field amplitude')
ax.set_title(f't = {t*1e9:.2f} ns')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, L)
ax.set_ylim(-1.5, 1.5)
plt.suptitle('1D Electromagnetic Wave Propagation (TEM Mode)', fontsize=14)
plt.tight_layout()
plt.savefig('em_wave_1d.png', dpi=150, bbox_inches='tight')
plt.show()
# electromagnetic_wave_1d()
2. Discretization of Maxwell's Equations¶
2.1 Finite Difference Approach¶
1D TEM wave (Ey, Hz components):
Faraday's Law (z component):
βEy/βx = -βBz/βt = -ΞΌβ βHz/βt
Ampère's Law (y component):
βHz/βx = -Ξ΅β βEy/βt
Discretization (central difference):
βEy/βx β (Ey[i+1] - Ey[i]) / Ξx
βHz/βt β (Hz[n+1/2] - Hz[n-1/2]) / Ξt
Temporal staggering (Leapfrog):
- E at integer time steps: t = nΞt
- H at half-integer time steps: t = (n+1/2)Ξt
def fdtd_discretization_concept():
"""Visualization of FDTD discretization concept"""
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# (1) Spatial discretization
ax1 = axes[0]
# Staggered positions of E and H
n_points = 6
for i in range(n_points):
# E points (integer positions)
ax1.plot(i, 0.5, 'bo', markersize=15)
ax1.text(i, 0.7, f'$E_y^n[{i}]$', ha='center', fontsize=10, color='blue')
# H points (half-integer positions)
if i < n_points - 1:
ax1.plot(i + 0.5, 0.5, 'r^', markersize=12)
ax1.text(i + 0.5, 0.3, f'$H_z^{{n+1/2}}[{i}]$', ha='center', fontsize=9, color='red')
# Grid lines
for i in range(n_points):
ax1.axvline(x=i, color='gray', linestyle=':', alpha=0.5)
ax1.axhline(y=0.5, color='black', linestyle='-', alpha=0.3)
ax1.set_xlim(-0.5, n_points - 0.5)
ax1.set_ylim(0, 1)
ax1.set_xlabel('Position index i')
ax1.set_title('Spatial Staggering')
ax1.set_yticks([])
# Legend
ax1.plot([], [], 'bo', markersize=10, label='E field')
ax1.plot([], [], 'r^', markersize=10, label='H field')
ax1.legend(loc='upper right')
# (2) Temporal discretization (Leapfrog)
ax2 = axes[1]
n_steps = 5
for n in range(n_steps):
# E points (integer time)
ax2.plot(0.3, n, 'bo', markersize=15)
ax2.text(0.5, n, f'$E^{n}$', ha='left', fontsize=12, color='blue')
# H points (half-integer time)
ax2.plot(0.3, n + 0.5, 'r^', markersize=12)
ax2.text(0.5, n + 0.5, f'$H^{{{n}+1/2}}$', ha='left', fontsize=11, color='red')
# Arrows (update sequence)
for n in range(n_steps - 1):
# E -> H
ax2.annotate('', xy=(0.3, n + 0.5), xytext=(0.3, n),
arrowprops=dict(arrowstyle='->', color='green', lw=1.5))
# H -> E
ax2.annotate('', xy=(0.3, n + 1), xytext=(0.3, n + 0.5),
arrowprops=dict(arrowstyle='->', color='purple', lw=1.5))
ax2.set_xlim(0, 1)
ax2.set_ylim(-0.5, n_steps)
ax2.set_ylabel('Time step n')
ax2.set_title('Temporal Staggering (Leapfrog)')
ax2.set_xticks([])
plt.tight_layout()
plt.savefig('fdtd_discretization.png', dpi=150, bbox_inches='tight')
plt.show()
# fdtd_discretization_concept()
2.2 Update Equations¶
1D FDTD update equations:
1. H update (n β n+1/2):
Hz^(n+1/2)[i] = Hz^(n-1/2)[i] - (Ξt/ΞΌβΞx)(Ey^n[i+1] - Ey^n[i])
2. E update (n+1/2 β n+1):
Ey^(n+1)[i] = Ey^n[i] - (Ξt/Ξ΅βΞx)(Hz^(n+1/2)[i] - Hz^(n+1/2)[i-1])
Parameterization:
C_a = Ξt/(ΞΌβΞx) (H coefficient)
C_b = Ξt/(Ξ΅βΞx) (E coefficient)
Hz^(n+1/2)[i] = Hz^(n-1/2)[i] - C_a (Ey^n[i+1] - Ey^n[i])
Ey^(n+1)[i] = Ey^n[i] - C_b (Hz^(n+1/2)[i] - Hz^(n+1/2)[i-1])
def simple_1d_fdtd():
"""Simple 1D FDTD simulation"""
# Grid setup
Nx = 200
dx = 1e-3 # 1 mm
# Time setup
dt = dx / (2 * c0) # Satisfy Courant condition
n_steps = 500
# Initialize arrays
Ey = np.zeros(Nx)
Hz = np.zeros(Nx)
# Coefficients
Ca = dt / (mu0 * dx)
Cb = dt / (eps0 * dx)
# Courant number
S = c0 * dt / dx
print(f"Courant number S = {S:.4f}")
# Source position
source_pos = Nx // 4
# Recording
Ey_history = []
times_to_record = [0, 100, 200, 300, 400]
# Main loop
for n in range(n_steps):
# H update
Hz[:-1] = Hz[:-1] - Ca * (Ey[1:] - Ey[:-1])
# Source (Gaussian pulse)
t = n * dt
t0 = 30 * dt
tau = 10 * dt
source = np.exp(-((t - t0) / tau) ** 2)
Ey[source_pos] += source
# E update
Ey[1:] = Ey[1:] - Cb * (Hz[1:] - Hz[:-1])
# Boundary conditions (simple absorption)
Ey[0] = 0
Ey[-1] = 0
# Recording
if n in times_to_record:
Ey_history.append((n, Ey.copy()))
# Visualization
fig, axes = plt.subplots(2, 3, figsize=(15, 8))
x = np.arange(Nx) * dx * 1000 # mm
for idx, (n, Ey_snap) in enumerate(Ey_history):
ax = axes[idx // 3, idx % 3]
ax.plot(x, Ey_snap, 'b-', linewidth=1.5)
ax.axvline(x=source_pos * dx * 1000, color='red', linestyle='--', alpha=0.5)
ax.set_xlabel('x [mm]')
ax.set_ylabel('Ey')
ax.set_title(f'Step n = {n}')
ax.grid(True, alpha=0.3)
ax.set_ylim(-1.2, 1.2)
# Last subplot with description
ax = axes[1, 2]
ax.text(0.5, 0.8, 'FDTD Parameters:', fontsize=12, ha='center', transform=ax.transAxes)
ax.text(0.5, 0.6, f'Nx = {Nx}', fontsize=10, ha='center', transform=ax.transAxes)
ax.text(0.5, 0.5, f'dx = {dx*1000:.2f} mm', fontsize=10, ha='center', transform=ax.transAxes)
ax.text(0.5, 0.4, f'dt = {dt*1e12:.2f} ps', fontsize=10, ha='center', transform=ax.transAxes)
ax.text(0.5, 0.3, f'Courant S = {S:.2f}', fontsize=10, ha='center', transform=ax.transAxes)
ax.axis('off')
plt.suptitle('1D FDTD Simulation - Gaussian Pulse Propagation', fontsize=14)
plt.tight_layout()
plt.savefig('fdtd_1d_simple.png', dpi=150, bbox_inches='tight')
plt.show()
return Ey, Hz
# Ey, Hz = simple_1d_fdtd()
3. Yee Lattice¶
3.1 3D Yee Cell¶
Yee Lattice (1966):
- E and H components are spatially staggered
- Also staggered by half a step in time
- Natural discretization of Maxwell's equations
3D Yee cell structure:
- Ex: Edge centers of y-z faces
- Ey: Edge centers of x-z faces
- Ez: Edge centers of x-y faces
- Hx: Centers of x-perpendicular faces
- Hy: Centers of y-perpendicular faces
- Hz: Centers of z-perpendicular faces
Each E component is surrounded by 4 H components
Each H component is surrounded by 4 E components
def yee_cell_visualization():
"""3D Yee cell visualization"""
fig = plt.figure(figsize=(14, 6))
# (1) 3D view
ax1 = fig.add_subplot(121, projection='3d')
# Cube vertices
a = 1.0 # Cell size
# E field positions (edge centers)
# Ex: (a, a/2, 0), (a, a/2, a), (0, a/2, 0), (0, a/2, a)
Ex_pos = [(a, a/2, 0), (a, a/2, a), (0, a/2, 0), (0, a/2, a)]
# Ey positions
Ey_pos = [(a/2, a, 0), (a/2, a, a), (a/2, 0, 0), (a/2, 0, a)]
# Ez positions
Ez_pos = [(0, 0, a/2), (a, 0, a/2), (0, a, a/2), (a, a, a/2)]
# H field positions (face centers)
Hx_pos = [(a/2, a/2, 0), (a/2, a/2, a)]
Hy_pos = [(a/2, 0, a/2), (a/2, a, a/2)]
Hz_pos = [(0, a/2, a/2), (a, a/2, a/2)]
# Draw cube
vertices = [
[0, 0, 0], [a, 0, 0], [a, a, 0], [0, a, 0],
[0, 0, a], [a, 0, a], [a, a, a], [0, a, a]
]
# Edges
edges = [
[0, 1], [1, 2], [2, 3], [3, 0], # Bottom
[4, 5], [5, 6], [6, 7], [7, 4], # Top
[0, 4], [1, 5], [2, 6], [3, 7] # Vertical
]
for edge in edges:
points = [vertices[edge[0]], vertices[edge[1]]]
ax1.plot3D(*zip(*points), 'k-', alpha=0.3)
# E field (arrows)
for pos in Ex_pos:
ax1.quiver(pos[0]-0.15, pos[1], pos[2], 0.3, 0, 0, color='blue', arrow_length_ratio=0.3)
for pos in Ey_pos[:2]:
ax1.quiver(pos[0], pos[1]-0.15, pos[2], 0, 0.3, 0, color='blue', arrow_length_ratio=0.3)
for pos in Ez_pos[:2]:
ax1.quiver(pos[0], pos[1], pos[2]-0.15, 0, 0, 0.3, color='blue', arrow_length_ratio=0.3)
# H field (arrows)
for pos in Hx_pos:
ax1.quiver(pos[0]-0.15, pos[1], pos[2], 0.3, 0, 0, color='red', arrow_length_ratio=0.3)
for pos in Hy_pos:
ax1.quiver(pos[0], pos[1]-0.15, pos[2], 0, 0.3, 0, color='red', arrow_length_ratio=0.3)
for pos in Hz_pos:
ax1.quiver(pos[0], pos[1], pos[2]-0.15, 0, 0, 0.3, color='red', arrow_length_ratio=0.3)
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_zlabel('z')
ax1.set_title('3D Yee Cell')
# Legend
ax1.plot([], [], 'b-', linewidth=2, label='E field')
ax1.plot([], [], 'r-', linewidth=2, label='H field')
ax1.legend()
# (2) 2D view (x-y plane)
ax2 = fig.add_subplot(122)
# Grid
for i in range(3):
ax2.axhline(y=i, color='gray', linestyle='-', linewidth=0.5)
ax2.axvline(x=i, color='gray', linestyle='-', linewidth=0.5)
# Ez (cell vertices)
for i in range(3):
for j in range(3):
ax2.plot(i, j, 'bo', markersize=12)
ax2.text(i+0.1, j+0.1, f'Ez({i},{j})', fontsize=8, color='blue')
# Hx (horizontal edge centers)
for i in range(3):
for j in range(2):
ax2.plot(i, j+0.5, 'r>', markersize=10)
# Hy (vertical edge centers)
for i in range(2):
for j in range(3):
ax2.plot(i+0.5, j, 'r^', markersize=10)
# Hz (cell centers)
for i in range(2):
for j in range(2):
ax2.plot(i+0.5, j+0.5, 'rs', markersize=10)
ax2.text(i+0.6, j+0.5, f'Hz', fontsize=8, color='red')
ax2.set_xlabel('i')
ax2.set_ylabel('j')
ax2.set_title('2D Yee Grid (TM Mode)')
ax2.set_aspect('equal')
ax2.set_xlim(-0.3, 2.5)
ax2.set_ylim(-0.3, 2.5)
# Legend
ax2.plot([], [], 'bo', markersize=10, label='Ez')
ax2.plot([], [], 'r>', markersize=10, label='Hx')
ax2.plot([], [], 'r^', markersize=10, label='Hy')
ax2.plot([], [], 'rs', markersize=10, label='Hz')
ax2.legend(loc='upper right')
plt.tight_layout()
plt.savefig('yee_cell.png', dpi=150, bbox_inches='tight')
plt.show()
# yee_cell_visualization()
3.2 2D FDTD Modes¶
2D FDTD decomposition:
TM Mode (Transverse Magnetic):
- Components: Ez, Hx, Hy
- Ez in z direction, H in x-y plane
Equations:
βHx/βt = -(1/ΞΌ) βEz/βy
βHy/βt = (1/ΞΌ) βEz/βx
βEz/βt = (1/Ξ΅)(βHy/βx - βHx/βy) - Ο/Ξ΅ Ez
TE Mode (Transverse Electric):
- Components: Hz, Ex, Ey
- Hz in z direction, E in x-y plane
Equations:
βEx/βt = (1/Ξ΅) βHz/βy - Ο/Ξ΅ Ex
βEy/βt = -(1/Ξ΅) βHz/βx - Ο/Ξ΅ Ey
βHz/βt = (1/ΞΌ)(βEx/βy - βEy/βx)
def fdtd_2d_tm_mode():
"""2D FDTD TM mode simulation"""
# Grid setup
Nx, Ny = 100, 100
dx = dy = 1e-3 # 1 mm
# Time setup
dt = 1 / (c0 * np.sqrt(1/dx**2 + 1/dy**2)) * 0.99 # Courant
n_steps = 300
# Material properties
eps = eps0 * np.ones((Ny, Nx)) # Permittivity
mu = mu0 * np.ones((Ny, Nx)) # Permeability
sigma = np.zeros((Ny, Nx)) # Conductivity
# Initialize arrays
Ez = np.zeros((Ny, Nx))
Hx = np.zeros((Ny, Nx))
Hy = np.zeros((Ny, Nx))
# Coefficients
Ca = (1 - sigma * dt / (2 * eps)) / (1 + sigma * dt / (2 * eps))
Cb = (dt / eps) / (1 + sigma * dt / (2 * eps))
# Source position
source_x, source_y = Nx // 2, Ny // 2
# Courant number
S = c0 * dt * np.sqrt(1/dx**2 + 1/dy**2)
print(f"2D Courant number S = {S:.4f}")
# Snapshot recording
snapshots = []
record_steps = [50, 100, 150, 200, 250]
# Main loop
for n in range(n_steps):
# H update
Hx[:, :-1] = Hx[:, :-1] - dt / (mu[:, :-1] * dy) * (Ez[:, 1:] - Ez[:, :-1])
Hy[:-1, :] = Hy[:-1, :] + dt / (mu[:-1, :] * dx) * (Ez[1:, :] - Ez[:-1, :])
# Source (Gaussian pulse)
t = n * dt
t0 = 50 * dt
tau = 20 * dt
source = np.exp(-((t - t0) / tau) ** 2)
# E update
Ez[1:-1, 1:-1] = (Ca[1:-1, 1:-1] * Ez[1:-1, 1:-1] +
Cb[1:-1, 1:-1] * (
(Hy[1:-1, 1:-1] - Hy[:-2, 1:-1]) / dx -
(Hx[1:-1, 1:-1] - Hx[1:-1, :-2]) / dy
))
# Source injection
Ez[source_y, source_x] += source
# Snapshot recording
if n in record_steps:
snapshots.append((n, Ez.copy()))
# Visualization
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
x = np.arange(Nx) * dx * 1000
y = np.arange(Ny) * dy * 1000
X, Y = np.meshgrid(x, y)
for idx, (step, Ez_snap) in enumerate(snapshots):
ax = axes[idx // 3, idx % 3]
vmax = np.max(np.abs(Ez_snap)) * 0.8
if vmax == 0:
vmax = 1
im = ax.pcolormesh(X, Y, Ez_snap, cmap='RdBu_r', shading='auto',
vmin=-vmax, vmax=vmax)
ax.plot(source_x * dx * 1000, source_y * dy * 1000, 'k*', markersize=10)
ax.set_xlabel('x [mm]')
ax.set_ylabel('y [mm]')
ax.set_title(f'Step n = {step}')
ax.set_aspect('equal')
plt.colorbar(im, ax=ax, label='Ez')
# Last subplot
ax = axes[1, 2]
ax.text(0.5, 0.7, '2D FDTD TM Mode', fontsize=14, ha='center', transform=ax.transAxes)
ax.text(0.5, 0.5, f'Grid: {Nx} x {Ny}', fontsize=11, ha='center', transform=ax.transAxes)
ax.text(0.5, 0.4, f'dx = dy = {dx*1000:.1f} mm', fontsize=11, ha='center', transform=ax.transAxes)
ax.text(0.5, 0.3, f'Courant S = {S:.2f}', fontsize=11, ha='center', transform=ax.transAxes)
ax.axis('off')
plt.suptitle('2D FDTD Simulation - Circular Wave Propagation from Point Source', fontsize=14)
plt.tight_layout()
plt.savefig('fdtd_2d_tm.png', dpi=150, bbox_inches='tight')
plt.show()
return Ez, Hx, Hy
# Ez, Hx, Hy = fdtd_2d_tm_mode()
4. Courant Condition¶
4.1 Stability Condition for Electromagnetic Waves¶
Courant-Friedrichs-Lewy (CFL) Condition:
1D:
cΒ·Ξt/Ξx β€ 1
2D:
cΒ·ΞtΒ·β(1/ΞxΒ² + 1/ΞyΒ²) β€ 1
3D:
cΒ·ΞtΒ·β(1/ΞxΒ² + 1/ΞyΒ² + 1/ΞzΒ²) β€ 1
Physical meaning:
- Wave cannot travel more than one cell in one time step
- Numerical information propagation speed β₯ Physical wave speed
Isotropic grid (Ξx = Ξy = Ξz = Ξ):
1D: Ξt β€ Ξ/c
2D: Ξt β€ Ξ/(cβ2)
3D: Ξt β€ Ξ/(cβ3)
def courant_condition_analysis():
"""Courant condition analysis"""
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# (1) Stability vs Courant number
ax1 = axes[0]
# 1D FDTD simulation (various Courant numbers)
def run_1d_fdtd(S, n_steps=200):
Nx = 100
dx = 1e-3
dt = S * dx / c0
Ey = np.zeros(Nx)
Hz = np.zeros(Nx)
Ca = dt / (mu0 * dx)
Cb = dt / (eps0 * dx)
max_values = []
for n in range(n_steps):
Hz[:-1] = Hz[:-1] - Ca * (Ey[1:] - Ey[:-1])
if n == 0:
Ey[Nx//4] = 1.0 # Initial pulse
Ey[1:] = Ey[1:] - Cb * (Hz[1:] - Hz[:-1])
Ey[0] = Ey[-1] = 0
max_values.append(np.max(np.abs(Ey)))
return max_values
courant_numbers = [0.5, 0.9, 1.0, 1.01, 1.1]
colors = ['green', 'blue', 'orange', 'red', 'darkred']
for S, color in zip(courant_numbers, colors):
max_vals = run_1d_fdtd(S)
label = f'S = {S}' + (' (stable)' if S <= 1 else ' (unstable)')
ax1.semilogy(max_vals, color=color, linewidth=1.5, label=label)
ax1.axhline(y=1, color='gray', linestyle='--', alpha=0.5)
ax1.set_xlabel('Time step')
ax1.set_ylabel('max|Ey|')
ax1.set_title('1D FDTD: Stability vs Courant Number')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_ylim(1e-3, 1e10)
# (2) Dispersion relation
ax2 = axes[1]
# FDTD numerical dispersion relation
# Ο_num = (2/Ξt) arcsin(S sin(kΞx/2))
# Analytical: Ο = ck
k_norm = np.linspace(0, np.pi, 100) # kΞx
for S in [0.5, 0.8, 1.0]:
omega_exact = k_norm / S # Normalized ΟΞt
omega_fdtd = 2 * np.arcsin(S * np.sin(k_norm / 2))
# Phase velocity ratio
vp_ratio = omega_fdtd / omega_exact
ax2.plot(k_norm / np.pi, vp_ratio, linewidth=2, label=f'S = {S}')
ax2.axhline(y=1, color='black', linestyle='--', alpha=0.5, label='Exact')
ax2.set_xlabel(r'$k\Delta x / \pi$')
ax2.set_ylabel(r'$v_p^{num} / c$')
ax2.set_title('FDTD Numerical Dispersion (Phase Velocity Ratio)')
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0, 1)
ax2.set_ylim(0.8, 1.05)
plt.tight_layout()
plt.savefig('courant_condition.png', dpi=150, bbox_inches='tight')
plt.show()
# courant_condition_analysis()
4.2 Numerical Dispersion¶
FDTD numerical dispersion:
1D dispersion relation:
sin(ΟΞt/2) = SΒ·sin(kΞx/2)
Where S = cΞt/Ξx (Courant number)
Analytical solution: Ο = ck (no dispersion)
FDTD: Ο β ck (numerical dispersion occurs)
Problems:
- Severe dispersion at short wavelengths (high frequency)
- Waveform distortion, phase error
Solutions:
- Use minimum 10-20 cells per wavelength
- Use higher-order difference schemes
- Apply dispersion correction techniques
5. Material Modeling¶
5.1 Dielectrics and Conductors¶
Considering material properties:
Dielectrics (Ξ΅ > Ξ΅β):
- Wave velocity decreases: v = c/β(Ξ΅α΅£)
- Wavelength decreases: Ξ» = Ξ»β/β(Ξ΅α΅£)
- Grid resolution adjustment needed
Conductors (Ο > 0):
- Current induction: J = ΟE
- Wave attenuation
- Skin depth: Ξ΄ = β(2/(ΟΞΌΟ))
Lossy media:
Add -ΟE/Ξ΅ to βE/βt term
def material_modeling_fdtd():
"""2D FDTD with material properties"""
# Grid setup
Nx, Ny = 150, 100
dx = dy = 1e-3
# Time setup
dt = 1 / (c0 * np.sqrt(1/dx**2 + 1/dy**2)) * 0.99
n_steps = 400
# Initialize material property arrays
eps_r = np.ones((Ny, Nx)) # Relative permittivity
sigma = np.zeros((Ny, Nx)) # Conductivity
# Add dielectric block (Ξ΅α΅£ = 4)
eps_r[30:70, 80:100] = 4.0
# Add conductor block (PEC approximation)
sigma[40:60, 40:55] = 1e7 # High conductivity
# Actual permittivity
eps = eps0 * eps_r
# Initialize arrays
Ez = np.zeros((Ny, Nx))
Hx = np.zeros((Ny, Nx))
Hy = np.zeros((Ny, Nx))
# Coefficients (including loss)
Ca = (1 - sigma * dt / (2 * eps)) / (1 + sigma * dt / (2 * eps))
Cb = (dt / eps) / (1 + sigma * dt / (2 * eps))
# Source position
source_x, source_y = 20, Ny // 2
# Snapshots
snapshots = []
record_steps = [50, 100, 200, 300]
for n in range(n_steps):
# H update
Hx[:, :-1] = Hx[:, :-1] - dt / (mu0 * dy) * (Ez[:, 1:] - Ez[:, :-1])
Hy[:-1, :] = Hy[:-1, :] + dt / (mu0 * dx) * (Ez[1:, :] - Ez[:-1, :])
# Source
t = n * dt
t0 = 50 * dt
tau = 15 * dt
freq = 5e9 # 5 GHz
source = np.exp(-((t - t0) / tau) ** 2) * np.sin(2 * np.pi * freq * t)
# E update
Ez[1:-1, 1:-1] = (Ca[1:-1, 1:-1] * Ez[1:-1, 1:-1] +
Cb[1:-1, 1:-1] * (
(Hy[1:-1, 1:-1] - Hy[:-2, 1:-1]) / dx -
(Hx[1:-1, 1:-1] - Hx[1:-1, :-2]) / dy
))
Ez[source_y, source_x] += source
if n in record_steps:
snapshots.append((n, Ez.copy()))
# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
x = np.arange(Nx) * dx * 1000
y = np.arange(Ny) * dy * 1000
X, Y = np.meshgrid(x, y)
for idx, (step, Ez_snap) in enumerate(snapshots):
ax = axes[idx // 2, idx % 2]
vmax = np.max(np.abs(Ez_snap)) * 0.5
if vmax == 0:
vmax = 1
im = ax.pcolormesh(X, Y, Ez_snap, cmap='RdBu_r', shading='auto',
vmin=-vmax, vmax=vmax)
# Display material regions
ax.contour(X, Y, eps_r, levels=[2], colors='green', linewidths=2)
ax.contour(X, Y, sigma, levels=[1e6], colors='black', linewidths=2)
ax.plot(source_x * dx * 1000, source_y * dy * 1000, 'r*', markersize=10)
ax.set_xlabel('x [mm]')
ax.set_ylabel('y [mm]')
ax.set_title(f'Step n = {step}')
ax.set_aspect('equal')
plt.colorbar(im, ax=ax, label='Ez')
# Legend
axes[0, 0].plot([], [], 'g-', linewidth=2, label=r'Dielectric ($\epsilon_r=4$)')
axes[0, 0].plot([], [], 'k-', linewidth=2, label='Conductor (PEC)')
axes[0, 0].legend(loc='upper right')
plt.suptitle('2D FDTD with Material Properties', fontsize=14)
plt.tight_layout()
plt.savefig('fdtd_materials.png', dpi=150, bbox_inches='tight')
plt.show()
return Ez
# Ez = material_modeling_fdtd()
6. CEM Method Comparison¶
6.1 Major Numerical Techniques¶
Computational Electromagnetics (CEM) methods:
1. FDTD (Finite-Difference Time-Domain):
- Direct time-domain analysis
- Broadband characteristics in one simulation
- Simple implementation, easy parallelization
2. FEM (Finite Element Method):
- Advantageous for complex geometries
- Uses unstructured grids
- Higher accuracy with higher-order elements
3. MoM (Method of Moments):
- Based on integral equations
- Suitable for open-region problems
- Antenna, scattering problems
4. FIT (Finite Integration Technique):
- Integral form Maxwell's equations
- Excellent energy conservation
- Commercial software (CST)
5. FDFD (Finite-Difference Frequency-Domain):
- Frequency-domain steady state
- Efficient for specific frequency analysis
def cem_methods_comparison():
"""CEM methods comparison"""
print("=" * 70)
print("Computational Electromagnetics (CEM) Method Comparison")
print("=" * 70)
comparison = """
βββββββββββββββ¬ββββββββββββββ¬ββββββββββββββ¬ββββββββββββββ¬ββββββββββββββ
β Method β FDTD β FEM β MoM β FDFD β
βββββββββββββββΌββββββββββββββΌββββββββββββββΌββββββββββββββΌββββββββββββββ€
β Domain β Time β Time/Freq β Frequency β Frequency β
β Grid β Structured β Unstructuredβ Surface β Structured β
β Matrix β None β Sparse β Dense β Sparse β
β Broadband β Efficient β Multi-calc β Multi-calc β Multi-calc β
β Inhomog. β Easy β Easy β Difficult β Easy β
β Open region β ABC/PML β Inf. elem. β Automatic β ABC β
β Parallel β Very easy β Easy β Difficult β Easy β
β Nonlinear β Possible β Possible β Difficult β Difficult β
βββββββββββββββ΄ββββββββββββββ΄ββββββββββββββ΄ββββββββββββββ΄ββββββββββββββ
Application areas:
FDTD:
- Optical devices, antennas, EMC/EMI
- Bioelectromagnetics
- Radar cross section (RCS)
FEM:
- Waveguides with complex geometries
- Microwave circuits
- Antenna feed structures
MoM:
- Wire antennas
- Planar microstrip
- Scattering problems
Commercial software:
- FDTD: Lumerical, XFdtd
- FEM: ANSYS HFSS, COMSOL
- MoM: FEKO, NEC
- FIT: CST Studio
"""
print(comparison)
cem_methods_comparison()
7. Exercises¶
Exercise 1: Maxwell's Equations¶
Derive the wave equation for the magnetic field B by combining Faraday's law and Ampère's law.
Exercise 2: 1D FDTD¶
Simulate reflection and transmission at a material interface (Ξ΅β β Ξ΅β) in 1D FDTD code. Compare the reflection and transmission coefficients with Fresnel formulas.
Exercise 3: Courant Condition¶
Compare the numerical dispersion for Courant numbers S = 0.5 and S = 1.0 in a 2D isotropic grid. Which case is more accurate?
Exercise 4: Yee Lattice¶
Draw the positions of H components needed for updating Ex in a 3D Yee lattice and write the update equation.
8. References¶
Key Papers¶
- Yee (1966) "Numerical Solution of Initial Boundary Value Problems Involving Maxwell's Equations in Isotropic Media" - Original FDTD paper
- Taflove & Brodwin (1975) - Absorbing boundary conditions
Textbooks¶
- Taflove & Hagness, "Computational Electrodynamics: The Finite-Difference Time-Domain Method"
- Sullivan, "Electromagnetic Simulation Using the FDTD Method"
- Jin, "The Finite Element Method in Electromagnetics"
Open Source Tools¶
- MEEP (MIT, FDTD)
- gprMax (Ground Penetrating Radar)
- OpenEMS (FDTD + circuits)
Summary¶
Computational Electromagnetics Essentials:
1. Maxwell's Equations:
- βΓE = -βB/βt (Faraday)
- βΓH = J + βD/βt (AmpΓ¨re)
- βΒ·D = Ο (Gauss electric)
- βΒ·B = 0 (Gauss magnetic)
2. FDTD Essentials:
- E and H spatially/temporally staggered
- Yee lattice structure
- Leapfrog time advancement
3. Courant Condition:
1D: cΞt/Ξx β€ 1
2D: cΞtβ(1/ΞxΒ² + 1/ΞyΒ²) β€ 1
3D: cΞtβ(1/ΞxΒ² + 1/ΞyΒ² + 1/ΞzΒ²) β€ 1
4. Numerical Dispersion:
- Severe at short wavelengths
- Recommended 10-20 cells per wavelength
5. Material Modeling:
- Dielectrics: Ξ΅ = Ξ΅βΞ΅α΅£
- Conductors: Ο > 0
- Lossy media: Modified Ca, Cb coefficients
In the next lesson, we will cover detailed FDTD implementation and absorbing boundary conditions.