Thrust calculation of electric solar wind sail by particle-in-cell simulation

Abstract

In this study, thrust characteristics of an electric solar wind sail were numerically evaluated using full three-dimensional particle-in-cell (PIC) simulation. The thrust obtained from the PIC simulation was lower than the thrust estimations obtained in previous studies. The PIC simulation indicated that ambient electrons strongly shield the electrostatic potential of the tether of the sail, and the strong shield effect causes a greater thrust reduction than has been obtained in previous studies. Additionally, previous expressions of the thrust estimation were modified by using the shielded potential structure derived from the present simulation results. The modified thrust estimation agreed very well with the thrust obtained from the PIC simulation.

Ann. Geophys., 34, 845855, 2016 www.ann-geophys.net/34/845/2016/ doi:10.5194/angeo-34-845-2016 Author(s) 2016. CC Attribution 3.0 License.

Thrust calculation of electric solar wind sail by particle-in-cell simulation

Kento Hoshi1, Hirotsugu Kojima2, Takanobu Muranaka3, and Hiroshi Yamakawa2

1Department of Electrical Engineering, Graduate School of Engineering, Kyoto University, Kyoto, Japan

2Research Institute for Sustainable Humanosphere, Kyoto University, Kyoto, Japan

3Department of Electrical Engineering, Chukyo University, Nagoya, Japan

Correspondence to: Kento Hoshi ([email protected])

Received: 1 July 2016 Revised: 25 August 2016 Accepted: 7 September 2016 Published: 26 September 2016

Abstract. In this study, thrust characteristics of an electric solar wind sail were numerically evaluated using full three-dimensional particle-in-cell (PIC) simulation. The thrust obtained from the PIC simulation was lower than the thrust estimations obtained in previous studies. The PIC simulation indicated that ambient electrons strongly shield the electro-static potential of the tether of the sail, and the strong shield effect causes a greater thrust reduction than has been obtained in previous studies. Additionally, previous expressions of the thrust estimation were modied by using the shielded potential structure derived from the present simulation results. The modied thrust estimation agreed very well with the thrust obtained from the PIC simulation.

Keywords. General or miscellaneous (instruments useful in three or more elds; new elds (not classiable under other headings); techniques applicable in three or more elds)

1 Introduction

An electric solar wind sail, called the E-sail, is a recently proposed propulsion device that consists of 50100 conductive tethers with lengths of 1020 km and thicknesses of 0.1 1 m. The E-sail was rst proposed by Janhunen (2004). The main body of the spacecraft expands the tethers to form a sail-like structure. The E-sail has electron guns to maintain a positive surface potential on the order of several kilovolts, in order to deect solar wind protons. The tethers obtain the momentum of these deected protons via Coulomb scattering and use it as their propulsive force. The system requires electron sources and electrical power for the electron guns to

produce thrust. The E-sail is expected to be used as a new propellantless space propulsion device.

The thrust characteristics of the E-sail were rst investigated by Janhunen and Sandroos (2007). They performed a one-dimensional (1-D) particle-in-cell (PIC) simulation of a conductive tether with a radius of 1.0 m. They found that an ansatz of the electrostatic potential structure around the tether, which is expressed as

V (r) =

V0 2

1 + (2 D/r)2 ln(2 D/rw) , (1)

agreed very well with the result of their PIC simulation, where r is the distance from the tether, V0 is the surface potential of the tether, D =

ln

p"0kBTe/e2ne is the Debye length of the electron, and rw is the radius of the tether. For evaluating the performance of an E-sail, the thrust per unit length is often used. The total thrust can be calculated from the thrust per unit length by multiplying the number of tethers by the length of one tether. Janhunen and Sandroos (2007) also conducted a two-dimensional (2-D) PIC simulation and suggested that the thrust per unit length acting on the tether is

dF dL =

Kmpn0v2dr0

sexp

, (2)

where L is the length of the tether, K is the coefcient of proportionality (K 3.09 was obtained from their Monte Carlo

simulation), mp is the mass of a proton, n0 is the ambient plasma density, and vd is the drift velocity of the solar wind.Moreover, Janhunen (2009) suggested the thrust of the E-sail

Published by Copernicus Publications on behalf of the European Geosciences Union.

mpv2deV0 ln( r0rw )

[bracketrightbigg] 1

846 K. Hoshi et al.: Thrust calculation of electric solar wind sail

may increase because of a lack of electrons around the tether.

He proposed the use of multiple tethers to collect ambient electrons so that the electron density around each tether decreases and the ambient electrons cannot completely shield the potential of the tether. According to Janhunen (2009), the thrust with ambient electron removal is 5 times larger than that obtained by Janhunen and Sandroos (2007) without ambient electron removal.

Sanchez-Torres (2014) also investigated the thrust of the E-sail considering Coulomb scattering, assuming the absence of trapped electrons around the tether. He expressed the potential structure as

