Richtmyer-Meshkov Instability at high Mach Number: Non-Newtonian Effects

The Richtmyer-Meshkov instability (RMI) occurs when a shock wave passes through an interface between fluids of different densities, a phenomenon prevalent in a variety of scenarios including supersonic combustion, supernovae, and inertial confinement fusion. In the most advanced current numerical modelling of RMI, a multitude of secondary physical phenomena are typically neglected that may crucially change in silico predictions. In this study, we investigate the effects of shear-thinning behaviour of a fluid on the RMI at negative Atwood numbers via numerical simulations. A parametric study is carried out over a wide range of Atwood and Mach numbers that probes the flow dynamics following the impact on the interface of the initial shock wave and subsequent, reflected shocks. We demonstrate agreement between our numerical results and analytical predictions, which are valid during the early stages of the flow, and examine the effect of the system parameters on the vorticity distribution near the interface. We also carry out an analysis of the rate of vorticity production and dissipation budget which pinpoints the physical mechanisms leading to instability due to the initial and reflected shocks. Our findings indicate that the shear-thinning effects have a significant impact on instability growth and the development of secondary instabilities, which manifest themselves through the formation of Kelvin-Helmholtz waves. Specifically, we demonstrate that these effects influence vorticity generation and damping, which, in turn, affect the RMI growth. These insights have important implications for a range of applications, including inertial confinement fusion and bubble collapse within non-Newtonian materials.


