Proceedings of International Forum for Green Asia 2017
Special issue of Evergreen - Joint Journal of Novel Carbon Resource Sciences & Green Asia StrategyP-31

3D FDTD-PIC Simulation of Particles in a Miniature Ion Thruster
Ali Yousefian,1 Naoji Yamamoto1 and Takahiko Miyazaki 1,2
1Interdisciplinary Graduate School of Engineering, Kyushu University
2
Green Asia Education Center, Kyushu University
Kasuga-koen, Kasuga, Fukuoka 816-8580, Japan,

Abstract
Our research focuses on simulation of particles inside an Electric Ion thrusters. An Ion engine works bycreating Plasma through ECR (Electron Cyclotron Resonance) method and accelerating the resultant ionsout of the engine thus producing thrust. It is, however, difficult to measure inner properties of the plasma,since the ion engine is so small that we cannot insert measurement equipment without perturbations.
Numerical simulations are effective tool to understand the phenomena occurring in a miniature microwave
ion thruster. This research focuses on development of a three dimensional FDTD (Finite Difference TimeDomain) approximation to Maxwell’s equations simulation for propagation of TEM (TransverseElectro-Magnetic) microwaves in a Coaxial Cable and development of PIC (Particle In Cell) simulation forparticle interactions inside the miniature ion thruster.

1. Introduction
Current research involves the development of an
all-inclusive full PIC-FDTD code in order tosimulate the particles inside a miniature ionthruster. As mentioned earlier it is very challenging
to measure inner properties of the plasma inside
the very small thruster, a fully developed
Numerical simulations code which utilizes the PIC
method as well as FDTD would shed light on theinner workings of our engines. FDTD (FiniteDifference Time Domain) simulation method isextensively used in plasma simulation as a way to
calculate the fluctuations of electric and magnetic
field values with respect to space and time.
Therefor this method is utilized to simulate
propagation of electromagnetic waves and charged
particles. Coaxial cables are used in order to feed
the microwave power to the engine and an accurate
simulation of microwave propagation inside a
coaxial cable leads to a more realistic simulation ofparticles for the entire engine. Figure 1 shows one
such engine developed in our research group inKyushu University with the coaxial cable attachedto the thruster. Our current research involves
simulation of microwaves in the coaxial cable aswell as simulation of neutral particles in theengine.

Fig. 1: A miniature Ion thruster developed in
Yamamoto lab In Kyushu University

2. FDTD simulation
2.1 Electromagnetic Fields
Since the three dimensional Cartesian model isused for this simulation, each grid is practically aperfect cube. One such grid along with the locationof the field points is illustrated in Figure 4. Since
none of the field points fits on the grid, a 3D
FDTD simulation is very challenging in its
calculation and indexing. Adjusting the original
Maxwell formulas to the simulation model inFigure 4 is a challenging task. Considering thebasic scalar Maxwell equations in threedimensions:85Proceedings of International Forum for Green Asia 2017
Special issue of Evergreen - Joint Journal of Novel Carbon Resource Sciences & Green Asia Strategy

P-31𝑥𝑟𝑒𝑞𝑢. 9𝑦𝑟𝑒𝑞𝑢. 10𝐸𝑥(𝑖,𝑗,𝑘) = 𝐸𝑥(𝑖,𝑗,𝑘) + 𝐸𝑟(𝑖,𝑗)𝜕𝐵𝑥 𝜕𝐸𝑦 𝜕𝐸𝑧=−𝜕𝑡𝜕𝑧𝜕𝑦𝑒𝑞. 1𝐸𝑦(𝑖,𝑗,𝑘) = 𝐸𝑦(𝑖,𝑗,𝑘) + 𝐸𝑟(𝑖,𝑗)𝜕𝐸𝑥 1 𝜕𝐵𝑧 𝜕𝐵𝑦𝜀= (−)
𝜕𝑡𝜇 𝜕𝑦𝜕𝑧𝑒𝑞. 2

By transforming these equations to the grid point
system in Figure 2, Maxwell’s equations can beutilized in the simulation. Equations 3 and 4
indicate such transformation for 𝐸𝑥 and 𝐵𝑥 :𝐸𝑛
𝜀𝑑𝑖− 𝐸 𝑛−1 11𝑥(𝑖+2,𝑗,𝑘)𝑥(𝑖+2,𝑗,𝑘)
𝑑𝑡
𝐵−B

