A conformal symplectic multi-resolution time-domain method with ADE-PML

Abstract

In this paper, the conformal finite-difference time-domain (CFDTD) technique and the concept of effective dielectric constant are applied to symplectic multi-resolution time-domain (SMRTD) algorithm based on Daubechies scaling functions, resulting in a conformal SMRTD (C-SMRTD) method. Moreover, a novel implementation of perfect matched layer with auxiliary differential equation (ADE-PML) is presented for the SMRTD method. Numerical results of the radar cross section (RCS) of a dielectric sphere verify the effectiveness of the C-SMRTD method, which shows good agreement with the results derived by analytical methods. It is also shown that the accuracy of the C-SMRTD method has been significantly improved compared with the MRTD method with step approximation, and the proposed method has an advantage over CFDTD at saving the CPU time and memory. Furthermore, a vault-top tunnel model is established with the proposed method, and the characteristic of the field cross-section distribution of the electromagnetic pulse (EMP) propagation in the tunnel model is also presented and discussed.

Keywords

Conformal technique symplectic multi-resolution time-domain coarse meshes RCS PML

1. Introduction

With highly-linear dispersion performance and accuracy, the SMRTD method implies that a low sampling rate in space can still provide for a relatively small phase error in the numerical simulation of a wave propagation problem, so it becomes possible that larger targets can be simulated without sacrificing accuracy [1]. However, the high numerical accuracy can only be maintained in the homogeneous medium. In order to accurately simulate the irregular shaped target, it is necessary to divide the mesh into very fine, which greatly reduces the advantages of the SMRTD method.

For the electromagnetic modeling of non-uniform boundary, subgridding algorithm [2], conformal gird method [3–10] have been used in the traditional FDTD method, and are increasingly used in grid generation, stability control, adaptive segmentation and so on. However, due to the limitation of dispersion and stability of FDTD algorithm, the computing resources needed for the simulation of large-scale targets are amazing. In order to solve this problem, this paper introduces the conformal technique in FDTD to the SMRTD algorithm, resulting in the C-SMRTD method. According to the local conformal technique, only the electric field on the deformed grid near the target surface needs modified. The magnetic field on the deformed grid are still calculated as normal MRTD. The modified electric field on the deformed grid is usually realized by equivalent electromagnetic parameters.

In this paper, the ADE-PML for the SMRTD is also described in detail. Numerical results show that the technique is more efficient at numerical reflection than that of the widely used anisotropic perfect matched layer (APML) [11–13]. Moreover, a side-wall vault-top tunnel model is established with the proposed conformal SMRTD method, and the characteristic of field cross-section distribution of electromagnetic pulse (EMP) propagation excited by TE₁₀ mode is studied. It should be notes that the compact support wavelet-Daubechies with two vanishing moments (D2) [14–17] is employed in this work.

2. C-SMRTD method

2.1. Multi domain decomposition of electric field component

Without loss of generality, here we take the electric field component E _x as an example. According to the symmetry relation of the connection coefficient a (l) [16], the semi discrete iterative formula of the E _x in SMRTD can be written as follows in the case of uniform dispersion (Δx = Δy = Δz = Δ). $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}\frac{\partial E_{i+1/2,j,k}^{{\phi}x}(t)}{\partial t}+{\sigma}E_{i+1/2,j,k}^{{\phi}x}(t)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \quad =a(0)\frac{1}{{\rm\Delta}}(H_{i+1/2,j+1/2,k}^{{\phi}z}(t)-H_{i+1/2,j-1/2,k}^{{\phi}z}(t)-H_{i+1/2,j,k+1/2}^{{\phi}y}(t)+H_{i+1/2,j,k-1/2}^{{\phi}y}(t))\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \qquad +∼3a(1)\frac{1}{3{\rm\Delta}}(H_{i+1/2,j+3/2,k}^{{\phi}z}(t)-H_{i+1/2,j-3/2,k}^{{\phi}z}(t)-H_{i+1/2,j,k+3/2}^{{\phi}y}(t)+H_{i+1/2,j,k-3/2}^{{\phi}y}(t))\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \qquad +\cdots +(2L_{s}-1)a(2L_{s}-1)\frac{1}{(2L_{s}-1){\rm\Delta}}\left(\begin{array}{@{}l@{}}H_{i+1/2,j+L_{s}-1/2,k}^{{\phi}z}(t)-H_{i+1/2,j-L_{s}+1/2,k}^{{\phi}z}(t)\\[5.0pt] -H_{i+1/2,j,k+L_{s}-1/2}^{{\phi}y}(t)+H_{i+1/2,j,k-L_{s}+1/2}^{{\phi}y}(t)\end{array}\right)\end{eqnarray}$ (1) Where L _s is the effective support size of Daubechies scale function. For Daubechies scale function, in order to ensure that the numerical results converge to the correct solution when the space grid size tends to 0, the connection coefficient a must satisfy [17]: $\begin{eqnarray}\displaystyle \mathop{\sum }_{l=0}^{L_{s}-1}a(l)\cdot (2l+1)=1. & & \displaystyle\end{eqnarray}$ (2) Therefore, Eq. (1) can be divided into the following subequations: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle a(0){\varepsilon}\frac{\partial E_{i+1/2,j,k}^{{\phi}x}(t)}{\partial t}+a(0){\sigma}E_{i+1/2,j,k}^{{\phi}x}(t)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \quad =a(0)\frac{1}{{\rm\Delta}}(H_{i+1/2,j+1/2,k}^{{\phi}z}(t)-H_{i+1/2,j-1/2,k}^{{\phi}z}(t)-H_{i+1/2,j,k+1/2}^{{\phi}y}(t)+H_{i+1/2,j,k-1/2}^{{\phi}y}(t))\end{eqnarray}$ (3) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle 3a(1){\varepsilon}\frac{\partial E_{i+1/2,j,k}^{{\phi}x}(t)}{\partial t}+3a(1){\sigma}E_{i+1/2,j,k}^{{\phi}x}(t)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \quad =3a(1)\frac{1}{3{\rm\Delta}}(H_{i+1/2,j+3/2,k}^{{\phi}z}(t)-H_{i+1/2,j-3/2,k}^{{\phi}z}(t)-H_{i+1/2,j,k+3/2}^{{\phi}y}(t)+H_{i+1/2,j,k-3/2}^{{\phi}y}(t))\end{eqnarray}$ (4) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \vdots \nonumber\\ \displaystyle \text{} & \text{} & \displaystyle (2L_{s}-1)a(L_{s}-1){\varepsilon}\frac{\partial E_{i+1/2,j,k}^{{\phi}x}(t)}{\partial t}+(2L_{s}-1)a(L_{s}-1){\sigma}E_{i+1/2,j,k}^{{\phi}x}(t)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \quad =(2L_{s}-1)a(L_{s}-1)\frac{1}{(2L_{s}-1){\rm\Delta}}\left(\begin{array}{@{}l@{}}H_{i+1/2,j+L_{s}-1/2,k}^{{\phi}z}(t)-H_{i+1/2,j-L_{s}+1/2,k}^{{\phi}z}(t)\\[5.0pt] -H_{i+1/2,j,k+L_{s}-1/2}^{{\phi}y}(t)+H_{i+1/2,j,k-L_{s}+1/2}^{{\phi}y}(t)\end{array}\right).\end{eqnarray}$ (5)

