The Richtmyer-Meshkov instability of a double-layer interface in convergent geometry with magnetohydrodynamics

The Richtmyer-Meshkov instability (RMI) refers to the instability of an interface between two fluids that is impulsively accelerated, usually by a shock wave [1,2]. It occurs over wide ranges of length and time scales in technological applications and natural phenomena, such as supernova [3], combustion [4], hypersonic air breathing engines [5] and inertial confinement fusion (ICF), a promising approach for fusion energy generation [6]. In ICF, a small target filled with deuterium-tritium fuel mixture is heated by high-power lasers that drive an imploding shock into the target, compressing the fuel to a hotspot of sufficiently high temperature and pressure to initiate fusion reactions. In this high temperature and high energy-density scenario, the materials are expected to be in the plasma state, and thus may be influenced by a magnetic field. ICF experiments conducted on the Omega Laser showed that an external strong magnetic field might enhance the implosion performance by increasing the hotspot ion temperature and neutron yield [7,8]. A two dimensional radiation-hydrodynamics numerical investigation found that the temperature and pressure of hotspots for ignition decreased under the influence of a strong magnetic field; it also found that the field might suppress the growth of hydrodynamic (HD) instabilities [9], such as the RMI and Rayleigh-Taylor instability (RTI) [10,11]. The rapid growth of these instabilities on the target surface was responsible for the reduction of the energy production by breaking the spherical symmetry of the flow and severely degrading the final compression of the target [6].

In this investigation, we mainly focus on the evolution of RMI in converging flows under an external magnetic field within the framework of single-fluid magnetohydrodynamics (MHD), a fluid description of plasma dynamics. In planar geometry, the effects of an initial seed magnetic field on the RMI has been explored extensively. Samtaney [12] numerically studied the interaction between shock and inclined density interface in MHD, and found that the RMI was suppressed in the presence of a magnetic field. Wheatley et al. [13,14] investigated the case of the magnetic field perpendicular to the interface and showed that the mechanism of the suppression was attributed to the transport of the baroclinic voricity by MHD waves away from the density interface. For the case where the magnetic field was parallel [15] and oblique [16] to the interface, the RMI was also suppressed by the field. Samtaney [17] performed linear simulations of the RMI of double-interfaces in the presence of the magnetic field and showed that the growth rate of both interfaces decayed and oscillated around zero. In converging geometries, Bakhsh et al. [18,19] examined the evolution of instabilities via linear simulations in cylindrical geometry in the presence of normal and azimuthal fields and observed a transition from the early RMI phase to the RTI dominated phase (the distinction between RMI and RTI was most clearly seen in hydrodynamic cases). In addition, a significant suppression of the instability by a sufficiently strong magnetic field was observed. Mostert et al. [20] investigated the MHD RMI under the influence of two seed magnetic field configurations (uniform field and saddle field) in cylindrical and spherical converging flows. It showed that the extent of RMI suppression was not strongly dependent on field configuration, but that the saddle field resulted in a lower degree of implosion distortion or asymmetry. These prior nonlinear single fluid MHD investigations of the converging RMI culminated in a proposed octahedrally symmetric magnetic field configuration in 3D simulations [21]. The octahedrally symmetric field suppresses the instability comparably to the other previously considered seed field configurations for light-heavy interface accelerations while results in a higher degree of symmetry of the underlying flow even at high field strengths. These results reveal that the applied field of higher symmetry degree helps maintain the symmetry of imploding flow of single interface, which is very possibly for double interfaces.

