# An Efficient Laguerre-Based FDTD Iterative Algorithm in 3D Cylindrical Coordinate System

Authors: Da-Wei Zhu; Hai-Lin Chen; Zi-Yu Pan; Bo-Ao Xu; Bin Chen
Journal: International Journal of Antennas and Propagation
 (2019)

---

## Abstract

Here an efficient Laguerre-based finite-difference time-domain iterative algorithm is proposed. Different from the previously developed iterative procedure used in the efficient FDTD algorithm, a new perturbation term combined with the Gauss–Seidel iterative procedure is introduced to form the new Laguerre-based FDTD algorithm in the 3D cylindrical coordinate system. Such a treatment scheme can reduce the splitting error to a low level and obtain a good convergence; in other words, it can improve the efficiency and accuracy than other algorithms. To verify the performance of the proposed algorithm, two scattering numerical examples are given. The computation results show that the proposed algorithm can be better than the ADI-FDTD algorithm in terms of efficiency and accuracy. Meanwhile, the proposed algorithm is extremely useful for the problems with fine structures in the 3D cylindrical coordinate system.


---

## Body

## 1. Introduction

In recent years, the unconditionally stable Laguerre-based finite-difference time-domain (FDTD) algorithm has been applied to simulate transient electromagnetic problems in the Cartesian coordinate. By transforming the time-domain problem to the Laguerre domain using the temporal Galerkin’s testing procedure, the transient solution is independent of time discretization. Thus, Laguerre-based FDTD formulation has the advantage of less numerical error when a larger time step is used.The main drawback of the conventional Laguerre-based FDTD algorithm is that it requires solving a very large sparse matrix. To overcome this problem, a factorization-splitting scheme [1] was used to resolve the huge sparse matrix into six tridiagonal matrixes, and then the chasing method was used to solve them [2]. This provided excellent computational accuracy and can be efficiently parallel-processed on a computing cluster. To further reduce the splitting error, the modified perturbation term was introduced in [3, 4]; meanwhile, the Gauss–Seidel iterative procedure was applied, which made the new Laguerre-based FDTD scheme more efficient and accurate. Obviously, the good feasibility of the efficient Laguerre scheme had been conformed in the rectangular coordinate system.However, in many applications, we have to deal with 3D cylindrical structures such as in optical fiber communication, integrated optics, and defense industry. Moreover, the geometry of interest may consist of fine structures. If we adopt the conventional FDTD method to discretize the cylindrical structure with the Cartesian grid, a significant staircasing error appears. In fact, due to the existence of the1/ρ term in the 3D cylindrical coordinate system, the decomposition of the fields along the ϕ direction will lead the splitting error. How to reduce this kind of error in the iterative FDTD calculation has become a difficult problem. In addition, how to conduct ϕ direction transformation of the fields on the z-axis is also a research difficult in the 3D cylindrical coordinate system. Perhaps this is why, since 2000, the FDTD method has done little applications in the 3D cylindrical coordinate system. There is not even an article on the Laguerre technology.Therefore, in order to expand the research field of the FDTD method, we propose an efficient Laguerre-based finite-difference time-domain iterative algorithm in the 3D cylindrical coordinate system. Firstly, according to the ideology of the weighted Laguerre polynomial (WLP) FDTD scheme, the WLP-FDTD equations of the proposed algorithm in the 3D cylindrical coordinate system are deduced by introducing the new perturbation term and the nonphysical intermediate variables. Secondly, the field components on thez-axis are amended with special treatment. Finally, to verify the proposed algorithm, two scattering numerical examples are given. Numerical results show that the proposed algorithm can be better than the ADI-FDTD algorithm [5] in terms of efficiency and accuracy.


## 2. Iterative Theory and Formulations

### 2.1. WLP-FDTD Equations in 3D Cylindrical Coordinate System

Introducing the WLP technology [6] into Maxwell’s equations of the 3D cylindrical coordinate [7], one can obtain the equations of the conventional WLP-FDTD in the 3D cylindrical coordinate system, directly:(1)eρq+2∑k=0q−1eρk=2sε1ρ∂ϕhzq−∂zhϕq,(2)eϕq+2∑k=0q−1eϕk=2sε∂zhρq−∂ρhzq,(3)ezq+2∑k=0q−1ezk=2sε1ρ∂ρhϕq∂ρ−1ρ∂ϕhρq,(4)hρq+2∑k=0q−1hρk=2sμ∂zeϕq−1ρ∂ϕezq,(5)hϕq+2∑k=0q−1hϕk=2sμ∂ρezq−∂zeρq,(6)hzq+2∑k=0q−1hzk=2sμ1ρ∂ϕeρq−1ρ∂ρeϕq∂ρ,where s>0 is a time-scale factor, q is the order of weighted Laguerre polynomials, and ε and μ are the electric and magnetic permeability, respectively.


