Orbital maneuvering of electric solar wind sail based on an advanced solar wind force model

Journal Pre-proof Orbital Maneuvering of Electric Solar Wind Sail Based on an Advanced Solar Wind Force Model Kohei Yamaguchi, Kikuko Miyata PII:

S0094-5765(19)31307-4

DOI:

https://doi.org/10.1016/j.actaastro.2019.10.001

Reference:

AA 7686

To appear in:

Acta Astronautica

Received Date: 2 August 2019 Revised Date:

24 September 2019

Accepted Date: 2 October 2019

Please cite this article as: K. Yamaguchi, K. Miyata, Orbital Maneuvering of Electric Solar Wind Sail Based on an Advanced Solar Wind Force Model, Acta Astronautica, https://doi.org/10.1016/ j.actaastro.2019.10.001. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 IAA. Published by Elsevier Ltd. All rights reserved.

Orbital Maneuvering of Electric Solar Wind Sail Based on an Advanced Solar Wind Force Model Kohei Yamaguchia,1,∗, Kikuko Miyataa,2 a

Nagoya University, Furo-cho, Chikusa-ku, Nagoya-city, Aichi 464-8603, Japan

Abstract The electric solar wind sail is a propulsion system that extracts the solar wind momentum for the thrust force of a spacecraft by using an interaction between solar wind protons and the electric potential structure around charged long thin conducting tethers. The system enables a spacecraft to generate a thrust force without consuming reaction mass. This paper investigates the capability of the electric solar wind sail as a propulsion system for deep space exploration missions. The shape of the conducting tether that is determined by the equilibrium of the solar wind force and centrifugal force is numerically calculated for formulating an advanced solar wind force model. The conducting tethers deviate from the ideal sail spin plane, and the maximum value of the thrust direction varies from 13◦ to 19◦ . To estimate the spacecraft thrust vector, which is calculated as the sum of solar wind force vectors exerted on each tether, best-fit polynomial equations are proposed. We performed numerical simulations for a two-dimensional orbital transfer mission to investigate the capability of the electric solar wind sail. Results of numerical simulations show that the electric solar wind sail enables spacecraft to perform Earth–Venus, Earth–Mars, and Earth–Itokawa transfer missions. Additionally, this paper performs three-dimensional simulations for an Earth–Ryugu transfer mission. The electric solar wind sail ∗

Corresponding author Email addresses: [email protected] (Kohei Yamaguchi), [email protected] (Kikuko Miyata) 1 Assistant Professor, Graduate School of Engineering, Department of Aerospace Engineering (+8152-789-3287) 2 Assistant Professor, Graduate School of Engineering, Department of Aerospace Engineering (+8152-789-3289)

Preprint submitted to Acta Astronautica

October 10, 2019

achieves a more complicated orbital transfer in a reasonable mission time. Keywords: electric solar wind sail, propellantless propulsion, trajectory optimization 1. Introduction The electric solar wind sail (E-sail) is a propulsion system that extracts the solar wind momentum for a spacecraft propulsion force. The fundamental concept of the E-sail was devised in 2004 by Janhunen [1] and has been studied by many researchers. The basic idea consists of spinning the spacecraft main body around a spin axis. A number of conducting tethers are deployed by the rotational motion. The main tethers are connected to the auxiliary tether to keep the tether rig dynamically stable [2]. Remote units are placed at the tips of each tether to deploy the auxiliary tether and to control the angular velocity of the main tethers. Conducting tethers are held at a high positive electric potential through the loss of electrons caused by the work of electron guns mounted on the main body. The electron guns are usually powered by solar panels [3, 4]. The specific acceleration of the Esail-based spacecraft is estimated to be sufficient for deep space exploration, and even up to 10 mm/s2 [2]. A method to estimate the electric potential structure around a conducting tether was proposed based on the results of three-dimensional particle-in-cell (3D PIC) simulations [5]. Researchers also proposed a method to estimate the characteristics of the E-sail thrust force based on a 3D PIC simulation [6]. Additionally, a feasible mass budget model of the E-sail was proposed [7]. A tether manufacturing technique was also developed. In 2013, automatic production of a 1 km tether composed of 25 and 50 µm diameter wires was demonstrated [8]. A satellite called ESTCube-1 was launched in 2013 to deploy a 10 m tether and to measure the thrust force [9]. Although the experiment failed, it was still significant that this was the first demonstration of the E-sail in space. An experiment to deploy a 100 m tether and measure the electrostatic drag force was planned as Aalto-1 [10, 11]. The E-sail uses the natural solar wind plasma stream to generate a thrust force, and it can enable a spacecraft to perform various missions efficiently. The performance of the E-sail was numerically investigated in minimum-time orbital transfer scenarios to Mars and Venus [12]. The E-sail was also studied to bring a spacecraft out of the solar system [13]. In addition, the E-sail 2

was studied as a propulsion system for several kinds of asteroidal missions. These were rendezvous to potentially hazardous asteroids [14], sample return missions from an asteroid [15], and flyby missions [16]. Even in a planetary defense field, the E-sail is an attractive system for kinetic impactors [17] and gravity tractors [18] to deflect an asteroid with an unacceptable Earth collision risk. Moreover, since a feasibility study showed that the E-sail can accelerate a spacecraft effectively, the NASA Marshall Space Flight Center studied the E-sail for a rapid heliopause mission3 . To investigate the capability of the E-sail, a solar wind force model was developed. In the first force model, the corresponding thrust angle was assumed to be equal to half the sail plane’s inclination angle. In addition, the thrust acceleration magnitude was assumed to be unaffected by the sail plane’s inclination angle [12]. An advanced solar wind force model in which the thrust angle and magnitude were expressed as a form of the polynomial function of the sail inclination was proposed [19]. Even assuming this advanced force model, the E-sail still can conduct missions such as accelerating an impactor [17] and minimum-time circle–to–circle transfers [19, 20]. In Ref. [21], the E-sail force model was refined with an analytical/geometrical approach. We also proposed an orbital maneuvering technique that fixes the E-sail attitude to an inertial frame. Unlike the traditional thrust vectoring technique that actively changes the E-sail attitude, the tether voltage is only changed at appropriate timings [17, 19]. In this paper, we incorporate the effects of changes in the shapes of the conducting tethers into the solar wind force model of the E-sail. The relationship between the tether shape and E-sail force vector was studied. Toivanen and Janhunen introduced a coning angle that defines the deviation of a straight E-sail tether from an ideal sail plane [22]. They found that the deviation of tethers affects the E-sail force vector even when assuming straight tethers. They also showed a numerical method to calculate the shape of warped tethers by separating the tether into many uniform segments and introducing local coning angles for each segment [23]. We expand the proposed method and formulate a solar wind force model by summing force vectors exerted on all segments of the conducting tethers. For sun-facing E-sail, 3

Heliopause Electrostatic Rapid Transit System http://www.nasa.gov/content/heliopause-electrostatic-rapid-transit-systemherts/#.VTsZdpOo1c4 (Accessed on Aug. 7th, 2017)

3

(HERTS),

Bassetto et al. investigated the thrust force and torque acting on the E-sail tether based on the tether shape analysis [24]. They described the shape of the tether through continuously differentiable functions and formulated analytical expressions of the total force and torque. For near-sun-facing Esail, in which the sail plane’s inclination is sufficiently small, they also derive thrust and torque vectors [25]. However, the orbital maneuvering technique we investigate in this paper needs to determine the thrust force for larger sail plane’s inclination. In this paper, we calculate the shape of E-sail tethers numerically for the range of the inclination angle 0◦ ≤ α ≤ 90◦ to determine the thrust force of the E-sail. Instead of the analytical expression, we introduce a polynomial form that provides a good approximation for the numerically calculated result. This paper is organized as follows. Section 2 provides a method to calculate the shapes of charged conducting tethers in three-dimensional space. A coordinate system that is fixed to a given conducting tether is defined, and equations of the local coning angle are numerically integrated. In addition, some examples of the E-sail shape are provided by varying the ideal E-sail spin plane’s inclination angle. A new solar wind force model is developed in Sec. 3. Since the E-sail tether shape and force depend on the tether voltage and solar wind proton density, the direction and magnitude of the thrust force are expressed as a function of the sail plane’s inclination angle, thrust modulation coefficient, and distance between the spacecraft and the sun. The radial, transverse, and out-of-plane thrust components are also provided. In Sec. 4, we investigate the capability of the E-sail-based spacecraft for orbital transfer missions by exploiting the proposed force model. Two-dimensional trajectories for rendezvous missions from Earth to Venus, Mars, and Itokawa are calculated in an optimal framework. In addition, as a more challenging case, a three-dimensional orbital transfer mission targeting an inclined orbit is studied. Section 5 concludes this paper with a brief discussion. 2. Tether Shape Analysis 2.1. Tether shape in 2D coordinate system In this paper, the attitude of the E-sail is determined in the coordinate systems defined in Ref. [23], as shown in Fig. 1. The reference frame moves together with the spacecraft center of mass. Since the E-sail uses the solar wind coming from the sun, it is useful to define the attitude of the E-sail in a rotating reference frame such that one axis always points to the sun. This 4