V (r) = V0

8 >

>

<

>

>

:

Table 1. Simulation parameters (BG: background).

BG plasma species Electron, proton

BG plasma density n0 107 m3 BG electron temperature Te 100 eVBG proton temperature Tp 12.0 eV Drift velocity vd of BG plasma 400 km s1

Electron Debye length D 23.5 m Time step width [Delta1]t 20 ns

Cell spacing [Delta1]x 0.30 mParticle number per cell pcell 60 particles cell1

(30 electrons, 30 protons)

System size 512 512 512 cells

= (153.6 m)

3

ln

rsheb

r

(r rsheb)

0 (r > rsheb)

ln

rsheb

rw

present simulation results. The estimated thrust and the PIC results are found to be in good agreement with each other and lower than those estimated in previous studies.

2 Full PIC simulation of E-sail

2.1 Simulation settings

This section describes the PIC simulation congurations of a positively charged tether in the solar wind environment.The simulation code HiPIC, which was developed by the Japan Aerospace Exploration Agencys Engineering Digital Innovation Center (JEDI) (Muranaka et al., 2011), was used to perform this simulation. HiPIC is an electrostatic code that models 3-D rectangular cells in space and uses the full PIC method to calculate collisionless kinetic plasma. HiPIC solves Newtons equations of motion for each particle using the BunemanBoris method and solves Poissons equations to obtain the electric potential structure in the computational domain using a discrete sine transformation. HiPIC can be used to calculate the interaction between plasmas and the spacecraft, which is modeled with rectangular internal boundaries. A detailed description of the performance of the code is given in Muranaka et al. (2011).

The simulation and physical parameters are given in Table 1. The electron density ne and the proton density np are ne = np = n0 = 1.0107 m3, and their temperatures are

kBTe = 100 eV and kBTp = 12 eV, respectively. The back

ground plasmas have a solar wind drift velocity of vd =

400 km s1 along the x axis. All of the edges of the simulation domain were xed to V = 0V (Dirichlet boundary con

dition). The background plasma particles were injected from all the domain boundaries in each time step as many times as the number of outgoing particles in previous time step.

A tether-like rectangular model was set in the center of the computational domain, as shown in Fig. 1. The size of the model is 12001 cells, which is equivalent to a tether with

dimensions of rw = 15 cm and L = 60 m. Tethers with L =

Ann. Geophys., 34, 845855, 2016 www.ann-geophys.net/34/845/2016/

, (3)

where the approximated parameter b is 0.65 for a potential bias of V0 =1040 kV and rsh is the sheath radius for a highly

positive bias tether and can be calculated from the ambient plasma parameters and V0 (Sanmartn et al., 2008). Using the solar wind parameters at 1 AU, the thrust per unit length is 407 nN m1 for a 20 kV charged tether of rw = 20 m. This is

lower than the thrust per unit length of 500 nN m1 estimated by Janhunen (2009) but higher than the value of 100 nN m1 estimated by Janhunen and Sandroos (2007) because Eq. (3) yields greater values of the potential than Eq.(1) under almost all conditions. The study by Sanchez-Torres (2014) was purely analytical; no plasma simulations were performed.

However, Hoshi et al. (2016) showed that the actual electrostatic potential structure around the tether was lower than those given by Eqs. (1) and (3) based on the results of a full three-dimensional (3-D) PIC simulation with V0 = 240 V.

The thrust of the E-sail is generated from the deection of solar wind protons by the electrostatic potential. If the electrostatic potential derived from the tethers is greatly shielded by ambient electrons, the actual thrust is lower than that estimated in previous studies. Hoshi et al. (2016) did not consider the thrust because a potential of 240 V was not sufcient to deect solar wind protons with a drift velocity of approximately 400 km s1 ( 0.8 keV).

In the present paper, we performed 3-D full PIC simulations to simulate a transient of the thrust of the E-sail. The thrust found in this paper is lower than that obtained in previous studies, with a sufciently high potential to deect ambient protons (V0 > 1 kV). Section 2 discusses the 3-D PIC simulation with V0 4.0 kV that was performed to conrm

that the potential was lower than those obtained from Eqs. (1) and (3). The propulsive force that acts on the tether is also calculated in the PIC simulation. In Sect. 3, the thrust is numerically estimated and compared with the PIC results. Two estimation procedures employed in previous studies (Janhunen and Sandroos, 2007; Sanchez-Torres, 2014) are modied to contain a shielded potential structure derived from the

K. Hoshi et al.: Thrust calculation of electric solar wind sail 847

the ambient particles had not yet begun to respond to the electrostatic potential of the tether. Fx then increased with time and converged to a specic value. At the end of the simulation, Fx was 0.35, 0.76, and 0.96 N for L = 30, 60, and