### 2.2. Iterative Equations for the Off-Axis Fields

Equations (1)–(6) can be written in the following matrix form:(7)Wq−Vq−1=A+BWq,where Wq and Vq−1 are the general form of q‐order WLP fields and the summation of q−1 order WLP fields, respectively, and A and B are the coefficient matrices.Obviously, solving equation (7) must involve large sparse matrix problems, and it is impossible to solve it. Therefore, the efficient algorithm must be adopted to avoid the direct solution of the large sparse matrix.Referring to the brief form of the 3-D FDTD equation in [4], when the excitation sources JM are added, the equations of the efficient Laguerre-based FDTD in the 3D cylindrical coordinate can be rewritten as(8)I−AI−BWq=−ABWq−1+Vq−1+JMq,(9)I−AI−BWp+1q=ABWpq+Vq−1+JMq,where p is the pth iteration and equations (8) and (9) are the initial and iterative computation, respectively. Clearly, (9) have the same form with (8). So, one can combine them as follows:(10a)Wpq=−Wq−1,(10b)I−AI−BWp+1q=ABWpq+Vq−1+JMq,where(11)Wq=EϑqHϑqT,Vq−1=Eϑq−1Hϑq−1T,Eϑq=eρqeϕqezqT,Hϑq=hρqhϕqhzqT,Eϑq−1=−2∑k=0q−1eρk−2∑k=0q−1eϕk−2∑k=0q−1ezkT,Hϑq−1=−2∑k=0q−1hρk−2∑k=0q−1hϕk−2∑k=0q−1hzkT,JMq=JMEqJMHqT,JMEq=−aJeρq−aJeϕq−aJezqT,JMHq=bMhρqbMhϕqbMhzqT.In addition, theA and B in (10b) should be(12)A=0A12A210,B=0B12B210,A12=0−a∂z000−a∂ρ−a∂ϕρ00,A21=00−b∂ϕρ−b∂z000−b1ρ∂∂ρρ0,B12=00a∂ϕρa∂z000a1ρ∂∂ρρ0,B21=0b∂z000b∂ρb∂ϕρ00,where a=2/sε and b=2/sμ. ∂ρ, ∂ϕ, and ∂z are the first-order central difference operators along ρ, ϕ, and z axes, respectively.In fact, equation (10a) can be implemented by simple programming, and we only need to implement the design of the iterative calculation. Therefore, we actually only need to introduce the perturbation term ABWp+1q−Wpq in [4] into equation (7) to obtain (10b), and the advantage of this is that it greatly simplifies the programming process and improves the computational efficiency.In accordance with the design theory of the proposed algorithm in this paper, we introduce the nonphysical intermediate variablesW∗q=Eϑ∗qHϑ∗qT=eρ∗qeϕ∗qez∗qhρ∗qhϕ∗qhz∗qT, and it can divide (10b) into(13a)I−AWp+1∗q=BWpq+Vq−1+JMq,(13b)I−BWp+1q=Wp+1∗q−BWpq.Expanding (13a) and (13b), one can obtain(14a)Eϑ,p+1∗q−A12Hϑ,p+1∗q=B12Hϑ,pq+Eϑq−1+JMEq,(14b)Hϑ,p+1∗q−A21Eϑ,p+1∗q=B21Eϑ,pq+Hϑq−1+JMHq,(15a)Eϑ,p+1q−B12Hϑ,p+1q=Eϑ,p+1∗q−B12Hϑ,pq,(15b)Hϑ,p+1q−B21Eϑ,p+1q=Hϑ,p+1∗q−B21Eϑ,pq.Substituting (14b) into (14a), one can obtain(16)I−A12A21Eϑ,p+1∗q=A12B21Eϑ,pq+B12Hϑ,pq+A12Hϑq−1+Eϑq−1+A12JMHq+JMEq.Expanding (16), one can obtain three implicit equations eρ,p+1∗q, eϕ,p+1∗q, and ez,p+1∗q. Here, the Gauss–Seidel iterative procedure should be used, and the equations after the Gauss–Seidel procedure are(17a)1−ab∂z2eρ,p+1∗q=−ab∂z∂ρez,pq+a∂ϕρhz,pq+V1q−1+JM1q,(17b)1−ab∂ρ1ρ∂∂ρρeϕ,p+1∗q=−ab∂ρ∂ϕρeρ,p+1∗q+a∂zhρ,pq+V2q−1+JM2q,(17c)1−ab∂ϕρ2ez,p+1∗q=−ab∂ϕρ∂zeϕ,p+1∗q+a1ρ∂∂ρhϕ,pq+V3q−1+JM3q,(17a‐1)V1q−1=2a∂z∑k=0q−1hϕk−2∑k=0q−1eρk,(17b‐1)V2q−1=2a∂ρ∑k=0q−1hzk−2∑k=0q−1eϕk,(17c‐1)V3q−1=2a∂ϕρ∑k=0q−1hρk−2∑k=0q−1ezk,(17a‐2)JM1q=−ab∂zMhϕq−aJeρq,(17b‐2)JM2q=−ab∂ρMhzq−aJeϕq,(17c‐2)JM3q=−ab∂ϕρMhρq−aJezq.Obviously, the necessary nonphysical variableseρ,p+1∗q, eϕ,p+1∗q, and ez,p+1∗q for solving (15a) have been obtained through solving (17a)–(17c).Next is how to solve the fields of the nonphysical variablesHϑ,p+1∗q for solving (15b). Here, it is not necessary because substituting (14b) into (15b), equation (15b) can be rewritten as(18)Hϑ,p+1q=B21Eϑ,p+1q+A21Eϑ,p+1∗q+Hϑq−1+JMHq.Clearly, the nonphysical variablesHϑ,p+1∗q are eliminated already. Meaning that, the nonphysical variables Hϑ,p+1∗q do not have to be solved in the whole iteration.Substituting (18) into (15a), one can obtain a new equation which is similar with (16):(19)I−B12B21Eϑ,p+1q=Eϑ,p+1∗q+B12A21Eϑ,p+1∗q−B12Hϑ,pq+B12Hϑq−1+B12JMHq.Expanding (18) and (19), one can obtain the computational equations of the electric and magnetic fields of the proposed algorithm:(20)1−ab∂ϕρ2eρ,p+1q=eρ,p+1∗q−ab∂ϕρ1ρ∂∂ρρeϕ,p+1∗q−a∂ϕρhz,pq+V4q−1+JM4q,(21)1−ab∂z2eϕ,p+1q=eϕ,p+1∗q−ab∂z∂ϕρez,p+1∗q−a∂zhρ,pq+V5q−1+JM5q,(22)1−ab1ρ∂∂ρρ∂ρez,p+1q=ez,p+1∗q−ab1ρ∂∂ρρ∂zeρ,p+1∗q−a1ρ∂∂ρρhϕ,pq+V6q−1+JM6q,where(20‐1)V4q−1=−2a∂ϕρ∑k=0q−1hzk,(20‐2)JM4q=ab∂ϕρMhzq,(21‐1)V5q−1=−2a∂z∑k=0q−1hρk,(21‐2)JM5q=ab∂zMhρq,(22‐1)V6q−1=−2a1ρ∂∂ρρ∑k=0q−1hϕk,(22‐2)JM6q=ab1ρ∂∂ρρMhϕq,(23)hρ,p+1q=b∂zeϕ,p+1q−b∂ϕρez,p+1∗q−2∑k=0q−1hρk+bMhρq,(24)hϕ,p+1q=b∂ρez,p+1q−b∂zeρ,p+1∗q−2∑k=0q−1hϕk+bMhϕq,(25)hz,p+1q=b∂ϕρeρ,p+1q−b1ρ∂∂ρρeϕ,p+1∗q−2∑k=0q−1hzk+bMhzq.At this point, the implementation scheme of the proposed algorithm has been elaborated, and the overall execution process is described in Figure1.Figure 1
Overall execution process of the proposed algorithm.