rotating reference frame is defined by (x∗ , y ∗, z ∗ ), where the z ∗ -axis points to the sun. The coordinate system (x, y, z), which shares its origin with (x∗ , y ∗, z ∗ ), is fixed to the centroid of the spacecraft main body. The z-axis is a spin axis of the main body, and the angle between the z ∗ - and z-axes is defined as the sail plane’s inclination angle α. As can be seen in Fig. 1, the solar wind is perpendicular to the y-axis. Note that the y- and y ∗-axes point in the same direction. If there is no solar wind force, then the tethers are in the xy-plane. This is defined as an ideal sail plane. Consider the ρ-axis that is in the xy-plane and crosses the z-axis at the origin of the (x, y, z) frame. If we assume that an E-sail tether remains in the ρz-plane, forces exerted on a tether are defined as shown in Fig. 2. A single tether shape is expressed not with a straight line but by a smooth warped line consisting of many short segments. A local coning angle λ, which expresses the angle between the ρ-axis and a tether segment, is introduced. Unlike the coning angle, which defines the angle between a straight tether and an ideal sail plane, the angle λ is useful in modeling a more realistic warped tether shape. Since the solar wind force is parallel to the component of the solar wind velocity perpendicular to the tether, v ⊥ , the force exerted on a single tether, F , can be written in an integration form as Z ρL √ (1) F = ξv ⊥ 1 + u2 dρ ρ

where ρL is the position of the remote unit on ρ-axis and u = cos λ. The E-sail force factor ξ is the magnitude of the solar wind force per tether length [26]. Equation (1) can be separated into two components as follows: Z ρL (cα − sα u)u √ Fρ = −ξvsw dρ, (2) 1 + u2 ρ Z ρL cα − sα u √ Fz = ξvsw dρ (3) 1 + u2 ρ where cα = cos α and sα = sin α. In addition, the total centrifugal force G is given by Z ρL √ 2 G = µω ρ 1 + u2 dρ + mru ω 2ρL (4) ρ

where µ is the linear mass density of the tether, ω is the E-sail angular velocity, and mru is the mass of the remote unit placed at the tip of the 5

tether. Note that all integrals in Eqs. (2)–(4) can be rearranged as forms of summation by assuming that the problem can be locally iscretized. Local values of Fρ , Fz , and G for the i-th tether segment are given by Z ρL (cα − sα ui−1)ui−1 i p dρ + Fρi−1 , (5) Fρ = −ξvsw i−1 2 i−1 1 + (u ) ρ Z ρL c − sα ui−1 i pα Fz = −ξvsw dρ + Fzi−1 , (6) i−1 2 i−1 1 + (u ) ρ Z ρL p i 2 G = µω ρ 1 + (ui−1 )2 dρ + Gi−1 . (7) ρi−1

In addition, an equation for a local tether segment can be written as u=

Fz . G + Fρ

(8)

Substituting Eqs. (2)–(4) into Eq. (8), and after some algebra, the local value of u for the i th tether segment ui is given by ξvsw cα ∆L + Fzi−1 u = (ξvsw sα + µω 2 ρi−1 )∆L + Gi−1 + mru ω 2 ρL − Fρi−1 i

(9)

where ∆L is the unit length of the tether. To find the shape of the tether, Eq. (9) is numerically integrated from a tip of the tether (ρ = ρL ) to the tether attachment point on the spacecraft main body (ρ = ρ0 ). The integration is (0) begun from an initial value of ρL = L over the tether length. When the (i−1) result of the i − 1-th iteration, ρL , is known, the initial value of ρL at the (i) (i−1) (i−1) i-th iteration is updated as ρL = ρL − ρ0 . This process is iteratively (i) conducted until the condition |ρ0 | < ǫ is satisfied. 2.2. Tether shape in 3D coordinate system To consider the 3D shape of a single conducting tether, a coordinate system (ρ, ψ, z) is introduced by rotating the body-fixed frame (x, y, z) around the rotation axis z by the angle β, as shown in Fig. 3. Note that the two coordinate systems share the z-axis. If the angle α is not 0◦ , then the solar wind force pushes the tether in the ψ direction, and the tether deviates from the ρz-plane. Note that the centrifugal force G is assumed to be parallel to the ρ-axis. Here, a new coordinate system (ρ, ψ ′ , z ′ ) is introduced to calculate 6

the conducting tether shape in the ρz ′ -plane. Figure 4 shows the relationship between the two coordinate systems (ρ, ψ, z) and (ρ, ψ ′ , z ′ ). By rotating the (ρ, ψ, z) coordinate system around the ρ-axis to make the ρz ′ -plane parallel to the solar wind velocity vector, the tether shape can be calculated in this plane. The rotation angle around the ρ-axis is defined as φroll . In the ρz ′ plane, the local coning angle λ can also be introduced. In addition, the solar wind coning angle in the ρz ′ -plane, αρz ′ , is defined. Now, we can calculate the shape of a given E-sail tether by using the same scheme proposed in Ref. [23]. As Fig. 4 shows, the solar wind direction in the ρz ′ -plane, αρz , can be calculated with some trigonometric identities as αρz ′ = sin−1 (sin α cos β).

(10)

Similarly, the rotation angle φroll is given by cos α

−1

φroll = cos

p

1 − sin2 α cos2 β

!

.

(11)

To avoid numerical issues that arise when β = 0◦ , a little ingenuity is required for Eq. (11). After calculating the tether shape in the ρz ′ -plane, it is rotated in order to visualize the shape of each tether in the (x, y, z) coordinate system. In addition, by summing all force vectors exerted on all tethers, the thrust force vector of the E-sail can be investigated. 2.3. Numerical solution The shape of the E-sail tethers is calculated with the proposed method. Table 1 lists the E-sail parameters assumed here. We assume an E-sail contains 100 tethers (N = 100) with a length L of 20 km each. By using a constant solar wind bulk speed vsw = 400 km/s, the solar wind force exerted on a unit length of the tether ξvsw = 0.4724 mN/km is determined. The rotation period ∆t is assumed to be 125 min [23]. Figure 5 shows three examples of the E-sail tether shape. Figure 5a is the result when α = 0◦ , b is when α = 45◦ , and c is when α = 90◦ . Relations between the normalized thrust acceleration vector asw /a0 and vector n, which is normal to the ideal spin plane, are also indicated. As can be seen in Fig. 5a, the solar wind is exerted uniformly on the tethers, and the tethers deviate from the ideal spin plane (xy-plane) by more than 5 km. This means that if we take the actual tether shape into account, the thrust components parallel to the xy-plane 7

Table 1: Simulation condition.

Items L N vsw ξ Vmax mru ∆t (ω

Values 20 km 100 400 km/s 1.176 × 10−9 mN · s · m−2 20 kV 1 kg 125 min [23] 2.88 deg/min)