I. INTRODUCTION
The Richtmyer-Meshkov instability (RMI) results from the misalignment of pressure and density gradients at the interface of two materials with different densities when a shock wave passes through.This instability generates vorticity, i.e., a rotational eddy in the local fluid flow, which leads to the growth of interface perturbations.Immediately following the passage of the shock, in the linear growth phase, the interface perturbations grow quickly, driven by their initial amplitude, wave number and Atwood number, whilst roughly maintaining their shape.Later on, however, the perturbations grow nonlinearly and couple to additional hydrodynamic instabilities, most notably the Kelvin-Helmholtz instability (KHI), resulting in the onset and development of turbulence.The Atwood number At, a dimensionless parameter that represents the ratio of the density difference to the sum of the densities of two fluids, significantly influences the progression of the RMI.This parameter quantifies the density contrast between the interacting fluids, and can be experimentally adjusted to access different growth regimes.
The RMI was first proposed by Richtmyer 1 , and the theory was later confirmed by the shock-tube experiments of Meshkov 2 and the numerical simulation of Meyer and Blewett 3 .Vandenboomgaerde et al. 4 derived a general formula for the growth rate of the RMI within the framework of the impulsive model.This formula enables the prediction of the growth rate in both 'heavy-to-light' and 'light-to-heavy' configurations.The RMI has important implications for a wide range of natural and engineering problems, including supersonic combustion, supernovae, and inertial confinement fusion 5 .For example, the occurrence of RMIs in supersonic combustion triggers better mixing and thus enhanced combustion efficiency 5 .In supernova explosions 5 , RMIs break the spherical symmetry of the progenitor star, and create conditions for the synthesis of heavy elements.In inertial confinement fusion, RMIs initiate mixing between the hotter and colder layers, resulting in a reduction of temperature that inhibits the thermonuclear burn 5 .
The impulsive model proposed by Richtmyer 1 only analyzed the light-to-heavy flow configuration, where the shock travels in the direction of the denser fluid.However, Meshkov 2 ran shock-tube experiments with both flow configurations, suggesting that the impulsive model could be applied for both light-to-heavy and heavy-to-light flow configurations.The Richtmyer impulsive model, grounded in a linear framework, posits linear growth for the perturbation amplitude, reflecting its reliance on linear theory assumptions.This model is only valid for early times when the perturbation amplitude is considerably smaller than its wavelength.As the amplitude exceeds the wavelength, nonlinear effects become significant rendering the impulsive model invalid.In response to these limitations, Yang et al. 6 introduced a refined version of Richtmyer's model to incorporate compressibility effects, which provides a more comprehensive description of perturbation evolution under shock compression.
Additionally, the Zhang and Sohn 7 model offers an alternative approach to simulate the late-time behaviour of the perturbation, distinguished by its utilization of Padé approximants.However, this model has its limitation as well, notably its tendency to underestimate the late-time average growth rate of the perturbation.To address this issue, Sadot et al. 8 introduced a model developed using data from shock-tube experiments FIG.1: Schematic of the flow configuration in which a perturbation of wavelength λ is applied to an interface separating two phases at pressure p 0 of initial densities ρ h,0 and ρ l,0 shown by the light grey and white regions, respectively, where ρ h,0 > ρ l,0 .A shock is initiated in the heavier phase (dark grey region) with density ρ h,s and pressure p h,s , which travels in the direction of the lighter phase at constant speed U s .
with air and sulfur hexafluoride.Another significant model for the late-time stage is the Dimonte and Schneider 9 model, which, unlike that of Richtmyer 1 , uses a drag coefficient to account for the reduced growth rate at the late-time stage.The predictions obtained from both these models will be compared to those provided by our numerical computations.The Nova laser experiments by Dimonte and Remington 10 , Dimonte et al. 11 and Farley et al. 12 showed reduced growth rates as compared to the impulsive model, particularly at high Mach numbers.Building upon the foundational work in this field, Thornber et al. 13 studied the critical influence of three-dimensional multi-mode initial conditions on the growth of turbulent mixing layers in the RMI.Their use of high-resolution, large-eddy simulations underscored the significant impact of initial conditions on the evolution of these layers.
Further advancing our understanding of the RMI, Thornber et al. 14 also explored the reshock phase, i.e., the period where the initial shock wave is reflected back towards the evolving interface, leading to further, more complex shockinterface interactions.Through high-order accurate Implicit Large-Eddy-Simulations, they extended theoretical models to accurately predict the behaviour of re-shocked mixing layers, offering invaluable insights for scenarios involving broadband initial perturbations.Numerical simulations have also been used to study the effects of viscosity on the RMI.Mikaelian 15 proposed a theoretical model in which viscosity dampens the growth of the instability, while Carlès and Popinet 16 used direct numerical simulations as a benchmark to derive a theoretical model for the RMI of viscous flows.Movahed and Johnsen 17 performed numerical simulations using a solution-adaptive method, which revealed that viscosity dampens the growth of the instability in agreement with the work of Mikaelian 15 .The study of vorticity evolution within RMI has also undergone considerable progress with Kotelnikov, Ray, and Zabusky 18 shedding light on the dynamics of secondary baroclinic vorticity deposition.
Building upon this work, Peng, Zabusky, and Zhang 19 investigated the generation of secondary baroclinic vorticity accelerated by vortices.Their study focused on the effects of this phenomenon on the flow evolution and features of the RMI at various stages, which include amplitude, neck, span, and stream.These aspects are integral to the validation studies presented later in this paper.Baroclinic vorticity has traditionally been viewed as the exclusive source of vorticity evo-lution, often overlooking the impact of viscosity.To investigate the effects of viscosity Wang et al. 20 studied a shockbubble-interaction problem.They introduced a dimensionless number for the baroclinic-to-viscous source ratio in the vorticity evolution.Liu et al. 21demonstrated that viscosity gradients, which emerge due to temperature variations within strong shock waves, make a positive contribution to the generation of vorticity, particularly at high Mach numbers.
Despite the extensive research that has been conducted on RMI and associated phenomena, there are still many unresolved questions, particularly regarding the non-Newtonian effects on the instability.While previous studies have mainly focused on the effects of viscosity on the instability, non-Newtonian fluids exhibit additional complexities, such as shear-thinning and elasticity, which could have a significant impact on the instability growth.This gap in the literature will be addressed in the present work.We investigate the influence of non-Newtonian effects on the generation and damping of vorticity that influence the dynamics of the RMI.By incorporating these effects into our numerical simulations, we aim to provide insights into the dynamics associated with shock impact wherein at least one of the fluids is a non-Newtonian material; this can be relevant to nuclear fusion-type applications such as in the shock-driven fusion application of Derentowicz et al. 22 .
In this study, we use numerical simulations to investigate how shear-thinning properties in the more dense phase (with the lighter one maintained as Newtonian) impact the development of vorticity and subsequently influence the progression of both the RMI and KHI.The non-Newtonian fluid is modelled as a Bird-Carreau-Yasuda fluid; elasto/plastic effects are neglected in this study.We carry out a parametric study to examine the effect of the Mach and Atwood numbers on the emergent dynamics and provide a breakdown of all the terms that contribute to vorticity production and damping to pinpoint the destabilising mechanisms unambiguously over the various stages of the flow; these include the dynamics following the passage of the initial shock as well as subsequent reflected shocks.The results are compared with those obtained when both fluids are either inviscid or Newtonian.
The remainder of this paper is organised as follows.In Sec. 2, we provide details of the problem formulation and the numerical methods used to carry out the computations; a breakdown of the terms contributing to the rate of change of vor-ticity is also included in this section.In Sec. 3, we present the numerical results of our parameteric study, which include comparisons with theoretical predictions, and a discussion of the destabilising mechanisms highlighting the role of shearthinning.Finally, we provide concluding remarks in Sec. 4.

II. PROBLEM FORMULATION
We consider an interface perturbation of wavelength λ separating two fluids of different densities corresponding to a more dense and a less dense fluid, to which we herein refer as 'heavy' and 'light', respectively.The initial densities and viscosities of the light and heavy fluids are ρ l,0 and ρ h,0 , and µ l,0 and µ h,0 , respectively, and both fluids are at a pressure p 0 .As shown in Figure 1, a shock wave resulting from a jump in pressure and density within the heavy phase, propagates towards the lighter phase and impacts the interface, leading to the development of an RMI on the interface.
The pressure and density ratios associated with the shock are denoted by p h,s /p 0 and ρ h,s /ρ h,0 , where the subscripts s and 0 designate the shocked left state (away from the interface) and the right state (close to the interface), respectively; interfacial tension effects are neglected.This assumption is motivated by practical considerations, as seen in the interactions at plastic-gas interfaces in ICF experiments such as those conducted by Derentowicz et al. 22 .The gas phase rheology is considered to have minimal impact on the dynamics as compared to that of the plastic phase, where non-Newtonian effects are more pronounced and relevant.We consider here that the heavy phase can be either inviscid or viscous but with a non-Newtonian rheology featuring a shear rate-dependent apparent viscosity.In contrast, the lighter phase is either considered to be inviscid or Newtonian, with a constant viscosity.