75 m, respectively. These values represent the total force acting on the tether, including the sum of the particles hitting the whole tether and the force calculated from the Maxwell stress tensor. However, the thrust per unit length indicates the performance of the E-sail; thus, Fx/L was calculated, as shown in Fig. 2b. Fx/L was 11.7, 12.7, and 12.8 nN m1 for

L = 30, 60, and 75 m, respectively. The thrust per unit length

for L = 30 m was slightly smaller than that for 60 and 75 m.

This lower value of Fx/L may be due to the end effect or the effect which arises when L is small in comparison with D. At L = 60 and 75 m, the values of Fx/L were almost

equal; thus, L = 60 m is considered to be sufciently long to

simulate an innite tether.The kinks in the thrust between t = 10 and 20 s shown

in Fig. 2 correspond to collections of ambient electrons. Figure 3 shows the time history of the current on the surface of the tether with L = 60 m. The ambient electron current (pur

ple line) varied dramatically between t = 10 and 20 s. This

is an initial response of the ambient electrons to the potential of the tether. The electron plasma frequency !p was approximately 178 kHz so the response time of ambient electrons were approximately 5.6 s. Collected electrons do not directly contribute to the thrust, because their momentum is small. Instead, they temporarily shield the potential structure more strongly than in a steady state, causing the thrust to stop increasing with time, as shown in Fig. 2.

Figure 3 also reveals that the dominant source of the force acting on the tether was the force from the Maxwell stress tensor, not from the protons hitting the tether because the proton current was almost zero. This fact shows that the computation successfully simulates the thrust generation by proton deection.

To compare the present results with previous thrust estimations, Fx/L was calculated at various V0 values, as shown in Fig. 4. The blue line in Fig. 4 is the result of our PIC simulation with L = 60 m and rw = 15 cm. The black

and red dotted lines show the thrust estimated by Janhunen and Sandroos (2007) and Sanchez-Torres (2014), respectively, with rw = 15 cm. In the present simulation, Fx/L was

67.1 nN m1 at V0 = 4.0 kV; in contrast, Janhunen and San

droos (2007) and Sanchez-Torres (2014) estimated Fx/L to be 123 and 246 nN m1, respectively. These results indicate that the thrust characteristics of the E-sail are different from those of conventional estimations. This difference is likely due to the difference in the potential structure around the tether.

Figure 5 compares the electrostatic potential structures obtained in the present study and two previous studies at V0 = 1.0, 2.0, 3.0, and 4.0 kV. The purple lines show the PIC

simulation results for the potential at the center of the tether (y = 76.8 m, z = 76.65 m). At all potentials, the potential ob-

www.ann-geophys.net/34/845/2016/ Ann. Geophys., 34, 845855, 2016

v = 400 km

x = 30 cm

Figure 1. Denition of the tether model. A tether-like rectangle is located in the center of the computational domain. The solar wind originates from x = 0.

30 and 75 m were also simulated to conrm the inuence of the length of the tether on the thrust per unit length.

Because of the limitation of the calculation resources, we cannot include the emission of the electron beam from the tethers edge and simulate the self-charging of the tether. Instead of emitting the electron, the surface potential V0 of the tether was xed to an inputted value. Hoshi et al. (2016) showed that effects of emitted electrons on the potential structure were small, so we consider that the absence of emitted electrons do not cause signicant differences in the force acting on the tether.

V0 was varied from 0 to 4.0 kV, and the thrust acting on the tether Fx was calculated. Fx is the x component of the thrust and was calculated as the sum of the total momentum of the particles impinging on the tether during each time step divided by [Delta1]t and the Coulomb force calculated from the Maxwell stress tensor. The simulation progressed with time steps of [Delta1]t = 20 ns until the time variation of the external

force became zero. In almost all the cases, the total iteration was 10 000 steps (= 0.2 ms). An additional 2000 steps were

calculated for the V0 = 4.0 kV case.

To perform the computation, we applied an MPI parallelization and an OpenMP parallelization. Each case used 2048 cores on Cray XE4 for calculation (1024 processes for MPI, two threads for OpenMP), requiring approximately 12 h to run.

2.2 Simulation results

Figure 2a shows the time variation of Fx at V0 = 1.3 kV. At

the beginning of the simulation, Fx was almost zero because

848 K. Hoshi et al.: Thrust calculation of electric solar wind sail

1000

14

a) Thrust

Fx

75 m

60 m

30 m

L = 30 m L = 60 m L = 75 m

b) Thrust per unit length

Fx / L

75 m

30 m

60 m

L = 30 m L = 60 m L = 75 m

900

12

800

-1

Thrust per unit length [nN P ]

700

10

600

Thrust [nN]

8

500

6

400

300

4

200

2

100

0 0 50 100 150 200

0 0 50 100 150 200

Time [s]

Time [s]

Figure 2. Time history of the thrust with V0 = 1.3 kV: (a) x component of the total external force; (b) x component of the thrust per unit