Equation (3)–Eq. (5) show that the semi discrete iterative formula of the SMRTD is actually the linear combination of traditional FDTD algorithms with mesh sizes of Δ, 3Δ… and (2L _s −1)Δ, respectively. For the SMRTD with support size L _s = 3, it is actually the linear combination of the traditional FDTD with three mesh sizes of Δ, 3Δ and 5Δ. Figure 1 shows the schematic diagram of multi-region decomposition of electric field component, and the three rectangles from the inside to the outside are the integration paths of traditional FDTD with grid sizes Δ, 3Δ and 5Δ.

Fig. 1.

Multi region decomposition of electric field component.

2.2. Solution of equivalent electromagnetic parameters

In CFDTD algorithm, volume weighted average method [9], area weighted average method [10] and length weighted average method [7] are usually used to solve the equivalent electromagnetic parameters. Compared with the volume and area weighted methods, length weighted method has advantages in both the realization and the calculation accuracy. As shown in Fig. 2, the dielectric parameter is modified on the corresponding edge by the length weighted average method, while the parameter of the edge that does not intersect with the medium surface remains unchanged. And ${\varepsilon}_{x}^{eq},{\varepsilon}_{y}^{eq}$ are shown in Eq. (6) and Eq. (7) ( ${\sigma}_{x}^{eq}$ and ${\sigma}_{y}^{eq}$ have similar expressions) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}_{x}^{eq}(i,j,k)=[{\rm\Delta}x_{1}(i,j,k)\times {\varepsilon}_{1}+({\rm\Delta}x(i,j,k)-{\rm\Delta}x_{1}(i,j,k))\times {\varepsilon}_{2}]/{\rm\Delta}x(i,j,k)\end{eqnarray}$ (6) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}_{y}^{eq}(i,j,k)=[{\rm\Delta}y_{1}(i,j,k)\times {\varepsilon}_{1}+({\rm\Delta}y(i,j,k)-{\rm\Delta}y_{1}(i,j,k))\times {\varepsilon}_{2}]/{\rm\Delta}y(i,j,k).\end{eqnarray}$ (7)

For symplectic SMRTD algorithm, the scope of deformable mesh should be extended. For the grid near the dielectric target surface, if at least one of the FDTD decomposition integral surfaces of the electric field component multi region decomposition of SMRTD on the grid is partially filled by the dielectric target, then it is called the extended deformation grid, as shown in Fig. 3(a)–Fig. 3(c). Figure 3(a) shows that only one decomposition integral surface is filled with media, Fig. 3(b) shows that part of decomposition integral surfaces are filled with media, and Fig. 3(c) shows that all decomposition integral surfaces are filled with media.

Fig. 2.

Schematic diagram of length weighted average method.

Fig. 3.

Extended deformation grid. (a) Only one decomposition integral surface is filled with medium. (b) Part of the decomposition integral surfaces are filled with medium. (c) All decomposition integral surfaces are filled with medium.

Taking the E _x(i +1∕2, j, k) shown in Fig. 3 as an example, the shaded part in the figure is the dielectric target. The SMRTD iterative formula based on D2 scale function can be decomposed into three traditional FDTD iterative formulas with grid sizes of Δ, 3Δ and 5Δ. The electromagnetic parameters ϵ, σ in the SMRTD iterative formula of E _x(i +1∕2, j, k) are replaced by the equivalent permittivity ϵ_eq and the equivalent conductivity σ_eq respectively. According to the symmetry relationship of the connection coefficient a (l), we can get: $\begin{eqnarray}\displaystyle {\varepsilon}_{eq}\frac{\partial E_{i+1/2,j,k}^{{\phi}x}(t)}{\partial t}+{\sigma}_{eq}E_{i+1/2,j,k}^{{\phi}x}(t)=\mathop{\sum }_{l=0}^{L_{s}-1}a(l)\frac{1}{{\rm\Delta}}\left(\begin{array}{@{}l@{}}H_{i+1/2,j+(l+1)-1/2,k}^{{\phi}z}(t)-H_{i+1/2,j-(l+1)+1/2,k}^{{\phi}z}(t)\\[5.0pt] -H_{i+1/2,j,k+(l+1)-1/2}^{{\phi}y}(t)+H_{i+1/2,j,k-(l+1)+1/2}^{{\phi}y}(t)\end{array}\right). & & \displaystyle\end{eqnarray}$ (8)

