We investigate the shock-induced turbulent mixing between a light and heavy gas, where a Richtmyer-Meshkov instability (RMI) is initiated by a shock wave. The prescribed initial conditions define a deterministic multimode interface perturbation between the gases, which can be imposed exactly for different simulation codes and resolutions to allow for quantitative comparison. Well-resolved large-eddy simulations are performed using two different and independently developed numerical methods with the objective of assessing turbulence structures, prediction uncertainties and convergence behaviour. The two numerical methods differ fundamentally with respect to the employed subgrid-scale regularisation, each representing state-of-the-art approaches to RMI. Unlike previous studies the focus of the present investigation is to quantify uncertainties introduced by the numerical method, as there is strong evidence that subgrid-scale regularisation and truncation errors may have a significant effect on the linear and non-linear stages of the RMI evolution. Fourier diagnostics reveal that the larger energy containing scales converge rapidly with increasing mesh resolution and thus are in excellent agreement for the two numerical methods. Spectra of gradient dependent quantities, such as enstrophy and scalar dissipation rate, show stronger dependencies on the small-scale flow field structures in consequence of truncation error effects, which for one numerical method are dominantly dissipative and for the other dominantly dispersive. Additionally, the study reveals details of various stages of RMI, as the flow transitions from large scale, non-linear entrainment, to fully developed turbulent mixing. The growth rates of the mixing zone widths as obtained by the two numerical methods are before re-shock, and long after re-shock. The decay rate of turbulence kinetic energy is consistently at late times, where the molecular mixing fraction approaches an asymptotic limit . The anisotropy measure approaches an asymptotic limit of implying that no full recovery of isotropy within the mixing zone is obtained, even after re-shock. Spectra of density, turbulence kinetic energy, scalar dissipation rate and enstrophy are presented and show excellent agreement for the resolved scales. Probability density function of the heavy-gas mass fraction and vorticity reveal that the light-heavy gas composition within the mixing zone is accurately predicted, whereas it is more difficult to capture the long-term behaviour of the vorticity.
On the Richtmyer-Meshkov instability]On the Richtmyer-Meshkov instability evolving from a deterministic multimode planar interface V. K. Tritschler, B. J. Olson, S. K. Lele, S. Hickel, X. Y. Hu and N. A. Adams]V. K. Tritschler††thanks: Email address for correspondence: email@example.com,B. J. Olson ††thanks: Email address for correspondence: firstname.lastname@example.org,S. K. Lele, S. Hickel, X. Y. Hu and N. A. Adams
Key words: Richtmyer-Meshkov instability, turbulent mixing, compressible turbulence, shock waves, large-eddy simulation
The Richtmyer-Meshkov instability (Richtmyer 1960; Meshkov 1969) is a hydrodynamic instability that occurs at the interface separating two fluids of different densities. It shows similarities with the Rayleigh-Taylor instability (Rayleigh 1883; Taylor 1950) where initial perturbations at the interface grow and eventually evolve into a turbulent flow field through the transfer of potential to kinetic energy. In the limit of an impulsive acceleration of the interface, e.g., by a shock wave, the instability is referred to as Richtmyer-Meshkov instability (RMI). In RMI, baroclinic vorticity production at the interface is caused by the misalignment of the pressure gradient () associated with the shock wave and the density gradient () of the material interface. The baroclinic vorticity production term is the initial driving force of RMI. See Brouillette (2002) and Zabusky (1999) for a comprehensive review.
RMI occurs on enormous scales in astrophysics (Arnett et al. 1989; Arnett 2000; Almgren et al. 2006), on intermediate scales in combustion (Yang et al. 1993; Khokhlov et al. 1999) and on very small scales in inertial confinement fusion (Lindl et al. 1992; Taccetti et al. 2005; Aglitskiy et al. 2010). Due to the fast time scales associated with RMI, laboratory experimental measurements have difficulties to characterise quantitatively initial perturbations of the material interface and to capture the evolution of the mixing zone. General insight into the flow physics of RMI relies to a considerable extent on numerical investigations, where large-eddy simulations (LES) have become an accepted tool during the past decade.
Hill et al. (2006) performed a rigorous numerical investigation of RMI with re-shock. The authors used an improved version of the TCD-WENO hybrid method of Hill & Pullin (2004). The method employs a switch to blend explicitly between a tuned centered-difference (TCD) stencil in smooth flow regions and a weighted essentially non-oscillatory (WENO) shock capturing stencil at discontinuities. The TCD-WENO hybrid method is used together with the stretched-vortex model (Pullin 2000; Kosović et al. 2002) for explicitly modelling the subgrid interaction terms. This approach was also used by Lombardini et al. (2011) to study systematically the impact of the Atwood number for a canonical three-dimensional numerical setup, and for LES of single-shock (i.e. without re-shock) RMI (Lombardini et al. 2012).
Thornber et al. (2010) studied the influence of different three-dimensional broad- and narrowband multimode initial conditions on the growth rate of a turbulent multicomponent mixing zone developing from RMI. In a later study (Thornber et al. 2011) the same authors presented a numerical study of a re-shocked turbulent mixing zone, and extended the theory of Mikaelian and Youngs to predict the behaviour of a multicomponent mixing zone before and after re-shock, c.f. Mikaelian (1989) and Thornber et al. (2010). They used an implicit LES (Drikakis 2003; Thornber et al. 2008; Drikakis et al. 2009) approach based on a finite-volume Godunov-type method to solve the Euler equations with the same specific heat ratio for both fluids.
In a recent investigation Weber et al. (2013) derived a growth-rate model for the single-shock RMI based on the net mass flux through the centre plane of the mixing zone. Here, the compressible Navier-Stokes equation were solved by a -order compact difference scheme for spatial differentiation. Artificial grid dependent fluid properties, proposed by Cook (2007), were used for shock and material-interface capturing as well as for subgrid-scale modelling.
Grid-resolution independent statistical quantities of the single-shock RMI were presented by Tritschler et al. (2013a). The kinetic energy spectra exhibit a Kolmogorov inertial range with scaling. The spatial flux discretisation was performed in characteristic space by an adaptive central-upwind -order accurate WENO scheme (Hu et al. 2010) in the low-dissipation version of Hu & Adams (2011).
LES relies on scale separation where the energy containing large scales are resolved and the effect of non-resolved scales is modelled either explicitly or implicitly. However, turbulent mixing initiated by RMI for typical LES mainly occurs at the marginally resolved or non-resolved scales. The interaction of non-resolved small scales with the resolved scales as well as the effect of the interaction of non-resolved scales with themselves is modelled by the employed subgrid-scale model. Moreover, discontinuities such as shock waves and material interfaces need to be captured by the numerical scheme. Due to the broad range of scales coarse-grained numerical simulations of RMI strongly rely on the resolution capabilities for the different types of subgrid scales (turbulent small scales, shocks, interfaces) of the underlying numerical scheme.
So far, research mainly focused on the identification and quantification of parameters that affect the evolution of Richtmyer-Meshkov unstable flows. The influence of the Atwood number (Lombardini et al. 2011), the Mach number (Lombardini et al. 2012) as well as the specific initial interface perturbations (Thornber et al. 2010; Grinstein et al. 2011; Schilling & Latini 2010) on the temporal evolution of the instability have been investigated. Results from numerical simulations have been compared to experiments (Schilling & Latini 2010; Hill et al. 2006; Tritschler et al. 2013b) and theoretical models have been derived (Thornber et al. 2011; Weber et al. 2013). These investigations have assumed, based on standard arguments such as empirical resolution criteria, that the marginally and non-resolved scales have a negligible effect on the resolved scales, and therefore on the evolution of the instability. Uncertainties introduced by the numerical method, i.e., the subgrid-scale regularisation and truncation errors, have not yet been investigated systematically. There is, however, strong evidence that numerical model uncertainty can significantly affect the linear and non-linear stages of evolution, and in particular the mixing measures. In fact, it is unclear how subgrid-scale regularisation and dispersive or dissipative truncation errors can affect the resolved scales and turbulent mixing measures.
In the present investigation, two independently developed and essentially different numerical methods are employed to study the prediction uncertainties of RMI simulations. The first method has a dominantly dissipative truncation error at the non-resolved scales, whereas the second one exhibits a more dispersive behaviour. At the marginally resolved scales the numerical truncation error is not small and the particular character of the truncation error is essential for the implicit modelling capabilities of the method, and thus also affects the resolved-scale solution. For the purpose of investigating this effect integral and spectral mixing metrics as well as probability density functions are analysed on four computational grids with resolutions ranging from to . The simulations employing two different numerical methods on a very fine grid-resolution of provide a data set with high confidence in the results.
We emphasise that the purpose of this study is (i) to present RMI results with a clear identification of the resolved-scale range by systematic grid refinement, (ii) to assess the physical effects of numerical subgrid-scale regularisations on the marginally resolved and on the non-resolved scale range. We do not intend to propose or improve a certain subgrid-scale model or regularisation scheme.
The paper is structured as follows: The governing equations along with the employed numerical models are described in section 2. Details about the computational domain and the exact generic initial conditions are given in section 3. Results are presented in section 4 and the key findings of the present study are discussed in section 5.
2 Numerical model
2.1 Governing equations
We solve the three-dimensional multicomponent Navier-Stokes equations
In (2.1) is the velocity vector, is the pressure, is the total energy, is the mixture density, is the mass fraction and is the diffusive mass flux of species with as the total number of species. The identity matrix is .
The viscous stress tensor for a Newtonian fluid is
with the mixture viscosity and the strain rate tensor . According to Fourier’s law we define the heat flux as
and the interspecies diffusional heat flux (Cook 2009) as
indicates the effective binary diffusion coefficient of species , and is the individual species enthalpy. The equations are closed with the equation of state for an ideal gas
where is the ratio of specific heats of the mixture and is the internal energy
2.2 Numerical methods
2.2.1 The Miranda simulation code
The Miranda simulation code has been used extensively for simulating turbulent flows with high Reynolds numbers and multi-species mixing (Cabot & Cook 2006; Cook et al. 2004; Olson & Cook 2007; Olson et al. 2011; Weber et al. 2013). Miranda employs a -order compact difference scheme (Lele 1992) for spatial differentiation and a five stage, -order Runge-Kutta scheme (Kennedy et al. 2000) for temporal integration of the compressible multicomponent Navier-Stokes equations. Full details of the numerical method are given by Cook (2007) which includes an -order compact filter which is applied to the conserved variables each time step and smoothly removes the top 10% of wave numbers to ensure numerical stability. For numerical regularisation of non-resolved steep flow gradients, artificial fluid properties are used to damp locally structures which exist on the length scales of the computational mesh. In this approach, artificial diffusion terms are added to the physical ones which appear in equations. (2.0), (2.0) and (2.0) as
This LES method employing artificial fluid properties was originally proposed by Cook (2007) but has been altered by replacing the (magnitude of the strain rate tensor) with in the equation for . Mani et al. (2009) showed that this modification substantially decreases the dissipation error of the method. Here we give the explicit formulation of the artificial terms on a Cartesian grid
where is the magnitude of the strain rate tensor, is the local grid spacing, is the sound speed and is the time-step size. The polyharmonic operator, , denotes a series of Laplacians, e.g., corresponds to the biharmonic operator, . The overbar denotes a truncated-Gaussian filter applied along each grid direction as in Cook (2007) to smooth out sharp cusps introduced by the absolute value operator. In LES of RMI, acts as the shock capturing scheme. The is primarily used as a numerical stabilisation mechanism rather than an SGS model. The artificial shear viscosity is found to not be needed to maintain numerical stability in the current calculations and its inclusion has a small impact on the solution. The dissipation of the vortical motion is primarily dependent on the -order compact filter.
2.2.2 The INCA simulation code
The INCA simulation code is a multi-physics simulation method for single- and multicomponent turbulent flows. With respect to the objective in this paper it has been tested and validated for shock induced turbulent multi-species mixing problems at finite Reynolds numbers (Tritschler et al. 2013b, a, 2014).
For all simulations presented in this paper we use a discretisation scheme that employs for the hyperbolic part in (2.1) a flux projection on local characteristics. The Roe-averaged matrix required for the projection is calculated for the full multi-species system (Roe 1981; Larouturou & Fezoui 1989; Fedkiw et al. 1997). The numerical fluxes at the cell faces are reconstructed from cell averages by the adaptive central-upwind -order weighted essentially non-oscillatory (WENO-CU6) scheme (Hu et al. 2010) in its scale-separation formulation by Hu & Adams (2011).
The fundamental idea of the WENO-CU6 scheme is to use a non-dissipative -order central stencil in smooth flow regions and a non-linear convex combination of -order stencils in regions with steep gradients. The reconstructed numerical flux at the cell boundaries is computed from
where is the weight assigned to stencil with the -degree reconstruction polynomial approximation for . In the WENO-CU6 framework the weights are given by
with being a small positive number . The optimal weights are defined such that the method recovers the -order central scheme in smooth flow regions. The constant parameters in (2.0) are set to and , see Hu & Adams (2011) is a reference smoothness indicator that is calculated from a linear combination of the other smoothness measures with
is also calculated from (2.0) but with the -degree reconstruction polynomial approximation of the flux which gives the 6-point stencil for the -order interpolation.
After reconstruction of the numerical fluxes at the cell boundaries the fluxes are projected back onto the physical field. A local switch to a Lax-Friedrichs flux is used as entropy fix, see, e.g., Toro (1999). A positivity-preserving flux limiter (Hu et al. 2013) is employed in regions with low pressure or density, maintaining the overall accuracy of the -order WENO scheme. It has been verified that the flux limiter has negligible effect on the results, and avoids excessively small time step sizes. Temporal integration is performed by a -order total variation diminishing Runge-Kutta scheme (Gottlieb & Shu 1998).
3 Numerical setup
3.1 Computational domain
We consider a shock tube with constant square cross section. The fine-grid domain extends in the - and -direction symmetrically from to and from to in the x-direction. An inflow boundary condition is imposed far away from the fine-grid domain in order to avoid shock reflections. To reduce computational costs a hyperbolic mesh stretching is applied between the inflow boundary and . is set to , and . At the boundaries normal to the - and -direction, periodic boundary conditions are imposed and an adiabatic-wall boundary at the end of the shock tube at is used. A schematic of the computational domain is shown in figure 1.
The fine-grid domain is discretised by four different homogeneous Cartesian grids with , , and cells in the y- and z-direction and , , and cells in the x-direction resulting in cubic cells of size . The total number of cells in the fine-grid domain amounts to for the coarsest resolution and to for the finest resolution.
3.2 Initial conditions
We consider air as mixture of nitrogen () and oxygen () with (in terms of volume fraction) and . The equivalent mass fractions on the air side give and , i.e., . The heavy gas is modelled as a mixture of and acetone () with mass fractions and , i.e., . The material interface between light (air) and heavy gas is accelerated by a shock wave that is initialised at propagating in the positive -direction. The pre-shock state is defined by the stagnation condition and . The corresponding post-shock thermodynamic state is obtained from the Rankine-Hugoniot conditions
with . The initial data of the post-shock state of the light gas as well as the pre-shock state of the light and heavy gas are given in table 1.
|Quantity||Post-shock||Pre-shock light-gas side||Pre-shock heavy-gas side|
Tritschler et al. (2013a) introduced a generic initial perturbation of the material interface that resembles a stochastic random perturbation being, however, deterministic and thus exactly reproducible for different simulation runs. This multimode perturbation is given by the following function
with the constant amplitudes and and wavenumbers , and . The amplitudes and the phase shifts and are given by
For facilitating a grid sensitivity study we impose an initial length scale by prescribing a finite initial interface thickness in the mass fraction field as
with being the characteristic initial thickness. The individual species mass fractions are set as
Figure 2 shows the initial condition in terms of the power spectrum of density for Miranda and INCA at all grid resolutions. The initial perturbation given in (3.0) and shown in figure 2 has been designed with the objective to obtain a reproducible and representative data set. Nevertheless, we cannot exclude that some of the observations presented in this paper do not apply to very different initial perturbations.
For exploring the effect of the finite truncation error arising from grid resolution and numerical method, four meshes were used to compute the temporal evolution of RMI with both Miranda and INCA. The simulation reaches which is well beyond the occurrence of re-shock at . At this stage, the effects of reflected shock waves and expansion waves on the shock location has become small as the shock wave are attenuated with each subsequent reflection. The space-time () diagram shown in figure 3 depicts the propagation of the shock wave and interface during the simulation.
The initial conditions described in the previous section are entirely deterministic and due to their band-limited representations are identically imposed at the different grid resolutions and for the two numerical methods. Therefore, the obtained results exhibit uncertainties only due to the numerical method and due to grid resolution, but exclude initial-data uncertainties.
For illustration we show the three-dimensional contour plots of species mass fraction of the heavy gas obtained with Miranda and INCA, respectively, in figure 4. Similarities at the large scales are clearly visible after re-shock, but also differences exist at the fine scales, more clearly visible from the inset.
4.1 Integral quantities
Integral measures of the mixing zone are presented here for both numerical models and all resolutions. Often, these time dependent integral measures are the only metrics available for comparison with experiment and are therefore of primary importance for validation.
Figure 5 shows the transition process predicted by the reference grid with a resolution of cells in the transverse directions. The numerical challenge, prior to re-shock, is to predict the large-scale non-linear entrainment and the associated interface steepening. The interface eventually becomes under-resolved when its thickness reaches the resolution limit of the numerical scheme and further steepening is prevented by numerical diffusion. The equilibrium between interface steepening and numerical diffusion occurs later in time as the grid is refined. The accurate prediction of the interface steepening phenomenon is one of the main challenges in modelling pre-transitional RMI where large-scale flow structures are still regular. This is because the numerical model largely determines the time when mixing transition occurs. In nature, mixing transition is due to the presence of small-scale perturbations whereas in numerical simulation, the transition is triggered by back-scatter from the under-resolved scales as predicted by the particular numerical model. Hence, details of mixing transition of the material interface evolve differently for the two codes.
Nevertheless, similarities before re-shock are striking and large-scale similarities in the resolved wavenumber range even persist throughout the entire simulation time. Following re-shock the large interfacial scales break down into smaller scales and develop a turbulent mixing zone as can be seen in figure 4 and figure 5. By visual inspection of figure 5 one finds that the post-re-shock turbulent structures are very similar, whereas the long term evolution of the small scales appears to be different between the codes. Differences in the observed flow field at may indicate slightly different effective Reynolds numbers for the two numerical methods and therefore they also exhibit different decay rates of enstrophy (Dimotakis 2000; Lombardini et al. 2012), as can be seen in figure 8 after re-shock.
The mixing width is a length scale that approximates the large-scale temporal evolution of the turbulent mixing zone. It is defined as an integral measure by
where denotes the ensemble average in the cross-stream yz-plane. For a quantity it is defined by
The mixing width plotted in figure 6(a) shows that data from both numerical methods converge to a single solution throughout the entire simulation time. Furthermore, it is observed that even with very-high-order models a minimum resolution of appears to be necessary for an accurate prediction of the mixing-zone-width. As will be shown later, coarser grids not only tend to overpredict the growth of the mixing zone but also molecular mixing.
Figure 6(b) shows the mixing zone width time evolution on a log-log scale. The (bubble) growth rate model of Zhou (2001) predicts accurately the pre-re-shock mixing zone growth rate that is consistently recovered by both numerical methods as . However, this is, according to Zhou (2001), the growth rate that is associated with turbulence of Batchelor type (Batchelor & Proudman 1956) with as . The kinetic energy spectra in the present investigation are of Saffman type (Saffman 1967b, a) with as (Tritschler et al. 2013a) for which Zhou (2001) predicts a growth that scales with . The present growth rates are also in good agreement with the experimental and numerical results of Dimonte et al. (1995) with and and their model predictions , where accounts for the time the shock needs to traverse the interface. As the mixing zone has not yet reached self-similar evolution, the initial growth rate depends on the specific initial conditions.
Llor (2006) found that the self-similar growth rate of the energy-containing eddies, i.e., the integral length scale, for incompressible RMI at vanishing Atwood number should scale as with , if the turbulence kinetic energy decays as . These growth rates slightly differ from the growth rate prediction for homogeneous isotropic turbulence by the same author. The predictions of Llor (2006), however, are at odds with Kolmogorov’s classical decay law (Kolmogorov 1941) for turbulence kinetic energy and more recent investigations of decaying isotropic turbulence by Ishida et al. (2006) and Wilczek et al. (2011), which found Kolmogorov’s decay law to hold if the Loitsyansky integral is constant and if the Taylor-scale Reynolds number exceeds . Based on Rayleigh-Taylor experiments driven by either sustained or impulsive acceleration at various Atwood numbers, Dimonte & Schneider (2000) found scaling laws for the bubble and spike growth rate. For the present density ratio the exponents become for the former and for the latter. The late-time mixing zone growth rate is therefore expected to correlate with the spike growth rate. The late-time growth rate prediction of the present work is , i.e., , once the turbulent mixing zone is fully established. is a virtual time origin set to . This is consistent with the mixing zone width growth rate predictions of Llor (2006) and the late time growth rate predictions of Dimonte & Schneider (2000), but underestimates the predictions of Zhou (2001), with a scaling of with long after re-shock, once the non-linear time scale has become the dominant time scale. In the numerical investigation of Lombardini et al. (2012) the authors found the mixing zone width to grow as . Before re-shock the infrared part of the kinetic energy spectrum, see figure 14, exhibits a range, for which a post-re-shock growth rate of is predicted by the model of Youngs (2004), which is in good agreement with the present data.
and quantifies the amount of mixed fluid within the mixing zone. It can be interpreted as the ratio of molecular mixing to large-scale entrainment by convective motion.
As bubbles of light air and spikes of heavy gas begin to interfuse, the initially mixed interface between the fluids steepens and the fluids become more segregated on the molecular level, see figure 7(a). The molecular mixing fraction reaches its minimum at before Kelvin-Helmholtz instabilities lead to an increase of molecular mixing. The onset of secondary instabilities is very sensitive to the numerical method as the numerical scheme determines how sharp the material interface can be represented or whether numerical diffusion or dispersion effects lead to an early mixing transition.
After re-shock molecular mixing is strongly enhanced and reaches its maximum of by the end of the simulation. This finding is consistent with Lombardini et al. (2012) who also found an asymptotic late-time mixing behaviour with independent of the shock Mach number but without re-shock. The asymptotic limit is already accurately calculated on grid-resolutions of . As the second shock wave compresses the mixing zone, the instability becomes less entrained yet equally diffused (at least in the y- and z-directions) and therefore causes a steep rise in . A gradual increase of the mixing fraction after the steep rise occurs as the mixing zone becomes more homogeneously distributed (Thornber et al. 2011) due to turbulent motion.
The temporal evolution of the scalar dissipation rate is plotted in figure 7(b) and is derived from the advection-diffusion equation for a scalar. The instantaneous scalar dissipation rate of the three-dimensional RMI is estimated from the concentration field as
which quantifies the rate at which mixing occurs. For consistency of post-processing, a -order central-difference scheme has been used for the calculation of the spatial derivatives in (4.0) and (4.0) for all simulation data sets. Note that the order of the finite-difference scheme with which the gradients in (4.0) and (4.0) are approximated affect their results.
The variation of the scalar dissipation rate with grid resolution before re-shock is largely due to the under-resolved material interface and the onset of mixing transition. Mixing is strongly enhanced after the second shock-interface interaction, but the mixing zone is also confined to a much smaller region which results in a decrease of the integral . Also, only represents the resolved part of the dissipation rate and therefore certainly underestimates the true value.
The turbulence kinetic energy (TKE) and the enstrophy () are integrated over cross-flow planes in the mixing zone that satisfy
This region is referred to as the inner mixing zone (IMZ) in the following.
Baroclinic vorticity is deposited at the material interface during shock passage. The amount of generated vorticity scales directly with the pressure gradient of the shock wave and the density gradient of the material interface. The enstrophy is calculated by
where is the vorticity.
As can be seen from figure 8, the enstrophy also exhibits a strong grid dependency. Fully grid-converged results are only obtained for times up to . As the interface steepens due to strain and shear the effective interface thickness is determined by numerical diffusion which appears to occur at . This is consistent with the evolution of shown in figure 7(a). Following Hahn et al. (2011) and Youngs (2007) integration of enstrophy with a theoretical scaling of up to the cut-off wavenumber yields a proportionality between enstrophy and grid resolution as . From this follows an increase of enstrophy by a factor of about from one grid resolution to the next finer, which is in good agreement with the present data.
The amount of turbulence kinetic energy created by the impulsive acceleration of the interface is calculated as
The fluctuating part of a quantity is calculated from
where is the Favre-average of .
Grid-converged turbulence kinetic energy TKE is obtained on grids with a minimum resolution of , see figure 8. This is consistent with the convergence rate of the mixing zone width. The total TKE deposited in the IMZ by the first shock-interface interaction can be seen in the inset of figure 8. The re-shock occurring at deposits approximately times more TKE than the initial shock wave. Hill et al. (2006) found a similar relative increase by the re-shock at the same shock Mach number. A significant decay in energy occurs immediately following re-shock. The material interface interacts with the first expansion fan, see figure 3, and results in a further increase in TKE between . The amount of energy deposited by the first expansion wave, however, is much weaker than that deposited by the reflected shock wave. Grinstein et al. (2011) and (Hill et al. 2006) found the amplification of TKE by the first rarefaction to be much stronger than for our data. Such differences are not surprising as Grinstein et al. (2011) reported a strong dependence of energy deposition on the respective initial interface perturbations. After the first expansion wave has interacted with the interface, TKE decays slowly and the pressure gradients associated with the subsequent rarefactions are too shallow to generate any further noticeable increase in TKE.
Lombardini et al. (2012) found the decay rate of TKE to be larger than approaching . In our data, the late time TKE decay is also approximately with being the virtual time origin, see figure 9. This scaling would be characteristic for Batchelor-type turbulence (Batchelor & Proudman 1956) with a constant Loitsyansky integral (Kolmogorov 1941; Ishida et al. 2006) in contrast to typical for turbulence of Saffman-type (Saffman 1967b, a). The similar TKE decay rates in homogeneous isotropic turbulence and Lombardini et al. (2012) and the present data suggest that the dissipation mechanism in RMI is quite similar to that in homogeneous isotropic turbulence.
In the limit of a self-similar quasi-isotropic state the temporal evolution of the integral length scale is related to the evolution of TKE in the mixing zone. From the growth rate of the integral scale follows as . Llor (2006) derived a maximum decay rate of turbulence kinetic energy that corresponds to a growth rate scaling of the energy-containing eddies of . These predictions are in excellent agreement with the growth rate predictions of the mixing zone width of the present investigation, see figure 6, and the decay rate of TKE, see figure 9.
The scalings indicated for the growth rate of the mixing zone and the decay rate of TKE in figures 6 and 9 were not fitted in a strict sense. They merely serve as reference for comparison with incompressible isotropic decaying turbulence. The narrow data range of only after re-shock for which the flow exhibits a self-similar regime precludes any precise estimates for decay and growth rate laws.
4.2 Anisotropy and inhomogeneity of the mixing zone
In the following the anisotropy in the mixing zone is investigated. We define the local anisotropy as
where corresponds to having all turbulence kinetic energy in the streamwise velocity component , whereas corresponds to having no energy in the streamwise component. In figure 10 (a) and (b) we show the -plane averaged anisotropy as a function of the dimensionless mixing zone coordinate and time from Miranda and INCA. The dimensionless mixing zone coordinate is defined as
with being the x-location where is maximal.
The light gas side of the mixing zone remains more anisotropic than the heavy gas side but with a homogeneous anisotropy distribution after re-shock on either side. The volume averaged anisotropy in the inner mixing zone is shown in figure 10 (c). No full recovery of isotropy of the mixing zone is achieved, and the re-shock does not significantly contribute in the sense of the volume averaged quantity , but leads to a stratified anisotropy distribution around the centre of the mixing zone. After an asymptotic limit of is reached, which temporally coincides with the onset of the self-similar decay of TKE, see figures 9 and 10. The positive value of implies that the streamwise component remains, despite re-shock, the dominant velocity component throughout the simulation time. Lombardini et al. (2012) also found a temporal asymptotic limit of the isotropisation process in their simulations. Grinstein et al. (2011) observed that the velocity fluctuations in the mixing zone are more isotropic when the initial interface perturbations also include short wavelengths in which case the authors nearly recovered full isotropy. When Grinstein et al. (2011) used long wavelength perturbations the mixing zone remained anisotropic except for a narrow range on the heavy gas side.
In order to quantify the homogeneity of mixing we calculate the density self-correlation (Besnard et al. 1992)
which is non-negative. corresponds to homogeneously mixed fluids with constant pressure and temperature. Large values indicate spatial inhomogeneities in the respective yz-plane. The density self-correlation has gained some attention in recent years and was subject in several experimental investigations of the RMI, see Balakumar et al. (2012); Balasubramanian et al. (2012, 2013); Orlicz et al. (2013); Tomkins et al. (2013); Weber et al. (2014).
Figure 11 (a) and (b) shows the density self-correlation normalised by the maximal value at time
as a function of the dimensionless mixing zone coordinate and time from Miranda and INCA. The largest values of are found around the centre of the mixing zone slightly shifted towards the heavy-gas side. peaks around the region where mixing between light and heavy gas occurs and tends to zero outside the mixing region, towards the respective pure gas side. Weber et al. (2014) observed in their experiment that the peak of the density self-correlation is initially shifted towards the light-gas side, but moves towards the center of the mixing zone with increasing time.
In contrast to the anisotropy, where the re-shock does not contribute to the isotropisation and which levels out after , the mixing zone becomes significantly more homogeneous after re-shock as can be observed from the temporal evolution of the volume average of the density self-correlation in the inner mixing zone . Following re-shock the fluids become more and more mixed, see figure 11 (c) with a value of at the latest time. The measured values of the density self-correlation in the single shock-interface interaction experiment of Weber et al. (2014) at are in good agreement with our simulated values at late times , whereas at the lower Mach number Weber et al. (2014) observed a more inhomogeneous mixing zone . These values are significantly larger than those measured for instance in the shock-gas-curtain experiments of Orlicz et al. (2013) and Tomkins et al. (2013).
4.3 Spectral quantities
From homogeneous isotropic turbulence it is well known that vorticity exhibits coherent worm-like structures with diameter on the order of the Kolmogorov length scale and of a length that scales with the integral scale of the flow. The work of Jiménez et al. (1993) suggests that theses structures are especially intense features of the background vorticity and independent of any particular forcing that generates the vorticity. In contrast to forced homogeneous isotropic turbulence where self-similar stationary statistics are achieved, shock-induced turbulent mixing is an inhomogeneous anisotropic unsteady decay phenomenon. Nevertheless, homogeneous isotropic turbulence is used as theoretical framework for most of the numerical analysis of RMI. However, it is unclear at what time and at what locations the mixing zone exhibits the appropriate features and if homogeneous isotropic turbulence is achieved at all. A fully isotropic mixing zone is never obtained, as the anisotropy, even though decreasing with time, reaches an asymptotic limit at .
The temporal evolution of the initial perturbation is depicted in figure 12. Before re-shock the dominant modes of the initial perturbation slowly break down. After re-shock, however, the additional vorticity deposited during the second shock-interface interaction rapidly destroys structures generated by the initial perturbation and initial shock, leading to a self-similar decay after .
Thornber et al. (2010) and Thornber et al. (2012) found, formally in the limit of infinite Reynolds numbers, a persistent scaling of the turbulence kinetic energy spectrum as well as a spectrum with a spectrum at high wavenumbers that covers more and more of the spectrum as time proceeds. Furthermore, the same authors (Thornber et al. 2011) found (depending on the initial conditions) a or a scaling range after re-shock. Long after re-shock however, these scalings return to a scaling at intermediate scales and to a scaling at high wavenumbers, close to the cut-off wavenumber. The authors evaluated the radial spectra either in the centre of the mixing zone or averaged over a fixed number of yz-planes within the mixing zone.
A different scaling behaviour was observed by Lombardini et al. (2012) and Hill et al. (2006) who found in their multicomponent LES at finite Reynolds numbers a scaling in the centre of the mixing zone. Whereas Cohen et al. (2002) found a scaling range for the single-shock RMI averaged over four transverse slices within the mixing zone. In a recent experimental investigation of a shock accelerated shear-layer Weber et al. (2012) showed a inertial range followed by an exponential decay in the dissipation range of the scalar spectrum. This result was numerically reproduced by Tritschler et al. (2013a). Here, the authors averaged over a pre-defined IMZ.
All spectra shown in this section are radial spectra with a radial wavenumber that is defined as . The radial spectra are averaged over all yz-planes within the IMZ in the x-direction that satisfy the condition in (4.0).
The radial power spectra of density are plotted in figure 13, where figures 13(a) and (b) show the spectra before and figures 13(c) and (d) after re-shock. Power spectra of density and mass fraction concentration (not shown) show a close correlation, even though they are not directly related as the mass fractions are constrained to be between zero and one.
Before re-shock, the dominant initial modes slowly break down and redistribute energy to smaller scales. Re-shock causes additional baroclinic vorticity production with inverse sign that results in a destruction process of the pre-shock structures, see also figure 12. This process in conjunction with a vorticity deposition that is one order of magnitude larger than the pre-shock deposition leads to rapid formation of complex disordered structures, which eliminates most of the memory of the initial interface perturbation as can be seen in figures 12, 13 and 14. Schilling et al. (2007) reported that during re-shock vorticity production is strongly enhanced along the interface where density gradients and misalignment of pressure and density gradients is largest. The vorticity deposited by the re-shock transforms bubbles into spikes and vice versa, which subsequently results in more complex and highly disordered structures.
At late times the power spectra of density appear to be more shallow than , and rather approach as was found by Cohen et al. (2002).
The smallest length scale in scalar turbulence is the Batchelor scale. For isotropic turbulence and Schmidt numbers of order unity it has the same order of magnitude as the Kolmogorov microscale . Therefore, the TKE spectra are closely correlated with the scalar power spectra. Figure 14 shows the spectra of TKE before and after re-shock. The significant increase in TKE is mainly due to the interaction of the enhanced small scale structures with comparatively steep density gradients and the reflected shock wave. The re-shock at leads to a self-similar lifting of the spectrum, see figure 14. The destruction process of the vortical structures initiated by the re-shock leads to the formation of small scales, which rapidly remove the memory of the initial condition. The intense fluctuating velocity gradients past re-shock are rapidly smoothed out by viscous stresses. This results in a fast decay of the turbulence kinetic energy following the first after re-shock, see figure 8 and figures 14(c) and (d).
The sharp drop-off of the spectral energy in figures 14 and 13 in the Miranda data at high wavenumbers is due to the filtering operator of the numerical method. Opposite behaviour, that is an increase of spectral energy at the highest wavenumbers, is observed for the less dissipative INCA code, where the spurious behaviour at the non-resolved scales is mainly dispersive.
The scaled turbulence kinetic energy spectra represent the effective energy contributed by each mode. Artifacts of the initial conditions still exist immediately before re-shock at as can be seen in figure 15(a), where most energy is contained at mode . At re-shock, baroclinic vorticity is deposited at the interface and the energy containing wavenumber range immediately widens as vortex stretching and tangling introduce new scales and higher vorticity. This broader profile is plotted in figure 15(b) which clearly shows that the relative difference between the imposed initial length scale and the remaining length scales (both larger and smaller) is vanishing. Indeed, as the mixing layer fully transitions to turbulence, the flow reaches a self-similar state where the memory of initial perturbations is lost.
The spectra of the scalar dissipation rate in figure 16 quickly builds up in the cut-off wavenumber range after the initial shock impact, see figure 16(b). After re-shock and at late time, see figure 16(c)-(d), the inertial subrange broadens to wavenumbers where numerical dissipation damps out structures. The inertial range is observed to scale with after re-shock, which is consistent with the scaling observed for and . For the resolved wavenumbers, there is good agreement between both codes at the finest two resolutions. Differences observed in figure 16 (b) are also reflected in figures 7 and 11. A sharper material interface and the associated segregation of the fluids lead to a higher scalar dissipation rate (), whereas at late times the difference in the scalar dissipation rate does not significantly influence the mixing measures and .
Larger quantitative differences are observed in the power spectra of enstrophy shown in figure 17. Immediately after either of the shock-interface interactions the quantitative agreement between the predicted enstrophy levels is excellent, see figures 17(a) and (c). The observed scalings of the inertial range following re-shock are predicted consistently and agree with the inertial range scalings for the scalar dissipation rate . However, the temporal decay of the small-scale enstrophy is significantly different for either code as can be seen immediately before re-shock and long after re-shock in figures 17(b) and 17(d), respectively.
In isotropic homogeneous turbulence, the scaled spectra of the enstrophy, see figure 18, has a single peak at the wavenumber where the dissipation range begins. Therefore, under grid refinement this peak will shift to higher wavenumbers and magnitudes as smaller scales are captured. The peak at is associated with the initial perturbation and disappears after re-shock as the flow becomes turbulent, see figure 18. Good agreement for lower wavenumbers is observed between codes and resolutions. Larger differences are observed at high wavenumbers where the dependence on numerical dissipation is greatest. At the peak in the scaled enstrophy spectra is at for both codes at the highest resolution. Later, at this peak has shifted to in INCA, whereas in Miranda, there is no apparent shift, although both have substantially decayed in magnitude.
As RMI is a pure decay process after re-shock, differences in the numerical approach become most apparent at late times. The numerical models of this study predict different turbulence decay rates as is evident from differences in the enstrophy spectrum, figure 17 and figure 18, and in TKE, figure 9. The differences in enstrophy () and scalar dissipation rate () have a qualitative effect that becomes apparent in the fine scale structures of figure 5 at . Although INCA resolved less scales with smaller enstrophy levels, it does resolve steeper mass fraction gradients, which is reflected in the higher and higher levels of . Although it is unclear which dissipation rate (scalar or kinetic) has most effect on the mixing process, both are important (Dimotakis 2000).
4.4 Probability density functions
The bin size for computing discrete probability density function (pdf) is defined as for a quantity . The number of bins for all quantities and all grid-resolutions is . Each discrete value of is distributed into the bins, yielding a frequency for each bin. The is then defined by
such that with as the total number of cells in the IMZ that fall within the range . The limits and are held constant for all resolutions and times.
The of the heavy-gas mass fraction is constrained to be . Figure 19 shows at times before re-shock (, ) and following re-shock (, ). From figure 19(a) it is evident that at early times following the initial shock-interface interaction the IMZ consists mostly of segregated fluid as the large peaks at the bounds indicate. Before re-shock, interspecies mixing is largely dominated by the inviscid linear and nonlinear entrainment. Molecular diffusion processes have not yet had enough time to act, see figure 19(b). Following re-shock, a fundamental change in the of () is observed, see figures 19(c) and (d). The additional vorticity deposited by the re-shock leads to rapid formation of small and very intense vortical structures that lead to very effective mixing and destruction of the initial interface perturbation. The takes a uni-modal form at as it was also reported by Hill et al. (2006). The peak value, however, is not as well correlated with the average value of the mixture mass fraction as it was reported by Hill et al. (2006). With our data the peak value is slightly shifted towards the heavy-gas side centered around . The degree of convergence between codes and resolutions is reassuring at . Note that is a very sensitive measure of the light-heavy gas mixing.
The rarefaction wave at does not significantly contribute to the mixing, as it is not as pronounced as found in comparable investigations (Grinstein et al. 2011; Hill et al. 2006). Long after re-shock the mixing process continues, which is reflected in narrower tails of . The peak value of predicted by Miranda now coincides with the average value of the mixture mass fraction. In the INCA results this value remains slightly shifted towards the heavy gas side. However, the bimodal character of reported by Hill et al. (2006), who used - as light-heavy gas, is not observed on the finest grid. Despite the strong mixing past re-shock the turbulent mixing zone remains inhomogeneous until the end of the simulation time, which makes the observed very sensitive to their location of evaluation within the mixing layer.
The of the normalised vorticity is constrained between with , where is the initial shock velocity and is a characteristic length scale of the perturbations taken as where is the width of the domain in the transverse direction and .
Figure 20 shows the of the normalised vorticity . Before re-shock, mixing is driven by weak large-scale vortices, see figures 20(a) and (b). Following re-shock, however, structures with very intense vorticity develop with a dual-mode shape in at on the finest grid. The early times after the second shock-interface interaction are again consistently predicted by both codes, figure 20(c). Nevertheless, we observe larger differences for at . The peak and distribution of from Miranda are shifted towards larger values of as compared to INCA. This supports the previous observation that the vorticity decay is affected by the numerical approach. The difference in the vorticity intensity observed in figure 20, however, does not lead to noticeable differences for the integral mixing measures shown in figures 7(a) and 11 or for the integral length scale in figure 6.
We have investigated the shock-induced turbulent mixing between a light (, ) and heavy (, acetone) gas in highly resolved numerical simulations. The mixing was initiated by the interaction of a shock wave with a deterministic multimode interface. After the initial baroclinic vorticity deposition, the shock wave is reflected at the opposite adiabatic wall boundary. The reflected shock wave impacts the interface (re-shock) and deposits additional vorticity with enstrophy that is more than two orders of magnitude larger than that of the initial vorticity deposition. The transformation of spike structures into bubbles and vice versa in conjunction with a large increase in vorticity results in the formation of disordered structures which eliminate most of the memory of the initial interface perturbation.
A proposed standardised initial condition for simulating the Richtmyer-Meshkov instability (RMI) has been assessed by two different numerical approaches, Miranda and INCA, over a range of grid resolutions. The deterministic interface definition allows for spectrally identical initial conditions for different numerical models and grid resolutions. A direct comparison shows that larger energy containing scales are in excellent agreement. Different subgrid-scale regularisation affect marginally resolved flow scales, but allow for a clear identification of a resolved scale-range that is unaffected by the subgrid-scale regularisation.
Mixing widths are nearly identical between the two approaches at the highest resolution. At lower resolutions, the solutions differ, and we found a minimum resolution of to be necessary in order to produce reasonable late-time results. The initial mixing zone growth rate scaled with , whereas long after re-shock the predicted growth rate was . The decay of turbulence kinetic energy was also found to be consistent and in good agreement between the approaches. The decay scaled with at late times, that corresponds to a growth rate scaling of the energy-containing eddies of . The agreement in the large scales of the solution between the two approaches is striking and has not been observed before.
Previous work on three-dimensional LES of RMI has examined numerical dependence only indirectly. E.g. Thornber et al. (2010) performed a code comparison of single-shock RMI. However, the initialisation was different for the two codes and the purpose was not to quantify the effect of different numerical methods. Accordingly, the comparison of results at different resolutions for mean, spectral and gradient based quantities was limited. With the current work we have presented for the first time a comprehensive quantitative analysis of numerical effects on RMI.
Results conclusively show that the large scales are in excellent agreement for the two methods. Differences are observed in the representation of the material interface. We conclude that the numerical challenge, prior to re-shock, is to predict the large-scale non-linear entrainment and the associated interface sharpening. Under shear and strain the interface steepens and eventually becomes under-resolved with a thickness defined by the resolution limit of the numerical scheme. Therefore, the saturation of the interface thickness by the numerical method occurs later as the grid is refined. The molecular mixing fraction reached an asymptotic limit as after re-shock, which was already correctly calculated on grid-resolutions of . The -plane averaged anisotropy revealed that the mixing zone exhibits a stratified anisotropy distribution with lower anisotropy on the heavy-gas side and higher anisotropy on the light-gas side. Moreover, the volume-averaged anisotropy approached an asymptotic limit of , implying that the fluctuating velocity component remains the dominant component even after re-shock and that no full recovery of isotropy of the mixing zone is obtained. The density-self correlation has been investigated in order to better understand the mixing inhomogeneity in the mixing zone. The volume-averaged density self-correlation showed that the re-shock significantly increases mixing homogeneity approaching a value of at the latest time.
The spectra demonstrate a broad range of resolved scales, which are in very good agreement. Data also show that differences exist in the small-scale range. The frequency dependence of the velocity and density fluctuations shows the existence of an inertial subrange and that the two approaches agree at lower frequencies. The observed spectral scalings were consistent among the methods with .
Quantities that are gradient dependent and therefore more sensitive to small scales, such as the scalar dissipation rate and enstrophy, exhibit stronger dependence on numerical method and grid resolution. The flow field shows visual differences for the fine-scale structures at late times. The -order compact scheme and the explicit filtering and artificial fluid properties used in Miranda resolved more small scales in turbulence kinetic energy and enstrophy, whereas the -order WENO based scheme used in INCA resolves more of the small-scale scalar flow features as observed in the spectra of density and scalar dissipation rate. This result is somewhat intuitive given the numerics of the two codes. High-order compact methods are capable of resolving higher modes than explicit finite difference methods, (Lele 1992). Given that the artificial shear viscosity in Miranda has only a small effect on the solution compared to the effect of the -order filter, the primary difference, we conclude, of the resolving power between methods is due to the difference in order of accuracy and modified wavenumber profiles between the schemes. The compact finite difference method with high-order filtering appears to capture a broader range of dynamic scales at late times.
The statistics of heavy gas mass fraction revealed that the IMZ remains inhomogeneous until the end of the simulation and that the peak probability is centered around and thus is slightly shifted towards the heavy gas side. Although the overall quantitative agreement was very good, the of the vorticity showed larger differences once intense small-scale vortical structures exist. The decay of vorticity differs accordingly between the numerical methods.
The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, www.lrz.de). This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract number DE-AC52-07NA27344. VKT gratefully acknowledges the support of the TUM Graduate School. BJO thanks A. Cook and W. Cabot for valuable insight and for use of the Miranda code.
A Multicomponent mixing rules
The specific heat capacity of species is found by
where is the ratio of specific heats. The ratio of specific heats of the mixture follows as
where is the mass fraction of species and is the specific gas constant of the mixture with . The molar mass of the mixture is given by
For the gas mixture Dalton’s law shall be valid with . The mixture viscosity and the mixture thermal conductivity is calculated from (Reid et al. 1987)
The effective binary diffusion coefficients (diffusion of species into all other species) are approximated as (Ramshaw 1990)
where is the mole fraction of species . (A 0) ensures that the inter-species diffusion fluxes balance to zero.
B Molecular mixing rules
The viscosity of a pure gas is calculated from the Chapman-Enskog model (Chapman & Cowling 1990)
where is the collision diameter and is the collision integral
with , , , , and and . is the Lennard-Jones energy parameter, where is the minimum of the Lennard-Jones potential and is the Boltzmann constant.
The thermal conductivity of species is defined by
with as the species specific Prandtl number.
The mass diffusion coefficient of a binary mixture is calculated from the empirical law (Reid et al. 1987)
with the collision integral for diffusion
where and , , , , , , , and
The molecular properties of all species in the present study are given in Table 2.
- Aglitskiy et al. (2010) Aglitskiy, Y., Velikovich, A. L., Karasik, M., Metzler, N., Zalesak, S. T., Schmitt, A. J., Phillips, L., Gardner, J. H., Serlin, V., Weaver, J. L. & Obenschain, S. P. 2010 Basic hydrodynamics of richtmyerâmeshkov-type growth and oscillations in the inertial confinement fusion-relevant conditions. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 368 (1916), 1739–1768.
- Almgren et al. (2006) Almgren, A. S., Bell, J. B., Rendleman, C. A. & Zingale, M. 2006 Low Mach numbermodeling of type Ia supernovae I. Hydrodynamics. The Astrophysical Journal 637, 922.
- Arnett (2000) Arnett, D. 2000 The role of mixing in astrophysics. Ap. J. Suppl. 127, 213.
- Arnett et al. (1989) Arnett, W. D., Bahcall, J. N., Kirshner, R. P. & Stanford, E. W. 1989 Supernova 1987a. Annu. Rev. Astrom. Astrophys. 27, 629.
- Balakumar et al. (2012) Balakumar, B J, Orlicz, G C, Ristorcelli, J R, Balasubramanian, S, Prestridge, K P & Tomkins, C D 2012 Turbulent Mixing in a Richtmyer-Meshkov Fluid Layer after Reshock: Velocity and Density Statistics. J. Fluid Mech. accepted.
- Balasubramanian et al. (2013) Balasubramanian, S, Orlicz, GC & Prestridge, KP 2013 Experimental study of initial condition dependence on turbulent mixing in shock-accelerated richtmyer–meshkov fluid layers. Journal of Turbulence 14 (3), 170–196.
- Balasubramanian et al. (2012) Balasubramanian, S., Orlicz, G. C., Prestridge, K. P. & Balakumar, B. J. 2012 Experimental study of initial condition dependence on Richtmyer-Meshkov instability in the presence of reshock. Physics of Fluids 24, 034103.
- Batchelor & Proudman (1956) Batchelor, G. K. & Proudman, I. 1956 The large-scale structure of homogeneous turbulence. Phil. Trans. R. Soc. Lond. A 248, 369.
- Besnard et al. (1992) Besnard, Didier, Harlow, FH, Rauenzahn, RM & Zemach, C 1992 Turbulence transport equations for variable-density turbulence and their relationship to two-field models. NASA STI/Recon Technical Report N 92, 33159.
- Brouillette (2002) Brouillette, M 2002 The Richtmyer-Meshkov instability. Annu. Rev. Fluid Mech. (34), 445–468.
- Cabot & Cook (2006) Cabot, W. H. & Cook, A. W. 2006 Reynolds number effects on rayleigh–taylor instability with possible implications for type ia supernovae. Nature Physics 2 (8), 562–568.
- Chapman & Cowling (1990) Chapman, S. & Cowling, T. G. 1990 The Mathematical Theory of Non- Uniform Gases: An Accou