# Approximate Crank–Nicolson Algorithm with Higher-Order PML Implementation for Plasma Simulation in Open Region Problems

Authors: Liqiang Niu; Yongjun Xie; Jie Gao; Peiyu Wu; Haolin Jiang
Journal: International Journal of Antennas and Propagation
 (2021)

---

## 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.


---

## Body

## 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 [3–5]. 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 [19–23]. 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 [24–27]. Recently, the unconditionally stable higher-order PMLs based on the ADI and CNAFS procedures are proposed [28–31].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 [37–39]. 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 as(1)jωDω=∇s×Hω,(2)−jωμ0Hω=∇s×Eω,where(3)D=ε0εrωE,where εrω is the complex relative permittivity which can be expressed as(4)εrω=1+ωp2−ω2+jωυ,where ωp is the resonance radian frequency and υ is the damping constant. Furthermore, the operator ∇s in (1) and (2)represents(5)∇s=x^1Sx∂∂x+y^1Sy∂∂y+z^1Sz∂∂z,where Sη,η=x,y,z is the higher-order stretched coordinate variable which can be defined as(6)Sη=κη1+ση1αη1+jωε0κη2+ση2αη2+jωε0,where αηn,n=1,2 and σηn are assumed to be positive real and κηn is assumed to be larger than one. By employing the partial fraction expansion, Sη−1 can be expressed as(7)Sη−1=kηjω+aη1jω+bη1jω+bη2jω+bη2,where the coefficients can be given as kη=1/κη1κη2, aηn=αηn/ε0, and bηn=aηn+σηn/κηnε0. For simplifying the demonstration, Dz and Hy are selected as examples. By rewriting (1) and (2), the equations can be rewritten as(8)jωDz=1Sx∂Hy∂x−1Sy∂Hx∂y,(9)−jωμ0Hy=1Sz∂Ex∂z−1Sx∂Ez∂x,By substituting (7) into (8) and (9) and rearranging the resultants, the equations can be obtained according to the following forms as(10)jωDz=kxjω+ax1jω+bx1jω+bx2jω+bx2∂Hy∂x−kyjω+ay1jω+by1jω+by2jω+by2∂Hx∂y,jωμ0Hy=kxjω+ax1jω+bx1jω+bx2jω+bx2∂Ez∂xkzjω+az1jω+bz1jω+bz2jω+bz2∂Ex∂z.Introducing the auxiliary variablesFzηn and Gyηn for updating, the equations can be simplified as(11)jωDz=jωFzx2+ax1Fzx2−jωFzy2−ay1Fzy2,(12)jωμ0Hy=jωGyz2+az1Gyz2−jωGyx2−ax1Gyx2,where auxiliary variables can be given as(13)Fzη1=kzη1jω+bη1∂Hη˜∂η⟹jωFzη1+bη1Fzη1=kη∂Hη˜∂η,(14)Fzη2=jω+aη2jω+bη2Fzη1⟹jωFzη2+bη2Fzη2=jωFzη1+aη2Fzη1,where η˜ is another component of η which can be further demonstrated by an example, i.e., when calculating Dz components, η=x, it can be observed that η˜=y according to the relationship. Due to the fact that Fzηn and Gyηn hold the similar forms, Fzηn is chosen as an example for demonstration. By substituting (13) and (14) into (11) and (12), the equations can be rewritten as(15)jωDz=ax1−bx2Fzx2+ax2−bx1Fzx1−ay1−by2Fzy2−ay2−by1Fzy1+kx∂Hy∂x−ky∂Hx∂y,(16)jωμ0Hy=az1−bz2Gyz2+az2−bz1Gyz1−ax1−bx2Gyx2−ax2−bx1Gyx1+kz∂Ex∂z−kx∂Ez∂x.By employing the transformation relationshipjω⟷∂/∂t, PLRC method, and CN scheme, equations (16) and (17) can be rewritten in the FDTD domain as(17)Ezn+1=a1Ezn+a2φzn+p1xFzx2n+p2xFzx1n−p1yFzy2n−p2yFzy1n+p3xδxHyn+1+Hyn−p3yδyHxn+1+Hxn,(18)Hyn+1=Hyn+p4zGyz2n+p5zGyz1n−p4xGyx2n−p5xGyx1n+p6zδzExn+1+Exn−p6xδxEzn+1+Ezn,(19)φzn+1=b1Ezn+1+Δξ0Ezn+b2φxn,where the coefficients can be given as follows(20)p1η=a2Δtaη1−bη2ε0,p2η=a2Δtaη2−bη1ε0,p3η=a2Δtkη2ε0,p4η=Δtaη1−bη2μ0,p5η=Δtaη2−bη1μ0,p6η=Δtkη2μ0,a1=a21−ξ0,a2=11+χ0−ξ0,χ0=ωp2Δtv−ωp21−e−vΔtv2,Δχ0=−ωp21−e−vΔt2v2,ξ0=ωp2Δt2v−ωp21−1+vΔte−vΔtv3Δt,Δξ0=−ωp21−1+vΔte−vΔt,1−e−vΔtv3Δt.The auxiliary variables can be given in the FDTD domain as the following forms:(21)Fzη1n+1=p7ηFzη1n+p8ηδηHη˜n+1+Hη˜n,Fzη2n+1=p9ηFzη2n+p10ηFzη1n+1−p11ηFzη1n,where the coefficients can be given as(22)p7η=1−bη1Δt/21+bη1Δt/2,p8η=Δtkη2Δη1+Δtbη1/2,p9η=1−bη2Δt/21+bη2Δt/2,p10η=1+aη2Δt/21+bη2Δt/2,p11η=1−aη2Δt/21+bη2Δt/2.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 as(23)I−D1−D2Φn+1=I1+D1+D2Φn+An,where I is the identity matrix, I1=a1I, Φ=Ex,Ey,Ez,Hx,Hy,HzT, and A is the other components at the n-th time step:(24)D1=0000−p3zδz000000−p3xδx000−p3yδy0000−p6yδy000−p6zδz000000−p6xδx0000,D2=00000p3yδy000p3zδz000000p3xδx00p6zδz000000p6xδx000p6yδy00000.According to the directly splitting procedure, equation (24) can be spitted into several sub-matrices as(25)I−D1Φ∗=I1+D1+2D2Φn+An,(26)I−D2Φn+1=Φ∗−D2Φn.Equations (25) and (26) can be written in the single field component as(27)Ez∗=a1Ezn+a2φzn−p3y∂yHx∗+Hxn+2p3x∂xHyn+p1xFzx2n+p2xFzx1n−p1yFzy2n−p2yFzy1n,(28)Hy∗=Hyn−p6z∂zEx∗+Exn+2p6x∂xEzn+p4xGyx2n+p5xGyx1n−p4zGyz2n−p5zGyz1n,(29)Ezn+1=Ez∗+p3x∂xHyn+1+Hyn,(30)Hyn+1=Hy∗+p6x∂xEzn+1+Ezn.By substituting (28) into (27) to eliminate Hx∗, the equation can be given as(31)1−p6yp3y∂y∂yEz∗=a1+p6yp3y∂y∂yEzn+a2φzn−2p3y∂yHxn+2p3x∂xHyn−2p6zp3y∂y∂zEyn−p4zp3y∂yGxz2n−p5zp3y∂yGxz1n+p4yp3y∂yGxy2n+p5yp3y∂yGxy1n+p1xFzx2n+p2xFzx1n−p1yFzy2n−p2yFzy1n.By substituting (30) into (29) to eliminate Hyn+1, the equation can be rewritten in the following forms:(32)1−p6xp3x∂x∂xEzn+1=Ez∗−p3hzp3x∂x∂zEx∗+Exn+2p3x∂xHyn+p6xp3x∂x∂xEzn+p4xp3x∂xGyx2n+p5xp3x∂xGyx1n−p4zp3x∂xGyz2n−p5zp3x∂xGyz1n.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 Figure1. 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 ωp=2π×28.7 Grad/s and ν=20Grad/s.Figure 1
The sketch picture and detail parameters of flight vehicle model: (a) whole model; (b) front view.
(a)(b)The target characteristics are implemented by employing different PML algorithms. Through the obtained resultants, the accuracy of various algorithms is demonstrated. Figure2 shows the sketch picture of FDTD computational domain with dimensions of 400Δx×400Δy×500Δz. The model is located at the centre of domain. The rest part is filled with vacuum.Figure 2
The sketch picture of the computational domain with the flight vehicle model.The incidence wave which is a 0.8 GHz Gaussian pulse propagates along positive side ofy-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 κη1=90, αη1=2.1, mη1=2, ση1_max=0.4ση1_opt, κη2=1, αη2=0.2, mη2=3, and ση2_max=0.06ση2_opt, where(33)σηn_opt=mηn+1150πΔη.The parameters of CFS-PML are chosen asκη=120, αη=1.2, mη=3, and ση_max=2.3ση_opt. The uniform mesh sizes are chosen as Δx=Δy=Δz=Δ=3mm. The maximum time step of the conventional FDTD method can be obtained as ΔtmaxFDTD=5.78ps. The CFL number (CFLN) is defined as CFLN=Δt/ΔtmaxFDTD, where Δt 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.Figure 3
The RCS of the flight vehicle model versus frequency obtained by (a) FDTD-PML, CNDS-PML, CNAFS-PML, HO-PML, and CNDS-HO-PML CFLN = 1; (b) CNDS-PML, CNAFS-PML, and CNDS-HO-PML CFLN = 10.
(a)(b)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. Figure4 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.Figure 4
The RCS of the flight vehicle model versus frequency with plasma sheath obtained by (a) FDTD-PML, CNDS-PML, CNAFS-PML, HO-PML, and CNDS-HO-PML CFLN = 1; (b) CNDS-PML, CNAFS-PML, and CNDS-HO-PML CFLN = 10.
(a)(b)In addition, it can be observed from Figure5 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.Figure 5
Comparison of the RCS versus frequency between flight vehicle with plasma sheath and flight vehicle without plasma sheath obtained by CNDS-HO-PML with CFLN = 10.Table 1
The memory consumption, CFLN, CPU time, iteration steps, and time reduction obtained by different PML algorithms.
PML algorithmsCFLNStepsMemory (GB)CPU time (min)Reduction (%)FDTD-PML13276824.139.2—CNDS-PML13276834.8180.33599.5CNAFS-PML13276839.6360.68199.0HO-PML13276832.989.21275.5Proposed13276851.5427.89913.3CNDS-PML10327734.812.967.1CNAFS-PML10327739.620.747.2Proposed10327721.523.141.1Here, the algorithm is implemented by 32768 and 3277 steps when CFLN = 1 and 10, respectively. As shown in Table1, 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. Figure6 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 ωp=2π×28.7Grad/s and ν=20Grad/s. 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 20×20×10mm. The PEC plate which has dimensions of 10×10×1mm 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.Figure 6
The antenna location and the sketch picture of the computational domain: (a) the monopole antenna on the surface of the flight; (b) the computational domain.
(a)(b)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κη1=11, αη1=0.5, mη1=4, ση1_max=2.7ση1_opt, κη2=3, αη2=1.1, mη2=1, and ση2_max=0.9ση2_opt. For comparison, the parameters of the other PML algorithms are chosen as κη=22, αη=2.1, mη=1, and ση_max=0.3ση_opt.In order to demonstrate the effectiveness of the algorithm, uniform mesh sizes are chosen asΔx=Δy=Δz=Δ=0.3mm. The time step ΔtmaxFDTD is chosen as 0.58 ps according to the CFL condition. Thus, the whole computational domain can be discretized as 67Δx×67Δy×34Δz. The absorption of the PML regions can be reflected by the relative reflection error in the time domain which can be defined as(34)RdBt=20log10Eztt−EzrtmaxEzrt,where Eztt is the test solution which can be measured directly from the observation point and Ezrt 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.Figure 7
The relative reflection error versus time obtained by (a) FDTD-PML, CNDS-PML, CNAFS-PML, HO-PML, and CNDS-HO-PML when CFLN = 1; (b) CNDS-PML, CNAFS-PML, and CNDS-HO-PML when CFLN = 16.
(a)(b)Table 2
The memory consumption, CFLN, CPU time, iteration steps, and time reduction obtained by different PML algorithms.
PML algorithmsCFLNStepsMemory (GB)CPU time (min)Reduction (%)FDTD-PML1327681.14.7‐CNDS-PML1327681.852.9−1025.5CNAFS-PML1327682.575.1−1497.9HO-PML1327681.922.3−377.5Proposed1327683.381.6−1636.2CNDS-PML1632771.82.253.2CNAFS-PML1632772.52.644.7Proposed1632773.33.721.3It can be observed from Table1 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.Figure 8
TheS11 versus frequency obtained by (a) FDTD-PML, CNDS-PML, CNAFS-PML, HO-PML, and CNDS-HO-PML when CFLN = 1 (b) CNDS-PML, CNAFS-PML, and CNDS-HO-PML when CFLN = 16.
(a)(b)In order to investigate the influence of plasma environment on itsS11 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.Figure 9
TheS11 versus frequency obtained by CNDS-HO-PML when CFLN = 16 in plasma and vacuum, respectively.