Replace the electromagnetic parameters ϵ, σ in Eqs (3)–(5) with the equivalent electromagnetic parameters, we can get $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle a(0){\varepsilon}_{eq}(0)\frac{\partial E_{i+1/2,j,k}^{{\phi}x}(t)}{\partial t}+a(0){\sigma}_{eq}(0)E_{i+1/2,j,k}^{{\phi}x}(t)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \quad =a(0)\frac{1}{{\rm\Delta}}(H_{i+1/2,j+1/2,k}^{{\phi}z}(t)-H_{i+1/2,j-1/2,k}^{{\phi}z}(t)-H_{i+1/2,j,k+1/2}^{{\phi}y}(t)+H_{i+1/2,j,k-1/2}^{{\phi}y}(t))\end{eqnarray}$ (9) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle 3a(1){\varepsilon}_{eq}(1)\frac{\partial E_{i+1/2,j,k}^{{\phi}x}(t)}{\partial t}+3a(1){\sigma}_{eq}(1)E_{i+1/2,j,k}^{{\phi}x}(t)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \quad =3a(1)\frac{1}{3{\rm\Delta}}(H_{i+1/2,j+3/2,k}^{{\phi}z}(t)-H_{i+1/2,j-3/2,k}^{{\phi}z}(t)-H_{i+1/2,j,k+3/2}^{{\phi}y}(t)+H_{i+1/2,j,k-3/2}^{{\phi}y}(t))\end{eqnarray}$ (10) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \vdots \nonumber\\ \displaystyle \text{} & \text{} & \displaystyle (2L_{s}-1)a(L_{s}-1){\varepsilon}_{eq}(L_{s}-1)\frac{\partial E_{i+1/2,j,k}^{{\phi}x}(t)}{\partial t}+(2L_{s}-1)a(L_{s}-1){\sigma}_{eq}(L_{s}-1)E_{i+1/2,j,k}^{{\phi}x}(t)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \quad =(2L_{s}-1)a(L_{s}-1)\frac{1}{(2L_{s}-1){\rm\Delta}}\left(\begin{array}{@{}l@{}}H_{i+1/2,j+L_{s}-1/2,k}^{{\phi}z}(t)-H_{i+1/2,j-L_{s}+1/2,k}^{{\phi}z}(t)\\[5.0pt] -H_{i+1/2,j,k+L_{s}-1/2}^{{\phi}y}(t)+H_{i+1/2,j,k-L_{s}+1/2}^{{\phi}y}(t)\end{array}\right)\end{eqnarray}$ (11) Where ϵ_eq(l), σ_eq(l) (l = 0,1,2, …, L _s −1) are the equivalent permittivity and equivalent conductivity with mesh sizes of Δ, 3Δ, … and (2L _s −1)Δ, respectively. Add Eqs (9)–(11), and the following formula can be obtained. $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \mathop{\sum }_{l=0}^{L_{s}-1}a(l)(2l+1){\varepsilon}_{eq}(l)\frac{\partial E_{i+1/2,j,k}^{{\phi}x}(t)}{\partial t}+\mathop{\sum }_{l=0}^{L_{s}-1}a(l)(2l+1){\sigma}_{eq}(l)E_{i+1/2,j,k}^{{\phi}x}(t)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \quad =\mathop{\sum }_{l=0}^{L_{s}-1}a(l)\frac{1}{{\rm\Delta}}\left(\begin{array}{@{}l@{}}H_{i+1/2,j+(l+1)-1/2,k}^{{\phi}z}(t)-H_{i+1/2,j-(l+1)+1/2,k}^{{\phi}z}(t)\\[5.0pt] -H_{i+1/2,j,k+(l+1)-1/2}^{{\phi}y}(t)+H_{i+1/2,j,k-(l+1)+1/2}^{{\phi}y}(t)\end{array}\right).\end{eqnarray}$ (12) Then we can get the equivalent electromagnetic parameters on the deformed mesh $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}_{eq}=\mathop{\sum }_{l=0}^{L_{s}-1}a(l)(2l+1){\varepsilon}_{eq}(l)\end{eqnarray}$ (13) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\sigma}_{eq}=\mathop{\sum }_{l=0}^{L_{s}-1}a(l)(2l+1){\sigma}_{eq}(l).\end{eqnarray}$ (14) By substituting the equivalent electromagnetic parameters into the iterative formula of E _x and discretizing them with symplectic algorithm in time, we can get the conformal SMRTD iterative formula of E _x $\begin{eqnarray}\displaystyle E_{i+1/2,j,k}^{{\phi}x,n+v/m} & = & \displaystyle \frac{2{\varepsilon}_{eq}-{\sigma}_{eq}d_{v}{\rm\Delta}t}{2{\varepsilon}_{eq}+{\sigma}_{eq}d_{v}{\rm\Delta}t}E_{i+1/2,j,k}^{{\phi}x,n+(v-1)/m}\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle +\frac{2d_{v}{\rm\Delta}t}{2{\varepsilon}_{eq}+{\sigma}_{eq}d_{v}{\rm\Delta}t}\mathop{\sum }_{l}a(l)\frac{1}{{\rm\Delta}}[H_{i+1/2,j+l+1/2,k}^{{\phi}z,n+v/m}-H_{i+1/2,j,k+l+1/2}^{{\phi}y,n+v/m}].\end{eqnarray}$ (15) The iterative formula of other electric field components can be obtained by the same method.

3. Implementation of the ADE-PML

3.1. Formulation