are generated. Such components are cancelled out because of the symmetric structure of the E-sail, and the total E-sail thrust acceleration magnitude is decreased compared with the model with which the tethers lie on the ideal sail spin plane. In this case, the magnitude of the thrust force is 96.96% of the previous model. The equilibrium shape for α = 0◦ can also be found by using the analytical model in Ref. [24] if ω is sufficiently large. When the sail plane’s inclination angle is increased, the deviation becomes small, as shown in Fig. 5b, and becomes zero when α = 90◦ is satisfied, as shown in Fig. 5c. Figure 5 also indicates that the characteristics of the thrust vector are affected by changes in the tether shape. For example, if all tethers are in the ideal sail plane, unlike Fig. 5c shows, then the value of |asw /a0 | is 0.5 [19]. That is, the direction and thrust acceleration magnitude vary depending on α. Variations of the thrust direction αT and magnitude γ with the sail plane’s inclination angle α are investigated by varying α in steps of 0.9◦ (= π/100). Here, results for three cases of the rotation period, ∆t = 100 min, ∆t = 125 min, and ∆t = 150 min are calculated. Comparisons of these results with that of the polynomial model proposed in the previous work [19] are plotted in Fig. 6. Note that the dashed lines in Fig. 6 represent the model proposed in the previous work. The proposed model indicates a similar tendency as the previous model. Since the tether shape approaches flat when the spin rate is large, the polynomial model agrees better with the ∆t = 100 min case. When the spin rate is small, as the case for ∆t = 150 min, the result shows a larger deviation from the result of the polynomial model. In this paper, the case for ∆t = 125 min is used for the latter numerical simulations. Figure 6a shows that the maximum value of the thrust direction of ∆t = 125 min case 8

decreases by about 2◦ from the polynomial model. Figure 6b shows that the normalized thrust acceleration magnitude γ is larger than that of the previous model at α = 90◦ because the acceleration is divided by the value when α = 0◦ . Figure 6 shows that the polynomial model can capture the qualitative features of the thrust direction and magnitude. Besides, it also implies that we need a better model to deal with the deviations caused by the warped tether shape. 3. Simulation Model 3.1. Modeling of thrust direction and magnitude of E-sail force As described in the preceding section, conducting tethers of the E-sail are warped because of forces exerted on the tethers. The shape of the tethers depends on the tether voltage, rotation period that defines the centrifugal force, and solar wind parameters. Here, we assume that the rotation period is maintained at a given constant value through the work of remote units. Factors that affect the tether shape are contained in the E-sail force factor ξ. We introduce a normalized force factor ξ/ξ0. Note that we assume the reference force factor value is ξ0 = 1.176 × 10−9 mN · s · m−2 . The relation between the sail plane’s inclination angle α, value of ξ, and thrust direction αT is numerically calculated by varying α in steps of 0.9◦ and ξ/ξ0 of 0.04. Figure 7 shows the calculation result. In Fig. 7, the black continuous lines are placed at 5◦ increments, and gray continuous lines are placed at 1◦ increments. Figure 7 shows that larger values of ξ/ξ0 decrease the maximum value of the thrust direction angle αT . When the value of ξ/ξ0 increases, the local coning angles also increase because of the increment of the solar wind force. On the other hand, when the value of ξ/ξ0 decreases, the tethers experience a smaller solar wind force, and their deviation from the ideal plane becomes small. Therefore, the value of the thrust direction angle approaches the value of the previous force model when ξ/ξ0 is small. If we assume that the minimum solar distance is 0.3 au, then the maximum value of the force factor ξ/ξ0 is 3.3. This means that the thrust vectoring capability of the E-sail αmax varies from 13◦ to 19◦ depending on the force factor. In Fig. 8, the normalized thrust acceleration magnitude γ is also shown as a function of the angle α and force factor ξ/ξ0. The thrust force is divided by the value of the force at α = 0◦ . The contour lines in Fig. 8 also warp as the value of ξ/ξ0 changes. Given the thrust force direction shown in Fig. 7, the magnitude γ approaches that of the previous force model when ξ/ξ0 becomes small. As 9

described, a larger value of ξ/ξ0 increases the local coning angle λ. When λ becomes large, thrust components that are cancelled out at α = 0◦ are also increased. Since the value of γ is normalized by the acceleration magnitude at α = 0◦ , γ increases when ξ/ξ0 is large. To improve our previous solar wind force model [19], the dependency of the E-sail force on the force factor ξ is taken into consideration. The E-sail force factor ξ is defined as follows [26]: √ ξ = 0.18max[0, V0 − Vsw ] ǫ0 np mp (12) where ǫ0 is the conductivity of vacuum, np is the number density of proton, and mp is the mass of proton, respectively. Since the voltage Vsw corresponding to the bulk kinetic energy of protons is relatively small (Vsw ∼ 1 kV) [2], it is usually neglected. Here, to introduce the concept of artificial tether voltage modulation, the net tether voltage V = V0 − Vsw is defined. To simplify the subsequent analysis, the voltage V is expressed with a switching parameter τ and maximum tether voltage Vmax as V = τ Vmax .

(13)

The range of the switching parameter is defined by 0 ≤ τ ≤ 1. Since the number density of protons np in Eq. (12) is proportional to the inverse square of the distance between the Sun and the spacecraft, r [27, 28], the E-sail thrust force is proportional to 1/r. The force factor ξ is proportional to a component constructed by such parameters, and it is ξ∝

τ . r

(14)

This means that the value of the E-sail force factor can be calculated by using the switching parameter τ and the distance from the sun r at a given time. We define this ratio τ /r as the force distance ratio, and it satisfies the following: ξ τ = . (15) ξ0 r We define the following 5th-order polynomial equation for the thrust direction αT to model characteristics of the E-sail force shown by numerical results: 5 X 5  τ X  τ j αT α, = Pi,j αi . (16) r r i=0 j=0 10

Table 2: Best-fit interpolation coefficients Pi,j for thrust direction αT . PP j 0 1 2 3 4 PP i P −2 0 1 2 3 4 5

−8.421 × 10 0.4858 1.065 × 10−3 −3.516 × 10−5 −6.077 × 10−7 2.081 × 10−9

0.9775 −5.251 × 10−2 −1.305 × 10−3 5.079 × 10−5 −3.274 × 10−7 0

−1.415 −1.260 × 10−2 −8.516 × 10−5 2.755 × 10−6 0 0

0.7404 4.476 × 10−3 −6.225 × 10−5 0 0 0

−0.1598 1.853 × 10−4 0 0 0 0

5

1.137 × 10−2 0 0 0 0 0

Table 3: Best-fit interpolation coefficients Qi,j for normalized thrust acceleration magnitude γ. PP j 0 1 2 3 4 5 PP i P 0 1 2 3 4 5

1.001 5.765 × 10−4 −1.718 × 10−4 1.992 × 10−6 −2.539 × 10−8 1.806 × 10−10

−1.224 × 10−2 −3.265 × 10−4 3.530 × 10−5 −2.999 × 10−7 −8.559 × 10−10 0

2.139 × 10−2 4.334 × 10−4 −7.470 × 10−7 3.602 × 10−8 0 0

1.447 × 10−2 −1.819 × 10−4 −2.589 × 10−7 0 0 0

4.093 × 10−4 2.058 × 10−5 0 0 0 0

−4.040 × 10−4 0 0 0 0 0

The polynomial expression for γ is defined as 5 X 5  τ j  τ X = Qi,j αi . γ α, r r i=0 j=0

(17)

By using a best-fit interpolation of numerical data shown in Figs. 7 and 8, coefficients Pi,j and Qi,j are found. Tables 2 and 3 list these coefficients. Note that the angles α and αT are assumed in degrees. Equation (16) provides an approximation of the numerical result with absolute errors smaller than 1◦ . The absolute errors in γ are at most 0.014. Hence, the polynomial expressions provide good approximations of the numerical results. As mentioned above, since the rotation period determines the E-sail tether shape, it also affects the force model. The thrust direction and magnitude as a function of sail plane’s inclination angle and E-sail force factor are calculated for the cases of ∆t = 75 min, ∆t = 100 min, and ∆t = 150 min and shown in Appendix A. Besides, by using a best-fit interpolation of numerical data, coefficients Pi,j and Qi,j are found for these three cases. 3.2. Equations of motion A solar wind force model is proposed here. To define the E-sail thrust ˆ are introduced as shown in Fig. 9. The unit vector, unit vectors rˆ and p vector rˆ points in the direction of sun–spacecraft line. Using the spacecraft ˆ is along the rˆ × v direction. By using heliocentric velocity vector v, p 11

the reference frame in Fig. 9, the radial, transverse, and out-of-plane thrust components (S, T, W ) are given by  τ 1  τ S =τ γ α, a0 cos αT α, , (18) r r r  τ 1  τ T =τ γ α, a0 sin αT α, sin δ, (19) r r r  τ  τ 1 W =τ γ α, a0 sin αT α, cos δ (20) r r r

