J Infrared Milli Terahz Waves (2010) 31:620–628 DOI 10.1007/s10762-010-9619-y
Finite-difference Time-domain Algorithm for Plasma Based on Trapezoidal Recursive Convolution Technique Song Liu & Sanqiu Liu & Shaobin Liu
Received: 29 November 2009 / Accepted: 21 January 2010 / Published online: 9 February 2010 # Springer Science+Business Media, LLC 2010
Abstract The electromagnetic wave propagation in plasma media is modeled using finitedifference time-domain (FDTD) method based on the trapezoidal recursive convolution (TRC) Technique. The TRC Technique requires single convolution integral in the formulation as in the recursive convolution (RC) method, while maintaining the accuracy comparable to the piecewise linear convolution integral (PLRC) method with two convolution integrals. The three dimensional (3-D) TRC-FDTD formulations for plasma are derived. The high accuracy and efficiency of the presented method is confirmed by computing the transmission and reflection coefficients for a unmagnetized collision plasma slab. The backward radar cross section (RCS) of perfectly conducting sphere covered by homogeneous and inhomogeneous plasma is calculated. Keywords Finite-difference time-domain (FDTD) . Trapezoidal recursive convolution (TRC) . Unmagnetized plasma . Radar cross section (RCS)
1 Introduction The finite-difference time-domain method (FDTD) has been widely used to simulate the transient solutions of electromagnetic problems involving the analysis and design of microwave structures, many other engineering applications and the electromagnetic wave propagation in various media. Over the past decade, there have been numerous investigations of FDTD dispersive media formulations. These include the recursive convolution (RC) methods [1, 2], frequency-dependent Z transform methods [3, 4], piecewise linear recursive convolution (PLRC) method [5], trapezoidal recursive S. Liu : S. Liu School of Sciences, Nanchang University, Nanchang, Jiangxi 330031, People’s Republic of China S. Liu (*) : S. Liu College of Information Science & Technology, Nanjing University of Aeronautics and Astronautics, Nanjing, Jiangsu 210016, People’s Republic of China e-mail:
[email protected]
J Infrared Milli Terahz Waves (2010) 31:620–628
621
convolution (TRC) method [6, 7], JE convolution (JEC) method [8], and piecewise linear current density recursive convolution (PLJERC) method [9], and so on. Although the PLRC-FDTD method attains higher accuracy than RC-FDTD method, it requires complicated formulations of two convolution integrals. The TRC method requires single convolution integral in the formulation as in the recursive convolution (RC) method, while maintaining the accuracy comparable to the piecewise linear convolution integral (PLRC) method. However, that the application of the TRC-FDTD method has been limited to the one dimensional of the Debye and Lorentz models electromagnetic problem [6, 7, 10]. In this letter, we introduce the TRC technique into the plasma media. The propagation of electromagnetic waves through unmagnetized collision plasma slab is also studied using the TRC-FDTD method. The high efficiency and accuracy of the method are confirmed by computing the transmission and reflection coefficients through a unmagnetized plasma layer. The numerical results indicate that the TRC-FDTD method not only improves accuracy over the RC-FDTD with the same computational efficiency but also spends less computational time than the PLRC-FDTD based on the same accuracy. The 3-D TRC-FDTD formulations for plasma are derived. The backward radar cross section (RCS) of conducting sphere covered by homogeneous and inhomogeneous plasma is calculated.
2 Methodology Considering an isotropic collision plasma medium. The Maxwell’s equations and constitutive relation for unmagnetized collision plasma are given by @D ¼ r H; @t
ð1Þ
@H 1 ¼ r E; @t m0
ð2Þ
DðwÞ ¼ "0 "r ðwÞEðwÞ;
ð3Þ
where, E is the electric field, H is the magnetic intensity, ε0 and μ0 are the permittivity and permeability of free space, respectively, εr is the relative permittivity of plasma. In the time domain we have Z t DðtÞ ¼ "0 EðtÞ þ "0 Eðt t Þcðt Þdt ð4Þ 0
Where c(τ) is the electric susceptibility. Using Yee’s notation, we let t = nΔt in (4), and each vector component of D can be written as Z Di ðtÞ ¼ Di ðnΔtÞ ¼ "0 Ein þ "0 Where i = x, y, z.
0
nΔt
Ei ðnΔt t Þcðt Þdt
ð5Þ
622
J Infrared Milli Terahz Waves (2010) 31:620–628
For the TRC technique [6], (5) is approximated using an average of the electric fields over two consecutive time steps in the following form Dni ¼ "0 Ein þ "0
Z n1 nm X E þ E nm1 i
i
2
m¼0
ðmþ1ÞΔt
c ðt Þdt
ð6Þ
mΔt
The value of D at the next time step is ¼ "0 Einþ1 þ "0 Dnþ1 i
Z n X E nþ1m þ Enm i
2
m¼0
ðmþ1ÞΔt
i
cðt Þdt
ð7Þ
mΔt
Following Yee, we quantize space, so that x = iΔx, y = jΔy and z = kΔz. Then we have from (1), using Yee’s notation and unit cell, Dnþ1 Dni nþ1=2 i ¼ ðr HÞi Δt
ð8Þ
We spatially quantize the susceptibility in each cell, so that cðt; x; y; zÞ ¼ c ðt; i; j; k Þ. When (6) and (7) are substituted in (8), we find Z n1 nm1 X E nþ1 þ Ein Δt Ei þ Einm "0 Einþ1 Ein þ "0 i c ðt; i; j; k Þdtþ"0 2 2 0 m¼0 "Z # Z ðmþ1ÞΔt ðmþ2ÞΔt nþ1=2 cðt; i; j; k Þdt cðt; i; j; k Þdt ¼ Δtðr HÞi ðmþ1ÞΔt
ð9Þ
mΔt
To simplify our notation, we let Z c ði; j; k Þ ¼ m
ðmþ1ÞΔt
cðt; i; j; k Þdt
ð10Þ
mΔt
Then (9) becomes n1 P Einm1 þEinm Enþ1 þEn "0 Einþ1 Ein þ "0 i 2 i c0 ði; j; k Þ "0 2 m¼0 nþ1=2
ð11Þ
½ cm ði; j; k Þ cmþ1 ði; j; k Þ ¼ Δtðr HÞi To further simplify our notation we define Δcm ¼ cm cmþ1
ð12Þ
Then (11) becomes n1 nm1 X E nþ1 þ Ein 0 Ei þ Einm m nþ1=2 "0 Einþ1 Ein þ "0 i ð13Þ c "0 Δc ¼ Δt ðr HÞi 2 2 m¼0
Thus we define y ni ¼
n1 nm1 X E þ E nm i
m¼0
i
2
Δcm
ð14Þ
J Infrared Milli Terahz Waves (2010) 31:620–628
623
Substituting (14) into (13) yields 0
Einþ1 ¼
1 c2
0
1 þ c2
Ein þ
1 0
1 þ c2
y ni þ
Δt "0 þ "02c
nþ1=2
0
ðr HÞi
ð15Þ
From(15), the discretized electric field components can be rewritten 0 1 c2 n 1 1 1 1 n ; j; k þ ; j; k Exnþ1 i þ ; j; k ¼ E i þ y i þ 0 0 x x 2 2 2 1 þ c2 1 þ c2 1 Δt 1 1 1 1 1 nþ 2 nþ 12 þ H ; j þ ; k H ; j ; k i þ i þ z z 0 Δy 2 2 2 2 "0 1 þ c2 1 1 1 1 1 nþ 1 nþ 1 Hy 2 i þ ; j; k þ Hy 2 i þ ; j; k Δz 2 2 2 2
ð16Þ 0 1 c2 n 1 1 1 1 n ; k þ ; k Eynþ1 i; j þ ; k ¼ E i; j þ y i; j þ 0 0 y y 2 2 2 1 þ c2 1 þ c2 1 Δt 1 1 1 1 1 nþ 2 nþ 12 þ H ; k þ ; k i; j þ i; j þ H x x 0 Δz 2 2 2 2 "0 1 þ c2 1 1 1 1 1 nþ 1 nþ 1 Hz 2 i þ ; j þ ; k Hz 2 i ; j þ ; k Δx 2 2 2 2
ð17Þ 0 1 c2 n 1 1 1 1 n Eznþ1 i; j; k þ E i; j; k þ y i; j; k þ ¼ þ 0 0 z z 2 2 2 1 þ c2 1 þ c2 Δt 1 1 1 1 1 nþ 12 nþ 12 H ; j; k þ ; j; k þ i þ i þ H y y 0 Δx 2 2 2 2 "0 1 þ c2 1 1 1 1 1 nþ 1 nþ 1 Hx 2 i; j þ ; k þ Hx 2 i; j ; k þ Δy 2 2 2 2
ð18Þ The complex permittivity ε(w) and susceptibility c(w) for an isotropic plasma are give by ! w2p ¼ "0 ð1 þ cðwÞÞ ð19Þ "ðwÞ ¼ "0 1 þ wð jn wÞ pffiffiffiffiffiffiffi Where v is the collision frequency, wp is the radian plasma frequency, and j ¼ 1.The time domain susceptibility function obtained by computing the inverse Fourier transform of c(w) is c ðt Þ ¼
w2p ½1 expðnt ÞU ðt Þ n
ð20Þ
624
J Infrared Milli Terahz Waves (2010) 31:620–628
Where U(τ) is the unit step function. Substituting (20) into (10) yields c0 ¼
w2p w2p Δt 2 ½1 expðnΔtÞ n n
w2p Δc ¼ n m
ð21Þ
!2 expðmnΔt Þ½1 expðnΔtÞ2
ð22Þ
The recursion relation is given by Δcmþ1 ¼ expðnΔt ÞΔcm y ni ¼
ð23Þ
Δc0 n Δc0 n1 E þ E þ expðnΔtÞy n1 i 2 i 2 i
ð24Þ
From(24), the discretized recursion relation can be rewritten 1 Δc 0 n 1 Δc0 n1 1 1 Ex i þ ; j; k þ Ex y nx i þ ; j; k ¼ i þ ; j; k þ enΔt y xn1 i þ ; j; k 2 2 2 2 2 2 ð25Þ 1 Δc 0 n 1 Δc0 n1 1 1 Ey i; j þ ; k þ Ey y ny i; j þ ; k ¼ i; j þ ; k þ enΔt y yn1 i; j þ ; k 2 2 2 2 2 2
ð26Þ 1 Δc0 n 1 Δc0 n1 1 1 ¼ þ þ enΔt y n1 Ez i; j; k þ Ez y nz i; j; k þ i; j; k þ i; j; k þ z 2 2 2 2 2 2 ð27Þ H field is updated as traditional FDTD algorithms.
3 Verification and accuracy To demonstrate the aforementioned TRC-FDTD formulation for unmagnetized plasmas, we will give two examples. [Case 1] We calculated the transmission and reflection coefficients of electromagnetic wave through an unmagnetized collision plasma slab with a thickness of 9 mm. In order to eliminate zero frequency incident energy, calculations were made for a Gaussian-derivative pulse normally incident to the plasma slab. "
ðt 5t Þ2 Ei ðtÞ ¼ ðt 5t Þ exp 2t 2 where t =15Δt.
# ð28Þ
J Infrared Milli Terahz Waves (2010) 31:620–628
625
The computational domain considered had 520δ in z direction, where δ is the mesh width taken to be δ=75 μm. The time step is Δt=0.125 ps. The plasma slab occupies cells 200 to 320, five cells perfectly matched layer (PML) [11] were used at the terminations of the space to eliminate unwanted reflections, and other is free space. The plasma parameters were wp ¼ 2p 28:7Grad=s, ν=20 GHz. The transmission and reflection coefficients are calculated by using the fast Fourier transform (FFT) method. The simulations are allowed to run for 10000 time steps and it is convergent. Figures 1 and 2 illustrate the magnitudes of the transmission coefficients and reflection coefficients computed using the TRC-FDTD, PLRC-FDTD, and RC-FDTD method with those of the analytical solution. Figures 1 and 2 depict that the TRC-FDTD and PLRCFDTD method are in good agreement with the analytical magnitude. But the RC-FDTD method exhibits slight error. Moreover, the memory and the CPU time used of computing by an AMD sempron (tm) processor 3500+ based computer for the three different methods is shown in Table 1. The numerical results indicate that the TRC-FDTD method not only improves accuracy over the RC-FDTD with the same computational efficiency but also spends less computational time than the PLRC-FDTD based on the same accuracy. [Case 2] We calculated the backward RCS of perfectly conducting sphere covered by homogeneous and inhomogeneous plasma. The distribution function of inhomogeneous plasma electron density " ne ðrÞ ¼
Nemax
ðr a d Þ2 2d ðr a d Þ d2
# ð29Þ
where ne(r) is the electron density, Nemax is the peak value electron density, d is the thickness of the plasma, r is the radial variable. The radius of the conductive sphere is a=0.1 m.
Fig. 1 Reflection coefficient magnitude versus frequency. Reflection coefficient magnitude (dB)
0
-10
-20
TRC method analytical solution PLRC method RC method
-30
-40
-50
0
10
20
30
40
Frequency (GHz)
50
60
70
626
J Infrared Milli Terahz Waves (2010) 31:620–628
Fig. 2 Transmission coefficient magnitude versus frequency. Transmission coefficient magnitude (dB)
0 -10 -20 -30
TRC method analytical solution PLRC method RC method
-40 -50 -60 -70 -80 0
10
20
30
40
50
60
70
Frequency (GHz)
Normally incidence on conducting sphere, being propagation in z-axis, the incident wave used in the simulation is a Gaussian pulsed wave " # ðt 5t Þ2 ð30Þ Ei ðtÞ ¼ exp 2t 2 where t =15Δt. The FDTD spatial discretisation was fixed to be d ¼ Δx ¼ Δy ¼ Δz ¼ 0:005m, the temporal step is Δt ¼ 0:5d=c, where c is the speed of light in vacuum. The FDTD total computation domain is 100×100×100δ3. The simulations are allowed to run for 1150 time steps. The plasma parameters were wp ¼ 4p 109 rad=s, v=20 GHz. Figure 3 shows the backward RCS of conducting sphere covered with homogeneous plasma. It is clear from Fig. 3 that when the electromagnetic (EM) wave frequency f<1.5GHz, the plasma can increase RCS. Physically, this can be explained on the basis that isotropic plasma acts as a “high-pass filter”. Namely, EM wave frequencies above the plasma frequency constitute a pass band, and frequencies below the plasma frequency constitute a stop band. Alternatively, when the EM wave frequency 1.5
