Delineation of the flow and mixing induced by Rayleigh-Taylor instability through tracers

by   Ge Zhang, et al.

Rayleigh-Taylor-instability(RTI) induced flow and mixing are of great importance in both nature and engineering scenarios. To capture the underpinning physics, tracers are introduced to make a supplement to discrete Boltzmann simulation of RTI in compressible flows. Via marking two types of tracers with different colors, the tracer distribution provides a clear boundary of two fluids during the RTI evolution. Fine structures of the flow and thermodynamic nonequilibrium behavior around the interface in a miscible two-fluid system are delineated. Distribution of tracers in its velocity phase space makes a charming pattern showing quite dense information on the flow behavior, which opens a new perspective for analyzing and accessing significantly deep insights into the flow system. RTI mixing is further investigated via tracer defined local mixedness. The appearance of Kelvin-Helmholtz instability is quantitatively captured by mixedness averaged align the direction of the pressure gradient. The role of compressibility and viscosity on mixing are investigated separately, both of which show two-stage effect. The underlying mechanism of the two-stage effect is interpreted as the development of large structures at the initial stage and the generation of small structures at the late stage. At the late stage, for a fixed time, a saturation phenomenon of viscosity is found that further increase of viscosity cannot see an evident decline in mixedness. The mixing statues of heavy and light fluids are not synchronous and the mixing of a RTI system is heterogenous. The results are helpful for understanding the mechanism of flow and mixing induced by RTI.



There are no comments yet.



Mixing Metaphors

Mixed metaphors have been neglected in recent metaphor research. This pa...

Geometrically reduced modelling of pulsatile flow in perivascular networks

Flow of cerebrospinal fluid in perivascular spaces is a key mechanism un...

Simulating dense granular flow using the μ(I)-rheology within a space-time framework

A space-time framework is applied to simulate dense granular flow. Two d...

Numerical Design of Distributive Mixing Elements

This paper presents a novel shape-optimization technique for the design ...

Gram Charlier and Edgeworth expansion for sample variance

In this paper, we derive a valid Edgeworth expansions for the Bessel cor...

The influence of initial perturbation power spectra on the growth of a turbulent mixing layer induced by Richtmyer-Meshkov instability

This paper investigates the influence of different broadband perturbatio...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Rayleigh-Taylor instability (RTI) occurs at the perturbed interface of two fluids when the pressure gradient points from the heavy to the light (Rayleigh, 1883; Taylor, 1950; Lewis, 1950). In inertial confinement fusion (ICF), which is expected to be a viable approach to fusion that relies on the inertia of fuel mass to provide confinement, control of the RTI development plays an important role on the successful ignition as pointed by Lindl et al. (2004); Edwards et al. (2013a, b); Wang et al. (2016). If it allows free development of RTI, the applied compression would be limited and in more serious cases, destroying the integrity of the shell and leading to self-ignition failure (Bradley, 2014). Additionally, RTI is also an important phenomenon in both scientific and engineering scenarios, like supernova evolution in astrophysics (Ribeyre et al., 2004), stratigraphic dynamics in geophysics (Lev & Hager, 2008), and material mixing (Vinningland et al., 2007). Therefore, it is necessary to understand and investigate the evolution of RTI both for fundamental science and engineering practices (Zhou, 2017a, b). The flow induced by RTI is time dependent, which undergoes linear and non-linear growing stages with fluids penetrating each other and mixing, and in some cases, it transits to fully turbulent flow. On different focuses, plentiful theoretical works have been made on the linear stage, weakly non-linear stage as well as the self-similar turbulence on the final stage (Goncharov, 2002; Mikaelian, 1998; Zhao et al., 2018; Clark & Zhou, 2003). However, it is rather hard to investigate the strong nonlinear and heterogeneous mixing processes induced by RTI analytically, so mechanisms of the nonlinear developing and fluids mixing processes are still subtle (Zhou et al., 2019). In this case, experiments and numerical simulations were utilized which provide direct and comprehensive images of RTI development.

Through experiments, the physical pattern of RTI was induced and firstly analyzed by Rayleigh (1883), Taylor (1950), Lewis (1950) and et al. The concept of RTI was then proposed and the coarse-grained description on the growth rate was derived (Kull, 1991a). In recent researches, different experimental apparatuses were designed to investigate RTI. Utilizing the experimental system consists of an accelerator to put the container of fluids down on lead rail with prescribed acceleration as well as the diagnostic instruments with planer laser-induced fluorescence (PLIF), Waddell et al. (2001) observed the development of RTI induced mixing layer in a miscible two-fluid system of low Atwood number, found that the initial growth rate agrees well with linear stability theory and the average of the late-time bubble and spike velocities is observed to be constant. Wilkinson & Jacobs (2007) designed similar experimental apparatus to create required acceleration and investigated the RTI triggered by a single-mode, three-dimensional perturbation. It was found that early time measured growth rate is consistent with linear stability theory, but the measured vertical interfacial velocity disagrees with many simplified theoretical models due to the formation of vortices is neglected in these models. Adopting a Z-fold schlieren system for visualizing, which illuminates the flow field by a mercury lamp light source and captures the flow by a high-speed video camera, Luo et al. (2018) uncovered an abating effect of RTI on Richtmyer-Meshkov instability (RMI). Although experimental results are delineated and persuasive, the drawback is that sometimes apparatuses are expensive and the results severely depends on the ability of flow-field capturing techniques such as PLIF and particle image velocimetry (PIV).

Therefore, numerical researches have also been conducted on RTI, and different methods, including level-set method (Chang et al., 1996), front tracking method (Glimm et al., 1986), volume-of-fluid method (Gerlach et al., 2006), smooth particle hydrodynamics (Shadloo et al., 2013), large-eddy simulation (Mellado et al., 2005), phase field method (Celani et al., 2009), direct numerical simulation (Mueschke & Schilling, 2009), and etc., were utilized for RTI simulation. With consideration of stabilizing and destabilizing effects on RTI, several factors, viscosity, compressibility, surface tension, and etc., were investigated in previous researches as well as delineated flow field information of RTI. With direct numerical simulation, Young et al. (2001) detailed pictures for both 2D and 3D development of the RTI induced mixing layer and distinguished three distinct evolutionary stages; Wei & Livescu (2012) investigated the effect of the initial perturbation shape on the growth of a two-dimensional single-mode RTI at low Atwood number. Liang et al. (2019) studied the multimode immiscible incompressible RTI under low Atwood number with an improved phase-field lattice Boltzmann method, and uncovered the effect of Reynolds number and initial condition (the perturbation wavelength and amplitude). Xie et al. (2017) optimized the initial mode of perturbation considering retarding effect in viscous RTI with numerical simulation based on volume of fluid method. Regarding the effect of viscosity, Hu et al. (2019) made a detailed numerical research in which they studied the growth rate of a 2D single-mode RT instability under different Reynolds numbers (viscosity) and formulated the relationship between the vorticity intensity, viscosity, and time at the reacceleration stage of RTI. Also, the effect of Kelvin-Helmholtz instability (KHI, growing from an unstable perturbed interface between two-fluids subject to a parallel shear flow) induced vortices on the later deceleration-acceleration stage was discussed in their work. However, as traditionally the simulation of RTI was based on macroscopic governing equations of fluid dynamics, mostly Euler and Navier-Stokes equations (NSE), and consequently processes on the fluent and complex nonequilibrium sides are little investigated. To reveal the interplay of various nonequilibrium behaviors in complex flow, the discrete Boltzmann method (DBM) has been developed and utilized as an effective physical modeling and simulation tool for nonequilibrium flows with various hydrodynamic instabilities (Xu et al., 2018; Gan et al., 2019; Lai et al., 2016; Chen et al., 2016, 2018; Lin et al., 2017; Zhang et al., 2019), including RTI, RMI, and KHI. Lai et al. (2016) investigated the effect of compressibility on RTI via DBM and found that the compressibility has a two-stage effect which delays RTI growing at initial stage and accelerates it at the later stage, meanwhile, thermodynamic nonequilibrium (TNE) characteristics can be utilized as physical criteria for tracking various interfaces. Based on a multi-relaxtion-time DBM, besides investigated viscosity, heat conductivity, and Prandtl number effects on RTI, the collaboration and competition between RMI and RTI, Chen et al. (2016, 2018) uncovered a series of high correlations between the nonuniformities of macroscopic quantity and nonequilibrium strength in the flow system. With a two-component DBM, Lin et al. (2017) studied the entropy of mixing in RTI evolution, and showed different stages of growth rate of mixing under low and high Reynolds numbers. Zhang et al. (2019) found that, via interfaces captured by appropriate TNE properties, the mixing of materials and that of internal energy resulted from hydrodynamic instability can be more conveniently probed. For example, the newly observed double-layer vortex shows the correlation and difference of density mixing and temperature mixing/unification.

Regarding a complex flow system induced by RTI, which though has been investigated for many years, a numerical research on it still faces dual challenges. The first challenge is how to simulate the internal structures and processes that we interested objectively. Noted that an accurate simulation is based on the assumption that the physical model is effective and has necessary functions on the corresponding spatial and time scale. Also, the error generated from physical modeling cannot be compensated for by pure attempts on the algorithm. Secondly, it is a challenge to extract valuable information from complex flow field. A simulation generally yields massive data containing extremely rich information. However, different approaches of analysis only delineate the system from its own perspective and physical scale. Consequently, how to design a method to extract more features, mechanisms and laws from a flow system remains a challenge.

