Quant Principle

ACKNOWLEDGMENTS

The acknowledgments have been removed from this online version. To see the Acknowledgements please contact q.boiler.

Contents

LIST OF TABLES

Table Page

LIST OF FIGURES

Figure Page

NOMENCLATURE

DT = deuteriun and tritium mixture
n_i = number density of species i
τ_E = energy confinement time
W_dt = energy out per deuterium tritium reaction
η_T = thermal conversion efficiency
η = amplitude of some interfacial perturbation
c_s = sound speed
ρ = mass density
U_DT = internal energy of DT fuel
E_F = fermi energy
AR = aspect ration
CR = convergence ratio
λ = wave length of a sinusoidal perturbation
γ = growth rate of linear Rayleigh-Taylor instability
γ = ratio of specific heats
C_p∕C_v (typically about 5/3)
ω = imaginary part of growth rate
D = d/dy
β = denstiy gradient factor
ϵ = specific internal energy of a fluid system
k = wave number
g = acceleration due to gravity
v_a = velocity of ablating material
α₁ = ρ₁∕(ρ₁ + ρ₂)
α₂ = ρ₂∕(ρ₁ + ρ₂)
A = Atwood number

ABSTRACT Alcock, Bryce L. M.S.N.E., Purdue University, May 2007. Eulerian Hydrodynamics modeling. Major Professor: Chan K. Choi. Alcock, Bryce Leander. Master of Science in Nuclear Engineering, Purdue University, August 1996. Eulerian Hydrodynamics Modeling the ICF Rayleigh-Taylor Instability. Major Professor: Prof. Chan K. Choi. The use of an Eulerian hydrodynamic code for the modeling of ICF type hydrodynamic instabilities is investigated. The Rayleigh-Taylor instability is the most significant limitation for the inertial confinement fusion (ICF) scheme. This instability is investigated by looking at the linear growth regime of superposed fluids. The RAGE (Radiation Adaptive Grid Eulerian) solution of this problem is compared to the small amplitude analytic solution. The computational solutions on the linear growth rates are found to compare favorably with the analytic results. The non-linear evolution of the Rayleigh-Taylor instability is then looked at to gain a qualitative knowledge of Eulerian hydrodynamics in a highly rotational flow. The formation of late-time bubbles and spikes are well represented by the Eulerian RAGE algorithm. The RAGE is used to model the hydrodynamic implosion process of an ICF target. The model is found to agree well with the experimental observations in that the time-dependent implosion of the inner radius of the target is matched remarkably well with the experimental measurements. A time-dependent temperature source with Planckian radiation diffusion was used to model the laser energy deposition. This was then compared to the results of an experiment and a one-dimensional Lagrangian hydrodynamics code with multigroup radiation diffusion. The calculational results from RAGE are shown to model correctly the laser driven implosion phase of an ICF target, by simulating the perturbation growth due to Rayleigh-Taylor instability at both linear and late-time non-linear regimes, and also by predicting the time-dependent implosion of the inner radius of the ICF target.

Chapter 1
INTRODUCTION

1.1

1.1 Historical Background

It is well known that as the world becomes more industrialized the need for larger and more vast sources of energy will become important. Limitations of the Earth’s resources of conventional fuels are becoming more apparent. It will be crucial to develop alternative energy sources before the current natural resources become exhausted. It has been known for several decades that nuclear fusion reactions are the major energy source in stars. Typical fusion fuel, e.g. deuterium, can be extracted from sea water in one part per 6700 of hydrogen atom. Because of this virtually limitless source of energy, scientists have turned their attention toward fusion for a new source of energy [1]. Nuclear fission has been known and used as a source of energy production for some time now. The draw backs to nuclear fission as an energy source are the limited fuel supply and radioactive waste disposal. On the other hand the concept of nuclear fusion does not have the problem of radioactive waste nor the limited fuel supply. This fact makes nuclear fusion an attractive clean source of energy for the future.

Fusion is the process of combining light atoms to make heavier ones. This process goes on in the sun every day. An example of a fusion reaction is the combining of a deuteron, D⁺ (²H) and a triton, T⁺ (³H) which is given as

(1.1)

In this reaction the deuteron and triton combine to form an alpha particle (α) and a neutron (n) and release 17.6 MeV of energy which is divided between the neutron and the alpha particle. Due to conservation of momentum the lions share (80%) of the released energy will go to the neutron. This is not the best scenario because the neutron is not a charged particle and hence will not follow magnetic field lines for direct energy conversion. Therefore the neutron cannot be directly converted to electrical energy. Because of this a thermodynamic heat cycle will need to be employed to convert the kinetic energy of the neutrons into electrical power which does not exceed the Carnot cycle efficiency. A more efficient way of converting the kinetic energies of fusion products is to use a direct energy conversion scheme in which the charged fusion products are collected directly via the magnetic field and the electric potential grids. The energy conversion efficiencies with this direct energy conversion scheme could be as high as 85-90%.

In order to fuse the ions one must meet three rather stringent requirements simultaneously. These requirements are temperature (T), number density (n_i), and confinement time ( τ_E). A large temperature is required in order to over come the repulsive Coulomb force seen by the similarly charged ions. A large product of number density and confinement time is needed to make the total number of fusion reactions significant. The product of the number density and confinement time is known as the Lawson criterion. The Lawson criterion is derived by balancing the energy out with the overall energy in, and is given as  [2]

(1.2)

for scientific break even. Where k is the Boltzmann’s constant, η_T is the reactor energy conversion efficiency, b is the bremsstrahlung coefficient, < σv > is the reaction rate parameter, which is dependent on the temperature of the ions, and W_DT is the fusion energy out of the reaction. For the D-T fusion reaction, this is 17.6 MeV.

There are different methods for confining fusion a plasma. Whichever method is used the Lawson criterion must be met. On the sun a gravitational confinement scheme is used. This would not be practical, however, in the laboratory due to the large mass of the reactor required to achieve self-gravitational attraction. Two other schemes are magnetic confinement fusion (MCF) and inertial confinement fusion (ICF). In the MCF scheme, charged particles are confined to follow magnetic field lines. The MCF schemes can be divided into two categories; open and closed magnetic confinement fusion systems. In the open magnetic confinement systems, field lines leave the plasma confinement area and allow the plasma to flow out of the confinement area. With the closed confinement system the magnetic field lines stay within the plasma confinement area. Most of the closed confinement systems take the form of a toroidal shape. In the ICF scheme, however lasers or particle beams are used to confine the plasma to super high densities. When the plasma is compressed to a large density significant fusion burn can occur prior to expansion of the pellet due to rarefaction waves. This scheme makes possible efficient thermonuclear burn and makes possible fusion power reactors with practical size lasers [3].

1.2

1.2 The ICF Scheme Overview

As noted earlier the ability to achieve fusion is governed by satisfying the Lawson criterion. In order to do this, one must have a high product of number density and confinement time. With the ICF concept, the idea is to have a very high number density, but short confinement time. The large number density will cause significant fusion burn such that in a short confinement time significant fusion burn can occur. The fusion reaction time for the ICF scheme is the burn time of the pellet. The disassembly time R∕c_s is about equal to the time it takes for a shock wave to travel through the center of the pellet. The burn time, however, is given by R∕4c_s where R is the radius of the compressed pellet, c_s is the sound speed, and the takes into account the fact that the pellet burn time is less than the disassembly time. Thus the nτ criterion can be expressed as an equivalent ρR criterion for the ICF scheme. One can easily see the relationship n = . Therefore the nτ relationship becomes

(1.3)

The ρR criterion is the figure of merit for the ICF scheme.

Figure 1.1 illustrates ICF process which is similar to the internal combustion engine in that it has the following steps  [1].

---

---

The first step is the injection, this is where the fuel or fuel pellet is injected into the combustion chamber. The second step is compression, in the ICF case this is where the pellet is compressed to many times the initial density. The lasers in the ICF scheme can be thought of as a spherical piston. This step is followed by the ignition due to high compression, here a micro core of the fuel is heated to ignition temperatures. The final stage is where the burn is propagated outward, followed by the pellet break up due to rarefaction waves. 1.2.1

1.2.1 ICF Hydrodynamics

In Figure 1.2 the stages of the implosion of an inertial confinement fusion target are shown.

---

---

In the first stage driver beams (lasers, electron beams, or ion beams) deposit energy on the surface of the capsule. Stage two sees the formation of a plasma atmosphere. Stage three has the driver beam energy being absorbed in the atmosphere, followed by stage four where the ablation driven imploding shocks are seen. These shocks compress the fuel core, which is followed by ignition. The last stage of the process is the burn propagation, which is terminated when the pellet breaks up due to the rarefaction waves.