length.

-2.0

present at V0 = 1.0 kV, it was formed by protons impinging

on the surface of the tether. At a very small rw, a zero-density region may appear for V0 = 1.0 kV.

This study then considered the high-density proton region in front of the tether. Figure 7 shows the electron and proton density structure along the x axis at V0 = 2.0 kV. As

with Fig. 3a in Janhunen and Sandroos (2007), Fig. 7b reveals a high-density region in front of the tether. The maximum proton density was 2.0 107 m3 at V0 = 2.0 kV. The

maximum electron density was approximately 1.0108 m3

at V0 = 2.0 kV. The maximum positive charge density of

the high-density region is one fth of the negative charge density of the electrons. The high-density region does not compensate for the potential; that is, it does not reduce the shielding effect of the electrons. Trapped electron removal, which was discussed by Janhunen (2009), was not observed in the present simulation. Thus, the thrust models given by Eqs. (1) and (3) are inappropriate for estimating the thrust of the E-sail, and a new model considering appropriate potential shielding must be developed.

3 Numerical estimation of the thrust of the E-sail

This section describes the proposed method of estimating the thrust of the E-sail and presents the estimation results.First, a semi-analytical solver of the electrostatic potential structure around an innite cylinder in plasma is introduced.This solver was used to obtain a realistic estimation of the electrostatic potential using a 2-D inverse fast Fourier transform (FFT). The two thrust estimation procedures used to evaluate the thrust characteristics of the E-sail are then described. The rst is the effective radius method by Janhunen and Sandroos (2007), and the second is the Coulomb scattering method by Sanchez-Torres (2014). These two methods are hereafter called Methods 1 and 2, respectively. These two evaluation procedures were then used in combination with

Ann. Geophys., 34, 845855, 2016 www.ann-geophys.net/34/845/2016/

0.2

0.0

Proton current

Electron current

Electron (L = 60 m)

Proton (L = 60 m)

-0.2

-0.4

-0.6

Current [mA]

-0.8

-1.0

-1.2

-1.4

-1.6

-1.8

0 50 100 150 200

Time [s]

Figure 3. Time history of the ambient electron and proton currents (V0 = 1.3 kV, L = 60 m).

tained from the present PIC simulations was lower than the potentials obtained using Eqs. (1) and (3). Figure 5 indicates that potential shielding by ambient electrons is not appropriately included in Eqs. (1) and (3). Figure 5 also shows that the sheath length assumed by Eq. (3) is consistent with that obtained by PIC simulation.

Figure 6 shows the proton density structure at t = 0.2 ms.

At V0 = 2.0, 3.0, and 4.0 kV, a zero-density (np = 0 m3) re

gion was present in front of the tether. No protons impinged on the surface of the tether. The momentum of the protons was transferred to the tether through the Coulomb force acting on the tether. In contrast, there was no zero-density region in front of the tether at V0 = 1.0 kV. This is because the

radius rw of the tether is relatively large. The drift velocity of the proton (400 km s1) is equivalent to 0.83keV, but the gure indicates that V0 = 1.0 kV is not sufcient to deect

all of the protons for rw = 15 cm. Although a wake region is

K. Hoshi et al.: Thrust calculation of electric solar wind sail 849

140

120

Estimated thrust per unit length [nN P ]

-1

100

80

60

40

20

0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

Potential [kV]

Figure 4. Comparison of present thrust simulation results with previous estimations. The blue line shows the thrust obtained from the present PIC simulation with L = 60 m and rw = 15 cm. The black and red dashed lines show the thrust estimated by Janhunen and Sandroos (2007)

and Sanchez-Torres (2014), respectively, with rw = 15 cm. The estimation procedures used in these two previous studies are described in

Sect. 3.

0.0

1.0

-80 -60 -40 -20 0 20 40 60 80

2.0

-80 -60 -40 -20 0 20 40 60 80

(a) V0

PIC simulation

Janhunen and Sandroos [Eq. (1)]

Sanchez-Torres [Eq. (3)]

= 1.0 kV

(b) V0

PIC simulation

Janhunen and Sandroos [Eq. (1)]

Sancehz-Torres [Eq. (3)]

= 2.0 kV

0.8

1.5

0.6

Potential [kV]

Potential [kV]

1.0

0.4

0.2

0.5

0.0

Distance from tether [m]

Distance from tether [m]

3.0

-80 -60 -40 -20 0 20 40 60 80

4.0

-80 -60 -40 -20 0 20 40 60 80

(c) V0

PIC simulation

Janhunen and Sandroos [Eq. (1)]

Sanchez-Torres [Eq. (3)]

= 3.0 kV

(d) V0

PIC simulation

Janhunen and Sandroos [Eq. (1)]

Sanchez-Torres [Eq. (3)]

= 4.0 kV