where the characteristic acceleration a0 is the E-sail acceleration when τ = 1, r = 1 au, α = 0◦ , and δ = 90◦ . The E-sail thrust direction and magnitude are defined by Eqs. (16) and (17). By introducing these components into heliocentric equations of motion, an improved numerical simulation model for the E-sail is formulated. Let the reference frame (r, θ, ϕ) be a heliocentric inertial spherical polar coordinate system, as shown in Fig. 10. The equations of motion in this system are given by 1 µ⊙ v˙ r = (vϕ 2 + vθ 2 ) + 2 + ar , r r 1 v˙ θ = (vθ vϕ tan ϕ − vr vθ ) + aθ , r 1 v˙ ϕ = − (vr vϕ + vθ2 tan ϕ) + aϕ r

(21) (22) (23)

where (vr , vθ , vϕ ) and (ar , aθ , aϕ ) are velocities and accelerations in the reference frame respectively and µ⊙ is the heliocentric gravitational constant. 3.3. Local switching laws To control the orbit of the E-sail, thrust vectoring by using an active modulation of the sail plane’s inclination angle was studied. Without any control, the sail angle slowly (∼1 deg/day) changes because the sail spin plane maintains its orientation. Toivanen and Janhunen showed that the appropriate tether voltage modulation achieves the desired change in the sail plane’s inclination [22]. For the orbital maneuvering method using the active control of the tether voltage and E-sail attitude, locally optimal steering laws that maximize the change in a given orbital element was formulated [29]. They were demonstrated for some cases including a complex mission. For a near-sun-facing E-sail, a voltage modulation law to deal with an undamped oscillatory motion was developed [29]. We proposed an orbital maneuvering 12

technique in which the ideal sail plane of the E-sail is not actively changed and is fixed to the inertial frame [19]. Instead of changing the E-sail attitude, the voltage of the tethers is changed only at appropriate timings. In this paper, the latter method is employed to control the E-sail orbit. Local switching laws (LSLs) that maximize or minimize the change in a given orbital element for a small change in the true anomaly were proposed and studied. Lagrange’s variational equations that account for the spacecraft semimajor axis a, eccentricity e, and inclination i are given by   da l 2lr 2 Se sin f + T , (24) = df r µ⊙ (1 − e2 )2  r2 n r r o de cos f + T e , = S sin f + T 1 + (25) df µ⊙ l l di r3 = cos (f + ω)W (26) df µ⊙ l where f is the spacecraft true anomaly and l is the semilatus rectum, respectively. Taking an inner product of thrust component vector (S, T, W ) and an appropriate thrust vector, the value of the switching parameter τ is determined based on the sign of the inner product. Appropriate thrust vectors are defined by Eqs. (24)–(26) as follows:  : for semimajor axis  (e sin f, rl , 0)
r r (sin f, 1 + l cos f + e l , 0) : for eccentricity (27)  (0, 0, cos (f + ω)) : for inclination To maximize the change in a given orbital element, the voltage is turned on when the inner product is positive. On the other hand, to minimize the change, the voltage is turned on when the inner product is negative. Since the proposed force model is expressed as polynomial functions, estimating the achievable change rate in the orbital elements by using LSLs is difficult. Here, we numerically integrate the equations of motion to evaluate the effectiveness of the LSLs. Figure 11 shows an example of a locally optimal trajectory of an E-sail using the LSLs to increase the semimajor axis. Note that the characteristic spacecraft acceleration a0 = 0.5 mm/s2 and initial sail attitude (α0 , δ0 ) = (0◦ , 90◦ ) are assumed. Since the attitude of the E-sail is fixed in the inertial reference frame, the inner product is positive when the spacecraft is in the second and fourth quadrants. The semimajor axis is gradually increased and reaches 1.9 au with a flight time of 10 years. Figure 12 13

lists the time histories of changes in the semimajor axis and eccentricity. Figure 12a assumes the LSL for increasing the semimajor axis, Fig. 12b is for a decreasing semimajor axis, and Fig. 12c is for increasing eccentricity. Black continuous lines are calculated assuming the force model proposed in this paper and the dashed lines are from the old force model [19]. We use the characteristic acceleration of a0 = 5.156 mm/s2 to calculate the trajectory of the E-sail with the previous force model because the degradation of the thrust force caused by the warped tether shape does not occur with the model. As the figure shows, a warped tether shape decreases the change rates of the orbital elements achieved by LSLs. Compared with the previous model, the increment in the semimajor axis was 90.51%, the decrement was 88.32%, and the increment in eccentricity was 97.07%. 4. Trajectory Optimization In this section, trajectories of the spacecraft are optimized to investigate the capability of the E-sail as a propulsion system for deep space exploration missions. 4.1. Mission planning Orbital transfer missions for the E-sail are explained here. We assume an E-sail with the characteristic acceleration of a0 = 0.5 mm/s2 . The attitude of the E-sail is fixed to the inertial frame, and the orbit is controlled by modulating the tether voltage. At the E-sail’s departure from Earth, it is assumed that the E-sail is in a 1 au circular heliocentric orbit. The flight time is separated into 100 uniform segments, and the switching parameters τi (i = 1, · · · , 100) are defined for each segment. The optimization variables are the flight time tF , values of τi , initial attitude α0 , and δ0 . Taking into account actual ephemeris constraints enables us to solve a problem for a given realistic transfer mission. Such constraints, on the other hand, derive local solutions and can make it difficult to investigate potentialities of the E-sail. Hence, to optimize the E-sail trajectory, the following objective function and constraints are defined: Minimize :w1 (af − a|t=tF )2 + w2 (ef − e|t=tF )2 + w3 (if − i|t=tF )2  ≤ τi ≤ 1 (i = 1, · · · , 100)  0 ◦ . Subject to : −90 ≤ α0 ≤ 90◦  0◦ ≤ δ0 ≤ 180◦ 14

(28) (29)

Table 4: Final boundary condition of target orbit.

Venus Mars Itokawa

af [au] ef 0.723332 0.00677323 1.52368 0.0934123 1.32411 0.280168

The boundary condition at the end of the flight time is defined as the function of orbital elements of the target: af , ef , and if . Weight coefficients w1 , w2 , and w3 are also introduced. The spacecraft transverse velocity and sun-spacecraft distance at the terminal t = tF were used as the final boundary condition for circle–to–circle rendezvous missions [12, 19, 20]. When the eccentricity of the target orbit is zero, the boundary condition enables us to leave the spacecraft angular position and velocity free. In this paper, on the other hand, defining the final position and velocity of the spacecraft does not enable us to leave the spacecraft angular position free because the eccentricity of the target orbit is not zero. Hence, in order to investigate the rendezvous capability of the E-sail, the final boundary is given as a set of the spacecraft semimajor axis, eccentricity, and inclination. To generate the initial trajectory for the E-sail, LSLs can be used. The problem is solved by the sequential quadratic programming method [30], and the trajectories are optimized for Eqs. (28) and (29). 4.2. Two-dimensional test cases Here, E-sail trajectories are optimized in a two-dimensional coordinate system. Hence, the orbital inclination is assumed to be zero, and the ascending node is not defined. The trajectory of the E-sail is optimized to inject a spacecraft into the orbits of Venus, Mars, and Itokawa. Table 4 summarizes the orbital elements of the target orbits. Note that the orbital elements of Itokawa were taken from the NASA Jet Propulsion Laboratory (JPL) website4 . Looking at Table 4, one can see that the values of the eccentricity of Mars and Venus are relatively small. On the other hand, Itokawa has a larger value of eccentricity at 0.28, and this can make the problem more difficult. Figure 13 summarizes details of the Earth–Venus mission. Figure 13a is the optimal trajectory for Earth–Venus transfer, b is the time history 4

https://cneos.jpl.nasa.gov/ (Accessed on Aug. 8, 2017)

15

