Computational electromagnetics – A review

S. M. Rao and N. Balakrishnan*

Department of EE, Auburn University, Auburn, AL 36849, USA

*Department of Aerospace, Indian Institute of Science, Bangalore 560 012, India

In this paper, we describe the important computational methods available to successfully predict the radar cross-section of a given object. Further, we evaluate these methods based on the applicability to general structures, efficiency and accuracy. We confine our attention to the integral equation (IE) based methods although some passing references were made to the differential equation (DE) methods.

COMPUTATIONAL Electromagnetics (CEM) has evolved rapidly during the past decade to a point where extremely accurate predictions can be made for very general scattering and antenna structures. In general, all the available methods may be classified broadly into two categories, viz. a) differential equation (DE) solution methods and b) integral equation (IE) methods.

Although the Maxwell curl equations are usually first encountered in the time domain (TD), i.e. with time as an explicit, independent variable, until relatively recently, most electromagnetic instruction and research has taken place in the frequency domain (FD) where time-harmonic behaviour is assumed. A principal reason for favouring the FD over the TD in the pre-computer era had been that a FD approach was generally more tractable analytically. Furthermore, the experimental hardware available for making measurements in past years was largely confined to the FD.

The inferior position of TD electromagnetics (EM) began to change with the arrival of the digital computer, which has not only profoundly affected what can be done numerically (or computationally), but also experimentally. Since the beginning of what has come to be called computational electromagnetics (CEM) in the early 1960s, there has been a steady growth in both TD and FD modeling. This growth, which began slowly at first, was primarily confined to integral-equation (IE) treatments, but it has become almost explosive over the last 10 years as TD differential-equation (DE) modeling has attracted wide attention. This presentation summarizes the status of computational electromagnetic modeling and highlights some of the current research areas.

Modeling choices in CEM

In discussing CEM, it is appropriate to consider two basic questions:

  1. What are the various alternative modeling approaches available for CEM?
  2. What are the advantages of one model relative to the other possibilities?

 

To answer both questions, we observe that there are four major, first-principles, models in CEM, given by,

 

  1. Time Domain Differential Equation (TDDE) models, the use of which has increased tremendously over the past several years, primarily as a result of much larger and faster computers.
  2. Time Domain Integral Equation (TDIE) models, although available for well over 30 years, have gained increased attention in the last decade. The recent advances in this area make these methods very attractive for a large variety of applications.
  3. Frequency Domain Integral Equation (FDIE) models which remain the most widely studied and used models, as they were the first to receive detailed development.
  4. Frequency Domain Differential Equation (FDDE) models whose use has also increased considerably in recent years, although most work to date has emphasized low frequency applications.

 

These four choices can actually be narrowed down to two choices, i.e. a) IE models and b) DE models, depending on the mathematical formulation. The well-known method of moments (MoM) in general, involves IE modeling whereas the well-known finite element method (FEM) uses DE formulation.

General aspects of CEM modeling

The formulation and numerical development of a CEM model, in general, involves a number of basic steps whether a DE or an IE approach is being followed. In the following, we discuss some of these considerations.

 

 

Model development. For any numerical solution, it is necessary to develop the required equations and solve them on a computer. The equations thus developed must include the physics of the problem as well as the geometrical features. The following four steps are carried out in CEM problems:

  1. Develop integral equations using potential theory along with appropriate boundary conditions or alternatively, begin with the time-dependent Maxwell curl equations or their equivalent to develop methods such as FD, TD or FEM.
  2. Sample these equations in space, and also in time if it is a time-dependent equation, utilizing an appropriate geometrical space grid and suitable basis and testing functions. Note that, depending on the choice of formulation, the space grid may cover the structure and/or the surrounding space.
  3. Develop a set of simultaneous equations relating known and unknown quantities. Generally, the known and unknown quantities are the excitation field or its derivatives and the radiated/scattered field or induced current and charge, respectively.
  4. Generate a computer solution of this system in space and time as an initial-value problem.

 

Comparison of DE and IE models

The two primary choices for CEM modeling, as indicated so far, are those based on DEs and IEs. Whatever the details of the specific approach, any numerical method in its most general form, provides a way of solving IEs and DEs, involving approximating integrals as finite sums, and derivatives as finite differences in the generic forms

 