3.5

2.5

3.0

2.0

Potential [kV]

2.5

1.5

Potential [kV]

2.0

1.0

1.5

1.0

0.5

0.5

0.0

0.0

Distance from tether [m]

Distance from tether [m]

Figure 5. Electrostatic potential structure around the tether. Purple lines are the results of the PIC simulations at the center of the tether (y = 76.8 m, z = 76.65 m). Green and blue dashed lines are plots of Eqs. (1) and (3), respectively.

www.ann-geophys.net/34/845/2016/ Ann. Geophys., 34, 845855, 2016

850 K. Hoshi et al.: Thrust calculation of electric solar wind sail

Figure 6. Proton density structure at the center of the tether (y = 76.8 m) at tether surface potentials of (a) V0 = 1.0 kV, (b) V0 = 2.0 kV,

(c) V0 = 3.0 kV, and (d) V0 = 4.0 kV. The ambient plasma conditions are n0 = 10 cm

3, Te = 100 eV, and Tp = 12 eV.

the realistic electrostatic potential given by the 2-D FFT. The thrust estimation script written by Python 2.7 is found in the Supplement (available online).

3.1 Semi-analytical solver of electrostatic potential structure

The Poisson equation in plasma is expressed as

r2V (r) =

(r)

"0

, (4)

where is the space charge density. From the velocity distribution function of electrons, we assume that the electrons density distribution becomes the Boltzmann distribution:

ne(r) = n0 exp

"0 . (7)

The 1-D solution of Eq. (7) is well known and is called the Yukawa potential. However, the 2-D analytical solution of Eq. (7), which would represent the shielded potential structure around an innite tether, remains unknown. This is why previous studies had to assume an articial potential structure, such as those given by Eqs. (1) and (3).

Hoshi et al. (2016) developed a numerical method of calculating the potential structure around a tether in plasma, adding the term h2k2D to the solution of the difference equation corresponding to the differential equation given by

Eq. (7) as follows

Vm,n =

h2

"0

e

"0 np(r) ne(r)

eV (r)

kBTe . (5)

Assuming that the normalized potential eV (r)/kBTe becomes small in the distance and also assuming that np = n0,

Eq. (4) becomes

r2V (r) =

(r)

"0 +

m,n

4 + h2k2D Wm Wm Wn Wn

,

W = e2i/N, (8)

where Vm,n is the electrostatic potential solution in k space at the grid point (m,n), N is the total number of cells in 2-D

Ann. Geophys., 34, 845855, 2016 www.ann-geophys.net/34/845/2016/

e2n0

"0kBTe V (r) (6)

with the rst-order approximation. Dening kD =

e2n0/"0kBTe, we obtain

(r2 k2D)V (r) =

(r)

K. Hoshi et al.: Thrust calculation of electric solar wind sail 851

0.0 100

1.2 108

(a) Electron density

t = 0 s t = 25 s t = 50 s t = 75 s t = 100 s t = 125 s t = 150 s t = 175 s t = 200 s

1.0 108

Density [m-3 ]

8.0 107

6.0 107

4.0 107

2.0 107

t = 0 s t = 25 s t = 50 s t = 75 s t = 100 s t = 125 s t = 150 s t = 175 s t = 200 s

2.0 107

(b) Proton density

1.5 107

Density [m-3 ]

1.0 107

5.0 106

0.0 100

(c) Space charge density

4.0 10-12

-80 -60 -40 -20 0 20 40 60 80

Charge density [C/m3 ]

0.0 100

-4.0 10-12

-8.0 10-12

t = 0 s t = 25 s t = 50 s t = 75 s t = 100 s t = 125 s t = 150 s t = 175 s t = 200 s

-1.2 10-11

-1.6 10-11

-2.0 10-11

Distance from tether [m]

Figure 7. Density structure along the x axis (y = 76.8 m, z = 76.65 m, V0 = 2.0 kV). (a) Electron density. (b) Proton density. (c) Space

charge density.

space, and h is the cell width used in the Fourier transformation. The 2-D inverse FFT of Eq. (8) was taken to obtain the shielded electrostatic potential in plasma, and the solution was found to be consistent with the potential given by the full PIC simulation.

The proposed estimation method was used to obtain a realistic potential without performing the full PIC simulation.To realize an equipotential within the radius of the tether, the capacity matrix method was also used (Hockney and Eastwood, 1981). Figure 8 compares the potential structures obtained using the 2-D FFT method with the PIC results. A cell width of h = 0.03 and N = 8192 cells in 2-D space were

used for the Fourier transformation. The potential structures estimated using the proposed method were consistent with those obtained from the PIC simulation results (Fig. 5). At V0 = 1 kV, there was a small difference (approximately 20 V)

between the potential structures estimated using the proposed

method and the PIC simulation results behind the tether, but this difference did not cause a difference in the thrust, because only the potential in front of the tether contributes to the thrust.

