# An Unconditionally Stable Cylindrical FDTD Method to Analyze the EM Ground Wave Propagation

Authors: Zheng-Yu Huang; Li-Hua Shi; Zheng Sun; Ying-Hui Zhou
Journal: International Journal of Antennas and Propagation
 (2018)

---

## Abstract

We extended the unconditionally stable associated Hermite (AH) FDTD method to the cylindrical coordinate system for analyzing the electromagnetic (EM) ground wave propagation. With AH domain differentiate operator technology and paralleling-in-order solution scheme, the cylindrical time-domain Maxwell equations with complex frequency shifted perfectly matched layer (CFS-PML) are transformed to AH domain generating a five-point banded equation, which also unified the formula from the central axis of cylinder. The numerical results represent the accuracy and efficiency of the proposed method when comparing for the conventional FDTD method.


---

## Body

## 1. Introduction

In the finite-difference time-domain (FDTD) method [1], the time step size is limited by the Courant-Friedrich-Levy (CFL) stability condition [2]. So for the problems with fine structures, it should choose a small time step and then a large number of marching-on-in-time steps for simulation, which is a time-consuming process. To eliminate the stability condition for fine structure, many unconditionally stable (US) methods were proposed, such as alternating direction implicit (ADI) method [3] and orthogonal function-based methods [4–9]. The orthogonal function-based methods can be classified into two main categories: the marching-on-in-order scheme and the paralleling-in-order scheme. The former is based on weighted Laguerre polynomials (WLP), and the latter is based on associated Hermite (AH) functions. For the latter, all of the works are based on in Cartesian grid previously.In many applications, such as EM ground wave propagation from a lightning channel, a cylindrical coordinate is preferred. Here, we extend the AH FDTD method to cylindrical coordinates system to analyze this case. The time-domain Maxwell equations in cylindrical coordinates are transformed to AH domain using a AH differential operator technology [6]. Then a five-point banded equation is derived for calculation of expanding coefficients. And a paralleling-in-order scheme [6] is used for indirectly but efficiently solution. For the CFS-PML is wildly used to effectively attenuate evanescent waves [10, 11], we apply it to AH FDTD method and make a unified formula including a treatment of central axis of cylinder when following the Ampere law. One should note that although CFS-PML is used in [6, 9], it is implemented in Cartesian coordinate system. In this paper, we firstly extended it to cylindrical coordinate system, which might give AH FDTD method much more convenient to analyze the scattering problems for the cylindrical shape-based objective. Finally, a numerical example is used to validate the accuracy and efficiency of the proposed method.


## 2. Formulation

With the lossy dielectric media, the 2D cylinder Maxwell equations with CFS-PML can be written as [11]
(1)εr∂Er∂t+σerEr=1sez∂Hz∂z−Jr,(2)εz∂Ez∂t+σezEz=−1ser∂rHφr∂r−Jz,(3)μφ∂Hφ∂t+σmφHφ=1smz∂Er∂z−1smr∂Ez∂r−Mφ.where μφ, σmz, εξ, and σeξξ=r,z are the permeability, the magnetic conductivity, the permittivity, and the electric conductivity of the media, respectively. seξ and smξ are the frequency-domain coordinate-stretching variables, which are defined as seξ=κeξ+σpeξ/ηeξ+jωε0 and smξ=κmξ+σpmξ/ηmξ+jωμ0, where σpeξ, σpmξ, ηeξ, ηmξ, κeξ, and κmξ are the respective PML parameters [11]. Similar to [6], a set of AH basis functions ϕnt,Tf,l is chosen to expand all of the field components (1), (2), and (3). l is the time-scaling parameter, and Tf is the time-translating parameter. Then, a temporal Galerkin testing procedure is used to eliminate the time variables. By applying central difference scheme, we can transform (1), (2), and (3) to AH domain as
(4)aeri,kEri,k=aSezi,k−1Hφi,k−Hφi,k−1Δz¯k−Jri,k,(5)aezi,kEzi,k=−aSeri,k−11+1/2i−1Hφi,k−1−1/2i−1Hφi−1,kΔr¯i−Jzi,k,(6)amφi,jHφi,k=aSmri,k−1Ezi+1,k−Ezi,kΔri−aSmzi,k−1Eri,k+1−Eri,kΔzk−Mφi,k,where
(7)a=22l1−12−2⋱⋱Q−1−Q−1Q×Q,(8)aeξi,k=εξi,ka+σeξi,kI,(9)amφi,k=μφi,ka+σmφi,kI,(10)aSeξi,k=κeξi,kI+σpeξi,kηeξi,kI+aε0−1,(11)aSmξi,k=κmξi,kI+σpmξi,kηmξi,kI+aμ0−1.The matrixa above is AH domain differential operator [6], which can be considered as the similar time-domain operator d/dt or frequency-domain operator jω and respective identity matrix I is for 1. By assembling (4), (5), and (6) and eliminating the electric field components, we can derive a five-point banded matrix equation with Hφ as
(12)ali,kHφi−1,k+ari+1,kHφi+1,k+ami,kHφi,k+adi,kHφi,k−1+aui,k+1Hφi,k+1=bi,k,where
(13)aui,k+1=aSmzi,k−1aSezi,k+1−1aeri,k+1−1Δz¯k+1/Δzk,adi,k=aSmyi,k−1aSeyi,k−1aeri,k−1Δz¯k/Δzk,ali,k=1−1/2i−1aSmri,k−1aSeri,k−1aezi,k−1Δri/Δri,ari+1,k=1+1/2iaSmri,k−1aSeri+1,k−1aezi+1,k−1Δri+1/Δri,ami,k=−2i−12i+1ari+1,k+2i−12i−3ali,k+aui,k+1+adi,k+amφi,k,bi,k=aSmzi,k−1amri,k+1−1Jri,k+1−amri,j−1Jri,k/Δz−aSmri,k−1amzi+1,k−1Jzi+1,k−amzi,j−1Jzi,k/Δr+Mφi,k.When for the boundary condition of central axis, it needs to be treated separately owing to singularities. According to Ampere’s law [12],
(14)∮lHφ⋅dl=∬S∂∂tεEz+σeEz⋅ds.Then, the differential form for (14) can be obtained as follows:
(15)2πrHφr,z,t=πr2ε∂Ezr,z,t∂t+σeEr,z,t.Converting it to the AH domain by using the AH domain differential operatora and performing discretization for i=1
(16)Ez1,j=2raezi,k−1Hφ1,j.By applying (16) to (5), the modified magnetic field (12) in central axis can be derived as follows:
(17)aezi+1,k−1Δr21+12iHφi+1,k+aeξi,k+1−1Δz2Hφi,k+1+aeri,k−1Δz2Hφi,k−1−aezi+1,k−1Δr21−12i+aeri,k+1−1Δz2+aeri,k−1Δz2+4εri,kαΔr2+amφi,kHφi,k=1Δraezi+1,k−1Jzi+1,k−aezi,j−1Jzi,k−1Δzaeri,k+1−1Jri,k+1−aeri,k−1Jri,k.After the central axis formula to modify (12) is updated, we can apply the lower-upper (LU) decomposition procedure to calculate this equation. Also, we can perform paralleling-in-order solution scheme [6] for (12). If every magnetic field variable Hφ∣i,j in the entire computational space is calculated, the results of electronic field variables Er∣i,j and Ez∣i,j can be gotten from (4) and (5). In addition, only the lossy medium is considered above, while for the case of dispersive medium, the formulation derivation is almost the same as before, except for a slight modification of (8) with εξ being replaced by its AH domain form from [7].


