Difference between revisions of "Basic Principles of The Finite Difference Time Domain Method"

Revision as of 15:02, 22 May 2017

An Overview of FDTD Modeling

In the Finite Difference Time Domain (FDTD) method, a discretized form of Maxwell’s equations is solved numerically and simultaneously in both the 3D space and time. During this process, the electric and magnetic fields are computed everywhere in the computational domain and as a function of time starting at t = 0. From knowledge of the primary fields in space and time, one can compute other secondary quantities including frequency domain characteristics like scattering parameters, input impedance, far-field radiation patterns, radar cross section, etc.

Since FDTD is a finite domain numerical technique, the computational domain of the problem must be truncated. At the boundaries of the computational domain, proper boundary conditions must be enforced. In a shielded structure, all objects are enclosed within a perfect electric (or magnetic) conductor box. In an open boundary problem like an antenna, some kind of absorbing boundary conditions such as a perfectly matched layer (PML) must be used to emulate the free space. The absorbing boundaries should act such that the incident fields and waves propagate through them without any back reflection. The FDTD simulation time depends directly on the size of the computational domain and on how close you can place the PML walls to the enclosed objects.

The FDTD computational domain must be discretized using an appropriate meshing scheme. EM.Tempo uses a non-uniform, adaptive, voxel-based Yee mesh with a mesh density that you can customize. A fixed-cell mesh generator is also available, where you can set constant cell dimensions along the three principal axes for the entire computational domain. The variable mesh density is specified in terms of the effective wavelength inside material media. As a result, the mesh resolution and average mesh cell size differ in regions that are filled with different types of material. EM.Tempo's adaptive mesher generates more cells in the areas that are occupied by dielectric materials, fewer cells in the free-space regions and no cells inside PEC regions. EM.Tempo's default adaptive mesh generator also refines the mesh around curved segments of lines, surface or solids to produce a far more accurate representation of your geometry.

The FDTD method provides a wideband simulation of your physical structure. In order to produce sufficient spectral information, an appropriate wideband temporal waveform is needed to excite the physical structure. The choice of the waveform, its bandwidth and time delay all affect the convergence behavior of the FDTD time marching loop. By default, EM.Tempo uses a modulated Gaussian waveform with optimal parameters. Another issue of concern is the numerical stability of the time marching scheme. You might expect to get better and more accurate results if you keep increasing the FDTD mesh resolution. However, in order to satisfy the Courant-Friedrichs-Levy (CFL) stability condition, the time step must be inversely proportional to the maximum grid cell size . A high resolution mesh requires a smaller time step. To let the fields in the computational domain fully evolve over time, a smaller time step will require a larger number of time steps to converge. EM.Tempo automatically chooses a time step that satisfies the CFL condition.

Differential Form of Maxwell's Equations & the Yee Cell

Maxwell’s equations for an isotropic, time-invariant and homogeneous medium are given by:

[math] \nabla \times \mathbf{E} = - \dfrac{\partial \mathbf{B}}{\partial t} - \mathbf{M} [/math]

[math] \nabla \times \mathbf{H} = \dfrac{\partial \mathbf{D}}{\partial t} + \mathbf{J} [/math]

where E and H are the electric and magnetic fields, respectively, D and B are the electric and magnetic flux densities, respectively, J and M are the volume electric and magnetic current densities, respectively, and the following constitutive relationships hold:

[math] \mathbf{D} = \epsilon \mathbf{E}, \quad \quad \mathbf{J} = \sigma \mathbf{E} [/math]

[math] \mathbf{B} = \epsilon \mathbf{H}, \quad \quad \mathbf{M} = \sigma_m \mathbf{H} [/math]

where ε is the permittivity, μ is the permeability, σ is the electric conductivity, and σ_m is the magnetic conductivity of the medium.

In a medium without electric or magnetic losses these equations can be rewritten as:

[math]\dfrac{\partial \mathbf{H}}{\partial t} = -\dfrac{1}{\mu} \nabla \times \mathbf{E}[/math]