The choice of driver is of important consequence in the field of ICF fusion. The main driver choices are the laser, relativistic electron beam, light ion beam and the heavy ion beam. If lasers are used it will most likely be the gas lasers, because of the need for high power and high efficiency. The CO₂ gas laser could have the high power and efficiency, but the problem is that the 10.6 μm wave length is too long for a good coupling with the target and preheat problems. Therefore the gas laser based on the krypton-fluorine mixture will be a better choice because of its 0.2 μm wave length. The relativistic electron beam has many problems that will make it a poor choice. The relativistic electrons produce hard x-rays from bremsstrahlung that penetrate the target and cause preheating of the target. It is also difficult to focus electron beams. The third difficulty with the relativistic electron beams is the long electron energy deposition range. These difficulties raise serious doubts about the ability of this driver. The light and heavy ion beams look good. The reasons that the ion beams look good is mainly two fold; first the ion beams can have a large beam energy density, and secondly the beams can deposit their energy in a very small region defined by the Bragg peak. The only real problems with the ion beams could be the cost, and the vacuum technology. With the ion beams the vacuum in the target chamber must be good to get good beam transport.

The National Ignition Facility (NIF) is a proposed laser facility that will demonstrate the ignition of the laser fusion target. The facility will be designed with 192 laser beams delivering 1.8 MJ of energy to an ICF target. The proposed NIF target design is shown in Figure 1.3  [4].

---

---

The implosion of a laser fusion pellet can be broken down into three distinct hydrodynamic phases. Two of the phases are clearly unstable to the Raleigh-Taylor (R-T) instability, which will be discussed in chapter 2. These phases are illustrated in Figure 1.4.

---

---

1. The shell is accelerated inward. During this acceleration the shell is Rayleigh-Taylor (R-T) unstable at the outer ablation front and at any interface where a lighter fluid is pushing on a heavier one.
2. The coasting phase occurs after the shell has reached its final velocity. In this phase one possibility is that the perturbation evolves such that the mass of material in the perturbed region is conserved.
3. Finally the shell decelerates as the gas converges on the central region. In this situation the inner surface is also Rayleigh-Taylor(R-T) unstable.

After the pellet has ignited a central core, a burn wave will propagate outward. The burn wave will cause significant heating to bring the rest of the pellet to thermonuclear burn temperatures. The pellet will continue to burn until it is blown apart. The time it takes for the pellet to disassemble, as mentioned earlier, is about equal to the time that it takes a sound wave to travel to the center of the pellet which is about 1 ns. The disassembly of the pellet marks the end of the cycle. The chamber is then cleared and the next pellet is injected. This can then be a continuous mode of operation. 1.3

1.3 Literature Review

The production of energy in stars is due almost entirely to the combination of four protons and two electrons to make an α particle. There are two ways this reaction can take place. The first is the direct combination of the hydrogen ions. The second is a method that uses carbon and nitrogen as catalysts. The reaction that uses nitrogen and carbon as catalysts is about 1000 times more likely than the direct combination. This fusion cycle of the sun has been known since the thirties [5]. One of the most challenging goals of this century has been to harness that energy of the sun in the laboratory. Hydrogen may be compressed to more than 10,000 times liquid density by an implosion system energized by a high energy laser. This scheme makes possible efficient thermonuclear burn of small pellets of heavy hydrogen isotopes, and makes feasible fusion power reactors using practical lasers [3].

The inertial confinement fusion scheme uses laser, ion beam or electron beam energy to heat up the outer shell. This causes an ablation pressure which drives the fuel or payload inward. Where the payload is the DT fuel plus some of the ablator material. In order to achieve large target gains two conditions will have to be met. First the majority of the mass will have to be compressed isentropically. This will minimizes the work needed to do the compression. Secondly a small central region will need to be brought to a large enough temperature to ignite. Once this small “micro core” has ignited, α particles from the fusion reactions will cause heating of the rest of the fuel [6]. By using Bernoulli’s equation we know PV = mv_f² where v _f² is the final velocity of the payload. Then mv_f² = U _DT ∝ ρ^2∕3 where U _DT is the internal energy of the DT calculated by using a Fermi-Dirac equation of state for kT ≪ E_F where E_F is the Fermi energy. It is clear that for a given final density, the driving pressure can be reduced by increasing the volume enclosed by the shell. If the volume is made too large the shell will become too thin to withstand the shell break up caused by the Rayleigh-Taylor instability [7]. If the perturbation converges a distance r₀(1 -) where r₀ is the initial radius, and CR is the convergence ratio, the most dangerous wavelength is λ_d = ϵd where ϵ is a constant of order unity, an estimate of the number of growth times given by Kothe et al. [8] is nt ≈ [4απ(AR)(CR - 1)∕ϵ]^1∕2 where AR is the initial fuel layer aspect ratio, defined as AR ≡. Evidently the most stable targets are those having thick shells not requiring high convergence for ignition [8]. This indicates that there is a trade off between high gain and stability, and for a successful demonstration of the inertial confinement fusion scheme this trade off must be optimized. As mentioned earlier there are three criteria that must be met; confinement time, number density, and temperature, simultaneously. In the ICF scheme the temperature is created by the entropy that is generated due to the converging shock waves. The plan is to have the shock waves heat a small central core called the hot spot, and then have the burn propagate outward causing self-heating significant enough to ignite the rest of the pellet. A substantial fraction of the energy absorbed by the ablator of an ICF target is expended to bring the central hot spot to ignition temperatures. If the hot spot convergence ratio is reduced by a factor of 2 the hot spot radius is 2 fold increased and 8 times more fuel must be heated to ignition temperatures [9]. Thus more energy is needed to ignite the pellet. For this reason only a small micro core is desired to be ignited. This will minimize the driver energy requirement.

The Rayleigh-Taylor (R-T) instability is the most significant limitation of inertial confinement fusion [10]. The acceleration stage and the deceleration stage are both unstable to the R-T instability. In the acceleration stage the low density ablator is accelerating the high density pusher. In the deceleration stage the high-pressure low-density fuel decelerates the shell, thus the inner surface is R-T unstable. A successful demonstration of the ICF will require understanding and minimizing the growth of the fluid instabilities that arise during the implosion of the fuel containing capsule [11].

1.3.1

1.3.1 Current Technology

Numerical simulation reproduces the correct analytic growth rates and eigenfunctions for classical linear problems, including feed through of perturbations from one interface to another. For laser-driven plastic targets at an intensity of 2 * 10¹⁴ growth rates were found to be 40-80 percent of the classical rates which are consistent with more recent calculations [12]. Unfortunately, a fully Lagrangian or semi-Eulerian code bowties when strong shear develops in the flow (A bowtie is a quadrilateral zone in which two opposite sides intersect each other) [12]. To over come this problem computer codes have been developed which use Eulerian coordinates. The problem with modeling the inertial ICF target on the computer is that there are no experiments to benchmark the calculations with because there is no direct way of measuring the growth of the instabilities. For this reason the problems of convergence effects have also been formulated in cylindrical geometry. Cylindrical geometry has the distinct advantage of creating an implosion which will have the deceleration phase while still allowing for a direct diagnostic line-of-sight [11]. Investigations of the R-T instability on the pusher-fuel contact surface in the deceleration phase of the implosion in cylindrical geometry using a two-dimensional fluid code IMPACT-2D [13] have been done. IMPACT-2D (IMPlosion Analysis Code with Tvd scheme) is fully Eulerian, which is necessitated from the requirement to model rotational flows accurately.

1.4

1.4 Goals

The goals of this thesis will be to demonstrate the use of an Eulerian hydrodynamics code to model the Rayleigh-Taylor (R-T) instability. This will be done by looking at the R-T instability for various Atwood numbers and various wave numbers, k. Then the implosion of a cylindrical capsule will be done. This will be compared to experimental data and to that of a Lagrangian hydrodynamics code.

The organization of the thesis is as follows: Chapter two will contain theoretical analysis of the Rayleigh-Taylor instability as pertaining to the inertial confinement fusion scheme. Chapter three will discuss the numerical methods that are used for the computational modeling of ICF, and will have a detailed discussion of the Godunov method. The Godunov method is the method employeed by the Eulerian hydrodynamics code used in this study. Chapter four will contain the computational results of this study and will have comparison to experimental data and a Lagrangian hydrodynamics code.

Chapter 2
analysis of hydrodynamic instabilities

Hydrodynamic instabilities play a major role in the current approaches to ICF. The direct drive approach has many lasers or particle beams focus on a target. This heats up an ablating layer which then compresses the fuel to high densities. In the indirect drive approach photons from high intensity lasers are absorbed by means of a high Z material, which then emmits x-rays. These x-rays drive the implosion. In both cases the low density plasma which is blown off accelerates a high density pusher at the ablation front. In the accelerating frame of the pusher a heavy pusher is being accelerated by the low density ablating plasma. This sets up the Rayleigh-Taylor (R-T) instability. Small imperfections on the capsule will grow and could cause the capsule to break up or seed perturbations onto the inner surface of the shell. The perturbations seeded to the inner surface would grow during the deceleration phase. This would degrade the capsule performance. The understanding of hydrodynamic instabilities is of much importance to the success of the inertial confinement fusion scheme. Because of this much work has been directed toward establishing the growth rate of the ablative R-T instability. In this chapter a development of the classical R-T instability will be presented. The relationship of the classical R-T instability to the ablatively accelerated R-T instability will be discussed. A review of the current research in the area of the ablative R-T instability will be presented. Finally a discussion of two other hydrodynamic instabilities, the Kelvin-Helmholtz (K-H) and Richtmyer-Meshkov (R-M), will be reviewed. Both the K-H and R-M instabilities will play a role in ICF. 2.1