Given the above considerations, this paper aims at a fine-grained physical delineation of RTI induced flow from two aspects, the modeling and the information capturing. More than conducting a DBM simulation for the compressible viscous flow in a two-miscible-fluids system with RTI features, one-way coupled tracer particles will also be introduced to provide a synchronous Lagrangian viewpoint on the flow field, which could also closely relate the simulation results to the experiments of RTI if PIV or similar flow imaging techniques is expected. Based on particle tracking, a new access for complex flow field analysis is desired, which would be introduced in this paper. Furthermore, the analysis based on tracer particles would bring some physical understanding of the RTI induced complex flow and mixing process.

Accordingly, the rest of this paper is organized as follows. The numerical modeling and method are presented in §2. In §3 we demonstrate particle tracking results as well as discussions on RTI flow and mixing. Specifically, effects of compressibility and viscosity are investigated in §3.4 and §3.5 separately. Finally, the conclusion is provided in §4.

2 DBM coupled with tracer particles

According to nonequilibrium statistical physics, state of the flow system is described by a distribution function (DF,

) or equivalently all its kinetic moments. In many practical cases, it is neither possible nor necessary to know all details of the DF. What one needs is to capture the main features of the DF, therefore delineating parts of flow behavior of the most concern. The main features of the DF can be captured/described by parts of its kinetic moments which are most relevant to the main features of flow behavior. At this side, a hydrodynamic model is a coarse-grained model of Boltzmann equation. For example, The NSE is for the case where the system deviates always slightly from its thermodynamic equilibrium state. The hydrodynamic quantities used to describe the flow system are just the initial several kinetic moments of the DF. It is clear that, with increasing of deviation from thermodynamic equilibrium state, more kinetic moments of the DF are needed to capture the main feature of the flow system. In the case without external force, from the Boltzmann equation to a DBM, there are three steps: (i) simplification of the collision term, (ii) discretization of the particle velocity space, (iii) characterize the nonequilibrium state and pick out meaningful nonequilibrium information. The first two steps belong to coarse-grained physical modeling. The third step is the core and main purpose of DBM. In the first two steps, the kinetic moments (or hydrodynamic quantities) used to describe the flow system must keep the same values before and after the simplification/discretization. In this way, the model is significantly simplified while the main features of described flow system are kept unchanged

(Gan et al., 2013).

The acceleration effect has been taken into account in DBM to simulate the RTI. As did in nearly all previous numerical researches, this study focus on the case where the nonequilibrium effects purely induced by the acceleration are negligible. Thus, the distribution function in the force term can be approximated by its corresponding equilibrium distribution function . Based on coarse-grained modeling of Boltzmann equation, including simplifications on collision term as Bhatnagar et al. (1954), force term as He et al. (1998), and discretization on the velocity space, a typical single-relaxation-time conservation equation of DBM is shown as below,


where is the discrete distribution function of fluid particles, r denotes the spatial position, is time, and indicates discrete velocity ( corresponds to different discrete velocities decided by velocity models). u is the fluid velocity, a is external body force, is the relaxation time, and is the equilibrium distribution function.

Through Chapman-Enskog expansion analysis, a discrete Boltzmann model with at least 16 different discrete velocity covers fully compressible NSE on the description of fluid dynamics (Xu et al., 2012, 2016; Gan et al., 2018). Thus, DBM simulation provides credible results of flow field comparing with numerical solvers of NSE, which has already been well tested in previous researches; on the other hand, DBM also provides a quantitative indicator of local thermodynamic nonequilibrium effects by defining nonequilibrium moments as (Xu et al., 2018; Gan et al., 2013; Lin et al., 2016),


where are the kinetic central moments, in which used the velocity as ; and

relates to the tensor order and specific component of the moment respectively.

Toward the delineation of flow and mixing from Lagrangian viewpoint, tracer particles will be introduced into the DBM simulation. Considering that the tracer particles should be small and light enough so that they have negligible influence on the flow field under investigation, one-way coupled particles are introduced into DBM simulation. The tracer particles are treated as points in the flow field and there is no collision between them (Dentz et al., 2016).

In order to start a RTI flow, initial interface of two fluids is described by a single sinusoidal function with known amplitude () and wavenumber (),


where wavenumber relates to wavelength as, .

Although the above single-mode RTI is a rather special case, in the long run, small irregularities will inevitably grow and break the repetitive symmetry of the initial conditions. Also, theoretically any perturbation is possible to be represented as a sum of different modes through Fourier expansion. Thus, single-mode RTI is still meaningful and the focus of many researches both experimentally and numerically.

In the simulation, hydrostatic initial state is guaranteed by an pre-set condition for the heavy fluid at upper side and the light fluid at bottom as follows (Scagliarini et al., 2010),

Thus, the pressure on the interface gained from each side is equal,


Atwood number indicates the degree of difference of two fluids which is defined as,


where the term on the right of second equal sign is valid for the case that ideal gas state equation is introduced. Following the above setting, the DBM code utilized for the following RTI simulation has been well tested in one of our previous research (Lai et al., 2016).

After introducing tracer particles, the Stokes number () indicates the working accuracy of tracer particles,


where is the relaxation time of a particle, is local flow velocity, is the representative length which is always selected as diameter of a particle. For , particles will detach from the flow especially where the flow decelerates abruptly; for , particles follow fluid streamlines closely. By adjusting to small enough, particles are capable for RTI tracking.

Following the trajectory of each particle, the description of each point particle is made on Lagrangian description. For the -th particle, its motion equation writes,


where is the spatial position of -th particle and indicates the particle velocity.

Assuming there is no delay on the momentum exchange between fluid and particle, the velocity of particles is directly determined by the local flow as,


where is a sharp Dirac function which would be implemented with a discrete form in the simulation. The velocity of the -th particle will change according to its current position and the local flow velocity, as from to during time period . To refresh the position of point particles, a 4th-order Runge-Kutta scheme is utilized for discretization of the particle motion equation according to Ramsden & Holloway (1991). Details of the algorithm are shown in appendix A.

Besides, in particle tracking experiments of fluid flow, the most desirable state is that tracer particles are less affected by factors such as inertia, diffusion, or Brownian motion, so that their movement fully conform to the streamline. Therefore, the assumption made on tracer particles is in line with ideal situation in experiments and the results can serve as a reference.

3 Results and discussion

A 2D simulation of compressible flow with RTI is conducted based on the above modeling and method. Combining physical resolution and computational cost, mesh size is chosen as and about 250 000 tracer particles are distributed randomly to the region initially and initial fluid parameters are selected as table 1. Many physical properties are changeable in a compressible fluid system both with time and space, including sound speed (), viscosity (dynamic or kinematic, or or ), heat conductivity (), Knudsen number (), and etc. Particles above and below the interface are marked differently (type-a above and type-b below) in order to distinguish mixing condition of the RTI system. Additionally, some particles marked as type-c are set at the interface on the purpose of monitoring deformation and physical changes.

Time step () 2.010
Grid size (=) 1.010
Relaxation time () 5.010
Acceleration () 1.0
Upper fluid temperature () 1.0
Interface pressure () 1.0
Atwood number () 0.6
Table 1: Initial physical parameters for RTI simulation

3.1 Fluid/particle evolution patterns of RTI

Previously, fluid density is always shown to indicate the evolution of RTI in literatures, which is undoubtedly clear in an incompressible immiscible fluids system. However, in widely existing cases of miscible and compressible fluids systems, the density pattern cannot provide a clear interface due to the formation of density transition layer (Gerya & Yuen, 2015; Kilkenny et al., 1994; Yang et al., 2002; Aziz & Chandra, 2000). Therefore, effective approaches are demanded for unequivocal visualization of flow field. Particle tracking technique is then utilized in this paper. Through identifying markers of tracer particles, particle tracking can provide a clear boundary of two fluids during RTI evolution, which is shown and compared with fluid density patterns in figure 1. At initial stages, from figures 1(a) to 1(d), fluid density patterns (the left half) and tracer particle patterns (the right half) are nearly identical, although the interface in the right half is sharper. The similarity of two patterns verifies the validity of particle tracking approach. However, the advantage of particle tracking is demonstrated in the late stage of RTI growing and mixing, from figures 1(e) to 1(h). When fluids are mixed, it is extremely hard to figure out which part the fluid comes from initially through identifying the density pattern. Nevertheless, distinct boundaries between two kinds of particles are provided even in figure 1(h) which is a rather turbulent case that drops are dissociate from the main body of fluids. Clearly shown by the particle distribution pattern in figure 1, fluids are mixed in the late stages and tiny vortex structures are generated. These distinguished tiny flow structures are important to the RTI mixing and will be analyzed in the following sections of this paper. Noted that, due to the discrete nature of particles, there are some unavoidable white-noisy points on the particle distribution pattern, but theses points do not interfere with the recognition of the characteristic structures if enough number of particles are provided. Previously, a two-component DBM was employed to simulate RTI by Lin et al. (2017), which also delineated the mixing status. However, two-component DBM model nearly doubles the computational cost of the simulation and it is inevitable and hard to determine the interaction of the two components physically, so when used to simulate miscible fluids with RTI characters, the method employed here is more convenient and efficient.

Figure 1: Comparing density and particle tracking manifestation of RTI growing and mixing in miscible fluids; the left half in each part: density contours (the boundary is indistinct at late time); the right half in each part: tracer particles patterns (with distinct boundary at late time). (a) ; (b) ; (c) ; (d) ; (e) ; (f) ; (g) ; (h) .