[math]\dfrac{\partial \mathbf{E}}{\partial t} = \dfrac{1}{\epsilon} \nabla \times \mathbf{H}[/math]

The Yee cell geometry and definition of the electric and magnetic field components.

In the FDTD method, both time- and space-derivatives are approximated with central finite differences. This results in six differential equations, one for each field component, on a Yee cell define at a point in the 3D space as shown in the figure on the right:

[math] H_x^{n+\frac{1}{2}} (i,j,k) = H_x^{n-\frac{1}{2}} (i,j,k) + \frac{\Delta t}{\mu (i,j,k)} \left[ \frac{E_{y}^{n}(i,j,k+1) - E_{y}^{n}(i,j,k)}{\Delta z} - \frac{E_{z}^{n} (i,j+1,k)-E_{z}^{n} (i,j,k)}{\Delta y} \right] [/math]

[math] H_y^{n+\frac{1}{2}} (i,j,k) = H_y^{n-\frac{1}{2}} (i,j,k) + \frac{\Delta t}{\mu (i,j,k)} \left[ \frac{E_{z}^{n}(i+1,j,k) - E_{z}^{n}(i,j,k)}{\Delta x} - \frac{E_{x}^{n} (i,j,k+1)-E_{x}^{n} (i,j,k)}{\Delta z} \right] [/math]

[math] H_z^{n+\frac{1}{2}} (i,j,k) = H_z^{n-\frac{1}{2}} (i,j,k) + \frac{\Delta t}{\mu (i,j,k)} \left[ \frac{E_{x}^{n}(i,j+1,k) - E_{x}^{n}(i,j,k)}{\Delta y} - \frac{E_{y}^{n} (i+1,j,k)-E_{y}^{n} (i,j,k)}{\Delta x} \right] [/math]

[math] {E_{x}^{n+1} (i,j,k) = E_{x}^{n} (i,j,k)} + \frac{\Delta t}{\epsilon (i,j,k)} \left[ \frac{H_{z}^{n+\frac{1}{2} } (i,j,k) - H_{z}^{n+\frac{1}{2} } (i,j-1,k)}{\Delta y} - \frac{H_{y}^{n+\frac{1}{2} } (i,j,k) - H_{y}^{n+\frac{1}{2} } (i,j,k-1)}{\Delta z} \right] [/math]

[math] {E_{y}^{n+1} (i,j,k) = E_{y}^{n} (i,j,k)} + \frac{\Delta t}{\epsilon (i,j,k)} \left[ \frac{H_{x}^{n+\frac{1}{2} } (i,j,k) - H_{x}^{n+\frac{1}{2} } (i,j,k-1)}{\Delta z} - \frac{H_{z}^{n+\frac{1}{2} } (i,j,k) - H_{z}^{n+\frac{1}{2} } (i-1,j,k)}{\Delta x} \right] [/math]

[math] {E_{z}^{n+1} (i,j,k) = E_{z}^{n} (i,j,k)} + \frac{\Delta t}{\epsilon (i,j,k)} \left[ \frac{H_{y}^{n+\frac{1}{2} } (i,j,k) - H_{y}^{n+\frac{1}{2} } (i-1,j,k)}{\Delta x} - \frac{H_{x}^{n+\frac{1}{2} } (i,j,k) - H_{x}^{n+\frac{1}{2} } (i,j-1,k)}{\Delta y} \right] [/math]

where i, j, k are the grid position indices along the X, Y, Z axes, respectively, n is the current time step index and Δx, Δy, Δz are the grid position indices along the X, Y, Z axes, respectively. Similar expressions are obtained for the Y and Z components of the electric and magnetic fields. When your physical structure involves lossy materials with nonzero electric conductivity σ and/or nonzero electric conductivity σ_m, the above update equations become more complicated. In the case of anisotropic materials with tensorial constitutive parameters, the electric displacement vector D and magnetic induction vector B need to be involved in the update of Maxwell's equations at every time step. This results in a total of twelve update equations at every time step. In the case of dispersive materials with time-varying constitutive parameters, additional auxiliary differential equations are invoked and updated at every time step. Applying the proper boundary conditions for all the materials inside the computational domain and at the boundaries of the domain itself, EM.Cube calculates and "updates" all the necessary field components at every mesh node, at every time step.