2.1 The Rayleigh-Taylor Instability

The R-T instability is seen when a less dense fluid is accelerating a more dense fluid. This situation is of particular interest in the ICF area. The reason for this is that the R-T instability must be well understood and minimized in order to maximize the implosion efficiency, thus maximizing the gain of an ICF pellet. The ICF pellet will pass through three regimes as was discussed in Figure 1.4. The first regime is the acceleration of the dense pusher by the ablating plasma, this is R-T unstable due to the fact that the light ablating plasma is accelerating the more dense pusher. The second regime is the coasting regime where the instability may grow, not due to the acceleration of the fluid but instead due to the conservation of perturbation and geometry effects. The third regime is the deceleration phase when the low density D-T is decelerating the pusher.

The planar R-T instability driven under constant acceleration has three phases that are well understood at this point, and two phases that are not as well understood. A diagram of the first three phases is shown in Figure 2.1. In the first phase an unstable

---

---

situation is seen where a more dense fluid is sitting on top of a less dense fluid, or a more dense fluid is being accelerated by a less dense fluid. In this situation any small perturbation will tend to grow. It is interesting to note that if the interface between the more dense and less dense fluid is perfectly smooth the situation is stable until a perturbation is introduced. The second phase of the R-T instability is the linear growth regime. In this regime the perturbation grows such that η = η₀e^γt, where γ is a function of , k, Atwood number, and geometry, with being the acceleration and k the wave number defined such that λk = 2π. Once the amplitude of the perturbation has grown such that η ≈ the growth rate decreases and is governed by the bubble velocity and spike velocity as seen in (c) of Figure 2.1. Once this stage is reached the perturbation is said to have reached saturation. Following the bubble and spike regime, the Kelvin-Helmholtz (K-H) instability is seen in the region between the bubbles and spikes. This instability arises mainly due to the velocity shear and will be discussed in more detail later. Finally the two fluids will enter a regime of turbulent mixing.

The derivation of the classical growth rate for the Rayleigh-Taylor instability has been done by many authors; most notably, see  [14] [15]. One can derive the governing equation of the R-T instability by looking at the equations of motion, mass conservation, and incompressibility. If one ignores viscosity, surface tension, and other non-ideal effects the equations of motion, mass conservation, and incompressibility become

These equations can then be linearized such that ρ = ρ₀ + ρ₁,p = p₀ + p₁ and v = v₁. Where the initial fluid is considered to be at rest. Then assume that the perturbed quantities vary as f(x,z,t) = f(x)e^i(kz-wt), thus one can linearize the equations to get the following equations These four equations can be solved by eliminating ρ₁,p₁, and v_z. Doing this will lead to the equation that governs the Rayleigh-Taylor instability which is given as

(2.8)

This is the equation for the simple case of no viscosity, no surface tension, and incompressibility. There have been many authors that have looked at compressibility, surface tension, magnetic field, and other non-ideal effects see  [16] [17]. 2.1.1

2.1.1 Discontinuous Density Profile

The simplest case of the R-T instability is that of two constant density fluids separated by a discontinuity at the interface between them. This is the case that is shown in Figure 2.1 and will be solved analytically for a growth rate of the linear regime. Assuming that the density is constant in both the upper and lower fluids, the last term in equation 2.8 becomes zero. Therefore the equation that governs the motion is simplified to

(2.9)

where D=d/dy. This equation has two decaying solutions in the upper and lower half plane given as

(2.10)

Integrating equation 2.9 through the interface and using equation 2.10 one obtains

(2.11)

Thus if ρ₂ > ρ₁ we see that γ² > 0 and the system is unstable.

2.1.2

2.1.2 Density Gradient Profile

The ideal R-T instability as shown in Figure 2.1 has a growth rate not adequate for high gain ICF targets. One possibility would be to deposit the energy at various depths in the target. This would spread out the region of energy deposition and cause a density gradient. The density gradient growth rate will be derived for an exponential density gradient. This model takes into account the density gradient but doesn’t account for the time-dependent acceleration, compressibility, mass flow, and other non-ideal effects. The inertial confinement fusion scheme has the low-density high-pressure ablator pushing on the high-density fuel. The Atwood number for this situation is close to one  [10]. The gradient of density in the ablator can have a stabilizing effect on the growth rate. This can be seen by looking at the case of a density profile given by ρ = ρ₀e^βz. For this case the equation governing the growth rate is given as

(2.12)

This equation can then be reduced by defining D=d/dz and dividing through by ρ₀e^βz. The result is a second order ordinary differential equation given as

(2.13)

This equation then has solutions of the form

(2.14)

The solving for the solution in a box that extdends from z=0 to z=d, gives the fluid boundary conditions such that v=0 at z=0 and z=d one then gets the following two expressions

where in the second equation m is any integer. The general solution of the expression is then obtained to be

(2.17)

From this equation it can be seen that the stratification is stable if β is negative and unstable if β is positive [16].

The exponential density profile does not completely describe the situation for ICF. The ICF case will be approximated by an exponential density profile in the region where the energy is being deposited, joined by two regions of constant density outside this exponential region. The region in which the energy is being deposited will extend from x=0 to x=d. The mathematical expression for such a density profile is given as

(2.18)

Solving for the velocity in the three regions will give a solution of the following form

(2.19)

Where A,B,C, and D are constants and γ₁ and γ₂ are solutions of equation 2.13 for the region 0 < x ≤ d which are

(2.20)

Multiplying equation 2.8 by W^m and solving for one obtains

(2.21)

Then taking m=0 gives

(2.22)

Substitution of equations 2.18 and 2.19 into 2.22 gives

(2.23)

where A, B, C and D are constants from the eigenvalue solutions of equation 2.8, β indicates how fast the density profile increases, ρ₁ and ρ₂ are the densities in the regions outside the exponential density profile, k is the wave number and γ₁ and γ₂ are given by equation 2.20. This growth rate will be more important in the light and heavy ion driven fusion due to the fact that with these approaches the energy deposition region is more exactly controlled by the Bragg peak.

2.1.3

2.1.3 Effect of a Finite Geometry

The fact that the fluid interface is located within a finite geometry will have an effect on the growth rate of the instability. This effect can be seen as a slight decrease in the growth rate. For the case of two stratified fluids in a box of height extending from -b to a, the growth rate for two constant density superposed fluids becomes [18]

(2.24)

Looking at equation 2.24, it is clear that for ka = kb ≥ 2; 1.037 ≥ coth(ka) ≥ 1.0 thus coth(ka) ≈1 for this situation. Thus for large ka and kb the geometry is not an important effect and the system can be viewed as having the same growth rate as the infinite system while in the linear regime.

2.1.4

2.1.4 Nonlinear Growth of the R-T Instability

When the amplitude, η, of the R-T instability becomes large compared with the wave length, the growth rate is seen to decrease, until it reaches the asymptotic limit of the bubble rise velocity and spike freefall. The overall flow is then governed by the rise velocity of the bubble. Recall Figure 2.1, in this Figure the velocity of the bubble and the spike are not the same, but if one assumes incompressible flow then the volumetric flow rates must be equal. Therefore the velocity of the spike will be larger than the velocity of the bubble due to the fact that the bubble has a larger area.

In the case of the ICF pellet the growth of spikes tend to saturate. Their growth decreases significantly at large amplitude and in some cases the spike structure becomes approximately stationary except for the time-dependent mass removal effects [19]. For ICF ignition one would like to have a large aspect ratio, where for ICF aspect ratio is defined as AR ≡. The reason for this is that in general the performance of the ICF pellet increases with increasing aspect ratio. The problem with the large aspect ratio is that the pellet becomes more susceptible to break up caused by the R-T instability. Some of the typical designs require the aspect ratio AR to be in the range 25-50 [10], and due to fabrication limitations it is thought that the outer surface finishes of the pellet will be no better that 250 angstroms. Therefore the aspect ratio is limited to

(2.25)

where ϵ is a growth factor reduction from the classical R-T growth rate. And n is a factor that indicates the number of e foldings of the perturbation amplitude. For this situation n was taken to be six. It is then seen that ϵ must be about 0.4-0.5 for the ICF to be successful.

2.1.5

2.1.5 Ablatively Driven Rayleigh-Taylor Instability

The ICF community has proven that the R-T instability has a significantly reduced growth rate due to ablative stabilization. A widely accepted dispersion relation for the ablative R-T instability is give by [20] as

(2.26)

