INTRODUCTION
The behavior and performance of inertial fusion targets are calculated with complex numerical simulations. These simulations have evolved over the decades and have incorporated physical processes that are important, but are not always fully understood. The simplest one-dimensional (1D) (i.e., radial) fluid models incorporating only classical transport, generally show robust high gain for well designed targets. As complicating effects such as two-dimensional (2D) fluid instabilities (Rayleigh-Taylor and Richtmyer-Meshkov) (Weber et al., Reference Weber, Glendinning, Kalantar, Key, Remington, Rotheneberg, Wolfrum, Verdon and Knauer1997; Schmitt et al., Reference Schmitt, Colombant, Velikovich, Zalesak, Gardner, Fyfe and Metzler2004; Radha et al., Reference Radha, Goncharov, Collins, Delettrez, Elbaz, Epstein, Glebov, Goncharov, Keck, Knauer, Marozas, Marshall, McCrory, McKenty, Meyerhofer, Regan, Sangster, Seka, Shvarts, Skupsky, Srebro and Stoeckl2005), flux-limitation (Malone et al., Reference Malone, McCrory and Morse1975), nonlocal electron transport (Luciani et al., Reference Luciani, Mora and Virmont1983; Epperlein & Short, Reference Epperlein and Short1991; Sunahara et al., Reference Sunahara, Delettrez and Stoeckel2003; Manheimer & Colombant, Reference Manheimer and Colombant2004), and laser plasma instabilities (Kruer, Reference Kruer2000) are included, the problem becomes more complicated and the viable regions of parameter space generally shrink. One effect not included in the fluid simulations is ion viscosity.
However, a recent calculation (Li et al., Reference Li, Igumenshchev and Goncharov2006) included the effect of viscosity on the momentum equation (but in their original abstract, had not yet included it in the energy equation) and found that near the peak of the implosion, the effect of it could be important. We have examined this as well, using a somewhat different approach. We perform a fluid simulation on a laser implosion without ion viscosity, and then post-process the data to determine where viscosity could be important. We look not only at the dynamics as Li et al. (Reference Li, Igumenshchev and Goncharov2006) did, but we look at two other issues. First we confirm the standard result, namely that viscosity is nowhere near important enough that one can dispense with shock capturing algorithms in the fluid simulation. Second we examine the effect of the viscosity on the Rayleigh-Taylor instability. Using a simple slab model, we find that the effect on the growth rate is negligible at all times except right before the peak of the implosion, where it has a significant, but certainly not a dominant effect. In that sense, we confirm the results of Li et al. (Reference Li, Igumenshchev and Goncharov2006) that for times just near (but not after) the peak of the implosion, viscosity could be playing a significant role. However for all other times, it has virtually no effect on the dynamics.
EVALUATION OF THE EFFECTS OF ION VISCOSITY FOR A TYPICAL LASER IMPLOSION AND THE NEED FOR SHOCK CAPTURING ALGORITHMS
We post-process our simulation of the 0.5 MJ laser fusion target (Obenschain et al., Reference Obenschain, Colombant, Schmitt, Sethian and McGeoch2006), to see where viscosity could be playing a role. The target and laser pulse characteristics for this target are shown in Figure 1. The target consists of 153 µm of deuterium-tritium (DT) fuel surrounded by a 146 µm DT-wicked foam ablator and a 5 µm CH outer layer. The laser pulse consists of a 3.4 ns foot (3.1 TW) followed by a gradual ramp-up to full power (178 TW) at t = 7.33 ns that is maintained for 2.1 ns. Zooming also occurs in two stages as the target radius shrinks. The yield for this target is 30 MJ for incident laser energy of 480 kJ, leading to a gain slightly above 60. The implosion diagram for this pellet, as well as the inward velocity, as a function of time is shown in Figures 2a and 2b.
Fig. 1. (Color online) Schematics of the 1/2 MJ target and its corresponding laser pulse.
Fig. 2. (a) r-t diagram for target and laser pulse shown in Figure 1 and (b) average maximum implosion velocity for this target.
We use the simplest theory of bulk viscosity for the 1D spherical implosions, assuming that in the heat front, the temperature, and density gradient scale lengths are much less than the radius. We will see shortly that this is a very good approximation. Hence, the configuration is nearly planar, so we use for the divergence of the viscous plus ion pressure stress tensor the quantity (Braginskii, Reference Braginskii1965; Huba, Reference Huba2006)
Thus, the ratio of the viscous to ion pressure stress tensor is simply 1.28 τii(∂v/∂r). Here τii is the ion–ion collision time, given by
α is the ratio of ion mass to proton mass, and Λ is the Coulomb logarithm. Note that since these calculations apply to laser fusion targets, the Z's are rather small, usually unity in the fuel or just over unity in the foam ablator, and always less than 3.5 as appropriate for a CH layer. For other types of target at higher Z, ion viscosity will be less important generally.
We selected three times during the laser implosion, 5, 7, and 9.5 ns and plotted out as a function of radius, a variety of quantities shown in Figures 3–5. Panels a are the electron density, panels b are the velocities, panels c are both the electron and ion temperature, and panels d are the magnitude of the ratio of the viscous to ion pressure tensor. At 5 ns, this ratio is small everywhere except right at the sharp density drop at the inner edge of the pellet. However, everywhere else, the effect of the viscosity is on the order of 10−3. At 7 ns, the effect is similar, but now it is generally on the order of 10−2, except for the large value at the inner edge of the pellet. In these regions of large viscosity, its effect would be to somehow smooth out the transition in the velocity profile between the accelerating shell and the inner low density plasma. However, the viscosity is so large here, that in this small region, the entire concept of a fluid formulation is questionable. But this is not a very important region of the plasma at all; the density and velocity are both very small in the region of large viscosity (or more accurately, the region where a fluid formulation is suspect). At 9.5 ns, we find that the effect of viscosity maximizes. In the inner region, it is an order unity effect, while in the outer regions; it is still an order 10−2 effect. In order to get a better idea of the effect of the viscosity at this time, we plotted in Figure 6a, the total electron and ion pressure on a greatly expanded scale horizontally, and in Figure 6b, the electron plus ion pressure plus the viscous stress tensor. The force on the fluid is the negative gradient of this quantity. The quantities are plotted as a function of grid cell up to grid cell 100 (the simulation has about 500 grid cells, and grid cell 100 is at about r = 0.024 cm where the density is n e = 1.7×1023 cm−3 and the inward velocity is v = 4.4×107 cm/s). Clearly at around this time, and this time alone, the viscosity could be playing a non-negligible role in a significant portion of the plasma. Its effect will be to smooth out the velocity profile somewhat for this inner portion of the implosion. However at 10 ns, just after peak compression where all the flow is outward, the viscosity is everywhere a correction on the order of 10−2.
Fig. 3. (a) Electron density, (b) velocity, (c) ion and electron temperature profiles, and (d) ratio of ion viscosity to total pressure versus radius at t = 5 ns (during compression phase).
Fig. 4. Same as Figure 3 but at t = 7 ns, just slightly ahead of shock break-out. The discontinuity in ion temperature profile occurs at the ablator/CH interface and ensures continuity of the total pressure at this location.
Fig. 6. Differences between the total pressure (a) without and (b) with ion viscosity as a function of grid cell number at t = 9.5 ns. Effects of the ion viscosity are seen to take place very close to the center of the target where the ion temperature is larger than the electron temperature.
We now discuss whether viscosity can ever be so important that we do not need shock capturing algorithms in our fluid simulations. Since a strong viscous shock, of the type we utilize in a laser implosion, has a thickness on the order of the ion mean free path, we simply post-process the data of the implosion to calculate the ratio of Δx/λii as a function of x for a variety of times. Here Δx is the grid spacing and λii is the ion mean free path. Without displaying the graphs, we simply summarize the result, namely that this ratio is very large, on the order of hundreds or thousands in virtually all regions of the plasma and at all times. Thus, viscosity will never sufficiently smooth the profile that shock capturing algorithms are unnecessary.
EFFECT OF VISCOSITY ON THE RAYLEIGH-TAYLOR INSTABILITY
We consider a simple slab model of the Rayleigh Taylor instability. Gravity g points downward. A fluid with density ρ1 and kinematic viscosity η1 (i.e., η has dimension of a length squared divided by time), is in equilibrium on top of a lower density fluid with density ρ2, and viscosity η2. Looking at the plots of density as a function of space in Figures 3–5, the slab model seems reasonable due to the very large drop in density at the accelerating surface. It is also reasonable, because, as we will see, the viscous correction to the growth rate depends almost entirely on the parameters of the heavier fluid, the fluid for which the uniform slab is almost certainly a reasonable approximation.
The theory of the Rayleigh-Taylor instability of viscous fluids is generally considered in specialized text books only, such as Chandrasekhar's (Chandrasekhar, Reference Chandrasekhar1961). In the general case, the calculation of the instability requires the numerical solution of a fourth order polynomial. One way in which this can be simplified is the so-called Hide approximation (Mikaelian, Reference Mikaelian1993; Piriz et al., Reference Piriz, Cortazar, Cela and Tahir2006) where the vertical dependence of the eigenfunction is taken as that for the non-viscous fluid. Where valid, this approximation can be quite accurate. However, in some cases, the eigenfunction in the vertical direction can depart from the non-viscous case in a significant way. In these cases, the Hide approximation is not valid and one must solve the dispersion relation in a more accurate way. We find this to be the case for laser produced plasmas. In the appendix, we solve for the dispersion relation in the appropriate limit. There, we worked out the theory of the instability, but here, we simply give approximate analytic expressions for the dispersion relation in several limits. In all cases, we assume that the effect of viscous damping is small compared to the basic growth rate (kg)1/2 and that ρ2<<ρl. Expanding the dispersion relation from the appendix in lowest power of ηk2/(kg)1/2, we find
![\gamma^2=\hbox{kg}\left[{1 - {2\rho_2 \sqrt{\eta_1 \eta_2 } \hbox{k} \over \lpar \hbox{kg}\rpar ^{1/4} \lpar \rho_2 \sqrt{\eta_2}+\rho_1 \sqrt{\eta_1}\rpar }}\right].](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151021081604112-0734:S0263034607000663_eqn3.gif?pub-status=live){kind=small}
Notice that the correction to the growth rate goes as a fractional power of the ratio of viscous damping rate to growth rate. Other approximations find that the correction goes as the ratio of these two rates (Mikaelian, Reference Mikaelian1993; Piriz et al., Reference Piriz, Cortazar, Cela and Tahir2006). However these calculations use a variational approach where they have to assume an eigenfunction, and then, in terms of the assumed eigenfunction, calculate the correction to the growth rate. But generally these assumed eigenfunctions are chosen to be the same as the eigenfunctions in the absence of viscosity. As we see in the appendix, the actual eigenfunctions are the sum of two parts, one the non-viscous part, and another, a boundary layer effect, dominated by the viscosity. It is this boundary layer effect that provides the dominant dissipation.
Note that in the Rayleigh Taylor instability without viscosity, the velocities of the fluids parallel to the interface, at the interface, is equal and opposite, giving rises to a strong viscous stress. In the presence of viscosity, the fluids are governed by a no slip boundary condition. This means that there must be a great deal of dissipation in the viscous boundary layer, implying a strong effect on the growth rate.
However, what happens if the top fluid is supported by a massless bottom fluid (essentially a vacuum, but one that supports the top fluid's pressure)? Then ρ2 vanishes in Eq. (3) and there is no correction to the growth rate on this order. However if there is no bottom fluid, a no slip and a no separation boundary is no longer meaningful and the problem must be reformulated. This is also discussed in the appendix. There it is shown that the dispersion relation becomes approximately
Hence if there is no strong shear flow at the interface, the reduction in growth rate is much less.
Now let us use our knowledge of the physics of a laser plasma implosion to get a simpler expression for the growth rates, valid in all limits, as long as the correction to the growth rate is small. In the laser implosion, pressure is approximately constant across the interface, certainly the variation in pressure is always far less than the variation in density or temperature, so ρ1Ti1 ~ ρ2Ti2. Also for either fluid, η is proportional to λii[Ti/M]1/2, and λii ~ Ti2/ρ. Using these relations, we find that in the denominator of Eq. (3), the low density term dominates by roughly a factor of [ρ1/ρ2]3/4. As is apparent from Figures 3–5, density ratios are typically about 20–25, so the low density term dominates by about an order of magnitude. But it cancels the low density term in the numerator, so Eq. (3) reduces to
Notice that the viscous damping is governed entirely by the viscosity of the high density fluid. This makes the theory more credible in the inhomogeneous fluid because the uniform density approximation is much better satisfied in the high density fluid than it is in the low density fluid. Thus we expect the slab model to be a reasonable approximation.
In our calculations of laser implosion, the code estimates Rayleigh-Taylor growth using various dispersion relations for the growth rate. With an expression for the growth rate as a function of mode number l (k = l/R ab, with R ab being the ablation radius), the post- processor calculates the growth rate at various times of the implosion. The calculation of ablative stabilization of the Rayleigh-Taylor instability is complicated and different theories (Bodner, Reference Bodner1974; Takabe et al., Reference Takabe, Montierth and Morse1983; Betti et al., Reference Betti, Goncharov, McCrory and Verdon1998) give somewhat different results. In the Naval Research Laboratory (NRL) 1D simulations, we typically apply several of these theories to get several different graphs of γ(l) where l/r ab = k and r ab is the radius of the ablation surface. In this way, we have several estimates (without performing a much more difficult and time consuming 2D simulation) of the effects of the Rayleigh-Taylor instability on laser fusion targets. The theory we use here and which is widely used in the community is that of Takabe et al. (Reference Takabe, Montierth and Morse1983). While this is regarded as sufficient for present purposes, it is of course only an approximation to the effect of the stabilizing effect of the dynamic overpressure created at the unstable interface. We consider the growth rate
where v abl is the ablation velocity. Eq. (6), without the curly bracket is similar to the standard expression for growth rate used in our 1D simulation to estimate the Rayleigh-Taylor growth. We plot the growth rate in s−1 of mode number one from Eq. (6) in Figure 7a. To get an idea of the reduction of the growth rate due to viscosity, we plot out the square root of the curly brackets as a function of mode number at a variety of times ranging between 5 and 9.5 ns as shown in Figure 7b. All times except the last are bunched up to make effectively a single graph for the reduction. At 9.5 ns, the effect of viscosity on the growth rate is considerably larger, but at all times the dominant stabilization mechanism is the ablative stabilization. Even assuming a 10% reduction in growth rate at l = 200, for half an ns, we see that the integrated growth is only reduced by about 0.1.
Fig. 7. (a) Rayleigh-Taylor growth rate from dispersion relation without ion viscosity at t = 8.5 ns and (b) correction to the growth rate when including ion viscosity as a function of mode number for various times.
CONCLUSIONS
For almost the entire time of the implosion of the NRL 0.5 MJ pellet, viscosity plays essentially no role. At the places in the plasma where it is important, just to the back of the fuel, the plasma has essentially no density or velocity. In fact, in this tiny region, viscosity is so important that the viscous stress is by no means a perturbation, but is the dominant effect. This calls into question the very idea of a fluid formulation. However, this region of the plasma is so unimportant to the overall implosion, that for this alone, it is not worth including a viscosity in the fluid formulation, and even if one did, it is not clear that it would provide any more accuracy. Similarly, the effect of viscosity on the Rayleigh-Taylor instability is so small for these times (it reduces the growth rates by at most a couple of percent) that calculations neglecting the viscosity should be reasonably accurate.
However for a very small time before the peak of the implosion, the ion viscosity has a significant, although by no means dominant effect on the dynamics and on the Rayleigh-Taylor instability. After the peak of the implosion, at 10 and 10.5 ns in our simulation, the effect of viscosity is again down to a percent or less. Target designers must consider whether it is worth the effort to include viscosity so as to more accurately model these effects which become significant just before the peak of the implosion. Of course if viscous stress is included in the momentum equation, then viscous heating should also be included in the energy equation.
Acknowledgements
We would like to acknowledge discussions with Dr. S. Zalesak. This work was supported by US Department of Energy/NNSA.
Appendix: Solution of a Model Problem for Viscous Damping
Here we assume a heavy fluid with density ρ1 and kinematic viscosity η1 on top, a gravity g pointing down in the negative z direction, and a light fluid with ρ2, η2 on the bottom. We assume a simple viscous model so we take as the linearized equations for each fluid
and
Assuming that the fluid quantities have spatial and temporal and spatial dependence going as exp[γt+ikx], Eqs. (A1 and 2) can be combined into a single equation for vz:
![\gamma \left[1 - \hbox{k}^{ - 2} {\partial ^2 \over \partial \hbox{z}^2}\right]\hbox{v}_{\rm z}=\eta \left(-\hbox{k}^2+{\partial^2 \over \partial \hbox{z}^2}\right)\left(1 - \hbox{k}^{ - 2} {\partial^2 \over \partial \hbox{z}^2}\right)\hbox{v}_{\rm z}](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151021081604112-0734:S0263034607000663_eqn9.gif?pub-status=live)
and the expression for the perturbed pressure p in terms of the z component of velocity is
There are four solutions to Eq. (A3),
and
Clearly only the upper signs are used in the lower fluid, and only the lower signs are used for the upper fluid, so as to have solutions which do not diverge. Thus, the problem is specified by four coefficients which we specify as v1a, v1b, v2a, and v2b. The former, Eq. (A5a) is the same as the solution without viscosity, while the latter Eq. (A5b) is, in the limit of small viscosity, a thin boundary layer surface mode. In the inviscid Rayleigh-Taylor instability, there is a strong shear at the boundary. That is the fluid motions parallel to the surface, just across the surface from one another, are in opposite directions. Thus when viscosity is present, and a no slip boundary condition is imposed, we expect that the viscosity dominated mode will play an important role. We will see shortly that this is the case. To find the dispersion relation, we need four boundary conditions relating the four coefficients. Clearly one is that there is no separation of the fluids, or v1z = v2z. For viscous fluids, there is also a no slip boundary condition, so at the interface v1x = v2x.
The other two boundary conditions come from the fact that the stress tensor must be continuous across the perturbed boundary between the two fluids. The ambient gravity gives a contribution to the stress tensor at the perturbed surface, all other components of the stress tensor arise only from perturbed quantities. In our simple model however, we take the viscous stress tensor as −ηρ∂vi/∂xj. Continuity of both components of the stress tensor across the perturbed boundary gives the result that
and
are continuous across the interface. Imposing these four conditions leads to a 4×4 determinant of the coefficients. Setting this determinant equal to zero gives the dispersion relation. This relation is
where
Eq. (A7) could be solved numerically for γ2/kg in terms of three-dimensionless parameters, ρ1/ρ2, η1/η2, and kg/η12 k4. However at this point, we prefer to get analytic insight instead. The determinant has a η to a variety of powers in the numerators and denominators. Assuming small η's, we can keep only the dominant term and solve for the vb's (the coefficient of the boundary layer solutions) in terms of the va's. The result is
and
Assuming that both the correction to the classical growth rate from viscosity is small, and furthermore that ρ2 << ρ1 (which as Figs. 3–5 show, is clearly a good approximation), then we find the correction to the classical growth rate is as given in Eq. (3). Notice that the reduction in growth goes as the square root of the viscosity, rather than the viscosity. The reason is that the dominant dissipation arises from the boundary layer viscous mode, as we expect, due to the strong shear at the interface in the absence of viscosity. Any theory, which assumes a mode proportional to exp±kz will, of course miss this effect. Also, Eq. (3) shows that the reduction in growth rate scales as the density of the lighter fluid. This is reasonable, because the lighter the lower fluid is, the less important this shearing motion is. In the ultimate limit of no lighter fluid at all, ρ2 = 0, this shear obviously does not exist at all, so it cannot exert a stabilizing effect on the instability.
In the limit of no lighter fluid at all, that is, the heavier fluid supported only by the pressure of the vacuum, it is a simple matter to calculate the dispersion relation. Since there is no fluid underneath, there is no need for the no slip and no separation boundary condition. There are only two solutions, v 1a and v 1b, the standard Rayleigh-Taylor and the boundary layer modes of the upper fluid. The boundary condition is that both components of the stress tensor must be continuous from the fluid to the vacuum. Imposing this condition, it is a simple matter to see that the dispersion relation, in the limit of ρ2 -> 0 is given by Eq. (4). In this case, where there is no shear motion between the top and bottom fluid, viscous stabilization is a much weaker effect.
