Numerical Simulations of Shock-driven Heavy Fluid Layer
Satyvir Singh
Institute for Applied and Computational Mathematics, RWTH Aachen University, Germany
E-mail: singh@acom.rwth-aachen.de
Received 02 November 2024; Accepted 05 January 2025
The present study presents the numerical simulations for a shocked-heavy fluid layer with a stratified N/SF/N configuration. Simulations were conducted using a third-order modal discontinuous Galerkin method to solve the compressible two-component Euler equations. The results were validated against experimental data, confirming the accuracy of the computational approach. Dynamics of the heavy fluid layer were found to be strongly influenced by the shock Mach numbers M 1.15, 1.25, 1.5. At a lower Mach number M 1.15, the interface deformations remained smooth and relatively symmetric, with limited vorticity generation and weak perturbations. Baroclinic effects at this stage were minimal, and the instability growth remained linear. As the Mach number increased to M 1.25, the interaction became nonlinear, leading to the formation of small-scaled vortex structures driven by moderate baroclinic effects. Interface mixing intensified as rotational motion increased. At the highest Mach number M 1.5, the interface rapidly evolved into chaotic structures, characterized by significant vorticity amplification, vortex roll-up, and the onset of turbulence. The baroclinic vorticity, resulting from the misalignment of pressure and density gradients, dominated the vorticity production mechanism, particularly at higher Mach numbers. Quantitative analysis demonstrated that average vorticity, baroclinic vorticity, and enstrophy grew rapidly with increasing Mach numbers. Enstrophy, which quantifies turbulence intensity, exhibited pronounced growth at M 1.5, marking the transition to turbulent mixing.
Keywords: Shock wave, heavy fluid layer, vorticity, interface deformation.
The Richtmyer-Meshkov (RM) instability occurs when a shock wave or an impulsive acceleration disturbs the interface separating two fluids [1, 2]. This instability is triggered by the generation of baroclinic vorticity, which arises when the pressure gradient from the shock wave is misaligned with the density gradient at the interface. In its early stages, the perturbation amplitude grows linearly, but it eventually transitions into a nonlinear regime, forming prominent structures like spikes and bubbles [3]. Over time, this process may evolve further, leading to turbulent mixing. The RM instability plays a vital role in various engineering and scientific contexts. For instance, in inertial confinement fusion (ICF), the RM instability proves detrimental as it contaminates the fusion fuel within the capsule, reducing the thermonuclear yield [3]. Beyond engineering, the RM instability provides insights into phenomena related to supernovae and stellar evolution [4]. Both ICF capsules and supernova consist of multiple material layers with distinct initial perturbations at their interfaces. Essentially, the RM instability occurs at fluid layers with finite thickness and discrete boundaries, where the instability growth at one interface influences that of adjacent ones. This complex interaction, termed interface coupling, arises due to varying initial conditions at each interface [5, 6]. Additionally, the instability growth rates at the interfaces often differ, resulting in intricate mixing across the fluid layer. Therefore, understanding the instability behavior in fluid layers with finite thickness and two distinct interfaces is crucial for unravelling such phenomena.
Numerous studies have been conducted on the simplified fluid-layer scenario, in which the single-mode perturbations at both interfaces are identical [7]. Taylor initially examined this scenario theoretically in the context of Rayleigh-Taylor (RT) instability [8, 9]. After then, Mikaelian extended the investigation to fluid layers in RM instability and created an analytical model to account for the two interfaces’ linear amplitude rise [10, 11]. To shed light on the physics of RM instability mixing, Balakumar et al. [12] experimentally examined the characteristics of turbulent mixing in an RM instability caused by a fluid layer under the influence of a single shock and a reshock wave. Jacobs et al. [13] developed an alternative linear model by solving velocity potential functions. While the formulation differs from the approach of Mikaelian [11], it yields predictions that are nearly identical to his approach. Liang and Luo [14] expanded the model introduced by Jacobs et al. [13] to include a layer that divides three fluids with varying densities. Later, Liang et al. [15] investigated the impact of interface coupling and reverberating waves on perturbation growth by constructing several fluid layers with different gas mixtures and thicknesses using the soap-film approach. The fundamental mechanism of the RM instability under the influence of a reshock wave was recently revealed by Li et al. [16], which numerically examined the evolution of a shock-induced fluid layer. Recently, Guo et al. [17] studied the RM instability finger collisions in light fluid layers under reshock conditions through shock-tube experiments.
In this work, we present the numerical simulations for a shock-driven heavy fluid layer with two different interfaces. A fluid layer with composed stratified N/SF/N configuration is considered as a heavy fluid layer, which involves two different kinds of interfaces, such as V-shaped (I1) and fixed (I2). For numerical simulations, we apply a high-order modal discontinuous Galerkin solver to investigate the RM instability at the heavy fluid layer. The goal is to analyse the effect of shock Mach number on the evolution of flow structures, vorticity generation, and mixing dynamics as the shock interacts with this fluid layer. By leveraging a high-order DG scheme, the present study captures the detailed wave interactions, vortex dynamics, and instability growth that characterize the RM instability. The remainder of this investigation is organized as follows: Section 2 presents the initial setup and the governing equations used to simulate the RM instability. Section 3 details the numerical solver employed and the validation process. Section 4 discusses the numerical results and provides an in-depth analysis of the shocked heavy fluid layer. Finally, conclusions are drawn in Section 5.
Figure 1 shows the initial configurations of the shock wave interaction with the N/SF/N fluid layer used in this study. The numerical simulations are conducted within a rectangular domain measuring . I1 is the first V-shaped interface in this fluid layer, and I2 is the second fixed interface. While, mm denotes the fluid-layer thickness. The I1 interface is located downstream of the shock, with the left end positioned 5 mm away from the incident shock (IS) wave moving from left to right along the x-direction, and situated in N gas. Moreover, 35 mm is the IS wave’s beginning distance from the domain’s left boundary. The Mach number M characterizes the IS wave. In the present setup, the total length of the I1 interface, which is twice the standard amplitude of a single-mode interface, is referred to as the initial amplitude (). The wavelength of the I1 interface is determined by taking . The vertex angle of the I1 interface is . The relationship between , , and for this V-shaped geometry is /. The gives the location of the I1 interface. We take into account a vertex angle of for numerical simulations. For the shock wave strengths, three Mach numbers, M 1.15, 1.25, and 1.5, are taken into consideration. In this computational domain, the top, bottom, and right boundaries function as outlets, while the left boundary acts as the inflow boundary. Around the fluid layer, the initial pressure and temperature are specified as Pa and K, respectively.
Figure 1 Schematic description of the initial configurations for a shocked heavy fluid layer. Here d denotes the fluid-layer thickness, I1 denotes the initial first interface, I2 denotes the initial second interface. , and are the vortex angle, wavelength, amplitude of the I1 interface, respectively. The solid arrow illustrates the direction of shock wave propagation.
The compressible two-component Euler equations are numerically solved in the current Work. These equations can be re-written in using the conservative form as
(1) |
where is the density; and are the velocity components in the and directions, respectively; is the energy, is the mass fraction; and is the pressure which is evaluated from the mathematical expression as follows:
(2) |
In this context, represents the mixture’s specific heat ratio. The equation that describes the relationship between pressure (), density (), temperature (), and the mixture-specific gas constant (R) is given is . The assumption is made that both gas components are in thermal equilibrium and behave as calorically perfect gases. These gases are characterized by their specific heats at constant pressure (, ), specific heats at constant volume (, ), and specific heat ratios (, ). The specific heat ratio of a mixture can be evaluated by
(3) |
where subscripts 1 and 2 denotes for bubble and ambient gas, respectively. Remarkably, in numerical experiments, spurious oscillations can be caused by the jump in the specific heat ratio across an interface, especially in issues involving compressible flow. Adaptive strategies, artificial viscosity, boundary conditions, and careful numerical method selection are frequently needed to address this problem and reduce or eliminate oscillations while preserving simulation accuracy.
In this study, the two-dimensional system of a two-component compressible Euler equations is solved using in-house developed explicit modal discontinuous Galerkin solver [18]. Scaled Legendre polynomial functions that are divided into non-overlapping rectangular components are used in the computational domain. Using the Gauss–Legendre quadrature rule, volume and flux are both integrated [18]. Numerical fluxes at the elemental interfaces for two-component flows are calculated using the HLLC scheme. An accurate scaled Legendre polynomial expansion of third order is used to approximate the solutions in the finite element space. An explicit third-order accurate strong stability-preserving (SSP) Runge–Kutta approach is used to integrate time. Furthermore, the computational solutions reduce non-physical oscillations by using a high-order moment limiter [19].
To validate the current numerical solver, the numerical results are compared with those obtained from experiments [20]. In this validation, a V-shaped air/SF interface was taken into consideration. During the experimental examination, a shock Mach number of M was utilized, and a vertex angle of occurred. As shown in Figure 2(a), Schlieren image comparison between our current findings and the experimental results is presented. As the shock wave advances, our observations show that the V-shaped contact becomes increasingly distorted. This interface distortion is consistent with findings from a related experiment, as are the ensuing intricate wave patterns. Furthermore, Figure 2(b) shows the time evolution of the upstream interface displacement, represented by the symbol Ds. It can be seen that the numerical simulation precisely reproduces the positions of Ds and closely aligns with the experimental data.
Figure 2 Comparison of experiments [20] and numerical results for a shock wave interaction with V-shaped air/SF interface: (a) numerical Schlieren, (b) variation in the upstream interface displacement.
In this section, we present the numerical results for the shock-driven N/SF/N layer Emphasis is placed on the evolution of flow morphology, wave patterns, vorticity generation, and enstrophy. In the following numerical simulations, we use the normalized time to display flow morphology snapshots. The characteristic time is used to standardize the actual computational time, where is the initial wavelength of the V-shaped interface and is the IS wave speed. Thus we obtain dimensionless time scale as . We employ grid points for the numerical simulations in order to accurately represent the intricate flow field structure of the RM instability at heavy fluid layer.
Figure 3 Flow field evolution in the shock-driven N/SF/N layer: density contours with three different shock Mach numbers (M, 1.25, 1.5).
Figure 3 illustrates the flow field evolution of the N/SF/N layer at three different shock Mach numbers (M, 1.25, 1.5) through density contour plots at varying times (). This figure effectively demonstrates the influence of shock strength on the RMI showing the progression from smooth interface deformation to turbulent mixing as Mach number increases. At the initial shock Interaction (), the shock wave interacts with the interface of the heavier SF gas layer (initially stratified between N layers). A sharp interface deformation occurs, with a slight tilt toward the shock propagation direction. The interface shows a primary shock compression effect, visible as a small wedge-like deformation. As time goes (), the fluid layer deformation remains smooth and relatively symmetric at M. The amplitude of the disturbance increases, but the flow is less chaotic. More significant deformation occurs due to stronger compressibility effects. At M, vortex structure starts to emerge, especially at the corners of the interface. While, at high Mach number i.e. M, The fluid layer shows strong instability and vortex formation. The deformation is much larger, and small-scale structures (vorticity) appear around the interface. During intermediate development (), the fluid layer becomes more nonlinear, with pronounced finger-like protrusions. At high Mach number (M), the development of instabilities is much faster. The interface breaks down into fine turbulent structures. The flow structures at M and M show clear vortex roll-ups, indicating stronger growth of the RMI. Subsequently, later time evolution (), the fluid layer at M remains less turbulent, with elongated smooth structures. Interestingly, the fluid layer at M breaks into multiple vortices, and mixing intensifies. Turbulent structures are visible but less chaotic compared to M. The fluid layer at M becomes fully turbulent with significant mixing and small-scale structures dominating the flow.
Figure 4 Vorticity generation in the shock-driven N/SF/N layer: vorticity contours with three different shock Mach numbers (M, 1.25, 1.5).
Figure 4 illustrates the vorticity generation process in the shocked N/SF/N gas layer for three different shock Mach numbers (M, 1.25, 1.5) over time. Vorticity fields are visualized at different dimensionless times (). At all Mach numbers, the interface begins to distort due to the shock wave interacting with the density gradient (). Vorticity appears as thin red (positive) and blue (negative) bands, representing shear layers. It is found that the vorticity is weak and localized along the interface at M. While, the vorticity is stronger, with more prominent positive and negative regions forming at M. During the intermediate stages (), vorticity remains weak and the shear layers are stable with minimal deformation at M. The interface is stretched but not significantly perturbed. As Mach number increases, vorticity strengthens and small-scale structures begin to appear M. The interface shows slight perturbations, indicating the early stages of instability. At high Mach number i.e. M, vorticity generation becomes highly nonlinear. Small vortices start forming due to strong Baroclinic effects, and the interface is highly distorted. At later times i.e. , the interface remains mostly smooth, with vorticity confined to the shear layers for M. No significant secondary structures are observed. For M, vorticity intensifies, and small vortices appear along the interface. The interface shows noticeable deformation and mixing. Interestingly, the vorticity field becomes highly turbulent, with multiple vortices rolling up and interacting at M. The interface is completely distorted, showing significant mixing regions and complex vortex interactions.
Now, we investigate the vorticity generation in the shocked heavy fluid layer by analysing two spatially integrated fields: average and baroclinic vorticities. The spatially integrated field of average vorticity is defined as
(4) |
where denotes the vorticity. The spatial integrated field of baroclinic vorticity production is given by
(5) |
Figure 5 presents spatially integrated fields of average vorticity and baroclinic in the shocked N/SF/N layer for different shock Mach numbers (M, 1.25, 1.5). It can be seen from Figure 5(a) that higher Mach numbers (M) produce significantly higher average vorticity over time compared to lower Mach numbers (M) and (M). For M, the average vorticity grows rapidly and shows a nearly linear increase at later times, indicating sustained rotational motion. At M, the average vorticity increases more gradually. While, for M, the growth remains small, reflecting weaker shock-induced rotational motion. On other hand, baroclinic vorticity arises from the Baroclinic effect, which occurs when pressure and density gradients are misaligned. This effect is the primary source of vorticity generation during shock interaction with an interface, as illustrated in Figure 5(b). Similar to average vorticity, baroclinic vorticity grows much faster for higher Mach numbers (M) compared to lower Mach numbers. For M, baroclinic vorticity increases steadily and significantly over time, reaching the highest values. While, for M, the growth rate is slower but still notable, while for M, it remains quite small. It can be found that increasing the Mach number intensifies both average and baroclinic vorticity. This happens because stronger shocks generate larger pressure gradients, enhancing the baroclinic effect and consequently the vorticity. Both quantities show an initial rise immediately after shock interaction, followed by sustained growth, especially at higher Mach numbers (M). The growth rate is higher for baroclinic vorticity compared to average vorticity, highlighting the dominance of the baroclinic mechanism in driving vorticity production.
Figure 5 Spatially integrated fields of (a) average vorticity, and (b) baroclinic vorticity in the shocked N/SF/N layer with different shock Mach numbers (M, 1.25, 1.5).
Lastly, we examine the spatially integrated field of enstrophy, which aids in comprehending the physical mechanisms underlying the interaction’s vorticity creation or attenuation. It is characterized by
(6) |
Figure 6 illustrates the spatially integrated of enstrophy for the shock-driven N/SF/N layer with different shock Mach numbers (M, 1.25, 1.5). When a shock wave interacts with a heterogeneous fluid layer, it causes baroclinic vorticity generation due to the misalignment between pressure and density gradients. This vorticity evolves into turbulent mixing regions, increasing the enstrophy of the system. Higher shock Mach numbers produce stronger shocks, leading to larger pressure gradients and greater vorticity generation. Enstrophy increases rapidly at the start due to the immediate shock impact on the interface, generating vorticity. It can be seen that for M, the enstrophy grows significantly faster and reaches much higher values compared to lower Mach numbers. Also, the enstrophy shows oscillations and a significant long-term rise, likely due to secondary instabilities and turbulence. For lower Mach numbers (M), enstrophy saturates at a relatively low value because the shock strength is insufficient to trigger significant turbulence.
Figure 6 Spatially integrated fields of enstrophy in the shocked N/SF/N layer with different shock Mach numbers (M, 1.25, 1.5).
In this study, the evolution of the RM instability in a heavy fluid layer has been thoroughly examined by analysing the interaction of a shock wave with a stratified N/SF/N configuration. Key observations reveal that the instability dynamics are significantly influenced by the strength of the shock wave, characterized by different Mach numbers (M, 1.25, 1.5). A third-order modal discontinuous Galerkin approach was used to simulate unstable compressible two-component Euler equations, which produced high-resolution numerical simulations. The computational model was validated with existing experimental results. The numerical results show that interface deformations remain smooth and relatively symmetric at lower Mach numbers (M), with limited vorticity generation and weak perturbations. However, as the Mach number increases, the interaction becomes progressively nonlinear, leading to vortex roll-up and the onset of turbulence. At M, moderate baroclinic effects drive the formation of vortical structures and enhance interface mixing. For high Mach numbers (M), the interface evolves rapidly into highly chaotic structures, exhibiting significant vorticity amplification, turbulent mixing, and enstrophy growth. Quantitative analysis of average vorticity, baroclinic vorticity, and enstrophy confirms that stronger shocks intensify rotational motion and accelerate the instability growth. The baroclinic vorticity, arising from the misalignment of pressure and density gradients, dominates the vorticity production mechanism, especially at M. The enstrophy, which quantifies turbulence intensity, demonstrates rapid growth for higher Mach numbers, reflecting the transition to turbulent mixing.
The author acknowledges funding through the German Research Foundation within the Research Unit DFG–FOR5409.
[1] R.D. Richtmyer, ‘Taylor instability in shock acceleration of compressible fluids,’ Communications on Pure and Applied Mathematics, 13:297, 1960.
[2] E.E. Meshkov, ‘Instability of the interface of two gases accelerated by a shock wave,’ Fluid Dynamics, 4:101–104, 1969.
[3] S. Singh, ‘Role of Atwood number on flow morphology of a planar shock-accelerated square bubble: A numerical study,’ Physics of Fluids, 32: 126112, 2020.
[4] J. Lindl, O. Landen, J. Edwards, E. Moses, and NIC team, ‘Review of the national ignition campaign 2009-201,’ Physics of Plasmas, 21: 020501, 2014.
[5] D. Arnett, ‘The role of mixing in astrophysics,’ The Astrophysical Journal Supplement Series, 127: 2131, 2000.
[6] K.O. Mikaelian, ‘The role of mixing in astrophysics,’ Physical Review A, 28: 1637–1646, 1983.
[7] J.W. Jacobs, D.G. Jenkins, D.L. Klein, and R.F. Benjamin, ‘Nonlinear growth of the shock-accelerated instability of a thin fluid layer,’ Physical Review A, 28: 1637–1646, 1983.
[8] G.I. Taylor, ‘The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. I,’ Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 201: 192–196, 1950.
[9] L. Rayleigh, ‘Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density,’ Proceedings of the London Mathematical Society, 14: 170–177, 1883.
[10] K.O. Mikaelian, ‘Richtmyer-Meshkov instabilities in stratified fluids,’ Physical Review A, 31: 410–419, 1985.
[11] K.O. Mikaelian, ‘Rayleigh–Taylor and Richtmyer–Meshkov instabilities in finite thickness fluid layers,’ Physics of Fluids, 7: 888–890, 1995.
[12] B.J. Balakumar, G.C. Orlicz, J.R. Ristorcelli, S. Balasubramanian, K.P. Prestridge, and C.D. Tomkins, ‘Turbulent mixing in a Richtmyer–Meshkov fluid layer after reshock: velocity and density statistics,’ Journal of Fluid Mechanics, 696: 67–93, 2012.
[13] J.W. Jacobs, D.G. Jenkins, D.L. Klein, and R.F. Benjamin, ‘Turbulent mixing in a Richtmyer–Meshkov fluid layer after reshock: velocity and density statistics,’ Journal of Fluid Mechanics, 28: 23–42, 1995.
[14] Y. Liang, and X. Luo, ‘Review on hydrodynamic instabilities of a shocked gas layer,’ Science China Physics, Mechanics & Astronomy, 66: 104701, 2023.
[15] Y. Liang, and X. Luo, ‘On shock-induced heavy-fluid-layer evolution,’ Journal of Fluid Mechanics, 920: A131, 2021.
[16] L. Li, T. Jin, L. Zou, K. Luo, and J. Fan, ‘Numerical study of Richtmyer-Meshkov instability in finite thickness fluid layers with reshock,’ Physical Review E, 109: 055105, 2024.
[17] X. Guo, Z. Cong, T. Si, and X. Luo, ‘On Richtmyer–Meshkov finger collisions in a light fluid layer under reshock conditions,’ Journal of Fluid Mechanics, 1000: A87, 2024.
[18] S. Singh, A. Karchani, T. Chourushi, and R.S. Myong, ‘A three-dimensional modal discontinuous Galerkin method for second-order Boltzmann-Curtiss constitutive models of rarefied and microscale gas flow,’ Journal of Computational Physics, 457: 111052, 2022.
[19] L. Krivodonova, ‘Limiters for high-order discontinuous Galerkin methods,’ Journal of Computational Physics, 226: 879–896, 2007.
[20] X. Luo, P. Dong, T. Si, and Z. Zhai, ‘The Richtmyer–Meshkov instability of a ‘V’ shaped air/SF interface,’ Journal of Fluid Mechanics, 802: 186–202, 2016.
Satyvir Singh is currently working as a Research Associate Fellow in the Institute of Applied and Computational Mathematics at RWTH Aachen University, Germany. Dr. Singh earned his Ph.D. in Computational Fluid mechanics in the School of Mechanical and Aerospace Engineering at Gyeongsang National University, South Korea. Subsequently, he worked as a Senior Research Fellow at Research Center for Aircraft Parts Technology, Gyeongsang National University, South Korea in 2018. After then, Dr. Singh worked as a Research Fellow in School of Physical and Mathematical Sciences at Nanyang Technological University Singapore. Dr. Singh completed Master Degree M.Tech. in Industrial Mathematics & Scientific Computing at Indian Institute of Technology Madras, India. He qualified two highly competitive Indian examinations – Junior Research Fellowship and National Eligibility Test in Mathematical Sciences (2011) with All Indian Rank – 38, and Graduate Aptitude Test for Engineering in Mathematics (2012) with All Indian Rank – 244. As a teaching background, Dr. Singh worked in India as Assistant Professor in Galgotias College of Engineering & Technology Greater Noida, and IMS Engineering College Ghaziabad. Dr. Singh has a vast research area, including computational fluid dynamics, high-order numerical methods, hydrodynamic instability, gas kinetic theory, heat and mass transfer, and computational biology. His commitment to research is reflected in his more than 50 research articles (700+ citations) in reputable journals, including Physics of Fluids, J. Computational Physics, Computers & Fluids, Int. J. Heat and Mass Transfer, Physical Review Fluids, Applied Mathematics Computations, etc. Also, he has published one book. Besides it, he has attended many international conferences and presented his research work in USA, UK, Italy, South Korea, Singapore, Germany, Greece, China, Japan and India. As a Co-PI, Dr. Singh has received the research fund for the project “Mathematical Modelling and High-fidelity Simulations for Brain Tumor Dynamics” by the Deanship of Graduate Studies and Scientific Research Scheme, Jazan University, Saudi Arabia (2024).
Journal of Graphic Era University, Vol. 13_1, 75–90.
doi: 10.13052/jgeu0975-1416.1314
© 2025 River Publishers