The mean momentum and heavy mass fraction, turbulent kinetic energy, and heavy mass fraction variance fields, as well as the budgets of their transport equations are examined several times during the evolution of a narrowband Richtmyer-Meshkov instability initiated by a Mach 1.84 shock traversing a perturbed interface separating gases with a density ratio of 3. The results are computed using the “quarter scale” data from four algorithms presented in the θ-group study of Thornber et al. [“Late-time growth rate, mixing, and anisotropy in the multimode narrowband Richtmyer-Meshkov instability: The θ-group collaboration,” Phys. Fluids 29, 105107 (2017)]. The present study is inspired by a previous similar study of Rayleigh-Taylor instability and mixing using direct numerical simulation data by Schilling and Mueschke [“Analysis of turbulent transport and mixing in transitional Rayleigh-Taylor unstable flow using direct numerical simulation data,” Phys. Fluids 22, 105102 (2010)]. In addition to comparing the predictions of the data from four implicit large-eddy simulation codes, the budgets are used to quantify the relative importance of the terms in the transport equations, and the balance of the terms is employed to infer the numerical dissipation. Terms arising from the compressibility of the flow are examined, in particular the pressure-dilatation. The results are useful for validation of large-eddy simulation and Reynolds-averaged modeling of Richtmyer-Meshkov instability.

Richtmyer-Meshkov instability (RMI) occurs when a perturbed interface between two fluids of differing properties is impulsively accelerated.1,2 This acceleration imparts a vorticity field near the interface, causing any perturbations on the interface to grow in time, regardless of the direction of the impulse relative to the layer. For initial perturbations with amplitude over wavelength a/λ ≪ 1, the growth rate of the instability is at first linear.1 As the perturbation amplitude becomes nonlinear, the perturbation growth rate decreases and the shear layers between the peaks and troughs roll up to form vortices. Shear along the interface also triggers smaller scale Kelvin-Helmholtz instabilities, which at later times contribute to the development of a range of vortical structures of differing sizes. For multimode perturbations, this process occurs simultaneously for a range of length-scales, leading to the development of an inhomogeneous turbulent mixing layer at late times. As the deposition of vorticity is impulsive, the layer with linear perturbations evolves through several growth regimes: (i) linear, (ii) nonlinear, (iii) nonlinear transitional (range of length-scales developing), and (iv) decaying inhomogeneous variable density turbulent layer.

This instability occurs in applications ranging from inertial confinement fusion3 and astrophysical flows4 to augmented mixing in scramjets and jet exhaust plumes.5 For a comprehensive recent review of research in this field, see Refs. 6 and 7. RMI differs significantly from Kelvin-Helmholtz instability (KHI) and closely related Rayleigh-Taylor instability (RTI) in that the driving is impulsive. With a sufficiently strong impulse, time to develop, and high Reynolds number, linear growth will be followed by a transition to turbulence, where velocity fluctuations decay in time and the layer grows at a relatively slow rate ∝ tθ with θ ≈ 0.2–1 (possible if dominated by linear modes, although a fully turbulent layer is limited to8θ ≤ 2/3) with a strong dependence on initial conditions.9–15 As an example θ ≈ 0.29 for the configuration studied in this paper,16 compared to KHI (∝t) and RTI (∝t2). Despite the relatively small temporal exponent, the growth rate of RMI is sufficiently high that it significantly impacts the physics of the aforementioned applications, so there is substantial interest in the accurate modeling of turbulent transport and mixing in RMI.

In the context of current computational power, turbulent transport models for unsteady Reynolds-Averaged Navier-Stokes (URANS) or Large-Eddy Simulation (LES) are particularly important. URANS methods are in active development and use for design problems, principally due to the complexity of applied computations typically impacted by multiple physical phenomena, complex equations of state, and extreme conditions. Turbulence resolving computations such as Large-Eddy Simulations require resolution of the full dimensionality of the problem,3 whereas for rapid design evaluation and iteration, it is desirable to work with a problem reduced to one or two dimensions by appropriate averaging and application of modeling of the full three-dimensional effects. This is the key motivation for the development and use of URANS modeling.

Modeling approaches for decaying, inhomogeneous, and variable density turbulent flows where density variations do not satisfy the Boussinesq approximation are particularly complex. The majority of turbulence modeling applied to variable density flows employs a density-weighted or Favre-averaged approach.17 This simplifies the resulting equations describing the evolution of the density-weighted mean compared to Reynolds-averaging. Following the Favre decomposition, two principal additional complexities arise (i) additional terms associated with density fluctuations and (ii) a change in the definition of the “mean” evolved, perhaps necessitating a modification to the closures of terms from the Reynolds averaged formulation.18 URANS modeling approaches applicable to RMI range from single-fluid two-, three-, and four-equation19–24 and Reynolds stress25–28 to two fluid models.29–31 There is a strong need for data to validate and further advance these and other models.

Despite significant advances in facilities and diagnostics,11,32–40 experimental data are very challenging to acquire due to the relatively short time scale of RMI in terrestrial facilities, extreme conditions, or uncertainty in the initial perturbations. Such experimental research is critical and permits important validation of individual terms; however, to date, it is insufficient to fully determine the required model coefficients for URANS computations. While there is considerable literature examining turbulent transport in shear layers induced by KHI, and some prior studies of compressible RTI,41 to date, there has not been an equivalent study of all individual terms in the turbulent transport equations computed simultaneously for a turbulent mixing layer developing from a RMI. Thus, the objective of the present study is to address a gap in the identification of the key mechanisms responsible for variable density effects in the RMI and the terms that must be modeled to accurately represent mean transport.

Recently, a benchmark multimode, three-dimensional RMI case was computed with multiple independent computational fluid dynamics codes, establishing a well understood “θ-group” dataset.16 The principal goal of that study was to understand and quantify numerical uncertainty in the computation of such flows and to provide a foundation for basic quantities of interest. The present study computes the transport equation budgets for the mean momentum, mean heavy fluid mass fraction, heavy fluid mass fraction variance, and specific turbulent kinetic energy from the results of the four codes presented in that paper: Triclade, Turmoil, Flamenco, and Flash. The analysis of these data identifies the terms of greatest importance in the mean transport equations and provides a reliable benchmark for LES and URANS model development and validation.

This paper is organized as follows: Sec. II summarizes the numerical methods employed to generate the θ-group dataset, the test case examined, the transport equations examined, and numerical methods employed to compute the terms. Section III presents and analyses the turbulent transport budgets for four representative times ranging from just postshock to a weakly turbulent state. Finally, the key conclusions are summarized in Sec. IV.

This paper utilizes the “quarter scale” data from four algorithms presented in the θ-group study of Thornber et al.16 In that paper, the full details of the initial conditions, results, and grid convergence studies for each individual algorithm are presented. It established an initialization to explore the physics of narrowband Richtmyer-Meshkov instability formulated in a way, which could be simulated by a wide range of algorithms. The configuration consisted of a light and a heavy gas, with a density ratio of 3, separated by a perturbed interface. A shock wave of strength Mach 1.84 impinges on the perturbed layer, triggering Richtmyer-Meshkov instability. The perturbation on the interface is defined using a narrowband power spectrum with constant amplitude between a defined minimum and maximum wavenumber, where modal amplitudes and phases are defined using deterministic random numbers, thus generating an exactly reproducible initial condition.

The initial conditions in the shocked heavy fluid domain 0.0 m < x < 3.0 m are (ρ, u, p) = (6.375 kg/m3, −61.49 m/s, 400 kPa). Given a surface perturbation A(y, z), the unshocked heavy fluid lies between 3.0 m < x < 3.5 m + A(y, z) with initial properties (ρ, u, p) = (3.0 kg/m3, −291.6 m/s, 100 kPa). The third region contains the unshocked light fluid that lies between 3.5 m + A(y, z) < x < 2.8π m and is initialized with properties (ρ, u, p) = (1.0 kg/m3, −291.6 m/s, 100 kPa). The configuration is shown schematically in Fig. 1.

FIG. 1.

Schematic of the initial condition, where the arrows indicate the initial velocity field.

FIG. 1.

Schematic of the initial condition, where the arrows indicate the initial velocity field.

Close modal

The perturbation power spectrum is constant between length scales of L/32 (λmin) and L/16, where L is the cross section, and an initial diffuse interface of error function form and thickness δ = L/128 is applied as follows:

f 1 ( y , z ) = 1 2    e r f c π [ x S ( y , z ) ] δ ,
(1)

where f1 is the volume fraction of the heavy gas and S(y, z) = 3.5 + A(y, z). The standard deviation of the perturbation amplitude is σ0 = 0.1λmin which places the shock-interface interaction initially at the large-amplitude end of the linear regime at the highest wavenumber. The diffuse initial condition ensures that the initial condition was accessible for all numerical schemes, as some algorithms are unstable in the presence of sharp interfaces. Together, the addition of a diffuse interface and the larger amplitude will impact the initial growth rates, particularly for the highest wavenumbers, which could be expected to be ≈20% lower than that predicted from linear theory.42 Although diffusion is not explicitly represented in the computation, an implicit assumption is that the gases are miscible. Turbulent layers evolving in immiscible fluids may have different growth rates and mixing parameters (see, for example, Refs. 43–45).

