COMPUTATION OF SCATTERING USING THE FINITE

DIFFERENCE TIME DOMAIN (FDTD) METHOD

Introduction

The finite difference time domain (FDTD) method is a means of determining the scattered field from objects by solving Maxwell’s equations in the time domain. The spatial and temporal differential operators are approximated by finite differences. Update equations can be derived (see reference [1]) that give the field at a point in space and time as a function of the field at the same and neighboring points at previous times. Therefore the solution is said to “march in time.”

As with most numerical solutions, the computational region is discretized into appropriate subdomains. For example, in 3-dimensional space they might be cubes. The incident wave is introduced into the computational grid and the scattered fields computed throughout the grid as a function of time. The fields at the boundaries of the computational grid are used to compute equivalent currents, which, in turn, are used in the radiation integrals to compute the far field. The resulting far fields are a function of time. They can be Fourier transformed to obtain the RCS as a function of frequency, .

Reference [1] describes the limitations, restrictions and tradeoffs involved in the computation. Important computational parameters and tradeoff issues include:

  • grid size,
  • time step, : determines the highest frequency of the computed RCS,
  • time window, ( time steps): frequency resolution (i.e., spacing between frequency samples after the Fourier transform) depends on the length of the time window. In most cases the fast Fourier transform (FFT) algorithm is used, which requires , M an integer. The frequency resolution is
  • and cannot be chosen independently and the Nyquist sampling criterion must be satisfied
  • the effects of terminating the computational grid must be addressed

The incident waveform is a Gaussian pulse. The pulse parameters are:

  • effective time duration,
  • effective bandwidth,
  • dynamic range,
  • time duration,
  • signal-to-noise ratio,

A Gaussian pulse is shown below:

The computer code fdtdpanel.m solves for the field scattered by a panel of inifinite extent, but finite thickness. The codes rcs2dbir.m rcs2inc.m, rcs2far.m and rcs2scat.m are used sequentially in order to compute the field scattered by an infinitely long square cylinder.

One-dimensional Solution: fdtdpanel.m

The program fdtdpanel.m uses the FDTD method to compute the scattering from an infinite panel of thickness . The waveform is a Gaussian pulse and the computational parameters are assigned in the code as shown below. The panel extends from diel_start to diel_stop. The conductivity, relative dielectric constant, and relative permeability can be chosen by the user. For each time step the fields are plotted. Lines representing the faces of the panel are also plotted as a reference. The oversampling ratio is K (if K=1 the Nyquist rate is used).

Setup for the program fdtdpanel.m:

%%% INPUT PARAMETERS %%%

K = 4; % Oversampling Ratio

N_samples= 2048*K; % Time Samples

v_rat = 1; % Velocity Ratio (=c/ugrid)

T_max = 1024*10^(-9); % Time Window

dT = T_max/N_samples; % Time Step

dL = 0.15/K; % Spatial Step

N_nodes = K*101; % Number of E-nodes

L = (N_nodes-1)*dL; % Spatial Domain Length

T0 = K*16*dT; % Truncated Gaussian Pulse Duration

T_shift = T0/2; % Gaussian Pulse Time Shift

B_eff = 10^9; % Effective Bandwidth

f_max = 1/(2*dT); % Nyquist Frequency

N_obs = round(N_nodes/(K*2));% Observation Point, Total Field

c = 3e8;

eta0 = 377;

ugrid = dL/dT;

%%% ASSIGN INPUT VARIABLES %%%

eps_rel = 4; % Relative permittivity is fixed

mu_rel = 1; % Relative permeability is fixed

sigma = 0; % Conductivity inside of panel

diel_start = 6; % Left edge of dielectric

diel_stop = 7; % Right edge of dielectric

%%% SET UP PLOT RANGE AND DRAW THE WALL &

Lmax=15;

Amax=0;

Amin=-40;

Xw1 = [diel_start, diel_start];

Xw2 = [diel_stop, diel_stop];

A snapshot of the computation is shown below for the data listed above. Note that the reflections from the front and rear faces are traveling in the reverse direction, while the transmitted wave is traveling in the forward direction.