, (1)

 

leading, after some additional manipulation, to a linear system of equations, or ‘system’ matrix. The process of discretizing and quantifying DEs numerically is known by various names including finite-difference, finite-area (or volume), and finite-element procedures. The term finite-element is usually, but not necessarily, associated with a variational formulation, while use of a designation other than finite difference usually refers to the use of more general basis and testing functions. A numerical model based on an IE is also called a ‘boundary element’ method in structural dynamics and acoustics.

Most computer modeling involves replacing an infinite domain, first principles, analytical description of a problem by a finite domain, discretized, numerical one. The numerical model is finite in nature because only a limited number of unknowns of limited precision can be used in the solution process. An analytical model, however, entails, symbolically at least, an infinite dimensionality such as that exhibited by a series expansion for a sphere. However, it is worth noting that, from a practical viewpoint, the analytical model is finite in nature as well because the process of quantitatively evaluating any analytical model is automatically subject to limited precision and accuracy. This means that any observable of interest will exhibit no quantitative change over some specified dynamic range after an appropriate number of terms have been summed in its series solution. This will be the case whenever we deal with numerical answers as opposed to analytical solutions.

Some basic differences between DE and IE models are as follows:

 

  1. In general, the differential equation methods generate a sparse matrix, while the integral equation methods generate full matrices.
  2. Homogeneous/inhomogeneous/anisotropic materials can be handled in a relatively simple manner, while the level of complexity for the integral equation methods varies enormously for each of these cases.
  3. The code generation is straight forward for DE methods. This is usually not the case for integral equation methods.
  4. For DE methods, the solution space includes the object’s surroundings, the radiation condition is not enforced in exact sense, thus leading to certain error in the solution. For the IE solution, the solution space is confined to the object and the radiation condition is automatically enforced.
  5. The IE solutions are generally more accurate and efficient.
  6. Spurious solutions exist in DE methods whereas such solutions are absent in IE methods.
  7. For the DE solutions, developing numerical solution using parallel computer architecture is easy. However, for IE solutions, this is possible only in the time-domain. A lot of time and effort is needed to generate a parallel version for the frequency-domain IE solution.

 

Integral equation solutions in CEM

Mathematically speaking, an equation involving the integral of an unknown function of one or more variables is known as integral equation. One of the most common integral equations encountered in electrical engineering is the convolution integral given by

 

(2)

 

In eq. (2), we note that the response function Y(t) and the system function H(t, t ) is known and we need to determine the input X(t ). Of course, if X(t ) and H(t, t ) are known and we need to determine Y(t), then eq. (2) simply represents a integral relationship which can be performed in a straightforward manner. We further note that H(t, t ) is also commonly known as impulse response if eq. (2) represents the system response of a linear system. In general, in mathematics and in engineering literature, H(t, t ) is known as Green’s function or kernel function. We also acknowledge that, for some other physical systems Y(t) and X(t) may represent the driving force and response functions, respectively.

Next, we note that eq. (2) is known as integral equation of first kind. We also have another type of integral equation given by

 

(3)

 

where C1 and C2 are constants.

In eq. (3), we note that the unknown function X(t) appears both inside and outside the integral sign. Such equation is known as the integral equation of second kind. Further, we also see in electrical engineering yet another type of integral equation given by

 

(4)

 

which is known as integro-differential equation.

It may be noted that for a limited number of kernel and response functions, in eqs. (2–4), it is possible to obtain the solution using analytical methods. However, for a majority of practical problems, these equations can be solved using numerical methods only. Fortunately, in this day and age, we can obtain very accurate numerical solutions owing to the availability of fast digital computers. In the following section, we discuss a general numerical technique, popularly known as Method of Moments, to solve the integral equations (2)–(4).

Method of moments solution

The method of moments (MoM) solution procedure was first applied to electromagnetic scattering problems by Harrington1. Consider a linear operator equation given by

 

AX = Y, (5)

where A represents the integral operator, Y is the known excitation function and X is the unknown response function to be determined. Now, let X be represented by a set of known functions, termed as basis functions or expansion functions (p1, p2, p3,...,) in the domain of A as a linear combination:

