Contents

1 Introduction1

1.1

Deep space exploration . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

Space Transportation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.3

Laser Fusion Rocket (LFR) . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.4

Previous researches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.4.1

Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.4.2

Replication of experimental result . . . . . . . . . . . . . . . . . .

6

1.4.3

Thrust deflection . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

Aim of this study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

1.5

2 Methodology
2.1

2.2

12

Numerical simulation codes . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.1.1

Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.1.2

Radiation hydrodynamics code . . . . . . . . . . . . . . . . . . . .

13

2.1.3

Hybrid code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

2.1.4

Coupling code . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

Simulation models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

2.2.1

Replication of experimental result . . . . . . . . . . . . . . . . . .

20

2.2.2

Thrust deflection . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

3 Results

23

3.1

23

Replication of experimental result . . . . . . . . . . . . . . . . . . . . . .

1

3.2

Thrust deflection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4 Discussion

32
36

4.1

Replication of experimental result . . . . . . . . . . . . . . . . . . . . . .

36

4.2

Thrust deflection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

5 Summary43

List of figures

44

List of tables

46

Bibliography

46

Acknowledgment

502Chapter 1
Introduction
1.1 Deep space exploration
Mankind has planned and done many space exploration projects up to this day. The
farthest spacecraft, Voyager 1, escaped from the heliosphere [1]. Twelve astronauts walked
on the moon [2]. These spacecraft or spaceships were launched on sixties or seventies. Afterthese projects, there is no realized project that exceeds them for the distance from the Earth.
A game-changing ticket to cruise the deep space is a cornerstone for the continuing success of human beings. The resources on the Earth will be sure to dry up. A severe meteorite
impact might be occurred in the future. To brace for a planetary scale crisis, we should get a
way to move freely in space. Konstantin Tsiolkovsky, who was a Russian and Soviet rocket
scientist and pioneer of the astronautic theory, said that ”Earth is the cradle of humanity, but
one cannot live in a cradle forever”. Someday, we should go out from the comfortable home
planet.

1.2 Space Transportation
Current main space propulsion systems [3] are categorized as two types, chemical andelectric propulsion as shown in Fig. 1.1 and 1.2, respectively. Chemical propulsion can
produce large thrust (over MN), thus, all of launch vehicles use chemical propulsion sys-1tem. However, the specific impulse of chemical propulsion is not so high (several hundred
seconds). It means that chemical propulsion requires large amount of propellant to get large
acceleration, ∆V. On the other hand, electric propulsion can produce high specific impulse(several thousands seconds) but the thrust is quite low (up to N class) due to the restrictionof power generation in space. In general, large thrust and high specific impulse are incompatible.

Figure 1.1: Launch of the first Epsilon

Figure 1.2: Deep Space 1 with the NSTAR Ion

Launch Vehicle (Epsilon-1)[4]Engine[5]Manned mars exploration project is ongoing but it takes over three years to do a roundtrip with current technologies [6]. NASA said that ”Radiation exposure for human explorers
could exceed NASA’s career limit for astronauts if current propulsion systems are used [7]”during a mars exploration. The data was measured by radiation assessment detector thatwas mounted on Curiosity, NASA’s mars rover. Thick wall is required as radiation shield
while spacecraft should be light. Most suitable solution is shortening a time of exposure
to radiation. Thus, a new propulsion system that can reduce a mission period is strongly
required.

2

1.3 Laser Fusion Rocket (LFR)
LFR is one of the candidates for interplanetary manned flight in the future because it cangenerate both large thrust and high specific impulse simultaneously. The concept design of
the LFR was proposed by Lawrence Livermore National Laboratory[[8]].The image of theflight model of VISTA is shown in figure 1.3. It is estimated that we can do a round trip tomars in half year. In LFR system, fusion plasma is controlled by an external magnetic field.
The magnetic field controls the direction of plasma flow and converts its thermal energy tokinetic energy. This magnetic nozzle system is called

magnetic thrust chamber . The

mechanism of magnetic thrust chamber is shown in figures 1.4(a)-1.4(c). At first, plasma
produced by inertial confinement fusion expands in an external magnetic field produced byan electromagnetic coil. Then, the plasma compresses the magnetic field until the magnetic
pressure balances with the plasma pressure. At the same time, a magnetic cavity appears ina plasma due to a diamagnetic current that is produced by Larmor motion. As a result, theplasma is accelerated by a Lorentz force between the diamagnetic current and the external
magnetic field. Then, the coil is accelerated to opposite direction by a reaction force.

Figure 1.3: Flight model image of VISTA

3

Figure 1.4: The mechanism of a magnetic thrust chamber
Single solenoid coil system is the basic and simplest thrust generation system for a LFR.
This system is shown in figures 1.4.
Multi coil system, it was called ”square coil system” in the previous researches, is theadvanced coil system for a LFR that was proposed by Dr. Kajimura. Figure 1.5 shows its
concept. This system consists of multiple coils that are controlled independently. In thisfigure, when the upper coils and lower coils are driven by different current, the magnetic
field structure deflects from axial symmetry. As a result, plasma ejection vector deflects tothe weaker side of the magnetic field. Thus, the plasma ejection vector (i.e. thrust vector)
will be controllable using this system. When all coils are driven by same current, multi coilsystem works on a parity with a single solenoid coil system.

4

Figure 1.5: Concept image of multi coil system
This system will be used to control the attitude of space craft. On the acceleration phase,
LFR ejects plasma backwords. LFR gets large ∆V in this phase to trip in space with highspeed. On the other hand, LFR must decrease its speed on the rendezvous phase. Reverse
thrust is required on this phase because there are not external braking force in space. Thus,
LFR must ejects plasma to the direction of traveling. Therefore, LFR must make a turnbetween the acceleration or cruising phase and the deceleration phase as Fig. 1.6.

Figure 1.6: The image of the mission phase of LFR

5