In this section, the ADE-PML for SMRTD based on Daubechies scaling functions is discussed. For the sake of generality example, a lossy medium is assumed here. In the PML layer, the formulation is posed in the stretched coordinate space $\begin{eqnarray}\displaystyle j{\omega}{\varepsilon}E_{x}+{\sigma}E_{x}=\frac{1}{s_{y}}\frac{\partial }{\partial y}H_{z}-\frac{1}{s_{z}}\frac{\partial }{\partial z}H_{y} & & \displaystyle\end{eqnarray}$ (16) where s _i are the stretched-coordinate metric, which are proposed to be $\begin{eqnarray}\displaystyle s_{i}={\kappa}_{i}+\frac{{\sigma}_{i}}{{\alpha}_{i}+j{\omega}{\varepsilon}_{0}},\quad i=x,y,z & & \displaystyle\end{eqnarray}$ (17) and then we can get $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}{\displaystyle \frac{1}{s_{i}}}={\displaystyle \frac{1}{{\kappa}_{i}}}+{\displaystyle \frac{{\sigma}_{i}^{\prime }}{{\alpha}_{i}^{\prime }+j{\omega}{\varepsilon}_{i}^{\prime }}}\\[10.0pt] ({\sigma}_{i}^{\prime }=-{\sigma}_{i},∼{\varepsilon}_{i}^{\prime }={\varepsilon}_{0}{{\kappa}_{i}}^{2},∼{\alpha}_{i}^{\prime }={{\kappa}_{i}}^{2}{\alpha}_{i}+{\kappa}_{i}{\sigma}_{i})\end{array}\right. & & \displaystyle\end{eqnarray}$ (18) where σ_i, 𝜅_i and 𝛼_i are nonnegative reals, moreover, 𝜅_i ≥ 1.

Inserting (18) into (16), we obtain $\begin{eqnarray}\displaystyle j{\omega}{\varepsilon}E_{x}+{\sigma}E_{x} & = & \displaystyle \frac{1}{{\kappa}_{y}}\frac{\partial }{\partial y}H_{z}-\frac{1}{{\kappa}_{z}}\frac{\partial }{\partial z}H_{y}+\frac{{\sigma}_{y}^{\prime }}{{\alpha}_{y}^{\prime }+j{\omega}{\varepsilon}_{y}^{\prime }}\frac{\partial }{\partial y}H_{z}-\frac{{\sigma}_{z}^{\prime }}{{\alpha}_{z}^{\prime }+j{\omega}{\varepsilon}_{z}^{\prime }}\frac{\partial }{\partial z}H_{y}\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{{\kappa}_{y}}\frac{\partial }{\partial y}H_{z}-\frac{1}{{\kappa}_{z}}\frac{\partial }{\partial z}H_{y}+{\varphi}_{e_{xz}}-{\varphi}_{e_{xy}}\end{eqnarray}$ (19) where $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varphi}_{e_{xy}}=\frac{{\sigma}_{y}^{\prime }}{{\alpha}_{y}^{\prime }+j{\omega}{\varepsilon}_{y}^{\prime }}\frac{\partial }{\partial y}H_{z}\end{eqnarray}$ (20) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varphi}_{e_{xz}}=\frac{{\sigma}_{z}^{\prime }}{{\alpha}_{z}^{\prime }+j{\omega}{\varepsilon}_{z}^{\prime }}\frac{\partial }{\partial z}H_{y}.\end{eqnarray}$ (21) Here we rewrite (20)–(21) as $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle j{\omega}{\varepsilon}_{y}^{\prime }{\varphi}_{e_{xy}}+{\alpha}_{y}^{\prime }{\varphi}_{e_{xy}}={\sigma}_{y}^{\prime }\frac{\partial }{\partial y}H_{z}\end{eqnarray}$ (22) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle j{\omega}{\varepsilon}_{z}^{\prime }{\varphi}_{e_{xz}}+{\alpha}_{z}^{\prime }{\varphi}_{e_{xz}}={\sigma}_{z}^{\prime }\frac{\partial }{\partial z}H_{y}\end{eqnarray}$ (23) then transform (19), (22) and (23) into time domain $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}\frac{\partial }{\partial t}E_{x}+{\sigma}E_{x}=\frac{1}{{\kappa}_{y}}\frac{\partial }{\partial y}H_{z}-\frac{1}{{\kappa}_{z}}\frac{\partial }{\partial z}H_{y}+{\varphi}_{e_{xz}}-{\varphi}_{e_{xy}}\end{eqnarray}$ (24) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}_{y}^{\prime }\frac{\partial {\varphi}_{e_{xy}}}{\partial t}+{\alpha}_{y}^{\prime }{\varphi}_{e_{xy}}={\sigma}_{y}^{\prime }\frac{\partial }{\partial y}H_{z}\end{eqnarray}$ (25) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}_{z}^{\prime }\frac{\partial {\varphi}_{e_{xz}}}{\partial t}+{\alpha}_{z}^{\prime }{\varphi}_{e_{xz}}={\sigma}_{z}^{\prime }\frac{\partial }{\partial z}H_{y}.\end{eqnarray}$ (26)

Also for the sake of simplicity in the presentation and without loss of generality, the fields E _x is expanded in terms of scaling functions only in space domain and pulse functions in time domain. $\begin{eqnarray}\displaystyle E_{x}(r,t)=\mathop{\sum }_{i,j=-\infty }^{+\infty }E_{i+1/2,j,k}^{{\phi}x}(t){\phi}_{i+1/2}(x){\phi}_{j}(y){\phi}_{k}(z) & & \displaystyle\end{eqnarray}$ (27) where $E_{i+1/2,j,k}^{{\phi}x}$ is the coefficient for the field in terms of scaling functions which is equal to the corresponding field. The indexes i, j and k are the discrete space indices related to the space coordinates via x = iΔx, y = jΔy and z = kΔz, where Δx, Δy and Δz represent the space discretization intervals in x-, y- and z-direction.

