FDTD Method and Discretization

• Numerical Transmission Line Solutions will be explored in this module

The Finite-Difference Time-Domain Method (FDTD)

• FDTD is a popular computational electromagnetics method.
• It directly and intuitively solves Maxwell's equations.
• It is powerful and flexible, excelling in simulations of transients and broadband responses.
• Drawbacks include long computation times for resonant interactions and the need for large computational domains.
• Matlab will be used to program a one-dimensional FDTD algorithm to solve for voltages and currents along a transmission line.
• The same code will then be applied to electric and magnetic fields of plane waves propagating through space.

Discretization and Taylor Expansion

• Voltages and currents along a transmission line can be represented in Matlab by sampling them at discrete points.
• Telegrapher's equations relate these quantities through derivatives.
• A centered finite difference equation approximates the derivative of a discrete function based on the Taylor expansion:
  • f(z + Δz/2) - f(z - Δz/2) / Δz = df(z)/dz + O(Δz^2).
• This approximation has an error term of O(Δz^2), which approaches zero as Δz decreases.

Discretizing Telegrapher's Equations

• Telegrapher's equations are discretized using centered finite differences.
• The equations are evaluated at a point (z, t) using spatial discretization Δz and temporal discretization Δt.
• This discretization leads to second-order accuracy in Δz and Δt.
• The discretized equations are:
  • V(z + Δz/2, t) - V(z - Δz/2, t) / Δz = -L * (I(z, t + Δt/2) - I(z, t - Δt/2) / Δt)
  • I(z + Δz/2, t) - I(z - Δz/2, t) / Δz = -C * (V(z, t + Δt/2) - V(z, t - Δt/2) / Δt)
• The voltage and current values in each equation are offset in space by Δz/2 and in time by Δt/2.

Yee's Algorithm and Staggered Grids

• Yee proposed an algorithm in 1966 for solving Maxwell's equations using staggered electric and magnetic field lattices.
• This algorithm solves the fields by leap-frogging through time.
• A similar concept can be applied to solve for voltages and currents on a transmission line using grids staggered in space and time.
• Updating voltage and current on every half time step provides a natural and efficient solution to the Telegrapher's equations.
• Temporal finite differences contain field points located half a time step into the future.
• The equations are the FDTD method update equations.

Update Equations and Grid Placement

• Voltages are placed at integer space and time steps, while currents are placed at half-integer space and time steps.
• The indices k and n represent space and time, respectively.
• Update equations:
  • V_k^(n+1) = V_k^n - (Δt / (C * Δz)) * (I_(k+1/2)^(n+1/2) - I_(k-1/2)^(n+1/2))
  • I_(k+1/2)^(n+1/2) = I_(k+1/2)^(n-1/2) - (Δt / (L * Δz)) * (V_(k+1)^n - V_k^n)

Matlab Implementation

• To implement the 1D transmission line in Matlab, specify the per-unit-length inductance L and capacitance C of the line first; for example, L = 166.8 nH/m and C = 66.71 pF/m.
• The length of the line is discretized into N unit cells each with length Δz; for example, a 2.5 m line discretized into N = 250 unit cells, each of length Δz = 1 cm.
• Increasing the number of unit cells allows you to capture finer spatial variations in voltages and currents.
• When frequencies go up, smaller space steps are needed to capture smaller wavelengths.
• A good rule of thumb is to discretize the transmission line using at least ten unit cells per wavelength, calculated using the wavelength of the highest frequency signal on the line.

Time Step and Courant Condition

• The time step needs to be specified.
• Values in each unit cell are updated based on the values of the adjacent unit cells.
• The Courant condition is used to select the time step because it is the limit of stability.
  • Δt ≤ Δz / c
