Implementation of an Approximate Conformal UPML in 2-D DGTD

Using the numerical discrete technique with unstructured grids, conformal perfectly matched layer (PML) absorbing boundary in the discontinuous Galerkin time-domain (DGTD) can be set flexibly so as to save lots of computing resources. Based on the DGTD equations in an orthogonal curvilinear coordinate system, the processes of parameter transformation for 2-D UPML between the coordinate systems of elliptical and Cartesian are given; and the expressions of transition matrix are derived. The calculation scheme of conductivity distribution in elliptic cylinder absorbing layer is given, and the calculation coefficient of DGTD in elliptic UPML is calculated. Furthermore, the 2-D iterative formulas of DGTD and that of auxiliary equation in the elliptical cylinder UPML are derived; the conformal UPML calculation in DGTD is realized. Numerical results show that very good accuracy and computational efficiency are achieved by using the method in this paper. Compared to the rectangular computational region, both the memory and computation time of conformal UPML absorbing boundary are reduced by more than 20%.


Introduction
Truncated boundary condition is the key to ensure the accuracy for many electromagnetic numerical methods.In the past decades, different kinds of absorbing boundary conditions (ABC) have been proposed and successfully applied to the finite difference time-domain (FDTD) method and other methods.The perfectly matched layer (PML) ABC is presented by Berenger in 1994 [1].Later, the coordinate stretched perfectly matched layer (CPML) ABC was proposed by Chew and Weedon [2], and uniaxial anisotropic perfectly matched layer (UPML) ABC was presented by Sacks et al. [3].The conformal PML was analytically derived and applied to the hyperbolic model by Teixeira et al. in the 1990s [4][5][6].The corner-free truncation strategy was proposed and used in FDTD method by Zhang et al. [7].The discontinuous Galerkin time-domain (DGTD) [8][9][10][11][12] method, which was developed based on the finite volume element time-domain method, has advantages of mesh flexibility of finite element time-domain (FETD) and explicit iterative of FDTD.As a new kind of time-domain algorithm, the study of ABC for DGTD is a hot topic in the computational electromagnetic field.
The first-order Silver-Muller (SM) ABC is widely used owing to its easy realization [13].However, this kind of ABC only has a better absorption when the wave is perpendicular to the boundary.In order to avoid reflection, it is necessary to increase the distance between the scattering target and the absorption boundary; thus, the calculation amount is increased.Lu et al. introduced nonconforming UPML theory into DG calculations in the case of TM wave and compared the local relative error between PML and first-order SM-ABC [14].König et al. applied the stretched-coordinate PML technique to DGTD [15].Gedney and Zhao proposed a kind of auxiliary differential equation (ADE) multilevel PML theory by complex frequency domain shift, which is used to cut off the boundary of DGTD region and take good absorption [16].In the literature [17], the PML theory is used in 3-D DGTD algorithm to realize the absorption of electromagnetic wave through the conformal truncation face by Dosopoulos and Lee.Donderici and Teixeira incorporated conformal PML into the mixed FETD algorithm by utilizing PML constitutive tensors in 2008 [18].Based on a higherorder curvilinear finite-element and the concept of transformation electromagnetics, Smull et al. proposed anisotropic locally conformal PML; this new PML has favorable accuracy and efficiency [19].Modave et al. designed PMLs for transient acoustic wave propagation in generally shaped convex truncated domains in DGTD [20].Yang et al. applied the corner-free truncation strategy to DGTD calculation, which saved the resource effectively [21].Based on well-posed PML theory, Ren et al. proposed multiaxial full anisotropic media PML for subdomain and applied it to electromagnetic simulation of nonconformal mesh bicrystals; this method avoids the potential late-time instability found in classical PML [22].However, many physical problems can be solved successfully by two-dimensional models.In twodimensional models, it has advantages of small memory, fast calculation speed, easy handling the electromagnetic problem of electric large size by computer, and so on in three-dimensional models.Nevertheless, there is only little research focus on the study of DGTD method in twodimensional problems.Implementation and application of conformal UPML in 2-D are not yet reported.In this paper, the 2-D DGTD iterative formula for elliptic cylindrical UPML is derived and applied to the nested dielectric column scattering.
In order to apply the UPML technique to DGTD in two-dimensional cases, the wave equations of UPML in an orthogonal curvilinear coordinate system are considered with the theory of parameter transformation, the processes of parameters transformation for 2-D UPML, and the specific expressions of transition matrix between the coordinate systems of elliptical and Cartesian are derived.The penalty flux is used for field exchange between units.The 2-D DGTD iterative formula and auxiliary equation in elliptical cylindrical UPML are derived.Numerical results show that very good effectiveness is achieved by using the proposed algorithm.