Two-dimensional Scattering

The scattering from an infinitely long cylinder is computed using the four “rcs2xxx.m” codes in the following order:

rcs2dbir.m

rcs2inc.m (typical: 1-10 MHz for resolution; 0.5-1 oversampling ratio)

rcs2scat.m (a movie can be generated by using getframe with the array Mov)

rcs2far.m

A brief description of the setup (variable assignment) section of each code follows. More detailed explanations are given in reference [1], Chapter 4. Several helpful hints:

  1. Some of these programs were written under Matlab 4.x. Some warnings may be encountered, but the programs will run with Matlab 5.x. If the warnings are a nuisance they can be turned off by typing “warnings off”
  2. The programs may run a long time if a large number of FFT points are used.

The discretization of the cylinder is illustrated below:

Setup for rcs2dbir.m

N = 111;

Nin = 51;

Nout = Nin + 2;

h_dur = 50; % duration of impulse responses

vr = 1/sqrt(2); % grid "speed"

rowsh = 0.5*(N - Nin); % Row "shift" for the "inside" matrix

colsh = 0.5*(N - Nin); % Column "shift" for the "inside" matrix

irt = rowsh + 1; % Top row for the "inside" matrix

irb = irt + Nin - 1; % Bottom row for the "inside" matrix

icl = colsh + 1; % Leftmost column for the "inside" matrix

icr = icl + Nin - 1; % Rightmost column for the "inside" matrix

nodes_out = 4*(Nin+1); % Number of nodes on the boundary

nodes_in = 4*(Nin-1); % Number of nodes just inside the boundary

rcs2dbir.m writes the files dbir2par and dbir2dat.

Setup for rcs2inc.m:

load dbir2par

% constants

eps0 = 10^(-9)/(36*pi);

mu0 = 4*pi*10^(-7);

v0 = 1/sqrt(mu0*eps0);

% User Input

fmax = 10^9;

dt_Nyquist = 1/(2*fmax);

eps_rel = 1;

dt = dt_Nyquist/sqrt(eps_rel);

Df = input('Enter the desired frequency resolution in MHz:');

K_over = input('Enter the oversampling ratio:');

dt = dt/K_over;

Df = Df*10^6;

dl_air = sqrt(2)*v0*dt;

% Gaussian Pulse Parameters

B_eff = fmax;

To = 8/B_eff;

T_shift = To/2;

Gp_dur = ceil(To/dt);

% Time Samples

Tmax = 1/Df;

N_samples = 2^(ceil(log10(Tmax/dt)/log10(2)));

rcs2inc.m saves several variables which are used in rcs2scat.m

Setup for rcs2scat.m:

load dbir2dat

load dbir2par

load ez_inc

% constants

eps0 = 10^(-9)/(36*pi);

mu0 = 4*pi*10^(-7);

v0 = 1/sqrt(mu0*eps0);

% User Input

fmax = 10^9;

dt_Nyquist = 1/(2*fmax);

eps_rel = 1;

dt = dt_Nyquist/sqrt(eps_rel);

%K_over =input('Enter the oversampling ratio:');

dt = dt/K_over;

dl_air = sqrt(2)*v0*dt;

dl_diel = dl_air/sqrt(eps_rel);

%%%%% INPUT PARAMETERS %%%%%

N = Nin + 2;

nodes_out = 4*(N-1); % Number of nodes on the boundary

nodes_in = 4*(Nin-1); % Number of nodes just inside the boundary

% Time Samples

%N_samples = 512;

Tmax = N_samples*dt;

% Gaussian Pulse Parameters

B_eff = fmax;

To = 8/B_eff;

T_shift = To/2;

Example: Snapshot of the field for a square cylinder (computational parameters listed above). Top: mesh plot; bottom: contour plot. The plots are obtained from the program rcs2inc.m

The normalized scattered far-field pattern is shown above for a frequency of 100 MHz. The result is obtained by running the program rcs2far.m.

Reference:

[1] D. Jenn, Radar and Laser Cross Section Engineering, AIAA Education Series, 1995