A. Mathematical Model
The shock and interface dynamics are governed by the mass conservation of each phase and the one-fluid formulation of the compressible momentum, and energy equations, respectively expressed by where t, ρ h , ρ l , ρ, µ, u, p, g, e, T , and κ correspond to the time, heavy fluid density, light fluid density, local mixture density and viscosity, mixture velocity, mixture pressure, acceleration due to gravity, mixture total specific energy, temperature, and thermal conductivity, respectively.Using the wavelength of the interface perturbation λ as the length scale and the shock speed U S as the velocity scale, the corresponding time scale is t S = λ /U S ; in addition, we use the dynamic pressure scale ρ h,0 U 2 S , the energy per unit mass scale U 2 S , and the density and viscosity scales as ρ h,0 and µ h,0 , respectively, the above governing equations are written in dimensionless form as follows: where we have scaled T on T h,o − T l,o , the initial temperature difference between the heavy and light phases.Using mass conservation, Eq. ( 7) can also be re-written as In the above equations, C h and C l are the speeds of sound in the heavy and light phases, and Ma = U S /C h , Re = ρ h,0 U S λ /µ h,0 , and Fr = U S / g λ stand for the Mach, Reynolds, and Froude numbers, respectively.Furthermore, ) is the Brinkman number, and ReBr, which characterises the conduction term in Eq. ( 8) can be re-expressed as RePrEc = PeEc where Pe, Pr = C p h,0 µ h,0 /κ, and Ec = U 2 s /C p h,0 (T h,0 − T l,0 ) correspond to the Péclet, Prandtl, and Eckert numbers, respectively, wherein C p h,0 is the specific heat capacity of the heavy phase.In the context of melting plastic in ICF applications, such as those studied by Derentowicz et al. 22 , the high Prandtl number characterises a regime where viscous diffusion is dominant over thermal diffusion, leading us to concentrate on the viscous effects while neglecting conductive heat flow.Thus, in our study, we investigate cases where the influence of conduction is considered negligible corresponding to scenarios where the Péclet number is large, i.e., the Pe → ∞ limit.We also go one step further and consider the isothermal case to focus attention on the non-Newtonian effects on the RMI.
Additionally, the Atwood number, At, is defined as the ratio of the difference of the initial densities of the two phases to their sum, and expressed by: According to this definition, a negative (positive) At suggests the propagation of the shock from the heavier (lighter) to the lighter (heavier) fluid.We will focus below on flows parameterised by negative At.This study is primarily motivated by the impact of Non-Newtonian effects on RMI in ICF, as exemplified by the work of Derentowicz et al. 22 on thermonuclear fusion neutrons, which operates within the negative At range.Equations ( 5)-( 8) are solved using the open-source software blastFOAM 23 , which has been developed to model highlycompressible flows.blastFOAM also employs the volumeof-fluid (VOF) method to capture the interface by advecting the volume fraction of the heavy fluid, α, according to where α = 0 and α = 1 for cells occupied by the light and heavy fluids, respectively; the interface resides in cells wherein 0 < α < 1.Within the VOF framework, α is used to calculate the local density and viscosity as follows: We employ the stiffened gas equation-of-state to calculate the pressure in each phase pi and the mixture pressure p in Mie-Gruneisen form 24 where γ i and ãi are the compressibility and the reference pressures, respectively, which are (constant) material properties in each phase.The non-Newtonian behaviour of the heavier fluid is modelled using the Bird-Carreau-Yasuda (BCY) 25 model in which we take into account the shear-thinning dependence: here, μh,∞ and μh,0 represent the dimensionless infinite and zero-shear rate viscosities, respectively, γ is the dimensionless shear rate, the flow index is denoted by n, while τ represents the ratio of the relaxation time, τ, to the flow time scale, t S .
The EOS parameters are given in Table I as are the rheological properties used for the two phases.
We utilize the advantages of blastFOAM's explicit solution approach for conservative variables that offers both high accuracy and computational efficiency.We use a third-order accurate Runge-Kutta time integration method along with the SFCD interpolation scheme and HLLC flux scheme 26 .We conduct simulations of the RMI within a rectangular domain of size 10 × 1, where unity is the initial wavelength of RMI that is equal to the width of the domain.The shock is initiated in the heavier fluid, on the left side of the interface, as shown in Figure 1, and propagates towards the interface.
In scenarios similar to ICF, the shock travels from the heavier to the lighter fluid 22,27 .Given this, our simulation is designed to emulate conditions relevant to ICF thereby closely reflecting the RMI phenomena we are aiming to study.The top and bottom boundaries are set to be symmetric to the shock axis.We use a zero-gradient boundary condition for density, volume fraction, and pressure at both the left and right boundaries.Moreover, a zero-gradient boundary condition is employed for velocity at the left boundary, while a zero velocity is imposed at the right boundary to simulate a reflected shock.The interface at the start of our simulations is modelled as f (y) = α cos 2πy λ where y is the vertical coordinate, as shown in Figure 1.Furthermore, to mitigate the influence of small-scale perturbations originating from the VOF method and the mesh, a diffusion-type layer is introduced 19 .
We investigate the impact of non-Newtonian effects on the RMI and vorticity evolution for Mach numbers in the range of Ma = [2.5, 10] and Atwood numbers in the range of At = [−0.048,−0.66], respectively.Our focus is on studying how shear-thinning effects on the RMI and vorticity evolution vary with Mach numbers and Atwood numbers.The negative At values reflect the flow situations we will focus on exclusively in the present work, which are associated with shock propagation from the heavier to the lighter phase, typical of ICF applications.The results of these simulations will be compared with those generated for cases wherein both fluids are either inviscid or Newtonian to highlight the role played by the shear-thinning behaviour in the development of instabilities.A discussion of our numerical results is presented next.I.