The Expression of UPML Equations from Orthogonal Curvilinear Coordinate System to Cartesian Coordinate System
Suppose that the interface from free space to PML is the isosurface ξ 3 in the orthogonal curve coordinate system (as shown in Figure 1); if the parameters in UPML satisfy the similar condition in Cartesian coordinate system UPML, the incident electromagnetic waves from free space will pass through the interface without reflecting and attenuate rapidly inside the PML.

UPML Wave Equation in Orthogonal Curvilinear
Coordinate System.From the literature [17,23], we know that the form of diagonal matrix of the parameters in the orthogonal curve coordinate system is unchanged.So UPML matching matrix is written as follows: where u 1 , u 2 , u 3 are unit vectors.And there are When the interface is the isosurface ξ 3 (as shown in Figure 1), the no-reflect condition is s 1 = s 2 = 1, then (1) is rewritten as follows: UPML equations in an orthogonal curvilinear coordinate system are as follows: where E and H are electromagnetic fields; μ 1 and ε 1 are permittivity and permeability, respectively.a, b, and c are diagonal tensors.P h and P e are auxiliary variables, which are satisfied with the following relationship: where d and κ −1 are the same as a, b, and c.Their expressions are as follows: 2 International Journal of Antennas and Propagation In practical applications, ( 4) and ( 5) often need to be solved in Cartesian coordinates, thus, they are transformed from a curvilinear coordinate system by transformation matrix, which can be written as follows: where u 1x , u 2x , u 3x , … are components of unit vectors u 1 , u 2 , and u 3 in x, y, and z direction, respectively.Transformation matrix has a unitary nature [24].
Taking the example of tensor a, the relationship of a from a curvilinear coordinate system to Cartesian coordinates is as follows: It is noted that after the transformation, the a xyz is not the diagonal matrix anymore.

The Parameter Calculation for Elliptic Cylinder UPML.
Since the expression of the transformation matrix depends on the specific coordinate system, the process of calculational parameters in the elliptical cylinder UPML can be described as follows.
Let the ellipse cylindrical axis be z-axis, and the elliptic equation is as follows: where a and b are the semimajor axis and semiminor axis, respectively.The variables w and t are introduced, and let

11
where w ≥ 0, 0 ≤ t ≤ 2π, and c is the focal distance.Its trail is shown in Figure 2. w represents the distance sum between focal points in the ellipse.The transformation matrix between local coordinate of elliptical system and Cartesian coordinate system is decided by rotating angle β (as shown in Figure 2), which is defined as included angle between normal vector ŵ and the x-axis in the Cartesian coordinate system.From [25], we know that the tangent line equation of the observation point P x P , y P on the ellipse is as follows: The slope of the tangent line k is as follows: where α is the included angle between tangent line and normal vector ŵ.Then, the direction cosine and direction sine between tangential unit vectors and normal vector ŵ are as follows: The β is as follows: Then, the direction cosine and direction sine between normal unit vectors and x-and y-axes are as follows: UPML are distributed between the regions w 0 ≤ w ≤ w max .ξ 3 = w in elliptical cylinder coordinate systems; thus, (2) is rewritten as follows: where w is the isosurface of the ellipse in the barycentric coordinates of the element in UPML region.w 0 and w max are the isosurface of the inside and the outside, respectively.m is taken as an integer.Supposing that x p , y p , and focal distance c are known, we need to calculate the parameters of a, b, and w; the equation of the ellipse is Thus, the formulas of semimajor axis a and semiminor axis b are

