Late-time Instability for UPML and Periodic Boundary in Elongated Multilayer Thin Plate Simulations

In this article, we examine the late-time instability properties of hybrid boundary conditions in the discontinuous Galerkin time-domain (DGTD) simulations of an elongated multilayer thin plate. The hybrid boundary is combined by uniaxial perfectly matched layer (UPML) and periodic boundary condition (PBC). Herein, the PBC is employed to approximate an infinite long target. For the target studied, when implementing the UPML within the discrete DGTD domain, late-time instabilities would occur. This instable or spurious information can severely corrupt the solution of a physical problem in time domain. To suppress them, two effective ways are proposed, i.e., increasing the size of the air space (the distance away from the interface between the target studied and the UPML) and decreasing the conductivity of the UPML. The numerical experiments verify that the instability characteristics can be efficiently attenuated by two methods proposed in this paper.


Introduction
To mimic unbounded regions in time domain calculation of Maxwell's equations, an absorbing boundary condition (ABC) must be introduced to truncate the computational region.A popular kind of artificial absorbing boundary treatment is the perfectly matched layer (PML) invented by Berenger [1] in which each field of Maxwell's equations is split into two orthogonal tensor forms.The innovation of Berenger's PML is that plane waves of arbitrary incidence, polarization, and frequency could be absorbed perfectly and without reflections into the computational domain [2].However, due to the field splitting, the memory requirements will double.Fortunately, a lossy UPML suggested by Sacks and Ziolkowski does not require the splitting of the Maxwell's fields [3], [4].The formulations of the well-posed UPML are identical to the Maxwell's equations [5], therefore, only a few relatively straightforward complex-coordinate modifications are needed to implement in existing numerical methods [6].The lossy UPML employs the split conductivity to realize the anisotropic property of the PML and achieve the decay of the fields inside the UPML region [7].The UPML is a Maxwellian formula-based perfectly matched absorbing medium, it offers a number of significant advantages: more computational efficient, more flexible implementation, lower buffer spaces, and easier to extend to unstructured grid techniques.
DGTD is an excellent unstructured-grid-based numerical method [8], which is originated from the finite volume time domain (FVTD) and the finite element method (FEM) [9], [10].It retains spatial high-order convergence and adaptive unstructured meshes of the FEM and explicit iterative and easy parallelization of the finite difference time domain (FDTD) [2].Thus, the DGTD approach is more accurate, stable, flexible, and efficient than FDTD, FVTD and FEM.What's more, DGTD employs a Galerkin test procedure for each element to obtain the spatial discretization, and applies a unique numerical flux to provide the coupling information from neighbor elements [11], [12].It is a very powerful full-wave numerical approach to solve all kinds of time-dependent electromagnetic problems, especially suitable for handling radiation, scattering, and propagation problems which involve open regions.In an open region computations, it is common to surround artificial boundaries of a computational domain with the UPML of finite thickness to prevent artificially reflected waves from contaminating a numerical simulation [13].Theoretically, the UPML can ensure zero reflection from an interface with the computational domain.However, in the discretized implementation of the UPML equations some unwanted features appear, e.g., frequency-dependent reflection and even divergence [15].Using UPML ABC, a gradual long-term growth of the solution will destroy the accuracy of a numerical simulation everywhere.Eventually, the DGTD simulations would be subjected to some issues of late-time instabilities or spurious results [13], [18].This is mainly because waves inherently have long-term interactions with the UPML boundary or an irregularly shaped computational domain is inherently ill-posed.
In this paper, the issues pertaining to variational spatial distance and UPML conductivity are systematically studied, to verify how they affect the late-time stability properties of the hybrid boundary conditions.Herein, the high-order discontinuous Galerkin scheme in space accompanied by the explicit 4 th -order Runge-Kutta evolution in time is adopted.The PBC is employed to approximate an infinite long target and the UPML is applied to truncate the infinite width computational domain.The elongated multilayer thin plate could be regarded as a large aspectration target which is infinitely long in the y-axis direction but extremely thin in the x-axis direction, that is, the dimension of the computational domain in the x-axis direction is not sufficiently large compared to the one in the y-axis directions.For this special computational model, the existence of non-physical higher-order modes may cause instable or spurious characteristic [7].There are two ways to improve the late-time instability: one is to increase dramatically the size of the air space w to postpone the unstable wave shape; the other is by adjusting the value of the polynomial-graded conductivity σ to achieve a stable solution.

