An efficient one-dimensional time-domain algorithm for obliquely incident plane wave in BOR-SMRTD

Abstract

In this paper, the concept of body of revolution (BOR) finite-difference time-domain method time-domain (BOR-FDTD) method is applied to symplectic multi-resolution time-domain (SMRTD) algorithm to propose a BOR-SMRTD method. In order to solve the problem that the obliquely incident pulse wave cannot be directly introduced into BOR-SMRTD method, an efficient one-dimensional time-domain algorithm for propagating obliquely incident pulse plane wave in BOR-SMRTD method is also proposed. The numerical results show that the proposed algorithm runs 51 times faster than the algorithm based on fast Fourier transform (FFT) completely. Moreover, the scattering problem of a perfect electric conductor (PEC) cylinder with finite length is calculated to verify the effectiveness of the algorithm. The surface current density distribution along a generatrix of the cylinder at a giving frequency was extracted, and the results are in good agreement with that obtained with the methods of moment (MoM) and BOR-FDTD. Finally, the radar cross section (RCS) of a PEC sphere was calculated, and the results agree well with the theoretical values.

Keywords

Body of revolution symplectic multi-resolution time-domain FFT

1. Introduction

Body of revolution finite-difference time-domain (BOR-FDTD) method [1–10] is a special numerical method for the efficient analysis of electromagnetic field with rotationally symmetric structure. This method greatly reduces the CPU time and the requirement of computer memory, and is an effective method to analyze the problem of bodies of revolution. However, due to the limitation of stability conditions, BOR-FDTD is mainly used to the analysis of rotational symmetry and enclosed cavity in the case of small electrical size and small incident angle of plane wave. Compared with the traditional FDTD and MRTD, SMRTD has higher numerical stability and better numerical dispersion characteristics [11]. Therefore, a BOR-SMRTD method is proposed in this paper.

When analyzing the electromagnetic scattering problems, pulse wave is usually used as the excitation source in order to obtain the wide frequency band data through a time domain calculation. However, the obliquely incident pulse plane wave cannot be directly introduced into BOR-SMRTD because of its non-rotational symmetry. The plane wave must first be expanded into a series of cylindrical waves in the frequency domain [12,13], then the time domain component can be obtained by FFT, and the incident wave can be introduced by total field/scattered field surface (TFSF surface). However, the FFT is needed for each time-domain incident field component at the TFSF surface, which will consume a lot of CPU time. In order to avoid a large number of Fourier transforms, an effective one-dimensional time-domain algorithm for obliquely incident plane waves in BOR-SMRTD is proposed. With the proposed algorithm, the differential equation of two-dimensional BOR is transformed into one-dimensional equation by variable substitution according to the time delay characteristic of incident plane wave in the symmetry axis direction. The numerical results show that the algorithm is 50 times faster than the algorithm using FFT. To validate the implementation of the algorithm, the scattering problem of a metal cylinder with finite length is calculated. The surface current density distribution along a generatrix of the cylinder at a certain frequency is extracted and compared with the results obtained by MoM and BOR-FDTD, and the results are in good agreement. Then, the RCS of a metal sphere is calculated, and the results are in good agreement with the theoretical values.

2. BOR-SMRTD method

Firstly, the electric field and magnetic field components are expanded into Fourier series. $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \mathbf{E}=\mathop{\sum }_{m=0}^{\infty }(\mathbf{e}_{u}\cos m{\phi}+\mathbf{e}_{v}\sin m{\phi})\end{eqnarray}$ (1) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \mathbf{H}=\mathop{\sum }_{m=0}^{\infty }(\mathbf{h}_{u}\cos m{\phi}+\mathbf{h}_{v}\sin m{\phi}).\end{eqnarray}$ (2)

Assuming that the medium is homogeneous isotropic, and the coordinate variables in cylindrical coordinate system are 𝜌, 𝜙 and z, then the e _u, e _v, h _u and h _v are Fourier coefficients, which are only related to coordinates 𝜌, z and time t, and are required in BOR-SMRTD. 𝜙 is azimuth, m is azimuth modulus. Where subscript u denotes even function term, which is multiplied by cos mφ, and v denotes odd function term, which is multiplied by sin m𝜙.