In this equation k is the wave number, g is gravity, β is a factor in the range of 3-4, and v_a is the velocity of the ablating material. There have been many papers which confirm the finding of a significant reduction in the R-T growth rate. Some of these papers include theoretical  [21]  [22] [10],computational  [23] [12]  [24] [25], and experimental works [26] [27]  [28] [29] [30] [31]. Most of this work has been done in the planar geometry because of ease of diagnostics for the experimental work. The planar geometry allows a direct face on image of the experimental results which is not possible in the spherical geometry. The planar geometry also allows for a simpler dispersion relation.

2.2

2.2 Kelvin-Helmholtz Instability

The K-H instability will arise when two stratified fluids are in relative motion. An example of the K-H instability is the waves on the ocean. In this example the wind causes the waves to grow. The waves are stabilized by gravity and surface tension. The main aspect of the K-H instability is that in the absence of gravity and surface tension any perturbation that is observed would grow if there is relative velocity between the two fluids. This is independent of the relative density of the two fluids. The equation for the growth if the K-H instability is given as [16]

(2.27)

where α₁ = ρ₁∕(ρ₁ + ρ₂). This instability in the absence of surface tension and gravity reduces to an interesting equation which could be applied to inertial confinement fusion applications. The equation reduces to

(2.28)

This is clearly seen to be most unstable for α₁ = α₂. The reason this may be the most important growth rate for the ICF application is that in ICF the velocity shear is occurring at the interface between the bubble and spike, thus there may be no accelerations in the plane perpendicular to the interface, and since the fluids are likely to be ionized the surface tension and viscosity may be negligible. Then equation 2.28 would clearly indicate that the most K-H unstable situations will arise when the two densities are about equal.

2.3

2.3 Richtmyer-Meshkov Instability

The Richtmyer-Meshkov (R-M) instability is similar to the R-T instability in that it is seen to occur in situations where the density is stratified. The R-M instability is different than the R-T in that it is driven by an impulsive force, such as a shock, as opposed to a constant body force such as gravity [32]. The growth rate for the R-M instability is also different than the growth rate for the R-T instability. The growth rate for the R-M instability is given as [18]

(2.29)

Where n is the initial amplitude of the perturbation and γ = ΔvkA, where v is the velocity jump across the shock, k is the wave number, and A is the atwood number. From this growth rate equation it is seen that amplitudes of perturbations can increase or decrease depending on the sign of v. That means that the amplitude will not know whether it is going to grow or shrink until it finds out the direction of the shock wave. But if the amplitude does shrink it will not stop when it gets to zero, it will instead shift phase and continue to grow. It can be seen from equation 2.29 that if the growth rate is negative it will shrink to 0 and then continue to grow at the rate γ. The concept of freezing out a perturbation could have a great deal of positive applications for the ICF field. This concept is when one applies the R-M instability to cause a mode perturbation to freeze, or stop growing. It can be seen that if the correct shock hits a perturbed surface that is growing the perturbed surface could be frozen, thus it would quit growing. A more detailed explanation of the R-M instability as related to the ICF program is discussed in references [33] [34].

Chapter 3
Eulerian Hydrodynamics and the Godunov Method

3.1

3.1 Eulerian vs. Lagrangian Point of View

A problem must be formulated in a coordinate system in order to be solved. There are two main coordinate systems to choose from, the Eulerian and the Lagrangian systems. Both of these reference frames have some advantages and disadvantages. The comparison of these methods will be made in this section.

The Eulerian view point of hydrodynamics is based on looking at the fluid that passes through a control volume which is fixed in space. Therefore in this view point the independent variables are x, y, z, and t [35]. When solving problems in this frame of reference one must formulate the laws of conservation such that the net flow of mass, momentum and energy that flow through the control volume are equal to the time rate change of these quantities in the control volume. The equations of mass, momentum and energy in this formulation then become

where ϵ is the specific energy of a mass of the fluid.

The Lagrangian frame of reference for hydrodynamics is different in that one is now interested in looking at a fixed mass of fluid as it flows. This is a control mass approach to the problem. Therefore one is to look at a particular fluid particle as it flows. Then the principles of mass, momentum and energy conservation are applied to the particular fluid particle. In this reference frame x, y, z, and t are no longer independent variables. Now, if x₀, y₀, and z₀ are chosen then at any time t later the values of x, y, and z are uniquely determined for that fluid particle [35]. Therefore the fluid particle that is chosen at t₀ has the independent variables x₀, y₀, and z₀. The equations of mass, momentum and energy in Lagrangian coordinates for the case of one dimension then become

It can then be seen that in the Lagrangian frame work the total time derivative of some quantity, m, is simply = . This is not the case of the Eulerian formulation where the mass that is entering and leaving the control volume must be taken into account. It can also be seen that the Lagrangian formulation will give more information than the Eulerian formulation, that is in the Lagrangian formulation one knows where the fluid particle was at at time t=0 and every time thereafter, but in the Eulerian formulation this information is not available [32]. This fact makes it hard to track the material interfaces in the Eulerian formulation and thus makes the Eulerian formulation a more difficult formulation to use.

The Lagrangian formulation is the more popular of the two for computational work because of the quicker run time and less numerical diffusion. The Eulerian formulation does have some advantages, however, and theses advantages are starting to become very important in the area of computational fluid dynamics (cfd). The ICF pellet will experience the R-T instability which will grow into the non-linear regime where the bubble-spike formation is seen. This bubble-spike formation distorts the contact surface and large vorticity is seen. This cannot be modeled well at late times with the Lagrangian algorithms because the mesh will get tangled. Thus if one is going to use the Lagrangian mesh a rezoning will need to be done, but this will introduce an unacceptable amount of numerical diffusion [13].

The Eulerian hydrodynamic code, RAGE, which stands for Radiation Adaptive Grid Eulerian, is being developed to solve this problem of large rotational flow. Another code that is also being developed specifically for this problems is a 3-D code called PLATO. The PLATO is a 3-D spherical hydrodynamics code with a fixed Eulerian grid. With PLATO the fluid equations are closed with the ideal gas law and γ = . PLATO does not have thermal conduction, laser energy deposition, or α particle heating. RAGE on the other hand has non-equilibrium gray radiation diffusion and can have a time-dependent radiation temperature source or time-dependent energy source deposited in a material or region of the user’s choice.

3.2

3.2 RAGE and the Godunov Method

The computer code RAGE (Radiation Adaptive Grid Eulerian) is a multidimensional, multi-material, Eulerian hydrodynamics code with adaptive mesh refinement (AMR), second-order accurate numerical techniques, and non-equilibrium gray radiation diffusion with flux limiters [36] The code was developed by SAIC (Scientific Applications International Corporation). RAGE uses the higher order piecewise linear Godunov numerical method.

3.2.1

3.2.1 Survey of Numerical Methods for Discontinuous Flow

There are three major approaches to the representation of discontinuities. These three major approaches will be discussed from a theoretical point of view and the advantages and disadvantages of each will be pointed out. The first approach is that of introducing artificial viscosity. In this method a non-physical pressure term is introduced to spread the discontinuity over several mesh. This can be done by introducing an artificial viscous pressure such that this term is large near steep gradients of flow and zero everywhere else. This term will not affect the flow pattern, but will tend to spread the discontinuity over several cells. The most powerful advantage of this method is that it is very easy to implement and it is the most common method of handling flow discontinuities. The disadvantages of this method include large computational time and large computer memory requirements. The reason for the demand on computer time and memory is due to the fact that discontinuities are spread over several mesh, thus if one wants to get fine representation of the discontinuity very fine mesh must be used, and as higher resolution is desired the mesh must get smaller. The ultimate sacrifice is a prohibitively small Courant time step [37].

The second method is that of linear hybridization. This is a hybrid of a low order scheme and a high order scheme. The high order scheme is used in areas of smooth flow, and is very good at representing the flow in the smooth regions, but the high order scheme is not stable near discontinuous flow and will give non physical oscillating solutions. Thus a linear combination of the high and low order schemes are used in the areas of discontinuous flow. The low order scheme acts like artificial viscosity in that it introduces terms that are non-physical which help to stablize the solution in these regions. The disadvantages of this method include the difficulty in programming and the difficulty in determining a good weighting factor for the post shock oscillations. The advantage of this method is that the discontinuity can be represented with very good resolution [37].

The third approach, the Godunov method makes use of the knowledge of the propagation of flow discontinuities. This knowledge is built in to the method in the form of the Reimann solver. The Reimann solver determines the interaction of the two states and then tells what form of wave is going to emerge. From the interaction of a discontinuity three waves will emerge. Two of these will either be of the rarefaction or compression variety and the third will be a contact wave. The interaction of these waves is determined by the Reimann constants. This information is then used to calculate the other properties of the fluid. The major advantage of this scheme is the narrow structures with which flow discontinuities can be represented, and the fact that non-physical dissipative terms are not added to the equations. The disadvantage of this method includes the complexity of the computer algorithm for the Riemann solver. The Riemann solver becomes much more complex if the flow conditions cannot be represented by the ideal gas laws or the isothermal assumption. The equations are implicit in the case of strong shocks. This causes the need for an iterative solution of these equations. This method can take twice as much computational time per zone per time step as some artificial viscosity algorithms.

