Radiation hydrodynamics simulation of high-Z doped ICF targets

We conducted a comparative study among the Bi-CG family with incomplete LU factorization employed for solving diffusion-type radiative transfer equation in radiation hydrodynamics (RHD) simulations for laser-produced plasmas. The result suggests that Bi-CGSTAB is the most suitable method thanks to its high convergence stability and low computational cost. RHD simulations were performed for Rayleigh–Taylor instability experiments of an inertial confinement fusion target. We found that higher harmonics of static pressure perturbation become important at an intermediate wavelength.


Introduction
Suppression of Rayleigh-Taylor (RT) instability [1] at the ablation surface of a fuel shell is still a critical issue for inertial confinement fusion (ICF) [2] despite the fact that continuous efforts have been devoted to perform it for a long time. On a frame of the ablation surface, the low density plasma corona sustains massive fuel plasma in a gravitational field. Hence the surface becomes RT unstable while the ablated material flows through it. In classical RT instability which occurs at a material interface of an incompressible fluid, an analytical growth rate of the interface γ cl is obtained by linear theory: where ρ h (ρ l ) is the density of higher (lower) material, k is the wavenumber of the given perturbation, and g is the acceleration. Perturbation growth in an ICF target is suppressed due to the ablative flow through the unstable surface. The growth rate γ can be expressed by the modified Takabe formula [3,4] where α and β are the constants determined by the ablation structure, L is the density scale length at the ablation surface, and v a is the ablation velocity.
A concept of indirect-direct hybrid drive is preferable for achieving both target stability and high energy coupling efficiency. A high-Z doped target was proposed in the context of the hybrid drive by utilizing the double ablation structure formed by self-emission of high-Z dopant [5]. Radiation hydrodynamics (RHD) simulations demonstrated that perturbation growth in the high-Z doped target is suppressed because the radiative ablation (RA) surface is stabilized by its high ablation rate. However, the detailed property of such a target should be assessed for practical use; in particular the distance between the RA front and the electron-conduction ablation (EA) one may govern the behavior of a specific wavelength [6].
Although analysis of the perturbation growth can be performed by a multi-dimensional RHD code, reliable simulations should be conducted with a large number of computational cells for resolving the ablation structure (< µm) in the acceleration area (∼ mm) and so they have a large computational costs. Since the radiative transfer solver is the most time-consuming module in the RHD code, it is important to make it efficient. If an optically thick flowfield can be assumed, the radiative transfer equation is simplified to a diffusion-type equation. We usually employ the bi-conjugate gradient (Bi-CG) method [7] for obtaining a solution of such an equation because the coefficient matrix that appears in implicit differentiation with 9-point stencil becomes asymmetric. In this study, the coefficient matrix becomes symmetric since we use Cartesian grids with isotropic radiation. However, we investigate the convergence properties of the Bi-CG family to develop a versatile code.
In this paper, we first conducted a comparative study among methods of the Bi-CG family. The Bi-CG method, conjugate gradient squared (CGS) method, bi-conjugate gradient stabilized (Bi-CGSTAB) method, and Bi-CGSAFE method are examined for the comparison with test matrices obtained from radiation transfer problems of a practical RHD simulation. Then, we performed a set of RHD simulations corresponding to planar ICF target experiments and analyzed the flowfield to discuss a significant suppression appearing in intermediate wavelengths.

Numerical methods and conditions
Numerical simulations were carried out by a two-dimensional Eulerian RHD code, RAICHO [8]. The code solves the two-dimensional Euler equations including some energy sources with realistic equations of state (EOS) for electrons and ions. The EOS for ions is based on the Cowan model and that for electrons is based on the Thomas-Fermi model. In the case of low-density and high-temperature plasma, the electrons and ions are in thermal non-equilibrium. Hence the energy equation is divided into two energy equations for electrons and ions. One of the source term is estimated by solving the radiative transfer equation which is approximated in a diffusion-type form.
The governing equations were discretized in a cell-centered finite volume fashion. The AUSM-DV method [9] was employed for the numerical flux with the second-order MUSCL method [10]. The dependent variables were integrated with a first-order explicit Euler method with inviscid flux and with an implicit method with the source terms in an operator splitting manner.
The laser is assumed to be normally irradiated on the planar target, and the laser absorption of inverse-bremsstrahlung is calculated with one-dimensional ray-tracing. The thermal conduction is formulated by a diffusion approximation with the Spitzer-Härm conductivity and a flux limiter of 0.1 [11]. In this study, the number of energy groups is 32 for the radiative transfer equation under a multi-group approximation. wavelength, and 1.9-ns FWHM. The initial target density is 1.26 g/cm 3 , and its thickness is set to 16 µm for all the simulations. The critical density of the target plasma is 3.98 × 10 21 cm −3 . Figures 3 and 4 show examples of time history of relative residual norm for the Bi-CG, CGS, Bi-CGSTAB, and Bi-CGSAFE for two cases, in the radiation transfer computation. All methods are preconditioned by incomplete LU factorization vectorized in a hyperplane manner. Test matrices of the diffusion coefficient were obtained from the radiation transfer process of the practical RHD simulation. Case-A is a typical example for the early stage of laser incidence at which the radiation field drastically changes due to initial ionization. Case-B is for a matrix when the input laser becomes constant and the ablation structure is in a quasi-steady state. In Case-A, the matrix becomes stiffer than in Case-B, so a larger number of iterations is required to obtain a solution. Table 1 shows numbers of iterations and CPU time for solving Case-A and Case-B with each method. The numbers of iterations for the Bi-CGSTAB and the Bi-CGSAFE are much smaller than that of the others because of their high convergence stability. Although the Bi-CGSAFE is more state-of-the-art than the Bi-CGSTAB, the former needs more CPU time than the latter. This difference originates from the low computational cost of the Bi-CGSTAB.    Hence the Bi-CGSTAB is the most suitable method among the Bi-CG family especially for the radiative transfer problem. 4. Perturbation growth Figure 5 shows that a significant suppression of the RT growth rate in the intermediate wavelength was found as predicted in [6]. In order to identify the cause of such a suppression, we defined static pressure perturbation between two ablation fronts in an acceleration frame of the EA surface as follows:

Convergence property in radiative transfer solver
where p(x, y) is the static pressure in the EA acceleration frame, and λ is the perturbation wavelength. Thus the static pressure perturbation means the difference between static pressure    and its average value along y-axis at a given x-position in this paper. In the case of λ = 30 µm, at which the perturbation growth is significantly suppressed, higher harmonics of the static pressure perturbation appear between the RA front and EA front (figures 6-8). The fundamental mode is damped for this wavelength and then cannot grow at the RA front deformation. There is some possibility that a suppression effect emerges via pressure waves from the EA front to the RA front.

Conclusions
We conducted a comparative study among methods of family of conjugate gradient. The results suggest that the Bi-CGSTAB is the most suitable method among the Bi-CG family on RHD simulations. We also investigated the dynamics of ablation fronts in a high-Z doped ICF target using a two-dimensional RHD code. We found that higher harmonics of the static pressure perturbation between two ablation fronts develop at intermediate wavelengths, at which a significant suppression appears.