A Simple, Method of Moments Solution for the Integral Equations for Multiple Dielectric Bodies of Arbitrary Shape in Time Domain

Sadasiva M. Rao

Naval Research Laboratory
Washington DC 20375, USA.
sadasiva.rao@nrl.navy.mil

Submitted On: June 21, 2023; Accepted On: January 20, 2024

ABSTRACT

In this work, we present a straightforward and simple method of moments (MOM) solution procedure, with minimum mathematical manipulations, to solve the coupled integral equations for multiple, homogeneous and inhomogeneous, dielectric bodies of arbitrary shape directly in the time domain. The standard surface and volume integral equation formulations are used for homogeneous and inhomogeneous bodies, respectively. The numerical solution procedure does not involve a time-marching process as is usually adopted for time domain problems and seems to be one of the primary reasons for the late-time instabilities as a result of error accumulation. The present solution method is stable for a very long time as evidenced by several representative numerical examples presented for validation.

Index Terms: dielectric bodies, integral equations, method of moments, time domain.

I. INTRODUCTION

Many applications developed in recent times, such as short pulse radar and 5G cellular systems, require the calculation of wide-band signature for efficient design and to avoid multi-band interference. The traditional frequency domain techniques for such calculations are expensive and time consuming. As a result, direct time domain techniques received considerable attention to develop user-friendly, general purpose, stable algorithms that can be easily adopted to a wide variety of problems. It may be noted that the goal of the present work is developing an efficient algorithm in time domain to solve a variety of problems and does not target a specific application problem.

The general purpose numerical methods in the time domain can be broadly classified into two categories $viz.$ a) Differential Equation (DE) solution techniques and b) Integral Equation (IE) solution techniques. The Finite Difference Time Domain (FDTD) is a popular DE solution method that is applicable to a variety of problems [1–5]. Similarly, Method of Moments (MOM) is the popular IE based solution for solving time domain electromagnetic problems.

In this work, we deal with IE solution methods only. One advantage of the IE method over the DE method is that the solution space is confined to the space occupied by the object. In contrast, for DE based methods the surrounding space also needs to be included in thesolution.

Until recently, a time-stepping process based on MOM, popularly known as the Marching-on-in-Time (MOT) method, has been the preferred technique for solving the numerical solution of Time Domain Integral Equation (TDIE) for electromagnetic field problems [6–11]. However, the MOT procedure is prone to late-time instabilities. The primary source of instability seems to be the accumulation of error at each time step eventually resulting into rapidly growing oscillations commonly known as late-time instabilities.

Recently, a new type of algorithm was developed using MOM that did not employ time-stepping procedure and remained stable for a very long time signature [12, 13]. In the present work, we apply this new method to multiple dielectric bodies of arbitrary shape using the surface equivalence principle [14], resulting in a surface integral equation (SIE) formulation or the volume equivalence principle [15] resulting in a volume integral equation (VIE) formulation. We note that SIE formulation is applicable only to homogeneous dielectric bodies whereas VIE formulation is applicable to both homogeneous and inhomogeneous volumes.

This work is organized as follows: In sections II and III we develop the integral equations along with the numerical solution scheme for the VIE and SIE formulations, respectively. In sections IV and V, we present several representative numerical results for validation purposes. Finally, in section VI, we summarize the work and present a few conclusions.

II. VIE FORMULATION

Although we are developing the solution procedure for multiple dielectric bodies, for the sake of brevity, let us consider only a single inhomogeneous body, illuminated by a Gaussian Plane Wave (GPW). We note that extending the present solution technique to multiple bodies is straightforward and does not require any special treatment.

A. Integral equation formulation

Let $V$ denote an arbitrarily-shaped, loss-less, inhomogeneous volume, surrounded by a homogeneous material $({ϵ}_0,{\mu }_0)$ as shown in Fig. 1. The material parameters inside the volume, $V$, are continuously changing from point to point and denoted by permittivity, $ϵ(\mathbf{r})$, and permeability, $\mu (\mathbf{r})$.

Figure 1: Inhomogeneous volume surrounded by a homogeneous material (free space).

An electric field, ${\mathbf{E}}^{inc}(\mathbf{r},t)$, defined in the absence of the scatterer, is incident on and induces polarization volume currents, $\mathbf{J}(\mathbf{r},t)$, in the volume, $V$. These volume currents then generate the scattered field which is a function of space and time. The total electric field is the sum of incident and scattered fields generated by the induced volume currents, and we have the following basic equations:

$\mathbf{J}(\mathbf{r},t)$	$=$	${\displaystyle \frac{\partial \mathbf{P}(\mathbf{r},t)}{\partial t}}={\displaystyle \frac{\partial [\mathbf{D}(\mathbf{r},t)-{ϵ}_0\mathbf{E}(\mathbf{r},t)]}{\partial t}}$	(1)
		$=$	$\left[{ϵ}_r(\mathbf{r})-1\right]{ϵ}_0{\displaystyle \frac{\partial \mathbf{E}(\mathbf{r},t)}{\partial t}},$	(2)
$\mathbf{E}(\mathbf{r},t)$	$=$	${\mathbf{E}}^s(\mathbf{r},t)+{\mathbf{E}}^{inc}(\mathbf{r},t),$	(3)
${\mathbf{E}}^s(\mathbf{r},t)$	$=$	$-{\displaystyle \frac{\partial \mathbf{A}(\mathbf{r},t)}{\partial t}}-\nabla \mathrm{\Phi }(\mathbf{r},t),$	(4)

where $\mathbf{P}(\mathbf{r},t)$ and $\mathbf{D}(\mathbf{r},t)$ represent the polarization vector and electric flux density vector, respectively. The magnetic vector potential, $\mathbf{A}(\mathbf{r},t)$, and the electric scalar potential, $\mathrm{\Phi }(\mathbf{r},t)$, are:

$\mathbf{A}(\mathbf{r},t)$	$=$	${\mu }_0{\displaystyle {\int }_V}{\displaystyle \frac{\mathbf{J}({\mathbf{r}}^{\prime },t-\frac{R}{c})}{4\pi R}}\mathit{d}{V}^{\prime },$	(5)
$\mathrm{\Phi }(\mathbf{r},t)$	$=$	${\displaystyle \frac{1}{{ϵ}_0}}{\displaystyle {\int }_V}{\displaystyle \frac{q_v({\mathbf{r}}^{\prime },t-\frac{R}{c})}{4\pi R}}\mathit{d}{V}^{\prime },$	(6)
$R$	$=$	$|\mathbf{r}-{\mathbf{r}}^{\prime }|.$	(7)

In Equations (5) and (6), $\mathbf{r}$ and ${\mathbf{r}}^{\prime }$ are the locations of the observation and source points on the scatterer, respectively. The volume charge density, $q_v$, is related to the polarization current, $\mathbf{J}$, by:

${\displaystyle \frac{\partial q_v}{\partial t}}$

$=$

$-\nabla \cdot \mathbf{J}$

(8)

Next, we can re-write Equations (2), using Equations (3) and (4), as:

${\displaystyle \frac{\mathbf{J}(\mathbf{r},t)}{{ϵ}_r(\mathbf{r})-1}}$	$=$	${ϵ}_0{\displaystyle \frac{\partial {\mathbf{E}}^{inc}(\mathbf{r},t)}{\partial t}}$	(9)
		$-$	${ϵ}_0\left[{\displaystyle \frac{{\partial }^2\mathbf{A}(\mathbf{r},t)}{\partial t^2}}-\nabla {\displaystyle \frac{\partial \mathrm{\Phi }(\mathbf{r},t)}{\partial t}}\right].$

Now, we can write the time derivative of the scalar potential, using Equation (8), as:

${\displaystyle \frac{\partial \mathrm{\Phi }(\mathbf{r},t)}{\partial t}}$

$=$