3.2.2

3.2.2 The Method of S.K. Godunov

In this section an Eulerian formulation of the Godunov method will be discussed. This method can also be formulated in the Lagrangian coordinates. The Eulerian equations in conservation-law form are chosen. The equations are

These equations must then be solved in conjunction with the integral relationships valid on any contour of integration [38]. These relationships if integrated over a discontinuity will give the same results as the Rankine-Hugoniot relations. These integral relations are The procedure for solving the problems is to note that if x_m is the interface between two grid points, then processes in not too long time τ will allow the values of u,ρ, p, and ϵ to be regarded as unchanged [38]. Thus one can define u(x_m,t) = U_m, ρ(x_m,t) = R_m, p(x_m,t) = P_m and ϵ(x_m,t) = E_m as constant values for the time of interest. From this a closed form of the difference equations for one dimension follows [38] as In these equations the superscript represents the next time step and the subscript represents the current time step. In order to close these equations one must add the equation of state which is of the form p = p(ρ,ϵ). When this method is extended to two dimensions, the calculation of the decay of discontinuity becomes more complicated. The form for the major constants are not in a closed form and must be solved using an iterative procedure.

3.3

3.3 The RAGE Algorithm and Philosophy

3.3.1

3.3.1 Philosophy and History of RAGE Development

The main motivation for the development of the RAGE computer algorithm was for the study of strong shocks in water. There is interest in the propagation and structure of such shocks in shallow water, and to analyze the effects of these on the ocean bottom [36].

Because of the need to resolve the the structure of the shock very well, yet look at very late times, an adaptive mesh was needed. The need for high resolution of flow discontinuities also required the need for a second-order accurate numerical method. Thus the second-order accurate piecewise linear Godunov numerical method in Eulerian coordinates was chosen for this problem. Along with the second-order accurate numerical method, a unique adaptive mesh algorithm was implemented. “RAGE is as much a philosophy as it is a hydrocode” [36].

3.3.2

3.3.2 RAGE Problem Approach

RAGE has a powerful problem generator built into the code. The first step is to chose the computational grid. Because of the adaptive mesh, this can be defined with only a few large mesh. This area must be chosen such that the entire area of interest is included. Therefore, if late times are going to be looked at, a very large grid must be chosen, such that any cell of activity is included in the computational grid. Once the initial grid is chosen and the material regions are chosen RAGE configures the grid. This is done by looking at each individual cell, if a cell contains different material or the same material at different temperature, pressure or velocity that cell is divided into four cells. This grid refinement process is repeated until the desired level of resolution is achieved. The cells are then filled with the material defined at the centers of the cell. Initially all level one cells contain a pure material but the cell refinement allows material interfaces to be resolved to any level of desired resolution. If one is clever and chooses the material boundaries and grid size such that material boundaries and cell boundaries are falling on the same spot at relatively large cell sizes, it is not needed to subdivide the cells to the limiting number of cells to define the material interface for the initial cycle. Once the grid is initialized the hydrodynamic equations based on the Godunov method are solved. The time step is first calculated as dt_next = .9dx∕c_s where c_s is the sound speed, and the .9 is chosen because dt_next < dx∕c_s. The hydrodynamic phases are then calculated by splitting the equations into radial(x) and axial(y) sweeps. The sweeps are performed in alternating order. On even cycles the equations are swept in the radial direction, then in the axial direction. On odd cycles the sweeps are performed in the axial then radial directions. The sweeps are simply one-dimensional Lagragian hydrodynamic calculations which determine the impulsive work done by stresses and velocities. The solutions to these Lagrangian sweeps are then mapped back onto the original Eulerian grid. Following the two sweeps the radiation subroutine is run. This subroutine solves the non-equilibrium gray radiation diffusion equation. Once the entire hydrodynamic and radiation subroutines have been run, the entire grid is re-examined to determine if a mesh reconfiguration is needed. This is done in the subroutine RECON. There are several configuration criterion to determine if a cell is to be split or dezoned (unsplit). One requirement is that every cell can be no more than one level of refinement from its neighbor. No cell will be removed if the following cycle will force the region to be redivided. Also if a cell is divided due to an adverse gradient of pressure, temperature or velocity, then all of the neighbor cells will be divided as well.

3.4

3.4 Shock Tube Problem for RAGE Benchmark

The Godunov method and the adaptive mesh refinement features of RAGE will be demonstrated in this section. The analytic solution of a shock tube problem will be compared to the numerical solution from RAGE. The results will be compared both qualitatively and quantitatively.

There are many shock-tube type problems that can be solved analytically to compare to numerical calculations. The shock tube chosen for this comparison is quite simple. The setup is shown in Figure 3.1. There are two chambers separated by a diaphragm. The two chambers both contain an ideal gas with γ = 1.4. The initial temperature is 0.025 eV, the pressure to the left of the diaphragm is initially 400 × 10⁶ , and the pressure to the right of the diaphragm is initially 200 × 10⁶. At time t=0 the diaphragm is broken and three waves develop: one is a shock wave, another is a rarefaction wave, and the last one is a contact wave. These waves will propagate until they interact with a boundary.

The equations of the shock tube problem are easy to derive, but will be presented here with no derivation. The equations following from Amsden and Harlow  [39] are

In theses equations ρ₊ is the density in the region prior to the arrival of a shock, ρ_R is the density in the region just to the right of the contact wave in the shocked region, ρ_L is the density just to the left of the contact wave, and ρ_- is the density in the region to the left of the rarefaction wave. And P_- is the pressure to the left of the rarefaction wave. Using the same labeling for the regions, I is the specific internal energy. The velocity of the shock wave is v_s, the velocity of the contact wave is u_c, and P ≡P_c/P_-. The value of P is determined from the following implicit equation

(3.22)

Then P_- is assumed given and thus P_c can be determined once P is known. Table 3.1 contains the numerical results from the RAGE calculation which are compared to the results from the analytic solution.

The plot in Figure 3.1 shows the results of the pressure as a function x for time t=400E-6 sec. The two lines compare the effect of the adaptive mesh on the solution of this problem. The dotted line is with 1-cm mesh and no refinement allowed, and the solid line allows the mesh to refine to 0.125 cm. It is clear that the smaller mesh gives the better result.

---

---

---

---

The 1-cm fixed mesh run took 1.838 second on a Cray Y-MP C90, whereas the 0.125 cm-adaptive mesh took 10.525 seconds on the same computer. This indicates that the small mesh slows the time step down considerably, which is due to the fact that the time step for each cycle is controlled by the slowest time step on the grid. The slowest time step on the grid is often controlled by the smallest cell. The plot labeled Figure 3.2 shows the levels of mesh refinement and the density for the case in which the mesh was allowed to refine. It is seen that the areas of largest mesh refinement are those of the steepest gradients in density. The steep gradients in density, as well as steep gradients in other flow quantities, cause the mesh to refine. The most well refined areas are the shock front and the contact discontinuity. The rarefaction wave does not see as much refinement because the gradient of flow properties is not as steep.

Chapter 4
Modeling of ICF Phenomena

Modeling of the inertial confinement fusion pellet requires that the radiation driven implosion process and the hydrodynamic instability growth processes be accurately modeled. In the following sections these processes are looked at. The first section deals with the modeling of the Rayleigh-Taylor instability in the planar geometry. The second section deals with the modeling of the radiation driven implosion process. 4.1

4.1 Modeling of the Rayleigh-Taylor Instability

This section will be strongly coupled to the results of chapter 2. The results of chapter 2 are the theoretical predictions for the behavior of fluid instabilities. This section contains some of the computational results for parameters of the Rayleigh-Taylor instability.

There are three main areas of study for the comparison of RAGE to the theoretical analysis. These are growth rate (γ) vs. wave number (k) studies, growth rate (γ) vs. Atwood number (A) studies, and the nonlinear studies. Each of these three studies will be discussed in the following sections.

4.1.1

4.1.1 Growth Rate vs. Wave Number Studies

The study of the growth rate versus the wave number for the linear regime of the instability was done by comparing the theoretical value of γ to the computational value, which was calculated from the RAGE output.

In order to make this comparison a series of inputs of similar input decks with various initial perturbation wave numbers were run. The data were then analyzed and compared to theory. The computational problem setup is a 1.0 cm square region of interest. The grid had 256 level one mesh in both the axial and radial directions, and was allowed to refine such that the maximum number of mesh in either the axial or radial directions would be 4×256 = 1024 mesh. The amplitude, η, of the initial perturbation was 0.01 cm. The most refined mesh had a size of cm, thus there were about 10 mesh per perturbation at the initial time. In order simulate incompressible flow, RAGE has a feature called ‘rttype′. If this is set to one RAGE will force the density in two adjacent regions to stay constant while letting other variables change. This feature was used so that at initial times no compression (stressing) of the grid was seen. Also the specific internal energy and density were the two state properties defined. They were defined such that the pressure across the interface was zero.

