Abstract

By incorporating the higher-order concept with the perfectly matched later (PML) scheme, unconditionally stable approximate Crank–Nicolson algorithm is proposed for plasma simulation in open region problems. More precisely, the proposed implementation is based on the CN Direct-Splitting (CNDS) procedure for the finite-difference time-domain (FDTD) unmagnetized plasma simulation. The unmagnetized plasma can be regarded as frequency-dependent media which can be calculated by the piecewise linear recursive convolution (PLRC) method. The proposed implementation shows the advantages of higher-order concept, CNDS procedure, and PLRC method in terms of improved absorbing performance, enhanced computational efficiency, and outstanding calculation accuracy. Numerical examples are introduced to indicate the effectiveness and efficiency. It can be concluded from results that the proposed scheme shows considerable efficiency, accuracy, absorption, and unconditional stability.

1. Introduction

With its outstanding performance in broadband simulations, the finite-difference time-domain (FDTD) algorithm, proposed by Yee, has received considerable attention in modelling microwave devices, calculating Maxwell’s equations and analysing frequency-dependent dispersive materials [1]. As the mesh size and time step are limited by Courant–Friedrichs–Levy (CFL) condition, the simulation duration will become unacceptable in solving large amount of time-step problems [2]. One of the most efficient methods which can significantly improve computational efficiency is to remove the CFL condition. In order to realize such aim, series of unconditionally stable algorithms are proposed including alternating-direction-implicit (ADI), locally one-dimensional, and split-steps algorithms [35]. However, it should be noticed that these algorithms are implemented by substep procedure resulting in decrement of efficiency and accuracy. Fortunately, the Crank–Nicolson (CN) scheme which can update the equation within a single step is proposed. However, it should be noticed that the original CN scheme is merely efficient and accurate in one-dimensions. By directly employing the original CN scheme to multidimension problems, large sparse matrices will be formed resulting in the computation much more expensive. Such condition leads to the unrealized computation in some practical problems. To not only avoid the calculation of huge matrices but also improve the entire performance of the simulation, the approximate CN algorithms are proposed, including the two-dimensional approximate-decoupling (AD) and Douglas–Gunn (DG) procedures [6, 7]. However, it has been testified that both CNAD and CNDG algorithms cannot be extended directly to three-dimensions [8]. Thus, series of approximate CN algorithms are introduced in three-dimensions, including cycle-sweep (CS), direct-splitting (DS), and approximate-factorization-splitting (AFS) procedures [9, 10]. Most importantly, it should be noticed that the CNCS algorithm has been proven to be conditionally stable both in theory and experiments [11]. Although the CNAFS algorithm shows considerable accuracy, nine tridiagonal matrices must be calculated at each time step resulting in quite low efficiency [9, 12]. Thus, the CNDS algorithm with considerable accuracy and efficiency is regarded as one of the most potential algorithms in the practical engineering [13].

By simulating the infinite extension in finite space, absorbing boundary condition must be employed to terminate the unbounded lattice [2]. The perfectly matched layer (PML), proposed by Berenger, is regarded as one of the most powerful absorbing boundary conditions [14]. The Berenger’s PML is implemented by the split-field scheme resulting in degeneration of absorption and efficiency. The unsplit-field stretched coordinate (SC) PML is carried out which can simplify the implementations at corners and edges of PML regions [15]. The SC-PML scheme shows disadvantages in absorbing low-frequency evanescent waves and reducing late-time reflections [16]. To alleviate such phenomenon, the complex-frequency-shifted (CFS) PML is carried out [17]. However, it has been testified that the SC-PML and CFS-PML are both inefficiency in the low-frequency band. The reason is that low-frequency propagation waves cannot be absorbed [18]. Thus, the higher-order concept is incorporated with various PML schemes not only to solve such condition but also to enhance the absorbing performance during the whole simulation [1923]. It can be found that the efficiency within higher-order PML regions will be affected by increment of auxiliary variables and manipulation operators. To alleviate such problems, the unconditionally stable higher-order PML is regarded as one of the most powerful ways. During the past recent decades, the unconditionally stable higher-order PML schemes were mainly developed in two-dimensions [2427]. Recently, the unconditionally stable higher-order PMLs based on the ADI and CNAFS procedures are proposed [2831].