The heavy and light unshocked gases are at the same temperature in the undisturbed flow, and the specific heats for each component (Cv1 and Cv2) are chosen to ensure this. Finally, the problem description is completed by assuming an ideal gas equation of state with a ratio of specific heats γ = 5/3.

The Cartesian computational domain is Lx × Ly × Lz = 2.8π × 2π × 2π m3. The simulations were run from t = 0 to 5 s, which is equivalent to a dimensionless time τ = 246 (the nondimensionalization is detailed in Sec. II C). However, at times later than 2 s, individual spikes were observed to exit the domain; thus, only data up to 2 s are typically employed. Prior analysis has indicated that by the final time an asymptotic approximately self-similar state is approached—this assumption is examined in more detail here. Four time instants are examined, where the first time is immediately following shock passage (0.01 s), the second is shortly after this time where layer mixedness as measured by the mix parameter Θ has approximately doubled (0.025 s), and the third and fourth are approaching a self-similar state (1.0 and 2.0 s). See the discussion in Sec. III A and Fig. 4 for more details.

The configuration was chosen such that the low order statistics for the problem are just converged on the finest grid level using modern compressible algorithms. This introduces a degree of uncertainty on the influence of numerical dissipation on each transport term. This has been addressed in two ways: the first is through a grid convergence study and the second is by computing all transport terms using data produced from four independent codes/algorithms: Flamenco (Godunov method with the fifth order low-Mach-number-corrected MUSCL scheme46,47), Flash (compressible Euler with the Piecewise Parabolic Method and Monotonized Central limiter48–50), Triclade (conservative finite difference using the wave propagation algorithm of LeVeque51 with high order accurate corrections52), and Turmoil (Lagrange remap method with high-order artificial viscosity53–55). The boundary conditions are periodic in the homogeneous directions, but the best practice is adopted for each code in the shock-direction, which permits the transmission of the reflected rarefaction and transmitted shock through the boundaries with minimal reflections. The implicit numerical dissipation mechanisms in each of the algorithms differ substantially. Thus, agreement between the results either implies that for a given quantity (i) numerical dissipation has been minimized or (ii) numerical dissipation has acted to dissipate the quantity in a statistically similar manner. All simulations have no model for physical viscosity, diffusivity, or conduction.

This study explores the individual transport terms in the mean momentum and mass fraction, mass fraction variance, and fluctuating kinetic energy equations. A detailed prior study of Rayleigh-Taylor mixing of Schilling and Mueschke41 presented the required transport equations which are considered here. Referring to the equation numbering in Schilling and Mueschke,41 this paper examines the transport of mean momentum, Eq. (6), mean mass fraction, Eqs. (7) and (8), turbulent kinetic energy, Eqs. (14)–(15e), and mass fraction variance, Eqs. (19)–(20c). Averaging processes are defined in Eqs. (2) and (3) and accompanying text. In this contribution, the key differences are that the homogeneous directions are y and z, and the shock-direction is x, with a velocity component u, different from the Rayleigh-Taylor case where the principal direction and velocity are labeled z and w, respectively.

The Reynolds average of a field ϕ x , t over the statistically homogeneous plane is

ϕ ¯ ( x , t ) = 1 L y L z 0 L y 0 L z ϕ ( x , t ) d z d y ,
(2)

which enables a decomposition into mean and fluctuating components, where ϕ ( x , t ) = ϕ ¯ ( x , t ) + ϕ ( x , t ) and the fluctuating component is denoted by a prime. The Favre average is

ϕ ̃ ( x , t ) = ρ ϕ ¯ ( x , t ) ρ ¯ ( x , t ) = 0 L y 0 L z ρ ( x , t ) ϕ ( x , t ) d z d y 0 L y 0 L z ρ ( x , t ) d z d y ,
(3)

where ρ(x, t) is the density and ϕ ̃ ( x , t ) is defined such that ϕ ( x , t ) = ϕ ̃ ( x , t ) + ϕ ( x , t ) . This averaging reduces the fields to being one-dimensional functions of x only.

The individual terms which may be extracted from a single ILES are listed below, along with the notation used in the figures and analysis. For mean momentum transport,

ρ ¯ ũ t + u t r a n s ũ x t u + ρ ¯ ũ u t r a n s ũ x A u = p ¯ x F u τ 11 x R u + N D u ,
(4)

where τ 11 = ρ ¯ u 2 ̃ , u is the x-direction velocity component, p is the pressure, and NDu represents the effective numerical dissipation or diffusion for this equation. For mean mass fraction transport,

ρ ¯ m ̃ 1 t + u t r a n s m ̃ 1 x t m 1 + ρ ¯ ũ u t r a n s m ̃ 1 x A m 1 = ρ ¯ m 1 u ̃ x T m 1 + N D m 1 ,
(5)

where m1 is the mass fraction of the heavy fluid. For turbulent kinetic energy transport,

ρ ¯ E ̃ t + u t r a n s E ̃ x t E + ρ ¯ ũ u t r a n s E ̃ x A E = u ¯ p ¯ x P b E ρ ¯ u i u ̃ u i ̃ x P s E + p u k x k ¯ Π E x ρ ¯ E u ̃ + p u ¯ T E + N D E ,
(6)

where kinetic energy is E = u i u i / 2 . For mass fraction variance transport,

ρ ¯ m 1 2 ̃ t + u t r a n s m 1 2 ̃ x t m 1 2 + ρ ¯ ũ u t r a n s m 1 2 ̃ x A m 1 2 = 2 ρ ¯ m 1 u ̃ m ̃ 1 x P m 1 2 ρ ¯ m 1 2 u ̃ x T m 1 2 + N D m 1 .
(7)

The mean translational velocity utrans of the layer

u t r a n s ( t ) = 0 L x ū m 1 ̃ ( 1 m 1 ̃ ) d x 0 L x m 1 ̃ ( 1 m 1 ̃ ) d x
(8)

accounts for the small nonzero net velocity of the mixing layer in some simulations, for example, in the quarter scale case of Flamenco, there is a residual ũ ≈ 0.5 m/s attributed to imperfect boundary conditions (shock wave evacuation outside the computational domain). Each term is computed principally via sixth order central differences, but with biased stencils near the boundaries.

The nondimensionalization follows the analysis outlined by Thornber et al.16 These include estimates of the initial growth rate W ̇ 0 using a power-weighted mean wavelength λ ¯ = 2 π / k ¯ and mean postshock density ρ c + , which are defined as

W ̇ 0 = 0.564 7 12 k m a x A t + σ + Δ u , k ¯ = 7 12    k m a x ,     ρ c + = ρ 1 + + ρ 2 + 2 ,
(9)

where the postshock standard deviation of the perturbation amplitude is σ 0 + = C σ 0 , an approximate compression factor1 is C = (1 − Δu/Ui), where Ui = 434 m/s is the velocity of the incident shock, Δu = 291.575 m/s is the velocity impulse imparted by the shock; postshock densities were 5.22 and 1.80 kg/m3 for the heavy and light fluids, respectively, giving ρ c + = 3.51 kg/m3, uc = W ̇ 0 = 12.649 m/s, and λ ¯ = 0.2571 m. The layer center is defined as the location at which m ̃ 1 ( x , t ) = 0.5 .

Three measures of self-similarity may be employed. The first is to examine whether the layer properties at multiple time instants collapse when normalized by characteristic length, time, and mass scale. To examine this, some figures are plotted using axes scaled with the integral width

W ( t ) = 0 L x f ¯ 1 ( x , t ) [ 1 f ¯ 1 ( x , t ) ] d x ,
(10)

and W ̇ (t), the values of which are given in Table I. This table also gives a second measure, log 2 ( 3 W / λ ¯ ) , which is an indicator of the number of modal generations. Elbaz and Shvarts proposed that a mixing layer must achieve ≥3 generations to reach self-similarity.10 Here, approximately 2.29 generations are reached, an indication that the layer may not be in the fully self-similar regime. A third measure of self-similarity is an integral mixing measure such as the “molecular mixing”

Θ = f 1 f 2 ¯ d x f ¯ 1 f ¯ 2 d x ,
(11)

which varies by less than 1.5% in the last 3 s of the Flamenco simulation. The first measure will be explored further within this paper. For the rest of the paper, quantities are plotted against x/W, where x = 0 is the center of the mixing layer.

TABLE I.

Mixing layer width W and growth rate W ̇ as a function of time for the code-averaged data for the quarter scale problem.

t (s) τ W (m) log 2 ( 3 W / λ ¯ ) W ̇ (m/s)
0.01  0.49  0.076 05  −0.17  4.12 
0.025  1.23  0.110 7  0.37  1.436 
0.5  24.6  0.235  1.46  0.122 
49.7  0.283 3  1.73  0.077 9 
98.5  0.344 3  2.0  0.049 2 
197  0.419 5  2.29  0.029 56 
t (s) τ W (m) log 2 ( 3 W / λ ¯ ) W ̇ (m/s)
0.01  0.49  0.076 05  −0.17  4.12 
0.025  1.23  0.110 7  0.37  1.436 
0.5  24.6  0.235  1.46  0.122 
49.7  0.283 3  1.73  0.077 9 
98.5  0.344 3  2.0  0.049 2 
197  0.419 5  2.29  0.029 56 