By substituting Eqs (1) and (2) into Maxwell’s equations in cylindrical coordinates, the following equations can be obtained $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\mu}\frac{\partial h_{{\rho}}}{\partial t}=\pm {\displaystyle \frac{m}{{\rho}}}e_{z}+{\displaystyle \frac{\partial e_{{\phi}}}{\partial z}}\end{eqnarray}$ (3) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\mu}\frac{\partial h_{{\phi}}}{\partial t}=-{\displaystyle \frac{\partial e_{{\rho}}}{\partial z}}+\frac{\partial e_{z}}{\partial {\rho}}\end{eqnarray}$ (4) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\mu}{\displaystyle \frac{\partial h_{z}}{\partial t}}=-{\displaystyle \frac{1}{{\rho}}}{\displaystyle \frac{\partial ({\rho}e_{{\phi}})}{\partial {\rho}}}\mp {\displaystyle \frac{m}{{\rho}}}e_{{\rho}}\end{eqnarray}$ (5) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}{\displaystyle \frac{\partial e_{{\rho}}}{\partial t}}=\pm \frac{\text{m}}{{\rho}}h_{z}-\frac{\partial h_{{\phi}}}{\partial z}\end{eqnarray}$ (6) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}{\displaystyle \frac{\partial e_{{\phi}}}{\partial t}}=\frac{\partial h_{{\rho}}}{\partial z}-{\displaystyle \frac{\partial h_{z}}{\partial {\rho}}}\end{eqnarray}$ (7) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}{\displaystyle \frac{\partial e_{z}}{\partial t}}={\displaystyle \frac{1}{{\rho}}}\frac{\partial ({\rho}h_{{\phi}})}{\partial {\rho}}\mp {\displaystyle \frac{m}{{\rho}}}h_{{\rho}}.\end{eqnarray}$ (8)

The field discretization in space domain of the BOR-SMRTD method is consistent with that of the MRTD method, taking h _𝜌 and e _𝜌 as examples, the field expansion can be expressed as follows: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle h_{{\rho}}(\mathbf{r},t)=\mathop{\sum }_{i,j,k,n=-\infty }^{+\infty }h_{i,j+1/2,k+1/2}^{{\varphi}{\rho}}(t){\varphi}_{i}({\rho}){\varphi}_{j+1/2}({\phi}){\varphi}_{k+1/2}(z)\end{eqnarray}$ (9) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle e_{{\rho}}(\mathbf{r},t)=\mathop{\sum }_{i,j,k,n=-\infty }^{+\infty }e_{i+1/2,j,k}^{{\phi}{\rho}}(t){\varphi}_{i+1/2}({\rho}){\varphi}_{j}({\phi}){\varphi}_{k}(z)\end{eqnarray}$ (10) where $h_{i,j+1/2,k+1/2}^{{\varphi}{\rho}}$ and $e_{i+1/2,j,k}^{{\phi}{\rho}}$ are the coefficients for the fields expansions in terms of scaling functions. The indexes i, j and k are the discrete space indices related to the space coordinates via 𝜌 = iΔ𝜌, 𝜙 = jΔ𝜙 and z = kΔz, where, Δ𝜌, Δ𝜙 and Δz represent the space discretization intervals in 𝜌-, 𝜙- and z-direction. The function φ(v), where v = 𝜌, 𝜙, z, is defined as Daubechies’ scaling function. With the wavelet-Galerkin scheme based on Daubechies’ compactly supported wavelets [14], the following equations can be obtained: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\mu}\frac{\partial h_{i,k+1/2}^{{\varphi}{\rho}}(t)}{\partial t}=\pm {\displaystyle \frac{m}{i{\rm\Delta}{\rho}}}e_{i,k+1/2}^{{\varphi}z}(t)+\mathop{\sum }_{l}a(l){\displaystyle \frac{e_{i,k+l+1}^{{\varphi}{\phi}}(t)}{{\rm\Delta}z}}\end{eqnarray}$ (11) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}{\displaystyle \frac{\partial e_{i+1/2,k}^{{\varphi}{\rho}}(t)}{\partial t}}=\pm {\displaystyle \frac{m}{(i+\frac{1}{2}){\rm\Delta}{\rho}}}{\displaystyle \frac{h_{i+1/2,k}^{{\varphi}z}(t)}{{\rm\Delta}y}}-\mathop{\sum }_{l}a(l){\displaystyle \frac{h_{i+1/2,k+l+1/2}^{{\varphi}{\phi}}(t)}{{\rm\Delta}z}}.\end{eqnarray}$ (12)