III. RESULTS AND ANALYSIS
A. Perturbation growth and vorticity generation: inviscid and BCY fluids We begin by describing the development of the RMI at different stages for a situation wherein a shock propagates from the heavier to the lighter fluid (see Figure 2a) characterised by At = −0.66 and Ma = 5.The shock impacts the interface, as depicted in Figure 2b which leads to the development of counter-clockwise vorticity and, in turn, the protrusion of the heavier phase into the lighter one (see Figure 2c).The mechanisms giving rise to vorticity generation and damping, with a particular focus on the role of non-Newtonian effects, are discussed below.Following the propagation of the shock past the interface, the protrusion becomes elongated and the development of roll-up driven by a KHI is clearly evident, giving rise to a classic mushroom-like structure at the protrusion leading edge, as shown in Figure 2d.The shock is reflected by the right boundary and propagates towards the interface from the lighter fluid side of the domain (see Figure 2d).As shown in Figure 2e, upon impact, the reflected shock is repelled by the interface due to the higher density of the heavier fluid, which induces the formation of shock waves that travel towards the upper, lower, and right boundaries.These waves are then reflected from these boundaries and interact once again with the interface, only to be reflected back, and the process is then repeated.As shown in Figure 2f, the mushroom-like structure becomes larger, with a flattened leading edge; highly pronounced KHI-driven roll-up phenomena are also clearly evident in this figure .In Figure 3(a), we plot the temporal variation of the perturbation amplitude for three situations: one in which both phases are inviscid; the second in which both phases are Newtonian; and the third which involves a Newtonian light fluid and a heavier non-Newtonian fluid, the properties of which are given in Table I; the rest of the parameters remain unchanged from those used to generate Figure 2. In the Newtonian-Newtonian case, the impact time between the shock and the interface is denoted by T 1 and is followed in the other cases by the protrusion of the heavier phase into the lighter one at T 2 .The results presented in Figure 3(a) are compared to predictions obtained from the Richtmyer model 1 given by where h0 is the initial amplitude of the perturbation, k is its wavenumber, and ∆ ũ is the velocity jump in the direction of travel of the shock following its propagation through the interface, at which stage this model assumes that the flow becomes incompressible.The Nova laser experiments [10][11][12] showed a lower growth rate at high Ma as compared to the predictions by the Richtmyer incompressible model 1 ; the latter model is only valid for the early, linear stage of the RMI, prior to the development of secondary KHI whose vortices decelerate the RMI growth.At high Ma numbers strong shocks are formed; for this range of Ma, the behaviour of fluid dynamics significantly diverges from simpler, incompressible models.In the case of strong shocks, due to high compression, the transmitted shock recedes slowly and represents an effectively impenetrable boundary, which reduces the growth of the RMI.This phenomenon leads to a complex interaction between the shock waves and the interface perturbations, further complicating the fluid motion and the subsequent development of instabilities.It is therefore unsurprising that the agreement with the Richtmyer model 1 deteriorates after early times giving way to significantly larger growth rate estimates than those generated numerically for both the inviscid and viscous cases.We also compare the numerical predictions for the perturbation amplitude with those obtained from the work of Dimonte and Schneider 9 who introduced a late-stage model which accounts for the drag effects on the growth of RMI.The evolution of the amplitude h of the single-mode RMI is described in this model by the following equation: In Eq. ( 18), c d = 24 Re represents a drag coefficient for simple laminar flow around a sphere.A comparison of the Dimonte and Schneider model 9 predictions and those from the inviscid and viscous numerical models demonstrate an improvement over those obtained from the work of Richtmyer 1 .
We show x − y projections of the interface shape superimposed on vorticity contour plots in Figure 3(b)-(e) for t = 0.075, 0.18, 1.41, 1.48, respectively, for the inviscid and viscous cases; the variation of the z-component of the vorticity, w z , along the arc length is also shown for the same times.It is seen clearly that vorticity is generated in the regions immediately adjacent to the interface, and this is closely correlated to the interfacial deformation associated with the RMI and the KHI-related roll-up phenomena; the development of the latter, whose onset is at T 3 , leads to a deceleration of the perturbation amplitude.Furthermore, the dominant contribution is due to the heavier phase during the earliest stages of the flow; following the crossing of the shock into the lighter phase, the latter becomes the dominant contributor to vorticity generation.
The divergence in the evolving dynamics between the inviscid and viscous cases becomes more accentuated with increasing time as can be confirmed upon inspection of Figures 3(d), (e), (h) and (i); it is also seen that the vortical structures are more pronounced for the inviscid case, for which the interface becomes more stretched as revealed by the expanded arc length range associated with this case (see Figure 3(i)).Lastly, the impact of the reflected shock on the interface at t = T 4 leads to further development of the RMI and KHI and acceleration of perturbation growth.
In Figure 4, we plot the analogous variables to those shown in Figure 3 except Ma = 10, and the rest of the parameters remain unchanged.A comparison of the results presented in these figures reveals similar trends in terms of perturbation amplitude growth and the close association between the vorticity in the vicinity of the interface and the development of the RMI and KHI roll-up phenomena.An exception is provided by the Newtonian-Newtonian case where the strong influence of viscosity prevents any increase in perturbation amplitude, as can be seen from Figures 3(a) and 4(a).At the higher Ma, it is seen clearly that not only are these phenomena significantly more pronounced, the difference in vorticity generated between the inviscid and BCY cases is lower than that associated with the Ma = 5 case.This is due to the higher shear rates at Ma = 10 and therefore lower viscosity in the heavy phase in comparison to that at Ma = 5.In contrast, the Newtonian-Newtonian cases exhibits comparatively very little perturbation growth, even at the elevated Ma.
Following the passage of the incident shock, viscous damping associated with the residual viscosity becomes dominant  and acts to suppress the RMI and KHI roll-up.This is seen to be more effective for the lower Ma studied due to the comparatively smaller shear rates and therefore larger residual viscosity.There is also evidence for perturbation amplitude deceleration due to the development of KHI roll-up in the Ma = 10 case, as can be seen in Figure 4(a).Furthermore, the impact of the first reflected shock at t = T 3 on the KHI and RMI growth rate is similar in both the viscous and inviscid cases, as shown in Figure 4(d).The impact of the second reflected shock at t = T 4 leads to the development of small-scale roll-up structures.These small-scale structures lead to additional vorticity generation, the so-called vortexaccelerated secondary baroclinic vorticity, which has been investigated in detail by Peng et al. 19 .This mechanism involves the interaction of vortices generated by RMI and KHI, causing a further enhancement of baroclinic vorticity deposition, which contributes to an increase in the growth rate of RMI and a more developed KHI.In our simulations, we can observe the effects of vortex-accelerated secondary baroclinic vorticity on the RMI and KHI: the RMI growth rate is higher, and KHI is further developed in the inviscid case.In the BCY case, however, the span of the RMI, which is a low shear rate region, is characterised by a relatively high viscosity that leads to a lower growth rate of the RMI after the second reflected shock.
We now explore the effects of system parameters on the features of the mushroom-like structure discussed above, and show in Figure 5 an annotated schematic of a fully-developed such structure; in this figure, we highlight the dimensions of the various features: the 'span', 'stream', and 'neck'.In Figure 6 we demonstrate the influence of the non-Newtonian effects on these features.As shown in Figure 6(a), there is very little difference in stream development between the inviscid and BCY cases after the first reflected shock for Ma = 10.The second reflected shock leads to the development of smallscale structures, and the occurrence of these structures induces vortex-accelerated vorticity.This additional vorticity in the inviscid case results in a higher stream in comparison to the BCY case wherein the higher viscosity due to the comparatively lower shear rates dampens the small-scale structures, leading to a lower stream.
For Ma = 5, the stream does not develop until the first reflected shock for the BCY case in contrast to the inviscid case.The second reflected shock leads to the development of smallscale structures and vortex-accelerated vorticity, resulting in elevation of the BCY stream, which remains below that of the inviscid case due to larger viscous damping in the BCY case.Similar behaviour is observed in the development of the span and neck for both Ma = 5 and Ma = 10, as shown in Figures 6(b) and (c), respectively.
We have also carried out a parametric study of the effects of At and Ma on the development of RMI and associated phenomena.In Figures 7 and 8 7 for all Atwood and Mach number values investigated.Notably, however, the flow features exhibited by both cases are qualitatively similar due to the shear-induced viscosity reduction in the BCY fluid.This is in marked contrast to the Newtonian-Newtonian case (not shown) which is significantly more stable and is accompanied by virtually no perturbation growth.