$-{\displaystyle \frac{1}{{ϵ}_0}}{\displaystyle {\int }_V}{\displaystyle \frac{\nabla \cdot \mathbf{J}({\mathbf{r}}^{\prime },t-\frac{R}{c})}{4\pi R}}\mathit{d}{V}^{\prime }.$

(10)

Next, we define the following relationships:

$\kappa (\mathbf{r})$	$=$	${\displaystyle \frac{{ϵ}_r-1}{{ϵ}_r}},$	(11)
$\mathbf{D}(\mathbf{r},t)$	$=$	${ϵ}_r{ϵ}_0\mathbf{E}(\mathbf{r},t),$	(12)
$\mathbf{J}(\mathbf{r})$	$=$	$\kappa (\mathbf{r}){\displaystyle \frac{\partial \mathbf{D}(\mathbf{r},t)}{\partial t}},$	(13)

where $\kappa (\mathbf{r})$ is the contrast ratio.

Note that the normal component of $\mathbf{D}$ is continuous at media interfaces and, hence, provides a convenient way to solve for the unknown quantity, $\mathbf{J}(\mathbf{r},t)$.

Now using Equations (11) - (13), and carrying out a few simple mathematical steps, we can write Equation (1) as:

${\displaystyle \frac{\mathbf{D}(\mathbf{r},t)}{{ϵ}_r}}+{\mu }_0{\displaystyle \frac{{\partial }^2}{\partial t^2}}{\displaystyle {\int }_V}{\displaystyle \frac{\kappa ({\mathbf{r}}^{\prime })\mathbf{D}({\mathbf{r}}^{\prime },t-\frac{R}{c})}{4\pi R}}\mathit{d}{V}^{\prime }$
$-{\displaystyle \frac{1}{{ϵ}_0}}\nabla \left[{\displaystyle {\int }_V}{\displaystyle \frac{{\nabla }^{\prime }\cdot \kappa ({\mathbf{r}}^{\prime })\mathbf{D}({\mathbf{r}}^{\prime },t-\frac{R}{c})}{4\pi R}}\mathit{d}{V}^{\prime }\right]={\mathbf{E}}^{inc}(\mathbf{r},t).$	(14)

Equation (14) is the required integro-differential equation that needs to be solved by numerical methods to obtain the unknown quantity, $\mathbf{D}(\mathbf{r},t)$.

B. MOM solution

As a first step, we define an upper limit on the time variable $t=T$, where $T$ represents the time when the incident pulse becomes negligible. Then, we divide the time axis $0\to T$ into $N_t$ uniform time intervals given by $\mathrm{\Delta }t$ and denote $t_n=n\mathrm{\Delta }t$ for $n=1,2,\mathrm{\dots }.N_t$. We note that, initially, the MOM scheme is applied to a finite interval $0\to T$. We also note that extending the time interval to later times is trivial and simply repeats the same steps from $0\to T$.

Next, we define the triangle functions to approximate the time variable in the interval $0\to T$ as:

$$g_n(t)\equiv \{\begin{array}{cc}\hfill 1-\frac{|t-t_n|}{\mathrm{\Delta }t}\hfill & t\in (t_{n-1},t_{n+1})\hfill \\ \hfill 0\hfill & \mathrm{otherwise}.\hfill \end{array},$$

(15)

for $n=1,2,\mathrm{\cdots },N_t$.

Now, we define the Schaubert-Wilton-Glisson (SWG) basis functions to represent the spatial variation of the electric flux density vector $\mathbf{D}$ [18].

Assuming a suitable tetrahedral model for the scattering structure, the basis functions are described as follows:

Figure 2: Description of the basis function.

Figure 2 shows two tetrahedrons $T_n^p$ and $T_n^q$, associated with the $n^{th}$ triangular surface denoted by $S_n$. Points in $T_n^p$ may be designated either by the position vector $\mathbf{r}$, or by ${\mathit{\rho }}_n^p$ defined with respect to the free vertex of $T_n^p$ . Similar remarks apply to the position vector $\mathbf{r}$ in $T_n^q$ except that it is directed toward the free vertex of $T_n^q$. It is assumed that the reference direction for $\mathbf{D}$ associated with the $n^{th}$ triangle is from $T_n^p$ to $T_n^q$.

Referring to Fig. 2, the vector basis functions associated with the $n^{th}$ triangular surface is:

$${\mathbf{f}}_n(\mathbf{r})=\{\begin{array}{cc}\frac{a_n}{3V_n^p}(\mathbf{r}-{\mathbf{r}}_n^p)\hfill & \mathbf{r}\in T_n^p\hfill \\ \frac{a_n}{3V_n^q}(\mathbf{r}-{\mathbf{r}}_n^q)\hfill & \mathbf{r}\in T_n^q\hfill \\ 0.0\hfill & \mathrm{otherwise},\hfill \end{array}$$

(16)

where $a_n$ is the area of the triangular surface, $V_n^p$ and $V_n^q$ are the volumes of the $p$ and $q$ tetrahedrons attached to the triangular surface $S_n$, and ${\mathbf{r}}_n^p$ and ${\mathbf{r}}_n^q$ represent the position vectors to the free vertex of tetrahedrons $p$ and $q$, respectively.

Now, we approximate the induced electric flux density, $\mathbf{D}(\mathbf{r},t)$, as:

$$\mathbf{D}(\mathbf{r},t)\approx \sum _{m=1}^{N_S}\sum _{n=1}^{N_t}{I}_{m,n}{\mathbf{f}}_m(\mathbf{r})g_n(t),$$

(17)

where $N_S$ and $N_t$ represent the number of basis (expansion) functions in space and time, respectively.

The next step in applying the method of moments is to select the testing procedure. As testing functions, we choose the same functions described in [18]. Defining:

$<{\mathbf{f}}_m(\mathbf{r})g_n(t),\mathbf{F}(\mathbf{r},t)>=$
${\displaystyle {\int }_V}{\displaystyle {\int }_T}{\mathbf{f}}_m(\mathbf{r})g_n(t)\bullet \mathbf{F}(\mathbf{r},t)\mathit{d}V\mathit{d}t,$	(18)

we write Equation (14) as:

$<{\mathbf{f}}_m(\mathbf{r})g_n(t),{\displaystyle \frac{\mathbf{D}(\mathbf{r},t)}{{ϵ}_r}}>$
$+<{\mathbf{f}}_m(\mathbf{r})g_n(t),{\mu }_0{\displaystyle \frac{{\partial }^2}{\partial t^2}}{\displaystyle {\int }_V}{\displaystyle \frac{\kappa ({\mathbf{r}}^{\prime })\mathbf{D}({\mathbf{r}}^{\prime },t-\frac{R}{c})}{4\pi R}}d{V}^{\prime }>$
$-<{\mathbf{f}}_m(\mathbf{r})g_n(t),{\displaystyle \frac{1}{{ϵ}_0}}\nabla \left[{\displaystyle {\int }_V}{\displaystyle \frac{{\nabla }^{\prime }\cdot \kappa ({\mathbf{r}}^{\prime })\mathbf{D}({\mathbf{r}}^{\prime },t-\frac{R}{c})}{4\pi R}}d{V}^{\prime }\right]>$
$=<{\mathbf{f}}_m(\mathbf{r})g_n(t),{\mathbf{E}}^{inc}>,$	(19)

for $m=1,2,\mathrm{\cdots }{N}_S$ and $n=1,2,\mathrm{\cdots }{N}_t$.

The first term, $<{\mathbf{f}}_m(\mathbf{r})g_n(t),\frac{\mathbf{D}(\mathbf{r},t)}{{ϵ}_r}>$, in Equation (19) may be written as:

$⟨{\mathbf{f}}_m(\mathbf{r})g_n(t),{\displaystyle \frac{\mathbf{D}(\mathbf{r},t)}{{ϵ}_r}}⟩$
$=\mathrm{\Delta }t[{\displaystyle \frac{a_m}{3V_m^i}}{\displaystyle {\int }_{T_m^i}}(\mathbf{r}-{\mathbf{r}}_m^i)\cdot {\displaystyle \frac{\mathbf{D}(\mathbf{r},t_n)}{{ϵ}_r}}dV$
$+{\displaystyle \frac{a_m}{3V_m^j}}{\displaystyle {\int }_{T_m^j}}(\mathbf{r}-{\mathbf{r}}_m^j)\cdot {\displaystyle \frac{\mathbf{D}(\mathbf{r},t_n)}{{ϵ}_r}}dV],$	(20)

where $T_m^i$ and $T_m^j$ represent the tetrahedrons attached to the triangular surface, $S_m$. The integrals in Equation (20) may be analytically evaluated after substituting the expansion functions, (17), for $\mathbf{D}(\mathbf{r},t)$.

Next, let us consider the second term in Equation (19):

$⟨{\mathbf{f}}_m(\mathbf{r})g_n(t),{\mu }_0{\displaystyle \frac{{\partial }^2}{\partial t^2}}{\displaystyle {\int }_V}{\displaystyle \frac{\kappa ({\mathbf{r}}^{\prime })\mathbf{D}({\mathbf{r}}^{\prime },t-\frac{R}{c})}{4\pi R}}\mathit{d}{V}^{\prime }⟩,$

which is approximated as:

$\left[{\displaystyle \frac{g(t_n)-2g(t_{n-1})+g(t_{n-2})}{\mathrm{\Delta }t}}\right]$
$\times [\kappa ({\mathbf{r}}_m^i)\mathbf{D}({\mathbf{r}}_m^{ci},t-{\displaystyle \frac{R_m^{ci}}{c}})\cdot {\displaystyle {\int }_{T_m^i}}{\mathbf{f}}_{m}d{V}^{\prime }$
$+\kappa ({\mathbf{r}}_m^j)\mathbf{D}({\mathbf{r}}_m^{cj},t-{\displaystyle \frac{R_m^{cj}}{c}})\cdot {\displaystyle {\int }_{T_m^j}}{\mathbf{f}}_{m}d{V}^{\prime }].$	(21)

Note that in Equation (21), $\mathbf{D}(\mathbf{r},t)$ is evaluated at ${\mathbf{r}}_m^{ci}$ and ${\mathbf{r}}_m^{cj}$, which represent the position vectors to centroids of the tetrahedrons $T_m^i$ and $T_m^j$, respectively. The integrals in Equation (21) are trivial and may be carried out analytically.

Using similar mathematical steps, the right hand side of the Equation (19) is written as:

$\mathrm{\Delta }t[{\mathbf{E}}^{inc}({\mathbf{r}}_m^{ci},t-{\displaystyle \frac{R_m^{ci}}{c}})\cdot ({\mathbf{r}}_m^{ci}-{\mathbf{r}}_m^i)$
$+{\mathbf{E}}^{inc}({\mathbf{r}}_m^{cj},t-{\displaystyle \frac{R_m^{cj}}{c}})\cdot ({\mathbf{r}}_m^{cj}-{\mathbf{r}}_m^j)]$	(22)

Lastly, we consider the third term in Equation (19):

$$⟨{\mathbf{f}}_m(\mathbf{r})g_n(t),\nabla \left[\frac{1}{{ϵ}_0}{\int }_V\frac{{\nabla }^{\prime }\cdot \kappa ({\mathbf{r}}^{\prime })\mathbf{D}({\mathbf{r}}^{\prime },t-\frac{R}{c})}{4\pi R}\mathit{d}{V}^{\prime }\right]⟩.$$

(23)

Here, we first perform testing on the time variable and the result is:

$\mathrm{\Delta }t⟨{\mathbf{f}}_m(\mathbf{r}),\nabla \left[{\displaystyle \frac{1}{{ϵ}_0}}{\displaystyle {\int }_V}{\displaystyle \frac{{\nabla }^{\prime }\cdot \kappa ({\mathbf{r}}^{\prime })\mathbf{D}({\mathbf{r}}^{\prime },t-\frac{R}{c})}{4\pi R}}\mathit{d}{V}^{\prime }\right]⟩.$

(24)

Denoting the term in the square bracket as $\mathrm{{\rm Y}}$, we can write the previous expression, after using the identity $\nabla \cdot (\mathbf{A}\mathrm{{\rm Y}})=\nabla \mathrm{{\rm Y}}\cdot \mathbf{A}+\mathrm{{\rm Y}}\nabla \cdot \mathbf{A}$, as:

$⟨{\mathbf{f}}_m,\nabla \mathrm{{\rm Y}}⟩$

$=$

${\displaystyle {\int }_V}\nabla \cdot ({\mathbf{f}}_m\mathrm{{\rm Y}})\mathit{d}V-{\displaystyle {\int }_V}\mathrm{{\rm Y}}\nabla \cdot {\mathbf{f}}_m\mathit{d}V.$

(25)

We note that the first integral in Equation (25) vanishes because of the properties of the basis functions and we are left with the following:

$-\mathrm{\Delta }t{\displaystyle {\int }_V}\nabla \cdot {\mathbf{f}}_m\left[{\displaystyle {\int }_V}\nabla \cdot \kappa ({\mathbf{r}}^{\prime })\mathbf{D}({\mathbf{r}}^{\prime })G(\mathbf{r},{\mathbf{r}}^{\prime })\mathit{d}{V}^{\prime }\right]\mathit{d}V.$

(26)

Next, we consider the expansion procedure. Substitution of the expansion function, Equation (17), into Equations (20), (21), and (26) yields a $P\times P$ system of linear equations, where $P=N_S\times N_t$. These equations may be written in a matrix form as:

$$\mathbf{ZI}=\mathbf{V},$$

(27)

where $\mathbf{Z}=[Z_{mn}]$ is an $P\times P$ matrix and $\mathbf{I}=[I_n]$ and $\mathbf{V}=[V_m]$ are column vectors of length, $P$. Obviously, it is possible to obtain the unknown vector $\mathbf{I}=[I_n]$ by inverting the $\mathbf{Z}$-matrix and multiplying by $\mathbf{V}$. However, there is a better and efficient way as described in the following:

Here, we note that the $\mathbf{Z}$-matrix is not a full matrix, unlike in the frequency domain MOM procedure. In fact, it is a lower triangular, block-wise, Toeplitz matrix given by:

$$\mathbf{Z}=\left[\begin{array}{cccc}\hfill {\mathbf{Z}}_{1,1}\hfill & \hfill \mathbf{O}\hfill & \hfill \mathrm{\cdots }\hfill & \hfill \mathbf{O}\hfill \\ \hfill {\mathbf{Z}}_{2,1}\hfill & \hfill {\mathbf{Z}}_{2,2}\hfill & \hfill \mathrm{\cdots }\hfill & \hfill \mathbf{O}\hfill \\ \hfill \mathrm{⋮}\hfill & \hfill \mathrm{⋮}\hfill & \hfill \mathrm{⋮}\hfill & \hfill \mathrm{⋮}\hfill \\ \hfill {\mathbf{Z}}_{N_S,1}\hfill & \hfill {\mathbf{Z}}_{N_S,2}\hfill & \hfill \mathrm{\cdots }\hfill & \hfill {\mathbf{Z}}_{N_S,N_S}\hfill \end{array}\right],$$

(28)

where each ${\mathbf{Z}}_{p,q},p=1,2,\mathrm{\cdots }{N}_S$ and $q=1,2,\mathrm{\cdots },N_S$, is a matrix of dimension, $N_S$, representing the mutual interaction between the spatial basis functions for a given pair of testing time function, $g_n$, and source time function, $g_k$. Further, because of the Toeplitz property, we have ${\mathbf{Z}}_{p,q}={\mathbf{Z}}_{|p-q|+1,1}$. Hence, we only need to compute the first column of Equation (28) and distribute the elements accordingly. In other words, we only have to compute the matrix elements for the first source time function and the testing time functions $1,2,\mathrm{\cdots },N_T$. The solution of such a matrix equation is very efficient, involves inverting only once, a matrix of size $N_S\times N_S$, and solving the matrix equation.

Lastly, the elements of the right hand side of Equation (14) are give by:

$V_m={ϵ}_0[{\mathbf{E}}^{inc}({\mathbf{r}}_m^{ci},t_k)\cdot ({\mathbf{r}}_m^{ci}-{\mathbf{r}}_m^i)$
$+{\mathbf{E}}^{inc}({\mathbf{r}}_m^{cj},t_k)\cdot ({\mathbf{r}}_m^{cj}-{\mathbf{r}}_m^j)].$	(29)

For a plane wave incidence, we set:

$${\mathbf{E}}^{inc}(\mathbf{r},t)=E_{\theta }(\mathbf{r},t-\frac{\mathbf{r}\cdot \widehat{k}}{c})\widehat{\theta }+E_{\varphi }(\mathbf{r},t-\frac{\mathbf{r}\cdot \widehat{k}}{c})\widehat{\varphi },$$

(30)

where the propagation vector, $\widehat{k}$, is given by:

$$\widehat{k}=\sin {\theta }_0\cos {\varphi }_0\widehat{x}+\sin {\theta }_0\sin {\varphi }_0\widehat{y}+\cos {\theta }_0\widehat{z},$$

(31)

and $({\theta }_0,{\varphi }_0)$ defines the angles of arrival of the plane wave in the usual spherical coordinate system.

On the right hand side of the matrix equation, $\mathbf{V}$ is obtained by using Equation (29) and consists of $N_t$ blocks of vectors of dimension $N_S$. At this stage, we note that multiple incident pulses with varying frequency content can be easily accommodated by adding more column blocks to the $\mathbf{V}$-matrix. Also, we note that obtaining currents for $T$ to $2T$ and later instants is similar to solving the equation for $0$ to $T$ as presented in [12] and [13].

Lastly, note that the numerical procedure presented so far allows us to obtain the current distribution on the scattering structure as a function of time. Once an accurate current distribution is obtained, it is a simple process to obtain near-fields, far-fields, and any other required parameters. The mathematical details to obtain such parameters are well-known and available in [16] and hence not repeated here.

III. SIE FORMULATION

Once again, we consider a single, homogeneous dielectric body for developing the integral equations for the sake of brevity.

A. Integral equation formulation

Let $S_d$ denote the surface of the dielectric body, surrounded by free space and illuminated by an incident plane wave pulse. The regions exterior and interior to the dielectric body, denoted by “e” and “i”, are characterized by medium parameters (${\mu }_0$, ${ϵ}_0$) and (${\mu }_d$, ${ϵ}_d$), respectively. Using the standard equivalence principle [14], and defining the equivalent currents, ${\mathbf{J}}_d$ and $M_d$, we derive the following equations:

${\displaystyle \frac{\partial }{\partial t}}{\left[{\mathbf{E}}_e^s({\mathbf{J}}_d)+{\mathbf{E}}_e^s({\mathbf{M}}_d)+{\mathbf{E}}^{inc}\right]}_{tan}=0,$	(32)
${\displaystyle \frac{\partial }{\partial t}}{\left[{\mathbf{E}}_i^s({\mathbf{J}}_d)+{\mathbf{E}}_i^s({\mathbf{M}}_d)\right]}_{tan}=0.$	(33)

Next, the time derivative of the scattered electric fields radiated by the equivalent electric and magnetic currents are written, in terms of potential functions, as:

$$\frac{\partial }{\partial t}\left[{\mathbf{E}}_{\nu }^s({\mathbf{J}}_d,{\mathbf{M}}_d)\right]=-\frac{{\partial }^2{\mathbf{A}}_{\nu }}{\partial t^2}-\frac{\partial \nabla {\mathrm{\Phi }}_{\nu }}{\partial t}-\frac{\partial }{\partial t}\left[\frac{1}{{ϵ}_{\nu }}\nabla \times {\mathbf{F}}_{\nu }\right],$$

(34)

where ${\mathbf{A}}_{\nu }$ and ${\mathbf{F}}_{\nu }$ are the magnetic and electric vector potentials, respectively, and ${\mathrm{\Phi }}_{\nu }$ is the electric scalar potential, given by:

${\mathbf{A}}_{\nu }(\mathbf{r},t)$	$=$	${\mu }_{\nu }{\displaystyle {\int }_S}{\displaystyle \frac{{\mathbf{J}}_d({\mathbf{r}}^{\prime },t-\frac{R}{c_{\nu }})}{4\pi R}}\mathit{d}{s}^{\prime }\text{,}$	(35)
${\mathbf{F}}_{\nu }(\mathbf{r},t)$	$=$	${ϵ}_{\nu }{\displaystyle {\int }_S}{\displaystyle \frac{{\mathbf{M}}_d({\mathbf{r}}^{\prime },t-\frac{R}{c_{\nu }})}{4\pi R}}\mathit{d}{s}^{\prime }\text{,}$	(36)
${\mathrm{\Phi }}_{\nu }(\mathbf{r},t)$	$=$	${\displaystyle \frac{1}{{ϵ}_{\nu }}}{\displaystyle {\int }_S}{\displaystyle \frac{q_e({\mathbf{r}}^{\prime },t-\frac{R}{c_{\nu }})}{4\pi R}}\mathit{d}{s}^{\prime }\text{,}$	(37)

for $\nu =e$ or $\nu =i$. In Equations (35) - (37), $R=|\mathbf{r}-{\mathbf{r}}^{\prime }|$ is the distance from the field point, r, to the source point, ${\mathbf{r}}^{\prime }$. The electric surface charge densities, $q_e$ is related to ${\mathbf{J}}_d$, by the continuity equation:

$$\nabla \cdot {\mathbf{J}}_d=-\frac{\partial q_e}{\partial t}\Rightarrow q_e=-{\int }_{\tau =0}^t\nabla \cdot {\mathbf{J}}_d\mathit{d}\tau .$$

(38)

Using Equation (38), Equation (37) is re-written as:

${\mathrm{\Psi }}_{\nu }={\displaystyle \frac{\partial {\mathrm{\Phi }}_{\nu }(\mathbf{r},t)}{\partial t}}=-{\displaystyle \frac{1}{{ϵ}_{\nu }}}$

${\displaystyle {\int }_S}{\displaystyle {\int }_{\tau =0}^t}{\displaystyle \frac{\nabla \cdot {\mathbf{J}}_d({\mathbf{r}}^{\prime },\tau -\frac{R}{c_{\nu }})}{4\pi R}}\mathit{d}{s}^{\prime }\mathit{d}\tau .$

(39)

Finally, we have:

${\left[{\displaystyle \frac{{\partial }^2{\mathbf{A}}_e}{\partial t^2}}+\nabla {\mathrm{\Psi }}_e+{\displaystyle \frac{\partial }{\partial t}}\left\{{\displaystyle \frac{1}{{ϵ}_e}}\nabla \times {\mathbf{F}}_e\right\}\right]}_{tan}$	$=$	${\left[{\displaystyle \frac{\partial {\mathbf{E}}^{inc}}{\partial t}}\right]}_{tan},$	(40)
${\left[{\displaystyle \frac{{\partial }^2{\mathbf{A}}_i}{\partial t^2}}+\nabla {\mathrm{\Psi }}_i+{\displaystyle \frac{\partial }{\partial t}}\left\{{\displaystyle \frac{1}{{ϵ}_i}}\nabla \times {\mathbf{F}}_i\right\}\right]}_{tan}$	$=$	$0.$	(41)

The integral equations (40) and (41) are solved as described in the next subsection.

B. MOM solution

Assuming a suitable triangulation for the scattering structure, ${\mathbf{J}}_d$ and ${\mathbf{M}}_d$ are approximated as:

${\mathbf{J}}_d(\mathbf{r},t)$	$=$	${\displaystyle \sum _{m=1}^{N_d}}{\displaystyle \sum _{n=1}^{N_t}}{\alpha }_{m,n}{\mathbf{f}}_m(\mathbf{r})g_n(t),$	(42)
${\mathbf{M}}_d(\mathbf{r},t)$	$=$	${\displaystyle \sum _{m=1}^{N_d}}{\displaystyle \sum _{n=1}^{N_t}}{\beta }_{m,n}{\mathbf{a}}_n\times {\mathbf{f}}_m(\mathbf{r})g_n(t),$	(43)

where $N_d$ and $N_t$ represent the number of basis functions in space and time, respectively, ${\mathbf{a}}_n$ is the normal vector,

$$g_n(t)\equiv \{\begin{array}{cc}\hfill 1-\frac{|t-t_n|}{\mathrm{\Delta }t}\hfill & t\in (t_{n-1},t_{n+1})\hfill \\ \hfill 0\hfill & \mathrm{otherwise},,\hfill \end{array}$$

(44)

for $n=1,2,\mathrm{\cdots },N_t$, and, ${\mathbf{f}}_m(\mathbf{r})$ are the standard Rao-Wilton-Glisson (RWG) functions [17].

Using the symmetric product defined in Equation (18), we can write Equations (40) and (41) as:

$<{\mathbf{f}}_m(\mathbf{r})g_n(t),[{\displaystyle \frac{{\partial }^2{\mathbf{A}}_e}{\partial t^2}}+\nabla {\mathrm{\Psi }}_e+{\displaystyle \frac{\partial }{\partial t}}\{{\displaystyle \frac{1}{{ϵ}_e}}\nabla \times {\mathbf{F}}_e\}]>$
$=<{\mathbf{f}}_m(\mathbf{r})g_n(t),\left[{\displaystyle \frac{\partial {\mathbf{E}}^{inc}}{\partial t}}\right]>,$	(45)
$<{\mathbf{f}}_m(\mathbf{r})g_n(t),[{\displaystyle \frac{{\partial }^2{\mathbf{A}}_i}{\partial t^2}}+\nabla {\mathrm{\Psi }}_i+{\displaystyle \frac{\partial }{\partial t}}\{{\displaystyle \frac{1}{{ϵ}_i}}\nabla \times {\mathbf{F}}_i\}]>,$
$=0,$	(46)

for $m=1,2,\mathrm{\cdots }{N}_d$ and $n=1,2,\mathrm{\cdots }{N}_t$.

The first term, $<{\mathbf{f}}_m(\mathbf{r})g_n(t),\frac{{\partial }^2{\mathbf{A}}_e}{\partial t^2}>$, in Equation (45) may be written as:

$⟨{\mathbf{f}}_m(\mathbf{r})g_n(t),{\displaystyle \frac{{\partial }^2{\mathbf{A}}_e}{\partial t^2}}⟩={\displaystyle \frac{{\partial }^2}{\partial t^2}}\left[\mu {\displaystyle {\int }_S}{\displaystyle \frac{{\mathbf{J}}_d({\mathbf{r}}^{\prime },t-\frac{R}{c_e})}{4\pi R}}\right]$
$=\mu \left[{\displaystyle \frac{{\partial }^2{g}_n}{\partial t^2}}\right]{\displaystyle {\int }_S}{\displaystyle \frac{{\mathbf{f}}_m({\mathbf{r}}^{\prime },t-\frac{R}{c_e})}{4\pi R}}\mathit{d}{s}^{\prime }.$	(47)

Next, the second term in Equation (45) is written as,

$⟨{\mathbf{f}}_m(\mathbf{r})g_n(t),\nabla \mathrm{\Psi }⟩=\nabla \left[{\displaystyle {\int }_S}{\displaystyle \frac{\nabla \cdot {\mathbf{f}}_m({\mathbf{r}}^{\prime },t-\frac{R}{c_e})}{4\pi R}}\right]d{s}^{\prime }.$

(48)

The third term in Equation (45) is:

$⟨{\mathbf{f}}_m(\mathbf{r})g_n(t),{\displaystyle \frac{\partial }{\partial t}}\left\{{\displaystyle \frac{1}{{ϵ}_i}}\nabla \times {\mathbf{F}}_i\right\}⟩$
$={\displaystyle \frac{\partial g_n}{\partial t}}{\displaystyle {\int }_S}{\displaystyle \frac{\nabla \times {\mathbf{a}}_n\times {\mathbf{f}}_m({\mathbf{r}}^{\prime },t-\frac{R}{c_e})}{4\pi R}}\mathit{d}{s}^{\prime }.$	(49)

It is easy to see that similar expressions are valid for Equation (46) with $c_e$ replaced by $c_i$. The integrals present in Equations (46) - (49) may be carried out analytically using the procedures developed for the triangular domains [17].

Using the standard expansion procedure for MOM problems, it is possible to generate a matrix equation $\mathbf{ZX}=\mathbf{Y}$ of dimension $P=N_t\times 2N_d$. The matrix equation can be efficiently solved using the special procedure developed in [12, 13].

IV. NUMERICAL RESULTS - VIE FORMULATION

In this section, we present numerical results for several inhomogeneous objects modeled by tetrahedral elements in the TD.

For all the examples presented in this section, and also for the next section, the following statements apply:

1. The object is placed at the center of the right-handed coordinate system, with the origin approximately coinciding with the geometrical center of the object. For all examples, $\theta$ and $\varphi$ represent the angles measured with respect to $z$ and $x$ axes, respectively.

2. The incident field is a GPW, given by:

$${\mathbf{E}}^{inc}(\mathbf{r}),t={\mathbf{E}}_o\frac{4}{T_P\sqrt{\pi }}{\text{e}}^{-{\gamma }^2},$$

(50)

where:

$$\gamma =\frac{4}{T_P}(ct-c{t}_o-\mathbf{r}\cdot \widehat{k}).$$

(51)

In Equations (50) and (51), $\widehat{k}$ is the unit vector in the direction of propagation of the incident wave, $T_P$ is the pulse width of the Gaussian impulse, ${\mathbf{E}}_o\cdot \widehat{k}=0$, $\mathbf{r}$ is a position vector relative to the origin, $c$ is the velocity of propagation in the external medium, and $t_o$ is a time delay which represents the time at which the pulse peaks at the origin. It may be noted that the GPW represents a smoothed impulse and as a result the response obtained may be considered as the impulse response. It is well-known that once the impulse response is available, the response to any other incident time domain wave form can be obtained by performing straight-forward convolution.

3. The incident plane wave is traveling along the $z$-direction, and the electric field vector is linearly polarized along the $x$-axis.

In the next subsection, we consider only homogeneous dielectric bodies and compare the results with the frequency domain (FD) solution. We note that to generate the TD data, the FD solution must be performed at several frequencies and then take the inverse Fourier transform into the time domain.

A. Homogeneous dielectric objects

As a first example, consider a dielectric sphere, radius$=0.1m$, ${ϵ}_r=3.0$, located with the origin coinciding with the center of the Cartesian coordinate system. The dielectric sphere is illuminated by the GPW described in Equation (50) with $T_p=4.0$ LM and $t_0=6.0$ LM, where the unit “LM” implies Light-Meter ($1$ LM = $3.333\times {1.0}^{-9}$ s). The sphere is modeled by 119 tetrahedrons resulting in 268 triangle faces implying that we have 268 space basis functions. We modeled the time variable with 36 triangle functions with $\mathrm{\Delta }t=2t_0/36=0.33$ LM. We note that the number of triangle functions for time are dictated by the pulse width of GPW and is not very critical. The numerical results are presented in Fig. 3. Here, we compare the RCS obtained as a function of $\theta$ with $\varphi =0$ by the present method with the frequency domain SIE solution at 100 MHz, 200 MHz and 300 MHz. We note that for the sphere problem, it is possible to generate exact solution using Mie series. However, the SIE formulation was well tested for canonical shapes [7] and for the sake of uniformness we compared our results with numerical MOM solution only. We note that both solutions compare very well at the selected frequencies.

Figure 3: RCS vs $\theta$ at $\varphi =0$ of a dielectric sphere (radius=0.1m and ${ϵ}_r=3.0$) illuminated by a GPW.

Figure 4: RCS vs $\theta$ at $\varphi =0$ of a dielectric slab ($1.0m\times 1.0m\times 0.1m$, and ${ϵ}_r=2.0$) illuminated by a GPW.

Next, we consider a rectangular homogeneous dielectric slab, with dimensions $1.0m\times 1.0m\times 0.2m$ and ${ϵ}_r=2.0$, located in the Cartesian coordinate system as shown in the Fig. 4. The dielectric slab is illuminated by the GPW described in Equation (50) with $T_p=8.0$ LM and $t_0=12.0$ LM. The slab is modeled by 626 tetrahedrons resulting in 1386 triangles. Thus, we have 1386 basis functions for space, and we modeled the time variable with 36 triangle functions. The numerical results are presented in Fig. 4. Here, we compare the RCS obtained by the present method with the SIE solution at 100 MHz, and 200 MHz. We note that both solutions compare very well for this case also.

Next, we consider a thick dielectric cylinder, with radius and height equal to $0.2m$, and ${ϵ}_r=3.0$, located with the origin coinciding with the center of the Cartesian coordinate system. The dielectric cylinder is illuminated by the GPW described in Equation (50) with $T_p=8.0$ LM and $t_0=12.0$ LM. The cylinder is modeled by 136 tetrahedrons resulting in 314 basis functions for space. We modeled the time variable with 36 triangle functions. The numerical results are presented in Fig. 5. Here, we plot the normal component of the electric flux density as a function of time and compare that with the solution obtained by the frequency domain MOM solution and Inverse Discrete Fourier Transform (VIE_FD_IDFT) method. The normal component is sampled at the center of the top face. We note that both solutions compare very well and also notice the absence of any late-time instabilities in the direct time domain solution. Although not shown here, the direct time domain solution was obtained up to $300$ LM and remainedstable.

Figure 5: Dielectric cylinder (radius=$0.2m$, height =$0.1m$, and ${ϵ}_r=3.0$) illuminated by a GPW.

Next, we consider a few examples of composite bodies where two or more homogeneous bodies are joined to form the inhomogeneous, composite body.

B. Objects with multiple dielectric materials

As a first example, consider a composite dielectric sphere, formed by combining two homogeneous dielectric hemispheres, each with radius=$0.2m$ and with distinct dielectric materials ${ϵ}_r=3.0$ and ${ϵ}_r=5.0$, as shown in the inset of Figs. 6 and 7. The composite dielectric sphere is illuminated by the GPW described in Equation (50) with $T_p=8.0$ LM and $t_0=12.0$ LM. Each hemisphere is modeled by 124 tetrahedrons resulting in 286 triangles. Therefore, the total number of basis functions for this case is 572, and we modeled the time variable with 36 triangle functions. In Fig. 6, we present the normal component of the electric flux density as a function of time. The normal component is sampled at $\theta ={45}^0$ and $\varphi ={90}^0$. The present solution is compared with the solution obtained in the FD, performing the calculations for 128 frequency points between 2- 256 MHz at 2 MHz interval, and performing IDFT. We note a good comparison between the two solutions and the absence of any late-time oscillations. In Fig. 7, we plot the backscattered field ($\theta ={180}^0$ and $\varphi =0^0$) as a function of frequency for both solutions. We note that both solutions compare very well for this case in the frequency range set by the pulse-width of the incident field.

Figure 6: Composite dielectric sphere illuminated by a GPW.

Figure 7: Backscattering RCS vs Frequency of a composite dielectric sphere illuminated by a GPW.

Next, we consider the case of a composite dielectric slab, formed by placing two slabs, one on top of the other slab, each with a distinct dielectric material as shown in the inset of Fig. 8. The dimensions of each slab are $1.0m\times 1.0m\times 0.2m$ and the dielectric constants are ${ϵ}_r=2.0$ and ${ϵ}_r=3.0$, respectively. The dielectric slab is illuminated by the GPW described in Equation (50) with $T_p=8.0$ LM and $t_0=12.0$ LM. Each slab is modeled by 626 tetrahedrons resulting in 1386 triangle faces. Therefore, the total number of basis functions for this problem is 2772, and we modeled the time variable with 36 triangle functions. The comparison solution is obtained by developing the frequency domain solution using VIE at 128 frequency points and then performing the IDFT to obtain the TD solution. The numerical results are presented in Fig. 8. Here, we plot the electric flux density as a function of time. The normal component is sampled at the center of the top face. We note that both solutions compare reasonably well for this case also, and the direct time domain solution remains stable even at a very late time. The apparent difference in the peak value is due to the small number of frequency samples available for the IDFT solution.

Figure 8: Composite dielectric slab illuminated by a GPW.

Next, we consider a composite dielectric cylinder, formed by joining two homogeneous cylinders. The height and radius of each cylinder is, $0.5m$ and $0.2m$, respectively. The dielectric constant for each cylinder is ${ϵ}_r=2.0$ and ${ϵ}_r=3.0$. The composite cylinder is located with the origin coinciding with the center of the Cartesian coordinate system and the axis coinciding with the $z$-axis. The whole body is illuminated by the GPW described in Equation (50) with $T_p=8.0$ LM and $t_0=12.0$ LM. Each cylinder is modeled by 410 tetrahedrons resulting in 902 triangle faces, and the total number of unknowns for this problem is 1804. We modeled the time variable with 36 triangle functions. The numerical results are presented in Fig. 9. Here, we plot the normal component of the electric flux density as a function of time and compare with the solution obtained using the frequency domain MOM solution and IDFT method. The normal component is sampled at the center of the top face. We note that both solutions compare very well and also notice the absence of any late-time instabilities. We also note that the IDFT solution behaved in a strange fashion (not tending to zero at late time) that can be attributed to a low sampling rate employed in thesolution.

Figure 9: Composite dielectric cylinder illuminated by a GPW.

Figure 10: Composite dielectric disk illuminated by a GPW.

As a last example in this subsection, we consider the case of three homogeneous dielectric discs, joined together as shown in the inset of Fig. 10. The radius and thickness of each disk is equal to $0.2m$. The dielectric constants are 3.0, 4.0, and 5.0 as shown in the figure. The whole body is illuminated by the GPW described in Equation (50) with $T_p=4.0$ LM and $t_0=6.0$ LM. Each cylinder is modeled by 136 tetrahedrons resulting in 314 triangle faces. The total number of basis functions for this problem is 942. We modeled the time variable with 24 triangle functions. The numerical results are presented in Fig. 10. Here, we plot the backscattered field as a function of time for both the frequency domain solution and the present algorithm. The back-scattered field data for the present time domain solution is obtained by performing a straightforward Fourier transform. We note that both solutions compare very well for this case in the frequency range set by the pulse-width of the incident field.

V. NUMERICAL RESULTS - SIE FORMULATION

In this section, we present numerical results for a few representative single/multiple homogeneous objects using direct time domain SIE formulation. We only consider canonical shapes in this work. For all the examples presented in this section, the incident field is a GPW, given by Equations (50) and (51). We present comparisons for the equivalent electric and magnetic currents, $\mathbf{J}$ and $\mathbf{M}$, respectively, with other methods. Also, note that for all results presented in this section, $\mathbf{M}$ is normalized with respect to the free space impedance, $\eta$.

A. Homogeneous dielectric objects

As a first example, consider a dielectric sphere, radius=$1.0m$, ${ϵ}_r=10.0$, located with the origin coinciding with the center of the Cartesian coordinate system. The dielectric sphere is illuminated by the GPW described in Equation (50) with $T_p=20.0$ LM and $t_0=30.0$ LM. The sphere is modeled by 288 triangles resulting in 432 edges implying that we have 864 unknowns (432 unknowns each for $\mathbf{J}$ and $\mathbf{M}$, respectively) for the solution scheme. We modeled the time variable with 45 triangle functions. The numerical results are presented in Fig. 11. In the figure, we present the induced equivalent current components, $J_x$ and $M_y$, at a selected point, ($\theta ={90}^{\circ }$ and $\varphi =0^{\circ }$), as a function of time. Note that we are presenting the results only for a single point on the body although the data is available for any point on the sphere. The time domain results are compared with the results obtained by IDFT solution. We further note that the IDFT solution is presented for a relatively shorter duration compared to direct TD solution (120 LM vs 600 LM). It is because the IDFT solution is periodic by nature since the inverse Fourier transform is performed and hence dictated by the frequency interval between two successive samples ($\mathrm{\Delta }f$). For a longer time signature, one must sample the frequency scale more closely which dramatically increases the computationaltime.

Figure 11: Dielectric sphere (radius=1.0m and ${ϵ}_r=10.0$) illuminated by a GPW.

Figure 12: Dielectric cube (a=2.0m and ${ϵ}_r=10.0$) illuminated by a GPW.

Figure 13: Dielectric disk (radius 1.0m, thickness 0.2m, and ${ϵ}_r=3.0$) illuminated by a GPW.

Figure 14: Dielectric cylinder (radius 0.2m, length 1m, and ${ϵ}_r=2.0$) illuminated by a GPW.

Next, we consider a dielectric cube of side length $2.0m$, and ${ϵ}_r=10.0$, located with the origin coinciding with the center of the Cartesian coordinate system. The dielectric cube is illuminated by the GPW described in Equation (50) with $T_p=20.0$ LM and $t_0=30.0$ LM. The cube is modeled by 432 triangles resulting in 648 edges implying that we have 1296 unknowns. We modeled the time variable with 45 triangle functions. The numerical results are presented in Fig. 12. In the figure, we present the induced equivalent currents, $J_x$ and $M_y$, at a selected point (at the center of the top face) as a function of time. The comparison between the two results is reasonably well. We note that the currents induced are very weak and the small differences we notice are due to the numerical errors. Further, we notice small oscillations for the IDFT scheme and absence of them in the TD scheme. This is due to the large time step in the TD scheme to cover the time scale 0–600 LM.

Next, we consider a circular dielectric disk with $1.0m$ radius, $0.2m$ thickness, and ${ϵ}_r=3.0$, located with the origin coinciding with the center of the Cartesian coordinate system. The disk is illuminated by the GPW described in Equation (50) with $T_p=20.0$ LM and $t_0=30.0$ LM. The body is modeled by 48 triangles resulting in 72 edges implying that we have 144 unknowns. We modeled the time variable with 45 triangle functions. The numerical results are presented in Fig. 13. In the figure, we present the induced equivalent currents, $J_x$ and $M_y$, at a selected point (at $x=0.0$, $y=0.5$ on the top face) as a function of time. We note that the results remain stable for a very long time whereas the IDFT solution is terminated much earlier. This is because, to obtain a long time signature using the IDFT solution, the frequency range needs to be densely sampled making the solution very expensive.

Next, we consider a dielectric cylinder of length $1.0m$, $0.2m$ radius, and ${ϵ}_r=2.0$, located with the center of the body coinciding with the center of the Cartesian coordinate system. The dielectric cylinder is illuminated by the GPW described in Equation (50) with $T_p=20.0$ LM and $t_0=30.0$ LM. The body is modeled by 176 triangles resulting in 264 edges implying that we have 528 unknowns. We modeled the time variable with 20 triangle functions. The numerical results are presented in Fig. 14. In the figure, we present the induced equivalent currents, $J_x$ and $M_y$, at a selected point (at $x=0.0$, $y=0.1$ on the top face) as a function of time. The time domain results are compared with the results obtained by the IDFT solution. Although the results compare well, there is some difference, particularly related to the electric current. We feel that the IDFT solution is showing oscillations because of the loss-less, perfect dielectric material used in the simulation. A small amount of loss would dampen the oscillations but may also contribute to some loss in the peak values.

B. Objects with multiple dielectric materials

Now, we consider the case of two dielectric cylinders placed along the $z$-axis and placed one cylinder on the top of the other cylinder as shown in the Figure 15. Each dielectric cylinder is of length $1.0m$ and $1.0m$ radius. The dielectric constant of the top cylinder is ${ϵ}_r=5.0$, where as the bottom cylinder’s dielectric constant is ${ϵ}_r=1.0$ implying that it is an air-dielectric body. Because of the air-dielectric nature, the result for a single cylinder and the combination of a dielectric cylinder with air-dielectric cylinder should be identical. The dielectric cylinder combination is illuminated by the GPW described in Equation (50) with $T_p=12.0$ LM and $t_0=18.0$ LM. Each body is modeled by 48 triangles resulting in 72 edges implying that we have 288 unknowns for the whole system (144 unknowns per cylinder). We modeled the time variable with 24 triangle functions. The numerical results are presented in Fig. 15. In the figure, we present the induced equivalent currents, $J_x$ and $M_y$, at a selected point (at $x=0.0$, $y=0.5$ on the top face of the top cylinder) as a function of time. As expected, the air-dielectric cylinder did not contribute to the scattering phenomenon.

Figure 15: A system of two dielectric cylinders, radius $1.0m$, length $1.0m$, illuminated by a GPW.

Next, we consider the case of two dielectric cylinders placed along the $z$-axis and touching each other as shown in the Fig. 16. Each dielectric cylinder is of length $1.0m$ and radius =$1.0m$. The dielectric constants of the top and bottom cylinders are ${ϵ}_r=3.0$, and ${ϵ}_r=5.0$, respectively. The dielectric cylinders are illuminated by the GPW described in Euation. (50) with $T_p=12.0$ LM and $t_0=18.0$ LM. Each body is modeled by 48 triangles resulting in 72 edges implying that we have 288 unknowns for the whole system (144 unknowns per cylinder). We modeled the time variable with 24 triangle functions. The numerical results are presented in Fig. 16. In the figure, we present the induced equivalent currents, $J_x$ and $M_y$, at a selected point (at $x=0.0$, $y=0.5$ on the top face of the top cylinder) as a function of time. The time domain results are compared with the results obtained by the IDFT solution and the comparison is reasonable.

Figure 16: Dielectric cylinders illuminated by a GPW.

Next, we consider the case of two dielectric cubes touching each other as shown in the Fig. 17. Each dielectric cube is of side length $1.0m$. The dielectric constants of the top and bottom cubes are ${ϵ}_r=3.0$, and ${ϵ}_r=5.0$, respectively. The two-body system is illuminated by the GPW described in Equation (50) with $T_p=12.0$ LM and $t_0=18.0$ LM. Each cube is modeled by 108 triangles resulting in 162 edges implying that we have 648 unknowns for the whole system (324 unknowns per cube). We modeled the time variable with 24 triangle functions. The numerical results are presented in Fig. 17. In the figure, we present the induced equivalent currents, $J_x$ and $M_y$, at the center of the top face as a function of time. The time domain results are compared with the results obtained by the IDFT solution and the comparison is reasonable.

Figure 17: Dielectric cubes illuminated by a GPW.

As a last example, we consider the case of five dielectric cylinders, each with $0.2$ m radius and $10$ m length, joined together to form a $50$ m long cylinder. The dielectric constant of each cylinder is 16, 14, 12, 8, and 4. The whole structure is placed along the $z$-axis and illuminated by a Gaussian pulse with $20$ LM pulse width. The time domain result is transformed into frequency domain and the RCS is compared with the direct frequency domain results at $10$, $20$ and $25$ MHz as shown in the Fig. 18. Each cylinder is modeled by 132 triangles resulting in 198 edges implying that we have 990 unknowns for the whole system. We note that the comparison is excellent for 10 and 20 MHz cases whereas we note a different result for the 25 MHz case. We attribute this difference to the bandwidth of the incident pulse which drops off steeply after 20 MHz. This example also highlights the limitation of the time domain solution. It is often said that one TD simulation is sufficient to obtain the frequency response from DC to daylight via Fourier transform. It is only true if one uses a true impulse for the incident field, which is not possible in our method. Hence, one should note that the frequency response that can be obtained is limited by the bandwidth of the incident pulse.

Figure 18: RCS of a 5-section, $50$ m inhomogeneous cylinder.

VI. CONCLUSIONS

In this work, we presented direct time domain formulations for dielectric bodies of arbitrary shape using surface and volume integral equation formulations. The main objective of this work is to demonstrate that the solution remains stable by eliminating time marching as was done in the previous works. Further, the formulation and the solution methodology is simple, not requiring any complex mathematical manipulations. Lastly, the present method can easily handle multiple right hand sides efficiently as required for monostatic radar cross section (RCS) studies, thus preserving the advantages of MOM solution scheme. Unfortunately, the conventional MOT scheme and all the DE methods including FDTD is not capable of performing this task efficiently and for each right hand side the solution must be started from the beginning.

ACKNOWLEDGEMENTS

This work is sponsored by the Office of Naval Research via the NRL Base Program.

REFERENCES

[1] K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Transactions on Antennas and Propagation, vol. 14, pp. 302–307, 1966.

[2] K. S. Yee, J. S. Chen and A. H. Chang, “Conformal finite-difference time domain (FDTD) with overlapping grids,” IEEE Transactions on Antennas and Propagation, vol. 40, pp. 1068–1075, 1992.

[3] R. Gomez, J. A. Morente and A. Salinas, “Time domain analysis of an array of straight-wire coupled antennas,” IEE Electronic Letters, vol. 22, pp. 316–318, 1986.

[4] R. Gomez, A. Salinas, A. R. Bretones, J. Fornieles and M. Martin, “Time domain integral equations for EMP analysis,” International Journal of Numerical Modeling, vol. 4, pp. 153–162, 1991.

[5] R. Gomez, A. Salinas and A. R. Bretones, “Time domain integral equation methods for transient analysis,” IEEE AP-S Magazine, vol. 34, pp. 15–22, 1992.

[6] S. M. Rao and D. R. Wilton, “Transient scattering by conducting surfaces of arbitrary shape,” IEEE Transactions on Antennas and Propagation, vol. 39, pp. 56–61, 1991.

[7] D. A. Vechinski, S. M. Rao, and T. K. Sarkar, “Transient scattering from three-dimensional arbitrarily shaped dielectric bodies,” Journal of the Optical Society of America, vol. 11, pp. 1458–1470, 1994.

[8] G. Manara, A. Monorchio, and R. Reggiannini, “A space-time discretization criterion for a stable time-marching solution of the electric field integral equation,” IEEE Transactions on Antennas and Propagation, vol. 45, no. 3, pp. 527–533, Mar. 1997.

[9] N. Gres, A. Ergin, E. Michielssen and B. Shanker, “Volume integral equation based electromagnetic scattering from three-dimensional inhomogeneous dielectric objects,” Radio Science, vol. 36, pp. 379–386, 2001.

[10] B. Shanker, K. Aygun, and E. Michielssen, “Fast analysis of transient scattering from lossy inhomogeneous dielectric bodies,”Radio Science, vol. 39, RS2007, doi:10.1029/2003RS002877, 2004.

[11] G. Kobidze, J. Gao, B. Shanker, and E. Michielssen, “A fast time domain integral equation based scheme for analyzing scattering from dispersive objects,” IEEE Transactions on Antennas and Propagation, vol. 53, pp. 1215–1226, 2005.

[12] S. M. Rao, “A simple and efficient method of moments solution procedure for solving time-domain integral equation - Application to wire-grid model of perfect conducting objects,” IEEE Journal on multiscale and multiphysics computational techniques, vol. 4, pp. 57–63, Mar. 2019.

[13] S. M. Rao “A straight-forward method of moments procedure to solve the time domain integral equation to PEC bodies via triangular patch modeling,” ACES Journal, vol. 35, pp. 843–854, Aug. 2020.

[14] R. Harrington, Time Harmonic Electromagnetic Fields, New York: IEEE Press, 2001.

[15] R. Harrington, Field Computation by Moment Methods, New York: Macmillan, 1968.

[16] S. M. Rao, Time Domain Electromagnetics, London: Academic Press, 1999.

[17] S. M. Rao, D. R. Wilton, and A. W. Glisson, “Electromagnetic scattering by surfaces of arbitrary shape,” IEEE Transactions on Antennas and Propagation, vol. 30, pp. 409–418, 1982.

[18] D. H. Schaubert, D. R. Wilton, and A. W. Glisson, “Tetrahedral modeling method for electromagnetic scattering by arbitrarily shaped inhomogeneous dielectric bodies,” IEEE Transactions on Antennas and Propagation, vol. 32, pp. 77–85, 1984.

BIOGRAPHIES

Sadasiva M. Rao received the Bachelor’s degree in electrical communication engineering from Osmania University in 1974, Master’s degree in microwave engineering from Indian Institute of Sciences in 1976, and Ph.D. degree with specialization in electromagnetic theory from University of Mississippi in 1980.

Dr. Rao served as an Assistant Professor in the Department of Electrical Engineering, Rochester Institute of Technology from 1980 to 1985, Senior Scientist at Osmania University from 1985 to 1987, and as a Professor in the Department of Electrical and Computer Engineering, Auburn University, from 1988 to 2009. He also held visiting Professorships at University of Houston (1987–1988), Osmania University, and Indian Institute of Science. Presently, he is with the Radar Division, Naval Research Laboratory, Washington, DC.

Dr. Rao worked extensively in the area of numerical modeling techniques as applied to Electromagnetic/Acoustic Scattering. He and his team at the University of Mississippi, were the original researchers to develop the planar triangular patch model and to solve the problem of EM scattering by arbitrary shaped conducting bodies. For this work, he received the best paper award for the period 1979–1981 from SUMMA Foundation. He published/presented over 150 papers in international journals/conferences. For his contributions in numerical electromagnetic problems, he was awarded the status of Fellow of IEEE. Further, he was recognized as a Highly Cited Researcher by Thomson ISI in 2001. Dr. Rao’s research interests are in the area of numerical methods applied to antennas and scattering.

ACES JOURNAL, Vol. 39, No. 3, 201–214
doi: 10.13052/2024.ACES.J.390305
© 2024 River Publishers

ABSTRACT

I. INTRODUCTION

II. VIE FORMULATION

A. Integral equation formulation

B. MOM solution

III. SIE FORMULATION

A. Integral equation formulation

B. MOM solution

IV. NUMERICAL RESULTS - VIE FORMULATION

A. Homogeneous dielectric objects

B. Objects with multiple dielectric materials

V. NUMERICAL RESULTS - SIE FORMULATION

A. Homogeneous dielectric objects

B. Objects with multiple dielectric materials

VI. CONCLUSIONS

ACKNOWLEDGEMENTS

REFERENCES

BIOGRAPHIES