An Improved PML Implementation In The Transmission Line Method

A new framework for the development of the perfectly matched layer (PML) in the Transmission Line Modelling (TLM) method is demonstrated. The formulation is based on the stretched coordinate PML and involves the mapping of the TLM node to an analytically extended geometric domain. The implementation demonstrates the more stable second order discretization of the PML equations in TLM. The validity of the PML formulation is demonstrated on a canonical case of a metal waveguide and compared with an existing PML TLM scheme.


Introduction
The perfectly matched layer (PML), developed by Bérenger offers a powerful and effective way to model open boundary wave propagation problems [1]. Originally devised for the Finite-Difference Time Domain (FDTD) method, the PML is an artificial absorbing layer introduced at the boundaries of a computational domain and achieves a theoretical reflection free outward propagation by an impedance matching scheme. The success and simplicity of the PML has encouraged different revisions and extensions of Bérengers original theory and has resulted in various formulations beyond the initial FDTD technique [2]- [4]. To date, the performance of PMLs have been demonstrated in the transmission line model (TLM) method for the split-field [5,6] and stretched coordinate based (SCB) [7] PML interpretations, all of which have a first-order approach in their implementations. Whilst these have achieved some success, they also report such issues as temporal (exponential) growth in the solution domain and, a lower absorption performance when compared to their FDTD counterpart.
In this paper, we explore an alternative method to the implementation of PMLs based on the discretization of the second-order SCB PML equations. A motivation for this approach is reported by [8]- [9] as an improved stability performance of the discrete second-order PMLs compared to their first order complement. The proposed formulation follows a mapping of the TLM node to a complex stretched domain which is shown to transform the constituent LC transmission line components.
For simplicity, the theory is presented for the two-dimensional (2D) TLM node which can be straightforwardly extended to the three-dimensional (3D) case and as space precludes giving the details here, will be the subject of a future communication.
The results are demonstrated for the 3D canonical case of the metal waveguide where the formulation is validated by comparing the absorption performance and stability against an existing TLM PML scheme [6]. The framework presented can be easily extended to advanced TLM nodes thus encouraging a review of the PML in TLM.

Mapping the TLM node to a complex stretched domain
TLM is a well-established numerical technique developed based on an analogy between the wave equation and passive circuits. It offers an intuitive approach to the simulation of complex electromagnetic wave propagation problems by a translation to an equivalent transmission line network. Therefore, fundamental to the TLM method is establishing a field-circuit equivalence which allows the medium properties to be properly encapsulated into the LC components of the transmission lines [10].
A second-order 2D SCB PML wave equation for the transverse magnetic case is derived from [3] to obtain where and are magnetic permeability and electric permittivity, respectively and and represent the complex stretch factors given as where ≥ 0 is a conductivity function used to control the attenuation of the fields propagating in the i-direction, and 0 is the permittivity of free space. It can be shown [10] that there exists a circuit expression defined in a block of space analogous to Equation (1) given as where ∆ represents the dimension of the block of space, and represent the total inductance in x and y-directed transmission lines, respectively, and and represent the total capacitance and voltage in z-directed transmission line, respectively. The block of space modelled by Equation (3) is referred to as a shunt node which is well defined in the literature for real space i.e. = 1 [10]. By comparing the form of Equations (1) and (3) the following equivalences can be made between the field and circuit parameters in the complex stretched space: Applying the transmission line equations relating the characteristic admittance = √ / and propagation delay ∆ = ∆ √ , the equivalent circuits of the shunt node mapped to a complex stretched domain is derived as illustrated in Fig. 1a,b. In the 2D TLM node, the reflected voltage at each port n is given as [10] = − (4) where is the voltage incident on the line. By a simple circuit analysis of (Fig. 1b) the following expression is obtained Arriving at a time domain solution to Equation (5) is nontrivial and requires an appropriate handling of the complex and frequency dependent line admittances. This requires a series of transformations to be deployed for the implementation of digital filters that enable a mapping of frequency dependent material parameters to the time domain and these are [11]: Upon solving Equation (5)  where is a scaling factor dependent on , and is a nonphysical voltage introduced by the complex stretching.
Finally, as previously highlighted, the propagation delay along each i-direction in the mapped shunt node is a complex stretched variable given as According to Equation (7) the stretching is only applied to the imaginary time component, therefore the TLM criteria for impulse synchronism is maintained. The complex time component ∆ 0 is mathematically interpreted, through an inverse Fourier transform, as an attenuation factor − ∆ 0 for the propagating pulses along each i-directed line. This direct attenuation is applied to both pulses incident through the PML and those reflected within the PML region and is applied as a scaling factor in the connection process. For example, the voltage reflected from port 3 in k time step is incident on port 1 of the adjacent node in k+1 time step by the following rule