B. Mechanisms for vorticity evolution
In the previous sections, we highlighted the significant role of vorticity in the dynamics associated with the KHI and RTI roll-up processes.It was particularly noted that vorticity generation is enhanced with increasing Atwood and Mach numbers.The non-Newtonian effects can slow down the mixing rate within the span of the KHI.In this region, the non-Newtonian effects influence the small-scale structures that play an important role in mixing.For example, in general ICF scenarios such as those studied by Derentowicz et al. 22 , shear-thinning effects can influence the development of smallscale structures, which in turn can slow down the mixing rate of hotter and colder fluids.FIG.9: Orientation of the pressure, density, and viscosity gradients following shock impact and protrusion of the heavier phase into the lighter one.
In this section, we identify the mechanisms that lead to vorticity generation and damping, paying particular attention to the role of non-Newtonian rheology.In the vicinity of an interface between two immiscible fluids of different densities, it is well known that the presence of a pressure gradient as a result of hydrostatic pressure or a shock leads to vorticity generation in RTIs and RMIs, respectively.However, when considering generalised Newtonian fluids, not only should the viscous dissipation be taken into account in the vorticity balance, but also the viscosity gradients.Taking the curl of Eq. ( 7), the two-dimensional vorticity equation can be written in the one-fluid formulation framework, expressed by Eqs. ( 5)-(8) (details are in Appendix A): Note that in the simultaneous limits Ma → 0 and Re → ∞, this equation reduces to the greatly simplified form of the vorticity equation associated with the inviscid case 28 Dω Dt = 1 which implies that the baroclinic term on the right-hand-side of this equation is responsible for vorticity generation resulting from misalignments of density and pressure at the interface.We illustrate in Figure 9 how the major gradients such TABLE II: A breakdown of the rate of change of vorticity into the terms that appear on the right-hand-side of equation (19).The terms are characterised as 'Stabilising' or 'Destabilising' depending on their time-averaged contribution to this equation, and are ranked according to the magnitude of their respective contributions, e.g, 'S1' and 'D1' correspond to the terms with the largest stabilising and destabilising contributions, respectively.Note: Some terms, like Term 5,6,9 and 10 have mixed effects (both stabilising and destabilising at different times), but are classified by their dominant time-averaged effect.
as ∇ρ, ∇µ, and ∇p are oriented after the passage of the shock and the creation of the protrusion of the heavier phase into the lighter one.This protrusion leads to large gradients in pressure, density, and viscosity across the interface, transitioning from high values in the heavier fluid to lower ones in the lighter fluid.The misalignment of ∇p and ∇ρ leads to a destabilising baroclinic mechanism, which is present in the inviscid case, that drives the RMI.For finite Re, other effects compete with this mechanism.The presence of non-Newtonian rheology in one of the phases gives rise to additional mechanisms related to ∇µ.
We list the various terms in Table II, Terms 1-12 that appear in the dimensionless vorticity Eq. ( 19) derived from the dimensionless momentum equations; time is nondimensionalised as specified in the earlier section.These dimensionless terms are integrated over half of the dimensionless wavelength for volume fractions α = 0.1 and α = 0.9 yielding integrals I 1 − I 12 , respectively, and the temporal evolution of these integrals is then plotted in Figures 10, 11 , respectively.We also obtain time-averages of I 1 − I 12 , which highlight the overall contributions of Terms 1-12 to the system stability over the duration of the flow investigated, and this is also indicated in Table II.We now examine the contribution of each term.Inspection of Figure 10(a) reveals that the primary shock impact leads to a rapid rise in I 3 , I 4 , and I 11 which exhibit positive values highlighting the destabilising contributions of Terms 3, 4, and 11; in contrast, I 8 is comparatively smaller in magnitude though it is also destabilising albeit marginally.It is also seen that I 4 associated with the inviscid case is virtually indistinguishable from that of the BCY case.The largest instantaneous contributor to instability is Term 4, which is related to the baroclinic mechanism discussed above, and this is followed by Term 11, represented by − 1 Re 1 ρ ∇µ ×∇×ω, a distinct source of vorticity generation specifically driven by the gradient in viscosity.This, in turn, is followed by the contribution from Term 3, which is represented by