Why Does FDTD Need Domain Termination?

The FDTD simulation time depends directly on the size of the computational domain. For free space radiation or scattering problems, the computational domain must be extended to infinity, which means an infinite number of cells in the computational domain. The solution to this problem is to truncate the domain by a set of artificial boundaries at a certain distance from the objects in the computational domain. The absorbing boundaries should be such that the field propagates through them without any back reflection. Different methods have been used to simulate an absorbing boundary condition in FDTD simulations. The most common ones are Mur, Liao, and the perfectly match layer (PML). The Mur boundary condition calculates the boundary field values from the three dimensional scalar wave equations, while the Liao boundary condition is based on extrapolation of the fields in space and time. In 1994, Berenger proposed a new boundary condition called the perfectly matched layer (PML), which provides a much better performance than the Liao and Mur boundary conditions. The PML medium properties surrounding the computational domain are chosen to effectively absorb all the outgoing waves propagating towards the boundaries. In PML regions, an artificial conductivity is introduced such that it starts with very small values at the free space-PML interfaces and gradually increases until it reaches its maximum value at the last layer of the PML region.

CPML vs. PML

The PML boundary condition is not effective in absorbing evanescent waves, and it suffers from late-time reflections when simulating fields with very long time signatures. This is partly due to the weakly causal nature of the PML. A strictly causal form of the PML, known as the complex frequency-shifted PML (CFS-PML), was later developed by simply shifting the frequency dependent pole off the real axis and into the negative-imaginary half of the complex plane. It has been shown that the CFS-PML is highly effective at absorbing evanescent waves and signals with a long time signature. Therefore, using the CFS-PML, the boundaries can be placed closer to the objects in the computational domain, and considerable time and memory savings can be achieved. The convolutional PML (CPML) is an efficient implementation of the CFS-PML based on a stretched coordinate formulation in conjunction with recursive convolution. It has been shown that CPML requires only two auxiliary variables per discrete field point and absorbs waves in isotropic, homogeneous, inhomogeneous, lossy, dispersive, anisotropic or non-linear media without any further generalization.

The boundary CPML cells placed outside the visible domain box.

The CPML Formulation

The CPML is formulated in the stretched coordinate space. The CPML layers are assumed to terminate the FDTD computational domain. The X components of Maxwell's frequency-domain curl equations can then be written in the following form:

[math] (j\omega\epsilon_x + \sigma_{ex})\tilde{E}_x = \frac{1}{s_{ey}} \frac{\partial \tilde{H}_z}{\partial y} - \frac{1}{s_{ez}} \frac{\partial \tilde{H}_y}{\partial z} [/math]

[math] (j\omega\mu_x + \sigma_{mx}) \tilde{H}_x = -\frac{1}{s_{my}} \frac{\partial \tilde{E}_z}{\partial y} + \frac{1}{s_{mz}} \frac{\partial \tilde{E}_y}{\partial z} [/math]

where s_ei and s_mi are the stretched coordinate metrics defined by:

[math] s_{ei} = \kappa_{ei} + \frac{\sigma_{ei}}{\alpha_{ei} + j\omega\varepsilon_0}, \quad i=x,y,z [/math]

[math] s_{mi} = \kappa_{mi} + \frac{\sigma_{mi}}{\alpha_{mi} + j\omega\mu_0}, \quad i=x,y,z [/math]

s_ei and s_mi are the anisotropic components of the synthesized electric and magnetic conductivities in the CPML region. κ_ei , κ_mi, α_ei and α_mi are all assumed to be positive real and κ_ei, κ_mi ≥ 1. Similar equations hold for the Y and Z components of the electric and magnetic fields in the CPML layers. The requirement for zero reflection at PML-PML interfaces imposes the following condition:

[math] s_{ei} = s_{mi} \quad \Rightarrow \quad \kappa_{ei} = \kappa_{mi} ,\quad \frac{\sigma_{ei}}{\varepsilon_0} = \frac{\sigma _{mi}}{\mu_0}, \quad \frac{\alpha_{ei}}{\varepsilon_0} = \frac{\alpha_{mi}}{\mu _0} [/math]

The tilde notation above denotes the Fourier transform of the field components in the frequency domain. Transforming the above equations back to the time domain, one encounters convolution on the right hand side due to the frequency dependence of the stretched coordinate metrics:

[math] \left(\varepsilon_x \frac{\partial E_x}{\partial t} + \sigma_{ex}E_x \right) = \breve{s}_{ey}(t) \frac{\partial H_z}{\partial y} - \breve{s}_{ez}(t)\frac{\partial H_y}{\partial z} [/math]

[math] \left(\mu_x \frac{\partial H_x}{\partial t} + \sigma_{mx}H_x \right) = -\breve{s}_{my}(t)\frac{\partial E_z}{\partial y} + \breve{s}_{mz}(t)\frac{\partial E_y}{\partial z} [/math]

where š_ei(t) and š_mi(t) denote functions of time that are indeed the inverse Laplace transform of s_ei^-1 and s_mi^-1,respectively, given by:

[math] \breve{s}_{ei}(t) = \frac{\delta(t)}{\kappa_{ei}} - \frac{\sigma_{ei}}{\kappa_{ei}^2} u(t)\exp \left[-\left(\frac{\sigma_{ei}}{\kappa_{ei}} + \alpha_{ei} \right)\frac{t}{\varepsilon_0} \right] [/math]

[math] \breve{s}_{mi}(t) = \frac{\delta(t)}{\kappa_{mi}} - \frac{\sigma_{mi}}{\kappa_{mi}^2} u(t)\exp \left[-\left(\frac{\sigma_{mi}}{\kappa_{mi}} + \alpha_{mi} \right)\frac{t}{\mu_0} \right] [/math]

where d(t) and u(t) denote the Dirac delta and unit step functions, respectively. The convolutions on the right hand side of the time domain equations can be accelerated by the use of the recursive convolution (RC) method.

The CPML parameters are chosen to be an increasing function of the distance from the boundaries of the computational domain. EM.Cube uses a polynomial profile of degree n_PML. Given the interrelationships among these parameters, one can write:

[math] \sigma_{ei}(r) = \frac{\varepsilon_0}{\mu_0}\sigma_{mi}(r) = \sigma_{max}\left( \dfrac{r}{\delta} \right)^{n_{PML}} [/math]

[math] \kappa_{ei}(r) = \kappa_{mi}(r) = 1 + \left( \kappa_{max} - 1 \right) \left( \dfrac{r}{\delta} \right)^{n_{PML} } [/math]

[math] \alpha_{ei}(r) = \frac{\varepsilon_0}{\mu_0} \alpha_{mi}(r) = \alpha _{min} + \left(\alpha_{max } - \alpha_{min } \right) \left( \frac{r}{\delta} \right)^{n_{PML}} [/math]

where r is the distance of field observation point inside the CPML layer from the edge of the computational domain. The parameters σ_max, κ_max, α_min and α_max as well as n_PML can be modified by the user.

EM.Tempo's Waveform Types

The accuracy of the FDTD simulation results depends on the right choice of temporal waveform. The excitation waveform is the temporal source function that sets the initial conditions at t = 0 and provides an excitation signal afterwards. EM.Tempo currently offers four different choices:

• Sinusoidal
• Gaussian Pulse
• Modulated Gaussian Pulse
• Custom

EM.Tempo's default waveform choice is a modulated Gaussian pulse. If you use a Gaussian pulse or a modulated Gaussian pulse waveform to drive your FDTD source, after a certain number of time steps, the total energy of the computational domain drops to very negligible levels. At that point, you can consider your solution to have converged. If you drive your FDTD source by a sinusoidal waveform, the total energy of the computational domain will oscillate indefinitely, and you have to force the time loop to terminate after a certain number of time steps assuming a steady state have been reached.

