Influence of initial particle configuration upon sphericity of shock loaded density interface in SPH modelling of Richtmyer-Meshkov instability

The Richtmyer-Meshkov Instability (RMI) occurs when a perturbed interface between two fluids of different densities undergoes impulsive acceleration by, for example, a shock wave. Its study is important for gaining better understanding of processes governing inertial confinement fusion where a spherical solid deuterium-tritium (DT) capsule filled with the DT gas is compressed by a powerful laser [1]. The imploding shock moving to the center of the spherical shell causes the RMI to form on its inner surface due to a large density difference between solid and gaseous DT phases, which leads to a rapid growth of initially present randomperturbations. Usually, the density interface consists of random perturbations (e.g., surface roughness, machining inaccuracies, etc.) which initially grow linearly before nonlinear effects emerge. To improve reliability of shock-driven flow simulations, it is necessary to understand the factors which inuence the RMI growth. Our work aims to find the best initial particle configuration that does not introduce artificial perturbations of the spherical density interface. The shocked interface must preserve spherical symmetry if no perturbations are applied. In order to reduce the number of simulated particles in calculation, it is necessary to adapt spatial resolution, i.e., to increase the concentration of particles in target regions.


Problem statement
The method of Smoothed Particle Hydrodynamics (SPH) was proposed by Lucy, Gingold and Monaghan in 1977 for astrophysical problems. Libersky and Petschek in 1990 extended it to study the dynamic response of materials. The method is used to investigate such processes as free surface flows, explosion phenomena, impact and penetration, and instability growth.
The problem statement is shown in figure 1. All materials are modelled with an ideal gas EOS by the SPH method implemented in the code package MOLOCH [2].
In this work we investigate two types of initial conditions (IC): the lattice IC generated with Face Centered Cubic lattice (figure 2, left), and the spherical IC generated with concentric spheres (figure 2, right). Particle positions in the spherical IC are found with an algorithm from [3].

Two templates
Let us do calculation with no adaptation for a sample with the lattice and spherical IC. In a 1/8 of a sphere we define 10 million particles.
After shock wave passing through the interface one can see obvious differences in its sphericity for the different IC ( figure 3 and figure 4). Spurious perturbations are much higher for the lattice IC than for the spherical IC.   traced and the lattice of particles gives almost a uniform distribution over radii. However after loading, the spherical IC preserve sphericity, while the lattice IC depart from uniform which is indicative of broken sphericity. So, the spherical IC perform best and in what follows we are only working with them.  Number of particles in spherical layers after loading for the lattice and spherical IC.

Adaptivity
Adaptive IC are needed in order to minimize the number of particles required for attaining desired spatial resolution near the boundaries of spherical layers. Particle splitting is a typical way of improving local resolution [4]. For each particle in the region of higher resolution, an icosahedral template of particles is generated ( figure 7). Thus, we increase the number of particles by a factor of 12, linear resolution is increased by a factor of 2.3. Figure 7. An icosahedral template for particle splitting with an initial particle (red) and the added ones (blue). Now do calculation with adaptation by splitting for the spherical IC. Define 10 million particles in a 1/8 of a sphere and split the particles, which are near the interface. We will simulate 10 + 1 = 11 million particles.
After shock wave passing through the interface one can easily see that sphericity is broken (figure 8). Systematic perturbations are clearly seen. Let us try to remove them by rotating the  particle template arbitrarily at every splitting. The result is shown in figure 9. The perturbations are seen to be random but their amplitudes are still rather high.
Differences in calculated results can be quantified by the number of particles in spherical layers after loading (figure 10). The spheres adapted by splitting with no rotation reduce perturbations but their amplitudes are still rather high. Let us try to find some other solution.
Let us add spheres by increasing their radius not at a constant step, but at a smooth decrease in the radius increment towards the boundary (figure 12, left) and then at a smooth increase in it (figure 12, right). As seen from figure 11, most particles are added just on the boundary to describe it with desired resolution.    Figure 12. Adapted resolution: a section of particles in condensing spheres (left) and radius increment versus distance to boundary (right) with no adaptation (red) and with adaptation (black). Now do calculation with radial adaptation. Define 2 million particles in a half of a sphere. Adapt the radius by a factor of 5. We will simulate 2 + 21 = 23 million particles. Figure 13 shows DT ice after shock loading. The outer and inner boundaries show good sphericity. One can see small systematic circular perturbations perpendicular to the vertical axis. These are the artifacts of spheres construction: the vertical axis is the axis of the Fibonacci spiral, which gives us the spherical arrangement of particles.
Rotate arbitrarily each sphere and repeat calculation. It is seen from figure 14 that systematic perturbations turn random.
Evaluate interface sphericity in these two calculations. figure 15 shows that compared to the non-rotated spheres, the rotated ones induce more perturbations, entrapping more particles in them.

3D perturbations
All the above was done in order to make more stable the numerical simulation of real interface instabilities with predefined perturbations. The effect of initial conditions on instability evolution can be evaluated from calculations with and without perturbations on the inner boundary of the ablator. figure 16 shows the boundary between the ablator and DT ice with 3D perturbations.  Figure 15. Number of particles in spherical layers after loading.
Define 4 million particles in a 1/8 of a sphere and perform radial adaptation by a factor of 10 near the ablator -DT-ice interface. We will simulate 4 + 126 = 130 million particles. Apply 3D perturbations. figure 17 shows DT ice with and without initial perturbations after shock loading. The perturbations induced by imperfections in the IC and numerical scheme are seen to have much lower amplitudes than the initially applied perturbations.  So, we have demonstrated that in the simulation of spherical samples compiled of spherical layers which strongly differ in density, it is important to use initial conditions that do not induce artificial perturbations. The lattice IC with splitting, which are often used in plane problems, showed poor results when applied to spherical compression. We suggest that the spherical IC with radial adaptation should be used in this case to preserve local particle ordering and induce the least spurious perturbation under shock loading.