Dp
Dt ∇µ, is marginally destabilising and signifies the interactions between the density and viscosity gradients, the temporal variation in pressure, and the Mach number, which provides a measure of the flow compressibility.
The ordering of the terms remains unchanged following the impact of the reflected shock, though, as shown in Figure 10(b), the magnitudes of I 3 , I 4 , I 8 , and I 11 are smaller than their early-times counterparts.We also note that whereas the contribution of the baroclinic term is short-lived, as can be dis-cerned from Figure 10(a), that of Term 11 is relatively longerlasting, and the time-averages of I 3 , I 4 , I 8 , and I 11 indicate that the largest contributor to instability is associated with Term 11 (see Table II), which exceeds Term 4 and further highlights the role of viscosity gradients, related to the non-Newtonian rheology of the heavier fluid, in driving instability.Lastly, in the insets in Fig. 10, we show x − y projections of the density field associated with the inviscid and BCY cases following the impact of the reflected shock; inspection of these Figures reveals that the interface is unstable due to the shock impact in the BCY case, albeit not to the extent of the instability observed in the absence of viscous effects.
In a similar vain to Figure 10, we show in Figure 11 the temporal evolution of the spatial integrals of Terms 5, 6, 9, and 10, I 5 , I 6 , I 9 , and I 10 , respectively, over the same time duration and parameter values.It is clearly seen that Terms 5 and 6, which are represented by 1 Re µ∇ × ω, are destabilising and stabilising, respectively, following the impact of the primary shock, as shown in Figure 11(a); the impact of the reflected shock appears to have virtually no effect on I 5 , while I 6 is slightly negative.This suggests that Terms 5 and 6 have a weakly stabilising influence on the late-stage dynamics.Although Term 9, represented by − 1 vorticity generation.The reflected shock impact also leads to an increase in I 10 .The overall contribution to the flow, however, is destabilising due to Term 5, and stabilising due to Terms 6, 9, and 10.In the case of the latter, it appears that its stabilising long-time influence outweighs its destabilising effect during the shock impact episodes at the early-and latestage dynamics.
In Figure 12, we assess the contributions due to Terms 1, 2, 7, and 12.  damping, vorticity-attenuating effect of Terms 1, 2, 7, and 12 is provided by the temporal evolution of the spatial integrals of these terms, which are shown in Figure 12 as well as their respective time-averages as highlighted in Table II.
In Figure 13, we show the temporal variation of the major vorticity-generating terms, Terms 4 and 11, integrated along half of the arc length, I (s) 4 and I (s) 11 , respectively.This analysis is performed when the shock wave interacts with the interface at Mach numbers 2.5, 5, and 10, with At = −0.66,both in the inviscid case and in situations where the denser and lighter phases are represented by BCY and Newtonian fluids, respectively; results were generated for α = 0.1 and α = 0.9.It is seen that the magnitude of Term 4 associated with inviscid flow is larger than that for the viscous case, as expected, for all Ma considered, and higher in the heavier than the lighter fluid.It can also be observed that although both I increases.This indicates that the baroclinic mechanism driving vorticity growth gains in significance with increasing Ma relative to the one associated with viscosity gradients arising from non-Newtonian effects in the heavier phase.