## 2.1. WLP-FDTD Equations in 3D Cylindrical Coordinate System

Introducing the WLP technology [6] into Maxwell’s equations of the 3D cylindrical coordinate [7], one can obtain the equations of the conventional WLP-FDTD in the 3D cylindrical coordinate system, directly:(1)eρq+2∑k=0q−1eρk=2sε1ρ∂ϕhzq−∂zhϕq,(2)eϕq+2∑k=0q−1eϕk=2sε∂zhρq−∂ρhzq,(3)ezq+2∑k=0q−1ezk=2sε1ρ∂ρhϕq∂ρ−1ρ∂ϕhρq,(4)hρq+2∑k=0q−1hρk=2sμ∂zeϕq−1ρ∂ϕezq,(5)hϕq+2∑k=0q−1hϕk=2sμ∂ρezq−∂zeρq,(6)hzq+2∑k=0q−1hzk=2sμ1ρ∂ϕeρq−1ρ∂ρeϕq∂ρ,where s>0 is a time-scale factor, q is the order of weighted Laguerre polynomials, and ε and μ are the electric and magnetic permeability, respectively.


## 2.2. Iterative Equations for the Off-Axis Fields

Equations (1)–(6) can be written in the following matrix form:(7)Wq−Vq−1=A+BWq,where Wq and Vq−1 are the general form of q‐order WLP fields and the summation of q−1 order WLP fields, respectively, and A and B are the coefficient matrices.Obviously, solving equation (7) must involve large sparse matrix problems, and it is impossible to solve it. Therefore, the efficient algorithm must be adopted to avoid the direct solution of the large sparse matrix.Referring to the brief form of the 3-D FDTD equation in [4], when the excitation sources JM are added, the equations of the efficient Laguerre-based FDTD in the 3D cylindrical coordinate can be rewritten as(8)I−AI−BWq=−ABWq−1+Vq−1+JMq,(9)I−AI−BWp+1q=ABWpq+Vq−1+JMq,where p is the pth iteration and equations (8) and (9) are the initial and iterative computation, respectively. Clearly, (9) have the same form with (8). So, one can combine them as follows:(10a)Wpq=−Wq−1,(10b)I−AI−BWp+1q=ABWpq+Vq−1+JMq,where(11)Wq=EϑqHϑqT,Vq−1=Eϑq−1Hϑq−1T,Eϑq=eρqeϕqezqT,Hϑq=hρqhϕqhzqT,Eϑq−1=−2∑k=0q−1eρk−2∑k=0q−1eϕk−2∑k=0q−1ezkT,Hϑq−1=−2∑k=0q−1hρk−2∑k=0q−1hϕk−2∑k=0q−1hzkT,JMq=JMEqJMHqT,JMEq=−aJeρq−aJeϕq−aJezqT,JMHq=bMhρqbMhϕqbMhzqT.In addition, theA and B in (10b) should be(12)A=0A12A210,B=0B12B210,A12=0−a∂z000−a∂ρ−a∂ϕρ00,A21=00−b∂ϕρ−b∂z000−b1ρ∂∂ρρ0,B12=00a∂ϕρa∂z000a1ρ∂∂ρρ0,B21=0b∂z000b∂ρb∂ϕρ00,where a=2/sε and b=2/sμ. ∂ρ, ∂ϕ, and ∂z are the first-order central difference operators along ρ, ϕ, and z axes, respectively.In fact, equation (10a) can be implemented by simple programming, and we only need to implement the design of the iterative calculation. Therefore, we actually only need to introduce the perturbation term ABWp+1q−Wpq in [4] into equation (7) to obtain (10b), and the advantage of this is that it greatly simplifies the programming process and improves the computational efficiency.In accordance with the design theory of the proposed algorithm in this paper, we introduce the nonphysical intermediate variablesW∗q=Eϑ∗qHϑ∗qT=eρ∗qeϕ∗qez∗qhρ∗qhϕ∗qhz∗qT, and it can divide (10b) into(13a)I−AWp+1∗q=BWpq+Vq−1+JMq,(13b)I−BWp+1q=Wp+1∗q−BWpq.Expanding (13a) and (13b), one can obtain(14a)Eϑ,p+1∗q−A12Hϑ,p+1∗q=B12Hϑ,pq+Eϑq−1+JMEq,(14b)Hϑ,p+1∗q−A21Eϑ,p+1∗q=B21Eϑ,pq+Hϑq−1+JMHq,(15a)Eϑ,p+1q−B12Hϑ,p+1q=Eϑ,p+1∗q−B12Hϑ,pq,(15b)Hϑ,p+1q−B21Eϑ,p+1q=Hϑ,p+1∗q−B21Eϑ,pq.Substituting (14b) into (14a), one can obtain(16)I−A12A21Eϑ,p+1∗q=A12B21Eϑ,pq+B12Hϑ,pq+A12Hϑq−1+Eϑq−1+A12JMHq+JMEq.Expanding (16), one can obtain three implicit equations eρ,p+1∗q, eϕ,p+1∗q, and ez,p+1∗q. Here, the Gauss–Seidel iterative procedure should be used, and the equations after the Gauss–Seidel procedure are(17a)1−ab∂z2eρ,p+1∗q=−ab∂z∂ρez,pq+a∂ϕρhz,pq+V1q−1+JM1q,(17b)1−ab∂ρ1ρ∂∂ρρeϕ,p+1∗q=−ab∂ρ∂ϕρeρ,p+1∗q+a∂zhρ,pq+V2q−1+JM2q,(17c)1−ab∂ϕρ2ez,p+1∗q=−ab∂ϕρ∂zeϕ,p+1∗q+a1ρ∂∂ρhϕ,pq+V3q−1+JM3q,(17a‐1)V1q−1=2a∂z∑k=0q−1hϕk−2∑k=0q−1eρk,(17b‐1)V2q−1=2a∂ρ∑k=0q−1hzk−2∑k=0q−1eϕk,(17c‐1)V3q−1=2a∂ϕρ∑k=0q−1hρk−2∑k=0q−1ezk,(17a‐2)JM1q=−ab∂zMhϕq−aJeρq,(17b‐2)JM2q=−ab∂ρMhzq−aJeϕq,(17c‐2)JM3q=−ab∂ϕρMhρq−aJezq.Obviously, the necessary nonphysical variableseρ,p+1∗q, eϕ,p+1∗q, and ez,p+1∗q for solving (15a) have been obtained through solving (17a)–(17c).Next is how to solve the fields of the nonphysical variablesHϑ,p+1∗q for solving (15b). Here, it is not necessary because substituting (14b) into (15b), equation (15b) can be rewritten as(18)Hϑ,p+1q=B21Eϑ,p+1q+A21Eϑ,p+1∗q+Hϑq−1+JMHq.Clearly, the nonphysical variablesHϑ,p+1∗q are eliminated already. Meaning that, the nonphysical variables Hϑ,p+1∗q do not have to be solved in the whole iteration.Substituting (18) into (15a), one can obtain a new equation which is similar with (16):(19)I−B12B21Eϑ,p+1q=Eϑ,p+1∗q+B12A21Eϑ,p+1∗q−B12Hϑ,pq+B12Hϑq−1+B12JMHq.Expanding (18) and (19), one can obtain the computational equations of the electric and magnetic fields of the proposed algorithm:(20)1−ab∂ϕρ2eρ,p+1q=eρ,p+1∗q−ab∂ϕρ1ρ∂∂ρρeϕ,p+1∗q−a∂ϕρhz,pq+V4q−1+JM4q,(21)1−ab∂z2eϕ,p+1q=eϕ,p+1∗q−ab∂z∂ϕρez,p+1∗q−a∂zhρ,pq+V5q−1+JM5q,(22)1−ab1ρ∂∂ρρ∂ρez,p+1q=ez,p+1∗q−ab1ρ∂∂ρρ∂zeρ,p+1∗q−a1ρ∂∂ρρhϕ,pq+V6q−1+JM6q,where(20‐1)V4q−1=−2a∂ϕρ∑k=0q−1hzk,(20‐2)JM4q=ab∂ϕρMhzq,(21‐1)V5q−1=−2a∂z∑k=0q−1hρk,(21‐2)JM5q=ab∂zMhρq,(22‐1)V6q−1=−2a1ρ∂∂ρρ∑k=0q−1hϕk,(22‐2)JM6q=ab1ρ∂∂ρρMhϕq,(23)hρ,p+1q=b∂zeϕ,p+1q−b∂ϕρez,p+1∗q−2∑k=0q−1hρk+bMhρq,(24)hϕ,p+1q=b∂ρez,p+1q−b∂zeρ,p+1∗q−2∑k=0q−1hϕk+bMhϕq,(25)hz,p+1q=b∂ϕρeρ,p+1q−b1ρ∂∂ρρeϕ,p+1∗q−2∑k=0q−1hzk+bMhzq.At this point, the implementation scheme of the proposed algorithm has been elaborated, and the overall execution process is described in Figure1.Figure 1
Overall execution process of the proposed algorithm.