Considering that both the fingers of light and heavy fluid penetrate the intermedia position of the field during RTI growing, in figure 2 we extract the changes of macroscopic physical quantities (fluid density, temperature, and velocity) at over an initial time period, coarsely corresponding to the formation of fingers. It can be pointed from figures 2(a) and 2(b) that the density and temperature curves gradually shrink horizontally toward the center indicating that the finger of heavy fluid tapers and the finger of light fluid widens, so the two fingers are visually named spike (the heavy one) and bubble (the light one). As shown in figures 2(c) and 2(d), the fluid velocities of the left and the right half are antisymmetric on direction () and symmetric on direction (). At , the direction of reverses five times looking from left to right or backward. Inner the light fluid, the direction of is back to back, indicating that the bubble is blowing. At the interface, the velocity directs from the light to the heavy and shows the peak value horizontally, so the interface moves to the heavy part accordingly. Inner the heavy fluid, the velocity is front to front, indicating the heavy fluid shrinks. It is found that initially at the center of heavy fluid, there is a repulsion velocity to reject the shrinkage so as in the center of the light fluid, but the repulsion velocity diminishes quickly in heavy fluid (at ) and gradually diminishes in light fluid (at ). Thus, the direction of velocity of whole field is unified pointing from the light fluid to the heavy fluid. Additionally, it can be pointed that the movement of fluids starts at the interface of two fluids from figure 2(d) and there is a fluctuation of velocity horizontally at initial time. The shrinkage of spike and expansion of bubble lead to the figure of curve gradually takes a hopper shape at . Also, during RTI growing, is smoother inner the light fluid and has the peak value at the center, but there is a delay of in the center of heavy fluid especially forming a valley, which means the flow inner the spike and bubble are totally different. Additionally, the spike and the bubble not only differs on physical quantity distribution but on the growing velocity that the of spike is nearly 0.3 greatly larger than the of bubble which is about 0.2.

Figure 2: Macroscopic quantities at the position of during the fluid fingers formation time. (a) density (), (b) temperature (), (c) and (d) fluid velocity ().

During RTI growing and mixing, tracer particles transport according to local velocity and collectively indicate the state of RTI system from their perspective. Based on statistic of fully 250 000 particles, a particle velocity distribution (PVD) pattern is provided in figure 3. Initially at , when the RTI starts, most particles locate at the origin point (0, 0), indicating that the disturbance of RTI has not impact the system on large scale, but there are also some particles have already gained velocity and the PVD of them is nearly circular in shape, especially noted that a ring shape distribution pattern is formed in figure 3(a1) whose center slightly offsets from (0, 0). Particles on a ring shape PVD pattern may have equal kinetic energy if the ring center is at the origin point. However, due to the gravity induced overall movement of heavy fluid, the center of ring deviates from the origin. Eliminating the overall movement, the ring delineates a homogeneous feature of the instability evolution of RTI initially. At different stages of RTI evolution, varying characteristic patterns of PVD are distinguished. In summary, induced by RTI, groups of particles gained specific velocity directly lead to the formation of these characteristic patterns. Therefore, it provide us a perspective for RTI analysis and what we need is to describe these characteristic PVDs morphologically and try to link them with RTI. At , particles concentrates at several layers and onion-like structures are formed both in PVDs of type-a and type-b particles, see figures 3(b1) and 3(b2) respectively, but inner structures of two onions are different. As shown in figure 3(b1) visually, a drip develops from the origin and encircled by a crown. This PVD pattern indicates that small part of type-a particles which aligns with the spike break through the surrounded type-b particles during the RTI growing. At the same time, two mushroom structures are formed vertically inside an envelope see figure 3(b2). The area of mushroom is large compared with the whole pattern and points are sparsely spilled in it, indicating that the growing of bubble is on relative large scale compared with spike. As the histograms in figures 3(a1), 3(a2), 3(b1), and 3(b2) show, both the horizontal velocity () of type-a and type-b particles distribute normally at = 0.5 and 1.5, and although with slightly difference, the PVDs of vertical velocity (

) are close to normal distribution. At

, two petunia-like patterns are shown in figures 3(c1) and 3(c2). The outside profiles of two PVDs look like petals and with a pistil in the center which indicates that many particles have zero horizontal velocity and the vertical velocity of particles ranges on large scope. The type-a pistil is slenderer than type-b pistil. Therefore, the spike maintains its flow direction when it penetrates the light fluid while the bubble disperses when it is impeded in the heavy fluid. From the histograms at , the distribution of of both type-a and type-b particles centralize at (0, 0) and are close to normal distribution. Only a small swell exists at of type-a histogram; however, a distinct double-peak distribution is observed from histogram of type-b particles. At the time, KHI develops and lead to the vortex motion of both spike and bubble, which generates reversed motion of particles. The small swell of type-a histogram indicates that the KHI vortex is consisted by a small portion of heavy fluids and the double-peak of type-b histogram shows that the light fluid constitutes large portion of the KHI vortex. At , when there is no distinct spike and bubble of RTI system and the mixing dominates, spiral patterns are formed in two PVDs, see figures 3(d1) and 3(d2). Both the histogram of type-a and type-b extends so that more particles are with horizontal velocity. The histogram of type-a also extends but with a main peak at (0, 0). Each column of histogram of type-b gets close and the two-peak disappears, forming a near normal distribution pattern with high deviation. From the difference between histogram of type-a and type-b, we can conclude that the mixing statues of heavy and light fluids are not synchronous and the mixing of a RTI system is heterogenous. Although we have demonstrated some analysis and results from characteristic patterns of PVD, further investigation of pattern dynamics can still be made and is needed to uncover more physical understanding.

Figure 3: Velocity distribution of two kinds of tracer particles at four selected time, , responding to different stages of RTI development in figure 1, initiation of RTI, initiation of KHI, late RTI growing, and RTI mixing; (a1)-(d1) type-a particles in red color; (a2)-(b2) type-b particles in blue color. Histograms on the top and on the right of each figure indicate the distribution of and respectively.

3.2 Interfacial tracking manifestation

In order to delineate characteristics of interface during RTI growing, type-c particles are distributed initially at the interface. The deformation of interface provided by tracking method is presented in figure 4. Initially, a start-up slight perturbation exists at the interface. At , fingers of two fluids (spike and bubble) are formed at interface, and then penetrate each other. However, the developing rates of the two parts are not consistent. At , the spike crosses the position of while the bubble just reaches the position of . Then, the spike still grows faster than bubble, indicated by figures of late times as , 0.6, and 0.8. Due to acceleration, the interface continually deforms, while at , the velocity difference on shear side between the interface initiates observable KHI, with two large vortices take their shapes at each side of the spike, which are in according with classical observation (Kull, 1991a). The development of KHI is a critical step toward the opening of turbulent mixing (Peltier & Caulfield, 2003) and is also the reason for the evolution of the mushroom structures during the RTI (Zabusky, 1999). Due to the interface is stretched and elongated during KHI formation, limited number of particles (totally 128 in figure 4) distribute sparsely at late time, so it is not possible to sketch the profile of interface precisely as decreased number of particles are at front, see the profile at . After that, the particle tracking description for the interface is unacceptable. Therefore, it is important to design a tracking approach with appropriate number or concentration of tracer particles afore a particle tracking RTI experiment.

Figure 4: The deformation of interface during RTI growing illustrated by tracer particles. The fingering morphology of fluids is well captured at early stage but failed at late stage due to KHI development.

Furthermore, as the main place for hydro-instability initiation and development, it is worthwhile to pay special attention to physical feathers of the interface, which may also be helpful for controlling instability in ICF (Betti et al., 1998; Yu et al., 2018)

. Also, the design and improvement of all the interface tracking techniques need better physical understanding of the interface. Fundamentally, thermodynamic nonequilibrium (TNE) is the driving force of fluid system evolution. The complexity of non-equilibrium behavior determines the diversity of the description of non-equilibrium degree. Knudsen number, viscosity, heat transfer, the macroscopic gradients of density, flow rate, temperature and pressure are commonly used to characterize the degree of non-equilibrium. They all describe the degree of nonequilibrium of the system from their respective perspectives. But they are also highly condensed and averaged coarse-grained descriptions of certain information. Many specific information about non-equilibrium states, specific allocation of internal energy in different degrees of freedom, etc., through which they are invisible and cannot be directly studied. Therefore, in addition to the description mentioned above, we need a more detailed description on TNE. In this respect, the non-conservative dynamic moments of

can be used for a more detailed description of the nonequilibrium behavior, and the concept of non-equilibrium phase space to origin distance can be used to define the amount of coarsening non-equilibrium. These descriptions are interrelated and constitute a relatively complete description of the system’s non-equilibrium states and behaviors. The collaboration of TNE quantities given by DBM presents a rough but quantitative indicator for the specific nonequilibrium state of fluid system. Each independent component of the TNE quantities describes the nonequilibrium state from its own side. For example, indicates viscous stress tensor and indicates heat flux. and are of higher-order kinetic indicators, which do not have direct counterparts in traditional hydrodynamics, but can be seen as delineations for the transport behavior of the flow rate of momentum flow and energy flow. We can open a high-dimensional phase space by using the independent components of all the TNE quantities given by DBM. In this phase space, the origin point corresponds to the equilibrium state. Any one of other points in this phase space corresponds to a specific TNE state. The distance between this point and the origin point can be defined as a TNE strength, see equation 11. It is obvious that the nonequilibrium strength is a very rough description for the TNE state because all TNE states in the same sphere around the origin point have the same TNE strength (noted that the complexity of the nonequilibrium behavior of flow systems determines the multi-angle descriptions of the non-equilibrium states. Besides the nonequilibrium strength, Knudsen number, viscosity, heat conductivity, gradients of macroscopic quantities (density, flow rate, temperature, pressure, etc.) are commonly used to characterize the degree of non-equilibrium, and they all describe the degree of non-equilibrium of the system from their respective perspectives).