With the wavelet-Galerkin scheme based on Daubechies’ compactly supported wavelets, the following equations can be obtained $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}\frac{\partial E_{i+\frac{1}{2},j,k}^{{\phi}x}(t)}{\partial t}+{\sigma}E_{i+\frac{1}{2},j,k}^{{\phi}x}(t)=\frac{1}{{\kappa}_{y}{\rm\Delta}y}\mathop{\sum }_{l}a(l)H_{i+\frac{1}{2},j+l+\frac{1}{2},k}^{{\phi}z}(t)-\frac{1}{{\kappa}_{z}{\rm\Delta}z}\mathop{\sum }_{l}a(l)H_{i+\frac{1}{2},j,k+l+\frac{1}{2}}^{{\phi}y}(t)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \hspace{132.0pt}+{\varphi}_{e_{xz}}-{\varphi}_{e_{xy}}\end{eqnarray}$ (28) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}_{y}^{\prime }\frac{{\varphi}_{e_{xy},i+1/2,j,k}^{{\phi}x}(t)}{\partial t}+{\alpha}^{\prime }{\varphi}_{e_{xy},i+1/2,j,k}^{{\phi}x}(t)=\frac{{\sigma}_{y}^{\prime }}{{\rm\Delta}y}\mathop{\sum }_{l}a(l)H_{i+\frac{1}{2},j+l+\frac{1}{2},k}^{{\phi}z}(t)\end{eqnarray}$ (29) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \frac{{\varepsilon}_{z}^{\prime }{\varphi}_{e_{xz},i+1/2,j,k}^{{\phi}x}(t)}{\partial t}+{\alpha}_{z}^{\prime }{\varphi}_{e_{xz},i+1/2,j,k}^{{\phi}x}(t)=\frac{{\sigma}_{z}^{\prime }}{{\rm\Delta}z}\mathop{\sum }_{l}a(l)H_{i+\frac{1}{2},j,k+l+\frac{1}{2}}^{{\phi}y}(t).\end{eqnarray}$ (30)

Discretizing by 5-stage and 4-order explicit symplectic integral in time domain, Eqs (28)–(30) can be written as $\begin{eqnarray}\displaystyle E_{i+1/2,j,k}^{{\phi}x,n+v/m} & = & \displaystyle \mathit{CA}_{i+1/2,j,k}E_{i+1/2,j,k}^{{\phi}x,n+(v-1)/m}+\mathit{CB}_{i+1/2,j,k}\left[\left(\frac{1}{{\kappa}_{y}{\rm\Delta}y}\mathop{\sum }_{l}a(l)H_{i+1/2,j+l+1/2,k}^{{\phi}z,n+v/m}\right.\right.\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \left.\left.-\frac{1}{{\kappa}_{z}{\rm\Delta}z}\mathop{\sum }_{l}a(l)H_{i+1/2,j,k+l+1/2}^{{\phi}y,n+v/m}\right)+{\psi}_{e_{xy},i+1/2,j,k}^{n+v/m}-{\psi}_{e_{xz},i+1/2,j,k}^{n+v/m}\right]\end{eqnarray}$ (31) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varphi}_{e_{xy},i+1/2,j,k}^{{\phi}x,n+v/m}=P_{j}^{y}{\varphi}_{e_{xy},i+1/2,j,k}^{{\phi}x,n+(v-1)/m}+\frac{Q_{j}^{y}}{{\rm\Delta}y}\mathop{\sum }_{l}a(l)H_{i+1/2,j+l+1/2,k}^{{\phi}z,n+v/m}\end{eqnarray}$ (32) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varphi}_{e_{_{xz}},i+1/2,j,k}^{{\phi}x,n+v/m}=P_{k}^{z}{\varphi}_{e_{_{xz}},i+1/2,j,k}^{{\phi}x,n+(v-1)/m}+\frac{Q_{k}^{z}}{{\rm\Delta}z}\mathop{\sum }_{l}a(l)H_{i,j+l+1/2}^{{\phi}y,n+v/m}\end{eqnarray}$ (33) where $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \mathit{CA}_{i+1/2,j,k}=\frac{2{\varepsilon}_{i+1/2,j,k}-{\sigma}_{i+1/2,j,k}d_{v}{\rm\Delta}t}{2{\varepsilon}_{i+1/2,j,k}+{\sigma}_{i+1/2,j,k}d_{v}{\rm\Delta}t}\end{eqnarray}$ (34) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \mathit{CB}_{i+1/2,j,k}=\frac{2d_{v}{\rm\Delta}t}{2{\varepsilon}_{i+1/2,j,k}+{\sigma}_{i+1/2,j,k}d_{v}{\rm\Delta}t}\end{eqnarray}$ (35) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle P_{j}^{y}=\frac{2{\varepsilon}_{yj}^{\prime }-{\alpha}_{yj}^{\prime }{\rm\Delta}t}{2{\varepsilon}_{yj}^{\prime }+{\alpha}_{yj}^{\prime }{\rm\Delta}t},\quad Q_{j}^{y}=\frac{2{\sigma}_{yj}^{\prime }{\rm\Delta}t}{2{\varepsilon}_{yj}^{\prime }+{\alpha}_{yj}^{\prime }{\rm\Delta}t}\end{eqnarray}$ (36) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle P_{k}^{z}=\frac{2{\varepsilon}_{zk}^{\prime }-{\alpha}_{zk}^{\prime }{\rm\Delta}t}{2{\varepsilon}_{zk}^{\prime }+{\alpha}_{zk}^{\prime }{\rm\Delta}t},\quad Q_{k}^{z}=\frac{2{\sigma}_{zk}^{\prime }{\rm\Delta}t}{2{\varepsilon}_{zk}^{\prime }+{\alpha}_{zk}^{\prime }{\rm\Delta}t}.\end{eqnarray}$ (37)