(6)

where a i’s are scalar constants to be determined. Substituting eq. (6) into eq. (5), and using the linearity of A, we have

(7)

 

where the equality is usually approximate. Let (q1, q2, q3,...) define a set of testing functions in the range of A. Now, multiplying eq. (7) with each qj and using the linearity property of the inner product, we obtain

 

(8)

 

for j = 1, 2, ... , N. The set of linear equations represented by eq. (8) may be solved using simple matrix methods to obtain the unknown coefficients a i.

The simplicity of the method lies in choosing the proper set of expansion and testing functions to solve the problem at hand. Further, the method provides a most accurate result if properly applied. While applying the method of moments to complex practical problems, the solution region, in general, is divided into triangular subdomains, as shown in Figure 1. Then, one can define suitable basis and testing functions and develop a general algorithm to solve the electromagnetic problem2.

In Figure 2, we present the current induced on the aircraft shown in Figure 1 using the numerical procedure described in ref. 2, when illuminated by a 300 MHz electromagnetic plane wave. The plane wave is polarized along the length of the aircraft and travelling perpendicular to the body. Efficient solutions have been obtained for very complex problems using these methods in electromagnetics and acoustics3. It is possible that these methods found applications in other areas of engineering. Lastly, using the method of moments, solutions have been obtained for initial value problems also4.

Fast multipole method

One major problem with MoM is the generation of a dense matrix and for certain problems, the dimension of this matrix can be prohibitively large. Usually, for electromagnetic and acoustic scattering problems, it is necessary to divide the solution region into small enough

1345.gif (10940 bytes) 

1346.jpg (12353 bytes)

subdomains in order to obtain accurate results. By ‘small enough’, we mean about 200–300 subdomains per square wavelength. In usual practice, we may typically solve for several thousand unknowns for large, complex problems. Quickly, this requirement becomes expensive, in terms of computational resources, and may even become impossible to handle. Hence, we look for alternate schemes to reduce the computational resources.

The fast multiple method (FMM) dramatically reduces the time and memory required to compute radar cross sections and antenna radiation patterns compared to dense matrix techniques5. It is fairly simple to implement the FMM in a method of moments program to compute the electromagnetic scattering from large bodies of arbitrary shape. In the following, we describe the essential steps involved in the FMM implementation.

In the method of moments solutions to boundary integral equations, one is faced with solving large systems of equations of the form

 ZI = V, (9)

where Z is a dense matrix and I and V are column vectors of length N, where N is the number of current expansion functions. Eq. (9) can be solved by a number of iterative schemes. These techniques involve the computation of the product of Z and a solution vector one or more times for each iteration. This operation takes O(N2) operations and usually dominates all other operations in the iterative loop.

In the FMM, one groups the N basis functions into M groups so that the basis functions in each group have neighbouring support. For the simplest single stage FMM, the optimum value for M is proportional to ÷ N. Let the index m run over the groups and the index a refer to a basis function within a particular group. The dense matrix Z is then replaced with the expression

Z =  + VTV+, (10)

where , V, and T are sparse matrices. are those components of the original Z matrix for interactions between nearby regions of the target (typically within about one wavelength). The approximation can be made arbitrarily precise by the appropriate choice of FMM parameters in the computation of V and T.

The components of V are given by

 

(11)

 

where f are the basis functions. V is evaluated at K = ÷ N angles of k needed for a quadrature over the surface of a sphere.

 

The sparse matrix T is

(12)

where RmmĘ is the distance between the centers of the group m and . The number of terms in the sum, L, is chosen to give the desired accuracy in the FMM expansions.

Recently, the FMM algorithm was implemented on massively parallel architectures such as Intel Paragon. These machines typically consist of several hundred fast RISC microprocessors interconnected by a communication network. Let T(1) be the time required to compute ZI with one processor and T(P) be the time required with P processors. Ideally, one should have T(P)/T(1) = O(P). This is often difficult or impossible to achieve due to inherently sequential portions of the algorithm and communication costs.

For the FMM, the essential problem is to find an optimal distribution of the data structures. This is done by assigning one group to each processor.