• c is the velocity in the Courant condition, it is the velocity of the transmission line (i.e. the velocity of the signals on the line) and depends on the transmission-line parameters. For example, when L and C are increased, this velocity will decrease.
• If instabilities are seen in the code, check that the time step is sufficiently small.
• For example, signals propagate at the speed of light in vacuum; given this, the maximum time step for the code is `Δt ≤ 0.01 m / (3 * 10^8 m/s) = 33 ps``

Programming the FDTD Kernel

• Voltages and currents need to be represented as vectors to program the FDTD kernel.
• Step through time, executing an update equation every half time step.
• Spatial offsets between voltage and current must be implicit in the code.
• Matlab will only accept integers when indexing arrays, so notes must be kept of which indices represent half steps and which represent whole steps.
• If the transmission line is discretized into N unit cells, the update equations are applied to the voltages over unit cells 1 to N-1 and to the currents over unit cells 2 to N.
• If V[N] and I[1] are not updated, they will remain constant over time.
• Signals on the transmission line can be visualized by plotting voltages and currents after each update.
• Update voltage data using: set(h, 'YData', V)
• Include a drawnow command after set to make sure Matlab updates the figure window right away.
• You need to supply energy to t-line (i.e. add a source!)

Adding a Source

• A Gaussian pulse can be programmed using the equation f(t) = e^(-(t-t0)^2 / τ^2).
  • t0 is the pulse center.
  • τ = FWHM / (2 * sqrt(ln(2))) ≈ 0.6 * FWHM.
• The voltage pulse satisfies the correct sign for power, need to have +V/I to the right & -V/I to the left.
• At right end V flips, I peaks; V = 0 (short circuit has Vmin, Imax).
• At left end: I flips, V peaks; I = 0 (open circuit has Imin, Vmax).

Source Behavior

• When the pulse returns to the source location, it may be offset for the reflected pulses to keep running into each other.
• Pulses reflected off a voltage source behave like a short circuit.
  • The internal resistance is zero (s.c.).
• Pulses reflected off a current source behave like an open circuit.
  • The internal resistance is infinite (o.c.).

Applying Pulse to Current

• If you apply the pulse to the current, rather than the voltage:
  • +V/-I to the right & -V/+I to the left.
• Pulses propagate at c ≈ 3*10^8 m/s.
• The voltage to current ratio was V/I ≈ 50 Ω.

Implementing a Soft Source

• A series current source or a shunt voltage source would produce a hard source.
• Series voltage source or shunt current source:
  • Add a line after the update equation to update using Vg & Ig.
• Pulses resembled the hard current source pulses (+V/+I right, -V/+I left) but with V0/2 peak amplitudes.
• Vo/2 amplitudes because Vo in the differential voltage applied relative to adjacent unit cell (rather than rel. to zero).

Terminating the Line

• A resistive loading can be implemented at the ends to represent open or short circuits, where previously the voltage and current values were unchanged
• Using careful time indexing is crucial because of spatial locality.
• In this example, the equations are:
  • V_0 = I_(1/2) R_s
  • I_(N+1/2) = V_N / R_L.
• With these equations:
  • You don't have V^(n+1/2) or I^(n+1).
    • The algorithm is using old / delayed data.
  • I_(N+1/2) & V_N and V_0 & I_(1/2) are offset by Δz/2.

Updating Equation

• It is crucial to use the correct differential equations to derive update equations for the first unit-cell current in the line I[1] and the last unit-cell voltage in the line V[N].
• End of transmission line with a resistor in series with the inductor:
  • −V_1 = I_(1/2)*R_s + L*dI_(1/2)/dt
  • Update the equation as follows: I_(1/2)^(n+1/2) = (L / Δt + R_s / 2)^(-1) * (I_(1/2)^(n-1/2) + (V_0^n / (L * Δz / Δt)))
• Equations for the right end (load) of the transmission line where the resistor is in parallel with the capacitor:
  • I_(N-1/2) = V_N / R_L + C*dV_N/dt
  • V_N^(n+1) = ((C * Δz) / Δt + 1 / (2 * R_L))^(-1) * (((C * Δz) / Δt - 1 / (2 * R_L)) *V_N^n + I_(N-1/2)^n

Summary of Updating

• Update I_(3/2) ... I_(N-1/2) from V_1 ... V_N.
• Then I_(1/2) from V_1 & R_s.
• Update V_1 to V_(N-1) from I_(3/2) ... I_(N-1/2).
• Then V_N from I_(N-1/2) & R_L.
• Small R flips V, not I (like a s.c.).
• Large R flips I, not V (like an o.c.).
• Reflected amplitude decreased as R approaches 50 Ω (but still observed some reflections!).
• Voltage and current should be co-localized at the terminals.

Model Corrections

• Must therefore add a factor of one-half to L and C in your update equations for I[1] and V[N].
• Effectively reduces the length of the transmission line by Δz/2.
• Can add a voltage source to the left end of the line (in series with the resistor) or a current source to the right end of the line (in shunt with the resistor) by modifying the V_1^(n+1/2) term in the I[1] update equation or the I_(N+1/2)^(n+1/2) term in the V[N] update equation.
• For voltage source: V_G - V_1^n = L * (dI_(1/2)/dt) + R_s * I_(1/2)^n.
• For current source: I_(N-1/2)^(n+1/2) - I_G = (C * Δz / 2) * (dV_N / dt) + V_N / R_L.

Other Sources and Inhomogeneities

• The most common way to model a transmission line is with a Thévenin source on the left and a passive load on the right.
• Implement a step or impulse function.
  • Saw distortion in signals (the ringing of HF components).
  • Different frequencies travel at slightly different speeds ("dispersion").
  • Impulse has a flat spectrum.
  • Step has a 1/f spectrum (less distortion with step vs. impulse).
• Cosine has discontinuity introduces HF components & distortion.
• Try smoothing functions which smoothly increase from f(0) = 0 to f(1) = 1 to eliminate ringing.
  • Example: f(t) = 3t^2 - 2t^3.
    • Continuous from 0 to 1 with continuous derivative also.

Inhomogeneities

• The final modification to allow for inhomogeneities in our transmission line is to allow our transmission-line parameters to change along its length.
• Replace L and C constants with vectors representing the parameter values in each unit cell and then using element-wise multiplication and division in our update equations (in Matlab, this is done using the .* and ./ commands).
• Smaller reflections at the interface, but no longer at load / Pulse slowed by V & I, were unchanged.
• If both L and C change by the same factor, there should not be a mismatch.
• To improve the example: Place the boundary at spatial index, k=100.
• Averaging helps improve these models.