The other Ampere’s and Faraday’s equations in the ADE-PML can also be similarly obtained. Moreover, observing (31)–(33), it is seen that the explicit time-marching schemes for the field and the auxiliary variables are obtained at the (n + v∕m) time step. And the schemes satisfied the stability condition for the SMRTD scheme. It is also obvious that the ADE-PML implementation requires no more than two auxiliary variables per field component, which is less than that reported by previous implementation of this method [11]. Therefore, the ADE-PML method is more straightforward and memory saving. Moreover, we use perfectly electric conductor (PEC) walls to terminate the PML regions.

3.2. Result of the reflection error

To demonstrate ADE-PML termination of the SMRTD lattice, a 3-D lattice of the dimension Nx × Ny × Nz = 40 × 40 × 26 was used, surrounding a computational domain of the dimension 24 × 24 × 10 with 8-cell thick ADE-PML. The media with constitutive parameters ϵ_r = 8.0 and σ = 0.27 is introduced into the computational domain. The excitation was applied to the electric field component E _x at the center of the computation domain as follows: $\begin{eqnarray}\displaystyle E_{i,j}^{{\phi}z,n+v/m}=E_{i,j}^{{\phi}z,n+(v-1)/m}-\frac{((n+(v-1)/m){\rm\Delta}t-3t_{0})}{t_{0}}\exp \left(-\frac{((n+(v-1)/m){\rm\Delta}t-3t_{0})^{2}}{{t_{0}}^{2}}\right) & & \displaystyle\end{eqnarray}$ (38) where t ₀ = 2.0 ns. The space is discretized with a mesh with Δx = Δy = Δz = 0.05 m, and the time step is Δt = 60 ps. Within the ADE-PML layer, the constitutive parameters σ_i and 𝜅_i are scaled using a pth order polynomial scaling $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\sigma}({\rho})={\sigma}_{\max }\left(\frac{{\rho}}{d}\right)^{p}\end{eqnarray}$ (39) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\kappa}({\rho})=1+({\kappa}_{\max }-1)\left(\frac{{\rho}}{d}\right)^{p}\end{eqnarray}$ (40) where 𝜌 denotes the distance from the interface between the computational domain and the ADE-PML into the PML layer, d is the depth of the ADE-PML, and p is the order of the polynomial. A choice for σ_max can be expressed as $\begin{eqnarray}\displaystyle {\sigma}_{opt}=\frac{(m+1)}{150{\pi}\sqrt{{\varepsilon}_{r}}{\rm\Delta}} & & \displaystyle\end{eqnarray}$ (41) where Δ is the grid spacing along the normal axis and there is no difference between x-, y-, and z-direction in all computations in this article. Another ADE-PML parameter 𝛼 is constant through the ADE-PML. In this article, the reflection error was computed at the sampling point A that corresponding to E _x a cell from the ADE-PML interface, where electromagnetic wave is incident normally. In order to isolate the error due to the ADE-PML from grid dispersion error, a reference problem was also simulated: the same mesh is extended 100 cells out in all dimensions, leading to a 240 × 240 × 226 cell lattice. The error relative to reference solution was then computed as $\begin{eqnarray}\displaystyle \mathit{error}_{\text{dB}}=20\log _{10}\frac{|E_{x}(t)-E_{x_{ref}}(t)|}{\max |E_{x_{ref}}(t)|} & & \displaystyle\end{eqnarray}$ (42) where E _x(t) represents the time-dependent discrete field computed within the working volume, E _{x
_ref}(t) represents the same discrete field computed by the reference problem.

Defining that k = σ_max∕σ_opt. It is instructive to observe the maximum reflection error experienced by the ADE-PML as a function of the constitutive parameters 𝜅_max, k and 𝛼. Figure 4 shows the contour plots of the maximum relative error over 1000 time steps at point A versus 𝜅_max and k with 𝛼 = 0.001. It can be demonstrated that as low as −120 dB maximum error is achieved.

Fig. 4.

Contour plots of the maximum relative error for the first 900 time steps (CFS-PML).

Figure 5 does the same with the same 𝜅_max and k as Fig. 4 but for the APML. It can be seen that the maximum error is on the order of 110 dB. Compared with the APML, an improvement of 10 dB is obtained with the ADE-PML. Moreover, the optimal error is realized over a much broader range of 𝜅_max and k, making these values easier to predict.

Fig. 5.

Contour plots of the maximum relative error for the first 900 time steps (APML).

Fig. 6.

RCS of the dielectric sphere.

4. Application of the C-SMRTD method

4.1. RCS of dielectric sphere

The radius of the sphere is 0.5 m, and the relative permittivity ϵ_r = 4, relative permeability μ_r = 1, conductivity σ = 0. The wavelength of the incident sine wave is 1 m, and the incident wave propagates along the z-direction, the electric field polarizes along the x-direction. The mesh size and the Courant-Friedrichs-Lewy (CFL) number is 0.05 m and 0.35 for C-SMRTD and MRTD method. The mesh size and the CFL number is 0.02 m and 0.5 for CFDTD method. In Fig. 6, the RCS of dielectric sphere calculated by Mie series, MRTD with step approximation, CFDTD and C-SMRTD are compared.

It can be seen from Fig. 6 that the accuracy of the C-SMRTD has been significantly improved compared with the MRTD with step approximation, and the results of C-SMRTD and CFDTD are in good agreement with the analytical curve, but as shown in Table 1, it is also obvious that the C-SMRTD has an advantage over CFDTD at saving the CPU time and memory.