The results of these runs are shown in Figure 4.1.

---

---

It is clear from these results that the values of large k, i.e., short wavelength, do not agree with theory. The results for small k, on the other hand, agree quite well with theory. The reasons for this disagreement are due to several factors. First, the the large values of k have less mesh in the radial (x) direction. Second, the larger values of k will reach the nonlinear regime faster. In order to model the linear regime of the R-T instability well for large values of k, one must use very fine mesh and very small initial perturbation amplitudes. The cost of making the mesh so small is a prohibitively slow time step. 4.1.2

4.1.2 Growth Rate vs. Atwood Number Studies

Studies of γ as a function of A, Atwood number, were done in order to get a qualitative feel for the ability of RAGE to model the R-T instability for various Atwood numbers. The studies were done by fixing gravity, wave number, and the density of the lower fluid and causing the density of the upper fluid to change for various runs. Several runs were made and it was found that for Atwood numbers between 0.2 and 0.5 there was good agreement with the theory. The resluts of this study are shown in Figure 4.2.

---

---

One problem that RAGE has in dealing with the Rayleigh-Taylor problem is that it is an Eulerian grid. This causes problems when material flows through the grid. The main problem is that there are cells that are composed of more than one material. These mixed cells then act as a density gradient. From the theoretical analysis of the R-T instability, it was seen that the density gradient acted as a stabilizing agent to decrease the growth rate of the instability. The growth rates of the instabilities in the γ vs. k section were lower than theory because of this effect, and because higher order modes were saturating and thus causing a decrease in the growth rate of the mode of interest. The higher order modes (see Appendix C) are seeded by the finite difference approximation of a cosine curve. These seeds would then grow in time and because of the large wave numbers these small modes would become nonlinear very fast and enter a nonlinear regime. The mode coupling of these high modes and the low modes changes the R-T growth and a nonlinear pattern evolves. At the initial times the modes grow independent of each other but as the higher wave number modes saturate they couple with the lower wave numbers and cause the funny flow patterns seen in Appendix C. 4.1.3

4.1.3 Nonlinear Evolution

The nonlinear evolution of the R-T instability is marked by the formation of rising bubbles of light fluid and falling spikes of heavy fluid. This motion is most easily characterized by the bubble rise velocity and the spike free-fall velocity. The thicknesses of the bubble region and spike region are determined by the Atwood number. The features of the nonlinear evolution of the bubble and spike regime are shown for two different Atwood numbers in Figures 4.3 and 4.4. The Atwood numbers for this flow are 0.5 and 0.2, respectively.

---

---

---

---

4.2

4.2 Modeling of Radiation Source

This section describes the experiment that is modeled along with results. In addition to this, the physical processes of laser interaction will be discussed and the method with which RAGE models these phenomena will be indicated.

4.2.1

4.2.1 Radiation Interaction with Matter and Plasma

The indirect-drive process of ICF sees a laser entering a hohlraum cavity and interacting with its wall. The wall material emits radiation, which then irradiates the ICF pellet thus driving an implosion indirectly. The emission of radiation from the wall of the hohlraum is seen to follow the behavior of radiation transport, whose equation is

(4.1)

Where I_ν is the beam intensity of frequency ν, c is the speed of light, the first term in the brackets on the right is the spontaneous emission, the second term in the brackets on the right is the stimulated emission, and the last term on the right is the absorption [40].

The radiation from the wall of the hohlraum will interact with the ICF pellet. The interaction process involves many mechanisms. These mechanisms are any of the photon interactions with a plasma, which include collisional absorption of laser light, resonance absorption, acoustic turbulence, and parametric processes. Each of these processes is responsible for the deposition of some radiation energy in the target. This energy is seen as an increase in the temperature of the electrons and ions or an increase in the plasma wave energy. The electrons are mainly responsible for the thermal conduction transport of energy. But the large temperature gradients are responsible for the radiation energy transport.

4.2.2

4.2.2 RAGE Model of Radiation Transport

The RAGE algorithm is a simplified model of the physical processes that are occurring in the ICF pellet. The laser interaction with the wall of the hohlraum, and subsequent emission of radiation from the wall is replaced by a time-dependent temperature profile in low density air surrounding the ICF pellet. The energy from this time dependent temperature will then heat the surface of the pellet through thermal conduction and non-equilibrium gray radiation diffusion. The thermal conduction follows Fouriers law of heat transfer with flux limiters implemented to exclude nonphysical transfer rates. The non-equilibrium gray radiation diffusion is approximated through the use of Rosseland opacities and is given as follows. If approximate local equilibrium is assumed, the total energy flux can be expressed as

(4.2)

Where σ is the Stefan-Boltzmann constant and l is the average value of the radiation mean free path. The value for the average mean free path, l, is

(4.3)

where the factor G(u) is the Rosseland mean free path and is given as

(4.4)

4.3

4.3 Description of the Experiment Modeled with RAGE

A group of scientists at Los Alamos National Laboratory  [41] are working with cylindrical implosions in order to analyze the growth rate of the R-T instability in this geometry. The reason for using cylinders is that there is a direct line of sight for ease of diagnostics while maintaining a convergent geometry. The geometry of the experimental setup is shown in Figure 4.5.

---

---

The cylinder is sitting in the hohlraum transversely for multiple reasons. First, if the cylinder was sitting coaxially the x-rays from the gold would not hit the cylinder as evenly and this would lead to a nonsymmetric implosion. Second, the coaxial arrangement would allow for x-rays to enter the axial(z) direction of the cylinder thus causing preheat and other non-ideal effect to the inner foam. Finally, the transverse arrangement allows for an easy diagnostic. The diagnostic is setup by allowing a laser beam to strike a silver backlighter which will produce x-rays in the 3-4 keV range. These x-rays will travel the axial length of the cylinder and enter a pinhole camera called a Gated X-ray Imager (GXI).

The LANL group has done experimental and computational analyses with cylindrical geometry. This has helped gain understanding of the R-T instability in converging geometry, and the ablative effects. In this section RAGE (Radiation Adaptive Grid Eulerian) will be used to calculate various properties of the implosion process. These data will then be compared to the experiment and calculations by HYADES. HYADES is a one-dimensional Lagrangian hydrodynamics code that uses multigroup diffusion for the radiation transport. RAGE, as mentioned above, uses a Planckian radiation source. The cylinder that was used in the experiments at Lawrence Livermore National Laboratory with the NOVA laser is shown in Figure 4.6.

---

---

The cylinder shown in Figure 4.6, as mentioned earlier, is mounted transversely in a hohlraum, and eight NOVA laser beams are pointed at the wall of the hohlraum. The laser beam power is shaped in time. The pulse shape that was used for this experiment is designated PS-22, which stands for pulse shape 22. The power as a function of time is shown in Figure 4.7.

---

---

The materials used in the design of the cylinder were chosen to serve very specific purposes. The design, simple as it may look, is actually quite complicated and precise. The materials that was chosen for the central core of the cylinder was polystyrene foam. This was chosen because it had a low opacity, thus the backlighter x-rays would be able to travel the length of the cylinder without many interactions. The foam is also easy to fabricate and can be generated in a wide range of densities. For this experiment a density of 60 was chosen for two reasons. First, the low density would allow a high convergence ratio which will allow the R-T instability to be modeled to the highly nonlinear regime. Second, the low density foam will have a low enough opacity even at high convergence that a significant enough number of x-rays from the backlighter will get through the material without interaction. The reason this material was chosen to be 800 μm, instead of the entire length of the cylinder, is so that the attenuation of x-rays would not be significant.

Outside the cylinder, a thin 4 μm marker layer of C₆H₈Cl₂ was placed. This layer was axially centered and had a total axial length of 160 μm. This material has the purpose of creating a dark spot for the GXI. This material has a density of 1.42 and is opaque to the backlighter. Thus a dark image will appear in the GXI image where x-rays have to travel through this material. The choice of using only 160 μm was so that any parallax caused by the placement of the GXI pinhole would not attenuate the x-rays in the entire region of interest. Parallax is seen when the camera is at a slight angle to the axial direction of the cylinder. It would be impossible to setup the GXI such that there is no parallax because the pinholes are fixed in space and the camera is place such that there is a 2.7 degree offest for the four corner pinholes of the GXI.

The layer outside of the marker layer is polystyrene plastic. This is at a density of 1.044 . This is the material that is ablating away during the implosion process. This material was chosen because of ease of fabrication and because the large density was designed so that most of the x-rays from the gold will be deposited in this material.

The plastic, as shown in Figure 4.6, is tappered out to a radius of 315 μm to prevent any material from collapsing into the vacuum region until late times. The regions axially past the tapper are covered with 5 μm of gold. The gold serves two useful purposes. These are to keep the hohlraum hotter by reflecting and emitting x-rays which also keep the implosion more symmetric, and also to prevent x-rays from the hohlraum wall from preheating the plastic and imploding into the vacuum until late times.