The Maxwell’s equations (MEs) can be viewed as a Hamiltonian system with an infinite dimensional, and the Hamiltonian formulation can be written as $\begin{eqnarray}\displaystyle {\Pi}(\mathbf{H},\mathbf{E})={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{{\varepsilon}}}\mathbf{H}\bullet {\nabla}\times \mathbf{H}+{\displaystyle \frac{1}{{\mu}}}\mathbf{E}\bullet {\nabla}\times \mathbf{E}\right). & & \displaystyle\end{eqnarray}$ (13) The Hamiltonian system corresponding to Eq. (13) is as follows $\begin{eqnarray}\displaystyle \frac{\partial }{\partial t}\left(\begin{array}{@{}c@{}}\mathbf{H}\\ \boldsymbol{ E}\end{array}\right)=(\mathbf{A}+\mathbf{B})\left(\begin{array}{@{}c@{}}\mathbf{H}\\ \boldsymbol{ E}\end{array}\right) & & \displaystyle\end{eqnarray}$ (14) where $\begin{eqnarray}\displaystyle \mathbf{A}=\left(\begin{array}{@{}cc@{}}\{0\}_{3\times 3} & {\mu}^{-1}\mathbf{U}_{3\times 3}\\ \{0\}_{3\times 3} & \{0\}_{3\times 3}\end{array}\right)_{6\times 6},\quad \mathbf{B}=\left(\begin{array}{@{}cc@{}}\{0\}_{3\times 3} & \{0\}_{3\times 3}\\ {\varepsilon}^{-1}\mathbf{V }_{3\times 3} & \{0\}_{3\times 3}\end{array}\right)_{6\times 6} & & \displaystyle\end{eqnarray}$ (15) and $\begin{eqnarray}\displaystyle \mathbf{U}=\left(\begin{array}{@{}ccc@{}}0 & {\displaystyle \frac{\partial }{\partial z}} & \pm {\displaystyle \frac{m}{{\rho}}}\\ -{\displaystyle \frac{\partial }{\partial z}} & 0 & {\displaystyle \frac{\partial }{\partial {\rho}}}\\ \mp {\displaystyle \frac{\text{m}}{{\rho}}} & -{\displaystyle \frac{1}{{\rho}}}{\displaystyle \frac{\partial ({\rho})}{\partial {\rho}}} & 0\end{array}\right),\quad \mathbf{V }=\left(\begin{array}{@{}ccc@{}}0 & -{\displaystyle \frac{\partial }{\partial z}} & \pm {\displaystyle \frac{m}{{\rho}}}\\ {\displaystyle \frac{\partial }{\partial z}} & 0 & -{\displaystyle \frac{\partial }{\partial {\rho}}}\\ \mp {\displaystyle \frac{m}{{\rho}}} & {\displaystyle \frac{1}{{\rho}}}{\displaystyle \frac{\partial ({\rho})}{\partial {\rho}}} & 0\end{array}\right). & & \displaystyle\end{eqnarray}$ (16) From time t = 0 to t = Δt, the operator approximates exp(Δt (A + B)) as [15 ] $\begin{eqnarray}\displaystyle \exp ({\rm\Delta}t(\mathbf{A}+\mathbf{B}))=\mathop{\prod }_{v=1}^{q}\exp (d_{v}{\rm\Delta}tB)\exp (c_{v}{\rm\Delta}tA)+O({\rm\Delta}t^{p+1}) & & \displaystyle\end{eqnarray}$ (17) where c _v and d _v are the symplectic integrator propagator time coefficients, q and p denote the stage number of the propagator, and the order of approximation, respectively. In this work, we adopt the five-stage and four-order symplectic integrator propagator as an example, the values of c _v and d _v are as follows $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle c_{1}=c_{5}=0.1786178958,\quad c_{2}=c_{4}=-0.0662645827,\quad c_{3}=0.7752933737,\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle d_{1}=d_{4}=0.7123418311,\quad d_{2}=d_{3}=-0.2123418311,\quad d_{5}=0.\nonumber\end{eqnarray}$