1.4 Previous researches
1.4.1 Overview
The flight model of LFR will use a fusion plasma as its energy source. However, it is really difficult to create and use laser fusion plasmas in an experiment on the current situation
because huge laser energy more than a few MJ is required. Thus, we have conducted numerical and experimental researches on a magnetic thrust chamber using laser-produced plasmas
to prove the concept of a magnetic thrust chamber. In previous experimental researches with
Gekko-XII (GXII) laser (maximum output energy of 1 kJ/beam) and Extreme Ultra Violet(EUV) database laser facilities (∼10 J) at the Institute of Laser Engineering, Osaka University, we have examined the behaviors of plasmas and magnetic fields with the observationof self-emission[9, 10, 11], interferometry, and Thomson scattering, and with B-dot probes,
respectively. In addition, thrust generation has been confirmed with magnetic thrust chamber
experimentally using pendulum thrust stand[12, 13]. In previous numerical researches, wehave simulated plasma generation process by a laser irradiation using one-dimensional (1D)
and two-dimensional (2D) radiation hydrodynamics codes[14, 15], and interaction between
plasmas and magnetic fields using three-dimensional (3D) hybrid code[16].

1.4.2 Replication of experimental result
Some differences between the experimental and the numerical results were pointed out
in previous researches. Figure 1.7 shows the experimental setup. A 500 µm-diameter spheretarget made of polystyrene was placed at target chamber center (TCC) and irradiated by alaser beam. An electromagnetic coil was placed 11 mm away from the target to generate themagnetic field of 1.1 T approximately at the target position. A∼F are the charge collectorsthat were placed 150 mm away from TCC to collect ions in a laser produced plasma. The
aperture diameter of each charge collector is 1 mm. The experimental conditions are shownin Table 1.4.2.

6

Table 1.1: Experimental parameters

Beam shape

Gaussian

Number of beams

1

Energy7.5 J

Pulse width
Laser

9.4 ns(FWHM)
Spot diameter
500 µm(FWHM)
Wave length1064 nm

Inner diameter

13 mmOuter diameter

25 mm

Length10 mm

Number of turns

96

Magnetic field strength

1.1 T(at TCC)(approx.)Shape

SphereDiameter

500 µmMaterial

Polystyrene

Coil

Figure 1.7: Setup of EUV experiment

TargetFigure 1.8 and 1.9 show the comparison between experimental and numerical results.
Figure 1.8 shows the comparison in ion currents densities. In this figure, the experimental result was measured by the charge corrector at 15 degrees. The numerical result is theconverted value from the number of ions that was counted at the same position of charge corrector to the ion current density. On the experimental result, the first peak is at 1 µs and theion current density at the first peak is 115 mA/mm2 . On the other hand, the first ion detection on the numerical result is 2.5 µs and the ion current density is only several mA/mm2 .

7

Figure 1.8: Comparison of numerical result with experimental one
Figure 1.9 shows the comparison between self emission observation in the experimentand the ion density distribution in the numerical simulation. Both results are at 0.5 µs after
the laser timing. The intensity and density are comparable because the self emission fromplasma is proportional to square value of its density. On the experimental result, the highintensity area is expend ± 5 mm on the z axis. On the other hand, the high density area on
the numerical result is smaller than the experimental result and it only stands on the negative
area of z.

8Figure 1.9: Comparison with (a)self-emission and (b)Ion density distribution at 0.5 µs

1.4.3 Thrust deflection
In the previous research, multi coil system was proposed [17]. In this research, deflectionof plasma flow was observed numerically as Figs. 1.10 and 1.11. The maximum value ofdeflection angle was 5.2 degrees. This research was conducted based on the flight model
parameters. It assumed the isotropic fusion plasma as the initial condition of plasmas.
A laser ablation plasma is not an isotropic plasma, generally. It’s necessary to simulate alaser ablation process to compare the numerical and experimental result.

9Figure 1.10: The ion flow behaviors with multi coil system [17]Figure 1.11: Previous results of thrust deflection system [17]10

1.5 Aim of this study
The first aim of this thesis is to develop a numerical simulation that can reproduce experimental results. The final goal is quantitative reproduction of plasma parameters and magnetic
field structures in a magnetic thrust chamber. In this thesis, the goal is set as a qualitative
reproduction of plasma behavior comparing with previous experimental results.
The second aim of this thesis is the validation of a thrust deflection system numerically. In
the previous researches, the numerical simulations with flight model parameters have been
conducted. In this thesis, the goal is to validate the thrust deflection system with realistic
experimental parameters.

11

Chapter 2
Methodology
2.1 Numerical simulation codes
2.1.1 Overview
Our numerical simulation consists of multiple codes. Figure 2.1 shows the numericalsimulation procedure.

Figure 2.1: Numerical simulation procedure
At first, one-dimensional (1D) or two-dimensional (2D) radiation hydrodynamics code
calculates the laser ablation process. After that, the result from radiation code is expanded tothree-dimensional (3D) data. 1D and 2D data are expanded as spherical and axial symmetry,12respectively. Thus, 1D radiation hydrodynamics code is suitable to calculate a sphericalsymmetry plasma expansion like a fusion plasma from a spherical symmetric irradiation.
On the other hand, 2D radiation code is suitable to calculate asymmetric plasma expansion
like an expand with a single beam. The 3D hybrid code receives the expanded data and the
3D initial magnetic field data. Using these data, 3D hybrid code calculates the interaction
between a plasma and an external magnetic field.

2.1.2 Radiation hydrodynamics code
The laser irradiation process is calculated by the radiation hydrodynamics code that wasdeveloped by Dr. Sunahara[14, 15]. Here, 2D radiation hydrodynamics code is used to
calculate one side irradiation condition. It is the same condition as the experiments on EUV
database laser facility at Osaka University. In this code, a laser beam is considered as single
and parallel beam with gaussian shape both in space and time. The calculation mesh typeis Eulerian and all of simulation space except the target is filled with low density gas. The
composition of background gas is assumed as same as target material. This code can handle
only one material. We use averaged parameters in case of multiple ion species.

basic equation
The basic equations of this code are shown in equations (2.1)-(2.4): equations of continuity, motion, ion energy, and electron energy, respectively. The nomenclature on this partis shown in table 2.1. As can be seen from these equations, two temperature model is usedin this code. In addition, it is assumed that the radiation heating affects only electrons.

13

Table 2.1: Nomenclature for the equations of radiation hydrodynamics code
t

time

cvi