The frequency-dependent dispersive materials play an important part during the development of astronavigation. Especially in the hypersonic flight, it has been testified that the blackout phenomenon will occur resulting in interruption of communication system [32]. During this period, plasma sheath is formed at the surface of the hypersonic flight [33]. The plasma sheath can be explained by magnetized plasma model and unmagnetized plasma model [34]. The magnetized plasma model is introduced mainly to investigate the microcosmic performance [35]. The unmagnetized plasma model has gained considerable attention in the investigation on macroscopical behaviour [36]. The investigation on unmagnetized plasma plays an important part in the series of blackout related problems. Several numerical methods have been employed into the simulation of frequency-dependent unmagnetized plasma [3739]. Among them, the piecewise linear recursive convolution (PLRC) method which assumes that the electric field has a piecewise linear functional dependence during each time step and then employs the recursive convolution to calculate time convolution and the dielectric susceptibility function is deemed as one of the most accurate ones [40, 41].

Here, unconditionally stable higher-order approximate CN-PML algorithm is proposed for the plasma simulation. More precisely, the proposed implementation is based on CNDS procedure, higher-order PML formulation, and PLRC method for the termination of unmagnetized plasma. Numerical examples are carried out to demonstrate the effectiveness and efficiency including the target chrematistics example and antenna performance demonstration. Through the resultants, it can be concluded that the proposed scheme takes advantages of higher-order PML, CNDS algorithm, and PLRC method in terms of improved absorption and considerable computational efficiency when the time step surpasses far beyond CFL limit. Meanwhile, it can be demonstrated that the forming of plasma sheath significantly affects its target chrematistics and antenna performance.

2. Formulation

In higher-order PML regions for terminating unmagnetized plasma, the modified Maxwell’s equations can be given aswherewhere is the complex relative permittivity which can be expressed aswhere is the resonance radian frequency and is the damping constant. Furthermore, the operator in (1) and (2)representswhere is the higher-order stretched coordinate variable which can be defined aswhere and are assumed to be positive real and is assumed to be larger than one. By employing the partial fraction expansion, can be expressed aswhere the coefficients can be given as , , and . For simplifying the demonstration, and are selected as examples. By rewriting (1) and (2), the equations can be rewritten as

By substituting (7) into (8) and (9) and rearranging the resultants, the equations can be obtained according to the following forms as

Introducing the auxiliary variables and for updating, the equations can be simplified aswhere auxiliary variables can be given aswhere is another component of which can be further demonstrated by an example, i.e., when calculating components, , it can be observed that according to the relationship. Due to the fact that and hold the similar forms, is chosen as an example for demonstration. By substituting (13) and (14) into (11) and (12), the equations can be rewritten as

By employing the transformation relationship , PLRC method, and CN scheme, equations (16) and (17) can be rewritten in the FDTD domain aswhere the coefficients can be given as follows

The auxiliary variables can be given in the FDTD domain as the following forms:where the coefficients can be given as

The other components can be obtained by employing the same approach. It can be observed that large sparse matrices are formed by employing (17)–(19) directly as an algorithm resulting in the computation much more expensive which becomes unpractical sometimes. To avoid the calculation of such huge matrices and improve the calculation, the CNDS procedure is employed during the implementation. According to (17) and (18), the components can be written in matrix form aswhere is the identity matrix, , , and is the other components at the n-th time step:

According to the directly splitting procedure, equation (24) can be spitted into several sub-matrices as

Equations (25) and (26) can be written in the single field component as

By substituting (28) into (27) to eliminate , the equation can be given as

By substituting (30) into (29) to eliminate , the equation can be rewritten in the following forms:

It can be observed from (31) and (32) that tridiagonal matrices are formed. Thus, they can be updated by employing the Thomas algorithm. It can be concluded that large number of operators are introduced during the calculation resulting in quite low efficiency. Thus, the proper auxiliary variables can be introduced during the calculation. The method is the same as that in [13].

3. Numerical Results and Discussion

To demonstrate the effectiveness of proposed implementation, numerical examples are carried out. Firstly, the target characteristics of the hypersonic flight vehicle model are investigated to verify the accuracy of calculation. Then, the performance of the antenna model which is located on surface of flight is investigated. For comparison, the CFS-PML based on conventional FDTD algorithm in [42] (FDTD-PML), CFS-PML based on CNDS algorithm in [13] (CNDS-PML), CFS-PML based on CNAFS algorithm in [12] (CNAFS-PML), and HO-PML based on conventional FDTD algorithm in [19] are chosen as examples. For simplify the demonstration, the proposed implementation can be denoted as CNDS-HO-PML.