Also taking h _𝜌 and e _𝜌 as examples, the BOR-SMRTD iteration equations at the v-stage are given by $\begin{eqnarray}\displaystyle h_{i,k+1/2}^{{\varphi}{\rho},n+v/q} & = & \displaystyle h_{i,k+1/2}^{{\varphi}{\rho},n+(v-1)/q}+{\displaystyle \frac{{\rm\Delta}t}{{\mu}}}\left[\mathop{\sum }_{l}a(l){\displaystyle \frac{e_{i,k+l+1}^{{\varphi}{\phi},n+(v-1)/q}}{{\rm\Delta}z}}\pm {\displaystyle \frac{m}{i{\rm\Delta}{\rho}}}e_{i,k+1/2}^{{\varphi}z,n+(v-1)/q}\right]\end{eqnarray}$ (18) $\begin{eqnarray}\displaystyle e_{i+1/2,k}^{{\varphi}{\rho},n+v/q} & = & \displaystyle e_{i+1/2,k}^{{\varphi}{\rho},n+(v-1)/q}-{\displaystyle \frac{{\rm\Delta}t}{{\varepsilon}}}\left[\mathop{\sum }_{l}a(l){\displaystyle \frac{h_{i+1/2,k+l+1/2}^{{\varphi}{\phi},n+v/q}}{{\rm\Delta}z}}\mp {\displaystyle \frac{m}{(i+{\textstyle \frac{1}{2}}){\rm\Delta}{\rho}}}{\displaystyle \frac{h_{i+1/2,k}^{{\varphi}z,n+v/q}}{{\rm\Delta}y}}\right].\end{eqnarray}$ (19)

It should be noted that the components on the axis need special treatment due to singularity [14,16].

3. Cylindrical wave expansion of obliquely incident plane waves

Fig. 1.

The schematic diagram of incident plane waves.