Grid convergence of the mean flow quantities has been demonstrated by Thornber et al.16 for a highly resolved simulation where all length scales are four times larger (relative to the domain size). In the current study, higher order statistics are examined; thus, it is important to revisit grid convergence to ensure that these metrics are also reasonably converged. The  Appendix plots all of the individual terms of each transport equation for the code Flamenco at grid resolutions ranging from 180 × 1282 to 720 × 5122, in Figs. 27 and 28. The overall agreement at the two highest grid resolutions for all transport terms is excellent, indicating that the results are reasonably grid converged.

The spikes of heavy fluid are intense vortical structures of relatively small size (compared to the grid spacing) and as such are the most challenging “large” scale features to converge. Given even small differences in the trajectory at early times, strong interactions of the spike’s ring-shaped vortices with other spikes and the mixing layer can cause large differences at late times. This can be seen as a region of relatively lower convergence at x/W ≈ 8 for the kinetic energy and kinetic energy balances in Figs. 28(c) and 29(c). Such an intense impact of the spikes on the kinetic energy budget has also been observed in an analysis of the RMI spectral kinetic energy budget.56 

The levels of noise in each of the quantities are a good indicator of statistical convergence, as no smoothing has been applied. However, the statistical noise does not mask the overall trends for the majority of the results presented. Based on the previous study,16 it is expected that results from Turmoil and Triclade would converge at a similar grid resolution to Flamenco. Flash, however, has historically demonstrated slower convergence compared to the other algorithms, but given the good agreement with the other algorithms presented here, and the fact that Flamenco is reasonably converged at the intermediate grid level (a factor of eight times fewer points than the highest resolution), the Flash results are also expected to be reasonably converged.

Compiling the data from each code highlighted numerical and physical issues at different times. Flamenco exhibits postshock oscillations at the earliest times which are clear in the plots of ũ, Triclade and Turmoil show highly fluctuating acoustic pressure correlations at late times (terms such as ΠE), and Flash has difficulties at the spike side boundary at late times (visible in the spike side E ̃ / u c 2 at t = 2 s, for example). An effort has been made here to present all results except where (i) the result is so noisy that it obscures interpretation of the rest of the data or (ii) it is clearly dominated by a nonphysical effect. There are some cases which are marginal (perhaps unphysical) but can fit on the plot without obscuring the data. Where data are deliberately excluded due to such effects, the omission and reason for omission are stated in the figure caption.

Visualization of the distribution of the volume fractions within the flow field is shown for the four time instants chosen using slices in the y-plane in Fig. 2 and three-dimensional isosurfaces of volume fraction in Fig. 3 from the Flamenco code.

FIG. 2.

Visualization of slices of volume fraction of the heavy fluid from the Flamenco simulation of the mixing layer at t = 0.01, 0.025, 1, and 2 s (from left to right). Blue indicates the heavy fluid and red the light fluid.

FIG. 2.

Visualization of slices of volume fraction of the heavy fluid from the Flamenco simulation of the mixing layer at t = 0.01, 0.025, 1, and 2 s (from left to right). Blue indicates the heavy fluid and red the light fluid.

Close modal
FIG. 3.

Visualization of volume fractions from the Flamenco simulation of the mixing layer at t = 0.01, 0.025, 1, and 2 s, looking at the spike side of the mixing layer. Red indicates the isosurface of heavy fluid volume fraction 0.001 and blue 0.999.

FIG. 3.

Visualization of volume fractions from the Flamenco simulation of the mixing layer at t = 0.01, 0.025, 1, and 2 s, looking at the spike side of the mixing layer. Red indicates the isosurface of heavy fluid volume fraction 0.001 and blue 0.999.

Close modal

At the initial time t = 0.01 s, the individual modes at the highest wavenumber have already formed characteristic mushroom shaped bubbles (light fluid penetrating the heavy) and spikes (heavy fluid penetrating the light). As expected at this Atwood number, several spikes have propagated substantially beyond the bulk of the mixing layer. Two of these are visible on the left-most image in Fig. 2; however, due to the perspective, they are not as clearly seen in the three-dimensional rendering.

At t = 0.025 s, the layer has begun to transition, where the spikes contain relatively well mixed fluid, but the bubble heads still contain some unmixed material. This provides a preliminary indication that turbulent fluctuations are initially stronger on the spike side of the mixing layer. There is visible evidence of the vortex projectiles propagating far from the layer on the spike side; however, their intensity and size have reduced. At the two latest times, there is a core of well mixed fluid, but the flow is still relatively nonuniform in the homogeneous direction. The layer width grows by ≈20% between t = 1 and 2 s, and low wavenumber modes dominate at t = 2 s. Strong three-dimensionality and a range of vortex length-scales are apparent at t = 1 and 2 s in Fig. 3; however, there is still evidence of coherent spike structures such as the ring vortex at the upper left of the three-dimensional visualization at t = 2 s.

Figure 4 plots the integral width and molecular mixing measure Θ at early and late times. Superimposed on this figure are the chosen sampling times. The first two sampling times are chosen to be just after the shock passes through the interface and at a short time later. At this stage, the layer is expected to be formed of coherent structures and a wide range of turbulent length scales is not expected. The two latest times are chosen to be at a time where the layer may be approaching self-similarity, but with a compromise between having the widest possible range of turbulent length scales and a moderate impact of a lack of statistical accuracy in the homogeneous direction and domain-size constraint.57 From t = 1 to 4 s, Θ changes by only 1.5% in Flamenco simulations, an indication that it is approaching a self-similar state (but not yet reached, as discussed earlier).

FIG. 4.

Integral width and molecular mixing measure Θ presented for (a) early and (b) late times, illustrating the times considered with a solid filled circle.

FIG. 4.

Integral width and molecular mixing measure Θ presented for (a) early and (b) late times, illustrating the times considered with a solid filled circle.

Close modal

Figure 5 plots the profiles of ũ, , and w ̃ normalized by uc immediately following shock interaction (0.01 s), a short time later where mixing is increasing rapidly (0.025 s), at the onset of approximate self-similarity (1.0 s), and at the latest time (2.0 s). The overall agreement for ũ is very good across the codes, particularly, at late times, with the main discrepancies at early times where the peak ũ varies by 14% from the average peak value. The Favre-averaged velocity is positive, representing the net flux of heavy material into the light material as expected. At early times, the profile of ũ is asymmetric, with a larger gradient on the bubble side (x < 0) and more gradual decrease on the spike side (x > 0); however, at late times, it is reasonably symmetric and the peak velocity is at the centerline (defined by mean mass fraction of one half). Asymmetries are expected at this moderate Atwood number.

FIG. 5.

Mean velocity components ũ, , and w ̃ at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s. All terms are nondimensionalized by uc.

FIG. 5.

Mean velocity components ũ, , and w ̃ at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s. All terms are nondimensionalized by uc.

Close modal

The Favre-averaged in-plane velocities and w ̃ are relatively large compared with ũ at late times, where it could be expected that these quantities should average to zero in the limit of an infinite domain. As the domain size is finite, these components may take on values roughly equivalent to the noise observed in the evaluation of ũ. This is clearly not the case as the profile of ũ is relatively smooth. The large magnitude of and w ̃ also extends well into the spikes (x/W > 4), where the plane-averaged heavy mass fraction is <1%, and these fluctuations are associated with the propagation of “vortex projectiles” exiting the mixing layer. As this is a turbulent mixing layer, it is not surprising to see differences in the locations of the peaks and troughs for the in-plane components between the algorithms, especially given that the mean and w ̃ are 0.6% of the peak values at a given x location.

The relatively large magnitude is most likely due to the mean gradient of ũ in the x direction being constrained to be close to zero due to approximate incompressibility and the corresponding divergence-free constraint. Due to numerical miscibility, the current simulations will not satisfy ∇ · u = 0, as ∇ · u = −D ln ρ/Dt, where D(·)/Dt is the material derivative. Although the flow at late times is at low Mach number and nearly incompressible, the advected density changes due to numerical mixing. The periodic boundary conditions in the y and z directions do not constrain the net flow in that direction, thus leading to the observed profiles.

Late time self-similarity is examined in Fig. 6, where the Flamenco data have been nondimensionalized by the integral width W(t) and the time rate of change of the integral width W ̇ (t). An additional point at t = 4 s has been added, where at that time, spikes have left the computational domain, but the bulk of the mixing layer may be assumed to be reasonably well resolved. The momentum profiles collapse well for the bulk of the layer under the chosen scalings for t > 1 s; however, the spike side is taking longer to achieve approximate self-similarity.

FIG. 6.

Mean velocity component ũ at t = 0.5, 1, 2, and 4 s (Flamenco results only) nondimensionalized by W(t) and W ̇ (t).

FIG. 6.

Mean velocity component ũ at t = 0.5, 1, 2, and 4 s (Flamenco results only) nondimensionalized by W(t) and W ̇ (t).

Close modal

Figure 7 shows the individual components of the mean momentum transport equation for all codes, Fig. 8 shows the same time instants but only the Flamenco results so that the relative magnitudes of the individual terms may be more clearly discerned, and Fig. 9 plots the Fu term alone, where the data have been smoothed by a simple five point moving average. The codes agree well for all terms except for the pressure gradient term Fu, which shows good agreement on the bubble side of the mixing layer but relatively poor agreement on the spike side. Although the bulk of the mixing layer is contained in −4 ≤ x/W ≤ 4, there is a substantial contribution to the mean momentum equation at x/W > 4 which may be attributed to highly energetic ejected spikes and features which combine both very high velocities and low pressure cores in a localized feature. Such features have been observed in previous analysis of spectral energy transfer in the RMI58 and RTI.59 