IV. CONCLUSIONS
We investigated the effects of non-Newtonian rheology on the development of a Richtmyer-Meshkov instability (RMI) in flows characterised by negative Atwood numbers wherein a shock is initiated in the heavier, non-Newtonian fluid, and travels towards the lighter, Newtonian one.In order to resolve the shock and interfacial dynamics, a one-fluid formulation for the mass, momentum, and energy equations for a compressible flow was used.A shear-thinning rheology was considered for the heavier phase, and a stiffened gas equationof-state was used to relate the pressure to the density.
A volume-of-fluid approach was utilised to capture the interface, and the open-source software blastFOAM was employed to obtain numerical solutions of the governing equations for a range of Atwood and Mach numbers and a fixed set of rheological parameters for the heavier phase.We also derived a two-dimensional vorticity transport equation, and using a budget analysis of this equation, we identified terms that generate and dampen vorticity, which provide valuable insights into the interplay between non-Newtonian effects, vorticity generation, and RMI growth in various regimes.
Our results demonstrated that the incident shock creates a viscosity gradient along the interface, which generates additional vorticity due to non-Newtonian effects compared to the inviscid case, which was used as a benchmark.This additional vorticity generation is more pronounced at moderate Mach numbers Ma (< Ma = 5) than at high Mach numbers (Ma = 10) because lower shear rates in the former Ma range result in a higher residual viscosity gradient.
After the passage of the shock, the damping effects of viscosity become dominant, leading to a reduction in the vorticity.The reflected shock then generates additional vorticity due to the non-Newtonian viscosity gradient.For the low-to-moderate Ma range, the damping effects of viscosity hinder the development of secondary instabilities, such as the Kelvin-Helmholtz instability (KHI), leading to lower RMI growth.This observation might seem counterintuitive.One might intuitively assume that suppressing secondary instabilities would lead to enhanced growth of RMI, given that such instabilities act as a mechanism to limit or saturate growth.However the dampening effects of viscosity are not only present in nonlinear phase, where the viscosity effects dampen the onset of KHI at low-moderate Mach numbers.The dampening effects of viscosity also play a crucial role in the linear phase.At high Ma, the residual viscosity is low due to the very high shear rates associated with high-Ma flows, resulting in low rates of vorticity generation associated with non-Newtonian effects.Furthermore, viscous effects in these situations damp out the small-scale structure, which could have given rise to vortex-accelerated vorticity deposition that would have enhanced the RMI and associated mixing.
Although neglected in the present work, incorporating temperature dependence into the viscosity model aligns the growth rate more closely with the inviscid scenario, attributed to the substantial reduction of viscosity at elevated temperatures.The significance of viscosity becomes especially apparent within the regime span, marked by the presence of smallscale structures.These structures are instrumental in augmenting the generation of baroclinic vorticity through vortex acceleration, a phenomenon thoroughly investigated by Peng et al. 19 .A reduced occurrence of these small-scale structures could lead to reduced mixing, presenting a potential area for further exploration as an extension of the present work.

FIG. 2 :
FIG. 2: Projections of the density field in the x − y plane depicting the development of a RMI in the inviscid case, and a subsequent KHI, driven by the left-to-right propagation of a shock wave from the heavier towards the lighter phase for t = 0, 0.06, 0.25, 1.35, 1.42, 1.80 shown in (a)-(f), respectively.The parameter values are At = −0.66,Ma = 5, and Fr = 48571; the remaining parameters are in TableI.

FIG. 3 :
FIG. 3: Temporal variation of the perturbation amplitude for At = −0.66 and Ma = 5, (a), where the black solid and dashed curves represent the predictions of the Richtmyer 1 and Dimonte and Schneider 9 models, the blue curves represent the inviscid case, while the red dashed and solid lines represent the numerical predictions obtained for the Newtonian-Newtonian and BCY cases, respectively; in the latter case, the light and heavy fluids correspond to Newtonian and BCY fluids, respectively; here, T 1 − T 4 correspond to t = 0.075, 0.18,1.41,and 1.48, respectively.Snapshots of the x − y projections of the density field at t = T 1 − T 4 depicting the interface shape superimposed on vorticity contours shown in (b)-(e), respectively; in each panel the predictions associated with the BCY and inviscid cases are shown on the left and right, respectively.Variation of the z-component of the vorticity with arc length for the inviscid and viscous cases represented by blue and red lines for t = T 1 − T 4 shown in (f)-(i), respectively, plotted on the heavy-(with α = 0.9) and light-side (with α = 0.1) of the interface using solid and dashed lines, respectively.The rest of the parameters remain unchanged from the previous Figure.