=

1𝑛−21
(𝜇01𝑛−− 𝐵 21 1
1 1𝑧(𝑖+2,𝑗+2,𝑘)𝑧(𝑖+2,𝑗−2,𝑘)𝐵

1𝑛−2𝑑𝑦1𝑛−211 −𝐵11𝑦(𝑖+2,𝑗,𝑘+2)𝑦(𝑖+2,𝑗,𝑘−2)

1n+2

11 −B11x(i,j+ ,k+ )x(i,j+ ,k+ )2
222dtEn−)

𝑑𝑧

1n−2

𝑒𝑞. 3

En=− En 11y(i,j+2,k+1)y(i,j+2,k)
dz

n1 −E
1z(i,j+1,k+2)z(i,j,k+2)

Fig. 2: a 3D FDTD grid with field positions

𝑒𝑞. 4

dy

And finally by adjusting the indexing of the above
equations for the simulation code we will have:𝐸𝑥(𝑖,𝑗,𝑘) = 𝐸𝑥(𝑖,𝑗,𝑘) + (𝐵𝑧(𝑖,𝑗,𝑘) − 𝐵𝑧(𝑖,𝑗−1,𝑘) )𝐶𝑦− (𝐵𝑦(𝑖,𝑗,𝑘) − 𝐵𝑦(𝑖,𝑗,𝑘−1) )𝐶𝑧
𝐵𝑥(𝑖,𝑗,𝑘) = 𝐵𝑥(𝑖,𝑗,𝑘) + (𝐸𝑦(𝑖,𝑗,𝑘+1) − 𝐸𝑦(𝑖,𝑗,𝑘) )Where 𝐶𝑦 =

1 𝑑𝑡

1 𝑑𝑡
2 𝑑𝑦

11and 𝐶𝑧 =1 𝑑𝑡

2 𝑑𝑧 𝜀𝑑𝑖 𝜇0𝑒𝑞. 6.

2.2 Microwave Input Method
The microwave input power in a coaxial cable iscalculated by:𝑝𝑜𝑤𝑒𝑟 =

2𝜋 2𝑟𝑜𝑢𝑡𝑒𝑟𝑣 ln ()𝑟𝑖𝑛𝑛𝑒𝑟𝜇0 0𝑏√𝜀𝑑𝑖𝑒𝑞. 7

Where 𝑟𝑜𝑢𝑡𝑒𝑟 and 𝑟𝑖𝑛𝑛𝑒𝑟 are the outer and inner
radius of the cable respectively and the 𝑣0𝑏 is thevoltage between the cable’s inner and outer layers.
Therefore input field values are calculated by:𝐵𝜃 = 𝑣0𝑏11𝑟 𝑐𝐷,𝐸𝑟 = 𝑣0𝑏1𝑟𝑦𝑟𝑒𝑞𝑢. 11𝐵𝑦(𝑖,𝑗,𝑘) = 𝐵𝑦(𝑖,𝑗,𝑘) + 𝐵𝜃(𝑖,𝑗)𝑥𝑟𝑒𝑞𝑢. 12The simulated waves are illustrated in figures 3through 8. The results in these figures are taken at
a mid-plane cross section of the cable where the yvalues are constant. Since the wave in thissimulation is supposed to act as a TEM wave the𝐸𝑧 and 𝐵𝑧 values are expected to be close to zero.
As it can be seen in figure 5 the comparative valueof 𝐸𝑧 with respect to the amplitude of 𝐸𝑥 wavein figure 3 is too insignificant (24000 timessmaller) which means that it can be ignored. Samething applies to the magnetic field values depicted
in figures 7 and 8. The 𝐵𝑧 value is half a million
times smaller than By which means it will not affectthe accuracy of the simulation as well. Since this
output data is taken at a plane in the center of the yaxis, Ey and 𝐵𝑥 values are also expected to be zero.
However as mentioned earlier since the simulated
geometry of the cable is not a perfect cylinder we
expect some imprecisions as well. According tofigures 4 and 6 since the amplitudes of Ey and Bx
waves are four orders of magnitude smaller thantheir respective wave amplitudes, theseimperfections can also be ignored as they will notcause any defects on the outcome of thesimulation.