where indicates local TNE strength by measuring the distance between a specific point in phase space and the origin point. , , , and are different components of a TNE state.

In this case, we extract the TNE information of type-c particles and hence demonstrate the TNE property of interface. Mean TNE strength is then calculated as equation 12,


where is the length of interface and is the number of type-c particles.

In the following analysis of interfacial TNE, the period of = 0.0 to 1.2 is selected during which the application of type-c tracking particles is roughly acceptable. Figure 5 shows the evolution of average TNE strength on interface. The mean TNE strength has its maximum value at initial stage as , and then gradually decreases. At the start of the simulation, TNE triggers the initial instability of the system, and then the system evolves spontaneously to a new equilibrium state, so there is a peak of the mean TNE strength at the beginning and gradual decreases subsequently.

Figure 5: The mean TNE strength of interface over time.

We demonstrate the averaged velocity of interface and TNE quantities in figures 6 and 7 respectively. Due to the symmetry of the system, we only use the information on the left half interface and the mean values are obtained by arithmetically averaging the quantities on each type-c particles. In figure 6, the horizontal velocity of interface () rises near the beginning () and gradually drops until . Then it slightly increases and decreases in short time. The value of is greater than 0 from to 1.2, indicating that during the growing of spike, the interface moves toward the center wholly. On the vertical velocity (), it increases from 0 at the beginning to about -0.6 at , but stabilizes at following time period of to 1.0 and rises distinctly to about -2.0 at . The interface endures three stages of moving downward, the slow acceleration stage, the stabilization stage, and re-acceleration stage. From each TNE component in figures 7(a)-(d), we observe that most curves have a peak value at bout , which corresponds to the time when the interface reaches its maximum mean TNE strength. At time , the interface has not developed into any characteristic patterns, so before the system accumulates the driving potential to drive the system to evolve. In this case, fields of TNE quantities at are provided by figures 20 to 23 in appendix B to make a compensation for figure 7. Different TNE components, , , , and are all specific manifestation of a TNE state. Figure 7(a) shows the changes of over time, which corresponds to viscosity stress during the development of interface. At , and reach their maximum value respectively and then decrease to bottom and finally increase. The curves of and are nearly symmetrical along value of 0.0. During to 1.2, fluctuates but slowly increases, indicating that the momentum exchange on shear direction is slightly enhanced during interface deformation. In figure 7(b), four independent components of are shown. and indicate the flux of component of kinetic energy on and direction respectively and they deviate from equilibrium state in the opposite direction. Similarly, and indicate the flux of component of kinetic energy on and direction respectively and also they deviate from equilibrium state in the opposite direction. Good symmetry is found between and along value of 0.0. However, small deviation is shown between the inversed curve of and . The value of is greater than the inversed curve of especially at . The value of and both decrease over time, but and increase firstly and then decline at about . According to figure 7(c), reaches its peak value at and monotonously drops from to 1.2. The curve declines slightly in short period of time from and rises again and finally drops. Initially at the interface, the energy flux indicated by has larger value on vertical direction () than horizontal direction (), but the surpasses at about , delineating that the vertical energy flux is not dominant over the horizontal flux in the system. From dimensional analysis, corresponds to the flux of energy flux, which indicates a more complex physical process. It is found that two components of , and , reach their peak values at , and decline until and 1.1, respectively. Also, both the two components rise during the following time. Analysis is made out of our most present knowledge, but TNE phenomena is so fruitful that cannot be fully investigated, so we leave it for a future study. In sum, a specific non-equilibrium state requires a comprehensive consideration of various non-equilibrium components, and their collaboration is combined to give a more completed description of the flow system.

Figure 6: The mean velocity of left half interface; and are horizontal and vertical velocity respectively.
Figure 7: The mean TNE quantities of left half interface. Each TNE component indicates the physical state of interface differently. (a) ; (b) ; (c) ; (d) .

3.3 Delineation of RTI mixing

On the description of mixing state, different metrics are employed in previous researches. A coarse-grained indicator is the width of mixing layer, which is defined as the distance between the amplitude of the spike and the bubble (Cook & Zhou, 2002). Although many researches has been made on predicting the width of mixing layer in RTI (Cabot & Cook, 2006; Olson & Jacobs, 2009; Dimonte et al., 2004), their results and conclusions may differ due to different definition of the edge of mixing layer (Zhou & Cabot, 2019). Toward a fine-grained description of mixing, the actual mixed mass is utilized by Zhou et al. (2016), which is directly decided by the mass fractions of the two fluids. In experiment, with advanced high-energy-density-physics instrument, the measurements on mixed mass has already leading to some meaningful results by Ma et al. (2017). In a previous two-component DBM research made by Lin et al. (2017), the mixing process of two fluids is studied with an indicator, mixing entropy, but it also needs the local mass fraction of each fluid which is not possible to be achieved in a single-component DBM system. However, as an effective tool to delineate complex flows, the particle tracking technique provides us an approach for visualizing and analyzing local mixing state. In order to quantify the status of mixing, a local mixedness (LM) metric is defined based on the spatial information of particles as,


where , are local particle density (LPD) of type-a and type-b particles respectively. When there is no type-a or type-b particles in the statistical cell, equals to the minimum value 0.0, while if , equals to maximum value 1.0. Thus, ranges from 0.0 to 1.0, indicating a no mixing state to complete mixing state.

Noted that there is a key step of mapping discrete spatial information of particles to continuous field for obtaining the particle density. The statistical cell bridges the discrete information to the continuous field and it should be selected with the appropriate size in order to make an accurate demonstration of mixing. The statistical cell with the size of times of the grid step ( and ) is employed, and ranges from 1.0 to 8.0 is selected for a test.

Firstly, field averaged mixedness (AM) field is calculated according to equation 14,


where is the area of the whole field.

The growth of AM over time under different statistical cells is shown in figure 8. Initially, AM is nearly zero, while after a short period, the mixing accelerates and AM grows rapidly. Comparing results from different statistical cell, it can be pointed that with small cells ( for example), the AM value is relatively low during the mixing process, and with large cell ( for example), perturbation on AM is evident during the process. However, a smooth and distinct growing line of AM is achieved with a medium statistical size. In this case, we choose the cell size of for main demonstration in the following parts, but similar results with other cell sizes are also available which is presented in the appendix C.

Figure 8: Comparison of AM growing curves under different statistical cell sizes.

Figure 9 provides the details of AM based on . On the full range, AM grows exponentially with time, indicating that the mixing has a long low-speed growing period and then accelerates in late stages. Four representative time (, 1.5, 2.5, and 3.5) of LM are chosen and field patterns are shown in figures 9(b1) to 9(b4). It is shown that the mixing occurs at the contacting area of two fluids primly. As the interface elongates and distorts, the contacting area increases accordingly and AM rises consequently. The exponential relationship not only indicates the growth rate of AM but also uncovers the developing features of the RTI. Similarly, based on statistical cell size of and , different patterns of LM are shown by figure 24 in appendix C. According to figure 24, areas with LM value is too rare under statistical cell; therefore, we cannot observe any distinct information from the pattern. With large statistical cells of , many details are blurred due to coarse resolution; thus, the width of interface indicated by LM is too thick and not adequately accurate.

Figure 9: Delineation of mixing characteristic through AM and LM with statistical cell size of . (a) AM growing curve; (b1-b4) contours of LM at four selected time, , 1.5, 2.5, and 3.5 respectively. Exponential fit is made for AM and four selected time is marked on AM curve.

Integrating the LM along axis and then divide the integral length, the vertical mean of mixedness (VM) is,


where is the height of the simulation system.

The evolution pattern of VM based on statistical cell is shown in figure 10. Initially, VM is almost zero and fluctuates along axis. With time, the flow system evolves and two peaks of LM appear along axis. Comparing the appearing time and positions of two peaks and flow field visualization in figure 1, we conclude that the peaks arise with the evolution of KHI, which occurs at the tangent of the two fluids on both sides of RTI. The two peaks of the VM gradually increase along with time, signifying the enhancement of KHI over time. Finally, the VM is relatively homogeneous with no obvious peaks and discontinuities. At this point, the large-scale KHI has been evolved, and the mixing happens on smaller scales, which indicates that the flow system enters another mixing stage. Therefore, VM provides a good indicator for KHI occurrence, which can directly provide the time and horizontal position of KHI as well as roughly indicate the intensity of KHI thorough the value of VM. During capturing the KHI through the VM indicator, there is no need of empirical judgement or time consuming picturing, so it suits for many different cases especially three-dimensional cases whose flow field data is very large. Also, the VM processing could be done separately both during and post the simulation, which is another benefit for RTI research. Additionally, it can be seen that during the development of KHI, the VM on both sides is not completely symmetrical, due to the growing of small irregularities on interface. In the real world, complete symmetry of a system is not possible to achieve for the limitation of processing technology; at the same time, a disturbance during the evolution will inevitably cause the destruction of symmetry, so the asymmetry in the simulation can partially reflect the real RTI process.