EM.Tempo's standard sinusoidal excitation waveform.

EM.Tempo's standard Gaussian pulse excitation waveform.

EM.Tempo's standard modulated Gaussian pulse excitation waveform.

Waveform, Bandwidth, Stability

The FDTD method provides a wideband simulation of your physical structure. Frequency domain techniques often require a tedious frequency sweep to calculate the port characteristics (S/Y/Z parameters). By contrast, EM.Tempo performs a discrete Fourier transform (DFT) of the time domain data to calculate these characteristics at the end of a single FDTD simulation run. In order to produce sufficient spectral information, an appropriate wideband temporal waveform is needed to excite the physical structure. The general form of EM.Tempo's default excitation waveform is a Modulated Gaussian Pulse given by:

[math] E_j^{inc}(r,t)=E_0(r) \exp \left(-\frac{(t-t_0)^2}{\tau ^2} \right) \cos \left(2\pi f_0 (t-t_0)-\Phi \right),\quad j=x,y,z [/math]

where f₀ is the center frequency, t₀ is the time delay, t is the Gaussian pulse width, and Φ is a constant phase. In the limits, the above waveform can be reduced either to a simple Gaussian pulse:

[math] E_j^{inc}(r,t) = E_0(r) \exp \left(-\frac{(t-t_0)^2}{\tau ^2} \right),\quad j=x,y,z [/math]

or to a continuous, single-tone, sinusoidal waveform with a frequency of f₀:

[math] E_j^{inc}(r,t) = E_{0} (r)\cos (2\pi f_0 (t-t_0)-\Phi),\quad j=x,y,z [/math]

The choice of the waveform, its bandwidth and time delay are important for the convergence behavior of the FDTD time marching loop. By default, EM.Tempo uses a modulated Gaussian waveform with optimal parameters: t = 0.966/Δf and t₀ = 4.5t, where Δf is the specified bandwidth of the simulation. The time delay t₀ is chosen so that the temporal waveform has an almost zero value at t = 0. Of the above waveforms, modulated Gaussian and sinusoidal waveforms are band pass with no DC content, while the Gaussian pulse is low pass with a frequency spectrum that is concentrated around f = 0. In a typical FDTD simulation, you set a center frequency for the structure of interest and then specify a bandwidth around this center frequency. These together determine the lowest and highest spectral contents of your FDTD waveform. Note that setting a bandwidth equal to 2f₀ sets the lowest frequency to DC (f_min = 0), which you may want to avoid in certain applications. On the other hand, using a Gaussian pulse waveform, you do want to set Δf = 2f₀. In contrast to the wideband, exponentially decaying, Gaussian pulse and modulated Gaussian waveforms, the sinusoidal waveform is extremely narrowband and single-frequency indeed. It does not decay over time and continues to oscillate indefinitely after reaching a steady state.

Another issue of concern in an FDTD simulation is the numerical stability of the time marching scheme. You can set the mesh grid cell size to any fraction of a wavelength. Normally, you would expect to get better and more accurate results if you increase the mesh resolution. However, the time step is inversely proportional to the maximum grid cell size in order to satisfy the Courant-Friedrichs-Levy (CFL) stability condition:

[math] \Delta t \le \frac{K_{CFL}} {c\sqrt{\left(\dfrac{1}{\Delta x_{min}} \right)^2 + \left(\dfrac{1}{\Delta y_{min}} \right)^2 + \left(\dfrac{1}{\Delta z_{min}} \right)^2 } } [/math]

where c is the speed of light, and K_CFL is a constant. EM.Tempo uses a default value of K_CFL = 0.9. For a uniform grid with equal cell dimensions along the X, Y and Z directions, i.e. Δx = Δy = Δz = Δ, and the CFL condition reduces to:

[math] \Delta t \le K_{CFL} \frac{\Delta}{\sqrt{3}c}[/math]