FIG. 7.

Individual terms in the ũ transport equation at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s. Flamenco time derivative at t = 0.025 s is not shown as it is dominated by the dissipation of the oscillations (see Fig. 8), and Turmoil Fu at t = 0.025, 1, and 2 s are not shown due to large fluctuations which mask the other curves (see Fig. 9 for smoothed Fu). All terms are nondimensionalized by ρ c u c 2 / λ ¯ .

FIG. 7.

Individual terms in the ũ transport equation at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s. Flamenco time derivative at t = 0.025 s is not shown as it is dominated by the dissipation of the oscillations (see Fig. 8), and Turmoil Fu at t = 0.025, 1, and 2 s are not shown due to large fluctuations which mask the other curves (see Fig. 9 for smoothed Fu). All terms are nondimensionalized by ρ c u c 2 / λ ¯ .

Close modal
FIG. 8.

Individual terms in the ũ transport equation at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s computed from the Flamenco results. All terms are nondimensionalized by ρ c u c 2 / λ ¯ .

FIG. 8.

Individual terms in the ũ transport equation at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s computed from the Flamenco results. All terms are nondimensionalized by ρ c u c 2 / λ ¯ .

Close modal
FIG. 9.

Smoothed Fu term at (a) t = 1 s and (b) 2 s. All terms are nondimensionalized by ρ c u c 2 / λ ¯ .

FIG. 9.

Smoothed Fu term at (a) t = 1 s and (b) 2 s. All terms are nondimensionalized by ρ c u c 2 / λ ¯ .

Close modal

As has been shown for a small Atwood number Rayleigh-Taylor instability,15 the time rate of change of mean momentum tu is mostly generated as the remainder of two larger terms, the mean pressure gradient Fu and the Reynolds stresses Ru, where Fu ≈ −Ru at all times. At early times, both Ru and Fu are approximately antisymmetric across the centerline, but at later times, the bubble side amplitudes are larger and closer to the core of the mixing layer than the spike side amplitudes, as the gradient of the Reynolds stress τ11 is much sharper on the bubble than the spike side (see Sec. III E). Mean advection is nearly negligible at all times, which is expected as the mean layer position is approximately stationary with respect to the chosen reference frame.

The mean momentum transport equation has contributions from numerical dissipation NDu representing the under-resolved components of the Reynolds stresses. As the grid resolution increases, the magnitude of these implicitly modeled effects should decrease. Figure 10 plots the nondimensional numerical dissipation term for Flash, Flamenco, and Triclade, computed as the residual of the mean momentum equation. Although the overall magnitude is similar to the other terms in the momentum equation, the magnitude is largely due to the presence of noise in the computation of the individual physical terms as shown clearly in Fig. 7. Examining only the Flamenco results (which have relatively low noise) shows that the peak magnitude is less than 3.5 × 10−5, which is ten times lower than the dominant physical terms, and the rms magnitude defined over −5 < x/W < 5 is 7.6 × 10−6. A similar conclusion can be reached if a simple moving average is applied to the results. It should be noted however that as the net change of mean momentum is a result of the difference between the two largest physical terms Fu and Ru, the numerical dissipation represents approximately one quarter of the magnitude of the overall rate of change of ũ and is not negligible.

FIG. 10.

Nondimensional numerical dissipation estimated by the residual in the mean momentum transport balance at t = 1 s.

FIG. 10.

Nondimensional numerical dissipation estimated by the residual in the mean momentum transport balance at t = 1 s.

Close modal

Figure 11 presents the profile of the mean heavy fluid mass fraction m ̃ 1 . The profile is approximately linear and at late times exhibits no discernable kinks indicating a region of well mixed flow, with very small mean diffusive transport. All codes are in good agreement barring some late time differences on the spike side. Figure 12 shows the scaled m ̃ 1 profile from Flamenco at four times, demonstrating an excellent collapse for the three latest times on the bubble side in particular. Figure 12 also shows a close-up of the spike side, where the size of the spikes relative to the mixing layer is still decreasing at late times, a further indication that a fully self-similar state is not yet attained. This is consistent with breakup and dissipation of these strong vortical structures.

FIG. 11.

Mean mass fraction m 1 ̃ at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s.

FIG. 11.

Mean mass fraction m 1 ̃ at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s.

Close modal
FIG. 12.

Mean mass fraction m 1 ̃ at t = 0.5, 1, 2, and 4 s (right), and a close up of the spike side (left) (Flamenco results only).

FIG. 12.

Mean mass fraction m 1 ̃ at t = 0.5, 1, 2, and 4 s (right), and a close up of the spike side (left) (Flamenco results only).

Close modal

Figure 13 plots all terms in the m ̃ 1 transport equation, where all code results collapse well. Turbulent transport Tm1 is approximately antisymmetric about the center of the mixing layer and dominates for x/W < 0; however, due to mean advection Am1 being negative throughout the layer, the time rate of change tm1 dominates for x/W > 0. Mean advection Am1 is non-negligible, symmetric at the two latest times, and negative throughout the layer due to the net transport from the heavy to light side, but is less than half the magnitude of turbulent transport. As a result, the time variation of m ̃ 1 largely follows Tm1 but is further skewed toward the heavy side due to the contribution from Am1.

FIG. 13.

Individual terms in the m 1 ̃ transport equation at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s. All terms are nondimensionalized by ρ c u c / λ ¯ .

FIG. 13.

Individual terms in the m 1 ̃ transport equation at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s. All terms are nondimensionalized by ρ c u c / λ ¯ .

Close modal

Similar to the mean momentum equation, the results presented here also include numerical dissipation NDm1, which is estimated as the remainder in the governing equations and shown in Fig. 14 for t = 1 s. Again, the result is quite noisy, and the peak of 9 × 10−5 is at least one order of magnitude lower than the peak of the dominant physical terms. There is no discernable structure to the numerical dissipation profile, exhibiting as symmetric noise about zero.

FIG. 14.

Nondimensional numerical dissipation estimated by the residual in the m ̃ 1 turbulent transport balance at t = 1 s.

FIG. 14.

Nondimensional numerical dissipation estimated by the residual in the m ̃ 1 turbulent transport balance at t = 1 s.

Close modal

Profiles of the heavy species mass fraction variance m 1 2 ̃ are shown in Fig. 15 for four representative times. At early times, the profiles peak moderately on the spike side of the mixing layer and are relatively symmetric. There is a substantial difference in the peak mass fraction variance from each code, bounded by Flamenco as the lowest and Turmoil as the highest. This is expected as at early times the layer is in a state of transition, where the initially smooth density variation across the interface is being sharpened to a grid-scale discontinuity by shear flows driven by the growth of RMI. The minimum thickness of the layer is a function of the dissipation within each code, hence the increased uncertainty. At the latest two times, the mean gradients within the layer are smoother and the codes are in better agreement, although the peaks are still only in agreement to within ±10%. As time progresses, the variance decreases toward a self-similar profile as stirring reduces the gradients through mixing.

FIG. 15.

Mass fraction variance m 1 2 ̃ at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s.

FIG. 15.

Mass fraction variance m 1 2 ̃ at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s.

Close modal

Figure 16 shows the collapse of the m 1 2 ̃ profiles at several times. Except for the earliest “transitional” time of 0.5 s, the profiles at the other three times collapse well, having a similar peak value of ≈0.046. Mass fraction variance gradients are steeper on the spike side of the layer compared to the bubble side.

FIG. 16.

Mass fraction variance m 1 2 ̃ at t = 0.5, 1, 2, and 4 s (Flamenco results only).

FIG. 16.

Mass fraction variance m 1 2 ̃ at t = 0.5, 1, 2, and 4 s (Flamenco results only).

Close modal

Individual terms in the transport equation for m 1 ̃ are plotted in Fig. 17 for all codes. At the two earliest times, all terms but mean advection are important. At the latest time, mean production dominates all other physically resolved terms by a factor of three, although it must be noted that mean numerical dissipation is of a similar magnitude. Turbulent transport increases variance at the bubble and spike fronts but reduces it within the core. Although this also occurs at the later times, it is dominated by mean production, which is asymmetrically weighted toward the bubble side of the layer.

FIG. 17.

Individual terms in the m 1 2 ̃ transport equation at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s. All terms are nondimensionalized by ρ c u c / λ ¯ . Flash time derivatives are not shown at 1 and 2 s as they follow the same trends as Flamenco and Triclade but contain substantially higher noise, which otherwise masks the underlying trend.

FIG. 17.

Individual terms in the m 1 2 ̃ transport equation at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s. All terms are nondimensionalized by ρ c u c / λ ¯ . Flash time derivatives are not shown at 1 and 2 s as they follow the same trends as Flamenco and Triclade but contain substantially higher noise, which otherwise masks the underlying trend.

Close modal