As shown in Fig. 1, the incident plane wave is assumed to propagate in the direction, θ = π −θ_iφ = 0, and to be polarized in a direction making an angle θ_p with the unit vector $-\hat{{\theta}}$ . To simplify the expansion of the plane wave, let us separate the wave into a wave polarized parallel to the x–z plane of incidence, and a wave polarized perpendicular to this plane, as follows: $\begin{eqnarray}\displaystyle \mathbf{E}_{i}=\mathbf{E}_{i}^{//}+\mathbf{E}_{i}^{\bot }=-E_{i}^{//}\hat{{\theta}}+E_{i}^{\bot }{\hat{y}}=-\cos {\theta}_{p}E_{i}\hat{{\theta}}+\sin {\theta}_{p}E_{i}{\hat{y}}. & & \displaystyle\end{eqnarray}$ (20)

Referring to the method in reference [13] and the transformation relationship between rectangular coordinates and cylindrical coordinates, taking E _𝜌i as example, the incident plane wave can be expanded into the following form:

Parallel polarization wave: $\begin{eqnarray}\displaystyle E_{{\rho}i}^{//}=-E_{0}\cos {\theta}_{p}\cos {\theta}_{i}G({\omega})\text{e}^{-jkz\cos {\theta}_{i}}\mathop{\sum }_{m=0}^{\infty }{\varepsilon}_{m}j^{m+1}J_{m}^{\prime }(k{\rho}\sin {\theta}_{i})\cos m{\phi}. & & \displaystyle\end{eqnarray}$ (21) Vertical polarization wave: $\begin{eqnarray}\displaystyle E_{{\rho}i}^{\bot }=-E_{0}\sin {\theta}_{p}G({\omega})\text{e}^{-jkz\cos {\theta}_{i}}\mathop{\sum }_{m=0}^{\infty }m{\varepsilon}_{m}j^{m+1}\frac{J_{m}(k{\rho}\sin {\theta}_{i})}{k{\rho}\sin {\theta}_{i}}\sin m{\phi} & & \displaystyle\end{eqnarray}$ (22) where $G({\omega})={\mathcal{F}}[g(t)]$ is the Fourier transform of incident wave g (t), and ϵ_m = 1 at m = 0, ϵ_m = 2 at m≠ 0, and J _m is the Bessel function of order m. In theory, m tends to be infinite in all the above formulas, but in practice, according to the asymptotic characteristics of Bessel function and the requirement of calculation accuracy, only the first M terms are taken as approximation values. Generally, the series are truncated at M = kR _max sin θ_i +6, and R _max is the maximum radius of calculation domain [12].

Taking parallel polarized wave as an example, the time-domain expressions of Eq. (21) can be obtained by inverse Fourier transform $\begin{eqnarray}\displaystyle E_{{\rho}}^{i}=-E_{0}\cos {\theta}_{p}\cos {\theta}_{i}\mathop{\sum }_{m=0}^{\infty }A_{m}\left({\rho},t-{\displaystyle \frac{z\cos {\theta}_{i}}{c}}\right)\cos m{\phi} & & \displaystyle\end{eqnarray}$ (23) where $\begin{eqnarray}\displaystyle A_{m}({\rho},t)={\mathcal{F}}^{-1}[G({\omega}){\varepsilon}_{m}j^{m+1}J_{m}^{\prime }(k{\rho}\sin {\theta}_{i})]. & & \displaystyle\end{eqnarray}$ (24)

4. One-dimensional time-domain algorithm for obliquely incident plane wave

It can be seen from Eqs (23) that the calculation of the incident field component at the TFSF surface requires a large number of FFTs, which will consume a lot of CPU time. It is noted that each component of the incident wave field has the same time term t − zcosθ_i∕c, which means that the waveform of the field component is the same at all points in the z direction, only with a time delay. Therefore, as shown in Fig. 2, as long as the incident wave at z = k ₁Δz is calculated by FFT, the field components at other points can be obtained by a corresponding time delay.

Fig. 2.

Scattering calculation diagram of BOR-FDTD.

Equations (23) can be further abbreviated as: $\begin{eqnarray}\displaystyle E_{{\rho}}^{i}=\mathop{\sum }_{m=0}^{\infty }e_{{\rho}}\left({\rho},t-{\displaystyle \frac{\cos {\theta}_{i}}{c}}z\right)\cos m{\phi}. & & \displaystyle\end{eqnarray}$ (25) Let t ^′ = t − zcos θ_i∕c = t −𝛼z, then for each field component f, we can get $\begin{eqnarray}\displaystyle {\displaystyle \frac{\partial f}{\partial t}}={\displaystyle \frac{\partial f}{\partial t^{\prime }}},\quad {\displaystyle \frac{\partial f}{\partial z}}=-{\alpha}{\displaystyle \frac{\partial f}{\partial t^{\prime }}}. & & \displaystyle\end{eqnarray}$ (26) Substituting Eq. (26) into Eqs (3)–(8), the following results can be obtained $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\mu}{\displaystyle \frac{\partial h_{{\rho}}}{\partial t}}=\frac{m}{{\rho}}e_{z}-{\alpha}\frac{\partial e_{{\phi}}}{\partial t}\end{eqnarray}$ (27) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\mu}\frac{\partial h_{{\phi}}}{\partial t}={\alpha}\frac{\partial e_{{\rho}}}{\partial t}+\frac{\partial e_{z}}{\partial {\rho}}\end{eqnarray}$ (28) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}\frac{\partial e_{{\rho}}}{\partial t}={\alpha}\frac{\partial h_{{\phi}}}{\partial t}\end{eqnarray}$ (29) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}\frac{\partial e_{{\phi}}}{\partial t}=-{\alpha}\frac{\partial h_{{\rho}}}{\partial t}\end{eqnarray}$ (30) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle {\varepsilon}\frac{\partial e_{z}}{\partial t}={\displaystyle \frac{1}{{\rho}}}\frac{\partial ({\rho}h_{{\phi}})}{\partial {\rho}}-{\displaystyle \frac{m}{{\rho}}}h_{{\rho}}.\end{eqnarray}$ (31) Obviously, Eqs (27)–(31) are one-dimensional partial differential equations about 𝜌. Since the h _z component of the parallel polarization wave is 0, there are only five equations.