Figure 10: VM profiles at different time. The area circled with dashed lines indicates the mixing area with KHI, which accompanies with the late growing stage of RTI and leads to the acceleration of mixing.

Furthermore, the time-dependent profiles of VM based on statistical cells with size of and are provided by figures 25(a) and 25(b) in appendix C. It is found that although the general patterns of VM are in accord with figure 10, details provided by these two figures are different. With a small statistical cell, the VM fluctuates violently and the distribution curve is very rough, while with a large statistical cell, the distribution curve is very smooth but lacks most of details. An intermediate-sized statistical cell can take advantages and avoid disadvantages of the both, providing a relatively comprehensive view of the VM information, better delineating the mixing process of RTI. Therefore, in the following parts of investigating some influencing factors of RTI mixing, results are given based on statistical size of .

3.4 The effect of compressibility

According to previous researches, compressibility of flow plays an important role on the growth rate of RTI, as well as the RTI mixing (Bernstein & Book, 1983; Kull, 1991b; Gauthier & Le Creurer, 2010). In this case, we intend to discuss the effect of compressibility on RTI mixing through the particle tracking approach and the afore-defined mixedness. According to the analysis in appendix D, the compression factor derives as Lai et al. (2016),


where is the wavenumber defined as ; is the wavelength of perturbation; is the sound speed. The sound speed at upper side of the interface is chosen to make a representation in the following analysis.

At the same time, in order to keep same the characteristic physical scale of different compressible systems, a necessary step is made to ensure the equality of Knudsen numbers (). According to the analysis in appendix D, the dimensionless factor number or is determined as,


By adjusting different and , it is possible to change the compression factor while keeping constant (in the current simulation, ). Six different characteristic compression factors are selected as 0.011596, 0.017393, 0.023191, 0.028989, 0.034787, and 0.040584. Noted that, due to changes on the temperature of fluids, sound speed adapts accordingly, so initial compression factors are selected as representation. Also, physical time is scaled based on the characteristic time and then the effect of compressibility on mixing is provided in figure 11. It is shown that AM curves of different compressibility are relatively close initially and the mixing accelerates later with the AM curves of strong compressibility obviously above other curves of relatively weak compressibility. Thus, we conclude that high compressibility may enhance the mixedness at late stage of RTI mixing. However, two selected periods of RTI mixing in figures 11(b) and 11(c) provide us with more details of the process. At the beginning (see figure 11(b)), there is an obvious decline of AM before under all compressibility conditions. As above patterns of particles already demonstrated, such as figure 9(b), the mixing occurs at the contacting area of two-kind particles. Therefore, a LPD pattern at is attached to make a detailed delineation. Clearly, there is a darker band near the interface, which means that the LPD is small, closely to 0. Therefore, in the early stage of instability development, there is a short-term pullback effect tending to smooth the perturbation and resulting in a break of LPD and a decline of AM. However, due to the effect of acceleration, the instability continues to develop. A decline-ascent pattern of AM is observed in figure 11. After this startup process, the AM of high compressibility becomes lower than that of low compressibility, which is opposite to the relationship between AM and compressibility at late stage. In the subsequent development of mixing, the AM curve of high compressibility catches up with the low compressibility AM curve at about shown by figure 11(c). Accordingly, we conclude that there are at least two stages in the mixing and note that the intersection of every two AM curves does not necessarily coincide. (I) In the first stage, when the fingers of two fluids insert into each other, high compressibility is not favorable for mixing; (II) after the instability develops for a certain period of time, the AM of high compressibility accelerates and ultimately exceed AM of low compressibility. As demonstrated in literatures, compressibility destabilizes RTI (Xue & Ye, 2010), so the AM at late state is certainly larger in a more compressible flow system; however, the detailed study above on AM shows a two-stage effect of compressibility which is in accord with our conclusion in a previous research (Lai et al., 2016).

Figure 11: The AM of fluid systems with different compressibility; where is the compression factor, the greater its value, the stronger the compressibility. Time is scaled as . (a) an overall profile of AM evolution; (b) a magnified view of AM with ranges from 0.0 to 2.5, which corresponds to the initial development of RTI mixing (the first stage); the attached LPD figure indicates the particle distribution at and the dark color band is area with low particle density; (c) a magnified view of AM with ranges from 3.5 to 5.5, in which an intersection of high and low compressibility curves exists (transition from the first to the second stage).

Detailed information of the RTI mixing under different compressibility is provided by figures 12 and 13. figure 12 shows the particle tracking patterns of fluid systems with different compressibility and figure 12 makes a compensation for mixing patterns from the distribution of local mixedness. Three representative times, , 4.0, and 7.0, are selected which correspond to the first stage, the transition and the late stage. Three compression factors, , 0.028989, and 0.040584, are chosen to make representation. The particle tracking patterns of RTI mixing are shown in figure 12. At , see figures 12(a1) to 12(a3), we can observe that the spike of low compressibility grows evidently faster and with larger deformation than that of high compressibility. The fast development of large structures in low compressibility fluid system accounts for its high AM at the first stage. At the transition time , see figures 12(b1) to 12(b3), regions of each kind of fluid in low compressibility system are morphologically integral and not highly expand compared with that of high compressibility system. The elongation and deformation of flow structures in high compressibility system lead to further generation of small structures, which contribute to the increase of mixedness. Finally, distinct difference is found between patterns ranging from low to high compressibility at . The state of high compressibility system is of highly turbulent with many dissociated drops and generated small KHI vortex structures; however, the pattern of low compressibility is regular and without evident dissociation of fluids. Therefore, the mixing status of high compressibility system fully transcends the low compressibility system at late stage. From the perspective of mixedness distribution, figure 13 shows similar results to figure 12 and clearly provides the difference on mixednees of the whole field. Based on the above analysis, we can attribute the two-stage effect of compressibility on RTI mixing to the development (penetration or expansion) of large structures and small structures at two stages respectively. Intuitively, a fluid system with low compressibility is “hard” and one with high compressibility is “soft”. At first stage, the “hard” one has relatively high penetrating ability and leads to high mixedness, while the “soft” one has low penetrating ability and leads to low mixedness. However, at late stage, more small structures are tend to generate in the “soft” system due to it is more delicate than the “hard” one, which then leads to high mixedness.

Figure 12: Particle tracking manifestation of RTI growing and mixing in different compressibility systems provided by two-type particles; (1)-(3) indicate different compressibility, , 0.028989, and 0.040584 respectively; (a)-(c) indicate different time of the system, , 4.0, and 7.0 respectively.

Figure 13: Local mixedness distribution in different compressibility systems; (1)-(3) indicate different compressibility, , 0.028989, and 0.040584 respectively; (a)-(c) indicate different time of the system, , 4.0, and 7.0 respectively.

Through fitting different curves in figure 11 with exponent formula expressed as: , where , , and are fitting parameters, each mode of mixing in different compressibility cases are extracted. In the fitting formula, is the characteristic time scale of the mixing process, and is a characteristic mixedness which also relates to the amplification of mixedness at a specific time. We make a comparison of and under different compressibility in figure 14. Exponential fits are made on the two parameters over compression factor , which both declines with . In this case, higher compressibility would lead to higher rate of mixing () but lower characteristic mixedness, which two together decide the value of AM. The underpinning physics of the two-stage effect of compressibility is obtained. Initially, in short of time, the exponent parts of all compressibility cases are small and AM is mainly decided by the characteristic mixedness; thus, compressibility decelerates the mixing. However, as time increases, the exponent part dominates, so compressibility enhances the mixing. Also, with this relationship, it is possible to normalize the AM feature of different compressibility.

Figure 14: Characteristic time scale () and characteristic mixedness () of AM differs with compressibility.

3.5 The effect of viscosity

The effect of viscosity is worth discussing due to the extreme high temperature in ICF would lead to significant changes on kinematic viscosity as pointed by Weber et al. (2014). Generally, viscosity is considered as a stabilizing effect on RTI at both linear and nonlinear growing stage (Sohn, 2009). During reacceleration stage of RTI, viscosity inhibits the development of KHI which then inhibits the RTI (Chen et al., 2018). Furthermore, after the reacceleration stage, the viscosity effect can lead to different complex behaviors of RTI (Hu et al., 2019). However, as both the development of RTI and KHI may affect the RTI mixing, mixedness at specific time can vary according to different viscosity. Therefore, this section intends to delineate the effect of viscosity on mixing with particle tracking approach. In the single-relaxation-time discrete Boltzmann model (SRT-DBM) with ideal gas state equation, the relaxation time and temperature decides the kinematic viscosity and the dynamic viscosity is then determined by relaxation time and fluid pressure ,


As pressure changes with time and space, the initial pressure () at interface is chosen for representation in calculating dynamic viscosity. In this case, different relaxation time are employed to change the kinematic/dynamic viscosity. Noted that, due to Prandtl number is fixed as 1 in SRT-DBM, the heat conductivity is also aligned with as,


Additionally, as an indicator of intrinsic physical characteristics of fluid, the relaxation time links with Knudsen number , which indicates the non-equilibrium state of the fluid, shown by equation 17.