of the semimajor axis, and c is the time history of the eccentricity. After 1464 days of flight time, the E-sail is injected into the target orbit. Since the semimajor axis of Venus is smaller than 1 au, the LSL for decreasing the semimajor axis is used to generate the initial trajectories of the E-sail. As Fig. 13b shows, the value of the semimajor axis gradually decreases during the flight phase. The value of the eccentricity, on the other hand, fluctuates between 0 and 0.15, as shown in Fig. 13c. The error in the semimajor axis is about 1200 km. This is sufficiently small compared with one astronomical unit (∼ 1.5 × 108 km). The error in the eccentricity is also sufficiently small, and the order is smaller than 10−7 . Two orbital elements successfully converged to the target values. Figure 14 shows the time history of the switching parameter τ . Ten areas filled with gray occur if the semimajor axis of the spacecraft is decreased when the value of switching parameter τ is higher than zero. These areas are defined by the LSL for decreasing the semimajor axis. As can be seen in Fig. 14, the E-sail thrust is turned on mainly in these gray areas in order to decrease the semimajor axis until it converges to the value of Venus, 0.7233 au. The time history of switching parameter τ is well accorded with the tendency of the LSL for decreasing the semimajor axis. Details of the Earth–Mars mission are shown in Fig. 15. Unlike the Earth– Venus mission, the semimajor axis of Mars is greater than 1 au, and the LSL for increasing the semimajor axis is used to generate the initial trajectories of the E-sail. The E-sail increases the value of the semimajor axis through the flight phase, as shown in Figs. 15a and b. After 2558 days of flight time, the orbital elements converge to the target values, as shown in Figs. 15b and c. The error in the semimajor axis is on the order of 100 km, and the eccentricity is on the order of 10−7 . Figure 16 shows the time history of switching parameter τ for the Earth–Mars transfer mission. The time history of the switching parameter τ is well accorded with the tendency of the LSL to increase the semimajor axis. Since the target value of the eccentricity is small, the LSL for the semimajor axis well represents the optimized trajectory. The optimization results of the Earth-Venus and Earth-Mars missions are compared with missions optimized in the previous work [19]. In that work, trajectories were optimized assuming an old solar wind force model. In addition, we approximated the orbits of Venus and Mars as circular coplanar orbits. It is useful to estimate how much the E-sail capability for these missions decreases when we consider the actual force model and target orbits. In the previous work, the Earth-Venus mission was completed in 1425 days, 16

Table 5: Final boundary conditions of Ryugu orbit.

Ryugu

af [au] ef if [deg] 1.18956 0.190284 5.88390

and the Earth-Mars mission was completed in 2162 days. The extensions are still acceptable values. The result of the Earth–Itokawa transfer is shown in Fig. 17. Since the semimajor axis of the Itokawa orbit is greater than 1 au, the LSL for increasing the semimajor axis is used to generate the initial trajectories of the E-sail as in the case of the Earth–Mars mission. After 2054 days of flight time, the E-sail is injected into the target orbit, as shown in Fig. 17a. The semimajor axis converges to the target value at the end of the flight time, as shown in Fig. 17b. To put the spacecraft into Itokawa orbit, unlike the preceding two missions, the eccentricity of the orbit should be increased from zero to the target value ef = 0.2802. As Fig. 17c shows, the value of the eccentricity increases and converges to the target value at the end of the flight phase. Figure 18 shows the time history of switching parameter τ . Figure 18a shows the time history with gray areas that occur if the semimajor axis of the E-sail decreases when τ > 0 is satisfied. Figure 18b shows the time history with gray areas that occur if the eccentricity of the E-sail decreases when τ > 0 is satisfied. Since both the semimajor axis and eccentricity should converge to the target values, the time history of switching parameter τ in Fig. 18 can no longer be explained simply with LSLs like the preceding two missions. The E-sail with a maneuvering technique can achieve the orbit transfer targeting even in a more complicated orbit. 4.3. Three-dimensional test case We perform trajectory optimizations for the E-sail in a three-dimensional coordinate system. To this end, sail clock angle δ is considered here. Since the attitude of the E-sail is fixed in the inertial reference frame, the thrust vectoring capability is strongly limited. Thus, the orbital transfer problem in three-dimensional space is clearly a more challenging case. To investigate the capability of the E-sail to inject a spacecraft into inclined orbits, the orbit of the asteroid Ryugu is chosen as a target. Table 5 lists the orbital elements of Ryugu. Note that these orbital elements were also taken from JPL’s website. The final boundary of the ascending node is left free for this problem. We neglect the ascending node because the misalignment in the 17

ascending node can be removed by optimizing the launch date. To complete this mission, appropriate changes in the semimajor axis, eccentricity, and inclination should be achieved simultaneously. The optimal trajectory is shown in Fig. 19. Circles with dotted lines in Fig. 19 are every 0.25 au. Figure 20 provides time histories of the orbital elements. Figure 19 shows that the orbit of the E-sail is gradually inclined and arrived at the target orbit after 1833 days of flight time. Three orbital elements successfully converged to the target values at the end of the flight phase, as shown in Fig. 20. Errors in orbital elements at the end of the flight phase are sufficiently small as those in preceding missions. Figure 21 shows the time history of switching parameter τ . Note that in Fig. 21, areas filled with gray occur if a given orbital element of the spacecraft orbit decreases when the switching parameter τ is higher than zero. Figure 21a is the semimajor axis, b is the eccentricity, and c is the inclination. Since the three orbital elements should converge to the target values at the end of the flight time, changes in the switching parameter are more complicated than in the Earth–Venus, Earth–Mars, and even Earth–Itokawa missions solved in two–dimensional space. The LSLs for increasing the semimajor axis and the eccentricity are unable to represent the changes in τ . The change rate of the inclination by the LSL depends on the initial attitude of the E-sail defined by the sail plane’s inclination angles α and δ, and thus the areas in which the inclination decreases at τ > 0 are smaller, as shown in Fig. 21 c). The optimized initial sail attitude angles are α = −42.21◦ and δ = 123.8◦ , and the thrust force along the out-of-plane direction is generated. The E-sail can inject a spacecraft into such a realistic orbit with a reasonable flight time. 5. Conclusion To investigate the capability of the electric solar wind sail as a propulsion system for deep space exploration missions, we formulated an advanced solar wind force model. A conducting tether was separated into uniform segments, and we numerically calculated the forces exerted on each segment. The results showed that the thrust force vector of the electric solar wind sail is represented by a complex function of the sail plane’s inclination angle and distance from the Sun because the solar wind and centrifugal forces warp the tether. We also presented trajectory optimization results. The twodimensional trajectories calculated within this paper demonstrated that the electric solar wind sail can inject a spacecraft into realistic orbits such as 18

Venus, Mars, and Itokawa in acceptable required flight times. Moreover, we showed the three-dimensional trajectory to the Ryugu orbit. Although details of how to deal with a mechanical instability in the spacecraft at a high sail plane inclination angle and misalignment in the spacecraft ascending node are out of the scope of this paper, such issues have to be considered. The real degradation of the solar wind force and changes in the shape of the tethers are investigated prior to missions. Acknowledgement This work was supported by JSPS KAKENHI Grant Number JP19K15210. Appendix A. Force model for different cases of E-sail rotation period In this paper, the E-sail force model is developed based on the tether rotation period of ∆t = 125 min [23]. Based on the model, numerically optimized trajectories for several missions were shown. However, the E-sail tether shape changes depending on the rotation period, because the centrifugal force changes. The change in the tether shape can affect the thrust direction and magnitude. To provide simulation models for the E-sail with different rotation periods, thrust direction and acceleration magnitude of the E-sail are calculated for some cases of the rotation period. Calculation results are shown in Figs. A.22-A.24. Figure A.22 is for the ∆t = 75 min case, Fig. A.23 is for the ∆t = 100 min case, and Fig. A.24 is the ∆t = 150 min case, respectively. Since warped tethers deviate from the ideal flat sail plane, the thrust capability is degraded as ∆t becomes large. When the rotation period ∆t is large, the peak of the thrust direction αT becomes small. The relations between the sail plane’s inclination angle, force coefficient, thrust direction, and thrust magnitude show same kinds of tendency with the ∆t = 125 min case. The best-fit interpolation of the numerical data is used to find the coefficients Pi,j and Qi,j . They are listed in tables A.6–A.11. Substituting coefficients corresponding to a given rotation period ∆t into Eqs. (16) and (17), we can obtain the E-sail force model.

19

Table A.6: Best-fit interpolation coefficients Pi,j for thrust direction αT (∆t = 75 min case). PP j 0 1 2 3 4 5 PP i P 0 1 2 3 4 5

−7.066 × 10−2 0.5390 −4.648 × 10−3 1.547 × 10−4 −3.154 × 10−6 1.409 × 10−8

1.993 × 10−2 1.372 × 10−2 −1.169 × 10−3 2.602 × 10−5 −1.689 × 10−7 0

−6.817 × 10−2 −1.194 × 10−2 1.253 × 10−4 4.144 × 10−7 0 0

1.774 × 10−2 1.083 × 10−3 −1.894 × 10−5 0 0 0

6.666 × 10−4 4.896 × 10−5 0 0 0 0