Iterative Formulas of 2-D DGTD in Elliptical
Cylinder UPML In the case of two-dimensional TM wave, E x = E y = 0 and H z = 0. Substituting ( 18) into (4), we can get There are nondiagonal elements in ( 24) and ( 25), such as a xy , a yx , b xy , b yx , and so on.This will lead to a complex form of matrix equations.In order to facilitate the acquisition of the time-domain iterative formula, ∂H y /∂t and ∂H x /∂t in (24) need to be eliminated.Suppose that After using the penalty flux [11] as well as expending it by basic functions in (27), we can get After ( 28) is discrete in time axis, then there are

z,n 29
The iterative formulas of auxiliary equations are

Numerical Examples
In the elliptical cylindrical computing region, the semimajor and semishort axes of the outer PML boundary are 4 m and 3 m, respectively.And the inner PML boundary has a semimajor axis of 3.5 m.Both the inner and outer PML boundaries have the same focus length of 2.29 m.The region is divided into 46,122 triangular units and 23,298 nodes with the discrete scale of 0.05 m.The time discrete interval is taken as dt = 0 16 × 10 −10 s.The radiant source of Gaussian impulse which has the width and peak of τ = 200 dt and t 0 = 100 dt is set at the center of the region (0 m, 0 m).Time-domain waveform and spectrum at the monitoring point of (0, 2.5 m) are shown in Figure 3(a).The solid and dotted line are the impulse waveform in time-domain and the spectrum, 5 International Journal of Antennas and Propagation respectively.And "o" is the analytic solution obtained by the second kind of zero order Hankel function.Refer to [18]; the resulting reflection error is shown in Figure 3(b).It can be seen from the figure that the reflection error is below −37 dB.
The PEC elliptic column has a semimajor axis and a semishort axis of 1.5 mm and 0.5 mm, respectively.Incident plane wave has a wavelength of 2 mm.The computational region is truncated with conformal PML layer.The result of bistatic scattering for metal elliptical column is shown in Figure 4 ("o").For easy comparison, the result of method of moment (solid line) is also shown in figure.As illustrated, the results are consistent.
Two nested dielectric elliptical columns target truncated by conformal PML.Sparse discrete conformal model is shown in Figure 5     As can be seen, the calculation results of the three coincide.
In order to verify the benefits of conformal UPML, the numbers of unit and node in two different regions by three discrete scales of 0.1 mm, 0.08 mm, and 0.04 mm are given in Table 1.Compared to the elliptical region, their numbers in rectangular computational region (the region surrounded by rectangular boundary which is shown in Figure 5(a)) increase by about 27%.In order to further compare the time and memory consumption, the following calculation is given.Take the case of discrete scale of 0.08 mm example (shown in Table 1), the numbers of unit and node in rectangular region are 121,474 and 61,175, respectively, which increase by 28.3%.And the memory and time increase by 25.5% and 32.5%, respectively.It is clear that conformal UPML has obvious advantages in terms of computational resources and computational time savings.

Conclusions
According to the DGTD equations and the parameter transformation theory in orthogonal curvilinear coordinate systems, the processes of parameter transformation for 2-D UPML between the coordinate systems of elliptical and Cartesian are presented, and expressions of transition matrix are derived.The 2-D iterative formulas of DGTD and that of auxiliary equation in elliptical cylinder UPML are derived, and the conformal UPML calculation in DGTD are realized.Numerical results show that good absorption for the outgoing wave can be realized by using conformal UPML.Because of the discrete characteristics of the unstructured grids, the absorption layer can be set up flexibly according to the shape of the target.Compared with the computational cost in rectangular PML layer, the memory and computation time for conformal structure decrease by 20.3% and 24.5%, respectively.Thus, the method presented in this paper will be very beneficial to the engineering computation since lots of computing resources can be saved.
(a).The inner elliptical target has a relative dielectric coefficient of 4, and the semimajor axis and semiminor axis are 5 mm and 2 mm, respectively.The outer has a relative dielectric coefficient of 2, which are 6 mm and 3 mm.The rest of conformal boundaries from inside to outside are connected boundary, extrapolated boundary, and PML boundary.Both the inner and outer PML boundaries have the same focus distance.The calculating region is

Figure 3 :
Figure 3: Radiant field and reflection error.
Bistatic RCS for two nested elliptical columns

Table 1 :
Comparison of the number of unit and node in different calculated regions.