Fusion Engineering and Design 44 Ž1999. 163]169
2D simulation of hydrodynamic instability in ICF stagnation phase A. SunaharaU , H. Takabe, K. Mima Institute of Laser Engineering, Osaka Uni¨ ersity, Osaka 565, Japan
Abstract In the stagnation phase of inertial confinement fusion, hydrodynamic instabilities prevent fuel heating, and fusion neutron yield is reduced. We calculated the neutron yield, fuel temperature, and the energy transport in the stagnation phase. The measured neutron yield of the implosion experiment carried out at this institute ŽILE. is compared with the calculated value. The trend of calculation values agrees with that of the experiment. Q 1999 Elsevier Science S.A. All rights reserved. Keywords: Inertial confinement fusion; Hydrodynamic instability
1. Introduction In inertial confinement fusion ŽICF., a capsule containing deuterium and tritium is imploded by the high pressure induced by laser irradiation, and high temperature and high density are achieved to ignite thermonuclear reaction in the compressed fuel. This scheme is called central ignition. The process to achieve the central ignition can be divided into four phases; the acceleration, inertial, stagnation, and expansion phases. In the stagnation phase, the pusher is decelerated by the fuel pressure, and the kinetic energy of the fuel and pusher is mainly converted into the fuel
U
Corresponding author.
internal energy. Through this stagnation heating, the thermonuclear burn begins. However, in the experiments carried out so far, the measured neutron yields are lower by an order of magnitude than the values calculated using the one dimensional implosion codes, such as ILESTA ] 1D w1x. This discrepancy seems mostly due to growth of the hydrodynamic instability seeded by non-uniformities of laser irradiation and targets w2x. Such reduction of implosion performance caused by the hydrodynamic instability is the most important issue to be studied in inertial confinement fusion. The hydrodynamic instability affects the energy transport in the fuel corresponding to the heating or cooling processes of the fuel. Therefore, we evaluated effects of the hydrodynamic instability on the energy transport and the implo-
0920-3796r99r$ - see front matter Q 1999 Elsevier Science S.A. All rights reserved. P I I: S 0 9 2 0 - 3 7 9 6 Ž 9 8 . 0 0 3 6 2 - 7
A. Sunahara et al. r Fusion Engineering and Design 44 (1999) 163]169
164
sion performance in the stagnation phase, and compared the calculated value with the experimental values in order to investigate a physical mechanism of the reduction of neutron yield. 2. Numerical scheme and simulation condition As an initial attempt to understand the time evolution of the implosion process, we have carried out two-dimensional simulations using the second order Godunov scheme w3x, which solves the shock tube problem between two adjacent zones at each time step using the Riemann solver. We assume that the fuel is an ideal fluid with a constant specific heat ratio of 5r3. Then, Euler equations in the conservative form are given by
U 1 xF x F y q q s 0, t x x y r ru Us , F xs r¨ rE
ru
Ž 1a .
r¨ r ¨u
pq r u , F ys , r u¨ pq r ¨ 2 r uEq up r ¨E q¨p Ž 1b .
0 0 0 2
and Es
p q 1 Ž u2 q¨ 2 . , Ž g y 1. ? r 2
Ž 1c .
where r is the density, and u and ¨ are the x and y components of the flow velocity, respectively. The total energy density and pressure are denoted by E and p, respectively. The density vector U consists of the following physical variables; the density r , the x and y components of the momentum density r u and r ¨ and the energy density r E. The x and y components of the physical flux are given by F x and F y , respectively. In the present calculation, the physical quantities are assumed to have linear profiles in each zone, namely, the so-called ‘piecewise linear approximation’ is adopted. Consequently, the numerical accuracy of calculation is of the second order in space. In evaluating the physical variables at each mesh, we used the van-Leer’s flux limiter w4x in
order to guarantee monotonicity of the variables. Euler equations ŽEqs. 1a]1c. are solved iteratively by using Newton’s method for the D tr2 time step. The calculation is carried out in the xand y-direction separately. This is called ‘operator splitting’ and as a result the calculation becomes essentially one dimensional. The time derivative is calculated by the center finite difference scheme, and numerical accuracy of calculation is of second order in time. In the present simulation, we assume that the ion and electron temperatures are the same, since the temperature relaxation time is short enough compared to the typical time scale of the dynamics. The temperature T is evaluated by equation of state for an ideal gas, Ts
mi p , Ž Z q 1. r
Ž2.
where Z and m i are the charge number and the mass of the ion. In the present calculation, we assume plasmas with Z s 1 and m i s 1.45m p for hydrogen-deuterium ŽH 3 D 2 . fuel and Z s 3.5 and m i s 7.5m p for CH Žpusher., where m p is the proton mass. We traced the fuelrpusher contact surface at every time step and distinguished the materials by the surface tracking code w5x. In the numerical calculation, the temperature is evaluated at every time step using Eq. Ž2.. Then, the equation for thermal conduction given by
r Cn
T s = ? Ž k =T . , t
Ž3.
is solved. In Eq. Ž3., C¨ is the heat capacity given by C¨y1 s Žg y 1. m i and the thermal conductivity k is assumed to be that of an ideal plasma given in Ref. w6x,
k s 15.8 te s 0.3
Z n e Tete , Z q 4 me
3r2 m1r2 e Te e 4 ² Z 2 : n i lnL
Ž 4a . Ž 4b .
where n e and n i are electron and ion number densities given to be n e s Zn i s Z rrm i . In Eq.
A. Sunahara et al. r Fusion Engineering and Design 44 (1999) 163]169
Ž4a. and Ž4b., lnL is the Coulomb logarithm given, for example, in Ref. w6x, and the electron temperature Te is equal to T of Eq. Ž2.. In numerically solving Eq. Ž3., we used the ADI algorithm. After solving Eq. Ž3., the pressure p is modified using Eq. Ž2. with the resultant temperature T. Then the total energy E given by Eq. Ž1c. is evaluated again using this pressure. In order to evaluate the radiation cooling and transport effect of the compressed targets, we have solved the equation for radiation transport w7x.
r
d dt
En q = ? w F n y uE n x s 4ph n y c x ¨ E n r
ž /
Ž5.
where n is the photon frequency and En and Fn are the photon energy density and flux, respectively. Emissivity hn and opacity xn which depend on density and temperature, are determined from a table w8x. The c is the speed of light. We assume that the photon energy flux is defined by Fnsy
c
=E n 3x q< n < E n
=E n
Ž6.
165
which is based on the flux limited diffusion approximation model w7x. We calculate transport of the photon energy region from 50 eV to 10 keV by dividing into 100 energy bins. After carrying out energy transport calculation, we correct the internal energy r E of the fluid cells. The 2D simulation starts at 1.93 ns when the shock wave collides with the center of the target. We adopted the simulation condition corresponding to the present implosion experiment Žshot a19714. as shown in Fig. 1. We obtain the initial condition at 1.93 ns from the ILESTA ] 1D, and perturb the fluid velocity as the 2D initial condition. 3. Comparison of calculated neutron yield and ion temperature with experiment 3.1. Neutron yield First, we calculated the deuterium-deuterium ŽDD. reaction neutron yields from the fuel for various velocity perturbation amplitudes of mode l , l s 0 Ž1D., 2, 6, and 12, as shown in Fig. 2, using the thermonuclear reactivity formula of Ref. w9x, where mode number l is perturbation wave
Fig. 1. Initial condition for shot a19714.
166
A. Sunahara et al. r Fusion Engineering and Design 44 (1999) 163]169
Fig. 2. The DD neutron yield for the various velocity perturbation amplitudes at mode number 0 Ž1D., 2, 6, and 12. Experimental value is 7.95= 10 5, with the systematic error bar of "35% Žshaded region..
number divided by 2p . The neutron yield decreases rapidly with increase of the mode number from 0 to 2, and then increases with increase of the mode number. The dependence of the neutron yield on the perturbation amplitude is significant near the mode number 2. However, at larger mode numbers, the neutron yield becomes less sensitive to the perturbation amplitude. At the mode number 12, the neutron yield becomes larger than that of 1D. This is because at higher modes with shorter wavelength the spikes generated by the hydrodynamic instability confine the fuel and the temperature of the compressed fuel increases. As a result, neutron yield increases. The measured DD neutron yield is 7.95= 10 5 with "35% systematic error w10x. However, dominant perturbation mode number is not determined experimentally. A calculation of laser irradiation non-uniformity suggests that the mode 6 is dominant w11x. We compared the calculated value at mode number 6 with the experiment. If the velocity perturbation amplitude is 20]30% at mode 6, the calculated neutron yield agrees with the measurement. In this paper, we considered the hydrodynamic instability caused by velocity perturbation only. However, in the present experiment, contact surface displacement is measured.
Fig. 3. Calculated ion temperature, Ti TOF ŽeV., of the fuel weighted in the neutron yield. Experimental value is 807 eV with the statistical error bar of "7% Žshaded region..
We will calculate the hydrodynamic instability caused by contact surface displacement in future. 3.2. Ion temperature The ion temperature of the fuel weighted by the neutron emission rate, which corresponds to the measured ion temperature by the neutron time of flight method w10x, is shown in Fig. 3. In the figure, it is seen that the temperature is minimum near the mode number 2, as observed in the neutron yield vs. mode number curves in Fig. 2. However, when the mode number is larger than 6, the ion temperature increases significantly as the velocity perturbation increases. This is in contrast to the case of the neutron yield in Fig. 2 where the larger perturbation tends to the decrease of the neutron yield and less sensitive at higher mode numbers. This can explain the experimental result w10x that inspite of decrease of the neutron yield, the ion temperature increased. The measured temperature is 807 eV with "7% statistical error w10x. From Figs. 2 and 3, both the neutron yield and the ion temperature measured by the experiment w10x agree with the calculation results if the mode number 6 is dominant. From the above investigation, we conclude that in the present implosion experiment, the pertur-
A. Sunahara et al. r Fusion Engineering and Design 44 (1999) 163]169
167
bation mode number 6 is dominant, and the implosion dynamics can be explained by the present simulation if we set the mode number of 6 with the velocity perturbation 20]30% is assumed as the initial conditions for the stagnation phase. The spike grown by the hydrodynamic instability confines the fuel and increases the temperature at the center region of the fuel. 3.3. Effect of multi modes The previous study has been carried out for the case of a single mode. In order to investigate effects of multi-modes on both the neutron yield and the ion temperature, we repeated the simulation by assuming multi-modes Žmode 6, 2 and 12. with different amplitude ratios. In Fig. 4, the horizontal axis shows the ratio of mode number 6 to mode number 2 Žin the left., and to mode number 12 Žin the right half.. The full line represents the neutron yield Žthe value at the left. and the dotted line the ion temperature Žthe value at the right.. The perturbation amplitude is 20% Žrms.. We see that both the neutron yield and ion temperature decrease when the mode number 2 dominates, as observed in the single mode simulation. Also, the ion temperature does not change considerably when mode number 12 is mixed to mode number 6. Although the neutron yield is relatively sensitive to mixing of mode number 12, the value remains within the range bounded by calculated values for the single mode 6 and single mode number 12. Therefore, we conclude that the effect of multi-modes does not change the previous conclusions obtained from the single mode simulations. 4. Energy transport In Fig. 5, time development of the kinetic energy and the internal energy of the fuel and the neutron yield are shown. We compared the 1D simulation and the 2D simulation with mode number 6 and the velocity perturbation amplitude 10%. We see that the internal energy of the fuel increases with time, and around the maximum of the internal energy, the neutron yield rapidly increases and at about 2.4 ns, the neutron pro-
Fig. 4. DD neutron yield and Ti TOF Žion temperature of the fuel. for combination of mode 6}mode 2, mode 6}mode 12. The total velocity perturbation amplitude is 20% Žr.m.s...
duction terminates. The kinetic energy decreases from the beginning of the stagnation phase in both 1D and 2D simulations. Therefore, conversion of the kinetic energy to the internal energy is seen both in 1D and 2D simulations at the beginning. However the kinetic energy of the fuel increases later in 2D, though it remains small in 1D. This indicates that the hydrodynamic instability seen in 2D tends to reduce the internal energy in the later phase of implosion. In order to understand the reduction of the internal energy in 2D simulation comparing to the 1D case as shown in Fig. 5, we plot in Fig. 6 the time developments of the heating, electron thermal conduction loss and radiation loss in both 1D and 2D simulations. Since neutron production terminates at ; 2.4 ns as seen in Fig. 5, we focus our attention to the phenomena around 2.4 ns in Fig. 6. In 2D simulation, the net loss by electron thermal conduction and radiation is larger than that in 1D by about 0.5 J. However, heating in 2D is less than that in 1D by 4 J. This explains the difference in the internal energy in 1D and 2D simulations shown in Fig. 5. Therefore, we conclude that in the present case of the implosion, the reduction of the internal energy in 2D compared to 1D is mostly due to reduction of fuel heating by hydrodynamic compression, and not
168
A. Sunahara et al. r Fusion Engineering and Design 44 (1999) 163]169
Fig. 5. Time integrated DD neutron yields and time developments of the internal energy and kinetic energy of the fuel for the ID and 2D simulations. The mode number 6 with perturbation amplitude 10% is assumed in the 2D simulations.
due to the losses by electron thermal conduction and radiation.
5. Conclusion
We have calculated the hydrodynamic instability in the stagnation phase, and its effects on the energy transport and neutron reduction with use of 2D simulation code. Parameters the same as the present implosion experiment Žshot a19714. are used for the simulation. By comparing the simulation with the experiment, it is concluded that mode number 6 dominates with the perturbation amplitude of 20]30% required to explain the experimental results. It was demonstrated that the spikes generated by the growth of the hydrodynamic instability heats the center of the fuel. This explains the experimental finding of decreasing neutron yields when ion temperature increases, whose trend can be seen in other experimental results. Simulation including multi-modes showed that combination of modes 6 and 2 reduces the neutron yield and the ion temperature. On the other hand, combination of modes 6 and 12 increases the neutron yield, but the ion tem-
Fig. 6. Time integrated values of heating and cooling energy of the fuel for 1D and 2D simulations. PdV work by the pusher and shock heating contribute to the heating of the fuel. The cooling of the fuel is caused by the radiation energy loss and electron thermal conduction loss. The mode number is 6 and the perturbation amplitude 10% case.
perature is remaining constant. We conclude that in the energy transport, hydrodynamic heating by PdV work is reduced by the growth of the hydrodynamic instability and neutron yield is reduced. However, ion temperature increases for the case of mode 6 or 12 are dominant for the growth of the hydrodynamic instability. Acknowledgements We are grateful to Dr N. Izumi, Institute of Laser Engineering, Osaka University for discussions. References w1x H. Takabe, M. Yamanaka, K. Mima, et al., Phys. Fluids 31 Ž1988. 2884. w2x K. Mima, Y. Kato, H. Azechi, et al., Phys. Plasmas 3 Ž1996. 2077. w3x P. Colella, J. Comput. Phys. 54 Ž1984. 174. w4x B. van Leer, J. Comput. Phys. 32 Ž1979. 101. w5x P. Lotstedt, J. Comput. Phys 47 Ž1982. 211. w6x S.I. Braginskii, in: Leontovich ŽEd.., Review of Plasma Physics, vol. 1, Consultant Bureau, New York, 1965 p. 205.
A. Sunahara et al. r Fusion Engineering and Design 44 (1999) 163]169 w7x G.C. Pomraning, The Equation of Radiation Hydrodynamics, Pergamon Press, Oxford, 1973. w8x H. Takabe, T. Nishikawa, J. Quant. Spectrosc. Radiant. Transfer 51 Ž1994. 379. w9x H.-S. Bosch, G.M. Hale, Nucl. Fus. 32 Ž1992. 61. w10x N. Izumi, private communications. w11x N. Miyanaga, M. Nakatsuka, T. Kanabe, T. Jitsuno, K.
169
Tsubakimoto, N. Nishi, S. Matsuoka, A. Ando, S. Amano, S. Urushihara, N. Morio, T. Kawasaki, K. Suzuki, S. Matsuo, H. Asahara, K. Kawakami, M. Murakami, H. Azechi, K. Nishihara, S. Nakai, Inst. Phys. Proc. of Laser Interact. with Matter Conf., Conf. Ser. No 140, Section 4, IOP Publishing, Oxford, 1995, p. 81.