In this work, we continue the thread of these previous investigations in 2D and hence exclude the investigation of octahedral 3D fields. In our work, motivated by the presence of multiple density interfaces in ICF, we numerically investigate the RMI when a converging cylindrical shock interacts with two interfaces separating fluids of three different densities (referred to as a double density layer) in the presence of a saddle-topology seed magnetic field. The nonlinear interactions present in the double density layer are significantly more complex than that for a single density interface, and there is nonlinear coupling between the two interfaces. We believe it is important to determine whether the presence of a seed magnetic field still effectively suppresses the interfacial growth and how the essential physical suppression mechanism, vorticity transport via MHD waves, is influenced by the double density layer. It's worth mentioning that Mikaelian [22,23] has linearly investigated RTI and RMI in stratified cylindrical and spherical concentric shells with incompressible hydrodynamic models. For double interface case, he noted that perturbations fed through from one interface to another when the shell was thin in the early stages of its evolution, while the coupling between interfaces decreased as the implosion thickened the shell, causing perturbations at each surface to grow independently. We observe the feedthrough phenomenon in the early stages when the shell is very thin, however we haven't observed the shell thickening by implosions. This difference may due to nonlinear effects which were not considered in Mikaelian's model. Wheatley et al. [13,15] developed analytical MHD incompressible models in Cartesian geometry for impulsively accelerated interfaces. Incompressible models for RMI/RTI in converging geometry with double interfaces in MHD entail several technical difficulties and are outside the main scope of this paper. The remainder of this paper is organized as follows: In Section 2, MHD equations and initial setup of the problem are introduced, along with a brief description of the numerical method. In Section 3, simulation results and their interpretation are presented. Conclusions are presented in Section 4.

The initial physical setup of the problem is shown in Fig. 1: two density interfaces (DI1 & DI2) centered at the origin with mean radii r₁ = 1.0 and r₂ = 0.5, respectively, are perturbed with a single-mode of azimuthal wavenumber k = 32, and an amplitude equal to 4% of their wavelength λ. The perturbed interface radii ζ_i for DI1 and DI2 are given as,

where we define ϕ = arctan y/x and consider the domain 0 < ϕ < π/2. Solution symmetry is imposed along ϕ = 0 and ϕ = π/2. A cylindrical driven shock results from a Riemann problem initialized across a Riemann interface (RI) located radially upstream at radius r_b = 1.5. The region r < r_b is separated into three sub-domains by the two DIs. The initially quiescent states in each sub-domain are initialized as p_I = p₀, ρ_I = 3ρ₀; p_II = p₀, ρ_II = 5ρ₀ (case A) or 2ρ₀ (case B) and p_III = p₀, ρ_III = ρ₀, respectively. Thus, for Case A, DI1 is a light-to-heavy interface (the incident shock propagates from the light fluid to the heavy) and DI2 is a heavy-to-light interface, while for Case B, both DI1 and DI2 are light-to-heavy interfaces. The state outside the RI is set as p_b = 12.1p₀, ρ_b = 3ρ₀. This Riemann problem leads to a primary converging shock of initial acoustic Mach number M = 2 at ϕ = 0. To investigate the influence of thickness d₁₂ of region II on the flow, r₂ is varied to be 0.5, 0.8 and 0.9, corresponding to the cases where thickness d₁₂ ≈ 2.5λ₁, λ₁ and 0.5λ₁.

A saddle-topology seed magnetic field is applied initially, expressed as [20],

with $αi={+α0,−α0,+α0,−α0}$ a signed scaling parameter that sets $|B(r)|=B0$ at r = r₁, and (x_i,y_i) = {(5,5),(−5,5), (−5,−5),(5,−5)}. The nondimensional parameter $β0=2p0/B02$ is used to describe the strength of the initial seed field. For both Case A and Case B, the effect of the field strength is discussed by varying β₀ to be ∞, 32 and 8, where β₀ = ∞ denotes the hydrodynamic (HD) case. For simplicity, all cases are described by the abbreviations outlined in Table 1.

Following several of the previous investigations (e.g. Refs. [15,21]), we employ the ideal MHD model for this study. The dimensionless variables are defined as:

where ρ, p, u and B are the density, pressure, velocity and magnetic field, respectively; and μ₀ is the permeability of free space. Neglecting the effect of body forces and dissipation, the nondimensionalized ideal MHD equations with above notations can be written as follows [24], with the carets omitted for simplicity,

where the specific heat ratio is fixed as γ = 5/3 throughout this study. In addition, we have divergence free constraint of the magnetic field, i.e.,

Because of the discretization errors, ∇⋅B can be non-zero and may increase with time leading to unphysical results [25]. Thus, this constraint should be numerically satisfied all the time during the simulation. A second-order nonlinear compressible finite volume code developed by Samtaney [26] is applied to solve the ideal MHD equations expressed in strong conservation form, using an unsplit upwinding scheme with a Roe flux solver. A projection method is used to enforce the divergence free constraint of the magnetic field [25]. In addition, our code has adaptive mesh refinement capability using the Chombo framework [27]. From symmetry considerations, the computation is performed on a quarter-domain, for 0 < x,y < 2 . All the simulations use a coarsest mesh at resolution 256² with two levels of refinement, with the refinement ratio of 4 in each direction for an effective resolution of 4096². The criterion for refinement, based on the local density gradient, is $|∇ρ|>0.02ρ$⁠. This mesh is sufficiently refined to resolve the perturbation amplitude according to the study by Mostert et al. [20]. Hence we omit details of other convergence tests here.