3.1. Target Characteristic of Flight Vehicle

The target characteristics of flight vehicle model is investigated not only to verify accuracy of the calculation but also to study the target characteristic of the formed plasma sheath on its surface. The sketch of flight vehicle model and its detailed parameters is shown in Figure 1. As shown in Figure 1(a), it can be regarded as the combination of half-sphere and pyramis models. The detailed parameters are shown in Figure 1(b). It can be observed that the radius of sphere is 0.2 m. In the pyramis model, the radiuses of upper and bottom are 0.2 m and 0.4 m, respectively. The material of flight vehicle model can be expressed by perfectly electronic conductor (PEC). The plasma sheath with the thickness of 0.1 m is formed on the surface of the model. It holds the parameters of Grad/s and .

The target characteristics are implemented by employing different PML algorithms. Through the obtained resultants, the accuracy of various algorithms is demonstrated. Figure 2 shows the sketch picture of FDTD computational domain with dimensions of . The model is located at the centre of domain. The rest part is filled with vacuum.

The incidence wave which is a 0.8 GHz Gaussian pulse propagates along positive side of y-direction. At the boundaries of domain, they are terminated by 10 cells PML region to absorb outgoing waves and reduce wave reflections. The parameters of PML regions are chosen to obtain the best absorbing performance both in time domain and frequency domain. The parameters of higher-order PMLs are chosen as , , , , , , , and , where

The parameters of CFS-PML are chosen as , , , and . The uniform mesh sizes are chosen as . The maximum time step of the conventional FDTD method can be obtained as . The CFL number (CFLN) is defined as , where is the time step of the unconditionally stable algorithm. The target characteristic can be reflected by the radar cross section (RCS) parameters. Figure 3 shows the RCS of the flight vehicle model versus frequency obtained by different PML algorithms and CFLNs.

It can be observed that the RCS parameters obtained by different PML algorithms are almost overlapped with CFLN = 1 indicating that they almost hold the similar accuracy. When CFLN = 10, it can be observed that CNDS-HO-PML holds the best simulation accuracy which is better than CNDS- and CNAFS-PML implementations. Figure 4 shows RCS of flight vehicle model with plasma sheath versus frequency obtained by different PML algorithms and CFLNs. One can obtain the same conclusion as previous that the CNDS-HO-PML holds the best calculation accuracy among unconditionally stable algorithms.

In addition, it can be observed from Figure 5 that the RCS with the plasma sheath changes significantly especially after the frequency of 0.4 GHz. The memory consumption, CPU time, CFLN, and iteration steps obtained by different PML algorithms are shown in Table 1.

Here, the algorithm is implemented by 32768 and 3277 steps when CFLN = 1 and 10, respectively. As shown in Table 1, it can be observed that the CPU time of unconditionally stable algorithm increases significantly compared with FDTD-PML. The reason is that the unconditionally stable algorithms solve tridiagonal matrices at each time step resulting in such conditions. It can be concluded that the CPU time can be decreased by employing larger CFLNs. The reason is that the total iteration steps can be decreased by employing larger time step. The proposed scheme with CFLN = 10 can save 41.1% compared with FDTD-PML which is close to CNAFS-PML.

3.2. Antenna Performance and Treatment in Plasma Sheath Open Region Problem

For the further demonstration of absorption, efficiency, and accuracy, the antenna performance inside plasma open region problem is simulated. Figure 6 shows the antenna model and its detail parameters on the surface of flight vehicle with plasma sheath. More precisely, antenna is located inside plasma sheath with parameters of and . Figure 6(a) demonstrates the location of antenna on the surface of flight. Such conditions can be simplified as a more intuitionistic model for simulation. The sketch picture of computational domain is shown in Figure 6(b). It is shown that the structure is composed of PEC plate and monopole antenna model. The whole computational domain takes the direction of . The PEC plate which has dimensions of is located at one-third height of z-direction. The monopole PEC antenna with the height of 1.75 mm is situated at the centre of xOy plane. The source which is a modulated Gaussian pulse with the centre frequency and maximum frequency of 4 GHz and 5 GHz is employed at the bottom of antenna for excitation. The rest part is filled with unmagnetized plasma. At the boundaries of domain, PML regions with 10 cells are employed to terminate the unbounded lattice and absorb outgoing waves.

