Two-dimensional and three-dimensional computational fluid dynamics studies of a spherical bubble impacted by a supersonic shock wave (Mach 1.25) have been performed to fully understand the complex process involved in shock–bubble interaction (SBI). The unsteady Reynolds-averaged Navier–Stokes computational approach with a coupled level set and volume of fluid method has been employed in the present study. The predicted velocities of refracted wave, transmitted wave, upstream interface, downstream interface, jet, and vortex ring agree very well with the relevant available experimental data. The predicted non-dimensional bubble and vortex velocities are also in much better agreement with the experiment data than values computed from a simple model of shock-induced Rayleigh–Taylor instability (the Richtmyer–Meshkov instability). Comprehensive flow visualization has been presented and analyzed to elucidate the SBI process from the beginning of bubble compression (continuous reflection and refraction of the acoustic wave fronts as well as the location of the incident, refracted and transmitted waves at the bubble compression stage) up to the formation of vortex rings as well as the production and distribution of vorticity. Furthermore, it is demonstrated that turbulence is generated with some small flow structures formed and more intensive mixing, i.e., turbulent mixing of helium with air starts to develop at the later stage of SBI.
I. INTRODUCTION
The flow field created by shock interaction with a geometrically distinct density inhomogeneity entails the intense coupling of various types of fluid dynamic phenomena, such as shock wave refraction and reflection, generation and transport of vorticity, as well as turbulence mixing. This interaction leads to a complicated configuration of shock and rarefaction waves mainly through focusing and scattering with the concurrent production of characteristic vortices and frequent mixing of the surrounding gas. This physical phenomenon is essential for the studies and investigation of complex challenges related to shock travel through an arbitrary material with varying density, temperature, and other thermodynamic state variables. Ranjan et al.1 suggested that the simplest configuration through which all these processes can be studied in detail is the interaction of a shock wave with a cylindrical or spherical bubble and this study is focused on the interaction of a shock wave with a spherical bubble.
Considering the specific situation where a planar incident shock wave interacts with a discrete sharply well-defined spherical gas bubble surrounded by ambient gas, assuming that the surrounding un-shocked ambient gas has a density and speed of sound of ρ1 and c1 while the density and speed of sound of the un-shocked bubble gas are ρ2 and c2, a parameter called the Atwood number, A, is defined as2
From Eq. (1), two scenarios may be defined corresponding to the light bubble case where A < 0 and the heavy bubble case where A > 0. The present study focuses on the light bubble case.
Early experimental study on the interaction of a planar shock wave with light- or heavy-gas spherical or cylindrical inhomogeneities was carried out by Rudinger and Somers,3 and they also introduced a simple model to calculate the relative velocity of the shocked bubble. Experimental studies of Haas and Sturtevant4 and Jacobs5 showed that vorticity is mainly generated on the bubble gas/ambient gas interface by the baroclinic mechanism (measure of misalignment between the pressure gradient of the interacting shock waves and density gradient in the bubble gas) throughout the propagation of the shock wave, which leads to the overturning of the interface’s upper section and the subsequent production of a vortex ring. This process of vorticity production is analogous to the Richtmyer–Meshkov instability.6 The shock refraction pattern is determined by the impedance mismatch at the interface, which exists in the early phase of the flow. Subsequent to several shock-passages through the bubble, the characteristics noticed in the flow field are dominated by the vortical movement. Layes et al.7,8 conducted experimental studies on the interaction of a planar shock wave with a spherical gas inhomogeneity and provided qualitative characterization of the flow evolution for SBIs. Shock–bubble interactions at a higher Mach number of 2.88 than those studied by Haas and Sturtevant4 and Layes et al.7,8 were studied experimentally by Ranjan et al.9 and distinct secondary vortex rings formed at later times, which had not been observed previously in experiments at lower Mach numbers.
Experimental studies have limitations in providing detailed spatial/temporal information due to restriction by instrumentation and experimental conditions, etc. Hence, there have been many numerical/theoretical studies carried out to study SBIs. Picone and Boris10 presented a new theoretical model that showed that after a shock passed a bubble, the formation of vortex structures and continual distortion of the bubble were due to the generation of long-lived vorticity at the edge of the discrete inhomogeneity. They also carried out a 2D numerical study on the interaction of a weak shock with a bubble to verify their theoretical model, and the simulated results were similar to those observed in the experiments of Haas and Sturtevant.4 Klein et al.11,12 carried out 2D and 3D experimental and numerical studies of the interaction between a planar shock wave and a spherical inhomogeneity. Their study provided a good description of the shocked cloud: the size and shape of the cloud, the mean density, the mean velocity, the total circulation, and so on. Quirk and Karni13 conducted a detailed numerical study of shock–bubble interaction (SBI) and presented the shock refraction pattern plus vortical features observed by Haas and Sturtevant.4 Samtaney and Zabusky14 described a detailed analysis of the baroclinic vorticity deposition on a planar interface utilizing shock polar analysis. This finding was applied to the spherical scenario using a near-normality condition. Their model revealed that for large density ratio across the interface, the circulation does not rely on the density ratio. In addition, for high Ma number flows, they showed that circulation is directly proportional to Ma number. Zhai et al.15 carried out a 2D numerical study of the interaction between a planar shock wave and a spherical bubble. Detailed flow field structures were obtained in their study, and it was found that the generation and distribution of vorticity were the dominant factors for the interface deformation and turbulent mixing. Levy et al.16 applied an interface-tracking 2D arbitrary Lagrangian-Eulerian (ALE) hydrodynamic code to simulate a SBI, and their numerical model was an extension of Samteney and Zabusky’s14 circulation model to the velocity field scaling, which revealed that the bubble velocity does not rely on the radius of the bubble and that the velocity scaling failed for M > 2. Zhu et al.17 investigated effects of the Atwood number (At) on the evolution of the shock wave and gas bubble through 2D numerical simulations. Their main conclusion is that the Atwood number has a non-monotonic influence on the evolution of mixedness, average vorticity, and circulation.
The above numerical studies are mainly 2D, and so far not many full 3D numerical studies have been carried out. Niederhaus et al.2 carried out comprehensive 3D Euler simulations of SBI across a wide range of Mach and Atwood numbers. They proposed a new model for the total circulation and showed primary/secondary vorticity production in sequenced visualizations of the density and vorticity fields. Their work also revealed that the mean density and bubble velocity could be scaled over a broad range of Mach numbers at a fixed Atwood number using quantities predetermined from the one-dimensional gas dynamics. Yang et al.18 performed a 3D numerical study on the interaction of a shock wave with a spherical helium bubble at different Mach numbers. Their results showed that at a fixed Atwood number, the compression ratio is greatly amplified with increasing Mach numbers and the time-dependent mean bubble velocity is also affected. Hejazialhosseini et al.19 conducted high resolution 3D simulations of shock–bubble interaction for several Mach and Atwood numbers. Their results show the evolution of the flow structures and demonstrate that a number of analytical models cannot predict accurately the late time circulation of shock–bubble interactions at high Mach numbers. Rybakin and Goryachev20 studied the process of a shock wave interaction with low-density gas bubble through a 3D Euler simulation, visualizing the complex SBI process of formation of a system of reflected, refracted, and passed waves, as well as the formation of vortex rings.
Despite many previous experimental and numerical studies, our current understanding of the SBI is far from complete. In the present study, a detailed 3D simulation of SBI has been carried out to advance our current understanding of the complex shock/bubble interaction process. In particular, as far as we know that no previous studies have addressed the generation and development of turbulence at the later stage of SBI, which will be elucidated in the present study. Furthermore, as the SBI is inherently 3D in many aspects and hence it is necessary to perform 3D simulations, but so far a majority of numerical studies are 2D due mainly to its computational efficiency. However, there is always a question mark on the accuracy of 2D simulations. Hence, a 2D simulation of SBI has also been conducted in the present study to assess the accuracy or inaccuracy of 2D simulation in such a case through a direct comparison against the 3D simulation.
This paper is structured as follows: methodology, including grid independence study, is described in Sec. II. Section III presents quantitative and qualitative comparison between the predictions and experimental data, plus comprehensive analysis of the instantaneous flow fields to elucidate the process of SBI. Concluding remarks are drawn in Sec. IV.
II. METHODOLOGY
A. Governing equations and numerical method
The fundamental governing equations for fluid flows, called Navier–Stokes equations, are based on fundamental conservation laws for mass, momentum, and energy, which can be written in a generic form as follows:
and
where ρ, p, and E denote the density, pressure, and total energy, respectively, while u,v, and w represent the velocity component in the x, y, and z directions, respectively. H, which is zero in the continuity equation, represents viscous force term in the momentum equation and viscous work term in the energy equation, which will not be presented here as those are standard terms given in many textbooks.
SBI is predominantly unsteady, and turbulence is usually generated after the interaction at the later stage of SBI so that the Unsteady Reynolds-Averaged Navier–Stokes (URANS) approach is adopted in the present study. Two other more accurate numerical approaches (large-eddy simulation and direct numerical simulation) for simulating turbulent flows have not been selected in the present study mainly because that turbulence is generated only at the later stage of SBI and the flow consists of mainly unsteady large scale flow structures that URANS can capture very well at a much lower computational cost. The URANS equations are derived by averaging the above instantaneous Navier–Stokes equations, and the averaging procedure leads to extra terms called Reynolds stresses that need to be modeled using a turbulence model. There have been many turbulence models developed, and the selection of a suitable turbulence model in the present study will be given in Sec. II D.
It is crucial to capture the bubble deformation in SBI studies, and hence an interface tracking method is needed. There are many methods available to track the interface, but some of those methods could fail when handling large interface deformations, such as disintegration and fusing of fronts.21 It has been demonstrated that a Coupled Level Set and Volume of Fluids (CLSVOF) model boosts a higher accuracy than the standalone level set or volume of fluid method as it combines their merits and overcomes their deficiencies (Olsson and Kreiss).22 The CLSVOF approach has proven to be a potent tool to solve huge interface deformations and was applied successfully by Niederhaus et al.2 to study SBI. Therefore, this approach has been employed in the present study.
The computer code used is ANSYS Fluent 19.0, and a third-order Monotone Upstream-centered Schemes for Conservation Laws (MUSCL) scheme is used for the spatial discretization. The scheme is selected as it has the capabilities to boost spatial accuracy for all grid types by decreasing numerical diffusion particularly for complex 3D flows. A first-order implicit scheme is employed for temporal discretization.
B. Computational details
The computational setup is based on the experiment carried out by Haas and Sturtevant4 with the main parameters listed in Table I.
Bubble gas . | Ambient gas . | Ma number . | Atwood number . |
---|---|---|---|
Helium | Air | 1.25 | −0.715 |
Bubble gas . | Ambient gas . | Ma number . | Atwood number . |
---|---|---|---|
Helium | Air | 1.25 | −0.715 |
Figures 1(a) and 1(b) present the computational domains and the boundary/initial conditions for the 2D (half of the domain/axisymmetric model) and 3D cases. The flow is assumed axis-symmetrical in the 2D case with the horizontal x-axis is in the direction of the shock tube axis and the y-axis in the vertical direction. The domain length in the 2D case is Lu + 0.44d (d = 45 mm is the bubble diameter), where Lu is variable to ensure enough room that will allow the shocked bubble travel throughout the monitored time. The domain height is 0.99d, as shown in Fig. 1(a). A no-slip wall boundary condition is applied at upper boundary (wall of shock tube). The computational setups are matching those in the experiment by Haas and Sturtevant4 with the right boundary being the inlet as the incident shock propagated from right to left in the experiment. At inlet, density, pressure, and velocity are specified, and at outlet, a zero gradient boundary condition is applied.
The computational domain length in the 3D case is the same as that in the 2D case while the domain height and width are both 1.98d, as shown in Fig. 1(b). The same boundary conditions used in the 2D case are applied in the 3D case, and a no-slip wall boundary condition is applied on the side walls.
An unstructured mesh is employed to capture the bubble surface accurately, which would pose a great challenge if a structured mesh were used. For the 2D case, a hybrid mesh having a structured uniform grid in most of the region and an unstructured grid in the vicinity of the bubble is generated, as shown in Fig. 2(a). For the 3D case, a more or less uniform tetrahedral volume mesh is generated in most of the flow region, and a triangular surface mesh is employed at the gas bubble boundary, as shown in Fig. 2(b). Adaptive mesh refinement method has been used in the present study to concentrate more cells in the bubble, its vicinity, and the primary shock wave, guaranteeing that fine cells surround the bubble to support an effective simulation of the strong interaction/mixing of the planar shock wave and the bubble gas.
C. Mesh independence study
A mesh independence study has been carried out with three meshes for both 2D case (0.16, 0.24, and 0.28 × 106 cells) and 3D case (3.6, 6.2, and 9.2 × 106 cells). Figure 3 presents the bubble compression and changes against time in terms of three representative points (upstream, jet, and downstream locations as indicated in the figure), and it can be seen that for the 2D case, the results are hardly changing when the mesh is refined from 0.24 to 0.28 × 106 cells, and hence there is need to refine the mesh further, and the rest of the 2D simulations have been performed using 0.28 × 106 cells. Similarly for the 3D case, the results obtained using the mesh with 6.2 × 106 cells are very close to the results obtained using the mesh with 9.2 × 106 cells, and hence there is no need to refine the mesh further, and the rest of the 3D simulations have been carried out using 9.2 × 106 cells.
D. Turbulence model selection
There are many turbulence models, and there is no general consensus which model is the best as their performances are very much case dependent. There is hardly any knowledge gained so far about the performance of turbulence models in the simulation of SBI. Hence in the present study, three widely used and high rated turbulence models, the realizable k-ε, the shear stress transport (SST) k-ω, and a Reynolds Stress Model (RSM), have been tested to assess their performance in the 3D case. Figure 4 presents the comparison between the predicted three location changes of the bubble against the experimental data of Haas and Sturtevant.4 It is evident from Fig. 4 that the predictions using all three turbulence models are very close to each other and the results obtained using the realizable k-ε are slightly closer to the experimental data. Therefore, the realizable k-ε is employed in the present study.
III. RESULTS AND DISCUSSION
A. Comparison between the measured and predicted velocities
Table II presents the comparison between the measured and predicted velocities of incident (vI), refracted (vR) and transmitted waves (vT), initial upstream interface (vIU), final upstream interface (vFU), initial downstream interface (vID), air-jet head (vAJ), and vortex ring (vVR). The average velocities of these acoustic wave characteristics are derived by utilizing minimum squares line adjustments to approximate their approximate locations, at different time, along the advected bubble trajectory. Generally, the calculation of these velocities entails the measurement of the changes in their respective positions over 50-time steps and dividing by the elapsed duration, i.e., . As such, the values provided in Table II have been computed by averaging over a number of such intervals.
. | vI . | vR . | vT . | vIU . | vFU . | vID . | vAJ . | vVR . |
---|---|---|---|---|---|---|---|---|
Experimental data | 420 | 960 | 365 | 190 | 125 | 145 | 335 | 165 |
Predictions (2D) | 406.7 | 932.1 | 351.5 | 178.7 | 115 | 134.6 | 318.9 | 155.8 |
Predictions (3D) | 416 | 951.4 | 361.3 | 184.7 | 122 | 142.4 | 331.1 | 161.8 |
. | vI . | vR . | vT . | vIU . | vFU . | vID . | vAJ . | vVR . |
---|---|---|---|---|---|---|---|---|
Experimental data | 420 | 960 | 365 | 190 | 125 | 145 | 335 | 165 |
Predictions (2D) | 406.7 | 932.1 | 351.5 | 178.7 | 115 | 134.6 | 318.9 | 155.8 |
Predictions (3D) | 416 | 951.4 | 361.3 | 184.7 | 122 | 142.4 | 331.1 | 161.8 |
As can be seen from Table II that an excellent agreement between the 3D predictions and experimental data has been obtained for all velocities with the maximum error being 2.8% for the initial upstream interface velocity. The 2D predictions also agree well with the measured values with the maximum error being 8% for the final upstream interface velocity.
B. Bubble acceleration and vortex formation
Air is accelerated in the shock tube as the shock wave propagates from a state of rest to a uniform velocity, . As helium is less dense than air, the major flow velocity trailing the incident shock functions as a piston on the bubble, thus contributing to bubble deformation and compression.23 The sphere will then be accelerated to a higher velocity than , which means that the light helium gas travels ahead of the ambient air as it is transported down the shock tube. Rudinger and Somers3 developed a basic two-stage model for this: in the first stage and during the early transient, the helium bubble accelerates as a solid body to velocity ; in the second stage, it changes into a vortex ring with velocity based on the Taylor mechanism.24 Within the first stage, the traveling undistorted bubble acts as Taylor’s “dissolved” vortex-generating disk and as such, , denotes the disk velocity. Rudinger and Somers3 proposed that in the first stage, the impulse per unit volume experienced by the bubble, Ib would be the same as that underwent by the ambient air, i.e., product of ρair and . This is given mathematically in the following equation:
k denotes the apparent mass fraction for a sphere (=0.5), and ρb is the gas density in the bubble. From Eq. (4), the initial non-dimensional velocity of the bubble is then computed as
where σ = ρb/ρair.
The conversion of the helium bubble into a vortex indicates a drop in the relative velocity. This is shown mathematically as
β has a numerical value of 0.436. Similarly, from Eq. (6), the non-dimensional vortex velocity can be computed as
The and obtained from Eqs. (5) and (7), as well as the predicted and measured values, are presented in Table III.
. | . | . |
---|---|---|
Theoretical | 2.199 | 1.52 |
Experimental data | 1.3 | 1.28 |
Predictions (2D) | 1.177 | 1.181 |
Predictions (3D) | 1.26 | 1.246 |
. | . | . |
---|---|---|
Theoretical | 2.199 | 1.52 |
Experimental data | 1.3 | 1.28 |
Predictions (2D) | 1.177 | 1.181 |
Predictions (3D) | 1.26 | 1.246 |
It can be seen from Table III that the predictions and the experimental data agree very well, especially that the 3D predictions are very close to the experimental data while the theoretical values are much larger with the non-dimensional bubble velocity being almost twice the prediction and experimental values. This strongly indicates that the two-stage model proposed by Rudinger and Somers3 does not quite represent the real situation of bubble acceleration and vortex formation in the process of SBI.
Further comparison of the 2D and 3D predictions against the experimental data4 is shown in Fig. 5, which presents the predicted and measured three location changes of the bubble against the time. It can be clearly seen that the predictions agree well with the experimental data, and as expected, better agreement has been obtained between the 3D predictions and the experimental data.
It is evident from the above quantitative comparison between the predictions and experimental data that the present numerical simulations, especially the 3D simulations, have captured the complex process of SBI very well. Hence, all the analysis below will be based on the 3D results.
C. Visualization of the shock–bubble interaction process
The intricate and key processes involved in bubble distortion, vorticity production, air-jet development, and the overall late-time progression of the SBI process will be further analyzed in this section. Figure 6 presents six snapshots of the simulated images (left) and shadow-photographs (right) taken in the experiment. Time, t, is normalized using velocity of the shock and the bubble diameter.
Overall, the predictions capture the SBI process, especially the key features, very well as the simulated images are quite similar to the shadow-photographs at all those six times. Figure 6(a) shows the incident shock wave shortly after its collision with the helium bubble, and the refracted wave, traveling ahead of the incident shock, can be clearly seen from the simulated image while it is barely observable from the shadow-photograph. The shadow-photograph in Fig. 6(b) reveals the spherical transmitted shock wave at the left trailed by a flat front that represents the torus-shaped secondary transmitted wave’s projection. All these features are well captured in the simulation as shown in the simulated image [left side of Fig. 6(b)]. Figure 6(b) also shows that the bubble is flattened in the direction of shock travel. As time progresses, a developing air-jet has started to penetrate the upstream interface of the helium volume, which appears to be sinking inwards as seen from Fig. 6(c), thus making the distorted bubble appear kidney shaped due to the re-entrant jet’s formation at the upstream end. The re-entrant jet is formed due the impact of the ambient air on the upstream interface driven by the generated vorticity17 and pierces much further into the bubble as shown in Fig. 6(d), accompanied by back-scattered waves to the right, which are more clearly observable in the simulated images. In the meantime, diffracted waves are observed in both the simulated images and shadow-photographs, more obvious in the shadow-photographs. During bubble deformation and the formation of the air-jet, which is shaped as a convergence nozzle, clearly shown in Fig. 6(f), also observed in the study by Yoo and Sung,25 there is a gradual rise in the velocity magnitude arriving at the bubble center. Due to the interfacial density mismatch, the pressure preceding the upstream region of the bubble increases with the air-jet causing the lighter helium gas with the bubble to move to the right and the heavier air to move to the left. These movements are responsible for the bubble hollowing from the upstream interface and through the center. Similarly, as can be seen from Figs. 6(a) to 6(d), the pressure upstream of the helium volume is considerably higher than that on the downstream region, and this pressure differential means that the heavier air accelerates the lighter helium, thus subjecting the deformed bubble to the Richtmyer–Meshkov (RM) instability (shock-induced Rayleigh–Taylor instability). It is understood that vortex rings in the post-shock are formed due to the RM instability at the bubble interface. The simulated images in Figs. 6(e) and 6(f), in the present study, show a visible appearance and the formation of such vortex rings, demonstrating the existence of the RM instability. Both the simulated images and shadow-photographs in Figs. 6(e) and 6(f) also show that the air-jet impinges on the downstream air–helium interface and eventually go through the bubble.
D. Visualization of vorticity
Vorticity production and transport arguably represents the most essential component of SBI, and Fig. 7 presents the predicted contours of vorticity at different stages of SBI process. Vorticity is so important in SBI as it, alongside aerodynamic forces, predominantly controls the bubble’s motion and shape23 so that the visualization of vorticity can help to understand the SBI process. It is generally accepted that vorticity is generated due to a misalignment of the local pressure and density gradients, which is then deposited and transported in the flow when the shock waves (incident, refracted, diffracted, and focused waves) travel through the bubble. It is interesting to note that even at the early stage, vorticity is generated, as shown in Fig. 7(a), around the air–helium interface (bubble surface) where pressure and density gradient are located. The vorticity produced by the baroclinic mechanism initiates the deformation of the upstream surface and leads to the caving in/inversion of the upstream bubble surface toward the downstream bubble surface, as shown in Figs. 7(b)–7(d), ensuring that the related rotational fluid motion creates an air-jet that pierces the downstream bubble interface, as shown in Fig. 7(e). The magnitude of the deposited vorticity is determined by the non-collinearity of the local pressure and density gradients. As the maximum misalignment is located at the diametral plane, the maximum vorticity is thus accumulated at this location, thus a justification for the formation of a pair of vortex rings close to the diametral plane. During the process of shock wave interaction with the helium bubble, the RM instability causes the bubble collapse and deformation after which the fluids, i.e., helium and air, rotate and continue to develop to form the vortical structures as shown in both Figs. 6(f) and 7(f). As these perturbations develop, vortices produced at the interface rolls up and drags helium into a distinctive downstream (primary) vortex ring. Figure 7(f) also shows that most of the generated vorticity is situated in the downstream vortex ring.
E. 3D flow visualization
3D flow visualization will be presented in this section to further illustrate the SBI process. The notations used in the analysis are defined in Fig. 8. ui, di, and aj in Figs. 8(a) and 8(b) represent the upstream interface, downstream interface, and air-jet (at early stage of the SBI process), respectively. is, rw, and tw in Fig. 8(a) represent the incident shock, refracted wave, and transmitted wave, respectively. df, uf, and hν in Fig. 8(c) denote downstream, upstream interface, and helium vortex head (at a later stage of the SBI process). pvr and hl in Fig. 8(c) represent the primary vortex ring and helium lobe, respectively.
Snapshots of the predicted pressure iso-surfaces showing a quarter bubble (better suited to show certain feature) are presented below to elucidate the change in position of the helium bubble, generation of the various wave patterns as well as the formation of the vortical structures. There are two points, a point on the upstream edge of the bubble—point C and the highest point of the bubble—point D, as shown in Fig. 9(a), which are very relevant in understanding how the shocked bubble travel as well as the formation of the vortical structures. At an early stage as shown in Fig. 9(a), the planar incident shock has collided with the bubble face, generating a refracted shock inside the bubble. The shocked bubble is subsequently deformed, and at t = 0.84, the incident shock has propagated more than half the length of the bubble and lagged the refracted shock, as shown in Fig. 9(a). At t = 1.33, the incident shock travels past the entire bubble length, as shown in Fig. 9(b), and lags the transmitted shock culminating from the encroachment of the refracted shock on the downstream bubble edge. At this stage, points C and D have moved a considerable distance from their original positions. As previously discussed, the resulting baroclinic vorticity generation from shock motion across the bubble ensures that the associated rotational motion pulls a jet of surrounding air through the bubble center. This air-jet is first detected t = 1.33, as shown in Fig. 9(b), which also shows how far the transmitted wave has traveled ahead of the shocked bubble with the generated vorticity causing the upstream surface to continuously cave-in in the direction of the impinging air-jet and toward the downstream surface. This air-jet that penetrates the helium volume and pierces through the downstream edge of the bubble is analogous to the so-called spike of the RM instability at a perturbed gaseous air–helium interface. This caving-in/inversion continues, as shown in Figs. 9(c) and 9(d), until point C is no longer visible in Fig. 9(e). At the later stage (t = 3.73 and t = 5.76), the inverted part of the upstream surface has impinged on the downstream surface with vorticity spinning up and dragging the bubble fluid (helium) into a characteristic vortex ring at the downstream bubble end.
Figure 10 presents snapshots of the predicted pressure iso-surfaces showing more clearly certain features of the full bubble compression process and shock propagation. It can be seen from Fig. 10(a) that at t = 0.84 the transmitted wave is already behind the bubble while the incident shock is still propagating through the bubble, and distortion and movement of the spherical inhomogeneity has already taken place. The upstream region of the bubble has gradually become flat. The transmitted wave travels further behind the bubble, and significant distortion of the bubble is observable in Fig. 10(b). The air-jet has become evident, as shown in Fig. 10(b), and the piercing effects of the air-jet are pronounced in Figs. 10(c) and 10(d) where the upstream interface approaches the downstream interface. This is then followed by the formation of a well-defined vortex ring, as shown in Fig. 10(e).
F. Vortex ring evolution
The primary vortex ring (pvr) involves highly complex, distorted rotational motion and mixing, which have formed from the nonlinear coupling of shock compression and acceleration, nonlinear acoustic impacts, as well as generation and transport of vorticity.1,7,26,27 Hence, the development of pvr at the later stage of the SBI in terms of its formation and evolution will be further elucidated via flow visualization in this section. It can be seen from Figs. 11(a-i) and 11(b-i) that the initial spherical bubble has evolved into three distinct volumes at t = 3.73, with the smallest volume at the downstream end. This smallest volume grows bigger with time, as shown in Figs. 11(a-ii), 11(a-iii), 11(b-ii), and 11(b-iii), and as time goes by helium is drawn in and rolled up by the generated vorticity induced by shock propagation within this smallest volume with onset of the pvr’s formation observed in Figs. 11(a-iii) and 11(b-iii) at t = 4.69. This process keeps on, and at t = 5.28, a distinctive pvr has been formed, which is clearly observable in Figs. 11(a-iv) and 11(b-iv). This pvr becomes more and more distinctive, as shown in Figs. 11(a-v)–11(a-viii) and Figs. 11(b-v)–11(b-viii). Those figures also clearly show that the hl extends from the upstream edge to where the pvr exists and carries more than 60% of the helium volume, i.e., the hl is much larger than the pvr in volume.
G. Onset and development of turbulent mixing
The RM instability, which appears when the shock wave accelerates the perturbed interface separating the heavier air and the lighter helium, initially starts with perturbations that develop linearly followed by a nonlinear stage where these perturbations grow into complex structures shaped as a converging spike as seen previously in Figs. 6(c) and 6(d), showing the air-jet penetrating the helium bubble. As time progresses, the downstream end of the helium bubble is ruptured by the impinging air-jet, as shown in Figs. 6(e) and 6(f). Subsequently, more intensive mixing of helium with air occurs with some small flow structures formed, as shown in Fig. 12, indicating that turbulent mixing starts to develop at the later stage of the SBI in the present study. It can be seen from Fig. 12 that the mixing is mainly concentrated inside the pvr, the bridge region that connects the hl to the pvr as well as the outer interface of these areas. This is concordant with the findings of Tomkins et al.28 who observed three main regions of mixing, i.e., the vortex core, the outer boundary, and the bridge area joining the two major vortices. Furthermore, it is evident from Fig. 12 that large-scale entrainment occurs around the pvr, leading to the increase of pvr and reduction of hl as time passes by.
H. Turbulence generation and development
It is well known that SBI involves a broad range of complicated features, from shock wave refraction and reflection, generation and transport of vorticity, to the onset of turbulent mixing, but hardly any previous experimental and numerical studies have addressed one very important feature—turbulence generation and development. This section focuses on the turbulence generation and development of the SBI process. Figure 13 presents contours of turbulence intensity, and it can be seen from Fig. 13(a) that turbulence starts to be generated initially around the bubble interface at the early stage when t = 0.84. As times passes by, the bubble interface is severely distorted, and more turbulence is generated around the tip of the air-jet, as shown in Fig. 13(b), with the maximum turbulence intensity reaching about 11%. After the air-jet pierces through the downstream bubble interface turbulence is mainly generated around the pvr region, and maximum turbulence intensity level increases gradually to about 14% when a distinctive pvr has already been formed at t = 5.76, as shown in Fig. 13(d). Afterward the turbulence field looks similar, and the turbulence level increases gradually, as shown in Figs. 13(e)–13(g), to a maximum turbulence intensity about 20% at t = 9.9, as shown in Fig. 13(h).
IV. CONCLUSION
The mechanism of a spherical bubble deformation and compression from the interaction with an incident shock at Ma = 1.25 has been examined through a numerical study. The URANS computation approach is employed with a coupled level set and VOF method to capture the helium bubble and air interface accurately. Both 2D and 3D simulations have been carried out and it is demonstrated that the 3D predictions are much closer to the measured values, with a very good agreement between the predicted velocities of refracted wave, transmitted wave, upstream interface, downstream interface, jet, vortex ring, and the corresponding measured values. This suggests that 3D simulations are necessary for accurate predictions of SBI that is inherently 3D, especially at the later stage of SBI.
Flow visualization shows that the predictions have captured many salient flow features observed experimentally very well as the simulated images generally conform to the experimental shadowgraphs. The simulated images have clearly shown the initial distortion and flattening of the upstream edge of helium-filled spherical bubble, followed by the formation of an air-jet that grows toward and impinges on the downstream bubble edge. The evolution of wave patterns involved in the SBI, i.e., incident wave, transmitted wave, and refracted wave, is also clearly illustrated. As the SBI progresses at the air–helium interface that has already been deformed by little undulating perturbations of long wavelengths that grew linearly, the formation of the primary vortex ring induced by the RM instability is observed.
The detailed evolution process of the primary vortex ring from its onset, development to a full distinctive vortex ring, has been revealed via the flow visualization. In particular, the present study shows clearly that turbulence is generated at relatively early stage of the SBI before the formation of primary vortex ring (pvr) and after the formation of a distinctive pvr turbulence is mainly generated around the pvr region with the maximum turbulence intensity reaching around 20%.
ACKNOWLEDGMENTS
The Ph.D. studentship for Mr. Solomon Onwuegbu and the research facilities/support provided by the University of Derby are gratefully acknowledged.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.