Before we discuss the RMI of the perturbed interfaces, it is instructive to examine the complex wave structures that arise from various nonlinear interactions. According to the initial setup, the incident shock (IS) generated from the RI propagates radially inwards and interacts with the two DIs. Fig. 2 shows various wave interactions with an x − t diagram wherein we plot contours of $log10(|∇ρ|+1)$ along line $ϕ=arctan12$ for both cases. It shows that contact discontinuity (CD) and expansion waves arising from the RI have no interaction with the DIs over the duration of the simulation. For HD (β = ∞) case (Fig. 2(a) and (c)), the IS converging processes can be summarized as follows: (1) the initially generated IS interacts with DI1 and the impulse lead to a radially inward motion of DI1; (2) the transmitted shock (TS) produced from step 1 interacts with the DI2, resulting in converging motion of DI2; (2) the reflected shock/rarefaction (shock for case B, rarefaction for case A) generated in step 2 propagates outward and interacts with DI1, while the reflected shock (RS) formed after the TS generated from step 2 reflects at the origin and travels radially outwards interacting first with DI2 and then with DI1, changing the direction of it's motion (this latter process is sometimes referred to as ”reshock”).

For the MHD (finite β) cases, the wave structures (Fig. 2(b) and (d)) are more complicated than the HD ones since more waves are formed during the primary MHD shock refraction at the interface, which in turn imply a higher number of shock-interface interactions. Initially, two incoming shocks are generated: the incident fast shock (IFS) and incident slow shock (ISS). For the IFS, the converging processes of MHD cases are similar to the HD ones except more waves are generated. For instance, during the IFS-DI1 interaction of case AL-8, a transmitted fast shock (TFS), a reflected fast shock (RFS), a transmitted slow shock (TSS) and a reflected slow shock (RSS) are produced, as seen more clearly in Fig. 3, which shows the vorticity field near DI1 for two cases with perturbed interfaces (these slow waves are too weak to be seen in the x − t diagram in Fig. 2). Unlike the IFS, the ISS has little influence on DI2 since the waves generated during the ISS-DI1 interaction are too weak to produce any discernible interactions, but it shows an obvious impact on the banded structure which arises from the perturbations on the interfaces. Since the density in region II of case AL is higher than that of case BL, the TS travels faster in this region of case BL and converges to the origin at an earlier time.

The density plots of case AL (t = 0.67) and case BL (t = 0.51) are shown in Fig. 4. At these times, the fast shocks have traveled across the two interfaces while the reflected waves generated from the TFS-DI2 interaction have not yet interacted with the converging DI1. A visual inspection of the density images make it apparent that the amplitude of the perturbations of the MHD case is smaller than that of the HD case, which suggests that the magnetic field does suppress the RMI, although not completely. Due to the lack of axisymmetry of the saddle field, the extent of the suppression varies on the interface, which leads to the non-axisymmetric growth of the perturbations. It appears that the perturbations are most diminished at $ϕ≈π/4$ where the field is nearly parallel to the interface and becomes the largest near region ϕ = 0 where the field is nearly perpendicular to the layer. We draw attention to the ”bulge” formed on the CD at ϕ = π/4 under the seed field, which results from the non-axisymmetry of the ISS. Under the effect of the seed field, the radial velocity of the ISS is minimum at ϕ = π/4, which leads to the formation of a curvature singularity on the ISS at this point resulting in two reflected shocks. The high pressure behind these two shocks pushes the CD out and forms the ”bulge” to balance the pressure [28]. A stronger field strength results in a more pronounced ”bulge”.