## 3. Special Treatment for the On-Axis Fields

### 3.1. Theory and Analysis of the Special Treatment for the On-Axis Field Components

Because of the strangeness of the cylindrical coordinate system, drastic change of the field nearby the axis will lead the splitting error to be great [8]. Therefore, on-axis fields need to be specially treated in the 3D cylindrical coordinate system.The positions of the electric/magnetic fields of the 3D cylindrical coordinate are shown in Figure2, deducting the fields nearby the axis, and the ρoϕ section is a partial circular ring. Seeing the black shadow sector in Figure 3, clearly, the ρoϕ section can be added by many sectors. That is to say, all the electric/magnetic fields on the axis exhibit mathematical strangeness. As shown in Figure 3, when ρ=0, there are Ez and Hz field components on the axis; therefore, we should treat Ez0,j,k+1/2 and Hz1/2,j+1/2,k specially.Figure 2
Position of field in cylindrical coordinates.Figure 3
Top view of the field nearby the axis.According Stokes’s theorem, we start from the following integral form of Maxwell’s equations in the time domain:(26)ε∬s∂E⇀∂t⋅ds⇀=∮lH⇀⋅dl⇀,(27)μ∬s∂H⇀∂t⋅ds⇀=−∮lE⇀⋅dl⇀,where l⇀ is a closed curve surrounding the ρ=0 axis and s⇀ is the surface bounded by the curve l⇀. Equations (26) and (27) are actually the integral form of the Ampere loop theorem.