Implementation of UPML in DGTD
When one maps the UPML into the discrete DGTD space, to avoid the convolution-type constitutive relation in the time domain, an auxiliary differential equation (ADE) approach is used to acquire an efficient formulation [2].Consider the normalized time-domain form of Maxwell's equations for TM z (in 2D case, H z = 0 and E x = E y = 0) waves in physical and UPML regions [4], [5]  where E z is the z-component of the normalized electric field, H x and H y are the x-and y-component of the normalized magnetic field, respectively. r and  r are separately the relative permeability and permittivity.Q x , Q y and P z are auxiliary variables introduced by the UPML.These auxiliary fields exist only in the UPML region and thus the additional computational cost is little.σ x and σ y are the normalized relative conductivity for the UPML in the xand y-direction.Due to the discrete approximation of electric and magnetic fields, the material parameters at the interface of the UPML will result in a spurious impedance loading.Usually, the conductivity σ is set to a polynomialgraded profile along the normal axes of the UPML [2], [20], [21], as shown , , where σ i,max is the maximum conductivity in i-axis direction, d is the thickness of the UPML, l is the distance from the air-UPML interface to the internal UPML, which is varied in the range [0, d] and m is order of the polynomial variation, which represents a finer spatial discretization.As far as (2), it is shown that by choosing appropriate profile one can control the feature of the PML reflectance decrease with the increase of the thickness d [14], [15].
Supposing the computational domain is decomposed into K non-overlapping triangles in 2D space.For given an arbitrary element D k , the unknown fields can be expanded by interpolating multivariate Lagrange polynomial where N p = (N + 1)(N + 2)/2 stands for the minimum number of nodal points, N signifies the maximum order of the polynomial, and x = (x, y) is the position vector.On account of the fact that correctly choosing interpolation nodes can bring about good numerical behaviors, this work employs the Legendre-Gauss-Lobatto (LGL) interpolating nodes as x i [8].
Choose the same multivariate Lagrange interpolation polynomials as test functions.Multiplied by the test functions, integrated over each element, and followed by integration by parts, the fully explicit semi-discrete scheme are obtained as Here, k D  is the boundary of D and n denotes the local outward pointing normal vector.The matrices J, D, and M represent the local transformation Jacobian, the differentiation matrix and the mass-integration matrix, respectively (see Ref. [8] for details).Using of the notation x y z H ,H ,E  q and d     q q q are simple notations.
To facilitate the coupling with the neighbor elements, we use a pure upwind flux which could strongly damp unphysical modes [8], as follows,  (5) where Z ± and Y ± are the local impedance and admittance, respectively, and defined as 1 The characteristic impedance and admittance of free space is given by . The superscript minus signs refer to the local element and the superscript plus signs to its adjacent element across the interface [5].The numerical fluxes play a role in performing the communication be-tween adjacent elements.For the special case of a PBC boundary, the boundary fluxes can be acquired by setting Y    because the PBC boundary behaves as a material with an infinite admittance.