Table 1
CPU time and memory for different methods

Methods Cell lattice CPU time/s (3000Δt)

C-SMRTD 64 × 64 × 64 1364.23

CFDTD 80 × 80 × 80 1786.56

Methods	Cell lattice	CPU time/s (3000Δt)
C-SMRTD	64 × 64 × 64	1364.23
CFDTD	80 × 80 × 80	1786.56

4.2. Field cross-section distribution of the TE₁₁ propagation model in vaulted tunnel

The compendious model of the vaulted tunnel is illustrated in Fig. 7. The dimensions of the tunnel cross-section are shown in Fig. 7(a). The source is placed near the PML. For the purposes of this study, constitutive parameters for soil were assumed, giving σ_s = 0.035, ϵ_s = 8.5. The constitutive parameters for the wall and arris of the straight part can be defined as $\begin{eqnarray}\displaystyle \text{wall:}∼{\varepsilon}_{eq}=\frac{{\varepsilon}_{s}+{\varepsilon}_{0}}{2},{\sigma}_{eq}=\frac{{\sigma}_{s}}{2},\quad \text{arris:}∼{\varepsilon}_{eq}=\frac{3{\varepsilon}_{0}{\varepsilon}_{r}+{\varepsilon}_{0}}{4},{\sigma}_{eq}=\frac{3{\sigma}_{s}}{4}. & & \displaystyle \nonumber\end{eqnarray}$

The constitutive parameters for the wall and arris of the crooked part can be defined via the proposed conformal technique, and Fig. 8 shows the conformed vault of the tunnel.

Fig. 7.

The computational model of the tunnel.

Fig. 8.

The conformed vault of the tunnel.

In the waveguide system [18], the excitation source is usually introduced robustly according to propagation model such as TE₁₁, TM₁₁, etc. The way that the excitation sources induced in the waveguide system can still be employed here as follows $\begin{eqnarray}\displaystyle E_{tan}^{n+1}(i,j,k_{s})=E_{tan}^{n}(i,j,k_{s})+f(i,j,k_{s})g(t) & & \displaystyle\end{eqnarray}$ (43) where the subscript ‘tan’ denotes the E-field distributed in a transverse cross section at z = k _sΔz of the tunnel structure in Fig. 7, f (i, j, k _s) is the function of the field distribution and g (t) refer to the time function determine the bandwidth of the sources. Here we set f (i, j, k _s) as the model of TE₁₁ propagation in waveguide approximately, and suppose that the long side of the tunnel is a and the short side is b, then the TE₁₁ propagation model is defined as $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}E_{x}={\displaystyle \frac{j{\omega}{\pi}{\mu}}{k_{c}^{2}b}}H_{0}\cos \left({\displaystyle \frac{{\pi}}{a}}x\right)\sin \left({\displaystyle \frac{{\pi}}{b}}y\right)e^{-jk_{z}z}\\[15.0pt] E_{y}=-{\displaystyle \frac{j{\omega}{\pi}{\mu}}{k_{c}^{2}a}}H_{0}\sin \left({\displaystyle \frac{{\pi}}{a}}x\right)\cos \left({\displaystyle \frac{{\pi}}{b}}y\right)e^{-jk_{z}z}\\[8.0pt] E_{z}=0\\[3.0pt] H_{x}={\displaystyle \frac{jk_{z}{\pi}}{k_{c}^{2}a}}H_{0}\sin \left({\displaystyle \frac{{\pi}}{a}}x\right)\cos \left({\displaystyle \frac{{\pi}}{b}}y\right)e^{-jk_{z}z}\\[15.0pt] H_{y}={\displaystyle \frac{jk_{z}{\pi}}{k_{c}^{2}b}}H_{0}\cos \left({\displaystyle \frac{{\pi}}{a}}x\right)\sin \left({\displaystyle \frac{{\pi}}{b}}y\right)e^{-jk_{z}z}\\[12.0pt] H_{z}=H_{0}\cos \left({\displaystyle \frac{{\pi}}{a}}x\right)\cos \left({\displaystyle \frac{{\pi}}{b}}y\right)e^{-jk_{z}z}\end{array}\right. & & \displaystyle\end{eqnarray}$ (44) Where g (t) is set to be a differential Gaussian electric pulse, g (t) = E ₀(t − t ₀)exp(4π(t − t ₀)²∕τ²), τ = 3.0 ns, E ₀ = 1000 V/m, t ₀ = τ and a = 6 m, b = 4 m.

Fig. 9.

The cross-section distributions of the fields.

Figure 9(a)–(f) denote the cross-section distributions of the fields of E and H in the tunnel. Observing the results, it is seen that the mainly characteristics of the fields’ distributions are generally similar with that in the perfect waveguide. Though the field E _z is not zero, it is very small compared with the other electric fields. And the distribution characteristics for E _x and H _y, E _y and H _x are similar, which can be validated from Eq. (44). Taking E _x as example, E _x has two peaks along the x-direction, one of which is positive, the other one is negative. Moreover, the E _z, field basically distributes as the stationary wave along the x and y coordinates, and the energy of E _z field is mainly located in the middle area of the cross-section.

5. Conclusion

In this paper, the C-SMRTD method is derived in detail based on the concept of equivalent electromagnetic parameters and conformal grid technology. The proposed ADE-PML for SMRTD is more straightforward to implement compared with conventional APML. Numerical examples show that maximum errors on the order of −120 dB were recorded for the ADE-PML, compared to −110 dB for the APML. The RCS of a dielectric sphere is calculated to verify the effectiveness of the C-SMRTD method. The proposed method shows good agreement with the results derived by analytical methods, and has higher accuracy compared with the MRTD method with step approximation. It is also shown that the proposed method has an advantage over CFDTD at saving the CPU time and memory. Furthermore, the calculation model of a vault tunnel is established by using the C-SMRTD method, the characteristics of the field cross-section distributions of the TE₁₁ propagation model in the vaulted tunnel are obtained and analysed.