### 3.2. Differentia of the Special Treatment for the On-Axis Field Components in 3D Cylindrical Coordinate System

The FDTD difference equation ofEz0,j,k+1/2 can be derived from the formula (26), its integral area is the dotted line circle, and the radius is Δρ/2, as shown in Figure 3:(28)πΔρ22ε∂Ez0,j,k+1/2∂t=∑j1=1jmHϕ12,j1,k+12Δρ2Δϕ,noting j1≠j.Obviously, the form of the right side of equation (28) is different from equations (17a)–(17c-2) in [9]. Supposing Hϕ1/2,j1,k+1/2 in equation (28) is calculated in a body of revolution structure and each Hϕ1/2,j1,k+1/2 is equal for all the different j1 values, equation (28) can be modified as(29)ε∂Ez0,j,k+1/2∂t=2jmΔϕπΔρHϕq12,j1,k+12=4ΔρHϕq12,j1,k+12,where jmΔϕ=2π.In this way, the form of the right side of equation (29) is similar with (17a)–(17c-2) in [9]. By comparing equations (28) and (29), we can explain the different characteristics of the special treatment scheme on the axis in the 3D cylindrical coordinate structure, which is clearly distinguished from the scheme in the body of the revolution structure.


### 3.3. Special Treatment ofEzq0,j,k+1/2 and Hzq1/2,j+1/2,k

