# The CFS-PML for 2D Auxiliary Differential Equation FDTD Method Using Associated Hermite Orthogonal Functions

Authors: Feng Jiang; Xiao-Ping Miao; Feng Lu; Li-Yuan Su; Yao Ma
Journal: International Journal of Antennas and Propagation
 (2017)

---

## Abstract

The complex frequency shifted (CFS) perfectly matched layer (PML) is proposed for the two-dimensional auxiliary differential equation (ADE) finite-difference time-domain (FDTD) method combined with Associated Hermite (AH) orthogonal functions. According to the property of constitutive parameters of CFS-PML (CPML) absorbing boundary conditions (ABCs), the auxiliary differential variables are introduced. And one relationship between field components and auxiliary differential variables is derived. Substituting auxiliary differential variables into CPML ABCs, the other relationship between field components and auxiliary differential variables is derived. Then the matrix equations are obtained, which can be unified with Berenger’s PML (BPML) and free space. The electric field expansion coefficients can thus be obtained, respectively. In order to validate the efficiency of the proposed method, one example of wave propagation in two-dimensional free space is calculated using BPML, UPML, and CPML. Moreover, the absorbing effectiveness of the BPML, UPML, and CPML is discussed in a two-dimensional (2D) case, and the numerical simulations verify the accuracy and efficiency of the proposed method.


---

## Body

## 1. Introduction

In the conventional finite-difference time-domain (FDTD) method [1, 2], the time step is constrained by the Courant-Friedrichs-Lewy (CFL) stability condition. When fine structures such as thin material and slot are simulated, more computer memory and computation time are required. In order to eliminate the CFL stability condition, some unconditionally stable FDTD methods have been developed such as alternating-direction implicit (ADI) method [3], Crank-Nicolson method [4], locally one-dimensional method [5], and Weighted Laguerre Polynomials (WLP) FDTD method [6]. Recently, the Associated Hermite (AH) FDTD methods to model wave propagation have been introduced in [7, 8].The perfectly matched layer (PML) absorbing boundary conditions (ABCs) introduced by Berenger [9] have been widely used for truncating FDTD domains. Some PML variants have been proposed to improve the absorbing effectiveness for various electromagnetic waves, like modified PML (MPML), uniaxial anisotropic PML (UPML), generalized PML (GPML), and so forth. Recently, UPML-ABC for AH-FDTD method [10] is used in conductive medium. According to Kuzuoglu and Mittra in [11], the complex frequency shifted (CFS) PML [12, 13] is the most accurate in all the PML. The frequency-domain coordinate-stretching variable of CFS-PML (CPML) has three adjustable variables, while there are two adjustable variables for UPML-ABC. The reflection error of CPML can be greatly reduced from adjusting variables, and it has been implemented in the Cartesian coordinate, periodic structures, cylindrical coordinates, dispersive materials, WLP-FDTD, and so on.In this paper, a 2D CPML for auxiliary differential equation (ADE) FDTD method using AH orthogonal functions is proposed. It is shown that the constitutive parameters of CPML are complicated, and if the method in [7] is used directly here, the second derivative field components would be involved and the final matrix equations will also be complex. In order to get a simple and easy formula, the auxiliary differential variables are introduced. Based on the ADE technique [14, 15] and Galerkin temporal testing procedure, one relationship between field components and auxiliary differential variables is derived. According to auxiliary differential variables and Maxwell’s equations of CPML, the other relationship between field components and auxiliary differential variables is derived. Then the formulations of all orders of AH functions are obtained to calculate the magnetic field expansion coefficients. According to [7], the matrix equations of CPML, Berenger’s PML (BPML), and free space can be unified. At last, the electric field expansion coefficients can also be obtained, respectively. To validate the efficiency of the proposed method, a 2D case is calculated. And the efficiency of the proposed method is verified through the comparison with BPML and UPML ABCs.


## 2. Formulation