The numerical dissipation NDm1″ has been computed for t = 1 s and is plotted in Fig. 18 for Flamenco, Flash, and Triclade, alongside the production over dissipation for Flamenco. The agreement between the codes is excellent, particularly, given the different algorithmic approaches. As described above, numerical dissipation is of similar magnitude to the largest physical terms and is purely dissipative. The ratio of production to dissipation peaks at ≈1.5, somewhat larger than the Rayleigh-Taylor instability value of 1.25 reported by Schilling and Mueschke41 and takes a volume weighted average of 0.92 over the region −5 < x/W < 5.

FIG. 18.

Nondimensional numerical dissipation estimated by the residual in the mass fraction variance transport balance at t = 1 s (left) and the production rate over dissipation rate (right).

FIG. 18.

Nondimensional numerical dissipation estimated by the residual in the mass fraction variance transport balance at t = 1 s (left) and the production rate over dissipation rate (right).

Close modal

The dimensionless turbulent kinetic energy E ̃ is shown in Fig. 19 for four times. At the earliest time, all codes show a “double peak” structure, where the peaks are focused on the bubble and spike heads, respectively. The highly energetic spikes persist at late times, causing a broadening of the profile into the light fluid (x/W > 0). The bulk of the turbulent kinetic energy content in the bubbles is restricted to the region x/W > −6 and increases steeply on the bubble side. This is consistent with the visualizations in Fig. 2, which shows an earlier onset of mixing on the spike side (an indicator of stirring by a range of length scales) than on the bubble side. The advection of spikes away from the layer at early times smooths the spatial contribution of the spike to the turbulent kinetic energy profiles, and turbulent transport acts to merge the early “double bump” profile into a single peak. The peak is slightly biased to the spike side of the mixing layer at late times, where most of the kinetic energy is located. For example, at t = 1 s and x/W = 5.5, there is substantial E ̃ even though m ̃ 1 < 0.1 %. At late times, the profile is asymmetric. In comparing the codes, there are consistent peak positions in Triclade, Turmoil, and Flamenco. Flash shows large values on the spike side at the latest time, which is attributed to unphysical interactions with the boundaries.

FIG. 19.

Normalized turbulent kinetic energy E ̃ / u c 2 at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s.

FIG. 19.

Normalized turbulent kinetic energy E ̃ / u c 2 at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s.

Close modal

The scaling of the profiles of E ̃ with the mixing layer width and its growth rate is plotted in Fig. 20. Overall, the collapse is excellent for the core of the layer, but notable differences can be seen in the prediction of the peaks in E ̃ in the spikes. As observed previously, the approximate vortex ring structures do not scale with the same power law as the rest of the mixing layer and generate isolated peaks and troughs in the E ̃ profile which are not consistent from one time to the next when scaled using the integral layer properties. However, the overall shape of the profile is in good agreement given the factor of 8 reduction in dimensional kinetic energy between the first time and the last shown in Fig. 20.

FIG. 20.

Normalized turbulent kinetic energy E ̃ / u c 2 at t = 0.5 s, 1, 2, and 4 s (Flamenco results only).

FIG. 20.

Normalized turbulent kinetic energy E ̃ / u c 2 at t = 0.5 s, 1, 2, and 4 s (Flamenco results only).

Close modal

There are six terms in the kinetic energy budget, of which the time variation tE, pressure-dilatation correlation ΠE, and turbulent transport TE terms are relatively large, and mean advection AE, buoyancy production P b E , and shear production P s E are relatively small. The large terms are plotted together in Fig. 21 and the small terms in Fig. 22. To aid a clear interpretation of the relative magnitude of each of the terms, the full budget is plotted for the Flamenco results only in Fig. 23. Finally, P b E and TE have strong fluctuations; thus, Fig. 24 plots those two terms independently at t = 1 and 2 s where a five point moving average has been applied to smooth the data.

FIG. 21.

Large magnitude individual terms in the E ̃ transport equation at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s. All terms are nondimensionalized by ρ c u c 3 / λ ¯ . Note that t = 0.025 s does not include TE for Turmoil and t = 1 s does not include TE for Turmoil and ΠE for Triclade. Turmoil data are not shown as fluctuations mask the other lines (see Fig. 24 for full data for TE).

FIG. 21.

Large magnitude individual terms in the E ̃ transport equation at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s. All terms are nondimensionalized by ρ c u c 3 / λ ¯ . Note that t = 0.025 s does not include TE for Turmoil and t = 1 s does not include TE for Turmoil and ΠE for Triclade. Turmoil data are not shown as fluctuations mask the other lines (see Fig. 24 for full data for TE).

Close modal
FIG. 22.

Small magnitude individual terms in the E ̃ transport equation at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s. All terms are nondimensionalized by ρ c u c 3 / λ ¯ . Note that the t = 0.01, t = 0.025, and 2 s does not include Turmoil P b E , 0.025 s also does not include Turmoil P s E , and t = 2 s does not show the Flash P b E term. Turmoil data are not shown as fluctuations mask the other curves (see Fig. 24 for full data for P b E ).

FIG. 22.

Small magnitude individual terms in the E ̃ transport equation at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s. All terms are nondimensionalized by ρ c u c 3 / λ ¯ . Note that the t = 0.01, t = 0.025, and 2 s does not include Turmoil P b E , 0.025 s also does not include Turmoil P s E , and t = 2 s does not show the Flash P b E term. Turmoil data are not shown as fluctuations mask the other curves (see Fig. 24 for full data for P b E ).

Close modal
FIG. 23.

Individual terms in the E ̃ transport equation at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s computed from the Flamenco results. All terms are nondimensionalized by ρ c u c 3 / λ ¯ .

FIG. 23.

Individual terms in the E ̃ transport equation at (a) t = 0.01 s, (b) 0.025 s, (c) 1 s, and (d) 2 s computed from the Flamenco results. All terms are nondimensionalized by ρ c u c 3 / λ ¯ .

Close modal
FIG. 24.

Smoothed P b E [(a) and (b)] and TE [(c) and (d)] plotted for all codes at 1 s [(a) and (c)] and 2 s [(b) and (d)]. All terms are nondimensionalized by ρ c u c 3 / λ ¯ .

FIG. 24.

Smoothed P b E [(a) and (b)] and TE [(c) and (d)] plotted for all codes at 1 s [(a) and (c)] and 2 s [(b) and (d)]. All terms are nondimensionalized by ρ c u c 3 / λ ¯ .

Close modal

First, a note on some numerical issues in computing specific terms. The overall level of fluctuations in the results is substantially larger than in the transport equations for m ̃ 1 and ũ, which is partly expected as E ̃ is a second order quantity; however, there were also specific difficulties in the computation of terms involving correlations of pressure fluctuations. This led to a large variation in results from code to code for the pressure-dilatation ΠE and turbulent transport TE terms.

Several different stencils were employed to compute the pressure-dilatation term, ranging from second to sixth order to explore the convergence of the numerical results. For Flamenco, it was noted that both the magnitude and sign changed for the pressure-dilatation correlation from second order estimation of derivatives to sixth order. It is encouraging that the third and fourth order stencils are converging toward the sixth order result.

Now, consider the individual terms in Fig. 21. At the first and second times, turbulent transport TE forms two symmetric peaks at the bubble and spike fronts, indicating a transport/redistribution of turbulent kinetic energy from the core of the layer to the edges. Turbulent transport is the largest term contributing to the time rate of change of E ̃ . The peak in TE at the bubble side remains at later time; however, the peak on the spike side advects away from the layer, becoming flatter and indicating a net transport of E ̃ from the spike to bubble side.

The pressure-dilatation term is symmetric, peaking at the center of the layer, and increases E ̃ . As pointed out by Schilling and Mueschke,41 although the Mach number is low, mean dilatation is non-negligible due to mixing of the fluids; thus, ∇ · u = −∇ · (D∇ ln ρ), with D being an effective numerical diffusion, the form of which may be more complex than Fickian.

Both Turmoil and Triclade results exhibited large fluctuations in pressure fluctuation-based terms which precluded the extraction of clear trends from their data. This is due to the treatment of acoustic waves in both of those algorithms, which both dissipate acoustic fluctuations less strongly than Flamenco and Flash. Acoustic waves add a “fast” component to the kinetic energy budget, but are a mostly isentropic process, and thus should not impact the overall time rate of change of E ̃ . Figure 19 shows that this is the case: the peak kinetic energies have reduced by a factor of 1000 over the simulation time, yet all codes remain in good agreement; thus, net dissipation must be of a similar order of magnitude for all algorithms and very likely to be scaling relative to a physical (rather than numerical) dissipation rate. This division of the pressure-dilatation into compressible and incompressible components is in agreement with previous observations for homogeneous decaying turbulence and shear layers,18,60,61 where it was noted that although the fast component has a large magnitude, its global impact is small compared with the “incompressible” correlation.

Figure 25 plots the spatial distribution of the rms pressure variance p 2 at several times using Flamenco and the temporal variation at the center of the mixing layer [defined by x such that m ̃ 1 ( x , t ) = 0.5 ] and the integral average in the heavy bulk fluid [x such that m ̃ 1 ( x , t ) > 0.99 ] and light fluid [x such that m ̃ 1 ( x , t ) < 0.01 ]. At late times, the pressure fluctuations within the mixing layer in Flamenco clearly follow an incompressible scaling ≈21.5ρav W ̇ 2 with ρ a v = ( ρ H + + ρ L + ) / 2 . This is a reasonable result noting that the integral width is substantially smaller than the actual width of the layer which leads to the large constant of proportionality. However, for Triclade and Turmoil, the pressure fluctuations are substantially larger.