Introducing the WLP technology [6] into equation (28), one can obtain(30)Ezq0,j,k+12+2∑k=0q−1Ezk0,j,k+12=2aΔϕπΔρ∑j1=1jmHϕq12,j1,k+12.In order to solve equation (30), we present the first-order FDTD central difference expansion of equation (24) at ρ=0:(31)hϕ,p+1q12,j,k+12=bΔρez,p+1q1,j,k+12−ez,p+1q0,j,k+12−bΔzeρ,p+1∗q12,j,k+1−eρ,p+1∗q12,j,k−2∑k=0q−1hϕk12,j,k+12+bMhϕq12,j,k+12.Substituting (31) into (30), one can obtain the special treatment equation of Ez0,j,k+1/2 in the 3D cylindrical coordinate system:(32)1+2jmabΔϕπΔρ2ez,p+1q0,j,k+12−2abΔϕπΔρ2∑j1=1jmez,p+1q1,j1,k+12=−2abΔϕπΔρΔz∑j1=1jmeρ,p+1∗q12,j1,k+1−eρ,p+1∗q12,j1,k−4aΔϕπΔρ∑j1=1jm∑k=0q−1hϕk12,j1,k+12−2∑k=0q−1ezk0,j,k+12+2abΔϕπΔρ∑j1=1jmMhϕq12,j1,k+12.In a similar way, the FDTD difference equation ofHz1/2,j+1/2,k can be derived from formula (27), its integral area is the dotted line circle, and the radius is Δρ, as shown in Figure 3. Its method of obtaining can refer to equations (30) to (32).Obviously, whenρ=0, the expansion form of equation (25) is consistent with Hz1/2,j+1/2,k; therefore, the special treatment of Hz1/2,j+1/2,k is not necessary.


## 3.1. Theory and Analysis of the Special Treatment for the On-Axis Field Components

Because of the strangeness of the cylindrical coordinate system, drastic change of the field nearby the axis will lead the splitting error to be great [8]. Therefore, on-axis fields need to be specially treated in the 3D cylindrical coordinate system.The positions of the electric/magnetic fields of the 3D cylindrical coordinate are shown in Figure2, deducting the fields nearby the axis, and the ρoϕ section is a partial circular ring. Seeing the black shadow sector in Figure 3, clearly, the ρoϕ section can be added by many sectors. That is to say, all the electric/magnetic fields on the axis exhibit mathematical strangeness. As shown in Figure 3, when ρ=0, there are Ez and Hz field components on the axis; therefore, we should treat Ez0,j,k+1/2 and Hz1/2,j+1/2,k specially.Figure 2
Position of field in cylindrical coordinates.Figure 3
Top view of the field nearby the axis.According Stokes’s theorem, we start from the following integral form of Maxwell’s equations in the time domain:(26)ε∬s∂E⇀∂t⋅ds⇀=∮lH⇀⋅dl⇀,(27)μ∬s∂H⇀∂t⋅ds⇀=−∮lE⇀⋅dl⇀,where l⇀ is a closed curve surrounding the ρ=0 axis and s⇀ is the surface bounded by the curve l⇀. Equations (26) and (27) are actually the integral form of the Ampere loop theorem.


## 3.2. Differentia of the Special Treatment for the On-Axis Field Components in 3D Cylindrical Coordinate System