As can be seen from the above criterion, a high resolution mesh requires a smaller time step. Since you need to let the fields in the computational domain fully evolve over time, a smaller time step will require a larger number of time steps to achieve convergence. EM.Tempo automatically chooses a time step that satisfies the CFL condition.

The Relationship Between Excitation Waveform and Frequency-Domain Characteristics

The accuracy of the FDTD simulation results depends on the right choice of temporal waveform. At the end of an FDTD simulation, the time domain field data are transformed into the frequency domain at your specified frequency or bandwidth to produce the desired observables. The frequency at which the frequency-domain characteristics are computed is usually the same as the center frequency of your project, but it doesn't have to be as long as it falls within the project's bandwidth. Different waveforms have different spectral behaviors. It is critical that you understand the time-domain-to-frequency-domain transformations in order to correctly interpret some simulation data.

If you open up the property dialog of any source type in EM.Tempo, you will see an Excitation Waveform... button located in the "Source Properties" section of the dialog. Clicking this button opens up EM.Tempo's Excitation Waveform dialog. From this dialog, you can override EM.Tempo's default waveform and customize your own temporal waveform.

The sinusoidal waveform has a frequency f₀ equal to the project's center frequency and does not have any parameter that you can modify. The delay in this case is set to t₀ = 1/(4f₀) so that the waveform vanishes at t = 0:

where sinc(x) = sin(πx) / (πx), and T is the length of the temporal window from t = 0 to the exit time t = NΔt, and Δt is the time step of the time marching loop.

For Gaussian and modulated Gaussian waveforms you can set the pulse width and time delay. The pulse width should be chosen such that the Fourier transform of your selected waveform encompasses your project bandwidth. Rather than setting an arbitrary pulse width, EM.Cube lets you set a parameter called Spectral Truncation Level or delta (δ), which is the normalized truncation threshold for the spectral content of the excitation pulse. This is explained as follows. Recall that the Fourier transform of a Gaussian pulse is also a Gaussian pulse. The Fourier transform of a modulated Gaussian pulse consists of two Gaussian pulses located at the two opposite sides of f = 0 (DC) shifted by the modulation frequency f₀:

Since a Gaussian pulse waveform has considerable DC content, it must be used to model lowpass structures. In that case, the center frequency and bandwidth of the project must be set such that Δf = 2f₀ and hence, f_min = f₀ - Δf/2 = 0, and f_max = f₀ + Δf/2 = 2f₀. The width τ of the temporal Gaussian pulse is then determined such that at f_max the spectral Gaussian pulse drops to the δ-level from its maximum value of 1. With a bandwidth of Δf, the pulse width must satisfy the following equation:

[math] \exp[-(\pi f_{\delta} \tau)^2] = \exp[-(\pi \Delta f \tau)^2] = \delta \quad \Rightarrow \quad \tau = \frac{\sqrt{-\ln\delta}}{\pi f_{max}} [/math]

When you set δ = 0.1 (the default value), it means that the Fourier transform of your excitation waveform drops to 10% of its peak at the upper edge of your specified frequency range. For a Modulated Gaussian pulse waveform with a bandwidth of Δf, the pulse width must satisfy the following equation:

[math] \exp[ -[\pi(f_{\delta} - f_0)\tau] ^2 ] = \exp \left[ -\left(\pi \frac{\Delta f}{2} \tau\right)^2 \right] = \delta \quad \Rightarrow \quad \tau = \frac{2\sqrt{-\ln \delta}}{\pi \Delta f} [/math]

For the default value of δ = 0.1, you get the standard relation: τ = 0.966 / BW, which is typically used for modulated Gaussian waveforms. The source waveform requires a time delay t₀ so that its value drops to almost zero at t = 0. The time delay is expressed in terms of a multiple of the pulse width: t0 / tau, and its default value is 4.5.