−4.161 × 10−4 0 0 0 0 0

Table A.7: Best-fit interpolation coefficients Qi,j for normalized thrust acceleration magnitude γ (∆t = 75 min case). PP j 0 1 2 3 4 5 PP i P 0 1 2 3 4 5

0.9984 5.140 × 10−4 −1.517 × 10−4 1.110 × 10−6 −1.044 × 10−8 9.525 × 10−11

1.734 × 10−3 −4.258 × 10−4 1.974 × 10−5 −3.101 × 10−7 1.541 × 10−9 0

−7.910 × 10−4 1.928 × 10−4 8.433 × 10−7 −2.067 × 10−8 0 0

1.946 × 10−4 −4.309 × 10−5 3.031 × 10−7 0 0 0

−4.796 × 10−5 6.323 × 10−7 0 0 0 0

9.999 × 10−6 0 0 0 0 0

Appendix B. Nomenclature a a0 (ar , aθ , aϕ ) asw asw cα e F F f G G γ i L ∆L l mp mru N np

: : : : : : : : : : : : : : : : : : : : :

semimajor axis characteristic acceleration acceleration in inertial spherical polar coordinate system solar wind acceleration vector solar wind acceleration cosine factor eccentricity solar wind force vector solar wind force true anomaly centrifugal force vector centrifugal force normalized thrust acceleration magnitude inclination length of main tether unit length of main tether semilatus rectum mass of proton mass of remote unit number of main tethers number density of proton 20

Table A.8: Best-fit interpolation coefficients Pi,j for thrust direction αT (∆t = 100 min case). PP j 0 1 2 3 4 5 PP i P 0 1 2 3 4 5

−5.029 × 10−2 0.5076 −1.622 × 10−3 5.909 × 10−5 −1.929 × 10−6 8.571 × 10−9

0.3788 −4.661 × 10−3 −1.517 × 10−3 4.280 × 10−5 −2.903 × 10−7 0

−0.5303 −1.749 × 10−2 6.970 × 10−5 2.147 × 10−6 0 0

0.2313 2.952 × 10−3 −5.394 × 10−5 0 0 0

−3.768 × 10−2 2.151 × 10−4 0 0 0 0

1.346 × 10−3 0 0 0 0 0

Table A.9: Best-fit interpolation coefficients Qi,j for normalized thrust acceleration magnitude γ (∆t = 100 min case). PP j 0 1 2 3 4 5 PP i P 0 1 2 3 4 5

0.9988 6.948 × 10−4 −1.721 × 10−4 1.862 × 10−6 −2.180 × 10−8 1.551 × 10−10

n ˆ p Pi,j Qi,j rˆ (r, θ, ϕ) sα (S, T, W ) t tF ∆t T u V V0 Vmax Vsw v v v⊥ vsw

: : : : : : : : : : : : : : : : : : : : :

−2.587 × 10−3 −5.158 × 10−4 2.748 × 10−5 −3.297 × 10−7 7.442 × 10−10 0

6.051 × 10−3 4.064 × 10−4 3.724 × 10−7 −9.310 × 10−9 0 0

−4.277 × 10−3 −1.237 × 10−4 3.442 × 10−7 0 0 0

1.152 × 10−3 7.751 × 10−6 0 0 0 0

−9.924 × 10−5 0 0 0 0 0

ideal sail plane normal vector unit vector along out-of-orbit plane best-fit coefficients for thrust direction best-fit coefficients for thrust acceleration magnitude unit vector along Sun-spacecraft direction heliocentric inertial spherical polar coordinate sine factor radial, transverse, and out-of-plane thrust components time flight time rotation period of the spacecraft main body total force exerted on tether local tether tangent tan λ net tether voltage voltage of main tethers maximum tether voltage bulk kinetic energy of protons spacecraft velocity vector spacecraft velocity solar wind velocity vector perpendicular to the tether direction solar wind velocity

21

Table A.10: Best-fit interpolation coefficients Pi,j for thrust direction αT (∆t = 150 min case). PP j 0 1 2 3 4 5 PP i P 0 1 2 3 4 5

−0.1168 0.4751 2.876 × 10−3 −1.039 × 10−4 4.177 × 10−7 −3.253 × 10−9

1.454 −0.1200 −8.960 × 10−4 5.553 × 10−5 −3.170 × 10−7 0

−2.220 6.487 × 10−3 −3.364 × 10−4 1.903 × 10−6

1.279 4.307 × 10−3 −2.712 × 10−5 0

−0.3114 −1.177 × 10−4 0 0

2.640 × 10−2 0 0 0

0

0

0

0

Table A.11: Best-fit interpolation coefficients Qi,j for normalized thrust acceleration magnitude γ (∆t = 150 min case). PP j 0 1 2 3 4 5 PP i P 0 1 2 3 4 5

1.002 3.905 × 10−4 −1.661 × 10−4 1.903 × 10−6 −2.514 × 10−8 1.830 × 10−10

(vr , vθ , vϕ ) w1 , w2 , w3 (x, y, z) (x∗ , y ∗, z ∗ ) α α0 αρz ′ αT β δ δ0 ǫ ǫ0 λ µ µ⊙ φroll (ρ, ψ, z) (ρ, ψ ′ , z ′ ) ρ0 ρL τ ω ξ ξ0

: : : : : : : : : : : : : : : : : : : : : : : : :

−2.039 × 10−2 9.624 × 10−5 4.559 × 10−5 −3.268 × 10−7 −1.962 × 10−9 0

3.397 × 10−2 1.846 × 10−4 −1.830 × 10−6 8.583 × 10−8 0 0

−2.297 × 10−2 −1.677 × 10−4 −1.098 × 10−6 0 0 0

6.698 × 10−3 2.986 × 10−5 0 0 0 0

−6.965 × 10−4 0 0 0 0 0

velocity components in the spherical polar coordinate system weight coefficients rotating reference frame body fixed reference frame sail plane’s inclination angle initial sail plane’s inclination angle solar wind coming angle in the ρz ′ -plane thrust direction angle between x- and ρ-axes sail plane’s clock angle initial sail plane’s clock angle sufficient small value conductivity of vacuum local tether coning angle linear mass density of tether heliocentric gravitational constant rotation angle around ρ-axis rotating reference frame rotating reference frame root of a tether position in ρ-axis position of remote unit in ρ-axis switching parameter rotation period of spacecraft main body E-sail force factor E-sail reference force22factor

Subscripts f ρ z

: target value : component along ρ-axis : component along z-axis

Superscripts (i) : value after i-th iteration i : value of i-th tether segment References [1] P. Janhunen, Electric Sail for Spacecraft Propulsion, journal of proulsion and power 20 (4) (2004) 763–764, doi:10.2514/1.8580. [2] P. Janhunen, P. K. Toivanen, J. Polkko, S. Merikallio, P. Salminen, E. Heggstorm, H. Sepp¨anen, R. Kurppa, J. Ukkonen, S. Kiprich, G. Thornell, H. Kratz, L. Tichter, O. Kr¨omer, R. Rosta, M. Noorma, J. Envall, S. L¨att, G. Mengali, A. A. Quarlta, H. Koivisto, O. Travainen, T. Kalvas, J. Kauppinen, A. Nuottaj¨arvi, A. Obraztsov, Electric Solar Wind Sail: Towards Test Missions, Review of Scientific Instruments 81 (11) (2010) 1–11, doi:10.1063/1.3514548, paper111301. [3] P. Janhunen, A. Sandroos, Simulation Study of Solar Wind Push on a Charged Wire: Basis of Solar Wind Electric Sail Propulsion, Annales Geophysicae 25 (3) (2007) 755–768, doi:10.5194/angeo-25-755-2007. [4] P. Janhunen, The Electric Sail - a New Propulsion Method Which May Enable Fast Missions to the Outer Solar System, Journal of the British Interplanetary Society 61 (2008) 322–325. [5] K. Hoshi, H. Kojima, H. Yamakawa, Numerical Analysis of Potential Structure around Electric Solar Wind Sail, The 30th ISTS Special Issue of Transactions of JSASS, Aerospace Technology Japan 14 (ists30) (2016) 83–89, doi:10.2322/tastj.14.Pb 83, URL http://doi.org/10.2322/tastj.14.Pb 83. [6] K. Hoshi, H. Kojima, T. Muranaka, H. Yamakawa, Thrust Calculation of Electric Solar Wind Sail by Particle-In-Cell simulation, Annales Geophysicae 34 (9) (2016) 845–855, doi:10.5194/angeo-34-845-2016, URL http://doi.org/10.5194/angeo-34-845-2016. 23