Following the procedure of Section 2, Eqs (27)–(31) can be discretized as $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle h_{i}^{{\varphi}{\rho},n+v/q}(i)=h_{i}^{{\varphi}{\rho},n+(v-1)/q}+\frac{m{\rm\Delta}t}{\text{CA}i{\rm\Delta}{\rho}}e_{i}^{{\varphi}z,n+(v-1)/q}\end{eqnarray}$ (32) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle h_{i+1/2}^{{\varphi}{\phi},n+v/q}=h_{i+1/2}^{{\varphi}{\phi},n+(v-1)/q}+\frac{{\rm\Delta}t}{\text{CA}{\rm\Delta}{\rho}}\mathop{\sum }_{l}a(l)e_{i+l+1}^{{\varphi}z,n+(v-1)/q}\end{eqnarray}$ (33) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle e_{i+1/2}^{{\varphi}{\rho},n+v/q}=e_{i+1/2}^{{\varphi}{\rho},n+(v-1)/q}+\frac{{\alpha}{\rm\Delta}t}{2\text{CA}{\varepsilon}{\rm\Delta}{\rho}}\mathop{\sum }_{l}a(l)\left(e_{i+l+1}^{{\varphi}z,n+v/q}+e_{i+l+1}^{{\varphi}z,n+(v-1)/q}\right)\end{eqnarray}$ (34) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle e_{i}^{{\varphi}{\phi},n+v/q}=e_{i}^{{\varphi}{\phi},n+(v-1)/q}-\frac{m{\alpha}{\rm\Delta}t}{2\text{CA}{\varepsilon}i{\rm\Delta}{\rho}}\left(e_{i}^{{\varphi}z,n+v/q}+e_{i}^{{\varphi}z,n+(v-1)/q}\right)\end{eqnarray}$ (35) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle e_{i}^{{\varphi}z,n+v/q}=e_{i}^{{\varphi}z,n+(v-1)/q}+{\displaystyle \frac{{\rm\Delta}t}{{\varepsilon}i{\rm\Delta}{\rho}}}\mathop{\sum }_{l}a(l)\left((i+l+1/2)h_{i+l+1/2}^{{\varphi}{\phi},n+(v-1)/q}-mh_{i}^{{\varphi}{\rho},n+(v-1)/q}\right)\end{eqnarray}$ (36) where CA = μ −𝛼²∕ϵ.

It can be seen that the Eqs (34) and (35) are implicit and cannot be calculated directly. But as long as each component is calculated in h _𝜌 → h _𝜙 → e _z → e _𝜌 → e _𝜙 order, this problem can be easily solved. The whole calculation process only needs two FFTs, which greatly improves the computational efficiency.

5. Numerical results

5.1. Validation of the proposed algorithm

Fig. 3.

Calculation results of the pulse wave.

Firstly, in order to verify the effectiveness of the proposed one-dimensional time-domain algorithm, the incident electric field component at z = k ₁Δz is calculated by using the proposed algorithm and FFT method respectively. A Gaussian pulse excitation is used here $\begin{eqnarray}\displaystyle E_{i}(t)=E_{0}\exp \left(-{\displaystyle \frac{4{\pi}(t-t_{0})^{2}}{{\tau}^{2}}}\right) & & \displaystyle\end{eqnarray}$ (37) where τ = 5 ×10⁻⁸ s, t ₀ = 5 ×10⁻⁸ s. The incident wave is parallel-polarized, the incident angle is θ_i = 3π∕4, the maximum grid number in 𝜌 direction is i _max = 110, and the calculation time steps are 10000. Figure 3 shows the E _z field at the grid point (i = 55, k = k ₁, 𝜙 = 0). The results obtained by the two methods are in good agreement with the theoretical values. The time of the proposed algorithm is 8 s, while that of FFT is 409 s, which is much larger than the former. With the increase of computational space, the advantage of the proposed algorithm will be more and more significant.

5.2. Calculation of current distribution on the surface of metal cylinder