Keep in mind that FDTD is a time domain algorithm. At the end of an FDTD simulation, a Discrete Fourier Transform (DFT) is performed on the time domain data to calculate frequency domain characteristics such as near fields, far field radiation patterns, RCS, S/Y/Z parameters, etc.:

[math]\tilde{F}(f) = \int_{-\infty}^{\infty} f(t) e^{-j2\pi f t} \, dt \quad \approx \quad \Delta t \sum_{n=0}^N f(n\Delta t) e^{-j2 \pi n f \Delta t} [/math]

Of EM.Tempo's observables, the near fields, far fields and all of their associated parameters like directivity, RCS, etc., are calculated at a certain frequency that is specified as part of the definition of the observable. On the other hand, port characteristics like S/Y/Z parameters, VSWR and periodic characteristics like reflection and transmission coefficients, are calculated over the entire specified bandwidth of your project. In other words, you get a wideband frequency response at the end of a single time domain FDTD simulation run. The number of frequency point data over this bandwidth is equal to the number of DFT samples that are generated during the time marching loop. This number is set to 200 by default. In other words, 200 frequency data are generated at the end of an FDTD simulation. You can change this number through the box labeled No. DFT Samples in the "Discrete Fourier Transform" section of the FDTD Simulation Engine Settings dialog.

Normalization of Frequency-Domain Field Data

From the above equations, it can be seen that the discrete Fourier transform multiplies the samples of time-domain field quantities by the time step Δt. This means that the resulting Fourier transforms of electric and magnetic field components now have units of V/m/Hz or A/m/Hz, respectively. Moreover, the Fourier transforms of the three waveform types have different spectral values at the observation frequency f₀. This makes it difficult to compare the FDTD simulation results with the results from EM.Cube's other computational modules. For example, in the Planar, MoM3D and Physical Optics Modules, a plane wave source typically has a complex-valued functional form of exp(-jk₀k.r), which has a unit magnitude. In EM.Tempo, the time domain plane wave source has a functional dependence of the following form:

[math] \mathbf{E^{inc}}(r,t) = (E_{\theta}^{inc} \hat{\theta} + E_{\phi}^{inc} \hat{\phi}) f \left[ (t-t_0) - \frac{\mathbf{\hat{k} \cdot r} - l_0}{c} \right] [/math]

where f(t) is the temporal waveform, t₀ is the time delay, l₀ is a spatial shift, and c is the speed of light in the free space. The Fourier transform of the above temporal function evaluated at f = f₀ has an amplitude different than 1. For this purpose, EM.Tempo normalizes the temporal waveform by the magnitude of its Fourier transform at the observation frequency f₀. As a result, a temporal plane source with any of the three waveform types will always create spectral incident source with |E^inc(r, f₀)| = 1. The temporal waveform normalization factors for the three waveform types are given below:

[math] \begin{align} &(N.F.)_{Sinusoidal} = \frac{2}{T} = \frac{2}{N \Delta t} \\ &(N.F.)_{Gaussian} = \frac{ e^{(\pi f_0 \tau)^2} }{ \sqrt{\pi} \tau } = \frac{1}{\sqrt{\pi} \tau \delta^{1/4}} \\ &(N.F.)_{Modulated} = \frac{2}{\sqrt{\pi} \tau} \end{align} [/math]

Time Domain Simulation of Periodic Structures

A periodic structure is one that repeats itself infinitely in one, two or three directions. EM.Tempo allows you to simulate doubly periodic structures with periodicities along the X and Y directions. Many interesting structures such as frequency selective surfaces (FSS), electromagnetic band-gap (EBG) structures and metamaterial structures can be modeled using periodic geometries. In the case of an infinitely extended periodic structure, it is sufficient to analyze only a unit cell. In the FDTD method, this is accomplished by applying periodic boundary conditions (PBC) at the side walls of the computational domain. The application of the PBC is straightforward for the case of a normally incident plane wave source since the fields do not experience any delay as they travel across the unit cell. Obliquely incident plane waves, on the other hand, cause a time delay in the transverse plane. This delay requires knowledge of the future values of the fields at any time step.

