Abstract
Based on polarizability in the form of a complex quadratic rational function, a novel finite-difference time-domain (FDTD) approach combined with the Newmark algorithm is presented for dealing with a complex dispersive medium. In this paper, the time-stepping equation of the polarization vector is derived by applying simultaneously the Newmark algorithm to the two sides of a second-order time-domain differential equation obtained from the relation between the polarization vector and electric field intensity in the frequency domain by the inverse Fourier transform. Then, its accuracy and stability are discussed from the two aspects of theoretical analysis and numerical computation. It is observed that this method possesses the advantages of high accuracy, high stability, and a wide application scope and can thus be applied to the treatment of many complex dispersion models, including the complex conjugate pole residue model, critical point model, modified Lorentz model, and complex quadratic rational function.
1. Introduction
The finite-difference time-domain (FDTD) method has been widely used to study electromagnetic (EM) wave interaction with a wide variety of materials due to its robustness and its ability to calculate accurate broadband response via a single simulation [1, 2]. In the analysis of dispersive materials using the FDTD approach, it is required that the variation of dielectric parameters with frequency is modeled efficiently and accurately. In recent years, some novel dispersive models have been introduced, for example, the complex-conjugate pole-residue (CCPR) model [3], critical point (CP) model [4], modified Lorentz (m-Lo) model [5], and quadratic complex rational function (QCRF) model [6]. Some hybrid dispersion models, such as the Drude-Lorentz model [7, 8] and Drude-CP model [4, 9], were also usually adopted. However, the introduction of these new dispersion models and their hybrid models undoubtedly challenges the versatility of the existing FDTD methods. For the typical dispersive FDTD approaches, the recursive convolution (RC) method and Z transform (ZT) method are not convenient to deal with the high-order models such as m-Lo and QCRF. The auxiliary differential equation (ADE) method has more stringent stability conditions than the Courant-Friedrichs-Lewy (CFL) condition of Yee’s FDTD scheme [10–12].
The Newmark algorithm was originally used for the numerical solution of structural dynamics and was introduced into the finite-element time-domain (FETD) method for electromagnetic simulations [13–15]. Then, it is introduced into the analysis of the dispersive FDTD computation in 2014 for the unified treatment of the Debye, Drude, and Lorentz medium [16–18]. In this paper, the Newmark-FDTD method is further extended to the high-order model so that it can be used as a general approach for dispersive media, not only to treat with Debye, Drude, and Lorentz media, but also to deal with complex dispersion models including CCPR, CP, m-Lo, and QCRF media. The scope of application of the algorithm is extended greatly. Then, the computational accuracy and stability are investigated from the two aspects of theory and numerical computation, and higher accuracy and stability are shown.
2. Extended Newmark-FDTD Method
For the linear electric dispersive medium, the Maxwell’s curl equation can be written as where , , , and are the electric field strength, electric displacement, magnetic field strength, and magnetic flux density, respectively. is the electric conductivity. For simplicity, consider only the electric dispersion medium. The constitutive relation in a frequency domain can be characterized as where is the magnetic permeability. and are the relative electric permittivity in free space and at the infinite frequency, respectively. and are the electric permittivity and electric susceptibility in frequency domain, respectively. In this paper, the electric susceptibility is adopted in the following rational function form: where and (i = 0, 1, 2) are both real constants. In this way, the three typical kinds of dispersive models, that is, Debye, Drude, and Lorentz, and some of the complex dispersive models, that is, CCPR, CP, m-Lo, and QCRF, can be viewed as special cases of (5). The higher degrees of freedom and wider range of application are possessed.
The time-stepping equation of can be simply obtained by applying the Yee discrete scheme to (2). However, to obtain the time-stepping equation of from (1), we must firstly solve (4) in a time-stepping manner. For this reason, the polarization vector in the frequency domain is introduced and let ; thus, (4) can further be represented as Substituting (5) into , we can obtain Converting (8) into a time domain by the relation , we arrive at a two-order differential equation It can be observed that the right side of (9) not only contains but also contains the first-order derivative and second-order derivative of so that it cannot be solved directly using the Newmark algorithm employed in [16, 17]. For this reason, a temporary variable W is introduced and (9) is further represented as Then, the Newmark algorithm is applied to (10) with a unified time step and to obtain higher stability, and are adopted [13, 14]. Thus, we can obtain where , , and are similar to , , and . It only needs replacing , , and in (13), (14), and (15) with , , and .
Due to be a unified , in (11) and (12), the value of at , , and must be the same, so we can obtain the update formula Finally, substituting (16) and (7) into the discretized form of (1), the FDTD-updating formula of E is given In this way, a complete FDTD time-domain iterative computation is constituted by (17) and (16) and the standard FDTD formula of corresponding to (2). The FDTD computation can be realized.
3. Analysis of Accuracy and Stability
In the FDTD method, Maxwell’s curl equations are replaced with a different equation and the dispersion relation in the discrete form is thus quite different from that in the continuous form. For the dispersive FDTD method, the accuracy of the FDTD method is related also to a discrete scheme applied to the dispersive model of the media in the time domain. So, the accuracy of FDTD method can clearly be understood by investigating the errors between the numerical dispersive model of media and analytical model of media. The smaller the errors, the higher the computational accuracy [19]. The accuracy of the extended Newmark-FDTD algorithm is discussed below. Moreover, for simplicity and without loss of generality, let and ; the permittivity of dispersive media becomes a rational fraction form similar to (5).
To investigate the accuracy of the extended FDTD method, the z-domain electric susceptibility is given by applying the z transform to (16) Substituting , , , , , , and into (18) and then rearranging the resulting equation, we can arrive at To derive a more compact form, the following double-angle formula of trigonometric identity is applied to (19): The numerical electric susceptibility corresponding to (5) is derived as follows: where Comparing (21) with (5), it is shown that the forms of numerical electric susceptibility and electric susceptibility are exactly the same, apart from replacing with .
By expanding (22) into the power series, we can obtain The error between and is not higher than . For practical FDTD computation, restricted by CFL stable condition is always very small. Hence, is very close to . In comparison with the ADE method given in [10, 19], the extended Newmark-FDTD method possesses higher accuracy.
Then, the stability of the extended FDTD method is investigated, which is of great importance to fully employ this method for general dispersive media. The stability of the dispersive FDTD method depends not only on the dispersive model but also on the discretization scheme [11]. It can be carried out by the Von Neumann method combined with the Hurwitz-Routh (R-H) criterion [20]. Consider a time-dependent wave equation in a source-free homogeneous dispersive medium, we have where In (24) and (25), and are complex amplitudes of field quantities, is a numerical wavenumber in s direction, is the space step in s direction, and . Substituting (18) into (4), we have The stability polynomial in z-plane can be given as follows from (24) and (26): It is known from the knowledge of signal processing that the root of must lie in the unit circle on the z-plane to make the FDTD algorithm stable. In order to avoid solving directly the root of , we can transform (27) into the s-plane by the relation , then substituting , , , , , and , the following stability polynomial in the s-plane is obtained: In this way, the problem of judgment of whether the root of is in the unit circle on the z-plane is converted into the problem of whether the root of lies in the left of the s-plane, which can be calculated according to the R-H criterion with the Routh table. Finally, the following effective stability conditions are obtained: Taking (32) into account, for , the last numerical stability condition of the extended Newmark-FDTD method can be written in terms of : A comparison of the extended Newmark-FDTD algorithm with the double average scheme (DAS) of the ADE method as well as the Mobius transformation (MT) method indicates that a similar time-domain update formula is presented except for the coefficient difference with a constant in the three dispersive FDTD algorithms, although these three formulations are based on different approximation principles [11, 21–23]. The aforementioned numerical stability condition expressed by (33) is also the same as the stability condition of DAS and MT and is more relaxed than the stability condition of the direct scheme (DS) of the ADE algorithm [12]. Therefore, the extended Newmark-FDTD method also has the best stability [24, 25]. However, the derivation of the extended Newmark scheme is much simpler than that of the ADE and MT method, as the problem of the discrete scheme and time synchronism in the ADE algorithm as well as the tedious derivation in the MT algorithm can be avoided [22, 26].
4. Numerical Verification
To verify the advantages of extended Newmark-FDTD method, the numerical analysis of accuracy and stability of this method is further conducted.
Firstly, to verify the computational accuracy of the extended Newmark-FDTD algorithm, consider a 1-D human dry skin tissue modeled in the frequency range of 400 MHz–3 GHz by a QCRF with the coefficients of , , , , , , and [19]. The error is measured using the relative error between and defined as
The results at and are shown in Figure 1, where and are one tenths of the minimum wavelength in incident wave. As a comparison, the relative error of the direct scheme (DS) of the ADE method is also given. It is seen from Figure 1 that both the relative errors of the extended Newmark algorithm and ADE algorithm are very small and decrease with the decrease of the discrete interval . A comparison of the two methods shows that the extended Newmark algorithm has a smaller magnitude of relative error than the ADE method, so it has higher computational accuracy.
The stability of the extended Newmark-FDTD method is also investigated by the root locus of the stability polynomial in z-plane as (27). Consider a human muscle tissue modeled in the frequency range of 400 MHz–3 GHz by QCRF as (5) with the following parameters [21]: , , , , , , , and . Let m and , then the maximum time step satisfying the CFL stability limit of 1-D FDTD in free space is ps and the time step of 1-D FDTD in human muscle tissue is , where , is computed from (33). Then, the root locus of the stability polynomial of an extended Newmark-FDTD implementation is computed and shown in Figure 2(a) for . At the same time, the root locus of the stability polynomial of the direct scheme of the ADE algorithm is given in Figure 2(b).
(a) Extended Newmark-FDTD method
(b) Direct scheme of ADE method
It is seen from Figure 2(a) that for 0 < v < 5.1177, all the roots of the stability polynomial of the extended Newmark-FDTD implementation are inside the unit circle that shows the extended Newmark-FDTD method is stable for satisfied (33). This is consistent with the theoretical analysis. On the other hand, it is seen from Figure 2(b) that one root of the stability polynomial of the DS-ADE method escapes outside the unit circle for . This indicates that the direct scheme of the ADE method becomes unstable. This result is also consistent with the results from [21]. The comparison of Figures 2(a) and 2(b) shows that the stability condition of the extended Newmark scheme is more relaxed than that of the ADE direct implementation scheme.
To validate further the theoretical stability results, the FDTD simulation of electromagnetic wave propagation in a 1-D domain is considered. The simulation domain is composed of a human muscle tissue slab with the same parameters as described above and sandwiched between two semi-infinite dielectrics with . The size of the simulation domain is limited to 8000 space cells and terminated by the convolutional perfectly matched layer (CPML) with a thickness of 10 cells. The muscle slab is located between the nodes 4000 and 5000 of the domain. A modulated Gaussian pulse given by is introduced at the S point that is 10 away from the left boundary of the simulation domain (where GHz and ns). The electric field can be observed at point O that is 3000 away from the left boundary of the simulation domain [21]. The simulation domain geometry is shown in Figure 3.
If the spatial step is m, the maximum time step of the 1-D FDTD computation in free space that satisfied the Courant stability condition is ps and the upper bound of the time step of 1-D FDTD in the human muscle tissue for the discretization of extended Newmark FDTD is ps obtained from (33). The maximum time step of the 1-D FDTD computation in the surrounding semi-infinite dielectrics is ps. In the first test, the FDTD simulation is performed with a time step limit . This time step selection meets the stability limit in the human muscle tissue. Figure 4(a) shows the time-domain evolution of as recorded at the observation point O. It can be seen that the field remains stable during the whole period of time. The simulation is repeated by using ps that is bigger than the stability limit in the human muscle tissue. Figure 4(b) shows the time-domain evolution of as recorded at observation point O. It can be seen that the field begins to increase without being bound in the late-time response, and therefore, it is unstable.
(a)
(b)
5. Conclusion
In this paper, an extended Newmark-FDTD approach for the complex dispersive medium is presented in terms of polarizability in the form of a quadratic complex rational function, in which the time-stepping equation of the polarization vector is obtained by the Newmark algorithm used widely in finite-element time-domain (FETD) computation. The high accuracy and stability of this method are demonstrated from the two aspects of theoretical analyses and numerical results. It is also indicated that the application scope of this approach is extended greatly. In addition, this approach can be further used for the analysis of some complex dispersive models and some hybrid models, for example, CCPR, CP, m-Lo, QCRF, Drude-Lorentz, and Drude-CP.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This project is supported by the National Natural Scientific Foundation of China (Grant no. 61701428) and the Foundation of Construction of High-Level University Project in Shaanxi Province, China (Grant no. 2015SXTS02).