In this case, a metal cylinder is used, and the height and radius of the metal cylinder are both 60 m. The incident wave is applied pulse form, E _i(t) = E ₀ exp[−4π(t − t ₀)²∕τ²], which is parallel polarization, and the incident angle is θ_i = 3π∕4, τ = t ₀ = 1 × 10⁻⁷ s. In this section, the computation domain is discretized with a mesh Δ𝜌 × Δz = 1.5 m × 1.5 m, Δt = 0.417 ns, and M = 11.

The surface current density of the cylinder is decomposed into J _t (tangential component) and J _𝜙 (𝜙-direction component), and the following surface current density formulas can be obtained according to the boundary condition:

Bottom edge: $\quad J_{t}^{n}(i,k_{1})=h_{{\phi}}^{n}(i,k_{1}-1),\quad J_{{\phi}}^{n}(i,k_{1})=-h_{{\rho}}^{n}(i,k_{1}-1)$

Side: $\quad J_{t}^{n}(i_{1},k)=h_{{\phi}}^{n}(i_{1},k),\quad J_{{\phi}}^{n}(i_{1},k)=h_{z}(i_{1},k)$

Top edge: $\quad J_{t}^{n}(i,k_{2})=h_{{\phi}}^{n}(i,k_{2}-1),\quad J_{{\phi}}^{n}(i,k_{2})=-h_{{\rho}}^{n}(i,k_{2}-1)$ .

Figure 4 shows the time-domain waveform of current density at the middle point of the cylinder surface. J _t represents the tangential component of current density and J _𝜙 represents the 𝜙 directional component.

As shown in Fig. 5, the current density distribution of the cylindrical busbar with frequency of 10 MHz is extracted by Fourier transform and compared with the results of MoM, and the two results are also in good agreement.

Fig. 4.

Time-domain waveform of current density on the surface of the cylinder.

Fig. 5.

Current density distribution of cylindrical bus with 10 MHz.

Moreover, as shown in Table 1, it is obvious that the BOR-SMRTD method has an advantage over BOR-FDTD at saving the CPU time and memory. The proposed method uses 47.6% of the CPU memory and consumes 67.5% of the CPU time by the BOR-FDTD method.

Table 1

CPU time and memory for different schemes

Schemes	Mesh/Δ𝜌 ×Δz	CPU time/s	CPU memory/MB
		(20000Δt)
BOR-FDTD	1 m × 1 m	785.5	3.47
BOR-SMRTD	1.5 m × 1.5 m	530.2	1.65

Fig. 6.

RCS of metal sphere with radius 1 m.

Fig. 7.

RCS of metal sphere with radius 1 m.

5.3. RCS Calculation of metal ball

The radius of the metal sphere is 1 m, and the incident wave is still a Gauss pulse, τ = 3 × 10⁻⁹ s, t ₀ = 3 × 10⁻⁹ s. The computation domain is discretized with a mesh Δ𝜌 ×Δz = 1.5 m × 1.5 m, Δt = 6.67 ps, and M = 9. Figure 6 shows the far-field backscattering electric field of the metal ball. As shown in Fig. 7, Fourier transform is made for the electric field and incident pulse to obtain the RCS of the metal ball. Mie series solutions are also given as a comparison. It can be seen that the results obtained by BOR-FDTD are in good agreement with Mie series solutions, which also proves the accuracy of the proposed method.

6. Conclusion

The obliquely incident plane wave cannot be directly introduced into BOR-SMRTD due to its non-rotational symmetry. In this paper, an efficient one-dimensional time-domain algorithm for obliquely incident plane waves in BOR-SMRTD is proposed, which is based on the characteristic that the incident wave has only time delay in the symmetrical axis direction in the cylindrical coordinate system and the cylindrical expansion of plane waves. Numerical results show that the speed of the proposed algorithm is 51 times faster than that of the algorithm using FFT. In order to verify the effectiveness of the algorithm, the scattering problem of a metal cylinder with finite length is calculated. The surface current density distribution along the cylindrical generatrix at a given frequency is extracted, and the results are consistent with those obtained by MoM and BOR-FDTD methods. Furthermore, it is shown that the proposed algorithm uses only 47.6% of the CPU memory and consumes 67.5% of the CPU time by the BOR-FDTD method. Then the RCS of a metal spheres is calculated, and the numerical result is also in good agreement with the theoretical values.

