WU Peiyu ,YU Han ,HU Yenan,* ,XIE Yongjun ,and JIANG Haolin

1.School of Electronic and Information Engineering,Beihang University,Beijing 100191,China;2.Beijing Institute of Aerospace Systems Engineering,Beijing 100076,China;3.School of Electronic &Information Engineering,Nanjing University of Information Science and Technology,Nanjing 210096,China

Abstract:In order to simulate metamaterial rotational symmetric open region problems,unconditionally stable perfectly match layer (PML) implementation is proposed in the body of revolution (BOR) finite-difference time-domain (FDTD) lattice.More precisely,the proposed algorithm is implemented by the Crank-Nicolson (CN) Douglas-Gunn (DG) procedure for BOR metamaterial simulation.The constitutive relationship of metamaterial can be expressed by the Drude model and calculated by the piecewise linear recursive convolution (PLRC) approach.The effectiveness including absorption,efficiency,and accuracy is demonstrated through the numerical example.It can be concluded that the proposed implementation is to take the advantages of the CNDG-PML procedure,PLRC approach,and BORFDTD algorithm in terms of considerable accuracy,enhanced absorption and remarkable efficiency.Meanwhile,it can be demonstrated that the proposed scheme can maintain its unconditional stability when the time step exceeds the Courant-Friedrichs-Levy (CFL) condition.

Keywords:body of revolution (BOR),Crank-Nicolson (CN),finitedifference time-domain (FDTD),perfectly match layer (PML),metamaterial.

1.Introduction

With the characteristics of directly solving the Maxwell’s equations in time domain,the finite-difference time-domain(FDTD) algorithm which can obtain the responses over a broadband has become more prevalent than ever [1-3].The whole computational domain can be discretized according to the regulation of Yee’s grid.However,by applying the conventional FDTD algorithm to rotational symmetric structures,such cube grid generation affects the entire calculation accuracy especially for curve boundaries [4].The body of revolution (BOR) FDTD is proposed to alleviate such condition.In the BOR-FDTD algorithm,three-dimensional problems can be converted to two-dimensions resulting in the significant improvement of efficiency and accuracy [5].As a time-explicit algorithm,the stability is limited by the Courant-Friedrichs-Levy (CFL) condition which means the established relationship between mesh size and time step [6].By applying them directly to the large amount of time step problems,the simulation duration will become unacceptable.Unconditionally stable algorithms were proposed to remove the CFL condition.Among them,the Crank-Nicolson (CN) procedure which can solve Maxwell’s equation during a single iteration has aroused much attention [7].The original CN algorithm is efficient merely in one-dimension situation.By extending CN procedure to multi-dimensions,large sparse matrices are formed resulting in expensive computation.To alleviate such condition,approximate CN procedure is introduced including the approximate-decoupling (AD) and Douglas-Gunn(DG) schemes [8,9].It has been testified that the CNDG procedure is much more accurate through adding the distribution terms together at both sides of equations [10].

To simulate problems in open regions,absorbing boundary condition must be employed at the boundaries[6].The perfectly matched layer (PML) is the most prevalent one [11].As the original split-field PML shows quite low efficiency and absorption,the unsplit-field stretched coordinate PML and complex-frequency-shifted(CFS) PML are proposed to simplify the implementation,reduce late-time reflections and attenuate evanescent waves,respectively [12-15].Investigations on metamaterial employing CN procedure with double negative characteristic are mainly focused on the rectangular coordinate system [16-18].In addition,the constitutive relationship can be expressed by the Drude model for the double negative unique property [19,20].For the considerable calculation accuracy,the piecewise linear recursive convolution (PLRC) approach is employed for the simulation of the Drude model [21-23].

Here,by incorporating CNDG procedure,PLRC approach,unconditionally stable PML implementation is proposed for rotational symmetric metamaterial in the BOR-FDTD algorithm.The proposed implementation is to take advantages of the above-mentioned methods in terms of considerable absorption,remarkable accuracy and enhanced efficiency with the increment of time steps which exceeds the CFL condition.

2.Formulation