[7] P. Janhunen, A. A. Quarta, G. Mengali, Electric Solar Wind Sail Mass Budget Model, Geocentric Instrumentation Methods and Data Systems 2 (2013) 85–95, doi:10.5194/gi-2-85-2013. [8] H. Sepp¨anen, T. Rauhala, S. Kiprich, J. Ukkonen, M. Simonsson, R. Kurppa, P. Janhunen, E. Hæggstr¨om, One Kilometer (1 km) Electric Solar Wind Sail Tether Produced Automatically, Review of Scientific Instruments 84, doi:10.1063/1.4819795, URL http://aip.scitation.org/toc/rsi/84/9. [9] A. Slavinskis, H. Ehrpais, K. Zalite, R. Rosta, M. Pajusalu, P. Liias, K. Kahn, T. Kalvas, H. Kuuste, E. Kulu, S. L¨att, R. Vendt, E. Ilbis, J. Viru, J. Envall, V. Allik, T. onis Eenm¨ae, J. Kalde, P. Toivanen, M. Noorma, I. S¨ unter, U. Kvell, J. Polkko, K. Laizans, J. K¨ utt, P. Janhunen, ESTCube-1 in-orbit experience and lessons learned, IEEE Aerospace and Electronics Systems Magazine 30 (8) (2015) 13–22, doi: 10.1109/MAES.2015.150034. [10] O. Khurshid, T. Tikka, J. Praks, M. Hallikainen, Accomodating the Plasma Brake Experimant On-board the Aalto-1 Satellite, in: Estonian Academy of Sciences, vol. 63, Aalto University, 02150 Espoo, Finland, 258–266, doi:10.3176/proc.2014.2S.07, 2014. [11] A. Kestil¨a, T. Tikka, P. Peitso, J. Rantanen, A. N¨asil¨a, K. Nordling, H. Saari, R. Vainio, P. Janhunen, J. Praks, M. Hallikainen, Aalto-1 nanosatellite–technical description and mission objectives, Geoscientific Instrumentation, Methods and Data Systems 2 (1) (2013) 121–130, doi:10.5194/gi-2-121-2013, URL https://www.geosci-instrum-method-data-syst.net/2/121/2013/. [12] G. Mengali, A. A. Quarlta, P. Janhunen, Electric Sail Performance Analysis, Journal of Spacecraft and Rockets 45 (1) (2007) 122–129, doi: 10.2514/1.31769. [13] A. A. Quarta, G. Mengali, Electric Sail Mission Analysis for Outer Solar Exploration, Journal of Spacecraft and Rockets 33 (3) (2008) 740–755, doi:10.2514/1.47006. [14] A. A. Quarta, G. Mengali, Electric Sail Missions to Potentially Hazardous Asteroids, Acta Astronautica 66 (9–10) (2013) 1506–1519, doi: 10.1016/j.actaastro.2009.11.021. 24

[15] A. A. Quarta, G. Mengali, P. Janhunen, Electric Sail for a Near-Earth Asteroid Sample Return Mission: Case 1998 KY26, Journal of Aerospace Engineering 27 (6), doi:10.1061/(ASCE)AS.1943-5525.0000285. [16] G. Mengali, A. A. Quarta, Optimal Nodal Flyby with Near-Earth Asteroids Using Electric Sail, Acta Astronautica 104 (2014) 450–457, doi: 10.1016/j.actaastro.2014.02.012. [17] K. Yamaguchi, H. Yamakawa, Electric Solar Wind Sail Kinetic Energy Impactor for Asteroid Deflection Missions, The Journal of the Astronautical Sciences 63 (1) (2016) 1–22, doi:10.1007/s40295-015-0081-x. [18] S. Merikallio, P. Janhunen, Moving an Asteroid with Electric Solar Wind Sail, Astrophysics and Space Sciences Transactions 6 (2010) 41–48, doi: 10.5194/astra-6-41-2010. [19] K. Yamaguchi, H. Yamakawa, Orbital Maneuver of Electric Sail Using ON/OFF Thrust Control, Koku-Uchu-Gijyutsu(Aerospace Technology Japan, the Japan Society for Aeronautical and Space Sciences) 12 (2013) 79–88, doi:org/10.2322/astj.12.79. [20] A. A. Quarta, G. Mengali, Minimum-time trajectories of electric sail with advanced thrust model, Aerospace Science and Technology 55 (2016) 419–430, doi:10.1016/j.ast.2016.06.020. [21] M. Huo, G. Mengali, A. A. Quarta, Electric Sail Thrust Model from a Geometrical Perspective, Journal of Guidance, Control and Dynamics 41 (3) (2018) 734–740, doi:10.2514/1.G003169, URL https://arc.aiaa.org/doi/pdf/10.2514/1.G003169. [22] P. K. Toivanen, P. Janhunen, Spin Plane Control and Thrust Vectoring of Electric Solar Wind Sail, Journal of Propulsion and Power 29 (1) (2013) 178–185, doi:10.2514/1.B34330. [23] P. K. Toivanen, P. Janhunen, Thrust vectoring of an electric solar wind sail with a realistic sail shape, Acta Astronautica 131 (2017) 145–151, doi:10.1016/j.actaastro.2016.11.027. [24] M. Bassetto, G. Mengali, A. A. Quarta, Thrust and torque vector characteristics of axially-symmetric E-sail, Acta Astronautica 146 (2018) 134–143, doi:10.1016/j.actaastro.2018.02.035, URL https://doi.org/10.1016/j.actaastro.2018.02.035. 25

[25] M. Bassetto, G. Mengali, A. A. Quarta, Stability and Control of Spinning Electric Solar Wind Sail in Heliostationary Orbit, Journal of Guidance, Control, and Dynamics 42 (2) (2019) 425–431, doi: 10.2514/1.G003788. [26] P. Janhunen, Increased electric sail thrust through removal of trapped shielding electrons by orbit chaotisation due to spacecraft body, Annales Geophysicae 27 (8) (2009) 3089–3100, doi:10.5194/angeo-27-3089-2009. [27] E. C. Sittler, J. D. Scudder, An Empirical Polytrope Law for Solar Wind Thermal Electrons Between 0.45 and 4.76 AU: Voyager and Mariner 10, Journal of Geophysical Research 85 (A10) (1980) 5131–5137, doi: 10.1029/JA085iA10p05131. [28] J. A. Slavin, R. E. Holzer, Solar Wind Flow about the Terrestrial Planets, 1. Modeling Bow Shock Position and Shape, Journal of Geophysical Research 86 (11) (1981) 11401–11418, doi:10.1029/JA086iA13p11401. [29] M. Bassetto, G. Mengali, A. Quarta, Attitude Dynamicsof an Electric Sail Model with a Realistic Shape, Acta Astronautica 159 (2019) 250– 257, doi:10.1016/j.actaastro.2019.03.064. [30] T. Ibaraki, M. Fukushima, Optimization Programming, Iwanami Computer Science, Iwanami, Tokyo, Japan, 1991.

26

E-sail

Figure 1: Two coordinate systems to determine E-sail attitude.

27

‫ݖ‬௅

‫ܨ‬ఘ௜

ߣ

ܶఘ Tether

‫ܩ‬ Remote unit

‫ܨ‬௜

‫ܨ‬௭௜

Figure 2: E-sail tether in 2D coordinate system.

28

ܶ௜

‫ݖ‬

Solar wind

ߙ

ߚ ‫ݔ‬ ߩ

Figure 3: Relationship between body-fixed frame (x, y, z) and tether-fixed frame (ρ, ψ, z).

29

-plane

Tether Remote unit

Figure 4: Attitude of coordinate system (ρ, ψ ′ , z ′ ) relative to (ρ, ψ, z) and conducting tether.

30

z* [km]

z [km]

0 -10 -20-20

asw/a0 -10

-10

0 y* [km]

10

2020

10

-10 0 x* [km]

-20 -20 20

10 x* [km]

a) Tether shape when α = 0o.

20 20 0.5 |

10 z* [km]

z* [km]

10 0 -10 -20-20

asw/a0 n

0 -10

-10

0 y* [km]

10

2020

10

-10 0 x* [km]

-20 -20 20

10 x* [km]

b) Tether shape when α = 45o.