FIG. 25.

The variation of (a) p 2 ( t ) as a function of time and space in Flamenco, and as a function of time in (b) the heavy fluid, (c) the light fluid, and (d) the layer center for Flamenco, Triclade, and Turmoil.

FIG. 25.

The variation of (a) p 2 ( t ) as a function of time and space in Flamenco, and as a function of time in (b) the heavy fluid, (c) the light fluid, and (d) the layer center for Flamenco, Triclade, and Turmoil.

Close modal

Focusing on the fluctuations in the pure heavy fluid in Fig. 25(b) and pure light fluid in Fig. 25(d), the fluctuations are modeled as a simple spherical source deposited at shock-interaction, where given conservation of acoustic energy, the pressure fluctuations should scale with p 2 ( t ) p 2 ( t 0 ) r 0 ( t 0 ) / [ r 0 ( t 0 ) + a H , L t ] . Here, aH and aL are the speed of sound in the shocked heavy/light fluid, respectively, t0 is a sampling time measured from shock interaction, and the initial acoustic radius is modeled as r 0 ( t 0 ) = λ ¯ / 2 + a H t 0 .

The overall agreement is very good from t = 0.01 to 0.1 s between Flamenco, Triclade, and Turmoil, despite the complexities of the acoustic waves and the influence of the exiting rarefaction at early times (t < 0.1 s). At later times, numerical dissipation in Flamenco reduces the magnitude of the acoustic waves in the heavy and light fluids, to below the pressure fluctuations produced by turbulent motion in the layer. This is not seen in Turmoil and Triclade, which preserve acoustic fluctuations better. One-dimensional tests indicate that Turmoil damps acoustic waves less strongly than Flamenco, which is consistent with this result, especially noting that by t = 1 s the wave has traveled through the cross section more than 400 times and so numerical dissipation has had ample opportunity to cause the observed differences. This explains the discrepancy between the three codes in the layer center, where pressure fluctuations in the heavy fluid penetrate the layer and are at a larger magnitude than pressure fluctuations generated through vortical motion alone.

At even later times, acoustic pressure fluctuations will dominate the vortical fluctuations by orders of magnitude. This also explains the difficulty in obtaining agreement among the codes for transport terms including pressure gradients, and correlations including pressure fluctuations and gradients. Despite this, the impact on the overall transport balance is small.

The time derivative tE is plotted in Fig. 21 and shows an increase in turbulent kinetic energy at both the bubble and spikes at the earliest times. However, at t = 1 and t = 2 s, the time rate of change is negative throughout the core of the mixing layer, with positive variations only at the isolated vortex “projectiles.” Note that the overall balance is impacted by dissipation provided by the implicit numerical model, which is discussed at the end of this section.

Examining the smaller terms plotted in Fig. 22, mean advection AE shows a transport of kinetic energy from the bubble to the spike side of the layer. Both buoyancy and shear production are negative on the bubble side and positive on the spike side, and all three terms are only nonzero within the core of the mixing layer. Some oscillations in P s E may be seen at x/W > 8 at t = 2 s for Flash, due to the unphysical interaction of the spikes with the boundary condition.

The relative magnitude of all six terms is shown clearly in Fig. 23 for the Flamenco results only. It is clear that the time evolution of E ̃ is dominated by two physical effects, the turbulent transport and the pressure-dilatation correlation, and by numerical dissipation with a mean magnitude of the same order as the resolved physical terms.

The numerical dissipation NDE has been computed for t = 1 s and is plotted in Fig. 26. This figure compares the numerical dissipation computed for Flamenco and Flash, where the Triclade results are not shown as they are dominated by the fluctuations in the pressure-dilatation term due to the acoustic field. As mentioned in the previous paragraph, the numerical dissipation of turbulent kinetic energy is of a similar order to the dominant terms in the transport balance. This large order of magnitude is to be expected. Turbulent kinetic energy is conserved until it reaches the viscous scales, where it is dissipated. In an implicit LES, turbulent kinetic energy is dissipated by numerical dissipation as it reaches the grid scale. The mechanism of turbulent dissipation is highly code-dependent, with some algorithms having dissipation in the remap phase for low Mach number flows (e.g., Turmoil), and others using upwind schemes to dissipate on the individual wave strengths (e.g., Flamenco, Flash, and Triclade).

FIG. 26.

Nondimensional numerical dissipation estimated by the residual in the mean kinetic energy transport balance at t = 1 s.

FIG. 26.

Nondimensional numerical dissipation estimated by the residual in the mean kinetic energy transport balance at t = 1 s.

Close modal

However, implicit LES relies on the concept that the actual dissipation rate is driven by the largest scales of turbulence, and independent of the detailed mechanism of dissipation as long as the dissipation does not impact the large scales. Thus, when a sufficiently wide range of large scales is resolved by a numerical algorithm, the flow would evolve in a statistically identical manner regardless of whether the smallest scales of the flow are resolved or not. The word “statistical” is used since turbulent flows are very sensitive to small perturbations; hence, an exact match is not expected.

That the four codes used in the current study have such good agreement in the turbulent kinetic energy profiles at a very late time, where peak values have reduced by a factor of 1000, indicates that the decay rates are very similar from algorithm to algorithm. Thus, although only results from two codes are plotted in Fig. 26, there is reasonable confidence that a similar magnitude would be observed for the other two codes. Finally, the grid convergence of NDE for Flamenco as shown in Fig. 27 also indicates that the simulations have achieved the necessary scale separation such that the total dissipation rate is determined by processes at resolved scales, rather than by numerical dissipation.

FIG. 27.

Convergence of nondimensional numerical dissipation for Flamenco estimated by the residual in the turbulent transport balances: (a) ũ, (b) m 1 ̃ , (c) E ̃ , and (d) m 1 ̃ .

FIG. 27.

Convergence of nondimensional numerical dissipation for Flamenco estimated by the residual in the turbulent transport balances: (a) ũ, (b) m 1 ̃ , (c) E ̃ , and (d) m 1 ̃ .

Close modal

An analysis of turbulent transport in an inhomogeneous compressible turbulent mixing layer induced by Richtmyer-Meshkov instability has been presented. The analysis used data from four independent ILES codes which were applied to the quarter-scale problem presented in the θ-group collaboration.16 All terms in the high-Reynolds number (inviscid) limit of the Reynolds-averaged transport equations for momentum, species mass fraction, turbulent kinetic energy, and mass fraction variance were computed in the inhomogeneous direction. Numerical dissipation provided an implicit subgrid model in all algorithms, which was computed from the balance of the transport equations. The objective of this study was to provide insight into the dominant physical mechanisms influencing turbulent transport and to facilitate the future development and validation of reduced order modeling of those quantities. The inclusion of four numerical methods also provided an understanding of the confidence in predicting each of the terms, and to compare and contrast the impact of numerical dissipation.

For the momentum transport equation, there is high confidence in the mean momentum in the turbulent stage, with codes agreeing to within ≈±4% in the shock-direction. The time rate of change of momentum is dominated by the balance of the pressure gradient and Reynolds stresses, and the numerical dissipation is very low. The pressure gradient term was particularly noisy due to acoustic modes. Advecting spikes which escape the mixing layer cause large variations in the individual transport terms, the position and magnitude of which were not easily captured numerically. The results here may also be impacted by the small sample size of such spikes which escape the developing mixing layer at an early time. The individual transport terms have a spread of ≈±10% due to statistical errors and acoustic fluctuations, but the agreement between codes is good and could be employed to validate URANS and LES modeling approaches.

For mass fraction transport, the mean profiles are nearly identical at all times; thus, there is a high level of confidence. The time variation of mass fraction distribution largely follows the turbulent transport term, moderated by mean advection, and numerical dissipation is small. Here, there is very good agreement between all algorithms, and thus a high level of confidence in the interpretation of the dominant terms. The mass fraction variance shows a higher variation at early times, but again agreement is very good when the layer transitions to turbulence (≈±5%). All terms in the transport equation are important, but mean advection is lower by a factor of approximately two. Numerical dissipation is purely negative as expected and is of a similar magnitude to the largest resolved terms. All terms considered are in sufficient agreement to be useful for validation of lower order modeling results.

The analysis of the turbulent kinetic energy transport equation was more challenging. The largest error in turbulent kinetic energy occurred at the earliest, nonturbulent times (≈±10%), but with excellent agreement at later times (≈±3%). Temporal evolution is dominated by turbulent transport and numerical dissipation (implicit subgrid dissipation), with pressure-dilatation the third largest. The computed pressure-dilatation term was strongly impacted by acoustic modes, rendering the confidence in the results for that term relatively low compared to the rest of the data, given the numerically challenging problem of evolving sound waves for relatively long periods. Despite this uncertainty, the impact on the overall development of kinetic energy in the layer is concluded to be low, as evidenced by the excellent agreement in late time kinetic energy for algorithms which have substantially different mean pressure fluctuations at late times.