3.2 Effective radius method (Method 1)

Method 1 considers the effective radius rs that satises the following equation:

eV (rs) =

1

2mpv2d, (9)

where rs is the distance from the tether at which the electro-static potential energy is equal to the kinetic energy of the drifting proton. Janhunen and Sandroos (2007) assumed that the scattering cross section is proportional to rs with a coefcient of proportionality of K, meaning the thrust per unit

www.ann-geophys.net/34/845/2016/ Ann. Geophys., 34, 845855, 2016

852 K. Hoshi et al.: Thrust calculation of electric solar wind sail

0.0

1.0

-80 -60 -40 -20 0 20 40 60 80

2.0

-80 -60 -40 -20 0 20 40 60 80

(a) V0

PIC simulation 2D FFT

= 1.0 kV

(b) V0

PIC simulation 2D FFT

= 2.0 kV

0.8

1.5

Potential [kV]

Potential [kV]

0.6

1.0

0.4

0.2

0.5

0.0

Distance from tether [m]

Distance from tether [m]

3.0

-80 -60 -40 -20 0 20 40 60 80

4.0

-80 -60 -40 -20 0 20 40 60 80

(c) V0

PIC simulation 2D FFT

= 3.0 kV

(d) V0

PIC simulation 2D FFT

= 4.0 kV

3.5

2.5

3.0

2.0

2.5

Potential [kV]

1.5

Potential [kV]

2.0

1.0

1.5

1.0

0.5

0.5

0.0

0.0

Distance from tether [m]

Distance from tether [m]

Figure 8. Comparison of potential structures at surface potentials of (a) V0 = 1 kV, (b) V0 = 2 kV, (c) V0 = 3 kV, and (d) V0 = 4 kV. Blue

lines show the potential structure obtained using the proposed method. Purple lines show the results of the PIC simulations.

length dF/dL can be expressed as

dFdL = KrsPdyn

= Krsn0mpv2d, (10)

where Pdyn = mpn0v2d is the dynamic pressure of a solar

wind proton.

In Janhunen and Sandroos (2007), V (r) is given by Eq. (1). In the present study, V (r) was replaced with the numerical solution obtained using the 2-D FFT method, which considers potential shielding by ambient electrons. For the numerical calculation, the value of r that minimizes the difference between the potential energy and the kinetic energy f1, which is expressed as

f1(r) = eV (r)

1

2mpv2d, (12)

where Lm = mpvd is the angular momentum and is the

impact parameter. Dening Esw = 12mpv2d yields

eV (r)

Esw +

L2m

2mpr2 =

2r2 = 1. (13)

Additionally, l and Ueff were dened as

l

r (14)

Ueff(l)

eV (r)

Esw + l2. (15)

Ueff(l) is equal to the left-hand side of Eq. (13), meaning the distance rmin that minimizes f2(r ) = 1 Ueff(r ) is equiv

alent to the distance rs obtained using the effective radius method. When the scattering angle is dened as [notdef]() =

(), as shown in Fig. 4 of Sanchez-Torres (2014), the

following equations give [notdef]():

1

2mpv2d, (11)

is used as rs instead the value of r. The coefcient of proportionality was set to K = 3.09, as obtained by Janhunen and

Sandroos (2007).

3.3 Coulomb scattering method (Method 2)

The thrust modeling method by Sanchez-Torres (2014) is based on Coulomb scattering. The following equation is formed by adding an angular momentum term to Eq. (9):

eV (r) +

Ann. Geophys., 34, 845855, 2016 www.ann-geophys.net/34/845/2016/

K. Hoshi et al.: Thrust calculation of electric solar wind sail 853

140

120

-1

Estimated thrust per unit length [nN m ]

100

80

60

40

20

0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

Potential [kV]

Figure 9. Comparison of thrust estimations. The purple line shows the modied estimation obtained using the effective radius method. The green line shows the modied estimation based on the Coulomb scattering method. The blue and black lines are the same as those in Fig. 4.

() =

In the present estimation, V (r) in Eq. (15) was replaced with the numerical solution obtained using the 2-D FFT method, and the maximum affection length rmax = 2.0 D

was used instead of the value of rmax = rsheb used by

Sanchez-Torres (2014).

3.4 Estimation results

The thrust per unit length was obtained for V0 = 0 to 4.0 kV

with rw = 15 cm. Figure 9 shows the thrust estimated using

the effective radius method (purple line) and the Coulomb scattering method (green line). These modied thrust estimations are similar to the thrust obtained from the PIC simulation, meaning they are also signicantly lower than estimations from previous studies. The modied estimated thrusts obtained using Methods 1 and 2 are 79.1 and 62.0 nN m1, respectively. The thrust estimated using Method 2 at 4.0 kV is in very good agreement with the PIC result of 67.1 nN m1.