In the CFS-PML regions for TEφwave,the Maxwell’s equations can be given in the frequency domain as

whereDη,(η=r,φ,z) is the electric flux density andBηis the magnetic flux density along η direction,respectively.They can be obtained asDη(ω)=ε0εr(ω)Eη(ω) andBη(ω)=μ0μr(ω)Hη(ω).Sηis the stretched coordinate variable.In the metamaterial regions,constitutive relationship is expressed by Drude model.Meanwhile,εr(ω)and μr(ω) are assumed to be identical,given as

where ωpis the plasma frequency and ν is the damping constant.Within (1b),is the complex spatial coordinate stretching variable,given as

wherer1is interface between boundaries and PML interface alongr-direction,andSris the stretch coordinate variable alongr-direction.Sηis the stretched coordinate variable with CFS factor,defined as

By transforming Maxwell’s equations intoz-domain,one obtains

By introducing the auxiliary variables and rewriting the resultants,(7a)-(7c) can be reorganized as

whereFφz,Fzr,Pzr,Gφz,andGφrare the auxiliary variables.Due to the fact thatFφzandGφzhold almost similar forms,FzrandPzrare selected as examples for the further demonstration which can be given as

Within (8b) and (9b),coefficients λrand θrcan be given as

By substituting auxiliary variables (9a) and (9b) into(8a)-(8c) and employing thez-transform characteristic,the results can be discretized after employing the PLRC method and the CN scheme as

The operator δηis the first-order finite difference form,for example,

The auxiliary variables can be given as

It can be observed that the group equations of (10) are coupled.Huge sparse matrices must be calculated at each time step by directly employing them as an algorithm.To alleviate such condition,the approximate CN algorithm is introduced which is the CNDG procedure.According to the CNDG procedure,(10a) and (10b) are substituted into(10c),the results can be given as

whereAnand theBn-1represent the terms at thenth and the (n-1)th time step,respectively.It can be observed that magnetic component can be updated implicitly by solving tri-diagonal matrices.OnceHφis calculated,the other components and auxiliary variables can be calculated explicitly.

3.Numerical examples

The effectiveness of the proposed implementation including accuracy,efficiency,and absorption is demonstrated through the full-filled structure.Different PML implementations are performed on a PC with Intel Core i7-10700K 3.80 GHz and DDR4 (128 GB,3 200 MHz).The sketch picture of computational domain is shown in Fig.1.

Fig.1 Sketch of the computational domain of the metamaterial structure in the BOR algorithm

It can be demonstrated that the whole computational domain takes dimensions of 70Δr×80Δz.The structure can be revolved alongz-axis.Inside the computational domain,the perfectly electronic conductor (PEC) bulk with the size of 30Δr×30Δzis employed.Four vertexes are located at the positions of (28,33),(53,33),(28,58),and (53,58),respectively.The rest part of the structure is filled with metamaterial which has the parameters of ωp=rad/s and ν=1012rad/s,wheref0is the center frequency of source [24].The plane source incident from the positive side alongr-direction.The source locates at 20 cells distance in verticalz-direction.The modulated Gaussian pulse with the maximum frequency of 3 THz and bandwidth of 1 THz is employed.Three sides of the computational domain are terminated by a 10-cell-PML to absorb outgoing waves and reduce wave reflections.The receiver point is employed to receive propagation waves and evaluate wave reflection.The receiver point is located at the right bottom corner of the entire domain with the distances of one grid from both sides of the PML regions.Inside the PML regions,the parameters are selected to obtain the best absorbing performance both in the time domain and frequency domain.For comparison,FDTD-PML in [25,26] and ADI-PML in[27-29] are employed for the further demonstration.Meanwhile,the commercial software high frequency simulate software (HFSS) uniaxial PML (UPML) and CST convolutional PML (CPML) are also introduced for the comparison [30,31].The parameters of them are chosen as κη=140,αη=0.6,mη=3,and ση_max=1.9ση_opt,where

For the demonstration of accuracy and effectiveness,the uniform mesh size is employed as Δr=Δz=Δ=0.5 μm.The maximum time step of the conventional FDTD algorithm (=0.95Δ/(2c0)) can be obtained as 79.1 fs wherec0is the speed of light.The absorbing performance can be reflected by the relative reflection error in the time domain which can be defined as