Numerical Experiments
The developments presented in section 2 have been extended to the symmetrical condensed node (SCN) to demonstrate the validity of the formulation in 3D. We compared the reflection performance and numerical stability of the mapped SCN PML against an existing SCN-PML scheme [6] for the canonical metal waveguide test case.

Reflection performance
The reflection performance was investigated for an empty rectangular waveguide WR90 with dimensions 22.86mm by 10.16mm by 40mm (∆ = 0.5mm) truncated by a 20layer PML with parabolic conductivity grading with an end conductivity = 5 / . The PML layer was backed by a perfect electrical conductor (PEC). The dominant 10 mode was excited by an amplitude modulated Gaussian pulse with a centre frequency of 10GHz and a bandwidth of 4GHz. The plane of excitation was located at one node from the PML-waveguide interface and the observation node was set at the other end of the waveguideone node away from the PML. An illustration of the simulated waveguide is shown in Fig.2. The reflection coefficient shown in Fig. 3 was computed using the reference solution of a longer waveguide terminated with PMLs. As observed, the mapped SCN demonstrates good agreement with the PML formulation proposed in [6].

Stability
We investigated the stability of our formulation by loading the waveguide simulated above with a capacitive iris. The iris was placed close to the PML (5 nodes away) to ensure high evanescent energy is present at the PML interface. The PML was terminated by a matched boundary condition and the end conductivity was set to = 28 / . Observation and excitation were made at the same node at the opposite end of the waveguide. The time domain plot of the component is presented in Fig. 4, where the mapped SCN PML is shown to exhibit no instability even with the iris placed at such close proximity to the PML. This demonstrates an improved performance over existing PML schemes simulated in a similar context [5]- [7]. The appearance of oscillations in the time domain arise due to the inability of the PML equations in the continuous form to properly absorb evanescent modes. This is well reported and not related to the numerical method employed.
To corroborate our claim of improved stability, other investigations have been carried out involving computational domains containing strong discontinuities which were simulated for time steps of five to six orders of magnitude and operated at the maximum time step. In these contexts, the stability of the mapped SCN PML also demonstrated an improved late time stability over existing TLM-PML schemes. A possible explanation for the stability improvement can be found in the approach employed by this formulation. Unlike the previous TLM-PML schemes, the mapping of the TLM node elaborates the transformation of the LC line components which are then accounted for by modifications in the TLM scatter-connect processes. In addition to this, a truly secondorder implementation is maintained by employing an equivalent network model [12], [13] in the derivation of its scatter equations. This is different to existing TLM-PML schemes which have been developed based on a direct firstorder discretization of the modified Maxwell's equation.

Conclusion
Perfectly matched layers are a critical enabler for large-scale electromagnetic simulations and it is imperative to optimise their performance. It is essential that PML implementations guarantee stability as well as high quality absorbing properties. A framework for the development of second-order SCB PMLs in TLM was proposed. The formulation is shown to modify the standard TLM scatter-connect processes whilst maintaining the fundamental TLM criteria of synchronism and connection. Conceptually the formulation is stable since each unit operation can be proved to be stable. This agrees with results obtained in which no instabilities have been found. The reflection performance is also shown to be in agreement with existing TLM-PML schemes. Based on its simple approach extensions will be made to advanced TLM nodes in the future.