Abstract
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%.
1. 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–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–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 higher-order 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 two-dimensional 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 two-dimensional 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.
2. 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 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.
2.1. 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 are unit vectors. And there are
When the interface is the isosurface (as shown in Figure 1), the no-reflect condition is , then (1) is rewritten as follows:
UPML equations in an orthogonal curvilinear coordinate system are as follows: where and are electromagnetic fields; and are permittivity and permeability, respectively. a, b, and c are diagonal tensors. and are auxiliary variables, which are satisfied with the following relationship: where d and are the same as a, b, and c. Their expressions are as follows:
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 are components of unit vectors , , and 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 is not the diagonal matrix anymore.
2.2. 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 where , , and c is the focal distance. Its trail is shown in Figure 2. 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 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:
Thus, the transformation matrix in (7) can be written as follows:
Substituting (17) into (8), then
UPML are distributed between the regions . in elliptical cylinder coordinate systems; thus, (2) is rewritten as follows: where where is the isosurface of the ellipse in the barycentric coordinates of the element in UPML region. and are the isosurface of the inside and the outside, respectively. m is taken as an integer.
Supposing that , , and focal distance c are known, we need to calculate the parameters of a, b, and ; the equation of the ellipse is
Thus, the formulas of semimajor axis a and semiminor axis b are where .
Substituting (22) into (11), we can get
Inserting (23) into (20), we can obtain the value of conductivity .
3. Iterative Formulas of 2-D DGTD in Elliptical Cylinder UPML
In the case of two-dimensional TM wave, and . Substituting (18) into (4), we can get
The auxiliary equations in a two-dimensional case are
There are nondiagonal elements in (24) and (25), such as , 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, and in (24) need to be eliminated. Suppose that Thus
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
The iterative formulas of auxiliary equations are where , , , , , , , , , and are iterative coefficients.
4. 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 . The radiant source of Gaussian impulse which has the width and peak of and 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, 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.
(a) Time-domain waveform and spectrum
(b) Reflection error
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(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 divided into 94,634 triangular units and 47,662 nodes with the discrete scale of 0.08 mm. The incident plane wave has a wavelength of 3 mm. Bistatic RCS for coated dielectric elliptical column target is shown in Figure 5(b) (“o”). The results of FDTD (solid line) and DGTD (“▲”) (its region is truncated with rectangular PML) are also shows in the figure. As can be seen, the calculation results of the three coincide.
(a) Sparse discrete model
(b) Bistatic RCS for two nested elliptical columns
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.
5. 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.
Conflicts of Interest
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research is supported in part by the National Natural Science Foundation of China under Grant nos. 61571348, 61231003, and 61401344.