## 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 Figure1. 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 ωp=2π×28.7 Grad/s and ν=20Grad/s.Figure 1
The sketch picture and detail parameters of flight vehicle model: (a) whole model; (b) front view.
(a)(b)The target characteristics are implemented by employing different PML algorithms. Through the obtained resultants, the accuracy of various algorithms is demonstrated. Figure2 shows the sketch picture of FDTD computational domain with dimensions of 400Δx×400Δy×500Δz. The model is located at the centre of domain. The rest part is filled with vacuum.Figure 2
The sketch picture of the computational domain with the flight vehicle model.The incidence wave which is a 0.8 GHz Gaussian pulse propagates along positive side ofy-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 κη1=90, αη1=2.1, mη1=2, ση1_max=0.4ση1_opt, κη2=1, αη2=0.2, mη2=3, and ση2_max=0.06ση2_opt, where(33)σηn_opt=mηn+1150πΔη.The parameters of CFS-PML are chosen asκη=120, αη=1.2, mη=3, and ση_max=2.3ση_opt. The uniform mesh sizes are chosen as Δx=Δy=Δz=Δ=3mm. The maximum time step of the conventional FDTD method can be obtained as ΔtmaxFDTD=5.78ps. The CFL number (CFLN) is defined as CFLN=Δt/ΔtmaxFDTD, where Δt 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.Figure 3
The RCS of the flight vehicle model versus frequency obtained by (a) FDTD-PML, CNDS-PML, CNAFS-PML, HO-PML, and CNDS-HO-PML CFLN = 1; (b) CNDS-PML, CNAFS-PML, and CNDS-HO-PML CFLN = 10.
(a)(b)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. Figure4 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.Figure 4
The RCS of the flight vehicle model versus frequency with plasma sheath obtained by (a) FDTD-PML, CNDS-PML, CNAFS-PML, HO-PML, and CNDS-HO-PML CFLN = 1; (b) CNDS-PML, CNAFS-PML, and CNDS-HO-PML CFLN = 10.
(a)(b)In addition, it can be observed from Figure5 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.Figure 5
Comparison of the RCS versus frequency between flight vehicle with plasma sheath and flight vehicle without plasma sheath obtained by CNDS-HO-PML with CFLN = 10.Table 1
The memory consumption, CFLN, CPU time, iteration steps, and time reduction obtained by different PML algorithms.
PML algorithmsCFLNStepsMemory (GB)CPU time (min)Reduction (%)FDTD-PML13276824.139.2—CNDS-PML13276834.8180.33599.5CNAFS-PML13276839.6360.68199.0HO-PML13276832.989.21275.5Proposed13276851.5427.89913.3CNDS-PML10327734.812.967.1CNAFS-PML10327739.620.747.2Proposed10327721.523.141.1Here, the algorithm is implemented by 32768 and 3277 steps when CFLN = 1 and 10, respectively. As shown in Table1, 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. Figure6 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 ωp=2π×28.7Grad/s and ν=20Grad/s. 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 20×20×10mm. The PEC plate which has dimensions of 10×10×1mm 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.Figure 6
The antenna location and the sketch picture of the computational domain: (a) the monopole antenna on the surface of the flight; (b) the computational domain.
(a)(b)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κη1=11, αη1=0.5, mη1=4, ση1_max=2.7ση1_opt, κη2=3, αη2=1.1, mη2=1, and ση2_max=0.9ση2_opt. For comparison, the parameters of the other PML algorithms are chosen as κη=22, αη=2.1, mη=1, and ση_max=0.3ση_opt.In order to demonstrate the effectiveness of the algorithm, uniform mesh sizes are chosen asΔx=Δy=Δz=Δ=0.3mm. The time step ΔtmaxFDTD is chosen as 0.58 ps according to the CFL condition. Thus, the whole computational domain can be discretized as 67Δx×67Δy×34Δz. The absorption of the PML regions can be reflected by the relative reflection error in the time domain which can be defined as(34)RdBt=20log10Eztt−EzrtmaxEzrt,where Eztt is the test solution which can be measured directly from the observation point and Ezrt 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.Figure 7
The relative reflection error versus time obtained by (a) FDTD-PML, CNDS-PML, CNAFS-PML, HO-PML, and CNDS-HO-PML when CFLN = 1; (b) CNDS-PML, CNAFS-PML, and CNDS-HO-PML when CFLN = 16.
(a)(b)Table 2
The memory consumption, CFLN, CPU time, iteration steps, and time reduction obtained by different PML algorithms.
PML algorithmsCFLNStepsMemory (GB)CPU time (min)Reduction (%)FDTD-PML1327681.14.7‐CNDS-PML1327681.852.9−1025.5CNAFS-PML1327682.575.1−1497.9HO-PML1327681.922.3−377.5Proposed1327683.381.6−1636.2CNDS-PML1632771.82.253.2CNAFS-PML1632772.52.644.7Proposed1632773.33.721.3It can be observed from Table1 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.Figure 8
TheS11 versus frequency obtained by (a) FDTD-PML, CNDS-PML, CNAFS-PML, HO-PML, and CNDS-HO-PML when CFLN = 1 (b) CNDS-PML, CNAFS-PML, and CNDS-HO-PML when CFLN = 16.
(a)(b)In order to investigate the influence of plasma environment on itsS11 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.Figure 9
TheS11 versus frequency obtained by CNDS-HO-PML when CFLN = 16 in plasma and vacuum, respectively.


## 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.



---

*Source: 6618492-2021-02-10.xml*