20 20 n

10 z* [km]

z* [km]

10 0 -10 -20-20

asw/a0

0 -10

-10

0 y [km] *

10

2020

10

-10 0 * x [km]

-20 -20 20

10 x* [km]

c) Tether shape when α = 90o.

31 Figure 5: Shape of E-sail tethers for three patterns of sail plane’s inclination angle α.

Proposed tether shape analysis  = 125 min  ∆t = 100 min Old polynomial model Normalized thrust acceleration magnitude γ

Thrust direction αT [deg]

20

15

10

5

1 0.9 0.8 0.7 0.6 0.5

0 0

10

20 30 40 50 60 70 80 Sail plane’s inclination α [deg]

90

0

a) Thrust direction over sail plane’s inclination angle.

10

20 30 40 50 60 70 80 Sail plane’s inclination α [deg]

90

b) Thrust acceleration magnitude over sail plane’s inclination angle.

Figure 6: Thrust direction and acceleration magnitude as a function of sail plane’s inclination angle α.

32

Normalized force factor ξ / ξ

3.5

5

3

10°

2.5

10° 15°

2 1.5 1 0.5 0 0 10 20 30 40 50 60 70 80 Sail plane’s inclination angle α [deg]

Figure 7: Numerical result of thrust direction as a function of the sail plane’s inclination angle and E-sail force factor.

33

Thrust magnitude

0.9

0.8

0.8 0.7

0.7 0.6

0.6

0.5

4 3

r 2 cto 0 10 a f 20 30 e 1 40 50 orc Sail plan f 60 d e’s inclin 0 ation ang 70 80 90 lize a le α [deg rm ] No

Figure 8: Numerical result of thrust acceleration magnitude as a function of sail plane’s inclination angle and E-sail force factor.

࢔ : Ideal spin plane normal E-sail center of mass

࢖ ෝ ൈ ࢘ො ࢇ ୱ୵

ߜ

ߙ ߙ் Sun Orbit plane

Figure 9: Sail plane’s inclination angle α, thrust direction αT , and clock angle δ.

34

࢘ො





Sun
 0∘






Figure 10: Heliocentric inertial spherical polar coordinate system.

35

150

o

Sun

180o

0.

5

1 1. 5

210

2

o

2.

5

au

au

au

au

au

240o 270

Figure 11: Example of locally optimal trajectory of E-sail using LSLs to increase semimajor axis.

36

Semimajor axis a [au]

Proposed forcemodel Old force model

2 1.75 1.5 1.25 1 0

500

1000 1500 2000 2500 3000

Semimajor axis a [au]

Time t [d] a) LSL for increasing semimajor axis. 1 0.75 0.5 0.25 0 0

500

1000 1500 2000 2500 3000 Time t [d] b) LSL for decreasing semimajor axis.

0

500

Eccentricity e

1 0.75 0.5 0.25 0 1000 1500 2000 2500 3000 Time t [d] c) LSL for increasing eccentricity.

Figure 12: Time histories of orbital elements using LSLs.

37

Semimajor axis a

Venus 120

150

90

o

o

60

o

o

o

30

0.9 0.8 a = 0.723332 au 0.7 0

200

400

600

800

1000 1200 1400

Time t [d]

0.5

0. 7

5

1

au o 210 1 .2 5 a

b) Time history of semimajor axis.

Departure o 0

Sun

o

0.15 au

au

Arrival 330

Eccentricity e

180

o

u

240

o

270

o

300

o

0.1 0.05 e = 0.00677323

0 0

200

400

600

800

1000 1200 1400

Time t [d] c) Time history of eccentricity.

a) Optimal Earth - Venus trajectory for E-sail.

Switching parameter

Figure 13: Optimal Earth–Venus transfer for a0 = 0.5 mm/s2 E-sail.

0 0

200

400

600

800

1000

Time t [d]

Figure 14: Time history of change in switching parameter.

38

1200

1400

Semimajor axis a

Mars 90o 120o

60o 30o

150o

1.5 1.25 1 0

500

1000

1500

2000

2500

2000

2500

Time t [d] Sun

180o 0. 5

1

75

au 25 210 1 5 au au .7 5 au o

1.

1.

au

au

Eccentricity e

0.

b) Time history of semi-major axis.

0o

Departure

330o

Arrival

240o

270o

300o

0.1

e = 0.0934123

0.05

0 0

a) Optimal Earth - Mars trajectory for E-sail.

500

1000 1500 Time t [d]

c) Time history of eccentricity.

Switching parameter

Figure 15: Optimal Earth–Mars transfer for a0 = 0.5 mm/s2 E-sail.

0 0

500

1000

1500 Time t [d]

Figure 16: Time history of change in switching parameter.

39

2000

Semimajor axis a [au]

Itokawa 90o

120o

60o

o

o

30

150

1.3 1.2 1.1 1 0.9 0

500

Arrival

210o

5

75

au

1500

2000

0.4

au au

1. 25 5 a au 75 u 2 au au

330o

1.

240o

2000

0o

Eccentricity e

1 1.

Departure

0.

0.

1500

b) Time history of semi-major axis.

Sun

180o

1000 Time t [d]

270o

300o

e = 0.280168

0.3 0.2 0.1 0 0

a) Optimal Earth - Itokawa trajectory for E-sail.

500

1000 Time t [d]

c) Time history of eccentricity.

Figure 17: Optimal Earth–Itokawa transfer for a0 = 0.5 mm/s2 E-sail.

40

Switching parameter

1 0.75 0.5 0.25 0 0

500

1000

1500

Time t [d]

Switching parameter τ

a) Areas in which the semimajor axis is decreased when τ > 0 (filled with gray).

1 0.75 0.5 0.25 0 0

500

1000 Time t [d]

1500

b) Areas in which the eccentricity is decreased when τ > 0 (filled with gray). Figure 18: Time history of change in switching parameter.

41

Figure 19: Optimal Earth–Ryugu transfer for a0 = 0.5 mm/s2 E-sail.

42

Semimajor axis a

a = 1.18956 au

1.2 1.1 1 0

200

400

600

800 1000 Time t [d]

1200

1400

1200

1400

1200

1400

Eccentricity e

a) Time history of semimajor axis.

e = 0.190284

0.2 0.1 0 0

200

400

600

800 1000 Time t [d]

Inclination i [deg]

b) Time history of eccentricity.

i = 5.88390 deg

6 4 2 0 0

200

400

600

800 1000 Time t [d]

c) Time history of inclination.

Figure 20: Time histories of changes in orbital elements of Earth–Ryugu transfer.

43

Switching parameter

1 0.75 0.5 0.25 0 0

200

400

600

800 1000 1200 Time t [d]

1400

1600

1800

Switching parameter τ

a) Areas in which the semi-major axis is decreased when τ > 0 (filled with gray). 1 0.75 0.5 0.25 0 0

200

400

600

800

1000

1200

1400

1600

1800

Time t [d]

Switching parameter τ

b) Areas in which the eccentricity is decreased when τ > 0 (filled with gray). 1 0.75 0.5 0.25 0 0

200

400

600

800 1000 1200 Time t [d]

1400

1600

c) Areas in which the inclination is decreased when τ > 0 (filled with gray). Figure 21: Time history of change in switching parameter.

44

1800

a) Thrust direction over sail plane’s inclination and normalized force factor

b) Normalized thrust magnitude over sail plane’s inclination and normalized force factor

Figure A.22: Thrust direction and acceleration magnitude as a function of sail plane’s inclination angle and normalized force factor (∆t = 75 min).

a) Thrust direction over sail plane’s inclination and normalized force factor

b) Normalized thrust magnitude over sail plane’s inclination and normalized force factor

Figure A.23: Thrust direction and acceleration magnitude as a function of sail plane’s inclination angle and normalized force factor (∆t = 100 min).

45

a) Thrust direction over sail plane’s inclination and normalized force factor

b) Normalized thrust magnitude over sail plane’s inclination and normalized force factor

Figure A.24: Thrust direction and acceleration magnitude as a function of sail plane’s inclination angle and normalized force factor (∆t = 150 min).

46

Electric solar wind sail generates thrust by solar wind and tether interaction By investigating the change in the tether shape, advanced force model is developed Orbital maneuvering technique that doesn’t change the spacecraft attitude is proposed Earth-Mars, Venus, and asteroids missions are numerically demonstrated Results show that the system can complete the missions in reasonable flight times