The RTI mixing of different viscosity is simulated and the AM during the RTI evolution is delineated in figures 15(a) to 15(c). Observing from figure 15(a), all AM curves well-nigh overlap with each other during a long time from beginning, indicating that the effect of viscosity is not very prominent initially. However, at late time, curves separate with each other and the AM of the low-viscosity fluid is much higher than that of the high-viscosity fluid. The result aligns with the knowledge that viscosity acts as an inhibiting effect on RTI and the RTI mixing is also restrained in high viscosity fluids. However, an interesting phenomena is found at late stage of RTI mixing and shown by figure 15(b). AM profiles at three representative time, , 3.5, and 4.0, of late stage are chosen to demonstrate the relationship between the mixing and viscosity. We found that the AM decreases with viscosity but slowly drops with further increment of viscosity and tends to a constant value. Performing optimal exponential fitting of the curve, the AM has a saturation feature over viscosity, indicating that enough large viscosity can sufficiently suppress the RTI mixing. Therefore, in different fluid systems, if only with viscosity change, a saturation AM exists at different time. When magnifying AM curves at initial stage, a similar decline-ascent phenomenon of AM also exists which is in accord with results in the previous section. Besides, the AM of high-viscosity fluid is slightly higher than that of the low-viscosity fluid after the decline-ascent period. Therefore, we found an intersection between high viscosity curve and low viscosity curve form figure 15(c) with to 3.0. The result indicates that the viscosity is also a two-stage effect on RTI mixing. At initial stage, the viscosity slightly enhances the mixing; however, in the late stage, viscosity inhibits the mixing, and with the viscosity increase the AM tends to be constant.

Figure 15: The AM of fluid systems with different viscosity. The dynamic viscosity changes from to ; (a) an overall profile of AM evolution; (b) a cross-section profile of AM at , 3.5, 4.0 with different viscosity; the AM tends to be a constant value with viscosity increase; (c) a magnified view of AM with ranges from 1.3 to 3.0, in which the curves of high and low viscosity intersect (transition from the first to the second stage).

Both figures 16 and 17 are made for an illustration of viscosity effect on mixing at different states (four representative time, , 2.5, 3.0, and 3.5 are selected which correspond to the state at initial stage, the transition state of two stages, the further mixing state and more developed state at late stage respectively; also, results are based on three representative viscosities, , , and ). At , see figure 16(a), although the pattern of low viscosity has relatively elongated inner structure compared with the high viscosity system, but the AM of high viscosity system is greater than the low viscosity system. We can observe that the high viscosity system has an integral spike which directly penetrates into the light fluid of long distance, so this large scale development contributes to mixedness at first stage mostly. Also, slight difference is on the two high viscosity cases of and , which is in accord with the close AM of the two. However, significant differences are observed at both on the large structure growing and small structure developing, the low viscosity system fully transcends the high viscosity system, see figure 16(b). At the time, AM of low viscosity system also exceeds the AM of high viscosity system. At , shown by figure 16(c), as more tiny structures are generated in the low viscosity system, the AM of low viscosity system is significantly higher than that of high viscosity system. At , see figure 16(d), small flow structures break down from main large structure and drops dissociate from main fluid body in the low viscosity system, so the AM of low viscosity system constantly increases. By contrast, in high viscosity system, large structures are maintained and small structures are not developed. Also, little difference on the AM is found between the two high viscosity examples ( and ). The results of RTI mixing at late stage is in accord with conclusion made by Liang et al. (2019) that relatively high Reynolds number (corresponding to low viscosity in this paper) leads to high growing rate of RTI and generation of dissociated tiny fluid structures but no dissociation of drops in a low Reynolds number (high viscosity) system. However, detailed results presented above compensate for the understanding, which delineate a two-stage effect of viscosity on RTI mixing. We can also attribute the two-stage effect of viscosity to the development of large structures and small structures at two stages respectively. A fluid system with high viscosity is considered “hard” and one with low viscosity is considered “soft”. At first stage, the “hard” one has relatively high penetrating ability and leads to high mixedness, while the “soft” one has low penetrating ability and leads to low mixedness. In this case, the development of small structures at first stage caused by viscosity not prominently affects the mixing but the large scale penetration decides the mixedness. However, at late stage, more tiny structures are tend to generate in the “soft” system due to it is more delicate than the “hard” one; thus, generation of small structures leads to high mixedness. Additionally, a “hard” enough system can successfully suppress the generation of tiny structures, which causes a saturation effect of viscosity on mixedness.

Figure 16: Particle tracking manifestation of RTI growing and mixing in different viscosity systems provided by two-type particles; , , and . (a)-(d) indicate different time of the system, , 2.5, 3.0, and 3.5 respectively.

Figure 17: Local mixedness distribution in different viscosity systems; , , and . (a)-(d) indicate different time of the system, , 2.5, 3.0, and 3.5 respectively.

Furthermore, a scale analysis based on optimal exponential fitting of AM over viscosity is conducted. In figure 18, it is found that both the time scale () and characteristic mixedness () can be separated to two stages, which is totally different from the compressibility. Firstly, time scale increases with viscosity and then becomes nearly constant or with slight decrease with viscosity. The relationship between characteristic mixedness and viscosity is similar. In this case, the two stage effect of viscosity is more complex because it includes not only the transition of dominant role between the time scale and characteristic mixedness but also the self-transition of two parameters. At late stage, the overall effect of two factors which both not change distinctly decides the phenomenon that the AM saturates with viscosity increase.

Figure 18: Characteristic time scale () and characteristic mixedness () of AM differs with viscosity.

4 Conclusions

In this paper, tracer particles are introduced into DBM simulation of RTI in compressible flow. The tracer particles makes possible to observe finer structures with clear interfaces in the late mixing stage. One can clearly observe the mixing process and distinguish tracer particles from the heavy and light fluids. Various nonequilibrium behaviors in the RTI and material mixing processes are systematically investigated via the tracer particles. In particle velocity space the distribution of tracer particles show interesting patterns. These patterns bring quite dense information for the RTI evolution, which opens a new way for analysis and accesses significantly deep insights into the flow system. A mixedness is further defined via the tracer particles. The distribution of vertically averaged mixedness along the horizontal direction shows an interesting perspective to study the stage behavior in RTI evolution. The occurrence of double peak distribution can be roughly regarded as a criterion for the start of KHI.

Two physical factors, compressibility and viscosity, are explored regarding their effects on RTI mixing. Both of them show two-stage effects. Specifically, the compressibility initially inhibits and later enhances the RTI mixing, while the viscosity firstly enhances and then inhibits. The phenomena can be roughly understood via the following physical picture. At the initial period, the mixing is mainly induced by evolution of large structures, while at the late period the mixing is mainly induced by generation of small structures. Low compressibility and/or high viscosity make the system “hard” which favors the initial evolution of large structures. High compressibility and/or low viscosity make the system “soft” which favors the generation of small structures in the late stage.

The field averaged mixedness increases nearly exponentially with time, which indicates the existence of a characteristic time scale . The time scale show different dependence on compressibility and viscosity. It monotonically decreases with compressibility, while shows a two-stage dependence on viscosity. It firstly increases and then slightly decreases with viscosity, which means there is an appropriate viscosity for the system to get minimum mixedness. It is indicated that one can inhibit the mixedness via adjusting system viscosity. At the late stage, for a fixed time, the field averaged mixedness shows an exponential decrease with the viscosity, which indicates two points, (i) there exists a characteristic viscosity and (ii) the field averaged mixedness will not show meaningful decease with viscosity for high viscosities.

Many thanks to Dr. Yudong Zhang, Dr. Kaige Zhao, Prof. Lihao Zhao, Prof. Wen Song, Mr. Yiming Shan, Mr. Jiahui Song, and Mr. Dejia Zhang for their help. We acknowledge the support from the National Natural Science Foundation of China (Grant Nos. 11772064 and 11574390), CAEP Foundation (Grant No. CX2019033), Science Challenge Project (Grant No. JCKY2016212A501), foundation of computational physics, the opening project of State Key Laboratory of Explosion Science and Technology (Beijing Institute of Technology, Grant No. KFJJ19-01M) and Natural Science Foundation of Fujian Provinces (Grant No. 2018J01654).

Appendix A Implement of tracer particles

After the spatial discretization of fluid domain, point particles seldom locate at the grid node. Thus, the velocity of a point particle should be determined from the velocity of its nearby nodes (see figure 19), and the process is realized through adopting discrete Dirac function, denoted by which suits for discrete grids,


Figure 19: A schematic for the determination of particle velocity from grid nodes. The shade indicates the extent of velocity contribution made by each grid.

Two-dimensionally, function can be decomposed as the product of two parts,


The function utilized in this paper is,


In this way, the velocity information of the flow field is transferred to particles.

Appendix B Details of the Thermodynamic Noneequilibrium

At the time of , TNE properties of the fluid system are shown as follows, including four components, , , , and (see figures 20 to 23).

Figure 20: Field details of TNE component at . (a) ; (b) ; (c) .

Figure 21: Field details of TNE component at . (a) ; (b) ; (c) ; (d) .

Figure 22: Field details of TNE component at . (a) ; (b) .

Figure 23: Field details of TNE component at . (a) ; (b) ; (c) .

Appendix C Mixedness based on different statistical cells

The defined local mixedness depends on the size of statistical cell. Thus, the following figures demonstrate different results from two statistical cells, and 8.0 respectively (see figures 24 and 25).