Footnotes

Acknowledgements

This work was supported by the Key Science and Technology Projects of the Education Department of Jilin Province (Grant No. JJKH20210178KJ); Scientific and Technological Project of Jilin Engineering Normal University (Grant No. BSKJ202008); Youth Fund of National Natural Science Foundation of China (Grant no. 61801511) and the Youth Fund of Jiangsu Natural Science Foundation (Grant no. BK20180580).

References

Wei

Huang

Z.X.

and Wang

, A novel symplectic multi-resolution time-domain scheme for electromagnetic simulations, IEEE Microwave and Wireless Components Letters 23(4) (2013), 175–177, doi:10.1109/lmwc.2013.2247990.

Jurgens

T.G.

Taflove

Umashankar

, Finite-difference time-domain modeling of curved surface, IEEE Transactions on Antennas & Propagation 40(4) (1992), 357–366, doi:10.1109/8.138836.

Zheng

Zhang

, Wide-band modeling on-chip spiral inductors using frequency-dependent conformal ADI-FDTD method, IEEE Access 7 (2019), 184940–184949, doi:10.1109/ACCESS.2019.2960284.

Cabello

M.R.

Angulo

L.D.

Alvarez

, A new efficient and stable 3D conformal FDTD, IEEE Microwave and Wireless Components Letters 26(8) (2016), 553–555, doi:10.1109/LMWC.2016.2588579.

Muraoka

Inoue

Sekine

, A hybrid implicit-explicit and conformal (HIE/C) FDTD method for efficient electromagnetic simulation of nonorthogonally aligned thin structures, IEEE Transactions on Electromagnetic Compatibility 57(3) (2015), 505–512, doi:10.1109/TEMC.2013.2273852.

Zagorodnov

I.A.

Schuhmann

and Weiland

, A uniformly stable conformal FDTD-method in Cartesian grids, International Journal of Numerical Modelling Electronic Networks Devices and Fields 16(2) (2003), 127–141, doi:10.1002/jnm.488.

W.H.

and Mittra

, A conformal finite difference time domain technique for modeling curved dielectric surfaces, IEEE Microwave and Wireless Components Letters 11(1) (2002), 25–27, doi:10.1109/7260.905957.

Okoniewski

Okoniewska

and Stuchly

M.A.

, Three-dimensional subgridding algorithm for FDTD, IEEE Transactions on Antennas and Propagation 45(3) (2002), 422–429, doi:10.1109/8.558657.

Dey

and Mittra

, A conformal finite-difference time-domain technique for modeling cylindrical dielectric resonators, IEEE Transactions on Microwave Theory and Techniques 47(9) (1999), 1737–1739, doi:10.1109/22.788616.

10.

Wang

Yin

W.Y.

and Liu

Q.H.

, FDTD (2, 4)-compatible conformal technique for treatment of dielectric surfaces, Electronics Letters 45(3) (2009), 146–147, doi:10.1049/el:20093215.

11.

Cao

and Chen

, Application of an anisotropic perfectly matched layer absorber for open boundary truncation in the multiresolution time domain scheme, IEEE Trans. Antennas Propagat. 51(2) (2003), 350–357, doi:10.1109/TAP.2003.809068.

12.

Mukherjee

, Numerical solution in FDTD for absorbing boundary condition over dielectric surfaces, Journal of Advance Research in Scientific Computing, IASR 4(1) (2012), 13–23.

13.

Mukherjee

and Vishwakarma

D.K.

, Application of finite difference time domain to calculate the transmission coefficient of an electromagnetic wave impinging perpendicularly on a dielectric interface with modified MUR-I ABC, Defence Science Journal 62(4) (2012), 228–235, doi:10.14429/dsj.62.792.

14.

Cheong

Y.W.

and Lee

Y.M.

, Wavelet-Galerkin scheme of time-dependent inhomogeneous electromagnetic problems, IEEE Microwave and Guided Wave Letters 9(8) (1999), 297–299, doi:10.1109/75.779907.

15.

Fujii

and Hoefer

W.J.R.

, Dispersion of time domain wavelet Galerkin method based on Daubechies’ compactly supported scaling functions with three and four vanishing moments, IEEE Microwave and Guided Wave Letters 10(4) (2000), 125–127, doi:10.1109/75.846920.

16.

Liu

Y.W.

Chen

Y.W.

Chen

and Xu

, A cylindrical MRTD algorithm with PML and Quasi-PML, IEEE Trans. Microwave Theory Tech. 61(3) (2013), 1006–1017, doi:10.1109/TMTT.2013.2243745.

17.

Daubechies

, Ten Lectures on Wavelets, SIAM, Philadelphia, PA, 1992.

18.

Taflove

, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech House, Norwood, MA, 1995.

A conformal symplectic multi-resolution time-domain method with ADE-PML

Abstract

Keywords

1. Introduction

2. C-SMRTD method

2.1. Multi domain decomposition of electric field component

3.1. Formulation

4.1. RCS of dielectric sphere

Table 1 CPU time and memory for different methods Methods Cell lattice CPU time/s (3000Δt) C-SMRTD 64 × 64 × 64 1364.23 CFDTD 80 × 80 × 80 1786.56

Footnotes

Acknowledgements

References

Table 1
CPU time and memory for different methods

Methods Cell lattice CPU time/s (3000Δt)

C-SMRTD 64 × 64 × 64 1364.23

CFDTD 80 × 80 × 80 1786.56