where(t) is the test solution which can be calculated by Maxwell’s equations directly,(t) is the reference solution.The reference solution can be simulated by the extended domain with the size of 700Δr×800Δzand terminated by a 128-cell-PML without changing the relative position between the source and the receiver.Owing to such procedure,the reflection wave cannot be observed at the receiver during the calculation.Fig.2 shows the relative reflection error with different PML algorithms when CFLN=1 and 16,respectively.

Fig.2 Relative reflection error obtained by different PML algorithms and CFLNs in the time domain

For clarifying the demonstration,0.8 ns is employed in Fig.2.The absorbing performance can be evaluated by the maximum relative reflection error (MRRE) and latetime reflections.The MEER obtained by (i) FDTD-PML,(ii) HFSS-HPML,(iii) CST-CPML,(iv) ADI-PML with CFLN=1 and 16,(v) CNDS-PML with CFLN=1 and 16 are -86.5 dB,-60.8 dB,-63.7 dB,-69.8 dB,-73.9 dB,-46.2 dB,and -56.5 dB,respectively.The proposed scheme shows better performance compared with the commercial software including HFSS-UPML and CSTCPML.Such phenomenon indicates the effectiveness of the proposed implementation.Compared with the published ADI-PML,it can be observed that the proposed scheme can further decrease MRRE and late-time reflections indicating the improvement of the entire absorption with different CFLNs [28,29,32-34].

It can be observed that the relative reflection error of unconditionally stable algorithms degenerates slightly compared with FDTD-PML.This means the absorption decrement with the increment of CFLNs.The reason is that the numerical dispersion becomes larger in such circumstance.Furthermore,MRRE can be decreased by 10.3 dB in such condition.The effectiveness can be further demonstrated by employing the computational efficiency and resources.Table 1 shows the memory consumption,iteration steps,CFLNs,CPU time,and time reduction with different PML algorithms.

Table 1 Comparison of CPU time, iteration steps, memory, reduction, and MRRE of different PML algorithms

As can be illustrated in Table 1 that the commercial software occupies much more computational resource compared with the proposed scheme including the memory and calculation time.Moreover,the memory consumption of CNDG-PML increases slightly compared with ADI-PML.The reason is that disturbance terms must be calculated at both sides of equations.In addition,it can be observed that the CPU time of CNDS-and ADIPMLs increases significantly compared with FDTD-PML when CFLN=1.The reason is that tri-diagonal matrices must be calculated during each time step.The efficiency can be further improved by employing larger CFLNs to decrease the total iteration steps.It should be noticed that although the absorbing performance degenerates by employing larger CFLNs,the efficiency can be improved by 76.5 % and 78.3 %,respectively.

The absorbing performance can be further illustrated by employing the reflection coefficient in the frequency,given as

where operator FFT{·} denotes the Fourier transformation.Fig.3 shows the reflection coefficient obtained by different PML algorithms in the frequency domain.It can be observed that the proposed scheme can further improve the absorption within the entire frequency band compared with commercial software and ADI-PML.Meanwhile,compared with FDTD-PML,the proposed scheme shows acceptable effectiveness especially in efficiency.

Fig.3 Reflection coefficient obtained by different PML algorithms and CFLNs in the frequency domain

4.Conclusions

The PML formulation is a medium dependent scheme whose implementation should be modified according different materials.As the previous work cannot be employed in the metamaterial simulation with double negative property,based upon the BOR-FDTD algorithm and CNDG procedure,unconditionally stable PML implementation is proposed for rotational symmetric structure in metamaterial open region problems.The proposed formulation shows significant difference compared with the previous work.The proposed scheme takes advantages of them in terms of considerable accuracy,enhanced absorption and remarkable efficiency.A numerical example is carried out for the demonstration of absorption,efficiency and accuracy when the time step exceeds CFL condition.Through the results,the effectiveness and efficiency of the proposed CNDG-PML implementation are further indicated and illustrated.