FIG. 4 :
FIG. 4: Temporal variation of the perturbation amplitude for At = −0.66 and Ma = 10, (a), where the black solid and dashed lines represent the predictions of the Richtmyer 1 and Dimonte and Schneider 9 models, the blue curves represent the inviscid case, while the red dashed and solid lines represent the numerical predictions obtained for the Newtonian-Newtonian and BCY cases, respectively; in the latter case, the light and heavy fluids correspond to Newtonian and BCY fluids, respectively; here, T 1 − T 4 correspond to t = 0.037, 0.085, 0.815, and 1.175, respectively.Snapshots of the x − y projections of the density field at t = T 1 − T 4 depicting the interface shape superimposed on vorticity contours shown in (b)-(e), respectively; in each panel the predictions associated with the BCY and inviscid cases are shown on the left and right, respectively.Variation of the z-component of the vorticity with arc length for the inviscid and viscous cases represented by blue and red lines for t = T 1 − T 4 shown in (f)-(i), respectively, plotted on the heavy-(with α = 0.9) and light-side (with α = 0.1) of the interface using solid and dashed lines, respectively.The rest of the parameters remain unchanged from the previous Figure.
FIG.5: Schematic of a fully-developed mushroom-like structure formed following the passage of a shock travelling from a heavy phase towards a lighter one through the separating interface.The features stream, span, neck and amplitude19 characterise the structure are shown together with their defined dimensions.

Figure 8
Figure 8 highlights the combined effects of At, Ma, and the non-Newtonian rheology in the heavy phase on RMI development.Inspection of this Figure reveals that viscous damping decelerates the development of the RMI in comparison to the inviscid case shown in Figure7for all Atwood and Mach number values investigated.Notably, however, the flow features exhibited by both cases are qualitatively similar due to the shear-induced viscosity reduction in the BCY fluid.This is in marked contrast to the Newtonian-Newtonian case (not shown) which is significantly more stable and is accompanied by virtually no perturbation growth.

FIG. 6 :
FIG. 6: Temporal evolution of the length of the stream, span, and neck (see Figure 5 for the definitions of these features) At = −0.66shown in (a)-(c), respectively.In each panel, the red and blue solid lines designate the inviscid and BCY cases for Ma = 10, and their dashed counterparts the Ma = 5 case, respectively.The solid and dashed vertical lines at t = 0.71 and t = 0.95, and t = 1.41 and t = 1.89 highlight the times at which the interface is impacted by the first and second reflected shocks for Ma = 10 and Ma = 5, respectively.The rest of the parameters are unchanged from previous Figures.
FIG. 7: Fully-developed x − y projections of the density field for the inviscid case in the (At, Ma) space for At = (−0.048,−0.230, −0.66) and Ma = (5, 7.5, 10); the rest of the parameters remain unchanged from previous Figures.

8 :
FIG. 8: Fully-developed x − y projections of the density field for the case wherein the lighter and heavier phases correspond to Newtonian and BCY fluids, respectively, in the (At, Ma) space for At = (−0.048,−0.230, −0.66) and Ma = (5, 7.5, 10); the rest of the parameters remain unchanged from previous Figures.

FIG. 10 :
FIG. 10: Temporal evolution of the spatial integrals of Terms 3, 4, 8, and 11, I 3 , I 4 , I 8 , and I 11 shown between t = 0.06 − 0.16 and t = 1.4 − 1.5 in (a) and (b), respectively, highlighting the effect of impact by the primary shock at t ≈ 0.7 in (a), and reflected shock at t ≈ 1.415 in (b), for At = −0.66 and Ma = 5.The green solid and dotted lines represent I 4 and I 11 , the solid and dashed khaki lines represent I 3 and I 8 for there case wherein the heavier and lighter phases correspond to BCY and Newtonian fluids, respectively.The dashed black line represents I 4 for the case wherein both fluids are treated as inviscid.The insets in panel (b) represent x − y projections of the density field for the inviscid and BCY cases following the impact of the reflected shock.All other parameters mirror those from preceding Figures.

1 1 ρFIG. 11 :FIG. 12 :
FIG. 11: Temporal evolution of the spatial integrals of Terms 5, 6, 9, and 10, I 5 , I 6 , I 9 , and I 10 shown between t = 0.06 − 0.16 and t = 1.4 − 1.5 in (a) and (b), respectively, highlighting the effect of impact by the primary shock at t ≈ 0.7 in (a), and reflected shock at t ≈ 1.415 in (b), for At = −0.66 and Ma = 5.The blue solid and dashed lines represent I 5 and I 6 , and the gray solid and dashed lines represent I 10 and I 9 , respectively.All other parameters mirror those from preceding Figures.

1 1 ρFIG. 13 :
FIG. 13: Variation of the arc length integrals of terms 4 and 11, I 11 , respectively, when the primary shock impacts the interface with At = -0.66 for Ma = 2.5, Ma = 5, and Ma = 10 shown in (a) and (b), (c) and (d), and (e) and (f), respectively.The definition of Terms 4 and 11 is provided in Table II.The blue and red lines represent the inviscid case and that wherein the heavier and lighter phases correspond to BCY and Newtonian fluids, respectively.The solid and dashed lines represent the results associated with α = 0.9 and α = 0.1, respectively.
Ma, the ratio of the magnitudes of I

TABLE I :
Dimensionless stiffened gas equation of state parameters which appear in equation(14).The initial densities are scaled on ρ h,0 .Rheological parameters for the light and heavy phases which appear in equations (16). 1