A number of techniques have been proposed to solve this problem. EM.Cube uses a recently developed novel technique that is known as Direct Spectral FDTD or Constant Transverse Wavenumber method. In this technique, the components of the transverse (horizontal) wavenumber are kept constant in the direction of periodicity. This technique shows a significant advantage over the other methods for simulation of the incident illuminations close to the grazing angles.

Diagram of a periodic structure illuminated by an obliquely incident plane wave in EM.Tempo.

The figure above shows a doubly periodic structure with periods S_x and S_y along the X and Y directions, respectively. The computational domain is terminated with PBC in both X and Y directions. Along the positive and negative Z directions, it is terminated with CPML layers. Bear in mind that the PBC is also applied to the CPML layers. The computational domain is excited by a TM_z or TE_z plane wave incident at z = z₀. The plane wave incidence angles are denoted by θ (elevation) and φ (azimuth) in the spherical coordinate system. The constant wavenumber components k_x and k_y in this case are defined as:

[math]k_x = k_0 \sin\theta\cos\phi[/math]

[math]k_y = k_0 \sin\theta \sin\phi[/math]

where [math]k_0 = \omega/c = 2\pi f/c = 2\pi/\lambda_0[/math] is the free space propagation constant, f is the operational frequency, ω is the angular frequency, λ₀ is the free space wavelength, c is the speed of light in the free space. The constant transverse wavenumber k_l is then given by:

[math] k_l = \sqrt{k_x^2 + k_y^2} = k_0\sin\theta [/math]

which depends only on θ and not on φ. On the excitation plane, the incident field adopts a modulated Gaussian waveform and a complex phase delay along the periodicity direction with the following form:

[math] H_x^{inc}(x,y,t) = -\frac{1}{\eta_0} \sin\phi \; \exp \left(-\frac{(t-t_0)^2}{\tau^2} \right) \exp(j2\pi f_0 t) \exp(-jk_x x) \exp(-jk_y y)[/math]

[math] H_y^{inc}(x,y,t) = \frac{1}{\eta_0} \cos\phi \; \exp \left(-\frac{(t-t_0)^2}{\tau^2} \right) \exp(j2\pi f_0 t) \exp(-jk_x x) \exp(-jk_y y)[/math]

for TM_z polarization and

[math] E_x^{inc}(x,y,t) = \sin\phi \; \exp \left(-\frac{(t-t_0)^2}{\tau^2} \right) \exp(j2\pi f_0 t) \exp(-jk_x x) \exp(-jk_y y)[/math]

[math] E_y^{inc}(x,y,t) = -\cos\phi \; \exp \left(-\frac{(t-t_0)^2}{\tau^2} \right) \exp(j2\pi f_0 t) \exp(-jk_x x) \exp(-jk_y y)[/math]

for TE_z polarization. Here, f₀ is the center frequency of the modulated Gaussian pulse waveform, t₀ is the time delay, and τ is the Gaussian pulse width. The choices of the Gaussian waveform parameters are very critical in order to avoid possible resonances. For a fixed value of k_l, the horizontal resonance occurs at:

[math] f_{res} = \frac{k_l c}{2 \pi} [/math]

For a fixed frequency [math]f_0[/math] and a fixed incidence angle [math]\theta_0[/math], the resonant frequency is reduced to:

[math] f_{res} = f_0 \mathbf{\sin} \theta_0 [/math]

The modulated Gaussian waveform must be chosen such that its effective bandwidth avoids the horizontal resonant frequency. Otherwise, the temporal response of the structure starts to oscillate, and the time marching loop will not converge. To avoid this problem, the modulation frequency and bandwidth of the waveform are chosen to satisfy the following condition:

[math] f_{mod} \ge f_{res} + \dfrac{1}{2}\Delta f = \dfrac{k_{l,fixed}\;c}{2\pi} + \dfrac{1}{2}\Delta f [/math]

Back to the Top of the Page

Back to Maxwell's Equations Page

Back to EM.Tempo Manual

Back to EM.Cube Main Page