Inside PML regions, the parameters are chosen to obtain the best absorbing performance both in the time domain and in the frequency domain. The parameters of higher-order PMLs are chosen as , , , , , , , and . For comparison, the parameters of the other PML algorithms are chosen as , , , and .

In order to demonstrate the effectiveness of the algorithm, uniform mesh sizes are chosen as . The time step is chosen as 0.58 ps according to the CFL condition. Thus, the whole computational domain can be discretized as . The absorption of the PML regions can be reflected by the relative reflection error in the time domain which can be defined aswhere is the test solution which can be measured directly from the observation point and is the reference solution which can be obtained by enlarging the computational domain by 20 times and terminating by 64-cell PML. During the reference solution calculation, the reflection wave can be ignored due to the thick PML regions and enlargement computational domain. Figure 7 shows the relative reflection error in the time domain obtained by different PML algorithms with CFLN = 1 and 16. The absorption can be reflected by the maximum relative reflection error (MRRE) and late-time reflections. As is shown in Figure 7(a), both MRRE and late-time reflections can be decreased by employing the higher-order formulation. The HO-PML based on the conventional FDTD algorithm holds the best absorption. The proposed scheme shows its inferior performance. Furthermore, compared with FDTD-PML, CNDS-PML, and CNAFS-PML, the proposed scheme shows the improved absorbing performance. Such conditions become much more significant with larger CFLNs. As shown in Figure 7(b), the entire absorbing performance becomes worse. The reason is that the numerical dispersion increases with the enlargement of CFLNs resulting in decrement of the computational accuracy. However, the absorption can be improved significantly with the proposed scheme. In conclusion, the absorption can be further enhanced by employing the higher-order PML formulation. The effectiveness of the proposed can not only be reflected by the absorbing performance but also the efficiency. Table 2 shows memory consumption, CPU time, CFLN, and iteration steps obtained by different PML algorithms in the antenna simulation.

It can be observed from Table 1 that the unconditionally stable algorithms occupy much computational resources and simulation duration. The reason is that several tridiagonal matrices are formed at left sides of the equations which can be solved by the Thomas decomposition resulting in such phenomenon. The CNAFS algorithm occupies much resources and duration compared with CNDS-PML. The reason is that nine matrices should be calculated in the CNAFS algorithm while the CNDS algorithm calculates six. The higher-order PML occupies much computational resources and simulation duration. The reason is that auxiliary variables and coefficients increases with the increment of PML order. The simulation duration can be decreased by employing larger CFLNs resulting in the decrement of iteration steps. Most importantly, the proposed scheme with CFLN = 16 shows much considerable absorption and efficiency compared with FDTD-PML. Such a condition indicates the effectiveness of the algorithm. The return loss (S11) is one of the most important parameters in antenna simulation. Thus, the absorption and accuracy can also be evaluated by the S11 in the frequency domain. Figure 8 shows S11 parameters obtained by different PML algorithms and CFLNs in the plasma environment. It can be shown in Figure 8(a) that the curves are almost overlapped indicating that they almost hold the same accuracy. It shown in Figure 8(b) that the proposed scheme shows the least shifting indicating that the proposed scheme shows the best performance. Compared with the CNDS-PML scheme, the CNAFS-PML scheme shows better accuracy.

In order to investigate the influence of plasma environment on its S11 parameter, such simulation is performed in the vacuum environment again. Figure 9 shows the S11 parameters in vacuum and plasma environments obtained by the proposed scheme when CFLN = 16. It can be observed that S11 parameters shift and decrease in the plasma environment compared with those inside vacuum. The shifting and the decrement of return loss result in the communication interruption in the blackout.

4. Conclusions

By incorporating the approximate CN procedure, higher-order formulation, and PLRC method, unconditionally stable PML implementation is proposed for plasma simulation. More precisely, the higher-order CNDS-PML implementation is proposed for the termination of unmagnetized plasma. The proposed scheme is to take advantages of these approaches in terms of improved absorbing performance, enhanced computational efficiency, and outstanding calculation accuracy. Numerical examples are carried out for further demonstration. The results indicate that the target characteristics of flight vehicle and the antenna performance inside the plasma sheath changes significantly compared with that without the sheath. Meanwhile, the results also demonstrate that the proposed implementation can receive considerable effectiveness, efficiency, and accuracy during the calculation.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (61571022 and 61971022) and National Key Laboratory Foundation (HTKJ2019KL504013 and 61424020305).