Computational Setup of Model
Assume the elongated multilayer thin plate is composed of two skin layers and one core layer.For the sake of simplicity, the skin thickness is set to 0.1λ (where the wavelength λ = 10 mm) and the core thickness is set to 1.4λ, thus the total thickness of the thin plate is 1.6λ.The thin plate is assumed to be lossless and has a dielectric constant of ε r,skin = 3.3 and ε r,core = 1.08.The origin is assigned to the center of the bottom edge of the thin plate.The air space which is signified by w is variable as illustrated in Fig. 1.Furthermore suppose a height of 1.0λ small unit is backed by the PBC boundaries in the y-direction, which is to mimic an infinite long shape.To truncate the computational domain in the x-direction we add two additional layers whose thickness is 1.0λ, symmetrically located at two sides of the thin plate.Moreover, the UPML is terminated by PEC boundaries, to enable propagating waves are sufficiently suppressed, resulting to provide just an insignificant contribution to the studied EM fields.The source is located at (-1.5λ, 0.5λ) and the monitoring point is located at (1.5λ, 0.5λ).
The UPML conductivity σ i has the polynomial-graded profile as given in (2).However, there are the PBC boundaries in the vertical direction of our model, thus (2) is modified as  x (x) =  x,max (l/d) m ,  y (y) = 0.In all subsequent computations we use this absorption profile.For the order of the profile m, typical values are 2～4 and we shall choose m = 3 unless stated otherwise.The computational model is illuminated by a vertically directed electric current source with a Gaussian time-signature.In all computations x 0 = -1.5λ,y 0 = 0.5λ, δ x = 0.03λ，δ y = 0.06λ.The initial magnetic field components are zero, as shown, , ,0 0, , ,0 0, , ,0 exp . 2 2 x y z x y

H x y H x y x x y y E x y
To calculate the scattering parameters on the elongated multilayer thin plate, DGTD with 3 rd -order basis functions is applied in space and the 4 th -order low-storage explicit Runge-Kutta (LSERK) scheme in time is used.
Here, taking the air space w = 0.8λ, the UPML conductivity σ x,max = 50, and FinalTime = 500, i.e., a total of 224,317 steps of numerical calculations is performed.The shorttime and long-time propagation behavior are illustrated in Fig. 2. A view of the late-time response of E z field at the observation point in Fig. 2(c) reveals any instable or spurious value once generated, no matter how small, will grow up to be significant in the very late time.The gradual values will pollute the entire computational domain and lead to late-time instability [22].

Mitigate Late-Time Instability
In this section, we examine the issues pertaining to the air space w and the UPML conductivity σ x,max , to see how they impact the late-time instability property for the elongated multilayer thin plate with UPML and PBC boundaries.

Enlarge the Air Space
The existence of non-physical higher-order modes generated in the elongated domain may cause a late-time instability characteristic [7].In this case, we fix the UPML conductivity σ x,max , i.e., taking σ x,max = 50.And let the width of the air space w to be the only free parameter, i.e., from 0.8λ to 5.2λ.The computational domain is meshed with unstructured triangular elements, and much finer triangles are used near the thin plate to better approximate the thin plate surface.The shortest edge of element is 0.05λ.According to the stability conditions of the DGTD method, the smallest time discrete interval t = 0.0031.After extensive trials, it is interesting to observe that increasing the width of the air space w dramatically improves the late-time instability.In physics and engineering, the envelope of an oscillating signal is a smooth curve outlining its extremes [23].Owing to the electric field at the monitoring point E ztest is an approximate oscillating signal, we employ its root mean square (RMS) of the upper envelope to represent the instable characteristic of the long-time propagation behavior, as shown in Fig. 3. From Fig. 3 we can see that when w  1.2λ, solutions of a longtime interaction become unstable, i.e., narrower w is easier to render significant late-time instability.Compared to the case of w = 0.8λ, if the UPML is moved to 2.2λ, the upper envelope of E ztest will drop by about 5:1.When w  3.2λ, we assert the solution to be stable even up to about 230,000 time steps.

Decrease the UPML Conductivity
Owing to increasing the width of the air space shall aggravate cost in system memory and computational time.To still generate a stable solution, we fix the width of the air space, i.e., taking w = 0.8λ and change the value of σ x,max in the UPML formula from 50 to 1.The RMS of the upper envelope of E z at the monitoring point is presented in Fig. 4. It is interesting to observe that the larger the conductivity σ x,max , the worse the solution after long time running, that is, a larger value of σ x,max will bring about more significant late-time instability.By inspection of Fig. 4 it is evident that taking 5  σ x,max  15 would prevent the late-time instability.However, using a too small value of conductivity i.e., σ x,max = 1, the mitigation of the latetime instability is a little poorly effective, comparing with σ x,max = 10 and σ x,max = 5.After extensive trials, it is interesting to observe that the first maximum of E ztest 0.06757 is obtained at Timestep = 794.Changing the UPML conductivity σ x,max from 50 to 1, the different effects relative to the different σ x,max are shown in Tab. 2. One can clearly see the late-time instability problem can be mitigated by using the smaller value of σ x,max , especially, when σ x,max = 10, the RMS of the E ztest reduces to about 3% relative to the first maximum value.