The FDTD difference equation ofEz0,j,k+1/2 can be derived from the formula (26), its integral area is the dotted line circle, and the radius is Δρ/2, as shown in Figure 3:(28)πΔρ22ε∂Ez0,j,k+1/2∂t=∑j1=1jmHϕ12,j1,k+12Δρ2Δϕ,noting j1≠j.Obviously, the form of the right side of equation (28) is different from equations (17a)–(17c-2) in [9]. Supposing Hϕ1/2,j1,k+1/2 in equation (28) is calculated in a body of revolution structure and each Hϕ1/2,j1,k+1/2 is equal for all the different j1 values, equation (28) can be modified as(29)ε∂Ez0,j,k+1/2∂t=2jmΔϕπΔρHϕq12,j1,k+12=4ΔρHϕq12,j1,k+12,where jmΔϕ=2π.In this way, the form of the right side of equation (29) is similar with (17a)–(17c-2) in [9]. By comparing equations (28) and (29), we can explain the different characteristics of the special treatment scheme on the axis in the 3D cylindrical coordinate structure, which is clearly distinguished from the scheme in the body of the revolution structure.


## 3.3. Special Treatment ofEzq0,j,k+1/2 and Hzq1/2,j+1/2,k

Introducing the WLP technology [6] into equation (28), one can obtain(30)Ezq0,j,k+12+2∑k=0q−1Ezk0,j,k+12=2aΔϕπΔρ∑j1=1jmHϕq12,j1,k+12.In order to solve equation (30), we present the first-order FDTD central difference expansion of equation (24) at ρ=0:(31)hϕ,p+1q12,j,k+12=bΔρez,p+1q1,j,k+12−ez,p+1q0,j,k+12−bΔzeρ,p+1∗q12,j,k+1−eρ,p+1∗q12,j,k−2∑k=0q−1hϕk12,j,k+12+bMhϕq12,j,k+12.Substituting (31) into (30), one can obtain the special treatment equation of Ez0,j,k+1/2 in the 3D cylindrical coordinate system:(32)1+2jmabΔϕπΔρ2ez,p+1q0,j,k+12−2abΔϕπΔρ2∑j1=1jmez,p+1q1,j1,k+12=−2abΔϕπΔρΔz∑j1=1jmeρ,p+1∗q12,j1,k+1−eρ,p+1∗q12,j1,k−4aΔϕπΔρ∑j1=1jm∑k=0q−1hϕk12,j1,k+12−2∑k=0q−1ezk0,j,k+12+2abΔϕπΔρ∑j1=1jmMhϕq12,j1,k+12.In a similar way, the FDTD difference equation ofHz1/2,j+1/2,k can be derived from formula (27), its integral area is the dotted line circle, and the radius is Δρ, as shown in Figure 3. Its method of obtaining can refer to equations (30) to (32).Obviously, whenρ=0, the expansion form of equation (25) is consistent with Hz1/2,j+1/2,k; therefore, the special treatment of Hz1/2,j+1/2,k is not necessary.


## 4. Numerical Results