The suppression of the instability is attributed to two related effects:one is the transport of baroclinically generated vorticity, created at the interface during the shock-interface interaction, away from the interface by MHD waves. This is illustrated by the vorticity distributions shown in Fig. 5 of case AL-8 at different times. An examination of region $ϕ≠π/2$ shows that the baroclinic vorticity deposited on DI1 during shock interaction is transported away by the slow waves. The other suppression mechanism is observed at ϕ = π/4, where the vorticity remains in the vicinity of the interface, but according to Ref. [15], the initial distribution breaks up into waves traveling parallel and antiparallel to the interface. As these waves propagate, the induced velocity at a given interface location (and the perturbation growth rate) oscillates in time.

Fig. 5(a) shows the vorticity plot after the IFS-DI1 interaction at t = 0.30. At this point, the ISS has not yet interacted with DI1. The two fast shocks, TFS and RFS, from the IFS-DI1 interaction rapidly travel away from DI1, while the two slow shocks, TSS and RSS, are moving more closely tied to DI1, nearly co-moving with the interface at ϕ = π/4. The strength of ISS reaches its maximum at ϕ = π/4 and vanishes at ϕ = 0. The trajectory of ISS changes only slightly after the RFS passes and the waves produced during the interaction are too weak to discern in the plot. At t = 0.57, the ISS-DI1 interaction is in progress and the TFS-DI2 interaction has occurred (see Fig. 5(b)). Before interacting with DI1, the ISS is somewhat weakened by the RSS from the IFS-DI1 interaction. Unlike the IFS-DI1 interaction, the ISS-DI1 interaction is a long-term one. In the low ϕ region, the strength of ISS is low and does not interact strongly with DI1; while near ϕ = π/4, the ISS is sufficiently strong to distort DI1. In Fig. 5(c), the reflected TFS from the origin has propagated across DI2 and causes it to move radially outwards. At this time, the ISS has completed its interaction with DI1 and the reflected shocks behind ISS begin to influence DI1. Near ϕ = π/4, a concave-deformation of the interface forms. The reflected TFS keeps propagating outwards and interacts with all the slow waves from DI1. These slow waves are too weak to visibly affect the TFS. At t = 0.99, the reflected TFS has traveled across the DI1 and becomes distinctly non-axisymmetric after this interaction. Some refracted shocks generated by the perturbations appear behind the TFS. Along the ray ϕ = π/4, phase inversion of the perturbation amplitude occurs for both DI1 and DI2 interfaces. Also apparent at this time are the formation of additional transmitted and reflected slow shocks during the reshock of DI1 by the outgoing TFS. These perturbed slow shocks transport the reshock generated vorticity in the same manner described for the initial shock interaction. Thus the reshock induced RMI should be similarly mitigated.

Fig. 6 shows the vorticity plots of case BL-8 at different time. In this case, both DI1 and DI2 are light-to-heavy interfaces, thus, the behavior of DI2 is similar to that of DI1 if the influence of the ISS-DI1 interaction is not taken into account. At t = 0.51, the incoming fast waves have interacted with both DI1 and DI2. The slow shocks attached to each DI transport the vorticity away at angles where the field is not tangential to the interface. After the outgoing reflected fast shock (from the origin) has traversed over all DIs at t = 0.81, a second generation of slow waves travels with DIs together with slow waves from the previous interaction.

Fig. 7 plots the evolution of DI1 and DI2 of case A in (ϕ,r) space. For DI1, after being impacted by the incident shock, the RMI occurs, and the perturbations grow rapidly for the HD case while showing partial suppression by the magnetic field in the MHD case. It shows that the extent of suppression achieves its maximum in the region around $ϕ≈π/4$⁠, where the field is locally parallel to the interface, and becomes minimum around region $ϕ≈0$⁠. In Fig. 7(b), we note that the amplitude of perturbations decreases after reaching a certain value. This is due to the effect of the RTI stemming from the continuous acceleration of the interface as it converges. The effect of RTI is most evident near region ϕ = π/4, where the phase inversion happens around t = 0.62. Since DI2 is heavy to light, after its interaction with incident shock, the phase inversion occurs immediately due to the generated vorticity (with the opposite sign of vorticity generated to DI1), then the perturbations grow driven by RMI, as shown in Fig. 7(d).