Experimental results from this cylinder design are included in the Appendix A. The pictures shown there are the 12 shots from the GXI. From these shots both the inner and outer edge of the marker layer are quite clear. This gives a good image to process and see the behavior of the implosion process. It can be seen that the the marker layer expands as a function of time, for the region of time the shots were taken in. This is due to the fact that the marker layer will absorb a large amount of radiation energy, thus causing this region to heat and expand. It is thought that the marker layer will also expand at the early times, but there is no experimental data from these times. The ideal is based on the fact that the marker layer will absorb significant amounts of radiation that is not attenuated in the outer polystyrene plastic. This will cause a temperature rise and an expansion. This theory is backed up by computational data from HYADES. 4.3.1

4.3.1 Analysis of Computational Results for the Original Cylinder Design

A comparison of the RAGE model of the cylinder is shown with the experimental cylinder in Figure 4.8. The Figure shows the r-θ plane of the two cylinders. In the RAGE model the cylinder is not two dimensional but is instead one dimensional in the r-direction with θ-symmetry imposed.

---

---

The hohlraum and time-dependent laser energy deposition were simulated by depositing a time-dependent temperature profile in a low density air outside the outer layer of polystyrene. The outer air has a density of 1 × 10^-4. The time-dependent temperature profile that is deposited in the low density air follows that of pulse shape 22, which is shown in Figure 4.9. The temperature is scaled such that the peak temperature is 165 eV.

---

---

The results from this computation were looked at in both a qualitative and quantitative manner. The qualitative results of the analysis are shown in Appendix A, Original Cylinder Design. Quantitatively, the most important feature of the computation is the minimum inner radius of the marker layer and the time at which that minimum occurred. The results from the RAGE calculations are compared to both experimental data and a Lagrangian hydrodynamics code HYADES. The results of this comparison can be seen in Figures 4.10-4.13.

---

---

---

---

---

---

---

---

The radiation diffusion that is used by RAGE is not adequate to model the implosion process. This is because for optically thick, opaque, materials no preheating of the target will be experienced. There are ways to over come this, e.g., by using a multigroup diffusion. This would allow the preheat phenomena. This is not able to be done with RAGE, so instead one must define an artificially higher temperature in the marker layer to simulate the preheating. Notice that the marker layer is more dense than the polystyrene, and the marker layer contains chlorine. The chlorine will make the marker layer opaque to x-rays coming from the gold, thus a significant portion of the x-rays that made it through the outer polystyrene will deposit their energy in the marker layer. This is the reason that one must choose to simulate preheating of the marker layer by adding an artificially high temperature in the marker layer at time t=0. This second choice is what was chosen for the RAGE model of this cylinder. It is easy to see from the plot of material radius vs. time that the preheat simulation of the marker layer is effective in causing the marker layer to spread. The unfortunate aspect is that without the multigroup diffusion it is not possible to get the marker layer to spread again at later times. This is the most significant limitation of using the Planckian radiation diffusion. This later time spreading of the marker layer will be important in the modeling of the ablative R-T instability. The feature of the implosion that looks great is the inner marker layer as a function of time as illustrated in Figure 4.11. From the two plots it is seen that this follows both the experimental data and HYADES, a Lagrangian multigroup radiation diffusion hydrodynamics code. In Figure 4.12 the thickness of the marker layer as a function of time is shown. This shows that at early time the simulation of the preheat is doing a good job of causing the marker layer to spread. But at about 1 ns the first weak shock is seen to pass through the marker layer, and at this point the thickness of the marker layer decreases quickly for the RAGE calculation but stays constant for the HYADES calculation. At about 2 ns the marker layer is seen to expand due to converging effects and some heating. This expansion from the RAGE data is not as significant as that seen with the HYADES code, nor does it compare well with the experimental data. The conclusion that this result draws upon is that in order to model the ablatively driven implosion accurately with a hydrodynamics code one will need to implement multigroup radiation diffusion.

A unique feature of the RAGE algorithm, as mentioned in Chapter 3, is the AMR. The importance of this feature is most notable in large computations, and its effectiveness is seen in this calculation. The computational grid for this calculation starts out as 250 level 1 mesh of size equal to 0.0008 cm. The mesh are allowed to refine to a smallest size of 0.00005 cm, which will be 4 levels of refinement. Figure 4.14 shows this AMR feature of RAGE by displaying the levels of mesh refinement as a function of time and radius. The refinement is seen to follow the line of steep gradients in flow properties. In Appendix B this Figure is repeated on a larger plot and the other flow properties of specific internal energy, velocity, temperature, pressure, and density are included. Thus one is able to compare the gradients of these properties to see that in fact the AMR feature is working quit well.

---

---

Appendix B contains mesh plots of the various properties that were calculated by RAGE. The plots are three dimensional plots of the hydrodynamic property vs. radius and time. The properties that are looked at are temperature, pressure, specific internal energy, density, grid refinement, sound speed, and material speed. The maximum values of these properties are listed in the following table.

All of the peak values of fluid properties were found to be at time t ≈ 2.5 ns and at radius r=0. The exception was the velocity which was found to be maximum at time t ≈ 2.5ns but the radius was slightly larger than 0. The reason for this is that the wall boundary condition requires that the velocity goes to zero at the boundaries of the computational grid.

Chapter 5
Conclusions and Recommendations

The primary goal of this thesis was to establish the need to use Eulerian hydrodynamic codes to model inertial confinement fusion (ICF) phenomena, then to verify the fact that an Eulerian hydrodynamic code will in fact be able to adequately model the ICF implosion phenomena. In this chapter conclusions of the work performed is presented, followed by recommendations and future work that may need to continue in this area. 5.1

5.1 Conclusions

The computational capability to model the ICF phenomena will be crucial to understand the hydrodynamic instabilities in ICF. It was shown that a Lagrangian hydrodynamic model was not adequate to describe the strong rotational flow patterns in ICF. Because of this difficulty, an Eulerian grid which will be able to model the strong rotational flows without the grid lines crossing each other has been developed.

There have been Eulerian hydrodynamic codes used to model the ICF phenomena, but for the most part they have come short of the goal in predicting exact hydrodynamic responses in ICF plasmas. For this reason there is still a need for an Eulerian hydrodynamic code that can model the implosion phenomena of ICF accurately. The algorithm that Radiation Adaptive Grid Eulerian (RAGE) uses to deal with discontinuous flows, which are often seen in ICF, is better than some of the algorithms developed in the past. The Godunov method was shown to deal quite well with the shock tube problem. It was shown in chapter 3 that RAGE reproduced the same results as the analytic solution from the simple shock tube. Also the Godunov method had no post shock oscillations which are experienced by some hydrodynamics algorithms.

The Adaptive Mesh Refinement (AMR) feature of RAGE was shown to work well in reducing the computational grid size, and thus the memory needed by the computer. Because the time step for each cycle is determined by the smallest time step on the grid, each cycle has the same time step as a fixed grid. However, the adaptive mesh has a smaller number of total cells so each cycle takes less computational time and, thus, the problem runs faster than the fixed grid. For the Eulerian hydrodynamics codes an adaptive mesh is imperative.

The RAGE was used to model an experiment performed to determine the growth rate of the R-T instability in cylindrical coordinates. It was found that RAGE did well for smaller wave numbers and large Atwood numbers, but was not as good for the larger wave numbers. It was indicated that the reason for this discrepancy in the large wave numbers was because the sinusoidal perturbations, which are used to model this instability, are approximated by a square grid. In a Lagrangian code the material boundaries and mesh boundaries are the same but in an Eulerian code this is not true at all times. Thus the Eulerian code was less accurate at early times. The RAGE model was a one-dimensional model of the radial direction. The results from this model showed that a Planckian radiation source could be used to model the radiation diffusion for the overall implosion process. The RAGE predicted the radial implosions which are in good agreement with experimental observation. The code, however, fell short in modeling the thickness of the marker layer. The inadequate results from the marker layer were due to the fact that multigroup radiation diffusion was not used.

It is understood that the Lagrangian hydrodynamic code is not adequate for modeling ICF phenomena because of the inability to deal with rotational flow. Also the Eulerian hydrodynamic code alone may not be a practical tool because of the large computational resources needed by this approach. The solution to this dilemma may be the coupling of the better aspects of each of these methods. This could be done by allowing the hydrodynamics code to run in Lagrangian mode until the mesh starts to tangle. Once the mesh have bowtied, one can back up the computation a few cycles to the point at which the grid was still good in order to run the code in Eulerian mode. Doing this will produce the most efficient computer mode for the flow pattern in both early and late times.

5.2

5.2 Recommendations

The main goal of RAGE as a hydrodynamics code is to be able to model the three-dimensional hydrodynamic phenomena with full physics. As of this date both numerical models and large enough computers are unavailable. The Advanced Strategic Computing Initiative (ASCI) computer being built by Intel will be a 1 tera flops machine. This machine will be used to solve some of the complicated three-dimensional problems.

The future work of RAGE will be to develop a two-dimensional (R-θ) model for the purpose of modeling the growth of perturbations on the material interface. This model will then be extended to three dimensions.