The authors greatly acknowledge Dr. Karnig Mikaelian who suggested the original “θ-group” collaboration at the Fourteenth International Workshop on the Physics of Compressible Turbulent Mixing in 2014. This research was supported under the Australian Research Council’s Discovery Projects funding scheme (Project No. DP150101059). The lead author would like to acknowledge the computational resources at the National Computational Infrastructure through the National Computational Merit Allocation Scheme which were employed for all Flamenco cases presented here. Flash was developed by the DOE-sponsored ASC/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago. This work was performed under the auspices of the U.S. Department of Energy by the Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344.

Due to the large number of individual terms computed, convergence is demonstrated here through the examination of the Flamenco results alone, as detailed in Sec. II D, using simulations ranging from 180 × 1282 to 720 × 5122. Figure 28 shows the convergence of mean properties at t = 1 s. The two highest grid levels are in very good agreement. Figure 29 shows the grid convergence of the individual terms in each transport equation for Flamenco. Again, the two highest resolutions agree well, particularly considering the higher-order nature of the metrics under examination. Finally, Fig. 27 shows the convergence of computed numerical dissipation for Flamenco. For the mean momentum and mass fraction, numerical dissipation is small and dominated by noise. For the mass fraction variance and kinetic energy, numerical dissipation is significant and the difference between the two highest grid resolutions is again small.

FIG. 28.

Convergence of mean properties at 1 s for Flamenco, (a) ũutrans, , w ̃ , (b) m 1 ̃ , (c) E ̃ , and (d) m 1 ̃ .

FIG. 28.

Convergence of mean properties at 1 s for Flamenco, (a) ũutrans, , w ̃ , (b) m 1 ̃ , (c) E ̃ , and (d) m 1 ̃ .

Close modal
FIG. 29.

Convergence of individual terms in the (a) ũ, (b) m 1 ̃ , (c) E ̃ , and (d) m 1 ̃ transport equations at t = 1 s.

FIG. 29.

Convergence of individual terms in the (a) ũ, (b) m 1 ̃ , (c) E ̃ , and (d) m 1 ̃ transport equations at t = 1 s.

Close modal
1.
R. D.
Richtmyer
, “
Taylor instability in shock acceleration of compressible fluids
,”
Commun. Pure Appl. Math.
13
,
297
319
(
1960
).
2.
E. E.
Meshkov
, “
Instability of the interface of two gases accelerated by a shock wave
,”
Izv. Akad. Nauk. SSSR, Mekh. Zhidk. Gaza
4
,
151
157
(
1969
).
3.
D. S.
Clark
 et al, “