To quantitatively study the effect of magnetic field on the growth of amplitude of perturbations, it is convenient to consider the amplitude of perturbations in different azimuthal sectors demarcated by different intervals of ϕ. Fig. 8 shows the evolution of the normalized amplitude $η/η0$ of DI1 and DI2 in two sectors for all cases. The results are compared with two reference amplitude curves: one is from a linear simulation of cylindrical RMI, the other is the initial growth rate calculated through Lombardini-Pullin (LP) model [29]. To aid the quantification of the effect of the magnetic field and dependence on the angle, we denote the sector where $ϕ∈[0,π/16]$ as ”low-ϕ” region and the section where $ϕ∈[3π/16,π/4]$ as ”high-ϕ” region. Labels ”a”, ”b”, ”c”, ”d” and ”e” denote specially chosen time instances for both HD cases, also approximately coinciding with times for MHD cases, in the high-ϕ region. Label ”a” represents the time when the IS interacts with DI1, ”b” is the time when the reflected wave from TS-DI2 interaction hits DI1, ”c” is the time when the reflected shock from origin interacts with DI1, ”d” is the time when TS-DI2 interaction happens, and ”e” represents the time when the reflected shock from origin interacts with DI2. In all cases, it shows that the amplitude of the perturbations is decreased by the magnetic field.

Because interface DI1 in case AL, and interfaces DI1 and DI2 of case BL are light-to-heavy interfaces, the evolutions of the perturbations appear similar. After time instance ”a” or ”d”, the amplitude of the perturbation grows rapidly to a peak followed by its decrease. This is attributed to the competition mechanism between RMI and RTI. The initial impulse due to the shock wave leads to the RMI but as the shock interface decelerates radially, the RTI manifests itself, driving growth in the opposite direction. In the same sector, a stronger field strength leads to a lower peak amplitude. For the same field strength, the peak amplitude value in the high ϕ region is less than that in the low ϕ region, since RMI is more highly suppressed where the field is close to parallel to the interface [20]. In addition, the oscillation of the amplitude in high ϕ region indicates the existence of phase inversion, which is observed in Fig. 7(b). The small change of the amplitude at ”b” comes from the interaction between DI1 and the reflected wave, which turns out to be the expansion wave for case AL and shock wave for case BL. Table 2 shows the maximum normalized amplitude of perturbations on DIs before time slot ”c” or ”e” in the low ϕ region. It shows that the extent of the suppression of case AL is smaller than that of case BL, especially for DI1.

The extent of asymmetry is an important factor considered in ICF. In order to investigate the effect of magnetic field on the DI asymmetry, the mean radial position of each DI, $r¯$ is measured and the asymmetry degree ζ(r) of interface is quantified by comparing amplitudes in the four sectors where $ϕ=[nπ/16,(n+1)π/16]$⁠, $n∈[0,3]$⁠. We divide each section into N (32 in this study) equivalent segments and calculate standard deviation of mean radial position of corresponding pieces in all sections. Thus, the asymmetry degree ζ(r) is the mean of these N standard deviations normalized by the initial perturbation amplitude. Fig. 9 compares ζ(r) for cases AL and BL with different magnetic field strengths before DI interaction with reflected shocks from origin (reshock). For the HD cases, the flow is symmetric, thus ζ(r) = 0 should be satisfied all the time. In our HD simulations, ζ(r) = 0 is only approximately true due to a mild asymmetry in the adaptive meshes. For MHD cases, the initial saddle magnetic field causes the IFS speed to vary depending on the local field conditions, breaking the axisymmetry. Thus during the IFS-DI interaction, ζ(r) increases rapidly due to the interaction of the asymmetric IFS with the initially symmetric interface, while the mean radial position of DI almost remains unchanged. Consequently, the stronger magnetic field is, the larger asymmetry degree becomes, as indicated in Fig. 9. The degree of asymmetry is large under strong field strength. In such situations, ISS-DI interaction becomes strong, resulting in a severe distortion of the DI, as can be seen in Fig. 5(c).