It should be noted that the difference between the modied and previous estimation methods is simply the electrostatic potential structure. The effective proton deection area with shielding is smaller than that without the shielding. Thus, the inclusion of potential shielding results in reduced thrust.

These two modied estimation methods yielded similar estimated values of the thrust per unit length. Thus, the two estimation procedures were shown to be essentially similar, and the difference between the original estimated thrust values was revealed to have been caused by differences in the potential structures. The effective radius method contains an approximation coefcient K, whereas the Coulomb scattering method does not have contain an approximation coef-

www.ann-geophys.net/34/845/2016/ Ann. Geophys., 34, 845855, 2016

8 <

:

rmax

[integraldisplay]

rmin

1

[integraldisplay]

rmax

dr

dr

r2

(16)

r2

q1 2r2 eV (r)Esw+

q1 2r2

9 =

;

= 1() + 2() (17)

1() =

lmax

[integraldisplay]

(18)

lmin

dlp1 Ueff(l)

2() =

lmin

[integraldisplay]0

dlp1 l2 =

sin(lmin), (19)

where lmax = /rmin and lmin = /rmax. rmax denotes the

maximum length that scattered particles can be affected by an electrostatic potential. rmax was dened as rsheb by

Sanchez-Torres (2014) and has been calculated as follows by Sanmartn et al. (2008):

1.53"1 2.56

D rsh

rshrw

[parenrightbigg] =

eV0kBTe . (20)

Then, the thrust per unit length can be expressed as

dFdL = n0vd 2

rmax

[integraldisplay]0

4/5

[bracketrightBigg][parenleftbigg]

rsh

D

4/3ln

2mpvd sin

[notdef]() 2

cos

[notdef]() 2

d (21)

= 4Pdyn

rmax

[integraldisplay]0

sin2

[notdef]() 2

d. (22)

854 K. Hoshi et al.: Thrust calculation of electric solar wind sail

Table 2. Differences between the present PIC simulation and those by Janhunen and Sandroos (2007).

Present simulation 1-D by Janhunen and Sandroos (2007) 2-D by Janhunen and Sandroos (2007)

Method 3-D electrostatic PIC 1-D electrostatic PIC 2-D electrostatic PIC Plasma species electron, proton electron (proton density is constant) electron, proton, helium ions Plasma density n0 = 1.0 10

7 m3 n0 = 0.73 10

7 m3 n0 = 0.73 10

7 m3

Plasma temperature Te = 100 eV, Tp = 12 eV Te = Tp = 6 to 24 eV Te = Tp =6 to 24 eV

Cell width 0.3 m 2.0 m 1.25 mTime step 20 ns not described 15.625 nsSimulation duration 0.2 ms not described 40 to 100 msSystem length 152 m (5123 cells) 100 m (50 cells) 320 m 160 m (256 128 cells)

Tether radius 0.15 m 1.0 m 1.25 mParticles per cell 60 (30 of each species) not described 40 (20 of each species)

cient; thus, the Coulomb scattering method was considered to be more consistent with the PIC results.

4 Discussion

In this study, a full PIC simulation of the E-sail was conducted. The most signicant difference between the results obtained in this paper and those obtained in previous studies is that the electrostatic potential structure around the tether was considered, which yielded different values of the thrust. In the PIC simulation by Janhunen and Sandroos (2007), the potential shielding effect of the ambient electrons was not signicant. In the analytical estimation by Sanchez-Torres (2014), they assumed the absence of ambient electrons so their estimations of the potential (Eq. 3) were higher than those obtained using Eq. (1) at high positive potentials. In the present PIC simulation, the ambient electrons shield the potential of the tether more effectively than those in Janhunen and Sandroos (2007) and Sanchez-Torres (2014).

This study focused on the source of the differences between the present PIC simulation and the PIC simulations performed by Janhunen and Sandroos (2007), which are shown in Table 2.

There are several differences among the simulations, particularly regarding the tether dimensions, plasma conditions, simulation duration, and cell width. We consider that the differences caused the disagreement of the simulation result. Note that the electron temperature kBTe = 100 eV adopted

in this paper is the typical value in solar wind at 0.5 AU and is about 1 order of magnitude higher than the typical value at 1 AU (applied in previous studies).

To the best of the authors knowledge, this study is the rst to perform a full 3-D PIC simulation of the E-sail without any approximations. The radius of the tether (rw = 15 cm)

was large in comparison with several tens of micrometers; thus, the E-sail was not completely simulated, but the difference is negligible because 15 cm is much smaller than D (approximately several meters to a few tens of meters in in-

terplanetary space), so the ambient electron collection is not signicantly different.