Three-dimensional simulations of low foot and high foot implosion experiments on the National Ignition Facility
,”
Phys. Plasmas
23
,
056302
(
2016
).
4.
A.
Burrows
, “
Supernova explosions in the Universe
,”
Nature
403
,
723
733
(
2000
).
5.
J.
Yang
,
T.
Kubota
, and
E. E.
Zukoski
, “
Applications of shock-induced mixing to supersonic combustion
,”
AIAA J.
31
,
854
862
(
1993
).
6.
Y.
Zhou
, “
Rayleigh-Taylor and Richtmyer-Meshkov instability induced flow, turbulence, and mixing. I
,”
Phys. Rep.
720-722
,
1
136
(
2017
).
7.
Y.
Zhou
, “
Rayleigh-Taylor and Richtmyer-Meshkov instability induced flow, turbulence, and mixing. II
,”
Phys. Rep.
723-725
,
1
160
(
2017
).
8.
G. I.
Barenblatt
,
G.
Looss
, and
D. D.
Joseph
,
Nonlinear Dynamics and Turbulence
(
Pitman Publishing
,
1983
).
9.
B.
Thornber
,
D.
Youngs
,
D.
Drikakis
, and
R. J. R.
Williams
, “
The influence of initial conditions on turbulent mixing due to Richtmyer-Meshkov instability
,”
J. Fluid Mech.
654
,
99
139
(
2010
).
10.
Y.
Elbaz
and
D.
Shvarts
, “
Modal model mean field self-similar solutions to the asymptotic evolution of Rayleigh-Taylor and Richtmyer-Meshkov instabilities and its dependence on the initial conditions
,”
Phys. Plasmas
25
,
062126
(
2018
).
11.
C. D.
Stefano
,
G.
Malamud
,
C.
Kuranz
,
S.
Klein
, and
R.
Drake
, “
Measurement of Richtmyer–Meshkov mode coupling under steady shock conditions and at high energy density
,”
High Energy Density Phys.
17
,
263
269
(
2015
).
12.
K. A.
Flippo
,
F. W.
Doss
,
J. L.
Kline
,
E. C.
Merritt
,
D.
Capelli
,
T.
Cardenas
,
B.
DeVolder
,
F.
Fierro
,
C. M.
Huntington
,
L.
Kot
,
E. N.
Loomis
,
S. A.
MacLaren
,
T. J.
Murphy
,
S. R.
Nagel
,
T. S.
Perry
,
R. B.
Randolph
,
G.
Rivera
, and
D. W.
Schmidt
, “
Late-time mixing sensitivity to initial broadband surface roughness in high-energy-density shear layers
,”
Phys. Rev. Lett.
117
,
225001
(
2016
).
13.
V. V.
Krivets
,
K. J.
Ferguson
, and
J. W.
Jacobs
, “
Turbulent mixing induced by Richtmyer-Meshkov instability
,”
AIP Conf. Proc.
1793
,
150003
(
2017
).
14.
O.
Soulard
,
F.
Guillois
,
J.
Griffond
,
V.
Sabelnikov
, and
S.
Simoëns
, “
Permanence of large eddies in Richtmyer-Meshkov turbulence with a small Atwood number
,”
Phys. Rev. Fluids
3
,
104603
(
2018
).
15.
O.
Schilling
and
M.
Latini
, “
High-order WENO simulations of three-dimensional reshocked Richtmyer-Meshkov instability to late times: Dynamics, dependence on initial conditions, and comparisons to experimental data
,”
Acta Math. Sci.
30
,
595
620
(
2010
).
16.
B.
Thornber
,
J.
Griffond
,
O.
Poujade
,
N.
Attal
,
H.
Varshochi
,
P.
Bigdelou
,
P.
Ramaprabhu
,
B.
Olson
,
J.
Greenough
,
Y.
Zhou
,
O.
Schilling
,
K. A.
Garside
,
R. J. R.
Williams
,
C. A.
Batha
,
P. A.
Kuchugov
,
M. E.
Ladonkina
,
V. F.
Tishkin
,
N. V.
Zmitrenko
,
V. B.
Rozanov
, and
D. L.
Youngs
, “
Late-time growth rate, mixing, and anisotropy in the multimode narrowband Richtmyer-Meshkov instability: The θ-group collaboration
,”
Phys. Fluids
29
,
105107
(
2017
).
17.
A.
Favre
, “
Équations statistiques aux fluctuations d’entropie, de concentration, de rotationnel dans les écoulements compressible
,”
C.R. Acd. Sci. Paris
273
,
1289
1294
(
1958
).
18.
P.
Chassaing
, “
The modeling of variable density turbulent flows. A review of first-order closure schemes
,”
Flow, Turbul. Combust.
66
,
293
332
(
2001
).
19.
G.
Dimonte
and
R.
Tipton
, “
K-l turbulence model for the self-similar growth of the Rayleigh-Taylor and Richtmyer-Meshkov instabilities
,”
Phys. Fluids
18
,
085101
(
2006
).
20.
S.
Gauthier
and
M.
Bonnet
, “
A k–ε model for turbulent mixing in shock-tube flows induced by Rayleigh-Taylor instability
,”
Phys. Fluids A
2
,
1685
1694
(
1990
).
21.
J. T.
Morán-López
and
O.
Schilling
, “
Multicomponent Reynolds–averaged Navier–Stokes simulations of reshocked Richtmyer–Meshkov instability-induced mixing
,”
High Energy Density Phys.
9
,
112
121
(
2013
).
22.
D.
Besnard
,
F.
Harlow
,
R.
Rauenzahn
, and
C.
Zemach
, “
Turbulence transport equations for variable-density turbulence and their relationship to two-field models
,” Technical Report No. LA-12303-MS,
Los Alamos
,
1992
.
23.
A.
Banerjee
,
R.
Gore
, and
M.
Andrews
, “
Development and validation of a turbulent-mix model for variable-density and compressible flows
,”
Phys. Rev. E
82
,
046309
(
2010
).
24.
K.
Stalsberg-Zarling
and
R.
Gore
, “
The BHR2 turbulence model incompressible isotropic decay, Rayleigh-Taylor, Kelvin-Helmholtz and homogeneous variable density turbulence
,” Technical Report No. LA-UR 11-04773,
Los Alamos
,
2011
.
25.
J. D.
Schwarzkopf
,
D.
Livescu
,
R. A.
Gore
,
R. M.
Rauenzahn
, and
J. R.
Ristorcelli
, “
Application of a second-moment closure model to mixing processes involving multicomponent miscible fluids
,”
J. Turbul.
12
,
N49
(
2011
).
26.
J. D.
Schwarzkopf
,
D.
Livescu
,
J. R.
Baltzer
,
R. A.
Gore
, and
J.
Ristorcelli
, “
A two-length scale turbulence model for single-phase multi-fluid mixing
,”
Flow, Turbul. Combust.
96
,
1
43
(
2016
).
27.
O.
Grégoire
,
D.
Souffland
, and
S.
Gauthier
, “
A second-order turbulence model for gaseous mixtures induced by Richtmyer-Meshkov instability
,”
J. Turbul.
6
,
1
20
(
2005
).
28.
D.
Souffland
,
O.
Soulard
, and
J.
Griffond
, “
Modeling of Reynolds stress models for diffusion fluxes inside shock waves
,”
ASME J. Fluids Eng.
136
,
091102
(
2014
).
29.
D. L.
Youngs
, “
Numerical simulation of mixing by Rayleigh-Taylor and Richtmyer-Meshkov instabilities
,”
Laser Part. Beams
12
,
725
750
(
1994
).
30.
A.
Llor
,
Statistical Hydrodynamic Models for Developing Mixing Instability Flows
(
Springer
,
Berlin
,
2005
).
31.
A.
Llor
and
P.
Bailly
, “
A new turbulent two-field concept for modeling Rayleigh–Taylor, Richtmyer–Meshkov, and Kelvin–Helmholtz mixing layers
,”
Laser Part. Beams
21
,
311
315
(
2003
).
32.
E.
Leinov
,
G.
Malamud
,
Y.
Elbaz
,
L.
Levin
,
G.
Ben-Dor
,
D.
Shvarts
, and
O.
Sadot
, “
Experimental and numerical investigation of the Richtmyer-Meshkov instability under re-shock conditions
,”
J. Fluid Mech.
626
,
449
475
(
2009
).
33.
G.
Malamud
,
E.
Leinov
,
O.
Sadot
,
Y.
Elbaz
,
G.
Ben-Dor
, and
D.
Shvarts
, “
Reshocked Richtmyer-Meshkov instability: Numerical study and modeling of random multi-mode experiments
,”
Phys. Fluids
26
,
084107
(
2014
).
34.
M.
Mohaghar
,
J.
Carter
,
G.
Pathikonda
, and
D.
Ranjan
, “
The transition to turbulence in shock-driven mixing: Effects of Mach number and initial conditions
,”
J. Fluid Mech.
871
,
595
635
(
2019
).
35.
C.
Weber
,
N.
Haehn
,
J.
Oakley
,
D.
Rothamer
, and
R.
Bonazza
, “
An experimental investigation of the turbulent mixing transition in the Richtmyer–Meshkov instability
,”
J. Fluid Mech.
748
,
457
487
(
2014
).
36.
D. T.
Reese
,
A. M.
Ames
,
C. D.
Noble
,
J. G.
Oakley
,
D. A.
Rothamer
, and
R.
Bonazza
, “
Simultaneous direct measurements of concentration and velocity in the Richtmyer–Meshkov instability
,”
J. Fluid Mech.
849
,
541
575
(
2018
).
37.
D. S.
Clark
,
A. L.
Kritcher
,
S. A.
Yi
,
A. B.
Zylstra
,
S. W.
Haan
, and
C. R.
Weber
, “
Capsule physics comparison of National Ignition Facility implosion designs using plastic, high density carbon, and beryllium ablators
,”
Phys. Plasmas
25
,
032703
(
2018
).
38.
S. R.
Nagel
,
K. S.
Raman
,
C. M.
Huntington
,
S. A.
MacLaren
,
P.
Wang
,
M. A.
Barrios
,
T.
Baumann
,
J. D.
Bender
,
L. R.
Benedetti
,
D. M.
Doane
,
S.
Felker
,
P.
Fitzsimmons
,
K. A.
Flippo
,
J. P.
Holder
,
D. N.
Kaczala
,
T. S.
Perry
,
R. M.
Seugling
,
L.
Savage
, and
Y.
Zhou
, “
A platform for studying the Rayleigh–Taylor and Richtmyer–Meshkov instabilities in a planar geometry at high energy density at the National Ignition Facility
,”
Phys. Plasmas
24
,
072704
(
2017
).
39.
W. C.
Wan
,
G.
Malamud
,
A.
Shimony
,
C. A.
Di Stefano
,
M. R.
Trantham
,
S. R.
Klein
,
D.
Shvarts
,
C. C.
Kuranz
, and
R. P.
Drake
, “
Observation of single-mode, Kelvin-Helmholtz instability in a supersonic flow
,”
Phys. Rev. Lett.
115
,
145001
(
2015
).
40.
M. D.
Slessor
,
C. L.
Bond
, and
P. E.
Dimotakis
, “
Turbulent shear-layer mixing at high Reynolds numbers: Effects of inflow conditions
,”
J. Fluid Mech.
376
,
115
138
(
1998
).
41.
O.
Schilling
and
N. J.
Mueschke
, “
Analysis of turbulent transport and mixing in transitional Rayleigh-Taylor unstable flow using direct numerical simulation data
,”
Phys. Fluids
22
,
105102
(
2010
).
42.
R. E.
Duff
,
F. H.
Harlow
, and
C. W.
Hirt
, “
Effects of diffusion on interface instability between gases
,”
Phys. Fluids
5
,
417
425
(
1962
).
43.
G.
Dimonte
and
M.
Schneider
, “
Density ratio dependence of Rayleigh-Taylor mixing for sustained and impulsive acceleration histories
,”
Phys. Fluids
12
,
304
321
(
2000
).
44.
A.
Shimony
,
G.
Malamud
, and
D.
Shvarts
, “
Density ratio and entrainment effects on asymptotic Rayleigh–Taylor instability
,”
J. Fluids Eng.
140
,
050906
(
2018
).
45.
K. O.
Mikaelian
, “
Rayleigh-Taylor and Richtmyer-Meshkov instabilities in multilayer fluids with surface tension
,”
Phys. Rev. A
42
,
7211
7225
(
1990
).
46.
B.
Thornber
,
A.
Mosedale
,
D.
Drikakis
,
D.
Youngs
, and
R.
Williams
, “
An improved reconstruction method for compressible flows with low Mach number features
,”
J. Comput. Phys.
227
,
4873
4894
(
2008
).
47.
A.
Garcia-Uceda Juarez
,
A.
Raimo
,
E.
Shapiro
, and
B.
Thornber
, “
Steady turbulent flow computations using a low Mach fully compressible scheme
,”
AIAA J.
52
,
2559
2575
(
2014
).
48.
B.
Fryxell
,
K.
Olson
,
P.
Ricker
,
F. X.
Timmes
,
M.
Zingale
,
D. Q.
Lamb
,
P.
MacNeice
,
R.
Rosner
,
J. W.
Truran
, and
H.
Tufo
, “
Flash: An adaptive mesh hydrodynamics code for modeling astrophysical thermonuclear flashes
,”
Astrophys. J. Suppl.
131
,
273
(
2000
).
49.
P.
Colella
and
P. R.
Woodward
, “
The piecewise parabolic method (PPM) for gas-dynamical simulations
,”
J. Comput. Phys.
54
,
174
201
(
1984
).
50.
B.
van Leer
, “
Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow
,”
J. Comput. Phys.
23
,
263
275
(
1977
).
51.
R.
LeVeque
,
Finite Volume Methods for Hyperbolic Problems
, Cambridge Texts in Applied Mathematics Vol. 31 (
Cambridge University Press
,
2002
).
52.
V.
Daru
and
C.
Tenaud
, “
High order one-step monotonicity-preserving schemes for unsteady compressible flow calculations
,”
J. Comput. Phys.
193
,
563
594
(
2004
).
53.
Implicit Large Eddy Simulation: Computing Turbulent Fluid Dynamics
, edited by
F. F.
Grinstein
,
L. G.
Margolin
, and
W. J.
Rider
(
Cambridge University Press
,
Cambridge
,
2007
).
54.
D. L.
Youngs
, “
Three-dimensional numerical simulation of turbulent mixing by Rayleigh-Taylor instability
,”
Phys. Fluids A
3
,
1312
1320
(
1991
).
55.
R. J. R.
Williams
, “
Sub-grid properties and artificial viscous stresses in staggered-mesh schemes
,”
J. Comput. Phys.
374
,
413
443
(
2018
).
56.
B.
Thornber
and
Y.
Zhou
, “
Energy transfer in the Richtmyer-Meshkov instability
,”
Phys. Rev. E
86
,
056302
(
2012
).
57.
B.
Thornber
, “
Impact of domain size and statistical errors in simulations of homogeneous decaying turbulence and the Richtmyer-Meshkov instability
,”
Phys. Fluids
28
,
045106
(
2016
).
58.
B.
Thornber
,
D.
Drikakis
,
D.
Youngs
, and
R. J. R.
Williams
, “
Physics of the single-shocked and reshocked Richtmyer-Meshkov instability
,”
J. Turbul.
13
,
1
17
(
2012
).
59.
A. W.
Cook
and
Y.
Zhou
, “
Energy transfer in Rayleigh-Taylor instability
,”
Phys. Rev. E
66
,
026312
(
2002
).
60.
S.
Sarkar
, “
The pressure-dilatation correlation in compressible flows
,”
Phys. Fluids A
4
,
2674
2682
(
1992
).
61.
J.
Ristorcelli
, “
A pseudo-sound constitutive relationship for the dilatational covariances in compressible turbulence
,”
J. Fluid Mech.
347
,
37
70
(
1997
).