After shock interactions have occurred in the MHD cases, the perturbation amplitudes are different in each sector due to the asymmetric RMI suppression induced by the magnetic field. We introduce another measure of the perturbation amplitude $η¯$ defined as the root-mean-square of the perturbation amplitude in each sector. Fig. 10 compares the growth of relative root-mean-square amplitude of DIs for cases with different thicknesses. It shows that the magnetic field has a suppression effect on the perturbation growth for all cases. For DI1 of case A, the reflected rarefaction wave from DI2 accelerates the perturbation growth, enhancing the RMI. In Fig. 10(a), this is seen as the change in slope at $τ≈0.1,0.2,0.5$ for the small, medium and large layer thickness, respectively. Eventually the perturbation amplitude $η¯$ reaches a peak and then decreases due to the dominance of the RTI, which drives the perturbations towards a phase inversion. We note that the duration of the RMI phase is positively correlated to the layer thickness, i.e., as d₁₂ increases the peak amplitude and subsequent decrease in amplitude occurs later in time. The competition between the RTI and the RMI phases results in a non-monotonic behavior, i.e., the peak amplitude is reached by the layer of medium thickness. When d₁₂ is small, the amplitude growth reinforcement due to the reflected rarefaction driven RMI from DI2 is also small. Although the extent of this RMI reinforcement for AL cases is stronger than that of AM cases, the duration over which RMI dominates RTI is reduced (since the RMI reinforcement happens), and this eventually makes the maximum amplitude of AL cases smaller than that of AM cases (see Table 3).

For finite β, the RMI growth is suppressed due to the magnetic field and the suppression mechanisms discussed above. The effect of the layer thickness follows the same trend as in the HD cases, i.e., the medium layer thickness exhibits the largest growth. For DI2 of case A in HD, the shock interacts with a heavy-to-light interface and leads to a phase inversion of the amplitude. The second interface being closer to the origin experiences a stronger shock (due to convergence effects) for the larger layer thickness cases. In Fig. 10(b), it is seen that the medium and large layer thickness cases show a similar growth in DI2 perturbation amplitude until $τ≈0.28$⁠, while the small layer case exhibits substantially small amplitude. The largest layer width case experiences the highest early DI2 perturbation growth rate, but quickly experiences a reshock at $τ≈0.28$ (the reshock is seen as a sharp decrease in the amplitude in these plots). Consequently, the largest overall growth is experienced for the medium layer interface, as shown in Table 3. A similar behavior occurs for DI2 in the finite β case as seen in Fig. 10(b).

For interface DI1 of case B in HD, the reflected wave from DI2 is a shock that interacts with DI1 as if DI1 is a heavy-to-light interface, and weakens the RMI. Overall, the thickness d₁₂ shows a positive correlation with the weakened strength (see Fig. 10(c)). Although the extent of this RMI for BM cases is weaker than that for BL cases, the duration over which RMI dominates RTI is increased (since the RMI weakening happens), thus the maximum amplitude of BM cases is smaller than that of BL cases. For case BL-32 (finite β), since amplitude achieves maximum before the RMI dominance abates, the maximum amplitude of case BL-32 will be greater than that of case BM-32. For case B, DI2 is a light-to-heavy interface and upon interaction with the shock, shows a positive growth of the interface perturbations. As the layer width increases, the second interface is closer to the origin, and experiences RTI more rapidly for higher layer thickness cases (See Fig. 10(d)). The sharp change in amplitude is due to reshock from the shock reflection off the origin.

The isnteraction between a converging cylindrical shock and double interfaces in the presence of a saddle-topology seed magnetic field was numerically investigated. Two cases (light-heavy-light and light-heavy-heaviest) with various field strengths and layer thicknesses were computed and the results showed that the magnetic field did suppress the RMI, while the extent of the suppression varied with ϕ and led to azimuthal symmetry breaking of the DIs. The growth of the amplitude of perturbations was affected by both RMI and RTI. For light-to-heavy interfaces, the amplitude increased when RMI dominated over RTI and decreased when RTI dominated over RMI. In the region where the interface was nearly parallel to the magnetic field, the perturbations oscillated due to phase inversion. For heavy-to-light interfaces, amplitudes were increased subsequently after the RMI induced perturbation phase inversion. In addition, for a given field strength, the extent of suppression of case B (both interfaces were light-to-heavy) was larger than that of case A (light-to-heavy for the first interface and heavy-to-light for the second). The saddle field increased the asymmetry degree of the interface by distorting it through the ISS-DI interaction. The effect varying the layer thickness was also examined and the case of medium layer thickness generally experienced a larger growth compared with smaller or larger layer thickness cases due to the competing effects of stronger initial shock interactions, as the inner density interface was moved closer to the origin, reflected wave driven RMI, and RTI onset.

The authors have no conflict of interest.

This work was supported by the KAUST Office of Sponsored Research under Award No. URF/1/2162-01.