Compare the Two Strategies
Figure 5 shows the RMS of the upper envelope of long time response of the electric field E ztest for a small σ x,max value and a large w.In this experiment, for the magenta solid line, the value of σ x,max is fixed, e.g., taking σ x,max = 100 and the larger width of the air space, i.e., taking w = 5.2λ; for the blue dot line, the value of w is fixed, e.g., taking w = 0.8λ and the smaller value of σ x,max , i.e., taking σ x,max = 5.In contrast to the large air space case, reducing the σ x,max value neither increases the simulation time nor the memory requirement.Evidently, there are some advantages to the small σ x,max as compared to the large w in terms of mitigating late-time instability.

Conclusions
This work successfully performs a systematic study of the late-time instability characteristic of the elongated multilayer thin plate with the UPML and PBC boundary.The unstructured triangular finite element is employed to model the elongated domain, the 2D DGTD-UPML scheme is applied for the solution of Maxwell's equations.To mitigate the late-time instability caused by the ill-posed elongated thin target, there are two ways to be adopted: 1) increasing the width of the air space and 2) decreasing the UPML conductivity σ x,max .A large number of numerical trials have shown that both ways can decrease the gradual growth in the UPML region, ultimately, the late-time instability can be effective to be mitigated and reduced.Furthermore, numerical results illustrate using the small σ x,max is superior to using the large air space w.
About the Authors ...

Fig. 3 .
Fig. 3. RMS of the upper envelope of long time response of the electric field E ztest for different air space w.

Fig.
Fig. 4. RMS of the upper envelope of long time response of the electric field E ztest for different σ x,max .

Fig. 5 .
Fig. 5. RMS of the upper envelope of long time response of E ztest for large w and smallσ x,max .

Table 1
records the grid information for different width values.
Tab. 1. Grid information versus various widths of w.

4 .
RMS of the upper envelope of long time response of the electric field E ztest for different σ x,max .Late-time response of E ztest for different σ x,max .
Huiping LI received the master's degree in Communication and Information Systems from the Guangxi University, Guangxi, China, in 2007, and she is currently pursuing the Ph.D. degree in Electromagnetic Field and Microwave Technology from Nanjing University of Aeronautics and Astronautics, Nanjing, China, since 2014.She has been a lecturer with the School of Physics and Electronics, Henan University, Kaifeng, China, since 2008.Her current research interests include computational electromagnetism, wireless communication and circuit analysis.Yi WANG received the B.S. and Ph.D. degrees in Communication and Information System from Nanjing University of Aeronautics and Astronautics (NUAA), Nanjing, China in 2006 and 2012.After that, Yi Wang joined the College of Electronic and Information Engineering, NUAA, as an assistant Professor.His research interests include computational electromagnetics, especially the finite-difference time-domain (FDTD) method, the FDTD modeling of the entire Earth-ionosphere system, the earthquake electromagnetics and radome design.Qunsheng CAO received the Ph.D. degree in Electronic Engineering from The Hong Kong Polytechnic University, Hong Kong, in 2000.From 2000 to 2005, he was a Research Associate with the Department of Electrical Engineering, University of Illinois at Urbana-Champaign, and with the Army High Performance Computing Research Center (AHPCRC), University of Minnesota, respectively.In 2006, he joined the University of Aeronautics and Astronautics (NUAA), Nanjing, China, as a Professor of Electronic Engineering.He has authored or co-authored over 70 papers in refereed journals and conference proceedings.He has co-authored Multi-resolution Time Domain Scheme for Electromagnetic Engineering (Wiley, 2005).His current research interests are in computational electromagnetism and antennas designs, particularly in time domain numerical techniques for the study of microwave devices and scattering applications.