In order to verify the performance of the proposed algorithm, two scattering numerical examples are given.First example: an irregular scatter with two grooves is calculated, and the angle of the groove isα=π/4. The radius of the groove is 19 cm, and the radius and height of the cylinder are 20 cm and 40 cm, respectively. The oblique incident plane wave is added through Huygen’s surface at ρ=22.5cm and z=±22.5cm, and parallel to the ϕ=0 plane, the incident angle is θi=45°, as shown in Figure 4. The incident electric field is(33)Eit=E0cos2πf0texp−4πt−t0τ2,where E0=103V/m, f0=10GHz, τ=1.0/f0, and t0=1.25/f0. And s=7×1011 and q=180 are chosen for the computation of the field values. The computational domain is meshed using a uniform grid with Δρ=Δz=1mm, leading to the total mesh of 40ρ×40ϕ×80z, and the first-order Mur absorbing boundary condition is used to truncate the boundary [10].Figure 4
Illustration of the irregular cylinder with two grooves under the oblique incident wave.The observation points of the simulated electric fields are set atP129,22,0, P229,2,19, and P325,25,−15, and the iteration of the whole computational domain is denoted Nt. By using the conventional FDTD algorithm, the ADI-FDTD algorithm [5], and the proposed algorithm, the time-domain field waveforms are shown in Figures 5–7. The numerical results show an excellent agreement between the conventional FDTD and the proposed algorithms.Figure 5
Computational results ofEϕ at point P1 of the conventional FDTD algorithm, the ADI-FDTD algorithm [5], and the proposed algorithm.Figure 6
Computational results ofEϕ at point P2 of the conventional FDTD algorithm, the ADI-FDTD algorithm [5], and the proposed algorithm.Figure 7
Computational results ofEϕ at point P3 of the conventional FDTD algorithm, the ADI-FDTD algorithm [5], and the proposed algorithm.In addition, to demonstrate the higher accuracy of the proposed algorithm, here two error formulations are provided to describe the accuracy, which are defined as(34‐1)Eϕ,error=Eϕt−Eϕ,FDTDt,(34‐2)Eϕ,error=Eϕt−Eϕ,FDTDtEϕtmax%,where Eϕt are the simulation results of the electric fields of the ADI-FDTD algorithm and the proposed algorithm and Eϕ,FDTDt is the simulation result of the conventional FDTD algorithm.Figures8–10 afford the error about E of the proposed algorithm (with Nt = 3). Clearly, the error of the proposed algorithm is much smaller than the ADI-FDTD algorithm (with cfln=4, cfln=5, and cfln=6, cfln=Δt/Δtcfl). The CPU time of the ADI-FDTD algorithm and the proposed algorithm is shown in Table 1.Figure 8
Errors obtained byEϕ at point P1 of the ADI-FDTD algorithm [5] and the proposed algorithm.Figure 9
Errors obtained byEϕ at point P2 of the ADI-FDTD algorithm [5] and the proposed algorithm.Figure 10
Errors obtained byEϕ at point P3 of the ADI-FDTD algorithm [5] and the proposed algorithm.Table 1
Simulation results of the two algorithms in the first example.
ADI-FDTD algorithm [5]
Proposed algorithm
cfln = 4
cfln = 5
cfln = 6
Nt = 3
CPU time (s)
337.27
270.14
224.26
218.95Observing Figures8–10 and Table 1, when the proposed algorithm Nt = 3 and the ADI-FDTD algorithm cfln=6, the CPU time of both algorithms is almost equal; nevertheless, the error of the proposed algorithm can be reduced almost 64% than the ADI-FDTD algorithm at least. Even though, when the ADI-FDTD cfln=4, the comparing results show that the error of the proposed algorithm can be reduced almost 44% at least and the CPU time can be saved almost 35% than the ADI-FDTD algorithm. Obviously, the proposed algorithm can provide better improvement in terms of efficiency and accuracy.Second example: another irregular scatter with one wedge-shaped bulge is shown in Figure11: the angle of the bulge is π/4, the radius of the bulge is 22 cm, and the radius and height of the cylinder are 20 cm and 40 cm, respectively. The other conditions are same with the first example.Figure 11
Illustration of the irregular cylinder with one wedge-shaped bulge under the oblique incident wave.The observation points of simulated electric fields are set atP129,10,15, P229,20,0, and P315,38,−29, respectively. The waveforms of the proposed algorithm are shown in Figures 12–14, which also form an excellent agreement with the conventional FDTD algorithm.Figure 12
Computational results ofEϕ at point P1 of the conventional FDTD algorithm, the ADI-FDTD algorithm [5], and the proposed algorithm.Figure 13
Computational results ofEϕ at point P2 of the conventional FDTD algorithm, the ADI-FDTD algorithm [5], and the proposed algorithm.Figure 14
Computational results ofEϕ at point P3 of the conventional FDTD algorithm, the ADI-FDTD algorithm [5], and the proposed algorithm.To demonstrate the accuracy of the proposed algorithm again, formulations (34-1) and (34-2) are used again here. Figures 15–17 afford the error about E of the proposed algorithm (with Nt = 3) and the ADI-FDTD algorithm (with cfln=4, cfln=5, and cfln=6). The CPU time of the both algorithms is shown in Table 2.Figure 15
Errors obtained byEϕ at point P1 of the ADI-FDTD algorithm [5] and the proposed algorithm.Figure 16
Errors obtained byEϕ at point P2 of the ADI-FDTD algorithm [5] and the proposed algorithm.Figure 17
Errors obtained byEϕ at point P3 of the ADI-FDTD algorithm [5] and the proposed algorithm.Table 2
Simulation results of the two algorithms in the second example.
ADI-FDTD algorithm [5]
Proposed algorithm
cfln = 4
cfln = 5
cfln = 6
Nt = 3
CPU time (s)
335.48
271.93
226.06
222.99Observing Figures15–17 and Table 2, when the proposed algorithm Nt = 3 and the ADI-FDTD algorithm (with cfln=4, cfln=5 and cfln=6), the comparison results are similar with the first example, namely, the efficiency and the accuracy of the proposed algorithm are still better than the ADI-FDTD algorithm. Obviously, the second example numerical results verify that the proposed algorithm can provide better improvement than the ADI-FDTD algorithm [5] in terms of efficiency and accuracy again.


## 5. Conclusion

In this paper, an efficient Laguerre-based finite-difference time-domain iterative algorithm in the 3D cylindrical coordinate system is proposed. By adopting a new perturbation termABWp+1q−Wpq and using the Gauss–Seidel iterative procedure method in the iterative algorithm of this paper, this kind of treatment can get better convergence; meanwhile, it also can ensure the enhancing of the efficiency and accuracy. It can be proved by computing the two scattering numerical examples that the proposed algorithm can provide better efficiency and accuracy in high-frequency range than the ADI-FDTD algorithm.



---

*Source: 9542976-2019-09-25.xml*