With free space, the frequency-domain Maxwell’s equations of CPML for 2DTEz model case are(1)jωε0Ex=1Ry∂Hz∂y,jωε0Ey=-1Rx∂Hz∂x,jωμ0Hzx=-1Rx∂Ey∂x,jωμ0Hzy=1Ry∂Ex∂y,where ε0 is the electric permittivity and μ0 is the magnetic permeability of free space; RΦ=κΦ+σΦ/(αΦ+jωε0), Φ=x,y. For CPML, αΦ and σΦ are assumed to be positive real and κΦ is real and ≥1.An orthonormal set of basis functionsφ0,φ1,φ2,φ3,… can be defined as(2)ϕnt~=2nn!π1/2-1/2e-t~2/2Hnt~,where Hnt~=-1net~2dn/dtne-t~2 are Hermite polynomials, t~=t-Tf/γ, Tf is a time-translating parameter, and γ is a time-scaling parameter. Choosing a finite Q order of basis functions [7] and proper parameters for Tf and γ, and then using these transformed basis functions, the causal field components, taking u for example, can be expanded as(3)ur,t=∑n=0Qunrϕ¯nt~.The time derivative of thenth order AH function is(4)ddtϕ¯nt~=-1γ12ϕ¯1t~,n=0,1γn2ϕ¯n-1t~-1γn+12ϕ¯n+1t~,n≥1.According to (4), the first derivative of field components ur,t with respect to t is [7](5)∂∂tur,t=1γ∑n=0Qun+1rn+12-un-1rn2ϕ¯nt~,where u-1r=0.Using ADE scheme [13, 14], four auxiliary differential variables are introduced:(6)ψEx=jωExαy+jωε0,(7)ψEy=jωEyαx+jωε0,(8)ψHzx=jωHzxαx+jωε0,(9)ψHzy=jωHzyαy+jωε0.With the transition relationship from frequency domain to time domain, (6) can be written as(10)αyψEx+ε0∂ψEx∂t=∂Ex∂t.Applying (3) and (5) to (10), the field functions in (10) can be expanded by(11)αyψExpr+ε01γ∑n=0QψExn+1rn+12-ψExn-1rn2ϕ¯nt~=1γ∑n=0QExn+1rn+12-Exn-1rn2ϕ¯nt~.Multiplying both sides by ϕ¯nt~ and integrating over (-∞,∞), we get(12)n+12ε0γψExq+1r+αyψExqr-n2ε0γψExq-1r=n+121γExq+1r-n21γExq-1r.Applying the same above procedure with (7)–(9), we have(13)n+12ε0γψEyq+1r+αxψEyqr-n2ε0γψEyq-1r=n+121γEyq+1r-n21γEyq-1r,n+12ε0γψHzxq+1r+αxψHzxqr-n2ε0γψHzxq-1r=n+121γHzxq+1r-n21γHzxq-1r,n+12ε0γψHzyq+1r+αyψHzyqr-n2ε0γψHzyq-1r=n+121γHzyq+1r-n21γHzyq-1r.Similar to [7, 8], we can rewrite (12)-(13) in a matrix form:(14)ξyψEx=ςEx,ξxψEy=ςEy,ξxψHzx=ςHzx,ξyψHzy=ςHzy,where(15)ξΦ=αΦ12ε0γ-12ε0γαΦ1+12ε0γ⋱⋱-Q-22ε0γαΦQ-2+12ε0γ-Q-12ε0γαΦΦ=x,yς=121γ-121γ1+121γ⋱⋱-Q-221γQ-121γ-Q-121γ.Equation (14) can be rewritten in a new matrix form:(16)ψEx=ξy-1ςEx,ψEy=ξx-1ςEy,ψHzx=ξx-1ςHzx,ψHzy=ξy-1ςHzy.Applying (6)–(9) to (1) and transferring from frequency domain to time domain, we have(17)ε0κy∂Ex∂t+σyψEx=∂Hz∂y,ε0κx∂Ey∂t+σxψEy=-∂Hz∂x,μ0κx∂Hzx∂t+σxHzx=-∂Ey∂x,μ0κy∂Hzy∂t+σyHzy=∂Ex∂y.Applying the same above procedure with (17), we have(18)κjαExqr+γσjβψExqr=γε0β∂Hzqr∂y,κiαEyqr+γσiβψEyqr=-γε0β∂Hzqr∂x,κi+0.5αHzxqr+γσi+0.5βψHzxqr=-γμ0β∂Eyqr∂x,κj+0.5αHzyqr+γσj+0.5βψHzyqr=γμ0β∂Exqr∂y,where(19)α=1-12⋱⋱1-Q-2Q-11-Q-1Q,β=2122⋱2Q-12Q.Using central difference scheme, substituting (16) into (18), grafting (18), and using Hz=Hzx+Hzy, we have (20)Exi+0.5,j=C¯yEi,jηj-1βHzi+0.5,j+0.5-Hzi+0.5,j-0.5,(21)Eyi,j+0.5=-C¯xEi,jηi-1βHzi+0.5,j+0.5-Hzi-0.5,j+0.5,(22)Hzi+0.5,j+0.5=C¯yHi,jηj+0.5-1βExi+0.5,j+1-Exi+0.5,j-C¯xHi,jηi+0.5-1βEyi+1,j+0.5-Eyi,j+0.5,where ηΦ=κΦα+γσΦβ[ξ]Φ-1[ς](Φ=x,y), C¯xEi,j=γ/ε0Δx¯i, C¯yEi,j=γ/ε0Δy¯j, C¯xHi,j=γ/μ0Δxi, C¯yHi,j=γ/μ0Δyj, and Δxi and Δyj are the lengths of the cell edge where the electric fields are located. Δx-i and Δy-j are the distances between the center nodes where the magnetic fields are located. For the above (i,j) subscripts for Cx and Cy in (20)–(22), (i,j) is not a real position but an array index of each field variable, as shown in Figure  1 of [16].In order to form the matrix equations only containing magnetic field, we apply (20)-(21) to (22): (23)blHzi-0.5,j+0.5+brHzi+1.5,j+0.5+aHzi+0.5,j+0.5+buHzi+0.5,j+1.5+bdHzi+0.5,j-0.5=0,where (24)bl=C¯xHi,jC¯xEi,jηi+0.5-1βηi-1β,br=C¯xHi,jC¯xEi+1,jηi+0.5-1βηi+1-1β,bu=C¯yHi,jC¯yEi,j+1ηj+0.5-1βηj+1-1β,bd=C¯yHi,jC¯yEi,jηj+0.5-1βηj-1β,a=-bu+bd+bl+br+I.I is unit matrix.From (23), only magnetic field vector variables remained. And it can also be found that each magnetic field variable is related to the four adjacent magnetic fields. If R=1 in (1), the final matrix equations of CPML will degrade to the matrix equations of free space in AH-FDTD method. If κ=1 and αΦ=0 in (1), the final matrix equations of CPML will degrade to the matrix equations of BPML in AH-FDTD method. So the matrix equations of CPML and BPML and free space can be unified. And this greatly facilitates the programming. Finally, a banded sparse coefficient matrix is obtained like [7], which contains all the points in free space and CPML ABC. And the paralleling-in-order solution scheme in [8] is applied to indirectly but efficiently calculate all of the expansion coefficients. If all of the expansion coefficients of the magnetic field are calculated, the expansion coefficients of the electric field can be obtained from (20)-(21). Finally, Exr,t, Eyr,t, and Hzr,t can be reconstructed by using (3).