## 3. Numerical Verification

To validate the performance of the proposed AH FDTD method, a numerical example with an EM ground wave propagation problem is presented, as shown in Figure1. The computational domain is discretized into 50 × 50 cells with Δr=Δz=0.1m, including an axial symmetric boundary condition for the left boundary and the other three boundaries terminated by 10 PML layers. The background material is set to be earth with medium parameters: εξ=12, σeξ=0.05, μφ=1, and σmz=0. The PML parameters are scaled following the expressions in [11].
(18)κξ=1+κξmax−1ξ−ξ0mdm,σξ=σξmaxξ−ξ0mdm,σopt=m+1150πΔξ,where ξ0 represents the interface between FDTD and PML grids, d is the thickness of the PML, and m=4 is the order of the polynomial. A Gaussian pulse source Jz modulated by sinusoidal signal is located at the grid (0, 27).
(19)Jzt=exp−t−TcTd2sin2πfct−Tc,where Td=0.5fc, Tc=3Td, and fc=0.05MHz. Two observation points p1 (20, 23) and p2 (38, 38) are chosen to record the results.Figure 1
The configuration of the computational domain.To verify the performance of the proposed method, the result from conventional FDTD method is used as reference. The time step for conventional FDTD method is chosen asΔt=0.21ns to meet the CFL condition, and the time duration is chosen as 127.4 μs, which leads to the total time-marching number of 600,000. While for AH FDTD method, Δt can be chosen as 21 ns to evaluate expanding coefficients of source. AH parameters are selected as the number of AH functions Q=48, the time-translating parameter Tf=63.7μs, and time-scaling parameter l=6.8×10−6, respectively. PML parameters are selected to meet the best performances in this example as ηξ=0.005, κmax=7, and σmax=1.2σopt.The magnetic fieldsHφ at p1 and p2 are shown in Figure 2, where we can see the good agreement between the proposed method and the conventional FDTD method.Figure 2
Comparison of the magnetic fields calculated by proposed method and conventional FDTD at two observation points.The comparison of computational resources between the AH FDTD method and the conventional FDTD method is listed in Table1. Although the proposed method has much more memory consumption, it gives much higher efficiency in CPU time and there is about 95% reduction of computing time in contrast with conventional FDTD method.Table 1
Comparison of the computational resources.
△t (ns)
Memory (Mb)
CPU times (s)
Conventional FDTD
0.21
1.78
587.4
Proposed method
21
24.9
31.2


## 4. Conclusion

The AH FDTD method is developed to analyze the ground wave propagation problem under the cylindrical coordinates in this paper. With the CFS-PML absorbing boundary, a unified banded equation, including the central axis formula, is derived in AH domain by using the paralleling-in-order solution scheme and AH differential operator technology. The numerical example with a cylindrical ground wave propagation case is given, which also makes a good verification of the performance for the proposed method. Therefore, for its good performance of numerical dispersion, the AH FDTD method has a potential to analyze the propagation of waves on long distances such as ground penetrating radar (GPR) applications, where segmentation technology should be explored to deal with a probable long-time response (e.g., signals as biexponentials) in the future.



---

*Source: 3046838-2018-05-08.xml*