It is the recommendation of the author that RAGE be used in conjunction with a Lagrangian hydrodynamic code to develop accurate modeling capabilities. The application of coupling the two methods can also be implemented in other areas of Computational Fluid Dynamics (CFD) such as atmospheric modeling, and two phase flow.

Chapter A
Computational Results From RAGE for the Original Cylinder Design

---

---

---

---

---

---

---

---

---

---

---

---

---

---

Chapter B
Hierarchical Data Format (HDF) Raster-Images of the Rayleigh-Taylor Instability Problem

---

---

---

---

---

---

---

---

---

---

---

---

Chapter C
Images from the Experiment

---

---

Chapter D
Code Description of the RAGE

The RAGE algorithm has a steep learning curve, and the following comments have been generated to ease the difficulty. The input deck can be divided into three main categories:

• Control; these are the first group of commands and are used to control RAGE during run time.
• Setup; these are what defines your material, regions, eos, and other problem generation commands.
• Output; these are the commands that control the output of RAGE. There are two main output data: HDF files, and tek files.

The sample input deck below has most of the features of RAGE in it. Since RAGE is in a high state of development there is no implication that all or any of the features of this input deck will be available on the code in the future.

The RAGE flow will be discussed to give the reader a more clear understanding of the RAGE problem solving approach. Because the RAGE flow is linear with only a few decision points, an ordered proceeder, will be presented.

1. The input file is read in and checked for correctness; if the input is correct RAGE continues. If the input is not correct RAGE will exit at this point.
2. Problem generation; once the input is read in, RAGE will start to generate the problem by reading in eos files, forming a level 1 grid and filling materials into there respective regions. If all definitions are correct to this point, RAGE will continue on, if not RAGE will exit gracefully at this point.
3. Grid adjustment and refinement; at this point RAGE will examine the grid and refine it if needed. On the first cycle RAGE will also deal with cells that have mixed materials.
4. Calculation of the next time step for the entire grid is now performed. The time step for the entire grid is determined by the slowest individual time step on the grid. If the time step is less than dtkill, RAGE will exit at this point, if not RAGE will continue.
5. Hydrodynamic equations are now solved.
6. Radiation diffusion is now solved if the ‘dorad flag’ is set to true.
7. At this point the time and cycle numbers are now checked to see if any output files need to be generated. Most of the output files can be dumped on cycle numbers or time. RAGE supports HDF format and TEK plots as well as ascii format output dmp-files.
8. The grid is now reexamined to see if any adjustment will be needed. If any adjustment is needed it is done now.
9. At this point RAGE checks to see if any flags are set, if so they are dealt with by exiting if needed. If no flags are set, RAGE will cycle back to step 4, and continue until the calculation is complete or RAGE is force to exit.

LIST OF REFERENCES

[1] J J Duderstadt and G A Moses. Inertial Confinement Fusion. John Wiley and Sons, New York, 1982.
[2] Tom Dolan. Fusion Research. Pergamon Press, New York, NY, 1982.
[3] J Nuckolls, L Wood, A Thiessen, and G Zimmerman. Laser compression of matter to super-high densities: Thermonuclear (ctr) applications. Nature, 239:139–142, 1972.
[4] W J Krauser et al. Ignition target design and robustness studies for the national ignition facility. Physics of Plasmas, 3(5):2084–2093, 1996.
[5] H A Bethe. Energy production in stars. Physical Review, 55:434–456, 1939.
[6] R E Kidder. Laser-driven isentropic hollow-shell implosions: The problem of ignition. Nuclear Fusion, 19(2):223–234, 1979.
[7] N A Tahir and K A Long. Numerical simulation and theoretical analysis of implosion, ignition and burn of heavy-ion-beam reactor size icf targets. Nuclear Fusion, 23(7):887–916.
[8] D B Kothe, J U Brackbill, and C K Choi. Implosion symmetry of heavy-ion-driven inertial confinement fusion targets. Physics of Fluids B, 2(8):1898–1906, 1990.
[9] Yuli Pan. Reducing the convergence ratio of high gain icf targets. Nuclear Fusion, 28(3):363–368, 1988.
[10] J D Kilkenny and et al. A review of the ablative stabilization of the rayleigh-taylor instability in regimes relevant to inertial confinement fusion. Physics of Plasmas, 1(5):1379–1389, 1994.
[11] J B Beck, W W Hsing, N M Hoffman, and C K Choi. Experimental analyses of rayleigh-taylor growth in cylindrical implosions. LA-UR-95, 1995.
[12] K O Mikaelian. Lasnex simulations of the classical and laser-driven rayleigh-taylor instability. Physical Review A, 42(8):4944–4951, 1990.
[13] H Sakagami and K Nishihara. Rayleigh-taylor instability on the pusher-fuel contact surface of stagnating targets. 2(11):2715–2730, 1990.
[14] Lord Rayleigh. Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density. Proceedings of the London Mathematical Society, 14:170–177, 1883.
[15] Sir Geoffrey Taylor. The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. Proc. Roy. Soc., 201:192–196, 1950.
[16] S Chandrasekhar. Hydrodynamic and Hydromagnetic Stability. Dover, New York, 1961.
[17] G Bateman. MHD Instabilities. MIT Press, Cambridge Massachusetts, 19.
[18] C K Choi. Class Notes Nuclear 697M Inertial Confinement Fusion. Unpublished, Spring 1996.
[19] R L McCrory, L Montierth, R L Morse, and C P Verdon. Nonlinear evolution of ablative-driven rayleigh-taylor instabilities. 46(5):336–339, 1981.
[20] H Takabe et al. Self-consistent growth rate of the rayleigh-taylor instability in an ablatively accelerating plasma. Physics of Fluids, 28(12):3676–3682, 1983.
[21] H J Kull and S I Anisimov. Ablative stabilization in the incompressible rayleigh-taylor instability. Physics of Fluids, 29(7):2067–2075, 1986.
[22] V V Byshkov et al. Self-consistent model of the rayleigh-taylor instability in ablatively accelerated laser plasma. Physics of Plasmas, 1(9):2976–2986, 1994.
[23] C P Verdon et al. Nonlinear effects of multifrequency hydrodynamic instabilities on ablatively accelerated thin shells. Physics of Fluids, 25(9):1653–1674, 1982.
[24] M H Emery, J H Gardner, and S E Bodner. Stronly inhibited rayleigh-taylor growth with 0.25-μm lasers. Physical Review Letters, 57:703–706, 1986.
[25] M H Emery, J P Dahlburg, and J H Gardner. The rayleigh-taylor instability in ablatively accelerated targets with 1, 1/2 and 1/4 mum laser light. Physics of Fluids, 31(5):1007–1016, 1988.
[26] B A Remington et al. Laser-driven hydrodynamic instability experiment. Physics of Fluids B, 7(5):2589–2595, 1993.
[27] B A Remington et al. Large growth rayleigh-taylor experiments using shaped laser pulses. Physical Review Letters, 67(23):3259–3262, 1991.
[28] S G Glendinning et al. Laser-driven planar rayleigh-taylor instability experiments. Physical Review Letters, 69(8):1201–1204.
[29] B A Remington et al. Large growth, planar rayleigh-taylor experiments on nova. Physics of Fluids B, 4(4), 1992.
[30] M Desselberger and O Willi. Measurement and analysis of rayleigh-taylor instability in targets driven by incoherent laser radiation. Physics of Fluids B, 5(3):896–909, 1993.
[31] M Desselberger, O Willi, M Savage, and M J Lamb. Measurement of the rayleigh-taylor instability in targets driven by optically smoothed laser beams. Physical Review Letters, 65(24):2997–3000, 1990.
[32] R D Richtmyer and K W Morton. Difference Methods for Initial-Value Problems. John Wiley and Sons, New York, NY, 1967.
[33] K O Mikaelian. Freeze-out and the effect of compressibility in the richtmyer-meshkov instability. Physics of Fluids, 6(1):356–368, 1994.
[34] K O Mikaelian. Richtmyer-meshkov instabiliteis in stratified fluids. 31(1):410–419, 1985.
[35] I G Currie. Fundamental Mechanics of Fluids. McGraw-Hill, New York, 1993.
[36] M G Lumsden. RAGE User Manual. SAIC, 1994.
[37] P Woodward and P Colella. The numerical simulation of two-dimensional fluid flow with strong shocks. Journal of Computational Physics, 54:115–173, 1984.
[38] S K Godunov, A V Zabrodin, and GP Prokopov. A computational scheme for two-dimensional non-stationary problems of gas dynamics and calculations of the flow from a shock wave approaching a stationary state. Zh. Vych. Mat., 1(6):1187 – 1218, 1961.
[39] F H Harlow and A A Amsden. Fluid dynamics. LA-4700, pages 58–71, 1971.
[40] Y B Zel’dovich and Y P Raizer. Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena. Academic Press, New York, 1966.
[41] Nelson M Hoffman. Personal Comunications. Los Alamos, NM, Summer 1995.