specific ion heatρmass densityTeelectron temperature⃗v

velocityTi

ion temperaturepeelectron pressureα

Spitzer relaxation time [18]pi

ion pressure

κeelectron conductivity [19]ptotaltotal pressure

κiion conductivity [20]qartificial viscosity

QL

laser heating

cve

specific electron heat

Qr

radiation heating

Dρ
Dt
Dv
ρ
Dt
DTi
ρcvi
Dt
DTe
ρcve
Dt

= −ρ∇ · ⃗v

(2.1)

= −∇(ptotal + q)

(2.2)

= −(pT Hi + q)∇ · ⃗v + ∇ · (κi ∇Ti ) + α(Te − Ti )(2.3)

= −pT He ∇ · ⃗v + ∇ · (κe ∇Te ) − α(Te − Ti ) + QL + Qr

(2.4)

Here, ptotal is defined as
ptotal = pe + pi

(2.5)

And, pT Hi and pT Hi are defined as,
∂pe
)ρ
∂Te
∂pi
= Ti ()ρ∂Ti

pT He = Te (

(2.6)

pT Hi

(2.7)

The ion and electron heat conductions in equation 2.3 and 2.4 are calculated simultaneously. In this calculation, flux-limited Spitzer-Harm model [21] with the flux-limiter of 0.1
is applied.

14

Laser absorption
The laser beam is injected along the meshes from the boundary of calculation area. On
the laser absorption process, inverse-bremsstrahlung is assumed as its mechanism [22]. Theenergy damping rates on each points, νabs , are given by equation (2.8). The nomenclature on
this part is shown in table 2.2.
Table 2.2: Nomenclature for the equation of laser absorption
νabs

laser damping ratekB

Boltzmann constantZ∗effective chargeTeelectron temperaturelnΛ

Coulomb logarithm [23]zeffective chargefLangdon factor [24]neelectron density

nenormalized electron density

ncritical

critical density of laser

λµm

laser wavelength<>

averaging over the ion spaciesνabs = 1.195 × 10−3Z ∗ lnΛf ne2
3

λ2µm (kB Te ) 2(2.8)

Here, Z ∗ is defined as,Z∗ =< z2 >
<z>

n
e=

ne

(2.9)

And, ne is defined as,
ncritical(2.10)

2.1.3 Hybrid code
The interaction between magnetic fields and plasmas is calculated by 3D hybrid code [16].
In this code, ions are treated as super-particles, and electrons are treated as an inertialesscold fluid. The charge density of plasma is also assumed as quasi-neutral. In addition, thecollisions between a plasma and coils are ignored. In this simulation, the number of superparticles is 106 , meaning 1010 ions are included in a super-particle, approximately.

15

basic equation
The nomenclature on this part is shown in table 2.3.
Table 2.3: Nomenclature for the equations of 3D hybrid codeneelectron densityZeffective chargeni

ion density

νeicollison frequencymeelectron massη

resistivity

mi

ion massj⃗i

ion current densityveelectron velocitykB

Boltzmann constantvi

ion velocity

Teelectron temperaturepeelectron pressureµ0

permeability of vacuum⃗Eelectric fieldTeelectron temperature⃗Bmagnetic field

The equations of motion for ions, particles, and electrons, fluid, are shown in equation (2.12) and equation (2.11).
dv⃗e= −ene (E + ve × B) + R − ∇pedtd⃗vi
n i mi= Zeni (E + vi × B) − Rdtn e me(2.11)(2.12)The term R on these equations are the friction term for the momentum of electron-ion
collision [25, 26]. This term is given by following equation.
R = ne me (vi − ve )νei
= ne ηe (ji − ene ve )(2.13)The term ∇pe on equation (2.11) is the electron pressure. This term is given by following

16equation.
∇pe = ∇(ne kB Te )= ne ηe(ji − ene ve )(2.14)

Here, the electrons are treated as a inertialess fluid. It means that the mass of electron
is assumed to be zero. From this assumption, the equation (2.11) is rewritten as following
equation (2.15)
0 = −ene (E + ve × B) − ∇pe(2.15)

Equation 2.17 is the Ampere’s law with Darwin approximation. This approximationneglects high frequency phenomena in a pasma such as electromagnetic radiations and waves
that are associated with electron inertia.
⃗ = µ0⃗j∇×B(2.16)= µ0 (j⃗i − ene v⃗e )

v⃗e =

⃗
1 ⃗ ∇×B
(ji −
)
ene
µ0

(2.17)

Here, we can get equation 2.18 to plug in the equation 2.17 to the equation 2.15.
(
)
}
⃗
1
∇
×
B
⃗+
⃗ − η⃗j + kB Te ∇ne
0 = ene E
j⃗i −
×B
ene
µ0
}
{(
)
⃗
1
∇
×
B
⃗ + kB Te Z∇ni + η⃗j
⃗ = −
×B
E
j⃗i −
Zeni
µ0
{(
)
}
⃗
1
∇
×
B
⃗ + kB Te Z∇ni + η ∇ × B
⃗
= −
j⃗i −
×B
Zeni
µ0
µ0
{

(2.18)

Using Maxwell-Faraday’s law as shown in equation 2.19, we can get equation 2.20.
⃗⃗ = − ∂B
∇×E
∂t

17

(2.19)

⃗
∂B
1
∇×
=
∂t
Zeni

{(

⃗
∇B
j⃗i −
µ0

)

}
⃗ + kB Te Z∇ni
×B

+

η 2⃗
∇B
µ0

(2.20)

Here, η and Te is assumed to be zero for simplicity. Thus, the equation 2.20 is rewritten
as the equation 2.21.
⃗
∂B
1
=
∇×
∂t
Zeni

{(

⃗
∇B
j⃗i −
µ0

)

}
⃗
×B

(2.21)

2.1.4 Coupling code
In this code, 2D data from radiation hydrodynamics code are extended to 3D data thatare usable on the 3D hybrid code. The 2D data are extended as axial symmetry to the laser
axis as shown in figure 2.2

Figure 2.2: The simulation model for the 2D radiation hydrodynamics code

The radiation hydrodynamics code calculates laser ablation process based on the CGS
system of units. On the other hand, the hybrid code calculates interaction between plasma
and magnetic field based on the MKS system of units. Thus, the coupling code converts CGS
into MKF system of units.
Secondly, the ion fluid is converted to super-particles. The number of super-particles ineach mesh are estimated based on the ion density of each mesh and the its volume as the18following equations.
N umber of ions in each mesh = (Ion density) × (mesh volume)(2.22)

N umber of total ions = Σ(N umber of ions in each mesh)(2.23)

N umber of ions/super-particle =N umber of super-particle in each mesh =N umber of total ions
N umber of super-particles(2.24)

N umber of ions in each mesh(2.25)
N umber of ions/super-particle

Here, the volume for each mesh is ring-like shape as shown in Fig. 2.3.

Figure 2.3: The simulation model for the 2D radiation hydrodynamics code

The volume is calculated from equation 2.26.
22)(zmax − zmin )− rminV olume = π(rmax(2.26)

Thirdly, the positions and velocities are converted to 3D Cartesian coordinate system.
The x and y components are decided with a random rotation angle around the z axis, the axis
of axial symmetry. Here, θn is the rotation angle and randn is a random number.
19θ1 = 2πrand1(2.27)

x = rmin + rand2 (rmax − rmin ) cos θ1(2.28)

y = rmin + rand2 (rmax − rmin ) sin θ1(2.29)z = zmin + rand3 (zmax − zmin )(2.30)

θ2 = 2πrand4(2.31)

Vx = Vr cos θ2(2.32)

Vy = Vr sin θ2(2.33)

Vz = Vz(2.34)

2.2 Simulation models
2.2.1 Replication of experimental result
Figure 2.4 shows the simulation model of the 2D radiation hydrodynamics code. The
simulation area is rectangle of ± 2.5mm × 3mm. Each axis is divided to 100 meshes. The
target is placed at r=z=0 . The laser beam comes from -z direction. The calculation area
except target is filled with the low density background gas because the code uses Eulerian
mesh.

20

Figure 2.4: The simulation model for the 2D radiation hydrodynamics code

Figure 2.5 shows the simulation model of the 3D hybrid code for the single solenoid coilsystem. The simulation area is a cube of the ranges of x, y, and z are ± 150mm, respectively.
Each axis is divided to 100 meshes equally. The target is placed at origin. The coil is placed
at 18 mm apart from the center. The coil is a 96 turn solenoid coil. The inner radius, outerradius, and length are 13 mm, 25 mm, and 10 mm, respectively. It is driven by 90 kA·turncurrent to produce 1.1 T magnetic field at the target position.

Figure 2.5: The simulation model for the 3D hybrid code with a single solenoid coil21

2.2.2 Thrust deflection
Figure 2.6 and 2.7 shows the setup image of multi coil system and the simulation modelof the 3D hybrid code for the multi coil system, respectively. he simulation area is a cubeof the ranges of x, y, and z are ± 150mm, respectively. Each axis is divided to 100 meshes
equally. The target is placed at origin. The coils are placed at 40 mm apart from the center.
The multi coil system consists of four coils. All coils have same 42 mm square shape andplaced toward the laser axis. Each coil is 10 turn solenoid coil. These coils are controlled
independently, 0∼50 kA·turn.

Figure 2.6: The setup image of multi coil system

Figure 2.7: The simulation model for the 3D hybrid code with multi coil system

22

Chapter 3
Results
3.1 Replication of experimental result
At first, I calculated using the existing codes to find the problems. In the previous method,
the radiation hydrodynamics code and the hybrid code were connected at 75 ns. The time
evolution of laser intensity is shown in Figure 3.1. As shown in this figure, the intensity has apeak at 25 ns. Thus, the coupling timing is after 50 ns from the peak timing of laser intensity.
This time was decided based on the heat relaxation [27] because the plasma temperatures are
ignored in the hybrid code calculation. Figure 3.2 shows the time evolutions of radiation,internal, and kinetic energy. The plasma temperatures should be low enough on the coupling
timing. In other words, kinetic energy should be dominant compared to thermal energy.

23

Figure 3.1: Time evolution of the laser intensity

Figure 3.2: Time evolution of the energy balance [27]Figure 3.3 shows the density distribution at the coupling timing (t=75 ns). This figure
shows that the high density area is still remaining in the calculation area. On the other hand,24

Figure 3.4 shows the velocity distribution at the coupling timing. The previous experimental
result, Figure 1.8, shows that there were high speed ions over 150 km/s. However, there are
only low speed ions up to 70 km/s on the numerical result.

Figure 3.3: Ion density distribution at 75 ns

Figure 3.4: Plasma velocity distribution at 75 ns

25

To solve this problem, I checked the time evolution of density and velocity distributionof plasma. Then, the high speed plasma reaches the boundary of calculation area at 33
ns. Figure 3.5 and 3.6 shows the ion density and velocity distributions at 33 ns, respectively.
Thus, the coupling time of 75 ns is too late to include fast ions (v ≤ 250 km/s), and this is oneof the reasons why the previous simulation did not express the experimental ion distribution
early in time.

Figure 3.5: Ion density distribution at 33 ns

Figure 3.6: Plasma velocity distribution at 33 ns
26

Figure 3.7 shows the time evolution of the ion current densities that are measured on theexperiment and the numerical simulation in which the radiation hydrodynamics and hybrid
codes are connected at 33 ns. On this result, the high speed ions appear early in time att ∼1-2µs as the experimental result indicated. However, the amount of ion current density isquite low.

Figure 3.7: Comparison of numerical result with experimental one
I found this difference comes from the method how to allocate super-particles with thecoupling code. In the coupling code, the super-particles are allocated based on the massdensity of each mesh. The ion density of center area of target is quite high compared withouter area. Thus, most super-particles are allocated on center area of the target. Therefore,
the number of super-particles in the outer area is small, meaning the number of high-speedions is too small.
To solve this problem, I installed the solid removing scheme to the coupling process.
Here, the threshold value between solid and plasma is set to 1027 m−3 . This value is the 10% of the original solid density of the target. As a result, the high density area was removed
as figure 3.8.
27

Figure 3.8: Ion density distribution from the improved coupling code
Figure 3.9 shows the comparison of ion current density from numerical and experimental
results using the solid-removed data as the initial plasma condition of 3D hybrid code. As canbe seen in this figure, two peaks as seen in the experimental result appeared on the numericalresult. In addition, the time scale is almost same. The amount of ion current density isfive times larger than the experimental result on the first peak. However, the discrepancy is
reduced significantly from previous calculations.

28

Figure 3.9: Comparison of numerical result with experimental one
Figures 3.10 and 3.11 are the super-particle and ion density distributions at 1µs in which
the initial plasma for the hybrid code is Fig. 3.5 and 3.6. On the other hand, Figs. 3.12 and3.13 are the solid-removed case, Fig.3.8. As can be seen in these figures, the amount ofhigh-speed super-particles is increased significantly.

29

Figure 3.10: Super-particle distribution at 1 µs Figure 3.11: Ion density distribution at 1 µs

Figure 3.12: Super-particle distribution at 1 µs Figure 3.13: Ion density distribution at 1 µs

Figures 3.14 shows the comparison of a the self-emission from the experiment and theion density distribution from the simulation with the magnetic field strength of 1.1 T at t =0.5 µs. These figure are comparable because the intensity of self emission is proportional to
the square of plasma density. This result shows that plasma expansion in a magnetic field ismore or less the same in both simulation and experiment.
30

Figure 3.14: Comparison with (a)self-emission and (b)Ion density distribution at 0.5 µs

31

3.2 Thrust deflection
Here, I simulated with three different magnetic field configurations. The initial conditions of the magnetic field structure is shown in Fig. 3.15 and 3.16. First, case(a), all coils
are driven by 50 kA-turn. Second, case(b), three coils on the ± y and -x side are driven by
50 kA·turn, and the coil on +x side is driven by 25 kA·turn.Third, case(c), three coils on the± y and -x side are driven by 50 kA·turn, and the coil on +x side is off (zero current). As

can be seen on these figures, the magnetic field strengths of the +x side become weaker from(a) to (c).

Figure 3.15: Initial magnetic fields on x-y plane at z=0(a) All coils are driven by 50 kA·turn(b) Three coils except on +x side by 50 kA·turn, and on +x side is driven by 25 kA·turn(c) Three coils except on +x side by 50 kA·turn, and on +x side is off (zero current)·turn32

Figure 3.16: Initial magnetic fields on x-z plane at y=0(a) All coils are driven by 5 kA(b) Three coils except on +x side by 50 kA·turn, and on +x side is driven by 25 kA·turn(c) Three coils except on +x side by 50 kA·turn, and on +x side is off (zero current)
The super-particle distributions and the density distributions are shown in Figs. 3.17 and
3.18, respectively. As shown is these figures, the deflection angle become greater with decreasing magnetic field from (a) to (c). In addition, the plasma flow that goes through the-x side coil increases from (a) to (c). It shows that the plasma flows propagate along themagnetic field.

Figure 3.17: Super-particle distributions(a) All coils are driven by 5 kA(b) Three coils except on +x side by 50 kA·turn, and on +x side is driven by 25 kA·turn(c) Three coils except on +x side by 50 kA·turn, and on +x side is off (zero current)33

Figure 3.18: Ion density distributions(a) All coils are driven by 5 kA(b) Three coils except on +x side by 50 kA·turn, and on +x side is driven by 25 kA·turn(c) Three coils except on +x side by 50 kA·turn, and on +x side is off (zero current)
Figures 3.20 and 3.21 show the time evolution of thrust deflection angle. The angle isdefined as Fig.3.19. As shown in these figures, the thrust angle on the initial condition is 180
degrees on every case and it decrease rapidly. It shows that the initial plasma flow is headed
to upstream of the laser beam, and the plasma flow is changed its direction by the interaction
with the magnetic field. After that, the thrust angle becomes almost constant at 10 µs. The
maximum deflection angle was 3 degrees on the case (c).

Figure 3.19: The definition of a deflection angle

34

Figure 3.20: Thrust deflection angle

Figure 3.21: Thrust deflection angle, 8-10 ns

35

Chapter 4
Discussion
4.1 Replication of experimental result
The final result still shows differences in the absolute value of ion current density andin the timing of the peak compared with experimental result. The differences are shown in
Fig. 4.1

Figure 4.1: The differences between the experimental and the numerical result
36

(A) is the temporal difference. The waveform from numerical simulation has delay of∼0.4µs from the experimental result on the first peak, approximately. The first peak on theexperimental result and the numerical one is 150 km/s and 107 km/s, approximately. I will
itemize the possible reasons for this difference.
Part of laser energy is ignored
The laser pulse is gaussian shape in time. The laser intensity on the coupling
timing, 33 ns, is 35.5 M W/cm2 as shown in Fig. 4.2. It is 13.4 % of peak intensity.

Figure 4.2: The time evolution of the laser intensity around 33 ns

The FWHM of this pulse is 9.4 ns. The σ of this gaussian pulse is,
F W HM√
2 2 ln 2≈ 4 ns

σ =

(4.1)
(4.2)

Thus, 33 ns is the point of +2σ. It means that 2.5 % of laser energy has not incident
at the coupling timing.
To solve this problem, longer calculation using radiation hydrodynamics code isrequired. However, the plasma will expand without the interaction with a external
magnetic field in the radiation hydrodynamics simulation. It will be an another problem.
37

Low ion density regions are ignored in the coupling process
In the coupling process, the ions are converted to super-particles. A super-particle
contain 1010 ions. When a number of ions in a mesh is less than one super-particles, theions are ignored. In this case, fast ions may be ignored and ion current density maybe
underestimated at t ∼ 1 µsṗar

To solve this problem, increasing the number of thesuper-particles is required. It is directly linked to increasing the calculation cost. Thus,
speed-up or/and parallelization of the simulation code are required simultaneously.
Temperatures are ignored
The plasma temperatures, ion and electron temperatures, are ignored to simplify
the calculation process in the hybrid code. Nevertheless, the laser ablation plasma is
not a cold plasma at 33 ns. The temperature distributions are shown in Fig. 4.3.

Figure 4.3: The distributions of the plasma temperatures

Here, the energy ratio between kinetic and thermal is calculated as following equations. The velocity of ion is assumed as 100 km/s. The mass of ion is assumed as
6.5 mproton , averaged value of Hydrogen and Carbon.
12
mion vion2(4.3)≈ 5.44 × 10−17 J(4.4)≈ 340 eV

(4.5)

Kinetic energy of ion =38

The thermal energy is up to 50 eV approximately. When the thermal and kinetic
energy of an ion is 340 eV and 50 eV, the thermal energy is around 12 %. If thethermal energy is converted to kinetic energy in a magnetic thrust chamber, the ion get
additional 38 km/s. It is not negligible value.
To solve this problem, installation of the plasma temperatures in the hybrid code isrequired.
The external magnetic field is ignored in the radiation hydrodynamics code
The magnetic field is ignored in the radiation hydrodynamics code. On the other
hand, there are quasistationary external magnetic field on the experimental condition. Thus, to replicate the experimental condition accurately, installation of an external magnetic field in the simulation of a laser ablation process is required, meaningradiation-MHD simulation is required.

The velocity difference between a proton and a carbon ion
Current simulation codes can treat only one material. Thus, the averaged parameters of hydrogen and carbon are used as the parameters of ion. However, there should
be differences between hydrogen and carbon plasma. If a hydrogen plasma get higher
velocity than a carbon plasma, the first peak consist of hydrogen plasma and it should
appear earlier time.
To solve this problem, simulation with two materials is required.

39

(B) is the difference of magnitude of the ion current density on the first peak. The ion
current density on the numerical result is 4 - 5 times larger than the experiment. One of thepossible reason is the incidence angle to the charge corrector. In the numerical simulation,
the ion current density is calculated from the ion flow rate at the inspection surface without
considering the incidence angle of ions. However, the actual charge corrector is covered bymetal cap as shown in Fig. 4.4. In this condition, ions that have large incidence angle cannot
reach the Faraday cup. Thus, the numerical value of ion current density should be larger thanexperimental one. To solve this problem, the filtering process by the incidence angle of ion
is required on the numerical simulation.

Figure 4.4: The incidence angle of ion for a charge corrector(C) is the difference of magnitude of the ion current density on the second peak and later.
The time evolution of ion current density on the numerical result has long tail compared withexperiment. One of the possible reason is electron-ion re-combination. Re-combination
processes is not consider in the hybrid code and charge collectors cannot measure the neutral
particles. To solve this problem, to validate the possibility of recombination in the plasma
conditions and to install the recombination processes to the hybrid code are required.

40

4.2 Thrust deflection
The coil conditions and the thrust angles are summarized in table 4.1.
Table 4.1: The results of the thrust deflection
Driven current of coils [kA-turn] Deflection angle [deg.]
+x

-x

+y

-y

at 10 µs

50

50

50

50

0.8

25

50

50

50

1.2

0

50

50

50

3.0

The maximum deflection angle is 3 degrees on the case (c). It is smaller than the previous
result [17], 5.2 degrees. To increase the deflection angle, increasing the ratio of the magnetic
field on the upper (+x) and lower (-x) side is required. Increasing current is the simplest way
to make high ratio but it is directly linked to increasing mass of power supply systems. Here,I propose the 45 degrees rotation setup of multi coil system. Figure 4.5 shows the layout ofmultiple coils. In this layout, only one coil is off to deflect plasma flow to +x side. On theother hand, Fig. 4.6 shows the proposed layout of multiple coils. In this layout, two coils are
off to deflect plasma flow to +x side. Thus, the ratio of the magnetic fields on the upper andlower side should increase.

Figure 4.6: Proposed multi coil system

Figure 4.5: Current multi coil system

Ideally, the deflection angle of case (a) should be zero because the initial conditions of41

plasma and magnetic field are completely plane symmetry to the y-z plane. However, theangle is not zero. To solve this problem, I checked the initial conditions of the plasma. As
a result, the initial momentum of plasma is little biased to +x and -y side. It is caused by
the coupling process. In the coupling process, 2D data are extended to 3D data as axial
symmetry. Then, the angles around the symmetrical axis are set randomly. This random
process is considered the reason of the bias and the bias is a possible reason of deflection on
the case (a).
From this result, the thrust deflection system with 4 coils worked well numerically.
As a next step, it is required to validate the efficiency of this system as an application to
LFR.

42

Chapter 5
Summary
The high-speed ions and the two peaks as seen in the experimental result are appearedon the numerical result by the improvements. The time scale is almost same but the waveform from numerical simulation has delay of ∼0.4µs from the experimental result on the firstpeak, approximately. In addition, the amount of ion current density is five times larger thanthe experimental result on the first peak. There are long tail after second peak. These differences are come from the approximation in the radiation hydrodynamics code and/or hybrid
code. To improve the reproducibility of the experimental result with the numerical simulation and to understand the physical mechanism of the magnetic thrust chamber, improvementof numerical simulation is required: validation of applicability of simulation procedure withcoupling process, installing magnetic field to radiation hydrodynamics code, and installing
temperatures to hybrid code.
The thrust deflection system with four coils worked well numerically. The maximum
deflection angle was 3 degrees. It is smaller than the previous result: 5.2 degree [17]. It is required to optimize this system such as number of coils, coil position, and current to optimize
the thrust deflection system. In addition, experimental validation is required simultaneously.
On the experiment, self-emission observation, arrayed charge correctors, and thrust stand areviewed as way to measure the deflection angle.

43

List of Figures
1.1 Launch of the first Epsilon Launch Vehicle (Epsilon-1)[4] . . . . . . . . . .
1.2 Deep Space 1 with the NSTAR Ion Engine[5] . . . . . . . . . . . . . . . .
1.3 Flight model image of VISTA . . . . . . . . . . . . . . . . . . . . . . . .
1.4 The mechanism of a magnetic thrust chamber . . . . . . . . . . . . . . . .
1.5 Concept image of multi coil system . . . . . . . . . . . . . . . . . . . . .
1.6 The image of the mission phase of LFR . . . . . . . . . . . . . . . . . . .
1.7 Setup of EUV experiment . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.8 Comparison of numerical result with experimental one . . . . . . . . . . .
1.9 Comparison with (a)self-emission and (b)Ion density distribution at 0.5 µs .
1.10 The ion flow behaviors with multi coil system [17] . . . . . . . . . . . . .
1.11 Previous results of thrust deflection system [17] . . . . . . . . . . . . . . .

2
2
3
4
5
5
7
8
9
10
10

2.1
2.2
2.3
2.4
2.5
2.6
2.7

Numerical simulation procedure . . . . . . . . . . . . . . . . . . . . . . .
The simulation model for the 2D radiation hydrodynamics code . . . . . .
The simulation model for the 2D radiation hydrodynamics code . . . . . .
The simulation model for the 2D radiation hydrodynamics code . . . . . .
The simulation model for the 3D hybrid code with a single solenoid coil . .
The setup image of multi coil system . . . . . . . . . . . . . . . . . . . . .
The simulation model for the 3D hybrid code with multi coil system . . . .

12
18
19
21
21
22
22

3.1 Time evolution of the laser intensity . . . . . . . . . . . . . . . . . . . . .
3.2 Time evolution of the energy balance [27] . . . . . . . . . . . . . . . . . .
3.3 Ion density distribution at 75 ns . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Plasma velocity distribution at 75 ns . . . . . . . . . . . . . . . . . . . . .
3.5 Ion density distribution at 33 ns . . . . . . . . . . . . . . . . . . . . . . . .
3.6 Plasma velocity distribution at 33 ns . . . . . . . . . . . . . . . . . . . . .
3.7 Comparison of numerical result with experimental one . . . . . . . . . . .
3.8 Ion density distribution from the improved coupling code . . . . . . . . . .
3.9 Comparison of numerical result with experimental one . . . . . . . . . . .
3.10 Super-particle distribution at 1 µs . . . . . . . . . . . . . . . . . . . . . . .
3.11 Ion density distribution at 1 µs . . . . . . . . . . . . . . . . . . . . . . . .
3.12 Super-particle distribution at 1 µs . . . . . . . . . . . . . . . . . . . . . . .
3.13 Ion density distribution at 1 µs . . . . . . . . . . . . . . . . . . . . . . . .
3.14 Comparison with (a)self-emission and (b)Ion density distribution at 0.5 µs .

24
24
25
25
26
26
27
28
29
30
30
30
30
31

44

3.15 Initial magnetic fields on x-y plane at z=0(a) All coils are driven by 50 kA·turn(b) Three coils except on +x side by 50 kA·turn, and on +x side is driven by
25 kA·turn(c) Three coils except on +x side by 50 kA·turn, and on +x side is off (zerocurrent)·turn . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.16 Initial magnetic fields on x-z plane at y=0(a) All coils are driven by 5 kA(b) Three coils except on +x side by 50 kA·turn, and on +x side is driven by
25 kA·turn(c) Three coils except on +x side by 50 kA·turn, and on +x side is off (zerocurrent) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.17 Super-particle distributions(a) All coils are driven by 5 kA(b) Three coils except on +x side by 50 kA·turn, and on +x side is driven by
25 kA·turn(c) Three coils except on +x side by 50 kA·turn, and on +x side is off (zerocurrent) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.18 Ion density distributions(a) All coils are driven by 5 kA(b) Three coils except on +x side by 50 kA·turn, and on +x side is driven by
25 kA·turn(c) Three coils except on +x side by 50 kA·turn, and on +x side is off (zerocurrent) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.19 The definition of a deflection angle . . . . . . . . . . . . . . . . . . . . . .
3.20 Thrust deflection angle . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.21 Thrust deflection angle, 8-10 ns . . . . . . . . . . . . . . . . . . . . . . . .

34
34
35
35

4.1
4.2
4.3
4.4
4.5
4.6

36
37
38
40
41
41

The differences between the experimental and the numerical result . . . . .
The time evolution of the laser intensity around 33 ns . . . . . . . . . . . .
The distributions of the plasma temperatures . . . . . . . . . . . . . . . . .
The incidence angle of ion for a charge corrector . . . . . . . . . . . . . .
Current multi coil system . . . . . . . . . . . . . . . . . . . . . . . . . . .
Proposed multi coil system . . . . . . . . . . . . . . . . . . . . . . . . . .

45

32

33

33

List of Tables
1.1

Experimental parameters . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.1
2.2
2.3

Nomenclature for the equations of radiation hydrodynamics code . . . . . .
Nomenclature for the equation of laser absorption . . . . . . . . . . . . . .
Nomenclature for the equations of 3D hybrid code . . . . . . . . . . . . . .

14
15
16

4.1

The results of the thrust deflection . . . . . . . . . . . . . . . . . . . . . .

41

46

Bibliography[1] WR Webber and FB McDonald. Recent voyager 1 data indicate that on 25 august 2012
at a distance of 121.7 au from the sun, sudden and unprecedented intensity changes
were observed in anomalous and galactic cosmic rays. Geophysical Research Letters,
40(9):1665–1668, 2013.
[2] B Lyndon. Apollo program summary report. Johnson Space Center, NASA, TM X,
68725, 1975.
[3] Dan M Goebel and Ira Katz. Fundamentals of electric propulsion: ion and Hall
thrusters, volume 1. John Wiley & Sons, 2008.
[4] JAXA Digital Archives. Launch of the first epsilon launch vehicle (epsilon1). http://jda-strm.tksc.jaxa.jp/archive/photo/P100008529/
ba76785908182f9302a4df6138fa10ac.jpg.
[5] NASA Image and Video Library.
Deep space 1 ion engine.
https://
images-assets.nasa.gov/image/PIA04247/PIA04247˜orig.jpg.
[6] National Aeronautics and Space Administration. Nasa’s journey to mars : pioneering
next steps in space exploration. NP-2015-08-2018-HQ, 2015.
[7] NASA News Releases. Radiation measured by nasa’s curiosity on voyage to mars has
implications for future human missions, 2013. https://www.nasa.gov/home/
hqnews/2013/may/HQ_13-165_MSL_Radiation_Findings.html.
[8] CD Orth et al. Vista-a vehicle for interplanetary space transport application poweredby inertial confinement fusion. University of California, Lawrence Livermore National
Laboratory, Report UCRL-LR-110500, in final preparation, 2003.
[9] N Saito, R Kawashima, T Morita, N Yamamoto, H Nakashima, S Fujioka, A Yogo,
H Nishimura, A Sunahara, Y Mori, and T Johzaki. Optimizing a magnetic field strengthof a magnetic thrust chamber. In 30th International Symposium on Space Technology
and Science, 34th International Electric Propulsion Conference, number IEPC-2015465p/ISTS-2015-b-465p in Conference paper, Kobe, Japan, 2015.
[10] R Kawashima, T Morita, N Yamamoto, N Saito, S Fujioka, H Nishimura, H Matsukuma, A Sunahara, Y Mori, T Johzaki, and H Nakashima. The measurementof plasma structure in a magnetic thrust chamber. Plasma and Fusion Research,
11(3406012), 2016.

47

[11] T Morita, N Yamamoto, R Kawashima, N Saito, M Edamoto, S Fujioka, Y Itadani,
T Johzaki, S Miura, Y Mori, et al. Plasma structure and energy dependence in a magnetic thrust chamber system. Journal of Physics: Conference Series, 717(1):012071,
2016.
[12] M Yasunaga, A Maeno, N Yamamoto, H Nakashoma, S Fujioka, A Sunahara, andJohzaki T. Experimental demonstration of magnetic thrust chamber for a laser fusion rocket. Transactions of the Japan Society for Aeronautical and Space Sciences,
Aerospace Technology Japan, 10(ists28):Pb 109–Pb 114, 2012.
[13] A Maeno, N Yamamoto, S Fujioka, A Sunahara, T Johzaki, and H Nakashima. Analysisof laser wavelength and energy dependences of the impulse in a magnetic thrust chamber system for a laser fusion rocket. Transactions of the Japan Society for Aeronauticaland Space Sciences, 56(3):170–172, 2013.
[14] A Sunahara, K Nishihara, and A Sasaki. Optimization of extreme ultraviolet emission from laser-produced tin plasmas based on radiation hydrodynamics simulations.
Plasma and Fusion Research, 3:043–043, 2008.
[15] A Sunahara, T Johzaki, H Nagatomo, and K Mima. Generation of pre-formed plasma
and its reduction for fast-ignition. Laser and particle beams, 30(01):95–102, 2012.
[16] Y Nagamine and H Nakashima. Analysis of plasma behavior in a magnetic thrustchamber of a laser fusion rocket. Fusion Science and Technology, 35(1):62–70, 1999.
[17] Y Kajimura, R Kawabuchi, N Matsuda, and H Nakashima. Numerical simulation offusion plasma behavior in a magnetic nozzle for laser fusion rocket. Fusion Scienceand Technology, 51(2T):229–231, 2007.
[18] L Spitzer Jr. Physics of fully ionized gases:Second Revised Edition. Dover Publications,
2006.
[19] L Spitzer Jr and R Härm. Transport phenomena in a completely ionized gas. PhysicalReview, 89(5):977, 1953.
[20] SI Braginskii. Transport processes in a plasma. Reviews of plasma physics, 1:205,
1965.
[21] R. C. Malone, R. L. McCrory, and R. L. Morse. Indications of strongly flux-limited
electron thermal conduction in laser-target experiments. Physical Review Letters,
34(12):721, 1975.
[22] W. L. Kruer. The physics of laser plasma interactions. Addison-Wesley PublishingCo., 1988.
[23] S Skupsky. coulomb logarithm for inverse-bremsstrahlung laser absorption. Physical Review A, 36(12):5701, 1987.
[24] A B. Langdon. Nonlinear inverse bremsstrahlung and heated-electron distributions.
Physical Review Letters, 44(9):575, 1980.

48

[25] A Mankofsky, R. N. Sudan, and J Denavit. Hybrid simulation of ion beams in background plasma. Journal of Computational Physics, 70(1):89–116, 1987.
[26] T Bagdonat and U Motschmann. 3d hybrid simulation code using curvilinear coordinates. Journal of Computational Physics, 183(2):470–485, 2002.
[27] A Maeno. Study on Experimental Demonstration and Numerical Simulation of Magnetic Thrust Chamber System. PhD thesis, Kyushu University, 2014.

49

Acknowledgments
I would like to express the deepest appreciation to my supervisors; Professor Naoji Yamamoto, and Assistant Professor Taichi Morita, and Professor Hideki Nakashima.
I would like to express the appreciation to my advisers and collaborative researchers; Dr.
Atsushi Sunahara of Institute for Laser Technology, Professor Shinsuke Fujioka of Osaka
University, Dr. Yoshitaka Mori of The Graduate School for the Creation of New Photonics
Industries, Dr. Tomoyuki Johzaki of Hiroshima University, Dr. Akifumi Yogo of Osaka University, and Professor Hiroaki Nishimura of Osaka University. I would also like to express
the appreciation to Mr. Eiji Sato and GEKKO technical crews for their exceptional support
during the experiments.
I would like to appreciate Mrs. Takako Mawatari who take care of my paperwork.
Advice and comments given by Akihiro Maeno has been a great help in my research.
I would like to thank the LFR group; Naoya Saito, Ryosuke Kawashima, Satoshi Miura,
Yutaro Itadani, and Tomihiko Kojima. In addition, I thank my labmates; Lee Hou-Yi, Koichi
Ushio, Yu Okuda, Kohei Takase, Kazuhiko Nakano, Atsushi Yamaguchi, Kensuke Iijima,
Tomohiro Ichimaru, Bunsuke Ueno, Yusuke Egawa, Masataka Iwamoto, Yutaro Kawahara,
Hirotaka Fuchigami, Yu Fukumoto, Tomoki Koyama, and Ippei Takesue.
This work was supported by the Japan Society for the Promotion of Science (JSPS)
KAKENHI (Grant Numbers 25420852 and 15K18283) and by the joint research projectof the Institute of Laser Engineering, Osaka University (2014B1-31, 2015B1-15, 2016A1YAMAMOTO).
The computation was mainly carried out using the computer facilities at Reaearch Institute for Information Technology, Kyushu University.
Finally, I would like to express the deepest gratitude to my family for their kind supportand help.

50