Figure 24: Contours of local mixedness at four selected times obtained from two different statistical cells; (a) ; (b) ; (1)-(4) =0.5, 1.5, 2.5, and 3.5.
Figure 25: Profiles of vertically averaged mixedness distribution at different time (obtained from two different statistical cells). (a) ; (b) . The area circled by dashed lines indicates the occurrence of KHI.

Appendix D Deviation of Compression factor and Knudsen number

The Euler equation with ideal gas equation of state reads,


Choosing the unit of length as , the unit of time as , and the unit of velocity as , Euler equation is transferred to the following,


In the above equation, is the square of the ratio of the representative speed to the sound speed, defined as compression factor,


Following a similar procedure, BGK-Boltzmann equation can be transferred as,


where is the ratio of relaxation time to representative time, defined as Knudsen number,



  • Aziz & Chandra (2000) Aziz, Shiraz D. & Chandra, Sanjeev 2000 Impact, recoil and splashing of molten metal droplets. International Journal of Heat and Mass Transfer 43 (16), 2841–2857.
  • Bernstein & Book (1983) Bernstein, Ira B & Book, David L 1983 Effect of compressibility on the rayleigh ctaylor instability. Physics of Fluids 26 (2), 453–458.
  • Betti et al. (1998) Betti, R., Goncharov, V. N., Mccrory, R. L. & Verdon, C. P. 1998 Growth rates of the ablative rayleigh ctaylor instability in inertial confinement fusion. Physics of Plasmas 5 (5), 1446–1454.
  • Bhatnagar et al. (1954) Bhatnagar, P. L., Gross, E. P. & Krook, M. K. 1954 A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems. Physical Review 94 (3), 511–525.
  • Bradley (2014) Bradley, PA 2014 The effect of mix on capsule yields as a function of shell thickness and gas fill. Physics of Plasmas 21 (6), 062703.
  • Cabot & Cook (2006) Cabot, William H & Cook, Andrew W 2006 Reynolds number effects on rayleigh ctaylor instability with possible implications for type ia supernovae. Nature Physics 2 (8), 562.
  • Celani et al. (2009) Celani, Antonio, Mazzino, Andrea, Muratore-Ginanneschi, Paolo & Vozella, Lara 2009 Phase-field model for the rayleigh ctaylor instability of immiscible fluids. Journal of Fluid Mechanics 622, 115–134.
  • Chang et al. (1996) Chang, Yu-Chung, Hou, TY, Merriman, B & Osher, Stanley 1996 A level set formulation of eulerian interface capturing methods for incompressible fluid flows. Journal of Computational Physics 124 (2), 449–464.
  • Chen et al. (2016) Chen, Feng, Xu, Aiguo & Zhang, Guangcai 2016 Viscosity, heat conductivity, and prandtl number effects in the rayleigh ctaylor instability. Frontiers of Physics 11 (6), 114703.
  • Chen et al. (2018) Chen, Feng, Xu, Aiguo & Zhang, Guangcai 2018 Collaboration and competition between richtmyer-meshkov instability and rayleigh-taylor instability. Physics of Fluids 30 (10), 102105.
  • Clark & Zhou (2003) Clark, Timothy T & Zhou, Ye 2003 Self-similarity of two flows induced by instabilities. Physical Review E 68 (6), 066305.
  • Cook & Zhou (2002) Cook, Andrew W & Zhou, Ye 2002 Energy transfer in rayleigh-taylor instability. Physical Review E 66 (2), 026312.
  • Dentz et al. (2016) Dentz, Marco, Kang, Peter K, Comolli, Alessandro, Le Borgne, Tanguy & Lester, Daniel R 2016 Continuous time random walks for the evolution of lagrangian velocities. Physical Review Fluids 1 (7), 074004.
  • Dimonte et al. (2004) Dimonte, Guy, Youngs, DL, Dimits, A, Weber, S, Marinak, M, Wunsch, S, Garasi, C, Robinson, A, Andrews, MJ & Ramaprabhu, P 2004 A comparative study of the turbulent rayleigh ctaylor instability using high-resolution three-dimensional numerical simulations: the alpha-group collaboration. Physics of Fluids 16 (5), 1668–1693.
  • Edwards et al. (2013a) Edwards, MJ, Patel, PK, Lindl, JD, Atherton, LJ, Glenzer, SH, Haan, SW, Kilkenny, JD, Landen, OL, Moses, EI & Nikroo, A 2013a Progress towards ignition on the national ignition facility. Physics of Plasmas 20 (7), 070501.
  • Edwards et al. (2013b) Edwards, MJ, Patel, PK, Lindl, JD, Atherton, LJ, Glenzer, SH, Haan, SW, Kilkenny, JD, Landen, OL, Moses, EI, Nikroo, A & others 2013b Progress towards ignition on the national ignition facility. Physics of Plasmas 20 (7), 070501.
  • Gan et al. (2013) Gan, Yanbiao, Xu, Aiguo, Zhang, Guangcai & Yang, Yang 2013 Lattice bgk kinetic model for high-speed compressible flows: Hydrodynamic and nonequilibrium behaviors. EPL (Europhysics Letters) 103 (2), 24003.
  • Gan et al. (2018) Gan, Yanbiao, Xu, Aiguo, Zhang, Guangcai, Zhang, Yudong & Succi, Sauro 2018 Discrete boltzmann trans-scale modeling of high-speed compressible flows. Physical Review E 97 (5), 053312.
  • Gan et al. (2019) Gan, Yan Biao, Xu, Ai Guo, Zhang, Guang Cai, Lin, Chuan Dong, Lai, Hui Lin & Liu, Zhi Peng 2019 Nonequilibrium and morphological characterizations of kelvin chelmholtz instability in compressible flows. Frontiers of Physics 14 (4).
  • Gauthier & Le Creurer (2010) Gauthier, S & Le Creurer, B 2010 Compressibility effects in rayleigh ctaylor instability-induced flows. Philosophical Transactions of The Royal Society A: Mathematical, Physical and Engineering Sciences 368 (1916), 1681–1704.
  • Gerlach et al. (2006) Gerlach, D, Tomar, G, Biswas, G & Durst, F 2006 Comparison of volume-of-fluid methods for surface tension-dominant two-phase flows. International Journal of Heat and Mass Transfer 49 (3-4), 740–754.
  • Gerya & Yuen (2015) Gerya, Taras V & Yuen, David A 2015 Rayleigh ctaylor instabilities from hydration and melting propel cold plumes at subduction zones. Earth and Planetary Science Letters 212 (1), 47–62.
  • Glimm et al. (1986) Glimm, J, McBryan, O, Menikoff, R & Sharp, DH 1986 Front tracking applied to rayleigh ctaylor instability. SIAM Journal on Scientific and Statistical Computing 7 (1), 230–251.
  • Goncharov (2002) Goncharov, VN 2002 Analytical model of nonlinear, single-mode, classical rayleigh-taylor instability at arbitrary atwood numbers. Physical Review Letters 88 (13), 134502.
  • He et al. (1998) He, Xiaoyi, Shan, Xiaowen & Doolen, Gary D. 1998 Discrete boltzmann equation model for nonideal gases. Physical Review E 57 (1), R13.
  • Hu et al. (2019) Hu, Ze-Xi, Zhang, You-Sheng, Tian, Baolin, He, Zhiwei & Li, Li 2019 Effect of viscosity on two-dimensional single-mode rayleigh-taylor instability during and after the reacceleration stage. Physics of Fluids 31 (10), 104108.
  • Kilkenny et al. (1994) Kilkenny, J. D., Glendinning, S. G., Haan, S. W., Hammel, B. A., Lindl, J. D., Munro, D., Remington, B. A., Weber, S. V., Knauer, J. P. & Verdon, C. P. 1994 A review of the ablative stabilization of the rayleigh-taylor instability in regimes relevant to icf. Physics of Plasmas 1 (5), 1379.
  • Kull (1991a) Kull, H. J. 1991a Theory of the rayleigh-taylor instability. Physics Reports 206 (5), 197–325.
  • Kull (1991b) Kull, H. J. 1991b Theory of the rayleigh-taylor instability. Physics Reports 206 (5), 197–325.
  • Lai et al. (2016) Lai, Huilin, Xu, Aiguo, Zhang, Guangcai, Gan, Yanbiao, Ying, Yangjun & Succi, Sauro 2016 Nonequilibrium thermohydrodynamic effects on the rayleigh-taylor instability in compressible flows. Physical Review E 94 (2), 023106.
  • Lev & Hager (2008) Lev, Einat & Hager, Bradford H 2008 Rayleigh ctaylor instabilities with anisotropic lithospheric viscosity. Geophysical Journal International 173 (3), 806–814.
  • Lewis (1950) Lewis, D J 1950 The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. ii. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 202 (1068), 81–96.
  • Liang et al. (2019) Liang, Hong, Hu, Xiaoliang, Huang, Xuefeng & Xu, Jiangrong 2019 Direct numerical simulations of multi-mode immiscible rayleigh-taylor instability with high reynolds numbers. Physics of Fluids 31 (11), 112104.
  • Lin et al. (2016) Lin, Chuandong, Xu, Aiguo, Zhang, Guangcai & Li, Yingjun 2016 Double-distribution-function discrete boltzmann model for combustion. Combustion and Flame 164, 137–151.
  • Lin et al. (2017) Lin, C., Xu, A., Zhang, G., Luo, K. H. & Li, Y. 2017 Discrete boltzmann modeling of rayleigh-taylor instability in two-component compressible flows. Physical Review E 96 (5), 053305.
  • Lindl et al. (2004) Lindl, John D, Amendt, Peter, Berger, Richard L, Glendinning, S Gail, Glenzer, Siegfried H, Haan, Steven W, Kauffman, Robert L, Landen, Otto L & Suter, Laurence J 2004 The physics basis for ignition using indirect-drive targets on the national ignition facility. Physics of Plasmas 11 (2), 339–491.
  • Luo et al. (2018) Luo, Xisheng, Zhang, Fu, Ding, Juchun, Si, Ting, Yang, Jiming, Zhai, Zhigang & Wen, Chih Yung 2018 Long-term effect of rayleigh ctaylor stabilization on converging richtmyer cmeshkov instability. Journal of Fluid Mechanics 849, 231–244.
  • Ma et al. (2017) Ma, T, Patel, PK, Izumi, N, Springer, PT, Key, MH, Atherton, LJ, Barrios, MA, Benedetti, LR, Bionta, R & Bond, E 2017 The role of hot spot mix in the low-foot and high-foot implosions on the nif. Physics of Plasmas 24 (5), 056311.
  • Mellado et al. (2005) Mellado, Juan Pedro, Sarkar, Sutanu & Zhou, Ye 2005 Large-eddy simulation of rayleigh-taylor turbulence with compressible miscible fluids. Physics of Fluids 17 (7), 076101.
  • Mikaelian (1998) Mikaelian, Karnig O 1998 Analytic approach to nonlinear rayleigh-taylor and richtmyer-meshkov instabilities. Physical Review Letters 80 (3), 508.
  • Mueschke & Schilling (2009) Mueschke, Nicholas J & Schilling, Oleg 2009 Investigation of rayleigh ctaylor turbulence and mixing using direct numerical simulation with experimentally measured initial conditions. i. comparison to experimental data. Physics of Fluids 21 (1), 014106.
  • Olson & Jacobs (2009) Olson, DH & Jacobs, JW 2009 Experimental study of rayleigh ctaylor instability with a complex initial perturbation. Physics of Fluids 21 (3), 034103.
  • Peltier & Caulfield (2003) Peltier, WR & Caulfield, CP 2003 Mixing efficiency in stratified shear flows. Annual review of fluid mechanics 35 (1), 135–167.
  • Ramsden & Holloway (1991) Ramsden, Dave & Holloway, Greg 1991 Timestepping lagrangian particles in two dimensional eulerian flow fields. Journal of Computational Physics 95 (1), 101–116.
  • Rayleigh (1883) Rayleigh 1883 Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density. Proceedings of the London Mathematical Society 14 (1), 170–177.
  • Ribeyre et al. (2004) Ribeyre, X, Tikhonchuk, VT & Bouquet, S 2004 Compressible rayleigh ctaylor instabilities in supernova remnants. Physics of Fluids 16 (12), 4661–4670.
  • Scagliarini et al. (2010) Scagliarini, Andrea, Biferale, Luca, Sbragaglia, Mauro, Sugiyama, Kazuyasu & Toschi, Federico 2010 Lattice boltzmann methods for thermal flows: Continuum limit and applications to compressible rayleigh ctaylor systems. Physics of Fluids 22 (5), 055101.
  • Shadloo et al. (2013) Shadloo, MS, Zainali, Amir & Yildiz, M 2013 Simulation of single mode rayleigh ctaylor instability by sph method. Computational Mechanics 51 (5), 699–715.
  • Sohn (2009) Sohn, Sung-Ik 2009 Effects of surface tension and viscosity on the growth rates of rayleigh-taylor and richtmyer-meshkov instabilities. Physical Review E 80 (5), 055302.
  • Taylor (1950) Taylor, Geoffrey Ingram 1950 The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. i. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 201 (1065), 192–196.
  • Vinningland et al. (2007) Vinningland, Jan Ludvig, Johnsen, Øistein, Flekkøy, Eirik G, Toussaint, Renaud & Måløy, Knut Jørgen 2007 Granular rayleigh-taylor instability: Experiments and simulations. Physical review letters 99 (4), 048001.
  • Waddell et al. (2001) Waddell, JT, Niederhaus, CE & Jacobs, Jeffrey W 2001 Experimental study of rayleigh ctaylor instability: low atwood number liquid systems with single-mode initial perturbations. Physics of Fluids 13 (5), 1263–1273.
  • Wang et al. (2016) Wang, LF, Ye, WH, Wu, JF, Liu, Jie, Zhang, WY & He, XT 2016 A scheme for reducing deceleration-phase rayleigh ctaylor growth in inertial confinement fusion implosions. Physics of Plasmas 23 (5), 052713.
  • Weber et al. (2014) Weber, CR, Clark, DS, Cook, AW, Busby, LE & Robey, HF 2014 Inhibition of turbulence in inertial-confinement-fusion hot spots by viscous dissipation. Physical Review E 89 (5), 053106.
  • Wei & Livescu (2012) Wei, Tie & Livescu, Daniel 2012 Late-time quadratic growth in single-mode rayleigh-taylor instability. Physical Review E 86 (4), 046405.
  • Wilkinson & Jacobs (2007) Wilkinson, JP & Jacobs, Jeffrey W 2007 Experimental study of the single-mode three-dimensional rayleigh-taylor instability. Physics of fluids 19 (12), 124102.
  • Xie et al. (2017) Xie, C. Y., Tao, J. J., Sun, Z. L. & Li, J. 2017 Retarding viscous rayleigh-taylor mixing by an optimized additional mode. Physical Review E 95 (2), 023109.
  • Xu et al. (2012) Xu, Aiguo, Zhang, Guangcai, Gan, Yanbiao, Chen, Feng & Yu, Xijun 2012 Lattice boltzmann modeling and simulation of compressible flows. Frontiers of Physics 7 (5), 582–600.
  • Xu et al. (2016) Xu, AiGuo, Zhang, GuangCai, Ying, YangJun & Wang, Cheng 2016 Complex fields in heterogeneous materials under shock: modeling, simulation and analysis. SCIENCE CHINA Physics, Mechanics & Astronomy 59 (5), 650501.
  • Xu et al. (2018) Xu, Aiguo, Zhang, Guangcai & Zhang, Yudong 2018 Discrete Boltzmann modeling of compressible flows, pp. 450–458. IntecOpen.
  • Xue & Ye (2010) Xue, Chuang & Ye, Wenhua 2010 Destabilizing effect of compressibility on rayleigh ctaylor instability for fluids with fixed density profile. Physics of Plasmas 17 (4), 042705.
  • Yang et al. (2002) Yang, J., D Onofrio, A., Kalliadasis, S. & Wit, A. De 2002 Rayleigh ctaylor instability of reaction-diffusion acidity fronts. Journal of Chemical Physics 117 (20), 9395–9408.
  • Young et al. (2001) Young, Y-N, Tufo, H, Dubey, A & Rosner, R 2001 On the miscible rayleigh ctaylor instability: two and three dimensions. Journal of Fluid Mechanics 447, 377–408.
  • Yu et al. (2018) Yu, C. X., Xue, C., Liu, J., Hu, X. Y., Liu, Y. Y., Ye, W. H., Wang, L. F., Wu, J. F. & Fan, Z. F. 2018 Multiple eigenmodes of the rayleigh-taylor instability observed for a fluid interface with smoothly varying density. Physical Review E 97 (1), 013102.
  • Zabusky (1999) Zabusky, Norman J 1999 Vortex paradigm for accelerated inhomogeneous flows: Visiometrics for the rayleigh-taylor and richtmyer-meshkov environments. Annual review of fluid mechanics 31 (1), 495–536.
  • Zhang et al. (2019) Zhang, Yudong, Xu, Aiguo, Zhang, Guangcai, Chen, Zhihua & Pei, Wang 2019 Discrete boltzmann method for non-equilibrium flows: based on shakhov model. Computer Physics Communications 238, 50–65.
  • Zhao et al. (2018) Zhao, KG, Wang, LF, Xue, C, Ye, WH, Wu, JF, Ding, YK & Zhang, WY 2018 Thin layer model for nonlinear evolution of the rayleigh-taylor instability. Physics of Plasmas 25 (3), 032708.
  • Zhou (2017a) Zhou, Ye 2017a Rayleigh ctaylor and richtmyer cmeshkov instability induced flow, turbulence, and mixing. i. Physics Reports 720-722, 1–136.
  • Zhou (2017b) Zhou, Ye 2017b Rayleigh ctaylor and richtmyer cmeshkov instability induced flow, turbulence, and mixing. ii. Physics Reports 723-725, 1–160.
  • Zhou & Cabot (2019) Zhou, Ye & Cabot, William H 2019 Time-dependent study of anisotropy in rayleigh-taylor instability induced turbulent flows with a variety of density ratios. Physics of Fluids 31 (8), 084106.
  • Zhou et al. (2016) Zhou, Ye, Cabot, William H & Thornber, Ben 2016 Asymptotic behavior of the mixed mass in rayleigh ctaylor and richtmyer cmeshkov instability induced flows. Physics of Plasmas 23 (5), 052712.
  • Zhou et al. (2019) Zhou, Ye, Clark, Timothy T, Clark, Daniel S, Gail Glendinning, S, Aaron Skinner, M, Huntington, Channing M, Hurricane, Omar A, Dimits, Andris M & Remington, Bruce A 2019 Turbulent mixing and transition criteria of flows induced by hydrodynamic instabilities. Physics of Plasmas 26 (8), 080901.