A logical extension of FMM is the development of multilevel FMM algorithm. Further, we have yet another technique to improve the computational speed known as adaptive integral method (AIM). The details of these algorithms may be obtained from refs 6, 7.

Sparse matrix methods

The generation of a sparse matrix in the method of moment solution procedure may be achieved in two ways, viz. a) by defining a special set of basis functions to represent the unknown quantity or b) by handling the influence of the kernel function in a novel way. The usage of well-known wavelet-type basis functions to provide the required sparsity belong to the former category8 and the application of fast multipole method (FMM) belongs to the latter category5. So far, the wavelet-type basis functions have been applied to integral equations with one variable only and it remains to be seen how these functions can be utilized for two or more variable case. In contrast, in the FMM scheme, the matrix-vector product is carried out in a novel way and seem to work well for more complex problems.

There is also yet another scheme, known as impedance matrix localization scheme (IML) which achieves modest sparsity for simple problems9. Notice that the kernel function is, in general, a decaying function with respect to the distance between the source and observation points. Thus, with increasing distances, the influence of a given source becomes negligible at a sufficiently distant observation point and may be actually set to zero. The IML scheme cleverly exploits this fact. However, there is a certain degree of arbitrariness in this scheme and seems to work for simple problems only.

Recently a new method, known as generalized sparse matrix reduction scheme (GSMR), is proposed which seem to improve on the IML method10. The basic concept utilized in the GSMR technique may be qualitatively illustrated as follows. Following similar procedures of the MoM, a moment matrix is also generated in the GSMR method. However, in contrast to the conventional moment method where interaction is computed from each and every cell on other cells, only the interaction from the self-cell and few neighbouring cells is computed in the GSMR technique. In fact, for single variable problems (wire scatterers and two-dimensional infinite cylinders) only the self-term and two neighbouring terms on either side of the self-cell are generated in this technique. However, this technique, although appears promising, needs to be validated for more complex geometries.

Conclusions

In this paper, we have described the use of integral equation based methods in computational electromagnetics. In general, these methods are applicable to scatterers whose characteristic dimensions are of the order of a wavelength. The iterative methods and FMM are then called in to extend the IE methods to scatterers of larger dimension. Since the integral equation methods are global in nature, these methods work very well for scatterers with smooth geometries. The local nature of the DE formulation may be exploited to cater for the sharp variations, thus creating hybrid techniques. Both DE and IE, even with the present day super computers, are suited for computing the RCS of a fighter aircraft at best up to 1 GHz. For higher frequencies, the recourse is often taken to asymptotic methods like the geometrical theory of diffraction (GTD), uniform theory of diffraction (UTD), and uniform asymptotic theory (UAT). In essence, the CEM has grown today to be a mature tool for predicting precisely the scattering and radiation characteristics of most complex structures encountered in real life from human body to aircraft.


  1. Harrington, R. F., Field Computation by Moment Methods, Macmillan, New York, 1968.
  2. Rao, S. M., Wilton, D. R. and Glisson, A. W., IEEE Trans. Antenna Propagation, 1982, 30, 409–418.
  3. Miller, E. K., Medgyesi–Mitschang, L. and Newman, E. H., Computational Electromagnetics – Frequency-Domain Method of Moments, IEEE Press, New York, 1992.
  4. Rao, S. M., Time Domain Electromagnetics, Academic Press, New York, 1999.
  5. Coifman, R., Rokhlin, V. and Wandzura, S., IEEE AP-S Mag., 1993, 35, 7–12.
  6. Song, J. M. and Chew, W. C., Microwave Optical Technol. Lett., 1995, 10, 14–19.
  7. Bleszynski, E., Bleszynski, M. and Jaroszewicz, T., IEEE AP-S International Symposium Digest. Seattle, WA, June 1994, pp. 416–419.
  8. Steinberg, B. Z. and Leviatan, Y., IEEE Trans. Antenna Propagation, 1993, 41, 610–619.
  9. Canning, F. X., IEEE Trans. Antenna Propagation, 1993, 41, 659–667.
  10. Rao, S. M. and Gothard, G. K., Microwave Optical Technol. Lett., 1998, 19, 271–274.