The increase in thrust caused by the removal of electrons, which was discussed by Janhunen (2009), was not investigated in this study; the present simulation did not consider multiple tethers, and the simulation duration (0.2 ms) was not long enough to describe such a effect. If any efcient trapped electron removal mechanisms exist, the thrust of E-sail may increase asymptotically, so we must remark that the thrust characteristics obtained in this paper do not consider the effect of the trapped electron removal. However, the presented simulation successfully described the response of ambient electrons as shown in Fig. 3. Hence, our results reveal at least the minimum thrust characteristics of the E-sail.

A modied thrust estimation proposed in this paper, which is obtained by replacing the electrostatic potential structures used in the estimations in the previous studies, is a better reference model of the minimum thrust of the E-sail. The proposed estimation method can be easily used to calculate the minimum thrust of E-sail.

5 Conclusions

In this study, the rst full 3-D PIC simulation of the tether of the E-sail was performed, and the transient of its thrust was numerically calculated. At V0 4.0 kV, the thrust ob

tained by PIC simulation was almost half of the thrust estimated in previous studies. This difference is caused by the electrostatic potential structure around the tether. The potential structure in the present simulation differed greatly from the structure used in previous estimations due to the strong potential shielding by ambient electrons.

Additionally, a modied thrust estimation method with a shielded potential structure was proposed. In this new method, the potential structure employed in previous estimations was replaced with the potential structure derived from our simulation result. The estimated thrust obtained using the modied method agreed very well with the PIC simulation

Ann. Geophys., 34, 845855, 2016 www.ann-geophys.net/34/845/2016/

K. Hoshi et al.: Thrust calculation of electric solar wind sail 855

results. The proposed method can be easily used to calculate the minimum thrust of E-sail.

In future work, we will perform the long-duration simulation and investigate an asymptotic thrust characteristics. We also plan the PIC simulation of a much thinner tether using various simulation techniques, such as the ctitious surface method.

6 Data availability

The PIC simulation data (approximately 1TB), which include potential structure, electric eld, and density structure, are available upon requests.

The Supplement related to this article is available online at http://dx.doi.org/10.5194/angeo-34-845-2016-supplement

Web End =doi:10.5194/angeo-34-845-2016-supplement .

Author contributions. K. Hoshi improved the performance of the PIC code, performed the PIC calculation, analyzed the results, wrote the thrust estimation code, and wrote the paper. H. Kojima and H. Yamakawa directed the study and discussed the interpretation of the results. T. Muranaka rst developed the full PIC code and helped run the code on the supercomputer system.

Acknowledgements. The computations in the present study were performed using the Kyoto-daigaku Denpa-kagaku Keisanki-jikken (KDK) system at the Research Institute for Sustainable Humanosphere (RISH) at Kyoto University. The present study was supported by JSPS KAKENHI Grant-in-Aid for JSPS Fellows Number 15J08941.

The topical editor, E. Roussos, thanks T. Laeur and one anonymous referee for help in evaluating this paper.

References

Hockney, R. W. and Eastwood, J. W.: Computer Simulation Using

Particles, McGraw-Hill International Book Co., New York, 1981. Hoshi, K., Kojima, H., and Yamakawa, H.: Numerical Analysis of Potential Structure around Electric Solar Wind Sail Tether, Transactions of JSASS Aerospace Technology, in press, 2016. Janhunen, P.: Electric sail for spacecraft propulsion, J. Propul.

Power, 20, 763764, 2004.

Janhunen, P.: Increased electric sail thrust through removal of trapped shielding electrons by orbit chaotisation due to spacecraft body, Ann. Geophys., 27, 30893100, doi:http://dx.doi.org/10.5194/angeo-27-3089-2009

Web End =10.5194/angeo- http://dx.doi.org/10.5194/angeo-27-3089-2009

Web End =27-3089-2009 , 2009.

Janhunen, P. and Sandroos, A.: Simulation study of solar wind push on a charged wire: basis of solar wind electric sail propulsion, Ann. Geophys., 25, 755767, doi:http://dx.doi.org/10.5194/angeo-25-755-2007

Web End =10.5194/angeo-25-755-2007 http://dx.doi.org/10.5194/angeo-25-755-2007

Web End = , 2007.

Muranaka, T., Shinohara, I., Funaki, I., Kajimura, Y., Nakano, M., and Tasaki, R.: Research and development of plasma simulation tools in JEDI/JAXA, Journal of Space Technology and Science, 25, 118, 2011.

Sanchez-Torres, A.: Propulsive Force in an Electric Solar Sail, Contrib. Plasma Phys., 54, 314319, 2014.

Sanmartn, J. R., Choinire, E., Gilchrist, B. E., Ferry, J., and Martnez-snchez, M.: Bare-Tether Sheath and Current: Comparison of Asymptotic Theory and Kinetic Simulations in Stationary Plasma, IEEE Transactions of Plasma Science, 36, 28512858, 2008.

www.ann-geophys.net/34/845/2016/ Ann. Geophys., 34, 845855, 2016