Footnotes

Acknowledgements

This work was supported by the 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

D.Y.

Borges

B.H.V.

and Teixeira

F.L.

, Finite element time-domain body-of-revolution Maxwell solver based on discrete exterior calculus, Journal of Computational Physics 376 (2018), doi:10.1016/j.jcp.2018.09.024.

Gentili

G.G.

Khosronejad

Pelosi

and Stefano

, Analysis of elliptical structures with constant axial ratio by body-of-revolution finite element method and transformation optics, International Journal of Microwave and Wireless Technologies 11 (2019), 1–8, doi:10.1017/S1759078719000217.

Warecka

Lech

and Kowalczyk

, Efficient finite element analysis of axially symmetrical waveguides and waveguide discontinuities, IEEE Transactions on Microwave Theory and Techniques 67(11) (2019), 4291–4297, doi:10.1109/TMTT.2019.2940021.

Chen

Y.W.

J.Z.

Liu

Y.W.

Dong

X.H.

and Zhang

, The alternating direction implicit body of revolution multiresolution time domain method with convolution perfect matched layer, Mathematical Problems in Engineering 2017(2) (2017), 1–8, doi:10.1155/2017/2073563.

Hong

Huang

W.F.

Chen

and Liu

, Novel and stable formulations for the response of horizontal-coil eccentric antennas in a cylindrically multilayered medium, IEEE Transactions on Antennas and Propagation (2017), 1–1, doi:10.1109/TAP.2017.2670360.

Davidson

D.B.

and Ziolkowski

R.W.

, Body-of-revolution finite-difference tim-domain modeling of space-time focusing by a three-dimensional lens, J. Opt. Soc. Am. A 11(4) (1994), 1471–1490, doi:10.1364/JOSAA.11.001471.

Van

D.V.M.J.M.

and De Maagt

P.J.I.

, Efficient body of revolution finite-difference time-domain modeling of integrated lens antennas, IEEE Microwave and Wireless Components Letters 12(7) (2002), 258–260, doi:10.1109/LMWC.2002.801138.

Shi

S.Y.

and Prather

D.W.

, Electromagnetic analysis of axially symmetric diffractive optical elements illuminated by oblique incident plane waves, J. Opt. Soc. Am. A 18(1) (2001), 2901–2907, doi:10.1364/JOSAA.18.002901.

Chen

Y.C.

Mittra

and Harms

, Finite- difference time-domain algorithm for solving Maxwell’s equations in rotationally symmetric geometries, IEEE Trans. Microwave Theory Tech. 44(6) (1996), 832–839, doi:10.1109/22.506441.

10.

W.H.

Arakaki

and Mittra

, On the solution of a class of large body problems with full or partial circular symmetry by using the finite-difference time-domain (FDTD) method, IEEE Trans. Antennas Propag. 48(12) (2000), 1810–1817, doi:10.1109/8.901269.

11.

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.

12.

Sullivan

D.M.

, Three-dimensional computer simulation in deep regional hyperthermia using the finite-difference time-domain method, IEEE Trans. Microwave Theory and Techniques 38(2) (1990), 204–211, doi:10.1109/22.4643.

13.

Piket-May

M.J.

Taflove

Lin

W.C.

, Initial results for automated computational modeling of patient-specific electromagnetic hyperthermia, IEEE Trans. Biomedical Engineering 39(3) (1992), 226–237, doi:10.1109/10.125007.

14.

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.

15.

Hirono

Lui

Seki

and Yoshikuni

, A three-dimensional fourth-order finite-difference time-domain scheme using a symplectic integrator propagator, IEEE Transactions on Microwave Theory and Techniques 49(9) (2001), 1640–1648, doi:10.1109/22.942578.

16.

Prather

D.W.

and Shi

S.Y.

, Formulation and application of the finite-difference time-domain method for the analysis of axially symmetric diffractive optical elements, J. Opt. Soc. Am. A 16(5) (1999), 1131–1142, doi:10.1364/JOSAA.16.001131.