1 𝑑𝑡
2 𝑑𝑧

− (𝐸𝑧(𝑖,𝑗+1,𝑘) − 𝐸𝑧(𝑖,𝑗,𝑘) )2 𝑑𝑦 𝜀𝑑𝑖 𝜇0𝑒𝑞. 5𝐵𝑥(𝑖,𝑗,𝑘) = 𝐵𝑥(𝑖,𝑗,𝑘) − 𝐵𝜃(𝑖,𝑗)𝑒𝑞𝑢. 8Where 𝐶𝑑 is the speed of light in the dielectric
region. The field components 𝐸𝑟 and 𝐵𝜃 shouldbe implemented on the cables cross section (Figure3) and should be divided into their 𝑥𝑦
components. In order to get the right 𝐸 × 𝐵
vector and the right values, therefore the following
formulation is used:86Proceedings of International Forum for Green Asia 2017
Special issue of Evergreen - Joint Journal of Novel Carbon Resource Sciences & Green Asia StrategyP-31

Fig. 3: 𝐸𝑥 field values through the coaxial cable

Fig. 6: 𝐵𝑥 field values through the coaxial cable

Fig. 4: 𝐸𝑦 field values through the coaxial cable

Fig. 7: 𝐵𝑦 field values through the coaxial cable

Fig. 5: 𝐸𝑧 field values through the coaxial cable

Fig. 8: 𝐵𝑧 field values through the coaxial cable87Proceedings of International Forum for Green Asia 2017
Special issue of Evergreen - Joint Journal of Novel Carbon Resource Sciences & Green Asia Strategy

P-31

Reference

3. PIC simulation
For the analysis of plasma behavior in thedischarge chamber, the PIC method is used. Effectsof electromagnetic field on plasma particles arecalculated by linear interpolation in the PIC
method. The PIC method can treat collision
process and distribution function since in thismethod the simulation space is divided into meshesin which individual particles are tracked along
these meshes while variables such as density andcurrent are calculated at these mesh points. Using
this method the particle interaction with the engine
and each other such as reflection on the engine
walls and collision of particles is also simulated. In[1]. Francis F.Chen, in Introduction to Plasma
Physics and Controlled Fusion, 2nd ed. byPlenum Publishing Corporation, (1974).
[2]. C.K.Birdsall and A.B.Langdon, in Plasma
Physics via Computer Simulation, by IoP
Publishing (1985).
[3]. Yee K 1966 “Numerical solution of initialboundary
valueproblems
involvingmaxwell’s equations in isotropic media”
IEEE Trans. Antennas Propagat.

order to achieve that, first we have simulated theinteraction of neutral particles in the engine. Figures
9 shows the particles distribution of neutral particles
in a 3D engine. In order to simplify the outcome
demonstration, these results are taken at a central
plane rather than the whole engine.

[4]. Kousaka.H, Ono.K 2002 “ Numerical
Analysis of the electromagnetic fields in a
Microwave Plasma Source Excited by
Azimuthally Surface Waves” Jpn J. Phys.
Vol 41[5]. Takao Y, Koizumy H, Eriguchi K,Komurasaki K, and Ono K 2012
“Three-dimensionalParticle
inCellSimulation of Miniature Plasma Source for aMicrowave Discharge Ion Thruster ” Plasma
Source Sci. Technol Vol.23

4. Conclusion
The simulation for a microwave propagation insidea coaxial cable was discussed using FDTD method
in a 3D model. The simulation results werediscussed and their validation was investigated.
Furthermore the simulation results of a 3D PIC
simulation of neutral particles in miniature ionthruster was discussed. These two codes will beutilized in future works for a more realistic
simulation of an Ion thruster. This new simulation
model will replace the previous inaccurate method
for future simulations in our research group andThis in turn will result in development of better
performing engines in future.

[6]. Koizumi H and Kuninaka H 2010
“Miniature microwave discharge ion thruster
driven by 1 watt microwave power” J.
Propul. Power 26 601-4
Email: 3ES17308W@s.kyushu-u.ac.jp

Fig. 9: Neutral particle distribution in a miniature Ion thruster
developed in Yamamoto lab In Kyushu University

88