## 3. Numerical Demonstration

In order to validate the efficiency of the presented method, wave scatter from a PEC medium is simulated. The BPML, UPML, and CPML ABCs are used to truncate the FDTD domains. The PML constitutive parameters are scaled usingmth order polynomial scaling [17]:(25)σρΦ=σΦ,maxρΦDΦmΦ=x,y,σopt=m+1150πΔρ,κρΦ=1+κΦ,max-1ρΦDΦmΦ=x,y,where ρ indicates the distance from the dispersive medium-PML interface into the PML, D is the depth of the PML, and m=4 is the order of polynomial.As illustrated in Figure1, the computational domain is truncated by 10 additional PML in x-direction and y-direction, respectively. And we choose b1 = b5 = a1 = a6 = 10 mm, a2 = 12 mm, a3 = 6 mm, a4 = 10 mm, a5 = 2 mm, b2 = 2 mm, b3 = 13 mm, and b4 = 15 mm (in Figure 1), which results in a 50×50 cells lattice with Δx=Δy=1 mm. For the PEC rectangle’s scatter, we choose c1 = 16 mm, c2 = 8 mm, and c3 = 10 mm.Figure 1
Computational domain.A sinusoidal-modulated Gaussian pulse is chosen as they-direction excitation source: (26)Jyt=e-t-tc2/td2sin⁡2πfct-tc,where td=1/2fc, tc=3td, and fc=0.8 GHz. The time duration of interest for the analyzed fields is chosen as Ts=50 ns and the bandwidth is limited up to the frequency Ws=3 GHz. We choose the order Q=300 and the range of γ = 1.3 × 10−9 from the condition above. The time-translating parameter Tf = 25 ns.In order to compare the performance of the proposed CPML with that of BPML, a reference solution was also simulated with no reflection coming from the boundary. And the same mesh is extended by 100 cells inx-direction and y-direction, leading to a 250×250 cells lattice. And the reference field is calculated using AH-FDTD method. The reflection error is calculated as follows:(27)RdB=20log10ExBt-ExTtmax⁡ExBt,where ExTt is the field computed in the test domain and ExBt is the reference field.In Figure2, the waveform of electric field at points P is graphed by reference field, BPML, UPML, and CPML. And it should be pointed out that parameters with different PML are different. Using (27), the reflection error is computed at the measurement point P. In Figure 3, it is noted that the maximum relative error is −43 dB, −58 dB, and −70 dB for the BPML, UPML, and CPML, respectively. And it is obvious that, with the CFS factor, the CPML is superior to BPML and UPML [10]. In Figures 2 and 3, the computation time is much higher than the time-support of source wave. And late-time reflections of CPML are much lower than those of UPML and BPML. The CFL stability condition of this model is Δt≤2.33 ps. The time step size for the proposed method is 23.3 ps. In Table 1, comparison of CPU resource for four methods is presented. In order to guarantee computation accuracy, the CPML [11, 18] is also adopted in WLP-FDTD and conventional FDTD method. The total memory storage for the proposed method is increased to 150.2 Mb, about 57.8 times of the conventional FDTD method and 3.8 times of WLP-FDTD method, while the total CPU time for the proposed method can be reduced to about 1.1% of the conventional FDTD method and 25.9% of WLP-FDTD method.Table 1
Comparison of CPU resource for various methods.
Methods
Δ
t (ps)
Memory (Mb)
CPU time (s)
Conventional FDTD
2.33
2.6
183
UPML of [10]
23.3
150.2
2.1
WLP-FDTD
23.3
39.5
8.1
CPML of proposed method
23.3
150.2
2.1Figure 2
Transient electric field ofx component at the observation point.Figure 3
Reflection error in the electric field intensity relative to the field’s maximum amplitude versus time for BPML (withD = 10 cells, σmax/σopt = 1, and m = 4), UPML (with D = 10 cells, κmax = 22, σmax/σopt = 1.2, and m = 4), and CPML (with D = 10 cells, α = 0.003, κmax = 72, σmax/σopt = 1.12, and m = 4) at point p.


## 4. Conclusion

Using ADE scheme, the CPML for 2D AH-FDTD method has been presented in this paper. This method is free from CFL stability condition for it has eliminated the time variable in calculation. The final matrix equations of CPML and BPML and free space can be unified. By solving the banded sparse coefficient matrix, magnetic field expansion coefficients of all orders of AH functions can be obtained. Then the electric field expansion coefficients can also be obtained, respectively. Numerical results show that this implementation is very effective in absorbing the electromagnetic waves, which means that the proposed method can save more computation time and computer memory.



---

*Source: 5817380-2017-05-18.xml*